CN101793927A - 步进应力加速退化试验优化设计方法 - Google Patents

Abstract

本发明公开了一种步进应力加速退化试验优化设计方法，包括以下几个步骤，步骤一、确定产品性能退化模型，加速模型和可靠度函数；步骤二、建立目标函数；步骤三、确定约束条件和优化变量；步骤四、实现优化算法；本发明以p分位寿命可靠度的渐进方差最小为目标。计算p分位寿命可靠度的渐进方差时，采用差分方式对可靠度函数求偏导，该方法解决了复杂函数求导困难，容易出错等问题；根据工程经验提出对每一应力下监测次数的约束条件，使试验优化结果更具有工程实际意义，M₁≥M₂≥……≥M_K≥10的约束也大大减少了计算次数，提高了试验优化速率；将试验应力作为优化变量，并根据实际情况给出应力水平的约束条件，对试验方案更系统更全面地进行优化。

Description

步进应力加速退化试验优化设计方法

技术领域

本发明涉及一种步进应力加速退化试验的优化设计方法，属于加速退化试验技术领域。

背景技术

随着现代科技的发展，高可靠、长寿命产品越来越多，为了评估这些产品的可靠性与寿命，美国罗姆航空中心在20世纪60年代就正式提出加速寿命试验(ALT)，目前ALT技术的研究已经取得了大量研究成果，但长寿命和高可靠性产品在ALT中只有少量失效出现或根本没有失效出现，这就对产品的寿命与可靠性评估带来了极大的困难。为解决这一问题，加速退化试验(ADT)应运而生。ADT克服了ALT只记录产品失效时间而不考虑产品性能退化情况的不足，通过对加速退化数据处理可以对高可靠长寿命产品的可靠性及寿命进行较好的评估。

ADT需要在了解产品退化失效机理并建立加速模型之后实施，通过监控并记录产品的性能退化参数，了解其性能退化趋势，从而获得产品的退化数据，基于退化数据对产品进行可靠性与寿命的评估。因此，在不改变产品失效机理的前提下，如何对ADT方案进行科学合理的设计，以便快速有效的获得产品的性能退化数据，提高产品可靠性评估精度是当前加速退化试验研究的一个热点。

恒定应力加速退化试验(CSADT)是指选择一组高于正常应力的加速应力水平，将一定数量的试件分为N组，每组在一个应力下实施加速退化试验；而步进应力加速退化试验(SSADT)是指选择一组高于正常应力的加速应力水平，在试验开始时将一定数量的样品都放置在试验箱内，所选的加速应力呈台阶式递增实施。CSADT的研究和应用相对成熟，但由于SSADT比CSADT具有节省样本、试验时间、监测仪器和试验箱的使用，大大减少试验费用等优势，因此目前SSADT受到更为广泛的关注，SSADT的相关研究也越来越多。SSADT方案主要包括试验样本量、试验应力、性能参数的监测记录频率、每一应力下的监测次数等内容。

加速退化试验方案设计的研究文献大约从1994年开始，Boulanger&Escobar两位学者针对一类经典的退化模型对ADT进行了方案设计的研究，讨论了试验应力、总试验样本量、每一应力下样本量比例以及测试时间的确定问题。之后Tseng ST和Yu H F两位学者在1997和1998年分别对退化试验和加速退化试验的截尾时间进行了探讨。2002年Shuo-Jye Wu&Chun-Tao Chang在费用约束条件下以p阶分位寿命的均方误差最小为目标，提出了退化率服从指数分布的退化试验优化设计方法。2004年LC Tang，GY Yang&M Xie针对两应力水平的SSADT提出以试验费用最少为目标，以一定可靠度估计精度为约束，退化量服从随机过程的ADT优化设计方法，确定试验样本量和各应力水平的监测次数。同年，Yu H F&Tseng ST以p阶分位寿命方差最小为目标，在试验费用约束条件下，对退化率服从倒数Weibull分布的退化试验进行了优化设计研究，得到了最优的样本量、监测频率和试验监测次数。2006年Chen-Mao Liao&Tseng ST以试验费用不超过总预算为目标，以产品p分位寿命估计的渐进方差最小为约束，用随机扩散过程对典型SSADT问题建模，来确定试验样本量，监测频率和试验截止时间。2009年Yashun Wang；Chunhua Zhang&Xun Chen提出基于Monte Carlo仿真的SSADT优化方法，以p阶分位寿命渐进方差的局部估计最小为目标，以试验费用为约束，利用Monte Carlo仿真方法模拟退化试验过程生成数据，通过统计分析计算p分位寿命渐进方差，选取渐进方差最小的方案为最优方案，方案给出试验样本量，监测频率和监测次数。Tseng，S T&Balakrishnan，N&Tsai，C C等学者对服从gamma过程的疲劳失效产品进行SSADT优化设计，以试验费用不超过预算为约束，以产品MTTF估计的渐进方差最小为目标，进行试验优化，最终确定试验样本量，监测频率和试验截止时间。Xiaoyang Li&Tongmin Jiang讨论了在多种退化失效模式下制订SSADT方案的方法，根据竞争失效机理，以漂移布朗运动建立多态退化系统的可靠性模型。以试验费用不超过预算为约束，以p分位寿命估计的渐进方差最小为目标，进行试验优化，最优方案给出试验样本，监测频率和监测次数。

到目前，对SSADT方案优化设计的研究中，优化设计变量都没有包含试验应力水平，试验优化在试验应力水平给定的情况下进行，最优方案中只提供试验样本，监测频率，试验截尾时间或监测次数。但是试验应力水平是试验方案中很重要的一项，只有在合适的应力水平下才能更好地激发产品性能退化，获得更为有效的退化数据，从而提高产品可靠性评估精度。

因此，在加速退化试验中将试验应力水平也作为一项优化变量，对试验样本量、试验应力水平、监测频率、各应力下监测次数进行系统的优化设计是目前加速退化试验优化设计领域的当务之急。

发明内容

本发明的目的是为了解决上述问题，提出了将试验应力水平作为优化变量的步进应力加速退化试验(SSADT)优化方法。将试验方案中所涉及到的所有变量，包括试验样本量、试验应力水平、监测频率、各应力下监测次数都考虑在内，以p分位寿命可靠度的渐进方差最小为目标，以试验总费用不超过预算为约束，进行系统的优化设计，使SSADT优化设计方法更为全面完善，优化结果更具有工程实际意义和应用价值。

本发明是步进应力加速退化试验优化设计方法，包括以下几个步骤：

步骤一、确定产品性能退化模型，加速模型和可靠度函数；；

步骤二、建立目标函数；

步骤三、确定约束条件和优化变量；

步骤四、实现优化算法。

本发明的优点在于：

(1)建立目标函数：以p分位寿命可靠度的渐进方差最小为目标。计算p分位寿命可靠度的渐进方差时，采用差分方式对可靠度函数求偏导，该方法解决了复杂函数求导困难，容易出错等问题。

(2)根据工程经验提出对每一应力下监测次数的约束条件，使试验优化结果更具有工程实际意义。由于低应力下产品性能退化慢，高应力下退化快，故低应力下的监测次数应大于高应力下的监测次数，而且根据工程实际情况，各应力下的监测次数应不小于10次，使每一应力下都能采集到足够的性能退化数据，即M₁≥M₂≥……≥M_K≥10，这一约束也大大减少了计算次数，提高了试验优化速率。

(3)将试验应力作为优化变量，并根据实际情况给出应力水平的约束条件，对试验方案更系统更全面地进行优化。在实施SSADT前，应了解到产品的极限应力，即已知正常应力和工作极限应力S_max的情况下，试验应力约束条件可设为S₀＜S₁＜S₂＜……＜S_K≤S_max。试验的优化模型中，目标函数min AsVar(R(ξ_p))是一个关于样本量，监测间隔，监测次数，以及试验应力的函数，将试验应力作为优化变量，可将试验方案进行综合的全面的优化设计。

附图说明

图1是本发明步进应力加速退化试验优化设计方法的流程图；

图2是本发明中步骤四的流程图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种步进应力加速退化试验优化设计方法，流程如图1所示，包括以下几个步骤：

步骤一、确定产品性能退化模型、加速模型和可靠度函数；

根据产品特点、敏感应力和性能参数退化趋势等，确定产品性能退化模型和产品的加速模型，进而确定退化模型的概率密度分布函数f(t，C)和可靠性函数R(t)；

产品的性能退化模型是指对产品性能参数随时间的退化趋势进行拟合的数学描述。目前，性能参数退化趋势拟合方法主要采用对性能参数随时间变化的确定函数中添加随机过程分量描述。性能退化模型主要有布朗(brown)运动、伽马(gamma)过程等，其中漂移布朗运动(带有漂移系数的布朗运动)的应用最多。

加速模型用来描述产品性能参数退化速率与应力之间的关系，是一个仅与应力相关的确定性函数。常用的加速模型有阿伦尼斯(Arrhenius)模型、逆幂率模型、艾琳(Egring)模型等，其形式都可表示为：

其中，

是应力s的某一已知函数，d(s)为性能退化率，A、B为常数。

若产品性能退化模型可采用对性能参数随时间变化的确定函数中添加随机过程分量来描述，加速模型通过变形可表示为

则可采用本发明的方法进行试验优化设计。

例如：当选取漂移布朗运动模型为产品性能退化模型，描述产品的性能退化过程，如式(1)所示：

Y(t)＝σB(t)+d(s)·t+y₀                            (1)

其中：y₀-漂移布朗运动的起始点，即产品性能初始值；

Y(t)-产品性能退化过程，是一个漂移布朗运动；

B(t)-均值为0，方差为时间t的标准布朗运动，B(t)～N(0，t)；

σ-扩散系数，σ＞0，不随应力和时间而改变，是常数；

d(s)-漂移系数，也称为性能退化率；

本发明中选取Arrhenius模型为加速模型，加速应力为温度T，则d(s)为：

d(T_k)＝exp(A+B/T_k)                            (2)

A为常数，B＝-Ea/k，k是波尔兹曼常数8.6171×10-5eV/K，Ea是激活能，单位为eV，T_k是绝对温度，单位为K。

漂移布朗运动模型的概率密度分布函数为：

$f(t,C)=\frac{C}{\mathrm{\σ}\sqrt{2\mathrm{\πt}}}\exp \{-\frac{{[C-d\left(s\right)\·t-y_0]}^2}{2{\mathrm{\σ}}^{2}t}\}---\left(3\right)$

其中：C为产品失效临界值；σ为扩散系数；d(s)为应力s下的性能退化率；y₀为性能初始值。具有形如上式的概率密度函数的分布，通常称为逆高斯分布，则产品的可靠度函数为：

$R\left(t\right)=\mathrm{\Φ}\left[\frac{C-y_0-d\left(s\right)\·t}{\mathrm{\σ}\sqrt{t}}\right]-\exp \left(\frac{2d\left(s\right)\·(C-y_0)}{{\mathrm{\σ}}^2}\right)\mathrm{\Φ}[-\frac{C-y_0+d\left(s\right)\·t}{\mathrm{\σ}\sqrt{t}}]---\left(4\right)$

Φ表示标准正态分布的累积概率分布函数，R(t)即为基于漂移布朗运动的可靠性评估模型。

步骤二、建立目标函数；

选取p阶分位寿命估计值的可靠度的渐进方差最小作为目标函数；即：

min  AsVar(R(ξ_p))

其中，ξ_p是产品在正常条件下的p阶分位寿命的估计值；

具体为：

(1)建立产品性能退化模型的对数似然函数；

设K个应力水平的步进应力加速退化试验，样本量为n，产品性能监测时间间隔为Δt，试验应力为S_k，各应力水平的性能监测次数为M_k，试验中共监测M次，k＝1，2，…，K，每一应力水平的试验时间为M_k·Δt，总试验时间为M·Δt；每次进行监控的时间为t_kij，k＝1，...，K；i＝1，...，n，j＝1，...，M_k，监控到的性能值为y_ikj，性能增量ΔD_ikj＝y_ikj-y_ik(j-1)。由采用性能增量表示概率密度分布函数f(Δt，ΔD)，得到极大似然函数L，进而得到对数似然函数ln L；

(2)求解产品在正常条件下的p阶分位寿命的估计值ξ_p；

由步骤一确定的可靠性函数R(t)，对给定的分位点p，则

是产品在正常条件下的p分位寿命的估计值；

(3)求解渐进方差AsVar(R(ξ_p))；

ξ_p的估计值

在n→∞时，服从均值为ξ_p，方差为h^TI^-1(θ)h，记为AsVar(R(ξ_p))的渐进正态分布：

即，渐进方差为AsVar(R(ξ_p))＝h^TI^-1(θ)h；

其中：θ为可靠度函数R(t)中的未知参数，假设存在，n个未知参数，则θ＝[θ₁，θ₂，…，θ_n]；

$h^T=(\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\θ}}_1},\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\θ}}_2},\·\·\·,\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\θ}}_n})---\left(5\right)$

I(θ)为θ的Fisher信息矩阵，是一个三阶正定方阵；

①采用差分方法求h；

采用差分方法对可靠度函数求偏导，得：

$\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\θ}}_i}=\underset{\mathrm{\Δ}{\mathrm{\θ}}_i\→0}{\lim }\frac{R({\mathrm{\ξ}}_p,{\mathrm{\θ}}_i+\mathrm{\Δ}{\mathrm{\θ}}_i)-R({\mathrm{\ξ}}_p,{\mathrm{\θ}}_i)}{\mathrm{\Δ}{\mathrm{\θ}}_i}---\left(7\right)$

式中：θ_i∈[θ₁，θ₂，…，θ_n，]，Δθ_i趋于0，取值为θ_i的10^-6或更小；

②推导Fisher信息矩阵I(0)；

根据对数似然函数ln L，对参数θ_i求一阶偏导和二阶偏导，求解式(6)，得到信息矩阵I(θ)；

③得到p阶分位寿命可靠度估计的渐进方差，建立目标函数；

AsVar(R(ξ_p))是为关于样本量n，监测间隔Δt，监测次数[M₁，M₂，......，M_K]和试验应力[S₁，S₂，......，S_K]的函数，最后将p阶分位寿命可靠度的渐进方差最小作为目标函数，即min AsVar(R(ξ_p))；

例如，针对步骤一中的漂移布朗运动模型，步骤二具体为：

(1)建立产品性能退化模型的对数似然函数；

将加速模型式(2)带入概率密度函数式(3)，用性能增量表示为

$f(\mathrm{\Δt},\mathrm{\ΔD})=\frac{1}{\mathrm{\σ}\sqrt{2\mathrm{\π\Δt}}}\exp \{-\frac{{[\mathrm{\ΔD}-(A+B/T_k)\·\mathrm{\Δt}]}^2}{2{\mathrm{\σ}}^2\mathrm{\Δt}}---\left(8\right)$

其极大似然函数为：

$L\∝\underset{i=1}{\overset{n}{\mathrm{\Π}}}\underset{k=1}{\overset{K}{\mathrm{\Π}}}\underset{j=1}{\overset{M_k}{\mathrm{\Π}}}\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}^2\mathrm{\Δt}}}\exp \{-\frac{{[{\mathrm{\ΔD}}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}^2}{2{\mathrm{\σ}}^2\mathrm{\Δt}}\}---\left(9\right)$

对数似然函数为：

$\ln L\∝-\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\Π}}}\underset{k=1}{\overset{K}{\mathrm{\Π}}}\underset{j=1}{\overset{M_k}{\mathrm{\Π}}}\left\{\right[\ln \left(2\mathrm{\π\Δt}\right)+\ln \left({\mathrm{\σ}}^2\right)]+\frac{{[{\mathrm{\ΔD}}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}^2}{{\mathrm{\σ}}^2\mathrm{\Δt}}\}---\left(10\right)$

(2)求解产品在正常条件下的p阶分位寿命的估计值ξ_p；

由式(4)，对给定的分位点p，则

是产品在正常条件下的p分位寿命的估计值。

(3)求解渐进方差AsVar(R(ξ_p))；

ξ_p的估计值

在n→∞时，服从均值为ξ_p，方差为h^TI^-1(θ)h(记为AsVar(R_S(ξ_p))的渐进正态分布：

即，渐进方差为AsVar(R(ξ_p))＝h^TI^-1(θ)h；

其中：                              θ＝(A，B，σ²)

$h^T=(\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂A},\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂B},\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\σ}}^2})---\left(11\right)$

$I\left(\mathrm{\θ}\right)=\left[\begin{array}{ccc}E\left(\frac{{\∂}^2\ln L}{\∂A^2}\right)& E(-\frac{{\∂}^2\ln L}{\∂A\∂B})& E\left(\frac{{\∂}^2\ln L}{\∂A\∂{\mathrm{\σ}}^2}\right)\\ & E(-\frac{{\∂}^2\ln L}{\∂B^2})& E(-\frac{{\∂}^2\ln L}{\∂B\∂{\mathrm{\σ}}^2})\\ \mathrm{symmetrical}& & E(-\frac{{\∂}^2\ln L}{\∂{\left({\mathrm{\σ}}^2\right)}^2})\end{array}\right]---\left(12\right)$

I(θ)为θ的Fisher信息矩阵，是一个三阶正定方阵。

①采用差分方法求h；

采用差分方法对可靠度函数求偏导，得：

$\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂A}=\underset{\mathrm{\ΔA}\→0}{\lim }\frac{R({\mathrm{\ξ}}_p,A+\mathrm{\ΔA})-R({\mathrm{\ξ}}_p,A)}{\mathrm{\ΔA}}$

$\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂B}=\underset{\mathrm{\ΔB}\→0}{\lim }\frac{R({\mathrm{\ξ}}_p,B+\mathrm{\ΔB})-R({\mathrm{\ξ}}_p,B)}{\mathrm{\ΔB}}---\left(13\right)$

$\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\σ}}^2}=\underset{\mathrm{\Δ}{\mathrm{\σ}}^2\→0}{\lim }\frac{R({\mathrm{\ξ}}_p,{\mathrm{\σ}}^2+\mathrm{\Δ}{\mathrm{\σ}}^2)-R({\mathrm{\ξ}}_p,{\mathrm{\σ}}^2)}{\mathrm{\Δ}{\mathrm{\σ}}^2}$

由于可靠度函数R(t)的表达式复杂，故对其参数求偏导的解析式难以给出，推导过程极易出错，Chen-Mao Liao&Tseng ST对SSADT优化设计方法的研究中，附录里给出了对产品寿命分布函数求偏导的解析式，在对模型参数a求导时直接套用φ(t，a)即令其他参数b，a，σ²也是如此，但Ф为标准正态分布，对t求导

其他参数不能直接套用φ(t)，故该文献中计算渐进方差有误。

本发明中采用差分方式对可靠度函数求偏导，只要ΔA，ΔB，Δσ²趋于0，取值为参数A，B，σ²的10^-6或更小，则该方法适用。

②推导Fisher信息矩阵I(θ)；

在Δt时间内的性能增量ΔD服从均值为d(s)·Δt，方差为σ²·Δt的正态分布。可知

$E\left(\frac{\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}}{\sqrt{{\mathrm{\σ}}^2\mathrm{\Δt}}}\right)=0---\left(14\right)$

$\mathrm{Var}\left(\frac{\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}}{\sqrt{{\mathrm{\σ}}^2\mathrm{\Δt}}}\right)=1---\left(15\right)$

利用公式(14)、公式(15)求解Fisher信息矩阵，公式推导如下：

$\frac{\∂\ln L}{\∂A}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\{\frac{[\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}{{\mathrm{\σ}}^2}\·\exp (A+B/T_k)\},$

$\frac{\∂\ln L}{\∂B}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\{\frac{[\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}{{\mathrm{\σ}}^2}\·\frac{\exp (A+B/T_k)}{T_k}\},$

$\frac{\∂\ln L}{\∂{\mathrm{\σ}}^2}=-\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\Π}}}\underset{k=1}{\overset{K}{\mathrm{\Π}}}\underset{j=1}{\overset{M_k}{\mathrm{\Π}}}\{\frac{1}{{\mathrm{\σ}}^2}-\frac{{[{\mathrm{\ΔD}}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}^2}{{\mathrm{\σ}}^4\mathrm{\Δt}}.\}$

1)

的求解为：

$\frac{{\∂}^2\ln L}{\∂A^2}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\{\frac{[\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}{{\mathrm{\σ}}^2\mathrm{\Δt}}\·\exp (A+B/T_k)-\frac{\exp (2A+2B/T_k)\·\mathrm{\Δt}}{{\mathrm{\σ}}^2}\},$

$\⇒E(-\frac{{\∂}^2\ln L}{\∂A^2})=n\·\frac{\mathrm{\Δt}}{{\mathrm{\σ}}^2}\·\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\·\exp \left[2(A+B/T_k)\right];$

2)的求解为：

$\frac{{\∂}^2\ln L}{\∂A\∂B}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\left\{\frac{\exp (A+B/T_k)}{T_k}\right[\frac{\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}}{{\mathrm{\σ}}^2}]-\frac{\exp (2A+2B/T_k)\·\mathrm{\Δt}}{{\mathrm{\σ}}^2{T}_k}\},$

$\⇒E(-\frac{{\∂}^2\ln L}{\∂A\∂B})=n\·\frac{\mathrm{\Δt}}{{\mathrm{\σ}}^2}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\·\frac{\exp (2A+2B/T_k)}{T_k};$

3)的求解为：

$\frac{{\∂}^2\ln L}{\∂B^2}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\{\frac{[\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}{{\mathrm{\σ}}^2\mathrm{\Δt}}\·\frac{\exp (A+B/T_k)}{T_k^2}-\frac{\exp (2A+B/T_k)\·\mathrm{\Δt}}{{\mathrm{\σ}}^2{T}_k^2}\},$

$\⇒E(-\frac{{\∂}^2\ln L}{\∂B^2})=n\·\frac{\mathrm{\Δt}}{{\mathrm{\σ}}^2}\·\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\·\frac{\exp (2A+2B/T_k)}{T_k^2};$

4)

的求解为：

$E(-\frac{{\∂}^2\ln L}{\∂B\∂A})=E(-\frac{{\∂}^2\ln L}{\∂A\∂B});$

5)

的求解为：

$\frac{{\∂}^2\ln L}{\∂{\left({\mathrm{\σ}}^2\right)}^2}=-\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\{-\frac{1}{{\mathrm{\σ}}^4}+\frac{2\×{[\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}^2}{{\mathrm{\σ}}^6\mathrm{\Δt}}\}=-\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\{-\frac{1}{{\mathrm{\σ}}^4}+0\},$

$\⇒E(-\frac{{\∂}^2\ln L}{\∂{\left({\mathrm{\σ}}^2\right)}^2})=\frac{n}{2{\mathrm{\σ}}^4}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k;$

6)

的求解为：

$\frac{{\∂}^2\ln L}{\∂A\∂{\mathrm{\σ}}^2}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\{-\frac{[{\mathrm{\ΔD}}_{1\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}{{\mathrm{\σ}}^2}\·\frac{1}{{\mathrm{\σ}}^2}\exp (A+B/T_k)\},$

$\⇒E(-\frac{{\∂}^2\ln L}{\∂A\∂{\mathrm{\σ}}^2})=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\{\frac{\mathrm{\Δt}}{{\mathrm{\σ}}^2}\·\exp (A+B/T_k)\·0\}=0;$

由上式可见，除了方差自身的二阶偏导数，一切与方差相关的偏导数均等于零，故有：

$E(-\frac{{\∂}^2\ln L}{\∂{\mathrm{\σ}}^2\∂A})=0,$ $E(-\frac{{\∂}^2\ln L}{\∂B\∂{\mathrm{\σ}}^2})=0,$ $E(-\frac{{\∂}^2\ln L}{\∂{\mathrm{\σ}}^2\∂B})=0.$

由以上可知，则Fisher信息矩阵为

$I\left(\mathrm{\θ}\right)=$

$\left[\begin{array}{ccc}n\·\frac{\mathrm{\Δt}}{{\mathrm{\σ}}^2}\·\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\·\exp (2A+2B/T_k)& n\·\frac{\mathrm{\Δt}}{{\mathrm{\σ}}^2}\·\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\·\frac{\exp (2A+2B/T_k)}{T_k}& 0\\ n\·\frac{\mathrm{\Δt}}{{\mathrm{\σ}}^2}\·\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\·\frac{\exp (2A+2B/T_k)}{T_k}& n\·\frac{\mathrm{\Δt}}{{\mathrm{\σ}}^2}\·\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\·\frac{\exp (2A+2B/T_k)}{T_k^2}& 0\\ 0& 0& \frac{n}{2{\mathrm{\σ}}^4}\·\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\end{array}\right],---\left(16\right)$

③得到p阶分位寿命可靠度估计的渐进方差，建立目标函数；

$\mathrm{AsVar}\left(R\left({\mathrm{\ξ}}_p\right)\right)=h^T{I}^{-1}\left(\mathrm{\θ}\right)h=$

$(\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂A},\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂B},\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\σ}}^2})\×{\left[\begin{array}{ccc}E(-\frac{{\∂}^2\ln L}{\∂A^2})& E(-\frac{{\∂}^2\ln L}{\∂A\∂B})& 0\\ E(-\frac{{\∂}^2\ln L}{\∂B\∂A})& E(-\frac{{\∂}^2\ln L}{\∂B^2})& 0\\ 0& 0& E(-\frac{{\∂}^2\ln L}{\∂{\left({\mathrm{\σ}}^2\right)}^2})\end{array}\right]}^{-1}\×\left(\begin{array}{c}\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂A}\\ \frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂B}\\ \frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\σ}}^2}\end{array}\right)---\left(17\right)$

将p阶分位寿命可靠度的渐进方差最小作为目标函数，即min AsVar(R(ξ_P))。

步骤三、确定约束条件和优化变量；

①确定约束条件：

(1)试验总费用C_t；试验费用的主要因素包括两部分：

a)单位时间内的试验费用C_o，包括试验设备折旧、产品性能测试仪器的使用、试验人员的工资以及资源(水、电等)消耗等。单位时间：元/小时；

b)试验样本单价C_d。

因此试验总费用可以表示为：

$C_t=n\·C_d+\mathrm{\Δt}\·\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\·C_o$

(2)监测次数M₁≥M₂≥……≥M_K≥10；

根据工程经验每一应力下的性能检测次数应不小于10次，若某一应力下监测次数过少，则在评估时该应力下退化率的回归拟合精度难以确保；在失效机理不变的前提下，产品退化率随应力的升高而加快，为使每一应力下能够检测到足够的退化量来外推产品的寿命与可靠性，应使低应力水平下的监测次数应高于高应力水平。

(3)应力水平：S₀＜S₁＜S₂＜……＜S_K≤S_max；

S₁应高于产品使用时的正常应力S₀，S_max为产品的工作应力极限，通过预试验或者产品的设计和使用说明等相关信息可以获得S_max，进而根据实际情况确定最高应力S_K，各应力下产品的失效机理不变。

②确定优化变量；

(1)样本量n；

(2)监测时间间隔Δt；

(3)监测次数M_k，(k＝1，2，…K)；

(4)应力水平S_k，(k＝1，2，…K)。

即步进应力加速退化试验共投入n个样本，对每一样本每隔Δt小时监测一次，进行K步应力水平试验，每应力水平S_k下，监测次数为M_k。

步骤四、实现优化算法；

由步骤一至步骤三，确定SSADT试验方案的优化问题，用数学描述为：

min  AsVar(R(ξ_p))

s.t. $n\·C_d+\mathrm{\Δt}\·\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\·C_o\≤C_t---\left(18\right)$

M₁≥M₂≥……≥M_K≥10

S₀＜S₁＜S₂＜……＜S_K≤S_max

上式为本发明中的优化模型。其中，需要优化的决策变量为n，Δt，M_k，S_k，它们的取值均为正整数。

优化算法具体为：

首先根据工程经验确定退化模型式中的参数[θ₁，θ₂，…，θ_n，]、产品性能初值y₀、产品退化临界值C、单位时间内的试验费用C_o、试验样本单价C_d、试验总费用C_t、产品试验正常应力S₀、最高应力S_K、试验应力水平数K；然后根据步骤三所述的约束条件，确定所有满足约束条件的试验应力[S₁，S₂，…，S_K]和监测次数[M₁，M₂，…，M_K]，计算每个满足约束条件时相应的p阶分位寿命可靠度的渐进方差，最后选取方差最小的试验方案作为最优方案。优化算法的流程如图2所示，包括以下几个步骤：

Step1：根据历史数据、相似产品信息、工程经验、手册、标准等与试验样本相关的信息确定初始参数，包括：产品性能退化模型的参数[θ₁，θ₂，…，θ_n，]；产品性能初值y₀；产品退化临界值C；单位时间内的试验费用S_o；试验样本单价C_d；试验总费用C_t；产品试验正常应力S₀、最高应力S_K、试验应力水平数K。

Step2：计算p阶分位寿命ξ_p；

根据试验最高应力S_K和试验应力水平数K，以及应力水平的约束条件S₀＜S₁＜S₂＜……＜S_K≤S_max，求解[S₁，S₂，…，S_K]的所有组合，记为解空间Ω_S。

例如，当S_K表示温度时，令S₁从50℃开始(正常应力S₀＝25℃)，并令各应力之间间隔为5的倍数，当S_K确定后，S₁，S₂，…，S_K-1的取值空间为{50，55，60，......，S_K-5}，以S₁＜S₂＜……＜S_K-1＜S_K为约束，可得到[S₁，S₂，…，S_K-1，S_K]的解空间Ω_S。

以K＝3，最高应力S₃＝70为例，S₁，S₂的取值空间为{50，55，60，65}，以S₁＜S₂＜S₃为约束，[S₁，S₂，S₃]的解空间为{[50，55，70]，[50，60，70]，[50，65，70]，[55，60，70]，[55，65，70]，[60，65，70]}共6个解。

Step3：令样本量n从3开始(一般样本量不宜少于3)，由试验费用约束条件

为求n的最大值，若令M_k取最小值10，Δt＝1，则n的上限为n_max＝(C_t-C_o·K·10)/C_d。

Step4：若n≤n_max继续下一步，否则，跳到Step14。

Step5：令监测间隔Δt＝1，根据试验费用约束条件

若令M_k取最小值，即M_k＝10(k＝1，…，K)则可得到Δt的最大值，其上限

Step6：若Δt≤Δt_max继续下一步，否则，跳到Step12。

Step7：计算

根据约束条件M₁≥M₂≥……≥M_K≥10，则若

继续下一步，否则跳到Step12。

Step8：由监测次数的约束条件M₁≥M₂≥……≥M_K≥10和Step7计算得到的

求解[M₁，M₂，……，M_K]的所有组合，记为解空间Ω_M。

以K＝3，

为例，根据约束条件M₁≥M₂≥M₃≥10，[M₁，M₂，M₃]的解空间为：{[15，10，10]，[14，11，10]，[13，12，10]，[13，11，11]，[12，12，11]}共5个解。

Step9：由Step8求得的解空间Ω_M和Step2得到的Ω_S，两个解空间的所有组合Ω_M×Ω_S，按Step2和Step8中的例子，则监测次数[M₁，M₂，......，M_K]和试验应力[S₁，S₂，…，S_K]的所有组合共6×5＝30个。根据步骤二目标函数的建立过程，AsVar(R(ξ_p))是一个关于样本量n，监测间隔Δt，监测次数[M₁，M₂，......，M_K]和试验应力[S₁，S₂，…，S_K]的函数，将所有组合带入Fisher信息矩阵I(θ)，计算渐进方差AsVar(R(ξ_P))＝h^TI^-1(θ)h，记为：

g(AsVarR，n，Δt，[M₁，M₂，……，M_K]，[S₁，S₂，…，S_K])(其中AsVarR为求得的渐进方差)；

Step10：取g中渐进方差最小，记为

G(Δt，：)＝ming(AsVarR，n，Δt，[M₁，M₂，……，M_K]，[S₁，S₂，…，S_K])；

Step11：Δt＝Δt+1，重复循环Step6～Step10；

Step12：若Δt＞Δt_max，取G(Δt，：)中渐进方差最小，记为Gn(n，：)＝minG(Δt，：)；

Step13：令n＝n+1，重复循环Step4～Step12；

Step14：若n＞n_max，取Gn(n，：)中渐进方差最小，记为minGn＝minGn(n，：)，即为达到目标函数渐进方差最小的最优试验。

实施例：

步骤一、确定产品性能退化模型，加速模型和可靠度函数；

根据产品特点、敏感应力和性能参数退化趋势，选用漂移布朗运动模型描述产品性能退化过程，试验应力为温度T，对应的加速模型选用Arrhenius模型。则确定产品性能退化模型为式(1)。

产品在正常应力下的可靠度函数为：

$R\left(t\right)=\mathrm{\Φ}\left[\frac{C-y_0-(A+B/T_0)\·t}{\mathrm{\σ}\sqrt{t}}\right]-\exp \left(\frac{2(A+B/T_0)\·(C-y_0)}{{\mathrm{\σ}}^2}\right)\mathrm{\Φ}[-\frac{C-y_0+(A+B/T_0)\·t}{\mathrm{\σ}\sqrt{t}}]$

步骤二、建立目标函数；

建立目标函数：min  AsVar(R(ξ_p))

p阶分位寿命可靠度渐进方差AsVar(R(ξ_p))的计算过程见具体实施方式中的步骤二。

步骤三、确定约束条件和优化变量；

①确定约束条件：

(1)试验总费用C_t＝100×10³元；试验费用的主要因素包括两部分：

a)单位时间内的试验费用C_o，包括试验设备折旧、产品性能测试仪器的使用、试验人员的工资以及资源(水、电等)消耗等。单位时间：元/小时；

C_o＝0.6×10³元/小时

b)试验样本单价C_d

C_d＝3×10³元/小时。

试验总费用表示为：

$\mathrm{Ct}=n\·C_d+\mathrm{\Δt}\·\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\·C_o$

(2)监测次数M₁≥M₂≥……≥M_K≥10。

(3)应力水平T₀＜T₁＜T₂＜……＜T_K≤T_max。

②确定优化变量

(1)样本量n；

(2)监测时间间隔Δt；

(3)监测次数M_k，(k＝1，2，…K)；

(4)应力水平为温度T_k，(k＝1，2，…K)。

步骤四、实现优化算法。

Step1：确定初始参数：

表1退化模型参数信息表

模型参数	A	激活能Ea	扩散系数σ	性能初始值y0	性能临界值C
						估计值	11.44	0.5	0.07	100	95

B＝-E_a/k＝-5802.4；

表2试验费用表

费用构成	试验样本的单价Cd×103元/个	试验时间单价Co×103元/小时
			单价	3	0.6

试验总费用C_t＝100×10³元；

试验应力水平	正常应力	最高应力
			摄氏温度℃	25	105
绝对温度K	298.15	378.15

选应力水平数K＝4；

Step2：计算p阶分位寿命ξ_p：p＝0.5，ξ_0.5＝R^-1(0.5)＝6349；

寻找所有符合要求的[T₁，T₂，……，T_K-1，T_K]的解空间Ω_T，经计算Ω_T中共165组解。

Step3：令样本量n从3开始，其上限为n_max＝(C_t-C_o·K·10)/C_d＝25；

Step4到Step14按照优化算法的具体步骤实施，求出产品p分位寿命可靠度渐进方差最小时对应的最优试验。

运算结果如下：

试验总成本上限Ct×103	样本量n	监测间隔Δt(h)	试验应力[T1，T2，T3，T4]	监测次数[M1，M2，M3，M4]	试验实际成本Ct产103	AsVar(R(ξ0.5))％
							100	17	1	[70，75，80，105]	[28，28，13，12]	99.6	3.0687

Claims (4)

1.步进应力加速退化试验优化设计方法，其特征在于，包括以下几个步骤：

步骤一、确定产品性能退化模型，加速模型和可靠度函数；

确定产品性能退化模型和产品的加速模型，进而确定退化模型的概率密度分布函数f(t，C)和可靠性函数R(t)；

步骤二、建立目标函数；

选取p阶分位寿命估计值的可靠度的渐进方差最小作为目标函数；即：

min  AsVar(R(ξ_p))(1)

其中，ξ_p是产品在正常条件下的p阶分位寿命的估计值；

具体为：

(1)建立产品性能退化模型的对数似然函数；

设K个应力水平的步进应力加速退化试验，样本量为n，产品性能监测时间间隔为Δt，试验应力为S_k，各应力水平的性能监测次数为M_k，试验中共监测M次，

k＝1，2，…，K，每一应力水平的试验时间为M_k·Δt，总试验时间为M·Δt，每次进行监控的时间为t_kij，k＝1，...，K；i＝1，...，n，j＝1，...，M_k，监控到的性能值为y_ikj，性能增量ΔD_ikj＝y_ikj-y_ik(j-1)，由采用性能增量表示概率密度分布函数f(Δt，ΔD)，得到极大似然函数L，进而得到对数似然函数lnL；

(2)求解产品在正常条件下的p阶分位寿命的估计值ξ_p；

由步骤一确定的可靠性函数R(t)，对给定的分位点p，则

是产品在正常条件下的p分位寿命的估计值；

(3)求解渐进方差AsVar(R(ξ_p))；

ξ_p的估计值

在n→∞时，服从均值为ξ_p，方差为h^TI^-1(θ)h，记为AsVar(R(ξ_p))的渐进正态分布：

即，渐进方差为AsVar(R(ξ_p))＝h^TI^-1(θ)h；

其中：θ为可靠度函数R(t)中的未知参数，假设存在n个未知参数，则θ＝[θ₁，θ₂，…，θ_n，]；

$h^T=(\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\θ}}_1},\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\θ}}_2},...,\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\θ}}_n})---\left(2\right)$

I(θ)为θ的Fisher信息矩阵，是一个三阶正定方阵；

①采用差分方法求h；

采用差分方法对可靠度函数求偏导，得：

$\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\θ}}_i}=\underset{\mathrm{\Δ}{\mathrm{\θ}}_i\→0}{\lim }\frac{R({\mathrm{\ξ}}_p,{\mathrm{\θ}}_i+{\mathrm{\Δ\θ}}_i)-R({\mathrm{\ξ}}_p,{\mathrm{\θ}}_i)}{{\mathrm{\Δ\θ}}_i}---\left(4\right)$

式中：θ_i∈[θ₁，θ₂，…，θ_n，]，Δθ_i趋于0，取值为θ_i的10^-6或更小；

②推导Fisher信息矩阵I(θ)；

根据对数似然函数lnL，对参数θ_i求一阶偏导和二阶偏导，求解式(3)，得到信息矩阵I(θ)；

③得到p阶分位寿命可靠度估计的渐进方差，建立目标函数；

AsVar(R(ξ_p))是为关于样本量n，监测间隔Δt，监测次数[M₁，M₂，......，M_K]和试验应力[S₁，S₂，......，S_K]的函数，最后将p阶分位寿命可靠度的渐进方差最小作为目标函数，即min AsVar(R(ξ_p))；

步骤三、确定约束条件和优化变量；

①确定约束条件：

(1)试验总费用C_t包括：

a)单位时间内的试验费用C_o，单位时间：元/小时；

b)试验样本单价C_d；

试验总费用C_t表示为：

$C_t=n\·C_d+\mathrm{\Δt}\·\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\·C_o---\left(5\right)$

(2)监测次数：M₁≥M₂≥……≥M_K≥10；

每一个应力下的性能检测次数大于等于10次；

(3)应力水平：S₀＜S₁＜S₂＜……＜S_K≤S_max；

其中，S₁至S_K表示K个应力水平，S₀为产品使用时的正常应力，S_max为产品的工作应力极限；

②确定优化变量；

(1)样本量n；

(2)监测时间间隔Δt；

(3)监测次数M_k，(k＝1，2，…K)；

(4)应力水平S_k，(k＝1，2，…K)；

即步进应力加速退化试验共投入n个样本，对每一样本每隔Δt小时监测一次，进行K步应力水平试验，每应力水平S_k下，监测次数为M_k；

步骤四、实现优化算法；

由步骤一至步骤三，确定步进应力加速退化试验方案的优化问题，用数学描述为：

min  AsVar(R(ξ_p))

$s.t.n\·C_d+\mathrm{\Δt}\·\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{M}_k\·C_o\≤C_t---\left(6\right)$

M₁≥M₂≥……≥M_K≥10

S₀＜S₁＜S₂＜……＜S_K≤S_max

上式为本发明中的优化模型，其中，需要优化的决策变量为n，Δt，M_k，S_k，它们的取值均为正整数；

具体为：

首先根据工程经验确定退化模型式中的参数[θ₁，θ₂，…，θ_n]、产品性能初值y₀、产品退化临界值C、单位时间内的试验费用C_o、试验样本单价C_d、试验总费用C_t、产品试验正常应力S₀、最高应力S_K、试验应力水平数K；然后根据步骤三所述的约束条件，确定所有满足约束条件的试验应力[S₁，S₂，…，S_K]和监测次数[M₁，M₂，…，M_K]，计算每个满足约束条件时相应的p阶分位寿命可靠度的渐进方差，最后选取方差最小的试验方案作为最优方案。

2.根据权利要求1所述的步进应力加速退化试验优化设计方法，其特征在于，当选取漂移布朗运动模型为产品性能退化模型，所述的步骤一具体为：

产品的退化模型，如式(7)所示：

Y(t)＝σB(t)+d(s)·t+y₀(7)

其中：y₀-漂移布朗运动的起始点，即产品性能初始值；

Y(t)-产品性能退化过程，是一个漂移布朗运动；

B(t)-均值为0，方差为时间t的标准布朗运动，B(t)～N(0，t)；

σ-扩散系数，σ＞0，不随应力和时间而改变，是常数；

d(s)-漂移系数，也称为性能退化率，它是一个仅与应力相关的确定性函数，是加速模型；

假定d(s)为Arrhenius模型，加速应力为温度T，则d(s)为：

d(T_k)＝exp(A+B/T_k)(8)

A为常数，B＝-Ea/k，k是波尔兹曼常数8.6171×10-5eV/K，Ea是激活能，单位为eV，T_k是绝对温度，单位为K；

漂移布朗运动模型的概率密度分布函数为：

$f(t,C)=\frac{C}{\mathrm{\σ}\sqrt{2\mathrm{\πt}}}\exp \{-\frac{{[C-d\left(s\right)\·t-y_0]}^2}{2{\mathrm{\σ}}^{2}t}\}---\left(9\right)$

其中：C为产品失效临界值，σ为扩散系数，d(s)为应力s下的性能退化率，y₀为性能初始值，具有形如上式的概率密度函数的分布，称为逆高斯分布，则产品的可靠度函数为：

$R\left(t\right)=\mathrm{\Φ}\left[\frac{C-y_0-d\left(s\right)\·t}{\mathrm{\σ}\sqrt{t}}\right]-\exp \left(\frac{2d\left(s\right)\·(C-y_0)}{{\mathrm{\σ}}^2}\right)\mathrm{\Φ}[-\frac{C-y_0+d\left(s\right)\·t}{\mathrm{\σ}\sqrt{t}}]---\left(10\right)$

Φ表示标准正态分布的累积概率分布函数，R(t)即为基于漂移布朗运动的可靠性评估模型。

3.根据权利要求1所述的步进应力加速退化试验优化设计方法，其特征在于，当选取漂移布朗运动模型为产品性能退化模型，所述的步骤二为：

根据步骤一中的可靠性函数R(t)，选取p阶分位寿命估计值的可靠度的渐进方差最小作为目标函数即式(1)；

具体为：

(1)建立产品性能退化模型的对数似然函数；

将加速模型式(8)带入概率密度函数式(9)，用性能增量表示为

$f(\mathrm{\Δt},\mathrm{\ΔD})=\frac{1}{\mathrm{\σ}\sqrt{2\mathrm{\π\Δt}}}\exp \{-\frac{{[\mathrm{\ΔD}-(A+B/T_k)\·\mathrm{\Δt}]}^2}{2{\mathrm{\σ}}^2\·\mathrm{\Δt}}\}---\left(11\right)$

其极大似然函数为：

$L\∝\underset{i=1}{\overset{n}{\mathrm{\Π}}}\underset{k=1}{\overset{K}{\mathrm{\Π}}}\underset{j=1}{\overset{M_k}{\mathrm{\Π}}}\frac{1}{\sqrt{{2\mathrm{\π\σ}}^2\mathrm{\Δt}}}\exp \{-\frac{{[\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}^2}{2{\mathrm{\σ}}^2\mathrm{\Δt}}\}---\left(12\right)$

对数似然函数为：

$\ln L\∝-\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\left\{\right[ln\left(2\mathrm{\π\Δt}\right)+\ln \left({\mathrm{\σ}}^2\right)]+\frac{{[\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}^2}{{\mathrm{\σ}}^2\mathrm{\Δt}}\}---\left(13\right)$

(2)求解产品在正常条件下的p阶分位寿命的估计值ξ_p；

由式(10)，对给定的分位点p，则

是产品在正常条件下的p分位寿命的估计值；

(3)求解渐进方差AsVar(R(ξ_p))；

ξ_p的估计值

在n→∞时，服从均值为ξ_p，方差为h^TI^-1(θ)h，记为AsVar(R_S(ξ_p))的渐进正态分布：

则渐进方差为AsVar(R(ξ_p))＝h^TI^-1(θ)h；

其中：θ＝(A，B，σ²)

$h^T=(\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂A},\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂B},\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\σ}}^2})---\left(14\right)$

I(θ)为θ的Fisher信息矩阵，是一个三阶正定方阵；

①采用差分方法求h；

采用差分方式对可靠度函数求偏导，得：

$\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂A}=\underset{\mathrm{\ΔA}\→0}{\lim }\frac{R({\mathrm{\ξ}}_p,A+\mathrm{\ΔA})-R({\mathrm{\ξ}}_p,A)}{\mathrm{\ΔA}}$

$\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂B}=\underset{\mathrm{\ΔB}\→0}{\lim }\frac{R({\mathrm{\ξ}}_p,B+\mathrm{\ΔB})-R({\mathrm{\ξ}}_p,B)}{\mathrm{\ΔB}}---\left(16\right)$

$\frac{\∂R\left({\mathrm{\ξ}}_p\right)}{\∂{\mathrm{\σ}}^2}=\underset{\mathrm{\Δ}{\mathrm{\σ}}^2\→0}{\lim }\frac{R({\mathrm{\ξ}}_p,{\mathrm{\σ}}^2+\mathrm{\Δ}{\mathrm{\σ}}^2)-R({\mathrm{\ξ}}_p,{\mathrm{\σ}}^2)}{\mathrm{\Δ}{\mathrm{\σ}}^2}$

②推导Fisher信息矩阵I(θ)；

在Δt时间内的性能增量ΔD服从均值为d(s)·Δt，方差为σ²·Δt的正态分布，则

$E\left(\frac{{\mathrm{\ΔD}}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}}{\sqrt{{\mathrm{\σ}}^2\mathrm{\Δt}}}\right)=0---\left(17\right)$

$\mathrm{Var}\left(\frac{{\mathrm{\ΔD}}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}}{\sqrt{{\mathrm{\σ}}^2\mathrm{\Δt}}}\right)=1---\left(18\right)$

利用公式(17)、公式(18)求解Fisher信息矩阵，具体过程中如下：

$\frac{\∂\ln L}{\∂A}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\{\frac{[\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}{{\mathrm{\σ}}^2}\·exp(A+B/T_k)\},$

$\frac{\∂\ln L}{\∂B}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\{\frac{[\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}{{\mathrm{\σ}}^2}\·\frac{\exp (A+B/T_k)}{T_k}\},$

$\frac{\∂\ln L}{\∂{\mathrm{\σ}}^2}=-\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{M_k}{\mathrm{\Σ}}}\{\frac{1}{{\mathrm{\σ}}^2}-\frac{{[\mathrm{\Δ}{D}_{\mathrm{ikj}}-\exp (A+B/T_k)\·\mathrm{\Δt}]}^2}{{\mathrm{\σ}}^4\mathrm{\Δt}}\·\}$

1)

的求解为：

2)

的求解为：

的求解为：

的求解为：

的求解为：

的求解为：

由上式，除了方差自身的二阶偏导数，一切与方差相关的偏导数均等于零，故有：

$E(-\frac{{\∂}^2\ln L}{\∂{\mathrm{\σ}}^2\∂A})=0,E(-\frac{{\∂}^2\ln L}{\∂B\∂{\mathrm{\σ}}^2})=0,E(-\frac{{\∂}^2\ln L}{\∂{\mathrm{\σ}}^2\∂B})=0;$

由以上，则Fisher信息矩阵为

③得到p阶分位寿命可靠度估计的渐进方差，建立目标函数；

将p阶分位寿命可靠度的渐进方差最小作为目标函数，即min AsVar(R(ξ_p))。

4.根据权利要求1所述的步进应力加速退化试验优化设计方法，其特征在于，所述的步骤四具体包括以下几个步骤：

Step1：根据工程经验确定初始参数，包括：产品性能退化模型的参数[θ₁，θ₂，…，θ_n，]、产品性能初值y₀、产品退化临界值C、单位时间内的试验费用C_o、试验样本单价C_d、试验总费用C_t、产品试验正常应力S₀、最高应力S_K、试验应力水平数K；

Step2：计算p阶分位寿命ξ_p；

根据试验最高应力S_K和试验应力水平数K，以及应力水平的约束条件S₀＜S₁＜S₂＜……＜S_K≤S_max，得到[S₁，S₂，…，S_K]的所有组合，记为解空间Ω_S；

Step3：令样本量n从3开始，由试验费用约束条件

求n的最大值，令M_k＝10，Δt＝1，则n的上限为n_max＝(C_t-C_o·K·10)/C_d；

Step4：若n≤n_max继续下一步，否则，跳到Step14，其中n_max为总样本量；

Step5：令监测间隔Δt＝1，根据试验费用约束条件若令M_k取最小值10，得到Δt的最大值，

Step6：若Δt≤Δt_max继续下一步，否则，跳到Step12；

Step7：计算

根据约束条件M₁≥M₂≥……≥M_K≥10，若

则继续下一步，否则跳到Step12；

Step8：由监测次数的约束条件M₁≥M₂≥……≥M_K≥10和Step7计算得到的

求解[M₁，M₂，……，M_K]的所有组合，记为解空间Ω_M；

Step9：由Step8求得的解空间Ω_M和Step2得到的Ω_S，两个解空间的所有组合Ω_M×Ω_S，将所有组合带入Fisher信息矩阵I(θ)，计算渐进方差AsVar(R(ξ_p))＝h^TI^-1(θ)h，记为：g(AsVarR，n，Δt，[M₁，M₂，……，M_K]，[S₁，S₂，…，S_K])；

Step10：取g中渐进方差最小，记为G(Δt，：)＝ming(AsVarR，n，Δt，[M₁，M₂，……，M_K]，[S₁，S₂，…，S_K])；

Step11：Δt＝Δt+1，重复循环Step6～Step10；

Step12：若Δt＞Δt_max，取G(Δt，：)中渐进方差最小，记为Gn(n，：)＝minG(Δt，：)；

Step13：令n＝n+1，重复循环Step4～Step12；

Step14：若n＞n_max，取Gn(n，：)中渐进方差最小，记为minGn＝minGn(n，：)，则Gn(n，：)即为达到目标函数渐进方差最小的最优试验。

Images