最近在学习玻尔兹曼机，里面用到了模拟退火算法，经过一天的实验，总算顺利完成，本文打算记录这一过程，以作备忘。
本文内容如下：
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<-scan("tsp.txt",sep=",") ##读取城市距离数据
tsp<-na.omit(tsp) ##去掉空值
label<-c("北京","上海","天津","石家庄","太原","呼和浩特","沈阳","长春","哈尔滨","济南","南京","合肥","杭州","南昌","福州","台北","郑州","武汉","长沙","广州","海口","南宁","西安","银川","兰州","西宁","乌鲁木齐","成都","贵阳","昆明","拉萨") ##加上城市标签
3、转成距离矩阵
##城市距离矩阵
m_dis<-dis(tsp,label)
dis<-function(tsp,label) ##还原城市距离矩阵函数
{
n<-length(label);
m_dis<-matrix(0,n,n);
colnames(m_dis)<-label;
rownames(m_dis)<-label;
k<-1;
for(i in 1:(n-1))
{
m_dis[i,(i+1):n]<-tsp[k:(k+n-i-1)];
k<-k+n-i;
}
for(i in 1:n)
{
for(j in 1:n)
{
if(j<i)
{
m_dis[i,j]<-m_dis[j,i]
}
}
}
return(m_dis);
}
4、随机产生新的路径
genseq <- function(sq) { ##随机产生一个路径函数
n<-length(sq);
if(runif(1,0,1)<0.25)
{
index <- sample(2:(n-1), size=2, replace=FALSE);
tmp <- sq[index[1]];
sq[index[1]] <- sq[index[2]];
sq[index[2]] <- tmp;
}else
{
index <- sort(sample(2:(n-1), size=3, replace=FALSE));
old_sq<-sq;
sq[index[1]:index[3]]<-c(old_sq[index[2]:index[3]],old_sq[index[1]:(index[2]-1)]);
}
return(sq)
}
5、模拟退火算法主程序
流程图：
sann<-function(m_dis,sq,m_iter=100,t_iter=50,t=100000,method="exp",d_rate=0.97)
{
require(stringr);
##加载计算距离函数
source("distance.txt");
##加载随机产生一个路径函数
source("genseq.txt");
##初始路径
if(is.null(sq))
{
sq <- c(1,2:NROW(m_dis),1);
}
##初始距离
jl<-distance(m_dis,sq);
##优化过程记录
log<-data.frame(iter=0:(m_iter*t_iter),dis=0,t=0,sq="");
log$sq<-as.character(log$sq);
log$dis[1]<-jl;
log$t[1]<-t;
log$sq[1]<-str_c(sq,collapse=",");
k<-1;##循环次数
t0<-t;
##外循环：降温
for(i in 1:m_iter)
{
##内循环：跳坑
for(j in 1:t_iter)
{
##基于之前的路径随机变动两个节点生成下一条路径
tmp_sq<-genseq(sq);
##计算该路径距离
tmp_jl<-distance(m_dis,tmp_sq);
if(tmp_jl<jl)
{
jl<-tmp_jl;
sq<-tmp_sq;
}else
{
if(exp((jl-tmp_jl)/t)>runif(1,0,1))
{
jl<-tmp_jl;
sq<-tmp_sq;
}
}
log$dis[k+1]<-jl;
log$t[k+1]<-t;
log$sq[k+1]<-str_c(sq,collapse=",");
k<-k+1;
}
##降温
if(method=="exp")
{
t<-t*d_rate;
}else if(method=="log")
{
t<-t0/log(i+1);
}else if(method=="rec")
{
t<-t0/(i+1);
}else
{
t<-t*d_rate;
}
}
return(list(sq=sq,jl=jl,log=log));
}
四、运行结果解析
4.1 运行代码
T<-c(500,1000,5000) ##初始温度
Method<-c("exp","log","rec") ##降温方式
Iter<-c(10,30,50) ##随机跳坑次数
count<-10 ##为了降低随机性影响，评估算法稳定性，每类参数运行的次数
res<-data.frame(idx=1:(27*count),M="",T=0,Iter=0,max_iter=0,min_jl=0,best_iter=0,best_t=0);
res$M<-as.character(res$M);
k<-1;##循环次数
for(m in Method)
{
for(t in T)
{
for(iter in Iter)
{
for(i in 1:count)
{
x<-sann(m_dis=m_dis,sq=sq,m_iter=1000,t_iter=iter,t=t,method=m,d_rate=0.97);
res$M[k]<-m;
res$T[k]<-t;
res$Iter[k]<-iter;
res$max_iter[k]<-nrow(x$log)-1;
res$min_jl[k]<-x$jl;
best_iter<-min(which(x$log$dis==x$jl));
res$best_iter[k]<-best_iter;
res$best_t[k]<-x$log$t[best_iter];
k<-k+1;
}
}
}
}
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.txt4.3 结果分析
4.3.1 最佳距离
从上两图来看，log算法表现不如exp和rec，主要原因是其后期温度下降缓慢，如果外循环迭代次数不够，无法收敛，而且温度越高，收敛越困难，另外从整体获得最优解看，初始温度高一些更有机会获得最优解，但是所需的计算迭代代价也很高，在合理的计算成本下，选择合适的温度加上适当的内循环探索次数也能达到相同的效果，而且所需的计算成本也不大。
4.3.2 最佳终止温度
从上图结果来看，本案例需要温度降低到0度左右，才能收敛到最佳距离，而且最佳距离与终止温度相关性较大，这也解释了为何log算法表现不好的原因。
4.3.3 最佳循环比
从上图可以看出，内循环探索次数越高，达到最佳循环比更低，因此设置合适的内循环次数（小球在周边自由晃动的次数）对于找到最优解很有帮助。
五、算法的优缺点分析
通过上述的原理介绍和运行结果，大致可以总结出模拟退火算法的一些优缺点
5.1 优点：有大概率在复杂、非凸、表面粗糙的问题中获得最优解，或者说即使得不到最优解，也能获得一个不错的近似解
5.2 缺点：应用算法需要对求解问题非常熟悉，才能选择到合适的参数，同时算法迭代次数较多，对于大数据量的场景需要慎重考虑
六、关于模拟退火算法改进的一些想法
通过上述实验，发现迭代更新策略对算法的影响很大，也就是说小球如何晃动，本文实现的代码非常随机，可以说毫无方向性可言，也许加入遗传算法的思想，让小球向着进化的方向晃动，也许能够更快收敛，从而提升算法的运行效率，后续有时间在试一试这个想法。
