CN106383946A - 一种加速退化模型参数估计方法 - Google Patents

Abstract

本发明涉及一种实用的加速退化模型参数估计方法，该方法主要包括以下步骤：（1）基于加速退化模型建立似然函数；（2）定义一个matlab函数，实现对数似然函数求和；（3）利用最小二乘法确定出待估参数的初值；（4）调用matlab软件中的fminsearch函数，求取待估参数的极大似然估计值；（5）判断所得极大似然估计值是否正确。该方法可避免传统极大似然法估计加速退化模型参数时适应性差的弱点，建立一种适用性强、易于使用的参数估计方法。

Description

一种加速退化模型参数估计方法

技术领域

本发明属于可靠性工程领域，涉及一种加速退化模型参数估计方法。

背景技术

利用加速退化试验技术评估产品可靠性指标或寿命信息的必要环节是对加速退化数据进行准确的统计分析，核心工作是根据加速退化数据估计出加速退化模型的各参数值。目前的参数估计方法主要有极大似然法、最小二乘法、矩估计法、最佳线性估计法等，由于加速退化数据具有大样本量特点，采用极大似然法的参数估计效果最好。

利用传统极大似然法估计加速退化模型参数的一般步骤为，首先建立包含模型各参数的对数似然函数，然后由对数似然函数获取各参数的偏导，最后解出各偏导为0时对应的参数值。然而，对于由加速退化模型建立的对数似然函数，未知参数的偏导表达式普遍繁琐，导致求解未知参数的工作量很大甚至无法解得未知参数值，因此传统极大似然法的适用性不强。为避免传统方法的不足，介绍一种利用Matlab软件中的fminsearch函数配合最小二乘法实现极大似然估计的参数求解方法。

发明内容

本发明的目的在于建立一种适用性强、易于使用的加速退化模型参数估计方法，该方法的具体技术方案为：

步骤一：基于加速退化模型建立似然函数。

设某产品开展了温度应力加速退化试验，T_k为第k个加速温度，t_ijk为T_k下第j个产品的第i次测量时间,y_ijk为对应的退化测量值，Λ(t)＝t^Λ为时间函数，代表时间增量,Δy_ijk＝y_ijk-y_(i-1)jk代表退化增量，其中i＝1,2,···,H_jk；j＝1,2,···,N_k；k＝1,2,···,M。如果产品的性能退化服从Wiener随机过程，称性能退化模型为Wiener退化模型，根据Wiener加速退化模型的独立增量特性

${\mathrm{\Δy}}_{ijk}~N\left(\exp \right({\mathrm{\γ}}_1-{\mathrm{\γ}}_2/T_k){\mathrm{\Δ\Λ}}_{ijk},\exp (2{\mathrm{\γ}}_3-{\mathrm{\γ}}_2/T_k\left){\mathrm{\Δ\Λ}}_{ijk}^2\right)---\left(1\right)$

式中，N(·)表示正态分布函数，γ₁,γ₂,γ₃,Λ为加速退化模型的4个待估参数。

根据式(1)，建立如下似然函数，

$L({\mathrm{\γ}}_1,{\mathrm{\γ}}_2,{\mathrm{\γ}}_3,\mathrm{\Λ})=\underset{k=1}{\overset{M}{\Π}}\underset{j=1}{\overset{N_k}{\Π}}\underset{i=1}{\overset{H_{jk}}{\Π}}\frac{1}{\sqrt{2\mathrm{\π}\exp (2{\mathrm{\γ}}_3-{\mathrm{\γ}}_2/T_k){\mathrm{\Δ\Λ}}_{ijk}}}\exp \[-\frac{{({\mathrm{\Δy}}_{ijk}-\exp ({\mathrm{\γ}}_1-{\mathrm{\γ}}_2/T_k\left){\mathrm{\Δ\Λ}}_{ijk}\right)}^2}{2\exp (2{\mathrm{\γ}}_3-{\mathrm{\γ}}_2/T_k){\mathrm{\Δ\Λ}}_{ijk}}\]---\left(2\right)$

步骤二：定义一个matlab函数，实现对数似然函数求和。

定义此matlab函数的步骤为，

步骤21：建立matlab函数如f＝MLEfun(x,data)，令x＝(x(1),x(2),x(3),x(4))分别对应待估参数(γ₁,γ₂,γ₃,Λ)，data为所有样品的加速退化数据(t_ijk,y_ijk,T_k)，f为MLEfun的输出；

步骤22：在函数MLEfun内实现对数似然函数求和

步骤23：令从而将求极大值问题转换为求极小值问题。

步骤三：利用最小二乘法确定出待估参数的初值。

相比于传统估计方法，基于fminsearch函数的极大似然估计方法不受复杂的似然函数所局限，适用性较强。然而此方法对参数初值选取的要求较高，初值设置不当无法获取最优解，因此通常需要多次试设参数初值，这在未知参数数量较多时会造成较大工作量，例如采用多应力加速模型时，可利用以下最小二乘法确定参数初值。实现步骤为，

步骤31：根据Wiener退化模型建立如下似然函数，

$L({\mathrm{\μ}}_k,{\mathrm{\σ}}_k,{\mathrm{\Λ}}_k)=\underset{j=1}{\overset{N_k}{\Π}}\underset{i=1}{\overset{H_{jk}}{\Π}}\frac{1}{\sqrt{2{\mathrm{\π\σ}}_k^2{\mathrm{\Δ\Λ}}_{ijk}}}\exp \{-\frac{{({\mathrm{\Δy}}_{ijk}-{\mathrm{\μ}}_k{\mathrm{\Δ\Λ}}_{ijk})}^2}{2{\mathrm{\σ}}_k^2{\mathrm{\Δ\Λ}}_{ijk}}\}---\left(3\right)$

式中，μ_k,σ_k,Λ_k分别表示Wiener退化模型在加速应力T_k下的漂移参数，扩散参数和时间参数；

步骤32：分别将各加速应力下的加速退化数据代入式(3)后极大化，得到参数估计值向量

步骤33：利用Arrhenius加速模型将表示为

${\widehat{\mathrm{\μ}}}_k=\exp ({\mathrm{\γ}}_1-{\mathrm{\γ}}_2/T_k)---\left(4\right)$

将表示为

${\widehat{\mathrm{\σ}}}_k=\exp ({\mathrm{\γ}}_2-0.5{\mathrm{\γ}}_2/T_k)---\left(5\right)$

步骤34：分别对式(4)，式(5)两边求对数，得到

$\ln {\widehat{\mathrm{\μ}}}_k={\mathrm{\γ}}_1-{\mathrm{\γ}}_2/T_k---\left(6\right)$

$\ln {\widehat{\mathrm{\σ}}}_k={\mathrm{\γ}}_3-0.5{\mathrm{\γ}}_2/T_k---\left(7\right)$

步骤35：建立如下矩阵方程

$\left[\begin{array}{c}\ln {\widehat{\mathrm{\μ}}}_1\\ \mathrm{...}\\ \ln {\widehat{\mathrm{\μ}}}_M\\ \ln {\widehat{\mathrm{\σ}}}_1\\ \mathrm{...}\\ \ln {\widehat{\mathrm{\σ}}}_M\end{array}\right]=\left[\begin{array}{ccc}1& -1/T_1& 0\\ ...& ...& ...\\ 1& -1/T_M& 0\\ 0& -0.5/T_1& 1\\ ...& ...& ...\\ 0& -0.5/T_M& 1\end{array}\right]\left[\begin{array}{c}{\mathrm{\γ}}_1\\ {\mathrm{\γ}}_2\\ {\mathrm{\γ}}_3\end{array}\right]---\left(8\right)$

解矩阵方程得到最小二乘估计令将作为参数初值。

步骤四：调用matlab软件中的fminsearch函数，求取待估参数的极大似然估计值。

调用形式为θ＝(θ(1),θ(2),θ(3),θ(4))分别为待估参数(γ₁,γ₂,γ₃,Λ)的极大似然估计值，fval为MLEfun(x,data)的极小值。

步骤五：判断所得极大似然估计值是否正确。

将每项参数初值分别增大5％或减小5％，按照步骤四中的方法求取极大似然估计值，如果每次所得结果一致，说明所得极大似然估计值正确；如果所得结果不一致，说明所得极大似然结果估计值有误，错误原因很可能为加速退化模型不合理，解决方法是重新建立加速退化模型后利用本方法进行参数估计。

附图说明

图1一种加速退化模型参数估计方法的步骤流程示意图。

具体实施方式

下面结合附图对本发明实现步骤进行进一步说明。

实施例：为了掌握某型宝石轴承支撑摆式加速度计的寿命信息，对此型加速度计开展了恒定应力加速退化试验。加速应力为温度T，分别为T₁＝338.16K，T₂＝348.16K，T₃＝358.16K，每个加速应力水平下的加速度计样本量为6。性能退化量y为一次项标度因素x相对于初始值x₀的变化y＝x-x₀，待样品冷却到室温25℃时测量y值，失效阈值为D＝0.006。具体的加速退化数据见表1.

表1加速度计加速退化数据

步骤一，建立如下似然函数，

$L({\mathrm{\γ}}_1,{\mathrm{\γ}}_2,{\mathrm{\γ}}_3,\mathrm{\Λ})=\underset{k=1}{\overset{3}{\Π}}\underset{j=1}{\overset{6}{\Π}}\underset{i=1}{\overset{H_{jk}}{\Π}}\frac{1}{\sqrt{2\mathrm{\π}\exp (2{\mathrm{\γ}}_3-{\mathrm{\γ}}_2/T_k){\mathrm{\Δ\Λ}}_{ijk}}}\exp \[-\frac{{({\mathrm{\Δy}}_{ijk}-\exp ({\mathrm{\γ}}_1-{\mathrm{\γ}}_2/T_k\left){\mathrm{\Δ\Λ}}_{ijk}\right)}^2}{2\exp (2{\mathrm{\γ}}_3-{\mathrm{\γ}}_2/T_k){\mathrm{\Δ\Λ}}_{ijk}}\]$

其中，H_j1＝11,H_j2＝11,H_j3＝10；

步骤二，定义一个matlab函数f＝MLEfun(x,data)，令x＝(x(1),x(2),x(3),x(4))分别对应待估参数(γ₁,γ₂,γ₃,Λ)，data为所有样品的加速退化数据(t_ijk,y_ijk,T_k)，为MLEfun的输出。

步骤三，利用最小二乘法确定出待估参数的初值。

求取Wiener退化模型在T_k(k＝1,2,···,M)下的参数估计值如表2所示

表2 Wiener退化模型在各加速应力下的参数估计值

采用最小二乘法解得

步骤四，调用matlab软件中的fminsearch函数，将作为初值代入，解得极大似然估计值

步骤五，分别将中每项增大5％、减小5％后执行步骤四，所得极大似然估计值都为

通过以上案例说明了本文方法的创造性和实用性，本案例中无法使用传统极大似然法解出极大似然估计值。

Claims (3)

1.一种加速退化模型参数估计方法，其特征在于，包括以下步骤：

步骤一，基于加速退化模型建立似然函数；

步骤二，定义一个matlab函数，实现对数似然函数求和；

步骤三，利用最小二乘法确定出待估参数的初值；

步骤四，调用matlab软件中的fminsearch函数，求取待估参数的极大似然估计值；

步骤五，判断所得极大似然函数是否正确，将每项参数初值分别增大5％或减小5％，按照步骤四中的方法求取极大似然估计值，如果每次所得结果一致，说明所得极大似然估计值正确；如果所得结果不一致，说明所得极大似然结果估计值有误，错误原因很可能为加速退化模型不合理，解决方法是重新建立加速退化模型后进行参数估计。

2.如权利要求1所述的加速退化模型参数估计方法，其特征在于，步骤二的具体方法为：

步骤21，建立matlab函数如f＝MLEfun(x,data)，令x＝(x(1),x(2),x(3),x(4))分别对应待估参数(γ₁,γ₂,γ₃,Λ)，data为所有样品的加速退化数据(t_ijk,y_ijk,T_k)，f为MLEfun的输出；

步骤22，在函数MLEfun内实现对数似然函数求和

步骤23，令从而将求极大值问题转换为求极小值问题。

3.如权利要求1所述的加速退化模型参数估计方法，其特征在于，步骤三的具体方法为：

步骤31，根据Wiener退化模型建立似然函数；

步骤32，分别将各加速应力下的加速退化数据代入似然函数后极大化，得到参数估计值向量；

步骤33，利用Arrhenius加速模型将表示为将表示为

步骤34，分别对式(4)，式(5)两边求对数，得到

步骤35，建立如下矩阵方程

$\left[\begin{array}{c}\ln {\widehat{\mathrm{\μ}}}_1\\ \mathrm{...}\\ \ln {\widehat{\mathrm{\μ}}}_M\\ \ln {\widehat{\mathrm{\σ}}}_1\\ \mathrm{...}\\ \ln {\widehat{\mathrm{\σ}}}_M\end{array}\right]=\left[\begin{array}{ccc}1& -1/T_1& 0\\ \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 1& -1/T_M& 0\\ 0& -0.5/T_1& 1\\ \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 0& -0.5/T_M& 1\end{array}\right]\left[\begin{array}{c}{\mathrm{\γ}}_1\\ {\mathrm{\γ}}_2\\ {\mathrm{\γ}}_3\end{array}\right]$

解矩阵方程得到最小二乘估计令将作为参数初值。