模拟退火算法

模拟退火算法——脱胎于自然界物理过程，与优化问题相结合

模拟退火算法

工程中许多实际优化问题都是非凸的，存在许多局部最优解，特别是随着优化规模的增大，局部最优解的数目将会迅速增加。因此，有效地求一般的随机优化问题的解至今仍是一个难题。目前求解一全局最优化问题的方法一般分为两类，一类是确定性方法，另一类是随机性方法。

模拟退火算法是一种通用的概率算法，用来在一个很大的空间内寻找问题的最优解。由于其能有效解决NP问题、避免陷入局部最优、初值依赖性较弱等特点，已在控制工程、生产调度、机器学习、图像处理等领域获得了广泛的应用。

算法基本思想

现代的模拟退火算法思想源于固体退火过程，即将固体热至足够高的温度，再缓慢冷却；升温的过程中，固体内部的粒子随温度升高变为无序状，内能增大，而缓慢冷却时粒子又逐渐趋于有序。理论上，如果冷却的速度足够慢，那么冷却中任一温度时固体都能达到热平衡，而冷却到低温时，将达到这一低温状态的内能最小状态。

在这一过程中，任一恒定温度都能达到热平衡是个很关键的步骤，这一点可以用Monte Carlo算法模拟。但因为物理系统总是趋向于能量最低，而分子热运动则趋向于破坏这种低能量的状态，故而只需要着重取贡献比较大的状态即可达到比较好的效果。

设当前状态i生成新状态j，若新状态的内能小于状态i的内能（即Ej<Ei），则接受新状态j作为新的当前状态；否则，以概率exp[ -( Ei-Ej )/kT ]接受状态j，其中k为Boltzmann常数，这就是通常所说的Metropolis准则。

后来，人们把模拟退火算法思想与组合优化的相似点进行类比，将模拟退火算法应用到了组合优化领域。在模拟退火算法应用的时候，一般可以将温度T当做控制参数，目标函数值f视为内能E，而固体在某温度T时的一个状态对应一个解xi。然后算法试图随着控制参数的降低是目标函数值也逐渐降低，直至区域全局最小值。

冷却进度表

1. 控制参数的初值T0：冷却开始的温度
2. 控制参数T的衰减函数：由于计算机只能处理离散的数据，因此需要把连续的降温过程离散化为一系列温度点，衰减函数即计算一系列温度的表达式
3. 控制参数T的终值Tf ：冷却结束的温度
4. Markov链的长度Lk：任一温度T的迭代次数

算法基本步骤

1、令T=T0，随机生成一个初始解x0，并计算其的目标函数的值E(x0)

2、令T等于冷却进度表中的下一个值Ti

3、根据当前解xi进行扰动，产生一个新的解xj，计算相应的目标函数值E(xj)，得到δE=E(xj)-E(xi)

4、若δE＜0，则接受新解xj作为当前解；若δE＞0，则新解xj按概率exp(-δ/Ti)接受，Ti为当前温度

5、在温度Ti下，重复Lk次的扰动和接受过程，即执行步骤(3)和步骤(4)

6、判断T是否达到Tf，若达到，则终止算法；若未达到，转到步骤2继续执行

算法的实质是两重循环，在任一温度下随机扰动产生新解，并计算目标函数值的变化，决定是否被接受。由于算法的初始温度比较高，这样，使E增大的新解在初始时也能被接受，因而能跳出局部最小值，然后通过缓慢地降低温度，算法就最终可能收敛到全局的最优解。因为低温时虽然接收函数很小，但也有可能碰到更差的解，因此一般会把历史最优解记录下来，与终止算法前最后一个接受解一起输出

模拟退火算法的python实现（解决背包问题）

import numpy as np
from math import *

Inf=float('inf') #定义正无穷

a=0.95 #衰减函数的参数
k=np.array([5,10,13,4,3,11,13,10,8,16,7,4]) #物品的价值
k=-k #模拟退火算法求的是最小值，故取负数
d=np.array([2,5,18,3,2,5,10,4,11,7,14,6]) #物品的质量
restriction=46 #背包的最大容量
num=d.shape[0] #物品的数量

"""
    E_new是新解的目标函数值
    E_current是当前解对应的目标函数值
    E_best是最优解
    初始温度T0为97
    终止温度Tf为3
    控制温度为T
"""
sol_new=np.ones((1,num)) 
E_current ,E_best = Inf ,Inf
sol_current ,sol_best = sol_new ,sol_new
T0 ,Tf = 97 ,3
T = T0
p = True

while T>=Tf:
    for i in range(0,100):
        #产生随机扰动
        tmp = round( np.random.rand()*num - 1)
        sol_new[0,tmp] = 1-sol_new[0,tmp]
        
        #检查是否满足约束
        while True:
            q = np.sum(sol_new*d) <= restriction
    
            if not q:
                p = not p
                tmp = np.where(sol_new == 1 )
                if p : sol_new[sol_new ==1 ] = 0
                else : sol_new[0,tmp[-1][-1]] = 0
            
            else : break

        
        #计算背包中的物体价值
        E_new = np.sum(sol_new*k)
        if E_new < E_current :
            E_current = E_new
            sol_current = sol_new
            if E_new < E_best:
                E_best = E_new
        else:
            if np.random.rand()< exp(-(E_new - E_current)/T ):
                E_current = E_new
                sol_current = sol_new
            else:
                sol_new = sol_current
        
    T = T*a

print(sol_best)#最优解
print(-E_best)#最优总价值

欢迎关注公众号 ： 数学算法实验室