周 鑫 趙 眾 孫 康
(1.北京化工大學(xué)信息科學(xué)與技術(shù)學(xué)院,北京 100029;2.中國(guó)石油化工股份有限公司北京化工研究院,北京 100013)
在氣相聚乙烯[1]的工業(yè)生產(chǎn)過(guò)程中,有兩個(gè)非常重要的產(chǎn)品指標(biāo)——熔融指數(shù)和密度,它們能實(shí)現(xiàn)對(duì)產(chǎn)品質(zhì)量指標(biāo)的精確估計(jì),有助于提高聚乙烯工業(yè)裝置優(yōu)化控制的結(jié)果。目前工業(yè)應(yīng)用的在線測(cè)量?jī)x表存在不準(zhǔn)確的問(wèn)題,聚乙烯產(chǎn)品質(zhì)量主要通過(guò)定時(shí)采樣得到分析數(shù)據(jù),再根據(jù)分析數(shù)據(jù)來(lái)調(diào)控以實(shí)現(xiàn)對(duì)質(zhì)量指標(biāo)的控制。這樣的操控方式使得聚乙烯生產(chǎn)過(guò)程中產(chǎn)品質(zhì)量會(huì)出現(xiàn)波動(dòng),嚴(yán)重時(shí)可能造成停車。在實(shí)際的氣相聚乙烯工業(yè)生產(chǎn)中,氣相法生產(chǎn)聚乙烯的牌號(hào)有多種,以此來(lái)滿足市場(chǎng)對(duì)不同產(chǎn)品的需求,在實(shí)際生產(chǎn)過(guò)程中常常需要進(jìn)行多種牌號(hào)間的相互切換。牌號(hào)切換會(huì)造成聚乙烯產(chǎn)品質(zhì)量的波動(dòng),同時(shí)也會(huì)造成定時(shí)采樣分析數(shù)據(jù)過(guò)程不準(zhǔn)確。筆者針對(duì)該生產(chǎn)過(guò)程進(jìn)行質(zhì)量指標(biāo)的在線估計(jì),有助于提高企業(yè)的產(chǎn)品質(zhì)量和經(jīng)濟(jì)效益。
在Unipol工藝的聚乙烯工業(yè)裝置中,通常用熔融指數(shù)MI代表聚合物的分子量,密度Den代表共聚組成。實(shí)際生產(chǎn)過(guò)程中,聚合物的分子量受到諸多因素的影響,當(dāng)催化劑確定后,產(chǎn)品的質(zhì)量主要受以下因素影響:聚合溫度、氫氣/乙烯比、共聚單體/乙烯比。一般提高溫度,可以加快鏈轉(zhuǎn)移反應(yīng),使聚合物分子量下降,熔融指數(shù)升高;氫氣作為鏈轉(zhuǎn)移劑加入反應(yīng)體系中主要用來(lái)控制熔融指數(shù),氫烯比的提高,使得樹(shù)脂分子量降低而樹(shù)脂的熔融指數(shù)隨之提高。同時(shí),熔融指數(shù)還與聚合溫度和共聚單體的濃度有關(guān)。聚合溫度關(guān)系到催化劑的活性,共聚單體的加入,可以增加聚乙烯分子鏈的支鏈數(shù)目。因此,熔融指數(shù)主要由氫氣濃度進(jìn)行調(diào)節(jié),其他參數(shù)可作為熔融指數(shù)MI的微調(diào)參數(shù)。建立熔融指數(shù)機(jī)理模型如下[2]:
lnMIc(i+1)=fr(i)lnMIc(i)+[1-fr(i)][θ0+
(1)
聚合物的密度取決于聚合物鏈的結(jié)構(gòu)參數(shù),可用一般的經(jīng)驗(yàn)公式表示[3]:
Denc(i+1)=fr(i)Denc(i)+[1-fr(i)]·
(2)
根據(jù)實(shí)際經(jīng)驗(yàn),式(2)可進(jìn)一步簡(jiǎn)化為:
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為采樣時(shí)間,τ為顆粒在反應(yīng)器中的停留時(shí)間。
用下面的非線性動(dòng)態(tài)模型來(lái)描述式(1)、(3)的質(zhì)量指標(biāo)預(yù)測(cè)模型:
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為辨識(shí)得到的非線性穩(wěn)態(tài)函數(shù),F(xiàn)∶Rq→Rn;C為測(cè)量矩陣,C∶Rm×n。
實(shí)際非線性動(dòng)態(tài)過(guò)程可表示為:
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)計(jì)性質(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為模型噪聲方差,是正定對(duì)稱矩陣;E為數(shù)學(xué)期望符號(hào);δk,j為Kronecker函數(shù);R(Δ)為測(cè)量噪聲方差,具有如下不確定性:
R-1(Δ)=R-1+MΔN+NTΔMT
(8)
其中Δ∈{Δ∶‖Δ‖≤1},表征多牌號(hào)切換時(shí)由于工藝變動(dòng)引起的不確定性。
濾波估計(jì)如下[4]:
(9)
其中K(k)為濾波器增益矩陣,K(k)∈Rn×m。
定義 預(yù)測(cè)模型的真實(shí)最大逼近偏差為:
ef≡‖f(x(k),u(k))-A(k)x(k)-B(k)F(u(k))‖∞
(10)
估計(jì)誤差偏差為:
K(k)ξ2(k)
(11)
預(yù)測(cè)模型誤差為:
(12)
估計(jì)誤差協(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)的上限估計(jì),其中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)。
通過(guò)極小化誤差協(xié)方差上限大小,優(yōu)化得到狀態(tài)估計(jì)器,其中估計(jì)器增益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)‖∞為估計(jì)誤差協(xié)方差上限的無(wú)窮范數(shù)。
(19)
引理3[6]給定具有適當(dāng)維數(shù)的實(shí)矩陣F=FT、M、N,當(dāng)且僅當(dāng)存在λ使得下式成立:
(20)
則對(duì)于所有的Δ,‖Δ‖≤1,F(xiàn)+MΔN+NTΔTMT>0成立。
定理1 優(yōu)化問(wèn)題式(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時(shí)刻,求解式(21)、(22),得到最壞情況下的遞推矩陣P*(k+1)和K(k);
c. 根據(jù)協(xié)方差上限遞推矩陣(17),計(jì)算k+1時(shí)刻誤差協(xié)方差上限矩陣P(k+1),重復(fù)步驟b。
將筆者提出的魯棒濾波方法應(yīng)用在中石化某氣相法聚乙烯工業(yè)生產(chǎn)裝置上,生產(chǎn)工藝流程如圖1所示。乙烯、丁烯(或己烯)、氫氣、氮?dú)夂推渌栊詺怏w一起被加入反應(yīng)系統(tǒng)組成循環(huán)氣體。循環(huán)氣體通過(guò)反應(yīng)器,主要是流化聚乙烯顆粒,同時(shí)把部分乙烯聚合(轉(zhuǎn)化率為1%~2%)。循環(huán)氣體從反應(yīng)器出來(lái)之后先進(jìn)入循環(huán)氣體壓縮機(jī),再經(jīng)過(guò)換熱器,最后再重新回到循環(huán)反應(yīng)器中。精制乙烯、氫氣、丁烯(或己烯)、三乙基鋁助催化劑從流化床循環(huán)管路進(jìn)入聚合系統(tǒng)。把催化劑直接注入到流化床反應(yīng)器內(nèi),反應(yīng)器的出料機(jī)構(gòu)把聚合物產(chǎn)品輸出,并且保持反應(yīng)器的料位高度恒定。
圖1 氣相聚乙烯生產(chǎn)工藝流程簡(jiǎn)圖
設(shè)定質(zhì)量指標(biāo)預(yù)測(cè)模型辨識(shí)目標(biāo)函數(shù)為:
(27)
運(yùn)用粒子群優(yōu)化算法[7],選擇牌號(hào)DGM-1815的數(shù)據(jù),通過(guò)求解式(27),得到質(zhì)量指標(biāo)預(yù)測(cè)模型參數(shù)(表1),模型最大偏差和均方根誤差見(jiàn)表2,模型辨識(shí)效果如圖2所示。由圖2可以看出,所得到的模型變化趨勢(shì)與實(shí)際生產(chǎn)過(guò)程中質(zhì)量指標(biāo)的變化趨勢(shì)相符。
表1 質(zhì)量指標(biāo)預(yù)測(cè)模型參數(shù)
表2 模型辨識(shí)與泛化效果
圖2 DGM-1815熔融指數(shù)和密度的預(yù)測(cè)模型與實(shí)際數(shù)據(jù)比較
獲得模型參數(shù)后,筆者對(duì)不同牌號(hào)的DGM-1830進(jìn)行了模型參數(shù)泛化驗(yàn)證,驗(yàn)證效果如圖3所示。模型最大偏差和均方根誤差見(jiàn)表2。從圖3和表2可以看出,所得動(dòng)態(tài)模型反映了實(shí)際生產(chǎn)過(guò)程中熔融指數(shù)和密度的變化趨勢(shì),但對(duì)于多牌號(hào)的泛化會(huì)引起模型預(yù)測(cè)誤差的顯著增加,因而根據(jù)實(shí)驗(yàn)室分析數(shù)據(jù)進(jìn)行濾波修正是非常必要的。
圖3 DGM-1830熔融指數(shù)和密度的預(yù)測(cè)模型與實(shí)際數(shù)據(jù)比較
美國(guó)聯(lián)碳校正算法為:
(28)
EMXf=λEXX+(1-λ)EMXfold
(29)
其中EMX為瞬時(shí)模型偏差;EMXfold為上一時(shí)刻濾波后的模型偏差;λ為模型偏差的濾波系數(shù)(缺省取0.4~0.8)。此外還需要一些條件來(lái)對(duì)校正進(jìn)行選擇:RX>8或EMX>3σX。所提方法與美國(guó)聯(lián)碳應(yīng)用比較如圖4所示。最大方差和均根誤差比較結(jié)果見(jiàn)表3。從圖4可以看出,所提方法的質(zhì)量指標(biāo)曲線比美國(guó)聯(lián)碳校正方法的質(zhì)量指標(biāo)曲線更加貼近實(shí)際值,有效地修正了多牌號(hào)模型預(yù)測(cè)的偏差。
表3 校正效果比較
所提方法已在中石化某氣相聚乙烯工業(yè)裝置上成功投用,DCS上實(shí)際投用效果曲線如圖5所示。從投用效果圖可以看出,所提方法真實(shí)反映了多牌號(hào)的質(zhì)量指標(biāo)變化趨勢(shì),達(dá)到了實(shí)時(shí)估計(jì)的目的,為質(zhì)量指標(biāo)閉環(huán)優(yōu)化控制提供了良好的條件。
圖5 熔融指數(shù)和密度的DCS實(shí)際投用效果
針對(duì)氣相聚乙烯工業(yè)裝置多牌號(hào)質(zhì)量指標(biāo)實(shí)時(shí)估計(jì)的復(fù)雜性,在聚乙烯工業(yè)裝置質(zhì)量指標(biāo)實(shí)時(shí)預(yù)測(cè)模型的基礎(chǔ)上,提出了一種氣相聚乙烯工業(yè)裝置多牌號(hào)質(zhì)量指標(biāo)魯棒濾波器設(shè)計(jì)方法。利用該方法根據(jù)實(shí)驗(yàn)室分析數(shù)據(jù)反饋修正模型預(yù)測(cè)并實(shí)時(shí)估計(jì)質(zhì)量指標(biāo)。該方法將模型逼近偏差顯式地引入到濾波器設(shè)計(jì)中,并考慮了由多牌號(hào)切換引起的工藝不確定性,利用LMI和最小方差估計(jì)原理進(jìn)行濾波器增益的實(shí)時(shí)求解。該方法在氣相聚乙烯工業(yè)裝置上的應(yīng)用結(jié)果和與美國(guó)聯(lián)碳技術(shù)的應(yīng)用對(duì)比證實(shí)了其有效性和可行性,為實(shí)現(xiàn)氣相聚乙烯工業(yè)裝置多牌號(hào)質(zhì)量指標(biāo)的閉環(huán)優(yōu)化控制創(chuàng)造了良好的條件。
[1] 鄧哲文.氣相流化床法聚乙烯工藝技術(shù)比較[J].化工設(shè)計(jì),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)強(qiáng)跟蹤濾波估計(jì)[J].化工學(xué)報(bào),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.