項(xiàng)亞南，潘 豐

(輕工過(guò)程先進(jìn)控制教育部重點(diǎn)實(shí)驗(yàn)室，江南大學(xué)，江蘇無(wú)錫 214122)

?

基于多元統(tǒng)計(jì)時(shí)間序列模型的間歇過(guò)程故障預(yù)測(cè)方法研究

項(xiàng)亞南，潘 豐

(輕工過(guò)程先進(jìn)控制教育部重點(diǎn)實(shí)驗(yàn)室，江南大學(xué)，江蘇無(wú)錫 214122)

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

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

0 引言

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

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

1 時(shí)間序列預(yù)測(cè)方法與MPCA

1.1 ARMA與MPCA

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

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

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

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

X=I·J·K

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

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

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

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

1．2 時(shí)間序列模型模型的確定

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

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

(1)

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

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

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

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

均值：

(2)

方差：

(3)

協(xié)方差：

(4)

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

(5)

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

φ11=ρ1

(6)

(7)

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

(8)

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

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

(9)

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

1.2.2 參數(shù)的估計(jì)

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

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

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

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

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

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

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

(6) 在線應(yīng)用，對(duì)將來(lái)某一時(shí)刻的數(shù)據(jù)作出預(yù)測(cè)。

圖1 時(shí)間序列建模流程圖

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

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

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

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

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

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

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

(10)

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

表2 T2實(shí)際值和預(yù)測(cè)值

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

圖2 T2實(shí)際值與預(yù)測(cè)值

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

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

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

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

(11)

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

表4 SPE統(tǒng)計(jì)量實(shí)際值和預(yù)測(cè)值

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

圖3 SPE實(shí)際值與預(yù)測(cè)值

3 仿真研究

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

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

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

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

4 結(jié)論

本文基于ARMA時(shí)間序列預(yù)測(cè)模型，結(jié)合MPCA建立間歇過(guò)程的和SPE統(tǒng)計(jì)量的預(yù)測(cè)模型來(lái)預(yù)測(cè)緩變故障的發(fā)生，采用青霉素發(fā)酵Pensim2．0仿真平臺(tái)對(duì)該方法進(jìn)行驗(yàn)證，仿真結(jié)果表明該方法是有效的，對(duì)實(shí)際的生產(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í)行機(jī)構(gòu)故障診斷方法．計(jì)算機(jī)仿真,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)計(jì)模量分析間歇過(guò)程故障檢測(cè)方法研究．儀器儀表學(xué)報(bào),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] 汪志鋒，袁景淇．動(dòng)態(tài)MPCA在發(fā)酵過(guò)程監(jiān)測(cè)與故障診斷中的應(yīng)用．生物工程學(xué)報(bào),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àn)楣I(yè)自動(dòng)化和物聯(lián)網(wǎng)。E-mail：lisx@sia．cn 王宏(1963—)，研究員，碩士，主要研究領(lǐng)域?yàn)榍度胧较到y(tǒng)和計(jì)算機(jī)網(wǎng)絡(luò)。E-mail：wang@sia．cn

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

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

TP273

A

1002-1841(2015)06-0135-03