智能优化算法——蝙蝠算法（PythonMatlab实现）
⽬录
1 前⾔
蝙蝠算法是2010年杨教授基于群体智能提出的启发式搜索算法，是⼀种搜索全局最优解的有效⽅法。该算法基于迭代优化，初始化为⼀组随机解，然后迭代搜寻最优解，且在最优解周围通过随机飞⾏产⽣局部新解，加强局部搜索速度。该算法具有实现简单、参数少等特点。
该算法主要⽤于⽬标函数寻优，基于蝙蝠种群利⽤产⽣的声波搜索猎物和控制飞⾏⽅向的特征来实现函数的寻优。以⼀只蝙蝠作为基本单元，且每只蝙蝠都有⼀个适应值来对函数解空间进⾏优化。每只蝙蝠可以调
整⾃⾝发射声波的响度、频率等对空间进⾏搜索，使整个种群的活动逐步由⽆序变为有序。但蝙蝠算法在寻优末段容易陷⼊局部的极值，本⽂引⼊变速度权重因⼦修正系数(和粒⼦群参数调整的⾃适应过程相似)，尽可能避免局部极值的困境，从⽽达到全局寻优的效果。
2 蝙蝠算法原理细讲
⾸先在⾃变量范围内产⽣蝙蝠的随机位置；
然后每个蝙蝠向周围某处移动，飞⾏靠声波反馈来更变⽅向，⼜向周围移动是随机的，蝙蝠下⼀时刻移动到某⼀位置有⼀定的出现频率，所以本⽂在运动公式中加了声波和频率两个因素。每⼀时刻的位移看作是⼀次飞⾏，每次飞⾏的距离很短（距离长短反映出搜素精度）。
每只蝙蝠初始位置是随机的，在初始位置时刻其中⼀只蝙蝠对应的函数最值，算法开始会记录该位置，然后所有蝙蝠逐渐向该位置靠近，飞⾏⽅向⼤致向当前最值⽅向（即该⽅向的位置该蝙蝠的出现频率更⾼），但是在飞⾏⽅向也是随机飞⾏的，相当于逐步搜索过去。蝙蝠飞⾏过程中不会因为当前飞⾏出现的最⼤值位置⽽改变种群的飞⾏趋势，这样可以尽量避免算法陷⼊局部极值。
算法每飞⾏⼀次就记录⼀次种群所处位置的最值，最后找出记录中的最值。
算法有两个参数可以影响最终的结果：种群数量和飞⾏次数。其中种群数量的影响是最⼤的。
3 详细步骤
3.1 初始化相关参数
蝙蝠的位置为Xi，飞⾏速度Vi，声⾳响度为Ai，频率yi范围，设有⽬标函数为 ：
w(t)其中为时刻变速惯性权重因⼦，作⽤是使蝙蝠的前期搜索对后期搜索提供参照，wmax为w(t)的最⼤值、wmin为w(t)的最⼩值；
，⼀般取2，Tmax为最⼤迭代次数；为当前位置最优解；y(i)为频率满⾜正态均匀分布的⼀个随机数，β是⼀个随机变量，且。开始运⾏时，蝙蝠在随机进⾏频率分配。
3.4 通过蝙蝠多次飞⾏产⽣多个新解，进⾏全局搜索，若得到的新解
3.5 排列所有蝙蝠的位置，并找出当前最优值及对应的位置
4 Python实现
4.1 代码
#=========导⼊相关库===============
import numpy as np
from numpy.random import random as rand
#========参数设置==============
# objfun:⽬标函数
# N_pop: 种群规模，通常为10到40
# N_gen: 迭代数
# A: 响度（恒定或降低）
# r: 脉冲率（恒定或减⼩）
# 此频率范围决定范围
# 如有必要，应更改这些值
# Qmin: 频率最⼩值
# Qmax: 频率最⼤值
# d: 维度
# lower: 下界
# upper: 上界
def bat_algorithm(objfun, N_pop=20, N_gen=1000, A=0.5, r=0.5,
Qmin=0, Qmax=2, d=10, lower=-2, upper=2):
N_iter = 0 # Total number of function evaluations
#=====速度上下限================
Lower_bound = lower * np.ones((1,d))
Upper_bound = upper * np.ones((1,d))
Q = np.zeros((N_pop, 1)) # 频率
v = np.zeros((N_pop, d)) # 速度
S = np.zeros((N_pop, d))
#=====初始化种群、初始解=======
# Sol = np.random.uniform(Lower_bound, Upper_bound, (N_pop, d)) # Fitness = objfun(Sol)
Sol = np.zeros((N_pop, d))
Fitness = np.zeros((N_pop, 1))
for i in range(N_pop):
Sol[i] = np.random.uniform(Lower_bound, Upper_bound, (1, d)) Fitness[i] = objfun(Sol[i])
#====找出初始最优解===========
fmin = min(Fitness)
Index = list(Fitness).index(fmin)
best = Sol[Index]
#======开始迭代=======
for t in range(N_gen):
#====对所有蝙蝠/解决⽅案进⾏循环 ======
for i in range(N_pop):
# Q[i] = Qmin + (Qmin - Qmax) * np.random.rand
Q[i] = np.random.uniform(Qmin, Qmax)
v[i] = v[i] + (Sol[i] - best) * Q[i]
S[i] = Sol[i] + v[i]
#===应⽤简单的界限/限制====
Sol[i] = simplebounds(Sol[i], Lower_bound, Upper_bound)
# Pul rate
if rand() > r:
# The factor 0.001 limits the step sizes of random walks
S[i] = best + 0.001*np.random.randn(1, d)
#====评估新的解决⽅案 ===========
# print(i)
Fnew = objfun(S[i])
#====如果解决⽅案有所改进，或者声⾳不太⼤，请更新====
if (Fnew <= Fitness[i]) and (rand() < A):
Sol[i] = S[i]
Fitness[i] = Fnew
Fitness[i] = Fnew
#====更新当前的最佳解决⽅案======
if Fnew <= fmin:
best = S[i]
fmin = Fnew
N_iter = N_iter + N_pop
print('Number of evaluations: ', N_iter)
print("Best = ", best, '\n fmin = ', fmin)
return best
def simplebounds(s, Lower_bound, Upper_bound):
Index = s > Lower_bound
s = Index * s + ~Index * Lower_bound
Index = s < Upper_bound
s = Index * s + ~Index * Upper_bound
return s
#====⽬标函数=============
def test_function(u):
a = u ** 2
return a.sum(axis=0)
if __name__ == '__main__':
# print(bat_algorithm(test_function))
bat_algorithm(test_function)
4.2 结果
5 Matlab实现
5.1 代码
clear
wmax=0.9;%惯性权重最⼤值
wmin=0.4;%惯性权重最⼩值
n=10000; % 群体⼤⼩
A=rand(1,n); % 声⾳响度 (不变或者减⼩)
%% 频率范围
Qmin=0; % 最低频率
Qmax=1; % 最⾼频率
d=2;% 搜索变量的维数（即频率和速度）
%% 初始矩阵
Q=zeros(n,1); % 频率矩阵初始化
v=zeros(n,d); % 速度矩阵初始化，初始化意义就是产⽣⼀个初始矩阵
%% x⾃变量范围
u=-3;
o=12.1;
% y⾃变量范围
p=4.1;
l=5.8;
%% 初始化群体/解
for i=1:n
Sol(i,1)=-3+(12.1+3)*rand(1,1);%x⾃变量范围【-3，12.1】
Sol(i,2)=4.1+(5.8-4.1)*rand(1,1);%y⾃变量【4.1，5.8范围】
%将随机⽣成的两个⾃变量带⼊函数式
Fitness(i)=Fun(Sol(i,:));%函数值
end
%% 寻找当前最优解
[fmax,I]=max(Fitness);
best=Sol(I,:);
T=100;%飞⾏次数
%% 开始飞⾏
for t=1:T
for i=1:n,
Q(i)=Qmin+(Qmin-Qmax)*rand;%rand均匀分布的随机数
%v(i,:)=v(i,:)+(Sol(i,:)-best)*Q(i);（原速度）
w=(wmax-wmin)*exp(-2*(t/T)^2)+wmin;%惯性权重因⼦
v(i,:)=w*v(i,:)+(Sol(i,:)-best)*A(i)*Q(i);%更改后的速度
S(i,:)=Sol(i,:)+v(i,:);%位置移动
%% 边界问题，如果下次飞⾏超出⾃变量范围外了，那么下次飞⾏的位置为投影在的边界上的位置 %x轴
if S(i,1)>o
S(i,1)=o;
end
if S(i,1)<u
S(i,1)=u;
end
%y轴
if S(i,2)>l
S(i,2)=l;
end
if S(i,2)<p
S(i,2)=p;
end
%% 评估该次飞⾏后产⽣的新解
Fnew(i)=Fun(S(i,:));
end
[Fmax,Z]=max(Fnew);%找出该次飞⾏后产⽣的最⼤值
C(t,:)=S(Z,:);
FFnew(t)=Fmax;
end
[Ffmax,N]=max(FFnew);%找出整个飞⾏过程中的最⼤值
M=C(N,:)
Ffmax
