模拟退火算法
基本思路
爬山算法（Hill Climbing）：兔子朝着比现在高的地方跳去。它找到了不远处的最高山峰。但是这座山不一定是珠穆朗玛峰。这就是爬山算法，它不能保证局部最优值就是全局最优值。
模拟退火算法（Simulated Annealing）：兔子喝醉了。它随机地跳了很长时间。这期间，它可能走向高处，也可能踏入平地。但是，它渐渐清醒了并朝最高方向跳去。这就是模拟退火。
玄学算法。
模拟退火的主要步骤有几个：
- 设置初始温度 \(T\)，初始符合条件的答案；
- 通过某种神奇的方式，找到另一个符合条件的新状态；
- 分别将两个状态的答案计算出来，并作差得到 \(ΔE\)；
- 根据题目要求，贪心的决定是否更换答案（即 选择最优解）；
- 如果无法替换答案，则根据一定概率替换答案，即运用到平衡概率 \(\displaystyle e^{\frac{ΔE}{T}}\) 随机的决定是否替换；
- 每一次操作后，进行降温操作，即：将温度 \(T\) 乘上某一个系数，一般是 0.985−0.999 随具体题目（
随缘）定。（源：模拟退火 by peng-ym）
const double eps = 1e-15;
const double t0 = 0.985~0.999;
inline void mnth() { // (需要玄学调试)
double T = 初始温度;
while (T>eps) {
double now = 新状态 如 ans + (rand()*2 - RAND_MAX) * T;
double delta = f(now) - f(ans);
if (delta 更优) ans = now;
else if (exp(±delta / T) * RAND_MAX > rand()) ans = now; // 从局部最优解中跳出
T *= t0; // 降温
}
}
ans = 初始状态;
mnth();
print(ans); // 退火后的近似最优解
[JSOI2004] 平衡点
有 \(n\) 个重物，每个重物系在一条足够长的绳子上。每条绳子自上而下穿过桌面上的洞，然后系在一起。记 \(X\) 处为公共的绳结。假设绳子是完全弹性的（不会造成能量损失），桌子足够高（因而重物不会垂到地上），且忽略所有的摩擦。问绳结 \(X\) 最终平衡于何处。注意：桌面上的洞都比绳结 \(X\) 小得多，所以即使某个重物特别重，绳结 \(X\) 也不可能穿过桌面上的洞掉下来，最多是卡在某个洞口处。
如题，经过物理学分析，易得本题求解使 \(\sum m_is_{xi}\) 取最小值的 \(x\) 点位置。然后开始编写程序，然后面向数据编程。（玄学调参）
/* [JSOI2004] 平衡点
* Au: GG
*/
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <ctime>
using namespace std;
#define db long double // double 精度不够，改成 long double
const db eps=1e-15;
const db t0=0.98; // 0.99 会 TLE，改成 0.98
int n, px[1003], py[1003], w[1003];
db ansx, ansy;
inline db dis(db ax, db ay, db bx, db by) {
return sqrt((ax-bx)*(ax-bx) + (ay-by)*(ay-by));
}
inline db calc(db x, db y) {
db res=0.0;
for (int i=1; i<=n; ++i) res+=dis(x, y, px[i], py[i])*w[i];
return res;
}
inline void mnth() {
db T=200.0; // 温度调得挺暖和
while (T>eps) {
db nowx=ansx + (rand()*2-RAND_MAX) * T;
db nowy=ansy + (rand()*2-RAND_MAX) * T;
db delta = calc(nowx, nowy) - calc(ansx, ansy);
if (delta<0) ansx=nowx, ansy=nowy;
else if (exp(-delta/T) * RAND_MAX > rand()) ansx=nowx, ansy=nowy;
T*=t0;
}
}
int main() {
srand(time(0)); // 设定固定种子容易被卡（其实素数 19491001 比较少有人卡吧……）
scanf("%d", &n);
for (int i=1; i<=n; ++i)
scanf("%d%d%d", &px[i], &py[i], &w[i]), ansx+=px[i], ansy+=py[i];
ansx/=n; ansy/=n;
mnth();
printf("%.3Lf %.3Lf\n", ansx, ansy);
return 0;
}
[HAOI2006] 均分数据
已知 \(N\) 个正整数：\(A_1,A_2,\ldots,A_n\)。今要将它们分成 \(M\) 组，使得各组数据的数值和最平均，即各组的数据的数值和均方差最小。
均方差 $\displaystyle \sigma=\sqrt{\frac{\sum_{i=1}n(x_i-\overline{x})2}{n}} $.
算术平均值 $\displaystyle \overline{x}=\frac{\sum_{i=1}^nx_i}{n} $.
我们可以贪心地把 \(N\) 个数放入 \(M\) 组中，即每次放入数值和最小的组。但是这样不一定是最优解。
模拟退火算法。贪心的顺序很重要，于是不断修改贪心顺序，直到找到近似最优解。
/* [HAOI2006] 均分数据
* Au: GG
*/
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <ctime>
using namespace std;
const double eps=1e-15;
const double t0=0.998;
int n, m, p[23], sum[9];
double ave, ans, now;
inline void calc() {
now=0.0, memset(sum, 0, sizeof sum);
for (int i=1; i<=n; ++i) {
sum[1]+=p[i];
sort(sum+1, sum+m+1);
}
for (int i=1; i<=m; ++i) now+=(sum[i]-ave)*(sum[i]-ave);
now=sqrt(now/m);
}
inline void mnth() {
double T=10000.0; // 高炉炼铁
while (T>eps) {
int x=rand()%n+1, y=rand()%n+1; while (y==x) y=rand()%n+1;
swap(p[x], p[y]); calc(); // 随机替换贪心序列中任意两个元素，构造新状态
double delta=now-ans;
if (delta<0) ans=now;
else if (exp(-delta/T)*RAND_MAX>rand()) ans=now;
else swap(p[x], p[y]); // 交换回来，即不改变状态
T*=t0;
}
}
int main() {
srand(time(0)), srand(rand());
scanf("%d%d", &n, &m);
for (int i=1; i<=n; ++i) scanf("%d", &p[i]), ave+=p[i];
ave/=m;
calc(); ans=now;
mnth();
printf("%.2lf\n", ans);
return 0;
}
[USACO2013 Open] Haywire
Farmer John 有 \(N\) 只奶牛, (\(4 \le N \le 12\), 其中 \(N\) 是偶数). 他们建立了一套原生的系统，使得奶牛与他的朋友可以通过由干草保护的线路来进行对话交流. 每一头奶牛在这个牧场中正好有 3 个朋友，并且他们必须把自己安排在一排干草堆中. 一条长 \(L\) 的线路要占用刚好 \(N\) 堆干草来保护线路. 比如说，如果有两头奶牛分别在草堆 4 与草堆 7 中，并且他们是朋友关系，那么我们就需要用 3 堆干草来建造线路，使他们之间能够联系. 假设每一对作为朋友的奶牛都必须用一条单独的线来连接，并且我们可以随便地改变奶牛的位置，请计算出我们建造线路所需要的最少的干草堆.
此题也是随机序列的模拟退火。方法同上一题。
/* Haywire
* Au: GG
*/
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <ctime>
using namespace std;
const double eps=1e-15;
const double t0=0.998;
int n, G[13][3], f[13], ans, now;
inline int ABS(int x) {return x<0?-x:x; } // 卡常
inline void calc() {
now=0;
for (int i=1; i<=n; ++i) for (int j=0; j<3; ++j)
now+= ABS(f[i]-f[G[i][j]]);
}
inline void mnth() {
double T=200.0; // 数据范围小，不需要高温
while (T>eps) {
int x=rand()%n+1, y=rand()%n+1; if (x==y) y=rand()%n+1;
swap(f[x], f[y]), calc();
double delta=now-ans;
if (delta<0) ans=now;
else if (exp(-delta/T)*RAND_MAX>rand()) ans=now;
else swap(f[x], f[y]);
T*=t0;
}
}
int main() {
srand(time(0)), srand(rand()); // 玄之又玄，众妙之门
scanf("%d", &n);
for (int i=1; i<=n; ++i) for (int j=0; j<3; ++j) scanf("%d", &G[i][j]);
for (int i=1; i<=n; ++i) f[i]=i;
calc(); ans=now;
mnth();
printf("%d\n", ans/2);
return 0;
}
[NOIP2017TG] 宝藏
本题标算 状压DP！
观察题目描述和数据范围，显然，本题可以使用模拟退火算法。（所以是不是所有 状压DP/记忆化搜索 的题目都可以模拟退火乱搞？hhh……）
贪心地把序列中每一个点连到生成树的代价最小的节点后面，即局部最优解。此时本题转化为序列顺序的问题。模拟退火！
注意到本人提交此题 9 次才 AC，可见本题调参是比较有代表性的。（上面几题几乎一遍过……）
/* [NOIP2017TG] 宝藏
* Au: GG
*/
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <ctime>
using namespace std;
#define inf 0x3f3f3f3f
const double eps=1e-15;
const double t0=0.998;
int n, m;
int G[15][15], f[15], dep[15];
int ans, now;
inline void calc() { // 贪心方法是关键
now=0; memset(dep, 0, sizeof dep);
dep[f[1]]=1;
for (int i=2; i<=n; ++i) {
int w=inf;
for (int j=1, ww; j<i; ++j) if (G[f[j]][f[i]]<inf)
if ((ww=dep[f[j]]*G[f[j]][f[i]])<w) w=ww, dep[f[i]]=dep[f[j]]+1;
if (w<inf) now+=w; else {now=inf; return; }
}
}
inline void mnth() {
double T=10000.0; // 本题对精确度要求比较高，所以要加热到上万度
while (T>eps) {
int x=rand()%n+1, y=rand()%n+1; if (x==y) y=rand()%n+1;
swap(f[x], f[y]); calc();
double delta=now-ans;
if (delta<0) ans=now;
else if (exp(-delta/T)*RAND_MAX > rand()) ans=now;
else swap(f[x], f[y]);
T*=t0;
}
}
int main() {
srand(time(0)), srand(rand()+19260817); // 迷信一下
scanf("%d%d", &n, &m); memset(G, inf, sizeof G);
for (int i=1, a, b, c; i<=m; ++i)
scanf("%d%d%d", &a, &b, &c), G[a][b]=G[b][a]=min(G[a][b], c);
for (int i=1; i<=n; ++i) f[i]=i;
calc(); ans=now;
for (int i=1; i<=20; ++i) mnth(); // 连续退火多次保证正确率，当然次数太多也不行，
printf("%d\n", ans); // 100 次以上会 TLE（因为精度已经挺高了）。
return 0; // 目前此代码跑官方数据用时平均 146ms，极限 297ms，速度刚刚好。
}
推荐习题：[NOIP2016TG] 愤怒的小鸟（标算 搜索。试试看，你的乱搞能 AC 行吗？）
拓展
在不会 TLE 的情况下尽量多地跑 SA：
我们知道，有一个
clock()函数，返回程序运行时间。那么这样即可：while ((double)clock()/CLOCKS_PER_SEC<MAX_TIME) SA();其中
MAX_TIME是一个自定义的略小于 1 的正数，可以取 0.7~0.8。– M-sea
例：三角形牧场
本题标算 DP。
因为可以贪心，所以可以使用随机算法。
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <ctime>
using namespace std;
const double eps=1e-10; // 答案只要 2 位精度，所以 eps 可以调小
const double t0=0.998;
int n, L[43], sum, now, ans;
double heron(int x, int y, int z) {
double p=(x+y+z)/2.0;
return sqrt(p*(p-x)*(p-y)*(p-z));
}
inline void calc() { // 贪心
register int a=0.0, b=0.0, c=0.0, i=1;
for (; i<=n; ++i) if ((a+=L[i])*3.0>=sum) break;
for (++i; i<=n; ++i) if ((b+=L[i])*2.0>=sum-a) break;
c=sum-a-b;
if(a>=b+c||b>=a+c||c>=a+b||a<=0||b<=0||c<=0) now=-1;
now=100.0*heron(a, b, c);
}
inline void mnth() {
double T=5000.0;
int nocow=ans; // 缓存一下，解决 -1 的问题
while (T>eps) {
int x=rand()%n+1, y=rand()%n+1; while (y==x) y=rand()%n+1;
swap(L[x], L[y]); calc();
double delta=now-nocow;
if (now>0 && delta>=0) nocow=now; // 更优解
else if (exp(-delta/T)*RAND_MAX>rand()) nocow=now;
else swap(L[x], L[y]);
if (ans<nocow) ans=nocow; // 更新答案
T*=t0;
}
}
int main() {
srand(time(0)+19260817), srand(rand()^rand());
scanf("%d", &n);
for (int i=1; i<=n; ++i) scanf("%d", &L[i]), sum+=L[i];
random_shuffle(L+1, L+n+1);
calc(); ans=now;
while ((double)clock()/CLOCKS_PER_SEC<.9) mnth();
if (ans==(-2147483648)) printf("-1\n");
else printf("%d\n", ans);
return 0;
}
玄学算法拓展：
遗传算法（Genetic）：兔子们吃了失忆药片，并被发射到太空，然后随机落到了地球上的某些地方。他们不知道自己的使命是什么。但是，如果你过几年就杀死一部分海拔低的兔子，多产的兔子们自己就会找到珠穆朗玛峰。这就是遗传算法。
禁忌搜索算法（Tabu Search）：兔子们知道一个兔的力量是渺小的。他们互相转告着，哪里的山已经找过，并且找过的每一座山他们都留下一只兔子做记号。他们制定了下一步去哪里寻找的策略。这就是禁忌搜索。
Post author 作者: Grey
Copyright Notice 版权说明: Except where otherwise noted, all content of this blog is licensed under a CC BY-NC-SA 4.0 International license. 除非另有说明，本博客上的所有文章均受 知识共享署名 - 非商业性使用 - 相同方式共享 4.0 国际许可协议 保护。
