無跡卡爾曼濾波算法在污水處理BSM1過程中的應(yīng)用研究

曾 靜 桑耀凱 李 元

(沈陽化工大學(xué)信息工程學(xué)院 遼寧 沈陽 110142)

0 引 言

污水處理是水循環(huán)的重要環(huán)節(jié)，涉及到復(fù)雜的生物和化學(xué)反應(yīng)。污水處理過程(Wastewater Treatment Process,WWTP)是一個典型的非線性系統(tǒng)，其出水水質(zhì)與環(huán)境的可持續(xù)性密切相關(guān)，通常由環(huán)境立法進(jìn)行管理，在滿足環(huán)境法規(guī)的嚴(yán)格要求的同時，還要求確保工藝安全并將運行成本降至最低。污水入水流量和成分的顯著變化導(dǎo)致設(shè)計污水處理廠的先進(jìn)控制和監(jiān)測方案的復(fù)雜性增加。污水處理過程的典型仿真模型BSM1(Benchmark Simulation Model no.1)由活性污泥1號模型(ASMI)和雙指數(shù)沉淀速度模型組成，作為標(biāo)準(zhǔn)仿真平臺廣泛應(yīng)用于污水過程的控制、優(yōu)化、系統(tǒng)監(jiān)視及故障診斷等研究中[1-2]。

關(guān)于污水處理過程的控制算法有許多研究成果。早期利用比例積分(Proportional-Integral，PI)控制進(jìn)行污水處理廠營養(yǎng)物去除，PI控制可以達(dá)到運行中的穩(wěn)定性要求，但不能有效處理約束[3]。近年來模型預(yù)測控制(Model Predictive Control，MPC)等先進(jìn)控制方法也引起了人們的廣泛關(guān)注。英國蘭卡斯特的污水處理廠的應(yīng)用證明了模型預(yù)測控制的有效性[4]。Francisco等[5]為BSM1基準(zhǔn)仿真模型1號開發(fā)了MPC控制器，文獻(xiàn)[6-8]考慮了污水處理廠MPC的同步設(shè)計和控制問題。Zeng等[9]將經(jīng)濟(jì)的MPC算法(Economic Model Predictive Control，EMPC)應(yīng)用于BSM1，設(shè)計考慮終端代價(Terminal Cost)的經(jīng)濟(jì)模型預(yù)測控制器，并將EMPC性能與常規(guī)MPC、PI進(jìn)行了分析比較，在保證系統(tǒng)穩(wěn)定性的基礎(chǔ)上顯著提高了出水質(zhì)量，同時降低了全局能耗指標(biāo)。

雖然關(guān)于污水處理廠控制系統(tǒng)設(shè)計的研究成果很多，但對污水處理廠狀態(tài)估計的關(guān)注相對較少。狀態(tài)估計是基于輸出測量和系統(tǒng)模型重構(gòu)系統(tǒng)狀態(tài)的過程。狀態(tài)估計與污水處理廠的控制和監(jiān)測密切相關(guān)，是污水處理廠運行中的一個重要課題。污水處理中的許多相關(guān)狀態(tài)是不可測量的或受到嚴(yán)重的噪聲污染。擴(kuò)展卡爾曼濾波(Extended Kalman Filter，EKF)曾被應(yīng)用于BSM1的狀態(tài)估計[10]。然而EKF濾波器線性化存在局限性,例如對于較強非線性模型需要計算的步長盡可能小。Julier等[11-13]提出了一種新的非線性濾波方案來替代EKF在目標(biāo)跟蹤領(lǐng)域的應(yīng)用，后來被命名為無跡卡爾曼濾波(Unscented Kalman Filter，UKF)。

基于上述考慮，本文針對污水處理的典型模型BSM1設(shè)計了無跡卡爾曼濾波器，并與擴(kuò)展卡爾曼濾波的狀態(tài)估計結(jié)果進(jìn)行比較。主要工作包括：(1) 詳細(xì)地描述EKF和UKF的算法理論及二者的優(yōu)缺點，并描述無跡卡爾曼濾波的理論依據(jù)。(2) 將EKF與UKF在同一噪聲下應(yīng)用于BSM1仿真，將仿真結(jié)果進(jìn)行對比分析，得出二者方案各自的優(yōu)缺點和最好的估計方案。(3) 在不同天氣條件不同噪聲大小狀況下，對EKF、UKF兩種方案在BSM1應(yīng)用進(jìn)行比較,分析比較結(jié)果。

1 EKF算法描述

擴(kuò)展卡爾曼濾波算法理論在于對原系統(tǒng)方程和測量方程采取一階泰勒展開來近似逼近實際系統(tǒng)非線性模型。通過求偏導(dǎo)計算所需的雅可比矩陣確定出當(dāng)前的系統(tǒng)矩陣F(k)和輸出矩陣H(k),為了不失一般性，我們設(shè)系統(tǒng)模型的狀態(tài)方程和輸出方程如下：

x(k)=f(x(k-1))+w(k-1)

(1)

z(k)=h(x(k))+v(k)

(2)

式中：f(·)為系統(tǒng)狀態(tài)方程、h(·)為系統(tǒng)輸出方程，二者為非線性函數(shù)；x為系統(tǒng)的狀態(tài)變量；w(k-1)為系統(tǒng)k-1時刻過程噪聲；v(k)為系統(tǒng)k時刻測量噪聲；z(k)為系統(tǒng)在k時刻測量值。其中噪聲均服從零均值的高斯分布。為了確定系統(tǒng)矩陣和輸出矩陣，對狀態(tài)方程和測量方程求偏導(dǎo)，得：

(3)

(4)

(5)

P(k|k-1)=F(k-1)P(k-1|k-1)FT(k-1)+Q

(6)

式中：P(k-1|k-1)為k-1時刻的誤差的協(xié)方差矩陣；P(k|k-1)為k-1時刻預(yù)測k時刻的誤差協(xié)方差矩陣；Q為過程噪聲的協(xié)方差矩陣。測量更新如下：

K(k)=P(k|k-1)HT(k)[HT(k)P(k|k-1)H(k)+R]-1

(7)

(8)

P(k|k)=P(k|k-1)-K(k)H(k)P(k|k-1)

(9)

2 UKF濾波算法描述

(1) 利用初始估計值，生成初始的L個Sigma點：

(10)

(2) 利用狀態(tài)方程計算Sigma點的一步預(yù)測：

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

(11)

(3) 對經(jīng)過非線性傳遞之后Sigma點加權(quán)計算：

(12)

(4) 計算協(xié)方差矩陣的一步預(yù)測：

(13)

(5) 利用測量方程計算一步Sigma點的測量值：

Z(i)(k|k-1)=h(χ(i)(k|k-1))

(14)

(6) 計算測量值的預(yù)測加權(quán)求和：

(15)

(7) 計算新息的協(xié)方差：

(16)

(8) 計算協(xié)方差：

(17)

(9) 計算Kalman增益：

(18)

(10) 更新誤差協(xié)方差矩陣：

PX(k|k)=PX(k|k-1)-K(k)PzzK(k)T

(19)

(11) 利用新息更新狀態(tài)，z(k)為當(dāng)前實測輸出：

(20)

3 仿真實驗

將上述介紹的UKF算法應(yīng)用于BSM1模型，并與EKF在不同天氣和不同大小噪聲條件下的估計性能進(jìn)行比較分析。BSM1的模型描述詳見文獻(xiàn)[10，14]。三種不同天氣條件數(shù)據(jù)文件分別為干燥天氣數(shù)據(jù)文件、降雨天的數(shù)據(jù)文件和暴風(fēng)雨天的數(shù)據(jù)文件，詳見國際水協(xié)會的官網(wǎng)(http://www.benchmarkwwtp.org)。干燥天氣包含了兩周相同的動態(tài)干燥條件數(shù)據(jù)；雨天數(shù)據(jù)文件中包含一周的動態(tài)干燥條件數(shù)據(jù)和發(fā)生在第二周的一次長降雨事件；暴風(fēng)雨天的數(shù)據(jù)文件中包含了一周的動態(tài)干燥條件和發(fā)生于第二周期間的兩次暴風(fēng)雨事件。生物反應(yīng)器和理想分離器的初始值通過模擬工廠100天的輸入，并在干燥條件下模擬14天計算得到的。由于三種不同天氣狀況的輸入數(shù)據(jù)前一周是完全相同的，因此仿真直接從第二周開始。在仿真中，假設(shè)對整個過程測量進(jìn)行采樣是同步和周期性的，采樣周期取15 min，系統(tǒng)的狀態(tài)變量包括5個反應(yīng)器及二次沉降池變量(SALK、XBH、XI、XS、XND、XP、XBA、SI、SS、SO、SNO、SNH、SND)見表1，反應(yīng)器及沉降池可測量變量(SO，SALK，SNH，SNO，COD，CODf，BOD，TSS)見表2。在每組的仿真中，將相同的測量噪聲、過程噪聲序列和相同的初始值分別應(yīng)用到EKF、UKF方案上，過程噪聲和測量噪聲均為有界高斯白噪聲。

表1 反應(yīng)器中狀態(tài)變量

表2 測量輸出量

首先比較了UKF、EFK方案在干燥天氣下的性能，在本組仿真中78個狀態(tài)變量分別產(chǎn)生零均值，標(biāo)準(zhǔn)差為0.2x0的高斯白噪聲，并有界于[-0.2x0,0.2x0],測量傳感器自帶噪聲為零均值，標(biāo)準(zhǔn)差為0.2y0的高斯白噪聲，并有界于[-0.2y0,0.2y0]，y0=Cx0，其中：x0為初始狀態(tài)值；y0為初始的測量值；C為輸出測量矩陣。濾波方案所采用的初始估計值xinit均設(shè)置為1.02x0，根據(jù)協(xié)方差計算公式得出過程噪聲和測量噪聲的協(xié)方差矩陣分別為Q=diag((0.2x0)2)和R=diag((0.2y0)2)，初始的誤差方差矩陣設(shè)為P(0)=Q，設(shè)過程噪聲和測量噪聲以及狀態(tài)誤差之間均不相關(guān)。圖1為反應(yīng)器3中實際過程狀態(tài)和UKF和EKF的濾波軌跡。其中：虛線為EKF軌跡；點線為UKF軌跡；實線為實際過程狀態(tài)軌跡。其他反應(yīng)器和沉淀池具有與其類似的特點，可以看出EKF、UKF都能很好實現(xiàn)跟蹤，大多狀態(tài)跟蹤效果都挺好，只有個別狀態(tài)估計偏差較大如XI、XP、XBH。為了更加細(xì)致地比較EKF、UKF的性能優(yōu)劣，我們需要定量計算一下二者在每個時間點的誤差均值并繪成圖像。比較方案將每組誤差統(tǒng)一標(biāo)準(zhǔn)化之后采用RMS(Root Mean Square)。將整個時間段[t0,ts]兩種方案每個狀態(tài)的最大估計誤差值作分母來標(biāo)準(zhǔn)化每個狀態(tài)誤差，t0、ts分別為初始時間和終止時間，狀態(tài)誤差如式(21)所示。

(21)

式中：ei(tk)為一個標(biāo)準(zhǔn)化的狀態(tài)變量的誤差序列，k為仿真時刻。計算78個狀態(tài)在時間段[t0,ts]內(nèi)每一個時刻的整體誤差如式(22)所示。

(22)

采用誤差均值及最大值來比較兩種方案估計性能的優(yōu)劣，計算公式如下：

(23)

Max|e|=Max(e(tk))

(24)

式中：K為[t0,ts]時間段計算的點數(shù)。根據(jù)式(21)-式(24)，通過定量計算EKF、UKF的誤差均值分別為2.94、2.39，最大標(biāo)準(zhǔn)化誤差分別為5.04、4.15，可以看出EKF給出的估計性能沒有UKF好。干燥天氣數(shù)據(jù)下UKF對比EKF性能大概提升了12.7%，圖2為二者的誤差軌跡，實線為EKF的誤差軌跡，虛線為UKF誤差軌跡，可以看到實線始終高于虛線，一開始二者比較相近，隨著輸入數(shù)據(jù)和仿真時長的增加，二者的曲線也在不斷波動變化。

圖1 干燥天氣下反應(yīng)器3中EKF軌跡、UKF軌跡、 實際過程狀態(tài)軌跡

圖2 EKF誤差和UKF誤差軌跡

在另一組的仿真中，在不同天氣狀況分別對兩種估計方案的性能進(jìn)行了測試。除了系統(tǒng)輸入數(shù)據(jù)外，所有的條件都與上一組仿真完全相同。過程噪聲取標(biāo)準(zhǔn)差為0.2x0的高斯白噪聲，測量噪聲取標(biāo)準(zhǔn)差為0.2y0的高斯白噪聲，圖3給出了雨天狀況下反應(yīng)器3中的13種成分的實際狀態(tài)和兩種估計方案的軌跡。其中：虛線為EKF軌跡；點線為UKF軌跡；實線為實際過程狀態(tài)軌跡。其他反應(yīng)器和沉淀池具有與反應(yīng)器3類似的特點。此天氣狀況和當(dāng)前噪聲環(huán)境下，EKF、UKF的估計誤差均值分別為2.6、2.26。UKF對比EKF估計精度提升了13.08%。暴風(fēng)雨天狀況下，EKF、UKF的估計誤差均值分別為2.9、2.43，UKF較EKF估計精度提升16.2%。這表明了盡管惡劣天氣下輸入數(shù)據(jù)波動比較大，UKF算法仍然比EKF濾波算法的估計性能好。

圖3 雨天狀況下反應(yīng)器3中EKF軌跡、UKF軌跡、 實際過程狀態(tài)軌跡

在最后一組仿真中，通過改變過程噪聲和測量噪聲的大小來測試EKF、UKF性能的優(yōu)劣，初始條件xinit仍設(shè)為1.02x0，通過參數(shù)wQ和wR來調(diào)整噪聲的幅度，假設(shè)過程噪聲大小為wQx0，并有界于[-wQx0,wQx0]，則Q=diag((wQx0)2)，測量噪聲的大小為wRy0,并有界于[-wRy0,wRy0]，均為有界的高斯白噪聲，有y0=Cx0，R=diag((wRy0)2)。初始誤差方差陣P(0)=Q，不同大小過程噪聲和測量噪聲下EKF和UKF狀態(tài)估計的誤差均值如表3所示?？梢钥闯鲈诟稍锾臁⒂晏旌捅┯晏烨闆r下，無論過程噪聲及測量噪聲大小如何改變，UKF估計性能始終優(yōu)于EKF估計效果。

表3 不同噪聲、天氣條件下估計誤差均值比較

4 結(jié) 語

無跡卡爾曼濾波算法通過對非線性函數(shù)的概率密度分布近似來完成估計過程，與擴(kuò)展卡爾曼算法相比不需要對Jacobin矩陣進(jìn)行求導(dǎo)，也沒有忽略高階項，因此對非線性系統(tǒng)的狀態(tài)估計有較高的計算精度。本文將無跡卡爾曼濾波算法引入到污水處理過程帶的狀態(tài)估計問題中，并與擴(kuò)展卡爾曼濾波算法相比較，結(jié)果表明無論天氣狀況如何，在相同模型噪聲與測量噪聲下的UKF估計精度始終高于EKF的估計精度，并且魯棒性好，UKF較EKF的估計精度提升了10%～20%左右，顯著提高了對復(fù)雜污水處理過程的狀態(tài)估計效果。