自适应SA-PSO优化的威布尔混合分布参数估计方法及应用

郭 森，王大为，张绍伟，姚永超

(上海机电工程研究所，上海 201109)

0 引言

随着装备复杂化程度不断提升，综合保障成为影响其作战使用的重要因素，而通过装备全寿命周期数据的应用分析挖掘出故障及使用维修规律，确定预防性维修间隔，是提升其综合保障能力的关键环节[1-2]。其中，根据装备失效维护数据确定其寿命分布是重要途径。威布尔分布是一种连续型概率分布模型，它能够描述各种类型机械零部件失效数据的分布规律，在寿命数据分析、可靠评估和预测等方面得到了一定程度的应用[3]。对于复杂装备而言，失效数据往往不会对各失效模式进行区分，用简单威布尔分布并不能有效地对设备的寿命进行描述[4]。多重威布尔混合模型具有良好的拟合性能，更适合于复杂装备失效规律的描述，但也会引入更多的未知参数，增加了模型参数估计算法的复杂度[5]。

常见的威布尔分布参数估计方法包括图解法[6]、期望极大化(EM)法、极大似然估计(MLE)法[7]和贝叶斯法[8]等。其中：图解法是粗略的估计，精度较低；EM算法和MLE算法需求解超越方程，在大样本下估计精度较高，但是求解过程比较复杂，有时不易收敛，且面对截尾小样本时偏差较大；贝叶斯法可以充分利用专家经验、设计参数等信息提高参数估计的准确性，并且随着使用过程中新信息的出现不断更新结果，是解决小样本参数估计的有效方法，但对设备先验信息要求较高。

为了取得更精确的寿命分布模型，并减少对先验信息的依赖性，许多学者引入了非线性拟合优化方法来对威布尔模型参数进行估计。文献[9]将两重二参数威布尔混合分布参数估计转化为非线性最小二乘优化问题，并采用Quasi-Newton法进行迭代求解，能够得到更加精确的结果，但理论要求比较高；文献[10]采用灰色模型与支持向量机相结合的方式，解决了小样本情况下三参数威布尔分布参数估计精度较低的问题；文献[11]采用混合粒子群(PSO)算法对模型似然函数进行求解，避免了传统极大似然估计方法的短板，一定程度上解决了威布尔混合分布模型的参数估计问题；文献[12]引入模拟退火(SA)算法，直接对威布尔分布模型进行参数估计，取得了较好的结果。

针对威布尔混合分布模型复杂、参数估计难度大的问题，本文提出采用基于自适应改进模拟退火粒子群(SA-PSO)算法加以解决。将威布尔混合分布模型参数估计视为一个非线性优化问题，采用粒子群算法求解，并对粒子速度与位置更新进行自适应改进，在全局最优值的选择上引入模拟退火机制，进一步改善粒子群算法的全局搜索能力。将改进的算法应用到某型柴油机喷油器的失效数据分析中，利用图解法对参数初值进行粗估，用于自适应SA-PSO优化算法的种群初始化，对威布尔混合分布进行参数估计，结果表明该算法能够有效提高参数估计的精度和效率。

1 模拟退火粒子群算法及其改进

1.1 PSO算法基本原理

PSO算法是一种用于非线性函数优化的集群智能优化算法[13]，它在解决大规模非线性问题中具有更高的优化效率和更好的优化结果。该算法将各个备选解称为“粒子”，全部的粒子集合称作“种群”。PSO算法首先生成初始种群，即在可行解空间中随机初始化生成一群粒子，并根据要优化的目标函数确定与之相对应的适应度值。优化过程就是初始种群在可行解空间内部飞行，并寻找最大或最小适应度值的过程，粒子飞行的方向与距离由速度决定。

其数学过程描述如下所述。

vid(t+1)=vid(t)+c1r1(pid(t)-xid(t))+

c2r2(gd(t)-xid(t))

(1)

xid(t+1)=xid(t)+vid(t+1)

(2)

1≤i≤N,1≤d≤D

c1、c2为正数，称作加速因子,c1调节粒子飞向自身最好位置方向上的步长，c2调节粒子飞向全局最好位置方向上的步长;r1、r2为[0,1]的随机数。

通过比较更新前后各粒子的适应度值，对单个粒子最优位置和全局最优位置进行更新，直到满足迭代终止条件。

1.2 PSO算法的自适应改进

与其他群体智能算法相比，PSO算法概念简单，易于编程实现求解，但在优化过程中同样会面临易于陷入局部极值点、搜索精度不高、收敛速度慢的问题。为此，引入惯性权重因子ω，对粒子飞行速度进行调节，提高PSO算法的收敛速度及全局搜索能力，使其能够快速得到全局最优结果。惯性权重因子ω对速度的调节原理为

vid(t+1)=ωvid(t)+c1r1(pid(t)-xid(t))+

c2r2(gd(t)-xid(t))

(3)

由式(3)可以看出，引入的ω代表的含义主要是上一代粒子状态对当前粒子状态的影响程度，ω取值越大，则上一代粒子状态对与当前粒子的状态影响越大；反之，影响则越小。同时，算法的全局搜索能力也会随着ω的增大而提高。因此，通过控制惯性权重因子ω的大小，可以改进粒子群算法的全局搜索能力，当ω值取为常数1时，式(3)即退化为基本粒子群算法的速度更新公式。

为了进一步提升优化性能，并加快粒子群算法的收敛速度，引入距离控制因子C(t)对惯性权重因子及加速系数进行自适应调整[14]，原理为

(4)

D为各粒子到全局最优值的平均距离；max(D)为至今平均距离的最大值；C(t)为距离控制因子，即各粒子到全局最优点的平均距离与平均距离最大值的比值。

1.3 基于模拟退火的自适应粒子群优化策略

模拟退火的思想来自常用于金属材料提升性能的热处理工艺退火，通过模拟高温下金属逐渐降温直至热平衡的过程，引入概率接收准则，并适当地控制温度参数的下降，反复迭代寻找能量最小值点，最后得到问题近似的全局最优解[12]。

将模拟退火思想引入PSO算法的优化过程，改善算法在训练过程精度降低容易发散的问题。在最优粒子的取舍中加入Metropolis准则，在接收更优粒子的同时以一定的概率接收较差粒子，能够从局部极值中跳出，随着温度的降低逐渐收敛到全局最优位置。

基于模拟退火的自适应粒子群优化过程如下：

a.初始化种群中各粒子的位置和速度。

b.计算各粒子的适应度值，将当前各粒子位置和适应度值存储在各粒子的pi中，将所有Pb中最优个体位置和适应度值存储在pb中。

c.根据式(5)确定初始温度t0，即

t0=f(pb)/ln5

(5)

d.根据式(6)确定当前温度下各pi的适配值，即

(6)

f.计算各粒子新的适应度值，更新各粒子的pi值及群体的pb值。

g.采用衰减系数法进行退温操作，即

tk+1=λtk

(7)

h.若满足终止条件搜索停止，输出结果，否则返回步骤d继续搜索。

2 威布尔分布及其参数估计方法

2.1 二参数威布尔分布

由于威布尔分布既能描述失效率上升的产品寿命，又能描述失效率下降的产品寿命，因此威布尔分布是可靠性工程中适于描述寿命分布规律的统计模型之一。

若设备寿命t服从二参数威布尔分布，则其累积失效概率密度函数、可靠度和失效率为：

(8)

(9)

(10)

β、θ分别为威布尔模型的形状参数、尺度参数。

2.2 威布尔混合分布模型

复杂机械系统包含多个零部件，每个零部件的故障可能是在多种失效机理共同作用下发生的，而且在不同的寿命阶段，不同的失效机理对系统的失效起主导作用。上述情况下，采集到的可靠性试验数据在威布尔概率图(WWP)上表现为曲线，如果用标准的单威布尔分布模型来描述时，会出现较大的误差，威布尔混合分布提供了更符合实际的寿命数据拟合模型[3]。

如果1个总体由n个子体组成，设各子体的概率密度函数分别为f1(t)，f2(t),…,fn(t)，各子体的混合权重分别为p1(t),p2(t),…,pn(t)，权重之和为1，则总体概率密度函数f(t)为

(11)

若各子体都服从威布尔分布，即

(12)

βi为第i个威布尔分布的形状参数；θi为第i个威布尔分布的尺度参数；pi为第i个威布尔分布的权重，它表示第i种失效原因或失效模式的失效数据在产品总的失效数据中所占的比例。

2.3 基于自适应SA-PSO算法的参数估计优化

对于单一威布尔分布而言，可以经过简单的数学变换将样本数据的非线性关系转换线性关系，并采取线性回归方法对数据进行拟合。但威布尔混合分布模型较为复杂，一般的数学变换无法将模型向线性问题进行转换，因此采用非线性优化方式对其参数估计问题进行求解[15]。针对两重威布尔混合分布模型，基于误差平方和最小的思想，建立多参数威布尔混合分布的非线性最小二乘(NLS)优化模型。

将样本容量为m的失效数据按寿命大小进行顺序排列(t1,t2,…,tm)，采用中位秩公式计算累计失效概率的观测值[3]，如式(13)：

(13)

(14)

根据误差平方和最小思想构建威布尔混合模型参数非线性最小二乘优化模型，如式(15)：

(15)

采用自适应改进的模拟退火粒子群算法对威布尔混合分布参数优化模型进行求解，流程如图1所示。在威布尔分布的参数估计方法中，图解法参数估计值的精度较低，往往不能满足实际使用需求，但是可以作为其他方法的迭代初值。因此，先对设备失效数据用图解法进行求解，得到迭代初始值，之后按照基于模拟退火的自适应粒子群优化策略进行迭代优化，直到满足终止条件。

图1 基于自适应SA-PSO的威布尔混合分布参数估计流程

3 实例分析

为了验证本文所提方法的准确性，采用文献[3]中某型柴油机喷油器服从两重威布尔混合分布的故障里程进行拟合优化，其数据如表1所示。根据式(13)构造样本集，其可靠度函数为

(16)

不可靠度为

F(t)=1-R(t)

(17)

优化目标函数为

(18)

图解法粗估值为p=0.620 5，θ1=9 996.6，β1=1.408 5，θ2=2 250.3，β2=11.007 9。在此结果基础上构造自适应SA-PSO优化算法的初始种群，以优化目标函数式(18)作为适应度函数按图1所示流程进行优化。同时，为了验证自适应SA-PSO方法的优势，将其结果与标准PSO算法、自适应PSO算法结果进行对比。

表1 某型柴油机喷油器的故障里程

3.1 自适应SA-PSO参数优化性能分析

不同优化算法对模型解决参数估计问题的优化过程如图2所示。从图2可以看出，标准PSO算法在迭代704步之后适应度值收敛至0.031 70，自适应PSO算法在迭代至250步时适应度值收敛至0.030 60，说明二者在收敛精度上相差不大，但后者收敛速度明显提高；与前2种算法相比，自适应SA-PSO算法在迭代过程中，适应度值于第113步便快速达到局部极值点0.033 25，并随后在第527步、第780步产生突变，最终收敛到0.021 00，这表明，引入模拟退火机制并经权重自适应更新的粒子群算法能够有效跳出局部极值点收敛到更低的适应度值，具有更优的全局搜索能力以及更快的收敛速度。综合分析，3种算法中，自适应SA-PSO算法优化性能最优，自适应PSO算法次之，标准PSO算法由于未经过任何改进，收敛速度较慢且易于陷入局部最优，优化性能最弱。

图2 不同算法参数优化过程

3.2 威布尔分布参数估计结果分析

采用图解法、基于Levenberg-Marquardt (L-M)的非线性最小二乘法(NLS)、标准PSO算法、自适应PSO算法和自适应SA-PSO算法，得到的模型参数估计结果及相对均方根误差如表2所示。采用基于L-M的NLS方法得到的模型相对均方根误差为0.089 1，通过智能优化算法拟合的模型误差都要更小，而且随着算法的改进，误差也变得越来越小，其中经过自适应改进的SA-PSO算法结果最好，拟合模型的相对均方根误差仅为0.024 8，与优化性能分析结果相符。

表2 不同方法参数估计结果

根据各方法参数估计结果画出相应的模型曲线，如图3所示。图中能直观看出喷油器失效数据中位秩点比较均匀的分布于相应拟合曲线两边。

结果表明，与传统方法相比，基于群体智能优化算法拟合的模型具有更好的精度，其中自适应SA-PSO算法拟合精度最好。根据该算法结果得到某型柴油机喷油器的可靠度函数为

(19)

图3 不同方法拟合曲线结果对比

4 结束语

本文对威布尔混合分布模型传统参数估计方法进行了分析。针对模型复杂、求解难度较大的问题，在残差最小化思想基础上建立了多重两参数威布尔混合分布的非线性最小二乘优化模型，并提出了一种基于自适应SA-PSO算法加以求解，最后采用某型柴油机喷油器失效数据，对本文所构造的参数估计模型及求解方法进行了验证，可以得到如下结论：

a.与标准PSO算法及经过自适应改进的PSO算法相比，通过引入模拟退火机制对优化过程全局最优值选取进行控制，可以有效避免陷入局部最优，同时可以提高算法的收敛速度。

b.与传统用于解决参数估计的图解法、非线性最小二乘法相比，采用基于自适应SA-PSO算法求解威布尔混合分布模型参数的最小二乘估计值，能够最大程度上减少先验信息缺乏的影响，有效简化了计算的复杂度，并且具有更高的求解精度。