CN105487524A - 具有多工况特性的超高速小盒包装机状态监测与诊断方法 - Google Patents

Abstract

本发明公开了一种具有多工况特性的超高速小盒包装机状态监测与诊断方法，该方法针对运行过程多工况特性及工况普遍未知的问题，基于滑动窗口计算稳定度因子识别不同工况，采用自适应聚类方法对稳定工况进行划分形成若干稳定工况数据簇，利用PCA方法为每个稳定工况数据簇建立监测模型。多工况过程状态监测时，根据当前滑动窗口内数据的稳定度因子判断工况类型，若为过渡工况，则将监测统计量赋值为0；若为稳定工况，计算当前有效数据与各个聚类中心的距离，获得匹配的PCA监测模型进行实时监测，任一统计量超限时采用贡献图确定引起故障的主要过程变量。

Description

具有多工况特性的超高速小盒包装机状态监测与诊断方法

技术领域

本发明涉及意大利G.D公司超高速小盒包装机H1000的统计监测建模、状态监测与诊断技术。

背景技术

随着近年来烟草行业全球化、新的法律法规以及外部环境的重大变化，中国烟草面临愈来愈严峻的竞争压力和社会压力。在保证产品质量的前提下，有效提升设备的智能化水平和高效运行能力成为卷烟工厂关注的焦点。智能化，即生产设备具有感知、分析、推理和决策等功能；高效运行，即在规定时间内消耗更少的资源生产更多、更好的产品。卷烟工厂设备的智能化水平和高效运行能力与其自动化程度密不可分，设备状态监测与诊断技术是提升智能化水平、保证高效运行的重要手段。

意大利G.D公司超高速包装设备是中国烟草总公司“十二五”期间重点引进的目前世界上最高速、最先进的卷烟包装设备，其设计生产能力为1000包香烟/分钟。超高速包装设备由小盒包装机H1000、小盒薄膜封装机W1000和条盒包装机BV三部分组成。小盒包装机H1000是一台高速连续运动硬盒翻盖包装机，利用烟支、内衬纸、内框纸、商标纸将烟支按7-6-7共20支的形式进行翻盖式小盒包装。小盒包装机H1000是超高速包装设备最关键、最核心的单元，也是生产现场故障最集中的部位。由于G.D公司超高速包装设备底层PLC开放的程序源代码和库文件非常有限，生产现场主要通过每班一组的数据报表、生产时的人工目测、轮保时的人工检测、设备故障的计划维修和事后抢修等方式保证小盒包装机H1000的设备完好性。杭州卷烟厂2014年超高速包装设备的运行效率仅为70.08％，目前采用的设备状态监测方式已经不能满足卷烟生产的实际需要，亟需一种切实有效的超高速小盒包装机H1000状态监测与诊断方法。

为解决G.D公司超高速包装设备底层PLC参数数据难以获取的问题，杭州卷烟厂(朱立明,钱杰,王琪,李钰靓,李勇.一种超高速包装机数据采集系统:中国,CN204056360U[P].2014-12-31.)提出了一种基于Beckoff软PLC的超高速包装设备数据采集方法，通过对设备底层PLC源代码和库文件的分析，能够采集超高速小盒包装机H1000的部分参数数据，从而为基于多元统计分析方法实施具有多工况特性的超高速小盒包装机H1000状态监测与诊断提供了丰富的运行数据。

发明内容

超高速小盒包装机H1000生产过程中，操作人员根据目测的设备运行状态频繁调整车速设定值，使得生产过程具有多个运行工况。由于小盒包装机H1000的高速运行，从一个稳定运行工况到另外一个稳定运行工况的时间非常短(如从车速为940包/分钟降低到50包/分钟仅需20秒，即在两个采样周期内完成)，因此本发明仅考虑各个不同的稳定运行工况。不同稳定运行工况之间数据的均值、方差、相关关系等特征变量会有明显的变化，直接应用PCA方法进行统计建模和状态监测会产生大量的误报警或漏报警。

为满足超高速小盒包装机H1000状态监测与诊断的实际需要，本发明面向设备生产过程的多工况特性，提出一种包含离线工况识别和划分、稳定工况统计监测建模、多工况过程状态监测、故障原因变量分离等关键技术的超高速小盒包装机H1000状态监测与诊断方法。该方法基于滑动时间窗口离线计算稳定度因子识别稳定工况和过渡工况，采用自适应k-means聚类方法对稳定工况进行划分形成若干个稳定工况数据簇，再利用PCA方法对每个稳定工况数据簇独立建立统计监测模型。多工况过程状态监测时，根据当前滑动时间窗口内数据的稳定度因子判断工况类型，若为过渡工况，则将当前时刻T²、SPE两个监测统计量赋值为0；若为稳定工况，计算当前滑动时间窗口内有效数据与各个聚类中心的欧式距离，利用距离最小的聚类所对应的PCA监测模型进行实时监测，任一统计量超限时利用贡献图方法确定引起故障的主要过程变量，最终为超高速小盒包装机H1000的状态监测与诊断提供了一种有效可行的方法。

本发明所采用的具体技术方案如下：

一种具有多工况特性的超高速小盒包装机H1000状态监测与诊断方法，包括：

1)根据超高速小盒包装机H1000的运行机理确定统计模型的输入变量，采集某一机台底层PLC中的设备运行数据，获得二维原始数据D(K×J)，K为采样点个数，J为监测变量个数；

2)基于车速检测值，采用滑动时间窗口H方法计算稳定度因子SF，当其值较小时(小于等于阈值ζ)表明过程处于稳定工况，反之表明过程处于过渡工况。采用平滑滤波方法对处于稳定工况的滑动时间窗口H内原始数据进行算术平均，获得一组移动平均数据，将移动平均数据中“车速”变量值大于阈值ψ的数据作为稳定工况的有效数据D＝{D₁,D₂,…,D_N}，N为有效数据的个数；

3)稳定工况的有效数据被识别出以后，采用自适应k-means聚类方法对有效数据的相似度进行聚类分析，获得每一种稳定工况的建模数据C表示将稳定工况划分成C个聚类，每个聚类中具有相似过程特性的建模数据可以用一个统计监测模型来描述；

4)对每一种稳定工况的建模数据(r＝1,2,…,C)进行均值中心化和方差归一化的数据预处理，获得处理后的建模数据X_r(K_rm×J)。对建模数据X_r进行PCA分解，得到对应的主元个数A_r、负载矩阵P_r、控制限和建立对应第r个稳定工况的PCA监测模型；

5)多工况过程状态监测时，计算当前滑动时间窗口内数据的稳定度因子SF_new，通过与给定阈值ζ的比较确定当前所处的工况类型。若为过渡工况，则将当前时刻和SPE_new两个监测统计量赋值为0；否则采用平滑滤波方法计算当前滑动时间窗口内的移动平均数据，并计算其与各个聚类中心的欧式距离，调用距离最小的聚类所对应的PCA监测模型在线计算和SPE_new两个监测统计量。如果两个统计量都位于控制限以内，表明设备状态正常，如果至少其中一个统计量超出控制限，表明设备状态异常；

6)当检测到设备出现有异常时，根据超限统计量处于T²的主元子空间或是SPE的残差子空间，计算每个过程变量对超限统计量的贡献率，贡献率较大的变量被初步确定为造成设备异常的原因变量。

所述步骤2)中，采用滑动时间窗口H方法计算稳定度因子SF，计算公式如下：

$SF=\frac{V_{Hmax}-V_{H\min }}{V_{max}}$

其中V_Hmax表示滑动时间窗口H内车速数据的最大值，V_Hmin表示滑动窗口H内车速数据的最小值，V_max表示设备允许的车速最大值。H表示滑动时间窗口的长度。

在步骤3)中，采用自适应k-means聚类方法对有效数据的相似度进行聚类分析，通过定义两个聚类中心的最小距离的阈值θ，同时最小化同一聚类数据之间的局部距离平方和以及所有聚类之间的全局距离平方和，获得最优的聚类数目和聚类中心，实现对稳定工况的自适应划分。

在步骤4)中，对建模数据预处理包括依次进行的减均值、除标准差处理获得X_r(K_rm×J)，计算公式如下：

$x_{rk,j}=\frac{{\tilde{x}}_{rk,j}-{\stackrel{\‾}{x}}_{r,j}}{{\tilde{s}}_{r,j}},rk=1,...K_{rm},;j=1,...,J$

其中表示建模数据中第rk行第j列的数据，下标r代表工况、k代表样本、j代表变量，表示第r个稳定工况第j个变量的均值，表示第r个稳定工况第j个变量的标准差。

将数据X_r进行PCA分解，PCA分解的计算公式为：

$X_r=T_r{P_r}^T+E_r=\underset{a=1}{\overset{A_r}{\Σ}}{t}_{ra}{p}_{ra}^T+E_r$

其中T_r为主元子空间中(K_rm×A_r)维的得分矩阵；P_r为主元子空间中(J×A_r)维的负载矩阵；E_r为残差子空间中(K_rm×J)维的残差矩阵。A_r表示主元个数，由累计贡献率方法确定。

离线计算主成分子空间的Hotelling-T²和残差子空间的SPE(SquaredPredictionError)监测统计量。监测统计量的计算公式为：

$T_{rk}^2=t_{rk}{S_r}^{-1}{t_{rk}}^T,k=1,2,...,K_{rm}$

其中t_rk＝x_rkP_r为(1×A_r)维的主元得分向量，x_rk为建模数据X_r的第k行，P_r为主元子空间的负载矩阵，对角矩阵由建模数据X_r的协方差矩阵∑X_r ^TX_r的前A_r个特征值所构成。

SPE_r监测统计量的计算公式为：

${\mathrm{SPE}}_{rk}=e_{rk}{e_{rk}}^T=\underset{j=1}{\overset{J}{\Σ}}{(x_{rk,j}-{\hat{x}}_{rk,j})}^2,k=1,2,...,K_{rm}$

其中表示重构得到的(1×J)维估计向量。

根据F分布计算监测统计量的控制限为：

${\mathrm{ctr}}_{T_r^2}=\frac{A_r(K_{rm}-1)}{K_{rm}-A_r}{F}_{\mathrm{\α}}(A_r,K_{rm}-A_r)$

其中α为置信度，F_α(A_r,K_rm-A_r)是对应于置信度为α，自由度为A_r、K_rm-A_r条件下的F分布临界值。

根据χ²分布计算SPE_r监测统计量的控制限为：

${\mathrm{ctr}}_{{\mathrm{SPE}}_r}={\mathrm{g\χ}}_{h,\mathrm{\α}}^2$

其中g＝v_r/(2n_r)，n_r、v_r分别为SPE_r监测统计量的均值和方差。

在步骤5)中，计算处于稳定工况滑动时间窗口的移动平均数据与各个聚类中心的欧式距离，根据距离最小判断该时刻对应的稳定工况监测模型。欧式距离的计算公式为：

${\mathrm{dist}}_{new,i}=\underset{j=1}{\overset{J}{\Σ}}(x_{new,j}-{\mathrm{Cen}}_{i,j}),i=1,2,...,C$

其中x_new,j表示移动平均数据的第j个变量的检测值，Cen_i,j表示第i个聚类中心第j个变量的值，C表示聚类中心的个数。

设当前移动平均数据对应第L(1≤L≤C)个稳定工况监测模型，利用该监测模型建模数据的均值和标准差对移动平均数据进行标准化处理，获得x_Lnew(1×J)。将其投影到第L个监测模型上，计算和SPE_Lnew监测统计量指标。

$\left\{\begin{array}{c}{T}_{Lnew}^2=x_{Lnew}{P}_L{S_L}^{-1}{P_L}^T{x_{Lnew}}^T\\ {\mathrm{SPE}}_{Lnew}=e_{Lnew}{e_{Lnew}}^T=\underset{j=1}{\overset{J}{\Σ}}{(x_{Lnew,j}-{\hat{x}}_{Lnew,j})}^2\end{array}\right.$

其中P_L表示第L个监测模型主元子空间的负载矩阵，对角矩阵S_L是由建模数据X_L的协方差矩阵的前A_L个特征值所构成，表示重构得到的(1×J)维估计向量。

如果两个监测统计量都位于控制限内，表明超高速小盒包装机处于正常状态，如果其中至少一个监测统计量超出控制限，表明设备处于异常状态。

在步骤6)中，当任一监测统计量超出控制限时，计算各个过程变量对超限统计量的贡献值，确定引起异常的原因变量。主成分子空间监测统计量的超出正常控制限时，主成分得分变量t_La对的贡献率可以计算如下：

$C_{T_{Lnew}^2,t_{La}}=\frac{t_{La}^2}{{\mathrm{\λ}}_{La}}/T_{Lnew}^2,a=1,...,A_L$

其中t_La表示对应第L个监测模型的第a个主成分得分变量，λ_La表示对应第L个监测模型的第a个特征值。

过程变量x_Lnew,j对t_La的贡献率可以计算如下：

$C_{t_{La},x_{Lnew,j}}=x_{Lnew,j}{p}_{La,j}/t_{La},a=1,\mathrm{...},A_L;j=1,\mathrm{...},J$

其中p_La,j表示对应第L个监测模型的第a个负载向量的第j个变量。

当残差子空间SPE_Lnew监测统计量的超出正常控制限时，过程变量x_Lnew,j对SPE_Lnew的贡献率可以计算如下：

$C_{{\mathrm{SPE}}_{Lnew},x_{Lnew,j}}=sign(x_{Lnew,j}-{\hat{x}}_{Lnew,j})\·\frac{{(x_{Lnew,j}-{\hat{x}}_{Lnew,j})}^2}{{\mathrm{SPE}}_{Lnew}}$

其中表示重构估计向量的第j个变量，表示残差的正负信息。

本发明的有益效果是：

本发明主要涉及离线工况识别和划分、稳定工况统计监测建模、多工况过程状态监测、故障原因变量分离等关键技术。基于海量历史数据，实现了多工况过程稳定工况、过渡工况的有效识别以及稳定工况的自适应划分，进而自动完成稳定工况的离线监测建模及控制限的求取。多工况过程状态监测时，匹配当前稳定工况的移动平均数据与各个聚类中心的相似度判断该时刻对应的稳定工况监测模型，对当前的生产过程进行有效监测，监测异常发生时采用变量贡献图追溯和确定引起故障的原因变量。最终形成一套集工况识别、统计建模、状态监测和故障诊断于一体的切实可行的实施方案，能有效提高设备诊断结果的可靠性和准确性，从而为超高速小盒包装机H1000的保养和维护提供科学决策和有效指导。

附图说明

图1为本发明的具有多工况特性超高速小盒包装机H1000状态监测与诊断流程图；

图2为工况识别前1～1000采样点的H1000车速曲线图；

图3为1～1000采样点H1000车速的稳定度因子SF曲线图；

图4为经过平滑滤波后稳定工况的H1000车速曲线图；

图5为稳定工况有效数据的H1000车速曲线图；

图6为多工况过程正常状态测试数据的T²和SPE过程监测图；

图7为多工况过程故障状态数据1的T²和SPE过程监测图；

图8为多工况过程故障状态数据2的T²和SPE过程监测图；

图9为多工况过程故障状态数据3的T²和SPE过程监测图；

图10为故障状态1下首次故障报警时刻的T²变量贡献图；

图11为故障状态2下首次故障报警时刻的T²变量贡献图；

图12为故障状态3下首次故障报警时刻的T²变量贡献图。

具体实施方式

为了更好的理解本发明的技术方案，以下结合说明书附图对本发明的实施方式作进一步描述。本实施是用于意大利G.D公司超高速小盒包装机H1000的状态监测与诊断方法。G.D公司超高速包装设备由小盒包装机H1000、小盒薄膜封装机W1000和条盒包装机BV三部分组成。小盒包装机H1000是一台高速连续运动硬盒翻盖包装机，其运行以间歇式往复运动和公转加自传的行星驱动为主，利用烟支、内衬纸、内框纸、商标纸将烟支按7-6-7共20支的形式进行翻盖式小盒包装。本发明的具有多工况特性的超高速小盒包装机状态监测与诊断方法的实现框图如图1所示，方法主要分为以下几步：

(1)历史运行数据的获取

根据超高速小盒包装机H1000的运行机理确定统计模型的输入变量，采集某一机台底层PLC中的设备运行数据，获得二维原始数据D(K×J)，K为采样点个数，J为监测变量个数。

本实例中，选取A4机台的超高速小盒包装机，采集的设备运行数据中包含以下变量：伺服电机电流、伺服电机温度、热熔胶温度、车速等44个监测变量，如表1所示。监测变量的采样频率为10秒/次，获得的原始数据位D(64394×44)。

表1超高速小盒包装机的监测变量表

(2)离线工况的识别

基于车速检测值，采用滑动时间窗口H方法计算稳定度因子SF，当其值较小时(小于等于阈值ζ)表明过程处于稳定工况，反之表明过程处于过渡工况。采用平滑滤波方法对处于稳定工况的滑动时间窗口H内原始数据进行算术平均，获得一组移动平均数据，将移动平均数据中“车速”变量值大于阈值ψ的数据作为稳定工况的有效数据D＝{D₁,D₂,…,D_N}，N为有效数据的个数。

采用滑动时间窗口H方法计算稳定度因子SF，计算公式如下：

$SF=\frac{V_{Hmax}-V_{Hmin}}{V_{max}}$

其中V_Hmax表示滑动时间窗口H内车速数据的最大值，V_Hmin表示滑动窗口H内车速数据的最小值，V_max表示设备允许的车速最大值。H表示滑动时间窗口的长度。

本实例中，滑动时间窗口的长度H取为6，允许的车速最大值V_max取为1000，阈值ζ取为0.02，阈值ψ取为800，获得的有效数据为D(30283×44)。工况识别前1～1000采样点的H1000车速如图2所示，相应的稳定度因子如图3所示，经过平滑滤波后稳定工况的H1000车速如图4所示，稳定工况有效数据的H1000车速如图5所示。

(3)稳定工况的划分

稳定工况的有效数据被识别出以后，采用自适应k-means聚类方法对有效数据的相似度进行聚类分析，获得每一种稳定工况的建模数据C表示将稳定工况划分成C个聚类，每个聚类中具有相似过程特性的建模数据可以用一个统计模型来描述。自适应k-means聚类方法的求取过程如下：

步骤1：设定初始聚类中心个数C₀，两个聚类中心的最小距离阈值θ，每个聚类俘获的最小数据量阈值minL。初始化C₀个聚类中心W_i,c(i＝1,…,J；c＝1,…,C₀)；

步骤2：计算两个聚类中心的欧式距离dist(W_i,c1，W_i,c2)，若其小于预定的阈值θ，剔除其中一个聚类中心；

步骤3：计算每个稳定工况的有效数据D_k(k＝1,2,…,N)到所有聚类中心的欧式距离dist(D_k,W_i,c)，若D_k与第c^*类的中心的距离最小，则将D_k的隶属关系定义为m(k)＝c^*；

步骤4：I_num次迭代后，若某聚类中心俘获的数据样本个数小于预定的阈值minL，则剔除该类。

步骤5：更新聚类个数为C_i+1，并根据数据的隶属关系重新计算新的聚类中心W_i+1,c(c＝1,2,…,C_i+1)。

步骤6：如果算法满足收敛条件则结束，否则返回步骤2，进行下一次迭代计算。收敛条件有：两次迭代中的聚类中心距离的变化小于一个很小的阈值ε，或者每个聚类中数据D_k到聚类中心的距离平方和以及聚类之间的距离平方和达到最小。

本实例中，初始聚类中心个数C₀取为15，最小距离阈值θ取为300，最小数据量阈值minL取为100，阈值ε取为0.0001。通过自适应聚类，稳定工况有效数据被划分为7个聚类(即C＝7)，获得每一种稳定工况的建模数据

(4)稳定工况统计监测建模

对每一种稳定工况的建模数据(r＝1,2,…,C)进行均值中心化和方差归一化的数据预处理，获得处理后的建模数据X_r(K_rm×J)。对建模数据X_r进行PCA分解，得到对应的主元个数A_r、负载矩阵P_r、控制限和建立对应第r个稳定工况的PCA监测模型。

对建模数据预处理包括依次进行的减均值、除标准差处理获得X_r(K_rm×J)，计算公式如下：

$x_{rk,j}=\frac{{\tilde{x}}_{rk,j}-{\stackrel{\‾}{x}}_{r,j}}{{\tilde{s}}_{r,j}},rk=1,...K_{rm},;j=1,...,J$

其中表示建模数据中第rk行第j列的数据，下标r代表工况数、k代表样本数、j代表变量数，表示第r个稳定工况第j个变量的均值，表示第r个稳定工况第j个变量的标准差。

将数据X_r进行PCA分解，PCA分解的计算公式为：

$X_r=T_r{P_r}^T+E_r=\underset{a=1}{\overset{A_r}{\Σ}}{t}_{ra}{p}_{ra}^T+E_r$

其中T_r为主元子空间中(K_rm×A_r)维的得分矩阵；P_r为主元子空间中(J×A_r)维的负载矩阵；E_r为残差子空间中(K_rm×J)维的残差矩阵。A_r表示主元个数，由累计贡献率方法确定。

离线计算主成分子空间的Hotelling-T²和残差子空间的SPE(SquaredPredictionError)监测统计量。监测统计量的计算公式为：

$T_{rk}^2=t_{rk}{S_r}^{-1}{t_{rk}}^T,k=1,2,\mathrm{...},K_{rm}$

其中t_rk＝x_rkP_r为(1×A_r)维的主元得分向量，x_rk为建模数据X_r的第k行，P_r为主元子空间的负载矩阵，对角矩阵由建模数据X_r的协方差矩阵∑X_r ^TX_r的前A_r个特征值所构成。

SPE_r监测统计量的计算公式为：

${\mathrm{SPE}}_{rk}=e_{rk}{e_{rk}}^T=\underset{j=1}{\overset{J}{\Σ}}{(x_{rk,j}-{\hat{x}}_{rk,j})}^2,k=1,2,...,K_{rm}$

其中表示重构得到的(1×J)维估计向量。

根据F分布计算监测统计量的控制限为：

${\mathrm{ctr}}_{T_r^2}=\frac{A_r(K_{rm}-1)}{K_{rm}-A_r}{F}_{\mathrm{\α}}(A_r,K_{rm}-A_r)$

其中α为置信度，F_α(A_r,K_rm-A_r)是对应于置信度为α，自由度为A_r、K_rm-A_r条件下的F分布临界值。

根据χ²分布计算SPE_r监测统计量的控制限为：

${\mathrm{ctr}}_{{\mathrm{SPE}}_r}={\mathrm{g\χ}}_{h,\mathrm{\α}}^2$

其中g＝v_r/(2n_r)，n_r、v_r分别为SPE_r监测统计量的均值和方差。

本实例中，7个PCA监测模型的主元个数依次为A₁＝26、A₂＝21、A₃＝26、A₄＝26、A₅＝27、A₆＝26、A₇＝26，则P₁为(44×26)负载矩阵、P₂为(44×21)负载矩阵、P₃为(44×26)负载矩阵、P₄为(44×26)负载矩阵、P₅为(44×27)负载矩阵、P₆为(44×26)负载矩阵、P₇为(44×26)负载矩阵。统计控制限的置信度α取为0.99，监测模型1的控制限和监测模型2的控制限和监测模型3的控制限和监测模型4的控制限和监测模型5的控制限和监测模型6的控制限和监测模型7的控制限和

(5)多工况过程状态监测

多工况过程状态监测时，计算当前滑动时间窗口内数据的稳定度因子SF_new，通过与给定阈值ζ的比较确定当前所处的工况类型。若为过渡工况，则将当前时刻和SPE_new两个监测统计量赋值为0；否则采用平滑滤波方法计算当前滑动时间窗口内的移动平均数据，并计算其与各个聚类中心的欧式距离，调用距离最小的聚类所对应的PCA监测模型在线计算和SPE_new两个监测统计量。如果两个统计量都位于控制限以内，表明设备状态正常，如果至少其中一个统计量超出控制限，表明设备状态异常。

计算处于稳定工况滑动时间窗口的移动平均数据与各个聚类中心的欧式距离，根据距离最小判断该时刻对应的稳定工况监测模型。欧式距离的计算公式为：

${\mathrm{dist}}_{new,i}=\underset{j=1}{\overset{J}{\Σ}}(x_{new,j}-{\mathrm{Cen}}_{i,j}),i=1,2,...,C$

其中x_new,j表示移动平均数据的第j个变量的检测值，Cen_i,j表示第i个聚类中心第j个变量的值，C表示聚类中心的个数。

设当前移动平均数据对应第L(1≤L≤C)个稳定工况监测模型，利用该监测模型建模数据的均值和标准差对移动平均数据进行标准化处理，获得x_Lnew(1×J)。将其投影到第L个监测模型上，计算和SPE_Lnew监测统计量指标。

$\left\{\begin{array}{c}{T}_{Lnew}^2=x_{Lnew}{P}_L{S_L}^{-1}{P_L}^T{x_{Lnew}}^T\\ {\mathrm{SPE}}_{Lnew}=e_{Lnew}{e_{Lnew}}^T=\underset{j=1}{\overset{J}{\Σ}}{(x_{Lnew,j}-{\hat{x}}_{Lnew,j})}^2\end{array}\right.$

其中P_L表示第L个监测模型主元子空间的负载矩阵，对角矩阵S_L是由建模数据X_L的协方差矩阵的前A_L个特征值所构成，表示重构得到的(1×J)维估计向量。

如果两个监测统计量都位于控制限内，表明超高速小盒包装机处于正常状态，如果其中至少一个监测统计量超出控制限，表明设备处于异常状态。

本实例中，选取多工况过程8816组采样数据D _new(8816×44)进行分析，该测试数据T²和SPE统计量的监测结果如图6所示。引入首次故障报警时间，定义为连续9个监测统计量超过对应控制限的采样时间点。在两个监测子空间中有个别样本点的监测统计量超出控制限，根据首次故障报警的定义，两个监测统计量均未发生故障报警，对当前过程的正常状态做出了正确的指示，表明所提状态监测与诊断方法具有准确监测多工况过程正常状态的能力。

作为比较，重新选取多工况过程4505组采样数据D _new(4505×44)进行状态监测与诊断的分析。在第1001个样本点的时候，分别令热熔胶温度阶跃增长1％(突变故障)、A1电机电流阶跃降低30％(突变故障)、A5A电机温度以0.03的斜率增加(渐变故障)，一直持续到第2000个样本点的时候恢复正常，形成故障状态数据1、2和3。故障状态数据T²和SPE统计量的监测结果分别如图7至图9所示。

对于故障1，主元子空间的T²统计量出现连续超限，主成分子空间出现异常，检测出的首次故障报警时间为第1030个样本点时刻，即监测模型延后290秒准确检测出故障的发生；对于故障2，主元子空间的T²统计量出现连续超限，主成分子空间出现异常，检测出的首次故障报警时间为第1030个样本点时刻，即监测模型延后290秒准确检测出故障的发生；对于故障3，主元子空间的T²统计量出现连续超限，主成分子空间出现异常，检测出的首次故障报警时间为第1318个样本点时刻，即监测模型延后3170秒准确检测出故障的发生。说明所提方法能对多工况过程的突变和渐变故障做出正确指示，表明该方法具有准确监测多工况过程故障状态的能力。

(6)基于贡献图的故障诊断

当任一监测统计量超出控制限时，计算各个过程变量对超限统计量的贡献值，确定引起异常的原因变量。主成分子空间监测统计量的超出正常控制限时，主成分得分变量t_La对的贡献率可以计算如下：

$C_{T_{Lnew}^2,t_{La}}=\frac{t_{La}^2}{{\mathrm{\λ}}_{La}}/T_{Lnew}^2,a=1,...,A_L$

其中t_La表示对应第L个监测模型的第a个主成分得分变量，λ_La表示对应第L个监测模型的第a个特征值。

过程变量x_Lnew,j对t_La的贡献率可以计算如下：

$C_{t_{La},x_{Lnew,j}}=x_{Lnew,j}{p}_{La,j}/t_{La},a=1,\mathrm{...},A_L;j=1,\mathrm{...},J$

其中p_La,j表示对应第L个监测模型的第a个负载向量的第j个变量。

当残差子空间SPE_Lnew监测统计量的超出正常控制限时，过程变量x_Lnew,j对SPE_Lnew的贡献率可以计算如下：

$C_{{\mathrm{SPE}}_{Lnew},x_{Lnew,j}}=sign(x_{Lnew,j}-{\hat{x}}_{Lnew,j})\·\frac{{(x_{Lnew,j}-{\hat{x}}_{Lnew,j})}^2}{{\mathrm{SPE}}_{Lnew}}$

其中表示重构估计向量的第j个变量，表示残差的正负信息。

本实例中，根据图7的监测结果，故障发生在主成分子空间中。在首次故障报警的第1030个样本点时刻，过程变量对监测统计量超限的贡献率如图10所示，确定引起故障的原因变量为第43个变量，即H1000热熔胶温度。

根据图8的监测结果，故障发生在主成分子空间中。在首次故障报警的第1030个样本点时刻，过程变量对监测统计量超限的贡献率如图11所示，确定引起故障的原因变量为第1个变量，即A1电机电流。

根据图9的监测结果，故障发生在主成分子空间中。在首次故障报警的第1318个样本点时刻，过程变量对监测统计量超限的贡献率如图12所示，确定引起故障的原因变量为第16个变量，即A5A电机温度。

以上所述仅为本发明的较佳实施举例，并不用于限制本发明，凡在本发明精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种具有多工况特性的超高速小盒包装机状态监测与诊断方法，其特征在于，包括步骤：

1)根据超高速小盒包装机H1000的运行机理确定统计模型的输入变量，采集相应的运行数据，获得二维原始数据K为采样点个数，J为监测变量个数；

2)基于车速检测值，采用滑动时间窗口H方法计算稳定度因子SF，根据稳定度因子SF与阈值ζ的关系判断是否处于稳定工况，采用平滑滤波方法对处于稳定工况的滑动时间窗口H内原始数据进行算术平均，获得一组移动平均数据，将移动平均数据中“车速”变量值大于阈值ψ的数据作为稳定工况的有效数据D＝{D₁,D₂,…,D_N}，N为有效数据的个数；

3)采用自适应k-means聚类方法对有效数据的相似度进行聚类分析，获得每一种稳定工况的建模数据C表示将稳定工况划分成C个聚类，每个聚类中具有相似过程特性的建模数据用一统计监测模型来描述；

4)对每一种稳定工况的建模数据进行均值中心化和方差归一化的数据预处理，获得处理后的建模数据X_r(K_rm×J)，对建模数据X_r进行PCA分解，得到对应的主元个数A_r、负载矩阵P_r、控制限和建立对应第r个稳定工况的PCA监测模型；

5)实时计算当前滑动时间窗口内数据的稳定度因子SF_new，通过与给定阈值ζ的比较确定当前所处的工况类型；若为过渡工况，则将当前时刻和SPE_new两个监测统计量赋值为0；否则采用平滑滤波方法计算当前滑动时间窗口内的移动平均数据，并计算其与各个聚类中心的欧式距离，调用距离最小的聚类所对应的PCA监测模型在线计算和SPE_new两个监测统计量，若一个监测统计量超出控制限，表明设备状态异常；

6)当检测到设备出现有异常时，根据超限的监测统计量处于T²的主元子空间或是SPE的残差子空间，计算每个过程变量对超限统计量的贡献率，贡献率较大的变量被初步确定为造成设备异常的原因变量。

2.如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其特征在于，在步骤2)中采用滑动时间窗口H方法计算稳定度因子SF进行工况的识别，计算公式如下：

$SF=\frac{V_{Hmax}-V_{Hmin}}{V_{max}}$

其中V_Hmax表示滑动时间窗口H内车速数据的最大值，V_Hmin表示滑动窗口H内车速数据的最小值，V_max表示设备允许的车速最大值，H表示滑动时间窗口的长度。

3.如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其特征在于，在步骤3)中采用自适应k-means聚类方法进行工况的划分，对有效数据的相似度进行聚类分析，通过定义两个聚类中心的最小距离的阈值θ，同时最小化同一聚类数据之间的局部距离平方和以及所有聚类之间的全局距离平方和，获得最优的聚类数目和聚类中心，实现对稳定工况的自适应划分。

4.如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其特征在于，在步骤4)中对建模数据预处理包括依次进行的减均值、除标准差处理获得X_r(K_rm×J)，计算公式如下：

$x_{rk,j}=\frac{{\tilde{x}}_{rk,j}-{\stackrel{\‾}{x}}_{r,j}}{{\tilde{s}}_{r,j}},rk=1,...K_{rm},;j=1,...,J$

其中，表示建模数据中第rk行第j列的数据，下标r代表工况、k代表样本、j代表变量，表示第r个稳定工况第j个变量的均值，表示第r个稳定工况第j个变量的标准差。

5.如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其特征在于，所述步骤4)中PCA分解的计算公式为：

$X_r=T_r{P_r}^T+E_r=\underset{a=1}{\overset{A_r}{\Σ}}{t}_{ra}{p}_{ra}^T+E_r$

其中，T_r为主元子空间中(K_rm×A_r)维的得分矩阵；P_r为主元子空间中(J×A_r)维的负载矩阵；E_r为残差子空间中(K_rm×J)维的残差矩阵；A_r表示主元个数，由累计贡献率方法确定；

离线计算主成分子空间的Hotelling-T²和残差子空间的SPE(SquaredPredictionError)监测统计量；监测统计量的计算公式为：

$T_{rk}^2=t_{rk}{S_r}^{-1}{t_{rk}}^T,k=1,2,...,K_{rm}$

其中t_rk＝x_rkP_r为(1×A_r)维的主元得分向量，x_rk为建模数据X_r的第k行，P_r为主元子空间的负载矩阵，对角矩阵由建模数据X_r的协方差矩阵的前A_r个特征值所构成；

SPE_r监测统计量的计算公式为：

${\mathrm{SPE}}_{rk}=e_{rk}{e_{rk}}^T=\underset{j=1}{\overset{J}{\Σ}}{(x_{rk,j}-{\hat{x}}_{rk,j})}^2,k=1,2,...,K_{rm}$

其中表示重构得到的(1×J)维估计向量；

根据F分布计算监测统计量的控制限为：

${\mathrm{ctr}}_{T_r^2}=\frac{A_r(K_{rm}-1)}{K_{rm}-A_r}{F}_{\mathrm{\α}}(A_r,K_{rm}-A_r)$

其中α为置信度，F_α(A_r,K_rm-A_r)是对应于置信度为α，自由度为A_r、K_rm-A_r条件下的F分布临界值；

根据χ²分布计算SPE_r监测统计量的控制限

${\mathrm{ctr}}_{{\mathrm{SPE}}_r}={\mathrm{g\χ}}_{h,\mathrm{\α}}^2$

其中g＝v_r/(2n_r)，分别为SPE_r监测统计量的均值和方差。

6.如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其特征在于，在步骤5)中，欧式距离的计算公式为：

${\mathrm{dist}}_{new,i}=\underset{j=1}{\overset{J}{\Σ}}(x_{new,j}-{\mathrm{Cen}}_{i,j}),i=1,2,...,C$

其中x_new,j表示移动平均数据的第j个变量的检测值，Cen_i,j表示第i个聚类中心第j个变量的值，C表示聚类中心的个数。

7.如权利要求6所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其特征在于，设当前移动平均数据对应第L(1≤L≤C)个稳定工况监测模型，利用该监测模型建模数据的均值和标准差对移动平均数据进行标准化处理，获得x_Lnew(1×J)，将其投影到第L个监测模型上，计算和SPE_Lnew监测统计量指标；

$\left\{\begin{array}{c}{T}_{Lnew}^2=x_{Lnew}{P}_L{S_L}^{-1}{P_L}^T{x_{Lnew}}^T\\ {\mathrm{SPE}}_{Lnew}=e_{Lnew}{e_{Lnew}}^T=\underset{j=1}{\overset{J}{\Σ}}{\left(x_{Lnew,j}-{\hat{x}}_{Lnew,j}\right)}^2\end{array}\right.$

其中P_L表示第L个监测模型主元子空间的负载矩阵，对角矩阵S_L是由建模数据X_L的协方差矩阵的前A_L个特征值所构成，表示重构得到的(1×J)维估计向量。

8.如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其特征在于，在步骤6)中，主成分子空间监测统计量的超出正常控制限时，主成分得分变量t_La对的贡献率可以计算如下：

$C_{T_{Lnew,}^2{t}_{La}}=\frac{t_{La}^2}{{\mathrm{\λ}}_{La}}/T_{Lnew}^2,a=1,\mathrm{...},A_L$

其中t_La表示对应第L个监测模型的第a个主成分得分变量，λ_La表示对应第L个监测模型的第a个特征值；

过程变量x_Lnew,j对t_La的贡献率可以计算如下：

$C_{t_{La},x_{Lnew,j}}=x_{Lnew,j}{p}_{La,j}/t_{La},a=1,...,A_L;j=1,...,J$

其中p_La,j表示对应第L个监测模型的第a个负载向量的第j个变量；

当残差子空间SPE_Lnew监测统计量的超出正常控制限时，过程变量x_Lnew,j对SPE_Lnew的贡献率可以计算如下：

$C_{{\mathrm{SPE}}_{Lnew,}{x}_{Lnew,j}}=sign(x_{Lnew,j}-{\hat{x}}_{Lnew,j})\·\frac{{(x_{Lnew,j}-{\hat{x}}_{Lnew,j})}^2}{{\mathrm{SPE}}_{Lnew}}$

其中表示重构估计向量的第j个变量，表示残差的正负信息。