基于高階Radau方法的高壓輸電線路電磁暫態快速計算方法
【技術領域】
[0001] 本發明涉及一種電力系統電磁暫態數值計算方法，具體是涉及一種基于多級高階 隱式Radau IIA方法及V-變換的高壓輸電線路電磁暫態快速數值計算方法。
【背景技術】
[0002] 高壓輸電線路是電力系統的重要元件之一。利用數值計算方法分析高壓輸電線路 的電磁暫態過程，可解決電力系統諸多領域的問題:可對輸電線路故障進行仿真，得到短路 電流以及諧波分量等;可用于輸電線路故障測距算法的研究以及輸電線路空載合閘過電壓 計算研究等。
[0003] 概括起來，高壓輸電線路的電磁暫態數值計算方法大致包括頻域法和時域法這2 大類方法。頻域類方法通常是先利用快速傅立葉變換或拉氏變換求得輸電線的頻率響應， 然后利用相應的反變換以得到其時域解。頻域類方法概念比較簡單，但涉及的計算極為繁 瑣，而且此類方法的求解過程主要是基于線性疊加原理，因而很難處理非線性問題。
[0004]時域法的基本思路是采用時間離散或空間離散方法將偏微分方程轉化為常微分 方程，并在此基礎上利用各類數值積分方法求解各狀態變量的值。目前使用較多的時域類 方法主要包括:時域有限差分法、基于隱式梯形積分方法和隱式歐拉法相結合的臨界阻尼 調整法（Critical Damping Adjustment，CDA)等。
[0005] 時域有限差分法的主要思路是將描述高壓輸電線路電磁暫態計算的電報方程在 空間和時間上同時進行離散，用差分方程代替一階偏微分方程，再通過求解線性方程組以 計算出各狀態變量的值。時域有限差分法算法簡單，但該方法中時間步長的選取受穩定性 條件的限制，當空間步長很小時亦須選取很小的時間步長，從而導致計算效率低下。
[0006] 目前，在商用化的電磁暫態計算程序EMTP(Electromagnetic Transient Program，EMTP)中主要是采用⑶A法來求解高壓輸電線路的電磁暫態計算問題。CDA方法計 算過程簡單，從理論上講可以避免傳統的電磁暫態計算方法中的數值振蕩問題，但對一些 無法準確判斷的突變現象，使用該方法進行計算仍將產生"虛假的"數值振蕩現象。
[0007] 迄今為止，在進行電力系統電磁暫態數值計算時，主要是采用低階的數值積分方 法，例如，應用最早的隱式梯形積分法是2階的數值方法，而CDA方法中所包含的隱式歐拉法 只是1階的數值方法。熟知:在采用低階的數值積分方法時，為保證計算精度通常只能采用 很小的計算步長，這自然導致極低的計算效率。為提高電磁暫態數值計算的計算效率，一個 可行的思路就是采用多級、高階的數值積分方法，因為在保持相同的計算精度的前提下，若 數值積分方法的階數愈高，則可選用的積分步長自然可以更大。然而，現有的高階數值積分 方法通常是多級的；對一個含有m個狀態變量的常微分初值問題，若采用s級的數值積分方 法對其進行求解，則在每一步的積分過程中需要求解一個s Xm維的線性或非線性方程組， 當m和s較大時，其巨大的計算量會"抵消"采用大步長所帶來的"好處"，有時甚至是"得不償 失"的。因此，在采用多級、高階數值積分方法以便采用較大的積分步長的同時，須解決高維 方程組的高效求解問題，只有合理解決這一問題才能真正提高電磁暫態計算過程的計算效 率。
【發明內容】
[0008] 本發明所要解決的技術問題，就是提供一種基于多級、高階隱式Radau IIA方法的 高壓輸電線路電磁暫態數值計算方法。由于Radau IIA方法既是線性L-穩定的也是非線性 B-穩定的，因而該方法可以有效避免數值振蕩問題;對高壓輸電線路電磁暫態計算所涉及 的非齊次常微分方程組的求解，該方法能利用Radau IIA方法的系數矩陣滿足V-變換這一 特性，將多級高階方法所導致的高維線性方程組的求解問題解耦并形成多級分塊遞推求解 方法，因而能通過采用大步長來提高高壓輸電線路電磁暫態數值計算的計算效率。
[0009] 解決上述技術問題，本發明采取如下的技術方案：
[00?0] -種基于高階Radau方法的高壓輸電線路電磁暫態快速計算方法，其特征在于包 括以下步驟:利用Radau IIA方法的系數矩陣滿足V-變換這一特性，將多級高階數值積分方 法所導致的高維線性方程組的求解問題解耦并形成多級分塊遞推求解方法。
[0011] 具體步驟包括：
[0012] 1)將高壓輸電線路進行空間離散化，建立相應的電磁暫態數值計算的數學模型， 即：
[0013]
(1);
[0014]式（1)中，X是高壓輸電線路電磁暫態計算中的狀態變量;B是與高壓輸電線路參數 相關的定常系數矩陣;w(t)即是電磁暫態計算的輸入信號或激勵信號；
[0015] 2)電磁暫態數值計算初始化
[0016] 置t = 0.0s，積分步數n = 0;確定Radau IIA方法的級數s(s 2 3);
[0017] 確定數值積分步長h、電磁暫態數值計算時程T;
[0018] 確定各狀態變量的初值，即X(t = 0)=XO;
[0019] 輸入電磁暫態數值計算的故障或操作；
[0020] 3)故障或操作判斷
[0021 ]依據時刻t判斷系統有無故障或操作：
[0022]若無故障或操作，則直接轉至步驟4);
[0023] 若有故障或操作，則修改系數矩陣B以及激勵信號w(t)并重新形成方程(1);
[0024] 同時，若涉及狀態變量的突變，則還需依據具體情況修改相應的狀態變量值Xn ⑴；
[0025] 4)數值積分
[0026]已知狀態變量奴〇在1時刻的值知，采用s級（s 2 3)的Radau IIA方法求解狀態變 量X(t)在tn+Ι時刻的值Xn+l;
[0027] 5)t = t+h；n = n+l
[0028] 6)數值計算是否終止判斷 [0029] 若t<T，則返回步驟3):
[0030] 若t 2 T，則轉至步驟7);
[0031] 7)數值計算結果輸出。
[0032]所述步驟4)數值積分是本發明的核心部分，其中所用到的s級的Radau IIA方法可 用以下Butcher表來描述：
[0033]
(2)；
[0034] 上式（2)中：A= [aij] ERsXs，bT= (bi …bs)，c= (ci …cs)T;其中，A即是s級的 Radau IIA方法的系數矩陣，當級數一定時，它是一個定常矩陣；Ci，i e (1，s)稱為Radau ΠΑ方法的內節點，它是由右Radau多項式
[0035]
(3)；
[0036]的s個零點組成。
[0037] 研究人員已證明：s級的Radau IIA方法滿足簡化階條件B(2s_l)，C(s)，D(s-Ι)，它 是A-穩定、L-穩定以及非線性B-穩定的2s-1階的數值積分方法；s級的Radau IIA方法滿足 簡化階條件C( s )，因此，其系數矩陣A滿足以下的V-變換：
[0038] A = VDV-1 (4)；
[0039] 上式(4)中：
[0040]
[0041] (6)；
[0042]事實上，隱式歐拉法即是單級（s = l)的Radau IIA方法。作為示例，下面給出幾個 多級的Radau IIA方法。
[0043] s = 3時，Radau IIA方法為：
[0044] :(〒):;
[0045] (8)；
[0046]
(9):;
[0047] s = 5時，Radau IIA方法為：
[0048] c = [0.0571 0.2768 0.5836 0.8602 1]τ (10)；
[0049] bT=[0.1437 0.2814 0.3118 0.2231 0.0400] (11)；
[0050]
(12)·?
[0051 ]上式（12)中，
[0052]
[0053] 所述的步驟4)數值積分也即對常微分方程初值問題(1)的求解步驟具體如下：
[0054] 第一步:將微分方程(1)在時間域離散，形成差分代數方程組
[0055] 利用數值積分方法(2)求解常微分方程(1)，可得：
[0056]
(13)；
[0057] 上式（13)中：h為積分步長;m為狀態變量X的維數;&是所有分量均為1的s維列向 量；
[0058] Im表示mXm維的單位矩陣；
[0059]
(14)?
[0060] 顯然，可以將方程(13)整理成以下形式的線性代數方程組：
[0061] JX=F (15)；
[0062] 上式（15)中，
[0063] (16)?
[0064] (17)；
[0065]第二步:利用V-變換對系數矩陣J進行分塊三角分解
[0066] Radau IIA方法的系數矩陣A滿足V-變換即方程(4)，因此，有：
[0067]
[0068]
[0069] (19)?
[0070]顯然，依據方程(6)，J的表達式為：
[0071]
(20)；
[0072] 對J進行分塊三角分解，可得：