一种基于时序数据分析的连续搅拌反应釜运行状态监测方法

本发明涉及一种化工过程监测方法，特别涉及一种基于时序数据分析的连续搅拌反应釜运行状态监测方法。

化工过程实时采样数据是化工“大数据”建设与应用的基础，利用采样数据来监测化工过程的运行状态已成为化工智能系统的重要组成部分。数据驱动的化工过程监测经过多年的研究发展，出现了许多以主成分分析(Principal Component Analysis，缩写：PCA)与独立成分分析(Independent Component Analysis，缩写：ICA)算法为主的过程监测方法。这些主流的过程监测方法实施的核心主要关注于数据潜在特征的挖掘。换句话说，所建立的数据驱动模型都是旨在提取数据特征。近十几年来，无论是学术界还是工业界，都投入了大量的人力与物力研究数据驱动的故障检测方法与技术。连续搅拌反应釜(ContinuousStirred Tank Reactor，简称：CSTR)是化工生产中进行各种物理变化和化学反应广泛使用的设备，在反应装置中占有重要地位。由于CSTR设备在实际生产过程中的广泛应用和重要性，CSTR的运行状态监测一直受到化工过程监测领域专业技术人员的关注。

由于先进测量仪表的广泛应用，CSTR运行过程中的采样数据采样间隔短，不可避免地在时间序列上存在自相关性。因此，通过监测采样数据时间序列的异常变化情况是可以反映出CSTR的运行是否进入异常状态。在现有文献与专利材料中，对化工过程对象实施动态过程监测大多依赖于为各个采样数据引入延时测量数据，即将多个在采样时间上连续的样本数据当成一个样本，然后再实施建模与监测。这种方法技术的典型代表就是动态PCA与动态ICA，都是将时间序列上的自相关性与交叉相关性混合在一起同时提取。所提取的特征不具备对时间序列相关与交叉相关特征的可解释性。最近，也有研究工作提出通过最大化潜在特征成分间协方差的方式，引导采样数据的潜在特征的挖掘，而不依赖于使用增广向量或矩阵，从而实现对时序相关特征的挖掘。典型的代表方式主要是基于动态潜变量(Dynamic Latent Variable，缩写：DLV)模型的过程监测方法。

然而，协方差虽然能在一定程度上体现相关性，但是数据间的共线性问题同样会导致协方差指标的最大化。因此，提取采样数据中的时间序列相关特征应该考虑使用典型相关系数指标。众所周知，两个随机变量间的典型相关系数是在区间[-1，1]之间取值，正负符号表示的是方向，而数值则表示两个随机变量间的相关性大小。因此，典型相关系数接近于-1的两个随机变量同样存在较大的典型相关性。此外，由于时间序列采样数据通常考虑的是多个连续时间上的样本数据，为了提取时间序列相关特征，需要充分考虑各个时间序列样本数据间的相关性问题。因此，挖掘采样数据的时间序列相关特征不仅需要使用典型相关系数指标，而且还需要考虑多个采样时间上连续的多个样本间的时序性。

本发明所要解决的主要技术问题是：推理出一种时序数据分析的算法，通过监测连续搅拌反应釜的实时采样数据，来实现对连续搅拌反应釜的运行状态监测。具体来讲，本发明方法首先根据时序数据分析算法两个要求：其一，典型相关系数量化相关性大小；其二，样本数据的采样时序特性，来设计相应的优化目标；其次，本发明方法提出一种迭代算法来优化求解各个载荷向量；再次，本发明方法使用自回归模型描述时序相关特征的动态性；最后，基于此模型实施对连续搅拌反应釜运行状态的在线实时监测。

本发明方法解决上述问题所采用的技术方案为：一种基于时序数据分析的连续搅拌反应釜运行状态监测方法，包括以下步骤：

步骤(1)：按照采样时间先后顺序，采集连续搅拌反应釜设备在正常运行状态下的n个数据向量x1，x2，...，xn，从而组成训练数据矩阵X＝[x1，x2，...，xn]∈Rm×n，并计算均值向量μ＝(x1+x2+...+xn)/n以及标准差向量δ∈Rm×1：

其中，⊙表示向量(xi-μ)与(xi-μ)中的对应元素相乘、m为连续搅拌反应釜设备所涉及的测量变量的总个数、i＝1，2，...，n、R为实数集、Rm×n表示m×n维的实数矩阵、xi∈Rm×1表示第i个数据向量。

值得注意的是，数据向量xi中m个元素分别由对应的测量仪表测量得到，分别包括反应器进料流量、反应器压力、反应器液位、反应器温度、反应器进料阀门开度、反应器冷凝水流量，以及冷凝器冷却水流量，这7个测量数据。因此，m＝7。

此外，由于各个测量变量的变化范围不可能一致，也就导致各个测量变量之间存在量纲的差异影响。因此，需要使用标准化处理的方式，将各个测量变量的采样数据皆变换成均值为0，标准差为1的数据，如步骤(2)所示。

步骤(2)：根据公式对样本数据x1，x2，...，xn分别实施标准化处理，从而得到矩阵其中，表示向量(xi-μ)与标准差向量δ中的对应元素相除，为标准化后的数据向量。

步骤(3)：设置时间序列相关阶数为D(一般可取D＝3或4)后，根据如下所示公式构造时间序列矩阵X1，X2，...，XD：

上式中，d＝1，2，...，D，N＝n-D+1。

步骤(2)中的标准化为对采样数据实施预处理的阶段，是几乎所有数据驱动的过程监测方法都必须实施的基本过程；接下来，需要以矩阵为训练数据集对其实施时序相关特征分析算法。

首先，需要量化出时序相关特征分析算法的目标函数。考虑到该算法有两个量化要求：其一，典型相关系数的大小量化相关性；其二，样本数据的采样时序特性，因此设计出如下所示的优化目标函数：

上式中，c＝1，2，...，D-1、w1，w2，...，wD分别为对应于间序列矩阵X1，X2，...，XD的载荷向量、d＝1，2，...，D、上标号T表示矩阵或向量的转置、arg max表示最大化目标函数、s.t.为单词Subject To的首字母缩写，表示约束条件的意思。

从上式③中可以看出，该优化问题中的目标函数对得分向量与之间的典型相关系数实施平方和操作；目标函数的设计只考虑到了典型相关系数的大小而不会受其正负方向的影响，且是按照时间先后的时序特性计算相关性大小的。

其次，上式③中带约束条件的优化问题求解需要使用拉格朗日乘子法。在此之前，若是令则上式③可等价转换成如下所示的问题：

上式④中，vd∈Rm×1表示伪载荷向量、d＝1，2，...，D。

然后，针对上式④所定义的最大化问题，构造如下所示的拉格朗日函数L：

上式⑤中，λ1，λ2，...，λD分别表示拉格朗日乘子。

按照拉格朗日乘子法的求解思路，需将拉格朗日函数L分别对v1，v2，...，vD求偏导数，可得如下结果：

再分别令各个偏导数等于零，即可得到如下所示的等式关系：

若是在上式⑦中各个等式两边依次左乘v1T，v2T，...，vDT，并将的约束条件考虑进来，即可得到如下所示的等式关系：

由此可见，λD＝λ1+λ2+...+λD-1即等价于上式④中需要最大化的目标函数。

此外，将上式⑦中前D-1个等式代入上式⑦中最后一个等式中，可得到如下所示的特征值问题：

因此，最大化λD等价于求解上式⑨中最大的特征值问题，并且vD为最大特征值所对应的特征向量。若vD已知，根据上式⑦可以发现vc的方向是与向量的方向是一致的，因此可直接根据式⑦中的等式计算得到v1，v2，...，vD-1。

如此一来，利用时间序列相关特征分析算法优化求解载荷向量的实施步骤就简化成了简单的特征值与特征向量的求解问题。而且，由于矩阵是对称的，其对应的特征向量则是相互正交的。

步骤(4)：求解矩阵Θ最大的k个特征值所对应的特征向量此步骤要求所有特征向量的长度都为1，具体的实施过程如下所示：

步骤(4.1)：根据公式计算矩阵Θ∈Rm×m，并设置j＝1。

步骤(4.2)：初始化伪载荷向量vD∈Rm×1为任意实数向量后，根据公式vD＝vD/||vD||更新vD，其中||vD||表示计算vD的长度。

步骤(4.3)：若j＜2，则依次根据公式vD＝ΘvD与vD＝vD/||vD||更新vD；若j≥2，则依次根据公式vD＝(Im-VDVDT)ΘvD与vD＝vD/||vD||更新vD，其中Im为m×m维的单位矩阵，伪载荷矩阵

步骤(4.4)：判断vD是否收敛，收敛的标准为vD中各元素不再发生变化，若否，则返回步骤(4.3)；若是，则设置第j个特征向量为后执行步骤(4.5)。

步骤(4.5)：依次根据公式与vc＝vc/||vc||计算对应于第c个时间序列矩阵Xc的第j个伪载荷向量其中c＝1，2，...，D-1。

步骤(4.6)：根据公式计算第j个特征值后，判断是否小于0.1，若否，则设置j＝j+1后返回步骤(4.2)；若是，则设置k＝j-1后，得到各个时间序列矩阵所对应的伪载荷矩阵其中d＝1，2，...，D。

步骤(5)：根据公式分别计算得到载荷矩阵W1，W2，...，WD后，再根据公式计算得到得分矩阵S1，S2，...，SD，其中矩阵

步骤(6)：分别根据公式FD＝SD-ZDθD与公式计算残差矩阵ED与FD，其中ZD＝[S1，S2，...，SD-1]、

步骤(7)：根据公式与分别计算得到监测指标向量与Q，并分别将与Q中的元素最大值记录为与Qmax，其中diag{ }表示将大括号中矩阵对角线的元素变成列向量的操作运算，矩阵

步骤(8)：根据公式计算综合监测指标向量φ后，再利用核密度估计(Kernel Density Estimation)法确定出在置信限α＝99％条件下的控制上限φlim。

上述步骤(1)至步骤(8)为本发明方法的离线建模阶段，在离线建模阶段完成后，即可利用实时测量的样本数据实施对连续搅拌反应釜的运行状态监测。此外，由于连续搅拌反应釜的在线采样数据有可能存在个别测量变量的数据缺失，因此还需对在线采样数据做数据校正处理。

步骤(9)：在线采集连续搅拌反应釜最新采样时刻的数据向量xt∈Rm×1，判断数据向量xt中是否存在缺失数据，若是，则将均值向量μ中对应的元素补充至数据向量xt中，下标号t表示最新采样时刻。

步骤(10)：根据公式对其实施标准化处理得到向量

步骤(11)：根据公式计算得到得分向量sD，并根据公式计算将t-D+1至t-1采样时刻数据向量所对应的得分向量s1，s2，...，sD-1，其中c＝1，2，...，D-1。

步骤(12)：根据公式ξD＝sD-zDθD与分别计算残差向量ξd与ed。

步骤(13)：根据公式与计算t采样时刻的监测指标与Qt后，再根据公式计算综合监测指标

步骤(14)：判断是否满足条件：φt≤φlim？若是，则当前采样时刻连续搅拌反应釜的运行状态正常，返回步骤(9)继续对下一新时刻的采样数据实施监测；若否，则当前采样时刻连续搅拌反应釜的运行状态出现异常，触发故障警报并返回步骤(9)继续对下一新时刻的采样数据实施监测。

与传统数据驱动化工过程监测方法(尤其是针对时序相关性问题的动态过程监测方法)相比，本发明方法的优势在于：

首先，本发明方法所涉及的时间序列相关特征分析算法能够有效地提取采样数据体现在时间序列上典型相关的潜在特征成分，对采样数据时序相关动态特征挖掘相比于以方差或协方差为量化指标的其他传统方法更充分。其次，本发明方法通过自回归模型进一步描述时序相关特征的时序特性，并通过监测误差的变化来实现对连续搅拌反应釜过程采样数据时序异常变化的监测。最后，通过具体的实施案例，对比验证了本发明方法在监测连续搅拌反应釜运行状态上的优越性与有效性。

图1为本发明方法的实施流程示意图。

图2为连续搅拌反应釜的流程图及其相应的测量仪表。

图3为本发明方法与其他传统方法在监测连续搅拌反应釜运行状态时的对比图。

下面结合附图和具体实施方式对本发明进行详细说明。

如图1所示，本发明公开了一种基于时序数据分析的连续搅拌反应釜运行状态监测方法，下面结合一个具体应用实例来说明本发明方法的具体实施方式。

如图2所示，连续搅拌反应釜的流程图及其相应的测量仪表。本实施案例中的应用对象是一个涉及放热反应过程的CSTR设备，因此该CSTR设备需配备一个冷凝器对反应物的出口温度进行降温处理。从图2中可以发现，该连续搅拌反应釜所涉及的测量数据包括7个，分别是：反应器进料流量、反应器压力、反应器液位、反应器温度、反应器进料阀门开度、反应器冷凝水流量，以及冷凝器冷却水流量。

首先，利用CSTR过程对象在正常工况下采样的n＝960个样本数据实施本发明方法的离线建模阶段，包括以下步骤：

步骤(1)：按照采样时间先后顺序，采集连续搅拌反应釜设备在正常运行状态下的n＝960个数据向量x1，x2，...，x960组成训练数据矩阵X＝[x1，x2，...，xn]∈R7×960，并计算均值向量μ＝(x1+x2+...+x960)/960以及标准差向量δ∈R7×1。

步骤(2)：根据公式对样本数据x1，x2，...，x960分别实施标准化处理，从而得到矩阵

步骤(3)：设置时间序列相关阶数为D＝3，根据上述公式②构造时间序列矩阵X1，X2，X3。

步骤(4)：求解矩阵最大的k＝3个特征值所对应的特征向量并按照上述步骤(4.1)至步骤(4.6)得到各个时间序列矩阵所对应的伪载荷矩阵

步骤(5)：根据公式分别计算得到载荷矩阵W1，W2，W3后，再根据公式计算得到得分矩阵S1，S2，S3，其中矩阵

步骤(6)：分别根据公式F3＝S3-Z3θ3与公式计算残差矩阵E3与F3，其中Z3＝[S1，S2]；

步骤(7)：根据公式与分别计算得到监测指标向量与Q，并分别将与Q中的元素最大值记录为与Qmax，其中diag{ }表示将大括号中矩阵对角线的元素变成列向量的操作运算，矩阵

步骤(8)：根据公式计算综合监测指标向量φ后，再利用核密度估计法确定出在置信限α＝99％条件下的控制上限φlim。

离线建模阶段至此完成，接下来进入在线动态过程监测阶段，需要实时利用CSTR化工过程对象的在线采样数据。在本实施案例中，CSTR的在线采样数据的前160个数据采集自正常运行状态，故障工况从第161个样本点开始TE过程才进入故障工况。

步骤(9)：在线采集CSTR过程对象最新采样时刻的数据向量xt∈R7×1，判断数据向量xt中是否存在缺失数据，若是，则将均值向量μ中相应位置的元素补充至数据向量xt中，下标号t表示最新采样时刻。

步骤(10)：根据公式对其实施标准化处理得到向量

步骤(11)：根据公式计算得到得分向量sD，并根据公式计算将t-D+1至t-1采样时刻数据向量所对应的得分向量s1，s2，其中c＝1，2。

步骤(12)：根据公式ξD＝sD-zDθD与分别计算残差向量ξd与ed。

步骤(13)：根据公式与计算t采样时刻的监测指标与Qt后，再根据公式计算综合监测指标

步骤(14)：判断是否满足条件：φt≤φlim？若是，则当前采样时刻运行正常，返回步骤(9)继续对下一新时刻的采样数据实施监测；若否，则当前采样时刻运行出现异常，触发故障警报并返回步骤(9)继续对下一新时刻的采样数据实施监测。

如图3所示，本发明方法与传统DPCA与DLV方法的在监测故障工况数据的监测图。从图3中对比可以很明显地发现，本发明方法在故障检测成功率明显优越于其他动态过程监测方法。因此，可以说本发明方法具有更可靠的过程监测性能。

上述实施案例只用来解释说明本发明的具体实施，而不是对本发明进行限制。在本发明的精神和权利要求的保护范围内，对本发明做出的任何修改，都落入本发明的保护范围。