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

沧山牧水

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

0 ^9 W+ I5 G* \7 L" \, O+ w%简单潮流计算的小程序，相关的原始数据数据数据输入格式如下：
%B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写
%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点
8 g6 k! q7 B- o. k/ I%编号，将变压器的串联阻抗置于低压侧处理。
%第三列为支路的串列阻抗参数。
  z3 O* k; a  q0 P+ \%第四列为支路的对地导纳参数。
9 P/ `" v! b# D2 F%第五烈为含变压器支路的变压器的变比
%第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，
%“0”为不含有变压器。
%B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点
%负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中
! ^9 C9 Z. V' l* w& U6 X%“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。
% [# s/ z2 O( L
) S( J* e' T5 g$ ^; \' K& @%X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地
+ `6 n  y2 i8 P6 b2 u3 g%参数。
# t" _6 F9 n2 i  On=input('请输入节点数:n=');
, [. v: g2 X7 ~, T% Mn1=input('请输入支路数:n1=');
isb=input('请输入平衡节点号:isb=');
pr=input('请输入误差精度:pr=');
B1=input('请输入支路参数:B1=');
' N' X2 E* c# Y& [: P) `6 FB2=input('请输入节点参数:B2=');
X=input('节点号和对地参数:X=');
Y=zeros(n);
3 g" B6 R7 H; \1 T: ]Times=1;                %置迭代次数为初始值
7 u( H: b( ]& D%创建节点导纳矩阵
for i=1:n1
    if B1(i,6)==0       %不含变压器的支路
        p=B1(i,1);
        q=B1(i,2);
2 q( {6 p! P  p6 Z        Y(p,q)=Y(p,q)-1/B1(i,3);
* g8 J6 Q3 t3 H2 s4 H        Y(q,p)=Y(p,q);
        Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);
        Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);
% s8 w' j) e0 T# g& v* a    else                %含有变压器的支路
        p=B1(i,1);
        q=B1(i,2);
1 `' O( ^% Y( N- s7 L* O* `% G0 j        Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5));
9 J+ Z& \9 T* W# o7 `        Y(q,p)=Y(p,q);
        Y(p,p)=Y(p,p)+1/B1(i,3);
" q( h: l# f" f        Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));
, p( S( i8 [3 X, ]    end
end
Y
OrgS=zeros(2*n-2,1);
. Y) I3 t" C% X6 u- _DetaS=zeros(2*n-2,1);   %将OrgS、DetaS初始化
%创建OrgS，用于存储初始功率参数
h=0;
j=0;
for i=1:n            %对PQ节点的处理
8 ~% Z6 p. {# v: ~4 G, M- j    if i~=isb&B2(i,6)==2   
: u3 y! u$ q4 V) O3 L; Q        h=h+1;
        for j=1:n
$ f1 T& c5 \7 L$ u            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           %对PV节点的处理，注意这时不可再将h初始化为0
& a6 u6 u- r$ m' h+ K; P- X    if i~=isb&B2(i,6)==3
( W5 c; _& k( l% x# x        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)));
$ u+ F# T3 n) ~            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 q8 z6 e" }. [        end
- }1 V5 ?! V  o0 G3 g3 A    end
end
OrgS
%创建PVU 用于存储PV节点的初始电压
8 x4 i- c6 |3 h$ P9 sPVU=zeros(n-h-1,1);
& {' R3 x: y, w. At=0;
) w3 l4 K( U" Y. sfor i=1:n
+ K7 y* ^$ B( F, M    if B2(i,6)==3
* }, P% p& _7 ]. q+ Z: h7 M# c* V        t=t+1;
: I- H  O. Q8 @! O9 Z& A        PVU(t,1)=B2(i,3);
    end
# U: ^2 F' D* Q( I" P" kend
" e& @8 ~/ x9 O3 y  r2 f, V7 dPVU
9 f+ U* i0 C8 L; l%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量
% o3 R6 \; L% Yh=0;
5 m2 D) t+ @( `for i=1:n           %对PQ节点的处理
    if i~=isb&B2(i,6)==2
        h=h+1;
        DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
7 e" s7 s& P- a/ }$ V0 X2 ~" i        DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);
    end
end
t=0;
# a: ~3 N% ?0 Y; T- @. c8 bfor i=1:n           %对PV节点的处理，注意这时不可再将h初始化为0
% S2 h0 |# p. u9 {3 {2 G* X( r8 @5 H( x    if i~=isb&B2(i,6)==3
' t, l% R  j, [6 f. j" Y8 A        h=h+1;
3 T2 a, ^& ~& u' @. c1 f        t=t+1;
' I1 w8 l3 y6 s1 \# P        DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
0 |7 ^( ]( d& C; d. S, U) L, G6 O        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
% m6 Y/ b) S, Z3 gend
6 d: _7 N; s/ P3 C* Z. i# j. xDetaS
) Z2 s/ [2 Q' Z# o% x- ^8 V%创建I，用于存储节点电流参数
, |5 i; X$ R0 wi=zeros(n-1,1);
( m( p+ [7 {1 g" B2 bh=0;
for i=1:n
    if i~=isb
& t8 `. W" z  o        h=h+1;
        I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));
* d" \! h! F& e9 t# t, N: t" n0 R    end
( E# K. ?: a4 M5 j, send
I
%创建Jacbi(雅可比矩阵)
Jacbi=zeros(2*n-2);
h=0;
6 k* f, v4 D5 S1 r' |7 A: a# C2 Xk=0;
' n7 I0 j# p( `for i=1:n       %对PQ节点的处理
# q1 m1 n+ {8 V" R6 W    if B2(i,6)==2
        h=h+1;
4 g, V* E2 D. S; u5 Y        for j=1:n
2 a$ U. T- @% h9 N; a            if j~=isb
8 @7 a; K3 ~& Z& v                k=k+1;
" Q5 u+ g8 t8 P7 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));
/ j$ z! k4 U# @                    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));
5 O9 x, Q% {! F1 a# h                    Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));
4 i. v" ~+ _9 o# o/ c                else        %非对角元素的处理
2 z# z3 E  N. g7 q                    Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));
9 F. x& }9 C* f$ T  L: j& x4 T                    Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));
# Q( d, c  a; j) m7 {                    Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);
: T' \# i4 k# D+ P5 E# [/ [                    Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);
                end
. b* R, o# b$ l                if k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行
' ^1 a2 z9 W5 O8 Y( o/ ^* k4 \7 ?                    k=0;
                end
            end
        end
    end
end
3 J% @  B1 c( t4 Z! [/ c* y' ek=0;
for i=1:n       %对PV节点的处理
+ A( ?8 s& a, ~' f6 B0 O3 \. G    if B2(i,6)==3
        h=h+1;
        for j=1:n
            if j~=isb
/ K- a+ X5 n8 g% B- G4 J                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));
' I! p3 \. n# A9 x                    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+ b& s6 p4 K1 v! j                    Jacbi(2*h,2*k-1)=2*imag(B2(i,3));
                    Jacbi(2*h,2*k)=2*real(B2(i,3));
9 _- M% ?) Y" ~, k- B- u                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;
                    Jacbi(2*h,2*k)=0;
                end
' V  o# |6 U$ I3 @                if k==(n-1)     %将用于内循环的指针置于初始值，以确保雅可比矩阵换行
                    k=0;
. N) C- P) a4 e0 P2 p( a                end
( S0 A, S! t* M- k6 o$ u8 W            end
: B( m% _" r+ j1 U+ Y& D        end
: I8 [- }- K* ^7 ~- r    end
end
2 |5 E: \8 D9 V: j9 s8 ]Jacbi
/ ]. `; c7 A6 @0 k; J. I%求解修正方程，获取节点电压的不平衡量
% h. v/ t: o' J) V  T; pDetaU=zeros(2*n-2,1);
( ?3 c' d% x: j8 @* DDetaU=inv(Jacbi)*DetaS;
DetaU
%修正节点电压
j=0;
for i=1:n       %对PQ节点处理
5 A' `. J$ U1 C4 }    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);
7 Z. O1 @3 m" Z    end
- g& C3 U" T6 b8 w/ i- Y; Vend
$ u0 E& y! A7 m  j" R  @6 Rfor i=1:n       %对PV节点的处理
    if B2(i,6)==3
! v& J$ F1 t# p- i! L; Z( t0 c        j=j+1;
        B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
7 Q- E) r) h2 u7 Y, |( h, Q' ^. K0 B- K    end
end
0 @  z, P  X. tB2
2 G6 ?3 \$ c! _+ X) Z! K" \%开始循环**********************************************************************
7 B9 t- v4 _$ Fwhile abs(max(DetaU))>pr
- w$ O' K$ A% `OrgS=zeros(2*n-2,1);        %!!!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵
h=0;
j=0;
8 W2 y+ K$ Q7 ffor i=1:n
    if i~=isb&B2(i,6)==2
        h=h+1;
- `4 R" ?* c5 ]+ t        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)));
, n. c+ Y0 S$ e: ~: z' j            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)));
# H- X# c' o' W8 S" a* h        end
- O# i' G: R7 Y7 \, i    end
4 }! W* F8 O& R6 Y5 E2 |7 ^end
( w) X  |& Z, m, C3 e- Ffor i=1:n
) \( F/ _9 ~5 m* L& u) t1 n9 }    if i~=isb&B2(i,6)==3
. e0 R' T# S& N7 T, j, D        h=h+1;
6 K" n. Q8 D# `  @5 z8 G        for j=1:n
4 g' h& }! M/ |# s            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)));
+ x. O% p' @& n# j        end
* B' W; U/ a" D; c3 \& R    end
# ?7 v/ U! T# J+ B/ a/ o# M3 D! Iend
OrgS
%创建DetaS
7 ]% f8 P- L( j) g! i* xh=0;
for i=1:n
    if i~=isb&B2(i,6)==2
        h=h+1;
8 y( o; c9 P2 l/ U# }. A0 H        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);
( f  d  t4 |' r  p- W+ m1 k0 t    end
; K% j# [. Z4 |* x+ [# {end
t=0;
/ a( S. M* _! c! H, i/ O/ b2 Sfor i=1:n
    if i~=isb&B2(i,6)==3
        h=h+1;
/ J- L3 u/ d. `  \0 r4 M        t=t+1;
        DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
, E1 z6 k7 g% S: @+ ?        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);
' ^/ x9 }  s/ |  V5 x2 e4 \0 ?# qh=0;
for i=1:n
& ^  d& c. U% A2 Z% V4 R    if i~=isb
& `( R) ]  r$ D  v        h=h+1;
4 \# K) d0 X# `1 }7 e  L        I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));
    end
8 u0 q9 O- Y( f- G0 F: {' k" Pend
I
, l; t3 ]' r3 I* {9 o: A9 L%创建Jacbi
. ^. q; L! y2 H# MJacbi=zeros(2*n-2);
h=0;
k=0;
9 S  D, {+ J8 p9 O8 F( kfor i=1:n
    if B2(i,6)==2
        h=h+1;
) a' P. W% O6 W+ q4 ?5 C& a        for j=1:n
4 Y/ s) I4 N- b) t+ e6 _            if j~=isb
                k=k+1;
7 z# n8 L6 v4 a5 }                if i==j
7 ?$ W- `5 ~5 j4 B                    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));
% a  e) w+ q0 y+ z. e4 Q3 t                    Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));
+ R4 ]  d& D; y                    Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));
                else
' e8 t9 J8 F2 T4 {+ o) p2 h2 n                    Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));
& N# W( t! g$ Y# \0 i! o                    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);
4 e7 z7 a9 n) g5 P0 S                    Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);
7 M- G$ Y; S/ k                end
+ Y( \0 q  t% T% Z4 {                if k==(n-1)
) f( E; q! i. w( h                    k=0;
. s. e: s3 A7 M. j6 `' O1 V: y/ O                end
- x/ s5 q1 ~2 N/ B. ^4 I' h* T            end
. l2 Z! o; T- d, M% Y        end
; h% q5 _% ~$ ~# [4 ]9 G    end
2 D" W5 H: Z3 k) e: \end
  [. E$ k% L) ]  u, H: ]  Lk=0;
for i=1:n
    if B2(i,6)==3
+ K6 l* d5 }6 Q0 }' O        h=h+1;
        for j=1:n
; J0 f( o7 f, P( c5 b! o( e            if j~=isb
                k=k+1;
8 C2 X3 {- R2 M# E" l% }                if i==j
: z, o1 ~5 T. v                    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));
& I( d. D/ V! Q% I( f                    Jacbi(2*h,2*k-1)=2*imag(B2(i,3));
                    Jacbi(2*h,2*k)=2*real(B2(i,3));
                else
& E& b/ u( p8 ?2 _1 _                    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));
+ w% q' T5 i% A8 _                    Jacbi(2*h,2*k-1)=0;
                    Jacbi(2*h,2*k)=0;
                end
# t8 K2 ~6 l, c4 {$ d: d2 m- l                if k==(n-1)
                    k=0;
4 P6 n2 {) b3 \  Y3 Y                end
            end
+ ~0 V2 v/ ?1 t# y4 |1 n* t" I        end
/ n) ^( E* ]- i# ~& u4 V6 t    end
* ^+ b: @& J9 }- T- g: H1 h7 Pend
Jacbi
DetaU=zeros(2*n-2,1);
DetaU=inv(Jacbi)*DetaS;
% ]7 r( R+ Q7 k9 n" e) s9 WDetaU
%修正节点电压
" T) t8 A1 T) Z" e4 Y3 {j=0;
for i=1:n
% ^1 {, O9 H( ^2 q0 s    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);
% d( m7 p, b1 n7 _( Q    end
/ T0 q+ O& }$ F$ tend
for i=1:n
3 p1 Y* ^# P2 L7 y/ Y/ ~' h    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);
    end
end
! s/ t( c& a5 x0 V8 {& XB2
Times=Times+1;      %迭代次数加1
8 K. [0 u. Y5 `" e; tend
Times
+ ?) m2 M( }& a& g; @/ w5 d
一个原始数据的例子
节点数 5
支路数 5
9 D3 ~! w+ l6 C平衡节点编号 5
精度pr 0.000001
B1(支路参数矩阵)
% S6 p" F1 ^* `. S[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]
5 k. d) N& p' R& [8 r, ?; W  \  XB2(节点参数矩阵)
# N$ a( s) [4 t' ^. v" K[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]
1 G8 T5 @% z. N" @* PX(节点号和对地参数)
7 k" d2 _5 G% o0 j[1 0;2 0;3 0;4 0;5 0]

沧山牧水

程序2：
%本程序采用牛顿-拉夫讯法对一个五节点电力系统进行潮流计算
% n=input('请输入系统节点数目:n=');
% nl=input('请输入系统支路数目:nl=');
% ph=input('请输入平衡节点号:ph=');
5 |) u, V& ]/ R# z9 J% jd=input('请输入误差精度：jd=');
% B=input('请输入由支路参数形成的矩阵:B=');
% A=input('请输入各节点参数形成的矩阵:A=');
clear all
/ }# ?$ U, p$ O+ Vn=5;
: S4 j, [3 d0 h7 G" N% Y3 P& Y( vnl=5;
ph=5;
4 K2 l: ^6 z* F- r2 Ojd=1E-6;
%支路参数矩阵
Br=[1 2 0.04+0.25*j 0.25j 1;
/ _4 s- e; I  h0 O   1 3 0.1+0.35*j  0     1;
$ x4 k; K+ ^  a. |. y+ B   2 3 0.08+0.30*j 0.25*j 1;
9 U% o7 @( U8 m2 {% a   2 4  0.015*j    0     1.05;
   3 5  0.03*j     0     1.05];
. a6 S8 J$ g4 M) I5 Y% g$ f%节点参数矩阵:2-PQ节点，3-PV节点，4-平衡节点
( {' l! H$ u+ `2 `A=[-1.6-0.8*j     1     0    2;
    -2-1*j        1     0    2;
    -3.7-1.3*j    1     0    2;
    5            1.05  1.05  3;
    0            1.05  0     4];

" E0 n" u- j5 W; P, B+ |! t%雅克比矩阵形成
Y=zeros(n);
4 e# F+ D; E- }, a- |+ \* ?e=zeros(1,n);
f=zeros(1,n);
  u. k. o' J9 g. k- n7 Y/ [V=zeros(1,n);
  }: i+ a8 p  W9 n7 s4 Ffor i=1:nl
    p=Br(i,1);
  D% z  K( C* @4 i7 }) o) u5 f: p    q=Br(i,2);
    Y(p,q)=Y(p,q)-1./(Br(i,3)*Br(i,5));
    Y(q,p)=Y(p,q);
, @9 R' A2 I% T. L$ [: g3 g) U    Y(q,q)=Y(q,q)+1./(Br(i,3))+Br(i,4);
5 E" b) v+ W; v" U6 c# V# O    Y(p,p)=Y(p,p)+1./(Br(i,3)*Br(i,5)^2)+Br(i,4);
  `$ q3 W+ k: ?; s6 x' c7 i# }end         
+ F9 M7 x; c' ^# g# }3 K" Y2 Z) Zdisp('节点导纳矩阵为：Y=');
; F8 _, d/ t$ U/ O2 C$ Hdisp(Y);
G=real(Y);
B=imag(Y);
) L/ p6 ?" u- }for i=1:n
. T, m" P9 F. z9 l" c- w    e(i)=real(A(i,2));
/ s1 @2 P* M' o+ G1 [% z  T; W    f(i)=imag(A(i,2));
    V(i)=A(i,3);
end
for i=1:n
$ |% e6 c- f& E5 R2 @    S(i)=A(i);
end
P=real(S);
Q=imag(S);
%雅克比矩阵求取
. }* X( `, R1 M9 P5 v7 |Ci=0;
a=1;
* W5 G, u/ G$ @2 q, L' l/ }" q9 v5 hNO=2*n;
5 @7 Y+ ^; f% J, g' R) U5 |( ~N=NO-1;
1 E% Q$ |. o, Rwhile a~=0
3 R( n1 M: ?- X    a=0;
    for i=1:n
        if i~=ph
            C(i)=0;
" c; M; d2 j& i. {5 b            D(i)=0;
            for p=1:n
0 V, [5 [& P2 Y& w                C(i)=C(i)+G(i,p)*e(p)-B(i,p)*f(p);
6 q1 y" r% _) p, N* {. Y                D(i)=D(i)+G(i,p)*f(p)+B(i,p)*e(p);
            end
            P1=e(i)*C(i)+f(i)*D(i);
& z4 b" S1 A# ~! q, \            Q1=f(i)*C(i)-e(i)*D(i);
            V2=e(i)^2+f(i)^2;
" ~% P: M1 @; e; s$ ]& P            if A(i,4)~=3
/ w) i& O- G6 j( u0 y! r& V  u                DP=P(i)-P1;
                DQ=Q(i)-Q1;
/ h; D; e- \: k0 w& \& ~( [( ^                for k=1:n
%非平衡节点时，非对角线元素
                    if k~=ph & k~=i                       
9 d7 ?" l8 E1 F! ^9 X7 O                        X1=-G(i,k)*f(i)+B(i,k)*e(i);
                        X2=-G(i,k)*e(i)-B(i,k)*f(i);
                        X3=-X2;
; u6 M9 i/ ]& A  {( l2 k. A/ A                        X4=X1;
7 O& H  F6 c3 s, l. t" o1 ~2 p                        p=2*i-1;
                        q=2*k-1;
8 K; C! {& k  n. O8 E                        J(p,q)=X1;
                        J(p,N)=DP;
                        m=p+1;
" j7 V: n  |/ q7 G8 w8 n                        J(m,q)=X3;
                        J(m,N)=DQ;
                        q=q+1;
: N$ [6 f' o# W0 v                        J(p,q)=X2;
                        J(m,q)=X4;                       
& N+ X( N8 M5 F+ G- V$ E/ j% h%非平衡节点时对角线元素
                    elseif k~=ph & k==i
, p/ u; n: z; g$ b5 t0 a                            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);                           
" H( c6 E, ]/ A6 ^; S# q                            X3=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);
) z) O9 Y! Q) {7 @8 I6 G                            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;
                            m=p+1;
7 P3 f6 F1 ~- P1 Q+ b! j                            J(m,q)=X3;
; q6 E# T4 X3 Y- x2 Z7 C! o                            J(m,N)=DQ;
                            q=q+1;
( R5 x9 f& v; T$ D5 t/ N                            J(p,q)=X2;
* I$ L+ @9 @* y& K/ L                            J(m,q)=X4;
                        end
                    end
                    else
! p# O% t) i6 [                        DP=P(i)-P1;
                        DV=V(i)^2-V2;
                        for s=1:n
                            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;
" X: B9 |; j/ X* G, G! q                                X6=0;
                                p=2*i-1;
' t2 K- X+ k0 X9 ]/ Y$ I4 m                                q=2*s-1;
                                J(p,q)=X1;
8 }" ?2 C. L; Q0 z; r: `: L5 `                                J(p,N)=DP;
; t$ _! r0 \! [6 c" L  Z# q" w                                m=p+1;
                                J(m,q)=X5;
                                J(m,N)=DV;
                                q=q+1;
0 H3 b% I( \5 O: ~( u$ h                                J(p,q)=X2;
                                J(m,q)=X6;
* M3 _& s2 }4 J8 A  _                            elseif s~=ph & s==i
                                    X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);
6 P4 |% m" j& C; }                                    X2=-C(i)-G(i,s)*e(i)-B(i,s)*f(i);
6 l; P% h2 H- `" d                                    X5=-2*f(i);
                                    X6=-2*e(i);                                 
                                    p=2*i-1;
9 g$ x" _0 a; N2 \- I                                    q=2*s-1;
                                    J(p,q)=X1;
. {2 A6 ^) P/ F/ c* ~                                    J(p,N)=DP;
                                    m=p+1;
                                    J(m,q)=X5;
1 ~3 `4 N- U8 y' ]; |; I                                    J(m,N)=DV;
                                    q=q+1;
$ Q2 i# _9 {8 R' l+ J' N                                    J(p,q)=X2;
                                    J(m,q)=X6;
                                end
                            end
$ g  w  F; z- c# ^                        end
0 D  F0 F1 U2 E, T# @                end
) B0 c' Y, o3 L7 z+ [/ U- T            end
1 ~/ E( i; b% J4 Q# i     disp('系统雅克比矩阵为：J=');
; e) K& ]( J* g7 s     disp(J);
%对于高压电力系统，最大元素出现在对角线位置，简单计算如下
- w* F+ Z7 J) }; R* \     DJ=J(:,[1:NO-2]);
     DP=J([1:NO-2],N);
     DU=DJ\DP;
     disp('************************************');
     disp('     ')
" y1 W$ i; ~7 H* \     disp('第几次的修正值DU？');
6 ~8 b, z  S' e9 |* ?0 D6 K     disp(Ci);
2 X" T7 `8 M7 K/ V  H- T/ J6 t     disp('***********************************');
  G4 q- M7 A8 ^4 v3 s9 S# l     disp(DU);
. s8 U+ _- [$ z     for i=1:NO-2
4 n$ t4 T! o, f; E6 C2 o         eps=abs(DP(i));
; g6 a  N% \# ], I6 B4 ?         if eps>=jd;
4 Q1 \" \4 ^' p7 f8 A% @8 G             a=a+1;
+ h1 j$ w2 v: Y, ^         end
     end
7 g+ x% T8 g. A6 l2 U7 n     Ci=Ci+1;
( j$ Q# h9 r1 q/ v3 E     for i=1:n-1
         e(i)=e(i)-DU(2*i);
1 `+ j! E* Y, R         f(i)=f(i)-DU(2*i-1);
! }$ q7 {- G$ Q1 Y! E0 t     end
# p3 X3 p# |9 X) Y% X1 U2 n     E=e+j*f;
, V$ l+ t7 f# ^  M" n     disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
6 h5 N- `; n0 u  z0 N+ z4 C     disp('    ');
     disp('第几次迭代后的近似计算值');
8 g; s* C/ Y, Y1 |  d, ]     disp(Ci);
     disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
     for i1=1:n
         V(i1)=sqrt(e(i1)^2+f(i1)^2);
         ANGLE(i1)=atan(f(i1)./e(i1))*(180./pi);
, ]7 c2 p4 k& |5 A$ k( u$ [     end
%下列输出结果均按节点号从小到大排列
     disp('各个节点的电压实际值为：E=');
     disp(E);
% m) j9 l  m( M$ ]     disp('各个节点电压幅值为：V=');
     disp(V);
; n# w: N# H. |8 [% f     disp('各个节点的电压相交为：ANGLE=');
     disp(ANGLE);
     disp('各个节点的功率误差变化');
     disp(J([1:8],9)');
     MDU(Ci)=max(abs(J([1:8],9)'));
     for p=1:n
/ I( C, Z( x3 _0 u* e3 T% x2 n! S4 S         Sm(p)=0;
         for q=1:n
             Sm(p)=Sm(p)+conj(Y(p,q))*conj(E(q));
) Z5 f' N6 p) n2 x2 R1 U0 n4 b1 c         end
* u+ m( h' C" R( e" [         S(i)=E(p)*Sm(p);
; F" h7 G$ @& k6 |1 O, p" h     end
6 ?0 p. N1 P5 u  R2 pend
disp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');
7 b$ c, \, i: d+ Z4 D& E1 ]! `disp(' ');
disp('计算的最终功率近似值如下各项值');
5 U, f2 x+ ^4 f+ K" o% L/ }disp(' ');
% X* y8 O+ l' r) F6 _9 ?disp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');
$ P" i! L. @8 u" Kdisp('各个节点功率为:S=');
disp(S);
disp('各条支路的首端功率Si为:');
for i=1:nl
    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);      
end
( G0 @" ~# \( q1 s9 T0 Pdisp(SBS);
disp('各条支路的末端功率Sj为:');
for i=1:nl
    p=Br(i,1);q=Br(i,2);
    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);
disp('各条支路的功率损耗为:');
for i=1:nl
$ ^+ h7 j: |$ X4 Y  d  v   p=Br(i,1);q=Br(i,2);
' ]' Z9 K9 C4 j/ a9 q# x( }   SBL(i)=SBS(i)+SBE(i);
end
. x) c' H4 F& |: B4 Udisp(SBL);
%绘制功率误差曲线，先用最小二乘法进行曲线拟合
for i=1:Ci
% U7 u7 i6 r: T6 V2 e6 L    CSH(i)=i;
0 c$ ~) S! d7 S% o" i7 [8 [2 g1 w1 A5 dend
P=polyfit(CSH,MDU,3);
plot(CSH,MDU,'-O');

tj612

谢谢楼主分享！！！！！

201314

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

^_^

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

windlot

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

liusen19880123

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

jimxi_6

程序可以运行？
4 H! ~- o  Y- E& b1 K??? Undefined function or variable 'm'.

sailor0305

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

沧山牧水

上面的程序可以运行