一种多应力多退化量步进加速退化试验方案优化设计方法

技术领域

本发明涉及多应力多退化量步进加速退化试验方案优化设计方法，属于可靠性工程技术 领域。

背景技术

对于高可靠长寿命产品而言，其寿命与可靠性如果通过传统的寿命试验技术或自然条件 试验技术进行预测，则往往难以在可行的时间内完成。即使采用加速寿命试验技术，也很可 能出现零失效的情况，给产品寿命预测带来困难。由于大部分产品在工作或贮存过程中性能 指标(其函数被定义为退化量)会随时间逐渐退化，如果充分合理的利用这些性能退化数据， 产品的寿命预测会更加高效与准确。加速退化试验对产品施加超出使用应力水平的加速应力， 在退化机理不变的条件下，分析产品在加速应力水平下的性能退化数据，外推得出产品在使 用条件下的寿命与可靠度。如果加速应力随时间成台阶状逐步提高，则称为步进加速退化试 验，它具有试验效率高的优点。加速退化试验技术在可靠性试验工程领域得到越来越广泛的 应用。

如何设计加速退化试验方案使产品寿命预测结果最准确、代价最小，是加速退化试验工 程应用面临的核心问题之一，即加速退化试验方案优化设计问题。针对这一问题，目前出现 的解决方案要么仅适用于单一应力加速退化试验，要么仅适用于单一退化量产品的加速退化 试验。然而，在工程实际中，产品正常工作通常受到多种应力的作用，包括工作应力(如电 流、电压等)和环境应力(如温度、湿度、振动等)，因此单一应力不能真实体现产品实际的 工作应力与环境应力特征。此外，对于高可靠长寿命产品而言，必须将多种应力作为加速应 力才能得到更大的加速系数，同时保证加速退化机理的不变性，这是单一应力无法满足的。

同时，表征产品性能的退化量通常有多个，要完整衡量产品的性能状态需要借助多个退 化量。例如，陀螺在长期贮存过程中的退化量有X方向漂移量、Y方向漂移量、Z方向漂移 量等；火箭管路安全阀的退化量有其导阀调整弹簧松弛应力、主阀复位弹簧松弛应力、密封 圈的永久压缩变形、壳体裂纹长度等。需要将多个退化量综合考虑，才能对其寿命及可靠性 进行正确的建模和分析。

如果利用多应力步进加速退化试验对此类产品进行寿命预测，必然会面临多应力多退化 量场合的复杂步进加速退化试验方案优化设计问题，其设计变量多、类型多样，如样本量(离 散型)、各加速应力的应力水平(连续型)、每一应力水平下的监测时间间隔(离散型)、监测 次数(离散型)等，各设计变量与多退化量所对应的多种退化过程相互作用，对产品寿命或 可靠性估计精度的影响非常复杂，是可靠性工程领域亟待解决的难题。

发明内容

本发明的目的是，提供一种多应力多退化量步进加速退化试验方案优化设计方法，能够 由产品相关先验信息得到费用约束条件下的优化试验方案，使得产品寿命预测结果精度最好。 本发明针对优化目标函数的解析形式难以推导的难题，基于蒙特卡罗统计仿真理论，建立多 应力多退化量的步进加速退化试验方案优化设计优化模型，并提出相应的优化算法解决此优 化问题，最终提出多应力多退化量步进加速退化试验方案优化设计方法。本发明所提出的方 法易于流程化、便于工程应用，可为多应力多退化量场合产品寿命预测提供优化的试验方案 支撑，以最小的试验代价实现最准确的寿命预测。

为实现上述目的，本发明所采取的技术方案为：

一种多应力多退化量步进加速退化试验方案优化设计方法，包括以下步骤：

步骤1、获取产品加速退化试验相关信息

1-1)产品的退化量及失效阈值信息

产品在工作或贮存过程中有m个退化量Y_i(i＝1,2,...,m)随时间逐渐退化，一旦某个退化 量Y_i超过失效阈值D_i(i＝1,2,...,m)，产品就会发生失效。

1-2)产品退化量的联合概率密度函数信息

时刻t产品退化量Y＝(Y₁,Y₂,…,Y_m)^T服从多维正态分布，其联合概率密度函数可表示为

$f\left(Y\right)=\frac{1}{{\left(2\pi \right)}^{m/2}|\Sigma {|}^{1/2}}\exp (-\frac{1}{2}{\left(Y-\mu \right)}^T{\Sigma }^{-1}(Y-\mu \left)\right)---\left(1\right)$

其中，μ＝(μ₁,μ₂,…,μ_m)^T为均值向量，Σ为方差协方差矩阵

σ_ij(i＝1,...,m，j＝1,...,m)为退化量Y_i与Y_j的协方差。当i＝j时，σ_ij为退化量Y_i的方差。|Σ|为Σ 的行列式值。如前所述，步进加速退化试验要求产品退化失效机理不发生改变，此时Σ一般 不随应力水平组合变化。

1-3)产品退化模型及加速模型信息

在不同应力水平组合下，其中表示第i种加速应力的第li个水平， i＝1,...,s，li＝1,...,L，L为应力水平的个数，s为加速应力的个数；产品的m维正态分布均值向 量μ的第j维元素与试验时间的关系满足如下退化模型

μ_j＝b_j+a_jtj＝1,2,…,m(3)

式中，b_j为截距参数，a_j为退化速率参数，t为试验时间。

退化速率参数a_j与不同应力水平组合α之间满足如下多应力加速模型

$a_j={\eta }_{j0}\underset{i=1}{\overset{s}{\Pi }}{\left\{T_i\left(S_i\right)\right\}}^{{\eta }_{ji}}---\left(4\right)$

其中η_j0、η_ji为加速模型的系数，T_i(·)为任意的单调函数，S_i为第i种加速应力。常用的单应 力加速模型如下：当s＝1，T(S)＝exp(1/S)，式(4)即为Arrhenius加速模型a_j＝η_j0exp(η_j1/S)； 当s＝1，T(S)＝S，式(4)即为幂律模型a_j＝η_j0S^ηj1。式(4)两边取自然对数

即

$A_j\left(x\right)={\gamma }_{j0}^{\prime }+\underset{i=1}{\overset{s}{\Sigma }}{\gamma }_{ji}^{\prime }{x}_i---\left(6\right)$

式中，$A_j\left(x\right)=\ln a_j,{\gamma }_{j0}^{\prime }={\mathrm{ln\eta }}_{j0},{\gamma }_{ji}^{\prime }={\eta }_{ji},x_i=ln\left[T\left(S_i\right)\right],i=1,2,...,s.$x_i称为等效应 力水平。

将x_i其进行标准化，得到取值范围为[0,1]的标准化等效应力水平ξ_i

${\xi }_i=\frac{x_{iH}-x_i}{x_{iH}-x_{iL}}---\left(7\right)$

其中，x_iL和x_iH分别为应力x_i的最低水平和最高水平。因此，式(6)可重写为

其中，

$\left(\begin{array}{cc}{\gamma }_{j0}={\gamma }_{j0}^{\prime }+\underset{i=1}{\overset{s}{\Sigma }}{\gamma }_{ji}^{\prime }{x}_{iH},& {\gamma }_{ji}=-(x_{iH}-x_{iL}){\gamma }_{ji}^{\prime }\end{array}\right)$

1-4)产品性能退化的累积损伤模型信息

令ν_i表示第i个应力水平组合α_i下退化轨迹的起始时间，且此时的退化量与第i-1个应力 水平组合α_i-1结束时的退化量相等，则ν₁是以下方程的解

μ_j(ν₁|α₂)＝μ_j(t₁|α₁)(9)

类似地，ν_i满足

μ_j(ν_i|α_i+1)＝μ_j(ν_i-1+t_i-t_i-1|α_i)(10)

其中，i＝1,...,K-1。因此，μ_j(t)可表示为

${\mu }_j\left(t\right)=\left(\begin{array}{ccc}{\mu }_j\left(t\right|{\alpha }_1)& if& 0\le t<t_1\\ {\mu }_j(t-t_1+v_1|{\alpha }_2)& if& t_1\le t<t_2\\ ·& ·& \\ ·& ·& \\ ·& ·& \\ {\mu }_j(t-t_{K-1}+v_{K-1}|{\alpha }_K)& if& t_{K-1}\le t<t_K\end{array}\right)---\left(11\right)$

其中，j＝1,...,m。

因此，获取的产品加速退化试验模型参数先验信息可描述为

I＝(Σ,b_j,γ_ji,D_j),j＝1,…,m；i＝0,1,…,s(12)

步骤2、设计产品多应力多退化量步进加速退化试验基本方案。

Y_i的退化受到S₁,S₂,...,S_s种应力的影响，高于使用条件或贮存条件的s种应力组合能加 速Y_i退化过程。在进行多应力步进加速退化试验时，这s种加速应力的应力水平数均取为L。 s种加速应力的最高应力水平设置应不使加速退化试验过程中产品的退化机理发生改变，即产 品在这s种加速应力的加速退化试验中的退化机理与正常使用过程中的退化机理保持一致。

令表示一种应力水平组合。根据现有技术中的均匀设计及正交设计 原则(首先确定每种应力的应力水平数，然后选取相应的正交表，最后将应力及其水平按正 交表排列形成试验方案)，选取一系列应力水平组合α₁,α₂,…,α_K形成试验方案，其中K为应 力水平组合的个数。如果试验方案是分式析因设计方案，则K＝L^s-1。例如：当s＝3、L＝2时， 试验应力水平及其组合为α₁,α₂,α₃,α₄，如表1所示。

表1三种应力两水平的正交试验方案

注：“1”表示低应力水平，“2”表示高应力水平；此时，s＝3，L＝2，K＝4，$\left(\begin{array}{ccc}{\alpha }_2=(S_1^{\left(1\right)},S_2^{\left(2\right)},S_3^{\left(2\right)}),& {\alpha }_3=(S_1^{\left(2\right)},S_2^{\left(1\right)},S_3^{\left(2\right)}),& {\alpha }_4=(S_1^{\left(2\right)},S_2^{\left(2\right)},S_3^{\left(1\right)}),\end{array}\right)$这个步进加速退化试验方案有 四个应力水平组合。

在开展多应力多退化量步进加速退化试验时，随机抽取N个样品在应力水平组合α₁下进 行试验，每隔F单位时间测试一次性能参数(监测频率)，一共监测M₁次(监测次数)。当 试验进行到时间τ₁时，应力水平组合由α₁变为α₂，继续进行试验，监测频率为F，监测次数 为M₂。当试验进行到时间τ₂时，应力水平组合由α₂变为α₃，继续进行试验，监测频率为F， 监测次数为M₃。试验按如此方式进行，直到预定的时间结束。也就是说，应力水平组合最终 变为τ_K，监测频率为F，监测次数为M_K，试验到时间τ_K时试验全部结束。步进加速退化试 验每一应力水平组合下的试验时间为τ_i(i＝1,2,...,K)，且τ_i＝F·M_i·t_u，其中t_u为单位时间， 为1天或1小时。因此，总试验时间τ可表示为

$\begin{array}{c}\tau =\underset{i=1}{\overset{K}{\Sigma }}{\tau }_i\\ =\underset{i=1}{\overset{K}{\Sigma }}F·M_i·t_u\end{array}---\left(13\right)$

于是步进加速退化试验的应力水平组合随时间的变化规律可表示为

$\alpha =\left(\begin{array}{ccc}{\alpha }_1& if& 0\le t<t_1\\ {\alpha }_2& if& t_1\le t<t_2\\ ·& ·& \\ ·& ·& \\ ·& ·& \\ {\alpha }_K& if& t_{K-1}\le t<t_K\end{array}\right)---\left(14\right)$

其中，例如，当s＝3，L＝2，K＝4时，步进加速退化试验应力 水平组合剖面如图1及图2所示，其中图1为三应力水平的设置，图2为应力水平组合随时 间的变化规律。

步骤3、建立多应力多退化量步进加速退化试验方案优化模型。

3-1)确定优化模型的目标函数

将产品在使用应力水平组合下的p阶分位寿命估计的均方误差平方 根RMSE作为优化的目标函数：

$U=RMSE\left({\hat{\tau }}_{p0}\right)=\sqrt{E\left[{({\hat{\tau }}_{p0}-{\tau }_{p0})}^2\right]}---\left(15\right)$

其中，E[·]表示数学期望；τ_p0为p阶分位寿命；为p阶分位寿命估计。及τ_p0的 求解方法可描述如下。

p阶分位寿命τ_p0指产品在时刻τ_p0时的失效概率为p，此时产品的可靠度为1-p，即

R₀(τ_p0)＝1-p(16) 而产品在使用应力水平组合下时刻τ_p0的可靠度为

$\begin{array}{c}{R}_0\left({\tau }_{p0}\right)=P\{Y_{01}\left(t\right)\le D_1,Y_{02}\left(t\right)\le D_2,...,Y_{0m}\left(t\right)\le D_m\}\\ ={\int }_0^{D_m}...{\int }_0^{D_2}{\int }_0^{D_1}f(y_{01},y_{02},...,y_{0m}|{\tau }_{p0}){\mathrm{dy}}_{01}{\mathrm{dy}}_{02}...{\mathrm{dy}}_{0m}\end{array}---\left(17\right)$

通过联立方程(16)、(17)求解τ_p0。式(17)中，τ_p0为p阶分位寿命，R₀(τ_p0)为τ_p0时 刻产品的可靠度，Y₀₁(t),...,Y_0m(t)为产品在使用应力水平组合下时刻t的m 个退化量，D₁,...,D_m为产品m个退化量阈值，f(y₀₁,y₀₂,…,y_0m)为产品m个退化量的联合分 布密度函数。而f(y₀₁,y₀₂,…,y_0m)由步骤1中方程(12)确定，具体方法如下：

①将公式(12)中的γ_ji代入式(8)，可计算得到a_j；

②将a_j及公式(12)中的b_j代入式(3)，可计算得到μ_j；

③由μ_j及公式(12)中的Σ即可确定f(y₀₁,y₀₂,…,y_0m)的均值向量参数及方差协方差 矩阵参数，从而确定f(y₀₁,y₀₂,…,y_0m)。

若f(y₀₁,y₀₂,…,y_0m)的参数由试验数据或仿真数据估计出，则可用与上述步骤相同的方法 求解

3-2)确定优化模型的设计变量

多应力多退化量步进加速退化试验的每一要素都可以作为设计变量：

①试验应力S₁,S₂,…,S_s；

②试验应力水平

$\begin{array}{c}{S}_1^{\left(1\right)},S_1^{\left(2\right)},...,S_1^{\left(L\right)}\\ S_2^{\left(1\right)},S_2^{\left(2\right)},...,S_2^{\left(L\right)}\\ ·\\ ·\\ ·\\ S_s^{\left(1\right)},S_s^{\left(2\right)},...,S_s^{\left(L\right)}\end{array}---\left(18\right)$

及其组合α₁,…,α_K；

③试验样品个数N；

④监测频率F；

⑤应力水平组合α_j下的监测次数M_j。

因此，试验方案可表示为d＝(s,S_i,L,K,α_j,N,F,M_j)，i＝1,…,s，j＝1,…K。经过方 案的基本设计，多应力多变量步进加速退化试验的设计变量可简化为d＝(N,F,M_j)。

3-3)确定优化模型的约束条件

优化模型的约束条件如下：

①试验总费用C_T不超过试验预算C_b，C_T≤C_b；

②试验的应力数不低于2，s≥2；

③每一应力的水平数不少于2，L≥2；

④应力水平的组合数与应力水平数满足K＝L^s-1；

⑤试验样品个数不少于5，N≥5；

⑥监测频率不小于1个时间单位，F≥1；

⑦监测次数不少于3次，M_j≥3；

⑧

总的试验费用C_T由试验运行费用、测量费用、试验样品费用组成，由下式计算：

$C_T=C_{op}·F·\underset{j=1}{\overset{K}{\Sigma }}{M}_j+C_M·N·\underset{j=1}{\overset{K}{\Sigma }}{M}_j+C_d·N---\left(19\right)$

其中，C_op表示单位时间试验运行费用，C_M表示单次测量的费用，C_d表示样品单价。

综上所述，建立多应力多退化量步进加速退化试验方案优化模型可描述为

步骤4、按如图3所示流程对优化模型进行优化，具体如下：

4-1)根据约束条件构造可行试验方案集D，输入选取的试验方案数Z及蒙特卡罗仿真次 数N_mc，令z＝1；

4-2)在D中选取一个方案d_z＝(s,S_i,L,K,α_j,N,F,M_j)，i＝1,…,s，j＝1,…,K，z＝1,…,Z， 令n_mc＝1；

4-3)根据方案d_z及前述获取的先验参数信息I＝(Σ,b_j,γ_ji,β_j,D_j),j＝1,…,m；i＝0,1,…,s 计算仿真参数：

4-3-1)将α₁,α₂,…,α_K中的标准化等效应力水平代入加速模型(8)计算a_ji(j＝1,...,m， i＝1,...,K)；

4-3-2)根据式(3)、(9)、(10)、(11)可推导出

${\mu }_j\left(t\right)=\left(\begin{array}{ccc}{b}_j+a_{j1}t& if& 0\le t<t_1\\ b_j+a_{j2}(t-t_1)+a_{j1}{t}_1& if& t_1\le t<t_2\\ b_j+a_{j3}(t-t_2)+a_{j1}{t}_1+a_{j2}(t_2-t_1)& if& t_2\le t<t_3\\ ·& ·& \\ ·& ·& \\ ·& ·& \\ b_j+a_{jK}(t-t_{K-1})+{\Sigma }_{i=1}^{K-1}{a}_{ji}(t_i-t_{i-1})& if& t_{K-1}\le t<t_K\end{array}\right)---\left(21\right)$

其中j＝1,2,...,m。

4-3-3)令t＝F,2F,…,MF，其中根据式(21)计算 μ(t)＝[μ₁(t),μ₂(t),…,μ_m(t)]^T。

4-4)根据μ(t)及先验信息I中的Σ参数，对于t＝F,2F,…,MF生成N个m维正态分布 N(μ(t),Σ)向量如下

其中，Y_n(t_j)＝[Y_n1(t_j),Y_n2(t_j),…,Y_nm(t_j)]^T且n＝1,…,N，j＝1,…,M，t₁＝F,t₂＝2F,…,t_M…＝MF。 Y_n(t_j)从m维正态分布N(μ(t_j),Σ)中抽样获得。

4-5)分析仿真数据(22)，按如下步骤计算

4-5-1)对t＝F,2F,…,MF估计均值向量和方差协方差矩阵

$\hat{\mu }\left(t\right)={[{\mu }_1\left(t\right),{\mu }_2\left(t\right),...,{\mu }_m\left(t\right)]}^T---\left(23\right)$

其中

${\hat{\mu }}_j\left(t\right)=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{Y}_{ij}\left(t\right),j=1,...,m---\left(25\right)$

${\hat{\sigma }}_{kl}\left(t\right)=\frac{1}{N}\underset{q=1}{\overset{N}{\Sigma }}[Y_{qk}\left(t\right)-{\hat{\mu }}_k\left(t\right)][Y_{ql}\left(t\right)-{\hat{\mu }}_l\left(t\right)],k,l=1,...,m---\left(26\right)$

因此，方差协方差矩阵可用下式估计

$\hat{\Sigma }=\frac{1}{M}\underset{k=1}{\overset{M}{\Sigma }}{\hat{\Sigma }}_k---\left(27\right)$

其中可令t＝kF，由式(24)计算。

4-5-2)用模型(21)拟合4-5-1)得到的估计模 型(21)的参数${\hat{b}}_j,{\hat{\gamma }}_{ji},i=0,1,...,s,j=1,...,m.$

4-5-3)令即$\begin{array}{cccc}{S}_1=S_1^{\left(0\right)},& S_2=S_2^{\left(0\right)},& ...,& S_s=S_s^{\left(0\right)}\end{array}.$根据式(7)将 它们转化为标准应力水平并代入式(8)计算使用应力水平组合下的

4-5-4)将和代入式(3)确定使用应力水平组合下和t的关系，然后联立方程(16)、 (17)求解令n_mc＝n_mc+1。

4-5-5)如果n_mc≤N_mc则返回步骤4-3)，重复步骤4-3)～4-5-4)。否则，由以上步骤可得 根据先验信息I中的参数b_j,γ_ji,i＝0,1,…,s,j＝1,…,m，再由步骤4-5-3)、 4-5-4)求解τ_p0。由下式计算方案d_z对应的优化目标函数值U_z

$U_z=\sqrt{\frac{1}{N_{mc}}\underset{i=1}{\overset{N_{mc}}{\Sigma }}\left[{({\hat{\tau }}_{p0i}-{\tau }_{p0})}^2\right]}---\left(28\right)$

令z＝z+1。

4-6)如果z≤Z，则返回步骤4-2)选取另一方案，重复步骤4-3)～4-5)，否则转到步骤 4-7)。

4-7)选取使U(d)最小的试验方案作为最优试验方案d^*。

进一步地，步骤1-2)中，如果产品的退化量Y服从多维非正态分布，一般有两种方法 处理：一是将其变换为多维正态分布，二是将正态分布的情况扩展为非正态分布的情况。

进一步地，步骤1-3)中，可对试验时间t进行变换，使变换后的试验时间与μ_j成线性关 系，例如t＝t′^βj，其中t′为实际的日历时间，单位为小时或天，β_j则是待估计的系数。

进一步地，步骤4-5)中，仿真次数N_mc越多，蒙特卡罗仿真误差将越小，但较大的N_mc会增加计算量，可通过初步仿真分析确定。

本发明提供了一种多应力多退化量步进加速退化试验方案优化设计方法。该方法以蒙特 卡罗仿真优化理论为基础，将使用应力水平组合下产品p阶分位寿命误差的均方根作为优化 目标，将试验应力及其水平组合、试验样品个数、监测频率、每一应力水平组合下的监测次 数作为优化变量，将试验总费用及相关优化变量的整数限制作为约束条件，通过仿真、分析、 寻优得到试验方案的最优解。该方法突破了多应力多退化量加速退化试验方案优化目标函数 表达式难以求解的问题，易于流程化，便于工程应用，能够在满足试验费用约束条件的同时， 得到产品寿命指标的最精确估计，可为产品基于多应力多退化量加速退化试验的寿命预测提 供最优试验方案支持。

本发明提供的方法目前已经成功应用于某型安全阀橡胶密封圈的多应力多退化量加速退 化试验方案优化设计，为此型安全阀基于多应力多退化量加速退化试验的寿命预测提供最优 试验方案支持。

附图说明

图1三应力步进加速退化试验三应力水平的设置

图2三应力步进加速退化试验应力水平组合随时间的变化规律

图3多应力多退化量步进加速退化试验方案优化设计方法流程图

图4实施例中该型橡胶密封圈两应力两退化量步进加速退化试验设置

图5实施例中该型橡胶密封圈两应力两退化量步进加速退化试验应力水平组合随时间的 变化规律

图6实施例中得到的U(d)-d(RMSE-z)关系图

具体实施方式

下面以某安全阀橡胶密封圈为例进一步说明本发明所述方法的具体实施方式。该型橡胶 密封圈广泛应用于安全阀，而此型安全阀则是火箭发动机姿态控制的关键部件。正常贮存环 境的温度为25℃，温度为50％RH。工程经验表明，该型橡胶密封圈是此安全阀的寿命薄弱 环节。该型橡胶密封圈在贮存环境中会逐渐发生老化，直到其某些性能参数不能满足要求时 发生失效。该型橡胶密封圈老化决定了此安全阀的贮存寿命。为利用加速退化试验预测该型 橡胶密封圈的贮存寿命，需要设计其多应力多退化量步进加速退化试验的最优方案。

需要特别指出的是，以下实施仅用于说明目的，而非用于限定本发明的范围。

实施例1、某安全阀橡胶密封圈多应力多退化量步进加速退化试验优化设计

步骤1、获取该橡胶密封圈的加速退化试验相关信息。

1-1)产品的退化量及失效阈值信息

该橡胶密封圈在贮存过程中压缩永久变形量Y₁和压缩应力松弛系数Y₂随时间会逐渐变 大，当Y₁>0.3567或Y₁>0.3567时判定该橡胶密封圈贮存失效。即该橡胶密封圈有2个退化量 Y_i(i＝1,2)随时间逐渐退化，一旦某个退化量Y_i超过失效阈值D_i(i＝1,2)，D＝(0.3567,0.3567)^T， 该橡胶密封圈就会发生失效。高于正常贮存环境的温度S₁和湿度S₂能加速的退化速度。

1-2)该橡胶密封圈退化量的联合概率密度函数信息

时刻t该橡胶密封圈退化量Y＝(Y₁,Y₂)^T服从二维正态分布，其联合概率密度函数可表示 为

$f\left(Y\right)=\frac{1}{{\left(2\pi \right)}^{m/2}|\Sigma {|}^{1/2}}\exp (-\frac{1}{2}{\left(Y-\mu \right)}^T{\Sigma }^{-1}(Y-\mu \left)\right)---\left(29\right)$

其中，μ＝(μ₁,μ₂)^T为均值向量，Σ为方差协方差矩阵

$\Sigma =\left(\begin{array}{cc}{\sigma }_{11}& {\sigma }_{12}\\ {\sigma }_{12}& {\sigma }_{22}\end{array}\right)---\left(30\right)$

|Σ|为行列式值，|Σ|不随应力水平组合变化。通过预试验分析估计σ₁₁＝9.332×10^-5、 σ₁₂＝1.9064×10^-5、σ₂₂＝4.9297×10^-5。

1-3)该型橡胶密封圈退化模型及加速模型信息

在不同应力水平组合α下，样品二维正态分布均值向量μ的第j维元素与试验时间的关 系满足如下退化模型

$\left(\begin{array}{c}{\mu }_1=b_1+a_1{t}^{{\beta }_1}\\ {\mu }_2=b_2+a_2{t}^{{\beta }_2}\end{array}\right)---\left(31\right)$

式中，b_1、b₂为截距参数，a_1、a₂为退化速率参数，t为试验时间单位为天，β_1、β₂是时间变换系 数。b₁＝-0.08865、b₂＝-0.134793、β₁＝0.7、β₂＝0.4。

退化速率参数a_1、a₂与不同应力水平组合α之间满足如下加速模型

$\left(\begin{array}{c}{a}_1={\eta }_{10}·{\left\{T_1\left(S_1\right)\right\}}^{{\eta }_{11}}·{\left\{T_2\left(S_2\right)\right\}}^{{\eta }_{12}}\\ a_2={\eta }_{20}·{\left\{T\left(S_1\right)\right\}}^{{\eta }_{21}}·{\left\{T_2\left(S_2\right)\right\}}^{{\eta }_{22}}\end{array}\right)---\left(32\right)$

其中，由于加速应力分别为温度S₁(℃)、湿度S₂(％RH)，因此采用如下应力转换函数

$\left(\begin{array}{c}{T}_1\left(S_1\right)=\exp \left(\frac{1000}{S_1+273.16}\right)\\ T_2\left(S_2\right)=S_2\end{array}\right)---\left(33\right)$

式(32)两边取自然对数

$\left(\begin{array}{c}{\mathrm{lna}}_1={\mathrm{ln\eta }}_{10}+{\eta }_{11}ln\left\{T_1\left(S_1\right)\right\}+{\eta }_{12}ln\left\{T_2\left(S_2\right)\right\}\\ {\mathrm{lna}}_2={\mathrm{ln\eta }}_{20}+{\eta }_{21}ln\left\{T_1\left(S_1\right)\right\}+{\eta }_{22}ln\left\{T_2\left(S_2\right)\right\}\end{array}\right)---\left(34\right)$

即

$\left(\begin{array}{c}{A}_1\left(x\right)={\gamma }_{10}^{\prime }+{\gamma }_{11}^{\prime }{x}_1+{\gamma }_{12}^{\prime }{x}_2\\ A_2\left(x\right)={\gamma }_{20}^{\prime }+{\gamma }_{21}^{\prime }{x}_1+{\gamma }_{22}^{\prime }{x}_2\end{array}\right)---\left(35\right)$

式中，将x_i标准化，使加速 应力水平的取值范围为[0,1]

$\begin{array}{cc}{\xi }_1=\frac{x_{1H}-x_1}{x_{1H}-x_{1L}},& {\xi }_2=\frac{x_{2H}-x_2}{x_{2H}-x_{2L}}\end{array}---\left(36\right)$

其中，x_1L、x_2L和x_1H、x_2H分别为应力x₁、x₂的最低水平和最高水平

$\left(\begin{array}{c}{x}_{1L}=ln\left[T_1\left(S_{1L}\right)\right]=\frac{1000}{S_{1L}+273.16}=\frac{1000}{50+273.16}=3.0945\\ x_{1H}=ln\left[T_1\left(S_{1H}\right)\right]=\frac{1000}{S_{1H}+273.16}=\frac{1000}{115+273.16}=2.5763\\ x_{2L}=ln\left[T_2\left(S_{2L}\right)\right]=ln\left(S_{2L}\right)=ln\left(0.70\right)=-0.3567\\ x_{2H}=ln\left[T_2\left(S_{2H}\right)\right]=ln\left(S_{2H}\right)=ln\left(0.95\right)=-0.0513\end{array}\right)---\left(37\right)$

同时，正常贮存时环境应力水平

$\left(\begin{array}{c}{x}_1^{\left(0\right)}=\ln \left[T_1\left(S_1^{\left(0\right)}\right)\right]=\frac{1000}{S_1^{\left(0\right)}+273.16}=\frac{1000}{25+273.16}=3.3540\\ x_2^{\left(0\right)}=\ln \left[T_2\left(S_2^{\left(0\right)}\right)\right]=\ln \left(S_2^{\left(0\right)}\right)=\ln \left(0.50\right)=-0.6931\end{array}\right)---\left(38\right)$

因此，正常贮存应力、最低及最高加速应力如表2所示。

表2正常贮存应力、最低及最高加速应力

因此，式(35)可重写为

$\left(\begin{array}{c}{\mathrm{lna}}_1=A_1\left(x\right)={\gamma }_{10}+{\gamma }_{11}{\xi }_1+{\gamma }_{12}{\xi }_2\\ {\mathrm{lna}}_2=A_2\left(x\right)={\gamma }_{20}+{\gamma }_{21}{\xi }_1+{\gamma }_{22}{\xi }_2\end{array}\right)---\left(39\right)$

γ₁₀＝-2.9246、γ₁₁＝-0.93626、γ₁₂＝-0.47946、γ₂₀＝-2.2378、γ₂₁＝-0.52551、γ₂₂＝-0.31084。

1-4)该型橡胶密封圈性能退化的累积损伤模型信息

基于方程(3)及累积损伤理论，在应力水平组合α_i下，多应力多退化量步进加速退化试 验的μ_j(t)和相应恒定加速退化试验的μ_αij(t)，i＝1,...,K；j＝1,2如图4所示。

令ν_i表示组合α_i下退化轨迹的起始时间，且此时的退化量与组合α_i-1结束时的退化量相 等，则ν₁是以下方程的解

μ_j(ν₁|α₂)＝μ_j(t₁|α₁)

类似地，ν_i满足

μ_j(ν_i|α_i+1)＝μ_j(ν_i-1+t_i-t_i-1|α_i)

其中，i＝1,...,K-1。因此，μ_j(t)可表示为μ_αij(t)的形式

${\mu }_j\left(t\right)=\left(\begin{array}{ccc}{\mu }_j\left(t\right|{\alpha }_1)& if& 0\le t<t_1\\ {\mu }_j(t-t_1+v_1|{\alpha }_2)& if& t_1\le t<t_2\\ ·& ·& \\ ·& ·& \\ ·& ·& \\ {\mu }_j(t-t_{K-1}+v_{K-1}|{\alpha }_K)& if& t_{K-1}\le t<t_K\end{array}\right)$

其中，j＝1,2。

因此，获取的该型橡胶密封圈加速退化试验模型参数先验相关信息可描述为

I＝(Σ,b_j,γ_ji,D_j),j＝1,2；i＝0,1,2 ＝(σ₁₁,σ₁₂,σ₂₂,b₁,b₂,γ₁₀,γ₁₁,γ₁₂,γ₂₀,γ₂₁,γ₂₂) （40） ＝(9.332×10^-5,1.9064×10^-5,4.9297×10^-5,-0.08865,-0.134793,-2.9246,-0.93626, -0.47946,-2.2378,-0.52551,-0.31084)

步骤2、设计该型橡胶密封圈多应力多退化量步进加速退化试验基本方案。

如前所述，压缩永久变形量Y₁和压缩应力松弛系数Y₂的退化受到温度S₁和湿度S₂两种 应力的影响，高于贮存条件的两(s＝2)种应力的水平组合能加速Y₁和Y₂的退化过程。在进 行多应力步进加速退化试验时，这两种加速应力的应力水平数均取为L＝3。两种加速应力的 最高应力水平设置不使加速退化试验过程中产品的退化机理发生改变，即产品在这两种加速 应力的加速退化试验中的退化机理与正常使用过程中的退化机理保持一致。按此原则，设置 多应力多退化量步进应力加速退化试验的最高应力水平最 低应力水平$s_{Low}=(S_1^{\left(Low\right)},S_2^{\left(Low\right)})=(50,70\%),$中间应力水平$s_{Middle}=(S_1^{\left(Middle\right)},S_2^{\left(Middle\right)})$$=(80,82\%).$

令表示一种应力水平组合。根据均匀设计及正交设计原则，选取如表3所 示的应力水平组合α₁,α₂,α₃形成试验方案，此方案可更直观的表示为表4及图4。

表3两应力两退化量步进加速退化试验的均匀正交试验方案

表4两应力两退化量步进加速退化试验的均匀正交试验方案直观表示

随机抽取N个样品在应力水平组合α₁下进行试验，每隔F单位时间测试一次性能参数(监 测频率)，一共监测M₁次(监测次数)。当试验进行到时间τ₁时，应力水平组合由α₁变为α₂， 继续进行试验，监测频率为F，监测次数为M₂。当试验进行到时间τ₂时，应力水平组合由α₂变为α₃，继续进行试验，监测频率为F，一直试验到监测次数达到M₃时试验结束。步进加速 退化试验每一应力水平组合下的试验时间为τ_i(i＝1,2,3)，且τ_i＝F·M_i·t_u，其中t_u＝1天。因 此，总试验时间τ可表示为

τ＝τ₁+τ₂+τ₃＝F·M₁·t_u+F·M₂·t_u+F·M₃·t_u

于是步进加速退化试验的应力水平组合随时间的变化规律可表示为

$\alpha =\left(\begin{array}{ccc}{\alpha }_1& if& 0\le t<t_1\\ {\alpha }_2& if& t_1\le t<t_2\\ {\alpha }_3& if& t_2\le t<t_3\end{array}\right)$

其中，t₁＝τ₁,t₂＝τ₁+τ₂,t₃＝τ＝τ₁+τ₂+τ₃。图5为应力水平组合随时间变化规律示意。

步骤3、建立多应力多退化量步进加速退化试验方案优化模型。

3-1)确定优化模型的目标函数

将产品在使用应力水平组合下的p阶分位寿命估计的均方误差平方根

RMSE作为优化的目标函数：

$U=RMSE\left({\hat{\tau }}_{p0}\right)=\sqrt{E\left[{({\hat{\tau }}_{p0}-{\tau }_{p0})}^2\right]}---\left(41\right)$

其中，E[·]表示数学期望。及的求解方法可描述如下。

p阶分位寿命τ_p0指产品在时刻τ_p0时的失效概率为p，此时产品的可靠度为1-p，即

R₀(τ_p0)＝1-p(42) 而产品在使用应力水平组合下时刻τ_p0的可靠度为

$\begin{array}{c}{R}_0\left({\tau }_{p0}\right)=P\{Y_{01}\left(t\right)\le D_1,Y_{02}\left(t\right)\le D_2\}\\ ={\int }_0^{D_2}{\int }_0^{D_1}f(y_{01},y_{02},|{\tau }_{p0}){\mathrm{dy}}_{01}{\mathrm{dy}}_{02}\end{array}---\left(43\right)$

通过联立方程(35)、(36)求解τ_p0。其中D₁、D₂、f(y₀₁,y₀₂)的参数由步骤1中的相关先验 信息确定。若f(y₀₁,y₀₂)的参数由试验数据或仿真数据估计出，则可用类似的方法求解

3-2)确定优化模型的设计变量

由于经过了基本设计，因此该型橡胶密封圈多应力多退化量步进加速退化试验的设计变 量为：

①试验样品个数N；

②监测频率F；

③应力水平组合α₁、α₂、α₃下的监测次数M₁、M₂、M₃。

因此，试验方案可表示为d＝(N,F,M₁,M₂,M₃)。

3-3)确定优化模型的约束条件

优化模型的约束条件如下：

①试验总费用C_T不超过试验预算C_b＝60000，C_T≤C_b；

②试验样品个数不少于5，N≥5；

③监测频率不小于1个时间单位，F≥1；

④监测次数不少于3次，M₁,M₂,M₃≥3；

⑤

总的试验费用C_T由试验运行费用、测量费用、试验样品费用组成，由下式计算：

$C_T=C_{op}·F·\underset{j=1}{\overset{3}{\Sigma }}{M}_j+C_M·N·\underset{j=1}{\overset{3}{\Sigma }}{M}_j+C_d·N$

其中，C_op＝20元/小时表示单位时间试验运行费用，C_M＝5元/次表示单次测量的费用，C_d＝850 元/个表示样品单价。

令x＝(x₁,x₂,x₃,x₄,x₅)且x₁＝N、x₂＝F、x₃＝M₁、x₄＝M₂、x₅＝M₃。建立多应力多退

化量步进加速退化试验方案优化模型可描述为

步骤4、按以下流程对优化模型进行优化。

4-1)根据约束条件构造可行方案集D，输入选取的试验方案数Z＝2594，以及蒙特卡罗 仿真次数N_mc＝1000，令z＝1；

4-2)从D中选取一个方案d₁＝(N,F,M₁,M₂,M₃)＝(20,1,15,15,44)，令n_mc＝1；

4-3)根据方案d₁及前述获取的式(40)所示的先验参数信息I＝(Σ,b_j,γ_ji,β_j,D_j),j＝1,2, i＝0,1,2计算仿真参数：

4-3-1)将如表3所示的α₁,α₂,α₃标准化等效应力水平代入加速模型(32)计算a_ji(j＝1,2， i＝1,2,3)：

$\left(\begin{array}{c}{a}_{11}=\exp ({\gamma }_{10}+{\gamma }_{11}{\xi }_1^{\left(1\right)}+{\gamma }_{12}{\xi }_2^{\left(1\right)})=\exp (-2.9246-0.93626\times 0-0.47946\times 1.0)=0.0332\\ a_{12}=\exp ({\gamma }_{10}+{\gamma }_{11}{\xi }_1^{\left(2\right)}+{\gamma }_{12}{\xi }_2^{\left(2\right)})=\exp (-2.9246-0.93626\times 0.5-0.47946\times 0)=0.0336\\ a_{13}=\exp ({\gamma }_{10}+{\gamma }_{11}{\xi }_1^{\left(3\right)}+{\gamma }_{12}{\xi }_2^{\left(3\right)})=\exp (-2.9246-0.93626\times 1.0-0.47946\times 0.5)=0.0166\\ a_{21}=\exp ({\gamma }_{20}+{\gamma }_{\mathrm{21}}{\xi }_1^{\left(1\right)}+{\gamma }_{\mathrm{22}}{\xi }_2^{\left(1\right)})=\exp (-2.2378-0.52551\times 0-0.31084\times 1.0)=0.0782\\ a_{22}=\exp ({\gamma }_{20}+{\gamma }_{\mathrm{21}}{\xi }_1^{\left(2\right)}+{\gamma }_{\mathrm{22}}{\xi }_2^{\left(2\right)})=\exp (-2.2378-0.52551\times 0.5-0.31084\times 0)=0.0820\\ a_{23}=\exp ({\gamma }_{20}+{\gamma }_{\mathrm{21}}{\xi }_1^{\left(3\right)}+{\gamma }_{\mathrm{22}}{\xi }_2^{\left(3\right)})=\exp (-2.2378-0.52551\times 1.0-0.31084\times 0.5)=0.0540\end{array}\right)$

4-3-2)根据式(3)、(9)、(10)、(11)及b₁、b₂、a₁₁、a₁₂、a₁₃、a₂₁、a₂₂、a₂₃可推导出

${\mu }_1\left(t\right)=\left(\begin{array}{ccc}-0.08865+0.0332t^{0.7}& if& 0\le t<15\\ -0.08865+0.0336(t^{0.7}-6.6568)+0.0332\times 6.6568& if& 15\le t<30\\ -0.08865+0.0166(t^{0.7}-10.8140)+0.0332\times 6.6568+0.0336(10.8140-6.6568)& if& 30\le t<74\end{array}\right)$

${\mu }_2\left(t\right)=\left(\begin{array}{ccc}-0.134793+0.0782t^{0.4}& if& 0\le t<15\\ -0.134793+0.0820(t^{0.4}-2.9542)0.0782\times 2.9542& if& 15\le t<30\\ -0.134793+0.0540(t^{0.4}-3.8981)+0.0782\times 2\mathrm{.9542+0}\mathrm{.0820}(3.8981-2.9542)& if& 30\le t<74\end{array}\right)$

化简，得

${\mu }_1\left(t\right)=\left(\begin{array}{ccc}-0.08865+0.0332t^{0.7}& if& 0\le t<15\\ 0.1324+0.0336(t^{0.7}-6.6568)& if& 15\le t<30\\ 0.2720+0.0166(t^{0.7}-10.8140)& if& 30\le t<74\end{array}\right)---\left(45\right)$

${\mu }_2\left(t\right)=\left(\begin{array}{ccc}-0.134793+0.0782t^{0.4}& if& 0\le t<15\\ 0.0962+0.0820(t^{0.4}-2.9542)& if& 15\le t<30\\ 0.1736+0.0540(t^{0.4}-3.8981)& if& 30\le t<74\end{array}\right)---\left(46\right)$

4-3-3)令t＝F,2F,…,MF＝1,2,…,74，其中M＝M₁+M₂+M₃＝15+15+44＝74，将t 代入式(45)、(46)计算μ(t)＝[μ₁(t),μ₂(t)]^T。μ₁(t)、μ₂(t)均为1行74列的数组，μ(t)为2 行74列矩阵，此处省略。

4-4)根据μ(t)及式(40)所示先验信息I中的Σ参数，对于t＝1,2,…,74生成N个2维 正态分布N(μ(t),Σ)向量如下

其中，Y_n(t_j)＝[Y_n1(t_j),Y_n2(t_j)]^T且n＝1,…,20，j＝1,…,74，t₁＝1,t₂＝2,…,t_M＝74。Y_n(t_j)从2 维正态分布N(μ(t_j),Σ)中抽样获得。生成的仿真数据为40行74列矩阵，此处省略。

4-5)按如下步骤计算

4-5-1)对t＝1,2,…,74估计均值向量和方差协方差矩阵

$\hat{\mu }\left(t\right)={[{\mu }_1\left(t\right),{\mu }_2\left(t\right)]}^T$

$\hat{\Sigma }\left(t\right)=\left(\begin{array}{cc}{\hat{\sigma }}_{11}\left(t\right)& {\hat{\sigma }}_{12}\left(t\right)\\ {\hat{\sigma }}_{12}\left(t\right)& {\hat{\sigma }}_{22}\left(t\right)\end{array}\right)$

其中

${\hat{\mu }}_j\left(t\right)=\frac{1}{20}\underset{i=1}{\overset{20}{\Sigma }}{Y}_{ij}\left(t\right),j=1,2$

${\hat{\sigma }}_{kl}\left(t\right)=\frac{1}{20}\underset{q=1}{\overset{20}{\Sigma }}[Y_{qk}\left(t\right)-{\hat{\mu }}_k\left(t\right)][Y_{ql}\left(t\right)-{\hat{\mu }}_l\left(t\right)],k,l=1,2$

通过计算，得到的为2行74列矩阵，此处省略。方差协方差矩阵可用下式估计

$\hat{\Sigma }=\frac{1}{74}\underset{k=1}{\overset{74}{\Sigma }}{\hat{\Sigma }}_k$

其中，${\hat{\Sigma }}_k=\hat{\Sigma }\left(t_k\right)=\hat{\Sigma }\left(kF\right),$即

$\hat{\Sigma }=\left(\begin{array}{cc}8.882\times {10}^{-5}& 2.136\times {10}^{-5}\\ 2.136\times {10}^{-5}& 5.201\times {10}^{-5}\end{array}\right)$

4-5-2)用模型(21)拟合4-5-1)得到的及其中t＝1,2,…,74，估计模 型(21)的参数计算结果为：${\hat{b}}_1=-0.0884,{\hat{b}}_2=-0.1349,{\hat{\gamma }}_{10}-2.9234,$${\hat{\gamma }}_{11}=-0.9429,{\hat{\gamma }}_{12}=-0.4796,{\hat{\gamma }}_{20}=-2.2265,{\hat{\gamma }}_{21}=-0.5458,{\hat{\gamma }}_{20}=-\mathrm{0.3202.}$

4-5-3)令即S₁＝25，S₂＝50％。根据式(36)将它们转化为 标准应力水平并代入式(39)计算使用应力水平组合下的${\hat{a}}_j;$${\hat{a}}_1=0.0048,{\hat{a}}_2=\mathrm{0.0243.}$

4-5-4)将和代入式(31)确定使用应力水平组合下和t的关系，然后联立方程(42)、 (43)求解令n_mc＝n_mc+1。

4-5-5)如果n_mc≤N_mc则返回步骤4-3)，重复步骤4-3)～4-5-4)。否则，由以上步骤可得 (具体值省略)。根据先验信息I中的参数b_j,γ_ji,i＝0,1,2,j＝1,2，再由步骤 4-5-3)、4-5-4)求解τ_p0＝612.3302。由下式计算方案d₁对应的优化目标函数值U₁

$U_1=\sqrt{\frac{1}{1000}\underset{i=1}{\overset{1000}{\Sigma }}\left[{({\hat{\tau }}_{p0i}-{\tau }_{p0})}^2\right]}---\left(48\right)$

令z＝z+1。

4-6)如果z≤Z，则返回步骤4-2)选取另一方案，重复步骤4-3)～4-5)，否则转到步骤 4-7)。

4-7)得到使U(d)最小的最优试验方案d^*＝(21,1,15,18,39)，此时U(d^*)＝27。计算得到 的U(d)-d关系如图6所示，其中横坐标z为方案d的序号，纵坐标RMSE为方案d对应的 目标函数值U(d)，即估计误差的均方根值，d^*为最优试验方案。

在本发明的上述实例中，通过优化过程，得到在约束条件下的两应力两退化量步进加速 退化试验的最优试验方案。投入21个该型橡胶密封圈样品，以1天/次的监测频率，在 α₁＝(115℃,70％)下试验监测15次，然后将应力变化为α₂＝(80℃,95％)继续进行试验监测18 次，最后将应力变化为α₃＝(50℃,82％)继续进行试验监测39次，可得到正常贮存条件 α₀＝(25℃,50％)下95％可靠贮存寿命最精确估计，其误差均方根(RMSE)达到最小值27。