含重复点的蚁群算法

背景：以论文《汉中市城市生活垃圾收运路线优化研究》为背景。共37个位置，一个是车库，一个是垃圾处理中心，剩下35个是垃圾收集站。每次都是垃圾搬运车从车库出发，经过7个垃圾收集站，运到到垃圾处理中心，重复5次，直到35个垃圾收集站的垃圾都收集完。

1.先画出各个位置观察。

clear all
close all;
clc;


%%1.input data
C = [4.8 0.2;         %车库位置
    0.9 5.0;     %垃圾处理中心位置
    8.6 7.8;
    9.9 3.9;
    0.6 5.4;
    0.1 2.2;
    9.7 7.6;
    5.1 2.1;
    2.4 1.0;
    3.7 2.7;
    5.7 3.0;
    9.5 5.1;
    5.6 3.2;
    6.9 6.2;
    4.4 6.5;
    2.4 1.1;
    5.6 3.7;
    9.0 1.1;
    5.9 3.8;
    4.3 2.4;
    7.1 5.4;
    6.7 9.7;
    2.7 7.1;
    8.7 7.2;
    3.7 6.5;
    7.3 3.4;
    8.1 8.2;
    0.2 2.3;
    6.7 5.1;
    4.0 3.4;
    9.1 2.3;
    4.7 9.0;
    7.8 2.7;
    7.6 9.2;
    1.2 7.2;
    2.5 1.5;
    8.8 2.8];
figure(1);
subplot(2,2,1)
x1 = C(:,1);
y1 = C(:,2);
scatter(x1,y1,'.')
title('垃圾收集点位置')

得到图。

2.计算两点间距离

%%2.计算两两垃圾收集站的距离
[M,N] = size(C);
distance = zeros(M,M);     % 用来记录任意两点的距离
% for m=1:M
%     for n=1:M
%         distance(m,n) = sqrt(sum((C(m,:)-C(n,:)).^2));
%     end
% end
for i = 1:M
    for j = 1:M
        if i ~= j
            distance(i,j) = sqrt(sum((C(i,:) - C(j,:)).^2));
        else
            distance(i,j) = 1e-4;  %对角矩阵 距离为0，我将对角矩阵赋值一个很小的距离，不影响计算。因为前文图中公式里有距离的倒数   
        end
    end    
end
for m=1:M
    for n=1:M
        x2(37*(m-1)+n,1)=38-m;
        y2(37*(m-1)+n,1)=n;
        distance_norm(37*(m-1)+n,1)=distance(n,m);
    end
end
subplot(2,2,2);
scatter(x2,y2,50,distance_norm,'s','filled');
title('距离矩阵');

把距离矩阵可视化，越黄表明距离越大，越蓝表示距离越小。得到图。

3.进行每次迭代。
基本思想是，每次有m只蚂蚁，在m只蚂蚁里面找到最优的。再进行下次迭代，每次迭代的路径是根据信息素和轮盘赌。车库是第一行数据，垃圾处理中心是第二行数据，我们需要求的最优路径可理解为，从车库（经过7个垃圾收集站）到垃圾处理中心，再从垃圾处理中心（经过7个垃圾收集站）到车库（因为这个距离和从车库经7个垃圾收集站到处理中心的距离一样，转化成这样更有连续性。），以此往返，完成35个垃圾收集站的采集。

进行参数定义。

%%3.初始化参数
p=M;     
m = 65;                              % 蚂蚁数量
alpha = 1;                           % 信息素重要程度因子
beta = 5;                            % 启发函数重要程度因子
rho = 0.5;                           % 信息素挥发因子
Q = 1;                               % 常系数
Eta = 1./distance;                   % 启发函数
Tau = ones(p,p);                     % 信息素矩阵
Table = zeros(m,p+4);                  % 路径记录表
%%Table_all=zeros(m,p+5);              %完整的路线记录表
iter = 1;                            % 迭代次数初值
iter_max = 1000;                      % 最大迭代次数 
Route_best = zeros(iter_max,p);      % 各代最佳路径 
Route_best_all = zeros(iter_max,p+4);      % 各代最佳路径
Length_best = zeros(iter_max,1);     % 各代最佳路径的长度  
Length_ave = zeros(iter_max,1);      % 各代路径的平均长度 
a=ones(m,1);
b=2*ones(m,1);

进行每次迭代。

%%开始迭代找最佳路径
while iter<=iter_max
    
    %随机产生各个蚂蚁的起点城市
    start =zeros(m,1);
    for i=1:m  %遍历每一只蚂蚁
        temp=randperm(1); %产生随机数
        start(i)=temp(1);   %将每次产生的随机数进行赋值，得到每次蚂蚁的起点
    end
    
    Table(:,1)=start;  %路径记录表
    Table(:,9)=b;
    Table(:,17)=a;
    Table(:,25)=b;
    Table(:,33)=a;
    Table(:,41)=b;

    citys_index=1:p;
    
    
    %逐个蚂蚁路径选择
    for i=1:m
        j=2;
        while(j<=p+3)   %从第二个城市开始
            yy=1;
            for w=1:p+3
                %num=ones(50,39);
               if (Table(i,w) ~=0)
                  num( i,yy)=Table(i,w);
                  yy_tabu=yy;
                  yy=yy+1;
               end
            end
            %yy_tabu=yy;
            yy=1;
            %tabu=Table(i,1:(j-1));  %已访问的禁忌表
            tabu=num(i,1:(yy_tabu));  %已访问的禁忌表
            allow_index=~ismember(citys_index,tabu);  %取出未访问的索引
            allow = citys_index(allow_index);  % 待访问的城市集合
            pp=allow;
            
            

                %计算城市间的转移概率
                 for k=1:length(allow)
                     pp(k) = Tau(tabu(end),allow(k))^alpha * Eta(tabu(end),allow(k))^beta;
                 end
                 pp=pp/sum(pp);
            
                %轮盘赌法选择下一个访问城市
                if(j==9||j==17||j==25||j==33||j==41)
                    Table(i,j)=Table(i,j);
                    j=j+1; 
                    else
                        pc=cumsum(pp);
                        target_index=find(pc>=0.5*rand);
                        target=allow(target_index(1));
                        Table(i,j)=target;
                        j=j+1;
                    %end
                end
  
        end
    end
    
    %计算各个蚂蚁的路径距离
    Length=zeros(m,1);
    for i=1:m
        Route=Table(i,:);
        for j=1:(p+3)
           Length(i)= Length(i)+distance(Route(j),Route(j+1));
        end
    end
    
    %计算最短路径距离和平均距离
    if iter == 1
          [min_Length,min_index] = min(Length);
          Length_best(iter) = min_Length;  
          Length_ave(iter) = mean(Length);
          Route_best_all(iter,:) = Table(min_index,:);
      else
          [min_Length,min_index] = min(Length);
          Length_best(iter) = min(Length_best(iter - 1),min_Length);
          Length_ave(iter) = mean(Length);
          if Length_best(iter) == min_Length
              Route_best_all(iter,:) = Table(min_index,:);
          else
              Route_best_all(iter,:) = Route_best_all((iter-1),:);
          end
    end
    
    Delta_Tau = zeros(p,p);
      
      %逐个蚂蚁计算
      for u=1:m
          %逐个城市更新
          for w=1:p+3
               Delta_Tau(Table(u,w),Table(u,w+1)) = Delta_Tau(Table(u,w),Table(u,w+1)) + Q/Length(u);   
          end
      end
       Tau = (1-rho) * Tau + Delta_Tau;
       
     %画图
     subplot(2,2,3);
      for i=1:40
        if(i<=8)
            plot([C( Route_best_all(iter,i),1),C( Route_best_all(iter,i+1),1)],[C( Route_best_all(iter,i),2),C( Route_best_all(iter,i+1),2)],'-*y');
            hold on;
        else if(i<=16)
                 plot([C( Route_best_all(iter,i),1),C( Route_best_all(iter,i+1),1)],[C( Route_best_all(iter,i),2),C( Route_best_all(iter,i+1),2)],'-*m');
                 hold on;
            else if(i<=24)
                     plot([C( Route_best_all(iter,i),1),C( Route_best_all(iter,i+1),1)],[C( Route_best_all(iter,i),2),C( Route_best_all(iter,i+1),2)],'-*b');
                     hold on;
                else if(i<=32)
                         plot([C( Route_best_all(iter,i),1),C( Route_best_all(iter,i+1),1)],[C( Route_best_all(iter,i),2),C( Route_best_all(iter,i+1),2)],'-*k');
                         hold on;
                    else
                         plot([C( Route_best_all(iter,i),1),C( Route_best_all(iter,i+1),1)],[C( Route_best_all(iter,i),2),C( Route_best_all(iter,i+1),2)],'-*g');
                         hold on;
                    end
                end
            end
        end          
    end
    %plot([C( Route_best_all(iter,41),1),C( Route_best_all(iter,1),1)],[C( Route_best_all(iter,41),2),C( Route_best_all(iter,1),2)],'ro-');
    title(['最短路径：',num2str(Length_best(iter))]);
    hold off;
    pause(0.05);
    
     %迭代次数+1，清空记录路线
      iter = iter + 1
      Table = zeros(m,p+4);                  % 路径记录表
      Table_all=zeros(m,p+4);              %完整的路线记录表
      num=ones(m,39);
    
end

在重要位置预设了车库、垃圾处理中心的位置，运行到这部分的时候，要注意不能让其变换掉。并且每一段路用不同颜色进行了绘制，方便观察。得到这样的图。

在这里进行了结果显示。

%% 5. 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best_all(index,:);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str([Shortest_Route Shortest_Route(1)])]);

把每次迭代的路径距离进行作图。

figure(1)
subplot(2,2,4);
%plot(1:iter_max,Length_best,'b',1:iter_max,Length_ave,'r:')
plot(1:iter_max,Length_best,'b')
legend('最短距离')
%xlabel('迭代次数')
%ylabel('距离')
title('最优历史路径长度')

得到这样的图。

再多迭代几次看。

观察到第三个图的那个路线，个人感觉最短距离还可以更短，本想再引入“随机”考虑的方向，觉得问题不是出在这里，能够得到更优化的代码后，再在这里更新。