基于最大似然集合濾波器的電力系統狀態估計方法

施永豪 董志誠

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

0 引 言

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

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

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

1 最大似然系綜濾波器

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

1.1 預測步驟

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

(1)

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

(2)

分析誤差協方差矩陣的平方根如下：

(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)是對應于最可能出現動態的一步控制估計，并由最大似然方法獲得。

1.2 分析步驟

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

(6)

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

(7)

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

(8)

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

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

(9)

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

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

(10)

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

(11)

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

2 基準模型

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

圖1 帶分接開關的基準電力系統模型

圖2 系統相空間軌跡模型

本文模型可視為與較大外部網絡的局域網連接的等效電路。該模型由一個負載節點和兩個發電機節點組成，形成一個三節點系統[20]和一個用于電壓控制的變換器[19]。其中一個發電機節點是表示與外網相連接的節點。該系統考慮發電機的雙軸模型，系統負載由動態感應電動機來表示，并聯時則以基于負載電壓和頻率的恒定PQ節點表示。負載節點還包括一個電容器，以保持電壓大小在時刻位于適當的范圍。

模型詳細數學表達如下：

(12)

(13)

(14)

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

(15)

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

(16)

(17)

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

(18)

(19)

(20)

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

(21)

(22)

20Vcos(0.087 3-δ)+

(23)

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

0.087 3)]}

(24)

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

xk=M(xk-1)+Ωk

(25)

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

測量模型如下：

yk=H(xk-1)+Rk

(26)

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

3 算例分析

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

3.1 3節點測試系統仿真

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

圖3 各狀態下MLEF估計的仿真結果

圖4 不同狀態下MLEF與PF估計仿真結果的比較

表1 各狀態一次模擬運行的兩種不同集合數方法的RMSE估計

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

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

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

圖5 各狀態的MLEF估計的仿真結果

3.2 68節點測試系統仿真

進一步將MLEF方法也在IEEE68測試系統上進行了測試[22]。系統包含68個節點和16臺發電機。發電機模型采用次暫態模型和DC1型勵磁機。

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

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

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

圖6 68節點系統1號和4號發電機端電壓和電壓角的MLEF估計仿真結果

表4 68節點系統單次10集合下MLEF估計的RMSE

4 結 語

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