Loading

随机化算法入门——模拟退火

随机化算法入门——模拟退火

1 算法简介

模拟退火算法 (Simulated Annealing，SA) 最早的思想是由 N. Metropolis 等人于1953 年提出。1983， S. Kirkpatrick 等成功地将退火思想引入到组合优化领域。它是基于 Monte-Carlo 迭代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降，结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是一种通用的优化算法，理论上算法具有概率的全局优化性能,目前已在工程中得到了广泛应用，诸如 VLSI 、生产调度、控制工程、机器学习、神经网络、信号处理等领域。

模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。

就像这个图一样：

通过图我们不难发现，一开始温度升高这个指针跳的很厉害，但是随着温度降低，我们会逐渐逼近全局最优解。

也就是说，温度越高，对于一个不优的解，能被接受的概率也就越高。

2 算法实现

因为是模拟退火过程，所以有以下几个参数：

  • 初温 \(t_0\)
  • 温度变化量 \(\Delta\)
  • 末温度 \(eps\)

在打模拟退火的时候基本上是一个模板，下面我们以例题为例子，讲解模拟退火算法。

2.1 算法初实现

例题

  • 简述题意：给你 \(n\) 点，求出费马点，需要输出费马点到 \(n\) 个点的距离之和。

我们的模拟退火过程如下：

struct point{
    dd x,y;
    point(){}
    point(dd x,dd y) : x(x),y(y) {}
};
point p[N];

int T,n;
dd sumx,sumy,ans,ansx,ansy;

inline dd dis(point a,point b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

inline dd all_dis(dd x,dd y){
    dd res=0;
    for(int i=1;i<=n;i++) res+=dis(point(x,y),p[i]);
    return res;
}

inline void simulate(){
    dd x=ansx,y=ansy;
    dd t=t0;
    while(t>eps){
        dd dx=x+((rand()<<1)-RAND_MAX)*t;
        dd dy=y+((rand()<<1)-RAND_MAX)*t;
        dd now=all_dis(dx,dy);dd del=now-ans;
        if(del<0){
            x=dx;y=dy;ans=now;ansx=dx;ansy=dy;
        }
        else if(exp(-del/t)*RAND_MAX>rand()) x=dx,y=dy;
        t*=delta;
    }
}

inline void work(){
    ansx=sumx/n;ansy=sumy/n;
    for(int i=1;i<=5;i++) simulate();
}

稍后会展示全部程序。程序中 \(sumx,sumy\) 分别是所有点横坐标，纵坐标之和。

我们关注到 \(simulate\) 函数，也就是我们模拟退火算法的主函数，而 \(work\) 函数决定了我们要跑多少次模拟退火，\(ansx,ansy\) 存储的是我们目前最优解所对应的点。注意，及时模拟退火算法求得是一个值，我们一般还是要保存取到这个值的那个状态是多少，这个所谓“状态“可以是一个点，可以是一个序列，也可以是一个矩阵。

那么每一次我们用 \(x,y\) 存一下当前最优解，然后随机一个新的点，算一下这个点的值，如果比当前最优解优，我们就用这个解去更新最优解。——这是程序第 \(25\) 到 \(30\) 行所干的事情。

但是，如果这个解不优我们就不接受它了吗？显然不行。因为如果这样我们可能就会陷入局部最优解，所以对于一个不优的解我们还是需要一定的概率接受它，但是我们同时还需要保证随着温度的降低，我们接受的可能性就越小。

所以第 \(31\) 行就是我们看是否需要接受不优的解。我们下面解释一下为什么需要这么写。

注意这里的 \(del\) 是一个正数，然后 \(-del\) 是一个负数，当 \(t\) 越小时 \(-del/t\) 就越小，所以左边越小，大于右边的概率就越小，所以我们接受新解的概率也就越小。

现在挂上全部的程序：

#include<bits/stdc++.h>
#define dd double
#define ld long double
#define ll long long
#define uint unsigned int
#define ull unsigned long long
#define N 50010
#define M number
using namespace std;

const int INF=0x3f3f3f3f;
const dd delta=0.996;
const dd t0=3000;
const dd eps=1e-15;

template<typename T>  inline void read(T &x) {
    x=0; int f=1;
    char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c == '-') f=-f;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0';
    x*=f;
}

struct point{
    dd x,y;
    point(){}
    point(dd x,dd y) : x(x),y(y) {}
};
point p[N];

int T,n;
dd sumx,sumy,ans,ansx,ansy;

inline dd dis(point a,point b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

inline dd all_dis(dd x,dd y){
    dd res=0;
    for(int i=1;i<=n;i++) res+=dis(point(x,y),p[i]);
    return res;
}

inline void simulate(){
    dd x=ansx,y=ansy;
    dd t=t0;
    while(t>eps){
        dd dx=x+((rand()<<1)-RAND_MAX)*t;
        dd dy=y+((rand()<<1)-RAND_MAX)*t;
        dd now=all_dis(dx,dy);dd del=now-ans;
        if(del<0){
            x=dx;y=dy;ans=now;ansx=dx;ansy=dy;
        }
        else if(exp(-del/t)*RAND_MAX>rand()) x=dx,y=dy;
        t*=delta;
    }
}

inline void work(){
    ansx=sumx/n;ansy=sumy/n;
    for(int i=1;i<=5;i++) simulate();
}

int main(){
    srand(time(0));
    read(T);
    while(T--){
        ans=INF;sumx=0;sumy=0;
        read(n);
        for(int i=1;i<=n;i++){
            read(p[i].x);read(p[i].y);
            sumx+=p[i].x;sumy+=p[i].y;
        }
        work();
//        printf("%lf\n",ans);
        printf("%.0lf\n",round(ans));
        if(T) printf("\n");
    }
    return 0;
}

2.2 最大值

例题

\(1.1\) 中的例题是找最小值，那么我们怎么找最大值呢？

其实需要改变的只有我们接受不优解的概率，我们这样来找我们的最大值：

if(exp(-(dd)del/t)*RAND_MAX<rand())

你会发现这个时候 \(-del\) 是一个正数，那么 \(t\) 越小，\(-del/t\) 越大，左边就越大，接受不优解的概率就越低。

其余都很好写，我们放上总代码：

#include<bits/stdc++.h>
#define dd double
#define ld long double
#define ll long long
#define uint unsigned int
#define ull unsigned long long
#define N 15
#define M 3010
using namespace std;

const int INF=0x3f3f3f3f;
const dd eps=1e-10;
const dd t0=4000;
const dd delta=0.9972;

template<typename T>  inline void read(T &x) {
    x=0; int f=1;
    char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c == '-') f=-f;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0';
    x*=f;
}

int n,m,r;
dd ans;
dd ansx,ansy;

struct point{
    dd x,y;
    point() {}
    point(dd x,dd y) : x(x),y(y) {}
};
point b[M];
struct Round{
    dd x,y,r;
};
Round a[N];

inline dd dis(point a,point b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

inline dd compeat(dd x,dd y){
    dd res=0;
    dd min_dis=INF;
    for(int i=1;i<=n;i++) min_dis=min(min_dis,dis(point(x,y),point(a[i].x,a[i].y))-a[i].r);
    min_dis=min(min_dis,(dd)r);
    for(int i=1;i<=m;i++){
        dd d=dis(b[i],point(x,y));
        if(d<=min_dis) res++;
    }
    return res;
}

inline void simulate(){
    dd x=ansx,y=ansy;
    dd t=t0;
    while(t>eps){
        dd dx=x+((rand()<<1)-RAND_MAX)*t;
        dd dy=y+((rand()<<1)-RAND_MAX)*t;
        int now=compeat(dx,dy);int del=now-ans;
        if(del>0){
            ans=now;
            x=dx;y=dy;
            ansx=dx;ansy=dy;
        }
        else if(exp(-(dd)del/t)*RAND_MAX<rand()) x=dx,y=dy;
        t*=delta;
    }
}

inline void work(){
    ansx=0,ansy=0;
    for(int i=1;i<=6;i++) simulate();
}

int main(){
    srand(time(0));
    read(n);read(m);read(r);
    for(int i=1;i<=n;i++) scanf("%lf%lf%lf",&a[i].x,&a[i].y,&a[i].r);
    for(int i=1;i<=m;i++) scanf("%lf%lf",&b[i].x,&b[i].y);
    work();
    printf("%.0lf",ans);
    return 0;
}

但其实我们还可以这样写：

else if(exp((dd)del/t)*RAND_MAX>rand()) x=dx,y=dy;

2.3 如何调参

通常来说，我们初温定在 \(2500\) 左右，\(\Delta\) 约为 \(0.96\) 到 \(0.99\) ，然后末温约是 \(1e-10\) 到 \(1e-15\) 。下面我会不定期的更新一些其他的规律。

  • 如果需要一开始波动大一点，就把初温调小，\(\Delta\) 接近 \(1\)。这样写可以让初始状态不是很劣，从而逼近最优解。

2.4 技巧

  • 卡时：while((dd)clock()/CLOCKS_PER_SEC<=0.6) simulate();

3 习题

#include<bits/stdc++.h>
#define dd double
#define ld long double
#define ll long long
#define uint unsigned int
#define ull unsigned long long
#define N 51
#define M number
using namespace std;

const int INF=0x3f3f3f3f;
const dd eps=1e-10;
const dd delta=0.99;
const dd t0=4000;

inline void Swap(int &x,int &y){
    x^=y;y^=x;x^=y;
}

inline int Max(int a,int b){
    return a>b?a:b;
}

inline int Min(int a,int b){
    return a<b?a:b;
}

template<typename T>  inline void read(T &x) {
    x=0; int f=1;
    char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c == '-') f=-f;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0';
    x*=f;
}

int n,m,k,dis[N][N],r[N],d[N];
int ansxu[N],ans,tail,a[N],b[N];
bool vis[N];

inline int compeat(int *c){
    int res=-INF;
    for(int i=1;i<=k;i++) vis[c[i]]=1;
    for(int i=k+1;i<=tail;i++){
        int nowres=INF;
        for(int j=0;j<n;j++){
            if(!vis[j]||j==c[i]) continue;
            nowres=Min(nowres,dis[c[i]][j]);
        }
        res=Max(res,nowres);
    }
    for(int i=1;i<=k;i++) vis[c[i]]=0;
    return res;
}

inline void simulate(){
    for(int i=1;i<=tail;i++) a[i]=ansxu[i];
    dd t=t0;
    while(t>=eps){
        for(int i=1;i<=tail;i++) b[i]=a[i];
        int x=rand()%tail+1,y=rand()%tail+1;
        swap(b[x],b[y]);
        int now=compeat(b);int del=now-ans;
        if(del<0){
            ans=now;
            for(int i=1;i<=tail;i++) ansxu[i]=a[i]=b[i];
        }
        else if(exp(-(dd)del*t)*RAND_MAX>rand()) for(int i=1;i<=tail;i++) a[i]=b[i];
        t*=delta;
    } 
}

inline void work(){
    ans=INF;
    for(int i=1;i<=6;i++) simulate();
}

int main(){
    srand(time(0));
    read(n);read(m);read(k);
    for(int i=0;i<n;i++) read(r[i]);
    for(int i=0;i<n;i++) read(d[i]);
    for(int i=1;i<=m;i++){
        int x;read(x);
        vis[x]=1;
    }
    memset(dis,INF,sizeof(dis));
    for(int i=0;i<n;i++){
        if(!vis[i]) ansxu[++tail]=i;
        dis[i][r[i]]=dis[r[i]][i]=Min(dis[r[i]][i],d[i]);;
        dis[i][i]=0;
    }
//    printf("tail:%d\n",tail);
//    for(int i=1;i<=tail;i++) printf("%d ",ansxu[i]);
    for(int i=0;i<n;i++) for(int j=0;j<n;j++) for(int k=0;k<n;k++)
        dis[j][k]=Min(dis[j][k],dis[j][i]+dis[i][k]);
//    while(1){
//        int a,b;read(a);read(b);
//        printf("dis:(%d %d) is %d\n",a,b,dis[a][b]);
//    }
    if(m+k==n) return printf("0\n"),0;
    work();
    printf("%d\n",ans);
    return 0;
}

#include<bits/stdc++.h>
#define dd double
#define ld long double
#define ll long long
#define uint unsigned int
#define ull unsigned long long
#define N 21
#define M number
using namespace std;

const int INF=0x3f3f3f3f;
const dd eps=1e-15;
const dd delta=0.9999;
const dd t0=0.0005;

template<typename T>  inline void read(T &x) {
    x=0; int f=1;
    char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c == '-') f=-f;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0';
    x*=f;
}

int n,m,c,p[N],ans,ansxu[N][N],a[N][N],b[N][N];
int P[N],cnt[N];

#define For(i,n) for(int i=1;i<=n;i++)
inline int compeat(){
    int res=0;
    For(i,n) For(j,m){
        if(i!=n) if(b[i][j]!=b[i+1][j]) res++;
        if(j!=m) if(b[i][j]!=b[i][j+1]) res++;
    }   
    return res;
}

inline void simulate(){
    dd t=t0;
    For(i,n) For(j,m) a[i][j]=ansxu[i][j];
    while(t>=eps){
        For(i,n) For(j,m) b[i][j]=a[i][j];
//        if(!check()){
//            printf("here1\n");
//            while(1);
//        }
        int x1=rand()%n+1,x2=rand()%n+1;
        int y1=rand()%m+1,y2=rand()%m+1;
//        printf("%d %d\n",b[x1][y1],b[x2][y2]);
        swap(b[x1][y1],b[x2][y2]);
//        printf("%d %d\n",b[x1][y1],b[x2][y2]);
//        for(int i=1;i<=n;i++){
//            for(int j=1;j<=m;j++) printf("%d ",ansxu[i][j]);
//            putchar('\n');
//        }
//        putchar('\n');
//        if(!check()){
//            for(int i=1;i<=c;i++) printf("%d ",cnt[i]);
//            while(1); 
//        }
        int now=compeat();int del=now-ans;
        if(del<0){
            ans=now;
            For(i,n) For(j,m) ansxu[i][j]=a[i][j]=b[i][j];
        }
        else if(exp(-(dd)del/t)*(RAND_MAX)>rand()) For(i,n) For(j,m) a[i][j]=b[i][j];
        t*=delta;
    }
}

inline void work(){
    ans=INF;
    while((dd)clock()/CLOCKS_PER_SEC<=4.2)
//    for(int i=1;i<=5;i++)
        simulate();
}

int main(){
    srand(time(0));
    read(n);read(m);read(c);
    for(int i=1;i<=c;i++) read(p[i]),P[i]=p[i];
    int k=1;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=m;j++){
            while(!P[k]) k++;
            ansxu[i][j]=k;P[k]--;
        }
//    for(int i=1;i<=n;i++){
//        for(int j=1;j<=m;j++) printf("%d ",ansxu[i][j]);
//        putchar('\n');
//    }
//    putchar('\n');

    work();
//    if(!check()){
//        printf("here\n");
//        while(1);
//    }
//    printf("%d\n",ans);
    for(int i=1;i<=n;i++){
        for(int j=1;j<=m;j++) printf("%d ",ansxu[i][j]);
        putchar('\n');
    }
    return 0;
}
posted @ 2021-06-22 11:57  hyl天梦  阅读(438)  评论(0)  编辑  收藏  举报