背景

这个 TP 的背景是这样的：有一个商务旅行团需要去 n 个城市，怎样安排他们的行程可以使行程变得更短？行程的起点和终点都在巴黎。为了解决这个问题，可以使用两种算法：第一种就是最简单的求解 n 个城市 n！种行程的各个距离，然后选取最短的行程；第二种使用la méthode du recuit simulé 模拟退火算法。

关于模拟退火算法

模拟退火算法的思想借鉴于固体的退火原理，当固体的温度很高的时候，内能比较大，内部粒子处于快速无序运动。当温度慢慢降低的过程中，固体的内能减小，慢慢趋于有序。最终，当固体处于常温时，内能达到最小。此时，粒子最为稳定。

模拟退火算法从某一较高的温度出发，这个温度称为初始温度，伴随着温度参数的不断下降，算法中的解趋于稳定。但是，可能这样的稳定解是一个局部最优解。此时，模拟退火算法中会以一定的概率跳出这样的局部最优解，以寻找目标函数的全局最优解。

导入所需的模块

import numpy as np
import time
import itertools
import matplotlib.pyplot as plt

单个路程的总距离计算

# 两个输入
# city_list : 需要去的城市列表
# city_dict : 城市的坐标字典
# 一个输出
# 路程的总距离
    distance = 0 # 初始化距离
    distance+=np.sqrt(city_dict[city_list[0]]['x']**2+city_dict[city_list[0]]['y']**2)  # 就是求两个坐标点之间的距离，先求第一个（巴黎的坐标为（0，0））
    # 循环求并累加
    for i in range(len(city_list)-1):
        x1=city_dict[city_list[i]]['x']
        y1=city_dict[city_list[i]]['y']
        x2=city_dict[city_list[i+1]]['x']
        y2=city_dict[city_list[i+1]]['y']
        distance=distance+np.sqrt((x1-x2)**2+(y1-y2)**2)
    # 求最后一个（巴黎的坐标为（0，0））
    distance+=np.sqrt(city_dict[city_list[len(city_list)-1]]['x']**2+city_dict[city_list[len(city_list)-1]]['y']**2)
    return distance  # 返回总距离

随机打乱原有路程

# 一个输入：list_in 当前旅行团的行程列表
# 一个输出：尝试路程列表
# np.random.permutation：输入一个数或者数组，生成一个随机序列，对多维数组来说是多维随机打乱
# np.asarray：从已有的数组创建数组
    return np.random.permutation(np.asarray(list_in).tolist())

计算退火新的温度（降温）

# 输入：h>0 计算参数，浮点数。h越小，算法越有风险；越大算法收敛时间越长。
# 输入：k 算法参数，整数。
# 输入：ind 算法的迭代次数，整数。
# 输入：Temp 当前算法的温度
# 输出：新的参数k
# 输出：新的温度Temp
    while ind <= np.exp((k-1)*h) or ind>np.exp(k*h):
         k+=1
         Temp=1.0/k

    return k,Temp

绘图函数

    chemin=['Paris']+chemin+['Paris']
    plt.plot([dico[ville]['x'] for ville in chemin],[dico[ville]['y'] for ville in chemin],'bo-') # 其他城市用蓝色表示
    plt.plot(dico['Paris']['x'],dico['Paris']['y'],'ro') # 巴黎用红色表示

    plt.xlim([-500.0,500.0]) # 坐标范围
    plt.ylim([-800.0,300.0])
    plt.title("Trajet")
    ax=plt.gca()
    ax.set_aspect('equal')

def plot_temp(Temp_list):
    plt.figure()
    plt.plot(Temp_list)
    plt.xlabel('$n$')
    plt.ylabel('$T$')
    plt.title(u'Profil de température')
    plt.grid()

基本参数

# 城市列表，字典
dico=dict()
dico['Lille']      ={'x':52.0, 'y':197.0}
dico['Orléans']    ={'x':-33.0, 'y':-105.0}
dico['Lyon']       ={'x':185.0, 'y':-343.0}
dico['Paris']      ={'x':0.0, 'y':0.0}
dico['Marseille']  ={'x':225.0, 'y':-617.0}
dico['Strasbourg'] ={'x':403.0, 'y':-31.0}
dico['Rennes']     ={'x':-300.0, 'y':-88.0}
dico['Metz']       ={'x':285.0, 'y':30.0}
dico['Bordeaux']   ={'x':-213.0, 'y':-448.0}
dico['Perpignan']  ={'x':40.0, 'y':-688.0}
dico['Cherbourg']  ={'x':-289.0, 'y':86.0}

parcours=np.random.permutation(['Marseille',
                                'Lyon',
                                'Rennes',
                                'Lille',
                                'Perpignan',
                                'Cherbourg']).tolist()
# 城市越多，传统的算法将花费更多时间。</pre>

用传统解决办法

# 计算所有的可能路线，如果是5个城市就是5！种情况；如果是10个城市就是10！种情况。
t1=time.time()   #计时开始

permutations=list()
for i in itertools.permutations(parcours):  #排列顺序可变，元素不重复:每一项用长度为r的元组表示
    permutations.append(i)      #所有的可能情况，一种种情况往上叠。
print('Nombre de trajets étudiés : ', len(permutations))

cost=list()
for i in range(0,len(permutations)-1):
    parcours_chaque=list(permutations[i])
    parcours_chaque=['Paris']+parcours_chaque+['Paris']
    print(parcours_chaque)
    print(cost_function(parcours_chaque,dico)) #传到上方的cost_function算法里
    cost.append(cost_function(parcours_chaque,dico)) #所有的可能情况，一种种情况往上叠。
t2=time.time() # 计时结束

# 打印结果
cost=np.asarray(cost)
print('Trajet le plus court :')
for ind, item in zip(range(1,len(parcours)+1),permutations[cost.argmin()]):
    print(ind,' ',item)
print('Distance totale : ', cost.min())
print('Temps de calcul : ',t2-t1)

用模拟退火法

candidate=parcours
# 先设置参数
itermax=10000
hpar=3.0
kpar=1
Temp=1.0/kpar
Temp_list=[Temp]

t1=time.time()   #计时开始
for ind in range(itermax):
    #计算新的行程尝试（就是根据现有行程随机生成一个新的）
    new_candidate=compute_new_candidate(candidate)
    #计算新的行程和旧的行程的路程差距
    delta=cost_function(new_candidate,dico)-cost_function(candidate,dico)
    #如果新的候选路线更长，则仍可以一定的概率接受它。
    if delta <= 0:
        candidate = new_candidate
    else:
        u=np.random.uniform(0,1) # 产生一个0到1之间的随机数
        if np.exp(-delta/Temp) >= u :
            candidate = new_candidate
    # 降温
    kpar, Temp = compute_Temp(hpar,kpar,ind+2,Temp)
    Temp_list.append(Temp)
t2=time.time()

# 打印结果
print('Trajet le plus court :')
for ind, item in zip(range(1,len(parcours)+1),candidate):
    print(ind,' ',item)
print('Distance totale : ', cost_function(candidate,dico))
print('Temps de calcul : ',t2-t1)

# 绘制图表
plot_chemin(parcours,dico)
plot_chemin(candidate,dico)

plot_temp(Temp_list)
plt.show()

参考资料

https://blog.csdn.net/google19890102/article/details/45395257