基于蚁群算法的MTSP问题

蚁群算法

一、简介

20世纪50年代中期创立了仿生学，人们从生物进化的机理中受到启发。提出许多用以解决复杂优化问题的新方法，如进化规划、进化策略、遗传算法等，这些算法成功地解决了一些实际问题。而同时一种对自然界蚂蚁的寻径方式进行模拟而得出的仿生算法——蚁群算法也应运而生。

二、思想与简易模型

现考察如下图所示的简单情况：

我们可以清晰的看到，在蚁穴和食物之间存在两条路径，其中上路径较短，而下路径较长。在对蚂蚁的实际考察中，发现蚂蚁总会以一条最优的路径在蚁穴与食物之间进行运输，而蚂蚁又是如何判断出最优路径的呢？
研究发现，蚂蚁在寻找路径时会存在着释放信息素和受信息素导引的现象。所谓释放信息素是指蚂蚁在一条路径上行走时，会在该路径上释放信息素来指引其他蚂蚁的行径;而受信息素导引是指，蚂蚁在走过路口时，选择信息素浓度较高的路径概率会相对较大，下面我会以最简单的情况描述为什么短路径上的信息素浓度较大。

假设：
1.初始时，蚂蚁从A点出发，以等概率选择路径ABD或ACD到达食物点D，然后原路返回。
2.每只蚂蚁的行走速率相同，即每单位时间行走一步。
3.假设ABD路径距离为8步，ACD的路径距离为16步。
4.假设蚂蚁每走一步释放一个信息素单位，且不考虑信息素随时间的流失。

由假设1，我们可以做出简单等效，初始时，每条路径上派一只蚂蚁。（蚂蚁数据量较大时，数量分布近似概率分布，我们可以进行归一化，假设初始每条路径只有一只蚂蚁）。上图显示的是第8个时间单位的情况，可以看到，选择路径ABD的蚂蚁已经达到食物源，而选择路径ACD的蚂蚁到达路径的中点C。此时的信息素分布如下图所示：

经过16个单位时间后，选择路径ABD的蚂蚁已往返了一趟，而选择ACD的蚂蚁刚刚到达食物源D。此时信息素分布如下图所示：

经过32个时间单位后，选择路径ABD与路径ACD的蚂蚁都已经回到了蚁穴，此时选择路径ABD的蚂蚁已经往返了两趟，而选择路径ACD的蚂蚁往返了一趟，此时信息素分布如下：

寻找食物的过程继续，此时两只蚂蚁已经取得了食物而且回到了蚁穴，但此时两条路线上信息素分布产生了差异，由上图可知ABD路线上每一处的信息素数量与ACD路线上每一处的信息素数量之比为2:1。这意味着蚂蚁选择ABD路径的概率大于选择ACD路径概率，我们姑且假设选择路径概率与信息素浓度呈线性关系，则蚂蚁选择ABD路径的概率与选择ACD路径概率之比为2:1。
等效的做法是，我们在ABD路径上增派一只蚂蚁，则在第64个时间单位时的信息素分布如下：

由上图可见，第64个时间单位时，ABD路线上每一处的信息素数量与ACD路线上每一处的信息素数量之比为3:1。假设我们将此过程无线进行下去，则按照信息素的指导，两条路径上信息素浓度之比将越来越大，且变化越来越快，形成一种正反馈，最终将导致所有的蚂蚁都会放弃ACD路线，而都选择ABD路线。

三、TSP问题的算法思路

参考链接：
蚁群算法(含matlab代码)
蚁群算法
matlab—基于蚁群算法求解Tsp旅行商问题
运用蚁群算法求解TSP问题时的基本原理是：将m个蚂蚁随机地放在多个城市，让这些蚂蚁从所在的城市出发，n步（一个蚂蚁从一个城市到另外一个城市为1步）之后返回到出发的城市。如果m个蚂蚁所走出的m条路经对应的中最短者不是TSP问题的最短路程，则重复这一过程，直至寻找到满意的TSP问题的最短路径为止。为了说明这一个算法下面用一个算法流程图来表示一下：

（一）算法的相关参数
在算法中用到的变量比较多，而每个变量的对算法最终收敛的结果影响程度也不一样，在这里把一些重要的参数说明一下。

τij\tau _{ij}τij：城市i和城市j之间边eije_{ij}eij上信息素的残留强度

Δτij\Delta\tau _{ij}Δτij：一次循环后边eije_{ij}eij上信息素的增量

Δτijk\Delta \tau_{ij}^{k}Δτijk：一次循环之后蚂蚁k对边eije_{ij}eij上信息素的贡献量

ηij\eta _{ij}ηij：城市i，j之间的能见度，反映了由城市i转移到城市j的启发程度

dijd_{ij}dij：城市i到城市j之间的距离

pijkp_{ij}^{k}pijk：蚂蚁k从当前所在的城市到下一个城市去的概率

tabu(k)tabu(k)tabu(k)：禁忌表，用于存放第k只蚂蚁已经走过的城市

QQQ：信息素总量，信息素总信息量为蚂蚁循环一周后向经过路径释放信息素的总量

ρ\rhoρ：信息素残留系数，由于蚂蚁释放的信息量回随着时间的转移而逐渐挥发，以至于路径上的信息素不能无限递增，该系数太小时会降低算法的全局搜素能力，过大时容易使算法陷入局部最优，影响全局搜素能力。

mmm：蚂蚁总数，在TSP问题中，每次循环当中，每只蚂蚁所走出的每条路径为TSP问题的候选解，m只蚂蚁一次循环所走出来的m条路经为TSP问题的一个解子集，所以这个解子集越大则算法的全局搜索能力越强，但是过大会使算法的收敛速度降低。如果m太小的话，算法也很容易就陷入局部最优，过早的出现停滞现象，所以选择合理的蚂蚁总数是非常重要的。

α\alphaα：信息启发因子，反映了蚂蚁在从城市i向城市j移动时，这两个城市之间道路上所累积的信息素在指导蚂蚁选择城市j的程度，即蚁群在路径搜素中随即性因素作用的强度。

β\betaβ：期望值启发式因子，反映了蚂蚁在从城市i向城市j转移时候期望值ηij\eta _{ij}ηij在指导蚁群搜素中的相对重要程度。其大小反映了蚁群在道路搜素中的先验性、确定性等因素的强弱，α\alphaα、β\betaβ的大小也会影响算法的收敛性。

(二) 算法过程中的关键步骤
1、上面讲了蚂蚁在选择下一个要转移的城市时候是基于概率选择的，当然这个概率不是随机概率，而是其中的一些参数来决定。假定在t时刻，蚂蚁从目前所在的i城市要选择转移去的下一个j城市，概率大小为pijkp_{ij}^{k}pijk：
pijk(t)={τijkηijβ∑s∈allowedkτisk(t)ηisβ(t)j∈allowedk0otherwisep_{ij}^{k}(t)=\left\{\begin{matrix} & \frac{\tau_{ij}^{k}\eta _{ij}^{\beta }}{\sum_{s\in allowed_{k}}\tau_{is}^{k}(t)\eta_{is}^{\beta}(t)} &j\in allowed_{k} \\ \\ &0 &otherwise \end{matrix}\right.pijk(t)=⎩⎪⎨⎪⎧∑s∈allowedkτisk(t)ηisβ(t)τijkηijβ0j∈allowedkotherwise

其中allowedkallowed_{k}allowedk表示允许蚂蚁k下一步可容许去的城市的集合，τijα\tau_{ij}^{\alpha}τijα为边eije_{ij}eij上的信息素因数，ηijβ\eta_{ij}^{\beta}ηijβ为城市i，j间能见度因数。至于这个过程具体是怎么实现的，在程序中会有相关的源码。

2、对任意两个城市i，j之间道路对应的边eije_{ij}eij信息素增量按照下式进行：
{τij(t+1)=ρτij(t)+Δτij(t,t+1)Δτij(t,t+1)=Σk=1mΔτijk(t,t+1)\left\{\begin{matrix} &\tau_{ij}(t+1)=\rho\tau_{ij}(t)+\Delta \tau_{ij}(t,t+1)\\ \\&\Delta\tau_{ij}(t,t+1)=\Sigma_{k=1}^{m}\Delta\tau_{ij}^{k}(t,t+1) \end{matrix}\right.⎩⎨⎧τij(t+1)=ρτij(t)+Δτij(t,t+1)Δτij(t,t+1)=Σk=1mΔτijk(t,t+1)

其中，Δτijk(t,t+1)\Delta\tau_{ij}^{k}(t,t+1)Δτijk(t,t+1)为蚂蚁k对边eije_{ij}eij上所贡献的信息素增量，Δτij(t,t+1)\Delta\tau_{ij}(t,t+1)Δτij(t,t+1)是经过边eije_{ij}eij的所有蚂蚁对边eije_{ij}eij的信息素量贡献，ρ\rhoρ为信息素残留系数。对于Δτijk(t,t+1)\Delta\tau_{ij}^{k}(t,t+1)Δτijk(t,t+1)的调整方式不同，可以将蚁群算法分为三种模型：蚁密模型、蚁量模型和蚁周模型。

2.1 蚁密模型
在蚁密模型当中，每一只蚂蚁在经过城市i,j时，对该边eije_{ij}eij所贡献的信息素增量为常量，每个单位长度是Q：

2.2 蚁量模型
在蚁量模型中，一只蚂蚁k在经过城市i,j之间边eije_{ij}eij时，对该边所贡献的信息素增量为变量，即每单位长度上为Q/dijd_{ij}dij，它与城市i,j之间的路径长度dijd_{ij}dij有关，具体为：

2.3 蚁周模型
上述两种模型，对两城市之间边eije_{ij}eij上信息素贡献的增量在蚂蚁k经过边的同时完成，而蚁周模型对边eije_{ij}eij信息素的增量是在本次循环结束时才进行更新调整。一只蚂蚁在经过城市i，j时对边eije_{ij}eij上贡献的增量为每单位长度Q/LkQ/L_{k}Q/Lk，LkL_{k}Lk为蚂蚁在本次循环走出路径的长度。

本文章就是基于蚁周模型来实现的，介绍完上述的几个关键步骤以后，就要用matlab来实现MTSP问题了，另外，关于TSP文体的算法主体框架可参考上述链接中的三篇文章，本文不再详述。

四、MTSP算法实现

本文的重点是运用蚁群算法解决MTSP问题，而MTSP问题的算法框架与TSP问题的算法框架基本相似，但关于蚂蚁单位的处理有一些需要注意的细节，首先让我们考察一个具体的问题。附件中列出了陕西曲线公路里程表，现在有两个旅行者从碑林区出发，求合计拜访每个区县一次且仅一次且最后都回到碑林区的最短路径。
和单纯的TSP问题不同的是，我们的单位不再是一只蚂蚁，而是一对蚂蚁，我们通过一对蚂蚁进行全局搜索。这里便涉及到几个TSP问题不需要考虑到的问题，
问题一：如何让两只蚂蚁协同搜索路径，即每进行一个时间单位，应该是哪一只蚂蚁进行搜索，亦或是两只蚂蚁一起搜索。
问题二：两只蚂蚁释放的信息素是否为同一种，即确定下一步搜索路径时，是要将两种蚂蚁释放的信息素综合起来考虑，还是应该分开考虑两种信息素的影响.
解决方法：
1.若每次让两只蚂蚁一起搜素，则每只蚂蚁会固定的走访一半的城市，而往往最优的路径并不是每只蚂蚁一定要走访一半的城市。所以我们选择对两只蚂蚁即将走访的两条路径，选择路径较短的蚂蚁进行下一步的路径搜索。
2.经过试验发现，信息素选择为一种与两种的实验结果影响不大，故将两种蚂蚁释放的信息素综合起来考虑。

多旅行商问题的代码如下所示:
(具体的实现细节详见代码注释)

D = P3;
n = 107;
m = 50;                              % 蚂蚁对数，每对蚂蚁分为蚂蚁1与蚂蚁2
alpha = 1;                           % 信息素重要程度因子
beta = 5;                            % 启发函数重要程度因子
rho = 0.1;                           % 信息素挥发因子
Q = 1;                               % 常系数（信息素释放量）
Eta = 1./D;                          % 启发函数
Tau = ones(n,n);                     % 信息素矩阵
Table1 = zeros(m,n);                 % 蚂蚁1路径记录表
Table2 = zeros(m,n);                 % 蚂蚁2路径记录表
iter = 1;                            % 迭代次数初值
iter_max = 150;                      % 最大迭代次数
Route_best1 = zeros(iter_max,n);     % 蚂蚁1各代最佳路径
Route_best2 = zeros(iter_max,n);     % 蚂蚁2各代最佳路径
Length_best = zeros(iter_max,1);     % 各代最佳路径的长度
Length_ave  = zeros(iter_max,1);     % 各代路径的平均长度while iter <= iter_max  start = ones(m,1)+1;           %初始出发城市为2城市，即碑林区Table1(:,1) = start;Table2(:,1) = start;citys_index = 1:n;% 逐个蚂蚁路径选择for i = 1:m% 逐个城市路径选择for j = 2:ntabu1 = Table1(i,1:(j - 1));       %蚂蚁1已访问的城市集合tabu1(tabu1==0) = [];              %把表中的0元素舍去，因为无意义tabu2 = Table2(i,1:(j - 1));       %蚂蚁2已访问的城市集合tabu2(tabu2==0) = [];              %把表中的0元素舍去，因为无意义tabu = [tabu1 tabu2];              allow_index = ~ismember(citys_index,tabu);allow = citys_index(allow_index);       %待访问的城市集合P1 = allow;P2 = allow;% 计算城市间转移概率for k = 1:length(allow)P1(k) = Tau(tabu1(end),allow(k))^alpha * Eta(tabu1(end),allow(k))^beta;P2(k) = Tau(tabu2(end),allow(k))^alpha * Eta(tabu2(end),allow(k))^beta;endP1 = P1/sum(P1);P2 = P2/sum(P2);% 轮盘赌法选择下一个访问城市Pc1 = cumsum(P1);Pc2 = cumsum(P2);target_index1 = find(Pc1 >= rand);target_index2 = find(Pc2 >= rand);target1 = allow(target_index1(1));target2 = allow(target_index2(1));% 计算蚂蚁1与蚂蚁2到下一个城市的路径长度，选择路径长度较少的蚂蚁进行行走% 假设蚂蚁1行走，则蚂蚁1的路径表将添加下一个城市的标号，而蚂蚁2的路径表将添加0lengthant1 = D(tabu1(end), target1);lengthant2 = D(tabu2(end), target2);if( lengthant1 <= lengthant2)Table1(i,j) = target1;elseTable2(i,j) = target2;endendend% 计算各个蚂蚁对的路径距离Length1 = zeros(m,1);    % 蚂蚁1路径距离之和Length2 = zeros(m,1);    % 蚂蚁2路径距离之和for i = 1:mRoute1 = Table1(i,:);Route1(Route1 == 0) = [];    %舍去表中0元素Route2 = Table2(i,:);Route2(Route2 == 0) = [];    %舍去表中0元素for j = 1:length(Route1)-1Length1(i) = Length1(i) + D(Route1(j),Route1(j + 1));endLength1(i) = Length1(i) + D(Route1(end),Route1(1));for j = 1:length(Route2)-1Length2(i) = Length2(i) + D(Route2(j),Route2(j + 1));endLength2(i) = Length2(i) + D(Route2(end),Route2(1));endLength = Length1 +  Length2;% 计算最短路径距离及平均距离if iter == 1[min_Length,min_index] = min(Length);Length_best(iter) = min_Length;Length_ave(iter) = mean(Length);       Route_best1(iter,:) =  Table1(min_index,:);     Route_best2(iter,:) =  Table2(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_best1(iter,:) =  Table1(min_index,:);        Route_best2(iter,:) =  Table2(min_index,:);elseRoute_best1(iter,:) = Route_best1((iter-1),:);Route_best2(iter,:) = Route_best2((iter-1),:);endend% 更新信息素Delta_Tau = zeros(n,n);% 逐个蚂蚁计算for i = 1:mRoute1 = Table1(i,:);Route1(Route1 == 0) = [];     %舍去表中0元素Route2 = Table2(i,:);Route2(Route2 == 0) = [];     %舍去表中0元素% 逐个城市计算for j = 1:length(Route1)-1Delta_Tau(Route1(j),Route1(j+1)) = Delta_Tau(Route1(j),Route1(j+1)) + Q/Length1(i);endfor j = 1:length(Route2)-1Delta_Tau(Route2(j),Route2(j+1)) = Delta_Tau(Route2(j),Route2(j+1)) + Q/Length2(i);endDelta_Tau(Route1(end),Route1(1)) = Delta_Tau(Route1(end),Route1(1)) + Q/Length1(i);Delta_Tau(Route2(end),Route1(1)) = Delta_Tau(Route2(end),Route1(1)) + Q/Length2(i);endTau = (1-rho) * Tau + Delta_Tau;% 迭代次数加1，清空路径记录表iter = iter + 1;Table1 = zeros(m,n);Table2 = zeros(m,n);
end[Shortest_Length,index] = min(Length_best);
Shortest_Route1 = Route_best1(index,:);
Shortest_Route1(Shortest_Route1 == 0) = [];
Shortest_Route2 = Route_best2(index,:);
Shortest_Route2(Shortest_Route2 == 0) = [];
disp(['最短距离: ' num2str(Shortest_Length)]);
disp(['最短路径1: ' num2str([Shortest_Route1 Shortest_Route1(1)])]);
disp(['最短路径2: ' num2str([Shortest_Route2 Shortest_Route2(1)])]);
figure(2)
plot(1:iter_max,Length_best,'b',1:iter_max,Length_ave,'r:')
legend('最短距离','平均距离')
xlabel('迭代次数')
ylabel('距离')
title('各代最短距离与平均距离对比')
% r = TSP;
% sum = 0;
% for  i = 1:107
%     sum = sum+D(r(i),r(i+1));
% end
% sum

实验结果：

最短距离: 7192.7
最短路径1: 2    6    9   12   35   34   36   37   39   40   51    7   16   14   15   17   66   63   62   61   54   58   57   59   60   83   82   81   78   79   80   89   84   85   86   87   88   56   55   64   65   50   49   48   47   52   46   53   45   44   43   10  100  102  103  101  104  105  106   93   92   91   90   98   99   96   97   95   94   75   71   69   70   68   67   72   73   74   76   28   77    2
最短路径2: 2    3    5    1    4    8   13   32   33   31   30  107   42   41   11   23   22   21   20   18   19   26   25   29   24   27   38    2
>>

陕西公路里程表数据
提取码：1234
才疏学浅，欢迎大家指正~