CN108446714A - 一种多工况下的非马尔科夫退化系统剩余寿命预测方法 - Google Patents

Abstract

本发明公开了一种多工况下的非马尔科夫退化系统剩余寿命预测方法，属于预测与健康管理中剩余寿命预测领域，包括以下步骤：读入系统状态监测数据，初始化退化模型参数；利用野生二值分割方法检测变点的发生时刻，并在此基础上通过小波离散二阶导数算法估计模型的赫斯特指数；基于分层极大似然方法进一步估计模型的群体参数和个体参数；利用一种分形布朗运动的弱收敛准则，给出系统剩余寿命的概率密度函数的解析解；最后，结合施瓦茨信息准则和均方误差分析模型对数据的拟合程度以及寿命预测的准确性。对比仿真表明，本方法较传统方法具有更高的精度。

Description

一种多工况下的非马尔科夫退化系统剩余寿命预测方法

技术领域

本发明属于预测与健康管理中剩余寿命预测领域，具体涉及一种多工况下的非马尔科夫 退化系统剩余寿命预测方法。

背景技术

预测与健康管理是系统可靠性分析领域下的一项新兴技术，能够为长寿命操作设备的决 策过程提供大量有效信息。对于在役设备，剩余使用寿命是反映其健康状态的一个重要指标。 一般地，寿命预测的核心目标为在设备失效前准确地估计出当前时刻下剩余寿命的期望或概 率密度函数。

经过近些年的发展，可将现有的剩余寿命预测方法归纳为三大类，包括基于物理机理的 方法，基于经验的方法，以及数据驱动的方法。事实上，还有一系列融合不同类型思想的混 合方法。注意到数据驱动方法仅仅需要获取监测数据即可完成预测，而不需要其它的物理机 理或是专家知识，因此具有更为广泛的适用范围。通常来说，该方法的核心思想在于建立可 靠的随机过程模型，例如马尔科夫链、布朗运动、伽马过程以及逆高斯过程。

由于实际工业系统正趋向于大型化和复杂化的方向快速发展，对应的退化过程往往会历 经多个不同工况。以大型高炉为例，在生产过程中可能发生多种异常工况，包括悬料、管道 结瘤等等，会导致退化速率产生较为明显的变化。此外，高炉数据大多为非平稳、非马尔科 夫的，即退化过程中可能存在长程相关性之类的记忆效应。因此，在这种情况下，最为关键 的问题在于如何准确地预测含有多个工况的非马尔可夫退化系统的剩余寿命。

目前，针对多工况系统的寿命预测研究已经取得了一定的成果。现有方法主要可以分为 两种，一种是将变点或同阶段持续时间看作是随机变量，并假定其服从某种指定分布，例如， 高斯分布、对数正态分布以及高斯逆伽马分布等等，另一种则是直接进行变点检测，常用的 方法包括赤池准则、极大似然算法以及残差平方和等等。考虑到任意指定参数的先验分布的 可信度相对较低，第二类方法一般更为有效或是鲁棒。

发明内容

针对现有技术中存在的上述技术问题，本发明提出了一种多工况下的非马尔科夫退化系 统剩余寿命预测方法，设计合理，克服了现有技术的不足，具有良好的效果。

为了实现上述目的，本发明采用如下技术方案：

一种多工况下的非马尔科夫退化系统剩余寿命预测方法，包括如下步骤：

步骤1：读入M组状态监测数据，记为X₁,X₂,…,X_m,…,X_M；其中，

步骤2：根据多阶段退化模型，初始化参数集

其中，x_m(t)为第m个设备在t时刻的退化数据，j＝1,2,…,D，D为变点个数，为第 m个设备的第j个变点，N为监测时刻的总数，δ为采样间隔，表示跳变，为漂移系数，σ^(j+1)为扩散系数，B_H(t)为标准分数布朗运动，对于第一阶段 有

步骤3：利用野生二值分割方法检测变点时刻

步骤4：通过小波离散二阶导数算法估计模型的赫斯特指数H；

步骤5：利用分层极大似然方法分别估计模型的群体参数和个体参数

步骤6：基于分形布朗运动的弱收敛准则求取剩余寿命的概率密度函数f_m,k(Δ_m,k)；

步骤7：分别计算性能指标SIC_pop，SIC_ind和MSE，即

SIC_pop＝-2lΩ_pop|X+γ_popln N (15)；

其中，lΩ_pop|X为群体参数的对数似然函数，为第m 个设备的个体参数的对数似然函数，γ_pop和γ_ind分别为群体参数和个体参数的数量，为 Δ_m,k时刻剩余寿命的真值；

步骤8：根据三个性能指标值检验模型的预测效果，并最终输出概率密度函数f_m,kΔ_m,k。

优选地，在步骤3中，具体包括如下步骤：

步骤3.1：针对第m个设备，给出如下累积和统计量

其中，κ＝e-s+1表示区间间隔，e是终止时刻，s是起始时刻，τ的取值在e与s之间；

步骤3.2：任意选取P个子区间其中，p＝1,2,…,P，同时，[s₁,e₁]＝[s,e]， 则有如下检测逻辑

其中，为第p个区间内变点的估计值，δ为采样间隔，为步骤3.1中的累积 和统计量在第p个区间下各个监测时刻的取值；

将阈值设定为其中，MAD为x_m(t)的绝对中位差，C取1或1.3，则当时，认定为系统的一个变点，然后分别在和上重复步骤3.2的变点检测过程，即在拆分后的子区间内利用公式(3)进一步估计变点的位置，变点的最终取值为P次估计的均值。

优选地，在步骤4中，具体包括如下步骤：

步骤4.1：构造如下形式的二次变分

其中，a_i表示由Symlet小波滤波器生成的第i个选定序列，q为各个序列的最大长度；

步骤4.2：对步骤4.1中的公式(4)进行线性回归求取赫斯特指数H；或者利用MATLAB 提供的wfbmesti工具函数求解此问题，输出向量的第二个元素即为赫斯特指数H的估计值。

优选地，在步骤5中，具体包括如下步骤：

步骤5.1：记群体参数集为根据前面给出的退化模型的定 义，进一步记 重写分段观测向量为显然，其中，表示的协方差矩阵；则Ω_pop的对数似然函数为：

步骤5.2：对求偏导并将导数置为零，可得

步骤5.3：将式(6)代回式(5)，对应的剖面对数似然函数为

其中，为第m个设备在第j阶段的长度；

步骤5.4：利用序列二次规划方法求解式(7)的最大化问题，即可得到两类群体参数和σ^(1:D+1)的估计值，分别记作和再将其代入式(6)即可求得剩余的群体参数的最终结果，记作

步骤5.5：记个体参数集为针对第m个设备，给出该 设备对应的个体参数的对数似然函数

步骤5.6：令m＝1,2,…,M，采用序列二次规划方法求解式(8)的最大化问题，即可求 得各个设备的个体参数的极大似然估计值，记作

优选地，在步骤6中，令表示独立同分布的随机变量，且有构造如下随机过程

其中，[·]为向下取整符号，当i>0，b_i＝(H-0.5)i^H1.5，当i≤0，b_i＝0；

当r→∞时，y_H,r(t)弱收敛于且有

则B_H(t)可近似为B(h(t))，因此，RUL_m,k的 概率密度函数可表示为

其中，

表示t_k时刻下的阈值，且有

结合多阶段模型的特性对RUL_m,k的概率密度函数，即式(11)进行修正：当令来计算Δ_m,k等价于的求和，i＝j+1,j+2,…,D1；假定相互独立，则f_m,k(Δ_m,k)的修正值可由依次 卷积进一步给出。

本发明所带来的有益技术效果：

考虑到绝大多数传统的基于马尔科夫过程的退化建模方法不能很好地解决上述问题，本 发明建立了一种改进的基于分数布朗运动的多阶段退化模型，同时考虑了单个系统的长程相 关性，多个系统之间的不确定性，以及工况切换时刻下可能存在的跳变；为了更好地辨识模 型参数，提出了一种四阶段参数估计方法，逐步估计变点、赫斯特指数、群体参数以及个体 参数；利用分数布朗运动的弱收敛准则，给出了任意时刻下系统剩余寿命的概率密度函数的 解析形式；相比传统方法，本发明提供的剩余寿命预测算法具有较高的预测精度和较小的不 确定性。

本发明方法能够处理含有多个未知工况(包括异常工况)的复杂退化系统，且剩余寿命 的预测精度有所提升；本发明方法相对传统的多阶段剩余寿命预测算法实用性更好，准确度 更高，主要原因在于所说模型将无记忆效应的马尔科夫过程升级为基于分数布朗运动的退化 过程，考虑了长期数据中实际存在的长程相关性，体现在分数布朗运动B_H(t)当中，具体来说： 当赫斯特指数0.5＜H＜1时，未来的退化状态趋向于沿着过去的方向继续向前发展；当 0＜H＜0.5时，未来的退化状态则会呈现出一个负相关的趋势；而当H＝0.5时，B_H(t)退化 为传统的标准布朗运动；当赫斯特指数0.5＜H＜1时，退化过程存在长期的记忆效应；对于 多个同批次设备，本发明能够同时实现各个系统的剩余寿命预测，同时通过结合不同退化信 息之间的相关关系，有效提升方法的准确性。

本发明可以应用到以下场合：

I高速列车制动系统，航天器推进装置等设备的性能退化过程；

II高炉炼铁，火力发电等工业生产过程。

附图说明

图1是本发明进行剩余寿命预测的流程图；

图2是实例一的两阶段退化轨迹示意图；

图3是实例一的变点检测结果示意图；

图4是实例一的所提模型的剩余寿命预测结果示意图；

图5是实例一的单阶段模型的剩余寿命预测结果示意图；

图6是实例一的布朗运动模型的剩余寿命预测结果示意图。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

一种多工况下的非马尔科夫退化系统剩余寿命预测方法，其流程如图1所示，包括如下 步骤：

步骤1：读入M组状态监测数据，记为X₁,X₂,…,X_m,…,X_M；其中，

步骤2：根据多阶段退化模型，初始化参数集

其中，x_m(t)为第m个设备在t时刻的退化数据，j＝1,2,…,D，D为变点个数，为第m个设备的第j个变点，N为监测时刻的总数，δ为采样间隔，表示跳变，为漂移系数，σ^(j+1)为扩散系数，B_H(t)为标准分数布朗运动，对于第一阶段有

步骤3：利用野生二值分割方法检测变点时刻具体包括如下步骤：

步骤3.1：针对第m个设备，给出如下累积和统计量

其中，κ＝e-s+1表示区间间隔，e是终止时刻，s是起始时刻，τ的取值在e与s之间；

步骤3.2：任意选取P个子区间其中，p＝1,2,…,P，同时，[s₁,e₁]＝[s,e]， 则有如下检测逻辑

其中，为第p个区间内变点的估计值，δ为采样间隔，为步骤3.1中的累积 和统计量在第p个区间下各个监测时刻的取值；

将阈值设定为其中，MAD为x_m(t)的绝对中位差，C取1或1.3，则当时，认定为系统的一个变点，然后分别在和上重复步骤3.2的变点检测过程，即在拆分后的子区间内利用公式(3)进一步估计变点的位置，变点的最终取值为P次估计的均值。

步骤4：通过小波离散二阶导数算法估计模型的赫斯特指数H；具体包括如下步骤：

步骤4.1：构造如下形式的二次变分

其中，a_i表示由Symlet小波滤波器生成的第i个选定序列，q为各个序列的最大长度；

步骤4.2：对步骤4.1中的公式(4)进行线性回归求取赫斯特指数H；或者利用MATLAB 提供的wfbmesti工具函数求解此问题，输出向量的第二个元素即为赫斯特指数H的估计值

步骤5：利用分层极大似然方法分别估计模型的群体参数和个体参数 具体包括如下步骤：

步骤5.1：记群体参数集为根据前面给出的退化模型的定 义，进一步记 重写分段观测 向量为显然，其中，表示的协方差矩阵；则Ω_pop的对数似然函数为：

步骤5.2：对求偏导并将导数置为零，可得

步骤5.3：将式(6)代回式(5)，对应的剖面对数似然函数为

其中，为第m个设备在第j阶段的长度；

步骤5.4：利用序列二次规划方法求解式(7)的最大化问题，即可得到两类群体参数和σ^(1:D+1)的估计值，分别记作和再将其代入式(6)即可求得剩余的群体参数的最终结果，记作

步骤5.5：记个体参数集为针对第m个设备，给出该 设备对应的个体参数的对数似然函数

步骤5.6：令m＝1,2,…,M，采用序列二次规划方法求解式(8)的最大化问题，即可求 得各个设备的个体参数的极大似然估计值，记作

步骤6：基于分形布朗运动的弱收敛准则求取剩余寿命的概率密度函数f_m,k(Δ_m,k)；

令表示独立同分布的随机变量，且有构造如下随机过 程

其中，[·]为向下取整符号，当i>0，b_i＝(H-0.5)i^H1.5，当i≤0，b_i＝0；

当r→∞时，y_H,r(t)弱收敛于且有

则B_H(t)可近似为B(h(t))，因此，RUL_m,k的 概率密度函数可表示为

其中，

表示t_k时刻下的阈值，且有

结合多阶段模型的特性对RUL_m,k的概率密度函数，即式(11)进行修正：当令来计算Δ_m,k等价于的求 和，i＝j+1,j+2,…,D1；假定相互独立，则f_m,k(Δ_m,k)的修正值可由依次卷积进一步给出。

步骤7：分别计算性能指标SIC_pop，SIC_ind和MSE；

SIC_pop＝-2lΩ_pop|X+γ_popln N (15)；

其中，lΩ_pop|X为群体参数的对数似然函数，为第m 个设备的个体参数的对数似然函数，γ_pop和_ind分别为群体参数和个体参数的数量，为 Δ_m,k时刻剩余寿命的真值；

步骤8：根据三个性能指标值检验模型的预测效果，并最终输出概率密度函数。

2、实例寿命预测结果

下面结合实例对该发明的具体实施步骤做进一步的具体说明。

下述实例的仿真环境如下：

机型：Intel Core i5-5200(CPU 2.20Ghz，8.00GBRAM)；

操作系统：Windows 7；

软件：Matlab R2013b。

a、实例一：通过一组数值仿真，给出实验步骤的详细说明以及剩余寿命预测结果：

1)考虑一类包含两种工况的退化系统，换言之，每个设备的退化过程存在一个变点，即 D＝1。根据模型(1)生成退化数据，其中，参数设置为M＝20，N＝250，δ＝0.1， σ⁽¹⁾＝0.6，σ⁽²⁾＝0.4，H＝0.6。个体参数由指定的 高斯分布随机生成，而变点均取在中点附近，系统寿命为25，具体的退化轨迹如图2所示；

2)初始化模型参数集

3)利用式(2)和式(3)检测变点时刻以第一个设备为例，5000次变点检测的结果如图3所示，可以看出，绝大多数次检测结果和真值相吻合，事实上，通过取中位数的操作，20个设备的变点均能成功地由野生二值分割方法检测得出；

4)利用式(4)或MATLAB提供的wfbmesti函数估计系统的赫斯特指数H，结果为与其真值十分接近，说明该方法能够准确地评估非马尔科夫过程当中存在的长 期记忆效应；

5)利用式(5)-式(7)估计模型的群体参数并在此基础上利用式 (8)估计个体参数

6)根据式(9)-式(14)求取剩余寿命的概率密度函数f_m,k(Δ_m,k)，为分析算法性能，同时选用所提模型(1)，传统的单阶段模型，以及布朗运动模型进行对比仿真，对应的剩余寿命预测结果如图4-6所示，可以看出，单阶段模型的预测结果不够理想，特别是在第一阶段，误差相对较大，这是由于仿真数据在变点处存在跳变，而单阶段模型无法实现变点的检测；粗略地看，布朗运动模型和所提模型均能得到不错的预测结果，即各个监测时刻下剩余寿命估计的期望值和其真值基本吻合，事实上，所提模型的概率密度函数具有更小的方差，说明预测的不确定性更低；

7)利用式(15)-式(17)计算三个不同模型的SIC_po_p，SIC_ind，以及MSE，结果见表1，可以看出，所提模型(1)的各个性能指标均为最小，说明该模型对仿真数据的拟合能力最强， 同时具有最高的预测精度。

表1.模型性能对比

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的 技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护 范围。

Claims (5)

1.一种多工况下的非马尔科夫退化系统剩余寿命预测方法，其特征在于：包括如下步骤：

步骤1：读入M组状态监测数据，记为X₁,X₂,…,X_m,…,X_M；其中，

步骤2：根据多阶段退化模型，初始化参数集

其中，x_m(t)为第m个设备在t时刻的退化数据，j＝1,2,…,D，D为变点个数，为第m个设备的第j个变点，N为监测时刻的总数，δ为采样间隔，表示跳变，为漂移系数，σ^(j+1)为扩散系数，B_H(t)为标准分数布朗运动，对于第一阶段有

步骤3：利用野生二值分割方法检测变点时刻

步骤4：通过小波离散二阶导数算法估计模型的赫斯特指数H；

步骤5：利用分层极大似然方法分别估计模型的群体参数和个体参数

步骤6：基于分形布朗运动的弱收敛准则求取剩余寿命的概率密度函数f_m,k(Δ_m,k)；

步骤7：分别计算性能指标SIC_pop，SIC_ind和MSE，即

SIC_pop＝-2l(Ω_pop|X)+Υ_popln N (15)；

其中，l(Ω_pop|X)为群体参数的对数似然函数，为第m个设备的个体参数的对数似然函数，Υ_pop和Υ_ind分别为群体参数和个体参数的数量，为Δ_m,k时刻剩余寿命的真值；

步骤8：根据三个性能指标值检验模型的预测效果，并最终输出概率密度函数f_m,k(Δ_m,k)。

2.根据权利要求1所述的多工况下的非马尔科夫退化系统剩余寿命预测方法，其特征在于：在步骤3中，具体包括如下步骤：

步骤3.1：针对第m个设备，给出如下累积和统计量

其中，κ＝e-s+1表示区间间隔，e是终止时刻，s是起始时刻，τ的取值在e与s之间；

步骤3.2：任意选取P个子区间其中，p＝1,2,…,P，同时，[s₁,e₁]＝[s,e]，则有如下检测逻辑

其中，为第p个区间内变点的估计值，δ为采样间隔，为步骤3.1中的累积和统计量在第p个区间下各个监测时刻的取值；

将阈值设定为其中，MAD为x_m(t)的绝对中位差，C取1或1.3，则当时，认定为系统的一个变点，然后分别在和上重复步骤3.2的变点检测过程，即在拆分后的子区间内利用公式(3)进一步估计变点的位置，变点的最终取值为P次估计的均值。

3.根据权利要求1所述的多工况下的非马尔科夫退化系统剩余寿命预测方法，其特征在于：在步骤4中，具体包括如下步骤：

步骤4.1：构造如下形式的二次变分

其中，a_i表示由Symlet小波滤波器生成的第i个选定序列，q为各个序列的最大长度；

步骤4.2：对步骤4.1中的公式(4)进行线性回归求取赫斯特指数H；或者利用MATLAB提供的wfbmesti工具函数求解此问题，输出向量的第二个元素即为赫斯特指数H的估计值。

4.根据权利要求1所述的多工况下的非马尔科夫退化系统剩余寿命预测方法，其特征在于：在步骤5中，具体包括如下步骤：

步骤5.1：记群体参数集为根据前面给出的退化模型的定义，进一步记 重写分段观测向量为显然，其中，表示的协方差矩阵；则Ω_pop的对数似然函数为：

步骤5.2：对求偏导并将导数置为零，可得

步骤5.3：将式(6)代回式(5)，对应的剖面对数似然函数为

其中，为第m个设备在第j阶段的长度；

步骤5.4：利用序列二次规划方法求解式(7)的最大化问题，即可得到两类群体参数和σ^(1:D+1)的估计值，分别记作和再将其代入式(6)即可求得剩余的群体参数的最终结果，记作

步骤5.5：记个体参数集为针对第m个设备，给出该设备对应的个体参数的对数似然函数

步骤5.6：令m＝1,2,…,M，采用序列二次规划方法求解式(8)的最大化问题，即可求得各个设备的个体参数的极大似然估计值，记作

5.根据权利要求1所述的多工况下的非马尔科夫退化系统剩余寿命预测方法，其特征在于：在步骤6中，令表示独立同分布的随机变量，且有E(ρ₀)＝0，构造如下随机过程

其中，[·]为向下取整符号，当i>0，b_i＝(H-0.5)i^H-1.5，当i≤0，b_i＝0；

当r→∞时，y_H,r(t)弱收敛于且有

则B_H(t)可近似为B(h(t))，因此，RUL_m,k的概率密度函数可表示为

其中，

表示t_k时刻下的阈值，且有

结合多阶段模型的特性对RUL_m,k的概率密度函数，即式(11)进行修正：当令来计算Δ_m,k等价于的求和，i＝j+1,j+2,…,D+1；假定相互独立，则f_m,k(Δ_m,k)的修正值可由依次卷积进一步给出。