前言
模拟退火算法（Simulate Anneal，SA）是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的。V.Černý在1985年也独立发明此演算法。模拟退火算法是解决TSP问题的有效方法之一。
模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法是一种通用的优化算法，其物理退火过程由加温过程、等温过程、冷却过程这三部分组成
——百度百科
原来我是学打铁来了
抛开上面这一段，首先SA是个随机化算法。它能够解决什么样的问题呢？
另：友情提示，由于是随机化算法，所以需要一个随机数种子，本文所有代码的随机数种子已经得到了前辈的加成，可以放心使用。
正文
1.在下面这样的函数中找到最小值
由于不能通过枚举实数来模拟，那么就先考虑一下一个很显然的思路：贪心。
加入我们从最左端开始跳，逐渐往下跳到了这个位置：
接下来往左右跳都发现不对劲：都不能取得更优，所以就不会继续跳下去了。
但是肉眼可见，这段函数的最小值并不在这里。那么怎么办呢？模拟退火就来了。
首先要认识什么是退火。
这个退火，是指物理上的退火。退火是一种金属热处理工艺，指的是将金属缓慢加热到一定温度，保持足够时间，然后以适宜速度冷却。适宜的速度能使最终金属的内能最小。于是在处理最值问题的时候，我们能不能也找到这样适宜的速度呢？
再认识一下几个模拟退火中的参数：
(T_0)：模拟退火的初温，这个数值较大，比如(2000)。
(T_1)：模拟退火的末温，是一个略大于1的正数，比如(1e-14)。
(Δ)：降温系数，即温度变化率，是一个略小于1的正数，一般在([0.95,1))之间，比如(0.99)，这个数模拟的是缓缓降温。
我们就是通过一个变量(T)，先让它等于(T_0)，然后每次操作后(T*=Δ)，当(T<=T_1)的时候退火结束。
如何把这个温度和上面的问题结合起来呢？
在贪心的过程中，我们只单调地向目前看来的最优解跳，但是这样会被局限在一个局部最优解。所以我们需要一种手段来跳出这个局部最优解。但是呢又不能乱跳，否则快到最优解的时候又跳回局部最优解就很不合理。
来看一张演示图：（寻找函数最大值，摘自mediawiki）
可以看到，温度越低，波动幅度越小。并且在经过了更优解之后，依然会跳回不是那么优的解，但并非始终在乱跳。
那么就可以知道，模拟退火时可能会跳向目前看来并不是那么优的答案，随着温度的降低，答案也逐渐固定在了最优解周围。
对应到OI中，就是每次随机生成一个新的解，如果它比当前的答案更优，那么接受它。否则以一定的概率接受它。
而这个概率是多少呢？设新解与当前最优解的差是(Δx)，当前温度为(T)，生成随机数为k((Δx<0,T>0)）
这个概率就是(e^{frac{Δx}{k*T}})。但查了好久百度也只看到了这个。
到底为什么呢？让我们看下(e^x)的图像:
可以看到，当x<0,y∈(0,1)，并且这个x越大，y也越大。
对应到这个概率中，也就是在高温下，可以接受和当前最优解相差较大的解，而较低温度时则不行，所以跳到局部最优解也能跳出来，最终会逐渐稳定在最优解上。
这也解释了上面的“温度越低，波动幅度越小”，其实是因为在温度低时，更难接受与当前最优解相差大的解，反映到变化量上就是越来越小。
所以大概的代码框架是这样的：
const double delta=0.99;
double t;
int ans;
void SA()
{
t=2000;
while(t>1e-14)
{
(随机产生一组新解)
int now=计算这组解的价值
int del=now-ans;
if(del<0)接受这组解
else if(eps(-del/t)*RAND_MAX>rand())接受这组解
t*=delta;
}
}
当然，一次SA不一定，或者说在大数据下基本不能跑出最优解。那么我们就——多跑几次。
通过调用(double)clock()/CLOCKS_PER_SEC可以得到程序已经运行的秒数。那么...
while((double)clock()/CLOCKS_PER_SEC<=MAX_T)SA();这样就是最保险的了。MAX_T是自行设定的一个程序运行时间的上界。
另外，由于要生成随机数，所以在程序开始的时候可以srand(一个你喜欢的数)，比如114514,1919810,19****17之类的。
例题
前言
在例题的讲解之前，先解决一个问题：如何随机生成一组新解。下面是常见的几种方法：
平面系内：随机生成一对坐标。
序列问题：(random\_shuffle)，写法和sort无排序函数时一样，作用是随机重新排列序列。或随机交换两个数。
网格问题：随机生成两对坐标。
大致题意：把(n)枚价值分别为(w_i)的金币分成(n/2)和(n-(n/2))两堆，并使这两堆价值和的差最小。
(Tleq 20)组数据，(nleq 30,v_ileq 2^{30})。
正解貌似折半状压，但是退火永远不会畏惧。
这是我认为最适合新手练习模拟退火的一道题了。(虽然是道序列问题
初始就按照题意分为([1,n/2]),((n/2,n])两个区间。累加两个区间的和值到(s1,s2)两个变量中。
然后在SA中，每次降温时随机出两个位置x和y，并直接把它们交换并修改s1,s2。
如果修改后更优，那么当然接受，接受后如果答案更优就更新答案，否则按概率接受，但如果脸黑不接受的话换回去就好了。
一个小细节就是，如果n=1，需要特判一下，否则前一段长度是0，直接换炸。
至于SA的次数，这道题我是选用了5次。但明显的，随机数种子影响了得分。
所以，越臭的种子得分越高（暴论）
另外顺便说一下，(Δ)会对答案会有较大的影响，而(T_0)和(T_1)不会有较大的影响。不过也需要根据题目微调亿下。
于是一道紫色的dp题就变成了一道暴力题。
#include<bits/stdc++.h>
#define int long long
using namespace std;
inline int random(int a,int b){return (rand()%(b-a+1)+a);}
int a[50];
int ans=1e9,k,n,s1,s2,now_ans=1e9;
double t;
const double delta=0.995;
inline void exg(int x,int y)
{
s1=s1-a[x]+a[y];
s2=s2-a[y]+a[x];
swap(a[x],a[y]);
}
void SA()
{
now_ans=ans;
t=3000;
while(t>1e-14)
{
int x=random(1,k),y=random(k+1,n);
exg(x,y);
int now=abs(s1-s2);
int del=now-ans;
if(del<0){ans=min(ans,now);now_ans=now;}
else if(exp(-del/t)*RAND_MAX>rand()){now_ans=now;}
else exg(x,y);
t*=delta;
}
}
signed main()
{
int t;
scanf("%lld",&t);
while(t--)
{
ans=1e9;
srand(114514+1919810);
scanf("%lld",&n);
for(int i=1;i<=n;i++)
scanf("%lld",&a[i]);
k=n>>1;
s1=0,s2=0;
for(int i=1;i<=k;i++)s1+=a[i];
for(int i=k+1;i<=n;i++)s2+=a[i];
ans=abs(s1-s2);
if(n==1){printf("%lld
",ans);continue;}
SA();SA();SA();SA();SA();
printf("%lld
",ans);
}
return 0;
}
这道题被各路大佬都拿来当模拟退火入门题，但看见实数就难受的我只能把它放在第二位。
大致题意：给定(n)个点的(x_i,y_i,w_i)（横纵坐标，重量），求找到一个点(A(x,y))使得(sum^{n}_{i=1}w_i*dis_{A,i})最小。保留三位小数。
(1leq n,wleq 1000)，坐标在([-10000,10000])内。
对于给定的((x,y))，计算贡献时我们直接暴力的(O(n))计算即可。
那么依然随机的生成一对点再计算答案，和之前的降温过程没有区别。
#include<bits/stdc++.h>
using namespace std;
struct node{int x,y,w;}a[1010];
double ansx,ansy;
double ans=1e18,t;
const double delta=0.995;
int n;
inline double count(double x,double y)
{
double res=0;
for(int i=1;i<=n;i++)
{
double dx=x-a[i].x,dy=y-a[i].y;
res+=sqrt(dx*dx+dy*dy)*a[i].w;
}
return res;
}
inline void SA()
{
double x=ansx,y=ansy;
t=2000;
while(t>1e-14)
{
double dx=x+((rand()<<1)-RAND_MAX)*t;
double dy=y+((rand()<<1)-RAND_MAX)*t;
double now=count(dx,dy);
double del=now-ans;
if(del<0){x=dx,y=dy,ans=now;ansx=x,ansy=y;}
else if(exp(-del/t)*RAND_MAX>rand()){x=dx,y=dy;}
t*=delta;
}
}
int main()
{
int sx=0,sy=0;
srand(114514+1919810);
scanf("%d",&n);
for(int i=1;i<=n;i++)
{
scanf("%d%d%d",&a[i].x,&a[i].y,&a[i].w);
sx+=a[i].x,sy+=a[i].y;
}
ansx=(double)sx/n,ansy=(double)sy/n;
SA();SA();
printf("%.3lf %.3lf
",ansx,ansy);
return 0;
}
大致题意：
(nleq 50)个城市由(n)条边(u_i,v_i,d_i)连通，其中已有(mleq n)个城市建立了城堡，还需要再建(k)座城堡在未建立城堡的城市。令(dis(c))表示城市(c)到最近的城堡的距离。求确定这(k)座城堡的建立位置使得(max{dis(c)})最小。
看起来变成了图论，但是(nleq 50)就说明这不是一般的图论，那么考虑每次用dij跑出所有(dis(i))，再更新答案就可以了。思路和第一题一模一样，也只需要随机把一个选入这k个城市中的一个换出去，再随便换一个进来即可。
所以依然是暴力啊。。
特别的，如果m+k=n，没得选，直接输出一次dij后的答案即可。
#include<bits/stdc++.h>
#define pi pair<int,int>
#define v first
#define u second
using namespace std;
const int N=110;
int head[N],to[N],nxt[N],val[N],cnt;
int g[N],c[N],pos[N];
int d[N],vis[N],add[N];
int n,m,k,di,p,k0=0,ans=0;
void dij()
{
memset(d,0x7f,sizeof d);
memset(vis,0,sizeof vis);
priority_queue<pi> q;
for(int i=1;i<=m;i++)d[pos[i]]=0,q.push(pi(0,pos[i]));
for(int i=1;i<=k;i++)d[add[i]]=0,q.push(pi(0,add[i]));
while(!q.empty())
{
int u=q.top().u;
q.pop();
if(vis[u])continue;
vis[u]=1;
for(int i=head[u];i;i=nxt[i])
{
int v=to[i];
if(d[v]>d[u]+val[i])
{
d[v]=d[u]+val[i];
q.push(pi(-d[v],v));//大根堆->小根堆
}
}
}
}
inline int count(){int mx=0;dij();for(int i=1;i<=n;i++)mx=max(mx,d[i]);return mx;}
const double MAX_T=0.6;
const double delta=0.996;
double t;
inline int random(int a,int b){return rand()%(b-a+1)+a;}
inline void SA()
{
t=2000;
while(t>1e-14)
{
int x=random(1,k),y=random(k+1,k0);
if(add[x]==add[y])continue;
swap(add[x],add[y]);
int now=count();
int del=now-ans;
if(del<0)ans=now;
else if(exp(-del/t)*RAND_MAX<=rand())swap(add[x],add[y]);
t*=delta;
}
}
void adde(int u,int v,int w)
{
cnt++;
to[cnt]=v;
val[cnt]=w;
nxt[cnt]=head[u];
head[u]=cnt;
}
int main()
{
srand(114514+1919810);
scanf("%d%d%d",&n,&m,&k);
for(int i=1;i<=n;i++){scanf("%d",&g[i]);g[i]++;}
for(int i=1;i<=n;i++){scanf("%d",&di);adde(i,g[i],di);adde(g[i],i,di);}
for(int i=1;i<=m;i++){scanf("%d",&pos[i]);pos[i]++;c[pos[i]]=1;}
dij();for(int i=1;i<=n;i++)ans=max(ans,d[i]);
for(int i=1;i<=n;i++)if(!c[i])add[++k0]=i;
if(m+k==n){printf("%d
",count());return 0;}
while((double)clock()/CLOCKS_PER_SEC<=MAX_T)SA();
printf("%d
",ans);
return 0;
}
第二道黑题
update:2020.10.7
无聊又做了一道模拟退火的题目，来更新。
大致题意：将(1leq nleq 20)个数(1leq a_ileq 50)分成(2leq mleq 6)组，使得每组的和值(x_i)的均方差最小。
即(ar x=frac{1}{n}sum^n_{i=1}a_i)，求(sqrt{frac{1}{n}sum^m_{i-1}(ar x-x_i)^2})最小值。
首先考虑给定一个序列a,如何把它分组使得这个均方差最小。思路很简单，我们可以通过贪心地把1~n里的数一次放进目前和值最小的组里，由于数据太小所以这个操作直接(O(nm))就好。
然后就是一个序列问题模拟退火的板子，每次随机交换两个数并计算新答案，更优则接受，否则则以一定概率接受。随机数种子114514+1919810即可。
#include<bits/stdc++.h>
#define p(x) ((x)*(x))
using namespace std;
int x[10];
int a[30];
int n,m,s;
double xj;
inline double count()
{
memset(x,0,sizeof x);
for(int i=1;i<=n;i++)
{
int k=1;
for(int j=1;j<=m;j++)
if(x[j]<x[k])k=j;
x[k]+=a[i];
}
double now=0;
for(int i=1;i<=m;i++)now+=(x[i]-xj)*(x[i]-xj);
return sqrt(now/(double)m);
}
const double delta=0.99;
double t;
double ans=0;
void SA()
{
t=2000;
while(t>1e-14)
{
int x=rand()%n+1,y=rand()%n+1;
while(x==y)y=rand()%n+1;
swap(a[x],a[y]);
double now=count();
double del=now-ans;
if(del<0)ans=now;
else if(exp(-del/t)*RAND_MAX<=rand())swap(a[x],a[y]);
t*=delta;
}
}
double MAX_T=0.8;
int main()
{
srand(114514+1919810);
scanf("%d%d",&n,&m);
for(int i=1;i<=n;i++)
{
scanf("%d",&a[i]);
s+=a[i];
}
xj=(double)s/m;
ans=count();
if(n==2){printf("%.2lf
",ans);return 0;}
while((double)clock()/CLOCKS_PER_SEC<=MAX_T)SA();
printf("%.2lf
",ans);
return 0;
}
最后是一道不属于模拟退火的例题，但是..可以用随机化贪心水过。
什么是随机化贪心？
题意自己感性理解一下。具体操作就是random_shuffle之后统计答案即可。
为什么不枚举全排列？当然是因为时间不够。
为什么不能过？我/你是非酋。
至于是否能跑出正解呢？跑够了肯定是可以的，所以我跑了1e6次就过了。建议先洗脸后提交。
#include<bits/stdc++.h>
using namespace std;
int a[13][3],b[13],ans=1e9,n;
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++)
scanf("%d%d%d",&a[i][0],&a[i][1],&a[i][2]),b[i]=i;
for(int i=1;i<=1000000;i++)
{
random_shuffle(b+1,b+1+n);
int res=0;
for(int i=1;i<=n;i++)
for(int j=0;j<3;j++)
res+=(abs(b[i]-b[a[i][j]]));
ans=min(ans,res>>1);
}
cout<<ans<<endl;
return 0;
}
遇到其他类型的题目再丕定更新。