CN102982208A - 一种基于Bayes因子优化的动态可靠性模型更新方法 - Google Patents

Abstract

本发明涉及一种基于Bayes因子优化的动态可靠性模型更新方法，包括如下步骤：步骤1：分析产品的失效物理，建立产品的失效物理可靠性模型，估计在某些给定时刻点处的可靠度R₀(t_i)；步骤2：现场收集产品相关的可靠性信息，对可靠性信息进行统计和物理特性分析，建立统计可靠性模型，估计在某些给定时刻点处的可靠度；步骤3：将失效物理可靠性模型作为选择模型，结合现场可靠性数据和可靠性试验数据，建立不同时间区间的Bayes因子表达式；本发明的有益效果是：实现小样本下的动态可靠性模型更新，同时可为新技术、新材料下新一代产品的可靠性设计优化，具有较高的实用价值和较强的工程意义。

Description

一种基于Bayes因子优化的动态可靠性模型更新方法

技术领域

本发明属于可靠性分析技术领域，尤其涉及机电产品的可靠性分析技术领域。

背景技术

目前产品可靠性分析多是基于统计数据或者通过经验确定的失效模型来实现的，因此也称为统计可靠性分析方法。这是一种自上而下的可靠性分析方法，需要大量的现场或试验数据做支撑。同时该可靠性分析方法无法建立统计数据与产品失效物理之间的联系，因此无法满足新材料和新技术应用下的可靠性评估要求。

基于失效物理的可靠性分析方法主要是从影响材料和系统退化和失效的机理出发，充分考虑材料性能和状态的变化规律、来自外部和内部各种能量的作用以及影响这些过程的不确定性因素，实现基于失效物理的可靠性分析。这是一种自下而上的可靠性分析方法，它能够从根本上处理产品的退化、失效过程，且同时将对产品性能产生影响的因素考虑进去，从而可以实现新材料、新技术下的产品可靠性分析与设计。

随着科学技术的不断发展，产品的高可靠、长寿命的特点表现的越来越明显，导致可用的现场数据越来越少，这给数据统计可靠性分析工作带来了很大困难。因此，为了克服对现场数据的依赖性，基于失效物理的可靠性分析方法逐渐成为研究的热点。但在对失效机理的推理分析及计算机仿真的过程中经常采用简化模型或近似加载来替代产品的真实状况。为了能够减少基于失效物理的可靠性模型与真实状况下产品可靠性模型之间的差别，许多研究者开展了可靠性模型更新方面的研究工作。目前普遍使用的可靠性模型更新方法多是基于贝叶斯方法或极大似然估计方法。为了能够表征更新后的可靠性模型与真实状况下的可靠性模型之间的吻合程度，业已提出了各种验证指标，包括：传统假设检验、贝叶斯因子指标、Frequentist指标和Area指标。传统假设检验和Frequentist指标是采用置信区间进行验证的；贝叶斯因子指标是通过建立在贝叶斯假设检验基础上的贝叶斯因子的取值大小进行验证的；Area指标是通过预测分布函数与观察分布函数之间的干涉面积进行验证的。目前的研究工作主要关注在可靠性模型更新方法和验证指标的表达上，而尚未开展更新信息需求与可靠性预测精度之间的关系问题的研究。

发明内容

本发明的目的是针对目前可靠性模型更新技术存在的缺陷，提出了一种基于Bayes因子优化的动态可靠性模型更新方法，本方法通过Bayes因子优化方法建立信息更新与可靠性预测精度的关系，从而获取满足精度要求的失效物理可靠性模型参数，为产品可靠性设计、监测及维修策略的制定奠定理论基础和技术支撑。

本发明的技术方案是：一种基于Bayes因子优化的动态可靠性模型更新方法，包括如下步骤：

步骤1：分析产品的失效物理，建立产品的失效物理可靠性模型（全称基于失效物理的动态可靠性近似模型），利用蒙特卡洛仿真对失效物理可靠性模型进行动态可靠性分析，估计在某些给定时刻t_i(i＝1,2,...,l)点处的可靠度R₀(t_i)；

步骤2：现场收集产品相关的可靠性信息，对可靠性信息进行统计和物理特性分析，建立统计可靠性模型(也称动态可靠性模型），估计在某些给定时刻t_i(i＝1,2,...,l)点处的可靠度R₁(t_i)；

步骤3：将失效物理可靠性模型作为选择模型，结合现场可靠性数据和可靠性试验数据，建立不同时间区间的Bayes因子表达式，并作为Bayes因子优化模型中的约束条件；

步骤4:以失效物理可靠性模型中的模型参数为设计变量，模型参数主要指表征系统性能的相关参数，失效物理可靠度和统计可靠度之差平方和最小为目标，各时间区间内Bayes因子为约束条件，建立基于Bayes因子的可靠性模型参数更新的优化模型；

步骤5：利用Matlab中的fmincon函数步骤4中的优化模型进行求解，若得到最优解，则此最优解为满足要求的失效物理模型更新参数；若无解，则继续制定可靠性试验方案，通过可靠性试验增加现场可靠性数据作为原有现场可靠性数据的补充，再重复步骤2到步骤4，直至Bayes因子优化模型找到最优解。

本发明的有益效果是：本发明充分考虑同一产品的失效物理可靠性模型与统计可靠性模型的统一性以及现场数据和可靠性试验数据对产品寿命反映的真实性，建立现场数据（包括可靠性试验数据）与失效物理可靠性模型更新精度之间的关系，实现小样本下的动态可靠性模型更新，同时可为新技术、新材料下新一代产品的可靠性设计优化、维修策略制定提供技术支撑，从而缩短研制周期、减少研制和使用过程中的维修费用，具有较高的实用价值和较强的工程意义。

附图说明

图1是本发明动态可靠性模型参数更新方法的流程图。

图2是具体实施例中机床主轴的结构示意图。

图3是模型参数更新前失效物理可靠性模型的可靠度函数。

图4是模型参数更新后失效物理可靠性模型的可靠度函数。

具体实施方式

下面结合附图对本发明的通用方法进行说明：

如图1所示，一种基于Bayes因子优化的动态可靠性模型更新方法，包括如下步骤：

步骤1：分析产品的失效物理，建立产品的失效物理可靠性模型（全称基于失效物理的动态可靠性近似模型），利用蒙特卡洛仿真对失效物理可靠性模型进行动态可靠性分析，估计在某些给定时刻t_i(i＝1,2,...,l)点处的可靠度R₀(t_i)，l为自然数；

可靠度R₀(t_i)的计算公式如公式（1）所示：

$R_0\left(t_i\right)=\mathrm{Pr}\{g(d\left(t_i\right),X\left(t_i\right),P\left(t_i\right);t_i)>0\}$ 公式（1）

$=\underset{g(d\left(t_i\right),X\left(t_i\right),P\left(t_i\right);t_i)>0}{\∫}{f}_{X,P}(X\left(t_i\right),P\left(t_i\right))\mathrm{dxdp}$

公式（1）中：d(t_i)表示与时间相关的确定性设计变量；X(t_i)表示与时间相关的随机设计变量；P(t_i)表示与时间相关的随机设计参数；Pr{□}表示概率；f_X,P(X(t_i),P(t_i))表示在某些给定时刻t_i点时X(t_i)与P(t_i)的联合概率密度函数；g(d(t_i),X(t_i),P(t_i);t_i)>0表示设计安全区。

本步骤中，建立产品的失效物理可靠性模型是通过产品动态性能分析软件实现的，此过程被认为是本领域的惯用手段而不再详细描述。

步骤2：现场收集产品相关的可靠性信息，对可靠性信息进行统计和物理特性分析，建立统计可靠性模型(也称动态可靠性模型），估计在某些给定时刻t_i(i＝1,2,...,l)点处的可靠度R₁(t_i)；

利用中位秩公式计算给定时刻t_i(i＝1,2,...,l)点处的可靠度R₁(t_i)：

$R_1\left(t_i\right)=1-F_n\left(t_i\right)=1-\frac{i-0.3}{n+0.4}$ 公式（2）

公式(2)中,F_n(t_i)表示t_i时刻的失效概率；i表示t_i时刻失效样本寿命数量；n表示从现场和可靠性试验数据中获取的样本寿命总量；公式（2）主要是用来评估样本量比较小情况下的可靠度。

步骤3：将失效物理可靠性模型作为选择模型，结合现场可靠性数据和可靠性试验数据，建立不同时间区间的Bayes因子B_i表达式如公式（3）所示，并作为Bayes因子优化模型中的约束条件；

$B_i=(n+1)P\left(k_i\right|x_{0i}^{k_i},n)$

公式（3）

$=x_{0i}^{k_i}{(1-x_{0i}^{k_i})}^{n-k_i}(n+1)\left(\begin{array}{c}n\\ k_i\end{array}\right)$

公式(3)中：表示基于失效物理可靠性模型得到的产品寿命落在第i个时间区间[t_i-1，t_i]内的概率；n表示从现场和可靠性试验数据中获取的产品寿命样本总量，其中k_i表示从现场和可靠性试验数据中获取的产品寿命落在第i个时间区间[t_i-1，t_i]内的样本量。

步骤4:以失效物理可靠性模型中的模型参数为设计变量，模型参数主要指表征系统性能的相关参数，失效物理可靠度和统计可靠度之差平方和最小为目标，各时间区间内Bayes因子为约束条件，建立基于Bayes因子的可靠性模型参数更新的优化模型：

$\min \underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(R_0(P_{\mathrm{FM}};t_i)-R_1\left(t_i\right))}^2$

$s.t.B_i\≥B_i^{T}i=1,\·\·\·,m;$ 公式（4）

$P_{\mathrm{FM}}^U\≥P\≥P_{\mathrm{FM}}^L$

公式(4)中，m表示选取的时间区间个数；P_FM表示需要更新的失效物理可靠性模型的关键参数；B_i如公式（3）所示；

和

分别表示关键参数P_FM可能取值的上限和下限；

表示Bayes因子需要满足的目标值，根据实际情况进行选取，

越大说明更新后的失效物理模型与现实情况越吻合。

步骤5：利用Matlab中的fmincon函数对公式（4）中的优化模型进行求解，若得到最优解，则此最优解为满足要求的失效物理模型更新参数；若无解，则继续制定可靠性试验方案，通过可靠性试验增加现场可靠性数据作为原有现场可靠性数据的补充，再重复步骤2到步骤4，直至Bayes因子优化模型找到最优解。

本步骤中Matlab为一款商业数学软件，fmincon函数为非线性优化函数，两者均为本零用的惯用软件和函数，故不再详细描述，这并不会影响本步骤的实施。

为了进一步的对本发明做进一步详细说明，下面再结合一具体实施例对本发明的方案做一个说明。本实施例以考虑性能退化的机床主轴失效物理可靠性模型更新为实施例，在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

如图2所示，为本发明的方法所针对的机床主轴的结构示意图。图中，R1、R2和R3分别表示阶梯轴的外径；R4和R5分别表示阶梯轴的内径；L1、L2和L3分别表示阶梯轴的长度。其中，R2=R1+0.0025，R3=R1+0.0050，L2=L3。

按照图1所示的流程进行实施。一种基于Bayes因子优化的动态可靠性模型更新方法，针对如图2所示的机床主轴，其具体步骤为：

（1）.机床主轴过度的弯曲变形将直接影响所加工零件的精度，因此将主轴的变形作为考虑的主要失效模式。利用试验设计、有限元分析以及响应面模型近似方法，构建关于机床主轴最大挠度与设计变量之间的函数。其中采用的响应面的基本形式为：

$y_f={\mathrm{\&beta;}}_0+\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_i{x}_i+\underset{i<j}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_{\mathrm{ij}}{x}_i{x}_j$

公式（5）

公式（5）中，y_f表示挠度；x_i、x_j表示尺寸参数和载荷；k表示尺寸参数和载荷个数；β₀、β_i、β_ij(i＝1,…,k;j=1,…,k)表示待定模型参数。

经过数据拟合，各待定模型参数分别为β₀=0.0050、β₁=0.0451、β₂=0.4546、β₃=-0.1519、β₄=-0.0764、β₅=0.0535、β₆=-2.7101×10^-8、β₁₂=-4.7026、β₁₃=-3.2264、β₁₄=1.7895、β₁₅=0.0906、β₁₆=0.5844×10^-6、β₂₃=0.2318、β₂₄=-0.2923、β₂₅=-2.3735、β₂₆=-0.8465×10^-5、β₃₄=2.5071、β₃₅=1.7403、β₃₆=1.1281×10^-6、β₄₅=-0.7215、β₄₆=2.8713×10^-6、β₅₆=1.3697×10^-6。

将估计得到的待定模型参数带入公式（5），则数据拟合后的挠度表达式为：

y_f=0.0050+0.0451x₁+0.4546x₂-0.1519x₃-0.0764x₄+0.0535x₅-2.7101×10^-8x₆-4.7026x₁x₂-3.2264x₁x₃+1.7895x₁x₄+0.0906x₁x₅+0.5844×10^-6x₁x₆+0.2318x₂x₃-0.2923x₂x₄-2.3735x₂x₅-0.8465×10^-5x₂x₆+2.5071x₃x₄+1.7403x₃x₅+1.1281×10^-6x₃x₆-0.7215x₄x₅+2.8713×10^-6x₄x₆+1.3697×10^-6x₅x₆

其中：[x₁,x₂,x₃,x₄,x₅,x₆]＝[R₁,R₂,R₃,L₁,L₂,F_t]，R₁,R₂,R₃,L₁,L₂为尺寸参数，其含义参照附图2；F_t是主轴所受的切向力。R₁,R₂,R₃,L₁,L₂,F_t皆为随机变量，相关分布信息如表1所示。

表1随机变量的分布信息

变量	均值	方差	分布类型
				x1	0.0345(m)	5×10-5(m)	正态分布
x2	0.0280(m)	5×10-5(m)	正态分布

x3	0.0216(m)	5×10-5(m)	正态分布
				x4	0.0540(m)	5×10-5(m)	正态分布
x5	0.0701(m)	5×10-5(m)	正态分布
				x6	2241(N)	224.1(N)	正态分布

主轴使用过程中的腐蚀、磨损等将导致主轴尺寸的改变，从而使得主轴挠度随着使用时间不断增加。主轴挠度的增加量近似服从Gamma分布，

其中，D表示挠度随机增加量，Gama(□)表示Gamma分布，t表示时间，α₁,α₂,α₃表示Gamma分布的相关参数。根据已有的信息近似得到D~Gama(0.0015t^2.2,2.2)/200小时。考虑腐蚀、磨损等性能退化因素时，与主轴挠度相关的性能函数为：

g_f＝[n_f]-y_f-D；

其中，[n_f]为机床主轴加工过程中允许的最大挠度值，在本实施例中[n_f]=10^-4m。利用蒙特卡洛仿真计算各个时间点的可靠度，计算结果如表2所示。则模型参数更新前失效物理可靠性模型的可靠度函数如图3所示。

表2模型参数更新前基于失效物理可靠性近似模型的不同时刻点可靠度

	t＝0	t＝600	t＝1200	t＝1800	t＝2400	t＝3000
							R0(ti)	1.000	0.872	0.532	0.210	0.0572	0.0073

（2）.从现场收集到的满足主轴弯曲失效模式的寿命数据为：458,782,1021,1396,1651,2243;。因此，寿命样本总量n=6，落在时间区间[0,600]，[600,1200]，[1200,1800]，[1800,2400]，[2400,3000]内的寿命样本量分别为1，2，2，1，0。利用公式（2）计算各个时刻点处的可靠度如表3所示。

表3基于统计可靠性近似模型的不同时刻点可靠度

	t＝0	t＝600	t＝1200	t＝1800	t＝2400	t＝3000
							R0(ti)	1.000	0.891	0.578	0.266	0.109	0.109

（3）.将性能退化参数α₁,α₂,α₃作为需要更新的关键模型参数，则将步骤（2）中的相关参数代入公式（3），则各个时间区间的Bayes因子为：

$B_1=\left(x_{01}({\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,{\mathrm{\&alpha;}}_3)\right){(1-\left(x_{01}({\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,{\mathrm{\&alpha;}}_3)\right))}^5\&times;7\&times;\begin{array}{}\left(\begin{array}{c}6\\ 1\end{array}\right)\end{array}$

$B_2={\left(x_{02}({\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,{\mathrm{\&alpha;}}_3)\right)}^2{(1-\left(x_{02}({\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,{\mathrm{\&alpha;}}_3)\right))}^4\&times;7\&times;\begin{array}{c}\left(\begin{array}{c}6\\ 2\end{array}\right)\end{array}$

$B_3={\left(x_{03}({\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,{\mathrm{\&alpha;}}_3)\right)}^2{(1-\left(x_{03}({\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,{\mathrm{\&alpha;}}_3)\right))}^4\&times;7\&times;\begin{array}{c}\left(\begin{array}{c}6\\ 2\end{array}\right)\end{array}$

$B_4=\left(x_{04}({\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,{\mathrm{\&alpha;}}_3)\right){(1-\left(x_{04}({\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,{\mathrm{\&alpha;}}_3)\right))}^5\&times;7\&times;\begin{array}{c}\left(\begin{array}{c}6\\ 1\end{array}\right)\end{array}$

$B_5={(1-\left(x_{01}({\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,{\mathrm{\&alpha;}}_3)\right))}^6\&times;7\&times;\begin{array}{c}\left(\begin{array}{c}6\\ 0\end{array}\right)\end{array};$

（4）.根据工程实际需求，基于Bayes因子可靠性模型参数更新的优化模型为

$\min \underset{i=1}{\overset{5}{\mathrm{\&Sigma;}}}(R_0({\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,{\mathrm{\&alpha;}}_3;t_i)-R_1\left(t_i\right))$

s.t.B_i≥3.6  i=1,…,5;

0.0001≤α₁≤0.01;

0.5≤α₂≤4;

0.1≤α₃≤5;

（5）.若步骤（4）中的优化模型无解，需要可靠性试验增加可靠性数据，在本实施例中共增加2次可靠性试验，得到的寿命数据分别为932,1452。重复步骤（2）~（4），得到步骤（4）中的优化模型的最优解为α₁=0.0165,α₂=2.19,α₃=2.0。则更新后失效物理可靠性模型的可靠度函数如图4所示。通过对比图3和图4，可以看出更新后失效物理可靠性模型与工程实际更加吻合。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (5)

1.一种基于Bayes因子优化的动态可靠性模型更新方法，包括如下步骤：

步骤1：分析产品的失效物理，建立产品的失效物理可靠性模型（全称基于失效物理的动态可靠性近似模型），利用蒙特卡洛仿真对失效物理可靠性模型进行动态可靠性分析，估计在某些给定时刻t_i(i＝1,2,...,l)点处的可靠度R₀(t_i)；

步骤2：现场收集产品相关的可靠性信息，对可靠性信息进行统计和物理特性分析，建立统计可靠性模型(也称动态可靠性模型），估计在某些给定时刻t_i(i＝1,2,...,l)点处的可靠度R₁(t_i)；

步骤3：将失效物理可靠性模型作为选择模型，结合现场可靠性数据和可靠性试验数据，建立不同时间区间的Bayes因子表达式，并作为Bayes因子优化模型中的约束条件；

步骤4:以失效物理可靠性模型中的模型参数为设计变量，模型参数主要指表征系统性能的相关参数，失效物理可靠度和统计可靠度之差平方和最小为目标，各时间区间内Bayes因子为约束条件，建立基于Bayes因子的可靠性模型参数更新的优化模型；

步骤5：利用Matlab中的fmincon函数步骤4中的优化模型进行求解，若得到最优解，则此最优解为满足要求的失效物理模型更新参数；若无解，则继续制定可靠性试验方案，通过可靠性试验增加现场可靠性数据作为原有现场可靠性数据的补充，再重复步骤2到步骤4，直至Bayes因子优化模型找到最优解。

2.根据权利要求1所述的一种基于Bayes因子优化的动态可靠性模型更新方法，其特征在于，步骤1中可靠度R0₍t_i)的计算公式如公式（1）所示： $R_0\left(t_i\right)=\mathrm{Pr}\{g(d\left(t_i\right),X\left(t_i\right),P\left(t_i\right);t_i)>0\}$ 公式（1）

$=\underset{g(d\left(t_i\right),X\left(t_i\right),P\left(t_i\right);t_i)>0}{\&Integral;}{f}_{X,P}(X\left(t_i\right),P\left(t_i\right))\mathrm{dxdp}$

公式（1）中：d(t_i)表示与时间相关的确定性设计变量；X(t_i)表示与时间相关的随机设计变量；P(t_i)表示与时间相关的随机设计参数；Pr{□}表示概率；f_X,P(X(t_i),P(t_i))表示在某些给定时刻t_i点时X(t_i)与P(t_i)的联合概率密度函数；g(d(t_i),X(t_i),P(t_i);t_i)>0表示设计安全区。

3.根据权利要求1所述的一种基于Bayes因子优化的动态可靠性模型更新方法，其特征在于，步骤2中对统计可靠性模型的计算是利用中位秩公式计算给定时刻t_i(i＝1,2,...,l)点处的可靠度R₁(t_i)： $R_1\left(t_i\right)=1-F_n\left(t_i\right)=1-\frac{i-0.3}{n+0.4}$ 公式（2）

公式(2)中,F_n(t_i)表示t_i时刻的失效概率；i表示t_i时刻失效样本寿命数量；n表示从现场和可靠性试验数据中获取的样本寿命总量；公式（2）主要是用来评估样本量比较小情况下的可靠度。

4.根据权利要求1所述的一种基于Bayes因子优化的动态可靠性模型更新方法，其特征在于，步骤3中建立不同时间区间的Bayes因子表达式的具体过程为：

$B_i=(n+1)P\left(k_i\right|x_{0i}^{k_i},n)$

公式（3）

$=x_{0i}^{k_i}{(1-x_{0i}^{k_i})}^{n-k_i}(n+1)\left(\begin{array}{c}n\\ k_i\end{array}\right)$

公式(3)中：

表示基于失效物理可靠性模型得到的产品寿命落在第i个时间区间[t_i-1，t_i]内的概率；n表示从现场和可靠性试验数据中获取的产品寿命样本总量，其中k_i表示从现场和可靠性试验数据中获取的产品寿命落在第i个时间区间[t_i-1，t_i]内的样本量。

5.根据权利要求4所述的一种基于Bayes因子优化的动态可靠性模型更新方法，其特征在于，步骤4中建立基于Bayes因子的可靠性模型参数更新的优化模型：

$\min \underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(R_0(P_{\mathrm{FM}};t_i)-R_1\left(t_i\right))}^2$ $s.t.B_i\&GreaterEqual;B_i^{T}i=1,\&CenterDot;\&CenterDot;\&CenterDot;,m;$ 公式（4）

$P_{\mathrm{FM}}^U\&GreaterEqual;P\&GreaterEqual;P_{\mathrm{FM}}^L$

公式(4)中，m表示选取的时间区间个数；P_FM表示需要更新的失效物理可靠性模型的关键参数；B_i如公式（3）所示；

和

分别表示关键参数P_FM可能取值的上限和下限；

表示Bayes因子需要满足的目标值，根据实际情况进行选取，

越大说明更新后的失效物理模型与现实情况越吻合。

Images