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