%% 第22章 蚁群算法的优化计算——旅行商问题(TSP)优化
% <html>
% <table border="0" width="600px" id="table1"> <tr> <td><b><font size="2">该案例作者申明：</font></b></td> </tr> <tr><td><span class="comment"><font size="2">1：本人长期驻扎在此<a target="_blank" href="http://www.matlabsky.com/forum-78-1.html"><font color="#0000FF">板块</font></a>里，对该案例提问，做到有问必答。</font></span></td></tr><tr> <td><span class="comment"><font size="2">2</font><font size="2">：此案例有配套的教学视频，视频下载请点击<a href="http://www.matlabsky.com/forum-91-1.html">http://www.matlabsky.com/forum-91-1.html</a></font><font size="2">。 </font></span></td> </tr> <tr> <td><span class="comment"><font size="2"> 3：此案例为原创案例，转载请注明出处（《MATLAB智能算法30个案例分析》）。</font></span></td> </tr> <tr> <td><span class="comment"><font size="2"> 4：若此案例碰巧与您的研究有关联，我们欢迎您提意见，要求等，我们考虑后可以加在案例里。</font></span></td> </tr> <tr> <td><span class="comment"><font size="2"> 5：以下内容为初稿，与实际发行的书籍内容略有出入，请以书籍中的内容为准。</font></span></td> </tr> </table>
% </html>
%% 清空环境变量
clear all
clc
%% 导入数据
load citys_data.mat
%% 计算城市间相互距离
n = size(citys,1); %citys 是一个31行2列的data数据，是matlab中的一个数据存储文件，类似于excel。这里的size(citys,1)表示的是只取citys中第一列的数据，n=size(citys,1)中的n为size(citys,1)，31行1列中的列数。
D = zeros(n,n);
for i = 1:n
for j = 1:n
if i ~= j
D(i,j) = sqrt(sum((citys(i,:) - citys(j,:)).^2));
else
D(i,j) = 1e-4; %因为对角线元素都是0.0001，然后矩阵D又表示成1e+3，即乘以1000，所以对角线上的数字就变成1e-7，即0.0000001。多以四舍五入，取小数点后四位，就变成了0.000
end
end
end
%% 初始化参数
% m = 50; % 蚂蚁数量
m = 5;
alpha = 1; % 信息素重要程度因子
beta = 5; % 启发函数重要程度因子
rho = 0.1; % 信息素挥发因子
Q = 1; % 常系数
Eta = 1./D; % 启发函数 %1./D表示1除以D中的所有数字，因为D中除了对角线的元素，都是1.0e+03*,即除了对角线以外的元素都是几千，1除以几点得到的是0.0000几，取小数点后四位四舍五入，得到的就是0.0000。而对角线上的元素本来就是1e-4，即0.00001，用1除以0.00001得到的就是就是10000，所以对角线元素就是1，然后乘1e+4。
Tau = ones(n,n); % 信息素矩阵
Table = zeros(m,n); % 路径记录表
iter = 1; % 迭代次数初值
% iter_max = 200; % 最大迭代次数
iter_max = 10;
Route_best = zeros(iter_max,n); % 各代最佳路径
Length_best = zeros(iter_max,1); % 各代最佳路径的长度
Length_ave = zeros(iter_max,1); % 各代路径的平均长度
%% 迭代寻找最佳路径
while iter <= iter_max
% 随机产生各个蚂蚁的起点城市
start = zeros(m,1);
for i = 1:m
temp = randperm(n);
start(i) = temp(1);
end
Table(:,1) = start;
% 构建解空间
citys_index = 1:n;
% 逐个蚂蚁路径选择
for i = 1:m
% 逐个城市路径选择
for j = 2:n
tabu = Table(i,1:(j - 1)); % 已访问的城市集合(禁忌表)
allow_index = ~ismember(citys_index,tabu);
allow = citys_index(allow_index); % 待访问的城市集合
P = allow;
% 计算城市间转移概率
for k = 1:length(allow)
P(k) = Tau(tabu(end),allow(k))^alpha * Eta(tabu(end),allow(k))^beta; %P(k)的等号右边是一个公式，用来计算转移概率。其中Tau的行和列表示的都是城市，以第一行为例，表示第一个城市到1,2,3城市的信息素。
end
P = P/sum(P);
% 轮盘赌法选择下一个访问城市
Pc = cumsum(P);
target_index = find(Pc >= rand);
target = allow(target_index(1)); % Pc中可能产生两个概率，都大于rand的情况，取得是第一个大于rand所对应的城市，所以target_index(1)
Table(i,j) = target;
end
end
% 计算各个蚂蚁的路径距离
Length = zeros(m,1);
for i = 1:m
Route = Table(i,:);
for j = 1:(n - 1)
Length(i) = Length(i) + D(Route(j),Route(j + 1));
end
Length(i) = Length(i) + D(Route(n),Route(1));
end
% 计算最短路径距离及平均距离
if iter == 1
[min_Length,min_index] = min(Length);
Length_best(iter) = min_Length;
Length_ave(iter) = mean(Length);
Route_best(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(iter,:) = Table(min_index,:); %如果本次的最短距离就是本次的最优最短距离，则本次最优的最短路径就是本次最短距离所对应的路径
else
Route_best(iter,:) = Route_best((iter-1),:); %如果本次的最短距离不是本次的足有最优最短距离，则本次的最优的最短路径就是上一次迭代所得到的最短距离所对应的路径
end
end
% 更新信息素
Delta_Tau = zeros(n,n);
% 逐个蚂蚁计算
for i = 1:m
% 逐个城市计算
for j = 1:(n - 1)
Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i);
end
Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i);
end
Tau = (1-rho) * Tau + Delta_Tau;
% 迭代次数加1，清空路径记录表
iter = iter + 1;
Table = zeros(m,n);
end
%% 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best(index,:);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str([Shortest_Route Shortest_Route(1)])]); %括号里之所以有Shortest_Route(1),表示的是最终要回到原点
%% 绘图
figure(1)
plot([citys(Shortest_Route,1);citys(Shortest_Route(1),1)],...
[citys(Shortest_Route,2);citys(Shortest_Route(1),2)],'o-');
grid on
for i = 1:size(citys,1)
text(citys(i,1),citys(i,2),[' ' num2str(i)]);
end
text(citys(Shortest_Route(1),1),citys(Shortest_Route(1),2),' 起点');
text(citys(Shortest_Route(end),1),citys(Shortest_Route(end),2),' 终点');
xlabel('城市位置横坐标')
ylabel('城市位置纵坐标')
title(['蚁群算法优化路径(最短距离:' num2str(Shortest_Length) ')'])
figure(2)
plot(1:iter_max,Length_best,'b',1:iter_max,Length_ave,'r:')
legend('最短距离','平均距离')
xlabel('迭代次数')
ylabel('距离')
title('各代最短距离与平均距离对比')
%%
% <html>
% <table width="656" align="left" > <tr><td align="center"><p align="left"><font size="2">相关论坛：</font></p><p align="left"><font size="2">Matlab技术论坛：<a href="http://www.matlabsky.com">www.matlabsky.com</a></font></p><p align="left"><font size="2">M</font><font size="2">atlab函数百科：<a href="http://www.mfun.la">www.mfun.la</a></font></p></td> </tr></table>
% </html>
