简单步加寿命试验在定时截尾逐级删失模型下的优化设计
杨春燕, 姜万民
云南大学 数学与统计学院,云南 昆明 650091

作者简介:杨春燕(1971-),女,云南人,博士,教授,主要从事概率论与数理统计方面的研究.E-mail:yangchy@ynu.edu.cn.

摘要

研究了简单步进应力加速寿命试验在带有随机移走的定时截尾逐级删失模型下应力转换时间的优化问题.假定产品寿命服从Weibull分布,推导了似然函数和Fisher信息阵,得到了产品平均寿命估计量的渐进方差的表达式,利用Monte-Carlo方法对不同程度的移走和不同的失效模型进行模拟,通过极小化渐进方差得到了双应力步进寿命试验中应力改变时间的最优值.研究结果对实验者进行该类实验具有指导意义.

关键词: 步进应力加速寿命试验; Weibull分布; 随机移走; 定时截尾; 逐级删失
中图分类号:O213.2 文献标志码:A 文章编号:0258-7971(2018)06-1063-08
Optimal design of simple step stress ALT under type I progressive censoring
YANG Chun-yan, JIANG Wan-min
School of Mathematics and Statistics,Yunnan University,Kunming 650091,China
Abstract

We investigate the optimal plans of step stress accelerated life tests (ALT) under type I progressive censoring with random removals.Assume that the lifetime of the product follows a Weibull distribution,the maximum likelihood function and the Fisher information matrix are derived.Then the expression of asymptotic variance of the estimator of the mean lifetime is obtained.Monte-Carlo simulation is conducted to find the optimal time changing points which minimizes the asymptotic variance for different removal probabilities and lifetime parameters.The results of this study could provide the guidance for the experimenter to conduct the test.

Keyword: step stress accelerated life test; Weibull distribution; random removals; type I censoring; progressive censoring

可靠性是评价产品质量的一个重要指标, 它反映了产品在规定的时间和条件(温度、湿度、压强等)下能够正常工作的可能性.为评估产品的可靠性, 寿命试验通常被用来收集数据, 之后通过统计推断的方法得出产品寿命分布的规律性, 从而对产品的质量作出评价.随着科学技术的不断提高, 高可靠、长寿命的产品越来越多, 通过普通寿命试验很难收集到数据.加速寿命试验应运而生, 并且自出现以来就得到了快速发展.加速寿命试验是将产品置于高应力水平(高温、高湿、高压等)下进行试验, 以促使产品尽快失效, 从而缩短试验周期, 提高试验效率.对于加速寿命试验的统计分析, 文献中有详细的讨论和研究[1, 2, 3].

根据应力变化方式的不同, 加速寿命试验可分为3种:恒定应力加速(恒加)、步进应力加速(步加)和序进应力加速(序加).其中, 对恒加和步加试验的研究比较成熟, 而步加试验相对于恒加试验又具有样品失效较快、所需样本量较小的优点, 因此在工程领域中被广泛使用.

在试验过程中, 由于经费的削减、其它试验的需要或者试验温度过高等各种原因, 通常会在试验的各个阶段将还未失效的产品移走, 这就导致了逐级删失的发生.关于逐级删失模型的讨论可参见文献[4].逐级删失情况下步加试验的优化设计在国外已有很多研究, 而在国内几乎还是空白.假定产品寿命服从指数分布, 文献[5, 6, 7]研究了逐级删失定时截尾情况下的优化设计, 文献[8]研究了逐级删失定数截尾情况下的优化设计.在这些研究中均假定在试验进行的每个阶段被移走的产品数目是已知的, 在进行试验之前就被确定.但事实上, 由产生逐级删失的原因可知被移走的产品数目只能是随机的, 在试验之前不可能预先确定.文献[9]首次对带有随机移走的逐级删失模型在普通寿命试验条件下的参数估计问题进行了研究, 之后该模型在普通寿命试验和加速寿命试验中引起了广泛讨论.其中, 文献[10, 11]分别假设产品寿命服从指数分布和对数正态分布, 对恒加试验进行了优化设计; 文献[12, 13]假定产品寿命服从指数分布, 从不同的角度考虑了带有随机移走的定时截尾步加试验的优化设计.文献[14, 15]分别讨论了产品寿命服从Frechet分布和广义Rayleigh分布时带有随机移走的定数截尾步加试验的优化设计.文献中对步加寿命试验优化设计的讨论有很多, 最新的研究可参看文献.

Weibull分布是可靠性研究中最常用的一种分布, 用途非常广泛.本文假定产品的寿命服从Weibull分布, 在带有随机移走的定时截尾逐级删失模型下, 对简单步加试验中应力的转换时间进行了优化设计, 并讨论了试验参数的改变对最优转换时间产生的影响, 研究结果可为试验者进行该类型的寿命试验提供有价值的参考.

1 模型的建立

本文中使用的符号:G(· ):CE模型累积函数; α , β :Weibull分布的形状和尺度参数; F* (· ):残存函数(生存函数); xi(i=0, 1, 2):步进应力水平; p:随机移走的概率; τ :应力转换时间; T:定时截尾时间.

基本假设:

(1) 产品在任何应力水平下都是服从Weibull分布的;

(2) 形状参数α 是一个常数;

(3) 尺度参数β 相对于应力水平xi满足lnβ i=a+bxi(b< 0, i=0, 1, 2), 其中a和b是由产品自身属性决定的未知参数, x0是一般应力条件;

(4) 样品之间是相互独立的;

(5) 在双应力试验下产品的分布符合累计失效模型(CE模型).

试验过程:在这里我们只考虑2个应力水平的情形.将n个试验样品放在应力水平x1下, 直到转换时间τ 为止, 观察到的失效样品个数记为n1, 失效时间记为t1j(j=1, 2, …, n1); 此时随机移走r个样品, 接着应力水平变为x2, 一直持续到时间T为止, 这时观测到的样品失效个数记为n2, 失效时间记为t2j(j=1, 2, …, n2), 剩余个数为nc=n-n1-r-n2.

根据假设, 低应力水平x1和高应力水平x2下产品的寿命分布分别为

F1(t)=1-exp -tβ1α和F2(t)=1-exp -tβ2α,

相对应的步进应力试验下的累积寿命函数为

G(y)= F1(y), 0yτ, F2(y-τ+t), y> τ,

由累计失效模型原理我们可知F1(τ )=F2(t), 其中t= β2β1τ 为解.相对应的密度函数为

f(y)= f1(y)=αβ1yβ1α-1exp-yβ1α, 0yτ, f2(y)=αβ2y-τβ2+τβ1α-1exp-y-τβ2+τβ1α, τyT.(1)

试验中时间τ 之前失效产品数记为n1, 在到达τ 的时候移走产品数目记为R, 从时间τ 到截尾时间T之间的失效数目记为n2.2次失效产品所有的观察值为tij(i=1, 2; j=1, 2, …, ni), 在给定移走数目R=r的条件下, 可得似然函数:

L(β 1, β 2|R=r)= j=1n1αβ1t1jβ1α-1exp -t1jβ1αj=1n2αβ2t2j-τβ2+τβ1α-1·

exp -t2j-τβ2+τβ1αexp -rτβ1α-ncT-τβ2+τβ1α. (2)

假设移走的产品数目服从二项分布, 且每一个产品被移走的概率相同, 都为p, 即R~binomial(n-n1, p), 则在n1个产品失效后, 从剩余未失效的n-n1个产品中随机移走r个产品的概率P(R=r)可表示为:

P(R=r)= (n-n1)!r!(n-n1-r)pr(1-p )n-n1-r,

其中r是介于0到n-n1间的整数.

因此, 考虑整个试验过程, 得到似然函数的表达式为:

L(β 1, β 2, p)=L(β 1, β 2|R=r)P(R=r),

具体表示为:

L(β 1, β 2, p)=C j=1n1αβ1t1jβ1α-1exp -t1jβ1αj=1n2αβ2t2j-τβ2+τβ1α-1·

exp -t2j-τβ2+τβ1αexp -rτβ1α-ncT-τβ2+τβ1α·

pr(1-p )n-r-n1, (3)

其中:C= (n-n1)!r!(n-r-n1)!.

lnβ i=a+bxi代入(3)式整理, 可得关于参数a, b和p的似然函数, 取对数可得:

lnL(a, b, p)=lnC+(n1+n2)lnα -n1α (a+bx1)-n2α (a+bx2)+(α -1) j=1n1lnt1j+

(α -1) j=1n2ln t2j-τexp(a+bx2)+τexp(a+bx1)- j=1n1(t1jexp(-a-bx1))α -

α exp(-(a+bx1)α )- j=1n2((t2j-τ )exp(-a-bx2)+exp(-a-bx1)τ )α -

nc((T-τ )exp(-a-bx2)+exp(-a-bx1)τ )α +rlnp+(n-r-n1)ln(1-p). (4)

通过对数似然函数分别对a, b和p求偏导, 可得:

lnLa=n2-n1α -2n2α +α j=1n1(t1jexp(-a-bx1))α +a j=1n2((t2j-τ )exp(-a-bx2)+

exp(-a-bx1)τ )α +α nc((T-τ )exp(-a-bx2)+exp(-a-bx1)τ )α +α rτ α exp(-(a+bx1)α ),

lnLb=n2-n1x1α -2n2x2α +α x1 j=1n1(t1jexp(-a-bx1))α +α x1 j=1n2((t2j-τ )exp(-a-bx2)+

exp(-a-bx1)τ )α +x1nc((T-τ )exp(-a-bx2)+exp(-a-bx1)τ )α +α x1α exp(-(a+bx1)α ),

lnLp= rp+ n-r-n11-p,

令以上3式均等于零可得到似然方程组, 利用Newton-Raphson方法, 可得到参数a, b和p的极大似然估计的数值解.

2 Fisher信息阵和渐进方差

以一般应力水平x0下对数平均寿命极大似然估计(MLE)的渐进方差(AV)达到最小作为优化准则寻找最优转换时间τ * .首先推导渐进方差的表达式.

在Weibull分布场合下给定α , 一般应力水平x0下平均寿命θ 的极大似然估计 θ˙可以表示为:

θ˙= β˙Γ 1+1α1+1αexp( a˙+ b˙x0), (5)

由于α 是已知的, 则(5)式是关于参数ab的函数, 我们可以通过Fisher信息矩阵的逆求出平均寿命估计量的渐进方差的表达式.

关于变量a, bp的Fisher信息矩阵可以表示为:

F=E -2lnLa2-2lnLab0-2lnLba-2lnLb2000-2lnLp2=n A11A120A21A22000A33, (6)

利用似然函数分别对a, bp求二阶偏导数, 并分别令:

B1j=exp(-a-bx1)t1j, B2j=exp(-a-bx2)(t2j-τ ),

C1=exp(-a-bx1)τ , C2=exp(-a-bx2)(T-τ ).

整理可得:

A11= 1nE j=1n1(-α2B1jα)+ 1nE j=1n2-α2(B2j+C1)α+E ncn[-α 2(C2+C1)]+E rn(-α 2 C1α),

A12=A21= 1nE j=1n1(-α2x1B1jα)+ 1nE j=1n2-α2(B2j+C1)α(x2B2j+x1C1)+

E ncn[-α 2(C2+C1)α -1(x2C2+x1C1)]+E rn(-α 2x1 C1α),

A22= 1nE{ j=1n2{-(α -1)( x22B2jα+ x12C1)(B2j+C1)-1+(α -1)(x2B2j+x1C1)2(B2j+C1)-2-

α (α -1)(B2j+C1)α -2(x2B2j+x1C1)2-α (B2j+C1)α -1( x22B2j+ x12C1)}}+

E ncn[-α (α -1)(C2+C1)α -2(x2C2+x1C1)2-α (C2+C1)α -1( x22C2+ x12C1)]+

1nE j=1n1(-α2x12B1jα)+E rn(-α 2 x12C1α),

A33=E rnp2+n-r-n1n(1-p)2.

为计算以上各项期望, 记A, B, C分别为相应公式中的指数或者系数.根据试验过程知道, 相同样品在第一个应力水平下失效的个数服从二项分布B(n1, F1(τ )), 因此有:

E j=1n1CB1jα=n 0τCtα f(t)dt.

同理我们可以得到:

E j=1n2[C(B2j+C1)A(x2B2j+x1C1)B]=n 0τC((t-τ )+C1)A(x2(t-τ )+x1C1)Bf2(t)dt,

E rn=p τTF1(t)dt=p(1-F1(τ )),

E ncn=(1-p)(1-F1(τ ))- τTf2(t)dt,

E rnp2+n-r-n1n(1-p)2= 1p(1-F1(τ ))+ 11-p(1-F1(τ )).

F1=n A11 A12A21 A22, 在一般应力水平下, 对数平均寿命lnθ 的极大似然估计ln θ˙的渐进方差可以表示为:AV(ln θ˙)=B F1-1BT, 其中B= lnθ˙a, lnθ˙b.

3 主要数值结果及分析

使渐进方差达到最小的转换时间τ * 无法根据表达式得到解析式, 所以我们利用数值方法来计算最优转换时间.按照试验方法, 假设α 已知, 给定应力水平x0, x1, x2和试验总截尾时间T, 对给定的试验样品数n, 移走概率p以及参数ab, 用Monte-Carlo方法产生来自Weibull分布的容量为n的一个样本, 按照试验方法进行试验, 得到试验观测值n1, n2, r, tij, i=1, 2; j=1, 2, …, nj.利用似然方程和Newton-Raphson方法得到模型参数ab的估计值 a˙, b˙.将所有给定值、观测值和估计值带入渐进方差表达式就可以得到关于应力转换时间τ 的函数, 由于时间τ 所属的时间区间为(0, T), 将此区间作满足一定精度的分割, 并将每个分割点的值代入渐进方差的表达式, 即可得到不同的渐进方差的值, 对应最小渐进方差的分隔点τ 即为最优解τ * .

α =0.5, 1, 2, 分别对应前期失效、偶发失效和磨损失效3种不同的失效模式; 应力水平x0=1, x1=1.2, 1.4, 1.6, 1.8, x2=2, 2.1, 2.3, 2.5, 移走概率p=0.01, 0.05, 0.1, 0.2, 0.3, 0.5, 0.7分别代表不同程度的移走.给定模型参数a=10和b=-1, 产生一个n=50的来自Weibull分布的样本, 得到极大似然估计的值 a˙=9.518645, b˙=-1.080143.不失一般性, 设试验的总截尾时间T=1, 将τ 所在区间(0, 1)以0.001为单位进行分割, 可以获取999个格点.将已知的参数值、观测值代入渐进方差表达式得到关于τ 的函数, 再将999个分割点的数值代入可得999个渐进方差的值, 这些值中的最小值对应的τ 即为优化值τ * .对各参数的不同取值进行模拟, 经过计算, 我们可以得到如下结果(见表1~3).

表1 α =0.5时的最优解τ * Tab.1 The optimal plans τ * for α =0.5
表2 α =1时的最优解τ * Tab.2 The optimal plans τ * for α =1
表3 α =2时的最优解τ * Tab.3 The optimal plans τ * for α =2

由上面的数值结果发现:

(1) 当移走概率增大的时候, 转换时间的最优值会随之减小;

(2) 固定应力水平x1, 转换时间τ 的优化值τ * 随着第二个应力水平x2的增大而增大;

(3) 固定应力水平x2, 转换时间τ 的优化值τ * 会随着第一个应力水平x1的增大而减小.

此外, 比较表1~3中形状参数α 对优化值的影响发现:

(1) α =0.5(前期失效模式), 此时转换时间τ 的优化值τ * 较大;

(2) α =1(偶发失效模式), 此时转换时间τ 的优化值τ * 相对α =0.5时的结果来说较小;

(3) α =2(磨损失效模式), 此时转换时间τ 的优化值τ * 更小.

当形状参数α 增大的时候, 转换时间τ 的优化值τ * 会随之减小, 而渐进方差的值却会相应增大.

4 结 论

假定产品寿命服从Weibull分布, 本文研究了双应力步进加速寿命试验在带有随机移走的定时截尾逐级删失模型下应力转换时间的优化问题.对给定的应力水平, 通过数值模拟给出了3种失效模型在不同程度的移走情况下应力转换时间的最优值, 为实验者进行该类试验提供了有价值的参考.

致谢 该研究的部分工作于作者访问英国肯特大学(University of Kent, UK)期间完成, 在此对肯特大学以及导师张健教授提供的帮助表示衷心感谢!

The authors have declared that no competing interests exist.

参考文献
[1] BAGDONAVICIUS V, NIKULIN M. Accelerated life models: modelling and statistical analysis[M]. Boca Raton, FL: Chapman & Hall, 2002. [本文引用:1]
[2] MEEKER W Q, ESCOBAR L A. Statistical methods for reliability data[M]. New York: Wiley, 1998. [本文引用:1]
[3] NELSON W. Accelerated testing: statistical models, test plans, and data analysis[M]. New York: Wiley, 1990. [本文引用:1]
[4] BALAKRISHNAN N, AGGARWALA R. Progressive censoring: theory, methods and applications[M]. Boston: Birkhauser, 2000. [本文引用:1]
[5] GOUNO E, SEN A, BALAKRISHNAN N. Optimal step-stress test under progressive type-I censoring[J]. IEEE Transactions on Reliability, 2004, 53(3): 388-393. [本文引用:1]
[6] WU S J, LIN Y P, CHEN Y J. Planning step-stress life test with progressively type I group-censored exponential data[J]. Statistica Neerland ica, 2006, 60(1): 46-56. [本文引用:1]
[7] BALAKRISHNAN N, HAN D. Optimal step-stress testing for progressively Type-I censored data from exponential distribution[J]. Journal of Statistical Planning and Inference, 2009, 139(5): 1782-1798. [本文引用:1]
[8] ISMAIL A, SARHAN A M. Optimal design of step-stress life test with progressively type-II censored exponential data[J]. International Mathematical Forum, 2009, 40: 1963-1976. [本文引用:1]
[9] YUEN H K, TSE S K. Parameters estimation for Weibull distributed lifetime under progressive censoring with rand om removals[J]. Journal of Statistical Computation and Simulation, 1996, 55(1): 57-71. [本文引用:1]
[10] YANG C Y, TSE S K. Planning accelerated life tests under progressive type I interval censoring with rand om removals[J]. Communications in Statistics: Simulation and Computation, 2005, 34(4): 1001-1025. [本文引用:1]
[11] 樊爱霞, 杨春燕. 对数正态分布恒加寿命试验在带有随机移走的定时截尾模型下的优化设计[J]. 云南大学学报: 自然科学版, 2009, 31(5): 444-448.
FAN A X, YANG C Y. The optimal design of accelerated life tests for lognormal distribution under type-I interval censoring with progressively rand om removals[J]. Journal of Yunnan University: Natural Sciences Edition, 2009, 31(5): 444-448. [本文引用:1]
[12] WU S J, LIN Y P, CHEN S T. Optimal step-stress test under type I progressive group-censoring with rand om removals[J]. Journal of Statistical Planning and Inference, 2008, 138(4): 817-826. [本文引用:1]
[13] SHEN K F, SHEN Y J, LEU L Y. Design of optimal step-stress accelerated life tests under progressive type I censoring with rand om removals[J]. Quality and Quantity, 2011, 45(3): 587-597. [本文引用:1]
[14] SHAHAB S. Optimal design of step stress partially accelerated life test under progressive type-II censored data with rand om removal for Frechet distribution[J]. Reliability: Theory & Applications, 2015, 10(2): 12-26. [本文引用:1]
[15] ISMAIL A. Planning step-stress life tests for the generalized Rayleigh distribution under progressive type-II censoring with binomial removals[J]. Strength of Materials, 2017, 49(2): 292-306. [本文引用:1]
[16] HAN D. Time and cost constrained optimal designs of constant-stress and step-stress accelerated life tests[J]. Reliability Engineering and System Safety, 2015, 140: 1-14. [本文引用:1]
[17] XU X J. Constrained optimal designs for step-stress accelerated life testing experiments[C]. IEEE 3rd International Conference on Cloud Computing and Big Data Analysis, 2018: 633-638. [本文引用:1]
[18] HAMDY M S. Parameter estimations and optimal design of simple step-stress model for Gamma dual Weibull distribution[J]. Pakistan Journal of Statistics and Operation Research, 2018, 14(1): 109-120. [本文引用:1]