覃溢波, 柳天虹, 任春雨, 齊勝利

(揚(yáng)州大學(xué)信息工程學(xué)院, 江蘇 揚(yáng)州 225127)

連續(xù)攪拌釜式反應(yīng)器(continuous stirred tank reactor, CSTR)是化工生產(chǎn)過(guò)程中進(jìn)行各種物理變化和化學(xué)反應(yīng)廣泛使用的設(shè)備, 在反應(yīng)裝置中占有重要地位．釜內(nèi)攪拌裝置的主要作用是使反應(yīng)釜內(nèi)物料混合均勻且物料溫度處處相等,以保證物料的質(zhì)量和反應(yīng)釜的安全性．近年來(lái), 關(guān)于CSTR的研究主要集中在反應(yīng)釜預(yù)測(cè)控制[1-2]及狀態(tài)估計(jì)[3]等方面．由于釜式反應(yīng)器中化學(xué)反應(yīng)機(jī)理較復(fù)雜,所采集的數(shù)據(jù)通常表現(xiàn)出非線性強(qiáng)、樣本維數(shù)高以及存在較大噪聲和孤立點(diǎn)等特性．目前,針對(duì)反應(yīng)釜輸出量的預(yù)測(cè)主要采用神經(jīng)網(wǎng)絡(luò)模型[4-5],該模型以經(jīng)驗(yàn)風(fēng)險(xiǎn)最小化為原則,模型建立在樣本數(shù)趨于無(wú)窮的基礎(chǔ)上, 但是實(shí)際訓(xùn)練樣本往往是有限的, 導(dǎo)致訓(xùn)練過(guò)程中模型過(guò)于復(fù)雜而產(chǎn)生過(guò)擬合現(xiàn)象．高斯過(guò)程回歸(Gaussian process regression, GPR)[6]因結(jié)構(gòu)靈活性高而成為主流的機(jī)器學(xué)習(xí)模型之一．GPR是一種非參數(shù)建模方法, 可表征為隨機(jī)變量的集合, 其中的任意隨機(jī)變量組合均服從聯(lián)合高斯分布．GPR也可等效為一種特殊結(jié)構(gòu)的神經(jīng)網(wǎng)絡(luò), 即隱含層含有無(wú)限節(jié)點(diǎn)的神經(jīng)網(wǎng)絡(luò)．相較于神經(jīng)網(wǎng)絡(luò)等方法, GPR因包含的參數(shù)較少, 故參數(shù)優(yōu)化過(guò)程簡(jiǎn)單且更易收斂,建模時(shí)可將先驗(yàn)知識(shí)表征為先驗(yàn)概率融入模型,通過(guò)選取不同類型的協(xié)方差函數(shù)增強(qiáng)模型對(duì)實(shí)際過(guò)程的解釋能力．GPR被廣泛用于鋰電池功率預(yù)測(cè)和健康狀態(tài)估計(jì)[7-8]、風(fēng)電功率預(yù)測(cè)[9]、柴油機(jī)氮氧化物預(yù)測(cè)[10]以及大型核動(dòng)力渦輪發(fā)電機(jī)的渦流損耗預(yù)測(cè)[11]等領(lǐng)域．在已有的研究中, 模型多為靜態(tài)非線性映射模型,即基于過(guò)程處于穩(wěn)定工況的假設(shè)．然而, 工業(yè)過(guò)程中存在多種動(dòng)態(tài)特性, 故有必要建立動(dòng)態(tài)GPR模型, 不僅考慮當(dāng)前輸入對(duì)輸出的影響, 而且須考慮過(guò)去的輸入和輸出對(duì)系統(tǒng)當(dāng)前輸出是否產(chǎn)生影響．GPR模型的精度取決于協(xié)方差函數(shù)中參數(shù)的優(yōu)化程度．通常采用共軛梯度算法優(yōu)化GPR模型參數(shù), 但由于共軛梯度算法搜索參數(shù)最優(yōu)解過(guò)程受初值影響較大, 且存在迭代次數(shù)難確定和局部?jī)?yōu)化的缺點(diǎn), 故本文選擇粒子群優(yōu)化(particle swarm optimization, PSO)算法[12-14]優(yōu)化核函數(shù)參數(shù)．本文擬提出一種基于粒子群參數(shù)優(yōu)化的動(dòng)態(tài)GPR模型(particle swarm optimization-dynamic Gaussian process regression, PSO-DGPR)預(yù)測(cè)釜式反應(yīng)器的出料濃度, 在建模過(guò)程中首先考慮模型的階次, 確定動(dòng)態(tài)模型的輸入, 進(jìn)而建立基于粒子群參數(shù)優(yōu)化的動(dòng)態(tài)高斯過(guò)程模型, 并分別與靜態(tài)模型和未參數(shù)優(yōu)化的動(dòng)態(tài)模型進(jìn)行對(duì)比．

1 高斯過(guò)程回歸

(1)

(2)

(3)

采用最大似然估計(jì)方法確定協(xié)方差函數(shù)參數(shù)．對(duì)數(shù)似然函數(shù)可表示為

(4)

2 標(biāo)準(zhǔn)粒子群算法

粒子群算法中種群由若干個(gè)粒子構(gòu)成, 粒子的個(gè)數(shù)即種群大小或規(guī)模．每個(gè)粒子都有一個(gè)位置矢量Xi和速度矢量Vi, 其維數(shù)代表問(wèn)題解空間的維數(shù)．假設(shè)問(wèn)題的搜索空間是一個(gè)d維空間,則第i個(gè)粒子的位置矢量和速度矢量分別為Xi=[xi1,xi2,…,xid]和Vi=[vi1,vi2,…,vid]．初始化一群隨機(jī)粒子,每個(gè)粒子在搜索空間中以一定速度飛行, 并通過(guò)迭代找到最優(yōu)解．在迭代過(guò)程中, 各粒子通過(guò)跟蹤粒子本身所找到的最優(yōu)解Pbest和整個(gè)種群目前找到的最優(yōu)解Gbest這兩個(gè)“極值”不斷進(jìn)行位置更新．各粒子的速度和位置更新公式如下:

Vid(t+1)=Vid(t)+c1r1(Pbestid(t)-Xid(t))+c2r2(Gbestid(t)-Xid(t)),

(5)

Xid(t+1)=xid(t)+Vid(t+1),

(6)

式中加速因子c1,c2為正常數(shù), 目的是使粒子具有自我總結(jié)和向群體中優(yōu)秀個(gè)體學(xué)習(xí)的能力, 從而向自己的歷史最優(yōu)點(diǎn)及群體的全局最優(yōu)點(diǎn)靠近,c1,c2用于表示個(gè)體經(jīng)驗(yàn)和群體信息在飛行中的權(quán)重時(shí)須根據(jù)問(wèn)題分析取值, 通常c1=c2=2;r1,r2為[0,1]之間的隨機(jī)數(shù)．

3 PSO-DGPR建模

3.1 模型定階

單純的靜態(tài)映射無(wú)法滿足動(dòng)態(tài)過(guò)程中模型對(duì)輸出的預(yù)測(cè)要求, 故采用動(dòng)態(tài)預(yù)測(cè)模型進(jìn)行建模, 同時(shí)考慮當(dāng)前輸入對(duì)輸出的影響以及過(guò)去輸入和輸出對(duì)系統(tǒng)當(dāng)前輸出的影響．動(dòng)態(tài)模型在建模時(shí)的首要問(wèn)題是確定模型的階次, 假設(shè)系統(tǒng)的輸出y(t)與過(guò)去l個(gè)輸入和k個(gè)輸出之間服從非線性關(guān)系, 則l和k的值即待確定的模型階次．輸入與輸出服從如下非線性關(guān)系:

y(t)=f(y(t-1),…,y(t-k+1),x(t),…,x(t-l+1))+ε(t),

(7)

其中f(·)為高斯過(guò)程擬合的非線性函數(shù),ε(t)為白噪聲．

常用的模型選擇準(zhǔn)則有基于信息論的Akaike信息準(zhǔn)則(Akaike information criterion, AIC)[15]和基于Bayes方法的Bayes信息準(zhǔn)則(Bayesian information criterion, BIC)[16]等．使用這類準(zhǔn)則時(shí), 通過(guò)計(jì)算不同參數(shù)數(shù)量模型的AIC值和BIC值來(lái)確定最終模型．當(dāng)模型參數(shù)數(shù)量增多時(shí), 易導(dǎo)致AIC值過(guò)大且模型過(guò)于復(fù)雜, 從而增大似然函數(shù), 產(chǎn)生過(guò)擬合現(xiàn)象．然而, BIC因考慮樣本數(shù)量, 有效避免了模型精度高所致模型復(fù)雜度過(guò)高的問(wèn)題．BIC準(zhǔn)則的計(jì)算式如下:

(8)

3.2 模型建立

基于粒子群算法參數(shù)優(yōu)化的動(dòng)態(tài)高斯過(guò)程回歸預(yù)測(cè)建模步驟如下:

1) 對(duì)反應(yīng)釜內(nèi)冷卻液流量和出料濃度數(shù)據(jù)進(jìn)行歸一化處理, 并將數(shù)據(jù)集劃分為訓(xùn)練集、驗(yàn)證集和測(cè)試集,其中訓(xùn)練集用于訓(xùn)練模型, 驗(yàn)證集用于確定模型階次,測(cè)試集用于檢驗(yàn)最終模型的預(yù)測(cè)性能;

2) 利用BIC確定動(dòng)態(tài)高斯過(guò)程回歸模型的階數(shù), 計(jì)算輸出y(t)分別與過(guò)去l個(gè)輸入和k個(gè)輸出之間BIC值取最小時(shí)對(duì)應(yīng)的階數(shù)l和k值, 即得到最優(yōu)階數(shù);

圖1 動(dòng)態(tài)GPR模型結(jié)構(gòu)Fig.1 Dynamic GPR model structure

3) 采用粒子群算法優(yōu)化GPR核函數(shù)超參數(shù)．初始化粒子群位置和速度,計(jì)算每個(gè)粒子的適應(yīng)度．適應(yīng)度函數(shù)取高斯模型中的對(duì)數(shù)似然函數(shù)(4)．根據(jù)適應(yīng)度更新每個(gè)粒子的歷史最佳位置Pbest和全局最佳位置Gbest, 并根據(jù)式(5)(6)更新各粒子的速度和位置．通過(guò)判斷算法是否達(dá)到最大迭代次數(shù)或全局最優(yōu)位置是否滿足最小界限來(lái)確定算法是否停止;

4) 按照?qǐng)D1所示的模型結(jié)構(gòu), 利用參數(shù)優(yōu)化的動(dòng)態(tài)高斯過(guò)程回歸模型對(duì)測(cè)試樣本{y(t-1),…,y(t-k+1),x(t),x(t-1),…,x(t-l+1)}進(jìn)行預(yù)測(cè)．

4 實(shí)驗(yàn)驗(yàn)證

連續(xù)攪拌反應(yīng)釜數(shù)據(jù)來(lái)源于DaISy數(shù)據(jù)集．輸入為反應(yīng)釜內(nèi)冷卻液流量q/(L·min-1), 對(duì)應(yīng)的輸出分別為出料濃度ca/(mol·L-1)和出料溫度T/K, 本文以出料濃度為例作為預(yù)測(cè)模型的輸出．樣本的采樣時(shí)間為0.1 min, 樣本個(gè)數(shù)為7 500, 其中50%樣本作為訓(xùn)練集, 剩余的50%樣本劃分為驗(yàn)證集和測(cè)試集．為了對(duì)模型的預(yù)測(cè)性能進(jìn)行評(píng)估, 采用均方根誤差(root mean square error, RMSE)、平均絕對(duì)誤差(mean absolute error, MAE)和平均絕對(duì)百分比誤差(mean absolute percentage error, MAPE)作為評(píng)價(jià)標(biāo)準(zhǔn):

(9)

(10)

(11)

計(jì)算GPR模型的階次, 對(duì)應(yīng)不同階次的BIC值如圖2所示．由圖2可見(jiàn), 當(dāng)k=3,l=3時(shí), 對(duì)應(yīng)的BIC值均為最小值．由此確定動(dòng)態(tài)模型的輸入y(t)=f(y(t-1),y(t-2),y(t-3),x(t-1),x(t-2),x(t-3))+ε(t).

PSO-DGPR模型中種群規(guī)模取80, 最大迭代步數(shù)為50, 不同階次對(duì)應(yīng)的預(yù)測(cè)誤差如表1所示．由表1可見(jiàn), 動(dòng)態(tài)模型的預(yù)測(cè)結(jié)果均優(yōu)于靜態(tài)映射模型, 當(dāng)取3階歷史輸入和3階歷史輸出作為當(dāng)前輸入(即k=3,l=3)時(shí)預(yù)測(cè)誤差最?。o態(tài)模型及三階動(dòng)態(tài)高斯過(guò)程模型的預(yù)測(cè)結(jié)果分別如圖3～4所示．由圖3～4可見(jiàn), 由于觀測(cè)樣本值較小, 故模型預(yù)測(cè)誤差值較小, 但動(dòng)態(tài)模型的預(yù)測(cè)誤差明顯低于靜態(tài)預(yù)測(cè)模型的．

圖2 階次k(a), l(b)的BIC值Fig.2 The identification of order k (a) and l (b)

表1 不同階次下的預(yù)測(cè)誤差

圖3 靜態(tài)映射預(yù)測(cè)結(jié)果Fig.3 Static mapping prediction results

圖4 三階動(dòng)態(tài)高斯過(guò)程模型預(yù)測(cè)結(jié)果Fig.4 Prediction results of third-order dynamic Gaussian process model

為了驗(yàn)證本文方法的有效性, 將PSO-DGPR與GPR進(jìn)行對(duì)比實(shí)驗(yàn), 采用梯度下降法求取GPR模型參數(shù), 并選擇三階動(dòng)態(tài)輸入作為模型的輸入, 預(yù)測(cè)結(jié)果如圖5所示．由圖5可見(jiàn), 在數(shù)據(jù)平滑性較好的情況下,采用動(dòng)態(tài)輸入的GPR模型和PSO-DGPR模型均能較好地反映數(shù)據(jù)變化趨勢(shì); 在數(shù)據(jù)波動(dòng)較大的情況下, 基于參數(shù)優(yōu)化的PSO-DGPR模型能夠更好地應(yīng)對(duì)數(shù)據(jù)的非線性特性,其靈活性更高,預(yù)測(cè)效果更佳．

5 結(jié)語(yǔ)

本文綜合考慮過(guò)去的輸入和輸出對(duì)當(dāng)前輸出的影響, 提出了一種基于粒子群算法優(yōu)化的動(dòng)態(tài)高斯過(guò)程回歸模型預(yù)測(cè)CSTR出料濃度．通過(guò)BIC方法確定模型階次, 并構(gòu)建動(dòng)態(tài)輸入．利用粒子群算法代替共軛梯度法優(yōu)化GPR模型的核函數(shù)超參數(shù),構(gòu)建了基于PSO-DGPR的預(yù)測(cè)模型．實(shí)驗(yàn)結(jié)果表明, 所提方法能有效提高CSTR出料濃度預(yù)測(cè)模型的精度, 對(duì)于其他工業(yè)時(shí)間序列的預(yù)測(cè)起到一定的啟發(fā)作用．