項亞南，潘 豐

(輕工過程先進控制教育部重點實驗室，江南大學(xué)，江蘇無錫 214122)

?

基于多元統(tǒng)計時間序列模型的間歇過程故障預(yù)測方法研究

項亞南，潘 豐

(輕工過程先進控制教育部重點實驗室，江南大學(xué)，江蘇無錫 214122)

青霉素發(fā)酵過程是一個典型的間歇過程并且生產(chǎn)過程伴隨著非線性，多階段，動態(tài)性能等特點，一旦過程中發(fā)生故障，發(fā)酵底物可能會被破壞，造成巨大的財產(chǎn)損失，因此提前預(yù)測間歇過程的故障是非常有意義的。針對間歇過程多元統(tǒng)計(MPCA)故障診斷的研究方法，文中提出一種時間序列模型與MPCA相結(jié)合的方法來預(yù)測緩變故障，在建立和SPE統(tǒng)計量的自回歸模型之后，再與各自建立的控制限對比，判斷是否超限來預(yù)測漸變故障。文中采用的是Pensim 2．0仿真平臺對該方法進行仿真研究，仿真結(jié)果驗證了該方法的有效性。

青霉素發(fā)酵；間歇過程；多元統(tǒng)計；故障診斷；時間序列模型；控制限

0 引言

青霉素發(fā)酵是一個典型的間歇過程，由于初始條件(如菌體濃度、初始底物的濃度、菌體的初始活性等)和操作條件(例如補料的數(shù)量、速率)的不同，各批次之間的過程變量往往表現(xiàn)出不同的特性。另外即使在相同批次發(fā)酵過程的各個階段(菌體生長階段、產(chǎn)物合成階段)之間，過程變量的自身特性和變量之間的相關(guān)關(guān)系也都會變化[1-2]，發(fā)酵過程中一旦發(fā)生故障將產(chǎn)生很嚴重的后果，但是如果可以根據(jù)在線的數(shù)據(jù)，對未來一段時間的故障作出早期的預(yù)測，將給操作員爭取更多的時間來排除故障，保證安全生產(chǎn)。在以往的故障預(yù)測研究中：基于模型的卡爾曼濾波器方法[3]需要建立間歇過程的數(shù)學(xué)模型，間歇過程復(fù)雜的機理很難建立準確的模型；神經(jīng)網(wǎng)絡(luò)的故障診斷適用于有大量生產(chǎn)的經(jīng)驗和專家知識的場合[4]，而間歇過程在這方面的知識和專家知識是有限的。

通過上述分析，考慮到前后數(shù)據(jù)的相關(guān)性關(guān)系，本文提出以當(dāng)前時段變量的變化趨勢為基礎(chǔ)的間歇過程監(jiān)控方法，采用時間序列(ARMA)與多元統(tǒng)計方法MPCA[5]相結(jié)合的方法來預(yù)測間歇過程的緩變故障。

1 時間序列預(yù)測方法與MPCA

1.1 ARMA與MPCA

ARMA是在20世紀70年代提出的[6]，開始是用于對市場經(jīng)濟的預(yù)測，隨著研究的不斷深入，廣泛的應(yīng)用在各個不同的領(lǐng)域之中。

過程的數(shù)據(jù)看做是一個隨機的序列，但是相鄰的數(shù)據(jù)具有相關(guān)性和依賴性，根據(jù)輸入輸出數(shù)據(jù)，建立相互之間的映射關(guān)系，以線性自回歸滑動平均模型為基礎(chǔ)建立數(shù)據(jù)的模型，按照時間建立一個隨時間變化變量之間的相關(guān)關(guān)系，借以預(yù)測下一個時間段或者以后若干時間內(nèi)數(shù)據(jù)可能達到的水平。

間歇過程與連續(xù)生產(chǎn)過程不同在于特殊的運行機制和工業(yè)要求，并且產(chǎn)品質(zhì)量更易受到環(huán)境、投放的材料、設(shè)備運行的狀況等因素的影響。針對間歇過程的特點，本文采用的是基于數(shù)據(jù)的主元分析方法來分析，這種監(jiān)控方法不需要建立過程的數(shù)學(xué)模型，通過挖掘輸入輸出數(shù)據(jù)中的信息以此來建立過程的統(tǒng)計量指標，就可以進行故障診斷。

MPCA方法在間歇生產(chǎn)中得到了較為廣泛的研究與應(yīng)用[7-9]，它的基本原理是：

X=I·J·K

式中：I表示批次；J表示變量個數(shù)；K表示采樣時刻序列。

首先把正常的三維數(shù)據(jù)按照一個維度分解為二維的空間，再用一般的PCA進行主元分析，分解為主元空間和殘差空間。在正常的工況數(shù)據(jù)下分別建立T2統(tǒng)計量和平方預(yù)測誤差SPE (Squared prediction error)統(tǒng)計量控制限，其中T2是根據(jù)主元空間建立的，表示的某時刻主元模型描述原模型的程度，SPE是根據(jù)殘差空間建立的，表示的是某時刻數(shù)據(jù)在殘差空間偏離主元模型的程度。

基于MPCA的時間序列的預(yù)測方法來預(yù)測故障，即是預(yù)測下一批次數(shù)據(jù)的T2統(tǒng)計量與SPE統(tǒng)計量是否在控制限之內(nèi)。

采集到數(shù)據(jù)之后，把過程的數(shù)據(jù)均值化處理之后再轉(zhuǎn)化為一組平穩(wěn)的時間序列數(shù)據(jù)。利用當(dāng)前的數(shù)據(jù)預(yù)測未來的趨勢，首先選擇合適的擬合模型，即確定模型的階數(shù)和辨識出未知的參數(shù)，檢驗預(yù)測。

1．2 時間序列模型模型的確定

時間序列模型，首先是要確定選用AR、MA或ARMA模型，本文利用的是自相關(guān)系數(shù)和偏自相關(guān)系數(shù)來確定模型的類型。

ARMA是一種由AR(自回歸模型)、MA(滑動平均模型)相結(jié)合的模型預(yù)測方法，預(yù)測模型的表達式如下：

(1)

式中：p、q表示的分別是自回歸和滑動平均的階數(shù)，簡記為ARMA(p,q)；φi為自回歸系數(shù)；θj為滑動平均系數(shù)，這兩項是待定系數(shù)。

當(dāng)p=0時，模型為MA模型，當(dāng)q=0時模型為AR模型。對于ARMA，要找到與其擬合最好的模型。

1.2.1 模型階數(shù)的確定

對時間序列{xt|t=1,2,…,n}，其中統(tǒng)計分析如下

均值：

(2)

方差：

(3)

協(xié)方差：

(4)

自相關(guān)系數(shù)ACF(ρk)：

(5)

偏相關(guān)系數(shù)PACF(φkk)就是在已知序列值的情況下，對未知序列的一種量度:

φ11=ρ1

(6)

(7)

φki=φk-1,i-φkkφk-1,k-1；i=1,2…

(8)

模型類別的判斷依據(jù)是：對于AR模型偏相關(guān)系數(shù)φkk是截尾的，即若偏相關(guān)系數(shù)在n>p后，φkk=0即可判斷AR模型的階數(shù)p=n-1，對于MA模型，自相關(guān)系數(shù)ρk是截尾的，當(dāng)m>p時，求得ρk=0時，即可判斷MA模型的階數(shù)是q=m-1。

對于ARMA模型，采用的是準則函數(shù)方法來確定模型的階數(shù)。采用一種最小AIC準則函數(shù)作為定階方法[10]，即可以利用其數(shù)據(jù)的自相關(guān)函數(shù)來判斷階數(shù)，AIC的表達式為：

(9)

從低階到高階分別對p,q取值分別建立模型，在進行參數(shù)估計，比較各模型的AIC值，當(dāng)AIC取得極小值的時對應(yīng)的p,q確定為最佳的模型參數(shù)。

1.2.2 參數(shù)的估計

在確定了模型的階數(shù)之后，接下來的問題是確定未知的參數(shù)，文中采用的是最小二乘的方法來辨識未知的參數(shù)。

采用模型預(yù)測的算法流程總結(jié)如下，流程圖如圖1所示。

(1) 零均值化處理，0-1均值化處理；

(2) 平穩(wěn)信號處理，ARMA要求時間序列是平穩(wěn)的信號，既非平穩(wěn)信號做差分處理；

(3) 分別計算在不同p和q取值下自相關(guān)系數(shù)和偏相關(guān)系數(shù)，確定模型的階數(shù)p和q，選擇合適的ARMA模型；

(4) 采用最小二乘的方法辨識參數(shù)φi和θj；

(5) 白噪聲對模型進行測試，并對模型不斷地改進與優(yōu)化，重復(fù)步驟(3)直到滿足要求；

(6) 在線應(yīng)用，對將來某一時刻的數(shù)據(jù)作出預(yù)測。

圖1 時間序列建模流程圖

2 基于統(tǒng)計量T2和SPE時間序列模型對故障的預(yù)測

根據(jù)正常數(shù)據(jù)分別建立T2統(tǒng)計量與SPE統(tǒng)計量的控制限，再根據(jù)在線預(yù)測出兩個統(tǒng)計量的的未來值，若兩個統(tǒng)計量都在在控制限范圍內(nèi)認為過程是正常的，當(dāng)任一統(tǒng)計量超過控制限時，就判斷過程發(fā)生了故障。

利用Pensim2．0仿真平臺得到30批次正常工況下的數(shù)據(jù)，前24組用來建立預(yù)測模型，后6組用來檢測預(yù)測模型，在青霉素發(fā)酵18個變量中選取10個起主要作用的變量，共80個樣本，選擇的數(shù)據(jù)樣本是X(30×10×80)，利用MPCA方法沿著采樣時刻方向分解為二維矩陣X(30×800)，建立正常數(shù)據(jù)的T2和SPE統(tǒng)計量控制限。

利用AIC準則函數(shù)和最小二乘的方法步驟，分別計算出自相關(guān)系數(shù)ACF(ρk)：和偏相關(guān)系數(shù)PACF(φffkk)，發(fā)現(xiàn)偏相關(guān)系數(shù)是截尾的，而自相關(guān)系數(shù)是沒有截尾的，下面分別測試3階、4階5階和6階預(yù)測模型并計算出對應(yīng)的AIC值。

表1 T2的ARMA模型階數(shù)和AIC值

由表1可知，由AIC最小值原則，選擇模型的階數(shù)為4階。那么通過時間序列模型建立的T2的表達式為

Tt=0.523Tt-1+0.296Tt-2+0.127Tt-3+0.054Tt-4

(10)

通過上述的模型表達式，用正常數(shù)據(jù)后6組的數(shù)據(jù)對模型進行測試，實際值和預(yù)測值如表2所示。

表2 T2實際值和預(yù)測值

如圖2所示，T2的時間序列模型得到的預(yù)測值對現(xiàn)場的實際值有很好的擬合效果。

圖2 T2實際值與預(yù)測值

同樣的方法建立SPE的數(shù)據(jù)模型，按照上面類似的做法，這里選擇的依然是AR模型，分別采用2階、3階、4階和5階做測試，各參數(shù)如表3所示。

表3 SPE的ARMA模型階數(shù)和AIC值

表3可知，當(dāng)階數(shù)選為5階時，AIC的值取最小值此時得到SPE統(tǒng)計量的時間序列模型表達式為

SPEt=0.425SPEt-1+0.317SPEt-2+0.205SPEt-3+0.037SPEt-4+0.016SPEt-5

(11)

類似的，把正常數(shù)據(jù)和SPE的時間序列模型預(yù)測的數(shù)據(jù)作比較，實際值和預(yù)測值如表4所示。

表4 SPE統(tǒng)計量實際值和預(yù)測值

同樣的，如圖3所示，統(tǒng)計量SPE的時間序列模型，對實際值有很好的預(yù)測的效果，擬合的效果較好。

圖3 SPE實際值與預(yù)測值

3 仿真研究

采用青霉素發(fā)酵Pensim2．0仿真平臺對該方法仿真，在建立合適的ARMA和SPE的時間序列模型之后，為了驗證該方法在整個過程中的有效性，再選取30個批次，10個變量，125個采樣時刻的數(shù)據(jù)，組成X(30×10×125)，在80 h的時候加入一個底物流加速度+5%的階躍擾動作為測試。和SPE仿真預(yù)測圖，如圖4、圖5所示。

圖4 T2仿真預(yù)測圖

圖5 SPE仿真預(yù)測圖

在第80個采樣時刻的時刻發(fā)生故障，有圖4、圖5的仿真圖可知，T2到95個采樣時刻的時候超過控制限，此時檢測出故障發(fā)生，延時15 h左右，預(yù)測的模型的預(yù)測值才有變化，故障的預(yù)測值越來越大，SPE在90 h左右檢測出發(fā)生了故障。故障的發(fā)生有延時的，這是因為底物是菌體生長和產(chǎn)物合成所必須的物質(zhì)，但是它造成的影響是間接的，通過如培養(yǎng)基的體積，氧氣含量，CO2的濃度等表現(xiàn)出來。符合實際的生產(chǎn)情況，達到了故障檢測的效果。

4 結(jié)論

本文基于ARMA時間序列預(yù)測模型，結(jié)合MPCA建立間歇過程的和SPE統(tǒng)計量的預(yù)測模型來預(yù)測緩變故障的發(fā)生，采用青霉素發(fā)酵Pensim2．0仿真平臺對該方法進行驗證，仿真結(jié)果表明該方法是有效的，對實際的生產(chǎn)具有很好的指導(dǎo)作用。

[1] SEBASTIEN R,ANDREI D,JACKY D．Detection and characterization of physiological states in bioprocesses based on holder exponent．Knowledge-Based Systems,2008,21(1)：70-79．

[2] WEI S N,YAN M G,AHMET P ,et al．A method for multiphase batch process monitoring based on auto phase identification ．Journal of Process Control,2011,21(8)：627-638．

[3] 韓冰潔，緱林峰．，基于卡爾曼濾波器組的執(zhí)行機構(gòu)故障診斷方法．計算機仿真,2013,30(2)：293-296．

[4] WANG H Q,CHEN P．Intelligent diagnosis method for rolling element bearing faults using possibility theory and neural network ．Computers&IndustrialEngineering,2011,60(4),511-518．

[5] NOMIKOS P，MACGREGOR J F．Monitoring batch process using multi—way principle component analysis．American Institute of Chemical Engineers,1994,40(8)：136-1375．

[6] BOX G E,JENKINS G M．Time Series Analysis Forecasting and Contro(academic thesis)l．San Francisco：Holden-Day Inc,1970．

[7] 張成,李元．基于統(tǒng)計模量分析間歇過程故障檢測方法研究．儀器儀表學(xué)報,2013,9(6)：2103-2110．

[8] CHUANG W L ,CHEN C H,YEN J Y,et al．Using MPCA of spectra model for fault detection in a hot strip mill．Journal of Materials Processing Technology,2009,209(8)：4162-4168．

[9] 汪志鋒，袁景淇．動態(tài)MPCA在發(fā)酵過程監(jiān)測與故障診斷中的應(yīng)用．生物工程學(xué)報,2006,22(3)：483-488．

[10] EKONOMOU S P L,KARAMOUSANTAS D C,CHATZARAKIS G E,et al．Liatsis．Electricity demand loads modeling using AutoRegressive Moving Average (ARMA) models．Energy,2008,33(9)：1353-1360．

Batch Processes Fault Prediction Method Research Based onMultivariate Statistical Analysis Time Series Model

XIANG Ya-nan,PAN Feng

(Key Laboratory of Advanced Process Control for Light Industry,Jiangnan University,Wuxi 214122,China)

Penicillin fermentation process is a typical batch process along with nonlinear,multi-stage,dynamic performance and other characters．If a fault occurs in the process of fermentation,substrate will be destroyed,which will lead to a huge property loss ．So it is meaningful to predict the faults during the batch process．This article put forward a kind of time series model with MPCA to predict graded fault,and compared with the control limit of each statistic after establishing predict autoregressive model of T2 and SPE．This article simulated the method on Pensim 2．0 penicillin fermentation platform．The results confirm the effectiveness of this method．

penicillin fermentation；batch process；MPCA；fault diagnosis；time series mode；control limit

李世興(1986—)，博士研究生，主要研究領(lǐng)域為工業(yè)自動化和物聯(lián)網(wǎng)。E-mail：lisx@sia．cn 王宏(1963—)，研究員，碩士，主要研究領(lǐng)域為嵌入式系統(tǒng)和計算機網(wǎng)絡(luò)。E-mail：wang@sia．cn

國家自然科學(xué)基金資助項目(61273131)；江蘇省產(chǎn)學(xué)研聯(lián)合創(chuàng)新資金項目(BY2013015-39)

2014-06-19 收修改稿日期：2015-03-06

TP273

A

1002-1841(2015)06-0135-03