牛拉法潮流程序

星夜飘雪

希望对大家有所帮助哦

1. clc
   
2. disp('此程序为牛拉法解潮流')
   
3. nPQ=input('请输入PQ节点的个数：');
   2 @  P( K6 ?3 a' o1 K8 M- n; f  ^  Y
4. nPV=input('请输入PV节点的个数：');
   * L8 q( }/ Q! e, q3 L0 C* q( o
5. n=nPQ+nPV+1;
   " f; J8 r9 U8 l: u8 [
6. Ps=[0;-0.5;0.2];
   & i& l0 \# N0 c8 X4 H1 L0 q6 {
7. Qs=[0;-0.3];
   
8. Us=[1.0+j*0;1.0+j*0;1.05+j*0;1.05+j*0];
   
9. % nl nr R X Bl Br
   
10. zdata=[1 2 0 0.1880 -0.6815 0.6040;
    
11. 1 3 0.1302 0.2479 0.0129 0.0129;
    % C; z# R( z/ F) q7 Z3 f
12. 1 4 0.1736 0.3306 0.0172 0.0172;
    * ^, }+ I, q3 ^2 W0 R. ]
13. 3 4 0.2603 0.4959 0.0259 0.0259];
    
14. % nx B
    
15. xdata=[2 0.05];
    7 W8 w4 \9 N6 G: j" C
16. dPQU=1;
    
17. %计算导纳矩阵
    0 N9 L- [% Z7 y9 [' b/ ~, k
18. nl=zdata(:,1);
    ; }# o( S6 Z1 Q% u3 Z! m) Z
19. nr=zdata(:,2);
    
20. R=zdata(:,3);
    
21. X=zdata(:,4);
    ( ?4 ?4 X5 j5 F6 ~& j
22. Bl=zdata(:,5);
    : O0 }" g/ v: Q9 \4 _9 D4 D
23. Br=zdata(:,6);
    9 ]! [: Q  s) V
24. nx=xdata(:,1);
    
25. Bx=xdata(:,2);
    
26. nbr=length(zdata(:,1));
    
27. nbrx=length(xdata(:,1));
    4 b" Y1 S/ L6 [9 B
28. Z=R+j*X;
    - u+ M) O9 w" }5 }! E; v( X6 J9 V
29. y=ones(nbr,1)./Z;
    
30. Y=zeros(n,n);
    ' I( v1 M4 p- e& V' j$ ]2 Q2 H2 Q
31. %计算非对角元素
    * N) u8 [$ q1 j+ Y
32. for ii=1:nbr
    ( E) a9 m$ w3 e0 ^
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. %计算对角元素
    1 h; a$ \6 B+ u; [9 |: z& N
37. for ii=1:n
    " R) a- m% S& f5 m# N. x5 S+ Y8 |
38. for jj=1:nbr
    
39. if nl(jj)==ii|nr(jj)==ii
    
40. Y(ii,ii)=Y(ii,ii)+y(jj);
    
41. end
    
42. end
    
43. end
    " D9 _  Z# s2 v- i
44. for ii=1:nbr
    
45. Y(nl(ii),nl(ii))= Y(nl(ii),nl(ii))+j*Bl(ii);
    ' r6 H: a- [2 @
46. Y(nr(ii),nr(ii))= Y(nr(ii),nr(ii))+j*Br(ii);
    & [/ e* {9 f* R* k: S3 a8 A
47. end
    ! o4 u9 o( @$ y; X. r2 n, S
48. for ii=1:nbrx
    
49. Y(nx(ii),nx(ii))=Y(nx(ii),nx(ii))+j*Bx(ii);
    
50. end
    
51. %分离G、B
    4 w) ^% a( C$ k( V6 p4 {: |8 y
52. G=real(Y);
    
53. B=imag(Y);
    7 h( J/ s; Q& L
54. disp('导纳矩阵：');
    
55. Y
    
56. e=real(Us);
      e; _' D5 A8 G( ?: K* t2 V; A
57. f=imag(Us);
    ; f" q6 w- C7 a* k
58. k=0;
    ' L0 ^4 d7 r5 `9 x/ e; D3 v' R
59. while dPQU>0.00001
    
60. %求dP
    
61. dP=zeros(n-1,1);
    3 I1 o, _+ l, a4 X/ A( Q7 c
62. for ii=1:n-1
    
63. t=0;
    8 K' y/ ^6 d3 h/ s' [( Q- J1 K5 {
64. for jj=1:n
    
65. t=t+conj(Y(ii,jj))*(e(jj)-j*f(jj));
    
66. end
    
67. dP(ii)=Ps(ii)-real(t*(e(ii)+j*f(ii)));
    5 {$ Y" U3 Z. ]- {* x" s4 G
68. end
    
69. %求dQ
    4 v9 H# a( |; E
70. dQ=zeros(nPQ,1);
    " A: w& \9 E8 d, ?1 G2 c' ?
71. for ii=1:nPQ
    9 z/ Q$ I  ?8 L5 A7 w
72. t=0;
    " i$ I7 G2 E: j/ q. }: m; S
73. for jj=1:n
    
74. t=t+conj(Y(ii,jj))*(e(jj)-j*f(jj));
    
75. end
    - ~$ M: `. E  L! ~  p/ S
76. dQ(ii)=Qs(ii)-imag(t*(e(ii)+j*f(ii)));
    
77. end
    3 d+ ~* U% f1 B# n5 A( i
78. %求dU^2
    
79. dU2=zeros(nPV,1);
    
80. ii=1:nPV;
    
81. dU2(ii)=abs(Us(ii+nPQ))^2-abs(e(ii+nPQ)+j*f(ii+nPQ))^2;
    1 j1 ]' I5 r. Y% m) U$ g
82. dS=[dP;dQ;dU2];
    
83. dPQU=max(abs(dS));
    
84. if(dPQU>0.00001)
    3 l3 K. P0 C+ {- k0 Y
85. k=k+1
    
86. %形成雅克比行列式
    5 U) x* O3 o' E+ Z
87. Jacob=zeros(2*(n-1),2*(n-1));
    ( ?" A5 Q4 i: u2 k5 l$ @- O
88. %P部分
    
89. for ii=1:n-1
    
90. mid1=0;
    # \* r5 i% W6 a! Q* m
91. mid2=0;
    
92. for jj=1:n
    ; G% \: o# E' v/ f4 I9 P9 f
93. if ii~=jj&&jj<n
    - R" k, ^! z- u/ k
94. Jacob(ii,2*jj-1)=-(G(ii,jj)*e(ii)+B(ii,jj)*f(ii));
    
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);
    ! m0 i% ^! h4 F& O& K
98. mid2=mid2+G(ii,jj)*e(jj)-B(ii,jj)*f(jj);
    / u% `/ ?9 z% N8 S4 Y/ o6 h
99. end
    % A& X. y" {& V& y
100. Jacob(ii,2*ii-1)=-mid2-G(ii,ii)*e(ii)-B(ii,ii)*f(ii);
     
101. Jacob(ii,2*ii)=-mid1+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);
     ) D' r# x- p" ?; h  r& ^
102. end
     ' X8 I* f! B1 [
103. %Q部分
     0 i9 `0 j7 f4 f6 @0 j
104. for ii=1:nPQ
     ; N( O; v; A, A9 a* J, I
105. mid1=0;
     5 [) `) q  d1 w( g" {5 y
106. mid2=0;
     
107. for jj=1:n
     
108. if ii~=jj&&jj<n
     
109. Jacob(ii+n-1,2*jj-1)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);
     ; H$ F" i* g6 A( `" @
110. Jacob(ii+n-1,2*jj)=G(ii,jj)*e(ii)+B(ii,jj)*f(ii);
     ) j+ ^4 u' J( D8 H2 S
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);
     
114. end
     $ {9 O" q8 q' }; d; s& k$ B
115. Jacob(ii+n-1,2*ii-1)=mid1+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);
     ; w  s- p, b' c  a; e* l" a8 J6 T) i2 A
116. Jacob(ii+n-1,2*ii)=-mid2+G(ii,ii)*e(ii)+B(ii,ii)*f(ii);
     
117. end
     
118. %U2部分
     $ r( |! m' @! Y" _- D' B% J, ~
119. for ii=nPQ+1:n-1
     5 N- j: M. Q) b+ D
120. Jacob(ii+n-1,2*ii-1)=-2*e(ii);
     4 H: @; `! `" k' N0 t( f/ U# P
121. Jacob(ii+n-1,2*ii)=-2*f(ii);
     
122. end
     
123. dU=-inv(Jacob)*dS;
     
124. de=zeros(n-1,1);
     ; I" {. T6 g" v5 }; j. z
125. df=zeros(n-1,1);
     
126. ii=1:n-1;
     
127. de(ii)=dU(2*ii-1);
     
128. df(ii)=dU(2*ii);
     
129. e(ii)=e(ii)+de(ii);
     / F7 t8 X, W# N4 f& E1 |- s3 h& ?4 f4 [; R
130. f(ii)=f(ii)+df(ii);
     9 O6 _1 k- G) ~4 m! j  ^$ Q
131. end
     9 C! A4 F& d0 e2 [# Q3 `( n  U
132. end%迭代结束
     1 S) X1 j7 f- W0 }
133. U=e+j*f;
     : t! m6 `* W: P; D4 w! g; ?5 V
134. %计算PV节点的Q
     
135. P=zeros(n,1);
     . s6 R: J$ }* D8 j& I" F
136. Q=zeros(n,1);
     
137. for ii=1:nPV
     
138. t=0;
     0 n2 E6 F/ r  v8 X) w. s. j
139. for jj=1:n
     
140. t=t+conj(Y(ii+nPQ,jj))*(e(jj)-j*f(jj));
     5 O3 {. Q+ X6 v  J3 _
141. end
     . f* ~9 l! R4 f- M7 U2 ?- Y
142. Q(ii+nPQ)=imag(t*(e(ii+nPQ)+j*f(ii+nPQ)));
     9 t8 }8 T+ c& d8 I3 U4 G3 |
143. end
     
144. %计算平衡节点
     
145. t=0;
     * C9 }8 N/ T7 f: |; o2 H2 W9 Q
146. for jj=1:n
     
147. t=t+conj(Y(n,jj))*(e(jj)-j*f(jj));
     
148. end
     
149. P(n)=real(t*(e(n)+j*f(n)));
     / Z9 [1 X8 p) \- Y7 n( L' j
150. Q(n)=imag(t*(e(n)+j*f(n)));
     ! j9 q8 F' P6 G: Q+ {
151. ii=1:n-1;
     
152. P(ii)=Ps(ii);
     - u* Q3 m) P$ D2 A' _7 ^. W/ i( q4 G
153. ii=1:nPQ;
     
154. Q(ii)=Qs(ii);
     6 [, w* w! [( C& W
155. %计算线路潮流
     ) k: M& S% l$ I! P9 G5 ]
156. Sij=zeros(nbr,1);
     
157. Sji=zeros(nbr,1);
     $ @% X5 r* J0 A% _8 y' d
158. dSij=zeros(nbr,1);
       p2 _+ x% v! N! n: {
159. for ii=1:nbr
     5 d/ c# V$ \+ H* |; l: b
160. Sij(ii)=U(nl(ii))*(conj(U(nl(ii)))*(-j*Bl(ii))+(conj(U((nl(ii))))-conj(U((nr(ii)))))*conj(y(ii)));
     
161. Sji(ii)=U(nr(ii))*(conj(U(nr(ii)))*(-j*Br(ii))+(conj(U((nr(ii))))-conj(U((nl(ii)))))*conj(y(ii)));
     , e8 O# c8 t; |, y8 O1 L
162. dSij(ii)=Sij(ii)+Sji(ii);
     
163. end
     
164. nn=[1:n]';
     
165. disp(' n e f P Q');
     
166. Display1=[nn e f P Q]
     9 @/ S8 h0 H3 N' D) r
167. disp(' nl nr Sij Sji dSij');
     $ I* m8 [/ r- N0 B& c- @$ B* R
168. Display2=[nl nr Sij Sji dSij]
     . v" H) \8 m* [% h
169. 
     ( B7 ?3 K, s: V% i) k& b6 B) z

复制代码

leichen

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

wangwangmoon

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

花尹儿

正在学习中

云来123

谢谢分享

风晨

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