，，，，,2

(1.三峽大學(xué) 電氣與新能源學(xué)院，湖北 宜昌 443002；2.國(guó)家電網(wǎng)直流建設(shè)分公司 ，北京 100052)

1 引言

卡爾曼濾波算法是R.E.Kalman于1960年提出的適合數(shù)字計(jì)算機(jī)的遞推濾波方法，它采用狀態(tài)方程和量測(cè)方程組成的狀態(tài)空間來(lái)描述濾波器，該方法所研究的信號(hào)過(guò)程除了可以是平穩(wěn)的純量隨機(jī)過(guò)程，也可包括非平穩(wěn)的向量隨機(jī)過(guò)程。本文將所有引起負(fù)荷變化的因素(如氣候、降雨、溫度)歸為隨機(jī)系統(tǒng)的隨機(jī)因素來(lái)處理；為了對(duì)噪聲協(xié)方差有比較準(zhǔn)確的估計(jì)，運(yùn)用時(shí)變次優(yōu)無(wú)偏噪聲估值器對(duì)噪聲協(xié)方差進(jìn)行自適應(yīng)估計(jì)，提高了濾波預(yù)測(cè)精度。用實(shí)際數(shù)據(jù)進(jìn)行Matlab仿真驗(yàn)證的結(jié)果表明，該方法使誤差得到了有效控制。

卡爾曼濾波是基于最小均方誤差地遞推預(yù)測(cè)算法，通過(guò)對(duì)模型參數(shù)的辨識(shí)，實(shí)現(xiàn)對(duì)觀測(cè)序列的預(yù)測(cè)。其用于負(fù)荷預(yù)測(cè)算法的指導(dǎo)思想是：將負(fù)荷分為確定分量和隨機(jī)分量，確定分量一般采用線性回歸模型預(yù)測(cè)，隨機(jī)分量采用卡爾曼濾波算法預(yù)測(cè)。文獻(xiàn)[2]較早地將隨機(jī)系統(tǒng)狀態(tài)模型辨識(shí)技術(shù)用于電力系統(tǒng)負(fù)荷預(yù)報(bào)；文獻(xiàn)[3]提出了用卡爾曼濾波結(jié)合最小二乘法線性擬合及3次樣條插值的數(shù)學(xué)方法；文獻(xiàn)[4]提出了自校正移動(dòng)窗卡爾曼濾波-最小二乘算法，該算法通過(guò)移動(dòng)窗最小二乘算法得到較精確的參數(shù)估計(jì)；文獻(xiàn)[5]則根據(jù)對(duì)測(cè)量新息做出估計(jì)，實(shí)現(xiàn)了一步預(yù)測(cè)基礎(chǔ)上的二次修正；文獻(xiàn)[6]分別以負(fù)荷需求序列和時(shí)變參數(shù)序列為狀態(tài)變量建立了兩個(gè)狀態(tài)空間模型，采用兩段自適應(yīng)卡爾曼濾波方法進(jìn)行負(fù)荷預(yù)測(cè)；文獻(xiàn)[7]提出了在預(yù)測(cè)值的基礎(chǔ)上加上溫度修正值的負(fù)荷預(yù)測(cè)方法。

2 卡爾曼濾波理論簡(jiǎn)介

假設(shè)有如下的自由狀態(tài)定常線性系統(tǒng)模型。

(1)

z(k)=Hx(k)+v(k)

(2)

當(dāng)系統(tǒng)的初始狀態(tài)x(k0)及系統(tǒng)干擾ω(k)和測(cè)量噪聲v(k)是隨機(jī)過(guò)程，為了使問(wèn)題簡(jiǎn)化，需要對(duì)系統(tǒng)干擾、測(cè)量噪聲及系統(tǒng)的初始狀態(tài)作一定的假設(shè)。系統(tǒng)干擾和測(cè)量噪聲是零均值或非零均值的白噪聲過(guò)程或高斯白噪聲過(guò)程，即

(3)

式中，Q(k)非負(fù)定矩陣，是ω(k)的強(qiáng)度函數(shù)；R(k)為對(duì)稱正定矩陣，是v(k)的方差強(qiáng)度矩陣；d(k)是狄拉克函數(shù)；A是n×n維狀態(tài)轉(zhuǎn)移矩陣；H是1×n維測(cè)量矩陣。系統(tǒng)干擾和測(cè)量噪聲互不相關(guān)，且與初始狀態(tài)x(k0)無(wú)關(guān)

Cov[ω(k),v(k)]=0

Cov[x(k0),ω(k)]=0

Cov[x(k0),v(k)]=0

(4)

利用正交投影定理推導(dǎo)卡爾曼濾波遞推公式如下：

狀態(tài)預(yù)測(cè)方程：

(5)

誤差協(xié)方差預(yù)

Pk|k-1=APk-1AT+Qk-1

(6)

狀態(tài)估計(jì)校正：

(7)

誤差協(xié)方差估計(jì)校正：

Pk|k=(I-KkHk)Pk|k-1

(8)

卡爾曼增益：

(9)

詳細(xì)推導(dǎo)過(guò)程見(jiàn)文獻(xiàn)[8]。

3 次優(yōu)無(wú)偏估值器的應(yīng)用改進(jìn)

為盡量減少隨機(jī)噪聲因素對(duì)預(yù)測(cè)結(jié)果的影響，在利用觀測(cè)數(shù)據(jù)進(jìn)行遞推濾波的同時(shí)，需要不斷了不斷地修正模型中噪聲協(xié)方差陣，以提高濾波精度。

確定系統(tǒng)干擾ω(k)、v(k)統(tǒng)計(jì)特性，也就是噪聲強(qiáng)度矩陣Q(k)、R(k)的確定。Q(k)、R(k)選取好壞對(duì)濾波精度有直接的影響，動(dòng)態(tài)模型越不精確，這種影響就越大。但在實(shí)際應(yīng)用中，通過(guò)已知信息很難得到比較準(zhǔn)確的Q、R值，況且隨著新息的不斷加入，Q、R值的不準(zhǔn)確性越來(lái)越明顯，這樣就直接導(dǎo)致了預(yù)測(cè)結(jié)果偏差比較大。

因此，在對(duì)遞推方程(6)、(9)進(jìn)行計(jì)算時(shí)，有必要對(duì)噪聲協(xié)方差采取自適應(yīng)的方法，引入次優(yōu)無(wú)偏噪聲估值器[9]。

(10)

(11)

以上運(yùn)用次優(yōu)無(wú)偏估值器計(jì)算簡(jiǎn)單，可同時(shí)在線估計(jì)出系統(tǒng)噪聲協(xié)方差和測(cè)量噪聲協(xié)方差，但是在ε(k+1)出現(xiàn)較大值或系統(tǒng)階數(shù)較高時(shí)，濾波器易發(fā)散。

(12)

(13)

4 算例分析

本文以某實(shí)際電網(wǎng)2011年夏季8月7日至15日每天24小時(shí)的負(fù)荷作為歷史負(fù)荷。歷史負(fù)荷曲線見(jiàn)圖1，抽離平均值后的隨機(jī)負(fù)荷頻數(shù)直方圖見(jiàn)圖2。從圖1可見(jiàn)日負(fù)荷曲線雖然呈較強(qiáng)的非線性，但是也有明顯的以天為周期的周期性，且峰谷負(fù)荷出現(xiàn)在固定的時(shí)間段。從圖2可見(jiàn)隨機(jī)負(fù)荷服從正太分布。

在進(jìn)行數(shù)據(jù)分析前先要對(duì)異常數(shù)據(jù)進(jìn)行處理，再用上文所述方法進(jìn)行仿真預(yù)測(cè)。

圖1 歷史負(fù)荷曲線圖

圖2 隨機(jī)負(fù)荷頻數(shù)直方圖

4.1 模型參數(shù)辨識(shí)

4.1.1 模型階數(shù)辨識(shí)

設(shè)已知被辨識(shí)系統(tǒng)的隨機(jī)負(fù)荷序列l(wèi)0,l1,…,lN。并按下式構(gòu)造Hankel矩陣(簡(jiǎn)稱H矩陣)。

(14)

其維數(shù)為m×m。在弱噪聲情況下，可分別令m=1,2,3,…，求出每個(gè)m值下H陣行列式值，當(dāng)行列式值達(dá)到極大時(shí)的m值，即為系統(tǒng)的階數(shù)n；在強(qiáng)噪聲情況下，不能直接用隨機(jī)負(fù)荷序列，而是采用隨機(jī)負(fù)荷序列的自相關(guān)系數(shù)構(gòu)造H矩陣，即用隨機(jī)負(fù)荷序列的自相關(guān)系數(shù)ρi代替li，然后計(jì)算H陣的行陣式的值，當(dāng)行陣式值取極小值時(shí)的m即為系統(tǒng)階數(shù)[9]。

4.1.2 模型轉(zhuǎn)移矩陣和測(cè)量矩陣的辨識(shí)

本文采用相關(guān)函數(shù)——最小二乘相結(jié)合的方法辨識(shí)系統(tǒng)模型的轉(zhuǎn)移矩陣?；舅枷胧前驯孀R(shí)分兩步進(jìn)行。第一步，利用相關(guān)函數(shù)法對(duì)隨機(jī)負(fù)荷分量進(jìn)行一次相關(guān)分析，獲得被辨識(shí)對(duì)象的相關(guān)函數(shù)；第二步利用最小二乘法進(jìn)一步估計(jì)模型的參數(shù)。

狀態(tài)轉(zhuǎn)移矩陣A、測(cè)量矩陣H采用與其等價(jià)的，辨識(shí)參數(shù)最少的典范性[2]。

H=[1 0 … 0]1×n

輸出隨機(jī)負(fù)荷序列的相關(guān)函數(shù)為：

R(n)=E{Z(k+n)Z(k)}

(15)

根據(jù)平穩(wěn)過(guò)程遍歷性，統(tǒng)計(jì)均值可以用時(shí)間平均代替，(15)式可改寫(xiě)為：

(16)

可以推得(16)式的遞推式為：

(17)

式子中ak為A中第一行參數(shù)，利用(16)、(17)式，采用最小二乘法即可計(jì)算出模型參數(shù)ak的值。

4.2 算法步驟

(1)分析歷史負(fù)荷數(shù)據(jù)，處理異常數(shù)據(jù)點(diǎn)。

(2)根據(jù)4.1求模型參數(shù)，模型參數(shù)及遞推初始值見(jiàn)表1。

表1 模型參數(shù)

4.3 仿真結(jié)果

下面給出利用上述方法得到的預(yù)測(cè)誤差和預(yù)測(cè)結(jié)果，預(yù)測(cè)相對(duì)誤差見(jiàn)圖3，預(yù)測(cè)結(jié)果對(duì)比見(jiàn)圖4。為了衡量負(fù)荷預(yù)測(cè)的有效性用百分相對(duì)誤差(Relative Percentage Error ，RPE)來(lái)計(jì)算某一時(shí)刻的誤差，用平均絕對(duì)百分誤差(Mean Absolute Percentage Error,MAPE)來(lái)反映誤差總體大小。其計(jì)算式如下：

(18)

(19)

圖3 相對(duì)誤差圖

圖4 負(fù)荷對(duì)比圖

其平均相對(duì)誤差為-0.38%，平均絕對(duì)誤差為1.24%。由此可見(jiàn)對(duì)次優(yōu)無(wú)偏噪聲估值器的應(yīng)用改進(jìn)是有效的。

5 結(jié)論

本文運(yùn)用帶次優(yōu)無(wú)偏噪聲估值器的卡爾曼濾波理論建立了短期負(fù)荷預(yù)測(cè)模型，并進(jìn)行短期負(fù)荷預(yù)測(cè)。通過(guò)Matlab仿真計(jì)算證明了該預(yù)測(cè)模型的可行性。同時(shí)針對(duì)負(fù)荷預(yù)測(cè)的特點(diǎn)，通過(guò)對(duì)次優(yōu)無(wú)偏噪聲估值器的改進(jìn)，提高了預(yù)測(cè)精度。

參數(shù)辨識(shí)后的帶次優(yōu)無(wú)偏估值器的卡爾曼濾波器始終處于穩(wěn)定狀態(tài)，說(shuō)明卡爾曼濾波模型參數(shù)辨識(shí)方法是可行的。

另外，將其他影響負(fù)荷的因素(如溫度，降雨，體感溫度)等與本文的方法綜合考慮，可期望得到更好的預(yù)測(cè)結(jié)果。

[1] P.C.Gupta,K.Yamada.Adaptive short-term forecasting of hourly loads usingweather information[J],IEEE Transactions powerApparatusSystem,PAS-91,1972:2085-2094.

[2] 李洪心，易允文.隨機(jī)辨識(shí)系統(tǒng)在電力系統(tǒng)負(fù)荷預(yù)報(bào)中的應(yīng)用[J].信息與控制，1985(3):31-34.

[3] 張民，鮑海，晏玲，等.基于卡爾曼濾波的短期負(fù)荷預(yù)測(cè)方法的研究[J].電網(wǎng)技術(shù)，2003,27(10)：39-42.

[4] 王正歐，張建平，等.電力系統(tǒng)超短期負(fù)荷預(yù)報(bào)的一種方法[J].天津大學(xué)學(xué)報(bào)，1992(2)：125-130.

[5] 孟思齊，楊洪耕.基于卡爾曼濾波二次修正的短期負(fù)荷預(yù)測(cè)[J]，電網(wǎng)與水力發(fā)電進(jìn)展,2008,24(2)：78-82.

[6] 馬靜波，楊洪耕.自適應(yīng)卡爾曼濾波在電力系統(tǒng)短期負(fù)荷預(yù)測(cè)中的應(yīng)用[J].電網(wǎng)技術(shù)，2005,29(1)：75-79.

[7] 李明干，孫健利，劉沛.基于卡爾曼濾波的電力系統(tǒng)短期負(fù)荷預(yù)測(cè)[J].繼電器，2004,32(4)：9-12.

[8] 王志賢.最優(yōu)狀態(tài)估計(jì)與系統(tǒng)辨識(shí)[M].西安：西北工業(yè)大學(xué)出版社，2004.

[9] 劉勝，張紅梅.最優(yōu)估計(jì)理論[M].北京：科學(xué)出版社，2011.