CN102750445B - 基于复化Simpson公式改进多变量灰色模型的故障预测方法 - Google Patents

Abstract

基于复化Simpson公式改进多变量灰色模型的故障预测方法，属于信号处理领域，本发明为解决现有故障预测设备系统的组成具有复杂性、结构关系的模糊性以及特征参数获取的不确定性、不完整性的问题。本发明方法包括以下步骤：步骤一、输入待测设备的故障原始特征数据序列和故障相关数据序列，计算其累加生成数；步骤二、建立多变量灰色模型，计算灰色系数；步骤三、利用步骤二求得的灰色系数获取的预测值步骤四、获得的预测值步骤五、计算多变量灰色模型的预测精度；步骤六、判断待测设备是否有故障。

Description

基于复化Simpson公式改进多变量灰色模型的故障预测方法

技术领域

本发明涉及基于复化Simpson公式改进多变量灰色模型的故障预测方法，属于信号处理领域。

背景技术

故障预测为制定精准的设备维修保障计划，减少维修保障费用，提高设备完好率和任务成功率，避免因故障造成巨大损失等具有重大意义。然而，实际设备系统组成的复杂性、结构关系的模糊性以及特征参数获取的不确定性、不完整性，成为故障预测的难点。

发明内容

本发明目的是为了解决现有故障预测设备系统的组成具有复杂性、结构关系的模糊性以及特征参数获取的不确定性、不完整性的问题，提供了一种基于复化Simpson公式改进多变量灰色模型的故障预测方法。

本发明所述基于复化Simpson公式改进多变量灰色模型的故障预测方法，该方法包括以下步骤：

步骤一、输入的待测设备的故障原始特征数据序列为

$X_1^{\left(0\right)}=\{x_1^{\left(0\right)}\left(1\right),x_1^{\left(0\right)}\left(2\right),...,x_1^{\left(0\right)}\left(k\right),...,x_1^{\left(0\right)}\left(n\right)\},$ k＝1，2，...，n    (1)

且 $x_1^{\left(0\right)}\left(k\right)\≥0;$

输入的待测设备的故障相关数据序列为

$X_i^{\left(0\right)}=\{x_i^{\left(0\right)}\left(1\right),x_i^{\left(0\right)}\left(2\right),...,x_i^{\left(0\right)}\left(k\right),...,x_i^{\left(0\right)}(n+m)\},$ i＝2，3，...，N，t＝1，2，...，n+m    (2)

且 $x_i^{\left(0\right)}\left(t\right)\≥0,$

其中，m为待预测数据个数，

由i＝2，3，...，N决定，

通过公式 $x_1^{\left(1\right)}\left(k\right)=\underset{l=1}{\overset{t}{\mathrm{\Σ}}}{x}_1^{\left(0\right)}\left(l\right)---\left(3\right)$

和

公式 $x_i^{\left(1\right)}\left(t\right)=\underset{l=1}{\overset{k}{\mathrm{\Σ}}}{x}_i^{\left(0\right)}\left(l\right)---\left(4\right)$

计算的一阶累加生成数

步骤二、以为初值，建立如下多变量灰色模型S-GMC(1，N)：

$\frac{{\mathrm{dx}}_1^{\left(1\right)}\left(t\right)}{\mathrm{dt}}+{\mathrm{\αx}}_1^{\left(1\right)}\left(t\right)=\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(t\right)+{\mathrm{\β}}_N---\left(5\right)$

由κ＝[α β₁ β₂ ...β_N]^T＝(B^TB)^-1B^TY，计算S-GMC(1，N)的灰色系数κ，其中

$Y={\left[\begin{array}{cccc}{x}_1^{\left(0\right)}\left(2\right)+x_1^{\left(0\right)}\left(3\right)& x_1^{\left(0\right)}\left(3\right)+x_1^{\left(0\right)}\left(4\right)& ...& x_1^{\left(0\right)}(n-1)+x_1^{\left(0\right)}\left(n\right)\end{array}\right]}^T---\left(6\right)$

步骤三、利用步骤二求得的灰色系数，就可以利用下式对进行预测，获取的预测值

步骤四、由步骤三的并按下式获得的预测值

${\hat{x}}_1^{\left(0\right)}\left(t\right)={\hat{x}}_1^{\left(1\right)}\left(t\right)-{\hat{x}}_1^{\left(1\right)}(t-1)---\left(9\right)$

其中，初始化 ${\hat{x}}_1^{\left(0\right)}\left(1\right)=x_1^{\left(0\right)}\left(1\right);$

步骤五、根据公式

$\mathrm{RMSPEPR}=\sqrt{\frac{1}{n}\underset{t=1}{\overset{n}{\mathrm{\Σ}}}{[{\hat{x}}_1^{\left(0\right)}\left(t\right)-x_1^{\left(0\right)}\left(t\right)]}^2/{\left[x_1^{\left(0\right)}\left(t\right)\right]}^2}\×100\%---\left(10\right)$

$\mathrm{RMSPEPO}=\sqrt{\frac{1}{m}\underset{t=n+1}{\overset{n+m}{\mathrm{\Σ}}}{[{\hat{x}}_1^{\left(0\right)}\left(t\right)-x_1^{\left(0\right)}\left(t\right)]}^2/{\left[x_1^{\left(0\right)}\left(t\right)\right]}^2}\×100\%---\left(11\right)$

计算S-GMC(1，N)的预测精度，

其中，RMSPEPR为由前部分样本，即t＝1，2，...，n组成的根平均平方误差百分比，而RMSPEPO为由后部分样本，即t＝n+1，n+2，...，n+m组成的根平均平方误差百分比；

步骤六、判断条件是否满足，其中ε₁为故障阈值上界，ε₂为故障阈值下界，

若满足，则说明待测设备无故障；否则，判定待测设备有故障。

本发明的优点：

1)在计算灰色系数的过程中，本发明用代数精度较高的Simpson公式代替梯形公式，对微分方程进行有效的离散化，最终得到准确的灰色系数；

2)在计算卷积积分的过程中，本发明方法用复化Simpson公式代替复化梯形公式，并把t≥5的数据分为奇数和偶数单独处理，最终得到精确的故障预测值；

3)本发明不改变原有多变量灰色模型故障预测的框架，计算过程简单，易于实现，有利于推广应用。

附图说明

图1是本发明所述一种基于复化Simpson公式改进多变量灰色模型的故障预测流程图；

图2是本发明方法与GMC(1，N)求得的t＝1，2，...，10的故障预测值及原始值对比图；

图3是本发明方法与GMC(1，N)求得的t＝1，2，...，10的故障预测值及原始值对比图；

图4是本发明方法与GMC(1，N)的故障预测精度(包括RMSPEPR和RMSPEPO)对比图。

具体实施方式

具体实施方式一：下面结合图1至图4说明本实施方式。

灰色系统理论以“部分信息已知、部分信息未知”的“贫信息”、“小样本”不确定性系统为研究对象，通过对“部分”已知信息的生成、开发、提取有价值的信息，实现对系统运行行为、演化规律的正确描述和有效监控，为故障预测提供了一种有效途径。

作为灰色系统理论的一个重要内容，灰色模型(Grey Model，以下简称GM)具有很多优良特性，如：1)不改变原始信号特性；2)有效预测非线性或非平稳信号；3)计算量小，速度快；4)允许少量数据预测；5)短期预测精度高等，这些特性促使其在经济、社会、农业、工程等领域尤其是故障预测方面，得到日益广泛的应用。需要指出的是，GM是一种特殊的微分方程，可以统称为GM(P，Q)，其中G代表Grey，M代表Model，P代表GM是P阶方程，Q代表GM有Q个变量。最常用的是一维情形的GM(1，1)(P＝1，Q＝1)，而本发明要探讨的多变量GM可以表示为GM(1，N)(P＝1，Q＝N)。

遗憾的是，目前GM的处理对象一般都是用一维状态变量描述，而在二维乃至多维状态变量描述的领域中，GM的应用却很少见。究其原因，主要是GM(1，N)的发展还很不完善，预测能力远不及GM(1，1)。在已有的GM(1，N)方法中，效果最好的要数台湾国立虎尾科技大学的Tzu Li Tien教授研究的一种基于卷积积分(Convolution Integral)的多变量GM(以下简称GMC(1，N))。通过计算具有特殊定义的卷积积分，GMC(1，N)不再把原始GM(1，N)(Original GM(1，N)，以下简称OGM(1，N))中难于求积分的非常数项笼统地设定为常数，从而有效地提高了GM(1，N)的预测精度。然而，OGM(1，N)在计算灰色系数时采用了代数精度较低的梯形公式，而GMC(1，N)在计算灰色系数和卷积积分的过程中，分别采用了代数精度较低的梯形公式和复化梯形公式，阻碍了预测精度的进一步提高。

GM(1，N)在理论上的缺陷阻碍了其在故障预测方面的深入应用，虽然很多学者已经认识到了该问题的严重性，但目前尚无一套良好的解决办法，这就更加凸显出提高GM(1，N)预测精度的必要性。

本实施方式所述基于复化Simpson公式改进多变量灰色模型的故障预测方法，该方法包括以下步骤：

步骤一、输入的待测设备的故障原始特征数据序列为

$X_1^{\left(0\right)}=\{x_1^{\left(0\right)}\left(1\right),x_1^{\left(0\right)}\left(2\right),...,x_1^{\left(0\right)}\left(k\right),...,x_1^{\left(0\right)}\left(n\right)\},$ k＝1，2，...，n    (1)

且 $x_1^{\left(0\right)}\left(k\right)\≥0;$

输入的待测设备的故障相关数据序列为

$X_i^{\left(0\right)}=\{x_i^{\left(0\right)}\left(1\right),x_i^{\left(0\right)}\left(2\right),...,x_i^{\left(0\right)}\left(k\right),...,x_i^{\left(0\right)}(n+m)\},$ i＝2，3，...，N，t＝1，2，...，n+m    (2)

且 $x_i^{\left(0\right)}\left(t\right)\≥0,$

其中，m为待预测数据个数，

由i＝2，3，...，N决定，

通过公式 $x_1^{\left(1\right)}\left(k\right)=\underset{l=1}{\overset{t}{\mathrm{\Σ}}}{x}_1^{\left(0\right)}\left(l\right)---\left(3\right)$

和

公式 $x_i^{\left(1\right)}\left(t\right)=\underset{l=1}{\overset{k}{\mathrm{\Σ}}}{x}_i^{\left(0\right)}\left(l\right)---\left(4\right)$

计算的一阶累加生成数

步骤二、以为初值，建立如下多变量灰色模型S-GMC(1，N)：

$\frac{{\mathrm{dx}}_1^{\left(1\right)}\left(t\right)}{\mathrm{dt}}+{\mathrm{\αx}}_1^{\left(1\right)}\left(t\right)=\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(t\right)+{\mathrm{\β}}_N---\left(5\right)$

由κ＝[α β₁ β₂ ...β_N]^T＝(B^TB)^-1B^TY，计算S-GMC(1，N)的灰色系数κ，其中

$Y={\left[\begin{array}{cccc}{x}_1^{\left(0\right)}\left(2\right)+x_1^{\left(0\right)}\left(3\right)& x_1^{\left(0\right)}\left(3\right)+x_1^{\left(0\right)}\left(4\right)& ...& x_1^{\left(0\right)}(n-1)+x_1^{\left(0\right)}\left(n\right)\end{array}\right]}^T---\left(6\right)$

步骤三、利用步骤二求得的灰色系数，就可以利用下式对进行预测，获取的预测值

步骤四、由步骤三的并按下式获得的预测值

${\hat{x}}_1^{\left(0\right)}\left(t\right)={\hat{x}}_1^{\left(1\right)}\left(t\right)-{\hat{x}}_1^{\left(1\right)}(t-1)---\left(9\right)$

其中，初始化 ${\hat{x}}_1^{\left(0\right)}\left(1\right)=x_1^{\left(0\right)}\left(1\right);$

步骤五、根据公式

$\mathrm{RMSPEPR}=\sqrt{\frac{1}{n}\underset{t=1}{\overset{n}{\mathrm{\Σ}}}{[{\hat{x}}_1^{\left(0\right)}\left(t\right)-x_1^{\left(0\right)}\left(t\right)]}^2/{\left[x_1^{\left(0\right)}\left(t\right)\right]}^2}\×100\%---\left(10\right)$

$\mathrm{RMSPEPO}=\sqrt{\frac{1}{m}\underset{t=n+1}{\overset{n+m}{\mathrm{\Σ}}}{[{\hat{x}}_1^{\left(0\right)}\left(t\right)-x_1^{\left(0\right)}\left(t\right)]}^2/{\left[x_1^{\left(0\right)}\left(t\right)\right]}^2}\×100\%---\left(11\right)$

计算S-GMC(1，N)的预测精度，

其中，RMSPEPR为由前部分样本，即t＝1，2，...，n组成的根平均平方误差百分比，而RMSPEPO为由后部分样本，即t＝n+1，n+2，...，n+m组成的根平均平方误差百分比；

步骤六、判断条件是否满足，其中ε₁为故障阈值上界，ε₂为故障阈值下界，

若满足，则说明待测设备无故障；否则，判定待测设备有故障。

步骤二中κ的推导获取过程为：

存在公式(5)，

令 $\mathrm{\ψ}\left(t\right)=\frac{{\mathrm{dx}}_1^{\left(1\right)}\left(t\right)}{\mathrm{dt}}=-\mathrm{\α}{x}_1^{\left(1\right)}\left(t\right)+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(t\right)+{\mathrm{\β}}_N---\left(12\right)$

则有

$\underset{t}{\overset{t+2}{\∫}}\frac{{\mathrm{dx}}_1^{\left(1\right)}\left(t\right)}{\mathrm{dt}}\mathrm{dt}=\underset{t}{\overset{t+2}{\∫}}\mathrm{\ψ}\left(t\right)\mathrm{dt}---\left(13\right)$

于是

$x_1^{\left(1\right)}(t+2)-x_1^{\left(1\right)}\left(t\right)=x_1^{\left(0\right)}(t+1)+x_1^{\left(0\right)}(t+2)=\underset{t}{\overset{t+2}{\∫}}\mathrm{\ψ}\left(t\right)\mathrm{dt}---\left(14\right)$

采用代数精度为3的Simpson公式来代替梯形公式，并得

$x_1^{\left(0\right)}(t+1)+x_1^{\left(0\right)}(t+2)=\frac{t+2-t}{6}[\mathrm{\ψ}\left(t\right)+4\mathrm{\ψ}(t+1)+\mathrm{\ψ}(t+2)]$ (15)

$=\frac{1}{3}[\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i(x_{i+1}^{\left(1\right)}\left(t\right)+4x_{i+1}^{\left(1\right)}(t+1)+x_{i+1}^{\left(1\right)}(t+2))-\mathrm{\α}(x_1^{\left(1\right)}\left(t\right)+{4x}_1^{\left(1\right)}(t+1)+x_1^{\left(1\right)}(t+2))+6{\mathrm{\β}}_N]$

简化为Bκ＝Y，其中κ＝[α β₁ β₂ ... β_N]^T，Y和B分别由公式(6)和公式(7)求得

因此，由最小二乘法可以推导得出κ＝(B^TB)^-1B^TY。

在原始GM(1，N)(Original GM(1，N)，以下简称OGM(1，N))和基于卷积积分(Convolution Integral)的GM(1，N)(以下简称GMC(1，N))中，对公式(14)的离散化采用的是梯形公式，其代数精度仅为1，本发明采用代数精度为3的Simpson公式来代替梯形公式，提高了灰色系数的精度。

步骤三中公式(8)的推导获取过程为：

求解公式(5)所表示的S-GMC(1，N)，得

${\hat{x}}_1^{\left(1\right)}\left(t\right)=x_1^{\left(1\right)}{\left(1\right)e}^{-\mathrm{\α}(t-1)}+\underset{1}{\overset{t}{\∫}}{e}^{-\mathrm{\α}(t-\mathrm{\τ})}(\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(\mathrm{\τ}\right)+{\mathrm{\β}}_N)\mathrm{d\τ}$

(16)

$=x_1^{\left(1\right)}\left(1\right)e^{-\mathrm{\α}(t-1)}+\frac{{\mathrm{\β}}_N}{\mathrm{\α}}(1-e^{-\mathrm{\α}(t-1)})+\underset{1}{\overset{t}{\∫}}{e}^{-\mathrm{\α}(t-\mathrm{\τ})}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(\mathrm{\τ}\right)\mathrm{d\τ}$

定义利用复化Simpson公式可以简化I(t)，则

$I\left(1\right)=\underset{1}{\overset{1}{\∫}}g\left(\mathrm{\τ}\right)\mathrm{d\τ}=\underset{1}{\overset{1}{\∫}}{e}^{-\mathrm{\α}(1-\mathrm{\τ})}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(\mathrm{\τ}\right)\mathrm{d\τ}=0---\left(17\right)$

$I\left(2\right)=\underset{1}{\overset{2}{\∫}}g\left(\mathrm{\τ}\right)\mathrm{d\τ}=\frac{1}{2}[g\left(1\right)+g\left(2\right)]$

$=\frac{1}{2}[e^{-\mathrm{\α}(2-1)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(1\right)+e^{-\mathrm{\α}(2-2)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(2\right)]---\left(18\right)$

$=\frac{1}{2}[e^{-\mathrm{\α}}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(1\right)+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(2\right)]$

$I\left(3\right)=\underset{1}{\overset{3}{\text{\∫}}}g\left(\mathrm{\τ}\right)\mathrm{d\τ}=\frac{3-1}{6}[g\left(1\right)+4g\left(2\right)+g\left(3\right)]$

$=\frac{1}{3}[e^{-\mathrm{\α}(3-1)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(1\right)+4e^{-\mathrm{\α}(3-2)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(2\right)+e^{-\mathrm{\α}(3-3)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(3\right)]---\left(19\right)$

$=\frac{1}{3}[e^{-2\mathrm{\α}}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(1\right)+4e^{-\mathrm{\α}}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(2\right)+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(3\right)]$

$I\left(4\right)=\underset{1}{\overset{4}{\∫}}g\left(\mathrm{\τ}\right)\mathrm{d\τ}=\underset{1}{\overset{2}{\∫}}g\left(\mathrm{\τ}\right)\mathrm{d\τ}+\underset{2}{\overset{4}{\∫}}g\left(\mathrm{\τ}\right)\mathrm{d\τ}$

$=\frac{1}{2}[g\left(1\right)+g\left(2\right)]+\frac{4-2}{6}[g\left(2\right)+4g\left(3\right)+g\left(4\right)]$

$=\frac{1}{2}[e^{-\mathrm{\α}(4-1)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(1\right)+e^{-\mathrm{\α}(4-2)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(2\right)]---\left(20\right)$

$+\frac{1}{3}[e^{-\mathrm{\α}(4-2)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(2\right)+4e^{-\mathrm{\α}(4-3)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(3\right)+e^{-\mathrm{\α}(4-4)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(4\right)]$

$=\frac{1}{2}{e}^{-3\mathrm{\α}}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(1\right)+\frac{5}{6}{e}^{-2\mathrm{\α}}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(2\right)+\frac{1}{3}[4e^{-\mathrm{\α}}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(3\right)+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(4\right)]$

对于t≥5，可以将t分成两类：奇数和偶数，设t≥5，t＝2l+1， $l=\frac{t-1}{2},$ $h=\frac{t-1}{l}=\frac{t-1}{t-1/2}=2,$ 那么有

$I\left(t\right)=\underset{1}{\overset{t}{\∫}}g\left(\mathrm{\τ}\right)\mathrm{d\τ}=\frac{h}{6}\underset{j=0}{\overset{l-1}{\mathrm{\Σ}}}[g\left({\mathrm{\τ}}_j\right)+4g\left({\mathrm{\τ}}_{j+1/2}\right)+g\left({\mathrm{\τ}}_{j+1}\right)]$

$=\frac{1}{3}\underset{j=0}{\overset{l-1}{\mathrm{\Σ}}}[g(2j+1)+4g(2j+2)+g(2j+3)]$

$=\frac{1}{3}[g\left(1\right)+\underset{j=1}{\overset{l-1}{\mathrm{\Σ}}}g(2j+1)+4\underset{j=0}{\overset{l-1}{\mathrm{\Σ}}}g(2j+2)+\underset{j=1}{\overset{l}{\mathrm{\Σ}}}g(2j+3)]$

$=\frac{1}{3}[g\left(1\right)+2\underset{j=1}{\overset{l-1}{\mathrm{\Σ}}}g(2j+1)+4\underset{j=0}{\overset{l-1}{\mathrm{\Σ}}}g(2j+2)+g(2l+1)]$

$=\frac{1}{3}[g\left(1\right)+2\underset{j=1}{\overset{t-3/2}{\mathrm{\Σ}}}g(2j+1)+4\underset{j=0}{\overset{t-3/2}{\mathrm{\Σ}}}g(2j+2)+g\left(t\right)]$

$=\frac{1}{3}[e^{-\mathrm{\α}(t-1)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(1\right)+2\underset{j=1}{\overset{t-3/2}{\mathrm{\Σ}}}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{e}^{-\mathrm{\α}[t-(2j+1)]}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}(2j+1)]$

$+4\underset{j=0}{\overset{t-3/2}{\mathrm{\Σ}}}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{e}^{-\mathrm{\α}[t-(2j+2)]}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}(2j+2)+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(t\right)]---\left(21\right)$

其中，τ_j＝2j+1，且τ_j+1/2为τ_j和τ_j+1的中值。

同理，当t为偶数时，假设t≥6，t＝2l+2， 则

$I\left(t\right)=\underset{1}{\overset{t}{\∫}}g\left(\mathrm{\τ}\right)\mathrm{d\τ}=\underset{1}{\overset{2}{\∫}}g\left(\mathrm{\τ}\right)\mathrm{d\τ}+\underset{2}{\overset{t}{\∫}}g\left(\mathrm{\τ}\right)\mathrm{d\τ}$

$=\frac{1}{2}[g\left(1\right)+g\left(2\right)]+\frac{h}{6}\underset{j=0}{\overset{l-1}{\mathrm{\Σ}}}[g\left({\widetilde{\mathrm{\τ}}}_j\right)+4g\left({\widetilde{\mathrm{\τ}}}_{j+1/2}\right)+g\left({\widetilde{\mathrm{\τ}}}_{j+1}\right)]$

$=\frac{1}{2}[g\left(1\right)+g\left(2\right)]+\frac{1}{3}\underset{j=0}{\overset{l-1}{\mathrm{\Σ}}}[g(2j+2)+4g(2j+3)+g(2i+4)]$

$=\frac{1}{2}[g\left(1\right)+g\left(2\right)]+\frac{1}{3}[g\left(2\right)+\underset{j=1}{\overset{l-1}{\mathrm{\Σ}}}g(2j+2)+4\underset{j=0}{\overset{l-1}{\mathrm{\Σ}}}g(2j+3)+\underset{j=0}{\overset{l-1}{\mathrm{\Σ}}}g(2j+4)]$

$=\frac{1}{2}g\left(1\right)+\frac{5}{6}g\left(2\right)+\frac{1}{3}[2\underset{j=1}{\overset{l-1}{\mathrm{\Σ}}}g(2j+2)+4\underset{j=0}{\overset{l-1}{\mathrm{\Σ}}}g(2j+3)+g(2l+2)]$

$=\frac{1}{2}g\left(1\right)+\frac{5}{6}g\left(2\right)+\frac{1}{3}[2\underset{j=1}{\overset{t-4/2}{\mathrm{\Σ}}}g(2j+2)+4\underset{j=0}{\overset{t-4/2}{\mathrm{\Σ}}}g(2j+3)+g\left(t\right)]$

$=\frac{1}{2}{e}^{-\mathrm{\α}(t-1)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(1\right)+\frac{5}{6}{e}^{-\mathrm{\α}(t-2)}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(2\right)+\frac{1}{3}[2\underset{j=1}{\overset{t-4/2}{\mathrm{\Σ}}}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{e}^{-\mathrm{\α}[t-(2j+2)]}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}(2j+2)$

$+4\underset{j=0}{\overset{t-4/2}{\mathrm{\Σ}}}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{e}^{-\mathrm{\α}[t-(2j+3)]}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}(2j+3)+\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(t\right)]---\left(22\right)$

其中， ${\widetilde{\mathrm{\τ}}}_j=2j+2,$ 且为和的中值，

将(17)-(22)代入(16)，即可推导得到公式(8)。

具体实施方式二：本实施方式给出一个具体实施例：

下面对表1所示某机械设备的运行指标测量值(其中，前5个值用于建模，后5个值用于求预测精度)，建立基于复化Simpson公式改进多变量灰色模型(以下简称S-GMC(1，N)，本实施例中取S-GMC(1，N)，N＝2)，来阐述本发明的具体实施方式。

执行步骤一、假设输入的故障原始特征数据序列和相关数据序列分别为 $X_1^{\left(0\right)}=\{x_1^{\left(0\right)}\left(1\right),x_1^{\left(0\right)}\left(2\right),...,x_1^{\left(0\right)}\left(5\right)\}$ 及 $X_2^{\left(0\right)}=\{X_2^{\left(0\right)}\left(1\right),X_2^{\left(0\right)}\left(2\right),...,x_2^{\left(0\right)}\left(5\right),...,x_2^{\left(0\right)}\left(10\right)\},$ 具体数值如表1所示，其中由决定，待预测数据个数为5，k＝1，2，...，5且t＝1，2，...，10，为的原始值，用于计算预测精度。通过k＝1，2，...，5和t＝1，2，...，10，计算i＝1，2的一阶累加生成数。

执行步骤二、以为初值，建立如下S-GMC(1，2)

$\frac{{\mathrm{dx}}_1^{\left(1\right)}\left(1\right)}{\mathrm{dt}}+\mathrm{\α}{x}_1^{\left(1\right)}\left(t\right)=\underset{i=1}{\overset{2-1}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_{i+1}^{\left(1\right)}\left(t\right)+{\mathrm{\β}}_2,$ t＝1，2，...，10

由κ＝[α β₁ β₂]^T＝(B^TB)^-1B^TY，计算S-GMC(1，2)的灰色系数，其中

$Y={\left[\begin{array}{ccc}{x}_1^{\left(0\right)}\left(2\right)+x_1^{\left(0\right)}\left(3\right)& x_1^{\left(0\right)}\left(3\right)+x_1^{\left(0\right)}\left(4\right)& x_1^{\left(0\right)}\left(4\right)+x_1^{\left(0\right)}\left(5\right)\end{array}\right]}^T$

$B=\frac{1}{3}\left[\begin{array}{ccc}-(x_1^{\left(1\right)}\left(1\right)+4_{x_1}^{\left(1\right)}\left(2\right)+x_1^{\left(1\right)}\left(3\right))& (x_2^{\left(1\right)}\left(1\right)+{4x}_2^{\left(1\right)}\left(2\right)+x_2^{\left(1\right)}\left(3\right))& 6\\ -(x_1^{\left(1\right)}\left(2\right)+4x_1^{\left(1\right)}\left(3\right)+x_1^{\left(1\right)}\left(4\right))& (x_2^{\left(1\right)}\left(2\right)+4x_2^{\left(1\right)}\left(3\right)+x_2^{\left(1\right)}\left(4\right))& 6\\ -(x_1^{\left(1\right)}\left(3\right)+{4x}_1^{\left(1\right)}\left(4\right)+x_1^{\left(1\right)}\left(5\right))& (x_2^{\left(1\right)}\left(3\right)+{4x}_2^{\left(1\right)}\left(4\right)+x_2^{\left(1\right)}\left(5\right))& 6\end{array}\right]$

其实，将上述Y和B公式代入κ＝(B^TB)^-1B^TY，得

κ＝[α β₁ β₂]^T＝[-0.5965 -0.7915 4.0849]^T。

执行步骤三、利用步骤二求得的灰色系数κ＝[-0.5965 -0.7915 4.0849]^T，就可以通过下式对t＝1，2，...，10进行预测

其中，t＝1，2，...，10为t＝1，2，...，10的预测值，图2所示为本发明方法与基于卷积积分的多变量模型(Convolution Integral based GM(1，N)，以下简称GMC(1，N)，本实施例中取GMC(1，2)，N＝2)求得的t＝1，2，...，10的预测值及原始值对比图。可以看出，本发明得到的t＝1，2，...，10与原始值较为贴近，预测结果较为理想。

执行步骤四、由步骤三得到的t＝1，2，...，10，可求得

${\hat{x}}_1^{\left(0\right)}\left(t\right)={\hat{x}}_1^{\left(1\right)}\left(t\right)-{\hat{x}}_1^{\left(1\right)}(t-1),$ t＝2，3，...，10

其中， t＝1，2，...，10为t＝1，2，...，10的预测值，图3显示了本发明方法与GMC(1，2)求得的t＝1，2，...，10的预测值及原始值对比图。与图2类似，可以发现本发明得到的t＝1，2，...，10与原始值较为贴近，预测结果较为理想。

执行步骤五、由步骤一至四得到了t＝1，2，...，10，则可以通过下式计算S-GMC(1，2)的预测精度

$\mathrm{RMSPEPR}=\sqrt{\frac{1}{5}\underset{t=1}{\overset{5}{\mathrm{\Σ}}}{[{\hat{x}}_1^{\left(0\right)}\left(t\right)-x_1^{\left(0\right)}\left(t\right)]}^2/{\left[x_1^{\left(0\right)}\left(t\right)\right]}^2}\×100\%$

$\mathrm{RMSPEPO}=\sqrt{\frac{1}{5}\underset{t=6}{\overset{10}{\mathrm{\Σ}}}{[{\hat{x}}_1^{\left(0\right)}\left(t\right)-x_1^{\left(0\right)}\left(t\right)]}^2/{\left[x_1^{\left(0\right)}\left(t\right)\right]}^2}\×100\%$

其中，RMSPEPR为由前部分样本(即t＝1，2，...，5)组成的根平均平方误差百分比(RootMean Squared Percentage Error of PRiori-sample period)，而RMSPEPO为由后部分样本(即t＝6，7，...，10)组成的根平均平方误差百分比(Root Mean Squared Percentage Error ofPOst-sample period)。图4所示为本发明方法与GMC(1，2)的预测精度(包括RMSPEPR和RMSPEPO)对比图。从图4可以得出，本发明方法得到的RMSPEPR和RMSPEPO比GMC(1，2)分别降低了3.45％和20.98％，预测精度大大提高。

步骤六、判断条件t＝6，7，...，10是否满足，其中ε₁＝0，ε₂＝45分别为故障上界和下界。从图3知，GMC(1，2)在序列6之后开始严重偏离原始值，并在序列9至序列10判断出故障，产生了虚警，而S-GMC(1，2)与原始值吻合较好，能准确地对设备故障进行预测。

总之，图2至图4充分说明了本发明S-GMC(1，N)方法能有效提高故障预测精度，准确反映原始数据变化趋势，是对原有方法的巧妙改进。

表1 实施例用到的数据序列

表1给出的是实施例用到的数据序列。

Claims (1)

1.基于复化Simpson公式改进多变量灰色模型的故障预测方法，其特征在于，该方法包括以下步骤：

步骤一、输入的待测设备的故障原始特征数据序列为

且

输入的待测设备的故障相关数据序列为

且

其中，m为待预测数据个数，

由决定，

通过公式

和

公式

计算的一阶累加生成数

步骤二、以为初值，建立如下多变量灰色模型S-GMC(1,N)：

由κ＝[α β₁ β₂ ... β_N]^T＝(B^TB)^-1B^TY，计算S-GMC(1,N)的灰色系数κ，其中

步骤三、利用步骤二求得的灰色系数，就可以利用下式对进行预测，获取的预测值

步骤四、由步骤三的并按下式获得的预测值

其中，初始化

步骤五、根据公式

计算S-GMC(1,N)的预测精度，

其中，RMSPEPR为由前部分样本，即t＝1,2,...,n组成的根平均平方误差百分比，而RMSPEPO为由后部分样本，即t＝n+1,n+2,...,n+m组成的根平均平方误差百分比；

步骤六、判断条件是否满足，其中ε₁为故障阈值上界，ε₂为故障阈值下 界，

若满足，则说明待测设备无故障；否则，判定待测设备有故障；

步骤二中κ的推导获取过程为：

存在公式(5)，

令

则有

于是

采用代数精度为3的Simpson公式来代替梯形公式，并得

简化为Bκ＝Y，其中κ＝[α β₁ β₂ ... β_N]^T，Y和B分别由公式(6)和公式(7)求得

因此，由最小二乘法可以推导得出κ＝(B^TB)^-1B^TY；

步骤三中公式(8)的推导获取过程为：

求解公式(5)所表示的S-GMC(1,N)，得

定义利用复化Simpson公式可以简化I(t)，则

对于t≥5，可以将t分成两类：奇数和偶数，设 那么有

其中，τ_j＝2j+1，且τ_j+1/2为τ_j和τ_j+1的中值 ，

同理，当t为偶数时，假设则

其中，且为和的中值，

将(17)-(22)代入(16)，即可推导得到公式(8)。