模拟退火学习笔记

网上的绝大部分博客，包括洛谷题解原来都是错的写法……具体看转移部分QWQ

终于成为了少数派，这是好的

模拟退火的由来

模拟退火算法（Simulate Anneal，SA）是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。

模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的。V.Černý在1985年也独立发明此演算法。模拟退火算法是解决TSP问题的有效方法之一。

模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。

模拟退火算法是一种通用的优化算法，其物理退火过程由加温过程、等温过程、冷却过程这三部分组成。

——百度百科

简单说，模拟退火是一种随机化算法。当一个问题的方案数量极大（甚至无穷）而且不是单峰函数时，常使用模拟退火求解

我其实一直觉得模拟退火挺玄学的。（后来发现还真就靠RP吃饭

Metropolis接受原则

每次随机出一个新解，如果这个解更优，则接受它，否则以一个与温度和与当前最优解的差相关的概率接受它。

设这个新的解与最优解的差为 (Delta E)，温度为 (T) ， (k) 为一个随机数，那么这个概率为 (e^{frac{Delta E}{kT}})

具体来说：设 (P) 为接受当前这个解的概率，那么

若 (Delta E >0) ，那么 (P=1) ，也就是必定接受更优解；

若 (Delta E <0)，那么 $P=exp( -frac{Delta E}{kT} ) $，也就是以一个和差相关的概率接受这个不优的解。

对第二种情况的具体做法将在后面进行详细解释。

SA的简要流程

所需参数： 初始温度 (Tb)，降温的系数 (D)，最终温度 (Te)

（注：一般来讲 (Tb) 会比较大，(D) 是一个小于 (1) 且接近 (1) 的 (double) ，比如 (0.97) , (Te) 是一个接近 (0) 的正值，如 (1e-10) ，这些都是程序中自定的常数，但是在实际情况中这些参数需要进行不断调优，也就是所谓的调参。这或许也是模拟退火很难应用到 OI 赛制中的原因，毕竟不能实时查看结果的话就真的很靠 RP ，当然调参也是有迹可循的，比如 (TLE) 的时候可以考虑把初温调小，(WA) 的时候调大）

流程如下：

1. 令 (T=Tb)
2. 进行转移
3. (T*=D)
4. 重复步骤 (2、3) 直到 (T<Te)

转移（重点）

设估价函数为 (f(x)) （其实就是对你随机出来的方案计算答案的一个函数），(x)为原先得到的解，(x1=x+T0*rand()) 为新得到的随机出来的解（ (-1<=R<=1),且(R!=0) ），(Delta E =f(x1)- f(x))

虽然我觉得到这里还是一头雾水，但是看到代码应该就比较清楚

放个伪代码

if  ( DelE>0 )  更新解；
else if ( exp(DelE/sqrt(T))>(double)rand()/1215 ) 更新解；

然后发现一个奇怪的东西：if ( exp(DelE/sqrt(T))>(double)rand()/1215 ) 这是干嘛用的？

exp(DelE/sqrt(T)) ，由于前面的 if 语句，DelE/sqrt(T) 显然小于 (0) ，(exp) 之后也就是一个 (0sim 1) 的小数。下面来分析为什么这句话能达到 “以和温度、最优解相关的概率接受这个不优的解”。

根据简单的数学知识我们知道，DelE/sqrt(T) 越大，(exp) 之后的结果越大。而 DelE/sqrt(T) 的结果，当 DelE 较大或者 T 较大的时候都会变大，也就是说，当退火过程刚刚开始，或者当前结果和目前的最优解差距较小时，exp(DelE/sqrt(T)) 越大，也就更容易使 >rand()/1215 成立，这个解更容易被接受，反之亦然。

注意 ：if(exp(-del / ts) < (double)rand() / RAND_MAX) 这样的写法是错误的！ 显然 exp(-del/ts) 大于 (1) ，而 rand()/RAND_MAX 小于 (1) ，一定不会成立。但是不知道为什么就是能过去，反正是错的 感谢 @cjlworld 指出这个写法的问题。

一些细节

1. 随机的初始化真的很重要

2. 调参这个东西有时候很考经验

3. 其实模拟退火可以解决基本所有最优解问题 (所以快去，用SA全做一遍

4. 有时候为了求得最优解，需要多次SA，这里有个不TLE的好办法：

while ( (double)clock()/CLOCKS_PER_SEC<=0.8 ) SA();

(0.8) 根据时限而定，一般比时限小一点点

模板

luogu P5544 炸弹攻击1

//Author:RingweEH
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <ctime>
#define ll long long
#define db double
using namespace std;
int min( int a,int b ) { return (a<b) ? a : b; }
int max( int a,int b ) { return (a>b) ? a : b; }
void bmin( int &a,int b ) { a=(a<b) ? a : b; }
void bmax( int &a,int b ) { a=(a>b) ? a : b; }
int read()
{
    int x=0,w=1; char ch=getchar();
    while ( ch>'9' || ch<'0' ) { if ( ch=='-' ) w=-1; ch=getchar(); }
    while ( ch<='9' && ch>='0' ) { x=x*10+ch-'0'; ch=getchar(); }
    return x*w;
}

const int N=11,M=1010;
const double delta=0.996,Te=1e-10;          //D和最终温度
int n,m,x[N],y[N],r[N],p[M],q[M],R,ansout=0;
double ansx,ansy;
 
double dis( double ax,double ay,double bx,double by )		//求两点间距离
{
    double res=(bx-ax)*(bx-ax)+(by-ay)*(by-ay);
    res=sqrt(res);
    return res;
}
 
int calc( double xx,double yy )			//估价函数
{
    double nowr=R;
    for ( int i=1; i<=n; i++ )
        nowr=min( nowr,dis(xx,yy,x[i],y[i])-(double)r[i] );
    int cnt=0;
    for ( int i=1; i<=m; i++ )
        if ( dis(xx,yy,p[i],q[i])<=nowr ) cnt++;
    return cnt;
}
 
void SA()
{
    double T=6000;		//初始温度
    int ans=0;
    while ( T>Te )
    {
        double nowx=ansx+(rand()*2-RAND_MAX)*T;
        double nowy=ansy+(rand()*2-RAND_MAX)*T;
        int nowans=calc( nowx,nowy );
        int DelE=nowans-ans;
        if ( DelE>0 ) ansx=nowx,ansy=nowy,ans=nowans,ansout=max( ansout,ans );
        //更优的解，必然接受
        else if ( exp(DelE/sqrt(T))>(double)rand()/1215 ) ansx=nowx,ansy=nowy,ans=nowans;
        //按概率接受，保证了越到后期，和最优解的差距越大，越难被接受
        T*=delta;
    }
}
 
int main()
{
    srand(time(0)); ansx=ansy=0;
    n=read(); m=read(); R=read();
    for ( int i=1; i<=n; i++ )
        x[i]=read(),y[i]=read(),r[i]=read();
    for ( int i=1; i<=m; i++ )
        p[i]=read(),q[i]=read(),ansx+=x[i],ansy+=y[i];
    
    ansx/=m; ansy/=m; ansout=calc( ansx,ansy );
    while ( (double)clock()/CLOCKS_PER_SEC<=0.85 ) SA();
 
    printf( "%d",ansout );
}

习题

UVA10228-A Star not a Tree?

给定一个 (N) 边形所有顶点坐标 (x,y) ，求其费马点到所有顶点距离和。费马点是指到多边形所有顶点距离和最小的点

Solution

也是比较板子的题，但是注意多组数据，计时的 CLOCK 是总时间。

如果用 while ( (double)clock()/CLOCKS_PER_SEC<=0.8 ) ，这样的话就会出问题。所以要改成固定次数。

//Author:RingweEH
const int N=110;
const double delta=0.998,Te=1e-10;
int n,m,x[N],y[N];
double ansx,ansy;

double dis( double ax,double ay,double bx,double by )
{
    double res=(bx-ax)*(bx-ax)+(by-ay)*(by-ay);
    res=sqrt(res);
    return res;
}

double calc( double xx,double yy )
{
    double res=0;
    for ( int i=1; i<=n; i++ )
        res+=dis( xx,yy,x[i],y[i] );
    return res;
}

void SA()
{
    double T=1e7; 
    while ( T>Te )
    {
        double nowx=ansx+(rand()*2-RAND_MAX)*T;
        double nowy=ansy+(rand()*2-RAND_MAX)*T;
        double nowans=calc( nowx,nowy ),ans=calc(ansx,ansy);
        double DelE=nowans-ans;
        if ( DelE<0 ) ansx=nowx,ansy=nowy;
        else if ( exp(-DelE/T)>(double)rand()/RAND_MAX ) ansx=nowx,ansy=nowy;
        T*=delta;
    }
}

int main()
{
    srand(time(0));  int Case=read();
    for ( int cas=0; cas<Case; cas++ )
    {
        if ( cas ) printf( "
" );
        ansx=ansy=0; n=read(); 
        for ( int i=1; i<=n; i++ )
            x[i]=read(),y[i]=read(),ansx+=x[i],ansy+=y[i];
        
        ansx/=n; ansy/=n;
        for ( int i=1; i<=100; i++ ) SA();
        printf( "%.0lf
",calc( ansx,ansy ) );
    }

    return 0;
}

后记

蒟蒻初学SA，若有误希望指正QWQ