非高斯噪聲下基于CKMC-CKF的發(fā)電機動態(tài)狀態(tài)估計

趙志雷，黃繼東，王義，楊志偉

(1.國能神東煤炭集團(tuán)有限責(zé)任公司供電中心，陜西 榆林 719300；2.許繼電氣股份有限公司，河南 許昌 461000；3.鄭州大學(xué) 電氣與信息工程學(xué)院，河南 鄭州 450001；4. 國網(wǎng)三門峽供電公司，河南 三門峽 472000)

隨著相量測量單元(phasor measurement unit，PMU)的廣泛應(yīng)用，不同類型的狀態(tài)估計(state estimation，SE)技術(shù)被用于跟蹤電力系統(tǒng)的隱藏狀態(tài)變量[1-3]。靜態(tài)SE假定系統(tǒng)處于準(zhǔn)穩(wěn)態(tài)狀態(tài)下，可以監(jiān)測某時刻母線電壓幅值和相位角。動態(tài)SE(dynamic state estimation，DSE)能實時監(jiān)測系統(tǒng)的動態(tài)變化[4-5]。因此，DSE被廣泛應(yīng)用于電氣設(shè)備控制方案的實現(xiàn)，如同步發(fā)電機的機電暫態(tài)跟蹤。

對于DSE，近年來已提出了諸多經(jīng)典卡爾曼濾波器(Kalman filter，KF)，如擴展KF(extended KF，EKF)[6]、無跡KF(unscented KF，UKF)[7]、容積KF(cubature KF，CKF)[8]及其擴展形式[9-10]，旨在準(zhǔn)確跟蹤系統(tǒng)的內(nèi)部數(shù)據(jù)。鑒于非高斯過程的不確定性可能導(dǎo)致的偏差估計，文獻(xiàn)[11]提出了一種簡化的高斯綜合EKF方法來處理非線性狀態(tài)估計。考慮到線性化近似引起的EKF截斷誤差，研究人員又開發(fā)了幾種無導(dǎo)數(shù)DSE方法。文獻(xiàn)[12]基于平方根UKF和作用于測量狀態(tài)估計的加權(quán)因子提出一種改進(jìn)的方法，提高了狀態(tài)估計的穩(wěn)定性和魯棒性。文獻(xiàn)[13]引入了自適應(yīng)UKF，同時監(jiān)控集成電機變速器的控制輸入和狀態(tài)信息。文獻(xiàn)[14]開發(fā)了基于變分貝葉斯的自適應(yīng)CKF，可用于跟蹤隨機發(fā)生的注入攻擊的狀態(tài)。上述所提出的狀態(tài)估計方法大多基于高斯噪聲分布假設(shè)，在存在非高斯噪聲的情況下，這些方法的估計性能可能會顯著下降。

為了減輕測量函數(shù)中異常值的不利影響，如非高斯噪聲、負(fù)載波動[15]等，研究者也提出了很多改進(jìn)措施。為了降低惡意網(wǎng)絡(luò)攻擊導(dǎo)致的估計誤差，文獻(xiàn)[16]中構(gòu)建了基于分布式壓縮感知的估計策略。此外，作為信息理論學(xué)習(xí)(information-theoretic learning，ITL)的局部相似函數(shù)，熵包含了概率密度函數(shù)的高階矩，在處理非高斯噪聲假設(shè)時具有優(yōu)異的特性。其中，基于最大熵準(zhǔn)則的DSE算法被融合到EKF[17]、UKF[18]中，大幅增強了卡爾曼濾波器在非高斯噪聲下的魯棒性。然而，基于高斯核相關(guān)熵的最優(yōu)核帶寬(kernel bandwidth，KB)往往是通過大量計算和測試獲得的，這在實際應(yīng)用中限制了該方法的通用性。此外，高斯核函數(shù)在矩陣Cholesky分解中會出現(xiàn)奇異值[19]，嚴(yán)重影響了DSE的計算精度。

為了解決上述問題，提高傳統(tǒng)DSE對非高斯噪聲問題的魯棒性，本文基于柯西核相關(guān)熵準(zhǔn)則，提出一種基于柯西核最大相關(guān)熵(Cauchy kernel maximum correntropy，CKMC)的CKF算法(簡稱CKMC-CKF算法)。CKMC-CKF算法通過引入柯西內(nèi)核優(yōu)化準(zhǔn)則，將CKF算法和CKMC準(zhǔn)則相結(jié)合，可以有效抑制同步發(fā)電機模型中非高斯噪聲對測量數(shù)據(jù)的不利影響。在IEEE 39節(jié)點系統(tǒng)上進(jìn)行仿真分析，以驗證CKMC-CKF算法在非高斯噪聲條件下的跟蹤性能。

1 發(fā)電機動態(tài)狀態(tài)估計模型

實際電力系統(tǒng)的離散微分方程由狀態(tài)變量和測量值2個時變離散非線性函數(shù)組成，采用改進(jìn)歐拉法[19]實現(xiàn)，非線性函數(shù)的具體表達(dá)式為：

(1)

式中：φ(·)、h(·)分別為狀態(tài)傳播函數(shù)和狀態(tài)測量函數(shù)；Xk、Zk分別為時刻k的狀態(tài)向量和觀測向量；uk為時刻k的輸入向量；wk、vk分別為時刻k的狀態(tài)噪聲和測量噪聲。

同步發(fā)電機是電力系統(tǒng)中關(guān)鍵設(shè)備之一，對其準(zhǔn)確建模對電力系統(tǒng)的動態(tài)跟蹤具有重要意義。本文采用同步發(fā)電機的四階模型，比傳統(tǒng)的二階發(fā)電機模型更接近發(fā)電機的實際運行模式。發(fā)電機狀態(tài)函數(shù)建立如下[20]：

δ=ω-ω0，

(2)

(3)

(4)

(5)

式中：δ為發(fā)電機的功率角；ω為電角速度，ω0為其初始值；ed、eq分別為d、q軸上的瞬時電動勢(下標(biāo)d、q分別表示d、q軸，下同)；KD為阻尼系數(shù)；Aj、Tm、Te分別為慣性常數(shù)、機械功率、電磁功率；Ufd為定子的勵磁電壓；xd、xq為同步電抗；x′d、x′q為瞬態(tài)電抗；Td0、Tq0為時間常數(shù)；id、iq為定子電流。

本研究中設(shè)置狀態(tài)向量Xk=(δ，ω，eq，ed)T，輸入向量uk=(Tm，Ufd，iR，iI)T，其中iR、iI分別為定子R軸和I軸上的電流。為了提高同步發(fā)電機的估計精度，采用四維量測量，發(fā)電機的測量函數(shù)設(shè)置為Zk=(δ，ω，eR，eI)T，其中：

eR=(ed+idxq)sinδ+(eq-idxd)cosδ，

(6)

eI=(eq-idxd)sinδ-(ed+iqxq)cosδ.

(7)

為便于求解，將式(6)、(7)表示為狀態(tài)變量和輸入變量的函數(shù)：

id=iRsinδ-iIcosδ，

(8)

iq=iIsinδ+iRcosδ.

(9)

2 基于CKMC準(zhǔn)則的CKF魯棒動態(tài)估計

2.1 CKMC

作為統(tǒng)計相似性的度量，相關(guān)熵可以由1對隨機變量A和B定義[17]：

V(A，B)=E[κσ(A，B)]=

(10)

式中：V(A，B)為變量A和B的均值；E(·)為期望算子，κσ(A，B)為對應(yīng)的Mercer核，σ為KB；fA，B(A，B)為A和B之間的聯(lián)合概率密度函數(shù)。

由于樣本的有限性，V(A，B)的近似估計

(11)

式中：N為樣本總數(shù)；Ai、Bi為第i對隨機變量。

在信息理論學(xué)習(xí)(information theory learning，ITL)中，高斯核通常被認(rèn)為是相關(guān)熵的核函數(shù)。本文采用柯西核作為熵的核函數(shù)，與高斯核[16]相比，柯西核結(jié)構(gòu)簡單，對KB不敏感，避免了最優(yōu)KB選擇的大量測試，提高了計算效率[17]?？挛骱撕瘮?shù)定義為

(12)

式中：e=A-B，Cσ(e)為對應(yīng)的柯西核。

通過最大化式(11)并結(jié)合式(12)，可以獲得基于CKMC的目標(biāo)函數(shù)

(13)

2.2 CKMC-CKF方法的工作機理

本文所提的CKMC-CKF是在傳統(tǒng)CKF框架上進(jìn)行改進(jìn)，通過線性回歸方程推導(dǎo)出2種加權(quán)局部相似度函數(shù)，并結(jié)合CKMC準(zhǔn)則，利用固定點迭代法不斷更新誤差協(xié)方差矩陣，從而降低不良數(shù)據(jù)的權(quán)重，最終提高動態(tài)狀態(tài)估計的精度和魯棒性。具體原理包含如下過程：

(14)

(15)

當(dāng)k>0時，首先，根據(jù)球面徑向容積法則，生成一組容積Sigma點：

(16)

j=1，2，…，2n.

(17)

(18)

式中：Xj，k-1|k-1為k-1時刻第j個Sigma點；Ij為單位矩陣；n為狀態(tài)維度；Sk-1|k-1表示對矩陣Pk-1|k-1Cholesky分解。

(19)

(20)

(21)

(22)

(23)

Zi，k∣k-1=h(Xk，Xi，k∣k-1)；

(24)

(25)

(26)

式中：Xi，k∣k-1為Sigma點；Zi，k∣k-1為經(jīng)量測方程傳播后的Sigma點。

c)線性回歸方程的推導(dǎo)。為了迭代計算出最佳的狀態(tài)估計值，首先定義偽量測向量

(27)

將量測函數(shù)線性化近似，得到近似后的量測值

(28)

由式(1)、(28)，可得線性回歸方程組

(29)

其中，

(30)

(31)

式(29)—(31)中：ηk為誤差矩陣；Bp，k∣k-1、Br，k分別為Pk∣k-1、Rk進(jìn)行Cholesky分解所得的分解因子。

基于式(13)，所提CKMC-CKF方法的目標(biāo)函數(shù)可進(jìn)一步表示為：

(32)

(33)

(34)

式中：eX，i、eZ，i為加權(quán)矩陣eX、eZ的第i個變量；m為量測信息的維度數(shù)。

柯西核相關(guān)熵采用2種加權(quán)相似函數(shù)對狀態(tài)誤差和量測誤差進(jìn)行修正，從而降低不良數(shù)據(jù)的權(quán)重。對式(32)求偏導(dǎo)，并且令?JCKMC/?Xk=0，可得：

(35)

進(jìn)一步整理式(35)，利用量測量對預(yù)測后的狀態(tài)值進(jìn)行修正，得到狀態(tài)量的最優(yōu)估計形式為

(36)

其中：

(37)

(38)

(39)

(40)

(41)

最后，狀態(tài)量后驗誤差協(xié)方差Pk更新為

(42)

基于CKMC-CKF的同步發(fā)電機DSE流程如圖1所示，其中，t為內(nèi)循環(huán)的時間變量，k為外循環(huán)的時間變量，ε為判斷閾值。

圖1 基于CKMC-CKF的同步發(fā)電機DSE流程Fig.1 Synchronous generator DSE flowchart based on CKMC-CKF

3 算例分析

3.1 算例設(shè)計與評價指標(biāo)

在本節(jié)中，使用IEEE 10機39節(jié)點標(biāo)準(zhǔn)系統(tǒng)模擬電力系統(tǒng)的網(wǎng)絡(luò)結(jié)構(gòu)，其拓?fù)淙鐖D2所示。由于噪聲分布的隨機性，使用蒙特卡洛技術(shù)來模擬噪聲的統(tǒng)計信息，并設(shè)置采樣數(shù)N=200。假設(shè)系統(tǒng)在16-21線路處發(fā)生三相金屬性接地短路，且當(dāng)仿真進(jìn)行到0.5 s時出現(xiàn)系統(tǒng)故障，經(jīng)0.2 s后系統(tǒng)過渡至新穩(wěn)態(tài)。其中，以穩(wěn)態(tài)運行值作為狀態(tài)變量的初始值，并且初始誤差協(xié)方差設(shè)置為P0∣0=10-5I。此外，過程和測量的噪聲協(xié)方差矩陣分別設(shè)置為10-5I和10-6I。本文以同步發(fā)電機Gen2為例，將CKF[8]、MCC-CKF[22]與所提CKMC-CKF方法進(jìn)行比較，驗證CKMC-CKF方法的有效性。

圖2 IEEE 39節(jié)點系統(tǒng)Fig.2 IEEE 39 bus system

為了評估各種濾波方法對于發(fā)電機狀態(tài)的跟蹤效果，本文將平均絕對誤差M和總體性能誤差J作為評價指標(biāo)對不同方法的估計結(jié)果進(jìn)行量化分析，定義如下[17]：

(43)

(44)

3.2 KB測試

如文獻(xiàn)[16]所述，相關(guān)熵中KB對狀態(tài)估計的準(zhǔn)確性具有決定性作用。具體而言，過大的KB可能無法抑制異常值，過小的KB可能導(dǎo)致收斂緩慢，因此KB的取值對于DSE的估計結(jié)果至關(guān)重要。為了驗證所提方法在不同KB下的性能，本文所提出的CKMC-CKF算法在不同KB下的誤差指標(biāo)見表1。

表1 不同KB下CKMC-CKF的誤差Tab.1 Errors of CKMC-CKF at different KBs

從表1數(shù)據(jù)可知，CKMC-CKF算法在不同KB下的估計誤差沒有太大差異?；诳挛鲀?nèi)核損失準(zhǔn)則的最大熵對KB的選擇不敏感，在很大程度上避免了對最優(yōu)KB的大量實驗，從而提高了算法的通用性。為了便于后續(xù)研究，CKMC-CKF的KB設(shè)置為50。

3.3 噪聲不確定性測試

由于實際電網(wǎng)中運行環(huán)境的時變性，噪聲的統(tǒng)計特性在實際系統(tǒng)中會隨著時間而變化[21]，因此，實際仿真中的系統(tǒng)噪聲通常是未知的。為了驗證所提方法在該工況下的有效性，假設(shè)系統(tǒng)的過程噪聲協(xié)方差設(shè)為Qk=diag(10-3，10，10-1，10-1)，測量噪聲設(shè)為Rk=diag(10-5，10-5，10-3，10-5)。

在該工況下，CKF、MCC-CKF、CKMC-CKF 3種方法的狀態(tài)變量估計結(jié)果如圖3所示?？梢钥闯觯寒?dāng)故障出現(xiàn)后，CKF對于同步發(fā)電機的狀態(tài)跟蹤出現(xiàn)了較大的估計偏差，精度下降嚴(yán)重；MCC-CKF對不確定的高斯噪聲具有一定的抵抗力，但也因其對KB的依賴而受到限制。相比而言，CKMC-CKF具有良好的數(shù)值穩(wěn)定性和對不同KB的不敏感性，因此其估計曲線最貼近系統(tǒng)的真實狀態(tài)，狀態(tài)估計的精度更高。

圖3 噪聲不確定工況下的狀態(tài)估計Fig.3 State estimation under noise uncertainty

噪聲不確定工況下CKF、MCC-CKF和CKMC-CKF方法的估計誤差指標(biāo)對比情況見表2?？梢钥闯?，與傳統(tǒng)CKF和MCC-CKF相比，本文所提的CKMC-CKF方法誤差最小，估計精度更高，對于不確定噪聲的魯棒性更強。

表2 噪聲不確定工況下的估計誤差指標(biāo)Tab.2 Estimation error metrics under noise uncertainties

3.4 非高斯噪聲測試

發(fā)電機動態(tài)估計模型中通常假設(shè)噪聲為高斯白噪聲，即噪聲協(xié)方差是一個常數(shù)矩陣。在實際電力系統(tǒng)中，PMU采集的量測信息在傳送至控制中心時，通信信道易遭受噪聲污染，導(dǎo)致出現(xiàn)非高斯噪聲，這會影響動態(tài)狀態(tài)估計性能，降低濾波精度[22]。

為了驗證所提方法在非高斯噪聲下的濾波性能，假設(shè)量測方程中的噪聲信號rk是由混合高斯函數(shù)表示的重尾噪聲，即

(45)

圖4展示了v1=10-4、v2=10-6、θ=0.98時的狀態(tài)估計情況?？梢姡?dāng)系統(tǒng)出現(xiàn)非高斯噪聲時，基于高斯噪聲假設(shè)的CKF方法跟蹤速度大幅下降，與狀態(tài)真實值具有明顯的誤差，跟蹤精度較低，已經(jīng)不滿足動態(tài)估計的準(zhǔn)確性和實時性要求；相比于CKF，MCC-CKF方法對非高斯噪聲具有一定的魯棒性，但是由于MCC準(zhǔn)則基于高斯內(nèi)核，對于非高斯噪聲情形的估計誤差仍然較大；CKMC-CKF方法與真實值更加接近，相比于CKF和MCC-CKF，估計精度分別提高了約78%和35%。這是由于CKMC-CKF方法將柯西內(nèi)核作為相關(guān)熵的內(nèi)核函數(shù)，既克服了KB選擇復(fù)雜性的困難，也極大降低了非高斯噪聲情形下發(fā)電機的估計誤差。因此，該工況下CKMC-CKF方法具有最高的估計精度和最強的魯棒性能，能夠準(zhǔn)確地跟蹤發(fā)電機的狀態(tài)。

圖4 非高斯噪聲工況下的狀態(tài)估計Fig.4 State estimation under non-Gaussian noise

不同濾波器在混合程度為θ=0.98、θ=0.95下的誤差指標(biāo)分別見表3、表4。可以看出，本文方法在非高斯噪聲下估計誤差最小，濾波性能最好。

表3 非高斯噪聲混合度為0.98時的誤差指標(biāo)Tab.3 Error metrics for non-Gaussian noise with mixing degrees θ=0.98

表4 非高斯噪聲混合度為0.95時的誤差指標(biāo)Tab.4 Error metrics for non-Gaussian noise with mixing degrees θ=0.95

4 結(jié)論

針對電力系統(tǒng)中同步發(fā)電機動態(tài)狀態(tài)估計中存在的非高斯噪聲問題，本文提出了一種基于CKMC-CKF的發(fā)電機動態(tài)狀態(tài)估計方法，主要結(jié)論如下：

a)利用柯西核作為相關(guān)熵的核函數(shù)，避免了最優(yōu)KB選擇的大量測試，提高了動態(tài)狀態(tài)估計的計算效率；利用熵增益動態(tài)修正誤差協(xié)方差，增強了狀態(tài)估計的精確性。

b)通過在IEEE 39節(jié)點系統(tǒng)進(jìn)行仿真分析可知，相比于CKF和MCC-CKF，所提CKMC-CKF方法在非高斯噪聲環(huán)境下具有更準(zhǔn)確的估計精度和更強的魯棒性，可以有效地跟蹤同步發(fā)電機的狀態(tài)。

c)后續(xù)將進(jìn)一步研究電力系統(tǒng)遭受網(wǎng)絡(luò)攻擊情況時，同步發(fā)電機的量測信息和狀態(tài)估計問題。