show下我自己用matlab编写的潮流程序

小草1987

本文的的案例是在福州大学编写的《电力系统计算程序及其实现》书中的第二页，所使用的方法也是本书中的第1，2，3章，考虑了节点优化和稀疏导纳矩阵储存问题，以下各子程序都是存在M文件中，只需把各个子程序拷贝在M文件中，分别调用即可。

5 k: \+ t: N5 M$ K7 oLGP子程序（导入数据）：

SJ1=xlsread('branch.xls')
+ V# L7 _% k1 m4 v0 U" ~3 MIZA=SJ1(:,1)       %存支路状态数%
IZ1=SJ1(:,2)       %存支路一端的节点号%
IZ2=SJ1(:,3)       %存支路另一端的节点号%
Z1=SJ1(:,4)        %存支路正序电阻%
Z2=SJ1(:,5)        %存支路正序电抗%
/ f4 B" v5 p2 B5 OZ3=SJ1(:,6)        %存支路正序电纳或非标准变比%   
" n. B2 m0 K/ g" r/ bN=SJ1(1,7)         %存网络节点数%
9 ~1 \' a8 \% r% M3 t! u" T5 sM=SJ1(1,8)         %存网络支路数%
' K+ z1 m, |# NSJ2=xlsread('generate.xls')
: u- [7 i0 A$ A' q$ lIWGA=SJ2(:,1)      %存发电机状态数%
0 U1 B) W/ d2 d; nIWG=SJ2(:,2)       %存发电机节点数%
2 ?+ k- s# w" o3 S3 B- @' N" G! bW1=SJ2(:,3)        %存发电机有功功率%
4 {+ I% H  {3 A' HW2=SJ2(:,4)        %存发电机无功功率%
. Z) Q# {# e# O: J4 V1 _+ lIQ=SJ2(1,5)        %存发电机台数%
( X. t& x- M- d. c, o; T3 s6 ySJ3=xlsread('load.xls')
ILP=SJ3(:,1)       %存负荷状态数%
8 P7 ^* {  L4 mILD=SJ3(:,2)       %存负荷节点号%
WL1=SJ3(:,3)       %存负荷有功功率%
* u9 j9 b* F! SWL2=SJ3(:,4)       %存负荷无功功率%
IP=SJ3(1,5)        %存负荷个数%
+ u5 r4 i5 K9 J$ @  b  HSJ4=xlsread('jd.xls')
  u1 E" u( v3 s# a, B$ q3 s- jN0=SJ4(1,1)        %存平衡节点的节点号%
0 t9 J2 J$ Y  _3 A5 `' R( tU0=SJ4(1,2)        %存平衡节点的给定电压值%
" O& N& I0 h: q6 B2 t6 F2 dIPV=SJ4(:,3)       %存PV节点的节点号%
- ~) e8 x: D" ?3 G0 tPV=SJ4(:,4)        %存PV节点的给定电压值%
N1=SJ4(1,5)        %存PV节点数%
% i( l0 K& u4 z! aUP=SJ4(1,6)        %存PQ节点的电压初值%
2 p7 y$ M- D* ^
LOP子程序（节点优化编号）：

for I=1:N                %寄存各节点连接支路数的ID数组要预先充零%
5 T7 Q5 u3 q3 i, e8 P( Y9 j    ID(I)=0;
end
0 _2 M( v$ n# W, S* z. R. `for K=1:M
3 }2 a/ E5 k' [( J$ R) a" e   if IZA(K)==0|IZA(K)==4  %停运支路IZA(K)=0和对地支路IZA(K)=4不属于统计范围之内%
      continue
   end
   I=IZ1(K);             %先将当前支路两端的节点号从储存它们的IZ1和IZ2数组取出，存于I,J单元%
   J=IZ2(K);
   I1=ID(I);             %将I,J节点已连接的支路数从ID数组取出，存于I1,J1单元%
   J1=ID(J);
   X=0;
   Y=0;
- t: l. w; E- B0 G6 N+ j9 @& M   for K1=1:I1            %判别当前支路是否与前面支路并联%
, l4 [$ m3 Q* B/ q      if IX(I,K1)==J
( N, M+ D2 k1 e! P/ [         X=1;
$ t% r1 I+ ?; \# I         break
      end
   end
   if X~= 0
      for K1=1:J1
          if IX(J,K1)==I
             Y=1;
6 Y1 i. {3 Q# |3 I- a9 W8 T. I( [  o          break
          end
      end
8 W3 C- t4 k( _# E5 V   end
   if X==0|Y==0            %如果是非并联支路，先将I,J的支路数各自增一，然后各自寄存对端的节点号与IX数组%
$ I+ B+ ^" X2 t9 A3 z   I1=I1+1;
   J1=J1+1;
   ID(I)=I1;
4 Y3 ^' d% d5 k8 K0 ?6 J9 H3 g   ID(J)=J1;
# {9 G# V0 e+ x. t- Q3 B; r   IX(I,I1)=J;
   IX(J,J1)=I;
   end
  c" H" [4 U6 O8 I) j' @/ }end
; f4 D; T/ |  Q8 Hfor I5=1:N                  %在尚未编号的网络中查找连接支路最少的节点%
; D  i# F) G  I% u    I1=1000;
    for I=1:N
        I2=ID(I);
& I+ O+ w& ~! y$ }# R        if I2<=1
8 Q& _8 _% J- Z1 K           I1=I2;
           J=I;
           break
        end
        if I2<I1
# `; |* d6 g& b7 `           I1=I2;
; |1 h0 |7 \1 I8 K) G% F6 q$ P/ K# _0 h           J=I;
        end
: s1 [# `' B3 `: y! L: V    end
    INA(I5)=J;
6 R% r2 U: a5 E* X    INB(J)=I5;
    ID(J)=1000;
    for J1=1:I1                         %挨个取出一个与J节点相连接的节点号，与J相连接的支路数有I1条，即与J节点相连接的节点数有I1个%                     
        I3=IX(J,J1);                        %先将IX数组取出一个与J节点相连接的节点号，存于I3单元%
        J2=ID(I3);                          %再从ID数组取出I3节点连接的支路数存于J2中%
        for J3=1:J2                         %挨个取出一个与I3节点相连接的节点号%
  X% `1 [6 @; \            J4=IX(I3,J3);                        
* S, Z6 R. ~4 w- }* C3 k            if J4==J                            %判断J4是否等于J%
               break
# _0 ^+ j1 @1 g4 y. Y            end
3 }/ h% c6 [& ]1 i- S+ F1 I! }        end
7 G! [5 y$ |4 A  E4 d; w3 J        IX(I3,J3)=IX(I3,J2);                %去掉与I3节点相连接的节点号J%
& Q0 W" z  t* Y# x        ID(I3)=J2-1;
8 i5 M8 i5 a4 i. o( E" g' |    end
    X=0;
    for J1=1:I1-1                        %去掉J节点后，使原于J节点相连接的所有节点每两个之间增加一条新支路%
        K1=IX(J,J1);                         %挨个从IX数组取出一个原与J节点相连接的节点，存于K1单元%
        K3=ID(K1);                           %从ID数组取出K1节点所连接的支路数，，存于K3单元%
5 C+ E) v- M; \4 t6 T; k/ R8 e        for J2=J1+1:I1                       %挨个取出原与J节点相连接的K1节点之后的节点，存于K2%
2 b2 [0 q: A. s3 X* X4 G4 r            K2=IX(J,J2);                        
) Z/ a( z* T  ~# \2 V            for J3=1:K3                        
                if IX(K1,J3)==K2                     %从IX数组挨个取出与K1节点相连接的节点号IX(K1,J3)，与K2进行比较，如果均不等于K2，则K1与K2节点之间无支路关系%
1 p3 g/ h8 X4 _& r                   X=1;                                 %如果IX(K1,J3)等于K2,则表明K1与K2节点之间已有支路连接，不必增加新支路%
                   break
8 d7 M, u+ x$ x% M                end
8 `! ^. m9 M- `9 U5 }7 B            end
            if X==0                              %K1与K2节点之间增加一条新支路%
                K3=K3+1;                             %K1节点连接支路数增加一%
$ X9 o8 y7 z9 r% f+ F* d4 `* B                IX(K1,K3)=K2;                        %寄存对端节点%
. a( Y- Z) Y6 L( L+ l% V2 u" k7 W                ID(K1)=K3;
                K4=ID(K2)+1;                         %K2节点连接支路数增加一%
) q5 x. E9 X9 F% }                IX(K2,K4)=K1;                        %寄存对端的节点%
                ID(K2)=K4;
            end
        end
    end
* \9 V, x3 r, b7 }end

LKP子程序（优化编号后的新节点储存）：

+ q$ R0 ~( v: s* F7 _. k8 }% Efor K=1:M             %将支路原有旧节点号换成新的节点号%
    I=IZ1(K);
    J=IZ2(K);
    IZ1(K)=INB(I);
    if J==0
        continue
    end
    IZ2(K)=INB(J);
end
for K=1:IQ            %将发电机旧节点号换成新的节点号%
* k4 I/ a0 c8 ], O    I=IWG(K);
    IWG(K)=INB(I);
" C6 J5 @7 l2 C7 l5 L* Send
/ d' x- j" |7 m3 u: |. afor K=1:IP            %将负荷旧的节点号换成新的节点号%
3 a: q, L3 d, L    I=ILD(K);
( O% Z. W7 Q, w0 H    ILD(K)=INB(I);
, w- \1 G  L' w- D; c* Tend
8 [4 y1 e0 c& k5 r" b  v' `for K=1:N1            %将平衡节点的旧节点号换成新的节点号%
    I=IPV(K);
    IPV(K)=INB(I);
end
3 w' T- K6 G" x1 zN0=INB(N0);           %将平衡节点旧的换成新的%
+ d. M0 _5 e5 @  w
: Y  f( z& g# E! k9 P/ XLDP子程序（形成导纳矩阵）：
" }. F0 j+ f  s0 @  \% d0 l; I
8 q$ E0 |6 c4 d' Ffor I=1:N                 %储存自导纳的D11,D12数组要预先置零%
    D11(I)=0;
    D12(I)=0;
end
L=0;                      %非零互导纳元素的计数单元，开始置零%
for K=1:M
    IG=IZA(K);            %IG是当前支路状态数的临时寄存单元%
  R1 l2 x: n3 Q" c4 J    if IG==0
        continue
    end
: M' C5 ], j' m9 e/ w    I=IZ1(K);
) M' R/ c' g) @# p+ b! P8 e/ j    J=IZ2(K);
    R=Z1(K);
    X=Z2(K);
% n, Z8 O6 Q% \1 F8 n2 P6 k, o    B=Z3(K);
    A=R*R+X*X;
& |7 e+ ?6 u% s* x' G2 v    if A==0
        continue
    end
   GIJ=R/A;              %计算支路导纳%
+ u8 @+ R- [0 D   BIJ=-X/A;
   Y=0;
   if IG==1   
       Y=1;
+ a5 L" R5 V! O3 {! W3 f       GI=GIJ;           %计算输电线支路I节点自导纳，J节点自导纳和它们之间的互导纳%
3 I% M1 p( @* M! f* H3 k9 l       GJ=GIJ;
3 c! F+ v1 Z# E  h4 j& x       BI=BIJ+B/2;
       BJ=BI;
9 q) k  ?3 X; |   end
   if IG==2|IG==3
       Y=1;
       GJ=GIJ;            %计算变压器支路I节点自导纳，J节点自导纳和它们之间的互导纳%
       BJ=BIJ;
) g: G, d( [9 b# r# b       GIJ=GIJ/B;
       BIJ=BIJ/B;
       GI=GIJ/B;
       BI=BIJ/B;
& h, m! h7 m0 W+ ~- H   end
; y* |  g6 c/ F7 R5 _& }3 f2 G   if Y==0
       D11(I)=D11(I)+GIJ; %对地支路时，只将当前支路的导纳累加到I节点的自导纳上%
       D12(I)=D12(I)+BIJ;
! h. t$ n( W4 Y9 z( _' G# r# f       continue
1 `6 `& u9 y3 j* n; N, r2 e   end
, o) j2 C, U4 z5 M( j   if Y==1
. S& \  w0 l8 C/ j$ U& [       D11(I)=D11(I)+GI;  %非对地支路时，将当前非对地支路导纳累加到I,J节点的自导纳上%
5 b+ E7 D; o' o4 C. y' ?, E       D12(I)=D12(I)+BI;
# ~- Q8 T  G# j, F2 }       D11(J)=D11(J)+GJ;
1 a5 i# p+ ]% `7 G% f5 \% W       D12(J)=D12(J)+BJ;
       L=L+1;             %L是非零互导纳元素的计数单元%
       YZ1(L)=-GIJ;
0 E! R% l  H4 ?: r/ \       YZ2(L)=-BIJ;
       IY1(L)=I;
6 |8 v# J6 l* {1 |( r$ @; G( f' K       IY2(L)=J;
! I7 s( r* u9 Y" q7 i" g4 r3 N! H       L=L+1;             %L是非零互导纳元素的计数单元%
       YZ1(L)=-GIJ;
       YZ2(L)=-BIJ;
0 U# m  w. o9 a& m1 {3 m) @# B       IY1(L)=J;
* Z$ j+ Q2 Y$ ^4 r: g       IY2(L)=I;
   end
6 z4 ]& q5 u: z3 vend
7 U# k4 f; Q+ {% Z% U6 LJ=0;                   %J是有规则非零互导纳元素的计数单元，挨个累计%
: a; }8 u/ e+ O" {K0=0;                  %K0是有规则非零互导纳元素的计数单元，挨行累计%
for I=1:N              %I循环实现按行号由小到大将非零互导纳元素排列在Y1,Y2数组中%
    J1=0;              %J1是当前行I非零互导纳元素的计数单元，开始置零%
    for K=1:L          %K循环挨个检查不规则非零互导纳%
        if IY1(K)~=I
$ @+ C! a4 b$ g* b' U$ f0 b  U           continue
% Z) I: N9 a0 K% f6 k8 _        end
        J3=IY2(K);     %IY1(K)如果等于I，则表明该非零互导纳YZ1(K),YZ2(K)是第I行的元素，将该非零互导纳的列号从IY2(K)单元取出，存于J3单元%
; N# k* l, M* I2 `5 j+ y        Y=0;
        for K1=1:J1    %K1循环对当前行I已有规则排列在Y1,Y2数组中的互导纳元素进行挨个排查，是否有列号等于J3的元素%
            K2=K0+K1;
            if J3==IY(K2)
: |: y/ T( x$ y3 R               Y=1;
4 B: G% f& A6 ^' {, K; t; y$ y               break
            end
( ^" T9 y0 }+ h0 t5 _        end
) X0 m) ]$ r0 @: j; }        if Y==0        %不存在列号等于J3,即非并联支路%
           J=J+1;
           Y1(J)=YZ1(K);
0 y  t/ w* j  x: _6 o+ z" }# O  W           Y2(J)=YZ2(K);
9 p/ ?( o1 q4 d. F0 z; E2 y4 A% Z           IY(J)=J3;
8 c3 R" {: E' R: \3 S$ |0 _* _           J1=J1+1;
           continue
        end
- Q/ P7 D4 E0 a! A: c. T$ r        if Y==1        %存在列号等于J3,即并联支路，将非零互导纳YZ1(K)，YZ2(K)累加到Y1(K2),Y2(K2)单元中%
$ g" R% q' W* S) t/ U           Y1(K2)=Y1(K2)+YZ1(K);
0 M7 _0 A: V# b. U/ O- a) N           Y2(K2)=Y2(K2)+YZ2(K);
        end
    end   
    IN(I)=J1;           %将当前行I非零互导纳元素的个数J1存于IN数组，IN数组是用来存放正序导纳矩阵每行非零互导纳元素个数的%
% [3 R7 _# r! W    K0=K0+J1;           
0 r5 V/ y2 ~; j  d& l  q! a& G' K! iend
LIP子程序（PV,PQ，平衡节点赋予初值）：

for I=1:N           
    U1(I)=UP;
    U2(I)=0;
! W- b( l& N* D7 Z' g    PD(I)=0;
1 |" Y+ m, D7 Z  x& J* v9 |& }    QD(I)=0;
' R4 }( h# j. e    PF(I)=0;
6 I: Y+ o, t0 {/ W: u. n  p7 a    QF(I)=0;
0 [! u7 |" Z2 n5 m0 b( t4 l4 w    IVI(I)=0;
( \, o. P# Z+ {# S. lend
2 }- z+ r7 d& [. Vfor I=1:IP                  %PD,QD数组中将接有负荷的节点填上给定的负荷功率%
    if ILP(I)==0
5 V8 w/ D: b" _) H1 _- Q& M3 M        continue
    end
    J=ILD(I);
3 R9 b( W3 H' O* ]6 e! n0 {    PD(J)=WL1(I);
' j! e: e* I" h& J    QD(J)=WL2(I);
) B) L, E  x: m0 j6 ]end
3 {& X; a& x' tfor I=1:IQ                  %在PF,QF数组中将接有发电机的节点填上给定的发电功率%
8 m  W3 l+ P3 c! }    if IWGA(I)==0
+ _+ |5 J1 ~/ V, a        continue
    end
    J=IWG(I);
8 m+ G9 ]) @8 J    PF(J)=W1(I);
; n$ J3 U+ Y1 F3 z9 A0 W, A    QF(J)=W2(I);
: }' S, T& o# s7 |/ W1 aend
for I=1:N1                  %给PV节点加标志1%
5 ?4 E: n9 F* n+ f8 T    J=IPV(I);
( E8 `) s% u5 i, L9 s    U1(J)=PV(I);
    IVI(J)=1;
end
U1(N0)=U0;

LJP子程序（牛顿拉夫逊法解雅可比矩阵）：
5 Q& {, i- b$ G% E
0 I; A- b9 x6 |8 W% P+ ~for IT=1:20                   %IT是迭代次数储存单元%
. L2 ^% f1 |# a    AM=0;                      %AM是用来寄存节点功率误差的最大绝对值%
    K0=1;                      %K0是导纳矩阵非零互导纳元素的指针，开始置1%
    for I=1:N                  %每次形成雅可比矩阵的一行元素及其相应的常数项，然后进行消去和规格化运算%
        A=D11(I)*U1(I);
0 e2 ?/ w# w. N1 N9 e; ]% h        B=D12(I)*U1(I);
- I' ^2 |. H- h  @+ ]9 F! F/ a        R=D11(I)*U2(I);
        X=D12(I)*U2(I);
7 [3 t0 L  Z+ G6 A& E( c, D# H3 H. R        A3=A-X;                %A3，B3单元寄存第I节点电流的实部和虚部%
- a8 i6 s' |. U8 ?. n  w$ t5 U        B3=R+B;
- Z: t/ I3 m* @$ _) j        J5=1;                  %J5是雅可比矩阵第I行非零元素的计数单元，开始置1%
9 T5 b; ]0 m- R  ]8 ?        for IG=1:IN(I)         %IG循环表示挨次取出导纳矩阵第I行的一个非零互导纳元素，形成雅可比矩阵相应的一个非零元素%
            J=IY(K0);
            A3=A3+Y1(K0)*U1(J)-Y2(K0)*U2(J);
& c8 ~9 h" ^4 i. D$ @9 U# d            B3=B3+Y1(K0)*U2(J)+Y2(K0)*U1(J);
  ^& K% I# Z- d: Q            if I~=N0&J~=N0
                J5=J5+1;
9 o, p; A4 v) e                JK(J5)=IY(K0);
- F7 ^: n. T. A" y5 R                AK1(J5)=Y1(K0)*U2(I)-Y2(K0)*U1(I);
7 y" i5 W4 |4 ~. i# f- E. y  ^. r                AK3(J5)=Y1(K0)*U1(I)+Y2(K0)*U2(I);
) c3 F3 M0 O8 S" P. p                AK2(J5)=-AK3(J5);
7 A* A$ U% U* t' ]  W! v# M                AK4(J5)=AK1(J5);
            end
            K0=K0+1;
        end
        if I==N0               %第I行如果是平衡节点%
            GQ(I)=A3*U2(I)-B3*U1(I);
            GP(I)=A3*U1(I)+B3*U2(I);
            CK1(I)=0;
1 f7 J* y( O$ _3 A% k" N6 a5 X4 r' U" F            CK2(I)=0;
4 J2 _& M% }( x5 F4 t            JF(I)=0;
            continue
( L( T  ]# n. S2 o+ u& u9 m        end
        P2=A3*U1(I)+B3*U2(I)+PD(I);
8 k' M. T) K3 @$ |+ f) V* r  F        Q2=A3*U2(I)-B3*U1(I)+QD(I);
        GP(I)=P2;
1 x' h2 j1 B! c7 H4 g        GQ(I)=Q2;
        P2=PF(I)-P2;
        Q2=QF(I)-Q2;
- j% g) y; S' {" R        P3=P2;
        Q3=Q2;
        if IVI(I)==1           
5 d* K2 o& M8 A+ H9 s            Q(I)=Q2;
            if Q(I)>0|IT-1<5
! x2 n  T! i% Y1 P$ P& ~% }( Z                Q3=0;
            end
4 Z$ v: @% T) @7 N% A7 P; I        end
        AK3(1)=A3+A+X;
6 X/ H; o8 L4 `' [0 k$ \        AK4(1)=B3-B+R;
5 j" s. f% d1 k* w) q% a2 n        JK(1)=I;
! ~  G$ Y  T: Y2 A0 q  u        CK2(I)=P2;
        Y=0;
. H& W2 V& B8 ]1 m        if IVI(I)~=1
            Y=1;
; M8 B$ o3 H$ u0 `        end
2 h% W9 z8 ~# q        if Y==0
            Q(I)=Q2;
; c% ~6 M8 w1 L! I! y9 H1 V            if IT>=5
                if Q2<0
4 ~1 E- p- U- c/ ]- x# J+ I                    Y=1;
  h$ k, Z0 d5 H( X- d) n$ n* u                end
- D6 y% X$ U5 B: w  h4 J& S5 q            end
            if Y==0
                for J=2:J5
                    AK1(J)=0;
1 D" u$ U4 \4 o" S( @                    AK2(J)=0;
                end
                AK1(1)=2*U1(I);
- d% Q& S. s& u9 T8 _; [9 S                AK2(1)=2*U2(I);
  P) {7 z0 J% D0 D                Q2=0;
                CK1(I)=Q2;
5 S! H+ {6 F5 Z( t  g' W9 @            end
! W6 W# J% Z" V4 S6 g! I% ]        end
        if Y==1
           AK1(1)=R-B-B3;
           AK2(1)=A3-A-X;
           CK1(I)=Q2;
        end
6 h  X4 k9 s! P9 f) J        C5=abs(P3);                %判别节点功率误差的绝对值是否大于AM%
- r9 H: L5 s- x& Y% C; |        D5=abs(Q3);
        if C5>AM
           I0=I;
           AM=C5;
0 e3 }+ \$ `) M1 L% }        end
" I. k, Z# W' p: x4 Z/ ~1 ]! i8 f+ f5 P        if D5>AM
           I0=I;
; [/ g. X& n6 p! C$ s7 Z+ d           AM=D5;
9 {: H% v* ~; `, o* u2 V        end
* t# S6 q4 |: i2 X. D5 F        K=1;                  
        for I1=1:I-1           %进行消去运算%
4 _$ P# z  c2 ~" f/ G6 A! x            X=0;
3 \5 g, `8 N9 M: d            for I3=2:J5
                if JK(I3)==I1
# H# _/ a% g! l                   X=1;
                   break
                end
            end
            if X==0
# \3 P+ Z2 E8 {               K=K+JF(I1);
% {$ h0 P5 R( m               continue
+ \# n# k7 A1 J$ a" ~            end
            for IG=1:JF(I1)
/ e$ H7 Q6 R7 |6 ]                 Y=0;
/ x' f  l2 q6 j0 j1 M                 for I2=1:J5
) \. S9 x- q- b' C& Q- Q' \1 Z                     if JK(I2)==IJ(K)
! Y6 L5 ]2 ]: Z                        Y=1;
; X- ^' x& c' w' k9 n                        break
                     end
                 end
                 if Y==0
                    J5=J5+1;
                    AK1(J5)=-AK1(I3)*AJ1(K)-AK2(I3)*AJ3(K);
                    AK2(J5)=-AK1(I3)*AJ2(K)-AK2(I3)*AJ4(K);
                    AK3(J5)=-AK3(I3)*AJ1(K)-AK4(I3)*AJ3(K);
2 S5 p; K! b% e" {/ [4 D                    AK4(J5)=-AK3(I3)*AJ2(K)-AK4(I3)*AJ4(K);
                    JK(J5)=IJ(K);
                 end
                 if Y==1
                   AK1(I2)=AK1(I2)-AK1(I3)*AJ1(K)-AK2(I3)*AJ3(K);
                   AK2(I2)=AK2(I2)-AK1(I3)*AJ2(K)-AK2(I3)*AJ4(K);
& E1 _2 Q7 A; `                   AK3(I2)=AK3(I2)-AK3(I3)*AJ1(K)-AK4(I3)*AJ3(K);
                   AK4(I2)=AK4(I2)-AK3(I3)*AJ2(K)-AK4(I3)*AJ4(K);
# N0 A8 l: H' P, E; f$ l. l8 j% ?0 P                 end
0 `) f) `5 I$ Y                 K=K+1;
- E9 c! R" s  y, x  _& V              end
+ P: u8 v' y0 t& V              CK1(I)=CK1(I)-AK1(I3)*CK1(I1)-AK2(I3)*CK2(I1);
3 U# j& R7 ~$ i: h7 I: y              CK2(I)=CK2(I)-AK3(I3)*CK1(I1)-AK4(I3)*CK2(I1);
9 s2 |4 A4 d6 G* O           end
           K5=K;
           A=AK1(1)*AK4(1)-AK2(1)*AK3(1);
           if A==0
! Z! G' U' R, e2 ?, H  A7 l' u                CK1(I)=0;
                CK2(I)=0;
( }& e9 x# E) {1 M9 }' Z: g& _, O4 P+ H                JF(I)=0;
                continue
! H2 u. v4 k1 l+ k  `1 V           end
           A3=AK4(1)/A;
           B3=-AK2(1)/A;
           C3=-AK3(1)/A;
           D3=AK1(1)/A;
        for J=2:J5
            if JK(J)>I
               IJ(K)=JK(J);
, I% t+ n" S6 D               AJ1(K)=A3*AK1(J)+B3*AK3(J);
* `/ H; a% ]) J. B& m               AJ2(K)=A3*AK2(J)+B3*AK4(J);
               AJ3(K)=C3*AK1(J)+D3*AK3(J);
1 ^- U) |2 H* u  F               AJ4(K)=C3*AK2(J)+D3*AK4(J);
& M3 {2 e# Y2 f; r: l               K=K+1;
  p/ v: _2 S8 \. A            end
        end
  \/ v4 X& P% t& s. l# k) @2 V        A5=CK1(I);
        B5=CK2(I);
        CK1(I)=A3*A5+B3*B5;
1 f8 y- b5 w& i2 r2 p1 P        CK2(I)=C3*A5+D3*B5;
8 T: K# u3 _* F# h% K6 y# R% X        JF(I)=K-K5;
  ^# @5 R/ q1 h# G8 }/ t    end
5 p' R/ K2 ?  F) [/ B* ~6 Z4 F        fprintf('%6.2',AM)
5 s: E9 a. ~) ?        fprintf('%6.2',IT)
        fprintf('%6.2',I0)
        if AM<=1.0E-4
           break
; l6 ~0 N# Z% }1 I        end
+ X* a) r6 S8 t& q/ u* ^, |% X" F        if IT>20
7 V! T: k& `" Y6 d+ Y6 H+ j           break
$ j$ j! V' Q; H        end
        for I=N-1:-1:1                    %I循环表示回代运算从n-1开始，倒推到第1行，每次计算一个节点电压的修正量%
            for IG=1:JF(I)
( `  Z1 k. l. o3 a( [& M( I# a                K=K-1;
                J=IJ(K);
                CK1(I)=CK1(I)-AJ1(K)*CK1(J)-AJ2(K)*CK2(J);
                CK2(I)=CK2(I)-AJ3(K)*CK1(J)-AJ4(K)*CK2(J);
            end
7 \, [) y. B+ f* K" a' r        end
& r3 P9 e+ ~% x9 a. D& Z, a        for I=1:N                         %计算各节点的新电压%
            U1(I)=U1(I)+CK1(I);
            U2(I)=U2(I)+CK2(I);            
       end
' ]. q4 ~- Z# e) A1 V( B       for I=1:N1
           J=IPV(I);
6 x* m3 `# X3 W& d5 I0 `0 E& w           if Q(J)<0&IT>=5
              continue
           end
           C=PV(I)/sqrt(U1(J)*U1(J)+U2(J)*U2(J));
           U1(J)=C*U1(J);
           U2(J)=C*U2(J);
      end
: b: E& [+ Y3 S  end               

LRP子程序（输出节点信息和支路信息）：
* _' O) ?% n7 q2 S" C3 k; S
% G& L9 j& S! p* b. `GP(N0)=GP(N0)+PD(N0);
GQ(N0)=GQ(N0)+QD(N0);
& e5 y) L, L$ E& IPG=0;
9 T; q: a+ f2 v5 jQG=0;
PL=0;
9 a7 M! S3 O( y" u1 Z0 y( R' y! bQL=0;
PI=180/3.14159;
for I=1:N                    %I循环实现每次输出一个节点的信息，输出是按旧的节点号为顺序号%
0 G) V8 L& a+ b$ `  V  }4 A" O    I5=INB(I);
    E=U1(I5);
/ ?5 p' k% _% D3 O* M5 H1 s: s8 @    F=U2(I5);
    A=sqrt(E*E+F*F);
; e* Q' p( r: @& [    B=PI*atan(F/E);
0 y4 e+ X( r$ b5 I& u    PG=PG+GP(I5);
    QG=QG+GQ(I5);
    PL=PL+PD(I5);
    QL=QL+QD(I5);
    JD(I,1)=I;
% y$ b: ]  k1 M0 O0 D  \  M7 H    JD(I,2)=A;
/ A; S/ n% T6 R1 k( l2 ?0 j    JD(I,3)=B;
    JD(I,4)=GP(I5);
* j* W& d8 A: e+ ]% L    JD(I,5)=GQ(I5);
( S7 {) I' g4 G; i+ y( c    JD(I,6)=PD(I5);
    JD(I,7)=QD(I5);
end
3 p! ^% H, m! x1 A6 I( s; ~  |PG
  `3 o1 q' Z7 d8 ]/ f, J. d1 v# KQG
, i4 R- O1 n" b1 o; o7 [( iPL
QL
PLOSS=0;
QLOSS=0;
QB=0;
for K=1:M                     %K循环表示每次输出一条支路的信息%
$ S& n: V4 \& a$ m7 P9 e- N    IG=IZA(K);
    I=IZ1(K);
    J=IZ2(K);
    R=Z1(K);
    X=Z2(K);
    B1=Z3(K);
    if IG==0
        ZL(K,1)=K;
        I=INA(I);
- c3 E( E% p6 F1 m$ X        J=INA(J);
        ZL(K,2)=I;
        ZL(K,3)=I;
% o6 a- v9 S& X6 l$ O       continue
    end
    if R==0&X==0
' V5 [* {+ S8 E; o5 F  j# j& z# `        I=INA(I);
        J=INA(J);
1 o9 v  U. _+ q2 R- P& {  n- Z        ZL(K,1)=K;
        ZL(K,2)=I
" O" [) |7 c' i6 l9 ~9 ~* `% c        ZL(K,3)=I;
        continue
+ k6 z9 {" t  R1 Y8 z& d    end
    E=U1(I);
4 q# o+ F" Q( z& I% f+ T5 J; I+ b    F=U2(I);
3 D5 v  a8 N0 N8 W    if IG==4
  ?" a( G% ?2 X9 O3 k$ e        A=0;
/ _% D, G7 \4 F) W( T. ~& ]3 @/ L        B=0;
. U# L! W# i) k/ s" D1 r    end
    if IG~=4
$ G1 N8 n6 X$ V; x        A=U1(J);
9 F7 ^! e7 G% L" Z0 G+ L2 M& a        B=U2(J);
    end
    if IG==2|IG==3
        E=E/B1;
1 U& D6 a6 P( v* _        F=F/B1;
        B1=0;
$ v5 H4 n1 |) W% I    end
' l$ X! V; y0 R, O+ V    A3=R*R+X*X;
    C3=E-A;
    D3=F-B;
    A5=(C3*R+D3*X)/A3;
    B5=(D3*R-C3*X)/A3;
- }: `, M) |2 D    C3=(E*E+F*F)*B1/2;
    D3=(A*A+B*B)*B1/2;
    P1=E*A5+F*B5;
    Q1=F*A5-E*B5-C3;
    P2=-A*A5-B*B5;
    Q2=-B*A5+A*B5-D3;
    if IG~=4
        PLOSS=PLOSS+P1+P2;
        QLOSS=QLOSS+Q1+Q2;
! s" J) f0 L5 g, w7 u* B0 E, V        QB=QB-C3-D3;
    end
3 ]' F% K- m0 H% ?4 f: l7 R    I=INA(I);
    if J~=0
        J=INA(J);
    end
    ZL(K,1)=K;
    ZL(K,2)=I;
    ZL(K,3)=J;
    ZL(K,4)=P1;
    ZL(K,5)=P2;
    ZL(K,6)=Q1;
    ZL(K,7)=Q2;
end
) ^6 j8 C* k$ a7 u( s* TPLOSS
6 Z% d+ P4 ?- l% n. UQLOSS
QB
6 V5 G' [/ i9 K; X$ J. @( T7 Vxlswrite('outext.xls',JD);        
3 C4 J+ x* D. ?: Axlswrite('outext1.xls',ZL);

gxuyang

我想看看，学习一下

zgk001

那个excel文件读入是怎么处理的？对于节点数有没有要求？

lhw606

不错，支持原创~~~~~~~~~~~~

klein

有些复杂，没弄明白excel是怎么调用的？

小草1987

回复 5# klein


4 d: M; h1 L6 R  j: T2 v3 b    福州大学编写的《电力系统计算程序及其实现》书中的第二页给了个表格，你把它输入到excel文件中，matlab中xlsread就可读取

原随风

这样的帖子要顶

fzhang1978

学习一下，谢谢

rena1121

好,力顶！！！

wz0430131

看一下，呵呵