CN108427398A - 一种基于分散式ar-pls模型的动态过程监测方法 - Google Patents

Abstract

本发明公开一种基于分散式AR‑PLS模型的动态过程监测方法，本发明方法首先将各测量变量引入在采样时间上的多个延时测量值，然后针对每个测量变量为之建立该变量与其他测量变量及延时测量值之间的AR‑PLS模型，最后利用模型的预测误差实施对动态过程的在线监测。相比于传统动态过程监测方法，本发明方法首先在建模阶段体现出了分散式建模的特点，因此具备分散式多模型的优势。其次，本发明方法只利用AR模型的预测误差作为被监测对象，而误差是AR模型中剔除了自相关性与交叉相关性后的结果，通过这种思路巧妙地避免了动态过程监测中自相关性问题。可以说，本发明方法是一种更为优选的动态过程监测方法。

Description

一种基于分散式AR-PLS模型的动态过程监测方法

技术领域

本发明涉及一种数据驱动的故障检测方法，尤其涉及一种基于分散式AR-PLS模型的动态过程监测方法。

背景技术

过程监测的目的在于及时而准确地发现故障，这对于保证安全生产与维持产品质量稳定是具有重要意义的。当前，过程监测的主流实施技术手段是数据驱动的方法，这主要得益于现代化工过程的大型化建设以及先进仪表与计算机技术的广泛应用，生产过程可以采集海量的数据。采样数据容易获取而机理模型难以获取的发展方式使得传统基于机理模型的故障检测方法逐渐没落。相比之下，数据驱动的故障检测方法不需要机理模型只需要采样数据，比较适合于现代工业过程运行状态的监测。一般来讲，特征挖掘是数据驱动的故障检测方法的基本出发点，相应故障检测的性能直接取决于对正常工况下的采样数据是否进行了有效的特征挖掘。实施在线故障检测时，通常是将在线采样数据经投影变换得到相应的潜在特征，然后依据潜在特征的统计量指标大小判断是否为故障样本。

从本质上讲，数据驱动的故障检测方法与基于机理模型的故障检测方法是存在显著差异的，前者旨在挖掘潜在特征而后者旨在生成误差。发展至今，数据驱动的故障检测研究领域涌现出了很多特征挖掘算法以及各式各样的建模思路。就特征挖掘算法而言，主成分分析(Principal Component Analysis，PCA)、独立成分分析(Independent ComponentAnalysis，ICA)、偏最小二乘(Partial Least Squares，PLS)、以及它们的非线性扩展形式是被广泛研究的算法。就建模思路而言，对大规模过程对象可采用多变量块的分散式建模策略，对多工况过程对象可使用多阶段建模等等。可以说，数据驱动的故障检测方法已如雨后春笋般蓬勃发展。

对现代工业过程对象而言，过程数据的自相关性(或称动态性)是一个很常见的问题。最经典的处理方式莫过于采用增广型矩阵的建模方式，通过为各个测量变量引入延时测量值，在建立特征挖掘模型从而将自相关性考虑进来。可以说，动态过程监测方法的有效性除取决于后续特征挖掘算法外，还取决于前面如何描述测量变量的自相关特性。现有专利与科研文献中，有部分研究人员对变量的自相关特征进行更为细致的描述，有采用自回归(Auto Regression，AR)模型过滤自相关特性，也有对每个变量的自相关性实施细致分析以选择体现再不同采样时刻上的数据。一般而言，自相关性主要体现在当前时刻的数据可以由之前采样时刻的数据预测出来。但是每个测量变量的自相关性一般是不一致的，同一而论的话可能会弱化相应动态过程监测方法的故障检测效果。因此，针对动态过程监测的研究还需进一步深入。

发明内容

本发明所要解决的主要技术问题是：如何针对每个测量变量利用PLS算法建立其各自的AR模型，然后对AR模型的预测误差作为被监测对象实施动态过程监测。具体来讲，本发明方法首先将各测量变量引入在采样时间上的多个延时测量值，然后针对每个测量变量为之建立该变量与其他测量变量及延时测量值之间的AR-PLS模型，最后利用模型的预测误差实施对动态过程的在线监测。

本发明解决上述技术问题所采用的技术方案为：一种基于分散式AR-PLS模型的动态过程监测方法，包括以下步骤：

(1)采集生产过程正常运行状态下的样本，组成训练数据集X∈R^N×m，并按照如下形式构造增广型矩阵X_a∈R^{(N-d)×m(d+1)}：

其中，N为训练样本数，m为过程测量变量数，R为实数集，R^N×m表示N×m维的实数矩阵，x_n∈R^1×m为第n个采样时刻的样本数据，下标号n＝1，2，…，N，d为引入的延时测量值的个数(一般取d＝1或2)。

(2)对矩阵X_a中每列实施标准化处理，得到均值为0，标准差为1的新数据矩阵后初始化i＝1，并将表示成如下形式：

其中，上标号1，2，…，m为测量变量标号，下标号0，1，…，d为当前时刻(即0)采样值及其d个延时测量值(即1，2，…，d)的标号，表示延时d个时刻测量到的m个N-d维的列向量。

(3)将矩阵中第i列数据取出作为AR模型的输出，矩阵中剩余的列则作为AR模型的输入矩阵再利用PLS算法训练得到与输出之间的AR模型，具体的实施过程如下所示：

①置k＝1与后，初始化向量

②根据公式w_k＝Z_iu_k/||Z_iu_k||、s_k＝Z_iw_k、分别计算得到系数向量w_k、得分向量s_k、和系数q_k，其中||Z_iu_k||表示计算向量Z_iu_k的长度；

③根据公式计算向量u_new；

④判断是否满足条件||u_k-u_new||＜10^-6？若否，则置u_k＝u_new后返回步骤②；若是，则执行⑤；

⑤根据公式p_k＝Z_i ^Ts_k/(s_k ^Ts_k)计算得到第k个投影向量p_k∈R^{[m(d+1)-1]×1}，并保留向量p_k、向量w_k、和系数q_k；

⑥判断矩阵Y_k＝s_kp_k ^T中的最大元素是否大于0.01？若是，根据公式Z_i＝Z_i-s_kp_k ^T更新矩阵Z_i后执行步骤⑦；若否，则得到投影矩阵P_i＝[p₁，p₂，…，p_k]、系数矩阵W_i＝[w₁，w₂，…，w_k]、和列向量Q＝[q₁，q₂，…，q_k]^T，并执行步骤⑧；

⑦判断k＜m(d+1)-1？若是，则置k＝k+1后，返回步骤②；若否，则得到最终的投影矩阵P_i＝[p₁，p₂，…，p_k]与系数矩阵W_i＝[w₁，w₂，…，w_k]；

⑧计算回归向量那么输入与输出x_i之间的AR模型为：其中，e_i表示第i个变量的预测误差。

(4)判断是否满足条件i＜m？若是，则置i＝i+1后返回步骤(3)；若否，则将得到的估计误差向量组成矩阵E＝[e₁，e₂，…，e_m]∈R^(N-d)×m后继续执行下一步骤(5)。

(5)将估计误差矩阵E作为新的训练数据矩阵，对其中的每一列实施标准化处理，得到均值为0，标准差为1的新数据矩阵

(6)对矩阵实施奇异值分解(Singular Value Decomposition，SVD)，即：具体的实施过程如下所示：

①设置r＝1与后，初始化列向量t_r为矩阵中的第一列；

②根据公式v_r＝F^Tt_r/(t_r ^Tt_r)计算得到向量v_r；

③根据公式t_new＝Fv_r/(v_r ^Tv_r)计算向量t_new；

④判断是否满足条件||t_new-t_r||＜10^-6？若是，则执行下一步骤⑤；若否，则置t_r＝t_new后返回步骤②；

⑤根据公式与μ_r＝t_rλ_r ^-1分别计算得到第r个奇异值λ_r与向量μ_r，并根据公式F＝F-t_rv_r ^T更新F；

⑥判断是否满足条件λ_r≤10^-3？若否，则设置r＝r+1与向量t_r为矩阵F中的第一列后返回步骤②；若是，则将所有得到的奇异值λ₁，λ₂，…，λ_r组成对角矩阵Λ∈R^r×r，将所有得到的向量μ₁，μ₂，…，μ_r组成矩阵U＝[μ₁，μ₂，…，μ_r]，再将所有的向量v₁，v₂，…，v_r组成矩阵V＝[v₁，v₂，…，v_r]。

(7)根据公式计算监测统计量的控制上限D_lim，并保留参数集Θ＝{V，Λ，D_lim}以备在线监测时调用，其中表示自由度为r的卡方分布在置信度α＝99％下的取值。

(8)收集新采样时刻的数据样本x_t∈R^1×m，引入其前d个采样时刻的样本得到增广向量x_a＝[x_t，x_t-1，…，x_t-d]，其中t表示当前采样时刻。

(9)对x_a实施与步骤(2)同样的标准化处理得到并初始化i＝1。

(10)将行向量中的第i个元素y_i取出后作为AR模型输出，而将向量中其余元素组成AR模型输入并根据步骤(3)中训练得到的AR模型计算y_i的预测误差f_i：

(11)判断是否满足条件i＜m？若是，则置i＝i+1后返回步骤(10)；若否，则将得到的误差组成向量f＝[f₁，f₂，…，f_m]并继续执行下一步骤(12)。

(12)对f实施与步骤(5)中相同的标准化处理得到新向量调用步骤(7)中保存的参数集Θ对进行转换得到向量μ，即：

(13)根据公式D＝μμ^T计算监测统计指标D，判断是否满足条件D＞D_lim？若否，则当前样本为正常工况采样；若是，则当前采样数据有可能来自故障工况。

与传统方法相比，本发明方法的优势在于：

首先，本发明方法利用PLS算法为每个测量变量建立其各自的AR模型，体现出了分散式建模的特点，因此具备分散式多模型的优势。其次，本发明方法只利用AR模型的预测误差作为被监测对象，而误差是AR模型中剔除了自相关性与交叉相关性后的结果，通过这种思路巧妙地避免了动态过程监测中自相关性问题。再者，本发明方法为误差实施奇异值分解后只建立一个统计量指标，与传统计算两个指标相比，本发明方法给出的监测结果更单一而不会因组合出现多种过程监测结果。

附图说明

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

图2为TE过程反应器冷却水进口温度故障的监测详情对比图。

具体实施方式

下面结合附图与具体的实施案例对本发明方法进行详细的说明。

如图1所示，本发明公开一种基于分散式AR-PLS模型的动态过程监测方法。下面结合一个具体的工业过程的例子来说明本发明方法的具体实施过程，以及相对于现有方法的优越性。

应用对象是来自于美国田纳西-伊斯曼(TE)化工过程实验，原型是伊斯曼化工生产车间的一个实际工艺流程。目前，TE过程因其流程的复杂性，已作为一个标准实验平台被广泛用于故障检测研究。整个TE过程包括22个测量变量、12个操作变量、和19个成分测量变量。所采集的数据分为22组，其中包括1组正常工况下的数据集与21组故障数据。而在这些故障数据中，有16个是已知故障类型，如冷却水入口温度或进料成分的变化、阀门粘滞、反应动力学漂移等，还有5个故障类型是未知的。为了对该过程进行监测，选取如表1所示的33个过程变量，接下来结合该TE过程对本发明具体实施步骤进行详细的阐述。

表1：TE过程监测变量。

序号	变量描述	序号	变量描述	序号	变量描述
						1	物料A流量	12	分离器液位	23	D进料阀门位置
2	物料D流量	13	分离器压力	24	E进料阀门位置
						3	物料E流量	14	分离器塔底流量	25	A进料阀门位置
4	总进料流量	15	汽提塔等级	26	A和C进料阀门位置
						5	循环流量	16	汽提塔压力	27	压缩机循环阀门位置
6	反应器进料	17	汽提塔底部流量	28	排空阀门位置
						7	反应器压力	18	汽提塔温度	29	分离器液相阀门位置
8	反应器等级	19	汽提塔上部蒸汽	30	汽提塔液相阀门位置
						9	反应器温度	20	压缩机功率	31	汽提塔蒸汽阀门位置
10	排空速率	21	反应器冷却水出口温度	32	反应器冷凝水流量
						11	分离器温度	22	分离器冷却水出口温度	33	冷凝器冷却水流量

首先，利用TE过程正常工况下的采样数据建立故障检测模型，包括以下步骤：

(1)收集生产过程正常运行状态下的数据样本，组成训练数据集X∈R^960×33，并构造增广型矩阵X_a∈R^958×99；

(2)对矩阵X_a中每列实施标准化处理，得到均值为0，标准差为1的新数据矩阵后初始化i＝1；

(3)将矩阵中第i列数据取出作为AR模型的输出，矩阵中剩余的列则作为AR模型的输入矩阵并训练得到与输出之间的AR模型；

(4)判断是否满足条件i＜33？若是，则置i＝i+1后返回步骤(3)；若否，则将得到的估计误差向量组成矩阵E＝[e₁，e₂，…，e_m]∈R^958×33后继续执行下一步骤(5)；

(5)将估计误差矩阵E作为新的训练数据矩阵，对其中的每一列实施标准化处理，得到均值为0，标准差为1的新数据矩阵

(6)对矩阵实施奇异值分解，即：

(7)根据公式计算监测统计量的控制上限D_lim，并保留参数集Θ＝{V，Λ，D_lim}以备在线监测时调用；

其次，采集TE过程在反应器冷却水进口温度故障条件下的测试数据集，实施在线过程监测。

(8)收集新采样时刻的数据样本x_t∈R^1×33，引入其前d＝2个采样时刻的样本得到增广向量x_a＝[x_t，x_t-1，x_t-2]；

(9)对x_a实施与步骤(2)同样的标准化处理得到并初始化i＝1；

(10)：将行向量中的第i个元素y_i取出后作为AR模型输出，而将向量中其余元素组成AR模型输入并根据步骤(3)中训练得到的AR模型计算y_i的预测误差f_i：

(11)判断是否满足条件i＜m？若是，则置i＝i+1后返回步骤(10)；若否，则将得到的误差组成向量f＝[f₁，f₂，…，f_m]并继续执行下一步骤(12)；

(12)对f实施与步骤(5)中相同的标准化处理得到新向量调用步骤(7)中保存的参数集Θ对进行转换得到向量

(13)根据公式D＝μμ^T计算监测统计指标D，判断是否满足条件D＞D_lim？若否，则当前样本为正常工况采样；若是，则当前采样数据有可能来自故障工况。

最后，将本发明方法与传统动态PCA方法的过程监测详情对比于如图2中。从图2中可以发现，本发明方法的监测效果要优越于传统动态PCA方法。

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

Claims (3)

1.一种基于分散式AR-PLS模型的动态过程监测方法，其特征在于，包括以下步骤：

离线建模阶段的实施过程如下所示：

步骤(1)：采集生产过程正常运行状态下的样本，组成训练数据集X∈R^N×m，并按照如下形式构造增广型矩阵X_a∈R^{(N-d)×m(d+1)}：

其中，N为训练样本数，m为过程测量变量数，R为实数集，R^N×m表示N×m维的实数矩阵，x_n∈R^1×m为第n个采样时刻的样本数据，下标号n＝1，2，…，N，d为引入的延时测量值的个数(一般取d＝1或2)；

步骤(2)：对矩阵X_a中每列实施标准化处理，得到均值为0，标准差为1的新数据矩阵后初始化i＝1，并将表示成如下形式：

其中，上标号1，2，…，m为测量变量标号，下标号0，1，…，d为当前时刻(即0)采样值及其d个延时测量值(即1，2，…，d)的标号，表示延时d个时刻测量到的m个N-d维的列向量；

步骤(3)：将矩阵中第i列数据取出作为自回归(Auto Regression，AR)模型的输出，矩阵中剩余的列则作为AR模型的输入矩阵再利用偏最小二乘(Partial Least Squares，PLS)算法训练得到与输出之间的AR模型：

其中，b_i为输入-输出回归向量，e_i表示第i个变量的预测误差；

步骤(4)：判断是否满足条件i＜m？若是，则置i＝i+1后返回步骤(3)；若否，则将得到的估计误差向量组成矩阵E＝[e₁，e₂，…，e_m]∈R^(N-d)×m后继续执行下一步骤(5)；

步骤(5)：将估计误差矩阵E作为新的训练数据矩阵，对其中的每一列实施标准化处理，得到均值为0，标准差为1的新数据矩阵

步骤(6)：对矩阵实施奇异值分解(Singular Value Decomposition，SVD)，即：

步骤(7)：根据公式计算监测统计量的控制上限D_lim，并保留参数集Θ＝{V，Λ，D_lim}以备在线监测时调用，其中表示自由度为r的卡方分布在置信度α＝99％下的取值；

在线过程监测阶段的实施过程如下所示：

步骤(8)：收集新采样时刻的数据样本x_i∈R^1×m，引入其前d个采样时刻的样本得到增广向量x_a＝[x_i，x_i-1，…，x_i-d]，其中t表示当前采样时刻；

步骤(9)：对x_a实施与步骤(2)同样的标准化处理得到并初始化i＝1；

步骤(10)：将行向量中的第i个元素y_i取出后作为AR模型输出，而将向量中其余元素组成AR模型输入并根据步骤(3)中训练得到的AR模型计算y_i的预测误差f_i：

步骤(11)：判断是否满足条件i＜m？若是，则置i＝i+1后返回步骤(10)；若否，则将得到的误差组成向量f＝[f₁，f₂，…，f_m]并继续执行下一步骤(12)；

步骤(12)：对f实施与步骤(5)中相同的标准化处理得到新向量调用步骤(7)中保存的参数集Θ对进行转换得到向量μ，即：

步骤(13)根据公式D＝μμ^T计算监测统计指标D，判断是否满足条件D＞D_lim？若否，则当前样本为正常工况采样；若是，则当前采样数据有可能来自故障工况。

2.根据权利要求1所述的一种基于分散式AR-PLS模型的动态过程监测方法，其特征在于，所述步骤(3)中利用PLS算法训练得到与输出之间的AR模型的具体实施过程如下所示：

①置k＝1与后，初始化向量

②根据公式w_k＝Z_iu_k/||Z_iu_k||、s_k＝Z_iw_k、分别计算得到系数向量w_k、得分向量s_k、和系数q_k，其中||Z_iu_k||表示计算向量Z_iu_k的长度；

③根据公式计算向量u_new；

④判断是否满足条件||u_k-u_new||＜10^-6？若否，则置u_k＝u_new后返回步骤②；若是，则执行⑤；

⑤根据公式p_k＝Z_i ^Ts_k/(s_k ^Ts_k)计算得到第k个投影向量p_k∈R^{[m(d+1)-1]×1}，并保留向量p_k、向量w_k、和系数q_k；

⑥判断矩阵Y_k＝s_kp_k ^T中的最大元素是否大于0.01？若是，根据公式Z_i＝Z_i-s_kp_k ^T更新矩阵Z_i后执行步骤⑦；若否，则得到投影矩阵P_i＝[p₁，p₂，…，p_k]、系数矩阵W_i＝[w₁，w₂，…，w_k]、和列向量Q＝[q₁，q₂，…，q_k]^T，并执行步骤⑧；

⑦判断k＜m(d+1)-1？若是，则置k＝k+1后，返回步骤②；若否，则得到最终的投影矩阵P_i＝[p₁，p₂，…，p_k]与系数矩阵W_i＝[w₁，w₂，…，w_k]；

⑧计算回归向量b_i＝W_i(P_i ^TW_i)^-1Q，那么输入与输出x_i之间的AR模型为：其中，e_i表示第i个变量的预测误差。

3.根据权利要求1所述的一种基于分散式AR-PLS模型的动态过程监测方法，其特征在于，所述步骤(6)中对实施奇异值分解的具体实施过程如下所示：

①设置r＝1与后，初始化列向量t_r为矩阵中的第一列；

②根据公式v_r＝F^Tt_r/(t_r ^Tt_r)计算得到向量v_r；

③根据公式t_new＝Fv_r/(v_r ^Tv_r)计算向量t_new；

④判断是否满足条件||t_new-t_r||＜10^-6？若是，则执行下一步骤⑤；若否，则置t_r＝t_new后返回步骤②；

⑤根据公式与μ_r＝t_rλ_r ^-1分别计算得到第r个奇异值λ_r与向量μ_r，并根据公式F＝F-t_rv_r ^T更新F；

⑥判断是否满足条件λ_r≤10^-3？若否，则设置r＝r+1与向量t_r为矩阵F中的第一列后返回步骤②；若是，则将所有得到的奇异值λ₁，λ₂，…，λ_r组成对角矩阵Λ∈R^r×r，将所有得到的向量μ₁，μ₂，…，μ_r组成矩阵U＝[μ₁，μ₂，…，μ_r]，再将所有的向量v₁，v₂，…，v_r组成矩阵V＝[v₁，v₂，…，v_r]。