施永豪 董志誠

1(中原工學(xué)院信息商務(wù)學(xué)院 河南 鄭州 450000) 2(西藏大學(xué) 西藏 拉薩 850000)

0 引 言

隨著能源需求的不斷增長，新能源不斷并入電網(wǎng)，使能源管理變得更加復(fù)雜。因此，對(duì)能源管理系統(tǒng)(Energy management systems，EMS)或監(jiān)測(cè)控制與數(shù)據(jù)采集(Supervisory control and data acquisition，SCADA)系統(tǒng)進(jìn)行適當(dāng)?shù)谋O(jiān)測(cè)、優(yōu)化和控制，以維持電網(wǎng)的可靠和持續(xù)運(yùn)行是十分必要的。狀態(tài)估計(jì)是EMS/SCADA系統(tǒng)中輔助這些功能的關(guān)鍵應(yīng)用之一。

自狀態(tài)估計(jì)引入電力系統(tǒng)以來，一直是跟蹤電力系統(tǒng)狀態(tài)的有效途徑。電力系統(tǒng)狀態(tài)估計(jì)主要分為靜態(tài)估計(jì)和動(dòng)態(tài)估計(jì)。靜態(tài)估計(jì)(Static state estimation，SSE)是基于穩(wěn)態(tài)模型的估計(jì)。突發(fā)事件或電能質(zhì)量要求不匹配將導(dǎo)致電網(wǎng)拓?fù)浣Y(jié)構(gòu)快速變化，供電間斷問題時(shí)有發(fā)生，操作人員的應(yīng)對(duì)可能也不足以進(jìn)行解決。因此，監(jiān)控系統(tǒng)應(yīng)配備主動(dòng)控制功能，從而實(shí)現(xiàn)自動(dòng)運(yùn)行[1-2]。為了對(duì)電力系統(tǒng)進(jìn)行連續(xù)監(jiān)測(cè)，必須在較短的時(shí)間間隔內(nèi)進(jìn)行狀態(tài)估計(jì)。此外，隨著具有高測(cè)量頻率和高數(shù)據(jù)交換能力的相量測(cè)量單元(Phasor measurement unit，PMU)的出現(xiàn)，對(duì)資源高效利用的工具需求也進(jìn)一步增大。然而，SSE基于穩(wěn)態(tài)模型，計(jì)算量大，數(shù)據(jù)更新速度慢，通常無法準(zhǔn)確捕捉電力系統(tǒng)動(dòng)態(tài)[1，3]。由此出現(xiàn)了另一個(gè)概念，即動(dòng)態(tài)狀態(tài)估計(jì)(Dynamic state estimation，DSE)。DSE的思想是基于預(yù)測(cè)工具對(duì)估計(jì)值進(jìn)行遞歸更新，其能夠在下一個(gè)時(shí)間步長內(nèi)預(yù)測(cè)系統(tǒng)狀態(tài)。利用DSE可以更準(zhǔn)確地實(shí)現(xiàn)實(shí)時(shí)監(jiān)測(cè)，并且DSE能夠代替丟失的測(cè)量數(shù)據(jù)。因此，這為EMS/SCADA系統(tǒng)通過各種濾波或平滑方法進(jìn)行系統(tǒng)分析和對(duì)電力系統(tǒng)采取控制措施提供了顯著的時(shí)序優(yōu)勢(shì)[1，3]。本文采用最大似然集合濾波器(Maximum likelihood ensemble filter，MLEF)技術(shù)對(duì)一個(gè)基準(zhǔn)電力系統(tǒng)的動(dòng)態(tài)狀態(tài)進(jìn)行估計(jì)，該方法可以幫助電網(wǎng)實(shí)現(xiàn)自治運(yùn)行[2]。

DSE中使用的方法有很多，大多數(shù)都是基于Kalman和Bucy引入的Kalman濾波理論[4]。Kalman濾波器(Kalman filter，KF)被廣泛應(yīng)用于線性高斯系統(tǒng)的狀態(tài)估計(jì)，然而它并不適合非高斯和非線性系統(tǒng)[5-6]。為了將卡爾曼濾波過渡到非線性和非高斯系統(tǒng)，以往研究改進(jìn)得到擴(kuò)展KF(Extended KF，EKF)、集合KF(Ensemble KF，EnKF)、無跡KF(Unscented KF，UKF)和粒子濾波(Particle filter，PF)算法[5-6]，并已有較多文獻(xiàn)將上述算法應(yīng)用于從低維到高維系統(tǒng)的廣泛?jiǎn)栴}中[7-12]。利用DSE的雅可比矩陣對(duì)非線性進(jìn)行線性化，實(shí)現(xiàn)了EKF方法。然而，生成雅可比矩陣以及狀態(tài)估計(jì)和協(xié)方差預(yù)測(cè)所需的計(jì)算量通常很大，且當(dāng)泰勒級(jí)數(shù)展開中的高階項(xiàng)不可忽略時(shí)，線性化近似的誤差較大[13]。EnKF使用蒙特卡羅方法，幫助估計(jì)誤差協(xié)方差，獲得卡爾曼[CD*2]布西濾波器的近似值，并生成可用于預(yù)測(cè)的初始條件集合[11，14-15]。EnKF將非線性嵌入到原始線性KF解中，并使用額外的協(xié)方差來考慮非線性情況[16]。UKF使用近似平均值和協(xié)方差作為若干傳播點(diǎn)(稱為sigma點(diǎn))的線性組合[10，13]，然而，它也包含文獻(xiàn)[18]中所述的缺點(diǎn)。PF使用基于采樣方法的蒙特卡羅模擬來逼近狀態(tài)向量的后驗(yàn)密度，因此它能夠模擬非線性和非高斯性情況。盡管PF是估計(jì)非線性和非高斯系統(tǒng)(包括電力系統(tǒng))的有效工具[12，17]，但它對(duì)采樣的要求較高。MLEF不同于EnKF和PF，它是在狀態(tài)空間而不是樣本空間中起作用的，它通過最大似然方法優(yōu)化非線性函數(shù)，從而減少了計(jì)算時(shí)間。本質(zhì)上，它是SSE中進(jìn)行優(yōu)化的動(dòng)態(tài)版本，解決了隨機(jī)性和不連續(xù)性，同時(shí)它利用了低維空間中的采樣，并使用了文獻(xiàn)[16]和文獻(xiàn)[18]中所提的Hessian信息。在過去的一些研究中，在從低維到高維系統(tǒng)的數(shù)據(jù)同化背景下，MLEF已經(jīng)證明了它在不同領(lǐng)域與其他方法相比的優(yōu)越性[16]。本文首次將MLEF應(yīng)用在電力系統(tǒng)狀態(tài)估計(jì)問題上。

1 最大似然系綜濾波器

最大似然狀態(tài)估計(jì)采用迭代最小化算法進(jìn)行，因此MLEF方法與迭代KF緊密相關(guān)[5]。成本函數(shù)定義為一個(gè)非線性問題，該優(yōu)化方法將MLEF與數(shù)據(jù)同步、控制理論聯(lián)系起來。MLEF采用迭代最小化算法得到最大似然狀態(tài)解。最大似然估計(jì)主要包括預(yù)測(cè)和分析兩個(gè)不同的階段。

1.1 預(yù)測(cè)步驟

預(yù)測(cè)步驟與關(guān)于離散KF的預(yù)測(cè)誤差協(xié)方差變化有關(guān)，表示如下：

(1)

式中：Pf(k)、Mk-1，k和Ω(k-1)分別表示預(yù)測(cè)誤差協(xié)方差矩陣、非線性模型演化矩陣和模型誤差協(xié)方差矩陣；k是時(shí)間參數(shù)；Pa(k-1)表示分析協(xié)方差矩陣。本文忽略了模型誤差。通過因子分解，無時(shí)間指標(biāo)的預(yù)測(cè)誤差協(xié)方差表示如下：

(2)

分析誤差協(xié)方差矩陣的平方根如下：

(3)

Pf(k)1/2=[b1b2…bS]

(4)

bi=M(Xk-1+pi)-M(Xk-1)≈Mpi

(5)

式中：Xk-1是上一步的分析變量；Xk-1+Pi是非線性集合的一步估計(jì)；M(Xk-1)是對(duì)應(yīng)于最可能出現(xiàn)動(dòng)態(tài)的一步控制估計(jì)，并由最大似然方法獲得。

1.2 分析步驟

最大似然法最后可導(dǎo)出最大化后驗(yàn)概率分布值Xk-1。在MLEF的分析步驟中，將問題轉(zhuǎn)化為成本函數(shù)的最小化，其形式如下：

(6)

式中：x表示狀態(tài)向量；y表示觀測(cè)向量；R為觀測(cè)誤差協(xié)方差矩陣；H為非線性觀測(cè)算子；xb為前一步得到的起始狀態(tài)?？紤]到Pf的定義僅限于集合空間，相應(yīng)的Pf的逆也保持在同一子空間范圍內(nèi)。因此，成本函數(shù)也存在于Pf的有效范圍內(nèi)。為了最小化集成子空間中的成本函數(shù)，通過變量的變化來實(shí)現(xiàn)Hessian預(yù)處理。

(7)

式中：ξ表示子空間中控制變量的向量。

(8)

式中：C是一個(gè)信息矩陣，H是雅可比矩陣形式。為了避免可能出現(xiàn)的問題，以獲得矩陣C，它通過已知的平方根預(yù)測(cè)誤差協(xié)方差矩陣來進(jìn)行近似處理如下：

R-1/2H(x+bi)-R-1/2H(x)

(9)

在通過式(7)預(yù)處理之后，非線性優(yōu)化問題式(6)可以通過子空間中基于梯度計(jì)算的迭代最小化來處理。需要注意的是，借助梯度計(jì)算的有限差分近似，避免了矩陣轉(zhuǎn)置的使用[16]。如果所使用的算子是線性的，則通過預(yù)處理最陡下降進(jìn)行的最小化迭代等價(jià)于基于集合的降階KF或基于蒙特卡羅的EnKF。

因?yàn)閦i與觀測(cè)空間具有相同的維數(shù)，如果矩陣Z的定義為Z=[z1z2…zs]，即可將C表示為矩陣Z乘積形式如下：C=ZTZ。通過共軛梯度算法獲得最小解，然后即可更新C，最后計(jì)算平方根分析誤差協(xié)方差如下：

(10)

(I+C)-T/2=E(I+Λ)-1/2ET

(11)

MLEF通過矩陣C的定義與集合變換Kalman濾波器(Ensemble transform Kalman filter，ETKF)具備一定相關(guān)性，而C的特征值分解等價(jià)于ETKF的矩陣變換。

2 基準(zhǔn)模型

本文所用模型如圖1所示。該模型因?yàn)槠滹@示非線性的動(dòng)態(tài)特性，被廣泛應(yīng)用于較多研究。圖2為以文獻(xiàn)[19]中參數(shù)系統(tǒng)進(jìn)行模擬所得的負(fù)載電壓和發(fā)電機(jī)角狀態(tài)相圖。

圖1 帶分接開關(guān)的基準(zhǔn)電力系統(tǒng)模型

圖2 系統(tǒng)相空間軌跡模型

本文模型可視為與較大外部網(wǎng)絡(luò)的局域網(wǎng)連接的等效電路。該模型由一個(gè)負(fù)載節(jié)點(diǎn)和兩個(gè)發(fā)電機(jī)節(jié)點(diǎn)組成，形成一個(gè)三節(jié)點(diǎn)系統(tǒng)[20]和一個(gè)用于電壓控制的變換器[19]。其中一個(gè)發(fā)電機(jī)節(jié)點(diǎn)是表示與外網(wǎng)相連接的節(jié)點(diǎn)。該系統(tǒng)考慮發(fā)電機(jī)的雙軸模型，系統(tǒng)負(fù)載由動(dòng)態(tài)感應(yīng)電動(dòng)機(jī)來表示，并聯(lián)時(shí)則以基于負(fù)載電壓和頻率的恒定PQ節(jié)點(diǎn)表示。負(fù)載節(jié)點(diǎn)還包括一個(gè)電容器，以保持電壓大小在時(shí)刻位于適當(dāng)?shù)姆秶?/p>

模型詳細(xì)數(shù)學(xué)表達(dá)如下：

(12)

(13)

(14)

Kqω(P-P0-P1)-KPω(Q-Q0-Q1)

(15)

式中：M、dm和Pm分別表示發(fā)電機(jī)慣性、阻尼和機(jī)械功率；δm為發(fā)電機(jī)角度(單位為弧度)；ω為發(fā)電機(jī)角速度(單位為rad/s)；δ表示負(fù)載角(單位為rad);V表示負(fù)載電壓(單位為p.u.)；參數(shù)P0和Q0表示電機(jī)的恒定實(shí)際功率和無功功率；P1和Q1表示PQ負(fù)載。K系數(shù)來自感應(yīng)電動(dòng)機(jī)負(fù)荷模型?；鶞?zhǔn)模型的靜態(tài)部分如下:

(16)

(17)

計(jì)算電容的戴維南等效電路的調(diào)整值如下：

(18)

(19)

(20)

將式(16)和式(17)代入式(14)和式(15)獲得δ和V，由此將微分代數(shù)方程形式的系統(tǒng)重新排列成常微分方程(ODE)的形式。此外，將文獻(xiàn)[20]中模型所述的常數(shù)參數(shù)代入方程中，得到的各項(xiàng)系統(tǒng)表達(dá)式如下：

(21)

(22)

20Vcos(0.087 3-δ)+

(23)

0.087 3)]-0.03[-0.6-P1+

0.087 3)]}

(24)

式中：Q1、P1和n為求導(dǎo)后剩下的參數(shù)。因?yàn)橄到y(tǒng)根據(jù)其值表現(xiàn)出不同的混沌行為，因此可稱為分岔參數(shù)。進(jìn)一步，對(duì)上述連續(xù)節(jié)點(diǎn)進(jìn)行離散化后，離散電力系統(tǒng)狀態(tài)空間模型可表示為：

xk=M(xk-1)+Ωk

(25)

式中：狀態(tài)向量xk-1=[δk-1ωk-1δk-1Vk-1]T;Ωk為模型噪聲。

測(cè)量模型如下：

yk=H(xk-1)+Rk

(26)

式中：Rk是測(cè)量噪聲。根據(jù)文獻(xiàn)[12,21]，在此可以用均值為零、方差為σ2的高斯正態(tài)分布進(jìn)行建模。此外，本文對(duì)基準(zhǔn)模型的所有狀態(tài)進(jìn)行了觀測(cè)。

3 算例分析

本節(jié)在基準(zhǔn)模型算例研究中，將使用均方根誤差(RMSE)值進(jìn)行評(píng)估。

3.1 3節(jié)點(diǎn)測(cè)試系統(tǒng)仿真

首先，在圖中用1 000個(gè)集合描述在3節(jié)點(diǎn)電力系統(tǒng)狀態(tài)下一次仿真運(yùn)行的MLEF估計(jì)收斂性。圖3和圖4所示為10個(gè)集合的MLEF和PF之間的估值對(duì)比。表1給出了兩種方法在各個(gè)系統(tǒng)狀態(tài)下模擬的RMSE值。表2給出了兩種方法在各個(gè)系統(tǒng)狀態(tài)下的1 000次蒙特卡洛模擬的平均RMSE值。結(jié)果表明，平均一步估計(jì)時(shí)間為0.039 9 s的10個(gè)集合MLEF與平均一步估計(jì)時(shí)間為2.716 7 s的1 000個(gè)集合PF具有幾乎相同的RMSE結(jié)果。由于仿真的可優(yōu)化性，可以通過采取一定手段減少運(yùn)行時(shí)間。在不同同化時(shí)間下，10個(gè)集合的MLEF估計(jì)RMSE結(jié)果如表3所示。可見，更加頻繁的測(cè)量可以提高RMSE。

圖3 各狀態(tài)下MLEF估計(jì)的仿真結(jié)果

圖4 不同狀態(tài)下MLEF與PF估計(jì)仿真結(jié)果的比較

表1 各狀態(tài)一次模擬運(yùn)行的兩種不同集合數(shù)方法的RMSE估計(jì)

表2 每種狀態(tài)10個(gè)集合下兩種方法100次蒙特卡羅運(yùn)行的平均估計(jì)RMSE

表3 單次仿真運(yùn)行MLEF估計(jì)的RMSE

為了解決干擾問題，系統(tǒng)參數(shù)P1和n分別取為0和1 p.u.的情況下，用10個(gè)同化時(shí)間為0.25的集合進(jìn)行重復(fù)的MLEF估計(jì)，參數(shù)Q1對(duì)系統(tǒng)的干擾在第4秒從11 p.u.增加到11.3 p.u.。圖5展示了針對(duì)該場(chǎng)景的一次模擬運(yùn)行的MLEF估計(jì)性能。

圖5 各狀態(tài)的MLEF估計(jì)的仿真結(jié)果

3.2 68節(jié)點(diǎn)測(cè)試系統(tǒng)仿真

進(jìn)一步將MLEF方法也在IEEE68測(cè)試系統(tǒng)上進(jìn)行了測(cè)試[22]。系統(tǒng)包含68個(gè)節(jié)點(diǎn)和16臺(tái)發(fā)電機(jī)。發(fā)電機(jī)模型采用次暫態(tài)模型和DC1型勵(lì)磁機(jī)。

假設(shè)所有發(fā)電機(jī)的終端總線上安裝了16個(gè)PMU，而對(duì)其他類型的測(cè)量設(shè)備(如遠(yuǎn)程終端裝置)則沒有限制。同時(shí)，PMU的配置位置是隨機(jī)的。PMU測(cè)量的同化時(shí)間選擇為0.04。模擬的集合數(shù)為10。模擬是在節(jié)點(diǎn)23失去負(fù)載的情況下運(yùn)行的，測(cè)量過程考慮了發(fā)電機(jī)的電壓和相角。

穩(wěn)態(tài)下，同步發(fā)電機(jī)的轉(zhuǎn)子角可以從其相量圖中通過端電壓和電流值獲得。暫態(tài)過程中，電機(jī)的電抗會(huì)改變其有效值，通常通過分析法來估計(jì)轉(zhuǎn)子角。對(duì)于實(shí)際測(cè)量值，可以計(jì)算瞬態(tài)電抗值，然后通過求解描述同步電機(jī)動(dòng)力學(xué)方程或使用終端測(cè)量來估計(jì)轉(zhuǎn)子角度。此外，還可以使用合適的估計(jì)方法估計(jì)暫態(tài)電抗及轉(zhuǎn)子角。狀態(tài)估計(jì)將同時(shí)使用模型狀態(tài)和稀疏度量來進(jìn)行。由于問題的性質(zhì)，與完整模型狀態(tài)相比，該度量總是稀疏的，并且該方法具備一定的普適性。

為了在沒有測(cè)量的情況下觀察MLEF的性能，在仿真的第150步后停止，并且在未觀測(cè)的情況下繼續(xù)仿真。對(duì)于發(fā)電機(jī)1和發(fā)電機(jī)4的電壓和相角，圖6展示了系統(tǒng)上一次模擬運(yùn)行的MLEF估計(jì)的收斂性。所有發(fā)電機(jī)的電壓和相角的一次模擬運(yùn)行RMSE結(jié)果如表4所示。

圖6 68節(jié)點(diǎn)系統(tǒng)1號(hào)和4號(hào)發(fā)電機(jī)端電壓和電壓角的MLEF估計(jì)仿真結(jié)果

表4 68節(jié)點(diǎn)系統(tǒng)單次10集合下MLEF估計(jì)的RMSE

4 結(jié) 語

針對(duì)非線性電力系統(tǒng)，本文引入了極大似然函數(shù)。所采用的MLEF是一種將優(yōu)化與低維預(yù)處理相結(jié)合的空間算法，它使MLEF成為低維集成狀態(tài)估計(jì)工具，可以解決沒有定義導(dǎo)數(shù)的非連續(xù)問題。本文通過針對(duì)典型測(cè)試系統(tǒng)的狀態(tài)估計(jì)驗(yàn)證了MLEF的優(yōu)越性，并且在較小的集合規(guī)模下實(shí)現(xiàn)了算法的收斂。在典型模型上的成功試驗(yàn)為今后更復(fù)雜電力系統(tǒng)中狀態(tài)估計(jì)的MLEF性能研究提供了參考。