鐘 浩，姚 丹

（三峽大學 梯級水電站運行與控制湖北省重點實驗室，湖北 宜昌 443002）

0 引言

電力系統(tǒng)在當前運行狀態(tài)下的電壓穩(wěn)定程度[1-2]是運行管理人員最為關(guān)心的問題之一。負荷裕度指標[3]在體現(xiàn)電壓穩(wěn)定程度上具有直觀性強和線性度好等優(yōu)點，因此得到了廣泛的應用。負荷裕度指標是指在功率注入空間中，從當前運行點出發(fā)，按給定方向增長負荷直至電壓崩潰時，當前運行點與電壓分岔點之間的距離。因此，要獲得負荷裕度指標關(guān)鍵是要計算得到系統(tǒng)的電壓穩(wěn)定分岔點[4-6]。

若只考慮電力系統(tǒng)的靜態(tài)特性，電壓靜態(tài)分岔類型[7]主要有鞍結(jié)分岔 SNB（Saddle Node Bifurcation）和極限誘導分岔LIB（Limit-Induced Bifurcation）。目前大多數(shù)靜態(tài)電壓分岔點的計算都是針對SNB點。但在某些情況下，系統(tǒng)可能在SNB點之前就發(fā)生了LIB，導致電壓崩潰。只計算系統(tǒng)的SNB點可能會使指標過于樂觀，給系統(tǒng)電壓穩(wěn)定帶來隱患。因此，研究LIB的電壓穩(wěn)定裕度計算方法有重要意義[8-9]。

隨著負荷參數(shù)或其他參數(shù)的逐漸增大，系統(tǒng)的無功源發(fā)生無功越限，可能會伴隨LIB發(fā)生[10-12]?，F(xiàn)階段求取誘導分岔點的方法主要有連續(xù)潮流方法、直接法和優(yōu)化方法。連續(xù)潮流方法[13]可追蹤到LIB點，但是步長難以控制且追蹤時間較長；直接法[14]的應用不如連續(xù)潮流方法普遍，因為直接法求解的非線性方程組維數(shù)高、占用內(nèi)存空間大且不易利用稀疏技術(shù)；優(yōu)化方法[15]易于加入各種約束，但隨著約束條件的增多，模型的收斂性難以保證，計算時間越來越長。

針對LIB點追蹤的計算精度和計算速度問題，本文提出了一種LIB點的快速追蹤方法。首先給出了無功越限導致系統(tǒng)電壓失穩(wěn)的判據(jù)，在遠離分岔點時利用局部曲線擬合技術(shù)，自適應確定負荷增長步長，快速逼近分岔點，并將無功越限的節(jié)點由PV節(jié)點轉(zhuǎn)換成PQ節(jié)點；在靠近分岔點處采用二階靈敏度方法對增長步長內(nèi)PV節(jié)點到PQ節(jié)點轉(zhuǎn)換的PV節(jié)點集的觸發(fā)無功上限排序，校正增長步長，再結(jié)合步長折半搜索方法和分岔點判據(jù)易追蹤到系統(tǒng)的LIB點。

1 負荷裕度求解模型

負荷的有功和無功以及發(fā)電機PV節(jié)點的有功沿某一方向連續(xù)增長，則表征電力系統(tǒng)運行情況的代數(shù)方程可表示為：

其中，x為狀態(tài)變量；ρ為發(fā)電機功率分配因子；λ為負荷參數(shù)。

可將式（1）詳細表示為：

其中，PiG0、QiG0分別為節(jié)點i的初始發(fā)電機有功、無功功率；PiL0、QiL0分別為初始負荷有功、無功功率；Pi、Qi分別為所計算出來的節(jié)點i有功、無功功率；ρiG為發(fā)電機i的分配因子；nL為參與增長的負荷個數(shù)。

負荷參數(shù)λ的變動范圍為：

λ=1對應基態(tài)負荷，λ=λcir對應最大負荷。本文假設負荷按原功率因數(shù)增長，同時也考慮PV節(jié)點的無功極限問題。若某一PV節(jié)點的無功超出其限值，則將其設置成PQ節(jié)點；若在減小步長后其計算無功在極限值之內(nèi)，則將其重新置成PV節(jié)點。在這種PV節(jié)點與PQ節(jié)點轉(zhuǎn)換的過程中，可能會伴隨LIB發(fā)生。

2 LIB點的識別

發(fā)電機無功越限時，發(fā)電機的運行狀態(tài)將發(fā)生變化，節(jié)點類型由PV轉(zhuǎn)換成PQ，此時會有2種情況出現(xiàn)，如圖1和圖2所示。圖1中發(fā)電機無功越限減小了電壓穩(wěn)定裕度但系統(tǒng)未發(fā)生電壓失穩(wěn)，運行點位于節(jié)點轉(zhuǎn)換后λ-U曲線上半支，此時dUj/dλ＜0（j為電壓弱節(jié)點號）；圖2中發(fā)電機無功越限后的LIB導致了電壓崩潰，運行點位于節(jié)點轉(zhuǎn)換后λ-U曲線下半支，此時dUj/dλ＞0。本文重點分析圖2情況。

圖1 極限誘導轉(zhuǎn)換過程Fig.1 Limit-induced transition

圖2 極限誘導分岔Fig.2 Limit-induced bifurcation

在圖2中PV節(jié)點集中有轉(zhuǎn)換節(jié)點a和b，預測點i+1已經(jīng)在可行域之外，使得潮流不收斂，因此需通過式（15）對PV節(jié)點集中的觸發(fā)無功上限排序，將運行點拉回轉(zhuǎn)換點a。在轉(zhuǎn)換節(jié)點a，發(fā)電機a達到無功極限，即QGa=QGliam，Ua=Ua0。隨著負荷參數(shù)增大，其電壓也將增大，此時運行點位于λ-U曲線（QGa=QlGima）下半支，誘發(fā)了系統(tǒng)的突然電壓崩潰，分岔類型為LIB，滿足QGa=QGlima、dUj/dλ＞0 判據(jù)。 若分岔類型為SNB，則可在逼近點斜率絕對值大于某一值情況下基于局部曲線擬合技術(shù)求取SNB點[4]。

將式（1）寫成潮流修正方程的形式，表示為：

式（4）對負荷參數(shù)求導：

在穩(wěn)定運行點處雅可比矩陣J為已知，因此很容易得到dU/dλ。可快速判斷運行點是在λ-U曲線的上半支還是下半支。

3 LIB點追蹤方法

3.1 局部曲線擬合確定增長步長

節(jié)點電壓與負荷參數(shù)呈非線性，為了確保增長步長預測的相對準確性，有必要計算節(jié)點電壓對負荷參數(shù)的2階導數(shù)，通過局部曲線擬合獲得增長步長。繼續(xù)對式（5）求導，得到：

獲得當前運行點的弱節(jié)點電壓對負荷參數(shù)二階靈敏度信息之后，可用局部曲線擬合技術(shù)求取當前運行點到預測點的負荷增長步長。局部曲線擬合式為：

其中，α、β為常系數(shù)；Δλ和ΔUc分別為當前運行點到預測點處的負荷增長步長和電壓弱節(jié)點的電壓幅值變化量。只要計算出α、β和ΔUc即可得到負荷增長步長Δλ。從式（7）可知，越靠近分岔點確定的增長步長越小，具有自適應性。

式（7）對負荷增長步長Δλ求1階導：

繼續(xù)對式（8）求2階導：

由于運行點 i的 dUc/dλ、d2Uc/dλ2已知，聯(lián)立式（8）、（9）可以解出參數(shù) α、β。

然后式（7）對節(jié)點電壓變化量求導，可得：

3.2 PV節(jié)點集觸發(fā)無功上限排序

LIB往往由無功越限所引起，因此在增長步長內(nèi)須計算出發(fā)電機輸出無功，以確定發(fā)電機無功是否越限。如有越限，將其納入PV節(jié)點集，并對節(jié)點集觸發(fā)無功上限排序，識別LIB點。

發(fā)電機PV節(jié)點發(fā)出的無功功率可表示為：

由式（11）可得到在運行點i處發(fā)電機PV節(jié)點的無功出力QPV對負荷參數(shù)λ的各階導數(shù)。其中1階導數(shù)為：

由式（12）再次對負荷參數(shù)λ求導，可得到QPV對負荷系數(shù)λ的2階導數(shù)：

根據(jù)式（14）很容易確定增長步長內(nèi)的PV節(jié)點集：

LIB點最終是通過某臺發(fā)電機無功極限約束來確定的，為了追蹤LIB點，需要對PV節(jié)點集觸發(fā)無功上限順序進行排序，然后取最先觸發(fā)的節(jié)點對應的增長步長Δλ，其求解式如式（15）所示：

其中，ΩPV為當前運行點與預測點之間PV節(jié)點到PQ節(jié)點轉(zhuǎn)換的PV節(jié)點集；nG為PV節(jié)點到PQ節(jié)點轉(zhuǎn)換的節(jié)點個數(shù)；QGlimj為發(fā)電機j無功出力上限；對應中相應值；由于在當前運行點已經(jīng)求出發(fā)電機PV節(jié)點的無功出力QGj，所以根據(jù)發(fā)電機無功出力對負荷參數(shù)導數(shù)很容易確定節(jié)點集中觸發(fā)無功上限的順序。由于Q與λ的非線性，采用二階靈敏度法預測的增長步長可能并不準確，后面將通過負荷增長步長折半搜索方法和設定閾值來重新識別。

系統(tǒng)的無功約束轉(zhuǎn)換點數(shù)目眾多，在遠離分岔點時，無功約束轉(zhuǎn)換點不會使系統(tǒng)發(fā)生電壓崩潰，因此只需將無功越限的PV節(jié)點轉(zhuǎn)換成PQ節(jié)點即可。當離分岔點較近（潮流不收斂情況），才需對預測步長內(nèi)PV節(jié)點集中節(jié)點進行識別。因此，避免了對所有轉(zhuǎn)換點識別所帶來的巨大計算量。

4 算法步驟

下面給出系統(tǒng)電壓穩(wěn)定LIB點追蹤基本步驟。

a.初始化。初始負荷系數(shù)λ0=1.0 p.u.，設置斜率常數(shù) r、斜率比較值 k、最大限值 Δλmax、閾值 ε，運行點斜率 S=0（S=dUj/dλ），標志位 F=0。

b.潮流計算。如果潮流收斂，則令F=0，計算運行點斜率S，跳到步驟c；如果潮流不收斂則跳到步驟e。

c.判斷斜率S正負：如果S為正，則為LIB點，結(jié)束程序；如果S為負且絕對值小于某一設定值，則繼續(xù)向下執(zhí)行；如果S為負且絕對值大于某一設定值，則根據(jù)曲線擬合技術(shù)計算SNB點值，程序結(jié)束。

e.F=0，則取PV節(jié)點集中最先無功越界發(fā)電機的無功/電壓約束轉(zhuǎn)換點作為預測點，并根據(jù)式（15）校正負荷增長步長Δλi，其他轉(zhuǎn)換節(jié)點置回PV節(jié)點，令 F=1，回到步驟 b；F=1 且 Δλi＞ε，則將上步所取的最先轉(zhuǎn)換點還原為PV節(jié)點且負荷增長步長折半，即 Δλi=Δλi/2，回到步驟 b；F=1 且 Δλi＜ε，則該點為LIB點，程序結(jié)束。

靠近分岔點處，PV節(jié)點集中節(jié)點的轉(zhuǎn)換使系統(tǒng)結(jié)構(gòu)發(fā)生突變，可能步長過大會使潮流不收斂，此時用折半搜索方法確保數(shù)據(jù)具有一定連續(xù)性，使本算法具有很好的穩(wěn)定性。另外，在整個過程中，求取λ-U曲線上運行點時具有一定連續(xù)性，對系統(tǒng)計算數(shù)據(jù)隨時存儲和隨時讀取，使得潮流數(shù)據(jù)具有一定的連續(xù)性。因此，算法特點決定了本方案是可行的。

5 算例

采用IEEE 118節(jié)點系統(tǒng)作為仿真算例驗證本文所提方法的快速性和有效性。表1給出了斜率常數(shù)r為1.15和1.20時追蹤LIB點的計算過程，表中n為潮流計算步數(shù)，粗體數(shù)字即為LIB點。

表1 不同參數(shù)逼近過程比較Table 1 Comparison of approaching process between different parameters

當r=1.15、n=8時，預測步長使得潮流不收斂，斜率標記為“—”；在n=9時，通過式（15）取PV節(jié)點集中最先轉(zhuǎn)換的節(jié)點對應的增長步長，潮流不收斂；在n=10時，通過折半搜索方法將轉(zhuǎn)換點置回PV節(jié)點重新識別；當n=11時，將轉(zhuǎn)換點置成PQ節(jié)點，從而S=0.3852為正，說明運行點已位于λ-U曲線下半支，轉(zhuǎn)換點觸發(fā)了LIB，整個追蹤過程如圖3所示。當r=1.20、n=6時，同樣運行點超出可行域，潮流不收斂，此時根據(jù)PV節(jié)點集，通過式（15）確定n=7時的負荷增長步長，將系統(tǒng)運行點拉回可行域；當n=9時，因S=0.0292為正，說明運行點已位于λ-U曲線下半支，轉(zhuǎn)換點觸發(fā)了LIB，整個追蹤過程見圖4。

表2為不同方法的計算值比較，其中NP和NCPF分別為本文所提方法和連續(xù)潮流法的潮流計算次數(shù)，λP和λCPF分別為本文所提方法和連續(xù)潮流法的最大負荷值。本文所提方法計算IEEE118節(jié)點系統(tǒng)，在選取不同斜率常數(shù)r時與連續(xù)潮流法所得最大負荷值[16]的誤差都在2%以內(nèi)，但本文方法所用潮流計算次數(shù)NP分別為11次和8次，而連續(xù)潮流法計算次數(shù)為18次。從算例可見，本文所提方法在計算精度相當?shù)那闆r下效率要高出連續(xù)潮流法1倍左右。

從表2還可以看出，斜率常數(shù)r的取值對潮流計算次數(shù)有一定的影響。常數(shù)r的取值大體思路是：r取值須大于1；當前運行點斜率較小且離分岔點較遠（可根據(jù)算法中增長步長后潮流計算是否收斂來確定當前運行點離分岔點遠近）的時候選擇較大的r值；當前運行點斜率較大且離分岔點較近的時候選擇較小的r值；特別是接近分岔點處，r值越小則計算越精確；在后續(xù)研究中需考慮選擇最佳斜率常數(shù)r的方法。

圖3 r為1.15時逼近LIB點示意圖Fig.3 Schematic diagram of LIB point approaching，when r=1.15

圖4 r為1.20時逼近LIB點示意圖Fig.4 Schematic diagram of LIB point approaching，when r=1.20

表2 不同方法計算值比較Table 2 Comparison of calculative results between different methods

6 結(jié)論

本文通過二階靈敏度確定負荷增長步長和增長步長內(nèi)PV節(jié)點到PQ節(jié)點轉(zhuǎn)換的PV節(jié)點集。由于采用步長限制策略且步長預測具有自適應能力，因此二階靈敏度具有較高的準確性。二階靈敏度只需在潮流計算基礎上增加非常少的計算量即可獲得。此外，本文只需在靠近分岔點處（潮流不收斂情況）才對預測步長內(nèi)PV節(jié)點集中節(jié)點進行識別，避免了對所有轉(zhuǎn)換點識別所帶來的巨大計算量，因此具有非常高的計算效率。本文方法的潮流計算次數(shù)明顯少于連續(xù)潮流方法，且計算精度較高，適用于大規(guī)模電力系統(tǒng)的在線求解。