基于多級高階微分求積法的非線性電磁暫態快速仿真研究

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

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

0 引 言

隨著現代化工業技術的發展，分布式能源的引入，電網結構變得龐大和復雜，而電網不僅可能會受到操作過電壓，雷電過電壓的沖擊，也會面臨不同類型的故障帶來的影響。為了確保電力系統的可靠運行，為保護整定提供參數，也為了滿足電網對實時仿真的要求，對電磁暫態仿真計算速度的研究顯得尤為重要。電磁暫態過程一般由一組微分方程表示，極小的仿真步長使得計算效率低下；精度和速度之間的矛盾制約著電磁暫態仿真效率的提高。為了在保證計算精度的同時加快電磁暫態計算速度，研究者在不同方向展開探討，并取得了可觀的結果[1-4]。

并行技術作為一種提高電磁暫態計算速度的有效手段活躍在電磁暫態仿真中。圖形處理器(GPU)結合所選算法的特點將電磁暫態計算中硬件設備的優點充分發揮出來[5-6],數據的并行處理極大提高了電磁暫態計算速度。而算法上的并行實現更是為電磁暫態計算效率的提高做出了卓越的貢獻[7-9],得到了可觀的加速比。

另一方面，大步長仿真也是提高電磁暫態仿真的重要手段。文獻[10]成功實現了基于時間尺度變化下的大步長電磁暫態仿真；而動態相量的引入也帶來更多電磁暫態大步長仿真的可能性[11-13]。動態相量建模法可以實現大步長仿真的本質是因為作了降頻處理，使快變化變成慢變化，但動態相量也暴露出很多不足，動態相量對非線性元件建模依然復雜，并且動態相量將時域信號轉變為復數信號，造成系數矩陣增維，也增加了計算量。總的說來，對動態相量的研究還有待增加。

電磁暫態計算中數值方法的選擇也為大步長仿真提供了選擇性。經典的隱式梯形法由于具有A穩定性，且擁有二階計算精度，廣泛應用在電磁暫態計算中。但隱式梯形法不具有L穩定性，即無法有效抑制突變情況下的振蕩特性，且其仿真步長小，因此更多多級多階的方法被應用于電磁暫態計算中[14-15]。其中微分求積法原理簡單，易于實現。具有A穩定性且可以有效抑制震蕩。文獻[16]已經證明了微分求積法與龍格庫塔(RK)方法的等值性，其計算結果為s級s階精度。和隱式梯形法相比可以采取更大的仿真步長。但對于s級微分求積法來說隨之而來的是內點的引入，由于需要在每個步長里引入s個內點，造成系數矩陣維數增加到原來的s倍，因此在矩陣求逆上會花費大量時間，使得大步長計算失去意義。文中算法在微分求積法的基礎上，結合V變換的特點，根據增維后的矩陣特點進行分塊求解，提高矩陣求逆的速度，以此來提高基于微分求積法的電磁暫態計算速度。

1 微分求積法簡介

若函數f(x)在區間上光滑可導，則f(x)在網格點ci,i∈(1,s)上的導數可由該區間上所有網格點處的函數值的線性加權和表示，即：

(1)

式中gij為加權系數;s為整個區間進行劃分的網格點數，不包括起始端點。將式(1)寫為矩陣形式為:

f(1)(c)=G0f(c0)+Gf(c)

(2)

其中：

(3)

(4)

且：G0=-Ge，e為s維單位列向量。微分求積法加權系數矩陣G一般有兩種求解方式，一種是將拉格朗日插值基函數作為試函數帶入式(1)求解，叫顯示表達式；另一種是將一般多項式函數帶入式(1)求解，其結果由范德蒙德矩陣及其逆矩陣的乘積來表示，稱為隱式表達式，其結果如下：

(5)

其中：

(6)

(7)

αs=[α1,α2,…,αs]T=(1/s)V-1cs

(8)

(9)

對于網格點的選取，常見的有勒讓德網格、切比雪夫網格、均勻網格等，采取的是切比雪夫網格，其表達式為：

(10)

電磁暫態過程是由磁場和電場的變化引起電壓、電流的變換過程，由于變換速度快，一般由一組帶有初值的微分方程來表示：

(11)

式中x∈Rm×1，x為系統中狀態變量的集合;g(t)只與時間變量有關;x0則為t0(t=0)時刻狀態變量x的初值。

將積分區間[tn,tn+1]正則化到c∈(0,1)上，并將式(11)帶入式(2)，可得：

(12)

(13)

令：A=G-1≡[aij]，?示直積計算，則可由下式來表示微分求積法求解步驟：

(14)

2 基于微分求積法的電磁暫態計算

將式(14)寫為矩陣格式則可得：

(15)

式中：

(16)

式中Im為m維單位矩陣，利用A=G-1，將式(15)變換得：

(17)

當式(12)描述為非線性問題時，定義：

(18)

用牛頓迭代法對式(17)進行整體求解，ξ為迭代次數。帶入式(18)定義的初值，則：

-JΔx(ξ)=?R(x(ξ))

(19)

式中：

(20)

(21)

(22)

式中Is為s維單位陣。令：

(23)

帶入式(5)，則式(19)變為：

(24)

式中：

(25)

(26)

(27)

式(19)可以寫為：

(28)

如果直接求解式(28)，在每一次迭代中都會對增維矩陣(s×m維)進行三角分解，若m很大將會增加很大的計算量，使得微分求積法大步長的優勢失去意義。觀察式(28)，由于As矩陣具有特殊結構，因此其逆矩陣也具有特殊結構，表達式為：

(29)

將式(29)帶入式(28)，并且展開分塊可得：

(30)

(31)

令：

(32)

(33)

(34)

(35)

(36)

將式(32)～式(36)帶入式(30)和式(31)可得：

(37)

觀察可得，H2為下三角矩陣，因此只需要進行簡單的前代運算就可以將式(37)第一式中的y2表示成y1的線性表達式，而不需要進行復雜的三角分解，即：

y2=φ(y1)

(38)

再將式(38)帶入式(37)的第二式，得：

(39)

Q∈Rm×m，利用LU三角分解，解出y1，帶入式(38)，解出y2，在通過式(23)反解出Δx即可。通過判斷是否小于收斂精度，小于收斂精度則計算下一個步長，大于收斂精度，則需要按照內點平均值作為初值方式進行下一次迭代，直到收斂到精度范圍內則開始進行下一個步長的計算。

當初值問題為線性時變微分方程時，即：

f(x)=Utx

(40)

Ut只與時間變量有關。將式(40)帶入式(17)，則有：

(41)

式中：

(42)

G′(t)=(G?Im)(e?xn)+hG(t)

(43)

令：

(44)

做近似處理：

(45)

帶入式(42)，則式(42)中對角矩陣塊變為相等的矩陣塊。式(41)變為：

(46)

(47)

將式(46)按照非線性部分算法處理，進行分塊處理。即可獲得可觀的加速比。

3 算例仿真

上述算法是在微分求積法基礎上提出來的，旨在解決由于內點引入，使得矩陣增維降低計算速度的問題。此部分中引入三個算例，分別為非線性和線性時變系統算例，并從精度和速度兩個方面進行討論，證明文中算法的有效性。算例仿真部分的算法都是在三級微分求積法上形成的。

3.1 加速比定義

定義文中算法仿真時間與其他算法仿真時間之比為加速比，加速比為仿真時間之比，沒有量綱。具體規定為：加速比K1表示傳統串行微分求積法與文中算法的仿真時間之比(步長一致)；加速比K2表示小步長隱式梯形法與文中算法的仿真時間之比；加速比K3表示小步長隱式梯形法與傳統串行微分求積法的仿真時間之比。大步長與小步長在每個算例里均有具體的定義。

3.2 含非線性負荷的高壓輸電線路仿真

含非線性負荷的高壓輸電線路系統如圖1所示。在系統合閘時，輸電線路首端和末端將會產生過電壓，需要評估過電壓對輸電線路電力設備絕緣的影響程度，為了提高過電壓電磁暫態過程的數值計算效率，需要在更短的時間內得到其仿真結果。

圖1 含非線性負荷的高壓輸電線路仿真圖Fig.1 Simulation diagram of high voltage transmission line with nonlinear load

圖1中，輸電線路L首端與電源，電阻Rs，電感Ls串聯，輸電線路末端接有一非線性負載，由定電阻RL和非線性電感LL組合表示，LL中磁鏈φ與iL的數學關系為：φ=atanhbiL。輸電線路全長100 km。其余參數如表1所示。

表1 參數表Tab.1 Parameter table

本算例中采用分段π模型電路對輸電線路進行建模。每段電路模型如圖2所示。

圖2 π型等值電路Fig.2 π model equivalent circuit

本算例里將輸電線路分為30段，每段輸電線路的電磁暫態過程由式(48)表示：

(48)

仿真總時長為0.004 s，t=0時刻合閘。輸電線路首端末端約束方程為：

(49)

將輸電線路方程和首尾兩端約束方程聯立起來，最終含非線性負載的輸電線路電磁暫態數學模型為：

(50)

式中x為系統狀態變量；U為系統線性部分系數矩陣；f(x,t)為系統非線性部分。

分別用隱式梯形法、串行的微分求積法對式(50)進行離散，牛頓-拉夫遜迭代求解離散后的非線性方程組,再用文中算法對式(50)進行離散求解。隱式梯形法仿真步長1 μs(后文稱為小步長)，串行的微分求積法以及文中算法仿真步長10 μs(后文稱為大步長)。設置收斂精度為10-4，三種算法每個步長內只需要迭代兩次就能完成收斂。

圖3為采用文中算法仿真得出的輸電線路首端和末端的電壓，圖4為文中算法(大步長)和隱式梯形法(小步長)仿真出的輸電線路首端電壓和末端電壓波形的誤差曲線。需要注意的是由于串行微分求積法和文中算法在離散方法和步長的選擇上是一致的，兩種算法的精度誤差在10-5之內，故不再給出誤差曲線。

由圖3和圖4可知，文中算法步長是隱式梯形法步長的10倍，但兩者仿真精度誤差如圖4所示，說明微分求積法能在大步長仿真下依舊保持仿真精度。加速比以及計算時間以表2給出。

圖3 輸電線路首端末端電壓Fig.3 Voltage at the front end of the transmission line

圖4 文中算法和隱式梯形法精度誤差圖Fig.4 Algorithm proposed in this paper and the implicit trapezoidal method precision error graph

由于傳統串行方式下，每一次迭代求解過程都會進行s×m維矩陣的LU分解，即使采用大步長計算，增維矩陣的三角分解也使得計算效率低下，反映為表2中的K3值，其加速比不明顯。對于文中算法來說，在采用大步長的計算基礎上，每一次迭代過程中都會將s×m維矩陣的三角分解轉換為(s-1)×m維矩陣的前代運算以及一個m維矩陣的三角分解，且相對于傳統串行微分求積法來說，這種加速優勢將會隨著一次又一次的迭代累積，最終形成可觀的加速比，反映為表2中的K1、K2值。

表2 加速比及計算時間Tab.2 Acceleration ratio and calculation time

3.3 柔性高壓直流系統算例

圖5為柔性高壓直流系統圖(VSC-HVDC)，下標為1代表系統整流側，下標為2代表系統逆變側。對上述系統有兩點假設：

圖5 柔性高壓直流系統圖Fig.5 Flexible HVDC system diagram

(1)換流橋均由理想閥元件組成，即其正向漏電流為零；

(2)系統電壓、電流滿足三相平衡條件且為工頻正弦波。

整流側、逆變側時域方程為：

(51)

直流側時域方程為：

(52)

式中：

將式(51)和式(52)整合成矩陣形式方程：

(53)

式中:

VSC-HVDC系統參數：兩端電壓源電壓分別為3 kV和2.8 kV,系統頻率為50 Hz，兩側電氣參數一致，即：

t=0.06 s時，VSC其中一側調制比由0.78下降到0.65，t=0.12 s時恢復到初始值。總的仿真時長為0.3 s。

為了滿足電網對實時仿真的要求，提高電磁暫態仿真速度是極其必要的。下面分別用隱式梯形法對式(53)進行離散仿真，步長為0.0001 s(小步長)；微分求積法對式(53)進行離散仿真，串行求解，仿真步長為0.001 s(大步長)；文中算法對式(53)離散仿真，步長為0.001 s(大步長)。圖6給出了由文中算法仿真出來的直流電壓波形圖以及與小步長隱式梯形法仿真結果的精度誤差曲線。由于微分求積法串行算法在離散算法和步長選擇上與文中一致，兩者之間的仿真誤差在10-8之內，故不再給出其誤差曲線。

圖6 文中方法計算結果Fig.6 Calculation results of the proposed method

由圖6可以看出，即使微分求積法仿真步長為隱式梯形法仿真步長的10倍，卻依然可以達到很高的計算精度。加速比由表3示出。

表3 加速比及計算時間Tab.3 Acceleration ratio and calculation time

由表3可以看出：雖然相比隱式梯形法，微分求積法增大了仿真步長，但由于每個步長引入s個內點使得計算量增加，因此微分求積法串行方式下其加速比K3只是略大于1，即加速優勢不明顯，而文中算法在采用大步長的前提下，針對引入內點增加計算量這點提出改進方法，在每一個計算步長里將s×m維矩陣的求逆計算分解為(s-1)×m維矩陣的前代運算和m維矩陣的三角分解，減少計算量，得到了較為可觀的加速比。

3.4 VFTO算例仿真

氣體絕緣開關設備中的隔離開關投切空載短母線時會產生特快速暫態過電壓(VFTO)，具有頻率，幅值，陡度高的特點，對與其相連設備造成安全威脅。為了給保護整定提供電氣參考值，維護設備正常運行，對VFTO進行高效率電磁暫態仿真是必要的。圖7是一個簡化的VFTO計算模型。

圖7 VFTO系統圖Fig.7 VFTO system diagram

圖7中，電源電壓幅值為550 kV，LT，CT組合表示變壓器等效電感和對地電容，L1，L2為連接短母線，由無損短傳輸線表示，L1長度為10 m，L2長度為3.5 m。CR表示隔離開關對地電容。隔離開關合閘過程的電弧模型由以下數學表達式來模擬：

R(t)=R1e-t/τ1+R2e-t/τ2

(54)

無損傳輸線單位長度的電感和對地電容分別由L0，C0表示。所有參數值由表4提供。

表4 參數表Tab.4 Parameter table

隔離開關投切空載短母線時，短母線上的殘余電荷電壓對VFTO有重要影響。考慮最嚴重的VFTO，即電源側與殘余電荷電壓差達到2.0 p.u.(電源側為1.0 p.u.,負荷側殘余電荷電壓取-1.0 p.u.)。

輸電線路采用分段π型等值電路，L1均分為20段，L2均分為7段。仿真時長為1.2 μs。VFTO數學模型表達式為：

(55)

Ut∈R57×57,且為時變系數矩陣，分別用隱式梯形法(步長為0.1 ns，小步長)，文中算法(步長為0.6 ns,大步長)，傳統的串行微分求積法(步長為0.6 ns)進行離散仿真。仿真結果如圖8所示，u(t)表示用隱式梯形法仿真出來的母線末端電壓，v(t)表示由文中算法仿真出來的末端電壓，Δuv(t)表示兩者之間的精度誤差。需要說明的是，由于文中算法和串行微分求積法所采用的離散方法和步長都是一致的，兩者之間的計算精度誤差在10-8之內，故不再給出誤差曲線。

圖8 VFTO結果對比圖Fig.8 VFTO result comparison diagram

由圖8可以看出，微分求積法采用大步長計算，其仿真結果依然保持一定精度。

加速比及計算時間在表5中給出。

表5 加速比及計算時間Tab.5 Acceleration ratio and calculation time

Ut為時變系數矩陣，對于每一個步長來說，Ut依然為定系數矩陣。因此每一個步長的計算都需要求解s×m維矩陣的逆。文中算法依舊將s×m維矩陣的求逆計算分解為(s-1)×m維矩陣的前代運算和m維矩陣的三角分解，因此會得到可觀的加速比。

3.5 結果分析

通過上述三個算例的分析，可以發現，文中算法對非線性和線性時變系統是有明顯加速效果的，其加速比是通過每次迭代或者是每個步長的快速求逆累積得到的。對于非線性電力系統算例來講，s的取值不會影響非線性部分算法的迭代次數。所以無論是文中算法，串行微分求積法，還是隱式梯形法，每個步長的迭代次數均為兩次。3.3節以及3.4節給出了兩個線性時變算例，雖然都獲得了較好的加速效果，但從這兩個算例可以看出，矩陣規模越大其加速效果是更加明顯的。而對于線性常微分算例來講，由于系數矩陣為定常數矩陣，整個計算仿真過程只需要進行一次求逆計算。只需要不斷更新初值和時間函數變量。加速優勢無法通過累積得到，因此文中算法對此類算例效果不明顯。

需要說明的是：針對非線性電力系統，無論其拓撲結構復雜程度如何，其電磁暫態過程最終可寫為式(50)的形式；而對于線性時變系統來說，其電磁暫態過程可描述為式(53)的形式，利用所提算法對式(50)或式(53)進行離散仿真，便可以在保證精度的前提下得到可觀的加速比。因此所提算法針對不同拓撲結構的非線性系統以及線性時變系統均具有有效性。

4 結束語

文中算法是在微分求積法基礎上形成的，利用微分求積法V變換的特點，對增維矩陣進行分塊處理，避免直接對增維矩陣進行三角分解，使得需要三角分解的矩陣由串行模式下的s×m維降低到m維，因此大大提高了基于微分求積法的電磁暫態計算速度。且算例表明，所提算法對非線性算例以及線性時變算例均有理想的加速效果，且矩陣規模越大其加速效果更加明顯。文中算法可以基于任何級的微分求積法。因此，文中算法可以有效提高基于微分求積法的電磁暫態仿真效率。