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

咻咻永远

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

1. %程序名：qiantuihuitui_I_3.m
   
2. %功能：支路电流前推回推法求解潮流
   
3. clc
   
4. clear all;
   " u3 y1 q7 `2 V0 Q
5. %--------------输入网络参数--------------
   
6. %1-支路编号，2-首节点，3-尾节点，4-自阻抗，5-尾节点复功率，6-支路性质（1-馈线段支路，2-变压器支路）,7-尾节点是否带负荷
   ) f& D7 f. f! ~" q5 v$ U/ P
7. DB=[1 1 2 0.000167+j*0.000208 0.42+j*0.31 1 1
   0 M, ~8 E. q( c$ W
8.     2 2 3 0.000151+j*0.000188 6.15 1 0
   + r2 k, }6 K* ^  _( [; p- A5 I
9.     3 2 4 0.000066+j*0.000082 0.38+j*0.29 1 1
   & ~1 k5 u. y3 f& E+ s9 d1 K5 Q
10.     4 2 5 0.000249+j*0.000310 0 1 0
    : L# E- o% e7 N" a3 E; [
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
    
14.     8 4 9 0.000345+j*0.000430 0 1 0
    ) o; \- j+ |& {" I& V/ l0 w
15.     9 4 10 0.000287+j*0.000358 0 1 0
    7 v: g* u3 a3 Z* ~3 C* J
16.     10  5 11 0.020563+j*0.321594 0 2 0
    
17.     11  6 12 0.020563+j*0.321594 0 2 0
    , g. p, v2 n6 ]/ m, s8 \8 `& K  P
18.     12  9 13 0.020563+j*0.321594 0 2 0
    , i* i4 u5 r% l! B# G" C/ v) Q
19.     13  10 14 0.020563+j*0.321594 0 2 0  
    + w( C/ k! A& f- a& Y
20.     14  11 15 0.000237+j*0.000408 5.72+j*0.12 1 1
    
21.     15  12 16 0.000237+j*0.000408 5.76+j*0.09 1 1
    
22.     16  13 17 0.000292+j*0.000502 5.86+j*0.11 1 1
    
23.     17  14 18 0.000274+j*0.000470 5.81+j*0.14 1 1];
    
24. [n,m]=size(DB);  
    / k$ l2 P6 x$ |1 u- W) o+ |
25. B=[1 sin(2*pi/3)+j*cos(2*pi/3) sin(4*pi/3)+j*cos(4*pi/3)];
    
26. C=[1.02 1 1.02 1 1 1 1.02 1.02 1 1 1 1 1 1 1 1 1 1];
    0 Z% j8 M! s( a! z
27. U(:,1)=B(1,1)*C';
    , u" M2 h3 I& i& g3 ]8 C; r+ Y
28. U(:,2)=B(1,2)*C';
    
29. U(:,3)=B(1,3)*C';
    ; T) s3 }2 y: K' ]# m( K  M( S2 |
30. %-------------------求解潮流-----------------
    & J5 s; s9 [' q6 p- a
31. for k=1:15
    ( I! d5 z& D5 o- j- z
32. %   I(:,k)=((DB(:,5).')*(diag(1./U(:,k))))';
    
33. for i=n:-1:1
    ( V+ D0 r/ g; M  w- s) k
34.   %如果尾节点带恒功率负荷，需计算节点上负荷注入的电流
    : m! [6 ^2 X0 D1 F3 D2 d/ p' o
35. if DB(i,7)==1
    
36. c=DB(i,5)/3;
    . M0 Y/ E) E( f1 U
37. d=c/U(i+1,(3*k-2));
    
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);
    2 R( N9 x5 S9 \: R
42. IL(i,3*k)=conj(d);
    0 l5 @! P6 h/ G+ z9 G6 ]5 N4 [( g
43. else
    ( E& x) I2 \( C1 n, S  S. @, i
44. IL(i,3*k-2)=0;
    ! d; E: n2 S, z( X2 \
45. IL(i,3*k-1)=0;
    
46. IL(i,3*k)=0;
      p$ w4 D! ~3 \! w
47. end
    
48. %找出所有与尾节点相连的支路,计算进支电流（末端电流），存IKj
    , G6 g& S# Z! ?9 n2 @+ ?% b
49. A=(find(DB(:,2)==(i+1)));
    * b% I4 t6 q* T- }4 K0 ?$ [
50. if isempty(A)~=1
    
51. IKj(i,3*k-2)=IL(i,3*k-2)-sum(IKj(A,3*k-2));
    
52. IKj(i,3*k-1)=IL(i,3*k-1)-sum(IKj(A,3*k-1));
    - {# s6 ~5 k7 [3 ~1 u1 m* I
53. IKj(i,3*k)=IL(i,3*k)-sum(IKj(A,3*k));
    
54. else
    ' W- H; e% h5 S; u/ w' i7 V' A! s1 `
55. IKj(i,3*k-2)=IL(i,3*k-2);
    - L4 n$ S5 R# e6 P9 F# G( f
56. IKj(i,3*k-1)=IL(i,3*k-1);
    1 r3 A. P. z! ]3 A
57. IKj(i,3*k)=IL(i,3*k);
    
58. end
    ' K: L+ o$ l. Z; Z6 Z0 J
59. %计算出支电流（始端电流），存IKi
    2 S% T+ a8 o9 W1 P) N. D$ z
60. a=DB(i,2);
    - `3 {- e! @  l  H7 b
61. b=DB(i,3);
    8 p) n2 e5 H+ {/ E0 n/ ^" R( |. F9 j
62. Yi=1/DB(i,4)*eye(3);
    ' _6 p3 G8 M* X8 T# v
63. if DB(A,6)==1
    
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)]';
      N. H" J) j0 g, }  V& b1 a
65. IKi(i,3*k-2)=F(1,1);
    1 _& b  J' e3 A0 l5 O
66. IKi(i,3*k-1)=F(2,1);
    ; C+ {3 g1 Q# ~6 [
67. IKi(i,3*k)=F(3,1);
    
68. else
    + t0 S# X6 t' F* s! @/ c8 X
69. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    
70. YI=YT*eye(3);
    
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);
    % N' O( ?2 F# f  a& f4 Q7 W
76. U(a,3*k+3)=D(3,1);
    
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)]';
    " l, X- P& i6 a: B* r- K" R8 O* V
78. IKi(i,3*k-2)=E(1,1);
    
79. IKi(i,3*k-1)=E(2,1);
    
80. IKi(i,3*k)=E(3,1);
    
81. end
    : ^* n! L% d& U" X1 z+ j% N0 N
82. end
    
83. %前推电压
    + ]% d+ z& L! b$ U* i# i4 D
84. for j=2:n
    
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));
    2 q- r$ @6 r% S
87. U(1,3*k)=1.02*(sin(4*pi/3)+j*cos(4*pi/3));
    % Z$ N" z# u1 y
88. a=DB(j-1,2);
    2 F' `( {0 u7 m6 b% w" ^, u
89. if DB(j-1,6)==1
    
90. Yi=1/DB(j-1,4)*eye(3);
    7 y8 C  P# U2 `( K( @: @
91. G1=[IKi(a,3*k-2) IKi(a,3*k-1) IKi(a,3*k)]';
    0 `$ J: h# a) c7 M
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)]');
    
93. U(j,3*k+1)=G(1,1);
    
94. U(j,3*k+2)=G(2,1);
    9 R. u. _, k: P* \6 P. s
95. U(j,3*k+3)=G(3,1);
    
96. else
    
97. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    ; |" L* m2 @) m+ W: M; p: r6 r
98. YI=YT*eye(3);
    
99. YII=YT*eye(3);
    , F& h$ }3 E( e  n1 P2 u5 n0 t% D1 o- w
100. YIII=-YT*eye(3);
     / O! O2 ]! y4 T
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)]');
     . U$ Z1 }+ D4 F) W
102. U(j,3*k+1)=H(1,1);
     
103. U(j,3*k+2)=H(2,1);
     
104. U(j,3*k+3)=H(3,1);
     . g! M3 w$ o% t
105. end
     
106. end
     
107. dU1=U(:,3*k+3)-U(:,3*k);
     + m5 a. t$ w$ G. E
108. dU2=U(:,3*k+2)-U(:,3*k-1);
     ) r1 j( P% |' H4 G# T7 t
109. dU3=U(:,3*k+1)-U(:,3*k-2);
       V# |5 s: V; [7 W' z4 p
110. if (max(abs(dU1))<0.001&&max(abs(dU2))<0.001&&max(abs(dU2))<0.001)
     
111. break;
     7 {# r' Y) C. i; h, _, d1 X1 {
112. else            
     
113. k=k+1;
     
114. end
     
115. end
     
116. U
     
117. IL
     8 k& p: C. |7 D6 _
118. IKj
     6 i2 ^( b, ~- t/ C% E$ j
119. IKi
     
120. k
     8 _) Y) H7 b* L" k

复制代码

redplum

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

咻咻永远

回复 2# redplum
" P6 d9 ?3 @! r# }3 B

, e9 H+ V4 ^3 v( b/ C" n! t    不应该是这样的，这是一个迭代计算，结果不收敛，所以还是哪个环节有问题，但是我还没有找到。

redplum

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

咻咻永远

回复 4# redplum
9 ^# |& u2 h& [7 [就是即使我弄成迭代次数为2，结果也不啊~~

魏丹

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