蔣德瓏 尹淑萍 王克文 代旭光
(鄭州大學(xué)電氣工程學(xué)院,河南 鄭州 450001)
經(jīng)濟(jì)和技術(shù)的進(jìn)步使電力系統(tǒng)得到了迅速發(fā)展,但也對(duì)電力系統(tǒng)穩(wěn)定性提出了更高的要求[1]。多年來,大量文獻(xiàn)研究了電力系統(tǒng)的專業(yè)化和自動(dòng)化,在電力系統(tǒng)暫態(tài)穩(wěn)定自動(dòng)化分析中更是取得了長足的進(jìn)步[2-11]。如采用計(jì)算出擾動(dòng)后各發(fā)電機(jī)轉(zhuǎn)子間相對(duì)角度隨時(shí)間的變化曲線(即搖擺曲線)來判斷系統(tǒng)的穩(wěn)定性[6-9];對(duì)暫態(tài)穩(wěn)定程序進(jìn)行模塊化處理,使程序具有較好的擴(kuò)展性和靈活性[10-11]。
相對(duì)顯式解法,隱式解法有其明顯的優(yōu)越性。在總結(jié)前人研究的基礎(chǔ)上,本文將經(jīng)典的改進(jìn)歐拉法與隱式解法相結(jié)合,用于電力系統(tǒng)的暫態(tài)穩(wěn)定計(jì)算,以有效地改善計(jì)算的精度和速度,提高暫態(tài)穩(wěn)定自動(dòng)化分析的水平,為電力系統(tǒng)的安全穩(wěn)定運(yùn)行提供更為有力的保障。
牛頓法是數(shù)學(xué)中解決非線性方程式的典型方法,它通過對(duì)非線性方程逐次線性化求解方程的根,并能保持較好的收斂性。潮流計(jì)算是研究電力系統(tǒng)暫態(tài)穩(wěn)定的基礎(chǔ),它根據(jù)給定的發(fā)電運(yùn)行方式和系統(tǒng)接線方式來確定整個(gè)電力系統(tǒng)各部分的運(yùn)行狀態(tài),包括各母線的電壓、各元件流過的功率和系統(tǒng)的功率損耗等,為暫態(tài)穩(wěn)定分析提供初始運(yùn)行狀態(tài)。當(dāng)前,牛頓法被廣泛應(yīng)用于電力系統(tǒng)潮流計(jì)算,其計(jì)算的核心問題是修正方程式的建立和求解[12]。本文采用經(jīng)典的牛頓法潮流計(jì)算方法,其求解過程大致分為以下幾步:首先給定節(jié)點(diǎn)導(dǎo)納矩陣YB和節(jié)點(diǎn)電壓初值e(0)、f(0),代入功率偏差量的表達(dá)式中,求解 ΔP(0)、ΔQ(0)、(ΔV2)(0),再求出修正方程式的系數(shù)矩陣(雅可比矩陣)各元素;接著通過解修正方程式得到修正量Δe(0)、Δf(0),然后對(duì)節(jié)點(diǎn)電壓 e(1)=e(0)-Δe(0)、f(1)=f(0)-Δf(0)進(jìn)行修正,進(jìn)而利用 e(1)、f(1)求得 ΔP(1)、ΔQ(1)、(ΔV2)(1);最后判斷結(jié)果是否收斂,如收斂,則求出各支路潮流并輸出結(jié)果,否則以e(1)、f(1)為初值重新迭代。
整個(gè)電力系統(tǒng)的模型,在數(shù)學(xué)上可以統(tǒng)一描述成如下一般形式的微分-代數(shù)方程組:式中:x為微分方程組中描述系統(tǒng)動(dòng)態(tài)特性的狀態(tài)變量;y為代數(shù)方程組中系統(tǒng)的運(yùn)行參量。本文重點(diǎn)關(guān)注發(fā)電機(jī)微分方程的求解。
考慮到暫態(tài)穩(wěn)定的計(jì)算精度和系統(tǒng)的復(fù)雜性,本文采用忽略定子電磁暫態(tài)但考慮轉(zhuǎn)子阻尼繞組作用的五階發(fā)電機(jī)模型(E'q、E″q、E'd、δ、ω),即考慮 f繞組、D繞組、Q繞組的電磁暫態(tài)以及轉(zhuǎn)子運(yùn)動(dòng)的機(jī)電暫態(tài)。
轉(zhuǎn)子電壓平衡方程為:
轉(zhuǎn)子運(yùn)動(dòng)方程為:
式中:T'd為勵(lì)磁繞組本身的時(shí)間常數(shù);T″d為勵(lì)磁繞組短路、定子繞組開路時(shí)d軸阻尼繞組的時(shí)間常數(shù);T″q為定子繞組開路時(shí)q軸阻尼繞組的時(shí)間常數(shù);TJ為發(fā)電機(jī)組的慣性時(shí)間常數(shù);E'q為發(fā)電機(jī)q軸的暫態(tài)電動(dòng)勢(shì);E″q為發(fā)電機(jī)q軸的次暫態(tài)電動(dòng)勢(shì);E″d為發(fā)電機(jī)d軸的次暫態(tài)電動(dòng)勢(shì);δ為功角;ω為轉(zhuǎn)子角頻率。將式(2)與式(3)聯(lián)立,即可得到同步發(fā)電機(jī)的五階數(shù)學(xué)模型。
對(duì)于一般的非線性微分方程式:
確定時(shí)間步長 Δt=h,由(tn,xn)點(diǎn)推算(tn+1,xn+1)點(diǎn)時(shí),改進(jìn)歐拉法的遞推公式即為:
改進(jìn)歐拉法的核心思想是取某時(shí)段始點(diǎn)導(dǎo)數(shù)值與終點(diǎn)導(dǎo)數(shù)值的平均值作為該時(shí)刻折線段的斜率,進(jìn)而得到更加精確的計(jì)算結(jié)果。而用f(xn,tn)和f(xn+1,tn+1)代替積分曲線[tn,tn+1]區(qū)間始點(diǎn)和終點(diǎn)的斜率,即可將式(5)進(jìn)一步變換為:
式(6)中只有xn+1為未知數(shù),因而利用求解代數(shù)方程式的方法就可以求得xn+1,這就得到了隱式改進(jìn)歐拉法。
傳統(tǒng)的改進(jìn)歐拉法是利用已知量,根據(jù)遞推公式逐次遞推計(jì)算出相應(yīng)時(shí)段終點(diǎn)的函數(shù)值,因此是一種顯式解法。隱式改進(jìn)歐拉法不是給出遞推公式,而是首先把微分方程化為差分方程,然后利用求解差分方程的方法確定函數(shù)值,其特點(diǎn)就是把微分方程的求解問題轉(zhuǎn)換成一系列代數(shù)方程的求解過程。因此,隱式解法有兩個(gè)明顯的優(yōu)點(diǎn):一是當(dāng)微分方程和代數(shù)方程需要聯(lián)立求解時(shí),利用隱式解法便于消除交接誤差;二是隱式解法可以采取較大的步長。
對(duì)于電力系統(tǒng)暫態(tài)計(jì)算涉及到的發(fā)電機(jī)、電動(dòng)機(jī)和控制器等數(shù)學(xué)模型中的微分方程,皆可以將其化為差分方程[13-14]。下面以同步發(fā)電機(jī)的數(shù)學(xué)模型為例,利用隱式改進(jìn)歐拉法,將五階發(fā)電機(jī)模型中的微分方程化為差分方程。
由此得到的用于描述各發(fā)電機(jī)、電動(dòng)機(jī)、控制器等差分方程組,可與用于描述電力網(wǎng)絡(luò)的代數(shù)方程組聯(lián)立求解,發(fā)揮了牛頓法所具有的數(shù)值穩(wěn)定、無交接誤差和收斂速度快等特點(diǎn)。
根據(jù)上述潮流計(jì)算法,可以計(jì)算出電力系統(tǒng)遭受擾動(dòng)以前的正常運(yùn)行狀態(tài),得到電力網(wǎng)絡(luò)的運(yùn)行參數(shù),并由此求出電力系統(tǒng)中有關(guān)元件的狀態(tài)參數(shù)初始值。這些都是進(jìn)行電力系統(tǒng)暫態(tài)穩(wěn)定分析的先決條件。
在進(jìn)行暫態(tài)穩(wěn)定分析之前,必須求出暫態(tài)過程計(jì)算所需要的初值條件,即首先要確定求解微分方程所需要的初值。
在簡化的暫態(tài)穩(wěn)定計(jì)算程序中,初值的計(jì)算包括求系統(tǒng)擾動(dòng)瞬間發(fā)電機(jī)的暫態(tài)電勢(shì)、轉(zhuǎn)子角度、原動(dòng)機(jī)的機(jī)械功率和等值導(dǎo)納等。這些參數(shù)在電力系統(tǒng)擾動(dòng)的瞬間是不會(huì)突變的[10]。因此,可以由擾動(dòng)前的正常運(yùn)行狀態(tài)計(jì)算得到。
在電力系統(tǒng)中,電力網(wǎng)絡(luò)將系統(tǒng)中看起來相互獨(dú)立的所有暫態(tài)元件聯(lián)系在一起,在暫態(tài)過程中的任一時(shí)刻,各暫態(tài)元件注入網(wǎng)絡(luò)的電流不但由其本身的特性決定,且整個(gè)電力網(wǎng)絡(luò)必須滿足基爾霍夫定律。其中,前者由各暫態(tài)元件自身的代數(shù)方程描述,后者反映在電力網(wǎng)絡(luò)方程中。因此,為了求解網(wǎng)絡(luò)方程,需要列出各暫態(tài)元件自身的代數(shù)方程,并對(duì)其進(jìn)行處理,從而與網(wǎng)絡(luò)方程聯(lián)立求解。將發(fā)電機(jī)作為電流源,其節(jié)點(diǎn)注入電流的表達(dá)式為:
對(duì)于負(fù)荷模型,只需對(duì)導(dǎo)納矩陣的相應(yīng)對(duì)角塊進(jìn)行變化,即可得到新的網(wǎng)絡(luò)方程。
利用上述方法編寫暫態(tài)穩(wěn)定計(jì)算程序,對(duì)典型的9節(jié)點(diǎn)三機(jī)系統(tǒng)進(jìn)行仿真計(jì)算,三機(jī)系統(tǒng)示意圖如圖1所示。
圖1 三機(jī)系統(tǒng)示意圖Fig.1 Three-machine system
圖1中,系統(tǒng)有3臺(tái)發(fā)電機(jī)、3個(gè)負(fù)荷和9條支路,頻率為60 Hz,支路數(shù)據(jù)和發(fā)電機(jī)參數(shù)見文獻(xiàn)[14]。
假設(shè)系統(tǒng)在0 s受到擾動(dòng),在線路5-7靠近母線7處發(fā)生三相接地短路,則故障在5個(gè)周波(約0.08333 s)后由斷開線路5-7而被消除。
本文采用牛頓法潮流計(jì)算得出正常運(yùn)行情況下的系統(tǒng)潮流,如表1所示。以其作為初始運(yùn)行條件,由隱式改進(jìn)歐拉法暫態(tài)穩(wěn)定分析程序可計(jì)算各發(fā)電機(jī)功角在故障過程中的變化情況。
表1 潮流計(jì)算結(jié)果Tab.1 Results of load flow calculation
在暫態(tài)穩(wěn)定計(jì)算程序中,各發(fā)電機(jī)用五階模型模擬,各負(fù)荷用恒定阻抗模擬,電力網(wǎng)絡(luò)用導(dǎo)納矩陣描述,微分方程用隱式改進(jìn)歐拉法求解,網(wǎng)絡(luò)方程用直接法求解。程序計(jì)算了從短路故障開始后,系統(tǒng)在2 s內(nèi)的暫態(tài)過程,得出計(jì)及發(fā)電機(jī)凸極效應(yīng)時(shí)各發(fā)電機(jī)的功角和相對(duì)搖擺角,其中數(shù)值積分采用0.001 s的步長。
系統(tǒng)故障過程中各發(fā)電機(jī)的功角數(shù)據(jù)如表2所示,計(jì)算時(shí)間步長取 0.01 s。表 2 中,1-1、2-2、3-3分別表示發(fā)電機(jī)①、②、③的功角。限于篇幅,表2中只截取了以0.2 s為間隔的數(shù)據(jù)。
表2 功角數(shù)據(jù)Tab.2 Data of power angle
系統(tǒng)故障過程中各發(fā)電機(jī)功角的相對(duì)角度如表3所示,即發(fā)電機(jī)相對(duì)搖擺角,計(jì)算時(shí)間步長取0.01 s,其中,下標(biāo)2-1表示發(fā)電機(jī)②相對(duì)于發(fā)電機(jī)①的角度,3-1表示發(fā)電機(jī)③相對(duì)于發(fā)電機(jī)①的角度。限于篇幅,表3中只截取了以0.4 s為間隔的數(shù)據(jù),“*”表示極值時(shí)刻。
表3 相對(duì)搖擺角數(shù)據(jù)Tab.3 Data of relative swing angle
由表2和表3可知,雖然發(fā)電機(jī)的功角不斷增大,但其相對(duì)搖擺角的變化并不大,第一擺時(shí)的最大搖擺角度 為 δ2-1=126.47259°(t=0.59s)、δ3-1=94.32671°(t=0.64 s),第二擺時(shí)的最大搖擺角度為δ2-1=126.74058°(t=1.88 s)、δ3-1=95.87463°(t=1.92 s),比第一擺時(shí)的角度稍大,2 s內(nèi)各發(fā)電機(jī)間的最大相對(duì)搖擺角均小于180°。
為了便于分析,利用Matlab7.1對(duì)2 s內(nèi)電力系統(tǒng)發(fā)生三相短路故障的功角暫態(tài)過程進(jìn)行仿真,發(fā)電機(jī)功角和相對(duì)搖擺角仿真結(jié)果如圖2和圖3所示。圖2表示系統(tǒng)故障過程中各發(fā)電機(jī)的功角變化曲線,反映出3臺(tái)發(fā)電機(jī)的功角逐漸增大,說明有功出力逐漸增大。圖3表示系統(tǒng)故障過程中各發(fā)電機(jī)的功角的相對(duì)變化,反映出在計(jì)及發(fā)電機(jī)凸極效應(yīng)時(shí),各發(fā)電機(jī)間的最大相對(duì)搖擺角均小于180°。通過與圖2的比較可以發(fā)現(xiàn),雖然3臺(tái)發(fā)電機(jī)的功角逐漸增大,但其相對(duì)搖擺角的變化并沒有超過180°,說明系統(tǒng)在發(fā)生該故障時(shí)仍能保持穩(wěn)定,即系統(tǒng)是暫態(tài)穩(wěn)定的,與表2和表3中數(shù)據(jù)反映的情況一致。
算例結(jié)果證明,該方法能夠較好地應(yīng)用于電力系統(tǒng)暫態(tài)穩(wěn)定的判定過程,同時(shí)用隱式解法對(duì)微分方程的處理,也更有利于在暫態(tài)計(jì)算過程中對(duì)各變量進(jìn)行修正,使發(fā)電機(jī)功角的計(jì)算結(jié)果更加接近真實(shí)值。
本文采用隱式改進(jìn)歐拉法對(duì)電力系統(tǒng)暫態(tài)計(jì)算中各元件的數(shù)學(xué)模型進(jìn)行處理,通過聯(lián)立求解的方法進(jìn)行暫態(tài)穩(wěn)定計(jì)算,最后對(duì)電力系統(tǒng)多機(jī)運(yùn)行方式下發(fā)電機(jī)功角的暫態(tài)過程進(jìn)行了分析。
與傳統(tǒng)的改進(jìn)歐拉法相比,該方法具有積分精度高和無交接誤差的優(yōu)點(diǎn),可以滿足長期仿真計(jì)算對(duì)數(shù)值精度的要求。同時(shí),本文介紹的方法保持了電力系統(tǒng)內(nèi)在的模塊性,在程序設(shè)計(jì)上具有良好的可擴(kuò)展性與靈活性。由于選擇了較為理想的系統(tǒng),忽略了勵(lì)磁系統(tǒng)和調(diào)相機(jī)對(duì)電力系統(tǒng)暫態(tài)的影響,采用線性化模型得到的功角搖擺曲線還存在一定的誤差,因此,下一步可以采用保留高階項(xiàng)的模型,以求進(jìn)一步改善計(jì)算精度。
[1]倪以信,陳壽孫,張寶霖.動(dòng)態(tài)電力系統(tǒng)的理論和分析[M].北京:清華大學(xué)出版社,2002.
[2]孫暉,趙菁.建立專業(yè)電力自動(dòng)化監(jiān)控系統(tǒng)的探討[J].自動(dòng)化儀表,2003,24(6):15-18.
[3]唐雯,羅安,唐杰,等.電力監(jiān)控系統(tǒng)圖形繪制及數(shù)據(jù)處理一體化平臺(tái)[J].自動(dòng)化儀表,2007,28(8):56-58.
[4]張海生,范征宇.S7-300 PLC和VC++.NET在電力監(jiān)控系統(tǒng)中的應(yīng)用[J].自動(dòng)化儀表,2006,27(8):50-52.
[5]程若發(fā),馮士芬,江曉舟.基于ARM的同步發(fā)電機(jī)故障錄波系統(tǒng)的設(shè)計(jì)[J].自動(dòng)化儀表,2007,28(9):11-14.
[6]西安交通大學(xué),清華大學(xué).電力系統(tǒng)計(jì)算[M].北京:水利電力出版社,1978.
[7]Al-Taee M A,Al-Azzawi F J,Al-Taee A A,et al.Real-time assessment of power system transient stability using rate of change of kinetic energy method[C]//IEEE Proceedings of Generation Transmission and Distribution,2001,148(6):505-510.
[8]吳穎,趙巖,蔣傳文,等.風(fēng)電場(chǎng)穿透功率極限計(jì)算方法及發(fā)展[J].自動(dòng)化儀表,2008,29(11):7-11.
[9]趙艷雷,童建忠.改進(jìn)歐拉法在電力系統(tǒng)暫態(tài)分析中的應(yīng)用和軟件設(shè)計(jì)[J].微計(jì)算機(jī)信息,2004,20(5):98-99.
[10]Kasztenny B,Kezunovic M.A method for linking different modeling techniques for accurate and efficient simulation[J].IEEE Transactions on Power Systems,2000,15(1):65-72.
[11]Suyono H,Nor K M,Yusof S.Component based development for transient stability power system simulation software[C]//Proceedings of Power Engineering Conference,PECon,2003:16-21.
[12]陳珩.電力系統(tǒng)穩(wěn)態(tài)分析[M].2版.北京:中國電力出版社,1995.
[13]武東亞,王克文,馬潔.基于概率潮流的多運(yùn)行方式下發(fā)電機(jī)功角的暫態(tài)過程分析[J].電氣應(yīng)用,2009,28(6):64-67.
[14]王錫凡,方萬良,杜正春.現(xiàn)代電力系統(tǒng)分析[M].北京:科學(xué)出版社,2003.