面向高密度聚乙烯生產(chǎn)過程的最優(yōu)采樣策略

尹 鵬, 陳偉鋒

(浙江工業(yè)大學 信息工程學院, 浙江 杭州 310023)

1 前 言

高密度聚乙烯(high-density polyethylene，HDPE)是5 大通用合成樹脂之一，具有優(yōu)異的化學穩(wěn)定性和耐低溫性。我國是HPDE 最大消費國，也是最大進口國[1]，2019 年消費量為15 478 kt，其中進口7 997 kt，產(chǎn)品自足率為49.39%。目前我國HDPE 產(chǎn)品結(jié)構(gòu)矛盾突出，中低端產(chǎn)能過剩，高端產(chǎn)品仍以進口為主。高密度聚乙烯實際上是混合物，由不同鏈長的高分子聚合物組成，組成成分不同會帶來不同的材料特性，為對聚合物材料特性進行表征，質(zhì)量分數(shù)的分布常常作為質(zhì)量指標。質(zhì)量分數(shù)的分布為所有不同鏈長分子質(zhì)量占聚合物總質(zhì)量的分數(shù)，它可以反映產(chǎn)品的微觀性質(zhì)[2]。高密度聚乙烯模型的開發(fā)對提升產(chǎn)品的質(zhì)量有極為重要的意義，而參數(shù)估計是實現(xiàn)基于質(zhì)量分數(shù)分布的高密度聚乙烯建模的核心技術(shù)手段之一。

機理模型在反應過程的設計、優(yōu)化和控制中起著重要的作用，參數(shù)估計是建立過程模型的關(guān)鍵組成部分。然而，由于參數(shù)可識別性和可估計性的問題，高質(zhì)量參數(shù)值的估計并不容易[3]。從另一個角度來看，參數(shù)估計問題往往是不定的，即不是所有的參數(shù)都可以唯一確定。常用的方法是通過將其他參數(shù)固定在合理的值上來估計參數(shù)的一個子集，而難點就在于確定可估計性高的參數(shù)子集。Yao 等[4]對參數(shù)子集選擇方法進行了研究，其中正交化方法具有明顯的幾何含義并且易于實現(xiàn)，它將靈敏度矩陣S正交化排序，可以看作是每一步驟的D 準則最大化的順序方法，有效性被Chu 等[5]證明。正交化方法是參數(shù)子集選擇中應用最廣泛的技術(shù)[6-9]，缺點是很難設置一個適當?shù)耐队跋蛄糠稊?shù)閾值來分離可估參數(shù)和不可估參數(shù)。Machado 等[10]基于費雪信息矩陣(Fisher information matrix, FIM)準則進行參數(shù)子集選擇，通過枚舉法評估每個參數(shù)子集，此外還有行列式準則和一個修正的E 準則之間的比率 (a ratio between a determinant criterion and a modified E criterion，RDE)準則，但是隨著參數(shù)數(shù)量上升，將會面臨組合爆炸的問題。Brun等[11]提出了相關(guān)性和共線性方法選擇參數(shù)子集，通過靈敏度矩陣S的列中的線性關(guān)系來確定不可估計的參數(shù)。Jolliffe 等[12]提出主成分分析法(principal component analysis，PCA)，利用STS或者FIM 的特征值來選擇參數(shù)子集，缺點在于很難將特征值與特定的參數(shù)相匹配。此外Chen 等[13]提出用簡約海森矩陣(reduced Hessian,Hr)代替FIM，并且證明了等效性，避免了靈敏度矩陣S的計算，并基于系統(tǒng)輸出對參數(shù)進行排序，同時引入了關(guān)于待估參數(shù)的標準差與估計值比率(ratio of standard deviation to estimated value，RSDEV)準則用于可估和不可估參數(shù)的劃分。

基于模型的實驗設計(design of experiments，DOE)的目標是設計最佳的系統(tǒng)輸入，并以經(jīng)濟的方式使收集到的輸出中包含盡可能多的系統(tǒng)信息以便于后期的參數(shù)估計，采樣時間的選擇是DOE 的一個重要組成部分，特別是在實驗和樣本分析成本高昂的情況下。有很多文獻的最優(yōu)采樣使用了FIM 的A 準則、D準則，E 準則[14～16]，但是在參數(shù)之間存在高相關(guān)性的情況下，A 準則是不可靠的，而D 準則往往過于重視敏感參數(shù)[17]。Kutalik 等[18]提出了一種次優(yōu)采樣時間選擇方法，該方法是假設在選擇采樣時間之前，通過一些時間點進行參數(shù)估計，基于多重打靶法最大化FIM 的行列式來選擇采樣時間。結(jié)合基于擴展靈敏度狀態(tài)方程系統(tǒng)的優(yōu)化采樣設計和有限元的變化，Hoang 等[19]提出了一種先進的采樣策略，優(yōu)點是直接計算控制變量的參數(shù)靈敏度的精確高階導數(shù)，并且使用聯(lián)立求解。Chen 等[20]基于RSDEV 準則，提出一種迭代增強的近似最優(yōu)的時間序列采樣策略。

HDPE 模型十分復雜，使用自適應步長計算靈敏度矩陣S[21]時，無法使用同一步長計算多個點的靈敏度信息，這使得靈敏度矩陣S的計算成本過于巨大，導致無法使用與S相關(guān)的方法去選擇參數(shù)子集，同時也無法對模型中參數(shù)排序。因此本研究采用了Chen 等[13]提出的Hr矩陣代替FIM，并且結(jié)合RSDEV準則選出可估計參數(shù)子集，繼而提出了一種基于拉丁超立方體采樣的最優(yōu)采樣策略，結(jié)合RSDEV 準則與Hr的關(guān)系建立了混合整數(shù)非線性規(guī)劃(mixed integer nonlinear programming，MINLP)模型，并使用了一種新的迭代增強近似最優(yōu)采樣策略。數(shù)值結(jié)果表明，在RSDEV 準則約束條件下，可以顯著降低采樣個數(shù)，同時得到可以接受的待估參數(shù)的相對誤差。

2 高密度聚乙烯生產(chǎn)過程模型

HDPE 由于生產(chǎn)成本低、機械和流變性能廣泛，是應用最廣泛的合成商品聚合物之一[22]。多活性位Ziegler-Natta 催化下的HDPE 聚合反應是HDPE 生產(chǎn)的重要工藝，聚合物在催化劑的活性位點生長，直到鏈轉(zhuǎn)移發(fā)生，從而形成死聚合物鏈[23]。

考慮模型基于在連續(xù)攪拌釜反應器(continuous stirred tank reactor，CSTR)中進行的HDPE 漿工藝，由乙烯(單體)、氫氣(鏈轉(zhuǎn)移劑)、正己烷(稀釋劑)和Ziegler-Natta 催化劑系統(tǒng)構(gòu)成，其中四氯化鈦(TiCl4)作為催化劑，三乙基鋁(Al(C2H5)3)作為助催化劑，連續(xù)注入CSTR，并且將漿液產(chǎn)物連續(xù)從反應器中去除。假設各個活性位吸附單體的能力相同，則鏈增長、鏈轉(zhuǎn)移和鏈失活的動力學常數(shù)與鏈的長度無關(guān)[24]。Ziegler-Natta 催化劑的不同活性位點具有不同的傳播速率常數(shù)，會產(chǎn)生很多不同鏈長的聚合物，從而導致一個非常寬廣的質(zhì)量分數(shù)的分布。建立一個由大量微分和代數(shù)方程(differential-algebraic equation，DAE)組成的聯(lián)立方程模型，以描述該過程的動力學和質(zhì)量分數(shù)的分布。

從過程控制的角度來看，可用的操縱變量為催化劑、助催化劑、單體、氫、六烷的供給量和反應器溫度，而被控變量為單體轉(zhuǎn)化率和質(zhì)量分數(shù)的分布。

2.1 乙烯聚合反應機理

乙烯聚合反應采用的催化劑是Ziegler-Natta 催化劑，是一種復合催化劑。乙烯聚合反應機理包括催化劑活化、鏈引發(fā)反應、鏈增長反應、鏈轉(zhuǎn)移反應以及鏈失活反應[25]。

1) 催化劑活化。通過Al(C2H5)3對鈦進行烷基化形成活性位點，反應如下：

式中：CP(j)為具有潛在活性位點的催化劑TiCl4；A 為助催化劑Al(C2H5)3；kaA(j)為催化劑活化反應的動力學常數(shù)，L·(mol·s)-1；j為活性位編號；P0(j)為具有還未聚合任何單位的空活性位點的催化劑。

2) 鏈引發(fā)反應。單體分子與活性中心相互反應形成鏈長為1 的聚乙烯鏈，反應如下：

式中：M 為單體；ki(j)為鏈引發(fā)反應的動力學常數(shù)，L·(mol·s)-1；P1(j)為鏈長為1 且仍具有活性位點的乙烯聚合物鏈。

3) 鏈增長反應。繼續(xù)插入乙烯單元使鏈長不斷增長，反應如下：

式中：kp(j)為鏈增長反應的動力學常數(shù)，L·(mol·s)-1；Pn(j)為鏈長為n的活聚合物鏈。

4) 鏈轉(zhuǎn)移反應?；罹酆衔镦湻謩e與單體、氫氣、助催化劑反應及自反應形成死鏈，反應如下：

式中：ktM(j）為向單體轉(zhuǎn)移鏈反應的動力學常數(shù)，L·(mol·s)-1；ktH(j）為向氫氣轉(zhuǎn)移鏈反應的動力學常數(shù)，L0.5·(mol0.5·s-1)；ktA(j）為向助催化劑轉(zhuǎn)移鏈反應的動力學常數(shù)，L·(mol·s)-1；kt(j）為自轉(zhuǎn)移反應的動力學常數(shù)，s-1；Dn(j)為鏈長為n的死聚合物鏈。

5) 鏈失活反應?；罹酆衔镦溑c活性中心全部失活，反應如下：

式中：kd(j）為自失活反應的動力學常數(shù)，s-1；Cd(j)為具有不能再聚合任何單體的死活性位點的催化劑。

在這個模型中，活性位點的數(shù)量被設為2，根據(jù)乙烯聚合反應機理，轉(zhuǎn)移終止的偽動力學速率系數(shù)定義為

式中：c[·]為括號中某成分的濃度，mol·L-1。

2.2 物料衡算

為了簡潔起見，模型中將在第j個活性位點的活和死聚合物的m階矩分別定義為

活聚合物鏈的零階矩表示為

對于長度為n=1 的活鏈，反應凈生產(chǎn)率為

對于長度n≥2 的活鏈來說，反應凈生產(chǎn)率為

同樣地，長度為n=1 和n≥2 的死鏈的反應凈生產(chǎn)率分別為

從式(14)～(17)可得，長度n=1,2,3,…,nmax的活鏈和死鏈的粒數(shù)衡算方程分別為

式中：Qout為體積流出率，%；V為反應堆的液體體積，L。類似的單體、氫氣、己烷、催化劑潛在活性位點、催化劑無活性位點、活性點、共催化劑以及活聚合物鏈的零階矩和一階矩的物料衡算方程也可以寫為

式中：Ns為活性位點個數(shù)，Ns=2；Fin(M)、Fin(H2)、Fin(C6H14)、Fin(CP)、Fin(A)分別為M、H2、C6H14、CP、A 的摩爾進料流率，mol·h-1。

Ma 等[26]提出一種利用二維正交配置的質(zhì)量分數(shù)分布重構(gòu)方法，以捕捉質(zhì)量分數(shù)分布的時間動態(tài)特征和鏈長度的分布特征。采用聯(lián)立配置法離散化模型，將模型最終離散成具有86 116 個變量，86 100 條等式約束的大型非線性規(guī)劃(nonlinear programming，NLP)問題。16 個待估參數(shù)分別為kp(1)、kp(2)、ktH(1)、ktH(2)、ktM(1)、ktM(2)、ktA(1)、ktA(2)、kaA(1)、kaA(2)、ki(1)、ki(2)、kt(1)、kt(2)、kd(1)、kd(2)。模型采用以下目標函數(shù)進行參數(shù)估計：

式中：2σ為測量噪聲的方差，tf為時間采樣點，wcal為模型質(zhì)量分數(shù)的分布仿真值，wrel為模型質(zhì)量分數(shù)的分布實際值。

3 最優(yōu)采樣

由于HDPE 采樣后的分析成本過高，故盡可能減少采樣數(shù)量，同時盡量獲得多的信息，本研究將結(jié)合RSDEV 準則和FIM 的A 準則來選擇參數(shù)子集，然后基于拉丁超立方體采樣，使用Chen 等[20]提出的一種基于Hr的近似最優(yōu)的時間序列采樣策略。

3.1 參數(shù)子集選擇

在參數(shù)子集選擇方法上，Chen等[13]證明了Hr與FIM的等效性，故采用Hr代替FIM，并且將基于FIM準則選擇參數(shù)子集。常用的FIM準則包括了A準則、D準則、修正后的E準則，此外還有RDE準則，如果參數(shù)數(shù)量太多，則會面臨組合爆炸問題，為避免此問題，采用Chen等[13]提出的RSDEV準則，即

式中：δq為第q個待估參數(shù)標準差；qθ為第q個待估參數(shù)的估計值；ε為比率；δq=aq,q1/2，aq,q為相應參數(shù)的協(xié)方差矩陣A的對角元素的值，即A矩陣的第q個對角元素，A=FIM-1=Hr-1。

將FIM 的A 準則和RSDEV 準則結(jié)合在一起形成新的準則：

min (tr (FIM-1))，且滿足max (ε)<εexp。式中：εexp為期望的比率閾值。

3.2 最優(yōu)采樣策略

將HDPE 模型表達成微分代數(shù)方程組，即

式中：K為插值的階次，也是配置點的數(shù)量，xik為在第i個有限元的第k個配置點的狀態(tài)變量向量，yik為在第i個有限元的第k個配置點的代數(shù)變量向量，τ為有限元中t的相對值，kτ、jτ分別為有限元中第k、j個配置點t的相對值。在第l個采樣時間點tl處的預測可以用以下方程來近似：

式中：zl,cal為采樣時間點tl處的模擬預測向量。

假設一個{1,…,ntp}的時間序列，它可以作為潛在的采樣點，相應在tl的測量是zl，l=1, …,ntp。基于似然原理，假設z的每個分量的測量噪聲呈正態(tài)分布，均值為0、方差為σp2，p=1,…,nz，通過求解問題(35)可以估計θ：

式中：zl,rel為采樣時間點tl處的實際預測向量，Σ為對角線方差矩陣，其中第p個元素為σp2；l j′(kτ)為導數(shù)。

設S={t1,t2,…,tn}是采樣點集合，是時間序列{1,2,…,ntp}的一個子集。為了便于最優(yōu)抽樣分析，將S={1,2,…,ntp}，并且將非線性規(guī)劃問題(35)中的離散變量和方程被重寫為

式中：I為單位矩陣。

確定δq/θq的值，即參數(shù)q的標準差除以其估計值，建立一個MINLP 問題尋找輸出測量時間序列的最小子集，從而可以找到一個滿足上述比率的A矩陣。

首先定義了Zr，它由Z的最后一個nz×ntp行組成：

用二元變量b=[b1,b2,…,bntp]T來確定哪個采樣時刻tl(l=1,2,…,ntp)從時間序列中被選擇。因此，關(guān)于時間序列S的Hr可以定義為

在當前參數(shù)θ下的協(xié)方差矩陣A=(Hr(b))-1，通過求解以下MINLP 問題來確定b，時間序列子集和合適的A矩陣元素am,s：

然而，使用問題式(40)的前提是，真實的估計參數(shù)值是先驗已知的，但事實并非如此。相反，該策略通過在參數(shù)空間上進行隨機采樣來選擇θ的候選值，并用每個候選值來解決問題(40)。從每個解中，記錄從b中選擇了哪些時間序列點，并觀察從整個采樣中選擇時間序列點的頻率，然后選擇大頻率的時間序列點進行參數(shù)估計。

并且Chen 等[20]證明了輸出采樣數(shù)量的增加會使得RSDEV 準則ε的值減小或者保持不變。因此，MINLP(40)的解決方案提供了一個最小的時間序列點集，并且同時保證了一個適當?shù)谋嚷师拧?/p>

3.3 拉丁超立方體采樣

采樣策略多次在參數(shù)空間中隨機采樣選擇θ的候選值，并用每個候選值來解決問題式(40)，然而隨機采樣使結(jié)果充滿不確定性。建立精確的代理模型要求采樣點集可以描述模型在輸入范圍上的整體趨勢變化，因而需要采樣點足夠均勻，數(shù)量也足夠代表模型的變化，為了達到這一點，本研究采用了拉丁超立方體采樣(Latin hypercube sampling，LHS)取代隨機采樣。LHS 是分層采樣抽樣技術(shù)中的一種，由Mckay等[27]于1979 年提出。分層抽樣法將總體按一定特征或者規(guī)則劃分為不同的層，然后從不同的層內(nèi)進行獨立、隨機的抽樣，相比于簡單隨機抽樣，分層采樣獲得的樣本代表性較好，抽樣誤差較小。

在LHS 方法中在每個維度上都進行了均勻的劃分，但是采樣點的數(shù)目不受限制。假設采樣空間的維度為n，想要獲得m個采樣點，那么要將n個維度中的每個維度上的取值空間劃分為m個等份，然后在每個維度的每個等分內(nèi)都隨機產(chǎn)生一個數(shù)值，由此就獲得了n組數(shù)量為m的數(shù)值集合，之后將這n組數(shù)值集合內(nèi)的點隨機組合，就獲得了所需要的m個采樣點[28]。

在此策略中，維度由參數(shù)子集選擇技術(shù)決定，設可估參數(shù)子集個數(shù)為nq，采樣次數(shù)為nrs，將參數(shù)子集的參數(shù)空間分成nrs 個等份，在每個等份中隨機取一個值，將nq 組數(shù)值集合內(nèi)的點隨機組合，獲得輸入采樣。

3.4 算法步驟

步驟1：通過模型式(32)確定可估計參數(shù)子集，可估參數(shù)子集個數(shù)為nq。

步驟2：設定式(40)中的εexp，乘子η>1，采樣次數(shù)為nrs。

步驟3：設定可估計參數(shù)的參數(shù)空間，用拉丁超立方體采樣在參數(shù)空間里采樣nrs 次生成候選參數(shù)，并解決nrs 次候選參數(shù)的問題式(40)，然后根據(jù)所有采樣的參數(shù)值計算時間序列點的頻率。

步驟4：設最終采樣集合S由nq×η個時間序列點組成，選擇頻率最高的nq×η個時間點組成S。

步驟5：用集合S中的采樣點時刻的測量值來求解問題(36)，即參數(shù)估計；計算參數(shù)協(xié)方差矩陣A，并表示最大比率δq/θq，q=1,…,nq。

步驟6：如果δq/θq<εexp，轉(zhuǎn)至步驟9。

步驟7：如果δq/θq≥εexp，設Sl={1, 2,…, ntp}-S，對Sl中的每個元素{v}，令St=S+{v}，并計算其Hr和協(xié)方差矩陣A，并保留其δq*/θq*的值。

步驟8：選擇δq*/θq*的最大值對應的元素vmax，令S=S+{vmax}，返回步驟5。

步驟9：結(jié)束。

算法流程圖如圖(1)所示：

圖1 最優(yōu)策略算法流程圖Fig.1 Flow chart of the optimal strategy algorithm

4 仿真實驗

以質(zhì)量分數(shù)的分布作為測量輸出，通過加入方差1×10-4來模擬測量值，有限元個數(shù)為30 個，每個有限元配置點K=3，共有90 個潛在采樣點。根據(jù)第3 節(jié)的參數(shù)子集選擇方法，設定εexp=0.05，被估計參數(shù)子集θ=[kp(2),ktH(1),kaA(2)]，nq=3，其參數(shù)空間的范圍分別為3 554.041≤kp(2)≤3 871.933，L·(mol·s)-1；182.361≤ktH(1)≤216.283，L0.5·mol-0.5·s-1；204.329≤kaA(2) ≤242.987，L·(mol·s)-1；設定式(40)中的εexp=0.05，乘子η=2，采樣次數(shù)nrs 為100。在最優(yōu)采樣策略的步驟3 中，使用隨機采樣，得出所有時間點的頻率，如圖2 所示。

圖2 采用隨機采樣策略得到的時間序列點選擇頻率Fig.2 Frequency of time series points for estimating the parameter range of HDPE using random sampling

根據(jù)頻率最高選擇采樣時間點S={9, 15, 43, 69, 77, 90}，采樣時間點數(shù)量為6。使用選擇的采樣時間點估計參數(shù)θ=[kp(2),ktH(1),kaA(2)]，其值分別為[3 423.207, 221.58, 244.6]，其δq/θq分別為[0.038 8, 0.045 5,0.048 4]，誤差分別為6.23%、6.27%、7.23%。

用拉丁超立方體采樣代替步驟3中的隨機采樣，同樣得到所有采樣的參數(shù)值計算時間序列點的頻率，如圖3所示。

圖3 采用拉丁超立方體采樣策略得到的時間序列點選擇頻率Fig.3 Frequency of time series points for estimating the parameter range of HDPE using Latin hypercube sampling

根據(jù)頻率最高選擇采樣時間點S={9, 12, 15, 43, 77, 90}，采樣時間點數(shù)量為6。使用選擇的采樣時間點估計參數(shù)θ=[kp(2),ktH(1),kaA(2)]，其值分別為[3 542.722, 213.84, 237.2]，其δq/θq分別為[0.033 3, 0.036 9,0.040 5]，誤差分別為2.96%、2.56%、3.97%。

為對比在步驟3中使用拉丁超立方體采樣與使用隨機采樣對最終參數(shù)估計的影響，分別使用這2種采樣方法仿真了10次，并且取其誤差平均值作為對比。使用隨機采樣的參數(shù)估計的誤差平均值為5.82%、6.12%、6.99%，其方差分別為0.000 325，0.000 436，0.000 387，使用拉丁超立方體采樣的參數(shù)估計的誤差平均值為3.63%、3.49%、4.11%，其方差分別為0.000 377、0.000 631、0.000 563。

5 結(jié) 語

針對HDPE 模型的樣本分析成本過高且參數(shù)估計困難問題，提出使用簡約Hessian 矩陣代替FIM 并結(jié)合RSDEV 準則作為新的參數(shù)子集選擇方法，其優(yōu)勢在于避免了HDPE 模型靈敏度矩陣S計算成本過大和組合爆炸問題。此外，在保證參數(shù)估計精度的情況下，為了減少采樣點個數(shù)，降低樣本分析成本，提出一種基于拉丁超立方體采樣的近似最優(yōu)采樣策略。通過拉丁超立方體采樣獲得參數(shù)空間候選值，并建立和求解最優(yōu)采樣子問題，構(gòu)建了各潛在采樣點的選擇頻率，使用一種啟發(fā)式策略和迭代增強策略選擇最終采樣點，并實現(xiàn)了參數(shù)估計。為了驗證采樣策略結(jié)果的可靠性，使用了基于隨機采樣的最優(yōu)采樣策略進行對比。使用新的參數(shù)子集選擇方法選出的待估參數(shù)為kp(2)、ktH(1)、kaA(2)，同時仿真結(jié)果表明，最終在90 個潛在采樣點中僅采樣其中6 個，大大減少了采樣點個數(shù)，降低了樣本分析成本，使用隨機采樣的參數(shù)估計相對誤差平均值分別為5.82%、6.12%、6.99%，使用拉丁超立方體采樣的參數(shù)估計的相對誤差平均值分別為3.63%、3.49%、4.11%。從仿真結(jié)果可以看出參數(shù)估計精度在可接受范圍，而且基于拉丁超立方體采樣比基于隨機采樣的參數(shù)估計結(jié)果更加精確、穩(wěn)定。此外，該最優(yōu)采樣策略還適用于不能在線測量和對樣本離線分析成本很高的情況。