模拟退火的进一步学习（P5544 [JSOI2016]炸弹攻击1）

P5544 [JSOI2016]炸弹攻击1

题目描述

给定平面上$n$个圆与$m$个点，找到平面上某个点，选择某个不大于$R$的长度，满足以这个点为圆心、这个长度为半径的圆与给定圆不相交（可以相切）且覆盖的（$m$个点里的）点数最大。输出这个最大值.

数据范围：$n\le 10$，$m\le 1000$，坐标绝对值和半径$\le 10^3$.

题目分析

我们找不到这道题有什么显然的性质或者推导式，而搜索大概率复杂度爆表，所以采用模拟退火算法.

首先，模拟退火算法一定要有一个初始状态，这里我们假设初始状态是圆在恰好是各个敌人的坐标和的平均值的位置.

然后就是用模拟退火瞎搞一气了.

主要注意坐标在更新时，新状态改变的值可以这样来设：

		double sx=(rand()*2-RAND_MAX)*T,sy=(rand()*2-RAND_MAX)*T;//随机变化坐标，变化幅度为 T
		//rand()*2的范围为0到65534（所有偶数），减去RAND_MAX（32767）后范围为-32767到32767（所有奇数），实现了正负覆盖

这算是一个套路了，很多涉及坐标的模拟退火题目（比如P1337）都是这样子表示状态的转移的.

---

另外，这里我还想简单说一下上次没说清楚的这一个语句：

		else if(exp(-delta/T)*RAND_MAX<rand())

模拟退火的这个地方，温度越低，与最优解的差值越大，越难被接纳.

也就是说，这一句的不等号是大于号还是小于号（当然也可以是$delta$符号的改变，为了统一，这里都讨论不等号），是与前一个$if$语句$delta$是正好还是负好，是有关系的.

比如说，这是本题的$SA$.

void SA(){
	double T=1200;
	while(T>esp){
		double sx=(rand()*2-RAND_MAX)*T,sy=(rand()*2-RAND_MAX)*T;
		int now=deal(sumx+sx,sumy+sy);
		int delta=now-ans1;//delta统一为新状态减去后状态，显然在这道题中，delta>0比较优
		if(delta>0){
			sumx+=sx;
			sumy+=sy;
			ans1=now;
			ans2=max(ans1,ans2);
		}
		else if(exp(-delta/T)*RAND_MAX<rand()){//那么这里剩下的情况就是delta<0了，明显delta越小，越难被选到，而delta符号统一为-，所以用小于号
			sumx+=sx;
			sumy+=sy;
			ans1=now;
		}
		T*=dt;
	} 
}
//rand()/RAND_MAX可以看做一个固定的概率

而在另外一道题P3878分金币中，我们的$delta$却是越小越好.

那道题的$SA$如下：

void SA(){
	double T=1200;
	while(T>esp){
		int x=(int )(rand()%len1+1),y=(int )(rand()%len2+len1+1);
		ll delta=abs((tot1+a[y]-a[x])-(tot2+a[x]-a[y]))-ans;//随机新状态
		if(delta<0){
			tot1=tot1+a[y]-a[x];
			tot2=tot2+a[x]-a[y];
			ans=abs(tot1-tot2);
			swap(a[x],a[y]);
		}
		else if(exp((-delta)/T)*RAND_MAX>rand()){//RAND_MAX表示rand函数能够取得到的最大值，这是对是否更新状态的估值；exp表示自然对数的底的k次方
			tot1=tot1+a[y]-a[x];
			tot2=tot2+a[x]-a[y];
			swap(a[x],a[y]);//改变的是状态，不是用来记录最优解的ans 
		}
		T*=dt;
	}
}

由于这道题中是delta<0肯定能取到，delta>0的话越大越不好，也就是-delta越小越不好，所以符号是大于号.

以上，就是关于$delta$概率语句的详细解释. 分金币的代码引自我的另外一篇博客模拟退火学习报告，有兴趣的同学可以看一下.

程序实现

#include<bits/stdc++.h>
#define esp 1e-10
#define dt 0.998
using namespace std;
struct building{
	int x,y,r;
}bld[20];
struct enemy{
	int x,y;
}enm[1010];
int n,m,r,ans1,ans2;
double sumx,sumy;
double dist(double x1,double y1,double x2,double y2){
	return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}
int deal(double x,double y){
	int sum=0;
	double R=(double)r;
	for(int i=1;i<=n;i++){
		double dis=dist(x,y,bld[i].x ,bld[i].y )-(double)bld[i].r ;
		R=min(R,dis);
	}
	for(int i=1;i<=m;i++){
		double dis=dist(x,y,enm[i].x ,enm[i].y );
		if(dis<=R)sum++;
	}
	return sum;
}//O(n+m)求被杀死的敌人数量
void SA(){
	double T=1200;
	while(T>esp){
		double sx=(rand()*2-RAND_MAX)*T,sy=(rand()*2-RAND_MAX)*T;//随机更新新状态
		int now=deal(sumx+sx,sumy+sy);
		int delta=now-ans1;
		if(delta>0){
			sumx+=sx;
			sumy+=sy;
			ans1=now;
			ans2=max(ans1,ans2);//ans2专门记录最优解，不是状态，不用改变
		}
		else if(exp(-delta/T)*RAND_MAX<rand()){
			sumx+=sx;
			sumy+=sy;
			ans1=now;//随机改变状态
		}
		T*=dt;
	} 
}
int main(){
	srand((int )time(0));//随机时间种子
	scanf("%d%d%d",&n,&m,&r);
	for(int i=1;i<=n;i++){
		scanf("%d%d%d",&bld[i].x ,&bld[i].y ,&bld[i].r );
	}
	for(int i=1;i<=m;i++){
		scanf("%d%d",&enm[i].x ,&enm[i].y );
		sumx+=enm[i].x ,sumy+=enm[i].y ;
	}
	sumx/=m,sumy/=m;//确定初始位置
	for(int i=1;i<=6;i++)SA();//为保证准确性，多次模拟退火
	printf("%d\n",ans2);
	return 0;
}