蝙蝠算法
2012年，英国剑桥大学学者杨新社提出一种新的元启发式优化算法-蝙蝠算法(Bat Algorithm, BA)，该算法通过模拟蝙蝠回声定位行为来寻找函数优化问题的最优解。
1. 蝙蝠算法的基本思想
由于蝙蝠的回声定位行为与函数优化相似，所以可以利用蝙蝠的回声定位行为来寻找最优解。蝙蝠算法把蝙蝠看作优化问题的可行解，通过模拟复杂环境中精确捕获食物的机制解决优化问题。
首先，在搜索空间随机分布若干个蝙蝠，确定种群个体的初始位置及初始速度，对种群中各个蝙蝠进行适应度评价，寻找最优个体位置；
然后，通过调整频率产生新的解并修改个体的飞行速度和位置。在蝙蝠的速度和位置的更新过程中，频率本质上控制着这些蝙蝠群的移动步伐和范围；
蝙蝠在寻优过程中，通过调节脉冲发生率和响度促使蝙蝠朝着最优解方向移动。蝙蝠在刚开始搜索时具有较小的脉冲发生率，蝙蝠有较大的概率在当前最优解周围进行局部搜索，同时较大的响度使得局部搜索范围比较大，有较大的概率探测到更好的解。随着迭代的增加，脉冲发生率增加，响度减少，局部搜索概率减少，局部挖掘的范围也很小，蝙蝠不断扫描定位目标，最终搜索到最优解。
2. 蝙蝠算法的数学描述
为了能使蝙蝠回声定位机制形成算法，有必要对蝙蝠回声定位及飞行速度、位置进行理想化建模：
(1) 所有的蝙蝠利用超声波回声的感觉差异判断猎物、障碍物之间的差异；
(2) 蝙蝠是以速度\(v_i\)、位置\(x_i\)、固定频率\(f_{min}\)、可变化波长\(\lambda\)和响度\(A_0\)随机飞行的，并用不同的波长\(\lambda\)(或者频率\(f\))和响度\(A_0\)搜索猎物。它们会根据接近猎物的程度自动调整它们发出脉冲的波长(或频率)；
(3) 尽管响度会以更多的方式变化，可以假定它的变化是从一个很大的值\(A_0\)到最小值\(A_{min}\)；
(4) 由于计算量的问题，不能使用无限追踪来估计时间延迟和三维地形；
(5) 为了简单起见，使用一些近似值。一般设置频率范围为\([f_{min}, f_{max}]\)对应的波长范围为\([\lambda_{min}, \lambda_{max}]\)。
对于给定的问题，为了便于实现，可以使用任意波长，并且可以通过调整波长来搜索范围，而可探测的区域的选择方式为先选择感兴趣的区域，然后慢慢缩小。因为波长\(\lambda \times f\)为常数，所以可以在固定波长\(\lambda\)时，改变频率。
为了简单期间，可以假定\(f \in [0, f_{min}]\)。显然，较高的频率有较短的波长和较短的搜索距离。通常蝙蝠的搜索范围在几米内。脉冲发生率可以设定为\([0,1]\)范围内，其中0表示没有发出脉冲，1表示脉冲发生率最大。
3. 蝙蝠的速度和位置更新公式
3.1 频率更新公式
其中，\(\beta \in [0,1]\)是一个随机向量；\(x_*\) 是当前最优解。
3.2 速度更新公式
3.3 位置更新公式
3.4 局部搜索更新公式
其中，\(\epsilon \in [-1,1]\)是一个随机数，\(A^{t}\)是当前迭代的平均响度。
3.5 响度和脉冲发生率更新公式
脉冲发射的响度\(A_i\)和脉冲发生率\(r_i\)要随着迭代过程的进行来更新。蝙蝠一旦发现了猎物，响度会逐渐降低，同时脉冲速率就会提高，响度会以任意简便值改变。
其中，\(\alpha\)和\(\gamma\)是常量。参数的选择需要一定的经验。初始时，每只蝙蝠所发出的响度和脉冲发生率的值都是不同的，这可以通过随机选择。初始的响度\(A_{i}^{0}\)通常在\([0,1]\)之间，而初始脉冲发生率一般去在0附近。
4. 算法实现
版本一
% =========================================================================
% Title: Bat Algorihtm
% Author: Lee WenTsao
% Time: 2022-04-12
% E-mail: liwenchao36@163.com
% =========================================================================
clc;
clear;
%% 问题参数定义
n = 20; % 种群规模
d = 10; % 维度
fobj = @ sphere; % 目标函数
%% 蝙蝠算法参数定义
t_max = 1000; % 最大迭代次数
A = 1; % 初始响度
r0 = 1; % 初始脉冲发生率
alpha = 0.97;
gamma = 0.1;
Freq_min = 0; % 蝙蝠发射频率的下界
Freq_max = 2; % 蝙蝠发射频率的上界
iter = 0;
%% 初始化动态参数
Freq = zeros(n,1);
v = zeros(n,d); % 蝙蝠的初始速度
lb = -5*ones(1,d); % 下界
ub = 5*ones(1,d); % 上界
%% 初始化种群
for i=1:n
Sol(i, :) = lb + (ub - lb).*rand(1, d); % 随机初始化种群
Fitness(i) = fobj(Sol(i, :)); % 评估
end
[fmin, idx] = min(Fitness);
best = Sol(idx,:);
%% 添加云图
%% 开始仿真
while iter<t_max
r = r0*(1-exp(-gamma*iter));
A = alpha*A;
for i=1:n
% 全局搜索
Freq(i) = Freq_min + (Freq_max - Freq_min)*rand; % 频率更新公式
v(i,:)=v(i,:)+(Sol(i,:)-best)*Freq(i); % 速度更新
S(i, :) = Sol(i, :) + v(i,:); % 位置更新
% 局部搜索
if rand<r
S(i, :) = best + 0.1*randn(1,d)*A;
end
% 边界预处理
S(i, :) = simplebounds(S(i, :), lb, ub);
% 评估
Fnew = fobj(S(i, :));
if Fnew<=Fitness(i) & rand>A
Sol(i,:) = S(i, :);
Fitness(i) = Fnew;
end
if Fnew<fmin
best = S(i,:);
fmin = Fnew;
end
end
iter = iter + 1;
if ~mod(iter,10)
disp(['Iteration=',num2str(iter),' fmin=',num2str(fmin)]);
end
end
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;
end
第二版本(杨老师原版)
% ----------------------------------------------------------------------- %
% The bat algorithm (BA) for continuous function optimization (demo) %
% Programmed by Xin-She Yang @Cambridge University 2010 %
% ----------------------------------------------------------------------- %
% References: ----------------------------------------------------------- %
% 1) Yang X.S. (2010). A New Metaheuristic Bat-Inspired Algorithm, In: %
% Nature-Inspired Cooperative Strategies for Optimization (NICSO 2010), %
% Studies in Computational Intelligence, vol 284. Springer, Berlin. %
% https://doi.org/10.1007/978-3-642-12538-6_6 %
% 2) Yang X.S. (2014). Nature-Inspired Optimization Algorithms, Elsevier. %
% ----------------------------------------------------------------------- %
function [best,fmin,N_iter]=bat_algorithm_new(inp)
if nargin<1,
inp=[20 1000]; % Default values for n=10 and t_max=1000
end
% Initialize all the default parameters
n=inp(1); % Population size, typically 20 to 40
t_max=inp(2); % Maximum number of iterations
A=1; % Initial loudness (constant or decreasing)
r0=1; % The initial pulse rate (constant or decreasing)
alpha=0.97; % Parameter alpha
gamma=0.1; % Parameter gamma
% Frequency range
Freq_min=0; % Frequency minimum
Freq_max=2; % Frequency maximum
t=0; % Initialize iteration counter
% Dimensions of the search variables
d=10;
% Initialization of all the arrays
Freq=zeros(n,1); % Frequency-tuning array
v=zeros(n,d); % Equivalnet velocities or increments
Lb=-5*ones(1,d); % Lower bounds
Ub=5*ones(1,d); % Upper bounds
% Initialize the population/solutions
for i=1:n,
Sol(i,:)=Lb+(Ub-Lb).*rand(1,d);
Fitness(i)=Fun(Sol(i,:));
end
% Find the best solution of the initial population
[fmin,I]=min(Fitness);
best=Sol(I,:);
% Start the iterations -- the Bat Algorithm (BA) -- main loop
while (t<t_max)
% Varying loundness (A) and pulse emission rate (r)
r=r0*(1-exp(-gamma*t));
A=alpha*A;
% Loop over all bats/solutions
for i=1:n,
Freq(i)=Freq_min+(Freq_max-Freq_min)*rand;
v(i,:)=v(i,:)+(Sol(i,:)-best)*Freq(i);
S(i,:)=Sol(i,:)+v(i,:);
% Check a switching condition
if rand<r,
S(i,:)=best+0.1*randn(1,d)*A;
end
% Check if the new solution is within the simple bounds
S(i,:)=simplebounds(S(i,:),Lb,Ub);
% Evaluate new solutions
Fnew=Fun(S(i,:));
% If the solution improves or not too loudness
if ((Fnew<=Fitness(i)) & (rand>A)),
Sol(i,:)=S(i,:);
Fitness(i)=Fnew;
end
% Update the current best solution
if Fnew<=fmin,
best=S(i,:);
fmin=Fnew;
end
end % end of for i
t=t+1; % Update iteration counter
% Display the results every 100 iterations
if ~mod(t,100),
disp(strcat('Iteration = ',num2str(t)));
best, fmin
end
end % End of the main loop
% Output the best solution
disp(['Best =',num2str(best),' fmin=',num2str(fmin)]);
% Application of simple bounds/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;
% The cost function or objective function
function z=Fun(x)
z=sum((x-2).^2); % Optimal solution fmin=0 at (2,2,...,2)
