周 鑫 趙 眾 孫 康

(1.北京化工大學(xué)信息科學(xué)與技術(shù)學(xué)院，北京 100029；2.中國石油化工股份有限公司北京化工研究院，北京 100013)

在氣相聚乙烯[1]的工業(yè)生產(chǎn)過程中，有兩個非常重要的產(chǎn)品指標(biāo)——熔融指數(shù)和密度，它們能實現(xiàn)對產(chǎn)品質(zhì)量指標(biāo)的精確估計，有助于提高聚乙烯工業(yè)裝置優(yōu)化控制的結(jié)果。目前工業(yè)應(yīng)用的在線測量儀表存在不準(zhǔn)確的問題，聚乙烯產(chǎn)品質(zhì)量主要通過定時采樣得到分析數(shù)據(jù)，再根據(jù)分析數(shù)據(jù)來調(diào)控以實現(xiàn)對質(zhì)量指標(biāo)的控制。這樣的操控方式使得聚乙烯生產(chǎn)過程中產(chǎn)品質(zhì)量會出現(xiàn)波動，嚴(yán)重時可能造成停車。在實際的氣相聚乙烯工業(yè)生產(chǎn)中，氣相法生產(chǎn)聚乙烯的牌號有多種，以此來滿足市場對不同產(chǎn)品的需求，在實際生產(chǎn)過程中常常需要進行多種牌號間的相互切換。牌號切換會造成聚乙烯產(chǎn)品質(zhì)量的波動，同時也會造成定時采樣分析數(shù)據(jù)過程不準(zhǔn)確。筆者針對該生產(chǎn)過程進行質(zhì)量指標(biāo)的在線估計，有助于提高企業(yè)的產(chǎn)品質(zhì)量和經(jīng)濟效益。

1 聚乙烯質(zhì)量指標(biāo)動態(tài)模型①

在Unipol工藝的聚乙烯工業(yè)裝置中，通常用熔融指數(shù)MI代表聚合物的分子量，密度Den代表共聚組成。實際生產(chǎn)過程中，聚合物的分子量受到諸多因素的影響，當(dāng)催化劑確定后，產(chǎn)品的質(zhì)量主要受以下因素影響：聚合溫度、氫氣/乙烯比、共聚單體/乙烯比。一般提高溫度，可以加快鏈轉(zhuǎn)移反應(yīng)，使聚合物分子量下降，熔融指數(shù)升高；氫氣作為鏈轉(zhuǎn)移劑加入反應(yīng)體系中主要用來控制熔融指數(shù)，氫烯比的提高，使得樹脂分子量降低而樹脂的熔融指數(shù)隨之提高。同時，熔融指數(shù)還與聚合溫度和共聚單體的濃度有關(guān)。聚合溫度關(guān)系到催化劑的活性，共聚單體的加入，可以增加聚乙烯分子鏈的支鏈數(shù)目。因此，熔融指數(shù)主要由氫氣濃度進行調(diào)節(jié)，其他參數(shù)可作為熔融指數(shù)MI的微調(diào)參數(shù)。建立熔融指數(shù)機理模型如下[2]：

lnMIc(i+1)=fr(i)lnMIc(i)+[1-fr(i)][θ0+

(1)

聚合物的密度取決于聚合物鏈的結(jié)構(gòu)參數(shù)，可用一般的經(jīng)驗公式表示[3]：

Denc(i+1)=fr(i)Denc(i)+[1-fr(i)]·

(2)

根據(jù)實際經(jīng)驗，式(2)可進一步簡化為：

Denc(i+1)=fr(i)Denc(i)+[1-fr(i)]·

(3)

其中，θ0～θ8為模型參數(shù)，[C2]、[Cx]、[H2]分別為乙烯濃度、共聚單體濃度和氫氣濃度，fr=exp(-Δt/τ)，lnMIc、Denc分別為lnMI、Den累積值，Δt為采樣時間，τ為顆粒在反應(yīng)器中的停留時間。

2 質(zhì)量指標(biāo)在線魯棒濾波估計

用下面的非線性動態(tài)模型來描述式(1)、(3)的質(zhì)量指標(biāo)預(yù)測模型：

x(k+1)=A(k)x(k)+B(k)F(u(k))y(k)=Cx(k)

(4)

其中x為狀態(tài)向量，x∈Rn；u為控制輸入向量，u∈Rq；y為輸出向量，y∈Rm；F為辨識得到的非線性穩(wěn)態(tài)函數(shù)，F(xiàn)∶Rq→Rn；C為測量矩陣，C∶Rm×n。

實際非線性動態(tài)過程可表示為：

x(k+1)=f(x(k),u(k))+ξ1(k)y(k)=Cx(k)+ξ2(k)

(5)

其中f∶Rn×Rq→Rn為未知非線性函數(shù)，ξ1(k)、ξ2(k)分別為n維和m維高斯白噪聲向量，且具有如下的統(tǒng)計性質(zhì)：

Eξ1(k)=Eξ2(k)=0，E[ξ1(k)ξ2(k)T]=0

(6)

E[ξ1(k)ξ1(j)T]=Qδk,j，E[ξ2(k)ξ2(j)T]=R(Δ)δk,j

(7)

其中Q為模型噪聲方差，是正定對稱矩陣；E為數(shù)學(xué)期望符號；δk,j為Kronecker函數(shù)；R(Δ)為測量噪聲方差，具有如下不確定性：

R-1(Δ)=R-1+MΔN+NTΔMT

(8)

其中Δ∈{Δ∶‖Δ‖≤1}，表征多牌號切換時由于工藝變動引起的不確定性。

濾波估計如下[4]：

(9)

其中K(k)為濾波器增益矩陣,K(k)∈Rn×m。

定義 預(yù)測模型的真實最大逼近偏差為：

ef≡‖f(x(k),u(k))-A(k)x(k)-B(k)F(u(k))‖∞

(10)

估計誤差偏差為：

K(k)ξ2(k)

(11)

預(yù)測模型誤差為：

(12)

估計誤差協(xié)方差陣為：

[A(k)-K(k)C]ρ(k)[A(k)-K(k)C]T+QI+K(k)R(Δ)KT(k)

(13)

引理1d∈Rn為一列向量，I為單位陣，有ddT≥0；dTdI≥ddT。

由引理1，則以下的矩陣不等式成立：

(14)

≤[A(k)-K(k)C]ρ(k)[A(k)-K(k)C]T+nef2I

(15)

把式(13)化為：

ρ(k+1)≤2nef2I+2[A(k)-K(k)C]P(k)[A(k)-

K(k)C]T+QI+K(k)R(Δ)K(k)T

(16)

取P(k+1)為ρ(k+1)的上限估計，其中P(0)≥ρ(0)，且P(k+1)滿足遞推關(guān)系：

P(k+1)=2nef2I+2[A(k)-K(k)C]P(k)[A(k)-

K(k)C]T+QI+K(k)R(Δ)K(k)T

(17)

易證：P(k+1)≥ρ(k+1)。

通過極小化誤差協(xié)方差上限大小，優(yōu)化得到狀態(tài)估計器，其中估計器增益K(k)和誤差協(xié)方差上限矩陣P(k+1)滿足：

(18a)

2nef2I+2[A(k)-K(k)C]P(k)[A(k)-K(k)C]T+

QI+K(k)R(Δ)K(k)T

(18b)

其中‖P(k+1)‖∞為估計誤差協(xié)方差上限的無窮范數(shù)。

(19)

引理3[6]給定具有適當(dāng)維數(shù)的實矩陣F=FT、M、N，當(dāng)且僅當(dāng)存在λ使得下式成立：

(20)

則對于所有的Δ，‖Δ‖≤1，F(xiàn)+MΔN+NTΔTMT>0成立。

定理1 優(yōu)化問題式(18)可用如下線性矩陣不等式(LMI)求解：

(21)

(22)

證明

式(18b)可轉(zhuǎn)化為：

(23)

由引理3，式(23)可轉(zhuǎn)化為：

(24)

將R-1(Δ)不確定形式代入，得：

(25)

由引理3，得：

(26)

綜上，基于LMI的濾波器增益矩陣K(k)的遞推算法如下：

a. 給定初始誤差協(xié)方差矩陣P(0)、噪聲方差Q、R(Δ)、初始狀態(tài)x(0)；

b. 在采樣k時刻，求解式(21)、(22)，得到最壞情況下的遞推矩陣P*(k+1)和K(k)；

c. 根據(jù)協(xié)方差上限遞推矩陣(17)，計算k+1時刻誤差協(xié)方差上限矩陣P(k+1)，重復(fù)步驟b。

3 實際應(yīng)用

將筆者提出的魯棒濾波方法應(yīng)用在中石化某氣相法聚乙烯工業(yè)生產(chǎn)裝置上，生產(chǎn)工藝流程如圖1所示。乙烯、丁烯(或己烯)、氫氣、氮氣和其他惰性氣體一起被加入反應(yīng)系統(tǒng)組成循環(huán)氣體。循環(huán)氣體通過反應(yīng)器，主要是流化聚乙烯顆粒，同時把部分乙烯聚合(轉(zhuǎn)化率為1%～2%)。循環(huán)氣體從反應(yīng)器出來之后先進入循環(huán)氣體壓縮機，再經(jīng)過換熱器，最后再重新回到循環(huán)反應(yīng)器中。精制乙烯、氫氣、丁烯(或己烯)、三乙基鋁助催化劑從流化床循環(huán)管路進入聚合系統(tǒng)。把催化劑直接注入到流化床反應(yīng)器內(nèi)，反應(yīng)器的出料機構(gòu)把聚合物產(chǎn)品輸出，并且保持反應(yīng)器的料位高度恒定。

圖1 氣相聚乙烯生產(chǎn)工藝流程簡圖

設(shè)定質(zhì)量指標(biāo)預(yù)測模型辨識目標(biāo)函數(shù)為：

(27)

運用粒子群優(yōu)化算法[7]，選擇牌號DGM-1815的數(shù)據(jù)，通過求解式(27)，得到質(zhì)量指標(biāo)預(yù)測模型參數(shù)(表1)，模型最大偏差和均方根誤差見表2，模型辨識效果如圖2所示。由圖2可以看出，所得到的模型變化趨勢與實際生產(chǎn)過程中質(zhì)量指標(biāo)的變化趨勢相符。

表1 質(zhì)量指標(biāo)預(yù)測模型參數(shù)

表2 模型辨識與泛化效果

圖2 DGM-1815熔融指數(shù)和密度的預(yù)測模型與實際數(shù)據(jù)比較

獲得模型參數(shù)后，筆者對不同牌號的DGM-1830進行了模型參數(shù)泛化驗證，驗證效果如圖3所示。模型最大偏差和均方根誤差見表2。從圖3和表2可以看出，所得動態(tài)模型反映了實際生產(chǎn)過程中熔融指數(shù)和密度的變化趨勢，但對于多牌號的泛化會引起模型預(yù)測誤差的顯著增加，因而根據(jù)實驗室分析數(shù)據(jù)進行濾波修正是非常必要的。

圖3 DGM-1830熔融指數(shù)和密度的預(yù)測模型與實際數(shù)據(jù)比較

美國聯(lián)碳校正算法為：

(28)

EMXf=λEXX+(1-λ)EMXfold

(29)

其中EMX為瞬時模型偏差；EMXfold為上一時刻濾波后的模型偏差；λ為模型偏差的濾波系數(shù)(缺省取0.4～0.8)。此外還需要一些條件來對校正進行選擇：RX>8或EMX>3σX。所提方法與美國聯(lián)碳應(yīng)用比較如圖4所示。最大方差和均根誤差比較結(jié)果見表3。從圖4可以看出，所提方法的質(zhì)量指標(biāo)曲線比美國聯(lián)碳校正方法的質(zhì)量指標(biāo)曲線更加貼近實際值，有效地修正了多牌號模型預(yù)測的偏差。

表3 校正效果比較

所提方法已在中石化某氣相聚乙烯工業(yè)裝置上成功投用，DCS上實際投用效果曲線如圖5所示。從投用效果圖可以看出，所提方法真實反映了多牌號的質(zhì)量指標(biāo)變化趨勢，達到了實時估計的目的，為質(zhì)量指標(biāo)閉環(huán)優(yōu)化控制提供了良好的條件。

圖5 熔融指數(shù)和密度的DCS實際投用效果

4 結(jié)束語

針對氣相聚乙烯工業(yè)裝置多牌號質(zhì)量指標(biāo)實時估計的復(fù)雜性，在聚乙烯工業(yè)裝置質(zhì)量指標(biāo)實時預(yù)測模型的基礎(chǔ)上，提出了一種氣相聚乙烯工業(yè)裝置多牌號質(zhì)量指標(biāo)魯棒濾波器設(shè)計方法。利用該方法根據(jù)實驗室分析數(shù)據(jù)反饋修正模型預(yù)測并實時估計質(zhì)量指標(biāo)。該方法將模型逼近偏差顯式地引入到濾波器設(shè)計中，并考慮了由多牌號切換引起的工藝不確定性，利用LMI和最小方差估計原理進行濾波器增益的實時求解。該方法在氣相聚乙烯工業(yè)裝置上的應(yīng)用結(jié)果和與美國聯(lián)碳技術(shù)的應(yīng)用對比證實了其有效性和可行性，為實現(xiàn)氣相聚乙烯工業(yè)裝置多牌號質(zhì)量指標(biāo)的閉環(huán)優(yōu)化控制創(chuàng)造了良好的條件。

[1] 鄧哲文．氣相流化床法聚乙烯工藝技術(shù)比較[J].化工設(shè)計,2006,16(4):9～13．

[2] McAuley K B,MacGregor J F．A Kinetic Model for Industrial Gas Phase Ethylene Copolymerization[J]．AIChE Journal,1990,36(6):837～844．

[3] Debling J,Kuijpers F,Ray W H．Dynamic Modeling of Product Grade Transition for Olefin Polymerization Process[J]．AIChE Journal,1994,40(3):506～513．

[4] 趙眾，馬博．大型聚乙烯工業(yè)裝置質(zhì)量指標(biāo)的次優(yōu)強跟蹤濾波估計[J]．化工學(xué)報,2008,7(15):1635～1639．

[5] Boyd S,Ghaoui L EI,Feron E．Linear Matrix Inequalities in Systems and Control Theory[M]．Philadelphia,PA:SIAM,1994．

[6] Calafiore G,Ghaoui L EI．Minimum Variance Estimation with Uncertain Statistical Model[C]．Proceedings of the 40th IEEE Conference on Decision and Control.Orlando, Florida USA:IEEE,2001:3497～3499．

[7] Ratnaweera A,Halgamuge S K,Watson H C．Self-organizing Hierarchical Particle Swarm Optimizer with Time-varying Acceleration Coefficients [J]．IEEE Transactions on Evolutionary Computation,2004,8(3):240～255．