電力系統(tǒng)靜態(tài)電壓穩(wěn)定臨界點(diǎn)的快速計(jì)算方法

黎 曉， 劉崇茹， 蔡 暉， 辛蜀駿， 鄭 樂

(1.新能源電力系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室(華北電力大學(xué))， 北京 102206；2.國網(wǎng)江蘇省電力有限公司經(jīng)濟(jì)技術(shù)研究院， 江蘇 南京 210008；3.國網(wǎng)經(jīng)濟(jì)技術(shù)研究院有限公司， 北京 102209)

0 引 言

20世紀(jì)70年代以來，國內(nèi)外發(fā)生過多次大的電壓崩潰事故，學(xué)者們對(duì)電壓穩(wěn)定的評(píng)估、預(yù)警與預(yù)防控制開展了大量的研究工作[1～3]。電壓穩(wěn)定臨界點(diǎn)的快速、準(zhǔn)確計(jì)算對(duì)于電壓穩(wěn)定評(píng)估以及電壓穩(wěn)定控制具有重要意義。分岔點(diǎn)對(duì)應(yīng)的負(fù)荷大小可表示當(dāng)前運(yùn)行點(diǎn)的負(fù)荷裕度，可用于電壓穩(wěn)定的預(yù)警與預(yù)防控制。臨界點(diǎn)的靈敏度信息，可以為配置無功補(bǔ)償[4]，切負(fù)荷[5]等措施提供重要信息。

常用的臨界點(diǎn)計(jì)算方法有連續(xù)潮流法和分岔點(diǎn)直接計(jì)算法。連續(xù)潮流法[6]具有魯棒性強(qiáng)、能同時(shí)考慮發(fā)電機(jī)的無功功率限制等約束條件[7]的優(yōu)勢(shì)，很容易處理極限誘導(dǎo)分岔。但是準(zhǔn)確地計(jì)算鞍結(jié)型分岔點(diǎn)，現(xiàn)有的連續(xù)潮流計(jì)算方法計(jì)算量大，難以滿足在線計(jì)算的需求。主要表現(xiàn)在以下兩個(gè)方面：(1)由于PV曲線在鞍結(jié)型分岔點(diǎn)附近的非線性特性明顯，預(yù)測(cè)步長過大會(huì)導(dǎo)致預(yù)測(cè)點(diǎn)距離PV曲線較遠(yuǎn)，使得校正環(huán)節(jié)中牛頓迭代法的迭代次數(shù)顯著增加，甚至導(dǎo)致牛頓法不收斂[8]或擴(kuò)展潮流方程無解[7]。通過步長調(diào)整策略和參數(shù)化過程能夠明顯改善連續(xù)潮流的收斂性，但會(huì)使得連續(xù)潮流的步長在遠(yuǎn)離鞍結(jié)型分岔點(diǎn)時(shí)使用較大的步長，在靠近分岔點(diǎn)時(shí)使用較小的步長[9，10]，這種策略需要較多的計(jì)算量來逼近負(fù)荷極限點(diǎn)。(2)“斜銳角型”PV曲線，可能會(huì)使得校正環(huán)節(jié)中的牛頓法不收斂[11]。這種不收斂現(xiàn)象在全局性參數(shù)化方法與不合適的局部參數(shù)化方法中都可能出現(xiàn)。出現(xiàn)這種情況時(shí)，計(jì)算過程中需要調(diào)整參數(shù)化方法[9]，重新計(jì)算校正步，或者采用線性搜索等求解方式[12]，會(huì)顯著增加總的迭代計(jì)算次數(shù)，降低計(jì)算速度。

另外一類方法基于分岔點(diǎn)的判別條件，直接計(jì)算分岔點(diǎn)[13～17]，被稱為直接法。直接法通常采用牛頓迭代求解，在合適的初值選擇下能夠快速準(zhǔn)確地計(jì)算出鞍結(jié)型分岔點(diǎn)。但是牛頓法對(duì)初值較為敏感，不合適的初值可能導(dǎo)致計(jì)算發(fā)散。文獻(xiàn)[18-20]提出的直接計(jì)算方法，不需要設(shè)置初值，提高了直接迭代方法的收斂性，但是計(jì)算節(jié)點(diǎn)的選擇直接影響到該方法的收斂特性。這類直接法對(duì)于PV曲線呈現(xiàn)“拋物線”形狀的薄弱節(jié)點(diǎn)具有較好的收斂特性，但是當(dāng)初始潮流解遠(yuǎn)離負(fù)荷極限點(diǎn)時(shí)，薄弱節(jié)點(diǎn)的選擇并不容易，因?yàn)镻V曲線的特性在逼近臨界點(diǎn)處才顯現(xiàn)出來[21，22]。文獻(xiàn)[23,24]嘗試將連續(xù)潮流法和直接法的優(yōu)勢(shì)綜合起來進(jìn)行計(jì)算，通過連續(xù)潮流快速穿過鞍結(jié)型分岔點(diǎn)，尋找鞍結(jié)型分岔點(diǎn)的近似解作為直接法的迭代初值，然后用直接法計(jì)算準(zhǔn)確的鞍結(jié)型分岔點(diǎn)。但是，由于前述分析可知，連續(xù)潮流方法在PV曲線上穿過鞍結(jié)型分岔點(diǎn)時(shí)，為了解決收斂性的問題而采用的步長調(diào)整策略或參數(shù)化方法均需要較多的計(jì)算量，限制了該方法的計(jì)算效率。

本文結(jié)合連續(xù)潮流和直接法兩者的優(yōu)勢(shì)，提出了一種靜態(tài)電壓穩(wěn)定臨界點(diǎn)的快速計(jì)算方法。該方法不需要連續(xù)潮流計(jì)算穿過分岔點(diǎn)，避免了連續(xù)潮流為解決收斂問題而增加的計(jì)算量，同時(shí)又能夠?yàn)橹苯臃ㄕ业浇咏?fù)荷極限點(diǎn)的合理初值，從而保證直接法的收斂性。首先，通過大步長的連續(xù)潮流快速靠近臨界點(diǎn)，當(dāng)判斷臨界點(diǎn)為極限誘導(dǎo)分岔時(shí)，則結(jié)束臨界點(diǎn)計(jì)算；當(dāng)判斷臨界點(diǎn)是鞍結(jié)型分岔點(diǎn)時(shí)，以最后一步連續(xù)潮流的解為初值，通過直接法快速準(zhǔn)確地計(jì)算出鞍結(jié)型分岔點(diǎn)。

本文的創(chuàng)新點(diǎn)主要體現(xiàn)在：(1)提出了一種線性程度較好的鞍結(jié)型分岔點(diǎn)逼近指標(biāo)，用以評(píng)估系統(tǒng)逼近鞍結(jié)型分岔點(diǎn)的程度，并根據(jù)這一指標(biāo)確定當(dāng)前連續(xù)潮流計(jì)算結(jié)果能夠滿足直接法初值的要求；(2)提出了連續(xù)潮流與直接法的切換策略，顯著降低了連續(xù)潮流的計(jì)算時(shí)間，同時(shí)保證直接法的初始值接近負(fù)荷臨界點(diǎn)；(3)提出了一種海森矩陣的隱式計(jì)算方法，提高了直接法的計(jì)算效率。

1 鞍結(jié)型分岔點(diǎn)逼近指標(biāo)

1.1 指標(biāo)的理論基礎(chǔ)

鞍結(jié)型分岔點(diǎn)指標(biāo)在第i個(gè)負(fù)荷節(jié)點(diǎn)處的體現(xiàn)記為ci，其在數(shù)值上應(yīng)具備兩個(gè)特征：

(1)在鞍結(jié)點(diǎn)分岔點(diǎn)處ci等于0。保證在鞍結(jié)點(diǎn)分岔點(diǎn)處ci=0能夠保證在計(jì)算過程中通過ci接近于0的程度而快速判斷是否已接近鞍結(jié)型分岔點(diǎn)。

(2)指標(biāo)ci數(shù)值隨著參數(shù)λ增大應(yīng)呈現(xiàn)近似線性的變化趨勢(shì)。由于某一個(gè)參數(shù)λ所對(duì)應(yīng)的指標(biāo)ci并不具有預(yù)測(cè)的能力，此指標(biāo)的關(guān)鍵特征在于其數(shù)值隨著參數(shù)λ增大的線性變化趨勢(shì)。

定義：鞍結(jié)型分岔點(diǎn)逼近指標(biāo)在第i個(gè)負(fù)荷節(jié)點(diǎn)的計(jì)算值ci：

(1)

式中：dVi表示節(jié)點(diǎn)i的電壓幅值的微分；dλ/dVi可以通過連續(xù)潮流計(jì)算過程直接獲??；Ii表示負(fù)荷節(jié)點(diǎn)i的負(fù)荷電流。

指標(biāo)ci還可以寫成式 (2)的形式

(2)

式中：w和τj為方程 (3)的解，w表示PV曲線上切向量的方向；w(i)表示向量w的第i個(gè)元素。fx表示潮流方程的雅各比矩陣，fλ表示潮流方程對(duì)參數(shù)λ的偏導(dǎo)數(shù)，ejT表示除了第j個(gè)元素外其余元素全為0的行向量。τj的上標(biāo)j和ej的下標(biāo)j相對(duì)應(yīng)，與ci的下標(biāo)i無關(guān)，所以可以任取j算出所有負(fù)荷節(jié)點(diǎn)的ci值。n等于除平衡節(jié)點(diǎn)之外的節(jié)點(diǎn)數(shù)量。式 (3)中的τj與Abbott給出的鞍結(jié)型分岔點(diǎn)測(cè)試函數(shù)相同，因此此方法對(duì)應(yīng)相同的理論基礎(chǔ)[25]。已經(jīng)證明[26]，在鞍結(jié)型分岔點(diǎn)處，τj/(w)i=0，因此對(duì)應(yīng)ci=0。

(3)

1.2 指標(biāo)的物理意義

由于視在功率與電壓電流幅值的關(guān)系S=VI，其導(dǎo)數(shù)可以表示為式 (4)。

dS=IdV+VdI

(4)

因此

(5)

式中：dSi/dλ表示節(jié)點(diǎn)i的視在功率的導(dǎo)數(shù)，也是負(fù)荷功率的增長方向在節(jié)點(diǎn)i上的投影，可以視為已知常量，那么容易使得其滿足dSi/dλ0。ZLi=Vi/Ii，表示負(fù)荷等值阻抗模。ZTHi=-dVi/dIi，在其他文獻(xiàn)中被稱為動(dòng)態(tài)等值阻抗模[27]，與通常所采用的戴維南等值阻抗模不同。通常所采用的戴維南等值阻抗采用電壓和電流相量變化量的比值，這里僅僅采用幅值變化量的比值。在物理意義上，式 (5)的形式類似于阻抗模值匹配指標(biāo)，可以同樣解釋為在鞍結(jié)型分岔處，負(fù)荷等值阻抗模等于動(dòng)態(tài)等值阻抗模。

1.3 保持指標(biāo)具有連續(xù)特征的修正方法

與其他鞍結(jié)型分岔點(diǎn)逼近指標(biāo)一樣，在發(fā)電機(jī)節(jié)點(diǎn)類型發(fā)生變化時(shí)，指標(biāo)ci可能會(huì)發(fā)生突變。

因此需要對(duì)式 (1)的定義式進(jìn)行修正，使得該指標(biāo)能夠抵消潮流方程中節(jié)點(diǎn)類型變化帶來的影響。該指標(biāo)在突變點(diǎn)前后具有相似的變化斜率，所以只要消除突變量，即使在節(jié)點(diǎn)電壓類型發(fā)生變化的情況下，保證指標(biāo)隨參數(shù)λ變化的光滑程度仍然較好。

(6)

圖1 鞍結(jié)型分岔點(diǎn)逼近指標(biāo)突變處理示意圖Fig. 1 Schematic diagram of sudden change treatment of saddle-node bifurcation point approximation index

1.4 指標(biāo)的線性化方法

(7)

其中參數(shù)βi可以通過最小方差的方式辨識(shí)，辨識(shí)方法如下。

步驟一：

A=[λ-E(λ),-c+E(c),-λ·c+E(λ·c)]

其中，λ=[λ1,λ2, …,λl]T；c=[ci1,ci2,…,cil]T；λ·c=[λ1·ci1,λ2·ci2,…,λk·cil]T；E(*)表示求平均數(shù)。

步驟二：

式中：Σi表示樣本的方差；Aj,:表示A矩陣的第j行。

步驟三：求解式 (8)中的z向量。

(8)

其中第二個(gè)約束是為了使得βi與ci的符號(hào)異號(hào)，避免式 (7)中分母等于0。z(1)與z(3)分別表示z向量的第1個(gè)元素和第3個(gè)元素。步驟三是一個(gè)帶約束的二次規(guī)劃問題，可以利用現(xiàn)有的工具很容易求解。

步驟四：βi=z(1)

2 連續(xù)潮流與直接法的切換決策

2.1 節(jié)點(diǎn)類型轉(zhuǎn)化識(shí)別的補(bǔ)邏輯

將發(fā)電機(jī)無功功率方程寫成QGi(λ)-fqi(x)=0，其中QGi(λ)表示發(fā)電機(jī)i發(fā)出的無功功率，fqi(x)表示通過節(jié)點(diǎn)電壓計(jì)算發(fā)電機(jī)節(jié)點(diǎn)i的無功功率的表達(dá)式。

(9)

式中：xj表示潮流方程中的變量，包括PQ節(jié)點(diǎn)的電壓幅值和相角、PV節(jié)點(diǎn)的電壓幅值。dxj/dλ由連續(xù)潮流方法中的預(yù)測(cè)環(huán)節(jié)得到。根據(jù)發(fā)電機(jī)無功的上限值和下限值，計(jì)算每個(gè)發(fā)電機(jī)節(jié)點(diǎn)的λ增量的估計(jì)值，如式 (10)所示。

(10)

(11)

式中：GPV表示發(fā)電機(jī)節(jié)點(diǎn)中PV節(jié)點(diǎn)的集合；GPQ表示發(fā)電機(jī)節(jié)點(diǎn)中PQ節(jié)點(diǎn)的集合；dVGi/dλ由連續(xù)潮流方法中的預(yù)測(cè)環(huán)節(jié)得到。然后取λ增量大于0中的最小值作為節(jié)點(diǎn)類型發(fā)生轉(zhuǎn)化的估計(jì)點(diǎn)，如式 (12)和式 (13)所示。

(12)

(13)

但是當(dāng)兩個(gè)節(jié)點(diǎn)幾乎同時(shí)發(fā)生類型轉(zhuǎn)換時(shí)，提前預(yù)判節(jié)點(diǎn)類型轉(zhuǎn)化難以識(shí)別其先后順序，可能收斂到不滿足發(fā)電機(jī)無功限制的解。為了避免這種情況下的多次縮小步長再次計(jì)算校正環(huán)節(jié)，本文增加以下兩個(gè)補(bǔ)邏輯。

(1)當(dāng)預(yù)判節(jié)點(diǎn)i將發(fā)生節(jié)點(diǎn)類型轉(zhuǎn)化，通過計(jì)算收斂之后，節(jié)點(diǎn)j出現(xiàn)了越限，那么將節(jié)點(diǎn)j的類型轉(zhuǎn)換之后，重新計(jì)算節(jié)點(diǎn)i的節(jié)點(diǎn)類型轉(zhuǎn)化臨界點(diǎn)。

(2)如果預(yù)判節(jié)點(diǎn)i和節(jié)點(diǎn)j發(fā)生節(jié)點(diǎn)類型轉(zhuǎn)換時(shí)刻的λ值之差小于，且按照補(bǔ)邏輯1轉(zhuǎn)換節(jié)點(diǎn)i的類型后仍然有新的發(fā)電機(jī)節(jié)點(diǎn)出現(xiàn)越界，那么就恢復(fù)節(jié)點(diǎn)i的節(jié)點(diǎn)類型，再以節(jié)點(diǎn)j作為類型轉(zhuǎn)換的關(guān)鍵節(jié)點(diǎn)計(jì)算其節(jié)點(diǎn)類型轉(zhuǎn)化的臨界值。其中可以取一個(gè)較小的數(shù)，例如=0.001。

2.2 失穩(wěn)類型的識(shí)別

失穩(wěn)類型的正確識(shí)別依賴于節(jié)點(diǎn)類型的識(shí)別和鞍結(jié)型分岔點(diǎn)逼近指標(biāo)的失穩(wěn)判斷。所以在計(jì)算流程上先識(shí)別節(jié)點(diǎn)類型轉(zhuǎn)化，然后判斷失穩(wěn)類型。

每次計(jì)算得到節(jié)點(diǎn)類型轉(zhuǎn)化臨界點(diǎn)之后，通過節(jié)點(diǎn)類型轉(zhuǎn)化前和轉(zhuǎn)化后的切向量的符號(hào)變化可以判斷該點(diǎn)是否為極限誘導(dǎo)分岔點(diǎn)[29，30]。如果是極限誘導(dǎo)分岔，則不需要進(jìn)入鞍結(jié)型分岔計(jì)算的直接法，直接輸出結(jié)果。

當(dāng)鞍結(jié)型分岔點(diǎn)逼近指標(biāo)判斷鞍結(jié)點(diǎn)時(shí)，就不用再進(jìn)入連續(xù)潮流的校正環(huán)節(jié)，轉(zhuǎn)為采用鞍結(jié)型分岔點(diǎn)的直接計(jì)算方法，能夠快速收斂到分岔點(diǎn)。判斷的方法如下：

(1)在每一步連續(xù)潮流法的校正環(huán)節(jié)求解前，綜合突變處理前的指標(biāo)ci與母線電壓幅值選擇薄弱母線k。比較所有負(fù)荷節(jié)點(diǎn)的鞍結(jié)型分岔點(diǎn)逼近指標(biāo)ci值，將其絕對(duì)值從小到大排序，記節(jié)點(diǎn)i的序號(hào)為Ord(ci)。相似地，將負(fù)荷電壓母線的電壓幅值也從小到大排序，記節(jié)點(diǎn)i地序號(hào)為Ord(Vi)，然后計(jì)算Ord(ci)+ Ord(Vi)，取結(jié)果最小的母線編號(hào)為k。

(14)

(15)

整體的具體決策過程見圖2。為了提高算法的魯棒性，保留了回到預(yù)測(cè)環(huán)節(jié)縮短步長再次計(jì)算校正環(huán)節(jié)的處理方式。

圖2 總體計(jì)算流程Fig. 2 Overall calculation process

3 鞍結(jié)型分岔點(diǎn)計(jì)算的改進(jìn)直接法

3.1 直流法的計(jì)算步驟

靜態(tài)電壓穩(wěn)定的鞍結(jié)型分岔點(diǎn)計(jì)算等同于求解以下擴(kuò)展潮流方程中的(x,λ)。

(16)

式中：rank(*)表示矩陣的秩；fx表示潮流方程的雅各比矩陣；g(x,λ)=0是鞍結(jié)型分岔點(diǎn)的判據(jù)，函數(shù)g也被稱為測(cè)試函數(shù)[31]。

方程(16)的解即為鞍結(jié)型分岔點(diǎn)。利用牛頓法求解時(shí)需要形成方程(16)的增廣雅各比矩陣J～，如式 (17)所示。fλ表示潮流方程對(duì)參數(shù)λ的偏導(dǎo)數(shù)。當(dāng)f(x,λ)=f1(x)+f2(λ)時(shí)，fλ與x無關(guān)。gλ=0，gx表示g函數(shù)對(duì)代數(shù)變量x的偏導(dǎo)數(shù)。

(17)

函數(shù)g的表達(dá)式具有不同的等價(jià)形式。本文采用式 (18)形式的函數(shù)g，與式 (16)中的函數(shù)g等價(jià)[31，32]。

g=(d-cTA-1b)

(18)

式中：矩陣A表示fx矩陣中去掉第k行第k列后剩下的部分。向量b表示fx中第k列中除k行之外的元素。cT表示fx中的第k行中除第k列之外的元素。向量d表示fx矩陣中第k行第k列的元素。k是一個(gè)需要設(shè)定的參數(shù)，對(duì)直接法的收斂性有直接影響，本文選擇k值的方法與鞍結(jié)型逼近指標(biāo)預(yù)測(cè)時(shí)選擇薄弱母線k的方法相同。

偏導(dǎo)數(shù)gx需要求解式 (18)中矩陣A的導(dǎo)數(shù)，也是f(x)在x點(diǎn)處展開的海森矩陣的一部分[33]。由于海森矩陣屬于三維矩陣，其維度較高，不利于算法快速求解。

3.2 海森矩陣的快速計(jì)算方法

3.2.1 引理

引理1：一個(gè)矩陣，如果fx，R，T都是連續(xù)函數(shù)或者常數(shù)，那么一定存在唯一的向量u,W, 和實(shí)數(shù)g使得式 (19)和式 (20)成立，且fx矩陣的秩等于2n-r-1的充分必要條件是g=0[34]。

(19)

(20)

本文在引理1的基礎(chǔ)上表示gx，并推導(dǎo)海森矩陣的隱式計(jì)算方法。將式(19)等號(hào)兩邊左乘(uT,g)，與式(20)聯(lián)立得到：

g=uTfxW

(21)

取R=[01×n,1]T，TT=[01×n,1]，uT=(-cTA-1,1)P，W=P-1(-(-A-1b)T,1)T。其中矩陣P為一個(gè)置換矩陣，用于將雅各比矩陣fx中的第k行k列移動(dòng)到最后一行一列。置換矩陣與矩陣fx相乘可以得到如式 (22)所示的分塊形式，其中A∈R(2 n-r-1)×(2 n-r-1),b∈R(2 n-r-1)×1,c∈R(2 n-r-1)×1,d∈R。r表示PV節(jié)點(diǎn)(給定注入有功功率與電壓幅值的節(jié)點(diǎn))的數(shù)量。那么，將R，T，u和W的表達(dá)式代入式(21)中就可以得到式 (18)。說明在以上條件下，式 (18)中的函數(shù)g適用于引理1。

(22)

進(jìn)一步通過式(19)和式(20)求解g對(duì)潮流方程變量x的導(dǎo)數(shù)的表達(dá)式。記g對(duì)任意變量的導(dǎo)數(shù)為g′，對(duì)式(21)等號(hào)兩邊同時(shí)求導(dǎo)得到式(23)。

g′=(uT)′(fxW)+uT(fx)′W+uTfxW′

(23)

根據(jù)式(19)和式(20)中fxW=Rg和uTfx=gTT，式(23)可以改寫為以下形式。

g′=(uT)′(Rg)+uT(fx)′W+gTTW′

(24)

根據(jù)式(19)和式(20)中，TTW=1和uTR=1，得到

TTW′=-(T′)TW

(25)

(uT)′R=-uTR′

(26)

由于向量T向量R為常數(shù)向量，所以式(25)和(26)中的T′和R′都等于全0的向量。因此，將式(25)和式 GOTO(26)代入式(24)，得到g對(duì)任意參數(shù)的導(dǎo)數(shù)，如式(27)所示。

g′=uT(fx)′W

(27)

將式(27)應(yīng)用到對(duì)所有代數(shù)變量求導(dǎo)，得到gx的表達(dá)式。

gx=uTfxxW

(28)

3.2.2 海森矩陣的隱式計(jì)算方法

式(28)中的fxx是潮流方程的海森矩陣，在大規(guī)模電力系統(tǒng)中，海森矩陣維度非常高。盡管對(duì)于電力系統(tǒng)的潮流方程，其海森矩陣高度稀疏，但稀疏矩陣計(jì)算技術(shù)很難在三維矩陣上直接利用。

根據(jù)方向?qū)?shù)的定義，沿著向量W的方向?qū)?shù)等于函數(shù)梯度與方向向量W的內(nèi)積，即

(29)

式中：▽表示梯度算子，有▽f=fx；下標(biāo)(i,j)表示矩陣中第i行第j列的元素；單個(gè)下標(biāo)i或者j表示向量中第i個(gè)或者第j個(gè)元素。

同理可得

(30)

式中：(▽f)x=fxx，下標(biāo)(:,j)表示矩陣的第j列。

那么，

(31)

以上分析過程說明，只需要把向量W當(dāng)作x的一個(gè)增量，在這個(gè)增量方向上兩次形成雅各比矩陣矩陣，就可以很方便的計(jì)算出gx。該計(jì)算過程計(jì)算量小，很容易在現(xiàn)有的潮流計(jì)算程序上實(shí)現(xiàn)。ε可以取一個(gè)較小的數(shù)值，理論上ε取得越小，gx的計(jì)算結(jié)果越精確，但是還需要考慮數(shù)值計(jì)算過程中的舍入誤差，所以不宜取得過小。對(duì)于雙精度浮點(diǎn)數(shù)，可以設(shè)置ε=10-7。

本文提出的海森矩陣隱式計(jì)算方法，顯著提高了直接法的計(jì)算效率。

4 算例分析

在配置為i7-8700 處理器、8 G 內(nèi)存的計(jì)算機(jī)上，用軟件 Matlab 編制本文提出的計(jì)算程序，以算例case2736sp, case2869pagase, case6468rte[35], case9241pegase[36]為例進(jìn)行分析，其網(wǎng)絡(luò)規(guī)模分別為2 736條母線，2 869條母線，6 468條母線，和9 241條母線。本節(jié)選擇了兩種負(fù)荷功率增長方式，增長模式1表示全系統(tǒng)負(fù)荷功率增長的情況，假設(shè)全系統(tǒng)的負(fù)荷有功功率、無功功率等比例增長，全系統(tǒng)的發(fā)電機(jī)(平衡節(jié)點(diǎn)除外)發(fā)出的有功功率也等比例增長以分擔(dān)負(fù)荷功率。增長模式2表示局部區(qū)域負(fù)荷功率增長的情況，假設(shè)編號(hào)最小的10個(gè)負(fù)荷節(jié)點(diǎn)的有功功率與無功功率等比例增長，全系統(tǒng)的發(fā)電機(jī)(平衡節(jié)點(diǎn)除外)等比例增長有功功率以分擔(dān)負(fù)荷功率。大步長連續(xù)潮流采用擬弧長的參數(shù)化方法，直接法中采用隱式海森矩陣計(jì)算方法。

4.1 鞍結(jié)型分岔點(diǎn)逼近指標(biāo)的性質(zhì)

盡管節(jié)點(diǎn)類型轉(zhuǎn)化會(huì)造成鞍結(jié)型分岔點(diǎn)逼近指標(biāo)的突變，圖3表明所提出的式 (6)的處理方式可以得到整體光滑性較好的指標(biāo)變化曲線，鞍結(jié)型分岔點(diǎn)逼近指標(biāo)的突變點(diǎn)的不連續(xù)特性得到了很好的修正，修正后的連續(xù)性是該指標(biāo)在非線性變換后呈現(xiàn)線性特征的基礎(chǔ)。

圖3 鞍結(jié)型分岔點(diǎn)逼近指標(biāo)(非線性變換前，突變修正前后)模式2Fig. 3 Saddle-node bifurcation point approximation index (before nonlinear transformation, whether jump corrected or not)

基于修正后的結(jié)果，通過參數(shù)β的辨識(shí)和式 (7)的非線性變換，得到的鞍結(jié)型分岔點(diǎn)逼近指標(biāo)就具有了比較理想的線性趨勢(shì)，如圖4所示。

圖4 鞍結(jié)型分岔點(diǎn)逼近指標(biāo)(非線性變換后)Fig. 4 Saddle-node bifurcation point approximation index (after nonlinear transformation)

以上結(jié)果驗(yàn)證了所提出的鞍結(jié)型分岔點(diǎn)逼近指標(biāo)的良好特性：在經(jīng)過節(jié)點(diǎn)類型轉(zhuǎn)化造成的突變修正和非線性變換之后，具有較好的線性特征。

4.2 補(bǔ)邏輯的有效性

考慮發(fā)電機(jī)節(jié)點(diǎn)類型轉(zhuǎn)化，擬弧長的預(yù)測(cè)步長σ=1.0，計(jì)算不同算例、不同負(fù)荷增長模式下的靜態(tài)電壓失穩(wěn)分岔點(diǎn)。計(jì)算過程中所觸發(fā)的補(bǔ)邏輯次數(shù)如表 1所示。

表1 補(bǔ)邏輯觸發(fā)次數(shù)比較Tab.1 Comparison of trigger times of complementary logic

表 1中SNB表示鞍結(jié)型分岔，LIB表示極限誘導(dǎo)分岔。結(jié)果表明，在絕大多數(shù)情況下，補(bǔ)邏輯1就足夠有效地避免節(jié)點(diǎn)類型轉(zhuǎn)化造成的縮小步長和再次計(jì)算校正環(huán)節(jié)。僅在case2736sp系統(tǒng)的負(fù)荷增長模式2的計(jì)算過程中觸發(fā)了補(bǔ)邏輯2。

在算例case6468rte中，母線6 372與母線1 366幾乎同時(shí)發(fā)生節(jié)點(diǎn)類型轉(zhuǎn)化。雖然在預(yù)測(cè)步長比較大的情況下(即σ=1.0)通過式 (10)和式 (11)錯(cuò)誤地識(shí)別到母線1 366先發(fā)生節(jié)點(diǎn)類型轉(zhuǎn)化，從而造成母線6 372的無功越界。但是通過補(bǔ)邏輯1，將母線6 372轉(zhuǎn)化為PQ節(jié)點(diǎn)后，計(jì)算結(jié)果正確的收斂到了預(yù)測(cè)步長σ取0.1時(shí)的正確解。

4.3 計(jì)算速度比較

為了比較步長大小對(duì)于鞍結(jié)點(diǎn)附近牛頓法計(jì)算收斂性以及計(jì)算速度的影響，不考慮發(fā)電機(jī)無功功率限制，比較以上4個(gè)算例在不同步長下的計(jì)算時(shí)間，如圖5所示。

圖5 計(jì)算時(shí)間比較Fig. 5 Comparison of computation time

圖5中左側(cè)線表示本文所采用的修正方法：通過鞍結(jié)型分岔點(diǎn)逼近指標(biāo)評(píng)估運(yùn)行點(diǎn)是否接近運(yùn)行極限，然后通過直接法計(jì)算分岔點(diǎn)。右側(cè)線表示采用擬弧長參數(shù)化方法的連續(xù)潮流計(jì)算方法。當(dāng)連續(xù)潮流法計(jì)算過程中校正環(huán)節(jié)的牛頓法計(jì)算發(fā)散時(shí)，將步長縮小為原來的一半，再次計(jì)算校正環(huán)節(jié)[7]。

在連續(xù)潮流法的計(jì)算中，如果逼近臨界點(diǎn)且步長過大，例如圖5 中步長等于20的情況，因?yàn)榕ｎD法計(jì)算不收斂或者收斂慢，造成更多的校正環(huán)節(jié)重復(fù)計(jì)算以及更多的牛頓法迭代次數(shù)，反而降低了計(jì)算速度。

圖5中的計(jì)算結(jié)果表明，在相同預(yù)測(cè)步長的情況下，本文方法比連續(xù)潮流法計(jì)算速度快，這說明通過本文提出的鞍結(jié)型分岔點(diǎn)逼近指標(biāo)預(yù)測(cè)逼近鞍結(jié)型分岔的程度，判斷當(dāng)系統(tǒng)運(yùn)行點(diǎn)接近負(fù)荷增長極限時(shí)，轉(zhuǎn)為直接法計(jì)算鞍結(jié)型分岔具有計(jì)算速度上的優(yōu)勢(shì)。

4.4 計(jì)算精度比較

由于計(jì)算結(jié)果都是收斂后的潮流值，計(jì)算誤差都在10-10以下，因此不論是本文方法還是連續(xù)潮流計(jì)算方法，所得到的計(jì)算結(jié)果都在PV曲線上。

圖6比較了不同步長下的分岔點(diǎn)計(jì)算結(jié)果?？梢钥闯觯?1)本文的方法通過直接法計(jì)算分岔點(diǎn)，計(jì)算精度不受連續(xù)潮流步長的影響。而采用連續(xù)潮流法計(jì)算分岔點(diǎn)的結(jié)果會(huì)受到預(yù)測(cè)步長影響；(2)由于所得到的計(jì)算結(jié)果都在PV曲線上，且本文方法的計(jì)算結(jié)果都大于等于傳統(tǒng)連續(xù)潮流的計(jì)算結(jié)果，說明本文的計(jì)算結(jié)果更加接近鞍結(jié)型分岔點(diǎn)。

圖6 臨界點(diǎn)計(jì)算精度的比較Fig. 6 Comparison of calculation accuracy of critical point

4.5 海森矩陣隱式計(jì)算效率

目前的稀疏矩陣技術(shù)針對(duì)二維矩陣之間或二維矩陣與向量之間的運(yùn)算進(jìn)行了優(yōu)化，且本文的隱式形成海森矩陣能夠比較方便地利用現(xiàn)有的稀疏矩陣技術(shù)。但是顯式形成的海森矩陣屬于三維矩陣，不方便直接利用現(xiàn)有稀疏矩陣技術(shù)，所以兩者之間的計(jì)算時(shí)間難以進(jìn)行有效的比較。因此，本文通過比較計(jì)算的非零元素個(gè)數(shù)與乘法次數(shù)來進(jìn)行說明。

顯式計(jì)算式 (28)中的海森矩陣fxxW，乘法計(jì)算的次數(shù)約等于fxx中非零元素的個(gè)數(shù)，因?yàn)橄蛄縒是稠密的；隱式的海森矩陣形成方法需要形成一次雅各比矩陣的計(jì)算量。根據(jù)fxx矩陣的元素組成，通過估計(jì)非零元素個(gè)數(shù)，從而估計(jì)形成過程中的乘法次數(shù)，如表 2所示。表格中的單位1等于導(dǎo)納矩陣的非0元素個(gè)數(shù)，分析結(jié)果表明，在各個(gè)算例下，隱式形式海森矩陣的計(jì)算量都小于顯式形成海森矩陣的計(jì)算量，而且隱式方法更容易在現(xiàn)有的潮流計(jì)算程序上實(shí)現(xiàn)。

表2 隱式形成海森矩陣與顯式形成海森矩陣的計(jì)算量比較

5 結(jié) 論

本文提出了一種靜態(tài)電壓穩(wěn)定臨界點(diǎn)的快速計(jì)算方法。該方法具有以下3個(gè)特點(diǎn)：

(1)提出了一種線性性質(zhì)較好的鞍結(jié)型分岔點(diǎn)逼近指標(biāo)，利用該指標(biāo)隨負(fù)荷參數(shù)λ變化的線性變化優(yōu)勢(shì)，實(shí)現(xiàn)了對(duì)潮流解逼近鞍結(jié)型分岔點(diǎn)的有效估計(jì)。

(2)提出了連續(xù)潮流法與鞍結(jié)型分岔點(diǎn)計(jì)算直接法的切換策略，通過鞍結(jié)型分岔點(diǎn)逼近指標(biāo)的判斷，在逼近鞍結(jié)型分岔點(diǎn)時(shí)從連續(xù)潮流法切換到直接法，避免了連續(xù)潮流法難以快速準(zhǔn)確地計(jì)算鞍結(jié)型分岔點(diǎn)的問題，同時(shí)為直接法提供了較好的計(jì)算初值，實(shí)現(xiàn)了兩種計(jì)算方法的高效結(jié)合。

(3)提出了一種海森矩陣的隱式計(jì)算方法，具有容易實(shí)現(xiàn)、計(jì)算快的優(yōu)點(diǎn)，提高了此類直接法的計(jì)算效率。在直角坐標(biāo)系下，潮流方程的海森矩陣為常數(shù)矩陣，在本文的基礎(chǔ)上可以進(jìn)一步加速臨界點(diǎn)的計(jì)算。

由于所提出的方法可以方便地處理節(jié)點(diǎn)類型的轉(zhuǎn)化，在本文的基礎(chǔ)上還可以考慮無功補(bǔ)償設(shè)備的其他約束，以及增加新能源機(jī)組和直流輸電系統(tǒng)的靜態(tài)模型，這是本文下一步的工作。