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

咻咻永远

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

1. %程序名：qiantuihuitui_I_3.m
   % T9 I" V  B; r0 \& S3 e
2. %功能：支路电流前推回推法求解潮流
   
3. clc
   
4. clear all;
   
5. %--------------输入网络参数--------------
   , y& B, P' ^7 V* B7 |: C
6. %1-支路编号，2-首节点，3-尾节点，4-自阻抗，5-尾节点复功率，6-支路性质（1-馈线段支路，2-变压器支路）,7-尾节点是否带负荷
   
7. DB=[1 1 2 0.000167+j*0.000208 0.42+j*0.31 1 1
   $ ^( s0 Z* x7 R3 A& a
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
    , x5 _' f* a3 o
11.     5 2 6 0.000172+j*0.000215 0 1 0
    / n% Y  H+ L- D, r% k
12.     6 4 7 0.000156+j*0.000195 6.06 1 0
    % H4 r" i4 X3 h% o2 _6 C
13.     7 4 8 0.000162+j*0.000202 6.04 1 0
    6 {& G% w) K3 q
14.     8 4 9 0.000345+j*0.000430 0 1 0
    - g0 t2 }. f' A- C1 l7 K7 k
15.     9 4 10 0.000287+j*0.000358 0 1 0
    * q7 d# r7 m. L  K. O; _
16.     10  5 11 0.020563+j*0.321594 0 2 0
    . u" G# m3 ]. C, p- Z( x" k1 I
17.     11  6 12 0.020563+j*0.321594 0 2 0
    - U* |  X/ F: c/ _, w- D
18.     12  9 13 0.020563+j*0.321594 0 2 0
    8 _& Z7 a% I8 R. l
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
    " N& W% i' Y& M3 G$ Y  i
21.     15  12 16 0.000237+j*0.000408 5.76+j*0.09 1 1
    4 I  _& q; N2 N2 m- M. z* W
22.     16  13 17 0.000292+j*0.000502 5.86+j*0.11 1 1
    ' O5 |* k' T$ E0 I3 g/ V
23.     17  14 18 0.000274+j*0.000470 5.81+j*0.14 1 1];
    
24. [n,m]=size(DB);  
    
25. B=[1 sin(2*pi/3)+j*cos(2*pi/3) sin(4*pi/3)+j*cos(4*pi/3)];
    " n, j1 E% H  g; f
26. C=[1.02 1 1.02 1 1 1 1.02 1.02 1 1 1 1 1 1 1 1 1 1];
    5 ^' y9 f1 K$ t1 z6 J1 z6 a, ?. k# L
27. U(:,1)=B(1,1)*C';
    
28. U(:,2)=B(1,2)*C';
    
29. U(:,3)=B(1,3)*C';
    
30. %-------------------求解潮流-----------------
    
31. for k=1:15
    
32. %   I(:,k)=((DB(:,5).')*(diag(1./U(:,k))))';
    
33. for i=n:-1:1
    
34.   %如果尾节点带恒功率负荷，需计算节点上负荷注入的电流
    ' w+ M5 t" C3 Q: T* |6 I; @+ G8 i, }
35. if DB(i,7)==1
    4 Z5 ?) v8 }! G9 f+ Z
36. c=DB(i,5)/3;
    
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);
    - f7 ?* ^& X7 L
40. IL(i,3*k-1)=conj(d);
    # e4 _, }- k* P* g! \
41. d=c/U(i+1,3*k);
    
42. IL(i,3*k)=conj(d);
    1 _0 O+ w  w' t6 @1 H6 O
43. else
    - s# w) }# k9 N6 M
44. IL(i,3*k-2)=0;
    7 ^- a; o) W9 {$ B- C, b
45. IL(i,3*k-1)=0;
    
46. IL(i,3*k)=0;
      `4 m! D3 |8 w" x
47. end
    
48. %找出所有与尾节点相连的支路,计算进支电流（末端电流），存IKj
    , Q- |2 P8 h3 Z+ G; x9 r
49. A=(find(DB(:,2)==(i+1)));
    7 m( c9 D3 a0 y' X5 C
50. if isempty(A)~=1
    
51. IKj(i,3*k-2)=IL(i,3*k-2)-sum(IKj(A,3*k-2));
    0 M/ f9 y7 X# C. }+ T
52. IKj(i,3*k-1)=IL(i,3*k-1)-sum(IKj(A,3*k-1));
    ) G$ P* ?0 r" x4 a! w: s
53. IKj(i,3*k)=IL(i,3*k)-sum(IKj(A,3*k));
    
54. else
    
55. IKj(i,3*k-2)=IL(i,3*k-2);
    
56. IKj(i,3*k-1)=IL(i,3*k-1);
    % w* R6 M/ y/ Y* P6 \! M! Q& Z9 b
57. IKj(i,3*k)=IL(i,3*k);
    
58. end
    + ?5 s; r) a& B% |' y- g  A0 Z
59. %计算出支电流（始端电流），存IKi
    
60. a=DB(i,2);
    ' o- F% y& A: w* X
61. b=DB(i,3);
    & @0 m- e  y7 n
62. Yi=1/DB(i,4)*eye(3);
    , T1 o5 g9 g2 B  C
63. if DB(A,6)==1
    1 O; g2 g# q! M/ Z7 A4 R$ S( ^
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)]';
    2 _: E* U) @* e! _! p
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);
    
68. else
    
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);
    ! n4 `" l( ]9 x! J" p- |
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);
    & J3 _1 P9 l2 y9 _( e( \$ d, Y1 N, K' w
75. U(a,3*k+2)=D(2,1);
    
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)]';
    
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);
      f& n9 u1 Z8 ~( U* K! i' l) q' G' }
81. end
    # x- [$ k: B, J+ a/ d! l' s- a
82. end
    
83. %前推电压
    
84. for j=2:n
    
85. U(1,3*k-2)=1.02;
    ) i% d* G/ I4 r8 E( W) [: M+ a
86. U(1,3*k-1)=1.02*(sin(2*pi/3)+j*cos(2*pi/3));
    6 b) r, a8 q2 g( O5 w$ d' w& K
87. U(1,3*k)=1.02*(sin(4*pi/3)+j*cos(4*pi/3));
    
88. a=DB(j-1,2);
    4 A: _, m4 m) Q' Q7 A+ A9 q! A
89. if DB(j-1,6)==1
    # |1 d  d4 y* f. O- ]
90. Yi=1/DB(j-1,4)*eye(3);
    3 q3 X" [* ~" c5 H
91. G1=[IKi(a,3*k-2) IKi(a,3*k-1) IKi(a,3*k)]';
    
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);
    
95. U(j,3*k+3)=G(3,1);
    $ \# I$ i3 D+ J6 |
96. else
    * B  u9 Q; o! B* b2 i
97. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    
98. YI=YT*eye(3);
    0 J% U- `1 X2 p/ f8 |: H/ }
99. YII=YT*eye(3);
    7 S$ \2 S. V+ Z+ e) _
100. YIII=-YT*eye(3);
     
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);
     2 }7 L. |, q/ Z% a
103. U(j,3*k+2)=H(2,1);
     
104. U(j,3*k+3)=H(3,1);
       o9 Y5 V, w0 X4 U; h4 R
105. end
     
106. end
     
107. dU1=U(:,3*k+3)-U(:,3*k);
     
108. dU2=U(:,3*k+2)-U(:,3*k-1);
     ) g$ |! C8 U: r/ E* G
109. dU3=U(:,3*k+1)-U(:,3*k-2);
     " B0 ^( k8 x* r) v
110. if (max(abs(dU1))<0.001&&max(abs(dU2))<0.001&&max(abs(dU2))<0.001)
     
111. break;
     4 t, T3 X$ `' W2 L9 y% H( Q. N
112. else            
     0 A4 u& f% r& r0 `4 [
113. k=k+1;
     8 q: G7 A* D( Z, _, Y. L$ ?
114. end
     
115. end
     
116. U
     
117. IL
     # \; z# U5 k1 n/ c  X
118. IKj
     % [$ t" U, \- r* p% Y
119. IKi
     % j; ^+ P  _& U
120. k
     

复制代码

redplum

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

咻咻永远

回复 2# redplum

" k8 n2 J6 G4 r" a& Q# c
3 e* q& Y6 T2 w  a; A4 }    不应该是这样的，这是一个迭代计算，结果不收敛，所以还是哪个环节有问题，但是我还没有找到。

redplum

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

咻咻永远

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

魏丹

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