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

咻咻永远

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

1. %程序名：qiantuihuitui_I_3.m
   % ^1 m: p6 Y# @8 ^% ]( p: |
2. %功能：支路电流前推回推法求解潮流
   " X& T. `. Y6 [, ~
3. clc
   % W- }8 T' S  e& G0 Z1 D
4. clear all;
   
5. %--------------输入网络参数--------------
   
6. %1-支路编号，2-首节点，3-尾节点，4-自阻抗，5-尾节点复功率，6-支路性质（1-馈线段支路，2-变压器支路）,7-尾节点是否带负荷
   . Y" c* \9 U: X* Z5 e+ H! C* a
7. DB=[1 1 2 0.000167+j*0.000208 0.42+j*0.31 1 1
   
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
   # p0 t7 u! I) I* H/ q! f
10.     4 2 5 0.000249+j*0.000310 0 1 0
    
11.     5 2 6 0.000172+j*0.000215 0 1 0
    . r' T/ ~! f; z& }
12.     6 4 7 0.000156+j*0.000195 6.06 1 0
    ' R* C8 P* [0 A3 R6 l
13.     7 4 8 0.000162+j*0.000202 6.04 1 0
    3 k9 B9 u) h  R* w& b" Q6 G& ?
14.     8 4 9 0.000345+j*0.000430 0 1 0
    & ^7 C+ n  ^! q# t! T
15.     9 4 10 0.000287+j*0.000358 0 1 0
    6 `  l, y! R+ O  @1 R. }$ m
16.     10  5 11 0.020563+j*0.321594 0 2 0
    
17.     11  6 12 0.020563+j*0.321594 0 2 0
    5 G' a9 j, m, T& |% B+ A+ c* x* j
18.     12  9 13 0.020563+j*0.321594 0 2 0
    1 |9 `7 L0 V2 i3 @" ]. ~5 c
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
    
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);  
    % ^2 u3 _1 M1 n- E
25. B=[1 sin(2*pi/3)+j*cos(2*pi/3) sin(4*pi/3)+j*cos(4*pi/3)];
      ]: p! w  y2 U4 }& b  P6 \1 f  g
26. C=[1.02 1 1.02 1 1 1 1.02 1.02 1 1 1 1 1 1 1 1 1 1];
    ; h: e  n6 Y* R. b% b9 O% A
27. U(:,1)=B(1,1)*C';
    8 y! v$ t+ z' q0 J+ ]
28. U(:,2)=B(1,2)*C';
    
29. U(:,3)=B(1,3)*C';
    
30. %-------------------求解潮流-----------------
    2 ^' j3 D. ^# k2 V, S: o
31. for k=1:15
    
32. %   I(:,k)=((DB(:,5).')*(diag(1./U(:,k))))';
    0 E6 v7 I# S7 \  ~. T+ B
33. for i=n:-1:1
    
34.   %如果尾节点带恒功率负荷，需计算节点上负荷注入的电流
    : p# V0 z' m  o: h
35. if DB(i,7)==1
    3 Y8 @! n8 q. b, l2 V8 e, B
36. c=DB(i,5)/3;
    4 O! L9 j9 V: X6 P# g
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);
    + \' y8 p$ F: l; U) ^. q$ j, J% b0 ]4 A
41. d=c/U(i+1,3*k);
    , S3 j% B8 @, U  v+ k( I
42. IL(i,3*k)=conj(d);
    
43. else
    
44. IL(i,3*k-2)=0;
    ) ^1 x" R2 A5 f8 F, }# m3 j' W
45. IL(i,3*k-1)=0;
    & U& |  l+ _' H* E% d1 o
46. IL(i,3*k)=0;
    7 ^. U( r$ M1 k+ I, A
47. end
    # C5 Z" d5 L( t' L0 o
48. %找出所有与尾节点相连的支路,计算进支电流（末端电流），存IKj
    
49. A=(find(DB(:,2)==(i+1)));
    ' D2 }& z! o$ v! [
50. if isempty(A)~=1
    4 j, E- e3 A4 k- }+ L! H
51. IKj(i,3*k-2)=IL(i,3*k-2)-sum(IKj(A,3*k-2));
    ( r& V6 K: X' ~) i& t
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));
    ' ~1 L$ W4 b1 N6 y
54. else
    * V4 ]/ X2 G+ d7 t+ j- {" e
55. IKj(i,3*k-2)=IL(i,3*k-2);
    2 x! p5 f3 ?4 |9 T  v8 T6 a
56. IKj(i,3*k-1)=IL(i,3*k-1);
    / d$ e. Q1 x3 I+ e3 B
57. IKj(i,3*k)=IL(i,3*k);
    
58. end
    
59. %计算出支电流（始端电流），存IKi
    
60. a=DB(i,2);
    
61. b=DB(i,3);
    4 r5 d1 b) M+ g) s' i
62. Yi=1/DB(i,4)*eye(3);
    $ u7 C$ y* I* h
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)]';
    9 F3 a: D" ]& U& v$ G' B
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 g, c" f# `$ l# i% {! _
68. else
    % s- b/ e% p3 x8 Q# S6 D) y
69. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    
70. YI=YT*eye(3);
    2 }9 P8 Y* Z: P  V' j* s5 a
71. YII=YT*eye(3);
    & {% u4 Y, h# v: v% y5 }3 L
72. YIII=-YT*eye(3);
    " m" r% F% z% W# |; I" L7 u8 ]8 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);
    - A( q4 q8 z/ @) O1 |
75. U(a,3*k+2)=D(2,1);
    
76. U(a,3*k+3)=D(3,1);
    $ I+ V  A, c6 G3 r% S0 J
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)]';
    * V+ }. |2 L1 I* c4 Y' H# J
78. IKi(i,3*k-2)=E(1,1);
    
79. IKi(i,3*k-1)=E(2,1);
    , q/ N3 x# D5 g
80. IKi(i,3*k)=E(3,1);
    
81. end
    
82. end
    
83. %前推电压
    % N# g) Z3 ]5 o) p# P. T; k
84. for j=2:n
    % w( Q. N6 ]5 }  w9 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));
    * \' c8 a( C4 y. }
87. U(1,3*k)=1.02*(sin(4*pi/3)+j*cos(4*pi/3));
    3 w0 k) @0 ]- p7 X0 v" `9 o3 J9 |
88. a=DB(j-1,2);
    / ?+ [  v9 G) Q: g8 Y
89. if DB(j-1,6)==1
    4 W1 K& |2 [$ s/ x
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)]');
    9 k9 S7 W' R0 h$ `3 o/ x
93. U(j,3*k+1)=G(1,1);
    ' b" X& @- B5 E6 M, x: l
94. U(j,3*k+2)=G(2,1);
    
95. U(j,3*k+3)=G(3,1);
      j* O* a9 [0 y) p$ Z! M
96. else
    ; F% f0 b  y( T! I( a
97. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    
98. YI=YT*eye(3);
    # c8 s' H$ e5 O, [9 p3 n
99. YII=YT*eye(3);
    ) V/ B& x  t! i3 J+ M7 a% ?
100. YIII=-YT*eye(3);
     : ~- \- C$ \) V
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);
     + S1 ^& f  d4 i0 K; M8 M+ Q6 ^
103. U(j,3*k+2)=H(2,1);
     
104. U(j,3*k+3)=H(3,1);
     ; q5 |! l0 m. {3 J" {
105. end
     5 s3 n. p* T& H0 ^+ z/ s3 u+ L& W
106. end
     
107. dU1=U(:,3*k+3)-U(:,3*k);
     " q" }$ Y( J. k5 J1 m9 c
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)
     9 ^/ z. u# q2 R1 Q' W2 A
111. break;
     
112. else            
     2 y' B( L( k, ?" ]/ o
113. k=k+1;
     
114. end
     1 L) H' P  x+ ]& E* y
115. end
     
116. U
     ; m4 q1 S# T: a0 D! [
117. IL
     4 `5 k% @# J, H
118. IKj
     
119. IKi
     
120. k
     2 g; q  o; v- z' A$ o! P* ^/ R2 L

复制代码

redplum

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

咻咻永远

回复 2# redplum


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

redplum

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

咻咻永远

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

魏丹

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