zoukankan      html  css  js  c++  java
  • 「笔记」模拟退火

    写在前面

    感谢 caq 的倾情讲解 qq_emoji: bx

    模拟退火是个随机化算法,正确性有一定保证,但如果你想我一样脸黑的话......

    实测模拟退火做多了 rp 会掉

    正文

    简介

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

    什么是退火?

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

    qq_emoji: fad

    扯远了。

    这个算法就是在温度不断降低的过程中,不断地从当前位置寻找别的位置进行计算,温度越低,也就是它的动能越小时,位置就会变化的越小,最后逐渐停留在最优解(或者附近)

    算法流程

    每次随机一个新的状态,如果状态更优就更新答案,否则以一定概率接受这个状态。

    Metropolis准则

    以求最小值为例。

    • 如果 (Delta E < 0),说明当前解更优,直接更新即可

    • 否则,如果

    [e^{frac{-Delta E}{T}} > frac{ ext{rand()}}{ ext{RAND_MAX}} ]

    就接受这个状态。

    • 否则 直接跳过。

    为什么?
    第一步因为是最优解所以一定选择更新答案
    第二步后边的是一个随机值我们暂且不论。
    考虑整个退火过程,
    假设温度 (T) 不变,新的解越劣,(Delta E) 越大,左项的值越小,接受的概率也越小。
    假设 (Delta E) 不变,随着温度的下降,求解的范围也趋于稳定,(T) 越小,左项得值也越小,接受的概率也越小

    扔一张图可能会更好理解:

    听上去很扯 ,但它还是有一定的正确性的。

    SA 函数

    通常降温系数 (d) 是一个很接近 (1) 的数,终止温度 (T_0) 是一个很接近 (0) 的数

    这里给一个伪代码:

    const double lim = ... // 温度最小值,通常为 1e-10 左右
    const double d = ... // 变化系数,通常为 0.996 左右
    void SA() {
        double T = ... // 初始温度,通常为 2021 左右
        while(T > lim) {
            ... // 获取一个随机的位置
            now = calc(); // 计算当前位置的答案 
            del = now - ans; // 计算 变化量
            if(del < 0) { // 以最小值为例
                ans = now; // 更新答案
                ...  // 更新答案和中间量的状态
            } else if(exp(-del/T) > (double)rand()/RAND_MAX) {
                ...  // 一定概率选择当前当前状态
            } 
            T *= d; // 降温
        }
    }
    

    计算函数 calc

    依据题目而定,这里不给出

    一些技巧/思想

    如果想要随机一个无限大平面内的一个点,可以这样:

    double nowx = limx + ((rand() << 1) - RAND_MAX) * T;
    double nowy = limy + ((rand() << 1) - RAND_MAX) * T;
    

    其中 nowx,nowy 是我们随机的位置, limx, limy 是我们一个中间状态的位置(注意不是答案的位置),
    后面的那一坨刚好对应着温度越小变化越小的实际情况。


    我们有时为了使得到的解更有质量,会在模拟退火结束后,以当前温度在得到的解附近多次随机状态,尝试得到更优的解(其过程与模拟退火相似)。


    模拟退火是个随机的算法,执行次数越多获得的解越有可能更优,所以我们可以执行多遍 SA 函数。至于如何控制时间?

    while((double)clock()/CLOCKS_PER_SEC < 0.90) SA();
    

    上面这个代码控制时间在 (0.90s) 左右,如果时间限制为 (1s),而每次 SA 函数执行时间略长时,就要小心可能会 ( ext{T}) 掉了。


    如果一个代码不行,就考虑换个种子吧。

    srand(...);
    

    为了获得更精确的解,也可以把 (d)(T_0) 调的更精准一点

    const double d = 0.996 -> 0.99996;
    const double lim = 1e-10 -> 1e-15;
    

    还有,随机乱搞一些 初温,终温,降温系数 也是可以的。


    随机微调

    对于序列分组这类问题,每次可以随机两个位置交换计算新解的贡献。
    但要注意如果没有接受这个状态时要把两个位置在交换回去。


    打乱排序

    对于一些统计贡献与序列顺序无关的题可以考虑。

    函数 random_shuffle() 可以在 (O(n)) 的时间内随机打乱一个序列
    包含在 ctime 中,使用方法和 sort 类似。
    可以用下面那道 P2503 来练练手。

    Tips

    模拟退火通常用来解决一些数据范围比较小,但对正解的计算不好算的问题。
    通过随机一些状态来不断逼近正解。
    其计算过程通常是 简单粗暴 的。
    其正确性可能是 玄学 的。

    这些例题中都有体现,注意体会。

    例题

    UVA10228 A Star not a Tree?

    题目传送

    初始点的位置考虑选择所有点的平均值(人类本质?)

    按照上面说的每次随机一个点即可。

    如何计算一个点的贡献?(n) 的数据范围只有 (100),暴力计算即可。

    看一下代码可能会更好理解。

    /*
    Work by: Suzt_ilymics
    Problem: 不知名屑题
    Knowledge: 垃圾算法
    Time: O(能过)
    CAQ yyds !!!!
    */
    #include<iostream>
    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<cmath>
    #include<queue>
    #include<ctime>
    #define LL long long
    #define orz cout<<"lkp AK IOI!"<<endl
    
    using namespace std;
    const int MAXN = 1e5+5;
    const int INF = 1e9+7;
    const int mod = 1e9+7;
    const double d = 0.996;
    const double lim = 1e-10;
    
    int n;
    double ax[MAXN], ay[MAXN];
    double limx, limy; // 上一次选择保留的位置 
    double ansx, ansy; // 存答案的位置; 
    double ans;
    
    int read(){
        int s = 0, f = 0;
        char ch = getchar();
        while(!isdigit(ch))  f |= (ch == '-'), ch = getchar();
        while(isdigit(ch)) s = (s << 1) + (s << 3) + ch - '0' , ch = getchar();
        return f ? -s : s;
    }
    
    double dis(double sx, double sy, double ex, double ey) {
        return sqrt((sx - ex) * (sx - ex) + (sy - ey) * (sy - ey));
    }
    
    double calc(double sx, double sy) {
        double sum = 0.0;
        for(int i = 1; i <= n; ++i) {
            sum += dis(sx, sy, ax[i], ay[i]);
        }
        return sum;
    }
    
    void SA() {
        double T = 2021; // 设置初始温度 
        while(T > lim) { // 当没降完温继续降温 
            double nowx = limx + ((rand() << 1) - RAND_MAX) * T; // 独特的随机一个点的公式
            double nowy = limy + ((rand() << 1) - RAND_MAX) * T;
            double now = calc(nowx, nowy);//计算贡献
            double del = now - ans; // 计算 改变量
            if(del < 0) {
                ansx = nowx, ansy = nowy;
                limx = nowx, limy = nowy;
                ans = now;
            } else if(exp(-del / T) > (double)rand() / RAND_MAX) { // 概率的选择接受这个解 
                limx = nowx, limy = nowy;
            }
            T *= d; // 降温 
        }
    }
    
    void work() {
        limx /= (1.0 * n); // 取一个平均值
        limy /= (1.0 * n);
        ansx = limx, ansy = limy;
        ans = calc(limx, limy);
        for(int i = 1; i <= 10; ++i) SA(); // 多执行几遍 SA 函数
    }
    
    int main()
    {
        int t = read();
        for(int i = 1; i < t; ++i) { // UVA 独特的数据测试方式
            limx = limy = ansx = ansy = 0.0;
            ans = 0;
            n = read();
            for(int j = 1; j <= n; ++j) {
                scanf("%lf%lf", &ax[j], &ay[j]);
                limx += ax[j], limy += ay[j];
            }
            work();
            printf("%.0lf
    
    ", ans);
            
        }
        n = read();
        for(int j = 1; j <= n; ++j) {
            scanf("%lf%lf", &ax[j], &ay[j]);
            limx += ax[j], limy += ay[j];
        }
        work();
        printf("%.0lf
    ", ans);
        return 0;
    }
    

    P5544 [JSOI2016]炸弹攻击1

    题目传送

    肯定要先随机一个点,考虑如何确定 (R)

    随机吗?

    注意随机一个量的时候只能有一个参照量,也就是说这里不能同时随机点的位置和爆炸半径。
    因为我们无法简单的确定一个决策依据。

    人类本质?

    能炸多少炸多少啊!

    所以暴力找到最大爆炸范围统计贡献就行啊。

    /*
    Work by: Suzt_ilymics
    Problem: 不知名屑题
    Knowledge: 垃圾算法
    Time: O(能过)
    */
    #include<iostream>
    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<cmath>
    #include<queue>
    #include<cstdlib>
    #include<ctime>
    #define LL long long
    #define orz cout<<"lkp AK IOI!"<<endl
    
    using namespace std;
    const int MAXN = 1e5+5;
    const int INF = 1e9+7;
    const int mod = 1e9+7;
    const double d = 0.996;
    const double lim = 1e-10;
    
    struct node {
        int x, y, r;
    }a[MAXN];
    
    struct node2 {
        int x, y;
    }b[MAXN];
    
    int n, m, R;
    double limx, limy, ansx, ansy;
    int ans;
    
    int read(){
        int s = 0, f = 0;
        char ch = getchar();
        while(!isdigit(ch))  f |= (ch == '-'), ch = getchar();
        while(isdigit(ch)) s = (s << 1) + (s << 3) + ch - '0' , ch = getchar();
        return f ? -s : s;
    }
    
    double GetR(double sx, double sy, int x) {
        return sqrt((sx - a[x].x) * (sx - a[x].x) + (sy - a[x].y) * (sy - a[x].y)) - a[x].r;
    }
    
    double Dis(double sx, double sy, int x) {
        return sqrt((sx - b[x].x) * (sx - b[x].x) + (sy - b[x].y) * (sy - b[x].y));
    }
    
    int calc(double sx, double sy) {
        double maxR = 1.0 * R;
        for(int i = 1; i <= n; ++i) {
            maxR = min(maxR, GetR(sx, sy, i));
        }
        int cnt = 0;
        for(int i = 1; i <= m; ++i) {
            if(Dis(sx, sy, i) <= maxR) cnt++;
        }
        return cnt;
    }
    
    void SA() {
        double T = 2021;
        while(T > lim) {
            double nowx = limx + ((rand() << 1) - RAND_MAX) * T;
            double nowy = limy + ((rand() << 1) - RAND_MAX) * T;
            int now = calc(nowx, nowy);
            int del = now - ans;
            if(del > 0) {
                ansx = nowx, ansy = nowy;
                limx = nowx, limy = nowy;
                ans = now;
            } else if(exp(-del/T) < (double) rand() / RAND_MAX) {
                limx = nowx, limy = nowy;
            }
            T *= d;
        }
    }
    
    void work() {
        ansx = ansy = limx = limy = 0.0;
        ans = calc(limx, limy);
        while((double)clock() / CLOCKS_PER_SEC < 0.70) SA(); // CLOCKS_PER_SEC
    }
    
    int main()
    {
        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();
        work();
        printf("%d", ans);
        return 0;
    }
    

    P2503 [HAOI2006]均分数据

    题目传送

    原题面中的式子可能有点错误,应该把 (m) 组当成 (m) 个元素来计算方差

    介绍一个函数 random_shuffle(),可以在 (O(n)) 的时间内,将一段区间随机打乱。

    考虑最简单的想法,

    每次随机打乱,贪心的分组。

    每次加进一个元素,将它放进当前 (m) 组中最小的那一组中。

    然后计算贡献。

    多随机几遍,总有一遍我们会碰到最优解的。

    /*
    Work by: Suzt_ilymics
    Problem: 不知名屑题
    Knowledge: 垃圾算法
    Time: O(能过)
    */
    #include<iostream>
    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<cmath>
    #include<queue>
    #include<ctime>
    #define LL long long
    #define orz cout<<"lkp AK IOI!"<<endl
    
    using namespace std;
    const int MAXN = 1e5+5;
    const int INF = 1e9+7;
    const int mod = 1e9+7;
    const double d = 0.996;
    const double lim = 1e-10;
    
    int n, m;
    int a[22];
    int b[7];
    int sum;
    double ave, ans = 10000000.0;
    
    int read(){
        int s = 0, f = 0;
        char ch = getchar();
        while(!isdigit(ch))  f |= (ch == '-'), ch = getchar();
        while(isdigit(ch)) s = (s << 1) + (s << 3) + ch - '0' , ch = getchar();
        return f ? -s : s;
    }
    
    double calc() {
        memset(b, false, sizeof b);
        for(int i = 1; i <= n; ++i) {
            int wz = 1;
            for(int j = 1; j <= m; ++j) if(b[wz] > b[j]) wz = j;
            b[wz] += a[i];
        }
        double Sum = 0;
        for(int i = 1; i <= m; ++i) {
            Sum += (ave - 1.0 * b[i]) * (ave - 1.0 * b[i]);
        }
        Sum /= (1.0 * m);
        return sqrt(Sum);
    }
    
    void SA() {
        double T = 2021;
        while(T > lim) {
            random_shuffle(a + 1, a + n + 1);
            double now = calc();
            ans = min(ans, now);
    //        double del = now - ans;
    ////        cout<<now<<endl;
    //        if(del < 0) {
    //            ans = now;
    //        } 
            T *= d;
        }
    }
    
    void work() {
        ave = 1.0 * sum / m;
        for(int i = 1; i <= 20; ++i) SA();
    }
    
    int main()
    {
        n = read(), m = read();
        for(int i = 1; i <= n; ++i) a[i] = read(), sum += a[i];
        work();
        printf("%.2lf", ans);
        return 0;
    }
    

    P3878 [TJOI2010]分金币

    题目传送

    前一半分一组,后一半分一组。

    然后每次随机两个位置交换一下计算贡献。

    nowx = rand() % n + 1; // 整数位置的随机直接 % n + 1 即可
    

    这里与板子略有不同

    如果更优就更新,
    否则考虑概率接受,
    否则你要把元素再换回去!

    计算函数应该很好想,依旧是暴力,如果还不会就看代码吧。

    /*
    Work by: Suzt_ilymics
    Problem: 不知名屑题
    Knowledge: 垃圾算法
    Time: O(能过)
    */
    #include<iostream>
    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<cmath>
    #include<queue>
    #include<ctime>
    #define int long long
    #define LL long long
    #define orz cout<<"lkp AK IOI!"<<endl
    
    using namespace std;
    const int MAXN = 1e5+5;
    const int INF = 1e12+7;
    const int mod = 1e9+7;
    const double d = 0.996;
    const double lim = 1e-10;
    
    int t, n, ans = INF;
    int a[40];
    
    int read(){
        int s = 0, f = 0;
        char ch = getchar();
        while(!isdigit(ch))  f |= (ch == '-'), ch = getchar();
        while(isdigit(ch)) s = (s << 1) + (s << 3) + ch - '0' , ch = getchar();
        return f ? -s : s;
    }
    
    int calc() {
        int sum1 = 0, sum2 = 0;
        for(int i = 1; i <= n / 2; ++i) sum1 += a[i];
        for(int i = n / 2 + 1; i <= n; ++i) sum2 += a[i];
        return abs(sum1 - sum2);
    }
    
    void SA() {
        double T = 2021;
        while(T > lim) {
            int nowx = rand() % n + 1;
            int nowy = rand() % n + 1;
            swap(a[nowx], a[nowy]);
            int now = calc();
            int del = now - ans;
            if(del < 0) {
                ans = now;
            } else if(exp(-del/T) > (double)rand() / RAND_MAX) {
            } else swap(a[nowx], a[nowy]);
            T *= d;
        }
    }
    
    void work() {
        for(int i = 1; i <= 100; ++i) SA();
    }
    
    signed main()
    {
        t = read();
        while(t--) {
            n = read();
            ans = INF;
            for(int i = 1; i <= n; ++i) a[i] = read();
            random_shuffle(a + 1, a + n + 1);
            work();
            printf("%lld
    ", ans);
        }
        return 0;
    }
    

    P3936 Coloring

    题目传送

    不要害怕,加油,奥利给!

    选择按顺序染色作为初始状态。
    然后开始微调。

    如何更快的计算贡献?
    在原来的基础上减去两个位置原来位置的贡献,然后加上交换后的贡献。

    推荐参数:

    d = 0.9999lim = 1e-15T = 500srand(20041011)

    这样的话分数就可以在 (97pts) 以上了

    有个疑问?
    我计算 (del) 的时候是当前值与交换前的值作差,而不是和最优值作差。可前面几道题都是和最优值作差,可这道题这样做却死活过不去。有没有大佬可以解释一下啊,也可以到我这个帖子看一看。

    /*
    Work by: Suzt_ilymics
    Problem: 不知名屑题
    Knowledge: 垃圾算法
    Time: O(能过)
    */
    #include<iostream>
    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<cmath>
    #include<queue>
    #include<ctime>
    #define LL long long
    #define orz cout<<"lkp AK IOI!"<<endl
    
    using namespace std;
    const int MAXN = 1e5+5;
    const int INF = 1e9+7;
    const int mod = 1e9+7;
    const double d = 0.99998;
    const double lim = 1e-15;
    int dx[] = {1, -1, 0, 0};
    int dy[] = {0, 0, 1, -1};
    
    int n, m, c, ans, now_;
    int p[55], cnt[55], top = 1;
    int col[22][22], acol[22][22];
    
    int read(){
        int s = 0, f = 0;
        char ch = getchar();
        while(!isdigit(ch))  f |= (ch == '-'), ch = getchar();
        while(isdigit(ch)) s = (s << 1) + (s << 3) + ch - '0' , ch = getchar();
        return f ? -s : s;
    }
    
    bool check(int x, int y) { return x < 1 || y < 1 || x > n || y > m; }
    
    int Dev(int sx, int sy) {
        int Cnt = 0;
        for(int i = 0; i < 4; ++i) {
            int x = sx + dx[i], y = sy + dy[i]; // 找相邻的四个格子计算贡献 
            if(check(x, y)) continue;
            if(col[sx][sy] != col[x][y]) Cnt++; 
        }
        return Cnt;
    }
    
    //int Dev(int x, int y) { return (col[x][y] != col[x-1][y] && x!=1) + (col[x][y] != col[x][y-1] && y!=1) + (col[x][y] != col[x+1][y] && x != n) + (col[x][y] != col[x][y+1] && y != m); }
    
    int calc(int sx, int sy, int ex, int ey) {
        int res = 0;
    //    orz;
        res = res - Dev(sx, sy) - Dev(ex, ey); // 减去原来的贡献 
        swap(col[sx][sy], col[ex][ey]);
        res = res + Dev(sx, sy) + Dev(ex, ey); // 加上现在的贡献 
        return res;
    }
    
    void SA() {
    //    orz;
    //    for(int i = 1; i <= n; ++i) for(int j = 1; j <= m; ++j) col[i][j] = acol[i][j]; 
    //    now_ = ans;
        double T = 500;
        while(T > lim) {
            int nowsx = rand() % n + 1;
            int nowsy = rand() % m + 1;
            int nowex = rand() % n + 1;
            int nowey = rand() % m + 1;
    //        int pre = Calc();
    //        swap(col[nowsx][nowsy], col[nowex][nowey]); // 先交换 
    //        cout<<calc(nowsx, nowsy, nowex, nowey)<<"
    ";
            int now = now_ + calc(nowsx, nowsy, nowex, nowey); // 再计算 
    //        cout<<ans<<" "<<now_<<" "<<now<<"
    ";
            int del = now - now_; // now_ 存的是当前 col 矩阵的答案 
    //        int del = now - ans; // ans 存的是最优解 
            if(del < 0) { // 更新这个解
                ans = now; 
                now_ = now;
                for(int i = 1; i <= n; ++i) for(int j = 1; j <= m; ++j) acol[i][j] = col[i][j];  
            } else if(exp(-del/T) > (double) rand() / RAND_MAX) {
                now_ = now; // 一定概率接受当前这个状态 
            } else {
                swap(col[nowsx][nowsy], col[nowex][nowey]); // 交换回去 
            }
            T *= d;
        }
    }
    
    void work() {
    //    for(int i = 1; i <= n; ++i) {
    //        for(int j = 1; j <= m; ++j) {
    //            ans += Dev(i, j);
    //        }
    //    }
    //    ans /= 2;
        for(int i = 2; i <= n; ++i) 
            for(int j = 1;j <= m; ++j) 
                if(col[i][j] != col[i - 1][j]) ans++;
        for(int i = 1; i <= n; ++i) 
            for(int j = 2;j <= m; ++j) 
                if(col[i][j] != col[i][j - 1]) ans++;
        now_ = ans;
    //    cout<<ans<<" "<<now_<<endl;
    //    cout<<ans<<endl;
    //    for(int i = 1; i <= 30; ++i) SA();
        while((double)clock() / CLOCKS_PER_SEC < 3.90) SA();//, cout<<ans<<" "<<now_<<endl;
    }
    
    int main()
    {
        srand(20041011);
        n = read(), m = read(), c = read();
        for(int i = 1; i <= c; ++i) p[i] = read();
        for(int i = 1; i <= n; ++i) {
            for(int j = 1; j <= m; ++j) {
                if(cnt[top] == p[top]) ++top;
                acol[i][j] = col[i][j] = top;
                ++cnt[top];
            }
        }
        work();
    //    printf("%d
    ", ans);
        for(int i = 1; i <= n; ++i) {
            for(int j = 1; j <= m; ++j) {
                printf("%d ", acol[i][j]);
            }puts("");
        }
        return 0;
    }
    
    
    /*
    int dev(int sx, int sy, int ex, int ey) {
        int cnt = 0;
        for(int i = 0; i < 4; ++i) {
            int x = sx + dx[i], y = sy + dy[i]; // 与起点相邻的点 
            if(check(x, y)) continue;
            if(col[ex][ey] != col[x][y]) cnt++; 
        }
        for(int i = 0; i < 4; ++i) {
            int x = ex + dx[i], y = ey + dy[i]; // 与终点相邻的点 
            if(check(x, y)) continue;
            if(col[sx][sy] != col[x][y]) cnt++; 
        }
    //    cout<<"cnt1 :"<< cnt<<" ";
        return cnt;
    } 
    
    int Calc() {
        int cnt = 0;
        for(int i = 2; i <= n; ++i) 
            for(int j = 1;j <= m; ++j) 
                if(col[i][j] != col[i - 1][j]) cnt++;
        for(int i = 1; i <= n; ++i) 
            for(int j = 2;j <= m; ++j) 
                if(col[i][j] != col[i][j - 1]) cnt++;
        return cnt;
    }
    */
    

    几个练习题

    P2210 Haywire

    P2538 [SCOI2008]城堡

    SP4587 FENCE3 - Electric Fences

  • 相关阅读:
    第一次MVC记录
    Treeview绑定以及添加多选功能
    BindingSource的使用
    WPF实现环(圆)形进度条
    WPF实现环(圆)形菜单
    WPF实现音乐字幕动画
    WPF加载高德地图
    WPF实现Android(3D)菜单翻转动画
    WPF实现头像裁剪
    WPF PointAnimationUsingKeyFrames 动画
  • 原文地址:https://www.cnblogs.com/Silymtics/p/14880620.html
Copyright © 2011-2022 走看看