CN106990768A - 基于Limited‑DTW的MKPCA间歇过程故障监测方法 - Google Patents

Abstract

本发明公开一种基于Limited‑DTW的MKPCA间歇过程故障监测方法，针对间歇过程固有的批次不等长特性，也为了克服传统解决批次间同步问题方法数据浪费、扭曲原始过程变量的自相关及变量之间的交叉相关关系的严重缺陷，本发明方法引入了全局路径限制和失真度阈值限制对动态时间规整(Dynamic Time Warping,DTW)方法进行改进，避免算法长时间运行造成的监控弊端，并且解决了其处理过程的复杂性与其离线性导致其实际应用的困难。Limited‑DTW和基于MKPCA的监测方法的组合用于离线和在线实施；实验设计由青霉素发酵仿真平台和重组大肠杆菌实际生产过程完成，结果显示了本文方法的可行性和有效性。

Description

基于Limited-DTW的MKPCA间歇过程故障监测方法

技术领域

本发明涉及基于数据驱动的多元统计过程监控(Multivariate StatisticalProcess Monitoring,MSPM)技术领域，尤其涉及一种基于Limited-DTW的MKPCA间歇过程故障监测方法。

背景技术

基于数据驱动的多元统计方法在用于间歇过程的过程监控时，不需要考虑复杂的过程机理特性，通过对历史数据的建模和分析，判断生产过程的运行状态是否出现异常。间歇过程作为现代流程工业中的重要生产方式之一，被广泛用于生物医药、食品及生化产品的制备，特别是青霉素(Penicillin，或音译盘尼西林)药品的制备。与连续生产过程相比，间歇过程有着明显的区别，其中批次数据不等长是其固有特征之一。

为了运用传统的统计分析方法，假设批次数据是等长的，因此操作事件在批次上是同步化的。然而，这种假设通常是不成立的。另外，在进行间歇过程建模时，一般都基于一个默认假设，同一个时刻的数据属于同一阶段，即同一个阶段在所有的生产批次中都是等长的。但是实际生产过程中并不是如此，不同批次的生产由于种种条件的不同，导致不同批次间的阶段长度会有所差异。这也是导致批次不等长的原因之一。现有的解决批次不等长问题的方法，最常用的是最短长度法，即直接按最短批次数据轨迹的长度切割其余批次，方法简单，但使数据轨迹的过程大量丢失，并且使点对点数据的相关性降低，导致数据可靠性降低。Neogi和Undey均提出采用指示变量法解决批次间的同步问题，但是指示变量法可能会扭曲原始过程变量的自相关及变量之间的交叉相关关系。另外，并不是每个工业过程中都存在这样的“指示”变量。所以，在进行在线监控前，解决批次数据不等长问题是十分重要的。

发明内容

本发明针对间歇过程批次数据不等长问题，提供一种基于Limited-DTW的MKPCA间歇过程故障监测方法，通过使原始批次按照轨迹中点与点的模式进行动态的方法匹配解决批次不等长问题。

为实现上述目的，本发明采用如下的技术方案：

一种基于Limited-DTW的MKPCA间歇过程故障监测方法，包括以下步骤：

步骤1、针对青霉素的发酵过程数据进行数据预处理

将选取的青霉素发酵过程的10个批次在数据预处理时采用了沿批次展开方法，之后将沿批次对数据进行按列标准化处理，其中，10个批次数据为矩阵X^，；

步骤2、将提取的过程数据的平均轨迹进行缩放

选取B_RAW,i,i＝1,...,I为参考轨迹，它包含来自I个正常批次的原始测量值，I等于10，通过平均每个批次的范围来找到每个变量平均范围，然后存储这些值，将所有批次中的每个变量除以其平均范围，得到B_i，定义B_i,i＝1,...,I作为生成的缩放批次轨迹；

步骤3、利用Limited-DTW方法对不同批次轨迹进行同步处理

具体过程如下：在所有轨迹中选择一个轨迹B_k作为规整轨迹，令B_REF＝B_k，令W为Limited-DTW算法中的权重矩阵，设W等于单位矩阵，对指定的最大迭代次数执行以下步骤，

(1)在B_i,i＝1,...,I和B_REF之间进行Limited-DTW同步，同步时，在DTW算法的基础上，搜索空间减少为上一步所提出的图中菱形和两条平行线之间交叉的部分，

令为与B_REF同步后的批次轨迹，此时的持续时间与B_REF相等。

(2)计算平均轨迹

(3)根据计算每个变量方差的和，的倒数将是下一次迭代的特定变量的权重，即W将是一个对角矩阵：

其中，k代表迭代空间的横坐标，j代表迭代空间的纵坐标，

然后归一化W，使得权重的和等于变量的数量，即令W等于

(4)对于前三次迭代，保持相同的参考轨迹：B_REF＝B_k，

对于后续迭代，将参考设置为等于平均轨迹：

(5)规整后的各批次数据形成新的矩阵X，此时，各批次数据等长；

步骤4、发酵过程故障监测

对同步后的批次进行MKPCA建模，对于完成上述数据处理的数据X(I×kJ)进行MKPCA建模，其中，MKPCA模型表达式如下：

X＝TP^T+E

其中，P(J×R)为负载矩阵，T(I×kR)为得分矩阵，E(I×kJ)为残差矩阵，R为由累计方差贡献率保留的主元数目；

最终由MKPCA模型确定得分矩阵T，进而求得X的非线性主元。

作为优选，对DTW算法增加全局路径限制和失真度阈值限制，得到Limited-DTW方法。

本发明的基于Limited-DTW的MKPCA间歇过程故障监测方法，针对间歇过程固有的批次不等长特性，也为了克服传统解决批次间同步问题方法数据浪费、扭曲原始过程变量的自相关及变量之间的交叉相关关系的严重缺陷，本发明方法引入了全局路径限制和失真度阈值限制对动态时间规整(Dynamic Time Warping,DTW)方法进行改进，避免算法长时间运行造成的监控弊端，并且解决了其处理过程的复杂性与其离线性导致其实际应用的困难。Limited-DTW和基于MKPCA的监测方法的组合用于离线和在线实施；实验设计由青霉素发酵仿真平台和重组大肠杆菌实际生产过程完成，结果显示了本文方法的可行性和有效性。

有益效果

针对DTW在处理批次不等长问题中存在的缺陷，对其增加全局路径限制和失真度阈值限制，通过在全局最优路径的查找过程中约束路径的斜率的方法，可以大大减少全局路径的查找范围，进而减少约三分之二的计算量。而且节省存储空间，降低运算复杂度，进而提升算法的实时性，提高了监控性能。

附图说明

图1为基于Limited-DTW的MKPCA间歇过程故障监测方法流程图；

图2为DTW全局路径限制的示意图；

图3为Limited-DTW算法的示意图，其中，包括DTW失真度阈值限制；

图4a为最短长度法MKPCA对故障批次的监控结果T²统计量；

图4b为最短长度法MKPCA对故障批次的监控结果SPE统计量；

图5a为DTW MKPCA对故障批次的监控结果T²统计量；

图5b为DTW MKPCA对故障批次的监控结果SPE统计量；

图6a为Limited-DTW MKPCA对正常批次的监控结果T²统计量；

图6b为Limited-DTW MKPCA对正常批次的监控结果SPE统计量。

具体实施方式

动态时间规整(Dynamic Time Warping，DTW)算法最早用于语音识别领域，是结合时间规整和间距测量计算的非线性规整技术，是一种计算时间矢量序列间相似度的常用方法。动态时间规整算法是一种柔性模式匹配算法，能够对存在全局或局部扩展、压缩或变形的模式进行匹配，解决动态模式的相似度量和分类问题，其实质是运用动态规划思想，按局部最优自动寻找一条路径。这种方法虽然按照轨迹中点与点的模式进行动态匹配，但其处理过程的复杂性与其离线性导致其实际应用的困难。所以本发明提出了基于Limited-DTW的MKPCA间歇过程故障监测方法，该方法通过对DTW算法增加全局路径限制和设定失真度阈值限制，减少算法运算量，使其可以快速有效的用于批次不等长现实问题的解决。

青霉素(Penicillin，或音译盘尼西林)是一种常见的临床抗菌药品，其生产制备过程具有典型的批次间不等长特性。本文基于美国Illinois州立理工学院Cinar教授研究开发的具有广泛国际影响力的Pensim仿真平台对间歇过程进行在线监测仿真研究。青霉素发酵过程的生产批次持续时间为400h左右，采样间隔为1h。在线监控时，选取10个主要的过程变量用于监控过程的运行状况，选取的过程变量如表1所示。本文选取了青霉素发酵过程的10个批次进行实验仿真。

表1过程变量

Tab.1Process Variables

基于以上描述，按照发明内容，将具体过程在Matlab中实现如下，如图1所示：

Ⅰ进行数据预处理。MKPCA作为MSPM技术的核心方法，在将其用于发酵过程监测时，面对三维形式的数据需要进行必要的预处理操作。本文将选取的青霉素发酵过程的10个批次在数据预处理时采用了沿批次展开方法，之后将沿批次对数据进行按列标准化处理，提取过程数据的平均轨迹。其中，10个批次数据为矩阵X^，。

Ⅱ将提取的过程数据的平均轨迹进行缩放。选取B_RAW,i,i＝1,...,I为参考轨迹，它包含来自I个正常批次的原始测量值，本文中I等于10。通过平均每个批次的范围来找到每个变量平均范围，然后存储这些值，将所有批次中的每个变量除以其平均范围，得到B_i，定义B_i,i＝1,...,I作为生成的缩放批次轨迹。

Ⅲ对DTW算法增加全局路径限制和失真度阈值限制，即本文提出的Limited-DTW方法。Limited-DTW方法对DTW算法改进如下:

(1)对DTW算法增加全局路径限制。根据动态规划条件，增加全局路径限制。DTW算法虽然简洁，对资源的硬件要求也较小，但运算量很大，能否减少运算量又不降低识别率，对DTW算法的实际应用非常关键。随着索引序列库的不断增大，会影响动态识别的效率，这将是一个严重缺点。针对上述问题，本文提出了对DTW算法的全局路径限制。如图2所示，实线为全局路径限制线，在算法进行匹配路径搜索时，则只需要对中间菱形的部分进行搜索。根据发酵过程的数据特性，通过交叉验证法，选取约束斜率k满足此时在减少计算量的同时能最大限度的保证数据的完整性。

增加全局路径之后，算法在寻找匹配的最优路径时，只计算菱形之内某点和其之前点的匹配距离，比较过程如下：

当X_A＝X_B时，比较分为两段：

当X_A＜X_B时：

当X_a＞X_b时，比较方法与上式类似。

其中，M和N分别为参与匹配的两组特征向量序列的长度。A为测试特征向量序列。

(2)对DTW算法增加失真度阈值限制。随着特征向量序列长度变大，该算法的复杂度也会随着增大。假设测试特征向量序列为A，则该特征向量序列与A的模之间的失真度应最小。最优路径长度均处于M、N中较大值与M+N之间，最优路径中不匹配的次数与最优路径长度成正比关系，选用α×(M+N)(α为正比例系数)作为失真度阈值。本文综合考虑上述因素以及通过实验结果的验证及分析，选取失真度阈值为0.25×(M+N)，如图3所示。

全局路径限制和失真度阈值限制结合使用，如图3所示，在这种情况下，搜索空间将是图中菱形和两条平行线之间交叉的部分。

Ⅳ利用Limited-DTW方法对不同批次轨迹进行同步处理，具体算法如下。在所有轨迹中选择一个轨迹B_k作为规整轨迹，令B_REF＝B_k。

令W为Limited-DTW算法中的权重矩阵，设W等于单位矩阵，对指定的最大迭代次数执行以下步骤。

(1)在B_i,i＝1,...,I和B_REF之间进行Limited-DTW同步。同步时，在DTW算法的基础上，搜索空间减少为上一步所提出的图中菱形和两条平行线之间交叉的部分。

令为与B_REF同步后的批次轨迹，此时的持续时间与B_REF相等。

(2)计算平均轨迹

(3)根据计算每个变量方差的和。的倒数将是下一次迭代的特定变量的权重，即W将是一个对角矩阵：

其中，k代表迭代空间的横坐标，j代表迭代空间的纵坐标。

然后归一化W，使得权重的和等于变量的数量，即令W等于

(4)对于前三次迭代，保持相同的参考轨迹：B_REF＝B_k。

对于后续迭代，将参考设置为等于平均轨迹：

迭代过程结束时的同步轨迹的长度将等于最初用作参考批次的轨迹的长度。或者，可以从初始轨迹估计平均持续时间，并且其持续时间最接近平均持续时间的轨迹可以用作前三次迭代的B_REF。通过这样做，结束时的同步轨迹的持续时间将为平均持续时间。初始参考轨迹的选择是用户偏好的问题。最大迭代次数是用户设置的方法的另一个参数。还可以监测权重矩阵W从一次迭代到下一次迭代的变化，并将其用作趋同的指标。

(5)规整后的各批次数据形成新的矩阵X，此时，各批次数据等长。

Ⅴ发酵过程故障监测。

对同步后的批次进行MKPCA建模；对于完成上述数据处理的数据X(I×kJ)进行MKPCA建模，

MKPCA模型表达式如下：

X＝TP^T+E (7)

其中P(J×R)为负载矩阵，T(I×kR)为得分矩阵，E(I×kJ)为残差矩阵，R为由累计方差贡献率保留的主元数目。

最终由MKPCA模型确定得分矩阵T，进而求得X的非线性主元。

发酵仿真使用在相近的初始条件下产生的30个批次的正常过程数据，分别按照T²统计量和SPE统计量所遵循的概率分布求取控制限。

本文选用的待监测的故障批次是由底物流加速率在200时刻引入1％的斜坡故障产生，延续至反应结束。为了保证仿真环境与实际生产现场环境的一致性，对训练样本集加入了一定的白噪声干扰。

为了验证本文方法用于过程监控的有效性，文中分别本文引入了传统DTW方法和Limited-DTW方法对30批数据进行批次间不等长处理建立正常批次模型后对故障批次进行监测对比分析验证。为进一步表明本文方法的优越性，本文还引入了最短长度法方法处理批次数据与本文方法进行对比分析。从图4a、4b、5a、5b、6a、6b中可以看出，最短长度法处理数据后建立的MKPCA模型的T²统计量在208时刻检测到故障，SPE统计量在209时刻左右检测到故障，但其存在较高的误报率和漏报率；传统DTW在故障检测的准确率上有所提高，但对故障的检测存在较大的延迟，其T²统计量在222时刻检测到故障，SPE统计量在223时刻检测到故障。而本文提出的Limited-DTW方法对DTW方法在寻优过程中增加了全局路径限制和失真度阈值限制，减少了对数据处理的计算量，因此较传统DTW方法在故障检测的实时性上有了很大提高，同时减小了误报率和漏报率。

Claims (2)

1.一种基于Limited-DTW的MKPCA间歇过程故障监测方法，其特征在于，包括以下步骤：

步骤1、针对青霉素的发酵过程数据进行数据预处理

将选取的青霉素发酵过程的10个批次在数据预处理时采用了沿批次展开方法，之后将沿批次对数据进行按列标准化处理，其中，10个批次数据为矩阵X，；

步骤2、将提取的过程数据的平均轨迹进行缩放

选取B_RAW,i,i＝1,...,I为参考轨迹，它包含来自I个正常批次的原始测量值，I等于10，通过平均每个批次的范围来找到每个变量平均范围，然后存储这些值，将所有批次中的每个变量除以其平均范围，得到B_i，定义B_i,i＝1,...,I作为生成的缩放批次轨迹；

步骤3、利用Limited-DTW方法对不同批次轨迹进行同步处理

具体过程如下：在所有轨迹中选择一个轨迹B_k作为规整轨迹，令B_REF＝B_k，令W为Limited-DTW算法中的权重矩阵，设W等于单位矩阵，对指定的最大迭代次数执行以下步骤，

(1)在B_i,i＝1,...,I和B_REF之间进行Limited-DTW同步，同步时，在DTW算法的基础上，搜索空间减少为上一步所提出的图中菱形和两条平行线之间交叉的部分，

令为与B_REF同步后的批次轨迹，此时的持续时间与B_REF相等。

(2)计算平均轨迹

(3)根据计算每个变量方差的和，的倒数将是下一次迭代的特定变量的权重，即W将是一个对角矩阵：

$W(j,j)={\[\underset{i=1}{\overset{I}{\Σ}}\underset{k=1}{\overset{b_{REF}}{\Σ}}{\[{\tilde{B}}_i(k,j)-\stackrel{\‾}{B}(k,j)\]}^2\]}^{-1}$

其中，k代表迭代空间的横坐标，j代表迭代空间的纵坐标，

然后归一化W，使得权重的和等于变量的数量，即令W等于

(4)对于前三次迭代，保持相同的参考轨迹：B_REF＝B_k，

对于后续迭代，将参考设置为等于平均轨迹：

(5)规整后的各批次数据形成新的矩阵X，此时，各批次数据等长；

步骤4、发酵过程故障监测

对同步后的批次进行MKPCA建模，对于完成上述数据处理的数据X(I×kJ)进行MKPCA建模，其中，MKPCA模型表达式如下：

X＝TP^T+E

其中，P(J×R)为负载矩阵，T(I×kR)为得分矩阵，E(I×kJ)为残差矩阵，R为由累计方差贡献率保留的主元数目；

最终由MKPCA模型确定得分矩阵T，进而求得X的非线性主元。

2.如权利要求1所述的基于Limited-DTW的MKPCA间歇过程故障监测方法，其特征在于，对DTW算法增加全局路径限制和失真度阈值限制，得到Limited-DTW方法。