基于聚类多目标分布估计算法的齿轮减速器优化设计方法

摘要

基于聚类多目标分布估计算法的齿轮减速器优化设计方法，本发明涉及基于聚类多目标分布估计算法的齿轮减速器优化设计方法。解决现有的多目标分布估计算法在求解多目标优化问题的过程中存在没有充分利用算法的局部搜索能力，求解过程中直接丢弃异常解，种群多样性容易丢失，过多的计算开销用于构建最优概率模型的问题。本发明首先利用凝聚层次聚类算法将种群划分为若干个局部类，从每一个局部类中随机选择一个个体构成一个全局类，然后为每个个体构建一个高斯模型去逼近种群结构，并抽样产生新个体；此高斯模型的均值为个体本身，协方差矩阵为个体所在局部类的协方差矩阵或者是全局类的协方差矩阵。本发明用于航天领域。

说明书

技术领域

本发明涉及齿轮减速器优化设计方法。

背景技术

实际工程中存在着大量的具有多约束、多变量以及非线性等性质的复杂多目标优化问题(Multiobjective Optimization Problem,MOP)。典型的约束MOP表达如下(王勇,蔡自兴,周育人,等.约束优化进化算法[J].软件学报,2009,20(1):11-29)：

minF(x)＝[f₁(x),f₂(x),...,f_m(x)]^T

x＝(x₁,x₂,…,x_n)^T∈Ω

其中，x为n维决策变量向量；F(x)为m维目标函数向量；p为不等式约束条件g_i(x)的个数；q为等式约束条件h_j(x)的个数。Ω为决策空间。

由于大多数情况下MOP中各子目标之间相互冲突，不存在一个最优解使所有子目标同时达到最优。因此，不同于单目标优化问题只有一个或者若干个孤立的最优解，MOP具有大量的对于所有目标都可以接受的折衷解，即Pareto最优解。所有Pareto最优解组成的集合称为Pareto解集(Pareto Set,PS)，Pareto解集投影到目标空间获得的目标向量的集合称为Pareto前沿(Pareto Front,PF)。并且连续MOP的PS和PF的结构具有规则特性，即根据Karush-Kuhn-Tucker条件，在宽松的条件下，具有m个目标的连续MOP的PS(或PF)的结构是一个m-1维的分段连续的流型。对于一个MOP，由于不可能求解出其所有的Pareto最优解，因此在求解过程中，决策者往往希望获得一个有限数目的逼近解的集合(逼近解集)，其对应的目标向量(构成逼近前沿)越靠近PF越好(收敛性)，并且沿着PF分布越广泛以及越均匀越好(多样性)。

由于传统的确定性优化技术不能较好地对复杂的MOP进行求解，因此基于自然启发搜索的全局优化算法——演化算法(Evolutionary Algorithm,EA)成为了解决MOP的流行的方法。多目标演化算法(Multiobjective Evolutionary Algorithm,MOEA)具有良好的并行性、鲁棒性，而且其求解不依赖于问题特性、通用性强，并且单次运行就可获得MOP的Pareto解集的逼近，近年来得到了蓬勃发展(Zhou A,Qu B Y,Li H,et al.Multiobjectiveevolutionary algorithms:A survey of the state of the art[J].Swarm&Evolutionary Computation,2011,1(1):32-49)。

在EA当中，包含有个体重组和环境选择两个重要组成部分。个体重组用于产生新解，而环境选择则负责为下一代挑选有效的新解。在MOEA中，根据新解产生的方式，重组算子可以粗略地划分为两大类，即基于遗传(Genetic-Based)的重组算子和基于模型(Model-Based)的重组算子。基于遗传的MOEA应用传统的重组算子(例如：模拟二进制交叉(Deb K,Beyer H G.Self-adaptive genetic algorithms with simulated binary crossover[J].Evolutionary Computation,2001,9(2):197-221)、多项式变异(Schaer JD.Multipleobjective optimization with vector evaluated genetic algorithms[C].Proceedings of the 1st International Conference on Genetic Algorithms andTheir Applications,Lawrence Erlbaum Associates,1985,93-100.)等)产生新解。基于模型的MOEA采用概率模型描述群体中个体的分布，并通过建立的模型采样产生新个体，多元高斯模型、Bayesian网络、流型学习、密度估计等常用的机器学习方法被广泛应用于群体建模(MartíL,Grimme C,Kerschke P,et al.Averaged Hausdorff approximations ofpareto fronts based on multiobjective estimation of distribution algorithms[C].Proceedings of the Companion Publication of the 2015Annual Conference onGenetic and Evolutionary Computation.ACM,2015:1427-1428)。目前已有的MOEA大多数采用的是基于遗传的新解产生方法，但是基于模型的MOEA也得到了学者们越来越多的关注，近些年变得流行的多目标分布估计算法(Multiobjective Estimation ofDistribution Algorithm,MEDA)就是一个重要的代表(Pelikan M,Sastry K,Goldberg DE.Multiobjective estimation of distribution algorithms[C].ScalableOptimization via Probabilistic Modeling.Heidelberg,Berlin,Germany:Springer-Verlag,2006:223-248)。

分布估计算法(Estimation of Distribution Algorithm，EDA)(P,Lozano J A.Estimation of distribution algorithms:A new tool for evolutionarycomputation[M].SpringerScience&Business Media,2002)是EA中一类特别的方法。EDA并不采用传统的交叉变异等遗传操作，取而代之，它从所挑选的有效解中明确地提取全局统计信息，基于提取的统计信息建立一个有效解后验概率分布模型，进而从建立的模型中抽样产生新解。在基于遗传的MOEA中，遗传操作可能会破坏种群强模式的建立，因此种群个体朝向最优解的移动方向非常难以预测。然而，MEDA能够预测PF的位置或模式，或者预测在搜索空间中有效的搜索方向。通过调整搜索使其沿着发掘或预测的有效的搜索方向，就能够较好地产生有效解。学者们已经提出了各种MEDAs，并且这些算法显示出了良好的性能。

虽然MEDA已经得到了越来越多学者的关注和研究，但是依然存在着不足，典型的有：算法中没有充分考虑MOP的规则特性，种群中的异常解处理不正确，种群多样性容易丢失，以及过多的计算开销用于构建最优的种群模型(MartíL,Grimme C,Kerschke P,etal.Averaged Hausdorff approximations of pareto fronts based on multiobjectiveestimation of distribution algorithms[C].Proceedings of the CompanionPublication of the 2015Annual Conference on Genetic and EvolutionaryComputation.ACM,2015:1427-1428)(Zhang Q,Zhou A,Jin Y.RM-MEDA:A regularitymodel-based multiobjective estimation of distribution algorithm[J].IEEETransactions on Evolutionary Computation,2008,12(1):41-63)。针对上述不足，本发明提出一种基于聚类的新型多目标分布估计算法(Clustering Based MEDA,CEDA)。在CEDA中的每一代，首先利用聚类算法发掘种群中个体的分布结构，然后基于结构信息，为每一个个体构建一个多元高斯模型(Multivariate Gaussian Model,MGM)，基于此模型，抽样产生新解。

发明内容

本发明的目的是为了解决现有的多目标分布估计算法在求解多目标优化问题的过程中存在没有充分利用算法的局部搜索能力，求解过程中直接丢弃异常解，种群多样性容易丢失，过多的计算开销用于构建最优概率模型的问题，而提出基于聚类多目标分布估计算法的齿轮减速器优化设计方法。

基于聚类多目标分布估计算法的齿轮减速器优化设计方法的具体步骤为：

步骤一：初始化种群P＝{x¹,x²,…,x^N}和控制概率β，设置演化代数t＝0；x¹,x²,…,x^N为种群中的个体；

步骤二：进行主循环；

步骤二一：设置一个空的外部文档A＝φ；

步骤二二：对种群P进行聚类，{LC¹,…,LC^K}＝AHC(P,K)；AHC为凝聚层次聚类算法，K为AHC中定义的最大聚类个数，LC¹,…,LC^K为聚类得到的K个局部类；

步骤二三：构建一个全局类GC；

步骤二四：分别计算局部类LC^k和全局类GC的协方差矩阵Σ^k(k＝1,…,K)和Σ^GC；

步骤二五：新解产生；

步骤二六：环境选择：更新种群P＝EnvSel(Α∪P)；

步骤二七：令t＝t+1；

步骤二八：如果t>T算法结束，输出P；否则转向步骤二；所述T为最大演化代数；

步骤三：停机，输出P。

关于CEDA，有如下说明：

(1)Jin等(Jin Y,Sendhoff B.Connectedness,regularity and the success oflocal search in evolutionary multi-objective optimization[C].Proceedings ofIEEE Congress on Evolutionary Computation.IEEE,2003,3:1910-1917)指出在MOEA中，相似个体重组，能提高产生新解的质量。产生此效果的原因是增强了算法的局部搜索，暗含地利用了MOP的规则特性。与此类似，本文的CEDA采用邻近个体为每个当前个体构建高斯模型逼近种群结构进而抽样产生新解，这一机制同样增强了算法的局部搜索，充分地考虑了MOP的规则特性，理应也能更好地产生高质量的新解。

(2)与RM-MEDA中利用局部主成分分析方法提取流型结构，然后抽样产生新解的方式相比，CEDA中的基于聚类建立高斯模型抽样产生新解的方式更简单易用。并且，在进化的早期阶段，PS的流型结构尚未显现，种群需要良好多样性，但是RM-MEDA的新解产生方式限制了新解的产生方向，不利于产生多样化的解，而CEDA利用完全协方差矩阵抽样产生新解，能从各个方向生成新解，更好地维护新解的多样性。

(3)MIEDA(Bosman PA,Thierens D.Multi-objective optimization withdiversity preserving mixture-based iterated density estimation evolutionaryalgorithms[J].International Journal of Approximate Reasoning,2002,31(3):259-289)等传统的为每一个类构建一个高斯模型进行抽样的方式，其产生的新解大量地分布在均值向量附近，新解的多样性不够，而CEDA为每个种群个体以自身为均值建立一个高斯模型抽样产生新解，实际上是为每个个体添加一个高斯扰动，此种方式能产生更为多样化的解。

(4)在为个体构建高斯模型的时候，如果对于每个个体都计算协方差矩阵，则需要大量的建模计算开销。为了解决此问题，CEDA中使同一类中的个体共享相同的协方差矩阵进行建模从而大大降低建模计算开销。之所以能够进行此策略是因为相似的个体理应具有相近的高斯模型，并且近似的高斯模型就已经满足算法的要求，没有必要花费大量的计算开销建立精确的模型。

(5)不同于以往丢弃异常解的建模方式，CEDA中为每个个体建立一个高斯模型进行抽样，实际上就是充分地考虑了异常解代表的解空间区域，因此能更好地加强对解空间的搜索。

本发明的有益效果为：

本发明设计了一种基于聚类的新型多目标分布估计算法(CEDA)。在CEDA中，首先利用凝聚层次聚类算法将种群划分为若干个局部类，从每一个局部类中随机选择一个个体构成一个全局类，然后为每个个体构建一个高斯模型(此高斯模型的均值为个体本身，协方差矩阵为个体所在局部类的协方差矩阵或者是全局类的协方差矩阵)去逼近种群结构，并抽样产生新个体。此新解产生方法充分地考虑了多目标优化问题的规则特性，其本质是为每个个体添加了一个外部扰动，能改善已有的大部分分布估计算法中存在的异常个体处理不合理，种群多样性容易丢失的问题。并且处于同类中的个体共享协方差矩阵用于建模，极大地降低了建模的计算开销。

以具有复杂Pareto前沿和复杂Pareto解集形状的多目标优化测试题为求解对象，将CEDA与典型的MOEAs进行了对比实验。实验结果表明，CEDA对于此类问题具有优异的求解性能。将CEDA算法应用于齿轮减速器优化设计中，结果表明，CEDA算法同样能够快速有效地求解此类复杂的实际工程问题。

附图说明

图1为对GLT1测试中获得的平均IGD指标值进化曲线图；图中1,2,3,4,5,分别代表CEDA,NSGA-II,SMS-EMOA,RM-MEDA,TMOEA/D五种算法；GLT为标准测试题名称；

图2为对GLT2测试中获得的平均IGD指标值进化曲线图；图中1,2,3,4,5,分别代表CEDA,NSGA-II,SMS-EMOA,RM-MEDA,TMOEA/D五种算法；

图3为对GLT3测试中获得的平均IGD指标值进化曲线图；图中1,2,3,4,5,分别代表CEDA,NSGA-II,SMS-EMOA,RM-MEDA,TMOEA/D五种算法；

图4为对GLT4测试中获得的平均IGD指标值进化曲线图；图中1,2,3,4,5,分别代表CEDA,NSGA-II,SMS-EMOA,RM-MEDA,TMOEA/D五种算法；

图5为对GLT5测试中获得的平均IGD指标值进化曲线图；图中1,2,3,4,5,分别代表CEDA,NSGA-II,SMS-EMOA,RM-MEDA,TMOEA/D五种算法；

图6为对GLT6测试中获得的平均IGD指标值进化曲线图；图中1,2,3,4,5,分别代表CEDA,NSGA-II,SMS-EMOA,RM-MEDA,TMOEA/D五种算法；

图7为对GLT1测试中，TMOEA/D获得的全部逼近前沿图；图中横坐标为目标1值，纵坐标为目标2值；

图8为对GLT1测试中，CEDA获得的全部逼近前沿图；图中横坐标为目标1值，纵坐标为目标2值；

图9为对GLT2测试中，TMOEA/D获得的全部逼近前沿图；

图10为对GLT2测试中，CEDA获得的全部逼近前沿图；

图11为对GLT3测试中，TMOEA/D获得的全部逼近前沿图；

图12为对GLT3测试中，CEDA获得的全部逼近前沿图；

图13为对GLT4测试中，TMOEA/D获得的全部逼近前沿图；

图14为对GLT4测试中，CEDA获得的全部逼近前沿图；

图15为对GLT5测试中，TMOEA/D获得的全部逼近前沿图；图中三个坐标分别代表目标1,2,3值；

图16为对GLT5测试中，CEDA获得的全部逼近前沿图；

图17为对GLT6测试中，TMOEA/D获得的全部逼近前沿图；

图18为对GLT6测试中，CEDA获得的全部逼近前沿图；

图19为对GLT1测试中，TMOEA/D获得的代表性逼近前沿图；

图20为对GLT1测试中，CEDA获得的代表性逼近前沿图；

图21为对GLT2测试中，TMOEA/D获得的代表性逼近前沿图；

图22为对GLT2测试中，CEDA获得的代表性逼近前沿图；

图23为对GLT3测试中，TMOEA/D获得的代表性逼近前沿图；

图24为对GLT3测试中，CEDA获得的代表性逼近前沿图；

图25为对GLT4测试中，TMOEA/D获得的代表性逼近前沿图；

图26为对GLT4测试中，CEDA获得的代表性逼近前沿图；

图27为对GLT5测试中，TMOEA/D获得的代表性逼近前沿图；

图28为对GLT5测试中，CEDA获得的代表性逼近前沿图；

图29为对GLT6测试中，TMOEA/D获得的代表性逼近前沿图；

图30为对GLT6测试中，CEDA获得的代表性逼近前沿图；

图31为重组控制概率(β)分析；

图32为聚类数目(K)分析；

图33为齿轮减速器结构简图；

图34为NSGA-II、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA算法对齿轮减速器优化设计模型分别独立运算33次获得的HV指标值的箱型图；图中横坐标1，2，3,4,5分别代表NSGA-II、SMS-EMOA、RM-MEDA、TMOEA/D、CEDA算法；

图35为NSGA-II、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA算法对齿轮减速器优化设计模型分别独立运算33次获得的HV指标值的箱型图的局部放大图；

图36为平均HV指标值演化曲线；图中1,2,3,4,5,分别代表CEDA,NSGA-II,SMS-EMOA,RM-MEDA,TMOEA/D五种算法；

图37为对齿轮减速器测试中，NSGA-II获得的全部逼近前沿图；图中横坐标为目标1的值，纵坐标为目标2的值；

图38为对齿轮减速器测试中，NSGA-II获得的代表性逼近前沿图；

图39为对齿轮减速器测试中，CEDA获得的全部逼近前沿图；图中横坐标为目标1的值，纵坐标为目标2的值；

图40为对齿轮减速器测试中，CEDA获得的代表性逼近前沿图。

具体实施方式

具体实施方式一：基于聚类多目标分布估计算法的齿轮减速器优化设计方法的具体过程为：

EDA已经被大量地应用于MOP的求解。Bosman和Thierens(Bosman PA,ThierensD.Multi-objective optimization with diversity preserving mixture-basediterated density estimation evolutionary algorithms[J].International Journalof Approximate Reasoning,2002,31(3):259-289)提出了一种基于混合的多目标迭代密度估计演化算法(MIEDA)，用于求解连续和离散优化问题，MIEDA被认为是其它MEDAs的基准算法。Pelikan等(Pelikan M,Sastry K,Goldberg D E.Multiobjective hBOA,clustering,and scalability[C].Proceedings of the 7th Annual Conference onGenetic and Evolutionary Computation.ACM,2005:663-670)采用基于支配的框架并且使用K-means聚类算法进行建模，设计了一种多目标分层Bayesian优化算法(mohBOA)。Sastry等(Sastry K,Goldberg D E,Pelikan M.Limits of scalability ofmultiobjective estimation of distribution algorithms[C].Proceedings of IEEECongress on Evolutionary Computation.IEEE,2005,3:2217-2224)提出了一种延伸的紧凑遗传算法(ECGA)以解决可扩展的欺骗问题。Shim等(Shim V A,Tan K C,Cheong C Y.Ahybrid estimation of distribution algorithm with decomposition for solvingthe multiobjective multiple traveling salesman problem[J].IEEE Transactionson Systems,Man,and Cybernetics,Part C(Applications and Reviews),2012,42(5):682-691)将EDA集成到基于分解的MOEA框架中，提出了一种混合局部搜索元启发式方法的基于分解的EDA。Cheng等(Cheng R,Jin Y,Narukawa K,et al.A multiobjectiveevolutionary algorithm using Gaussian process-based inverse modeling[J].IEEETransactions on Evolutionary Computation,2015,19(6):838-856)构建基于高斯过程的逆模型(Inverse Models)将所有已发掘的非支配解从目标空间映射到决策空间，通过对目标空间抽样产生新解，提出了基于逆模型的MOEA(IM-MOEA)。

为了利用连续MOP的规则特性提高MEDA的求解性能，学者们又提出了多种基于规则特性的MEDA。Zhou等提出了一种基于规则模型的MEDA，即(RM-MEDA)，其使用局部主成分分析方法对有效解建立概率分布模型，并通过概率模型抽样产生新个体。在此之后，又设计了一种基于概率模型的MOEA(MMEA)(Zhou A,Zhang Q,Jin Y.Approximating the set ofPareto-optimal solutions in both the decision and objective spaces by anestimation of distribution algorithm[J].IEEE Transactions on EvolutionaryComputation,2009,13(5):1167-1189)，在决策空间和目标空间同时建立概率模型逼近PS和PF。受到RM-MEDA思想的启发，出现了一系列变种的RM-MEDA，例如，基于冗余类消减的MEDA(Wang Y,Xiang J,Cai Z.A regularity model-based multiobjective estimationof distribution algorithm with reducing redundant cluster operator[J].AppliedSoft Computing,2012,12(11):3526-3538)，带有局部学习策略的MEDA(Li Y,Xu X,Li P,et al.Improved RM-MEDA with local learning[J].Soft Computing,2014,18(7):1383-1397)，基于流型学习的演化多目标优化算法等(Li K,Kwong S.A general framework forevolutionary multiobjective optimization via manifold learning[J].Neurocomputing,2014,146:65-74)。

目前为止，在已有的MEDAs中，大多数在设计过程中并没有充分地考虑MOP的规则特性，并且建模过程中，只是运用较少个数的高斯模型描述有效解的分布，往往丢弃了异常解。实际上，运用较少个数的高斯模型抽样产生新解，新解大量的集中在模型的中心位置(均值)附近，多样性不足，并且异常解可能代表着新的有效区域，需要开展搜索。另外，在已有的基于规则特性的MEDAs中，大部分借鉴RM-MEDA的思想，而RM-MEDA建模复杂，在搜索的早期阶段种群的多样性保持不好，并且难以设定主成分的个数。为了改善前述问题，本发明充分考虑MOP的规则特性，计划基于聚类运用较多数目的简单的高斯模型逼近种群结构，进而抽样产生新解，从而降低算法结构的复杂性，增强算法的易用性，并提高算法产生多样解的能力。

CEDA采用凝聚层次聚类(Agglomerative Hierarchical Clustering,AHC)算法(Xu R,Wunsch D.Clustering[M].John Wiley&Sons,Hokoben,New Jersey,2008)发掘种群结构。CEDA的基本框架为以下步骤。

步骤一：初始化种群P＝{x¹,x²,…,x^N}和控制概率β，设置演化代数t＝0；x¹,x²,…,x^N为种群中的个体；

步骤二：进行主循环；

步骤二一：设置一个空的外部文档A＝φ；

步骤二二：对种群P进行聚类，{LC¹,…,LC^K}＝AHC(P,K)；AHC为凝聚层次聚类算法，K为AHC中定义的最大聚类个数，LC¹,…,LC^K为聚类得到的K个局部类；

步骤二三：构建一个全局类GC；

步骤二四：分别计算局部类LC^k和全局类GC的协方差矩阵Σ^k(k＝1,…,K)和Σ^GC；

步骤二五：新解产生；

步骤二六：环境选择：更新种群P＝EnvSel(Α∪P)；

步骤二七：令t＝t+1；

步骤二八：如果t>T算法结束，输出P；否则转向步骤二；所述T为最大演化代数；

步骤三：停机，输出P。

在算法，N代表种群大小，K为AHC中定义的最大聚类个数，T为最大演化代数，GC和LC^k分别代表全局类和第K个局部类，为xⁱ所在的局部类的协方差矩阵，β表示利用LC^k建立抽样池的概率(称之为重组控制概率)，rand()生成一个在[0,1]区间内均匀分布的随机数。

在CEDA算法的每一代中，首先利用AHC将种群划分为K个局部类(步骤二二)，并从每一个局部类中随机抽取一个个体共同构建一个全局类(步骤二三)。然后计算出全局类和所有局部类的协方差矩阵Σ^GC和Σ^k(k＝1,…,K)(步骤二四)。紧接着为每个个体xⁱ确定一个协方差矩阵Σⁱ，该协方差矩阵分别以β和1-β的概率设置为Σ^k或Σ^GC(步骤二五一)，并且由xⁱ和Σⁱ构成高斯模型抽样产生新个体yⁱ(步骤二五二)，将yⁱ存于外部档案A中(步骤二五三)。最后基于A和P利用环境选择机制更新种群P(步骤二六)。以下内容对CEDA的细节进行详细介绍。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤二二中AHC(P,K)具体为：

利用AHC算法将种群P中包含的N个个体，即P＝{x¹,x²,…,x^N}，划分到K个类中的原理为以下步骤。

(1)将种群P中每个个体作为一个类；

(2)进行循环：

(2.1)计算每两个不同的类之间的欧氏距离；

(2.2)找出距离最小的两个类合并成新的类；

(2.3)判断是否满足终止条件，即聚类个数是否大于K，若满足则终止，输出最终的聚类结果，否则转至步骤(2.1)。

AHC首先将每一个个体视为一个类，然后利用一系列机制合并不同类，直至种群聚类个数不大于K。在CEDA的AHC算法中利用组间平均联接算法(Group average linkagealgorithm)定义两个类之间的距离。关于AHC算法的细节内容参考文献(Xu R,WunschD.Clustering[M].John Wiley&Sons,Hokoben,New Jersey,2008)。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述步骤二五中新解产生的具体过程为：

对于每一个个体xⁱ∈P,i＝1,…,N进行如下步骤：

步骤二五一：为个体xⁱ选择一个协方差矩阵Σⁱ；

其中所述为个体xⁱ所在的局部类的协方差矩阵，Σ^GC为全局类的协方差矩阵；

步骤二五二：产生新个体yⁱ＝SolGen(Σⁱ,xⁱ)；

步骤二五三：保留新解A＝A∪{yⁱ}。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述SolGen(Σⁱ,xⁱ)具体为：

步骤二五二产生新个体，此过程包含基于多元高斯模型的抽样和多项式变异，其具体为以下步骤。

(1)利用平方根法(Cholekey分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解)分解协方差矩阵Σⁱ得到一个下三角矩阵Λ，并且Σⁱ＝ΛΛ^T；

(2)产生向量v＝(v₁,…,v_n)^T，其中v_j～N(0,I),j＝1,…,n服从高斯分布；

(3)产生试验解y'＝xⁱ+Λv,y'＝(y'₁,…,y'_n)^T；

(4)修复该试验解：

a_j和b_j代表第j个变量的上下边界；

(5)对试验解进行变异：

其中

p_m为变异概率，η_m为变异指数，r＝rand()；rand()是随机产生一个0-1之间的数，r＝rand()就是随机产生一个0-1之间的数赋值给r。

(6)修复个体j＝1,2,...,n；j代表每一个个体中变量的个数，因为每一个个体相当于是一个多维向量，具有j个数组成。

(7)返回新解

对于种群中的每个个体，首先基于协方差矩阵为其建立多元高斯模型并抽样产生一个初始试验解(步骤(1)-步骤(3))。然后对试验解进行修补，保证其可行性(步骤(4))，紧接着对试验解进行变异操作以增强解的多样性(步骤(5))，最后再次对试验解进行边界修补确保其可行(步骤(6))。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤二六中EnvSel(A∪P)的具体过程为：

每一代当新解产生完毕之后，运用EnvSel(A∪P)从A∪P中选出优秀个体进入下一代演化过程。CEDA采用SMS-MOEA(Beume N,Naujoks B,Emmerich M.SMS-EMOA:Multiobjective selection based on dominated hypervolume[J].European Journalof Operational Research,2007,181(3):1653-1669)中提出的基于超体积指标的环境选择方法。超体积指标是已知的唯一一个为“Pareto兼容(Pareto compliant)”的一元指标(Zitzler E,Thiele L,Laumanns M,et al.Performance assessment of multiobjectiveoptimizers:an analysis and review[J].IEEE Transactions on EvolutionaryComputation,2003,7(2):117-132)，基于超体积指标的环境选择方法在求解具有复杂PF的MOP时展现出了良好的性能(Zhang H,Zhou A,Song S,Zhang Q,Gao X.Z.,Zhang J.Aself-organizing multiobjective evolutionary algorithm[J],IEEE Transactions onEvolutionary Computation,2016,in press)。EnvSel(A∪P)具体为：

(1)利用快速非支配排序方法对A∪P中的个体进行排序；

{B¹,…,B^L}＝Fast_Nondominated_Sort(A∪P)；

B¹,…,B^L为L个不同的非支配前沿；Fast_Nondominated_Sort为快速非支配排序方法，为一种现有算法。

(2)复制较优的个体到辅助种群P'中

(3)如果l>1，进行循环；当|P'|>N时，进行以下步骤：

(3.1)挑选出x^*，其中d(x,P')指的是在P'中支配x的点的个数；

(3.2)将x^*从P'中移除，P'＝P'{x^*}；

(4)如果l＝1(l的值为大于等于1的正整数，等于1的时候为步骤(4))，进行循环：当|P'|>N时，进行以下步骤：

(4.1)挑选出x^*，其中为x的超体积贡献度；

(4.1)将x^*从P'中移除，P'＝P'{x^*}；

(5)将P'赋给P，P＝P'；

(6)输出P。

首先将当前种群P和外部文档A合并成一个新的种群，并且利用NSGA-II(Deb K,Pratap A,Agarwal S,et al.A fast and elitist multiobjective genetic algorithm:NSGA-II[J].IEEE Transactions on Evolutionary Computation,2002,6(2):182-197)中提出的快速非支配排序方法将新种群中的个体划分到L个不同的非支配前沿{B¹,…,B^L}当中。然后根据排序的结果，将新种群中的较优个体复制到一个辅助种群P'当中，直到|P'|≥N。如果P'当中包含有多个非支配前沿(即l>1)_，则逐个移除第l前沿中d(x,P')最大的个体，直到|P'|＝N，d(x,P')指的是在P'中支配x的点的个数；否则如果l＝1，逐个移除P'中超体积贡献度最小的个体，直到|P'|＝N，超体积贡献度的计算方法参考文献(Beume N,Naujoks B,Emmerich M.SMS-EMOA:Multiobjective selection based ondominated hypervolume[J].European Journal of Operational Research,2007,181(3):1653-1669)。最后，将P'赋给P，作为下一代的种群。

其它步骤及参数与具体实施方式一至四之一相同。

采用以下实施例验证本发明的有益效果：

实施例一：

1、实验分析

标准测试实例和性能度量指标

为了测试CEDA的性能，首先利用标准测试题对其进行测试。大多数工程中的MOP具有复杂的PF结构，因此CEDA算法理应适用于求解此类具有复杂PF结构的MOPs。本文利用具有复杂PF和PS结构的6道标准测试题GLT1-GLT6对CEDA进行测试。其中，GLT1-GLT4为双目标测试问题，GLT5-GLT6为三目标测试问题。GLT测试题的具体细节参考文献(Zhang H,ZhouA,Song S,Zhang Q,Gao X.Z.,Zhang J.A self-organizing multiobjectiveevolutionary algorithm[J],IEEE Transactions on Evolutionary Computation,2016,in press)。

为了评估算法的性能，运用两个常用的性能指标，即反世代距离IGD(Zhang Q,Zhou A,Jin Y.RM-MEDA:A regularity model-based multiobjective estimation ofdistribution algorithm[J].IEEE Transactions on Evolutionary Computation,2008,12(1):41-63)(Zhou A,Zhang Q,Jin Y,et al.A model-based evolutionary algorithmfor bi-objective optimization[C].Proceedings of IEEE Congress on EvolutionaryComputation.IEEE,2005,3:2568-2575)和超体积HV(Zitzler E,ThieleL.Multiobjective evolutionary algorithms:a comparative case study and thestrength Pareto approach[J].IEEE Transactions on Evolutionary Computation,1999,3(4):257-271)，度量算法获得的逼近前沿的效果。IGD和HV是两个均能够综合评价获得的逼近前沿的收敛性与多样性的性能指标。并且IGD值越小、HV值越大代表算法所求得的逼近前沿的收敛性和多样性越好。

在接下来的实验中，计算HV指标值时，各测试题的参考点取值为：GLT1取r＝(2,2)^T,GLT2取r＝(2,11)^T,GLT3取r＝(2,2)^T,GLT4取r＝(2,3)^T,GLT5-GLT6取r＝(2,2,2)^T。

对比算法及参数设置

选取四种典型的MOEAs，即NSGA-II(Deb K,Pratap A,Agarwal S,et al.A fastand elitist multiobjective genetic algorithm:NSGA-II[J].IEEE Transactions onEvolutionary Computation,2002,6(2):182-19)、SMS-EMOA(Beume N,Naujoks B,Emmerich M.SMS-EMOA:Multiobjective selection based on dominated hypervolume[J].European Journal of Operational Research,2007,181(3):1653-1669)、RM-MEDA(Zhang Q,Zhou A,Jin Y.RM-MEDA:A regularity model-based multiobjectiveestimation of distribution algorithm[J].IEEE Transactions on EvolutionaryComputation,2008,12(1):41-63)和TMOEA/D(Liu H L,Gu F,Cheung Y.T-MOEA/D:MOEA/Dwith objective transform in multi-objective problems[C].Proceedings of2010International Conference on Information Science and ManagementEngineering(ISME).IEEE,2010,2:282-285)，与CEDA开展对比实验。NSGA-II是一种基于支配的MOEA，SMS-EMOA是一种基于指标的MOEA，TMOEA/D是一种用于解决具有复杂PF形状的MOPs的基于分解的MOEA，RM-MEDA是一种基于规则特性的MEDA，这几种算法涵盖了当前主流的MOEA类型。为了保证对比的公平性，所有对比算法的参数均通过前期的实验进行了系统地优化，对比实验中采用最佳的参数组合。所有算法均由Matlab进行实现，并且在同一台计算机上运行，具体的算法参数设置如下：

公共参数：

-种群大小N：在TMOEA/D中，种群的大小由权重向量的个数所决定，即(m为目标维数，D为预先设定的整数)。因此在TMOEA/D中将求解双目标(D＝65)和三目标(D＝10)MOPs的种群大小设置为N＝66。其它算法与TMOEA/D设置相同的种群大小；

变量维数：n＝10；

最大演化代数：T＝300.

NSGA-II参数：

模拟二进制交叉：P_c＝0.9,η_c＝20；

多项式变异算子控制参数：p_m＝1/n,η_m＝20.

SMS-EMOA参数：

模拟二进制交叉：P_c＝0.9,η_c＝20；

多项式变异算子控制参数：p_m＝1/n,η_m＝20.

RM-MEDA参数

聚类个数PCA：5；

局部主成分分析最大迭代次数：50；

扩展采样率：0.25.

TMOEA/D参数：

邻居大小：NS＝30；

第一搜索阶段演化代数：T1＝T/10；

第二搜索阶段演化代数：

T2＝αT，α＝{0.01,0.02,…,0.1,0.1,0.1,0.15}；

差分演化交叉算子控制参数：F＝0.5，CR＝1.

CEDA参数：

重组控制概率：β＝0.9；

最大聚类数目：K＝5；

多项式变异算子控制参数：p_m＝1/n,η_m＝20.

为了在实验中获得统计置信的结论，每种算法对每道测试题独立运行33次，基于统计指标值(均值和标准差)进行算法性能的比较。在比较表格中，关于某一道测试题，各算法对其统计运算获得的指标值的均值进行升序(IGD指标)或降序(HV指标)排序，排序结果在表格的方括号中显示，并且每种算法对GLT测试集的计算性能排序的平均值(平均秩)也列在表格中。对于每道测试题，各算法获得的平均指标值中最优值用深灰色背景表示，次优值用浅灰色背景表示。另外，当CEDA与任一算法进行比较的时候，执行在5％显著性水平的Wilcoxon秩和检验观测差异的显著性。“§”和“≈”表示CEDA求解某问题的性能在5％显著水平上是优于、劣于以及相似于比较算法对于该问题的求解能力。

实验结果

一、质量指标

表1给出了NSGA-II、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA算法分别独立计算GLT测试集33次各自获得的逼近前沿的HV和IGD值的统计结果。

从表中可以看出，通过演化300代，与对比算法进行比较，在12个指标值中，CEDA获得了8个最优和2个次优平均指标值。根据统计显著性检验，相对于NSGA-II、SMS-EMOA、RM-MEDA和TMOEA/D，在与每种算法的12次比较中，CEDA分别获得了12、11、10和7个明显较优的平均指标值。另外，平均秩的值表明在求解GLT测试集时，性能从最优到最劣的算法分别是CEDA、TMOEA/D、RM-MEDA、SMS-EMOA、NSGA-II。

二、搜索效率

图1-图6绘制出了NSGA-II、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA算法分别独立计算GLT测试集33次各自获得的平均IGD指标值的演化曲线。从图中可以看出对于GLT2-GLT3和GLT5-GLT6的求解中，CEDA均在最少的演化代数内获得了最低的平均IGD指标值。对于GLT1，CEDA的求解性能劣于RM-MEDA和TMOEA/D。对于GLT4，CEDA劣于TMOEA/D获得了次优的求解效果。总体而言，与其它四种相比，CEDA在求解GLT测试集的演化过程中收敛速度最快并且能够维持最好的种群多样性。

表1NSGA-II、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA分别独立计算GLT测试集33次所得的逼近前沿的IGD和HV指标值的统计结果(均值(标准差)[排序])

三、结果可视化

图7-图30绘制出了统计比较中性能最好的两种算法CEDA和TMOEA/D分别独立计算GLT测试集33次各自获得的全部最终逼近前沿(如图7-图18)，以及分别获得中位IGD指标值时对应的代表性逼近前沿(如图19-图30)。从图7-图18可以看出，在33次独立运算中，在求解GLT1和GLT4时，CEDA获得的逼近前沿还有部分并未收敛到PFs上，但是在求解GLT2、GLT3、GLT5、GLT6时，CEDA获得的全部逼近前沿都能稳定的收敛到PFs上，而且覆盖整个PFs。然而，TMOEA/D求解GLT4获得的逼近前沿并没有全部收敛到PF上，并且求解GLT5和GLT6时，其获得的逼近前沿并未完全覆盖整个的PFs。从图19-图30可以观察到，TMOEA/D求解GLT3和GLT4时，获得的代表性的前沿虽然最终都能收敛到PFs，但并不能完全覆盖PFs，求解GLT5和GLT6时，得到的代表性前沿中仍存在一些个体没有完全收敛到PFs，且前沿分布的均匀性也并不理想。与TMOEA/D相比，CEDA对于GLT2-GLT6获得的代表性的前沿均具有更好的收敛性和多样性。

根据上述的质量指标、搜索效率以及结果可视化，可以推断出相对于NSGA-II、SMS-EMOA、RM-MEDA和TMOEA/D，CEDA算法对于GLT测试集具有最佳的求解性能。

四、参数灵敏度分析

交配限制概率

在CEDA中，采用重组控制概率β维护算法演化过程中勘探(exploration)与开发(exploitation)之间的平衡。为了分析β对算法性能的影响，采用不同β值(β＝0.5,0.6,0.7,0.8,0.9)构造CEDA算法对GLT测试集进行求解，算法的其他参数与公共参数设置相同。每种带有不同β值的算法对每道测试题进行22次独立运算，获得的逼近前沿的IGD指标值的均值和标准差如图31所示。

由图31可以看出，当求解GLT1，GLT3和GLT4时，不同的β值获得的平均IGD值有明显的差异，然而对其它测试题求解时，不同的β值却得到了相似的平均IGD值。但是总体来说，当β＝0.9时，CEDA对于GLT1-GLT3以及GLT5-GLT6均具有较好的求解效果，因此说明算法的性能对于β的取值并不十分敏感。

聚类数目

CEDA中采用AHC方法发掘种群结构。为了分析AHC中的最大聚类数目K对CEDA性能的影响，采用不同K值(K＝4,5,7,10,20)构造CEDA算法对GLT测试集进行求解，算法中的其他参数与公共参数设置相同。每种带有不同K值的算法对每道测试题进行22次独立运算，获得的逼近前沿的IGD指标值的均值和标准差如图32所示。

由图32可以看出，当求解GLT1-GLT4时，不同K值的CEDA获得的平均IGD值有显著的差异，然而对于GLT5-GLT6求解时，不同的K值却得到了相近的平均IGD值。总体来说，当K＝5，CEDA对于不同的测试题均能获得较小的平均IGD值，因此说明CEDA的性能对于最大聚类数目的取值也并不十分敏感。

2、工程应用

优化模型

齿轮减速器是原动机和工作机之间的独立的闭式传动装置，用来降低转速和增大转矩，以满足工作需要。其无须联轴器和适配器，结构紧凑。负载分布在行星齿轮上，因而承载能力比一般斜齿轮减速机高，满足小空间高扭矩输出的需要，广泛应用于大型矿山、钢铁、化工、港口、环保等领域。虽然齿轮减速器应用广泛，但长期以来减速器的设计仅由设计人员依靠相关的资料、文献，以及多年的经验完成，因而不仅效率低，而且可能造成人力、物力和财力的浪费，因此目前需要找到一种快速有效的方法来优化设计齿轮减速器。齿轮减速器的优化设计实际上是一个多峰多目标优化问题，普通的算法难以较好地求解此问题(Farhang-Mehr A,Azarm S.Entropy-based multi-objective genetic algorithm fordesign optimization[J].Structural&Multidisciplinary Optimization,2002,24(5):351-361)本文以此问题为实例检验CEDA求解复杂工程优化问题的效果。齿轮减速器简易模型如图33所示。

该MOP的设计目标是使减速器的体积和轴2所承受的应力最小，并且满足轮齿的弯曲应力、接触应力、轴的扭转变形以及应力等约束。该问题的数学模型描述为：

s.t.:

g₆:x₁/x₂-12≤0>7:5-x₁/x₂≤0>8:1.9-x₄+1.5x₆≤0

g₉:1.9-x₅+1.1x₇≤0>stress≤1300

g_14,15:0.7≤x₂≤0.8>16,17:17≤x₃≤28>18,19:7.3≤x₄≤8.3

g_20,21:7.3≤x₅≤8.3>22,23:2.9≤x₆≤3.9>24,25:5.0≤x₁≤5.5

式中：x₁为齿宽；x₂为齿轮模数；x₃为小齿轮齿数；x₄为轴承1之间的距离；x₅为轴承2之间的距离；x₆为轴1的直径；x₇为轴2的直径；g₁为齿的弯曲应力约束；g₂为齿的接触应力约束；g₃、g₄为轴的变形约束；g₅、g₆、g₇为基于空间尺寸限制和经验约束；g₈、g₉为由经验得到的设计轴的要求；g₁₀、g₁₁为轴的应力约束；g₁₂至g₂₅为7个变量的上下边界。

对于步骤一中P＝{x¹,x²,…,x^N}，其中：

x¹＝{x₁,x₂,x₃,x₄,x₅,x₆,x₇}

x²＝{x₁,x₂,x₃,x₄,x₅,x₆,x₇}

x³＝{x₁,x₂,x₃,x₄,x₅,x₆,x₇}

.............................................

x^N＝{x₁,x₂,x₃,x₄,x₅,x₆,x₇}

x¹…x^N都代表x₁,x₂,x₃,x₄,x₅,x₆,x₇，但取值不同。

实验设计与结果分析

利用NSGA-II、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA对齿轮减速器优化设计模型进行求解。经过参数优化，运算过程中的参数设置如表2所示，其余设计与实施例1中参数设置相同。每种算法对模型独立运算33遍，利用超体积HV指标值度量每一次获得的逼近前沿的效果。其中，计算HV值时取参考点r＝[6600,1600]^T。

表2NSGA-II、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA算法求解齿轮减速器优化设计模型时参数设置

五种算法对齿轮减速器优化设计模型独立运行33次获得的HV指标值的箱型图对比结果如图34和图35(34图为原图，35图为局部放大图)。从图中可以看出CEDA获得了最大的中位HV指标值和最小的四分位距，从而说明CEDA对于齿轮减速器优化设计模型能稳定地求解出具有良好多样性和收敛性的解。

图36绘制出了NSGA-II、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA算法分别独立计算齿轮减速器优化设计模型33次各自获得的平均HV指标值的演化曲线。从图中可以看出CEDA在最小的演化代数内获得了最高的平均HV指标值。也就是说与其它四种相比，CEDA在演化过程中收敛速度最快并且能够维持最好的种群多样性。

图37-图40为分别利用NSGA-II和CEDA对齿轮减速器优化设计模型求解时，独立运算33次各自获得的全部逼近前沿以及中位IGD指标值对应的代表性逼近前沿。从图38可以看出，CEDA获得的全部逼近前沿均能够稳定地收敛，并且与NSGA-II获得的逼近前沿相比，其覆盖面更广。从图39可以看出，相对于NSGA-II，CEDA获得了更为宽广和均匀的代表性的逼近前沿。从对图34-图40的分析中可以推断出CEDA算法对于齿轮减速器优化设计模型具有优异的求解性能。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。