配电系统三相潮流计算 matlab仿真 结果不正确 求解答

咻咻永远

各位~~我根据下图的一个系统用MATLAB编写了一段程序，是电流型前推回代法的三相配电系统潮流计算，但是现在结果不是我想要的，不收敛。请各位帮我看一下，提提建议。图片见附件了。还有个问题：已知的线路都是三相对称线路的阻抗值，那么线路的阻抗矩阵中的互阻抗应该如何计算，是什么样的形式，应用怎样的公式？这里我不是很清楚。

1. %程序名：qiantuihuitui_I_3.m
   
2. %功能：支路电流前推回推法求解潮流
   
3. clc
   0 S$ `1 e" T+ Q# X
4. clear all;
   4 ^! I* g* Y. Q" g" J
5. %--------------输入网络参数--------------
   0 l* `$ Q$ v3 A4 b
6. %1-支路编号，2-首节点，3-尾节点，4-自阻抗，5-尾节点复功率，6-支路性质（1-馈线段支路，2-变压器支路）,7-尾节点是否带负荷
   0 Q' ?2 w/ B) |+ s8 o0 ?2 d3 g6 R/ d
7. DB=[1 1 2 0.000167+j*0.000208 0.42+j*0.31 1 1
   & Z* v/ T5 M" i' Y9 P. d+ G) X
8.     2 2 3 0.000151+j*0.000188 6.15 1 0
   
9.     3 2 4 0.000066+j*0.000082 0.38+j*0.29 1 1
   
10.     4 2 5 0.000249+j*0.000310 0 1 0
    
11.     5 2 6 0.000172+j*0.000215 0 1 0
    
12.     6 4 7 0.000156+j*0.000195 6.06 1 0
    
13.     7 4 8 0.000162+j*0.000202 6.04 1 0
    6 w6 |% ^3 K* Q/ {3 o
14.     8 4 9 0.000345+j*0.000430 0 1 0
    
15.     9 4 10 0.000287+j*0.000358 0 1 0
    
16.     10  5 11 0.020563+j*0.321594 0 2 0
    8 [, z5 h' F) _, {
17.     11  6 12 0.020563+j*0.321594 0 2 0
      z1 i: O8 o" E& V) Y) g8 Y: G9 z5 r
18.     12  9 13 0.020563+j*0.321594 0 2 0
    
19.     13  10 14 0.020563+j*0.321594 0 2 0  
    
20.     14  11 15 0.000237+j*0.000408 5.72+j*0.12 1 1
    4 V8 n2 E8 L( P" K
21.     15  12 16 0.000237+j*0.000408 5.76+j*0.09 1 1
    " p1 R. h" }, V  y" D+ O
22.     16  13 17 0.000292+j*0.000502 5.86+j*0.11 1 1
    ) e2 h. Z7 h7 E2 r, N4 v
23.     17  14 18 0.000274+j*0.000470 5.81+j*0.14 1 1];
    
24. [n,m]=size(DB);  
    # Y3 l1 k. C3 ]; v' z  |+ v' D
25. B=[1 sin(2*pi/3)+j*cos(2*pi/3) sin(4*pi/3)+j*cos(4*pi/3)];
    , u) @" J/ C! M% d  L& ]
26. C=[1.02 1 1.02 1 1 1 1.02 1.02 1 1 1 1 1 1 1 1 1 1];
    6 M1 F( c  L# Q5 ]
27. U(:,1)=B(1,1)*C';
    
28. U(:,2)=B(1,2)*C';
    
29. U(:,3)=B(1,3)*C';
    
30. %-------------------求解潮流-----------------
    : N$ D* _' i% r
31. for k=1:15
    ' o6 I' ^  @. U
32. %   I(:,k)=((DB(:,5).')*(diag(1./U(:,k))))';
    3 A9 P( e$ S  L4 O. o2 A
33. for i=n:-1:1
    
34.   %如果尾节点带恒功率负荷，需计算节点上负荷注入的电流
    . c$ @8 }. a' W8 b& d9 F4 y% n
35. if DB(i,7)==1
    # @4 a; l/ t* J  t& e
36. c=DB(i,5)/3;
    
37. d=c/U(i+1,(3*k-2));
    : @; V# v' Y0 p, }" y+ N; Q
38. IL(i,3*k-2)=conj(d);
    
39. d=c/U(i+1,3*k-1);
    
40. IL(i,3*k-1)=conj(d);
    
41. d=c/U(i+1,3*k);
    + \4 H' I. A$ m
42. IL(i,3*k)=conj(d);
    
43. else
    
44. IL(i,3*k-2)=0;
    # H- X/ s9 N8 q0 A" B: w3 A4 R
45. IL(i,3*k-1)=0;
    * H. r6 L7 q5 D- S; `0 A! b
46. IL(i,3*k)=0;
    : w; s+ F  H4 O$ y6 G
47. end
    
48. %找出所有与尾节点相连的支路,计算进支电流（末端电流），存IKj
      V# v" Q# U' ^  {% B
49. A=(find(DB(:,2)==(i+1)));
    - r# N7 @, W3 Z7 X
50. if isempty(A)~=1
    * h  x& ?# m) D" x1 c% ?
51. IKj(i,3*k-2)=IL(i,3*k-2)-sum(IKj(A,3*k-2));
    8 [( h2 Z9 W0 ?8 s
52. IKj(i,3*k-1)=IL(i,3*k-1)-sum(IKj(A,3*k-1));
    $ M0 Z4 c- R0 C" q5 B6 f
53. IKj(i,3*k)=IL(i,3*k)-sum(IKj(A,3*k));
    $ L7 r. h2 B: a; C+ Z# K% n3 `7 S
54. else
      y+ c2 F+ Y. V
55. IKj(i,3*k-2)=IL(i,3*k-2);
    8 m( N+ H' @- y6 H# ]2 [
56. IKj(i,3*k-1)=IL(i,3*k-1);
    4 |8 j: Z, }% J# ]% M+ o! b1 H6 n) J
57. IKj(i,3*k)=IL(i,3*k);
    
58. end
    2 N% ?6 R" W7 @7 \6 S" a$ X
59. %计算出支电流（始端电流），存IKi
    
60. a=DB(i,2);
    
61. b=DB(i,3);
    
62. Yi=1/DB(i,4)*eye(3);
    * M7 G; O0 G0 `* C2 i5 A
63. if DB(A,6)==1
    : F- s+ z5 a7 W& m( U7 C& B
64. F=0.5*Yi*[(U(a,3*k-2)+U(b,3*k-2)) (U(a,3*k-1)+U(b,3*k-1)) (U(a,3*k)+U(b,3*k))]'+[IKj(i,3*k-2) IKj(i,3*k-1) IKj(i,3*k)]';
    
65. IKi(i,3*k-2)=F(1,1);
    
66. IKi(i,3*k-1)=F(2,1);
    
67. IKi(i,3*k)=F(3,1);
    8 `' q( I2 d) e" X
68. else
    
69. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    2 }9 x% E& @# Q0 M
70. YI=YT*eye(3);
    2 t* A3 _' C6 j1 `
71. YII=YT*eye(3);
    
72. YIII=-YT*eye(3);
    
73. D=inv(YIII)*(-[IKj(i,3*k-2) IKj(i,3*k-1) IKj(i,3*k)]'-YII*([U(b,3*k-2) U(b,3*k-1) U(b,3*k)]'));
    
74. U(a,3*k+1)=D(1,1);
    
75. U(a,3*k+2)=D(2,1);
    4 d1 k9 x6 l# n9 a; i; v0 A
76. U(a,3*k+3)=D(3,1);
    % R% \: p& ]6 b* I7 r; Y
77. E=YI*[U(a,3*k+1) U(a,3*k+2) U(a,3*k+3)]'-YIII*[U(b,3*k-2) U(b,3*k-1) U(b,3*k)]';
    
78. IKi(i,3*k-2)=E(1,1);
    & Q# |# W% J- [. S
79. IKi(i,3*k-1)=E(2,1);
    
80. IKi(i,3*k)=E(3,1);
    " m$ Q- T( B# ~
81. end
    
82. end
    
83. %前推电压
    
84. for j=2:n
    # f$ p1 F$ A+ j: P
85. U(1,3*k-2)=1.02;
    
86. U(1,3*k-1)=1.02*(sin(2*pi/3)+j*cos(2*pi/3));
    
87. U(1,3*k)=1.02*(sin(4*pi/3)+j*cos(4*pi/3));
    
88. a=DB(j-1,2);
    " e& D/ A) ]' x
89. if DB(j-1,6)==1
    : A1 A9 b8 D5 e
90. Yi=1/DB(j-1,4)*eye(3);
    
91. G1=[IKi(a,3*k-2) IKi(a,3*k-1) IKi(a,3*k)]';
    & V1 c8 H2 c! I) a
92. G=[U(a,3*k-2) U(a,3*k-1) U(a,3*k)]'-DB(j-1,4)*eye(3)*(G1-0.5*Yi*[U(a,3*k-2) U(a,3*k-1) U(a,3*k)]');
    . q1 ]$ _! S9 _+ V6 @4 T1 `
93. U(j,3*k+1)=G(1,1);
    
94. U(j,3*k+2)=G(2,1);
    ' \, E# Q$ g5 q* v/ @# I3 p
95. U(j,3*k+3)=G(3,1);
    
96. else
    
97. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    
98. YI=YT*eye(3);
    ! _* {: M3 d. N
99. YII=YT*eye(3);
    6 V5 s1 N/ t, O$ R1 a
100. YIII=-YT*eye(3);
     ; Y/ L/ v" v0 X
101. H=inv(YIII)*([IKi(a,3*k-2) IKi(a,3*k-1) IKi(a,3*k)]'-YI*[U(a,3*k+1) U(a,3*k+2) U(a,3*k+3)]');
     
102. U(j,3*k+1)=H(1,1);
     
103. U(j,3*k+2)=H(2,1);
     ! Y; _1 d' A$ ^$ @, T
104. U(j,3*k+3)=H(3,1);
     ) C: S. d0 z! W# Z+ X
105. end
     
106. end
     
107. dU1=U(:,3*k+3)-U(:,3*k);
     % r( B! F% j3 @" R! [/ d2 F
108. dU2=U(:,3*k+2)-U(:,3*k-1);
     % G8 G4 q5 L( s7 @  l: h7 r  S/ ~( u
109. dU3=U(:,3*k+1)-U(:,3*k-2);
     8 q1 Z8 c& y" m$ v( @
110. if (max(abs(dU1))<0.001&&max(abs(dU2))<0.001&&max(abs(dU2))<0.001)
       }) U  h0 O, C7 y/ N, \
111. break;
     ( D7 C$ s2 Q; `6 T4 i
112. else            
     
113. k=k+1;
     
114. end
     5 E' h0 W4 y( ]2 X: \
115. end
     
116. U
     . \# V- ~' n! z+ x
117. IL
     * I0 Z9 n& k; N8 D
118. IKj
     
119. IKi
     
120. k
     / _' }2 Z4 |: l6 R# ^

复制代码

redplum

我运行了楼主的这段代码，发现最后结果矩阵中很多元素都是Nan或Inf，是不是潮流计算阻抗矩阵中多数情况是这样的，还是不应该很多元素是NaN或Inf

咻咻永远

回复 2# redplum
" z; X3 W, S7 J- S0 a4 l# a
. }1 ^# K: h9 i
3 N6 I3 V! x7 p) x7 `! F, Q0 Z    不应该是这样的，这是一个迭代计算，结果不收敛，所以还是哪个环节有问题，但是我还没有找到。

redplum

这段代码的原始数据有17行7列，K=1:15计算了15次循环。建议将数据弄少一些，这样才能更好判断到底是哪里不收敛

咻咻永远

回复 4# redplum
就是即使我弄成迭代次数为2，结果也不啊~~

魏丹

这个问题解决了吗 求答案~