【代码】模拟退火(SA)-优化TSP

TSP问题是智能算法中的经典问题，这里我们先看看模拟退火如何解决TSP问题。

一、模拟退火-解决TSP问题

1.1. TSP问题

我们简单回顾一下TSP(旅行商)问题，TSP问题就是有N个城市 ，旅行商需求穿过每个城市，最终回到出发城市：

现在要寻找一条路径，使旅行总距离最小。

那么，假设有5个城市，TSP问题实际就是求一个序列X，例如，使得按X所代表的路径去走时，总距离最小。

1.2. 模拟退火解决TSP问题

用模拟退火解决该问题时，主要需要确定以下两点：

1. 能量函数是什么

在TSP问题中，按路线x逐个行走的总距离就是我们所要优化的能量函数

2. 邻域中获得新解的方式

这里我们采用如下的方式来获得新解。

如图所示，通过对原解x序列上的任意两个位置进行交换，作为新解的生成方法，例如图中的[2,4,6,1,3,5]，交换了4和1的位置，就得到了[2,1,6,4,3,5]。

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

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

% 代码说明：模拟退火求解TSP旅行商问题
% 来自《老饼讲解神经网络》www.bbblearn.com ,matlab版本：2018a
clc;clear all;close all 
setdemorandstream(888);                                                   % 为方便复现结果,设定随机种子
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);                                                      % 城市个数
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                                                             % 迭代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                                                              % 先hold off
    plot(cityLoc(:,1),cityLoc(:,2),'o','MarkerEdgeColor',...            
	     'k','MarkerFaceColor','y','MarkerSize',10,'LineWidth',2)         % 画出城市
    hold on                                                               % 先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！多种 ，是非常庞大的，但经模拟退火求解，可以看到，在有限的步数内，已经能自动找到一条不错的路径了，这展示了模拟退火的简单与强大！

总结

这里我们展示了模拟退火如何解决TSP问题，可以看到，智能优化算法并不像梯度下降那样，只局限于连续函数的优化，它可优化任何目标函数！只要能根据X计算出y的值就行了！这使得它们特别的万金油，啥问题都能用它来优化。