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

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

0 引言

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

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

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

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

1 邊界值方法

1.1 計(jì)算格式

首先需要解決一階常微分方程組的初值問題，其標(biāo)準(zhǔn)形式為[5]：

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

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

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

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

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

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

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

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

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

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

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

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

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

2）當(dāng)α∈(-∞,-1]?(1,+∞)時(shí)，具有（1，1）邊界條件的式（4）是三階、01，1-穩(wěn)定及A1，1-穩(wěn)定的數(shù)值積分方法。當(dāng)參數(shù)α取值[-9，15]時(shí)的穩(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)定性計(jì)算的分析流程

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

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

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

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

Pei(δi)的表達(dá)式為：

式中：Cij，Dij均為常量；i，j均為常量，下標(biāo)i，j表示待研究電力系統(tǒng)中第i臺(tái)發(fā)電機(jī)節(jié)點(diǎn)與其它網(wǎng)絡(luò)節(jié)點(diǎn)j之間拓?fù)溥B接關(guān)系。

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

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

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

得到COI 坐標(biāo)系下的發(fā)電機(jī)轉(zhuǎn)子運(yùn)動(dòng)方程：

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

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

重新定義：

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

式（12）可簡(jiǎn)記為：

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

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

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

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

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

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

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

式中：y，w，Z，d和q均為新的計(jì)算過(guò)程的中間變量。

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

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

定義牛頓迭代法過(guò)程新變量為：

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

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

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

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

3 算例分析

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

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

圖4 步長(zhǎng)h=0.06 s 時(shí)的誤差曲線Fig.4 Error curves when h is 0.06 s

圖5 步長(zhǎng)h=0.1 s 時(shí)的誤差曲線Fig.5 Error curves for when h is 0.1 s

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

圖6 步長(zhǎng)h=0.06 s 時(shí)的二步二階塊BVM計(jì)算的相對(duì)功角曲線Fig.6 Relative power angle curves for generators computed by two-step two-order BVM when h is 0.06 s

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

對(duì)于部分比較靠近故障7 號(hào)母線的發(fā)電機(jī)機(jī)組的相對(duì)角速度曲線，如第20，24 和36 臺(tái)發(fā)電機(jī)組的轉(zhuǎn)子角速度φ20，φ24和φ36變化曲線如圖8 所示?？梢钥闯霎?dāng)?shù)?0 臺(tái)發(fā)電機(jī)相對(duì)于COI 在減速時(shí)，第24 臺(tái)發(fā)電機(jī)則相對(duì)于COI 在加速。這說(shuō)明在故障擾動(dòng)下，發(fā)電機(jī)轉(zhuǎn)子運(yùn)動(dòng)具有分群特性[25-28]。

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