隨機化演算法屬於省選芝士體系

0x01 前置芝士

你只需要會 rand 就可以啦！

當然如果你想理解的更透徹也可以先看看 爬山演算法

0x02 關於退火

退火是一種金屬熱處理工藝，指的是將金屬緩慢加熱到一定溫度，保持足夠時間，然後以適宜速度冷卻。目的是降低硬度，改善切削加工性；消除殘餘應力，穩定尺寸，減少變形與裂紋傾向；細化晶粒，調整組織，消除組織缺陷。準確的說，退火是一種對材料的熱處理工藝，包括金屬材料、非金屬材料。而且新材料的退火目的也與傳統金屬退火存在異同。

也就是你可以理解為，金屬在逐漸降溫的過程中，它含的某種物質會趨於恆定。

換成 OI 的話說，如果每個溫度對應一個答案，則隨著溫度越來越小，答案波動的範圍就會越來越小，最後答案趨於恆定。例如下面這個經典圖片，便是利用這個思想，最後得到了一個多峰函式的最值。模擬退火就是用於求解多峰函式的最值問題。

不過是不是看起來很玄學？現在我們來細化一下模擬退火思想。

0x03 模擬退火

你如果仔細看，剛剛那個圖中答案好像和溫度沒有直接聯絡，也就是說答案是隨機的。這個說法不太完全，答案和溫度確實沒有太大聯絡，但溫度卻可以限定答案所在的範圍。

首先我們引入幾個引數：當前最優解 \(A_0\)，新的解 \(A\)，上一個被接納的解 \(A_1\)，解的“變動量” \(ΔA\)（即 \(A\) 與 \(A_0\) 的差值，這裡規定為 \(A - A_0\)），開始的溫度 \(T_0\)，當前溫度 \(T\)，最終的溫度 \(T_{esp}\)。溫度變動量 \(q\)。接納：指對於一個解，我們認為它有可能是正確答案，或者說與正確答案有聯絡。

以下均以求多峰函式最大值為例來模擬整個過程。

首先在一開始，我們找到一個你認為接近正確答案的解 \(A_1\)，將 \(A_0\) 初始置為 \(A_1\)。然後你需要在溫度限制的範圍內，去合理的再取到一個解，即 \(A\)。至於這個溫度限制的範圍其實很簡單，你只需要求到一個滿足條件的隨機數即可。例如：你有一個橫座標，如果你要再求到一個橫座標，你就可以寫成。

double cx = now.x + ((rand() << 1) - RAND_MAX) * t; 
// RAND_MAX 是 <stdlib.h> 中的一個巨集，也就是隨機數的最大值。
// 通過 (rand << 1) - RAND_MAX 我們就得到了一個可能為正可能為負的隨機數
 

這樣的話，你就會發現隨著 \(T\) 的不斷減小，新的解與上一個解的差距就越小，即達到了我們的目的。

現在我們有了 \(A\)，自然就可以求到所謂 \(ΔA\)。

得到 \(ΔA\)後：

• \(ΔA > 0\) ：也就是說當前解大於當前最優解，因為我們是求最大值，所以此時顯然更新 \(A_0\)，並更新 \(A_1 = A\)。

• \(ΔA \leq 0\) ：這就比較麻煩了，首先這樣的情況是肯定不能更新 \(A_0\) 的。但我們考慮一下 \(A_1\)，如果接納 \(A\)，那麼我們下次就是從 \(A\) 開始擴充套件新的點，這很明顯有可能不是最優的（你從多峰函式的一個峰的半山腰跌到了一個山谷）。但也有可能是更優的，因為你現在接納的這個點，它可能不在最終答案的那一座峰上，所以你需要跳往另一座峰，也就是你需要接納 \(A\)。這下該怎麼辦呢？有一個非常玄學的東西，叫：Metropolis接受準則。它告訴我們有一定概率去接納 \(A\)，而這個概率只需比較 \(\exp(-delta / t) * RAND\_MAX\) 和 \(rand()\) 的大小即可。

在做完這些操作後，我們讓 \(T = T \times q\)。

好了，我們現在又有了一個 \(A_1\)，這個 \(A_1\) 可能等於 \(A_0\) 也可能等於 \(A\)。我們將這個 \(A_1\) 再去執行上述操作，反覆執行，隨著溫度的減小，我們就可以求到一個趨於恆定的答案了！

當然還有幾個未處理的地方。

首先：

1. \(T_0\) 的取值：根據題目而定，我通常使用 \(2000\) 到 \(5000\) 的一個整數。
2. \(q\) 的取值：在前人的總結中，\(q\) 近似於 \(0.996\)。
3. \(T_{esp}\) 的取值：這也需要視題目而定，如果這個值越小，按理說答案精度就越大，所以如果你不怕超時就可以開的小一點。

最後，關於時間複雜度，因為此演算法使用了大量隨機數，所以其時間複雜度近似於 \(O(rp)\)。

0x04 例題

LINK

一句話題意：給定 \(n\) 個座標，求這 \(n\) 個座標的 費馬點。

定義 \(sum\) 表示一個座標到 \(n\) 個點的距離和。

首先，你會發現答案是一個函式，因為一個座標只對應一個 \(sum\)。於是我們其實就是要求這個多峰函式 \(sum\) 的最小值。

那麼就是模擬退火嘛。讀者可以在根據下面這個程式碼理解一下實現。

值得注意的是，這道題函式 \(sum\) 的“橫座標”是一個座標。

#include <bits/stdc++.h>
using namespace std;

int read() {
    int k = 1, x = 0;
    char s = getchar();
    while (s < '0' || s > '9') {
        if (s == '-')
            k = -1;
        s = getchar();
    }
    while (s >= '0' && s <= '9') {
        x = (x << 3) + (x << 1) + s - '0';
        s = getchar();
    }
    return x * k;
}

const int MAXN = 5e3 + 5;
const double q = 0.996;
// 溫度變動量
struct node {
    double x, y;
    node() {}
    node(double X, double Y) {
        x = X;
        y = Y;
    }
} a[MAXN], now, anst;
// now 初值為 A1，在這裡是 (0, 0)
double ans = 1e18;
// 答案初值

int n;
double f(double x, double y) {
    double ret = 0;
    for (int i = 1; i <= n; i++) ret += sqrt((x - a[i].x) * (x - a[i].x) + (y - a[i].y) * (y - a[i].y));
    return ret;
}

void Make_Fire() {
    double t = 3000;
    while (t > 1e-16) {
        double cx = now.x + ((rand() << 1) - RAND_MAX) * t;
        double cy = now.y + ((rand() << 1) - RAND_MAX) * t;
        // 求到 A 的橫座標，在這裡是 (cx, cy)
        double cnt = f(cx, cy);
        // 求到 A
        double delta = cnt - ans;
        // 求到差
        if (delta < 0) { // 如果差是小於 0 的
            now = node(cx, cy);
            // 接納它
            anst = node(cx, cy);
            // 根據題目需要，更新答案的座標
            ans = cnt;
            // 跟新答案
        } else if (exp(-delta / t) * RAND_MAX > rand())
        // 玄學接受準則
            now = node(cx, cy);
            // 接納它
        t *= q;
        // 降溫
    }
}

int main() {
    srand(998244353);
    n = read();
    ans = 1e18;
    for (int i = 1; i <= n; i++) 
        scanf("%lf %lf", &a[i].x, &a[i].y);
    for (int i = 1; i <= 5; i++) Make_Fire();
    printf("%.2lf %.2lf\n", anst.x, anst.y);
    return 0;
}