【数值优化】启发搜索之模拟退火法

在热力学上,退火(annealing)现象指物体逐渐降温的物理现象,温度越低,物体的能量状态会低;够低后,液体开始冷凝与结晶,在结晶状态时,系统的能量状态最低。大自然在缓慢降温(退火)时,可”找到”最低能量状态:结晶。但是如果过程过急过快,快速降温(亦称”淬炼”,quenching)时,会导致不是最低能态的非晶形。

加温时,固体内部粒子随温升变为无序状,内能增大,而徐徐冷却时粒子渐趋有序,在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小。此时,粒子最为稳定。模拟退火算法(Simulate Anneal,SA)便是基于这样的原理设计而成,主要用于组合优化问题的求解。

模拟退火算法过程

模拟退火算法从某一较高的温度出发,这个温度对应的是A点,伴随着温度参数的不断下降,算法中的解趋于稳定,但是,可能这样的稳定解是一个局部最优解,此时,模拟退火算法中会以一定的概率跳出这样的局部最优解,以寻找目标函数的全局最优解。如上图中所示,若此时寻找到了B点处的解,模拟退火算法会以一定的概率跳出这个解,如跳到了C点重新寻找,这样在一定程度上增加了寻找到全局最优解D点的可能性。
该算法可以分解为如下流程。
STEP1: 初始化。初始温度 T (充分大),初始解状态 S (是算法迭代的起点), 约定迭代次数 L。
STEP2: 对 k=1,……,L 做第 3 至第 6 步:
STEP3: 产生新解 S’
STEP4: 计算增量 ΔE =E(S′)-E(S) ,其中 E(S) 为评价函数
STEP5: 若 Δt′<0 则接受 S′ 作为新的当前解,否则以概率 exp(-ΔE/kT) 接受 S′ 作为新的当前解.

Metropolis 准则,粒子在温度T时趋于平衡的概率为 exp(-ΔE/kT),其中E为温度T时的内能,ΔE为其改变量,k为Boltzmann常数。

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

模拟退火法的代码实现

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
# coding=utf-8
import numpy as np
import matplotlib.pyplot as plt
import math

#########################################################
# 定义目标函数
#########################################################

def aimFunction(x):
# 目标函数
y = x**3-80*x**2-4*x+6
return y

x = [i/10 for i in range(1000)]
y = []

for i in range(1000):
y.append(aimFunction(x[i]))

plt.plot(x, y)
plt.show()

#########################################################
# 模拟退火算法
#########################################################

def SA(T_init, T_min, parameter_space):
'''
模拟退火算法
INPUT -> 初始温度, 最低温度(即退出循环条件), 参数空间
'''
import random
T = T_init # 初始温度
x = random.choice(parameter_space) # 随机选取初解

alpha = 0.95 # 降温系数
iterL = 100 # 搜索次数

# 外层循环, 是否满足终止条件
while T >= T_min:

# 内层循环, 是否达到迭代次数
for i in range(iterL): # 一般来说,同一温度下的"充分"搜索是相当必要的

# 当前解的内能
E = aimFunction(x)
# 当前解变化量
delta_x = random.random() - 0.5
# 新解仍要求在参数空间之内
if int(x + delta_x) in parameter_space:
xNew = x + delta_x
# 计算新解的内能
ENew = aimFunction(xNew)
# 内能变化量
delta_E = ENew - E
if delta_E < 0: # 无条件更新
x = xNew
else: # 按照概率更新
p = math.exp(-delta_E/T)
if random.random() < p:
x = xNew

# 降到新的温度
T = T*alpha
print(x, T)
return x

x = SA(1000, 0.001, range(100))
print(x, aimFunction(x))

用模拟退火算法解决实际问题

旅行商问题(TravelingSalesmanProblem,TSP)是一个经典的组合优化问题。经典的TSP可以描述为:一个商品推销员要去若干个城市推销商品,该推销员从一个城市出发,需要经过所有城市后,回到出发地。应如何选择行进路线,以使总的行程最短。

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
# coding:utf-8
import math
import time
import random

# 初始温度
T = 50000
# 最低温度
T_end = 1e-8
# 在每个温度下的迭代次数
L = 100
# 退火系数
delta = 0.98

# 31个城市的坐标
citys = [[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]]
# 存储两个城市之间的距离
d = [[0 for i in range(31)] for j in range(31)]
# 存储一条路径
ans = []

# 计算降温次数
cnt=0

# 计算两个城市之间的距离
def get_city_distance():
for i in range(len(citys)):
for j in range(i, len(citys)):
d[i][j] = d[j][i] = math.sqrt((citys[i][0]-citys[j][0])**2 + (citys[i][1]-citys[j][1])**2)

# 使用随机交换路径中两个城市的位置来产生一条新路径
def create_new(a):
i = random.randint(0, len(a)-1)
j = random.randint(0, len(a)-1)
a[i], a[j] = a[j], a[i]
return a

# 获取路径的长度
def get_route_distance(a):
dist = 0
for i in range(len(a)-1):
dist+=d[a[i]][a[i+1]]
return dist

def saa():
get_city_distance()
cnt = 0
ans = range(0, len(citys))
t = T
result = 0
while t >= T_end:
for i in range(0, L):
ans_new = create_new(list(ans))
d1, d2 = get_route_distance(ans), get_route_distance(ans_new)
de = d2 - d1
result = d1
if de<0:
ans = ans_new
result = d2
else:
if(math.e**(-de/T) > random.random()):
ans = ans_new
result = d2
t = t * delta
cnt+=1
print("路径如下:")
print(ans)
print("路径长度:")
print(result)
print("降温次数:")
print(cnt)

start = time.time()
saa()
end = time.time()
print("运行耗时:" + str(end-start) + "s")
0%