林沁庭，王 永，郝躍東，吳秀海，張 磊，靳生鵬

（1.三峽大學電氣與新能源學院，湖北宜昌 443002；2.國網(wǎng)上海市電力公司特高壓換流站分公司，上海 201314；3.國網(wǎng)青海省電力公司超高壓公司，青海西寧 810008）

0 引言

電力系統(tǒng)暫態(tài)穩(wěn)定性的數(shù)值計算方法眾多，但基本上可以歸納為顯式分離求解法（Partitioned Explicit Method，PE）和隱式聯(lián)立求解法（Simultaneous Implicit Method，SI）兩大類[1]。其中，隱式聯(lián)立求解法由于不存在交接誤差，而被廣泛采用[1]。

但是，隱式聯(lián)立求解法存在的問題是采用嚴格的Newton 法在聯(lián)立求解過程需要多次形成雅克比矩陣并對其進行三角分解，計算速度較慢。眾所周知，隱式梯形積分法因具有A-穩(wěn)定性且計算精度為二階而被廣泛用于電力系統(tǒng)暫態(tài)穩(wěn)定性分析中。為此，有必要尋找一種新的具有良好數(shù)值穩(wěn)定性的隱式積分法，并結(jié)合與數(shù)值積分算法相適應的高效線性方程組求解方法，則在一定程度上可以滿足暫態(tài)穩(wěn)定性計算的實時性要求。

隱式數(shù)值積分法可大致分為2 類[2]: 第一類是Runge-Kutta 方法，該方法使用較少；第二類是線性多步法，常用于求解微分方程初值問題，包括向后差分方法（Backward Differentiation Formulae，BDF）方法和Adams 方法。除了這2 類方法，還有1 種在隱式線性多步法的基礎上擴展出的被稱為邊界值方法的第三類方法。然而，隱式線性多步法存在著Dahlquist 限制，在現(xiàn)階段眾多學者還無法克服這個問題，因此在求解剛性微分方程問題上受到限制。學者L.Lopez 在基于線性二步法的基礎上，成功地構(gòu)造出了二步邊界值方法[3]（Boundary Value Methods，BVM）。隨后，L.Brugnano 等學者在前人研究的基礎上提出了一系列BVM 方法，包括廣義BDF 方法（Generalized BDF，GBDF）、廣義Adams 方法和擴展的隱式梯形積分法（Extended Trapezoidal Rules，ETRs）。BVM 方法不存在任何階障礙，具有良好的數(shù)值穩(wěn)定性，因此被廣泛應用于求解微分動力系統(tǒng)[4]。

本文將帶參數(shù)二步邊界值方法用于電力系統(tǒng)暫態(tài)穩(wěn)定性計算，并根據(jù)單步塊格式的基本原理得到二步邊界值方法的新的計算格式，減小了牛頓迭代的方程組的維數(shù)，提高了計算的效率。首先，在二步邊界值方法的計算格式的基礎上分析了其數(shù)值穩(wěn)定性。其次，為降低離散后方程組的維數(shù)并提高計算效率，采用了單步塊格式，并推導了基于帶參數(shù)多步塊方法的電力系統(tǒng)暫態(tài)穩(wěn)定性計算公式。最后利用典型算例系統(tǒng)，并與電力系統(tǒng)分析綜合程序（Power System Analysis Software Package，PSASP）軟件作對比測試。數(shù)值結(jié)果表明，文中方法在計算精度、通用性和數(shù)值穩(wěn)定性上效果明顯，取得了較好的結(jié)果，可用于快速判別大擾動后發(fā)電機轉(zhuǎn)子首次搖擺的穩(wěn)定性。

1 邊界值方法

1.1 計算格式

首先需要解決一階常微分方程組的初值問題，其標準形式為[5]：

式中：t為時間變量；y0為狀態(tài)變量在t=0 s 時刻的數(shù)值；y為系統(tǒng)狀態(tài)變量所構(gòu)成的n維向量，y=(y1y2…yn)T∈Rn，Rn為狀態(tài)變量y是n維實數(shù)域內(nèi)的列向量；y′為狀態(tài)變量y的一階導數(shù)；函數(shù)f(t,y(t))為關于t和狀態(tài)變量y(t)的變量；T為待求解問題中t的取值區(qū)間的右端點。

式（1）可改寫成關于t的函數(shù)：

考慮帶參數(shù)β的二步二階邊界值方法[3]離散格式求解式（2），此時式（2）可以離散為：

式中：M為對時間區(qū)間劃分的網(wǎng)格個數(shù)；m為式（3）中代數(shù)方程的個數(shù)；h為時間步長，有h≡tm+1-tm=T/M；fm為函數(shù)f(t)在t=tm=m?h時的函數(shù)值，即fm=f(tm)；ym為狀態(tài)變量y(t) 在t=tm=mh時的數(shù)值。

考慮另外一種帶實數(shù)參數(shù)α的二步三階邊界值方法[3]對式（2）進行數(shù)值積分，有：

聯(lián)立求解式（3）或式（4），需要補充2 個關于狀態(tài)變量的邊界值方程。對于式（2）來說，僅只有首端邊值是提前給定的。為此，需要在式（3）或式（4）的兩端添加二階隱式梯形數(shù)值積分公式[3]：

式中：fM為函數(shù)f(t)在t=tM=Mh時的函數(shù)值，即fM=f(tM)；yM為狀態(tài)變量y(t)在t=tM=M?h時的數(shù)值。

1.2 算法的數(shù)值穩(wěn)定性分析

為驗證二步邊界值方法的數(shù)值穩(wěn)定性，采用的檢驗方程[2]：

式中：λ為實部為負數(shù)的虛數(shù)。

采用單步法以固定步長求解式（6），令變量ξ=λh，構(gòu)造以ξ為變量的帶參數(shù)二步邊界值方法特征多項式，結(jié)合數(shù)值算法的穩(wěn)定性定義[1]，L.Lopez詳細證明了關于帶參數(shù)二步邊界值方法穩(wěn)定性的結(jié)論[3]：

1）當β≤1時，對于具有（1，1）邊界條件的式（3）是二階、01，1-穩(wěn)定及A1，1-穩(wěn)定的數(shù)值積分方法。當參數(shù)β取值[-2，-6]時的穩(wěn)定域見圖1 所示。圖1中穩(wěn)定域定義為封閉曲線以外的區(qū)域。

圖1 二步二階邊界值方法的穩(wěn)定域Fig.1 Stability region of two-step two-order boundary value method

2）當α∈(-∞,-1]?(1,+∞)時，具有（1，1）邊界條件的式（4）是三階、01，1-穩(wěn)定及A1，1-穩(wěn)定的數(shù)值積分方法。當參數(shù)α取值[-9，15]時的穩(wěn)定域見圖2所示。圖2 中穩(wěn)定域定義為封閉曲線以外的區(qū)域。

圖2 二步三階邊界值方法的穩(wěn)定域Fig.2 Stability region of two-step three-order boundary value method

2 帶參數(shù)單步塊方法用于暫態(tài)穩(wěn)定性計算的分析流程

電力系統(tǒng)暫態(tài)過程可用微分-代數(shù)方程組描述為：

式中：x為狀態(tài)變量；u為網(wǎng)絡節(jié)點電壓；Y為電力網(wǎng)絡的節(jié)點導納矩陣；g為系統(tǒng)網(wǎng)絡節(jié)點的注入電流函數(shù)。

為對電力系統(tǒng)暫態(tài)穩(wěn)定性進行實時在線分析，采用經(jīng)典暫態(tài)穩(wěn)定性計算模型：

式中：下標i表示第i臺發(fā)電機；Pei為電磁功率；Pmi為機械功率；Mi為慣性常數(shù)；ωi為轉(zhuǎn)速；δi為轉(zhuǎn)角。

Pei(δi)的表達式為：

式中：Cij，Dij均為常量；i，j均為常量，下標i，j表示待研究電力系統(tǒng)中第i臺發(fā)電機節(jié)點與其它網(wǎng)絡節(jié)點j之間拓撲連接關系。

為了方便分析各發(fā)電機間的相對轉(zhuǎn)子角和角速度，采用慣性中心（Center of Inertia，COI）參考坐標系[12]。首先，定義系統(tǒng)慣性中心下的等值轉(zhuǎn)子角δCOI和等值角速度ωCOI為：

式中：MT為待分析電力系統(tǒng)中總的發(fā)電機慣性時間常數(shù)，。

定義COI 坐標系下的發(fā)電機轉(zhuǎn)子角θi和角速度φi為：

得到COI 坐標系下的發(fā)電機轉(zhuǎn)子運動方程：

式中：PCOI為慣性中心轉(zhuǎn)子運動方程的總功率差額。

Pei(θi)的表達式同式（9）。

重新定義：

式中：z(y,u)為由COI 坐標系下的發(fā)電機轉(zhuǎn)子角和角速度以及網(wǎng)絡節(jié)點電壓表示的狀態(tài)變量y的一階導函數(shù)。

式（12）可簡記為：

式中：為狀態(tài)變量y的一階導數(shù)；Rm×1為狀態(tài)變量y取值為m維實數(shù)域，且m=2n為系統(tǒng)狀態(tài)變量y的維數(shù)。

按照二步二階邊界值方法（參數(shù)β=-2）離散計算式（16），式（5）作為附加方法，如圖3 所示。

圖3 s級單步塊格式的內(nèi)在時間并行性Fig.3 Intrinsic time parallelism of s-level single step block format

s=3 表示將二步邊界值方法轉(zhuǎn)為為單步塊方法后的級數(shù)，h為二步邊界值方法的計算步長，為二步邊界值方法求解過程中內(nèi)點的計算步長，h=s且h=t3-t0。此外，若令tk=t0，tk+1=t3則在時間區(qū)間[tk,tk+1]內(nèi)將s個時刻處狀態(tài)變量按列采用該方法求解式（16），并仿照隱式Runge-Kutta 方法的計算格式表示如下：

當s=3 時，且有相當于時間區(qū)間[tk,tk+1]內(nèi)在s個時間點上同時求解才能得到y(tǒng)k+1，因而具有較好的內(nèi)在時間并行性。根據(jù)式（17）和式（18）可以得到Butcher 表格形式，即，c，和bT均為二步邊界值方法相關的系數(shù)矩陣，其中的表達式如下：

同理，采用二步三階塊邊界值方法（Block Boundary Value Methods，BBVM）（參數(shù)α=-3）離散計算式（16），隱式梯形積分公式（5）作為附加方法，則按照式（17）-式（18）的有如下的系數(shù)矩陣：

通過式（17）—式（20），可以看出此時的2 步塊邊界值方法可以看作是一種單步塊方法，不妨定義新的計算過程的中間變量如式（21）所示：

式中：y，w，Z，d和q均為新的計算過程的中間變量。

則將式（18）和式（7）中的網(wǎng)絡方程聯(lián)立求解可得：

式中：g為系統(tǒng)網(wǎng)絡節(jié)點的注入電流函數(shù)；=-yk；Im為m維的單位矩陣；?為矩陣或向量的直積。

定義牛頓迭代法過程新變量為：

記Δg構(gòu)成新的向量ΔG為：

則利用牛頓法聯(lián)立求解式（22）可得：

按照式（25）的排列方式，式（24）中雅克比矩陣中各分塊矩陣J11，J12，J21和J22的表達式參見文獻[12]。

本文采用廣義極小殘余方法（Generalized Minimal Residual，GMRES）[12]求解方程（25），并采用預處理矩陣提高算法的收斂性。

3 算例分析

算例采用IEEE 145 節(jié)點系統(tǒng)，發(fā)電機模型采用經(jīng)典模型[13]，負荷采用恒阻抗模型[14-20]。故障設定為在7 號母線處發(fā)生三相短路。在故障發(fā)生0.1 s 后切除，整體積分時間為1.50 s。由于104 號母線與7號母線距離最近[21-24]，具體的網(wǎng)絡拓撲圖參考文獻[21]，本文主要分析比較104 號母線上的發(fā)電機功角的計算結(jié)果誤差曲線。數(shù)值計算中以步長取h=0.001 s 時的PSASP 軟件[18]計算結(jié)果（數(shù)值方法為隱式梯形積分法，其中θ表示發(fā)電機之間相對功角，φ表示發(fā)電機之間相對角速度）[8-11]為基準，跟蹤觀察本文方法（s級單步塊邊界值方法BBVM）和同階Runge-Kutta 方法的數(shù)值計算誤差曲線。牛頓法的收斂精度（無量綱）設置為ε=1×10-5。

圖4 和圖5 分別為h等于0.06 s，0.1 s 時的誤差曲線。如圖4 和圖5 所示，本文方法即使采用60倍或10 倍于隱式梯形積分法的計算步長仍可與PSASP 計算結(jié)果保持一致（誤差控制在0.1 之內(nèi)）。此外，與同階的Runge-Kutta 方法相比，保持計算步長不變，本文方法的計算精度更高。

圖4 步長h=0.06 s 時的誤差曲線Fig.4 Error curves when h is 0.06 s

圖5 步長h=0.1 s 時的誤差曲線Fig.5 Error curves for when h is 0.1 s

對于部分比較靠近故障7 號母線的發(fā)電機的相對轉(zhuǎn)子角的搖擺曲線，如第20，24，26，28 和36臺發(fā)電機組的功角θ20，θ24，θ26，θ28和θ36變化曲線如圖6 和圖7 所示?？梢钥闯龉收习l(fā)生后到0.1 s后故障被切除后的1.5 s 內(nèi)，各發(fā)電機的相對轉(zhuǎn)子角均未超過180°，這表明系統(tǒng)在首次搖擺時是暫態(tài)穩(wěn)定的。

圖6 步長h=0.06 s 時的二步二階塊BVM計算的相對功角曲線Fig.6 Relative power angle curves for generators computed by two-step two-order BVM when h is 0.06 s

圖7 步長h=0.06 s 時的二步三階塊BVM計算的相對功角曲線Fig.7 Relative power angle curves for generators computed by two-step three-order BVM when h is 0.06 s

對于部分比較靠近故障7 號母線的發(fā)電機機組的相對角速度曲線，如第20，24 和36 臺發(fā)電機組的轉(zhuǎn)子角速度φ20，φ24和φ36變化曲線如圖8 所示。可以看出當?shù)?0 臺發(fā)電機相對于COI 在減速時，第24 臺發(fā)電機則相對于COI 在加速。這說明在故障擾動下，發(fā)電機轉(zhuǎn)子運動具有分群特性[25-28]。

圖8 步長h=0.06 s 時的二步三階塊BVM計算的相對角速度曲線Fig.8 Relative angular velocity curves for generators computed by two-step three-order BVM when h is 0.06 s

4 結(jié)論

本文將帶參數(shù)的二步塊邊界值方法巧妙地用于電力系統(tǒng)暫態(tài)穩(wěn)定性數(shù)值計算中。帶參數(shù)二步邊界值方法具有二到三階精度，且是A-穩(wěn)定的單步塊方法。針對二步塊邊界值方法差分離散后經(jīng)牛頓法得到的代數(shù)方程組是一個高維的分塊矩陣，采用廣義極小殘余方法求解方程經(jīng)單步塊方法離散后的代數(shù)方程組，并采用預處理矩陣提高算法的收斂性。相對于隱式梯形法及同階的Runge-Kutta 法而言，本文在計算效率及精度上具有明顯優(yōu)勢。通過IEEE145 節(jié)點算例系統(tǒng)仿真表明：本文方法可以快速求解電力系統(tǒng)暫態(tài)穩(wěn)定性，對于機組首次搖擺時的暫態(tài)穩(wěn)定性可以準確判別。此外，文中所提算法可方便實現(xiàn)暫態(tài)穩(wěn)定性的并行計算，相關成果可以進一步推廣至其它類似的計算領域。