基于收斂交叉映射的化工過程擾動因果分析方法

程非凡，趙勁松

?

基于收斂交叉映射的化工過程擾動因果分析方法

程非凡，趙勁松

（清華大學(xué)化學(xué)工程系，北京100084）

當(dāng)化工裝置出現(xiàn)異常情況，操作工人往往無法及時準(zhǔn)確地定位故障發(fā)生的原因?；跀?shù)據(jù)的方法能夠通過化工過程數(shù)據(jù)，分析異常工況中的擾動傳播路徑，確定異常工況出現(xiàn)的根本原因。針對化工動態(tài)系統(tǒng)，提出了具有時間特性的收斂交叉映射方法(CCM)，和基于赤池信息準(zhǔn)則的維度選擇方法。為了驗證提出算法的有效性，在簡單的生態(tài)系統(tǒng)，因果檢測基準(zhǔn)系統(tǒng)和全混釜反應(yīng)器（CSTR）中進行驗證，并與原有的收斂交叉映射算法進行對比，體現(xiàn)出具有時間特性的收斂交叉映射算法的優(yōu)越性。

擾動分析；因果分析；收斂交叉映射；時間特性；安全；系統(tǒng)工程；模擬

引 言

化工生產(chǎn)過程是一個復(fù)雜過程。為了使經(jīng)濟效益最大化，化工生產(chǎn)過程往往要求“安、穩(wěn)、長、滿、優(yōu)”，要求裝置能夠長周期、平穩(wěn)地運行。而在實際的生產(chǎn)過程中，操作條件的改變可能會導(dǎo)致異常工況的出現(xiàn)，影響裝置的平穩(wěn)運行，甚至導(dǎo)致非計劃停車或事故，所以當(dāng)異常工況出現(xiàn)之后，就需要對異常工況進行快速診斷，采取有效的控制措施。隨著自動控制技術(shù)的不斷發(fā)展，裝置的控制與狀態(tài)變量緊緊地耦合在一起，當(dāng)一個主要控制變量發(fā)生異常的時候，可能會導(dǎo)致整套裝置都發(fā)生波動，引起報警泛濫，此時，報警的根原因難以判斷。異常工況出現(xiàn)時，現(xiàn)場的操作工人會根據(jù)自己的經(jīng)驗來進行操作，如果他們不知道裝置中所有變量之間的相關(guān)性，就可能出現(xiàn)判斷失誤和錯誤的決策，進而可能導(dǎo)致重大事故。為輔助操作人員應(yīng)對異常工況，有關(guān)異常工況管理的研究已經(jīng)有近40年的歷史，國際上有關(guān)研究人員提出了多種模型和算法, 這些方法可以分為3類：基于定性模型的方法、基于定量模型的方法和基于數(shù)據(jù)過程數(shù)據(jù)驅(qū)動的方法，而基于數(shù)據(jù)驅(qū)動的方法已經(jīng)展現(xiàn)出明顯的優(yōu)勢[1], 因為這類方法確實能夠幫助操作人員縮小甚至確定異常工況出現(xiàn)的原因[2]。

異常工況是指由于流程內(nèi)部或外部的擾動而導(dǎo)致系統(tǒng)的某些狀態(tài)參數(shù)超出了控制系統(tǒng)所能夠控制的正常范圍，進而需要操作人員進行手動干預(yù)的那些工況。擾動一般指系統(tǒng)引入一個未知或者不可控的輸入[3]。分析擾動的傳播路徑，可以確定異常工況出現(xiàn)的原因。而分析擾動的傳播路徑其核心在于尋找變量之間的因果關(guān)系[4]。對于工程師來說，通過因果關(guān)系，就能知道如何調(diào)節(jié)系統(tǒng)的變量得到想要的結(jié)果。因果關(guān)系可以從兩個角度進行考慮：①概率上的因果性；②時間上的因果性。從概率角度上來說，事件A的發(fā)生會提高事件B的發(fā)生概率，就意味著事件A和事件B具有因果關(guān)系，這種因果關(guān)系可以通過統(tǒng)計檢驗或者信息熵最低的方法來檢驗。從時間角度來看，如果事件A是事件B發(fā)生的原因，那么事件A一定發(fā)生在事件B之前，分析信號在頻域的相似性和時間上的順序，判斷因果關(guān)系。

目前有許多方法進行因果性的判斷。格蘭杰檢驗由Granger[5]在1969年提出，被廣泛應(yīng)用于經(jīng)濟領(lǐng)域。該方法假設(shè)有兩個時間序列{}和{}，如果同時利用序列{}和序列{}的過去數(shù)據(jù)來預(yù)測{}，要優(yōu)于單獨用{}過去的數(shù)據(jù)來預(yù)測{}，那么就意味著存在從{}到{}的因果關(guān)系，格蘭杰檢驗不適用于耦合的動態(tài)系統(tǒng)[6]?；ハ嚓P(guān)函數(shù)通過兩個信號的相似性和時間特征，來判斷信號之間的因果關(guān)系，互相關(guān)函數(shù)適用于線性系統(tǒng)[7]。Schreiber等[8-9]提出傳遞熵算法來判斷因果關(guān)系，從信息熵的概念出發(fā)，如果一個變量能夠減小另外一個變量的不確定性，那么就意味著存在從變量到變量的因果性，但是傳遞熵要求數(shù)據(jù)具有時間穩(wěn)定性，也就是均值和方差不隨時間變化[10]。2012年，Sugihara等[11]提出了一種收斂交叉映射（CCM）的方法，其基本思想是，如果變量對變量有影響，那么變量將含有變量的信息，通過測量和的重構(gòu)流型的相關(guān)性，來檢測因果性，不適用于化工動態(tài)系統(tǒng)，因為原有的CCM算法并沒有對因果關(guān)系的時間特性進行討論[12-13]。針對以上問題，本文提出了具有時間特征的收斂交叉映射算法和確定最佳嵌入維度的方法，在簡單的生態(tài)系統(tǒng)，因果檢測基準(zhǔn)系統(tǒng)和全混釜反應(yīng)器（CSTR）中驗證所提算法的優(yōu)越性。

1 收斂交叉映射算法的介紹

假設(shè)維空間上有維（<=）隨時間變化的流形，{}是流形投影于一維空間產(chǎn)生的序列，{}是流形投影于另一維空間產(chǎn)生相同長度的序列[14]。對于兩個長度為時間序列[]和[]，設(shè)重構(gòu)流形的維度為，采樣間隔為，那么在時刻重構(gòu)流形的坐標(biāo)為

由式(1)和式(2)得到重構(gòu)流形M={()}，M={()}，根據(jù)Takens嵌入定理[15]，流形M、M和是微分同胚的。定義[]M從流形M通過交叉映射式(3)～式(5)得到[]。

從流形M找到距離()最近的+1個點，其對應(yīng)序列[]上的[t]

d[(),(t)]代表著流形上()和(t)之間的歐氏距離，并計算和[]的相關(guān)系數(shù)。相關(guān)系數(shù)的計算公式如下

2 算法改進

2.1 最佳嵌入維度的選擇

流形的嵌入維度對于CCM的因果關(guān)系分析結(jié)果有著重要的影響。流形維度選擇太小，重構(gòu)的流形無法涵蓋所有的信息。流形維度選擇太大，計算量增大。本文提出了用赤池信息準(zhǔn)則來判斷最佳的嵌入維度[16-17]。赤池信息準(zhǔn)則是建立在熵的概念基礎(chǔ)上[18]，可以權(quán)衡所估計模型的復(fù)雜度和此模型擬合數(shù)據(jù)的優(yōu)良性[19-20]。模型復(fù)雜度定義為，定義損失函數(shù)，樣本數(shù)量為，赤池信息準(zhǔn)則定義如式(7)所示。

在本文，通過[]的自回歸模型來判斷最佳的嵌入維度，自回歸模型如式（8）所示。

損失函數(shù)=SSRauto，模型復(fù)雜度=，當(dāng)AIC()取到最小值，所對應(yīng)的即為最佳嵌入維度。

2.2 具有時間特性的CCM算法

為了使CCM具有時間特征，在原有的收斂交叉映射算法上，在流形M重構(gòu)過程中，加入代表因果關(guān)系時間特征的時滯，此時流形重構(gòu)的坐標(biāo)如式(9)所示。

由式(1)和式(7)得重構(gòu)的流形M={()}，={(-)}通過交叉映射式(3)～式(5)來估計，通過計算估計值和[]的殘差平方和來判斷估計效果，殘差平方和計算公式如式(10)所示。

SSR越小，意味著估計效果越好。SSR是一個關(guān)于時滯的函數(shù)，記SSR取得最小值所對應(yīng)的時滯為min，若min小于0，則意味著不存在從到的因果關(guān)系。當(dāng)min大于0，通過自回歸模型得到殘差平方和SSRauto作為檢測因果關(guān)系的閾值。當(dāng)min>0SSR小于SSRauto，認(rèn)為存在從[]到[]的因果關(guān)系。為表達因果關(guān)系強度，本文定義了如式(11)所示的影響強度，當(dāng)CT值越接近于1，說明從[]到[]的因果關(guān)系越強。

具有時間特征的收斂交叉映射算法流程如下：

（1）對進行自回歸，根據(jù)赤池信息準(zhǔn)則，確定最佳的嵌入維度，其殘差平方和記作SSRauto；

（2）根據(jù)嵌入維度，重構(gòu)影子流形M和；

（3）通過影子流形來估計，得到；

（5）如果min>0SSR小于SSRauto，則存在從到的因果關(guān)系，計算影響強度CT。

3 案例研究

為了驗證本文所提的算法的有效性，分別在簡單的生態(tài)系統(tǒng)、因果檢測基準(zhǔn)系統(tǒng)和全混釜反應(yīng)器（CSTR）等10個系統(tǒng)進行應(yīng)用驗證。

3.1 簡單的生態(tài)系統(tǒng)

在Sugihara等[11]提出CCM算法的文章中，有一個關(guān)于沙丁魚和鳳尾魚的二元生態(tài)系統(tǒng)。其中代表沙丁魚的數(shù)量，代表著鳳尾魚的數(shù)量，其模型可以用式(12)和式(13)表示

(+1)=()[r-rX()-,y()] (12)

(+1)=()[r-rY()] (13)

其中，r=3.8,,y=0.375,r=3.6，初值為0.4，初值為0.2。下面對嵌入維度對結(jié)果影響進行研究。

從模型公式[式(12)和式(13)]，可以看出存在從到的因果關(guān)系。在Sugihara提出CCM算法的文章中，并沒有對嵌入維度選擇進行討論。圖1顯示了不同嵌入維度對于因果關(guān)系檢測影響，其中虛線代表著從到的因果關(guān)系隨長度的變化曲線。當(dāng)嵌入維度為1，由于重構(gòu)的流形并未包含所有的信息，所以無法檢測從到的因果關(guān)系。當(dāng)嵌入維度為2，可以發(fā)現(xiàn)存在從到的因果關(guān)系，當(dāng)嵌入維度為3，計算結(jié)果與嵌入維度為2計算結(jié)果相似，因此最佳嵌入維度為2。對于進行自回歸，得到AIC()的分布圖，如圖2所示，在=2時，AIC()取得最小值，這與上述結(jié)果相符。

3.2 因果檢測基準(zhǔn)系統(tǒng)

對如圖3所示的8個因果檢測基準(zhǔn)系統(tǒng)，用上述算法進行測試。8個基準(zhǔn)系統(tǒng)依次是：(Ⅰ)一階系統(tǒng)；(Ⅱ)一階平方系統(tǒng)；(Ⅲ)具有反饋回路的一階系統(tǒng)；(Ⅳ)引入噪聲的一階系統(tǒng)；(Ⅴ)引入階躍的一階系統(tǒng)；(Ⅵ)饋通一階系統(tǒng)；(Ⅶ)引入死區(qū)時間的一階系統(tǒng)；(Ⅷ)二階系統(tǒng)。其中輸入信號()為滿足(0,1)分布的白噪聲，采樣間隔為0.1 s。所有的標(biāo)準(zhǔn)系統(tǒng)均采用放大因子為1和時間常數(shù)為0.5 s的傳遞函數(shù)。

圖4為用原有CCM算法的計算結(jié)果，其中實線代表著從到y(tǒng)的因果關(guān)系判斷，可以看出除了基準(zhǔn)系統(tǒng)(Ⅵ)以外，其他基準(zhǔn)系統(tǒng)的實線都在橫軸附近，也就意味著無法檢測出從到的因果關(guān)系。

圖5顯示對于8個簡單系統(tǒng)收斂交叉映射計算結(jié)果，其中(Ⅰ)、(Ⅱ)、(Ⅲ)、(Ⅴ)、(Ⅵ)和(Ⅷ)均在時滯大于0出現(xiàn)低谷。

表1為原有的CCM算法和具有時間特性的CCM算法的對比結(jié)果，其中對號代表著成功檢測出因果關(guān)系。原有CCM很難檢測出化工動態(tài)系統(tǒng)的因果關(guān)系，或者檢測出一個相反的因果關(guān)系。具有時間特性CCM能夠檢測出8個系統(tǒng)中的6個系統(tǒng)的因果關(guān)系，體現(xiàn)具有時間特性的CCM對于因果關(guān)系判斷有著較好的效果。

表1 因果檢驗計算結(jié)果

3.3 全混釜反應(yīng)器（CSTR）

圖6(a)為CSTR的示意圖，圖6(b)為機理推導(dǎo)出的因果關(guān)系。CSTR中反應(yīng)物A反應(yīng)先生成中間產(chǎn)物B，再反應(yīng)生成最終產(chǎn)物C。

反應(yīng)模型如式(14)～式(16)所示。

其中，代表CSTR的體積，代表進入CSTR的體積流量，1和2分別代表著反應(yīng)(1)和反應(yīng)(2)的活化能，A()、B()和C()分別代表著反應(yīng)器出口物料中物質(zhì)A、物質(zhì)B和物質(zhì)C的摩爾濃度，in為全混釜進料中物質(zhì)A的濃度，in=1 mol·L-1，在in加入高斯白噪聲(0,0.1)。fl為反應(yīng)溫度，fl=350 K，在fl加入高斯白噪聲(0,32)。其他參數(shù)如表2所示。

表2 CSTR參數(shù)

用Matlab中的simulink進行模擬，運行時間1000 s，采樣間隔為0.1 s。

表3 具有時間特征的收斂交叉映射計算強度結(jié)果

表3為具有時間特征的收斂交叉映射計算出來的變量之間的因果強度，其中表中的值代表著行變量對列變量的影響。從表3的結(jié)果看出，具有時間特征的收斂交叉映射都可以檢測出大部分的因果關(guān)系，但是從B()到C()因果關(guān)系沒有檢測出來。

4 結(jié) 論

本文提出了具有時間特征的收斂交叉映射算法，并在簡單生態(tài)系統(tǒng)，因果檢測基準(zhǔn)系統(tǒng)和CSTR模型進行驗證，得到以下結(jié)論。

（1）嵌入維度對于因果關(guān)系判斷有著重要的影響，通過赤池信息準(zhǔn)則能夠有效選擇合適的嵌入維度，優(yōu)化計算過程。

（2）對于化工動態(tài)系統(tǒng)來說，時間特性是判斷因果關(guān)系的重要依據(jù)，通過在構(gòu)建流形過程中，加入代表著時間特性的時滯，能夠有效地判斷化工動態(tài)系統(tǒng)的因果關(guān)系。

具有時間特征的CCM并沒有將系統(tǒng)中所有的因果關(guān)系判斷出來，未來還需要對其進行改進。

[1] SHU Y, MING L, CHENG F,. Abnormal situation management: challenges and opportunities in the big data era[J]. Computers & Chemical Engineering, 2016, 91: 104-113.

[2] KUHNERT C. Data-driven Methods for Fault Localization in Process Technology[M]. KIT Scientific Publishing, 2013: 12-17.

[3] IERMANN R. Fault-diagnosis Systems: An Introduction from Fault Detection to Fault Tolerance[M]. Springer Science & Business Media, 2006: 34-40.

[4] KUHNERT C, BEYERER J. Data-driven methods for the detection of causal structures in process technology[J]. Machines, 2014, 2(4): 255-274.

[5] GRANGER C W. Investigating causal relations by econometric models and cross-spectral methods[J]. Econometrica: Journal of the Econometric Society, 1969, 37(3): 424-438.

[6] MCCRACKEN J M, WEIGEL R S. Convergent cross-mapping and pairwise asymmetric inference[J]. Physical Review E, 2014, 90(6): 062903.

[7] BILLINGS S A. Nonlinear System Identification: NARMAX Methods in the Time, Frequency, and Spatio-Temporal Domains[M]. John Wiley & Sons, 2013: 78-85.

[8] SCHREIBER T. Measuring information transfer[J]. Physical Review Letters, 2000, 85(2): 461.

[9] SHU Y D, ZHAO J S. Data driven causal inference based on a modified transfer entropy[J]. Computers & Chemical Engineering, 2013, 57(10): 173-180

[10] BAUER M, COX J W, CAVENESS M H,. Finding the direction of disturbance propagation in a chemical process using transfer entropy[J]. Control Systems Technology, IEEE Transactions on, 2007, 15(1): 12-21.

[11] SUGIHARA G, MAY R, YE H,. Detecting causality in complex ecosystems[J]. Science, 2012, 338(6106): 496-500.

[12] MA H, AIHARA K, CHEN L. Detecting causality from nonlinear dynamics with short-term time series[J]. Scientific Reports, 2014, 4: 7464.

[13] CLARK A T, YE H, ISBELL F,. Spatial convergent cross mapping to detect causal relationships from short time series[J]. Ecology, 2015, 96(5): 1174-1181.

[14] WEI W W S. Time Series Analysis[M]. Addison-Wesley, 1994: 77-82.

[15] TAKENS F. Detecting strange attractors in turbulence[M]//Dynamical Systems and Turbulence, Warwick 1980. Springer Berlin Heidelberg, 1981: 366-381.

[16] AKAIKE H. Akaike’s information criterion[M] // International Encyclopedia of Statistical Science. Springer Berlin Heidelberg, 2011: 25.

[17] WISMüLLER A, ABIDIN A Z, D'SOUZA A M,. Nonlinear functional connectivity network recovery in the human brain with mutual connectivity analysis (MCA): convergent cross-mapping and non-metric clustering[C]//SPIE Medical Imaging. International Society for Optics and Photonics, 2015: 94170M-94170M-9.

[18] YAMAOKA K, NAKAGAWA T, UNO T. Application of Akaike’s information criterion (AIC) in the evaluation of linear pharmacokinetic equations[J]. Journal of pharmacokinetics and biopharmaceutics, 1978, 6(2): 165-175.

[19] BROOKS D G. Akaike information criterion statistics[J]. Technometrics, 1989, 31(2): 270-271.

[20] BOZDOGAN H. Model selection and Akaike’s information criterion (AIC): the general theory and its analytical extensions[J]. Psychometrika, 1987, 52(3): 345-370.

Convergent cross mapping method in analysis of disturbances in chemical processes

CHENG Feifan, ZHAO Jinsong

(Department of Chemical Engineering, Tsinghua University, Beijing 100084, China)

When there is a fault in the chemical plant, the operators may not be able to find the root cause of the fault. Data-driven method can be a great help to reduce the uncertainty of root cause, even to locate the root cause. By analysis of chemical process data, the disturbance propagation way can be detected, which can help to locate the root cause. In this article, convergent cross mapping (CCM) with the characteristic of time is proposed to detect the causality in dynamic chemical process. Furthermore, the usage of Akaike information criterion is proposed to determine the most proper embedding dimension. To prove the effectiveness of the method, the new method is applied to the ecosystem examples, causality testing benchmark model and CSTR model. Comparing the result calculated by the original CCM, the effectiveness of the new method is found.

disturbance analysis；causality analysis；convergent cross mapping；characteristic of time；safety；systems engineering；simulation

date: 2016-09-21.

Prof. ZHAO Jinsong, jinsongzhao@tsinghua. edu. cn

10.11949/j.issn.0438-1157.20161318

TQ 015.9

A

0438—1157（2016）12—5082—07

國家自然科學(xué)基金項目（61433001）。

supported by the National Natural Science Foundation of China (61433001).

2016-09-21收到初稿，2016-09-26收到修改稿。

聯(lián)系人：趙勁松。第一作者：程非凡（1993—），男，博士研究生。