MATLAB 电力系统潮流计算程序---连发

沧山牧水

首先声明一下，这些是从网站上转载的，不是本人上编写的
基于MATLAB的电力系统潮流计算

%简单潮流计算的小程序，相关的原始数据数据数据输入格式如下：
8 T$ l. N& L+ H7 W7 w+ Q%B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写
8 R) X7 B' H1 }9 x%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点
%编号，将变压器的串联阻抗置于低压侧处理。
- y  k8 c' z' u/ p6 ~& Y: e9 N%第三列为支路的串列阻抗参数。
) l0 I3 r8 P+ ?% I! p%第四列为支路的对地导纳参数。
%第五烈为含变压器支路的变压器的变比
% e) k, E6 m# `% j; ~%第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，
# }, e0 y& ?; `%“0”为不含有变压器。
* i  }8 u3 r1 P! {8 Z" {* a' d) Y%B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点
%负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中
  l1 B& r0 ]5 B( V: }, O& s. R7 S%“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。

' C& W& w; k3 ~2 _# I3 t%X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地
%参数。
* Q7 ^5 b. G  o  [n=input('请输入节点数:n=');
  o9 b( @* Z. o3 ^7 V  h/ kn1=input('请输入支路数:n1=');
8 G) A+ V4 j6 T9 Risb=input('请输入平衡节点号:isb=');
" G5 x8 R# w- q. Y2 _pr=input('请输入误差精度:pr=');
B1=input('请输入支路参数:B1=');
B2=input('请输入节点参数:B2=');
! k9 V: g$ X% X& @, P+ mX=input('节点号和对地参数:X=');
Y=zeros(n);
Times=1;                %置迭代次数为初始值
%创建节点导纳矩阵
for i=1:n1
    if B1(i,6)==0       %不含变压器的支路
        p=B1(i,1);
8 `1 g8 u5 h9 Y# V- N5 A" ]        q=B1(i,2);
        Y(p,q)=Y(p,q)-1/B1(i,3);
        Y(q,p)=Y(p,q);
0 K. g  u" H6 v7 m5 X        Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);
+ A( c% r4 J- p( _1 x        Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);
4 B3 N! ^: P# X* P8 V( n8 ~    else                %含有变压器的支路
1 j7 q- N4 m) H5 u: f        p=B1(i,1);
        q=B1(i,2);
        Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5));
7 J' Z& W! p  y: C+ D6 c        Y(q,p)=Y(p,q);
        Y(p,p)=Y(p,p)+1/B1(i,3);
        Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));
    end
end
* h9 n+ N! W* n$ P  L9 qY
! `9 G( [! ]' |/ X! J0 l1 F) p) WOrgS=zeros(2*n-2,1);
, X" V* V, J* |/ f( A) O7 g1 Y8 v9 vDetaS=zeros(2*n-2,1);   %将OrgS、DetaS初始化
" X9 Y1 ^9 C2 s1 p1 n! Z%创建OrgS，用于存储初始功率参数
- ^" F7 v8 C3 `$ wh=0;
j=0;
9 r0 ~& K' ~9 G7 pfor i=1:n            %对PQ节点的处理
    if i~=isb&B2(i,6)==2   
9 ^' n% ~* \" N- C        h=h+1;
3 q9 c4 k3 w; ]: `7 M        for j=1:n
            OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
            OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
3 Z1 c* w2 h6 q& }' l( L/ A. ^6 \        end
    end
end
for i=1:n           %对PV节点的处理，注意这时不可再将h初始化为0
    if i~=isb&B2(i,6)==3
        h=h+1;
9 Q6 E! o  {2 k4 a! D! p. U  @        for j=1:n
            OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
7 ]; i/ z/ W7 j* A3 n; G4 l- m1 M            OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
        end
    end
8 N4 b' @% N' i  Z" Tend
OrgS
- [" l8 A* A7 Y%创建PVU 用于存储PV节点的初始电压
PVU=zeros(n-h-1,1);
5 V% v% k- P' q8 Lt=0;
# S/ J) ^+ z- ^  Q: @* ?. a! Pfor i=1:n
3 \) [1 T+ x/ q% j    if B2(i,6)==3
        t=t+1;
        PVU(t,1)=B2(i,3);
0 ^$ _* |/ L2 R( h3 S    end
" b% w% ]5 W" R% Qend
8 z& C. S/ ^* F% f) ZPVU
7 R) C' o9 G+ ]8 h! t: U7 L%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量
h=0;
/ x) v+ S2 V! e3 s" |for i=1:n           %对PQ节点的处理
* |+ J- }% r2 k1 K3 o, K' Q+ W    if i~=isb&B2(i,6)==2
        h=h+1;
8 j1 K' |3 u, {) V2 w/ Q* }        DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
        DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);
& c6 U5 U7 l# y% M$ s& _6 Y    end
end
t=0;
for i=1:n           %对PV节点的处理，注意这时不可再将h初始化为0
* T+ P) }" m4 s( n- j: u  l/ ~    if i~=isb&B2(i,6)==3
! y- l' H3 @5 o        h=h+1;
        t=t+1;
" [0 M9 F, g$ |6 T  @' x        DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
        DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;
    end
end
3 ~8 E* P3 M3 j0 d, k9 \DetaS
%创建I，用于存储节点电流参数
0 J- @8 ^/ ~9 T( `* Z$ Qi=zeros(n-1,1);
h=0;
1 T2 ^% C  r4 b) q, k2 t# dfor i=1:n
# U& q. X4 A" H! r/ ~    if i~=isb
, K- l1 n  V) n5 o7 i! R        h=h+1;
        I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));
2 e1 x9 f, S& z    end
end
I
%创建Jacbi(雅可比矩阵)
1 i7 }5 Y. x! A  ~7 C& vJacbi=zeros(2*n-2);
9 r. t9 M8 t5 ch=0;
. ~4 e# F; ]8 x6 b# }' `k=0;
7 n5 Y% l$ z& Tfor i=1:n       %对PQ节点的处理
: O! q+ f% {! j) R( x! _3 f2 h+ ]" f    if B2(i,6)==2
2 B: u1 S4 s6 ]: \& ?" u1 W        h=h+1;
- }; k3 q3 _  r8 d! J        for j=1:n
* t: A+ g* @& E# @4 O7 F/ Q            if j~=isb
9 C& w& t2 k, h' ?' D                k=k+1;
                if i==j     %对角元素的处理
  [/ v! g# U5 ]3 Y% a% \                    Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));
                    Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));
                    Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));
                    Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));
6 `+ z0 I( ?4 h% b8 }                else        %非对角元素的处理
                    Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));
/ Y: j' \/ P8 c; }7 ?                    Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));
2 C( d7 O0 r. a; b1 O$ ]) T0 N2 N. T                    Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);
/ r/ b2 O  k) `5 c+ I1 t$ z; \                    Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);
                end
" C* d2 P# @6 x- X8 ]3 y$ t                if k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行
/ B6 p! X7 W1 ^$ N- Q' q                    k=0;
, {* w: C. b- `9 j5 s' A! B                end
            end
        end
0 \) l4 W) ~+ z/ R% \    end
) J3 D9 ?8 l( R) S- M) Z' F- H! k  b- u& o' Xend
$ _  K; l7 n- B; A' K0 H& |0 ^8 zk=0;
8 W! h& G: A7 G$ [8 `for i=1:n       %对PV节点的处理
    if B2(i,6)==3
* \; Z+ K4 C( P& U/ Q* A7 V; ?  \        h=h+1;
# R0 j( `9 d. R9 B- K# a% S        for j=1:n
; P6 X1 s0 e# u, z. b            if j~=isb
  u; j7 E5 R/ E- l7 q2 z                k=k+1;
                if i==j     %对角元素的处理
                    Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));
9 G5 \! N3 \% v9 X1 U3 W' {                    Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));
, F4 B9 t- y+ _7 _                    Jacbi(2*h,2*k-1)=2*imag(B2(i,3));
  U9 x4 M0 t* X                    Jacbi(2*h,2*k)=2*real(B2(i,3));
                else        %非对角元素的处理
                    Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));
                    Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));
) j$ H+ a0 S) I8 Z4 \/ e. ?                    Jacbi(2*h,2*k-1)=0;
                    Jacbi(2*h,2*k)=0;
, R5 s* F, G+ h5 @                end
                if k==(n-1)     %将用于内循环的指针置于初始值，以确保雅可比矩阵换行
                    k=0;
                end
: P  N! l4 Q7 K+ O6 u            end
        end
    end
end
Jacbi
2 a' V6 n1 [! {# p%求解修正方程，获取节点电压的不平衡量
8 v, _, j0 T9 q$ f8 ]! RDetaU=zeros(2*n-2,1);
% W5 ~' d( ]5 O2 ^; j4 TDetaU=inv(Jacbi)*DetaS;
  }+ |! z6 s1 n1 C/ }DetaU
%修正节点电压
8 H" R; |+ M0 P* mj=0;
for i=1:n       %对PQ节点处理
    if B2(i,6)==2
& ~) h" ]% M$ n' `        j=j+1;
        B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
4 ]* N$ u( ~4 n& w6 G1 q) K    end
; O; |- l% T# M  d. {' z" E7 Gend
( U6 B8 X: D' Dfor i=1:n       %对PV节点的处理
    if B2(i,6)==3
. [5 T# ]" o+ x% h4 ~0 A        j=j+1;
        B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
    end
! X+ N8 f! V! Z. ?7 zend
, O0 y; V. _) }6 B4 N+ O" ~& FB2
6 P' y( f2 S5 x* O  J8 H0 H5 y1 h%开始循环**********************************************************************
' |& `& W5 h+ \+ twhile abs(max(DetaU))>pr
0 J! Z4 t5 U# N; z# h9 SOrgS=zeros(2*n-2,1);        %!!!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵
h=0;
* d( C7 [6 k. K, m+ I6 Q" r* Sj=0;
% r1 {4 ]# e9 \& c! T- ~( kfor i=1:n
: L* d" x4 R! @    if i~=isb&B2(i,6)==2
' r2 a& u$ M( [        h=h+1;
        for j=1:n
            OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
            OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
$ Z) Q/ B0 u4 f, x* Z        end
    end
1 h6 \8 B5 d4 t1 r) l4 C) tend
" p+ i9 ^  `% Q9 G: }for i=1:n
$ [5 h3 M3 W  T. [- T% }    if i~=isb&B2(i,6)==3
        h=h+1;
        for j=1:n
1 c3 l# R2 k/ d5 z- J/ R$ Q4 d+ h            OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
            OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
5 d! h4 n! T% s# z& L; @7 K        end
6 D, n+ ~/ q. L+ [, B    end
end
OrgS
9 ^" p9 ]: m6 h+ Y/ n) s, m%创建DetaS
5 c% E& }. ^6 J: {$ ^9 @  U- M" ^' Wh=0;
2 Y8 ]' n( b9 f3 xfor i=1:n
/ B) h, f6 C& E0 ~1 A2 R* W" I    if i~=isb&B2(i,6)==2
        h=h+1;
        DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
. C5 B/ T: Q$ Z1 V* |        DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);
    end
' R  n6 W' J* t( x# [% t2 `end
t=0;
for i=1:n
4 i7 ^; l5 L! A8 F4 X1 U    if i~=isb&B2(i,6)==3
        h=h+1;
* {9 t6 J2 j2 \" Q, K( g        t=t+1;
        DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
; R6 K% ~2 N1 S, h4 F* N        DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;
    end
3 A% u+ m- M  W% vend
DetaS
%创建I
i=zeros(n-1,1);
% t8 f% h' M/ l7 `& {h=0;
for i=1:n
: x& l/ v7 R* L    if i~=isb
/ J: _* X5 A5 ^4 G' ?        h=h+1;
, V9 ^: h' C- g8 q: T        I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));
  ~! J3 X/ W) a6 l    end
end
I
# _: |; z$ U0 j' ?5 d6 C3 E& a2 e%创建Jacbi
& M- F4 C3 i. y8 Y9 c7 nJacbi=zeros(2*n-2);
1 W/ }' i5 t5 A, Kh=0;
k=0;
. N' Z7 v) |. M4 C: H. I' j+ afor i=1:n
2 x; B9 u  ?$ Y    if B2(i,6)==2
        h=h+1;
        for j=1:n
            if j~=isb
8 u* Q& N8 {6 F3 M; e                k=k+1;
/ E( d/ F9 }( f# e                if i==j
                    Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));
( K. n6 B0 L8 A" Y                    Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));
+ _. p( }: Q) u6 {$ T1 h! @& Z/ ]/ a: D                    Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));
% x" j5 @0 W. ?7 a3 p                    Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));
* B- G) e6 W. S2 q- r9 g                else
                    Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));
                    Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));
                    Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);
                    Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);
+ b1 E. l. e8 D% ~' Q% v5 q                end
                if k==(n-1)
                    k=0;
                end
1 [7 J; B* V! I6 ]            end
$ M" R- `, W, d# ~  k" |        end
) ~& T, ~# \- E1 `    end
3 [/ k7 m+ l% t# @7 w( M1 lend
2 Y8 ^! V* b1 `/ B( t9 Jk=0;
- W6 W5 M' b+ K2 [for i=1:n
0 T: Q$ D3 w- n    if B2(i,6)==3
/ o5 H9 M7 K* U7 ?2 I7 W        h=h+1;
        for j=1:n
3 C: p& g; G# k            if j~=isb
4 q5 k% g; Y6 _$ C$ G2 m( S                k=k+1;
                if i==j
                    Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));
                    Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));
                    Jacbi(2*h,2*k-1)=2*imag(B2(i,3));
8 W' E4 u, p3 [0 s8 f- x- `% O                    Jacbi(2*h,2*k)=2*real(B2(i,3));
6 z3 l9 ?/ x+ g% L2 `                else
  n# s" K( l- L6 K) _) |                    Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));
% p2 S1 o" r6 t& m& x1 d                    Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));
- B6 I9 L( ]; ?' ]                    Jacbi(2*h,2*k-1)=0;
& s! R$ N7 ^+ w                    Jacbi(2*h,2*k)=0;
* j3 ?: _! N" K' x) j                end
                if k==(n-1)
                    k=0;
                end
            end
! l& e% t( P# ?4 p0 b3 h        end
& R: y" d2 h: t    end
! ?2 ^% o- w7 G9 t) H5 D1 L7 ?end
# t. V1 r% P' N) ?Jacbi
DetaU=zeros(2*n-2,1);
3 z9 K- z3 e- K9 ~* F  `DetaU=inv(Jacbi)*DetaS;
DetaU
%修正节点电压
j=0;
for i=1:n
    if B2(i,6)==2
        j=j+1;
8 n6 Y2 g5 f; W$ G5 O# C) }4 p        B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
    end
9 e, `. ]: L. W7 F4 aend
for i=1:n
: {, F" b% t: k8 A    if B2(i,6)==3
! G; D8 z5 D$ O        j=j+1;
        B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
" Z4 `$ s( l, k    end
$ R8 o- e4 P* O) Gend
B2
Times=Times+1;      %迭代次数加1
& \- F4 z8 V/ Y$ |' P7 @% }end
# i% `9 H/ G4 q5 H6 `Times

$ k6 @. }, u$ T/ ?6 ~一个原始数据的例子
节点数 5
支路数 5
平衡节点编号 5
精度pr 0.000001
B1(支路参数矩阵)
[1 2 0.04+0.25i 0.5i 1 0;1 3 0.1+0.35i 0 1 0;2 3 0.08+0.30i 0.5i 1 0;4 2 0.015i 0 1.05 1;5 3 0.03i 0 1.05 1]
B2(节点参数矩阵)
" K6 x8 R/ [& S, }2 ~8 W1 M[0 -1.6-0.8i 1 0 0 2;0 -2-1i 1 0 0 2;0 -3.7-1.3i 1 0 0 2;0 5+0i 1.05 1.05 0 3;0 0 1.05 1.05 0 1]
X(节点号和对地参数)
[1 0;2 0;3 0;4 0;5 0]

沧山牧水

程序2：
# [! `! J/ H& D; W7 e- |%本程序采用牛顿-拉夫讯法对一个五节点电力系统进行潮流计算
% n=input('请输入系统节点数目:n=');
/ s! q/ ]8 J: J3 b; A% nl=input('请输入系统支路数目:nl=');
, m% Z$ r1 k( w( m2 y, R. B% ph=input('请输入平衡节点号:ph=');
% jd=input('请输入误差精度：jd=');
7 O, @8 ^/ ^+ d5 U9 F5 `& }% B=input('请输入由支路参数形成的矩阵:B=');
/ k% u$ `2 X' ]5 D! o% A=input('请输入各节点参数形成的矩阵:A=');
clear all
n=5;
/ F) W6 O) p, G2 A' Onl=5;
  k: h0 s7 S2 |- Y9 ]+ E7 @! f5 ^/ Gph=5;
. \& ?: F# l' C% x$ d" t* h8 U9 ]9 j* Hjd=1E-6;
3 J5 @" w  c( u( z( V' h%支路参数矩阵
Br=[1 2 0.04+0.25*j 0.25j 1;
   1 3 0.1+0.35*j  0     1;
1 ]8 x, O% ]' G! c& X   2 3 0.08+0.30*j 0.25*j 1;
$ {) w( t$ T' f' Y: f   2 4  0.015*j    0     1.05;
5 n2 F% Q& ^  q) z- ~- B; O( _% h   3 5  0.03*j     0     1.05];
%节点参数矩阵:2-PQ节点，3-PV节点，4-平衡节点
) D3 w/ o6 y! ZA=[-1.6-0.8*j     1     0    2;
' z& F6 Y/ `. l5 E9 c. h; s$ X    -2-1*j        1     0    2;
2 @" j" x* E! s8 s9 k% {( }9 k* f    -3.7-1.3*j    1     0    2;
0 |1 A+ p8 H* J" N0 y    5            1.05  1.05  3;
' g: R* d" u; D8 U1 [! T/ d# \: C; E; I    0            1.05  0     4];
, d7 k& z, K- o& }8 \% s
%雅克比矩阵形成
6 ]6 S) ]- s* M0 D. }Y=zeros(n);
: Y' X4 s/ l8 O7 g& I% s, E* ^e=zeros(1,n);
f=zeros(1,n);
V=zeros(1,n);
$ b) I! `% C& q) l% r" ~! Jfor i=1:nl
" c- a, I- s( h  A) B3 l    p=Br(i,1);
    q=Br(i,2);
    Y(p,q)=Y(p,q)-1./(Br(i,3)*Br(i,5));
* ?$ ~8 X/ L. I: }5 c    Y(q,p)=Y(p,q);
    Y(q,q)=Y(q,q)+1./(Br(i,3))+Br(i,4);
    Y(p,p)=Y(p,p)+1./(Br(i,3)*Br(i,5)^2)+Br(i,4);
end         
! F, Z6 f8 H2 @5 Edisp('节点导纳矩阵为：Y=');
- z# w$ z9 v( rdisp(Y);
1 t0 p+ J! {9 w5 a* w6 pG=real(Y);
) f' \+ s, F' E# Z4 jB=imag(Y);
for i=1:n
* N7 E3 K& `' E% h2 b$ _    e(i)=real(A(i,2));
+ M2 s$ j1 j% {* S    f(i)=imag(A(i,2));
( T; T8 d0 P' G& z" c9 b0 d8 @    V(i)=A(i,3);
2 T* t4 c8 j" o; p" _end
for i=1:n
    S(i)=A(i);
end
P=real(S);
4 Q: }# j* ^: {+ H# }' GQ=imag(S);
/ Y% Y1 l; A( y3 p$ C%雅克比矩阵求取
Ci=0;
8 c4 N4 ^$ @$ H& G% s* pa=1;
NO=2*n;
# W2 }% A$ d; U; p2 ON=NO-1;
while a~=0
  P" V) A8 P# }  H1 b    a=0;
) U6 Y1 A& S, f0 q" f    for i=1:n
! D( `+ r1 d4 f: E" r: A" B' ~0 J        if i~=ph
& {4 [# s3 Z7 _, n. I: O+ _0 h- d            C(i)=0;
            D(i)=0;
            for p=1:n
                C(i)=C(i)+G(i,p)*e(p)-B(i,p)*f(p);
  B/ C. @" ^$ k9 a8 ~. I/ N                D(i)=D(i)+G(i,p)*f(p)+B(i,p)*e(p);
8 Z% F0 M" F* ^/ A2 b% _            end
            P1=e(i)*C(i)+f(i)*D(i);
            Q1=f(i)*C(i)-e(i)*D(i);
            V2=e(i)^2+f(i)^2;
            if A(i,4)~=3
                DP=P(i)-P1;
' j. `6 I' ^* n2 u0 a9 \                DQ=Q(i)-Q1;
1 o# a# R8 I. t( R, `% b  c% d# Q                for k=1:n
6 }7 P& I, M$ L1 {% D8 @$ d+ y%非平衡节点时，非对角线元素
( @: I; M' n4 b) ]                    if k~=ph & k~=i                       
" h+ L( a- ~9 p$ }2 }                        X1=-G(i,k)*f(i)+B(i,k)*e(i);
                        X2=-G(i,k)*e(i)-B(i,k)*f(i);
                        X3=-X2;
                        X4=X1;
: t( c3 N( o- }" |- K- {1 F                        p=2*i-1;
                        q=2*k-1;
# v5 e( n- C* }9 ]                        J(p,q)=X1;
                        J(p,N)=DP;
                        m=p+1;
: N+ e, E6 V$ r                        J(m,q)=X3;
. Q2 d& s7 j6 E0 u4 H                        J(m,N)=DQ;
9 h$ i& c/ T, }% \" @! ?                        q=q+1;
                        J(p,q)=X2;
                        J(m,q)=X4;                       
%非平衡节点时对角线元素
% ?2 n8 \; U- u9 H                    elseif k~=ph & k==i
                            X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);
                            X2=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);                           
                            X3=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);
1 J. `  h2 C4 I                            X4=D(i)+B(i,i)*e(i)-G(i,i)*f(i);
                            p=2*i-1;
                            q=2*k-1;
                            J(p,q)=X1;
                            J(p,N)=DP;
3 c% E0 h$ \9 \, G$ ?0 Y8 u8 [% F                            m=p+1;
* }: u8 C: V  i' W9 s8 O) F                            J(m,q)=X3;
                            J(m,N)=DQ;
                            q=q+1;
4 O5 w: k+ m' T9 k8 n                            J(p,q)=X2;
8 R& J  A3 f, e2 a                            J(m,q)=X4;
                        end
& z- E3 h0 Z8 E) o                    end
4 d* Z7 i, g) I. Q/ m$ D                    else
                        DP=P(i)-P1;
                        DV=V(i)^2-V2;
% R# v5 {( K1 q. a9 e+ D8 b6 Y                        for s=1:n
% L8 L7 h0 {  a2 t5 @5 w; x4 U; U                            if s~=ph & s~=i
                                X1=-G(i,s)*f(i)+B(i,s)*e(i);
                                X2=-G(i,s)*e(i)-B(i,s)*f(i);                           
                                X5=0;
                                X6=0;
                                p=2*i-1;
                                q=2*s-1;
                                J(p,q)=X1;
# o7 R2 r6 ^8 I+ L4 T                                J(p,N)=DP;
                                m=p+1;
4 J# L. D, ^, L                                J(m,q)=X5;
                                J(m,N)=DV;
                                q=q+1;
                                J(p,q)=X2;
                                J(m,q)=X6;
) ~$ p1 i& m9 t1 |) N1 C                            elseif s~=ph & s==i
                                    X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);
0 v( P2 R: _5 {+ E+ \8 f                                    X2=-C(i)-G(i,s)*e(i)-B(i,s)*f(i);
                                    X5=-2*f(i);
                                    X6=-2*e(i);                                 
8 o1 `2 T; h+ V9 {! a9 g! l3 t4 z+ f                                    p=2*i-1;
& Y. @6 P2 i2 H- j/ f                                    q=2*s-1;
" F9 R; Q- x; l# z                                    J(p,q)=X1;
                                    J(p,N)=DP;
$ }/ `  K, M4 g  w* s                                    m=p+1;
                                    J(m,q)=X5;
                                    J(m,N)=DV;
& ~& u: B6 U  |0 [                                    q=q+1;
  ?8 A5 ^# w2 X( z: Q5 N8 F* q( ]                                    J(p,q)=X2;
                                    J(m,q)=X6;
                                end
8 R1 |4 I1 ]9 e4 j2 E                            end
) W  ^" _0 J; N: W% V5 s. q' [                        end
                end
            end
     disp('系统雅克比矩阵为：J=');
     disp(J);
%对于高压电力系统，最大元素出现在对角线位置，简单计算如下
     DJ=J(:,[1:NO-2]);
" j4 |6 n& P% J5 ~& @4 X     DP=J([1:NO-2],N);
9 j& A" {2 W+ G1 s     DU=DJ\DP;
: p5 ]4 v7 @% r     disp('************************************');
     disp('     ')
1 s4 `' I1 a3 N9 J* v     disp('第几次的修正值DU？');
     disp(Ci);
     disp('***********************************');
     disp(DU);
     for i=1:NO-2
         eps=abs(DP(i));
         if eps>=jd;
             a=a+1;
         end
2 A4 u$ S  a. ]7 w/ `% ~- l- g0 d6 @     end
     Ci=Ci+1;
     for i=1:n-1
6 |( K0 s. H& e) }* H& K         e(i)=e(i)-DU(2*i);
! z1 a2 Z1 U7 X, [) d# H% z( L         f(i)=f(i)-DU(2*i-1);
. H" q$ u8 ?: _! [     end
8 C8 _4 T0 p, m     E=e+j*f;
     disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
     disp('    ');
% |- \  K% \& c7 X, F8 [6 z     disp('第几次迭代后的近似计算值');
     disp(Ci);
" y( s- ~9 m. C; r/ o     disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
" t) l( r) l" I( E9 Y. P5 z- D     for i1=1:n
. I0 j( C* G7 P4 X+ v9 ^8 u         V(i1)=sqrt(e(i1)^2+f(i1)^2);
         ANGLE(i1)=atan(f(i1)./e(i1))*(180./pi);
     end
" R) m' y! u+ H! `7 B' x! }%下列输出结果均按节点号从小到大排列
1 Y! |, c' ^/ u7 _' c     disp('各个节点的电压实际值为：E=');
: _2 G, W( @7 I; p5 B/ x     disp(E);
$ m( D1 G4 c/ S8 U' i     disp('各个节点电压幅值为：V=');
& R2 l" K& K" y2 a) r4 C9 p8 r4 I     disp(V);
- f8 a5 v4 t7 X3 \; f     disp('各个节点的电压相交为：ANGLE=');
     disp(ANGLE);
; |" T! f* z- L     disp('各个节点的功率误差变化');
& q" E6 b+ b) N/ v% N7 h     disp(J([1:8],9)');
$ q; V" U8 B( `4 z6 c6 V     MDU(Ci)=max(abs(J([1:8],9)'));
     for p=1:n
8 m9 P; t1 p$ W6 c; W         Sm(p)=0;
         for q=1:n
             Sm(p)=Sm(p)+conj(Y(p,q))*conj(E(q));
' |0 y4 X* G' s/ [  k  A% z! r+ j         end
         S(i)=E(p)*Sm(p);
9 c0 V; \  W( U) M, A+ {     end
end
+ n* ^/ M* W( c# h$ C: fdisp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');
disp(' ');
disp('计算的最终功率近似值如下各项值');
& \" J( q6 I4 G, U# z7 Adisp(' ');
disp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');
; t) ~. Y8 Q3 y* u/ Ndisp('各个节点功率为:S=');
disp(S);
disp('各条支路的首端功率Si为:');
  Y. W# t6 D7 Yfor i=1:nl
3 v- F: L/ p! A' r0 Q    p=Br(i,1);q=Br(i,2);
    A(i)=E(p)*[conj(E(p))-conj(E(q))]*conj(1./(Br(i,3)));
3 ^9 [' Q3 p6 Y- f, n, Y; s    SBS(i)=[abs(E(p))]^2*conj(Br(i,4))+A(i);      
end
disp(SBS);
  J: t5 b, l7 {disp('各条支路的末端功率Sj为:');
for i=1:nl
    p=Br(i,1);q=Br(i,2);
9 o: M/ @- C8 y1 _    A2=E(q)*[conj(E(q))-conj(E(p))]*conj(1./(Br(i,3)));
* ?" E) N5 ~% G9 l* ]8 Q    SBE(i)=[abs(E(q))]^2*conj(Br(i,4))+A2;
* K2 F- D2 C) e- h; q+ fend
2 Q$ s+ g9 b, G  Gdisp(SBE);
disp('各条支路的功率损耗为:');
! T1 k' K2 I+ T* ofor i=1:nl
   p=Br(i,1);q=Br(i,2);
: h& E$ e0 V" c5 {& J" Z9 w   SBL(i)=SBS(i)+SBE(i);
6 M" |0 u& X- B5 y" g1 ^end
disp(SBL);
%绘制功率误差曲线，先用最小二乘法进行曲线拟合
2 f+ m" `, g4 S9 s3 _for i=1:Ci
    CSH(i)=i;
: h1 c/ K8 V  R- \1 ~1 H9 Dend
P=polyfit(CSH,MDU,3);
plot(CSH,MDU,'-O');

tj612

谢谢楼主分享！！！！！

201314

哦，这么发源程序似乎很累人哦！

^_^

参见 张伯明老师的  高等电力网络分析~~~~~~~~~~~~~

windlot

这么发源码倒是省了下载附件的钱了,呵呵

liusen19880123

楼主这样发源代码的方法值得大家学习啊

jimxi_6

程序可以运行？
??? Undefined function or variable 'm'.

sailor0305

请问楼主有没有关于人工免疫算法最有潮流的程序啊！急切盼望中.....

沧山牧水

上面的程序可以运行