均勻?qū)嶒?yàn)設(shè)計(jì)在化學(xué)反應(yīng)動(dòng)力學(xué)建模中的應(yīng)用研究I.非線性代數(shù)動(dòng)力學(xué)模型

胡葛林，呂 慧，曹貴平

聯(lián)合化學(xué)反應(yīng)工程研究所華東理工大學(xué)分所，華東理工大學(xué)化工學(xué)院，上海 200237

化學(xué)反應(yīng)動(dòng)力學(xué)是化學(xué)工程重要內(nèi)容之一。根據(jù)反應(yīng)動(dòng)力學(xué)不僅可以解析反應(yīng)過(guò)程機(jī)理、探究反應(yīng)路徑，還可以進(jìn)行工業(yè)化反應(yīng)器選型、設(shè)計(jì)、放大、優(yōu)化、控制和預(yù)測(cè)。建立化學(xué)反應(yīng)動(dòng)力學(xué)模型的通常方法是：固定反應(yīng)物初始濃度，測(cè)定不同反應(yīng)溫度、反應(yīng)時(shí)間下反應(yīng)物和產(chǎn)物的濃度，得到濃度和時(shí)間的關(guān)系數(shù)據(jù)；再利用非線性模型參數(shù)估計(jì)方法，迭代求解得到各溫度水平下的速率常數(shù)和反應(yīng)級(jí)數(shù)；進(jìn)一步利用Arrhenius方程對(duì)溫度進(jìn)行擬合，得到頻率因子和活化能。這種求取模型參數(shù)的方法需要開展大量的實(shí)驗(yàn)工作，存在工作量大、取樣困難等問(wèn)題。

Fang等[1-2]基于擬蒙特卡洛方法、數(shù)論方法和多變量統(tǒng)計(jì)分析提出了均勻?qū)嶒?yàn)設(shè)計(jì)思想，即在實(shí)驗(yàn)范圍內(nèi)選取一系列均勻分布的代表性點(diǎn)，以較少的實(shí)驗(yàn)次數(shù)解決多因素、多水平工藝優(yōu)化和配方設(shè)計(jì)問(wèn)題。

相較于全面實(shí)驗(yàn)、正交實(shí)驗(yàn)等傳統(tǒng)實(shí)驗(yàn)設(shè)計(jì)方法，均勻?qū)嶒?yàn)設(shè)計(jì)的優(yōu)勢(shì)在于：（1）能夠在實(shí)驗(yàn)域內(nèi)選取一系列具有高度代表性且在空間內(nèi)均勻分布的實(shí)驗(yàn)點(diǎn)，可近似替代全面實(shí)驗(yàn)設(shè)計(jì)；（2）不需要對(duì)模型有較深的認(rèn)識(shí)，不對(duì)模型施加較強(qiáng)的假設(shè)；（3）可在合理的實(shí)驗(yàn)次數(shù)內(nèi)為因素安排較多的實(shí)驗(yàn)水平，增強(qiáng)因素的感知度。因此，均勻?qū)嶒?yàn)設(shè)計(jì)在紡織、冶金、醫(yī)藥、化學(xué)化工和農(nóng)業(yè)等領(lǐng)域得到了廣泛應(yīng)用[3-4]，其相關(guān)理論和應(yīng)用還在不斷發(fā)展中。Xin等[5]基于均勻?qū)嶒?yàn)設(shè)計(jì)進(jìn)行瀝青粘合劑配方優(yōu)化，研究了硬質(zhì)改性瀝青性能和其組分配比之間的關(guān)系，以提升瀝青高溫性能。Li等[6]采用均勻?qū)嶒?yàn)設(shè)計(jì)系統(tǒng)研究了等離子體噴涂過(guò)程中電弧電流、Ar流率、H2流率、噴涂距離和粉末進(jìn)料速率5個(gè)工藝參數(shù)對(duì)TiN涂層沉積效率、孔隙率、氧化物含量、顯微硬度和斷裂韌性的影響。此外，均勻設(shè)計(jì)能夠在空間中生成均勻分布實(shí)驗(yàn)點(diǎn)的特點(diǎn)也引起了人們的極大興趣[7]。在機(jī)器學(xué)習(xí)中，由于“更均勻”、“空間填充更多”的優(yōu)勢(shì)，均勻設(shè)計(jì)被用于模型選擇，以提升支持向量機(jī)的一般性能，均勻設(shè)計(jì)的運(yùn)用不僅大大減小了參數(shù)實(shí)驗(yàn)的次數(shù)，也提供了在計(jì)算時(shí)間限制下調(diào)整候選集大小的靈活性[8]。根據(jù)均勻?qū)嶒?yàn)設(shè)計(jì)的基本原理，也可將均勻?qū)嶒?yàn)設(shè)計(jì)用于高效求取反應(yīng)動(dòng)力學(xué)模型，可望大幅減少動(dòng)力學(xué)實(shí)驗(yàn)的次數(shù)?；诖耍竟ぷ魇状翁岢隼镁鶆?qū)嶒?yàn)設(shè)計(jì)原理進(jìn)行反應(yīng)動(dòng)力學(xué)模型建立的實(shí)驗(yàn)設(shè)計(jì)，將反應(yīng)溫度和反應(yīng)時(shí)間在其變化范圍內(nèi)均勻布置，每個(gè)溫度下僅在一個(gè)時(shí)間點(diǎn)取樣分析得到樣品濃度，將參數(shù)對(duì)溫度的依賴關(guān)系代入反應(yīng)動(dòng)力學(xué)模型，用非線性模型參數(shù)估計(jì)方法直接得到頻率因子和反應(yīng)活化能，對(duì)參數(shù)估計(jì)結(jié)果作統(tǒng)計(jì)學(xué)分析評(píng)估，較為系統(tǒng)地研究影響模型參數(shù)準(zhǔn)度和精度的因素。

1 化學(xué)反應(yīng)動(dòng)力學(xué)模型

式中：Ci為物質(zhì)i的濃度，mol/L；t為反應(yīng)時(shí)間，min；k為反應(yīng)速率常數(shù)，min-1。反應(yīng)初始條件為：t=0時(shí)，CA0=1 mol/L，CB0=CC0=0 mol/L。根據(jù)Arrhenius方程，k與溫度的關(guān)系為：

式中：k0為頻率因子，min-1；E為活化能，kJ/mol；R為氣體常數(shù)，J/(mol·K)；T為溫度，K。將式（2）代入式（1），得到任意溫度下動(dòng)力學(xué)模型：

式（3）所示的動(dòng)力學(xué)模型中待求參數(shù)有k10，E1，k20和E2。

控制反應(yīng)物A的轉(zhuǎn)化率在95%以內(nèi)，基于前期探索，選擇反應(yīng)溫度因素的區(qū)間為135～160 ℃，水平數(shù)為6個(gè)，步長(zhǎng)為5 ℃；反應(yīng)時(shí)間因素的區(qū)間為14～44 min，水平數(shù)為6個(gè)，步長(zhǎng)為6 min，溫度和時(shí)間的水平見表1。選用U6*(64)規(guī)劃實(shí)驗(yàn)[10]，如表2所示，按表2中第1列和第3列安排實(shí)驗(yàn)。

表1 因素和水平Table 1 Factors and levels

表2 均勻表U6＊(64)Table 2 Uniform Table U6＊(64)

2 模擬實(shí)驗(yàn)數(shù)據(jù)的獲取

采用隨機(jī)數(shù)函數(shù)方法獲得不同溫度和不同時(shí)間下物質(zhì)i的濃度實(shí)驗(yàn)值誤差，進(jìn)而得到實(shí)驗(yàn)值，如式（4）所示：

式中：Cij,exp表示第j組實(shí)驗(yàn)中物質(zhì)i的濃度實(shí)驗(yàn)值，mol/L；Cij,mod表示第j組中物質(zhì)i的濃度模型值，mol/L，通過(guò)式（3）計(jì)算得到；Tj表示第j組實(shí)驗(yàn)的反應(yīng)溫度，℃；tj表示第j組實(shí)驗(yàn)的反應(yīng)時(shí)間，min；k*表示模型參數(shù)真值；eij表示第j組實(shí)驗(yàn)中物質(zhì)i濃度的隨機(jī)誤差，eij服從數(shù)學(xué)期望為0、方差為2σ的正態(tài)分布；M表示實(shí)驗(yàn)組數(shù)，即均勻?qū)嶒?yàn)設(shè)計(jì)中溫度和時(shí)間的水平數(shù)。

不同溫度、時(shí)間下A，B和C濃度的模型值CA,mod，CB,mod和CC,mod是通過(guò)將模型參數(shù)的真值、溫度和時(shí)間代入式（3）計(jì)算得到，為不失一般性，這里模型參數(shù)真值按照文獻(xiàn)[9]方法計(jì)算給出：k10=6.856×1013min-1，E1=1.233 1×102kJ/mol，k20=1.500×1016min-1，E2=1.446 3×102kJ/mol。根據(jù)式（4）計(jì)算得到模型值Cij,mod，同時(shí)加上誤差隨機(jī)數(shù)eij得到模擬實(shí)驗(yàn)數(shù)據(jù)值CA,exp，CB,exp和CC,exp，隨機(jī)誤差eji由正態(tài)隨機(jī)數(shù)函數(shù)計(jì)算產(chǎn)生。控制隨機(jī)誤差的標(biāo)準(zhǔn)差（σ）為0.001，得到的實(shí)驗(yàn)數(shù)據(jù)如表3所示，其中也列出關(guān)鍵反應(yīng)物A的轉(zhuǎn)化率（xA）。

表3 U6＊(64)實(shí)驗(yàn)設(shè)計(jì)及實(shí)驗(yàn)數(shù)據(jù)Table 3 U6＊(64) uniform design scheme and experimental data

由于各濃度的權(quán)值未知，參數(shù)估計(jì)時(shí)采用最大似然估計(jì)（MLE）目標(biāo)函數(shù)[9]，如式（5），迭代算法為改進(jìn)單純形法（ISA）。

3 結(jié)果與討論

3.1 模型顯著性檢驗(yàn)

利用式（5）目標(biāo)函數(shù)，采用改進(jìn)單純形方法回歸得到的模型參數(shù)估計(jì)值：k10=7.092 3×1013min-1，k20=1.565 2×1016min-1，E1=1.234 3×102kJ/mol，E2=1.447 8×102kJ/mol。相關(guān)系數(shù)R2=0.999 9，表明99.99%的觀測(cè)值可以用此模型來(lái)解釋[11]。將模型參數(shù)估計(jì)值代入反應(yīng)動(dòng)力學(xué)模型，在反應(yīng)溫度和時(shí)間范圍內(nèi)，濃度計(jì)算值以空間網(wǎng)格曲面表示，按均勻?qū)嶒?yàn)設(shè)計(jì)表U6*(64)安排的溫度和時(shí)間下得到的CA,exp，CB,exp和CC,exp以圓點(diǎn)表示，如圖1～圖3所示。由圖1～圖3可知，實(shí)驗(yàn)點(diǎn)均勻地分布在實(shí)驗(yàn)范圍內(nèi)，與模型曲面有緊密的貼合，吻合很好，表明參數(shù)估計(jì)效果較好[12]。方差分析如表4所示。

圖1 CA隨溫度和時(shí)間的變化Fig.1 Surface plots of concentrations as a function of temperatures and time

圖2 CB隨溫度和時(shí)間的變化Fig.2 Surface plots of concentrations as a function of temperatures and time

圖3 CC隨溫度和時(shí)間的變化Fig.3 Surface plots of concentrations as a function of temperatures and time

表4 方差分析Table 4 Analysis of variance

表4中，DF表示自由度，F(xiàn)-value為回歸均方和與殘差均方和的比值，P-value為假設(shè)概率。對(duì)于非線性模型，當(dāng)F-value大于10倍Fα，可以認(rèn)為模型在α置信水平下是顯著的[13]。由表4可知，反應(yīng)動(dòng)力學(xué)模型的F-value=7.14×105F0.01(4,14)=50.4，表明所得到的動(dòng)力學(xué)模型具有99%的置信度。P-value=5.45×10-370.01，表明模型在統(tǒng)計(jì)學(xué)上是高度顯著的[14]。

為了驗(yàn)證和方便比較模型參數(shù)估計(jì)效果，計(jì)算了參數(shù)估計(jì)值與真實(shí)值的相對(duì)誤差（ε），結(jié)果如表5所示。由表5可知，各模型參數(shù)估計(jì)值與真值的相對(duì)誤差均在5%以內(nèi)?？梢钥闯觯镁鶆?qū)嶒?yàn)設(shè)計(jì)將溫度和時(shí)間因素在變化范圍內(nèi)均勻分布，每個(gè)溫度下僅在一個(gè)時(shí)間點(diǎn)取樣分析，僅安排6次實(shí)驗(yàn)就能獲得反應(yīng)動(dòng)力學(xué)模型參數(shù)。通常，為了建立動(dòng)力學(xué)模型，對(duì)于溫度和時(shí)間兩個(gè)因素的實(shí)驗(yàn)，溫度需要安排至少5水平，時(shí)間需要安排至少10水平，至少需要50次實(shí)驗(yàn)，實(shí)驗(yàn)周期很長(zhǎng)，特別是在高溫、高壓等苛刻條件下在線取樣困難，而且多次在線取樣會(huì)造成間隙反應(yīng)器中的物料損失，影響后續(xù)反應(yīng)動(dòng)力學(xué)行為，帶來(lái)不可低估的實(shí)驗(yàn)誤差，影響動(dòng)力學(xué)模型的準(zhǔn)確性。而通過(guò)均勻?qū)嶒?yàn)設(shè)計(jì)，實(shí)驗(yàn)次數(shù)減少到6次，實(shí)驗(yàn)工作量大幅減小，用少量的實(shí)驗(yàn)即可獲得實(shí)驗(yàn)數(shù)據(jù)，因?qū)嶒?yàn)數(shù)據(jù)空間的均勻性，具有動(dòng)力學(xué)規(guī)律的空間代表性，因而得到的模型參數(shù)值具有很好的準(zhǔn)確性和精確性，表明均勻?qū)嶒?yàn)設(shè)計(jì)在多因素、多水平動(dòng)力學(xué)模型建立中有著顯著優(yōu)勢(shì)[15]。

表5 參數(shù)估計(jì)值及其與真實(shí)值的相對(duì)誤差Table 5 Parameter estimation values and relative errors between estimated values and real values

3.2 實(shí)驗(yàn)隨機(jī)誤差的影響

在3.1節(jié)中，控制A，B和C濃度的實(shí)驗(yàn)值與模型值之間的隨機(jī)誤差標(biāo)準(zhǔn)差（σ）為0.001，但在實(shí)際問(wèn)題中的隨機(jī)誤差可能會(huì)更大。保持溫度、時(shí)間變化范圍和水平、參數(shù)初值不變，探究了實(shí)驗(yàn)隨機(jī)誤差對(duì)參數(shù)估計(jì)精度的影響，結(jié)果見表6。由表6可知，實(shí)驗(yàn)誤差標(biāo)準(zhǔn)差為0.001～0.1時(shí)，對(duì)參數(shù)估計(jì)結(jié)果不會(huì)產(chǎn)生影響；若進(jìn)一步增大σ至0.25或0.5，出現(xiàn)濃度計(jì)算值為負(fù)數(shù)或大于初始值的異常結(jié)果?？梢姦摇?.1可以得到滿意的結(jié)果。

表6 不同實(shí)驗(yàn)隨機(jī)誤差標(biāo)準(zhǔn)差下的參數(shù)估計(jì)值和真值的相對(duì)誤差Table 6 The relative errors between estimated values and real values under different σ

3.3 因素水平的影響

由于均勻?qū)嶒?yàn)設(shè)計(jì)所需的實(shí)驗(yàn)組數(shù)等于因素的水平數(shù)，可以改變因素在變化范圍內(nèi)的水平數(shù)來(lái)增大或減少實(shí)驗(yàn)次數(shù)，以滿足不同實(shí)驗(yàn)條件要求和參數(shù)估計(jì)準(zhǔn)度和精度的需要[16]。表7匯總了不同水平數(shù)對(duì)參數(shù)估計(jì)值的影響。由表7可知，在各個(gè)水平數(shù)下，4個(gè)模型參數(shù)估計(jì)值和真值的相對(duì)誤差均可控制在6%以內(nèi)，E1和E2的參數(shù)估計(jì)精度明顯優(yōu)于k10和k20，E1和E2的參數(shù)估計(jì)值與真值的相對(duì)誤差保持在1%以內(nèi)?？傮w來(lái)說(shuō)，不同因數(shù)水平下的參數(shù)估計(jì)結(jié)果均令人滿意。

表7 不同因素水平下的參數(shù)估計(jì)值和真值的相對(duì)誤差Table 7 The relative errors between estimated values and real values under different levels

3.4 擬水平的影響

采用擬水平的方法，仍保持溫度的變化為135～160 ℃，將溫度的水平數(shù)減小到3個(gè)，分別為135.0，147.5和160.0 ℃；時(shí)間的變化范圍和水平數(shù)不變，即14～44 min，水平數(shù)為6個(gè)。將表2均勻表U6*(64)中第1列的水平數(shù)進(jìn)行合并：[1,2]→[1]，[3,4]→[2]，[5,6]→[3]，得到擬水平的均勻?qū)嶒?yàn)設(shè)計(jì)表U6(3×6)，如表8所示。參數(shù)估計(jì)結(jié)果和估計(jì)值與真實(shí)值的相對(duì)誤差如表9所示。

表8 U6(3×6)實(shí)驗(yàn)設(shè)計(jì)及實(shí)驗(yàn)數(shù)據(jù)Table 8 U6(3×6) uniform design and experimental data

表9 擬水平下參數(shù)值及其與真實(shí)值的相對(duì)誤差Table 9 Estimated values of parameters and relative errors between estimated values and real values under pseudo level

擬水平實(shí)驗(yàn)設(shè)計(jì)參數(shù)估計(jì)結(jié)果與均勻?qū)嶒?yàn)設(shè)計(jì)參數(shù)估計(jì)的比較結(jié)果見表10。由表10可知，擬水平實(shí)驗(yàn)設(shè)計(jì)參數(shù)估計(jì)結(jié)果與非擬水平的均勻?qū)嶒?yàn)設(shè)計(jì)無(wú)顯著差異，最大參數(shù)估計(jì)值與真值的相對(duì)誤差為k10時(shí)的3.47%，擬水平實(shí)驗(yàn)設(shè)計(jì)的應(yīng)用在保證參數(shù)估計(jì)精度的同時(shí)，將溫度水平數(shù)從6個(gè)減小到了3個(gè)，提高了實(shí)驗(yàn)效率。

表10 擬水平均勻?qū)嶒?yàn)設(shè)計(jì)與均勻?qū)嶒?yàn)設(shè)計(jì)的參數(shù)估計(jì)結(jié)果比較Table 10 Comparison about parameter estimation results between Pseudo-level and uniform level

4 結(jié) 論

根據(jù)均勻?qū)嶒?yàn)設(shè)計(jì)原理和方法，首次提出在溫度和時(shí)間二維空間內(nèi)均勻布局實(shí)驗(yàn)點(diǎn)，獲得建立化學(xué)反應(yīng)動(dòng)力學(xué)模型的“濃度-溫度-時(shí)間”數(shù)據(jù)，僅用6組實(shí)驗(yàn)數(shù)據(jù)，便可一步完成模型參數(shù)的估計(jì)，為高效建立動(dòng)力學(xué)模型提供一種新思路和新方法。參數(shù)估計(jì)值與真值的相對(duì)誤差可控制在5%以內(nèi)，統(tǒng)計(jì)分析表明模型參數(shù)估計(jì)結(jié)果具有良好的準(zhǔn)確性、可靠性和顯著性。不同實(shí)驗(yàn)誤差方差下得到的參數(shù)估計(jì)結(jié)果基本一致，可控的實(shí)驗(yàn)誤差對(duì)參數(shù)估計(jì)結(jié)果的影響可以忽略。不同水平數(shù)下，模型參數(shù)估值和真值的相對(duì)誤差均可控制在6%以內(nèi)，參數(shù)估計(jì)結(jié)果并不會(huì)因因素水平數(shù)的改變而發(fā)生明顯變化，可以用最少的實(shí)驗(yàn)獲得動(dòng)力學(xué)模型參數(shù)。擬水平均勻?qū)嶒?yàn)設(shè)計(jì)的運(yùn)用可以將溫度因素的水平數(shù)進(jìn)一步減小到3個(gè)，同時(shí)保證參數(shù)估計(jì)值有較高的精度。