模拟退火
1.算法描述
1.1 概述
简单说,模拟退火是一种随机化算法,用于求函数的极值。当一个问题的方案数量极大(甚至是无穷的)而且不是一个单峰函数时,我们常使用模拟退火求解。
它与爬山算法最大的不同是,在寻找到一个局部最优解时,赋予了它一个跳出去的概率,也就有更大的机会能找到全局最优解。
在 OI 领域,对应的,每次随机出一个新解,如果这个解更优,则接受它,否则以一个与温度和与最优解的差相关的概率接受它。
以下图片来自:https://oi-wiki.org/misc/simulated-annealing/
1.2 相关参数
初始温度: T 0 T_0 T0
结束温度: T s T_s Ts
降温系数: Δ t \Delta t Δt。这样每次温度就是上次的温度乘上 Δ t \Delta t Δt
能量差: Δ E = f ( t n e w ) − f ( t n o w ) \Delta E=f(t_{new})-f(t_{now}) ΔE=f(tnew)−f(tnow),即新点的能量减去当前的能量(能量也就是函数值)
接受概率: P ( Δ E ) = e − Δ E t P(\Delta E)=e^{\frac{-\Delta E}{t}} P(ΔE)=et−ΔE,t为当前温度,这样保证当能量差小于0时,概率P是大于1的,也就是必然接受,当能量差大于0时,能量差越大越不容易接受,t越大越容易接受。针对具体问题为:
- 如果是求凹函数的最值,那么判断条件写为:
if (exp(-dt / t) > rand(0, 1)) cur = new; // 如果概率大于(0 ~ 1),那么就跳到新的点
else continue; // 否则不跳
- 如果是求凸函数的最值,那么判断条件写为:
if (exp(dt / t) > rand(0, 1)) cur = new; // 如果概率大于(0 ~ 1),那么就跳到新的点
else continue; // 否则不跳
1.3 技巧
1.由于较为玄学,所以需要多跑几次模拟退火:
for (int i = 0; i < 100; i ++ ) simulate_anneal();
2.更改随机种子
3.卡时间,例如小于0.8秒就一直跑模拟退火,充分利用测评时间:
while ((double)clock() / CLOCKS_PER_SEC < 0.8) simulate_anneal();
4.超时说明降温系数太大了,降温过程太慢,方法:可以把 Δ t \Delta t Δt 从0.999改为0.99
5.错误可以是因为降温太快,直接跳过了正解,此时可以把:可以把 Δ t \Delta t Δt 从0.99改为0.999
6.因为这个维数较小的时候,按照经验提高ΔT是要比跑多次效果好的
2.模板
2.1 模板一
#include <bits/stdc++.h>
using namespace std;
typedef pair<double, double> PDD;
const int N = 110;
int n;
PDD q[N];
double ans = 1e8;
double rand(double l, double r) {
return (double)rand() / RAND_MAX * (r - l) + l;
}
double get_dist(PDD a, PDD b) {
double dx = a.first - b.first;
double dy = a.second - b.second;
return sqrt(dx * dx + dy * dy);
}
double calc(PDD p) { // 计算取当前点时,到所有的点的距离
double res = 0;
for (int i = 0; i < n; i++) res += get_dist(p, q[i]);
ans = min(ans, res);
return res;
}
void simulate_anneal() { // 模拟退火
PDD cur(rand(0, 10000), rand(0, 10000)); // 先随机一个点
for (double t = 1e4; t > 1e-4; t *= 0.9) { // 按照降温系数进行降温
PDD new_point(rand(cur.first - t, cur.first + t), rand(cur.second - t, cur.second + t)); // 得到一个新的点
double dt = calc(new_point) - calc(cur); // 计算新的点和旧的点的差值
if (exp(-dt / t) > rand(0, 1)) cur = new_point; // 如果新的点更优,那么跳到新的点去
}
}
int main() {
scanf("%d", &n);
for (int i = 0; i < n; i++) scanf("%lf%lf", &q[i].first, &q[i].second);
for (int i = 0; i < 100; i++) simulate_anneal(); // 跑100次
printf("%.0lf\n", ans);
return 0;
}
2.2 模板二
这个模板利用了mt19937,较为精确,以后最好都用这个
#include <bits/stdc++.h>
using namespace std;
#define urd uniform_real_distribution
mt19937 rng((unsigned)time(NULL));
const int N = 1e5 + 5;
int n, m;
double c[N];
struct node {
double x, y;
} a[N];
double sq(double x) { return x * x; }
double getd(node p1, node p2) {
double dx = p1.x - p2.x;
double dy = p1.y - p2.y;
return sqrt(sq(p1.x - p2.x) + sq(p1.y - p2.y));
}
double getsum(node x) {
double re=0;
for (int i = 1; i <= n; i++) re+=getd(x, a[i]);
return re;
}
double rand(double l, double r) { //计算一个l到r的随机值
return (urd<>(l, r)(rng) - urd<>(l, r)(rng));
}
double sa() { //模拟退火
node p{0, 0};
double res = getsum(p);
for (double t = 1e4; t >= 1e-4; t *= 0.99) {
node np = {p.x + rand(0,10000)* t,p.y + rand(0,10000)* t};
double nw = getsum(np);
if (nw < res || urd<>(0, 1)(rng) <= exp((res - nw) / t))
p=np, res = nw;
}
return res;
}
int main() {
cin >> n;
for (int i = 1; i <= n; i++) cin >> a[i].x >> a[i].y ;
double ans = 1e18;
while ((double)clock() / CLOCKS_PER_SEC < 0.8)
ans = min(ans, sa());
printf("%.0lf", ans);
return 0;
}
3.典型例题
acwing3167. 星星还是树
题意: 在二维平面上有 n n n 个点,第 i i i 个点的坐标为 ( x i , y i ) (xi,yi) (xi,yi)。请你找出一个点,使得该点到这 n 个点的距离之和最小。该点可以选择在平面中的任意位置,甚至与这 n 个点的位置重合 1 ≤ n ≤ 100 , 0 ≤ x i , y i ≤ 10000 1≤n≤100,0≤xi,yi≤10000 1≤n≤100,0≤xi,yi≤10000
题解: 模拟退火。每次先随机化一个点,然后按照降温系数降低温度,每次跳转到一个新的点,如果新的点更优,那么就跳到这个点去。
代码:
#include <bits/stdc++.h>
using namespace std;
typedef pair<double, double> PDD;
const int N = 110;
int n;
PDD q[N];
double ans = 1e8;
double rand(double l, double r) {
return (double)rand() / RAND_MAX * (r - l) + l;
}
double get_dist(PDD a, PDD b) {
double dx = a.first - b.first;
double dy = a.second - b.second;
return sqrt(dx * dx + dy * dy);
}
double calc(PDD p) { // 计算取当前点时,到所有的点的距离
double res = 0;
for (int i = 0; i < n; i++) res += get_dist(p, q[i]);
ans = min(ans, res);
return res;
}
void simulate_anneal() { // 模拟退火
PDD cur(rand(0, 10000), rand(0, 10000)); // 先随机一个点
for (double t = 1e4; t > 1e-4; t *= 0.9) { // 按照降温系数进行降温
PDD new_point(rand(cur.first - t, cur.first + t), rand(cur.second - t, cur.second + t)); // 得到一个新的点
double dt = calc(new_point) - calc(cur); // 计算新的点和旧的点的差值
if (exp(-dt / t) > rand(0, 1)) cur = new_point; // 如果新的点更优,那么跳到新的点去
}
}
int main() {
scanf("%d", &n);
for (int i = 0; i < n; i++) scanf("%lf%lf", &q[i].first, &q[i].second);
for (int i = 0; i < 100; i++) simulate_anneal(); // 跑100次
printf("%.0lf\n", ans);
return 0;
}
acwing2424. 保龄球
题意: 给定n个击倒保龄球的操作,每次给定a和b,表示第一个保龄球击倒瓶子数目为a,第二个保龄球击倒的瓶子数目为b。计分方案如下:如果第i轮a=10,那么奖励下一轮的a和b得分;如果第i论的a + b = 10,那么奖励下一轮的a的得分。现在要求重新排列这n个击倒保龄球的操作,要求最终得分最高,求最高的得分是多少?
题解: 模拟退火。先随机得到的一段序列的操作顺序,然后每次随机选择a和b,表示交换a和b,如果交换完更优,那么进行交换。
代码:
#include <bits/stdc++.h>
using namespace std;
typedef pair<int, int> PII;
const int N = 55;
int n, m;
PII q[N];
int ans;
int calc() { // 计算按照当前排列的得分
int res = 0;
for (int i = 0; i < m; i++) {
res += q[i].first + q[i].second;
if (i < n) {
if (q[i].first == 10)
res += q[i + 1].first + q[i + 1].second;
else if (q[i].first + q[i].second == 10)
res += q[i + 1].first;
}
}
ans = max(ans, res);
return res;
}
void simulate_anneal() {
for (double t = 1e4; t > 1e-4; t *= 0.99) {
int a = rand() % m, b = rand() % m; // 随机找到2个交换
int x = calc(); // 计算原来的值
swap(q[a], q[b]); // 交换
if (n + (q[n - 1].first == 10) == m) { // 比如满足这个条件才能交换
int y = calc(); // 计算新的值
int delta = y - x; // 计算差值
if (exp(delta / t) < (double)rand() / RAND_MAX) swap(q[a], q[b]); // 如果不满足条件,再次交换,即跳到新的点
} else
swap(q[a], q[b]);
}
}
int main() {
cin >> n;
for (int i = 0; i < n; i++) cin >> q[i].first >> q[i].second;
if (q[n - 1].first == 10)
m = n + 1, cin >> q[n].first >> q[n].second;
else
m = n;
for (int i = 0; i < 100; i++) simulate_anneal(); // 进行100次模拟退火
cout << ans << endl;
return 0;
}
acwing2680. 均分数据
题意: 有n个正整数 a [ i ] a[i] a[i],要求将这n个数字分为m组,使得每组数据的平均数值和最平均,打印划分完后最小的均方差。 M < = N < = 20 , 2 < = M < = 6 , 1 < = A i < = 50 M <= N <= 20, 2 <= M <= 6, 1 <= A_i <= 50 M<=N<=20,2<=M<=6,1<=Ai<=50
题解: 模拟退火+贪心。每次随机将所有数字排序,然后随机选择2个数字,进行交换,计算交换后的值哪个更优来决定是否交换。calc时使用贪心的思路:维护m组数字的和,然后每次将当前数字放入最小的那组数字内,从而计算出最小值。
代码:
#include <bits/stdc++.h>
using namespace std;
const int N = 25, M = 10;
int n, m;
int w[N], s[M];
double ans = 1e8;
double calc() { // 贪心计算当前每组的均方差
memset(s, 0, sizeof s);
for (int i = 0; i < n; i++) { // 每次把当前值放入分组的最小值内
int k = 0;
for (int j = 0; j < m; j++)
if (s[j] < s[k]) k = j;
s[k] += w[i];
}
double avg = 0;
for (int i = 0; i < m; i++) avg += (double)s[i] / m; // 计算每组的均值
double res = 0;
for (int i = 0; i < m; i++) res += (s[i] - avg) * (s[i] - avg); // 计算每组的均方差
res = sqrt(res / m);
ans = min(ans, res);
return res;
}
void simulate_anneal() {
random_shuffle(w, w + n); // 先随机初始序列
for (double t = 1e6; t > 1e-6; t *= 0.95) {
int a = rand() % n, b = rand() % n; // 每次随机选择两个位置
double x = calc(); // 计算没有交换时的初始值
swap(w[a], w[b]); // 交换
double y = calc(); // 计算交换完后的值
double delta = y - x; // 计算差值
if (exp(-delta / t) < (double)rand() / RAND_MAX) swap(w[a], w[b]); // 如果不满足条件,不跳,那么就再交换一次,等价于不叫唤
}
}
int main() {
cin >> n >> m;
for (int i = 0; i < n; i++) cin >> w[i];
for (int i = 0; i < 100; i++) simulate_anneal();
printf("%.2lf\n", ans);
return 0;
}
ICPC 2018 Nanjing D. Country Meow
题意: 给出空间上的n个点,求出一个点的位置,使得到这个点距离最远的点的距离最小
题解: 模拟退火。本题和acwing3167. 星星还是树基本一样。
代码:
#include <bits/stdc++.h>
using namespace std;
typedef pair<double, double> PDD;
const int N = 110;
struct Point {
double x, y, z;
}q[N];
int n;
double ans = 1e8;
double rand(double l, double r) {
return (double)rand() / RAND_MAX * (r - l) + l;
}
double get_dist(Point a, Point b) {
double dx = a.x - b.x;
double dy = a.y - b.y;
double dz = a.z - b.z;
return sqrt(dx * dx + dy * dy + dz * dz);
}
double calc(Point p) { // 计算取当前点时,到所有的点的距离
double res = 0;
double tmp_max = 0;
for (int i = 0; i < n; i++) tmp_max = max(tmp_max, get_dist(p, q[i]));
ans = min(ans, tmp_max);
return tmp_max;
}
void simulate_anneal() { // 模拟退火
Point cur = {rand(-1e5, 1e5), rand(-1e5, 1e5), rand(-1e5, 1e5)}; // 先随机一个点
for (double t = 2e5; t > 1e-4; t *= 0.99) { // 按照降温系数进行降温
Point new_point= {rand(cur.x - t, cur.x + t), rand(cur.y - t, cur.y + t), rand(cur.z - t, cur.z + t)}; // 得到一个新的点
double dt = calc(new_point) - calc(cur); // 计算新的点和旧的点的差值
if (exp(-dt / t) > rand(0, 1)) cur = new_point; // 如果新的点更优,那么跳到新的点去
}
}
int main() {
scanf("%d", &n);
for (int i = 0; i < n; i++) scanf("%lf%lf%lf", &q[i].x, &q[i].y, &q[i].z);
while ((double)clock() / CLOCKS_PER_SEC < 0.8) simulate_anneal(); // 跑100次
printf("%.10lf\n", ans);
return 0;
}
NC17248 I、Expected Size of Random Convex Hull
大意:
给出一个三角形三个点的坐标,要求从中随机选择n个点,这n个点中在凸包上的点的个数的期望
思路:
由于拉伸三角形是不影响结果的,所以直接随机化n个点打表求即可
#include <algorithm>
#include <cmath>
#include <cstdio>
using namespace std;
#define maxn 15
const double eps = 1e-8;
const double PI = acos(-1.0);
struct Point {
double x, y;
double k;
Point() {}
Point(double _x, double _y) {
x = _x;
y = _y;
}
Point operator-(const Point &b) const { return Point(x - b.x, y - b.y); }
double operator^(const Point &b) const { return x * b.y - y * b.x; }
double operator*(const Point &b) const { return x * b.x + y * b.y; }
};
int sgn(double x) {
if (fabs(x) < eps) return 0;
if (x < 0)
return -1;
else
return 1;
}
double dist(Point a, Point b) {
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
Point list[maxn];
int Stack[maxn], top;
//相对于list[0]的极角排序
bool _cmp(Point p1, Point p2) {
double tmp = (p1 - list[0]) ^ (p2 - list[0]);
if (sgn(tmp) > 0)
return true;
else if (sgn(tmp) == 0 && sgn(dist(p1, list[0]) - dist(p2, list[0])) <= 0)
return true;
else
return false;
}
int Graham(int n) {
Point p0;
int k = 0;
p0 = list[0];
//找最下边的一个点
for (int i = 1; i < n; i++) {
if ((p0.y > list[i].y) || (p0.y == list[i].y && p0.x > list[i].x)) {
p0 = list[i];
k = i;
}
}
swap(list[k], list[0]);
sort(list + 1, list + n, _cmp);
if (n == 1) {
top = 1;
Stack[0] = 0;
return 1;
}
if (n == 2) {
top = 2;
Stack[0] = 0;
Stack[1] = 1;
return 2;
}
Stack[0] = 0;
Stack[1] = 1;
top = 2;
for (int i = 2; i < n; i++) {
while (top > 1 && sgn((list[Stack[top - 1]] - list[Stack[top - 2]]) ^
(list[i] - list[Stack[top - 2]])) <= 0)
top--;
Stack[top++] = i;
}
return top;
}
void random(int n) {
for (int i = 0; i < n; i++) {
double x = rand(), y = rand();
x /= RAND_MAX, y /= RAND_MAX;
if (x < y) swap(x, y);
list[i].x = x, list[i].y = y;
}
}
int main() {
/*
for(int n=3;n<=10;n++)
{
srand(time(0));
int cnt=1e7;
double p=cnt,ans=0;
while(cnt--)
{
random(n);
ans+=Graham(n);
}
printf("n=%d %.5f\n",n,ans/p);
}*/
int a;
for (int i = 1; i <= 6; i++) scanf("%d", &a);
int n;
scanf("%d", &n);
double ans = 3;
for (int i = 3; i < n; i++) ans += 2.0 / i;
printf("%.5f\n", ans);
return 0;
}
2018 ACM 上海大都会赛A Fruit Ninja
大意:
给出n个点。问是否存在至少有M个点在一条直线上,使得M/N>=x,X为只有一位小数的实数。n为1e4
思路:
如果枚举两个点,那么就是n^2的算法,会超时, 所以可以想到利用随机化,随机选取直线,看直线上是否存在大于等于m个点
#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 5;
typedef long long LL;
int t;
int n, flag;
double a[N], b[N], x;
const double eps = 1e-8;
const double pi = acos(-1.0);
// 和0做比较
int sgn(double x) {
if (fabs(x) < eps) return 0; // =0
if (x < 0)
return -1; // < 0
else
return 1; // > 0
}
struct Point {
double x, y;
Point() {}
Point(double _x, double _y) {
x = _x;
y = _y;
}
bool operator==(Point b) const {
return sgn(x - b.x) == 0 && sgn(y - b.y) == 0;
}
bool operator<(Point b) const {
return sgn(x - b.x) == 0 ? sgn(y - b.y) < 0 : x < b.x;
}
Point operator-(const Point &b) const { return Point(x - b.x, y - b.y); }
//叉积
double operator^(const Point &b) const { return x * b.y - y * b.x; }
//点积
double operator*(const Point &b) const { return x * b.x + y * b.y; }
Point operator+(const Point &b) const { return Point(x + b.x, y + b.y); }
Point operator*(const double &k) const { return Point(x * k, y * k); }
Point operator/(const double &k) const { return Point(x / k, y / k); }
};
struct Line {
Point s, e;
Line() {}
// 两点确定一条线段
Line(Point _s, Point _e) {
s = _s;
e = _e;
}
//返回点和直线关系
// 1 在左侧; 2 在右侧; 3 在直线上
int relation(Point p) {
int c = sgn((p - s) ^ (e - s));
if (c < 0)
return 1;
else if (c > 0)
return 2;
else
return 3;
}
};
void work() {
int u = rand() % n, v = rand() % n;
if (u == v) v = (v + 1) % n;
Line l = Line(Point(a[u], b[u]), Point(a[v], b[v]));
int num = 0;
for (int i = 0; i < n; i++) {
Point p = Point(a[i], b[i]);
if (l.relation(p) == 3) num++;
}
if (num >= n * x) flag = 1;
}
int main() {
cin >> t;
while (t--) {
flag = 0;
cin >> n >> x;
for (int i = 0; i < n; i++) {
cin >> a[i] >> b[i];
}
while ((double)clock() / CLOCKS_PER_SEC < 0.8) work();
if (flag) cout << "Yes" << endl;
else cout << "No" << endl;
}
return 0;
}