matlab蚁群算法的优化计算——旅行商问题（TSP）优化
蚁群算法（ ant colony algorithn,ACA）是由意大利学者M. dorigo等人于20世纪90年代初提出的一种新的模拟进化算法，其真实地模拟了自然界蚂蚁群体的觅食行为。M. Dorigo等人将其用于解决旅行商问题（ traveling salesman problem,TSP），并取得了较好的实验结果。
近年来，许多专家学者致力于蚁群算法的研究，并将其应用于交通、通信、化工、电力等领域，成功解决了许多组合优化问题，如调度问题（job- shop scheduling problem）、指派问题（ quadratic assignment problem）、旅行商问题（ traveling salesman problem）等。
将详细阐述蚁群算法的基本思想及原理，并以实例的形式介绍其应用于解决中国旅行商问题（ chinese tsp,CTSP）的情况。
蚁群算法
基本思想生物学家研究发现，自然界中的蚂蚁觅食是一个群体性行为，并非单只蚂蚁自行寻找食物源。蚂蚁在寻找食物源时，会在其经过的路径上释放一种信息素，并能够感知其他蚂蚁释放的信息素。信息素浓度的大小表征路径的远近，信息素浓度越高，表示对应的路径距离越短。通常，蚂蚁会以较大的概率优先选择信息素浓度较高的路径，并释放一定量的信息素，以增强该条路径上的信息素浓度，这样会形成一个正反馈。最终，蚂蚁能够找到一条从巢穴到食物源的最佳路径，即最短距离。值得一提的是，生物学家同时发现，路径上的信息素浓度会随着时间的推进而逐渐衰减。
将蚁群算法应用于解决优化问题的基本思路为：用蚂蚁的行走路径表示待优化问题的可行解，整个蚂蚁群体的所有路径构成待优化问题的解空间。路径较短的蚂蚁释放的信息素量较多，随着时间的推进，较短的路径上累积的信息素浓度逐渐增高，选择该路径的蚂蚁个数也愈来愈多。最终，整个蚂蚁会在正反馈的作用下集中到最佳的路径上，此时对应的便是待优化问题的最优解。
蚁群算法解决TSP问题基本原理
本节将用数学语言对上述蚁群算法的基本思想进行抽象描述，并详细阐述蚁群算法用于解决TSP问题的基本原理。
不失一般性，设整个蚂蚁群体中蚂蚁的数量为m，城市的数量为n，城市i与城市j之间的距离为d3（i，=1,2，…，n），t时刻城市i与城市j连接路径上的信息素浓度为τ（t）。初始时刻，各个城市间连接路径上的信息素浓度相同，不妨设为τ（0）=τo。
蚂蚁k（k=1,2，…，m）根据各个城市间连接路径上的信息素浓度决定其下一个访问城市，设P（t）表示t时刻蚂蚁k从城市i转移到城市j的概率，其计算公式为
其中，v（t）为启发函数，v（t）=1/d，表示蚂蚁从城市i转移到城市j的期望程度； allow（k=1,2，…，m）为蚂蚁k待访问城市的集合，开始时，allw中有（n-1）个元素，即包括除了蚂蚁k出发城市的其他所有城市，随着时间的推进， allowr中的元素不断减少，直至为空，即表示所有的城市均访问完毕；a为信息素重要程度因子，其值越大，表示信息素的浓度在转移中起的作用越大；B为启发函数重要程度因子，其值越大，表示启发函数在转移中的作用越大，即蚂蚁会以较大的概率转移到距离短的城市。
如前文所述，在蚂蚁释放信息素的同时，各个城市间连接路径上的信息素逐渐消失，设参数p（0<ρ<1）表示信息素的挥发程度。因此，当所有蚂蚁完成一次循环后，各个城市间连接路径上的信息素浓度需进行实时更新，即
其中，Δτ表示第k只蚂蚁在城市i与城市j连接路径上释放的信息素浓度；△r；表示所有蚂蚁在城市i与城市j连接路径上释放的信息素浓度之和针对蚂蚁释放信息素问题，M. dorigo等人曾给出3种不同的模型，分别称之为 ant cycle system、 ant quantity system和 ant density system，
基本步骤
1.初始化参数在计算之初，需要对相关的参数进行初始化，如蚁群规模（蚂蚁数量）m、信息素重要程度因子a、启发函数重要程度因子、信息素挥发因子p信息素释放总量Q、最大迭代次数Iter max迭代次数初值 I Iter=1。
2.构建解空间将各个蚂蚁随机地置于不同出发点，对每个蚂蚁k（k=1,2，…，m），按照式（22-1）计算其下一个待访问的城市，直到所有蚂蚁访问完所有的城市。
3.更新信息素计算各个蚂蚁经过的路径长度L4（k=1,2，…，m），记录当前迭代次数中的最优解（最短路径）。同时，根据式（22-2）和式（22-3）对各个城市连接路径上的信息素浓度进行更新。
4.判断是否终止若iter<iter_max，则令iter=iter+1，清空蚂蚁经过路径的记录表，并返回步骤2；否则，终止计算，输出最优解。
蚁群算法的特点
基于蚁群算法的基本思想及解决TSP问题的基本原理，不难发现，与其他优化算法相比，蚁群算法具有以下几个特点：
（1）采用正反馈机制，使得搜索过程不断收敛，最终逼近最优解。
（2）每个个体可以通过释放信息素来改变周围的环境，且每个个体能够感知周围环境的实时变化，个体间通过环境进行间接地通讯。
（3）搜索过程采用分布式计算方式，多个个体同时进行并行计算，大大提高了算法的计算能力和运行效率。
（4）启发式的概率搜索方式不容易陷入局部最优，易于寻找到全局最优解。
基本步骤
1.计算城市间相互距离根据城市的位置坐标，计算两两城市间的相互距离，从而得到对称的距离矩阵（维数为31的方阵）。需要说明的是，计算出的矩阵对角线上的元素为0，然而如前文所述，由于启发函数v（D）=1d，因此，为了保证分母不为零，将对角线上的元素修正为一个非常小的正数（如10-4或10-5等）。
2.初始化参数如前文所述，在计算之前，需要对相关的参数进行初始化，此处不再赘述。具体请参考本章22.3.4节程序实现部分。
3.迭代寻找最佳路径首先构建解空间，即各个蚂蚁根据转移概率公式（22-1）访问所有的城市。然后计算各个蚂蚁经过路径的长度，并在每次迭代后根据式（22-2）和式（22-3）实时更新各个城市连接路径上的信息素浓度。经过循环迭代，记录下最优的路径及其长度。
4.结果分析找到最优路径后，可以将之与其他方法得出的结果进行比较，从而对蚁群算法的性能进行评价。同时，也可以探究不同取值的参数对优化结果的影响，从而找到一组最佳或者较佳的参数组合。
主程序展示
%% 第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);
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;
end
end
end
%% 初始化参数
m = 50; % 蚂蚁数量
alpha = 1; % 信息素重要程度因子
beta = 5; % 启发函数重要程度因子
rho = 0.1; % 信息素挥发因子
Q = 1; % 常系数
Eta = 1./D; % 启发函数
Tau = ones(n,n); % 信息素矩阵
Table = zeros(m,n); % 路径记录表
iter = 1; % 迭代次数初值
iter_max = 200; % 最大迭代次数
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;
end
P = P/sum(P);
% 轮盘赌法选择下一个访问城市
Pc = cumsum(P);
target_index = find(Pc >= rand);
target = allow(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)])]);
%% 绘图
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>数据集
1304 2312
3639 1315
4177 2244
3712 1399
3488 1535
3326 1556
3238 1229
4196 1004
4312 790
4386 570
3007 1970
2562 1756
2788 1491
2381 1676
1332 695
3715 1678
3918 2179
4061 2370
3780 2212
3676 2578
4029 2838
4263 2931
3429 1908
3507 2367
3394 2643
3439 3201
2935 3240
3140 3550
2545 2357
2778 2826
2370 2975结果展示
与运行结果对应的路径如图22-4所示。从图中可以清晰地看到，自起点出发，每个城市访问一次，遍历所有城市后，返回起点。寻找到的最短路径为15601.9195km。
各代的最短距离与平均距离如图所示，从图中不难发现，随着迭代次数的增加，最短距离与平均距离均呈现不断下降的趋势。当迭代次数大于112时，最短距离已不再变化，表示已经寻找到最佳路径。
最新研究成果表明，中国TSP问题的最优解为15377km，因此，这里寻找到的最佳路径是局部最优解，而并非全局最优解。
