利用Matlab以蚁群算法（Ant Colony Algorithm）求解不闭合的旅行商问题（Travelling Salesman Problem, TSP）并进行可视化

旅行商问题是一种十分重要的问题，在生产生活中都有重要应用。其变式也很多。这里讨论一下不闭合的旅行商问题的Matlab表示，并用相应函数来表示其产生的路径。

不闭合的旅行商问题，也就是起点和终点不需要重合，寻找连接所有点的最短路径的问题。

做理论研究时可能就要另起一个话题了，但是对于编程求解来说，因为实现环路是把原始的城市编号集（例如50个城市，那么现在编号集就有50个数据）直接在尾部添加编号集的第一个数据，这样就成为51个数据了，所以只要去掉这一步，就成为不闭合的TSP了。这是对原理的解释，不同人写的程序的具体实现不同。有的是采用编号集增加数据，有的是采用计算路径长度时算上首尾的距离等等。

 

以下是对百度知道里面找到的已有算法（点击打开链接）进行的整理和补充优化（原算法用于求解闭合的，这里修改为不闭合的）：

数据样本：

C=[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  
];

处理：

clc
clear all
close all

C=[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  
];
m=31;Alpha=1;Beta=5;Rho=0.1;NC_max=200;Q=100;
[R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha,Beta,Rho,Q)
function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha,Beta,Rho,Q)
%% 主要符号说明
%C n个城市的坐标，n×2的矩阵
%NC_max 最大迭代次数
%m 蚂蚁个数
% Alpha 表征信息素重要程度的参数
% Beta 表征启发式因子重要程度的参数
% Rho 信息素蒸发系数
% Q 信息素增加强度系数
% R_best 各代最佳路线
% L_best 各代最佳路线的长度
%% 第一步：变量初始化
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，但后面的启发因子要取倒数，用eps（浮点相对精度）表示
        end
        D(j,i)=D(i,j);   %对称矩阵
    end
end
Eta=1./D;          %Eta为启发因子，这里设为距离的倒数
Tau=ones(n,n);     %Tau为信息素矩阵
Tabu=zeros(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)];
    end
    Tabu(:,1)=(Randpos(1,1:m))';    %此句不太理解？
    %% 第三步：m只蚂蚁按概率函数选择下一座城市，完成各自的周游
    for j=2:n     %所在城市不计算
        for i=1:m     
            visited=Tabu(i,1:(j-1)); %记录已访问的城市，避免重复访问
            J=zeros(1,(n-j+1));       %待访问的城市
            P=J;                      %待访问城市的选择概率分布
            Jc=1;
            for k=1:n
                if length(find(visited==k))==0   %开始时置0
                    J(Jc)=k;
                    Jc=Jc+1;                         %访问的城市个数自加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);     %cumsum，元素累加即求和
            Select=find(Pcum>=rand); %若计算的概率大于原来的就选择这条路线
            to_visit=J(Select(1));
            Tabu(i,j)=to_visit;
        end
    end
    if NC>=2
        Tabu(1,:)=R_best(NC-1,:);
    end
    %% 第四步：记录本次迭代最佳路线
    L=zeros(m,1);     %开始距离为0，m*1的列向量
    for i=1:m
        R=Tabu(i,:);
        for j=1:(n-1)
            L(i)=L(i)+D(R(j),R(j+1));    %原距离加上第j个城市到第j+1个城市的距离
        end
        %L(i)=L(i)+D(R(1),R(n));      %一轮下来后走过的距离
    end
    L_best(NC)=min(L);           %最佳距离取最小
    pos=find(L==L_best(NC));
    R_best(NC,:)=Tabu(pos(1),:); %此轮迭代后的最佳路线
    L_ave(NC)=mean(L);           %此轮迭代后的平均距离
    NC=NC+1                      %迭代继续
    %% 第五步：更新信息素
    Delta_Tau=zeros(n,n);        %开始时信息素为n*n的0矩阵
    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; %考虑信息素挥发，更新后的信息素
    %% 第六步：禁忌表清零
    Tabu=zeros(m,n);             %%直到最大迭代次数
end
%% 第七步：输出结果
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('Average Distance and Shortest Distance')     %标题
function DrawRoute(C,R)
%% 画路线图的子函数
% C Coordinate 节点坐标，由一个N×2的矩阵存储
% R Route 路线
N=length(R);
scatter(C(:,1),C(:,2));
hold on
%plot([C(R(1),1),C(R(N),1)],[C(R(1),2),C(R(N),2)],'g')
%hold on
for ii=2:N
    plot([C(R(ii-1),1),C(R(ii),1)],[C(R(ii-1),2),C(R(ii),2)],'g')
    hold on
end
title('TSP Result')

结果：

代码解释：

函数ACATSP.m的第74行注释掉了，这样在计算路径长度时就不计首与尾的连接了，也就是成为了不闭合的TSP问题。

函数DrawRoute.m的第8和9行注释掉了，所以作出的图也就是不闭合的了。

这两个关键点保证计算和图示都是不闭合的TSP问题。