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

咻咻永远

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

1. %程序名：qiantuihuitui_I_3.m
   2 N4 P7 ]) h; b' l$ l
2. %功能：支路电流前推回推法求解潮流
   ) b# G* x0 w- q
3. clc
   
4. clear all;
   
5. %--------------输入网络参数--------------
   
6. %1-支路编号，2-首节点，3-尾节点，4-自阻抗，5-尾节点复功率，6-支路性质（1-馈线段支路，2-变压器支路）,7-尾节点是否带负荷
   9 L9 l# L2 o6 Q/ U9 w8 k
7. DB=[1 1 2 0.000167+j*0.000208 0.42+j*0.31 1 1
   5 n. I6 R* ]3 P  O/ h3 K( W" M
8.     2 2 3 0.000151+j*0.000188 6.15 1 0
   ; t0 l* C1 B# _; @0 s
9.     3 2 4 0.000066+j*0.000082 0.38+j*0.29 1 1
   1 D& i$ P7 G  K9 A% f
10.     4 2 5 0.000249+j*0.000310 0 1 0
    % F: Q: s9 v% }+ B/ ]
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
    ! W" [& ^; k( V
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
    
17.     11  6 12 0.020563+j*0.321594 0 2 0
    
18.     12  9 13 0.020563+j*0.321594 0 2 0
    * B! P, v) p: E1 a! v8 T. T4 I2 b
19.     13  10 14 0.020563+j*0.321594 0 2 0  
    0 a2 U5 }7 z6 W2 M2 C; S! j
20.     14  11 15 0.000237+j*0.000408 5.72+j*0.12 1 1
    & J( A7 B$ p. Q5 @2 t
21.     15  12 16 0.000237+j*0.000408 5.76+j*0.09 1 1
    & X0 |- S' r  ?6 w! f
22.     16  13 17 0.000292+j*0.000502 5.86+j*0.11 1 1
    ; g3 n! m+ }1 J0 H7 W$ O
23.     17  14 18 0.000274+j*0.000470 5.81+j*0.14 1 1];
    
24. [n,m]=size(DB);  
    % J7 {- b/ F6 x# ^& i$ p
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];
    
27. U(:,1)=B(1,1)*C';
    
28. U(:,2)=B(1,2)*C';
    
29. U(:,3)=B(1,3)*C';
    
30. %-------------------求解潮流-----------------
    . R& O! B" @& z1 r! u
31. for k=1:15
    & M) y/ g# e; i4 ?* P. X
32. %   I(:,k)=((DB(:,5).')*(diag(1./U(:,k))))';
    
33. for i=n:-1:1
    - Q" S+ l$ K6 q% P9 `. D$ e
34.   %如果尾节点带恒功率负荷，需计算节点上负荷注入的电流
    " U! k- u9 h3 }, R
35. if DB(i,7)==1
    ! h+ a4 g5 ]# F) @! l- U
36. c=DB(i,5)/3;
    
37. d=c/U(i+1,(3*k-2));
    1 t9 j8 D) a$ b3 C8 m
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);
    
42. IL(i,3*k)=conj(d);
    . l0 {9 n+ C  n/ \+ n
43. else
    
44. IL(i,3*k-2)=0;
    ! C( |# J4 B& b7 \4 K
45. IL(i,3*k-1)=0;
    6 j1 Y, R2 h* P0 y, a
46. IL(i,3*k)=0;
    - d# V) T% ^: u5 D% u
47. end
    
48. %找出所有与尾节点相连的支路,计算进支电流（末端电流），存IKj
    
49. A=(find(DB(:,2)==(i+1)));
    
50. if isempty(A)~=1
    : i+ `+ z! ]  W$ Q7 i& U# I2 i/ z/ r
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));
    5 |" H. E6 h1 y5 Y6 h
53. IKj(i,3*k)=IL(i,3*k)-sum(IKj(A,3*k));
    : ]/ r5 q3 i' p+ S1 ?  S) X% t  i
54. else
    9 h- O1 h6 e2 \
55. IKj(i,3*k-2)=IL(i,3*k-2);
      Z2 q, u1 m: [" I1 c' e
56. IKj(i,3*k-1)=IL(i,3*k-1);
    7 x# k2 @$ W, @" E6 X' p
57. IKj(i,3*k)=IL(i,3*k);
    
58. end
    
59. %计算出支电流（始端电流），存IKi
    
60. a=DB(i,2);
    8 w: W2 Q. ?; ]% O* h6 {  A, _- ~
61. b=DB(i,3);
    4 Y! T. ~+ j! a
62. Yi=1/DB(i,4)*eye(3);
    
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)]';
    , d/ [- K, U* L" }& P3 \( _
65. IKi(i,3*k-2)=F(1,1);
    ; T, r  d9 e- O1 R
66. IKi(i,3*k-1)=F(2,1);
    ! ~2 V. a9 o2 K7 a' G
67. IKi(i,3*k)=F(3,1);
    
68. else
    " x. \+ i( K; n+ k7 k+ V
69. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    ! F8 I0 S& d0 C9 K3 Y2 M
70. YI=YT*eye(3);
    3 D8 F1 Z5 i& S& s# H: [4 Z
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)]'));
    + |8 J" |+ a; \- q& L: _. r( c
74. U(a,3*k+1)=D(1,1);
    % O( k5 Z5 p: B4 \2 i) Z
75. U(a,3*k+2)=D(2,1);
    
76. U(a,3*k+3)=D(3,1);
    . h9 z2 }$ O9 S, m; c- k
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);
    5 l' h/ f2 v% k4 Y0 ?$ r7 u
80. IKi(i,3*k)=E(3,1);
    * k& C8 ]3 s  X
81. end
    9 {6 t* [6 w6 u- J# L1 I1 P8 r
82. end
    
83. %前推电压
    8 f1 Z$ G/ ]. O: H
84. for j=2:n
    , a0 M5 b" B) Z; S; H" X
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));
    4 ~: K+ g4 M4 q
88. a=DB(j-1,2);
    " B/ }( z5 P8 S" u! E
89. if DB(j-1,6)==1
    2 ^# h0 v/ e% l2 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)]';
    8 z/ S/ K0 z. g" g
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)]');
    3 r6 i7 P. Z' J# I+ o+ p+ K9 @
93. U(j,3*k+1)=G(1,1);
    ) n6 j$ Y- c! N4 u
94. U(j,3*k+2)=G(2,1);
    
95. U(j,3*k+3)=G(3,1);
    9 a; _$ a  w5 v4 {
96. else
    
97. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    
98. YI=YT*eye(3);
    
99. YII=YT*eye(3);
    3 M3 S1 x7 ~- h8 H  m% ?
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)]');
     ; \) ?- j0 C, |/ X
102. U(j,3*k+1)=H(1,1);
     2 U0 |8 v0 T$ K
103. U(j,3*k+2)=H(2,1);
     
104. U(j,3*k+3)=H(3,1);
     
105. end
     4 L) n; _! U6 T* A0 \; q( V
106. end
     
107. dU1=U(:,3*k+3)-U(:,3*k);
     $ [, t0 H! M; J" z! f7 O/ Q
108. dU2=U(:,3*k+2)-U(:,3*k-1);
     3 Z0 g7 a* U% A& c, S8 v
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)
     # k" Q; B0 \- r2 ]; p# `
111. break;
     ( k4 n; t, X( u9 |" m2 K5 |2 ]# v+ }) u
112. else            
     
113. k=k+1;
     ( _- r5 }% O6 k. j8 {3 p% W, `8 e: ~
114. end
     $ U+ _' U, x0 H; e- C) n
115. end
     
116. U
     ; D/ l0 }3 ~+ D$ m: y
117. IL
     1 a8 Y  }+ p4 Z0 W
118. IKj
     
119. IKi
     $ m  J1 P/ d) V
120. k
     # e$ G) g" n, V* T6 o

复制代码

redplum

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

咻咻永远

回复 2# redplum


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

redplum

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

咻咻永远

回复 4# redplum
4 E* q; n. a# _3 ?8 |. ^9 W就是即使我弄成迭代次数为2，结果也不啊~~

魏丹

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