基于RKNd算法的暫態穩定性快速數值計算方法
2012-04-29 00:00:00張磊汪芳宗胡佳怡
計算技術與自動化 2012年3期
摘要：RKNd方法是一類新的數值積分方法。在相同級數條件下，RKNd方法可達到的最高代數階比傳統的RungeKutta方法以及RungeKuttNystrm方法均高，而且具有更高的計算效率。將RKNd方法引入電力系統暫態穩定性數值計算。以IEEE145節點電力系統為例，通過數值實驗將新方法與電力系統分析中常用的傳統數值計算方法進行了對比分析。數值實驗結果表明，RKNd方法在計算精度和計算效率等方面均具有明顯的優勢，因而更適合于電力系統暫態穩定性及相似問題的數值計算。
關鍵詞：暫態穩定性；數值積分方法；RungeKutta方法；RungeKuttNystrm方法；RKNd方法
中圖分類號：TM744 文獻標識碼：A
Fast Numerical Simulation of Power System Transient Stability by RKNd Methods
ZHANG Lei，WANG Fangzong，HU jiayi
（Electrical Engineering Renewable Energy School，China Three Gorges University，Yichang443002，China）
Abstract：The RKNd method is a new kind of numerical integration methods， of which the order is higher than that of the traditional RungeKutta methods and RungeKuttaNystrm methods for the same stage. In this paper， the RKNd method is introduced to the numerical simulation of power system transient stability， and then a fast numerical simulation method has been proposed. The proposed method has been compared to both the traditional numerical integration method and the symplectic Gauss method using IEEE145bus power system， and the tested results show that the implicit RKNd method has the advantages both in calculation accuracy and in computational efficiency respectively over the symplectic Gauss method and the implicit trapezoidal rule. Therefore the proposed methods should be more suitable to numerical analysis of transient stability and other likewise problems.
Key words：transient stability；numerical integration method；RungeKutta method；RungeKuttNystrm method；RKNd method
1引言
數值積分方法是電力系統暫態穩定性分析計算的基本方法。最常用的數值積分方法大致包括隱式梯形積分法以及RungeKutta方法（RK方法），前者是隱式積分類方法，后者是顯式積分類方法。
近年來，研究人員又提出了不少新的數值積分算法。文獻[1]和文獻[2]將辛Runge—Kutta算法（辛RK方法）、文獻[3]將可分Hamiltonian系統的顯辛算法、文獻[4]將辛代數動力學算法用于暫態穩定性的計算，并對這幾種新的數值積分方法進行了測試和對比分析。文獻[5]和文獻[6]分別將多級高階辛RK算法以及多級高階辛RungeKuttaNystrm算法用于暫態穩定性的并行計算。
文獻[7]利用一階常微分方程導出的二階方程，借鑒Nystrm方法，提出了一類新的數值積分方法，即RKNd方法。RKNd方法的最大優點是：在相同級數情況下，RKNd方法可達到的最高代數階比傳統的RK方法高。在傳統的RK系列方法中，s級的顯式RK方法可達到的最高階數是s階；s級的隱式RK方法可達到的最高階數是2s階。但RKNd方法不同，2級的顯式RKNd方法可以達到4階；2級的隱式RKNd方法可以達到5階。因此，與同級的RK方法相比，RKNd方法具有更高的計算精度；與同階或略低階的RK方法相比，RKNd方法具有更高的計算效率。
本文將RKNd方法引入電力系統暫態穩定性的數值計算。以IEEE 145節點系統為例，分別將2級4階顯式RKNd方法與傳統的4級4階顯式RK方法、2級5階隱式RKNd方法與2級4階隱式辛RK方法進行了對比測試。測試結果驗證了RKNd方法在計算效率方面具有明顯的優勢，因而可以推廣應用于電力系統暫態穩定性及其它領域的數值計算。
2RKNd方法簡介
RKNd方法既不同于傳統的RK方法，也與RKN方法有所不同。對給定的2階常微分方程初值問題
=f（t，x）=g（t，x），x（t0）=x0 （1）
其s級的RKNd方法的一般形式為
yi=xn+cihn+h2∑sj=1aijg（yj），i∈（1，s）
xn+1=xn+hn+h2∑si=1big（yi）
n+1=f（tn+h，xn+1） （2）
計算技術與自動化2012年9月
第31卷第3期張磊等：基于RKNd算法的暫態穩定性快速數值計算方法
上式中，h為積分步長，yi（i∈（1，s））為真解x（t）在ti=tn+cih時的逼近，通常被稱為方法的內級，即yi≈x（tn+cih）。該RKNd方法也可用Butcher表
cAbT
表示。其中：
b=（b1，b2，…，bs）T∈Rs
A=（aij）∈Rs×s，c=（c1，c2，…，cs）T∈Rs
傳統RK方法的內級yi僅在常數項和h項與x（tn+cih）一致；RKNd方法的內級可以精確到h2項。這是RKNd方法與傳統RK方法的主要不同之處，也是其優點所在。從上面的計算公式也可以看出，RKNd方法與RKN方法也有較大的不同。
RKNd方法亦是一個系列方法。迄今為止，導出的比較有代表性的RKNd方法主要是2級4階顯式方法和2級5階隱式方法。
2.12級4階顯式RKNd方法
2級4階顯式RKNd方法的Butcher表為下式（3），這個計算格式很簡單，但它是4階方法。同級的傳統顯式RK方法只能達到2階。只有4級的顯式RK方法才能達到4階，但4級的顯式RK方法的計算量差不多是方法（3）的2倍。
0001/21/801/61/3（3）
2.22級5階隱式RKNd方法
2級5階隱式RKNd方法的Butcher表為
4—61011—4610004+6102+36507—261009+6369—636（4）
這個計算格式是5階的。同級的傳統隱式RK方法最多只能達到4階，而且，2級4階隱式RK方法（即2級Gauss方法）是全隱的（其系數矩陣A是一個滿陣）。
方法（4）的系數矩陣A是一個下三角陣，這意味著2級之間的求解可以解耦。因此，2級5階隱式RKNd方法與同級的傳統隱式RK方法相比，不僅計算精度會更高，而且其計算效率也更高。
3暫態穩定性快速數值計算方法
RKNd方法只能用于求解式（1）或形如式（1）的微分動力系統。電力系統暫態穩定性數值計算通常是采用1次常微分方程組（或微分—代數方程組）這種形式來描述的。但是，無論采用經典模型還是復雜模型，其右端函數（即f（t，x））均是連續可導的。因此，將暫態穩定性數值計算的傳統模型轉換成形如式（1）的微分動力系統，是一件很容易的事情。換言之，RKNd方法完全適用于采用任何模型的電力系統暫態穩定性的數值計算。
為簡單起見，同時不失一般性，本文采用電力系統暫態穩定性計算的經典模型。在此情況下，描述電力系統暫態穩定性的微分方程為
i=ωiMii=Pmi—Pei，i∈（1，m）（5）
m為發電機臺數；Mi、Pmi、Pei分別為發電機的慣性時間常數、輸入機械功率、輸出電磁功率，δ為發電機功角，其中
Pei=∑m（Cijsin δij+Dijcos δij）（6）
很易理解，上述暫態穩定性計算的經典模型可以轉換為形如式（1）的微分動力系統，即
i=αi（δ）i=βi（δ，ω），i∈（1，m） （7）
式中，
αi（δ）=（Pmi—Pei）/Mi（8）
βi（δ，ω）=γi（δ，ω）/Mi（9）
γi（δ，ω）=∑m[ωij（Dijsin δij—Cijcos δij）]（10）
3.1基于顯式RKNd方法的暫態穩定性計算
基于顯式積分方式的電力系統暫態穩定性計算比較簡單。將公式（2）和表達式（3）應用于2階微分方程（7）和1階微分方程（5），可以導出基于2級4階顯式RKNd方法的暫態穩定性計算公式如下：
i=δi（tn）+12hωi（tn）+18h2αi（δ（tn）） （11）
i=ωi（tn）+12hαi（δ（tn））+
18h2βi（δ（tn），ω（tn）） （12）
δi（tn+1）=δi（tn）+hωi（tn）+
16h2αi（δ（tn））+13h2αi（） （13）
ωi（tn+1）=ωi（tn）+hαi（δ（tn））+
16h2βi（δ（tn），ω（tn））+13h2βi（，） （14）
式（11）和（12）可以看做是對狀態變量的預測過程；式（13）和（14）則是對狀態變量的校正過程。上述求解過程的主要運算量是αi（δ（tn））、αi（）、βi（δ（tn），ω（tn））、βi（，）等4項的求解計算。
很容易理解：上述求解過程的計算量與傳統4級4階顯式RK方法大致相當；與3級4階顯式RKN方法相比，上述求解過程的計算量似乎略大。將1階常微分方程組轉換為2階常微分方程組，是導致2級4階RKNd方法計算量倍增的主要原因。
3.2基于隱式RKNd方法的暫態穩定性計算
隱式RKNd法要將公式（2）和表達式（4）應用于2階微分方程（7）和1階微分方程（5），首先定義：
δ=（δ1，δ2，…，δm）T，ω=（ω1，ω2，…，ωm）T
y=[δΤ，ωΤ]Τ∈R2m×1
yi=y（tn+cih）∈R2m×1，i∈（1，2）
F（yi）=[α1（δ），…，αm（δ），β1（δ，ω），…，
βm（δ，ω）]Τ，F（yi）∈R2m×1
G（y（tn））=[ω1（tn），…，ωm（tn），α1（δ（tn）），…，αm（δ（tn））]Τ
G（y（tn））∈R2m×1
φi=y（tn）+cihG（y（tn））∈R2m×1，
ψ=[φΤ1，φΤ2]ΤY=（yΤ1，yΤ2）Τ，
Z=[FΤ（y1），FΤ（y2）]Τ
υ=y（tn）+hG（y（tn））∈R2m×1
此時公式（5）和（7）可寫為以下形式
=G（y）
=F（y）（15）
利用s=2級的隱式RKNd方法對式（15）進行差分求解后的公式為：
Y=ψ+h2（AI）Z（16）
y（tn+1）=υ+h（BI）Z （17）
（tn+1）=G（y（tn+1））（18）
上述方程中，I為m×m維單位矩陣，表示矩陣的Kronecker積（也稱直積或張量積）。B=diag（bi），i∈（1，2），很明顯，要求解公式（17）和（18）需要先求解公式（16），而式（16）的求解需要同時在s=2個時間點上計算yi，i∈（1，2）。
定義
Q=A—1=qij
方程（16）可以變換成：
（AI）—1Y=（AI）—1ψ+h2Z（19）
即
（QI）Y=（QI）ψ+h2Z （20）
采用牛頓法進行求解（20）可以導出：
JΔY=RYΔY=ΔyT1，ΔyT2T （21）
其中R（Y）為牛頓殘差向量，J為雅克比矩陣。由于J的維數s倍于傳統逐步串行計算時雅克比矩陣（即下面定義的gi的維數，因此，若直接對式（21）進行求解，在求解過程中將會耗費很多時間。為了解決這個問題，定義：
gi=F（yi）yΤiH=diag（gi），i∈（1，2），
R（Y）=[rΤ1，rΤ2]Τ
因為q12=0，使得系數矩陣A是一個下三角陣。則J可寫成分塊矩陣的形式，其表達式為：
Jij=[Jij]，Jij∈Rm×m，i，j∈（1，2）J12=0，J21=q21IJii=qiiI—gi
由雅克比矩陣分塊后的的特殊形式，我們可以最終導出：
J11Δy1=r1J22Δy2=r2—J21Δy1 （22）
至此我們得到了基于隱式2級5階RKNd方法的暫態穩定快速計算方法。
將式（16）轉換為式（20）并且利用系數矩陣A的特性實現矩陣分裂是算法的關鍵。由（22）可看出此種方法實現了在2級求解之間的“解耦”，與顯式方法相比這種方法將會有更好的計算精度和速度，是一種快速的數值計算方法。
4算例測試結果及對比分析
顯示RKNd方法在計算速度和計算量上并沒有多大的優勢，對其分析對比意義不大。本算例對隱式RKNd方法與傳統方法進行了比較，算例系統采用IEEE145節點電力系統。故障設定為7號母線發生三相接地短路，故障持續時間為Δtf=0.10s，整體積分時間為1.50秒，此時系統處于穩定狀態。K是牛頓求解過程中所需的最大迭代次數。收斂精度統一設定為10—4。
為便于對比分析，以h=0.001s時用隱式梯形積分法求解的發電機相對功角δ12=δ1—δ2的計算結果為基準（此時步長很小，不同算法的計算結果基本上是一樣的），跟蹤觀察不同算法的誤差曲線。
表1是算法性能測試結果，在相同的計算精度下，RKNd方法的計算速度高于隱式梯形積分法，其積分一步相當于隱式梯形積分方法連續積分3步，計算效率明顯提高。在相同的積分步長下，2級5階 RKNd方法和隱式辛RK方法有相同的迭代次數，且與同級的隱式辛RK方法的CPU運行時間相差很小。
表1各方法性能測試結果
方法
積分步長
迭代次數
CPU時間（S）
隱式梯形法
0.05
2
0.5156
隱式辛RK方法
0.15
3
0.2031
隱式RKNDd方法
0.15
3
0.2188
圖1、圖2和圖3分別是在不同積分步長情況下，隱式梯形積分法的誤差曲線與隱式RKNd方法的誤差曲線。從圖1、圖2、圖3可以看出：當RKNd方法的步長取到0.25s的時候其計算精度還明顯高于傳統的隱式梯形法。
time/s
圖1隱式RKNd方法（h=0.15 s）與
隱式梯形積分法（h=0.05 s）誤差曲線比較
time/s
圖2隱式RKNd法（h=0.20 s）與
隱式梯形積分法（h=0.05）誤差曲線比較
time/s
圖3隱式RKNd法（h=0.25s）與
隱式梯形積分法（h=0.05）誤差曲線比較
圖4是隱式辛RK方法與RKNd方法的誤差曲線，在相同積分步長且同階情況下，RKNd方法能達到5階，而辛RK方法只能達到4階。RKNd的計算精度比同階的隱式辛RungeKutta方法更高。5結束語
本文將RKNd方法引入電力系統暫態穩定性數值計算領域，以IEEE145節點電力系統為例，將新方法與傳統的算法進行了對比分析。根據初步的數值實驗結果，可以得出以下幾個結論：
time/s
圖4隱式RKNd法（h=0.15 s）與
隱式RK法（h=0.15）誤差曲線比較
1）顯式RKNd方法與傳統4級4階RK方法相比，在計算量上相當，在計算速度上沒有多大優勢，是一種可供選擇的顯示方法。
2）與同級的隱式RK方法相比，隱式RKNd方法具有更高的計算精度；與同階或略低階的RK方法相比，RKNd方法具有更高的計算效率。
3）2級5階隱式RKNd方法比傳統的隱式梯形積分法精度更好，且計算效率更高，在采用大步長同時又對計算精度要求較高的情況下，此方法是一種較好的新方法。
參考文獻
[1]Sun Geng.Construction of high order symplectic Runge—Kutta methods[J].Journal of Computational Mathematics，1993，11（3）：250—260.
[2]Sun Geng.A simple way constructing symplectic Runge—Kutta methods[J].Journal of Computational Mathematics，2000，18（1）：61—68.
[3]馮康，秦孟兆.哈密爾頓系統的辛幾何算法[M].杭州：浙江科學技術出版社，2003： 185—239.
[4]FENG Kang，QIN Meng—zhao.Symplectic Geometric Algorithm for Hamiltonian Systems[M].Hangzhou：Zhejiang Science Technology Press，2003： 185—239.
[5]Wang S J，Zhang H.Symplectic Algebraic Dynamics Algorithm[J].Science of China（Series G），2007，50（2）：133—143.
[6]汪芳宗，何一帆.基于多級高階辛Runge—Kutta方法的暫態穩定性并行計算方法[J].電力系統保護與控制，2011，4（35）：40—46.
[7]WANG Fang—zong，HE Yi—fan.Parallel computation oftransient stability by symplectic Runge—Kutta method[J].Power System Protection and Control，2011，4（35）：40—46.
[8]汪芳宗，何一帆.基于辛龍格—庫塔—奈斯通方法的電力系統暫態穩定性并行計算方法[J].電網技術，2011，11（39）：22—28.
[9]WANG Fangzong，HE Yi—fan.A parallel computation method for power system transient stability based on symplecit RungeKuttNystrm method[J].Power System Technology，2011，11（39）：22—28.
[10]陳丙振，游雄.求解初值問題的RKNd方法[J].計算數學，2010，4（10），399—413.
[11]Chen Bingzhen，Xiong You.RKNd Methods for Solving initial Value Problems[J]， Mathematica Numerical Sinica，2010，4（10），399—413.
[12]LIU J.The Multifrontal Method for Sparse Matrix Solution：Theory and Practice[J].SIAM Review，1992，34（1），82—109.
[13]李壽佛.剛性微分方程算法理論[M].長沙：湖南科學技術出版社，1997.
[14]倪以信，陳壽孫，張寶霖.動態電力系統的理論和分析[M].北京：清華大學出版社，2002.
