模拟退火略解
模拟退火
模拟退火 算法(Simulate Anneal,SA)是一种通用概率演算法,用来在一个大的搜寻空间内找寻命题的最优解。
模拟退火的出发点是基于物理中 固体物质的退火过程 与一般组合优化问题之间的相似性。
固体物质的退火过程
lz前天的市物理竞赛爆炸(被防AK了)……心情极度不爽
听说oi同志们物理都很好 (翘物理竞赛辅导课) ,还是提一下这东西吧。
根据牛顿冷却定律 (Newton’s law of cooling),当物体表面与周围存在温度差时,单位时间从单位面积散失的热量与温度差成正比,即
Δ τ = − k ( τ − τ 0 ) \Delta\tau=-k(\tau-\tau_0) Δτ=−k(τ−τ0)其中 k k k 为恒量。
式中 τ \tau τ 表示物体温度, τ 0 \tau_0 τ0 表示环境温度(不变量)。
换句话说,单位时间散失的热量 Δ τ ∝ ( τ − τ 0 ) \Delta\tau\propto(\tau-\tau_0) Δτ∝(τ−τ0)
下面我们以 luogu P1337 \text{luogu P1337} luogu P1337 为例,讲解模拟退火算法的原理和步骤。
题目描述 luoguP1337 \text{luoguP1337} luoguP1337
如图:有 n n n 个重物,每个重物系在一条足够长的绳子上。每条绳子自上而下穿过桌面上的洞,然后系在一起。图中 X X X 处就是公共的绳结。假设绳子是完全弹性的(不会造成能量损失),桌子足够高(因而重物不会垂到地上),且忽略所有的摩擦。
问绳结 X X X 最终平衡于何处。
注意:桌面上的洞都比绳结 X X X 小得多,所以即使某个重物特别重,绳结 X X X 也不可能穿过桌面上的洞掉下来,最多是卡在某个洞口处。
输入格式
文件的第一行为一个正整数 n n n( 1 ≤ n ≤ 1000 1≤n≤1000 1≤n≤1000),表示重物和洞的数目。
接下来的 n n n 行,每行是 3 3 3 个整数: X i . Y i . W i X_i.Y_i.W_i Xi.Yi.Wi,分别表示第 i i i 个洞的坐标以及第 i i i 个重物的重量。 ( − 10000 ≤ x , y ≤ 10000 , 0 < w ≤ 1000 ) (-10000≤x,y≤10000, 0<w≤1000 ) (−10000≤x,y≤10000,0<w≤1000)
输出格式
你的程序必须输出两个浮点数(保留小数点后三位),分别表示处于最终平衡状态时绳结 X X X 的横坐标和纵坐标。两个数以一个空格隔开。
输入样例
3
0 0 1
0 2 1
1 1 1
输出样例
0.577 1.000
Solution 1337 \text{Solution 1337} Solution 1337
定义 f ( x , y ) f(x,y) f(x,y) 表示 X ( x , y ) X(x,y) X(x,y) 时的答案为多少。
容易得到空间能量和 f ( x , y ) = ∑ k = 1 n ( ( X k − x ) 2 + ( Y k − y ) 2 ) × W k f(x,y)=\sum_{k=1}^{n}{\sqrt{((X_k-x)^2+(Y_k-y)^2)}\times W_k} f(x,y)=k=1∑n((Xk−x)2+(Yk−y)2) ×Wk
因为在一个空间中的能量和越小,该空间越稳定。所以 f min f_{\min} fmin 时的 X X X 即为答案。
观察模拟退火时当前最优解与温度的关系。
可以看出,当 τ → 0 \tau→0 τ→0 时,我们找到最优解的几率也越来越大。
对于每个温度 τ \tau τ,尝试找一个新解。
若新解更优,则接受;若新解次,则以一定概率接受,这个概率为
e Δ a n s k τ e^{\frac{\Delta ans}{k\tau}} ekτΔans
其中 k k k 是 0 0 0 到 1 1 1 之间的随机数。
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<ctime>#define reg registerint 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){ //计算fdouble 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 ansif(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)设计合适的算法终止准则。
题目描述 luogu P4360 \text{luogu P4360} luogu P4360
从山顶上到山底下沿着一条直线种植了 n n n 棵老树。当地的政府决定把他们砍下来。为了不浪费任何一棵木材,树被砍倒后要运送到锯木厂。
木材只能朝山下运。山脚下有一个锯木厂。另外两个锯木厂将新修建在山路上。你必须决定在哪里修建这两个锯木厂,使得运输的费用总和最小。假定运输每公斤木材每米需要一分钱。
你的任务是编写一个程序,从输入文件中读入树的个数和他们的重量与位置,计算最小运输费用。
输入格式
输入的第一行为一个正整数 n n n ——树的个数 ( 2 ≤ n ≤ 20000 ) (2≤n≤20000) (2≤n≤20000)。树从山顶到山脚按照 1 , 2 , . . . , n 1,2,...,n 1,2,...,n 标号。
接下来 n n n 行,每行有两个正整数(用空格分开)。
第 i + 1 i+1 i+1 行含有: w i w_i wi ——第 i i i 棵树的重量(公斤为单位)和 d i d_i di——第 i i i 棵树和第 i + 1 i+1 i+1 棵树之间的距离, 1 ≤ w i ≤ 10000 , 0 ≤ d i ≤ 10000 1≤w_i≤10000,0≤d_i≤10000 1≤wi≤10000,0≤di≤10000。
最后一颗树的 d n d_n dn 表示第 n n n 棵树到山脚的锯木厂的距离。保证所有树运到山脚的锯木厂所需要的费用小于 2 × 1 0 9 2×10^9 2×109 分。
输出格式
输出最小的运输费用。
输入样例
9
1 2
2 1
3 3
1 1
3 2
1 6
2 1
1 2
1 1
输出样例
26
说明
样例图示
黑点为锯木厂。
Solution P4360 \text{Solution P4360} Solution P4360
定义 f ( a , b ) f(a,b) f(a,b) 表示锯木厂在 a , b ( a < b ) a,b\ (a<b) a,b (a<b) 时的答案。
设第 i i i 棵树到山脚的距离为 D i D_i Di,则 D i = ∑ j = 1 i − 1 d j D_i=\sum_{j=1}^{i-1}{d_j} Di=j=1∑i−1dj
则 f ( a , b ) = ∑ i = 1 a w i ⋅ ( D a − D i ) + ∑ i = a + 1 b w i ⋅ ( D b − D i ) + ∑ i = b + 1 n D n + 1 − D i = D a ⋅ ∑ i = 1 a w i + D b ⋅ ∑ i = a + 1 b w i + D n + 1 ⋅ ∑ i = b + 1 n D i − ∑ i = 1 n w i ⋅ D i \begin{aligned}f(a,b)&=\sum_{i=1}^{a}{w_i·(D_a-D_i)}+\sum_{i=a+1}^{b}{w_i·(D_b-D_i)}+\sum_{i=b+1}^{n}{D_{n+1}-D_i}\\ &=D_a·\sum_{i=1}^{a}{w_i}+D_b·\sum_{i=a+1}^{b}{w_i}+D_{n+1}·\sum_{i=b+1}^{n}{D_i}-\sum_{i=1}^{n}{w_i·D_i} \end{aligned} f(a,b)=i=1∑awi⋅(Da−Di)+i=a+1∑bwi⋅(Db−Di)+i=b+1∑nDn+1−Di=Da⋅i=1∑awi+Db⋅i=a+1∑bwi+Dn+1⋅i=b+1∑nDi−i=1∑nwi⋅Di
设 W k = ∑ i = 1 k w i W_k=\sum_{i=1}^{k}{w_i} Wk=i=1∑kwi s k = ∑ i = 1 k w i ⋅ D i s_k=\sum_{i=1}^{k}{w_i·D_i} sk=i=1∑kwi⋅Di我们发现它们可以使用前缀和优化,所以我们可以用 O(1) \text{O(1)} O(1) 的时间复杂度求出 f f 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 1 1 | 浅谈玄学算法——模拟退火 |
2 2 2 | 题解 P4360 [CEOI2004]锯木厂选址 |
题表
序号 | 标题 | 题解 |
---|---|---|
1 1 1 | P2503 [HAOI2006]均分数据 | 未解决 |
2 2 2 | P4035 [JSOI2008]球形空间产生器 | 已解决 |
3 3 3 | P3878 [TJOI2010]分金币 | 已解决 |
4 4 4 | SP4587 FENCE3 - Electric Fences | 未解决 |
5 5 5 | CF2C Commentator problem | 未解决 |
6 6 6 | UVA10228 A Star not a Tree? | 已解决 |
7 7 7 | P3936 Coloring | 已解决 |
8 8 8 | P2210 Haywire | 已解决 |
更多推荐
模拟退火略解
发布评论