模拟退火

模拟退火

1.算法描述

1.1 概述

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

它与爬山算法最大的不同是,在寻找到一个局部最优解时,赋予了它一个跳出去的概率,也就有更大的机会能找到全局最优解。

在 OI 领域,对应的,每次随机出一个新解,如果这个解更优,则接受它,否则以一个与温度和与最优解的差相关的概率接受它。

以下图片来自:https://oi-wiki.org/misc/simulated-annealing/

img

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越大越容易接受。针对具体问题为:

  1. 如果是求凹函数的最值,那么判断条件写为:
if (exp(-dt / t) > rand(0, 1)) cur = new;  // 如果概率大于(0 ~ 1),那么就跳到新的点
else continue;  // 否则不跳
  1. 如果是求凸函数的最值,那么判断条件写为:
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 1n100,0xi,yi10000

题解: 模拟退火。每次先随机化一个点,然后按照降温系数降低温度,每次跳转到一个新的点,如果新的点更优,那么就跳到这个点去。

代码:

#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;
}
  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值