模拟退火
Tangninghaha
·
·
算法·理论
简介
模拟退火算法（Simulated Annealing）是一种基于概率论的随机化算法，是对爬山算法的改进，它通过模拟退火过程来寻找全局最优解。
一般 OI 比赛中可以使用模拟退火获得更高的部分分，在部分题目可能由于数据强度不够，甚至可能通过。
算法描述
模拟退火算法随机求得一个新的解，并与当前解进行比较，如果新的解更优，则接受，否则以一定概率接受较差的解。
该算法需要三个参数：温度参数 T_0，降温参数 \Delta t，终止温度 T_k。
我们记录当前的温度 T=T_0，在当前温度下进行一次 转移，得到新的温度 T=T\times \Delta t，当 T<T_k 时终止退火。
转移是更新答案的过程。具体而言，根据当前的温度，在当前状态的基础上随机得到一个新状态，求得两个状态的答案差值 \Delta E(\Delta E\ge 0)。如果新状态更优，那么将新状态作为当前状态，否则以概率：
P(\Delta E)=e^{\frac{-\Delta E}{T}}
将新状态作为当前状态。
这样描述可能很抽象，下面用几个例题来进行分析。
例题
P1337 [JSOI2004] 平衡点 / 吊打XXX
题目描述：
求给定 n 个点的带权费马点。具体地，若给定点点权为 v_i，到费马点的距离为 dis_i，那么代价为 \sum v_i\times dis_i，找到一个点最小化代价。输出点的坐标。
这题显然具有单调性，可以通过三分解决。但是同样可以使用模拟退火通过。
下面结合代码讲解如何使用模拟退火：
```c++
const int N = 1005;
struct Point {
int x, y, w;
} p[N];
const double D = 0.97;
const double MIN_ = -1e4, MAX_ = 1e4;
double calc(double x, double y) { // 对于一个坐标计算答案
double res = 0;
for (int i = 1; i <= n; ++i)
res += p[i].w *
sqrt((p[i].x - x) * (p[i].x - x) +
(p[i].y - y) * (p[i].y - y));
return res;
}
void simulateAnneal() {
double t = 1e10, t0 = 1e-10; // 定义初始温度和终止温度
random_device seed; // 使用 random_device 生成基于硬件的真随机数作为种子
mt19937 rad(seed()); // 定义随机数生成器
uniform_real_distribution<> random_coord(MIN_, MAX_); // 定义坐标随机生成器
uniform_real_distribution<> random_posb(0, 1); // 定义概率生成器
while (t >= t0) {
// 根据当前状态和温度随机新状态
double new_x = cur_x + random_coord(rad) * t, new_y = cur_y + random_coord(rad) * t;
// 计算新的答案
double new_ans = calc(new_x, new_y);
if (new_ans < cur_ans) {
// 如果新答案更优，直接更新
cur_ans = new_ans;
cur_x = new_x;
cur_y = new_y;
} else if (random_posb(rad) <= exp((cur_ans - new_ans) / t)) {
// 否则按照上述概率接受
cur_ans = new_ans;
cur_x = new_x;
cur_y = new_y;
}
t *= D; // 降温
}
}
```
其中关于随机数的部分是值得考究的：
`random_device` 是基于硬件的随机数发生器，生成的随机数常用作伪随机数种子。实际上如此并不是必要的，可以直接用你喜欢的数字作为种子。
`mt19937` 是基于 32 位梅森缠绕器的伪随机数发生器，生成随机数效果较好，常被使用。
`uniform_real_distribution(l,r)` 用于生成一个位于 $[l,r)$ 的随机数，是左闭右开的。
更详细的内容可见 [随机数与随机函数](./random) 一文。
### [P2503 [HAOI2006] 均分数据](https://www.luogu.com.cn/problem/P2503)
将 $n$ 个正整数分成 $m$ 组，每组的权值为该组的数字和，希望最小化权值的方差。求最小方差的值。
$m\le n\le 20,2\le m\le 6,a_i\in[1,50]$。
这题应该可以使用状压 DP 进行。考虑使用模拟退火，每一次随机得到一个数组，规定每一组在数组 $a$ 中都是连续的一段，使用简单 DP 可以在 $O(n^3)$ 的时间内计算出答案。每一次新状态由原状态随机交换位置得到。类似地进行模拟退火：
```c++
const int N = 21;
int a[N];
int n, m;
void upd_min(double &x, double y) {
x = min(x, y);
}
double calc(int *a) { // 使用 DP 计算最小方差
static int sum[N];
static double f[N][N];
for (int i = 0; i <= n; ++i)
for (int j = 0; j <= m; ++j)
f[i][j] = 1.0 / 0.0;
f[0][0] = 0;
for (int i = 1; i <= n; ++i)
sum[i] = sum[i - 1] + a[i];
double aver = (double)sum[n] / m;
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= m; ++j) {
for (int k = 0; k < i; ++k) {
upd_min(f[i][j], f[k][j - 1] + (sum[i] - sum[k] - aver) * (sum[i] - sum[k] - aver));
}
}
}
return sqrt(f[n][m] / m);
}
const double D = 0.915;
double ans;
void simulateAnnel() {
double t = 100000, t0 = 1e-10;
static int cur[N], b[N];
double cur_ans = ans;
for (int i = 1; i <= n; ++i) cur[i] = a[i];
random_device seed;
mt19937 rad(seed());
uniform_real_distribution<> posb(0, 1);
while (t >= t0) {
for (int i = 1; i <= n; ++i) b[i] = cur[i];
swap(b[rad() % n + 1], b[rad() % n + 1]); // 随机交换得到新状态
double new_ans = calc(b);
if (new_ans < cur_ans) {
for (int i = 1; i <= n; ++i) cur[i] = b[i];
cur_ans = new_ans;
if (new_ans < ans) {
for (int i = 1; i <= n; ++i) a[i] = cur[i];
ans = new_ans;
}
} else if (posb(rad) <= exp(cur_ans - new_ans) / t) {
for (int i = 1; i <= n; ++i) cur[i] = b[i];
cur_ans = new_ans;
}
t *= D;
}
}
```
注意到进行一次模拟退火不一定得到最优解，可以多次进行模拟退火提高正确率。为了避免超时，可以使用名为 **卡时** 的技巧，在规定时间内尽可能多次进行模拟退火。
```c++
while ((double)clock() / CLOCKS_PER_SEC <= 0.9) {
simulateAnnel();
}
```
这里需要包含 `ctime` 头文件，clock() 会返回从程序开始到当前的时间，`CLOCKS_PER_SEC` 是系统定义的一个宏，这个告诉你一秒包含多少个 CLOCK 的值。`0.9` 为一个略小于时限的数字。
注意实际用时是 `0.9` 加上单次运行退火的时间，因为时间不超过 0.9 都会进行退火，请估计单次退火时间，以免超时。
## 习题
### [P5544 [JSOI2016] 炸弹攻击1](https://www.luogu.com.cn/problem/P5544)
题目描述：
给定平面上 $N$ 个圆和 $M$ 个点。请你选择一个点和不超过 $R$ 的半径作圆，在这个圆与给定圆不相交的情况下，求最多能够覆盖多少个点。
$0\le N\le 10,0<M\le 10^{3}$，$1\le R,r_i\le 2\times 10^{4}$，所有坐标绝对值不超过 $2\times 10^{4}$。
```c++
#include <bits/stdc++.h>
using namespace std;
template <typename T>
void read(T &r) {
r = 0;
int f = 1;
char c = getchar();
for (; c < '0' || c > '9'; c = getchar())
if (c == '-') f = -1;
for (; c >= '0' && c <= '9'; c = getchar())
r = (r << 1) + (r << 3) + (c ^ 48);
r *= f;
}
template <typename First, typename... Rest>
void read(First &first, Rest &...rest) {
read(first);
read(rest...);
}
template <typename T>
void print(T v) {
if (v < 0) putchar('-'), v = -v;
static int sta[40], top;
top = 0;
do {
sta[++top] = v % 10;
v /= 10;
} while (v);
while (top) putchar(sta[top--] ^ 48);
}
const int N = 2e4 + 5;
struct Bulid {
int x, y, r;
}b[N];
struct Enemy {
int x, y;
}p[N];
double dist(double x1, double y1, double x2, double y2) {
return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
}
int n, m, r;
double calc(double x, double y) {
double rr = r;
for (int i = 1; i <= n; ++i)
rr = min(rr, dist(x, y, b[i].x, b[i].y) - b[i].r);
if (rr < 0) return 0;
int res = 0;
double minn = 1.0 / 0.0;
for (int i = 1; i <= m; ++i) {
double d = dist(x, y, p[i].x, p[i].y);
if (d > rr) minn = min(d - rr, minn);
res += (d <= rr + 1e-10);
}
if (res) return res;
return res - minn;
}
int calc2(double x, double y) {
double rr = r;
for (int i = 1; i <= n; ++i)
rr = min(rr, dist(x, y, b[i].x, b[i].y) - b[i].r);
if (rr < 0) return 0;
int res = 0;
for (int i = 1; i <= m; ++i)
res += (dist(x, y, p[i].x, p[i].y) <= rr + 1e-10);
return res;
}
double ans_x, ans_y, ans;
void simulateAnnel() {
double cur_x = ans_x, cur_y = ans_y, cur_ans = ans;
double t = 1e10, t0 = 1e-15;
const double D = 0.996;
random_device seed;
mt19937 rad(seed());
uniform_real_distribution<> diff(-100, 100);
uniform_real_distribution<> posb(0, 1);
while (t >= t0) {
double new_x = cur_x + diff(rad) * t;
double new_y = cur_y + diff(rad) * t;
double new_ans = calc(new_x, new_y);
if (new_ans > cur_ans) {
cur_ans = new_ans;
cur_x = new_x;
cur_y = new_y;
if (new_ans > ans) {
ans = new_ans;
ans_x = new_x;
ans_y = new_y;
}
} else if (posb(rad) <= exp((new_ans - ans) / t)) {
cur_ans = new_ans;
cur_x = new_x;
cur_y = new_y;
}
t *= D;
}
}
int main() {
read(n, m, r);
for (int i = 1; i <= n; ++i)
read(b[i].x, b[i].y, b[i].r);
for (int i = 1; i <= m; ++i)
read(p[i].x, p[i].y);
ans_x = 0;
ans_y = 0;
ans = calc(ans_x, ans_y);
while ((double)clock() / CLOCKS_PER_SEC <= 0.9) {
simulateAnnel();
}
print(calc2(ans_x, ans_y));
fprintf(stderr, "\n\nTime Used: %.3lfs\n", (double)clock() / CLOCKS_PER_SEC);
}
```
### [P3878 [TJOI2010] 分金币](https://www.luogu.com.cn/problem/P3878)
题目描述：
给定 $n$ 个数字，将 $n$ 个数字分成两组，两组的数字数量相差不超过 $1$，求两组金币价值差的最小值。
$1\le n\le 30$，每一个数字 $v_i$ 满足 $1\le v_i\le 2^{30}$。
```c++
#include <bits/stdc++.h>
using namespace std;
const int N = 31;
int v[N];
long long ans;
int n;
long long calc(int *a) {
long long res = 0;
for (int i = 1; i <= n / 2; ++i) res += a[i];
for (int i = n / 2 + 1; i <= n; ++i) res -= a[i];
if (res < 0) res = -res;
return res;
}
void simulated_annealing() {
double t = 100, delta_t = 0.965, tk = 1e-10;
random_device seed;
mt19937 Rad(seed());
uniform_int_distribution<> rad(1, n);
uniform_real_distribution<> pos(0, 1);
long long cur_ans = ans;
while (t >= tk) {
static int cur_v[N];
for (int i = 1; i <= n; ++i) cur_v[i] = v[i];
swap(cur_v[rad(Rad)], cur_v[rad(Rad)]);
long long nxt_ans = calc(cur_v);
if (nxt_ans < cur_ans || pos(Rad) <= exp((cur_ans - nxt_ans) / t)) {
cur_ans = nxt_ans;
ans = min(ans, nxt_ans);
for (int i = 1; i <= n; ++i) v[i] = cur_v[i];
}
t = t * delta_t;
}
}
int main() {
int tt;
scanf("%d", &tt);
double ti = 0.97 / tt;
for (int t = 1; t <= tt; ++t) {
scanf("%d", &n);
for (int i = 1; i <= n; ++i) scanf("%d", &v[i]);
ans = calc(v);
while ((double)clock() / CLOCKS_PER_SEC <= ti * t) {
simulated_annealing();
}
printf("%lld\n", ans);
}
}
```