最近在学习玻尔兹曼机，里面用到了模拟退火算法，经过一天的实验，总算顺利完成，本文打算记录这一过程，以作备忘。
本文内容如下：
1、实验环境
2、算法原理简介
3、TSP案例代码实现
4、运行结果解析
5、算法的优缺点分析
6、关于模拟退火算法改进的一些想法
一、实验环境
编程语言：R(R语言中求解模拟退火算法的函数有optim，sann等，其中sann在ConsPlan包中，本文基于模拟退火算法原理自行编写)
使用数据：中国31个城市的距离数据，详见附件
其他：无
二、算法原理简介
模拟退火算法是一种随机搜索算法，用来求解表面粗糙(含有非常多的局部最小点)，无法微分求导的最优值的问题，由于它结合了一定的贪心策略和启发策略，因此比纯粹的随机搜索算法更容易找到最优解，它的基本思想是模拟金属退火的过程，大家可以想象要将两种金属融合到一起的(比如炼钢，需要将铁和其他微量金属：铜、锰等融合在一起)，首先需要先将金属加热至高温融化，使其原子处于高速运动状态，此时各种金属原子会做布朗运动进行互相融合，然后在缓慢降温，最后达到稳定状态(收敛)，内能最低(最优解)，熵值最大。
这样解释，也许大家还是不太明白，我们回到函数求解场景上，假设有一个函数(非凸函数)，它的最优解无法用常规的微分求导(比如梯度下降、牛顿法等)来求出，且它的函数值域上有很多局部最小点，这时候有什么办法来尽量避免求解到局部最小点，而找到一个最优解或者近似最优解呢，我们想到了一个办法，随机丢一个小球到函数的值域空间上(即随机取一个初始值)，然后让小球随机晃动，初始的时候温度较高，小球会晃动剧烈，这时小球会更有机会掉入周围更低的坑(函数值更小)中，从而避免掉入一个局部的小坑中，然后在逐步降温，小球晃动减少，最终收敛到一个比较理想的坑中(坑比较深，因此晃动不出来)，从而找到最优解或者近似最优解。
从上述原理可以看出，模拟退火算法的关键在于温度的选取、降温过程的控制以及小球该如何做晃动运动(更新迭代策略)，下面总结一下我的个人经验
1.1 更新迭代策略
前面我们说过模拟退火算法是一种随机搜索算法，它的随机性体现更新迭代过程中，即在小球当前位置做随机扰动，如果小球晃动到一个比之前更低的点上(值比之前的小)，则将当前位置更新为最新晃动位置，如果小球晃动到一个比之前更高的点上(值比之前的大)，则以一定的概率更新为最新晃动位置(这个概率服从玻尔兹曼分布)，否则回到之前的点上继续晃动。
从上述描述可以看出，更新迭代策略具有一定贪心和启发性在里面。
1.2
上面更新策略提到一个概率分布，即玻尔兹曼分布，它的概率分布公式简化形式如下：P=exp(-(f(x1)-f(x0))/T) 其中f(x1)为晃动后的值，f(x0)为晃动前的值，T即为温度
> exp(-10/10)
[1] 0.3678794
> exp(-10/100)
[1] 0.9048374
> exp(-10/5)
[1] 0.1353353
> exp(-10/1)
[1] 4.539993e-05
> exp(-1/10)
[1] 0.9048374
> exp(-5/10)
[1] 0.6065307
> exp(-100/100)
[1] 0.3678794
> exp(-10/90)
[1] 0.8948393
> exp(-10/1)
[1] 4.539993e-05
>
从上可以看出，在差值不变的情况下，温度越高，更新迭代概率越大，在温度一定的情况下，差值绝对值越大，概率越小，这正是我们所期望的，因此关于温度的选择应该与差值的分布保持一致，温度过低会导致晃动能力太差，更容易收敛于局部最小点，温度过大，如果温度下降不快，则不容易收敛，始终处于随机跳动状态。
1.3 降温控制
目前有三种降温策略(指数exp,对数log,倒数rec)，详情如下：
很明显，指数法降温过程整体比较平稳，对数法前期降温比较明显，但后期下降非常缓慢，倒数法则降温迅速，至于选取哪种方式与具体的场景有关，比如古代炼剑，不同的淬火方式，对剑的任性影响很大，需要根据用途，不断的尝试才能找到合适的。
其他有关模拟退火的详细介绍，请读者自行百度，知乎上有很多大牛的解释很精辟。
三、TSP案例代码实现
1、TSP问题描述
本文需要求解中国31个城市的TSP问题，代码具体如下：
2、原始数据读取tsp
tsp
label
3、转成距离矩阵##城市距离矩阵
m_dis
dis
{
n
m_dis
colnames(m_dis)
rownames(m_dis)
k
for(i in 1:(n-1))
{
m_dis[i,(i+1):n]
k
}
for(i in 1:n)
{
for(j in 1:n)
{
if(j
{
m_dis[i,j]
}
}
}
return(m_dis);
}
4、随机产生新的路径genseq
n
if(runif(1,0,1)<0.25)
{
index
tmp
sq[index[1]]
sq[index[2]]
}else
{
index
old_sq
sq[index[1]:index[3]]
}
return(sq)
}
5、模拟退火算法主程序
流程图：
sann
{
require(stringr);
##加载计算距离函数
source("distance.txt");
##加载随机产生一个路径函数
source("genseq.txt");
##初始路径
if(is.null(sq))
{
sq
}
##初始距离
jl
##优化过程记录
log
log$sq
log$dis[1]
log$t[1]
log$sq[1]
k
t0
##外循环：降温
for(i in 1:m_iter)
{
##内循环：跳坑
for(j in 1:t_iter)
{
##基于之前的路径随机变动两个节点生成下一条路径
tmp_sq
##计算该路径距离
tmp_jl
if(tmp_jl
{
jl
sq
}else
{
if(exp((jl-tmp_jl)/t)>runif(1,0,1))
{
jl
sq
}
}
log$dis[k+1]
log$t[k+1]
log$sq[k+1]
k
}
##降温
if(method=="exp")
{
t
}else if(method=="log")
{
t
}else if(method=="rec")
{
t
}else
{
t
}
}
return(list(sq=sq,jl=jl,log=log));
}
四、运行结果解析
4.1 运行代码T
Method
Iter
count
res
res$M
k
for(m in Method)
{
for(t in T)
{
for(iter in Iter)
{
for(i in 1:count)
{
x
res$M[k]
res$T[k]
res$Iter[k]
res$max_iter[k]
res$min_jl[k]
best_iter
res$best_iter[k]
res$best_t[k]
k
}
}
}
}
4.2 运行结果> res
idx M T Iter max_iter min_jl best_iter best_t
1 1 exp 500 10 10000 16593.18 8577 2.303278e-09
2 2 exp 500 10 10000 16956.54 6910 3.727921e-07
3 3 exp 500 10 10000 15604.95 9665 8.326554e-11
4 4 exp 500 10 10000 15456.42 3927 3.262044e-03
5 5 exp 500 10 10000 16088.69 8298 5.404278e-09
6 6 exp 500 10 10000 15588.28 4628 3.868221e-04
7 7 exp 500 10 10000 16550.58 9725 6.935786e-11
8 8 exp 500 10 10000 15551.24 8575 2.303278e-09
9 9 exp 500 10 10000 16314.27 9966 3.339007e-11
10 10 exp 500 10 10000 15820.15 2915 7.072039e-02
11 11 exp 500 30 30000 16405.02 4953 3.283267e+00
12 12 exp 500 30 30000 15966.83 17112 1.441646e-05
13 13 exp 500 30 30000 15569.39 10276 1.495902e-02
14 14 exp 500 30 30000 15426.42 6131 1.000930e+00
15 15 exp 500 30 30000 15569.39 9155 4.616884e-02
16 16 exp 500 30 30000 16682.08 8669 7.748711e-02
17 17 exp 500 30 30000 15466.42 22742 4.698341e-08
18 18 exp 500 30 30000 16755.78 9383 3.730363e-02
19 19 exp 500 30 30000 15910.90 9304 3.964675e-02
20 20 exp 500 30 30000 15939.42 22075 9.466614e-08
21 21 exp 500 50 50000 15499.76 24318 1.862228e-04
22 22 exp 500 50 50000 15782.74 7857 4.189211e+00
23 23 exp 500 50 50000 15520.87 11644 4.265917e-01
24 24 exp 500 50 50000 15758.67 13757 1.151322e-01
25 25 exp 500 50 50000 16375.75 31689 2.115816e-06
26 26 exp 500 50 50000 16069.80 7797 4.452345e+00
27 27 exp 500 50 50000 15569.39 7030 7.030930e+00
...详见附件res.txt
4.3 结果分析
4.3.1 最佳距离
从上两图来看，log算法表现不如exp和rec，主要原因是其后期温度下降缓慢，如果外循环迭代次数不够，无法收敛，而且温度越高，收敛越困难，另外从整体获得最优解看，初始温度高一些更有机会获得最优解，但是所需的计算迭代代价也很高，在合理的计算成本下，选择合适的温度加上适当的内循环探索次数也能达到相同的效果，而且所需的计算成本也不大。
4.3.2 最佳终止温度
从上图结果来看，本案例需要温度降低到0度左右，才能收敛到最佳距离，而且最佳距离与终止温度相关性较大，这也解释了为何log算法表现不好的原因。
4.3.3 最佳循环比
从上图可以看出，内循环探索次数越高，达到最佳循环比更低，因此设置合适的内循环次数(小球在周边自由晃动的次数)对于找到最优解很有帮助。
五、算法的优缺点分析
通过上述的原理介绍和运行结果，大致可以总结出模拟退火算法的一些优缺点
5.1 优点：有大概率在复杂、非凸、表面粗糙的问题中获得最优解，或者说即使得不到最优解，也能获得一个不错的近似解
5.2 缺点：应用算法需要对求解问题非常熟悉，才能选择到合适的参数，同时算法迭代次数较多，对于大数据量的场景需要慎重考虑
六、关于模拟退火算法改进的一些想法
通过上述实验，发现迭代更新策略对算法的影响很大，也就是说小球如何晃动，本文实现的代码非常随机，可以说毫无方向性可言，也许加入遗传算法的思想，让小球向着进化的方向晃动，也许能够更快收敛，从而提升算法的运行效率，后续有时间在试一试这个想法。
625
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
