故障切除时间与三相快速重合闸对电力系统稳定性的影响仿真含MATLAB程序教材

1、影响电力系统暂态稳定性的因素1 研究对象1.1 系统模型研究影响系统暂态稳定性的因素，首先要确定一个系统模型，本例选取参数可靠的美国西部电网等值模型 SWCC-9，该系统为三机九节点系统，如图 1.1 所示100MW18kV 18/2300.0085 + j0.072B/2 = j0.07459 j0.0585 3230/13.8 13.8kV35MVar0.0119 + j0.1008B/2 = j0.10458230kV161.0j+230.351.0j=2/B071.0j+930.0.971.0230kV125MW50MVar880.0j=2/970.0j+710.90MW30MVar2

2、30kV6750.032/5.616.5kV 1图 1.1 WSCC-9 系统模型1.2 系统参数1.2.1 节点参数表 1.1 节点已知参数节点类型电压幅值电压角度发电机有功发电机无功负荷有功负荷无功1V1.04000.71640.2705002PV1.0251.63000.0665003PV1.0250.8500-0.1086004PQ00005PQ001.25000.50006PQ000.90000.30007PQ00008PQ001.00000.35009PQ0000上表中发电机有功、无功出力和负荷的有功无功功率均为以 100MVA 为基准时的标幺值。1.2.2 支路参数表 1.2 支

3、路参数首节点末节点电阻电抗电纳一半450.01000.08500.0880460.01700.09200.0790570.03200.16100.1530690.03900.17000.1790780.00850.07200.0745890.01190.10080.1045140.00000.05760.0000270.00000.06250.0000390.00000.05860.0000上表中所有的参数均为标幺值，对于变压器支路。最后三行表示三台变压器参数， 已经计算出变压器的等效电抗并直接在表格中给出。1.2.3 发电机参数 对于发电机，采用二阶经典模型，并对系统作如下假设：(1) 输入

4、的机械功率保持恒定；(2) 忽略阻尼效应；(3) 负荷采用恒阻抗模型。表 1.3 发电机参数发电机XdXdH10.14600.060823.6420.89580.11986.4031.31250.18133.01以上阻抗参数均以标幺值表示，额定转速下存储的能量（ H）也转化为以 100MVA 为基准的标幺值。2 待研究的影响因素及其仿真流程设计2.1 故障切除时间的影响由等面积法则可在理论上分析得出如下结论：系统发生短路故障后，故障切除的越 快，越有利于提高系统的暂态稳定性，且存在临界切除时间。假设系统 0 时刻以前处于稳定运行状态， 0 时刻在线路 5-7 上靠近 7 的母线出口处 发生三相

5、短路，故障切除时间为 tc，仿真三台发电机的转速与功角变化曲线，通过观察 曲线判断系统是否失去稳定。不断改变故障时间，求出恰好能使系统维持稳定的故障切 除时间，该时间就是临界切除时间。2.2 三相快速重合闸的影响由等面积法则分析可知，在系统发生短路故障后，在减速过程中自动重合闸动作将 会增大加速面积，可能使原来无法维持暂态稳定的系统保持稳定。但是重合闸要求的时 间比较苛刻，如果重合闸在系统已经失去稳定才动作则对系统的暂态稳定起不到作用， 因此要求时间通常很短。这时如果采用单相重合闸由于有潜供电流的影响，单相重合闸 的动作时间不能太快，因此本例只考虑三相重合闸。在 2.1 的基础上， tar 时

6、刻自动重合 闸动作，调整参数，做对比仿真实验，一次不投入重合闸，另一次投入重合闸，通过判 断系统是否失去稳定来得出重合闸对系统暂态稳定性是否有影响。然后，固定故障切除 时间 tc为某一个定值，改变 tar 的值并观察系统是否失去稳定来确定自动重合闸的最长动 作延时。2.3 仿真流程设计首先要求出系统的稳态运行参数，即系统潮流分布，这些参数作为暂态过程的初始值。潮流计算采用 Newton-Raphson 迭代法，求出各个母线的电压和角度。2.3.1 发电机初态在动态分析时，首先应将发电机和负荷用相应的模型等效。本文将发电机等效为二阶经典模型，将负荷等效为恒阻抗负荷。本系统中具有三台发电机， 因此

7、会引入三个内节点。 因此节点导纳矩阵 Y 将增广到12 阶。可表示为：Y12 12Y 3 3Y93Y39Y9 9在 Y 3 3 中：在 Y 3 9 和 Y 9 3 中：在 Y 9 9 中：Yii jx1 ,i 1,2,3 jxdiiiYiijx1 ,i 1,2,3jxdi1Yii Yii1 ,i 1,2,3jxdiYii Yii VPii2 jVQi2i ,i 5,6,8(2.1)(2.2)(2.3)(2.4)将网络等效在发电机内节点构成的网络中， 可以得到一个 3 阶的降阶节点导纳矩阵，该矩阵可由以下关系解出：Y3 3 Y3 9Y9 9 Y9 3(2.5)对于故障前和故障后的降阶节点导纳矩阵

8、，利用式 (2.5)可以计算出对应的降阶节点 导纳矩阵；对于故障中的节点导纳矩阵，在 12 阶的增广矩阵中去掉故障母线所在的那 一行和那一列，利用下式计算降阶矩阵。Y Y3 3 Y38Y818Y8 3(2.6)发电机初态中包含初始的电压幅值和功角，根据图 2.1，可以求解发电机 1、2、3 的初始状态。EiP jQ * Vi j( ii)* xdi,i 1,2,3Vi(2.7)Ei Ei i (2.8)2.3.2 列写发电机动态方程发电机采用经典模型，其动态方程为：didtdidtPm Per2H,i 1,2,3(2.9)三台发电机共有六个状态量，定义：3123得到：(2.10)x F t,x

9、2.3.3 求解发电机状态利用龙格库塔数值积分方法，对发电机状态方程积分，求解发电机的状态。认为在故障过程中，发电机内电压的幅值是不变的，只有功角改变；认为发电机的机械功率是不变的，其值等于故障前的发电机电磁功率3 仿真过程与结果3.1 潮流计算结果表 3.1 潮流计算结果节点类型电压幅值电压角度发电机有功发电机无功负荷有功负荷无功1V1.04000.00000.71640.2705002PV1.02509.28001.63000.0665003PV1.02504.66480.8500-0.1086004PQ1.0258-2.216800005PQ0.9956-3.9888001.25000.

10、50006PQ1.0127-3.6874000.90000.30007PQ1.02583.719700008PQ1.01590.7275001.00000.35009PQ1.03241.966700003.2 发电机初态表 3.2 发电机初态发电机内电压功角11.05662.271621.050219.731631.017013.16643.3 故障切除时间影响暂态稳定性仿真假设 0时刻发生故障，经过 5工频周期后切除故障， 即tc = 0.083 s，重合闸不投入得到三台发电机的转速和功角变化曲线。3883863843823803783760图 3.1 发电机转速变化曲线1400120010

11、00800600400200图 3.2 发电机功角变化曲线图 3.1 和图 3.2 分别为该过程中三台发电机的转速变化曲线和功角变化曲线，仿真总时间为 4 s。从曲线中可以看出系统没有失去暂态稳定。现在改变故障切除时间tc 为12个工频周期，即 0.2 s，三台发电机的动态曲线分别如下。图 3.3 发电机转速变化曲线图 3.4 发电机功角变化曲线图 3.3 和图 3.4 分别为该过程中三台发电机的转速变化曲线和功角变化曲线，仿真 总时间为 4 s。从曲线中可以明显看出系统已经失去暂态稳定。因此可以得出故障切除时间会影响系统的暂态稳定性，并且时间越短越有利于提高系统的暂态稳定性。通过改变故障切除

12、时间，可以求出临界故障切除时间 tcc = 0.16 s，这时的发电机暂 态曲线如下：)s/dar(,2,10.5 1 1.5 2 2.5 3 3.5 4 t (s)图 3.5 临界切除时间下的发电机转速变化曲线)geD(3,2,150000 0.5 11.52t (s)2.53 3.5 4图 3.6 临界切除时间下的发电机功角变化曲线3.4 重合闸影响暂态稳定性仿真)s/dar(0.51 1.5 2 2.5 3 3.5 4 t (s)3600图 3.7 tc = 0.17 s 未投入重合闸时发电机转速变化曲线图 3.8 tc = 0.17 s 未投入重合闸时发电机功角变化曲线取故障切除时间

13、tc = 0.17 s，由 3.3节可知，该时间大于临界故障切除时间，因此暂 态稳定不能维持，发电机暂态曲线如图 3.7 和图 3.8 所示。图 3.9 tc = 0.17 s 投入重合闸时发电机转速变化曲线( t图 3.10 tc = 0.17 s 投入重合闸时发电机功角变化曲线现将自动重合闸投入， 重合闸动作延时为 0.25 s，即动作时刻为 tar = 0.42 s，这时发 电机的转速和功角变化曲线分别如图 3.9 和图 3.10 所示。从图中课明显看出此时系统没 有失去暂态稳定性，说明三相快速重合闸能够提高系统的暂态稳定性。附录MATLAB 仿真程序代码：%生成节点导纳矩阵Ypofl

14、= zeros ( 9);for n = 1: 9Ypofl( BRANCH( n, 1), BRANCH( n, 1) = Ypofl ( BRANCH( n, 1), BRANCH( n, 1) +1j * BRANCH( n, 5) + 1/( BRANCH( n, 3) + 1j * BRANCH( n, 4);Ypofl( BRANCH( n, 2), BRANCH( n, 2) = Ypofl ( BRANCH( n, 2), BRANCH( n, 2) +1j * BRANCH( n, 5) + 1/( BRANCH( n, 3) + 1j * BRANCH( n, 4);Ypo

15、fl( BRANCH( n, 1), BRANCH( n, 2)= Ypofl( BRANCH( n, 1),BRANCH( n, 2) -1/( BRANCH( n, 3) + 1j * BRANCH( n, 4);Ypofl( BRANCH( n, 2), BRANCH( n, 1)= Ypofl( BRANCH( n, 2),BRANCH( n, 1) -1/( BRANCH( n, 3) + 1j * BRANCH( n, 4);endG = real( Ypofl );B = imag( Ypofl );tolerance= 1e-10 ;tol = 1;LoopCount= 0;d

16、elta_y= zeros ( 14 , 1);H = zeros( 8);N = zeros(8, 6);M = zeros(6, 8);L = zeros( 6);%计算潮流while tol tolerance% Delta_yfor m= 2: 9sum= 0;for n = 1: 9sum= sum + BUS ( m, 3) * BUS(n, 3) *( G( m, n) * cos( BUS( m, 4) - BUS ( n, 4)+ B ( m, n) * sin ( BUS( m, 4) - BUS ( n, 4);enddelta_y ( m- 1) = BUS ( m,

17、5) - BUS ( m, 7) - sum endfor m4: 9sumfor0;= 1: 9sum- B(m, n) *endcos= sum + BUS( m, 3) * BUS( ( BUS( m, 4) - BUS ( n, 4);n, 3) *(G(m,n) * sin ( BUS( m, 4) - BUS(n, 4)delta_yend( m+5) = BUS ( m, 6) - BUS( m, 8)sum ;% Hfor m = 2: 9for n = 2: 9 if m = n sum = 0;for k = 1: 9sum =sum +BUS(m,3) * BUS(k,3

18、) * (G(m,k) * sin ( BUS( m, 4) - BUS(k, 4) - B(m,k) * cos ( BUS( m, 4) - BUS ( k , 4);endH ( m- 1, m- 1) = BUS ( m, 3) 2 * B ( m, m) + sum ;BUS ( n, 4)elseH( m- 1, n- 1)= - BUS ( m, 3) * BUS(n,3) * (G(m, n) * sin ( BUS( m, 4)BUS( n, 4)- B( m, n) * cos( BUS( m, 4) - BUS ( n, 4);endendend% Nfor m =2:9

19、forn= 4: 9ifm = n%行列不等时N( m- 1, n- 3)= - BUS ( m, 3) * BUS(n,3) * (G(m, n) * cos ( BUS( m, 4)BUS( n, 4)+ B( m, n) * sin( BUS( m, 4) - BUS ( n, 4);endendendfor m =2:7sum=0;forn= 1: 9sum= sum+ BUS ( m, 3) * BUS(n, 3) * (G(m,n)cos ( BUS( m, 4)+ B ( m, n) * sin ( BUS( m, 4) - BUS ( n, 4);end( m- 1, m- 1

20、) = - BUS ( m, 3) 2 * G ( m, m) - sumendfor m4: 9for2: 9ifm = nsumfor k0;= 1: 9sumBUS( k, 4)( m, k) *end= sum + BUS ( m, 3) * BUS ( k, 3)( BUS( m, 4) - BUS ( k, 4);* ( G( m, k)* cos ( BUS( m, 4)elseBUS( n, 4)endsin( m- 3, m- 3)( m- 3, n- 1) + B ( m, n) * endsinendfor m4: 9for4: 9ifm = nsum0;sumfor k

21、1:BUS( k, 4)( m, k) *endelseBUS( n, 4)= BUS ( m, 3) 2 * G( m, m)= BUS ( m, 3) * BUS ( n, 3) * ( BUS( m, 4) - BUS ( n, 4);= sum + BUS ( m, 3) * BUS ( k, 3)( BUS( m, 4) - BUS ( k, 4);cos( m- 3, m- 3)= BUS ( m, 3) 2 * B ( m, m)( m- 3, n- 3)- B (m, n) * cos ( BUS( m, 4) - BUS (n, 4); end= - BUS ( m, 3)

22、* BUS ( n, 3)sum ;( G( m, n) * ( G( m, k)sum ;* ( G( m, n)cos ( BUS( m, 4)* sin ( BUS( m, 4)* sin ( BUS( m, 4)endend% 雅可比矩阵JACOB = H N ; M L ;% 修正量delta_x = - JACOB delta_y% 修正电压幅值和相位form = 4: 9BUS( m, 3) = BUS ( m, 3) + delta_x( m+5) * BUS (m,3);endform = 2: 9BUS( m, 4) = BUS ( m, 4) + delta_x( m-

23、1);end% 计算精度tol = abs ( max( delta_y ); LoopCount = LoopCount + 1; endBUS (:, 4) = BUS (:, 4) disp disp( POWERFLOW * 180 / pi ; , date );( BUS);故障分析 简化的 Y 阵 故障前 计算增广 Y 阵Y_pf33= zeros ( 3);Y_pf39= zeros ( 3, 9);Y_pf99= Ypofl ;Y_pf33(1, 1) = 1/( 1j * GEN( 1, 3);Y_pf33(2, 2) = 1/( 1j * GEN( 2, 3);Y_pf3

24、3(3, 3) = 1/( 1j * GEN( 3, 3);Y_pf39(1, 1) = -1/( 1j * GEN( 1, 3);Y_pf39(2, 2) = -1/( 1j * GEN( 2, 3);Y_pf39(3, 3) = -1/( 1j * GEN( 3, 3);Y_pf93= transpose ( Y_pf39 );Y_pf99(1, 1) = Y_pf99 (1,1) + 1/(1j * GEN( 1, 3);Y_pf99(2, 2) = Y_pf99 (2,2) + 1/(1j * GEN( 2, 3);Y_pf99(3, 3) = Y_pf99 (3,3) + 1/(1j

25、 * GEN( 3, 3);Y_pf99(5, 5) = Y_pf99 (5,5) + BUS( 5, 7)/ BUS( 5,Y_pf99(6, 6) = Y_pf99 (6,6) + BUS( 6, 7)/ BUS( 6,Y_pf99(8, 8) = Y_pf99 (8,8) + BUS( 8, 7)/ BUS( 8,Yex = Y_pf33 Y_pf39; Y_pf93 Y_pf99;Yrpf= Y_pf33 - Y_pf39 / Y_pf99* Y_pf93 ;% 故障中Y_df33= Y_pf33 ;Y_df38= Y_pf39 ;Y_df38(:, 7) = ;Y_df83= tra

26、nspose ( Y_df38 );Y_df88= Y_pf99 ;%3)3)3)1j1j1j* BUS( 5, 8)/* BUS( 6, 8)/* BUS( 8, 8)/BUS( 5, 3)BUS( 6, 3)BUS( 8, 3)2;2;2;Y_df88 ( 7,:) = ;Y_df88 (:, 7) = ;Yrdf = Y_df33 - Y_df38 / Y_df88 * Y_df83 % 故障后Y_af33 = Y_pf33Y_af39 = Y_pf39Y_af93 = Y_pf93Y_af99 = Y_pf99Y_af99(5,5) = Y_af99( 5, 5)- 1/(BRANCH

27、( 3, 3)+ 1j * BRANCH( 3, 4)- 1 j * BRANCH( 3 , 5);Y_af99(7,7) = Y_af99( 7, 7)- 1/(BRANCH( 3, 3)+ 1j * BRANCH( 3, 4)- 1 j * BRANCH( 3 , 5);Y_af99(5, 7) = 0;Y_af99(7, 5) = 0;Yraf= Y_af33 - Y_af39 /Y_af99* Y_af93 ;disp ( Before fault , date );disp ( Yrpf );disp ( Yrdf );disp ( Yraf );% 故障前发电机状态tempAng

28、= BUS ( 1: 3, 4)* pi / 180 ;tempV1 = BUS ( 1: 3, 3).* exp ( 1j .* tempAng );tempV2 = tempV1 + conj ( BUS(1:3, 5) + 1j .* BUS(1: 3, 6) ./ tempV1 ) .* 1j .* GEN( 1: 3, 3);E_pf = sqrt ( real ( tempV2 ). 2 + imag ( tempV2 ). 2);Delta_pf = atan ( imag ( tempV2 )./ real ( tempV2 ).* 180. / pi ;GEN (:, 5)

29、= real ( tempV2 .* conj ( Yrpf * tempV2 );disp( GenInit , datedisp( E );disp( E_pf );disp( delta );disp( Delta_pf );%故障切除时间%仿真结束时间%重合闸时间odetol options% 故障中发电机状态 t_cut = 0.17 ; t_end = 4; t_ar = 0.42 ;1e-10 ;= odeset ( RelTol , 1e-10 );Delta_init= Delta_pf.* pi ./180 ;x0_df = ones ( 3, 1).*2. *pi .*60 ;Delta_init Tdf_out, Xdf_out =ode45 ( t, x) Gen_Fun ( t, x,GEN,Yrdf, E_pfDelta_df= Xdf_out (:,4:6).*