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

咻咻永远

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

1. %程序名：qiantuihuitui_I_3.m
   , y  e/ W4 C" n& t+ ~* u
2. %功能：支路电流前推回推法求解潮流
   ( m) i: w: m4 Z. @% r6 {
3. clc
   , R: Z5 a. @% M2 m* h
4. clear all;
   1 Z) Q, M, d; p" V; n1 w( ]# T
5. %--------------输入网络参数--------------
   
6. %1-支路编号，2-首节点，3-尾节点，4-自阻抗，5-尾节点复功率，6-支路性质（1-馈线段支路，2-变压器支路）,7-尾节点是否带负荷
   ! b9 o' O) Y) W; o2 N0 v* T3 z
7. DB=[1 1 2 0.000167+j*0.000208 0.42+j*0.31 1 1
   & F9 C9 O  t+ J' F+ T. C
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
   7 t: J8 x- k6 g- Z* q8 Z/ s2 C
10.     4 2 5 0.000249+j*0.000310 0 1 0
    
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
    7 S! c5 g8 ~$ x9 \
14.     8 4 9 0.000345+j*0.000430 0 1 0
      n( W' }1 r2 w
15.     9 4 10 0.000287+j*0.000358 0 1 0
    ; s+ v( N0 W' R% O7 C8 B6 v
16.     10  5 11 0.020563+j*0.321594 0 2 0
    : H+ B! C2 q( G! }* S. T. }
17.     11  6 12 0.020563+j*0.321594 0 2 0
    ) B  J$ v: n+ [8 Q" ?
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
    
22.     16  13 17 0.000292+j*0.000502 5.86+j*0.11 1 1
    1 Z. t- V& \/ S8 `% h( ^
23.     17  14 18 0.000274+j*0.000470 5.81+j*0.14 1 1];
    2 x+ A4 x/ ]* p8 q) Z
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)];
    8 u5 ~* Z: o% M. U  u$ T
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';
    ; J1 b- D5 b& T& y: P  Z  @
28. U(:,2)=B(1,2)*C';
    ; G5 N' v5 }( G' g3 S5 i
29. U(:,3)=B(1,3)*C';
    * v' g; e( N% t" |5 N0 j
30. %-------------------求解潮流-----------------
    
31. for k=1:15
    2 Z! T5 o3 d; f& `) p. t
32. %   I(:,k)=((DB(:,5).')*(diag(1./U(:,k))))';
    0 N. H0 a' Z4 B) x& N% R
33. for i=n:-1:1
    3 I; n, I0 t& v" J0 I% l" J
34.   %如果尾节点带恒功率负荷，需计算节点上负荷注入的电流
    5 i8 z; S  t! G: `& h0 c
35. if DB(i,7)==1
    
36. c=DB(i,5)/3;
    
37. d=c/U(i+1,(3*k-2));
    / x; M5 b) t: i6 v- q) P! A
38. IL(i,3*k-2)=conj(d);
    . v9 I* e0 e2 s8 S* z
39. d=c/U(i+1,3*k-1);
    " u& r$ }9 z' i: e  F
40. IL(i,3*k-1)=conj(d);
    
41. d=c/U(i+1,3*k);
    
42. IL(i,3*k)=conj(d);
    " {' Y  f7 F" a' o3 O( Z/ I: u& D
43. else
    3 y2 u' `5 v( ]4 Y6 V, i9 Q4 s1 ~
44. IL(i,3*k-2)=0;
    
45. IL(i,3*k-1)=0;
    
46. IL(i,3*k)=0;
    
47. end
    
48. %找出所有与尾节点相连的支路,计算进支电流（末端电流），存IKj
    
49. A=(find(DB(:,2)==(i+1)));
    
50. if isempty(A)~=1
    5 B0 X  J1 n/ M1 t
51. IKj(i,3*k-2)=IL(i,3*k-2)-sum(IKj(A,3*k-2));
    / g7 y% U* c6 Y4 Q
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
    , R4 D& c1 {$ z5 K2 p( l1 _
55. IKj(i,3*k-2)=IL(i,3*k-2);
    3 i& p2 V- l3 H- Q( E
56. IKj(i,3*k-1)=IL(i,3*k-1);
    
57. IKj(i,3*k)=IL(i,3*k);
    9 ]5 W8 C; J6 M8 A- r8 Q
58. end
    ! Z& b' r, S5 _/ f; F- H+ r
59. %计算出支电流（始端电流），存IKi
    
60. a=DB(i,2);
    
61. b=DB(i,3);
    3 z! _& F" l1 A: N+ ?# u
62. Yi=1/DB(i,4)*eye(3);
    4 X8 _0 Z8 k7 [" j
63. if DB(A,6)==1
    ! }! r% Z( O2 W1 V3 f2 P
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)]';
    ' m" h5 L' [% L' T7 w
65. IKi(i,3*k-2)=F(1,1);
    ; }6 R% z" _6 q  F! Q
66. IKi(i,3*k-1)=F(2,1);
    
67. IKi(i,3*k)=F(3,1);
    
68. else
    " T' p8 }8 e3 h$ D! W
69. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    1 l: @+ s/ B1 p. d; j0 N% }# ~3 s
70. YI=YT*eye(3);
    
71. YII=YT*eye(3);
    
72. YIII=-YT*eye(3);
    9 |0 V/ r' t6 f0 h
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)]'));
    ; }! G% a: p* p% ~& a0 D
74. U(a,3*k+1)=D(1,1);
    - o- ~; z1 D5 K/ r
75. U(a,3*k+2)=D(2,1);
    : Q2 q3 c& |9 {- |9 I0 f( U
76. U(a,3*k+3)=D(3,1);
    $ Q# U& w! V) l" U$ V/ Q
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);
    6 e/ r+ B& S) ]- x) F1 o
79. IKi(i,3*k-1)=E(2,1);
    ' ^* y) b: n$ G0 C
80. IKi(i,3*k)=E(3,1);
    
81. end
    , [5 w8 `8 ?; g: T- }6 Q* u; X
82. end
    
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));
    
87. U(1,3*k)=1.02*(sin(4*pi/3)+j*cos(4*pi/3));
    & a% q2 ]* \9 L8 M/ N( y
88. a=DB(j-1,2);
    5 c- ]' D$ V$ g5 H( a
89. if DB(j-1,6)==1
    5 H* t2 y! p, {- _. [
90. Yi=1/DB(j-1,4)*eye(3);
    : |  Q/ H7 ]1 T; z- |8 q
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
    
97. YT=1/real(DB(i,4))+j*(1/imag(DB(i,4)));
    , W% a3 b7 O8 y. L
98. YI=YT*eye(3);
    
99. YII=YT*eye(3);
    
100. YIII=-YT*eye(3);
     , ^8 R! q* q) ?  ~, q# m3 O
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)]');
     + y0 M2 G5 e; v/ C2 K; v
102. U(j,3*k+1)=H(1,1);
     
103. U(j,3*k+2)=H(2,1);
     
104. U(j,3*k+3)=H(3,1);
     + c5 Y2 N# K  L4 A
105. end
     4 V, e7 m' Q# ~  @# U8 E
106. end
     7 f+ |( p. u8 E  r, i5 e2 y1 V7 t
107. dU1=U(:,3*k+3)-U(:,3*k);
     0 s3 s8 ?1 B& W6 k; f* k
108. dU2=U(:,3*k+2)-U(:,3*k-1);
     
109. dU3=U(:,3*k+1)-U(:,3*k-2);
     # z' J8 H9 P- A6 b# e
110. if (max(abs(dU1))<0.001&&max(abs(dU2))<0.001&&max(abs(dU2))<0.001)
     5 x$ s% ]6 j! q5 e$ q5 g6 k- w4 E
111. break;
     
112. else            
     
113. k=k+1;
     
114. end
     
115. end
     
116. U
     ) ^0 w* f' j5 Z) t9 x
117. IL
     
118. IKj
     - P- w5 P: a6 `0 G* F9 r) S1 l. x
119. IKi
     ( E$ J7 v1 I9 |+ o
120. k
     

复制代码

sun710504

学习学习