首先在前面的博客中有介绍到蝙蝠算法BA，这是一个在2010年新提出的算法，研究者也对这个算法进行了很多探索，BA算法在一些问题上的效果还是不错的。说明BA算法有它特定的使用范围。
而p-中位问题一个优化问题，np难。
%BA算法解决CPMP问题
%利用BA解决这个问题，主要还是如何设计编码的问题
clc;
clear;
close all;
% Problem Definition
model=CreateModel();
CostFunction=@(ba)MyCost(ba,model);
nba=20; %蝙蝠个体数
N_gen=100; %迭代次数
A=1+rand(nba,1); % 响度，按照p8要求产生[1,2]的随机数
r=rand(nba,1); % 脉冲率,设置为[0,1]的随机数
al = 0.85;
rr = 0.9;
r0 = r;
%f=zeros(1,nba); % 频率,离散蝙蝠算法里面没有频率设计
%生成初始的蝙蝠个体
Fitness = []; %适应度函数
for i = 1:nba
ba(i) = newBA(model);
Fitness(i) =CostFunction(ba(i)) ; %适应度函数
end
[fmin,I]=min(Fitness); % Best BA
BestSol_Cost = fmin ;
BestSol_ba = ba(I);
%总循环
for t = 1:N_gen
for i = 1:nba
%计算速度 用当前解和最优解的汉明距离作为上限，来随机生成速度
hanming = get_hanming( ba(i), BestSol_ba,model );
V.x = unidrnd(hanming.x);
V.y = unidrnd(hanming.y);
%实现x = x + v 进行蝙蝠进化，这里采用在x中随机取点做V.x次交换，在y中随机取点做V.y次交换
x_plus_v(ba(i),V,model);
Fitness(i) =CostFunction(ba(i)) ; %更新适应度
Fnew = inf;
%以一定概率产生一个解
if rand>r(i,1)
baNew =newBA(model);%这里的新的蝙蝠个体是由当前全局最好的个体产生的
Fnew=CostFunction(baNew);
end
if ((Fnew<Fitness(i)) && (rand<A(i,1)))
ba(i)=baNew;
Fitness(i)=Fnew;
A(i,1) = al*A(i,1); %对响度进行更新
r(i,1) = r0(i,1)*(1-exp(-1*rr*t )); %对脉冲率进行更新
end
% 更新当前最优解
if Fitness(i)<=BestSol_Cost
BestSol_ba = ba(i);
BestSol_Cost = Fitness(i);
end
end
disp(['Iteration ', num2str(t),'BestSol_Cost=',num2str(BestSol_Cost)]);
end
function baNew =newBA(model)
baNew.x = randperm(model.n * model.m); %随机化为一个n*p长的一维序列
baNew.visit_x = zeros(model.n,model.m);
baNew.real_x = zeros(model.n,model.m);
baNew.y = randperm(model.m); %随机化为一个m长的一维序列
baNew.visit_y = zeros(1,model.m);
baNew.real_y = zeros(1,model.m);
baNew = FFix(baNew,model); %把蝙蝠体解释成p-中位问题的一个可行解,即把x->real_x , y->real_y
end
%借鉴遗传算法对运输问题的编码方案
function ba = FFix(ba,model)
x = ba.x;
y = ba.y;
m = model.m;
n = model.n;
p = model.p;
s = model.s;
d = model.d;
%先把y解释为"m选p的设施选择"
py = p/n; %每个设施被选中的概率，初始化时每个设施被选中的概率是均等的
yf = 1;
seed_y = 1;
while(yf<=p) %m个设施中选出来p个，存放到ba.real_y中
if(rand<py)
poi_y = mod( seed_y-1, m )+1;
ba.real_y( y(poi_y) ) = 1;
yf = yf + 1;
end
seed_y = seed_y + 1;
end
%在设施已经选好的基础上，进一步把x解释为"客户的选择"
for c = 1:n*m
i = floor ((x(c)-1)/m )+1 ; %通过x(c)拿到行号
j = mod( x(c)-1, m)+1; %通过x(c)拿到列号
if(s(j)>=d(i)||ba.real_y(j)==1 ) %通过ba.real_y(j)==1 筛出上面一步选出的p个设施
ba.real_x(i,j) = 1;
s(j) = s(j) - d(i);
end
end
end
function hanming = get_hanming( ba, BestSol_ba,model )
%首先计算ba.real_x的汉明距离
hanming.x = 0;
hanming.y = 0;
for i = 1:model.n
for j = 1:model.m
hanming.x = hanming.x + abs( ba.real_x(i,j)-BestSol_ba.real_x(i,j) );
end
end
%计算ba.real_y的汉明距离
for j = 1:model.m
hanming.y = hanming.y + abs( ba.real_y(j)-BestSol_ba.real_y(j) );
end
end
%实现x = x + v 进行蝙蝠进化，这里采用在x中随机取点做V.x次交换，在y中随机取点做V.y次交换
function x_plus_v(ba,V,model)
time_x = V.x;
time_y = V.y;
for i = 1:time_x
poi = unidrnd(model.n * model.m);
if(poi==model.n*model.m)
poi_next = 1;
else
poi_next = poi + 1;
end
num_temp = ba.x(poi_next);
ba.x(poi_next) = ba.x(poi);
ba.x(poi) = num_temp;
end
for k = 1:time_y
poi_y = unidrnd(model.m);
if(poi_y==model.m)
poi_y_next = 1;
else
poi_y_next = poi_y + 1;
end
num_temp_y = ba.y(poi_y_next) ;
ba.y(poi_y_next) = ba.y(poi_y) ;
ba.y(poi_y) = num_temp_y;
end
end
参考文献：蝙蝠算法的改进与应用 何子旷 广东工业大学硕士学位论文 2016.5