[发明专利]基于贝叶斯理论的定时截尾加速验收抽样试验优化设计方法

权利要求书

1.一种基于贝叶斯理论的定时截尾加速验收抽样试验优化设计方法，其特征在于，包括以下几个步骤：

步骤一、确定产品的寿命分布，寿命验证指标参数，然后建立统计假设；

产品的寿命决定于设计与制造中对其功能、结构、原材料等的选择及质量控制过程中各种随机因素的影响。它是一个服从一定统计规律的随机变量，一般用寿命的分布函数(也称累积分布函数)来描述。

多数产品寿命服从连续型随机变量的概率分布，常用的有指数分布、正态分布、威布尔分布等。某些产品以工作次数、循环周期数作为其寿命度量单位，如开关的开关次数等，这时可以用离散型随机变量的概率分布来描述其寿命分布的规律，如二项分布、泊松分布等。

例如，在可靠性理论中，指数分布是最基本、最常用的分布，电子产品的寿命和复杂系统的故障时间均可用指数分布来描述。指数分布寿命的失效率为常数。很多电子产品在早期失效之后及损耗故障期之前，产品的失效率基本上是稳定的。

指数分布的密度函数有两种表达形式：

f(t)=λexp(-λt)与f(t)=1θexp(-t/θ)---(1)]]>

式中，λ为指数分布的失效率；θ为指数分布的平均寿命。

两个表达式实质相同，在使用条件下，参数λ_U和θ_U(产品使用条件下的失效率和平均寿命)的关系为λ_U=1/θ_U，相应的指数分布函数的形式为：

F(t)=1-exp(-t/θ)与F(t)=1-exp(-λt)                     (2)

假设试验的加速因子为AF，则根据指数分布场合，其加速应力条件下对使用条件下的加速因子定义为：

AF=θUθA=λAλU---(3)]]>

为了方便先验分布的选取，本发明选取失效率λ_A作为其寿命验证指标参数，根据协定的双方风险，在产品的抽样特性曲线(OC Curve)上选择对应的检验上下限λ₀(=1/θ₀)和(λ₁=1/θ₁)，然后建立统计假设如下：

H₀：λ_A≤λ₀·AF H₁：λ_A>λ₁·AF                     (4)

双方商定，当产品批的加速应力下失效率λ_A≤λ₀·AF时，以大概率接收这批产品，即规定生产方风险为α，则接收概率L(λ)≥1-α；当产品批的加速应力下失效率λ_A ≥λ₁·AF时，以小概率接收(高概率拒收)这批产品，即规定使用方风险为β，则接收概率L(λ)≤β。

步骤二、对选定的产品设计试验决策法则，确定其试验和接收流程；

确定了产品的分布及统计假设后，需要确定定时截尾方案的抽验规则，定时截尾试验是指对n个样品进行试验，事先规定试验截尾t₀，到了时刻t₀所有试验样品停止试验，利用试验数据评估产品的可靠性特征量。按试验过程中对发生故障的产品所采取的措施，又可分为有替换和无替换两种方案。有替换指的是试验中某产品发生了故障，立即用新产品代替，保持整个实验过程中样本数不变，而无替换是指当产品发生故障就立刻撤去，在试验过程中，随着故障产品的增加而使样本减少，在本发明中使用的是无替换定时截尾试验。

以指数分布型产品的定时截尾方案为例，方案通常记为(c，T)，其中T为截尾试验时间，c为接收拒绝数。在本发明中，我们使用恒定应力加速试验(CSALT)，寿命试验的决策法则即为：

（1）在产品批中选择n个样品进行CSALT，加速因子为AF，试验为无替换试验；

（2）试验进行到试验累计时间到达预定值T时停止试验，记录试验过程中的失效数；

（3）设在试验过程中出现了r故障，如果r≤c，则认为产品批合格，接收；如果r>c，则认为产品不合格，拒收批产品。

因此，定时截尾试验设计的主要任务就是选择合适的c和T。

步骤三、基于验后风险准则推导满足双方风险关于验证指标参数的约束条件；

本发明将验后风险准则应用在试验方案的求解过程中，验后风险准则借助于现场的试验数据，反向推导指标参数的先验分布中与现场试验得出的寿命水平不同的概率，其计算方法侧重对于参数的先验分布主观认可。

下面介绍验后风险准则(Posterior Risk Criteria)下的双方风险计算原理。

基于概率论与数理统计的原理，基于验后风险准则原理，弃真风险(生产方风险)α的计算公式如式(5)所示：

α=P(R∈θ0|Z∈D1)=∫θ0p(R|Z∈D1)dR=∫θ0P(Z∈D1|R)π(R)dR∫θP(Z∈D1|R)π(R)dR---(5)]]>

式(5)中,θ表示产品寿命参数R的全集，即取值范围p(R|Z∈D₁)表示在给定Z∈D₁的条件下R的概率密度函数，π(R)是寿命参数R的先验分布。

由(5)式所见，弃真风险(生产方风险)α的物理含义为：根据决策法则，依据抽样结果拒绝原假设的情况下，但是，产品总体的寿命却是满足要求的。其数学含义为：在拒绝原假设的前提下，产品寿命参数的验后分布满足要求的概率。

在验后风险准则中，采伪风险(使用方风险)β的定义如式(6)所示：

β=P(R∈θ1|Z∈D0)=∫θ1p(R|Z∈D0)dR=∫θ1P(Z∈D0|R)π(R)dR∫θP(Z∈D0|R)π(R)dR---(6)]]>

式(6)中，p(R|Z∈D₀)表示在现场试验数据得出Z∈D₀基础上的R的概率密度函数。

由式(6)可见，采伪风险β的物理含义为：根据决策法则，依据抽样结果接受原假设的情况下，但是，产品总体的寿命却是不满足要求的。其数学含义为：在接受原假设的前提下，产品寿命参数的验后分布不满足要求的概率。

以指数分布型产品为例，首先推导指数分布定时截尾下产品的接收概率L(λ)。

根据指数分布的累积分布函数F(t)=1-exp(-t/θ)可知,产品的可靠度R(t)=exp(-λt)，到时间t时，n个产品中出现r个故障的概率为

CnrF(t)rR(t)n-r---(7)]]>

到时间t时，产品的故障率r≤c，从而产品被接受的概率为

L(λ)=Σr=0cCnrF(t)rR(t)n-r---(8)]]>

由于λ的值一般都很小，故将R(t)=exp(-λt)泰勒展开可得

R(t)=e-λt=1-λt+12!λ2t2-···≈1-λt]]>

F(t)=1-R(t)=λt            (9)

即可得接受概率

L(λ)=Σr=0cCnr(λt)r(1-λt)n-r---(10)]]>

在nλt≤5，F(t)≤10%的条件下，二项概率可用泊松概率近似，于是得到：

L(λ)=Σr=0ce-nλt(nλt)rr!---(11)]]>

一般情况下n都较小，故T≈nt，从而

L(λ)=Σr=0ce-λT(λT)rr!---(12)]]>

对于验证指标失效率λ，根据Bayes理论，取其共轭先验分布为Gamma分布，记为Gamma(a,b),即：

π(λ|a,b)=baΓ(a)λa-1e-bλ---(13)]]>

其中，Γ(a)为Gamma函数，其定义为：

Γ(a)=∫0∞e-xxa-1dx---(14)]]>

根据验后风险准则和对指数型产品先验分布和接受概率的表达式的推导以及加速条件λ_A=λ_U·AF，对于生产方后验风险α(c,T)的计算公式为：

α(c,T)=P(λ≤AF·λ0|t≥T)=∫0AF·λ0p(λ|t≥T)dλ=∫0AF·λ0P(t<T|λ)π(λ)dλ∫0∞P(t<T|λ)π(λ)dλ=∫0AF·λ0(1-Σr=0ce-λT(λT)rr!)π(λ)dλ∫0∞(1-Σr=0ce-λT(λT)rr!)π(λ)dλ---(15)]]>

使用方后验风险β(c,T)的计算公式为：

β(r,T)=P(λ>AF·λ1|t≥T)=∫AF·λ1∞p(λ|t≥T)dλ=∫AF·λ1∞P(t≥T|λ)π(λ)dλ∫0∞P(t≥T|λ)π(λ)dλ=∫AF·λ1∞(1-Σr=0ce-λT(λT)rr!)π(λ)dλ∫0∞(1-Σr=0ce-λT(λT)rr!)π(λ)dλ---(16)]]>

为了获得验证试验方案，需要解下列满足双方风险的约束条件：

α(c,T)=P(λ≤AF·λ₀|t≥T)≥1-α           (17)

β(c,T)=P(λ>AF·λ₁|t≥T)≤β             (18)

步骤四、基于历史数据给出验证指标参数的先验分布，基于马尔科夫蒙特卡洛(Markov Chain Monte Carlo,MCMC)方法，利用WinBUGS14计算满足条件的方案集合；

贝叶斯试验设计一般很难获得后验分布的数学表达式，为了解决这个问题，本发明使用基于马尔科夫蒙特卡洛(Markov Chain Monte Carlo,MCMC)方法，并利用WinBUGS14计算满足条件的方案。

具体的，以指数分布型产品为例，由步骤三可知，指数分布型产品对于验证指标失效率λ，根据Bayes理论，一般取其共轭先验分布即Gamma分布，记为Gamma(a,b)，分别定义参数a和b均服从分层先验分布Gamma(η,κ).

为了描述加速因子的不确定性，将加速因子的波动以概率分布的形式进行描述，为了简单起见，定义加速因子服从的概率分布为均匀分布U(A,B).

则由MCMC方法并利用WinBUGS14可得到参数的后验均值，假设蒙特卡洛模拟的次数为N次，λ的N次后验抽样计为λ^(j)，先验分布两参数a，b的N次后验抽样分别计为(a^(j),b^(j))，加速因子AF的N次抽样计为AF^(j),令

g(j)(r)=(b(j))a(j)Trr!Γ(a(j))(T+b(j))r+a(j)---(19)]]>

则生产方的后验风险(15)可写为：

α(c,T)=Σj=1N[1-Σr=0ce-λ(j)T(λ(j)T)rr!]I(λ(j)≤λ0·AF(j))Σj=1N[1-Σr=0ce-λ(j)T(λ(j)T)rr!]=Σj=1N[Σr=0cg(j)(r)γ(r+a(j),(T+b(j))λ0AF(j))]Σj=1N[Σr=0cg(j)(r)Γ(r+a(j))]---(20)]]>

使用方的后验风险为：

β(r,T)=Σj=1N[1-Σr=0ce-λ(j)T(λ(j)T)rr!]I(λ(j)≥λ1·AF(j))Σj=1N[1-Σr=0ce-λ(j)(λ(j)T)rr!]=Σj=1N[Σr=0cg(j)(r)(Γ(r+a(j))-γ(r+a(j),(T+b(j))λ1AF(j)))]Σj=1N[Σr=0cg(j)(r)Γ(r+a(j))]---(21)]]>

式中的是不完全gamma函数。这时候求解关于双方风险的约束条件(17)和(18)，可以得到满足约束条件的方案集合(c,T)。

步骤五、确定试验的费用约束，计算得到最优方案。

由(17)和(18)可获得满足双方风险条件的试验方案(c,T)集合，以及实际的双方风险，在考虑试验费用的情况下，我们可以对方案进行进一步优化，考虑优化参数为c,T，α，β，假设下列参数的定义如下：

a₁：与试验时间相关的试验费用，包括试验过程中的电力，物力及人力损耗，试验时间越长损耗越高。

a₂₁：与试验样品有关的试验费用，失效的样品越多则损耗越高。

a₂₂：经过试验但未失效的样本造成的试验损失，这些样本虽未失效，但是已经过了一定时间的使用损耗，不具有最初的性能。

a₃：与生产方风险有关的试验费用，包括产品的重新设计与生产延误、销售市场的萎缩等。

a₄：与使用方风险有关的试验费用,括产品使用过程当中的维护保障、任务延迟损失等。

由以上参数可得考虑产品试验费用为：

f(c,T,α,β)=a₁·T+a₁₂·c+a₂₂·(n-c)+a₃·α+a₄·β     (22)

优化目标即使上述试验费用最小。

2.根据权利要求1所述的一种基于贝叶斯理论的定时截尾加速验收抽样试验优化设计方法，其特征在于，步骤四中，具体的抽样模拟过程采用下面的步骤：

子步骤1.在参数λ对应的先验分布Gamma分布中抽取i次参数λ_i，在抽样分布f(x|λ_i)中抽取j次仿真失效数据x_ij。

子步骤2.λ_i对应的Gamma分布的两参数分别为a_i和b_i，结合分层先验分布Gamma(η,κ)和仿真失效数据x_ij，利用Winbugs软件和MCMC方法可得到参数a_i和b_i的后验平均值E(π(a_i|x_ij))和E(π(b_i|x_ij))。

子步骤3.设定蒙特卡洛模拟次数，然后重复子步骤1和子步骤2仿真得到足够的数据。

利用仿真数据求解(20)和(21)方程组，得到满足条件的方案集合(c,T)。