模拟退火略解
模拟退火
模拟退火 算法（Simulate Anneal，SA）是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。
模拟退火的出发点是基于物理中 固体物质的退火过程 与一般组合优化问题之间的相似性。
固体物质的退火过程
lz前天的市物理竞赛爆炸（被防AK了）……心情极度不爽
听说oi同志们物理都很好 （翘物理竞赛辅导课） ，还是提一下这东西吧。
根据牛顿冷却定律 (Newton's law of cooling)，当物体表面与周围存在温度差时，单位时间从单位面积散失的热量与温度差成正比，即
下面我们以 \(\text{luogu P1337}\) 为例，讲解模拟退火算法的原理和步骤。
题目描述 \(\text{luoguP1337}\)
如图：有 \(n\) 个重物，每个重物系在一条足够长的绳子上。每条绳子自上而下穿过桌面上的洞，然后系在一起。图中 \(X\) 处就是公共的绳结。假设绳子是完全弹性的（不会造成能量损失），桌子足够高（因而重物不会垂到地上），且忽略所有的摩擦。
问绳结 \(X\) 最终平衡于何处。
注意：桌面上的洞都比绳结 \(X\) 小得多，所以即使某个重物特别重，绳结 \(X\) 也不可能穿过桌面上的洞掉下来，最多是卡在某个洞口处。
输入格式
文件的第一行为一个正整数 \(n\)（\(1≤n≤1000\)），表示重物和洞的数目。
接下来的 \(n\) 行，每行是 \(3\) 个整数：\(X_i.Y_i.W_i\)，分别表示第 \(i\) 个洞的坐标以及第 \(i\) 个重物的重量。\((-10000≤x,y≤10000, 0<w≤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→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≤n≤20000)\)。树从山顶到山脚按照 \(1,2,...,n\) 标号。
接下来 \(n\) 行，每行有两个正整数（用空格分开）。
第 \(i+1\) 行含有：\(w_i\) ——第 \(i\) 棵树的重量（公斤为单位）和 \(d_i\)——第 \(i\) 棵树和第 \(i+1\) 棵树之间的距离，\(1≤w_i≤10000,0≤d_i≤10000\)。
最后一颗树的 \(d_n\) 表示第 \(n\) 棵树到山脚的锯木厂的距离。保证所有树运到山脚的锯木厂所需要的费用小于 \(2×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{aligned}f(a,b)&=\sum_{i=1}{a}{w_i·(D_a-D_i)}+\sum_{i=a+1}+\sum_{i=b+1}^{n}{D_{n+1}-D_i}\
&=D_a·\sum_{i=1}{a}{w_i}+D_b·\sum_{i=a+1}+D_{n+1}·\sum_{i=b+1}{n}{D_i}-\sum_{i=1}
\end{aligned}$$
设$$W_k=\sum_{i=1}{k}{w_i}$$$$s_k=\sum_{i=1}$$我们发现它们可以使用前缀和优化，所以我们可以用 \(\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 | 已解决 |
浙公网安备 33010602011771号