CN112883550A - 一种考虑多重不确定性的退化设备剩余寿命预测方法 - Google Patents

Abstract

本发明针对长寿命和高可靠性设备缺乏充足退化数据的情况，提出一种考虑多重不确定的退化设备剩余寿命预测方法，包括基于非线性扩散过程的步进加速退化模型，该模型的优点在于仅需要较小的样本量和较短的测试时间。对于退化模型的固有特性、个体差异性、测量设备性能和退化过程中因人为偏差造成的多重不确定性，该模型考虑了时变不确定性、个体差异性以及性能退化和协变量的测量不确定性。为估计退化设备的剩余寿命，本发明推导得到了首达时间意义下非线性扩散过程穿过预定阈值的解析近似解。结合极大似然估计(MLE)方法和仿真外推(SIMEX)方法，得到一种MLE‑SIMEX方法，用于估计模型中的未知参数。通过仿真和实际案例证明了本发明所提出模型的有效性。结果表明，该方法具有更高的剩余寿命估计精度，具有一定的工程实用价值。

Description

一种考虑多重不确定性的退化设备剩余寿命预测方法

技术领域

本发明属于可靠性工程技术领域，具体涉及一种考虑多重不确定性的退化设备剩余寿命预测方法。

背景技术

随着现状设计水平和制造工艺的不断提高，人们对高质量产品的需求也与日俱增，出现了越来越多使用寿命长、可靠性高的设备，在航空航天和军事领域显得尤为突出。由于测试时间长、测试成本高，传统的剩余寿命预测方法难以有效的应用于这类设备。为获得更多的退化数据，工程师通过提高测试环境的严苛性来加速设备的退化，例如高温高压、随机振动等。因此，加速退化试验已成为收集产品退化数据的有效方法。

在获得退化数据后，可以利用随机过程模型对退化数据进行建模分析，以估计设备的剩余寿命。对于复杂的系统或设备而言，非线性和不确定性是影响退化过程建模精度的两个重要因素。在非线性加速退化过程模型中存在多重不确定性，包括时变不确定性，个体差异性以及性能退化和协变量中的测量不确定性。通过建立基于随机过程的退化模型解决设备退化的时变不确定性问题；考虑退化模型中漂移系数的随机效应以描述退化设备的个体差异性；此外，由于测量设备性能的限制以及测试数据读数的人为偏差，在大多数协变量和退化性能的测量过程中都会存在测量误差。目前，已有方法仅考虑了多重不确定性中的部分，缺乏考虑多重不确定性的非线性加速退化建模与剩余寿命预测方法的相关研究。

发明内容

本发明的目的针对长寿命和高可靠性设备缺乏充足退化数据的情况，提出一种基于非线性扩散过程的步进加速退化模型以估计设备的剩余寿命。对于退化模型的固有特性、个体差异性、测量设备性能和退化过程中因人为偏差造成的多重不确定性，该模型考虑了时变不确定性、个体差异性以及性能退化和协变量的测量不确定性。

本发明采用的技术方案是：

一种考虑多重不确定性的退化设备剩余寿命预测方法，包括以下步骤：

步骤1、针对退化设备，基于扩散过程建立一般的非线性加速退化过程模型，并利用标准布朗运动描述随机退化过程随时间变化的固有不确定性；

步骤2、利用阿伦尼斯模型来描述退化模型漂移系数与加速应力间的关系，同时考虑退化设备的个体差异性和加速应力作为协变量的测量不确定性，建立考虑多重不确定性的加速退化模型；

步骤3、考虑性能退化过程中的测量不确定性，基于首达时间的概念，推导退化设备剩余寿命分布的表达式；

步骤4、利用MLE-SIMEX方法估计虑多重不确定性的加速退化模型中的未知参数，将其代入退化设备剩余寿命分布的表达式，从而实现退化设备的寿命预测。

优选的，在步骤1中，所述非线性加速退化过程模型建立如下：

设X(t)表示样品在t时刻的性能退化量，则基于扩散过程的退化过程{X(t),t≥0}可以表示为：

其中，μ(t；θ)为漂移系数，是时间t的非线性函数，用以表示模型的非线性特征，其中参数向量θ＝(a,b)为未知参数，a表征不同设备间的差异性，b表征同类设备的共性特征；σ_B称为扩散系数，B(t)为标准Brownian运动，且B(t):N(0,t)；

对于式(1)，在首达时间概念下，设备的寿命T可定义为：

T＝inf{t:X(t)≥ω|X(0)＜ω} (2)

当不考虑参数a和σ_B的随机性时，设备寿命的概率密度函数可以近似为：

其中，

ω为失效阈值。

优选的，在步骤2中，考虑多重不确定性的加速退化模型的建立如下:

假设存在K步的步进加速退化试验数据，试验以最低的应力水平S₁开始，初始时间t＝0；每个设备都将在预定的时间经历预定的应力增加过程，当t＝T_k(k＝1,2,L,K)时，相应应力为S_k，其表达如下:

漂移系数a与加速应力相关，利用阿伦尼斯模型来描述a与加速应力S的关系如下：

a_k＝c_aexp(-d_a/S_k) (5)

其中，a_k表示漂移系数a在不同应力水平S_k下的值，c_a是关于失效模式的常数，d_a是关于激活能和玻尔兹曼常数的常数；

对于由标准布朗运动{B(t),t≥0}驱动的的非线性退化模型，随机过程退化模型随时间变化的固有属性可描述时变不确定性；由于不同的工作条件，退化轨迹将表现出不同的退化速率，可利用漂移系数可作为描述个体差异性的随机参数，漂移系数确定了设备在不同应力下的退化轨迹，考虑到a_k的随机效应，式(5)中，c_a可作为随机参数且假设

μ_c和

分别表示c_a的均值和方差；因此，能够得到

μ_ak和

分别表示a_k的均值和方差，即：

实际工程中，为描述测量不确定性的影响，令{Y(t),t≥0}表示设备的测量过程，用来描述隐藏退化状态与测量值之间的关系，具体如下：

Y(t)＝X(t)+ε (7)

其中，ε是随机测量误差，假设

表示ε的方差；

通常加速退化试验在预定的应力水平下进行，由于测量误差的存在，使得设备实际上在未知的加速应力水平下发生退化，加速应力的测量误差存在于协变量的测量过程中，可将协变量的测量误差用随机变量

表示，且

表示

的方差；因此，实际加速应力水平可表示为

即：

综上所述，考虑多重不确定性的加速退化模型可表示为：

优选的，在步骤3中，退化设备剩余寿命分布表达式推导如下：

在首达时间下，考虑多重不确定性的设备寿命T_e可定义为：

T_e＝inf{t:Y(t)≥ω|Y(0)＜ω} (10)

相应的，设备寿命的概率密度函数表示为

首先，在不考虑漂移系数随机效应的情况下推导

对于退化过程{Y(t),t≥0}，其寿命定义如式(10)所示，则寿命T_e的

可表示为

其中，B＝ω-a·h(t,θ)，C＝σ²t；

为了估计单个设备的剩余寿命，需要得到该设备在正常工作应力下的退化数据，Y(t_k)表示设备在t_k时刻的测量值，根据首达时间的概念，定义t_k时刻的剩余寿命L_k为

L_k＝inf{l_k>0:Y(l_k+t_k)≥ω} (12)

已知t_k时刻的退化状态测量值为Y(t_k)，则剩余寿命的概率密度函数可表示为

其中，

Ω(l_k)＝h(l_k+t_k,θ)-h(t_k,θ)；

利用漂移系数a的随机性表示退化设备的个体差异性，并且a服从式(6)的分布形式；

对于退化过程{Y(t),t≥0}，已知

并且t_k时刻的退化状态测量值为Y(t_k)，则t_k时刻剩余寿命的概率密度函数可表示为

其中，

L＝Ω(l_k)-l_kμ(l_k+t_k；θ)，M＝Ω(l_k)＝h(l_k+t_k,θ)-h(t_k,θ)，

优选的，在步骤4中，利用MLE-SIMEX方法估计模型参数如下:

(1)仿真步骤

给定一个正整数G和一个序列V＝{v₁,v₂,L,v_M}，其中，v₁＝0，v_M是给定的正数，并且v₂,L,v_M-1服从均匀分布U(0,v_M)；对于试验应力S_k，k＝1,2,L,K，令

g＝1,2,L,G，表示K×G个服从标准正态分布的样本；因此，可仿真得到第k个协变量的应力水平如下

其中，m＝1,2,L,M；

(2)极大似然估计步骤

假设试验应力

有N个试验样本，退化状态测量值y(t_ijk)表示第j个样本在第k个应力下的第i个测量数据，对应的测量时间表示为t_ijk，其中，i＝1,2,L,n_k，j＝1,2,L,N，k＝1,2,L,K，n_k表示在第k个应力下测量值的个数，实际的退化状态表示为x(t_ijk)；

因此，y(t_ijk)可表示为

y(t_ijk)＝c_ajΛ(t_ijk)+σB(t_ijk)+ε (16)

其中，

c_aj表示漂移参数a_j的系数，Λ(t_ijk)可表示为

为便于描述，令Y_j＝(y(t_1j1),y(t_2j2),L,y(t_nkjk))′表示第j个样本的退化测量矢量，T_j＝(Λ(t_1j1),Λ(t_2j2),L,Λ(t_nkjk))′，其中，(·)′表示矢量的转置，Y表示全部的加速退化数据，包括Y_j，j＝1,2,L,N；

根据式(16)可知，Y_j服从多维正态分布，其均值和协方差可表示为

其中，

I_l是l维的单位矩阵，

因此，关于未知参数Θ(g,m)＝(μ_c,σ_c,σ,σ_ε,θ,d_a)的对数似然函数可表示为

其中，

将对数似然函数(19)对μ_c求一阶偏导数，可得

令式(20)为0，可得

那么根据已估计的

参数σ_c，σ，σ_ε，θ和d_a的剖面似然函数可表示为

参数σ_c，σ，σ_ε，θ和d_a的极大似然估计值可利用多维搜索的方法，通过极大化剖面似然函数(22)得到，将

和

代入式(21)，即可得到

则极大似然估计值

即可得到；

(3)外推步骤

设备个体在工作应力下的退化数据已知，则可以根据贝叶斯方法更新退化模型的参数，从而得到该设备实时的剩余寿命预测值；

假设有κ个正常工作应力下的退化数据，表示为Y_κ＝[y₁,y₂,L,y_κ]，其中，y_κ表示设备在t_κ时刻的退化状态测量值；假设漂移系数a的先验估计分布服从正态分布，即

那么，联合后验分布p(a|Y_κ)可根据贝叶斯方法表示为

其中，p(Y_κ|a)为似然函数，p(a)为先验分布的PDF，Δy_g＝y_g-y_g-1，Δt_g＝t_g-t_g-1，

p(Y_κ|a)和p(a)均服从正态分布，那么后验分布p(a|Y_κ)与先验分布应属于同一分布族，即

其中，

将更新后的参数(25)代入(16)，即可得到设备个体剩余寿命的概率密度函数，实现剩余寿命的实时预测。

本发明的有益效果：本发明利用基于非线性扩散过程的步进加速退化试验模型描述设备的退化过程，该模型同时考虑了时变不确定性、个体差异性以及性能退化和协变量中的测量不确定性。该模型的优点在于仅需要较小的样本量和较短的测试时间。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为步进加速应力下陀螺仪退化曲线；

图2为工作应力下陀螺仪的退化曲线；

图3为模型参数更新曲线；

图4为陀螺仪剩余寿命概率密度函数。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明具体提供了一种考虑多重不确定性的退化设备剩余寿命预测方法，包括以下步骤：

步骤1、针对退化设备，基于扩散过程建立一般的非线性加速退化过程模型，并利用标准布朗运动描述随机退化过程随时间变化的固有不确定性；

非线性加速退化过程模型建立如下：

设X(t)表示样品在t时刻的性能退化量，则基于扩散过程的退化过程{X(t),t≥0}可以表示为：

其中，μ(t；θ)为漂移系数，是时间t的非线性函数，用以表示模型的非线性特征，其中参数向量θ＝(a,b)为未知参数，a表征不同设备间的差异性，b表征同类设备的共性特征；σ_B称为扩散系数，B(t)为标准Brownian运动，且B(t):N(0,t)。

对于不同函数形式的μ(t；θ)可以描述不同形式的非线性随机退化过程。显然，若μ(t；θ)＝μ，则式(1)转化为线性随机退化模型，即Wiener过程。因此，线性随机退化模型是本文所讨论的非线性随机退化模型的特例。不失一般性，本文主要考虑X(0)＝0的情况。

对于式(1)，在首达时间概念下，设备的寿命T可定义为：

T＝inf{t:X(t)≥ω|X(0)＜ω} (2)

当不考虑参数a和σ_B的随机性时，设备寿命的概率密度函数可以近似为：

其中，

ω为失效阈值。

步骤2、利用阿伦尼斯模型来描述退化模型漂移系数与加速应力间的关系，同时考虑退化设备的个体差异性和加速应力作为协变量的测量不确定性，建立考虑多重不确定性的加速退化模型

考虑多重不确定性的加速退化模型的建立如下:

在本发明中，假设存在K步的步进加速退化试验数据，试验以最低的应力水平S₁开始，初始时间t＝0；每个设备都将在预定的时间经历预定的应力增加过程，当t＝T_k(k＝1,2,L,K)时，相应应力为S_k，其表达如下:

漂移系数a与加速应力相关，利用阿伦尼斯模型来描述a与加速应力S的关系如下：

a_k＝c_aexp(-d_a/S_k) (5)

其中，a_k表示漂移系数a在不同应力水平S_k下的值，c_a是关于失效模式的常数，d_a是关于激活能和玻尔兹曼常数的常数；

对于由标准布朗运动{B(t),t≥0}驱动的的非线性退化模型，随机过程退化模型随时间变化的固有属性可描述时变不确定性；由于不同的工作条件，退化轨迹将表现出不同的退化速率，可利用漂移系数可作为描述个体差异性的随机参数，。漂移系数确定了设备在不同应力下的退化轨迹，考虑到a_k的随机效应，式(5)中，c_a可作为随机参数且假设

μ_c和

分别表示c_a的均值和方差；因此，能够得到

μ_ak和

分别表示a_k的均值和方差，即：

实际工程中，难以精确的测量设备的退化状态，不可避免的会受到噪声干扰和非理想仪器造成的测量不确定性影响，为描述测量不确定性的影响，令{Y(t),t≥0}表示设备的测量过程，用来描述隐藏退化状态与测量值之间的关系，具体如下：

Y(t)＝X(t)+ε (7)

其中，ε是随机测量误差，假设

表示ε的方差；

通常加速退化试验在预定的应力水平下进行，该应力在试验中测量并控制在规定的范围内，由于测量误差的存在，使得设备实际上在未知的加速应力水平下发生退化，加速应力的测量误差存在于协变量的测量过程中，可将协变量的测量误差用随机变量

表示，且

表示

的方差；因此，实际加速应力水平可表示为

即：

综上所述，考虑多重不确定性的加速退化模型可表示为：

步骤3、考虑性能退化过程中的测量不确定性，基于首达时间的概念，推导退化设备剩余寿命分布的表达式；

在首达时间下，考虑多重不确定性的设备寿命T_e可定义为：

T_e＝inf{t:Y(t)≥ω|Y(0)＜ω} (10)

相应的，设备寿命的概率密度函数表示为

首先，在不考虑漂移系数随机效应的情况下推导

引入如下引理。

引理1.存在Z:N(μ,σ²)，并且

可以得到

根据引理1，可以利用全概率公式得到首达时间意义下寿命T_e的概率密度函数(PDF)

具体结果如下表示：

对于退化过程{Y(t),t≥0}，其寿命定义如式(10)所示，则寿命T_e的

可表示为

其中，B＝ω-a·h(t,θ)，C＝σ²t。

证明：由式(7)可知，Y(t)＝X(t)+ε，其中

因此，寿命T_e即为{Y(t),t≥0}首达阈值ω-ε的时间。由于测量误差是随机的，则可通过全概率公式得到T_e的概率密度函数(PDF)。

令B＝ω-a·h(t,θ)，C＝σ²t，Z＝ε，根据引理1可得

证毕。

为了估计单个设备的剩余寿命，还需要得到该设备在正常工作应力下的退化数据，Y(t_k)表示设备在t_k时刻的测量值，根据首达时间的概念，定义t_k时刻的剩余寿命L_k为

L_k＝inf{l_k>0:Y(l_k+t_k)≥ω} (12)

已知t_k时刻的退化状态测量值为Y(t_k)，则剩余寿命的概率密度函数可表示为

其中，

Ω(l_k)＝h(l_k+t_k,θ)-h(t_k,θ)。

证明：已知t_k时刻的退化状态测量值为Y(t_k)，对于t≥t_k，退化过程可表示为

令l_k＝t-t_k，退化过程{Y(t),t≥t_k}可表示为

Y(l_k+t_k)-Y(t_k)＝a·(h(l_k+t_k,θ)-h(t_k,θ))+σB(l_k)

其中，

因此，t_k时刻设备的剩余寿命等价于首达时间意义下，退化过程{Z(l_k),l_k≥0}穿过阈值

的时间，其中，Z(l_k)＝Y(l_k+t_k)-Y(t_k)，并且Z(0)＝0。式可改写为

Z(l_k)＝a·(h(l_k+t_k,θ)-h(t_k,θ))+σB(l_k)

根据定理1的结论，即可得到剩余寿命PDF的表达式。

证毕。

利用漂移系数a的随机性表示退化设备的个体差异性，并且a服从式(6)的分布形式；

对于退化过程{Y(t),t≥0}，已知

并且t_k时刻的退化状态测量值为Y(t_k)，则t_k时刻剩余寿命的概率密度函数可表示为

其中，

证明：由全概率公式可得

令L＝Ω(l_k)-l_kμ(l_k+t_k；θ)，M＝Ω(l_k)＝h(l_k+t_k,θ)-h(t_k,θ)，

由引理2可得

证毕。

步骤4、利用MLE-SIMEX方法估计虑多重不确定性的加速退化模型中的未知参数，将其代入退化设备剩余寿命分布的表达式，从而实现退化设备的寿命预测

考虑到温度、湿度等试验应力通常由常用仪器进行测量，从独立样本中估计参数

是可行的，并且与其他未知参数不相关。因此，本发明假设

已知。本发明提出了一种改进的MLE-SIMEX方法来处理多变量非线性SSADT的参数估计问题。

利用MLE-SIMEX方法估计模型参数：

(1)仿真步骤

给定一个正整数G和一个序列V＝{v₁,v₂,L,v_M}，其中，v₁＝0，v_M是给定的正数，并且v₂,L,v_M-1服从均匀分布U(0,v_M)；对于试验应力S_k，k＝1,2,L,K，令

g＝1,2,L,G，表示K×G个服从标准正态分布的样本；因此，可仿真得到第k个协变量的应力水平如下

其中，m＝1,2,L,M；

(2)极大似然估计步骤

假设试验应力

有N个试验样本，退化状态测量值y(t_ijk)表示第j个样本在第k个应力下的第i个测量数据，对应的测量时间表示为t_ijk，其中，i＝1,2,L,n_k，j＝1,2,L,N，k＝1,2,L,K，n_k表示在第k个应力下测量值的个数，实际的退化状态表示为x(ti_jk)；

因此，y(ti_jk)可表示为

y(t_ijk)＝c_ajΛ(t_ijk)+σB(t_ijk)+ε (16)

其中，

c_aj表示漂移参数a_j的系数，Λ(t_ijk)可表示为

为便于描述，令Y_j＝(y(t_1j1),y(t_2j2),L,y(t_nkjk))′表示第j个样本的退化测量矢量，

其中，(·)′表示矢量的转置，Y表示全部的加速退化数据，包括Y_j，j＝1,2,L,N；

根据式(16)可知，Y_j服从多维正态分布，其均值和协方差可表示为

其中，

I_l是l维的单位矩阵，

因此，关于未知参数Θ(g,m)＝(μ_c,σ_c,σ_,σ_ε,θ,d_a)的对数似然函数可表示为

其中，

将对数似然函数(19)对μ_c求一阶偏导数，可得

令式(20)为0，可得

那么根据已估计的

参数σ_c，σ，σ_ε，θ和d_a的剖面似然函数可表示为

参数σ_c，σ，σ_ε，θ和d_a的极大似然估计值可利用多维搜索的方法，通过极大化剖面似然函数(22)得到，将

和

代入式(21)，即可得到

则极大似然估计值

即可得到；

(3)外推步骤

设备个体在工作应力下的退化数据已知，则可以根据贝叶斯方法更新退化模型的参数，从而得到该设备实时的剩余寿命预测值；

假设有κ个正常工作应力下的退化数据，表示为Y_κ＝[y₁,y₂,L,y_κ]，其中，y_κ表示设备在t_κ时刻的退化状态测量值；假设漂移系数a的先验估计分布服从正态分布，即

那么，联合后验分布p(a|Y_κ)可根据贝叶斯方法表示为

其中，p(Y_κ|a)为似然函数，p(a)为先验分布的PDF，Δy_g＝y_g-y_g-1，Δt_g＝t_g-t_g-1，

p(Y_κ|a)和p(a)均服从正态分布，那么后验分布p(a|Y_κ)与先验分布应属于同一分布族，即

其中，

将更新后的参数(25)代入(16)，即可得到设备个体剩余寿命的概率密度函数，实现剩余寿命的实时预测。

实施例应用分析

以某型惯性导航系统中的陀螺仪为例，验证本发明所提方法的有效性。考虑到温度是造成陀螺仪漂移的主要因素，漂移退化数据与温度应力的关系符合阿伦尼斯模型。在步进应力加速退化试验中，设定阈值为0.6°/h，步进应力为50℃，60℃和70℃。陀螺仪在实际工作环境中连续工作，每小时自动记录一次漂移退化数据，每个应力记录6次。六个样品的测试结果如图1所示。

根据步骤4中提出的参数估计方法，利用改进的MLE-SIMEX方法可以得到模型中未知参数的估计值，如表1所示。

表1退化模型参数估计值

对同类型同批次的陀螺仪个体在正常工作应力下进行测试，每1.5小时记录一次漂移退化数据，共得到14个数据，如图2所示。从图中可以看出，在第13个数据之后，退化轨迹首次达到失效阈值。因此，可认为样品的真实寿命L＝20.10小时。

在得到第i-th(1≤i≤13)个测量值后，可对样本进行剩余寿命(RUL)预测。在剩余寿命(RUL)预测过程中，一旦新的退化测量值可用，就可利用贝叶斯方法更新模型参数。具体的参数更新过程如图3所示。

为了说明该方法在剩余寿命(RUL)预测中的有效性，根据表1中模型参数的初始估计值和参数随退化数据的更新值，分别给出了基于实际观测退化数据在不同测量点下的RUL预测结果。剩余寿命(RUL)预测结果和相应相对误差见表2。

表2剩余寿命(RUL)预测值

从表2可以看出，随着模型参数的更新，剩余寿命(RUL)的预测结果更加接近于实际剩余寿命(RUL)。在第3次参数更新后，相对误差逐渐降低，稳定在2％左右。

此外，我们在图4中展示了剩余寿命(RUL)的概率密度函数(PDF)、预测的平均剩余寿命(RUL)和该时刻的实际剩余寿命(RUL)以进行比较。由图4可知，本发明方法所预测的剩余寿命与实际剩余寿命非常接近，并且所预测的剩余寿命的概率密度函数随着模型参数的更新变得更高更尖，说明预测结果的不确定性逐渐减小。

进一步验证了本发明提出的考虑多重不确定性的非线性加速退化建模与剩余寿命预测方法的有效性和优越性。

以上所述，仅用以说明本发明的技术方案而非限制，本领域普通技术人员对本发明的技术方案所做的其它修改或者等同替换，只要不脱离本发明技术方案的精神和范围，均应涵盖在本发明的权利要求范围当中。

Claims (5)

1.一种考虑多重不确定性的退化设备剩余寿命预测方法，其特征在于，包括以下步骤：

步骤1、针对退化设备，基于扩散过程建立一般的非线性加速退化过程模型，并利用标准布朗运动描述随机退化过程随时间变化的固有不确定性；

步骤2、利用阿伦尼斯模型来描述退化模型漂移系数与加速应力间的关系，同时考虑退化设备的个体差异性和加速应力作为协变量的测量不确定性，建立考虑多重不确定性的加速退化模型；

步骤3、考虑性能退化过程中的测量不确定性，基于首达时间的概念，推导退化设备剩余寿命分布的表达式；

步骤4、利用MLE-SIMEX方法估计虑多重不确定性的加速退化模型中的未知参数，将其代入退化设备剩余寿命分布的表达式，从而实现退化设备的寿命预测。

2.根据权利要求1所述的一种考虑多重不确定性的退化设备剩余寿命预测方法，其特征在于，在步骤1中，所述非线性加速退化过程模型建立如下：

设X(t)表示样品在t时刻的性能退化量，则基于扩散过程的退化过程{X(t),t≥0}可以表示为：

其中，μ(t；θ)为漂移系数，是时间t的非线性函数，用以表示模型的非线性特征，其中参数向量θ＝(a,b)为未知参数，a表征不同设备间的差异性，b表征同类设备的共性特征；σ_B称为扩散系数，B(t)为标准Brownian运动，且B(t):N(0,t)；

对于式(1)，在首达时间概念下，设备的寿命T可定义为：

T＝inf{t:X(t)≥ω|X(0)＜ω} (2)

当不考虑参数a和σ_B的随机性时，设备寿命的概率密度函数可以近似为：

其中，

ω为失效阈值。

3.根据权利要求1所述的一种考虑多重不确定性的退化设备剩余寿命预测方法，其特征在于，在步骤2中，考虑多重不确定性的加速退化模型的建立如下:

假设存在K步的步进加速退化试验数据，试验以最低的应力水平S₁开始，初始时间t＝0；每个设备都将在预定的时间经历预定的应力增加过程，当t＝T_k(k＝1,2,L,K)时，相应应力为S_k，其表达如下:

漂移系数a与加速应力相关，利用阿伦尼斯模型来描述a与加速应力S的关系如下：

a_k＝c_aexp(-d_a/S_k) (5)

其中，a_k表示漂移系数a在不同应力水平S_k下的值，c_a是关于失效模式的常数，d_a是关于激活能和玻尔兹曼常数的常数；

对于由标准布朗运动{B(t),t≥0}驱动的的非线性退化模型，随机过程退化模型随时间变化的固有属性可描述时变不确定性；由于不同的工作条件，退化轨迹将表现出不同的退化速率，可利用漂移系数可作为描述个体差异性的随机参数，漂移系数确定了设备在不同应力下的退化轨迹，考虑到a_k的随机效应，式(5)中，c_a可作为随机参数且假设c_a:

μ_c和

分别表示c_a的均值和方差；因此，能够得到a_k:

μ_ak和

分别表示a_k的均值和方差，即：

实际工程中，为描述测量不确定性的影响，令{Y(t),t≥0}表示设备的测量过程，用来描述隐藏退化状态与测量值之间的关系，具体如下：

Y(t)＝X(t)+ε (7)

其中，ε是随机测量误差，假设ε:

表示ε的方差；

通常加速退化试验在预定的应力水平下进行，由于测量误差的存在，使得设备实际上在未知的加速应力水平下发生退化，加速应力的测量误差存在于协变量的测量过程中，可将协变量的测量误差用随机变量

表示，且

表示

的方差，因此，实际加速应力水平可表示为

即：

综上所述，考虑多重不确定性的加速退化模型可表示为：

4.根据权利要求1所述的一种考虑多重不确定性的退化设备剩余寿命预测方法，其特征在于，在步骤3中，退化设备剩余寿命分布表达式推导如下：

在首达时间下，考虑多重不确定性的设备寿命T_e可定义为：

T_e＝inf{t:Y(t)≥ω|Y(0)＜ω} (10)

相应的，设备寿命的概率密度函数表示为

首先，在不考虑漂移系数随机效应的情况下推导

对于退化过程{Y(t),t≥0}，其寿命定义如式(10)所示，则寿命T_e的

可表示为

其中，B＝ω-a·h(t,θ)，C＝σ²t；

为了估计单个设备的剩余寿命，需要得到该设备在正常工作应力下的退化数据，Y(t_k)表示设备在t_k时刻的测量值，根据首达时间的概念，定义t_k时刻的剩余寿命L_k为

L_k＝inf{l_k>0:Y(l_k+t_k)≥ω} (12)

已知t_k时刻的退化状态测量值为Y(t_k)，则剩余寿命的概率密度函数可表示为

其中，

Ω(l_k)＝h(l_k+t_k,θ)-h(t_k,θ)；

利用漂移系数a的随机性表示退化设备的个体差异性，并且a服从式(6)的分布形式；

对于退化过程{Y(t),t≥0}，已知a:

并且t_k时刻的退化状态测量值为Y(t_k)，则t_k时刻剩余寿命的概率密度函数可表示为

其中，

L＝Ω(l_k)-l_kμ(l_k+t_k；θ)，M＝Ω(l_k)＝h(l_k+t_k,θ)-h(t_k,θ)，

5.根据权利要求1所述的一种考虑多重不确定性的退化设备剩余寿命预测方法，其特征在于，在步骤4中，利用MLE-SIMEX方法估计模型参数如下:

(1)仿真步骤

给定一个正整数G和一个序列V＝{v₁,v₂,L,v_M}，其中，v₁＝0，v_M是给定的正数，并且v₂,L,v_M-1服从均匀分布U(0,v_M)；对于试验应力S_k，k＝1,2,L,K，令u_kg，g＝1,2,L,G，表示K×G个服从标准正态分布的样本；因此，可仿真得到第k个协变量的应力水平如下

其中，m＝1,2,L,M；

(2)极大似然估计步骤

假设试验应力

有N个试验样本，退化状态测量值y(t_ijk)表示第j个样本在第k个应力下的第i个测量数据，对应的测量时间表示为t_ijk，其中，i＝1,2,L,n_k，j＝1,2,L,N，k＝1,2,L,K，n_k表示在第k个应力下测量值的个数，实际的退化状态表示为x(t_ijk)；

因此，y(t_ijk)可表示为

y(t_ijk)＝c_ajΛ(t_ijk)+σB(t_ijk)+ε (16)

其中，

c_aj表示漂移参数a_j的系数，Λ(t_ijk)可表示为

为便于描述，令Y_j＝(y(t_1j1),y(t_2j2),L,y(t_nkjk))′表示第j个样本的退化测量矢量，T_j＝(Λ(t_1j1),Λ(t_2j2),L,Λ(t_nkjk))′，其中，(·)′表示矢量的转置，Y表示全部的加速退化数据，包括Y_j，j＝1,2,L,N；

根据式(16)可知，Y_j服从多维正态分布，其均值和协方差可表示为

其中，

I_l是l维的单位矩阵，

因此，关于未知参数Θ(g,m)＝(μ_c,σ_c,σ_,σ_ε,θ,d_a)的对数似然函数可表示为

其中，

将对数似然函数(19)对μ_c求一阶偏导数，可得

令式(20)为0，可得

那么根据已估计的

参数σ_c，σ，σ_ε，θ和d_a的剖面似然函数可表示为

参数σ_c，σ，σ_ε，θ和d_a的极大似然估计值可利用多维搜索的方法，通过极大化剖面似然函数(22)得到，将

和

代入式(21)，即可得到

则极大似然估计值

即可得到；

(3)外推步骤

设备个体在工作应力下的退化数据已知，则可以根据贝叶斯方法更新退化模型的参数，从而得到该设备实时的剩余寿命预测值；

假设有κ个正常工作应力下的退化数据，表示为Y_κ＝[y₁,y₂,L,y_κ]，其中，y_κ表示设备在t_κ时刻的退化状态测量值；假设漂移系数a的先验估计分布服从正态分布，即a:

那么，联合后验分布p(a|Y_κ)可根据贝叶斯方法表示为

其中，p(Y_κ|a)为似然函数，p(a)为先验分布的PDF，Δy_g＝y_g-y_g-1，Δt_g＝t_g-t_g-1，

p(Y_κ|a)和p(a)均服从正态分布，那么后验分布p(a|Y_κ)与先验分布应属于同一分布族，即

其中，

将更新后的参数(25)代入(16)，即可得到设备个体剩余寿命的概率密度函数，实现剩余寿命的实时预测。

Images