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

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

1 引言

隨著現(xiàn)代電力電子技術(shù)的快速發(fā)展，越來(lái)越多的電力電子元件被運(yùn)用到電力系統(tǒng)中，使得電力系統(tǒng)運(yùn)行方式更加靈活、多變。但是隨之而來(lái)的是更加復(fù)雜的控制方式以及更加快速的動(dòng)態(tài)特性。針對(duì)于大規(guī)模多尺度的電力系統(tǒng)，為評(píng)估過(guò)電壓、過(guò)電流對(duì)電力系統(tǒng)穩(wěn)定性造成的影響，維持電網(wǎng)的安全平穩(wěn)運(yùn)行，電磁暫態(tài)過(guò)程的高效仿真顯得極其重要[1]。

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

另一方面，大步長(zhǎng)仿真也是提高電磁暫態(tài)仿真速度的重要手段。動(dòng)態(tài)相量法可以在保持精度的同時(shí)使用大步長(zhǎng)仿真，提高仿真速度[5]。而傳統(tǒng)的動(dòng)態(tài)相量由于存在諧波截?cái)嗾`差[6]，無(wú)法準(zhǔn)確描述含有多種高次諧波系統(tǒng)的動(dòng)態(tài)行為。因此文獻(xiàn)[7，8]采用多頻段動(dòng)態(tài)相量應(yīng)用在電磁暫態(tài)仿真中。

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

而在眾多的多級(jí)隱式方法中，2級(jí)3階Radau II A法具有3階計(jì)算精度，具有A穩(wěn)定性和L穩(wěn)定性[11]，屬于時(shí)域仿真，和動(dòng)態(tài)相量法相比不存在諧波截?cái)嗾`差。但若直接將Radau II A法應(yīng)用于電磁暫態(tài)計(jì)算中，由于引入內(nèi)點(diǎn)，會(huì)造成計(jì)算量倍增的問(wèn)題。因此，將從并行的途徑獲得加速比。

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

2 2級(jí)3階Radau II A方法簡(jiǎn)介

電力系統(tǒng)電磁暫態(tài)過(guò)程通?？捎梢唤M微分方程來(lái)描述

(1)

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

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

(2)

(3)

圖1 時(shí)間網(wǎng)格點(diǎn)圖

利用式(2)對(duì)式(1)在圖1所示的N個(gè)網(wǎng)格點(diǎn)上進(jìn)行連續(xù)離散，可得：

(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維單位矩陣。當(dāng)式(1)描述為非線性初值問(wèn)題時(shí)，利用嚴(yán)格的牛頓法對(duì)式(4)進(jìn)行求解可得

-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公式的電磁暫態(tài)并行計(jì)算

經(jīng)典的Sherman-Morrison公式可描述如下

(23)

式中：A0∈Rn×n；α，β∈Rn×1。對(duì)經(jīng)典的Sherman-Morrison公式進(jìn)行推廣，表達(dá)式如下

(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)

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

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

(29)

(30)

(31)

令

將式(14)變換得

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

(32)

(33)

(34)

(35)

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

(36)

(37)

(38)

(39)

(40)

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

(41)

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

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

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

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

6)利用式(36)求解出ΔX，判斷收斂性，若收斂則進(jìn)行下一個(gè)時(shí)間網(wǎng)格點(diǎn)的計(jì)算，若不收斂；利用ΔX對(duì)X進(jìn)行更新，轉(zhuǎn)至步驟(3)進(jìn)行迭代計(jì)算；

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

f(x)=Utx

(42)

式中，Ut為時(shí)變矩陣。

將式(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 算例仿真

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

4.1 輸電線路連接非線性負(fù)載系統(tǒng)仿真

圖2為輸電線路連接非線性負(fù)載系統(tǒng)圖。

圖2 輸電線路連接非線性負(fù)載系統(tǒng)圖

該系統(tǒng)的基本參數(shù)已在圖2中標(biāo)出。其中，輸電線路L長(zhǎng)為100km，其單位長(zhǎng)度電阻，電感，對(duì)地電容為

R0=0.07Ω/km

L0=2.08×10-3H/km

C0=12×10-9F/km

(47)

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

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

圖3 π型等值電路圖

其數(shù)學(xué)模型表示如下

(48)

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

(49)

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

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

圖4 輸電線路電壓曲線圖

圖5 仿真結(jié)果誤差圖

利用隱式梯形法小步長(zhǎng)仿真，其仿真時(shí)間為7.23秒，利用Radau II A法大步長(zhǎng)串行仿真，其仿真時(shí)間為7.56秒。每個(gè)時(shí)間步長(zhǎng)內(nèi)兩次迭代即完成收斂。

加速比由表1給出。

表1 加速比

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

4.2 VFTO算例仿真

隔離開(kāi)關(guān)投切空載短母線時(shí)會(huì)產(chǎn)生特快速暫態(tài)過(guò)電壓(VFTO)，會(huì)對(duì)與其相連的電氣設(shè)備造成損壞，因此對(duì)VFTO的快速仿真也顯得尤其重要。圖6為產(chǎn)生VFTO的系統(tǒng)圖。

圖6 VFTO系統(tǒng)圖

L1和L2為短母線，L1長(zhǎng)10米，L2長(zhǎng)3.5米，由隔離開(kāi)關(guān)連接，隔離開(kāi)關(guān)對(duì)地電容由CR表示，CR=240pF。而變壓器則由LT和CT組合表示，LT=20mH，CT=3000pF?？紤]最嚴(yán)重的VFTO，即短母線上的殘余電壓與電源電壓差高達(dá)兩倍.p.u(即電源側(cè)取1.0倍.p.u，短母線上取-1.0倍.p.u)。當(dāng)開(kāi)關(guān)合閘時(shí)，其變換過(guò)程由時(shí)變電阻Rt來(lái)表示，其數(shù)學(xué)表達(dá)式為

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

(50)

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

(51)

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

圖7 VFTO仿真結(jié)果和誤差圖

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

表2 加速比

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

5 總結(jié)

所提算法將RadauIIA法和Sherman-Morrison公式結(jié)合起來(lái)，提出了一類(lèi)基于RadauIIA法的電磁暫態(tài)并行計(jì)算方法。并通過(guò)兩個(gè)經(jīng)典算例在精度和加速比兩個(gè)方面的對(duì)比分析中得出了以下結(jié)論：

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

2)RadauIIA法因引入內(nèi)點(diǎn)，具有高階計(jì)算精度，因此可采用較隱式梯形法更大步長(zhǎng)仿真。

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

因此，本文算法可以提高基于2級(jí)3階RadauIIA法的電磁暫態(tài)計(jì)算速度。