基于ITLS和DE的加工中心三参数威布尔分布*

王晓峰 张英芝 申桂香 龙哲 张立敏

(1.吉林大学 建设工程学院，吉林 长春130022;2.吉林大学 机械科学与工程学院，吉林 长春130022)

近些年来，威布尔分布模型已经广泛应用于加工中心可靠性研究. 在进行加工中心的可靠性建模时，常假设数据服从位置参数γ=0 的二参数威布尔分布［1-2］.但加工中心并非总在t=0 时刻发生故障，并且二参数威布尔分布的处理过程中人为地丢掉了一个参数，增大了其参数估计结果的误差. 因此，选用位置参数γ≠0 的三参数威布尔分布更有利于反映加工中心可靠性的真实情况. 三参数威布尔分布的参数估计比较复杂，有关学者进行了许多相关研究［3-5］.在不断探讨其在加工中心可靠性建模应用的过程中发现以下问题［6］:①对加工中心三参数威布尔分布模型的参数优化多采用迭代法，即通过一定步长不断增大优化参数值，计算每次迭代得到的目标函数值，最后取最优目标函数所对应的参数值为估计值.这种迭代法如果步长过大可能会丢失最优解，步长过小则以牺牲时间为代价，有时迭代法可能会耗时很长.②在对三参数威布尔分布进行参数估计时多采用最小二乘法确定回归直线方程. 由于选择x 方向或者y 方向拟合所得的回归直线是不同的［7］，所以当在y 方向上拟合直线时，只考虑到y 坐标上的误差而没有考虑到x 坐标上的误差. 也就是说只考虑了故障发生概率F 的误差，没有考虑间隔时间t 的误差. 但是故障间隔时间t 的误差必然存在，即x 坐标误差必然存在，所以x 坐标的误差也应该予以充分考虑.③在三参数威布尔线性变换后，并非是尺度参数α 直接使得目标函数达到最佳，而且目标函数所考虑的误差是y 坐标上的误差，并没有直接考虑到故障发生概率F 的误差，但实际上两组误差并不相同.

针对上述问题，文中采用一种新兴的优化算法(即差分进化算法)进行参数优化，它可以准确快速地寻得目标函数全局最优解，克服了迭代法的缺点.其次，结合三参数威布尔分布特点，文中提出一种改进的整体最小二乘法. 该方法一方面具有整体最小二乘法的固有特点，即同时兼顾x 坐标和y 坐标误差，弥补了最小二乘法的不足，使所得结果更加真实合理;另一方面以故障发生概率F 的偏差平方和作为目标函数进行优化，确保最终所得的三参数威布尔曲线与原始数据点最大程度地接近.

1 试验数据

对23 台加工中心进行了现场试验，根据所得数据整理出加工中心整机的105个故障间隔时间数据，具体见表1.

表1 加工中心故障间隔时间1)Table 1 Time between failures of machining center h

2 基于改进的整体最小二乘法的参数估计

三参数威布尔分布函数为

式中:α 为尺度参数，α ＞0;β 为形状参数，β ＞0;γ 为位置参数，γ≥0.

对式(1)两次取对数后线性变换为

设一元线性回归模型为

最小二乘(LS)法的约束准则为

整体最小二乘(TLS)法的约束准则为

由此可见:①最小二乘法的实质是求一条直线，使所有已知点到这条直线的纵坐标的距离平方和最小，因此，LS 法只顾及了y 坐标的误差没有涉及x坐标的误差.②整体最小二乘法的实质是求一条直线，使所有的已知点到此直线的距离和的平方最小，所以TLS 法同时顾及了x 坐标和y 坐标的误差. 因此将整体最小二乘法代替最小二乘法引入加工中心可靠性建模后，既可以考虑到故障发生概率F 的误差，又可以考虑到故障间隔时间t 的误差，使所得三参数威布尔模型更加真实合理.

2.1 对β 的参数估计

整体最小二乘法最早由Golub 等［8］在1980年提出.Byungsoo［9］总结了TLS 法的重要理论结果和计算方法，发现在一些典型应用中用TLS 法相较于用LS 法可以使参数估计精度显著提高.本节主要通过TLS 法对威布尔分布中的形状参数β 进行估计.

根据式(5)取QTLS关于a 的偏导并令其为0:

由上式得到a 的表达式如下:

把式(6)代入式(5)中，得到

根据式(7)求QTLS关于b 的偏导并令其为0:

根据上式可得

因为经三参数威布尔分布线性变换后yi=ln ln所以

又因为b=β，所以此时β 的估计值可以表示成仅含有未知参数γ 的表达式β(γ).

2.2 对α 和γ 的参数估计

根据三参数威布尔线性变换过程中的特点，可以发现:首先，式(2)中a= -βlnα，b=β，其中β 是可以直接求的，但是求出的α 并非直接使QTLS达到最小，而是a 使QTLS达到最小;其次，以上方法考虑的是x和y 坐标上的误差，没有直接考虑到故障间隔时间t和故障发生概率F 的误差，但两组误差并不相同.

针对以上问题，在上节得到β 的估计值β(γ)的基础上，对α 和γ 的估计过程如下.

将β(γ)代入原三参数威布尔分布函数(式(1))中，此时分布函数F^中只含有未知参数α 与γ. 然后根据以下约束准则进行优化:

QITLS实际上是威布尔分布函数F 在节点ti处的偏差平方和.代入β 表达式β(γ)后QITLS为自变量α与γ 的函数，结合表1 中数据，应用差分进化算法对α 与γ 进行优化使QITLS达到最小值.

根据以上分析过程可知，改进的整体最小二乘法继承了整体最小二乘法的特点，即兼顾自变量和因变量的误差.根据三参数威布尔分布函数线性变换的特点，为给出参数α 的更佳估计值，改进的整体最小二乘法以F 的偏差平方和最小作为约束准则进行优化，考虑到了初始数据t 和F 的误差，使所得曲线在拟合效果上达到最佳.

3 差分进化算法

差分进化(DE)算法最早由Storn 等［10］在1995年提出.DE 可以动态跟踪当前的搜索情况以调整搜索策略，具有较好的全局收敛能力和鲁棒性，适用于一些常规的数学规划方法所无法求解的复杂环境中的优化问题.

3.1 DE 算法的基本步骤

在使用DE 算法之前需要设定以下几个参数:种群大小N，即群体中个体的数量;最大迭代代数G;变异因子M，通常取值范围为M∈［0，2］;交叉概率C，通常取值范围为C∈［0，1］. 实施DE 算法时包含以下4个基本步骤［11］:

(1)初始化.个体初始化公式为

式中，r1，r2，r3∈{1，2，…，N}，为互不相同的整数且r1，r2，r3与i 不同，因此种群规模N≥4.

rand(j)∈［0，1］，为均匀分布随机数;rand(i)∈［1，2，…，D］，是随机选择的维数变量索引.

(4)选择操作.以最小化优化为例，选择操作为

3.2 DE 对三参数威布尔模型的优化

根据2.2 节分析结果，结合表1 中数据，利用DE对参数α 与γ 进行优化，使式(9)达到最小值.文中的DE 优化过程是通过在Matlab 中编程实现的.其优化α 与γ 的基本流程图如图1 所示，具体步骤如下:

图1 DE 优化流程图Fig.1 Process of DE optimization

(1)初始化种群. 根据文献［12-15］研究结论，文中设置种群规模N =10;最大迭代代数G =150;变异因子M=0.5;交叉概率C=0.9;count=0;因为有α 与γ 两个参数，所以每个个体有两个分量，即种群维数D=2;γ 的取值应该为0≤γ≤tmin，根据该加工中心故障数据(见表1)，γ 的界限为γL=0，γU=18.29;αL=0，αU=500. 根据式(10)对α 与γ 进行个体初始化.

(2)变异、交叉.根据式(11)、(12)对种群中的每个个体进行变异和交叉操作.

(3)选择.选取式(9)为目标函数，根据式(13)进行下一子代的选择操作.

(4)判断. 如果达到规定的进化代数则停止程序输出结果，否则返回变异操作步骤继续运行.

根据以上流程和表1 中数据，计算得到该系列加工中心可靠性三参数威布尔分布模型在γ=13.7816、α=406.6531时QITLS达到最小，最小值为0.0402，根据式(8)有β=0.9720.故得到该加工中心可靠性三参数威布尔分布函数为

4 参数估计结果对比

根据普通最小二乘法和文中改进的整体最小二乘法分别得到两个不同的三参数威布尔分布模型，对这两个模型的相关参数进行对比分析，结果如表2所示.

表2 所得模型相关参数对比Table 2 Comparison of parameters of models

从表2 可以看出:对于尺度参数α，改进的整体最小二乘法得到的α 比最小二乘法得到的α 有所增大;对于形状参数β，由改进的整体最小二乘法得到的β 比最小二乘法有所增大，β 值更加接近1，说明该型号加工中心处于偶然故障期;对于位置参数γ，由改进的整体最小二乘法得到的γ 小于由最小二乘法得到的γ 值.

对于QITLS，改进的整体最小二乘法所得QITLS小于最小二乘法所得QITLS，这说明由改进的整体最小二乘法所得三参数威布尔分布函数F 的偏差平方和最小，所得分布函数F 的曲线拟合优度最佳.

5 结语

针对以往加工中心可靠性建模中存在的一些问题，提出了一种对三参数威布尔分布模型进行建模的新方法，并且结合某型号加工中心故障数据进行计算.最后通过对比分析可以发现文中所提出的方法具有以下特点:所采用的DE 算法具有快速收敛、全局收敛、指导记忆能力、稳定性、通用性等特点，较普通的迭代算法更适合解决这类复杂目标函数的寻优问题;改进的整体最小二乘法兼顾到了故障间隔时间t 和故障发生概率F 的误差，使分析结果更加合理.改进的整体最小二乘法以三参数威布尔分布函数F 在节点ti处的偏差平方和作为约束准则进行优化，所得曲线在整体上最大程度地接近原始数据点，对三参数威布尔分布曲线拟合效果最好. 另外，应用文中所提出的方法对该型号加工中心的三参数威布尔分布进行了参数估计，根据所得结果，可以发现形状参数β 近似等于1，说明该型号加工中心处于偶然故障期.

［1］申桂香，孟书，张英芝，等. 平均秩次法在子系统可靠性建模中的应用［J］.吉林大学学报:工学版，2014，44(1):101-105.Shen Gui-xiang，Meng Shu，Zhang Ying-zhi，et al. Application of average rank time method in reliability modeling for subsystems［J］. Journal of Jilin University:Engineering and Technology Edition，2014，44(1):101-105.

［2］王晓峰，申桂香，张英芝，等. 可靠性模型参数估计方法的对比［J］.华南理工大学学报:自然科学版，2011，39(6):47-52.Wang Xiao-feng，Shen Gui-xiang，Zhang Ying-zhi，et al.Comparison of parameter estimation methods for reliability model［J］.Journal of South China University of Technology:Natural Science Edition，2011，39(6):47-52.

［3］Hideki Nagatsuka，Toshinari Kamakura.A consistent method of estimation for the three-parameter Weibull distribution［J］. Computational Statistics and Data Analysis，2013，58:210-226.

［4］Moeini Asghar，Jenab Kouroush.Fitting the three parameter Weibull distribution with ross entropy ［J］. Applied Mathematical Modelling，2013，37:6354-6363.

［5］孙丽玢，汤银才. 在定时截尾样本下三参数威布尔分布的矩估计［J］. 数学的实践与认识，2010，40(22):156-161.Sun Li-bin，Tang Yin-cai.Moment estimation of the threeparameter Weibull distribution under type I censored samples［J］. Mathematics in Practice and Theory，2010，40(22):156-161.

［6］王晓峰.加工中心可信性影响度分析及增长技术研究［D］.长春:吉林大学机械科学与工程学院，2012.

［7］李雄军. 对X 和Y 方向最小二乘线性回归的讨论［J］.误差与数据处理，2005(1):50-52.Li Xiong-jun.The least squares discussion of X and Y direction ［J］. Modern Manufacturing Engineering，2005(1):50-52.

［8］Golub G H，Van Loan C F.An analysis of the total least squares problem［J］.SIAM Journal on Numerical Analysis，1980，17(6):883-893.

［9］Byungsoo Kim.Total least square method applied to rating curves［J］.Hydrol Process，2014，28(13):4057-4066.

［10］Storn R，Price K. Differential evolution—a simple and efficient adaptive scheme for global optimization over continuous spaces ［R］. Berkley:Technical Report International Computer Science Institute，1995.

［11］Neri Ferrante，Tirronen Ville. Recent advances in differential evolution:a survey and experimental analysis［J］.Artif Intell Rev，2010，33(1):61-106.

［12］吴亮红.差分进化算法及应用研究［D］. 长沙:湖南大学电气与信息工程学院，2007.

［13］Qin A K，Huang V L，Suganthan P N.Differential evolution algorithm with strategy adaptation for global numerical optimization［J］.IEEE Trans Evol Comput，2009，13(2):398-417.

［14］张宝珍，张尧，林凌雪. 基于改进差分进化算法的估计等值法［J］. 华南理工大学学报:自然科学版，2014，42(4):7-12.Zhang Bao-zhen，Zhang Yao，Lin Ling-xue. Estimation equivalence method based on modified differential evolution algorithm［J］. Journal of South China University of Technology:Natural Science Edition，2014，42(4):7-12.

［15］魏德敏，陈贵涛.自适应子群体米母算法及其在混凝土框架位移性能优化中的应用［J］.华南理工大学学报:自然科学版，2013，41(10):108-116.Wei De-min，Chen Gui-tao. Self-daptive subpopulationsbased memetic algorithm and its application to displacement performance optimization of concrete frames［J］.Journal of South China University of Technology:Natural Science Edition，2013，41(10):108-116.