電壓穩(wěn)定極限誘導(dǎo)分岔點(diǎn)快速追蹤方法

鐘 浩，姚 丹

（三峽大學(xué) 梯級(jí)水電站運(yùn)行與控制湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 宜昌 443002）

0 引言

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

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

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

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

1 負(fù)荷裕度求解模型

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

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

可將式（1）詳細(xì)表示為：

其中，PiG0、QiG0分別為節(jié)點(diǎn)i的初始發(fā)電機(jī)有功、無(wú)功功率；PiL0、QiL0分別為初始負(fù)荷有功、無(wú)功功率；Pi、Qi分別為所計(jì)算出來(lái)的節(jié)點(diǎn)i有功、無(wú)功功率；ρiG為發(fā)電機(jī)i的分配因子；nL為參與增長(zhǎng)的負(fù)荷個(gè)數(shù)。

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

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

2 LIB點(diǎn)的識(shí)別

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

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

圖2 極限誘導(dǎo)分岔Fig.2 Limit-induced bifurcation

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

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

式（4）對(duì)負(fù)荷參數(shù)求導(dǎo)：

在穩(wěn)定運(yùn)行點(diǎn)處雅可比矩陣J為已知，因此很容易得到dU/dλ?？煽焖倥袛噙\(yùn)行點(diǎn)是在λ-U曲線的上半支還是下半支。

3 LIB點(diǎn)追蹤方法

3.1 局部曲線擬合確定增長(zhǎng)步長(zhǎng)

節(jié)點(diǎn)電壓與負(fù)荷參數(shù)呈非線性，為了確保增長(zhǎng)步長(zhǎng)預(yù)測(cè)的相對(duì)準(zhǔn)確性，有必要計(jì)算節(jié)點(diǎn)電壓對(duì)負(fù)荷參數(shù)的2階導(dǎo)數(shù)，通過(guò)局部曲線擬合獲得增長(zhǎng)步長(zhǎng)。繼續(xù)對(duì)式（5）求導(dǎo)，得到：

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

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

式（7）對(duì)負(fù)荷增長(zhǎng)步長(zhǎng)Δλ求1階導(dǎo)：

繼續(xù)對(duì)式（8）求2階導(dǎo)：

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

然后式（7）對(duì)節(jié)點(diǎn)電壓變化量求導(dǎo)，可得：

3.2 PV節(jié)點(diǎn)集觸發(fā)無(wú)功上限排序

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

發(fā)電機(jī)PV節(jié)點(diǎn)發(fā)出的無(wú)功功率可表示為：

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

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

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

LIB點(diǎn)最終是通過(guò)某臺(tái)發(fā)電機(jī)無(wú)功極限約束來(lái)確定的，為了追蹤LIB點(diǎn)，需要對(duì)PV節(jié)點(diǎn)集觸發(fā)無(wú)功上限順序進(jìn)行排序，然后取最先觸發(fā)的節(jié)點(diǎn)對(duì)應(yīng)的增長(zhǎng)步長(zhǎng)Δλ，其求解式如式（15）所示：

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

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

4 算法步驟

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

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

b.潮流計(jì)算。如果潮流收斂，則令F=0，計(jì)算運(yùn)行點(diǎn)斜率S，跳到步驟c；如果潮流不收斂則跳到步驟e。

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

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

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

5 算例

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

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

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

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

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

圖3 r為1.15時(shí)逼近LIB點(diǎn)示意圖Fig.3 Schematic diagram of LIB point approaching，when r=1.15

圖4 r為1.20時(shí)逼近LIB點(diǎn)示意圖Fig.4 Schematic diagram of LIB point approaching，when r=1.20

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

6 結(jié)論

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