|
发表于 2010-7-28 12:24:23
|
显示全部楼层
%本程序采用牛顿-拉夫讯法对一个五节点电力系统进行潮流计算
5 I! }. M7 Y: z5 G) v( P% n=input('请输入系统节点数目:n=');3 f: a7 O( k6 K, b2 \
% nl=input('请输入系统支路数目:nl=');# J) [0 n0 \1 j# d( }& e
% ph=input('请输入平衡节点号:ph=');8 S! T6 v# `5 r4 C: H% ?, v
% jd=input('请输入误差精度:jd=');
" J. G4 E8 z- A$ B& N$ _3 \) d2 R& \% B=input('请输入由支路参数形成的矩阵:B=');
5 r3 l' q; T3 b/ `1 o* N( M% A=input('请输入各节点参数形成的矩阵:A=');7 @1 M- a$ X) d) K; `7 }7 \/ u
clear all: I# Z* h7 }! M
n=5;
" G( E; X! S; p7 f5 u" Q. B- t$ {nl=5;
: d4 W: V& b" m+ F9 N" H: Bph=5;' F: X- @% R& F( g/ m
jd=1E-6;
' k# Q- v- b" o( `%支路参数矩阵1 i, U! U+ ^. x/ P4 b
Br=[1 2 0.04+0.25*j 0.25j 1;
) `3 w' v) Y9 [+ S 1 3 0.1+0.35*j 0 1;
7 Q5 L- k0 ^) {& j+ w" d5 P* y N' z 2 3 0.08+0.30*j 0.25*j 1;" a3 y1 {: s9 i. r
2 4 0.015*j 0 1.05;
. Y/ S5 b" V# g, E7 k* m. l3 M" R 3 5 0.03*j 0 1.05]; O4 z9 f. d, t
%节点参数矩阵:2-PQ节点,3-PV节点,4-平衡节点
% T7 }1 G [6 g$ J% D3 c: \A=[-1.6-0.8*j 1 0 2;
; A* t G8 b/ D2 u/ J& V -2-1*j 1 0 2;
# m" q' T" M8 o: {% I -3.7-1.3*j 1 0 2;+ \( g/ x2 n' ^$ r2 U* ^
5 1.05 1.05 3;5 o1 l3 M+ o7 y9 [1 M7 x
0 1.05 0 4];$ z0 g* m) p2 p
# f r* f, i7 H; a%雅克比矩阵形成" {+ W: a4 { o% H$ G u8 P9 M
Y=zeros(n);
& G9 q7 E/ ~) q8 E9 re=zeros(1,n);
8 ]1 i& W/ B$ E7 O; }" s! kf=zeros(1,n);
9 b4 k" I* \2 u3 Q% E8 ~V=zeros(1,n);
' E3 M* L0 s0 v. lfor i=1:nl
& l( m: Q' N* E7 W) Q1 K9 v" ~$ e8 v p=Br(i,1);
& `) A- ~" W$ Z& P1 I. F q=Br(i,2);0 P" [7 ~5 {/ e( j6 ~% b e
Y(p,q)=Y(p,q)-1./(Br(i,3)*Br(i,5));5 S0 `/ [5 }6 W+ q* A7 P, b0 m
Y(q,p)=Y(p,q);
$ S( E2 t* B. k- a* A, u1 K Y(q,q)=Y(q,q)+1./(Br(i,3))+Br(i,4);$ g2 `, v1 ]4 D, k0 v4 k2 Q" F
Y(p,p)=Y(p,p)+1./(Br(i,3)*Br(i,5)^2)+Br(i,4);
# y: U0 x( j0 ?) S+ q% Mend 1 R# S7 k8 d4 F( v, Y* J
disp('节点导纳矩阵为:Y=');
+ y+ f8 r: a: T% `+ ]! Y* _7 |disp(Y);
/ M" v- O% I1 w% m% K h- iG=real(Y);8 e$ D) d% s+ i- v8 F+ ]5 s
B=imag(Y);: k1 M5 p& ~1 V
for i=1:n
. q. ~! k" S! B6 T9 J0 a( ` e(i)=real(A(i,2));- Y1 v, i) Q1 g2 x+ Q% g
f(i)=imag(A(i,2));
& D, o3 A' b! g V(i)=A(i,3);1 W$ O3 X" e$ g8 K
end( f$ ?9 u8 i2 B' {) z
for i=1:n; t, x, o) z' \. f
S(i)=A(i);* o4 ], V q0 l" R6 ^
end
6 G) u' L* V3 b& v+ y$ j, cP=real(S);
" q( o6 D! }# nQ=imag(S);
6 \ U9 J5 P+ B: M7 W' ^% `( L%雅克比矩阵求取
1 c+ a3 O/ d! c: D7 s6 YCi=0;
* l: f: A5 Y9 H7 F! A, L( ga=1;
: K5 Q) A; p( o$ X/ @2 NNO=2*n;+ N1 e/ L7 Q6 {' @
N=NO-1;8 p4 z0 s; t. p1 y( f1 f9 ?
while a~=0' b7 b( E2 I) b" e
a=0;& ?$ N2 k% ^5 L6 `, U( l
for i=1:n N1 W. z. O0 X. E6 r# u* V, ?. \
if i~=ph, a( q' X! D& d; X& x9 j, E
C(i)=0;
4 o$ R( S5 c( Q9 U D(i)=0;: T. j3 O6 X; E' N
for p=1:n
: P' v8 x* Y" Q$ I9 }+ b( l C(i)=C(i)+G(i,p)*e(p)-B(i,p)*f(p);
$ y; ]& k+ K, S D(i)=D(i)+G(i,p)*f(p)+B(i,p)*e(p);
- h4 i, F, ^, b3 x$ I7 w: o end5 H& \5 M( s0 G0 b. [+ D
P1=e(i)*C(i)+f(i)*D(i);
& }3 ]; m. k" r( | Q1=f(i)*C(i)-e(i)*D(i);
1 t, L& T) }0 [9 n D: a V2=e(i)^2+f(i)^2;
9 [) A( M3 i* j& ~ if A(i,4)~=3
# j5 R) W) n% X" K U4 r- I DP=P(i)-P1;& r9 X! V7 h0 L( S' b0 _
DQ=Q(i)-Q1;
$ x; x# F* s* O) K- L2 D* ? for k=1:n
3 E# Q+ c9 v' g4 R3 X9 E' ?" `%非平衡节点时,非对角线元素
) C# o0 ]; b6 S' h2 z5 B) q, [ if k~=ph & k~=i ( K2 \! W' e7 g6 b
X1=-G(i,k)*f(i)+B(i,k)*e(i); y4 f. `7 h( c# w+ o
X2=-G(i,k)*e(i)-B(i,k)*f(i);
2 P; _( b' k; j: I0 ?( B5 E X3=-X2;
1 q5 M8 j, J& b/ `* X" J0 `7 A X4=X1;
# J+ q& n0 X5 g8 l. L! {- w p=2*i-1;5 ], x c5 O. n! t" v! s _0 V' K
q=2*k-1;
, u# W+ d! p# N# U8 P) w J(p,q)=X1;# w3 a1 e; `2 E/ t8 p; [7 w* x
J(p,N)=DP;8 G* u, K5 l7 E3 w8 {
m=p+1;
0 H" y1 t& f# a1 ^ J(m,q)=X3;
3 `1 M# c" Y, O# d8 R6 W, n J(m,N)=DQ;* ]& Z! X4 n7 i% B1 p1 f! U+ D( Q4 V
q=q+1;
7 u8 z# a k9 ^( L: }. s* M6 T J(p,q)=X2;# Q: }9 O% j7 d
J(m,q)=X4; $ r* i+ \5 s* ^. d0 X
%非平衡节点时对角线元素2 f' M* a7 ~( @' c8 Q
elseif k~=ph & k==i
( E* h4 v8 \- Y$ C0 r4 O4 m X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);4 S+ }: a1 |6 K( Q
X2=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);
$ n- O1 p6 _/ {5 e X3=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);+ \# t- B. Q# W9 }
X4=D(i)+B(i,i)*e(i)-G(i,i)*f(i);: [3 U1 Q' u3 L& Q7 D# r& P! C8 J5 c
p=2*i-1;* v, T+ k. Q) o+ T3 Z7 J% A0 \( W1 ~4 a
q=2*k-1;" @# U% j, d' t n5 B$ v
J(p,q)=X1;
0 H" {7 F3 c$ C. L/ d J(p,N)=DP;2 |" u- A$ o- Q- D
m=p+1;
6 h4 n. O. m i3 Y J(m,q)=X3;9 P, l/ Y/ [: o3 }5 }- o! a& f* G
J(m,N)=DQ;
3 q" \& ?! @2 }9 e9 ? q=q+1;
8 E- A; W+ \4 A' z+ w- U J(p,q)=X2;, u: y0 t0 ^' i2 l- U$ s4 q
J(m,q)=X4;, i' t, I- Z1 x7 @# A; M w
end
& n9 p) A$ \6 @- N x9 Z4 V" | end% X6 h+ Y# T& [2 a% s1 C# r
else" p$ y z) c) h p. c, }/ r" D
DP=P(i)-P1;
% F- ]2 N& O B9 a$ {1 U DV=V(i)^2-V2;
- I0 s5 K/ M3 a1 K- a+ A for s=1:n
6 b. I& E8 ]: k# p" J4 { if s~=ph & s~=i# i( H$ h, V) E
X1=-G(i,s)*f(i)+B(i,s)*e(i);
- h) F4 T0 f X e) @+ {- o, H X2=-G(i,s)*e(i)-B(i,s)*f(i);
* o/ T2 U% j9 A$ H$ M X5=0;* {- G! }# S, q; g7 [+ i
X6=0;) u2 X7 `0 n! O8 X% I
p=2*i-1;6 u/ }6 Z* O4 P8 p
q=2*s-1;
" d& G' N# J. t+ y* S" i: b% W J(p,q)=X1;9 D5 B3 P3 S U$ U6 A
J(p,N)=DP;/ G) Q3 o1 ~/ r6 i9 w- M" y8 G4 l
m=p+1;" q0 J# d8 Z' J+ }3 y9 A
J(m,q)=X5;4 k2 O4 i) t- x$ E; s+ L2 q+ o2 Z8 y
J(m,N)=DV;$ j7 I) X$ u- k! c. @/ r
q=q+1;
$ V* @4 K3 z" Y9 c J(p,q)=X2;
+ p3 F8 t6 L( r2 p J(m,q)=X6;
$ F3 L! x/ s2 d elseif s~=ph & s==i
' f+ |% C* ^9 m; Z7 Z* O. V X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); # ]4 Z2 ~( u6 X& J3 O
X2=-C(i)-G(i,s)*e(i)-B(i,s)*f(i);
. G7 @! v u9 ?& m2 ~' b. z X5=-2*f(i);, m9 g2 E, X2 \2 w; S
X6=-2*e(i); 9 n4 K1 ]. C' n9 l- ]7 g5 @
p=2*i-1;
# o3 o& _0 U6 M! x. u$ M- q q=2*s-1;" [, y& j7 f. K8 n. N6 `
J(p,q)=X1;% h( l) Q( t. V" r. R: O( }
J(p,N)=DP;' d1 F: q2 P* F/ _
m=p+1;
5 V. f* A8 _% ^0 y5 H/ C3 g n J(m,q)=X5;
+ G0 H% b# \4 _$ `) w. t2 m5 s0 u J(m,N)=DV;
- H, T- Y9 f/ o4 j& H# [6 B q=q+1;2 A& E' A8 l/ C' H; e- w, p! B
J(p,q)=X2;
$ |1 r6 Y k1 j J(m,q)=X6;
( s$ E3 [( q; \) z/ C9 ~ end
1 p8 [4 S; c7 o+ J: L5 `+ ^ end+ {4 A$ {! `( x7 `+ u4 v3 ^! Y
end* w# Z1 X- r' B+ z/ M! C5 Q
end
& H, l6 g+ { ^9 s+ }* A) ]" J' X; ^6 N end/ K# R; v U9 y" C5 n. `
disp('系统雅克比矩阵为:J=');
& g5 J* f& R- G disp(J);
: S5 k/ O1 y+ F2 S6 z%对于高压电力系统,最大元素出现在对角线位置,简单计算如下9 a8 c- ?- U) O4 P- k
DJ=J(:,[1:NO-2]);# E: y7 ]. S/ [ N5 r
DP=J([1:NO-2],N);* J1 N4 E, e9 h, y
DU=DJ\DP;
# |: q/ j1 }% ?" K+ g disp('************************************');# p$ x! i4 U7 I4 l4 @" h& B+ c
disp(' ')
# A7 f* P* m2 r$ {' g n disp('第几次的修正值DU?');
+ c8 U& X# F5 v# X5 C disp(Ci);$ q7 V0 V6 I1 F
disp('***********************************');; s8 m# P+ W `; t! q
disp(DU);
8 {% K5 ?. t( i U for i=1:NO-2# x$ r: q* q2 z$ c, K4 I% c
eps=abs(DP(i));
; k: [- g& o8 A% X/ w if eps>=jd;
" v" L; x2 y3 p/ A: Z a=a+1; Y- ~5 `1 ?6 O- K L/ d' x. a5 k
end# T0 i( S/ z$ k6 m/ I
end0 b: o3 q3 ~5 P
Ci=Ci+1;# n$ I- A7 [/ r }
for i=1:n-1& [1 P4 S: \! w6 _+ }
e(i)=e(i)-DU(2*i);
6 H- g8 Q& S$ l5 t' L f(i)=f(i)-DU(2*i-1);
1 b. _0 N6 @% s) i- u. M! t. W end4 a- {4 a0 m: {" Q# R2 a" K
E=e+j*f;
0 A. \7 N( B/ x* ?) Z; v disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');/ v' f# P% p1 o$ U. b1 j" f3 v
disp(' ');4 U, H% r0 t" J0 G
disp('第几次迭代后的近似计算值');
% l, q2 S4 Y+ ]" [9 K$ B: [7 Z disp(Ci);
, m8 S1 v6 q1 h& n5 Y3 E; o disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
/ S# c; \9 G5 z5 {# T for i1=1:n U/ ^" D* J) \& y4 b
V(i1)=sqrt(e(i1)^2+f(i1)^2);0 m v: m* R0 _6 M' O5 h
ANGLE(i1)=atan(f(i1)./e(i1))*(180./pi);( N9 q0 q9 r+ o1 t
end
% ^ }) V) j1 p, a% U! l; h%下列输出结果均按节点号从小到大排列
9 y( }" i% }! u5 @4 { disp('各个节点的电压实际值为:E=');" R2 q9 A1 C i; Y0 |" x: k
disp(E);
5 j! ~: ^" l. d) g disp('各个节点电压幅值为:V=');0 u, g' J! j/ ~
disp(V);
4 J, m3 F" d8 a disp('各个节点的电压相交为:ANGLE=');
4 @- a' a0 e+ {, i1 @3 Z2 k disp(ANGLE);
3 R( \* N+ C: ?: W' p disp('各个节点的功率误差变化');6 N1 z; ~* V% q' i; h8 ?/ X
disp(J([1:8],9)');
) {4 N) ^$ f3 t' W4 C; D2 B MDU(Ci)=max(abs(J([1:8],9)'));9 r5 q- n" c, Q3 k* P
for p=1:n& x/ d$ W, G2 O% v+ {% Q6 W+ p* a" s
Sm(p)=0;( K, w6 C/ V5 ?, S( d) ?/ C! |
for q=1:n
# F6 B( O$ E: |4 u8 u! U, r Sm(p)=Sm(p)+conj(Y(p,q))*conj(E(q));
8 r' ]1 p1 n2 i7 T- f end$ o4 I; t0 M- h% L) s, q) u6 [2 p
S(i)=E(p)*Sm(p);& R" f; O9 x# D7 C% ~
end9 B7 ~ ^$ b0 B3 u
end
2 F5 C& _8 w/ e: d% C1 f6 o- c+ ]9 G" Hdisp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');
9 G3 `5 ]) k$ _disp(' ');
0 q4 `7 I) i9 c# `% Cdisp('计算的最终功率近似值如下各项值');
( l! X( f3 D2 x7 D kdisp(' ');# G/ d% {7 H+ Q. l$ c
disp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');
. d0 g) e: z% o0 @3 ?% hdisp('各个节点功率为:S=');* T$ P# x! a% s- O: ?% V
disp(S);- s3 i/ x" E l2 p
disp('各条支路的首端功率Si为:');
7 T2 D! W1 V/ \8 G4 w7 q, p! N5 _$ @for i=1:nl # @: v8 G. L; `
p=Br(i,1);q=Br(i,2);6 p3 Z" e. [; v* w$ B C- Y
A(i)=E(p)*[conj(E(p))-conj(E(q))]*conj(1./(Br(i,3)));( l0 O- N! E; j
SBS(i)=[abs(E(p))]^2*conj(Br(i,4))+A(i); * c* h8 p5 Y7 ~9 O# o
end
" E8 J8 y5 w6 Z- ndisp(SBS);
" O f+ H5 d1 [$ K, hdisp('各条支路的末端功率Sj为:');" h' c# D4 P9 r, X9 Q4 X
for i=1:nl
/ |5 E% X9 w# n) I. C7 ` p=Br(i,1);q=Br(i,2);0 P W7 M; h2 y
A2=E(q)*[conj(E(q))-conj(E(p))]*conj(1./(Br(i,3)));
, |9 e- L: P# u0 @ l8 r SBE(i)=[abs(E(q))]^2*conj(Br(i,4))+A2;
4 c- q- | [+ M" N5 Y! B. Gend7 Y: ^: b, b6 n, T- P
disp(SBE);
7 }- e ^. q# Rdisp('各条支路的功率损耗为:');
( e9 [$ J7 W) z' J* N& ?" l5 v' H; xfor i=1:nl
C) X7 d& m4 \% ^) q7 J, H; ^ p=Br(i,1);q=Br(i,2);0 A, H, X% A3 x; c5 g) o
SBL(i)=SBS(i)+SBE(i);" ~) _/ [) K1 I+ W$ Z* \; ^' Q
end
) }3 ^; }" R m! m2 M2 ]disp(SBL);
& k) _- Y( Q6 D, N* N. I%绘制功率误差曲线,先用最小二乘法进行曲线拟合
2 d( h0 P7 l' P+ J3 a* T5 qfor i=1:Ci$ s# Q i; \& |
CSH(i)=i;& C; E: @6 O$ t( C8 M
end" _ D. B( x4 X0 y+ Z
P=polyfit(CSH,MDU,3);
$ H4 ?$ Y9 Q2 I* dplot(CSH,MDU,'-O'); |
|