考慮量測(cè)相關(guān)性的容積卡爾曼濾波動(dòng)態(tài)狀態(tài)估計(jì)

盧慶春，張俊，許沛東，陳思遠(yuǎn)，徐箭，柯德平

(武漢大學(xué) 電氣與自動(dòng)化學(xué)院， 武漢 430072)

0 引 言

電力系統(tǒng)狀態(tài)估計(jì)用于獲得電力系統(tǒng)各節(jié)點(diǎn)母線的運(yùn)行狀態(tài)，是能量管理系統(tǒng)(Energy Management System，EMS)的重要組成部分[1]。電力系統(tǒng)狀態(tài)估計(jì)有兩種，包括靜態(tài)狀態(tài)估計(jì)及動(dòng)態(tài)狀態(tài)估計(jì)[2]。靜態(tài)狀態(tài)估計(jì)是對(duì)電力系統(tǒng)某一時(shí)刻的狀態(tài)進(jìn)行估計(jì)，而動(dòng)態(tài)狀態(tài)估計(jì)則利用當(dāng)前時(shí)刻的狀態(tài)估計(jì)結(jié)果對(duì)下一時(shí)刻狀態(tài)進(jìn)行預(yù)測(cè)，并根據(jù)下一時(shí)刻實(shí)際量測(cè)對(duì)狀態(tài)預(yù)測(cè)值進(jìn)行修正。因此，相比于靜態(tài)狀態(tài)估計(jì)，動(dòng)態(tài)狀態(tài)估計(jì)具有預(yù)測(cè)環(huán)節(jié)，其可以實(shí)時(shí)動(dòng)態(tài)地估計(jì)系統(tǒng)狀態(tài)，可以更好地對(duì)電力系統(tǒng)狀態(tài)進(jìn)行監(jiān)測(cè)。文中研究電力系統(tǒng)動(dòng)態(tài)狀態(tài)估計(jì)。

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

然而，量測(cè)量間是存在相關(guān)性的，且量測(cè)誤差協(xié)方差矩陣是非對(duì)角矩陣。文獻(xiàn)[14]通過(guò)分析變電站內(nèi)電力系統(tǒng)量測(cè)信號(hào)的采集過(guò)程驗(yàn)證了量測(cè)間存在相關(guān)性，并用點(diǎn)估計(jì)法計(jì)算量測(cè)相關(guān)性矩陣，提出了考慮量測(cè)相關(guān)性的最小二乘估計(jì)算法，仿真驗(yàn)證了考慮量測(cè)相關(guān)性的最小二乘估計(jì)算法相比傳統(tǒng)最小二乘估計(jì)算法能提高狀態(tài)估計(jì)的精度及計(jì)算效率。之后，一些文獻(xiàn)開(kāi)始研究考慮量測(cè)相關(guān)性的狀態(tài)估計(jì)。文獻(xiàn)[15]中考慮了變電站量測(cè)相關(guān)性，采用兩點(diǎn)估計(jì)法計(jì)算了變電站內(nèi)所有量測(cè)間協(xié)方差矩陣，并通過(guò)仿真驗(yàn)證考慮量測(cè)相關(guān)性能提高狀態(tài)估計(jì)的精度；文獻(xiàn)[16]分析了最小二乘算法中考慮量測(cè)相關(guān)性的重要性，通過(guò)仿真驗(yàn)證考慮量測(cè)相關(guān)性能提高狀態(tài)估計(jì)的精度。

目前量測(cè)相關(guān)性主要應(yīng)用在靜態(tài)狀態(tài)估計(jì)中，在動(dòng)態(tài)狀態(tài)估計(jì)中尚缺乏考慮量測(cè)相關(guān)性的研究。因此，文中提出了一種考慮量測(cè)相關(guān)性的容積卡爾曼濾波動(dòng)態(tài)狀態(tài)估計(jì)方法。首先對(duì)SCADA系統(tǒng)量測(cè)量存在的相關(guān)性進(jìn)行分析，然后基于雙指數(shù)平滑預(yù)測(cè)過(guò)程計(jì)算過(guò)程噪聲協(xié)方差矩陣，采用容積變換計(jì)算考慮SCADA系統(tǒng)量測(cè)相關(guān)性的量測(cè)誤差協(xié)方差矩陣，并提出考慮量測(cè)相關(guān)性的電力系統(tǒng)動(dòng)態(tài)狀態(tài)估計(jì)流程，最后在IEEE-39節(jié)點(diǎn)進(jìn)行仿真驗(yàn)證文中方法動(dòng)態(tài)估計(jì)效果。

1 電力系統(tǒng)動(dòng)態(tài)狀態(tài)估計(jì)

1.1 動(dòng)態(tài)狀態(tài)估計(jì)模型

電力系統(tǒng)動(dòng)態(tài)狀態(tài)估計(jì)以卡爾曼濾波為基礎(chǔ)， 假設(shè)非線性系統(tǒng)的狀態(tài)轉(zhuǎn)移及量測(cè)模型為：

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

(1)

zk=h(xk)+rk

(2)

式中xk為k時(shí)刻系統(tǒng)狀態(tài)向量;f(·)為狀態(tài)轉(zhuǎn)移函數(shù)；zk為k時(shí)刻系統(tǒng)量測(cè)向量；h(·)為系統(tǒng)量測(cè)函數(shù)；qk-1和rk分別為k-1時(shí)刻和k時(shí)刻的過(guò)程噪聲向量和量測(cè)誤差向量，方差分別為Qk-1和Rk。

(3)

1.2 CKF基本原理

由于CKF算法沒(méi)有EKF的線性化誤差及UKF中選擇無(wú)跡變換參數(shù)帶來(lái)的參數(shù)選取誤差，文中基于CKF進(jìn)行狀態(tài)估計(jì)。CKF狀態(tài)估計(jì)基礎(chǔ)是容積變換，容積變換指的是以一個(gè)狀態(tài)量的值作為采樣的均值點(diǎn)，依據(jù)三階球面—徑向法則[17]產(chǎn)生一些具有相同權(quán)值的容積點(diǎn)，然后將每個(gè)容積點(diǎn)按照非線性方程進(jìn)行變換，從而完成對(duì)非線性方程的近似。CKF包括預(yù)測(cè)步以及濾波步，如下所示:

(1)預(yù)測(cè)

假設(shè)k-1時(shí)刻系統(tǒng)狀態(tài)量估計(jì)誤差方差矩陣為Pk-1。對(duì)Pk-1進(jìn)行Cholesky分解，求得k-1時(shí)刻估計(jì)誤差方差陣的平方根矩陣：

(4)

(5)

(6)

根據(jù)式(5)產(chǎn)生的容積點(diǎn)，狀態(tài)量以及估計(jì)誤差方差矩陣的預(yù)測(cè)值可以按下式計(jì)算：

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

(7)

(8)

(9)

(2)濾波

zi,k=h(γi,k)

(10)

因此，量測(cè)預(yù)測(cè)值為：

(11)

預(yù)測(cè)量測(cè)自協(xié)方差矩陣為：

(12)

Tk=Sk+Rk

(13)

狀態(tài)預(yù)測(cè)及量測(cè)預(yù)測(cè)的互協(xié)方差矩陣為：

(14)

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

(15)

因此，k時(shí)刻的狀態(tài)估計(jì)值以及估計(jì)誤差方差矩陣可以表示為如下形式：

(16)

(17)

2 考慮量測(cè)相關(guān)性CKF動(dòng)態(tài)狀態(tài)估計(jì)

結(jié)合目前電網(wǎng)中廣泛配置的SCADA量測(cè)，文中研究SCADA量測(cè)間的相關(guān)性。電力系統(tǒng)動(dòng)態(tài)狀態(tài)估計(jì)的量測(cè)方程h(x)為：

(18)

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

2.1 SCADA系統(tǒng)量測(cè)相關(guān)性分析

對(duì)于SCADA系統(tǒng)，狀態(tài)估計(jì)所用的量測(cè)量包括節(jié)點(diǎn)電壓幅值、節(jié)點(diǎn)注入功率(有功/無(wú)功)、支路潮流功率(有功/無(wú)功)。除去節(jié)點(diǎn)電壓幅值，其余量測(cè)由電力系統(tǒng)中間裝置對(duì)互感器量測(cè)進(jìn)行中間計(jì)算得到。文中假設(shè)每個(gè)節(jié)點(diǎn)或線路單獨(dú)安裝了互感器。

假設(shè)互感器電壓幅值量測(cè)量為：

(19)

假設(shè)中間量測(cè)電壓相角差為：

(20)

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

(21)

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

因此，SCADA系統(tǒng)量測(cè)量間存在相關(guān)性。誤差協(xié)方差矩陣Rk是非對(duì)角陣，對(duì)角線元素表示量測(cè)量的誤差自方差，非對(duì)角線元素表示的是各量測(cè)量之間的誤差互方差。Rk可以表示為：

(22)

(23)

(24)

2.2 過(guò)程噪聲協(xié)方差矩陣Qk

過(guò)程噪聲方差矩陣Qk關(guān)系到狀態(tài)估計(jì)結(jié)果的精度，因此Qk的合適選取十分關(guān)鍵。第2.2節(jié)采用的過(guò)程如下：

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

(25)

對(duì)bk-1進(jìn)行展開(kāi)得到：

(26)

同理，也可對(duì)bk-2、bk-3及bk-n進(jìn)行展開(kāi)。

由狀態(tài)轉(zhuǎn)移方程可知，k時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)值與k時(shí)刻之前的所有時(shí)刻均有關(guān)，且受鄰近時(shí)刻的影響較大。因此，過(guò)程噪聲是動(dòng)態(tài)變化的，隨系統(tǒng)狀態(tài)的變化而時(shí)變。

由于α和β取值均為[0，1]，當(dāng)系數(shù)次數(shù)越高，取值越小。當(dāng)相鄰時(shí)刻系統(tǒng)狀態(tài)變化不大時(shí)，其差值很小，為加快求解效率，取狀態(tài)轉(zhuǎn)移方程為：

(27)

根據(jù)誤差傳遞規(guī)則，可以求得過(guò)程誤差協(xié)方差矩陣為：

(28)

2.3 考慮量測(cè)相關(guān)性的量測(cè)誤差協(xié)方差矩陣Rk

誤差協(xié)方差矩陣Rk中非對(duì)角元素體現(xiàn)著量測(cè)間的相關(guān)性，因此考慮量測(cè)相關(guān)性的狀態(tài)估計(jì)不能假設(shè)Rk為一個(gè)特定矩陣。由于容積變換方法能夠根據(jù)自變量的協(xié)方差矩陣容易得到因變量的協(xié)方差矩陣，因此，文中基于容積變換法計(jì)算Rk。

假設(shè)互感器量測(cè)的電壓幅值和相角差量測(cè)向量為：

(29)

向量維度為m1。求解互感器量測(cè)量的協(xié)方差矩陣為：

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

(30)

依照第1節(jié)中的容積變換思想，求取2m1個(gè)容積點(diǎn)ζi,k。

(31)

(32)

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

(33)

將容積點(diǎn)代入量測(cè)方程，得到SCADA系統(tǒng)的量測(cè)量如下：

Zi,k=h(ξi,k)

(34)

量測(cè)均值為：

(35)

取量測(cè)量誤差協(xié)方差矩陣近似為：

(36)

2.4 考慮量測(cè)相關(guān)性CKF動(dòng)態(tài)狀態(tài)估計(jì)流程

文中考慮量測(cè)相關(guān)性的容積卡爾曼濾波方法流程如下：

(4)根據(jù)k時(shí)刻互感器量測(cè)得到的互感器量測(cè)向量Lk，通過(guò)容積變換方法求得容積點(diǎn)，通過(guò)式(34)～式(36)計(jì)算量測(cè)誤差協(xié)方差矩陣Rk；

(6)計(jì)算狀態(tài)預(yù)測(cè)及量測(cè)預(yù)測(cè)的互協(xié)方差矩陣Ck；

(7)計(jì)算卡爾曼濾波增益Kk，k時(shí)刻狀態(tài)估計(jì)值，以及k時(shí)刻狀態(tài)協(xié)方差矩陣Pk；

(8)更新X_est，循環(huán)步驟(1)～ 步驟(8)，進(jìn)行動(dòng)態(tài)狀態(tài)估計(jì)。

3 算例分析

3.1 算例設(shè)計(jì)

為驗(yàn)證文中所提考慮量測(cè)相關(guān)性的動(dòng)態(tài)狀態(tài)估計(jì)算法的估計(jì)效果，第3.1節(jié)在IEEE-39節(jié)點(diǎn)系統(tǒng)(圖1)上進(jìn)行了仿真分析，并和不考慮量測(cè)相關(guān)性的CKF方法進(jìn)行對(duì)比。通過(guò)在DIgSILENT上進(jìn)行時(shí)域仿真產(chǎn)生一系列采樣數(shù)據(jù)，在Matlab上進(jìn)行動(dòng)態(tài)狀態(tài)估計(jì)。時(shí)域仿真時(shí)，對(duì)系統(tǒng)設(shè)置短路故障，在t=0.18 s時(shí)，節(jié)點(diǎn)16發(fā)生接地短路故障，并在t=0.36 s清除故障，時(shí)域仿真時(shí)間為10 s。

對(duì)于不考慮量測(cè)相關(guān)性的CKF方法(記為傳統(tǒng)CKF)，其量測(cè)誤差協(xié)方差矩陣的對(duì)角線元素取為文中量測(cè)誤差協(xié)方差矩陣Rk對(duì)角線元素，狀態(tài)轉(zhuǎn)移方程采用文中狀態(tài)轉(zhuǎn)移方程。

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

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

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

(37)

圖1 IEEE-39節(jié)點(diǎn)系統(tǒng)

3.2 結(jié)果分析

假設(shè)SCADA系統(tǒng)數(shù)據(jù)的刷新頻率為5幀/s，即1 s采樣5個(gè)數(shù)據(jù)，采樣間隔為0.2 s。在短路故障清除后，對(duì)仿真數(shù)據(jù)進(jìn)行采樣，跟蹤系統(tǒng)故障恢復(fù)的動(dòng)態(tài)過(guò)程。為了獲得具有統(tǒng)計(jì)意義的結(jié)果，進(jìn)行了50次蒙特卡羅模擬，取50次模擬的平均值作為最終估計(jì)結(jié)果。以節(jié)點(diǎn)16為例，分析動(dòng)態(tài)狀態(tài)估計(jì)結(jié)果。

(1)情況1：節(jié)點(diǎn)16上的電壓幅值和電壓相角的估計(jì)結(jié)果分別如圖2和圖3所示。各個(gè)時(shí)刻所有節(jié)點(diǎn)電壓幅值和相角的平均絕對(duì)誤差比較結(jié)果分別如圖4和圖5所示。

圖2 節(jié)點(diǎn)16上電壓幅值估計(jì)結(jié)果比較(情況1)

圖3 節(jié)點(diǎn)16上電壓相角估計(jì)結(jié)果(情況1)

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

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

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

(2)情況2：節(jié)點(diǎn)16上的電壓幅值和電壓相角的估計(jì)結(jié)果分別如圖6和圖7所示，各個(gè)時(shí)刻所有節(jié)點(diǎn)電壓幅值和相角的平均絕對(duì)誤差比較結(jié)果分別如圖8和圖9所示。

圖6 節(jié)點(diǎn)16上電壓幅值估計(jì)結(jié)果比較(情況2)

圖7 節(jié)點(diǎn)16上電壓相角估計(jì)結(jié)果(情況2)

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

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

當(dāng)互感器節(jié)點(diǎn)電壓幅值量測(cè)的誤差方差及中間量測(cè)節(jié)點(diǎn)電壓相角差的誤差方差由10-4變?yōu)?0-5時(shí)：從圖6和圖7可以看出，文中方法在節(jié)點(diǎn)16上的電壓幅值和電壓相角估計(jì)值仍能較為準(zhǔn)確地跟隨電壓幅值和相角的真值動(dòng)態(tài)變化。對(duì)比圖4和圖5的平均絕對(duì)誤差曲線，從圖8和圖9可以看出，兩種方法的平均絕對(duì)誤差均有所減少，說(shuō)明互感器量測(cè)值誤差大小會(huì)對(duì)估計(jì)結(jié)果產(chǎn)生影響，且互感器量測(cè)值誤差方差越小，估計(jì)結(jié)果越精確；在各時(shí)刻，文中考慮量測(cè)相關(guān)性方法的平均絕對(duì)誤差仍小于不考慮量測(cè)相關(guān)性的方法。

綜上可知，在電力系統(tǒng)動(dòng)態(tài)狀態(tài)估計(jì)中考慮量測(cè)相關(guān)性是必要的；且將考慮量測(cè)相關(guān)性的CKF方法與不考慮量測(cè)相關(guān)性的CKF方法相比，前者比后者更能提高狀態(tài)估計(jì)的精度。

4 結(jié)束語(yǔ)

基于工程實(shí)際中SCADA量測(cè)量間存在相關(guān)性的實(shí)際情況，文中提出了一種考慮量測(cè)相關(guān)性的容積卡爾曼濾波動(dòng)態(tài)狀態(tài)估計(jì)方法。在IEEE-39節(jié)點(diǎn)系統(tǒng)進(jìn)行的仿真驗(yàn)證了在電力系統(tǒng)動(dòng)態(tài)狀態(tài)估計(jì)中考慮量測(cè)相關(guān)性是必要的，文中方法能夠顯著提高電力系統(tǒng)動(dòng)態(tài)狀態(tài)估計(jì)的精度，且能為工程實(shí)際中進(jìn)行動(dòng)態(tài)狀態(tài)估計(jì)研究提供一定的借鑒意義。然而文中的研究?jī)H考慮了SCADA量測(cè)相關(guān)性，隨著電力系統(tǒng)廣域量測(cè)裝置PMU的增多，下一步將要研究考慮混合量測(cè)相關(guān)性的電力系統(tǒng)動(dòng)態(tài)狀態(tài)估計(jì)。