吳宏亮，趙忠蓋，劉 飛

(1.輕工過(guò)程先進(jìn)控制教育部重點(diǎn)實(shí)驗(yàn)室,江蘇 無(wú)錫 214122;2.江南大學(xué) 物自動(dòng)化研究所，江蘇 無(wú)錫 214122)

0 引言

間歇過(guò)程是一類重復(fù)生產(chǎn)過(guò)程，按照相同的工序以批次為單位進(jìn)行產(chǎn)品生產(chǎn)，越來(lái)越多地應(yīng)用于食品、生物制藥以及化工領(lǐng)域。間歇過(guò)程中很多關(guān)鍵參數(shù)直接與產(chǎn)品質(zhì)量相關(guān)，需要對(duì)其進(jìn)行實(shí)時(shí)監(jiān)測(cè)。但是，受限于傳感器技術(shù)的發(fā)展，間歇過(guò)程中的許多狀態(tài)變量難以測(cè)量或者測(cè)量成本很高，因此對(duì)間歇過(guò)程狀態(tài)變量的估計(jì)尤為重要，一直是工業(yè)界和學(xué)術(shù)界的關(guān)注焦點(diǎn)[1]?？柭鼮V波(Kalman filter,KF)是一種廣泛應(yīng)用的最優(yōu)狀態(tài)估計(jì)的方法，但它適用于線性系統(tǒng)，而間歇過(guò)程通常具有非線性特性。對(duì)于非線性系統(tǒng)，通常采用擴(kuò)展卡爾曼濾波(Extended Kalman filter,EKF)方法，無(wú)跡卡爾曼濾波和粒子濾波方法等[2-5]。但是，這些方法往往是連續(xù)過(guò)程狀態(tài)估計(jì)的簡(jiǎn)單復(fù)制，忽略了批次之間的相關(guān)性。

近年來(lái)，為了提高間歇過(guò)程狀態(tài)估計(jì)的性能，迭代學(xué)習(xí)策略被引入狀態(tài)估計(jì)，考慮了批次運(yùn)行之間的重復(fù)性和時(shí)間方向動(dòng)態(tài)特性。迭代學(xué)習(xí)的思想來(lái)自于重復(fù)過(guò)程的控制，Arimoto等[6]提出一種迭代學(xué)習(xí)控制(iterative learning control,ILC)方法，針對(duì)重復(fù)運(yùn)行的被控系統(tǒng)，通過(guò)前一次或前幾次的控制誤差信息修正當(dāng)前操作的控制輸入，最終實(shí)現(xiàn)在整個(gè)時(shí)間區(qū)間上系統(tǒng)的輸出完全跟蹤上期望軌跡。基于迭代學(xué)習(xí)的思想，文獻(xiàn)[7-9]設(shè)計(jì)了確定性重復(fù)系統(tǒng)的迭代學(xué)習(xí)觀測(cè)器。然而，實(shí)際間歇生產(chǎn)存在噪聲干擾。

Zhao等[10]考慮包含隨機(jī)噪聲的非線性間歇過(guò)程，提出了一種用于批處理的迭代學(xué)習(xí)狀態(tài)估計(jì)算法，在當(dāng)前批次粒子濾波的基礎(chǔ)上，引入上一批次輸出值的估計(jì)誤差作為修正項(xiàng)，包含了上一批次的信息，但修正項(xiàng)難以確定。Cao[11]等人提出了一種魯棒迭代學(xué)習(xí)卡爾曼濾波方法，用于狀態(tài)和輸出矩陣均具有范數(shù)有界不確定性的重復(fù)過(guò)程系統(tǒng)的狀態(tài)估計(jì)，但只給出了一種先驗(yàn)類型估計(jì)。在文獻(xiàn)[12]中，提出了線性系統(tǒng)的迭代學(xué)習(xí)卡爾曼濾波器(Iterative learning Kalman filter,ILKF)，并在重復(fù)注塑機(jī)過(guò)程中進(jìn)行了估計(jì)研究，驗(yàn)證了ILKF在間歇過(guò)程中優(yōu)越的估計(jì)性能。

但是，文獻(xiàn)[12]中的ILKF方法是一種線性狀態(tài)估計(jì)方法，而大多數(shù)間歇過(guò)程非線性嚴(yán)重。本文將基于標(biāo)稱模型的擬線性化方法引入到間歇過(guò)程狀態(tài)估計(jì)中來(lái)，提出一種迭代學(xué)習(xí)擬線性卡爾曼濾波(Iterative learning quasilinear Kalman filter,ILQKF)方法。其中，標(biāo)稱模型是間歇過(guò)程含噪聲模型的期望[13]，可以根據(jù)機(jī)理建模求得?？紤]間歇過(guò)程正常運(yùn)行下，每個(gè)時(shí)刻的狀態(tài)在標(biāo)稱軌跡周邊變化，而狀態(tài)的變化呈現(xiàn)線性特征，因此ILQKF以實(shí)際狀態(tài)與標(biāo)稱狀態(tài)之間的誤差為新狀態(tài)，首先基于標(biāo)稱軌跡的擬線性化方法建立了與誤差相關(guān)的線性化模型；然后，應(yīng)用ILKF方法對(duì)誤差模型進(jìn)行狀態(tài)估計(jì)，而狀態(tài)軌跡等于誤差軌跡與標(biāo)稱軌跡之和。ILQKF在狀態(tài)估計(jì)中充分利用了時(shí)間和批次雙維的特性，并消除了重復(fù)誤差。此外，相比ILKF方法，ILQKF還具有更寬松的初值條件。ILKF方法需要假設(shè)系統(tǒng)的初值滿足正態(tài)分布且均值為零，這在間歇過(guò)程中是一個(gè)不合理的假設(shè)，比如發(fā)酵過(guò)程的初始狀態(tài)可能是一個(gè)已知的非零量，如發(fā)酵過(guò)程中底物的初始濃度和發(fā)酵罐中的初值生物質(zhì)濃度通常是一個(gè)給定的值，ILQKF方法更具有實(shí)用性。本文通過(guò)啤酒發(fā)酵例子驗(yàn)證了ILQKF算法優(yōu)越的性能。

1 迭代學(xué)習(xí)卡爾曼濾波

ILKF方法結(jié)合了KF和ILC的思想，充分利用了間歇過(guò)程的多批次重復(fù)性。

對(duì)于線性間歇過(guò)程系統(tǒng)模型(1)：

(1)

其中:k是系統(tǒng)的批次數(shù)，k∈[1,+∞)，t是采樣時(shí)間，t∈[0,M]，xk(t)∈Rn是第k批、第t時(shí)刻的狀態(tài)變量，yk(t)∈Rl是第k批、第t時(shí)刻的輸出變量，l≤n，d(t)是系統(tǒng)每批的重復(fù)干擾，為了設(shè)計(jì)ILKF，作以下假設(shè)[12]：

假設(shè)1：ωk(t)和vk(t)分別是獨(dú)立的雙維高斯過(guò)程噪聲和測(cè)量噪聲：

E[ωk(t)]=E[vk(t)]=0

E[ωi(j)vkT(t)] = 0

E[wi(j)wkT(t)] =Qδj,t;i,k

E[vi(j)vkT(t)] =Rδj,t;i,k

δj,t;i,k是雙維克羅內(nèi)克函數(shù)，僅當(dāng)j=t,i=k時(shí)，δj,t;i,k=I，否則δj,t;i,k=0。

假設(shè)2：d(t)為未知的有界重復(fù)干擾，隨時(shí)間變化。

假設(shè)3：初始條件xk(0)服從的獨(dú)立正態(tài)分布N(0,ζ)。

以上3個(gè)假設(shè)中,假設(shè)1與假設(shè)2相當(dāng)于間歇生產(chǎn)過(guò)程中的兩種干擾，而假設(shè)3對(duì)初值的假設(shè)不符合絕大多數(shù)情況，本文第2節(jié)和第3節(jié)將給出合理解決假設(shè)3的方法。

定義δ為批次方向的后向差分算子：

δxk(t)=xk(t)-xk-1(t)

δyk(t)=yk(t)-yk-1(t)

(2)

同樣可以定義δxk(t),δyk(t),δωk(t),δvk(t)。

可以將系統(tǒng)改寫為兩個(gè)子系統(tǒng)，分別為時(shí)間方向上的子系統(tǒng)∑T和批次方向上的子系統(tǒng)∑B：

(3)

(4)

可以發(fā)現(xiàn)在∑B中，xk(t)與xk-1(t)之間的轉(zhuǎn)移矩陣是單位陣，即批次子模型∑B是個(gè)臨界穩(wěn)定系統(tǒng)，可能會(huì)導(dǎo)致xk(t)的值不斷上升或震蕩，可以在原∑B子系統(tǒng)中加入系數(shù)τ<1，使系統(tǒng)穩(wěn)定[12]：

(5)

(6)

在間歇過(guò)程運(yùn)行時(shí)，系統(tǒng)干擾分為重復(fù)干擾d(t)，隨機(jī)噪聲ωk(t)和vk(t)。通過(guò)把系統(tǒng)分解成時(shí)間方向和批次方向的子系統(tǒng)，并且設(shè)計(jì)子系統(tǒng)相應(yīng)的估計(jì)器，抑制了隨機(jī)噪聲和重復(fù)干擾，并提高了間歇過(guò)程的估計(jì)精度。

圖1是ILKF雙向動(dòng)態(tài)特性圖，其中橫軸Time代表時(shí)間方向，縱軸Batch代表批次方向，箭頭表示傳遞方向，此外在每個(gè)采樣時(shí)刻都有對(duì)應(yīng)的測(cè)量值。橫軸上x的轉(zhuǎn)移只是隨著時(shí)間轉(zhuǎn)移，縱軸的x則通過(guò)批次方向的卡爾曼濾波進(jìn)行更新，并且引入了矯正項(xiàng)δx，這個(gè)矯正項(xiàng)也可稱為學(xué)習(xí)項(xiàng)，它在時(shí)間方向上通過(guò)卡爾曼濾波進(jìn)行更新。

圖1 ILKF雙向動(dòng)態(tài)特性圖

2 非線性系統(tǒng)的擬線性化

2.1 動(dòng)機(jī)

針對(duì)文獻(xiàn)[12]中的ILKF存在的非線性問(wèn)題和不合理的初值假設(shè)等，本文將通過(guò)基于標(biāo)稱軌跡的擬線性化方法將ILKF推廣到非線性間歇過(guò)程，并構(gòu)造了一個(gè)誤差系統(tǒng)，通過(guò)對(duì)誤差系統(tǒng)的狀態(tài)估計(jì)來(lái)獲得系統(tǒng)的狀態(tài)估計(jì)，可以放寬初始值的假設(shè)。

考慮非線性間歇過(guò)程系統(tǒng)模型(7)：

(7)

第1節(jié)中針對(duì)干擾的合理假設(shè)1和假設(shè)2在本節(jié)仍然成立；針對(duì)假設(shè)3，本節(jié)系統(tǒng)初值可以分布在任意均值的高斯分布中。

2.2 基于標(biāo)稱軌跡的擬線性化方法

基于標(biāo)稱軌跡的擬線性化方法對(duì)非線性模型(7)進(jìn)行了線性化處理，在本節(jié)中，引入了一個(gè)誤差狀態(tài)系統(tǒng)，誤差系統(tǒng)的狀態(tài)為原系統(tǒng)狀態(tài)軌跡和標(biāo)稱軌跡的偏差，通過(guò)估計(jì)誤差系統(tǒng)的狀態(tài)來(lái)得到非線性系統(tǒng)的狀態(tài)估計(jì)。由于存在干擾和噪聲，每批次運(yùn)行所產(chǎn)生的狀態(tài)軌跡是不同的，但均是在標(biāo)稱軌跡附近的波動(dòng)[13]。

假設(shè)4：實(shí)際軌跡接近于標(biāo)稱軌跡。

(8)

并且定義Δ為狀態(tài)的真實(shí)軌跡和標(biāo)稱軌跡的誤差算子：

(9)

xk(t)=φ(xk(t-1))+d(t-1)+ωk(t-1)=

(10)

基于假設(shè)4，由(9)和(10)可得：

△xk(t)=Φ(t-1)△xk(t-1)+d(t-1)+ωk(t-1)

其中一階偏導(dǎo)近似系數(shù)為：

(11)

同理可得：

△zk(t)=Hk(t)△xk(t)+vk(t)

(12)

其中一階偏導(dǎo)近似系數(shù)Hk(t)為：

(13)

可以得到新的線性化系統(tǒng)：

(14)

Φ(t-1)和H(t)分別如式(11)和式(13)所示。

2.3 時(shí)間和批次方向的子系統(tǒng)

由式(2)中δ的定義易知：

δ△xk(t)=△xk(t)-△xk-1(t)

δ△zk(t)=△zk(t)-△zk-1(t)

工業(yè)生產(chǎn)中重復(fù)非線性過(guò)程一般每一批都是從相同的初始值開始，根據(jù)機(jī)理或者經(jīng)驗(yàn)?zāi)Ｐ?，每一批都具有相同的?biāo)稱軌跡：

(15)

則根據(jù)式(2)和式(9)可得：

δxk(t)=xk(t)-xk-1(t)=△xk(t)-△xk-1(t)

即δxk(t)=δ△xk(t)，結(jié)合式(14)有遞推式：

δ△xk(t)=△xk(t)-△xk-1(t)=

Φ(t-1)δ△xk(t)+δωk(t-1)

同理：

δ△zk(t)=△zk(t)-△zk-1(t)=

H(t)δ△xk(t)+δvk(t)

參考式(3)～(4)將式(14)改寫成時(shí)間和批次兩個(gè)子系統(tǒng)，時(shí)間方向上的子系統(tǒng)∑T：

(16)

批次方向上的子系統(tǒng)∑B：

(17)

(18)

2.4 相關(guān)量的計(jì)算

本小節(jié)計(jì)算迭代學(xué)習(xí)卡爾曼濾波中需要用的量，相關(guān)噪聲協(xié)方差計(jì)算如下：

E[δωk(t)δωkT(t)] =

E[(ωk(t)-ωk-1(t))(ωk(t)-ωk-1(t))T] =

E[ωk(t)ωkT(t)] +E[ωk-1(t)ωkT(t)]=2Q

E[(ωk(t)-ωk-1(t))(ωk-1(t)-ωk-2(t))T]=

所以有：

定義:

由式(14)得：

(19)

(20)

3 迭代學(xué)習(xí)擬線性化卡爾曼濾波

本節(jié)針對(duì)非線性間歇過(guò)程的狀態(tài)估計(jì)問(wèn)題，基于ILKF方法引入標(biāo)稱軌跡擬線性化方法，提出了一種迭代學(xué)習(xí)擬線性化卡爾曼濾波(ILQKF)方法。ILQKF與已有的傳統(tǒng)狀態(tài)估計(jì)方法相比有以下三點(diǎn)的改進(jìn)。

首先，與現(xiàn)有的非線性濾波方法比，現(xiàn)有的非線性濾波方法只考慮了時(shí)間方向的狀態(tài)轉(zhuǎn)移，而ILQKF參考ILKF的濾波框架可以利用之前所有批次的狀態(tài)信息，使得狀態(tài)估計(jì)誤差可以隨時(shí)間和批次兩個(gè)方向減小。

其次，與ILKF方法比，ILKF只針對(duì)線性系統(tǒng)，很難運(yùn)用到非線性間歇生產(chǎn)過(guò)程上，而ILQKF引入擬線性化方法可以解決非線性問(wèn)題，更適合運(yùn)用到實(shí)際生產(chǎn)過(guò)程。

最后，與ILKF方法比，有著更為寬松的初始條件，在ILKF中，由于中間量的計(jì)算需要，文獻(xiàn)[12]給出了假設(shè)3，即原系統(tǒng)的狀態(tài)初值滿足0均值的正態(tài)分布。事實(shí)上，如果不滿足假設(shè)3，ILKF的遞推式將會(huì)變得極復(fù)雜，并且不為0的初值可能會(huì)使得批次子系統(tǒng)中的誤差協(xié)方差變大，從而降低濾波器的性能。然而，ILQKF不必滿足假設(shè)3，其依據(jù)在于擬線性化時(shí)引入了標(biāo)稱模型作為中間變量。當(dāng)系統(tǒng)初值滿足任意高斯分布，將高斯分布的均值作為標(biāo)稱軌跡初值，則新狀態(tài)仍是0均值的正態(tài)分布，即ILQKF可以用于初值期望不為0的非線性系統(tǒng)的狀態(tài)估計(jì)。

3.1 ILQKF估計(jì)結(jié)構(gòu)

與ILKF類似，基于ILQKF設(shè)計(jì)狀態(tài)估計(jì)器分為兩個(gè)階段：時(shí)間方向狀態(tài)估計(jì)和批次方向狀態(tài)估計(jì)。時(shí)間方向狀態(tài)估計(jì)包括對(duì)時(shí)間子系統(tǒng)ΣT的時(shí)間更新和測(cè)量更新；批次方向狀態(tài)估計(jì)包括對(duì)批次子系統(tǒng)ΣB的批次更新和測(cè)量更新。ILQKF的雙向動(dòng)態(tài)特性與圖1不同點(diǎn)在于系統(tǒng)估計(jì)對(duì)象不是原系統(tǒng)狀態(tài)xk(t)，而是狀態(tài)與標(biāo)稱軌跡的誤差Δxk(t)。

時(shí)間子系統(tǒng)ΣT的狀態(tài)估計(jì)器如下：

(21)

批次子系統(tǒng)∑B的狀態(tài)估計(jì)器如下：

(22)

基于式(16)與式(21)得到∑T的估計(jì)誤差系統(tǒng)為：

(23)

基于式(18)與式(22)得到∑B的估計(jì)誤差系統(tǒng)為：

(24)

定義下列誤差協(xié)方差矩陣：

3.2 時(shí)間方向狀態(tài)估計(jì)

1)子系統(tǒng)∑T的時(shí)間更新：

(25)

(26)

2)子系統(tǒng)∑T的測(cè)量更新：

(27)

(28)

(29)

3.3 批次方向狀態(tài)估計(jì)

通過(guò)式(18)、式(22)，可得計(jì)算批次系統(tǒng)的估計(jì)誤差協(xié)方差Ck：

根據(jù)式(23)，式(24)和性質(zhì)1，我們可以計(jì)算得到下面這些相關(guān)協(xié)方差：

批次方向狀態(tài)估計(jì)流程：

1)子系統(tǒng)∑B的批次更新：

(30)

(31)

(32)

2)子系統(tǒng)∑B的測(cè)量更新：

(33)

(34)

(35)

(36)

下面給出ILQKF算法：

算法開始：

輸入：φ(xk(t-1)),h(xk(t)),Q,R,zk(t),M,K；

循環(huán)開始，t=1：

當(dāng)t=M，循環(huán)結(jié)束，否則t=t+1。

濾波開始，k=1：

循環(huán)開始，t=1：

當(dāng)t=M，循環(huán)結(jié)束，否則t=t+1；

當(dāng)k=K，濾波結(jié)束，否則k=k+1；

算法結(jié)束。

4 啤酒發(fā)酵過(guò)程仿真

啤酒發(fā)酵過(guò)程是一種輸入為生物質(zhì)(酵母)，基質(zhì)(糖)和水，輸出為酒精，二氧化碳和水的間歇反應(yīng)過(guò)程。發(fā)酵過(guò)程中，糖類和酵母的濃度都必須控制在合理的范圍之內(nèi)[15]，因此十分有必要對(duì)啤酒發(fā)酵過(guò)程中的糖濃度和酵母濃度進(jìn)行狀態(tài)估計(jì)。

將所提出的ILQKF方法應(yīng)用到啤酒發(fā)酵過(guò)程的步驟如下：

步驟1：根據(jù)已有數(shù)據(jù)或者機(jī)理確定啤酒發(fā)酵的非線性標(biāo)稱模型。

步驟2：初始批生產(chǎn)開始時(shí)，用EKF算法對(duì)糖濃度與酵母濃度進(jìn)行估計(jì)。

步驟3：從第二批生產(chǎn)開始，以上一批的估計(jì)數(shù)據(jù)和標(biāo)稱模型為基礎(chǔ)，用ILQKF算法進(jìn)行估計(jì)，直到生產(chǎn)結(jié)束。

啤酒發(fā)酵過(guò)程中的糖濃度和酵母濃度的狀態(tài)空間模型[16]如下所示：

d(t-1)+ωSk(t-1)

Xk(t-1)+d(t-1)+ωXk(t-1)

zSk(t)=Sk(t)+vSk(t)

zXk(t)=Xk(t)+vXk(t)

其中:Sk(t)，Xk(t)，zSk(t)，zXk(t)，ωk(t)，vk(t)分別代表第k批次、第t時(shí)刻糖濃度，酵母菌絲的濃度，糖濃度的測(cè)量值，酵母菌濃度的測(cè)量值，過(guò)程噪聲，測(cè)量噪聲，ωk(t)和vk(t)相互獨(dú)立。d(t)表示每批只隨時(shí)間變化的重復(fù)性干擾，tc表示采樣間隔。模型參數(shù)如表1所示。

表1 模型參數(shù)表

設(shè)定步長(zhǎng)為[0,3000]，初始值x(0)=xnom(0)=[S(0),X(0)]=[30;3]，白噪聲協(xié)方差矩陣Q=[0.00001,0;0,0.00001]，R=[0.25,0;0,0.0025]，設(shè)置τ=0.96，σ=10-12。

仿真時(shí)，初始批，即k=0批時(shí)，直接采用EKF進(jìn)行狀態(tài)估計(jì)；從k=1批開始，將第0批由EKF得到的狀態(tài)估計(jì)值與標(biāo)稱狀態(tài)之差作為ILQKF第1批的初始估計(jì)值，再采用ILQKF算法得到之后每一批的狀態(tài)估計(jì)值。

圖2顯示了啤酒發(fā)酵中基質(zhì)(糖)和生物質(zhì)(酵母)的狀態(tài)估計(jì)結(jié)果圖，實(shí)線代表實(shí)際值，虛線代表估計(jì)值。計(jì)算第20批基質(zhì)濃度和生物質(zhì)濃度平均估計(jì)誤差分別為0.082 203與0.017 091；計(jì)算第100批基質(zhì)濃度和生物質(zhì)濃度平均估計(jì)誤差分別為0.014 453與0.01 033?？梢园l(fā)現(xiàn)估計(jì)精度隨批次增加而增加。這是由于ILQKF建立了批次間的傳遞關(guān)系，可以過(guò)濾批次間的重復(fù)干擾，隨著批次的增加，濾波使用的信息越來(lái)越多，估計(jì)效果也越來(lái)越好。

圖2 啤酒發(fā)酵狀態(tài)估計(jì)結(jié)果

圖3比較了EKF和ILQKF的狀態(tài)估計(jì)誤差。由于EKF估計(jì)對(duì)于每一批次都是獨(dú)立的，估計(jì)精度只依賴于初始值與動(dòng)態(tài)噪聲的大小，隨著批次的增加，EKF的估計(jì)誤差沒(méi)有很大變化。而ILQKF可以利用上一批次的估計(jì)信息更新初始值減小初始誤差，還同時(shí)考慮了重復(fù)干擾與動(dòng)態(tài)噪聲，隨著批次的增加，ILQKF狀態(tài)估計(jì)誤差逐漸收斂。圖3結(jié)果也表明ILQKF估計(jì)誤差隨批次遞減，而EKF估計(jì)誤差隨批次變化不大，所以ILQKF在多批次運(yùn)行時(shí)的狀態(tài)估計(jì)效果明顯優(yōu)于EKF。

圖3 ILQKF和EKF的狀態(tài)估計(jì)誤差對(duì)比

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

本文基于迭代學(xué)習(xí)和標(biāo)稱軌跡線性化的思想，設(shè)計(jì)了一種用于間歇過(guò)程的迭代學(xué)習(xí)擬線性化卡爾曼濾波器。該ILQKF方法將文獻(xiàn)[12]的ILKF方法擴(kuò)展到非線性過(guò)程中，利用了間歇過(guò)程多個(gè)批次的信息，設(shè)計(jì)了時(shí)間方向和批次方向相應(yīng)的卡爾曼濾波器，得到當(dāng)前批次當(dāng)前時(shí)刻的狀態(tài)估計(jì)值。啤酒發(fā)酵過(guò)程的仿真結(jié)果表明，狀態(tài)估計(jì)誤差隨批次的增加而減小直至收斂，相對(duì)于傳統(tǒng)的估計(jì)方法有一定的優(yōu)勢(shì)。