模拟退火学习笔记

· · 个人记录

1 引入

1.1 算法相关

笔者认为模拟退火是一种以人类思维启发式的、不会回溯的搜索，通常这种算法都不能通过一道题，因为其正确性无法保证的基础上还引入了大量的随机化，但是在实际比赛中对于一些数据较小、用其它算法很不好写的题目模拟退火可以得到一个还不错的分数。

2 前置知识

2.1 爬山算法

爬山算法是不断地寻找当前可转移的局部最优解并进行转移，若所有可转移的局部解都劣于当前的解，则算法结束，并认定当前解为最优解。对于单峰函数而言，正确性显然，但如果全局最优解需要从当前局部最优解转移到一个更劣解再转移，这样的算法就会像贪心一样错误。

3 算法讲解

3.1 算法流程

模拟退火大致为爬山算法的改进版。首先会找到一个初始状态，接下来对于所有状态，我们会有概率的选择转移或不转移。设当前状态为 S，选择转移或否的状态为 S'，转移概率如下：

P= \begin{cases} 1 &\mathrm{S'\ is\ better\ than\ S} \\ e^{\frac{-\Delta E}{t}} &\mathrm{otherwise} \end{cases}

其中 \Delta E 是 S 和 S' 的差值，t 是温度，温度就是模拟退火的一个核心概念了。

一开始，温度通常非常大，则在上述转移的概率中，第二种情况的概率也比较大。随着时间的转移，温度将会逐渐缩小（即降温），那么第二种情况的概率也会减小。通常情况下，我们每一次转移后就会将 t 乘上一个略比 1 小的数（参考 [0.985,0.999]），当 t 小到一定程度（参考 10^{-10}）就会停止降温，故而停止算法。

没理解为什么？我们将这种算法转移到生活实际上：

Pretharp 在玩游戏，这个游戏是一个战略游戏，你可以选择保持当前的情况，也可以选择放手一搏，这样可能为 Pretharp 带来更好的情况，也有可能把 Pretharp 逼入绝境。

一开始，Pretharp 非常爱冒险，会随便挑选一个未来的局面作为目标。如果当前可以使局面变得更好，Pretharp 一定会这么干，但是如果会使局面更差，Pretharp 也会有很大的概率去选择这么做（温度高）。

随着时间的转移，Pretharp 开始变得谨慎，如果选择一个未来的局面作为目标会变得更好，Pretharp 还是会这么干，但是若会变得更差，Pretharp 会很谨慎的考虑去选择这么做（温度低）。> 当 Pretharp 已经不想再冒险时（温度到达很低值），游戏结束。

放一张图，图来自 Simulated annealing - Wikipedia，图链来自于 OI Wiki 的图：

可以看到，随着温度降低当前状态的跳跃程度越来越小，最终锁定一个最大值。

3.2 注意事项

因为模拟退火的正确性无法保证，因此，我们可以记录算法过程中遇到的最优答案作为答案。

其次，每一次运行算法的结果可能不一样，因此我们可以在不超时的情况下多次运行该算法以得到更高的概率来获得正确答案。考虑到确定一个确切的运行算法次数较难且有大概率超时，我们通常使用以下代码卡时：

while((double)clock() / CLOCKS_PER_SEC < TimeLimit) { // TimeLimit 是一个变量，是一个比本题实现略小的数（单位秒）。
    /*模拟退火算法本身*/
}

最后，对于求出 e^x，在 C++ 下可以直接使用 exp(x) 求出其值。

3.3 参考模板

考虑到模拟退火每一次算法的代码大体相似，这里提供一种模拟退火的参考实现方式（CPP）。

const double dcr = 0.988;
int ans, ans_pos, now;
double calc(type1 name1, ...) { // 用于求出当前局面的权值的函数。 
    code1;
    code2;
    ...
    return result;
}
#define ttSA SimulateAnneal
void SimulateAnneal() { // 模拟退火本体。 
    double t = 2000; // 初始化温度。 
    now = ans_pos; // 用于多次执行模拟退火，每一次开始时从最优状态开始转移。 
    while(t > 1e-14) { // 结束条件通常设置为温度过低。 
        double nxt = now + ...; // 随机出下一个状态。 
        double delta = calc(nxt) - ans; // 求出 delta E，即下一个状态和当前状态的差值。 
        if(delta > 0) { // delta 大于 0（也有可能是其它情况），说明下一个状态比当前状态更有，进行转移。 
            now = nxt, ans_pos = nxt; // 更新当前状态及最优状态。 
            ans = calc(ans_pos); // 求出当前最优答案。 
        } else if(exp(-delta / t) * RAND_MAX > rand()) { // 否则，概率转移。 
            now = nxt;
        }
        t *= dcr; // 降温。 
    }
}
signed main() {
    srand(time(NULL)); // 重置随机种子。 
    /* 输入部分 */
    ans_pos = ..., ans = calc(ans_pos); // 一般初始会将 ans_pos 设置为一个较为接近于平均值的值，也可以完全随机。 
    while((double)clock() / CLOCKS_PER_SEC < 0.95) { // 卡时执行模拟退火。 
        ttSA();
    }
    /* 输出部分 */ 
    return 0;
}

上述代码中很多参数可用于更改，代码仅供参考。

4 例题

4.1 例题 1~2

例题 1：[JSOI2004] 平衡点 / 吊打XXX

题目大意：

给定 n 个点的坐标 (x_i,y_i) 和权值 w_i，请求出一个点，使得 \sum\limits_{i=1}^n dist_i \times w_i 最小，其中 dist_i 是第 i 个点与你选择的点的欧几里得距离。

分析：

考虑模拟退火。不难发现，本题答案的函数跨度不大，我们可以理想出爬山算法的模型，在一个点附近有一个点答案很优（即山顶），所以我们可以使用模拟退火来寻找答案。在本题中，先找到一个初始状态，之后以这个点为中心向四周随机一个点作为下一个状态，且下一个状态和当前状态的距离会逐渐缩小。其余部分在之前的算法讲解基本都有介绍，可以理解为模板。

参考代码：

#include<bits/stdc++.h>
using namespace std;
const double dcr = 0.998;
const int N = 1010;
int n;
double nowx, nowy, ans, ansx, ansy;
struct Thing {
    int x, y, w;
} a[N];
double dist(double xa, double ya, double xb, double yb) {
    return sqrt((xa - xb) * (xa - xb) + (ya - yb) * (ya - yb));
}
double calc(double x, double y) {
    double res = 0;
    for(int i = 1; i <= n; i ++) {
        res += dist((double)a[i].x, (double)a[i].y, x, y) * (double)a[i].w;
    }
    return res;
}
#define ttSA SimulateAnneal
void SimulateAnneal() {
    double t = 2000;
    nowx = ansx, nowy = ansy;
    while(t > 1e-14) {
        double nxtx = nowx + t * (rand() * 2 - RAND_MAX);
        double nxty = nowy + t * (rand() * 2 - RAND_MAX);
        double delta = calc(nxtx, nxty) - ans;
        if(delta < 0) {
            nowx = nxtx, nowy = nxty;
            ansx = nowx, ansy = nowy, ans = calc(ansx, ansy);
        } else if(exp(-delta / t) * RAND_MAX > rand()) {
            nowx = nxtx, nowy = nxty;
        }
        t *= dcr;
    }
}
signed main() {
    srand(time(NULL));
    cin >> n;
    for(int i = 1; i <= n; i ++) {
        cin >> a[i].x >> a[i].y >> a[i].w;
        ansx += a[i].x, ansy += a[i].y;
    }
    ansx /= n, ansy /= n, ans = calc(ansx, ansy);
    while((double)clock() / CLOCKS_PER_SEC < 0.95) {
        ttSA();
    }
    return printf("%.3lf %.3lf\n", ansx, ansy), 0;
}

例题 2：[HAOI2006] 均分数据

题目大意：

给定 n 个元素，每个点有一个权值 a_i，你需要将这些元素分为 m 组，其中第 j 组中的元素和为 b_i，使得 \sqrt{\frac{1}{m} \sum\limits_{i=1}^m (avg - b_i)^2} 最小，其中 avg=\frac{1}{m} \sum\limits_{i=1}^{m} b_i。

分析：

还是同样的套路，初始状态可以将每个元素随机分组，之后转移每次选择一个随机的元素，将它移动到一个随机的组别。值得一提的是，模拟退火得不到满分也是很正常的事，对于每一道题可能需要适时修改各类参数。

参考代码：

#include<bits/stdc++.h>
using namespace std;
const int N = 25;
const double dcr = 0.996;
int n, m, a[N], sum[N], from[N];
double avg, ans;
double calc() {
    double res = 0;
    for(int i = 1; i <= m; i ++) {
        res += 1.0 * (avg - sum[i]) * (avg - sum[i]);
    }
    return res / m;
} 
void modify(int x, int y) {
    sum[from[x]] -= a[x];
    from[x] = y;
    sum[from[x]] += a[x];
} 
#define ttSA SimulateAnneal
void SimulateAnneal() {
    double t = 3000;
    while(t > 1e-11) {
        int p1 = rand() % n + 1, p2 = rand() % m + 1, p3 = from[p1];
        modify(p1, p2);
        double now = calc(), delta = now - ans;
        modify(p1, p3);
        if(delta < 0) {
            ans = now;
            modify(p1, p2);
        } else if(exp(-delta / t) * RAND_MAX > rand()) {
            modify(p1, p2);
        }
        t *= dcr;
    }
}
signed main() {
    srand(time(NULL));
    cin >> n >> m;
    for(int i = 1; i <= n; i ++) {
        cin >> a[i], avg += a[i];
        from[i] = rand() % m + 1, sum[from[i]] += a[i];
    }
    avg /= m, ans = calc();
    while((double)clock() / CLOCKS_PER_SEC < 0.9) {
        ttSA();
    }
    return printf("%.2lf", sqrt(ans)), 0;
}

5 习题

5.1 对于例题 1~2 的习题

  1. [SCOI2008] 城堡

  2. [JSOI2016] 炸弹攻击1

  3. [JSOI2008] 球形空间产生器

  4. [洛谷 P3936] Coloring

更多习题参见 [洛谷] 模拟退火系列。