CN102622473A - 基于贝叶斯理论的步进应力加速退化试验优化设计方法 - Google Patents

Abstract

本发明提出一种基于贝叶斯理论的步进应力加速退化试验优化设计方法，应用在加速退化试验技术领域。首先，确定产品性能退化模型、加速模型，并基于历史数据给出模型参数的先验分布；其次，确定优化设计空间，构成试验方案集合；然后，建立期望效用函数或期望损失函数，确定优化目标，并基于马尔科夫蒙特卡洛方法，确定试验方案集合中的各设计的优化目标值；最后，利用曲面拟合方法，找到最优试验方案。本发明方法避免了传统(局部)试验优化设计方法在假定模型参数取值已知的情况下进行易出现较大偏差的缺点，在给出模型参数先验分布的情况下进行试验优化设计，得到的优化方案更合理更符合实际。

Description

基于贝叶斯理论的步进应力加速退化试验优化设计方法

技术领域

本发明是一种基于贝叶斯理论的步进应力加速退化试验优化设计方法，属于加速退化试验技术领域，用于解决可靠性与系统工程领域的技术问题。

背景技术

随着科技的迅速发展，高可靠、长寿命产品越来越多，例如：航空航天类电子产品、光电产品等，为了评估这些产品的可靠性与寿命，加速寿命试验(Accelerated Life Testing，ALT)和加速退化试验(Accelerated Degradation Testing，ADT)相继发展起来。ALT需要搜集产品的失效数据，但有时由于高可靠产品的特性，难以获得足够的失效数据以满足ALT评估时的需求。ADT是基于产品的物理特性或退化特性，通过在短时间内搜集加速应力下的退化增量来外推评估产品正常条件下的寿命和可靠性。ADT克服了ALT的不足，成为加速试验领域的研究热点。

如何在有限的资源下对ADT进行优化设计，科学合理的安排试验样本、试验时间、性能监测间隔等，以获得更有效的试验数据，进而提高产品可靠性评估精度是ADT的一项重要研究内容。ADT优化设计根据不同原则可分为不同范畴，包括：

1)应力施加方式。主要包括恒定和步进；

虽然对恒定应力加速退化试验(Constant Stress ADT，CSADT)的研究比对步进应力加速退化试验(Step Stress ADT，SSADT)相对成熟，但SSADT具有节省样本和试验时间的优势，因此，近年来，对SSADT的研究也逐步增多。

2)退化模型。主要包括混合效应模型、伽马过程和布朗运动(wiener过程)；

基于这三种模型的ADT优化设计都有学者进行研究。

3)加速应力类型的数量。主要包括单应力和多应力；

应力类型主要包括温度应力、湿度应力、电应力、振动应力等。单应力的ADT优化设计，采用单应力加速模型，例如以温度为加速应力，采用阿伦尼斯模型；多应力的ADT优化设计研究较少，采用多应力加速模型。

4)测量方式。主要包括破坏性测量和非破坏性测量；

测量产品性能参数时，对产品本身是否造成破坏，是ADT优化设计需要考虑的因素之一。非破坏性测量，如电子、光电类产品的性能测量等；破坏性测量，如胶类产品拉伸剪切强度的测量等，不同测量方式，测量次数，样本分配等不同。

5)历史信息。主要分为传统(局部)设计和贝叶斯设计。

对于大多数产品而言，通常都会存在大量的历史信息。比如，相似产品的历史试验数据或实际使用数据；产品自身在ADT之前进行的性能试验、环境试验等。传统(局部)的ADT优化设计是在给定退化模型基础上，通过假定模型参数取值来进行优化。当“假定”的取值与真实值偏差较大时，按照该假定值进行优化设计的试验方案就不能提供最有效的试验数据，往往会出现过试验或欠试验的情况，从而失去了试验“优化设计”的意义。如果基于贝叶斯理论综合利用历史信息进行优化设计，则比假定参数取值进行优化得到的方案可信度要高很多。但是ADT的Bayes优化设计还很少有研究人员涉及。

传统(局部)ADT优化设计已有大量研究：

对CSADT，Qishan Li基于混合效应模型，研究了CSADT的优化设计(参考文献[1]：Qishan Li.Accelerated Degradation Test Planning and Optimization[D].The University of Arizona，2002)。Yang G.，Yang K基于线性退化模型，将失效阈值作为一项加速应力，以最小化正常临界值和正常应力下平均(对数)寿命的渐进方差为目标来优化CSADT(参考文献[2]：Yang G.，Yang K.Accelerated degradation-tests with tightened critical values.IEEE Transactions onReliability，2002，51(4)：463-468)。H-F Yu等研究了退化率分别服从倒数Weibull分布和对数正态分布情况下加速退化试验的优化设计方法(参考文献[3]：Yu H.F.，Chiao C.-H.Designing adegradation experiment by optimizing the interval estimation of the percentile.IndustrialEngineering Research，2000，2：33-48；参考文献[4]：Yu H.F.Designing a degradation experimentby minimizing the variance of estimating a product′s mean-time-to-failure.Pan-PacificManagement Review，2002，5：171-181；参考文献[5]：Yu H.F.Optimal selection of the mostreliable product with degradation data.Engineering Optimization，2002，34：579-590；参考文献[6]：Yu H.F.，Chiao C.-H.An optimal designed degradation experiment fir reliability improvement.IEEE Transactions on Reliability，2002，51：427-433；参考文献[7]：Yu H.F.，Tseng S.T.Designing ascreening experiment for highly reliable products.Naval Research Logistics，2002，9：514-526；参考文献[8]：Yu H.F.Designing an accelerated degradation experiment with a reciprocal Weibulldegradation rate.Journal of statistical planning and inference，2006，136：282-297；参考文献[9]：Yu H.F.，Tseng S.T.Designing a degradation experiment with a reciprocal Weibull degradation rate.Quality Technology and Quantitative Management，2004，1(1)：47-63)。Ying Shi，Luis A Escobarand W Q Meeker研究了加速破坏退化试验的优化设计，以最小化指定失效时间分位点极大似然估计的渐进方差为优化目标，采用一般等价定理(general equivalence theorem，GET)来验证不同方案的最优性(参考文献[10]：Ying Shi，Luis AEscobar，William Q Meeker.AcceleratedDestructive Degrad ation Test Planning.Technometrics，2009，51(1)：1-13)。Yashun Wang等基于混合效应模型，以最小化产品寿命分布p分位估计的均方误差为目标，以试验费用为约束，采用Monte Carlo仿真，研究了ADT的优化设计(参考文献[11]：Yashun Wang，Chunhua Zhang，Xun Chen，Yongqiang Mo.Simulation-based optimal design for accelerated degradation tests.ICRMS.July 2009：1302-1306)。对SSADT，Park等采用对数正态分布和恒定退化率对步进应力加速破坏退化试验建立了退化模型，以产品正常条件下寿命分布p分位极大似然估计渐进方差最小为目标，研究了步进应力加速破坏退化试验的优化设计(参考文献[12]：Park S.J.，YumB.J.，Balamurali S.Optimal design of step-stress degradation tests in the case of destructivemeasurement.Quality Technology and Quantitative Management，2004，1(1)：105-124)。L.C.Tang等以试验费用为优化目标，以正常应力下平均寿命极大似然估计的渐进方差为约束，研究了SSADT设计方法(参考文献[13]：Tang L.，Yang G.，Xie M.Planning Of Step-stress AcceleratedDegradation Test Poceedings of The Annual Reliability and Maintainability Symposium，2004：287-292)。S-T Tseng等分别基于Weiner过程和Gamma过程研究了SSADT的优化(参考文献[14]：Tseng S.T.，Peng C.Y Reliability improvement via degradation experiments.Conference on Mathematical Methods in Reliability，Santa Fe，New Mexico，2004；参考文献[15]：Tseng  S.T.，NarayanaswamyBalakrishnan，Chih-Chun Tsai.Optimal  Step-Stress AcceleratedDegradation Test Plan for Gamma Degradation Processes.IEEE Transactions on Reliability，2009，58(4)：611-618)。Xiaoyang Li以最小化正常应力下竞争可靠性模型p分位寿命估计的渐进方差为目标，研究了考虑竞争失效机制下的SSADT的优化方法(参考文献[16]：Xiaoyang Li，Tongmin Jiang.Optimal Design for Step-Stress Accelerated Degradation Testing with CompetingFailure Modes Proceedings Annual Reliability and Maintainability Symposium，2009：64-68)。

贝叶斯试验设计是一种决策论，是统计理论一个重要分支。通常要在试验数据(后验分布)获得之前做出决策，称为预后验分析。预后验分析的目的是选择可最小化总成本的试验方案。总成本包括决策损失和试验费用。降低决策损失必然要增加试验费用，因此，需要平衡两部分成本。

在《贝叶斯可靠性》一书中，作者Michael S.Hamada，Alyson G.Wilson等基于混合效用模型，用正态分布和逆伽马分布作为混合效应模型参数的先验，以可靠性估计精度为目标，对试验样本量和性能检测次数进行了优化设计(参考文献[17]：Michael S.Hamada，Alyson G.Wilson，C.Shane Reese，Harry F.Martz.Bayesian Reliability.Springer Science+Business Media，LLC，New York，USA，2008：330-331)。Kathryn Chaloner and Isabella Verdinelli，Merlise A.Clyde分别给出了贝叶斯试验设计的综述(参考文献[18]：Chaloner，K，Verdinelli，I..BayesianExperimental Design：A Review.Statistical Science，1995，10：273-304；参考文献[19]：Merlise A.Clyde.Experimental Design：A Bayesian Perspective.Social and Behavioral Sciences，2001：1-22)。贝叶斯试验设计一般很难获得后验分布的数学表达式，因此衍生了解决该问题的两种方法：一种是基于仿真，一种是基于大样本理论。基于贝叶斯理论的ALT设计已有一些学者做出研究。Alaattin Erkanli等人从决策论的角度研究了不同效用函数的恒定ALT的贝叶斯设计(参考文献[20]：AlaattinErkanli，R.Soyer.Simulation-based designs for accelerated life tests.Journalof Statistical Planning and Inference，2000，90：335-348)。Gladys D.C.Barriga等人提出指数-威布尔寿命分布和Arrhenius模型下，基于马尔可夫链蒙特卡罗(MCMC)的ALT的贝叶斯方法。Yao Zhang and W.Q.Meeker利用大样本近似获得贝叶斯准则，研究了对数位置尺度分布和位置参数的线性加速模型的CSALT的贝叶斯设计，并提出用一般等价定理(GET)寻找最优方案，以解决非线性问题(参考文献[21]：Yao Zhang and William Q.Meeker.Bayesian Methods forPlanning Accelerated Life Tests.Technometrics，2006，48(1)：49-60)。Xiao Liu和L-C Tang等提出了序贯CSALT设计方法，根据该方法，首先实施高应力水平的试验以迅速获得失效数据，然后，基于贝叶斯推断方法，利用这些失效数据，建立低应力水平下的先验信息。通过最小化产品正常应力水平下寿命特征分位估计的渐进后验方差的预后验期望，使试验设计得以优化(参考文献[22]：Liu X，and Tang LC.A Sequential Constant-stress Accelerated Life TestingScheme and Its Bayesian Inference.Quality and Reliability Engineering International，2009，25(1)：91-109；参考文献[23]：Tang LC，and Liu X Planning Sequential Constant-StressAccelerated Life Tests With Step Wise Loaded Auxiliary Acceleration Factor.Journal of StatisticalPlanning and Inference.140(2010)1968-1985)。进一步，在他们2010的论文中，针对CSADT提出贝叶斯优化设计方法。以平方损失最小为优化准则，基于仿真并结合曲面拟合方法给出优化结果。通过对比，证明贝叶斯设计方法可明显提高ADT设计的健壮性，减少不确定性(参考文献[24]：Tang LC，and Liu X.Planning and Inference of a Sequential Accelerated Life Test.Journal of Quality Technology，2010，42(1)：103-118)。

但是，到目前为止还没有基于贝叶斯理论的SSADT优化设计方法的研究。

发明内容

本发明针对目前传统(局部)加速退化试验设计存在的当假定的参数取值与真实值偏差较大时，则传统(局部)优化设计得到的方案与真实的最优方案相差较远，失去了“优化”的意义的问题，提出了一种基于贝叶斯理论的步进应力加速退化试验优化设计方法，该方法基于贝叶斯理论，充分利用产品的历史数据、相似产品信息等，给出模型参数先验分布的情况下进行试验优化设计，得到的优化方案更可信更合理。

本发明提出的基于贝叶斯理论的步进应力加速退化试验优化设计方法，具体步骤为：

步骤一、确定产品性能退化模型、加速模型，进而基于历史数据给出模型参数的先验分布；

步骤二、确定优化设计空间，构成试验方案集合；

步骤三、建立期望效用函数或损失函数，明确优化目标；

步骤四、对每一试验方案，基于马尔科夫蒙特卡洛(Markov Chain Monte Carlo，MCMC)方法，利用WinBUGS14计算优化目标；

步骤五、利用曲面拟合方法，找到最优试验方案，确定最优决策变量。

所述的步骤二的具体过程是：在确定试验应力水平范围，结合试验费用约束确定总样本量和总测量次数的情况下，需要优化的决策变量为各试验应力水平和各应力水平下的测量次数分配，令S表示试验应力向量，S＝(s₁，...，s_j，..，s_k)，s_j表示第j个应力水平，令M表示测量次数向量，M＝(m₁，...，m_j，..，m_k)，m_j表示第j个应力水平s_j下的测量次数；然后，引入应力比ξj对向量S进行标准化：ξ_j＝(s_j-s_min)/(s_max-s_min)，其中，s_min、s_max分别表示试验施加的加速应力的上下限；引入测量次数比π_j对向量M进行标准化：

最后确定的设计空间D为：D＝S×M，将设计空间沿应力比和测量次数比的方向，在包括边界的取值范围内，进行等分取值，得到有限个设计，组成试验方案集合。

所述的步骤三中具体给出两种优化目标：

(1)基于相对熵的优化目标为：

$\underset{\mathrm{\η}\∈D}{\max }{E}_x{E}_{\mathrm{\θ}}\left[\log \frac{p\left(x\right|\mathrm{\θ})}{p\left(x\right)}\right]$

其中，x表示产品性能退化增量；θ表示模型参数向量；E_θ表示关于θ的期望；p(θ|x)表示模型参数后验分布概率密度函数；E_x表示对样本空间的信息取数学期望；从试验方案集合中的设计η中获得的信息I(η，x，p(θ))＝I₁(x)-I₀，I₁(x)表示后验分布中获得的总信息量，I₁(x)＝∫p(θ|x)log p(θ|x)dθ，I₀表示先验分布中包含的信息，I₀＝∫p(θ)log p(θ)dθ＝E_θlog p(θ)；p(θ)表示模型参数先验分布概率密度函数；p(x)是边际似然函数，p(x)＝∫p(x|θ)p(θ)dθ；p(x|θ)表示参数θ已知条件下的似然函数。

(2)基于期望平方损失的优化目标为：

$\underset{\mathrm{\η}\∈D}{\min }{E}_x{E}_{\mathrm{\θ}}\left[{\left|\right|\mathrm{\θ}-E\left[\mathrm{\θ}\right|x,\mathrm{\η}\left]\right||}^2\right]$

其中，||·||为欧几里得范数，E[θ|x，η]表示在给定试验方案η和得到相应的性能退化增量数据x下的参数后验期望。E_x和E_θ分别表示对样本空间和模型参数θ空间取数学期望。

所述的步骤五具体是：将试验方案集合中的各设计，及步骤四得到的对应的优化目标值，构成数对

η_r表示试验方案集合中的第r个设计，分别采用参数多项式回归和非参数局部加权回归散点平滑法对数据

进行曲面拟合。对基于期望平方损失的优化目标，找到最小的优化目标值对应的设计就是最优试验方案，对基于相对熵的优化目标，找到最大的优化目标值对应的设计就是最优试验方案。

本发明方法的优点和积极效果在于：

(1)本发明方法首次将贝叶斯理论引入到SSADT的优化设计中，基于贝叶斯理论，可充分利用试验前的历史数据，相似产品信息等，在确定模型参数的先验分布的情况下对试验进行优化设计，避免了传统(局部)试验优化设计方法在假定模型参数取值已知的情况下进行易出现较大偏差的缺点；

(2)本发明方法将贝叶斯试验设计的两个优化准则：期望平方损失函数和相对熵，引入到SSADT优化设计中，建立了SSADT贝叶斯优化设计框架，该优化设计框架的构建是ADT优化设计研究领域的一个突破点，为基于贝叶斯理论从其他角度研究ADT优化设计提供了思路和依据。这两个优化准则比传统方法的优化准则更符合实际；

(3)将期望平方损失优化准则首次运用到SSADT中，该准则可依据模型参数对优化设计的重要度来选择最优方案，得到的方案具有侧重点；

(4)将相对熵优化准则首次运用到加速退化试验中，该准则以期望信息增益最大为目标，全面考虑了试验中可获得的信息量，得到的优化结果更贴合工程实际。

附图说明

图1是本发明的基于贝叶斯理论的SSADT优化设计方法的流程图；

图2是本发明的优化设计方法的步骤二的设计空间D划分及方案η示意图；

图3是本发明实施例中以期望平方损失为优化目标的曲面拟合示意图：a是采用二次多项式回归法拟合，b是采用局部加权回归散点平滑法拟合；

图4是本发明实施例以相对熵为优化目标的曲面拟合示意图：a是采用二次多项式回归法拟合，b是采用局部加权回归散点平滑法拟合。

具体实施方式

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

本发明将贝叶斯试验设计的两个经典优化准则：期望平方损失函数和相对熵，引入到SSADT优化设计中，建立基于贝叶斯理论的SSADT优化设计框架，给出SSADT的贝叶斯优化设计方法的具体步骤。在实施例中选用漂移布朗运动为退化模型、Arrhenius模型为加速模型来阐述本发明提出的贝叶斯优化设计方法。

本发明提出的一种基于贝叶斯理论的步进应力加速退化试验优化设计方法，如图1所示，包括以下几个步骤：

步骤一、确定产品性能退化模型和加速模型，进而基于历史数据给出模型参数的先验分布。

产品性能退化模型应用较多的主要包括混合效应模型、伽马过程和布朗运动(Wiener过程)三种模型。常用的加速模型有阿伦尼斯(Arrhenius)模型、逆幂率模型、艾琳(Eyring)模型等，其形式都可表示对数线性形式：

其中，

是应力s的某一已知函数，例如，对阿伦尼斯(Arrhenius)模型而言，

s＝T，T为绝对温度；对逆幂率模型而言，

s可表示电压、电流、功率等；d(s)为性能退化率；a、b为常数。根据产品自身特点、敏感应力和性能参数退化情况等，确定产品性能退化模型和加速模型，进而确定退化模型的概率密度分布函数和对数似然函数等。

根据产品历史数据、相似产品信息，性能退化量分布情况，结合贝叶斯共轭先验分布理论确定产品性能退化模型和加速模型中未知参数的先验分布。

步骤二、确定优化设计空间，构成试验方案集合。

对SSADT，在确定试验应力水平范围，结合试验费用约束确定总样本量和总测量次数的情况下，需要优化的决策变量为各试验应力水平和各应力水平下的测量次数分配。

令S表示试验应力向量，S中包括k个元素，k为正整数，S＝(s₁，...，s_j，..，s_k)，每个元素表示一个应力水平，s_j表示第j个应力水平，例如，试验中有3个应力水平，则S包括3个元素，S＝(s₁，s₂，s₃)。令M表示测量次数向量，则M＝(m₁，...，m_j，..，m_k)，m_j表示第j个应力水平s_j下的测量次数(j＝1，...，k)。对两个向量S和M标准化，对S，引入应力比ξ_j，表示第j个应力水平s_j的应力比，ξ_j∈[0，1]，具体ξ_j由式(1)得到：

ξ_j＝(s_j-s_min)/(s_max-s_min)                  (1)

s_min、s_max分别表示试验可施加的加速应力的上下限，需要根据产品的具体情况，由工程人员结合经验来选择，一般s_max为产品不改变失效机理前提下的极限应力，s_min为高于正常应力水平s₀的某一值，即s₀＜s_min≤s₁＜s₂…＜s_k≤s_max，一般取s_k＝s_max。

对M，引入测量次数比π_j，表示第j个应力水平sj所对应的测量次数m_j占总测量次数的比例，π_j∈(0，1)可表示为：

${\mathrm{\π}}_j=m_j/{\mathrm{\Σ}}_{j=1}^k{m}_j,$ ${\mathrm{\Σ}}_{j=1}^k{\mathrm{\π}}_j=1---\left(2\right)$

设计空间D：D＝S×M，可由应力比和测量次数比划分，如图2所示，(ξ_j，π_j)为设计空间内的某一设计(或称方案)η，设计区间内划分的所有方案构成试验方案集合。

应力比和测量次数比都是连续变量，可以划分出无数个方案。由于目标函数需要大量抽样仿真计算，常常需要运行几天甚至是几周。为此本发明采用曲面拟合方案来避免这一问题。沿应力比和测量次数比方向，在其取值范围内等分取值，将设计空间包括边界划分为有限个方案组成试验方案集合，例如将应力比和测量次数比在各自的取值范围内等分为9份，也即在每个方向各取10个值，若将取值范围[0.1，1]等分为9份，则取了10个值：0.1，0.2，0.3，0.4，0.5，0.6，0.7，0.8，0.9，1，这样共100个试验方案。对试验方案集合中每一个方案计算目标函数，利用曲面拟合找出目标函数最大或最小的区域，进而找出最优方案。

本发明采用三维曲面拟合技术，两个自变量为应力比ξ₁和测量次数比π₁，因变量为目标函数。也即将优化目标函数表示为应力比ξ₁和测量次数比π₁的函数。其他试验应力水平下的应力比ξ_j(1＜j≤k)和测量次数比π_j(1＜j≤k)根据实际情况给出或表示为ξ₁、π₁的函数。

对两应力水平(也称为2点设计)，试验方案的决策变量为：应力水平s₁、s₂，检测次数m₁、m₂，或者应力水平比ξ₁、ξ₂，测量次数比π₁、π₂。一般取s₂＝s_max，因此ξ₂＝1；s₂下的测量次数比π₂＝1-π₁。

对应力水平数k≥3的情况，试验方案的决策变量为：应力水平s₁，s₂，...，s_k，检测次数m₁，m₂，...，m_k，或者应力水平比ξ₁，ξ₂，...，ξ_k，测量次数比π₁，π₂，...，π_k，。令s_k＝s_max，因此其应力比ξ_k＝1，令中间应力水平比ξ_j(1＜j＜k)表示为ξ₁和ξ_k的函数，例如对温度应力，可采用倒数等间隔；对电应力可采用对数等间隔。同样，对测量次数比，中间的π_j(1＜j＜k)可用π₁的函数表示，则

本发明实施例中为k＝2时的优化设计。

步骤三、建立期望效用函数或期望损失函数，明确优化目标。

本发明方法中通过比较试验方案集合中的每个设计η的期望效用函数值，将得到的最大期望效用函数值对应的设计，作为最优试验方案；或者通过比较试验方案集合中的每个设计η的期望损失函数值，将得到的最小期望损失函数值对应的设计，作为最优试验方案。本步骤主要是建立期望效用函数或期望损失函数，明确优化目标。

贝叶斯试验设计是贝叶斯决策问题，试验设计优化目标是最大化期望效用或最小化期望损失。

在设计空间D内，对每一设计η，计算其期望效用函数

找出其最大值对应的设计，即为最优决策，也就是最优试验方案η_最优：

期望效用函数

如下：

其中，U(η，x，θ)表示某给定方案η的效用函数；θ表示模型参数向量；x表示产品性能退化增量；p(θ|x，η)表示在某给定方案η以及得到相应的性能退化增量数据x下的参数后验分布的概率密度函数，p(x|η)表示在某给定方案η下的似然函数；Θ表示模型参数θ的取值空间。

损失函数L(η，x，θ)和效用函数U(η，x，θ)之间的关系为：

L(η，x，θ)＝-U(η，x，θ)                    (5)

下面介绍期望平方损失和相对熵两种不同优化目标。

(1)期望平方损失(Expection Quadratic Loss)

典型损失函数一般有平方损失、绝对误差损失、0-1损失等。在这些损失函数中，平方损失因为易于计算且和经典最小二乘的关系，成为最受欢迎的函数。平方损失函数是凸函数，因此相应的决策是唯一的。在贝叶斯框架下，只有当损失函数是平方损失时，贝叶斯估计才是后验均值。平方损失函数可表示为：

L(η，x，θ)＝||θ-E[θ|x，η]||²                    (6)

式(6)中，||·||为欧几里得范数，E[θ|x，η]表示在给定方案η以及得到相应的性能退化增量数据x下的参数后验期望。

对样本空间和模型参数空间取数学期望，得到方案η的期望平方损失函数：

L(η)＝E_xE_θ[||θ-E[θ|x，η]||²]                    (7)

式(7)中，E_x和E_θ分别表示对样本空间和模型参数θ空间取数学期望。

因此，基于期望平方损失的优化目标为：

(2)相对熵(KL Divergence)

很多学者选择用渐进近似作为优化准则，但这不是一个好的选择。因为渐进近似只在样本量足够大的情况下才合理，若以其作为优化准则，则历史信息贡献在参数估计中的作用会很小。且大样本在实际中很难满足。

在贝叶斯理论中，信息熵是一种先验分布和后验分布之间距离的量度。从香农(Shannon)信息角度讲，信息熵还代表了通过试验获得的信息量，也就是众所周知的相对熵。

基于这一事实，本发明从最大化相对熵角度，也即从信息增益角度，对SSADT进行优化设计。

根据Lindley的研究，先验分布中包含的信息I₀为：

I₀＝∫p(θ)log p(θ)dθ＝E_θlog p(θ)               (8)

其中，p(θ)表示模型参数先验分布概率密度函数；E_θ表示关于θ的期望。

从后验分布中获得的总信息量I₁(x)为：

I₁(x)＝∫p(θ|x)log p(θ|x)dθ                      (9)

其中，p(θ|x)表示模型参数后验分布概率密度函数。

Lindley在其研究中定义从试验方案η中获得的信息为：

I(η，x，p(θ))＝I₁(x)-I₀                           (10)

试验方案设计应在数据x获得之前给出，因此需要对样本空间的信息取数学期望，则I(η，x，p(θ))的期望为：

I(η，p(θ))＝Ex[I₁-I₀]                             (11)

其中，E_x表示对样本空间的信息取数学期望，I(η，p(θ))也称为期望信息增益(ExpectedInformation Gain，EIG)，可作为期望效用函数

基于贝叶斯理论，试验信息的期望可表示为：

$I(\mathrm{\η},p\left(\mathrm{\θ}\right))=E_x{E}_{\mathrm{\θ}}\left[\log \frac{p\left(x\right|\mathrm{\θ})}{p\left(x\right)}\right]---\left(12\right)$

式中，p(x)是边际似然函数，也是标准常量，如下所示：

p(x)＝∫p(x|θ)p(θ)dθ                (13)

其中，p(x|θ)表示参数θ已知条件下的似然函数。

因此，基于相对熵的优化目标为：

步骤四、对每一试验方案，基于MCMC方法，利用软件WinBUGS14计算优化目标。

多数情况下，很难得到后验的解析表达式，因此数值仿真计算是贝叶斯理论中常用的解决手段。同样，本发明中涉及的期望平方损失(7)和期望信息收益(EIG)(12)也很难写出解析表达式，因此下面将从仿真的角度分别介绍优化求解步骤。

(1)基于期望平方损失的优化目标计算，具体采用下面步骤：

子步骤1.1，从步骤二得到的试验设计空间(共R_D个方案)内取出一个设计(方案)η_r，r＝1，...，R_D。

子步骤1.2，对方案η_r，从其相应的先验分布中仿真抽取模型参数θ；下面采用θ_r表示设计η_r对应的仿真抽取的模型参数。

子步骤1.3，对每次仿真的模型参数θ_r，从抽样分布f(x|θ_r，η_r)中生成产品性能退化增量x_r，共R₁次；f(x|θ_r，η_r)表示在设计η_r和其相应的模型参数θ_r下的产品性能退化增量的概率密度函数。

子步骤1.4，对生成的产品性能退化增量x_r，基于MCMC利用WinBUGS14计算对应的期望平方损失(即参数后验方差)，则方案η_r共有R₁个损失L_r，取均值即得到方案η_r的期望损失函数的值；

子步骤1.5，回到子步骤11，对设计空间中每一个方案重复子步骤1.2～子步骤1.4；得到所有试验方案集合内所有方案的期望损失函数的值。

(2)基于相对熵的优化目标计算，具体采用下面步骤：

式(11)所示的期望信息增益(EIG)也可写作：

$E_x{E}_{\mathrm{\θ}}\left[\log \frac{p\left(x\right|\mathrm{\θ})}{p\left(x\right)}\right]=E_x{E}_{\mathrm{\θ}}\left[\log p\left(x\right|\mathrm{\θ})\right]-E_x\left[\log p\left(x\right)\right]---\left(14\right)$

由于上式通常很难具有显示表达式，因此一般采用蒙特卡罗(Monte Carlo)仿真方法来计算。首先，式(13)中的p(x|θ)是似然函数，可以在参数空间和样本空间直接采用Monte Carlo仿真方法计算E_xE_θ[logp(x|θ)]，计算公式如下：

$E_x{E}_{\mathrm{\θ}}\left[\log p\left(x\right|\mathrm{\θ})\right]=\frac{1}{R_2}\·{\mathrm{\Σ}}_{h=1}^{R_2}\log p\left(x_h\right|{\mathrm{\θ}}_h)---\left(15\right)$

式(14)中R₂是仿真次数，通常取较大的整数，一般可选为100。

其次，对于边际似然函数p(x)而言，本发明方法采用Laplace-Metropolis算法来估计，计算公式如下：

$p\left(x\right)\≈{\left(2\mathrm{\π}\right)}^{d/2}{\left|{\mathrm{\Σ}}_{\mathrm{\θ}}\right|}^{1/2}p\left(x\right|\stackrel{\‾}{\mathrm{\θ}})p\left(\stackrel{\‾}{\mathrm{\θ}}\right)---\left(16\right)$

$\stackrel{\‾}{\mathrm{\θ}}=\frac{1}{R_{\mathrm{ML}}}\·\underset{g=1}{\overset{R_{\mathrm{ML}}}{\mathrm{\Σ}}}{\mathrm{\θ}}_g\mathrm{and}{\mathrm{\Σ}}_{\mathrm{\θ}}=\frac{1}{R_{\mathrm{ML}}-1}\·\underset{g=1}{\overset{R_{\mathrm{ML}}}{\mathrm{\Σ}}}({\mathrm{\θ}}_g-\stackrel{\‾}{\mathrm{\θ}}){({\mathrm{\θ}}_g-\stackrel{\‾}{\mathrm{\θ}})}^T---\left(17\right)$

其中，π是圆周率，π≈3.14；d为模型参数向量的维数；R_ML是仿真的性能退化增量x的个数，为样本量n与总检测次数m的乘积，即R_ML＝n×m；θ_g是基于MCMC采用软件WinBUGS14计算得到的第g个性能退化增量的参数后验均值，

是由N个仿真退化增量数据得到的θ_g的均值；∑_θ是参数后验的方差协方差矩阵。

基于上述公式(13)～(16)，基于相对熵的优化目标的求解步骤如下：

子步骤2.1，从试验设计空间(共R_D个方案)内取方案η_r，r＝1，...，R_D。

子步骤2.2，对方案η_r，从其相应的先验分布中仿真抽取R₂次参数θ_rh，h表示仿真次数，本步骤中h＝1，...，R₂，并利用每次仿真的参数θ_rh，从抽样分布f(x|θ_rh，η_r)中生成退化产品性能增量数据x_rh。

子步骤2.3，根据仿真生成的退化增量数据x_rh，计算方案η_r的对数似然函数logp(x_rh|θ_rh，η_r)，并根据式(14)计算方案η_r的E_xE_θ[logp(x_rh|θ_rh，η_r)]：

$E_x{E}_{\mathrm{\θ}}\left[\log p\left(x_{\mathrm{rh}}\right|{\mathrm{\θ}}_{\mathrm{rh}},{\mathrm{\η}}_r)\right]=\frac{1}{R_2}\·{\mathrm{\Σ}}_{h=1}^{R_2}\log p\left(x_{\mathrm{rh}}\right|{\mathrm{\θ}}_{\mathrm{rh}},{\mathrm{\η}}_r);$

子步骤2.4，对方案η_r，从其相应的先验分布中仿真抽取参数θ_r，针对θ_r仿真R₃次退化增量数据x_rh，h表示仿真次数，本步骤中h＝1，...，R₃；基于每次仿真得到的退化增量数据，结合MCMC方法和WinBUGS14软件计算第g个性能退化增量数据的参数后验均值θ_g，进而根据式(16)得到θ_g的均值

和参数后验的方差协方差矩阵∑_θ：

$\stackrel{\‾}{\mathrm{\θ}}=\frac{1}{R_{\mathrm{ML}}}\·\underset{g=1}{\overset{R_{\mathrm{ML}}}{\mathrm{\Σ}}}{\mathrm{\θ}}_g\mathrm{and}{\mathrm{\Σ}}_{\mathrm{\θ}}=\frac{1}{R_{\mathrm{ML}}-1}\·\underset{g=1}{\overset{R_{\mathrm{ML}}}{\mathrm{\Σ}}}({\mathrm{\θ}}_g-\stackrel{\‾}{\mathrm{\θ}}){({\mathrm{\θ}}_g-\stackrel{\‾}{\mathrm{\θ}})}^T;$

子步骤2.5，根据式(15)计算边际似然函数p(x_rh)：

$p\left(x_{\mathrm{rh}}\right)\≈{\left(2\mathrm{\π}\right)}^{d/2}{\left|{\mathrm{\Σ}}_{\mathrm{\θ}}\right|}^{1/2}p\left(x_{\mathrm{rh}}\right|\stackrel{\‾}{\mathrm{\θ}})p\left(\stackrel{\‾}{\mathrm{\θ}}\right);$

进而得到 $E_x\left[\log p\left(x\right)\right]=\frac{1}{R_3}\·{\mathrm{\Σ}}_{h=1}^{R_3}\log p\left(x_{\mathrm{rh}}\right);$ 其中，表示模型参数后验均值

已知条件下的似然函数，表示模型参数后验均值

带入先验分布得到的概率值。

子步骤2.6，根据式(13)计算EIG，即得到方案η_r的期望效用函数

子步骤2.7，回到子步骤2.1，对设计空间中每一个方案重复子步骤2.2～2.6；得到所有方案集内所有方案的优化目标。

步骤五、利用曲面拟合方法，比较所有方案的目标值，找到最优方案，确定最优决策变量。在有些情况下，曲面回归模型的确切结构很难确定。因此本发明选用非参数方法对所有方案的目标值进行回归拟合，例如，核平滑法。为比较不同模型的回归拟合效果，本发明方法分别用参数多项式回归和非参数局部加权回归散点平滑法(Locally Weighted ScatterplotSmoothing，LOWESS)对数据进行曲面拟合。

基于步骤四得到的所有方案目标值构成的数对

进行拟合曲面，对期望平方损失准则，找到目标值最小对应的方案即为最优试验方案；对相对熵准则，找到目标值最大对应的方案即为最优试验方案。

在得到最优试验方案后，也就得到了最优的应力比ξ_j和测量次数比π_j，进而确定出具体的加速应力水平及其相应的测量次数。

实施例：

基于贝叶斯理论对某超辐射发光二极管(super luminescent diode，SLD)进行步进应力加速退化试验(SSADT)优化设计，步骤如下：

步骤一、确定产品性能退化模型、加速模型，进而基于历史数据给出模型参数的先验分布。

首先给出模型及假设，如下，

(1)假设：

A1：退化趋势是单调不可逆的；

A2：退化失效机理不随应力而改变；

A3：在正常应力水平s₀和k个加速应力水平s₁＜s₂＜...＜s_k下，性能退化过程Y_l服从漂移布朗运动，漂移系数d(s_l)＞0，扩散系数σ_l＞0，l＝1，...，k：

Y_l(t)＝σ_lB(t)+d(s_l)·t+y₀                  (18)

A4：扩散系数σ_l不随应力量值改变，即，σ₀＝σ₁＝...＝σ_k＝σ。

A5：漂移系数d(s_l)也可称为退化率，是应力s的函数，也即漂移系数是加速模型

其中，s_l是第l个加速应力水平，

是应力s_l某种形式的函数，例如当加速应力是绝对温度时

A6：参数a，b和1/σ²的先验分布相互独立。

(2)模型：

以温度为加速应力，加速模型选为Arrhenius模型：

d(s_l)＝exp[a+b/s_l]                    (20)

1)退化模型为漂移布朗运动：

Y_l(t)＝σ_lB(t)+d(s_l)·t+y₀            (21)

其中，令初始值y₀＝0。则模型未知参数向量

2)似然函数：

设n个产品实施k水平的SSADT。假设试验中没有因为性能退化引起的失效。SSADT中，第l个应力水平下性能监测次数为m_l，则SSADT的累积监测次数为产品性能监测时间间隔为Δt，则第l个应力水平下的试验时间t_l＝m_l·Δt，总试验时间为t＝m·Δt。第l个应力水平下第i个受试产品的第j测量时间为t_ilj(i＝1，...，n，l＝1，...，k，j＝1，...m_l)，监控到的性能值为y_ilj。布朗运动是高斯过程，因此检测时间间隔Δt上的性能退化增量x独立且服从均值为d(s)Δt，方差为σ²Δt的正态分布，即x～N(d(s)Δt，σ²Δt)。独立增量x的概率密度函数为，

$f\left(x\right|\mathrm{\θ})=\frac{1}{\mathrm{\σ}\sqrt{2\mathrm{\π\Δt}}}\exp \{-\frac{{[x-d\left(s\right)\·\mathrm{\Δt}]}^2}{{2\mathrm{\σ}}^2\mathrm{\Δt}}\}---\left(22\right)$

根据累积损伤假设和式(21)，k个应力水平SSADT下n个样本所有退化增量的似然函数如下：

$p\left(x\right|\mathrm{\θ})=\underset{i=1}{\overset{n}{\mathrm{\Π}}}\underset{l=1}{\overset{k}{\mathrm{\Π}}}\underset{j=1}{\overset{m_l}{\mathrm{\Π}}}\frac{1}{\sqrt{{2\mathrm{\π\σ}}^2\mathrm{\Δt}}}\exp \{-\frac{{[x_{\mathrm{ilj}}-\exp (a+b/s_l)\·\mathrm{\Δt}]}^2}{{2\mathrm{\σ}}^2\mathrm{\Δt}}\}---\left(23\right)$

其对数似然函数为，

$\log p\left(x\right|\mathrm{\θ})=-\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{l=1}{\overset{k}{\mathrm{\Σ}}}\underset{j=1}{\overset{m_l}{\mathrm{\Σ}}}\left\{\right[\ln \left(2\mathrm{\π\Δt}\right)+\ln \left({\mathrm{\σ}}^2\right)]+\frac{{[x_{\mathrm{ilj}}-\exp (a+b/s_l)\·\mathrm{\Δt}]}^2}{{\mathrm{\σ}}^2\mathrm{\Δt}}\}---\left(24\right)$

然后，基于历史数据给出模型参数的先验分布：

根据性能退化增量x～N(d(s)Δt，σ²Δt)以及共轭先验分布理论，可假设参数a和b分别独立服从均值为μ_a，μ_b，方差为σ_a，σ_b为的正态分布。σ²的倒数服从尺度参数为α，形状参数为β的伽马分布。也即

$a~N({\mathrm{\μ}}_a,{\mathrm{\σ}}_a^2)---\left(25\right)$

$b~N({\mathrm{\μ}}_b,{\mathrm{\σ}}_b^2)---\left(26\right)$

$\frac{1}{{\mathrm{\σ}}^2}~\mathrm{\Γ}(\mathrm{\α},\mathrm{\β})---\left(27\right)$

根据历史数据等信息，确定模型参数的先验分布为，

$p\left(a\right)~N({\mathrm{\μ}}_a,{\mathrm{\σ}}_a^2)=N\left(\mathrm{16,1}\right)---\left(28\right)$

$p\left(b\right)~N({\mathrm{\μ}}_b,{\mathrm{\σ}}_b^2)=N(-\mathrm{7383,64})---\left(29\right)$

$p\left(\frac{1}{{\mathrm{\σ}}^2}\right)~\mathrm{\Γ}(\mathrm{\α},\mathrm{\β})=\mathrm{\Γ}\left(\mathrm{1,100}\right)---\left(30\right)$

步骤二、设置设计空间，构成试验方案集合。

本实例SSADT中温度应力水平数为2，即k＝2。令样本量n＝3。设计空间划分如下：

在2-水平设计中，加速应力向量S＝(s₁，s₂)。令最低温度水平为s_min＝40℃，最高为s_max＝80℃，则应力比ξ₁＝(s₁-40)/40，ξ₂＝1。需要优化的变量为ξ₁，ξ₁∈[0，1]，在其取值范围内等分9份，取10个值。

测量次数向量M＝(m₁，m₂)，令性能退化的总测量数m为100。测量次数分配比π₁＝m₁/m，π₂＝1-π₁。需要优化的变量为π₁，π₁∈(0，1)，为保证每一应力下采集到足够多的性能退化量，其测量次数不宜过小，因此令π₁∈[0.1，0.9]在其取值范围内等分9份，取10个值。因此，设计空间内共划分100个方案构成方案集合。

步骤三、建立期望效用函数或损失函数，明确优化目标。

按照具体实施方式中步骤二的过程，分别基于期望平方损失和相对熵两种优化准则，明确优化目标。

基于期望平方损失：

min||θ-E[θ|x，η]||²                    (31)

基于相对熵：

$\max E_x{E}_{\mathrm{\θ}}\left[\log \frac{p\left(x\right|\mathrm{\θ})}{p\left(x\right)}\right]---\left(32\right)$

步骤四、对每一试验方案，基于MCMC方法，利用WinBUGS14计算优化目标。

对每一试验方案，基于MCMC方法，利用WinBUGS14软件，按照具体实施方式中步骤四中两种优化准则的具体计算步骤计算优化目标。

步骤五、利用曲面拟合方法，比较所有方案的优化目标值，找到最优试验方案，确定最优决策变量。

按照具体实施方式中步骤五，对得到的所有方案的优化目标值构成数对采用多项式(二次)回归和局部加权(二次)回归散点平滑法(LOWESS)进行拟合曲面。

在曲面拟合前，需要对数据进行预处理。由于Monte Carlo仿真时容易从分布的尾端仿真出野点。另外，分布的均值和尾端之间的偏差会以平方的方式放大。这些异常值会导致决策(优化设计选择)偏差。因此，需要将这些异常值剔除掉。3σ(σ是数据的标准差)准则是一种很受欢迎的数据过滤方法，本发明方法采用3σ准则来剔除数据集的野点。对预处理后的数据进行曲面拟合，结果如下：

(1)基于期望平方损失的优化结果：

基于期望平方损失的优化结果见图3所示，z轴代表期望平方损失，用风险risk表示；图3中x轴代表应力水平s₁的应力比ξ₁，用xi₁表示；y轴代表应力水平s₁下的测量次数比π₁，用pi₁表示；用灰度表示z轴期望平方损失值大小，灰度越大(即颜色越重)的区域说明该区域内期望平方损失值越小，则最优方案在该区域内。从图3中可以看出，用二次多项式和LOWESS拟合，三个参数参数a、b和σ的优化结果都不一致，很难找出合适的最优解。一般根据实际情况位于等高线图四个角落(ξ₁＝0 or 1 and π_i＝0.1 or 0.9)的解不是最优结果。因此，应在图的中间区域寻找最优设计。从这点看用LOWESS拟合要优于用二次多项式拟合。

不同参数的优化结果不同，一般考虑以对试验影响最大的参数的优化结果作为最优解，由于参数b代表激活能E_a(b＝-E_a /k，k is the Boltzmann’s constant，k＝8.6171×10^-5eV/K)，对优化结果的影响较大，而参数a和σ的影响较小(参考文献[25]：Li，XY；Jiang，TM；Sun，FQ；Ma，J.Constant stress ADT for superluminescent diode and parameter sensitivity analysis.EKSPLOATACJA I

-Maintenance and Reliability.2010，46(2)：21-26)，因此，以b的结果为主，同时参考a和σ的结果来选择优化方案。在参数b的等高线图的中间区域，最小值大约在(ξ₁，π₁)＝(0.3，0.75)和(ξ₁，π₁)＝(0.75，0.55)附近，考虑到合理外推，加速应力的范围应加大，因此选择(ξ₁，π₁)＝(0.3，0.75)为优化方案。

将应力比和测量次数比转换为对应的应力和测量次数，结果见表1。

表1基于期望平方损失的优化结果

(2)基于相对熵的优化结果：

图4中，z轴代表相对熵，用KL divergence表示；x轴代表应力水平s₁的应力比ξ₁，用xi₁表示；y轴代表应力水平s₁下的测量次数比π₁，用pi₁表示；用灰度表示z轴相对熵的大小，灰度越小(即颜色越浅)的区域说明该区域内相对熵越大，则最优方案在该区域内。

从图4可以看出，设计空间内用LOWESS拟合的结果中相对熵最大的区域在(ξ₁，π₁)＝(0.05，0.45)附近，将应力比和测量次数比转换为应力和测量次数，结果如表2所示。

表1基于相对熵的优化结果

Claims (4)

1.一种基于贝叶斯理论的步进应力加速退化试验优化设计方法，其特征在于，具体包括如下步骤：

步骤一、确定产品性能退化模型、加速模型，进而基于历史数据给出模型参数的先验分布；

步骤二、确定优化设计空间，构成试验方案集合，具体过程是：

在确定试验应力水平范围，结合试验费用约束确定总样本量和总测量次数的情况下，需要优化的决策变量为各试验应力水平和各应力水平下的测量次数分配，令S表示试验应力向量，S＝(s₁，...，s_j，...，s_k)，s_j表示第j个应力水平，令M表示测量次数向量，M＝(m₁，...，m_j，...，m_k)，m_j表示第j个应力水平s_j下的测量次数；

然后，引入应力比ξ_j对向量S进行标准化：ξ_j＝(s_j-s_min)/(s_max-s_min)，其中，s_min、s_max分别表示试验施加的加速应力的上下限；

引入测量次数比π_j对向量M进行标准化：

最后确定的设计空间D为：D＝S×M，将设计空间沿应力比和测量次数比的方向，在包括边界的取值范围内，等分取值，得到有限个设计，组成试验方案集合；

步骤三、建立期望效用函数或期望损失函数，确定优化目标；

具体本步骤给出两种优化目标：

(1)基于相对熵的优化目标为：

$\underset{\mathrm{\η}\∈D}{\max }{E}_x{E}_{\mathrm{\θ}}\left[\log \frac{p\left(x\right|\mathrm{\θ})}{p\left(x\right)}\right]$

其中，x表示产品性能退化增量；θ表示模型参数向量；p(θ|x)表示模型参数后验分布概率密度函数；p(x)是边际似然函数，p(x)＝∫p(x|θ)p(θ)dθ；p(x|θ)表示模型参数θ已知条件下的似然函数，p(θ)表示模型参数先验分布概率密度函数；E_θ表示关于θ的期望；E_x表示对样本空间的信息取数学期望，从设计η中获得的信息I(η，x，p(θ))＝I₁x)-I₀，I₁(x)表示后验分布中获得的总信息量，I₁(x)＝∫p(θ|x)log p(θ|x)dθ，I₀表示先验分布中包含的信息，I₀＝∫p(θ)log p(θ)dθ＝E_θlog p(θ)；

(2)基于期望平方损失的优化目标：

$\underset{\mathrm{\η}\∈D}{\min }{E}_x{E}_{\mathrm{\θ}}\left[{\left|\right|\mathrm{\θ}-E\left[\mathrm{\θ}\right|x,\mathrm{\η}\left]\right||}^2\right]$

其中，||·||为欧几里得范数，E[θ|x，η]表示在给定试验方案η和得到相应的性能退化增量数据x下的参数后验期望，E_x和E_θ分别表示对样本空间和模型参数θ空间取数学期望；

步骤四、基于马尔科夫蒙特卡洛方法，利用软件WinBUGS14确定试验方案集合中的各设计的优化目标值；

步骤五、利用曲面拟合方法，找到最优试验方案，具体是：将试验方案集合中的各设计，及步骤四得到的对应的优化目标值，构成数对

η_r表示试验方案集合中的第r个设计，分别采用参数多项式回归和非参数局部加权回归散点平滑法对数据

进行曲面拟合，对基于期望平方损失的优化目标，找到最小的优化目标值对应的设计就是最优试验方案，对基于相对熵的优化目标，找到最大的优化目标值对应的设计就是最优试验方案。

2.根据权利要求1所述的一种基于贝叶斯理论的步进应力加速退化试验优化设计方法，其特征在于，步骤一中所述的退化模型为混合效应模型、伽马过程模型或者布朗运动随机过程模型，所述的加速模型为阿伦尼斯模型、逆幂率模型或者艾琳模型。

3.根据权利要求1所述的一种基于贝叶斯理论的步进应力加速退化试验优化设计方法，其特征在于，步骤四中所述的利用软件WinBUGS14确定试验方案集合中的各设计的优化目标值，当以期望平方损失最小为优化目标时，确定优化目标值的步骤为：

子步骤11：从试验方案集合中取出一个设计η_r，r＝1，...，R_D，R_D表示试验方案集合中设计的总个数，η_r表示第r个设计；

子步骤1.2：对方案η_r，从相应的先验分布中仿真抽取模型参数，标记该模型参数为θ_r；

子步骤1.3：对模型参数θ_r，从抽样分布f(x|θ_r，η_r)中生成产品性能退化增量x_r，共R₁次；f(x|θ_r，η_r)表示在设计η_r和f(x|θ_r，η_r)表示在设计η_r和其相应的模型参数θ_r下的产品性能退化增量的概率密度函数；

子步骤1.4：对生成的退化增量，基于马尔科夫蒙特卡洛方法利用软件WinBUGS14计算对应的期望平方损失，则设计η_r共有R₁个损失L_r，对R₁个损失L_r取均值，得到设计η_r的期望损失函数值；

子步骤1.5：回到子步骤1.1，对试验方案集合中每一个设计重复子步骤1.2～子步骤1.4，得到所有试验方案集合内所有设计的期望损失函数值。

4.根据权利要求1所述的一种基于贝叶斯理论的步进应力加速退化试验优化设计方法，其特征在于，步骤四中所述的利用软件WinBUGS14确定试验方案集合中的各设计的优化目标值，当以相对熵最大为优化目标时，确定优化目标值的步骤为：

子步骤2.1：从试验方案集合内取一个设计η_r，r＝1，...，R_D，R_D表示试验方案集合中设计的总个数，η_r表示第r个设计；

子步骤2.2：对设计η_r，从其相应的先验分布中仿真抽取R₂次模型参数θ_rh，h＝1，...，R₂，并利用第h(h＝1，...，R₂)次抽取的模型参数θ_rh，从抽样分布f(x|θ_rh，η_r)中生成产品性能退化增量x_rh；f(x|θ_rh，η_r)表示在设计η_r和f(x|θ_r，η_r)表示在设计η_r和其相应的模型参数θ_r下的产品性能退化增量的概率密度函数；

子步骤2.3：根据仿真生成的产品性能退化增量x_rh，得到方案η_r的对数似然函数log p(x_rh|θ_rh，η_r)，并进一步得到方案η_r的E_xE_θ[log p(x_rh|θ_rh，η_r)]：

$E_x{E}_{\mathrm{\θ}}\left[\log p\left(x_{\mathrm{rh}}\right|{\mathrm{\θ}}_{\mathrm{rh}},{\mathrm{\η}}_r)\right]=\frac{1}{R_2}\·{\mathrm{\Σ}}_{h=1}^{R_2}\log p\left(x_{\mathrm{rh}}\right|{\mathrm{\θ}}_{\mathrm{rh}},{\mathrm{\η}}_r)$

子步骤2.4：对设计η_r，从相应的先验分布中仿真抽取模型参数θ_r，针对模型参数θ_r仿真R₃次产品性能退化增量x_rh，h＝1，...，R₃，基于每次仿真得到的产品性能退化增量，基于马尔科夫蒙特卡洛方法利用软件WinBUGS14得到第g个产品性能退化增量的模型参数的后验均值θ_g，进而得到R₃个θ_g的均值

和模型参数的后验的方差协方差矩阵∑_θ：

$\stackrel{\‾}{\mathrm{\θ}}=\frac{1}{R_{\mathrm{ML}}}\·\underset{g=1}{\overset{R_{\mathrm{ML}}}{\mathrm{\Σ}}}{\mathrm{\θ}}_g\mathrm{and}{\mathrm{\Σ}}_{\mathrm{\θ}}=\frac{1}{R_{\mathrm{ML}}-1}\·\underset{g=1}{\overset{R_{\mathrm{ML}}}{\mathrm{\Σ}}}({\mathrm{\θ}}_g-\stackrel{\‾}{\mathrm{\θ}}){({\mathrm{\θ}}_g-\stackrel{\‾}{\mathrm{\θ}})}^T$

其中，R_ML是仿真的产品性能退化增量x_rh的个数，为样本量与总检测次数的乘积；

子步骤2.5：根据下式确定边际似然函数p(x_rh)：

$p\left(x_{\mathrm{rh}}\right)\≈{\left(2\mathrm{\π}\right)}^{d/2}{\left|{\mathrm{\Σ}}_{\mathrm{\θ}}\right|}^{1/2}p\left(x_{\mathrm{rh}}\right|\stackrel{\‾}{\mathrm{\θ}})p\left(\stackrel{\‾}{\mathrm{\θ}}\right)$

进而得到 $E_x\left[\log p\left(x\right)\right]=\frac{1}{R_3}\·{\mathrm{\Σ}}_{h=1}^{R_3}\log p\left(x_{\mathrm{rh}}\right);$ 其中，

表示模型参数后验的均值

已知条件下的似然函数，d为模型参数向量的维数，

表示模型参数后验均值

带入先验分布得到的概率值；

子步骤2.6：对于设计η_r，根据子步骤2.3与子步骤2.4得到的数据，确定设计η_r的期望效用函数

的值：

子步骤2.7：回到子步骤2.1，对试验方案集合中每一个设计重复步骤2.2～2.6，得到所有设计的期望效用函数的值。

Images