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

1、實(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)題，從而得出3 -t和3 -t的關(guān)系曲線(xiàn)。每臺(tái)發(fā)電機(jī)的轉(zhuǎn)子運(yùn)動(dòng)方程是兩個(gè)一階非線(xiàn)性的常微分方程。因此， 首先介紹常微分方程的初值問(wèn)題的數(shù)值解法。-、常微分方程的初值問(wèn)題(一)問(wèn)題及求解公式的構(gòu)造方法我們討論形如式(3-1 )的一階微分方程的初值問(wèn)題(3-1 )y (x) f(x, y), a x b y(xo) yo設(shè)初值問(wèn)題(3-1 )的解為y(x)，為了求其數(shù)值解而采取離散化方法，在求解區(qū)間a,b上取一組節(jié)點(diǎn)a xo x1 xi xi 1xn b稱(chēng)hixi 1 xi ( i 0,1, ,n 1 )為步長(zhǎng)。

2、在等步長(zhǎng)的情況下，步長(zhǎng)為b ahn用yi表示在節(jié)點(diǎn)Xi處解的準(zhǔn)確值y(Xi)的近似值。設(shè)法構(gòu)造序列yi所滿(mǎn)足的一個(gè)方程(稱(chēng)為差分方程)yi 1yi h (Xi,yi,h)(3-2 )作為求解公式，這是一個(gè)遞推公式，從(Xo, yo)出發(fā)，采用步進(jìn)方式，自左相右逐步算出y(x)在所有節(jié)點(diǎn)x上的近似值yi ( i 1,2, n )。在公式(3-2 )中，為求yi 1只用到前面一步的值 yi ,這種方法稱(chēng)為單步法。在公式(3-2 )中的yi 1由yi明顯表示出，稱(chēng)為顯式公式。而形如(3-3)yi 1yi h (Xi, yi, yi 1, h)(3-3 )的公式稱(chēng)為隱式公式，因?yàn)槠溆叶酥羞€包括y 1。

3、如果由公式求y 1時(shí)，不止用到前一個(gè)節(jié)點(diǎn)的值，則稱(chēng)為多步法。由式(3-1 )可得dy = f (x, y)dx(3-4 )兩邊在x , Xi i 上積分，得xi 1y(Xi 1) y(Xi)f (x, y(x)dx(3-5 )xi由此可以看出，如果想構(gòu)造求解公式，就要對(duì)右端的積分項(xiàng)作某種數(shù)值處理。這種求解公式的構(gòu)造方法叫做數(shù)值積分法。(二)一般的初值問(wèn)題的解法1.歐拉法和改進(jìn)歐拉法對(duì)于初值問(wèn)題(3-1 )，采用數(shù)值積分法，從而得到(3-5 )。對(duì)于(3-5 )右端的積分用矩形公式(取左端點(diǎn))，則得到xi 1f (x, y(x)dx h f(Xi,y(Xi)Xi進(jìn)而得到(3-1 )的求解公式(3

4、-2 )yi 1yi h f (Xi, yi)( i =0, 1, 2, n-1 )( 3-6)此公式稱(chēng)為歐拉(Euler )格式。如果對(duì)式(3-5 )右端的積分用梯形公式為1hf (x, y(x)dx(f(Xi ,y(Xi) f(x 1,y(Xi 1)Xi2則可以得到初值問(wèn)題(3-1 )的梯形求解公式如式(3-7 )hyi 1yif (Xi, yi) f (x 1, yi 1)(i =0, 1, 2, n-1 )(3-7 )2式(3-7 )是個(gè)隱式公式?？梢圆扇∠扔脷W拉格式求一個(gè)y(xi 1)的初步近似值，記作 yi 1 ,稱(chēng)之為預(yù)報(bào)值，然后用預(yù)報(bào)值 yi 1替代式(3-7 )右端的yi 1

5、，再計(jì)算得到y(tǒng)i 1，稱(chēng)之為校正 值，這樣建立起來(lái)的預(yù)報(bào)-校正方法稱(chēng)為改進(jìn)歐拉格式y(tǒng) 1 yi h f (Xi, yi)h-(3-8 )yi 1 yif (Xi, yi) f (x 1, yi 1)22.龍格一庫(kù)塔方法在單步法中，應(yīng)用最廣泛的是龍格庫(kù)塔（Runge-kutta ）法，簡(jiǎn)稱(chēng)R K法。下面直接給出一種四階的龍格一庫(kù)塔法的計(jì)算公式（3-9）yi iKi1K4)yi 丄(Ki 2K2 2K3 6(x'f,yi,yiyi加） 土心）K3）(3-9)h f(Xi,yJ它也稱(chēng)為標(biāo)準(zhǔn)（古典）龍格-庫(kù)塔法。例3-1 研究下列微分方程的初值問(wèn)題2y1y 21 Xy(0) 0解:這是一個(gè)特殊

6、的微分方程，其解的解析式可以給出，為應(yīng)用龍格庫(kù)塔法，取 h=,根據(jù)式（3-9）編寫(xiě)一段程序，由零開(kāi)始自左相右逐步算出y（x）在所有節(jié)點(diǎn)Xi上的近似值yi。計(jì)算結(jié)果見(jiàn)表3-1。計(jì)算結(jié)果表明，四階龍格庫(kù)塔方法的精 度是較高的。表3-1Xnyny(Xn ) yn0.0.0.0.實(shí)際上，MATLAB為常微分方程提供了很好的解題指令，使得求解常微分方程變得很容易，并且能將問(wèn)題及解答表現(xiàn)在圖形上。因此，我們可以不用根據(jù)式（3-9 ）編寫(xiě)較復(fù)雜的程序，而只需應(yīng)用MATLAB供的常微分方程解題器來(lái)解決問(wèn)題。下面給出用MATLAB編寫(xiě)的解題程序。首先編寫(xiě)描述常微分方程的ODE文件，文件名為/ myfun便于解題

7、器調(diào)用它。fun cti ondy = myfun( x,y)dy = zeros(1,1);dy=1/(1+xA2)-2*yA2;再編寫(xiě)利用解題器指令求解y的程序。clearx0=0;for i=1:4xm=2*i;y0=0;x,y = ode45('myfun' ,x0 xm,yO);format longy(le ngth(y)endplot(x,y,'-')運(yùn)行上述程序，在得到幾個(gè)點(diǎn)的函數(shù)值的同時(shí)，也得到函數(shù)y的曲線(xiàn)，如圖3-1所示。x圖3-1 根據(jù)運(yùn)算結(jié)果畫(huà)出y的曲線(xiàn)二 簡(jiǎn)單電力系統(tǒng)的暫態(tài)穩(wěn)定性(一)物理過(guò)程分析某簡(jiǎn)單電力系統(tǒng)如圖3-2(a)所示，正常

8、運(yùn)行時(shí)發(fā)電機(jī)經(jīng)過(guò)變壓器和雙回線(xiàn)路向無(wú)限大系統(tǒng)供電。發(fā)電機(jī)用電勢(shì) E作為其等值電勢(shì)，則電勢(shì) E與無(wú)限大系統(tǒng)間的電抗為XdXT1Xl2XT2(3-10)這時(shí)發(fā)電機(jī)發(fā)出的電磁功率可表示為EU .sinxPm sin如果突然在一回輸電線(xiàn)路始端發(fā)生不對(duì)稱(chēng)短路，如圖(3-11)3-2(b)所示。故障期間發(fā)電機(jī)電勢(shì)E與無(wú)限大系統(tǒng)之間的聯(lián)系電抗為XlXlX(XdXti ) (Xt2)2(XdXti )(XT2)(3-12)2在故障情況下發(fā)電機(jī)輸出的電磁功率為EUsinP m sin(3-13)在短路故障發(fā)生之后，線(xiàn)路繼電保護(hù)裝置將迅速斷開(kāi)故障線(xiàn)路兩端的斷路器，如圖3-2(c)所示。此時(shí)發(fā)電機(jī)電勢(shì)與無(wú)限大系統(tǒng)間

9、的聯(lián)系電抗為XdXt1XlXT2(3-14)發(fā)電機(jī)輸出的功率為E U sin Pxm sin(3-15)T2U= cLUGT1E jXdjXT1-jXL jXT2 UjXLjx(b)EjXdjXT1jXLjXT2(c)圖3-2簡(jiǎn)單電力系統(tǒng)及其等值電路(a)正常運(yùn)行方式及其等值電路；(b)故障情況及其等值電路；(c)故障切除后及其等值電路如果正常時(shí)發(fā)電機(jī)向無(wú)限大系統(tǒng)輸送的有功功率為P。，則原動(dòng)機(jī)輸出的機(jī)械功率Pt等Po。因此，可于Po。假定不計(jì)故障后幾秒種之內(nèi)調(diào)速器的作用，即認(rèn)為機(jī)械功率始終保持以得到此簡(jiǎn)單電力系統(tǒng)正常運(yùn)行、故障期間及故障切除后的功率特性曲線(xiàn)如圖3-3所示。PtPo圖3-3簡(jiǎn)單系

10、統(tǒng)正常運(yùn)行、故障期間及故障切除后的功率特性曲線(xiàn)對(duì)于上述簡(jiǎn)單電力系統(tǒng)，我們可以根據(jù)等面積定則求得極限切除角。 但是，實(shí)際工作需要知道在多少時(shí)間之內(nèi)切除故障線(xiàn)路， 也就是要知道與極限切除角對(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)方程可以得出3-t和3 -t的關(guān)系曲線(xiàn)。其中3 -t曲線(xiàn)一般稱(chēng)為搖擺曲線(xiàn)。在上述簡(jiǎn)單電力系統(tǒng)中故障期間的轉(zhuǎn)子運(yùn)動(dòng)方程為1)式中,d 1打 tj(PtPMSin功率角，其單位為弧度;(3-16 )轉(zhuǎn)子角速度，標(biāo)幺值;轉(zhuǎn)子的同步角速度，即1=2 f =，其單位為弧度/秒；Tj 發(fā)電機(jī)的慣性時(shí)間常數(shù)，其單位

11、為秒；Pt、Pm 分別為機(jī)械和電磁功率，標(biāo)幺值。這是兩個(gè)一階的非線(xiàn)性常微分方程，它的起始條件是已知的，即t =to=O;0=；= sin故障切除后，由于系統(tǒng)參數(shù)改變， 以致發(fā)電機(jī)功率特性發(fā)生變化，必須開(kāi)始求解另一組 微分方程:-(1) 1d1(3-17 )(Ptp m sin )dtTj式中變量含義冋前述，其中P M也為標(biāo)幺值。這組方程的起始條件為t =tc ；= c ；= c其中tc為給定的切除時(shí)間；c、c為與tc時(shí)刻對(duì)應(yīng)的和，它們可由故障期間的3 -t和3 -t的關(guān)系曲線(xiàn)求得（ 和 都是不突變的）。一般來(lái)說(shuō)，在計(jì)算故障發(fā)生后幾秒種的過(guò)程 中，如果3始終不超過(guò) 180 0 ,而且振蕩幅值越來(lái)

12、越小，則系統(tǒng)是暫態(tài)穩(wěn)定的。當(dāng)發(fā)電機(jī)與無(wú)限大系統(tǒng)之間發(fā)生振蕩或失去同步時(shí)，在發(fā)電機(jī)的轉(zhuǎn)子回路中， 特別是阻尼繞組中將有感應(yīng)電流而形成阻尼轉(zhuǎn)矩（也稱(chēng)為異步轉(zhuǎn)矩）。當(dāng)作微小振蕩時(shí)，阻尼功率可表達(dá)為：Pd = D =D（1）（3-18）式中，D稱(chēng)為阻尼功率系數(shù)；為轉(zhuǎn)子角速度的偏移量，標(biāo)幺值；為轉(zhuǎn)子角速度，標(biāo)幺值。阻尼功率系數(shù) d除了與發(fā)電機(jī)的參數(shù)有關(guān)外，還和原始功角、的振蕩頻率有關(guān)。在般情況下它是正數(shù)。 在原始功角較小，或者定子回路中有串聯(lián)電容使定子回路總電阻相對(duì)于總電抗較大時(shí)，d可能為負(fù)數(shù)。如果考慮阻尼功率的影響，則故障后的轉(zhuǎn)子運(yùn)動(dòng)方程又可 表達(dá)為dt d1)(3-19 )dtTJPtD(1) P

13、m sin電力系統(tǒng)暫態(tài)穩(wěn)定計(jì)算包括兩類(lèi)問(wèn)題，一類(lèi)是應(yīng)用數(shù)值計(jì)算法得出故障期間的曲線(xiàn)后,根據(jù)曲線(xiàn)找到與極限切除角對(duì)應(yīng)的極限切除時(shí)間，此時(shí)只需要求解微分方程（3-16 ）;另類(lèi)是已知故障切除時(shí)間，需要求出搖擺曲線(xiàn)來(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)值Sb=220MVA Ub=209KV。換算后的參數(shù)已經(jīng)標(biāo)在圖中，其中一回線(xiàn)的電抗Xl=, Tj =秒。設(shè)電力線(xiàn)路某一回的始端發(fā)生兩相接地短路。假定E =常數(shù)。(1)計(jì)算保

14、持暫態(tài)穩(wěn)定而要求的極限切除角。(2)計(jì)算極限切除時(shí)間,并且作出在秒切除故障時(shí)的3 -t曲線(xiàn)。Xd =0.295Xl =0.486Xlo4XlXt2=0.122Po=1.OQo=O.2X2 =0.432Xti =0.138圖3-4某簡(jiǎn)單電力系統(tǒng)的接線(xiàn)圖解：計(jì)算系統(tǒng)正常運(yùn)行方式，決定E和0。由3-3(a)的正序網(wǎng)絡(luò)可得，此時(shí)系統(tǒng)的總電抗為x =+=發(fā)電機(jī)的暫態(tài)電勢(shì)為:0.2 0.798)21.01(1.0 0.798)21.00.798=01.0 0.2 0.798(2)故障后的功率特性又由3-3(b)的負(fù)序、零序網(wǎng)絡(luò)可得故障點(diǎn)的負(fù)序、零序等值電抗為=(0.432 0.138) (0.243 0.

15、122)X2 = (0.432 0.138)(0.243 0.122)0.138(0.972 0.122) _Xo =0.138 (0.972 0.122)所以在正序網(wǎng)絡(luò)故障點(diǎn)上的附加電抗為:x于是故障時(shí)等值電路如圖°222 °1230.0790.222 0.1233-3(c)所示，則0.433 0.365x0.4330.3652.800.079因此，故障期間發(fā)電機(jī)的最大功率為:EUx1.41 1.02.80.504(3 )故障切除后的功率特性故障切除后的等值電路如圖3-3(d)所示01800.2950.138 0.486 0.1221.041此時(shí)最大功率為EU竺衛(wèi)1.35

16、1.0411800 sin1 型1.351322°j 0.138 j0.243 4 j 0.122I'(c)(a)正常運(yùn)行等值電路;(b)負(fù)序和零序等值電路；(C)故障時(shí)等值電路;(d)故障切除后等值電路計(jì)算極限切除角f(0)(b)U=1.01.41 j 0.295 j 0.138 j 0.243 j 0.122n_/"vv%j 0.079U=1.01.41j 0.295 j 0.138 j 0.486 j 0.122:-rr-r(d)圖3-5 例題7-12的等值電路COS cmPT( 0 h) P M COS h P M COS 0Pm P m00.504 cos

17、34.530° 面(1322 3453)35cos132.21.35 0.504cm62.740(5)找出極限切除時(shí)間tcm根據(jù)(3-16)，首先計(jì)算初值34.530.6027,1.0令y(1)=, y(2)=。編寫(xiě)描述故障期間轉(zhuǎn)子運(yùn)動(dòng)方程的ODE文件，文件名為/myequ'。fun cti ondy = myequ(t,y)dy = zeros(2,1);f=50;w1=2*pi*f;dy(1) = (y(2)-1)*w1;dy(2) = (1/*再編寫(xiě)利用解題器指令求解y的程序。cleart0=0;tm=;d0=180)*pi;w0=1;T,Y = ode45( '

18、;myequ' ,t0 tm,dO wO);plot(T,(Y(:,1)/pi)*180,'-','*')text,60, 'delta_cmax=circ' , 'FontSize' ,10)text,56, 't_cmax=', 'FontSize' ,10)圖3-6 例題7-12的S -t曲線(xiàn)圖3-6給出短路發(fā)生后0秒到秒期間的S -t計(jì)算曲線(xiàn)，根據(jù)最大切除角cm ( 62.740)找到極限切除時(shí)間tcm為秒。由圖3-6可見(jiàn)，如果故障切除時(shí)間大于秒，則發(fā)電機(jī)的功角將 不斷地增大，最終失

19、去暫態(tài)穩(wěn)定。在極限切除時(shí)間之前切除故障，發(fā)電機(jī)的搖擺曲線(xiàn)的狀況將在下面作計(jì)算、分析。(6 )不考慮阻尼轉(zhuǎn)矩的影響，當(dāng)故障切除時(shí)間為秒時(shí)通過(guò)計(jì)算得出S-t曲線(xiàn)首先編寫(xiě)描述故障期間轉(zhuǎn)子運(yùn)動(dòng)方程的ODE文件，文件名為”myfun01functiondy = myfun01(t,y)文件名為” myfun02 ”。'-' ,tc,(dc/pi)*180, '*' )f=50; w1=2*pi*f;TJ=; Pt=; P2m=;dy = zeros(2,1);dy(1) = (y(2)-1)*w1;dy(2) = (1/TJ)*(Pt-P2m*sin(y(1);再編寫(xiě)描述

20、故障切除后轉(zhuǎn)子運(yùn)動(dòng)方程的ODE文件,function dy = myfun02(t,y)f=50; w1=2*pi*f;TJ=; Pt=; P3m=;dy = zeros(2,1);dy(1) = (y(2)-1)*w1;dy(2) = (1/TJ)*(Pt-P3m*sin(y(1);編寫(xiě)利用解題器指令求解y的小程序。cleart0=0; tc=; tm=;d0=180)*pi; w0=;T1,Y1 = ode45( 'myfun01' ,t0 tc,d0 w0);dc=Y1(length(Y1),1);wc=Y1(length(Y1),2);T2,Y2 = ode45( &#

21、39;myfun02' ,tc tm,dc wc); plot(T1,(Y1(:,1)/pi)*180, '-' ,T2,(Y2(:,1)/pi)*180, text,50, 'itt_c=' , 'FontSize' ,8)text,43, 'itdelta_c=circ' , 'FontSize' ,8) xlabel( 'itt' )ylabel( 'itdelta' )計(jì)算結(jié)果表明， 功角 沿著故障切除后的功角特性曲線(xiàn)根據(jù)等面積定則作等幅振蕩，最終運(yùn)行在圖3-7所示。實(shí)際上，由于阻尼轉(zhuǎn)矩的影響，振蕩的幅度是逐漸衰減的，功角k=o。因此，發(fā)電機(jī)能夠保持暫態(tài)穩(wěn)定。圖3-7不考慮阻尼轉(zhuǎn)矩影響，當(dāng)秒切除故障時(shí)發(fā)電機(jī)的t曲線(xiàn)(7)考慮阻尼轉(zhuǎn)矩的影響，當(dāng)故障切除時(shí)間為秒時(shí)通過(guò)計(jì)算得出3-t曲線(xiàn)描述故障期間轉(zhuǎn)子運(yùn)動(dòng)方程的ODE文件與(6)相同，文件名也為” myfunOI ”。重新編寫(xiě)描述故障切除后