Dyd's Blog

He who has a strong enough why can bear almost any how.

模拟退火和爬山法

Posted on Edited on In
Symbols count in article: 7.3k Reading time ≈ 7 mins.

暴搜能过样例了！

模拟退火（Simulated Annealing）

思想

模拟退火是一种启发式的随机暴力算法，主要用于解决最优化问题。

基本概念

  • 温度：决定每次随机时我们随机的范围，温度越大，下一个点的随机区间越大，在模拟退火中，温度总是不断减小，一般把温度记为 $T$ ；
  • 初始温度：最初的温度值；
  • 终止温度：结束时的温度；
  • 衰减方式：每次温度的变化，可以是一次（如 $T_i=T_{i-1}-c$ )，也可以是指数级（如 $T_i=T_{i-1}*c,c\in (0,1)$ )等，根据实际情况而定，但一般要保证 $T_i<T_{i-1}$ ；
  • 能量值：当前答案的“权值”，或者说“优秀程度”，能量值越小，答案越优，个人喜欢记答案 $x$ 的能量值为 $E(X)$ ；

算法过程

  1. 确定初始点（初始答案），初始温度为 $T_0$ (一般保证足够大)，终止温度 $T_e$ （一般保证足够小）,衰减方式；
  2. 在当前点（设为 $now$ ）温度范围内随机选择一个值 $np$ ，计算答案 $E(np)$ ，记 $E(np)-E(now)=\Delta E$ ；
  3. 若 $\Delta E < 0$ ，说明 $np$ 比 $now$ 优，则跳转到 $np$ ;
  4. 若 $\Delta E > 0$ ，则以一定的概率跳转到 $np$ ，一般来说，这个概率取成 $e^{-\frac{\Delta E}{T}}$ ，这样取的好处是 $np$ 比 $now$ 越差，我们跳转过去的概率越小 ;

例题一

A Star not a Tree?

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#include<bits/stdc++.h>
using namespace std;
const int N=100+5,INF=0x3f3f3f3f;
const double T0=10000+5,TE=1e-4,C=0.99;
int n;
int T;
struct Node{
	double x,y;
	void inint(double _x,double _y){
		x=_x;
		y=_y;
	}
}q[N];
double ans=INF;
double rand_d(double l,double r){ //随机得到一个[l,r)的浮点数 
	return (double)rand()/RAND_MAX*(r-l)+l;
}
double get_dist(Node x,Node y){
	double dx=x.x-y.x;
	double dy=x.y-y.y;
	return sqrt(dx*dx+dy*dy);
}
double E(Node x){
	double res=0;
	for(int i=1;i<=n;++i) res+=get_dist(x,q[i]);
	ans=min(ans,res); //记得更新ans
	return res;
}
void simulate_anneal(){
	Node now;
	now.inint(rand_d(0,T0),rand_d(0,T0));
	for(double t=T0;t>TE;t*=C){
		Node np;
		np.inint(rand_d(now.x-t,now.x+t),rand_d(now.y-t,now.y+t));
		double dt=E(np)-E(now);
		if(exp(-dt/t)>rand_d(0,1)) now=np;
		/*
		这里exp(-dt/t)是求e的(-dt/t)次方,当dt<0时,其值恒大于1,一定会跳转到np
		当dt>0时,其值恒小于1,且dt越小,其值越大,跳转的概率也就越大 
		*/
	}
}
int main(){
	scanf("%d",&T);
	while(T--){
		ans=INF;
		scanf("%d",&n);
		for(int i=1;i<=n;++i) scanf("%lf%lf",&q[i].x,&q[i].y);
		for(int i=1;i<=100;++i) simulate_anneal(); //多次模拟退火使答案正确的概率变大
		printf("%.0lf\n",ans);
		if(T!=0) printf("\n"); 
	}
	return 0;
}

例题二

模拟退火不只能解决那些数学模型非常明显的题目，也可以解决一些其他最优化问题，只要保证答案具有“连续性”（即设 $\gamma$ 是大于 $0$ 且足够小的实数， $\forall x$ ，有 $|E(x)-E(x+ \gamma)|$ 足够小）。

看看这题保龄球

对于这种给数列 $A$ ，求 $A$ 的某种最优排列的问题，我们可以思考如何使用模拟退火：

  1. 温度可以直接像上题一样设置，衰减也可以自设；
  2. 初始点（初始排列）就用题目给的原排列，每次转移时在原排列上随机两个位置交换，得到下一个排列；
  3. 依题意打估价函数 $E(x)$ ，估价并概率转移；

由上，我们发现本题可以轻松用退火解决。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include<bits/stdc++.h>
using namespace std;
const int N=50+5,INF=0x3f3f3f3f;
const double T0=1e4+5,TE=1e-4,C=0.99;
int n,m;
int ans=-INF;
struct Try{
	int x,y;
}q[N];
int E(){
	int res=0;
	for(int i=1;i<=m;++i){
		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;
	/*
	"优"的答案,能量反而低
	当然也可以返回res,在simulate_anneal里使dt=x-y 
	*/ 
}
void simulate_anneal(){
	for(double t=T0;t>TE;t*=C){
		int a=rand()%m+1,b=rand()%m+1;
		int now=E();
		swap(q[a],q[b]);
		if(n+(q[n].x==10)==m){
			int np=E();
			int dt=np-now;
			if(exp(-dt/t)<(double)rand()/RAND_MAX) swap(q[a],q[b]);
			/*
			这里用小于号,代表若不跳转至np,就把a、b复原
			*/
		}
		else swap(q[a],q[b]);
	}
}
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;++i) scanf("%d%d",&q[i].x,&q[i].y);
	if(q[n].x==10) m=n+1,scanf("%d%d",&q[m].x,&q[m].y);
	else m=n;
	while((double)clock()/CLOCKS_PER_SEC<0.8) simulate_anneal();
	printf("%d\n",ans); 
	return 0;
}

例题三

一道类似于上题的排列退火，但加入一定的优化

均分数据

这道题直接模拟退火也能过，但略显的有点过于暴力，我们可以考虑模拟退火加一点贪心（当然也可以加DP，但那样还不如直接用DP，我就是做不出DP才模拟退火的好吧）：

先用退火得到数列的排列，计算当前答案时，对于当前数 $a_i$ ,我们每次都将 $a_i$ 放到当前和最小的集合中，贪心计算。可以证明，虽然我不会证，这样贪心是不会漏掉最优解的。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include<bits/stdc++.h>
using namespace std;
const int N=20+5,M=10,INF=0x3f3f3f3f;
const double T0=1e4+5,TE=1e-6,C=0.99;
int n,m;
int w[N],s[M];
double ans=INF;
double E(){
	for(int i=1;i<=m;++i) s[i]=0;
	for(int i=1,k=1;i<=n;++i){
		for(int j=1;j<=m;++j) if(s[j]<s[k]) k=j;
		s[k]+=w[i];
	}
	double avg=0;
	for(int i=1;i<=m;++i) avg+=(double)s[i]/m;
	double res=0;
	for(int i=1;i<=m;++i) res+=(s[i]-avg)*(s[i]-avg);
	res=sqrt(res/m);
	ans=min(ans,res);
	return res;
}
void simulate_anneal(){
	random_shuffle(w+1,w+1+n);
	for(double t=T0;t>TE;t*=C){
		int a=rand()%n+1,b=rand()%n+1;
		double x=E();
		swap(w[a],w[b]);
		double y=E();
		double dt=y-x;
		if(exp(-dt/t)<(double)rand()/RAND_MAX) swap(w[a],w[b]);
	}
}
int main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n;++i) scanf("%d",&w[i]);
	while((double)clock()/CLOCKS_PER_SEC<0.8) simulate_anneal();
	printf("%.2lf\n",ans);
	return 0;
}

以上两个题都用了 $while((double)clock()/CLOCKS_PER_SEC<0.8)$ 来控制退火次数，在实际做题时，一般计算好循环次数的代码会更稳定一些，但如果不好估算循环次数，这样也不失为一种选择。

爬山法

对于一个凸函数（单峰函数），我们可以用爬山法解决（当然，和退火一样，不是要求题目有明确的函数模型，只是需要答案对于每一维自变量你能感觉到一个单峰函数即可，有点玄学）。

爬山法的思路大体如图：

爬山

对于点 $now$ 我们并不关心它对应的具体值是多少，我们只关心它去向答案的“方向”（如图中应该是蓝色箭头方向），然后我们只需要向着“方向”走即可，图中是一个维度，而多个维度就分每个维度单独计算“方向”即可。当然，走的“距离”视具体情况而定。

需要注意的是，爬山法局限较多，不太常用，但仍不失为一种很好的启发式算法。

例题一

球形空间产生器

本题正解是高斯消元，建立方程联立求解即可，但如果你和我一样高斯消元和矩阵学的不好，也可以用爬山法搞定。

我们发现对于 $n$ 维坐标的每一维，都是一个单峰函数（球心到各点距离相等），所以对于一个可能答案的一维，求它到各点的距离的平均值，距离大于平均值的点表现出拉力，反之为推力，将所有力合成即可。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
#include<bits/stdc++.h>
using namespace std;
const int N=20+5;
int n;
double d[N][N];
double ans[N],dis[N],delta[N];
void calc(){
	double avg=0;
	for(int i=1;i<=n+1;++i){
		dis[i]=delta[i]=0;
		for(int j=1;j<=n;++j)
			dis[i]+=(d[i][j]-ans[j])*(d[i][j]-ans[j]);
		dis[i]=sqrt(dis[i]);
		avg+=dis[i]/(n+1);
	}
	for(int i=1;i<=n+1;++i)
		for(int j=0;j<=n;++j)
			delta[j]+=(dis[i]-avg)*(d[i][j]-ans[j])/avg;
} 
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n+1;++i)
		for(int j=1;j<=n;++j){
			scanf("%lf",&d[i][j]);
			ans[j]+=d[i][j]/(n+1); //先将初始答案定为重心 
		}
	for(double t=1e4;t>1e-6;t*=0.99997){
		calc();
		for(int i=1;i<=n;++i) ans[i]+=delta[i]*t;
	}
	for(int i=1;i<=n;++i) printf("%.3lf ",ans[i]);
	return 0;
}

当然，此题用模拟退火也行，代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#include<bits/stdc++.h>
using namespace std;
const int N=10+5,INF=0x3f3f3f3f;
const double T0=1e4+5,TE=1e-6,C=0.99995;
int n;
struct Node{
	double x[N];
	Node operator+(const Node &_t) {
	    Node res;
	    for(int i=1;i<=n;++i) res.x[i]=x[i]+_t.x[i];
	    return res;
  	}
}q[N],ans;
double as=INF;
double s[N];
Node now,np;
double rand_d(double l,double r){
	return (double)rand()/RAND_MAX*(r-l)+l;
}
double get_dist(Node x,Node y){
	double res=0;
	for(int i=1;i<=n;++i) res+=(x.x[i]-y.x[i])*(x.x[i]-y.x[i]);
	return sqrt(res);
}
double E(Node x){
	double avg=0;
	for(int i=1;i<=n+1;++i){
		s[i]=get_dist(x,q[i]);
		avg+=s[i];
	}
	avg=avg/(n+1);
	double res=0;
	for(int i=1;i<=n+1;++i) res+=(s[i]-avg)*(s[i]-avg);
	res*=avg;
	if(res<as){
		as=res;
		ans=x;
	}
	return res;
}
void simulate_anneal(){
	for(int i=1;i<=n+1;++i) now=now+q[i];
	for(int i=1;i<=n;++i) now.x[i]=now.x[i]/(n+1);
	for(double t=T0;t>TE;t*=C){
		for(int i=1;i<=n;++i) np.x[i]=now.x[i]+rand_d(-t,t);
		double dt=E(np)-E(now);
		if(exp(-dt/t)>rand_d(0,1)) now=np;
	}
}
int main(){
	srand(time(0));
	scanf("%d",&n);	
	for(int i=1;i<=n+1;++i)
		for(int j=1;j<=n;++j)
			scanf("%lf",&q[i].x[j]);
	simulate_anneal(); 
	/*
	时间有限,且答案对C要求较高,我们只退火一次,
	对应的,C调高
	*/ 
	for(int i=1;i<=n;++i) printf("%.3lf ",ans.x[i]);
	return 0;
}

牛顿迭代法

对于维度少的单峰函数，除了用三分法求最值外，也可用牛顿迭代法（又名牛顿爬山法）

牛顿迭代法的核心公式是
$$
x_{i+1}=x_i-\frac{f(x_i)}{f’(x_i)}
$$
但牛顿迭代法有一个巨大的缺陷——我这么弱根本求不出导函数来，笑死。