一种基于PSO-AR模型误差Kiviat图的故障检测方法

本发明涉及一种工业过程故障检测方法，特别涉及一种基于pso-ar模型误差kiviat图的故障检测方法。

背景技术：

1、工业过程的故障检测是工业生产的永恒话题，只有在工业过程运行于无故障的条件下，才能保证产品品质和生产安全。随着工业智能制造的大力推进，工业生产过程安装配套有先进的测量仪表和计算机辅助系统，能够实时的提供关于生产过程的各类状态数据。这些海量的状态数据中蕴藏有能反映故障的潜在有用信息，如何充分挖掘并有效利用这些信息实现对工业过程的故障实时检测是实现智能制造不可避免的问题。在已有的技术文献资料中，以多变量统计分析算法为代表的数据驱动建模技术已经在工业过程，特别是流程工业过程中得到的广泛的关注，其技术特点都是在数据特征分析提取的基础上，通过超椭球或超球体的方式定义出正常变化范围。

2、从现阶段工业过程中测量仪表采集数据的特点来看，采样数据间隔时间短势必导致海量的数据，并同时要求数据的单位处理时间更短。因此，复杂的数据运算与处理过程是不适合现代工业过程的故障检测的。此外，考虑到工业生产过程的持续平稳性是保证无故障运行的直接手段，采样数据体现在时间先后上的平稳变化是反映工业过程运行正常最直接的方式。从这个角度来讲，监测工业过程采样数据的时序关系的变化可实现故障检测目的。为挖掘这种采样数据的时序关系特征，建立自回归(auto-regressive，缩写：ar)模型是最有效的实现手段，而相应的模型误差即可反映采样数据的正常变化情况。

3、在建立ar模型时，传统的思路都是通过最小二乘或偏最小二乘算法来训练得到，原采样数据的耦合导致奇异性，潜特征个数的确定等都会对相应ar模型的精度造成影响。另一方面，传统的方法技术对模型误差的监测一般都是依赖于距离统计量，如平方马氏距离和平方欧式距离，这是无法有效的表征各个存在交互关系的误差变化特征。kiviat图作为一种能直接展示多维数据变化的技术手段，不仅能展示各个数据的变化，而且能兼顾反映多维数据之间的交互性变化特征。虽然，kiviat图已经在可视化数据分析领域得到了广泛的应用，但还未应用于工业过程的故障检测中。

技术实现思路

1、本发明所要解决的主要技术问题是：如何利用粒子群优化(particle swarmoptimization,缩写：pso)算法直接优化搜寻ar模型的回归系数矩阵，并利用kiviat图来监测模型误差的变化实现对工业过程的故障检测。具体来讲，本发明方法首先利用pso算法为工业过程的采样数据建立能体现采样时间先后关系的ar模型，然后利用kiviat图来实时监测工业过程的在线采样数据对应的模型误差的变化，从而完成故障检测任务。

2、本发明方法解决上述问题所采用的技术方案为：一种基于pso-ar模型误差kiviat图的故障检测方法，具体包括如下所示步骤。

3、步骤(一)，pso-ar模型训练前的准备，具体包括以下三个方面。

4、首先，设置pso-ar模型训练所需的参数，具体包括：工业过程每个采样时刻都可采集的测量变量个数m，时序相关阶数c，粒子总数n，最大迭代次数d，加速度因子α1和α2，惯性权重的初始值w1和终值w2。

5、其次，确定惯性权重w与迭代次数d之间的关系，具体是w＝w1-(w1-w2)(d/d)2。

6、最后，确定pso算法中粒子的表现形式，具体是：使用一个m×m维的实数矩阵表示粒子的位置，使用一个m×m维的实数矩阵表示粒子的速度；其中，m＝c·m。

7、步骤(二)，利用工业过程的采样数据训练pso-ar模型，保留回归系数矩阵θ，具体包括以下所示步骤(1)至步骤(3)。

8、步骤(1)：确定用于pso-ar模型训练的输入矩阵z和输出矩阵y，具体的实施过程如步骤(1.1)至步骤(1.3)所示。

9、步骤(1.1)：从工业过程的历史采样数据库中，获取其运行于正常状态下的n组采样数据，并将每组采样数据表示为一个m×1维的数据向量，从而对应得到n个m×1维的数据向量z1,z2,…,zn；其中，每组采样数据都是由m个测量变量对应的采样值组成。

10、步骤(1.2)：根据公式对z1,z2,…,zn中的各个元素分别实施标准化处理，对应得到n个m×1维的数据向量其中，和分别表示和zi中的第j个元素，μj表示z1,z2,…,zn中第j个元素的平均值，δj表示z1,z2,…,zn中第j个元素的标准差，j＝1,2,…,m，下标号i＝1,2,…,n。

11、步骤(1.3)：将合并成一个(n-c)×m维的输出矩阵并按照如下所示公式①组成一个(n-c)×m维的输入矩阵z：

12、

13、其中，上标号t表示矩阵或向量的转置。

14、步骤(2)：根据输入矩阵z和输出矩阵y确定出粒子的位置上限pmax和位置下限pmin，粒子的速度上限vmax和速度下限vmin，具体的确定过程如步骤(2.1)至步骤(2.3)所示。

15、步骤(2.1)：根据公式b0＝(ztz)-1zty计算最小二乘回归系数矩阵b0后，将b0中元素的最大值和最小值分别记录为bmax和bmin。

16、步骤(2.2)：确定粒子的位置上限pmax＝max{bmax,bmin|}和粒子的位置下限pmin＝min{bmin,-bmax}；其中，max{bmax,|bmin|}表示选择bmax和|bmin|两者中的最大值，|bmin|表示计算bmin的绝对值，min{bmin,-bmax}表示选择bmin和-bmax两者中的最小值。

17、步骤(2.3)：确定粒子的速度上限vmax＝(pmax-pmin)/10和粒子的速度下限vmin＝-vmax。

18、步骤(3)：利用pso算法确定回归系数矩阵θ，从而完成对pso-ar模型y＝zθ+e的训练；其中，e表示模型误差矩阵，具体的实施过程包括以下所示步骤(3.1)至步骤(3.7)。

19、步骤(3.1)：从区间[pmin,pmax]上随机生成n个m×m维的实数矩阵b1,b2,…,bn，分别对应表示n个粒子的位置，从区间[vmin,vmax]上随机生成n个m×m维的实数矩阵v1,v2,…,vn，分别对应表示n个粒子的的速度，再设置迭代次数d＝1。

20、步骤(3.2)：根据分别计算各个粒子的适应度值j1,j2,…,jn；其中，jk表示第k个粒子的适应度值，m×m维的实数矩阵bk表示第k个粒子的位置，表示计算y-zbk中所有元素的平方和，粒子编号k＝1,2,…,n。

21、步骤(3.3)：将j1,j2,…,jn中最小的适应度值及其对应的粒子编号分别标记为和后，再判断d是否等于1；若是，则分别设置全局最优位置和最佳适应度值后，再执行步骤(3.5)；若否，则执行步骤(3.4)；其中，表示第个粒子的位置。

22、步骤(3.4)：判断jg是否大于若是，则设置全局最优粒子及其对应的适应度值后，执行步骤(3.5)；若否，则直接执行步骤(3.5)。

23、步骤(3.5)：根据公式分别对速度矩阵v1,v2,…,vn实施更新后，再根据公式vk(g,j)＝max{min{vk(g,j),vmax},vmin}对速度矩阵中各个元素分别进行修正。

24、其中，λk和ηk是从区间[0,1]上生成的两个随机数，vk(g,j)表示vk中的第g行第j列的元素，max{bmax,|bmin|}n

25、步骤(3.6)：根据公式bk＝bk+vk分别对b1,b2,…,bn实施更新后，再根据公式bk(g,j)＝max{min{bk(g,j),pmax},pmin}对b1,b2,…,bn中的各个元素分别进行修正；其中，bk(g,j)表示bk中的第g行第j列的元素，j＝1,2,…,m。

26、步骤(3.7)：设置d＝d+1后，判断d是否大于d；若否，则返回步骤(3.2)；若是，则设置回归系数矩阵θ＝g后，计算模型误差矩阵e＝y-zθ，从而完成对pso-ar模型y＝zθ+e的训练。

27、步骤(三)，确定kiviat图监测误差变化的实施方式，具体包括以下所示步骤(4)至步骤(6)。

28、步骤(4)：确定1×m维的误差向量ε转变至kiviat图的转变方式，具体包括如下所示步骤(4.1)至步骤(4.2)。

29、步骤(4.1)：误差向量ε中的第1个至第m个元素依次标记为e1,e2,…,em后，分别根据公式xj＝|ej|×cos[(j-1)θ]和yj＝|ej|×sin[(j-1)θ]计算各个元素对应的kiviat图横坐标x1,x2,…,xm和kiviat图纵坐标y1,y2,…,ym；其中，夹角θ＝360°/m，ej表示误差向量ε中的第j个元素，cos和sin分别表示余弦和正弦函数，|ej|表示计算ej的绝对值，j＝1,2,…,m。

30、步骤(4.2)：根据公式和计算误差向量ε转变至kiviat图的中心点横坐标x和中心点纵坐标

31、步骤(5)：确定kiviat图监测误差变化的上限值具体包括以下所示步骤(5.1)至步骤(5.3)。

32、步骤(5.1)：依次设置误差向量ε等于模型误差矩阵e中的第1行向量至第n-c行向量，并同时根据步骤(4)中确定的转变方式，对应得到相应的中心点横坐标和中心点纵坐标，依次分别标记为和其中，和分别表示当ε等于模型误差矩阵e中的第1行向量时，转变至kiviat图的中心点横坐标和中心点纵坐标；和分别表示当ε等于模型误差矩阵e中的第2行向量时，转变至kiviat图的中心点横坐标和中心点纵坐标；和分别表示当ε等于模型误差矩阵e中的第n-c行向量时，转变至kiviat图的中心点横坐标和中心点纵坐标。

33、步骤(5.2)：按照如下所示公式②将和合并成一个(n-c)×2维的坐标值矩阵p后，再计算p中所有行向量的平均值向量u：

34、

35、步骤(5.3)：计算协方差矩阵s＝(p-u)t(p-u)/(n-c-1)，再计算对称矩阵φ＝(p-u)s-1(p-u)t后，将φ中对角线上最大的元素确定为上限值其中，(n-c)×2维的平均值矩阵u中的各行向量都等于u，上标号t表示矩阵或向量的转置。

36、步骤(四)，对工业过程实施在线故障检测，具体包括如下所示步骤(6)至步骤(10)。

37、步骤(6)：在工业过程的最新采样时刻t，获取由m个测量变量对应的采样值组成的一组采样数据，并将该组采样数据表示为一个m×1维的数据向量zt后，再根据公式对zt的各个元素分别实施标准化处理，对应得到一个m×1维的数据向量其中，和分别表示zt和中的第j个元素。

38、步骤(7)：获取t之前的c个采样时刻的数据向量zt-1,zt-2,…,zt-c，并对其实施与步骤(6)中相同的标准化处理，从而对应得到c个m×1维的数据向量再合并组成一个1×m维的输入向量

39、步骤(8)：根据公式生成误差向量ε后，执行步骤(4)中确定的转变方式，得到相应的中心点横坐标和中心点纵坐标

40、步骤(9)：将和组成一个1×2维的坐标值向量后，再根据公式计算得到最新采样时刻t的误差变化指标

41、步骤(10)：判断是否满足条件若是，则工业过程运行中未出现故障，返回步骤(6)继续利用最新采样时刻的一组采样数据实施在线故障检测；若否，则工业过程运行出现故障，触发故障警报后，再返回步骤(6)继续利用最新采样时刻的一组采样数据实施在线故障检测。

42、通过以上所述实施步骤，本发明方法的优势介绍如下。

43、本发明方法通过pso这种智能优化算法为工业过程的时序采样数据建立相应的ar模型，然后利用kiviat图监测通过这种pso-ar模型生成的误差的变化，实现对工业过程的故障检测目的。相比于传统方法，本发明方法能从最优化的角度完成对ar模型的训练而不必考虑奇异值、过拟合和欠拟合等问题。此外，本发明方法利用kiviat图监测误差之间的交互性变化情况，也比传统通过距离的方式检测故障更加全面和直观。因此，本发明方法是一种更为优选的故障检测方法。