程非凡，趙勁松

?

基于收斂交叉映射的化工過(guò)程擾動(dòng)因果分析方法

程非凡，趙勁松

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

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

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

引 言

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

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

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

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

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

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

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

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

2 算法改進(jìn)

2.1 最佳嵌入維度的選擇

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

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

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

2.2 具有時(shí)間特性的CCM算法

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

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

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

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

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

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

（3）通過(guò)影子流形來(lái)估計(jì)，得到；

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

3 案例研究

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

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

在Sugihara等[11]提出CCM算法的文章中，有一個(gè)關(guān)于沙丁魚(yú)和鳳尾魚(yú)的二元生態(tài)系統(tǒng)。其中代表沙丁魚(yú)的數(shù)量，代表著鳳尾魚(yú)的數(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。下面對(duì)嵌入維度對(duì)結(jié)果影響進(jìn)行研究。

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

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

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

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

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

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

表1 因果檢驗(yàn)計(jì)算結(jié)果

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

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

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

其中，代表CSTR的體積，代表進(jìn)入CSTR的體積流量，1和2分別代表著反應(yīng)(1)和反應(yīng)(2)的活化能，A()、B()和C()分別代表著反應(yīng)器出口物料中物質(zhì)A、物質(zhì)B和物質(zhì)C的摩爾濃度，in為全混釜進(jìn)料中物質(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進(jìn)行模擬，運(yùn)行時(shí)間1000 s，采樣間隔為0.1 s。

表3 具有時(shí)間特征的收斂交叉映射計(jì)算強(qiáng)度結(jié)果

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

4 結(jié) 論

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

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

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

具有時(shí)間特征的CCM并沒(méi)有將系統(tǒng)中所有的因果關(guān)系判斷出來(lái)，未來(lái)還需要對(duì)其進(jìn)行改進(jì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

國(guó)家自然科學(xué)基金項(xiàng)目（61433001）。

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

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

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