基于最大似然集合濾波器的電力系統(tǒng)狀態(tài)估計方法

施永豪 董志誠

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

0 引 言

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

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

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

1 最大似然系綜濾波器

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

1.1 預測步驟

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

(1)

式中：Pf(k)、Mk-1，k和Ω(k-1)分別表示預測誤差協(xié)方差矩陣、非線性模型演化矩陣和模型誤差協(xié)方差矩陣；k是時間參數(shù)；Pa(k-1)表示分析協(xié)方差矩陣。本文忽略了模型誤差。通過因子分解，無時間指標的預測誤差協(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是非線性集合的一步估計；M(Xk-1)是對應于最可能出現(xiàn)動態(tài)的一步控制估計，并由最大似然方法獲得。

1.2 分析步驟

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

(6)

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

(7)

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

(8)

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

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

(9)

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

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

(10)

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

(11)

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

2 基準模型

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

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

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

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

模型詳細數(shù)學表達如下：

(12)

(13)

(14)

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

(15)

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

(16)

(17)

計算電容的戴維南等效電路的調整值如下：

(18)

(19)

(20)

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

(21)

(22)

20Vcos(0.087 3-δ)+

(23)

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

0.087 3)]}

(24)

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

xk=M(xk-1)+Ωk

(25)

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

測量模型如下：

yk=H(xk-1)+Rk

(26)

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

3 算例分析

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

3.1 3節(jié)點測試系統(tǒng)仿真

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

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

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

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

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

表3 單次仿真運行MLEF估計的RMSE

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

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

3.2 68節(jié)點測試系統(tǒng)仿真

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

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

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

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

圖6 68節(jié)點系統(tǒng)1號和4號發(fā)電機端電壓和電壓角的MLEF估計仿真結果

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

4 結 語

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