考虑专家信息的威布尔型产品定时截尾可靠性验收试验方案设计

谭 尧, 赵 骞, 王文峰, 郭 波,*, 蒋 平

(1. 国防科技大学系统工程学院, 湖南 长沙 410005; 2. 国防科技大学信息通信学院,陕西 西安 710106; 3. 中国人民解放军93209部队, 北京 100085)

0 引 言

可靠性验收试验(reliability acceptance test, RAT)是检验产品可靠性是否达到要求水平的试验。根据可靠性验证试验的结果可对参试产品做出接收或拒绝的结论。验证试验主要是根据以下假设检验的结果来做出判断:

其中,是参数的可接受值,而是参数的拒绝值。如果根据试验结果计算判断是接受H,则接收该批产品;如果判断是接受,则拒绝该批产品。在检验过程中,当H成立而接受H,此类错误称为第一类错误;当H成立而接受H,此类错误称为第二类错误。发生第一类错误和第二类错误的概率分别称为“生产方风险”和“使用方风险”。只有当这两类风险都是受控的情况下,才能让试验方案被生产方和使用方接受。

在工程应用中,国内外都有标准来指导RAT。MIL-HDBK-781给出了较为完整的可靠性验收方案,为美军进行RAT提供了有力支持。GJB899A—2009《可靠性鉴定与验收试验》以定时截尾试验为例,当生产、接受双方确认验收关键参数,如寿命检验上下限、和生产方、使用方风险、后,可确定试验方案(,)。其中,为累计试验时间,为接收故障数。在此试验中对参加验收试验的产品进行定时截尾试验。当累计试验时间达到后试验结束。如果试验过程中故障数不大于,那么该批产品被接收,否则拒绝接收该批产品。

GJB899A给出的试验方案是基于产品寿命服从指数分布的假设前提下得到的。但在实际应用中,很多产品寿命并不服从指数分布。例如,雷达系统、汽车零部件的寿命往往服从威布尔分布。因此,本文基于威布尔分布,研究RAT设计问题。此外,GJB899A中定时截尾试验方案的试验时间都较长,短时试验方案又存在两类风险较大的问题,实际应用中颇受诟病。

对产品先验信息的利用是通用且有效的。在工程实践中,产品参加验收前,生产方可能会提供一些产品信息,如参数统计信息。Jun等提出寿命服从威布尔分布且形状参数已知条件下抽样验收方案。Tsai等提出了形状参数已知,基于定时截尾缺失数据,并且借助成本模型,完成了可靠性抽样验证方案的设计。Chen等利用贝叶斯决策理论和研制过程中产品寿命信息,制订了基于定时截尾数据的鉴定方案。若产品状态发生变化,验收鉴定方案也会调整。Balamurali等提出多重递延状态重复群抽样计划。

上述工作主要是对产品进行抽样检验以决定验收方案。在实际操作中,由于经费成本的限制,可投入的试验样本量往往较少,传统试验方法在解决小子样问题上存在局限性。目前正在使用的验收试验方案是GJB899A中提出的,但由于只利用了系统级试验数据来开展试验,导致试验时间过长,或者在短时间的试验方案的两类风险较大,实际操作性不强。因此,有学者研究了利用先验信息来开展可靠性试验和评估。专家信息是一类基于专家长期经验的主观信息,同时也是小子样条件下装备可靠性试验中先验信息的重要来源。专家信息的形式包括可靠度点估计、寿命点估计、可靠度置信下限等。明志茂等讨论了利用先验信息制定试验方案的问题。Zhang等讨论了专家评估是文本形式的情况,提出了四步结构的可靠性评估框架。Zhao等将专家信息转化为参数的分布,同时Zhao等还讨论了不同形式专家数据转化为参数分布的方法。杨军等在将专家信息转化为概率分布的同时研究了对其进行融合的问题。

上述研究讨论了如何利用专家信息确定参数的分布,但对如何利用其进行RAT的研究相对较少。为同时解决有效利用先验信息、小子样验收方案设计两个问题,本文提出了基于专家信息的威布尔型产品可靠性验收方案。

1 威布尔分布与参数分布

威布尔分布在可靠性工程中被广泛应用,尤其适用于描述机电类产品寿命分布。因而,本文基于威布尔型产品,研究验收试验方案设计方法。

1.1 威布尔分布

威布尔分布的累积分布函数(cumulative distribution function, CDF)表示为

(1)

式中：为形状参数；为尺度参数。

通常令=1,有

()=1-exp(-)

(2)

则威布尔分布时刻可靠度函数可表示为

()=1-()=exp(-)

(3)

1.2 参数的分布

当前验收试验样本量较小,传统可靠性验收方案由于未能利用先验信息,导致试验时间长、成本高。现利用融合专家信息的方式对传统方案进行改进。专家信息可用来确定参数的分布,在这之前,首先需要明确参数分布的形式。

可假设威布尔分布形状参数的分布服从均匀分布,即

(4)

式中：参数,可通过工程实践经验来确定。

Zhao等从产品可靠度的分布出发,推导寿命分布参数的分布函数,该方法便于操作且取得了较好的效果。根据共轭先验理论,系统在时刻可靠度先验分布可认为服从负对数伽马(negative log Gamma, NLG)分布,具体表示为

(5)

式中：,是NLG分布的两个参数。参数的分布已知,为了得到参数的分布,参考式(3),对关于求导可得

(6)

由式(6)可知,随的增大而减小,故有

(7)

对式(7)关于求导，可得

(8)

因此,参数的分布可表示为

(9)

结合式(4)可得和联合分布为

(10)

式中:∈[0,∞],∈[,]。由此可见,确定分布的核心便是根据先验信息计算分布参数和。

2 基于专家信息的参数分布确定

本节研究如何利用不同类型专家信息确定参数的分布,为制定更加合理的验收方案提供帮助。

2.1 基于可靠度点估计确定参数分布

当专家信息为“产品时刻的可靠度为”时，有

(11)

通过使用最大熵方法,可以将确定参数分布问题转化为非线性规划:

(12)

对式(12)进行化简可得

(13)

2.2 基于寿命点估计确定参数分布

当专家信息为“产品的寿命为”时,根据威布尔分布基本性质可知寿命点估计为

(14)

参考时刻可靠度点估计形式专家信息的处理方法,同样可利用极大熵方法确定参数分布,具体形式为

(15)

将式(10)和式(14)代入式(15),可化简为

(16)

进一步可得

(17)

2.3 基于可靠度置信下限确定参数分布

当专家信息为“产品在时刻置信度为100(1-)%的可靠度置信下限为”时，有

(18)

利用极大熵方法可得

(19)

对式(19)进一步简化为

max=-ln+ln()+-s.t(,-ln)=1-

(20)

其中,

(21)

为下不完全伽马函数。

3 两类风险计算与试验方案的确定

第1节介绍了威布尔分布性质以及参数分布的确定,本节将给出两类风险的计算过程。计算过程利用了第2节由专家信息确定的参数分布,并将其作为确定验收方案的重要依据。

假设在产品验收时,采用了试验时间为的定时截尾验收试验方案。在时间(0,)内发生次故障的概率可根据二项分布计算得到,即

(22)

式中：为参与试验样本数;()，()分别为产品的寿命CDF和可靠度函数。

3.1 两类风险的计算

311 生产方风险的计算

生产方风险表征产品寿命符合标准时产品被拒收的概率。根据生产方风险定义及威布尔分布假设,生产方风险可表示为

(23)

式中：为验收试验发生故障的产品个数；为接收产品的最大故障数；和分别为产品寿命检验的上限和下限。其中,当>时,拒绝接收产品。

结合式(10)及式(14),式(23)分母部分可表示为

(24)

式中：=((1+1))；(,)为下不完全伽马函数,表达式同式(21)。求解分子部分即为求解

(25)

将式(10)及式(22)代入式(25)可得分子部分解析式。

312 使用方风险的计算

使用方风险表征产品寿命不符合标准时产品被接收的概率。使用方风险的计算方式与生产方风险的计算方式相同。根据使用方风险的定义可知:

(26)

结合式(10)及式(14),式(26)分母部分可表示为

(27)

其中,

(28)

为上不完全伽马函数。参考生产方风险计算过程,式(26)分子部分可进一步表示为

(29)

将式(10)及式(22)代入式(29)可得分子解析式。

3.2 基于抽样的两类风险计算

由于上述计算两类风险的解析式比较复杂,本文拟使用基于抽样的算法,确定参数,的分布(,)以及仿真次数,该算法具体步骤如下所示。

算法 1 基于抽样的两类风险计算方法1. 根据先验分布π(λ,m)抽出容量为S的样本组(λi,mi)。2. 将步骤1中的样本代入式(14),根据判断条件θ>θ0,θ<θ1选出符合条件的两个样本组。3. 将两个样本组代入式(22)并分别求均值,可得式(23)和式(26),即生产方、接收方风险分子的估计值[19]。4. 将步骤2中筛选出的符合判断条件θ>θ0,θ<θ1的样本组分别除以总样本数,可得式(24)和式(27),即生产方、接收方风险分母的估计值。5. 根据步骤3和步骤4求得结果可计算两类风险值。

显然,算法1的核心问题在于如何从参数的分布中抽取随机样本。

参考式(10),从先验分布(,)中抽样得到随机样本组(,),具体过程如下:

(,)∝(|)()

(30)

式中：()参考式(4)；(|)参考式(9),同时

(31)

由式(31)可知,(|)服从伽马分布。故在算法1步骤1求得(,)后,可以从式(4)中随机抽样得到样本。将代入式(31)中可以对进行抽样,这样就得到了一组样本(,)。

3.3 试验方案(T,n,C)的确定

在GJB899A中,产品验收方案由两类风险,以及寿命检验上限、下限4个参数共同确定。本文基本思路与GJB899A相似,产品验收试验方式采取定时截尾试验。同时,由于所接受产品工作时长小于是不可接受的,故将截尾时间设为寿命检验下限。检验上限为产品寿命期望值。具体步骤参考算法2。

算法 2 产品验收试验方案(t0,n,C)的确定1. 根据产品设计参数及验收方案鉴别比确定寿命检验上下限θ0,θ1及截尾试验时间t0。2. 确定生产方提供的参与试验样本数n。3. C从1开始取值,逐一递增,到n为止。结合专家信息计算并记录C取不同值时的两类风险值。(实现过程可参考第4.1节,图1展示了仿真结果)4. 生产方、接收方协商选择合适的验收方案并实施,若故障数超过C则拒收产品。

4 案例分析

案例中产品验收试验采用定时截尾试验方案。假设产品寿命分布参数为=800,=15。检验上限为产品寿命期望值722 h。以GJB899A中方案12为例,其鉴别比为2,故检验下限为361 h,同时截尾时间为361 h。假设得到的专家信息为产品工作100 h时的可靠度=0.956 8。根据文献[20],对威布尔分布型产品进行定时截尾试验,若截尾时间较短(<05),可用指数分布代替威布尔分布。又361<05,故本文试验方案可与GJB899A试验方案进行对比。仿真试验步骤可参考算法2。

4.1 给定参加试验产品数n确定验收方案

当截尾时间及试验产品数确定后,两类风险值随的变化如图1所示。

图1 鉴别比d=2,测试产品数n=20,验收方案两类 风险值随故障数变化情况Fig.1 Variation of two types of risks given different settings of fault numbers with d=2, n=20

由图1可以发现,当投入试验产品数不变时,对两类风险计算的影响很大。同时,当取8时两类风险值最为接近。但是具体选择拒接收数可以根据实际情况由生产方、使用方专家共同决定。如要更加倾向于降低使用方风险,可选=7的方案。

4.2 给定拒绝接受产品故障数C确定验收方案

若确定,两类风险值会随参加验收试验产品数的变化而变化,如图2所示。假设=5,可以看出两类风险值在测试产品数=9时最为接近。

图2 鉴别比d=2,拒收数C=5,验收方案两类风险值随参与 试验产品数变化情况Fig.2 Variation of two types of risks given different settings of component numbers with d=2, C=5

取鉴别比=2,生产方、接收方商定要求两类风险,<01的情况下,在投入试验产品数取不同值时,可通过算法2依次进行仿真试验。表1展示了当试验产品数取4～20情况下的最佳试验方案。

表1 鉴别比d=2，两类风险α, β<0.1时RAT方案Table 1 Plan of the RAT when discrimination ratio d=2 and two types of risks α, β<0.1

表1中代表总试验时间。由表1可知,当融合专家信息的验收方案在试验时间与GJB899A对比方案接近时,可以明显降低两类风险。同时,在本文提出的验收方案下,两类风险值之和随试验样本数和总试验时间的减少而增加。此时,需要生产方及接收方根据实际情况决定采取何种方案:如果倾向于降低风险,可以考虑提高试验产品数,降低两类风险;如果倾向于控制成本,可以考虑在两类风险不超过要求的前提下减少试验样本数。

对于鉴别比=3,生产方、接收方商定要求两类风险,<0.1时,类比上述计算步骤,可以得到表2。由表2结果可知,融合专家信息的产品验收试验方案显著降低了两类风险及试验时间,进一步验证了所提方案的有效性。

表2 鉴别比d=3，两类风险α,β<0.1时RAT方案Table 2 Plan of the RAT when discrimination ratio d=3 and two types of risks α, β<0.1

4.3 专家信息与产品实际状态有出入时对验收方案设计的影响及解决

在工程应用中,专家信息由于存在一定的主观性,并不绝对可靠。本节讨论了当专家信息与产品实际情况有出入时对验收方案设计的影响。

4.3.1 专家信息变化对方案的影响

在第4节案例分析中,给出的专家信息为产品工作100 h时的可靠度=0.956 8。此为结合式(3)仿真计算产生的精确专家信息。但实际专家可能难以给出如此精确的数据。为验证专家信息与实际产品状态有偏差时对可靠性验收方案的影响,本节进行对比试验,分别计算专家信息为产品工作100 h时=097,09,085,08,07情况对可靠性验收方案的影响情况。

根据图1所示,当投入样本数=20,鉴别比=2时。若专家信息为=0956 8,验收试验方案选择标准为生产方、使用方风险尽量接近,使用方风险略小,那么应该选择=8的验收方案。在保持=20,=2不变的情况下,对上述提到不同可靠度取值对应的两类风险进行计算,结果如表3所示。

表3 d=2,n=20时,验收方案随专家信息变化情况Table 3 Variation of RAT with the changes of expert information while d=2, n=20

在表3中,当专家信息为=0.956 8,故障数取8～11时，计算所得两类风险最为准确,且风险都在可控范围内,可作为参考。其中，=8由于满足验收方案选择条件,是最合适的方案。当专家信息为=097,09,08,07,06时,根据两类风险计算结果,取9,10,偏离实际最优方案=8,但依然可接受。当专家信息为=055时取11,此时虽未超过<01,<01的限制,但生产方、使用方风险较悬殊,不宜采纳。

由上述结果可知,当专家信息与实际情况偏差在一定范围内(Δ=0356 8)时,对验收方案的设计有影响,但影响在可接受范围内。这也间接证明了本文所提出方法的稳定性。

4.3.2 利用融合的先验分布优化验收方案

专家信息是根据专家长期经验进行总结而得到的信息。可能会出现个别专家信息置信度不高,或与实际出入过大的情况。第4.3.1节证明了专家信息存在偏差时,验收方案两类风险值依然可以保持在限定范围内。若个别专家信息偏差较大,对多个专家信息进行融合较为合理。这样可以降低置信度较低的专家信息带来的影响。融合的先验分布形式为

(32)

式中：为通过第2节方法确定的先验分布；为其对应的权重。由此可知融合的核心为确定先验分布权重。具体方法可参考于春雨等应用基于D-S证据推理理论的权重求解方法。

以生产方风险为例,通过推导可知融合的先验分布条件下有

(33)

式中：为先验分布通过第311节所提出的算法计算的生产方风险,使用方风险同理。故融合先验分布两类风险计算相当于各先验分布独立计算两类风险后再进行加权。在利用第3.2节提出的算法1进行两类风险计算后,结合式(33)可得融合的两类风险。利用第3.3节提出的试验方案确定法可得到融合多专家信息的试验方案。

5 结 论

由于GJB899A的定时截尾试验方案的制订仅考虑产品寿命分布类型(如指数分布)信息,导致验收试验时间长、风险大,已经难以适用于大型复杂装备的可靠性鉴定验收。本文基于不同类型专家信息,对产品RAT方法进行了改进。通过将本文提出的验收方案与GJB899A中的方案对比,可以发现有效利用专家信息可以显著降低两类风险值、缩短试验时间。案例分析中列举的两种试验条件下进行的试验方案充分证明了所提出理论的有效性。

下一步准备在现有基础上继续深入,针对不同分布类型的产品,如正态分布型产品等,开展试验鉴定方案制定方法研究。

此外,为降低单个专家信息置信度不高所带来的不确定性,后续研究将对第4.3.2节进行更加细致的研究,提出更准确、更合理的融合多个专家信息的可靠性验收方案制订法,充分挖掘专家信息，以达到最大化利用的目的。