模拟退火
前言骚话
模拟退火其实不怎么难,我第一次学的时候也被那些可怕的物理知识吓到了,但是不需要担心,我们是oi考试,不是物竞,只需要关于oi的部分就可以了。(其实我就看了别人代码,就差不多知道了模拟退火是个什么东西了,逃)
爬山算法
爬山算法是模拟退火的前置算法,本质就是一个非常笨的贪心,如上图,如果我们问题是要找最高的点,那么我们从箭头的方向寻找答案,每次判断两边的是否都比当前的答案更劣,图中的答案就是A点,因为在这个区间内没有比A更高的节点。虽然这个算法很明显是错误的,但是却诞生出了一个怪物一般的模拟退火。
模拟退火
简介(来自百度百科)
模拟退火算法(Simulate Anneal,SA)是一种通用概率演算法,用来在一个大的搜寻空间内找寻命题的最优解。模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的。V.Černý在1985年也独立发明此演算法。模拟退火算法是解决TSP问题的有效方法之一。
模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法是一种通用的优化算法,其物理退火过程由加温过程、等温过程、冷却过程这三部分组成。
流程图
人话解释
我们还是爬山算法,还是判断当前答案是不是更优,但是我们加的操作是:如果不是更优,那么我们就来一个继续迭代的概率,这个概率会随着整体的程序越变越小。
做一个比喻:
爬山算法是一个瞎的兔子,要跳到最高的山峰,他跳着跳着发现路开始向下了,他就不跳了。
模拟退火是一个喝醉的兔子,要跳到最高的山峰,他一直在跳,越跳越清醒,最后调到了最高的地方。
物理和oi的关系(纯属瞎扯不想看的可以直接跳到下一部分)
模拟退火的原理也和金属退火的原理近似:将热力学的理论套用到统计学上,将搜寻空间内每一点想像成空气内的分子;分子的能量,就是它本身的动能;而搜寻空间内的每一点,也像空气分子一样带有“能量”,以表示该点对命题的合适程度。演算法先以搜寻空间内一个任意点作起始:每一步先选择一个“邻居”,然后再计算从现有位置到达“邻居”的概率。
原理
一开始一个处于高温的物体,现在要让物体降温,但是实际上,过快地降温使得物体来不及有序地收缩,难以形成结晶。而结晶态,才是物体真正内能降到最低的形态。
那么我们就要慢慢降温,也就是退火,物体的每一个粒子都有足够的时间找到自己的最佳位置并紧密有序地排列。开始温度高的时候,粒子活跃地运动并逐渐找到一个合适的状态。在这过程中温度也会越降越低,温度低下来了,那么粒子也渐渐稳定下来,相较于以前不那么活跃了。这时候就可以慢慢形成最终稳定的结晶态了。
退火的参数
\(T\)----表示当前的温度。
\(ΔT\)----表示温度的变化量,一般是0.95到0.99,让物体慢慢降温。
\(x\)----当前的解
\(Δx\)----解变动大小
\(x'\)----当前的解
\(Δf\)表示答案的差值
那么每次我们就比较差值随机出来的答案是不是在范围内,越到后面随机的答案就越小。
退火模板
提供一个比较好的固定随机种子,15346301,是ouhuang mod 66666666,很多人都用自己的生日,我之前也是用生日的。
例题
[JSOI2004]平衡点 / 吊打XXX
#include <bits/stdc++.h>
#define db double
#define N 1005
using namespace std;
struct node {
db x, y, w, tot;
}a[N], ans;
int n;
db sqr(db x) {
return x * x;
}
db dist(node a, node b) {
return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));
}
db calc(node p) {
db res = 0;
for (int i = 1; i <= n; i ++)
res += a[i].w * dist(p, a[i]);
return res;
}
db Rand() {
return rand() % 10000 / 10000.0;
}
void SA(db T) {
node now = ans, nw;
double delta;
while (T > 1e-4) {
nw.x = now.x + T * (2 * Rand() - 1);
nw.y = now.y + T * (2 * Rand() - 1);
nw.tot = calc(nw);
if (nw.tot < now.tot || exp((now.tot - nw.tot) / T) > Rand()) now = nw;
if (nw.tot < ans.tot) ans = nw;
T = T * 0.98;
}
for (int i = 1; i <= 1500; i ++) {
nw.x = ans.x + T * (2 * Rand() - 1);
nw.y = ans.y + T * (2 * Rand() - 1);
nw.tot = calc(nw);
if (nw.tot < ans.tot) ans = nw;
}
}
int main() {
srand(15346301);
scanf("%d", &n);
for (int i = 1; i <= n; i ++) {
scanf("%lf%lf%lf", &a[i].x, &a[i].y, &a[i].w);
ans.x += a[i].x;
ans.y += a[i].y;
}
ans.x /= n;
ans.y /= n;
ans.tot = calc(ans);
for (int i = 1; i <= 150; i ++) SA(100000);
printf("%.3lf %.3lf\n", ans.x, ans.y);
return 0;
}
[TJOI2010]分金币
#include <bits/stdc++.h>
#define ms(a, b) memset(a, b, sizeof(a))
#define inf 0x3f3f3f3f
#define db double
using namespace std;
inline char gc() {
static char buf[1 << 16], *S, *T;
if (S == T) {
T = (S = buf) + fread(buf, 1, 1 << 16, stdin);
if (T == S) return EOF;
}
return *S ++;
}
template <typename T>
inline void read(T &x) {
T w = 1;
x = 0;
char ch = gc();
while (ch < '0' || ch > '9') {
if (ch == '-') w = -1;
ch = gc();
}
while (ch >= '0' && ch <= '9') x = (x << 1) + (x << 3) + (ch ^ 48), ch = gc();
x = x * w;
}
template <typename T>
void write(T x) {
if (x < 0) putchar('-'), x = -x;
if (x > 9) write(x / 10);
putchar(x % 10 + 48);
}
#define N 105
int n, ans;
int a[N];
int calc() {
int res1 = 0, res2 = 0;
for (int i = 1; i <= n; i ++)
if (i <= (n + 1) / 2) res1 += a[i];
else res2 += a[i];
return abs(res1 - res2);
}
db Rand() {
return rand() % 10000 / 10000.0;
}
void SA(db T) {
while (T > 1e-3) {
int x = rand() % ((n + 1) / 2) + 1, y = rand() % ((n + 1) / 2) + ((n + 1) / 2);
if (x <= 0 || x > n || y <= 0 || y > n) continue;
swap(a[x], a[y]);
int res = calc();
if (ans > res) ans = res;
else if ((exp((1.0 * ans - 1.0 * res) / T)) <= Rand()) swap(a[x], a[y]);
T *= 0.98;
}
}
int main() {
srand(15346301);
int cas;
read(cas);
for (int _t = 1; _t <= cas; _t ++) {
read(n);
for (int i = 1; i <= n; i ++) read(a[i]);
ans = inf;
for (int i = 1; i <= 150; i ++) SA(10000);
printf("%d\n", ans);
}
return 0;
}
[HAOI2006]均分数据
洛谷传送门
题解:https://www.cnblogs.com/chhokmah/p/10529628.html
#include <bits/stdc++.h>
#define ms(a,b) memset(a, b, sizeof(a))
#define db double
using namespace std;
inline char gc() {
static char buf[1 << 16], *S, *T;
if (S == T) {
T = (S = buf) + fread(buf, 1, 1 << 16, stdin);
if (T == S) return EOF;
}
return *S ++;
}
template <typename T>
inline void read(T &x) {
T w = 1;
x = 0;
char ch = gc();
while (ch < '0' || ch > '9') {
if (ch == '-') w = -1;
ch = gc();
}
while (ch >= '0' && ch <= '9') x = (x << 1) + (x << 3) + (ch ^ 48), ch = gc();
x = x * w;
}
template <typename T>
void write(T x) {
if (x < 0) putchar('-'), x = -x;
if (x > 9) write(x / 10);
putchar(x % 10 + 48);
}
#define N 305
db ans = 1e30, ave = 0;
int sum[N], pos[N], a[N];
int n, m;
void SA(db T){
ms(sum, 0);
for (int i = 1; i <= n; i ++) {
pos[i] = rand() % m + 1;
sum[pos[i]] += a[i];
}
db res = 0;
for (int i = 1; i <= m; i ++)
res += (1.0 * sum[i] - ave) * (1.0 * sum[i] - ave);
while (T > 1e-4) {
int t = rand() % n + 1, x = pos[t], y;
if (T > 500) y = min_element(sum + 1, sum + 1 + m) - sum;
else y = rand() % m + 1;
if (x == y) continue;
db tmp = res;
res -= (sum[x] - ave) * (sum[x] - ave);
res -= (sum[y] - ave) * (sum[y] - ave);
sum[x] -= a[t], sum[y] += a[t];
res += (sum[x] - ave) * (sum[x] - ave);
res += (sum[y] - ave) * (sum[y] - ave);
if (res < tmp || rand() % 10000 <= T) pos[t] = y;
else sum[x] += a[t], sum[y] -= a[t], res = tmp;
ans = min(ans, res);
T *= 0.98;
}
}
int main() {
srand(15346301);
read(n); read(m);
for (int i = 1; i <= n; i ++) {
read(a[i]);
ave += 1.0 * a[i];
}
ave /= 1.0 * m;
for (int i = 1; i <= 1500; i ++) SA(10000);
printf("%.2lf\n", sqrt(ans / m));
return 0;
}