使用模拟退火算法求解TSP问题（python实现）

使用模拟退火算法求解TSP问题

案例和代码参考的是：华中科技大学创新团队WESharp新生培训系列讲座，主讲人：周航

模拟退火（Simulated annealing）是一种通用概率算法，用来在一定时间内寻找在一个很大搜寻空间内的近似最优解。

简介

模拟退火来自冶金学的专有名词退火。退火是将材料加热后再经特定速率冷却，目的是增大晶粒的体积，并且减少晶格中的缺陷。材料中的原子原来会停留在使内能有局部最小值的位置，加热使能量变大，原子会离开原来位置，而随机在其他位置中移动。退火冷却时速度较慢，使得原子有较多可能可以找到内能比原先更低的位置。

模拟退火的原理也和金属退火的原理近似：我们将热力学的理论套用到统计学上，将搜寻空间内每一点想像成空气内的分子；分子的能量，就是它本身的动能；而搜寻空间内的每一点，也像空气分子一样带有“能量”，以表示该点对命题的合适程度。算法先以搜寻空间内一个任意点作起始：每一步先选择一个“邻居”，然后再计算从现有位置到达“邻居”的概率。

可以证明，模拟退火算法所得解依概率收敛到全局最优解。
退火算法-中文维基百科

演算步骤

初始化
由一个产生函数从当前解产生一个位于解空间的新解，并定义一个足够大的数值作为初始温度。

迭代过程
迭代过程是模拟退火算法的核心步骤，分为新解的产生和接受新解两部分：

1、由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。

2、计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
3、判断新解是否被接受，判断的依据是一个接受准则，最常用的接受准则是Metropolis准则：若Δt′<0则接受S′作为新的当前解S，否则以概率exp（-Δt′/T）接受S′作为新的当前解S。
4、 当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。

模拟退火算法与初始值无关，算法求得的解与初始解状态S（是算法迭代的起点）无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率1收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。

停止准则
迭代过程的停止准则：温度T降至某最低值时，完成给定数量迭代中无法接受新解，停止迭代，接受当前寻找的最优解为最终解。

退火方案
在某个温度状态T下，当一定数量的迭代操作完成后，降低温度T，在新的温度状态下执行下一个批次的迭代操作。

TSP实例

#31个城市的坐标
city_loc = [(1304,2312),(3639,1315),(4177,2244),(3712,1399),(3488,1535),
                 (3326,1556),(3238,1229),(4196,1004),(4312,790),(4380,570),
                 (3007,1970),(2562,1756),(2788,1491),(2381,1676),(1332,695),
                 (3715,1678),(3918,2179),(4061,2370),(3780,2212),(3676,2578),
                 (4029,2838),(4263,2931),(3429,1908),(3507,2367),(3394,2643),
                 (3439,3201),(2935,3240),(3140,3550),(2545,2357),(2778,2826),(2370,2975)]
T0 = 50000#初始温度
T_end = 15#结束温度
q = 0.98#每次退火的比例
L = 1000#每个温度的迭代次数

#两个城市的距离
def dist(a, b):
    x1 = city_loc[a][0]
    x2 = city_loc[b][0]
    y1 = city_loc[a][1]
    y2 = city_loc[b][1]
    distance = ((x2 - x1)**2 + (y2 - y1)**2)**0.5
    return distance
#路程总长 即目标函数值
def totaldistance(a):
    value = 0
    for j in range(30):
        value += dist(a[j], a[j + 1])
    value += dist(a[30], a[0])
    return value

初始化 由一个产生函数从当前解产生一个位于解空间的新解，并定义一个足够大的数值作为初始温度。

#初始化一个解 [0,1,2,3..30]
def init_ans():
    ans = []
    for i in range(31):
        ans.append(i)
    return ans

迭代过程
迭代过程是模拟退火算法的核心步骤，分为新解的产生和接受新解两部分：
1、由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。

def creat_new(ans_before):#随机选取两个节点交换位置，得到新解
    ans_after = []
    for i in range(len(ans_before)):
        ans_after.append(ans_before[i])
    cuta = random.randint(0,30)
    cutb = random.randint(0,30)
    ans_after[cuta], ans_after[cutb] = ans_after[cutb], ans_after[cuta]
    return ans_after

2、计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
3、判断新解是否被接受，判断的依据是一个接受准则，最常用的接受准则是Metropolis准则：若Δt′<0则接受S′作为新的当前解S，否则以概率exp（-Δt′/T）接受S′作为新的当前解S。
4、 当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。
停止准则
迭代过程的停止准则：温度T降至某最低值时，完成给定数量迭代中无法接受新解，停止迭代，接受当前寻找的最优解为最终解。

退火方案
在某个温度状态T下，当一定数量的迭代操作完成后，降低温度T，在新的温度状态下执行下一个批次的迭代操作。

    while T > T_end:#直到温度降到指定结束温度时结束退火过程
        for i in range(L):#在每个温度迭代L次
            newans = creat_new(ans0)#产生函数产生新解
            old_dist = totaldistance(ans0)
            new_dist = totaldistance(newans)
            df = new_dist - old_dist#计算目标函数差
            if df >= 0:#当差>0时，以一定概率接受新解
                rand = random.uniform(0,1)
                if rand < 1/(exp(df / T)):
                    ans0 = newans
                #而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。
                #else:
                #	ans0 = ans0
            else:#新解更优时直接取代前一个解
                ans0 = newans

算法全部代码

import random
import numpy as np
from math import e
from math import exp
import matplotlib.pyplot as plt
#31个城市的坐标
city_loc = [(1304,2312),(3639,1315),(4177,2244),(3712,1399),(3488,1535),
                 (3326,1556),(3238,1229),(4196,1004),(4312,790),(4380,570),
                 (3007,1970),(2562,1756),(2788,1491),(2381,1676),(1332,695),
                 (3715,1678),(3918,2179),(4061,2370),(3780,2212),(3676,2578),
                 (4029,2838),(4263,2931),(3429,1908),(3507,2367),(3394,2643),
                 (3439,3201),(2935,3240),(3140,3550),(2545,2357),(2778,2826),(2370,2975)]

T0 = 50000
T_end = 15
q = 0.98
L = 1000

#两个城市的距离
def dist(a, b):
    x1 = city_loc[a][0]
    x2 = city_loc[b][0]
    y1 = city_loc[a][1]
    y2 = city_loc[b][1]
    distance = ((x2 - x1)**2 + (y2 - y1)**2)**0.5
    return distance
#路程总长
def totaldistance(a):
    value = 0
    for j in range(30):
        value += dist(a[j], a[j + 1])
    value += dist(a[30], a[0])
    return value

#初始化一个解 [0,1,2,3..30]
def init_ans():
    ans = []
    for i in range(31):
        ans.append(i)
    return ans
#产生新解
def creat_new(ans_before):
    ans_after = []
    for i in range(len(ans_before)):
        ans_after.append(ans_before[i])
    cuta = random.randint(0,30)
    cutb = random.randint(0,30)
    ans_after[cuta], ans_after[cutb] = ans_after[cutb], ans_after[cuta]
    return ans_after

if __name__ == '__main__':
    ans0 =init_ans()
    T = T0
    cnt = 0
    trend = []
    while T > T_end:
        for i in range(L):
            newans = creat_new(ans0)
            old_dist = totaldistance(ans0)
            new_dist = totaldistance(newans)
            df = new_dist - old_dist
            if df >= 0:
                rand = random.uniform(0,1)
                if rand < 1/(exp(df / T)):
                    ans0 = newans
            else:
                ans0 = newans
        T = T * q
        cnt += 1
        now_dist = totaldistance(ans0)
        trend.append(now_dist)
        #print(cnt,"次降温，温度为：",T," 路程长度为：", now_dist)
    distance = totaldistance(ans0)
    print(distance, ans0)
    plt.plot(trend)
    plt.show()