模拟退火

模拟退火

解决问题类型

对于最小化问题：

\[\min_{x \in \mathcal X} f(x) \]

其中，\(\mathcal X\) 表示可能的解空间

相应的

最大化问题：

\[\max_{x\in \mathcal X} f(x) \]

可以转换为最小化问题：

\[\begin{equation} \max_{x\in \mathcal X} f(x) = - \min_{x\in\mathcal X} -f(x) \end{equation} \]

实在无法理解这些符号，只用记住一句话：

SA(模拟退火) 解决的是一类求最值得问题

解决方案

模拟退火算法（Simulate Anneal，SA）是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。

模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的。V.Černý在1985年也独立发明此演算法。模拟退火算法是解决TSP问题的有效方法之一。

模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。

模拟退火算法是一种通用的优化算法，其物理退火过程由加温过程、等温过程、冷却过程这三部分组成。

——百度百科

明确一点，退火来自：

模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。

其实就是从生活中来的

方案大家都以明了，所以吾就阐明其正确性了

为了方便表达

记 \(\mathcal X^*\) 为最优解的集合

\(\text{Pr}\{·\}\) 表示一个事件的概率

有：

定理

对于一个温度序列 \(T_k\)，有 \(\lim_{k \to \infin} T_k = 0\)，则

\[\lim_{k\to \infin} \text{Pr}\{X_k \in \mathcal X^*\} = 1 \]

什么，看不懂，吾来带着汝等详尽分析：

首先，吾等应先明确此算法的一个特性：

当温度 \(T\) 越高，则算法搜索的空间（或者说范围）也就越大

温度越小，算法逐渐趋近于正解，最终在正解处停下：

这张图就很好的体现了此点

而定理说明了一件事：

当温度趋近于 \(0\) （一个极小的，大于 \(0\) 的数，例如1e-15）时，找到的解 \(X\) 是最优解的概率来到了 \(1\) ，没座儿，就是 \(1\)

对应的，还有证明，不过有 2635 个字，且晦涩难懂

好，既然以明晰原理（正确性保障）

吾等来详细分析做法

首先，要有一个初始温度 \(T_0\) ，一个降温系数 \(\alpha\) ，必要时可以在一次降温内做多次（iterations） 操作，要有一个终止温度 \(\text{EPS}\)

框架就是：

while (T >= EPS) {
    
    for (int iter = 1; iter <= iterations; iter++) {
        // ... 操作
    }
    
    T *= coolingRate; // 降温
}

当然，我们需要一个初始解 \(x_0\) ，一般可以用错误的贪心或者….（自己想）

像这样：

int x = init();

while (T >= EPS) {} // ......

好，来看转移

我们用随机数 rand() 来求出一个新的解 \(x^*\) ，如果更加优秀（\(f(x^*) < f(x_\text{best})\)），算法会无条件下接受新解

但是如果新解更劣（\(f(x^*)>f(x_\text{best})\)），注意，不同于普通贪心的点来了：

有一定概率接受这个更劣解

而这个概率由Metropolis准则给出：

\[\text{Pr}\{x_\text{best} \to x^*\} = \exp(-\frac{f(x^*) - f(x_\text{best})}{T}) \]

这说明了什么？

说明了随着温度的下降，接受更劣解的概率会越来越小，也就是越来越稳定

为什么敢这样子做？：

• 首先，这源自自然规律（物理中的一个模型还是定理，“我忘了，QWQ）
• 其次，正是因为随着温度降低，找到正确解的概率变大才这么做

而这正保证了SA的时间复杂度

所以正确性完全不用担心，doge

可能就是调参有点痛苦……

for (int iter = 1; iter <= iterations; iter++) {
    int/*或者 double*/ y = ...; // 从合理的范围内取值，称作 x 的邻域 N(x)。通常是均匀取值，即取到 y 的概率是 1 / |N(x)|
    
    int delta = f(y) - f(x);
    
    if (delta < 0 || exp(-delta / T) > (rand() / (double)RAND_MAX)) { // 更优解，或满足概率，更新解
        // ...
    }
}

然后就完了

整的模版：

int x = init();

while (T >= EPS) {
    
    for (int iter = 1; iter <= iterations; iter++) {
    int/*或者 double*/ y = ...; // 从合理的范围内取值，称作 x 的邻域 N(x)。通常是均匀取值，即取到 y 的概率是 1 / |N(x)|
    
    int delta = f(y) - f(x);
    
    if (delta < 0 || exp(-delta / T) > (rand() / (double)RAND_MAX)) { // 更优解，或满足概率，更新解
        // ...
    }
}
    
    T *= coolingRate; // 降温
}

*复杂度分析 （拓展内容）

壹.

由 ？？？ 给出：

会有：

\[\Omicron(\log_{\frac{1}{\alpha}} \frac{T}{T_\text{final}}) \]

次迭代

或者：

\[\Omicron(\log_{\alpha}\frac{T_\text{final}}{T}) \]

贰. 下界

复杂度下界由复杂度下界定理给出（很奇怪的名字）

前置：

能量障碍

从 \(x\to y\) 的能量障碍定义为：

\[\Gamma(x, y) = \min_{\text{path} \ x \to y} \max_{z\in \text{path}}f(z) \]

状态 \(x\) 的深度

\[d(x) = \min_{y \ \text{s.t.} \ f(y) < f(x)} [\Gamma(x, y) - f(x)] \]

若 \(x \in \mathcal X^*\) ，则 \(d(x) = \infin\)

全局最大深度

\[D = \max_{x\not \in \mathcal X^*}d(x) \]

SA复杂度下界定理

对于任意的冷却方案，要达到收敛条件，至少经过

\[k = \Omega (\exp(\frac{D}{\max_k T_k})) \]

次迭代

这说明模拟退火是一个亚指数时间算法，但在处理 NP-hard 类问题上仍然很优秀

*注：一般问题中的 \(D\) 并不大

---

例题

P1337 吊打XXX

最小化：

\[\sum_{i = 1}^n w_i\sqrt{(x_i - X)^2 + (y_i - Y)^2} \]

每次的邻域 \(N(x)\) 的选择：

随机游走:

lf X = result.x + ((lf)(rand() % 1000 - 500) / 10000.0);
lf Y = result.y + ((lf)(rand() % 1000 - 500) / 10000.0);

View Code

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

using lf = double;
using ll = long long;
using ull = unsigned long long;

lf centerX, centerY;
lf minX = INT_MAX, maxX = INT_MIN, minY = INT_MAX, maxY = INT_MIN;

const int maxn = 1e3 + 5;

int x[maxn], y[maxn];
int weight[maxn];

int n;

lf f(lf X, lf Y) {
	lf result = 0.0;
	for (int i = 1; i <= n; i++) result += weight[i] * sqrt((lf)(x[i] - X) * (x[i] - X) + (lf)(y[i] - Y) * (y[i] - Y));
	return result;
}

struct Point {
	lf x, y;
	Point() = default;
	Point(lf _x_, lf _y_): x(_x_), y(_y_) {}
};

bool outOfRange(lf x, lf y) {
	return (x < minX || x > maxX) || (y < minY || y > maxY);
}

const lf EPS = 1e-15;

Point SA(Point p) {
	Point result = {p.x, p.y};
	lf best = f(result.x, result.y);

	lf T = 13.0;
	lf coolingRate = 0.9996;

	srand(time(nullptr));

	while (T >= EPS) {

		lf X = result.x + ((lf)(rand() % 1000 - 500) / 10000.0);
		lf Y = result.y + ((lf)(rand() % 1000 - 500) / 10000.0);

		if (outOfRange(X, Y)) {
			T *= coolingRate;
			continue;
		}

		lf current = f(X, Y);

		lf delta = current - f(result.x, result.y);

		if (delta < 0 || exp(-delta / T) > (lf)(rand() / RAND_MAX)) {
			result.x = X;
			result.y = Y;
			if (delta < 0) best = current;
		}

		T *= coolingRate;
	}

	return result;
}

int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	cout.tie(nullptr);

	cin >> n;

	lf sumX = 0.0, sumY = 0.0;
	
	for (int i = 1; i <= n; i++) {
		cin >> x[i] >> y[i] >> weight[i];
		minX = min(minX, (lf)x[i]);
		maxX = max(maxX, (lf)x[i]);
		minY = min(minY, (lf)y[i]);
		maxY = max(maxY, (lf)y[i]);
		sumX += x[i];
		sumY += y[i];
	}
	
	centerX = sumX / n;
	centerY = sumY / n;

	Point result = SA((Point){
		centerX, centerY
	});
	for (int i = 1; i <= 2; i++) result = SA(result);
	
	cout << fixed << setprecision(3) << result.x << " " << result.y << "\n";

	return 0;
}

当然，调参是个很严重的问题！！！

交了3页，哭

（如图）