模拟退火解析

原文链接

前言

模拟退火 S i m u l a t e d Simulated Simulated A n n e a l i n g Annealing Annealing , 简称 S A SA SA ,最早在 1953 1953 1953 年由 N . M e t r o p o l i s N. Metropolis N.Metropolis 提出,后经优化得到现在广泛应用的算法,应用在很多领域当中。

本文题目链接

算法思想

模拟退火是随机化搜索的一种,若随机化搜索写得好,则可以实现高效率和答案的正确率高(虽说不是 100 % 100\% 100% )。很多时候在想不出解决办法,或方法的时间复杂度出现极大情况时,可使用模拟退火。所说是有较大几率正确,但还是有疏漏,那么可以多次试验来更加接近准确地求出这个值(还是要看运气)。

模拟退火,顾名思义,是模拟工业上固体降温的过程。先将固体加温到一定的温度后,在按照适当的温度进行冷却,冷却到改物体想要达到的状态。温度降低地越慢,则该物体的质量约高,因为分子在因热加速运动中找到了更加合适的位置。当温度逐渐降低,分子运动减缓,达成目的。

那么这一现象被科学家与计算机算法所联系起来,就成了现在的模拟退火。

网上的一张图,模拟退火可视化:
在这里插入图片描述

模拟退火有几个很关键的参数,这几个参数决定了模拟退火的优劣。

  1. 随机种子 s e e d seed seed ,可以使用 19260817 19260817 19260817 ,或是时间,不推荐使用其他参数,很可能会降低正确率。
  2. 初始温度 T e m p H i g h TempHigh TempHigh ,一般取 100 100 100 10000 10000 10000 不等,但作者更加倾向于 2000 2000 2000 3000 3000 3000 的数字。
  3. 目标温度 T e m p L o w TempLow TempLow ,一般取 1 − 10 1^{-10} 110 1 − 15 1^{-15} 115
  4. 温度变化率 T e m p L e s s TempLess TempLess ,一般取 0.99 0.99 0.99 0.9999 0.9999 0.9999 。建议不取太大,效率不高。
#define Seed 19260817
#define TempLess 0.9975
#define TempHigh 2879.0
#define TempLow 1e-12

来看看降火的主体部分。

void SA() {
	double temp = TempHigh;//初始化温度
	定义初始状态;
	while(temp > TempLow) {//打到降温条件
		double nowans = Get_Ans(当前状态);//更新最优解
		double diff = nowans - ans;//与当前答案的差值
		if(diff > 0) {//比当前答案更优
			转移状态;
			ans = nowans;//更新答案
		}
		else if(exp(-diff / temp) * RAND_MAX < rand()) {
//接受这个解,为什么这样写请见例题部分
			转移状态;
		}
		temp *= TempLess;//降温
	}
}

模拟退火查询的是多峰函数的最值。

以下曲线是解析式为 y = 0.05 x 3 − 0.5 x 2 y=0.05x^3-0.5x^2 y=0.05x30.5x2 的函数的图像:
在这里插入图片描述
先来考虑贪心的做法:
在这里插入图片描述
当到达点 A A A 时,程序会选择更高的一个点,那么会从 A A A 点到达 B B B 点,而再从 B B B 点俯瞰,看到了点 C C C ,由于 C C C 的纵坐标比 B B B 小,所以点 B B B 不会到达点 C C C 。换句话说,该程序 100 % 100\% 100% 不会接受点 C C C ,进而更不会到达点 D D D 。不难发现,这时候找到的局部最优解并不是全局最优解。

而模拟退火再次做出了改进。假设初始位置在点 A A A ,则会基于 A A A 做出左右摆动,经过数次摆动后到达 B B B 。再进一步摆动,假设摆动到了 C C C 点,但是 C C C 的纵坐标比 B B B 小,会以一定几率以 C C C 的纵坐标来接受 C C C 。进而在以同样的方式摆动到点 D D D ,找到更高点。

由于是该算法随机性较高,所以多跑几遍该函数。

下面结合一道例题更加深入地探究,题目链接已在上文给出。

题意

n n n 个圆, m m m 个点,请求出一个半径不超过 r r r 的圆,使得与这 n n n 个圆没有交集,且能够覆盖的点最大。

思路

此题的答案圆的圆心并不满足是整数,且由横纵坐标两个值来影响,并不具有规律。这样的问题通常使用模拟退火来解决。

void SA() {
	double temp = TempHigh;//初始化温度
	定义初始状态;
	while(temp > TempLow) {//打到降温条件
		double nowans = Get_Ans(当前状态);//更新最优解
		double diff = nowans - ans;//与当前答案的差值
		if(diff > 0) {//比当前答案更优
			转移状态;
			ans = nowans;//更新答案
		}
		else if(exp(-diff / temp) * RAND_MAX < rand()) {//接受这个解
			转移状态;
		}
		temp *= TempLess;//降温
	}
}

如上,初始状态包含了横坐标和纵坐标,为了提高正确率与效率,设为所有点的横纵坐标的平均值。

G e t A n s GetAns GetAns 函数也很简单,先确定半径,半径为这个点到各个圆的切线的距离的最小值,即两个圆心的距离减去当前枚举到的这个圆的半径。后枚举每个点,若这个点被覆盖则 r e s + + res++ res++ ,最后返回 r e s res res

double Get_Ans(double x, double y) {
	double res = 0;
	double rkill = r;
	for(int i = 1; i <= n; i++)//枚举圆
		rkill = Min(rkill, Dist_Cartesian(XC(i), YC(i), x, y) - RC(i));
	for(int i = 1; i <= m; i++)//枚举点
		if(Dist_Cartesian(XE(i), YE(i), x, y) <= rkill)
			res += 1.0;
	return res;
}

有了 G e t A n s GetAns GetAns 函数,主题部分也很快能出来。

void SA() {
	double temp = TempHigh, ansx = initx, ansy = inity;//降温前初始化
	while(temp > TempLow) {
		double nowx = ansx + ((rand() << 1) - RAND_MAX) * temp;
		double nowy = ansy + ((rand() << 1) - RAND_MAX) * temp;
		double nowans = Get_Ans(nowx, nowy);
		double diff = nowans - ans;
		if(diff > 0) {
			initx = nowx;
			inity = nowy;
			ansx = nowx;
			ansy = nowy;
			ans = nowans;
		}
		else if(exp(-diff / temp) * RAND_MAX < rand()) {
			ansx = nowx;
			ansy = nowy;
		}
		temp *= TempLess;
	}
}

首先来看这段代码

double nowx = ansx + ((rand() << 1) - RAND_MAX) * temp;
double nowy = ansy + ((rand() << 1) - RAND_MAX) * temp;

由答案左右摆动,生成新的当前状态 n o w x nowx nowx n o w y nowy nowy ,摆动幅度是随机的,应该是由分子做无规则运动而来。乘上 t e m p temp temp 当前温度是由分子在越热的环境中,运动得越快而得来。

紧接着两行就是求出当前状态的答案,在求出它与当前最优解的差值。

第一个 i f if if 是当前这个局部解大于当前最优解,则用当前最优的局部解来更新最优解。

重点是下一个 i f if if ,这行代码就是它与贪心的不同,以一定几率接受这个解,在用它更新当前状态,进行左右摆动,从而找到局部更优解,更加接近整体最优解。其条件的优越性由 M e t r o p o l i s Metropolis Metropolis 接受准则给出。也就是 e l s e else else i f if if 中的条件:

exp(-diff / temp) * RAND_MAX < rand()

思路整理完了,此题并没有多少思维难度,但是需要对上述几个参数进行调整,可以多总结一些正确率大的参数,以备下次使用

C++代码

#include <cmath>
#include <cstdio>
#include <cstdlib>
#define Min(a, b) ((a) < (b) ? (a) : (b)) 
#define Seed 19260817//随机种子
#define TempLess 0.9975//温度变化率
#define TempHigh 2879.0//初始温度
#define TempLow 1e-12//目标温度
void Quick_Read(double &N) {//double快速读入
	N = 0.0;
	double now, wei = 0.1;
	bool op = false;
	char c = getchar();
	while(c < '0' || c > '9') {
		if(c == '-')
			op = -1;
		c = getchar();
	}
	while(c >= '0' && c <= '9') {
		N = N * 10.0 + (c ^ 48) * 1.0;
		c = getchar();
	}
	if(c == '.') {
		c = getchar();
		while(c >= '0' && c <= '9') {
			N += (c ^ 48) * wei;
			wei /= 10.0;
			c = getchar();
		}
	}
	if(op)
		N = -N;
}
const int MAXN = 15;
const int MAXM = 1e3 + 5;
struct Circle {//题目中的圆
	double Abscissa_C, Ordinate_C, Radius_C;
	#define XC(x) buildings[x].Abscissa_C
	#define YC(x) buildings[x].Ordinate_C
	#define RC(x) buildings[x].Radius_C
};
Circle buildings[MAXN];
struct Enemy {//题目中的点
	double Abscissa_E, Ordinate_E;
	#define XE(x) foe[x].Abscissa_E
	#define YE(x) foe[x].Ordinate_E
};
Enemy foe[MAXM];
double n, m, r;
double initx, inity;//记录答案的横纵坐标
double ans;//答案
double Dist_Cartesian(double XA, double YA, double XB, double YB) {//两点间距离公式
	double frontx = (XA - XB) * (XA - XB);
	double fronty = (YA - YB) * (YA - YB);
	double dist = sqrt(frontx + fronty);
	return dist;
}
double Get_Ans(double x, double y) {//找到当前状态的答案
	double res = 0;
	double rkill = r;
	for(int i = 1; i <= n; i++)//求出最大半径
		rkill = Min(rkill, Dist_Cartesian(XC(i), YC(i), x, y) - RC(i));
	for(int i = 1; i <= m; i++)//求出被圆覆盖的点
		if(Dist_Cartesian(XE(i), YE(i), x, y) <= rkill)
			res += 1.0;
	return res;
}
void SA() {
	double temp = TempHigh, ansx = initx, ansy = inity;//初始化
	while(temp > TempLow) {//降温
		double nowx = ansx + ((rand() << 1) - RAND_MAX) * temp;//当前状态x
		double nowy = ansy + ((rand() << 1) - RAND_MAX) * temp;//当前状态y
		double nowans = Get_Ans(nowx, nowy);//当前局部答案
		double diff = nowans - ans;//当前答案与最优解的差值
		if(diff > 0) {//比当前最优解更优则更新最优解
			initx = nowx;
			inity = nowy;
			ansx = nowx;
			ansy = nowy;
			ans = nowans;
		}
		else if(exp(-diff / temp) * RAND_MAX < rand()) {
//按照Metropolis接受准则接受改状态
			ansx = nowx;
			ansy = nowy;
		}
		temp *= TempLess;//降温
	}
}
void Cool_Down() {
	int frequ = 6;
	while(frequ--)//随机化算法尽量多跑几次
		SA();
}
void Make_Seed() {//生成随机种子
    srand(Seed);
}
void Read() {//输入
    Quick_Read(n);
    Quick_Read(m);
    Quick_Read(r);
    for(int i = 1; i <= n; i++) {
    	Quick_Read(XC(i));
    	Quick_Read(YC(i));
    	Quick_Read(RC(i));
	}
	for(int i = 1; i <= m; i++) {
		Quick_Read(XE(i));
		Quick_Read(YE(i));
		initx += XE(i);
		inity += YE(i);
	}
	initx /= m;//以平均值开始提高效率与准确率
	inity /= m;
}
void Write() {//输出
	printf("%.0lf", ans);
}
int main() {
	Make_Seed();
	Read();
	Cool_Down();
	Write();
	return 0;
}
  • 2
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值