[发明专利]一种基于退化数据建模的服役设备剩余寿命自适应预测方法

摘要

本发明公开了一种基于退化数据建模的服役设备剩余寿命自适应预测方法，运用Bayesian方法和EM算法实现对服役设备的退化建模和剩余寿命预测，包括以下步骤1）随机退化建模；2）基于Bayesian方法的模型随机参数更新；3）剩余寿命预测；4）基于EM算法的模型参数估计。本发明在服役设备剩余寿命预测的指数随机退化模型构建方法中引进了基于EM算法的参数自适应更新机制，使得指数随机退化模型的所有参数都随着服役设备实时数据的积累不断更新，因而预测的结果更能反映设备的实际运行情况，达到减小预测不确定性的目的。本发明不需要多个同类设备历史数据来初始化退化模型，即可自适应地更新模型参数和剩余寿命分布。

主权项

一种基于退化数据建模的服役设备剩余寿命自适应预测方法，其特征在于，运用Bayesian方法和EM算法实现对服役设备的退化建模和剩余寿命预测，包括以下步骤：1)随机退化建模；2)基于Bayesian方法的模型随机参数更新；3)剩余寿命预测；4)基于EM算法的模型参数估计；所述步骤1)具体为：令S(t)表示t时刻的退化量，设备在离散时间监测点t1，t2，...进行退化监测，在指数随机退化模型下，设备在tk时刻的退化量为：S(tk)=φ+θexp(βtk+ϵ(tk)-σ22)---(1)]]>其中φ是固定常数；θ和β是随机变量，用来描述设备间的个体差异；ε(tk)是服从正态分布的随机误差项，且有ε(tk)～N(0，σ2)；对指数随机模型进行对数变换，得：L(tk)=ln(S(tk)-φ)=lnθ-σ22+βtk+ϵ(tk)=θ′+βtk+ϵ(tk),---(2)]]>其中θ′＝lnθ‑σ2/2；假设随机参数θ′和β的先验分布服从联合的二元正态分布，分别具有均值μ0′和μ1、方差和且相关系数为ρ0；所述步骤2)具体为：基于步骤1)的模型和参数假设，运用Bayesian理论和服役设备退化的监测数据来对先验分布中的参数进行更新，得到相应的后验分布；记Lk＝L(tk)，L1：k＝{L1，…，Lk}；假设在t1，…，tk时刻有退化观测量L1，…，Lk，由于误差项ε(ti)，i＝1，…，k是独立同分布的随机变量，在给定θ′和β的条件下，退化样本L1：k的条件联合密度函数为：f(L1:k|θ′,β)=(12πσ2)k×exp(-Σi=1k(Li-θ′-βti)22σ2)---(3)]]>由于θ′和β的先验分布是二元正态分布，与样本分布f(L1：k|θ′，β)是共轭的，故θ′和β关于L1：k的联合后验分布仍然是二元正态分布，即因此，有：f(θ′,β|L1:k)∝f(L1:k|θ′,β)~f(θ′,β)∝exp[-Σj=1k(Lj-θ′-βtj)22σ2]~exp[(θ′-μ0′)2σ02-2ρ02(θ′-μ0′)(β-μβ)σ0σ1+(β-μβ)2σ12]∝12πσθ′,kσβ,k1-ρk2exp[-12(1-ρk2)((θ′-μθ′,k)2σθ′,k2-2ρk(θ′-μθ′,k)(β-μβ,k)σθ′,kσβ,k+(β-μβ,k)2σβ,k2)]---(4)]]>上式对应的参数由下面命题给出；命题1：给定到tk时刻的所有退化观测数据L1：k＝{L1，…，Lk}，θ′，β关于L1：k的联合后验分布服从二元正态分布，即：具体参数表达式如下：μθ′,k=μ0′σ2σ1(Y+Zρ0)-μ1σ2σ0(Yρ0+Z)+(1-ρ02)σ0σ1(Yσ0Σi=1kLi-Zσ1Σi=1kLiti)σ1(XY-Z2)]]>μβ,k=μ1σ2σ0(X+Zρ0)-μ0′σ2σ1(Xρ0+Z)+(1-ρ02)σ0σ1(Xσ1Σi=1kLiti-Zσ0Σi=1kLi)σ0(XY-Z2)]]>σθ′,k2=[(1-ρ02)σ12Σi=1kti2+σ2](1-ρ02)σ2σ02[k(1-ρ02)σ02+σ2][(1-ρ02)σ12Σi=1kti2+σ2]-[(1-ρ02)σ0σ1Σi=1kti+ρ0σ2]2]]>σβ,k2=[k(1-ρ02)σ02+σ2](1-ρ02)σ2σ12[k(1-ρ02)σ02+σ2][(1-ρ02)σ12Σi=1kti2+σ2]-[(1-ρ02)σ0σ1Σi=1kti+ρ0σ2]2]]>ρk=-((1-ρ02)σ0σ1Σi=1kti-ρ0σ2[k(1-ρ02)σ02+σ2][(1-ρ02)σ12Σi=1kti2+σ2])---(5)]]>其中，由于是二元正态分布的，根据二元正态分布的性质，有如下推论成立；命题2：令～|L1：k表示条件随机变量，E(～L1：k)表示条件期望，则：①E(θ′β|L1：k)＝ρkσθ′，kσβ，k+μθ′，kμβ，k，E(θ′|L1：k)＝μθ′，k，②μβ|θ′，k＝μβ，k+ρkσβ，k(θ′‑μθ′，k)/σθ′，k，其中：所述步骤3)具体为：在得到θ′和β的后验估计后，对于给定的失效阈值ω，要实现剩余寿命的预测，首先要基于到当前时刻tk的退化观测数据L1：k预测t+tk时的退化量；对于将来时刻t+tk的退化量L(t+tk)，在给定L1：k的条件下，有以下命题成立；命题3：给定到tk时刻的所有退化观测数据L1：k，预测的t+tk时刻对应的退化量L(t+tk)服从正态分布，且具有的均值和方差分别为：μ~(t+tk)=μθ′,k+μβ,k(t+tk)-σ22,σ~2(t+tk)=σθ′,k2+σβ,k2(t+tk)+σ2+2ρc(t+tk)σθ′,kσβ,k.---(6)]]>根据以上结果，设T为设备在tk时刻的剩余寿命，即T满足L(t+tk)＝lnω；于是，在给定L1：k的条件下，预测的剩余寿命的条件概率分布函数为：FT|L1:k(t)=Pr(T≤t|LL1:k)=Pr(S(t+tk)≥lnω|LL1:k)=Pr(Z≥lnw-μ~(t+tk)σ~2(t+tk))=Φ(g(t)).---(7)]]>其中Z服从标准正态分布，Φ(·)为标准正态随机变量的累积分布函数，且有由于limt→‑∞g(t)＝‑μβ，k/σβ，k，所以T的值域为(‑∞，∞)；考虑到T表示设备的寿命是非负的实数，将在(0，∞)上截断的累积分布函数作为设备的剩余寿命预测结果，如下式：FT|L1:k,T≥0(t)=Pr(T≤t|L1:k,T≥0)=Pr(0≤T≤t|T1:k)Pr(T≥0|L1:k)=Φ(g(t))-Φ(g(0))1-Φ(g(0)).---(8)]]>基于以上的结果，剩余寿命预测的条件概率密度函数为：fT|L1:k,T≥0(t)=dFT|L1:k,T≥0(t)dt=φ(g(t))1-Φ(g(0))·g′(t).---(9)]]>其中φ(·)为标准正态随机变量的概率密度函数；通过上面的建模过程，在获得新的退化测量数据后，根据命题1重新计算随机参数θ′和β的后验分布，由此获得的新估计；然后，根据这些更新的参数，通过更新函数g(t)中和的值实现对预测的剩余寿命分布的更新；所述步骤4)具体为：在以上的退化建模和剩余寿命预测过程中，模型参数是未知的，令Θ＝[σ2，μ′0，μ1，σ02，σ12，ρ0]表示模型未知参数；基于极大似然估计的方法，在tk时刻的监测数据Lk得到后，计算关于L1：k的对数似然函数：其中p(L1：k|Θ)表示退化数据L1：k的联合密度函数；Θ的极大似然估计由下式得到：在退化模型(2)中，参数θ′和β是随机的，采用EM算法实现对Θ的估计；EM算法对Θ的估计通过迭代以下两步实现：Ⅰ E‑step:计算条件期望其中表示基于退化数据L1：k进行估计时第i次迭代的结果；Ⅱ M‑step：最大化条件期望为了表示估计的参数依赖于到当前时刻的所有退化监测数据，将基于到tk时刻的所有退化监测数据L1：k估计的参数表示为令EM算法中第i次迭代得到的估计为完全对数似然函数为lnp(L1:k,θ′,β|Θk)=lnp(L1:k|θ′,β,Θk)+lnp(θ′,β|Θk)=-k+22ln2π-k2lnσk2-Σj=1k(Lj-θ′-βtj)22σk2-lnσ0,k2-lnσ1,k2-ln1-ρ0,k2-12(1-ρ0,k2)[(θ′-μ0,k′)2σ0,k2-2ρ0,k2(θ′-μ0,k′)(β-μ1,k)σ0,kσ1,k+(β-μ1,k)2σ1,k2]---(15)]]>基于式(13)，得到如下：lk(Θ|Θ^k(i))=Eθ,β|L1:k,Θ^k(i){logp(L1:k,θ′,β|Θ)}=-k+22ln2π-k2lnσk2-lnσ0,k2-lnσ1,k2-ln1-ρ0,k2-Σj=1kLj2-2Lj(μθ′,k+μβ,ktj)+μθ′,k2+σθ′,k2+2tj(ρkσθ′,kσβ,k+μθ′,kμβ,k)+tj2(μβ,k2+σβ,k2)2σk2-12(1-ρ0,k2)[μθ′,k2+σθ′,k2-2μ0,k′μθ′,k+μ0,k′2σ0,k2-2ρ0,k2(ρkσθ′,kσβ,k+μθ′,kμβ,k-μ1,kμθ′,k-μ0,k′μβ,k+μ0,k′μ1,k)σ0,kσ1,k+μβ,k2+σβ,k2-2μ1,kμβ,k+μ1,k2σ1,k2]---(16)]]>令得到第i+1步的参数估计为：σ^k2(i+1)=1kΣj=1k(Lj2-2Lj(μθ′,k+μβ,ktj)+μθ′,k2+σθ′,k2+2tj(ρkσθ′,kσβ,k+μθ′,kμβ,k)+tj2(μβ,k2+σβ,k2))μ^0,k′(i+1)=2σ0,kσ1,k·ρ0,k2(μ1,k-μβ,k)+μθ′,kμ^1,k(i+1)=2σ1,kσ0,k·ρ0,k2(μ0,k-μθ′,k)+μβ,kσ^0,k2(i+1)=14(-12B+14B2+4A)2σ^1,k2(i+1)=14(-12D+14D2+4C)2ρ^0,k(i+1)=1-2E---(17)]]>其中：A=μθ′,k2+σθ′,k2-2μ0,k′μθ′,k+μ0,k′22(1-ρ0,k2)C=μβ,k2+σβ,k2-2μ1,kμβ,k+μ1,k22(1-ρ0,k2)]]>B=ρ0,k2(ρkσθ′,kσβ,k+μθ′,kμβ,k-μ1,kμθ′,k-μ0,k′μβ,k+μ0,k′μ1,k)σ1,k(1-ρ0,k2)D=ρ0,k2(ρkσθ′,kσβ,k+μθ′,kμβ,k-μ1,kμθ′,k-μ0,k′μβ,k+μ0,k′μ1,k)σ0,k(1-ρ0,k2)]]>E=ρkσθ′,kσβ,k+μθ′,kμβ,k-μ1,kμθ′,k-μ0,k′μβ,k+μ0,k′μ1,kσ0,kσ1,k.]]>