BA算法解决p-中位问题

首先在前面的博客中有介绍到蝙蝠算法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