模擬退火算法介紹(摘自08集訓隊 顧研《淺談隨機化思想在幾何問題中的應用》)
一.模擬退火算法的原理
模擬退火算法是一種元啟發(fā)式（Meta-Heuristics）算法，來源于固體退火原理，將固體加熱至充分高的溫度，再讓其徐徐冷卻。加熱時，固體內(nèi)部粒子隨溫升變?yōu)闊o序狀，內(nèi)能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達到平衡態(tài)，最后在常溫時達到基態(tài)，內(nèi)能減為最小。根據(jù)Metropolis準則，粒子在溫度T時趨于平衡的概率為e-ΔE/kt，其中E為溫度T時的內(nèi)能，ΔE為其改變量，k為Boltzmann常數(shù)。
二.模擬退火算法的模型
① 初始化：初始溫度T(充分大)，初始解狀態(tài)S(算法迭代的起點)， 每次迭代次數(shù)L
② for k=1 to L 做③至⑥
③ 產(chǎn)生新解S’
④ 計算增量Δt′=C(S′)-C(S)，其中C(S)為評價函數(shù)
⑤ 若Δt′<0則接受S’作為新的當前解，否則以概率e-Δt/k接受S’作為新的當前解
⑥ 如果滿足終止條件則輸出當前解作為最優(yōu)解，結(jié)束程序
⑦ T逐漸減少，然后轉(zhuǎn)②
回到此題，題意為：平面上有一個矩形，在矩形內(nèi)有一些陷阱。求得矩形內(nèi)一個點，該點離與它最近的已知陷阱最遠（點的個數(shù)≤1000）。精度要求：1/10。
這題的精度要求不高，一個很樸素的想法是枚舉平面中的一些網(wǎng)格點并進行判斷。當然這樣的點太多了，我們必須選一些比較有“希望”的點，同時還不能忽略任何平面上的任何一點。
我們使用類比的方法引入模擬退火算法，初始解狀態(tài)S可以在矩形內(nèi)隨便選取，初始溫度T對應于每次調(diào)整的距離D，產(chǎn)生新解的方式是在目前解為圓心、半徑為D的圓周上任取一點，評價函數(shù)取距離最近的陷阱的距離，終止條件為D足夠小。
由此可得本問題的模擬退火算法：由初始解S和控制參數(shù)初值D開始，對當前解重復“產(chǎn)生新解→計算目標函數(shù)差→接受或舍棄”的迭代，并逐步衰減D值，算法終止時可以得到一組解，這是基于蒙特卡羅迭代求解法的一種啟發(fā)式隨機搜索過程。退火過程由冷卻進度表（Cooling Schedule）控制，包括控制參數(shù)的初值D及其衰減因子D’、每個值D時的迭代次數(shù)L。
模擬退火算法還有一個特點：具有并行性。因此我們可以將初始解S改成初始解集，對于每個候選解都進行迭代，答案取最終解集的最優(yōu)解。
由于我們必須保證能覆蓋矩形上的每個位置，因此在①中確定p后（我們按網(wǎng)格狀放置），②中的delta可取max{矩形邊長}/sqrt(n)。設算法的迭代次數(shù)為t，則算法的復雜度O（P*L*t*n）。
POJ 1379
#include<iostream>
#include<time.h>
#include<cmath>
#define sqr(x) (x)*(x)
#define FOR(i,a,b) for(int i=a;i<b;i++)
#define FF(i,a) for(int i=0;i<a;i++)
using namespace std;
const double eps=1e-3;
const double INF=1e20;
const double pi=acos(-1.0);
struct Houxuan
{
double x,y,dist;
}pp[30];
struct PT
{
double x,y;
}hole[1000];
double dis(double x1,double y1,double x2,double y2)
{
return sqrt(sqr(x1-x2)+sqr(y1-y2));
}
int n,cas;
int X,Y;
int main()
{
int P=15,L=25;
scanf("%d",&cas);
srand((unsigned)time(NULL));
while(cas--)
{
scanf("%d%d%d",&X,&Y,&n);
double delta=double(max(X,Y))/(sqrt(1.0*n));
FF(i,n)
scanf("%lf%lf",&hole[i].x,&hole[i].y);
FF(i,P)
{
pp[i].x=double(rand()%1000+1)/1000.000*X;
pp[i].y=double(rand()%1000+1)/1000.000*Y;
pp[i].dist=INF;
FF(j,n)
pp[i].dist=min(pp[i].dist,dis(pp[i].x,pp[i].y,hole[j].x,hole[j].y));
}
while(delta>eps)
{
FF(i,P)
{
Houxuan tmp=pp[i];
FF(j,L)
{
double theta=double(rand()%1000+1)/1000.000*10*pi;
Houxuan t;
t.x=tmp.x+delta*cos(theta);
t.y=tmp.y+delta*sin(theta);
if(t.x>X||t.x<0||t.y>Y||t.y<0)
continue;
t.dist=INF;
FF(k,n)
{
t.dist=min(t.dist,dis(t.x,t.y,hole[k].x,hole[k].y));
}
if(t.dist>pp[i].dist)
pp[i]=t;
}
}
delta*=0.8;
}
int idx=0;
FOR(i,1,P)
{
if(pp[i].dist>pp[idx].dist)
idx=i;
}
printf("The safest point is (%.1lf, %.1lf).\n",pp[idx].x,pp[idx].y);
}
return 0;
}