基于收敛交叉映射的化工过程扰动因果分析方法

程非凡，赵劲松



基于收敛交叉映射的化工过程扰动因果分析方法

程非凡，赵劲松

（清华大学化学工程系，北京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—），男，博士研究生。