CN105426647B - 基于可靠度先验信息融合的冷备系统可靠度估计方法 - Google Patents

Abstract

本发明属于可靠性评估领域，具体涉及一种基于可靠度先验信息融合的冷备系统可靠度估计方法，包括以下步骤：(S1)获取部件的可靠度先验信息，将部件的可靠度先验信息转化为部件可靠度的验前分布；(S2)将部件可靠度的验前分布转化为分布参数的验前分布；(S3)根据分布参数的验前分布，求解分布参数的验后分布；(S4)基于分布参数的验后分布，对冷备系统的可靠度进行估计。本发明的目的在于提供一种方法，当部件寿命服从威布尔分布时，能够用于融合部件可靠度先验信息，并估计部件寿命的分布参数，进一步估计n个相同部件构成的n中取k冷备系统的可靠性，很好地解决了融合部件可靠度先验信息后，对部件构成的冷备系统的可靠度估计问题。

Description

基于可靠度先验信息融合的冷备系统可靠度估计方法

技术领域

本发明属于可靠度评估领域，具体涉及一种基于可靠度先验信息融合的冷备系统可靠度估计方法。

背景技术

可靠性是指产品在规定的条件下和规定的时间内完成规定功能的能力(具体参考文献：郭波，武小悦.系统可靠性分析[M].长沙:国防科技大学出版社,2002:5-6.)；它是产品的固有属性，是衡量产品质量好坏的重要指标。可靠性的概率度量称为可靠度，有时也常用产品的寿命这一指标进行衡量。相应地，产品在规定的条件下丧失规定的功能，则称之为故障。随着现代科学技术的发展，构成产品的元器件越来越多，产品的规模越来越庞大，研制和生产费用越来越高，这使得产品的可靠性问题变得越来越重要。工程上经常采用冗余技术来提高产品的可靠性，冷备是其中常见的一种方式。对于由n个相同部件组成的n中取k冷备系统，n，k为整数，任意时刻必须有k个部件工作，整个冷备系统才能正常工作，而剩余的n-k个部件则作为备份。当k个工作部件中有故障部件时，备份的部件立即顶替故障部件，直到有(n-k+1)个部件全部故障后，冷备系统才会故障。

对产品的可靠度进行准确的估计，有助于我们及时了解产品的运行情况，作出正确的决策。对可靠度的估计，通常是将产品的寿命视为随机变量，并认为产品的寿命服从某个特定分布，然后借助于数理统计理论进行分析。例如在理论分析和工程中，因为威布尔分布的良好特性，常用威布尔分布来拟合产品的寿命分布。威布尔分布的概率密度函数为

其中t为产品的寿命，m为威布尔分布的形状参数，η为威布尔分布的尺度参数，exp表示以自然对数e为底的指数函数。威布尔分布下的可靠度函数为:

由此可知，假如需要估计产品工作到τ时刻的可靠度，只要知道分布参数m和η的估计值和借助于式(2)即可求得可靠度R(τ)的估计值因此对可靠度的估计，关键在于对分布参数m和η的估计。

在实际中，往往需要首先利用一批试验样品进行可靠性寿命试验，收集试验样品的寿命数据，然后借助于统计分析理论，对分布参数和可靠度进行估计。假如收集得到的试验数据全部为故障数据，则称该组试验数据为完全样本，否则称之为截尾样本。Bayes理论在当前的可靠性分析中运用较多。Bayes理论将各种各样的其他可靠性信息视为先验信息，并转化为验前分布，进一步通过Bayes公式与收集到的试验数据相融合，然后估计分布参数和可靠度。因为在估计过程中，运用到了大量的可靠性信息，从而大大提高了估计的精度，所以Bayes理论受到了广泛运用。Bayes理论的核心是Bayes公式，如下式所示：

其中θ是需要运用Bayes理论估计的参数，π(θ)是参数θ的验前分布，D是试验数据构成的样本，L(D|θ)是根据样本求得的似然函数，∫_θπ(θ)L(D|θ)dθ是关于试验数据D的边缘分布，π(θ|D)是参数θ的验后分布。

冷备系统作为可靠性工程中常见的一种结构，对冷备系统的可靠度估计也有很多相应的方法，处理思路都是根据构成冷备系统的部件的寿命所服从的概率分布，结合冷备系统的结构特性和数理统计理论，对冷备系统的可靠度进行分析求解。目前针对冷备系统的研究大多都是假定部件寿命的概率分布的分布参数已知，从而单纯地研究冷备系统的可靠性。但在实际应用中，分布参数往往是未知的，需要首先进行估计，然后才能估计冷备系统的可靠度。而目前，将分布参数估计与冷备系统的可靠性结合起来的研究还相对较少。另外，在工程中，除了可靠性试验收集得到的部件寿命数据之外，还存在着一些关于部件可靠性的先验信息。这些先验信息可以与部件寿命数据一起用于估计分布参数。而目前极其缺乏在这种情况下，估计冷备系统的可靠度方法。本发明通过对现有技术的整合和改进，要解决的技术问题在于：(1)当部件的寿命服从威布尔分布时，如何结合Bayes理论对威布尔分布的分布参数(m,η)进行估计；当部件寿命服从威布尔分布时，按照常用作法，认为可靠度的验前分布服从负对数伽马分布。负对数伽马分布的概率密度函数为：

其中a，b是负对数伽马分布的分布参数，Γ(a)为Γ函数，Γ(a)＝∫₀ ^∞y^a-1e^-ydy。

(2)如何根据前一步得到的分布参数(m,η)，结合n中取k冷备系统的可靠度估计方法，对n中取k冷备系统的可靠度进行估计。

发明内容

为了解决上述技术问题，本发明主要基于蒙特卡罗马尔可夫(Monte CarloMarkov Chain，MCMC)算法对相关的分布函数进行抽样处理，具体技术方案为：

一种基于可靠度先验信息融合的冷备系统可靠度估计方法，包括如下步骤：

(S1)获取部件的可靠度先验信息，将部件的可靠度先验信息转化为部件可靠度的验前分布；

(S2)将部件可靠度的验前分布转化为分布参数的验前分布；

(S3)根据分布参数的验前分布，求解分布参数的验后分布；

(S4)基于分布参数的验后分布，对冷备系统的可靠度进行估计。

进一步地，所述步骤(S1)的具体过程为：

(S11)记获取部件的可靠度先验信息为部件在时刻处的可靠度真值R_i的估计值其中i＝1,2,…,M,M≥2，i为自然数，M为整数；将估计值视为验前分布π(R_i|a_i,b_i)的期望值，令

其中，Γ(a_i)为Γ函数,根据上式得到分布参数a_i和b_i的关系

(S12)根据最大熵原理，确定分布参数a_i和b_i的值，令熵H最大，记为max H：

(S13)根据分布参数a_i和b_i，求得对应R_i的验前分布π(R_i|a_i,b_i)。

进一步地，所述步骤(S2)的具体过程为：

(S21)根据时刻i＝1,2,…,M处的可靠度R_i的验前分布π(R_i)，对每个π(R_i)依次进行抽样得到抽样值序列i＝1,2,…,M；

(S22)从抽样值序列i＝1,2,…,M中随机选择抽样值和其中 分别为时刻 处的验前分布π(R_u)，π(R_v)的抽样值，u,v＝1,2,…,M，u≠v；若 和 满足下列关系

则按照下式计算得到：

其中m^p和η^p视为分布参数(m,η)的验前分布π(m,η)的抽样值；

(S23)除去抽样值和判断抽样值序列中剩余的抽样值个数是否大于2个，若是，从剩余的抽样值序列i＝1,2,…,M,i≠u,v中继续随机选择抽样值，重复步骤(S22)，继续求解分布参数(m,η)的验前分布π(m,η)的抽样值；否则，进入步骤(S24)。

(S24)重复步骤(S21)-(S23)，直到得到的抽样值(m^p,η^p)个数达到预先设定值l，记为j＝1,2,…,l。i,j,u,v,l为自然数，M为整数。

进一步地，所述步骤(S3)的具体过程为：

(S31)记针对部件进行可靠性寿命试验收集到的数据为t₁,t₂,…,t_N，设其中的故障数据构成的集合为F；则可根据数据t₁,t₂,…,t_N计算得该样本的似然函数为

其中f(t_i；m,η)和R(t_j；m,η)具体格式如下：

令j＝1，记验后分布抽样值初值和为任意正数；

(S32)从验前分布π(m,η)的抽样值序列j＝1,2,…,l中，依次选择按下式计算

其中和即为式(9)中分布参数取不同值时的似然函数；

(S34)重复步骤(S32)-(S33)，直到j＝l；由此得到分布参数(m,η)的验后分布抽样值序列其中j＝1,2,…,l。

进一步地，所述步骤(S4)的具体过程为：

(S41)记步骤(S3)中得到分布参数(m,η)的验后分布抽样值序列其中j＝1,2,…,l，令初始值j＝1。

(S42)根据分布参数抽样值和基于分布参数为和的威布尔分布，生成n个随机数序列T₁,T₂,…,T_n。

(S43)将T_n＝(T₁,T₂,…,T_n)及n、k作为函数输入参数，调用递归函数fun_T(T_n,n,k)，计算n中取k冷备系统的寿命T^c，其中T^c＝fun_T(T_n,n,k)的定义如下：

如果k＝n，令

否则，记tm为前k个随机数T_k＝(T₁,…,T_k),(k＜n)的最小值，并从前k个随机数T_k＝(T₁,…,T_k)中去掉tm，更新T_k中剩余的(k-1)个随机数为T_k-1＝(T₁-tm,…,T_k-tm)，即T_k中剩余的(k-1)个随机数分别减去tm；然后将k之后的(n-k)个随机数及更新后的T_k-1合并为T_n-1＝(T₁-tm,…,T_k-tm,T_k+1,…,T_n)；令T^c＝tm+fun_T(T_n-1,n-1,k)；

(S44)如果T^c≤τ，τ为冷备系统的任务时刻；令g＝1，否则令g＝0；

(S45)重复步骤(S42)-(S44)，直到循环次数达到预先设定的s，并得到数值序列g₁,g₂,…,g_s；

(S46)由此基于验后分布抽样值和通过仿真计算的方法得到n中取k冷备系统在τ处的可靠度估计值为：

(S47)令j＝j+1，返回至(S42)，直到j＞l后停止计算；

此时可获得l个冷备系统的可靠度估计值其中j＝1,2,…,l；根据蒙特卡罗马尔可夫算法，舍弃这l个可靠度估计值序列的初始部分，并对剩余估计值取平均，即可获得信息融合后冷备系统的可靠度的最终估计值为：

进一步地，至少要有2个针对部件的可靠度先验信息，即步骤(S11)中所述的M≥2。

采用本发明获得的有益效果：本发明提出的步骤简便易行，便于程序化处理。借助于计算机程序，可避免大量复杂的数学运算。本发明的目的在于提供一种方法，当部件寿命服从威布尔分布时，能够用于融合部件的可靠度先验信息，并估计部件寿命的分布参数，进一步估计n个相同部件构成的n中取k冷备系统的可靠性。本发明很好地解决了融合部件可靠度先验信息后，对部件构成的冷备系统的可靠度估计问题。

附图说明：

图1是本发明的流程图。

具体实施方式：

以下将结合具体实施例和附图对本发明做进一步详细说明。如图1所示，为本发明的流程图。

(S1)获取部件的可靠度先验信息，将部件的可靠度先验信息转化为部件可靠度的验前分布；

工程上存在的可靠度先验信息一般是部件在时刻处的可靠度真值R_i的估计值工程上获取的方式有很多，常用方式是专家判断或工程经验。因的获取方式不是本发明所要考虑的问题，因此在本发明中将视为已知值。若要根据Bayes公式来估计分布参数，则按照Bayes理论的原理，将可靠度真值R_i视为随机变量，同时需要将估计值对应转化为可靠度R_i的验前分布π(R_i)，i＝1,2,…,M。

当部件寿命服从威布尔分布时，按照常用作法，认为可靠度的验前分布服从负对数伽马分布。负对数伽马分布的概率密度函数为：

其中a，b是负对数伽马分布的分布参数，Γ(a)为Γ函数，Γ(a)＝∫₀ ^∞y^a-1e^-ydy。即在本发明中，认为π(R_i)为LΓ(R_i|a_i,b_i)。下面具体阐述如何根据已知的来确定验前分布π(R_i)中的分布参数a_i和b_i，其中i＝1,2,…,M。

将已知的估计值视为验前分布π(R_i)的期望值，即令

根据式(5)，可以得到分布参数a_i和b_i之间的关系进一步按照最大熵原理，确定分布参数a_i和b_i的值，即

本实施例的评估对象为3中取2的冷备系统，即n＝3，k＝2。针对部件，利用8个样品实施可靠性寿命试验，并收集到试验数据3,5,8,11,20,28,33,38，数据的单位是天。除了试验数据以外，还收集到专家对于部件可靠度的判断值，即部件在时刻为1天时的可靠度估计值为0.98，在时刻为10天时的可靠度估计值为0.79，它们被视为部件可靠度的先验信息。本实施例以估计3中取2冷备系统在任务时刻为5天时的可靠度估计值为例，阐述本发明的技术方案步骤。

取分布参数的验前分布形式为负对数伽马分布LΓ(x|a,b)。将部件在时刻为1天和10天时的可靠度估计值0.98和0.79分别转化为LΓ(x|a₁,b₁)和LΓ(x|a₂,b₂)，可得a₁＝0.9847,b₁＝48.2411及a₂＝0.8644,b₂＝3.1895。

(S2)将部件可靠度的验前分布转化为分布参数的验前分布；

为便于利用Bayes公式估计分布参数m和η，需要进一步将得到的M(M≥2)个验前分布π(R_i),(i＝1,2,…,M)转化为分布参数(m,η)的联合验前分布π(m,η)，然后才能得到分布参数(m,η)的联合验后分布π(m,η|D)。鉴于威布尔分布函数的复杂性，为了便于处理，结合MCMC方法的思想，通过抽样离散的方式，按照下列步骤进行求解转化。

本实施例中设定l＝1000，依次对LΓ(x|a₁,b₁)和LΓ(x|a₂,b₂)进行抽样，并依据式(7)和式(8)转化为分布参数(m,η)的验前分布π(m,η)的抽样值(m^p,η^p)，直到得到1000组抽样值(m^p,η^p)为止。

(S3)根据验前分布π(m,η)，求解分布参数(m,η)的验后分布

根据Bayes理论的原理，当获取分布参数的验前分布后，需要与可靠性寿命试验收集到的寿命数据进行融合，得到分布参数的验后分布，才能得到信息融合后分布参数的估计值。基于步骤(S2)已获得分布参数(m,η)的验前分布π(m,η)的抽样值序列j＝1,2,3,…,l，将直接根据抽样值序列j＝1,2,…,l获取分布参数(m,η)的验后分布抽样值，记为j＝1,2,3,…,l。

实施例中按照本发明方法，基于1000组抽样值(m^p,η^p)，结合试验数据3,5,8,11,20,28,33,38，按照步骤(S3)，得到1000组分布参数(m,η)的验后分布抽样值(m^f,η^f)。

(S4)基于得到的验后分布抽样值，对冷备系统的可靠度进行估计

针对n中取k冷备系统的可靠度估计，有很多方法。这里将继续根据上步中得到的验后分布抽样值j＝1,2,3,…,l，结合n中取k冷备系统的仿真估计方法，对冷备系统在τ处的可靠度进行估计。

实施例中，按照步骤(S4)获得1000个冷备系统的可靠度估计值j＝1,2,3,…,l，并最终依据公式(12)求得信息融合后，冷备系统在时刻为5天时的可靠度估计值实施例中l的取值一般都是1000以上的整数倍，本领域的默认做法若j＝0.1l+1出现小数，则从对该数取整数开始统计。

本发明首先根据部件的可靠度先验信息确定可靠度的验前分布，通过对部件可靠度的验前分布进行抽样，从而转化为分布参数的验前分布的抽样值，而后通过信息融合的方法，根据这些验前分布的抽样值得到分布参数的验后分布的抽样值，继而再与冷备系统的可靠度估计方法相结合，最终求得信息融合后冷备系统的可靠度估计值。通过上述步骤，本发明很好地解决了融合部件可靠度先验信息后，对部件构成的冷备系统的可靠度估计问题。

Claims (1)

1.一种基于可靠度先验信息融合的冷备系统可靠度估计方法，其特征在于，包括以下步骤：

(S1)获取部件的可靠度先验信息，将部件的可靠度先验信息转化为部件可靠度的验前分布；具体过程为：

(S11)记获取部件的可靠度先验信息为部件在时刻处的可靠度真值R_i的估计值将估计值视为验前分布π(R_i|a_i,b_i)的期望值，其中i＝1,2,…,M,M≥2，令

其中，Γ(a_i)为Γ函数，根据上式得到分布参数a_i和b_i的关系

(S12)根据最大熵原理，确定分布参数a_i和b_i的值，即令最大熵maxH：

其中将根据式(5)得到的关系式代入式(6)中，则确定分布参数a_i和b_i的问题就转化为单变量的优化问题，利用一维线性搜索方法求解分布参数a_i和b_i；

(S13)根据分布参数a_i和b_i，求得对应R_i的验前分布π(R_i|a_i,b_i)；

(S2)将部件可靠度的验前分布转化为分布参数的验前分布；具体过程为：

(S21)根据时刻处的可靠度R_i的验前分布π(R_i)，对每个π(R_i)依次进行抽样得到抽样值序列其中i＝1,2,…,M；

(S22)从抽样值序列中随机选择抽样值和其中分别为时刻处的验前分布π(R_u)，π(R_v)的抽样值，u,v＝1,2,…,M，u≠v；若和满足下列关系

则按照下式计算得到：

其中m^p和η^p视为分布参数(m,η)的验前分布π(m,η)的抽样值；

(S23)除去抽样值和判断抽样值序列中剩余的抽样值个数是否大于2个，若是，从剩余的抽样值序列中继续随机选择抽样值，其中i＝1,2,…,M,i≠u,v，重复步骤(S22)，继续求解分布参数(m,η)的验前分布π(m,η)的抽样值；否则，进入步骤(S24)；

(S24)重复步骤(S21)-(S23)，直到得到的抽样值(m^p,η^p)个数达到预先设定值l，记为

(S3)根据分布参数的验前分布，求解分布参数的验后分布；具体过程为：

(S31)记针对部件进行可靠性寿命试验收集到的数据为t₁,t₂,…,t_N，设其中的故障数据构成的集合为F；则可根据数据t₁,t₂,…,t_N计算样本的似然函数为

其中f(t_i；m,η)和R(t_j；m,η)具体格式如下：

其中m为威布尔分布的形状参数，η为威布尔分布的尺度参数；

令j＝1，记验后分布抽样值初值和为任意正数；

(S32)从验前分布π(m,η)的抽样值序列中依次选择其中j＝1,2,…,l；按下式计算

其中和即为式(3)中分布参数取不同值时的似然函数；

(S33)从均匀分布U(0,1)中生成随机数r，并令j＝j+1；如果则令否则令

(S34)重复步骤(S32)-(S33)，直到j＝l；由此得到分布参数(m,η)的验后分布抽样值序列其中j＝1,2,…,l；

(S4)基于分布参数的验后分布，对冷备系统的可靠度进行估计，具体过程为：

(S41)记步骤(S3)中得到分布参数(m,η)的验后分布抽样值序列其中j＝1,2,…,l，令初始值j＝1；

(S42)根据分布参数抽样值和基于分布参数为和的威布尔分布，生成n个随机数序列T₁,…,T_n；

(S43)将T_n＝(T₁,…,T_n)及n、k作为函数输入参数，调用递归函数fun_T(T_n,n,k)，计算n中取k冷备系统的寿命T^c，其中T^c＝fun_T(T_n,n,k)的定义如下：

如果k＝n，令

否则，记tm为前k个随机数T_k＝(T₁,…,T_k)的最小值，k＜n，并从前k个随机数T_k＝(T₁,…,T_k)中去掉tm，更新T_k中剩余的(k-1)个随机数为T_k-1＝(T₁-tm,…,T_k-tm)；然后将k之后的(n-k)个随机数及更新后的T_k-1合并为T_n-1＝(T₁-tm,…,T_k-tm,T_k+1,…,T_n)；令T^c＝tm+fun_T(T_n-1,n-1,k)；

(S44)如果T^c≤τ，τ为产品某个时刻点；令g＝1，否则令g＝0；

(S45)重复步骤(S42)-(S44)，直到循环次数达到预先设定的s，并得到数值序列g₁,g₂,…,g_s；

(S46)由此基于验后分布抽样值和通过仿真计算的方法得到n中取k冷备系统在τ处的可靠度估计值为：

(S47)令j＝j+1，返回至(S42)，直到j＞l后停止计算；

此时可获得l个冷备系统的可靠度估计值其中j＝1,2,…,l；根据蒙特卡罗马尔可夫算法，舍弃这l个可靠度估计值序列的初始部分，并对剩余估计值取平均，即可获得信息融合后冷备系统的可靠度的最终估计值为：