现实世界中的许多复杂优化问题多为动态的, 会随时间发生变化. 例如, 在调度问题中, 新任务不断到达, 产品原材料质量不断变化, 机器随时间增加不断磨损等^[1]. 演化算法是求解复杂优化问题效果显著的一种智能算法, 在求解优化问题时无需考虑函数本身的可微性、可导性、连续性等信息, 具有较强的自适应性和自组织性^[2], 已在求解静态优化问题上取得了很好的效果, 因此许多研究者致力于在演化算法中寻求解决动态优化问题的方法改进和创新. 动态优化问题的最优解是不断变化的, 求解动态优化问题的算法要能够在动态环境中及时跟踪到最优解^{[3, 4]}. 因此, 对动态环境进行检测与响应已成为动态优化问题的研究热点.

为了跟踪随环境变化而变化的最优解, 首先要及时准确地检测环境的变化. 其次在环境变化后, 优化算法要及时采取响应策略紧密地跟踪解的变化以获得最优解. 许多研究者对这两方面进行了研究, 主要工作有: Hu和Eberhart ^[5]通过监测全局最优解和次最优解的适应值变化来估计环境的变化, Carisle ^[6]提出随机选择"Sentry"粒子作为检测环境变化的感应粒子, 单世民和邓贵仕^[7]根据个体极值的适应值变化情况来判断环境是否变化. 在响应策略上, 他们选择重新初始化种群或部分个体. 姜立强等^[8]将种群分为跟踪和搜索两个种群, 也是通过监测跟踪种群的当前最优解和次优解来判断环境是否发生变化. 当发现环境发生变化时, 重新初始化搜索种群, 并对两个种群采取不同的变异策略, 提高搜索新最优解的效率. 胡静等^[9-10]通过监测相近粒子前后两代适应值与距离比值的变化率来检测环境变化, 并构造种群多样性函数, 通过粒子的扩散运动响应环境变化. 唐剑等^[11]通过激活种群中的停滞粒子适应环境变化. 赵传信等^[12]通过监测全局最优个体和距离其最远的个体的适应值变化来检测环境变化, 环境变化后进行以种群为中心的扩散运动. 成照乾等^[13]根据均匀分布的静态粒子的适应值变化来判断环境是否变化. 当静态粒子检测到某个区域环境变化后, 每个粒子都朝该区域运动以响应环境变化. 朱庆保等^[14]设置多组粒子进行搜索, 通过监测每组最优解的适应值变化来检测环境变化. 环境变化后, 令多样性较差的一组粒子向最优粒子周围扩散. 陈健等^[15]和袁亦川等^[16]都专门设置一个种群来检测环境变化, 种群中至少有一个个体的适应值变化时, 认为环境发生了变化. 环境变化后, 初始化种群以适应新的环境.

动态优化问题环境检测的实质是检测最优解是否变化及如何变化, 目前环境检测方法主要是通过监测某些点的适应值变化情况来检测环境的变化. 在已有检测方法中, 如果所选的监测点与最优解差异较大, 则会降低环境检测的准确性. 如: 最优解发生了较明显的变化, 但这些点的适应值变化却相对不明显; 或最优解没有发生明显变化, 而这些点的适应值的变化却相对较明显. 因此, 在环境检测方法设计中, 如何设置监测点十分重要. 在响应策略方面, 研究者们采用重新初始化种群或部分个体、扩散种群、激活粒子等方法. 这些响应方法在一定程度上提高了算法的空间搜索能力. 但是, 重新初始化方法没有很好地利用之前的信息, 而且不断地初始化也会降低算法性能. 此外, 没有检测到环境变化后的最优解所在区域就扩散种群或者激活粒子的做法具有一定的盲目性和随机性, 并不能快速准确地跟踪到变化后的最优解.

正交试验设计是一种多因素多水平的优化试验设计方法, 它从全面试验的样本点中挑选出部分有代表性的样本点进行试验, 只用较少的试验次数就能找出因素水平间的最优搭配^[17]. 一些研究者将正交试验设计思想引入到优化算法的交叉步骤中以获得更好的子代, 例如, 李志勇等^[18]将正交试验设计方法用到多目标优化算法中, 设计了正交交叉算子. 具体做法为: 对父代的两个个体进行量化和分解操作, 使其适应于四因素三水平的正交表, 通过正交试验产生子代. 此外, 正交试验设计也是确定输入变量(因素)关键区域的有效方法, 张晓丽等^[19]利用正交表区间收缩法求函数极值, 其方法大致为: 将自变量作为因子, 将自变量的取值区间二等分作为两个水平, 取区间中点为试验点, 截除实验结果“差”的试验点所在区间的一部分, 逐步缩小自变量的取值区间, 最终得到函数极值. 但是选取区间中点为试验点可能会导致极值点所在的区间被舍去, 进而导致结果不准确(见本文3.3节例子).

基于以上分析, 本文首先对文献[19]中提出的正交表区间收缩方法进行改进. 具体的, 将自变量每一维的取值区间等分为多个子区间, 在每个子区间中随机选取试验点, 分析结果; 重复多次, 挑选出现率高的子区间作为函数极值区间. 实验表明: 改进的方法能有效提高正交试验结果的准确性.

在环境检测方法的改进上, 本文在适应值备选函数的最优解区域内、外, 按一定比例选取勘探粒子, 并根据勘探粒子的适应值变化情况来检测环境变化. 而在环境响应方面, 针对群智能算法, 如粒子群算法、差分进化算法等, 引入记忆机制存储环境信息, 环境变化后根据环境信息和最优解变化情况, 在环境变化后的最优解区域内、外, 按一定比例选取个体以适应新的环境. 算法测试结果证明: 本文提出的环境检测与响应方法能够有效监测最优解的变化情况, 在环境变化后能够及时准确地跟踪变化后的最优解. 改进后的方法充分利用了过去的环境信息, 提高了算法快速收敛到最优解的能力.

在实际中, 动态优化问题的环境可能受到某些随机因素的影响, 从而导致动态优化问题对应的目标函数$ f(x, t) $ (环境)具有某种随机性. 从已查到的文献看, 针对随机目标函数动态优化问题(以下简称随机动态优化问题)的研究工作较少. 陈莉在文献[20]中研究了一类具有马氏性的随机动态优化问题, 其特点是: 已知环境发生变化的时刻, 并假设环境状态的改变规律可用有限状态的马氏链来描述. 基于马氏链预测和记忆机制, 该文设计了相应的遗传算法跟踪优化问题的最优解. 本文将考虑一类新的随机动态优化问题, 其目标函数$ f(x, t) $的最优解因受随机扰动而发生随机偏移. 本文将在随机扰动服从正态分布的假设下, 给出扰动前后环境检测无变化的随机扰动标准差取值上限, 并验证本文改进的环境检测与响应方法处理此类随机动态优化问题的有效性. 一般动态优化问题及两类随机动态优化问题的数学定义见本文2.1节的描述.

本文其余部分结构如下: 第2章介绍随机动态优化问题和正交试验设计; 第3章给出改进后的函数最优解区间收缩方法一般步骤, 进行模拟实验证明该方法的可行性, 并与其它方法比较; 第4章基于函数最优解区间收缩方法, 提出动态优化问题的环境检测与响应方法, 通过粒子群算法进行测试并与其它算法比较, 验证了方法的优良性; 第5章针对受随机扰动的动态优化问题, 给出扰动前后环境检测无变化对应的随机扰动标准差上限, 并检验本文环境检测与响应方法处理这类随机动态优化问题的有效性; 第6章为总结与展望.

本节介绍随机动态优化问题和正交试验设计.

动态优化问题(Dynamic Optimization Problem, DOP)是指目标函数、问题实例或限制条件随时间发生变化的优化问题^[21]. 根据目标函数的个数, 动态优化问题分为单目标动态优化问题和多目标动态优化问题. 当目标函数的个数为1时, 称为单目标动态优化问题, 根据Cruz ^[22]的介绍, 单目标动态优化问题定义如下

(2.1) $ \begin{equation} DOP=\left\{\begin{array}{ll} {\rm optimize} \ f(x, t) \\ {\rm s.t.} \ x\in F(t)\subseteq S, t\in T. \end{array}\right. \end{equation} $

其中$ S $表示搜索空间, $ t $表示时间, $ F(t) $表示$ t $时刻在搜索空间中可行解区域, $ f:S\times T\to R $, $ f $表示目标函数, $ f(x, t)\in R $, $ f(x, t) $表示在$ t $时刻可行解$ x $的目标函数值.

给定备选函数空间$ {\Bbb G}_{1} = \left \{ g_{1} (x), g_{2}(x), \cdots , g_{K} (x) \right \}, x\in S\subset R^{D} $, 则一般动态优化问题的适应值函数为

(2.2) $ \begin{equation} f(x, t)=\left\{\begin{array}{ll}g_{j_{1} } (x), t_{0} \le t< t_{1}, \\g_{j_{2}} (x), t_{1} \le t< t_{2}, \\g_{j_{3}} (x), t_{2} \le t< t_{3}, \\\cdots \\ g_{j_{n}} (x), t_{n-1} \le t\le t_{n}, \\\cdots. \end{array}\right. \end{equation} $

其中$ j_{1}, j_{2}, \cdots , j_{n}, \cdots $均取值于$ \{1, 2, \cdots, K\} $且满足: $ j_{k}\not= j_{k+1}, \forall k\geq 1 $.

随机动态优化问题没有统一定义, 依据随机性来源而定. 下面先给出文献[20]具有马氏性随机动态优化问题的数学定义.

假设已知动态优化问题的环境变化时刻为$ t_0, t_1, \cdots, t_n, \cdots $, 并给定概率分布$ \pi=(\pi_1, $$ \cdots, \pi_K) $与转移概率矩阵$ P=(p_{ij})_{1\leq i, j\leq K} $, 文献[20]具有马氏性的随机动态优化问题则可以描述为

(2.3) $ \begin{equation} f(x, t)=\left\{\begin{array}{ll} Z_1(x), t_{0} \le t< t_{1}, \\Z_2 (x), t_{1} \le t< t_{2}, \\Z_3 (x), t_{2} \le t< t_{3}, \\\cdots \\ Z_{n} (x), t_{n-1} \le t\le t_{n}, \\\cdots. \end{array}\right. \end{equation} $

其中, $ Z_1 $依概率分布$ \pi $取值于备选函数空间$ {\Bbb G}_{1} $; $ Z_2 $依概率分布$ \pi P $取值于备选函数空间$ {\Bbb G}_{1} $; 依此类推, 对任意的$ n> 2 $, $ Z_n $依概率分布$ \pi P^{n-1} $取值于备选函数空间$ {\Bbb G}_{1} $. 此外, 函数取值变化规律满足如下的马氏性

$ P\{Z_{k+1}=g_j|Z_1=g_{j_1}, Z_2=g_{j_2}, \cdots, Z_{k}=g_i\}=P\{Z_{k+1}=g_j| Z_{k}=g_i\}=p_{ij}. $

对于最优解受随机扰动发生随机偏移的随机动态优化问题, 其目标函数需要加上随机扰动. 设$ \varepsilon _{1} , \varepsilon _{2}, \cdots , \varepsilon_{K} $为$ K $个独立同分布的随机向量, 其共同分布为$ N(0, \sigma ^{2} I_{D\times D}) $, 其中$ I_{D\times D} $为$ D $阶单位矩阵. 扰动后的备选函数空间为

$ {\Bbb G}_{2} = \{ g_{1} (x+\varepsilon _{1}), g_{2}(x+\varepsilon _{2}), \cdots , g_{K} (x+\varepsilon_{K}) \}, x\in S\subset R^{D}, $

则本文受随机扰动的动态优化问题可描述为

(2.4) $ \begin{equation} f(x, t) =\left\{\begin{array}{ll}g_{j_{1} } (x+\varepsilon _{1}), t_{0} \le t< t_{1}, \\g_{j_{2}} (x+\varepsilon _{2}), t_{1} \le t< t_{2}, \\g_{j_{3}} (x+\varepsilon _{3}), t_{2} \le t< t_{3}, \\\cdots \\g_{j_{n}} (x+\varepsilon _{n}), t_{n-1} \le t\le t_{n}, \\\cdots. \end{array}\right. \end{equation} $

其中$ j_{1}, j_{2}, \cdots , j_{n}, \cdots $均取值于$ \{1, 2, \cdots, K\} $且满足: $ j_{k}\not= j_{k+1}, \forall k\geq 1 $.

在试验设计中, 假设需要考虑$ q $个因素, 每个因素有$ a $个水平, 若采取全面试验方案, 则需要做$ a^{q} $次试验, 当因素个数和水平数都比较大时, 显然全面试验较难实现, 此时希望用较少的试验次数找出因素水平间的最优搭配, 而正交试验设计正是解决这种问题的常用方法. 正交试验设计是使用按正交性排列好的正交表来安排试验的, 对于$ q $个因素、$ a $个水平的试验, 选择$ L_{b} (a^{q} ) $正交表, 其中$ b $表示试验次数, 一般来说$ b $远小于$ a^{q} $. 以四因素三水平的试验为例, 若依据正交表$ L_{9} (3^{4} ) $来安排试验, 则只需要做9次试验; 若实施全面试验, 则需要做$ 3^{4} =81 $次试验. 由此可见正交试验设计能够大幅提高效率, 对于因素和水平较多的试验, 该方法的优越性将更加明显.

目前动态优化问题的环境检测方法主要是通过监测某些点的适应值变化情况来检测环境的变化. 设置监测点的方法主要分为两种, 一是监测个体最优解、全局最优解等特殊点, 二是在整个搜索空间中随机设置一些监测点. 这两种方法都有一定的盲目性, 如果监测点与最优解的差异较大, 那么环境检测的准确性将大大降低. 因此要在最优解周围设置监测点, 需要先找出最优解所在的较小区域.

正交试验设计是确定输入变量(因素)关键区域的有效方法. 求光滑的多维函数在定义域内的最优解问题本质上是自变量在每一维上的分量与每一小段取值区间的组合问题, 因此可将其转化为统计问题, 基于正交试验设计思想寻找函数最优解所在的较小区域. 本节对文献[19]中的正交表区间收缩方法进行改进, 构造函数最优解区间收缩方法, 通过模拟实验来验证利用该方法确定动态优化问题最优解所在区域的可行性, 并与文献[19]中方法比较.

问题描述: $ f(x) $为某一$ D $元函数, 自变量$ x=(x_{1} , x_{2}, \cdots , x_{D} ) $, $ x_{i} \in\left [ a_{i} , b_{i}\right ], $$ i=1, 2, \cdots , D $. $ \exists x^{*} =(x_{1}^{*} , x_{2}^{*} , \cdots , x_{D}^{*}) $为函数$ f(x) $的最小值点, 即满足: $ f( x_{1}^{*} , x_{2}^{*} , \cdots , x_{D}^{*} ) $$ = \min\limits _{x_{i}\in [a_{i} , b_{i}]} f(x_{1} , x_{2} , \cdots , x_{D}) $. 求最优解$ x^{*} $所在的较小区域.

文献[19]正交表区间收缩方法过程如下.

将自变量作为因子, 将自变量的取值区间二等分作为两个水平, 取区间中点为试验点, 选择正交表安排试验, 并计算试验指标值. 若某个因子取第一个水平时, 对应的试验指标最大, 则该因子从左端去掉$ 1/4 $取值区间; 若某个因子取第二个水平时, 对应的试验指标最大, 则该因子从右端去掉$ 1/4 $取值区间. 重复此过程, 逐步缩小各个因子的取值区间.

本文改进后的函数最优解区间收缩方法步骤如下

(1) 问题转化, 把求$ f(x) $的最小值问题转化为多因子多水平的优化问题. $ x_{1} , x_{2}, x_{3}, \cdots $, $ x_{D} $为$ f(x) $的$ D $个影响因子, 把每个因子$ x_{i} $的取值区间$ \left [ a_{i} , b_{i}\right ] $平均划分为$ r $个子区间, 作为$ x_{i} $的$ r $个水平.

(2) 根据因子和水平个数选择合适的正交表, 在每个等长子区间内随机选取一点作为试验点, 根据正交表安排试验, 计算正交表上每个试验对应的指标值.

(3) 采用极差分析法分析结果, 找出函数$ f(x) $取得最小值时各因子$ x_{1} , x_{2}, \cdots , x_{D} $所在的子区间. 重复随机选取试验点做$ M(M>2) $轮正交试验, 统计结果, 找出各因子取值区间中出现次数最多的子区间, 作为各因子新的取值区间, 这些区间构成了该函数最优解的所在区域. 若区域较大, 将每个因子的取值区间再次划分重复上述过程即可缩小区域.

本节选用最优解已知的单峰函数和多峰函数进行模拟实验, 检验利用函数最优解区间收缩方法确定函数最优解所在区域的可行性.

实验1  选取单峰三元函数$ f(x_{1} , x_{2} , x_{3} )=(x_{1}-1.3)^{2} +(x_{2}-7.8)^{2}+(x_{3}-3.5)^{2}+1, $$ x_{i} \in [0, 9], i=1, 2, 3 $. 显然, 当$ x_{1}=1.3, x_{2}=7.8, x_{3}=3.5 $时, 此函数取得最小值1. 下面利用函数最优解区间收缩方法找出此函数最优解所在的较小区域.

(1) 将各因子$ x_{i} $的取值区间$ \left [ 0, 9 \right ] $平均划分为三个子区间: $ \left [ 0, 3 \right ), \left [ 3, 6 \right ), \left [ 6, 9 \right ] $, 作为$ x_{i} $的三个水平, 选择$ L_{9} (3^{4} ) $正交表, 将$ x_{1}, x_{2}, x_{3} $放在第1, 2, 3列上, 在$ \left [ 0, 3 \right ), \left [ 3, 6 \right ), \left [ 6, 9 \right ] $上分别随机选取一点作为试验点, 根据正交表安排试验, 计算每个试验对应的指标值.

(2) 采用极差分析法对结果进行分析, 找出三元函数$ f(x_{1} , x_{2} , x_{3} ) $取得最小值时各因子$ x_{1} , x_{2}, x_{3} $的取值子区间. 重复随机选取试验点做100轮正交试验, 统计结果, 找出$ x_{1} , x_{2}, x_{3} $各自取值区间中出现次数最多的子区间, 分别为$ \left [ 0, 3 \right ) $, $ \left [ 6, 9 \right ] $, $ \left [ 3, 6 \right ) $, 作为$ x_{1} , x_{2}, x_{3} $新的取值区间.

(3) 把$ \left [ 0, 3 \right ) $, $ \left [ 6, 9 \right ] $, $ \left [ 3, 6 \right ) $都再次平分为三个小区间, 多次重复上述过程缩小各因子的取值区间, 最后得出三元函数$ f(x_{1} , x_{2} , x_{3} ) $取得最小值时$ x_{1} , x_{2}, x_{3} $较小的取值区间分别为$ \left [ 1, 4/3 \right ) $, $ \left [ 23/3, 8 \right ] $, $ \left [ 10/3, 11/3 \right ) $, 这些区间构成了三元函数最优解所在的较小区域, 显然函数最优解$ (1.3, 7.8, 3.5) $在这个区域内.

实验2  取动态优化问题中常用的多峰Rastrigin函数来检验函数最优解区间收缩方法对于多峰函数的有效性. $ f(x)=\sum\limits_{i=1}^{D} [x_{i}^{2} -10\cos(2\pi x_{i})+10], x_{i} \in [-5, 5] $, 它的三维图如图 1所示, 由图 1可知, 此函数在定义域内有多个局部极小值, 并已知它在$ (0, 0, \cdots , 0) $处取得最小值0.

维数设置为3, 将$ x_{1} , x_{2} , x_{3} $看作多峰函数$ f(x) $的三个影响因子, 仿照实验1的步骤不断收缩各因子的取值区间, 结果为多峰函数$ f(x) $取得最小值时$ x_{1} , x_{2}, x_{3} $较小的取值区间都为$ \left [ -5/27, 5/27 \right ] $, 显然最优解$ (0, 0, 0) $在这个区域内.

对于单峰多元函数来说, 它在定义域内只有一个极小值点, 利用函数最优解区间收缩方法不断缩小各因子的取值区间, 最后能得到函数全局最小值点所在的较小区域. 对于多峰函数来说, 它在定义域内存在多个极小值点, 从逻辑上来说, 利用函数最优解区间收缩方法可能得到的是函数某个极小值点所在的区域. 但是该方法在每个区间选取试验点时是随机选取的, 而且重复100次实验, 如果函数足够光滑, 那么就能保证全局最小值点所在的区间不会被舍去, 因此最后得到的是函数全局最小值点所在的较小区域, 实验2的结果证明了该方法的准确性. 因此, 对于光滑的多维单峰和多峰函数在给定区域的最值问题, 可以利用函数最优解区间收缩方法寻找函数最优解所在的较小区域.

本节将分别利用文献[19]及本文改进后的函数最优解区间收缩方法, 求$ D=3 $时函数$ f(x_{1} , x_{2} , \dots, x_{D} )=\sum\limits_{i=1}^{D} [x_{i}\sin(x_{i})\cos(2x_{i})-2x_{i}\sin(3x_{i})] $, $ x_{i}\in [3, 7] $最小值点所在的较小区域. 当$ D=1 $时, 函数图像如图 2所示, 显然此时函数的最小值点位于区间$ \left [ 6.5, 7 \right ] $上. 因此, 当$ x_{i}\in \left [ 6.5, 7 \right ], i=1, 2, 3 $时, 函数$ f(x_{1} , x_{2} , x_{3} ) $取得最小值.

根据函数最优解区间收缩方法步骤, 选用$ L_{9} (3^{4}) $正交表, 设置$ r=3 $, $ M=100 $, 经过两次收缩, 结果为函数取得最小值时, $ x_{1} , x_{2} , x_{3} $的取值区间都为$ (59/9, 7] $, 与实际结果相符. 而按照文献[19]中的正交表区间收缩法步骤, 选用$ L_{4} (2^{3}) $正交表安排实验, 经过若干次收缩, $ x_{1} , x_{2} , x_{3} $的$ \left [ 6, 7 \right ] $部分取值区间都被舍去, 而函数的最小值点位于此区间. 因此, 文献$ [19] $中由正交表区间收缩法求得的最小值点所在区域不准确.

比较两种方法可知, 文献[19]中正交表区间收缩法将自变量原来的取值区间两等分, 并分别选取区间中点为试验点, 最小值点位于区间$ [5, 7] $上, 而该区间的中点6处对应的试验指标却最大, 从右端去掉$ 1/4 $区间, 导致最小值点被舍去, 最终结果不准确. 而本文的函数最优解区间收缩方法在取值区间内随机选取一个试验点, 重复$ M $轮, 统计每轮实验得到的最小值点所在的子区间, 出现次数最多的子区间为各个因子新的取值区间. 如果$ M $取值足够大(本文中$ M $取100), 那么就能以较大的概率保证最小值点所在的区间不会被舍去. 因此, 本文的函数最优解区间收缩方法求解函数最优解所在区域的准确性更高.

本节基于函数最优解区间收缩方法, 改进现有动态优化问题的环境检测与响应策略, 利用粒子群算法进行测试并与其它算法比较.

为了提高算法对动态环境的适应能力, 对环境变化进行检测是求解动态优化问题要完成的首要工作. 环境检测主要包括两个部分: 勘探粒子的选取以及衡量环境变化的标准设置. 现有勘探粒子的选取是随机的, 衡量环境变化的标准则依据勘探粒子适应值变化情况而定. 常用衡量变化的标准有如下三个

(1) 满足$ \left | f( x_{i}, t+1)- f( x_{i}, t) \right | >\beta $的勘探粒子个数超过勘探粒子总数的一定比例$ \eta $$ (\beta \ge 0, $$ 0<\eta <1) $;

(2) $ \frac{\sum\limits_{i=1}^{Km_{1} +m_{2}} \left | f( x_{i}, t+1)- f( x_{i}, t) \right | }{Km_{1} +m_{2}} >\beta (\beta \ge 0) $;

(3) $ \min\limits_{1\le i\le Km_{1} +m_{2}} \left | f( x_{i}, t+1)- f( x_{i}, t) \right | >\beta (\beta \ge 0) $.

在实际问题中, 依据最优解的精度要求选择合适的标准并设置阈值$ \beta $和$ \eta $.

环境检测的实质是检测最优解是否变化及如何变化. 本节将改进现有勘探粒子的选取方式, 在一定程度上避免现有方法的随机性和盲目性, 提高环境检测的准确性. 首先利用函数最优解区间收缩方法, 找出所有适应值备选函数的最优解所在区域(称这些区域为重点区域), 在重点区域设置一定数目的勘探粒子, 而在其他非重点区域设置少量勘探粒子, 然后根据勘探粒子的适应值变化情况来判断环境是否变化. 具体的, 由函数最优解区间收缩方法得出适应值备选函数$ g_{1} (x), g_{2} (x), \cdots , g_{K} (x) $的最优解所在区域分别记为$ A_{1} , A_{2} , \cdots , A_{K} $, 在这些区域分别随机选取$ m_{1} $个点, 在其他区域随机选取$ m_{2}(m_{2}<m_{1}) $个点, 共同组成规模为$ Km_{1} +m_{2} $的勘探种群$ x_{1}, x_{2}, \cdots , x_{Km_{1} +m_{2}} $, 计算前后两代勘探种群的适应值. 衡量环境变化的标准依旧选择第一种标准, 当前后两代勘探种群的适应值满足所选标准时, 认为环境发生变化; 否则, 认为环境未发生变化.

就粒子群算法和差分进化算法等群智能算法而言, 对环境变化的响应实质上是选取新的种群来适应变化后的环境, 使其能更好追踪变化后的最优解. 因此环境变化后, 新的种群要在最优解周围选取个体. 此外, 动态优化问题的环境虽然是变化的, 但某些环境也会在不同时刻再次出现. 因此, 算法设计中要充分利用前面的环境信息, 在响应环境变化的策略中引入记忆机制. 每次检测到环境发生变化后, 将上一环境信息以$ \left \langle W_{1}, W_{2} \right \rangle $的形式存储到记忆中, 其中$ W_{1} $中包含上一环境的编号, $ W_{2} $中包含上一环境中最后一代个体, 直到将有限个环境信息全部存储到记忆中为止. 当环境重复出现时, 直接利用记忆中存储的环境信息, 使得过去搜索到的有用信息传递到当前的种群中, 加快算法收敛速度.

检测到环境发生变化后的响应步骤如下

(1) 检测上一环境信息是否已存储, 将未存储的上一环境信息存储到记忆中, $ W_{1} $中存储上一环境编号, $ W_{2} $中存储上一环境中最后一代个体.

(2) 判断变化后的环境是否为新环境, 若是新出现的环境, 则在环境变化后的最优解区域内随机选取$ lN(0<l<1) $个个体, 在定义域内随机选取$ (1-l) N $个个体, 组成规模为$ N $的种群; 若是之前已经出现过的环境, 则从$ W_{2} $中直接调取对应环境的最后一代个体组成当前种群.

(3) 计算种群中所有个体的适应值, 并重新计算个体最优解的适应值.

(4) 将当前个体适应值与个体最优解的适应值比较, 选择较优个体的位置作为新的个体最优解, 从新的个体最优解中选择最优个体的位置作为新的全局最优解.

由于动态抛物线函数的动态变化易于控制, 领域内许多文献[5, 9-11, 14]都选用下式作为测试函数

(4.1) $ \begin{equation} f(x)=\sum\limits_{i=1}^{D} (x_{i}-offset)^{2} , x_{i} \in \left [ -50, 50 \right ], \end{equation} $

其中, $ offset $用于控制抛物线的动态变化, 每经历$ Q $代, $ offset $取值变化一次, 函数也就发生一次变化. 实验参数设置: 维数$ D=10 $, $ m _{1}=10, m _{2}=5, l=\frac{5}{6} $. 粒子群算法相关参数及其它参数设置与文献$ [10, 14] $相同.

比较本文算法与IPSO算法^[10]、LDPSO算法^[14]的收敛性. 初始时$ offset=0 $, 算法搜索最优解, 记录下适应值精度达到0.0001时算法的迭代次数. 然后函数动态变化, 即改变步长$ offset $, 继续搜索最优解, 记录下适应值再次达到精度要求时算法的迭代次数. 每次测试, 算法重复运行100次, 取平均值作为测试结果. 实验结果如表 1所示, 表 1为三种算法在环境变化前后收敛到最优解的平均迭代次数. "before"表示环境变化前算法收敛到最优解的平均迭代次数, "after"表示环境变化后算法再次收敛到最优解的平均迭代次数. 标$ * $的数据来源于文献[10], 标$ + $的数据来源于文献[14].

由表 1可知, 环境变化前后本文算法收敛到最优解的平均迭代次数普遍较少. 比较表 1中环境变化后三种算法再次收敛到最优解的平均迭代次数容易发现, 环境变化越剧烈, IPSO算法和LDPSO算法在环境变化后跟踪到最优解所需的迭代次数越多, 而本文算法的迭代次数与环境变化程度无关, 即使环境变化剧烈也能在较短时间内跟踪到新环境下的最优解. 因此本文算法在动态环境中的收敛性能优于其它两种算法. 而三种算法只有环境检测与响应方法不同, 其他设置完全相同, 由此说明, 本文的环境检测与响应方法提高了粒子群算法的收敛速度.

IPSO算法和LDPSO算法在环境检测时考虑了种群多样性, 多样性较差时才响应环境变化. 而本文算法在所有适应值备选函数的最优解区域设置了勘探粒子, 环境变化当代就准确检测到环境变化, 立即采取响应策略. 响应环境变化时, IPSO算法和LDPSO算法都是对种群进行某种扩散运动, 并没有准确找出变化后的最优解所在的区域, 因此粒子有“迷茫期”, 环境变化后不能及时跟踪新的最优解. 而本文算法在变化后的最优解周围选取个体以适应新的环境, 更具有针对性, 因此本文算法收敛速度更快, 对动态环境适应性更强.

本节针对受随机扰动的动态优化问题, 给出扰动前后环境检测无变化对应的随机扰动标准差的上限, 并检验本文环境检测与响应方法处理此类问题的有效性.

为了使随机扰动不影响动态优化问题的环境检测与响应, 需在一定的概率下保证扰动后的最优解仍位于原最优解区域. 本节将给出扰动前后环境检测无变化对应的随机扰动标准差上限.

设定最优解位置度量$ \gamma(0\le \gamma\le \frac{1}{2} ) $, $ x_{k}^{d} $为第$ k $个适应值备选函数的最优解在$ d $维上的分量, $ x_{k}^{d} $的最优解区间长度为$ L_{k}^{d} $, 则$ x_{k}^{d} $到最优解区间两端点的最短距离为$ \gamma L_{k}^{d} $. 要使扰动后的最优解仍位于原最优解区域, 只需各个方向上的扰动大小$ \left |\varepsilon _{k}^{d} \right | $满足

(5.1) $ \begin{equation} \left | \varepsilon _{k}^{d} \right | \le \gamma L_{k}^{d}, \forall 1\le k\le K, 1\le d\le D. \end{equation} $

给定置信水平$ \alpha (0.95\le \alpha <1) $, 使得$ P\left \{ \left | \varepsilon _{k}^{d} \right |\le \gamma L_{k}^{d} \right \} \ge \alpha $. 由于$ \varepsilon _{k}^{d} $服从正态分布$ N(0, \sigma ^{2} ) $, 要使$ P\left \{ \left | \varepsilon _{k}^{d} \right |\le \gamma L_{k}^{d} \right \} \ge \alpha $, 只需$ \frac{\gamma L_{k}^{d} }{\sigma } \ge u_{\frac{1+\alpha }{2} } $, 其中$ u_{\frac{1+\alpha }{2} } $为标准正态分布的$ \frac{1+\alpha }{2} $分位点. 从而标准差$ \sigma $与最优解位置度量$ \gamma $需满足如下关系

(5.2) $ \begin{equation} \sigma \le \frac{L_{k}^{d} }{u_{\frac{1+\alpha}{2} } } \gamma. \end{equation} $

由(5.2) 式可知, 要以$ \alpha (0.95\le \alpha <1) $的概率保证扰动后的最优解仍位于原最优解区域, $ \gamma $相同时, $ \sigma $的取值范围也相同, 即最优解位于关于区间中点对称的两个位置时, $ \sigma $的取值范围是相同的; $ \gamma $不同时, 对应扰动的标准差$ \sigma $的取值范围也不同. $ \gamma $越大, 即最优解距离最优解区域边界越远, 则对应的$ \sigma $能取到的值越大.

由于动态优化问题最优解的具体位置未知, 则$ \gamma $未知. 下面进一步分析扰动的标准差$ \sigma $与最优解位置度量$ \gamma $之间的关系, 进而确定扰动前后环境检测无变化对应$ \sigma $的上限.

命题5.1  设原目标函数的最优解$ x^{*} \in [0, 1] $, 加上服从$ N(0, \sigma ^{2} ) $分布的扰动$ \varepsilon $后, 最优解$ x^{*}+\varepsilon \in [0, 1] $的概率为$ P(\gamma , \sigma ) $, 则$ P(\gamma , \sigma ) $具有如下性质

(i) 固定$ \sigma $, 关于$ \gamma(0\le \gamma\le \frac{1}{2} ) $单调递增, 当$ \gamma=\frac{1}{2} $时取得最大值$ P(\frac{1}{2} , \sigma ) $;

(ii) 固定$ \gamma $, 关于$ \sigma $单调递减.

证

(5.3) $ \begin{eqnarray} P(\gamma , \sigma ) &\triangleq& P\left \{ x^{*}+\varepsilon \in [0, 1] \right \} =P\left \{0\le \gamma +\varepsilon \le 1 \right \}\\ &=&P\left \{- \frac{\gamma}{\sigma } \le \frac{\varepsilon }{\sigma } \le \frac{1 -\gamma}{\sigma } \right \} =\Phi (\frac{1 -\gamma}{\sigma })-\Phi (-\frac{ \gamma}{\sigma })\\ &=&2\Phi (\frac{ \gamma}{\sigma })+\left [ \Phi (\frac{1 -\gamma}{\sigma })-\Phi (\frac{ \gamma}{\sigma })\right ] -1. \end{eqnarray} $

(i) 对任意的$ 0<\gamma _{1} <\gamma _{2}\le \frac{1}{2} $, 由(5.3) 式可得

$ \begin{eqnarray*} P(\gamma _{2} , \sigma )-P(\gamma _{1} , \sigma ) & =&2\left [ \Phi (\frac{\gamma _{2} }{\sigma })- \Phi(\frac{\gamma _{1} }{\sigma }) \right ] + \left [ \Phi (\frac{1-\gamma _{2} }{\sigma })- \Phi(\frac{\gamma _{2} }{\sigma }) \right ]- \left [ \Phi (\frac{1-\gamma _{1} }{\sigma })- \Phi(\frac{\gamma _{1} }{\sigma }) \right ] \\ &=&2\left [ \Phi (\frac{\gamma _{2} }{\sigma })- \Phi(\frac{\gamma _{1} }{\sigma }) \right ] - \left [ \Phi (\frac{\gamma _{2} }{\sigma })- \Phi(\frac{\gamma _{1} }{\sigma }) \right ]- \left [ \Phi (\frac{1-\gamma _{1} }{\sigma })- \Phi(\frac{1-\gamma _{2} }{\sigma }) \right ]\\ & =&\left [ \Phi (\frac{\gamma _{2} }{\sigma })- \Phi(\frac{\gamma _{1} }{\sigma }) \right ] - \left [ \Phi (\frac{1-\gamma _{1} }{\sigma })- \Phi(\frac{1-\gamma _{2} }{\sigma }) \right ]. \end{eqnarray*} $

因为

$ \frac{1-\gamma _{1} }{\sigma } -\frac{1-\gamma _{2} }{\sigma }=\frac{\gamma _{2}-\gamma _{1} }{\sigma }, \frac{1-\gamma _{2} }{\sigma }\ge \frac{\gamma _{2} }{\sigma }(\gamma _{2} \le \frac{1}{2} ), $

由图 3容易得

$ \Phi (\frac{\gamma _{2} }{\sigma } )-\Phi (\frac{\gamma _{1} }{\sigma } )> \Phi (\frac{1-\gamma _{1} }{\sigma } )-\Phi (\frac{1-\gamma _{2} }{\sigma } ), $

所以$ P(\gamma _{2} , \sigma )-P(\gamma _{1} , \sigma )>0. $

(ii) 对任意的$ \sigma _{1} <\sigma _{2} $, 则$ -\frac{\gamma }{\sigma _{1} } <-\frac{\gamma }{\sigma _{2} }, \frac{1-\gamma }{\sigma _{2} } <\frac{1-\gamma }{\sigma _{1} } $, 从而

$ \left [ -\frac{\gamma }{\sigma _{2} } , \frac{1-\gamma }{\sigma _{2} } \right ] \subset \left [ -\frac{\gamma }{\sigma _{1} } , \frac{1-\gamma }{\sigma _{1} } \right ], $

由(5.3) 式可得$ P(\gamma , \sigma _{2} )-P(\gamma , \sigma _{1} )<0. $所以, 固定$ \sigma $, $ P(\gamma , \sigma ) $关于$ \gamma(0\le \gamma\le \frac{1}{2} ) $单调递增, 当$ \gamma=\frac{1}{2} $时取得最大值$ P(\frac{1}{2} , \sigma ) $; 固定$ \gamma $, $ P(\gamma , \sigma ) $关于$ \sigma $单调递减. 命题5.1证毕.

为验证上述结论, 设原最优解$ x^{*} \in [0, 1] $, 置信水平$ \alpha =0.95 $, $ \gamma $分别取$ 0.1, 0.2, \cdots , 0.5 $, 根据(5.2) 式, 计算$ \sigma $的取值范围并取最大值, 则$ \sigma $分别取$ 0.0510 $, $ 0.1020 $, $ 0.1531 $, $ 0.2041 $, $ 0.2551 $. 每组实验随机产生10000个扰动, 计算不同位置的最优解加上扰动后仍位于原最优解区域的概率, 称之为成功率$ P^{*} $, 重复10次实验, 取平均值作为实验结果. 实验所得结果如图 4所示.

从实验结果可以看出

(1) 关于区间中点对称的两个最优解的成功率曲线几乎完全重合, 因此扰动后的最优解仍位于原最优解区域所对应随机扰动标准差$ \sigma $的取值与最优解位置度量$ \gamma $有关, $ \gamma $相同时, $ \sigma $的取值范围也相同.

(2) 对于同一$ \sigma $值, 最优解距离最优解区域边界越远, 成功率越高; 对于同一位置的最优解, $ \sigma $越小, 成功率越高.

(3) 观察图 4容易发现

当$ \sigma=0.0510 $时, 若最优解$ x^{*} $满足$ 0.1\le x^{*}\le0.9 $, 则成功率$ P^{*} \ge 0.95 $;

当$ \sigma=0.1020 $时, 若最优解$ x^{*} $满足$ 0.2\le x^{*}\le0.8 $, 则成功率$ P^{*} \ge 0.95 $;

当$ \sigma=0.1531 $时, 若最优解$ x^{*} $满足$ 0.3\le x^{*}\le0.7 $, 则成功率$ P^{*} \ge 0.95 $;

当$ \sigma=0.2041 $时, 若最优解$ x^{*} $满足$ 0.4\le x^{*}\le0.6 $, 则成功率$ P^{*} \ge 0.95 $;

当$ \sigma=0.2551 $时, 若最优解$ x^{*} $满足$ x^{*}=0.5 $, 则成功率$ P^{*} \ge 0.95 $.

图像上其它点的成功率虽然没有达到0.95, 但是都大于0.65, 说明由$ \gamma $取$ 0.1 $, $ 0.2 $, $ 0.3 $, $ 0.4 $, $ 0.5 $得到的五个$ \sigma $的取值范围都能以0.65以上的概率保证扰动后的最优解仍位于原最优解区域. $ \sigma $越小, 成功率越高, 但是对于扰动的限制也越高, 因此要找到$ \sigma $能取到的上限.

下面从算法性能角度确定扰动前后环境检测无变化所对应的随机扰动标准差$ \sigma $的上限, 实验步骤如下

(1) 利用函数最优解区间收缩方法找出$ f(x)=\sum\limits_{i=1}^{10}x_{i}^{2} $, $ x_{i}\in [-50, 50] $的最优解所在区域为$ x_{i} \in [-50/81, 50/81], i=1, 2, 3, \cdots , 10 $;

(2) $ \gamma $分别取$ 0.1, 0.2, 0.3, 0.4, 0.5 $, 令$ \sigma = \frac{L_{k}^{d} }{u_{\frac{1+\alpha}{2} } } \gamma $, $ \alpha=0.95 $, 随机产生5组服从正态分布$ N(0, \sigma ^{2} ) $的扰动, 每组1000个数据, 即$ \varepsilon _{1} , \varepsilon _{2} , \cdots , \varepsilon _{1000} $;

(3) 设置算法运行60代, 前50代适应值函数为$ f_{1}(x)=\sum\limits_{i=1}^{10}x_{i}^{2} $, $ x_{i}\in [-50, 50] $, 第51代到60代适应值函数为$ f_{2}(x)=\sum\limits_{i=1}^{10}(x_{i}+\varepsilon _{j}) ^{2} $, $ j=1, 2, \cdots , 1000 $, $ x_{i}\in [-50, 50] $, 求适应值函数的最小值. 记录算法结束时每组最优解适应值的误差, 并计算每组平均误差.

第51代适应值函数发生变化后, 两种算法采取不同的策略. 在算法一中, 认为变化后的最优解位于原最优解附近, 因此认为此时环境不变, 利用第50代的粒子继续迭代. 算法二在原最优解区域重新随机选取粒子组成新的种群参与迭代. 除此之外, 两种算法完全相同. 实验结果见表 2, 表 2给出了$ \gamma $取不同的值时, 两种算法的平均误差.

由表 2可知, 当$ \gamma \le 0.3 $时, 算法一的误差小于算法二的误差, 算法一的性能更为优良. 说明此时扰动后的最优解位于原最优解附近, 可直接利用扰动前的粒子继续迭代, 寻找新的最优解, 这种做法充分利用了前面的信息, 在一定程度上提高了算法性能. 因此, $ \gamma $能取到的最大值为0.3, 由(5.2) 式可得扰动前后环境检测无变化所对应随机扰动标准差$ \sigma $的取值范围如下

(5.4) $ \begin{equation} \sigma \le \frac{L_{k}^{l} }{u_{\frac{1+\alpha}{2} } } \times 0.3. \end{equation} $

加上服从均值为0, 标准差满足(5.4) 式的正态分布的随机扰动之后, 最优解发生了一定的偏移, 但以一定的概率保证了扰动后的最优解仍位于原最优解附近, 因此认为环境无变化, 仍可利用扰动前的环境信息追踪新的最优解.

为验证本文改进的环境检测与响应方法处理受随机扰动的动态优化问题的有效性, 本节将通过两个常用测试函数^{[10, 15]}进行测试. 设置环境每隔一定代数连续发生变化, 最优解受服从均值为0, 标准差满足(5.4) 式正态分布的随机扰动而发生随机偏移.

测试函数1 (动态抛物线): 适应值函数设置如表 3所示, 其中$ f_{1} (x)=\sum\limits_{i=1}^{10} x_{i}^{2} $, $ f_{2}(x)=\sum\limits_{i=1}^{10} (x_{i}-10)^{2} , x_{i} \in \left [ -50, 50 \right ] $, $ \varepsilon _{1} $和$ \varepsilon _{2} $为随机扰动. 比较本文原算法与删去响应策略中记忆机制后的算法在相同环境中跟踪最优解的性能. 删去记忆机制后的算法不存储环境信息, 除此之外与原算法完全相同. 实验结果如图 5所示.

环境在第31、81、111代发生变化, 两种算法中环境检测方法相同, 由图 5可知, 每次环境发生变化, 两种算法总能在环境发生变化当代检测到环境变化并迅速响应, 但响应方法不同. 第81代和第111代时, 之前环境重复出现且最优解是由之前相同环境中的最优解发生随机偏移得到的, 由图 5(a)可知, 原算法直接利用记忆中存储的环境信息进行迭代, 加快了收敛速度, 在第120代左右算法就收敛到了最优解, 同时精度达到了0.00001. 删去记忆机制后的算法没有利用之前的环境信息, 由图 5(b)可以看出, 收敛速度较慢, 在算法结束时仍未达到0.0001的精度要求.

分别计算两种算法每代全局最优解适应值的误差并求平均, 可得原算法的平均误差为0.0313, 不含记忆机制的算法的平均误差为0.0696, 原算法的平均误差更小. 由此可知, 本文算法的响应策略中含有记忆机制, 在重复出现的环境中跟踪最优解的性能更为优良.

测试函数2 (移动峰函数): 动态优化问题目标函数如下

(5.5) $ \begin{equation} f(x, t)=\max\limits_{1\le i\le 2} \frac{h_{i}(t) }{1+\omega _{i}(t)\sum\limits_{j=1}^{5}[x^{j}-p_{i}^{j}(t)]^{2}}, \end{equation} $

其中, $ h_{i}(t) $、$ \omega _{i}(t) $分别为第$ i $个峰在$ t $时刻的高度和宽度, $ p_{i}^{j}(t) $为第$ i $个峰在$ t $时刻第$ j $维上的位置. 设置$ \omega _{i}(t) $为1, $ h_{i}(t) $在$ [1, 50] $上变化, $ p_{i}^{j}(t) $在$ [0, 12] $上变化. 每隔100代环境发生一次变化, 环境共变化5次, 最优解受服从均值为0, 标准差满足(5.4) 式正态分布的随机扰动, 发生随机偏移. 每当环境改变时, 记录当前环境编号, 得到当前环境下算法搜索到的最优解及对应的适应值, 算法重复运行100次, 计算环境检测的准确率, 并计算平均绝对误差和平均相对误差, 结果如表 4所示.

分析表 4可知, 环境变化时算法的绝对误差和相对误差均在合理的范围内, 算法在每次环境变化后都可以逼近真实最优解, 环境检测的准确率高达$ 96\% $, 第3次环境变化对应的误差稍大是仅有的几次环境检测失误导致的.

针对受随机扰动的动态优化问题, 本文提出的环境检测与响应方法不仅能够准确检测到环境变化, 对于最优解发生随机偏移的情况, 响应策略可直接调取记忆机制中存储的过去相同环境中的最后一代个体参与迭代, 充分利用了过去有用的环境信息, 同时也加快了算法收敛速度.

动态优化问题广泛存在于生产调度、人工智能、工程设计等诸多领域, 环境检测与响应是解决这类问题的前提. 然而在已有的环境检测与响应方法中, 还存在着环境检测不准确和过去的信息利用不充分等不足. 本文将正交试验设计思想用于环境检测, 给出了函数最优解区间收缩方法, 进而改进了动态优化问题的环境检测与响应方法. 之后, 确定了扰动前后环境检测无变化所对应的随机扰动的标准差上限, 并验证了本文环境检测与响应方法对于受随机扰动的动态优化问题的有效性.

本文提出的环境检测与响应方法通过正交试验设计对适应值备选函数预集成处理, 即利用函数最优解区间收缩方法找出所有适应值备选函数的最优解所在区域, 在此基础上设置勘探粒子检测环境变化, 提高了环境检测的准确性. 环境变化后, 根据最优解变化情况采取有针对性的响应策略, 并引入记忆机制存储环境信息, 当环境重复出现时, 直接利用之前环境中的最后一代个体进行迭代, 充分利用了过去有用的环境信息. 实验结果表明该方法能够准确检测并跟踪最优解的变化, 提高了算法对动态环境的适应能力. 因正交试验设计受因素个数的限制, 故本文提出的环境检测与响应方法难以处理高维自变量动态优化问题, 将在后续的工作中尝试利用灵敏度分析方法进行降维处理.