模拟退火略解

模拟退火

模拟退火 算法（Simulate Anneal，SA）是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。
模拟退火的出发点是基于物理中 固体物质的退火过程 与一般组合优化问题之间的相似性。

固体物质的退火过程

lz前天的市物理竞赛爆炸（被防AK了）……心情极度不爽
听说oi同志们物理都很好 （翘物理竞赛辅导课） ，还是提一下这东西吧。

根据牛顿冷却定律 (Newton’s law of cooling)，当物体表面与周围存在温度差时，单位时间从单位面积散失的热量与温度差成正比，即
$$\Delta \tau =-k(\tau -{\tau }_0)$$其中 $k$ 为恒量。
式中 $\tau$ 表示物体温度， ${\tau }_0$ 表示环境温度（不变量）。
换句话说，单位时间散失的热量$$\Delta \tau \propto (\tau -{\tau }_0)$$



下面我们以 $\text{luogu P1337}$ 为例，讲解模拟退火算法的原理和步骤。

题目描述 $\text{luoguP1337}$

如图：有 $n$ 个重物，每个重物系在一条足够长的绳子上。每条绳子自上而下穿过桌面上的洞，然后系在一起。图中 $X$ 处就是公共的绳结。假设绳子是完全弹性的（不会造成能量损失），桌子足够高（因而重物不会垂到地上），且忽略所有的摩擦。

问绳结 $X$ 最终平衡于何处。

注意：桌面上的洞都比绳结 $X$ 小得多，所以即使某个重物特别重，绳结 $X$ 也不可能穿过桌面上的洞掉下来，最多是卡在某个洞口处。

输入格式

文件的第一行为一个正整数 $n$（$1\le n\le 1000$），表示重物和洞的数目。

接下来的 $n$ 行，每行是 $3$ 个整数：$X_i.Y_i.W_i$，分别表示第 $i$ 个洞的坐标以及第 $i$ 个重物的重量。$(-10000\le x,y\le 10000,0<w\le 1000)$

输出格式

你的程序必须输出两个浮点数（保留小数点后三位），分别表示处于最终平衡状态时绳结 $X$ 的横坐标和纵坐标。两个数以一个空格隔开。

输入样例

3
0 0 1
0 2 1
1 1 1

输出样例

0.577 1.000

$\text{Solution 1337}$

定义 $f(x,y)$ 表示 $X(x,y)$ 时的答案为多少。
容易得到空间能量和$$f(x,y)=\sum _{k=1}^n\sqrt{((X_k-x{)}^2+(Y_k-y{)}^2)}\times W_k$$
因为在一个空间中的能量和越小，该空间越稳定。所以 $f_{\min }$ 时的 $X$ 即为答案。


观察模拟退火时当前最优解与温度的关系。

可以看出，当 $\tau \to 0$ 时，我们找到最优解的几率也越来越大。
对于每个温度 $\tau$，尝试找一个新解。
若新解更优，则接受；若新解次，则以一定概率接受，这个概率为

$$e^{\frac{\Delta ans}{k\tau }}$$

其中 $k$ 是 $0$ 到 $1$ 之间的随机数。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<ctime>

#define reg register

int n;
struct node{
	int x,y,w;
}e[1010];
int cx=0,cy=0;
double ansx,ansy,ans=1e18;

double calc(double x,double y){	//计算f
	double cnt=0;
	for(reg int i=1;i<=n;++i)
		cnt+=sqrt((x-e[i].x)*(x-e[i].x)+(y-e[i].y)*(y-e[i].y))*e[i].w;
	return cnt;
}
void SA(){
	double x=ansx,y=ansy;
	double t=2000.0;	//初温
	while(t>1e-14){		//末温
		//得到新解坐标
		double nx=x+((rand()*2)-32767)*t;
		double ny=y+((rand()*2)-32767)*t;
		double nw=calc(nx,ny);
		double delta=nw-ans;	//delta ans
		if(delta<0){	//新解更优，接受
			x=nx;y=ny;ans=nw;
			ansx=nx;ansy=ny;
		}
		//新解更次，以一定概率接受
		else if(exp(-delta/t)*32767>rand()) x=nx,y=ny;
		t*=0.993;	//根据牛顿冷却定律降温
	}
}
void work(){
	double px=(double)cx/n,py=(double)cy/n;
	ansx=px;ansy=py;
	//多次执行SA，得到正确答案的几率更大
	for(reg int i=1;i<=10;++i) SA();
}
int main(){
	srand(100007);
	scanf("%d",&n);
	for(reg int i=1;i<=n;++i){
		scanf("%d%d%d",&e[i].x,&e[i].y,&e[i].w);
		cx+=e[i].x;cy+=e[i].y;
	}
	//从平均值开始的话，与答案更接近
	work();
	printf("%.3lf %.3lf",ansx,ansy);
}

算法改进

（1）设计合适的状态产生函数，使其根据搜索进程的需要表现出状态的全空间分散性或局部区域性；

（2）设计高效的退火策略；

（3）避免状态的迂回搜索；

（4）采用并行搜索结构；

（5）为避免陷入局部极小，改进对温度的控制方式；

（6）选择合适的初始状态；

（7）设计合适的算法终止准则。

题目描述 $\text{luogu P4360}$

从山顶上到山底下沿着一条直线种植了 $n$ 棵老树。当地的政府决定把他们砍下来。为了不浪费任何一棵木材，树被砍倒后要运送到锯木厂。

木材只能朝山下运。山脚下有一个锯木厂。另外两个锯木厂将新修建在山路上。你必须决定在哪里修建这两个锯木厂，使得运输的费用总和最小。假定运输每公斤木材每米需要一分钱。

你的任务是编写一个程序，从输入文件中读入树的个数和他们的重量与位置，计算最小运输费用。

输入格式

输入的第一行为一个正整数 $n$ ——树的个数 $(2\le n\le 20000)$。树从山顶到山脚按照 $1,2,...,n$ 标号。

接下来 $n$ 行，每行有两个正整数（用空格分开）。
第 $i+1$ 行含有：$w_i$ ——第 $i$ 棵树的重量（公斤为单位）和 $d_i$——第 $i$ 棵树和第 $i+1$ 棵树之间的距离，$1\le w_i\le 10000,0\le d_i\le 10000$。

最后一颗树的 $d_n$ 表示第 $n$ 棵树到山脚的锯木厂的距离。保证所有树运到山脚的锯木厂所需要的费用小于 $2\times 10^9$ 分。

输出格式

输出最小的运输费用。

输入样例

9 
1 2 
2 1 
3 3 
1 1 
3 2 
1 6 
2 1 
1 2 
1 1

输出样例

26

说明

样例图示

黑点为锯木厂。

$\text{Solution P4360}$

定义 $f(a,b)$ 表示锯木厂在 $a,b(a<b)$ 时的答案。
设第 $i$ 棵树到山脚的距离为 $D_i$，则$$D_i=\sum _{j=1}^{i-1}{d}_j$$
则$$\begin{array}{rl}f(a,b)& =\sum _{i=1}^a{w}_i\cdot (D_a-D_i)+\sum _{i=a+1}^b{w}_i\cdot (D_b-D_i)+\sum _{i=b+1}^n{D}_{n+1}-D_i\\ & =D_a\cdot \sum _{i=1}^a{w}_i+D_b\cdot \sum _{i=a+1}^b{w}_i+D_{n+1}\cdot \sum _{i=b+1}^n{D}_i-\sum _{i=1}^n{w}_i\cdot D_i\end{array}$$
设$$W_k=\sum _{i=1}^k{w}_i$$$$s_k=\sum _{i=1}^k{w}_i\cdot D_i$$我们发现它们可以使用前缀和优化，所以我们可以用 $\text{O(1)}$ 的时间复杂度求出 $f$。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<ctime>

#define reg register
typedef long long ll;

ll w[200010],d[200010],s;
int n;
struct node{
	int a,b;
	double f(){
		if(a>b) a^=b^=a^=b;
		return w[a]*d[a]+(w[b]-w[a])*d[b]+(w[n]-w[b])*d[n+1]-s;
	}
}x;
int cx=0,cy=0;
int ansx,ansy;
double ans=1e17;
int s1,s2;

void SA(){
	double a=1,b=2;
	double t=(double)n;
	while(t>5e-1){
		int nx=round(a+((rand()*2)-32767)*t);
		int ny=round(b+((rand()*2)-32767)*t);
		x.a=(nx+n)%n;x.b=(ny+n)%n;
		double nw=x.f();
		double delta=nw-ans;
		if(delta<0){
			a=nx;b=ny;ans=nw;
			ansx=(nx+n)%n;ansy=(ny+n)%n;
		}
		else if(exp(-delta/t)*32767>rand()) a=(nx+n)%n,b=(ny+n)%n;
		t*=0.99;
	}
}
void work(){
	for(reg int i=1;i<=10;++i) SA();
}
int main(){
	srand(100007);
	scanf("%d",&n);
	w[0]=d[0]=s=0;
	for(reg int i=1;i<=n;++i){
		scanf("%d%d",&s1,&s2);
		w[i]=w[i-1]+s1;d[i+1]=d[i+1]+s2;
		s=s+s1*d[i];
	}
	work();
	x.a=ansx;x.b=ansy;
	printf("%lld",(ll)x.f());
}

参考资料

序号	标题
$1$	浅谈玄学算法——模拟退火
$2$	题解 P4360 [CEOI2004]锯木厂选址

题表

序号	标题	题解
$1$	P2503 [HAOI2006]均分数据	未解决
$2$	P4035 [JSOI2008]球形空间产生器	已解决
$3$	P3878 [TJOI2010]分金币	已解决
$4$	SP4587 FENCE3 - Electric Fences	未解决
$5$	CF2C Commentator problem	未解决
$6$	UVA10228 A Star not a Tree?	已解决
$7$	P3936 Coloring	已解决
$8$	P2210 Haywire	已解决