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

程非凡，趙勁松

?

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

程非凡，趙勁松

（清華大學化學工程系，北京100084）

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

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

引 言

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

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

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

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

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

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

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

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

2 算法改進

2.1 最佳嵌入維度的選擇

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

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

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

2.2 具有時間特性的CCM算法

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

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

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

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

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

（2）根據嵌入維度，重構影子流形M和；

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

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

3 案例研究

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

3.1 簡單的生態系統

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

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

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

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

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

3.2 因果檢測基準系統

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

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

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

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

表1 因果檢驗計算結果

3.3 全混釜反應器（CSTR）

圖6(a)為CSTR的示意圖，圖6(b)為機理推導出的因果關系。CSTR中反應物A反應先生成中間產物B，再反應生成最終產物C。

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

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

表2 CSTR參數

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

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

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

4 結 論

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

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

（2）對于化工動態系統來說，時間特性是判斷因果關系的重要依據，通過在構建流形過程中，加入代表著時間特性的時滯，能夠有效地判斷化工動態系統的因果關系。

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

[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

國家自然科學基金項目（61433001）。

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

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

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