什么是模拟退火
模拟退火算法(Simulate Anneal,SA)是一种通用概率演算法,用来在一个大的搜寻空间内找寻命题的最优解。模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的。V.Černý在1985年也独立发明此演算法。模拟退火算法是解决TSP问题的有效方法之一。
简单说,模拟退火是一种随机化算法。当一个问题的方案数量极大(甚至是无穷的)而且不是一个单峰函数时,我们常使用模拟退火求解。
降温
实现过程
-
T
—— 温度 -
Δx
—— 解变动值 -
x1
—— 当前的目标解,等于x+Δx
-
Δf
—— 当前解的函数值与目标解函数值之差,等于f(x)−f(x1)
我们给一个初始解\(x\),并让它不断变动。要模拟变动的大小随温度的降低而降低,我们每次的\(\Delta x\)应该在一个大小与T成正比的值域内随机取值。
这时候我们就要考虑是否将当前解变为目标解。因为我们还是需要贪心,所以如果\(f(x1)<f(x)\),那么接受目标解,\(x=x1\)
那如果\(f(x1)>f(x)\)呢?我们当然要以一定概率接受它啦,这样才能跳出局部的限制,去搜寻更优的解,弥补贪心的局限性。那么这个概率应该是多少呢?同样要模拟变动的大小随温度的降低而降低。科学家经过理论分析,得出这个概率应该是\(e^{\frac{\Delta f}{T}}\)
小例子
首先经过大幅波动,当前解由\(A->B->C\),找到了一个比较满意的解。
但还不能满足。由于温度还比较大,此时接受了一些不比当前解优的目标解,C->D->E,成功地爬了上去。而温度还在慢慢减小。
终于,发现了一个更优解\(F\),成功跳出了那个非常宽的局部凹函数
这时候,温度越来越小了,很难再次接受不比当前解优的目标解,再次翻出去。解终于渐渐趋于稳定,并最终到达了\(G\),找到了最优解
如此看来,基于随机的模拟退火能很大程度上提高正确率,但也不可能完全正确。上面的例子只是随意模拟出的一个情况。所以要多跑几遍,取每一次得到的解的最优值。
例题
最短Hamilton路径
题目大意
给定一张\(n\)个点的带权无向图,点从\(0∼n−1\)标号,求起点\(0\)到终点\(n−1\)的最短 Hamilton 路径。(Hamilton 路径的定义是从\(0\)到\(n−1\)不重不漏地经过每个点恰好一次。)
本题的正解是:状态压缩DP
但是这题数据规模较小我们可以跑多次模拟退火极高的概率求出正解。
题目解空间:所有城市排列的集合
1.首先以初始的排列为解。
2.跑模拟退火交换相邻的节点
#include<bits/stdc++.h>
using namespace std;
const double eps=1e-15,delta=0.95;
const int N=100;
int dis[N][N];
int n;
int ans,last;
int w[N];
int fun(){
int res=0;
for(int i=0;i<n-1;i++)
res+=dis[w[i]][w[i+1]];
return res;
}
void sa(){
random_shuffle(w+1,w+n-1);
for(double t=1000;t>eps;t*=delta){
int last=fun();
int x=rand()%(n-2)+1,y=rand()%(n-2)+1;
swap(w[x],w[y]);
int now=fun();
int de=now-last;
if(de<0){
last=now;
ans=min(ans,now);
}
else if(exp(-1.0*de/t)/RAND_MAX>rand()) //以一定的概率接受劣解
last=now;
else
swap(w[x],w[y]); //回溯
}
}
int main(){
cin>>n;
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
cin>>dis[i][j];
for(int i=0;i<n;i++) w[i]=i;
ans=last=fun();
for(int i=0;i<500;i++) sa();
cout<<ans<<endl;
return 0;
}
星星还是树
题目大意:
在二维平面上有 n 个点,第 i 个点的坐标为 (xi,yi)。
请你找出一个点,使得该点到这 n 个点的距离之和最小。
该点可以选择在平面中的任意位置,甚至与这 n 个点的位置重合。
选取各个坐标的平均值为可行解。
由于一次模拟退火时间复杂度很小故而进行多次模拟退火选取最小的距离作为答案
#include<bits/stdc++.h>
using namespace std;
const double eps=1e-18,delta=0.98;
const int N=110;
struct node{
double x,y;
}p[N];
int n;
double ansx,ansy,ans;
double fun(double x,double y){
double res=0;
for(int i=0;i<n;i++){
double dx=p[i].x-x;
double dy=p[i].y-y;
res+=sqrt(dx*dx+dy*dy);
}
return res;
}
void sa(){
for(double t=5000;t>eps;t*=delta){
double nx=ansx+(rand()*2-RAND_MAX)*t;
double ny=ansx+(rand()*2-RAND_MAX)*t;
double nt=fun(nx,ny),de=nt-ans;
//cout<<nt<<endl;
if(de<0){
ansx=nx,ansy=ny;
ans=nt;
}
else if(exp(-1.0*de/t)*RAND_MAX>rand())
ansx=nx,ansy=ny;
}
}
int main(){
cin>>n;
ansx=0,ansy=0;
for(int i=0;i<n;i++){
double x,y;
cin>>x>>y;
p[i]={x,y};
ansx+=x,ansy+=y;
}
ansx/=n,ansy/=n;
ans=fun(ansx,ansy);
for(int i=0;i<100;i++) sa();
cout<<(int)(ans+0.5)<<endl;
return 0;
}
[TJOI2010]分金币
题目大意:
现在有\(n\)枚金币,它们可能会有不同的价值,第\(i\)枚金币的价值为 \(v_i\)
现在要把它们分成两部分,要求这两部分金币数目之差不超过1,问这样分成的两部分金币的价值之差最小是多少?
因为本题是要求平分,所以我们先锁边乱搞一下,分成两个组别,
每次在把一个组和另一个组的金币交换
序列对半交换,用模拟退火求解就可以了。
#include<bits/stdc++.h>
using namespace std;
const int N=1010;
int a[N];
const double eps=1e-15,delta=0.999;
int n,m;
#define ll long long
ll ans=LONG_LONG_MAX;
ll ansx;
int fun(){
ll res=0;
for(int i=0;i<n;i++) {
if(i<n/2) res+=a[i];
else res-=a[i];
}
ans=min(abs(res),ans);
return abs(res);
}
int Q;
void sa(){
ansx=fun();
for(double t=400;t>eps;t*=delta){
int x=rand()%n,y=rand()%n;
swap(a[x],a[y]);
ll now=fun();
int de=now-ansx;
if(de<0){ansx=now;}
else if(exp(-1.0*de/t)*RAND_MAX>rand()) ansx=now;
else swap(a[x],a[y]);
}
}
int main(){
int t;
cin>>t;
while(t--){
cin>>n;
for(int i=0;i<n;i++) cin>>a[i];
ans=LONG_LONG_MAX;
for(int i=0;i<10;i++) sa();
cout<<ans<<endl;
}
return 0;
}
平衡点 / 吊打XXX
我们可以确定一个原点,将所有的力在这个原点上正交分解,最终我们可以得到所有的力的一个合力,而真正的平衡点一定在合力所指向的方向
每当分得到一个合力之后,将原点在合力的方向上位移一定的距离。每当原点位移的方向发生了改变的时候,缩小以后操作的位移距离。例如:上次操作是将原点像 x 轴正方向移动,而当前移动是将原点像 x 轴负方向上移动,这说明原点的横坐标一定在这两次假设的原点的横坐标中间,因此我们缩小以后原点移动的距离。
#include<bits/stdc++.h>
using namespace std;
const int N=1e5+100;
struct node{
double x,y,w;
}p[N];
const double eps=1e-18,delta=0.9987;
int n;
double ansx,ansy,now;
double fun(double x,double y){
double sum=0,dx,dy;
for(int i=1;i<=n;i++){
dx=x-p[i].x,dy=y-p[i].y;
sum+=sqrt(dx*dx+dy*dy)*p[i].w;
}
return sum;
}
void sa(){
for(double t=4000;t>eps;t*=delta){
double nx=ansx+(rand()*2-RAND_MAX)*t;
double ny=ansy+(rand()*2-RAND_MAX)*t;
double nt=fun(nx,ny),de=nt-now;
if(de<0){
ansx=nx,ansy=ny,now=nt;
}
else if(exp(-de/t)*RAND_MAX>rand())
ansx=nx,ansy=ny;
}
}
int main(){
cin>>n;
for(int i=1;i<=n;i++){
double x,y,w;
cin>>x>>y>>w;
p[i]={x,y,w};
ansx+=w,ansy+=y;
}
ansx/=n,ansy=n;
now=fun(ansx,ansy);
sa();
printf("%.3lf %.3lf\n",ansx,ansy);
return 0;
}