模拟退火算法与旅行商问题

模拟退火算法的学习
模拟退火算法（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) ''])

代码中的难点：

1. 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。
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中的位置。
3. 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