许多工程实际优化问题, 例如动态车辆路径规划、动态生产调度和最优控制器的动态设计等, 其待优化目标的个数以及目标函数的变量维数和参数会跟随环境发生动态变化, 从而导致最优解随之发生改变.如果待优化目标≥2, 且多个动态目标之间存在冲突, 则称之为动态多目标优化问题(Dynamic multi-objective optimization problems, DMOPs).本文主要考虑一种目标函数参数随环境动态发生改变的动态多目标优化问题.不失一般性, 该问题描述为

其中, $x$是决策变量, $S\subset {\bf R}^{ N}$是决策空间; ${\pmb F}:({\pmb x}, \alpha(t))$ $\to {\bf R}^{M}$包含$M$个目标函数$f_i ({\pmb x}, \alpha(t))$, ${\bf R}^{M}$是目标空间. $\alpha(t)$为环境相关的动态参数.

通常认为, 环境变化依时间轴呈现离散特性, 即环境变化仅发生在某些非连续的时间点, 则相应的动态参数可以被离散化为$\alpha(k)$, $k=1, 2, \cdots$, $K$, 并且满足$\forall\, \alpha (k)$, $\alpha (k)\ne \alpha(k+1)$.在任一环境中, $\alpha(k)$保持不变, 相应的目标函数也不发生变化.由此, 上述动态多目标优化问题被转化为$K$个静态多目标优化问题, 记为, $\min {\pmb F}({\pmb x}, \alpha (K))\rangle$.针对任一环境下的静态多目标优化问题, 获得的Pareto最优解集记为PS($k$), $k=1$, $2, {\cdots}, K$, 相应的Pareto前沿记为PF($k$), $k=1, 2$, ${\cdots}$, $K$, 则动态多目标优化问题的求解目标就是要在环境发生动态改变之前, 即$\alpha(k+1)$出现之前, 尽快找到满足$\min {\pmb F}({\pmb x}, \alpha (k))$的Pareto最优解PS($k$).

为有效解决上述动态多目标优化问题, 设计一种快速有效的优化算法, 以适应持续的环境变化成为其关键技术.目前, 国内外研究学者已经提出了众多解决动态多目标优化问题的方法, 其解决思路可以被划分为两大类:跟踪环境动态变化后的最优解(Track the moving optima, TMO)^[1-11]和动态鲁棒最优解^[12-16].

TMO方法的本质是当检测到环境发生变化时, 重新初始化种群, 并触发相应的寻优过程.该方法的核心在于如何在优化过程中, 平衡对决策空间的探索和开发过程, 从而兼顾种群的多样性和收敛的快速性, 保证在$\alpha(k)\to \alpha (k+1)$的有限时间内, 找到满足$\min {\pmb F}({\pmb x}, \alpha(k))$的Pareto最优解.现有研究主要集中在环境变化的检测、历史信息的有效利用和种群多样性保持等三方面.

针对动态环境的检测问题, Farina等^[1]将每一代种群中随机选出的10{%}个体进行重新评价.如果重新评价前后的个体适应值不同, 则表示环境发生变化.为有效保持种群多样性, 研究人员往往在探测到环境变化后, 将随机或变异产生的新个体, 按一定比例替代原种群, 并在新环境下评价生成的新种群.

动态多目标优化问题中, 历史信息的有效利用主要是指, 如何利用不同环境下获得的最优解PS($k$)来提高新环境下的寻优效率.一方面, 采用记忆机制来保存不同动态环境下的历史最优解^[1].当检测到相似动态环境时, 通过历史信息的重用来提高算法性能.记忆增强型动态多目标进化算法^[2], 则通过对不同动态环境的串式记忆方法, 利用历史最优解来加速对新环境的跟踪能力.记忆机制能有效解决具有周期特性的动态多目标优化问题.另一方面, 根据连续动态环境下获得的历史信息, 采用预测机制来估计未来动态环境中的较优解, 并采用该预测解集构成或部分替代新环境下初始种群. Hatzakis等^[3]提出的前馈法预测模型中, 采用自回归预测模型估计相邻动态环境下的PS, 进而生成新环境下的初始种群. Zhou等^[4]采用时间序列预测模型, 通过加入白噪声等方法来生成新环境下的初始种群.进而, 又提出基于流形的种群预测方法^[5], 将任一环境下的Pareto解集, 转化为中心点和流形两部分加以描述.通过对这两部分分别加以预测, 得到下一时刻的初始种群.虽然这种预测方法具有较小的时间和空间复杂度, 但是仅适用于${\rm PS}(k)$, $k=1, 2, {\cdots}, K $具有相似流形的动态优化问题.

多种群优化方法的引入, 为提高动态多目标优化问题的求解效率提供了一条新途径. Goh等^[6]采用协同进化算法, 通过处于子目标空间中各个子种群之间的竞争, 以及获胜者之间的合作进化, 实现动态环境下最优Pareto前沿的跟踪.基于多种群协同的动态多目标粒子群优化算法^[7]中, 通过种群之间的竞争实现对解空间的探索, 当竞争失效后, 种群之间自适应切换到协作模式, 对解空间进行搜索.通过种群竞争和种群协作两种模式之间的合理配合, 实现动态多目标优化问题的高效求解.

可见, TMO方法的本质在于跟踪环境的变化, 在有限时间内尽快寻找到适应于每一环境的最优解.虽然TMO为动态优化问题提供了有效的解决思路, 但是仍存在以下困难. 1)当动态优化问题的目标函数比较复杂以及目标空间或者决策空间具有较高维数时, 会导致每一环境下的最优解很难在有限时间内获得. 2)当环境变化较快时, 若PS($k$)与PS$(k$ -1)存在显著差别, 面向某些实际工程优化问题, 采用PS($k$)虽然更逼近该环境下的最优解, 但是可能导致较高的人员操作代价和资源成本, 甚至于无法在有限时间内执行该优化解.因此, 当每次环境变换时, 都去跟踪一个新的解往往是不可行的.

基于以上考虑, Yu等^[12]首次提出了基于时间的动态鲁棒最优解(Robust optimization over time, ROOT)的概念.其核心思想在于寻找一组鲁棒解集, 每个鲁棒解可以以一定满意度适用于多个连续动态环境. Jin等^[13]进一步给出了求解ROOT问题的一个通用框架, 包含基于种群的优化算法、数据库、适应度近似策略和适应度预测策略.其中, 解的鲁棒性由其过去和将来的适应度性能来估计. Fu等^[14-15]为寻找随时间变化的动态鲁棒优化解, 定义了生存时间和平均适应度两个可计算性能指标来描述解的鲁棒性.以上关于求解动态鲁棒最优解的进化方法都是针对动态单目标优化问题提出的.面向动态多目标优化问题, Chen等^[16]提出了描述Pareto解环境适应性的新型鲁棒Pareto最优解(Robust Pareto optima over time, RPOOT)的定义, 给出了鲁棒Pareto解的生存时间和平均适应度描述.

求解RPOOT和已有的ROOT有很大的不同. 1）由于目标个数从单一目标扩展成了多个目标, 而且多个目标之间往往存在矛盾性.满足单一目标最优的解往往为边界点.兼顾多个目标性能的折中解往往不能在所有目标上具有最优性能.因此, 某一较优解的鲁棒性不能仅考虑其在连续动态环境下某一目标的适应度值, 而要度量其所有子目标的相对变化程度. 2）求解鲁棒解的过程中往往需要预测未来连续动态环境下的适应度. ROOT仅仅是单一目标值构成的时间序列.由于多目标优化问题中, 彼此非占优解都会成为最优解构成Pareto最优解集, 所以RPOOT中, 多个非占优解的目标时间序列构成了一个多时间序列预测问题.针对上述问题, 本文基于Chen等^[16]提出的时间鲁棒性和性能鲁棒性评价指标, 将原有目标函数转化为兼顾多个连续动态环境下解的平均适应度函数.通过构造新型学习预测策略, 实现动态多目标鲁棒最优解集的求解.本文所提算法适用于解决动态多目标优化问题, 特别是环境参数动态变化频繁但不剧烈, 或者相邻动态环境下解的切换代价较大的动态优化问题.因为对于前者, 如果用传统动态多目标优化方法求解, 通常需要在每个新环境下, 通过重新激发寻优过程, 获得适应该环境的Pareto最优解.这可能导致较高的计算代价; 对于后者, 切换相邻动态环境下的最优解需要付出较大的资源成本, 甚至可能无法在有限时间内执行该优化解.本文给出的动态多目标鲁棒进化优化方法的核心思想是寻找一组鲁棒Pareto解集, 每个鲁棒Pareto解可以以一定满意度适用于多个连续动态环境, 这样可以避免频繁地求解和执行每个动态环境下的Pareto最优解.

本文第1节讨论了动态多目标优化问题的相关工作, 并且给出了动态多目标鲁棒优化问题的定义; 在第2节详细描述了本文给出的算法框架, 然后将所提出的算法框架应用到现有的基准函数问题中; 第3节的实验结果说明算法的有效性; 第4节给出结论和将来的工作.

1.   动态鲁棒Pareto最优解的性能度量

动态多目标鲁棒优化问题区别于式(1)所示动态多目标优化问题, 其核心在于寻找一组鲁棒Pareto最优解集RPS = ${{\rm RPS}(1), {\cdots}, {{\rm RPS} (LS)}}$, $LS\le K$, 其中, 每个鲁棒Pareto最优解RPS($i$)适用于多个连续动态环境.而传统动态多目标优化问题中, 评价Pareto解优劣仅仅依赖于其在某一动态环境下的收敛性和分布性, 未考虑该Pareto解在未来连续动态环境中的适用性.在动态多目标鲁棒优化问题中, 如何评价每个Pareto解的环境适应性, 以及其逼近每个动态环境下真实Pareto前沿的收敛性至关重要.由此, 面向式(1)所示动态多目标优化问题, 给出其动态多目标鲁棒Pareto解的时间鲁棒性和性能鲁棒性两类新型性能指标.

1.1   时间鲁棒性

定义$L_{i}$为鲁棒Pareto最优解集RPS中任意鲁棒Pareto解RPS($i$)的生存时间.它表示在未来$L_{i}$个连续动态环境下, RPS($i$)是依一定稳定性阈值$\eta$的满意解, 反映了RPS($i$)在连续动态环境下的时间鲁棒性.这里

其中, $\Delta (l(x_i^j))$表示鲁棒Pareto解RPS($i$)中某个个体$x_i^j$在相邻动态环境中的相对适应度差值.差值越小, 表示解的稳定性越好; 反之, 则解的稳定性越差. $\hat{{\pmb F}}(x_i^j, \alpha (k+l(x_i^j)))$表示$x_i^j $在未来动态环境$\alpha(k+l(x_i^j))$中的适应度预测值. $\eta$为稳定性阈值, 反映了决策者对$x_i^j $在相邻动态环境中的相对适应度差值$\Delta(l(x_i^j))$的最小容忍程度.通常, 稳定性阈值取决于决策者的偏好.

显然, 动态鲁棒Pareto解的时间鲁棒性取决于其在相邻动态环境中的相对适应值和稳定性阈值.它反映了任一鲁棒Pareto解RPS($i$)对连续动态环境的适应能力.

1.2   性能鲁棒性

虽然动态多目标鲁棒优化问题的寻优目的是要找到适用于多个连续动态环境的鲁棒Pareto解, 但是鲁棒Pareto解逼近其所适用的多个连续动态环境下多目标优化问题的真实Pareto前沿性能也是其重要的性能指标.假设任一鲁棒Pareto解RPS($i$)至少适用于$T$个连续动态环境.记$T$为固定时间窗, 反映了决策者对鲁棒Pareto解的环境适应能力的期望.定义鲁棒Pareto解RPS($i$)中的某个个体$x_i^j$在$T$时间窗内的平均适应度值

其中

鲁棒Pareto解的平均适应度反映了其所在时间窗内解的性能, 即性能鲁棒性.寻找鲁棒Pareto最优解RPS($i$)的目标就是最小化该平均适应度值, 即$\mathop{\min }_{x_i^j \in {\rm RPS}(i)} {\pmb F}^{\rm ave}(x_i^j, \alpha(k))$.

下面举例说明鲁棒Pareto解的时间鲁棒性和性能鲁棒性含义.考虑文献[1]给出的动态多目标测试函数FDA1和FDA3, 其中, FDA1的Pareto最优解PS($k$)随时间$k$而变化, 但相应的Pareto前沿PF($k$)不随时间变化; FDA3的PS($k$)和PF($k$)都随时间$k$而变化.以$t=0.1, 0.2, 0.3$三个连续动态环境为例, 相应的真实Pareto最优解PS($k$)和Pareto前沿PF($k$) ($k=1, 2, 3)$如图 1和图 2所示.如果固定时间窗$T=3$, 将$k=1$时刻的PS(1)作为鲁棒Pareto解RPS(1), 即RPS(1) = PS(1), 则其在时间窗内, 两个函数的鲁棒Pareto解RPS(1)的平均适应度曲线如图 3所示.

图 1  函数FDA1在三个时刻的真实Pareto解与前沿的比较

Fig. 1  True Pareto solutions and Pareto fronts of FDA1 in three environments

下载: 全尺寸图片 幻灯片

图 2  函数FDA3在三个时刻的真实Pareto解与前沿的比较

Fig. 2  True Pareto solutions and Pareto fronts of FDA3 in three environments

下载: 全尺寸图片 幻灯片

图 3  函数的平均适应度值与时刻$k=1, 2, 3$的真实Pareto前沿的比较

Fig. 3  Average fitness and true PF($k$) with $k=1, 2, 3$

下载: 全尺寸图片 幻灯片

鲁棒Pareto解中某个个体的生存时间取决于其在相邻动态环境中的相对适应度差值, 即式(4)中定义的$\Delta(l(x_i^j))$.本例中所取鲁棒Pareto解RPS(1)与相邻的$k=2, 3$环境时刻的相对适应度差值记为$\Delta(1)$和$\Delta(2)$, 如图 4和表 1所示.其中, $\hat {{\pmb F}}(x_1^j, \alpha (2))$, $\hat{{\pmb F}}(x_1^j, \alpha (3))$取其真实的适应度值.设定稳定性阈值$\eta=0.4$, 对于函数FDA1, 可得鲁棒Pareto解RPS(1)的生存时间为$L_1=2$; 而对于函数FDA3, 其鲁棒Pareto解RPS(1)的生存时间仅为$L_1=1$.如果增大稳定性阈值$\eta=0.7$, 则对于函数FDA1和FDA3, 可得鲁棒Pareto解RPS(1)的生存时间为都为$L_1=2$.显然, 稳定性阈值$\eta$体现了鲁棒Pareto解对相邻连续动态环境的适用程度. $\eta$越大, 鲁棒Pareto解的生存时间$L$越大; 反之, $\eta$越小, 解的生存时间$L$越小.

图 4  函数的平均适应度值与鲁棒解RPS在$t=1, 2, 3$时刻的鲁棒Pareto前沿的比较

Fig. 4  Average fitness and robust Pareto fronts with $k=1, 2, 3$

下载: 全尺寸图片 幻灯片

表 1  RPS(1)与相邻$k=2, 3$环境下解的相对适应度差值

Table 1  The relative fitness errors of RPS(1) and adjacent two environments

相对适应度差	$\Delta (1)$	$\Delta (2)$
FDA1	0.301	1.2177
FDA3	0.6728	1.7687

下载: 导出CSV 

| 显示表格

综上所述, 生存时间和平均适应度分别反映了任一鲁棒Pareto解在连续动态环境下的时间鲁棒性和性能鲁棒性.

2.   动态多目标鲁棒进化优化算法

2.1   模型转化

在动态多目标鲁棒优化中, RPOOT并不是寻找适用于每个动态环境的Pareto最优解PS($k$), 而是获得满足第1节定义鲁棒性的鲁棒Pareto最优解RPS($k$).它不仅要近似收敛于当前动态环境下的真实Pareto最优解, 而且要尽可能多地逼近后续连续动态环境下的真实Pareto最优解.这就要在进化优化过程中了解Pareto解PS($k$)的当前适应度, 及其在未来动态环境下的适应度值.基于第1节给出的时间鲁棒性和性能鲁棒性定义, 构建如下两种求解动态鲁棒Pareto最优解的模型.

模型1.  基于固定时间窗的模型转化

其中, ${\pmb F}^{\rm ave}(x, \alpha(k))$为动态环境时刻$k$条件下, 解$x$在固定时间窗$T$内的平均适应度值.该模型通过最小化固定时间窗内, 连续动态环境的平均适应度值, 来保证解$x$对相应连续动态环境下真实Pareto前沿的逼近能力.

模型2.  基于约束的模型转化

该模型将$x$对未来连续动态环境的适用程度转化成时间窗内的$T$个鲁棒约束.其中, 解$x$的当前适应度函数与将来相邻的$T-1$个环境时刻的适应度值的标准差小于决策者的稳定性阈值$\eta$.其中, 当$q=0$时, $\hat {{\pmb F}}({\pmb x}, \alpha(k+q))$表示当前时刻$k$的真实适应度值, 即$\hat {{\pmb F}}({\pmb x}, \alpha (k+q))={\pmb F}({\pmb x}, \alpha(k))$.

2.2   基于分解的RPOOT进化优化算法

上述两类转化模型中, 模型1将动态多目标鲁棒优化问题中的两类鲁棒性度量转化为固定时间鲁棒性下的解性能鲁棒性, 即某个解在当前环境及其未来$T-1$个连续动态环境下的平均适应度值.因此, 针对该类转化模型, 优化求解过程的核心是, 寻求在时间窗内同时对当前和未来动态环境具有良好收敛性的鲁棒Pareto最优解.不同于传统的DMOPs, 一方面, 优化目标是各个解在时间窗$T$内的平均适应值, 另一方面, 未来连续动态环境下解的适应度值是未知的, 需要采取时间序列预测方法来估计.预测方法见第2.3节.针对模型1的求解，采用传统的基于分解的多目标进化算法(Multi-objective evolutionary algorithm based on decomposition, MOEA/D).其中, 分解方法采用基于罚值的边界交集法(Penalty-based boundary intersection, PBI)^[17].具体伪代码见算法1.

不同于模型1, 模型2将解的时间鲁棒性条件作为约束, 限制其平均适应度范围, 从而将原有优化问题转化为一个动态多目标约束优化问题.针对该约束条件, 本文基于传统MOEA/D, 采用一种基于分解的无参数惩罚约束处理方法^[18].

算法1.基于分解的动态多目标鲁棒进化优化算法

1:开始;

2:初始化环境时刻$k=1$, 获得环境参数$\alpha (k)$;

3:随机初始化种群$P$, 估计种群中每个解${\pmb x}^i$的当前平均适应度值${\pmb F}^{\rm ave}({\pmb x}^i, \alpha (k))$, $i=1, 2, {\cdots}, N_{\rm pop}$;

4:初始化参考点${\pmb z}(k)=(z_1(k), {\cdots}, z_{ M} (k))^{\rm T}$, 其中, $z_j (k)$ = $\min \{f_j^{\rm ave} ({\pmb x}^1, \alpha (k)), {\cdots}, f_j^{\rm ave} ({\pmb x}^ {N_{\rm pop}}, \alpha (k))\}$;

5:根据均匀分布产生权重向量, $i$ = $1, {\cdots}, N_{\rm pop} $, 利用动态PBI分解法将目标分解为$N_{\rm pop}$个动态标量优化子问题, 其动态标量目标函数记为$f^{\rm PBI}({\pmb x}(k)|$ $\lambda, z)$;

6:计算任意两个权向量的欧氏距离, 为每个权向量选出最近的$r$个向量作为它的邻域, 记为${ B}(\lambda^i)=$ $(\lambda^{i1}$, ${\cdots}$, $\lambda^{ir})$, $i$ $=$ $1, {\cdots}, N_{\rm pop} $;

7:设进化代数$t=1$, 总进化代数Gen;

8: while $t $＜Gen do

9:   for $i =1$ to $N_{\rm pop}$ do

10:     根据差分进化和高斯变异产生子代$P_Z (t)$;

11:     对每一个$j=1, {\cdots}, M$更新参考点${\pmb z}(k)$;

12:     更新邻域解;

13:   end for

14: end while

15:输出$({\pmb x}^1(k), {\cdots}, {\pmb x}^{N_{\rm pop}}(k))$和$\{{\pmb F}^{\rm ave}({\pmb x}^1(k), \alpha (k)), $ $ {\cdots}$, ${\pmb F}^{\rm ave}({\pmb x}^{N_{\rm pop} }(k), \alpha (k))\}, $作为第$k$个动态环境下的鲁棒Pareto最优解;

16:如果$k>K$, 转步骤17, 否则$k=k+1$, 转步骤3;

17:结束.

在进化过程中, 约束处理的核心在于约束满意度的确定.由于在式(8)中, 受到预测误差及环境动态变化的影响, 随着$T$的增加, 同时满足$T$个约束条件的可行解将大大减少.为此，将上述约束条件转化为约束满足程度, 采用转化后的弱化约束条件来对个体加以比较.

基于上述转化后的约束条件, 定义约束满意度函数为

如果$V_s ({\pmb x}, \alpha(k))=0$, 说明个体$x$是满足式(8)中前$s$个约束条件的弱可行解; 否则$x$不是弱可行解.

模型2的主体算法框架采用算法1所示的MOEA/D.不同在于:步骤3的初始化和步骤10和更新种群操作, 增加了个体的约束满意度$V_s({\pmb x}, \alpha(k))$计算, 用来判断个体为时间窗内的一个弱可行域解.步骤12的更新邻域解操作中, 采用算法2所示的弱约束处理方法, 判断邻域解是否需要更新.

算法2.无参数惩罚约束方法

1:取$\forall\, {\pmb x}^i$, ${\pmb x}^j\in P(t)$;

2: for $s=T-1$ to 0 then

3:   if $V_s ({\pmb x}^i, \alpha (k))=0$, 且$V_s ({\pmb x}^j, \alpha (k))>0$ then

4:      ${\pmb x}^i$替代${\pmb x}^j$, 转步骤15;

5:   elseif $V_s ({\pmb x}^i, \alpha (k))=V_s ({\pmb x}^j, \alpha (k))=0$,

      且$f^{\rm PBI}({\pmb x}^i, k| \lambda, {\pmb z}) < f^{\rm PBI}({\pmb x}^j, k| \lambda, {\pmb z})$ then

6:     $x^i$替代$x^j$, 转步骤15;

7:   elseif $V_s ({\pmb x}^i, \alpha (k))>0, V_s ({\pmb x}^j, \alpha (k))>0$,

      且$s>0$ then

8:       $s=s-1$;

9:     elseif $s=0$ then

10:     令$f^{\rm{PBI}}({\pmb x}, k| {\pmb \lambda }, {\pmb z})=f^{\rm PBI}({\pmb x}, k| {\pmb \lambda }, {\pmb z})\, +$

       $V_{T-1} ({\pmb x}, \alpha (k))$;

11:     if $f^{\rm PBI}({\pmb x}^i, k| {\pmb \lambda }, {\pmb z})<f^{\rm PBI}({\pmb x}^j, k| {\pmb \lambda }, {\pmb z})_{ }$ then

12:       ${\pmb x}^i$替代${\pmb x}^j$, 转步骤15;

13:     end if

14:   end if

15: end for

2.3   预测机制

在动态多目标优化问题的求解中, 预测机制可以有效改善算法进化效率.在传统的TMO求解方法中, 预测的目的是在下一相邻动态环境下, 能够充分利用历史和当前最优Pareto解和其他已知信息, 来估计其动态环境参数的变化趋势, 构建其新环境下的初始种群.

与TMO方法中的预测目的不同, 在RPOOT中, 由于要考虑鲁棒Pareto解的鲁棒性, 确定其对未来动态环境适用情况.所以, 需要使用未来动态环境下该鲁棒Pareto解的适应度值.对于在线的实际动态优化问题, 未来动态环境下解的适应度值是未知的, 需要通过预测的方法来估计.因此, RPOOT中预测的本质是根据历史信息和动态环境(参数)变化趋势, 估计当前解在未来连续环境下的适应度值.在模型1和模型2中，优化目标均为平均适应度函数值${\pmb F}^{\rm ave}({\pmb x}, \alpha(k))$.根据式(6)可知, $x$在时间窗内连续动态环境下的平均适应度依赖于未来连续动态环境下$x$的适应度值.因此, 在动态多目标鲁棒优化过程中, 需要预测每代种群中所有解在多个目标上的未来适应度值, 从而使RPOOT的预测问题转化为一个多维时间序列预测问题, 其维数为$N \times M$.为降低计算复杂度, 预测方法的设计需考虑如下两方面, 1)预测方法的空间复杂度不能太高; 2)预测方法的时间复杂度不能太大, 否则会大大影响动态多目标优化算法的实时性.

本文采用移动平均(Moving average, MA)^[19]预测模型. MA是一种简单有效的时间序列预测模型.首先, 假设在第$k$个动态环境下, 针对种群中的每个解$x$构造一个长度为$m$的时间序列.其中, ${\pmb F}({\pmb x}, \alpha (i))$, $k$ -$m+1\le i\le k$是由$f_j ({\pmb x}, \alpha (k))$, $1\le j \le M$组成的该环境下解$x$的适应度向量.预测的目的就是根据上述时间序列来估计未来动态环境下的${\pmb F}({\pmb x}, \alpha(k$ $+$ $i))$, $1\le i\le T-1$. MA模型定义如下:

其中, $\varepsilon (k)\sim {\rm U}(0, \sigma ^2)$, 即方差为$\sigma^2$的高斯白噪声. $b$是对时间序列中相邻的前$m$个数据点的平均估计.

方差$\sigma ^2$通过下式估计

可见, $\hat {{\pmb F}}({\pmb x}, \alpha(k+i))$在未来连续动态环境下的适应度值可被预测为均值为$b$, 方差为$\hat{\sigma }^2$的值.

3.   实验分析与比较

实验仿真的主要任务是测试本文所提出算法的有效性.首先, 为有效度量和评价算法的鲁棒性、收敛性和分布性, 提出了两个性能指标.其次, 针对8个常用动态多目标优化测试函数, 设计了以下3组实验: 1)为分析说明本文所提模型对Pareto解集的性能影响, 通过实验对比分析了不同时间窗和不同稳定性阈值下的解性能; 2)为进一步分析本文算法对不同环境变化快慢及剧烈程度的适应能力, 基于提出的两类模型, 通过实验对比分析了不同动态环境频率和强度参数下鲁棒Pareto最优解集的性能; 3)为验证本文所提鲁棒模型转化的合理性, 将基于本文模型所得动态多目标鲁棒Pareto解集与文献[16]所得解集进行对比分析.

3.1   测试函数

本文采用8个常用动态多目标测试函数组, 包括由Farina和Deb^[1]构建的Fun 1~Fun 5 (FDA1-5), 以及由Goh与Tan^[6]构建的Fun 6~Fun 8 (DMOP1-3).具体函数定义如表 2所示.

表 2  测试函数

Table 2  Benchmark functions

函数	测试函数	变量维数	变量范围	说明
Fun 1	$ f_1 (X_{Ⅰ})=x_1, \ f_2 =g\times h$ $g(X_{Ⅱ})=1+\sum_{x_i \in X_{Ⅱ}} {(x_i -G(k))^2}$ $ h(f_1, g)=1-\sqrt {{{f_1 }/ g}}, \ G(k)=\sin (0.5\pi k)$	$ N=10$ $| X_{Ⅰ} | =1$ $| X_{Ⅱ} | $=9	$ X_{Ⅰ} =(x_1)\in [0, 1]$ $X_{Ⅱ} =(x_2, {\cdots}, x_N)\in [-1, 1]$	FDA1 in [6] PS随环境变化 PF不随环境变化
Fun 2	$f_1 (X_{Ⅰ})=x_1, f_2 =g\times h$ $g(X_{Ⅱ})=1+\sum_{x_i \in X_{Ⅱ} } {x_i^2 }$ $h(f_1, g)=1-\left({{{f_1 }/g}} \right)^{H(k)^{-1}}$ $H(k)=0.75+0.7\sin (0.5\pi k)$	$ N=10$ $| X_{Ⅰ} | =1$ $| X_{Ⅱ} | =9$	$X_\rm{1} =(x_1)\in [0, 1],$ $X_\rm{2} \in [-1, 1]$	FDA2 in [6] PS不随环境变化 PF随环境变化
Fun 3	$ f_1 (X_{Ⅰ})=\sum_{x_i \in X_\rm{1} }^n {x_i^{F(k)} }, \quad f_2 =g\times h$ $g(X_{Ⅱ})=1+G(k)+\sum_{x_i \in X_{Ⅱ}}{(x_i -G(k))^2}$ $h(f_1, g)=1-\sqrt {{{f_1 }/g}}, \quad G(k)=\left| {\sin (0.5\pi k)} \right|$ $F(k)=10^{2\sin (0.5\pi k)}$	$ N=10$ $| X_{Ⅰ} | =1$ $| X_{Ⅱ} |$ =9	$X_{Ⅰ}=(x_1)\in [0, 1]$ $X_{Ⅱ} =(x_2, {\cdots}, x_N)\in [-1, 1]$	FDA3 in [6] PS随环境变化 PF随环境变化
Fun 4	$ f_1 (x)=(1+g)\cos (0.5\pi x_1)\cos (0.5\pi x_2)$ $f_2 (x)=(1+g)\cos (0.5\pi x_1)\sin (0.5\pi x_2)$ $f_3 (x)=(1+g)\sin (0.5\pi x_1)$ $g(x)=\sum_{i=3}^n {(x_i -G(k))^2}$ $G(k)=\left| {\sin (0.5\pi t)} \right|$	$N=12$	$ x_i \in [0, 1]$ $\forall\, i=1, 2, {\cdots}, N $	FDA4 in [6] PS随环境变化 PF不随环境变化
Fun 5	$ f_1 (x)=(1+g)\cos (0.5\pi y_1)\cos (0.5\pi y_2)$ $f_2 (x)=(1+g)\cos (0.5\pi y_1)\sin (0.5\pi y_2)$ $f_3 (x)=(1+g)\sin (0.5\pi y_1)$ $g(x)=G(k)+\sum_{i=3}^n {(x_i -G(k))^2}$ $y_i =x_i^{F(k)}$ $G(k)=\left| {\sin (0.5\pi k)} \right|$ $F(k)=1+100\sin ^4(0.5\pi k)$	$N=12$	$ x_i \in [0, 1]$ $\forall\, i=1, 2, {\cdots}, N $	FDA5 in [6] PS随环境变化 PF随环境变化
Fun 6	$ f_1 (x_1)=x_1, f_2 =g\times h$ $g(x_2, {\cdots}, x_n, k)=1+9\sum_{i=2}^n {x_i^2 }$ $h(f_1, g)=1-\left({{{f_1 } /g}} \right)^{H(k)}$ $H(k)=0.75\sin (0.5\pi k)$+1.25	$N=10$	$ x_i \in [0, 1]$ $\forall\, i=1, 2, {\cdots}, N $	DMOP1 in [8] PS不随环境变化 PF随环境变化
Fun 7	$ f_1 (x_1)=x_1, f_2 =g\times h$ $g(x_2, {\cdots}, x_n, k)=1+\sum_{i=2}^n {(x_i -G(k))^2}$ $h(f_1, g)=1-\left({{{f_1 }/g}} \right)^{H(k)}$ $H(k)=0.75\sin (0.5\pi k)$+1.25	$N=10$	$ x_i \in [0, 1]$ $\forall\, i=1, 2, {\cdots}, N $	DMOP2 in [8] PS随环境变化 PF随环境变化
Fun 8	$f_1 (x_1)=x_1, f_2 =g\times h$ $g(x_2, {\cdots}, x_n, k)=1+9\sum_{i=2}^n {(x_i -G(k))^2}$ $h(f_1, g)=1-\sqrt {{{f_1 } / g}}$ $H(t)=\sin (0.5\pi k)$	$N=10$	$ x_i \in [0, 1]$ $\forall\, i=1, 2, {\cdots}, N $	DMOP3 in [8] PS随环境变化 PF不随环境变化

下载: 导出CSV 

| 显示表格

本文核心研究上述动态多目标优化问题所获得的Pareto鲁棒解对连续动态环境的适应能力.因此, 首先要确定其动态环境变化时刻$k$^{[1, 6]}.

其中, $\tau $是进化代数. $\tau_{ d}$表示环境保持不变时(即$k$不变), 所对应的进化代数规模, 用来控制环境的变化频率.环境的变化强度用$n_{ d}$来表示.

3.2   性能评价指标

RPOOT本质上是在进化过程中不断寻找对未来连续动态环境具有鲁棒性的Pareto解集.也就是在每个动态环境下, 搜索当前环境下的较优Pareto解, 而且该解能最大程度适用于未来多个动态环境.因此, 鲁棒Pareto解集不仅要考虑解的鲁棒性, 而且要考虑解在其适用环境中的收敛性.首先, 在时间尺度上定义平均生存时间评价解的鲁棒性能; 其次, 采用平均鲁棒逆向世代距离(Robust inverted generational distance, RIGD)度量鲁棒解的收敛性和分布性.

3.2.1   平均生存时间

动态多目标优化问题中, 鲁棒Pareto最优解的鲁棒性能往往体现为其可以适用的连续动态环境数目, 即Pareto解的生存时间.在此基础上, 为度量所有Pareto解集中解的鲁棒性, 定义其平均生存时间作为Pareto解集的时间鲁棒性测度.

其中, $K$是动态环境出现次数, $L_k$是鲁棒Pareto最优解RPS($k$)的生存时间.显然, 平均生存时间越长, 鲁棒Pareto解集的鲁棒性越好.

3.2.2   平均鲁棒逆向世代距离

鲁棒Pareto解集要兼顾其时间尺度和性能尺度下的鲁棒性能.后者主要体现为鲁棒Pareto解集对连续动态环境下的真实PF($k$)的逼近程度, 即其收敛性和分布性.传统多目标优化问题中, 逆向世代距离(Inverted generational distance, IGD)通常用来度量Pareto最优解的收敛性和分布性.因此, 将其推广到动态多目标鲁棒优化问题中, 定义新型平均鲁棒逆向世代距离来度量鲁棒Pareto解集的收敛性和分布性.

其中, ${\rm PF}(j)$和$|{\rm PF}(j)|$分别代表第$j$个动态环境下的真实Pareto解及其规模, 本文将${\rm PF}(j)$分别离散化为100个点(二目标)和300个点(三目标), MIGD$(k)$定义了第$k$个动态环境下鲁棒Pareto解RPF$(k)$与第$j$ 个动态环境下真实Pareto解PF$(j)$的IGD的最大值. RIGD$^{\rm RPOOT}$反映RPF$(k)$对其适用环境下的各个时刻MIGD测度的平均值.

3.3   实验参数设置

设计的全部仿真实验结果均采用MATLAB R2011b实现.针对所有测试函数, 解决两种模型的MOEA/D算法均采用实数编码、差分进化交叉和高斯变异.具体参数设置如下:

1) 群体规模$N_{\rm pop}$与权向量的设置:传统MOEA/D算法中, 权向量的数目等于种群规模, 选取为$N_{\rm{pop}} =C_{H+M-1}^{M-1} $.其中, $M$为目标函数个数, $H$取决于决策者.在实验中, 目标个数分别为2和3时, 选取$H$为99和23, 确定种群规模$N_{\rm pop}$为100和300.另外, 选取权向量的邻域个数为20.

2) 进化算子中的控制参数:差分进化操作中, 交叉率CR $=1.0$, 尺度因子F $=0.5$; 高斯变异概率$P_m$ = $1/N$, $N$为决策变量的维数; PBI分解方法中的惩罚参数$\theta =5$.

3) 算法运行次数与停止准则:针对任一测试问题, 各个算法独立运行15次; 最大进化代数$\tau=$ $3\, 000$; 动态环境变化频率$\tau_d =30$; 动态环境变化强度$n_d =10$; 动态环境变化次数$K=100$.

4) 动态鲁棒优化模型中, 目标函数值的预测采用MA方法实现, 其构造的时间序列长度$m=10$.

3.4   实验结果与分析

考虑到目前仅有文献[16]与本文所提动态多目标鲁棒优化方法相关, 因此, 仅对上述二者进行性能对比分析.根据文献[16]的鲁棒Pareto解集定义, 采用MOEA/D方法求解动态多目标优化问题, 获得其在100个动态环境下的鲁棒Pareto解集, 并生成一组真正的鲁棒Pareto解集.在如下实验分析中, 方法1代表文献[16]所提鲁棒优化方法; 方法2为本文所提鲁棒动态优化模型1;方法3为本文所给出的鲁棒动态优化模型2.

3.4.1   鲁棒Pareto解集对参数的敏感性分析

表 3给出了不同稳定性阈值下, 基于两种模型的动态多目标鲁棒优化方法所获得Pareto解集的平均生存时间.当时间窗$T=2$, 稳定性阈值$\eta$分别为0.2, 0.4和0.6时, 所获得鲁棒Pareto解集的平均生存时间, 表明当稳定性阈值越大时, 所得Pareto解集的平均生存时间越大.这是因为阈值反映决策者对平均适应度偏离程度的容忍程度. $\eta$越大, 决策者可以接受未来连续环境下适应值的波动程度越大, Pareto解的时间鲁棒性越好, 但每个动态环境下的解性能变弱.

表 3  不同稳定性阈值下两种模型的平均生存时间比较($T=2$, $n_{ d} =10$, $\tau_{ d} =30$)

Table 3  Average survival time of methods 2 and 3 in different stability thresholds ($T=2$, $n_{ d} =10$, $\tau_{ d} =30$)

方法	参数	函数1	函数2	函数3	函数4	函数5	函数6	函数7	函数8
方法2	$\eta =0.2$	2.4480±0.0314	7.2233±0.0082	2.3433±0.0342	2.6867±0.0049	1.6120±0.0041	9.9680±0.0041	2.2813±0.0106	1.2407±0.0110
	$\eta=0.4$	3.4640±0.0099	21.658±0.0086	3.9428±0.0208	4.2027±0.0096	2.4093±0.0026	30.197±0.5440	3.6000±9.1E-16	1.7000±7.3E-15
	$\eta=0.6$	${\bf 4.5720\pm0.0077}$	${\bf 49.51\pm0} $	$ {\bf 5.9933\pm0.0049} $	${\bf 6.1293\pm0.0080}$	$ {\bf 4.3727\pm0.0059} $	$ {\bf 49.510\pm0} $	$ {\bf 4.8247\pm0.0210} $	$ {\bf 2.2867\pm0.0105}$
方法3	$\eta=0.2$	2.5067±0.0315	7.4120±0.2402	2.5500±0.0746	2.7800±9.1E-16	1.7200±2.2E-16	9.9453±0.0074	2.3347±0.0106	1.2600±2.2E-16
	$\eta=0.4$	3.5133±0.0082	20.941±0.0083	4.3147±0.0280	4.3187±0.0106	2.5500±0.0318	30.007±0.4951	3.6520±0.0068	1.7800±4.6E-16
	$\eta=0.6$	${\bf 4.5940\pm0.0019}$	$ {\bf 49.510\pm0} $	$ {\bf 6.1227\pm0.0046} $	$ {\bf 6.2727\pm0.0110} $	$ {\bf 4.4787\pm0.0035} $	$ {\bf 49.510\pm0} $	$ {\bf 4.8607\pm0.0026} $	$ {\bf 2.4033\pm0.0082}$

下载: 导出CSV 

| 显示表格

表 4是在不同时间窗$T$=2, 4, 6下, 本文所提两种模型所获得鲁棒Pareto最优解集的平均生存时间随着时间窗的增大, 解集的平均生存时间呈现出一定的减小趋势.其原因在于, 随着时间窗的增大, 满足可行约束条件的Pareto解在不同动态环境下的目标值波动增强, 解的可适用性降低, 时间鲁棒性变差.

表 4  不同时间窗下两种模型的平均生存时间比较($\eta =0.4$, $n_{ d} =10$, $\tau_{ d} =30$)

Table 4  Average survival time of methods 2 and 3 in different time windows ($\eta =0.4$, $n_{ d} =10$, $\tau_{ d} =30$)

方法	参数	函数1	函数2	函数3	函数4	函数5	函数6	函数7	函数8
方法2	$T=2 $	$ {\bf 3.4640\pm0.0099} $	$ {\bf 21.658\pm0.0086} $	3.9428±0.0208	$ {\bf 4.2027\pm0.0096}$	2.4093±0.0026	$ {\bf 30.197\pm0.5440} $	$ {\bf 3.6000\pm9.1}$E-${\bf 16}$	$ {\bf 1.7000\pm7.3}$E-${\bf 15}$
	$T=4 $	3.1460±0.0051	21.606±0.1788	3.6447±0.0168	3.9440±0.0203	$ {\bf 2.4500\pm0.0038} $	29.448±0.7324	3.1707±0.0406	1.5220±0.0056
	$T=6 $	$ 3.1180\pm0.0115 $	$ 21.650\pm3.6$E-15	$ {\bf 4.0373\pm0.0059} $	$ 3.9007\pm0.0167$	$ {{\bf 2.4500\pm0}} $	29.085±0.6836	3.2027±0.0448	1.5220±0.0056
方法3	$T=2 $	$ {\bf 3.5133\pm0.0082} $	$ {\bf 20.941\pm0.0083} $	$ {\bf 4.3147\pm0.0280} $	$ {\bf 4.3187\pm0.0106}$	2.5500±0.0318	$ {\bf 30.007\pm0.4951} $	$ {\bf 3.6520\pm0.0068} $	$ {\bf 1.7800\pm4.6}$E$-{\bf 16}$
	$T=4 $	3.4033±0.0082	20.900±0	3.9533±0.0082	4.1893±0.0110	2.7100±4.5E-16	28.068±0.2995	3.3927±0.0546	1.5600±2.2E-16
	$T=6 $	3.4120±0.0108	20.644±0.3445	4.1467±0.0374	4.2313±0.0125	$ {\bf 2.8813\pm0.0376} $	28.340±0.4003	3.4820±0.0441	1.5940±0.0356

下载: 导出CSV 

| 显示表格

表 5是在不同时间窗下, 两类算法所获得鲁棒Pareto解集的鲁棒逆向世代距离RIGD的统计结果.该性能指标体现了鲁棒解集的收敛性和分布性. RIGD值越小, 表明解集的收敛性越好.数据对比可见, 当$T$取较小值时, 鲁棒Pareto解的平均鲁棒世代距离较小.这是因为随着时间窗$T$的增大, 期望的时间鲁棒性更好, 但是受到预测性能影响, 鲁棒Pareto解集在多个连续动态环境下的适应值会随着$T$的增加偏离真实PF$(k)$, 从而导致收敛性下降.

表 5  不同时间窗下两种模型的RIGD测度比较($\eta =0.4$, $n_{ d} =10$, $\tau_{ d} =30$)

Table 5  RIGD of methods 2 and 3 in different time windows ($\eta =0.4$, $n_{ d} =10$, $\tau_{ d} =30$)

方法	参数	函数1	函数2	函数3	函数4	函数5	函数6	函数7	函数8
方法2	$T=2$	$ 0.2267\pm0.0092 $	$ {\bf 0.0816\pm0.0555} $	$ {\bf 0.3815\pm0.0230} $	$ 0.4009\pm0.0178$	$ {\bf 0.4638\pm0.0218} $	$ 0.1009\pm0.0754 $	$ 0.3476\pm0.0132 $	$ 2.5449\pm0.1079$
	$T=4$	$ {\bf 0.1944\pm0.0081} $	$ 0.1364\pm0.0911 $	$ 0.4142\pm0.0296 $	$ 0.4009\pm0.0177$	$ 0.5424\pm0.0250 $	$ {\bf 0.0707\pm0.0323} $	$ {\bf 0.3306\pm0.0171} $	$ 2.1396\pm0.0972$
	$T=6$	$ 0.2270\pm0.0077 $	$ 0.1552\pm0.1226 $	$ 0.8548\pm0.1841 $	$ {\bf 0.4008\pm0.0200} $	$ 0.7334\pm0.0571 $	$ 0.0928\pm0.0457 $	$ 0.3476\pm0.0181 $	$ {\bf 2.0528\pm0.0880}$
方法3	$T=2$	$ 0.2276\pm0.0063 $	$ {\bf 0.0849\pm0.0559} $	$ {\bf 0.3842\pm0.0215} $	$ 0.4520\pm0.0169$	$ {\bf 0.5907\pm0.0444} $	$ {\bf 0.0811\pm0.0543} $	$ 0.3477\pm0.0128 $	$ {\bf 2.5629\pm0.0999}$
	$T=4$	$ {\bf 0.1976\pm0.0066} $	$ 0.2238\pm0.1382 $	$ 0.5184\pm0.1738 $	$ {\bf 0.4057\pm0.0070}$	$ 0.7655\pm0.0278 $	$ 0.1168\pm0.1083 $	$ {\bf 0.3328\pm0.0089} $	$ 3.4348\pm0.1043$
	$T=6$	$ 0.5411\pm0.1200 $	$ 0.1432\pm0.0457 $	$ 1.2358\pm0.0851 $	$ 0.6615\pm0.0048 $	$ 1.2572\pm0.1641 $	$ 0.1812\pm0.0857 $	$ 0.9319\pm0.0289 $	$ 11.416\pm0.2446 $

下载: 导出CSV 

| 显示表格

3.4.2   环境变化程度对算法的影响

环境变化程度取决于动态变化频率和动态变化强度两方面.为此, 在固定时间窗$T=2$和稳定性阈值$\eta=0.4$的条件下, 分别针对不同动态变化频率和动态变化强度设计了两类实验.其一, 固定动态变化频率$\tau_{ d}=30$, 分别设置动态变化强度为$n_{ d}=5$, 10, 20, 对比方法2和方法3的实验结果, 如图 5所示; 其二, 固定动态变化强度$n_{ d}=10$, 分别设置动态变化频率为$\tau_{\rm d}=10, 20, 30$, 对比分析不同算法的实验结果, 如图 6所示.

图 5  不同变化强度$n_{ d} =5, 10, 20$下方法2和方法3的平均生存时间比较($\eta=0.4$, $T=2$, $\tau _{ d} =30$)

Fig. 5  Average survival time of methods 2 and 3 in different change severity $n_{ d} =5, 10, 20$ ($\eta=0.4$, $T=2$, $\tau _{ d} =30$)

下载: 全尺寸图片 幻灯片

图 6  不同变化频率$\tau _{ d}=10, 20, 30$下方法2和方法3的平均生存时间比较($\eta =0.4$, $T=2$, $n_{ d} =10$)

Fig. 6  Average survival time of methods 2 and 3 in different change frequency $\tau _{ d}=10, 20, 30$ ($\eta =0.4$, $T=2$, )

下载: 全尺寸图片 幻灯片

1) 图 5给出了不同动态变化强度下, 基于方法2和方法3的动态多目标鲁棒优化方法所获得鲁棒Pareto解集的平均生存时间黑箱图.由图可知, 随着参数$n_{ d}$增大, 获得的鲁棒Pareto解集平均生存时间越大, 并且方法3所得解集的平均生存时间略优于方法2.这是因为, 当参数$n_{ d}$增大时, 相邻动态环境下真实Pareto解的动态变化程度减弱, 在一定鲁棒阈值条件下, 每个鲁棒Pareto解的生存时间将会增加, 从而导致鲁棒Pareto解集的平均生存时间越大.

2) 动态变化频率$\tau_{ d}$的取值越小, 代表环境动态变化越频繁. 图 6针对本文所提两种算法, 在不同的动态变化频率$\tau_{ d}=10, 20, 30$下, 对比8个测试函数的鲁棒Pareto解集所包含的鲁棒Pareto解个数.从图中可以看出, 动态变化频率对鲁棒Pareto解集规模的影响较小, 表明其解集的平均生存时间波动较小, 时间鲁棒性较稳定.

3.4.3   不同算法及模型下的性能比较

针对本文所提出的两种鲁棒优化模型, 考虑在时间窗$T=2$, 且稳定性阈值$\eta=0.4$时, 通过基于分解的鲁棒进化优化方法获得不同动态环境下的鲁棒Pareto解.比较上述两种方法获得的鲁棒Pareto解集的平均生存时间如表 6、表 7和图 7所示.

表 6  三种方法获得鲁棒Pareto解集的平均生存时间比较($\eta =0.4$, $T=2$, $n_{ d}=10$, $\tau _{ d} =30$)

Table 6  Average survival time of robust Pareto solutions on three methods ($\eta =0.4$, $T=2$, $n_{ d}=10$, $\tau _{ d} =30$)

方法	函数1	函数2	函数3	函数4	函数5	函数6	函数7	函数8
方法1	2.5800±0.0086	10.090±0.1376	2.6200±0.3880	2.7300±0.6077	2.9900±0.0560	17.940±0.8734	2.2400±0.3525	1.3500±3.6E-15
方法2	3.4640±0.0099	$ {\bf 21.658\pm0.0086} $	3.9428±0.0208	4.2027±0.0096	$2.4093\pm0.0026 $	$ {\bf 30.197\pm0.5440} $	3.6000±9.1E-16	1.7000±7.3E-15
方法3	${\bf 4.5073\pm0.0133} $	20.941±0.0083	$ {\bf 4.3147\pm0.0280} $	$ {\bf 4.3187\pm0.0106}$	${\bf 2.5500\pm0.0318} $	30.007±0.4951	$ {\bf 3.6520\pm0.0068} $	$ {\bf 1.7800\pm4.6}$E-${\bf 16}$

下载: 导出CSV 

| 显示表格

表 7  三种方法获得鲁棒Pareto解集的RIGD测度比较($\eta =0.4$, $T=2$, $n_{ d}=10$, $\tau _{ d} =30$)

Table 7  RIGD of robust Pareto solutions on three methods ($\eta =0.4$, $T=2$, $n_{ d}=10$, $\tau _{ d} =30$)

方法	函数1	函数2	函数3	函数4	函数5	函数6	函数7	函数8
方法1	${\bf 0.0095\pm0.0021} $	0.2438±0.0316	$ {\bf 0.1544\pm0.0579} $	$ {\bf 0.0543\pm0.0043}$	${\bf 0.3619\pm0.1683} $	0.1564±1.04E-7	$ {\bf 0.1335\pm0.0743} $	$ {\bf 0.1866\pm1.40}$E-${\bf 6}$
方法2	0.2267±0.0092	$ {\bf 0.0816\pm0.0555} $	0.3815±0.0230	0.4009±0.0178	0.4638±0.0218	0.1009±0.0754	0.3476±0.0132	2.5449±0.1079
方法3	0.2276±0.0063	0.0849±0.0559	0.3842±0.0215	0.4520±0.0169	0.5907±0.0444	$ {\bf 0.0811\pm0.0543} $	0.3477±0.0128	2.5629±0.0999

下载: 导出CSV 

| 显示表格

图 7  100个动态环境下三种方法的鲁棒生存时间比较$\eta=0.4$, $T=2$, $n_{ d} =10$, $\tau _{ d} =30$

Fig. 7  Robust survival time of three methods in 100 environments ($\eta =0.4$, $T=2$, $n_{ d} =10$, $\tau _{ d} =30$)

下载: 全尺寸图片 幻灯片

1) 如表 6和图 7所示, 三种方法所得鲁棒Pareto最优解的生存时间都大于等于1, 代表其适用于至少一个或多个动态环境.此外, 采用本文所提方法2和方法3获得的鲁棒Pareto最优解集的平均生存时间明显优于方法1所得Pareto解的平均生存时间, 表明本文所提方法较前者更加合理, 获得的Pareto最优解集具有更强的时间鲁棒性.对比方法2和方法3所得鲁棒Pareto最优解的平均生存时间, 发现方法3所获得解的平均生存时间略优于方法2, 表明适应度满意度这一约束条件对时间鲁棒性改善具有一定促进作用.

2) 图 8给出了三种方法所获得鲁棒Pareto最优解集中, 平均鲁棒生存时间分别为$L_i =1$和$L_i$ = 2的解占解集规模的比例, 记为Per 1和Per 2.三种方法均选取参数$\eta =0.4$, $T=2$.可见, 方法2和方法3的Per 1明显小于方法1的Per 1, 而方法2和方法3的Per 2明显大于方法1的Per 2.表明本文所给的动态多目标鲁棒优化方法, 可以有效地提高动态鲁棒Pareto解集的时间鲁棒性能.

图 8  3种方法获得鲁棒Pareto解集的生存时间百分比比较($\eta =0.4$, $T=2$, $n_{ d}=10$, $\tau _{ d} =30$)

Fig. 8  Robust survival time percentage of three methods ($\eta =0.4$, $T=2$, $n_{ d}=10$, $\tau _{ d} =30$)

下载: 全尺寸图片 幻灯片

3) 表 7统计了三种方法在阈值$\eta=0.4$, $T=2$下的鲁棒Pareto解集的鲁棒逆向世代距离.可以看出, 本文所提两种方法的RIGD值大多数不如方法1, 主要原因是方法2和方法3所得鲁棒Pareto解集有更好的鲁棒性, 鲁棒性和收敛性在一定程度上存在冲突.

4) 分析8个测试函数所得鲁棒Pareto最优解集的平均生存时间, 发现函数2和函数6的结果明显大于其他测试函数.进一步分析可知, 这两个函数的PS$(k)$不随动态环境变化而变化, 但PF$(k)$却跟随动态环境变化, 并且变化不显著.因此, 在进化过程中很容易找到Pareto最优解, 且该解对未来连续动态环境的适应性较强, 可以被应用于更多动态变化环境, 导致其动态Pareto解集具有较好的鲁棒性.

5) 统一设定问题的动态变化参数为$\eta=0.4$, $T$ = 2, $n_{ d}=10$, $\tau_{ d}=30$, 且发生100次环境变化.实验硬件为IBM X3100M5 E3 3.10 GHz CPU, 8 GB内存的电脑, 软件采用MATLAB2011b.表 8统计了三种不同算法获得鲁棒Pareto解集所消耗的运行时间(s).比较可得, 采用本文所提两种模型, 可以在相同条件下, 以较少用时得到一组有效的鲁棒Pareto解集.

表 8  三种方法求解鲁棒Pareto最优解的平均运行时间(s)比较($\eta=0.4$, $T=2$, $n_{ d} =10$, $\tau _{ d} =30$)

Table 8  Average elapsed time of robust Pareto solutions on three methods ($\eta=0.4$, $T=2$, $n_{ d} =10$, $\tau _{ d} =30$)

方法	函数1	函数2	函数3	函数4	函数5	函数6	函数7	函数8
方法1	51.3681±3.4762	13.2158±4.0439	70.8508±1.2940	185.9928±10.9027	710.6224±70.2617	1.7005±0.0077	54.8422±2.6081	128.4283±2.3889
方法2	${\bf 46.0068\pm1.4818} $	$ {\bf 7.0358\pm4.2498} $	59.0582±1.1731	$ {\bf 172.3075\pm7.6934}$	${\bf 315.2275\pm8.4846} $	$ {\bf 0.8673\pm0.0140} $	$ {\bf 45.7206\pm2.3813} $	89.7115±0.9206
方法3	49.3816±1.8679	10.2509±3.2404	$ {\bf 55.8101\pm2.2671} $	190.0556±19.8364	340.7134±21.4715	1.0038±0.0046	47.8325±2.0709	$ {\bf 84.1739\pm0.8742}$

下载: 导出CSV 

| 显示表格

图 9进一步针对$\eta =0.4$, $T=2$, $n_{ d} =10$, $\tau _{ d}$ $=$ 100动态环境, 给出了8个测试函数的前4个鲁棒Pareto解与其所适用的$L_i$个连续动态环境下真实Pareto解的适应度值比较.每个图中, 虚线为相应方法所得鲁棒Pareto最优前沿RPF$(i)$, 实线为它所适用的$L_i$个相邻动态环境下的真实Pareto前沿PF$(i)$.显然, 鲁棒Pareto解能够以较好的收敛性逼近相应动态环境下的真实Pareto解.

图 9  8个测试函数的前四个鲁棒Pareto解与其所适应的动态环境下的真实Pareto解的比较图\\ ($\eta=0.4$, $T=2$, $n_d =10$, $\tau _d =100$)

Fig. 9  The first four robust Pareto fronts and true Pareto fronts in their adaptive environments ($\eta=0.4$, $T=2$, $n_d =10$, $\tau _d =100$)

下载: 全尺寸图片 幻灯片

4.   结论

本文提出了一类寻找动态鲁棒Pareto最优解集的进化优化方法.对于动态多目标优化问题, 动态鲁棒Pareto解集是指在某一时刻下的Pareto较优解可以以一定稳定性阈值逼近未来多个连续动态环境下的真实前沿, 从而直接作为这些环境下的Pareto解集, 以减小计算代价.首先, 分别通过构造平均生存时间和平均适应度函数来描述解的时间鲁棒性能和性能鲁棒性.其次, 将上述鲁棒性能评价转化为两类鲁棒优化模型, 采用基于分解的多目标进化优化方法和无惩罚约束处理方法, 构建了动态多目标分解鲁棒进化优化方法.仿真实验结果表明, 两种模型所获得的动态鲁棒Pareto最优解, 在满足适应度满意阈值条件下具有更好的时间鲁棒性.另外, 对于不同时间窗$T$和稳定性阈值$\eta$下的解性能分析表明, 较大的$\eta$和较小的时间窗有利于鲁棒Pareto解集性能.如何针对不同动态多目标优化问题选取合理的参数是未来有待深入的问题.此外, 对未来动态环境下的适应度预测在算法中具有显著影响, 因此预测方法的研究也是我们今后的研究工作.