隱式改進(jìn)歐拉法在暫態(tài)穩(wěn)定分析中的應(yīng)用

蔣德瓏 尹淑萍 王克文 代旭光

(鄭州大學(xué)電氣工程學(xué)院，河南 鄭州 450001)

0 引言

經(jīng)濟(jì)和技術(shù)的進(jìn)步使電力系統(tǒng)得到了迅速發(fā)展，但也對電力系統(tǒng)穩(wěn)定性提出了更高的要求［1］。多年來，大量文獻(xiàn)研究了電力系統(tǒng)的專業(yè)化和自動化，在電力系統(tǒng)暫態(tài)穩(wěn)定自動化分析中更是取得了長足的進(jìn)步［2-11］。如采用計(jì)算出擾動后各發(fā)電機(jī)轉(zhuǎn)子間相對角度隨時(shí)間的變化曲線(即搖擺曲線)來判斷系統(tǒng)的穩(wěn)定性［6-9］;對暫態(tài)穩(wěn)定程序進(jìn)行模塊化處理，使程序具有較好的擴(kuò)展性和靈活性［10-11］。

相對顯式解法，隱式解法有其明顯的優(yōu)越性。在總結(jié)前人研究的基礎(chǔ)上，本文將經(jīng)典的改進(jìn)歐拉法與隱式解法相結(jié)合，用于電力系統(tǒng)的暫態(tài)穩(wěn)定計(jì)算，以有效地改善計(jì)算的精度和速度，提高暫態(tài)穩(wěn)定自動化分析的水平，為電力系統(tǒng)的安全穩(wěn)定運(yùn)行提供更為有力的保障。

1 牛頓法潮流計(jì)算

牛頓法是數(shù)學(xué)中解決非線性方程式的典型方法，它通過對非線性方程逐次線性化求解方程的根，并能保持較好的收斂性。潮流計(jì)算是研究電力系統(tǒng)暫態(tài)穩(wěn)定的基礎(chǔ)，它根據(jù)給定的發(fā)電運(yùn)行方式和系統(tǒng)接線方式來確定整個電力系統(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)，然后對節(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)為初值重新迭代。

2 穩(wěn)定分析數(shù)學(xué)模型的確定

整個電力系統(tǒng)的模型，在數(shù)學(xué)上可以統(tǒng)一描述成如下一般形式的微分-代數(shù)方程組:式中:x為微分方程組中描述系統(tǒng)動態(tài)特性的狀態(tài)變量;y為代數(shù)方程組中系統(tǒng)的運(yùn)行參量。本文重點(diǎn)關(guān)注發(fā)電機(jī)微分方程的求解。

2.1 發(fā)電機(jī)數(shù)學(xué)模型

考慮到暫態(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)動的機(jī)電暫態(tài)。

轉(zhuǎn)子電壓平衡方程為:

轉(zhuǎn)子運(yùn)動方程為:

式中:T'd為勵磁繞組本身的時(shí)間常數(shù);T″d為勵磁繞組短路、定子繞組開路時(shí)d軸阻尼繞組的時(shí)間常數(shù);T″q為定子繞組開路時(shí)q軸阻尼繞組的時(shí)間常數(shù);TJ為發(fā)電機(jī)組的慣性時(shí)間常數(shù);E'q為發(fā)電機(jī)q軸的暫態(tài)電動勢;E″q為發(fā)電機(jī)q軸的次暫態(tài)電動勢;E″d為發(fā)電機(jī)d軸的次暫態(tài)電動勢;δ為功角;ω為轉(zhuǎn)子角頻率。將式(2)與式(3)聯(lián)立，即可得到同步發(fā)電機(jī)的五階數(shù)學(xué)模型。

2.2 隱式改進(jìn)歐拉法

對于一般的非線性微分方程式:

確定時(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ù)方程的求解過程。因此，隱式解法有兩個明顯的優(yōu)點(diǎn):一是當(dāng)微分方程和代數(shù)方程需要聯(lián)立求解時(shí)，利用隱式解法便于消除交接誤差;二是隱式解法可以采取較大的步長。

2.3 差分方程

對于電力系統(tǒng)暫態(tài)計(jì)算涉及到的發(fā)電機(jī)、電動機(jī)和控制器等數(shù)學(xué)模型中的微分方程，皆可以將其化為差分方程［13-14］。下面以同步發(fā)電機(jī)的數(shù)學(xué)模型為例，利用隱式改進(jìn)歐拉法，將五階發(fā)電機(jī)模型中的微分方程化為差分方程。

由此得到的用于描述各發(fā)電機(jī)、電動機(jī)、控制器等差分方程組，可與用于描述電力網(wǎng)絡(luò)的代數(shù)方程組聯(lián)立求解，發(fā)揮了牛頓法所具有的數(shù)值穩(wěn)定、無交接誤差和收斂速度快等特點(diǎn)。

3 初值計(jì)算及運(yùn)行參數(shù)的求解

根據(jù)上述潮流計(jì)算法，可以計(jì)算出電力系統(tǒ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)定分析的先決條件。

3.1 初值的計(jì)算

在進(jìn)行暫態(tài)穩(wěn)定分析之前，必須求出暫態(tài)過程計(jì)算所需要的初值條件，即首先要確定求解微分方程所需要的初值。

在簡化的暫態(tài)穩(wěn)定計(jì)算程序中，初值的計(jì)算包括求系統(tǒng)擾動瞬間發(fā)電機(jī)的暫態(tài)電勢、轉(zhuǎn)子角度、原動機(jī)的機(jī)械功率和等值導(dǎo)納等。這些參數(shù)在電力系統(tǒng)擾動的瞬間是不會突變的［10］。因此，可以由擾動前的正常運(yùn)行狀態(tài)計(jì)算得到。

3.2 與網(wǎng)絡(luò)方程聯(lián)立求解

在電力系統(tǒng)中，電力網(wǎng)絡(luò)將系統(tǒng)中看起來相互獨(dú)立的所有暫態(tài)元件聯(lián)系在一起，在暫態(tài)過程中的任一時(shí)刻，各暫態(tài)元件注入網(wǎng)絡(luò)的電流不但由其本身的特性決定，且整個電力網(wǎng)絡(luò)必須滿足基爾霍夫定律。其中，前者由各暫態(tài)元件自身的代數(shù)方程描述，后者反映在電力網(wǎng)絡(luò)方程中。因此，為了求解網(wǎng)絡(luò)方程，需要列出各暫態(tài)元件自身的代數(shù)方程，并對其進(jìn)行處理，從而與網(wǎng)絡(luò)方程聯(lián)立求解。將發(fā)電機(jī)作為電流源，其節(jié)點(diǎn)注入電流的表達(dá)式為:

對于負(fù)荷模型，只需對導(dǎo)納矩陣的相應(yīng)對角塊進(jìn)行變化，即可得到新的網(wǎng)絡(luò)方程。

4 算例分析

利用上述方法編寫暫態(tài)穩(wěn)定計(jì)算程序，對典型的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臺發(fā)電機(jī)、3個負(fù)荷和9條支路，頻率為60 Hz，支路數(shù)據(jù)和發(fā)電機(jī)參數(shù)見文獻(xiàn)［14］。

假設(shè)系統(tǒng)在0 s受到擾動，在線路5-7靠近母線7處發(fā)生三相接地短路，則故障在5個周波(約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ī)的功角和相對搖擺角，其中數(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ī)功角的相對角度如表3所示，即發(fā)電機(jī)相對搖擺角，計(jì)算時(shí)間步長取0.01 s，其中，下標(biāo)2-1表示發(fā)電機(jī)②相對于發(fā)電機(jī)①的角度，3-1表示發(fā)電機(jī)③相對于發(fā)電機(jī)①的角度。限于篇幅，表3中只截取了以0.4 s為間隔的數(shù)據(jù)，“*”表示極值時(shí)刻。

表3 相對搖擺角數(shù)據(jù)Tab.3 Data of relative swing angle

由表2和表3可知，雖然發(fā)電機(jī)的功角不斷增大，但其相對搖擺角的變化并不大，第一擺時(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ī)間的最大相對搖擺角均小于180°。

為了便于分析，利用Matlab7.1對2 s內(nèi)電力系統(tǒng)發(fā)生三相短路故障的功角暫態(tài)過程進(jìn)行仿真，發(fā)電機(jī)功角和相對搖擺角仿真結(jié)果如圖2和圖3所示。圖2表示系統(tǒng)故障過程中各發(fā)電機(jī)的功角變化曲線，反映出3臺發(fā)電機(jī)的功角逐漸增大，說明有功出力逐漸增大。圖3表示系統(tǒng)故障過程中各發(fā)電機(jī)的功角的相對變化，反映出在計(jì)及發(fā)電機(jī)凸極效應(yīng)時(shí)，各發(fā)電機(jī)間的最大相對搖擺角均小于180°。通過與圖2的比較可以發(fā)現(xiàn)，雖然3臺發(fā)電機(jī)的功角逐漸增大，但其相對搖擺角的變化并沒有超過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í)用隱式解法對微分方程的處理，也更有利于在暫態(tài)計(jì)算過程中對各變量進(jìn)行修正，使發(fā)電機(jī)功角的計(jì)算結(jié)果更加接近真實(shí)值。

5 結(jié)束語

本文采用隱式改進(jìn)歐拉法對電力系統(tǒng)暫態(tài)計(jì)算中各元件的數(shù)學(xué)模型進(jìn)行處理，通過聯(lián)立求解的方法進(jìn)行暫態(tài)穩(wěn)定計(jì)算，最后對電力系統(tǒng)多機(jī)運(yùn)行方式下發(fā)電機(jī)功角的暫態(tài)過程進(jìn)行了分析。

與傳統(tǒng)的改進(jìn)歐拉法相比，該方法具有積分精度高和無交接誤差的優(yōu)點(diǎn)，可以滿足長期仿真計(jì)算對數(shù)值精度的要求。同時(shí)，本文介紹的方法保持了電力系統(tǒng)內(nèi)在的模塊性，在程序設(shè)計(jì)上具有良好的可擴(kuò)展性與靈活性。由于選擇了較為理想的系統(tǒng)，忽略了勵磁系統(tǒng)和調(diào)相機(jī)對電力系統(tǒng)暫態(tài)的影響，采用線性化模型得到的功角搖擺曲線還存在一定的誤差，因此，下一步可以采用保留高階項(xiàng)的模型，以求進(jìn)一步改善計(jì)算精度。

［1］倪以信，陳壽孫，張寶霖.動態(tài)電力系統(tǒng)的理論和分析［M］.北京:清華大學(xué)出版社，2002.

［2］孫暉，趙菁.建立專業(yè)電力自動化監(jiān)控系統(tǒng)的探討［J］.自動化儀表，2003，24(6):15-18.

［3］唐雯，羅安，唐杰，等.電力監(jiān)控系統(tǒng)圖形繪制及數(shù)據(jù)處理一體化平臺［J］.自動化儀表，2007，28(8):56-58.

［4］張海生，范征宇.S7-300 PLC和VC++.NET在電力監(jiān)控系統(tǒng)中的應(yīng)用［J］.自動化儀表，2006，27(8):50-52.

［5］程若發(fā)，馮士芬，江曉舟.基于ARM的同步發(fā)電機(jī)故障錄波系統(tǒng)的設(shè)計(jì)［J］.自動化儀表，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)電場穿透功率極限計(jì)算方法及發(fā)展［J］.自動化儀表，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.