电力系统及其自动化--潮流计算--咻咻永远

咻咻永远

各位，我想请教一个问题，我现在做的船舶交流电力系统，想要仿真一个18节点破冰船舶电力系统，主要参考的文献为《船舶电力系统拓扑分析和潮流计算研究-周容华》，需要的是3相潮流计算，但是现在编写的程序不收敛，想请大家帮我看看问题主要出在哪里？源代码在下边。

1. %程序名：qiantuihuitui_I_3.m
   # I$ M9 {+ a5 ~. r+ C2 O
2. %功能：支路电流前推回推法求解潮流
   / \% D( [8 M! a2 h& L, i
3. clc
   
4. clear all;
   
5. %--------------输入网络参数--------------
   % [$ P% m* i7 e6 w3 Z6 L. [! {% |/ E
6. %1-支路编号，2-首节点，3-尾节点，4-自阻抗，5-尾节点复功率，6-支路性质（1-馈线段支路，2-变压器支路）,7-尾节点是否带负荷
   2 c4 q! q, o/ v# d- o
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
   - N% A/ m: l9 c4 |3 V
9.     3 2 4 0.000066+j*0.000082 0.38+j*0.29 1 1
   * b+ Q9 I* m. r! n7 Y2 G# S& e
10.     4 2 5 0.000249+j*0.000310 0 1 0
    * z, K0 {' \* G$ O: Y* x+ X. \
11.     5 2 6 0.000172+j*0.000215 0 1 0
    + `- q: V: n5 M( L2 b6 H
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
    0 T2 Y! ^. m# z
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
    
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
    5 k( p) c0 ^; y1 B; n- ?" G
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);  
    1 V' `0 u) w; P" {& M! h& q
25. B=[1 sin(2*pi/3)+j*cos(2*pi/3) sin(4*pi/3)+j*cos(4*pi/3)];
    / m1 X  g$ u5 Y% Z9 O
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';
    + I$ W# _4 |  n! d) s
28. U(:,2)=B(1,2)*C';
    2 o! Y; a, Q& U9 a
29. U(:,3)=B(1,3)*C';
    
30. %-------------------求解潮流-----------------
    # M- l/ g1 t/ w, X+ v
31. for k=1:15
    8 {+ K) }, k4 S' w& d0 V4 z* f2 y
32. %   I(:,k)=((DB(:,5).')*(diag(1./U(:,k))))';
    ! G3 X! e7 j8 \+ m0 j* H
33. for i=n:-1:1
    ) V7 p# t9 r* Y& H; C' J0 u( N
34.   %如果尾节点带恒功率负荷，需计算节点上负荷注入的电流
    6 }# g" b* J) X2 p% K
35. if DB(i,7)==1
    ! X' C* b+ x' O
36. c=DB(i,5)/3;
    4 `# ?5 q  c$ L' d+ W7 ?
37. d=c/U(i+1,(3*k-2));
    2 ~4 q6 ]# c$ a0 B4 W  ]
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);
      y/ z( `( G* w# j, p, ]) X) Z: A  d
41. d=c/U(i+1,3*k);
    ( b0 N) _3 o, Q5 _5 A
42. IL(i,3*k)=conj(d);
    
43. else
    
44. IL(i,3*k-2)=0;
    9 o/ Z& v  r# T9 W6 e( V
45. IL(i,3*k-1)=0;
    6 T, H* I& L* J; m4 c7 g; b
46. IL(i,3*k)=0;
    
47. end
    
48. %找出所有与尾节点相连的支路,计算进支电流（末端电流），存IKj
    ! s) L- t0 N  ]7 n9 \8 c
49. A=(find(DB(:,2)==(i+1)));
    
50. if isempty(A)~=1
    8 H& H/ ^) U# j  J2 T2 G
51. IKj(i,3*k-2)=IL(i,3*k-2)-sum(IKj(A,3*k-2));
    " {" M7 J" w: y7 r) A
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));
    
54. else
    
55. IKj(i,3*k-2)=IL(i,3*k-2);
    
56. IKj(i,3*k-1)=IL(i,3*k-1);
    - W8 x, w$ l5 R9 C
57. IKj(i,3*k)=IL(i,3*k);
    ) t" G" g$ M  q& S8 j
58. end
    ; E1 c. j" k/ a1 G
59. %计算出支电流（始端电流），存IKi
    
60. a=DB(i,2);
    6 _- G! G- ~0 \
61. b=DB(i,3);
    , M$ s) G" M7 R4 Q0 r5 n3 O
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)]';
    $ r$ |8 A( E* p2 ?
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);
    0 Z3 ]7 E" [5 z4 N$ |' L* |+ H, o
68. else
    + A9 a& n1 A. L" J5 I. ~
69. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    / l) \3 g5 |9 Z0 m
70. YI=YT*eye(3);
    $ m' }/ [  \8 X, u3 U
71. YII=YT*eye(3);
    
72. YIII=-YT*eye(3);
    ( J, N: ^8 _/ t5 W5 m/ w9 k
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);
    
76. U(a,3*k+3)=D(3,1);
    % d" [, q3 R; ]" O* I1 F/ E
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);
      G- F5 H; E+ \. }
79. IKi(i,3*k-1)=E(2,1);
    
80. IKi(i,3*k)=E(3,1);
    
81. end
    
82. end
    4 c; {3 A  x- l! h1 r. R; Z1 X2 p
83. %前推电压
    
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));
    ' ^. u# f2 o+ g% }
87. U(1,3*k)=1.02*(sin(4*pi/3)+j*cos(4*pi/3));
    8 j. v+ v4 {3 s
88. a=DB(j-1,2);
    
89. if DB(j-1,6)==1
    ) W! x8 K! s9 y9 _" V2 f/ B- B
90. Yi=1/DB(j-1,4)*eye(3);
    6 @" D; f7 p/ e/ [: ?
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);
    
96. else
    6 ?+ [" B3 ?0 j& L' v$ M, q! }9 W
97. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    
98. YI=YT*eye(3);
    
99. YII=YT*eye(3);
    ) Q  A  k* Z) m
100. YIII=-YT*eye(3);
     ' h$ A+ O6 a; J2 @
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 H" }' f, S, m2 b: R7 F2 U
103. U(j,3*k+2)=H(2,1);
     / @+ Q. e) Y8 F  L: R; O9 H
104. U(j,3*k+3)=H(3,1);
     : n0 ^6 ?7 j7 v, I7 ~9 `0 ?: }
105. end
     * s: o$ b& B9 `) R, u/ P! l
106. end
     - u! s7 ]0 Y3 `4 e
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)
     : N- Q, M. {# ?% @
111. break;
     7 D* Z( _; t! r) U
112. else            
     
113. k=k+1;
     8 G- D0 m+ |' {$ }
114. end
     - b& m7 I. i5 w, C% v
115. end
     
116. U
     & v  v% s/ `$ }. l* }  g
117. IL
     " \! ]! X: |" \! u* V: ]  f5 B$ c
118. IKj
     
119. IKi
     
120. k
     ' c6 l3 r1 P9 z" m

复制代码

sun710504

学习学习