故障切除時(shí)間與三相快速重合閘對電力系統(tǒng)穩(wěn)定性的影響仿真(含MATLAB程序)

1、影響電力系統(tǒng)暫態(tài)穩(wěn)定性的因素1 研究對象1.1系統(tǒng)模型研究影響系統(tǒng)暫態(tài)穩(wěn)定性的因素，首先要確定一個(gè)系統(tǒng)模型，本例選取參數(shù)可靠的美國西部電網(wǎng)等值模型SWCC-9，該系統(tǒng)為三機(jī)九節(jié)點(diǎn)系統(tǒng)，如圖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系統(tǒng)模型

2、1.2系統(tǒng)參數(shù)1.2.1節(jié)點(diǎn)參數(shù)表1.1節(jié)點(diǎn)已知參數(shù)節(jié)點(diǎn)類型電壓幅值電壓角度發(fā)電機(jī)有功發(fā)電機(jī)無功負(fù)荷有功負(fù)荷無功1V01.04000.71640.2705002PV1.0251.63000.0665003PV1.0250.8500-0.1086004PQ00005PQ001.25000.50006PQ000.90000.30007PQ00008PQ001.00000.35009PQ0000上表中發(fā)電機(jī)有功、無功出力和負(fù)荷的有功無功功率均為以100MVA為基準(zhǔn)時(shí)的標(biāo)幺值。1.2.2支路參數(shù)表1.2支路參數(shù)首節(jié)點(diǎn)末節(jié)點(diǎn)電阻電抗電納半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上表中所有的參數(shù)均為標(biāo)幺值，對于變壓器支路。最后三行表示三臺變壓器參數(shù)已經(jīng)計(jì)算出變壓器的等效電抗并直接在表格中給出。1.2.3發(fā)電機(jī)參數(shù)對于發(fā)電機(jī)，采用二階經(jīng)典模型，并對系統(tǒng)作如下假設(shè)(1) 輸入的機(jī)械功率保持恒定；(2) 忽略阻尼效應(yīng)；(3)負(fù)荷采用恒阻抗模型。表1.3發(fā)電機(jī)參數(shù)發(fā)電機(jī)XdXH10.14

4、600.060823.6420.89580.11986.4031.31250.18133.01以上阻抗參數(shù)均以標(biāo)幺值表示，額定轉(zhuǎn)速下存儲的能量（H）也轉(zhuǎn)化為以100MVA為基準(zhǔn)的標(biāo)幺值。2 待研究的影響因素及其仿真流程設(shè)計(jì)2.1 故障切除時(shí)間的影響由等面積法則可在理論上分析得出如下結(jié)論：系統(tǒng)發(fā)生短路故障后，故障切除的越快，越有利于提高系統(tǒng)的暫態(tài)穩(wěn)定性，且存在臨界切除時(shí)間。假設(shè)系統(tǒng)0時(shí)刻以前處于穩(wěn)定運(yùn)行狀態(tài)，0時(shí)刻在線路5-7上靠近7的母線出口處發(fā)生三相短路，故障切除時(shí)間為tc，仿真三臺發(fā)電機(jī)的轉(zhuǎn)速與功角變化曲線，通過觀察曲線判斷系統(tǒng)是否失去穩(wěn)定。不斷改變故障時(shí)間，求出恰好能使系統(tǒng)維持穩(wěn)定的故障

5、切除時(shí)間，該時(shí)間就是臨界切除時(shí)間。2.2 三相快速重合閘的影響由等面積法則分析可知，在系統(tǒng)發(fā)生短路故障后，在減速過程中自動重合閘動作將會增大加速面積，可能使原來無法維持暫態(tài)穩(wěn)定的系統(tǒng)保持穩(wěn)定。但是重合閘要求的時(shí)間比較苛刻，如果重合閘在系統(tǒng)已經(jīng)失去穩(wěn)定才動作則對系統(tǒng)的暫態(tài)穩(wěn)定起不到作用，因此要求時(shí)間通常很短。這時(shí)如果采用單相重合閘由于有潛供電流的影響，單相重合閘的動作時(shí)間不能太快，因此本例只考慮三相重合閘。在2.1的基礎(chǔ)上，tar時(shí)刻自動重合閘動作，調(diào)整參數(shù)，做對比仿真實(shí)驗(yàn)，一次不投入重合閘，另一次投入重合閘，通過判斷系統(tǒng)是否失去穩(wěn)定來得出重合閘對系統(tǒng)暫態(tài)穩(wěn)定性是否有影響。然后，固定故障切除時(shí)間

6、tc為某一個(gè)定值，改變tar的值并觀察系統(tǒng)是否失去穩(wěn)定來確定自動重合閘的最長動作延時(shí)。2.3 仿真流程設(shè)計(jì)首先要求出系統(tǒng)的穩(wěn)態(tài)運(yùn)行參數(shù)，即系統(tǒng)潮流分布，這些參數(shù)作為暫態(tài)過程的初始值。潮流計(jì)算采用Newton-Raphson迭代法，求出各個(gè)母線的電壓和角度。2.3.1發(fā)電機(jī)初態(tài)在動態(tài)分析時(shí)，首先應(yīng)將發(fā)電機(jī)和負(fù)荷用相應(yīng)的模型等效。本文將發(fā)電機(jī)等效為二階經(jīng)典模型，將負(fù)荷等效為恒阻抗負(fù)荷。圖2.1發(fā)電機(jī)等效模型本系統(tǒng)中具有三臺發(fā)電機(jī)，因此會引入三個(gè)內(nèi)節(jié)點(diǎn)。因此節(jié)點(diǎn)導(dǎo)納矩陣Y將增廣到12階?？杀硎緸椋篩3x9Y9x9YY=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)將網(wǎng)絡(luò)等效在發(fā)電機(jī)內(nèi)節(jié)點(diǎn)構(gòu)成的網(wǎng)絡(luò)中，可以得到一個(gè)3階的降階節(jié)點(diǎn)導(dǎo)納矩陣，該矩陣可由以下關(guān)系解出:Y=Y-YY-iY(2.5)3x33x99x99x3對于故障前和故障后的降階節(jié)點(diǎn)導(dǎo)納矩陣，利用式(2.5)可以計(jì)算出對應(yīng)的降階節(jié)點(diǎn)導(dǎo)納矩陣；對于故障中的節(jié)點(diǎn)導(dǎo)納矩陣，在12階的增廣矩陣中去掉故障母線所在的那一行和那一列，利用下式計(jì)算降階矩陣。Y=Y-YY-iY(2.6)3x33x88x88x3發(fā)電機(jī)初態(tài)

8、中包含初始的電壓幅值和功角，根據(jù)圖2.1，可以求解發(fā)電機(jī)1、2、3的初始狀態(tài)。Ei=V.+j(P+jQ)x,/=1,2,3V力iE=EZ5iii2.3.2列寫發(fā)電機(jī)動態(tài)方程發(fā)電機(jī)采用經(jīng)典模型，其動態(tài)方程為：”d(P-P)卜merdt2Hd8卜-、dtir,i=1,2,3(2.7)(2.8)(2.9)三臺發(fā)電機(jī)共有六個(gè)狀態(tài)量，定義：123818283得到：(2.10)2.3.3求解發(fā)電機(jī)狀態(tài)利用龍格庫塔數(shù)值積分方法，對發(fā)電機(jī)狀態(tài)方程積分，求解發(fā)電機(jī)的狀態(tài)。認(rèn)為在故障過程中，發(fā)電機(jī)內(nèi)電壓的幅值是不變的，只有功角改變；認(rèn)為發(fā)電機(jī)的機(jī)械功率是不變的，其值等于故障前的發(fā)電機(jī)電磁功率。3 仿真過程與結(jié)果3

9、.1潮流計(jì)算結(jié)果表3.1潮流計(jì)算結(jié)果節(jié)點(diǎn)類型電壓幅值電壓角度發(fā)電機(jī)有功發(fā)電機(jī)無功負(fù)荷有功負(fù)荷無功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發(fā)電機(jī)初態(tài)表3.2發(fā)電機(jī)初態(tài)

10、發(fā)電機(jī)內(nèi)電壓功角11.05662.271621.050219.731631.017013.16643.3故障切除時(shí)間影響暫態(tài)穩(wěn)定性仿真假設(shè)0時(shí)刻發(fā)生故障，經(jīng)過5工頻周期后切除故障，即tc=0.083S,重合閘不投入。得到三臺發(fā)電機(jī)的轉(zhuǎn)速和功角變化曲線。圖3.1發(fā)電機(jī)轉(zhuǎn)速變化曲線圖3.1和圖3.2分別為該過程中三臺發(fā)電機(jī)的轉(zhuǎn)速變化曲線和功角變化曲線，仿真總時(shí)間為4s。從曲線中可以看出系統(tǒng)沒有失去暫態(tài)穩(wěn)定。現(xiàn)在改變故障切除時(shí)間tc為12個(gè)工頻周期，即0.2s,三臺發(fā)電機(jī)的動態(tài)曲線分別如下。圖3.3和圖3.4分別為該過程中三臺發(fā)電機(jī)的轉(zhuǎn)速變化曲線和功角變化曲線，仿真總時(shí)間為4s。從曲線中可以明顯看出

11、系統(tǒng)已經(jīng)失去暫態(tài)穩(wěn)定。因此可以得出故障切除時(shí)間會影響系統(tǒng)的暫態(tài)穩(wěn)定性，并且時(shí)間越短越有利于提高系統(tǒng)的暫態(tài)穩(wěn)定性。通過改變故障切除時(shí)間，可以求出臨界故障切除時(shí)間tcc=0.16s，這時(shí)的發(fā)電機(jī)暫態(tài)曲線如下：圖3.6臨界切除時(shí)間下的發(fā)電機(jī)功角變化曲線3.4重合閘影響暫態(tài)穩(wěn)定性仿真圖3.7tc=0.17s未投入重合閘時(shí)發(fā)電機(jī)轉(zhuǎn)速變化曲線取故障切除時(shí)間tc=0.17s,由3.3節(jié)可知，該時(shí)間大于臨界故障切除時(shí)間，因此暫態(tài)穩(wěn)定不能維持，發(fā)電機(jī)暫態(tài)曲線如圖3.7和圖3.8所示。圖3.9tc=0.17s投入重合閘時(shí)發(fā)電機(jī)轉(zhuǎn)速變化曲線圖3.10tc=0.17s投入重合閘時(shí)發(fā)電機(jī)功角變化曲線現(xiàn)將自動重合閘投入，

12、重合閘動作延時(shí)為0.25s,即動作時(shí)刻為tar=0.42s,這時(shí)發(fā)電機(jī)的轉(zhuǎn)速和功角變化曲線分別如圖3.9和圖3.10所示。從圖中課明顯看出此時(shí)系統(tǒng)沒有失去暫態(tài)穩(wěn)定性，說明三相快速重合閘能夠提高系統(tǒng)的暫態(tài)穩(wěn)定性。附錄MATLAB仿真程序代碼：%生成節(jié)點(diǎn)導(dǎo)納矩陣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)%計(jì)算潮流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%行列不等時(shí)(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%計(jì)算精度tol=abs(max(delta_y)LoopCount=LoopCount

19、+1endBUS(:,4)=BUS(:,4)*180/pidisp(POWERFLOW,date)disp(BUS)%故障分析%簡化的Y陣%故障前%計(jì)算增廣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)%故障前發(fā)電機(jī)狀態(tài)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)%故障中發(fā)電機(jī)狀態(tài)t_cut=0.17t_end=4t_ar=0.42%故障切除時(shí)間%仿真結(jié)束時(shí)間%重合閘時(shí)間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,