秒懂

爬山算法（Hill Climbing）

为了了解退火算法，这里先介绍爬山算法作为对比来抛砖引玉。

爬山算法是一种简单的贪心搜索算法，该算法每次从当前解的临近解空间中选择一个最优解作为当前解,直到达到一个局部最优解。

如下图所示，红色曲线是我们的优化函数，现在我们的目标是找到该函数的最大值。假设初始位置是D，那么很显然根据爬山算法，我们的解会逐渐朝着E点移动，当到达E点时，由于在其周围邻域内并没有比E更高的点，算法即刻停止。但显然可以看到此时陷入了局部最优解当中，无法找到全局最优解。

退火算法

那么模拟退火算法是如何解决这一问题的呢？

模拟退火其实也是一种贪心算法，但是它的搜索过程引入了随机因素。

模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。以图1为例，模拟退火算法在搜索到局部最优解E后，会以一定的概率接受到F的移动。也许经过几次这样的不是局部最优的移动后会到达G点，进而上升到H点，于是就跳出了局部最大值E。

详解

算法来源

模拟退火算法来源于固体退火原理，是一种基于概率的算法，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。
模拟退火算法(Simulated Annealing，SA)最早的思想是由N. Metropolis等人于1953年提出。1983年，S. Kirkpatrick 等成功地将退火思想引入到组合优化领域。它是基于Monte-Carlo迭代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。

数学推导

来源：metropolis准则（蒙特卡洛准则）

本质思想为以概率接受新状态：

• 在温度为T时，设当前状态为i，新状态为j，若$E_j$<$E_i$，则接受 j 为当前状态；
• 否则，若概率 $p=e^{-(Ej-Ei)/KT}$大于[0,1)区间的随机数，则仍接受状态 j 为当前状态；若不成立，则保留状态 i 为当前状态。
  $p=e^{-(Ej-Ei)/KT}$：在高温下，可接受与当前状态能量差较大的新状态；在低温下，只接受与当前状态能量差较小的新状态。
• 公式中的k为玻尔兹曼常量，系热力学的一个基本常量，数值为$1.380649\times 10-23J/K$

算法流程

首先，我们要确定该算法一共有三个参数：$T_0$、$T_k$、$delta$。

• $T_0$表示初始的模拟温度值，在这里我们简化了蒙特卡洛准则中的概率公式为$p=e^{-(E_j-E_i)/T}$，每一轮迭代时以概率1或p接受下一状态。
• $T_k$表示终止的温度值，当温度达到该值时算法自动停止。
• $delta$表示每个轮次T值的更新权重，$T_i=delta\ast T_{i-1}$，这里是由于T如果过大，就会导致退火太快，达到局部最优值就会结束迭代，如果取值较小，则计算时间会增加，实际应用中采用退火温度表，在退火初期采用较大的T值，随着退火的进行，逐步降低。

接着是算法流程

1. 选择初始状态$S(i_0,E_{i_0})$，确定$T_0$、$T_k$、$delta$。
2. 由某种产生函数产生一个新的状态$S(i_1,E_{i_1})$
3. 按照Metropolis准则判断是否接受新状态
4. 更新Metropolis准则中的T值：令$T=T\ast delta$
5. 若$T<=T_k$，停止算法；反之重新进入步骤2

算法优劣

• 突破了爬山算法的局限性，能够获得全局最优解
• 初始解和最终解都是随机的，拥有很高的鲁棒性
• 迭代次数往往受delta的影响，delta越大，冷却速度越慢，获得最优解的搜索时间就越长；繁殖可能速度较快收敛，但有可能直接跳过最优解。

例题

TSP问题：一名商人要到若干城市去推销商品，已知城市个数和各城市间的路程（或旅费），要求找到一条从城市1出发，经过所有城市且每个城市只能访问一次，最后回到城市1的路线，使总的路程（或旅费）最小。

import numpy as np
# 位置矩阵
coordinates = np.array([[565.0, 575.0], [25.0, 185.0], [345.0, 750.0], [945.0, 685.0], [845.0, 655.0],
                        [880.0, 660.0], [25.0, 230.0], [525.0, 1000.0], [580.0, 1175.0], [650.0, 1130.0],
                        [1605.0, 620.0], [1220.0, 580.0], [1465.0, 200.0], [1530.0, 5.0], [845.0, 680.0],
                        [725.0, 370.0], [145.0, 665.0], [415.0, 635.0], [510.0, 875.0], [560.0, 365.0],
                        [300.0, 465.0], [520.0, 585.0], [480.0, 415.0], [835.0, 625.0], [975.0, 580.0],
                        [1215.0, 245.0], [1320.0, 315.0], [1250.0, 400.0], [660.0, 180.0], [410.0, 250.0],
                        [420.0, 555.0], [575.0, 665.0], [1150.0, 1160.0], [700.0, 580.0], [685.0, 595.0],
                        [685.0, 610.0], [770.0, 610.0], [795.0, 645.0], [720.0, 635.0], [760.0, 650.0],
                        [475.0, 960.0], [95.0, 260.0], [875.0, 920.0], [700.0, 500.0], [555.0, 815.0],
                        [830.0, 485.0], [1170.0, 65.0], [830.0, 610.0], [605.0, 625.0], [595.0, 360.0],
                        [1340.0, 725.0], [1740.0, 245.0]])
# 点的总数量
num = coordinates.shape[0]
# 两点之间的距离矩阵
distance = np.zeros((num, num))
for i in range(num):
    for j in range(i, num):
        distance[i][j] = distance[j][i] = np.linalg.norm(coordinates[i] - coordinates[j])
# 三大超参数
delta = 0.99
t = 100
tk = 1

solution_new = np.arange(num)
solution_current = solution_new.copy()
value_current = 99000  # np.max这样的源代码可能同样是因为版本问题被当做函数不能正确使用，应取一个较大值作为初始值
solution_best = solution_new.copy()
value_best = 99000  # np.max
result = []  # 记录迭代过程中的最优解
while t > tk:
    for i in np.arange(1000):
        # 下面的两交换和三角换是两种扰动方式，用于产生新解
        if np.random.rand() > 0.5:  # 交换路径中的这2个节点的顺序
            while True:  # 产生两个不同的随机数
                loc1 = np.int(np.ceil(np.random.rand() * (num - 1)))
                loc2 = np.int(np.ceil(np.random.rand() * (num - 1)))
                if loc1 != loc2:
                    break
            solution_new[loc1], solution_new[loc2] = solution_new[loc2], solution_new[loc1]
        else:  # 三交换
            while True:
                loc1 = np.int(np.ceil(np.random.rand() * (num - 1)))
                loc2 = np.int(np.ceil(np.random.rand() * (num - 1)))
                loc3 = np.int(np.ceil(np.random.rand() * (num - 1)))
                if (loc1 != loc2) & (loc2 != loc3) & (loc1 != loc3):
                    break
            # 下面的三个判断语句使得loc1<loc2<loc3
            if loc1 > loc2:
                loc1, loc2 = loc2, loc1
            if loc2 > loc3:
                loc2, loc3 = loc3, loc2
            if loc1 > loc2:
                loc1, loc2 = loc2, loc1
            # 下面的三行代码将[loc1, loc2)区间的数据插入到loc3之后
            tmplist = solution_new[loc1:loc2].copy()
            solution_new[loc1:loc3 - loc2 + 1 + loc1] = solution_new[loc2:loc3 + 1].copy()
            solution_new[loc3 - loc2 + 1 + loc1:loc3 + 1] = tmplist.copy()
        value_new = 0
        for j in range(num - 1):
            value_new += distance[solution_new[j]][solution_new[j + 1]]
        value_new += distance[solution_new[0]][solution_new[51]]
        if value_new < value_current:  # 接受该解
            value_current = value_new
            solution_current = solution_new.copy()

            if value_new < value_best:
                value_best = value_new
                solution_best = solution_new.copy()
        else:  # 按一定的概率接受该解
            if np.random.rand() < np.exp(-(value_new - value_current) / t):
                value_current = value_new
                solution_current = solution_new.copy()
            else:
                solution_new = solution_current.copy()
    t = delta * t
    result.append(value_best)
print(solution_best)