算法笔记（3）模拟退火
原发表于个人博客=
模拟退火的引入
假如我们有一个函数，要求它的极大值，怎么求呢？
如果这个函数满足单调性，可以用二分的方法。
如果这是一个单谷（或单峰）函数，可以用三分法。
那要是多峰函数怎么半呢？
这时就可以用随机化算法。
一种朴素的方法是：每次在当前找到的最优方案\(x\)附近寻找一个新方案。如果这个新的解\(x'\)更优，那么转移到 \(x'\)，否则就不转移。（这被称作爬山算法）
但是这样就容易陷入一种局部最优解，如图，如果使用爬山算法，找到的答案可能在一个范围内是最优解，但是在全局上并不是最优解。
所以就有了模拟退火算法。
模拟退火算法的基本思路是，在一定的概率下接受一个非最优解而跳出局部最优解。
为什么这个算法被称为模拟退火呢？因为物理中金属降温是随机的，这个金属降温的过程也叫退火。
基本概念
刚刚讲到了模拟退火接受非最优解的概率，这个概率被称为\(Metropolis\)准则，也就是
其中，\(\Delta E\)表示最优解和当前解的差，\(T\)是当前温度。
模拟退火一般有三个参数，初始温度\(T_0\)，降温系数\(\Delta T\)（一般取\(0.99\)左右），最终温度\(T_{last}\)（一般是个极小值）。
当前温度\(T\)从\(T_0\)开始，不断乘以\(\Delta T\)，在这个过程中不断接受新解，直到\(T\)降至\(T_{last}\)，算法结束。
模板
这里以UVA10228（费马点问题）为例。
#include <bits/stdc++.h>
using namespace std;
const double eps=1e-8;
struct xy{double x,y;};
int n;
xy point[1010],anst,now;
double ans;
double calc(xy k){
//计算当前解距各点的距离
double ans=0;
for(int i=1;i<=n;i++) ans+=sqrt((point[i].x-k.x)*(point[i].x-k.x)+(point[i].y-k.y)*(point[i].y-k.y));
return ans;
}
void sa(){
for(double t=3000;t>eps;t*=0.999){
//降温过程
double cx=now.x+(rand()*2-RAND_MAX)*t;
double cy=now.y+(rand()*2-RAND_MAX)*t;//随机出新解
double k=calc({cx,cy});//计算新解
double da=k-ans;//计算新解与最优解的差
if(da<0){
//如果新解比最优解还好，直接接受
now=anst={cx,cy};
ans=k;
}
else if(exp(-da/t)*RAND_MAX>rand()) now={cx,cy};//否则根据概率接受
}
}
int main(){
srand(19260817);//随机种子
int t;
cin>>t;
while(t--){
cin>>n;
for(int i=1;i<=n;i++){
cin>>point[i].x>>point[i].y;
anst.x+=point[i].x,anst.y+=point[i].y;
}
anst.x/=(double)n,anst.y/=(double)n;
ans=calc(anst);//初始解（所有坐标的平均值）
for(int i=1;i<=100;i++) sa();//跑模拟退火
printf("%.0lf\n",calc(anst));//输出解
if(t) cout<<endl;
}
return 0;
}
应用
在OI中，模拟退火一般有三种用途，一种是计算几何（比如费马点问题），一种是序
列问题（随机交换），还有一种是骗分。
这篇博客的最后会有几道例题，并附讲解。
技巧
模拟退火本质上是一种随机化算法，能否正确取决于参数和rp。
这里提供几种技巧，在考试时可以使用。
卡时技巧
一般来说，模拟退火要跑很多遍才能跑出最优解，但是没遍模拟退火的时间我们也不知道，所以可以用\(clock\)函数，判断剩余时间，来卡时。
while((double)clock()/CLOCKS_PER_SEC<0.8) sa();
注：clock函数返回从程序开始到当前的时间，CLOCKS_PER_SEC将clock()函数的结果转化为以秒为单位的量，每个系统都不一样，在windows系统里是1000。
参数调小
一般来说，参数过大会WA，参数过小会TLE。
但是TLE的可能性一般比WA小，所以调参时宁小毋大，如平衡点一题，如果参数设置为\(10^{-8}\)，只能得到十分，如果设为\(10^{-14}\)，就能得到满分。
坐标的随机
在费马点问题中，假如当前坐标是\(x,y\)，如果随机到下一个坐标？
如果这么些：
double cx=x+rand(),cy=y+rand();
rand的取值是一个正值，所以坐标会逐渐递增，这显然是不行的。
所以我们可以这样写：
double cx=now.x+(rand()*2-RAND_MAX)*t;
double cy=now.y+(rand()*2-RAND_MAX)*t;
RAND_MAX是随机数的最大值，所以(rand()*2-RAND_MAX)的取值范围就是[-RAND_MAX,RAND_MAX]，再乘上当前的温度\(t\)，就会得到一个可能正负的随机浮点数，并且随着温度的减小越来越小，满足了我们的需求。
例题
P1337 平衡点 / 吊打XXX
这题需要一定的物理知识。
根据能量最低原理，一个系统内能量越低就越稳定，所以题目要求的平衡情况能量肯定最小。
分析这道题，如果静止不动，则能量只有重力势能，所以直接令重力势能最小即可。
我们知道，重力势能的公式是
\(g\)是个常量，所以直接算每个点的质量和高度乘积即可。
简化题意后，就是寻找一个点，让这个点到每个点的距离乘以质量最小，也就是一个带权费马点问题，套用模板即可。
AC代码
#include <bits/stdc++.h>
using namespace std;
const double eps=1e-15;
struct xy{double x,y;};
int n;
xy point[1010],anst,now;
double w[1010];
double ans;
inline double calc(xy k){
double ans=0;
for(int i=1;i<=n;i++) ans+=w[i]*sqrt((point[i].x-k.x)*(point[i].x-k.x)+(point[i].y-k.y)*(point[i].y-k.y));
return ans;
}
void sa(){
for(double t=3000;t>eps;t*=0.999){
double cx=now.x+(rand()*2-RAND_MAX)*t;
double cy=now.y+(rand()*2-RAND_MAX)*t;//随机出新解
double k=calc({cx,cy});//比较新解和最优解
double da=k-ans;//计算差值
if(da<0){
//如果新解比最优解好，直接接受
now=anst={cx,cy};
ans=k;
}
else if(exp(-da/t)*RAND_MAX>rand()) now={cx,cy};//否则按照概率接受
}
}
int main(){
srand(19260817);
cin>>n;
for(int i=1;i<=n;i++){
cin>>point[i].x>>point[i].y>>w[i];
anst.x+=point[i].x,anst.y+=point[i].y;
}
anst.x/=(double)n,anst.y/=(double)n;
ans=calc(anst);
while((double)clock()/CLOCKS_PER_SEC<0.8) sa();//卡时
printf("%.3lf %.3lf",anst.x,anst.y);
}
P5544 炸弹攻击1
大意是找到一个点，以它为圆心的圆能覆盖最多的点。
这题的关键是计算，首先这个半径不能超过\(r\)，其次这个圆不能与别的圆重合，找到符合这两个条件的最小的圆，计算它覆盖的点即可。
其它的部分套就是个费马点问题了，套模板即可
AC代码
#include <bits/stdc++.h>
using namespace std;
struct build{
double x,y,r;
}b[20];
struct xy{
double x,y;
}a[1010],ans;
double n,m,r,ansn=1e8;
inline double dis(xy a,xy b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int calc(xy k){
double minn=r,ans=0;
for(int i=1;i<=n;i++) minn=min(minn,dis(k,{b[i].x,b[i].y})-b[i].r);
for(int i=1;i<=m;i++) if(dis(k,a[i])<=minn) ans++;
return ans;
}
void sa(){
for(double T=3000;T>1e-16;T*=0.99){
double cx=ans.x+(rand()*2-RAND_MAX)*T,cy=ans.y+(rand()*2-RAND_MAX)*T;
double da=calc({cx,cy})-ansn;
if(da>0){
ans.x=cx,ans.y=cy;
ansn=calc({cx,cy});
}
else if(exp(da/T)>rand()) ans={cx,cy};
}
}
int main(){
srand(19260817);
cin>>n>>m>>r;
for(int i=1;i<=n;i++) cin>>b[i].x>>b[i].y>>b[i].r;
for(int i=1;i<=m;i++){
cin>>a[i].x>>a[i].y;
ans.x+=a[i].x,ans.y+=a[i].y;
}
ans.x/=m,ans.y/=m,ansn=calc(ans);
for(int i=1;i<=100;i++) sa();
cout<<ansn;
return 0;
}
P3878 分金币
以这道题为例，我们介绍一种模拟退火的其它应用。
本题大意是把\(n\)个数分成两半，让它们的和的差值最小。
模拟退火可以解决这种序列问题，以随机概率交换两个数，计算当前的差值，如果更优，就接受，否则以一定概率接受（如果不接受就把它们换回去）
AC代码
#include <bits/stdc++.h>
using namespace std;
const double q=0.999,eps=1e-14;
int n,a[60],ans=1e9;
int calc(){
int s1=0,s2=0;
for(int i=1;i<=(n+1)/2;i++) s1+=a[i];
for(int i=(n+1)/2+1;i<=n;i++) s2+=a[i];
return abs(s1-s2);
}
void sa(){
for(double T=5000;T>eps;T*=q){
int l=rand()%n+1,r=rand()%n+1;//取数列中随机两个数
swap(a[l],a[r]);//交换它们
int da=ans-calc();
if(da>0) ans-=da;
else if(exp(da/T)*RAND_MAX<rand()) swap(a[l],a[r]);
//这个意思是，如果不满足那个概率，就给他回复原状
}
}
int main(){
srand(19260817);
int t;
cin>>t;
while(t--){
ans=1e9;
cin>>n;
for(int i=1;i<=n;i++) cin>>a[i];
for(int i=1;i<=20;i++) sa();
cout<<ans<<endl;
}
}
