模拟退火算法与旅行商问题
模拟退火算法的学习
模拟退火算法(Simulate Anneal,SA)是一种通用概率演算法,用来在一个大的搜寻空间内找寻命题的最优解。模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的。V.Černý在1985年也独立发明此演算法。模拟退火算法是解决TSP问题的有效方法之一
TSP问题:旅行商问题,即TSP问题(Traveling Salesman Problem)又译为旅行推销员问题、货郎担问题,是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市,他必须选择所要走的路径,路径的限制是每个城市只能拜访一次,而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。
![Alt](https://avatar.csdn.net/7/7/B/1_ralf_hx163co
// An highlighted block
%清空环境中的变量
clc,clear
tic
% 迭代次数初值
iter = 1;
%温度衰减系数
a=0.99;
%初始温度
t0=120;
%最后温度
tf=1;
t=t0;
%Markov链长度
Markov=5000;
%读入城市的坐标
load DATA.txt
city=DATA(:,2:3);
%城市距离初始化
n = size(city,1);
D = zeros(n,n);
for i = 1:n
for j = 1:n
D(i,j) = sqrt(sum((city(i,:) - city(j,:)).^2));
end
end
%初始化路径向量,一行n列,初始化的路径就是:1->2->3->4->5->6->...
route=1:n;
route_new=route;
best_length=Inf;
Length=Inf;
best_route=route;
%%
while t>=tf%循环是否进行的标准就是看是否冷却完毕
for j=1:Markov
%进行扰动,生成新的序列route_new;
if (rand<0.7)
%交换 route_new中两个数的的顺序
ind1=0;ind2=0;%仅仅是初始化,0没有意义。
while(ind1==ind2&&ind1>=ind2)
%ceil:取比它大的最小整数,即朝正无穷方向取整,如ceil(-1.3)=-1; ceil(1.3)=2;ceil(-1.8)=-1,ceil(1.8)=2
ind1=ceil(rand*n);
ind2=ceil(rand*n);
end %while(ind1==ind2&&ind1>=ind2)
temp=route_new(ind1);
route_new(ind1)=route_new(ind2);
route_new(ind2)=temp;
else%若:rand>=0.7
ind=zeros(3,1);%ind=[0 0 0 ]T
L_ind=length(unique(ind));%L_ind = ind的秩(向量的秩)
while (L_ind<3)%秩小于3,也就是ind中存在重复的数值
ind=ceil([rand*n rand*n rand*n]);
L_ind=length(unique(ind));
end%while (L_ind<3)
%本循环结束之时,L_ind的三个元素都是不一样的
%下面继续
ind0=sort(ind);%sort(X)对X的元素进行升序排序;当X是矩阵时,sort(X)对X的每一列进行升序排序;
%上面一行代码运行以后,ind0中存储的是升幂排列的ind向量;
a1=ind0(1);b1=ind0(2);c1=ind0(3);
%也就是说:ind = [a1,a2,a3]';
%下面的代码将route0的两个区间进行了对调
route0=route_new;
route0(a1:a1+c1-b1-1)=route_new(b1+1:c1);
route0(a1+c1-b1:c1)=route_new(a1:b1);
route_new=route0;
end %if (rand<0.7)
%计算路径的距离,Length_new
length_new = 0;
Route=[route_new route_new(1)];%这里是将第一个元素在结尾复现。所谓首尾相接
for j = 1:n
length_new = length_new+ D(Route(j),Route(j + 1));
%D是前面确定的一个n*n大小的矩阵。 D(i,j) = sqrt(sum((city(i,:) - city(j,:)).^2));
%存储着city i j的距离
end% for j = 1:n
if length_new<Length
Length=length_new;
route=route_new;
%对最优路线和距离更新
if length_new<best_length
iter = iter + 1;
best_length=length_new;
best_route=route_new;
end%if length_new<best_length
else
if rand<exp(-(length_new-Length)/t)%这里是说,如果没有达到最优解,仍然以小概率接受这个非最优解。进而更新route Length
%这里用到了评价函数
route=route_new;
Length=length_new;
end% if rand<exp(-(length_new-Length)/(t))
end%if length_new<Length
route_new=route;
end %for j=1:Markov
t=t*a;%降温。
t
end%while t>=tf
%--------------------------------------------------------------------------
%% 结果显示
toc
Route=[best_route best_route(1)];
plot([city(Route ,1)], [city(Route ,2)],'o-');
disp('最优解为:')
disp(best_route)
disp('最短距离:')
disp(best_length)
disp('最优解迭代次数:')
disp(iter)
for i = 1:n
%对每个城市进行标号
text(city(i,1),city(i,2),[' ' num2str(i)]);
end
xlabel('城市位置横坐标')
ylabel('城市位置纵坐标')
title(['理论值:629,实际值::' num2str(best_length) ''])
代码中的难点:
- ceil
Y = ceil(X)将X的每个元素四舍五入到大于或等于该元素的最接近整数。
Y = ceil(t)将duration数组t的每个元素四舍五入到大于或等于此元素的最接近的秒数。
Y = ceil(t,unit)将t的每个元素四舍五入到大于或等于此元素的最接近的数(使用指定的时间单位)。
Matlab取整函数有: fix, floor, ceil, round.具体应用方法如下:
fix朝零方向取整,如fix(-1.3)=-1; fix(1.3)=1;
floor,顾名思义,就是地板,所以是取比它小的整数,即朝负无穷方向取整,如floor(-1.3)=-2; floor(1.3)=1;floor(-1.8)=-2,floor(1.8)=1
ceil,与floor相反,它的意思是天花板,也就是取比它大的最小整数,即朝正无穷方向取整,如ceil(-1.3)=-1; ceil(1.3)=2;ceil(-1.8)=-1,ceil(1.8)=2
round四舍五入到最近的整数,如round(-1.3)=-1;round(-1.52)=-2;round(1.3)=1;round(1.52)=2。 - unique
去掉矩阵中重复的元素
(1)B = unique(A)
获取矩阵A 的不同元素构成的向量,其中B可能是行向量也可能是列向量。
(2)B = unique(A,‘rows’)
获取矩阵A的不同行向量构成的矩阵。
(3)[ C,IA,IC ] = unique(A)\unique(A,‘rows’)
IA为矩阵C中的元素在矩阵A中的位置,
IC为矩阵A中的元素在矩阵C中的位置。 - sort(X)对X的元素进行升序排序;当X是矩阵时,sort(X)对X的每一列进行升序排序;
运行结果:
下面4图为运行结果,本计算数据采用:eil101,理论最短距离:629。
数据下载地址:https://wwwproxy.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/tsp/#opennewwindow
参考最优解:https://blog.csdn.net/encaidx/article/details/8646740