CN106874634A - 基于逆高斯退化模型的剩余寿命贝叶斯预测方法 - Google Patents

Abstract

本发明涉及一种基于逆高斯退化模型的剩余寿命贝叶斯预测方法，该方法主要包括以下步骤：（1）建立逆高斯退化模型与剩余寿命预测模型；（2）设逆高斯退化模型的尺度参数与均值参数为随机参数，构建随机参数的共轭先验分布函数；（3）根据Bayes理论推导随机参数的后验分布函数；（4）设计EM算法估计随机参数先验分布函数的超参数值；（5）估计随机参数的后验期望值；（6）预测产品的剩余寿命。该方法可有效融合先验退化数据与现场退化数据，提高剩余寿命预测值的准确度与置信度。

Description

基于逆高斯退化模型的剩余寿命贝叶斯预测方法

技术领域

本发明属于可靠性工程领域，涉及一种基于逆高斯退化模型的剩余寿命贝叶斯预测方法。

背景技术

当代社会对产品的可靠性要求越来越高，产品的质量监控与健康管理也朝着精细化方向发展，因此，如何提前准确预测出产品的剩余寿命成为了目前的一个研究热点。逆高斯退化模型具有优良的统计特性，适合对单调性能退化过程建模，已经广泛用于退化失效型产品的剩余寿命预测。为了有效融合多源退化数据提高剩余寿命的预测准确度和可信度，发明一种基于逆高斯退化模型的剩余寿命贝叶斯预测方法。

发明内容

发明一种基于逆高斯退化模型的剩余寿命贝叶斯预测方法，该方法的具体技术方案为：

步骤一：建立逆高斯退化模型与剩余寿命预测模型

设产品的性能退化过程{Y(t),t≥0}服从逆高斯退化模型，则Y(t)～IG(μΛ(t),λΛ(t)²)，Y(t)的概率密度函数表示为

式中，μ为均值参数，λ为尺度参数，Λ(t)＝t^Λ为时间函数。设产品性能指标的失效阈值为D，则产品的寿命ξ为Y(t)首次到达D的时间，ξ的累积分布函数为

式中，Φ(·)为标准Normal分布的累积分布函数。产品的剩余寿命L(t)是指从时刻t的性能退化量Y(t)到首次超过D的时间，可以表示为L(t)＝inf{x|Y(t+x)≥D,x>0}，则产品寿命ξ与剩余寿命L(t)之间的关系为ξ＝L(t)+t。

步骤二：设逆高斯退化模型的尺度参数与均值参数为随机参数，构建随机参数的共轭先验分布函数

为了有效融合先验退化数据有现场退化数据，提高剩余寿命预测值的准确度与可信度，设μ,λ作为随机参数。并且为了便于统计分析，采用μ,λ的共轭先验分布：设λ服从Gamma分布λ～Ga(a,b)，其概率密度函数表示为

设δ＝1/μ服从条件正态分布δ|λ～N(c,d/λ)，其概率密度函数为

其中，a,b,c,d是随机参数的超参数。

步骤三：根据Bayes理论推导随机参数的后验分布函数

设y_i,j为第i个产品的第j次性能退化测量值，△y_i,j为退化增量，△Λ_i,j为时间增量，根据逆Gaussian过程的统计特性建立似然函数为

建立完全对数似然函数为

设f(δ_i,λ_i)是随机参数δ_i,λ_i的联合先验概率密度函数，则有f(δ_i,λ_i)＝f(λ_i)·f(δ_i|λ_i)，通过Bayesian公式f(δ_i,λ_i|△y_i)∝L(δ_i,λ_i)·f(δ_i,λ_i)推导出联合后验概率密度函数f(δ_i,λ_i|△y_i)，得到随机参数δ_i,λ_i的后验分布为

其中，

步骤四：设计EM算法估计随机参数先验分布函数的超参数值

设计EM算法估计超参数Ω＝(a,b,c,d)，EM算法每一步迭代包含E步和M步。E步的任务是求取隐含数据项的期望值，设Ω^(l)为第l次迭代后的估计值向量，则在第l+1次迭代中，隐含数据项λ_i,lnλ_i,λ_iδ_i,的期望值为

式中，ψ(·)为digamma函数。将完全似然函数式(1)中的各隐含数据项利用对应的期望值代替后，M步的任务是极大化式(1)，解得c^(l+1)，d^(l+1)，b^(l+1)及a^(l+1)的表达式如下

式中，ψ^-1(·)为逆digamma函数。

EM算法的执行过程为：

初始化:设l＝0,Ω⁽⁰⁾＝(1,1,1,1)；

第l+1次迭代:

E步:计算E(λ_i|y_i,Ω^(l)),E(lnλ_i|y_i,Ω^(l)),E(λ_iδ_i|y_i,Ω^(l))及

M步:解得c^(l+1)，d^(l+1)，b^(l+1)及a^(l+1)，将Ω^(l)更新为Ω^(l+1)；

结束条件:max(Ω^(l+1)-Ω^(l))<10^-3或l达到最大迭代数。

步骤五：估计随机参数的后验期望值

设产品在时间t_j的现场性能退化数据为△x_j，j＝1,2,…,k，首先将随机参数在时间t_j的后验分布函数更新为

其中，

然后，根据Gamma函数与Normal函数的统计特性估计随机参数的后验期望值为

步骤六：预测产品的剩余寿命

利用随机参数的后验期望值更新累积失效分布函数

利用F_ξ(t|△x)计算出E(ξ|△x)，产品在t时刻的剩余寿命预测值为L(t)＝E(ξ|△x)-t。

附图说明

图1基于逆高斯退化模型的剩余寿命贝叶斯预测方法。

具体实施方式

下面结合附图对本发明实现步骤进行进一步说明。

实施例：电连接器的主要失效模式有机械失效，电气失效，绝缘失效三种，机械失效主要由接插件应力松弛造成，性能退化量y为接插件应力值x相对于初始应力值x₀的百分比变化y＝(x-x₀)/x₀×％，每个样品在0时刻的性能退化量为0，失效阈值为D＝30％。电连接器的6组先验性能退化数据及测量时间如表1所示，1组现场性能退化数据如表2所示。

表1电连接器的先验性能退化数据

表2电连接器的现场性能退化数据

步骤一：建立逆高斯退化模型与剩余寿命预测模型，产品性能指标的失效阈值为D＝30％，则产品的寿命ξ的累积分布函数为

产品的剩余寿命L(t)是指从时刻t的性能退化量Y(t)到首次超过D的时间L(t)＝inf{x|Y(t+x)≥0.3,x>0}。

步骤二：设逆高斯退化模型的尺度参数与均值参数为随机参数，构建随机参数的共轭先验分布函数。设μ,λ作为随机参数。并且为了便于统计分析，采用μ,λ的共轭先验分布：设λ服从Gamma分布λ～Ga(a,b)，其概率密度函数表示为

设δ＝1/μ服从条件正态分布δ|λ～N(c,d/λ)，其概率密度函数为

其中，a,b,c,d是随机参数的超参数。

步骤三：根据Bayes理论推导随机参数的后验分布函数。设y_i,j为第i个产品的第j次性能退化测量值，△y_i，j为退化增量，△Λ_i，j为时间增量，根据逆Gaussian过程的统计特性建立似然函数为

建立完全对数似然函数为

设f(δ_i,λ_i)是随机参数δ_i,λ_i的联合先验概率密度函数，则有f(δ_i,λ_i)＝f(λ_i)·f(δ_i|λ_i)，通过Bayesian公式f(δ_i,λ_i|△y_i)∝L(δ_i,λ_i)·f(δ_i,λ_i)推导出联合后验概率密度函数f(δ_i,λ_i|△y_i)，得到随机参数δ_i,λ_i的后验分布为

其中，

步骤四：设计EM算法估计随机参数先验分布函数的超参数值，EM算法的执行过程为：

初始化:设l＝0,Ω⁽⁰⁾＝(1,1,1,1)；

第l+1次迭代:

E步:计算E(λ_i|y_i,Ω^(l)),E(lnλ_i|y_i,Ω^(l)),E(λ_iδ_i|y_i,Ω^(l))及

M步:解得c^(l+1)，d^(l+1)，b^(l+1)及a^(l+1)，将Ω^(l)更新为Ω^(l+1)；

结束条件:max(Ω^(l+1)-Ω^(l))<10^-3或l达到最大迭代数。

解得

步骤五：估计随机参数的后验期望值

利用产品在10个测量时刻的现场性能退化数据，结合先验参数估计值分别对超参数的后验估计值进行更新，更新结果如表3所示

表3超参数值更新结果

根据Gamma函数与Normal函数的统计特性估计出，产品在时刻1800h的随机参数后验期望值为E(λ|△x)＝0.647，E(μ|△x)＝0.601。

步骤六：预测产品的剩余寿命

利用随机参数的后验期望值更新累积失效分布函数

F_ξ(t|△x)＝Φ(0.147(t^0.509-49.917))-exp(2.153t^0.509)Φ(-0.147((t^0.509+49.917)))，

进而计算出E(ξ|△x)＝2169.906h，产品在1800h的剩余寿命预测值为L(t)＝369.906h，利用Bootstrap自助抽样法建立剩余寿命预测值的90％置信度区间为[158h,623h]。

如果不利用Bayes理论融合先验性能退化数据进行剩余寿命预测，只利用现场性能退化数据预测出

F_ξ(t|△x)＝Φ(0.153(t^0.482-61.004))-exp(2.869t^0.482)Φ(-0.153((t^0.482+61.004)))，

剩余寿命值为3259h，原大于实际剩余寿命560h；利用Bootstrap自助抽样法建立剩余寿命预测值的90％置信度区间为[626h,6755h]，区间长度远大于发明方法所建立的置信区间长度。可见，发明方法提高了剩余寿命预测值的准确度与置信度。

Claims (3)

1.基于逆高斯退化模型的剩余寿命贝叶斯预测方法，其特征在于，包括以下步骤：

步骤一，建立逆高斯退化模型与剩余寿命预测模型；

步骤二，设逆高斯退化模型的尺度参数λ与均值参数μ为随机参数，构建随机参数的共轭先验分布函数为：

设λ服从Gamma分布λ～Ga(a,b)，先验分布函数为

$f\left(\mathrm{\&lambda;}\right)=\frac{{\mathrm{\&lambda;}}^{a-1}}{\mathrm{\&Gamma;}\left(a\right)b^a}\exp (-\frac{\mathrm{\&lambda;}}{b}),$

设δ＝1/μ服从条件正态分布δ|λ～N(c,d/λ)，先验分布为

$f\left(\mathrm{\&delta;}\right|\mathrm{\&lambda;})=\sqrt{\frac{\mathrm{\&lambda;}}{2\mathrm{\&pi;}d}}\exp (-\frac{\mathrm{\&lambda;}{(\mathrm{\&delta;}-c)}^2}{2d}),$

式中，Γ(·)为Gamma函数，a,b,c,d为随机参数的超参数；

步骤三，根据Bayes理论推导随机参数的后验分布函数，得到随机参数δ_i,λ_i的后验分布函数为

${\mathrm{\&lambda;}}_i|{\mathrm{\&Delta;y}}_i~Ga(\frac{m_i}{2}+a,1/A_i)$

${\mathrm{\&delta;}}_i|{\mathrm{\&Delta;y}}_i,{\mathrm{\&lambda;}}_i~N(\frac{{\mathrm{\&Sigma;}}_{j=1}^{m_i}{\mathrm{\&Delta;\&Lambda;}}_{i,j}+c/d}{{\mathrm{\&Sigma;}}_{j=1}^{m_i}{\mathrm{\&Delta;y}}_{i,j}+1/d},\frac{1}{{\mathrm{\&lambda;}}_i({\mathrm{\&Sigma;}}_{j=1}^{m_i}{\mathrm{\&Delta;y}}_{i,j}+1/d)});$

式中，△y_i为性能退化增量，△Λ_i,j表示第i个产品的第j个时间增量，m_i第i个产品性能退化数据测量总数，

步骤四，设计EM算法估计随机参数先验分布函数的超参数值；

步骤五，估计随机参数的后验期望值；

步骤六，预测产品的剩余寿命。

2.如权利要求1所述的基于逆高斯退化模型的剩余寿命贝叶斯预测方法，其特征在于，步骤四中，利用EM算法估计超参数值的流程为：

初始化:设l＝0,Ω⁽⁰⁾＝(1,1,1,1)；

第l+1次迭代:

E步:计算E(λ_i|y_i,Ω^(l)),E(lnλ_i|y_i,Ω^(l)),E(λ_iδ_i|y_i,Ω^(l))及E(λ_iδ_i ²|y_i,Ω^(l))；

以上各项的表达式推导为：

$E\left({\mathrm{\&lambda;}}_i\right|{\mathrm{\&Delta;y}}_i,{\mathrm{\&Omega;}}^{\left(l\right)})=(a^{\left(l\right)}+m_i/2)/A_i^{\left(l\right)},$

$E\left({\mathrm{ln\&lambda;}}_i\right|{\mathrm{\&Delta;y}}_i,{\mathrm{\&Omega;}}^{\left(l\right)})=\mathrm{\&psi;}(a^{\left(l\right)}+\frac{m_i}{2})-{\mathrm{lnA}}_i^{\left(l\right)},$

$E\left({\mathrm{\&lambda;}}_i{\mathrm{\&delta;}}_i\right|{\mathrm{\&Delta;y}}_i,{\mathrm{\&Omega;}}^{\left(l\right)})=\frac{a^{\left(l\right)}+m_i/2}{A_i^{\left(l\right)}}\&CenterDot;\frac{{\mathrm{\&Sigma;}}_{j=1}^n{\mathrm{\&Delta;\&Lambda;}}_{i,j}+c^{\left(l\right)}/d^{\left(l\right)}}{{\mathrm{\&Sigma;}}_{j=1}^{m_i}{\mathrm{\&Delta;y}}_{i,j}+1/d^{\left(l\right)}},$

$E\left({\mathrm{\&lambda;}}_i{\mathrm{\&delta;}}_i^2\right|{\mathrm{\&Delta;y}}_i,{\mathrm{\&Omega;}}^{\left(l\right)})=\frac{a^{\left(l\right)}+m_i/2}{A_i^{\left(l\right)}}{\left(\frac{{\mathrm{\&Sigma;}}_{j=1}^n{\mathrm{\&Delta;\&Lambda;}}_{i,j}+c^{\left(l\right)}/d^{\left(l\right)}}{{\mathrm{\&Sigma;}}_{j=1}^{m_i}{\mathrm{\&Delta;y}}_{i,j}+1/d^{\left(l\right)}}\right)}^2+\frac{1}{{\mathrm{\&Sigma;}}_{j=1}^{m_i}{\mathrm{\&Delta;y}}_{i,j}+1/d^{\left(l\right)}};$

式中，Ω^(l)＝(a,b,c,d)^(l)为第l次迭代后得到的超参数向量，

M步:解得c^(l+1)，d^(l+1)，b^(l+1)及a^(l+1)，将Ω^(l)更新为Ω^(l+1)；

以上各项的表达式为：

$c^{(l+1)}=\frac{{\mathrm{\&Sigma;}}_{i=1}^n{m}_{i}E\left({\mathrm{\&lambda;}}_i{\mathrm{\&delta;}}_i\right|{\mathrm{\&Delta;y}}_i,{\mathrm{\&Omega;}}^{\left(l\right)})}{{\mathrm{\&Sigma;}}_{i=1}^n{m}_{i}E\left({\mathrm{\&lambda;}}_i\right|{\mathrm{\&Delta;y}}_i,{\mathrm{\&Omega;}}^{\left(l\right)})},$

$d^{(l+1)}=\frac{1}{{\mathrm{\&Sigma;}}_{i=1}^n{m}_i}{\mathrm{\&Sigma;}}_{i=1}^n{m}_i\left(E\right({\mathrm{\&lambda;}}_i{\mathrm{\&delta;}}_i^2|{\mathrm{\&Delta;y}}_i,{\mathrm{\&Omega;}}^{\left(l\right)})-2c^{(l+1)}E\left({\mathrm{\&lambda;}}_i{\mathrm{\&delta;}}_i\right|{\mathrm{\&Delta;y}}_i,{\mathrm{\&Omega;}}^{\left(l\right)})+c^{2(l+1)}E\left({\mathrm{\&lambda;}}_i\right|{\mathrm{\&Delta;y}}_i,{\mathrm{\&Omega;}}^{\left(l\right)})),$

$b^{(l+1)}=\frac{{\mathrm{\&Sigma;}}_{i=1}^n{m}_{i}E\left({\mathrm{\&lambda;}}_i\right|{\mathrm{\&Delta;y}}_i,{\mathrm{\&Omega;}}^{\left(l\right)})}{a^{\left(l\right)}{\mathrm{\&Sigma;}}_{i=1}^n{m}_i},$

$a^{(l+1)}={\mathrm{\&psi;}}^{-1}(\frac{{\mathrm{\&Sigma;}}_{i=1}^n{m}_{i}E\left({\mathrm{ln\&lambda;}}_i\right|{\mathrm{\&Delta;y}}_i,{\mathrm{\&Omega;}}^{\left(l\right)})}{{\mathrm{\&Sigma;}}_{i=1}^n{m}_i}-{\mathrm{lnb}}^{(l+1)}),$

结束条件:max(Ω^(l+1)-Ω^(l))<10^-3或l达到最大迭代数。

3.如权利要求1所述的基于逆高斯退化模型的剩余寿命贝叶斯预测方法，其特征在于，步骤五中，通过如下步骤获取随机参数的后验期望值：

首先，将随机参数在时间t_j的后验分布函数更新为

$\mathrm{\&lambda;}|\mathrm{\&Delta;}x~Ga(\frac{k}{2}+\hat{a},1/\hat{A})$

$\mathrm{\&delta;}|\mathrm{\&Delta;}x,\mathrm{\&lambda;}~N(\frac{{\mathrm{\&Sigma;}}_{j=1}^k{\mathrm{\&Delta;\&Lambda;}}_j+\hat{c}/\hat{d}}{{\mathrm{\&Sigma;}}_{j=1}^k{\mathrm{\&Delta;x}}_j+1/\hat{d}},\frac{1}{\mathrm{\&lambda;}({\mathrm{\&Sigma;}}_{j=1}^k{\mathrm{\&Delta;x}}_j+1/\hat{d})}),$

然后，根据Gamma函数与Normal函数的统计特性估计随机参数的后验期望值为