智能优化算法学习第一期：布谷鸟搜索算法（CS）【附MATLAB代码】

一、引言

本期我们介绍布谷鸟搜索(Cuckoo Search ,CS)算法。该算法由Yang等人于2009年提出，
主要模拟布谷鸟的繁衍自然行为。布谷鸟繁衍特点是寄生，将卵产在其他鸟类的巢穴中。
一般情况下，布谷鸟蛋的孵化速度比寄主蛋快。一旦第一只布谷鸟雏鸟孵化，它会迅速
推出其他蛋，增加寄主对它的食物供给。研究还指出，杜鹃雏鸟可以模仿寄主雏鸟的叫
声，以获取更多被喂食的机会。（适应度越高，蛋获得能量越高）

二、算法基本过程

流程图

①群体位置初始化

②局部探索

③全局探索

Lévy 飞行生成随机数应包括两个步骤：随机方向的选择和服从 Lévy 分布的步长的生成。方向服从均匀分布，步长由Mantegna算法给出：

三、算法代码

主程序 main.m

clear;clc;close all
fobj = @fun;        %适应度值
lb = -5;            %参数下限
ub = 20;            %参数上限
dim = 2;            %待求解参数维度
pop_num=100;        %种群数量
Max_iter=1000;      %最大迭代次数

t1=clock;
[Best_Score,Best_Pos,Convergence_curve] = CS(pop_num,Max_iter,lb,ub,dim,fobj); %寻找最小值
t2=clock;
time_CS=(t2(end)+t2(end-1)*60+t2(end-2)*3600-t1(end)-t1(end-1)*60-t1(end-2)*3600); %计算运行时间

%最优适应度变化曲线图
figure(1);  
plot(Convergence_curve,'r-','linewidth',2);
grid on;
title('布谷鸟算法迭代曲线');
xlabel('迭代次数');
ylabel('适应度值');

%…………………………………… 结果显示……………………………………
disp(['求解得到的x1，x2是:',num2str(Best_Pos(1)),' ',num2str(Best_Pos(2))]);
disp(['最优解:',num2str(Best_Score)]);

%待求解函数最小值(x2-x1^2*5.1/(4*pi^2)+5/pi*x1-6)^2+10*(1-1/(8*pi)*cos(x1))+10
function o = fun(x)
    o=(x(2)-(x(1)^2)*5.1/(4*(pi^2))+5/pi*x(1)-6)^2+10*(1-1/(8*pi))*cos(x(1))+10;
end

CS.m

% 布谷鸟优化算法
% 标准布谷鸟算法函数
% Leader_score：最优适应度值 Leader_pos：最优位置 Convergence_curve：最优适应度变化曲线
% SearchAgents_no：种群个数 Max_iter：进化次数 lb：变量下限 ub：变量上限 dim：待优化参数个数 fobj：适应度函数
function [Leader_score,Leader_pos,Convergence_curve] = CS(SearchAgents_no,Max_iter,lb,ub,dim,fobj)

%参数初始化
Pa = 0.25;                                                         %鸟蛋被发现的概率
alpha0 = 0.01;                                                     %比例因子
%种群初始化
Positions = initialization(SearchAgents_no,dim,ub,lb);
fitness = zeros(1,SearchAgents_no);

%计算适应度
for i = 1:SearchAgents_no
    fitness(i) = fobj(Positions(i,:));
end

%初始最优个体位置与适应度
[min_score, index]=min(fitness);
Leader_score = min_score;
Leader_pos = Positions(index, :);
Convergence_curve = zeros(1,Max_iter);                             %初始化最优适应度值变化曲线

%主函数(迭代更新)
for t = 1:Max_iter
    
    %莱维飞行种群更新（全局搜索）
    for i = 1:size(Positions,1)                                    %种群个数    
        alpha = alpha0.*(Positions(i,:)-Leader_pos);               %步长，等式4    
        newpop = Positions(i,:)+alpha.*Levy(dim);                  %等式1
        
        %超限判断
        Flag4ub = newpop>ub;
        Flag4lb = newpop<lb;
        newpop = (newpop.*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
        
        %计算新个体适应度值 
        newfit = fobj(newpop);
        %适应度判断
        if newfit < fitness(i)
            fitness(i) = newfit;
            Positions(i,:) = newpop;
        end
    end
    
    %随机游走种群更新（局部开发）
    for i = 1:size(Positions,1)
       
        if rand > Pa
        %随机生成解索引
        j = randi([1,SearchAgents_no],1); 
        k = randi([1,SearchAgents_no],1); 
        %偏好随机游走
        newpop = Positions(i,:)+rand.*(Positions(j,:)-Positions(k,:)); %等式6
        
        %超限判断
        Flag4ub = newpop>ub;
        Flag4lb = newpop<lb;
        newpop = (newpop.*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
        
        %计算新个体适应度值 
        newfit = fobj(newpop);
        %适应度判断
        if newfit < fitness(i)
            fitness(i) = newfit;
            Positions(i,:) = newpop;
        end
        end
        
        %全局最优解更新
        if fitness(i) < Leader_score 
            Leader_score = fitness(i);
            Leader_pos = Positions(i,:);
        end
    end
    
    Convergence_curve(t) = Leader_score;  %记录当前迭代的最优适应度值

end
end

%% 种群初始化函数
function Positions = initialization(SearchAgents_no,dim,ub,lb)

Boundary_no= size(ub,2); %限界个数

if Boundary_no == 1 %所有参数的边界相等
    Positions = rand(SearchAgents_no,dim).*(ub-lb)+lb;
end

if Boundary_no>1    %参数边界不等
    for i =1:dim
        ub_i=ub(i);
        lb_i=lb(i);
        Positions(:,i)=rand(SearchAgents_no,1).*(ub_i-lb_i)+lb_i;
    end
end
end

%% 执行 Levy飞行更新蜻蜓位置
function step = Levy(d)

beta = 3/2; 
%等式3
sigma =(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);   %gamma(X)返回在 X 的元素处计算的 gamma 函数值。
u = randn(1,d)*sigma;
v = randn(1,d);
%等式2
step = u./abs(v).^(1/beta);

end

四、输出结果

求解得到x1,x2是：3.1416 2.275
最优解：0.39789

各位阅读文章的朋友觉得文章可以给个好评或者点赞，大家觉得有问题可以在评论区指出来，感谢大家的阅读，祝大家sci多多。