浅谈模拟退火

萌萌随机算法。


简介

模拟退火,是一种随机化算法,当一个问题的方案数量极大(甚至是无穷的),但是却都不是一个单峰函数时,我们常使用模拟退火求解。

什么是退火?

退火是一种金属热处理工艺,指的是将金属缓慢加热到一定温度,保持足够时间,然后以适宜速度冷却。目的是降低硬度,改善切削加工性;消除残余应力,稳定尺寸,减少变形与裂纹倾向;细化晶粒,调整组织,消除组织缺陷。准确的说,退火是一种对材料的热处理工艺,包括金属材料、非金属材料。而且新材料的退火目的也与传统金属退火存在异同。

——节选自百度百科

由于退火的规律引入了更多随机因素,得到最优解的概率会大大增加。于是我们可以去模拟这个过程,将目标函数作为能量函数,这也是模拟退火算法的由来。

过程

省流:如果随机出的新状态更优则修改答案,否则以一定概率接受新状态。

定义当前温度为 \(T\),当前状态为 \(S\) 且新状态为 \(S'\),并定义 \(\Delta E = |S-S'|\) 表示两个状态之差。那么接受新状态的概率即为:

\[P(\Delta E)= \begin{cases} 1 & S'\text{ is better than }S\\ \mathrm{e}^{\frac{-\Delta E}{T} } & \text{otherwise.}\\ \end{cases}\]

在某些时候,我们为了让最终得到的解更加靠谱,会在模拟退火结束后,以当前温度在得到的解附近多次随机状态,尝试得到更优的解。

退火方式

退火时我们有三个参数:初始温度 \(T_{\text{start}}\),降温系数 \(D\),以及终止温度 \(T_{\text{end}}\)。其中 \(T_{\text{start}}\) 是一个比较大的数,\(D\) 是一个接近 \(1\) 但小于 \(1\) 的实数(如 \(0.99\)),\(T_{\text{end}}\) 是一个接近 \(0\) 的正数(如 \(0.0001\))。

首先让当前温度 \(T = T_{\text{start}}\),然后按照【过程】部分所述的步骤尝试进行一次转移,接着让 \(T = D \times T\)。当 \(T < T_{\text{end}}\) 时停止退火,得到的答案就是最终的最优解。

注意为了使解更精确,我们通常不直接取当前解作为答案,而是在退火过程中维护遇到的所有解的最优值。

盗用借用 OI-wiki 上的一张动图来形象地展示一下退火的过程。

随着温度 \(\text{Temperature}\) 的降低,跳跃越来越不随机,最优解也越来越稳定。

例题讲解

P1337 平衡点 / 吊打XXX

怎么是物理题啊,欺负小学生吗。

在这道题中,一切自然变化进行的方向都是使能量降低,因为能量较低的状态比较稳定。

因为物体的重量是固定的,绳子越短,重物越低,势能越小,势能又与物重成正比,所以,只要使得 \(\sum_{i=1}^{n} dist_i \times weigh_i\) 也就是总的重力势能最小,就可以使整体平衡。

#include<bits/stdc++.h>
#define LL long long
#define UInt unsigned int
#define ULL unsigned long long
#define LD long double
#define pii pair<int,int>
#define pLL pair<LL,LL>
#define pDD pair<LD,LD>
#define fr first
#define se second
#define pb push_back
#define isr insert
using namespace std;
const LD eps = 1e-3;
const LD D = 0.998;
struct node{LD x,y,w;}a[10005];
LL n;LD Ans,ansx,ansy;
LL read(){
    LL su=0,pp=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')pp=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){su=su*10+ch-'0';ch=getchar();}
    return su*pp;
}
LD calc(LD X,LD Y){
    LD res=0;
    for(int i=1;i<=n;i++){
        LD dx=a[i].x-X,dy=a[i].y-Y;
        res+=sqrt(dx*dx+dy*dy)*a[i].w;
    }if(res<Ans)Ans=res,ansx=X,ansy=Y;
    return res;
}
void Solve(){
    LD T=100000.0,nowx=ansx,nowy=ansy;
    while(T>eps){
        LD mx=nowx+T*((1.0*rand()/RAND_MAX)*2-1);
        LD my=nowy+T*((1.0*rand()/RAND_MAX)*2-1);
        LD Delta=calc(mx,my)-calc(nowx,nowy);
        if(exp(-Delta/T)>(1.0*rand()/RAND_MAX))
            nowx=mx,nowy=my;T*=D;
    }for(int i=1;i<=1000;i++){
        LD mx=ansx+T*((1.0*rand()/RAND_MAX)*2-1);
        LD my=ansy+T*((1.0*rand()/RAND_MAX)*2-1);
        calc(mx,my);
    }return;
}
int main(){
    srand(time(0));
    n=read();
    for(int i=1;i<=n;i++)
        a[i].x=read(),a[i].y=read(),a[i].w=read(),
        ansx+=a[i].x,ansy+=a[i].y;
    ansx/=n,ansy/=n,Ans=calc(ansx,ansy);
    Solve();printf("%.3Lf %.3Lf\n",ansx,ansy);
    return 0;
}

P5544 炸弹攻击1

考虑对于每个点,计算得到距离该点最近的建筑边界,然后寻找这个范围内的敌人数量。基本就是纯板子了。

但是你这么写过不了 Hack,是 unAccepted 的 \(100\)pts,因为某个【】洛谷用户在捣乱,故意卡模拟退火。但是你只要多乱搞一下就可以过了。

由于在考场上没人会想着用模拟退火冲正解,它本来就是用来捞分的,所以我就放过不了 Hack 的普通模拟退火代码了。

#include<bits/stdc++.h>
#define LL long long
#define UInt unsigned int
#define ULL unsigned long long
#define LD long double
#define pii pair<int,int>
#define pLL pair<LL,LL>
#define pDD pair<LD,LD>
#define fr first
#define se second
#define pb push_back
#define isr insert
using namespace std;
const LD eps = 1e-3;
const LD D = 0.98971;
struct buildings{LD x,y,r;}a[105];
struct boss_{LD x,y;}b[1005];
LL n,m,Ans;LD R,sx,sy;
LL read(){
    LL su=0,pp=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')pp=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){su=su*10+ch-'0';ch=getchar();}
    return su*pp;
}
LD Dist(LD x_1,LD y_1,LD x_2,LD y_2){
    return sqrt((x_1-x_2)*(x_1-x_2)+(y_1-y_2)*(y_1-y_2));
}
LL calc(LD X,LD Y){
    LD res=R;LL sum=0;
    for(int i=1;i<=n;i++)
        res=min(res,Dist(X,Y,a[i].x,a[i].y)-a[i].r);
    if(res<0)return 0;
    for(int i=1;i<=m;i++)
        if(Dist(X,Y,b[i].x,b[i].y)<=res)sum++;
    return sum;
}
LD Rand01(){
    return (1.0*rand()/RAND_MAX);
}
void Rchange(LD &X,LD &Y){
    X=2*R*Rand01()-R,Y=2*R*Rand01()-R;return;
}
void Solve(){
    LD T=R;LL now=calc(sx,sy);
    while(T>eps){
        LD tx=sx+2*T*Rand01()-T;
        LD ty=sy+2*T*Rand01()-T;
        LL tmp=calc(tx,ty);
        if(tmp>now||exp((LD)10000*(tmp-now)/T)>Rand01())
            now=tmp,sx=tx,sy=ty;
        Ans=max(Ans,now);T*=D;
    }return;
}
int main(){
    srand(time(0));
    n=read(),m=read(),R=read();
    for(int i=1;i<=n;i++)
        a[i].x=read(),a[i].y=read(),a[i].r=read();
    for(int i=1;i<=m;i++)
        b[i].x=read(),b[i].y=read();
    for(int opt=1;opt<=25;opt++)
        Rchange(sx,sy),Solve();
    cout<<Ans<<"\n";
    return 0;
}

P4035 球形空间产生器

这哪里叫模拟退火——这分明是爬山算法好吗!

随机一个点当圆心,然后和每个点比较,求出平均距离 \(r\) ,如果到这个点的距离大于 \(r\),说明离这个点远了,就给圆心施加一个向这个点的力;若小于 \(r\),说明近了,就施加一个远离这个点的力。所有点比较完后,把假设的圆心按合力方向移动一个距离,距离和当前温度有关。最后找到最优的圆心即可。

#include<bits/stdc++.h>
#define LL long long
#define UInt unsigned int
#define ULL unsigned long long
#define LD long double
#define pii pair<int,int>
#define pLL pair<LL,LL>
#define pDD pair<LD,LD>
#define fr first
#define se second
#define pb push_back
#define isr insert
using namespace std;
const int N = 20;
const LD eps = 1e-4;
const LD D = 0.99995;
int n;
LD a[N][N],dis[N],c[N],ans[N];
int read(){
    int su=0,pp=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')pp=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){su=su*10+ch-'0';ch=getchar();}
    return su*pp;
}
void SA(){
    LD T=10000;
    while(T>eps){
        LD res=0;
        for(int i=1;i<=n+1;i++){
            dis[i]=0,c[i]=0;
            for(int j=1;j<=n;j++)
                dis[i]+=(a[i][j]-ans[j])*(a[i][j]-ans[j]);
            dis[i]=sqrt(dis[i]),res+=dis[i];
        }res/=n+1;
        for(int i=1;i<=n+1;i++)for(int j=1;j<=n;j++)
            c[j]+=(dis[i]-res)*(a[i][j]-ans[j])/res;
        for(int i=1;i<=n;i++)ans[i]+=c[i]*T;
        T*=D;
    }return;
}
int main(){
    n=read();
    for(int i=1;i<=n+1;i++)
        for(int j=1;j<=n;j++)
            {cin>>a[i][j];ans[j]+=a[i][j];}
    for(int i=1;i<=n;i++)ans[i]/=n+1;
    SA();for(int i=1;i<=n;i++)printf("%.3Lf ",ans[i]);
    return 0;
}

P3878 分金币

非常板子的板子模拟退火,打乱数组顺序并从中间砍半分成两份,然后每次随机两个不同位置交换尝试转移。为了让答案更精确,可以多做几次模拟退火保证答案准确性。

#include<bits/stdc++.h>
#define LL long long
#define UInt unsigned int
#define ULL unsigned long long
#define LD long double
#define pii pair<int,int>
#define pLL pair<LL,LL>
#define pDD pair<LD,LD>
#define fr first
#define se second
#define pb push_back
#define isr insert
using namespace std;
const LD eps = 1e-8;
const LD D = 0.98971;
LL T,n,a[50],Ans;
LL read(){
    LL su=0,pp=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')pp=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){su=su*10+ch-'0';ch=getchar();}
    return su*pp;
}
LL calc(){
    LL x=0,y=0;
    for(int i=1;i<=n/2;i++)x+=a[i];
    for(int i=n;i>n/2;i--)y+=a[i];
    return abs(x-y);
}
void Solve(){
    LD T=100.0;LL now=calc();Ans=min(Ans,now);
    while(T>eps){
        int p1=rand()%n+1,p2=rand()%n+1;
        if(p1==p2)continue;
        swap(a[p1],a[p2]);
        LL tmp=calc();Ans=min(Ans,tmp);
        if(exp((LD)(now-tmp)/T)<(1.0*rand()/RAND_MAX))swap(a[p1],a[p2]);
        else now=tmp;T*=D;
    }return;
}
int main(){
    srand(time(0));
    T=read();
    while(T--){
        n=read(),Ans=0x3f3f3f3f3f3f3f3f;
        for(int i=1;i<=n;i++)a[i]=read();
        if(n==1){cout<<a[1]<<"\n";continue;}
        for(int opt=1;opt<=50;opt++)
            random_shuffle(a+1,a+n+1),Solve();
        cout<<Ans<<"\n";
    }
    return 0;
}

P3936 Coloring

说是有点卡人的题,参数需要设得比较极端,虽然我没被卡过来着

首先我们可以发现,给定一个矩阵要求计算这个所谓 \(q\) 值其实只需要遍历矩阵一遍,记录每个方格颜色不同的临边,最后 \(\div 2\) 即可。

考虑每次交换两个颜色不一样的格子来尝试转移,注意更新 \(q\) 值的时候没必要每次都去遍历,只需要计算交换这两个格子带来的贡献就可以了。这样能快速提高效率,也能让我们退火的次数更多,进而使答案更加精确。

注意这个题目是需要存储矩阵具体情况的喔!

#include<bits/stdc++.h>
#define LL long long
#define UInt unsigned int
#define ULL unsigned long long
#define LD long double
#define pii pair<int,int>
#define pLL pair<LL,LL>
#define pDD pair<LD,LD>
#define fr first
#define se second
#define pb push_back
#define isr insert
using namespace std;
const LD eps = 1e-15;
const LD D = 0.99999;
struct Node{int nq,h[25][25];}Ans;
int n,m,col,p[55],a[25][25];
int dx[]={0,1,0,-1},dy[]={1,0,-1,0};
int read(){
    int su=0,pp=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')pp=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){su=su*10+ch-'0';ch=getchar();}
    return su*pp;
}
void start_init(){
    int now=1,cnt=0;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=m;j++)
            if(cnt<p[now])a[i][j]=now,cnt++;
            else a[i][j]=(++now),cnt=1;
    return;
}
int calc_answer(){
    int cnt=0;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=m;j++)
            for(int k=0;k<4;k++){
                int x=i+dx[k],y=j+dy[k];
                if(x<1||x>n||y<1||y>m)continue;
                if(a[i][j]!=a[x][y])cnt++;
            }
    return cnt/2;
}
LD Probly(){return (1.0*rand()/RAND_MAX);}
void updAns(int x){
    Ans.nq=x;
    for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)
        Ans.h[i][j]=a[i][j];return;
}
void solSA(){
    LD T=1.0;int nowq=Ans.nq;
    while(T>eps){
        int x=rand()%n+1,y=rand()%m+1;
        int x2=rand()%n+1,y2=rand()%m+1;
        if(a[x][y]==a[x2][y2])continue;
        int ord=0,nwl=0;
        for(int k=0;k<4;k++)
            ord+=(a[x+dx[k]][y+dy[k]]!=a[x][y]),
            ord+=(a[x2+dx[k]][y2+dy[k]]!=a[x2][y2]);
        swap(a[x][y],a[x2][y2]);
        for(int k=0;k<4;k++)
            nwl+=(a[x+dx[k]][y+dy[k]]!=a[x][y]),
            nwl+=(a[x2+dx[k]][y2+dy[k]]!=a[x2][y2]);
        int newq=nowq-ord+nwl;
        if(newq<=nowq||exp((LD)(nowq-newq)/T)>Probly())nowq=newq;
        else swap(a[x][y],a[x2][y2]);
        if(nowq<Ans.nq)updAns(nowq);T*=D;
    }return;
}
int main(){
    srand(time(0));
    n=read(),m=read(),col=read();
    for(int i=1;i<=col;i++)p[i]=read();
    start_init();updAns(calc_answer());
    for(int i=1;i<=3;i++)solSA();
    for(int i=1;i<=n;cout<<"\n"&&i++)
        for(int j=1;j<=m;j++)cout<<Ans.h[i][j]<<" ";
    return 0;
}

总结

模拟退火,是一种特别好用的随机算法,它模拟了金属退火的过程,通过慢慢降温和随机参数的调整,可以得到一个非常近似最优解的答案,在某些题目上可以捞到不少分,也是比赛时的赚分工具。对了,它还有一个卡时的方法,就是不断判断是否满足 ((double)clock()/CLOCKS_PER_SEC < MAX_TIME,如果满足则继续跑模拟退火,否则终止程序直接输出答案,这样能保证模拟退火算法在不超时的情况下得到尽可能优的解。其中的 clock() 函数是返回程序运行时间的函数,而 MAX_TIME 则是一个略小于题目限制的时间,如限制 \(1\)s 的题目可以设为 \(0.9\)。总之,模拟退火是一种非常棒的随机算法哦!

posted @ 2026-01-02 14:49  嘎嘎喵  阅读(7)  评论(0)    收藏  举报