張治濤， 曾靜

(沈陽(yáng)化工大學(xué)信息工程學(xué)院， 沈陽(yáng) 110142)

非線性系統(tǒng)的狀態(tài)估計(jì)是安全操作的重要先決條件，是過(guò)程監(jiān)控、故障檢測(cè)、過(guò)程優(yōu)化和模型控制的重要組成部分。大型復(fù)雜化學(xué)過(guò)程日益出現(xiàn)在現(xiàn)代加工工業(yè)中，由于其經(jīng)濟(jì)效益，這樣一個(gè)大規(guī)模的復(fù)雜化學(xué)過(guò)程由幾個(gè)操作單元組成，它們通過(guò)物質(zhì)和能量流連接在一起，為降低計(jì)算復(fù)雜度通常需要對(duì)系統(tǒng)進(jìn)行分解[1]。由于過(guò)程規(guī)模的增加和不同子系統(tǒng)之間的顯著相互作用，給這一類(lèi)大規(guī)模復(fù)雜化學(xué)過(guò)程的狀態(tài)估計(jì)提出了巨大的挑戰(zhàn)。

目前，擴(kuò)展卡爾曼濾波(extended Kalman filter，EKF)[2]、無(wú)跡卡爾曼濾波(unscented Kalman filter，UKF)[3]和粒子濾波(particle filter，PF)[4]是解決非線性系統(tǒng)狀態(tài)估計(jì)的主要技術(shù)手段。EKF在現(xiàn)階段有了進(jìn)一步的發(fā)展，文獻(xiàn)[5]中設(shè)計(jì)了遞推最小二乘與自適應(yīng)擴(kuò)展卡爾曼濾波結(jié)合的狀態(tài)估計(jì)，比純EKF算法具有更好的精度與魯棒性。文獻(xiàn)[6]中在EKF的基礎(chǔ)上利用自適應(yīng)插值技術(shù)提高了EKF的估計(jì)精度。文獻(xiàn)[7]中設(shè)計(jì)的分布式EKF與傳統(tǒng)的EKF相比復(fù)雜程度要小很多。EKF對(duì)于非線性程度較強(qiáng)的系統(tǒng)，通過(guò)雅可比線性化的方式引入誤差較大，UKF通過(guò)非線性函數(shù)的概率密度進(jìn)行近似，相比于EKF有更高的計(jì)算精度[8]。PF不受系統(tǒng)模型的限制，具有EKF與UKF無(wú)法比擬的優(yōu)勢(shì)[9-10]，但其計(jì)算量龐大，實(shí)時(shí)性非常差，因此粒子濾波適用的場(chǎng)景非常有限。選擇UKF作為復(fù)雜非線性系統(tǒng)的狀態(tài)觀測(cè)器是不錯(cuò)的選擇。

文獻(xiàn)[11]在卡爾曼濾波理論下，利用分布式最大后驗(yàn)概率估計(jì)方法，與自適應(yīng)飽和度增益結(jié)合設(shè)計(jì)了有適應(yīng)能力的分布式狀態(tài)估計(jì)方法。這種分布式設(shè)計(jì)主要是在線性系統(tǒng)的背景下發(fā)展起來(lái)的，不適合強(qiáng)非線性系統(tǒng)。受到文獻(xiàn)[12-13]的啟發(fā)，對(duì)某一噪聲的突變導(dǎo)致乙苯化工過(guò)程中UKF狀態(tài)估計(jì)準(zhǔn)確性下降問(wèn)題設(shè)計(jì)了一種分布式無(wú)跡卡爾曼(distributed unscented Kalman filter，DUKF)狀態(tài)估計(jì)策略。

基于此，現(xiàn)描述苯與乙烯烷基化產(chǎn)生乙苯的工藝過(guò)程，介紹DUKF的設(shè)計(jì)方法和估計(jì)算法，并與UKF在不同噪聲下比較估計(jì)性能，最后驗(yàn)證存在較大噪聲時(shí)DUKF的優(yōu)勢(shì)。以期為復(fù)雜化工過(guò)程提供了有效的動(dòng)力學(xué)參數(shù)估計(jì)方法。

1 乙苯化工過(guò)程

乙苯是重要的化工原料，乙苯每年的生產(chǎn)總量在千萬(wàn)噸級(jí)別，大多數(shù)的乙苯用來(lái)生產(chǎn)苯乙烯，小部分乙苯被用作稀釋劑、溶劑等。通過(guò)苯和乙烯進(jìn)行烷基化反應(yīng)制得乙苯是最常見(jiàn)的方法[14]。該生產(chǎn)過(guò)程中通過(guò)苯與乙烯反應(yīng)、1,3-二乙苯與苯反應(yīng)產(chǎn)生乙苯的工藝如圖1所示，由4個(gè)連續(xù)攪拌反應(yīng)器(continuous stirred tank reactor，CSTR)和一個(gè)分離器組成，其微分動(dòng)力學(xué)方程詳見(jiàn)文獻(xiàn)[15]。

圖1 乙苯化工過(guò)程流程圖Fig.1 Flow chart of ethylbenzene chemical process

純苯(反應(yīng)物A)通過(guò)物質(zhì)流F1流入CSTR-1，純乙烯(反應(yīng)物B)通過(guò)物質(zhì)流F2、F4、F6分別流入CSTR-1、CSTR-2、CSTR-3，在3個(gè)串聯(lián)的連續(xù)攪拌釜式反應(yīng)器中發(fā)生兩個(gè)反應(yīng)。

反應(yīng)一：苯(反應(yīng)物A)與乙烯(反應(yīng)物B)反應(yīng)產(chǎn)生乙苯(反應(yīng)物C)。

反應(yīng)二：乙烯(反應(yīng)物B)與乙苯(反應(yīng)物C)反應(yīng)產(chǎn)生1,3-二乙基苯(反應(yīng)物D)。

CSTR-3的流出物進(jìn)入閃蒸罐分離器，未來(lái)得及反應(yīng)的苯經(jīng)過(guò)分離器在頂部分離，乙苯產(chǎn)物通過(guò)物質(zhì)流F8流出。頂部的苯通過(guò)Fr一部分進(jìn)入CSTR- 4中一部分重新進(jìn)入CSTR-1中。進(jìn)料流F10是蒸餾過(guò)程產(chǎn)生的1,3-二乙苯(反應(yīng)物D)。在CSTR-4中存在反應(yīng)三。

反應(yīng)三：1,3-二乙苯(反應(yīng)物D)與苯(反應(yīng)物A)反應(yīng)生成乙苯(反應(yīng)物C)。

CSTR-4中所有的產(chǎn)物通過(guò)物質(zhì)流F9進(jìn)入分離器。該工藝的控制輸入是注入或排出Q1、Q2、Q3、Q4和Q55個(gè)容器的熱量，以及CSTR-2和CSTR-3、F4和F6的進(jìn)料流流速。假設(shè)過(guò)程的狀態(tài)對(duì)控制器是連續(xù)可用的。考慮由穩(wěn)態(tài)輸入Q1 s、Q2s、Q3s、Q4s、Q5s、F4s和F6s定義的穩(wěn)定工作點(diǎn)xs，如表1所示。5個(gè)容器中的穩(wěn)態(tài)溫度為：T1s=477.24 K，T2s=476.97 K，T3s=473.47 K，T4s=470.60 K，T5s=478.28 K。

表1 xs的穩(wěn)態(tài)輸入值

該模型可以描述為一個(gè)簡(jiǎn)潔的形式：

(1)

y(t)=Cx(t)+v(t)

(2)

式中：x是過(guò)程的狀態(tài)向量；u為系統(tǒng)的輸入向量；w為過(guò)程擾動(dòng)；y為測(cè)量的輸出向量；C為測(cè)量矩陣；v為測(cè)量噪聲。

將乙苯化工過(guò)程進(jìn)行局部設(shè)計(jì)，分解為兩個(gè)子系統(tǒng)，如圖2所示。每個(gè)UKF都能夠直接獲取相關(guān)子系統(tǒng)的輸出測(cè)量值，并且能夠與其他子系統(tǒng)進(jìn)行數(shù)據(jù)交換，交換的數(shù)據(jù)用于補(bǔ)償子系統(tǒng)之間的相互作用，以改進(jìn)狀態(tài)估計(jì)值。

在進(jìn)行子系統(tǒng)劃分時(shí)除了應(yīng)保證每個(gè)子系統(tǒng)是可觀測(cè)的，還應(yīng)遵循以下準(zhǔn)則：①子系統(tǒng)的分解盡量與系統(tǒng)的物理拓?fù)浣Y(jié)構(gòu)相一致；②每個(gè)操作單元只對(duì)應(yīng)一個(gè)局部狀態(tài)觀測(cè)器。

由兩個(gè)子系統(tǒng)組成的非線性系統(tǒng)其中第i個(gè)子系統(tǒng)可用如下方程描述。

圖2 乙苯化工過(guò)程子系統(tǒng)分解Fig.2 Systematic decomposition of ethylbenzene chemical processes

(3)

yi(t)=Cixi(t)+vi(t)

(4)

2 分布式無(wú)跡卡爾曼濾波

針對(duì)一類(lèi)具有過(guò)程擾動(dòng)和測(cè)量噪聲的非線性系統(tǒng)，設(shè)計(jì)了一種分布式估計(jì)算法，該算法規(guī)定了子系統(tǒng)之間的信息交換方法與DUKF算法的實(shí)現(xiàn)策略。

2.1 DUKF算法步驟

所提出的DUKF通過(guò)以下步驟獲取整個(gè)系統(tǒng)的估計(jì)。

步驟1當(dāng)t=0時(shí)，用初始子系統(tǒng)的狀態(tài)xi(0)(子系統(tǒng)i的狀態(tài))，輸出測(cè)量yi(0)(子系統(tǒng)i的輸出測(cè)量)與協(xié)方差矩陣Pi(0)(子系統(tǒng)i的協(xié)方差)對(duì)所有的UKFi進(jìn)行初始化，Pi(0)可設(shè)置為n×n的正定對(duì)角陣，n為子系統(tǒng)的維數(shù)，以上i=1,2,…，m，其中m為子系統(tǒng)個(gè)數(shù)。

步驟2當(dāng)t>0時(shí)，執(zhí)行以下步驟。

步驟3在下一時(shí)刻進(jìn)行步驟2。

其原理如圖3所示，從上述算法可以看出，在一個(gè)采樣時(shí)間內(nèi)UKF只被并行估計(jì)一次，局部UKF能夠直接獲取相關(guān)系統(tǒng)的測(cè)量值，并且能夠與其他子系統(tǒng)通信，以交換輸出測(cè)量值、狀態(tài)估計(jì)值和方差，這些輸出測(cè)量值、狀態(tài)估計(jì)值和方差用于補(bǔ)償系統(tǒng)之間的相互作用。

圖3 DUKF原理Fig.3 The proposed DUKF design

2.2 局部UKF算法

對(duì)于式(3)和式(4)描述的子系統(tǒng)對(duì)應(yīng)的局部UKF算法步驟如下。

步驟1獲取2n+1個(gè)Sigma采樣點(diǎn)及其權(quán)值。

(5)

(6)

步驟2時(shí)間更新。

(7)

(8)

(9)

(10)

(11)

步驟3測(cè)量更新。

(12)

(13)

(14)

(15)

Pi(k+1|k+1)=Pi(k+1|k)-Ki(k+

1)Pzkzk,iKiT(k+1)

(16)

3 仿真分析

將UKF與DUKF估計(jì)策略應(yīng)用于乙苯化工過(guò)程。系統(tǒng)的25個(gè)狀態(tài)變量如表2所示，測(cè)量輸出變量如表3所示。在仿真過(guò)程中UKF與DUKF使用相同的過(guò)程擾動(dòng)、測(cè)量噪聲以及初始狀態(tài)。過(guò)程擾動(dòng)與觀測(cè)噪聲隨機(jī)產(chǎn)生。過(guò)程擾動(dòng)根據(jù)反應(yīng)物濃度和溫度變化范圍取值。反應(yīng)物A、B、C、D以及溫度的過(guò)程噪聲設(shè)置為均值為0標(biāo)準(zhǔn)差分別為1、0.1、1、0.2、1的正態(tài)分布值。觀測(cè)噪聲設(shè)置為0.01Cx0。x0為反應(yīng)器與分離器的初始值，其參數(shù)如表4所示。UKF協(xié)方差矩陣R和Q為分別設(shè)置為R=diag([1,1,1,0.1,1,1,1,1])，Q=0.01diag[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,100,1,1,1,1,100])。DUKF中R1=diag([1,1,1])，R2=diag([0.1,1,1,1,1])，Q1=0.01diag[1,1,1,1,1,1,1,1,1,1,1,1,1,1])，Q2=0.01diag[1,1,1,1,100,1,1,1,1,100])。仿真中所用輸入為系統(tǒng)穩(wěn)態(tài)輸入值：Q1=-4.4×106J/s，Q2=-4.6×106J/s，Q3=-4.7×106J/s，Q4=9.2×106J/s，Q5=5.9×106J/s，F(xiàn)4=8.697×10-4m3/s，F(xiàn)6=8.697×10-4m3/s。

各反應(yīng)器中分布式UKF和UKF對(duì)真實(shí)狀態(tài)的跟蹤效果如圖4～圖8所示。

表2 乙苯化工過(guò)程中的狀態(tài)變量

表3 乙苯化工過(guò)程中的測(cè)量輸出量

表4 反應(yīng)器與分離器的初始值

K為熱力學(xué)溫度單位開(kāi)爾文圖4 CSTR-1中反應(yīng)物A、B、C、D及溫度的狀態(tài)軌跡Fig.4 State trace of reactants A, B, C, D and temperature in CSTR-1

圖5 CSTR-2中反應(yīng)物A、B、C、D及溫度的狀態(tài)軌跡Fig.5 State trace of reactants A, B, C, D and temperature in CSTR-2

圖6 CSTR-3中反應(yīng)物A、B、C、D及溫度的狀態(tài)軌跡Fig.6 State trace of reactants A, B, C, D and temperature in CSTR-3

圖7 分離器中反應(yīng)物A、B、C、D及溫度的狀態(tài)軌跡Fig.7 State trajectories of reactants A, B, C, D and temperature in the separator

圖8 CSTR-4中反應(yīng)物A、B、C、D及溫度的狀態(tài)軌跡Fig.8 State trace of reactants A, B, C, D and temperature in CSTR-4

通過(guò)仿真實(shí)驗(yàn)得出UKF與DUKF都可以較好地估計(jì)出乙苯化工過(guò)程的真實(shí)狀態(tài)。為更準(zhǔn)確地評(píng)估兩者的性能，引入歐幾里得范數(shù)對(duì)兩者的誤差進(jìn)一步分析?？紤]不同狀態(tài)下估計(jì)誤差的大小，根據(jù)這兩種估計(jì)方案的最大估計(jì)誤差對(duì)每個(gè)狀態(tài)的誤差進(jìn)行歸一化處理。歸一化處理后的歐幾里得范數(shù)定義為

(17)

式(17)中：k為時(shí)刻；i為系統(tǒng)的25個(gè)狀態(tài)。兩種算法中每個(gè)狀態(tài)的最大誤差由整個(gè)采樣周期的狀態(tài)決定。對(duì)歐幾里得范數(shù)取平均值，即

(18)

式(18)中：T為仿真時(shí)長(zhǎng)。

UKF與DUKF估計(jì)誤差平均值分別用Umean與DUmean表示。在本組實(shí)驗(yàn)數(shù)據(jù)中UKF與DUKF的歐幾里得范數(shù)平均值分別為2.62、2.50，最大值方面UKF、DUKF分別為4.36、3.68。分布式的設(shè)計(jì)并沒(méi)有削弱估計(jì)性能，在當(dāng)前誤差水平下DUKF的估計(jì)準(zhǔn)確度比UKF更高并且DUKF有更穩(wěn)定的表現(xiàn)。圖9給出了兩種算法的歸一化處理后的歐幾里得范數(shù)軌跡。以本組過(guò)程擾動(dòng)標(biāo)準(zhǔn)差為基準(zhǔn)通過(guò)增加過(guò)程噪聲標(biāo)準(zhǔn)差系數(shù)來(lái)獲取多組不同的擾動(dòng)。由表5可知，無(wú)論是過(guò)程擾動(dòng)較大時(shí)(前3組數(shù)據(jù))，還是過(guò)程擾動(dòng)較小時(shí)(后2組)DUKF的平均估計(jì)誤差都要優(yōu)于UKF。

在子系統(tǒng)2中的分離器引入一個(gè)特別大的測(cè)量噪聲，使反應(yīng)物D的測(cè)量噪聲由原來(lái)的7.2突變?yōu)?20來(lái)模擬傳感器受到強(qiáng)烈影響的誤差變化。當(dāng)引入強(qiáng)烈噪聲后子系統(tǒng)1中各反應(yīng)物狀態(tài)估計(jì)情況如圖10～圖12所示。

圖9 歐幾里得范數(shù)軌跡Fig.9 The trajectory of the Euclidean norm of the normalized estimate error of UKF

表5 反應(yīng)物與溫度在不同噪聲下的誤差

當(dāng)傳感器受到強(qiáng)干擾時(shí)，UKF的整體估計(jì)精度受到了很大的影響，DUKF可以很好地處理強(qiáng)干擾所帶來(lái)的影響，這是由于系統(tǒng)的局部處理設(shè)計(jì)把噪聲限制在了當(dāng)前子系統(tǒng)中。值得注意的是，當(dāng)誤差水平接近于理想狀態(tài)下時(shí)DUKF的估計(jì)性能低于UKF，這主要是系統(tǒng)局部處理時(shí)產(chǎn)生的誤差導(dǎo)致。對(duì)局部處理后的系統(tǒng)越接近于原系統(tǒng)兩者在理想情況下的估計(jì)結(jié)果越接近。

4 結(jié)論

為解決非線性系統(tǒng)的狀態(tài)估計(jì)問(wèn)題提出了一種基于DUKF 的狀態(tài)估計(jì)策略，并在乙苯化工過(guò)程中得到了驗(yàn)證。仿真結(jié)果表明，分布式UKF雖然對(duì)系統(tǒng)進(jìn)行了分解，但是在復(fù)雜非線性系統(tǒng)中擁有比集中式UKF更好地估計(jì)精度。DUKF相對(duì)于UKF有更高的穩(wěn)定性。當(dāng)傳感器受到強(qiáng)干擾時(shí)依然有較好的估計(jì)性能，容錯(cuò)性更強(qiáng)。DUKF的子系統(tǒng)設(shè)計(jì)使得系統(tǒng)有動(dòng)態(tài)可變的特點(diǎn)，根據(jù)實(shí)際情況改變子系統(tǒng)分配實(shí)現(xiàn)工廠效益最大化。

在今后的工作中，將對(duì)濾波器傳輸頻率進(jìn)行改善，在保證估計(jì)準(zhǔn)確性的同時(shí)降低信息傳輸頻率。

圖10 加入較大測(cè)量噪聲后CSTR-1中反應(yīng)物A、B、C、D及溫度的狀態(tài)軌跡Fig.10 State trace of reactants A, B, C, D and temperature in CSTR-1 after adding large measurement noise

圖11 加入較大測(cè)量噪聲后CSTR-2中反應(yīng)物A、B、C、D及溫度的狀態(tài)軌跡Fig.11 State trace of reactants A, B, C, D and temperature in CSTR-2 after adding large measurement noise

圖12 加入較大測(cè)量噪聲后CSTR-3中反應(yīng)物A、B、C、D及溫度的狀態(tài)軌跡Fig.12 State trace of reactants A, B, C, D and temperature in CSTR-3 after adding large measurement noise