楊仁建
(中國石化石油化工科學研究院,北京100083)
現(xiàn)代工業(yè)技術的不斷發(fā)展對產品質量和安全生產的要求越來越嚴格,這就需要通過各種各樣的檢測儀表實現(xiàn)對生產過程參數(shù)的實時監(jiān)控[1]。然而,由于技術或經濟的原因,導致很多指標無法在線檢測或者檢測耗時長、滯后嚴重,造成產品質量不合格,帶來巨大的經濟損失和安全隱患[2]。因此,在線實時檢測過程參數(shù),具有重要的應用意義。
為了更好地利用時延變量優(yōu)化選取和多模型預測的優(yōu)勢,本文提出了一種基于多元模型融合和時延變量選擇的動態(tài)軟測量方法。該方法一方面根據(jù)信息熵定義聯(lián)合互信息指標,通過螺旋優(yōu)化算法搜索得到最優(yōu)的時延變量參數(shù),另一方面將高斯模型和正則極限學習機RELM (regular extreme learning machine)相結合,綜合兩者的長處,引入融合更新機制,實現(xiàn)軟測量,最終將提出的方法用于脫丁烷塔塔底的丁烷(C4H10)體積分數(shù)預測中。
針對工業(yè)過程參數(shù)檢測,常用的解決方法有兩種: 一是從儀器設備角度不斷地研發(fā)新型儀表,提高設備的應用場景和檢測精度,但是存在設備成本高、維護困難、分析周期長且滯后時間長的問題;二是通過測量與待檢測參數(shù)相關聯(lián)的變量,間接得到目標檢測變量[3]。軟測量技術通過檢測易測變量,發(fā)掘其與目標變量之間的關系,構造預測模型,間接地完成對目標變量的實時估計。由于其成本低、易于實現(xiàn),因此在工業(yè)過程中得到廣泛采用。
圖1 動態(tài)軟測量建模示意
然而,工業(yè)過程參數(shù)往往具有時延和非線性等特性,常規(guī)的軟測量方法難以滿足預測精度。由于操作條件的變化、催化劑活性降低、機械磨損、進料成分變化以及季節(jié)變化等因素,導致原有的軟測量模型將不再適應新的工況而出現(xiàn)模型老化現(xiàn)象,因此無法準確預測當前狀態(tài)信息[4]。常用的動態(tài)建模方法通過將輔助變量的歷史測量值(時延變量)有選擇性地加入到輸入中,將原始時變動態(tài)建模問題轉化為靜態(tài)問題,從而建立動態(tài)軟測量模型[5],常用的方法包括反饋網(wǎng)絡建模方法、多點輸入建模方法、多模型結構建模方法、動態(tài)加權輸入建模方法等[6-7]。Galicia等通過選取固定長度的時延變量,引入輸入建立模型,對軟測量方法進行改進,但是該方法難以保證不同工況下的建模精度[8]。Osorio等采用神經網(wǎng)絡改進軟測量方法,通過將輸出反饋回輸入增廣輸入矩陣,具有較好的預測效果[9]。文獻[10]提出了一種基于互信息和最小二乘支持向量機的軟測量建模方法,采用互信息描述變量之間的相關性,實現(xiàn)了對水泥生料細度的預測。文獻[11]提出了一種動態(tài)校正的多模型軟測量建模方法,通過將高斯模型和高斯過程回歸結合,采用自適應進行反饋校正,較好地實現(xiàn)了硫回收裝置的體積分數(shù)估計。關于時延變量的選取,目前主要采用以下兩種方法:
1)相關性分析法。文獻[11—12]基于互信息對變量的相關性進行分析,同歸最大互信息指標,確定時延變量的長度。
2)啟發(fā)式搜索。根據(jù)時延變量建立優(yōu)化模型,通過啟發(fā)式搜索算法,如遺傳算法[13]、粒子群算法[14]、差分進化算法[15]等,獲取誤差最小的時延變量維度,文獻[16]通過構造合適的適應度函數(shù),將時延變量參數(shù)估計問題轉化為一個多維非線性優(yōu)化問題,進而采用混合差分進化算法得到最優(yōu)的時延變量個數(shù),進行軟測量建模,在常壓塔航空煤油閃點預測中取得了較好的應用效果。
對于一個工業(yè)過程,假設系統(tǒng)的數(shù)據(jù)集{x(k),y(k)},k∈[1,n],其中x(k)=[x1(k),x2(k), …,xm(k)]∈Rm與y(k)∈R分別表示k時刻的輸入和輸出??紤]系統(tǒng)延時,對k時刻的輸入進行修正:
(1)
綜上所述,對于給定的數(shù)據(jù)集,只需針對每個變量從候選時延變量中選取從k-di時刻開始的連續(xù)的li個主導變量,得到用于建模的輸入輸出數(shù)據(jù);然后基于數(shù)據(jù)進行建模,即可得到軟測量模型。針對以上情況,本文一方面提出時延變量選擇策略,基于聯(lián)合互信息最大化,采用螺旋優(yōu)化算法尋優(yōu),得到最佳時延參數(shù)di和數(shù)據(jù)長度li,確定建模數(shù)據(jù)集;另一方面引入多元模型融合策略,將高斯模型和RELM相結合,建立軟測量模型。
互信息能夠充分地描述變量之間的相互關系,被廣泛地應用于系統(tǒng)分析、變量選擇、圖像處理等方面。然而,常規(guī)的互信息估計方法計算量大,算法過程復雜、精度低,不能有效地求解高維問題。k-近鄰互信息估計方法是由Kraskov等針對高維變量求解提出的互信息估計方法,具有計算簡單、精度高等優(yōu)點[17-18]。本文針對k-近鄰互信息估計方法對時延變量進行分析。
由文獻[17]可知,對于變量x,y構成的空間z=(x,y),按照zi=(xi,yi)與其他點的距離進行排序,則定義εi/2為該點到其本身的k-近鄰距離,其中εx(i)/2和εy(i)/2分別為到x軸和y軸相應點的距離。假設變量xi和yi的k-近鄰距離內分別有nx(i)和ny(i)個樣本點,則變量x和y的互信息可通過式(2)進行計算:
MI(x,y)=φ(k)-〈φ(nx+1)+φ(ny+1)〉+φ(N)
(2)
式中:φ(k)——雙Γ函數(shù),φ(k)=Γ(k)-1dΓ(k)/dk,其中dΓ(k)/dk表示微分,滿足迭代關系φ(k+1)=φ(k)+1/k,且φ(1)≈-0.577 215 6,符號〈g〉表示對所有變量i∈[1,n]取平均,即:
(3)
針對一般問題,將式(2)拓廣到m維變量,則變量(x1,x2, …,xm)的互信息可表示為
MI(x1,x2, …,xm)=φ(k)-
〈φ(nx1)+…+φ(nxm)〉+(m-1)φ(N)
(4)
k-近鄰互信息估計的計算結果受參數(shù)k的影響非常大,k越大,計算結果越精確,但是計算量會大幅增加;k越小,計算量越小,但是精度相應會變差,目前尚未有合理的方法對k進行選擇。本文基于k-近鄰互信息估計的基本思想,基于互信息定義聯(lián)合互信息計算公式:
MI(x1,x2, …,xm;y)=MI(x1,x2, …,xm,y)-MI(x1,x2, …,xm)
(5)
式中:MI(x1,x2, …,xm,y)——輸入和輸出變量的高維互信息;MI(x1,x2, …,xm)——各輸入變量的高維互信息,兩者均可由式(4)得到。
由式(5)可知,MI(x1,x2, …,xm,y)越大,參數(shù)反映的輸入和輸出信息越多,輸入和輸出之間的相關性越大;MI(x1,x2, …,xm)越小,各輸入變量的相關性越小,得到的聯(lián)合互信息參數(shù)MI(x1,x2, …,xm;y)越大。因此,只要求得滿足MI(x1,x2, …,xm;y)最大時各輸入變量的相關參數(shù),即可確定時延變量。結合圖1中給出的模型延時變量,得到如下最優(yōu)化問題:
(6)
螺旋優(yōu)化算法是在模仿自然界中螺旋現(xiàn)象的基礎上提出的,通過個體圍繞中心點旋轉搜索,逐步逼近最優(yōu)解。文獻[19]給出了一種改進的螺旋優(yōu)化算法,通過引入自適應柯西變異和拉丁超立方采樣,增強了標準螺旋優(yōu)化的全局搜索能力,能有效得到全局最優(yōu)解。
對于n維空間的優(yōu)化問題,假設中心點為x*,則圍繞x*的旋轉可以描述為
x(k+1)=γTn(θ)x(k)
(7)
(8)
式中:T(n)(θ)——旋轉矩陣;θ——圍繞x*的旋轉角度,且0<θ<2π;γ——收縮系數(shù),能反應旋轉前后與x*的距離變化,且0<γ<1。
假設種群大小為pop,定義每一步迭代過程中存儲的個體最佳歷史解Pbestij(k)和全局最優(yōu)解Gbestj(k),則自適應柯西變異定義如下:
xij(k+1)=xij(k)+αχjC(0, 1)
(9)
(10)
式中:α——校正因子;C(0, 1)——由柯西分布生成的隨機數(shù)。
改進螺旋優(yōu)化算法的執(zhí)行過程如圖2所示。
圖2 改進螺旋優(yōu)化算法流程示意
對于給定的輸入輸出軟測量建模數(shù)據(jù),采用第2節(jié)和第3節(jié)的相關理論確定了時延參數(shù)后,得到新的建模數(shù)據(jù)集{xi(k-di),xi(k-di-1), …,xi(k-di-li),y(k)},k∈[1,n],為了簡化描述,下文中統(tǒng)一用{xk,yk)表示。為了提高軟測量建模的精度和泛化能力,本文提出一種多元模型融合策略基于數(shù)據(jù)集進行建模預測。
文獻[20]給出了高斯模型的基本理論,對于數(shù)據(jù)集{xk,yk},k=1, …,n,若給定1個新的輸入xn+1,則預測輸出yn+1可寫作:
(11)
δ2(xn+1)=c(xn+1)-c(xn+1)TC-1c(xn+1)
(12)
(13)
(14)
超參數(shù)θ可以通過貝葉斯推理進行估計,將高斯分布轉化為p(y|θ,x)=G(0,C),最大化的對數(shù)似然函數(shù)如下:
(15)
針對式(15)采用共軛梯度法求偏導數(shù),依次求解得到各個超參數(shù),即可確定具體的高斯模型,根據(jù)式(11)對下一時刻的輸出給出預測值。
極限學習機的本質是一種單層前饋神經網(wǎng)絡,可以通過選取輸入權重和隱層偏置對網(wǎng)絡進行逼近,能夠在保證較高精度的前提下較大地加快學習速度。文獻[21]給出了基于果嶺回歸的RELM,對于數(shù)據(jù)集{xk,yk),k=1, …,n,可通過求解如下最小二乘解得到:
(16)
(17)
式中:H——隱層輸出矩陣,它的第i列是第i個隱層節(jié)點相對于x1,x2, …,xn的輸出;G(·)——激勵函數(shù);wi=[wi1, …,win]T——隱層節(jié)點和輸入節(jié)點的權向量;βi=[βi1, …,βin]T——隱層節(jié)點和輸出節(jié)點的權向量;bi——偏置。
求解式(16)所述的優(yōu)化問題,即可得到RELM的預測輸出模型:
(18)
(19)
對于實際的工業(yè)過程,往往受各種工況條件、人為因素、機械磨損等因素影響,導致基于一次模型預測的結果不夠準確,需要不斷更新模型參數(shù),甚至更新模型。本節(jié)中,引入了模型更新機制,選取對模型有利的參數(shù)進行更新,分別針對高斯模型和RELM介紹了數(shù)據(jù)選取準則,采用遞歸求解的方法提高更新效率。
4.3.1高斯模型更新
高斯模型通過預測方差評價預測輸出的置信水平,當方差δ(x)較小時,可認為預測輸出是準確的;當方差較大時,則認為輸出不準確。本節(jié)中,通過定義方差閾值δlimit(x),選取對預測有理的樣本數(shù)據(jù),若預測方差超過閾值,則認為樣本是有利的,否則,將樣本剔除。
對于已知數(shù)據(jù)集,在加入新樣本xn+1后,相應的輸出修正為
(20)
(21)
(22)
(23)
(24)
4.3.2RELM更新
由于RELM不引入預測方差,無法對預測結果進行直接評價,這里定義誤差函數(shù):
(25)
通過設置誤差閾值elimit對模型進行評價,從而確定是否對模型進行更新。若誤差大于閾值,則認為新樣本對模型改進有利,保留;否則,刪掉樣本數(shù)據(jù)。
(26)
式中:Hn——初始隱層輸出矩陣。
當有利樣本加入訓練集后,隱層輸出矩陣更新為
(27)
(28)
根據(jù)Woodbury矩陣恒等式[22],式(28)可寫作:
(29)
輸出權重βn+1更新為
(30)
為了提高模型的預測精度和泛化能力,本文在多模型估計的基礎上,提出了一種多元模型融合策略,通過權值將高斯模型和RELM進行融合。假設模型集M={mj,j=1, 2, …,r}包含有限個模型,其中模型mj是對不同模型的描述。對于任意如下形式離散非線性系統(tǒng):
(31)
式中:xk——系統(tǒng)狀態(tài);uk——控制變量;yk——系統(tǒng)輸出;σk,μk——過程噪聲和測量噪聲,且滿足σk∶N(0, Qk), μk∶N(0,Rk)。
針對上述系統(tǒng),采用多元模型融合策略進行預測的基本步驟如下:
1)條件初始化。假設k-1時刻的匹配模型是mi,k時刻的匹配模型是mj,則在yk-1條件下的融合概率為
(32)
當j=1, 2, …,r時,初始化狀態(tài)和協(xié)方差的融合估計為
(33)
(34)
2)模型條件無跡卡爾曼濾波。采用無跡卡爾曼濾波[23],根據(jù)步驟1)中初始化的狀態(tài)和協(xié)方差,以及輸出的yk更新狀態(tài)估計。基于無跡卡爾曼濾波的基本公式如式(35)所示:
a)狀態(tài)采樣
(35)
式中:n——狀態(tài)的維度;λ——系數(shù)。
b)時間更新
(36)
(37)
(38)
其中,各系數(shù)定義如下:
(39)
(40)
(41)
c)量測更新
(42)
3)概率更新。模型概率計算公式如下:
(43)
(44)
4)多元模型融合軟測量建模
(45)
(46)
圖3 多元模型融合預測過程示意
脫丁烷塔是石油煉制生產過程中脫硫和石腦油分離裝置的重要組成部分,其主要指標是減小脫丁烷塔塔底的C4H10體積分數(shù),具體工藝流程如圖4所示,主要包括7個輸入變量和1個輸出變量,具體參數(shù)見表1所列。
圖4 脫丁烷塔工藝流程示意
x1x2x3x4x5x6x7y塔頂溫度塔頂壓力塔頂回流量塔頂產品流出量第六層塔板溫度塔底溫度1塔底溫度2C4H10體積分數(shù)
對于該工業(yè)過程,C4H10體積分數(shù)無法在塔底物料流出處檢測,通常通過檢測塔頂?shù)漠愇焱楫a物間接獲得。該過程產物復雜,呈現(xiàn)明顯的非線性特性,另外化工過程的反應時間、物料傳遞等都加大了過程參數(shù)的滯后,無法有效地監(jiān)測C4H10體積分數(shù)。因此,為了提高脫丁烷塔的控制品質,需要基于歷史生產數(shù)據(jù)對C4H10體積分數(shù)進行動態(tài)軟測量建模。
為了方便處理,由于輸入變量x6,x7均為塔底溫度,將2個變量通過計算均值簡化為1個變量。對實際工業(yè)過程進行采樣,得到具有2 394個樣本的數(shù)據(jù)集,將所有樣本平均分成2份,前1份作為訓練集,后1份作為測試集,所有樣本采樣時間均為12 min。根據(jù)專家知識,脫丁烷塔C4H10體積分數(shù)檢測儀表均存在15 min的測量周期,獲得C4H10體積分數(shù)值存在30~75 min的滯后。因此,輸入變量與系統(tǒng)輸出C4H10體積分數(shù)存在45~90 min的時間滯后[25]。
采用本文提出的基于多元模型融合和時延變量選擇的動態(tài)軟測量建模方法對脫丁烷塔工業(yè)數(shù)據(jù)進行軟測量建模,首先針對訓練集采用第3節(jié)給出的時延變量選擇策略,基于聯(lián)合互信息和螺旋優(yōu)化對時延變量進行優(yōu)化,得到時延變量的最佳個數(shù)。然后采用第4節(jié)給出的多元模型融合策略,進行建模預測。為了突出本文方法的優(yōu)越性,分別引入偏最小二乘法(PLS)和最小二乘支持向量回歸機(LSSVM)對原始數(shù)據(jù)進行軟測量建模。為了評價預測結果,分別計算三種方法預測的均方根誤差和最大絕對誤差,軟測量建模得到的塔底C4H10預測體積分數(shù)與真實體積分數(shù)的預測誤差曲線如圖5~圖10所示,軟測量具體優(yōu)化結果見表2所列。
表2 軟測量誤差統(tǒng)計
通過對比表2和圖5~圖10可以發(fā)現(xiàn),本文方法得到的軟測量預測結果,無論是誤差均方根值還是誤差最大絕對值,均比其他兩種方法要小,說明本文提出的方法具有較高的預測精度和泛化能力。由于本文通過最大聯(lián)合互信息采用螺旋優(yōu)化選取最佳的時延變量個數(shù),能夠充分發(fā)掘對預測有利的變量數(shù)據(jù),剔除不利因素,且能在保證反映系統(tǒng)信息的情況下最大限度地減少變量個數(shù)。另外,多元模型融合策略的引入,充分利用高斯模型統(tǒng)計評價的特性和RELM更新速度快的優(yōu)點,通過無跡卡爾曼濾波進行融合預測,動態(tài)的更新輸出權重,能夠充分地描述系統(tǒng)信息,增強預測精度和模型泛化能力,充分保證了軟測量建模的性能。
圖5 PLS軟測量輸出對比
圖6 PLS軟測量誤差
圖7 LSSVM軟測量輸出對比
圖8 LSSVM軟測量誤差
圖9 本文方法軟測量輸出對比
圖10 本文方法軟測量誤差
綜上可知,本文提出的基于多元模型融合和時延變量選擇的動態(tài)軟測量建模方法,具有很好的建模精度和泛化能力,能夠有效地解決復雜非線性、帶時延的工業(yè)過程動態(tài)軟測量建模問題。
本文提出了一種基于多元模型融合和時延變量選擇的動態(tài)軟測量方法,能有效地解決工業(yè)過程中難以測量參數(shù)的軟測量問題,方法的主要結論如下:
1)引入聯(lián)合互信息評價指標,通過螺旋優(yōu)化對時延變量參數(shù)進行尋優(yōu),最大程度上反映了系統(tǒng)信息,減小了變量個數(shù)。
2)引入多元模型融合策略,充分發(fā)揮了高斯模型和RELM在建模方面的優(yōu)勢,將統(tǒng)計評價和最小范數(shù)求解相融合,通過無跡卡爾曼濾波進行參數(shù)更新,動態(tài)更新軟測量模型,提高了模型精度和泛化能力。
3)本文的軟測量建模方法對脫丁烷塔底C4H10體積分數(shù)的軟測量取得了很好效果,但是所提方法不僅適用于帶時延的工業(yè)過程,對于科研和生活中常見的建模問題,也具有較好的建模精度。
4)本文僅考慮了多入單出系統(tǒng)的軟測量建模,對于多入多出系統(tǒng),本文方法同樣適用,不過計算復雜度會增加,后續(xù)研究可以從多入多出系統(tǒng)建模效率展開。