模擬退火算法介紹(摘自08集訓隊 顧研《淺談隨機化思想在幾何問題中的應用》)

    一.模擬退火算法的原理

模擬退火算法是一種元啟發式（Meta-Heuristics）算法，來源于固體退火原理，將固體加熱至充分高的溫度，再讓其徐徐冷卻。加熱時，固體內部粒子隨溫升變為無序狀，內能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達到平衡態，最后在常溫時達到基態，內能減為最小。根據Metropolis準則，粒子在溫度T時趨于平衡的概率為e-ΔE/kt，其中E為溫度T時的內能，ΔE為其改變量，k為Boltzmann常數。

    二.模擬退火算法的模型

① 初始化：初始溫度T(充分大)，初始解狀態S(算法迭代的起點)， 每次迭代次數L

② for k=1 to L 做③至⑥

③ 產生新解S’

④ 計算增量Δt′=C(S′)-C(S)，其中C(S)為評價函數

⑤ 若Δt′<0則接受S’作為新的當前解，否則以概率e-Δt/k接受S’作為新的當前解

⑥ 如果滿足終止條件則輸出當前解作為最優解，結束程序

⑦ T逐漸減少，然后轉②

    回到此題，題意為：平面上有一個矩形，在矩形內有一些陷阱。求得矩形內一個點，該點離與它最近的已知陷阱最遠（點的個數≤1000）。精度要求：1/10。

這題的精度要求不高，一個很樸素的想法是枚舉平面中的一些網格點并進行判斷。當然這樣的點太多了，我們必須選一些比較有“希望”的點，同時還不能忽略任何平面上的任何一點。

我們使用類比的方法引入模擬退火算法，初始解狀態S可以在矩形內隨便選取，初始溫度T對應于每次調整的距離D，產生新解的方式是在目前解為圓心、半徑為D的圓周上任取一點，評價函數取距離最近的陷阱的距離，終止條件為D足夠小。

由此可得本問題的模擬退火算法：由初始解S和控制參數初值D開始，對當前解重復“產生新解→計算目標函數差→接受或舍棄”的迭代，并逐步衰減D值，算法終止時可以得到一組解，這是基于蒙特卡羅迭代求解法的一種啟發式隨機搜索過程。退火過程由冷卻進度表（Cooling Schedule）控制，包括控制參數的初值D及其衰減因子D’、每個值D時的迭代次數L。

模擬退火算法還有一個特點：具有并行性。因此我們可以將初始解S改成初始解集，對于每個候選解都進行迭代，答案取最終解集的最優解。

由于我們必須保證能覆蓋矩形上的每個位置，因此在①中確定p后（我們按網格狀放置），②中的delta可取max{矩形邊長}/sqrt(n)。設算法的迭代次數為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;
}