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

咻咻永远

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

1. %程序名：qiantuihuitui_I_3.m
   % d  A' o0 E  A" ?! o
2. %功能：支路电流前推回推法求解潮流
   
3. clc
   / i1 Q3 x) m. _7 e5 p
4. clear all;
   3 ^# k" L  L* ]" ^1 E' }, K
5. %--------------输入网络参数--------------
   
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
   : n7 u+ y; ?; m, c9 F: ]9 a& l& T
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
     w- z8 a! N5 F1 M
10.     4 2 5 0.000249+j*0.000310 0 1 0
    ( c3 z" [5 ~4 H0 O
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
    
15.     9 4 10 0.000287+j*0.000358 0 1 0
    ) N9 C2 z7 _3 ~% l# ?) z
16.     10  5 11 0.020563+j*0.321594 0 2 0
    6 \# t8 i1 M3 |; I& `# G# e& s
17.     11  6 12 0.020563+j*0.321594 0 2 0
    0 v% e! A( d, Z
18.     12  9 13 0.020563+j*0.321594 0 2 0
    
19.     13  10 14 0.020563+j*0.321594 0 2 0  
    9 J' ^4 R* O5 q+ V2 A- Y6 `4 e
20.     14  11 15 0.000237+j*0.000408 5.72+j*0.12 1 1
    1 \7 p9 a4 H: E# z7 D# u# h# T
21.     15  12 16 0.000237+j*0.000408 5.76+j*0.09 1 1
    ! o5 o) A2 ~6 W, i
22.     16  13 17 0.000292+j*0.000502 5.86+j*0.11 1 1
    4 s1 V$ @* o- s2 n+ _6 |
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)];
    
26. C=[1.02 1 1.02 1 1 1 1.02 1.02 1 1 1 1 1 1 1 1 1 1];
    ; p5 C* _* y9 W6 y) j
27. U(:,1)=B(1,1)*C';
    7 M, X& W, c. O7 s" u  y! O9 x
28. U(:,2)=B(1,2)*C';
    + X0 j9 l5 E) z
29. U(:,3)=B(1,3)*C';
    # U. M. W1 i$ F$ \
30. %-------------------求解潮流-----------------
    : P$ }$ `( C1 b+ G8 Y. w" l8 @
31. for k=1:15
    
32. %   I(:,k)=((DB(:,5).')*(diag(1./U(:,k))))';
    5 _# i% l4 o1 V! f  V/ d7 J
33. for i=n:-1:1
    
34.   %如果尾节点带恒功率负荷，需计算节点上负荷注入的电流
    , j; b) i+ p5 g/ q2 C  E( o8 E; |- m
35. if DB(i,7)==1
    : v9 F2 _: I  c/ R) x) _
36. c=DB(i,5)/3;
    
37. d=c/U(i+1,(3*k-2));
    
38. IL(i,3*k-2)=conj(d);
    " C& X; T" H/ R' N4 [
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);
    
42. IL(i,3*k)=conj(d);
    
43. else
    # y2 a' l9 |0 L( p: s
44. IL(i,3*k-2)=0;
    - h1 D% w7 ~! C/ L
45. IL(i,3*k-1)=0;
    ) a2 d7 L/ F( C3 d
46. IL(i,3*k)=0;
    # t: [2 U$ v6 i, P+ B& z
47. end
    
48. %找出所有与尾节点相连的支路,计算进支电流（末端电流），存IKj
    
49. A=(find(DB(:,2)==(i+1)));
    
50. if isempty(A)~=1
    - I% k" e9 z+ g: ^1 w  a
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));
    
53. IKj(i,3*k)=IL(i,3*k)-sum(IKj(A,3*k));
    % P! [& w3 ]7 y. ]& H; `
54. else
    , V) w8 t$ K3 b$ a" `3 P
55. IKj(i,3*k-2)=IL(i,3*k-2);
    8 ~5 G( K$ E' d5 y
56. IKj(i,3*k-1)=IL(i,3*k-1);
    * d8 b1 D9 M7 F8 N% h  A6 w) Y: P, L
57. IKj(i,3*k)=IL(i,3*k);
    7 ?% U- E3 d; }9 [, D
58. end
    . w1 H9 D! b* t
59. %计算出支电流（始端电流），存IKi
    
60. a=DB(i,2);
    
61. b=DB(i,3);
    
62. Yi=1/DB(i,4)*eye(3);
    & l  d1 D5 o9 o, g  b+ n
63. if DB(A,6)==1
    $ L; P; @+ F; w. U5 _
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);
    1 W0 ]+ E: f- {3 T( d) U1 \0 s$ o4 S* p
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)));
    1 V. ?" e! t( S7 z
70. YI=YT*eye(3);
    & }# c6 k7 a+ u2 A1 {% o/ P+ k/ t
71. YII=YT*eye(3);
    
72. YIII=-YT*eye(3);
    9 w, ~9 y' P* k; c4 m  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)]'));
    1 t: c& A5 Y; R0 K' P
74. U(a,3*k+1)=D(1,1);
    : z* J7 W: U) j, T6 y/ S
75. U(a,3*k+2)=D(2,1);
    
76. U(a,3*k+3)=D(3,1);
    ! J' U8 s4 F! J" y3 s( 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)]';
    ! _: r+ [9 o% _7 r- g: v
78. IKi(i,3*k-2)=E(1,1);
    / t8 w6 z" D8 P: F; n7 z
79. IKi(i,3*k-1)=E(2,1);
    % d6 y" x, |% X* e5 `
80. IKi(i,3*k)=E(3,1);
    2 W- V) s8 l; G4 x, r
81. end
    ; u9 p  c% ^/ R+ c0 E" X
82. end
    
83. %前推电压
    : V# I- y6 s! Z; i# l
84. for j=2:n
    * a/ j& }1 {) t2 R4 i1 }
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));
    : m- H" W, }) v( r0 C
88. a=DB(j-1,2);
    5 E0 N+ q1 ~: s
89. if DB(j-1,6)==1
    - i. N' U4 S' W5 h0 f  l
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)]';
    
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)]');
    - r- {8 H4 k1 o7 v. {$ K9 l  {
93. U(j,3*k+1)=G(1,1);
    
94. U(j,3*k+2)=G(2,1);
    + Z* V% X5 k- L+ X+ f
95. U(j,3*k+3)=G(3,1);
      m3 q$ O- J3 O. y" S
96. else
    
97. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    
98. YI=YT*eye(3);
    6 K9 ^8 Q" X% A" q0 v2 ?
99. YII=YT*eye(3);
    
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)]');
     1 n) S1 y3 G2 x  S
102. U(j,3*k+1)=H(1,1);
     # \* K% v- W& V5 z& Q
103. U(j,3*k+2)=H(2,1);
     ) @2 v- l. T9 H. M4 J1 a6 a
104. U(j,3*k+3)=H(3,1);
     
105. end
     
106. end
     3 c' @5 n' K' U5 t$ H% @& f% s
107. dU1=U(:,3*k+3)-U(:,3*k);
     
108. dU2=U(:,3*k+2)-U(:,3*k-1);
     
109. dU3=U(:,3*k+1)-U(:,3*k-2);
     
110. if (max(abs(dU1))<0.001&&max(abs(dU2))<0.001&&max(abs(dU2))<0.001)
     
111. break;
     ) ]/ y1 L6 f2 n) W0 [1 _/ A
112. else            
     
113. k=k+1;
     4 o& \4 {5 B" Q1 f$ P
114. end
     + G3 U  z  @5 s6 s. U9 S- {# p( B
115. end
     
116. U
     
117. IL
     
118. IKj
     
119. IKi
     1 i2 A: f8 ]# w9 C
120. k
     

复制代码

redplum

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

咻咻永远

回复 2# redplum


    不应该是这样的，这是一个迭代计算，结果不收敛，所以还是哪个环节有问题，但是我还没有找到。

redplum

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

咻咻永远

回复 4# redplum
( s1 w# |; e6 I% H: S1 r就是即使我弄成迭代次数为2，结果也不啊~~

魏丹

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