蝙蝠算法

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 频率更新公式

\[f_i = f_{min} + (f_{max} - f_{min})\beta \]

其中，\(\beta \in [0,1]\)是一个随机向量；\(x_*\) 是当前最优解。
3.2 速度更新公式

\[v_i^{t} = v_i^{t-1} + (x_{i}^{t} - x_*)f_i \]

3.3 位置更新公式

\[x_i^{t} = x_i^{t-1} + v_i^{t} \]

3.4 局部搜索更新公式

\[x_i^{new} = x_i^{old} + \epsilon A^{t} \\ \]

其中，\(\epsilon \in [-1,1]\)是一个随机数，\(A^{t}\)是当前迭代的平均响度。
3.5 响度和脉冲发生率更新公式
脉冲发射的响度\(A_i\)和脉冲发生率\(r_i\)要随着迭代过程的进行来更新。蝙蝠一旦发现了猎物，响度会逐渐降低，同时脉冲速率就会提高，响度会以任意简便值改变。

\[A_{i}^{t+1} = \alpha A_{i}^{t} \\ r_{i}^{t+1} = r_{i}^{0}[1 - exp(-\gamma t)] \]

其中，\(\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)