模擬退火演算法解析

Emissaryofnight發表於2021-01-01
演算法

前言

模擬退火 \(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;
}

相關文章