布谷鸟搜索算法

布谷鸟搜索(Cuckoo Search,CS)是由 Xin-She Yang 和 Suash Deb 于 2009 年开发的自然启发式算法。CS 基于布谷鸟的寄生性育雏(brood parasitism,又巢寄生)行为。该算法可以通过所谓的 Levy 飞行来增强,而不是简单的各向同性随机游走。研究表明,该算法可能比遗传算法、PSO 以及其他算法更有效。

布谷鸟育雏行为

布谷鸟(杜鹃)是一种神奇的鸟,不仅因为它们动听的啼鸣,还因它们的积极的繁殖策略。杜鹃科中的犀鹃(Ani Cuckoo)和圭拉鹃(Guira Cuckoo),将它们的蛋放在其他鸟的巢中,通过去除其他鸟(寄主)的蛋来增加自己蛋的孵化几率。有相当多种类的鸟都有将自己的蛋放在其他鸟的巢中这种寄生性育雏行为 [19]。

寄生性育雏分为三种:种内寄生(intraspecific brood parasitism)、合作养育(cooperative breeding)和巢占据(nest takeover)。一些寄主鸟会与入侵的布谷鸟发生直接冲突。如果一个寄主鸟发现这些蛋不是他们自己的,那么他们要么将这些外来蛋清除掉,要么就直接放弃这个巢,在别处建造一个新的巢。一些布谷鸟,例如 New World brood-parasitic Tapera,已经进化成这样一种方式,雌杜鹃通常非常善于模仿几种特定寄主的卵的颜色和纹理。这减少了它们蛋被遗弃的可能性,从而增加了它们的繁殖力。

此外,该物种对产蛋时机的把握也非常到位。布谷鸟通常会选择那些寄主刚刚产下自己蛋的巢。一般来说,布谷鸟蛋的孵化时间要比寄主蛋的孵化时间要早一些。一旦第一只布谷鸟雏鸟孵化出来,第一个本能的动作就是通过盲目地推动将其他蛋从巢中推出,从而增加寄主对布谷鸟雏鸟的食物供给。研究还表明,杜鹃雏鸟还可以模仿寄主雏鸟的叫声,以获得更多的被喂食机会。

Levy 飞行

另一方面,各种研究表明,许多动物和昆虫的飞行行为表现出了具有幂律规律的 Lévy 飞行的典型特征。Reynolds 和 Frye 最近的一项研究表明,果蝇(或 Drosophila melanogaster)利用一系列直线飞行路径和突然的 90° 转弯来探索景观,从而产生 Lévy飞行式的间歇无标度搜索模式 [21]。针对人类行为的研究也表明,如 Ju/’hoansi 狩猎采集觅食模式等也表现出了 Lévy 飞行的典型特征 [4]。即使是光线也与 Lévy 飞行有联系 [2]。另外,该行为已被应用于优化搜索,结果表明其具有潜力 [20]。

布谷鸟搜索

CS 是由 Xin-She Yang 和 Suash Deb 于 2009 年开发的自然启发式算法 [30,32,33]。CS 基于布谷鸟的寄生性育雏(brood parasitism,又巢寄生)行为。

此外, 该算法可以通过所谓的 Levy 飞行来增强,而不是简单的各向同性随机游走。研究表明,该算法可能比遗传算法、PSO 以及其他算法更有效 [30]。为了简化描述标准 CS,这里我们引入以下三条理想化的规则:

  • 每只布谷鸟每次下一个蛋,并将其放入随机选择的巢中。
  • 具有优质蛋的最佳巢会被带到下一代。
  • 可用的寄主巢数量是固定的,且寄主以概率 发现布谷鸟放的蛋。在这种情况下,寄主可以消灭该蛋或放弃旧巢另建新巢。

进一步地,对于最后一个假设,新巢可以通过替换 个宿主巢穴的 来近似。对于最大化问题,解质量或适应度可以简单地假设为与目标函数的值成比例。其他形式的适应度可以用与遗传算法中的适应度函数类似的方式来定义。

从实现的角度来看,我们可以用下面的简单规则:每个巢中的蛋代表一个解,每个布谷鸟只能下一个蛋。目的是使用新的和可能更好的解(布谷鸟)来取代巢中不太好的解。显然,这个算法可以扩展到更复杂的情况,也就是,每个巢有多个蛋来代表一组解。这里,我们只考虑最简单的情况,每个巢只有一个蛋。在这种情况下,蛋,巢或布谷鸟之间没有区别,因为每个巢对应一个鸡蛋,这也代表一只布谷鸟。

该算法使用由开关参数 控制的局部随机游走和的全局探索随机游走的平衡组合。局部随机游走可以写成

其中 是通过随机置换选择的两个不同的解, 是一个 Heaviside 函数(单位阶跃函数)。 是从均匀分布中抽取的随机数, 是步长。这里, 表示两个向量的 entry-wise 积(点乘)。

另一方面,全局随机行走使用 Levy 飞行

其中

这里 是步长缩放因子,与感兴趣问题的尺度相关。大多数情况,我们使用 ,而 是感兴趣问题的特征尺度,在某些情况下 可能会更有效并且能避免飞得太远。显然,这两个更新方程的 值可能不同,因此,导致产生两个不同的参数, 。为了简化,这里我们使用

基于上述三条规则,CS 的基本步骤可以总结为以下伪代码:

Figure 1 Pseudo code of the cuckoo search for a minimization problem.

式(2)本质上是一种随机行走的随机(stochastic)公式。事实上,随机行走是一个马尔科夫链,其下一个状态/位置仅取决于当前状态(上式的第一项)和转移概率(上式的第二项)。然而,新解的很大一部分应该由远场随机化产生,它们的位置应该离当前最佳解足够远,这将确保系统不会陷入局部最优 [30, 32]。

关于布谷鸟搜索的文献在快速增长。它得到了广泛地关注,最近不同领域的很多研究都用到了布谷鸟搜索 [6,7,9-11,13,36]. 例如,Walton 等改进了该算法提出了修改的 CS 算法 [26];Yang 和 Deb 扩展该算法到多目标优化 [33]。一个全面的综述可以从 Yang 编写的书上看到 [35]。

布谷鸟搜索的特殊情况

CS 作为一种元启发式算法有着令人惊讶的丰富特性。如果我们仔细观察更新公式(1)(2),我们就能发现丰富的细节。对于式(1)我们可以把因子放在一起,设 ,然后我们有 。于是,式(1)变成了差分进化的主更新公式。更进一步,我们将 替换为当前最佳解 并设 ,我们有

这本质上是没有个体历史最佳的 PSO 变种。在这种情况下,非常类似 Yang 等开发的 APSO [29]。

另一方面,对于式(2),这种随机行走是具有 Levy 飞行转移概率的模拟退火(simulated annealing,SA)。在这种情况下,SA 的随机冷却表由 控制。

因此,差分进化、PSO 以及 SA 都可以看作是 CS 的特殊情况。相反,我们也可说 CS 将 DE、PSO 和 SA 好的和有效的部分组合在一个算法中。因此,CS 非常有效。

如何执行 Levy 飞行

从实现的角度来看,用 Lévy 飞行生成随机数应包括两个步骤:随机方向的选择和服从 Lévy 分布的步长的生成。方向的生成应该服从均匀分布,而生成步长是相当棘手的。有几种方法可以实现,但是最有效且直接的方法就是使用所谓的 Mantegna 算法来实现对称的 Lévy 稳定分布 [15]。

然而,生成正确服从 Lévy 分布的伪随机步长并不简单。 在 Mantegna 算法中,步长 s 可以通过以下变换使用两个服从高斯分布的变量 U 和 V 来计算:

其中

这里 意味着样本服从均值为 0 方差为 的高斯正态分布。方差可以使用下式计算:

该分布,对 服从预期 Lévy 分布,其中 是最小步长。原则上,,但实际上 只要值合理就可取,如 到 1。

这些公式看起来很复杂,但 函数对于给定的 是常数。例如,当 时,我们有 以及

Mantegna 算法在数学上已被证明能够产生符合服从要求的分布的随机样本 [15]。

参数选择

CS 中有多个参数。除了种群规模 n 之外,还有切换概率 ,步长缩放因子 和 Lévy 指数 。然而,关键的参数是 和 n,因为我们可以把 作为常数。通过改变它们的值,我们发现对于大多数问题,我们可以设

对于关键参数,我们也尝试改变寄主巢的数量(或种群规模 n)和概率 。我们使用了 以及 。从我们的模拟中,我们发现 到 40 和 对于大多数优化问题是足够的。结果和分析也表明收敛速度在一定程度上对所使用的参数不敏感。这意味着对于任何问题都不需要进行仔细地调整。

让我们看一个简单的例子。我们使用的许多测试函数之一有双变量 Michalewicz 函数

其中 。该函数在 (2.20319, 1.57049)处有全局最小值 。使用 CS 可以很容易地找到这个全局最优值,结果如图 2 所示,其中巢的最终位置也在图中以 标出。这里我们使用了 n = 15 个巢,

Figure 2 Search paths of nests using CS. The final locations of the nests are marked with in the figure.

从图中我们可以看到,随着最优值的接近,大多数巢朝着全局最优的方向汇聚。 我们还注意到,在多模态函数的情况下,巢也分布在不同的(局部)最优点。 这意味着,如果巢的数量远远大于局部最优的数量,那么 CS 可以同时找到所有的最优解。当我们处理多模态和多目标优化问题时,这个优势可能变得更加重要。

CS 变种

在过去的几年里,CS 的许多变种已经被开发出来。 尽管该算法只有短暂的历史,但由于其简单、高效和灵活,CS 引起了大量的关注。 结果就是,相关文献显著地增长。标准 CS 非常强大和高效,但它是为连续优化而开发的。一个有用的扩展是开发离散 CS,以使它可以有效地解决调度问题和组合优化。已经有很多 CS 的变种, 综述请参考 [34,12,35]。这里我们只概述几个变种:

  • 修改的布谷鸟搜索(Modified cuckoo search,MCS)。Walton 等人开发了修改的的布谷鸟搜索 [26]。已被用来优化网格生成和其他应用。

  • 改进的布谷鸟搜索(Improved cuckoo search,ICS)。在人工智能应用上, Valian 等人 [23] 提出了用改进的 CS 训练前馈神经网络 [23,24]。与此同时,Vazquez [25]也在使用 CS 来训练尖峰神经模型 [25]。

  • 量子布谷鸟搜索(Quantum cuckoo search,QCS)。Layeb [14] 通过给算法增加量子行为提出了一个 CS 的变体,称之为量子 CS。 QCS 已被应用于解决背包问题 [14]。

  • 离散布谷鸟搜索(Discrete cuckoo search,DCS)。对于调度和组合问题等离散应用,也有几个变种。Ouaarab 等人 [18] 开发了一个解决旅行商问题(TSPs)的离散 CS [18]。Chandrasekaran 和S imon [5] 提出了多目标 CS 方法来解决多目标调度问题 [5]。

  • 多目标布谷鸟搜索(Multi-objective cuckoo search,MOCS)。Yang 和 Deb [33] 提出了一个求解多目标工程优化问题的多目标 CS 方法 [33]。

  • 离散多目标布谷鸟搜索(Discrete multi-objective cuckoo search,DMOCS)。在多目标和调度问题的背景下,Chandrasekaran 和 Simon [5] 开发了 CS 的变体来解决离散多目标调度问题 [5]。

  • 混合布谷鸟搜索(Hybrid cuckoo search,HCS)。有些变体试图将 CS 与其他算法结合起来。例如,Wang 等人 [28] 将 CS 与 PSO 相结合,取得了很好的改进 [28]。Salimi 等人 [22] 结合了修改的 CS 与共轭梯度法 [22]。

还有很多其他的变种,Fister 等人提供了详细的综述 [12]。Yang 和 Deb 供了一个概念性的综述 [34]。最近编写的书有更多关于杜鹃搜索和萤火虫算法的文献 [35]。

为什么 CS 这么有效

除了前面的分析表明,DE、PSO 和 SA 是 CS 的特例,最近的理论研究也表明 CS 具有全局收敛性 [27],如下一小节所述。

PSO 的理论研究表明,它可以快速收敛到当前最优解,但不一定是全局最优解。事实上,一些分析表明,PSO 更新方程不满足全局收敛条件,因此不能保证全局收敛。另一方面,已经证明布谷鸟搜索能够满足全局收敛的要求,从而保证了全局收敛性 [27]。这意味着对于多模态优化,PSO 可能过早地收敛到局部最优,而 CS 通常可以收敛到全局最优。

此外,CS 具有两种搜索能力:局部搜索和全局搜索,由切换/发现概率控制。正如前面提到的那样,局部搜索是非常密集的,搜索时间约为 1/4(),而全局搜索约占总搜索时间的 3/4。这使得可以在全局范围内更高效地探索搜索空间,从而可以以更高的概率发现全局最优。

CS 的另一个优势是它的全局搜索使用 Lévy 飞行,而不是标准的随机行走。由于 Lévy 飞行具有无限的均值和方差,CS 可以比使用标准高斯过程的算法更有效地探索搜索空间。这一优势,加上局部搜索能力,保证全局收敛,使 CS 非常高效。事实上,各种研究和应用已经证明 CS 是非常有效的 [32,13,26,8]。

全局收敛:数学简析

Wang 等人为标准 CS 提供了全局收敛的数学证明,他们的方法基于马尔可夫链理论 [27]。他们的证明可以概括如下:

由于更新公式中有两个分支,局部搜索主要用于局部细化,而主要的探索是通过全局搜索来完成的。为了简化分析并强调全局搜索能力,我们使用简化版的 CS。也就是说,与发现/切换概率 相比, 我们只使用随机数 的全局分支。现在我们有,

因为 CS 算法是一个随机搜索算法,我们可以总结为以下主要步骤:

  1. 初始种群由 n 个随机位置的巢构成, ,然后评估它们的目标函数值,以找到当前的全局最佳
  2. 由下式更新新解/位置
  3. 从均匀分布 [0,1] 中取随机数 r。如果 ,则更新 。然后评估新解,找到新的全局最佳
  4. 如果满足停止要求,则 是目前为止发现的最好的全局解。否则,返回步骤(2)。

算法的全局收敛性。如果 是可测的,且可行解空间 上的一个可测子集,算法 A 满足前两个条件,搜索序列 ,则

也就是说,算法 A 可以以概率 1 全局收敛。这里 是第 k 次迭代时,在 上的第 k 个解的概率测度。

状态和状态空间。搜索历史中布谷鸟/巢的位置及其全局最优解 g 形成了布谷鸟的状态: ,其中 。所有可能状态的集合形成了状态空间,用下式表示

布谷鸟群/种群的状态和空间。所有 n 个布谷鸟/巢的状态组成群的状态,由 表示。所有布谷鸟的所有状态组成了群的状态空间,记为

显然,Q 包含整个种群的历史全局最优解 和历史上所有个体最优解 。此外,整个种群的全局最优解是所有 中最好的,所以

CS 中从状态 的转移概率是

其中 是 CS 中第 2 步的转移概率, 是该步骤中历史全局最优的转移概率。 是第 3 步的转移概率,而 是历史全局最佳。

对于优化问题 的全局最优解 ,最优状态集定义为

对于优化问题 的全局最优解 ,最优群状态集可定义为

所有这些集合将保证收敛条件得到满足。进一步详细的数学分析证明,当迭代次数接近足够大时,群状态序列将收敛到最优状态/解集 。因此,CS 能够保证全局收敛。

应用

CS 已经被应用于许多优化和计算智能领域。例如,在工程设计应用中,CS 在一系列连续优化问题(例如弹簧设计和焊接梁设计)上面比其他算法具有更好的性能 [32,13]。

另外,由 Walton 等人修改的 CS [26] 已被证明可以非常有效地解决非线性问题,如网格生成。 Yildiz [36] 利用 CS 在铣削加工中选择最优的机床参数,提高了结果,Zheng 和 Zhou [37] 提供了一个使用高斯过程的CS变体。

在数据融合和无线传感器网络中,CS 已被证明是非常有效的 [9,10]。此外,基于量子的 CS 能够有效地解决背包问题 [14]。从算法分析的角度来看,Civicioglu 和 Desdo [8] 提出的 CS 与粒子群优化(PSO),差分进化(DE)和人工蜂群(ABC)的概念比较表明,CS 与差分进化算法比 PSO 与 ABC 提供了更健壮的结果。Gandomi 等人 [13] 为解决各种结构优化问题,提供了更为广泛的比较研究,得出的结论是,CS 与其他算法如 PSO 和遗传算法(GA)相比能够获得更好的结果。在各种应用中有一些有趣的性能提升,如 Valian 等人使用 CS 来训练神经网络 [23],以及可靠性优化问题 [24]。

对于复杂的相平衡应用,Bhargava 等人 [1] 表明,CS 为解决热力学计算提供了一个可靠的方法。与此同时,Bulatovi’c 等人 [3] 利用 CS 解决了六杆双闭锁联锁问题,Moravej 和 Akhlaghi [16] 以良好收敛速度和性能解决了配电网络中的 DG 分配问题。

作为进一步的扩展,Yang 和 Deb [33] 针对设计工程应用提出了多目标布谷鸟搜索(MOCS)。对于多目标调度问题,Chandrasekaran 和 Simon [5] 使用 CS 算法取得了很大的进展,这证明了他们提出的方法的优越性。最近的研究表明,CS 在许多应用中可以比其他算法表现得更好 [13,17,37,36]。更详细的综述,请参阅 Yang [35] 和 Yang 等人 [31]。

源代码

Version 1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
% -----------------------------------------------------------------
% Cuckoo Search (CS) algorithm by Xin-She Yang and Suash Deb %
% Programmed by Xin-She Yang at Cambridge University %
% Programming dates: Nov 2008 to June 2009 %
% Last revised: Dec 2009 (simplified version for demo only) %
% -----------------------------------------------------------------
% Papers -- Citation Details:
% 1) X.-S. Yang, S. Deb, Cuckoo search via Levy flights,
% in: Proc. of World Congress on Nature & Biologically Inspired
% Computing (NaBIC 2009), December 2009, India,
% IEEE Publications, USA, pp. 210-214 (2009).
% http://arxiv.org/PS_cache/arxiv/pdf/1003/1003.1594v1.pdf
% 2) X.-S. Yang, S. Deb, Engineering optimization by cuckoo search,
% Int. J. Mathematical Modelling and Numerical Optimisation,
% Vol. 1, No. 4, 330-343 (2010).
% http://arxiv.org/PS_cache/arxiv/pdf/1005/1005.2908v2.pdf
% ----------------------------------------------------------------%
% This demo program only implements a standard version of %
% Cuckoo Search (CS), as the Levy flights and generation of %
% new solutions may use slightly different methods. %
% The pseudo code was given sequentially (select a cuckoo etc), %
% but the implementation here uses Matlab's vector capability, %
% which results in neater/better codes and shorter running time. %
% This implementation is different and more efficient than the %
% the demo code provided in the book by
% "Yang X. S., Nature-Inspired Metaheuristic Algoirthms, %
% 2nd Edition, Luniver Press, (2010). " %
% --------------------------------------------------------------- %

% =============================================================== %
% Notes: %
% Different implementations may lead to slightly different %
% behavour and/or results, but there is nothing wrong with it, %
% as this is the nature of random walks and all metaheuristics. %
% -----------------------------------------------------------------

function [bestnest,fmin]=cuckoo_search(n)
if nargin<1,
% Number of nests (or different solutions)
n=25;
end

% Discovery rate of alien eggs/solutions
pa=0.25;

%% Change this if you want to get better results
% Tolerance
Tol=1.0e-5;
%% Simple bounds of the search domain
% Lower bounds
nd=15;
Lb=-5*ones(1,nd);
% Upper bounds
Ub=5*ones(1,nd);

% Random initial solutions
for i=1:n,
nest(i,:)=Lb+(Ub-Lb).*rand(size(Lb));
end

% Get the current best
fitness=10^10*ones(n,1);
[fmin,bestnest,nest,fitness]=get_best_nest(nest,nest,fitness);

N_iter=0;
%% Starting iterations
while (fmin>Tol),

% Generate new solutions (but keep the current best)
new_nest=get_cuckoos(nest,bestnest,Lb,Ub);
[fnew,best,nest,fitness]=get_best_nest(nest,new_nest,fitness);
% Update the counter
N_iter=N_iter+n;
% Discovery and randomization
new_nest=empty_nests(nest,Lb,Ub,pa) ;

% Evaluate this set of solutions
[fnew,best,nest,fitness]=get_best_nest(nest,new_nest,fitness);
% Update the counter again
N_iter=N_iter+n;
% Find the best objective so far
if fnew<fmin,
fmin=fnew;
bestnest=best;
end
end %% End of iterations

%% Post-optimization processing
%% Display all the nests
disp(strcat('Total number of iterations=',num2str(N_iter)));
fmin
bestnest

%% --------------- All subfunctions are list below ------------------
%% Get cuckoos by ramdom walk
function nest=get_cuckoos(nest,best,Lb,Ub)
% Levy flights
n=size(nest,1);
% Levy exponent and coefficient
% For details, see equation (2.21), Page 16 (chapter 2) of the book
% X. S. Yang, Nature-Inspired Metaheuristic Algorithms, 2nd Edition, Luniver Press, (2010).
beta=3/2;
sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);

for j=1:n,
s=nest(j,:);
% This is a simple way of implementing Levy flights
% For standard random walks, use step=1;
%% Levy flights by Mantegna's algorithm
u=randn(size(s))*sigma;
v=randn(size(s));
step=u./abs(v).^(1/beta);

% In the next equation, the difference factor (s-best) means that
% when the solution is the best solution, it remains unchanged.
stepsize=0.01*step.*(s-best);
% Here the factor 0.01 comes from the fact that L/100 should the typical
% step size of walks/flights where L is the typical lenghtscale;
% otherwise, Levy flights may become too aggresive/efficient,
% which makes new solutions (even) jump out side of the design domain
% (and thus wasting evaluations).
% Now the actual random walks or flights
s=s+stepsize.*randn(size(s));
% Apply simple bounds/limits
nest(j,:)=simplebounds(s,Lb,Ub);
end

%% Find the current best nest
function [fmin,best,nest,fitness]=get_best_nest(nest,newnest,fitness)
% Evaluating all new solutions
for j=1:size(nest,1),
fnew=fobj(newnest(j,:));
if fnew<=fitness(j),
fitness(j)=fnew;
nest(j,:)=newnest(j,:);
end
end
% Find the current best
[fmin,K]=min(fitness) ;
best=nest(K,:);

%% Replace some nests by constructing new solutions/nests
function new_nest=empty_nests(nest,Lb,Ub,pa)
% A fraction of worse nests are discovered with a probability pa
n=size(nest,1);
% Discovered or not -- a status vector
K=rand(size(nest))>pa;

% In the real world, if a cuckoo's egg is very similar to a host's eggs, then
% this cuckoo's egg is less likely to be discovered, thus the fitness should
% be related to the difference in solutions. Therefore, it is a good idea
% to do a random walk in a biased way with some random step sizes.
%% New solution by biased/selective random walks
stepsize=rand*(nest(randperm(n),:)-nest(randperm(n),:));
new_nest=nest+stepsize.*K;
for j=1:size(new_nest,1)
s=new_nest(j,:);
new_nest(j,:)=simplebounds(s,Lb,Ub);
end

% Application of simple constraints
function s=simplebounds(s,Lb,Ub)
% Apply the lower bound
ns_tmp=s;
I=ns_tmp<Lb;
ns_tmp(I)=Lb(I);

% Apply the upper bounds
J=ns_tmp>Ub;
ns_tmp(J)=Ub(J);
% Update this new move
s=ns_tmp;

%% You can replace the following by your own functions
% A d-dimensional objective function
function z=fobj(u)
%% d-dimensional sphere function sum_j=1^d (u_j-1)^2.
% with a minimum at (1,1, ...., 1);
z=sum((u-1).^2);

Version 2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
% -----------------------------------------------------------------
% Cuckoo Search (CS) algorithm by Xin-She Yang and Suash Deb %
% Programmed by Xin-She Yang at Cambridge University %
% Programming dates: Nov 2008 to June 2009 %
% Last revised: Dec 2009 (simplified version for demo only) %
% -----------------------------------------------------------------
% Papers -- Citation Details:
% 1) X.-S. Yang, S. Deb, Cuckoo search via Levy flights,
% in: Proc. of World Congress on Nature & Biologically Inspired
% Computing (NaBIC 2009), December 2009, India,
% IEEE Publications, USA, pp. 210-214 (2009).
% http://arxiv.org/PS_cache/arxiv/pdf/1003/1003.1594v1.pdf
% 2) X.-S. Yang, S. Deb, Engineering optimization by cuckoo search,
% Int. J. Mathematical Modelling and Numerical Optimisation,
% Vol. 1, No. 4, 330-343 (2010).
% http://arxiv.org/PS_cache/arxiv/pdf/1005/1005.2908v2.pdf
% ----------------------------------------------------------------%
% This demo program only implements a standard version of %
% Cuckoo Search (CS), as the Levy flights and generation of %
% new solutions may use slightly different methods. %
% The pseudo code was given sequentially (select a cuckoo etc), %
% but the implementation here uses Matlab's vector capability, %
% which results in neater/better codes and shorter running time. %
% This implementation is different and more efficient than the %
% the demo code provided in the book by
% "Yang X. S., Nature-Inspired Metaheuristic Algoirthms, %
% 2nd Edition, Luniver Press, (2010). " %
% --------------------------------------------------------------- %

% =============================================================== %
% Notes: %
% Different implementations may lead to slightly different %
% behavour and/or results, but there is nothing wrong with it, %
% as this is the nature of random walks and all metaheuristics. %
% -----------------------------------------------------------------

% Additional Note: This version uses a fixed number of generation %
% (not a given tolerance) because many readers asked me to add %
% or implement this option. Thanks.%
function [bestnest,fmin]=cuckoo_search_new(n)
if nargin<1,
% Number of nests (or different solutions)
n=25;
end

% Discovery rate of alien eggs/solutions
pa=0.25;

%% Change this if you want to get better results
N_IterTotal=1000;
%% Simple bounds of the search domain
% Lower bounds
nd=15;
Lb=-5*ones(1,nd);
% Upper bounds
Ub=5*ones(1,nd);

% Random initial solutions
for i=1:n,
nest(i,:)=Lb+(Ub-Lb).*rand(size(Lb));
end

% Get the current best
fitness=10^10*ones(n,1);
[fmin,bestnest,nest,fitness]=get_best_nest(nest,nest,fitness);

N_iter=0;
%% Starting iterations
for iter=1:N_IterTotal,
% Generate new solutions (but keep the current best)
new_nest=get_cuckoos(nest,bestnest,Lb,Ub);
[fnew,best,nest,fitness]=get_best_nest(nest,new_nest,fitness);
% Update the counter
N_iter=N_iter+n;
% Discovery and randomization
new_nest=empty_nests(nest,Lb,Ub,pa) ;

% Evaluate this set of solutions
[fnew,best,nest,fitness]=get_best_nest(nest,new_nest,fitness);
% Update the counter again
N_iter=N_iter+n;
% Find the best objective so far
if fnew<fmin,
fmin=fnew;
bestnest=best;
end
end %% End of iterations

%% Post-optimization processing
%% Display all the nests
disp(strcat('Total number of iterations=',num2str(N_iter)));
fmin
bestnest

%% --------------- All subfunctions are list below ------------------
%% Get cuckoos by ramdom walk
function nest=get_cuckoos(nest,best,Lb,Ub)
% Levy flights
n=size(nest,1);
% Levy exponent and coefficient
% For details, see equation (2.21), Page 16 (chapter 2) of the book
% X. S. Yang, Nature-Inspired Metaheuristic Algorithms, 2nd Edition, Luniver Press, (2010).
beta=3/2;
sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);

for j=1:n,
s=nest(j,:);
% This is a simple way of implementing Levy flights
% For standard random walks, use step=1;
%% Levy flights by Mantegna's algorithm
u=randn(size(s))*sigma;
v=randn(size(s));
step=u./abs(v).^(1/beta);

% In the next equation, the difference factor (s-best) means that
% when the solution is the best solution, it remains unchanged.
stepsize=0.01*step.*(s-best);
% Here the factor 0.01 comes from the fact that L/100 should the typical
% step size of walks/flights where L is the typical lenghtscale;
% otherwise, Levy flights may become too aggresive/efficient,
% which makes new solutions (even) jump out side of the design domain
% (and thus wasting evaluations).
% Now the actual random walks or flights
s=s+stepsize.*randn(size(s));
% Apply simple bounds/limits
nest(j,:)=simplebounds(s,Lb,Ub);
end

%% Find the current best nest
function [fmin,best,nest,fitness]=get_best_nest(nest,newnest,fitness)
% Evaluating all new solutions
for j=1:size(nest,1),
fnew=fobj(newnest(j,:));
if fnew<=fitness(j),
fitness(j)=fnew;
nest(j,:)=newnest(j,:);
end
end
% Find the current best
[fmin,K]=min(fitness) ;
best=nest(K,:);

%% Replace some nests by constructing new solutions/nests
function new_nest=empty_nests(nest,Lb,Ub,pa)
% A fraction of worse nests are discovered with a probability pa
n=size(nest,1);
% Discovered or not -- a status vector
K=rand(size(nest))>pa;

% In the real world, if a cuckoo's egg is very similar to a host's eggs, then
% this cuckoo's egg is less likely to be discovered, thus the fitness should
% be related to the difference in solutions. Therefore, it is a good idea
% to do a random walk in a biased way with some random step sizes.
%% New solution by biased/selective random walks
stepsize=rand*(nest(randperm(n),:)-nest(randperm(n),:));
new_nest=nest+stepsize.*K;
for j=1:size(new_nest,1)
s=new_nest(j,:);
new_nest(j,:)=simplebounds(s,Lb,Ub);
end

% Application of simple constraints
function s=simplebounds(s,Lb,Ub)
% Apply the lower bound
ns_tmp=s;
I=ns_tmp<Lb;
ns_tmp(I)=Lb(I);

% Apply the upper bounds
J=ns_tmp>Ub;
ns_tmp(J)=Ub(J);
% Update this new move
s=ns_tmp;

%% You can replace the following by your own functions
% A d-dimensional objective function
function z=fobj(u)
%% d-dimensional sphere function sum_j=1^d (u_j-1)^2.
% with a minimum at (1,1, ...., 1);
z=sum((u-1).^2);

Version 3

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
% ======================================================== % 
% Files of the Matlab programs included in the book: %
% Xin-She Yang, Nature-Inspired Metaheuristic Algorithms, %
% Second Edition, Luniver Press, (2010). www.luniver.com %
% ======================================================== %


% Cuckoo Search for nonlinear constrained optimization %
% Programmed by Xin-She Yang @ Cambridge University 2009 %
% Usage: cuckoo_spring(25000) or cuckoo_spring; %

function [bestsol,fval]=cuckoo_search_spring(time)
format long;
help cuckoo_search_spring.m
if nargin<1,
% Number of iteraions
time=2000;
end

disp('Computing ... it may take a few minutes.');

% Number of nests (or different solutions)
n=25;
% Discovery rate of alien eggs/solutions
pa=0.25;

% Simple bounds of the search domain
% Lower bounds and upper bounds
Lb=[0.05 0.25 2.0];
Ub=[2.0 1.3 15.0];


% Random initial solutions
for i=1:n,
nest(i,:)=Lb+(Ub-Lb).*rand(size(Lb));
end

% Get the current best
fitness=10^10*ones(n,1);
[fmin,bestnest,nest,fitness]=get_best_nest(nest,nest,fitness);

N_iter=0;
%% Starting iterations
for t=1:time,

% Generate new solutions (but keep the current best)
new_nest=get_cuckoos(nest,bestnest,Lb,Ub);
[fnew,best,nest,fitness]=get_best_nest(nest,new_nest,fitness);
% Update the counter
N_iter=N_iter+n;
% Discovery and randomization
new_nest=empty_nests(nest,Lb,Ub,pa) ;

% Evaluate this solution
[fnew,best,nest,fitness]=get_best_nest(nest,new_nest,fitness);
% Update the counter again
N_iter=N_iter+n;
% Find the best objective so far
if fnew<fmin,
fmin=fnew;
bestnest=best ;
end
end %% End of iterations

%% Post-optimization processing
%% Display all the nests
disp(strcat('Total number of iterations=',num2str(N_iter)));
fmin
bestnest

%% --------------- All subfunctions are list below ------------------
%% Get cuckoos by ramdom walk
function nest=get_cuckoos(nest,best,Lb,Ub)
% Levy flights
n=size(nest,1);
% Levy exponent and coefficient
% For details, see equation (2.21), Page 16 (chapter 2) of the book
% X. S. Yang, Nature-Inspired Metaheuristic Algorithms, 2nd Edition, Luniver Press, (2010).
beta=3/2;
sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);

for j=1:n,
s=nest(j,:);
% This is a simple way of implementing Levy flights
% For standard random walks, use step=1;
%% Levy flights by Mantegna's algorithm
u=randn(size(s))*sigma;
v=randn(size(s));
step=u./abs(v).^(1/beta);

% In the next equation, the difference factor (s-best) means that
% when the solution is the best solution, it remains unchanged.
stepsize=0.01*step.*(s-best);
% Here the factor 0.01 comes from the fact that L/100 should the typical
% step size of walks/flights where L is the typical lenghtscale;
% otherwise, Levy flights may become too aggresive/efficient,
% which makes new solutions (even) jump out side of the design domain
% (and thus wasting evaluations).
% Now the actual random walks or flights
s=s+stepsize.*randn(size(s));
% Apply simple bounds/limits
nest(j,:)=simplebounds(s,Lb,Ub);
end


%% Find the current best nest
function [fmin,best,nest,fitness]=get_best_nest(nest,newnest,fitness)
% Evaluating all new solutions
for j=1:size(nest,1),
fnew=fobj(newnest(j,:));
if fnew<=fitness(j),
fitness(j)=fnew;
nest(j,:)=newnest(j,:);
end
end
% Find the current best
[fmin,K]=min(fitness) ;
best=nest(K,:);

%% Replace some nests by constructing new solutions/nests
function new_nest=empty_nests(nest,Lb,Ub,pa)
% A fraction of worse nests are discovered with a probability pa
n=size(nest,1);
% Discovered or not -- a status vector
K=rand(size(nest))>pa;

% In real world, if a cuckoo's egg is very similar to a host's eggs, then
% this cuckoo's egg is less likely to be discovered, thus the fitness should
% be related to the difference in solutions. Therefore, it is a good idea
% to do a random walk in a biased way with some random step sizes.
nestn1=nest(randperm(n),:);
nestn2=nest(randperm(n),:);
%% New solution by biased/selective random walks
stepsize=rand*(nestn1-nestn2);
new_nest=nest+stepsize.*K;
for j=1:size(new_nest,1)
s=new_nest(j,:);
new_nest(j,:)=simplebounds(s,Lb,Ub);
end

% Application of simple constraints
function s=simplebounds(s,Lb,Ub)
% Apply the lower bound
ns_tmp=s;
I=ns_tmp<Lb;
ns_tmp(I)=Lb(I);

% Apply the upper bounds
J=ns_tmp>Ub;
ns_tmp(J)=Ub(J);
% Update this new move
s=ns_tmp;

%% Spring desgn optimization -- objective function
function z=fobj(u)
% The well-known spring design problem
z=(2+u(3))*u(1)^2*u(2);
z=z+getnonlinear(u);

function Z=getnonlinear(u)
Z=0;
% Penalty constant
lam=10^15;

% Inequality constraints
g(1)=1-u(2)^3*u(3)/(71785*u(1)^4);
gtmp=(4*u(2)^2-u(1)*u(2))/(12566*(u(2)*u(1)^3-u(1)^4));
g(2)=gtmp+1/(5108*u(1)^2)-1;
g(3)=1-140.45*u(1)/(u(2)^2*u(3));
g(4)=(u(1)+u(2))/1.5-1;

% No equality constraint in this problem, so empty;
geq=[];

% Apply inequality constraints
for k=1:length(g),
Z=Z+ lam*g(k)^2*getH(g(k));
end
% Apply equality constraints
for k=1:length(geq),
Z=Z+lam*geq(k)^2*getHeq(geq(k));
end

% Test if inequalities hold
% Index function H(g) for inequalities
function H=getH(g)
if g<=0,
H=0;
else
H=1;
end
% Index function for equalities
function H=getHeq(geq)
if geq==0,
H=0;
else
H=1;
end
% ----------------- end ------------------------------

参考文献

[1] Bhargava V, Fateen SEK, Bonilla-Petriciolet A. Cuckoo search: a new nature-inspired optimization method for phase equilibrium calculations. Fluid Phase Equilib 2013;337(1):191–200.

[2] Barthelemy P, Bertolotti J, Wiersma DS. A Lévy flight for light. Nature 2008;453(6948):495–8.

[3] Bulatovi´ c RR, Bordevi´ c SR, Dordevi´ c VS. Cuckoo search algorithm: a metaheuristic approach to solving the problem of optimum synthesis of a six-bar double dwell linkage. Mech Mach Theory 2013;61(1):1–3.

[4] Brown C, Liebovitch LS, Glendon R. Lévy flights in Dobe Ju/’hoansi foraging patterns. Human Ecol 2007;35(2):129–38.

[5] Chandrasekaran K, Simon SP. Multi-objective scheduling problem: hybrid appraoch using fuzzy assisted cuckoo search algorithm. Swarm Evol Comput 2012;5(1):1–6.

[6] Chifu VR, Pop CB, Salomie I, Suia DS, Niculici AN. Optimizing the semantic web service composition process using cuckoo search. In: Brazier FMT, Nieuwenhuis K, Palvin G, Warnier M, Badica C, editors. Intelligent distributed computing V. Studies in computational intelligence, vol. 382. Berlin, Germany: Springer; 2012. p. 93–102.

[7] Choudhary K, Purohit GN. A new testing approach using cuckoo search to achieve multiobjective genetic algorithm. J Comput 2011;3(4):117–9.

[8] Civicioglu P, Besdok E. A conception comparison of the cuckoo search, particle swarm optimization, differential evolution and artificial bee colony algorithms. Artif Intell Rev 2013;39(3):315–46.

[9] Dhivya M, Sundarambal M, Anand LN. Energy efficient computation of data fusion in wireless sensor networks using cuckoo-based particle approach (CBPA). Int J of Commun Netw Syst Sci 2011;4(4):249–55.

[10] Dhivya M, Sundarambal M. Cuckoo search for data gathering in wireless sensor networks. Int J Mobile Commun 2011;9(6):642–56.

[11] Durgun I, Yildiz AR. Structural design optimization of vehicle components using cuckoo search algorithm. Mater Test 2012;3(2):185–8.

[12] Fister JrI, Yang XS, Fister I, Fister D. Cuckoo search: a brief literature review. In: Yang XS, editor. Cuckoo search and firefly algorithm: theory and applications. Heidelberg, Germany: Springer; 2013 [chapter 3].

[13] Gandomi AH, Yang XS, Alavi AH. Cuckoo search algorithm: a meteheuristic approach to solve structural optimization problems. Eng Comput 2013;29(1):17–35.

[14] Layeb A. A novel quantum-inspired cuckoo search for Knapsack problems. Int J Bioinspired Comput 2011;3(5):297–305.

[15] Mantegna RN. Fast, accurate algorithm for numerical simulation of Lévy stable stochastic process. Phys Rev E 1994;49(5):4677–83.

[16] Moravej Z, Akhlaghi A. A novel approach based on cuckoo search for DG allocation in distribution network. Electr Power Energy Syst 2013;44(1):672–9.

[17] Noghrehabadi A, Ghalambaz M, Vosough A. A hybrid power series: cuckoo search optimization algorithm to electrostatic deflection of micro fixed-fixed actuators. Int J Multidisciplinary Sci Eng 2011;2(4):22–6.

[18] Ouaarab A, Ahiod B, Yang XS. Discrete cuckoo search algorithm for the travelling salesman problem, Neural Computing and Applications; in press. http://dx.doi.org/10.1007/s00521-013-1402-2.

[19] Payne RB, Sorenson MD, Klitz K. The cuckoos. Oxford, UK: Oxford University Press; 2005.

[20] Pavlyukevich I. Lévy flights, non-local search and simulated annealing. J Comput Phys 2007;226(2):1830–44.

[21] Reynolds AM, Frye MA. Free-flight odor tracking in Drosophila is consistent with an optimal intermittent scale-free search. PLoS One 2007;2(4):e354–63.

[22] Salimi H, Giveki D, Soltanshahi MA, Hatami J. Extended mixture of MLP experts by hybrid of conjugate gradient method and modified cuckoo search. Int J Artif Intell Appl 2012;3(1):1–3.

[23] Valian E, Mohanna S, Tavakoli S. Improved cuckoo search algorithm for feedforward neural network training. Int J Artif Intell Appl 2011;2(3):36–43.

[24] Valian E, Tavakoli S, Mohanna S, Haghi A. Improved cuckoo search for reliability optimization problems. Comput Ind Eng 2013;64(1):459–68.

[25] Vazquez RA. Training spiking neural models using cuckoo search algorithm. In: IEEE congress on evolutionary computation. New Orleans: IEEE publication; 2012.p. 679–86.

[26] Walton S, Hassan O, Morgan K, Brown MR. Modified cuckoo search: a new gradient free optimization algorithm. Chaos Solitons Fractals 2011;44(9):710–8.

[27] Wang F, He X-S, Wang Y, Yang SM. Markov model and convergence analysis based on cuckoo search algorithm. Comput Eng 2012;38(11):180–5.

[28] Wang F, Lou L, He X, Wang Y. Hybrid optimization algorithm of PSO and cuckoo search. In: Proceedings of the second international conference on artificial intelligence, management science and electronic commerce(AIMSEC’11). Zhengzhou:IEEE publication; 2011. p. 1172–5.

[29] Yang XS, Deb S, Fong S. Accelerated particle swarm optimization and support vector machine for business optimization and applications. In: Networked digital technologies 2011. Communications in computer and information science, vol. 136. Berlin, Germany: Springer; 2011. p. 53–66.

[30] Yang XS, Deb S. Cuckoo search via Lévy flights. In: Proceeings of world congress on nature & biologically inspired computing (NaBIC 2009). USA: IEEE Publications; 2009. p. 210–4.

[31] Yang XS, Cui ZH, Xiao RB, Gandomi AH, Karamanoglu M. Swarm intelligence and bio-inspired computation: theory and applications. Waltham, MA, USA: Elsevier; 2013.

[32] Yang XS, Deb S. Engineering optimization by cuckoo search. Int J Math Model Numer Optimisation 2010;1(4):330–43.

[33] Yang XS, Deb S. Multiobjective cuckoo search for design optimization. Comput Oper Res 2013;40(6):1616–24.

[34] Yang XS, Deb S. Cuckoo search: recent advances and applications, Neural Computing and Applications; in press. http://link.springer.com/article/10.1007.

[35] Yang XS. Cuckoo search and firefly algorithm: theory and applications. Studies in computational intelligence, vol. 516. Heidelberg, Germany: Springer; 2013.

[36] Yildiz AR. Cuckoo search algorithm for the selection of optimal machine parameters in milling operations. Int J Adv Manuf Technol. 2013;64(1–4):55–61.

[37] Zheng HQ, Zhou Y. A novel cuckoo search optimization algorithm based on Gauss distribution. J Comput Inform Syst 2012;8(10):4193–200.

本文翻译自 Yang X S. Nature-inspired optimization algorithms[M]. Elsevier, 2014. 有删节

GreatX wechat
订阅公众号,获取更多信息。