模拟退火学习笔记

虽然说考前不应该碰这些随机化算法，容易影响思考，但是还是写一写吧，对于一些问题还是很好用的。

概念

什么是模拟退火。一句话解释，我们从一个旧状态随机出一个新状态，要从旧状态转移到新状态，如果新状态小于旧状态那么就接受，否则以一定概率接受这个新状态。

注意这里的接受并不是指更新答案，是指保留当前状态。

感觉这个解释很好理解啊！我们直接来看例题吧。

应用——计算几何

P1337 [JSOI2004] 平衡点 / 吊打XXX

这道题题面描述写的比较阴间，大概意思就是给定你 \(n\) 个点，每个点有一个质量 \(m_i\)，然后你要找一个点 \(x\) 使得 \(\sum_{i=1}^n m_i \times \text{dist}_{i,x}\) 最小。 其中 \(\text{dis}_{i,x}\) 表示 \(i\) 和 \(x\) 之间的距离。 这你就可以直接退火了。

模拟退火你要考虑清楚两个问题，第一个是初始状态，第二个是计算代价。

那么这道题的初始状态是什么？其实初始状态你想设啥设啥，但是你为了让答案更加准确，需要设定一个尽量优的解。所以我们令初始状态是给定的坐标的平均数。

至于这个计算代价，这道题直接模拟一下就好了。

提供一个我觉得还不错的模板，后面变形要用。

#include<bits/stdc++.h>
using namespace std;
#define temp 1e5
#define cold 0.996 
const int N=1e5+10;
int n;
double ans,ansx,ansy;
struct node{
	int x,y,w;
}a[N];
double cost(double x,double y)//计算代价
{
	double energy=0.0;
	for(int i=1;i<=n;i++)
	{
		energy+=sqrt((x-a[i].x)*(x-a[i].x)+(y-a[i].y)*(y-a[i].y))*a[i].w;
	}
	return energy;
}
void SA()
{
	double t=temp;//初始温度
	while(t>1e-18)//退火
	{
		double tmpx=ansx+(rand()+rand()-RAND_MAX)*t,tmpy=ansy+(rand()+rand()-RAND_MAX)*t;
		double tmp=cost(tmpx,tmpy);
		double d=tmp-ans;
		if(d<0.0)ans=tmp,ansx=tmpx,ansy=tmpy;//如果更优直接接受
		else if(exp(-d/t)*RAND_MAX>rand())ansx=tmpx,ansy=tmpy;//否则以一定概率接受
		t*=cold;//降温
	}
}
void solve()//多次退火增加正确性
{
	ans=cost(ansx,ansy);
	SA();
	SA();
	SA();
	SA();
	SA();
	SA();
	
}
int main()
{
	srand(time(0));
	cin>>n;
	for(int i=1;i<=n;i++)
	{
		cin>>a[i].x;
		cin>>a[i].y;
		cin>>a[i].w;
		ansx+=a[i].x;
		ansy+=a[i].y;
		
	}
	ansx/=n;
	ansy/=n;//初始状态
	solve();
	printf("%.3lf %.3lf\n",ansx,ansy);
	return 0;
}

应用——方差/平均数

P2503 [HAOI2006]均分数据

不知道为什么，但这种分组题真的很适合退火。

首先老套路，初始状态，这题你可以一开先跑一个假贪心，把这个数组扫一变，每个数放到当前和最小的那个组里。然后如何计算代价，实际上对于每一个状态跑一遍假贪心就可以了。

注意这道题跟上面那题的区别，上面那个题状态是一个坐标，那么这个题的状态是什么呢？我们考虑这个题暴力是怎么写的，是不是枚举全排列然后算。那么模拟退火时我们可以类似这样，每次随机出来两个下标，然后交换这两个下标。

#include<bits/stdc++.h>
using namespace std;
#define cold 0.95
#define temp 3000
const int N=50;
int n;
int m;
int a[N];
int tmp[N];
double ans=1e9;
double sum;
double calc()//假贪心计算代价
{
	int w,mn;
	double summ=0;
	memset(tmp,0,sizeof(tmp));
	for(int i=1;i<=n;i++)
	{
		mn=1e9;
		for(int j=1;j<=m;j++)
		{
			if(tmp[j]<mn)
			{
				mn=tmp[j];
				w=j;
			}
		}
		tmp[w]+=a[i];
	}
	for(int i=1;i<=m;i++)
	{
		summ+=(sum-tmp[i])*(sum-tmp[i]);
	}
	return sqrt(summ/m);
}
void SA()
{
	double t=temp;
	while(t>1e-16)
	{
		int x=rand()%n+1;
		int y=rand()%n+1;
		swap(a[x],a[y]);//随机两个下标交换
		double cnt=calc();//初始状态，这个其实可以放到下面去的
		
		if(cnt<ans)ans=cnt;//接受
		else if(exp((cnt-ans)/t)<double(rand())/RAND_MAX)swap(a[x],a[y]);//随机接受
		t*=cold;
	}
}
int main()
{
	srand(time(0));
	cin>>n>>m;
	for(int i=1;i<=n;i++)
	{
		cin>>a[i];
		sum+=a[i];
	}
	sum/=m;
	while(clock()<CLOCKS_PER_SEC*0.9) SA();//算是模拟退火中一个比较重要的应用，注意有些题卡不了时就别卡了
	printf("%.2f",ans);
	return 0;
}

P3878 [TJOI2010]分金币

这道题又是一个分组的题，那模拟退火去吧！

老套路，初始状态是什么，我们直接对原数组对半砍，作为初始代价。转移状态？跟上面那题一样，每次随机两个下标交换。计算代价？按照题目模拟即可。

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define cold 0.981
#define temp 3000
const int N=50;
int n;

int a[N];
int ans;
int calc()//计算代价
{
	int mid=(1+n)>>1;
	int sum1=0;
	int sum2=0;
	for(int i=1;i<=mid;i++) sum1+=a[i];
	for(int i=mid+1;i<=n;i++) sum2+=a[i];
	return abs(sum1-sum2);
}
void SA()
{
	double t=temp;
	while(t>1e-16)
	{
		int x=rand()%n+1;
		int y=rand()%n+1;
		swap(a[x],a[y]);//随机两个下标
		int cnt=calc();
		if(cnt<ans)ans=cnt;
		else if(exp((ans-cnt)/t)<double(rand())/RAND_MAX)swap(a[x],a[y]);
		t*=cold;
	}
}
signed main()
{
	srand(time(0));
	int T;
	cin>>T;
	while(T--)
	{
		cin>>n;
		for(int i=1;i<=n;i++) cin>>a[i];
		ans=1e9;
		for(int i=1;i<=100;i++) SA();
		cout<<ans<<endl;
		
	}
	return 0;
}

P2210 Haywire

初始状态是什么？可以对于最开始的这个排列即 \(1 \to n\) 求一下。

如何转移？随机两个下标交换。

计算代价？直接对于每一头牛算一下他和他三个朋友的位置差即可。但这样会重复算两次所以最后答案还要除 \(2\) 。

#include<bits/stdc++.h>
using namespace std;
#define cold 0.996
#define temp 1e5
const int N=50;
int n;

int a[N][N];
int ans;
int pos[N];
int calc()//计算代价
{
	int ret=0;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=3;j++)
		{
			ret+=abs(pos[i]-pos[a[i][j]]);
		}
	}
	return ret;
}
void SA()
{
	double t=temp;
	while(t>1e-16)
	{
		int x=rand()%n+1;
		int y=rand()%n+1;
		swap(pos[x],pos[y]);//交换位置
		int cnt=calc();
		if(cnt<ans)ans=cnt;
		else if(exp((ans-cnt)/t)<double(rand())/RAND_MAX)swap(pos[x],pos[y]);
		t*=cold;
	}
}
void solve()
{
	SA();
	SA();
	SA();
	SA();
	SA();
	SA();
	SA();
	cout<<ans/2<<endl;//答案要除2
}
int main()
{
	srand(time(0));
	cin>>n;
	for(int i=1;i<=n;i++)
	{
		cin>>a[i][1];
		cin>>a[i][2];
		cin>>a[i][3];
		pos[i]=i;//记录一下位置
	}
	ans=calc();
	solve();
	return 0;
}

P7962 [NOIP2021] 方差

模拟退火典中典应用。要不是这道题我为什么会来学模拟退火。

预警：这道题的数据不算特别水，退火能拿到较高的分数但想要通过需要一些玄学操作和性质的发现。

这道题首先那道题有一个显然的事实，题目给的这个操作的本质就时交换两个差分数组，因为你考虑 \(a_i=a_{i-1}+a_{i+1}-a{i}\) 就是 \(a_i-a_{i-1}=a_{i+1}-a_i\) 。

那么我们就可以写出一个最基础的退火。

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define cold 0.94
#define temp 1000
const int N=1e5+10;
int n,m;
int sum[N];
int a[N];
int ans=1e9;
int calc()//计算代价
{
	int cnt=0;
	int ret=0;
	for(int i=2;i<=n;i++)
	{
		a[i]=a[i-1]+sum[i];//差分数组的前缀和就是原数组
	}
	for(int i=1;i<=n; i++) 
	{
		cnt+=a[i];
	}
	for(int i=1;i<=n;i++) 
	{
		ret+=(a[i]*n-cnt)*(a[i]*n-cnt);
	}
	return ret/n;
}
void SA()
{
	double t=temp;
	while(t>1e-16)
	{
		int x=rand()%(n-1)+2;
		int y=rand()%(n-1)+2;
		swap(sum[x],sum[y]);//交换差分数组
		int cnt=calc();
		if(cnt<ans)ans=cnt;
		else if(exp((ans-cnt)/t)<double(rand())/RAND_MAX)swap(sum[x],sum[y]);
		t*=cold;
	}
}
int cmp(int a,int b)
{
	return a>b;
}
signed main()
{
	srand(time(0));
	cin>>n;
	for(int i=1;i<=n;i++)
	{
		cin>>a[i];
		sum[i]=a[i]-a[i-1];
	}
	while(clock()<CLOCKS_PER_SEC*0.95) SA();
	cout<<ans<<endl;
	return 0;
}

然后你就需要稍微厉害一点的找性质能力了，我们模拟一下样例发现差分数组时单谷的，那我们大胆对于所有的最优解答案都是单谷的。具体如何证明我也不大会，题解里说了。

所以我们一开始先将差分数组对半砍，前半段从大到小排序，后半段从小到大排序构成一个单谷函数，这样找答案的时候会快很多。

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define cold 0.94
#define temp 1000
const int N=1e5+10;
int n,m;
int sum[N];
int a[N];
int ans=1e9;
int calc()
{
	int cnt=0;
	int ret=0;
	for(int i=2;i<=n;i++)
	{
		a[i]=a[i-1]+sum[i];
	}
	for(int i=1;i<=n; i++) 
	{
		cnt+=a[i];
	}
	for(int i=1;i<=n;i++) 
	{
		ret+=(a[i]*n-cnt)*(a[i]*n-cnt);
	}
	return ret/n;
}
void SA()
{
	double t=temp;
	while(t>1e-16)
	{
		int x=rand()%(n-1)+2;
		int y=rand()%(n-1)+2;
		swap(sum[x],sum[y]);
		int cnt=calc();
		if(cnt<ans)ans=cnt;
		else if(exp((ans-cnt)/t)<double(rand())/RAND_MAX)swap(sum[x],sum[y]);
		t*=cold;
	}
}
int cmp(int a,int b)
{
	return a>b;
}
signed main()
{
	srand(time(0));
	cin>>n;
	for(int i=1;i<=n;i++)
	{
		cin>>a[i];
		sum[i]=a[i]-a[i-1];
	}
 	sort(sum+2,sum+1+n/2,cmp);
 	sort(sum+n/2+1,sum+n+1);
	while(clock()<CLOCKS_PER_SEC*0.95) SA();
	cout<<ans<<endl;
	return 0;
}