IEEE14节点潮流程序

luojinwen

clear
r(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;
* m+ E7 F# i" Sr(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;
* J: Q7 O; f* s' x" p" G7 Pr(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;
+ e% d. }4 X" M7 a2 K% w! C6 T$ G! ar(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;
, v% W, X2 P0 D; f- [, Qr(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;
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;
r(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;

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;

& }8 c; Z$ p; J0 Nfor m=1:14
; C2 |2 e/ B5 r3 H$ Y    t(m)=0;
    for n=1:14
        t(m)=t(m)+y(m,n);
& B5 W0 z2 _) p# h    end
    if m==4
% j# M' M1 |2 G; J        t(m)=t(m)+0.190;
    end
end
for m=1:14
    d(m)=0;
    for n=1:14
        if m==n
            
2 E+ S5 k5 y! X- M        elseif (r(m,n)==0)&(x(m,n)==0)
           d(m)=d(m)+0;
; ~8 h! `! N) J- E0 R        else
           d(m)=d(m)+1/(r(m,n)+j*x(m,n));   
        end
0 b  _4 Q4 I! I  w    end
, o; R9 Q; o! l9 `7 d7 Hend
for m=1:14
    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;
        else
9 F- }9 c- X* N# ]           Y(m,n)=-1/(r(m,n)+j*x(m,n));
        end
    end
0 m4 e& A5 M7 d. X: Jend
Y;
G=real(Y);B=imag(Y);
4 M4 d; j. }3 [%给定节点电压初值和循环次数k=0
5 f* s% U' `- L! wdelt(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;
/ ^7 Z7 f+ W" {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;
2 X4 |6 S# z0 q9 H9 v" T: ^k=0;precision=1;N1=13;
& t3 R1 Z( a' Z! I, v9 h, l3 X" ~%求节点功率的不平衡量�6�2P和�6�2Q
8 u+ X3 W  Q3 l! r( ]while precision>0.00001
" u4 Y7 W+ [9 V9 y) \8 r. t" g3 `u(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     
         if m<=9
; h8 L9 t- ^  v) I; v             for n=1:N1+1
# o1 Y, }. ^1 T7 R" `1 t' @7 Fpt(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)));     
) m5 Y2 r8 \& D7 `$ o8 m2 F/ j             end
             pp(m)=p(m)-sum(pt);qq(m)=q(m)-sum(qt);
6 ^( q, O5 [: E3 \: k9 E  ?7 [         else
             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)));
             end
, w1 Z" Z+ ?+ w& ]# L# `7 w8 w             pp(m)=p(m)-sum(f);qq(m)=0;
5 p! V+ [  [- n5 f  J/ e         end      
     end
1 M/ r' M3 w: d' {, |9 x%计算雅克比矩阵各元素
# u8 A1 r1 b. T9 Lfor 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)));
- U! }$ H! O5 Q0 L% _         n0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
) e9 k8 c# ~1 j8 G         j0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
         L0(n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
2 i! [: U0 b" U( m& P# S    end
0 Q8 `# p. D" k2 _    if m<=9
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)));
( b7 ~9 m$ O, l* g- jN(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)));
* Q4 d' w  F  }/ ~9 ]J(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;
6 V* h) Q- O: V! _+ ~- ]L(m,m)=0;
    end
% M5 n8 o( g5 a7 {1 e: J; z# G: rend
, p$ Q. f. }/ U. bfor m=1:N1
    JJ(2*m-1,2*m-1)=H(m,m);JJ(2*m-1,2*m)=N(m,m);
1 h3 \! u( z; h. g    JJ(2*m,2*m-1)=J(m,m);JJ(2*m,2*m)=L(m,m);
4 }7 a8 d+ K- }end
7 o8 e- j1 ?  J' I& m7 Y/ J
for m=1:N1
    for n=1:N1   
        if     (m<=9)&(n<=9)&(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)));
+ f  Y, g% t" Z& kJ(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
; B! h3 E4 W5 [0 @N(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);
; G! A  d% P3 k8 W. jJJ(2*m,2*n-1)=J(m,n);JJ(2*m,2*n)=L(m,n);
        elseif ((m==10)&(n<10))|((m==11)&(n<10))|((m==12)&(n<10))|((m==13)&(n<10))
- R* w. \9 ^. L1 UH(m,n)= -u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
0 D7 B" M5 _. e7 K2 U' _; Q) cJ(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
* r# n1 a, S& P9 _9 EN(m,n)=-J(m,n);
. v- l0 Y; |) B% a1 I: pJJ(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))
H(m,n)= -u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
+ G- P9 B: C, J! G+ AJ(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=0;
. T; J1 T- R! y5 K. P: I9 Z6 p# xJJ(2*m,2*n-1)=J(m,n);JJ(2*m,2*n)=0;
3 z! M  n% _5 C( V. \+ {# D7 {        elseif (m>=10)&(n>=10)&(m~=n)
; A; Y: a! n6 s( O# X4 U& UH(m,n)= -u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=0;
JJ(2*m,2*n-1)=0;JJ(2*m,2*n)=0;
        end
    end
7 u6 c' ^9 y; m# L2 Z/ pend
5 S* M' G0 N! s5 |, F/ @% 对雅克比矩阵进行修正，使之变成非奇异阵，以便后续求逆计算
b=0;
  X0 _2 C% v3 O5 F: a1 p9 Z: zfor m=1:22
7 Y- O6 q  B# M; q  Q, G, u" {! W    for n=1:22
        if (m<=18)&(n<=18)
7 o1 R8 o* K$ e; s: c# b        A(m,n)=JJ(m,n);
        elseif (m>18)&(n<=18)
        A(m,n)=JJ(m+b,n);
5 g" E& O- g6 A. f6 B( O7 x        end   
    end
2 V. W! R4 r/ ^% L' k    if m>18
: p4 k* v  E  z$ k0 ?( q- S        b=b+1;
4 z& ?* H7 c  q8 f: x    end
end

* w$ L4 K  E; h8 w. T  lb=0;
for n=1:22
    for m=1:22
; n9 t* b5 A9 U4 o! e, c        if n<=18
7 P0 F  |. q+ p4 Z        A(m,n)=A(m,n);
3 |; G# b3 E" Y! Z' m; C) F' E( V        elseif (m<=18)&(n>18)
        A(m,n)=JJ(m,n+b);
4 E: V' u& m* h. f3 U7 g+ O        end
; u4 n1 A2 k; d7 ?5 ^$ Y& @    end
6 o/ D: c* r3 w    if n>18
) H/ o5 g0 p6 F' K  e. k        b=b+1;
    end
end
. X0 V, C0 Q* e. V& M
b=0;
8 [# n8 {+ t# f- }! ?; jfor m=1:22
    for n=1:26
5 \6 @- N  z3 E! X. P. }8 y        if (m>=19)&(n>=19)
$ S  N; q2 @( e4 Z" Q4 O  ^' C        D(m,n)=JJ(m+b,n);
        end
6 i7 H+ b. P# R" u. V2 A    end
    if m>18
        b=b+1;
2 F% e8 X! p' P    end
end

b=0;
for n=1:22
  j0 b  {. D  g2 n2 L    for m=1:22
# q3 V0 _. y& X9 \0 \8 m3 w        if (m>=19)&(n>=19)
) ]' Z& [$ Y, x5 v, L" X% @" @        D(m,n)=D(m,n+b);
        end
/ Y, V7 a+ A0 Q+ u    end
    if n>18
        b=b+1;
    end
end
+ g3 X+ s  ~1 D/ H
0 j- f7 t0 ?2 M- k" S& x1 zfor m=1:22
    for n=1:22
5 Q/ e3 p7 J5 R6 p  v4 s' l        if (m>=19)&(n>=19)
        A(m,n)=D(m,n);
! f- k" z& N! ~( u: t7 n        else
* z% E4 b9 X' k" F        A(m,n)=A(m,n);
        end
    end
end
3 S9 G1 `4 U' X6 j% 形成不平衡量的列矩阵
for m=1:N1
, F5 b" j: P6 O  Y* a. {2 u! u2 p    PP(2*m-1)=pp(m);PP(2*m)=qq(m);
end
C=PP(1:1:18);
* z6 k  F8 {8 H1 a" |$ S4 |% r$ S7 pb=0;
for m=1:22
    if m<=18
- C# A& P. X" d% ^5 m    C=C;
    else
    C(m)=PP(m+b);
    b=b+1;
    end
! n6 L- a/ ~6 B5 W" O+ q; n8 D+ Nend
+ O( _" Z; d& g* \% 解修正方程式，得到修正量
, A( O+ d8 U" |' w" R' Zuu=-inv(A)*C';
% z% o' f: F! W1 \. @7 R  gprecision=max(abs(uu));
$ [  V, p& J- K* x9 Q! d% 计算个节点电压新值，即修正后值
. c7 D& L6 z# dfor n=1:N1
    if n<=9
    delt(n)=delt(n)+uu(2*n-1);
    u(n)=u(n)+uu(2*n);
  X. w# l8 |$ @  ~; z$ a    else
6 b4 s( b! s3 {2 g. K    delt(n)=delt(n)+uu(n+9);
    u(n)=u(n);
    end
end
% 循环次数k加一，带入新值进入下一次迭代
. f7 \; r( P/ d6 x! M1 j" `& Y: zk=k+1;
end
6 I% |5 [' O$ t% o6 {4 V% 迭代完成后，求迭代次数、节点母线电压幅值与相角
" J* U. Z1 e1 qk-1,(delt*180/pi)',u'
* ]* v. D% y- Z* K( _% 计算各节点的功率和线路功率
for n=1:N1+1
  G6 F- v, M1 \( U# @( ]2 [    U(n)=u(n)*(cos(delt(n))+j*sin(delt(n)));
end
% k5 w" f! p: Ffor m=1:N1+1
    for n=1:N1+1
0 u( J' n0 D8 `9 y6 l- [    I(n)=Y(m,n)*U(n);
    end
    S(m)=U(m)*sum(conj(I));
end
S
for m=1:N1+1
    for n=1:N1+1
* W! u/ N0 d: U. S/ I+ c0 U7 J4 D3 S        T(m,n)=U(m)*conj((U(m)-U(n))*(-Y(m,n)))+U(m)*conj(U(m)*j*y(m,n));
    end
end
. n& o' W/ ^- v, X9 y/ Z. @T

原随风

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

luojinwen

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

lotus_cao

楼主的分享精神不错。

wangwangmoon

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

heming007

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

color

牛人啊，佩服佩服

liulina_1013

这个程序matlab论坛里就有

姗姗来迟156

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