一种基于多向主元分析法的多阶段批次过程阶段划分方法

摘要

本发明涉及一种基于多向主元分析法的多阶段批次过程阶段划分方法，包括模型数据采集、三维数据展开、二维矩阵标准化、PCA分解、主元个数选取、负载矩阵分解和K-means聚类分析来实现，本发明首次将多向主元分析法模型的负载矩阵按照时间片顺序经过分解变形来辨识过程阶段，打破了传统理论中多向主元分析法不适用于多阶段批次过程建模的认识误区，该方法不仅提取了各个时间片上变量的交叉相关性，同时还提取了测量变量在整个批次各个采样时刻的自身的动态变化特性来辨识阶段，提高了过程监测过程的监测效率和故障诊断结果的准确性，为无过程先验知识条件下的多阶段过程阶段划分提供了新的可能性。

说明书

技术领域

本发明涉及一种基于多向主元分析法的多阶段批次过程阶段划分方法，属于间歇过程多 变量监测与故障监测领域。

背景技术

间歇过程操作中各操作变量的相关关系并非随时间时刻变化，而是随着过程操作进程或 者过程机理特性的变化发生而变化，多阶段是间歇过程的一个显著特点，了解间隙过程每个 阶段的变量关系，有助于提高间隙过程的监测效率，增强诊断的可靠性。

目前针对带有阶段特性的批次过程的多变量过程阶段的划分方法主要有三种，第一种是 依靠过程专家的经验对阶段进行辨识，不过显然，这种方法过于依赖于人，判断结果过于主 观；第二种方法是通过某一关键过程变量的轨迹异常点来判断，但这种方法没有考虑过程的 变量的相关性，阶段辨识比较粗糙；第三种方法是Sub-PCA法，它通过提取过程的时间片矩 阵建模，按照建模后的每个时间片的负载矩阵来聚类分阶段，但是这种方法是提取过程每个 时间片矩阵的相关性信息，忽略了每个时间片矩阵整个过程的动态变化信息，阶段划分结果 比较片面。

发明内容

本发明的目的，就是克服现有技术的不足，提供一种基于多向主元分析法的多阶段批次 过程阶段划分方法，该方法提出一种无过程先验知识的过程阶段辨识方法，为多阶段批次过 程的建模以及监测提供了新的途径。

为了达到上述目的，采用如下技术方案：

一种基于多向主元分析法的多阶段批次过程阶段划分方法，包括以下步骤得到，

1)模型数据采集

设一个间歇操作具有J个测量变量和K个采样点，则每一个测量批次可得到一个J×K的 矩阵，重复I批次的测量步骤后，得到的数据可以表述为一个三维矩阵X(I×J×K)，其中测 量变量为温度、速度、压力、行程等批次运行过程中可被测量的状态参数；

2)三维数据展开

将三维矩阵X按照采集批次方向展开，即将一个操作批次内的各采样点上的变量按照时 间顺序排开得到二维矩阵显然矩阵为I行KJ列；

3)二维矩阵标准化

设二维矩阵内任意一点的变量为对该变量进行减均值、除方差的标准化处理，标 准化处理的计算公式如下：

$x_{\mathrm{ijk}}=\frac{{\underset{~}{x}}_{\mathrm{ijk}}-{\bar{x}}_{\mathrm{jk}}}{s_{\mathrm{jk}}};---\left(1\right)$

其中：是矩阵任一列的均值，s_jk是矩阵任一列的方差，

${\bar{x}}_{\mathrm{jk}}=\frac{1}{I}\underset{i=1}{\overset{I}{\Sigma }}{\underset{~}{x}}_{\mathrm{ijk}},$

$s_{\mathrm{jk}}=\sqrt{\underset{i=1}{\overset{I}{\Sigma }}{({\underset{~}{x}}_{\mathrm{ijk}}-{\bar{x}}_{\mathrm{jk}})}^2/(I-1)};$

4)多向主元分析法建模

对上一步经标准化后的二维矩阵(I×JK)执行PCA分解，完成多向主元分析法的建模， 其中PCA分解公式如下：

$X={\mathrm{TP}}^T=\underset{i=1}{\overset{\mathrm{JK}}{\Sigma }}{t}_i{p}_i;---\left(2\right)$

S＝trace(T^TT/(I-1))；                    (3)

其中：t_i为正交的主元向量，p_i为正交归一化的负载向量，S是主元的协方差矩阵的迹， 代表各个主元对于过程的解释度大小；

公式(2)X分解得到得分矩阵T(I×JK)及负载矩阵P(JK×JK)；

5)选取主元个数

将公式(2)重新表述成如下形式：

$X=\underset{r=1}{\overset{R}{\Sigma }}{t}_r{p}_r+\underset{i=R+1}{\overset{\mathrm{JK}}{\Sigma }}{t}_i{p}_i=T_r{P}_r+E---\left(4\right)$

其中：T_r(I×R)、P_r(JK×R)分别为保留R个主元后的得分矩阵和负载矩阵，E为残差矩 阵；

通过上述变换，多向主元分析法模型将原始数据空间分解为主元空间和残差空间，主元 空间变量高度相关，一般来说足以描述数据的变异性；

主元个数R一般可根据用户的经验设定或者采用Broken-Stick准则，Broken-Stick的内容 是当第r个主元的解释度S(r)占所有主元总贡献sum(S)的百分比大于G(r)的时候保留该主 元，否则终止，其中G(r)的计算公式如下：

$G\left(r\right)=100\frac{1}{z}\underset{i=r}{\overset{z}{\Sigma }}\frac{1}{i};---\left(5\right)$

其中：S(r)是第r个主元的解释度，sum(S)是所有主元的贡献和；

6)负载矩阵分解

按照步骤(1)的方式将包含了JK个变量的相关性的三维负载矩阵P_r(JK×R)按照时间 片批次展开，由于负载矩阵P_r(JK×R)每J行是一个时间片负载矩阵则将负 载矩阵P_r按照时间次序叠加起来即可得到三维矩阵

将沿时间K方向分解可构成二维矩阵其代表每个采样点各个变量之间的 相关关系，为方便说明，将称作相关模式，则可以表述成K个相关模式组合的结构形式：

$\hat{P}=\{{\hat{P}}_1^K,{\hat{P}}_2^K,···{\hat{P}}_k^K,···{\hat{P}}_K^K\};---\left(6\right)$

7)K-means聚类分析

本步骤选用下面定义的距离作为度量两个相关模式相似性程度的指标，对K个相关模式 进行K-means聚类分析，使得具有相同阶段特性的相关模式划分为一类，不同的 类别代表不同的阶段特性，上述距离由下式定义：

$\mathrm{dist}({\hat{p}}_1^K,{\hat{p}}_2^K)={\left(\underset{j=1}{\overset{J}{\Sigma }}{({\hat{p}}_{1,j}^K,-{\hat{p}}_{2,j}^K)}^T({\hat{p}}_{1,j}^K-{\hat{p}}_{2,j}^K)\right)}^{1/2};---\left(7\right)$

K-means算法的输入是K个相关模式集合以及两个子类中心的最 小距离阈值θ，算法的输出是子类数量C，子类中心设为{W₁，W₂，…，W_C}，以及每个相关模式 属于不同子类的隶属关系变量i是算法中迭代次数的索引，k是分 类模式的索引，而c则是聚类中心的索引，算法步骤如下：

a、从K个相关模式中，任意选择C₀个相关模式作为初始聚类中心W_i，c(c＝1，2，…，C₀)， 对于W_i，c的选取，常用方法是从被分类模式中均匀抽取C₀个相关模式，建议C₀在区间 (K/3～K/2)内取值；

b、若两个子类中心的距离dist(W_i，c1，W_i，c2)小于预定的阈值θ，则剔除其中一个聚类中心；

c、计算每个相关模式到所有聚类中心的距离若和第c^*类的中心的距离最小，则将的隶属关系定义为m(k)＝c^*；

d、I_num次迭代后，若某子类中心没有俘获一定数量的相关模式(例如没有超过5个相关 模式)，则剔除该奇类；

e、更新子类数量为C_i+1，并根据相关模式的隶属关系重新计算新的聚类中心 W_i+1，c(c＝1，2，…，C_i+1)；

如果算法满足收敛条件则结束，否则返回步骤b，进行下一次迭代计算，以上过程使具有 相同阶段特性的相关模式划分为一类，完成对多变阶段批次过程的阶段划分。

作为一种具体实施例，所述K-means聚类分析中，所述K-means聚类分析中，K-means 算法的收敛条件为两次迭代中的聚类中心距离的变化小于ε，其中

作为另一种具体实施例，所述K-means聚类分析中，K-means算法的收敛条件为每个子 类中相关模式到子类中心的距离平方和达到最小或子类之间的距离平方和达到最小。

与现有技术相比，本发明的有益效果在于：

本发明首次将多向主元分析法模型的负载矩阵按照时间片顺序经过分解变形来辨识过程 阶段，打破了传统理论中多向主元分析法不适用于多阶段批次过程建模的认识误区，该方法 不仅提取了各个时间片上变量的交叉相关性，同时还提取了测量变量在整个批次各个采样时 刻的自身的动态变化特性来辨识阶段，提高了过程监测过程的监测效率和故障诊断结果的准 确性，为无过程先验知识条件下的多阶段过程阶段划分提供了新的可能性。

附图说明

图1是本发明所述基于多向主元分析法的多阶段批次过程阶段划分方法的流程图。

图2是本发明所述基于多向主元分析法的多阶段批次过程阶段划分方法的三维数据展开 模型示意图。

图3是本发明具体实施例中注塑过程的阶段划分结果图。

具体实施方式

下面结合附图及具体实施例，对本发明做进一步说明：

注塑成型是典型的多阶段批次过程，其主要包含注射、保压、塑化、冷却四个阶段，注 塑过程的每个时段都有其特定的控制目标，不同的主导变量以及过程特性。具体来说，在注 射段，液压缸的高压推动螺杆向前将机桶内的熔融塑料推到模腔，当模腔被完全或者将近填 充满的时候，过程切换到保压阶段，在该阶段，高压将少量材料继续填充到模腔中以补充由 于冷却和固化带来的材料收缩；当胶口冷却，模腔中的材料不再被注射喷嘴影响的时候，保 压段结束。螺杆旋转并后退，将足够量的熔融塑料推到螺杆前端。螺杆后退同时开始容积计 算。头部熔料达到一定的注射量后，螺杆停止后退和转动，这段时间的过程状态称为塑化阶 段。在保压段结束，塑化过程进行的时候，冷却阶段也同时进行着，直到模具内材料达到能 被弹出的硬度，冷却阶段结束。

以上述注塑成型过程为例，参见图1，本发明所述的基于多向主元分析法(Multi-way  principal component analysis，MPCA)的多阶段批次过程阶段划分方法，包括以下步骤得到：

(1)模型数据采集

设一个间歇操作具有J个测量变量和K个采样点，则每一个测量批次可得到一个J×K的 矩阵，重复I批次的测量步骤后，得到的数据可以表述为一个三维矩阵X(I×J×K)。为了确 保检测数据涵盖足够长时间的工作范围，一般工业上用来建模的数据批次I的取值大于100， 测量变量为温度、速度、压力、行程等批次运行过程中可被测量的状态参数；基于过程时间 长短、过程变化的快慢程度以及模型负担是否在合理的范围，采样点K个数一般小于1000。

本实施例中，测量变量实验室注塑机工作过程可获得的变量为8个：压力阀门开度，流 量阀门开度，注射行程，注射速度，注射压力，机桶温度(3段)，操作批次I取100，每个批 次保留的采样点K为488。

(2)三维数据展开

参见图2，将三维矩阵X按照采集批次方向展开，即将一个操作批次内的各采样点上的 变量按照时间顺序排开得到二维矩阵显然矩阵为I行KJ列。

(3)二维矩阵标准化

设二维矩阵内任意一点的变量为对该变量进行先减均值后除方差的标准化处理， 标准化处理的计算公式如下：

$x_{\mathrm{ijk}}=\frac{{\underset{~}{x}}_{\mathrm{ijk}}-{\bar{x}}_{\mathrm{jk}}}{s_{\mathrm{jk}}};---\left(1\right)$

其中：是矩阵任一列的均值，s_jk是矩阵任一列的方差；

${\bar{x}}_{\mathrm{jk}}=\frac{1}{I}\underset{i=1}{\overset{I}{\Sigma }}{\underset{~}{x}}_{\mathrm{ijk}};$

$s_{\mathrm{jk}}=\sqrt{\underset{i=1}{\overset{I}{\Sigma }}{({\underset{~}{x}}_{\mathrm{ijk}}-{\bar{x}}_{\mathrm{jk}})}^2/(I-1)}.$

本步骤的标准化处理相当于抽取了间歇过程中一次操作的平均运行轨迹，突出了间歇过 程不同操作批次之间的一种正常随机波动。

(4)MPCA建模

所谓的MPCA建模就是先将三维矩阵展开成一个大的二维矩阵，再执行常规的PCA分解 的方法，本步骤对上一步经标准化处理后的二维矩阵(I×JK)执行PCA分解，其分解公式如 下：

$X={\mathrm{TP}}^T=\underset{i=1}{\overset{\mathrm{JK}}{\Sigma }}{t}_i{p}_i;---\left(2\right)$

S＝trace(T^TT/(I-1))；                        (3)

其中：t_i为正交的主元向量，p_i为正交归一化的负载向量，S是主元的协方差矩阵的迹， 代表各个主元对于过程的解释度大小。

公式(2)X分解得到得分矩阵T(I×JK)及负载矩阵P(JK×JK)。

(5)选取主元个数

一般来说，前几个主元一般包含着间隙过程的大部分变异信息，其他的主元可能主要包 含噪声信息，因此公式(2)可以被重新表述成如下形式：

$X=\underset{r=1}{\overset{R}{\Sigma }}{t}_r{p}_r+\underset{i=R+1}{\overset{\mathrm{JK}}{\Sigma }}{t}_i{p}_i=T_r{P}_r+E---\left(4\right)$

其中：T_r(I×R)、P_r(JK×R)分别为保留R个主元后的得分矩阵和负载矩阵，E为残差矩 阵；

通过上述变换，MPCA模型将原始数据空间分解为主元空间和残差空间，主元空间变量 高度相关，一般来说足以描述数据的变异性。

主元个数R一般可根据用户的经验设定或者采用Broken-Stick准则，Broken-Stick的内容 是当第r个主元的解释度S(r)占所有主元总贡献sum(S)的百分比大于G(r)的时候保留该主 元，否则终止，其中G(r)的计算公式如下：

$G\left(r\right)=100\frac{1}{z}\underset{i=r}{\overset{z}{\Sigma }}\frac{1}{i};---\left(5\right)$

其中：S(r)是第r个主元的解释度，sum(S)是所有主元的贡献和，在本实施例中，主元R 的个数选择5，对于过程的解释度为56.64。

(6)负载矩阵分解

按照步骤(1)的方式将包含了JK个变量的相关性的三维负载矩阵P_r(JK×R)按照时间 片批次展开，由于负载矩阵P_r(JK×R)每J行是一个时间片负载矩阵则将负 载矩阵P_r按照时间次序叠加起来即可得到三维矩阵

将沿时间K方向分解可构成二维矩阵其代表每个采样点各个变量之间的 相关关系，为方便说明，本说明书将称作相关模式，则可以表述成K个相关模式组合的 结构形式：

$\hat{P}=\{{\hat{P}}_1^K,{\hat{P}}_2^K,···{\hat{P}}_k^K,···{\hat{P}}_K^K\};---\left(6\right)$

(7)K-means聚类分析

对于一个具有多阶段特性的批次过程，不同的阶段过程关系结构会表现出不同，但是在 一个阶段内部，各个采样点上的过程特性基本相同，K-means算法是硬聚类算法，是典型的 局域原型的目标函数聚类方法的代表，它是数据点到原型的某种距离作为优化的目标函数， 利用函数求极值的方法得到迭代运算的调整规则，其是将n个数据对象划分为m个聚类以使 所获得的聚类满足同一聚类中的对象相似度较高而不同聚类中的对象相似度较小的一种算 法，聚类相似度是利用各聚类中对象与“中心对象”的距离来进行计算，

本方法选用下面定义的距离作为度量两个相关模式相似性程度的指标，对K个相关模式 进行K-means聚类分析，使得具有相同阶段特性的相关模式划分为一类，不同的 类别代表不同的阶段特性，上述距离由下式定义：

$\mathrm{dist}({\hat{p}}_1^K,{\hat{p}}_2^K)={\left(\underset{j=1}{\overset{J}{\Sigma }}{({\hat{p}}_{1,j}^K,-{\hat{p}}_{2,j}^K)}^T({\hat{p}}_{1,j}^K-{\hat{p}}_{2,j}^K)\right)}^{1/2};---\left(7\right)$

K-means算法的输入是K个相关模式集合以及两个子类中心的最 小距离阈值θ，算法的输出是子类数量C，子类中心设为{W₁，W₂，…，W_C}，以及每个相关模式 属于不同子类的隶属关系变量i是算法中迭代次数的索引，k是分 类模式的索引，而c则是聚类中心的索引，算法步骤如下：

a、从K个相关模式中，任意选择C₀个相关模式作为初始聚类中心W_i，c(c＝1，2，…，C₀)， 对于W_i，c的选取，常用方法是从被分类模式中均匀抽取C₀个相关模式，建议C₀在区间 (K/3～K/2)内取值；

b、若两个子类中心的距离dist(W_i，c1，W_i，c2)小于预定的阈值θ，则剔除其中一个聚类中心；

c、计算每个相关模式到所有聚类中心的距离若和第c^*类的中心的距离最小，则将的隶属关系定义为m(k)＝c^*；

d、I_num次迭代后，若某子类中心没有俘获一定数量的相关模式(例如没有超过5个相关 模式)，则剔除该奇类；

e、更新子类数量为C_i+1，并根据相关模式的隶属关系重新计算新的聚类中心 W_i+1，c(c＝1，2，…，C_i+1)；

如果算法满足收敛条件则结束，否则返回步骤b，进行下一次迭代计算，上述收敛条件有： 两次迭代中的聚类中心距离的变化小于一个很小的数ε，其中本实施例中 ε＝0.045；或者每个子类中模式到子类中心的距离平方达到最小或以及子类之间的距离平 方和达到最小。

以上过程使注塑过程中具有相同阶段特性的相关模式划分为一类，完成对多变阶段批次 过程的阶段划分，图3所示为对本实施例中注塑过程进行阶段划分的结果，这与实际的注塑 过程特性相吻合，该方式的应用使用户清楚地了解各变量在每个阶段的相互关系，提高了过 程监测效率并增加了诊断的有效性。

应该理解，本发明并不局限于上述具体实施例的注塑过程，凡是熟悉本领域的技术人员 在不违背本发明精神的前提下还可做出等同变形或替换，这些等同的变型或替换均包含在本 申请权利要求所限定的范围内。