|
发表于 2010-7-28 12:24:23
|
显示全部楼层
%本程序采用牛顿-拉夫讯法对一个五节点电力系统进行潮流计算
! y4 K0 Y H% N* a; C! C% n=input('请输入系统节点数目:n=');
" Y1 D5 ^# U* f5 C$ J3 c. |" N% nl=input('请输入系统支路数目:nl=');, L& a# Y. S6 W& H9 F" v2 U% t
% ph=input('请输入平衡节点号:ph=');
! `( l/ U+ |' s/ |; S/ G9 |% jd=input('请输入误差精度:jd=');
2 W" q- Q2 k& h6 X! z9 E# {% B=input('请输入由支路参数形成的矩阵:B=');& R- l6 A- g, C. j
% A=input('请输入各节点参数形成的矩阵:A=');
" ?8 E! E* p( t yclear all; \: M- S) R; k) y4 X4 `$ I& Y0 @
n=5;% f/ B3 r, `! J# ~9 ?+ B; x
nl=5;
% Y8 w6 j& L) t# N2 g4 ? Fph=5;
8 c) O; m5 _9 f9 J/ v, a6 Pjd=1E-6;4 d2 }& v+ j# p/ r% `
%支路参数矩阵
* r' Z$ n) ^4 g( F8 a& {Br=[1 2 0.04+0.25*j 0.25j 1;/ Q2 _& _% W9 S- l
1 3 0.1+0.35*j 0 1;
& g' N3 f5 m$ C I# m9 a' U x6 t 2 3 0.08+0.30*j 0.25*j 1;. X# d. ?/ s- s& j% N3 x" N9 H
2 4 0.015*j 0 1.05;1 ]. m: ^7 |0 U) a: W2 Y
3 5 0.03*j 0 1.05];
; a; A9 @+ q) V%节点参数矩阵:2-PQ节点,3-PV节点,4-平衡节点
1 F' C) ?1 ^4 f4 Z6 g6 VA=[-1.6-0.8*j 1 0 2;7 u" I0 j8 |8 S( ^2 z) ~% H
-2-1*j 1 0 2;0 N9 @6 U+ J7 T/ |9 X: ^# X
-3.7-1.3*j 1 0 2;
% Q, O7 e2 }* N. D2 S7 A! g 5 1.05 1.05 3;% f) t" Z$ A) T" }. b5 w
0 1.05 0 4];
w$ Z% `- W7 S5 p2 w# M) s- V) B& g( y) L2 O) F( x
%雅克比矩阵形成4 C- ^# @" f; D
Y=zeros(n);
- M7 N4 B6 M8 P8 |# v2 `! m7 Ie=zeros(1,n);
! p0 D3 V. W; [+ ~$ r9 \f=zeros(1,n);) q; T( C0 R5 |. `- }4 Y9 W
V=zeros(1,n);
0 q+ f5 w- r4 V# N E: gfor i=1:nl
, a; B; `; D; f @ p=Br(i,1);5 @1 R, W* L' [
q=Br(i,2);
s$ S2 y% Z/ L9 U' k+ ], e/ m2 F Y(p,q)=Y(p,q)-1./(Br(i,3)*Br(i,5));1 U# O, V5 P9 B9 ?" F( C
Y(q,p)=Y(p,q);: N; W! q: {+ o. M" M n% r
Y(q,q)=Y(q,q)+1./(Br(i,3))+Br(i,4);+ t7 p& V; T. }, j6 q, c2 W
Y(p,p)=Y(p,p)+1./(Br(i,3)*Br(i,5)^2)+Br(i,4);
2 H! h6 O3 j1 W: I" L& [end
7 a! U1 A% e1 y$ W b# Pdisp('节点导纳矩阵为:Y=');
# U" V, p4 F" {/ y% rdisp(Y);) R4 D( p8 h1 E: m9 c( p1 Y3 x
G=real(Y);5 {. R' `$ p. r! P! A
B=imag(Y);/ N0 h, z/ W! i. J2 A
for i=1:n% A" A% P6 R( x$ L' Z5 e4 }
e(i)=real(A(i,2));4 [& }% m+ T2 {# ~9 N
f(i)=imag(A(i,2));/ t' j. U4 R3 ?0 D$ u3 w, i
V(i)=A(i,3);
; M; `3 o, E% j; S6 _) Xend
6 f$ G1 ]! \6 g% h6 d, ^for i=1:n
; l% M, d2 A; U, u O5 E3 \6 I S(i)=A(i);% Y3 k5 p: l2 m7 {
end6 b: [2 E' u: ?
P=real(S);
& L" V+ U5 g: cQ=imag(S);, ]) {0 e/ P5 N: x1 b
%雅克比矩阵求取. K e6 T- f/ w" `8 y# C
Ci=0;2 p e2 V1 f! y% k+ r% Z; B
a=1;
- i& X* j) Q) h8 aNO=2*n;8 ^( ^( }) {, ]1 ~8 J! z4 H
N=NO-1; S* w+ |; ]! c3 m5 q+ h, U
while a~=0
. M! [& t* o' {, r6 J a=0;+ y/ v. r$ W, K; \, m6 z, z5 O
for i=1:n
* O0 e3 d9 D, v- o3 R/ c if i~=ph
( C! x8 O1 g8 a; Z C(i)=0;
' c& w5 ]2 X# W) y* P6 ~5 I5 i$ e D(i)=0;
8 ~. [. f( l, D% ?# p/ m A for p=1:n
% }- i+ E* `; O4 Z+ @ C(i)=C(i)+G(i,p)*e(p)-B(i,p)*f(p);
* C+ R; h& j3 L- f* L7 e2 O D(i)=D(i)+G(i,p)*f(p)+B(i,p)*e(p);
8 e; P) d/ P6 |, v2 o. A! Z8 s end
0 l3 I4 U. o) T P1=e(i)*C(i)+f(i)*D(i);
" M) l- B, Q0 ~" D' R, z Q1=f(i)*C(i)-e(i)*D(i);
5 f m2 |* M: a4 ^, L. `% k V2=e(i)^2+f(i)^2;. b1 a' u1 \/ n- I: ~
if A(i,4)~=34 i6 a' M* Q( _4 }1 ~5 V2 o0 V
DP=P(i)-P1;" |7 j H0 `! Q! Y2 [! A1 ?
DQ=Q(i)-Q1;
" C3 {$ C6 |7 K4 |+ c for k=1:n" N A: Y1 I; p
%非平衡节点时,非对角线元素1 x- Q* L+ x: Y6 N" S
if k~=ph & k~=i
; [% M( d# J- [3 O/ Z$ X$ m X1=-G(i,k)*f(i)+B(i,k)*e(i);
: k* X' ^; G% ~1 i: J7 @ X2=-G(i,k)*e(i)-B(i,k)*f(i);
9 q: C+ _! s; m& V5 @* ` X3=-X2;$ m! m0 X9 z6 c7 x9 s. b3 P5 L
X4=X1;
+ L; A4 S1 Y3 l# w/ y p=2*i-1;% r" X. P/ a. g; B
q=2*k-1;
: w: R7 ]% r8 P1 E7 } J(p,q)=X1;1 S6 |3 e" L, d. j4 H/ M
J(p,N)=DP;
1 H# l A2 y5 C! p m=p+1;
! z" n4 \/ A1 }/ w J(m,q)=X3;
* C7 M: N4 _6 } J(m,N)=DQ;
6 _, t+ h! A/ ?2 S' d q=q+1;
( |4 c# v7 s: h J(p,q)=X2;4 a3 d3 p5 G, T. |3 c2 W3 c! k3 J
J(m,q)=X4;
2 f9 i2 s) U8 S% R0 V5 w( C%非平衡节点时对角线元素
! K, C+ x& f# M0 G6 c3 v" j) v: r elseif k~=ph & k==i
8 i/ v0 y( D) V6 B" N" z" m8 x X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);" y6 [' U* p! @0 N0 r
X2=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);
2 p7 F& j5 c; @" Y$ Z; C9 I X3=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);1 n) W, I6 v" F' ^6 m
X4=D(i)+B(i,i)*e(i)-G(i,i)*f(i);
0 A/ t. R9 `& d: [2 a Q p=2*i-1;0 h0 B/ g0 F0 w! A8 u1 T( r- p- S
q=2*k-1;
, `+ X& O. Y1 V2 ] J(p,q)=X1;
/ b9 s9 G( d& Q7 r9 e2 l+ V( ], z J(p,N)=DP;- j A. j- y. T. o9 p% e
m=p+1;
* _3 B0 r! B# {4 D" I1 e- U. J J(m,q)=X3;
- ^0 {: x3 c. {- J2 f5 a2 v! g h J(m,N)=DQ;
& \3 K# Z" d3 Z1 m q=q+1;6 @: b8 p0 t: k( k/ \) {
J(p,q)=X2;) l: P( g* ^8 c2 Z* g. O* |
J(m,q)=X4;* G7 g* W, ^; Y% t# E7 v4 }3 w
end% v0 Q% P+ }9 Y4 {1 J
end1 |/ F7 v/ L9 G9 e: ]; Y( u
else* ]0 S1 F' p4 [
DP=P(i)-P1;
8 E5 j& S1 [' C7 A% D( q; J DV=V(i)^2-V2;
$ T! Z( N1 ^' S* m for s=1:n8 S5 P* g% z. g" d; o; C& y& d
if s~=ph & s~=i; C) D3 Y0 Z: g6 J9 K r% r0 b
X1=-G(i,s)*f(i)+B(i,s)*e(i);
( y) ~. c! W* Y9 ~9 F C, { X2=-G(i,s)*e(i)-B(i,s)*f(i); 7 o7 S0 N' Q/ V8 U9 d! j
X5=0;' J; u0 J2 F: t0 l: |
X6=0;6 Z2 H2 P* S! ^ I8 _% H
p=2*i-1;7 t# P* ` V& W% r
q=2*s-1;8 H8 R; k: x# `
J(p,q)=X1;. r9 k/ i( y& C. q
J(p,N)=DP;
* E0 D$ A: |: G+ S& z$ w m=p+1;) ^' g' A! H% s& i
J(m,q)=X5;
3 l: Z! W0 z/ D J(m,N)=DV;
6 ?: `3 s5 R9 z+ ?+ R7 S q=q+1;
4 }% J/ q9 P( Q" r7 d J(p,q)=X2;+ ^) J# b% Z! d4 f2 D
J(m,q)=X6;
- w3 _" h# {( s2 F d9 v) H! Z elseif s~=ph & s==i
3 P. }% p- U5 h: E0 t; q X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);
% Z( w w1 b* \1 M8 ^ m X2=-C(i)-G(i,s)*e(i)-B(i,s)*f(i);
% _# D) S% G0 x: ` X5=-2*f(i);5 o* X! q# R1 R) N: X- ~; X; A
X6=-2*e(i);
0 y! q3 @( E8 E% S p=2*i-1;8 N) g! B$ T+ i. U6 q4 W
q=2*s-1;9 W. r$ U5 O. {5 q
J(p,q)=X1;9 [6 S, Z. i7 W8 V5 v0 d
J(p,N)=DP;+ @) \& x' S6 r2 ^3 L6 Z
m=p+1;
! A4 K4 q7 k9 A6 M J(m,q)=X5;( T5 q0 V; H+ M; i3 y/ T: ^% L
J(m,N)=DV;* L+ I4 ~6 x( v1 A3 K
q=q+1;2 D$ x0 ^% v7 A: Q' F0 b
J(p,q)=X2;
+ U$ W7 G* `) _; _, e J(m,q)=X6;% F4 ~, q$ @7 V" s
end$ N# `/ m% h I0 T7 t/ i
end
$ S9 R# m* K7 _! f. h2 U end# _8 c, U& r* Q" w# k4 Z
end
3 Y' v+ ]& M9 ^+ _2 W end; W5 j+ f( [6 M& ^
disp('系统雅克比矩阵为:J=');' P9 ]' a2 E7 Z8 F) U
disp(J);
0 \# O, D0 s7 {* t u%对于高压电力系统,最大元素出现在对角线位置,简单计算如下
/ z7 q0 h* ]" _, t DJ=J(:,[1:NO-2]);; ^3 T- |# ?" ^7 X8 W
DP=J([1:NO-2],N);8 R4 n+ ]) F3 I+ R5 I9 s. n
DU=DJ\DP;
, c% P- d4 F0 ?! w; ?( g' g5 R9 X disp('************************************');
& h7 Z, P. V+ d$ F, ~5 ? disp(' ')3 }' S( C& U% T o8 m N/ \
disp('第几次的修正值DU?');' S( _; i8 K2 {5 o
disp(Ci);
+ O1 z5 P( z! u6 D2 ` disp('***********************************');
5 M8 ^- D {6 S disp(DU);9 b% e# k* N- E/ O5 ?. Q- e8 ^+ }1 l! ?
for i=1:NO-2! E9 X# e0 S/ i0 {/ n& |0 P. a
eps=abs(DP(i));; P' |7 F6 R% {' e" H, W
if eps>=jd;
5 w4 A( h9 |2 t1 [1 C& q) j a=a+1;; i& K1 x0 i0 {: i$ g
end
1 e$ x: C) Y* f( K) ^3 q3 a end
\! ^% r. \3 \8 A s; x6 r* S Ci=Ci+1;1 W$ U q6 M' r9 r: m4 ?, M7 [
for i=1:n-1
( v0 Y' F. p4 r, [! h& v e(i)=e(i)-DU(2*i);
% |# X1 p' Z% C f(i)=f(i)-DU(2*i-1);
0 W4 S+ s; z% x! c* ~9 e) D. W end
7 S6 ~, j7 C4 Q" O: K% \ E=e+j*f;/ q( U* K s* _6 {9 m6 F
disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');/ ~! J4 I1 l' O5 X( u- t
disp(' ');+ J$ m: g3 c& }+ T- y& v3 {" F
disp('第几次迭代后的近似计算值');
6 J X3 @: h; A disp(Ci);: E9 M" ~/ B) {7 T6 Z
disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
( C* D% o5 E* B& P/ N& S for i1=1:n
" u2 H3 D; X# ^ V(i1)=sqrt(e(i1)^2+f(i1)^2);
# z, W" L6 p. N, n+ e" ? ANGLE(i1)=atan(f(i1)./e(i1))*(180./pi);
- d3 ? Y! P2 A% \9 P end
5 c9 g; C: `7 h7 E+ I; E1 ~/ R%下列输出结果均按节点号从小到大排列
; b9 Z) F, T4 L/ {5 K disp('各个节点的电压实际值为:E=');
/ D$ }$ P/ ^/ C7 ]6 i. p disp(E);
' E6 ^! } R0 ^) I4 J1 ^3 @ disp('各个节点电压幅值为:V=');' g) U6 E1 g8 [# O' F
disp(V);0 q3 z+ B- S! F u
disp('各个节点的电压相交为:ANGLE=');/ S* v2 i. ~" t+ v. Z
disp(ANGLE);
5 }0 l# y. `6 Z2 H" r- }5 B disp('各个节点的功率误差变化');1 K! p, B; x; i8 Q& g! D* ]
disp(J([1:8],9)');5 P) r X" L; Q' ]5 ^5 [/ k* p- l4 a
MDU(Ci)=max(abs(J([1:8],9)'));
& h3 {% G2 K \) ?* R1 ]/ v for p=1:n
/ U( {7 e8 J) H) m- [, I: a1 D Sm(p)=0;
, l t- T1 M3 z: i for q=1:n t- U2 k! _' c- S
Sm(p)=Sm(p)+conj(Y(p,q))*conj(E(q));7 N; t1 }4 V% H0 |5 a4 P
end
% B7 L+ J2 L" K" ]6 }( L' z S(i)=E(p)*Sm(p);. _6 @" O2 U! d Q
end
* g+ Q* l$ P- G2 o7 u! Aend
0 x6 V" r: t- j/ g! \disp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');
; E! X4 R9 B' @' d! Hdisp(' ');3 N8 [& m4 a- |! p- n
disp('计算的最终功率近似值如下各项值');* G; O" z5 {" ?; P* X0 w. ?
disp(' ');; ^5 n5 N% S! a4 t" ]8 c0 D' O; P2 u
disp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @');* a/ Y# r& l. V, |
disp('各个节点功率为:S=');3 h4 E% Q, ]1 P
disp(S);
6 ^* U4 O3 P; q7 ndisp('各条支路的首端功率Si为:');
4 @' v3 l( a! K/ h1 Q A2 _for i=1:nl
: i: i, Q- z! e. w( k+ N p=Br(i,1);q=Br(i,2);4 K5 r% U6 T B
A(i)=E(p)*[conj(E(p))-conj(E(q))]*conj(1./(Br(i,3)));7 |+ V+ M) Q. C7 g3 X% h2 Z4 q
SBS(i)=[abs(E(p))]^2*conj(Br(i,4))+A(i); 6 b3 f' q2 Q; |
end. {+ g: z! _: G8 o7 ^
disp(SBS);0 U) U+ k' _/ C
disp('各条支路的末端功率Sj为:');
- y* ~& d% \8 {- W+ Kfor i=1:nl" N( V; N& n8 E! q/ x
p=Br(i,1);q=Br(i,2);+ \9 }. M6 M: b% P
A2=E(q)*[conj(E(q))-conj(E(p))]*conj(1./(Br(i,3)));7 n, X9 x: c( d0 ^+ Y" d) V( v0 T
SBE(i)=[abs(E(q))]^2*conj(Br(i,4))+A2;
, _+ j, \; E$ G5 d$ p3 |% Dend
9 ~# N; P' j3 {# b6 V) d& W b% U: Jdisp(SBE);
& I( b: x0 p6 A) ?/ N3 }8 Kdisp('各条支路的功率损耗为:');5 A1 ]9 a0 ?, H P* m% _" d x4 T
for i=1:nl
/ t7 ?: O7 f, v* w7 d p=Br(i,1);q=Br(i,2);2 ]( T% y7 _: r2 p
SBL(i)=SBS(i)+SBE(i);" Z2 _+ m+ D$ E& w9 e( Y) I
end9 O- d9 s4 Y t+ E
disp(SBL);
" ^& [' t; F$ H$ X; n2 w%绘制功率误差曲线,先用最小二乘法进行曲线拟合
5 K1 j; F5 Q8 a7 q( p' }for i=1:Ci
/ U( C8 @2 g8 A4 L" j1 Y) I$ G CSH(i)=i;! r( @& a2 V- i {0 v
end0 G: F0 x \6 ~% X; a
P=polyfit(CSH,MDU,3);1 R9 A1 ~8 j8 c9 p7 n( s& W
plot(CSH,MDU,'-O'); |
|