【matlab】蚁群算法详解

转载声明：感谢：解放军信息工程大学某老师
尊重原作者劳动：http://blog.sina.com.cn/s/blog_5013f7e30100aodx.html
同谢：网友wayjj. 博客 http://blog.csdn.net/wayjj/article/details/72809344

作者重新排版并重新注释，水平不到家，欢迎指正！画图部分博主原创，有所保留。

17/9/13 10.06  第1次补丁，修正前人错误，关于正反馈，补丁如下：

        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问题    endendEta   = 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