变领域搜索算法求解旅行商问题

《MATLAB智能优化算法》学习笔记二

一、TSP概述

旅行商问题，即TSP问题（Traveling Salesman Problem）又译为旅行推销员问题、货郎担问题，是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值（等价于求图的最短哈密尔顿回路问题）。

举个栗子，下表是四个城市的位置信息：

因为一共需要访问4个城市，所以全部的行走路线方案的数目为$A_4^4=24$。因为数据规模较小，所以把所有的行走路线方案全部列出来。如下表所示，其中黄色底色的方案表示行走距离最短的行走路线方案。

进一步观察发现，方案1、10、17、和19其实是一种方案，方案6、8、15和24其实也是同一种方案。虽然看似方案的起点和终点不同，但是旅行商人的行走路线实际上是一条闭环回路。

因为一条行走路线的起点和终点是同一个城市，所以在表示该路线时为了表现得更加简洁，将终点城市从行走路线中省略。以12341这条行走路线为例，在省略终点城市1后，该条行走路线变为1234。后面部分所有关于路线的表现形式全部省略终点城市。

二、算法简介

就TSP而言，最终的目的是为旅行商人找到条最短路线。因为一开始最短路线是未知的，所以只能先尝试随便为旅行商人设计出一条一个城市只能访问一次的行走路线。假设这条路线不是最短的那条路线，那么这条路线一定还可以通过调整旅行商人访问城市的顺序，从而达到缩短行走总距离的目的。因此，调整旅行商人访问城市顺序的方法就显得尤为重要，不同的调整方法对于缩短行走总距离的效果是不同的。

以上节中1243这条路线为例，由表可知，这条路线的行走总距离为41.4km。因为在上节已知最短总距离为33.3km，所以1243这条路线一定还有调整的空间。那究竟如何调整1243路线呢？二个容易想到的策略是对调相邻城市的访问顺序。因此，采用该策略对1243路线调整后的路线分别如下。

调整路线1：2143（将1和2对调）．此时总距离为33.3km。

调整路线2：1423（将2和4对调），此时总距离为33.9km。

调整路线3：1234（将4和3对调）．此时总距离为33.3km。

调整路线4：3241（将3和1对调），此时总距离为33.9km。

这里需要注意的一点是，3和1其实也是相邻城市，因此也可以进行对调。接下来从上述4条路线中选择最短的一条路线，因为调整路线1和调整路线3总距离是相等的，这里在经过随机选择后，假设选择的是调整路线3。

因为此时不知道1234这条路线是否为最短路线，所以还会继续对这条路线进行调整。此时用一种新的策略调整1234这条路线，即对调中间间隔一个城市的两个城市的访问顺序。采用该策略对1234路线调整后的路线分别如下。

调整路线5：3214（将1和3对调），此时总距离为33.3km。

调整路线6：1432（将2和4对调）．此时总距离为33.3km。

虽然这2条路线总距离都为33.3km，但是此时依然不确定是否找到最短路线。因此，还会从调整路线5和调整路线6这2条路线中随机选择一条路线，然后继续对该路线进行调整，一直达到初始设定的条件为止，才会停止继续调整与寻找更短的路线。假设最初设定只允许调整路线3次，那么在调整3次路线后，就会自动停止寻找更短路线，此时所找到的最短路线就记为最短路线。

综上所述，对调相邻城市的访问顺序的策略称为一个规则，记为规则1；对调中间间隔一个城市的两个城市的访问顺序也是一个规则，记为规则。第1、2、3、4条调整路线实际上表示一个集合，记为集合1；第5、6条调整路线实际上也表示一个集合，记为集合2只有按照规则对原始路线进行操作后才能形成集合，形成集合的目的是从集合中找出比原始路线更短的路线。

至此，引出“变邻域搜索”的概念。规则对应邻域动作，集合对应邻域，不同的集合体现出变邻域中的变．从集合中找出比原始路线更短的路线对应搜索。变邻域搜索就是在不同的领域中不断地搜索更优的解。 VNS 算法求解组合优化问题（以求最小化问题为例）的常规流程如图所示。

三、算法设计

（1）构造初始路线

就 TSP 而言，一条好的初始路线能够节省VNS的求解时间。在这采用常规且应用较为广泛的构造 TSP 初始路线的方法一贪婪算法。

假设 TSP 中城市数目为 N ，贪婪算法构造 TSP 初始路线的步骤如下。

STEP1：初始化已被访问的城市集合 visited 为空，初始化未被访问的城市集合 unvisited ={1,2…, N}。

STEP2：将 N 行 N 列距离矩阵 dist 的主对角线上的0全部赋值为无穷大。

STEP3：从距离矩阵 dist 中找出最小距离对应的行序号 row 和列序号 col ,如果存在多个最小距离，则选择行序号 row 中的第一个数 row (1）作为初始路线的起点 first 。

STEP4：更新 visited =[ visited , first ]，即将 first 添加到 visited 中，同时也更新 unvisited ( unvisited = first ) = [ ] 即将 first 从 unvisited 中删除，然后将起点 first 赋值给紧前点 pre_point 。

STEP5：首先在距离矩阵 dist 中找到紧前点 pre_point 对应的那一行距离 pre_dist ，即紧前点 pre_point 与其他城市之间的距离；其次将已被访问的城市排除在外，即在 pre_dist 中将 visited 对应的列全设为无穷大；然后找出 pre_dist 中最小值对应的列序号作为下一个紧前点 pre_point 。

STEP6：更新 visited =[ visited , pre_point ], unvisited ( unvisited = pre_point ) = [ ]。

STEP7：若 unvisited 非空，则转至STEP5，否则转至STEP8。

STEP8：将 visited 赋值给初始路线 init_route ，初始路线构造完毕。

（2）交换操作

交换前路线： 1 2 3 4 5 6

若选择的交换位置$i=2$和$j=5$，那么交换第二位和第五位的两座城市

交换后路线： 1 5 3 4 2 6

（3）逆转操作

交换前路线： 1 2 3 4 5 6

若选择的逆转位置$i=2$和$j=5$，那么将位置二和五之间的城市顺序逆转

交换后路线： 1 5 3 4 2 6

（4）插入操作

例1：

交换前路线： 1 2 3 4 5 6

若选择的插入位置$i=2$和$j=5$，那么将第二位的城市插到第五位城市的后面

交换后路线： 1 3 4 5 2 6

例2：

交换前路线： 1 2 3 4 5 6

若选择的插入位置$i=5$和$j=2$，那么将第五位的城市插到第二位城市的后面

交换后路线： 1 2 5 3 4 6

（5）扰动操作

当使用某个邻域操作（交换、逆转和插入）前就先对当前路线使用该邻域操作以获得一个“扰动解”，然后再在这个“扰动解”上进行后续一些列操作。

（6）邻域搜索策略

在对一个当前解$S_{curr}$使用某个邻域操作得到$S_{curr}$的邻域集合后，如何对所得到的邻域集合进行处理才能使$S_{curr}$向“更好”的方向变换？

既然已经得到$S_{curr}$的邻域集合，那么就可以首先求出这个邻域集合中所有解的总距离；其次找到总距离最短的那个解$S_{min}$，并替换当前解$S_{curr}$，即令$S_{curr}=S_{min}$；然后求出$S_{curr}$的邻域集合，以及这个邻域集合中所有解的总距离，同样再找到总距离最短的那个解$S_{min}$，并替换当前解$S_{curr}$。

一直按照上述方式迭代，直至迭代 M 次后，停止对$S_{curr}$在当前邻域的搜索。

（7）邻域变化策略

邻域搜索策略详细介绍了如何针对一个当前解$S_{curr}$在一个邻域中进行搜索，因为本节介绍了3个邻域动作，所以会有3个不同结构的邻域。现在将这3个邻域分别编号为 k = 1、k = 2和 k = 3，即第1个邻域为交换操作得到的邻域，第2个邻域为逆转操作得到的邻域，第3个邻域为插人操作得到的邻域。

那么在对一个邻域搜索结束后，如何能够在下一个邻域中继续搜索呢？假设此时当前解为$S_{curr}$, $S_{curr}$的总距离为为$L_{curr}$，令最优路线$S_{best}=S_{curr}$，最优路线总距离$L_{best}=L_{curr}$

具体步骤如下。

STEP1：设 k = 1。

STEP2：如果 k = 1，转至STEP3；如果 k = 2，转至STEP4；如果 k = 3，转至STEP5；否则，转至STEP7。

STEP3：对$S_{curr}$进行扰动操作得到$S_{swap}$ ,$S_{swap}$的总距离为$L_{swap}$，令$L_{curr}=L_{swap}$。如果 $L_{curr}<L_{best}$ 则 $S_{best}=S_{curr}=S_{swap}$, $L_{best}=L_{curr}$ , k = 0，转至STEP6。

STEP4：对$S_{curr}$进行扰动操作得到$S_{reversion}$, $S_{reversion}$的总距离为$L_{reversion}$，令$L_{curr}=L_{reversion}$ 如如果 $L_{reversion}$ < $L_{best}$ ，则$S_{best}=S_{curr}=S_{reversion}$ , $L_{best}=L_{curr}$ , k =0，转至STEP6。

STEP5：对$S_{curr}$进行扰动操作得到$S_{insertion}$ , $S_{insertion}$ 的总距离为$L_{insertion}$ ，令$L_{curr}=L_{insertion}$ 。如果$L_{curr}<L_{best}$ ，则$S_{best}=S_{curr}=S_{insertion}$ , $L_{best}=L_{curr}$ , k =0，转至STEP6。

STEP6：k = k +1，转至STEP2。

STEP7：终止循环，输出$S_{curr}、S_{best}、L_{curr}$ 和$L_{best}$

（8）VNS求解TSP流程

四、MATLAB程序实现

（1）构造初始路线函数：

%% 贪婪算法构造TSP的初始解
%输入dist：             距离矩阵   
%输出init_route：      贪婪算法构造的初始路线
%输出init_len：        init_route的总距离  
function [init_route,init_len]=construct_route(dist)
N=size(dist,1);                     %城市数目
%先将距离矩阵主对角线上的0赋值为无穷大
for i=1:N
    for j=1:N
        if i==j
            dist(i,j)=inf;
        end
    end
end

unvisited=1:N;                      %初始未被安排的城市集合
visited=[];                         %初始已被安排的城市集合

min_dist=min(min(dist));            %找出距离矩阵中的最小值
[row,col]=find(dist==min_dist);     %在dist中找出min_dist所对应的行和列
first=row(1);                       %将min_dist在dist中所对应行序号作为起点

unvisited(unvisited==first)=[];     %将first从unvisit中删除
visited=[visited,first];            %把first添加到visit中
pre_point=first;                    %将fisrt赋值给pre_point
while ~isempty(unvisited)
    pre_dist=dist(pre_point,:);     %pre_point与其它城市的距离
    pre_dist(visited)=inf;          %将pre_point与已经添加进来的城市之间的距离设位无穷大
    [~,pre_point]=min(pre_dist);    %找出pre_dist中的最小值
    unvisited(unvisited==pre_point)=[]; %将pre_point从unvisit中删除
    visited=[visited,pre_point];        %把pre_point添加到visit中
end
init_route=visited;
init_len=route_length(init_route,dist); %计算init_route的总距离
end

（2）路线总距离计算函数：

%% 计算一条路线总距离
%输入route：            一条路线
%输入dist：             距离矩阵
%输出len：              该条路线总距离
function len=route_length(route,dist)
    n=numel(route);
    route=[route route(1)];
    len=0;
    for k=1:n 
        i=route(k);
        j=route(k+1); 
        len=len+dist(i,j); 
    end
end

（3）交换操作函数：

%% 交换操作
%比如说有6个城市，当前解为123456，随机选择两个位置，然后将这两个位置上的元素进行交换。
%比如说，交换2和5两个位置上的元素，则交换后的解为153426。
%输入route1：          路线1
%输入i,j：             两个交换点
%输出route2：          经过交换操作变换后的路线2
function route2=swap(route1,i,j)
route2=route1;
route2([i j])=route1([j i]);
end

（4）逆转操作函数：

%% 逆转操作
%有6个城市，当前解为123456，我们随机选择两个位置，然后将这两个位置之间的元素进行逆序排列。
%比如说，逆转2和5之间的所有元素，则逆转后的解为154326。
%输入route1：          路线1
%输入i,j：             逆转点i,j
%输出route2：          经过逆转结构变换后的路线2
function route2=reversion(route1,i,j)
i1=min([i,j]);
i2=max([i,j]);
route2=route1;
route2(i1:i2)=route1(i2:-1:i1);
end

（5）插入操作函数：

%% 插入操作
%有6个城市，当前解为123456，我们随机选择两个位置，然后将这第一个位置上的元素插入到第二个元素后面。
%比如说，第一个选择5这个位置，第二个选择2这个位置，则插入后的解为125346。
%输入route1：          路线1
%输入i,j：             插入点i,j
%输出route2：          经过插入结构变换后的路线2
function route2=insertion(route1,i,j)
if i<j
    route2=route1([1:i-1 i+1:j i j+1:end]);
else
    route2=route1([1:j i j+1:i-1 i+1:end]);
end
end

（6）计算距离差值函数：

%% 计算swap操作后与操作前路线route的总距离的差值
%输入route：            一条路线
%输入dist：             距离矩阵
%输入i，j：             交换点i,j
%输出delta1：           swap后路线的总距离-swap前路线的总距离
function delta1=cal_delta1(route,dist,i,j)
N=numel(route);         %城市数目
if (i==1)&&(j==N)
    delta1=-(dist(route(i),route(i+1))+dist(route(j-1),route(j)))+...
        (dist(route(j),route(i+1))+dist(route(j-1),route(i)));
elseif (i==1)&&(j==2)
    delta1=-(dist(route(N),route(i))+dist(route(j),route(j+1)))+...
        (dist(route(N),route(j))+dist(route(i),route(j+1)));
elseif i==1
    delta1=-(dist(route(N),route(i))+dist(route(i),route(i+1))+...
        dist(route(j-1),route(j))+dist(route(j),route(j+1)))+...
        (dist(route(N),route(j))+dist(route(j),route(i+1))+...
        dist(route(j-1),route(i))+dist(route(i),route(j+1)));
elseif (i==N-1)&&(j==N)
    delta1=-(dist(route(i-1),route(i))+dist(route(j),route(1)))+...
        (dist(route(i-1),route(j))+dist(route(i),route(1)));
elseif j==N
    delta1=-(dist(route(i-1),route(i))+dist(route(i),route(i+1))+...
        dist(route(j-1),route(j))+dist(route(j),route(1)))+...
        (dist(route(i-1),route(j))+dist(route(j),route(i+1))+...
        dist(route(j-1),route(i))+dist(route(i),route(1)));
elseif abs(i-j)==1
    delta1=-(dist(route(i-1),route(i))+dist(route(j),route(j+1)))+...
        (dist(route(i-1),route(j))+dist(route(i),route(j+1)));
else
    delta1=-(dist(route(i-1),route(i))+dist(route(i),route(i+1))+...
        dist(route(j-1),route(j))+dist(route(j),route(j+1)))+...
        (dist(route(i-1),route(j))+dist(route(j),route(i+1))+...
        dist(route(j-1),route(i))+dist(route(i),route(j+1)));
end
end

%% 将给定的route序列在i和j位置之间进行逆序排列，然后计算转换序列前和转换序列后的路径距离的差值
%输入route：            一条路线
%输入dist：             距离矩阵
%输入i，j：             逆转点i,j
%输出delta2：           reversion后路线的总距离-reversion前路线的总距离
function delta2=cal_delta2(route,dist,i,j)
N=numel(route);                 %城市个数
if i==1
    if j==N
        delta2=0;
    else
        delta2=-dist(route(j),route(j+1))-dist(route(N),route(i))+...
            dist(route(i),route(j+1))+dist(route(N),route(j));
    end
else
    if j==N
        delta2=-dist(route(i-1),route(i))-dist(route(1),route(j))+...
            dist(route(i-1),route(j))+dist(route(i),route(1));
    else
        delta2=-dist(route(i-1),route(i))-dist(route(j),route(j+1))+...
            dist(route(i-1),route(j))+dist(route(i),route(j+1));
    end
end
end

%% 计算insertion操作后与操作前路线route的总距离的差值
%输入route：            一条路线
%输入dist：             距离矩阵
%输入i，j：             逆转点i,j
%输出delta1：           insertion后路线的总距离-insertion前路线的总距离
function delta3=cal_delta3(route,dist,i,j)
N=numel(route);         %城市数目
if i<j
    if (i==1) && (j==N)
        delta3=0;
    elseif (i==1) && (j==2)
        delta3=-(dist(route(N),route(i))+dist(route(j),route(j+1)))+...
            (dist(route(N),route(j))+dist(route(i),route(j+1)));
    elseif i==1
        delta3=-(dist(route(N),route(i))+dist(route(i),route(i+1))+dist(route(j),route(j+1)))+...
            (dist(route(N),route(i+1))+dist(route(j),route(i))+dist(route(i),route(j+1)));
    elseif (i==N-1)&&(j==N)
        delta3=-(dist(route(i-1),route(i))+dist(route(j),route(1)))+...
            (dist(route(i-1),route(j))+dist(route(i),route(1)));
    elseif j==N
        delta3=-(dist(route(i-1),route(i))+dist(route(i),route(i+1))+dist(route(j),route(1)))+...
            (dist(route(i-1),route(i+1))+dist(route(j),route(i))+dist(route(i),route(1)));
    elseif (j-i)==1
        delta3=-(dist(route(i-1),route(i))+dist(route(j),route(j+1)))+...
            (dist(route(i-1),route(j))+dist(route(i),route(j+1)));
    else
        delta3=-(dist(route(i-1),route(i))+dist(route(i),route(i+1))+dist(route(j),route(j+1)))+...
            (dist(route(i-1),route(i+1))+dist(route(j),route(i))+dist(route(i),route(j+1)));
    end
else
    if (i==N) && (j==1)
        delta3=-(dist(route(i-1),route(i))+dist(route(j),route(j+1)))+...
            (dist(route(i-1),route(j))+dist(route(i),route(j+1)));
    elseif (i-j)==1
        delta3=0;
    elseif i==N
        delta3=-(dist(route(i-1),route(i))+dist(route(i),route(1))+dist(route(j),route(j+1)))+...
            (dist(route(i-1),route(1))+dist(route(j),route(i))+dist(route(i),route(j+1)));
    else
        delta3=-(dist(route(i-1),route(i))+dist(route(i),route(i+1))+dist(route(j),route(j+1)))+...
            (dist(route(i-1),route(i+1))+dist(route(j),route(i))+dist(route(i),route(j+1)));
    end
end
end

（7）更新距离差值函数：

%% swap操作后生成新的距离差矩阵Delta
%输入route：            一条路线
%输入dist：             距离矩阵
%输入i，j：             交换点i,j
%输出Delta1：           swap操作后的距离差值的矩阵
function Delta1=Update1(route,dist,i,j)
N=numel(route);                     %城市个数
route2=swap(route,i,j);             %交换route上i和j两个位置上的城市
Delta1=zeros(N,N);                  %N行N列的Delta初始化，每个位置上的元素是距离差值
for i=1:N-1
    for j=i+1:N
        Delta1(i,j)=cal_delta1(route2,dist,i,j);
    end
end

%% reversion操作后生成新的距离差矩阵Delta
%输入route：            一条路线
%输入dist：             距离矩阵
%输入i，j：             逆转点i,j
%输出Delta2：           reversion操作后的距离差值的矩阵
function Delta2=Update2(route,dist,i,j)
N=numel(route);                     %城市个数
route2=reversion(route,i,j);        %逆转route上i和j两个位置上的城市
Delta2=zeros(N,N);                  %N行N列的Delta初始化，每个位置上的元素是距离差值
for i=1:N-1
    for j=i+1:N
        Delta2(i,j)=cal_delta2(route2,dist,i,j);
    end
end

%% insertion操作后生成新的距离差矩阵Delta
%输入route：            一条路线
%输入dist：             距离矩阵
%输入i，j：             插入点i,j
%输出Delta1：           insertion操作后的距离差值的矩阵
function Delta3=Update3(route,dist,i,j)
N=numel(route);                     %城市个数
route2=insertion(route,i,j);        %插入route上i和j两个位置上的城市
Delta3=zeros(N,N);                  %N行N列的Delta初始化，每个位置上的元素是距离差值
for i=1:N
    for j=1:N
        if i~=j
            Delta3(i,j)=cal_delta3(route2,dist,i,j);
        end
    end
end

（8）扰动操作函数：

%% 扰动，随机选择当前邻域中的一个解更新当前解
%输入route：           一条路线
%输入dist：            距离矩阵
%输入k：               当前邻域序号
%输出route_shake：     扰动操作后得到的路线
%输出len_shake：       该条路线的距离
function [route_shake,len_shake]=shaking(route,dist,k)
N=numel(route);                                                     %城市数目
select_no=randi([1,N],1,2);                                         %随机选择进行操作的两个点的序号
i=select_no(1);
j=select_no(2);
if k==1
    route_shake=swap(route,i,j);
elseif k==2
    route_shake=reversion(route,i,j);
else
    route_shake=insertion(route,i,j);
end
len_shake=route_length(route_shake,dist);
end

（9）交换邻域搜索函数：

%% 对route不断进行交换操作后所得到的路线以及所对应的总距离
%输入route：           一条路线
%输入dist：            距离矩阵
%输入M：               最多进行邻域操作的次数
%输出swap_route：      对route不断进行交换操作后所得到的路线
%输出swap_len：        swap_route的总距离
function [swap_route,swap_len]=swap_neighbor(route,dist,M)
N=numel(route);         %城市数目
Delta1=zeros(N,N);      %交换任意两个位置之间序列的元素所产的距离差的矩阵
for i=1:N-1
    for j=i+1:N
        Delta1(i,j)=cal_delta1(route,dist,i,j);
    end
end
cur_route=route;                                                    %初始化当前路线
m=1;                                                                %初始化计数器
while m<=M
    min_value=min(min(Delta1));                                     %找出距离差值矩阵中最小的距离差值
    %如果min_value小于0，才能更新当前路线和距离矩阵。否则，终止循环
    if min_value<0
        [min_row,min_col]=find(Delta1==min_value);                  %找出距离差值矩阵中最小的距离差值所对应的行和列
        Delta1=Update1(cur_route,dist,min_row(1),min_col(1));       %更新距离差值矩阵
        cur_route=swap(cur_route,min_row(1),min_col(1));            %更新当前路线
    else
        break
    end
    m=m+1;
end
swap_route=cur_route;                                               %将当前路线cur_route赋值给swap_route
swap_len=route_length(swap_route,dist);                             %swap_route的总距离
end

（10）逆转邻域搜索函数：

%% 对route不断进行逆转操作后所得到的路线以及所对应的总距离
%输入route：           一条路线
%输入dist：            距离矩阵
%输入M：               最多进行M次邻域操作
%输出reversion_route： 对route不断进行逆转操作后所得到的路线
%输出reversion_len：   reversion_route的总距离
function [reversion_route,reversion_len]=reversion_neighbor(route,dist,M)
N=numel(route);         %城市数目
Delta2=zeros(N,N);      %逆转任意两个位置之间序列的元素所产的距离差的矩阵
for i=1:N-1
    for j=i+1:N
        Delta2(i,j)=cal_delta2(route,dist,i,j);
    end
end
cur_route=route;                                                    %初始化当前路线
m=1;                                                                %初始化计数器
while m<=M
    min_value=min(min(Delta2));                                     %找出距离差值矩阵中最小的距离差值
    %如果min_value小于0，才能更新当前路线和距离矩阵。否则，终止循环
    if min_value<0
        [min_row,min_col]=find(Delta2==min_value);                  %找出距离差值矩阵中最小的距离差值所对应的行和列
        Delta2=Update2(cur_route,dist,min_row(1),min_col(1));       %更新距离差值矩阵
        cur_route=reversion(cur_route,min_row(1),min_col(1));       %更新当前路线
    else
        break
    end
    m=m+1;
end
reversion_route=cur_route;                                          %将当前路线cur_route赋值给reversion_route
reversion_len=route_length(reversion_route,dist);                   %reversion_route的总距离
end

（11）插入邻域搜索函数：

%% 对route不断进行插入操作后所得到的路线以及所对应的总距离
%输入route：           一条路线
%输入dist：            距离矩阵
%输入M：               最多进行M次邻域操作
%输出insertion_route： 对route不断进行插入操作后所得到的路线
%输出insertion_len：   insertion_route的总距离
function [insertion_route,insertion_len]=insertion_neighbor(route,dist,M)
N=numel(route);         %城市数目
Delta3=zeros(N,N);      %插入任意两个位置之间序列的元素所产的距离差的矩阵
for i=1:N-1
    for j=i+1:N
        Delta3(i,j)=cal_delta3(route,dist,i,j);
    end
end
cur_route=route;                                                    %初始化当前路线
m=1;                                                                %初始化计数器
while m<=M
    min_value=min(min(Delta3));                                     %找出距离差值矩阵中最小的距离差值
    %如果min_value小于0，才能更新当前路线和距离矩阵。否则，终止循环
    if min_value<0
        [min_row,min_col]=find(Delta3==min_value);                  %找出距离差值矩阵中最小的距离差值所对应的行和列
        Delta3=Update3(cur_route,dist,min_row(1),min_col(1));       %更新距离差值矩阵
        cur_route=insertion(cur_route,min_row(1),min_col(1));       %更新当前路线
    else
        break
    end
    m=m+1;
end
insertion_route=cur_route;                                          %将当前路线cur_route赋值给insertion_route
insertion_len=route_length(insertion_route,dist);                   %insertion_route的总距离
end

（12）TSP路线图函数：

%% TSP路线可视化
%输入route：           一条路线
%输入x,y：             x,y坐标
function plot_route(route,x,y)
figure
route=[route route(1)];
plot(x(route),y(route),'k-o','MarkerSize',10,'MarkerFaceColor','w','LineWidth',1.5);
xlabel('x');
ylabel('y');
end

（13）主函数：

tic
clear
clc
%% 输入数据
dataset=importdata('input.txt');                %数据中，每一列的含义分别为[序号，x坐标，y坐标]
x=dataset(:,2);                                 %x坐标
y=dataset(:,3);                                 %y坐标
vertexes=dataset(:,2:3);                        %提取各个城市的xy坐标
N=size(dataset,1);                              %城市数目
h=pdist(vertexes);                              %计算各个城市之间的距离，一共有1+2+......+(n-1)=n*(n-1)/2个
dist=squareform(h);                             %将各个城市之间的距离转换为n行n列的距离矩阵
%% 参数初始化
MAXGEN=50;                                      %外层最大迭代次数
M=50;                                           %最多进行M次邻域操作
n=3;                                            %邻域数目
%% 构造初始解
[init_route,init_len]=construct_route(dist);    %贪婪构造初始解
disp(['初始路线总距离为',num2str(init_len)]);
cur_route=init_route;
best_route=cur_route;
best_len=route_length(cur_route,dist);
BestL=zeros(MAXGEN,1);                          %记录每次迭代过程中全局最优个体的总距离
%% 主循环
gen=1;                                          %外层计数器
while gen<=MAXGEN
    k=1;
    while(1)
        switch k
            case 1
                cur_route=shaking(cur_route,dist,k);
                [swap_route,swap_len]=swap_neighbor(cur_route,dist,M);
                cur_len=swap_len;
                if cur_len<best_len
                    cur_route=swap_route;
                    best_len=cur_len;
                    best_route=swap_route;
                    k=0;
                end
            case 2
                cur_route=shaking(cur_route,dist,k);
                [reversion_route,reversion_len]=reversion_neighbor(cur_route,dist,M);
                cur_len=reversion_len;
                if cur_len<best_len
                    cur_route=reversion_route;
                    best_len=cur_len;
                    best_route=reversion_route;
                    k=0;
                end
            case 3
                cur_route=shaking(cur_route,dist,k);
                [insertion_route,insertion_len]=insertion_neighbor(cur_route,dist,M);
                cur_len=insertion_len;
                if cur_len<best_len
                    cur_route=insertion_route;
                    best_len=cur_len;
                    best_route=insertion_route;
                    k=0;
                end
            otherwise
                break;
        end
        k=k+1;
    end
    disp(['第',num2str(gen),'代最优路线总距离为',num2str(best_len)]);
    BestL(gen,1)=best_len;
    %% 计数器加1
    gen=gen+1;
end
%% 画出优化过程图
figure;
plot(BestL,'LineWidth',1);
title('优化过程')
xlabel('迭代次数');
ylabel('总距离');
%% 画出全局最优路线图
plot_route(best_route,x,y);
toc

五、实验结果展示

¹. 《MATLAB智能优化算法》（曹旺著） ↩