(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)， 每个T值的迭代次数L
(2) 对k=1，……，L做第(3)至第6步：
(3) 产生新解S′
(4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数
(5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/T)接受S′作为新的当前解.
(6) 如果满足终止条件则输出当前解作为最优解，结束程序。
终止条件通常取为连续若干个新解都没有被接受时终止算法。
(7) T逐渐减少，且T->0，然后转第2步。
1. 模拟退火原理原理
模拟退火：是一种随机算法，用于解决最优化问题。要求求解的问题对应的函数要有连续性。
模拟退火算法是模拟物理过程，有如下参数：
（1）温度t：即步长。分为初始温度和终止温度，对应代码中就是初始搜索范围和终止搜索的范围。
（2）衰减系数：每次搜索范围减小的比例，是(0, 1)中的一个数，可以取0.999，需要手动调节。
在每次迭代的过程中，我们在给定步长区间内随机一个新点，令dt = f(新点)-f(当前点)，如果求函数极小值的话，分为两种情况：
（1）dt<0，则跳到新点；
（2）dt>0，则以一定该概率跳到该点，且dt越大，跳过去的概率越低。
跳过去的概率值可以取为 e − d t / t e^{-dt/t} e−dt/t。
模拟退火的过程可能会收敛到局部最优解，但是这个过程我们可以做多次，这样收敛到局部最优解的概率就很小了。比如达到局部最优解的概率是0.99，则我们做1000次，达到局部最优解的概率是： 0.9 9 1000 ≈ 4.3 × 1 0 − 5 0.99^{1000} \approx 4.3 \times 10^{-5} 0.991000≈4.3×10−5。
2. AcWing上的模拟退火题目
AcWing 3167. 星星还是树
问题描述
问题链接：AcWing 3167. 星星还是树
分析
本题求解的这个点是费马点，即到所有点距离和最小的点。如果是一维的，排个序找中位数即可。
可以证明，这个函数是个凸函数，具有连续性。使用模拟退火求解即可。
代码
C++
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <ctime>
#define x first
#define y second
using namespace std
typedef pair<double, double>PDD
const int N = 110
int n
PDD q[N]
double ans = 1e8
// 返回[l, r]之间的随机小数
double rand(double l, double r) {
return (double)rand() / RAND_MAX * (r - l) + l
}
double get_dist(PDD a, PDD b) {
double dx = a.x - b.x, dy = a.y - b.y
return sqrt(dx * dx + dy * dy)
}
// 计算p到给定点的距离和
double calc(PDD p) {
double res = 0
for (int i = 0i <ni++)
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 = 1e4t >1e-4t *= 0.9) {
PDD np(rand(cur.x - t, cur.x + t), rand(cur.y - t, cur.y + t))
double dt = calc(np) - calc(cur)
if (exp(-dt / t) >rand(0, 1)) cur = np
}
}
int main() {
scanf("%d", &n)
for (int i = 0i <ni++) scanf("%lf%lf", &q[i].x, &q[i].y)
for (int i = 0i <100i++) simulate_anneal()
printf("%.0lf\n", ans)
return 0
}
登录后复制
AcWing 2424. 保龄球
问题描述
问题链接：AcWing 2424. 保龄球
分析
本题需要求解最大值，相当于求全排列中的最大值。
每次我们可以随机交换两个轮次，计算交换前后的差距，更新答案。
代码
C++
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <ctime>
#define x first
#define y second
using namespace std
typedef pair<int, int>PII
const int N = 55
int n, m // n: 规定的轮次 m: 实际的轮次
PII q[N]
int ans
int calc() {
int res = 0
for (int i = 0i <mi++) {
res += q[i].x + q[i].y
if (i <n) {
if (q[i].x == 10) res += q[i + 1].x + q[i + 1].y
else if (q[i].x + q[i].y == 10)
res += q[i + 1].x
}
}
ans = max(ans, res)
return res
}
void simulate_anneal() {
for (double t = 1e4t >1e-4t *= 0.99) {
auto a = rand() % m, b = rand() % m
int x = calc()
swap(q[a], q[b])
// 交换后进行的轮次 n + (q[n - 1].x == 10) 等于 实际轮次m
if (n + (q[n - 1].x == 10) == m) {
int y = calc()
int dt = y - x
// 如果dt>0, 则不用恢复原状
if (exp(dt / t) <(double)rand() / RAND_MAX)
swap(q[a], q[b])
} else swap(q[a], q[b])
}
}
int main() {
cin >>n
for (int i = 0i <ni++) cin >>q[i].x >>q[i].y
if (q[n - 1].x == 10) m = n + 1, cin >>q[n].x >>q[n].y
else m = n
// for (int i = 0i <100i++) simulate_anneal()
// 卡时写法: 卡0.1秒
while ((double)clock() / CLOCKS_PER_SEC <0.1) simulate_anneal()
cout <<ans <<endl
return 0
}
登录后复制
AcWing 2680. 均分数据
问题描述
问题链接：AcWing 2680. 均分数据
分析
这里可以随机将这些数据放置到某个组中，为了使得收敛的速度更快，可以采用贪心的策略将n个数据放置到m个组中。
每次找到和最小的组，将该数据放到该组中。
代码
C++
5.3.1 仿真退火与系统能量最小状态
退火是众所周知的一种金属加工工艺，指金属从高温逐渐冷却的一种热处理方式。高温时，分子相对运动的热活动性大，在冷却时热活动性逐渐消失，最后形成原子的有序排列，处于热动力系统的最小能量状态。快速的突然冷却常不能达到最小能量状态，只有在缓慢冷却的条件下，在原子失去热活动性的同时有足够的时间让它们重新分布，才能达到最小能量状态。
退火这一物理过程的计算机仿真就成为一种求泛函极小的最优化算法。与退火过程类似，这种算法要求温度缓慢下降，但这与降低成本相矛盾。
仿真退火算法将非线性优化问题与统计力学中的热平衡问题进行类比，找出了优化问题求解与固体退火过程之间的相似性，从而另辟优化问题数值求解的新途径。
对于非线性优化问题，设非负目标函数为φ（x），要求
φ（v*）=minφ（vi），vi∈V （5.3.1）
其中vi为某一个可能的解，所有vi的集合记为V=｛vi｝，i=1，2，…，N；v*为V中使φ取极小的一个解估计。在与退火过程类比时，φ视为系统的内能，要使它由高到低，最后取极小而终止。vi为解空间的一个点，对应为固体内的一种粒子状态，它的扰动范围随内能的减小而逐渐减小。温度T用控制参数T′=kT类比，其中k为波尔兹曼系数。仿真退火算法在求解优化问题时先取很大的T′，T′可称为仿真温度，然后逐渐减小T′，计算时不同的迭代步对应着时间与温度T′的变化。当系统温度高内能大时粒子的随机扰动大，相应于系统状态的vi扰动的自由度也大。随迭代进行温度降低，当前解vi不断更新而形成序列，构成了对解空间的搜索。
按照统计力学中关于退火过程的理论，仿真退火算法遵循的基本规律是：某种状态vi在当前出现的概率服从波尔兹曼分布
地球物理数据处理教程
由（5.3.2）式可见，当仿真温度T′=kT→0和目标函数φ取极小时（vi=v*），系统处于内能最低的基态，这时概率f*取极大而接近于1，说明系统最稳定，熵最小。因此，仿真温度降至某一阀值以下，而且系统的熵达到最小可作为停止迭代的准则。
仿真退火算法的关键在于：在给定仿真温度T′i和当前状态vi之下，如何确定下一个状态vi+1，即如何构造迭代修改序列以便使空间中的搜索合理地进行，而最后找到使φ取全局极小的v*。仿真退火算法的创始人N.Metropolis于1953年提出了一种抽样算法。他说，假定当前的状态为v，利用随机扰动产生一种新扰动v′，其对应系统内能的变化为Δφ=φ（v′）-φ（v）。根据（5.3.2）式，这种新状态可以被接受的概率为
p=exp（-Δφ/T′） （5.3.3）
如果φ（v′）比φ（v）减小或相同，即Δφ≤0，则有概率p≥1，这种新状态毫无疑问要被接受。如果Δφ＞0，由（5.3.3）式计算p＜1，但也不一定为零。这时新状态虽然使目标函数增大，但为保证解空间的搜索是全局性的，不至于落入局部极小而无机会逃出，给予p＞0的机会是正确的，它给予从局部极小逃出的机会。根据（5.3.3）式建立了Metropolis内循环抽样算法，将在下面详述。
根据仿真退火原理，可以把算法的实现步骤大致归结如下。
第一步，选取i=0时的初始模型及仿真温度。随机选取初始地球模型v0，并依照经验选取初始温度T′0。也可通过实验选取T′0，例如，随机抽样一组模型 ｛vi｝，计算对应的目标函数 ｛φj｝，取 ｛φj｝ 的方差或最大Δφj的若干倍为T′0。
第二步，在给定vi和仿真温度T′i之下用Metropolis抽样对解空间进行搜索，具体由以下内循环迭代过程构成：
（1）令k=0时当前解估计v（0）=vi，参数T=T′i。
（2）按某一半随机准则在当前解v（k）附近产生一个邻近子集Vn（v（k））⊂V，它不能为空集，与v（k）的距离（或半径）与参数T成比例。在子集中随机地取出一个新的扰动模型v′作为下一步迭代的候选解，计算
Δφ=φ（v′）-φ（v（k）） （5.3.4）
（3）如Δφ≤0，则接收v′为下步迭代的当前解，即令v（k+1）=v′。如若Δφ＞0，则按（5.3.3）式计算概率p′，并计算一个在（0，1）之间的随机数p*。如果p′＞p*，则接受v′为下一步迭代的当前解，否则令v（k+1）=v（k），即保留原来的当前解估计。
（4）令k=k+1，如果k+1不超过某个规定的次数，返回（5.3.2）继续抽样搜索，否则用当前解估计v（k+1）跳出内循环，转入下述主程序的第三步。
第三步，按一定准则降温。此时令 i=i+1，Ti+1=λiTi，λi＜1。例如，λi可在（0.2，0.99）之间变化，但随i增大而增大。
第四步，判断是否停止外循环迭代，如果不是，则返回第二步，以当前vi+1=v（k+1）及Ti+1重复Metropolis抽样搜索。外循环停止计算的准则是T′小于某一规定阀值和系统的熵达到极小，或者由（5.3.3）式计算当前解的概率f→1等。
仿真退火算法的收敛性可以用马尔可夫链来研究，因为非平稳可数马尔可夫链可恰当地描述解空间的搜索序列。总的来说，在仿真退火算法中，如果仿真温度T′的变化等参数满足一定的条件收敛是有保证的，但收敛速度较慢。如果T′0足够高而且下降充分慢，以及对每一次下降，Metropolis抽样都很稳定，则能达到φ（vi）的全局极小值。这些条件正好要求加大计算成本，这是仿真退火固有的缺点。
5.3.2 仿真退火地震反演
进行地震反演时，目标泛函可选类似式（5.3.1）的任何形式。例如，对不考虑正则化的情况，以拟解选
地球物理数据处理教程
为计算第k次迭代时的目标泛函值，其中r为地震道反射系数序列，s为地震道数据，sk为第k次迭代时用 k-1次迭代结果计算合成的地震道。
初始温度T0可选为10000～20000。每一次迭代，温度下降为Tk+1=0.92Tk，这与混沌反演λk下降一致，几何级数递减。每一次迭代随机试验拟解的修改最大次数为3000。设当前第j个模型参数为rj，随机扰动的公式为
r′j=rj+Rand（I）·βTk （5.3.6）
式中Rand（I）为一个［-0.5，0.5 ］之间取值的随机数，β为与当前迭代温度有关的权。温度Tk越低时，权β应越小才能提高成功率。
按式（5.3.6）修改模型参数后作地震合成并代入（5.3.5）式求φk（r′），如果φk减小则修改成功，否则按Metropolis算法确定是否用r′代替r。接下去再作扰动试验。
对每个温度级（即每个k），规定最大扰动成功次数（如500），在此级成功修改多次后即降低温度进入（k+1）次迭代。
迭代停止的准则是：
（1）在k个温度级修改r3000次，成功率为零；
（2）温度已降到室温；
（3）|φk+1-φk|＜ε，ε为规定的小数。
按以上参数准则，利用仿真退火进行了与混沌反演相似的地震道反演，无误差结果示意如图5.7，子波为11点零相位Ricker子波，频率为30 Hz。
由图可见，仿真退火反演的分辨率是不错的，如中段低阻抗薄层有清晰显示；但是，方差较大，即抖动得厉害。
从目标泛函的下降曲线看，它呈现阶梯形状，与混沌反演的阶段性类似。在迭代230次以后φk已接近全局最小，最终在k=265时达到全局极小“山谷”，φmin=4.92×10-5。
以上是用双精度计算的反演结果。当用单精度计算时，迭代次数和φ曲线下降台阶减少，迭代次数在238次停止，φmin=4.88×10-4。
图5.7 仿真退火地震反演试验（精确数据）
（A）地层波阻抗模型，虚线为初始模型；（B）反演最终取得的波阻抗序列；（C）目标泛函随迭代的下降曲线
如果在数据中出现微小高斯噪音，则迭代226次，φmin=1.59×10-3，没有明显的全局极小“山谷”出现，此时计算时仍比混沌反演多一倍。
总之，仿真退火反演是一种非常稳定的非线性反演方法，但与混沌反演相比，优点并不突出。
