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

沧山牧水

首先声明一下，这些是从网站上转载的，不是本人上编写的
# q" I4 y2 S  ?. W$ q! G6 G+ T% k基于MATLAB的电力系统潮流计算

%简单潮流计算的小程序，相关的原始数据数据数据输入格式如下：
8 D4 g% \' R/ Q$ l, V2 M1 U: a( T%B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写
& a# H* x( O2 y2 l%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点
%编号，将变压器的串联阻抗置于低压侧处理。
%第三列为支路的串列阻抗参数。
%第四列为支路的对地导纳参数。
7 r# E! J1 d/ J& q%第五烈为含变压器支路的变压器的变比
% a6 A; p2 M" Y. s3 u# B%第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，
5 N  u% O% {* [& n6 p5 k5 S%“0”为不含有变压器。
%B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点
%负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中
%“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。
# G" a; h5 \) M
%X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地
! g% J, U. n7 l# g( d% t; c%参数。
- r7 i8 K$ C5 }- Q: I" ]# a5 {n=input('请输入节点数:n=');
n1=input('请输入支路数:n1=');
' }! Y% d% P& Q" sisb=input('请输入平衡节点号:isb=');
pr=input('请输入误差精度:pr=');
B1=input('请输入支路参数:B1=');
$ Y3 z: W) H# B6 {! JB2=input('请输入节点参数:B2=');
. w" H  a6 K% q# [X=input('节点号和对地参数:X=');
Y=zeros(n);
. V) n7 R% K4 i! I" ?6 eTimes=1;                %置迭代次数为初始值
  X- v8 P/ Z, p" ~6 O+ J%创建节点导纳矩阵
for i=1:n1
    if B1(i,6)==0       %不含变压器的支路
( Q  Q6 K6 m: ?6 o( U        p=B1(i,1);
% h& D+ F* }, X8 [, n+ C        q=B1(i,2);
        Y(p,q)=Y(p,q)-1/B1(i,3);
        Y(q,p)=Y(p,q);
1 L% E+ e4 v+ k0 X        Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);
' Z9 {' i) V5 X        Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);
1 d" O! _7 H( [& j6 l    else                %含有变压器的支路
        p=B1(i,1);
& v- r- c) Y4 ?3 w+ r' Y) q% k        q=B1(i,2);
        Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5));
  \& E' C# Y7 h& p        Y(q,p)=Y(p,q);
        Y(p,p)=Y(p,p)+1/B1(i,3);
  f, C6 A3 d! C/ z        Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));
, _* @! `0 u! P' @& n* Q: i    end
end
Y
OrgS=zeros(2*n-2,1);
DetaS=zeros(2*n-2,1);   %将OrgS、DetaS初始化
& m  j1 a7 O, _; t%创建OrgS，用于存储初始功率参数
h=0;
4 f8 T+ C  C8 v4 Gj=0;
for i=1:n            %对PQ节点的处理
9 @7 q$ y9 P& X8 h    if i~=isb&B2(i,6)==2   
        h=h+1;
5 z8 s0 O  L: v5 O/ I. A+ y& ]5 F3 ~        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)));
2 I3 u2 V5 e7 S- V" g            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)));
2 S  `# \# s+ {        end
    end
end
for i=1:n           %对PV节点的处理，注意这时不可再将h初始化为0
    if i~=isb&B2(i,6)==3
2 I% g% G  C7 s1 x        h=h+1;
        for j=1:n
1 Q/ x  a0 E4 z% \  v+ @# `            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)));
9 o) c$ r6 g3 j3 j5 Y            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
end
+ `5 U: P, t* p' ]OrgS
%创建PVU 用于存储PV节点的初始电压
PVU=zeros(n-h-1,1);
t=0;
( o2 O. }4 V- N' P6 Zfor i=1:n
    if B2(i,6)==3
        t=t+1;
        PVU(t,1)=B2(i,3);
* C! a6 q- C5 B9 p, P! G    end
2 v: L0 F1 f* `- M3 B& |4 H/ nend
PVU
& d2 Z, ?% G2 x* F% y( i, P! s%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量
1 V/ }1 o1 I9 O1 ?4 y% Uh=0;
& Z5 I& G. s- K8 E8 d6 g6 n* ufor i=1:n           %对PQ节点的处理
    if i~=isb&B2(i,6)==2
, W" C6 c' z& R9 Y% T1 U: D        h=h+1;
        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);
  T% ]/ r9 p3 f. F: B$ X! v9 v! R& A    end
( O5 e4 Z( l6 `% I8 O& ]. v5 Kend
t=0;
for i=1:n           %对PV节点的处理，注意这时不可再将h初始化为0
    if i~=isb&B2(i,6)==3
        h=h+1;
        t=t+1;
- q; I3 ^, m. q) J2 M' B" ~        DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
4 }& {' V/ B& N) R7 E        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
DetaS
%创建I，用于存储节点电流参数
i=zeros(n-1,1);
h=0;
for i=1:n
    if i~=isb
        h=h+1;
        I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));
    end
end
I
3 [9 L9 }9 G; J) Y%创建Jacbi(雅可比矩阵)
Jacbi=zeros(2*n-2);
. n( V! i  ?; Uh=0;
0 K5 t4 a6 Z0 A& F+ w. L4 \* ]k=0;
for i=1:n       %对PQ节点的处理
7 b" J: h) k% s    if B2(i,6)==2
: k5 Y6 j, y9 r& |  P$ y9 ?        h=h+1;
        for j=1:n
' `0 Y' u/ _0 B1 s            if j~=isb
( p" J* ^: S7 A" j1 Z1 L  v' w% Y                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));
3 Z2 x& Z% V. Y' f4 O) P                    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));
9 `1 J' ?" @, h- y3 k! m                else        %非对角元素的处理
0 I' n  m* s1 R2 f$ @5 V. t                    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));
# u# C- p, q  ?' [/ L8 [$ v                    Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);
                    Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);
                end
; z6 q  l2 x* f, B, I8 y6 p                if k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行
                    k=0;
                end
& F- e& S8 g3 w. g( f2 t+ `7 @            end
; C5 Y- j2 |# a+ Y, T        end
2 s7 I8 s/ {+ f+ r$ W6 m6 @    end
end
* e/ ~; A3 S  w& l6 ok=0;
# z% B+ S1 \1 K% Y0 G( bfor i=1:n       %对PV节点的处理
    if B2(i,6)==3
# P, Z/ h( U' |" L6 H: b        h=h+1;
        for j=1:n
  K  ]/ o0 a* z) b! X8 T: F* T            if j~=isb
) O4 ^2 {( [- i% J3 N                k=k+1;
7 x- f) o, [4 P                if i==j     %对角元素的处理
- s- g) d; @: q% X0 o                    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));
; T7 V: I: S) v6 e7 d# y# f                    Jacbi(2*h,2*k-1)=2*imag(B2(i,3));
! O4 Q* |- q0 j                    Jacbi(2*h,2*k)=2*real(B2(i,3));
                else        %非对角元素的处理
2 A# [/ t: L: k0 M( w, ]# I                    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));
& T1 K9 G% E9 l                    Jacbi(2*h,2*k-1)=0;
( L7 z7 r( s) g: Z6 E* s                    Jacbi(2*h,2*k)=0;
                end
                if k==(n-1)     %将用于内循环的指针置于初始值，以确保雅可比矩阵换行
                    k=0;
4 Z" X/ V3 @1 i. T" Q1 _                end
            end
6 b5 r  E# H7 K" X, c' g        end
    end
' K8 [/ w3 Q. c) N4 B4 p. e+ Pend
8 U5 h6 Q5 l+ c. j7 r4 jJacbi
2 G; m2 Q1 K2 G# Q' l%求解修正方程，获取节点电压的不平衡量
DetaU=zeros(2*n-2,1);
DetaU=inv(Jacbi)*DetaS;
* M; r7 l+ X& `DetaU
, T) g7 w' Y; G* x6 s%修正节点电压
j=0;
. X( U. v" ]- X% |3 ~/ u  w# qfor i=1:n       %对PQ节点处理
    if B2(i,6)==2
) F- y  Q8 i7 z        j=j+1;
& p* V2 I( d5 ?! b2 A7 |0 p2 `& ^        B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
    end
end
1 Z& q4 p4 c) g* k( w8 p, @4 tfor i=1:n       %对PV节点的处理
8 ^8 @$ J! @3 A$ p    if B2(i,6)==3
# P9 T, V) M5 ]6 `. S6 i        j=j+1;
$ o  {5 T! t! Q. ?. u        B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
1 R, u8 c2 O$ g* ^$ Y& \* r    end
end
B2
%开始循环**********************************************************************
7 C2 i9 Z) Z1 g, V/ A: F3 L* d' Qwhile abs(max(DetaU))>pr
/ J) e/ M+ f, eOrgS=zeros(2*n-2,1);        %!!!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵
1 j( \, t7 j" q  qh=0;
* ]" {! ?; J7 s7 lj=0;
for i=1:n
    if i~=isb&B2(i,6)==2
        h=h+1;
% q3 N/ G# g8 i1 D' w3 |$ x        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)));
        end
    end
end
for i=1:n
1 `: ?" l: d' Q: }: X1 l% P    if i~=isb&B2(i,6)==3
        h=h+1;
        for j=1:n
2 j/ g' G6 U6 ?2 K2 ?3 E            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)));
: G  k# n" |: O+ e$ h            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
# T9 U- W5 R& H  z- M1 p    end
end
OrgS
%创建DetaS
. ]1 I  y/ O9 |0 ih=0;
7 |) x5 M; I( B1 H# v* y0 o8 Y) n: jfor i=1:n
    if i~=isb&B2(i,6)==2
        h=h+1;
5 j( w, M% {- ^# c5 U, x        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);
    end
end
t=0;
: W, p! t/ S' K% q# r' \6 ?4 Pfor i=1:n
    if i~=isb&B2(i,6)==3
        h=h+1;
        t=t+1;
3 R# D0 I. S  N4 ]* a. ^6 q% f, n        DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
/ e2 m. d2 \; F7 W! {$ J3 }        DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;
, b/ w  ^& T$ s/ X' _1 l    end
end
DetaS
. ~8 v8 q: l: n# O%创建I
- V# \3 z' A- J* d; T  F9 |i=zeros(n-1,1);
7 F$ ~; [, |* D/ i$ mh=0;
1 `# i$ f" h* mfor i=1:n
- z- E1 R$ z; z8 V( H9 B    if i~=isb
( N6 j' H" R/ i: i* J0 W        h=h+1;
        I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));
7 n7 D1 H; K, X5 L8 J! L9 ~$ j, i. e    end
! e/ h% Y) T% ^- Vend
I
6 |+ I% ]" r% U3 d( S# P3 [" F$ v%创建Jacbi
+ f% q$ G* c: v) [, F& g  [; iJacbi=zeros(2*n-2);
h=0;
k=0;
for i=1:n
    if B2(i,6)==2
        h=h+1;
        for j=1:n
8 U) C! C) }* b; O; H            if j~=isb
                k=k+1;
                if i==j
5 n5 e: S5 L' b, H2 L                    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));
; h, t  V0 b9 N/ N4 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));
! ?% [+ S7 w9 h4 S3 r. }                    Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));
* o& E4 ]! d) J* r) L; l- P                    Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));
                else
6 {5 M) c! K! K) M# V3 T                    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));
& `% X" O; J* Z7 |+ x9 T                    Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);
                    Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);
( P5 L" W" y; Z' `                end
5 m$ J9 j+ b6 w  Z! A' E5 I                if k==(n-1)
                    k=0;
                end
( ^' n' U' z6 H. B5 @& R            end
        end
1 Y6 B( Y- _$ b, ?% ?  [1 W; k    end
* ^( d8 X9 `+ S7 p; Wend
k=0;
for i=1:n
    if B2(i,6)==3
: K0 j3 p' U8 D- i" n* D) u        h=h+1;
, L" v, j  x& f: F        for j=1:n
            if j~=isb
$ j/ l9 z  N* _8 N                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));
: Q) J- ?. y! r8 U                    Jacbi(2*h,2*k-1)=2*imag(B2(i,3));
                    Jacbi(2*h,2*k)=2*real(B2(i,3));
8 S5 x( N- }; p. z: 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)=0;
5 H1 J. b, F. b                    Jacbi(2*h,2*k)=0;
) E% k" M; h: R0 b8 |1 u# }                end
; u  F6 H( g+ F+ i                if k==(n-1)
$ \* h8 _" i; ^3 ~: b                    k=0;
1 v+ n. w9 L7 E; `/ d  U                end
* c; d! n6 D+ E$ e. c' l( F% ?            end
        end
    end
, ?& r$ N2 q4 s9 Nend
Jacbi
DetaU=zeros(2*n-2,1);
DetaU=inv(Jacbi)*DetaS;
! e+ I1 e$ C# p3 d$ ~) jDetaU
%修正节点电压
j=0;
for i=1:n
( a' w* G+ `6 r$ h8 K( {  P4 g1 `    if B2(i,6)==2
        j=j+1;
        B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
    end
- n0 v0 j1 U4 u4 u' P9 aend
for i=1:n
" P6 w! D; n0 h; P7 f0 d8 X    if B2(i,6)==3
        j=j+1;
        B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
. d% m( N! Z% O# u' ~  `    end
! \' l; _# g& ?1 h0 ]  u7 L& xend
5 I( I4 @* U* X8 RB2
Times=Times+1;      %迭代次数加1
end
, ^5 _( Q6 `; S: ~: |9 A$ |" P0 ~Times
, H0 }' W9 ~& x+ k
9 b1 M1 y5 r7 ?/ d0 b一个原始数据的例子
! r8 c8 {' \/ f9 d节点数 5
# W7 r3 M! T" n支路数 5
2 r0 n* e9 X; s8 }, C平衡节点编号 5
7 W3 F% A2 |6 M' _: b8 V精度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]
; d* E* G, o- x5 Y6 j% ?B2(节点参数矩阵)
[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(节点号和对地参数)
/ g; R1 W, V; v9 e( L% r[1 0;2 0;3 0;4 0;5 0]

沧山牧水

程序2：
3 d" n; q8 C! g6 b' `6 e0 u7 u* U5 I5 Z5 p%本程序采用牛顿-拉夫讯法对一个五节点电力系统进行潮流计算
% e8 Z5 O9 |. h5 I/ d% d% T! B7 c% n=input('请输入系统节点数目:n=');
" d8 }7 M6 |, ^: Y" I2 @% nl=input('请输入系统支路数目:nl=');
% ph=input('请输入平衡节点号:ph=');
% jd=input('请输入误差精度：jd=');
% B=input('请输入由支路参数形成的矩阵:B=');
% A=input('请输入各节点参数形成的矩阵:A=');
8 u- Y6 S' H) P# _clear all
0 S' C/ D0 I  g9 \0 o  s9 n: ?n=5;
nl=5;
ph=5;
9 I" F+ T2 f0 D, Cjd=1E-6;
%支路参数矩阵
Br=[1 2 0.04+0.25*j 0.25j 1;
   1 3 0.1+0.35*j  0     1;
   2 3 0.08+0.30*j 0.25*j 1;
  M3 R; |+ V7 f! y) f   2 4  0.015*j    0     1.05;
  A! x7 r3 ~/ Y7 Q; D4 o5 U   3 5  0.03*j     0     1.05];
%节点参数矩阵:2-PQ节点，3-PV节点，4-平衡节点
& n/ i! R3 x/ {6 iA=[-1.6-0.8*j     1     0    2;
7 l' j' S! {% \/ i9 c6 @1 s    -2-1*j        1     0    2;
; L! t8 @, F! v$ H5 Z  {# ~    -3.7-1.3*j    1     0    2;
2 b8 _6 o8 h1 x' ?  W' t2 h! u/ S    5            1.05  1.05  3;
6 p% A! k4 @& T+ `3 v    0            1.05  0     4];

%雅克比矩阵形成
Y=zeros(n);
e=zeros(1,n);
f=zeros(1,n);
V=zeros(1,n);
for i=1:nl
5 j+ N7 t% ~! ]    p=Br(i,1);
! v6 ~3 v; k7 O    q=Br(i,2);
! @( M5 l- `$ o" p    Y(p,q)=Y(p,q)-1./(Br(i,3)*Br(i,5));
    Y(q,p)=Y(p,q);
- g' A2 Y8 |) y2 k, ^- x    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         
disp('节点导纳矩阵为：Y=');
- P; Y2 ~* z7 q( j8 [disp(Y);
G=real(Y);
B=imag(Y);
for i=1:n
: K9 a9 ^/ p: F! k+ B2 j    e(i)=real(A(i,2));
    f(i)=imag(A(i,2));
( r% v8 [6 ?, J0 w% [" ^, h3 m/ T1 m    V(i)=A(i,3);
  F% T& h; g* P) gend
5 F8 }# y& c, r7 p( J0 l; cfor i=1:n
% i  P3 L4 d& m9 D! a2 @" }) @    S(i)=A(i);
  K. C% }$ [1 Mend
P=real(S);
+ b# J- B) P" ^+ ^* U7 ]Q=imag(S);
%雅克比矩阵求取
" S- z  [. m' Y4 u% u/ z& tCi=0;
a=1;
* d9 N6 n  U4 y) ^7 U* zNO=2*n;
" E7 i* e$ N" }/ V  cN=NO-1;
while a~=0
$ `- ?) b; A$ ^/ ^2 _    a=0;
    for i=1:n
2 u/ r1 @  B/ o- W( \        if i~=ph
            C(i)=0;
            D(i)=0;
  i. C4 f, D. r            for p=1:n
                C(i)=C(i)+G(i,p)*e(p)-B(i,p)*f(p);
                D(i)=D(i)+G(i,p)*f(p)+B(i,p)*e(p);
            end
) }, O5 Y  y) }! a4 [5 m            P1=e(i)*C(i)+f(i)*D(i);
6 [' S/ p1 S2 l1 ?; }, u3 Q6 |4 T            Q1=f(i)*C(i)-e(i)*D(i);
            V2=e(i)^2+f(i)^2;
( c  D, ^" j% Y/ u9 Q            if A(i,4)~=3
                DP=P(i)-P1;
                DQ=Q(i)-Q1;
                for k=1:n
" r; V: Z' c, h: R" K4 D% Q% D%非平衡节点时，非对角线元素
                    if k~=ph & k~=i                       
  x! Y2 O, l5 U6 i. h1 {+ ^+ W. B                        X1=-G(i,k)*f(i)+B(i,k)*e(i);
" A% F- t. F" H; e7 s                        X2=-G(i,k)*e(i)-B(i,k)*f(i);
                        X3=-X2;
                        X4=X1;
                        p=2*i-1;
+ {# r( \# e6 t1 a( T! U0 J                        q=2*k-1;
. b$ S. _2 w, ]# x* S                        J(p,q)=X1;
4 C2 J1 u; Q5 @! N- D3 p                        J(p,N)=DP;
* G6 N1 p7 s* ]4 x# W0 K6 T                        m=p+1;
& d5 k# w3 t; q) W                        J(m,q)=X3;
                        J(m,N)=DQ;
                        q=q+1;
                        J(p,q)=X2;
+ Z  Z7 i2 S" g' j6 k; R                        J(m,q)=X4;                       
: L- E' p" p# A9 m8 R%非平衡节点时对角线元素
5 E4 D. f" D/ D# ^8 [+ C                    elseif k~=ph & k==i
                            X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);
6 x. r% i9 q8 z/ I4 |) Z, t                            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);
                            X4=D(i)+B(i,i)*e(i)-G(i,i)*f(i);
' ^( \, o0 o, ^; y. q4 e9 h1 \                            p=2*i-1;
8 u; @& D! g) J' O8 |4 b# A                            q=2*k-1;
                            J(p,q)=X1;
  d4 g" v' Y7 R2 G6 B, R, V                            J(p,N)=DP;
2 \) I3 c1 j  W: R- N" Q                            m=p+1;
2 ~- j# `" x7 t- q  \0 x" J# N                            J(m,q)=X3;
                            J(m,N)=DQ;
                            q=q+1;
                            J(p,q)=X2;
7 J" B: x3 U2 l: A& E                            J(m,q)=X4;
                        end
                    end
+ [, @) [5 |0 q6 U" ?" Z                    else
                        DP=P(i)-P1;
( H/ O$ u0 n  e1 U; W                        DV=V(i)^2-V2;
( L/ E* U: u! S1 A' F  f" d                        for s=1:n
                            if s~=ph & s~=i
                                X1=-G(i,s)*f(i)+B(i,s)*e(i);
" j( r$ A' i$ r                                X2=-G(i,s)*e(i)-B(i,s)*f(i);                           
                                X5=0;
& B- K" ~9 O$ K- L( N& S; E                                X6=0;
8 I7 Z  s1 `4 R" n' u                                p=2*i-1;
                                q=2*s-1;
7 K0 }( b( g, Y; v: ?                                J(p,q)=X1;
* @, E' u6 R' k                                J(p,N)=DP;
                                m=p+1;
" |9 f6 ~" H" \                                J(m,q)=X5;
                                J(m,N)=DV;
                                q=q+1;
" E6 U) w& f$ \' K; F4 s' W$ Q                                J(p,q)=X2;
$ X2 I0 Q5 f% g  }                                J(m,q)=X6;
                            elseif s~=ph & s==i
( O# z% @! X+ o                                    X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);
                                    X2=-C(i)-G(i,s)*e(i)-B(i,s)*f(i);
                                    X5=-2*f(i);
' S* q+ ?0 J, j% l% ~& ~                                    X6=-2*e(i);                                 
                                    p=2*i-1;
# M. F6 K. h( _8 t4 F                                    q=2*s-1;
                                    J(p,q)=X1;
                                    J(p,N)=DP;
                                    m=p+1;
                                    J(m,q)=X5;
                                    J(m,N)=DV;
                                    q=q+1;
% a( K/ I+ D/ k. R7 v: U$ u                                    J(p,q)=X2;
) ?' N2 r5 P/ I* O; v                                    J(m,q)=X6;
                                end
                            end
2 a" N! S2 e. Z4 M                        end
                end
: k9 y5 y' A$ e/ s            end
9 g; x$ F4 g8 G. y% }* C     disp('系统雅克比矩阵为：J=');
4 M8 e8 E4 K1 d     disp(J);
# |* ^  b5 A/ D%对于高压电力系统，最大元素出现在对角线位置，简单计算如下
     DJ=J(:,[1:NO-2]);
     DP=J([1:NO-2],N);
. `8 {$ _7 ?, T/ i% a     DU=DJ\DP;
1 n2 L1 P7 [( R! b7 U: u: e1 a( y     disp('************************************');
4 D" g9 m, f! X& x, z7 T/ s     disp('     ')
     disp('第几次的修正值DU？');
1 z5 o) ~/ M% ~' r& z+ b     disp(Ci);
$ G1 u; f% C: h- t" E/ p; E     disp('***********************************');
: O; ]9 i  V; z/ K8 `: k/ k( {     disp(DU);
     for i=1:NO-2
         eps=abs(DP(i));
2 Q# B0 H8 B. y0 T1 l" T         if eps>=jd;
             a=a+1;
         end
     end
- ]% W* U  r! N' o8 W" T' V     Ci=Ci+1;
     for i=1:n-1
: N; q- Y( W* x' d         e(i)=e(i)-DU(2*i);
  [7 ^, @# U4 I         f(i)=f(i)-DU(2*i-1);
. y3 z; [% E! F     end
     E=e+j*f;
1 J5 M$ B( m7 J4 R1 J* h     disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
     disp('    ');
1 B3 }8 m% t1 P+ I* [! w8 s     disp('第几次迭代后的近似计算值');
3 o$ b. ^. W8 }0 R% H     disp(Ci);
) V2 F7 Y, M/ Z     disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
2 p6 n. t- V! D5 |2 @     for i1=1:n
         V(i1)=sqrt(e(i1)^2+f(i1)^2);
9 ]; z" @5 j. ?: q& v0 m         ANGLE(i1)=atan(f(i1)./e(i1))*(180./pi);
     end
%下列输出结果均按节点号从小到大排列
/ g: L$ w/ O4 u2 M; o! e     disp('各个节点的电压实际值为：E=');
1 I% n9 p  i+ E: T' d* B# I     disp(E);
! G& T* V' `2 z" f* K. u     disp('各个节点电压幅值为：V=');
     disp(V);
     disp('各个节点的电压相交为：ANGLE=');
, S3 j- l' L* E; G9 l     disp(ANGLE);
     disp('各个节点的功率误差变化');
- d/ L+ p; W7 b3 Q7 i' V, I     disp(J([1:8],9)');
     MDU(Ci)=max(abs(J([1:8],9)'));
     for p=1:n
- F) B1 B$ V6 N  P         Sm(p)=0;
         for q=1:n
             Sm(p)=Sm(p)+conj(Y(p,q))*conj(E(q));
         end
         S(i)=E(p)*Sm(p);
  M  W5 F; S% g' J3 v     end
end
8 Y& ^) ~8 s/ O9 W2 g2 S, udisp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');
disp(' ');
disp('计算的最终功率近似值如下各项值');
disp(' ');
1 c" b3 P1 C( `  C/ }- q" Sdisp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');
& d% I  r7 ^7 S# c7 Jdisp('各个节点功率为:S=');
disp(S);
1 p9 i1 p" @+ b/ Y& q" l# [/ A' xdisp('各条支路的首端功率Si为:');
+ G8 Q* ~4 k4 Q' y; ?& Y4 Kfor i=1:nl
8 O- [: c' i% I6 \4 i    p=Br(i,1);q=Br(i,2);
    A(i)=E(p)*[conj(E(p))-conj(E(q))]*conj(1./(Br(i,3)));
    SBS(i)=[abs(E(p))]^2*conj(Br(i,4))+A(i);      
2 J( R" J2 J; Qend
disp(SBS);
disp('各条支路的末端功率Sj为:');
for i=1:nl
2 b) M( V, K+ u) i    p=Br(i,1);q=Br(i,2);
3 `; ?7 v, [" e, w5 c    A2=E(q)*[conj(E(q))-conj(E(p))]*conj(1./(Br(i,3)));
    SBE(i)=[abs(E(q))]^2*conj(Br(i,4))+A2;
end
disp(SBE);
+ |, V* K- g# \8 G# U; M$ x( idisp('各条支路的功率损耗为:');
for i=1:nl
2 ~- Y- g; ?5 r1 R   p=Br(i,1);q=Br(i,2);
   SBL(i)=SBS(i)+SBE(i);
end
& B1 ]7 \/ n, ^! Z8 C7 S6 mdisp(SBL);
%绘制功率误差曲线，先用最小二乘法进行曲线拟合
for i=1:Ci
    CSH(i)=i;
end
P=polyfit(CSH,MDU,3);
plot(CSH,MDU,'-O');

tj612

谢谢楼主分享！！！！！

201314

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

^_^

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

windlot

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

liusen19880123

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

jimxi_6

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

sailor0305

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

沧山牧水

上面的程序可以运行