变邻域搜索算法求解TSP问题（VNS）（自用部分非原创内容）
变邻域搜索算法（Variable Neighborhood Search，简称VNS）是一种用于解决组合优化问题的启发式搜索算法。该算法由Mladenović和Hansen于1997年提出，旨在通过不断变化问题的邻域结构来寻找更好的解决方案。VNS通常用于NP困难问题，如旅行商问题（TSP）、车辆路径问题（VRP）等。
以下是变邻域搜索算法的基本思想和步骤：
- 初始解：算法从一个初始解开始，可以是随机生成的或基于其他启发式算法得到的解。
- 邻域结构：VNS使用一系列不同的邻域结构来探索潜在的解空间。每个邻域结构定义了一种改变当前解的方式。这些邻域结构可以是不同的，包括交换两个元素、插入一个元素、反转一部分元素等等。
- 邻域搜索：在每个迭代中，算法选择一个当前邻域结构，并在该结构下搜索一个新的解。搜索可以使用局部搜索方法，如模拟退火、爬山法等。如果找到了一个更好的解，则更新当前解为新的解。
- 邻域切换：如果在当前邻域结构下无法找到更好的解，算法切换到下一个邻域结构，并继续搜索。这个过程可能会在一系列不同的邻域结构之间进行多次迭代。
- 停止条件：算法可以根据一定的停止条件来终止，例如达到最大迭代次数、运行时间超过限制或者找到了一个满足要求的解。
VNS的关键思想是通过变化邻域结构来增加搜索的多样性，从而有更大的机会找到更优的解决方案。通过不断地切换邻域结构，VNS可以跳出局部最优解，进一步搜索解空间。然而，VNS的性能高度依赖于如何选择和切换邻域结构，以及在每个邻域中进行的局部搜索策略。
构造初始路线
初始解是算法开始搜索的起点，它可以通过多种方式构造，包括以下几种常见方法：
- 随机生成：可以随机生成一个初始解，这通常是一种简单而快速的方法。但是，随机生成的初始解可能会比较差，因此需要后续的搜索来改进。
- 启发式算法：可以使用其他启发式算法（如最近邻算法、插入算法、交换算法等）来生成一个初始解。这些算法通常会根据一些启发式规则来构建初始解，以期望得到相对较好的起点。
- 基于已知信息的构建：如果问题的特性允许，可以根据已知的信息来构建初始解。例如，对于VRP，可以基于顾客之间的距离或需求来构建初始路线。
初始解的质量对于变邻域搜索算法的性能至关重要。一个好的初始解可以加速算法的收敛，而一个差的初始解可能需要更多的搜索时间才能找到更优的解。因此，在实际应用中，通常会尝试不同的初始解构造方法，以找到最好的初始解来启动变邻域搜索算法的搜索过程。现在用贪婪算法来举例子
初始化一个空的解，表示还没有选择任何元素或路径。
从未选择的元素中选择一个元素，通常是基于某种启发式规则来做出选择。这个规则可能会考虑元素之间的成本、距离、权重等信息，以选择最有希望的元素。
将选定的元素添加到当前解中。
重复步骤2和步骤3，直到满足某个终止条件，例如所有元素都被选择完毕。
function [init_route,init_len]=construct_route(dist) %% 贪婪算法构造TSP的初始解 %输入dist： 距离矩阵 %输出init_route： 贪婪算法构造的初始路线 %输出init_len： init_route的总距离 N=size(dist,1);%城市数目 for i = 1:N for j = 1:N if i == j dist(i,j)=inf; end end end univisited=1:N;%初始未被安排的城市集合 visited=[];%初始已被安排的城市集合 min_dist=min(min(dist));%找出距离矩阵中的最小值 [row,col]=find(dist==min_dist);%在dist中找出min_dist所对应的行和列 first=row(1);%将min_dist在dist中所对应行序号作为起点 univisited(univisited==first)=[];%将first从unvisit中删除 visited=[visited,first];%%把first添加到visit中 pre_point=first;%%将fisrt赋值给pre_point while ~isempty(univisited) pre_dist=dist(pre_point,:); %pre_point与其它城市的距离 pre_dist(visited)=inf;%将pre_point与已经添加进来的城市之间的距离设位无穷大 [~,pre_point]=min(pre_dist); %找出pre_dist中的最小值 univisited(univisited==pre_point)=[];%将pre_point从unvisit中删除 %%第一个pre_point只是为了找到最小距离的起点序号,后面用于储存pre_dist %%第二个pre_point储存了pre_dist中的最小值然后并将其赋予[]值 visited=[visited,pre_point]; %把pre_point添加到visit中 end init_route=visited; init_len=route_length(init_route,dist); %计算init_route的总距离 end
交换操作
function route2 = swap(route1, i, j)
% 交换操作
% 比如说有6个城市，当前解为123456，随机选择两个位置，然后将这两个位置上的元素进行交换。
% 比如说，交换2和5两个位置上的元素，则交换后的解为153426。
% 输入 route1：初始路线
% 输入 i, j：两个交换点的位置索引
% 输出 route2：经过交换操作变换后的路线2
route2 = route1; % 复制初始路线到新路线
% 交换位置 i 和 j 上的元素
route2([i j]) = route1([j i]);
end
逆转操作
function route2=reversion(route1,i,j)
%% 逆转操作
%有6个城市，当前解为123456，我们随机选择两个位置，然后将这两个位置之间的元素进行逆序排列。
%比如说，逆转2和5之间的所有元素，则逆转后的解为154326。
%输入route1： 路线1
%输入i,j： 逆转点i,j
%输出route2： 经过逆转结构变换后的路线2
i1=min([i,j]);% 计算逆转起始位置的索引
i2=max([i,j]); % 计算逆转结束位置的索引
route2=route1;% 复制初始路线到新路线
% 使用逆序排列操作，将位置 i1 到 i2 之间的元素逆序排列
route2(i1:i2) = route1(i2:-1:i1);
end
插入操作
插入操作要分 i i i和 j j j的大小（如图已经很明了了）
function route2=insertion(route1,i,j)
%% 插入操作
%有6个城市，当前解为123456，我们随机选择两个位置，然后将这第一个位置上的元素插入到第二个元素后面。
%比如说，第一个选择5这个位置，第二个选择2这个位置，则插入后的解为125346。
%输入route1： 路线1
%输入i,j： 插入点i,j
%输出route2： 经过插入结构变换后的路线2
if i<j
% 如果 i 小于 j，插入操作的位置在 i 和 j 之间
%i左边的保留，i-j之间提前面去，i放j后面，j后面保留
%i=2 j=5 123456→134526
route2=route1([1:i-1,i+1:j,i,j+1:end]);
else
% 如果 i 大于 j，插入操作的位置在 j 和 i 之间
%j前面保留，i后面保留，j后面与i前面的整体的前面插入i
%i=5 j=2 123456→125346
route2 = route1([1:j, i, j+1:i-1, i+1:end]);
end
end
扰动操作
扰动操作用于生成一个局部最优解的邻域内的随机解，以探索可能的更优解。扰动操作是VNS算法中的一个重要步骤，用于引入随机性以打破当前局部最优解的限制，从而更广泛地搜索解空间。
扰动操作通常包括以下步骤：
- 选择一个邻域（neighborhood）：首先，从算法定义的一组邻域中选择一个邻域，这个邻域可以是一种特定的局部搜索操作，例如交换操作、逆转操作、插入操作等。
- 选择一个解：从当前解中随机选择一个点或一组点，这些点将成为扰动操作的目标。选择的点数和方式取决于所选的邻域和问题的特性。
- 应用扰动：对所选择的解应用扰动操作，生成一个随机的解。扰动操作的具体实现方式取决于所选的邻域。例如，如果邻域是交换操作，那么可以随机选择两个点并交换它们的位置，从而生成一个新解。
- 返回扰动后的解：将生成的随机解作为扰动后的解返回，用于接下来的局部搜索或全局搜索。
function [route_shake, len_shake] = shaking(route, dist, k)
% 扰动，随机选择当前邻域中的一个解更新当前解
% 输入 route：一条路线
% 输入 dist：距离矩阵
% 输入 k：当前邻域序号
% 输出 route_shake：扰动操作后得到的路线
% 输出 len_shake：该条路线的距离
N = numel(route); % 城市数目
% 随机选择进行操作的两个点的序号，这两个点将被用于扰动操作
%在1-N中生成一个1×2的数组
select_no = randi([1, N], 1, 2);
i = select_no(1);
j = select_no(2);
if k == 1
route_shake = swap(route, i, j); % 如果 k 等于 1，执行交换操作
elseif k == 2
route_shake = reversion(route, i, j); % 如果 k 等于 2，执行逆转操作
else
route_shake = insertion(route, i, j); % 否则，执行插入操作
end
len_shake = route_length(route_shake, dist); % 计算扰动后路线的总距离
end
路线距离计算
注意：计算时路线起点和终点需连接起来，也就是123456➡1234561,而下面的dist就是距离矩阵，在这之前已经用matlab自带的pdist函数计算了，再用squareform函数进行format distance matrix。
function len = route_length(route, dist)
% 计算一条路线总距离
% 输入 route：一条路线（城市编号的数组）
% 输入 dist：距离矩阵（城市之间的距离）
% 输出 len：该条路线总距离
n = numel(route); % 获取路线中城市的数量
route = [route route(1)]; % 将路线首尾连接，形成闭环
len = 0; % 初始化总距离为0
% 遍历路线上的每一对相邻城市，计算它们之间的距离并累加到总距离上
for k = 1:n
i = route(k); % 当前城市的编号
j = route(k+1); % 下一个城市的编号
len = len + dist(i, j); % 累加距离
end
end
距离差值计算
距离差值计算需计算三种变换对应的结果，其实可以用route_length计算，但这样每次都会调用route_length函数，会增加计算时间，所以可以用下面这些函数来缩短计算时间
交换操作后的差值
交换的差值分7种情况（这里偷个懒 就不画这么多图了，哈哈）
%% 计算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
更新距离差值函数
跟上面的cal_delata函数可以写在一起，但这样只是为了方便
function Delta1=Update1(route,dist,i,j)
%% swap操作后生成新的距离差矩阵Delta
%输入route： 一条路线
%输入dist： 距离矩阵
%输入i，j： 交换点i,j
%输出Delta1： swap操作后的距离差值的矩阵
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
function Delta2=Update2(route,dist,i,j)
%% reversion操作后生成新的距离差矩阵Delta
%输入route： 一条路线
%输入dist： 距离矩阵
%输入i，j： 逆转点i,j
%输出Delta2： reversion操作后的距离差值的矩阵
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
交换邻域搜索
function [reversion_route,reversion_len]=reversion_neighbor(route,dist,M)
%% 对route不断进行逆转操作后所得到的路线以及所对应的总距离
%输入route： 一条路线
%输入dist： 距离矩阵
%输入M： 最多进行M次邻域操作
%输出reversion_route： 对route不断进行逆转操作后所得到的路线
%输出reversion_len： reversion_route的总距离
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
function [swap_route,swap_len]=swap_neighbor(route,dist,M)
%% 对route不断进行交换操作后所得到的路线以及所对应的总距离
%输入route： 一条路线
%输入dist： 距离矩阵
%输入M： 最多进行邻域操作的次数
%输出swap_route： 对route不断进行交换操作后所得到的路线
%输出swap_len： swap_route的总距离
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;
swap_len=route_length(swap_route,dist);
end
function [insertion_route,insertion_len]=insertion_neighbor(route,dist,M)
%% 对route不断进行插入操作后所得到的路线以及所对应的总距离
%输入route： 一条路线
%输入dist： 距离矩阵
%输入M： 最多进行M次邻域操作
%输出insertion_route： 对route不断进行插入操作后所得到的路线
%输出insertion_len： insertion_route的总距离
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
主函数
流程图后面补，今天脑子累
tic
clear
clc
close all
%% 输入数据
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=1000; %外层最大迭代次数
M=60; %最多进行M次邻域操作
n=7; %邻域数目
%% 构造初始解
[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
图我用figurebest修了下，有意者可以关注下图通道这个公众号，非常好用
`
图我用figurebest修了下，有意者可以关注下图通道这个公众号，非常好用
扫一扫
5094
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
