MATLAB實(shí)驗(yàn)電力系統(tǒng)暫態(tài)穩(wěn)定分析

實(shí)驗(yàn)三電力系統(tǒng)暫態(tài)穩(wěn)定分析電力系統(tǒng)暫態(tài)穩(wěn)定計(jì)算實(shí)際上就是求解發(fā)電機(jī)轉(zhuǎn)子運(yùn)動(dòng)方程的初值問(wèn)題，從而得出0-t和3-t的關(guān)系曲線。每臺(tái)發(fā)電機(jī)的轉(zhuǎn)子運(yùn)動(dòng)方程是兩個(gè)一階非線性的常微分方程。因此，首先介紹常微分方程的初值問(wèn)題的數(shù)值解法。-、常微分方程的初值問(wèn)題（一）問(wèn)題及求解公式的構(gòu)造方法我們討論形如式（3-1）的一階微分方程的初值問(wèn)題（3-1）礦(x)=f(x,y),a<x<by(X0)=y0設(shè)初值問(wèn)題（3-1）的解為y（x），為了求其數(shù)值解而采取離散化方法，在求解區(qū)間a,b］上取一組節(jié)點(diǎn)a=x<x<…<x<x<…<x=b

01ii+（3-1）稱h=xi+1-xi（i=0,1,…,n-1）為步長(zhǎng)。在等步長(zhǎng)的情況下，步長(zhǎng)為7b一a

h=n用七表示在節(jié)點(diǎn)工,處解的準(zhǔn)確值y（xi）的近似值。設(shè)法構(gòu)造序列所滿足的一個(gè)方程（稱為差分方程）y=y+h叩（x,y,h）（3-2）i+1iii作為求解公式，這是一個(gè)遞推公式，從（x0，y0）出發(fā)，采用步進(jìn)方式，自左相右逐步算出y（x）在所有節(jié)點(diǎn)xi上的近似值yi（i=1,2,…,n）。在公式（3-2）中，為求卜^只用到前面一步的值yi,這種方法稱為單步法。在公式（3-2）中的y.1由y.明顯表示出，稱為顯式公式。而形如（3-3）y=y+h叩（x,y,y,h）（3-3）i+1iiii+1的公式稱為隱式公式，因?yàn)槠溆叶舜缰羞€包括y,+1。如果由公式求yi+1時(shí)，不止用到前一個(gè)節(jié)點(diǎn)的值，則稱為多步法。由式（3-1）可得dy=f（x,y）dx（3-4）兩邊在［x.,x.」上積分，得y（x）=y（x）+J氣+1f（x,y（x））dxi+1ixi由此可以看出，如果想構(gòu)造求解公式，就要對(duì)右端的積分項(xiàng)作某種數(shù)值處理。這種求解公式(3-5)的構(gòu)造方法叫做數(shù)值積分法。（二）一般的初值問(wèn)題的解法1.歐拉法和改進(jìn)歐拉法對(duì)于初值問(wèn)題（3-1）,采用數(shù)值積分法，從而得到（3-5）。對(duì)于（3-5）右端的積分用矩形公式(取左端點(diǎn))，則得到+1f(x,y(x))dx?h-f(x,y(x))iiXi進(jìn)而得到（3-1）的求解公式（3-2）y.=y.+h-f(x.,y.)(i=0,1,2,n-1)(3-6)此公式稱為歐拉（Euler）格式。如果對(duì)式（3-5）右端的積分用梯形公式』x.+1f(x,y(x))dx牝2-(f(x.,y(x.))+f(x.回，y(x)))x.2l+則可以得到初值問(wèn)題(3-1)的梯形求解公式如式(3-7)i+1y=y+—-[.f（x,y）+f（x,y）]（i=0,1,2,n-1）（3-7）i+1i2iii+1i=1式（3-7）是個(gè)隱式公式。可以采取先用歐拉格式求一個(gè)y（xi+「的初步近似值，記作"+1,稱之為預(yù)報(bào)值，然后用預(yù)報(bào)值y替代式（3-7）右端的y，再計(jì)算得到y(tǒng)，稱之為校正i+1i+1i+1值，這樣建立起來(lái)的預(yù)報(bào)一校正方法稱為改進(jìn)歐拉格式y(tǒng),+1=y,+h-f(x.,y.)y=y+h?[/(x,y)+f(x,yJi+1i2iii+1i+1(3-8)2.龍格一庫(kù)塔方法在單步法中，應(yīng)用最廣泛的是龍格一庫(kù)塔（Runge-kutta）法，簡(jiǎn)稱R—K法。下面直接給出一種四階的龍格一庫(kù)塔法的計(jì)算公式（3-9）n=七+6(K1+2K疽2K3+K4)K,=h-/3.’y.)jK=h.f(X;h,y+1k)(3-9)i2.21K=h.f(x+h,y+-K)，i2i22K耳=h.f(x.+h,y+%)它也稱為標(biāo)準(zhǔn)(古典)龍格一庫(kù)塔法。例3-1研究下列微分方程的初值問(wèn)題y'=~^-2y2jl+x2〔y(0)=0解：這是一個(gè)特殊的微分方程，其解的解析式可以給出，為Xy=1+x2應(yīng)用龍格一庫(kù)塔法，取h=0.25,根據(jù)式(3-9)編寫一段程序，由零開(kāi)始自左相右逐步算出y(x)在所有節(jié)點(diǎn)x.上的近似值y.。計(jì)算結(jié)果見(jiàn)表3-1。計(jì)算結(jié)果表明，四階龍格一庫(kù)塔方法的精度是較高的。表3-1xy?y(x)-y2.00.399956994.3e-54.00.235291592.5e-66.00.162161793.7e-78.00.123076839.2e-8實(shí)際上，MATLAB為常微分方程提供了很好的解題指令，使得求解常微分方程變得很容易，并且能將問(wèn)題及解答表現(xiàn)在圖形上。因此，我們可以不用根據(jù)式(3-9)編寫較復(fù)雜的程序，而只需應(yīng)用MATLAB提供的常微分方程解題器來(lái)解決問(wèn)題。下面給出用MATLAB編寫的解題程序。首先編寫描述常微分方程的ODE文件，文件名為’myfun，便于解題器調(diào)用它。functiondy=myfun(x,y)dy=zeros(1,1);dy=1/(1+x八2)-2火y八2;再編寫利用解題器指令求解y的程序。clearx0=0;fori=1:4xm=2火i;y0=0;[x,y]=ode45('myfun',[x0xm],[y0]);formatlongy(length(y))endplot(x,y,'-')運(yùn)行上述程序，在得到幾個(gè)點(diǎn)的函數(shù)值的同時(shí)，也得到函數(shù)y的曲線，如圖3-1所示。x圖3-1根據(jù)運(yùn)算結(jié)果畫出y的曲線二、簡(jiǎn)單電力系統(tǒng)的暫態(tài)穩(wěn)定性(一)物理過(guò)程分析某簡(jiǎn)單電力系統(tǒng)如圖3-2(a)所示，正常運(yùn)行時(shí)發(fā)電機(jī)經(jīng)過(guò)變壓器和雙回線路向無(wú)限大系統(tǒng)供電。發(fā)電機(jī)用電勢(shì)廳作為其等值電勢(shì)，則電勢(shì)E'與無(wú)限大系統(tǒng)間的電抗為X=X'+X+—L+x(3-10)IdT12T2這時(shí)發(fā)電機(jī)發(fā)出的電磁功率可表示為P=E°sin8=Psin8(3-11)I—IMI如果突然在一回輸電線路始端發(fā)生不對(duì)稱短路，如圖3-2(b)所示。故障期間發(fā)電機(jī)電勢(shì)E'與無(wú)限大系統(tǒng)之間的聯(lián)系電抗為(—'+X)(、+—)TOC\o"1-5"\h\z—、dT12T2x=(x+x)+(―^+x)+(3-12)IIdT12T2XA在故障情況下發(fā)電機(jī)輸出的電磁功率為P=E°sin8=Psin8(3-13)IIxIIMII在短路故障發(fā)生之后，線路繼電保護(hù)裝置將迅速斷開(kāi)故障線路兩端的斷路器，如圖3-2(c)所示。此時(shí)發(fā)電機(jī)電勢(shì)E'與無(wú)限大系統(tǒng)間的聯(lián)系電抗為

發(fā)電機(jī)輸出的功率為Pxsin8=Psin8IIIIII(3-14)(3-15)(c)jXT2UUjjxE'*jxd初jx'jx發(fā)電機(jī)輸出的功率為Pxsin8=Psin8IIIIII(3-14)(3-15)(c)jXT2UUjjxE'*jxd初jx'jxjx如果正常時(shí)發(fā)電機(jī)向無(wú)限大系統(tǒng)輸送的有功功率為P0，則原動(dòng)機(jī)輸出的機(jī)械功率Pt等于P0。假定不計(jì)故障后幾秒種之內(nèi)調(diào)速器的作用，即認(rèn)為機(jī)械功率始終保持Po。因此，可以得到此簡(jiǎn)單電力系統(tǒng)正常運(yùn)行、故障期間及故障切除后的功率特性曲線如圖3-3所示。圖3-3簡(jiǎn)單系統(tǒng)正常運(yùn)行、故障期間及故障切除后的功率特性曲線圖3-3簡(jiǎn)單系統(tǒng)正常運(yùn)行、故障期間及故障切除后的功率特性曲線PT對(duì)于上述簡(jiǎn)單電力系統(tǒng)，我們可以根據(jù)等面積定則求得極限切除角。但是，實(shí)際工作需要知道在多少時(shí)間之內(nèi)切除故障線路，也就是要知道與極限切除角對(duì)應(yīng)的極限切除時(shí)間。要解決這個(gè)問(wèn)題，必須求解發(fā)電機(jī)的轉(zhuǎn)子運(yùn)動(dòng)方程。（二）求解發(fā)電機(jī)的轉(zhuǎn)子運(yùn)動(dòng)方程

求解發(fā)電機(jī)轉(zhuǎn)子運(yùn)動(dòng)方程可以得出0-t和3-t的關(guān)系曲線。其中0-t曲線一般稱為搖擺曲線。在上述簡(jiǎn)單電力系統(tǒng)中故障期間的轉(zhuǎn)子運(yùn)動(dòng)方程為=(①一1)①(3-16)dt1d?1/=—(P一Psino)j式中，o——功率角，其單位為弧度；少一一轉(zhuǎn)子角速度，標(biāo)幺值；氣一一轉(zhuǎn)子的同步角速度，即氣=2#=314.16，其單位為弧度秒；T——發(fā)電機(jī)的慣性時(shí)間常數(shù)，其單位為秒;P、P分別為機(jī)械和電磁功率，標(biāo)幺值。(3-16)這是兩個(gè)一階的非線性常微分方程，它的起始條件是已知的，即t=t=0；?=①=1.0；0=0=sin-1-t-000PIM故障切除后，由于系統(tǒng)參數(shù)改變，以致發(fā)電機(jī)功率特性發(fā)生變化，必須開(kāi)始求解另一組微分方程：~~=(?-1)?"[1(3-17)d?1/dD?e.M=廠(Pt"iIIMSin0)IJ式中變量含義同前述，其中PIIIM也為標(biāo)幺值。這組方程的起始條件為t=t；0=0；?=?其中t為給定的切除時(shí)間；0、?為與t時(shí)刻對(duì)應(yīng)的0和?，它們可由故障期間的0-t和CCCC3-t的關(guān)系曲線求得(0和?都是不突變的)。一般來(lái)說(shuō)，在計(jì)算故障發(fā)生后幾秒種的過(guò)程中，如果0始終不超過(guò)180°，而且振蕩幅值越來(lái)越小，則系統(tǒng)是暫態(tài)穩(wěn)定的。當(dāng)發(fā)電機(jī)與無(wú)限大系統(tǒng)之間發(fā)生振蕩或失去同步時(shí)，在發(fā)電機(jī)的轉(zhuǎn)子回路中，特別是阻尼繞組中將有感應(yīng)電流而形成阻尼轉(zhuǎn)矩(也稱為異步轉(zhuǎn)矩)。當(dāng)作微小振蕩時(shí)，阻尼功率可表達(dá)為：P=DA?=D(&-1)(3-18)D式中，D稱為阻尼功率系數(shù)；A?為轉(zhuǎn)子角速度的偏移量，標(biāo)幺值；?為轉(zhuǎn)子角速度，標(biāo)幺值。阻尼功率系數(shù)D除了與發(fā)電機(jī)的參數(shù)有關(guān)外，還和原始功角、A0的振蕩頻率有關(guān)。在一般情況下它是正數(shù)。在原始功角較小，或者定子回路中有串聯(lián)電容使定子回路總電阻相對(duì)于總電抗較大時(shí)，D可能為負(fù)數(shù)。如果考慮阻尼功率的影響，則故障后的轉(zhuǎn)子運(yùn)動(dòng)方程又可表達(dá)為

-—=（①一1）8觥11（3-19）芍=尸『Di%、靜］lJ電力系統(tǒng)暫態(tài)穩(wěn)定計(jì)算包括兩類問(wèn)題，一類是應(yīng)用數(shù)值計(jì)算法得出故障期間的曲線后，根據(jù)曲線找到與極限切除角對(duì)應(yīng)的極限切除時(shí)間，此時(shí)只需要求解微分方程（3-16）；另一類是已知故障切除時(shí)間，需要求出搖擺曲線來(lái)判斷系統(tǒng)的穩(wěn)定性，此時(shí)需要分段分別求解微分方程（3-16）和（3-17）。如果考慮阻尼轉(zhuǎn)矩的影響，則此時(shí)需要分段分別求解微分方程（3-16）和（3-19）。三、例題例3-2某簡(jiǎn)單電力系統(tǒng)如圖3-4所示，取基準(zhǔn)值七=220MVA，U；=209KV。換算后的參數(shù)已經(jīng)標(biāo)在圖中，其中一回線的電抗七=0.486,=8.18秒。設(shè)電力線路某一回的始端發(fā)生兩相接地短路。假定E'=常數(shù)。（1）計(jì)算保持暫態(tài)穩(wěn)定而要求的極限切除角。（2）計(jì)算極限切除時(shí)間，并且作出在0.15秒切除故障時(shí)的＜5-t曲線。圖3-4圖3-4某簡(jiǎn)單電力系統(tǒng)的接線圖解：計(jì)算系統(tǒng)正常運(yùn)行方式，決定礦和80。由3-3（a）的正序網(wǎng)絡(luò)可得，此時(shí)系統(tǒng)的總電抗為氣=0.295+0.138+0.243+0.122=0.798發(fā)電機(jī)的暫態(tài)電勢(shì)為:礦+*798)2+(*理)2=1.410.7981.01.00.79880=tgt1.0+0.2x0.798~34,53（2）故障后的功率特性又由3-3(b)的負(fù)序、零序網(wǎng)絡(luò)可得故障點(diǎn)的負(fù)序、零序等值電抗為(0.432+0.138)-(0.243+0.122)x==02222￡(0.432+0.138)+(0.243+0.122)'

0.138(0.972+0.122)x==01230￡0.138+(0.972+0.122),所以在正序網(wǎng)絡(luò)故障點(diǎn)上的附加電抗為:x=些空1癸=0.079a0.222+0.123于是故障時(shí)等值電路如圖3-3（c）所示，則x=0.433+0.365+°433+°365=2.80ii0.079因此，故障期間發(fā)電機(jī)的最大功率為:PIIPIIMxIIER1.羅=Q504（3）故障切除后的功率特性故障切除后的等值電路如圖3-3（d）所示xIII=0.295+0.138+0.486+0.122=1.041此時(shí)最大功率為P=空=14頃.0=135IIIMxiii1.0415h=1800-sin-1巖=132.20U=1.0E'=七41心95jQ438j熾43j心25=34.530fmP_100⑴p-L0(a)Q-0.2|_jM3j4^jf(2)j,0/38.j0.243.x4jj0.122Ihf(0)(b)U=1.0E'=%41j心5j心8j蜘3j也22(c)j0.079U=1.0E'=X41j35j蜘8j26j皿2(d)圖3-5例題7-12的等值電路(c)（a）正常運(yùn)行等值電路；（b）負(fù)序和零序等值電路；（c）故障時(shí)等值電路；（d）故障切除后等值電路（4）計(jì)算極限切除角cos5=PT00一8「+二Mcos>*cos80c'"PlW-PlM兀1.0X面(132.2-34.53)+1.35cos132.2o-0.504cos34.53o一1.35-0.504=0.4585=62.740cm找出極限切除時(shí)間，cm根據(jù)(3-16),首先計(jì)算初值53453XK=0.6027,W=1.001800令y(1)=5,y(2)=o。編寫描述故障期間轉(zhuǎn)子運(yùn)動(dòng)方程的ODE文件，文件名為’myequ‘。functiondy=myequ(t,y)dy=zeros(2,1);f=50;w1=2火pi*f;dy(1)=(y(2)-1)*w1;dy(2)=(1/8.18)*(1.0-0.504火sin(y(1)));再編寫利用解題器指令求解y的程序。cleart0=0;tm=0.25;d0=(34.53/180)*pi;w0=1;[T,Y]=ode45('myequ',[t0tm],[d0w0]);plot(T,(Y(:,1)/pi)*180,'-',0.194,62.76,'*')text(0.194,60,'\delta_{cmax}=62.76\circ','FontSize',10)text(0.194,56,'t{cmax}=0.194s','FontSize',10)圖3-6例題7-12的0-t曲線圖3-6給出短路發(fā)生后0秒到0.25秒期間的0-t計(jì)算曲線，根據(jù)最大切除角5cm(=62.740)找到極限切除時(shí)間'為0.194秒。由圖3-6可見(jiàn)，如果故障切除時(shí)間大于0.194秒，則發(fā)電機(jī)的功角將不斷地增大，最終失去暫態(tài)穩(wěn)定。在極限切除時(shí)間之前切除故障，發(fā)電機(jī)的搖擺曲線的狀況將在下面作計(jì)算、分析。不考慮阻尼轉(zhuǎn)矩的影響，當(dāng)故障切除時(shí)間為0.15秒時(shí)通過(guò)計(jì)算得出S-t曲線首先編寫描述故障期間轉(zhuǎn)子運(yùn)動(dòng)方程的ODE文件，文件名為"myfun01”。functiondy=myfun01(t,y)f=50;w1=2火pi*f;TJ=8.18;Pt=1.0;P2m=0.504;dy=zeros(2,1);dy(1)=(y(2)-1)*w1;dy(2)=(1/TJ)*(Pt-P2m*sin(y(1)));再編寫描述故障切除后轉(zhuǎn)子運(yùn)動(dòng)方程的ODE文件，文件名為”myfun02”。functiondy=myfun02(t,y)f=50;w1=2火pi*f;TJ=8.18;Pt=1.0;P3m=1.35;dy=zeros(2,1);dy(1)=(y(2)-1)*w1;dy(2)=(1/TJ)*(Pt-P3m*sin(y(1)));編寫利用解題器指令求解y的小程序。cleart0=0;tc=0.15;tm=2.0;d0=(34.53/180)*pi;w0=1.0;[T1,Y1]=ode45('myfun01',[t0tc],[d0w0]);dc=Y1(length(Y1),1);wc=Y1(length(Y1),2);[T2,Y2]=ode45('myfun02',[tctm],[dcwc]);plot(T1,(Y1(:,1)/pi)*180,'-',T2,(Y2(:,1)/pi)*180,'-',tc,(dc/pi)*180,'*')text(0.28,50,'\it{t}_{c}=0.15s','FontSize',8)text(0.28,43,'\it{\delta}_{c}=51.71\circ','FontSize',8)xlabel('\it{t}')ylabel('\it{\delta}')計(jì)算結(jié)果表明，功角8沿著故障切除后的功角特性曲線根據(jù)等面積定則作等幅振蕩，如圖3-7所示。實(shí)際上，由于阻尼轉(zhuǎn)矩的影響，振蕩的幅度是逐漸衰減的，功角8最終運(yùn)行在8廣47.8。。因此，發(fā)電機(jī)能夠保持暫態(tài)穩(wěn)定。圖3-7不考慮阻尼轉(zhuǎn)矩影響，當(dāng)0.15秒切除故障時(shí)發(fā)電機(jī)的5-1曲線考慮阻尼轉(zhuǎn)矩的影