考慮量測相關性的容積卡爾曼濾波動態狀態估計

盧慶春，張俊，許沛東，陳思遠，徐箭，柯德平

(武漢大學 電氣與自動化學院， 武漢 430072)

0 引 言

電力系統狀態估計用于獲得電力系統各節點母線的運行狀態，是能量管理系統(Energy Management System，EMS)的重要組成部分[1]。電力系統狀態估計有兩種，包括靜態狀態估計及動態狀態估計[2]。靜態狀態估計是對電力系統某一時刻的狀態進行估計，而動態狀態估計則利用當前時刻的狀態估計結果對下一時刻狀態進行預測，并根據下一時刻實際量測對狀態預測值進行修正。因此，相比于靜態狀態估計，動態狀態估計具有預測環節，其可以實時動態地估計系統狀態，可以更好地對電力系統狀態進行監測。文中研究電力系統動態狀態估計。

動態狀態估計主要是基于卡爾曼濾波方法進行的[3]，國內外對其已經有大量的研究。最開始有文獻采用擴展卡爾曼濾波(Extended Kalman Filter，EKF)[4-5]的估計方法，而EKF對電力系統中非線性方程進行展開的方法存在著線性化誤差[6]。為了解決EKF方法存在的問題，無跡卡爾曼濾波(Unscented Kalman Filter，UKF)[6-9]被提出，UKF通過對狀態量進行離散化采樣，然后將采樣點代入非線性方程中進行無跡變換，以此計算離散化采樣點的離散化函數結果，因而可以避免EKF線性化展開帶來的誤差。但所有基于UKF的狀態估計方法均需要進行無跡變換參數的選取，而參數選取的好壞關系著狀態估計結果的準確度。因此，為了克服UKF存在參數選取誤差的缺點，有學者提出了容積卡爾曼濾波方法(Cubature Kalman Filter，CKF)[10-13]。文獻[10]建立了四階發電機動態模型，并采用CKF進行發電機的狀態估計；文獻[11]中提出了可以修正狀態預測值的魯棒CKF方法，并用于估計發電機的狀態；文獻[12]提出了一種基于Huber M估計的魯棒CKF方法來應對量測噪聲非高斯分布的情況；文獻[13]將一種魯棒的CKF估計方法用于估計同步發電機的狀態，其對同步發電機狀態模型進行改進，并可以處理輸出和輸入量測中存在的不良數據，仿真驗證算法能明顯提高同步發電機的估計精度。上述文獻中的誤差協方差矩陣均為對角矩陣，通常假設各個量測的誤差標準差相同，難免對估計結果造成影響。

然而，量測量間是存在相關性的，且量測誤差協方差矩陣是非對角矩陣。文獻[14]通過分析變電站內電力系統量測信號的采集過程驗證了量測間存在相關性，并用點估計法計算量測相關性矩陣，提出了考慮量測相關性的最小二乘估計算法，仿真驗證了考慮量測相關性的最小二乘估計算法相比傳統最小二乘估計算法能提高狀態估計的精度及計算效率。之后，一些文獻開始研究考慮量測相關性的狀態估計。文獻[15]中考慮了變電站量測相關性，采用兩點估計法計算了變電站內所有量測間協方差矩陣，并通過仿真驗證考慮量測相關性能提高狀態估計的精度；文獻[16]分析了最小二乘算法中考慮量測相關性的重要性，通過仿真驗證考慮量測相關性能提高狀態估計的精度。

目前量測相關性主要應用在靜態狀態估計中，在動態狀態估計中尚缺乏考慮量測相關性的研究。因此，文中提出了一種考慮量測相關性的容積卡爾曼濾波動態狀態估計方法。首先對SCADA系統量測量存在的相關性進行分析，然后基于雙指數平滑預測過程計算過程噪聲協方差矩陣，采用容積變換計算考慮SCADA系統量測相關性的量測誤差協方差矩陣，并提出考慮量測相關性的電力系統動態狀態估計流程，最后在IEEE-39節點進行仿真驗證文中方法動態估計效果。

1 電力系統動態狀態估計

1.1 動態狀態估計模型

電力系統動態狀態估計以卡爾曼濾波為基礎， 假設非線性系統的狀態轉移及量測模型為：

xk=f(xk-1)+qk-1

(1)

zk=h(xk)+rk

(2)

式中xk為k時刻系統狀態向量;f(·)為狀態轉移函數；zk為k時刻系統量測向量；h(·)為系統量測函數；qk-1和rk分別為k-1時刻和k時刻的過程噪聲向量和量測誤差向量，方差分別為Qk-1和Rk。

(3)

1.2 CKF基本原理

由于CKF算法沒有EKF的線性化誤差及UKF中選擇無跡變換參數帶來的參數選取誤差，文中基于CKF進行狀態估計。CKF狀態估計基礎是容積變換，容積變換指的是以一個狀態量的值作為采樣的均值點，依據三階球面—徑向法則[17]產生一些具有相同權值的容積點，然后將每個容積點按照非線性方程進行變換，從而完成對非線性方程的近似。CKF包括預測步以及濾波步，如下所示:

(1)預測

假設k-1時刻系統狀態量估計誤差方差矩陣為Pk-1。對Pk-1進行Cholesky分解，求得k-1時刻估計誤差方差陣的平方根矩陣：

(4)

(5)

(6)

根據式(5)產生的容積點，狀態量以及估計誤差方差矩陣的預測值可以按下式計算：

xi,k|k-1=f(χi,k-1)

(7)

(8)

(9)

(2)濾波

zi,k=h(γi,k)

(10)

因此，量測預測值為：

(11)

預測量測自協方差矩陣為：

(12)

Tk=Sk+Rk

(13)

狀態預測及量測預測的互協方差矩陣為：

(14)

卡爾曼濾波增益矩陣可以按下式計算：

(15)

因此，k時刻的狀態估計值以及估計誤差方差矩陣可以表示為如下形式：

(16)

(17)

2 考慮量測相關性CKF動態狀態估計

結合目前電網中廣泛配置的SCADA量測，文中研究SCADA量測間的相關性。電力系統動態狀態估計的量測方程h(x)為：

(18)

式中Pi為節點注入有功功率；Qi為節點注入無功功率；Pij為支路有功功率；Qij為支路無功功率；Vi和Vj為節點電壓幅值；Gij為節點導納矩陣i行j列元素電導部分；Bij為節點導納矩陣i行j列元素電納部分；θij=θi-θj為支路兩端節點相角差；gij為節點導納矩陣i行j列元素互電導部分；bij為節點導納矩陣i行j列元素互電納部分；y為支路對地導納。

2.1 SCADA系統量測相關性分析

對于SCADA系統，狀態估計所用的量測量包括節點電壓幅值、節點注入功率(有功/無功)、支路潮流功率(有功/無功)。除去節點電壓幅值，其余量測由電力系統中間裝置對互感器量測進行中間計算得到。文中假設每個節點或線路單獨安裝了互感器。

假設互感器電壓幅值量測量為：

(19)

假設中間量測電壓相角差為：

(20)

因此，式(18)可以寫成：

(21)

由式(21)可知，Pi、Qi、Pij、Qij均受互感器量測誤差εi和εij影響，因此這些量測量間存在著相關性。Pi和Pj由于受εi和εij影響，兩個量測量間也存在相關性，其他量測間的相關性同理可知。

因此，SCADA系統量測量間存在相關性。誤差協方差矩陣Rk是非對角陣，對角線元素表示量測量的誤差自方差，非對角線元素表示的是各量測量之間的誤差互方差。Rk可以表示為：

(22)

(23)

(24)

2.2 過程噪聲協方差矩陣Qk

過程噪聲方差矩陣Qk關系到狀態估計結果的精度，因此Qk的合適選取十分關鍵。第2.2節采用的過程如下：

(1-α)3xk-3|k-4+bk-1+(1-α)bk-2+(1-α)2bk-3

(25)

對bk-1進行展開得到：

(26)

同理，也可對bk-2、bk-3及bk-n進行展開。

由狀態轉移方程可知，k時刻的系統狀態預測值與k時刻之前的所有時刻均有關，且受鄰近時刻的影響較大。因此，過程噪聲是動態變化的，隨系統狀態的變化而時變。

由于α和β取值均為[0，1]，當系數次數越高，取值越小。當相鄰時刻系統狀態變化不大時，其差值很小，為加快求解效率，取狀態轉移方程為：

(27)

根據誤差傳遞規則，可以求得過程誤差協方差矩陣為：

(28)

2.3 考慮量測相關性的量測誤差協方差矩陣Rk

誤差協方差矩陣Rk中非對角元素體現著量測間的相關性，因此考慮量測相關性的狀態估計不能假設Rk為一個特定矩陣。由于容積變換方法能夠根據自變量的協方差矩陣容易得到因變量的協方差矩陣，因此，文中基于容積變換法計算Rk。

假設互感器量測的電壓幅值和相角差量測向量為：

(29)

向量維度為m1。求解互感器量測量的協方差矩陣為：

Hk=E[(Lk-E(Lk))(Lk-E(Lk))T]

(30)

依照第1節中的容積變換思想，求取2m1個容積點ζi,k。

(31)

(32)

其中，{1}i為矩陣{1}的第i個列向量{1}如下所示：

(33)

將容積點代入量測方程，得到SCADA系統的量測量如下：

Zi,k=h(ξi,k)

(34)

量測均值為：

(35)

取量測量誤差協方差矩陣近似為：

(36)

2.4 考慮量測相關性CKF動態狀態估計流程

文中考慮量測相關性的容積卡爾曼濾波方法流程如下：

(4)根據k時刻互感器量測得到的互感器量測向量Lk，通過容積變換方法求得容積點，通過式(34)～式(36)計算量測誤差協方差矩陣Rk；

(6)計算狀態預測及量測預測的互協方差矩陣Ck；

(7)計算卡爾曼濾波增益Kk，k時刻狀態估計值，以及k時刻狀態協方差矩陣Pk；

(8)更新X_est，循環步驟(1)～ 步驟(8)，進行動態狀態估計。

3 算例分析

3.1 算例設計

為驗證文中所提考慮量測相關性的動態狀態估計算法的估計效果，第3.1節在IEEE-39節點系統(圖1)上進行了仿真分析，并和不考慮量測相關性的CKF方法進行對比。通過在DIgSILENT上進行時域仿真產生一系列采樣數據，在Matlab上進行動態狀態估計。時域仿真時，對系統設置短路故障，在t=0.18 s時，節點16發生接地短路故障，并在t=0.36 s清除故障，時域仿真時間為10 s。

對于不考慮量測相關性的CKF方法(記為傳統CKF)，其量測誤差協方差矩陣的對角線元素取為文中量測誤差協方差矩陣Rk對角線元素，狀態轉移方程采用文中狀態轉移方程。

仿真分為兩種情況：(1)各個互感器電壓幅值量測誤差及中間量測電壓相角差的誤差均服從均值為0、方差為10-4的正態分布，即εi～N(0,10-4)，εij～N(0,10-4)，i和j為不同的任意節點；(2)各個互感器電壓幅值量測誤差及中間量測電壓相角差的誤差均服從均值為0、方差為10-5的正態分布，即εi～N(0,10-5)，εij～N(0,10-5)，i和j為不同的任意節點。

互感器量測電壓幅值由節點電壓幅值真值疊加εi得到，中間量測節點電壓相角差量測由節點電壓相角真值做差得到相角差后，再疊加εij得到。SCADA系統的其它量測量Pi、Qi、Pij和Qij，由式(18)計算得到。

仿真采用平均絕對誤差反映估計值誤差的實際情況，平均絕對誤差(MAEk)定義如下：

(37)

圖1 IEEE-39節點系統

3.2 結果分析

假設SCADA系統數據的刷新頻率為5幀/s，即1 s采樣5個數據，采樣間隔為0.2 s。在短路故障清除后，對仿真數據進行采樣，跟蹤系統故障恢復的動態過程。為了獲得具有統計意義的結果，進行了50次蒙特卡羅模擬，取50次模擬的平均值作為最終估計結果。以節點16為例，分析動態狀態估計結果。

(1)情況1：節點16上的電壓幅值和電壓相角的估計結果分別如圖2和圖3所示。各個時刻所有節點電壓幅值和相角的平均絕對誤差比較結果分別如圖4和圖5所示。

圖2 節點16上電壓幅值估計結果比較(情況1)

圖3 節點16上電壓相角估計結果(情況1)

圖4 電壓幅值平均絕對誤差比較(情況1)

圖5 電壓相角平均絕對誤差比較(情況1)

從圖2可以看出，在故障恢復過程中，兩種方法在節點16的電壓幅值估計值均能跟隨電壓幅值的真值而變化，但文中方法的估計曲線更貼近真實曲線。從圖3可以看出，兩種方法在節點16的電壓相角估計值均能跟隨電壓相角的真值而變化，相比電壓幅值，由于相角值較大，兩種方法的估計曲線幾乎重疊，但從局部放大圖仍然可以看出文中方法相角估計曲線更貼近相角真實值。從圖4和圖5的平均絕對誤差曲線可以看出，在各個采樣時刻，無論是電壓幅值還是電壓相角，文中考慮量測相關性的方法平均絕對誤差均小于不考慮量測相關性的方法，這是由于量測誤差方差矩陣中的非對角線元素對狀態估計預測值具有修正作用。

(2)情況2：節點16上的電壓幅值和電壓相角的估計結果分別如圖6和圖7所示，各個時刻所有節點電壓幅值和相角的平均絕對誤差比較結果分別如圖8和圖9所示。

圖6 節點16上電壓幅值估計結果比較(情況2)

圖7 節點16上電壓相角估計結果(情況2)

圖8 電壓幅值平均絕對誤差比較(情況2)

圖9 電壓相角平均絕對誤差比較(情況2)

當互感器節點電壓幅值量測的誤差方差及中間量測節點電壓相角差的誤差方差由10-4變為10-5時：從圖6和圖7可以看出，文中方法在節點16上的電壓幅值和電壓相角估計值仍能較為準確地跟隨電壓幅值和相角的真值動態變化。對比圖4和圖5的平均絕對誤差曲線，從圖8和圖9可以看出，兩種方法的平均絕對誤差均有所減少，說明互感器量測值誤差大小會對估計結果產生影響，且互感器量測值誤差方差越小，估計結果越精確；在各時刻，文中考慮量測相關性方法的平均絕對誤差仍小于不考慮量測相關性的方法。

綜上可知，在電力系統動態狀態估計中考慮量測相關性是必要的；且將考慮量測相關性的CKF方法與不考慮量測相關性的CKF方法相比，前者比后者更能提高狀態估計的精度。

4 結束語

基于工程實際中SCADA量測量間存在相關性的實際情況，文中提出了一種考慮量測相關性的容積卡爾曼濾波動態狀態估計方法。在IEEE-39節點系統進行的仿真驗證了在電力系統動態狀態估計中考慮量測相關性是必要的，文中方法能夠顯著提高電力系統動態狀態估計的精度，且能為工程實際中進行動態狀態估計研究提供一定的借鑒意義。然而文中的研究僅考慮了SCADA量測相關性，隨著電力系統廣域量測裝置PMU的增多，下一步將要研究考慮混合量測相關性的電力系統動態狀態估計。