基于Radau II A法的電磁暫態并行計算

張 靜，殷 明，葉 婧，李博文

(1. 三峽大學電氣與新能源學院，湖北 宜昌 443002；2. 國網湖北省電力有限公司隨州供電公司，湖北 隨州 441300)

1 引言

隨著現代電力電子技術的快速發展，越來越多的電力電子元件被運用到電力系統中，使得電力系統運行方式更加靈活、多變。但是隨之而來的是更加復雜的控制方式以及更加快速的動態特性。針對于大規模多尺度的電力系統，為評估過電壓、過電流對電力系統穩定性造成的影響，維持電網的安全平穩運行，電磁暫態過程的高效仿真顯得極其重要[1]。

并行計算是電磁暫態高效快速仿真的有效手段之一。圖形處理器(GPU)具有多核心、高集成度的特點[2]，結合所選算法特征，GPU被應用于電磁暫態并行計算中[3，4]，提高計算速度。

另一方面，大步長仿真也是提高電磁暫態仿真速度的重要手段。動態相量法可以在保持精度的同時使用大步長仿真，提高仿真速度[5]。而傳統的動態相量由于存在諧波截斷誤差[6]，無法準確描述含有多種高次諧波系統的動態行為。因此文獻[7，8]采用多頻段動態相量應用在電磁暫態仿真中。

算法的多樣性也為大步長仿真提供了可能性。隱式梯形法作為一種經典的離散方法被廣泛的應用于電磁暫態計算中，具有二階計算精度。但其仿真步長小，仿真速度慢，近年來，更多多級隱式的方法憑借著大步長仿真的優勢被應用在電力系統計算中。但是多級就意味著雅克比矩陣增維，從而導致計算量增加，使得采用大步長仿真失去意義，因此多級隱式法并不適合直接應用于電磁暫態計算中。文獻[9]將多級的廣義向后差分方法應用于電磁暫態穩定性計算中，通過對大矩陣進行分塊處理，減少計算量，得到可觀的加速比。文獻[10]也將一種多級的方法應用于電磁暫態計算中，通過并行處理三對角塊矩陣，從而提高了計算速度。

而在眾多的多級隱式方法中，2級3階Radau II A法具有3階計算精度，具有A穩定性和L穩定性[11]，屬于時域仿真，和動態相量法相比不存在諧波截斷誤差。但若直接將Radau II A法應用于電磁暫態計算中，由于引入內點，會造成計算量倍增的問題。因此，將從并行的途徑獲得加速比。

概括起來，所提算法的創新點在于，利用Sherman-Morrison[12，13]公式求解代數方程，相比文獻[13]，其得到的高維雅克比矩陣的非對角塊為一個高度稀疏的常系數矩陣，因此避免了牛頓迭代中因更新雅克比矩陣而需要重新對其分解的過程。并且在求逆之前，先對大矩陣進行分塊處理，相比文獻[12]，其非對角塊矩陣元素減少，進一步降低串行計算量。

2 2級3階Radau II A方法簡介

電力系統電磁暫態過程通常可由一組微分方程來描述

(1)

式中，x∈Rm×1，為電力系統狀態變量的集合，g(t)僅與時間變量有關，t0為積分起始時刻，x0為起始時刻電力系統狀態變量的初值。

2級3階Radau II A法屬于多級隱式Runge-Kunttu(RK)系列方法。其積分格式為

(2)

(3)

圖1 時間網格點圖

利用式(2)對式(1)在圖1所示的N個網格點上進行連續離散，可得：

(H?Im)X=(B?Im)F(X)+(B?Im)G(t)+Const(t0)

(4)

式中

(5)

(6)

(7)

(8)

(9)

F(X)=[fT(x1)fT(x2)…fT(xN)]Ti∈(1，N)

(10)

G(t)=[gT(t1)gT(t2)…gT(tN)]Ti∈(1，N)

(11)

g(ti)=[gT(ti-1+c1h)gT(ti-1+c2h)]Ti∈(1，N)

(12)

(13)

Im為m維單位矩陣。當式(1)描述為非線性初值問題時，利用嚴格的牛頓法對式(4)進行求解可得

-JΔX=ΔF

(14)

(15)

(16)

(17)

(18)

(19)

ΔF=(H?Im)X-(B?Im)F(X)-

(B?Im)G(t)-Const(t0)

(20)

令

(21)

i∈(1，N)

(22)

其中：x0，2=x0

3 基于Sherman-Morrison公式的電磁暫態并行計算

經典的Sherman-Morrison公式可描述如下

(23)

式中：A0∈Rn×n；α，β∈Rn×1。對經典的Sherman-Morrison公式進行推廣，表達式如下

(24)

式中：αk，βk∈Rn×1，k∈(1，n)，且

P=[pk]，pk∈Rn×1，k∈(1，n)

(25)

Q=[qk]，qk∈Rn×1，k∈(1，n)

(26)

R=diag(rk)，k∈(1，n)

(27)

(28)

對于式(14)的求解，將從并行的途徑獲得加速比。由式(15)可見，J由對角矩陣和一個高度稀疏的常系數矩陣構成，因此可以巧妙地結合Sherman-Morrison公式求逆。

再觀察式(28)，每一個pk，k∈(1，n)都由前k-1個pi，i∈(1，k-1)線性遞推得到，當并行度越高，即n值越大，對式(28)的求解仍然具有一定的計算量，為了在保證并行度的同時，減少計算項，先對式(15)進行分塊處理。

(29)

(30)

(31)

令

將式(14)變換得

-(Jd+ΔJ)ΔX=ΔF

(32)

(33)

(34)

(35)

利用擴展的Sherman-Morrison公式求解式(32)，可得

(36)

(37)

(38)

(39)

(40)

對式(36)的求解，可以得到以下結論：

(41)

由式(41)可知，對于JJi的求逆計算轉換為了對Ji，i∈(1，N)的求逆計算，并且彼此解耦，沒有增加矩陣求逆的計算量。

所提算法的計算步驟大致可以概括如下：

1)利用Radau II A法在N個時間點上連續離散，設置收斂精度為ε，并帶入初值，嚴格牛頓法求解非線性方程組，得到式(14)；

2)對式(14)中形成的大雅可比矩陣進行分塊處理，得到式(29)，式(30)以及式(31)；

6)利用式(36)求解出ΔX，判斷收斂性，若收斂則進行下一個時間網格點的計算，若不收斂；利用ΔX對X進行更新，轉至步驟(3)進行迭代計算；

若式(1)描述為線性時變微分方程時，即

f(x)=Utx

(42)

式中，Ut為時變矩陣。

將式(42)帶入式(4)，并整理合并得：

JX=C

(43)

式中：

(44)

(45)

其中：Utij=U(ti-1+cjh)，i∈(1，N)，j∈(1，2)

C=(B?Im)G(t)+Const(t0)

(46)

按照非線性部分算法，將J分塊后分解為Jd+ΔJ，利用Sherman-Morrison公式求解，后文不再贅述。

4 算例仿真

由于本部分兩個算例均無法求出解析解，因此其仿真結果將會和小步長隱式梯形法仿真結果做對比。并且定義小步長隱式梯形法與本文算法的仿真時間之比為加速比k1，串行的大步長Rudau II A法與本文算法的仿真時間之比為k2，小步長和大步長在具體算例里都給出了定義。由于本文算法和串行Radau II A采用的步長和離散方法均是一致的，其仿真結果精度一致，因此只給出本文算法仿真結果圖。兩個算例的仿真均有16個計算單元參與并行計算。

4.1 輸電線路連接非線性負載系統仿真

圖2為輸電線路連接非線性負載系統圖。

圖2 輸電線路連接非線性負載系統圖

該系統的基本參數已在圖2中標出。其中，輸電線路L長為100km，其單位長度電阻，電感，對地電容為

R0=0.07Ω/km

L0=2.08×10-3H/km

C0=12×10-9F/km

(47)

和輸電線路相連的負載由定電阻RL和非線性電感LL組合表示。LL中的非線性關系表示為：φ=atanhbiL。a=8.40×102V·s，b=5.95×10-3A-1。

輸電線路由π型等值電路建模。將輸電線路L分成M段，每段由圖3所示線路表示。

圖3 π型等值電路圖

其數學模型表示如下

(48)

補充輸電線路L的首端尾端約束方程如下

(49)

聯立式(48)和式(49)，總的仿真時長為0.004秒。將輸電線路分為30段，收斂精度設置為10-4。整理可得式(1)表達式。其中x∈R63×1。利用隱式梯形法對本算例進行離散仿真，其仿真步長為10-6秒，后文稱為小步長；利用Radau II A法串行仿真；利用本文算法進行仿真，由于Radau II A法為多級方法，因此可以取大步長，將大步長取為隱式梯形法小步長的2s倍，因此本文將步長取為4×10-6秒，后文稱為大步長。

圖4給出了輸電線路首端末端電壓曲線，圖五給出了與隱式梯形算法的誤差圖。由圖5可見，即使采用了大步長仿真，其仿真結果依然保持了一定的精度。

圖4 輸電線路電壓曲線圖

圖5 仿真結果誤差圖

利用隱式梯形法小步長仿真，其仿真時間為7.23秒，利用Radau II A法大步長串行仿真，其仿真時間為7.56秒。每個時間步長內兩次迭代即完成收斂。

加速比由表1給出。

表1 加速比

由表1可以看出，采用RadauIIA法大步長仿真的時間和隱式梯形法小步長仿真的時間大致相等，這是由于引入內點增加計算量的原因。因此若直接將RadauIIA法運用到電磁暫態計算中，將不會產生明顯的加速效果。所提算法在N個時間點上實現了并行求逆，并且基于非對角塊矩陣高度稀疏和常系數的特點，使得串行過程只涉及簡單的向量計算，并且分塊處理后相比于分塊前減少一半計算量，進一步地加快了所提算法的計算速度。且當并行度越高所獲得的加速比也越高。

4.2 VFTO算例仿真

隔離開關投切空載短母線時會產生特快速暫態過電壓(VFTO)，會對與其相連的電氣設備造成損壞，因此對VFTO的快速仿真也顯得尤其重要。圖6為產生VFTO的系統圖。

圖6 VFTO系統圖

L1和L2為短母線，L1長10米，L2長3.5米，由隔離開關連接，隔離開關對地電容由CR表示，CR=240pF。而變壓器則由LT和CT組合表示，LT=20mH，CT=3000pF。考慮最嚴重的VFTO，即短母線上的殘余電壓與電源電壓差高達兩倍.p.u(即電源側取1.0倍.p.u，短母線上取-1.0倍.p.u)。當開關合閘時，其變換過程由時變電阻Rt來表示，其數學表達式為

Rt=R1e-t/τ1+R2e-t/τ2

(50)

式中：R1=1012Ω，R2=0.5Ω，τ1=1ns，τ2=1μs。本算例里將短母線L1分為20段，將短母線L2分為7段，單位長度電感L0=2.5×10-7H/m，單位長度電容C0=4.45×10-11F/m。采用算例1用的π型等值電路建模的方法對輸電線路L1，L2建模。得到其數學模型為

(51)

式中：Ut∈R57×57，為時變矩陣。仿真總時長為1.2微秒，利用隱式梯形法對式(51)進行離散仿真，步長取為0.1ns(小步長)；RadauIIA法串行仿真，和算例一同樣的原理，將步長取為2s倍隱式梯形法仿真步長，即0.4ns(大步長)；本文算法進行大步長仿真。圖7給出了輸電線路L2末端電壓波形以及與隱式梯形法仿真結果的誤差曲線圖，由圖可見，RadauIIA法具有大步長高精度仿真的優勢。

圖7 VFTO仿真結果和誤差圖

加速比如下表可見，k1，k2定義和算例1一致。隱式梯形法小步長仿真時間為5.871秒，RadauIIA法大步長串行仿真時間為6.023秒。

表2 加速比

和算例1同樣的分析方法，所提算法利用Sherman-Morrison公式并行求解對角塊矩陣，分塊處理減少了非對角塊元素，進一步降低計算量，并獲得了可觀加速比。

5 總結

所提算法將RadauIIA法和Sherman-Morrison公式結合起來，提出了一類基于RadauIIA法的電磁暫態并行計算方法。并通過兩個經典算例在精度和加速比兩個方面的對比分析中得出了以下結論：

1)所提算法均是在嚴格牛頓迭代法上形成的，因此具有嚴格牛頓法的收斂性。

2)RadauIIA法因引入內點，具有高階計算精度，因此可采用較隱式梯形法更大步長仿真。

3)通過對高維雅克比矩陣進行分塊處理，利用Sherman-Morrison公式并行求逆，獲得了較好的加速比，且并行度越高，加速比越大。

因此，本文算法可以提高基于2級3階RadauIIA法的電磁暫態計算速度。