具有多工况特性的超高速小盒包装机状态监测与诊断方法
具有多工况特性的超高速小盒包装机状态监测与诊断方法
【技术领域】
[0001] 本发明涉及意大利G.D公司超高速小盒包装机H1000的统计监测建模、状态监测与 诊断技术。
【背景技术】
[0002] 随着近年来烟草行业全球化、新的法律法规以及外部环境的重大变化，中国烟草 面临愈来愈严峻的竞争压力和社会压力。在保证产品质量的前提下，有效提升设备的智能 化水平和高效运行能力成为卷烟工厂关注的焦点。智能化，即生产设备具有感知、分析、推 理和决策等功能；高效运行，即在规定时间内消耗更少的资源生产更多、更好的产品。卷烟 工厂设备的智能化水平和高效运行能力与其自动化程度密不可分，设备状态监测与诊断技 术是提升智能化水平、保证高效运行的重要手段。
[0003] 意大利G.D公司超高速包装设备是中国烟草总公司"十二五"期间重点引进的目前 世界上最高速、最先进的卷烟包装设备，其设计生产能力为1〇〇〇包香烟/分钟。超高速包装 设备由小盒包装机H1000、小盒薄膜封装机W1000和条盒包装机BV三部分组成。小盒包装机 H1000是一台高速连续运动硬盒翻盖包装机，利用烟支、内衬纸、内框纸、商标纸将烟支按7-6-7共20支的形式进行翻盖式小盒包装。小盒包装机H1000是超高速包装设备最关键、最核 心的单元，也是生产现场故障最集中的部位。由于G.D公司超高速包装设备底层PLC开放的 程序源代码和库文件非常有限，生产现场主要通过每班一组的数据报表、生产时的人工目 测、轮保时的人工检测、设备故障的计划维修和事后抢修等方式保证小盒包装机H1000的设 备完好性。杭州卷烟厂2014年超高速包装设备的运行效率仅为70.08%，目前采用的设备状 态监测方式已经不能满足卷烟生产的实际需要，亟需一种切实有效的超高速小盒包装机 H1000状态监测与诊断方法。
[0004] 为解决G.D公司超高速包装设备底层PLC参数数据难以获取的问题，杭州卷烟厂 (朱立明，钱杰，王琪，李钰靓，李勇.一种超高速包装机数据采集系统：中国，CN204056360U [P]. 2014-12-31.)提出了一种基于Beckoff软PLC的超高速包装设备数据采集方法，通过对 设备底层PLC源代码和库文件的分析，能够采集超高速小盒包装机H1000的部分参数数据， 从而为基于多元统计分析方法实施具有多工况特性的超高速小盒包装机H1000状态监测与 诊断提供了丰富的运行数据。
【发明内容】
[0005] 超高速小盒包装机H1000生产过程中，操作人员根据目测的设备运行状态频繁调 整车速设定值，使得生产过程具有多个运行工况。由于小盒包装机H1000的高速运行，从一 个稳定运行工况到另外一个稳定运行工况的时间非常短(如从车速为940包/分钟降低到50 包/分钟仅需20秒，即在两个采样周期内完成），因此本发明仅考虑各个不同的稳定运行工 况。不同稳定运行工况之间数据的均值、方差、相关关系等特征变量会有明显的变化，直接 应用PCA方法进行统计建模和状态监测会产生大量的误报警或漏报警。
[0006] 为满足超高速小盒包装机H1000状态监测与诊断的实际需要，本发明面向设备生 产过程的多工况特性，提出一种包含离线工况识别和划分、稳定工况统计监测建模、多工况 过程状态监测、故障原因变量分离等关键技术的超高速小盒包装机H1000状态监测与诊断 方法。该方法基于滑动时间窗口离线计算稳定度因子识别稳定工况和过渡工况，采用自适 应k-means聚类方法对稳定工况进行划分形成若干个稳定工况数据簇，再利用PCA方法对每 个稳定工况数据簇独立建立统计监测模型。多工况过程状态监测时，根据当前滑动时间窗 口内数据的稳定度因子判断工况类型，若为过渡工况，则将当前时刻T 2、SPE两个监测统计 量赋值为〇;若为稳定工况，计算当前滑动时间窗口内有效数据与各个聚类中心的欧式距 离，利用距离最小的聚类所对应的PCA监测模型进行实时监测，任一统计量超限时利用贡献 图方法确定引起故障的主要过程变量，最终为超高速小盒包装机H1000的状态监测与诊断 提供了一种有效可行的方法。
[0007] 本发明所采用的具体技术方案如下：
[0008] -种具有多工况特性的超高速小盒包装机H1000状态监测与诊断方法，包括：
[0009] 1)根据超高速小盒包装机H1000的运行机理确定统计模型的输入变量，采集某一 机台底层PLC中的设备运行数据，获得二维原始数据MKX J)，K为采样点个数，J为监测变量 个数；
[0010] 2)基于车速检测值，采用滑动时间窗口 Η方法计算稳定度因子SF，当其值较小时 (小于等于阈值ζ)表明过程处于稳定工况，反之表明过程处于过渡工况。采用平滑滤波方法 对处于稳定工况的滑动时间窗口 Η内原始数据进行算术平均，获得一组移动平均数据，将移 动平均数据中"车速"变量值大于阈值Φ的数据作为稳定工况的有效数据D= {Di，D2，…，DN}， N为有效数据的个数；
[0011] 3)稳定工况的有效数据被识别出以后，采用自适应k-means聚类方法对有效数据 的相似度进行聚类分析，获得每一种稳定工况的建模数据.
〕表示将 稳定工况划分成C个聚类，每个聚类中具有相似过程特性的建模数据可以用一个统计监测 模型来描述；
[0012] 4)对每一种稳定工况的建模数据乂
_，C)进行均值中心化和 方差归一化的数据预处理，获得处理后的建模数据Xr(KrmXJ)。对建模数据Xr进行PCA分解， 得到对应的主元个数A r、负载矩阵Pr、控制限和&rSPt,，建立对应第r个稳定工况的PCA 监测模型；
[0013] 5)多工况过程状态监测时，计算当前滑动时间窗口内数据的稳定度因子SFnew，通 过与给定阈值ζ的比较确定当前所处的工况类型。若为过渡工况，则将当前时刻
印 SPEne3W两个监测统计量赋值为0;否则采用平滑滤波方法计算当前滑动时间窗口内的移动平 均数据，并计算其与各个聚类中心的欧式距离，调用距离最小的聚类所对应的PCA监测模型 在线计算和SPEne3W两个监测统计量。如果两个统计量都位于控制限以内，表明设备状态 正常，如果至少其中一个统计量超出控制限，表明设备状态异常；
[0014] 6)当检测到设备出现有异常时，根据超限统计量处于T2的主元子空间或是SPE的 残差子空间，计算每个过程变量对超限统计量的贡献率，贡献率较大的变量被初步确定为 造成设备异常的原因变量。
[0015]所述步骤2)中，采用滑动时间窗口Η方法计算稳定度因子SF，计算公式如下：
limx
[0017] 其中VH max表示滑动时间窗口Η内车速数据的最大值，VH min表示滑动窗口Η内车速 数据的最小值，￥_\表示设备允许的车速最大值。Η表示滑动时间窗口的长度。
[0018] 在步骤3)中，采用自适应k-means聚类方法对有效数据的相似度进行聚类分析，通 过定义两个聚类中心的最小距离的阈值Θ，同时最小化同一聚类数据之间的局部距离平方 和以及所有聚类之间的全局距离平方和，获得最优的聚类数目和聚类中心，实现对稳定工 况的自适应划分。
[0019] 在步骤4)中，对建模数据预处理包括依次进行的减均值、除标准差 处理获得Xr(KrmXJ)，计算公式如下：
[0021] 其中足.^·表示建模数据...... 1第rk行第j列的数据，下标r代表工况、k代
表样本、j代表变量，I.,表示第r个稳定工况第j个变量的均值f表示第r个稳定工况第j 个变量的标准差。
[0022] 将数据Xr进行PCA分解，PCA分解的计算公式为：
[0024] 其中Tr为主元子空间中(KrmXAr)维的得分矩阵;Pr为主元子空间中（JXA r)维的负 载矩阵;Er为残差子空间中（KrmXJ)维的残差矩阵。Ar表示主元个数，由累计贡献率方法确 定。
[0025] 离线计算主成分子空间的Hotelling-T2和残差子空间的SPE(Squared Prediction Error)监测统计量。I]2监测统计量的计算公式为：
[0027] 其中trk = XrkPr为（1 X Ar )维的主元得分向量，Xrk为建模数据Xr的第k行，Pr为主元 子空间的负载矩阵，对角矩阵\ ，· · ·，.也,）由建模数据Xr的协 方差矩阵Σ XrTXr的 前Ar个特征值所构成。
[0028] SPEr监测统计量的计算公式为：
[0030]其中
表示重构得到的（IX J)维估计向量。
[0031 ]根据F分布计算I；2监测统计量的控制限为：
[0033] 其中α为置信度，Fa(Ar，Krm-Ar)是对应于置信度为α，自由度为A r、Krm-Ar条件下的F 分布临界值。
[0034] 根据X2分布计算SPEr监测统计量的控制限:
[0036] 其中g = Vr/(2nr)
nr、vr分别为SPE r监测统计量的均值和方差。
[0037]在步骤5)中，计算处于稳定工况滑动时间窗口的移动平均数据与各个聚类中心的 欧式距离，根据距离最小判断该时刻对应的稳定工况监测模型。欧式距离的计算公式为： r二1
[0039] 其中Xne3W,j表示移动平均数据的第j个变量的检测值，Cem,」表示第i个聚类中心第j 个变量的值，C表示聚类中心的个数。
[0040] 设当前移动平均数据对应第L(1 < L < C)个稳定工况监测模型，利用该监测模型建 模数据的均值I ,和标准差对移动平均数据进行标准化处理，获得An?(lXJ)。将其投 影到第L个监测模型上，计算和SPEL_监测统计量指标。
[0042]其中Pl表示第L个监测模型主元子空间的负载矩阵，对角矩阵SL是由建模数据XL的 协方差矩阵的前Al个特征值所构成：
「表示重构得到的（1 XJ)维估计向 量。
[0043]如果两个监测统计量都位于控制限内，表明超高速小盒包装机处于正常状态，如 果其中至少一个监测统计量超出控制限，表明设备处于异常状态。
[0044] 在步骤6)中，当任一监测统计量超出控制限时，计算各个过程变量对超限统计量 的贡献值，确定引起异常的原因变量。主成分子空间
监测统计量的超出正常控制限 时，主成分得分变量1：1^对.
的贡献率可以计算如下：
[0046]其中tLa表示对应第L个监测模型的第a个主成分得分变量，表示对应第L个监测 模型的第a个特征值。
[0047]过程变量ane3W, j对tu的贡献率可以计算如下：
[0049]其中pLa,j表示对应第L个监测模型的第a个负载向量的第j个变量。
[0050] 当残差子空间SPELnew监测统计量的超出正常控制限时，过程变量xLnew,j对SPE Lnew的 贡献率可以计算如下：
[0052]其彳
5重构估计向1
1，的第j个变量
示残 差的正负信息。
[0053]本发明的有益效果是：
[0054] 本发明主要涉及离线工况识别和划分、稳定工况统计监测建模、多工况过程状态 监测、故障原因变量分离等关键技术。基于海量历史数据，实现了多工况过程稳定工况、过 渡工况的有效识别以及稳定工况的自适应划分，进而自动完成稳定工况的离线监测建模及 控制限的求取。多工况过程状态监测时，匹配当前稳定工况的移动平均数据与各个聚类中 心的相似度判断该时刻对应的稳定工况监测模型，对当前的生产过程进行有效监测，监测 异常发生时采用变量贡献图追溯和确定引起故障的原因变量。最终形成一套集工况识别、 统计建模、状态监测和故障诊断于一体的切实可行的实施方案，能有效提高设备诊断结果 的可靠性和准确性，从而为超高速小盒包装机H1000的保养和维护提供科学决策和有效指 导。
【附图说明】
[0055] 图1为本发明的具有多工况特性超高速小盒包装机H1000状态监测与诊断流程图；
[0056] 图2为工况识别前1~1000采样点的H1000车速曲线图；
[0057] 图3为1~1000采样点H1000车速的稳定度因子SF曲线图；
[0058] 图4为经过平滑滤波后稳定工况的H1000车速曲线图；
[0059] 图5为稳定工况有效数据的H1000车速曲线图；
[0060] 图6为多工况过程正常状态测试数据的T2和SPE过程监测图；
[0061]图7为多工况过程故障状态数据1的T2和SPE过程监测图；
[0062] 图8为多工况过程故障状态数据2的T2和SPE过程监测图；
[0063] 图9为多工况过程故障状态数据3的T2和SPE过程监测图；
[0064] 图10为故障状态1下首次故障报警时刻的T2变量贡献图；
[0065] 图11为故障状态2下首次故障报警时刻的T2变量贡献图；
[0066] 图12为故障状态3下首次故障报警时刻的T2变量贡献图。
【具体实施方式】
[0067]为了更好的理解本发明的技术方案，以下结合说明书附图对本发明的实施方式作 进一步描述。本实施是用于意大利G.D公司超高速小盒包装机H1000的状态监测与诊断方 法。G. D公司超高速包装设备由小盒包装机H1000、小盒薄膜封装机W1000和条盒包装机BV三 部分组成。小盒包装机H1000是一台高速连续运动硬盒翻盖包装机，其运行以间歇式往复运 动和公转加自传的行星驱动为主，利用烟支、内衬纸、内框纸、商标纸将烟支按7-6-7共20支 的形式进行翻盖式小盒包装。本发明的具有多工况特性的超高速小盒包装机状态监测与诊 断方法的实现框图如图1所示，方法主要分为以下几步：
[0068] (1)历史运行数据的获取
[0069] 根据超高速小盒包装机H1000的运行机理确定统计模型的输入变量，采集某一机 台底层PLC中的设备运行数据，获得二维原始数据p(K X J)，K为采样点个数，J为监测变量个 数。
[0070] 本实例中，选取Α4机台的超高速小盒包装机，采集的设备运行数据中包含以下变 量:伺服电机电流、伺服电机温度、热熔胶温度、车速等44个监测变量，如表1所示。监测变量 的采样频率为10秒/次，获得的原始数据位64394 X 44)。
[0071] 表1超高速小盒包装机的监测变量表
[0074] (2)离线工况的识别
[0075]基于车速检测值，采用滑动时间窗口Η方法计算稳定度因子SF，当其值较小时（小 于等于阈值ζ)表明过程处于稳定工况，反之表明过程处于过渡工况。采用平滑滤波方法对 处于稳定工况的滑动时间窗口 Η内原始数据进行算术平均，获得一组移动平均数据，将移动 平均数据中"车速"变量值大于阈值Φ的数据作为稳定工况的有效数据D= {D^D^···，DN}，Ν 为有效数据的个数。
[0076]采用滑动时间窗口Η方法计算稳定度因子SF，计算公式如下：
[0078] 其中VH max表示滑动时间窗口Η内车速数据的最大值，VH min表示滑动窗口Η内车速 数据的最小值，￥_\表示设备允许的车速最大值。Η表示滑动时间窗口的长度。
[0079]本实例中，滑动时间窗口的长度Η取为6，允许的车速最大值Vmax取为1000，阈值ζ取 为0.02,阈值Φ取为800,获得的有效数据为D(30283X 44)。工况识别前1~1000采样点的 H1000车速如图2所示，相应的稳定度因子如图3所示，经过平滑滤波后稳定工况的H1000车 速如图4所示，稳定工况有效数据的H1000车速如图5所示。
[0080] (3)稳定工况的划分
[0081]稳定工况的有效数据被识别出以后，采用自适应k-means聚类方法对有效数据的 相似度进行聚类分析，获得每一种稳定工况的建模数据
〕表示将稳定 工况划分成C个聚类，每个聚类中具有相似过程特性的建模数据可以用一个统计模型来描 述。自适应k-means聚类方法的求取过程如下：
[0082] 步骤1:设定初始聚类中心个数Co,两个聚类中心的最小距离阈值Θ，每个聚类俘获 的最小数据量阈值minL。初始化Co个聚类中心Wi,。(i = 1，…，J; c = 1，…，Co);
[0083] 步骤2:计算两个聚类中心的欧式距离虹#(11,。1，1 1,。2)，若其小于预定的阈值0，剔 除其中一个聚类中心；
[0084] 步骤3:计算每个稳定工况的有效数据01{仏=1，2，一，《到所有聚类中心的欧式距 离disi^Di^Wi,。），若Dk与第c*类的中心的距离最小，则将Dk的隶 属关系定义为m(k) = c*;
[0085] 步骤4:Inum次迭代后，若某聚类中心俘获的数据样本个数小于预定的阈值minL，则 剔除该类。
[0086] 步骤5:更新聚类个数为C1+1，并根据数据的隶属关系重新计算新的聚类中心W1+1,。 (c = l ,2, ··· ,Ci+i) 〇
[0087] 步骤6:如果算法满足收敛条件则结束，否则返回步骤2,进行下一次迭代计算。收 敛条件有:两次迭代中的聚类中心距离的变化小于一个很小的阈值ε，或者每个聚类中数据 Dk到聚类中心的距离平方和以及聚类之间的距离平方和达到最小。
[0088]本实例中，初始聚类中心个数Co取为15，最小距离阈值Θ取为300，最小数据量阈值minL 取为100,阈值ε取为0.0001。通过自适应聚类，稳定工况有效数据被划分为7个聚类(即C=7)，获得每 一种稳定工况的建模数据
[0089] (4)稳定工况统计监测建模
[0090]对每一种稳定工况的建模数#
进行均值中心化和方 差归一化的数据预处理，获得处理后的建模数据Xr(KrmXJ)。对建模数据Xr进行PCA分解，得 到对应的主元个数Ar、负载矩阵Pr、控制限和，建立对应第r个稳定工况的PCA监 测模型。
[0091]对建模数据预处理包括依次进行的减均值、除标准差处理获得Xr (KrmX J)，计算公式如下： ?
[0093] 其中足υ表示建模数?： 第rk行第j列的数据，下标r代表工况数、k 代表样本数、j代表变量数，.....-况第j个变量的均值表示第r个稳定 工况第j个变量的标准差。
[0094] 将数据Xr进行PCA分解，PCA分解的计算公式为：
[0096] 其中Tr为主元子空间中（KrmXAr)维的得分矩阵;Pr为主元子空间中（JXAr)维的负 载矩阵;Er为残差子空间中（K rmXJ)维的残差矩阵。Ar表示主元个数，由累计贡献率方法确 定。
[0097] 离线计算主成分子空间的Hotelling-T2和残差子空间的SPE(Squared Pred i c t i on Error)监测统计量。监测统计量的计算公式为：
[0099] 其中trk = XrkPr为（1 X Ar )维的主元得分向量，Xrk为建模数据Xr的第k行，Pr为主元 子空间的负载矩阵，对角矩阵& = …，\ )由建模数据Xr的协方差矩阵Σ XrTXr的 前Ar个特征值所构成。
[0100] SPEr监测统计量的计算公式为：
[0102] 其中
表示重构得到的（IX J)维估计向量。
[0103] 根据F分布计算f监测统计量的控制限#%为：
[0105] 其中α为置信度，Fa(Ar，Krm-Ar)是对应于置信度为α，自由度为A r、K?-Ar条件下的F 分布临界值。
[0106] 根据X2分布计算SPEr监测统计量的控制阳
[0108] 其：
，_，nr、vr分别为SPE r监测统计量的均值和方差。
[0109] 本实例中，7个PCA监测模型的主元个数依次为A! = 26、A2 = 21、A3 = 26、A4 = 26、A5 = 27、六6 = 26^7 = 26，则?1为(44\26)负载矩阵、？2为(44\21)负载矩阵、？3为(44\26)负载矩 阵、P 4为(44X26)负载矩阵、Ρ^(44 X 27)负载矩阵、P6为(44X26)负载矩阵、P7为(44X26) 负载矩阵。统计控制限的置信度α取为0.99,监测模型1的控制限和 c?rSPE( =.16:·3218、监测模型 2 的控制限 cir# =40.88:73:和 cirSPt: =8,5462、监测模型 3 的 控制限改45.9832和你啊=14·9丨56、监测模型4的控制限g%. = 45.9867和 cirSPE4 =23.0371、监测模型5的控制限冰> =47.2113和=51.4481.、监测模型6 的控制限= 46.0821和C%PE:6 = 12.6:8:84、监测模型7的控制限= 45.9559和 cthm7 = I B.5386
[0110] (5)多工况过程状态监测
[0111] 多工况过程状态监测时，计算当前滑动时间窗口内数据的稳定度因子SFnew，通过 与给定阈值ζ的比较确定当前所处的工况类型。若为过渡工况，则将当前时刻7；i和SPEnew 两个监测统计量赋值为0;否则采用平滑滤波方法计算当前滑动时间窗口内的移动平均数 据，并计算其与各个聚类中心的欧式距离，调用距离最小的聚类所对应的PCA监测模型在线 雅4邮PE-wW傾灘i十1：。娜W侦i十胃雛刊空·艮眺綱#?紅 常，如果至少其中一个统计量超出控制限，表明设备状态异常。
[0112] 计算处于稳定工况滑动时间窗口的移动平均数据与各个聚类中心的欧式距离，根 据距离最小判断该时刻对应的稳定工况监测模型。欧式距离的计算公式为：
[0114] 其中Xne3W,j表示移动平均数据的第j个变量的检测值，Cem.j表示第i个聚类中心第j 个变量的值，C表示聚类中心的个数。
[0115] 设当前移动平均数据对应第L(1<L<C)个稳定工况监测模型，利用该监测模型建 模数据的均值忑，_/和标准差毛，_/对移动平均数据进行标准化处理，获得xl_(1 X J)。将其投 影到第L个监测模型上，计算和SPEL_监测统计量指标。
[0117]其中Pl表示第L个监测模型主元子空间的负载矩阵，对角矩阵SL是由建模数据XL的 协方差矩阵的前Al个特征值所构成
表示重构得到的（1 XJ)维估计向 量。
[0118] 如果两个监测统计量都位于控制限内，表明超高速小盒包装机处于正常状态，如 果其中至少一个监测统计量超出控制限，表明设备处于异常状态。
[0119] 本实例中，选取多工况过程8816组采样数据％ew(8816X44)进行分析，该测试数据 T2和SPE统计量的监测结果如图6所示。引入首次故障i警时间，定义为连续9个监测统计量 超过对应控制限的采样时间点。在两个监测子空间中有个别样本点的监测统计量超出控制 限，根据首次故障报警的定义，两个监测统计量均未发生故障报警，对当前过程的正常状态 做出了正确的指示，表明所提状态监测与诊断方法具有准确监测多工况过程正常状态的能 力。
[0120]作为比较，重新选取多工况过程4505组采样数据％ew(4505X44)进行状态监测与 诊断的分析。在第1001个样本点的时候，分别令热熔胶温i阶跃增长1% (突变故障）、A1电 机电流阶跃降低30% (突变故障）、A5A电机温度以0.03的斜率增加(渐变故障），一直持续到 第2000个样本点的时候恢复正常，形成故障状态数据1、2和3。故障状态数据T 2和SPE统计量 的监测结果分别如图7至图9所示。
[0121] 对于故障1，主元子空间的T2统计量出现连续超限，主成分子空间出现异常，检测 出的首次故障报警时间为第1030个样本点时刻，即监测模型延后290秒准确检测出故障的 发生;对于故障2,主元子空间的T 2统计量出现连续超限，主成分子空间出现异常，检测出的 首次故障报警时间为第1030个样本点时刻，即监测模型延后290秒准确检测出故障的发生； 对于故障3,主元子空间的T 2统计量出现连续超限，主成分子空间出现异常，检测出的首次 故障报警时间为第1318个样本点时刻，即监测模型延后3170秒准确检测出故障的发生。说 明所提方法能对多工况过程的突变和渐变故障做出正确指示，表明该方法具有准确监测多 工况过程故障状态的能力。
[0122] (6)基于贡献图的故障诊断
[0123] 当任一监测统计量超出控制限时，计算各个过程变量对超限统计量的贡献值，确 定引起异常的原因变量。主成分子空间，监测统计量的超出正常控制限时，主成分得分 变量tu
贡献率可以计算如下：
[0125] 其中tu表示对应第L个监测模型的第a个主成分得分变量，表示对应第L个监测 模型的第a个特征值。
[0126] 过程变量对tu的贡献率可以计算如下：
[0128]其中plm表示对应第L个监测模型的第a个负载向量的第j个变量 。
[0129] 当残差子空间SPELnew监测统计量的超出正常控制限时，过程变量xLnew,j对SPE Lnew的 贡献率可以计算如下：
L,n(￡W
[0131] 其中
.示重构估计向I
:勺第j个变量，sign
I示残 差的正负信息。
[0132] 本实例中，根据图7的监测结果，故障发生在主成分子空间中。在首次故障报警的 第1030个样本点时刻，过程变量对7^0，监测统计量超限的贡献率如图10所示，确定引起故 障的原因变量为第43个变量，即H1000热熔胶温度。
[0133] 根据图8的监测结果，故障发生在主成分子空间中。在首次故障报警的第1030个样 本点时刻，过程变量对^^胃监测统计量超限的贡献率如图11所示，确定引起故障的原因变 量为第1个变量，即A1电机电流。
[0134] 根据图9的监测结果，故障发生在主成分子空间中。在首次故障报警的第1318个样 本点时刻，过程变量对监测统计量超限的贡献率如图12所示，确定引起故障的原因变 量为第16个变量，即A5A电机温度。
[0135] 以上所述仅为本发明的较佳实施举例，并不用于限制本发明，凡在本发明精神和 原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
【主权项】
1. 一种具有多工况特性的超高速小盒包装机状态监测与诊断方法，其特征在于，包括 步骤： 1) 根据超高速小盒包装机HlOOO的运行机理确定统计模型的输入变量，采集相应的运 行数据，获得二维原始数据旮(Kx·/)，K为采样点个数，J为监测变量个数； 2) 基于车速检测值，采用滑动时间窗口 H方法计算稳定度因子SF，根据稳定度因子SF与 阈值ζ的关系判断是否处于稳定工况，采用平滑滤波方法对处于稳定工况的滑动时间窗口 H 内原始数据进行算术平均，获得一组移动平均数据，将移动平均数据中"车速"变量值大于 阈值Φ的数据作为稳定工况的有效数据D= (D1 ,D2，…，DN}，Ν为有效数据的个数； 3) 采用自适应k-means聚类方法对有效数据的相似度进行聚类分析，获得每一种稳定 工况的建模数据C表示将稳定工况划分成C个聚类，每个聚类中具 有相似过程特性的建模数据用一统计监测模型来描述； 4) 对每一种稳定工况的建模数据，进行均值中心化和方差 归一化的数据预处理，获得处理后的建模数据Xr (Krm X J )，对建模数据Xr进行PCA分解，得到 对应的主元个数Ar、负载矩阵Pr、控制限和&胃，建立对应第 r个稳定工况的PCA监测 丰旲型； 5) 实时计算当前滑动时间窗口内数据的稳定度因子SFne3w,通过与给定阈值ζ的比较确 定当前所处的工况类型;若为过渡工况，则将当前时刻:ζ?,和SPE new两个监测统计量赋值为 〇;否则采用平滑滤波方法计算当前滑动时间窗口内的移动平均数据，并计算其与各个聚类 中心的欧式距离，调用距离最小的聚类所对应的PCA监测模型在线计算T^jPSPE new两个监 测统计量，若一个监测统计量超出控制限，表明设备状态异常； 6) 当检测到设备出现有异常时，根据超限的监测统计量处于T2的主元子空间或是SPE的 残差子空间，计算每个过程变量对超限统计量的贡献率，贡献率较大的变量被初步确定为 造成设备异常的原因变量。2. 如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，在步骤2)中采用滑动时间窗口 H方法计算稳定度因子SF进行工况的识别，计算公 式如下：其中Vh max表示滑动时间窗口H内车速数据的最大值，Vh min表示滑动窗口H内车速数据的 最小值，Vmax表示设备允许的车速最大值，H表示滑动时间窗口的长度。3. 如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，在步骤3)中采用自适应k-means聚类方法进行工况的划分，对有效数据的相似度 进行聚类分析，通过定义两个聚类中心的最小距离的阈值Θ，同时最小化同一聚类数据之间 的局部距离平方和以及所有聚类之间的全局距离平方和，获得最优的聚类数目和聚类中 心，实现对稳定工况的自适应划分。4. 如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，在步骤4)中对建模数据预处理包括依次进行的减均值、除标准差 处理获得以^】），计算公式如下：其中，表示建模数据中第rk行第j列的数据，下标r代表工况、k代表样 本、j代表变量，巧,7:表示第r个稳定工况第j个变量的均值表示第r个稳定工况第j个变 量的标准差。5. 如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，所述步骤4)中PCA分解的计算公式为：其中，Tr为主元子空间中（KrmXAr)维的得分矩阵;Pr为主元子空间中（JXAr)维的负载 矩阵;Er为残差子空间中(KrmX J)维的残差矩阵;Ar表示主元个数，由累计贡献率方法确定； 离线计算主成分子空间的HotelIing-T2和残差子空间的SPE(Squared Prediction Error)监测统计量；I；2监测统计量的计算公式为：其中trk = XrkPr为（I X Ar )维的主元得分向量，Xrk为建模数据Xr的第k行，Pr为主元子空间 的负载矩阵，对角矩阵由建模数据Xr的协方差矩阵的前Ar 个特征值所构成； SPEr监测统计量的计算公式为：其中表示重构得到的（I XJ)维估计向量； 根据F分布计算?；2监测统计量的控制限为：其中α为置信度，Fa(Ar，Krm-Ar)是对应于置信度为α，自由度为A r、Krm-Ar条件下的F分布 临界值； 根据X2分布计算SPEr监测统计量的控制限#rspi为··；/?,、Vi.分别为SPEr监测统计量的均值和方差。6. 如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，在步骤5)中，欧式距离的计算公式为：其中xn?,j表示移动平均数据的第j个变量的检测值，Cem, j表示第i个聚类中心第j个变 量的值，C表示聚类中心的个数。7. 如权利要求6所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，设当前移动平均数据对应第L(I^LSC)个稳定工况监测模型，利用该监测模型 建模数据的均值&,,和标准差^^对移动平均数据进行标准化处理，获得an?(l X J)，将其 投影到第L个监测模型上，计算J^jPSPEL_监测统计量指标；其中Pl表示第L个监测模型主元子空间的负载矩阵，对角矩阵Sl是由建模数据Xl的协方 差矩阵的前Ai个特征值所构成，表示重构得到的（I XJ)维估计向量。8. 如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，在步骤6)中，主成分子空间?^_监测统计量的超出正常控制限时，主成分得分变 散3对jS-颜鮮卿雅女吓：其中Ua表示对应第L个监测模型的第a个主成分得分变量，表示对应第L个监测模型 的第a个特征值； 过程变量XkW, j对tu的贡献率可以计算如下：其中PUd表示对应第L个监测模型的第a个负载向量的第j个变量； 当残差子空间SPELnew监测统计量的超出正常控制限时，过程变量xLnew,j对SPE Lnew的贡献 率可以计算如下：其中表示重构估计向量.的第j个变量，1表示残差的正 负信息。
【专利摘要】本发明公开了一种具有多工况特性的超高速小盒包装机状态监测与诊断方法，该方法针对运行过程多工况特性及工况普遍未知的问题，基于滑动窗口计算稳定度因子识别不同工况，采用自适应聚类方法对稳定工况进行划分形成若干稳定工况数据簇，利用PCA方法为每个稳定工况数据簇建立监测模型。多工况过程状态监测时，根据当前滑动窗口内数据的稳定度因子判断工况类型，若为过渡工况，则将监测统计量赋值为0；若为稳定工况，计算当前有效数据与各个聚类中心的距离，获得匹配的PCA监测模型进行实时监测，任一统计量超限时采用贡献图确定引起故障的主要过程变量。
【IPC分类】G05B23/02
【公开号】CN105487524
【申请号】CN201511007689
【发明人】王伟, 赵春晖, 楼卫东, 李钰靓, 张利宏, 熊月宏
【申请人】浙江中烟工业有限责任公司
【公开日】2016年4月13日
【申请日】2015年12月29日
【技术领域】
[0001] 本发明涉及意大利G.D公司超高速小盒包装机H1000的统计监测建模、状态监测与 诊断技术。
【背景技术】
[0002] 随着近年来烟草行业全球化、新的法律法规以及外部环境的重大变化，中国烟草 面临愈来愈严峻的竞争压力和社会压力。在保证产品质量的前提下，有效提升设备的智能 化水平和高效运行能力成为卷烟工厂关注的焦点。智能化，即生产设备具有感知、分析、推 理和决策等功能；高效运行，即在规定时间内消耗更少的资源生产更多、更好的产品。卷烟 工厂设备的智能化水平和高效运行能力与其自动化程度密不可分，设备状态监测与诊断技 术是提升智能化水平、保证高效运行的重要手段。
[0003] 意大利G.D公司超高速包装设备是中国烟草总公司"十二五"期间重点引进的目前 世界上最高速、最先进的卷烟包装设备，其设计生产能力为1〇〇〇包香烟/分钟。超高速包装 设备由小盒包装机H1000、小盒薄膜封装机W1000和条盒包装机BV三部分组成。小盒包装机 H1000是一台高速连续运动硬盒翻盖包装机，利用烟支、内衬纸、内框纸、商标纸将烟支按7-6-7共20支的形式进行翻盖式小盒包装。小盒包装机H1000是超高速包装设备最关键、最核 心的单元，也是生产现场故障最集中的部位。由于G.D公司超高速包装设备底层PLC开放的 程序源代码和库文件非常有限，生产现场主要通过每班一组的数据报表、生产时的人工目 测、轮保时的人工检测、设备故障的计划维修和事后抢修等方式保证小盒包装机H1000的设 备完好性。杭州卷烟厂2014年超高速包装设备的运行效率仅为70.08%，目前采用的设备状 态监测方式已经不能满足卷烟生产的实际需要，亟需一种切实有效的超高速小盒包装机 H1000状态监测与诊断方法。
[0004] 为解决G.D公司超高速包装设备底层PLC参数数据难以获取的问题，杭州卷烟厂 (朱立明，钱杰，王琪，李钰靓，李勇.一种超高速包装机数据采集系统：中国，CN204056360U [P]. 2014-12-31.)提出了一种基于Beckoff软PLC的超高速包装设备数据采集方法，通过对 设备底层PLC源代码和库文件的分析，能够采集超高速小盒包装机H1000的部分参数数据， 从而为基于多元统计分析方法实施具有多工况特性的超高速小盒包装机H1000状态监测与 诊断提供了丰富的运行数据。
【发明内容】
[0005] 超高速小盒包装机H1000生产过程中，操作人员根据目测的设备运行状态频繁调 整车速设定值，使得生产过程具有多个运行工况。由于小盒包装机H1000的高速运行，从一 个稳定运行工况到另外一个稳定运行工况的时间非常短(如从车速为940包/分钟降低到50 包/分钟仅需20秒，即在两个采样周期内完成），因此本发明仅考虑各个不同的稳定运行工 况。不同稳定运行工况之间数据的均值、方差、相关关系等特征变量会有明显的变化，直接 应用PCA方法进行统计建模和状态监测会产生大量的误报警或漏报警。
[0006] 为满足超高速小盒包装机H1000状态监测与诊断的实际需要，本发明面向设备生 产过程的多工况特性，提出一种包含离线工况识别和划分、稳定工况统计监测建模、多工况 过程状态监测、故障原因变量分离等关键技术的超高速小盒包装机H1000状态监测与诊断 方法。该方法基于滑动时间窗口离线计算稳定度因子识别稳定工况和过渡工况，采用自适 应k-means聚类方法对稳定工况进行划分形成若干个稳定工况数据簇，再利用PCA方法对每 个稳定工况数据簇独立建立统计监测模型。多工况过程状态监测时，根据当前滑动时间窗 口内数据的稳定度因子判断工况类型，若为过渡工况，则将当前时刻T 2、SPE两个监测统计 量赋值为〇;若为稳定工况，计算当前滑动时间窗口内有效数据与各个聚类中心的欧式距 离，利用距离最小的聚类所对应的PCA监测模型进行实时监测，任一统计量超限时利用贡献 图方法确定引起故障的主要过程变量，最终为超高速小盒包装机H1000的状态监测与诊断 提供了一种有效可行的方法。
[0007] 本发明所采用的具体技术方案如下：
[0008] -种具有多工况特性的超高速小盒包装机H1000状态监测与诊断方法，包括：
[0009] 1)根据超高速小盒包装机H1000的运行机理确定统计模型的输入变量，采集某一 机台底层PLC中的设备运行数据，获得二维原始数据MKX J)，K为采样点个数，J为监测变量 个数；
[0010] 2)基于车速检测值，采用滑动时间窗口 Η方法计算稳定度因子SF，当其值较小时 (小于等于阈值ζ)表明过程处于稳定工况，反之表明过程处于过渡工况。采用平滑滤波方法 对处于稳定工况的滑动时间窗口 Η内原始数据进行算术平均，获得一组移动平均数据，将移 动平均数据中"车速"变量值大于阈值Φ的数据作为稳定工况的有效数据D= {Di，D2，…，DN}， N为有效数据的个数；
[0011] 3)稳定工况的有效数据被识别出以后，采用自适应k-means聚类方法对有效数据 的相似度进行聚类分析，获得每一种稳定工况的建模数据.
〕表示将 稳定工况划分成C个聚类，每个聚类中具有相似过程特性的建模数据可以用一个统计监测 模型来描述；
[0012] 4)对每一种稳定工况的建模数据乂
_，C)进行均值中心化和 方差归一化的数据预处理，获得处理后的建模数据Xr(KrmXJ)。对建模数据Xr进行PCA分解， 得到对应的主元个数A r、负载矩阵Pr、控制限和&rSPt,，建立对应第r个稳定工况的PCA 监测模型；
[0013] 5)多工况过程状态监测时，计算当前滑动时间窗口内数据的稳定度因子SFnew，通 过与给定阈值ζ的比较确定当前所处的工况类型。若为过渡工况，则将当前时刻
印 SPEne3W两个监测统计量赋值为0;否则采用平滑滤波方法计算当前滑动时间窗口内的移动平 均数据，并计算其与各个聚类中心的欧式距离，调用距离最小的聚类所对应的PCA监测模型 在线计算和SPEne3W两个监测统计量。如果两个统计量都位于控制限以内，表明设备状态 正常，如果至少其中一个统计量超出控制限，表明设备状态异常；
[0014] 6)当检测到设备出现有异常时，根据超限统计量处于T2的主元子空间或是SPE的 残差子空间，计算每个过程变量对超限统计量的贡献率，贡献率较大的变量被初步确定为 造成设备异常的原因变量。
[0015]所述步骤2)中，采用滑动时间窗口Η方法计算稳定度因子SF，计算公式如下：
limx
[0017] 其中VH max表示滑动时间窗口Η内车速数据的最大值，VH min表示滑动窗口Η内车速 数据的最小值，￥_\表示设备允许的车速最大值。Η表示滑动时间窗口的长度。
[0018] 在步骤3)中，采用自适应k-means聚类方法对有效数据的相似度进行聚类分析，通 过定义两个聚类中心的最小距离的阈值Θ，同时最小化同一聚类数据之间的局部距离平方 和以及所有聚类之间的全局距离平方和，获得最优的聚类数目和聚类中心，实现对稳定工 况的自适应划分。
[0019] 在步骤4)中，对建模数据预处理包括依次进行的减均值、除标准差 处理获得Xr(KrmXJ)，计算公式如下：
[0021] 其中足.^·表示建模数据...... 1第rk行第j列的数据，下标r代表工况、k代
表样本、j代表变量，I.,表示第r个稳定工况第j个变量的均值f表示第r个稳定工况第j 个变量的标准差。
[0022] 将数据Xr进行PCA分解，PCA分解的计算公式为：
[0024] 其中Tr为主元子空间中(KrmXAr)维的得分矩阵;Pr为主元子空间中（JXA r)维的负 载矩阵;Er为残差子空间中（KrmXJ)维的残差矩阵。Ar表示主元个数，由累计贡献率方法确 定。
[0025] 离线计算主成分子空间的Hotelling-T2和残差子空间的SPE(Squared Prediction Error)监测统计量。I]2监测统计量的计算公式为：
[0027] 其中trk = XrkPr为（1 X Ar )维的主元得分向量，Xrk为建模数据Xr的第k行，Pr为主元 子空间的负载矩阵，对角矩阵\ ，· · ·，.也,）由建模数据Xr的协 方差矩阵Σ XrTXr的 前Ar个特征值所构成。
[0028] SPEr监测统计量的计算公式为：
[0030]其中
表示重构得到的（IX J)维估计向量。
[0031 ]根据F分布计算I；2监测统计量的控制限为：
[0033] 其中α为置信度，Fa(Ar，Krm-Ar)是对应于置信度为α，自由度为A r、Krm-Ar条件下的F 分布临界值。
[0034] 根据X2分布计算SPEr监测统计量的控制限:
[0036] 其中g = Vr/(2nr)
nr、vr分别为SPE r监测统计量的均值和方差。
[0037]在步骤5)中，计算处于稳定工况滑动时间窗口的移动平均数据与各个聚类中心的 欧式距离，根据距离最小判断该时刻对应的稳定工况监测模型。欧式距离的计算公式为： r二1
[0039] 其中Xne3W,j表示移动平均数据的第j个变量的检测值，Cem,」表示第i个聚类中心第j 个变量的值，C表示聚类中心的个数。
[0040] 设当前移动平均数据对应第L(1 < L < C)个稳定工况监测模型，利用该监测模型建 模数据的均值I ,和标准差对移动平均数据进行标准化处理，获得An?(lXJ)。将其投 影到第L个监测模型上，计算和SPEL_监测统计量指标。
[0042]其中Pl表示第L个监测模型主元子空间的负载矩阵，对角矩阵SL是由建模数据XL的 协方差矩阵的前Al个特征值所构成：
「表示重构得到的（1 XJ)维估计向 量。
[0043]如果两个监测统计量都位于控制限内，表明超高速小盒包装机处于正常状态，如 果其中至少一个监测统计量超出控制限，表明设备处于异常状态。
[0044] 在步骤6)中，当任一监测统计量超出控制限时，计算各个过程变量对超限统计量 的贡献值，确定引起异常的原因变量。主成分子空间
监测统计量的超出正常控制限 时，主成分得分变量1：1^对.
的贡献率可以计算如下：
[0046]其中tLa表示对应第L个监测模型的第a个主成分得分变量，表示对应第L个监测 模型的第a个特征值。
[0047]过程变量ane3W, j对tu的贡献率可以计算如下：
[0049]其中pLa,j表示对应第L个监测模型的第a个负载向量的第j个变量。
[0050] 当残差子空间SPELnew监测统计量的超出正常控制限时，过程变量xLnew,j对SPE Lnew的 贡献率可以计算如下：
[0052]其彳
5重构估计向1
1，的第j个变量
示残 差的正负信息。
[0053]本发明的有益效果是：
[0054] 本发明主要涉及离线工况识别和划分、稳定工况统计监测建模、多工况过程状态 监测、故障原因变量分离等关键技术。基于海量历史数据，实现了多工况过程稳定工况、过 渡工况的有效识别以及稳定工况的自适应划分，进而自动完成稳定工况的离线监测建模及 控制限的求取。多工况过程状态监测时，匹配当前稳定工况的移动平均数据与各个聚类中 心的相似度判断该时刻对应的稳定工况监测模型，对当前的生产过程进行有效监测，监测 异常发生时采用变量贡献图追溯和确定引起故障的原因变量。最终形成一套集工况识别、 统计建模、状态监测和故障诊断于一体的切实可行的实施方案，能有效提高设备诊断结果 的可靠性和准确性，从而为超高速小盒包装机H1000的保养和维护提供科学决策和有效指 导。
【附图说明】
[0055] 图1为本发明的具有多工况特性超高速小盒包装机H1000状态监测与诊断流程图；
[0056] 图2为工况识别前1~1000采样点的H1000车速曲线图；
[0057] 图3为1~1000采样点H1000车速的稳定度因子SF曲线图；
[0058] 图4为经过平滑滤波后稳定工况的H1000车速曲线图；
[0059] 图5为稳定工况有效数据的H1000车速曲线图；
[0060] 图6为多工况过程正常状态测试数据的T2和SPE过程监测图；
[0061]图7为多工况过程故障状态数据1的T2和SPE过程监测图；
[0062] 图8为多工况过程故障状态数据2的T2和SPE过程监测图；
[0063] 图9为多工况过程故障状态数据3的T2和SPE过程监测图；
[0064] 图10为故障状态1下首次故障报警时刻的T2变量贡献图；
[0065] 图11为故障状态2下首次故障报警时刻的T2变量贡献图；
[0066] 图12为故障状态3下首次故障报警时刻的T2变量贡献图。
【具体实施方式】
[0067]为了更好的理解本发明的技术方案，以下结合说明书附图对本发明的实施方式作 进一步描述。本实施是用于意大利G.D公司超高速小盒包装机H1000的状态监测与诊断方 法。G. D公司超高速包装设备由小盒包装机H1000、小盒薄膜封装机W1000和条盒包装机BV三 部分组成。小盒包装机H1000是一台高速连续运动硬盒翻盖包装机，其运行以间歇式往复运 动和公转加自传的行星驱动为主，利用烟支、内衬纸、内框纸、商标纸将烟支按7-6-7共20支 的形式进行翻盖式小盒包装。本发明的具有多工况特性的超高速小盒包装机状态监测与诊 断方法的实现框图如图1所示，方法主要分为以下几步：
[0068] (1)历史运行数据的获取
[0069] 根据超高速小盒包装机H1000的运行机理确定统计模型的输入变量，采集某一机 台底层PLC中的设备运行数据，获得二维原始数据p(K X J)，K为采样点个数，J为监测变量个 数。
[0070] 本实例中，选取Α4机台的超高速小盒包装机，采集的设备运行数据中包含以下变 量:伺服电机电流、伺服电机温度、热熔胶温度、车速等44个监测变量，如表1所示。监测变量 的采样频率为10秒/次，获得的原始数据位64394 X 44)。
[0071] 表1超高速小盒包装机的监测变量表
[0074] (2)离线工况的识别
[0075]基于车速检测值，采用滑动时间窗口Η方法计算稳定度因子SF，当其值较小时（小 于等于阈值ζ)表明过程处于稳定工况，反之表明过程处于过渡工况。采用平滑滤波方法对 处于稳定工况的滑动时间窗口 Η内原始数据进行算术平均，获得一组移动平均数据，将移动 平均数据中"车速"变量值大于阈值Φ的数据作为稳定工况的有效数据D= {D^D^···，DN}，Ν 为有效数据的个数。
[0076]采用滑动时间窗口Η方法计算稳定度因子SF，计算公式如下：
[0078] 其中VH max表示滑动时间窗口Η内车速数据的最大值，VH min表示滑动窗口Η内车速 数据的最小值，￥_\表示设备允许的车速最大值。Η表示滑动时间窗口的长度。
[0079]本实例中，滑动时间窗口的长度Η取为6，允许的车速最大值Vmax取为1000，阈值ζ取 为0.02,阈值Φ取为800,获得的有效数据为D(30283X 44)。工况识别前1~1000采样点的 H1000车速如图2所示，相应的稳定度因子如图3所示，经过平滑滤波后稳定工况的H1000车 速如图4所示，稳定工况有效数据的H1000车速如图5所示。
[0080] (3)稳定工况的划分
[0081]稳定工况的有效数据被识别出以后，采用自适应k-means聚类方法对有效数据的 相似度进行聚类分析，获得每一种稳定工况的建模数据
〕表示将稳定 工况划分成C个聚类，每个聚类中具有相似过程特性的建模数据可以用一个统计模型来描 述。自适应k-means聚类方法的求取过程如下：
[0082] 步骤1:设定初始聚类中心个数Co,两个聚类中心的最小距离阈值Θ，每个聚类俘获 的最小数据量阈值minL。初始化Co个聚类中心Wi,。(i = 1，…，J; c = 1，…，Co);
[0083] 步骤2:计算两个聚类中心的欧式距离虹#(11,。1，1 1,。2)，若其小于预定的阈值0，剔 除其中一个聚类中心；
[0084] 步骤3:计算每个稳定工况的有效数据01{仏=1，2，一，《到所有聚类中心的欧式距 离disi^Di^Wi,。），若Dk与第c*类的中心的距离最小，则将Dk的隶 属关系定义为m(k) = c*;
[0085] 步骤4:Inum次迭代后，若某聚类中心俘获的数据样本个数小于预定的阈值minL，则 剔除该类。
[0086] 步骤5:更新聚类个数为C1+1，并根据数据的隶属关系重新计算新的聚类中心W1+1,。 (c = l ,2, ··· ,Ci+i) 〇
[0087] 步骤6:如果算法满足收敛条件则结束，否则返回步骤2,进行下一次迭代计算。收 敛条件有:两次迭代中的聚类中心距离的变化小于一个很小的阈值ε，或者每个聚类中数据 Dk到聚类中心的距离平方和以及聚类之间的距离平方和达到最小。
[0088]本实例中，初始聚类中心个数Co取为15，最小距离阈值Θ取为300，最小数据量阈值minL 取为100,阈值ε取为0.0001。通过自适应聚类，稳定工况有效数据被划分为7个聚类(即C=7)，获得每 一种稳定工况的建模数据
[0089] (4)稳定工况统计监测建模
[0090]对每一种稳定工况的建模数#
进行均值中心化和方 差归一化的数据预处理，获得处理后的建模数据Xr(KrmXJ)。对建模数据Xr进行PCA分解，得 到对应的主元个数Ar、负载矩阵Pr、控制限和，建立对应第r个稳定工况的PCA监 测模型。
[0091]对建模数据预处理包括依次进行的减均值、除标准差处理获得Xr (KrmX J)，计算公式如下： ?
[0093] 其中足υ表示建模数?： 第rk行第j列的数据，下标r代表工况数、k 代表样本数、j代表变量数，.....-况第j个变量的均值表示第r个稳定 工况第j个变量的标准差。
[0094] 将数据Xr进行PCA分解，PCA分解的计算公式为：
[0096] 其中Tr为主元子空间中（KrmXAr)维的得分矩阵;Pr为主元子空间中（JXAr)维的负 载矩阵;Er为残差子空间中（K rmXJ)维的残差矩阵。Ar表示主元个数，由累计贡献率方法确 定。
[0097] 离线计算主成分子空间的Hotelling-T2和残差子空间的SPE(Squared Pred i c t i on Error)监测统计量。监测统计量的计算公式为：
[0099] 其中trk = XrkPr为（1 X Ar )维的主元得分向量，Xrk为建模数据Xr的第k行，Pr为主元 子空间的负载矩阵，对角矩阵& = …，\ )由建模数据Xr的协方差矩阵Σ XrTXr的 前Ar个特征值所构成。
[0100] SPEr监测统计量的计算公式为：
[0102] 其中
表示重构得到的（IX J)维估计向量。
[0103] 根据F分布计算f监测统计量的控制限#%为：
[0105] 其中α为置信度，Fa(Ar，Krm-Ar)是对应于置信度为α，自由度为A r、K?-Ar条件下的F 分布临界值。
[0106] 根据X2分布计算SPEr监测统计量的控制阳
[0108] 其：
，_，nr、vr分别为SPE r监测统计量的均值和方差。
[0109] 本实例中，7个PCA监测模型的主元个数依次为A! = 26、A2 = 21、A3 = 26、A4 = 26、A5 = 27、六6 = 26^7 = 26，则?1为(44\26)负载矩阵、？2为(44\21)负载矩阵、？3为(44\26)负载矩 阵、P 4为(44X26)负载矩阵、Ρ^(44 X 27)负载矩阵、P6为(44X26)负载矩阵、P7为(44X26) 负载矩阵。统计控制限的置信度α取为0.99,监测模型1的控制限和 c?rSPE( =.16:·3218、监测模型 2 的控制限 cir# =40.88:73:和 cirSPt: =8,5462、监测模型 3 的 控制限改45.9832和你啊=14·9丨56、监测模型4的控制限g%. = 45.9867和 cirSPE4 =23.0371、监测模型5的控制限冰> =47.2113和=51.4481.、监测模型6 的控制限= 46.0821和C%PE:6 = 12.6:8:84、监测模型7的控制限= 45.9559和 cthm7 = I B.5386
[0110] (5)多工况过程状态监测
[0111] 多工况过程状态监测时，计算当前滑动时间窗口内数据的稳定度因子SFnew，通过 与给定阈值ζ的比较确定当前所处的工况类型。若为过渡工况，则将当前时刻7；i和SPEnew 两个监测统计量赋值为0;否则采用平滑滤波方法计算当前滑动时间窗口内的移动平均数 据，并计算其与各个聚类中心的欧式距离，调用距离最小的聚类所对应的PCA监测模型在线 雅4邮PE-wW傾灘i十1：。娜W侦i十胃雛刊空·艮眺綱#?紅 常，如果至少其中一个统计量超出控制限，表明设备状态异常。
[0112] 计算处于稳定工况滑动时间窗口的移动平均数据与各个聚类中心的欧式距离，根 据距离最小判断该时刻对应的稳定工况监测模型。欧式距离的计算公式为：
[0114] 其中Xne3W,j表示移动平均数据的第j个变量的检测值，Cem.j表示第i个聚类中心第j 个变量的值，C表示聚类中心的个数。
[0115] 设当前移动平均数据对应第L(1<L<C)个稳定工况监测模型，利用该监测模型建 模数据的均值忑，_/和标准差毛，_/对移动平均数据进行标准化处理，获得xl_(1 X J)。将其投 影到第L个监测模型上，计算和SPEL_监测统计量指标。
[0117]其中Pl表示第L个监测模型主元子空间的负载矩阵，对角矩阵SL是由建模数据XL的 协方差矩阵的前Al个特征值所构成
表示重构得到的（1 XJ)维估计向 量。
[0118] 如果两个监测统计量都位于控制限内，表明超高速小盒包装机处于正常状态，如 果其中至少一个监测统计量超出控制限，表明设备处于异常状态。
[0119] 本实例中，选取多工况过程8816组采样数据％ew(8816X44)进行分析，该测试数据 T2和SPE统计量的监测结果如图6所示。引入首次故障i警时间，定义为连续9个监测统计量 超过对应控制限的采样时间点。在两个监测子空间中有个别样本点的监测统计量超出控制 限，根据首次故障报警的定义，两个监测统计量均未发生故障报警，对当前过程的正常状态 做出了正确的指示，表明所提状态监测与诊断方法具有准确监测多工况过程正常状态的能 力。
[0120]作为比较，重新选取多工况过程4505组采样数据％ew(4505X44)进行状态监测与 诊断的分析。在第1001个样本点的时候，分别令热熔胶温i阶跃增长1% (突变故障）、A1电 机电流阶跃降低30% (突变故障）、A5A电机温度以0.03的斜率增加(渐变故障），一直持续到 第2000个样本点的时候恢复正常，形成故障状态数据1、2和3。故障状态数据T 2和SPE统计量 的监测结果分别如图7至图9所示。
[0121] 对于故障1，主元子空间的T2统计量出现连续超限，主成分子空间出现异常，检测 出的首次故障报警时间为第1030个样本点时刻，即监测模型延后290秒准确检测出故障的 发生;对于故障2,主元子空间的T 2统计量出现连续超限，主成分子空间出现异常，检测出的 首次故障报警时间为第1030个样本点时刻，即监测模型延后290秒准确检测出故障的发生； 对于故障3,主元子空间的T 2统计量出现连续超限，主成分子空间出现异常，检测出的首次 故障报警时间为第1318个样本点时刻，即监测模型延后3170秒准确检测出故障的发生。说 明所提方法能对多工况过程的突变和渐变故障做出正确指示，表明该方法具有准确监测多 工况过程故障状态的能力。
[0122] (6)基于贡献图的故障诊断
[0123] 当任一监测统计量超出控制限时，计算各个过程变量对超限统计量的贡献值，确 定引起异常的原因变量。主成分子空间，监测统计量的超出正常控制限时，主成分得分 变量tu
贡献率可以计算如下：
[0125] 其中tu表示对应第L个监测模型的第a个主成分得分变量，表示对应第L个监测 模型的第a个特征值。
[0126] 过程变量对tu的贡献率可以计算如下：
[0128]其中plm表示对应第L个监测模型的第a个负载向量的第j个变量 。
[0129] 当残差子空间SPELnew监测统计量的超出正常控制限时，过程变量xLnew,j对SPE Lnew的 贡献率可以计算如下：
L,n(￡W
[0131] 其中
.示重构估计向I
:勺第j个变量，sign
I示残 差的正负信息。
[0132] 本实例中，根据图7的监测结果，故障发生在主成分子空间中。在首次故障报警的 第1030个样本点时刻，过程变量对7^0，监测统计量超限的贡献率如图10所示，确定引起故 障的原因变量为第43个变量，即H1000热熔胶温度。
[0133] 根据图8的监测结果，故障发生在主成分子空间中。在首次故障报警的第1030个样 本点时刻，过程变量对^^胃监测统计量超限的贡献率如图11所示，确定引起故障的原因变 量为第1个变量，即A1电机电流。
[0134] 根据图9的监测结果，故障发生在主成分子空间中。在首次故障报警的第1318个样 本点时刻，过程变量对监测统计量超限的贡献率如图12所示，确定引起故障的原因变 量为第16个变量，即A5A电机温度。
[0135] 以上所述仅为本发明的较佳实施举例，并不用于限制本发明，凡在本发明精神和 原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
【主权项】
1. 一种具有多工况特性的超高速小盒包装机状态监测与诊断方法，其特征在于，包括 步骤： 1) 根据超高速小盒包装机HlOOO的运行机理确定统计模型的输入变量，采集相应的运 行数据，获得二维原始数据旮(Kx·/)，K为采样点个数，J为监测变量个数； 2) 基于车速检测值，采用滑动时间窗口 H方法计算稳定度因子SF，根据稳定度因子SF与 阈值ζ的关系判断是否处于稳定工况，采用平滑滤波方法对处于稳定工况的滑动时间窗口 H 内原始数据进行算术平均，获得一组移动平均数据，将移动平均数据中"车速"变量值大于 阈值Φ的数据作为稳定工况的有效数据D= (D1 ,D2，…，DN}，Ν为有效数据的个数； 3) 采用自适应k-means聚类方法对有效数据的相似度进行聚类分析，获得每一种稳定 工况的建模数据C表示将稳定工况划分成C个聚类，每个聚类中具 有相似过程特性的建模数据用一统计监测模型来描述； 4) 对每一种稳定工况的建模数据，进行均值中心化和方差 归一化的数据预处理，获得处理后的建模数据Xr (Krm X J )，对建模数据Xr进行PCA分解，得到 对应的主元个数Ar、负载矩阵Pr、控制限和&胃，建立对应第 r个稳定工况的PCA监测 丰旲型； 5) 实时计算当前滑动时间窗口内数据的稳定度因子SFne3w,通过与给定阈值ζ的比较确 定当前所处的工况类型;若为过渡工况，则将当前时刻:ζ?,和SPE new两个监测统计量赋值为 〇;否则采用平滑滤波方法计算当前滑动时间窗口内的移动平均数据，并计算其与各个聚类 中心的欧式距离，调用距离最小的聚类所对应的PCA监测模型在线计算T^jPSPE new两个监 测统计量，若一个监测统计量超出控制限，表明设备状态异常； 6) 当检测到设备出现有异常时，根据超限的监测统计量处于T2的主元子空间或是SPE的 残差子空间，计算每个过程变量对超限统计量的贡献率，贡献率较大的变量被初步确定为 造成设备异常的原因变量。2. 如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，在步骤2)中采用滑动时间窗口 H方法计算稳定度因子SF进行工况的识别，计算公 式如下：其中Vh max表示滑动时间窗口H内车速数据的最大值，Vh min表示滑动窗口H内车速数据的 最小值，Vmax表示设备允许的车速最大值，H表示滑动时间窗口的长度。3. 如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，在步骤3)中采用自适应k-means聚类方法进行工况的划分，对有效数据的相似度 进行聚类分析，通过定义两个聚类中心的最小距离的阈值Θ，同时最小化同一聚类数据之间 的局部距离平方和以及所有聚类之间的全局距离平方和，获得最优的聚类数目和聚类中 心，实现对稳定工况的自适应划分。4. 如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，在步骤4)中对建模数据预处理包括依次进行的减均值、除标准差 处理获得以^】），计算公式如下：其中，表示建模数据中第rk行第j列的数据，下标r代表工况、k代表样 本、j代表变量，巧,7:表示第r个稳定工况第j个变量的均值表示第r个稳定工况第j个变 量的标准差。5. 如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，所述步骤4)中PCA分解的计算公式为：其中，Tr为主元子空间中（KrmXAr)维的得分矩阵;Pr为主元子空间中（JXAr)维的负载 矩阵;Er为残差子空间中(KrmX J)维的残差矩阵;Ar表示主元个数，由累计贡献率方法确定； 离线计算主成分子空间的HotelIing-T2和残差子空间的SPE(Squared Prediction Error)监测统计量；I；2监测统计量的计算公式为：其中trk = XrkPr为（I X Ar )维的主元得分向量，Xrk为建模数据Xr的第k行，Pr为主元子空间 的负载矩阵，对角矩阵由建模数据Xr的协方差矩阵的前Ar 个特征值所构成； SPEr监测统计量的计算公式为：其中表示重构得到的（I XJ)维估计向量； 根据F分布计算?；2监测统计量的控制限为：其中α为置信度，Fa(Ar，Krm-Ar)是对应于置信度为α，自由度为A r、Krm-Ar条件下的F分布 临界值； 根据X2分布计算SPEr监测统计量的控制限#rspi为··；/?,、Vi.分别为SPEr监测统计量的均值和方差。6. 如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，在步骤5)中，欧式距离的计算公式为：其中xn?,j表示移动平均数据的第j个变量的检测值，Cem, j表示第i个聚类中心第j个变 量的值，C表示聚类中心的个数。7. 如权利要求6所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，设当前移动平均数据对应第L(I^LSC)个稳定工况监测模型，利用该监测模型 建模数据的均值&,,和标准差^^对移动平均数据进行标准化处理，获得an?(l X J)，将其 投影到第L个监测模型上，计算J^jPSPEL_监测统计量指标；其中Pl表示第L个监测模型主元子空间的负载矩阵，对角矩阵Sl是由建模数据Xl的协方 差矩阵的前Ai个特征值所构成，表示重构得到的（I XJ)维估计向量。8. 如权利要求1所述的具有多工况特性的超高速小盒包装机状态监测与诊断方法，其 特征在于，在步骤6)中，主成分子空间?^_监测统计量的超出正常控制限时，主成分得分变 散3对jS-颜鮮卿雅女吓：其中Ua表示对应第L个监测模型的第a个主成分得分变量，表示对应第L个监测模型 的第a个特征值； 过程变量XkW, j对tu的贡献率可以计算如下：其中PUd表示对应第L个监测模型的第a个负载向量的第j个变量； 当残差子空间SPELnew监测统计量的超出正常控制限时，过程变量xLnew,j对SPE Lnew的贡献 率可以计算如下：其中表示重构估计向量.的第j个变量，1表示残差的正 负信息。
【专利摘要】本发明公开了一种具有多工况特性的超高速小盒包装机状态监测与诊断方法，该方法针对运行过程多工况特性及工况普遍未知的问题，基于滑动窗口计算稳定度因子识别不同工况，采用自适应聚类方法对稳定工况进行划分形成若干稳定工况数据簇，利用PCA方法为每个稳定工况数据簇建立监测模型。多工况过程状态监测时，根据当前滑动窗口内数据的稳定度因子判断工况类型，若为过渡工况，则将监测统计量赋值为0；若为稳定工况，计算当前有效数据与各个聚类中心的距离，获得匹配的PCA监测模型进行实时监测，任一统计量超限时采用贡献图确定引起故障的主要过程变量。
【IPC分类】G05B23/02
【公开号】CN105487524
【申请号】CN201511007689
【发明人】王伟, 赵春晖, 楼卫东, 李钰靓, 张利宏, 熊月宏
【申请人】浙江中烟工业有限责任公司
【公开日】2016年4月13日
【申请日】2015年12月29日
最新回复(0)