模拟退火算法解决TSP问题（含MATLAB&Python两种语言代码）

​1、现代优化算法的由来

在寻找最优解的过程中，我们常常想到最简单，最直接的办法是能不能把所有解全部求出，然后再从这些解中寻找最好的那一个。

这种思路通过计算机强大的运算能力能解决一部分问题，但随着问题的进一步复杂，要处理的数据也越来越多，计算量也越来越大，计算机求出所有解也变得越来越困难。如经典的旅行商问题，如果商人要通过的城市数为100个，他能选择的方案数就有100！，这个数比10的158次方还大。如果城市数为1000个，甚至是10000个，这种计算量简直难以想象！

可如果用贪婪算法来求解，得到的往往解往往只是局部最优，难以达到全局最优。

在这种基础上就有人提出，能不能通过降低解的精度来达到减少计算量，找到一个近似最优解。这就是现代优化算法的由来。

2、模拟退火算法

2.1 模拟退火算法的基本原理

模拟退火算法出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。

退火过程：固体物质的降温过程中，固体物质内部不断进行重新排列，并逐渐排列成最低能量状态的结构。在高温状态下，由于分子的扰动能力较强，对较差状态(远离最值所对应的状态)的容忍性高，因此可以在给定状态空间内进行全局的随机搜索，从而有较高概率跳出局部极值。

算法上的优化过程：则是当前解内部不断进行重新排列，并逐渐排列成实现目标函数最小值的解。在不断优化解的过程中需要摆脱贪婪算法的局限性，能有一定的概率跳出局部最优，达到全局最优。

2.2模拟退火过程

在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。

在某一个特定温度下，进行了充分的转换之后，材料将达到热平衡。这时材料处 于状态i 的概率满足波尔兹曼分布：

其中 x 表示材料当前状态的随机变量， S 表示状态空间集合，K 是物理学中的波尔兹曼常数，T 是材料温度。
可以证明：在高温下所有状态出现都具有相同的概率；当温度降至很低时，材料会以极大概率进入最小能量状态。

Metropolis 算法用一个简单的数学模型描述了退火过程。

假设材料在状态i 之下的能量为 E(i) ，那么材料在温度T 时从状态i 进入状 态 j 就遵循如下规律：

（1）如果 E( j) ≤ E(i) ，接受该状态被转换。

（2）如果 E( j) > E(i) ，则状态转换以如下概率被接受：

2.3算法流程图

3、TSP旅行商问题

3.1 问题描述

有个旅行商，他要去中心城市A，出发到另外n个城市（已给出城市的经纬度坐标）出售商品，最后在回到中心城市A。如果旅行商想在路上花费的时间最短，请你帮他规划路线。

3.2 问题分析

模拟退火算法的实现主要可分为：解空间、新解的产生和目标函数三部分。下面我们一步一步实现

3.2.1 解空间

解空间 S 可表为{1,2,3,……n，n+1 ，n+2}的所有固定起点和终点的循环排列集合，其中1和n+2都表示中心城市，{2~n+1}表示旅行商经过的城市
如：[1 2 3 4 5 6 7 8 9 10]和[1 9 8 7 6 1 2 3 4 5 10] 表示的就是通过8个城市旅行商的两条不同路径。

另外为了计算过程更方便，我们还应提前构建表示两两城市之间的距离矩阵d

如：Dij就表示城市 i 到 城市 j 的距离

3.2.2 新解的产生

任选序号 u,v （u < v ）交换u 与v 之间的顺序，此时的新路径为：

产生新解的路径差为：

3.2.3 目标函数

即求路径的最小值 min f(1……n+2)

在后面的代码中我们计算了23个城市的路线图

注：每次运行的结果都不相同，可多运行几次，找到近似最优解。
代码中的数据可关注公众号：“数模编程手”，并在后台回复 “TSP数据” 进行提取

欢迎关注我微信公众号 ： “数模编程手” 会不定期推送一些数学建模常用的算法，并附有MATLAB和Python两种语言的代码。

MATLAB代码：

clc,clear
x=xlsread(\'题目数据.xlsx\',\'sheet1\',\'B3:B31\');
y=xlsread(\'题目数据.xlsx\',\'sheet1\',\'C3:C31\');
sj=[x y]; 
dl =[120.7015202,36.37423];
sj=[dl;sj;dl];
sj=sj*pi/180; 
d=ones(31);
for i =1:31 
for j=1:31
 mp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2));
 d(i,j)=6370*acos(mp);
end
end
long=inf;
path=[];
for j=1:300 %求较好的初始解
path0=[1,1 + randperm(29),31]; temp=0;
for i =1:30
temp = temp + d(path0(i) ,path0(i+1));
end
 if temp < long
path = path0; long = temp;
end
end
long
t=0.1^ 30;%最低温度
L = 100*30;%Markov链长度
at =0.99;%控制参数T的衰减函数参数
T=1;%退火初始温度​
while t<T
    for i=1:L
        c=2 + floor(29* rand(1,2)); %产生新解
        c=sort(c); c1=c(1);c2=c(2);%计算代价函数值的增量
        df =d(path(c1 - 1),path(c2)) +d(path(c1),path(c2 +1)) -d(path(c1 - 1),path(c1)) - d( path( c2) ,path(c2 +1));
        if df<0 %接受准则
        long= long +df;
        path= [path(1:c1 -1) ,path(c2: -1:c1),path(c2 +1:31)];
        elseif exp( -df/T) >= rand(1)
        path= [path(1:c1 -1) ,path(c2: -1:c1),path(c2 +1:31)];
        long = long+df;
        end
        T=T*at;
    end
end
path;
long%输出巡航路径及路径长度
sj=sj*180/pi;
xx=sj(path,1);yy = sj( path,2);
plot(xx,yy,\'- *\') %画出巡航路径
text(120.7015202,36.37423,\'中心城市\',\'color\',\'k\')%标注数据中心点
for i=1:29
    text(x(i),y(i),num2str(i),\'color\',\'r\')%标注各点
end
​

Python代码：

import csv
import math
import numpy as np
import random
from matplotlib import pyplot as plt
​
\"导入数据\"
filename = \'题目数据.csv\'
with open(filename) as f:
    reader = csv.reader(f)
    header_row = next(reader)
    x = []
    y = []
    longtitudes=[]
    dimensionalitys=[]
    for row in reader:
        x.append(row[1])
        y.append((row[2]))
        longtitudes = [float(x) for x in x]
        dimensionalitys = [float(y) for y in y]
longtitudes.append(longtitudes[0])
dimensionalitys.append(dimensionalitys[0])
​
\"构建距离矩阵\"
d = []
rows = 31
columns = 31
pi = 3.1415926
for row in range(rows):
    d.append([])
    for column in range(columns):
        x1 = longtitudes[row] * pi/180
        y1 = dimensionalitys[row] * pi/180
        x2 = longtitudes[column]*pi/180
        y2 = dimensionalitys[column]*pi/180
        num = 6370 * math.acos((math.cos(x1-x2)*math.cos(y1)*math.cos(y2)+math.sin(y1)*math.sin(y2)))
        d[row].append(num)
​
\"产生较好的初始路径并计算初始长度\"
path = [0]
path1 = np.random.permutation(29)+1
path.extend(path1)
path.append(30)
temp = 0
for i in range(30):
    temp = temp + d[path[i]][path[i + 1]]
long = temp
​
​
\"设置初始值\"
t = 0.1**30
T = 1
at = 0.999
L = 10*300
a = []
​
\"模拟退火过程\"
while t < T:
    for i in range(L):
        a = random.sample(range(2, 29), 2)
        a.sort()
        a1 = a[0]
        a2 = a[1]
        r = np.random.uniform(0, 1)
        \"计算新路径的差值\"
        df = d[path[a1 - 2]][path[a2 - 1]] + d[path[a1 - 1]][path[a2]] - d[path[a1 - 2]][path[a1 - 1]] - d[path[a2 - 1]][path[a2]]
        if df < 0:
            long = long + df
            change_path = path[a2 - 1:a1 - 2:-1]
            path[a1 - 1:a2] = change_path
        elif math.exp(-df/T) > r:
            long = long + df
            change_path = path[a2 - 1:a1 - 2:-1]
            path[a1 - 1:a2] = change_path
        T = T*at
​
print(long)
print(path)
​
\"画路线图\"
plt.rcParams[\'font.sans-serif\']=[\'SimHei\']
xx = []
yy = []
for i in range(31):
    aa = path[i]
    aa1 = longtitudes[aa]
    aa2 = dimensionalitys[aa]
    if aa>0 & aa<30:
        plt.text(aa1, aa2, aa)
    elif aa == 0:
        plt.text(aa1,aa2,\"中心城市\")
    xx.append(aa1)
    yy.append(aa2)
plt.plot(xx, yy)
plt.show()