CN103488851A - 基于几何结构信息的多目标优化方法 - Google Patents

Abstract

本发明公开了基于几何结构信息的多目标优化方法，首先初始化种群，判断优化问题的优化目标函数个数，将多目标优化函数进行降维处理，降维到2维或3维，然后对当前种群的前沿解集进行参数化，并将其拟合为非均匀有利B样条方程，将前沿解集由优化目标函数空间映射到非均匀有理B样条参数空间，根据拟合的非均匀有理B样条方程求解牵引点集，并利用牵引点集对所有被支配解进行优化，使当前种群朝着优化问题的真实前沿逼近。本发明的方法在处理多目标优化问题时能更快地收敛于真实的Pareto前沿，解集中包含的非支配解的个数和分布均匀性优于传统方法，且对于优化问题本身具有很好的鲁棒性。

Description

基于几何结构信息的多目标优化方法

技术领域

本发明涉及工程设计领域，具体涉及对于多域的复杂工程产品，多领域的设计优化问题，以及基于启发式方法原理的多目标优化方法。

背景技术

在处理工程领域问题时，通常需要同时优化多个目标，意大利经济学家帕累托(Pareto)提出“Pareto最优”的思想常被用来处理多目标优化问题。目前，研究人员对多目标优化问题已经进行了大量研究，这些工作大致可以分为两类：第一类是将多个目标函数聚合为一个目标或者将多余的目标函数形式化为约束条件只保留一个目标函数，即可简化单目标的优化问题，因而可以利用丰富的单目标优化方法；第二类方法是确定一个Pareto前沿解集或者它的一个非常有代表性的子集合，集合中的解是互不支配的。针对多目标优化问题，也已经有一些智能方法，其中，遗传方法和微粒群是研究最多、应用最广的两种启发性方法。非均匀有理B样条NURBS(Non-uniform Rational B-Splines)是专门做曲线曲面物体的一种造型方法，对于大量数据散乱点用NUBRS能够拟合出其几何造型。

为了得到一个多样性的解集和提升优化方法的收敛速度，传统的演化方法在求解多目标优化问题时利用一些特殊的适应度函数和选择策略。现有的很多遗传方法的变种通过构造特定的适应度函数对种群进行排序，根据粒子的适应度值选择父代粒子生成后代，通过这样的策略保持优秀基因的延续性。微粒群是用于求解多目标优化问题的一种经典的生物智能启发式方法，该方法能得到非常丰富的Pareto前沿解集，因为微粒群方法求解问题的广泛性和有效性使得它越来越受人们的关注，为此研究人员已经提出超过30多个类似的方法。其变种则是通过利用粒子自己的局部最优记录和整个种群的全局最优记录来更新粒子的飞行速度，使得粒子朝着性能更优的方向飞行。虽然这些启发式方法在处理多目标优化问题时取得了一定的进展，但仍存在以下三个问题：第一，没有强有力的方法保证解集分布的均匀性；第二，粒子可行域和最优飞行方向的不可预见性，导致其收敛速度的不稳定；第三，不能逼近在不可行域外但感兴趣的解，通常在具体的设计问题中，工程师经常会探索为了得到某个极优的性能值需要在其它方面付出更多的代价。

发明内容

本发明的目的就是针对现有技术中存在的上述问题，提出一种基于几何结构信息的多目标优化方法，利用启发式方法中种群的几何结构指导种群的演化。

一种基于几何结构信息的多目标优化方法，包括步骤：

步骤1、初始化种群，设置进化代数计数器t=0，设置最大进化代数T，随机生成M个个体作为初始种群；

步骤2、判断优化问题的优化目标函数个数，如果大于3则进行降维处理后进入下一步骤，否则直接进入下一步骤；

步骤3、抽取当前种群的前沿解集，并对其进行参数化，根据参数化的结果，拟合其非均匀有理B样条方程；

步骤4、根据拟合的非均匀有理B样条方程求非支配解的单位法线向量，以非支配解的范数作为非支配解的前进步长，以非支配解的坐标与其前进步长之和得到的新点称之为非支配解的牵引点，求解牵引点集；

步骤5、利用牵引点集对所有被支配解进行优化；

步骤6、对前沿解集作均匀化处理或向外延拓处理；

步骤7、令t=t+1，判断t是否小于T，如果小于T则返回步骤3继续进行迭代，否则返回最终前沿解集。

进一步地，所述降维处理，包括如下步骤：

步骤2.1、对于有d个优化目标函数的多目标优化问题，随机产生d个测试种群并求解其优化目标函数值，每个测试种群中个体数量为n，计算这d个测试种群优化目标函数值的平均值作为最后的测试种群；

步骤2.2、将最后测试种群的优化目标函数空间看成是n*d的矩阵，将每个优化目标函数看作一个目标点族，求解该矩阵的相关系数矩阵d*d；

步骤2.3、根据相关系数矩阵计算所有目标点族中任意两者之间的目标点族距离；

步骤2.4、根据目标点族距离，选择目标点族距离最大的两个目标点族，将其合并为一个新的目标点族，此时目标点族数量减1；

步骤2.5、判断目标点族数量是否等于2或3，如果否，则返回步骤2.3，如果是则对于每个包含多个优化目标函数的目标点族，将多优化目标函数聚合为单目标函数。

进一步地，所述步骤2.5还包括步骤：

步骤2.5.1、若当前目标点族中存在两个优化目标函数之间相关系数小于给定的阈值的模糊依赖关系，则分别计算两个优化目标函数与其他优化目标函数的距离值，选择值小的那个目标优化函数进行调整，将其从当前目标点族中删除；若当前目标点族中存在一个优化目标函数与所有其他的优化目标函数都负相关的模糊依赖关系，则对当前目标点族中与所有其他的优化目标函数都负相关的优化目标函数进行调整，将其从当前目标点族中删除；

步骤2.5.2、判断是否迭代完毕，如果迭代完毕则停止调整，否则依次迭代其他目标点族，将迭代的目标点族设为当前目标点族，将需要调整的优化目标函数加入，若当前目标点族中仍然存模糊依赖关系，则返回步骤步骤2.5.1，否则停止调整；

步骤2.5.3、对于一个需要调整的优化目标函数，如果没有一个目标点族能够在放入该优化目标函数后，不存在模糊依赖关系，则不进行调整，将该优化目标函数放入原来的目标点族。

进一步地，所述将多优化目标函数聚合为单目标函数，包括步骤：

在测试种群中计算此目标点族中每个优化目标函数的数量级；

对每个优化目标函数除以其对应的数量级进行归一化；

目标点族最终的优化目标函数为其所包含的所有优化目标函数归一化后之和;

对于标记的优化目标函数，取其相反数相加，从而减少它对最后的优化目标函数的影响。

进一步地，所述根据相关系数矩阵计算所有目标点族中任意两者之间的目标点族距离，若设有优化目标点族G1，含有a个优化目标点，优化目标点族G2，含有b个优化目标点，r_kj表示G1中第k个优化目标点与G2中第j个优化目标点之间的相关系数，即在相关系数矩阵中第k行j列的数值，则定义目标点族距离cd为：

$\mathrm{cd}=\frac{\underset{k=1}{\overset{a}{\mathrm{\Σ}}}\underset{j=1}{\overset{b}{\mathrm{\Σ}}}{r}_{\mathrm{kj}}}{\mathrm{ab}}\·$

进一步地，所述抽取当前种群的前沿解集，并对其进行参数化，若优化目标函数空间为二维，则其前沿解集上每个粒子的详细参数化过程包括步骤：

对前沿解集在优化目标函数空间根据第一维进行升序排序；

求出前沿解集中相邻粒子在优化目标函数空间的欧氏距离，计算所有欧氏距离值之和得到总长度；

对于某个粒子，计算在其前面粒子的欧氏距离值之和，将得到的欧氏距离值除以总长度即为该粒子对应的参数值u。

进一步地，所述抽取当前种群的前沿解集，并对其进行参数化，若优化目标函数空间为三维，则其前沿解集上每个粒子的详细参数化过程包括步骤：

找到两个正整数(p,q)，使得p*q的值为前沿解集的大小，将一维的前沿解集向量看成是二维的p*q的矩阵，p为参数u方向上的个数，q为参数v方向上的个数；

使用二维优化目标函数空间的方法对u参数方向进行参数化；

在v方向上，即二维矩阵p*q的每一列q个粒子根据粒子在优化目标函数空间第二维进行升序排序，使用二维优化目标函数空间的方法对v参数方向进行参数化；

计算每一列u的平均值，该列所有粒子的u值即为此平均值；

计算每一行v的平均值，该行所有粒子的v值即为此平均值。

进一步地，所述根据参数化的结果，拟合其非均匀有理B样条方程，如果优化目标函数个数为2，则拟合成非均匀有理B样条曲线，优化目标函数个数为3则拟合为非均匀有理B样条曲面。

进一步地，所述利用牵引点集对所有被支配解进行优化，包括步骤：

步骤5.1、初始化被支配解集中当前被支配解的速度为0；

步骤5.2、计算当前被支配解与前沿解集中所有非支配解在优化目标函数空间的欧氏距离，选择两个欧氏距离最近的非支配解，找到这两个非支配解对应的牵引点，计算两个牵引点与其对应的非支配解之间的欧氏距离和d1；

步骤5.3、利用两个欧氏距离最近的非支配解，按照速度更新公式改变当前被支配解的飞行速度，按照位置更新公式改变其优化变量的值，最后重新求解其对应的优化目标函数值；

二维速度更新公式为：

$\stackrel{\→}{v(i+1)}=\mathrm{\ω}\stackrel{\→}{v\left(i\right)}+C_1(\stackrel{\→}{x_{p_1}}-\stackrel{\→}{x})+C_2(\stackrel{\→}{x_{p_2}}-\stackrel{\→}{x});$

其中，

是当前被支配解在第i次迭代的速度，

是当前被支配解的当前位置，

和

是它的两个最近的非支配解的位置，ω是被支配解速度的权重，C₁和C₂是分别受两个非支配解影响的权重，是当前被支配解在第（i+1）次迭代的速度；

三维速度更新公式为：

$\stackrel{\→}{v(i+1)}=\mathrm{\ω}\stackrel{\→}{v\left(i\right)}+C_1(\stackrel{\→}{x_{p_1}}-\stackrel{\→}{x})+C_2(\stackrel{\→}{x_{p_2}}-\stackrel{\→}{x})+C_3(\stackrel{\→}{x_{p_3}}-\stackrel{\→}{x});$

其中，和

是当前被支配解的三个最近的非支配解的位置，C₁、C₂和C₃是分别受三个非支配解影响的权重；

二维和三维位置更新公式为：

$\stackrel{\→}{x(i+1)}=\stackrel{\→}{x\left(i\right)}+\stackrel{\→}{v(i+1)};$

其中，

是当前被支配解在第i次迭代时的位置，

是当前被支配解在(i+1)次迭代时的速度，是当前被支配解在第(i+1)次迭代时的更新位置；

步骤5.4、计算当前被支配解与两个牵引点之间的欧氏距离之和d2；

步骤5.5、如果d1小于d2，说明被支配解未被牵引到更优区域内，此时返回步骤5.3，若d1不小于d2，则被支配解被牵引到更优区域，优化过程结束。

进一步地，所述对前沿解集作均匀化处理或向外延拓处理，包括步骤：

根据前沿解集中非支配解的参数值计算相邻非支配解间的距离值，找到距离值最大的两个非支配解，表明在整个前沿解集中它们之间最稀疏；

取这两个非支配解参数值的平均值作为即将插入点的参数值，根据非均匀有理B样条方程求出这个参数值对应的数据点，即在优化目标函数空间的目标函数值；

利用信赖域优化方法求解优化变量的值来逼近期望的优化目标值。

本发明提出的基于几何结构信息的多目标优化方法，首先将多目标优化函数进行降维处理，降维到2维或3维，然后对当前种群的前沿解集进行参数化，并将其拟合为非均匀有利B样条方程，将前沿解集由优化目标函数空间映射到非均匀有理B样条参数空间，根据拟合的非均匀有理B样条方程求解牵引点集，并利用牵引点集对所有被支配解进行优化，使当前种群朝着优化问题的真实前沿逼近。本发明提出的方法在处理多目标优化问题时相比较于其他的经典方法，比如NSGAII，更快地收敛于真实的Pareto前沿，解集中包含的非支配解的个数和分布均匀性优于其他方法，且对于优化问题本身具有很好的鲁棒性。

附图说明

图1为本发明基于几何结构信息的多目标优化方法流程图；

图2为本发明对多优化目标进行降维处理方法流程；

图3为本发明模糊依赖关系处理方法流程图；

图4为本发明被支配解优化方法流程；

图5a为本发明模糊依赖关系示意图一；

图5b为本发明模糊依赖关系示意图二；

图5c为本发明模糊依赖关系示意图三；

图5d为本发明模糊依赖关系示意图四；

图6为本发明被支配解牵引示意图。

具体实施方式

下面结合附图和实施例对本发明技术方案做进一步详细说明，以下实施例不构成对本发明的限定。

本发明利用启发式方法中种群的几何结构指导种群的演化，其主要思想是将当前种群前沿解集当作几何点集来处理，通过拟合求得其几何结构信息，然后利用求得的几何结构信息来指导粒子的飞行。在本实施例中，将当前种群前沿解集当作几何点集来处理，下文中将所有前沿解集中的非支配解和非前沿解集中的被支配解，都通称为粒子。

如图1所示，本发明基于几何结构信息的多目标优化方法包括步骤：

步骤S1、初始化种群，设置进化代数计数器t=0，设置最大进化代数T，随机生成M个个体作为初始种群P(0)。

步骤S2、判断优化问题的优化目标函数个数，如果大于3则进行降维处理后进入下一步骤，否则直接进入下一步骤。

对于多目标优化问题，通常需要将多目标降维到二维或三维来进行处理，对多目标进行降维处理方法流程如图2所示，包括如下步骤：

步骤201、对于有d个优化目标函数的多目标优化问题，随机产生d个测试种群并求解其优化目标函数值，每个测试种群中个体数量为n，计算这d个测试种群优化目标函数值的平均值作为最后的测试种群。

对于待求解的多目标优化问题，通过随机产生的与初始种群独立的测试种群，假设测试种群中个体数量为n，优化目标函数个数为d，随机产生d个测试种群并求解其优化目标函数值，为了获得测试种群优化目标函数值的分布一般情况，计算这d个测试种群优化目标函数值的平均值作为最后的测试种群。

步骤202、将最后测试种群的优化目标函数空间看成是n*d的矩阵，将每个优化目标函数看作一个目标点族，求解该矩阵的相关系数矩阵d*d。

矩阵的每一行对应于一个个体的优化目标函数值，每一列对应于一个优化目标函数，通过相关系数矩阵知道任意两个优化目标函数之间的距离，且相关系数值越大表示距离越近，通过相关系数矩阵得到每个优化目标函数之间的相关度。

步骤203、根据相关系数矩阵计算所有目标点族中任意两者之间的目标点族距离。

为了便于说明，将最后测试种群的优化目标函数空间n*d矩阵中每个优化目标函数认为是一个优化目标点，本发明将由一个或多个优化目标点形成的集合称为目标点族，设有优化目标点族G1，含有a个优化目标点，优化目标点族G2，含有b个优化目标点，r_kj表示G1中第k个优化目标点与G2中第j个优化目标点之间的相关系数，即在相关系数矩阵中第k行j列的数值，则定义目标点族距离cd为：

$\mathrm{cd}=\frac{\underset{k=1}{\overset{a}{\mathrm{\Σ}}}\underset{j=1}{\overset{b}{\mathrm{\Σ}}}{r}_{\mathrm{kj}}}{\mathrm{ab}}\·$

可见通过求解的相关系数矩阵可以计算目标点族距离。

步骤204、根据目标点族距离，选择目标点族距离最大的两个目标点族，将其合并为一个新的目标点族，此时目标点族数量减1。

本发明基于求得的相关系数矩阵进行多目标降维，为了将多个目标函数映射为两个或三个目标函数，首先在初始时将每个优化目标点可以认为是一个目标点族，计算所有目标点族中任意两者之间的目标点族距离，选择目标点族距离最大的两个目标点族，将其合并为一个新的目标点族，此时目标点族数量减1。

步骤205、判断目标点族数量是否等于2或3，如果否，则返回步骤203继续迭代，如果是则对于每个包含多个优化目标函数的目标点族，将多优化目标函数聚合为单目标函数。

重复步骤203和204，直到目标点族数量等于2或3，完成多目标降维。

在使用上述方法进行多目标降维后，每个目标点族包含一个至多个优化目标函数，目标点族中的多个优化目标函数互相之间的依赖关系强弱不定，这种模糊的依赖关系在一个目标点族中经常会发生。降维后每个目标点族理想的状态是其中的任何两个优化目标函数之间都是正相关，如图5a所示，点A1，B1和C1之间都是正相关。但通常情况是一个目标点族中经常会包含负相关的情形，为此，设定一个阈值σ=-0.2，两个优化目标函数之间如果相关系数在-0.2与0之间，则认为是弱相关，当两个优化目标函数之间的相关系数值小于0但是大于阈值时，认为它虽然不理想但属于可以接受的状态，如图5b中B2和C2之间相关系数为-0.18，认为C2为可接受情形；反之如果其值小于-0.2，如图5c中所示，B3和C3之间相关系数为-0.84，此时需要进行调整。

本发明在降维过程中，即在步骤205判断判断目标点族数量是否等于2或3之前，对带有模糊依赖关系的目标点族，进行调整。这里的模糊关系，是指若当前目标点族中存在两个优化目标函数之间相关系数小于给定的阈值，或当前目标点族中存在一个优化目标函数与所有其他的优化目标函数都负相关，具体调整策略流程如图3所示，包括如下步骤：

步骤301、若当前目标点族中存在两个优化目标函数之间相关系数小于给定的阈值的模糊依赖关系，则分别计算两个优化目标函数与其他优化目标函数的距离值，选择值小的那个目标优化函数进行调整，将其从当前目标点族中删除；若当前目标点族中存在一个优化目标函数与所有其他的优化目标函数都负相关的模糊依赖关系，则对当前目标点族中与所有其他的优化目标函数都负相关的优化目标函数进行调整，将其从当前目标点族中删除。

如遍历图5a、图5b、图5c、图5d四个目标点族，发现图5c所示的目标点族中B3和C3间相关系数小于给定的阈值-0.2，则将图5c所示的目标点族设为当前目标点族。

或者，发现图5d所示目标点族中，有个优化目标函数B4虽然与其他优化目标函数的相关系数在阈值内，但是它与目标点族内所有其他的优化目标函数都负相关，表明它的存在虽然对某个目标的影响有限，但是对于整个目标点族影响很大，为此也需要调整，B4与A4和C4都负相关，则将图5d所示的目标点族设为当前目标点族，B4需要进行调整。

如图5c所示，对不符合要求的两个优化目标函数B3和C3，分别计算两个优化目标函数与其他优化目标函数的距离值，分别为0.58和0.36，选择值小的那个目标优化点进行调整，将其从当前目标点族中删除，如图5c中删除点C3。

如图5d所示，B4与A4和C4都负相关，为此B4需要调整，删除B4。

步骤302、判断是否迭代完毕，如果迭代完毕则停止调整，否则依次迭代其他目标点族，将迭代的目标点族设为当前目标点族，将需要调整的优化目标函数加入，若当前目标点族中仍然存模糊依赖关系，则返回步骤301，否则停止调整。

例如，在步骤302中删除的C3，依次迭代其他三个目标点族，如迭代如图5a所示的目标点族，将C3加入到图5a所示的当前目标点族，判断是否存在模糊依赖关系，即判断加入了C3的当前目标点族中是否存在两个优化目标函数之间相关系数小于给定的阈值，或当前目标点族中存在一个优化目标函数与所有其他的优化目标函数都负相关的现象，如果有则根据步骤301-步骤302，继续进行调整。

步骤303、对于一个需要调整的优化目标函数，如果没有一个目标点族能够在放入该优化目标函数后，不存在模糊依赖关系，则不进行调整，将该优化目标函数放入原来的目标点族。

如果不存在模糊依赖关系了，即调整后的目标点族都如图5a或图5b所示，则停止调整。在停止调整时，对于一个需要调整的优化目标函数，如果没有一个目标点族能达到调整要求，则不强行调整，因为已经得到该优化目标函数与原目标点族距离最大，在此目标点族中标记此需要调整的优化目标函数，进行合成时需要处理；如果找到更适合的目标点族，则将该需要调整的优化目标函数放入新的目标点族。

对于每个包含多个优化目标函数的目标点族，将多优化目标函数聚合为单目标函数，具体聚合过程为：1)在测试种群中计算此目标点族中每个优化目标函数的数量级；2)对每个优化目标函数除以其对应的数量级进行归一化;3)目标点族最终的优化目标函数为其所包含的所有优化目标函数归一化后之和;4)对于标记的优化目标函数，取其相反数相加，从而减少它对最后的优化目标函数的影响。

步骤S3、抽取当前种群的前沿解集，并对其进行参数化，根据参数化的结果，拟合其非均匀有理B样条方程。

种群在每一代演化时，计算其Pareto最优前沿，即前沿解集，前沿解集中的解在优化目标函数上是互不支配的，又称为非支配解。参数化前沿解集的目的是将前沿解集由优化目标函数空间映射到非均匀有理B样条曲线（Non-Uniform Rational B-Splines）NURBS参数空间，求解出其对应的参数值，即将由前沿解集形成的散乱点轻松地映射到几何空间。在几何空间中，每个散乱点有一个参数u值，散乱点集合的参数u值集合称为对应的NURBS参数u序列，如果优化目标函数空间是三维，每个散乱点有一个参数u值和一个参数v值，散乱点集合的参数(u,v)值集合称为对应的NURBS参数(u,v)值序列。如果优化目标函数空间是二维则将其映射到NURBS曲线空间，求解出其对应的NURBS参数u值序列，如果优化目标函数空间是三维则将其映射到NURBS曲面空间，求解出其对应的NURBS参数(u,v)值序列。

具体地，优化目标函数空间为二维则其前沿解集上每个粒子的详细参数化过程如下：

对前沿解集在优化目标函数空间根据第一维进行升序排序；

求出相邻粒子在优化目标函数空间的欧氏距离，计算所有距离值之和得到总长度；

对于某个粒子，计算在其前面粒子的距离值之和，将得到的距离值除以总长度即为该粒子对应的参数值u。

具体地，优化目标空间为三维则其前沿解集上每个粒子的参数化过程如下：

找到两个正整数(p,q)，使得p*q的值为前沿解集的大小，将一维的前沿解集向量看成是二维的p*q的矩阵，p为参数u方向上的个数，q为参数v方向上的个数；

使用前面介绍的二维优化目标函数空间的方法对u参数方向进行参数化；

在v方向上，即二维矩阵p*q的每一列q个粒子根据粒子在优化目标函数空间第二维进行升序排序，使用前面介绍的二维优化目标函数空间的方法对v参数方向进行参数化；

计算每一列u的平均值，该列所有粒子的u值即为此平均值；

计算每一行v的平均值，该行所有粒子的v值即为此平均值。

根据前面的参数化的结果，拟合NURBS方程，如果优化目标函数个数为2，则拟合成NURBS曲线，优化目标函数个数为3则拟合为NURBS曲面。将前沿解集映射到几何参数空间后，基于NURBS曲线和曲面模型对前沿解集进行拟合，其中包括确定NURBS模型的节点向量和求解控制顶点坐标。

节点向量是关于参数一个数值集合，如果是NURBS曲面的话，则关于参数u和参数v分别有一个节点向量，其大小等于该参数方向上控制点的个数与基函数阶数之和，基函数次数选取3，所以其阶数为4。本文中采用累积弦长参数化构造节点向量，对于NURBS曲线其具体过程如下：

在节点向量中最前面阶数个元素赋值为0，最后面阶数个元素赋值为1。

抽取所有粒子的参数值u，构造成一个结构化的参数序列。

针对节点向量中剩余未知的第i个节点，在参数序列中以第i个参数为起始，选择连续阶数加1个参数，将其平均值赋值给第i个节点。

对于NURBS曲面，参数u方向和参数v方向上的向量选取与NURBS曲线相同。

前沿解集中的散乱点是拟合为NURBS结构的数据点，NUBRS模型中控制点和数据点的个数相同，在本实施例中默认控制点的权重为1，对于NURBS曲线求解控制点具体过程如下：

输入控制顶点个数，基函数的次数，节点向量，权重矩阵，数据点的参数值，求解基函数的值，构造线性方程组的系数矩阵；

求解系数矩阵的逆矩阵；

针对数据点的每一维乘以上步得到的逆矩阵，得到对应控制点那一维度的值。

对于NUBRS曲面求解控制点过程，在第一步中要分别输入参数u方向和参数v方向的节点向量。

得到控制点以后，NURBS模型中所有涉及到的元素都已确定，即确定了对于前沿解集的散乱点拟合出了其对应的NURBS方程。

步骤S4、根据拟合的NUBRS方程求非支配解的单位法线向量，以非支配解的范数作为非支配解的前进步长，以非支配解的坐标与其前进步长之和得到的新点称之为非支配解的牵引点，求解牵引点集。

本实施例主要用到前沿解集几何结构中粒子的法向，利用其代表的方向表示解更优的方向指导粒子飞行。前沿解集中每个粒子对应一个牵引点，故牵引点集大小等于前沿点集，其求解过程如下：

针对当前粒子，根据其参数值求出单位切向量，如果是NURBS曲面则分别求出参数u方向和参数v方向上的单位偏导向量且此时偏导向量是三维的；如果是NURBS曲线，利用法线向量和切向量相互垂直的关系，将单位切向量旋转90度得到对应的单位法向量。如果是NURBS曲面，两个单位偏导向量的外积是对应的单位法向量。

求解当前粒子优化目标函数空间的范数，将其乘以法向量得到该粒子的前进步长，当前粒子的坐标与其前进步长之和即为该粒子的牵引点。

分析比较当前粒子与其牵引点的支配关系，如果牵引点被当前粒子支配，则说明法线方向取反了，此时将步长置反，重新计算牵引点的值。

对于每个粒子，按照如上步骤计算其对应的牵引点，构成了最后的牵引点集。

牵引点的高效计算主要得益于NURBS曲线和NURBS曲面能够方便地求解在特定参数位置的切向量。

步骤S5、利用牵引点集对所有被支配解进行优化。

本实施例在每一次迭代对于处于前沿解集上的粒子不进行更新，只更新处于不在前沿解集上的粒子的位置，即调整被支配解集的位置，指导的基本思想是利用牵引点在优化目标函数空间位置的更优性牵引被支配粒子朝其前进。对于每个粒子，牵引的过程可以认为是一个局部优化的过程，通过不断地改变其速度来改变其位置，能够在很少的次数内使得该被支配的粒子比在前沿的粒子更优，且与牵引点越来越接近。对于任意一个被支配解，粒子的局部优化具体过程如图4所示，包括步骤：

步骤401、初始化被支配解集中当前被支配解的速度为0。

步骤402、计算当前被支配解与前沿解集中所有非支配解在优化目标函数空间的欧氏距离，选择两个欧氏距离最近的非支配解，找到这两个非支配解对应的牵引点，计算两个牵引点与其对应的非支配解之间的欧氏距离和d1。

如图6中，坐标轴分别表示优化目标函数f1和f2，A、B、C、D和E为前沿解集上的点，对于被支配解P，前沿解A和C作为离其最近的两个非支配解，A和C与其对应的牵引点Ag、Cg构成的四边形更优区域为被支配解P的优化目标区域。

步骤403、利用两个欧氏距离最近的非支配解，按照速度更新公式改变当前被支配解的飞行速度，按照位置更新公式改变其优化变量的值，最后重新求解其对应的优化目标函数值。

二维粒子速度更新公式为：

$\stackrel{\→}{v(i+1)}=\mathrm{\ω}\stackrel{\→}{v\left(i\right)}+C_1(\stackrel{\→}{x_{p_1}}-\stackrel{\→}{x})+C_2(\stackrel{\→}{x_{p_2}}-\stackrel{\→}{x})$

上述公式中，

是当前粒子在第i次迭代的速度，

是当前粒子的当前位置，

和

是它的两个最近的非支配解的位置，ω是粒子速度的权重，取值0.1，C₁和C₂是分别受两个非支配解影响的权重，取值0.45，

是当前粒子在第（i+1）次迭代的速度；

三维粒子速度更新公式为：

$\stackrel{\→}{v(i+1)}=\mathrm{\ω}\stackrel{\→}{v\left(i\right)}+C_1(\stackrel{\→}{x_{p_1}}-\stackrel{\→}{x})+C_2(\stackrel{\→}{x_{p_2}}-\stackrel{\→}{x})+C_3(\stackrel{\→}{x_{p_3}}-\stackrel{\→}{x})$

上述公式中，和

是当前被支配粒子的三个最近的非支配解的位置，C₁、C₂和C₃是分别受三个非支配解影响的权重，取值0.3；

二维和三维粒子位置更新公式为：

$\stackrel{\→}{x(i+1)}=\stackrel{\→}{x\left(i\right)}+\stackrel{\→}{v(i+1)}$

上述公式中，

是当前被支配粒子在第i次迭代时的位置，

是当前被支配粒子在(i+1)次迭代时的速度，是当前被支配粒子在第(i+1)次迭代时的更新位置。

步骤404、计算当前被支配解与两个牵引点之间的欧氏距离之和d2。

步骤405、如果d1小于d2，说明被支配解未被牵引到更优区域内，此时返回步骤403，若d1不小于d2，则被支配解被牵引到更优区域，优化过程结束。

如图6中点P在局部优化第i代的位置P_i，经比较发现此时P_i点未在更优区域内，在第(i+1)代中点P_i+1位于更优区域内，被牵引到理想位置，此时局部优化过程结束。

选择两个最近的前沿解能够保证被支配解最终会被牵引到大致它们中间的位置，一定程度上保证了种群的均匀性。当所有被支配粒子都执行局部优化过程后，种群整体上往更优的区域内迁移了一步，有些被支配解被优化为非支配解，此时Pareto前沿解集进行了更新，如此迭代下去，则整个种群朝着优化问题的真实Pareto前沿逼近。

步骤S6、对前沿解集作均匀化处理或向外延拓处理。

前沿解集的均匀化程度越高，有意义的解数量越多，对于具体的工程设计问题设计人员在进行决策时有更多的方案可以选择。本发明通过在前沿解集上插入点使其内部均匀和在其外部插入点使其延拓的方法，均匀化前沿解集。传统的多目标优化方法是计算通过目标空间的欧氏距离分析前沿解集的稀疏性质，且很难保证前沿解集具有很好的均匀性，而本发明提出的方法则通过分析前沿解集的参数化信息，能方便快速地得到其稀疏的分布情况，同时在计算距离时考虑的维度更少，通过直观的NUBRS几何参数比较，最后得到可信度更高的稀疏分布。如果优化目标函数空间是二维，在NURBS几何空间分析的时候只需要比较一维的u参数值；如果优化目标函数空间是三维，则在NUBRS几何空间分析的时候只需要比较两维的(u,v)参数值。另外，更为重要的是：针对当前前沿解集，设计人员可基于本方法，根据当前实际需要决定需要插入多少个粒子来均匀化前沿解集。

对于优化目标函数空间为二维，插入一个粒子其具体过程如下：

根据粒子的u参数值计算相邻粒子间的距离值，找到距离值最大的两个粒子，表明在整个前沿解集中它们之间最稀疏。

取这两个粒子u值的平均值作为即将插入的粒子u值，根据NURBS曲线方程求出这个u值对应的数据点，即在优化目标函数空间的两个目标函数值。

利用信赖域优化方法求解优化变量的值来逼近期望的优化目标值。

对于优化目标函数空间为三维的情况，与二维不同的是要考虑到粒子u和v两个参数。

步骤S7、令t=t+1，判断t是否小于T，如果小于T则返回步骤S3继续进行迭代，否则返回最终前沿解集。

以上实施例仅用以说明本发明的技术方案而非对其进行限制，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (10)

1.一种基于几何结构信息的多目标优化方法，应用于多目标优化问题，其特征在于，包括步骤：

步骤1、初始化种群，设置进化代数计数器t=0，设置最大进化代数T，随机生成M个个体作为初始种群；

步骤2、判断优化问题的优化目标函数个数，如果大于3则进行降维处理后进入下一步骤，否则直接进入下一步骤；

步骤3、抽取当前种群的前沿解集，并对其进行参数化，根据参数化的结果，拟合其非均匀有理B样条方程；

步骤4、根据拟合的非均匀有理B样条方程求非支配解的单位法线向量，以非支配解的范数作为非支配解的前进步长，以非支配解的坐标与其前进步长之和得到的新点称之为非支配解的牵引点，求解牵引点集；

步骤5、利用牵引点集对所有被支配解进行优化；

步骤6、对前沿解集作均匀化处理或向外延拓处理；

步骤7、令t=t+1，判断t是否小于T，如果小于T则返回步骤3继续进行迭代，否则返回最终前沿解集。

2.根据权利要求1所述的基于几何结构信息的多目标优化方法，其特征在于，所述降维处理，包括如下步骤：

步骤2.1、对于有d个优化目标函数的多目标优化问题，随机产生d个测试种群并求解其优化目标函数值，每个测试种群中个体数量为n，计算这d个测试种群优化目标函数值的平均值作为最后的测试种群；

步骤2.2、将最后测试种群的优化目标函数空间看成是n*d的矩阵，将每个优化目标函数看作一个目标点族，求解该矩阵的相关系数矩阵d*d；

步骤2.3、根据相关系数矩阵计算所有目标点族中任意两者之间的目标点族距离；

步骤2.4、根据目标点族距离，选择目标点族距离最大的两个目标点族，将其合并为一个新的目标点族，此时目标点族数量减1；

步骤2.5、判断目标点族数量是否等于2或3，如果否，则返回步骤2.3，如果是则对于每个包含多个优化目标函数的目标点族，将多优化目标函数聚合为单目标函数。

3.根据权利要求2所述的基于几何结构信息的多目标优化方法，其特征在于，所述步骤2.5还包括步骤：

步骤2.5.1、若当前目标点族中存在两个优化目标函数之间相关系数小于给定的阈值的模糊依赖关系，则分别计算两个优化目标函数与其他优化目标函数的距离值，选择值小的那个目标优化函数进行调整，将其从当前目标点族中删除；若当前目标点族中存在一个优化目标函数与所有其他的优化目标函数都负相关的模糊依赖关系，则对当前目标点族中与所有其他的优化目标函数都负相关的优化目标函数进行调整，将其从当前目标点族中删除；

步骤2.5.2、判断是否迭代完毕，如果迭代完毕则停止调整，否则依次迭代其他目标点族，将迭代的目标点族设为当前目标点族，将需要调整的优化目标函数加入，若当前目标点族中仍然存模糊依赖关系，则返回步骤步骤2.5.1，否则停止调整；

步骤2.5.3、对于一个需要调整的优化目标函数，如果没有一个目标点族能够在放入该优化目标函数后，不存在模糊依赖关系，则不进行调整，将该优化目标函数放入原来的目标点族。

4.根据权利要求2所述的基于几何结构信息的多目标优化方法，其特征在于，所述将多优化目标函数聚合为单目标函数，包括步骤：

在测试种群中计算此目标点族中每个优化目标函数的数量级；

对每个优化目标函数除以其对应的数量级进行归一化；

目标点族最终的优化目标函数为其所包含的所有优化目标函数归一化后之和;

对于标记的优化目标函数，取其相反数相加，从而减少它对最后的优化目标函数的影响。

5.根据权利要求2所述的基于几何结构信息的多目标优化方法，其特征在于，所述根据相关系数矩阵计算所有目标点族中任意两者之间的目标点族距离，若设有优化目标点族G1，含有a个优化目标点，优化目标点族G2，含有b个优化目标点，rkj表示G1中第k个优化目标点与G2中第j个优化目标点之间的相关系数，即在相关系数矩阵中第k行j列的数值，则定义目标点族距离cd为：

$\mathrm{cd}=\frac{\underset{k=1}{\overset{a}{\mathrm{\Σ}}}\underset{j=1}{\overset{b}{\mathrm{\Σ}}}{r}_{\mathrm{kj}}}{\mathrm{ab}}\·$

6.根据权利要求1所述的基于几何结构信息的多目标优化方法，其特征在于，所述抽取当前种群的前沿解集，并对其进行参数化，若优化目标函数空间为二维，则其前沿解集上每个粒子的详细参数化过程包括步骤：

对前沿解集在优化目标函数空间根据第一维进行升序排序；

求出前沿解集中相邻粒子在优化目标函数空间的欧氏距离，计算所有欧氏距离值之和得到总长度；

对于某个粒子，计算在其前面粒子的欧氏距离值之和，将得到的欧氏距离值除以总长度即为该粒子对应的参数值u。

7.根据权利要求6所述的基于几何结构信息的多目标优化方法，其特征在于，所述抽取当前种群的前沿解集，并对其进行参数化，若优化目标函数空间为三维，则其前沿解集上每个粒子的详细参数化过程包括步骤：

找到两个正整数(p,q)，使得p*q的值为前沿解集的大小，将一维的前沿解集向量看成是二维的p*q的矩阵，p为参数u方向上的个数，q为参数v方向上的个数；

使用二维优化目标函数空间的方法对u参数方向进行参数化；

在v方向上，即二维矩阵p*q的每一列q个粒子根据粒子在优化目标函数空间第二维进行升序排序，使用二维优化目标函数空间的方法对v参数方向进行参数化；

计算每一列u的平均值，该列所有粒子的u值即为此平均值；

计算每一行v的平均值，该行所有粒子的v值即为此平均值。

8.根据权利要求1所述的基于几何结构信息的多目标优化方法，其特征在于，所述根据参数化的结果，拟合其非均匀有理B样条方程，如果优化目标函数个数为2，则拟合成非均匀有理B样条曲线，优化目标函数个数为3则拟合为非均匀有理B样条曲面。

9.根据权利要求1所述的基于几何结构信息的多目标优化方法，其特征在于，所述利用牵引点集对所有被支配解进行优化，包括步骤：

步骤5.1、初始化被支配解集中当前被支配解的速度为0；

步骤5.2、计算当前被支配解与前沿解集中所有非支配解在优化目标函数空间的欧氏距离，选择两个欧氏距离最近的非支配解，找到这两个非支配解对应的牵引点，计算两个牵引点与其对应的非支配解之间的欧氏距离和d1；

步骤5.3、利用两个欧氏距离最近的非支配解，按照速度更新公式改变当前被支配解的飞行速度，按照位置更新公式改变其优化变量的值，最后重新求解其对应的优化目标函数值；

二维速度更新公式为：

$\stackrel{\→}{v(i+1)}=\mathrm{\ω}\stackrel{\→}{v\left(i\right)}+C_1(\stackrel{\→}{x_{p_1}}-\stackrel{\→}{x})+C_2(\stackrel{\→}{x_{p_2}}-\stackrel{\→}{x});$

其中，

是当前被支配解在第i次迭代的速度，

是当前被支配解的当前位置，

和

是它的两个最近的非支配解的位置，ω是被支配解速度的权重，C₁和C₂是分别受两个非支配解影响的权重，

是当前被支配解在第（i+1）次迭代的速度；

三维速度更新公式为：

$\stackrel{\→}{v(i+1)}=\mathrm{\ω}\stackrel{\→}{v\left(i\right)}+C_1(\stackrel{\→}{x_{p_1}}-\stackrel{\→}{x})+C_2(\stackrel{\→}{x_{p_2}}-\stackrel{\→}{x})+C_3(\stackrel{\→}{x_{p_3}}-\stackrel{\→}{x});$

其中，

和

是当前被支配解的三个最近的非支配解的位置，C₁、C₂和C₃是分别受三个非支配解影响的权重；

二维和三维位置更新公式为：

$\stackrel{\→}{x(i+1)}=\stackrel{\→}{x\left(i\right)}+\stackrel{\→}{v(i+1)};$

其中，

是当前被支配解在第i次迭代时的位置，

是当前被支配解在(i+1)次迭代时的速度，

是当前被支配解在第(i+1)次迭代时的更新位置；

步骤5.4、计算当前被支配解与两个牵引点之间的欧氏距离之和d2；

步骤5.5、如果d1小于d2，说明被支配解未被牵引到更优区域内，此时返回步骤5.3，若d1不小于d2，则被支配解被牵引到更优区域，优化过程结束。

10.根据权利要求1所述的基于几何结构信息的多目标优化方法，其特征在于，所述对前沿解集作均匀化处理或向外延拓处理，包括步骤：

根据前沿解集中非支配解的参数值计算相邻非支配解间的距离值，找到距离值最大的两个非支配解，表明在整个前沿解集中它们之间最稀疏；

取这两个非支配解参数值的平均值作为即将插入点的参数值，根据非均匀有理B样条方程求出这个参数值对应的数据点，即在优化目标函数空间的目标函数值；

利用信赖域优化方法求解优化变量的值来逼近期望的优化目标值。

Images