牛拉法潮流程序

星夜飘雪 · 发表于 2010-7-15 09:31:48

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

您需要登录才可以下载或查看，没有账号？立即加入

×

希望对大家有所帮助哦

clc0 l. C4 I+ k) [( Z! ?/ G
disp('此程序为牛拉法解潮流')/ w+ |; A3 E8 V1 V
nPQ=input('请输入PQ节点的个数：');
) F/ _% n. [7 d! L
nPV=input('请输入PV节点的个数：');
6 m3 M; M7 B( @) E) t4 m
n=nPQ+nPV+1;% w2 \$ E8 ?; l/ v0 @: k0 ~2 d6 h
Ps=[0;-0.5;0.2];. M4 A/ H/ R7 n7 m1 A& l
Qs=[0;-0.3];% a* z: Y/ N: Q$ R
Us=[1.0+j*0;1.0+j*0;1.05+j*0;1.05+j*0];
% {9 W* x3 e5 [, [ `" p7 E( X
% nl nr R X Bl Br& Z' ^1 P- z6 g( I, F; d+ |0 |3 E- y% I
zdata=[1 2 0 0.1880 -0.6815 0.6040;- W: ~1 v4 W3 X( s
1 3 0.1302 0.2479 0.0129 0.0129;$ j' [5 M6 R# Y
1 4 0.1736 0.3306 0.0172 0.0172;
: Q, `( }5 i1 o( d
3 4 0.2603 0.4959 0.0259 0.0259];
2 R" V9 o/ y y8 f! L# ]
% nx B ( r6 @$ K7 ?7 |) a. N( ^
xdata=[2 0.05]; : d4 w( P' E# g- [
dPQU=1;
) D9 o; ~/ ^: R+ V0 {
%计算导纳矩阵
2 n4 a2 w0 v$ b0 U5 S9 J2 {
nl=zdata(:,1);
0 t" M9 D0 Z( s% A2 a: T2 A
nr=zdata(:,2);
j6 K+ r" t% }4 m# N* T
R=zdata(:,3);* T! a6 R! H; i" Z/ `. q$ ?
X=zdata(:,4);5 }+ v0 h, v9 f/ r& ^
Bl=zdata(:,5);
|1 `& Q3 ]- O3 x6 n6 W0 Z3 j6 C
Br=zdata(:,6);
?1 i4 L Q5 h
nx=xdata(:,1);; c) z2 s3 [1 j1 g
Bx=xdata(:,2);' X9 r$ A: w. Z ]" ^* j9 L. d7 x
nbr=length(zdata(:,1));
nbrx=length(xdata(:,1));
5 X6 G' o# k+ A# W' p
Z=R+j*X;
2 x; W+ S s$ i2 S1 x, T* ]; N
y=ones(nbr,1)./Z;
0 M" w, }$ a+ t6 p2 ^; V
Y=zeros(n,n);
- q! x" |. w* ?5 N- H
%计算非对角元素
' ?6 z/ \$ ?$ W, x& T5 X g
for ii=1:nbr
5 z: W& D% d9 B' [$ a# S; z0 l2 C1 ~
Y(nl(ii),nr(ii))= Y(nl(ii),nr(ii))-y(ii);
4 H0 e8 ?/ b, X
Y(nr(ii),nl(ii))= Y(nl(ii),nr(ii));
5 E9 ]* S8 ?( A! @. s
end
6 O5 S! @; U* ?" H+ c% W
%计算对角元素' P8 [8 B# }/ J* B" M
for ii=1:n
* A; a9 Q. Z3 J& D
for jj=1:nbr
8 I5 `+ F8 \' g% Y* M8 s1 x
if nl(jj)==ii|nr(jj)==ii7 g* W2 |% E4 j; G" g& Y& Q
Y(ii,ii)=Y(ii,ii)+y(jj);2 [# b" r/ d j2 L
end1 ^; { v8 w2 P) ]5 j( l$ ?; u
end* R: W$ s5 x3 S+ v* o g
end
) [# M4 x# M* ~: t& T$ a) K3 [
for ii=1:nbr6 y. u, V6 f- f* e5 [7 _" w
Y(nl(ii),nl(ii))= Y(nl(ii),nl(ii))+j*Bl(ii);
3 V- F# `/ t" h0 a
Y(nr(ii),nr(ii))= Y(nr(ii),nr(ii))+j*Br(ii);
7 l% d6 \6 Z3 d* y! F7 D% d" z
end
. _( \/ t4 R1 A: j/ r
for ii=1:nbrx# e* n& L" [/ R! ]# H6 R6 S
Y(nx(ii),nx(ii))=Y(nx(ii),nx(ii))+j*Bx(ii);5 Z+ q7 v- b+ V
end- n; {) w7 `, t
%分离G、B
I% Z2 z3 P. G' q. P! x8 l& q
G=real(Y);. p3 @8 [; n' C8 P! |6 d0 ]6 F
B=imag(Y);! r) w3 u6 z# C: u
disp('导纳矩阵：');; U* h# b" i8 o% D" j/ Y
Y6 M/ G( a% i2 l e# |& L
e=real(Us);
7 n+ a) ]/ @ T, i! f0 U7 I5 c; @
f=imag(Us);
- h. h) O/ L) @. w7 L; C, ~
k=0;
4 W. Z/ w) o6 u' k$ A: f* E2 J
while dPQU>0.00001 ( F) b N$ Z" }: Q; W# B
%求dP9 d5 C/ ^( {9 f H1 R8 ^
dP=zeros(n-1,1);: p: e4 [, S/ l ?) g
for ii=1:n-17 k. Z* }4 w0 k! o
t=0; " R3 c M3 e4 {
for jj=1:n
! N. c; s0 F" \7 m/ B+ t5 d
t=t+conj(Y(ii,jj))*(e(jj)-j*f(jj)); 9 u/ G5 f6 z2 U* ]3 s' \; T; ?
end
9 |, I3 x" t, g' ]" B$ o
dP(ii)=Ps(ii)-real(t*(e(ii)+j*f(ii))); ; S9 w0 f" t: f% m. c
end
, I; g) x1 Z) t6 C* G
%求dQ
" |0 {* b5 W# W P; B9 s7 z
dQ=zeros(nPQ,1);
9 N4 |9 w% T$ Y4 b
for ii=1:nPQ
* \2 ]# b$ c( f2 K1 [
t=0;
2 ~" C( s4 N9 A; R- n* t9 G
for jj=1:n
/ \# R' O; R7 B# x6 g# t
t=t+conj(Y(ii,jj))*(e(jj)-j*f(jj));( \" @' h6 K) @6 M2 N& X. }
end" x9 A; s& \8 a5 I* w/ Q2 j- ^0 P
dQ(ii)=Qs(ii)-imag(t*(e(ii)+j*f(ii)));
1 ~, E# _& a$ V( Z) k9 u
end
2 h& V" o' M* g( e6 P; d7 S, Z
%求dU^2- K9 t. c- Z, ^* w8 d
dU2=zeros(nPV,1);
6 l9 `% b5 _4 _: o- n
ii=1:nPV;
5 i7 D: |7 L' @# t& |
dU2(ii)=abs(Us(ii+nPQ))^2-abs(e(ii+nPQ)+j*f(ii+nPQ))^2;8 m7 q) Q" M* o: Z- a. C2 h
dS=[dP;dQ;dU2];
" u8 d: Q' P9 Q( D& K" C
dPQU=max(abs(dS));& P: t( ^/ A6 o; _# _( b
if(dPQU>0.00001)
- b( Q `' m5 z' ~0 L+ q( q
k=k+1% E$ x j8 k% u- i# G* j
%形成雅克比行列式: j7 W( p2 z! u: E* r
Jacob=zeros(2*(n-1),2*(n-1));- e+ [' b5 i1 q+ n6 J/ `" S' C% L
%P部分
/ }7 I, W; W- H3 ^
for ii=1:n-14 ?0 n& s" @2 }: N4 z: o H
mid1=0;
: }# ]5 T& A( S. l* ^. u! N
mid2=0;2 N2 W3 }. r# J4 V' t m# `
for jj=1:n
& n) M: i! R4 F1 s- j3 j
if ii~=jj&&jj<n
K9 \! k# g: A1 t" V/ P$ ?
Jacob(ii,2*jj-1)=-(G(ii,jj)*e(ii)+B(ii,jj)*f(ii));* J9 i$ J6 u9 @% ]! r; H" ]% n8 h
Jacob(ii,2*jj)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);, v$ X) g' @" ~5 p' R* M# M
end
1 j/ ]8 r9 G" ?' q5 h
mid1=mid1+G(ii,jj)*f(jj)+B(ii,jj)*e(jj);! \( F, Q4 Q9 V7 F/ t
mid2=mid2+G(ii,jj)*e(jj)-B(ii,jj)*f(jj);
3 Y( O! T3 u+ M* Z) G( F
end: s( G: \! `- A4 ~% ]' U
Jacob(ii,2*ii-1)=-mid2-G(ii,ii)*e(ii)-B(ii,ii)*f(ii);& R! `& F1 M. x" @
Jacob(ii,2*ii)=-mid1+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);8 f- y6 Y* y. n
end
5 @3 U3 `9 _8 u
%Q部分
/ ]' @. e. L1 L$ M5 x8 l
for ii=1:nPQ- U& g. J; ` d4 A
mid1=0;4 o$ j' S% z+ h& Q1 x9 Q4 |* k7 p
mid2=0;* l z6 f: C, r0 @( i! S4 z3 P6 m
for jj=1:n( z$ Z. n J5 |9 J- r
if ii~=jj&&jj<n+ P6 x8 F( [1 \$ |" R' M3 u
Jacob(ii+n-1,2*jj-1)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);/ f/ W* f2 M3 @, `' G* d! N
Jacob(ii+n-1,2*jj)=G(ii,jj)*e(ii)+B(ii,jj)*f(ii);5 c8 k* N5 ~8 y+ {1 p
end
- N! k; N7 t0 W3 G i
mid1=mid1+G(ii,jj)*f(jj)+B(ii,jj)*e(jj);
( }# `1 d& h w9 Q
mid2=mid2+G(ii,jj)*e(jj)-B(ii,jj)*f(jj);
9 J6 K. ~8 h4 C4 A) s- U4 T
end7 t, d2 d" A* G9 ]; V( |
Jacob(ii+n-1,2*ii-1)=mid1+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);1 W4 @3 X1 }% i$ l1 B
Jacob(ii+n-1,2*ii)=-mid2+G(ii,ii)*e(ii)+B(ii,ii)*f(ii);! p( ]! a4 x# Q( o0 I0 j
end
7 d7 S% Q* M+ E) U0 a: o
%U2部分
$ s6 \6 f& u, N' w9 F3 A- B* t2 J
for ii=nPQ+1:n-1; q- a% u: @! E m
Jacob(ii+n-1,2*ii-1)=-2*e(ii);
( W1 S" M+ A" Q8 ^$ Q5 k" j/ Q
Jacob(ii+n-1,2*ii)=-2*f(ii);
2 N" {4 i+ ?4 B4 n
end; g+ }, x- n2 e2 R, f) ?& c+ N
dU=-inv(Jacob)*dS;" I2 J1 n- c/ U3 C9 s
de=zeros(n-1,1);+ w6 x/ O( q0 W, y# v' X' _$ v5 g
df=zeros(n-1,1);
. H$ i& s# \) v$ ]8 h3 l
ii=1:n-1;9 A* ~1 D; y$ V, @3 B
de(ii)=dU(2*ii-1);& _/ s& z, V9 @) ^6 |% o+ t
df(ii)=dU(2*ii);% S3 }, P' ?2 x) z
e(ii)=e(ii)+de(ii);
: G* J" ~' _) G5 J% O* \
f(ii)=f(ii)+df(ii);! @7 L5 g7 ]+ R
end" F3 h b; ~" K; ]% X6 F
end%迭代结束
5 c; J4 k# r9 B' B2 e
U=e+j*f;% \! t0 g5 x% g, m+ N5 q1 |
%计算PV节点的Q7 t4 B# C$ \. Z: v) e+ A" C$ O
P=zeros(n,1);* O- \: v0 J) R( A, V6 K
Q=zeros(n,1);: C. W, {" w& R- Y1 w
for ii=1:nPV
+ N$ _- t* p: ^ n
t=0;
! K) ~/ m1 i1 c' Q f8 T* M8 B% Y4 g
for jj=1:n( I. i/ t# C) V* m# r
t=t+conj(Y(ii+nPQ,jj))*(e(jj)-j*f(jj));
6 j5 x) z4 [( j5 ^1 F
end
J& r& V6 A* J0 K
Q(ii+nPQ)=imag(t*(e(ii+nPQ)+j*f(ii+nPQ)));2 N% G' C8 T# G- R1 D
end
' P/ P ?: }% G% [
%计算平衡节点
9 n' H$ Q3 k9 o& V" k# _. C; y
t=0; ( E# |8 a7 c+ [# R
for jj=1:n$ V5 r' |( L |$ J) h$ U
t=t+conj(Y(n,jj))*(e(jj)-j*f(jj));
, T* ^) f, [2 L+ w
end
3 G/ A2 B" R/ w# T2 B
P(n)=real(t*(e(n)+j*f(n)));4 B( j1 K0 b! `; U
Q(n)=imag(t*(e(n)+j*f(n)));
# N) }( w$ I2 p4 o1 J: ?
ii=1:n-1;/ s1 b/ d& u5 W3 W. E# h8 n
P(ii)=Ps(ii);* v0 x% Q; l# C0 ?# t; @$ l6 \, V* W
ii=1:nPQ;. L- t! L; s7 [- W
Q(ii)=Qs(ii);+ d% q: B' J, M- e/ D" y$ [3 [
%计算线路潮流. H3 ~5 E+ [6 K) V( _# {
Sij=zeros(nbr,1);& Z$ ]8 J4 a+ e8 z9 C. I4 o+ A( i6 z
Sji=zeros(nbr,1);. f% @+ l$ v% J; z6 R
dSij=zeros(nbr,1);+ I" Q+ p' a" a8 D
for ii=1:nbr2 t0 S9 }( \* B- q
Sij(ii)=U(nl(ii))*(conj(U(nl(ii)))*(-j*Bl(ii))+(conj(U((nl(ii))))-conj(U((nr(ii)))))*conj(y(ii)));
U" I3 D5 Q3 \0 n
Sji(ii)=U(nr(ii))*(conj(U(nr(ii)))*(-j*Br(ii))+(conj(U((nr(ii))))-conj(U((nl(ii)))))*conj(y(ii)));
. a% V3 E t# L% S
dSij(ii)=Sij(ii)+Sji(ii);9 n5 ]* A% I" Z
end
% N, o& \" J: B- ]
nn=[1:n]';
9 U+ k/ ], h9 V) y2 l3 m! _4 }3 `
disp(' n e f P Q');
* R' s- i t1 o" E4 R9 G
Display1=[nn e f P Q]6 Y4 z& o0 E, E! n
disp(' nl nr Sij Sji dSij');
2 ?! [! e% k: [' B7 Q0 i/ J: m
Display2=[nl nr Sij Sji dSij]
& n4 K3 x: M, [- C c
5 ^2 x1 W% f0 x. v4 a: u* r

复制代码

leichen · 发表于 2010-12-11 19:35:42

我正好要用，虽然不知道这个能不能行，先感谢楼主分享啦

wangwangmoon · 发表于 2010-12-24 11:19:59

不行吧，不能保证雅可比矩阵是非奇异矩阵，不能直接求逆吧

花尹儿 · 发表于 2016-3-29 16:46:35

正在学习中

云来123 · 发表于 2016-3-29 17:01:51

谢谢分享

风晨 · 发表于 2018-12-27 08:50:46

原创，永远支持~~~~~~~~~~~~

账号		自动登录	找回密码
密码			立即加入

[分享] 牛拉法潮流程序

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

相关帖子