ABC算法求解TSP问题

clear all;
close all;
clc
C=[1304,2312;3639,1315;4177,2244;3712,1399;3488,1535;
     3326,1556;3238,1229;4196,1044;4312,790;  4386,570  ;
     3007,1970;2562,1756;2788,1491;2381,1676;1332,695  ;
     3715,1678;3918,2179;4061,2370;3780,2212;3676,2578;
     4029,2838;4263,2931;3429,1908;3507,2376;3394,2643;
     3439,3201;2935,3240;3140,3550;2545,2357;2778,2826;
     2370,2975];     %城市坐标
D=size(C,1);         %种群规模
Dis=zeros(D,D);      %城市间距离
for i=1:D
    for j=1:D
        Dis(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;
    end
end
N=50;                %蜜源数量
swarm=zeros(N,D); 
for i=1:N       
    swarm(i,:)=randperm(D);      
    fit(i)=Dis(swarm(i,D),swarm(i,1));
    for j=1:D-1   %求累积距离
        fit(i)=fit(i)+Dis(swarm(i,j),swarm(i,j+1));
    end
end
X_best=zeros(1,D);  %全局最优解
fit_best=inf;       %全局最优值
for i=1:N
    if fit(i)<fit_best      %更新最优值
        fit_best=fit(i);
        X_best=swarm(i,:);
    end
end
limit= round(0.6*D*N);    %蜜源不能被更新的最大限制次数
t=0;
trial=zeros(N,1);         %记录限制次数，达到后就舍弃该蜜源，选择新的
pro=zeros(N,1);
trace=[];                 %记录每一代的最优值
Gmax=500;                 %最大迭代次数
%% 雇佣蜂阶段
while t<Gmax
    for i=1:N
        K = setdiff(1:N,i);  %选择除了i的个体
        k = K(randi(N-1));   %在k中随机选择一个个体
        V=swarm(i,:);
        set=[];
        temp_swarm1=swarm(i,:);
        temp_swarm2=swarm(k,:);
        for j=1:D    %邻域搜索选择蜜源
            index=find(temp_swarm1(j)==temp_swarm2);
            if index~=j
                set=[set;[j,index]];
                temp=temp_swarm1(j);
                temp_swarm1(j)=temp_swarm1(index);
                temp_swarm1(index)=temp;
            end
        end
  if isempty(set)
            set=randi(D,[1,2]);
        end
        j=randi(size(set,1));
        temp=V(set(j,1));
        V(set(j,1))=V(set(j,2));
        V(set(j,2))=temp;
        Vfit=Dis(V(D),V(1));
        for j=1:D-1
            Vfit=Vfit+Dis(V(j),V(j+1));
        end
 if Vfit<fit(i)  %贪婪算法更新蜜源位置，更新了停滞次数不变
            swarm(i,:)=V;
            fit(i)=Vfit;
            trial(i)=0;
        else
            trial(i)=trial(i)+1;%记录停滞次数
        end
    end
Fit=1-(fit-min(fit))/(max(fit)-min(fit)+0.01);  %标准化（最小化问题） 
    for i=1:N
        pro(i)=Fit(i)/sum(Fit);  %计算去往每个蜜源的概率
    end
    cum_pro=cumsum(pro);    %计算累积概率
    i=1;
    n=1;
%%观察蜂
 while n<=N
        if rand<cum_pro(i)
            %选择除去i的其他个体
            K = setdiff(1:N,i); 
            k = K(randi(N-1));set=[];%从k中随机选择一个个体
            temp_swarm1=swarm(i,:);
            temp_swarm2=swarm(k,:);
            for j=1:D
                index=find(temp_swarm1(j)==temp_swarm2);
                if index~=j
                    set=[set;[j,index]];
                    temp=temp_swarm1(j);
                    temp_swarm1(j)=temp_swarm1(index);
                    temp_swarm1(index)=temp;
                end
            end          
   if isempty(set)
                set=randi(D,[1,2]);
            end
            j=randi(size(set,1));
            temp=V(set(j,1));
            V(set(j,1))=V(set(j,2));
            V(set(j,2))=temp;
            Vfit=Dis(V(D),V(1));
            for j=1:D-1
                Vfit=Vfit+Dis(V(j),V(j+1));%更新蜜源位置
            end
    if Vfit<fit(i)     %贪婪算法更新蜜源位置，更新了停滞次数不变
                swarm(i,:)=V;
                fit(i)=Vfit;
                trial(i)=0;
            else
                trial(i)=trial(i)+1;   %记录此次蜜源的更新停滞次数
            end
            n=n+1;
        end
        i=i+1;
        if i==N
            i=1;
        end
 end
 %%侦查蜂
    for i=1:N
          if trial(i)>limit    %超过限制次数，则放弃该位置，重新搜索
            swarm(i,:)=randperm(D);
            fit(i)=Dis(swarm(i,D),swarm(i,1));
              for j=1:D-1
                  fit(i)=fit(i)+Dis(swarm(i,j),swarm(i,j+1));
              end
            trial(i)=0;
            end
       end
 for i=1:N
        if fit(i)<fit_best   %更新全局最优值
            fit_best=fit(i);
            X_best=swarm(i,:);
        end
    end
    t=t+1;
    trace=[trace;fit_best]  %记录每一代的最优值
end
fit_best
figure
for i=1:D-1
    plot([C(X_best(i),1),C(X_best(i+1),1)],[C(X_best(i),2),C(X_best(i+1),2)],'bo-','linewidth',1.5)
    text(C(X_best(i),1)+50,C(X_best(i),2)+50,[num2str(X_best(i))])
    hold on;
end
plot([C(X_best(D),1),C(X_best(1),1)],[C(X_best(D),2),C(X_best(1),2)],'co-','linewidth',1.5)
text(C(X_best(D),1)+50,C(X_best(i),2)+50,[num2str(X_best(D))])
title('最短距离');
hold on;
figure(2),
plot(trace)
xlabel('迭代次数')
ylabel('目标函数')
title('目标函数曲线')