3 模拟退火算法
3 模拟退火算法
3.1 概述
模拟退火算法(Simulated Annealing Algorithm,SAA)是一种模拟物理退火过程而设计出来的优化算法,退火过程实际上可以理解成一个物体的温度缓慢从高温(能量高,活跃状态)降低到低温(能量低,稳定状态)的过程。
假设当前我们要求解目标函数\(f(x)\)的最小值,如果我们把\(f(x)\)的函数值理解为物体处于状态\(x\)时的能量,那么我们的目标就变成了找到使得粒子能量最低的状态。为了完整模拟物体退火(降温)过程,我们还需要为当前时刻定义一个温度\(T\)。于是,退火算法可以大体上阐述为:
- 首先,设置一个较高的起始温度\(T_0\)并选定一个随机初始状态\(x_0\)(通常\(f(x_0)\)都会比较大)。
- 然后,我们控制温度\(T_0\)缓慢下降至某个终止温度,在温度下降期间通过某种策略产生一系列的状态\(x_1,x_2,\cdots,x_n\),并且这些状态原则上应该满足\(f(x_1)\ge f(x_2)\ge\cdots\ge f(x_n)\)。
- 最后,我们就得到了一个能量较低的状态\(x_n\)。
当然,这只是一种较为粗略的描述。
3.2 算法具体细节
初始温度的选取:
对于SSA,我们首先需要选定一个初始温度\(T_0\),初始温度选取的不同将直接影响到最终解的质量。如果初始温度太高,我们会有更高的概率可以获得到更高质量的解,但代价是算法需要花费更长的时间来计算并搜索解。因此,在实际应用的时候,可以先做几组实验看看效果如何,然后在效率和解质量这两者进行权衡。
这里介绍两种经验公式确定初温的方法:首先,随机产生一组状态\(S\)。然后,计算这组状态中两个状态间的最大差值(正数),即
最后,利用经验公式进行确定:
其中\(K\)应该是一个比较大的数。
退温策略:
这里介绍两种比较常用的降温策略:
(1)等差方式降温
其中\(k\)表示算法迭代的次数;\(\Delta T_0\)表示每次降温的大小。
(2)等比方式降温
理论上,温度降低的速度越慢,获得高质量解的概率越大,但如果温度降得太慢,反而影响算法的效率。一种折中的做法是温度高的时候下降得慢些,温度低时下降得快些。上面介绍的两种降温策略都是按照固定规则进行降温,可以在这基础之上进行改进。
新状态的产生:
如何通过初始状态产生新的状态也是算法的一种重点,应遵循的原则是所产生的候选解(状态)尽可能覆盖了全部的解空间(状态空间)。新状态的产生通常是在当前状态的邻域内,以一定的概率方式抽样出一个新的状态。可供选择的概率分布很多,如均匀分布、高斯分布、指数分布等等,需要自己选择。
需要注意的是,新状态的产生方式需要根据问题的性质来决定,并且如何定义当前状态的邻域也是一个问题。这里列举两个不同的示例:
(1)状态空间是离散的。
如果对于\(1,2,\cdots,n\)这\(n\)个数字的每种排列都分别对应着一个代价,我们的目标是找到使得代价最小的排列。新状态产生函数可以定义为随机交换两个数字的位置,例如
类似的可以定义为随机反转某个子序列等。
(2)状态空间是连续的。
如果初始状态是\(x_0\),那么我们可以从高斯分布中抽样出一个随机向量\(v\),然后得到新的状态\(x=x_0+v\)。
状态转移策略:
新状态产生之后,如果无条件接受,那么算法就变成了在状态空间随机”游走“,显然是不行的。因此,我们需要指定一个合适的策略使得算法能够向好的方向发展(即函数值减少的方向)。一个常用的准则就是Metropolis准则,公式表示如下:
其中\(T>0\)表示当前的温度;\(\Delta f=f(x_{new})-f(x_{old})\);\(p\)表示新状态被接受的概率。
实际上在编程的时候可以将Metropolis准则写成更紧凑的形式:
如果上述不等式成立,则接受新的状态,实际上这是一种轮盘赌的方式。
也就是如果产生的新状态比当前状态的能量更低(函数值更小)就发生状态转移,进入下一次迭代。如果新状态的能量更高,就以一定的概率接受更差的状态,而不是直接舍弃,这有助于跳出局部最优,发现更好的解。通过公式可知当\(T\)比较大的时候,更有可能接受较差的解。这意味着当温度较高时,算法倾向于在全局范围内随机搜索,搜索能力强。但当温度较低时,算法表现得较为保守,倾向于接受更好的解。
显然这里还有一个问题:如果产生了一个更差的状态,并且不被接受应该怎么办。具体的做法有以下几种:
- 不断循环抽样,直到产生的新状态被接受为止。
- 状态不发生改变。
- 设置一个抽样次数上限,如果达到次数上限后,新状态仍然未被接受,我们可以选择保留最后一次抽样的状态,或是不发生改变。
当然可能还有其他的做法。
停止条件:
- 设置算法的中止温度。
- 查看算法在连续多次迭代过程中,如果目标函数值已经不发生变化,或变化很小就停止迭代。
算法描述:
设置好初温\(t=t_0\),一个随机初始状态\(s=s_0\),迭代次数计数器\(k=0\)
while 终止条件未满足:
while True:
s_new = generate_new_state(s)
if 满足Metropolis准则:
s = s_new
break
t(k+1) = update(t(k))
k = k + 1
输出历史最优解
3.3 算法实现与测试
模拟退火算法具有一定的随机性,为了寻找到较好的解,最好限制搜索空间的大小。并且在运行过程中记录历史最优解。
import numpy as np
class SA:
"""
模拟退火算法:寻找函数的最小值
"""
def __init__(self, current_temp, gamma, final_temp, x0, st_x, func):
# 初始温度
self.current_temp = current_temp
# 温度采用线性递减策略
self.gamma = gamma
# 终止温度
self.final_temp = final_temp
# 空间维数
self.n = x0.size
# 约束x坐标的范围
self.st_x = st_x
# 目标函数
self.func = func
# 当前找到的解
self.solution = x0
# 当前解对应的函数值
self.value = func(x0)
# 历史最优解
self.best_solution = x0
self.best_value = func(x0)
def solve(self):
while self.current_temp > self.final_temp:
accept = False
offset = self.get_offset()
next_solution = self.solution + offset
next_solution[next_solution > self.st_x[1]] = self.st_x[1]
next_solution[next_solution < self.st_x[0]] = self.st_x[0]
if self.func(next_solution) < self.value:
accept = True
elif np.random.rand() < np.exp((self.value - self.func(next_solution)) / self.current_temp):
accept = True
if accept:
self.solution = next_solution
self.value = self.func(next_solution)
if self.value < self.best_value:
self.best_solution = self.solution
self.best_value = self.value
self.current_temp -= self.gamma
return self.best_solution, self.best_value
def get_offset(self):
# 产生随机偏移量
v = np.random.rand(self.n) - 0.5
v *= self.current_temp / 100
return v
def ackley(x):
return - 20 * np.exp(-0.2 * np.sqrt((x * x).mean())) - np.exp(np.cos(2 * np.pi * x).mean()) + 20 + np.exp(1)
if __name__ == '__main__':
# ackley函数的最小值
for i in range(10):
sa = SA(
current_temp=500,
gamma=0.1,
final_temp=1,
st_x=np.array([-5, 5]),
x0=np.random.randn(2),
func=ackley
)
solution, value = sa.solve()
print(solution, value)
测试结果:
![image-20200305220547484](https://img2020.cnblogs.com/blog/1530546/202004/1530546-20200423224513359-1661103593.png)