玄学乱搞算法——模拟退火，SA

\(\texttt{0x00:}\) 前言

在此之前只对模拟退火的大名有所耳闻，但并未在我的认知上激起太大的风浪，直到……

在外培的一场模拟赛上，队内大佬 yyc 在丝毫没有思路的情况下用 SA 骗了 70pts，赛后使得给我们上课的清华姚班老师惊掉下巴。

至此，在感叹 SA 的神力的同时，它也进入了我的学习计划中。

\(\texttt{0x01:}\) 介绍

退火原本是一个物理概念，是指一种热处理工艺，将固体材料加热到一定温度后，让其缓慢冷却的过程。

模拟退火算法（Simulate Anneal，SA）是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。模拟退火是由 S.Kirkpatrick，C.D.Gelatt 和 M.P.Vecchi 在1983年所发明的。V.Černý 在 1985 年也独立发明此演算法。模拟退火算法是解决 TSP 问题的有效方法之一。——百度百科

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

简单来讲，模拟退火是一种启发式搜索，它通过模拟物理降温的过程，先随机取一个起始点，然后随机地搜索其他点，若更优就跳过去，否则以一定的概率跳过去。由于模拟退火的过程中有概率接受较劣的解，所以它有一定概率找到全局最优解而不是受困于某一局部最优解。

\(\texttt{0x02:}\) 算法过程

首先设定几个参数：初温 \(T_0\)，末温 \(T_E\)，降温系数 \(k\)。

在算法开始时钦定温度 \(T = T_0\)，然后进行最优解的寻找，重复找最优解的过程并不断将温度乘上 \(k\) 直到 \(T = T_E\)。

具体来讲，在找最优解时的准则：如果新状态的解更优则修改答案，否则以一定概率接受新状态。

我们定义当前温度为 \(T\)，新状态 \(S'\) 与已知状态 \(S\)（新状态由已知状态通过随机的方式得到）之间的能量（值）差为 \(\Delta_E\)（\(\Delta_E\geqslant 0\)），则发生状态转移（修改最优解）的概率为：

\[\large P(\Delta_E) = \begin{cases} 1, \space\space\space\space\space\space\space\space\space\space\space\space\text{S' is better than S,}\\ \frac{-\Delta_E}{T},\space\space\space\space\space\space \text{otherwise.} \end{cases}\]

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

注意：模拟退火只能用来求解具有一定连续性的函数的最优解，这个连续性可以具体看题目中一个微小的策略改变是否导致了函数值的剧变，若导致了，则没多少连续性，否则可以认为有一定连续性。

模拟退火的过程如图所示：

\(\texttt{Code:}\)

void SA() {
    double T = 1e3;
    while(T > 1e-7) {
        //do something...
        ll now = calc(); //求当前的函数值
        ll dt = now - ans; //ans 是全局最优解
        if(dt < 0) ans = now; //更优则直接接受
        else if(exp(-(double)dt / T) * RAND_MAX < (double)rand()) 
			//不是更优则以一定概率接受
        T *= k;
    }
}

\(\texttt{0x03:}\) 一些例题

P1337 [JSOI2004] 平衡点 / 吊打XXX

题目大意

给定 \(n\) 的点表示每根绳子的固定点，每条绳子都系着一个重物，求出平衡点。

思路

这道题正解应该是三分套三分，但是不会，考虑乱搞。

根据基本的物理知识，平衡点会偏向重物所在的地方，事实上，假如钦定一个点 \(p\)，定义每条绳子的势能为点 \(p\) 到固定点的距离乘上所挂重物的重量，那么平衡点的势能应该是最小的，这个可以通过移动平衡点来证明确实很有道理。

同时在点 \(p\) 移动的距离很小的情况下，势能的变化也较小，所以函数具有连续性，可以用模拟退火来做。

这里我们顺便用当前温度代表当前的步长（随机的范围），便能轻松写出如下代码：

#include <cmath>
#include <ctime>
#include <iostream>

#define x first
#define y second

using namespace std;

const int N = 1010;
const double eps = 1e-9;
typedef pair<double, double> PDD;
int n;
struct node{
    int xx, yy, w;
}q[N];
double ans = 1e18;
int sumx, sumy;
PDD anss;

inline double rand(double l, double r) {
    return (double)rand() / RAND_MAX * (r - l) + l;
} 

inline double get_dist(PDD a, PDD b) {
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}

inline double calc(PDD now) {
    double res = 0;
    for(int i = 1; i <= n; i++)
        res += get_dist(now, {q[i].xx, q[i].yy}) * q[i].w;
    return res;
}

void SA() {
    PDD now(rand(-10000.0, 10000.0), rand(-10000.0, 10000.0));
    double T = 1e4;
    while(T > 1e-5) {
        PDD newp(now.x + (double)((rand() << 1) - RAND_MAX) * T, now.y + (double)((rand() << 1) - RAND_MAX) * T);
        double fut = calc(newp);
        double dt = fut - ans;
        if(dt < 0) ans = fut, anss = newp, now = newp;
        else if(exp(-dt / T) * RAND_MAX > (double)rand()) now = newp;
        T *= 0.95;
    }
}

int main() {
    srand((unsigned)time(0));
    scanf("%d", &n);
    for(int i = 1; i <= n; i++) {
        scanf("%d%d%d", &q[i].xx, &q[i].yy, &q[i].w);
    }
    SA();
    printf("%.3lf %.3lf", anss.x, anss.y);
    return 0;
}

虽然函数都写对了，但是这样连样例都过不了，怎么办？

其实有模拟退火的三大调整策略：

改进随机数的生成或是用贪心等策略优化找最优解的过程；
在时间允许的条件下，可以适当调高降温系数 \(k\)，调整初温或末温，还可以多跑几遍，甚至卡时；~~（调参是模拟退火最好玩的一环，不爽不要玩！）~~
采用精度更高的数据类型，比如 long double。（其实这条没什么用，主要看前两条）

由于温度和步长有关，所以这道题可以将末温调低一点，同时时间很充裕，所以可以调大降温系数和多跑几遍。

这样我们就 AC 了这道题。

其实还可以将起点设为固定点中的平均值，这样更容易跑出最优解。

P2503 [HAOI2006] 均分数据

题目大意

给定一个长为 \(n\) 的序列 \(\{a_n\}\)，先要将其分成 \(m\) 段，使得各组数字之和的标准差最小。

思路

因为数据范围比较小，所以其实各种乱搞都能过（

考虑模拟退火。

比较朴素的一种想法是，每次找最优解时就随机将这 \(n\) 个数分到 \(m\) 个组里去，然后再计算标准差，更新最优解。

但是这样不容易搜到最优解，因为可能将许多大数丢到同一个集合去，所以很多情况都没有用。

综上，我们可以加一个贪心优化，具体地，先随机打乱原序列 \(a\)，然后每次将一个数加到总和最小的那个组里去，这样就能更快搜到最优解。

\(\texttt{Code:}\)

#include <cmath>
#include <ctime>
#include <cstring>
#include <iostream>
#include <algorithm>

using namespace std;

const int N = 25, M = 10;

int n, m;
int a[N], s[M];
double ans = 1e9;

double calc() {
    memset(s, 0, sizeof s);
    for(int i = 1; i <= n; i++) {
        int pos = 1;
        for(int j = 2; j <= m; j++)
            if(s[j] < s[pos]) pos = j;
        s[pos] += a[i];
    }

    double avg = 0;
    for(int i = 1; i <= m; i++) avg += 1.0 * s[i] / m;
    double res = 0;
    for(int i = 1;  i <= m; i++) res += (s[i] - avg) * (s[i] - avg);
    res = sqrt(res / m);
    ans = min(ans, res);
    return res;
}

void SA() {
    random_shuffle(a + 1, a + n + 1);
    double T = 1e4;
    while(T > 1e-6) {
        int x = rand() % n + 1, y = rand() % n + 1;
        if(x == y) continue;
        // printf("%d %d\n", x, y);
        double now = calc();
        swap(a[x], a[y]);
        double fut = calc();
        double dt = fut - now;
        if(exp(-dt / T) * RAND_MAX < (double)rand())
            swap(a[x], a[y]);
        T *= 0.99;
    }
}

int main() {
    scanf("%d%d", &n, &m);
    for(int i = 1; i <= n; i++) scanf("%d", &a[i]);
    for(int i = 0; i < 10; i++) SA();
    printf("%.2lf", ans);
    return 0;
}

P4044 [AHOI2014/JSOI2014] 保龄球

题目大意

一场保龄球比赛共有 \(N\) 轮，每一轮都会有 \(10\) 个木瓶放置在木板道的另一端。

每一轮中，选手都有两次投球的机会来尝试击倒这 \(10\) 个木瓶。

对于每一次投球机会，选手投球的得分等于这一次投球所击倒的木瓶数量。

选手每一轮的得分是他两次机会击倒全部木瓶的数量。

对于每一个轮次，有如下三种情况：

“全中”：若选手第一次尝试就击倒了全部 \(10\) 个木瓶，则“全中”。在一个“全中”轮中，不需要再进行第二次投球尝试。作为奖励，选手的下一轮的得分将会被乘 \(2\) 计入总分。

“补中”：若选手两次尝试击倒了 \(10\) 个木瓶，则“补中”。作为奖励，选手的下一轮中的第一次尝试的得分将会被乘以 \(2\) 计入总分。

“失误”：若选手未能通过两次尝试击倒全部的木瓶，则“失误”。失误没有奖励。

此外，如果第 \(N\) 轮是“全中”，那么选手将一共进行 \(N+1\) 轮比赛。

附加轮的规则只执行一次。

最后，选手的总得分就是附加轮规则执行过，并且分数按上述规则加倍后的每一轮分数之和。

JYY 刚刚进行了一场 \(N\) 个轮次的保龄球比赛，但是，JYY 非常不满意他的得分。

JYY 希望通过更改每一轮成绩的顺序来最大化他的得分。

当然了，JYY 不希望做的太假，他希望保证重新排列之后，所需要进行的轮数和重排前所进行的轮数是一致的：

比如如果重排前 JYY 在第 \(N\) 轮打出了“全中”，那么重排之后，第 \(N\) 轮还得是“全中”以保证比赛一共进行 \(N+1\) 轮；同样的，如果 JYY 第 \(N\) 轮没有打出“全中”，那么重排过后第 \(N\) 轮也不能是全中。

求 JYY 可以得到的最高的分数。

（题目太长了，已经尽力简化了）

思路

这道题的数据范围同样较小，所以可以乱搞（

考虑模拟退火。

先看一下有没有 \(n + 1\) 轮，然后每次随机生成一个排列，退火过程中随机交换此排列中的两个位置，再按规则计算分数，更新最大值即可。

\(\texttt{Code:}\)

#include <cmath>
#include <ctime>
#include <iostream>

#define x first
#define y second

using namespace std;

const int N = 55;
typedef pair<int, int> PII;
int n, nn;
PII q[N];
int ans;

inline void swap(PII &x, PII &y) {
    PII t = y;
    y = x;
    x = t;
}

inline int calc() {
    int res = 0;
    for(int i = 1; i <= nn; i++) {
        res += q[i].x + q[i].y;
        if(i <= n) {
            if(q[i].x == 10) res += q[i + 1].x + q[i + 1].y;
            else if(q[i].x + q[i].y == 10) res += q[i + 1].x;
        }
    }
    ans = max(ans, res);
    return res;
}

void SA() {
    double T = 1e6;
    while(T > 1e-6) {
        int x = rand() % nn + 1, y = rand() % nn + 1;
        int now = calc();
        swap(q[x], q[y]);
        if(n + (q[n].x == 10) == nn) {
            int fut = calc();
            int dt = fut - now;
            if(exp(dt / T) < (double)rand() / RAND_MAX)
                swap(q[x], q[y]);
        }
        else swap(q[x], q[y]);
        T *= 0.996;
    }
}

int main() {
    srand((unsigned)(time(0)));
    srand((unsigned)rand() * rand()); //玄学种子
    scanf("%d", &n);
    for(int i = 1; i <= n; i++) scanf("%d%d", &q[i].x, &q[i].y);
    if(q[n].x == 10) nn = n + 1, scanf("%d%d", &q[nn].x, &q[nn].y);
    else nn = n;
    for(int i = 1; i < 3; i++) SA();
    printf("%d\n", ans);
    return 0;
}

posted @ 2024-08-27 21:28 Brilliant11001 阅读(5) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

Brilliant11001

一只会咬人的氯化氢~

玄学乱搞算法——模拟退火，SA

\(\texttt{0x00:}\) 前言

\(\texttt{0x01:}\) 介绍

\(\texttt{0x02:}\) 算法过程

\(\texttt{0x03:}\) 一些例题

P1337 [JSOI2004] 平衡点 / 吊打XXX

题目大意

思路

P2503 [HAOI2006] 均分数据

题目大意

思路

P4044 [AHOI2014/JSOI2014] 保龄球

题目大意

思路

公告