牛拉法潮流程序

星夜飘雪

希望对大家有所帮助哦
2 |( q/ n. \* C& G

1. clc
   
2. disp('此程序为牛拉法解潮流')
   
3. nPQ=input('请输入PQ节点的个数：');
   
4. nPV=input('请输入PV节点的个数：');
   
5. n=nPQ+nPV+1;
   
6. Ps=[0;-0.5;0.2];
   " g2 ^. A% Z8 t, w- E. h. y
7. Qs=[0;-0.3];
   $ V/ i% Z. {  e9 X6 j
8. Us=[1.0+j*0;1.0+j*0;1.05+j*0;1.05+j*0];
   3 r! d. H6 Z/ \% y
9. % nl nr R X Bl Br
   9 [. U% @6 H' n- R$ F, w
10. zdata=[1 2 0 0.1880 -0.6815 0.6040;
    
11. 1 3 0.1302 0.2479 0.0129 0.0129;
    ' F+ E: a3 ]7 b% g& k) l, ^
12. 1 4 0.1736 0.3306 0.0172 0.0172;
    
13. 3 4 0.2603 0.4959 0.0259 0.0259];
    0 e0 ~+ z; b" r4 V
14. % nx B
    
15. xdata=[2 0.05];
    3 e1 S2 M9 S; f! F
16. dPQU=1;
    
17. %计算导纳矩阵
    
18. nl=zdata(:,1);
    7 d" E" \% d; _
19. nr=zdata(:,2);
    
20. R=zdata(:,3);
    - z/ C4 ^1 ^% X# h% g; M
21. X=zdata(:,4);
    5 D- W, k  b5 V: X1 ]! I
22. Bl=zdata(:,5);
    $ A4 I# ]. e: I; K) d
23. Br=zdata(:,6);
    
24. nx=xdata(:,1);
    9 P5 p- S" R; i) [/ a
25. Bx=xdata(:,2);
    , M9 ^7 _" d0 S, g; A4 v
26. nbr=length(zdata(:,1));
    
27. nbrx=length(xdata(:,1));
    
28. Z=R+j*X;
    ) j+ @  H  y2 |1 f: \9 ?8 F
29. y=ones(nbr,1)./Z;
    
30. Y=zeros(n,n);
    ; K) \# i" P' b3 @' Y
31. %计算非对角元素
    
32. for ii=1:nbr
    4 K0 D2 ]7 L4 x1 v/ ?
33. Y(nl(ii),nr(ii))= Y(nl(ii),nr(ii))-y(ii);
    
34. Y(nr(ii),nl(ii))= Y(nl(ii),nr(ii));
    
35. end
    
36. %计算对角元素
    
37. for ii=1:n
    . x+ L1 ]0 D6 ~* n
38. for jj=1:nbr
    : Z- Z: J6 Y* m! q, U5 N
39. if nl(jj)==ii|nr(jj)==ii
    + G% D" A4 }" t" `; J
40. Y(ii,ii)=Y(ii,ii)+y(jj);
    7 s) y* u) e/ N& V. n) x* ~2 b4 E
41. end
    
42. end
    
43. end
    
44. for ii=1:nbr
    9 G' c; C, y  G9 }7 p* A# F& P; P9 `
45. Y(nl(ii),nl(ii))= Y(nl(ii),nl(ii))+j*Bl(ii);
    
46. Y(nr(ii),nr(ii))= Y(nr(ii),nr(ii))+j*Br(ii);
    
47. end
    
48. for ii=1:nbrx
    7 e% {, H. \/ z0 r2 m
49. Y(nx(ii),nx(ii))=Y(nx(ii),nx(ii))+j*Bx(ii);
    : g0 }" f0 J( @4 n- G/ M
50. end
    4 j9 |& u/ x9 w7 Z3 u) h
51. %分离G、B
    
52. G=real(Y);
    ( e! q, |: q+ |& t3 m$ X5 ]
53. B=imag(Y);
    
54. disp('导纳矩阵：');
    5 O3 |& V- [% }/ I! n" d0 g
55. Y
    0 |- m# Q1 G( P3 O
56. e=real(Us);
    $ M; o  A) K3 s
57. f=imag(Us);
    
58. k=0;
    
59. while dPQU>0.00001
    , q* o7 E) t+ L3 a
60. %求dP
    
61. dP=zeros(n-1,1);
    
62. for ii=1:n-1
    
63. t=0;
    
64. for jj=1:n
    ' k6 E; A$ \( W( x7 |6 p
65. t=t+conj(Y(ii,jj))*(e(jj)-j*f(jj));
    " Q4 R+ [8 e: v1 h" G
66. end
    ' x" h, x4 P: R# d+ k. j# g% q
67. dP(ii)=Ps(ii)-real(t*(e(ii)+j*f(ii)));
    
68. end
    
69. %求dQ
    7 r, `  n: v3 I( T) O% i& e" v
70. dQ=zeros(nPQ,1);
    
71. for ii=1:nPQ
    2 k) q# G- W; E5 Z$ v+ M
72. t=0;
    9 h" [& Z8 r7 J9 T) Y0 V+ x
73. for jj=1:n
    
74. t=t+conj(Y(ii,jj))*(e(jj)-j*f(jj));
    # b" a9 _7 b5 @) s) C8 b
75. end
    - Y+ S8 E" C" o+ @" i- f9 J
76. dQ(ii)=Qs(ii)-imag(t*(e(ii)+j*f(ii)));
    
77. end
    ) m+ v# U; M4 p/ |, P1 D
78. %求dU^2
    
79. dU2=zeros(nPV,1);
    # U1 {7 B. T; M: S1 o% u$ J
80. ii=1:nPV;
    
81. dU2(ii)=abs(Us(ii+nPQ))^2-abs(e(ii+nPQ)+j*f(ii+nPQ))^2;
    ; P" Z/ Y; N* V" B
82. dS=[dP;dQ;dU2];
    
83. dPQU=max(abs(dS));
    
84. if(dPQU>0.00001)
    7 i0 N* |* C+ ~3 M& U
85. k=k+1
    - C3 W. F* ?" C- \+ Y" _: B
86. %形成雅克比行列式
    ; a* [9 P& l6 V" g$ d( [
87. Jacob=zeros(2*(n-1),2*(n-1));
    6 b3 _6 S. w1 f( ~& }# {3 o
88. %P部分
    6 d/ @! ]: z! a, M, o
89. for ii=1:n-1
    
90. mid1=0;
    : d, G( r, T9 `; `/ s; ~
91. mid2=0;
    
92. for jj=1:n
    
93. if ii~=jj&&jj<n
    * l* a* D0 D, V5 A0 x
94. Jacob(ii,2*jj-1)=-(G(ii,jj)*e(ii)+B(ii,jj)*f(ii));
    7 i3 M! z6 q) \& M% M6 Z  y+ G
95. Jacob(ii,2*jj)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);
    
96. end
    
97. mid1=mid1+G(ii,jj)*f(jj)+B(ii,jj)*e(jj);
    
98. mid2=mid2+G(ii,jj)*e(jj)-B(ii,jj)*f(jj);
    ) @0 u! _6 M  ~) o
99. end
    8 K* o" f. x/ P4 I
100. Jacob(ii,2*ii-1)=-mid2-G(ii,ii)*e(ii)-B(ii,ii)*f(ii);
     2 ?; s3 k: N' p7 q2 H* G% i8 n1 d
101. Jacob(ii,2*ii)=-mid1+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);
     ' H% y! B% j$ a0 o
102. end
     6 K/ g8 G3 e1 D$ @& q: A, i4 z
103. %Q部分
     
104. for ii=1:nPQ
     
105. mid1=0;
     / g- X( K! T* a) Z% T
106. mid2=0;
     
107. for jj=1:n
     ! i: O- y! w1 M* T6 t, u1 G
108. if ii~=jj&&jj<n
     9 y1 z5 o5 Q( X7 ]
109. Jacob(ii+n-1,2*jj-1)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);
     
110. Jacob(ii+n-1,2*jj)=G(ii,jj)*e(ii)+B(ii,jj)*f(ii);
     : Z3 T: n" Y- e( Y
111. end
     
112. mid1=mid1+G(ii,jj)*f(jj)+B(ii,jj)*e(jj);
     
113. mid2=mid2+G(ii,jj)*e(jj)-B(ii,jj)*f(jj);
     5 i2 L3 T; i! b' x( Y) K
114. end
     
115. Jacob(ii+n-1,2*ii-1)=mid1+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);
     3 C% N" F% d, o, o2 u, K0 V. O
116. Jacob(ii+n-1,2*ii)=-mid2+G(ii,ii)*e(ii)+B(ii,ii)*f(ii);
     + h' g4 i. X# e8 Q  D& x( u" K7 _/ c
117. end
     ; m& ~+ O& q% \+ F6 ^
118. %U2部分
     
119. for ii=nPQ+1:n-1
     8 }+ x8 M+ }7 d
120. Jacob(ii+n-1,2*ii-1)=-2*e(ii);
     : _5 [+ T( g4 a4 a" K0 k
121. Jacob(ii+n-1,2*ii)=-2*f(ii);
     
122. end
     " j3 u; l2 E* f& R1 s# @
123. dU=-inv(Jacob)*dS;
     0 `8 p+ q0 ^* Y
124. de=zeros(n-1,1);
     * X, {+ @: C1 i, R: ?8 }+ R9 k
125. df=zeros(n-1,1);
     
126. ii=1:n-1;
     0 t9 z0 Z5 q* o: D, ]  X) I0 q+ e0 E8 W
127. de(ii)=dU(2*ii-1);
     % K  d. q9 F" p4 l; r9 _8 w
128. df(ii)=dU(2*ii);
     : a3 {( \5 X! o1 t4 t
129. e(ii)=e(ii)+de(ii);
     
130. f(ii)=f(ii)+df(ii);
     
131. end
     
132. end%迭代结束
     
133. U=e+j*f;
     , O+ _0 b. |( N
134. %计算PV节点的Q
     / ^! O+ _9 P4 j* A3 I2 D9 n) a& k- |
135. P=zeros(n,1);
     
136. Q=zeros(n,1);
     
137. for ii=1:nPV
     # p5 `8 A  u# @
138. t=0;
     7 ~8 C* ^( I8 T% F1 D: S# q
139. for jj=1:n
     
140. t=t+conj(Y(ii+nPQ,jj))*(e(jj)-j*f(jj));
     ; m; c5 z5 J1 v: x# Y/ ^
141. end
     
142. Q(ii+nPQ)=imag(t*(e(ii+nPQ)+j*f(ii+nPQ)));
     % s; r$ B) y. a( ?3 c8 ]
143. end
     0 Y2 D: d. G# E
144. %计算平衡节点
     " A9 U1 j. n: W6 h
145. t=0;
     1 u* g4 r  `% m) \6 T( f8 s! l" r
146. for jj=1:n
     
147. t=t+conj(Y(n,jj))*(e(jj)-j*f(jj));
     ! c* G, H" l/ \2 W, Z' V
148. end
     3 x! D/ f0 ^# k, G9 L
149. P(n)=real(t*(e(n)+j*f(n)));
     
150. Q(n)=imag(t*(e(n)+j*f(n)));
     
151. ii=1:n-1;
     
152. P(ii)=Ps(ii);
     
153. ii=1:nPQ;
     - Y$ m4 W5 E) u# Q7 U+ D
154. Q(ii)=Qs(ii);
     9 Q4 D. X& C+ f; V9 o: N$ A# S: a
155. %计算线路潮流
     8 Q! X2 x0 O- g3 O/ i
156. Sij=zeros(nbr,1);
     
157. Sji=zeros(nbr,1);
     * m) Z% @. r& u2 R" ]+ U
158. dSij=zeros(nbr,1);
     
159. for ii=1:nbr
     
160. Sij(ii)=U(nl(ii))*(conj(U(nl(ii)))*(-j*Bl(ii))+(conj(U((nl(ii))))-conj(U((nr(ii)))))*conj(y(ii)));
     - d  ]/ K" s$ \+ f, l# _) X
161. Sji(ii)=U(nr(ii))*(conj(U(nr(ii)))*(-j*Br(ii))+(conj(U((nr(ii))))-conj(U((nl(ii)))))*conj(y(ii)));
     
162. dSij(ii)=Sij(ii)+Sji(ii);
     8 {+ Y6 n0 U3 `# b* d
163. end
     0 }  m% F, v( @# s5 I- J
164. nn=[1:n]';
     ; i8 A' l5 T5 M$ ^0 H
165. disp(' n e f P Q');
     
166. Display1=[nn e f P Q]
     
167. disp(' nl nr Sij Sji dSij');
     3 A4 F9 Q6 }7 H5 e' f" r$ l2 i% v/ }
168. Display2=[nl nr Sij Sji dSij]
     
169. 
     

复制代码

leichen

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

wangwangmoon

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

花尹儿

正在学习中

云来123

谢谢分享

风晨

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