幾何數(shù)值積分方法在暫態(tài)穩(wěn)定性計算中的應(yīng)用

汪芳宗，郭夢芳，宋墩文，楊學(xué)濤，張磊

(1.三峽大學(xué)電氣與新能源學(xué)院，湖北宜昌443002；2.中國電力科學(xué)研究院， 北京100040)

0 引言

傳統(tǒng)的數(shù)值積分方法主要包括線性多步法(linear multistep formulae, LMF)[1-2]以及Runge-Kutta(RK)系列方法[2-3]。電力系統(tǒng)暫態(tài)仿真計算中應(yīng)用最為廣泛的隱式梯形積分，是RK系列中的一種單級、2階方法[4]；而另一種應(yīng)用較多的方法是2階的向后差分方法(backward differentiation formulae, BDF)，此方法屬于LMF類方法[2]。

關(guān)于傳統(tǒng)的數(shù)值積分方法，研究人員很早就建立了有關(guān)數(shù)值方法的穩(wěn)定性以及其計算精度分析的理論體系，并給出了明確、具體的結(jié)論[3]。從工程實際應(yīng)用的角度看，研究人員在選擇傳統(tǒng)的數(shù)值積分方法時主要是關(guān)注數(shù)值積分方法的數(shù)值穩(wěn)定性，其次就是計算精度[5-7]。事實上，隱式梯形積分方法之所以在電力系統(tǒng)暫態(tài)仿真計算中得到廣泛的應(yīng)用，主要依據(jù)是該方法為具有最小局部誤差截斷系數(shù)的A-穩(wěn)定的數(shù)值計算方法，而2階的BDF方法則是LMF系列方法中具有A-穩(wěn)定性以及L-穩(wěn)定性的最高階方法。

與傳統(tǒng)的數(shù)值積分方法不同，幾何積分方法主要是關(guān)注動力系統(tǒng)的物理/幾何屬性。如果數(shù)值積分方法能精確地保留動力系統(tǒng)的一個或多個物理/幾何屬性，則稱為幾何數(shù)值積分方法。幾何屬性的具體示例大致包括動力系統(tǒng)的對稱性以及反對稱性、辛結(jié)構(gòu)、相空間中的體積(phase-space volume)等，而物理屬性主要涉及動力系統(tǒng)的能量守恒性，例如首次積分、Hamilton能量函數(shù)、動量守恒等，因此，保幾何屬性的數(shù)值積分方法通常稱為保結(jié)構(gòu)計算方法(structure-preserving methods)[8,9]，而保物理屬性的數(shù)值積分方法通常稱為保能量守恒方法，簡稱為保能量方法(energy-preserving methods)[10]。

幾何積分方法的興起與發(fā)展得益于我國學(xué)者馮康先生的卓越貢獻(xiàn)。早在20世紀(jì)80年代中期，馮康先生率先提出并系統(tǒng)建立了Hamilton動力系統(tǒng)的辛幾何算法，也就是保Hamilton動力系統(tǒng)所具有的辛結(jié)構(gòu)的保辛方法[11-12]。由于Hamilton動力系統(tǒng)具有辛結(jié)構(gòu)與Hamilton能量守恒2大特征，因此，在馮康先生提出辛幾何方法以后，研究人員自然就會想到并開始研究保Hamilton能量的數(shù)值積分方法[13]。換言之，馮康先生不僅創(chuàng)始了保結(jié)構(gòu)計算方法，對保能量積分方法的研究也有重要的啟發(fā)和推動作用。迄今為止，辛幾何方法已在諸多領(lǐng)域得到了廣泛應(yīng)用。除辛幾何方法以外，馮康先生還率先提出了無源動力系統(tǒng)的保積方法(volume-preserving methods)[8,14]，亦稱為無散度或零散度動力系統(tǒng)。由于無源動力系統(tǒng)的真解也就是相流在相空間中是保體積的，因此，保積方法也是一類重要的保結(jié)構(gòu)數(shù)值方法[15-16]。

關(guān)于幾何積分方法的應(yīng)用，一種思路是用幾何積分方法來替換傳統(tǒng)的數(shù)值積分方法。需要說明的是，盡管幾何積分方法是基于動力系統(tǒng)的物理/幾何屬性而發(fā)展起來的，但這并不意味著幾何積分方法不能用于一些不具備物理/幾何屬性的一般動力系統(tǒng)。事實上，在辛幾何方法發(fā)明以前，RK系列方法中的Gauss類方法很早就在科學(xué)及工程計算中得到了應(yīng)用，而Gauss方法就是一類辛幾何方法[17]。另外一種思路就是：對具有特定的物理/幾何屬性微分動力系統(tǒng)，數(shù)值積分方法應(yīng)該也必須保持相應(yīng)的特征。換言之，對具有特定的物理/幾何屬性的微分動力系統(tǒng)，最好或必須采用相應(yīng)的幾何積分方法。這一思路是馮康先生提出的一個重要的學(xué)術(shù)思想[8]。

基于馮康先生的學(xué)術(shù)思想，本文嘗試將幾何積分方法應(yīng)用于電力系統(tǒng)暫態(tài)穩(wěn)定性計算。在暫態(tài)穩(wěn)定性計算中，單機—無窮大系統(tǒng)既是一個經(jīng)典的Hamilton系統(tǒng)，也是一個可分的無源動力系統(tǒng)。基于此思想，本文首先將St?rmer-Verlet方法(也稱為Leapfrog方法)、顯辛RKN(Runge-Kutta-Nystr?m)方法[8]、簡化的Takahashi-Imada方法(以下簡稱STI方法)[18]、平均向量場方法(averaged vector field method, AVF)[10]等應(yīng)用于單機—無窮大系統(tǒng)的暫態(tài)穩(wěn)定性計算中，并將上述方法與傳統(tǒng)的數(shù)值積分方法進行對比分析。對單機—無窮大系統(tǒng)，Leapfrog方法以及顯辛RKN方法既是保辛的，也是保積的；STI方法則只是保積的方法；而AVF方法則是保Hamilton能量或首次積分的數(shù)值方法。

在上述基礎(chǔ)上，本文嘗試將幾何積分方法應(yīng)用于多機系統(tǒng)的暫態(tài)穩(wěn)定性計算。在采用經(jīng)典模型的前提下，多機系統(tǒng)暫態(tài)穩(wěn)定性計算的數(shù)學(xué)模型是一個可分的無源動力系統(tǒng)，但不是一個Hamilton系統(tǒng)，也就不具備辛結(jié)構(gòu)及Hamilton能量守恒性。本文將Leapfrog方法、顯辛RKN方法以及STI方法應(yīng)用于多機系統(tǒng)的暫態(tài)穩(wěn)定性計算，并將其與傳統(tǒng)的隱式梯形積分方法進行對比。針對可分無源動力系統(tǒng)，這3種數(shù)值方法均是保積算法。

1 幾何積分方法

1.1 可分Hamilton系統(tǒng)與顯辛方法

設(shè)有2n個自變量p=[p1,p2,…,pn]T，q=[q1,q2,…,qn]T其中pi=pi(t)，qi=qi(t)。給定一個連續(xù)可微函數(shù)H=H(p,q)，稱其為Hamilton函數(shù)或能量，則由該函數(shù)所確定的Hamilton系統(tǒng)為

(1)

或?qū)懗伤暮啙嵭问?亦稱正則方程)：

(2)

(3)

其中，In為n維的單位陣，J通常稱為標(biāo)準(zhǔn)辛矩陣。若H=H(p,q)=U(p)+V(q)，則相應(yīng)的Hamilton系統(tǒng)稱為可分Hamilton系統(tǒng)。在可分Hamilton系統(tǒng)中，有一類特殊的Hamilton系統(tǒng)，具體可表述為

(4)

相應(yīng)的Hamilton函數(shù)為

(5)

(6)

在不知道具體的Hamilton函數(shù)H(p,q)的情況下，式(6)是判斷動力系統(tǒng)式(1)是否為Hamilton系統(tǒng)的準(zhǔn)則。

針對Hamilton系統(tǒng)，研究人員已構(gòu)造出了多種多樣的辛積分格式。對一般的Hamilton系統(tǒng)，不存在顯式的辛RK方法[8]，但對可分Hamilton系統(tǒng)，研究人員構(gòu)造并發(fā)現(xiàn)了幾類顯式辛格式。由于隱式方法涉及的計算量較大，因此，本文主要考慮顯式辛方法。以下介紹的顯式方法均是辛方法。

(1)Leapfrog方法，具體可表述為

(7)

對形如式(4)的可分Hamilton系統(tǒng)，上述方法是2階、對稱的辛方法，其共軛或伴隨方法也是2階的辛方法。

(2)顯辛RKN方法，其相應(yīng)的Butcher表述為

3+3603-362+31203+3603605-33243+3121+3243-23121123+2312

(8)

對形如式(4)的可分Hamilton系統(tǒng)，上述方法是3級4階的辛方法，其共軛或伴隨方法也是4階的辛方法。此類方法是由秦孟兆先生率先提出的[19]。

1.2 可分無源動力系統(tǒng)與保積方法

令n個自變量x?[x1,x2,…,xn]T，f(x)?[f1(x),f2(x),…,fn(x)]T。若有

(9)

無源動力系統(tǒng)具有一個重要的幾何屬性：其相流在相空間中是保體積的，即保持Ω=dx1∧dx2∧…∧dxn不變，因此，對無源動力系統(tǒng)，相應(yīng)的數(shù)值積分方法也應(yīng)該是保積的。對一個數(shù)值積分方法xn→xn+1，若它所生成的映射滿足

(10)

則該數(shù)值積分方法就是保積方法[14-15]。式中det (■)表示矩陣的行列式。

關(guān)于無源動力系統(tǒng)的保積方法，馮康等在文獻(xiàn)[14]中已證明：對一般的無源動力系統(tǒng)，傳統(tǒng)的數(shù)值方法例如LMF方法以及RK方法均不可能是保積的，對一般的無源動力系統(tǒng)，構(gòu)造保積算法是比較困難的一件事，但是，對一類特殊的無源動力系統(tǒng)也就是可分無源動力系統(tǒng)，研究人員發(fā)現(xiàn)并構(gòu)造了眾多的保積方法，而且均是顯式方法。若一個動力系統(tǒng)可表述為

(11)

則該系統(tǒng)是可分無源動力系統(tǒng)，并且可分無源動力系統(tǒng)滿足下列方程：

(12)

關(guān)于可分無源動力系統(tǒng)的保積方法，本文主要介紹以下幾種保積方法：

①對形如式(4)的動力系統(tǒng)，Leapfrog方法式(7)是2階的保積方法[18]。

②對形如式(4)的動力系統(tǒng)，RKN方法式(8)是4階的保積方法[20]。

③對形如式(4)的動力系統(tǒng)，STI方法是4階的保積方法[18]，其具體計算格式可表述為

(13)

式中α=1/12。這里需要說明的是：對可分Hamilton系統(tǒng)(也是可分無源動力系統(tǒng))，STI方法只是一種保積方法，但一般情況下它不是一種辛方法[18]。

1.3 平均向量場方法

針對Hamilton系統(tǒng)的能量守恒性，研究人員提出了保Hamilton能量的平均向量場方法[10]。以Hamilton正則式(2)為例，平均向量場方法可表述為

(14)

平均向量場方法是2階的隱式方法，理論上它可以精確保持Hamilton系統(tǒng)的能量守恒性。相較于顯式積分方法，平均向量場方法的計算量較大，并且，該方法只是保能量的數(shù)值方法，不是辛方法。

2 單機—無窮大系統(tǒng)暫態(tài)穩(wěn)定性計算

2.1 數(shù)學(xué)模型及其幾何/物理特性

圖1 單機—無窮大系統(tǒng)Fig.1 Single machine-infinity system

如圖1所示的單機—無窮大系統(tǒng)，在發(fā)電機采用經(jīng)典模型的情況下，相應(yīng)的暫態(tài)穩(wěn)定性計算的數(shù)學(xué)模型可表述為

(15)

式中各變量的物理意義或定義可參見有關(guān)電力系統(tǒng)分析的教科書，其中：M=11.28 s；Pm=1.0 p.u.；Ev/x∑=1.384 p.u.。

上述單機—無窮大系統(tǒng)結(jié)構(gòu)雖然簡單，但具備了一些重要的特征：此系統(tǒng)滿足式(6)，既是一個經(jīng)典的可分Hamilton系統(tǒng)，也是一個可分無源動力系統(tǒng)。其Hamilton能量函數(shù)表示為

(16)

式中δs是系統(tǒng)的穩(wěn)定平衡點。將δs代入式(16)，可知Hamilton能量函數(shù)是守恒的[21]。即

(17)

可知，上述Hamilton能量函數(shù)式(16)也是單機—無窮大系統(tǒng)的首次積分。

2.2 數(shù)值測試結(jié)果

由于單機—無窮大系統(tǒng)具備上述幾何/物理屬性，因此，本文首先將Leapfrog方法式(7)、顯辛RKN方法式(8)、STI方法式(13)、AVF方法式(14)分別應(yīng)用于該系統(tǒng)的暫態(tài)穩(wěn)定性計算，并將這4種方法的仿真結(jié)果與傳統(tǒng)的隱式梯形積分方法的仿真結(jié)果進行了對比分析。

圖2是3種不同的2階數(shù)值方法計算出的相空間軌跡圖。數(shù)值試驗中，步長設(shè)定為h=0.02 s，積分初始值設(shè)定為δ0=1.111 2 rad，ω0=0.012 7 p.u.，此時系統(tǒng)處于臨界穩(wěn)定狀態(tài)。理論上，若單機—無窮大系統(tǒng)是暫態(tài)穩(wěn)定的，則其相空間軌跡是一個封閉的曲線[21]。

為測試數(shù)值方法的穩(wěn)定性，對上述單機—無窮大系統(tǒng)進行長時間數(shù)值積分計算。得到的相空間軌跡如圖2所示。如圖2(a)所示，是隱式梯形積分法的相空間軌跡曲線，在經(jīng)過157.4 s的長時間積分計算后，系統(tǒng)相空間軌跡發(fā)生破損(相空間軌跡不再是封閉的)，這是由于長時間仿真計算導(dǎo)致了錯誤的結(jié)果。如圖2(b)所示，是利用AVF方法所計算出的相空間軌跡曲線，在經(jīng)過14 172.8 s的長時間積分計算后，系統(tǒng)相空間軌跡發(fā)生破損，系統(tǒng)從穩(wěn)定狀態(tài)變?yōu)椴环€(wěn)定狀態(tài)。如圖2(c)所示，是利用Leapfrog方法所計算出的相空間軌跡曲線，在經(jīng)過14 200 s的長時間積分計算后，系統(tǒng)仍然能夠保持穩(wěn)定狀態(tài)。因此，從長時間計算數(shù)值穩(wěn)定性的角度可知，盡管Leapfrog方法是顯式方法，但該方法在保結(jié)構(gòu)穩(wěn)定性方面效果最好，AVF在保結(jié)構(gòu)穩(wěn)定性方面方法效果次之，隱式梯形積分方法在保結(jié)構(gòu)穩(wěn)定性方面效果最差。

(a) 隱式梯形積分方法

(b) AVF方法

(c) Leapfrog方法

圖3 (2階方法)Hamilton函數(shù)差值曲線 Fig.3 (2nd order method) Hamilton function difference curve

圖3是3種不同方法計算出的Hamilton函數(shù)差值(用ΔH(t)表示)曲線，其中，ΔH(t)?H(δ(t),ω(t))-H0，H0?(δ0,ω0)。Hamilton函數(shù)差值曲線反映數(shù)值積分方法的保能量特性。如圖3所示，根據(jù)能量差值的波動幅度可知，在保能量守恒方面，AVF方法的效果最好，Leapfrog方法的效果次之，隱式梯形積分方法的效果最差。

圖4是4階的RKN方法、STI方法以及傳統(tǒng)的4階顯式RK方法所計算出的相空間軌跡圖。其中，步長設(shè)定為h=0.05 s。測試結(jié)果顯示，經(jīng)過3 000 s的長時間積分計算后，提出的3種方法仍能保持相空間軌跡的封閉性。如圖4(c)所示，傳統(tǒng)顯式RK方法的相空間軌跡在長時間積分后產(chǎn)生了“漂移”現(xiàn)象，可知其在保結(jié)構(gòu)穩(wěn)定性方面比保積的RKN方法以及STI方法效果要差。

圖5是3種不同方法計算出的Hamilton函數(shù)差值曲線。圖5(b)是RKN方法以及傳統(tǒng)顯式RK方法的Hamilton函數(shù)差值曲線。由圖5可知，根據(jù)能量差值的波動幅度可知，在保能量守恒方面，RKN方法效果最好，傳統(tǒng)RK方法效果次之，STI方法的能量差值變化幅度最大、其效果最差。

(a) RKN方法

(b) STI方法

(c) 傳統(tǒng)4階顯式RK方法

(a) 3種方法對比結(jié)果

與傳統(tǒng)的概念有所不同，通過對單機—無窮大系統(tǒng)的數(shù)值仿真分析，發(fā)現(xiàn)并驗證了可分Hamilton系統(tǒng)或可分無源動力系統(tǒng)，顯式保結(jié)構(gòu)計算方法具備良好數(shù)值穩(wěn)定性。

3 基于保積方法的多機系統(tǒng)暫態(tài)穩(wěn)定性計算

基于顯式保結(jié)構(gòu)計算方法的優(yōu)勢，本節(jié)將顯式保積方法應(yīng)用于多機系統(tǒng)暫態(tài)穩(wěn)定性數(shù)值計算，并將其與常用的隱式梯形積分方法進行對比。

3.1 經(jīng)典模型及其幾何屬性

在負(fù)荷采用恒阻抗模型、發(fā)電機采用經(jīng)典模型的前提下，基于收縮導(dǎo)納矩陣的電力系統(tǒng)暫態(tài)穩(wěn)定性計算的經(jīng)典模型為

(18)

式中：m表示發(fā)電機的數(shù)量；Pei為電磁功率，其表達(dá)式可以寫為

(19)

與單機—無窮大系統(tǒng)不同，上述多機系統(tǒng)的經(jīng)典模型不是一個Hamilton系統(tǒng)，經(jīng)驗證可知，式(18)并不滿足式(6)。多機系統(tǒng)的首次積分滿足暫態(tài)穩(wěn)定性分析中常用的暫態(tài)能量函數(shù)[21]，但不是一個守恒量。運用辛幾何方法求解多機系統(tǒng)，雖然不具備“保結(jié)構(gòu)”和“保能量”的特征，但上述經(jīng)典模型仍具有重要的幾何屬性，即該系統(tǒng)是一個經(jīng)典的可分無源動力系統(tǒng)。

3.2 數(shù)值測試結(jié)果

暫態(tài)穩(wěn)定性計算的經(jīng)典模型是一個可分的無源動力系統(tǒng)，本節(jié)將Leapfrog方法、RKN方法以及STI方法分別應(yīng)用于多機系統(tǒng)暫態(tài)穩(wěn)定性計算，并將計算出的結(jié)果與傳統(tǒng)的積分方法進行對比分析。為此，利用隱式梯形方法采用小步長h=0.001 s進行數(shù)值仿真計算的結(jié)果作為基準(zhǔn)值，對比所使用算法進行數(shù)值仿真計算的結(jié)果與基準(zhǔn)值的差值Δδ(t)。

選用含50臺發(fā)電機的IEEE 145節(jié)點系統(tǒng)[22]作為算例。暫態(tài)穩(wěn)定性計算中，設(shè)定7號母線處發(fā)生三相短路，經(jīng)0.1 s后切除故障。這種情況下，系統(tǒng)是暫態(tài)穩(wěn)定的，其中，20號發(fā)電機受擾動影響最大，因此，將該臺發(fā)電機的功角計算結(jié)果用于分析、比較不同的數(shù)值計算方法。其中隱式梯形方法設(shè)置小步長h=0.001 s。

圖6是分別采用步長為h=0.02 s以及h=0.03 s時Leapfrog方法的差值曲線，圖7是采用步長為h=0.05 s時RKN方法、STI方法以及傳統(tǒng)4階顯式RK方法的差值曲線，圖8是采用步長為h=0.06 s時3種方法的差值曲線。

圖6 2階Leapfrog方法差值曲線Fig.6 Difference curve of the second-order Leapfrog method

圖7 4階方法差值曲線(h=0.05 s)Fig.7 Difference curve of 4th order method (h=0.05 s)

圖8 4階方法差值曲線(h=0.06 s) Fig.8 Difference curve of 4th order method (h=0.06 s)

從圖6可以看出，2階的Leapfrog方法在采用較大步長(h=0.03 s)的情況下，仍然可以獲得足夠精確的、穩(wěn)定的數(shù)值結(jié)果。從圖7、8可以看出，在采用大步長的情況下，保積的RKN方法在計算精度上效果最好。

表1是關(guān)于計算效率的測試結(jié)果。根據(jù)表1中的計算結(jié)果可知，在采用同樣步長的情況下，Leapfrog方法比隱式梯形積分方法約快3.3倍；在采用步長為隱式梯形積分方法的2倍情況下，RKN方法與STI方法大致相當(dāng)，但兩者均比傳統(tǒng)顯式RK方法略快，比隱式梯形積分方法約快5.8倍。

表1 計算效率測試結(jié)果Tab.1 Calculation efficiency test results

4 結(jié)語

本文通過將幾何積分方法引入電力系統(tǒng)分析，在單機—無窮大系統(tǒng)和暫態(tài)穩(wěn)定性計算的經(jīng)典模型中進行仿真計算，得到以下結(jié)論：

①對具有幾何屬性的動力系統(tǒng)，顯式保結(jié)構(gòu)數(shù)值積分方法具有很好的保結(jié)構(gòu)數(shù)值穩(wěn)定性；將幾何積分方法特別是顯式幾何積分方法引入電力系統(tǒng)數(shù)值仿真計算中，可以得到更好的穩(wěn)定性。

②顯式保積方法相對于傳統(tǒng)數(shù)值方法具備精確度更高、穩(wěn)定性更好、計算效率更快等優(yōu)勢。

③對采用詳細(xì)模型的復(fù)雜電力系統(tǒng)，在一般情況下，我們很難將其轉(zhuǎn)換為一個可分的Hamilton系統(tǒng)，但可以嘗試將其轉(zhuǎn)化為可分無源動力系統(tǒng)，然后采用本文所述的顯式保積方法進行計算。