陶靜靜 潘明帥 汪芳宗

(三峽大學(xué) 電氣與新能源學(xué)院,湖北 宜昌 443002)

微分動力系統(tǒng)是研究隨時間變化的動力系統(tǒng)的整體性質(zhì)及其在擾動中的變化.對微分動力系統(tǒng)的研究通?？梢詺w為一階微分方程問題的求解,微分動力系統(tǒng)中高維非線性問題歸根結(jié)底是微分初值問題(initial value problem,IVP).有關(guān)微分初值問題的研究,已經(jīng)獲得了豐富的研究成果.其中,具有重要影響且具代表意義的數(shù)值方法大致包括:Runge-Kutta方法(RK methods)[1],微分求積法(differential quadrature methods,DQM)[2-4],線性多步法(linear multistep formulae,LMF)[1],邊界值方法(boundary valuemethods,BVM)[5-6]等.

在實際工程應(yīng)用中,隱式梯形法因具有單步A-穩(wěn)定性和二階精確度而被廣泛使用.雖然隱式梯形算法對于逐步積分計算很容易,但仍然存在計算效率問題.由于該算法是單步和低階數(shù)值積分方法,因此在實際應(yīng)用中不能選擇過大的積分步長,從而將不僅會增加積分步數(shù),而且會損失由太多積分步數(shù)帶來的累積誤差,降低計算精度.為了克服單步數(shù)值法的不足,人們開始考慮多級高階整合算法,其中RK方法被廣泛研究.雖然RK方法具有良好的數(shù)值穩(wěn)定性,但是卻很少應(yīng)用到實踐當(dāng)中.這是因為當(dāng)微分動力系統(tǒng)的維數(shù)增大時,容易引起階障礙,也就是所謂的Dahlquist障礙(Dahlquist barriers)[7],這個問題同樣也出現(xiàn)在微分積分方法的求解中.對一個m維的非線性微分初值問題,s級的隱式RK方法在每一步的積分過程中需要求解s×m維非線性方程聯(lián)立的方程組;當(dāng)m很大時,其巨大的計算量讓人難以接受.Brugnano等學(xué)者為了克服LMF類方法所存在的穩(wěn)定性及階障礙,提出一系列新的方法,稱之為邊界值方法[5].Brugnano將BVM分為3類[6],分別為:擴展的梯形積分方法(extended trapezoidal rules,ETRs);廣義BDF方法(generalized backward differentiation formulae,GBDF);廣義Adams方法(generalized Adams methods,GAMs).概括起來,BVM類方法在保持LMF類方法所具有的多級高階優(yōu)點的同時,避免了高階LMF方法所存在的穩(wěn)定性問題,因而具有潛在的應(yīng)用價值.然而,BVM的現(xiàn)有研究主要針對線性微分動力學(xué)系統(tǒng)提出,限制了這類方法的應(yīng)用.

基于上述思路,本文將ETR2方法應(yīng)用于非線性微分動力系統(tǒng)數(shù)值計算,提出了一種新的數(shù)值計算方法.所提方法的主要優(yōu)點:基于邊界值方法的高階特性,該方法可以采用更大的積分步長,由此可以減少數(shù)值積分的步數(shù);基于整體雅可比矩陣所具有的特殊帶狀結(jié)構(gòu)特征,采用適當(dāng)?shù)木仃嚪匠谭至鸭记?避免了對高維的整體雅可比矩陣或多個分塊子矩陣進行三角分解,因而提高了計算效率.

1 ETR2方法簡介

以下述一階常微分初值問題為例:

基于均勻網(wǎng)格的s級的ETR2方法可用下述公式來描述:

式(2)中βi,i∈(0,μ-1)的取值可參見文獻[6];

s級的ETR2方法(2)是A-穩(wěn)定、s+1階的數(shù)值方法,其穩(wěn)定域與傳統(tǒng)的隱式梯形積分方法的穩(wěn)定域完全一致[6,8].作為示例,當(dāng)s=3時,上述方法(2)的具體表達式為:

與傳統(tǒng)的“逐步積分”求解過程不同,將BVM方法應(yīng)用于初值問題時,只能利用積分規(guī)則在多個時間點上連續(xù)進行離散,然后對離散后的方程進行整體求解.另外,顧名思義,將邊界值方法應(yīng)用微分方程的數(shù)值計算時還需結(jié)合相應(yīng)的邊界條件.若將該方法應(yīng)用于方程(1)的數(shù)值求解,則需要已知x0,x1,…,xv-1以及xN-s+v+1,…,xN等s個邊界值.然而,對微分初值問題而言,通常只有x0是已知的.因此,將邊界值方法應(yīng)用于初值問題時,必須將邊界值方法與所謂的附加方法(additional methods)[5-6]聯(lián)合起來,以便形成與待求點數(shù)N相匹配的N個方程.為此,Brugnano等學(xué)者提出了附加方法的選擇策略[5],具體可概括如下:附加方法與采用的邊界值方法應(yīng)相互獨立,但兩者應(yīng)具有相同的階數(shù).依據(jù)上述策略選擇的附加方法,將不會影響邊界值方法作為主體方法的數(shù)值穩(wěn)定性及計算精度.

2 基于ETR2方法的微分初值問題數(shù)值計算方法

以3級4階的ETR2方法(4)為例.

采用該方法時,需要利用附加方法形成2個附加方程.由于ETR2方法(4)是對稱的方法(symmetric method),因此,不能選擇該方法的反射方法作為附加方法,但可以選擇同階的、非對稱ETR2方法[9]作為附加方法.

以方法(4)為主體方法,選取以下方法作為時域始端的附加方法:

同時選取以下方法作為時域末端的附加方法:

利用方法(4)、(5)及(6)對方程(1)進行連續(xù)的時域差分離散,然后利用牛頓法對離散后的方程進行整體求解,可得

式(7)中:

上述表達式中:

下面介紹牛頓方程(7)的快速求解方法.依據(jù)JE的結(jié)構(gòu)特征,可以將方程(7)分解為以下3個方程組:

上述3個方程中:

在導(dǎo)出關(guān)系式(25)后,將其代入方程(20)和(21),最終可以導(dǎo)出一個關(guān)于的矩陣方程組.為方便起見,將這個矩陣方程組用式(26)表示:

顯然,在涉及N個網(wǎng)格點的積分求解過程中,上述求解方法只需對一個2m維的方陣進行三角分解.概括起來,對s級(s為奇數(shù))的ETR2方法,上述求解方法也只需對一個(s-1)×m維的方陣進行三角分解.

3 數(shù)值算例及分析

大規(guī)模電力系統(tǒng)是典型的高維非線性微分動力系統(tǒng),以采用經(jīng)典模型的電力系統(tǒng)暫態(tài)穩(wěn)定性數(shù)值計算為例.在此情況下,電力系統(tǒng)暫態(tài)穩(wěn)定性計算的數(shù)學(xué)模型可用以下微分方程組來描述[8]:

上式(27)中:i∈(1,p),p為發(fā)電機的數(shù)量;δi,ωi分別為發(fā)電機的功角與角頻率;Mi為發(fā)電機的慣性時間常數(shù);Pmi為發(fā)電機的機械功率;Cij、Dij均為常量.

隱式梯形法是屬于線性多步法的一種,因為其出色的數(shù)值穩(wěn)定性在電力系統(tǒng)領(lǐng)域被廣泛應(yīng)用.顯然,其每一步的求解可用式(28)來描述:

式(28)中:上標ki表示該步積分的迭代次數(shù).

很易理解:若用隱式梯形法積分N步,則該方法總共需要對N×ki個m維的矩陣Hi,i∈(1,N)進行三角分解,這里的ki即是隱式梯形法在第i步的積分中所需的迭代次數(shù);若在每一步的積分求解中采用定雅可比牛頓法或隱式梯形法(very dishonest Newton method,VDHN),則該方法總共需要對N個m維的矩陣Hi∈(1,N)進行三角分解.若采用本文所述的方法積分N步,則總共需要對η個(s-1)×m維的矩陣進行三角分解,這里的η即是該方法所需的迭代次數(shù).因此,從理論上講,若η＜＜N,則本文所提出的計算方法在計算效率上將比隱式梯形法具有優(yōu)勢.

首先選用IEEE145節(jié)點系統(tǒng)[10]作為算例系統(tǒng)(簡稱算例1).該系統(tǒng)含50臺發(fā)電機(p=50,m=100).暫態(tài)穩(wěn)定性計算中,故障設(shè)定為在7號母線處發(fā)生三相短路,經(jīng)0.1 s后切除.利用該算例系統(tǒng),將本文所提出的方法與隱式梯形法進行對比測試.對比測試中,對2種方法設(shè)定相同的收斂精度,即‖?R(X)‖≤ε=10-5,‖?r~(xi)‖≤ε=10-5.為便于對計算精度進行對比分析,利用中國電科院的PSASP程序同時采用0.001 s的小步長進行暫穩(wěn)計算,將所得的計算結(jié)果作為基準,以此來跟蹤、觀察不同算法的計算誤差曲線.

圖1是利用隱式梯形法(h=0.02 s)和3級4階ETR2方法(N=40,h=0.04 s)所計算出的幾臺發(fā)電機功角搖擺曲線,其中采用嚴格牛頓法的隱式梯形法每一步的積分需迭代2～3次,即ki=2～3,3級4階ETR2方法需迭代5次,即η=5,圖2是2種方法所得結(jié)果相對于基準結(jié)果的誤差曲線(故障切除后).

從圖2可以看出:3級4階ETR2方法在采用較大的步長的情況下,其計算精度仍優(yōu)于采用較小步長的隱式梯形法.

另選一個大規(guī)模電力系統(tǒng)作為算例系統(tǒng)(簡稱算例2).該系統(tǒng)含327臺發(fā)電機(p=327,m=654)、2 383個節(jié)點[11].暫態(tài)穩(wěn)定性計算中,故障設(shè)定為在171號母線處發(fā)生三相短路,經(jīng)0.1 s后切除.設(shè)定如上所述的、相同的收斂精度,將本文所提出的方法與隱式梯形法進行對比測試.

圖1 兩種方法所得發(fā)電機功角搖擺曲線(算例1)

圖2 兩種方法所得結(jié)果的誤差曲線(算例1)

圖3 兩種方法所得發(fā)電機功角搖擺曲線(算例2)

圖3是分別利用隱式梯形法(h=0.02 s)、3級4階ETR2方法(N=40,h=0.04 s)所計算出的幾臺發(fā)電機功角搖擺曲線,其中(采用嚴格牛頓法的)隱式梯形法每一步的積分需迭代2～3次,即ki=2～3,3級4階ETR2方法需迭代6次,即η=6.

圖4是2種方法所得結(jié)果相對于基準結(jié)果的誤差曲線(故障切除后).從圖4可以看出:3級4階ETR2方法在采用較大的步長的情況下,其計算精度仍略優(yōu)于采用較小步長的隱式梯形法.

圖4 兩種方法所得結(jié)果的誤差曲線(算例2)

表1是對兩種不同方法的計算速度的測試結(jié)果.

表1 兩種不同方法的計算速度測試結(jié)果

從表1可以看出:本文所提出的計算方法在計算效率上明顯優(yōu)于隱式梯形法.

4 結(jié) 論

本文將ETR2方法應(yīng)用于求解非線性微分動力系統(tǒng),提出了一種新的快速數(shù)值計算方法.該方法作為邊界值方法新的應(yīng)用方式,具有高階、數(shù)值穩(wěn)定性好等優(yōu)點.同時,本文所提方法是利用ETR2方法進行求解,在數(shù)值計算的過程中,充分利用ETR2方法本身的結(jié)構(gòu)特點,采用適當(dāng)?shù)木仃嚪纸庵貥?gòu)技巧,避免了迭代過程中對大型雅克比矩陣分解產(chǎn)生的大量計算,因此顯著地提高了數(shù)值計算效率.

理論分析和仿真測試結(jié)果表明:本文所提出的基于ETR2方法的快速數(shù)值計算方法,在計算效率上比基于隱式梯形積分規(guī)則的數(shù)值計算方法具有明顯的優(yōu)勢.

[1] Butcher J C.Numerical Methods for Ordinary Differential Equations,Second Edition[M].New York:John Wiley&Sons,2008.

[2] Bellman R,Casti J.Differential Quadrature and Longterm Integration[J].Journal of Mathematical Analysis and Applications,1971,34(2):235-238.

[3] 汪芳宗,廖小兵.基于微分求積法及V-變換的大規(guī)模動力系統(tǒng)快速數(shù)值計算方法[J].振動與沖擊,2016,35(3):73-78,128.

[4] 汪芳宗,廖小兵,謝 雄.微分求積分的特性及其改進[J].計算力學(xué)學(xué)報,2015,32(6):765-771.

[5] Brugnano L,Trigiante D.Solving Differential Problems by Multistep Initial and Boundary Value Methods[M].Amsterdam:Gordon and Breach,1998.

[6] Brugnano L,Trigiante D.Boundary Value Methods:the Third Way between Linear Multistep and Runge-Kutta Methods[J].Computers&Mathematics with Applications,1998,36(10):269-284.

[7] Dahlquist G.Convergence and Stability in the Numerical Integration of Ordinary Differential Equations[J].Journal of Mathematica Scandinavica,1956,4(4):33-53.

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

[9] Amodio P,Brugnano L.Parallel ODE Solvers Based on Block BVMs[J].Advances in Computational Mathematics,1997,7(1):5-26.

[10]Milano F.An Open Source Power System Analysis Toolbox[J].IEEE Trans.on Power Systems,2005,20(3):1199-1206.

[11]Zimmerman R D,Mueillo S X.MATPOWER:Steadystate Operations,Planning,and Analysis Tools for Power Systems Research and Education[J].IEEE Transactions on Power Systems,2011,26(1):12-19.