模拟退火算法解析

前言

模拟退火 \(Simulated\) \(Annealing\) , 简称 \(SA\) ,最早在 \(1953\) 年由 \(N. Metropolis\) 提出,后经优化得到现在广泛应用的算法,应用在很多领域当中。

本文题目链接

算法思想

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

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

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

网上的一张图,模拟退火可视化:
模拟退火算法解析

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

  1. 随机种子 \(seed\) ,可以使用 \(19260817\) ,或是时间,不推荐使用其他参数,很可能会降低正确率。
  2. 初始温度 \(TempHigh\) ,一般取 \(100\)\(10000\) 不等,但作者更加倾向于 \(2000\)\(3000\) 的数字。
  3. 目标温度 \(TempLow\) ,一般取 \(1^{-10}\)\(1^{-15}\)
  4. 温度变化率 \(TempLess\) ,一般取 \(0.99\)\(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.05x^3-0.5x^2\) 的函数的图像:
模拟退火算法解析
先来考虑贪心的做法:
模拟退火算法解析
当到达点 \(A\) 时,程序会选择更高的一个点,那么会从 \(A\) 点到达 \(B\) 点,而再从 \(B\) 点俯瞰,看到了点 \(C\) ,由于 \(C\) 的纵坐标比 \(B\) 小,所以点 \(B\) 不会到达点 \(C\) 。换句话说,该程序 \(100\%\) 不会接受点 \(C\) ,进而更不会到达点 \(D\) 。不难发现,这时候找到的局部最优解并不是全局最优解。

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

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

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

题意

\(n\) 个圆, \(m\) 个点,请求出一个半径不超过 \(r\) 的圆,使得与这 \(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;//降温
	}
}

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

\(GetAns\) 函数也很简单,先确定半径,半径为这个点到各个圆的切线的距离的最小值,即两个圆心的距离减去当前枚举到的这个圆的半径。后枚举每个点,若这个点被覆盖则 \(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;
}

有了 \(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;

由答案左右摆动,生成新的当前状态 \(nowx\)\(nowy\) ,摆动幅度是随机的,应该是由分子做无规则运动而来。乘上 \(temp\) 当前温度是由分子在越热的环境中,运动得越快而得来。

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

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

重点是下一个 \(if\) ,这行代码就是它与贪心的不同,以一定几率接受这个解,在用它更新当前状态,进行左右摆动,从而找到局部更优解,更加接近整体最优解。其条件的优越性由 \(Metropolis\) 接受准则给出。也就是 \(else\) \(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;
}
发表评论

评论已关闭。

相关文章