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