模拟退火算法（Simulated Annealing)详解_鲁棒的模拟退火算法

问题提出

旅行商问题（Travelling Salesman Problem, 简记TSP，亦称货郎担问题)：设有n个城市和距离矩阵D=[dij]，其中dij表示城市i到城市j的距离（i，j=1，2 … n），则问题是要找出遍访每个城市恰好一次的一条回路并使其路径长度为最短。
求极值问题，人们经常遇到这样的问题：在某个定义域S内，求某个函数f(x)的最小值，形式化为Min f(x)，x属于S。这是一个优化问题，根据f(x)的形式不同，有很多的优化算法来解决这类问题，简单的有穷举法，图解法，数学分析法（求导数法）等精确算法，如果很难精确求得，还有很多的近似求解法，如贪心法（如爬山法，最速下降法，梯度下降法），随机模拟方法（MCMC等）。本文将介绍的模拟退火方法属于随机模拟方法，但是是可以求得精确解的（概率为1求得全局优化解）

模拟退火思想

为了解决局部最优解问题， 1983年，Kirkpatrick等提出了模拟退火算法（Simulated Annealing）能有效的解决局部最优解问题。我们知道在分子和原子的世界中，能量越大，意味着分子和原子越不稳定，当能量越低时，原子越稳定。‘退火’是物理学术语，指对物体加温在冷却的过程。模拟退火算法来源于晶体冷却的过程，如果固体不处于最低能量状态，给固体加热再冷却，随着温度缓慢下降，固体中的原子按照一定形状排列，形成高密度、低能量的有规则晶体，对应于算法中的全局最优解。而如果温度下降过快，可能导致原子缺少足够的时间排列成晶体的结构，结果产生了具有较高能量的非晶体，这就是局部最优解。因此就可以根据退火的过程，给其在增加一点能量，然后再冷却，如果增加能量，跳出了局部最优解，那么本次退火就是成功的，下面我们就详细讲讲他是如何在局部最优解跳出来到全局最优解的：
模拟退火算法包含两个部分即Metropolis算法和退火过程。Metropolis算法就是如何在局部最优解的情况下让其跳出来，是退火的基础。1953年Metropolis提出重要性采样方法，即以概率来接受新状态，而不是使用完全确定的规则，称为Metropolis准则，计算量较低。下面先形象的说一下，然后再导出数学公式：

D是全局最优点，假设实际在寻找最优解的时候是从A点开始，随着迭代次数更新到B的局部最优解，这时候发现更新到B时，能量比A要低，这水命接近最优解了，因此百分之百转移。状态到B后，发现B的周围能量又上升了，如果是普通的梯度下降一般是不允许继续前行的，而这里会以一定的概率跳出这个坑，这跳出的概率和当前的状态、能量都有关系。下面会详细说，如果B最终跳出来了到达C，又会继续以一定的概率跳出来，可能有人会迷惑会不会跳回之前的B呢？下面会解释，直到到达D后，就会稳定下来。所以说这个概率的设计是很重要的，下面从数学方面进行解释。

假设前一个状态为x(n),系统根据某一指标（梯度下降，上节的能量），状态变为x(n+1),相应的，系统的能量由E(n)变为E(n+1),定义系统由x(n)变为x(n+1)的接受概率P为：

从上式我们可以看到，如果E(n+1)相比较E(n)能量减小了，那么这种转移就被接受（概率为1），如果能量增大了，就说明系统偏离全局最优值位置更远了，此时算法不会立刻将其抛弃，而是进行概率操作：首先在区间【0,1】产生一个均匀分布的随机数ε，如果ε < P, 其中P = EXP^(-E(n+1)-E(n)/T)，则此种转移接受，否则拒绝转移，进入下一步，往复循环。其中以能量的变化量和T进行决定概率P的大小，所以这个值是动态的。

退火算法的参数控制：

Metropolis算法是模拟退火算法的基础，但是直接使用Metropolis算法 可能会导致寻优速度太慢，以至于无法实际使用，为了确保在有限的时间收敛，必须设定控制算法收敛的参数，在上面的公式中，可以调节的参数就是T，T如果过大，就会导致退火太快，达到局部最优值就会结束迭代，如果取值较小，则计算时间会增加，实际应用中采用退火温度表，在退火初期采用较大的T值，随着退火的进行，逐步降低，具体如下：
（1）初始的温度T(0)应选的足够高，使的所有转移状态都被接受。初始温度越高，获得高质量的解的概率越大，耗费的时间越长。
（2） 退火速率。 最简单的下降方式是指数式下降：

其中λ是小于1的正数，一般取值为0.8到0.99之间。使的对每一温度，有足够的转移尝试，指数式下降的收敛速度比较慢，其他下降方式如下：

（3）终止温度
如果在若干次迭代的情况下具有可以更新的新状态或者达到用户设定的阈值，则退火完成。

模拟退火的步骤：

(1) 初始化：初始温度T(保证充分大)，初始解状态S(是算法迭代的起点)，每个T值的迭代次数L

(2) 对k=1,2 …, L做第(3)至第6步：

(3) 对当前解进行变换（例如对某些解中元素进行互换，置换，3点交叉置换等）产生相邻近的新解S′

(4) 计算增量ΔT=C(S′)-C(S)，其中C(S)为综合评价函数

(5) 若ΔT<0则接受S′作为新的当前解，否则以概率exp(-ΔT/T)接受S′作为新的当前解.

(6) 如果满足终止条件则输出当前解作为最优解，结束程序。终止条件通常取为连续若干个新解都没有被接受时终止算法。
(7) T逐渐减少，且T->0，然后转第2步。

模拟退火算法相邻解产生和接收的主要步骤：

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

第二步是计算与新解所对应的目标函数差。因为目标函数差仅由当前解的变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。

第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若ΔT<0则接受S′作为新的当前解S，否则以概率exp(-ΔT/T)接受S′作为新的当前解S。

第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上在这里插入代码片继续下一轮试验。

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

退火算法程序流程图;

模拟退火算法的优缺点

模拟退火算法的应用很广泛，可以高效地求解NP完全问题，如货郎担问题(Travelling Salesman Problem，简记为TSP)、最大截问题(Max Cut Problem)、0-1背包问题(Zero One Knapsack Problem)、图着色问题(Graph Colouring Problem)等等，但其参数难以控制，不能保证一次就收敛到最优值，一般需要多次尝试才能获得（大部分情况下还是会陷入局部最优值）。观察模拟退火算法的过程，具有以下主要优势：
1 迭代搜索效率高，并且可以并行化；
2 算法中有一定概率接受比当前解较差的解，因此一定程度上可以跳出局部最优；
3 算法求得的解与初始解状态S无关，因此有一定的鲁棒性；
4 具有渐近收敛性，已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法。

模拟退火算法同样存在以下三个参数问题：
　　(1) 温度T的初始值设置问题
　　温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高，则搜索到全局最优解的可能性大，但因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。
　　(2) 退火速度问题，即每个T值的迭代次数
　　模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说，同一温度下的“充分”搜索是相当必要的，但这也需要计算时间。循环次数增加必定带来计算开销的增大。
　　(3) 温度管理问题
　　温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中，由于必须考虑计算复杂度的切实可行性等问题，常采用如下所示的降温方式：
T=α×T.α∈(0,1).
注：为了保证较大的搜索空间，α一般取接近于1的值，如0.95、0.9。
主要参考资料
https://blog.csdn.net/weixin_42398658/article/details/84031235
https://www.jianshu.com/p/f04fe7b58080
https://blog.csdn.net/xianlingmao/article/details/7798647
附Python代码

import numpy as np
import matplotlib.pyplot as plt 
import pdb
import math
import random

coordinates = [[1304,2312],[3639,1315],[4177,2244],[3712,1399],[3488,1535],[3326,1556],[3238,1229],[4196,1004],[4312,790],[4386,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]]

num = len(coordinates)
distmat = np.zeros((31,31))
for i in range(num):
    for j in range(i,num):
        distmat[i][j] = distmat[j][i] = math.sqrt((coordinates[i][0]-coordinates[j][0])**2 + (coordinates[i][1]-coordinates[j][1])**2)

#np.savetxt('distmat.txt',distmat,fmt="%d",delimiter=",")

#get the distance matrix
num = len(coordinates)
#inatiate the primary route
def initpara():
    alpha = 0.99                         #coolrate
    t = (1,100)                          #T_end and T_start 
    markovlen = 1000                     #number of iterations per temperature
    return alpha,t,markovlen

num = len(coordinates)                   #num = 31
solutionnew = np.arange(num)
#valuenew = np.max(num)
solutioncurrent = solutionnew.copy()
valuecurrent = 99000  #np.max
#print(valuecurrent)

solutionbest = solutionnew.copy()
valuebest = 99000 #np.max
 
alpha,t2,markovlen = initpara()
t = t2[1]                                      #T_start = 1000
result = []                                    #record the best result during the iterations

route = []
while t > t2[0]:                               #T_end = t2[0] = 1,
    for i in np.arange(markovlen):
        if np.random.rand() > 0.5:             #generate new route
            while True:                        #generate two different data and swap the data
                loc1 = np.int(np.ceil(np.random.rand()*(num-1)))
                loc2 = np.int(np.ceil(np.random.rand()*(num-1)))
                if loc1 != loc2:
                    break
            solutionnew[loc1],solutionnew[loc2] = solutionnew[loc2],solutionnew[loc1]
        else:                                  #triple swap
            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
					
            if loc1 > loc2:
                loc1,loc2 = loc2,loc1
            if loc2 > loc3:
                loc2,loc3 = loc3,loc2
            if loc1 > loc2:
                loc1,loc2 = loc2,loc1
 
            tmplist = solutionnew[loc1:loc2].copy()
            solutionnew[loc1:loc3-loc2+1+loc1] = solutionnew[loc2:loc3+1].copy()
            solutionnew[loc3-loc2+1+loc1:loc3+1] = tmplist.copy()
			
        valuenew = 0
        for i in range(num-1):
            valuenew += distmat[solutionnew[i]][solutionnew[i+1]]
        valuenew += distmat[solutionnew[0]][solutionnew[30]]
		
        if valuenew < valuecurrent:  
            valuecurrent = valuenew
            solutioncurrent = solutionnew.copy()
            if valuenew < valuebest:
                valuebest = valuenew
                solutionbest = solutionnew.copy()
        else:
            if np.random.rand() < np.exp(-(valuenew-valuecurrent)/t):
                valuecurrent = valuenew
                solutioncurrent = solutionnew.copy()
            else:
                solutionnew = solutioncurrent.copy()
    t = alpha*t
    result.append(valuebest)
    print(valuebest)
    print(t)
    print(solutionbest)
    print("\n")
	
plt.plot(np.array(result))
plt.ylabel("bestvalue")
plt.xlabel("t")
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98