IEEE14节点潮流程序

luojinwen

clear
. [+ s3 u, C% e! C5 K% _6 Jr(14,10)=0.01938;x(14,10)=0.05917;r(10,11)=0.04699;x(10,11)=0.19797;r(10,1)=0.05811;x(10,1)=0.17632;r(14,2)=0.05403;x(14,2)=0.22304;r(10,2)=0.05695;x(10,2)=0.17388;
r(10,14)=0.01938;x(10,14)=0.05917;r(11,10)=0.04699;x(11,10)=0.19797;r(1,10)=0.05811;x(1,10)=0.17632;r(2,14)=0.05403;x(2,14)=0.22304;r(2,10)=0.05695;x(2,10)=0.17388;
r(11,1)=0.06701;x(11,1)=0.17103;r(1,2)=0.01335;x(1,2)=0.04211;r(2,12)=0.00000;x(2,12)=0.23488;r(1,3)=0.00000;x(1,3)=0.20452;r(3,13)=0.00000;x(3,13)=0.17615;
4 @6 ^& E7 C7 `9 e1 U3 dr(1,11)=0.06701;x(1,11)=0.17103;r(2,1)=0.01335;x(2,1)=0.04211;r(12,2)=0.00000;x(12,2)=0.23488;r(3,1)=0.00000;x(3,1)=0.20452;r(13,3)=0.00000;x(13,3)=0.17615;
r(1,4)=0.00000;x(1,4)=0.53894;r(3,4)=0.00000;x(3,4)=0.11001;r(4,5)=0.03181;x(4,5)=0.08450;r(12,6)=0.09498;x(12,6)=0.19890;r(12,7)=0.12291;x(12,7)=0.25581;
r(4,1)=0.00000;x(4,1)=0.53894;r(4,3)=0.00000;x(4,3)=0.11001;r(5,4)=0.03181;x(5,4)=0.08450;r(6,12)=0.09498;x(6,12)=0.19890;r(7,12)=0.12291;x(7,12)=0.25581;
4 Z' y, a: D$ t5 O  ]& {r(12,8)=0.06615;x(12,8)=0.13027;r(4,9)=0.12711;x(4,9)=0.27038;r(5,6)=0.08205;x(5,6)=0.19207;r(7,8)=0.22092;x(7,8)=0.19988;r(8,9)=0.17093;x(8,9)=0.34802;
0 x( m& u5 L1 B5 L2 t. |  Gr(8,12)=0.06615;x(8,12)=0.13027;r(9,4)=0.12711;x(9,4)=0.27038;r(6,5)=0.08205;x(6,5)=0.19207;r(8,7)=0.22092;x(8,7)=0.19988;r(9,8)=0.17093;x(9,8)=0.34802;
( u, G; Q: U; n1 w
- H% H) ]! j" \y(14,10)=0.02640;y(10,11)=0.02190;y(10,1)=0.01870;y(14,2)=0.02460;y(10,2)=0.01700;y(11,1)=0.01730;y(1,2)=0.00640;y(2,12)=-0.31063;y(1,3)=-0.10999;y(1,4)=-0.05936;
y(10,14)=0.02640;y(11,10)=0.02190;y(1,10)=0.01870;y(2,14)=0.02460;y(2,10)=0.01700;y(1,11)=0.01730;y(2,1)=0.00640;y(12,2)=0.28951;y(3,1)=0.10757;y(4,1)=0.05752;
, h) W1 T5 T/ [# g. _
2 T7 w1 ]0 z0 _7 N9 |for m=1:14
    t(m)=0;
4 r4 Z& T+ X; A6 z    for n=1:14
        t(m)=t(m)+y(m,n);
4 K, E, C) e$ L9 w& n6 t! o    end
    if m==4
* ~3 m' ?, H  K) N( p# ~1 w7 O        t(m)=t(m)+0.190;
) G  y8 d2 b$ k) E8 U    end
end
4 W4 R* Z, R5 X& s1 Yfor m=1:14
    d(m)=0;
! }& D  H$ R- N2 v    for n=1:14
        if m==n
4 v9 h- a5 V+ Q" p. {            
" N+ H  [+ G: x  B7 b        elseif (r(m,n)==0)&(x(m,n)==0)
* H9 N' F/ T+ F- l/ v           d(m)=d(m)+0;
3 Q% e6 m6 F( z3 u! S) S8 \7 v. C        else
           d(m)=d(m)+1/(r(m,n)+j*x(m,n));   
        end
* i8 P; @  x& l7 t( Y: q2 ]    end
3 d; ]% |5 ]: `) zend
for m=1:14
; m: Z1 j! n% i, o    for n=1:14
        if m==n
           Y(m,n)=j*t(m)+d(m);
        elseif (r(m,n)==0)&(x(m,n)==0)
           Y(m,n)=0;
; E; a6 |% V, C- w5 {        else
3 Z  u. v6 j7 w$ t3 ^* c1 W           Y(m,n)=-1/(r(m,n)+j*x(m,n));
        end
    end
: _$ G3 Q+ L3 kend
Y;
8 u, |+ \" X& m/ F& O4 ZG=real(Y);B=imag(Y);
8 G4 B' O, V+ X; F0 Z%给定节点电压初值和循环次数k=0
delt(10)=0;delt(11)=0;u(1)=1.0;delt(1)=0;u(2)=1.0;delt(2)=0;delt(12)=0;u(3)=1.0;delt(3)=0;delt(13)=0;u(4)=1.0;delt(4)=0;u(5)=1.0;delt(5)=0;u(6)=1.0;delt(6)=0;u(7)=1.0;delt(7)=0;u(8)=1.0;delt(8)=0;u(9)=1.0;delt(9)=0;
p(10)=0.183;p(11)=-0.942;p(1)=-0.478;q(1)=0.039;p(2)=-0.076;q(2)=-0.016;p(12)=-0.112;p(3)=0;q(3)=0;p(13)=0;p(4)=-0.295;p(5)=-0.09;q(5)=-0.058;p(6)=-0.035;q(6)=-0.018;p(7)=-0.061;q(7)=-0.016;p(8)=-0.135;q(8)=-0.058;p(9)=-0.149;q(9)=-0.05;
" O+ q* P& N! N7 }k=0;precision=1;N1=13;
%求节点功率的不平衡量�6�2P和�6�2Q
while precision>0.00001
& L9 ?' M7 J( r& s( U1 T# g1 Ju(10)=1.0450;u(11)=1.01;u(12)=1.07;u(13)=1.09;u(14)=1.06;delt(14)=0;q(4)=-0.166+(u(4))^2*0.190;
     for m=1:N1     
) {( b, S2 T+ j8 @5 c         if m<=9
' V' g, f1 z1 H5 S6 O" C  V             for n=1:N1+1
pt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));     
% @! s/ a; P( ~# y4 ~5 ^% D5 P) O             end
) i5 D0 ], U% Y: j, c  K( }             pp(m)=p(m)-sum(pt);qq(m)=q(m)-sum(qt);
4 ~# W" y6 \' s7 S3 e2 k! }0 M- ]1 ]         else
" b7 b* D, \/ R. s2 [' K             for n=1:N1+1
f(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
) ]8 C  P. P! k; ]8 @. A2 L- p             end
             pp(m)=p(m)-sum(f);qq(m)=0;
7 A2 V, J9 C( E0 a9 b. @         end      
     end
%计算雅克比矩阵各元素
/ y! B) D; H) K7 Ffor m=1:N1
    for n=1:N1+1
         h0(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
         n0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
( b, E! j+ s/ Y% `2 y: u         j0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
( K1 X% ?$ J1 r* W         L0(n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
( [4 V9 K5 A# c$ P+ o    end
    if m<=9
$ p) b1 ?% ?0 q  _7 |/ n/ {H(m,m)=sum(h0)- u(m)*u(m)*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m)));
N(m,m)=sum(n0)-2*u(m)^2*G(m,m)+u(m)*u(m)*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m)));
# q) C& t0 T; f$ `8 u9 Y" xJ(m,m)=sum(j0)+u(m)*u(m)*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m)));
L(m,m)=sum(L0)+2*u(m)^2*B(m,m)+u(m)*u(m)*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m)));
    else
H(m,m)=sum(h0)- u(m)*u(m)*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m)));
N(m,m)=0;
J(m,m)=0;
5 s! u4 ]5 f! q9 m9 a7 {L(m,m)=0;
    end
4 L. w* d0 z: Z1 H1 A7 L# [" rend
for m=1:N1
    JJ(2*m-1,2*m-1)=H(m,m);JJ(2*m-1,2*m)=N(m,m);
    JJ(2*m,2*m-1)=J(m,m);JJ(2*m,2*m)=L(m,m);
- o; b3 ?7 a4 R2 bend
8 J) O$ [. ~' \, l
' X9 e, D, b- Z# K8 Ofor m=1:N1
8 J$ ~$ g; n5 Z# e8 V    for n=1:N1   
. h: E9 t1 t# B# I% p6 Q. {$ _        if     (m<=9)&(n<=9)&(m~=n)                       
' Z! W: t9 X# f5 IH(m,n)= -u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
J(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
( M3 X" \: y! ^/ `8 vN(m,n)=-J(m,n);L(m,n)=H(m,n);
JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n);
8 z) v! t' w" JJJ(2*m,2*n-1)=J(m,n);JJ(2*m,2*n)=L(m,n);
% F) |9 h* {- C' A% R& I& Y# r        elseif ((m==10)&(n<10))|((m==11)&(n<10))|((m==12)&(n<10))|((m==13)&(n<10))
H(m,n)= -u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
J(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
3 S* e& M' k4 @N(m,n)=-J(m,n);
JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n);
JJ(2*m,2*n-1)=0;JJ(2*m,2*n)=0;
        elseif ((n==10)&(m<10))|((n==11)&(m<10))|((n==12)&(m<10))|((n==13)&(m<10))
% q8 k! N4 m, mH(m,n)= -u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
/ |. g; Z- L/ w+ WJ(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
" k2 s8 ~$ ?) T; v7 A: x1 h; |, JJJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=0;
- s9 F1 Z9 A' [7 c2 a; mJJ(2*m,2*n-1)=J(m,n);JJ(2*m,2*n)=0;
5 o$ X" v5 e+ x# K        elseif (m>=10)&(n>=10)&(m~=n)
H(m,n)= -u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
! J" H; q+ B7 R; \" ?$ C% XJJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=0;
( c1 o+ s1 S  m$ MJJ(2*m,2*n-1)=0;JJ(2*m,2*n)=0;
        end
- d* h8 N: h; s7 Q    end
/ e/ A4 z7 N4 r1 ^7 Pend
% 对雅克比矩阵进行修正，使之变成非奇异阵，以便后续求逆计算
5 g& \  j, M/ V. u; ^4 ]5 ib=0;
for m=1:22
    for n=1:22
        if (m<=18)&(n<=18)
        A(m,n)=JJ(m,n);
        elseif (m>18)&(n<=18)
6 i! f9 U: K+ E, L1 H+ g        A(m,n)=JJ(m+b,n);
* V. ~8 s) d' F% X6 y) M        end   
    end
# C* z- U4 I$ i    if m>18
2 T! }7 N8 T' d6 y2 X" B- D5 O' ?        b=b+1;
    end
. m+ z' e$ c" F1 C* @/ Qend
4 S, _( E4 `( W* M
/ `9 N$ T2 L# _6 Bb=0;
4 w! o- w$ v' \/ q! zfor n=1:22
( T$ W6 i6 M, d. i# U% u) n    for m=1:22
9 c+ l# ^/ t8 P( d  ~, d        if n<=18
0 c/ B2 {4 ?/ E& a0 K4 D        A(m,n)=A(m,n);
        elseif (m<=18)&(n>18)
        A(m,n)=JJ(m,n+b);
4 c, q: d2 R* a2 p        end
    end
    if n>18
        b=b+1;
    end
end
2 g7 _  j! z5 y, t8 }$ M
b=0;
for m=1:22
" R6 k+ C( N& r4 Y/ `    for n=1:26
0 I  t7 {: m6 g( y& e- o0 [        if (m>=19)&(n>=19)
1 c1 E) h) Q% s) V: v" J        D(m,n)=JJ(m+b,n);
        end
- M, d8 _, {0 R( L    end
    if m>18
        b=b+1;
8 D* p! n3 x: Q! @2 [9 d. i- d, b4 g    end
end
9 o+ ^0 f- i5 {1 _; G" N
. J3 z5 _! T: y, p, i# zb=0;
! k3 [: ?2 f8 E8 l5 tfor n=1:22
) C5 ]8 g, G& r6 `2 v& a) [5 t    for m=1:22
4 E3 v( r; y4 @' ]; V! F        if (m>=19)&(n>=19)
; O) y" n3 R  O3 K. Y  G        D(m,n)=D(m,n+b);
        end
. k+ t% k# e9 k3 m* @# v+ C0 @    end
& d! A  ^$ W* _6 D    if n>18
, E; i7 _- p. N        b=b+1;
# J# T( M! [8 y2 ]+ m    end
7 a- Y% b$ P2 O4 mend

. o1 i' U4 m4 ?2 C+ b+ Sfor m=1:22
3 i- H: f1 V$ [4 S, {9 n0 k7 I    for n=1:22
6 r; u; I3 z% O8 g% i' g  `/ p- R3 {        if (m>=19)&(n>=19)
        A(m,n)=D(m,n);
4 {. }( k. l1 |: Y/ p+ Y! D$ \7 L        else
        A(m,n)=A(m,n);
) V! M6 N5 d3 b' d2 W+ R  }5 n# b' I        end
' s/ ]3 }  G* ?8 S- i6 O( S! B# R    end
9 {. O) N# ~4 F3 Cend
% 形成不平衡量的列矩阵
for m=1:N1
) Q7 [" V6 I) ~    PP(2*m-1)=pp(m);PP(2*m)=qq(m);
7 j8 u, l& S( b: k' ]$ K" Zend
C=PP(1:1:18);
b=0;
4 j1 n+ d" e! J8 p: o  Dfor m=1:22
. e# t1 I: U6 a1 |    if m<=18
6 s4 i/ l: u# \2 ]% [% c+ t2 _5 v    C=C;
    else
    C(m)=PP(m+b);
    b=b+1;
8 U1 I$ k7 R  ^$ |* T    end
; F) |) C7 R) k, \. T+ `end
% 解修正方程式，得到修正量
( x% c1 c. Z+ ^+ c* Q% _uu=-inv(A)*C';
2 T4 F: {! C3 n  K: V+ aprecision=max(abs(uu));
; {  \; w+ r- o- W" M* F% 计算个节点电压新值，即修正后值
0 p+ A5 c5 O+ @' c- Yfor n=1:N1
: A# D4 N( ]) }& M: S/ B5 d: s    if n<=9
$ i" c# n4 ]# C/ N/ m; Y+ k" D    delt(n)=delt(n)+uu(2*n-1);
    u(n)=u(n)+uu(2*n);
+ i, ^, U0 `/ ]7 ?" A# l* Y    else
! ?. O/ U. {" Z9 s, J2 D  D    delt(n)=delt(n)+uu(n+9);
    u(n)=u(n);
  l5 w1 `  }! ?: {9 ^0 L! T6 U; G    end
8 u4 X6 J3 h9 e1 l6 q- kend
% v. G" b5 r+ ]7 k' T% 循环次数k加一，带入新值进入下一次迭代
6 U; A2 Z; n5 O3 y, Kk=k+1;
end
% 迭代完成后，求迭代次数、节点母线电压幅值与相角
! \! L. c- G6 J# W% J* Kk-1,(delt*180/pi)',u'
2 |! _: |% p' S, |* \- Z% 计算各节点的功率和线路功率
for n=1:N1+1
    U(n)=u(n)*(cos(delt(n))+j*sin(delt(n)));
2 b( b. x+ {' ?1 R  \end
9 L- }* S' V7 rfor m=1:N1+1
  H1 x9 J. u2 B    for n=1:N1+1
# }- O8 _3 Y! z2 f    I(n)=Y(m,n)*U(n);
    end
* j' M7 e4 V" x6 z5 {7 `    S(m)=U(m)*sum(conj(I));
0 q/ O( [( c% \6 T" send
- x6 F' U; W( f0 n0 F7 c' |. Q6 T$ M+ XS
; k0 C7 H) l3 A# q2 ifor m=1:N1+1
. n3 E: c1 r  `) Y  G7 c1 V4 V( F    for n=1:N1+1
9 c2 ^1 T* l+ E& y        T(m,n)=U(m)*conj((U(m)-U(n))*(-Y(m,n)))+U(m)*conj(U(m)*j*y(m,n));
    end
end
T

原随风

多谢楼主分享。。。。。。。。。

luojinwen

这个好像算的结果不是很准确，但能学习下思想

lotus_cao

楼主的分享精神不错。

wangwangmoon

现在我们编程的缺点是不能用矩阵的思想全局的看待系统问题，会用很多for循环去实现功能

heming007

太复杂了，看不懂、有14节点编的算节点导纳矩阵的程序么

color

牛人啊，佩服佩服

liulina_1013

这个程序matlab论坛里就有

姗姗来迟156

多谢楼主，可以借鉴来学习啦