模拟退火算法

-
OI
随机化

简介 ​

模拟退火算法（Simulated Annealing）是一种基于概率论的随机化算法，是对爬山算法的改进，它通过模拟退火过程来寻找全局最优解。

一般 OI 比赛中可以使用模拟退火获得更高的部分分，在部分题目可能由于数据强度不够，甚至可能通过。

算法描述 ​

模拟退火算法随机求得一个新的解，并与当前解进行比较，如果新的解更优，则接受，否则以一定概率接受较差的解。

该算法需要三个参数：温度参数 T0T_0，降温参数 Δt\Delta t，终止温度 TkT_k。

我们记录当前的温度 T=T0T=T_0，在当前温度下进行一次 转移，得到新的温度 T=T×ΔtT=T\times \Delta t，当 T<TkT<T_k 时终止退火。

转移是更新答案的过程。具体而言，根据当前的温度，在当前状态的基础上随机得到一个新状态，求得两个状态的答案差值 ΔE(ΔE≥0)\Delta E(\Delta E\ge 0)。如果新状态更优，那么将新状态作为当前状态，否则以概率：

P(ΔE)=e−ΔETP(\Delta E)=e^{\frac{-\Delta E}{T}}

将新状态作为当前状态。

这样描述可能很抽象，下面用几个例题来进行分析。

例题 ​

P1337 [JSOI2004] 平衡点 / 吊打XXX ​

题目描述：

求给定 nn 个点的带权费马点。具体地，若给定点点权为 viv_i，到费马点的距离为 disidis_i，那么代价为 ∑vi×disi\sum v_i\times dis_i，找到一个点最小化代价。输出点的坐标。

n≤1000n\le 1000，∣xi∣,∣yi∣≤10000\lvert x_i\rvert,\lvert y_i\rvert\le 10000，0<wi≤10000<w_i\le 1000。

这题显然具有单调性，可以通过三分解决。但是同样可以使用模拟退火通过。

下面结合代码讲解如何使用模拟退火：

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)[l,r) 的随机数，是左闭右开的。

更详细的内容可见 随机数与随机函数 一文。

P2503 [HAOI2006] 均分数据 ​

将 nn 个正整数分成 mm 组，每组的权值为该组的数字和，希望最小化权值的方差。求最小方差的值。

m≤n≤20,2≤m≤6,ai∈[1,50]m\le n\le 20,2\le m\le 6,a_i\in[1,50]。

这题应该可以使用状压 DP 进行。考虑使用模拟退火，每一次随机得到一个数组，规定每一组在数组 aa 中都是连续的一段，使用简单 DP 可以在 O(n3)O(n^3) 的时间内计算出答案。每一次新状态由原状态随机交换位置得到。类似地进行模拟退火：

Sample Code (c++)
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 ​

题目描述：

给定平面上 NN 个圆和 MM 个点。请你选择一个点和不超过 RR 的半径作圆，在这个圆与给定圆不相交的情况下，求最多能够覆盖多少个点。

0≤N≤10,0<M≤1030\le N\le 10,0<M\le 10^{3}，1≤R,ri≤2×1041\le R,r_i\le 2\times 10^{4}，所有坐标绝对值不超过 2×1042\times 10^{4}。

Sample Code (C++)
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] 分金币 ​

题目描述：

给定 nn 个数字，将 nn 个数字分成两组，两组的数字数量相差不超过 11，求两组金币价值差的最小值。

1≤n≤301\le n\le 30，每一个数字 viv_i 满足 1≤vi≤2301\le v_i\le 2^{30}。

Sample Code (C++)
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);
  }
}
Kirchhoff 矩阵-树定理
行列式