[学习笔记] 模拟退火 (Simulated Annealing)

真没想到这东西真的在考场上用到了...顺便水篇blog以示诈尸好了(逃

模拟退火算法

模拟退火是一种随机化算法, 用于求函数的极值qwq

好比给出一个问题, 咱们要求最优解的值, 可是可能的方案数量极大, 直接搜索会T飞(或者方案是连续的总数无穷根本无法搜), 这种时候咱们通常会有两种选择:

登山算法

登山算法每次在当前找到的方案附近寻找一个新的方案(常见方式是随机一个差值), 而后若是这个解更优那么直接转移.

对于单峰函数来讲这显然能够直接找到最优解(不过你都知道它是单峰函数了为啥不三分呢?)

可是对于多数咱们求解的函数来讲, 它并不必定会长成这个样子...因而就极其有可能钻进一个局部的最优解出不来了

算法得出的最优解与初始解的位置以及搜寻的附近解的区域大小有关.

固然若是你寻找新方案的区间很大的话有几率跳出去, 可是太大的话又可能跳来跳去跳乱了从而找不到最优解...

欧皇专用最优化求解方式(@liu_runda)

 

然而并非全部人都是欧皇, 像博主这样的非酋要怎么办捏?

固然是求助于天然规律(大雾)

退火的理论部分

退火其实原本是冶金工业里的术语...大概过程是先把晶体加热到极高的温度再缓慢降温, 在这个过程当中减小晶体中的缺陷(达到能量最低的最稳定状态)

具体理论是这样的:

设 $E[\{x_i\}]$ 表示某一物质体系在微观状态 $\{x_i\}$ 下的内能, 对于给定温度 $T$, 若体系处于热平衡态时, $E[\{x_i\}]$ 服从 Boltzmann 分布:

$$ f=c(T)\exp\left(-\frac{E[\{x_i\}]}{kT}\right) $$

其中 $k$ 为 Boltzmann 常数.

$T$ 降低, $E$ 随之降低. 若 $T$ 降低的速度足够慢, 则体系总能够保持热平衡态. 当 $T\to 0$ 时 $E$ 趋近于最小. 这样的物质降温过程被称为退火过程.

而后机智的咱们发现这个过程最终和咱们的最优化过程相似!

因而咱们去模拟这个过程, 按照退火的规律引入更多随机因素, 这样咱们获得最优解的几率就会增长辣233

能够证实, 模拟退火获得的最终解依几率收敛于最优解.

emmmm...

等等模拟这个过程? 这是计算机又不是实验室你怎么模拟啊(╯°□°)╯︵┻━┻

拿出物理化学(伪装本身仍是个ChO党)...

根据热力学规律并结合计算机对离散数据的处理, 咱们定义: 若是当前温度为 $T$ , 当前状态与新状态之间的能量差为 $\Delta E$ , 则发生状态转移的几率为:

$$ P(\Delta E) = e^{\frac { \Delta E } { kT } } $$

显然若是 $ \Delta E$ 为正的话转移是必定会成功的, 可是对于 $\Delta E < 0$ 咱们则以上式中计算获得的几率接受这个新解.

而后咱们维护温度 $T$ 便可. 这里咱们有三个参数: 初温 $T_0$ , 降温系数 $d$ , 终温 $T_k$

通常 $T_0$ 是个比较大的数, $d$ 是个接近 $1$ 可是小于 $1$ 的值, $T_k$ 是个接近 $0$ 的正值.

首先让温度 $T=T_0$ , 而后进行一次转移尝试, 而后让 $T=dT$.

当 $T<T_k$ 时模拟退火过程结束, 当前解做为最优解.

看起来好像并非很难理解?

Wikipedia上的动图:

通常来说模拟退火在参数合适的状况下效果拔群, TSP随便跑(x

模拟退火的实际使用

实际使用里这函数可不必定是个单元函数...并且寻找新解好像是个很模糊的东西, 毕竟不少时候咱们会发现咱们要求解的问题的全部可能解并非离散的...

先拿道题说说...

3680: 吊打XXX

Time Limit: 10 Sec  Memory Limit: 128 MBSec  Special Judge
Submit: 4247  Solved: 1566
[Submit][Status][Discuss]

Description

gty又虐了一场比赛，被虐的蒟蒻们决定吊打gty。gty见大势很差机智的分出了n个分身，但仍是被人多势众的蒟蒻抓住了。蒟蒻们将
n个gty吊在n根绳子上，每根绳子穿过天台的一个洞。这n根绳子有一个公共的绳结x。吊好gty后蒟蒻们发现因为每一个gty重力不一样，绳
结x在移动。蒟蒻wangxz脑洞大开的决定计算出x最后停留处的坐标，因为他太弱了决定向你求助。
不计摩擦，不计能量损失，因为gty足够矮因此不会掉到地上。

Input

输入第一行为一个正整数n(1<=n<=10000)，表示gty的数目。
接下来n行,每行三个整数xi，yi，wi，表示第i个gty的横坐标，纵坐标和重力。
对于20%的数据，gty排列成一条直线。
对于50%的数据，1<=n<=1000。
对于100%的数据，1<=n<=10000,-100000<=xi,yi<=100000

Output

输出1行两个浮点数（保留到小数点后3位），表示最终x的横、纵坐标。

Sample Input

3
0 0 1
0 2 1
1 1 1

Sample Output

0.577 1.000

HINT

 

Source

By wangxz

在这个题中咱们为了到达稳定状态要让整个体系的总重力势能最低.

重力势能怎么求呢? 别忘了这绳子的总长度是不会变的...因而某个质点的重力势能和到绳结的水平距离成一次函数关系.

咱们为了简化问题, 能够将某个质点对最终答案产生的贡献计算为 $ dis(o,x)\times m_x $ . 而后咱们要让这个值最小化.

这个时候咱们能够考虑模拟退火. 首先随机一个点做为初始解(为了加速收敛, 咱们能够直接取各个点坐标的平均值所在的店). 而后随机两个值做为差值加到这个点的坐标上做为下一个解.

而后模拟退火直接往上套就能够了233

具体实现就是一个  while  循环, 循环内有4步:

1. 根据当前解找到下一个解
2. 计算下一个解的 "能量" (也就是价值)
3. 决定是否要接受这个新解
4. 降温

找下一个解的时候有一个提升精度的小技巧: 根据当前温度决定差值的范围. 这样在降温即将结束接近最优解的时候能够有更大的几率更精确地命中最优解.

固然若是是解是离散的就不能这样搞了. 以及生成下一个解的时候万万不能所有从新随机生成, 那就和瞎随没区别了...要随机做出一些相对小的修改.

具体作法就是使用一个产生 $[0,1]$ 随机实数的函数, 将随机区间转为 $[-1,1]$ 后乘上 $T$ 做为差值. (也就是生成一个 $[-T,T]$ 的随机值做为差值)

不过实际操做的时候咱们较少直接输出最终解, 而是选择在模拟退火的过程当中单独维护一个解, 只在遇到更优解的时候将其更新, 增长正确率.

另外一个提升正确率的方法是: 屡次进行模拟退火过程(或者说"从新烧热再退火一遍"), 每次取最优解.

还有就是最后烧完以后能够再在全局最优解的基础上进行登山.

本题的参考实现:

 1 #include <bits/stdc++.h>
 2  
 3 const int MAXN=1e4+10;
 4  
 5 struct Point{
 6     double x;
 7     double y;
 8     Point(double x=0,double y=0){
 9         this->x=x;
10         this->y=y;
11     }
12 };
13 Point P[MAXN];
14  
15 int n;
16 Point ans;
17 int g[MAXN];
18 double minAns=DBL_MAX;
19  
20 double Rand();
21 double Sqr(double);
22 double Calc(const Point&);
23 bool Accept(double,double);
24 double EucDis(const Point&,const Point&);
25 Point SimulatedAnnealing(Point,double,double,double);
26  
27 int main(){
28     scanf("%d",&n);
29     Point init;
30     for(int i=0;i<n;i++){
31         scanf("%lf%lf%d",&P[i].x,&P[i].y,g+i);
32         init.x+=P[i].x;
33         init.y+=P[i].y;
34     }
35     init.x/=n;
36     init.y/=n;
37     SimulatedAnnealing(init,1e5,1-7e-3,1e-3);
38     printf("%.3f %.3f\n",ans.x,ans.y);
39     return 0;
40 }
41  
42 inline double Calc(const Point& origin){
43     double ans=0;
44     for(int i=0;i<n;i++){
45         ans+=EucDis(origin,P[i])*g[i];
46     }
47     if(ans<minAns){
48         ::ans=origin;
49         minAns=ans;
50     }
51     return ans;
52 }
53  
54 bool Accept(double delta,double tmp){
55     return delta<0||Rand()<exp(-delta/tmp);
56 }
57  
58 Point SimulatedAnnealing(Point initAns,double initT,double dec,double end){
59     double tmp=initT;
60     Point now=initAns;
61     double nowAns=Calc(now);
62     while(tmp>end){
63         Point next=Point(now.x+tmp*(Rand()*2-1),now.y+tmp*(Rand()*2-1));
64         double ans=Calc(next);
65         if(Accept(ans-nowAns,tmp)){
66             now=next;
67             nowAns=ans;
68         }
69         tmp*=dec;
70     }
71     for(int i=0;i<1000;i++){
72         Point rnd=Point(ans.x+tmp*(Rand()*2-1),ans.y+tmp*(Rand()*2-1));
73         Calc(rnd);
74     }
75     return now;
76 }
77  
78 inline double Rand(){
79     return double(rand())/double(RAND_MAX);
80 }
81  
82 inline double Sqr(double x){
83     return x*x;
84 }
85  
86 inline double EucDis(const Point& a,const Point& b){
87     return sqrt(Sqr(a.x-b.x)+Sqr(a.y-b.y));
88 }

View Code

可是调参的过程仍是比较看脸的...究极非洲大酋长慎用(x

通常状况下咱们会在时间容许的状况下尽可能多地尝试新的解. 通常降温系数 $d$ 与 $1$ 的差减小一个数量级, 时间消耗大约多 $10$ 倍, $T_0$ 和 $T_k$ 变化一个数量级, 时间消耗不会变化很大.

这种时候咱们能够试着先本机跑跑自造数据看看精度怎么样. 若是发现常常陷入局部最优解的话考虑增大 $T_0$ 和 $d$ , 若是发现最终精度不够的话考虑减少 $T_k$.

至于模拟退火的正确率计算么...好像只有实验是最方便的了吧(x

今天上午考试的时候手调一波参数而后极限数据下测试 $100$ 次发现精度达标率有 $60\%$ 就交了...而后A了...

因而借此在高二那边水了个 $\texttt{rk4}$ (逃

 

模拟退火解旅行商问题

刚刚说模拟退火TSP随便跑...那么咱们就说说TSP怎么跑...

有人可能会问了这个东西怎么求下一个解?

其实仍是随机...

对于TSP, 咱们的一个方案其实就是一个遍历顺序(也就是一个排列)

这时咱们在生产新解的时候能够选择随机选取两个结点, 而后将它们在排列中的位置交换一下(好暴力啊)

然而事实证实效果拔群...

 

总结

模拟退火在OI中是一种在最优化问题中骗分的好方法

对于一些奇奇怪怪的多元函数也能够用这个方法来求解

其实在上面的例子中也能够体现出来, 这个算法的要点在于新解的选取以及参数的调整...

实际上利用退火过程的性质大胆随机再配合调参经验通常效果拔群OwO

可是做为一个随机化算法并不必定能找到最优解qwq...IOI赛制/ACM赛制的话可能骗分更容易些?(毕竟能够屡次提交233)