張 明,陸東亮,徐詩露,夏若平,何順帆
(1.武漢紡織大學(xué)電子與電氣工程學(xué)院,湖北武漢 430200;2.中南民族大學(xué)計(jì)算機(jī)科學(xué)學(xué)院,湖北武漢 430074)
諧波治理一直都是電力系統(tǒng)中的關(guān)鍵問題。隨著電網(wǎng)中大量可再生能源的接入,加上非線性負(fù)載的廣泛使用,電網(wǎng)結(jié)構(gòu)異常復(fù)雜并由此產(chǎn)生了大量諧波[1-2]。因此,如何有效治理諧波成為了供用電雙方重點(diǎn)關(guān)注的問題,需要不斷提升電力諧波信號(hào)的檢測(cè)精度[3]。
對(duì)電力系統(tǒng)諧波進(jìn)行快速準(zhǔn)確的檢測(cè),目前人們已經(jīng)提出了很多檢測(cè)諧波的方法[4],從諧波檢測(cè)的狀態(tài)分類,主要可分為穩(wěn)態(tài)和動(dòng)態(tài)諧波檢測(cè),其中穩(wěn)態(tài)諧波比較常見也比較容易進(jìn)行治理,而動(dòng)態(tài)諧波僅出現(xiàn)于某種運(yùn)行工況之下,檢測(cè)的難度較高[5-6],故動(dòng)態(tài)諧波的檢測(cè)算法是本文研究的重點(diǎn)。各種諧波檢測(cè)的文獻(xiàn)表明卡爾曼濾波(Kalman Filter,KF)及其擴(kuò)展算法都可以檢測(cè)動(dòng)態(tài)諧波,文獻(xiàn)[7]首次將局部集成變換的卡爾曼濾波器算法應(yīng)用于電力系統(tǒng)諧波檢測(cè),能夠提高動(dòng)態(tài)諧波的檢測(cè)精度,但是其要求處理的信號(hào)具有線性關(guān)系。文獻(xiàn)[8-9]提出了一種擴(kuò)展卡爾曼濾波(Extended Kalman Filter,EKF)算法,這是一種局部方法,能夠?yàn)榉蔷€性系統(tǒng)提供次優(yōu)解,但是這種方法可能會(huì)因?yàn)槟P偷姆蔷€性化而導(dǎo)致結(jié)果出現(xiàn)偏差。文獻(xiàn)[10]在EKF 的基礎(chǔ)上引入了無跡卡爾曼濾波(Unscented Kalman Filter,UKF)算法,它減少了預(yù)測(cè)狀態(tài)的線性負(fù)擔(dān),能夠?qū)⒕忍嵘炼A等級(jí),但是在處理二階函數(shù)時(shí)仍然不能給出正確的二階時(shí)刻。文獻(xiàn)[11]提出了一種容積卡爾曼(Cubature Kalman Filter,CKF)算法來提高諧波的檢測(cè)精度,其利用了三階球面徑向規(guī)則,能夠解決高維狀態(tài)空間問題,且由于其采用的非線性模型,能夠克服線性化問題,其但是建立在噪聲矩陣不變的情況下。在動(dòng)態(tài)諧波檢測(cè)中,不僅要考慮諧波信號(hào)的非線性化問題,還需要考慮高維問題以及因?yàn)橄到y(tǒng)噪聲引起的檢測(cè)誤差,因此,本文在CKF的基礎(chǔ)上通過引入噪聲估值的遺忘因子,提出了一種自適應(yīng)容積卡爾曼濾波(Adaptive Cubature Kalman Filter,ACKF)算法,仿真結(jié)果表明本文提出的算法在動(dòng)態(tài)諧波檢測(cè)中具有更高的精度。
卡爾曼濾波是一種基于時(shí)域的濾波方法,將某一時(shí)刻測(cè)量的數(shù)據(jù)作為量測(cè)量來預(yù)測(cè)下一時(shí)刻的狀態(tài)量。設(shè)k時(shí)刻一組隨機(jī)的離散線性方程為[12-13]:
式中:Xk為系統(tǒng)的狀態(tài)向量;Ak-1為k-1 時(shí)刻的狀態(tài)轉(zhuǎn)移矩陣;Wk-1為k-1 時(shí)刻的系統(tǒng)噪聲向量;Zk為系統(tǒng)的量測(cè)向量;Hk為系統(tǒng)的量測(cè)矩陣;Vk為系統(tǒng)的量測(cè)噪聲向量。
假設(shè)系統(tǒng)中出現(xiàn)的噪聲都是高斯白噪聲,則其均值滿足為0,系統(tǒng)噪聲方差矩陣為Q,量測(cè)噪聲方差矩陣為R。
由式(1)可知,當(dāng)前狀態(tài)的預(yù)測(cè)X′k-1是根據(jù)當(dāng)前時(shí)刻測(cè)量的數(shù)據(jù)得來的,可以寫成關(guān)系式(2),而下一時(shí)刻的預(yù)測(cè)又是在當(dāng)前預(yù)測(cè)的基礎(chǔ)上做出的,所以可以得到式(3)。
因?yàn)橄到y(tǒng)噪聲只對(duì)狀態(tài)向量產(chǎn)生影響,所以忽略E[Wk-1/Z1Z2…Zk-1],E表示單位向量,由式(2)和式(3)就可以得到預(yù)測(cè)狀態(tài)向量:
而在實(shí)際應(yīng)用中,預(yù)測(cè)的狀態(tài)量與狀態(tài)量的真實(shí)值總存在一定的誤差,由式(1)可知,狀態(tài)量產(chǎn)生誤差也將會(huì)影響量測(cè)值Zk,所以就有如下推導(dǎo):
式中:ΔX′k為預(yù)測(cè)狀態(tài)量與真實(shí)狀態(tài)量之間的誤差。
式中:Z′k為預(yù)測(cè)狀態(tài)向量出現(xiàn)誤差下的測(cè)量值;B為Zk與Z′k之間的殘差值。
由式(6)可以看出,通過處理殘差可以達(dá)到降低預(yù)測(cè)誤差的目的,引入濾波增益Kk,則可以得到修正后的狀態(tài)向量為:
在預(yù)測(cè)的過程中求得合適的濾波增益系數(shù)是本算法的關(guān)鍵。
卡爾曼濾波是一種最小均方誤差估計(jì),在修正后的狀態(tài)向量X″k的基礎(chǔ)上繼續(xù)求得差值記為ΔX″k=Xk-X″k。
則結(jié)合式(1)、式(5)和式(7)可得:
于是有:
因?yàn)榱繙y(cè)噪聲Vk不會(huì)對(duì)狀態(tài)向量產(chǎn)生影響,所以二者滿足不相關(guān),就有[14]:
式中:Pk為系統(tǒng)的偏差矩陣,使其達(dá)到最小值則可以求得Kk的值:
然后將式(13)代入式(11)中可得系統(tǒng)的偏差矩陣為:
最后由狀態(tài)向量誤差推導(dǎo)出預(yù)測(cè)偏差矩陣為:
式中:Qk-1為k-1時(shí)刻系統(tǒng)噪聲的方差矩陣。
傳統(tǒng)的卡爾曼濾波算法能夠預(yù)測(cè)出動(dòng)態(tài)諧波信號(hào)的相關(guān)信息,但它是建立在一種具有線性關(guān)系的基礎(chǔ)上,具有一定的局限性,需要進(jìn)一步的分析研究。
傳統(tǒng)的卡爾曼濾波只能對(duì)線性關(guān)系的狀態(tài)方程進(jìn)行處理,在實(shí)際中是無法運(yùn)用的,因此需要在傳統(tǒng)方法的基礎(chǔ)上進(jìn)行改進(jìn),針對(duì)此問題,提出了一種基于自適應(yīng)的容積卡爾曼濾波算法[15]。
對(duì)于高斯分布下的非線性濾波問題,實(shí)際上就是求解狀態(tài)量的后驗(yàn)概率密度函數(shù)乘積的積分值,CKF 可以看作是一種利用三階球面徑向規(guī)則生成一組等權(quán)值的容積點(diǎn)來求解后經(jīng)驗(yàn)期望的積分[16]。傳統(tǒng)的卡爾曼濾波的一般積分形式可以表示為:
式中:G(f)為待求的積分;Rn為n維的積分域;f(x)為非線性積分;x為狀態(tài)向量。
在容積卡爾曼濾波算法中需要將一般積分形式變換為球面徑向積分形式和三階球面徑向積分形式,經(jīng)過對(duì)式(16)的轉(zhuǎn)化,球面徑向積分形式可以表示為:
由式(17)可以發(fā)現(xiàn),該變換是將狀態(tài)向量x的積分分解為半徑r與方向向量y相關(guān)的二重積分,式(17)中x=ry,r∈[0,+∞),令yTy=1,因此可得xTx=r2,可以看出變換后積分中的球面單位為Un={y∈Rn|yTy=1},σ(·)表示積分域的積分微元。
對(duì)式(17)進(jìn)一步轉(zhuǎn)化可得:
式中:S(r)為球面積分,S(r)=∫Un f(ry)dσ(y)
在完成了球面徑向積分形式的轉(zhuǎn)化后,對(duì)式(18)和S(r)使用三階球面準(zhǔn)則,可得mr×ms的球面徑向容積準(zhǔn)則為:
式中:mr=1,ms=2n表示容積點(diǎn)的點(diǎn)數(shù);{ai,ri} 為計(jì)算徑向積分容積點(diǎn)集合;{bj,yj} 為計(jì)算球面積分容積點(diǎn)集合。
為了更好地說明容積卡爾曼濾波算法在非線性動(dòng)態(tài)諧波檢測(cè)上的應(yīng)用,設(shè)一個(gè)隨機(jī)的非線性離散方程為:
式中:xk為k時(shí)刻的狀態(tài)量;zk為k時(shí)刻的量測(cè)量;wk,vk分別為系統(tǒng)噪聲和量測(cè)噪聲,假設(shè)其為高斯白噪聲,則其方差也滿足高斯白噪聲,即wk~(0,Qk),vk~(0,Rk),R,Q之間互不相關(guān),且其均值為0。
容積卡爾曼濾波算法的實(shí)現(xiàn)可以分為初始化、狀態(tài)預(yù)測(cè)和狀態(tài)校正3 個(gè)步驟[17-19]。
2.1.1 初始化
在利用CKF 進(jìn)行動(dòng)態(tài)諧波狀態(tài)估計(jì)時(shí)需要給出前一時(shí)刻的量測(cè)量,本文給定k時(shí)刻的狀態(tài)量和估計(jì)偏差矩陣分別為xk和Pk|k。
當(dāng)k=0 時(shí)刻時(shí),假設(shè)x0~(xˉ0,P0),Q0,R0,則得濾波器的最優(yōu)初始化為:
2.1.2 狀態(tài)預(yù)測(cè)
狀態(tài)預(yù)測(cè)就是根據(jù)球面徑向規(guī)則生成的容積點(diǎn)集結(jié)合狀態(tài)方程來估計(jì)下一時(shí)刻的狀態(tài)量,并能夠計(jì)算狀態(tài)量的偏差矩陣,據(jù)此可以得到關(guān)于容積點(diǎn)的表達(dá)式如下:
式中:ξj為容積點(diǎn),,其中[L]為n維空間的點(diǎn)集矩陣如式(23)所示;[L]j為第j個(gè)容積點(diǎn);j取為權(quán)值;m為容積點(diǎn)的個(gè)數(shù),在三階球面徑向規(guī)則下滿足m=2n。
根據(jù)預(yù)測(cè)的本質(zhì),通過k時(shí)刻的估計(jì)偏差矩陣Pk|k計(jì)算諧波狀態(tài)的等權(quán)值容積點(diǎn)。
在得到了諧波的等權(quán)值容積點(diǎn)后,結(jié)合離散方程就可以估計(jì)出下一時(shí)刻系統(tǒng)的狀態(tài)向量,其表達(dá)式如下:
由容積點(diǎn)的非線性傳播和進(jìn)一步的轉(zhuǎn)化可以計(jì)算得到諧波在k+1 時(shí)刻的估計(jì)偏差矩陣:
2.1.3 狀態(tài)校正
1)結(jié)合式(27)可計(jì)算諧波狀態(tài)向量測(cè)量的等權(quán)值容積點(diǎn):
將其代入離散方程式(20)中就可以得到測(cè)量的量測(cè)向量。
再由式(26)就可以估計(jì)預(yù)測(cè)的諧波狀態(tài)測(cè)量值:
2)將估計(jì)得到的和實(shí)際測(cè)量得到的諧波狀態(tài)量對(duì)比,可以得到誤差系數(shù),求得修正后的諧波狀態(tài)向量和諧波估計(jì)偏差矩陣。
則可得:
式中:Pvv,k+1為更新后的量測(cè)量估計(jì)偏差矩陣;Pvz,k+1|k為估計(jì)的交叉協(xié)方差矩陣。
最后由更新后的式(31)和式(32)計(jì)算得到濾波增益系數(shù)Kk+1,就可以得到k+1 時(shí)刻估計(jì)的諧波狀態(tài)向量均值和偏差矩陣Pk+1。
容積卡爾曼濾波諧波動(dòng)態(tài)估計(jì)可以得到很好的估計(jì)效果,但其是建立在電力系統(tǒng)中系統(tǒng)噪聲和量測(cè)噪聲都保持不變的前提下,而實(shí)際中的協(xié)方差矩陣Q和R都是在不斷變化的,這也就使得容積卡爾曼濾波只能停留在理想狀況下,不能捕捉到估計(jì)噪聲協(xié)方差所需的整體噪聲動(dòng)態(tài)。因此,在容積卡爾曼濾波的基礎(chǔ)上加入了基于漸消記憶指數(shù)加權(quán)法的噪聲估值算法來得到一種自適應(yīng)容積卡爾曼濾波(ACKF)算法,以自適應(yīng)系統(tǒng)不斷變化的噪聲[20-22]。
基于漸消記憶指數(shù)加權(quán)法引入了噪聲估值遺忘因子α,這個(gè)因子決定了對(duì)過去的噪聲協(xié)方差估計(jì)和當(dāng)前的噪聲瞬時(shí)估計(jì)的權(quán)重大小,并且可以提供噪聲協(xié)方差的過去和現(xiàn)在估計(jì)的加權(quán)平均值。而當(dāng)遺忘因子為所有過去的估計(jì)提供恒定的權(quán)重時(shí),可能無法跟蹤系統(tǒng)噪聲指標(biāo)的變化[23]。因此,引入了衰落記憶的概念,通過不斷遺忘過去的陳舊數(shù)據(jù)而加大對(duì)現(xiàn)在估計(jì)生成的噪聲數(shù)據(jù),以指數(shù)衰減的形式加權(quán)過去的估計(jì),通過式(34)確定加權(quán)系數(shù),加權(quán)系數(shù)由遺忘因子用于確定每個(gè)延遲的權(quán)重估計(jì),過去的權(quán)重往往呈指數(shù)衰減從而產(chǎn)生衰減記憶。在動(dòng)態(tài)諧波狀態(tài)估計(jì)中,由于新生成的噪聲數(shù)據(jù)與下一時(shí)刻的噪聲數(shù)據(jù)具有相似性,因此,漸消記憶指數(shù)加權(quán)法就是通過遺忘因子確定前后時(shí)刻的加權(quán)系數(shù),由加權(quán)系數(shù)再結(jié)合前面的離散方程就可以得到新的系統(tǒng)噪聲和量測(cè)噪聲的協(xié)方差矩陣。
式中:dk為加權(quán)系數(shù);α取0.85;Q′k+1為k+1 時(shí)刻生成的系統(tǒng)噪聲的協(xié)方差矩陣;R′k+1為k+1 時(shí)刻生成的量測(cè)噪聲的協(xié)方差矩陣,其中εk=Zk-h(xk),表示殘差值。
在原始的算法中加入基于漸消記憶指數(shù)加權(quán)法的噪聲估值器,得到自適應(yīng)的容積卡爾曼濾波算法,其步驟為:
1)采用漸消記憶指數(shù)加權(quán)法生成系統(tǒng)噪聲和量測(cè)噪聲的時(shí)變協(xié)方差矩陣分別表示為Q′k+1和R′k+1;
2)用修正后的Q′k+1替換式(27)中的Qk+1得到k+1 時(shí)刻諧波狀態(tài)量的估計(jì)偏差矩陣,然后由式(28)和式(30)分別計(jì)算諧波量測(cè)量的等權(quán)容積點(diǎn)和諧波量測(cè)量,進(jìn)而就可以求得在時(shí)變?cè)肼曄碌南到y(tǒng)狀態(tài)向量和總的誤差偏差矩陣。
為了驗(yàn)證本文ACKF 算法用于動(dòng)態(tài)諧波檢測(cè)的有效性,選擇KF,CKF 和ACKF 3 種算法,通過仿真分析比較這3 種算法進(jìn)行動(dòng)態(tài)諧波檢測(cè)的性能。
為了評(píng)估算法性能,引入評(píng)價(jià)指標(biāo)幅值均方根誤差Ra和相角均方根誤差Rp分別衡量檢測(cè)諧波幅值和相角的準(zhǔn)確度,其計(jì)算式分別為:
首先將ACKF 算法與CKF 算法進(jìn)行比較驗(yàn)證。仿真分析采用如表1 所示的動(dòng)態(tài)諧波信號(hào)頻譜,其他仿真參數(shù)設(shè)置為采樣頻率4 000 Hz,采樣時(shí)間0.2 s,數(shù)據(jù)共800 個(gè)點(diǎn),其中在0.1 s 時(shí)幅值增加1.5 倍,同時(shí)加入高斯噪聲,信噪比為50 dB。
表1 動(dòng)態(tài)諧波信號(hào)頻譜表Table 1 Dynamic harmonic signal spectrum
CKF 算法用于動(dòng)態(tài)諧波檢測(cè)的仿真結(jié)果如圖1~圖6 所示。
圖1 原始諧波信號(hào)與CKF算法跟蹤信號(hào)Fig.1 Original harmonic signal and tracking signal with CKF algorithm
圖1 為對(duì)原始動(dòng)態(tài)諧波信號(hào)的動(dòng)態(tài)跟蹤效果,圖2—圖6 展示了該算法實(shí)時(shí)估計(jì)出的各次諧波的幅值和相位。由圖1 可以看出CKF 算法可以追蹤到動(dòng)態(tài)諧波信號(hào),但在初始時(shí)刻出現(xiàn)了一些偏差。由圖2—圖6 中的波形可以得出CKF 算法對(duì)基波和各次諧波的幅值估計(jì)在初始位置存在誤差,需要進(jìn)一步修正。
圖2 CKF算法估計(jì)的基波幅值與相角Fig.2 Fundamental amplitude and phase angle estimated by CKF algorithm
圖3 CKF算法估計(jì)的5次諧波幅值與相角Fig.3 Fifth harmonic amplitude and phase angle estimated by CKF algorithm
圖4 CKF算法估計(jì)的7次諧波幅值與相角Fig.4 Seventh harmonic amplitude and phase angle estimated by CKF algorithm
圖5 CKF算法估計(jì)的11次諧波幅值與相角Fig.5 11th harmonic amplitude and phase angle estimated by CKF algorithm
圖6 CKF算法估計(jì)的13次諧波幅值與相角Fig.6 13th harmonic amplitude and phase angle estimated by CKF algorithm
ACKF 算法用于動(dòng)態(tài)諧波檢測(cè)的仿真結(jié)果如圖7—圖12 所示。其中圖7 為對(duì)原始動(dòng)態(tài)諧波信號(hào)的動(dòng)態(tài)跟蹤效果,圖8—圖12 展示了算法實(shí)時(shí)估計(jì)出的各次諧波的幅值和相位。
圖7 原始諧波信號(hào)與ACKF算法跟蹤信號(hào)Fig.7 Original harmonic signal and tracking signal with ACKF algorithm
圖8 ACKF算法估計(jì)的基波幅值與相角Fig.8 Fundamental amplitude and phase angle estimated by ACKF algorithm
圖9 ACKF算法估計(jì)的5次諧波幅值與相角Fig.9 Fifth harmonic amplitude and phase angle estimated by ACKF algorithm
圖10 ACKF算法估計(jì)的7次諧波幅值與相角Fig.10 Seventh harmonic amplitude and phase angle estimated by ACKF algorithm
圖11 ACKF算法估計(jì)的11次諧波幅值與相角Fig.11 11th harmonic amplitude and phase angle estimated by ACKF algorithm
圖12 ACKF算法估計(jì)的13次諧波幅值與相角Fig.12 13th harmonic amplitude and phase angle estimated by ACKF algorithm
由圖7 中的跟蹤波形可以看出,ACKF 算法的準(zhǔn)確度更高,在初始時(shí)刻的波形抖動(dòng)很小。由圖8—圖12 可以看出ACKF 算法在幅值和相位上能準(zhǔn)確估計(jì)出真實(shí)值。
最后將KF 算法同樣應(yīng)用于對(duì)基波、5 次諧波、7次諧波、11 次諧波和13 次諧波的動(dòng)態(tài)諧波信號(hào)檢測(cè)中,比較結(jié)果見表2。從表2 中可以看出,在有噪聲干擾下,ACKF 算法的動(dòng)態(tài)諧波檢測(cè)精度在總體上要高于CKF 算法與KF 算法,能夠很好地跟蹤動(dòng)態(tài)諧波信號(hào),同時(shí)也能很精確地估計(jì)出各階次諧波的幅值和相位。
表2 3種算法的Ra,Rp值比較Table 2 Comparison of Ra and Rp among three algorithms
有源濾波器(Active Power Filter,APF)[24-25]主要用于抑制電力系統(tǒng)的動(dòng)態(tài)諧波,其高效工作的基礎(chǔ)是實(shí)時(shí)、準(zhǔn)確地檢測(cè)諧波,傳統(tǒng)的諧波檢測(cè)方法如傅立葉變換(Fast Fourier Transform,F(xiàn)FT),雖然有廣泛應(yīng)用,但一般會(huì)有1 個(gè)周期以上的時(shí)間延遲,同時(shí)存在檢測(cè)精度不高等問題。這里嘗試將ACKF 算法應(yīng)用到APF 中,其系統(tǒng)原理框圖如圖13 所示。
圖13 有源濾波器系統(tǒng)結(jié)構(gòu)框圖Fig.13 Structure diagram of active filter system
基于Simulink 搭建400 V 有源濾波器仿真系統(tǒng)。采樣頻率設(shè)置為4 000 Hz,仿真時(shí)間為0.2 s,采集電流數(shù)據(jù)共800 個(gè)點(diǎn),其中在0.1 s 時(shí)將負(fù)荷增加1.2 倍,同時(shí)加入高斯噪聲,信噪比為50 dB,其仿真結(jié)果如圖14—圖17 所示。
圖14 ACKF算法估計(jì)的5次諧波幅值與相角Fig.14 Fifth harmonic amplitude and phase angle estimated by ACKF algorithm in active filter simulation
圖15 ACKF算法估計(jì)的7 次諧波幅值與相角Fig.15 Seventh harmonic amplitude and phase angle estimated by ACKF algorithm in active filter simulation
圖16 ACKF算法估計(jì)的11次諧波幅值與相角Fig.16 11th harmonic amplitude and phase angle estimated by ACKF algorithm in active filter simulation
圖17 ACKF算法估計(jì)的13次諧波幅值與相角Fig.17 13th harmonic amplitude and phase angle estimated by ACKF algorithm in active filter simulation
圖14—圖17 展示了ACKF 算法實(shí)時(shí)估計(jì)出的5 次,7 次,11 次,13 次諧波的幅值和相位,由仿真結(jié)果可以看出,在有噪聲干擾下,ACKF 算法應(yīng)用于有源濾波器諧波電流檢測(cè),實(shí)時(shí)性更好、精度更高。
以電網(wǎng)某些運(yùn)行工況下出現(xiàn)動(dòng)態(tài)諧波信號(hào)為切入點(diǎn),介紹了卡爾曼濾波動(dòng)態(tài)諧波檢測(cè)算法的基本原理,針對(duì)傳統(tǒng)卡爾曼濾波只能處理線性關(guān)系下的動(dòng)態(tài)諧波信號(hào)等問題,采用了一種容積卡爾曼濾波算法,在描述了容積卡爾曼濾波的原理后詳細(xì)推導(dǎo)了其實(shí)現(xiàn)的3 個(gè)步驟,然后又針對(duì)容積卡爾曼濾波算法存在的缺點(diǎn)引入了基于漸消記憶指數(shù)加權(quán)法的噪聲估值算法,將兩者結(jié)合起來構(gòu)成了自適應(yīng)容積卡爾曼濾波(ACKF)算法,最后通過仿真驗(yàn)證了本文提出的ACKF 算法在動(dòng)態(tài)諧波檢測(cè)中的有效性。