李博文，張 磊，葉宇國，張 靜

（1.三峽大學(xué) 電氣與新能源學(xué)院，湖北 宜昌 443002；2.國網(wǎng)湖北省電力公司棗陽供電公司，湖北 棗陽 441200）

隨著電力系統(tǒng)的規(guī)模越來越大以及越來越多的電力電子設(shè)備應(yīng)用于電力系統(tǒng)中，電網(wǎng)的動態(tài)特性變得越來越復(fù)雜，對于電力系統(tǒng)仿真實(shí)時性的要求也變得越來越高［1-2］。目前已有許多商業(yè)化的電磁暫態(tài)仿真軟件。其中，基于EMTP類算法的電磁暫態(tài)程序目前已經(jīng)成為設(shè)計電力設(shè)備和電力系統(tǒng)的一個必備工具，在技術(shù)研究和程序開發(fā)中也扮演著越來越重要的角色。Power System Toolbox，PSCAD/EMTDC，Transient Performance Advisor（TPA），EPRI/DCG等眾多仿真軟件都是基于EMTP所開發(fā)［3］。

電磁暫態(tài)數(shù)值計算所采用的時間尺度往往非常小，通常是微秒級［4］。對于特快速暫態(tài)過電壓（very fast transient overvoltage，VFTO）的計算，所采用的積分步長會更小，一般是納秒級［5］。因此即使仿真時間很短，其對應(yīng)的計算量依舊很大，這是造成現(xiàn)有電磁暫態(tài)仿真程序?qū)崟r性不足的原因之一。同時現(xiàn)有的大多數(shù)商業(yè)電磁暫態(tài)仿真軟件，往往是基于傳統(tǒng)的串行計算模式，因此不可避免的面臨計算效率低的問題，這也是造成仿真實(shí)時性不足的原因之一。

為了滿足對電力系統(tǒng)仿真實(shí)時性的需求，就必須提高電磁暫態(tài)數(shù)值計算的計算速度，對此有2個途徑，一個是采用新的數(shù)值算法［6-7］，另外一個是并行計算技術(shù)。其中并行計算技術(shù)可以通過2種手段實(shí)現(xiàn)：一是空間并行，考慮通過長輸電線將電力系統(tǒng)分割得到完全解耦的子系統(tǒng)，各子系統(tǒng)可以獨(dú)立并行計算［8-9］。但是由于各個子系統(tǒng)的時間常數(shù)不同，若根據(jù)時間常數(shù)最小的子系統(tǒng)確定積分步長，則會嚴(yán)重浪費(fèi)計算資源，增加仿真時間，為了解決這個問題，文獻(xiàn)［10］結(jié)合多速率電磁暫態(tài)仿真技術(shù)，根據(jù)不同的子系統(tǒng)確定不同的積分步長，優(yōu)化了計算資源，提高了電磁暫態(tài)仿真效率。二是時間并行，考慮將數(shù)值算法并行化，文獻(xiàn)［11］在臨界阻尼調(diào)整法的基礎(chǔ)上，利用矩陣對角化技術(shù)，可實(shí)現(xiàn)精確的時間解耦并行計算。文獻(xiàn)［12］基于不動點(diǎn)迭代，推導(dǎo)出了一套收斂性更好的電磁暫態(tài)時間并行算法。

本文從并行化數(shù)值算法的角度出發(fā)，基于二步二階邊界值方法，對一階常微分方程進(jìn)行離散，通過牛頓法對離散后的方程進(jìn)行求解。離散后的雅克比矩陣是塊三對角矩陣，對于這種形式的矩陣，常用的方法是追趕法［13］，追趕法是一種快速的串行方法，本質(zhì)是高斯消元，無法實(shí)現(xiàn)并行化。本文將初始雅克比矩陣重新排列為等價矩陣，通過一個耦合矩陣，將大矩陣分解為各自獨(dú)立的小矩陣，推導(dǎo)出了二步二階邊界值并行算法（boundary value parallel methods，BVPM）。為了驗證本文算法性能，最后通過2個典型的非線性電磁暫態(tài)算例進(jìn)行數(shù)值測試，將本文算法與追趕法和EMTP類方法進(jìn)行比較。

1 二步二階邊界值法簡介

邊界值方法是L.Brugnano等學(xué)者為了克服傳統(tǒng)線性多步法存在的穩(wěn)定性以及階障礙所提出的一大類方法。按照L.Brugnano的思路，邊界值方法主要分為3類［14］：廣義向后差分法、廣義Adams以及拓展的梯形積分公式。這3類方法均突破了Dahlquist限制性定理，它們都是A-穩(wěn)定的數(shù)值方法。因此在數(shù)值穩(wěn)定性方面均要優(yōu)于傳統(tǒng)的線性多步法。除此之外，L.Lopez還提出了一種特殊的二步邊界方法，并研究了其穩(wěn)定性和收斂精度［15］。

考慮一階常微分方程初值問題為：

式中：t表示時間；f（t，x（t））表示時間t和變量x的函數(shù)；x0表示變量x在t=0時的初始值。

本文采用帶獨(dú)立參數(shù)的二步二階邊界值方法離散式為：

式中：h=tn+1-tn=T／N，T表示仿真時長，N表示時間離散網(wǎng)格點(diǎn)數(shù)。

文獻(xiàn)［15］已經(jīng)證明，當(dāng)θ＜-1時，具有（1，1）邊界條件的二步二階邊界值方法是一種A-穩(wěn)定的數(shù)值方法，穩(wěn)定域如圖1所示。圖1中邊界曲線外部為二步二階邊界值法的穩(wěn)定域。

圖1 二步二階邊界值法穩(wěn)定域示意圖

顯然，對于微分初值問題，通常只有x0是已知的，而xN是未知的，因此還需要添加一個附加方程作為邊界條件，本文選用隱式梯形法為時域末端的附加方程，即：

2 基于二步二階邊界值的并行算法

2.1 并行實(shí)現(xiàn)

利用二步二階BVM法作為主方法，隱式梯形積分方法作為附加方法對式（1）進(jìn)行連續(xù)的時域差分離散，然后再利用牛頓法對離散后的方程進(jìn)行整體求解，得到的雅克比矩陣方程為：

式中：J為整體雅克比矩陣：▽F表示牛頓殘差向量。

式中，Im為m維單位矩陣。

通過式（6）可以看出，由于二步二階BVM是對整個積分區(qū)間一次性求解出所有時間節(jié)點(diǎn)上的點(diǎn)，從而導(dǎo)致式（5）的維數(shù)會非常大，存在“維數(shù)災(zāi)”問題。

為了解決這一問題，常采用直接或迭代的方法求解。具體的，對于式（5）形式的分塊三對角矩陣，往往采用追趕法進(jìn)行求解［13］，但是追趕法本質(zhì)上是一種串行方法。本文針對式（5）的結(jié)構(gòu)特性，通過交換行塊和列塊，對雅克比矩陣進(jìn)行重排，重排后的方程記為：

式中：u=［Δx1，Δx2，Δx4，Δx5，…，ΔxN］T；v=［Δx3，Δx6，Δx9，…，ΔxN-2］T；

D=diag［P3，P6，…，PN-2］。A，B，C，D矩陣中，M為子系統(tǒng)數(shù)。

初始矩陣被簡化為一些新的對角獨(dú)立塊矩陣、2個輔助塊矩陣以及1個耦合塊矩陣。根據(jù)式（7），可以得到：

為了避免對矩陣A求逆，于是令：

基于上述推導(dǎo)，將二步二階BVM并行算法概述如下：

1）根據(jù)需要確定需要離散的時間點(diǎn)個數(shù)N以及時間積分步長h。

2）利用二步二階邊界值方法作為主方法，隱式梯形法作為時域末端附加方法，對微分方程進(jìn)行離散。

3）利用牛頓法對離散后的方程整體求解，確定雅克比矩陣，將大系統(tǒng)雅克比矩陣分解為M組獨(dú)立的小系統(tǒng)矩陣，由矩陣D耦合到大系統(tǒng)中。

4）根據(jù)式（11）計算出矩陣v，再根據(jù)式（12）并行計算出矩陣u。

5）收斂性判斷，若收斂，則計算結(jié)束，若不收斂，則轉(zhuǎn)到新一輪牛頓迭代中，直至收斂。

顯然，本文算法需要用到三角分解的有式（9）和式（11），但是由于該方法的并行結(jié)構(gòu)以及新塊矩陣的維數(shù)遠(yuǎn)遠(yuǎn)小于原矩陣的維數(shù)，因此本文算法理論上會比直接利用追趕法求解式（4）方法具有更高的計算效率。

2.2 加速比分析

為了評價本文算法的加速效果，將本文算法與追趕法的計算速度進(jìn)行對比分析。定義加速比為：

式中：λ表示加速比；tBVM、tBVPM分別表示追趕法和本文所提并行算法的計算時間。

通過考慮塊矩陣的矩陣運(yùn)算次數(shù)來估計算法的計算速度。由于矩陣的加法運(yùn)算所消耗的時間遠(yuǎn)比乘法運(yùn)算小，因此忽略加法運(yùn)算，只考慮乘法運(yùn)算。根據(jù)文獻(xiàn)［16］，假設(shè)n×n的塊矩陣做一次矩陣乘法和求逆運(yùn)算均需要n3次乘法運(yùn)算。

追趕法本質(zhì)還是高斯消元法，利用雅克比帶狀稀疏矩陣的非零系數(shù)分布特點(diǎn)，來提高運(yùn)算速度，但它是一種串行方法，不具備并行化特點(diǎn)。追趕法求解式（4）所需要的矩陣運(yùn)算次數(shù)為：

對于本文算法，使分解后的每一個子系統(tǒng)有著相同的大小，保證每一個計算單元的計算量大致相同，避免了因為等待某些計算量大的計算單元而延長求解時間。設(shè)Nk表示子矩陣的對角塊數(shù)（k∈［1，M］），并且有：

得到二步二階BVM并行算法的塊矩陣運(yùn)算次數(shù)為：

假設(shè)參與并行計算的計算單元的數(shù)量為P，各個子系統(tǒng)大小相等，且子系統(tǒng)數(shù)M=iP，每個計算單元都按照各自所對應(yīng)子系統(tǒng)隊列，順序執(zhí)行操作。我們可以得到加速比的近似值為：

圖2給出了在不同i值情況下，并行度（計算單元數(shù)量P）與理論加速比λ關(guān)系的曲線。

圖2 理論加速比

通過圖2可以發(fā)現(xiàn)，當(dāng)計算單元數(shù)量較小時，加速比與計算單元數(shù)量呈近似線性增長的關(guān)系，隨著計算單元的增加，本文算法非并行的部分會逐漸增大，加速比的增長速度會慢慢變緩，理論加速比會在達(dá)到最大值后開始降低。同時我們可以發(fā)現(xiàn)，隨著i值的增加，加速比的最大值變得越來越小。當(dāng)計算單元較少時，i值對于加速比的影響并不大，然而隨著計算單元的增加，i值越大，本文算法非并行計算部分增加得就越快，導(dǎo)致加速比的最大值急劇下降。因此，為了充分發(fā)揮本文算法的性能，除了要滿足各子系統(tǒng)大小要相同，還必須要滿足子系統(tǒng)數(shù)M等于計算單元的數(shù)量。計算單元數(shù)量通過圖2所示所能達(dá)到的最大加速比時的并行度來確定。

3 算例分析

3.1 算例1含非線性負(fù)荷的高壓輸電線路

圖3為一個含非線性負(fù)荷的高壓輸電線路示意圖。其中es（t）為激勵電壓源，開關(guān)在t=0 s時合閘。輸電線路送端內(nèi)電阻Rs=20Ω、內(nèi)電感Ls=0.06H。輸電線路受端的非線性負(fù)荷是由非線性電感和負(fù)荷電阻并聯(lián)組成。負(fù)荷電阻RL=96Ω，非線性電感的磁鏈φ與線路電流iL滿足φ=a tanh biL，其中a=8.4×10-2，b=5.95×10-3。

描述輸電線路的電磁暫態(tài)仿真過程，一般采用著名的電報方程。電報方程是一種雙曲偏微分方程，因此需要將其轉(zhuǎn)換為式（1）的常微分方程?？紤]到高壓輸電線路可以忽略電導(dǎo)，因此本文采用π型級聯(lián)的等效電路對輸電線路進(jìn)行建模［17］。

圖3 含非線性負(fù)荷高壓輸電線路

輸電線路分為30段，根據(jù)基爾霍夫電壓、電流定律，可以列寫出上述空間離散的π型等效電路的電路方程為：

式中，l、c、r分別表示離散后的每一段線路的電感、電容以及電阻。其值依次為：

當(dāng)開關(guān)閉合時，應(yīng)該滿足如下邊界條件：

為了驗證本文算法的計算性能，將本文算法與經(jīng)典的EMTP算法的計算結(jié)果進(jìn)行了對比。本文算法獨(dú)立參數(shù)θ=-2，時間積分步長均取h=10-6s，牛頓迭代的收斂精度設(shè)置為ε=10-4。圖4和圖5分別為2種方法計算的送端、受端電壓曲線。從圖4、圖5中可以看出，在相同的積分步長情況下，本文算法有著與EMTP類方法相當(dāng)?shù)挠嬎憔?。在一些變化波動較為劇烈的地方，本文算法也能與EMTP類方法的結(jié)果基本吻合。本文算法只需要迭代2次即可達(dá)到收斂精度要求，與EMTP類方法的迭代次數(shù)相同，顯然本文算法具有良好的收斂性。

圖4 送端電壓

圖5 受端電壓

在時間積分步長相同的情況下，將EMTP類方法與本文算法的計算時間比值定義為加速比λ1，追趕法與本文算法的計算時間比值定義為加速比λ2。在不同的并行度下，本文算法的計算時間與加速比如表1所示。

表1 算例1 BVPM計算時間及加速比

從表1可以看出，在并行度等于2時，BVPM的計算效率低于EMTP類方法，這是由于BVPM算法是對整個積分區(qū)間一次性求解出所有時間節(jié)點(diǎn)上的值，所得到的矩陣維數(shù)高，在低并行度下，計算機(jī)處理起來花費(fèi)的時間依舊比EMTP類方法更多，但是隨著并行度的增加，BVPM的加速效果逐漸體現(xiàn)出來，計算時間明顯小于EMTP類算法。

3.2 算例2 VFTO 含非線性負(fù)荷的高壓輸電線路

隨著超高壓和特高壓輸電技術(shù)的不斷發(fā)展，氣體絕緣開關(guān)設(shè)備（gas insulated switchgear，GIS）被廣泛運(yùn)用。當(dāng)GIS設(shè)備中隔離開關(guān)投切空載短母線時，開關(guān)觸頭間隙會被多次重復(fù)擊穿，從而在GIS設(shè)備中產(chǎn)生VFTO。因為VFTO具有頻率高、幅值大、陡度高等特點(diǎn)，所以會對電力系統(tǒng)一次設(shè)備和二次設(shè)備帶來一定的危害。因此關(guān)于VFTO的建模、數(shù)值計算等是研究特高壓電力系統(tǒng)電磁暫態(tài)的重要內(nèi)容。圖6是一個550kV GIS中的VFTO計算的簡化模型。其中，L1和L2表示連接的短母線，用無損的短傳輸線模擬，單位長度的電感L0=2.5×10-7H/m、電容C0=4.45×10-11F/m；LT和CT分別表示變壓器的等值電感和對地電容。DS表示隔離開關(guān)，CR表示隔離開關(guān)的對地電容。隔離開關(guān)合閘過程的電弧模型用一個非線性時變電阻來表示為：

式中：τ1=1 ns；τ2=1μs。

圖6 VFTO計算簡化模型示意圖

采用π型等值電路等值短母線，短母線L1和L2分別均分為20段和7段，則可以將VFTO的計算用式（1）微分方程表示。

利用本文算法和EMTP方法對VFTO計算模型進(jìn)行仿真，牛頓迭代的收斂精度設(shè)置為ε=10-4，時間積分步長均取h=0.1 ns。當(dāng)電源側(cè)電壓與短母線殘余電壓之差達(dá)到2.0 p.u.時，VFTO所產(chǎn)生影響最為嚴(yán)重，因此將電源側(cè)電壓設(shè)置為1.0 p.u.，短母線負(fù)荷側(cè)電壓設(shè)置為-1.0 p.u.。

圖7為本文算法計算出的L2線路末端電壓，其中差值曲線是本文算法與EMTP方法比較的結(jié)果。從圖7中可以看出，本文算法與EMTP方法計算結(jié)果差別不大，有著相同的計算精度。

圖7 VFTO計算結(jié)果

表2為本文算法計算時間以及與EMTP方法和追趕法的加速比結(jié)果。

表2 算例2 BVPM計算時間及加速比

從表1和表2可以看出，本文算法的加速效果與理論分析結(jié)果相比稍差。究其原因，一方面是理論分析只考慮了塊矩陣乘法運(yùn)算的時間，忽略了加法運(yùn)算的時間以及生成雅克比矩陣的時間；另一方面是理論分析時并未考慮內(nèi)存管理時間。正是這2個原因?qū)е聦?shí)際的加速效果與理論分析結(jié)果之間的差異。但總的來說，本文算法顯著提高了電磁暫態(tài)仿真的計算效率。

4 結(jié)論

本文將二步二階邊界值法應(yīng)用于電磁暫態(tài)數(shù)值計算，根據(jù)離散后的雅克比矩陣是塊三對角矩陣這一特點(diǎn)，通過一個耦合矩陣，將初始矩陣分解為獨(dú)立的子矩陣，提出了基于二步二階邊界值法的電磁暫態(tài)并行算法。

理論分析了本文算法和追趕法的矩陣計算成本，估計了其加速比。算例結(jié)果表明，本文算法收斂性好，有著與EMTP類方法相當(dāng)?shù)挠嬎憔取Ｓ嬎阈矢哂趥鹘y(tǒng)的追趕法和EMTP類方法，有著良好的加速效果，能夠極大減少電磁暫態(tài)仿真的計算時間。