CN103218533B - 一种指数分布的串联系统可靠度置信下限估计方法 - Google Patents

Abstract

本发明给出了一种指数分布的串联系统可靠度置信下限估计方法，它包括八个步骤。它首先基于信仰推断给出系统内每个单元失效率λ_i的信仰分布，再根据指数分布的特点，从失效率角度出发得到了系统的组成结构关系，最后通过蒙特卡罗法得到在给定置信水平下的系统失效率上限，进而得到工程上特别关注的系统可靠度置信下限。本发明提出一种求解单元寿命服从指数分布的串联系统可靠度置信下限的精确方法，计算简便，对于随机截尾数据给出了数据填充算法，扩展了使用范围。

Description

一种指数分布的串联系统可靠度置信下限估计方法

技术领域

本发明涉及一种指数分布的串联系统可靠度置信下限估计方法，它可以精确估计单元寿命服从指数分布的串联系统可靠度置信下限，适用于可靠性评估、可靠性验证等相关领域。

背景技术

电子产品的寿命一般服从指数分布，关于寿命服从指数分布的电子产品某一单元而言，工程上常表现为随机截尾样本，而对这类不完全样本估计可靠度函数较为困难，因此本发明针对随机截尾样本的可靠度函数估计具有非常重要的现实意义。

电子产品通常包含多个单元，目前评估含多个单元组成的电子产品的系统可靠性的方法有LM法、MML法，二者均属于近似方法，计算方便，在工程上得到广泛应用，但其评估结果都偏保守，评估结果的精度难以控制。

系统可靠性评估的精确方法得到的评估结果具有较好的统计性质，但这类方法计算困难，工程实用性不强，因此，探索计算方便的系统可靠性评估精确方法具有十分重要的意义，而目前关于这方面的研究较为薄弱，为此本发明提供了一种指数分布串联系统可靠度置信下限估计方法。

发明内容

本发明给出了寿命服从指数分布的随机截尾样本的系统可靠度置信下限估计方法，研究较为成熟的定时截尾样本与定数截尾样本是随机截尾样本的特殊形式。首先通过分位数填充算法将随机截尾数据补充成虚拟完全样本，然后基于信仰推断给出系统内每个单元失效率λ_i的信仰分布，最后根据指数分布的特点，从失效率角度出发，通过蒙特卡罗法得到在给定置信水平下的系统失效率上限，进而得到工程上特别关注的系统可靠度置信下限。

（1）发明目的

在实际研制生产过程中，寿命试验的时间过长会导致费用的增加，有时还会使试验过期而失去作用，因此缩短试验时间而获得的随机截尾样本在工程实际中较为常见，由于随机截尾样本属于不完全样本，故相关可靠度指标难以得到，其中可靠度置信下限是产品设计定型的重要指标之一，但考虑到计算方便，目前在工程上一般采用近似方法得到产品的可靠度置信下限， 基于此本发明提供了一种单元服从指数分布的随机截尾样本的串联系统可靠度置信下限估计方法，它是一种计算简便的精确方法。

（2）技术方案：

本发明一种指数分布的串联系统可靠度置信下限估计方法，其具体实施所涉及到的理论知识与前提条件如下：

①伽马分布及相关性质

形状参数为r,尺度参数为λ的Gamma分布Ga(r,λ)是一种连续概率函数，其概率密度函数为：

$f\left(x\right)=\left\{\begin{array}{cc}\frac{x^{r-1}{e}^{-x/\mathrm{\λ}}}{{\mathrm{\λ}}^r\mathrm{\Γ}\left(r\right)},& x\≥0\\ 0,& x＜0\end{array}\right.$

Gamma分布相关性质：

1.r=1时的伽马分布就是指数分布，即

Ga(1,λ)=Exp(λ)

2.的伽马分布是自由度为n的χ²分布，即

$\mathrm{Ga}(\frac{n}{2},\frac{1}{2})={\mathrm{\χ}}_n^2$

3.伽马分布的可加性。设随机变量X～Ga(r₁,λ),Y～Ga(r₂,λ),，且X与Y独立，则

Z=X+Y～Ga(r₁+r₂,λ)。

4.伽马分布的线性性质。设X_i～Ga(α_i,a_iλ)(α_i>0,a_i>0),i=1,2,…,n,X₁,X₂,…,X_n相互独立，则

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{a}_i{X}_i~\mathrm{Ga}({\mathrm{\α}}_1+{\mathrm{\α}}_2+\·\·\·+{\mathrm{\α}}_n,\mathrm{\λ})$

5.伽马分布与卡方分布的关系。设X₁,X₂,…,X_n是指数分布Exp(λ)的随机样本，则

$2\mathrm{\λ}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{X}_i~{{\mathrm{\χ}}^2}_{2n}$

②卡方分布及相关性质

若k个随机变量Z₁、……、Z_k是相互独立，符合标准正态分布的随机变量（数学期望为0、方差为1），则随机变量X的平方和

$X=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{Z}_i^2$

被称为服从自由度为k的卡方分布，记作

X～X²(k).

卡方分布的概率密度函数为：

$f_k\left(x\right)=\frac{{(1/2)}^{k/2}}{\mathrm{\Γ}(k/2)}{x}^{k/2-1}{e}^{-x/2}$

其中x≥0，当x≤0时f_k(x)=0。这里Γ代表Gamma函数。

（3）实施方案：

假设某复杂系统由m个单元串联组成串联系统，如图1所示：其寿命分别服从失效率为λ₁,...λ_m的指数分布，则该串联系统失效率λ_*为

${\mathrm{\λ}}_{*}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_i---\left(1\right)$

式中m代表串联系统中的单元个数，λ₁,...λ_m表示第1，……m个单元的失效率，λ_*为该串联系统的失效率；

从指数分布F(t)=1-e^-λt,(t≥0)的总体中随机抽取n个样本进行寿命试验，t为时间，在试验中观察其工作状态，若某一时刻开始某一样本不能继续工作，则该时刻记为该样本的失效时间。试验结束后，有r个样本失效，记录其失效时间依次为t₁,...t_i...,t_r，剩下的n-r个样本

未失效，分别于τ_j时间有m_j个撤离试验，且∑m_j=n-r，其中，τ_j记为第j次撤离样本的时间，称为删失时间。

本发明一种指数分布的串联系统可靠度置信下限估计方法，该方法具体步骤如下：

步骤一：根据已观测到的失效时间t₁,...,t_r及删失时间τ_k利用极大似然估计得到

对于参数为λ指数分布F(t)=1-e^-λt,(t≥0)而言，按下面公式可以得到上述指数法分布的参数λ的极大似然估计

$\widehat{\mathrm{\λ}}=\frac{r}{\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{t}_i+\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\τ}}_j\·m_j}---\left(2\right)$

其中，步骤一所述的失效时间t₁,...,t_r是指从指数分布F(t)=1-e^-λt,(t≥0)的总体中随机抽取n个样本进行寿命试验，在试验中观察其工作状态，若某一时刻开始某一样本不能继续工作，则该时刻记为该样本的失效时间，试验结束后，有r个样本失效，记录其失效时间依次为t₁,...,t_r，剩下的n-r个样本未失效，分别于τ_j时间有m_j个撤离试验，且∑m_j=n-r，其中，τ_j记为第j次撤离样本的时间，称为删失时间，k为τ_j的个数，m_j是在τ_j时刻撤离试验的个数。是指数分布F(t)=1-e^-λt,(t≥0)的参数λ的极大似然估计。

步骤二：通过分位数填充算法将随机截尾数据按下式补充成虚拟完全样本 其中t₁,...,t_r即在步骤一中表示的失效时间，τ_ij表示在τ_i时刻撤离的第j个样本经补充后的虚拟失效时间，按下式得到：

${\mathrm{\τ}}_{\mathrm{ij}}=F^{-1}(F\left({\mathrm{\τ}}_i\right|\widehat{\mathrm{\λ}})+[1-F\left({\mathrm{\τ}}_i\right|\widehat{\mathrm{\λ}})\left]\frac{j}{m_i+1}\right),(i=1,\·\·\·,k;j=1,\·\·\·,m_i)---\left(3\right)$

式中符号说明如下：由得到F(t)的反函数j为在τ_i时刻 撤离的m_i个样本中的第j个，其他符号与步骤一一致，即k为τ_j的个数，m_j是在τ_j时刻撤离试验的个数，是指数分布F(t)=1-e^-λt,(t≥0)的参数λ的极大似然估计；

步骤三：基于虚拟完全样本通过信仰推断得到每个单元的失效率λ_i的分布F_i。其中λ_i为第i个单元的失效率，F_i为第i个单元的失效率的分布，具体算法如下：

由步骤二得到的虚拟完全样本来自于失效率为λ的指数分布的总体中，则样本的总试验时间为：

$T_n=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{t}_i+\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\underset{j=1}{\overset{m_i}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ij}}~\mathrm{Ga}(n,\mathrm{\λ})---\left(4\right)$

故 $2\mathrm{\λ}{T}_n~\mathrm{Ga}(\frac{2n}{2},\frac{1}{2})={\mathrm{\χ}}^2\left(2n\right),$

进而根据信仰推断得到：

$\mathrm{\λ}~\frac{{\mathrm{\χ}}^2\left(2n\right)}{2T_n}---\left(5\right)$

式中符号说明如下：寿命服从指数分布F(t)=1-e^-λt,(t≥0)的单元的失效率分布，该n个试验单元的总试验时间，χ²(2n)是自由度为2n的卡方分布;

步骤四：对于每一个单元而言，利用步骤三的式子（5）可以得到各单元失效率的分布，从该失效率分布中随机抽取一个样本λ_i，根据式子（1）计算该串联系统失效率

${\mathrm{\λ}}_{*}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_i---\left(1\right)$

式中符号说明如下：λ_*为该串联系统失效率，λ_i是从第i个单元的失效率分布函数中随机抽取的一个失效率样本，m为该串联系统的单元总个数；

步骤五：重复步骤四N次，每次计算得到的该串联系统失效率分别记为式中符号说明如下：N为计算系统失效率的总循环次数；j是第j次计算系统失效率，故 j=1,2,...N；为第j次计算得到的系统失效率的值；

步骤六：对从小到大排列，分别记为

式中符号说明如下：是N个系统失效率的值从小到大排列在各位次上的值，例如代表N个系统失效率的值从小到大排在第2位次上的值；

步骤七：给定1-α置信水平下，该串联系统的失效率上界为

${\mathrm{\λ}}_{*}^U={\mathrm{\λ}}_{*}^{\left(\right[N(1-\mathrm{\α})\left]\right)}---\left(6\right)$

式中符号说明如下：^[N(1-α)]代表向下取整，为该串联系统的失效率上界；

步骤八：得到串联系统在1-α置信水平下的可靠度置信下限为

$R_{*}^L\left(t\right)=e^{-{\mathrm{\λ}}_{*}^{U}t}---\left(7\right)$

式中符号说明如下：为该串联系统在1-α置信水平下的可靠度置信下限，t为时间，其它与步骤七一致，为该串联系统的失效率上界。

在步骤二里使用的分位数填充算法将删失样本补充成虚拟完全样本，再进行参数估计，其算法的收敛性及相合性，现有公知技术能够证明，使得系统可靠度函数估计的准确性得以提高。

在步骤三中，关于信仰推断结果的精确性在公知技术中都能证明在给定置信水平1-α下，由信仰推断得到的单元可靠度函数，具有频率意义下的实际置信水平1-α。这与经典方法得到的精确置信水平的置信界一致。因此基于信仰推断得到的失效率分布函数是精确的。

（3）本发明的优点：

i本发明针对工程实际中常见又较难处理的随机截尾寿命试验样本，提出一种解决单元寿命服从指数分布的串联系统可靠度置信下限估计的方法，是基于传统意义下定时截尾样本与定数截尾样本进行可靠性评估的一种泛化推广。

ii本发明提出的基于信仰推断与蒙特卡罗方法的系统可靠度置信下限估计是一种精确方法，根据指数分布的特点，从失效率角度出发，得到串联系统可靠度置信下限估计，简化了传统精确方法的计算困难。

附图说明

图1本发明串联系统结构图

图2随机截尾试验单元A/B/C及系统可靠度函数示意图

图3随机截尾试验系统可靠度置信下限估计流程图

图中符号说明如下：

m:串联系统的单元个数（见图1）

x₁,x₂,...x_m：串联系统中第1,2…m个单元（见图1）

i:单元编号；（见图3）

j:循环次数编号；（见3）

N：给定需循环次数；（见图3）

F_i：单元i的失效率的信仰分布；（见图3）

λ_i：单元i的一个失效率样本；（见图3）

：第j次循环得到的系统失效率；（见图3）

：目前观测到的随机截尾样本下单元失效率；（见图3）

具体实施方式

见图1-图3，本发明针对工程应用较常见的随机截尾样本，提出一种单元服从指数分布的串联系统可靠度置信下限估计方法，具体实施方式如下：

已知电子产品由3个电子单元串联组成，记为单元A、B、C，分别对与3个单元同批生产出来的单元进行随机截尾试验，单元A随机抽取10个，单元B抽12个，单元C抽15个，对于单元A的寿命试验，在发生第一部失效的时候，撤下3部未失效的单元A，在第二部和第三部失效时，分别撤下2部未失效的单元，试验过程依次记录的失效及删失数据见表1-3，设三种单元寿命都服从指数分布，试计算该系统在95%的置信水平下的可靠度。

表1单元A寿命试验数据

表2单元B寿命试验数据

表3单元C寿命试验数据

解：

步骤一：根据已观测到的失效时间t₁,...,t_r利用极大似然估计得到

${\widehat{\mathrm{\λ}}}_{a0}=\frac{3}{34+113+169+34\×3+113\×2+169\×2}=\frac{3}{982}$

${\widehat{\mathrm{\λ}}}_{b0}=\frac{5}{543\×3+843\×2+1634\×2+1734\×3+1809\×2}=\frac{5}{15403}$

${\widehat{\mathrm{\λ}}}_{c0}=\frac{7}{2164\×2+2292\×5+2419\×1+3176\×2+3338\×2+3666\×3}=\frac{7}{42233}$

步骤二：根据式(2)将随机截尾数据补充成虚拟完全样本

单元A的初始分布函数估计(t≥0)由题意可得：

τ₁=34,τ₂=113,τ₃=169,m₁=3,m₂=2,m₃=2,

待补充数据为：τ₁₁,τ₁₂,τ₁₃,τ₂₁,τ₂₂,τ₃₁,τ₃₂,

其中， ${\mathrm{\τ}}_{1j}=F^{-1}(F\left({\mathrm{\τ}}_1\right|{\widehat{\mathrm{\λ}}}_{a0})+[1-F\left({\mathrm{\τ}}_1\right|{\widehat{\mathrm{\λ}}}_{a0})\left]\frac{j}{m_1+1}\right),$ (j=1,...,m₁)其他类似得到：

${\mathrm{\τ}}_{11}=-\frac{982}{3}\ln \{1-[(1-e^{-\frac{3}{982}\×34})+\frac{1}{3+1}(1-1+e^{-\frac{3}{982}\×34})\left]\right\}\≈128.2$

${\mathrm{\τ}}_{12}=-\frac{982}{3}\ln \{1-[(1-e^{-\frac{3}{982}\×34})+\frac{2}{3+1}(1-1+e^{-\frac{3}{982}\×34})\left]\right\}\≈260.9$

${\mathrm{\τ}}_{13}=-\frac{982}{3}\ln \{1-[(1-e^{-\frac{3}{982}\×34})+\frac{3}{3+1}(1-1+e^{-\frac{3}{982}\×34})\left]\right\}\≈487.8$

${\mathrm{\τ}}_{21}=-\frac{982}{3}\ln \{1-[(1-e^{-\frac{3}{982}\×113})+\frac{1}{2+1}(1-1+e^{-\frac{3}{982}\×113})\left]\right\}\≈245.7$

${\mathrm{\τ}}_{22}=-\frac{982}{3}\ln \{1-[(1-e^{-\frac{3}{982}\×113})+\frac{2}{2+1}(1-1+e^{-\frac{3}{982}\×113})\left]\right\}\≈472.6$

${\mathrm{\τ}}_{31}=-\frac{982}{3}\ln \{1-[(1-e^{-\frac{3}{982}\×169})+\frac{1}{2+1}(1-1+e^{-\frac{3}{982}\×169})\left]\right\}\≈301.7$

${\mathrm{\τ}}_{32}=-\frac{982}{3}\ln \{1-[(1-e^{-\frac{3}{982}\×169})+\frac{2}{2+1}(1-1+e^{-\frac{3}{982}\×169})\left]\right\}\≈528.6$

同理可以得到单元B/C的虚拟完全样本

单元B：

τ₁₁≈1972.1,τ₁₂≈3927.4,

τ₂₁≈2978.3,

τ₃₁≈4069.3,

τ₄₁≈2983.1,τ₄₂≈5118.4,

τ₅₁≈3944.3

单元C：

τ₁₁≈6346,

τ₂₁≈4027.7,τ₂₂≈6474,τ₂₃≈10655.9

τ₄₁≈7358,

τ₅₁≈7520,

τ₆₁≈6112.3,τ₆₂≈10294.2

步骤三：基于以上虚拟完全样本通过信仰推断得到每个单元的失效率λ_i的分布F_i

$T_{\mathrm{na}}=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{t}_i+\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\underset{j=1}{\overset{m_i}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ij}}$ $=316+128.2+260.9+487.8+245.7+472.6+301.7+528.6=2741.5$

${\mathrm{\λ}}_a~\frac{{\mathrm{\χ}}^2\left(2n\right)}{2T_{\mathrm{na}}}=\frac{{\mathrm{\χ}}^2\left(20\right)}{5483}$

$T_{\mathrm{nb}}=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{t}_i+\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\underset{j=1}{\overset{m_i}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ij}}$ $=6563+1792.1+3927.4+2978.3+4069.3+2983.1+5118.4+3944.3=31375.9$

${\mathrm{\λ}}_b~\frac{{\mathrm{\χ}}^2\left(2n\right)}{2T_{\mathrm{nb}}}=\frac{{\mathrm{\χ}}^2\left(24\right)}{62751.8}$

$T_{\mathrm{nc}}=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{t}_i+\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\underset{j=1}{\overset{m_i}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ij}}$ $=19347+6346+4027.7+6474+10655.9+7358+7520+6112.3+10294.2=78135.1$

${\mathrm{\λ}}_c~\frac{{\mathrm{\χ}}^2\left(2n\right)}{2T_{\mathrm{nc}}}=\frac{{\mathrm{\χ}}^2\left(30\right)}{156270.2}$

步骤四：分别从A/B/C三个分布λ_a,λ_b,λ_c中随机抽取一个样本计算系统失效率

${\mathrm{\λ}}_{*}^j={\mathrm{\λ}}_a^j+{\mathrm{\λ}}_b^j+{\mathrm{\λ}}_c^j$

步骤五：重复步骤四N=10000次，每次计算得到的系统失效率分别记为步骤六：对从小到大排列，分别记为

步骤七：根据式（6），取1-α=95%置信水平下，系统的失效率上界为：

${\mathrm{\λ}}_{*}^U=0.0063$

步骤八：故根据式（7）得到该电子产品在1-α=95%置信水平下的可靠度函数置信下限为：

$R_{*}^L\left(t\right)=e^{-{\mathrm{\λ}}_{*}^{U}t}=e^{-0.0063t}$ 。

Claims (1)

1.一种指数分布的串联系统可靠度置信下限估计方法，其特征在于：该方法具体步骤如下：

步骤一：根据第i个单元已观测到的失效时间t₁,t₂,…,t_r及删失时间τ₁,τ₂,…,τ_k利用极大似然估计得到

对于失效率为λ_i指数分布而言，按下面公式得到上述指数分布的失效率λ_i的极大似然估计

其中，所述的失效时间t₁,…,t_r是指从指数分布的总体中随机抽取n个样本进行寿命试验，在试验中观察其工作状态，若某一时刻开始某一样本不能继续工作，则该时刻记为该样本的失效时间，试验结束后，有r个样本失效，记录其失效时间依次为t₁,…,t_r，剩下的n-r个样本未失效，分别于τ_j时间有m_j个未失效样本撤离试验，且其中，τ_j记为第j次撤离样本的时间，称为第j次删失时间；k为不同删失时间τ_j的个数，m_j是在τ_j时刻撤离试验的个数，是指数分布的失效率λ_i的极大似然估计；

步骤二：通过分位数填充算法将随机截尾数据按下式补充成虚拟完全样本

其中t₁,…,t_r即在步骤一中表示的失效时间，τ_jv表示在τ_j时刻撤离的第v个样本经补充后的虚拟失效时间，按下式得到：

式中符号说明如下：由得到F(t)的反函数τ_jv表示在τ_j时刻撤离的m_j个样本中的第v个样本的虚拟失效时间，其他符号与步骤一含义相同，即k为不同删失时间τ_j的个数，m_j是在τ_j时刻撤离试验的样本个数，是指数分布的失效率λ_i的极大似然估计；

步骤三：基于虚拟完全样本通过信仰推断得到第i个单元的失效率λ_i和失效分布F_i；其中λ_i为第i个单元的失效率，F_i为第i个单元的失效分布，具体算法如下：

由步骤二得到的虚拟完全样本来自于失效率为λ_i的指数分布的总体中，则第i个单元的样本总试验时间为：

故

进而根据信仰推断得到：

式中符号说明如下：寿命服从指数分布的第i个单元的失效率分布，T_n表示第i个单元的全部n个试验样本的总试验时间，χ²(2n)是自由度为2n的卡方分布；

步骤四：对于每一个单元而言，利用步骤三的式子(5)得到各单元失效率的分布，从该失效率分布中随机抽取一个样本λ_i，根据式子(1)计算该串联系统失效率

式中符号说明如下：λ_*为该串联系统失效率，λ_i是从第i个单元的失效率分布函数中随机抽取的一个失效率样本，m为该串联系统的单元总个数；

步骤五：重复步骤四N次，每次计算得到的该串联系统失效率分别记为

式中符号说明如下：N为计算系统失效率的总循环次数；B是第B次计算系统失效率，故B＝1,2,…,N；为第B次计算得到的系统失效率的值；

步骤六：对从小到大排列，分别记为

式中符号说明如下：是N个系统失效率的值从小到大排列在各位次上的值；

步骤七：给定1-α置信水平下，该串联系统的失效率上界为

式中符号说明如下：[N(1-α)]代表向下取整，为该串联系统的失效率上界；

步骤八：得到串联系统在1-α置信水平下的可靠度置信下限为

式中符号说明如下：为该串联系统在1-α置信水平下的可靠度置信下限，t为时间，其它与步骤七一致，为该串联系统的失效率上界。