【建模算法】基于模拟退火算法求解TSP问题（matlab求解）

【建模算法】基于模拟退火算法求解TSP问题（matlab求解）

TSP (traveling salesman problem，旅行商问题)是典型的NP完全问题，即其最坏情况下的时间复杂度随着问题规模的增大按指数方式增长，到目前为止还未找到一个多项式时间的有效算法。本文探讨了使用matlab软件，基于模拟退火算法求解TSP问题。

一、问题描述

​ 本案例以31个城市为例，假定31个城市的位置坐标如表1所列。寻找出一条最短的遍历31个城市的路径。

二、解决思路及步骤

1、什么是模拟退火？

固体的退火过程指的是将固体加热至足够高的温度，再使其慢慢冷却的过程。在加热过程中，原本有序排列的内部粒子开始无序运动，此时固体的内能不断增大；而在降温过程中，粒子的排列又逐渐有序。

而物体系统总是趋向于保持能量最低的状态，理论上讲，在退火过程的降温过程中，如果降温过程足够的慢，可将每一时刻离散化，那么每一个温度下的固体都可能达到允许范围内的最小内能状态，而如何寻找该温度及状态即为模拟退火算法的核心问题。

模拟退火算法和遗传算法、粒子群算法等智能算法作为一种通用的概率算法，该类智能算法都可用来寻找一个较大搜索空间之内的全局最优解。搜索空间之内的每一个元素都可能是问题的全局最优解，那么如何确立搜索方式、如何判断某一温度下的固体粒子内能状态便成为区分其他智能算法的核心工作。一般地，可以将温度T作为搜索空间；将实际情况中的目标函数值f看做内能E；而固体在某一温度下的状态看做问题的一个解x，即E=f(x);算法在一定的规则之下，控制温度T的降低，使固体的内能也随着最优值（内能最低）靠近，直至找到了全局最优解，就像固体退火过程一样。

如何判断内能减少？

Metropolis准则：从当前i状态发展到下一状态j，若新状态的内能小于状态i的内能即（$E_j$<$E_i$）,则接受j为新的状态，其对应的内能作为目前的最优解；否则，以概率$exp[\frac{-(E_j-E_i)}{KT}]$ 接收该新解, 其中k为Boltzmann常数，T为目前的温度。在后续的判断过程中可加入计算机的蒙特卡罗实验模拟产生概率（随机数)。

算法设计为了寻找最小的内能，那可以对应到实际情况中的全局最小解，若要寻找全局最大值，在目标函数前添加“负号”即可。

2、算法基本步骤：

3、参数选择

三、模拟退火算法求解TSP问题

1. TSP问题：TSP(Traveling salesman problem)即旅行商问题，旅行商希望在N个城市进行一次巡回旅行，可以恰好访问每一个城市一次，并且最终回到出发城市。并且要使得这次巡回旅行的总消耗最小（总距离或总花销等等），如何求这个路线？

2. 目标函数
   

3. 新解的产生
   

4. Metropolis接收准则
   

四、运行结果

1、模拟退火过程动画演示

2、优化结果

​ 优化后的路线图：

最优解:
16->23->11->24->19->17->3->18->22->21->20->25->26->28->27->30->31->29->1->15->14->12->13->7->6->5->2->10->9->8->4->16
总距离:15.3992

五、完整源代码

%模拟退火算法求解TSP问题完整源代码
clear,clc;
close all;
% 读取个各城市的坐标
X=importdata('p_xy.xlsx');
data=X.data.Sheet1(:,2:3);   %各城市坐标
num_city = size(data,1);
Initial_temp = 1000;
res = 1e-3;       % 最低温限制
ratio = 0.9;       % 降温参数，控制温度的下降
temperature = Initial_temp;
Markov_length = 3000;                % 改变解的次数
Distance_matrix = pdist2(data,data); % 矩阵的形式存储城市之间的距离
route_new = randperm(num_city);      % 新产生的解路线
Energy_current = inf;                % 当前解的能量
Energy_best = inf;                   % 最优解的能量

route_current = route_new;           % 当前解路线
route_best = route_new;              % 最优解路线
pic_num = 1;

% 外层while循环控制降温过程，内层for循环控制新解的产生。
while temperature > res
    Energy1=Energy_best;                             % 用于控制循环的结束条件
    for i = 1: Markov_length
        
        % 产生新解(对当前解添加扰动)
        if rand >0.5
            % 两点交换
            a = 0;
            b = 0;
            while (a==b)
                a = ceil(rand*num_city);
                b = ceil(rand*num_city);
            end
            temp = route_new(a);
            route_new(a) = route_new(b);
            route_new(b) = temp;
        else
            % 三变化
            factor = randperm(num_city,3);
            factor = sort(factor);      % 对三个元素排序，a<b<c
            a = factor(1);
            b = factor(2);
            c = factor(3);
            temp = route_new(a:b);
            route_new(a:a+c-b-1) = route_new(b+1:c);
            route_new(a+c-b:c) = temp;
        end
        
        Energy_new = 0;   %计算该条路线的总距离
        for j=1:num_city-1
            Energy_new = Energy_new + Distance_matrix(route_new(j),route_new(j+1));
        end
        % 回到起始点，加上首尾两个城市的距离
        Energy_new = Energy_new+Distance_matrix(route_new(1),route_new(num_city));
        
        
        % 按照Metroplis准则接收新解
        if Energy_new<Energy_current
            % 更新局部最优
            Energy_current = Energy_new;
            route_current = route_new;
            
            % 更新全局最优
            if Energy_new<Energy_best
                Energy_best=Energy_new;
                route_best = route_new;
            end
            
        else
            if rand<exp(-(Energy_new-Energy_current)/temperature)
                Energy_current = Energy_new;
                route_current = route_new;
            else
                route_new = route_current; % 否则路线不更新，保存更改之前的路线
            end
            
        end
        
    end
    
    %画图
    items = [route_new,route_new(1)];                % 从最后一个城市回到起始城市

    scatter(data(route_new,1),data(route_new,2))
    for i=1:length(route_new)
        text(data(route_new(i),1),data(route_new(i),2),num2str(items(i)))
    end
    
    hold on
    plot(data(items,1),data(items,2),'b-')
    text(1.5,3.5,['温度=',num2str(temperature)]);
    text(1.5,3.3,['总距离=',num2str(Energy_best)]);

    drawnow;
    F=getframe(gcf);
    I=frame2im(F);
    [I,map]=rgb2ind(I,256);
    hold off
    title('模拟退火算法求解TSP问题')
    xlabel('横坐标')
    ylabel('纵坐标')
    pic_num = pic_num + 1;
    if Energy1==Energy_best&&Energy_current==Energy_best
        break
    else
        temperature = temperature*ratio; %降温过程
    end
    
end

%Energy_best:最短路径；route_best:最短路径对应的最优距离（不唯一）
disp('最优解:');
R = [route_best, route_best(1)];
N = length(R);
p = num2str(R(1));
for i = 2: N
    p = [p, '->', num2str(R(i))];
end
disp(p)
disp(['总距离:',num2str(round(Energy_best,8))]);

[参考书目]

数学建模与数学实验（第5版），高等教育出版社，赵静，但琦.