转载声明：感谢：解放军信息工程大学某老师
尊重原作者劳动：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问题
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