摘要
为获取精确的产品寿命预测结果,在可靠性模型中引入服从非共轭分布的随机参数以描述产品的个体差异。针对单调退化的长寿命产品,提出了采用非共轭分布随机参数Gamma过程的可靠性模型;将传统的随机参数共轭分布扩展到非共轭分布,得到不同参数分布下统一的寿命分布函数表达式,利用蒙特卡罗算法计算寿命预测值;利用马尔可夫链蒙特卡罗算法估计模型参数,提出不同分布下似然函数值计算方法及参数分布的选取准则,分析不同随机参数分布对产品可靠性的影响。仿真算例和应用实例表明,该方法能够选择合适的参数分布。同时,验证了该方法参数估计及寿命预测的精确性。
Abstract
To obtain accurate product lifetime prediction results, random parameters which obey non-conjugate distributions were introduced to describe the unit-to-unit variation of products in the reliability model. For long-life products with monotonic degradation, the reliability model was proposed based on the Gamma process with non-conjugate distribution random parameters. Extending the traditional conjugate distribution of random parameters to non-conjugate distribution, unified expression of life distribution function under different parameter distributions was derived, and the lifetime prediction values were calculated by using the Monte Carlo method. Markov chain Monte Carlo method was adopted for parameter estimation. The calculation algorithm for the likelihood function under different distributions and the selection criterion for the parameter distribution were presented, and the effect of the distribution of random parameters on the product reliability was analyzed. The simulation test and the application example show that this modeling method can select an appropriate distribution. Simultaneously, the accuracy of parameter estimation and lifetime prediction of the method is validated.
随着产品可靠性水平的不断提高,目前在进行寿命试验时,常常难以观测到产品的失效时间数据。对于这类长寿命产品,工业界广泛使用退化模型分析产品性能退化数据,在此基础上进行可靠性建模,进而预测产品的使用寿命及贮存寿命。由于退化模型对寿命预测的精确性至关重要,国内外关于退化模型的研究一直是可靠性领域的研究热点[1-6]。
基于随机过程的退化模型将产品的性能退化视为与时间相关的随机过程,在各类产品的可靠性建模中得到广泛应用,如激光器[7]、发光二极管[8-9]、输油管道[10]等。和Wiener随机过程相比,Gamma过程是单调的随机过程,常用于对单调退化产品的性能退化进行建模。由于同类产品的不同个体之间存在显著差异[11-15],目前基于Gamma过程的退化模型中,通常将Gamma过程的尺度参数视为随个体变化的随机参数,以描述这种个体间的差异。Lawless等[11]和Yang等[12]在假设尺度参数服从共轭Gamma分布的基础上,进行了可靠性建模并应用于某金属及伺服阀的寿命预测。Tsai等[13]在尺度参数共轭分布的假设基础上,利用期望最大化(expectation maximization,EM)算法估计模型参数并对退化试验方案进行优化设计。Giorgio等[14]与Rodríguez-Picón等[15]在建模过程中引入了更多的随机参数,但依然假设随机参数服从Gamma分布。
上述研究均使用尺度参数服从共轭Gamma分布的假设,利用共轭分布的统计特性使得EM算法求解模型参数简易可行。除共轭Gamma分布外,适合描述产品随机参数的分布还包括非共轭分布(包括Weibull分布、对数正态分布等)。实际上,Gamma分布的假设并不适用于所有产品[16-17],将非共轭分布误指定为共轭分布会使产品寿命的预测不够精确。Wang等[16]及杨志远等[17]均选取非共轭分布作为随机参数分布,并采用马尔可夫链蒙特卡罗(Markov chain Monte Carlo,MCMC)算法估计未知参数,但相关研究用于产品单一个体的剩余寿命预测,不能分析得出非共轭分布下产品总体的寿命分布,应用于预测产品总体寿命的场合受限。Ntzoufras[18]同样指出实际工程应用中存在无法获取共轭分布或共轭分布不符合实际的情况,引入非共轭分布并采用MCMC算法获取参数估计值,但相关研究仅用于寿命数据的贝叶斯建模分析中,没有考虑随机参数非共轭分布下针对单调退化产品的可靠性建模方法。
针对上述问题,为提高产品寿命预测精确性,提出采用非共轭分布随机参数Gamma过程的可靠性模型,采用Gamma随机过程描述产品退化过程,同时考虑非共轭分布与共轭分布两种情形,建立两种情形下统一的产品寿命分布表达式,并利用统一表达式得出共轭分布下产品寿命分布闭合解析式,针对非共轭分布提出参数估计方法及似然函数值计算方法,建立参数分布的选取准则,最后通过仿真算例和应用实例验证所提方法的有效性。
1 基于Gamma过程的可靠性模型
1.1 无随机参数的可靠性模型
假设产品退化量Y(t)服从Gamma过程,不考虑个体差异引起的参数随机分布,Y(t)服从Gamma分布Ga(α(t),β),概率密度函数为
(1)
其中:β>0为尺度参数;α(t)>0为形状参数,常见的函数式为α(t)=αtρ,α代表不受时间影响的形状参数大小,ρ=1时为线性退化,否则为非线性退化;为Gamma函数。
退化增量ΔY(t)=Y(t+Δt)-Y(t)服从Gamma分布Ga(α(t+Δt)-α(t),β)。不考虑随机参数分布,无交叉时间增量下的退化增量相互独立,并且由式(1)可知ΔY(t)>0。
考虑产品性能退化量随时间逐渐增长的情形,退化量阈值为ω,当退化量Y(t)首次达到退化阈值ω时就认为产品已经失效,产品寿命T定义为
(2)
依据Gamma过程的单调递增特性,不考虑随机参数分布,由式(2)可得产品的寿命分布函数为
(3)
其中,为下不完全Gamma函数。
产品的可靠度函数与寿命分布函数之和为1,因此无随机参数下的可靠度函数为
(4)
引入随机参数分布后同样可以利用寿命分布函数得出可靠度函数,后续不再给出可靠度表达式。
退化试验在恒定应力水平下测试n个试验样本的性能退化量,测试m次,测试的时间节点为tk(k=1,2,···,m),测得性能退化增量数据为Δyjk=yjk-yj(k-1)(j=1,2,···,n)。利用增量数据的相互独立特性,可得不考虑随机参数分布下的似然函数为
(5)
其中,。
1.2 非特定分布下的可靠性模型
考虑产品个体差异,认为尺度参数β随个体变化而变化,需要确定β的分布类型。Gamma过程参数β的共轭分布为Gamma分布,不指定β服从这一特定分步,认为β可能服从共轭Gamma分布或非共轭分布(如Weibull分布)。此时Y(t)的条件概率密度函数为
(6)
Y(t)的概率密度函数为
(7)
其中,f(β|δ,γ)代表参数β的概率密度函数,后续简化表示为f(β),δ为形状参数,γ为尺度参数。
不指定随机参数β服从特定分布,非共轭分布及共轭分布下统一的寿命分布函数为
(8)
由式(8)可知,产品寿命分布函数F2(t)实际是随机参数β的函数F1(t)的数学期望。β服从非共轭分布时,式(8)难以积分解析计算,可以采用蒙特卡罗(Monte Carlo,MC)算法计算非特定分布下的寿命预测值。
步骤1:选定β的具体分布类型,利用退化数据,计算模型未知参数向量θ=(α,γ,δ,ρ)的估计值。
步骤2:从f(β|δ,γ)分布中依次独立抽取的样本值(b=1,2,···,B),抽取次数B可设定为B=5 000。
步骤3:将代入式(3)计算的样本值。
步骤4:利用大样本平均值收敛于数学期望的性质,统一的寿命分布函数估计值为
(9)
考虑产品个体差异,给定参数β的条件下增量数据相互独立,因此似然函数为
(10)
1.3 共轭分布下的可靠性模型
共轭分布下式(7)~(10)可以积分求解得到闭合解析式,因此共轭分布下的可靠性模型是非特定分布下可靠性模型的特例。若考虑产品个体差异并指定β服从共轭Gamma分布,β的概率密度函数为
(11)
由式(7)可知Y(t)的概率密度函数为
(12)
取Z(t)=δY(t)/(γα(t)),Z(t)的概率密度函数为
(13)
由式(13)可知Z(t)服从自由度为(2α(t),2δ)的F分布,进一步由式(3)可得
(14)
利用式(10),得出共轭分布下似然函数为
(15)
其中: ;,一般有t0=0,yj0=0(j=1,2,···,n)。
2 参数估计方法及分布选取准则
2.1 参数估计方法
引入随机参数后,似然函数式(10)中含有无法直接观测到的未知参数β,难以解析求解未知参数估计值。如果将未知参数β作为隐变量,共轭分布下β的后验分布与先验分布类型一致,此时β的后验期望计算简易可行。因此共轭分布下的相关研究通常采用EM算法[13,19-20]进行参数估计。当β的分布为非共轭分布时,难以解析确定β的后验分布类型,此时β的后验期望计算困难,导致难以采用EM算法估计模型参数。
MCMC算法适用于共轭分布及非共轭分布两种情形下的参数估计[21-23],考虑β可能服从共轭分布或非共轭分布的两种情形,引入MCMC算法获取未知参数的点估计及区间估计。β非共轭分布下,参数向量θ及产品个体随机参数βj的估计步骤如下。
步骤1:选择某个非共轭分布(如Weibull分布)f(βj)作为随机参数βj分布,并将无信息分布π(θ)作为θ的先验分布,设定θ及βj的初始值θ(0)、β(0)j(j=1,2,···,n)。
步骤2:建立退化增量数据Δy的条件似然函数L(Δy|θ,β)=L1,其中β=(β1,···,βn)代表所有个体的随机参数,并且各随机参数βj均服从同一非共轭分布。
步骤3:利用非共轭分布乘积、先验分布与条件似然函数可以计算完全条件后验分布,其中θ\j=(θ1,···,θj-1,θj+1,···,θd)。
步骤4:选取不同的参数初始值组合,产生3条马尔可夫链,每条链可利用MCMC算法从完全条件后验分布中抽取60 000个样本,取这些样本中的后50 000个样本用于参数估计。
步骤5:检验样本序列收敛性,利用Gelman-Rubin方法[18]计算统计量GR,GR为合并后验方差V与样本序列内后验方差平均值WSS之比,有GR=V/WSS,当GR趋近于1时,生成的样本序列是收敛的。
步骤6:利用收敛的样本序列,计算后验均值作为参数的点估计,计算后验分位数获取参数的区间估计。
将上述步骤1~3中的非共轭分布替换为共轭分布后,就可以将上述估计步骤应用于共轭分布下的参数估计。
2.2 参数分布选取准则
不同类型的参数分布会产生不同的寿命预测结果,获取参数估计值后,需要结合实测数据选择合适的参数分布,包括共轭分布(Gamma分布)或非共轭分布(如Weibull分布)。工程实际中通常采用赤池信息准则(Akaike information criterion,AIC)评价参数分布的数据拟合效果,某一分布下AIC数值越小,该分布对数据的拟合效果越好。AIC定义为:
(16)
其中,p为模型未知参数的个数,lnL为对数似然函数。共轭分布下对数似然函数为lnL3,完成参数估计后可以利用式(16)直接计算。非共轭分布下似然函数为式(10),难以积分解析计算。
由式(10)可知,似然函数L2实际是随机参数β的函数f(Δyjk|β)乘积的数学期望。同时考虑共轭分布及非共轭分布的两种情形,采用MC算法计算对数似然函数lnL2。
步骤1:选定β的具体分布类型,利用退化数据,计算模型未知参数θ=(α,γ,δ,ρ)的估计值。
步骤2:从f(β|δ,γ)分布中依次独立抽取的样本值(b=1,2,···,B),抽取次数B可设定为B=5 000。
步骤3:将代入f(Δyjk|β),计算的样本值。
步骤4:利用大样本平均值收敛于数学期望的性质,可得
(17)
将式(17)代入式(16),可以得到不同分布下的AIC值,取最小AIC对应的参数分布,作为参数β的统计分布,并进一步对产品的寿命进行预测。
3 仿真算例
3.1 仿真算例1
GaAs激光器的性能可用工作电流表征,在激光器的使用过程中其工作电流不断增加。将工作电流相对增加量作为产品的退化量,同一时间不同个体的退化量显著不同,即激光器退化量存在明显的个体差异[7]。采用随机参数Gamma过程作为激光器的退化模型,本例中设定随机参数分布为共轭Gamma分布,生成仿真退化数据。
仿真退化试验样本量n=15,测量间隔为250 h,每个样本退化量测量次数m=16。已知该激光器为线性退化轨迹,即已知ρ=1。设定其余退化模型参数分别为α=0.038 8、δ=28.45、 γ=1.455,利用式(6)与式(11)生成仿真退化数据。
假设随机参数服从Gamma分布,利用MCMC算法得到共轭分布下未知参数的点估计与区间估计,参数估计稳定收敛,将估计值与未知参数真值进行比较,结果如表1所示,其中,2.5%和97.5%代表参数的2.5%分位数估计及97.5%分位数估计,可用于计算参数95%置信区间。
同时,为验证MC算法计算似然函数值的准确性,将解析式(15)直接计算的结果与MC算法计算结果进行比较,结果如表2所示。
表1与表2的Gamma-Gamma代表仿真数据的生成与可靠性建模均采用Gamma分布,结果表明,利用MCMC算法所得参数估计值与真值相差较小,并且MC算法计算出的似然函数值与直接计算值接近,验证了所提方法的有效性。
表1Gamma-Gamma下参数估计值
Tab.1 Parameter estimation values with the Gamma-Gamma case
表2Gamma-Gamma下对数似然函数计算值
Tab.2 Logarithmic likelihood function calculation values with the Gamma-Gamma case
假设随机参数服从Weibull分布,参数估计值如表3所示,其中Gamma-Weibull代表仿真数据的生成与可靠性建模分别采用Gamma分布及Weibull分布。比较表1与表3,结果表明不同分布下的同一参数α估计值接近。
表3Gamma-Weibull下参数估计值
Tab.3 Parameter estimation values with the Gamma-Weibull case
进一步算得Gamma-Gamma下的AIC值为-162.01,Gamma-Weibull下的AIC值为-161.41。由于Gamma分布下的AIC值更小,选择Gamma分布作为β的分布。这一选择结果与共轭分布的原始设定相符,验证了分布选取准则的正确性。
3.2 仿真算例2
同样考虑GaAs激光器产品的性能退化仿真,与算例1不同的是,本例中设定随机参数分布为非共轭Weibull分布。设定Weibull分布参数分别为δ=5.947、γ=19.66,其他参数与算例1一致,生成仿真退化数据。
假设随机参数服从Weibull分布,利用MCMC算法估计参数,参数估计值如表4所示,其中Weibull-Weibull代表仿真数据的生成与可靠性建模均采用Weibull分布。此时参数估计值与真值同样相差较小,所有参数的真值均在参数的区间估计范围内。
表4Weibull-Weibull下参数估计值
Tab.4 Parameter estimation values with the Weibull-Weibull case
假设随机参数服从Gamma分布,参数估计值如表5所示。结果再次表明,不同分布下的同一参数α估计值接近。
表5Weibull-Gamma下参数估计值
Tab.5 Parameter estimation values with the Weibull-Gamma case
进一步算得Weibull-Weibull下的AIC值为-162.20,Weibull-Gamma下的AIC值为-160.99。由于Weibull分布下的AIC值更小,选择Weibull分布作为β的分布。这一选择结果与非共轭分布的原始设定相符,再次验证了分布选取准则的正确性。
4 应用实例
4.1 应用实例1
以GaAs激光器实测性能退化数据[24]为应用实例开展可靠性建模及寿命预测。随着激光器性能的不断退化,工作电流不断增大,当工作电流相对增加量达到阈值ω=10时,认为激光器已经失效。试验样本量为15,所有样本退化量均在同一时间节点测量,测试间隔为250 h,实测的性能退化曲线如图1所示。
图1表明该激光器为线性退化轨迹,可认为ρ=1,并且各样本间存在明显的个体差异。分别假设随机参数服从Gamma分布及Weibull分布,利用MCMC算法得到两种分布下参数的估计值。Tsai等[13]利用EM算法得到了Gamma分布下的参数估计值,MCMC算法与EM算法参数估计值的比较如表6所示。
表6表明Gamma分布下MCMC算法得到的参数估计值与EM算法接近,再次验证了所提参数估计方法的有效性。进一步得到两种分布下β的概率密度函数如图2所示。
图1激光器性能退化曲线
Fig.1Degradation curves of laser performance
表6MCMC算法与EM算法参数估计值的比较
Tab.6 Comparison of parameter estimation values with MCMC method and EM method
图2激光器退化过程随机参数的概率密度函数
Fig.2Probability density functions of random parameters for laser degradation process
图2表明不同分布下随机参数的概率密度函数存在差异,由此可知不同分布下产品的寿命分布函数也存在差异。
将参数估计值代入式(8),利用MC算法,计算出两种分布下激光器的寿命分布函数如图3所示。图3表明两种分布下激光器的中位寿命预测会有一定的差异。
进一步算得Gamma分布下的AIC值为-181.78,Weibull分布下的AIC值为-183.0。由于Weibull分布下的AIC值更小,应选择Weibull分布作为激光器退化过程参数β的分布,并利用Weibull分布预测激光器的中位寿命为5 175 h。
图3两种分布下激光器寿命分布函数比较
Fig.3Comparison of laser lifetime distribution functions under two distributions
Tsai等[13]采用共轭Gamma分布作为Gamma过程的随机参数分布,在此基础上预测激光器的寿命分布,并利用伪失效寿命估计值作为参考,验证了带有随机参数可靠性模型下的寿命预测精度优于无随机参数模型下的寿命预测精度。
同样利用最小二乘法拟合各样本退化轨迹,计算得到各样本的伪失效寿命值T如表7所示。
表7激光器各样本伪失效寿命估计值
Tab.7 Pseudo-failure lifetime estimation values for laser samples
利用表7数据可得激光器产品中位伪失效寿命估计值为5 165 h,利用式(14)可得Tsai等[13]所用共轭分布模型下中位寿命预测值为5 045 h,共轭分布及非共轭分布下的寿命预测与伪失效寿命估计比较如表8所示。
由表8可知,利用伪失效寿命值作为参考,可以验证非共轭Weibull分布模型下的寿命预测精度高于Tsai等[13]所用共轭分布模型下的预测精度。
表8伪失效寿命估计与两种寿命预测
Tab.8 Pseudo-failure lifetime estimation and two types of life prediction
4.2 应用实例2
以2024-T351铝合金实测性能退化数据[25]为应用实例开展可靠性建模及寿命预测。随着铝合金疲劳性能的不断退化,疲劳裂纹尺寸不断增长,当疲劳裂纹尺寸增长量达到阈值ω=10 mm时,认为铝合金已经疲劳失效。试验样本量为30,所有样本退化量均在同一时间节点测量,测试节点为1×104 r、1.5×104 r、2×104 r、2.5×104 r、3×104 r、3.5×104 r、4×104 r,实测的性能退化曲线如图4所示。
图4铝合金性能退化曲线
Fig.4Degradation curves of aluminum alloy performance
图4表明铝合金的退化轨迹并不是线性的,ρ≠1为未知参数,各样本间同样存在个体差异。利用MCMC算法得到Gamma分布下与Weibull分布下参数估计值,Weibull分布下参数估计值如表9所示,Gamma分布下情况类似,不再给出。
表9Weibull分布下参数估计值
Tab.9 Parameter estimation values with the Weibull distribution
表9中ρ>1确认了铝合金退化轨迹的非线性,进一步得到两种分布下β的概率密度函数如图5所示。
图5铝合金退化过程随机参数的概率密度函数
Fig.5Probability density functions of random parameters for aluminum alloy degradation process
比较图2与图5可以发现,两类产品不同分布下的概率密度函数均存在差异,且两类产品不同概率密度函数间的差异也不同。
利用MC算法,计算出两种分布下铝合金的寿命分布函数如图6所示。
图6两种分布下铝合金寿命分布函数比较
Fig.6Comparison of aluminium alloy lifetime distribution functions under two distributions
图6表明两种分布下铝合金失效概率小于0.9的分位寿命都会有一定的差异。为验证MC算法计算寿命分布函数的准确性,将解析式(14)直接计算的Gamma分布下寿命结果,与MC算法计算结果进行比较,结果如图6所示,两者结果接近表明了寿命预测的准确性。
进一步算得Gamma分布下的AIC值为-29.15,Weibull分布下的AIC值为-27.63。由于Gamma分布下的AIC值更小,应选择Gamma分布作为铝合金退化过程参数β的分布,并利用Gamma分布预测铝合金的中位寿命为5.83万转。将实例2与实例1的分布选取结果对比,发现不同产品退化过程随机参数的分布类型并不相同,因此需要利用参数分布选取准则选择合适的分布类型。
5 结论
利用随机参数非共轭分布下的可靠性建模方法分析单调退化产品的性能退化数据,可以精确预测产品的使用寿命,具体如下:
1)长寿命产品的退化过程存在个体差异,需要引入随机参数描述个体差异,共轭Gamma分布的假设不适用于所有产品的随机参数,考虑随机参数非共轭分布下的可靠性建模,能够提高基于随机过程的退化模型的适用性。
2)将共轭分布扩展到非共轭分布,可以获得统一的寿命分布函数及似然函数表达式,利用MCMC算法可以较好地估计未知参数,利用MC算法可以精确地计算寿命分布函数与似然函数值。
3)不同随机参数分布下产品的寿命分布存在差异,依据AIC可以选择合适的随机参数分布类型,得到精确的产品可靠性模型及寿命预测结果。