【代码】模拟退火-求解TSP问题

本站原创文章，转载请说明来自《老饼讲解-BP神经网络》www.bbbdata.com

本节展示一个使用模拟退火求解TSP问题的代码示例，

通过本文可以更具体地了解模拟退火中的用途以及实现方式

   01. 模拟退火求解TSP问题-思路分析    

本节展示如何使用模拟退火算法来解决TSP问题

  问题：TSP旅行商问题  

什么是TSP问题
有N个城市 ，要求穿过每个城市，最终回到出发城市，
问：寻找一条路径，使旅行总距离最小

TSP问题的数学表述
假设有5个城市，用矩阵D表示城市之间的距离,
  代表第i个城市到第j个城市的距离      
路径X可以表示为序列形式： 
它代表先到城市1，再到城市5，再到3.....最后由4回到1.
而总距离就是按这顺序行走的总距离
即求一 个1到N的组合X,令总距离F最小

✍️TSP问题的经典性
TSP问题有两个特有的经典特点
1. 不能穷举                
它的解空间是非常大的，如果有N个城市，那么就会有N！个解
这就不能用一般的穷举法来寻找最优解
2. x与f(x)的非纯数学性
它的目标函数 f(x)  和解 x 而更倾向于一系列描述性计算
而不是一般问题那样具有严格的数值表达式
这导致它不能用梯度下降一类的算法来求解

     模拟退火解决TSP问题的思路     

用模拟退火解决该问题时，主要需要确定以下两点：
1.能量函数是什么           
 在TSP问题中，按路线x逐个行走的总距离就是我们所要优化的能量函数
   
2.邻域中获得新解的方式
 这里我们对原解x序列上的任意两个位置进行交换，作为新解的生成方法
示例如下：
 

   02. 模拟退火求解TSP问题-代码实现    

本节展示模拟退火求解TSP问题的代码实现

      模拟退火求解TSP问题-代码实现     

根据上述思路与模拟退火算法流程，
用模拟退火求解TSP问题的代码实现如下：

%代码说明：模拟退火求解TSP旅行商问题
%来自《老饼讲解神经网络》www.bbbdata.com ,matlab版本：2018a
clc;clear all;close all 
cityLoc = [3.64,2.68  ;...  %beijing
           4.18,1.76  ;...  %shanghai
           3.71,2.60  ;...  %tianjin
           1.33,3.30  ;...  %wulumuqi
           4.20,2.96  ;...  %shenyang
           3.92,1.82  ;...  %nanjing
           2.55,1.64  ;...  %chengdu
           2.37,1.02  ;...  %kunming
           3.43,2.09  ;...  %zhengzhou
           3.54,0.70  ;...  %xianggang
           3.51,1.62  ;...  %wuhan
           3.44,0.80  ;...  %guangzhou
           3.24,2.77  ;...  %huhehaote
           2.38,2.32  ;...  %xiling
           2.56,2.24  ;...  %lanzhou
           3.01,2.03  ;...  %xian
           2.79,2.51  ;...  %yinchuan
           4.03,1.16  ;...  %fuzhou
           3.47,0.70  ;...  %aomen
           1.30,1.69] ;     %lasha
distMat = dist(cityLoc') ;                                                % 计算城市距离
N = size(cityLoc,1);                                                      % 城市个数
rng(888)                                                                  % 为方便复现结果,设定随机种子
															             
x = randperm(N);                                                          % 随机初始化解
E = distMat(x(end),x(1));                                                 % 初始化解的能量
for i=1:length(x)-1                                                       
    E = E+ distMat(x(i),x(i+1));                                          % 累加各个城市之间的距离
end


%---------温度初始化可用此方法(取N次邻域能量差均值的10倍)-------------
dE     = 0;                                                               % 计算平均能量差
try_cn = 20;                                                              % 尝试次数
for k = 1:try_cn                                                          
    pair =  randperm(N);                                                  % 随机生成要交换的两个位置
    new_x = x;  % 将旧解copy过来                                          
    new_x([pair(1),pair(2)])= new_x([pair(2),pair(1)]);                   % 将选择的两个位置进行交换
    new_E = distMat(new_x(end),new_x(1));                                 % 初始化新解的能量
    for i = 1:length(x)-1                                                   
        new_E = new_E+ distMat(new_x(i),new_x(i+1));                      % 累加各个城市之间的距离
    end                                                                   
    dE = dE+ abs(E-new_E)/try_cn;                                         % 更新平均能量差
end                                                                       
teamp = dE*10 ;                                                           % 温度初始化
														                  
														                  
%-------------外循环-------------------------                             
for t=1:100                                                               
    minE = E;                                                             % 记录本轮的最小能量
    maxE = E;                                                             % 记录本轮的最大能量
    %--------内循环-----------------------                                
    for k  =1:500                                                         
        pair  =  randperm(N);                                             % 随机生成要交换的两个位置
        new_x = x;  % 将旧解copy过来                                      
        new_x([pair(1),pair(2)])= new_x([pair(2),pair(1)]);               % 将选择的两个位置进行交换
        new_E = distMat(new_x(end),new_x(1));                             % 初始化新解的能量
        for i = 1:length(x)-1                                             
            new_E = new_E+ distMat(new_x(i),new_x(i+1));                  % 累加各个城市之间的距离
        end                                                               
        %--------按概率接受新解-------------------                        
        if(exp((E-new_E)/teamp)>rand())                                   % 如果接受新解
            x = new_x;                                                    % 更新解
            E = new_E;                                                    % 更新能量
            minE = min(E,minE);                                           % 更新本轮的最小能量
            maxE = max(E,maxE);                                           % 更新本轮的最大能量
        end
    end
    disp(['第',num2str(t),'轮:E=',num2str(E),',teamp:',num2str(teamp)])   % 打印当前温度与能量
    disp(['       minE=',num2str(minE),',maxE:',num2str(maxE)])           % 打印最小最大能量
    hold off
    plot(cityLoc(:,1),cityLoc(:,2),'o','MarkerEdgeColor','k','MarkerFaceColor','y','MarkerSize',10,'LineWidth',2)
    hold on
    plot(cityLoc([x,x(1)],1),cityLoc([x,x(1)],2))
    drawnow
    %--------降温，并检查是否退出循环--------
    teamp = teamp*0.9;                                                    % 降温
    if((maxE-minE)/maxE<0.001)                                            % 如果最优最差解并没太大差异
        break                                                             % 退出循环
    end
end

     运行结果     

运行结果如下：

20个城市，路径的组合有20！多种 ，是非常庞大的
但经模拟退火求解，可以看到，在有限的步数内，已经能自动找到一条不错的路径了
这展示了模拟退火的简单与强大！

 End