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

1、影响电力系统暂态稳定性的因素1 研究对象1.1系统模型研究影响系统暂态稳定性的因素，首先要确定一个系统模型，本例选取参数可靠的美国西部电网等值模型SWCC-9，该系统为三机九节点系统，如图1.1所示。230kV125MW50MVar0.0085+j0.072B/2=j0.0745AuJ二ozzb58UoJ+OAUo8SUoJ二ozzb100MW35MVar.0.0119+j0.1008B/2=j0.10458nmnuj+Q3UoQzuoJ+71AUo9j0.05853230/13-813.8kVQ71nuj二2/_230kV90MW30MVar230kV16.5kV图1.1WSCC-9系统模型

2、1.2系统参数1.2.1节点参数表1.1节点已知参数节点类型电压幅值电压角度发电机有功发电机无功负荷有功负荷无功1V01.04000.71640.2705002PV1.0251.63000.0665003PV1.0250.8500-0.1086004PQ00005PQ001.25000.50006PQ000.90000.30007PQ00008PQ001.00000.35009PQ0000上表中发电机有功、无功出力和负荷的有功无功功率均为以100MVA为基准时的标幺值。1.2.2支路参数表1.2支路参数首节点末节点电阻电抗电纳半450.01000.08500.0880460.01700.092

3、00.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) 输入的机械功率保持恒定；(2) 忽略阻尼效应；(3)负荷采用恒阻抗模型。表1.3发电机参数发电机XdXH10.14

4、600.060823.6420.89580.11986.4031.31250.18133.01以上阻抗参数均以标幺值表示，额定转速下存储的能量（H）也转化为以100MVA为基准的标幺值。2 待研究的影响因素及其仿真流程设计2.1 故障切除时间的影响由等面积法则可在理论上分析得出如下结论：系统发生短路故障后，故障切除的越快，越有利于提高系统的暂态稳定性，且存在临界切除时间。假设系统0时刻以前处于稳定运行状态，0时刻在线路5-7上靠近7的母线出口处发生三相短路，故障切除时间为tc，仿真三台发电机的转速与功角变化曲线，通过观察曲线判断系统是否失去稳定。不断改变故障时间，求出恰好能使系统维持稳定的故障

5、切除时间，该时间就是临界切除时间。2.2 三相快速重合闸的影响由等面积法则分析可知，在系统发生短路故障后，在减速过程中自动重合闸动作将会增大加速面积，可能使原来无法维持暂态稳定的系统保持稳定。但是重合闸要求的时间比较苛刻，如果重合闸在系统已经失去稳定才动作则对系统的暂态稳定起不到作用，因此要求时间通常很短。这时如果采用单相重合闸由于有潜供电流的影响，单相重合闸的动作时间不能太快，因此本例只考虑三相重合闸。在2.1的基础上，tar时刻自动重合闸动作，调整参数，做对比仿真实验，一次不投入重合闸，另一次投入重合闸，通过判断系统是否失去稳定来得出重合闸对系统暂态稳定性是否有影响。然后，固定故障切除时间

6、tc为某一个定值，改变tar的值并观察系统是否失去稳定来确定自动重合闸的最长动作延时。2.3 仿真流程设计首先要求出系统的稳态运行参数，即系统潮流分布，这些参数作为暂态过程的初始值。潮流计算采用Newton-Raphson迭代法，求出各个母线的电压和角度。2.3.1发电机初态在动态分析时，首先应将发电机和负荷用相应的模型等效。本文将发电机等效为二阶经典模型，将负荷等效为恒阻抗负荷。图2.1发电机等效模型本系统中具有三台发电机，因此会引入三个内节点。因此节点导纳矩阵Y将增广到12阶。可表示为：Y3x9Y9x9YY=3x312x12Y9x3在Y中:3x3在Y和Y中:3x99x3在Y中:9x91Y=

7、,i=1,2,3iijxdi1Y一,i二1,2,3iijxdi1Y二Y+,i二1,2,3iiiijxdiPQY=Y+#j*,i=5,6,8iiiiV2V2ii(2.1)(2.2)(2.3)(2.4)将网络等效在发电机内节点构成的网络中，可以得到一个3阶的降阶节点导纳矩阵，该矩阵可由以下关系解出:Y=Y-YY-iY(2.5)3x33x99x99x3对于故障前和故障后的降阶节点导纳矩阵，利用式(2.5)可以计算出对应的降阶节点导纳矩阵；对于故障中的节点导纳矩阵，在12阶的增广矩阵中去掉故障母线所在的那一行和那一列，利用下式计算降阶矩阵。Y=Y-YY-iY(2.6)3x33x88x88x3发电机初态

8、中包含初始的电压幅值和功角，根据图2.1，可以求解发电机1、2、3的初始状态。Ei=V.+j(P+jQ)x,/=1,2,3V力iE=EZ5iii2.3.2列写发电机动态方程发电机采用经典模型，其动态方程为：”d(P-P)卜merdt2Hd8卜-、dtir,i=1,2,3(2.7)(2.8)(2.9)三台发电机共有六个状态量，定义：123818283得到：(2.10)2.3.3求解发电机状态利用龙格库塔数值积分方法，对发电机状态方程积分，求解发电机的状态。认为在故障过程中，发电机内电压的幅值是不变的，只有功角改变；认为发电机的机械功率是不变的，其值等于故障前的发电机电磁功率。3 仿真过程与结果3

9、.1潮流计算结果表3.1潮流计算结果节点类型电压幅值电压角度发电机有功发电机无功负荷有功负荷无功1V01.04000.00000.71640.2705002PV1.02509.28001.63000.0665003PV1.02504.66480.8500-0.1086004PQ1.0258-2.216800005PQ0.9956-3.9888001.25000.50006PQ1.0127-3.6874000.90000.30007PQ1.02583.719700008PQ1.01590.7275001.00000.35009PQ1.03241.966700003.2发电机初态表3.2发电机初态

10、发电机内电压功角11.05662.271621.050219.731631.017013.16643.3故障切除时间影响暂态稳定性仿真假设0时刻发生故障，经过5工频周期后切除故障，即tc=0.083S,重合闸不投入。得到三台发电机的转速和功角变化曲线。图3.1发电机转速变化曲线图3.1和图3.2分别为该过程中三台发电机的转速变化曲线和功角变化曲线，仿真总时间为4s。从曲线中可以看出系统没有失去暂态稳定。现在改变故障切除时间tc为12个工频周期，即0.2s,三台发电机的动态曲线分别如下。图3.3和图3.4分别为该过程中三台发电机的转速变化曲线和功角变化曲线，仿真总时间为4s。从曲线中可以明显看出

11、系统已经失去暂态稳定。因此可以得出故障切除时间会影响系统的暂态稳定性，并且时间越短越有利于提高系统的暂态稳定性。通过改变故障切除时间，可以求出临界故障切除时间tcc=0.16s，这时的发电机暂态曲线如下：图3.6临界切除时间下的发电机功角变化曲线3.4重合闸影响暂态稳定性仿真图3.7tc=0.17s未投入重合闸时发电机转速变化曲线取故障切除时间tc=0.17s,由3.3节可知，该时间大于临界故障切除时间，因此暂态稳定不能维持，发电机暂态曲线如图3.7和图3.8所示。图3.9tc=0.17s投入重合闸时发电机转速变化曲线图3.10tc=0.17s投入重合闸时发电机功角变化曲线现将自动重合闸投入，

12、重合闸动作延时为0.25s,即动作时刻为tar=0.42s,这时发电机的转速和功角变化曲线分别如图3.9和图3.10所示。从图中课明显看出此时系统没有失去暂态稳定性，说明三相快速重合闸能够提高系统的暂态稳定性。附录MATLAB仿真程序代码：%生成节点导纳矩阵Ypofl=zeros(9)forn=1:9Ypof(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)Ypof(BRANCH(n,2),BRANCH(n,2)=Ypofl(BRANCH(n,2)

13、,BRANCH(n,2)+1j*BRANCH(n,5)+1/(BRANCH(n,3)+1j*BRANCH(n,4)Ypof(BRANCH(n,1),BRANCH(n,2)=Ypofl(BRANCH(n,1),BRANCH(n,2)1/(BRANCH(n,3)+1j*BRANCH(n,4)Ypof(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-10tol=1LoopCount=0delta_y=z

14、eros(14,)H=zeros(8)N=zeros(8,6)M=zeros(6,8)L=zeros(6)%计算潮流whiletoltolerance%Delta_yform=2:9sun=0forn=1:9su=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_(m1)=BUS(m,5)BUS(m,7)sumendform=4:9sun=0forn=1:9su=sum+BUS(m,3)*BUS(n,3)*(G(m,n)*sin(BUS(m,4)BUS(n,4)B(m,n

15、)*cos(BUS(m,4)BUS(n,4)enddelta(m+5)=BUS(m,6)BUS(m,8)sumend%Hform=2:9forn=2:9ifm=ns=m0fork=1:9s=Lmum+BUS(m,3)*BUS(k,3)*(G(m,k)*sin(BUS(m,4)BUS(k,4)B(m,k)*cos(BUS(m,4)BUS(k,4)end(mH1,m1)=BUS(m,3)A2*B(m,m)+sumelse(mH1,n1)=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)endenden

16、d%Nform=2:9forn=4:9ifm=n%行列不等时(mN1,n3)=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)endendendform=2:7sum=0forn=1:9su=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)end(m1,m1)BUS(m,3)A2*G(m,m)sumendform=4:9forn=2:9ifm=ns=m0forBUS(k,4)+B(m,k)

17、ends=msum+BUS(m,3)*BUS(k,3)*sin(BUS(m,4)BUS(k,4)*(G(m,k)*cos(BUS(m,4)(mM3,m3)=BUS(m,3)A2*G(m,m)elsesum(mM3,n1)=BUS(m,3)*BUS(n,3)*BUS(n,4)+B(m,n)*sin(BUS(m,4)BUS(n,4)end(G(m,n)*cos(BUS(m,4)endendform=4:9forn=4:9ifm=ns=m0forBUS(k,4)B(m,k)ends=msum+BUS(m,3)*BUS(k,3)*cos(BUS(m,4)BUS(k,4)*sin(BUS(m,4)(mL

18、3,m3)=BUS(m,3)A2*B(m,m)elsesum(mL3,n3)=BUS(m,3)*BUS(n,3)BUS(n,4)B(m,n)*cos(BUS(m,4)BUS(n,4)end*(G(m,n)*sin(BUS(m,4)endend%雅可比矩阵JACOB=HNML%修正量delta_x=JACOBdelta_y%修正电压幅值和相位form=4:9BU(m,3)=BUS(m,)+delta_x(m+)*BUS(m,3)endform=2:9BU(m,4)=BUS(m,)+delta_x(m1)end%计算精度tol=abs(max(delta_y)LoopCount=LoopCount

19、+1endBUS(:,4)=BUS(:,4)*180/pidisp(POWERFLOW,date)disp(BUS)%故障分析%简化的Y阵%故障前%计算增广Y阵Y_pf33=zeros()Y_pf39=zeros(3,9)Y_pf99=YpoflY_pf33(1,)=Y pf33(2,)sY pf33(3,)=Y_pf39(1,)=Y pf39(2,)sY pf39(3,)=1/(1j*GEN(1,3)=1/(1j*GEN(2,3)=1/(1j*GEN(3,3)=1/(1j*GEN(1,3)=1/(1j*GEN(2,3)=1/(1j*GEN(3,3)Y_pf93=transpose(Y_pf3

20、9)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*GEN(3,3)Y_pf99(5,5):=Y_pf99(5,5)+BUS(5,7)/BUS(5,3)A21j*BUS(5,8)/BUS(5,3)A2Ypf99(6,6)s=Y_pf99(6,6)+BUS(6,7)/BUS(6,3)A21j*BUS(6,8)/BUS(6,3)A2Ypf99(8,8)s=Y_pf99(8,8)+BUS(8,7)/BUS(8,3)A21j*BUS(

21、8,8)/BUS(8,3)A2Yex=Y_pf33Y_pf39Y_pf93Y_pf99Yrpf=Y_pf33Y_pf39/Y_pf99*Y_pf93%故障中Y_df33=Y_pf33Y_df38=Y_pf39Y_df38(:,)=Y_df83=transpose(Y_df38)Y_df88=Y_pf99Y_df88(7,:)=Y_df88(:,)=Yrdf=Y_df33Y_df38/Y_df88*Y_df83%故障后Y_af33=Y_pf33Y_af39=Y_pf39Y_af93=Y_pf93Y_af99=Y_pf99Y_af99(5,)=Y_af99(5,5)Y_af99(7,)=Y_af9

22、9(7,7)1j*BRANCH(3,5)1j*BRANCH(3,5)1/(BRANCH(3,3)+1j*BRANCH(3,4)1/(BRANCH(3,3)+1j*BRANCH(3,4)Y_af99(5,)=0Y_af99(7,)=0Yraf=Y_af33Y_af39/Y_af99*Y_af93disp(Beforefault,date)disp(Yrpf)disp(Yrdf)disp(Yraf)%故障前发电机状态tempAng=BUS(1:3,4)*pi/180tempV1=BUS(1:3,3).*exp(1j.*tempAng)tempV2=tempV1+conj(BUS(1:3,5)+1j

23、.*BUS(1:3,6)./tempV1).*1j.*GEN(1:3,3)E_pf=sqrt(real(tempV2).A?+imag(tempV2).)Delta_pf=atan(imag(tempV2)./real(tempV2).*180./piGEN(:,5)=real(tempV2.*conj(Yrpf*tempV2)disp(GenInit,date)disp(E)disp(E_pf)disp(delta)disp(Delta_pf)%故障中发电机状态t_cut=0.17t_end=4t_ar=0.42%故障切除时间%仿真结束时间%重合闸时间odetol=1e-10options=odeset(RelTol,1e-10)tspan=0,t_cutDelta_init=Delta_pf.*pi./180x0_df=ones(,1).*2*pi.*60Delta_initTdf_out,Xdf_out=ode45(t,x)Gen_Fun(t,x,GEN,Yrdf,E_pf),tspan,