智能优化算法学习第二期：蝙蝠算法（BA）【附MATLAB代码】

一、引言

本期我们介绍蝙蝠算法(Bat Algorithm,BA)。该算法由Yang等人于2010年提出，主要模拟了蝙蝠利用回声定位捕食的行为。因其模型简单、算法参数少、通用性强，目前已在各个领域得到广泛应用。

二、算法原理

2.1 工作原理
①、声波频率和速度调整：每只蝙蝠代表一个解决方案，它们通过调整飞行速度和声波的频率来探索解空间。
②、随机飞行和位置更新：根据当前解决方案、速度（解的变化速度）和最好的解决方案来生成新的解决方案。
③、动态响应和避免障碍：根据当前解的质量来调整搜索范围和速度，优化搜索效率。
④、局部搜索和变异：模拟蝙蝠捕食时的随机和精确的行为，发现潜在的良好解决方案时会进行局部搜索，这可能涉及在当前最优解周围进行随机游走以探索更好的解。

2.2 计算过程
①蝙蝠个体的位置更新公式如下：

其中，β 为[0,1]内一随机值，x*为当前最优个体位置；fi​ 为蝙蝠发出的声波频率，其值介于[fmin,fmax]间，原文设定fmin、fmax 分别为0、100，不过其值也要根据问题域的大小来调整。
在局部搜索阶段, 算法将生成一个新的局部解：

式中ε为[-1,1]内一随机值，At 为第t 次迭代时种群内所有个体的平均响度。
②蝙蝠发出的脉冲响度A 和脉冲发射频率r 会随着迭代进行更新，脉冲响度逐渐减低的同时会提高脉冲发射速率，更新公式如下：

不难看出，当t 趋于无穷时，A、r 均趋于0，这也就表示蝙蝠捕获了猎物而不再发出脉冲。原文设定α=γ=0.9：

2.3 具体步骤
①初始化种群：初始化蝙蝠的位置和速度。
②计算蝙蝠的适应度值：根据蝙蝠的位置计算其适应度值。
③随机化搜索：根据当前的最优解和蝙蝠的位置，随机生成新的位置和速度。
④跟新位置和速度：通过当前位置和速度，更新蝙蝠的位置和速度。
⑤局部搜索：对于一部分蝙蝠，进行局部搜索，以加速收敛速度。
⑥更新最优解：更新当前的最优解。
⑦终止条件：当达到一定的迭代次数或者最优解已经收敛时，停止搜索。

三、算法代码

主程序 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] = BA(pop_num,Max_iter,lb,ub,dim,fobj); %寻找最小值
t2=clock;
time_BA=(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

BA.m

% 蝙蝠优化算法
% 定义函数
% Leader_score：最优适应度值 Leader_pos：最优位置 Convergence_curve：最优适应度变化曲线
% SearchAgents_no：种群个数 Max_iter：进化次数 lb：变量下限 ub：变量上限 dim：待优化参数个数 fobj：适应度函数

function [Leader_score, Leader_pos,Convergence_curve] = BA(SearchAgents_no,Max_iter,lb,ub,dim,fobj)

%参数初始化
A = 0.6.*ones(SearchAgents_no,1);    %初始响度 (不变或者减小)
r = 0.5.*ones(SearchAgents_no,1);        %初始脉冲率 (不变或增加))
r0 = 0.7;                            %r0为初始脉冲发射率 等式19参数
Alpha = 0.9;                         %调节系数  等式18参数
gamma = 0.9;                         %常数      等式19参数
Qmin = 0;                            %最小频率
Qmax = 1;                            %最大频率
 
%种群初始化 initialization：初始化函数
Positions = initialization(SearchAgents_no,dim,ub,lb);
%初始速度与频率
popv = zeros(SearchAgents_no,dim);   % 速度
Q = zeros(SearchAgents_no,1);   % 频率

%计算适应度
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);

%主函数
t = 1; 
while t <= Max_iter
    %种群位置更新
    for i = 1:SearchAgents_no
        Q(i) = Qmin+(Qmax-Qmin)*rand();                         %等式2-2
        popv(i,:) = popv(i,:)+(Positions(i,:)-Leader_pos)*Q(i); %等式2-3
        Stemp = Positions(i,:)+popv(i,:);                       %等式2-4
        
        %局部搜索
        if rand > r(i)
            Amean = mean(A);
            Stemp = Leader_pos+(2*rand(1,dim)-1).*Amean;
        end
        
        %超限判断
        Flag4ub = Stemp>ub;
        Flag4lb = Stemp<lb;
        Stemp = (Stemp.*(~(Flag4ub+Flag4lb)))+ ub.*Flag4ub+lb.*Flag4lb;
        %计算新个体适应度
        fitTemp = fobj(Stemp);
        
        if (fitTemp <= fitness(i))&&(rand<A(i))
            Positions(i,:) = Stemp;
            fitness(i) = fitTemp;
            A(i) = Alpha*A(i);                       %等式2-6
            r(i) = r0*(1-exp(-gamma*t));             %等式2-7
        end
        %全局最优解更新
        if fitTemp <= Leader_score
            Leader_score = fitTemp;
            Leader_pos = Stemp;
        end
    end
    
    Convergence_curve(t) = Leader_score;
    t = t +1;
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

end

四、输出结果

求解得到x1,x2是：-3.1426 12.2815
最优解：0.39791

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