大规模邻域搜索算法求解旅行商问题

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

一、TSP概述

见上个学习笔记。

二、算法简介

虽然VNS在搜索过程中使用若干个不同的邻域以扩大搜索范围，但是当TSP中城市数目增多时，VNS搜索范围有限这一缺点就会暴露出来。为了能够进一步扩大搜索范围，本文使用大规模邻域搜索算法（LNS）求解TSP。LNS的思想是先“破坏”解，然后将破坏后的解进行“修复”，最终获得更高质量的解。

还是以上个学习笔记为例，假设TSP中城市数目为4，初始路线为1234，现在分别使用交换操作，逆转操作和插入操作求出初始解1234对应的邻域，可以得到10种不同邻域解。因此，除去初始解1234外还有13种可能解没有通过上述3种操作获得。由此可见，VNS在使用这3种邻域操作时搜索范围的局限性。

（1）在LNS求解问题过程中，如何“破坏”解，以及如何“修复”解？

“破坏”解其实就是将若干城市从当前路线中移除。“修复”解就是将被移除的城市再重新插入到“破坏”解中去。

比如当前路线为1234，我要移除城市2、3、4，此时“破坏”后的解就是1。

“修复”的过程可以看下图：

每次插入的时候都有不同的位置可以选择。

（2）LNS的大规模是如何体现的呢？

上述“修复”解的顺序为先将3插回，然后将4插回，最后将2插回，即插回顺序为342。实际上也可以按照234、243、324、423和432这5种插回顺序“修复”解。每种插回顺序同样也能得到这24个不同的“修复解”。因为当 TSP 中城市数目为4时，一共有$A_4^4=24$种排序方式，所以采用任何一种插回顺序都能得到全部排序方式，这其实就体现了大规模的思想，即能搜索到更多不同类型的解。

综上所述，LNS求解TSP的流程如下图所示，其中$f(S)$表示解$S$的行走距离。

三、求解策略

（1）构造初始解

与VNS构造TSP初始解相同，即采用贪婪算法构造TSP的初始解。

（2）“破坏”解

假设TSP中城市数目为$N$，当前解为$route$，且要从当前解中最多移除相连接的$L_{max}$个城市，最少移除相连接的$L_{min}$个城市，则具体的移除步骤如下。

STEP1：从$L_{min}\sim L_{max}$中随机取一个整数作为要移除的城市数目。

STEP2：从当前解$route$中随机选择一个城市$visit$作为即将移除的相连接$L$个城市的参考城市。

STEP3：计算$visit$在$route$中的位置$findv$，以起始点$route(1)$为界限计算在$route$中$visit$左侧的城市数目$v_{LN}=findv-1$，同理也计算出在$route$中$visit$右侧的城市数目$v_{RN}=N-findv$。

STEP4：如果$v_{LN}\le v_{RN}$，则转至STEP5，否则转至STEP6。

STEP5：分以下3种情况计算$visit$右侧要移除城市的数目$n_R$和左侧要移除城市的数目$n_L$。

（1）如果$v_{LN}<L-1$且$v_{RN}<L-1$，则$n_R=L-1-v_{LN}+(0\sim (v_{RN}-L+1+v_{RN}))$之间的随机整数，$n_L=L-1-n_R$。

（2）如果$v_{LN}>L-1$且$v_{RN}>L-1$，则$n_R=(0\sim v_{LN})$之间的随机整数，$n_L=L-1-n_R$。

（3）如果$v_{LN}\le L-1$且$v_{RN}\ge L-1$，则$n_R=L-1-v_{LN}+(0\sim v_{LN})$之间的随机整数，$n_L=L-1-N_R$。

SREP6：分以下3种情况计算$visit$右侧要移除城市的数目$n_R$和左侧要移除城市的数目$n_L$。

（1）如果$v_{LN}<L-1$且$v_{RN}<L-1$，则$n_L=L-1-v_{RN}+(0\sim (v_{LN}-L+1+v_{RN}))$之间的随机整数，$n_R=L-1-n_L$。

（2）如果$v_{LN}>L-1$且$v_{RN}>L-1$，则$n_L=(0\sim v_{RN})$之间的随机整数，$n_R=L-1-n_L$。

（3）如果$v_{LN}\le L-1$且$v_{RN}\ge L-1$，则$n_L=L-1-v_{RN}+(0\sim v_{RN})$之间的随机整数，$n_R=L-1-N_L$。

STEP7：从$route$中提取被移除的城市集合$removed=route(findv-n_L:findv+n_R)$，并将removed中所有城市从$route$中删除，得到“破坏”后的解$S_{destroy}$。

（3）“修复”解

在讲解“修复”解的步骤前，先阐述“插入成本”这一概念。如果当前“破坏”后的解为$S_{destroy}$，那么在将removed中的一个城市插回到$S_{destroy}$中的某个插入位置后，此时“修复”后的路线总距离减去$S_{destroy}$的总距离即为该城市插入该位置的“插入成本”。

在介绍“插入成本”之后，我们再来了解一下“遗憾值”这一概念。如果$S_{destroy}$中城市的数目为$lr$，那么将在removed中的一个城市插回到$S_{destroy}$时共有$lr+1$个可能的插入位置，这$lr+1$个插入位置对应$lr+1$个“插入成本”。

接下来将这$lr+1$个“插入成本”从小到大进行排序，若排序后的“插入成本”为up_delta，则该城市插回到$S_{repair}$的“遗憾值”即为排序后排在第二位的“插入成本”减去排在第一位的“插入成本”，即up_delta(2)-up_delta(1)。

“修复”的步骤如下：

STEP1：初始为“修复后”的解$S_{repair}$赋值，即$S_{repair}=S_{destroy}$。

STEP2：如果removed非空，转至STEP3，否则转至STEP6。

STEP3：计算当前removed中的城市数目$nr$，计算将removed中各个城市插回到$S_{repair}$的“遗憾值”regret，即regret是$nr$行1列的矩阵。

STEP4：找出regret中最大“遗憾值”对应的序号max_index，确定出即将被插回的城市renisert_city=removed(max_index)，将renisert_city插回到$S_{repair}$中“插入成本”最小的位置。

STEP5：更新removed(max_index)=[ ]，转至STEP2。

STEP6：“修复”结束，返回“修复”后的解$S_{repair}$。

四、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）“破坏”函数

%% 破坏函数destroy根据从当前解中连续移除若干个城市
%输入route：                   当前解，一条路线
%输出Sdestroy：                移除removed中的城市后的route
%输出removed：                 被移除的城市集合
function [Sdestroy,removed]=destroy(route)
N=numel(route);                                 %当前解中城市数目
Lmin=1;                                         %一条路径中所允许移除最小的城市数目
Lmax=min(ceil(N/2),25);                         %一条路径中所允许移除最大的城市数目
visit=ceil(rand*N);                             %从当前解中随机选出要被移除的城市
L=Lmin+ceil((Lmax-Lmin)*rand);                  %计算在该条路径上移除的城市的数目

findv=find(route==visit,1,'first');             %找出visit在route中的位置
vLN=findv-1;                                    %visit左侧的城市个数
vRN=N-findv;                                    %visit右侧的城市个数
%如果vLN小
if vLN<=vRN
    if (vRN<L-1)&&(vLN<L-1)
        nR=L-1-vLN+round(rand*(vRN-L+1+vLN));
        nL=L-1-nR;                              %visit左侧要移除城市的数目
    elseif (vRN>L-1)&&(vLN>L-1)
        nR=round(rand*vLN);                     %visit右侧要移除城市的数目
        nL=L-1-nR;                              %visit左侧要移除城市的数目
    else
        nR=L-1-vLN+round(rand*vLN);
        nL=L-1-nR;                              %visit左侧要移除城市的数目
    end
else
    %如果vRN小
    if (vLN<L-1)&&(vRN<L-1)
        nL=L-1-vRN+round(rand*(vLN-L+1+vRN));
        nR=L-1-nL;                              %visit右侧要移除城市的数目
    elseif (vLN>L-1)&&(vRN>L-1)
        nL=round(rand*vRN);                     %visit左侧要移除城市的数目
        nR=L-1-nL;                              %visit右侧要移除城市的数目
    else
        nL=L-1-vRN+round(rand*vRN);
        nR=L-1-nL;                              %visit右侧要移除城市的数目
    end
end
removed=route(findv-nL:findv+nR);               %移除城市的集合，即包括visit在内的连续L个城市

Sdestroy=route;                                 %复制route
for i=1:L
    Sdestroy(Sdestroy==removed(i))=[];          %将removed中的所有城市从route中移除
end
end

（4）“插入成本”计算函数

%% 将visit插回到插入成本最小的位置后的路线，同时还计算出插入到各个插入位置的插入成本
%输入visit：               待插入的城市
%输入dist：                距离矩阵
%输入route：               被插入路径
%输出new_route：           将visit插入到route最小插入成本位置后的解
%输出up_delta：            将visit插入到route中各个插入位置后的插入成本从小到大排序后的结果
function [new_route,up_delta]=ins_route(visit,dist,route)
lr=numel(route);                            %当前路线城市数目
rc0=zeros(lr+1,lr+1);                       %记录插入城市后的路径
delta0=zeros(lr+1,1);                       %记录插入城市后的增量
for i=1:lr+1
    if i==lr+1
        rc=[route visit];
    elseif i==1
        rc=[visit route];
    else
        rc=[route(1:i-1) visit route(i:end)];
    end
    rc0(i,:)=rc;                                            %将合理路径存储到rc0，其中rc0与delta0对应
    dif=route_length(rc,dist)-route_length(route,dist);     %计算成本增量
    delta0(i,1)=dif;                                        %将成本增量存储到delta0
end
up_delta=sort(delta0);                                      %将插入成本从小到大排序
[~,ind]=min(delta0);                                        %计算最小插入成本所对应的序号
new_route=rc0(ind,:);                                        %最小插入成本所对应的插入后的路径
end

（5）“修复”函数

%% 修复函数repair依次将removed中的城市插回路径中
%先计算removed中各个城市插回当前解中所产生最小增量，然后再从上述各个最小增量的城市中
%找出一个(距离增量第2小-距离增量第1小)最大的城市插回，反复执行，直到全部插回
%输入removed：                 被移除城市的集合
%输入Sdestroy：                “破坏”后的解
%输入dist：                    距离矩阵
%输出Srepair：                 “修复”后的解
%输出repair_length：           Srepair的总距离
function [Srepair,repair_length]=repair(removed,Sdestroy,dist)
Srepair=Sdestroy;
%反复插回removed中的城市，直到全部城市插回
while ~isempty(removed)
    nr=numel(removed);                                  %移除集合中城市数目
    regret=zeros(nr,1);                                 %存储removed各城市插回最“好”插回路径后的遗憾值增量
    %逐个计算removed中的城市插回当前解中各路径的目标函数值增
    for i=1:nr
        visit=removed(i);                               %当前要插回的城市
        [~,up_delta]=ins_route(visit,dist,Srepair);     %将visit插回到插入成本最小的位置后的路线，同时还计算出插入到各个插入位置的插入成本
        de12=up_delta(2)-up_delta(1);                   %计算第2小成本增量与第1小成本增量差值
        regret(i)=de12;                                 %更新当前城市插回最“好”插回路径后的遗憾值
    end
    [~,max_index]=max(regret);                          %找出遗憾值最大的城市序号
    reinsert_city=removed(max_index);                   %removed中准备插回的城市
    Srepair=ins_route(reinsert_city,dist,Srepair);      %将reinsert_city插回到Srepair
    removed(max_index)=[];                              %将removed(firIns)城市从removed中移除
end
repair_length=route_length(Srepair,dist);               %计算Srepair的总距离
end

（6）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

（7）主函数

tic
clear
clc
%% 输入数据
dataset=importdata('input.txt');                %数据中，每一列的含义分别为[序号，x坐标，y坐标]
x=dataset(:,2);                                 %x坐标
y=dataset(:,3);                                 %y坐标
vertexes=dataset(:,2:3);                        %提取各个城市的xy坐标
h=pdist(vertexes);
dist=squareform(h);                             %距离矩阵
%% 参数初始化
MAXGEN=300;                                     %最大迭代次数
%% 构造初始解
[Sinit,init_len]=construct_route(dist);         %贪婪构造初始解
init_length=route_length(Sinit,dist);
str1=['初始总路线长度 =  ' num2str(init_length)];
disp(str1)
%% 初始化当前解和全局最优解
Scurr=Sinit;
curr_length=init_length;
Sbest=Sinit;
best_length=init_length;
%% 主循环
gen=1;
BestL=zeros(MAXGEN,1);                          %记录每次迭代过程中全局最优个体的总距离
while gen<=MAXGEN
    %% “破坏”解
    [Sdestroy,removed]=destroy(Scurr);
    %% “修复”解
    [Srepair,repair_length]=repair(removed,Sdestroy,dist);
    if repair_length<curr_length
        Scurr=Srepair;
        curr_length=repair_length;
    end
    if curr_length<best_length
        Sbest=Scurr;
        best_length=curr_length;  
    end
    %% 打印当前代全局最优解
    disp(['第',num2str(gen),'代最优路线总长度 =  ' num2str(best_length)])
    BestL(gen,1)=best_length;
    %% 计数器加1
    gen=gen+1;
end
str2=['搜索完成！ 最优路线总长度 =  ' num2str(best_length)];
disp(str2)
%% 画出优化过程图
figure;
plot(BestL,'LineWidth',1);
title('优化过程')
xlabel('迭代次数');
ylabel('总距离');
%% 画出全局最优路线图
plot_route(Sbest,x,y);
toc

五、实验结果展示

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