首页 > 编程语言 >算法笔记(3)模拟退火

算法笔记(3)模拟退火

时间:2023-10-23 18:25:32浏览次数:43  
标签:rand RAND int double 笔记 算法 模拟退火 ans

原发表于个人博客=

模拟退火的引入

假如我们有一个函数,要求它的极大值,怎么求呢?

如果这个函数满足单调性,可以用二分的方法。

如果这是一个单谷(或单峰)函数,可以用三分法。

那要是多峰函数怎么半呢?

这时就可以用随机化算法。

一种朴素的方法是:每次在当前找到的最优方案\(x\)附近寻找一个新方案。如果这个新的解\(x'\)更优,那么转移到 \(x'\),否则就不转移。(这被称作爬山算法)

但是这样就容易陷入一种局部最优解,如图,如果使用爬山算法,找到的答案可能在一个范围内是最优解,但是在全局上并不是最优解。
图片来自oiwiki
所以就有了模拟退火算法。

模拟退火算法的基本思路是,在一定的概率下接受一个非最优解而跳出局部最优解。
百年老图

为什么这个算法被称为模拟退火呢?因为物理中金属降温是随机的,这个金属降温的过程也叫退火。

基本概念

刚刚讲到了模拟退火接受非最优解的概率,这个概率被称为\(Metropolis\)准则,也就是

\[e^{\frac{-\Delta E}{T}} \]

其中,\(\Delta E\)表示最优解和当前解的差,\(T\)是当前温度。

模拟退火一般有三个参数,初始温度\(T_0\),降温系数\(\Delta T\)(一般取\(0.99\)左右),最终温度\(T_{last}\)(一般是个极小值)。

当前温度\(T\)从\(T_0\)开始,不断乘以\(\Delta T\),在这个过程中不断接受新解,直到\(T\)降至\(T_{last}\),算法结束。

模板

这里以UVA10228(费马点问题)为例。

#include <bits/stdc++.h>
using namespace std;
const double eps=1e-8;
struct xy{double x,y;};
int n;
xy point[1010],anst,now;
double ans;
double calc(xy k){
    //计算当前解距各点的距离
	double ans=0;
	for(int i=1;i<=n;i++)  ans+=sqrt((point[i].x-k.x)*(point[i].x-k.x)+(point[i].y-k.y)*(point[i].y-k.y));
	return ans;
}
void sa(){
	for(double t=3000;t>eps;t*=0.999){
		//降温过程
		double cx=now.x+(rand()*2-RAND_MAX)*t;
		double cy=now.y+(rand()*2-RAND_MAX)*t;//随机出新解
		double k=calc({cx,cy});//计算新解
		double da=k-ans;//计算新解与最优解的差
		if(da<0){
		    //如果新解比最优解还好,直接接受
			now=anst={cx,cy};
			ans=k;
		}
		else if(exp(-da/t)*RAND_MAX>rand()) now={cx,cy};//否则根据概率接受
	}
}
int main(){
	srand(19260817);//随机种子
	int t;
	cin>>t;
	while(t--){
    	cin>>n;
    	for(int i=1;i<=n;i++){
    		cin>>point[i].x>>point[i].y;
    		anst.x+=point[i].x,anst.y+=point[i].y;
    	}
    	anst.x/=(double)n,anst.y/=(double)n;
    	ans=calc(anst);//初始解(所有坐标的平均值)
    	for(int i=1;i<=100;i++) sa();//跑模拟退火
    	printf("%.0lf\n",calc(anst));//输出解
    	if(t) cout<<endl;
	}
	return 0;
}

应用

在OI中,模拟退火一般有三种用途,一种是计算几何(比如费马点问题),一种是序
列问题(随机交换),还有一种是骗分。

这篇博客的最后会有几道例题,并附讲解。

技巧

模拟退火本质上是一种随机化算法,能否正确取决于参数和rp。

这里提供几种技巧,在考试时可以使用。

卡时技巧

一般来说,模拟退火要跑很多遍才能跑出最优解,但是没遍模拟退火的时间我们也不知道,所以可以用\(clock\)函数,判断剩余时间,来卡时。

while((double)clock()/CLOCKS_PER_SEC<0.8) sa();

注:clock函数返回从程序开始到当前的时间,CLOCKS_PER_SECclock()函数的结果转化为以秒为单位的量,每个系统都不一样,在windows系统里是1000。

参数调小

一般来说,参数过大会WA,参数过小会TLE。

但是TLE的可能性一般比WA小,所以调参时宁小毋大,如平衡点一题,如果参数设置为\(10^{-8}\),只能得到十分,如果设为\(10^{-14}\),就能得到满分。

坐标的随机

在费马点问题中,假如当前坐标是\(x,y\),如果随机到下一个坐标?

如果这么些:

double cx=x+rand(),cy=y+rand();

rand的取值是一个正值,所以坐标会逐渐递增,这显然是不行的。

所以我们可以这样写:

double cx=now.x+(rand()*2-RAND_MAX)*t;
double cy=now.y+(rand()*2-RAND_MAX)*t;		

RAND_MAX是随机数的最大值,所以(rand()*2-RAND_MAX)的取值范围就是[-RAND_MAX,RAND_MAX],再乘上当前的温度\(t\),就会得到一个可能正负的随机浮点数,并且随着温度的减小越来越小,满足了我们的需求。

例题

P1337 平衡点 / 吊打XXX

原题

这题需要一定的物理知识。

根据能量最低原理,一个系统内能量越低就越稳定,所以题目要求的平衡情况能量肯定最小。

分析这道题,如果静止不动,则能量只有重力势能,所以直接令重力势能最小即可。

我们知道,重力势能的公式是

\[E_p=mgh \]

\(g\)是个常量,所以直接算每个点的质量和高度乘积即可。

简化题意后,就是寻找一个点,让这个点到每个点的距离乘以质量最小,也就是一个带权费马点问题,套用模板即可。

AC代码

#include <bits/stdc++.h>
using namespace std;
const double eps=1e-15;
struct xy{double x,y;};
int n;
xy point[1010],anst,now;
double w[1010];
double ans;
inline double calc(xy k){
	double ans=0;
	for(int i=1;i<=n;i++)  ans+=w[i]*sqrt((point[i].x-k.x)*(point[i].x-k.x)+(point[i].y-k.y)*(point[i].y-k.y));
	return ans;
}
void sa(){
    for(double t=3000;t>eps;t*=0.999){
		double cx=now.x+(rand()*2-RAND_MAX)*t;
		double cy=now.y+(rand()*2-RAND_MAX)*t;//随机出新解
		double k=calc({cx,cy});//比较新解和最优解
		double da=k-ans;//计算差值
		if(da<0){
		    //如果新解比最优解好,直接接受
			now=anst={cx,cy};
			ans=k;
		}
		else if(exp(-da/t)*RAND_MAX>rand()) now={cx,cy};//否则按照概率接受
	}
}
int main(){
	srand(19260817);
	cin>>n;
	for(int i=1;i<=n;i++){
	    cin>>point[i].x>>point[i].y>>w[i];
	    anst.x+=point[i].x,anst.y+=point[i].y;
	}
	anst.x/=(double)n,anst.y/=(double)n;
	ans=calc(anst);
	while((double)clock()/CLOCKS_PER_SEC<0.8) sa();//卡时
	printf("%.3lf %.3lf",anst.x,anst.y);
}

P5544 炸弹攻击1

原题

大意是找到一个点,以它为圆心的圆能覆盖最多的点。

这题的关键是计算,首先这个半径不能超过\(r\),其次这个圆不能与别的圆重合,找到符合这两个条件的最小的圆,计算它覆盖的点即可。

其它的部分套就是个费马点问题了,套模板即可

AC代码

#include <bits/stdc++.h>
using namespace std;
struct build{
    double x,y,r;
}b[20];
struct xy{
    double x,y;
}a[1010],ans;
double n,m,r,ansn=1e8;
inline double dis(xy a,xy b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int calc(xy k){
	double minn=r,ans=0;
	for(int i=1;i<=n;i++) minn=min(minn,dis(k,{b[i].x,b[i].y})-b[i].r);
	for(int i=1;i<=m;i++) if(dis(k,a[i])<=minn) ans++;
	return ans;
}
void sa(){
	for(double T=3000;T>1e-16;T*=0.99){
		double cx=ans.x+(rand()*2-RAND_MAX)*T,cy=ans.y+(rand()*2-RAND_MAX)*T;
		double da=calc({cx,cy})-ansn;
		if(da>0){
			ans.x=cx,ans.y=cy;
			ansn=calc({cx,cy});
		}
		else if(exp(da/T)>rand()) ans={cx,cy};
	}
}
int main(){
	srand(19260817);
	cin>>n>>m>>r;
	for(int i=1;i<=n;i++) cin>>b[i].x>>b[i].y>>b[i].r;
	for(int i=1;i<=m;i++){
		cin>>a[i].x>>a[i].y;
		ans.x+=a[i].x,ans.y+=a[i].y;
	}
	ans.x/=m,ans.y/=m,ansn=calc(ans);
	for(int i=1;i<=100;i++) sa();
	cout<<ansn;
	return 0;
}

P3878 分金币

以这道题为例,我们介绍一种模拟退火的其它应用。

本题大意是把\(n\)个数分成两半,让它们的和的差值最小。

模拟退火可以解决这种序列问题,以随机概率交换两个数,计算当前的差值,如果更优,就接受,否则以一定概率接受(如果不接受就把它们换回去)

AC代码

#include <bits/stdc++.h>
using namespace std;
const double q=0.999,eps=1e-14;
int n,a[60],ans=1e9;
int calc(){
	int s1=0,s2=0;
	for(int i=1;i<=(n+1)/2;i++) s1+=a[i];
	for(int i=(n+1)/2+1;i<=n;i++) s2+=a[i];
	return abs(s1-s2);
}
void sa(){
    for(double T=5000;T>eps;T*=q){
		int l=rand()%n+1,r=rand()%n+1;//取数列中随机两个数
		swap(a[l],a[r]);//交换它们
		int da=ans-calc();
		if(da>0) ans-=da;
		else if(exp(da/T)*RAND_MAX<rand()) swap(a[l],a[r]);
		//这个意思是,如果不满足那个概率,就给他回复原状
	}
}
int main(){
	srand(19260817);
	int t;
	cin>>t;
	while(t--){
		ans=1e9;
		cin>>n;
		for(int i=1;i<=n;i++) cin>>a[i];
		for(int i=1;i<=20;i++) sa();
		cout<<ans<<endl;
	}
}

练习题

  1. P4703 偷上网
  2. P2538 城堡
  3. P4035 球形空间产生器

标签:rand,RAND,int,double,笔记,算法,模拟退火,ans
From: https://www.cnblogs.com/luhaoren/p/SA.html

相关文章

  • 算法笔记(4)莫队算法入门
    原发表于我的博客前言本来想学完回滚莫队、树上莫队、二离莫队之后一起写一个博客,但是一直学不会/kk,只好把已会的普通莫队和带修莫队写了(以后会补上的)普通莫队莫队——优雅的暴力莫队算法的引入例题:给定一个数列和若干询问,每次询问询问一段区间内不同种类数字的个数。暴力......
  • 算法笔记(5)贪心算法
    原发表于我的博客贪心算法贪心与其说是一种算法,不如说一种思想。贪心思想,顾名思义,就是总是做出当前最好的选择,这种方式可能在全局上不是最好的结果,但是在一些题目中就可以直接用。最简单的例子就是“货比三家”,在生活中,我们买东西时都会挑性价比最优的,这就是一种贪心。贪心算......
  • 算法笔记(6)数列分块
    原发表于我的博客前言分块不能说是一种数据结构,它是一种思想,无论是数列分块,块状链表,还是数论分块,莫队算法,都应用了分块的思想。本文主要介绍狭义上的分块,即数列分块。数列分块的引入数列分块可以说是暴力,一种优美的暴力,它的基本思路是,把数列分成若干块(一般取\(\sqrtn\)),分块......
  • 算法笔记(2)FHQtreap
    原发布于我的个人博客前言FHQtreap绝对是平衡树里最好写,最实用的,他几乎能做所有splay或其它平衡树能做的事,还能可持久化!这篇文章将会介绍FHQtreap的基本操作和维护区间的操作,并附上例题。基本操作FHQtreap的基本操作只有两个,分裂和合并。有些读者可能会问,分裂和合并和平衡树......
  • K-medoids聚类算法
    发展:们每次选簇的平均值作为新的中心,迭代直到簇中对象分布不再变化。因此一个具有很大极端值的对象会扭曲数据分布,造成算法对极端值敏感在聚类分析中,异常值通常会引起问题,因为它们可能会被分配到一个独立的聚类,从而干扰正常的聚类结果。这可能导致聚类算法产生不合理或不稳定的......
  • 算法笔记(1)线段树
    原发表于个人博客。前言线段树,是数据结构皇冠上的明珠(我编的)。它用途广泛,被一代代的oier应用,改进,优化。本文介绍了线段树的基础知识和各种拓展(包括权值线段树,可持久化线段树),各种优化方式(包括zkw线段树,动态开点,离散化),希望能帮到更多的oier。在学习线段树前,默认你应该学会一下......
  • 【学习笔记】莫比乌斯反演
    前言/声明首先,本人的数论水平极低,目前莫反只是刚刚入门的水平,此博客的主要作用是用于记录本人的学习过程,真的想要深入了解莫反的话这边推荐cmd大佬的博客(点这里),应该对你有更大帮助。建议学习的时候能多理解些就多去理解,少硬记些结论,这样更不容易忘记。前置知识:最基础的数论。......
  • 【学习笔记】线段树合并
    前置知识:动态开点权值线段树。线段树合并,顾名思义,就是将两棵权值线段树合并在一起。为什么不把两棵普通的线段树合并呢?因为那样好像没啥用。我们知道,权值线段树支持着查询某个数的个数、查询第\(k\)大/小的数等操作,有了合并操作之后就可能会支持一些令人意想不到的操作。放张......
  • 【学习笔记】数论——同余相关
    0前言闲的没事的时候可能会摸鱼写一写,都是些非常基础的东西。最高大概会写到exCRT和exBSGS吧,阶和原根往后的我也不会了,但是前面的内容会时不时来补充。为了方便偷懒,许多定理不会给出证明。1基本概念\(\gcd(a,b)\)或者\((a,b)\):\(a,b\)的最大公约数。\(\rm{lcm}......
  • 算法-共识算法
    一、Paxos    基础的Paxos算法包括如下三种:BasicPaxos、MultiPaxos、FastPaxos     Paxos将系统中的角色分为提议者(Proposer),决策者(Acceptor),和最终决策学习者(Learner):    【Proposer】:提出提案(Proposal)。Proposal信息包括提案编号(ProposalID)......