模拟退火算法

一、固体退火原理
在热力学上，退火（annealing）现象指物体逐渐降温的物理现象，温度愈低，物体的能量状态会低；够低后，液体开始冷凝与结晶，在结晶状态时，系统的能量状态最低。大自然在缓慢降温（亦即，退火）时，可“找到”最低能量状态：结晶。但是，如果过程过急过快，快速降温（亦称「淬炼」，quenching）时，会导致不是最低能态的非晶形。

二、模拟退火
首先看一副图，假设有一个函数如下图，为了求他的最大值，可以从A点开始向右试探，如果获得值比他大则试探过程继续，当到达B点时试探结束，因为此时向左向右试探函数值都是在变小。

而模拟退火在搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。以上图为例，模拟退火算法在搜索到局部最优解B后，会以一定的概率接受向右继续移动。也许经过几次这样的不是局部最优的移动后会到达C，此时继续试探就可以达到最优解D，于是就跳出了局部最小值B。
根据Metropolis准则，粒子在温度T时趋于平衡的概率为exp(-ΔE/(kT))，其中E为温度T时的内能，ΔE为其改变数,k为Boltzmann常数。Metropolis准则常表示为
Metropolis准则表明，在温度为T时，出现能量差为dE的降温的概率为P(dE)，表示为：P(dE) = exp( dE/(kT) )。其中k是一个常数，exp表示自然指数，且dE<0。所以P和T正相关。这条公式就表示：温度越高，出现一次能量差为dE的降温的概率就越大；温度越低，则出现降温的概率就越小。又由于dE总是小于0（因为退火的过程是温度逐渐下降的过程），因此dE/kT < 0 ，所以P(dE)的函数取值范围是(0,1) 。随着温度T的降低，P(dE)会逐渐降低。

我们将一次向较差解的移动看做一次温度跳变过程，我们以概率P(dE)来接受这样的移动。也就是说，在用固体退火模拟组合优化问题，将内能E模拟为目标函数值 f，温度T演化成控制参数 t，即得到解组合优化问题的模拟退火演算法：由初始解 i 和控制参数初值 t 开始，对当前解重复“产生新解→计算目标函数差→接受或丢弃”的迭代，并逐步衰减 t 值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值 t 及其衰减因子Δt 、每个 t 值时的迭代次数L和停止条件S。
三、算法流程图

四、算法关键参数

1. 状态产生函数
   设计的出发点：尽可能保证产生的候选解遍布全部解空间。
   最常用的状态产生函数：
   为扰动幅度参数，为随机扰动变量。
   随机扰动可服从柯西分布、高斯分布、均匀分布。

2. 状态接受函数

是SA实现全局搜索的关键因素，一般以概率方式给出。
最常用的状态接受函数：
          min{1, exp[-(C(sj)-C(si))/tk]}>=random[0,1]?
设计状态接受函数，应该遵循以下原则： 

 1. 固定温度下，接受使目标函数值下降的候选解的概率要大于使目标函数值上升的候选解的概率；
 2. 随温度下降，接受使目标函数值上升解的概率要逐渐减小；
 3. 当温度趋于零时，只能接受目标函数值下降的解

3. 初始温度

实验表明，初温t0越大，获得高质量解的几率越大，但计算时间增加。 
初温的确定应折衷考虑优化质量和效率。
常用的确定初温的方法包括：  

 1. 均匀抽样一组状态，以各状态目标值的方差为初温。
 2. 随机产生一组状态，确定两两状态间的最大目标值差|∆max|，计算初温，t0=-∆max / lnpr 初始接受概率pr理论上接近于1。
 3. 设为某个较大的常数。

4. 温度更新函数

用于在外循环中修改温度值。

最常用的是指数退温函数：
				tk+1=atk
为退温速率，0<a<l，a大小可以不断变化。 

5. 内循环终止准则

 - 或称Metropolis抽样稳定准则
 
 - 常用的抽样稳定准则包括：
      检验目标函数的均值是否稳定；
      连续若干步的目标值变化较小；
      按一定的步数抽样。

6. 外循环终止准则

 - 用于决定算法何时结束。
 - 设置温度终值te是一种简单的方法，SA算法的收敛性理论中要求te趋于零，这显然不实际。
 - 通常的做法包括：

    设置终止温度的阈值；
    设置外循环迭代次数；
    算法搜索到的最优值连续若干步保持不变；  
    检验系统熵是否稳定。

五、伪代码

/*
* J(y)：在状态y时的评价函数值
* Y(i)：表示当前状态
* Y(i+1)：表示新的状态
* r： 用于控制降温的快慢
* T： 系统的温度，系统初始应该要处于一个高温的状态
* T_min ：温度的下限，若温度T达到T_min，则停止搜索
*/
while( T > T_min )
{
　　dE = J( Y(i+1) ) - J( Y(i) ) ; 

　　if ( dE >=0 ) //表达移动后得到更优解，则总是接受移动
       Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
　　else
　　{
       // 函数exp( dE/T )的取值范围是(0,1) ，dE/T越大，则exp( dE/T )也
      if ( exp( -dE/T ) > random( 0 , 1 ) )
      Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
　　}
　　T = r * T ; //降温退火 ，0<r<1 。r越大，降温越慢；r越小，降温越快
　　/*
　　* 若r过大，则搜索到全局最优解的可能会较高，但搜索的过程也就较长。若r过小，则搜索的过程会很快，但最终可能会达到一个局部最优值
　　*/
　　i ++ ;
}

六、使用模拟退火算法解决旅行商问题

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(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"

代码出处：https://www.cnblogs.com/sench/p/9427193.html