种群具有离散Leslie年龄结构的动力学优化算法

黄光球，陆秋琴

西安建筑科技大学管理学院，西安 710055

群智能优化算法[1]是启发式算法的一种，其基本求解策略不依赖优化问题的特征（如梯度、多次可微性等），而是依赖于试探解之间的信息交换能力。群智能优化算法之所以采取此策略，是因为优化问题在很多情况下无法求得其梯度或者不具备任何可微性。为了使得试探解之间的信息交换具有合理性，群智能优化算法经常依据一些特殊的自然现象构造而得。因此，一种自然现象的质量优劣对群智能优化算法的性能影响重大。

目前，具有代表性的新一代群智能算法有乌鸦算法[2]、海洋掠食者算法[3]、完全迁徙型生物地理学算法[4]、人工藻类算法[5]、磷虾群算法[6]、鲸鱼算法[7]、杂交蜂群算法[8]、反向学习粒子群算法[9]、花授粉算法[10]等。这些算法均来自一些特殊的自然现象，但这些自然现象有一个共同的缺陷就是无法用数学模型进行描述。因此，基于这些自然现象所构造出来的算法的性能设计无法奠定在坚实的数学分析基础之上，该缺陷也是传统群智能算法的缺陷。此外，因这些算法没有数学模型基础，故相关参数选择也相当困难。于是，若一个自然现象能够被数学模型恰当描述，基于该自然现象的群智能算法就可以克服上述缺陷。

一个自然现象能否成为构造一个群智能优化算法的生物活动场景，必须包含有两个条件：一是生物个体不断进化，符合“适者生存，优胜劣汰”的自然规律；二是生物个体之间存在相互作用关系。第一个条件用来构造算法的逻辑结构并确保该算法具有全局收敛特征，第二个条件用来构造算子，个体类型越多，个体之间相互作用关系越丰富，所能构造出来的算子也就越多，算子的性能也越强。此外，若该自然现象还能被数学模型恰当描述，则据此所构造出来的群智能优化算法的性能分析和参数选择将具有很好的科学依据[11]。

由具有不同性别、不同年龄的生物个体所组成的生物种群在一个生态环境中生存进化，是一种很常见的自然现象[3]，该现象就满足构造一个群智能优化算法的两个条件。不仅如此，该自然现象既能被数学模型很好描述，又能提供高度多样化的个体间相互作用关系。本文正是依据该自然现象提出了一种新的群智能算法，即种群具有离散Leslie 年龄结构的动力学优化算法，简称PDO-DLAS（population dynamic optimization algorithm with discrete Leslie age structure）算法，该种群动力学模型简称为Leslie模型[12]。Leslie 模型具有较好的应用价值，已得到较好的研究，目前研究进展如下：

（1）种群特征分析[13-18]。在该方面，通过构建Leslie模型，研究种群的特征，如扩散特征、稳定性特征、持久和湮灭特征等。

（2）种群控制研究[19-22]。在该方面，通过构建Leslie控制模型，研究种群规模的最优控制或收获问题。

PDO-DLAS 算法依赖于Leslie 模型，因此上述模型的研究可为PDO-DLAS 算法的构造奠定良好的理论基础。与文献[2-10]描述的算法相比，该算法具有如下特点：

（1）Leslie 年龄结构使得个体自动依据年龄分成多个类别，可大幅提升个体的多样性，从而降低个体陷入局部陷阱的概率。

（2）算法参数的选择依据Leslie 数学模型确定，因而具有坚实的理论依据。

（3）依据Leslie 模型描述的生态场景设计出来的算子具有充分的理论依据。

（4）个体进化满足达尔文进化论规律，既可确保算法具有全局收敛性，又可使得个体进化每次只涉及极少特征的改变，可大幅提升计算速度。

由于PDO-DLAS 算法具有独特的优势，对其进行研究很有意义。

1 PDO-DLAS 算法构建原理

考虑一类优化问题：

式中，Rn是n维欧氏空间，n为优化问题的维数；X=(x1,x2,…,xn)是一个n维解向量；f(X)为目标函数；H为搜索空间。

1.1 决策变量与生物个体的对应关系

假设某种群在一个生态系统E中生长，该种群由具有不同年龄的生物个体组成，为简单起见，下面将生物个体简称为个体。个体分为雌性（用性别变量u=1 表示）和雄性（用性别变量u=2 表示）两种类型；种群自然地分成A个年龄组，A为该种群的最高年龄，0

在时期t，该种群中的性别为u的个体总数为Pu(t),u∈{1,2}，即：

每个个体都用唯一编号表示，即性别为u、年龄组为a的个体编号集合为{1,2},a=1～A。种群中的每个个体都具有n个特征，且特征数n不随时间变化而改变。时期t，对于性别为u、年龄组为a的个体i来说，用其特征表示就是其中就是性别为u、年龄组为a的个体i的第j个特征，

优化问题式（1）中的决策变量与个体的对应关系如下所述。

时期t，在式（1）的搜索空间H中随机生成的Pu(t)个试探解是令搜索空间H与生态系统E相对应，则时期t该生态系统中Pu(t)个个体就与搜索空间H中式（1）的Pu(t)个试探解一一对应，即时期t性别为u、年龄组为a的个体i与试探解一一对应；更进一步，即时期t性别为u、年龄组为a的个体i的特征与试探解的分量相对应，u∈{1,2}。

在该生态系统中，个体的演化规律如下：

（1）个体的新生是由雌性个体和雄性个体交配产生的。

（2）个体的死亡是因为个体太虚弱所导致的。

（3）在年龄相近的个体中，普通个体会不分性别地向比它优秀的个体学习。

（4）年龄较大的个体能够不分性别地对年龄较小的施加影响。

（5）个体之间的相互作用表现在对其特征的影响上，且此相互作用及其所导致的影响是随时间变化的。

对任何一个个体来说，按优化问题式（1）计算，其所对应的目标函数值越小，则其适应度越高，反之亦然。个体的适应度越高，其生存的概率也越高；而适应度越低的个体，其被淘汰的概率会越高。也就是说，生态系统中的个体的演化规律符合达尔文进化论规律。时期t性别为u、年龄组为a的个体i的适应度用FEI（fitness evaluation index）指数来表示，其计算方法为：

1.2 Leslie种群离散模型

假设表明第a年龄组每一个个体在一个时段内出生的后代（或新生的雌性或雄性）的平均数为出生率函数，特别地设它是一个常数，即：

由所有年龄组的后代相加而得到初始年龄组为：

假设描述由第a年龄组存活到第a+1 年龄组的比率为存活率函数它只是第a年龄组个体数的线性函数。

假设式（4）和式（6）表明不考虑参数依赖于环境条件的变化，并且忽略种群大小对出生率和死亡率的影响。

令Yu(t)=，则由式（5）、式（7）可得到方程：

这里，A×A方阵Lu形如：

Lu称为Leslie 矩阵，其元素称为活率或种群统计参数。称式（8）、式（9）为Leslie模型。

向量Yu(t)以绝对值而非相对值显示了种群大小在各年龄组中的一个分布，相对年龄分布称为年龄普。矩阵Lu确定了A维欧氏空间的一个线性算子，把这种算子称为Leslie算子。因为变量代表第a年龄组的个体数，它们是非负的。矩阵Lu的特征方程为：

如果矩阵Lu有模等于r的h个特征值，即λ1=r,λ2,…,λh，则所有这些特征值彼此不同，但都是λhrh=0 的根，h称为矩阵Lu的非本原性指数，它等于那些出生率不为0 的年龄组的编号的最大公约数。若h=1，则矩阵Lu称为素矩阵。

矩阵Lu可表示为：

式中，m为矩阵Lu的不同特征值的个数；Jj为kj×kj（kj为特征值λj的重数，k1+k2+…+km=A）Jordan块，其形状为：

Cu(λ)和Ru(λ)的内积满足：

(Cu(λa),Ru(λa))=1,(Cu(λa),Ru(λb))=0,a≠b,a,b=1～A

引入极限向量函数：

式（10）、式（11）可以对Leslie 模型的轨线的渐近性质建立一些推断。

若Lu是素矩阵，则极限向量函数与t无关，即：

定理1（关于循环的平均）[12]设h为矩阵Lu的非本原性指数，则：

定理2（关于周期性）[12]设h为矩阵L的非本原性指数，则极限向量函数L(Yu,t)关于t是周期的，周期为：

式中，lcm 表示最小公倍数；l取遍值1,2,…,h;ul/vl为不可约分数（当l=1 时，设ul=0,vl=1）。

这样，T是非本原性指数h的因子，并依赖于向量Yu的初始分布。任何情况下，极限分布的总体大小W(Yu,t)的真正周期必须是T的因子，其中W(Yu,t)=真正周期等于T的条件，由下面的定理给出。

定理3（关于种群总大小的周期）[12]W(Yu,t)的周期一致的充要条件是：

当R的值不作为稳定性的某种度量时，是所有生物种群再生速率的广义参数；而作为一种稳定性度量时，如果R=1，那么r=1，种群大小不会作指数性的减少和增长，而趋于极限分布状态。R<1 和R>1 分别表明r<1 和r>1，根据式（10）和式（12），种群或灭亡或者无限制地增长。

由式（14）知，成长系数R由存活率、出生率函数和年龄组数A决定。当确定后，可确定能够确保R=1 的A。

定理4（循环周期）[12]设Lu是一具有单位主特征值和非本原性指数h的Leslie 矩阵，则在模型式（8）的轨线之间存在循环，如果循环有周期k，则k是h的因子。

1.3 演化算子设计

时期t，对当前性别为u、年龄组为a的个体i来说，其特征个体集合的生成策略如下：

（1）优势个体集合SIu,a(t)。从个个体中随机选出M个个体，这些个体的FEI 指数要比当前个体要高，形成优势个体集合SIu,a(t)，M称为特征个体数。

（2）普通个体集合GIu,a(t)。从个个体中随机选出M个个体，形成普通个体集合GIu,a(t)。

演化算子设计方法如下：

（1）新生算子。该算子描述种群内个体的新生。从雌性和雄性优势个体集合SIu,a(t)中各随机选择一个个体进行交配，产生新个体i，其性别随机确定，所处年龄组为第1 组。

式中，w在{1,2}随机选择；α=Rnd(0,1),Rnd(a,b)为在区间[a,b]内产生满足均匀分布的随机数；a=A0～A,A0为雌性个体生育的起始年龄；iu在中随机选择，u∈{1,2};J0为在1～n中随机选择M个特征所形成的集合。由式（15）形成的新个体为

（2）死亡算子。该算子描述种群内个体的死亡。时期t，从性别为u、年龄组为a的个体集合中，选出FEI指数最低的个体，将其删除。

（3）学习算子。该算子描述的是处于相近年龄的普通个体会向其优势个体学习，从而使得优势个体的一些特征传递给普通个体。与年龄组a相近的年龄组取为U={a-M,a-M+1,…,a,a+1,…,a+M}，对于每个年龄组b∈U，任意选择被学习的个体性别w，令Q(b)=w,w∈{1,2}；进一步令则对于性别为u、年龄组为a的当前个体i来说，有：

（4）影响算子。该算子描述的是年龄较大的个体的一些行为会对年龄小的个体造成影响。比年龄组a大的年龄组取为V={a+1,a+2,…,a+M}，对于每个年龄组b∈V，任意选择施加影响的个体性别w，令Z(b)=w,w∈{1,2}；进一步令则对于性别为u、年龄组为a的当前个体i来说，有：

式中，YP=GP1⋃GP2。

（5）进化算子。个体的进化须满足达尔文进化论规律。对于性别为u、年龄组为a的当前个体i，其进化算子的定义如下：

1.4 个体初始化

假设优化问题搜索空间的维数为n，每个变量的搜索区间为[li,ui],i=1～n，利用正交拉丁方生成方法产生W个初始解的正交表LM(Wn)的构造算法如算法INIT 所述。

算法1INIT(W)//产生W个初始解的正交表LM(Wn)的构造算法

步骤1计算每个变量的离散点yij：

yij=li+(j-1)(ui-li)/(W-1),i=1～n;j=1～W

步骤2根据正交拉丁方的生成方法计算初始解xij：

xij=yjk,i=1～W,j=1～n

式中，k=(i+j-1)modW；若k=0，则k=W。

算法INIT(W)所确定的W个初始解Xi=(xi1,xi2,…,xin),i=1～W，具有很好的均衡分散性和整齐可比性。

1.5 PDO-DLAS 算法构造

算法2PDO-DLAS

步骤1初始化

步骤2执行下列操作

1.6 算法特点分析

1.6.1 时间复杂度

PDO-DLAS 算法的时间复杂度计算如表1 所示，其中

1.6.2 PDO-DLAS 算法的全局收敛性证明

定理5PDO-DLAS 算法具有全局收敛性。

证明（1）从新生算子、学习算子、影响算子、生长算子的定义式（15）～式（18）可知，满足关系表明PDO-DLAS 算法的演变过程具有Markov 特性；（2）从生长算子的定义式（18）知，任一个体的FEI 指数总是保持递增态势。依据文献[23]可知，满足上述两个特性的PDO-DLAS算法具有全局收敛性，其相关证明可参见文献[23]，本文不再赘述。

Table 1 Time complexity表1 时间复杂度

2 参数确定

PDO-DLAS 算法的参数由Leslie 模型参数和运行控制参数两部分组成，前者是算法的内置参数，一旦设定，无需修改；后者需要根据所求优化问题的实际特征进行人工设置。

2.1 Leslie模型参数设置方法

依据1.2 节的结论，Leslie 模型参数设置应满足如下要求：

（2）从定理2 和定理3 知，轨线的性质依赖于时间尺度的选择，故应避免模型轨线的周期性。例如，可合并时间间隔或延长时段，使得有非零出生率的年龄组成为相邻的。

（3）从定理4 知，应尽量使得Lu的主特征值不等于1，从而确保模型式（8）的轨线之间存在循环。

依据上述要求，Leslie 模型式（4）、式（6）和式（9）中的参数可按下列方式设置：A=20,A0=A/2,h=A0+1=11,n0=30,n1=50;b0=3,b1=5;s0=0.65,s1=0.846。此时，个体总数随演化周期的变化如图1所示。由式（14）计算得R=1.000 378 ≈1;Lu的主特征值λ=1.000 3 ≈1，式（8）的轨线之间的循环周期为11 个演化周期。从图1 可知，各年龄组个体数趋向极限分布，且具有很好的变化特征。

Fig.1 Changing rule of total number of individuals with time图1 个体总数随时间的变化规律

2.2 运行控制参数设置方法

含有n个变量的Michalewicz 函数优化问题具有n!个局部极小点。例如，对于50 维优化问题，Michalewicz 函数有高达50!=3.041 41×1064个局部极小点。因此，搜索过程极易陷入局部最优解陷阱。目前，该函数的理论全局最优解还未可知，寻找该函数的全局最优解具有很大的挑战性。Michalewicz 函数优化问题如下：

本文以求解Michalewicz 函数优化问题为例来确定PDO-DLAS 算法中的运行控制参数的适当设置，这些参数包括如下几类：

（1）G和ε。这两个参数是互补参数，只要满足其中一个即可。通常取G=106～1010,ε=10-10～10-5。

（2）E0和M。这两个参数对PDO-DLAS 算法的性能影响较大，下面重点讨论这两个关键参数的取值规律。

令n=50,E0=0.01,G=108,PDO-DLAS 算法运行100 次，表2 描述了M与最优目标函数值的平均值（AvgOFV）和计算时间的平均值（AvgCT）之间的关系。表2 表明，当M=3～6 时，AvgOFV的精度达到最佳，而AvgCT递增不大。因此，建议M=3～6。

Table 2 Relationship of M with AvgOFV and AvgCT表2 M 与AvgOFV 和AvgCT 之间的关系

令n=50,M=3,G=108,PDO-DLAS 算法运行100 次，表3 描述了E0、AvgOFV和AvgCT之间的关系。结果表明，当E0=0.004～0.100 时，AvgOFV精度较高，但AvgCT增加不大；当E0>0.100 时，AvgCT增加很大，但AvgOFV精度也大大降低；尤其当E0=1.000时，无法得到最优解。因此，当E0=0.004～0.100 时，PDO-DLAS 算法性能最好。

Table 3 Relationship of E0 with AvgOFV and AvgCT表3 E0 与AvgOFV 和AvgCT 之间的关系

3 PDO-DLAS 算法与其他算法的比较

CEC2013[24]是一个国际上专门用来测试智能优化算法的测试包，该测试包中包含有28 个优化问题，每个优化问题由一些著名的传统优化问题组合而成。这些优化问题不但克服了传统优化问题的一些缺陷（如全局最优解关于原点对称、条件数过低等），而且其求解难度均大幅提高。本文在该测试包中选择6 个难度很大的优化问题来对PDO-DLAS 算法与其他算法进行比较，如表4 所示。

Table 4 6 optimization problems表4 6 个优化问题

表4 中，O是随机产生的n维向量。求解这些优化问题时，PDO-DLAS 算法的参数设置是n=50,G=108，ε=10-8,E0=0.01,M=3。与PDO-DLAS 算法进行比较的7 种智能优化算法为RCGA（real-coded genetic algorithm）[25]、DASA（differential ant-stigmergy algorithm）[26]、NP-PSO（non-parametric particle swarm optimization）[27]、MpBBO（metropolis biogeographybased optimization）[28]、MDE-LiGO（modified differential evolution with locality induced genetic operators）[29]、SLADE（self-adaptive differential evolution）[30]、ABC（artificial bee colony）[31]，这些算法各参数设置可参见其对应文献。

求解各个优化问题时，每个算法均独立求解51次。表5 给出了各算法的求解结果。其中，表5 的列是所有参与比较的优化算法，表5 的行是求解每个优化问题时所获得的关键对比参数，即各算法所获得的最优解的均值、中值、标准误、最小值、最大值、适应度评价次数。

从表5 可以看出，这8 个算法按最终排名1 进行排序所得的结果均如下：

PDO-DLAS>MDE-LiGO>SLADE>NP-PSO>RCGA>DASA=MpBBO>ABC

按最终排名2 进行排序所得的结果均如下：

PDO-DLAS>MDE-LiGO>SLADE>NP-PSO>RCGA=DASA>MpBBO>ABC

图2（a）～（f）给出了各算法求解6 个优化问题时的样本收敛曲线。各图样本曲线之间的差异可以用PDO-DLAS 算法的探索和求解过程来说明。

Table 5 Solution results of each algorithm表5 各算法的求解结果

表5（续）

（1）对于图2（a），PDO-DLAS 算法在区间[0，102]、[102，103]、[103，107]、[107，108]搜索期间，搜索分别处于求精状态、探索状态、再求精状态、快速探索状态转求精状态，而进行“求精→探索→再求精→快速探索”的状态转换等价于不断从局部最优解陷阱跳出的行为；其他算法在区间[0，102]、[102，106]、[106，108]处于初始化状态、求精状态、探索状态转再求精状态。PDO-DLAS、MDE-LiGO、SLADE 具有最终的相同求解精度，其他算法最终的求解精度均很低。

Fig.2 Sample convergence curves图2 样本收敛曲线

（2）对于图2（b），PDO-DLAS算法在区间[0，102]、[102，5×105]、[5×105，9×105]、[9×105，108]搜索期间，搜索分别处于求精状态、缓慢探索状态、快速探索状态、再求精状态，即进行“求精→缓慢探索→快速探索→再求精”的状态转换，实施相应的局部最优解陷阱逃逸操作；DASA、MpBBO、MDE-LiGO、SLADE 算法在区间[0，102]、[102，5×105]、[5×105，9×105]、[9×105，108]搜索期间，处于初始化状态、缓慢探索状态、快速探索状态、求精状态；RCGA、NP-PSO、ABC 算法在区间[0，102]、[102，9×107]搜索期间，搜索分别处于初始化状态、缓慢求精或缓慢探索状态，RCGA 在[9×107，108]期间进入快速探索状态。PDO-DLAS、MDELiGO、SLADE 具有最终的相同求解精度，其他算法最终的求解精度均很低。

采用上述类似的分析方法，可以分析图2（c）～图2（f），本文不再重复。

4 结束语

本文基于种群具有离散Leslie 年龄结构的动力学模型提出了一种具有全局收敛性的新型优化算法，与其他典型群智能算法相比，PDO-DLAS 算法具有如下特点：

（1）个体依据其年龄组和性别被自动划分成2A类，每类个体数依据Leslie 模型自动进行动态计算，这样既大幅增加了个体的多样性，又避免了人工确定个体数的困难。

（2）所有算子是通过利用Leslie 模型以及同龄组和不同龄组个体间的相互作用关系进行构造的，PDO-DLAS 算法与所求解的实际优化问题无关，故具有很好的普适性。

（3）PDO-DLAS 算法中的每个算子具有明确功能，其中学习算子可实现年龄组相近个体之间的信息交换；影响算子可实现不同龄组个体之间的信息交换；新生算子可增加强壮个体数；死亡算子可以减少虚弱个体数；进化算子可确保算法具有全局收敛性。

（4）采用正交拉丁方方法生成初始试探解可确保其具有均衡分散性和整齐可比性。

（5）采用Leslie 模型的机理确定PDO-DLAS 算法中的相关参数，大幅减少了PDO-DLAS 算法参数的人工确定个数。实际上，PDO-DLAS 算法需要人工确定的参数只有2 个。

（6）在进行迭代计算时，PDO-DLAS 算法每次只处理个体特征数的1/250～1/10，从而使时间复杂度大幅降低。因此，PDO-DLAS算法适于求解高维优化问题。

PDO-DLAS 算法的下一步改进方向如下：

（1）深入研究学习算子、影响算子、新生算子和死亡算子的动态特征。

（2）深入研究个体的动态特征。