粒子群算法自1995年被Kennedy和Eberhart提出后^[1]，由于其简单高效，收敛速度快，逐渐在优化算法中脱颖而出。随之粒子群算法被Coelloran应用到多目标优化问题上^[2-3]，得到了各界的广泛认可。后来更多学者对多目标粒子群算法(MOPSO)进行了更深入的研究：Raquel等^[4]使用了拥挤距离作为适应值对粒子进行排序；Li等^[5]通过整合粒子的GMR(global margin ranking)和粒子密度信息的方法来对粒子进行排序，能高效快速地选择出p_best和g_best粒子；除了对排序选择的方法进行优化外，有的利用聚合函数来对多个目标函数进行处理^[6]；有的使用了目标空间分解的方法来保障粒子最优解的多样性^[7]；以及将MOPSO和其他优化算法进行交叉使用，如教与学方法^[8]、差分方法^[9]等, 通过一定的比例调用MOPSO。

上述的这些方法不断地对MOPSO的多样性和全局收敛进行改进和优化，虽然取得了较好的研究成果，但是这些方法的改进策略只对MOPSO某一性能效果进行改进，其他的性能指标并没有同时得到很好的改进。例如：文献[7]中，多样性方面得到很大的提升，但是其收敛性方面仍有很大不足。所以，本文借鉴Pareto支配强度^[10]和连续变异的方法^[11]对基于目标空间分解方法的多目标粒子群算法^[7]进行优化，在给每个子区域分配粒子时，删除没有归属的粒子，给没分到粒子的区域初始化新的粒子，增加获得较优粒子的可能性。这样既从整体和局部两方面提升了粒子的收敛性，同时多样性也得到了一定的优化提升，这种新的算法被称为基于目标空间分解和连续变异的多目标粒子群优化算法(multi-objective particle swarm optimization based on decomposition and continuous mutation ，MOPSO/DC)。

1 多目标优化问题 1.1 基本概念

$\left\{ {\begin{array}{*{20}{l}} {\min y = F(x) = ({f_1}(x),{f_2}(x),\cdots,{f_m}(x))} \\ {{\rm{s.t}}.{g_i}(x) \leqslant 0,\;i = 0,1,2,\cdots,q} \\ {{h_j}(x) = 0,\;j = 1,2,\cdots,p} \end{array}} \right.$

(1)

式中： $x = ({x_1},{x_2},\cdots,{x_n})$ 为 $n$ 维决策变量； $m$ 为目标函数的个数； $g(x)$ 函数为目标函数的 $q$ 个不等式约束； $h(x)$ 为目标函数的 $p$ 个等式约束。所有这些满足条件的决策变量用集合 $\varOmega $ 表示， $Y = \{ F(x)|x \in \varOmega \} $ 为目标空间。接下来介绍4个关于多目标问题的重要定义。

定义1　Pareto支配。解 $d,e \in \varOmega $ ， $d$ 支配 $e$ ，记为 $d \succ e$ ，满足下面的两个关系式：

$\left\{ {\begin{array}{*{20}{c}} {\forall i \in \{ 1,2,\cdots,m\},\,\,\, {f_i}(d) \leqslant {f_i}(e)} \\ {\exists i \in \{ 1,2,\cdots,m\},\,\,\, {f_i}(d) < {f_i}(e)} \end{array}} \right.$

(2)

定义2　Pareto最优。如果 $x$ 是Pareto最优解，则在 $\varOmega $ 中， $\neg \exists z \in \varOmega $ 使 $z \succ x$ 成立。

定义3　Pareto最优解集( ${\rm PS}$ )：

${\rm PS} = \{ x \in \varOmega |\neg \exists z \in \varOmega ;z \succ x\} $

(3)

定义4　Pareto前沿( ${\rm{PF}}$ )：

${\rm{PF}} = \{ F(x)|x \in {\rm{PS}}\} $

(4)

1.2 粒子群优化

粒子群算法中的粒子由速度信息和位置信息组成，更新公式为

$V_k^{t + 1} = WK_k^t + {C_1}{R_1}(P - X_k^t) + {C_2}{R_2}(G - X_k^t)$

(5)

$X_k^{t + 1} = V_k^{t + 1} + X_k^t$

(6)

式中：k指粒子群中第k个粒子；t为当前迭代次数；W为权衡局部搜索和全局搜索的参数， $W \in [0.1,0.9]$ ； ${C_1}$ 和 ${C_2}$ 为学习因子，其大小都为2； ${R_1}$ 和 ${R_2}$ 都是[0，1]之间的随机数；P指当前粒子最好位置p_best；G指引导粒子g_best的位置； $V_k^t$ 指在第t次迭代中第k个粒子的速度； $X_k^t$ 指在第t次迭代中第k个粒子的位置。

2 MOPSO/DC算法 2.1 目标空间分解

把目标空间Y分成M个子区域Y₁，Y₂，···，Y_M。令 $j \in \{ 1,2,\cdots,M\} $ ，对任一给定的第j个子区域，每个目标函数在所有目标函数中所占有的权重 ${a_i}\left( {i = 1,2, \cdots ,m,}\right.$ 且 ${\displaystyle\sum\limits_{i = 1}^m { {{a_i}} = 1} } )$ 所组成的向量 $({a_1},{a_2},\cdots,{a_m})$ 定义为该子区域的中心向量 ${{{A}}_j}$ 。当目标函数个数 $m = 2$ 时，则第 $j$ 子区域的中心向量 ${{{A}}_j}$ 表示为 $\left[\displaystyle\frac{{j - 1}}{{M - 1}}\;\;\;\; 1 - \displaystyle\frac{{j - 1}}{{M - 1}}\right]$ ；当 $m = 3$ 时，进行两层循环，令 ${k_1}$ 为第1层(外层)循环变量， ${k_2}$ 为第2层(内层)循环变量， ${k_1}$ 从0取到h， ${k_2}$ 从0取到 $h - {k_1}$ ，每次循环令 $u = {k_2} + 1 + \displaystyle\sum\limits_{{k_1} = 1}^h {(h + 2 - {k_1})} $ ，其中h和M满足：当 $C_{h + 2}^2 \geqslant M$ 时h取最小值，则第u个子区域的中心向量为 $\left[\displaystyle\frac{{{k_1}}}{h}\;\;\;\;\displaystyle\frac{{{k_2}}}{h}\;\;\;\;1 - \displaystyle\frac{{{k_1}}}{h} -\displaystyle \frac{{{k_2}}}{h}\right]$ ；当 $m > 3$ 时，进行 $m - 1$ 层循环， ${k_1}$ 为第一层循环变量(最外层)，依次由外向内，令 ${k_i}$ 为第i层循环变量，k_m-1为第 $m - 1$ 层(最内层)循环变量， ${k_1}$ 从0取到h，k₂从0取到 $h - {k_1}$ ， ${k_i}$ 从0取到 $(h - {k_1} - {k_2} - \cdots{k_i})$ ，最里层循环变量k_m-1从0取到 $(h - {k_1} - {k_2} - \cdots{k_{m - 2}})$ ，每次循环令

$ \begin{split} u = & \displaystyle\sum\limits_{{k_1} = 0}^h {\displaystyle\sum\limits_{{k_2} = 0}^{h - {k_1}}\cdots}\\ & \displaystyle\sum\limits_{{k_i} = 0}^{h - {k_1} - \cdots - {k_{i - 1}}} \cdots\displaystyle\sum\limits_{{k_{m - 3}} = 0}^{h - {k_1} - \cdots - {k_{m - 4}}} {{{\left[\displaystyle\sum\limits_{{k_{m - 2}} = 1}^{h - {k_1} - \cdots - {k_{m - 3}}} {(h + 2 - \displaystyle\sum\limits_{l = 1}^{m - 2} {{k_l}} ) + {k_{m - 1}} + 1} \right]} } } \end{split} $

其中参数h、m和子区数M满足关系：h取满足 $C_{h + m - 1}^{m - 1} \geqslant M$ 时的最小值，则第u个子区域的中心向量为 $\left[\displaystyle\frac{{{k_1}}}{h}\;\displaystyle\frac{{{k_2}}}{h}\;\cdots\;1 - \displaystyle\frac{1}{h}\sum\limits_{l = 1}^{m - 1} {{k_l}} \right]$ 。当 $m \geqslant 3$ 时，一共产生 $C_{h + m - 1}^{m - 1}$ 个中心向量，这些向量的下标为 $u = 1,2, \cdots C_{h + m - 1}^{m - 1}$ 。若 $C_{h + m - 1}^{m - 1} > M$ ，则从下标 $u = 2$ 开始，以 $\displaystyle\frac{{M - 2}}{{C_{h + m - 1}^{m - 1} - M}}$ 为步长依次删除 $C_{h + m - 1}^{m - 1} - M$ 个下标对应的中心向量，然后将剩下的M个中心向量依次作为 ${{{A}}_j}$ ，j=1, 2, ···, M。通过每个子区域的中心向量找出该子区域的T个相邻的子区域, 参考指标为两个中心向量的余弦值。

2.2 粒子的分类与更新

在对粒子进行分配时，通过参考点确定每个粒子的方向向量^[7]。令参考点为 $R({r_1},{r_2},\cdots,{r_m})$ ，其中 ${r_i} = \min \{ {f_i}(x)|x \in \varOmega \} ,i = 1,2,\cdots,m$ ，m为目标函数个数。例如：如图1所示， ${f_1}$ 和 ${f_2}$ 表示两个目标函数，m=2，则图中粒子C的方向向量为从点 $R({r_1},{r_2})$ 指向点C的向量X，X=C−R。通过比较粒子方向向量和所有子区域中心向量的余弦值，确定该粒子属于最大余弦值所对应的中心，向量子区域。

在对粒子进行分配时会有以下两种特殊情况：

1)有的子区域粒子数大于子区域的容量Vol时，通过适应值进行取舍，适应值公式为

$f = aT + {\rm{CD}}$

(7)

式中：T是Pareto支配强度，其数值为当前粒子支配的粒子数，在适应值计算中加入了Pareto支配强度T，增强每个子区域中的粒子趋向真实PF的能力；参数a表示支配强度对适应值的影响程度， $a = \displaystyle\frac{2}{M}$ ，M指目标空间子区域数；CD为拥挤距离，通过每个目标函数值对粒子进行排序，序列两端粒子在当前目标函数中的拥挤距离设为5，其他粒子在当前目标函数中的拥挤距离为在序列中该粒子前后两粒子的目标函数值之差的绝对值，最后将求出的该粒子在每个目标函数中计算的拥挤距离之和作为当前该粒子的拥挤距离CD。

计算该子区域包含的粒子的适应值，从大到小排序，选择序列中前30%的粒子，再从这些粒子中选择离所在子区域中心向量最近的Vol个粒子，多余的劣质粒子删除。

2)当子区域粒子数少于该子区域容量Vol时，由于情况1)中已经对劣质的粒子进行删除，为了尽可能保证整体的优越性，所以重新初始化该子区域所缺数目的粒子作为该区域的粒子，增加获得较优粒子的可能性，增加粒子多样性。然后计算新粒子的目标函数值。

2.3 连续变异操作

当只采用一种方法进行变异时，不能兼顾全局粒子和局部粒子的特性。在MOPSO/DC中提出差分+柯西^[12]+高斯^[13]连续的变异策略，公式为

$\left\{ \begin{gathered} {X_{{\rm{diff}}}} = X + 0.5({X_1} - {X_2}) \hfill \\ {X_{{\rm{Cauchy}}}} = \frac{{c{{(t)}^2}}}{{\text{π} ({X^2} + c{{(t)}^2})}} \hfill \\ {X_{{\rm{Gauss}}}} = \frac{1}{{g(t)\sqrt {2\text{π} } }}\exp (\frac{{ - {X^2}}}{2}) \hfill \\ \end{gathered} \right.$

(8)

式中：X为当前引导粒子g_best的位置；X₁和X₂是从EPOP中随机选出的两个不同粒子的位置。这样在变异时可保留全局粒子的一部分性质，起到一定信息交流的作用，增加了多样性。t为当前迭代的次数，g_max为最大迭代次数， $g(1) = 2$ ， $g(t + 1) = $ $ \displaystyle\frac{{g(t) - 1}}{{g_{\rm{max}}}}$ ； $c(1) = 1$ ， $c(t + 1) = \displaystyle\frac{{c(t) - 1}}{{g_{\rm{max}}}}$ ； ${X_{{\rm{diff}}}}$ 、 ${X_{{\rm{Cauchy}}}}$ 和 ${X_{{\rm{Gauss}}}}$ 分别是对g_best的位置进行变异操作之后产生的新粒子。每次变异后让产生的新粒子和当前的g_best进行比较，选择支配权优先的作为g_best，然后进行后续变异操作，直到3次变异结束，最终确定引导粒子g_best。

高斯变异步长较短，能很好地吸取局部粒子性质，柯西变异具有相对大的步长，能进行较大范围的变异，具有全局范围变异的特性，这样能很好地产生全新的粒子。

下面对标准柯西分布 ${y_{\rm{c}}}$ 和标准高斯分布 ${y_{\rm{g}}}$ 进行对比分析，根据 ${y_{\rm{c}}}$ 和 ${y_{\rm{g}}}$ 的单调性和对称性，对其第一象限进行分析： ${y_{\rm{c}}} = \displaystyle\frac{1}{{{\rm{{\text{π}} }}({x^2} + 1)}}$ , ${y_{\rm{g}}} = \displaystyle\frac{1}{{\sqrt {2{\rm{{\text{π}} }}} }}{\rm{exp(}} - \displaystyle\frac{{{x^2}}}{2}{\rm{)}}$ ，求反函数 $x_{\rm{c}}^2 = \displaystyle\frac{1}{{{\rm{{\text{π}} }}{y_{\rm{c}}}}} - 1$ 和 $x_{\rm{g}}^2 = - 2\ln (\sqrt {2{\rm{{\text{π}} }}} {y_{\rm{g}}})$ ，求两反函数之差 $(x_{\rm{c}}^2 - x_{\rm{g}}^2)$ ，并用 $y$ 同时替代 ${y_{\rm{c}}}$ 和 ${y_{\rm{g}}}$ ，即 $T = \displaystyle\frac{1}{{{\rm{{\text{π}} }}y}} - 1 +2\ln$ $ (\sqrt {2{\rm{{\text{π}} }}} y)$ ；对T求关于 $y$ 的导数，即 $T' = - \displaystyle\frac{1}{{{\text{π}} {y^2}}} + \displaystyle\frac{2}{y}$ , T在 $(0,\displaystyle\frac{1}{{2{\rm{{\text{π}} }}}})$ 上单减，在 $(\displaystyle\frac{1}{{2{\rm{{\text{π}}}}}}, + \infty )$ 上单增，所以 $T(0) > T(\displaystyle\frac{1}{{5{\rm{{\text{π}} }}}}) > T(\displaystyle\frac{1}{{2{\rm{{\text{π}} }}}})$ ；又因为 $T(\displaystyle\frac{1}{{5{\rm{\text{π} }}}}) > 0$ , 所以 $T(0) > 0$ ；所以在 ${y_{\rm{c}}}$ 和 ${y_{\rm{g}}}$ 同时趋向于0时, $x_{\rm{c}}^2 > x_{\rm{g}}^2$ , 即 ${x_{\rm{c}}} > {x_{\rm{g}}}$ , 所以柯西分布在x轴上覆盖的面积比高斯分布大。从图2中可直观地看出，柯西分布具有较宽的变异范围，粒子具有较大范围的变异。

用这3种方法进行变异操作，既增强了g_best的引导能力，同时g_best吸取局部或全局粒子的性质，促进了粒子之间信息交流共享，增加了多样性。

2.4 MOSPO/DC步骤

1)初始化2N个粒子的位置，这些粒子的位置集合记为POP，初始化N个速度，这些速度集合记为V，设目标函数的个数为m，当前迭代次数为t，初始化最大迭代次数g_max、每个子区域容量 ${\rm{Vol}} \in [1,3]$ 、邻域个数T以及子区域数目M。

2)进行目标空间分解操作。

3)计算粒子群POP目标函数值，选出每个目标函数的最小值作为参考点R用于计算每个粒子的方向向量，方向向量的计算参考图1。

4)进行粒子的分类与更新操作，接着将当前所有子区域粒子的位置信息存入集合EPOP中，最后更新参考点R，并重新计算每个粒子的方向向量，清空POP备用。

5)g_best的选择：产生[0，1]之间的随机数r₁,当r₁>0.8时，随机从EPOP集合中选择一个粒子作为当前粒子的g_best；否则，在该粒子所在区域的邻域内操作。首先计算该粒子邻域中每个子区域的中心向量和该子区域中粒子方向向量的余弦值，然后选出余弦值最大的子区域中的粒子作为当前粒子的g_best。过多地使用变异会让算法变得更加随机，降低了算法的效率，所以产生[0,1]之间的随机数r₂,当r₂<0.6时，对选出的引导粒子位置进行连续变异操作，确定最终的引导粒子g_best。否则，不对g_best进行变异操作。

6)p_best的选择：产生[0,1]之间的随机数r₃，当r₃>0.8时，随机从EPOP中选择一个粒子；否则，从该粒子的所在区域的邻域中随机选择一个粒子。让当前粒子和这个随机选出的粒子进行比较，选出支配权优先的作为当前粒子的p_best。

7)通过上述选出的g_best和p_best，以及粒子速度更新式(5)和位置更新式(6)产生下一代新的粒子群体，这些新粒子位置集合记为NPOP，下一代新的速度集合重新覆盖集合V。

8)如果当前迭代次数t大于最大迭代次数g_max，则循环结束，输出EPOP作为最优解集，否则把NPOP和EPOP合并放入POP中，然后跳转至3)，继续循环。

3 实验分析

为了测试所提出的MOPSO/DC算法的性能，将其与目前较流行的MOPSO/D^[7]、MOPSOTL^[8]、MOPSODE^[9]和NNIA^[14]算法进行对比，实验中同样选择了文献[7]中的6个测试函数，依次为F1，F2···F6。为了定量比较5种算法的性能，采用以下3种广泛使用的性能指标：HV^[15]，测量了粒子收敛于真实PF的效果以及最优解集的多样性，所获得值越大越好；GD^[16]，指算法获得的PF到真实PF的距离，数值越小越接近真实PF，效果越好；IGD^[15]，显示真实PF到算法所得到PF的距离，值越小越能说明算法得到的PF上的点可以分散地均匀地收敛于真实PF上。

MOPSO/DC参数设置：当目标函数个数 $m = 2$ 时， $M = 500$ 。当 $m = 3$ 时， $M = 1\;035$ 。Vol=1， $N = M{\rm{Vol}}$ ， $T = [0.1N]$ ，粒子维度n=10，最大循环次数g_max=1 000，一共进行了30次实验求各性能指标的均值和标准差。其他算法参数请参考文献[7]。

MOPSO/DC与其他算法关于这6个测试函数的3种性能指标对比结果见表1，表中给出了通过MOPSO/DC和其他4个算法获得关于这6个问题的IGD、GD和HV的平均值和标准差。从表1中可以看出：1)MOPSO/DC获得的IGD平均值远小于其他4种算法；2)MOPSO/DC获得的GD平均值都是最佳的，除F3函数外，不过关于F3的GD均值，本算法也优于MOPSO/D算法^[7]；3)通过MOPSO/DC获得的HV值总体来说比其他算法获得的HV值好(除F1函数)，说明MOPSO/DC的收敛性以及解的多样性也相当好。整体来说，MOPSO/DC的这3个性能指标优于其他4个算法，对MOPSO/D算法进行了很好的改进。

MOPSO/DC算法的仿真图如图3，结果显示该算法得到的PF和真实PF完全重合。

通过这些对比实验得出结论：MOPSO/DC对原算法^[7]的收敛性和多样性进行了很好地优化，增强了粒子趋向真实PF的效果，并且具有很好的稳定性。

4 结束语

本文提出的算法采用目标空间分解与子区域粒子更新和对引导粒子g_best的位置进行连续变异的方法来提高粒子多样性和收敛性，对原算法的不足之处进行了很好地弥补和优化。通过和4个不同种类的优化算法对不同种测试函数进行多种性能的测试对比实验证明，MOPSO/DC整体效果最好，具有很高收敛性，多样性方面也有提升。总之，MOPSO/DC对多目标优化问题的多个性能指标同时进行了很好的优化。未来的研究重点是，将该算法更好地适应更多目标的优化问题，增强其通用性，用于解决实际问题。