【matlab】蚁群算法详解

    
    if  NC >= 2
        L_tmp = 0;
        for j = 1:(n-1)
            L_tmp =  L_tmp + D( Tabu(1,j), Tabu(1,j+1) );
        end
        L_tmp = L_tmp + D( Tabu(1,1), Tabu(1,n) );
        if L_tmp > L_best(NC-1,:)
            Tabu(1,:) = R_best(NC-1,:);     % 正反馈机制
        end
    end
    %{
        此处正反馈机制不可或缺，实际上，每次迭代的最短距离是变化的，
        并非最优解，为了最短距离和最佳路线保持全局收缩，将上一次的
        最佳路线赋予下一次尤为重要！
    %}
    

function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length] ...
        = ACATSP(C,NC_max,m,Alpha,Beta,Rho,Q)
%% -------------------------------------------------------------
%  ACATSP.m 【Ant Colony Algorithm for travellingh Salesman Problem】
%  基于蚁群算法的旅行商问题求解
%  基本思路：栅格， 轮盘法则 
%  感谢：解放军信息工程大学某老师的无私。
%        尊重原作者劳动：http://blog.sina.com.cn/s/blog_5013f7e30100aodx.html
%  同谢：网友wayjj. 博客 http://blog.csdn.net/wayjj/article/details/72809344
%% -------------------------------------------------------------
%  输入参数列表
%  C                n 个城市的坐标，n×2的矩阵
%  NC_max           最大迭代次数
%  m                蚂蚁个数
%  Alpha            表征信息素重要程度的参数
%  Beta             表征启发式因子重要程度的参数
%  Rho              信息素蒸发系数
%  Q                信息素增加强度系数
%   
%  输出参数列表
%  R_best           各代最佳路线
%  L_best           各代最佳路线的长度
%  L_ave            ？
%  Shortest_Route   最短路线
%  Shortest_Length  最短路径长度
%% --------------------变量初始化--------------------------------
tic;
n = size(C,1);                  % n 表示问题的规模（城市个数）
D = zeros(n,n);                 % D 表示完全图的赋权邻接矩阵
for i = 1:n
    for j = 1:n
        if i ~= j
            D(i,j) = ((C(i,1)-C(j,1))^2 + (C(i,2)-C(j,2))^2)^0.5;
        else
            D(i,j) = eps;       % i=j时不计算，应该为0，但后面的启发因子要...
        end                     % 取倒数，用eps（浮点相对精度）表示
        D(j,i) = D(i,j);        % 考虑对称TSP问题
    end
end

Eta   = 1./D;                   % Eta为启发信息矩阵，城市间直线距离的倒数
Tau   = ones(n,n);              % Tau为信息素矩阵
Tabu  = zeros(m,n);  % m只蚂蚁分别遍历 n座城市 存储并记录路径的生成，禁忌表
NC    = 1;                      % 迭代计数器，记录迭代次数
R_best = zeros(NC_max,n);       % 各代最佳路线存储
L_best = inf .* ones(NC_max,1); % 各代最佳路线的长度
L_ave  = zeros(NC_max,1);       % 各代路线的平均长度

%% -------------------  启动蚂蚁觅食活动 ------------------------
while NC <= NC_max                  % 停止条件之一：迭代达到最大值
%%  第二步： 将m只蚂蚁放到n个城市上
    Randpos = [];                   % 随机存取
    for i = 1:ceil(m/n)
        Randpos = [Randpos, randperm(n)];
        %{Randpos(n) = [1-n随机分布]%}
    end
    Tabu(:,1) = ( Randpos(1,1:m) )';% 取m个数据作为m只蚂蚁，均从1号城市出发
    
    
%%  第三步：m 只蚂蚁按概率函数选择下一座城市，完成各自的周游
    for j = 2:n                         % 第一站已经随机产生
        for i = 1:m                     % 蚂蚁编号
            visited = Tabu( i,1:(j-1) );% 记录第i只蚂蚁已访问的城市，避免重复访问
            J  = zeros(1,(n-j+1));      % 待访问的城市init
            P  = J;                     % 待访问城市的选择概率分布
            Jc = 1;                     % 访问的城市个数
            for k = 1:n                 % 存储未访问的城市
                if length( find(visited==k) )==0
                    J(Jc) = k;          % J 记录仍未访问的城市
                    Jc    = Jc+1;
                end
            end
            
            % 下面计算待选城市的概率分布
            for k = 1:length(J)
                P(k) = ( Tau(visited(end), J(k))^Alpha ) * ...
                    ( Eta(visited(end), J(k))^Beta  );
            end
            P = P/sum(P);                   % 状态转移公式
            
            % 按概率原则选取下一个城市【算法核心，大有文章可搞】
            %%     ★  ★  ★  ★
            Pcum = cumsum(P);               % example->cumsum([1 1 3]) = [1 2 5]
            %{
           有一点要特别说明，用到cumsum(P)，蚂蚁要选择的下一个城市不是按
           最大概率，使用轮盘法则，保证不影响全局收缩
            %}
            Select = find( Pcum >= rand );  % 若Pcum中某概率>=某随机数，返回其索引
            
            %要选择其中总概率大于等于某一个随机数，找到大于等于这个随机数的城市的在J中的位置
            to_visit = J( Select(1) );      % Select有可能有多个，默认取第1个
            Tabu(i,j) = to_visit;           % i 当前蚂蚁，j 下一站城市
            %{
            ★ 纠正一下轮盘法则，前面博主认为此处基于轮盘法则Select(1)，1 保证可以选到最大概率的城市，事实上并不能 ！
            ★ 个人举例：p=[0.1 0.2 0.3 0.4] 最后那个城市概率最大
            此时Pcum=[0.1 0.3 0.6 1], rand随机，特举= 0.223423；故Select =[2 3 4]; Select(1) = 2, 并不能选到最大概率的城市！
            【 欢迎其他兴趣读者提出 质疑 或 见解 ！ 】
            %}
            %%     ★  ★  ★  ★
        end
    end
    
    if NC >= 2                        
        Tabu(1,:) = R_best(NC-1,:);    
    end
    
%%  第四步：记录本次迭代最佳路线
    L = zeros(m,1);                     % 开始距离为0，m*1的列向量
    %{
    L=zeros(m,1)   记录本次迭代最佳路线的长度，即每只蚂蚁遍历所有城市
                   所走路径距离
    %}
    for i = 1:m                         % 蚂蚁编号
        R = Tabu(i,:);                  % R 存储本次迭代第i只蚂蚁的最佳路线
        for j = 1:(n-1)
            L(i) = L(i) + D(R(j), R(j+1)); % 计算第i只蚂蚁本次迭代路径长度
        end
        L(i) = L(i) + D( R(1),R(n) );      % 首尾相连后才遍历所有城市，总距离
    end
    
    L_best(NC) = min(L);                % 第 NC次迭代 最短距离
    pos = find( L == L_best(NC) );      % 第 NC次迭代 路径最短的蚂蚁编号，可能出现两只以上，选1即可
    R_best(NC,:) = Tabu(pos(1),:);      % 第 NC次迭代 的最佳路线图
    %{
    R_best(NC,:)=Tabu(pos(1),:)：找到路径最短的那条蚂蚁所在的城市先后顺序，
                                 pos(1)中1表示万一有长度一样的两条蚂蚁，
                                 那就选第一个
    %}
    L_ave(NC) = mean(L);                % 第 NC次迭代 的平均距离
    NC = NC+1;                          % 迭代继续
    
%%  第五步：更新信息素
    Delta_Tau = zeros(n,n);             % 开始时信息素为 n*n的0矩阵
    % n 座城市，n*n种可能连线，每条线来个信息素初始化
    
    for i = 1:m                         % 蚂蚁编号
        for j = 1:(n-1)                 %
            Delta_Tau( Tabu(i,j),Tabu(i,j+1) ) = ...
                Delta_Tau( Tabu(i,j),Tabu(i,j+1) ) + Q/L(i); % ★ 公式 ★
            % 此次循环在路径（i，j）上的信息素增量
        end
        Delta_Tau( Tabu(i,n),Tabu(i,1) ) = ....
            Delta_Tau( Tabu(i,n),Tabu(i,1) ) + Q/L(i); % 首尾，不要漏掉
        % 此次循环在整个路径上的信息素增量
    end
    Tau = (1-Rho) .* Tau + Delta_Tau;   % 考虑信息素挥发，更新后的信息素

%%  第六步：整个流程全部走完！禁忌表清零，Ready For Next Time
    Tabu = zeros(m,n);                  % 直到最大迭代次数
end

%% ★★★ 手工封装一个画图函数 ★★★
%{
    http://blog.sina.com.cn/s/blog_a74f6fe7010162bh.html
    M文件中将函数的调用直接写到m脚本文件中的情况是不允许的，必须也把调用
    写成函数的形式，或者将子函数都写成单独的m文件。
%}
    function DrawRoute(C,R)
        N_D = length(R);
        scatter( C(:,1),C(:,2) );
        hold on;
        plot([C(R(1),1),C(R(N_D),1)], [C(R(1),2),C(R(N_D),2)], 'r');
        hold on;
        for i_D = 2:N_D
            plot([C(R(i_D-1),1), C(R(i_D),1)], [C(R(i_D-1),2), C(R(i_D),2)], 'r');
            hold on;
        end
        title('旅行商问题优化结果');
    end

%%  第七步：输出 迭代了NC次后的最佳路径相关信息
    Pos = find( L_best == min(L_best) );    % 找到最佳路径（非0为真）
    Shortest_Route = R_best(Pos(1),:);      % 最大迭代次数后最佳路径图
    Shortest_Length = L_best(Pos(1));       % 最大迭代次数后最短距离
    subplot(1,2,1);  DrawRoute(C,Shortest_Route);  % 画路线图的子函数
    subplot(1,2,2);  
    plot(L_best);hold on;plot(L_ave,'r');          % 画平均距离和最短距离     
    title('平均距离和最短距离');
end