基于2級3階單對角隱式Runge-Kutta法的電磁暫態(tài)計算方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及一種電力系統(tǒng)電磁暫態(tài)數(shù)值計算方法，具體是涉及一種基于2級3階、 非線性代數(shù)穩(wěn)定的單對角隱式龍格-庫塔（Runge-Kutta，RK)方法的電磁暫態(tài)計算方法。
【背景技術(shù)】
[0002] 電磁暫態(tài)計算是研究電力系統(tǒng)各個元件中電場和磁場以及相應(yīng)的電壓和電流的 變化過程，其主要目的在于研究電力系統(tǒng)故障或操作過后可能出現(xiàn)的暫態(tài)過電壓和過電 流。
[0003] 在電磁暫態(tài)分析計算過程中，考慮到元件的非線性、電磁耦合、長線路波過程和線 路三相不對稱、線路參數(shù)的頻率特性等因素，需建立系統(tǒng)元件的微分方程或偏微分方程，并 借助一定的數(shù)值計算方法對這些方程進行離散化，得到代數(shù)形式的差分方程，進而求解各 時間點的物理量。
[0004] 目前，經(jīng)典的電磁暫態(tài)數(shù)值計算主要是采用隱式梯形積分方法。隱式梯形積分方 法具有2階精度以及A-穩(wěn)定性，但該方法不是L-穩(wěn)定的。在電磁暫態(tài)數(shù)值計算過程中， 當發(fā)生電感電流或電容電壓突變以及開關(guān)元件動作等情況時，使用隱式梯形法進行電磁 暫態(tài)數(shù)值計算時會產(chǎn)生一系列"虛假的"、持續(xù)的數(shù)值振蕩。為了解決此問題，加拿大學者 J.R.Marti和我國學者林集明等將隱式梯形法與具有強阻尼特性的隱式歐拉法相結(jié)合，提 出了臨界阻尼調(diào)整法（CriticalDampingAdjustment，CDA)，并將此方法運用到電磁暫態(tài) 仿真程序（ElectromagneticTransientProgram，EMTP)中。隱式歐拉法既是A-穩(wěn)定的也 是L-穩(wěn)定的數(shù)值方法，可以避免數(shù)值振蕩問題，但它只是1階的數(shù)值方法。在使用CDA方 法進行電磁暫態(tài)仿真時，依舊以隱式梯形法作為數(shù)值計算的主要方法，僅在突變發(fā)生時刻 利用隱式歐拉法來進行計算，以此避免數(shù)值振蕩情況的發(fā)生。CDA方法可以有效避免數(shù)值振 蕩問題，但前提條件是需要檢測出突變現(xiàn)象及其發(fā)生的時刻。在電磁暫態(tài)數(shù)值計算中，突變 現(xiàn)象主要包括開關(guān)元件的動作，電感電流和電容電壓的突變，以及非線性電感、電容元件的 運行方式由分段線性曲線轉(zhuǎn)折點的一邊躍變至另一邊等。在實際仿真過程中，突變的種類 眾多，對某些突變現(xiàn)象很難準確判斷出突變發(fā)生的時刻。例如，當傳輸線始端電壓或電流發(fā) 生突變時，很難準確地判定其末端電壓或電流發(fā)生突變的時刻，或者是控制系統(tǒng)中電壓源 和電流源因為限幅環(huán)節(jié)的影響發(fā)生難以檢測到的突變等。因此，對一些很難檢測到的突變 現(xiàn)象，CDA方法仍然無法避免數(shù)值振蕩問題，例如，當一個階躍突變信號從傳輸線始端傳輸 到末端時（如圖1所示），若采用CDA方法進行電磁暫態(tài)計算，因不能有效檢測末端的電壓 突變現(xiàn)象（如圖2 (a)所示），仍然會發(fā)生數(shù)值振蕩（如圖2 (b)所示）。
[0005] 為進一步解決CDA方法的問題，日本研究人員TakuNoda等將2級2階單對角隱 式RK方法應(yīng)用于電磁暫態(tài)數(shù)值計算。該方法可用Butcher表（ButcherTable)表示如下：
[0006]
[0007] 與隱式歐拉法一樣，上述方法（1)也是L-穩(wěn)定的數(shù)值方法。因此，該方法可以避 免數(shù)值振蕩問題。與CDA方法相比較，上述方法（1)無需對突變現(xiàn)象進行檢測或判斷，這是 該方法的主要優(yōu)點。然而，無論有無突變現(xiàn)象，方法（1)均須在每個時間步長內(nèi)計算兩個內(nèi) 點的變量值，其每一步的積分相當于隱式歐拉方法用更小的步長連續(xù)積分2步。因此，其計 算量約為隱式梯形法的2倍，也比CDA方法的計算效率更低。

【發(fā)明內(nèi)容】

[0008] 本發(fā)明所要解決的技術(shù)問題，就是提供一種基于2級3階單對角隱式Runge-Kutta 法的電磁暫態(tài)數(shù)值計算方法，其能在不降低經(jīng)典電磁暫態(tài)數(shù)值計算方法計算效率的前提 下，解決隱式梯形積分方法所存在的數(shù)值振蕩問題，且計算效率比CDA方法以及基于2級2 階單對角隱式Runge-Kutta方法的電磁暫態(tài)數(shù)值計算方法的效率更高。
[0009] 解決上述技術(shù)問題，本發(fā)明采取如下的技術(shù)方案：
[0010] -種基于2級3階單對角隱式Runge-Kutta法的電磁暫態(tài)數(shù)值分析方法，其特征 在于：通過建立電力系統(tǒng)電磁暫態(tài)數(shù)值計算的時域微分方程，采用B-穩(wěn)定的2級3階單對 角隱式Runge-Kutta法進行時域數(shù)值積分計算，逐步求解出各物理量隨時間的變化曲線；
[0011] 完整的步驟包括：
[0012] 1)輸入原始數(shù)據(jù)，建立各元件的微分方程，形成電磁暫態(tài)數(shù)值計算的基本數(shù)學模 型（土 = /(i，JC));
[0013] 2)電磁暫態(tài)數(shù)值計算初始化
[0014]置t= 0· 0s，積分步數(shù)η= 0;
[0015] 確定數(shù)值積分步長h、電磁暫態(tài)數(shù)值計算時程Τ;
[0016] 確定各狀態(tài)變量的初值，即x(t= 0) =X。；
[0017] 輸入電磁暫態(tài)數(shù)值計算的故障或操作；
[0018] 3)故障或操作判斷
[0019] 依據(jù)時刻t判斷系統(tǒng)有無故障或操作；
[0020] 若有故障或操作，則修改相應(yīng)的微分方程以及相應(yīng)的狀態(tài)變量值xn(t);
[0021] 4)數(shù)值積分
[0022] 采用B-穩(wěn)定的2級3階單對角隱式RK方法，計算出狀態(tài)變量在t=t+h處的值 -^-n+1?
[0023] 5)t=t+h；n=n+1 ；
[0024] 6)數(shù)值計算是否終止判斷
[0025]若t<T，則返回步驟3)，繼續(xù)下一步即下一時刻的數(shù)值計算；
[0026] 若t彡Τ，則轉(zhuǎn)至步驟7);
[0027] 7)數(shù)值計算結(jié)果輸出。
[0028] 所述的步驟4)數(shù)值積分中，用到的2級3階單對角隱式RK方法的Butcher表為：
[0029]
[0030] 與隱式梯形積分方法以及方法（1) 一樣，方法（2)也是A-穩(wěn)定的數(shù)值方法，但方 法（2)不是嚴格意義上的L-穩(wěn)定的數(shù)值方法，而是非線性B-穩(wěn)定亦即非線性代數(shù)穩(wěn)定的 數(shù)值方法；
[0031] 對常微分方程初值問題：
[0033] 所述的步驟4)數(shù)值積分也即對常微分方程初值問題（3)的求解步驟具體如下：
[0034] 從tjljtn+1時亥I」，已知狀態(tài)變量X⑴在t=t"時刻的值Xn，求解其在t=tn+1時 亥1J的值xn+1;令計算的時間步長為h=tn+1_tn;
[0035] 第一步：計算狀態(tài)變量在第一個內(nèi)點處的近似值：
[0037] 其中，充》奴可)是狀態(tài)變量在內(nèi)點?= 〖"+βΑ處的近似值；
[0038] 若f(t，X)是X的線性函數(shù)，則依據(jù)方程⑷可直接求解出萬，由此得/(?)的 值；
[0039] 若f(t，X)是X的非線性函數(shù)，則方程⑷的求解采用牛頓迭代法，同樣可求解出 萬及相應(yīng)的
[0040] 第二步：計算狀態(tài)變量在第二個內(nèi)點處的近似值：
[0042] 上式中，/(??)為已知量；?Χ(ξ)是狀態(tài)變量在內(nèi)點ξ:? 處的近似 值；
[0043] 同理，若f(t，x)是X的線性函數(shù)，則依據(jù)方程（5)直接求解出毛，由此得/(^?) 的值；若f(t，X)是X的非線性函數(shù)，則方程（5)的求解采用牛頓迭代法，同樣可求解出％及 相應(yīng)的/(1么）；
[0044] 第三步：計算狀態(tài)變量在t=tn+1時刻的值：
[0046] 本發(fā)明的理論基礎(chǔ)：
[0047] 眾所周知，L-穩(wěn)定性是線性穩(wěn)定性的范疇。從理論上講，對線性微分動力系統(tǒng)， L-穩(wěn)定的數(shù)值方法可以避免數(shù)值振蕩問題，這是隱式歐拉法以及方法（1)能夠避免數(shù)值振 蕩問題的主要數(shù)學機理。但對非線性微分動力系統(tǒng)，L-穩(wěn)定的數(shù)值方法并不一定能完全避 免數(shù)值振蕩問題。為此，研究人員已建立了非線性穩(wěn)定性分析的概念及相關(guān)理論體系。關(guān) 于非線性穩(wěn)定性，一個重要的結(jié)論就是所謂的B-穩(wěn)定性以及非線性代數(shù)穩(wěn)定性。研究人員 已證明：對一個非退化的Runge-Kutta方法，非線性代數(shù)穩(wěn)定性等價于B-穩(wěn)定性。利用非 線性代數(shù)穩(wěn)定性的定義，可以驗證并得出以下結(jié)論：
[0048] 2級3階單對角隱式RK方法（即方法⑵）的Μ矩陣的表達式為：

[0050] 上式（7)中，B=diag(b)。顯然，Μ矩陣的特征值為：λ丨=〇， 此，數(shù)值方法（2)的Μ矩陣是非負定的。由于有匕=132= 1/2 >0,因此，上述2級3階單 對角隱式RK方法（2)是非線性代數(shù)穩(wěn)定的，也是Β-穩(wěn)定的。
[0051] 眾所周知，對常微分方程初值問題（方程（3))，Β-穩(wěn)定的數(shù)值方法滿足單邊 Lipschitz條件，即：
[0052]<f(t1;Xj)-f(t2,x2),X!-x2> ^ 〇 (8)；
[0053] 上式⑶中，〈·，·>表示內(nèi)積。因此，B-穩(wěn)定的數(shù)值方法具有能量耗散性。從 物理概念上講，這里的能量耗散性也就是非線性阻尼特性。因此，在系統(tǒng)發(fā)生突變時，B-穩(wěn) 定的數(shù)值方法不會產(chǎn)生數(shù)值振蕩問題。換言之，在電磁暫態(tài)數(shù)值計算中，若系統(tǒng)發(fā)生突變現(xiàn) 象，2級3階單對角隱式RK方法（2)可以避免數(shù)值振蕩問題。這就是本發(fā)明的理論基礎(chǔ)。
[0054] 下面給出2級3階單對角隱式RK方法（2)不會產(chǎn)生數(shù)值振蕩的幾個具體實例。
[0055] 圖3a是一個基本的線性R-L串聯(lián)電路，圖3b是施加的電流源；其中，開關(guān)K在施 加的電流i(t)降至零時（t= 0. 01秒）突然打開。圖4a是利用隱式梯形法（計算步長h =0. 05ms)對該測試電路進行數(shù)值計算的結(jié)果（產(chǎn)生嚴重的數(shù)值振蕩）；圖4b是利用2級 3階單對角隱式RK方法（計算步長h= 0.lms)進行數(shù)值計算的結(jié)果。顯然，從圖4b可以 看出：當突變發(fā)生時，2級3階單對角隱式RK方法沒有產(chǎn)生數(shù)值振蕩。
[0056] 圖5是一個R-L串聯(lián)電路，其中，電感部分由一個線性電感外加一個飽和電抗器組 成；開關(guān)K在t= 0秒時突然合閘。圖6a是利用隱式梯形法（計算步長h= 0. 05ms)對該 測試電路進行數(shù)值計算的結(jié)果（產(chǎn)生數(shù)值振蕩）；圖6b是利用2級3階單對角隱式RK方 法（計算步長h= 0.lms)進行數(shù)值計算的結(jié)果。顯然，從圖6b可以看出：當突變發(fā)生時， 2級3階單對角隱式RK方法沒有產(chǎn)生數(shù)值振蕩。
[0057] 由于隱式梯形法以及方法（1)均是2階的數(shù)值方法，而2級3階單對角隱式RK方 法則是3階的數(shù)值方法。因此，在采用相同的計算步長的情況下，2級3階單對角隱式RK方 法的計算精度要比隱式梯形積分方法以及方法（1)更高。很易理解，在滿足相同的計算精 度的前提下，2級3階單對角隱式RK方法可以采用比隱式梯形法更大的步長；在2級3階 單對角隱式RK方法采用2倍于隱式梯形積分法的步長的情況下，若兩者的計算精度大致相 同，則它們的計算效率大致相當。為此，以圖7所示的基本測試電路（可以獲得該測試電路 的精確解即解析解）為例，對2級3階單對角隱式RK方法相對于隱式梯形法的計算效率進 行了測試和評估。圖8是分別利用隱式梯形法（步長h= 0. 01ms)和2級3階單對角隱式 RK方法（步長h= 0.02ms)對圖7所示測試電路進行數(shù)值計算的誤差對比曲線。顯然，從 圖8可以看出：