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

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

2、步長的情況下，步長為hban用 yi 表示在節(jié)點(diǎn) xi 處解的準(zhǔn)確值 y( xi ) 的近似值。設(shè)法構(gòu)造序列 yi 所滿足的一個(gè)方程（稱為差分方程）yi 1yi h (xi , yi , h)（3-2 ）作為求解公式，這是一個(gè)遞推公式，從（x0 ， y0）出發(fā)，采用步進(jìn)方式，自左相右逐步算出 y( x) 在所有節(jié)點(diǎn) xi 上的近似值 yi （ i1,2, , n ）。在公式（3-2 ）中，為求 yi 1 只用到前面一步的值yi，這種方法稱為單步法。 在公式（ 3-2 ）中的 yi 1 由 yi 明顯表示出，稱為顯式公式。而形如（3-3 ）yi 1yi h (xi , yi , yi 1 , h

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

4、 ，則得到xi 1f ( x, y( x) dxh f (xi , y( xi)xi進(jìn)而得到（ 3-1 ）的求解公式（ 3-2 ）yi 1 yih f ( xi , yi )（ i=0， 1，2， n-1 ）（3-6）此公式稱為歐拉（ Euler ）格式。如果對(duì)式（ 3-5 ）右端的積分用梯形公式xi 1hf (x, y(x)dx( f ( xi , y(xi ) f ( xi1 , y(xi 1 )xi2則可以得到初值問題（ 3-1 ）的梯形求解公式如式（3-7 ）hf (xi , yi ) f (xi 1 , yi 1 )（ i =0， 1， 2，n-1 ） （3-7 ）y i 1y i2

5、式（ 3-7 ）是個(gè)隱式公式?？梢圆扇∠扔脷W拉格式求一個(gè)y( xi 1 ) 的初步近似值，記作yi 1 ，稱之為預(yù)報(bào)值，然后用預(yù)報(bào)值yi 1 替代式（ 3-7 ）右端的 yi 1，再計(jì)算得到 yi 1 ，稱之為校正值，這樣建立起來的預(yù)報(bào)校正方法稱為改進(jìn)歐拉格式y(tǒng) i1y ih f ( xi , y i )（3-8 ）y iy ih)1f (xi , yif (xi 1 , y i 1 )22 龍格庫塔方法在單步法中，應(yīng)用最廣泛的是龍格庫塔（Runge-kutta）法，簡稱R K 法。下面直接給出一種四階的龍格庫塔法的計(jì)算公式（3-9 ）yi 1yi1 (K12K2 2K 3K 4 )K 1h6

6、f (xi , yi )K 2h f (xih1（3-9 ）, yiK1 )22K 3h f (xih , yi1K2)K 4h f (xi22h, yiK 3 )它也稱為標(biāo)準(zhǔn)（古典）龍格庫塔法。例 3-1研究下列微分方程的初值問題y12 2 y 21xy (0)0解：這是一個(gè)特殊的微分方程，其解的解析式可以給出，為yx1 x 2應(yīng)用龍格庫塔法，取h =0.25 ，根據(jù)式（ 3-9 ）編寫一段程序，由零開始自左相右逐步算出y( x) 在所有節(jié)點(diǎn) xi 上的近似值 yi 。計(jì)算結(jié)果見表3-1 。計(jì)算結(jié)果表明， 四階龍格庫塔方法的精度是較高的。表 3-1xnyny( xn ) yn2.00.399

7、956994.3e-54.00.235291592.5e-66.00.162161793.7e-78.00.123076839.2e-8實(shí)際上， MATLAB為常微分方程提供了很好的解題指令，使得求解常微分方程變得很容易，并且能將問題及解答表現(xiàn)在圖形上。因此，我們可以不用根據(jù)式（3-9 ）編寫較復(fù)雜的程序，而只需應(yīng)用MATLAB提供的常微分方程解題器來解決問題。下面給出用MATLAB編寫的解題程序。首先編寫描述常微分方程的ODE文件，文件名為myfun，便于解題器調(diào)用它。functiondy = myfun(x,y)dy = zeros(1,1);dy=1/(1+x2)-2*y2;再編寫利用解

8、題器指令求解y的程序。clearx0=0;fori=1:4xm=2*i;y0=0;x,y = ode45(myfun,x0 xm,y0);format longy(length(y)endplot(x,y,-)運(yùn)行上述程序，在得到幾個(gè)點(diǎn)的函數(shù)值的同時(shí)，也得到函數(shù)y 的曲線，如圖3-1 所示。0.50.40.3y0.20.1002468x圖 3-1根據(jù)運(yùn)算結(jié)果畫出y 的曲線二、 簡單電力系統(tǒng)的暫態(tài)穩(wěn)定性（一）物理過程分析某簡單電力系統(tǒng)如圖3-2(a)所示，正常運(yùn)行時(shí)發(fā)電機(jī)經(jīng)過變壓器和雙回線路向無限大系統(tǒng)供電。發(fā)電機(jī)用電勢E 作為其等值電勢，則電勢E 與無限大系統(tǒng)間的電抗為xxdxT 1x LxT

9、 2（ 3-10 ）2這時(shí)發(fā)電機(jī)發(fā)出的電磁功率可表示為PE U sinP M sin（ 3-11 ）x如果突然在一回輸電線路始端發(fā)生不對(duì)稱短路，如圖3-2(b) 所示。故障期間發(fā)電機(jī)電勢 E 與無限大系統(tǒng)之間的聯(lián)系電抗為(xd xT 1 ) ( x L( xdxT 1 )( x LxT 2 )xxT 2 )2（ 3-12 ）2x在故障情況下發(fā)電機(jī)輸出的電磁功率為PE U sinP M sin（ 3-13 ）x在短路故障發(fā)生之后，線路繼電保護(hù)裝置將迅速斷開故障線路兩端的斷路器，如圖3-2(c)所示。此時(shí)發(fā)電機(jī)電勢E 與無限大系統(tǒng)間的聯(lián)系電抗為xxdxT1 xLxT 2（ 3-14 ）發(fā)電機(jī)輸出的

10、功率為PE U sinP M sin（ 3-15 ）xGT1LT2U = cEj xdjx T 1jx Ljx T 2UEjxdjxT 1jx LjxT 2 Ujx Ljx Ljx( a)( b)Ej xdjx T1jxLjx T 2U( c)圖 3-2 簡單電力系統(tǒng)及其等值電路（ a）正常運(yùn)行方式及其等值電路； （b）故障情況及其等值電路； （ c）故障切除后及其等值電路如果正常時(shí)發(fā)電機(jī)向無限大系統(tǒng)輸送的有功功率為P0 ，則原動(dòng)機(jī)輸出的機(jī)械功率PT 等于 P0 。假定不計(jì)故障后幾秒種之調(diào)速器的作用，即認(rèn)為機(jī)械功率始終保持P0 。因此，可以得到此簡單電力系統(tǒng)正常運(yùn)行、故障期間及故障切除后的功率

11、特性曲線如圖3-3 所示。PPPakhPTP0P0kcmh圖 3-3 簡單系統(tǒng)正常運(yùn)行、故障期間及故障切除后的功率特性曲線對(duì)于上述簡單電力系統(tǒng)， 我們可以根據(jù)等面積定則求得極限切除角。但是，實(shí)際工作需要知道在多少時(shí)間之切除故障線路，也就是要知道與極限切除角對(duì)應(yīng)的極限切除時(shí)間。要解決這個(gè)問題，必須求解發(fā)電機(jī)的轉(zhuǎn)子運(yùn)動(dòng)方程。（二）求解發(fā)電機(jī)的轉(zhuǎn)子運(yùn)動(dòng)方程求解發(fā)電機(jī)轉(zhuǎn)子運(yùn)動(dòng)方程可以得出-t和 -t 的關(guān)系曲線。 其中 -t曲線一般稱為搖擺曲線。在上述簡單電力系統(tǒng)中故障期間的轉(zhuǎn)子運(yùn)動(dòng)方程為d1)(1dt（ 3-16 ）d1(PTP M sin )dtT J式中，功率角，其單位為弧度；轉(zhuǎn)子角速度，標(biāo)幺值

12、；1 轉(zhuǎn)子的同步角速度，即1 = 2 f =314.16 ，其單位為弧度 / 秒； TJ 發(fā)電機(jī)的慣性時(shí)間常數(shù)，其單位為秒；PT 、 P M 分別為機(jī)械和電磁功率，標(biāo)幺值。這是兩個(gè)一階的非線性常微分方程，它的起始條件是已知的，即t = t 0 =0；= 0=1.0； =0 = sin 1 PTP M故障切除后， 由于系統(tǒng)參數(shù)改變，以致發(fā)電機(jī)功率特性發(fā)生變化，必須開始求解另一組微分方程：d1)(1（3-17 ）dtd1(PTP M sin )dtT J式中變量含義同前述，其中P M也為標(biāo)幺值。這組方程的起始條件為t = t c ；= c ； = c其中 t c 為給定的切除時(shí)間；c 、c 為與

13、t c 時(shí)刻對(duì)應(yīng)的和 ，它們可由故障期間的 -t和-t 的關(guān)系曲線求得（和都是不突變的） 。一般來說，在計(jì)算故障發(fā)生后幾秒種的過程中，如果始終不超過180o ，而且振蕩幅值越來越小，則系統(tǒng)是暫態(tài)穩(wěn)定的。當(dāng)發(fā)電機(jī)與無限大系統(tǒng)之間發(fā)生振蕩或失去同步時(shí)，在發(fā)電機(jī)的轉(zhuǎn)子回路中， 特別是阻尼繞組中將有感應(yīng)電流而形成阻尼轉(zhuǎn)矩（也稱為異步轉(zhuǎn)矩）。當(dāng)作微小振蕩時(shí)，阻尼功率可表達(dá)為：PD=D=D( 1)（ 3-18 ）式中， D 稱為阻尼功率系數(shù)；為轉(zhuǎn)子角速度的偏移量，標(biāo)幺值；為轉(zhuǎn)子角速度，標(biāo)幺值。阻尼功率系數(shù) D 除了與發(fā)電機(jī)的參數(shù)有關(guān)外，還和原始功角、的振蕩頻率有關(guān)。在一般情況下它是正數(shù)。在原始功角較小，

14、或者定子回路中有串聯(lián)電容使定子回路總電阻相對(duì)于總電抗較大時(shí)， D可能為負(fù)數(shù)。如果考慮阻尼功率的影響，則故障后的轉(zhuǎn)子運(yùn)動(dòng)方程又可表達(dá)為d(1)1dt（3-19 ）d1 PTD ( 1) P M sin dtT J電力系統(tǒng)暫態(tài)穩(wěn)定計(jì)算包括兩類問題，一類是應(yīng)用數(shù)值計(jì)算法得出故障期間的曲線后，根據(jù)曲線找到與極限切除角對(duì)應(yīng)的極限切除時(shí)間，此時(shí)只需要求解微分方程（3-16 ）；另一類是已知故障切除時(shí)間， 需要求出搖擺曲線來判斷系統(tǒng)的穩(wěn)定性， 此時(shí)需要分段分別求解微分方程（ 3-16 ）和（ 3-17 ）。如果考慮阻尼轉(zhuǎn)矩的影響，則此時(shí)需要分段分別求解微分方程（ 3-16 ）和（ 3-19 ）。三、 例題例

15、 3-2 某簡單電力系統(tǒng)如圖3-4 所示， 取基準(zhǔn)值 SB =220MVA，U B =209KV。換算后的參數(shù)已經(jīng)標(biāo)在圖中，其中一回線的電抗xL =0.486 ， TJ =8.18秒。設(shè)電力線路某一回的始端發(fā)生兩相接地短路。假定E =常數(shù)。（ 1）計(jì)算保持暫態(tài)穩(wěn)定而要求的極限切除角。（ 2）計(jì)算極限切除時(shí)間，并且作出在0.15 秒切除故障時(shí)的 -t曲線。GT1LT2U =1.0kx2=0.432xT 1 =0.138xL =0.486xT 2 =0.122xd=0.295x L04x L=1.0=0.2Q0P0圖 3-4某簡單電力系統(tǒng)的接線圖解：計(jì)算系統(tǒng)正常運(yùn)行方式，決定E和 0 。由 3-3

16、(a) 的正序網(wǎng)絡(luò)可得，此時(shí)系統(tǒng)的總電抗為x =0.295+0.138+0.243+0.122=0.798發(fā)電機(jī)的暫態(tài)電勢為：E =(1.00.20.798 ) 2(1.0 0.798 )2 =1.411.01.010.79834.53 o0 =tg0.20.798=1.0（2）故障后的功率特性又由 3-3(b) 的負(fù)序、零序網(wǎng)絡(luò)可得故障點(diǎn)的負(fù)序、零序等值電抗為x2 = (0.4320.138) (0.2430.122) =0.222(0.4320.138)(0.2430.122)x0 = 0.138(0.9720.122)=0.1230.138(0.9720.122)所以在正序網(wǎng)絡(luò)故障點(diǎn)上的

17、附加電抗為：x0.2220.1230.0790.2220.123于是故障時(shí)等值電路如圖3-3(c)所示，則x0.4330.3650.4330.3652.800.079因此，故障期間發(fā)電機(jī)的最大功率為：E U1.41 1.0P M0.504x2.8（3）故障切除后的功率特性故障切除后的等值電路如圖3-3(d) 所示x0.2950.1380.4860.1221.041此時(shí)最大功率為E U1.41 1.0P M1.35x1.041h 1800 sin 1 1.0 132.201.35E 1.41 j0.295j0.138j0.243j0.122U =1.0034.530f(1)P01.0( a)Q0

18、0.2j0.432j0.138j0.243j0.122f(2 )j0.138j 0.2434 j0.122f(0 )( b)U =1.0E1.41 j0.295j 0.138j 0.243j 0.122j 0.079( c)U =1.0E1.41 j0.295j 0.138j 0.486j 0.122(d )圖 3-5 例題 7-12 的等值電路（ a）正常運(yùn)行等值電路； （ b）負(fù)序和零序等值電路； （ c）故障時(shí)等值電路； （ d）故障切除后等值電路（4）計(jì)算極限切除角cosPT( 0h )P McoshP M cos 0cmP MP M1.0(132.2 34.53)1.35 cos13

19、2 .200.504 cos34.53 0= 1801.350.504=0.458cm62.740（5）找出極限切除時(shí)間t cm根據(jù)（ 3-16 ），首先計(jì)算初值034.530.6027 ， 01.0180令 y(1)=，y(2)= 。編寫描述故障期間轉(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=

20、(34.53/180)*pi;w0=1;T,Y = ode45(myequ,t0 tm,d0 w0);plot(T,(Y(:,1)/pi)*180,-,0.194,62.76,*)text(0.194,60,delta_cmax=62.76circ, FontSize,10)text(0.194,56,t_cmax=0.194s, FontSize,10)9080cmax =62.74 ?70tcmax =0.194s605040300.050.10.150.20.250t圖 3-6例題 7-12 的 -t曲線圖 3-6給出短路發(fā)生后0 秒到0.25秒期間的 -t計(jì)算曲線，根據(jù)最大切除角cm（

21、62.740 ）找到極限切除時(shí)間tcm 為0.194秒。由圖3-6可見，如果故障切除時(shí)間大于0.194秒，則發(fā)電機(jī)的功角將不斷地增大，最終失去暫態(tài)穩(wěn)定。在極限切除時(shí)間之前切除故障，發(fā)電機(jī)的搖擺曲線的狀況將在下面作計(jì)算、分析。（6）不考慮阻尼轉(zhuǎn)矩的影響，當(dāng)故障切除時(shí)間為0.15秒時(shí)通過計(jì)算得出-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

22、/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,t0 tc,d0 w0);dc=

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