原發(fā)表于個人博客=

模擬退火的引入

假如我們有一個函數(shù)，要求它的極大值，怎么求呢？

如果這個函數(shù)滿足單調(diào)性，可以用二分的方法。

如果這是一個單谷（或單峰）函數(shù)，可以用三分法。

那要是多峰函數(shù)怎么半呢？

這時就可以用隨機化算法。

一種樸素的方法是：每次在當前找到的最優(yōu)方案\(x\)附近尋找一個新方案。如果這個新的解\(x'\)更優(yōu)，那么轉(zhuǎn)移到 \(x'\)，否則就不轉(zhuǎn)移。（這被稱作爬山算法）

但是這樣就容易陷入一種局部最優(yōu)解，如圖，如果使用爬山算法，找到的答案可能在一個范圍內(nèi)是最優(yōu)解，但是在全局上并不是最優(yōu)解。

所以就有了模擬退火算法。

模擬退火算法的基本思路是，在一定的概率下接受一個非最優(yōu)解而跳出局部最優(yōu)解。

為什么這個算法被稱為模擬退火呢？因為物理中金屬降溫是隨機的，這個金屬降溫的過程也叫退火。

基本概念

剛剛講到了模擬退火接受非最優(yōu)解的概率，這個概率被稱為\(Metropolis\)準則，也就是

其中，\(\Delta E\)表示最優(yōu)解和當前解的差，\(T\)是當前溫度。

模擬退火一般有三個參數(shù)，初始溫度\(T_0\)，降溫系數(shù)\(\Delta T\)（一般取\(0.99\)左右），最終溫度\(T_{last}\)（一般是個極小值）。

當前溫度\(T\)從\(T_0\)開始，不斷乘以\(\Delta T\)，在這個過程中不斷接受新解，直到\(T\)降至\(T_{last}\)，算法結(jié)束。

模板

這里以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;//計算新解與最優(yōu)解的差
		if(da<0){
		    //如果新解比最優(yōu)解還好，直接接受
			now=anst={cx,cy};
			ans=k;
		}
		else if(exp(-da/t)*RAND_MAX>rand()) now={cx,cy};//否則根據(jù)概率接受
	}
}
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中，模擬退火一般有三種用途，一種是計算幾何（比如費馬點問題），一種是序
列問題（隨機交換），還有一種是騙分。

這篇博客的最后會有幾道例題，并附講解。

技巧

模擬退火本質(zhì)上是一種隨機化算法，能否正確取決于參數(shù)和rp。

這里提供幾種技巧，在考試時可以使用。

卡時技巧

一般來說，模擬退火要跑很多遍才能跑出最優(yōu)解，但是沒遍模擬退火的時間我們也不知道，所以可以用\(clock\)函數(shù)，判斷剩余時間，來卡時。

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

注：clock函數(shù)返回從程序開始到當前的時間，CLOCKS_PER_SEC將clock()函數(shù)的結(jié)果轉(zhuǎn)化為以秒為單位的量，每個系統(tǒng)都不一樣，在windows系統(tǒng)里是1000。

參數(shù)調(diào)小

一般來說，參數(shù)過大會WA，參數(shù)過小會TLE。

但是TLE的可能性一般比WA小，所以調(diào)參時寧小毋大，如平衡點一題，如果參數(shù)設置為\(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是隨機數(shù)的最大值，所以(rand()*2-RAND_MAX)的取值范圍就是[-RAND_MAX,RAND_MAX]，再乘上當前的溫度\(t\)，就會得到一個可能正負的隨機浮點數(shù)，并且隨著溫度的減小越來越小，滿足了我們的需求。

例題

P1337 平衡點 / 吊打XXX

原題

這題需要一定的物理知識。

根據(jù)能量最低原理，一個系統(tǒng)內(nèi)能量越低就越穩(wěn)定，所以題目要求的平衡情況能量肯定最小。

分析這道題，如果靜止不動，則能量只有重力勢能，所以直接令重力勢能最小即可。

我們知道，重力勢能的公式是

\(g\)是個常量，所以直接算每個點的質(zhì)量和高度乘積即可。

簡化題意后，就是尋找一個點，讓這個點到每個點的距離乘以質(zhì)量最小，也就是一個帶權費馬點問題，套用模板即可。

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});//比較新解和最優(yōu)解
		double da=k-ans;//計算差值
		if(da<0){
		    //如果新解比最優(yōu)解好，直接接受
			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\)個數(shù)分成兩半，讓它們的和的差值最小。

模擬退火可以解決這種序列問題，以隨機概率交換兩個數(shù)，計算當前的差值，如果更優(yōu)，就接受，否則以一定概率接受（如果不接受就把它們換回去）

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;//取數(shù)列中隨機兩個數(shù)
		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 偷上網(wǎng)
2. P2538 城堡
3. P4035 球形空間產(chǎn)生器