今天我要讲的内容是布谷鸟算法，英文叫做Cuckoo search (CS algorithm)。

首先还是同样，介绍一下这个算法的英文含义，Cuckoo是布谷鸟的意思，啥是布谷鸟呢，是一种叫做布谷的鸟，o(∩_∩)o ，这种鸟她妈很懒，自己生蛋自己不养，一般把它的宝宝扔到别的种类鸟的鸟巢去。

但是呢，当孵化后，遇到聪明的鸟妈妈，一看就知道不是亲生的，直接就被鸟妈妈给杀了。

于是这群布谷鸟宝宝为了保命，它们就模仿别的种类的鸟叫，让智商或者情商极低的鸟妈妈误认为是自己的亲宝宝，这样它就活下来了。

Search指的是搜索，这搜索可不是谷歌一下，你就知道。

而是搜索最优值，举个简单的例子，y=(x-0.5)^2+1，它的最小值是1，位置是(0.5,1)，我们要搜索的就是这个位置。

现在我们应该清楚它是干嘛的了吧，它就是为了寻找最小值而产生的一种算法，有些好装X的人会说，你傻X啊，最小值不是-2a/b吗，用你找啊? 说的不错，确实是，但是要是我们的函数变成y=sin(x^3+x^2)+e^cos(x^3)+log(tan(x) +10，你怎么办吶？你解不了，就算你求导数，但是你知道怎么解导数等于0吗？所以我们就得引入先进的东西来求最小值。

为了使大家容易理解，我还是用y=(x-0.5)^2+1来举例子，例如我们有4个布谷鸟蛋（也就是4个x坐标），鸟妈妈发现不是自己的宝宝的概率是0.25，我们x的取值范围是[0,1]之间，于是我们就可以开始计算了。

目标：求x在[0,1]之内的函数y=(x-0.5)^2+1最小值（1）初始化x的位置，随机生成4个x坐标，x1=0.4，x2=0.6，x3=0.8，x4 =0.3 ——> X=[0.4, 0.6 ,0.8, 0.3]（2）求出y1~y4，把x1~x4带入函数，求得Y=[1,31, 1.46, 1.69, 1.265]，并选取当前最小值ymin= y4=1.265（3）开始定出一个y的最大值为Y_global=INF（无穷大），然后与ymin比较，把Y中最小的位置和值保留，例如Y_global=INF>ymin=1.265，所以令Y _global=1.265（4）记录Y_global的位置，（0.3，1.265）。

（5）按概率0.25，随机地把X中的值过塞子，选出被发现的蛋。

例如第二个蛋被发现x2=0.6，那么他就要随机地变换位子，生成一个随机数，例如0.02，然后把x2=x2+0.02=0.62，之后求出y2=1.4794。

那么X就变为了X=[0.4, 0.6 2 ,0.8, 0.3]，Y=[1,31, 1.4794, 1.69, 1.265]。

（6）进行莱维飞行，这名字听起来挺高大上，说白了，就是把X的位置给随机地改变了。

怎么变？有一个公式x=x+alpha*L。

L=S*(X-Y_global)*rand3S=[rand1*sigma/|rand2|]^(1/beta)sigma=0.6966beta=1.5alpha=0.01rand1~rand3为正态分布的随机数然后我们把X=[0.4, 0.6 ,0.8, 0.3]中的x1带入公式，首先随机生成rand1=-1. 2371，rand2=-2.1935，rand3=-0.3209，接下来带入公式中，获得x1=0.39 85之后同理计算:x2=0.6172x3=0.7889x4=0.3030（7）更新矩阵X，X=[0.3985, 0.6172, 0.7889, 0.3030]（8）计算Y=[1.3092, 1.4766, 1.6751, 1.2661]，并选取当前最小值ymin= y4=1.2661，然后与ymin比较，把Y中最小的位置和值保留，例如Y_global =1.265<ymin=1.2661，所以令Y_global=1.265（9）返回步骤（5）用更新的X去循环执行，经过多次计算即可获得y的最优值和的最值位置(x,y)1.% -----------------------------------------------------------------2.% Cuckoo Search (CS) algorithm by Xin-She Yang and Suash Deb %3.% Programmed by Xin-She Yang at Cambridge University %4.% Programming dates: Nov 2008 to June 2009 %5.% Last revised: Dec 2009 (simplified version for demo only) %6.% -----------------------------------------------------------------7.% Papers -- Citation Details:8.% 1) X.-S. Yang, S. Deb, Cuckoo search via Levy flights,9.% in: Proc. of World Congress on Nature & Biologically Inspired10.% Computing (NaBIC 2009), December 2009, India,11.% IEEE Publications, USA, pp. 210-214 (2009).12.% /PS_cache/arxiv/pdf/1003/1003.1594v1.pdf13.% 2) X.-S. Yang, S. Deb, Engineering optimization by cuckoo search,14.% Int. J. Mathematical Modelling and Numerical Optimisation,15.% Vol. 1, No. 4, 330-343 (2010).16.% /PS_cache/arxiv/pdf/1005/1005.2908v2.pdf17.% ----------------------------------------------------------------%18.% This demo program only implements a standard version of %19.% Cuckoo Search (CS), as the Levy flights and generation of %20.% new solutions may use slightly different methods. %21.% The pseudo code was given sequentially (select a cuckoo etc), %22.% but the implementation here uses Matlab's vector capability, %23.% which results in neater/better codes and shorter running time. %24.% This implementation is different and more efficient than the %25.% the demo code provided in the book by26.% "Yang X. S., Nature-Inspired Metaheuristic Algoirthms, %27.% 2nd Edition, Luniver Press, (2010). " %28.% --------------------------------------------------------------- %29.30.% =============================================================== %31.% Notes: %32.% Different implementations may lead to slightly different %33.% behavour and/or results, but there is nothing wrong with it, %34.% as this is the nature of random walks and all metaheuristics. %35.% -----------------------------------------------------------------36.37.% Additional Note: This version uses a fixed number of generation %38.% (not a given tolerance) because many readers asked me to add %39.% or implement this option. Thanks.%40.function [bestnest,fmin]=cuckoo_search_new(n)41.if nargin<1,42.% Number of nests (or different solutions)43.n=25;44.end45.46.% Discovery rate of alien eggs/solutions47.pa=0.25;48.49.%% Change this if you want to get better results50.N_IterTotal=1000;51.%% Simple bounds of the search domain52.% Lower bounds53.nd=15;54.Lb=-5*ones(1,nd);55.% Upper bounds56.Ub=5*ones(1,nd);57.58.% Random initial solutions59.for i=1:n,60.nest(i,:)=Lb+(Ub-Lb).*rand(size(Lb));61.end62.63.% Get the current best64.fitness=10^10*ones(n,1);65.[fmin,bestnest,nest,fitness]=get_best_nest(nest,nest,fitness);66.67.N_iter=0;68.%% Starting iterations69.for iter=1:N_IterTotal,70. % Generate new solutions (but keep the current best)71. new_nest=get_cuckoos(nest,bestnest,Lb,Ub);72. [fnew,best,nest,fitness]=get_best_nest(nest,new_nest,fitness);73. % Update the counter74. N_iter=N_iter+n;75. % Discovery and randomization76. new_nest=empty_nests(nest,Lb,Ub,pa) ;77.78. % Evaluate this set of solutions79. [fnew,best,nest,fitness]=get_best_nest(nest,new_nest,fitness);80. % Update the counter again81. N_iter=N_iter+n;82. % Find the best objective so far83. if fnew<fmin,84. fmin=fnew;85. bestnest=best;86. end87.end %% End of iterations88.89.%% Post-optimization processing90.%% Display all the nests91.disp(strcat('Total number of iterations=',num2str(N_iter)));92.fmin93.bestnest94.95.%% --------------- All subfunctions are list below ------------------96.%% Get cuckoos by ramdom walk97.function nest=get_cuckoos(nest,best,Lb,Ub)98.% Levy flights99.n=size(nest,1);100.% Levy exponent and coefficient101.% For details, see equation (2.21), Page 16 (chapter 2) of the book 102.% X. S. Yang, Nature-Inspired Metaheuristic Algorithms, 2nd Edition, Lunive r Press, (2010).103.beta=3/2;104.sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2))) ^(1/beta);105.106.for j=1:n,107. s=nest(j,:);108. % This is a simple way of implementing Levy flights109. % For standard random walks, use step=1;110. %% Levy flights by Mantegna's algorithm111. u=randn(size(s))*sigma;112. v=randn(size(s));113. step=u./abs(v).^(1/beta);114.115. % In the next equation, the difference factor (s-best) means that 116. % when the solution is the best solution, it remains unchanged. 117. stepsize=0.01*step.*(s-best);118. % Here the factor 0.01 comes from the fact that L/100 should the typica l119. % step size of walks/flights where L is the typical lenghtscale; 120. % otherwise, Levy flights may become too aggresive/efficient,121. % which makes new solutions (even) jump out side of the design domain 122. % (and thus wasting evaluations).123. % Now the actual random walks or flights124. s=s+stepsize.*randn(size(s));125. % Apply simple bounds/limits126. nest(j,:)=simplebounds(s,Lb,Ub);127.end128.129.%% Find the current best nest130.function [fmin,best,nest,fitness]=get_best_nest(nest,newnest,fitness) 131.% Evaluating all new solutions132.for j=1:size(nest,1),133. fnew=fobj(newnest(j,:));134. if fnew<=fitness(j),135. fitness(j)=fnew;136. nest(j,:)=newnest(j,:);137. end138.end139.% Find the current best140.[fmin,K]=min(fitness) ;141.best=nest(K,:);142.143.%% Replace some nests by constructing new solutions/nests144.function new_nest=empty_nests(nest,Lb,Ub,pa)145.% A fraction of worse nests are discovered with a probability pa146.n=size(nest,1);147.% Discovered or not -- a status vector148.K=rand(size(nest))>pa;149.150.% In the real world, if a cuckoo's egg is very similar to a host's eggs, the n151.% this cuckoo's egg is less likely to be discovered, thus the fitness should152.% be related to the difference in solutions. Therefore, it is a good idea153.% to do a random walk in a biased way with some random step sizes. 154.%% New solution by biased/selective random walks155.stepsize=rand*(nest(randperm(n),:)-nest(randperm(n),:));156.new_nest=nest+stepsize.*K;157.for j=1:size(new_nest,1)158. s=new_nest(j,:);159. new_nest(j,:)=simplebounds(s,Lb,Ub);160.end161.162.% Application of simple constraints163.function s=simplebounds(s,Lb,Ub)164. % Apply the lower bound165. ns_tmp=s;166. I=ns_tmp<Lb;167. ns_tmp(I)=Lb(I);168.169. % Apply the upper bounds170. J=ns_tmp>Ub;171. ns_tmp(J)=Ub(J);172. % Update this new move173. s=ns_tmp;174.175.%% You can replace the following by your own functions 176.% A d-dimensional objective function177.function z=fobj(u)178.%% d-dimensional sphere function sum_j=1^d (u_j-1)^2. 179.% with a minimum at (1,1, ...., 1);180.z=sum((u-1).^2);。