|
|
发表于 2010-7-28 12:24:23
|
显示全部楼层
%本程序采用牛顿-拉夫讯法对一个五节点电力系统进行潮流计算
/ w# m6 _2 I" e) N% j% n=input('请输入系统节点数目:n=');
, b1 N, l! F& i& p4 n% nl=input('请输入系统支路数目:nl=');
1 p0 q6 B' m: H/ a1 ^% ph=input('请输入平衡节点号:ph=');
% f% u0 E/ c& [( s4 G: G O% jd=input('请输入误差精度:jd=');8 h6 ]1 x; z9 K& n9 R# M7 B
% B=input('请输入由支路参数形成的矩阵:B=');1 q8 i2 C$ M0 Z' k
% A=input('请输入各节点参数形成的矩阵:A=');
% `( T8 y; u7 q( ^3 T4 Vclear all
{& q$ ]; a' [" i8 O9 an=5;
# P; j5 c% f0 ?$ O$ Y3 r# O7 pnl=5;
0 P3 Q6 |. w% D+ C9 |ph=5;6 L2 r& E3 e F& c
jd=1E-6;/ D+ l' d. f$ J' _! d2 m
%支路参数矩阵- A, B. _' Q ?# j4 ?/ c# t! Z# A
Br=[1 2 0.04+0.25*j 0.25j 1;4 ?+ Q+ P# E5 j1 C
1 3 0.1+0.35*j 0 1;
7 f" h# n5 H1 ?7 J7 H 2 3 0.08+0.30*j 0.25*j 1;/ ~: h, T9 a8 u' J6 s
2 4 0.015*j 0 1.05;
) k9 l2 t; O: M" Y) r 3 5 0.03*j 0 1.05];
" i9 E; f7 \, |& q%节点参数矩阵:2-PQ节点,3-PV节点,4-平衡节点% S, v+ x' B0 S3 B' h8 z
A=[-1.6-0.8*j 1 0 2;. r& W$ R" a6 r" z2 U1 m
-2-1*j 1 0 2;8 I4 A5 m0 y6 B, G
-3.7-1.3*j 1 0 2;4 a$ ?1 x5 e7 l) m' U0 _5 k8 u
5 1.05 1.05 3;
# h& {) z& A. q$ J9 u# } 0 1.05 0 4];
+ F/ F! m3 A; E$ `3 m: i
6 N* f4 n& f/ ~; b5 Z1 p%雅克比矩阵形成
0 i( w0 V* ^8 u, R3 {. @Y=zeros(n);( o! `+ `! h+ x7 B4 @" S+ g- y
e=zeros(1,n);
% G9 o4 K; d% K' b% Hf=zeros(1,n);
3 a8 \! V, w6 X& u# ^+ C5 rV=zeros(1,n);
3 f- p) \; o5 G9 n( {2 F9 ]for i=1:nl
7 ?0 r/ v1 E: ^- F" e% u p=Br(i,1);8 I9 Y! `4 j- f1 N7 z, B
q=Br(i,2);2 K, x, r4 l# Z3 L
Y(p,q)=Y(p,q)-1./(Br(i,3)*Br(i,5));3 z4 w8 c, {: X) Q7 g
Y(q,p)=Y(p,q);: W9 I) B; e7 u! i" w5 W1 {# P
Y(q,q)=Y(q,q)+1./(Br(i,3))+Br(i,4);; n+ { t2 S! B, Z, H4 n
Y(p,p)=Y(p,p)+1./(Br(i,3)*Br(i,5)^2)+Br(i,4);
1 |8 B' h, g* A4 K0 x9 N- b& d- z: l0 ^end . c+ |, L2 E) O
disp('节点导纳矩阵为:Y=');
6 U( F* n' ~) R" \# w+ ?& L0 G8 _" \disp(Y);+ P1 L& P4 y. J. }
G=real(Y);1 B/ \+ h6 X& ^! m! b1 z* R! D
B=imag(Y);7 C! N# f# ?4 A
for i=1:n$ b* c$ k! [+ M& A% `4 Y3 U- l
e(i)=real(A(i,2));8 {! }- Z$ f& O+ {
f(i)=imag(A(i,2));
* N4 {3 n! I; x" t, s V(i)=A(i,3);& W2 ]+ d( i0 b A2 s
end
' C" J8 E/ c. C; T# o4 o5 w+ w+ C' i& jfor i=1:n
) h9 k C" h X: o1 J S(i)=A(i);3 f2 L- }+ g& O9 d% m
end
/ W, c. i; _: N% h5 g! \P=real(S);
; U4 B, s! F7 @( B6 qQ=imag(S);
: c% |( z& Z2 n" ~: j. a' w%雅克比矩阵求取
/ K2 M9 i! k$ hCi=0;! b+ c0 n+ Y y
a=1;
0 U$ e1 g! h/ G% E# p8 D' fNO=2*n;: \/ f0 M Q7 t6 Z( A% P
N=NO-1;
3 ~! d7 Y' L9 l. F3 | |4 k) Hwhile a~=0
1 Q/ ^' q& S% Y5 J' \" o- g a=0;5 J4 j. F4 N& ~+ l, q) F4 ?
for i=1:n
% O0 W L. t. F5 K% D4 r& p G! K( W if i~=ph
5 S. }4 t' J6 ~- T C(i)=0;
]/ Y5 }* P6 z" y3 Y D(i)=0;
5 n* ^* W! \: n, M! F8 [ @ for p=1:n
R+ W6 C/ s' n( `9 l9 c! |# a C(i)=C(i)+G(i,p)*e(p)-B(i,p)*f(p);
! k! H% x% p6 Y$ T) F7 \9 r# d& I D(i)=D(i)+G(i,p)*f(p)+B(i,p)*e(p);$ T9 |. H8 t/ Q6 R9 ]
end' C3 b8 H! g8 R1 n+ e# f
P1=e(i)*C(i)+f(i)*D(i);8 L/ c8 F0 V1 t6 H9 x( X4 Q
Q1=f(i)*C(i)-e(i)*D(i);
# W9 c. Q1 L: G6 {4 @% g$ J V2=e(i)^2+f(i)^2;
! r3 K' x' L i! m; n: {1 l if A(i,4)~=3
/ F) F8 H% K3 N3 q/ w DP=P(i)-P1;
# N3 F" @ R, a' F9 j DQ=Q(i)-Q1;9 D' c# q* l" ]$ K; e" r- ^
for k=1:n1 h# j; {1 y7 i1 a* W# J3 Q: w
%非平衡节点时,非对角线元素; \7 L! |7 r Q5 T% N
if k~=ph & k~=i ; `6 A, s7 l1 U; H
X1=-G(i,k)*f(i)+B(i,k)*e(i);
6 c7 C5 z3 @! m+ {$ T X2=-G(i,k)*e(i)-B(i,k)*f(i);5 E3 c. {2 W7 X8 }' Q
X3=-X2;0 a. v0 o/ l$ d
X4=X1;7 @9 Y6 v- v& c5 z9 @
p=2*i-1;
8 ]3 H/ O9 R% M- `1 Q' l q=2*k-1;/ G2 k7 C s% _ z' `# q0 w
J(p,q)=X1;: i' |1 u. K- e
J(p,N)=DP;2 M9 J! O7 }' u
m=p+1;
! Y0 H& G+ O v/ N J(m,q)=X3;9 ]" ~( ?7 t3 x `% P
J(m,N)=DQ;4 ]) n1 [# M" |0 Y
q=q+1;
1 N5 R8 [7 F4 W# f9 S# w$ S9 T+ K* H J(p,q)=X2;
( i' U' G) ~$ d+ J& i; Z1 s0 B J(m,q)=X4; ! F% Q; t3 J5 c+ g) |
%非平衡节点时对角线元素* V' @5 T6 a0 o, k' A
elseif k~=ph & k==i
# n" W2 n" V+ x% d; D) g X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);5 k7 M% \! W( P, P4 _" S
X2=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); + P- e( e" x6 ~
X3=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);& H0 `$ l9 I8 H. m/ Y
X4=D(i)+B(i,i)*e(i)-G(i,i)*f(i);
6 ?% t2 t$ f+ X X5 E( S6 y$ y8 E" D p=2*i-1;% ]( U* D- K( M3 f2 H, m
q=2*k-1;
% `% _7 I& y" `+ X8 Q J(p,q)=X1;. |. a. {* p# l4 P' E
J(p,N)=DP;
% z& b* W9 x' `" C m=p+1;
6 }) l" L; X# ~( v; [5 l) ]. l J(m,q)=X3;
; _. k7 K4 ? k" g4 l& k( }/ W& Y J(m,N)=DQ;
: S$ ]. \, ^2 X/ i6 C) n# v, g: } q=q+1;
8 _% f* {, X, M) \* g6 i# O$ ~+ K J(p,q)=X2;
- G" @& ~3 e5 g: Z. ^ J(m,q)=X4;" O) O8 \5 ?( Z) V3 V
end9 U. i+ y- C4 l- t% J4 v! R& i9 e
end
4 g. a1 J( \. i) ]( ?: T0 G! f else) J! e" u+ t0 e4 p
DP=P(i)-P1;. [. ?2 r s+ X: {2 |# h: w
DV=V(i)^2-V2;
3 q a4 h7 i" a& t1 } for s=1:n' a' @* Q/ H, i2 u4 \, d; e
if s~=ph & s~=i( T! `3 V, Y8 u5 G& j' }
X1=-G(i,s)*f(i)+B(i,s)*e(i);
8 |9 z0 {" O0 F. D* ^5 c: L X2=-G(i,s)*e(i)-B(i,s)*f(i);
! J6 k8 ? n# G3 z- d" } X5=0;
$ I1 x- a. r9 e6 B, u X6=0;
) J- W6 H2 f' w- X; X0 w3 D p=2*i-1;
! c6 Y/ Q1 E# Y8 X; [& V/ k q=2*s-1;/ p- j9 u# F* g4 q& |6 V
J(p,q)=X1;9 s {4 T7 [6 l+ `& l+ [8 e! K
J(p,N)=DP;
* [+ G- L) }/ h: x! s m=p+1;
9 w* k. ?) u S5 i/ W J(m,q)=X5;* V: z/ z5 ~2 b4 Q; k5 i
J(m,N)=DV;
2 Y+ z) R! `0 `2 O q=q+1;
* w/ S, L3 W( ?" [ c4 y) I" l3 _ J(p,q)=X2;. G3 k& B! p! h' e1 w
J(m,q)=X6;
0 e% I) y, S5 n. Z. x& s elseif s~=ph & s==i5 e0 }0 E* C3 ]1 J9 y$ n/ w
X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);
. N0 x! W" _: _* Y! r X2=-C(i)-G(i,s)*e(i)-B(i,s)*f(i);
+ ]9 f2 {/ E a+ U$ Y: G7 F X5=-2*f(i);, i3 P1 e5 Y+ o3 V9 a! I
X6=-2*e(i); * z% a y/ W7 N9 W& y3 D
p=2*i-1;9 J9 G$ h5 b+ Y8 E
q=2*s-1;
5 q0 M2 J- \, N' ?0 w( x J(p,q)=X1;, s* W' c6 a# Y I4 R" i/ e2 F
J(p,N)=DP;
8 W, o4 E: T9 J6 B- u$ d) _ m=p+1; |( m4 t9 T# U' ~4 p
J(m,q)=X5;
: X0 |1 ?: \+ A- I J(m,N)=DV;3 o+ w" A. s$ H1 h4 [
q=q+1;
* K+ @6 ]% T7 s3 k- E/ A J(p,q)=X2;
1 x$ A2 N) V8 s( p% @! {9 i J(m,q)=X6;
+ [8 G4 _6 w7 P2 F, r. { end+ U# I/ }! S1 ^
end
' S7 B( Y6 c3 \ end
: v6 j- G8 ]+ b1 s) {# N5 {4 | end# }4 r% h( w& Q% I# e7 e
end( E3 n0 ]7 f& b' L; w( H, c$ g
disp('系统雅克比矩阵为:J=');
' C4 P8 {4 h% j b$ ~ disp(J);
' t9 K+ u+ i' K8 a; g% c%对于高压电力系统,最大元素出现在对角线位置,简单计算如下6 s Y/ _2 l: y7 s N9 Z3 z( i
DJ=J(:,[1:NO-2]);
: }( {0 \- P! @2 ~& a9 B+ @+ n DP=J([1:NO-2],N);
1 o) W- c4 C( r5 t! p DU=DJ\DP;
+ x$ i8 {) v& ]) x# b disp('************************************');5 t' }2 u8 b' F6 q) O* }
disp(' '): s0 D3 \( d3 z" C6 D1 j0 O
disp('第几次的修正值DU?');
2 A& s. ~! { S5 _ disp(Ci);" ~2 M) P& F7 M
disp('***********************************');
+ G( w: {& g- U5 C6 a% `+ S disp(DU);
, n" j8 o k, Y0 U9 J for i=1:NO-2
F" o# ]5 [4 X! _7 M( O7 I8 `: _ eps=abs(DP(i));
' I0 z1 h' L( x. ? if eps>=jd;
! J3 W3 `4 J( Q a=a+1;' a' b0 I4 K4 X
end1 G/ i/ n3 N) Q X/ y
end
6 ^6 b$ T, X9 D3 a; v, T+ l Ci=Ci+1;
- R- B& B* F) v% G m' h/ t O for i=1:n-1
" j4 x( n0 g. m e(i)=e(i)-DU(2*i);
$ ], g* G) p7 N7 n6 I% m& J3 _ f(i)=f(i)-DU(2*i-1);
. B1 }. l& |' w( B' }& w% T4 n$ z end! B, [+ r7 h" d* T4 L$ b$ f
E=e+j*f;" Z2 B5 ^" G1 ~& b( Q# Q2 A- _. M6 k
disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');/ l# m7 E6 b7 Z* W( _* ?4 p M
disp(' ');0 v0 @3 b1 e% z( j* \3 k. n4 N
disp('第几次迭代后的近似计算值');/ o9 @% N- V- _( P
disp(Ci);: J9 H! f) V9 ^
disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
5 O f! G+ y/ j# [2 v for i1=1:n
, ]# S9 X5 ]: G! w! |3 D6 p! y V(i1)=sqrt(e(i1)^2+f(i1)^2);
; A# z7 p! r ~. S9 K3 Q% F ANGLE(i1)=atan(f(i1)./e(i1))*(180./pi);" A# J) G+ A7 B3 } j
end; Y' @2 c0 e1 W1 m
%下列输出结果均按节点号从小到大排列# o k ^' |5 k4 s6 I: U
disp('各个节点的电压实际值为:E=');
0 m) w" p# i9 b9 X# X disp(E);% {2 R1 |- @+ A6 n
disp('各个节点电压幅值为:V=');
) J$ h0 p- E. U! z# b2 B( m disp(V);! |- H9 [2 N# E
disp('各个节点的电压相交为:ANGLE=');
% @1 B L& d+ U* c% e disp(ANGLE);7 k0 n0 T' o3 `5 s2 i* i3 r
disp('各个节点的功率误差变化');
. z- h1 m) q, ~# m0 Z; Y% N' h( d disp(J([1:8],9)');' A4 v* |5 _$ w- q3 m. k: h" _( r1 K {
MDU(Ci)=max(abs(J([1:8],9)'));: A3 L7 @0 H* ]/ T! F+ f1 n* o" K3 |
for p=1:n c1 U7 n: v8 m* P# w8 u* y
Sm(p)=0;
. O) E- g8 l8 \8 }! ? for q=1:n+ z; V* Y, v3 ]4 j* K0 |
Sm(p)=Sm(p)+conj(Y(p,q))*conj(E(q));
# s3 O- g+ }$ S" P. n end9 w4 ]# ?& n) L. B5 h: H1 R
S(i)=E(p)*Sm(p);& l* T4 c9 j3 g7 h
end9 v+ F/ a/ e" h8 U+ g( i) t
end
, ~; d' X* m1 s# d5 Tdisp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');5 h; S5 S9 b/ ?5 E+ l# _+ |# N
disp(' ');
W0 c. u$ c0 |7 g1 r1 c; Ndisp('计算的最终功率近似值如下各项值');( w8 O" p4 ]. N
disp(' ');
. } Y+ D! h8 ~+ ]! J4 p: qdisp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');1 {- t* j/ y9 F8 f, Q
disp('各个节点功率为:S=');8 \" l5 w/ Q7 B2 ?6 f
disp(S);
z [/ s S7 A. Wdisp('各条支路的首端功率Si为:');
) a# I& b( x, ?6 ?4 L7 s9 c; H( ofor i=1:nl # n$ j. N" Q `$ I- N" N- m8 R
p=Br(i,1);q=Br(i,2);9 w( [$ J9 ~7 B0 H
A(i)=E(p)*[conj(E(p))-conj(E(q))]*conj(1./(Br(i,3)));. e" _1 h2 l% l5 G+ A
SBS(i)=[abs(E(p))]^2*conj(Br(i,4))+A(i);
3 h) B3 H7 w4 V1 @6 z% kend
) J& ?+ i2 E5 Hdisp(SBS);. h3 y# o5 X6 u, H' ?
disp('各条支路的末端功率Sj为:');( W$ `0 A5 \! x% i5 F
for i=1:nl/ x+ O d# P) l2 `) g1 G
p=Br(i,1);q=Br(i,2);- }# D @3 [/ Q2 e9 J8 [$ L) f
A2=E(q)*[conj(E(q))-conj(E(p))]*conj(1./(Br(i,3)));7 X4 P. x! G' o
SBE(i)=[abs(E(q))]^2*conj(Br(i,4))+A2; s5 E5 p% j! [$ g$ ]
end
* g* g& s" q% u Adisp(SBE);: f/ H c7 O/ I& p
disp('各条支路的功率损耗为:');+ p" g' S5 L9 N" R$ j8 N5 v$ D
for i=1:nl0 i% |5 ]* j0 v- t4 L
p=Br(i,1);q=Br(i,2);
: u& [2 u# R1 J5 M% u SBL(i)=SBS(i)+SBE(i);! O" }$ J& R. w$ S: z9 _- t
end
, Y& X' I2 w9 K/ v, P0 vdisp(SBL);
+ B s* v: M% {1 r$ p& t%绘制功率误差曲线,先用最小二乘法进行曲线拟合( H* A5 S6 E; G$ |* I# i+ ~. H4 u
for i=1:Ci! G* D7 d4 K# P# B
CSH(i)=i;
9 w0 n1 L' q) w: z( I- s$ q2 }* Pend
$ N+ m; I5 ~8 tP=polyfit(CSH,MDU,3);
; z+ T. |! d& x( A; hplot(CSH,MDU,'-O'); |
|