模擬退火筆記

阿新 • • 發佈：2020-11-02

0-模擬退火的由來

模擬退火演算法（Simulate Anneal，SA）是一種通用概率演演算法，用來在一個大的搜尋空間內找尋命題的最優解。模擬退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所發明的。V.Černý在1985年也獨立發明此演演算法。模擬退火演算法是解決TSP問題的有效方法之一。

模擬退火的出發點是基於物理中固體物質的退火過程與一般組合優化問題之間的相似性。模擬退火演算法是一種通用的優化演算法，其物理退火過程由加溫過程、等溫過程、冷卻過程這三部分組成。

——百度百科

簡單說，模擬退火是一種隨機化演算法。當一個問題的方案數量極大（甚至無窮）而且不是單峰函式時，常使用模擬退火求解

我其實一直覺得模擬退火挺玄學的。~~（後來發現還真就靠RP吃飯~~

Metropolis接受原則

每次隨機出一個新解，如果這個解更優，則接受它，否則以一個與溫度和與當前最優解的差相關的概率接受它。

設這個新的解與最優解的差為 $\Delta E$，溫度為 $T$ ， $k$ 為一個隨機數，那麼這個概率為 $e^{\frac{\Delta E}{kT}}$

具體來說：設 $P$ 為接受當前這個解的概率，那麼

若 $\Delta E >0$ ，那麼 $P=1$ ，也就是必定接受更優解；

若 $\Delta E <0$，那麼 $P=exp( -\frac{\Delta E}{kT} ) $，也就是以一個和差相關的概率接受這個不優的解。

對第二種情況的具體做法將在後面進行詳細解釋。

SA的簡要流程

所需引數：初始溫度 $Tb$，降溫的係數 $D$，最終溫度 $Te$

（注：一般來講 $Tb$會比較大，$D$是一個小於$1$且接近$1$的$double$，比如$0.97$,$Te$是一個接近$0$的正值，

如$1e-10$ ，這些都是程式中自定的常數，但是在實際情況中這些引數需要進行不斷調優，也就是所謂的調參，

這或許也是模擬退火很難應用到$OI$賽制中的原因，畢竟不能實時檢視結果的話就真的很靠$RP$，

當然調參也是有跡可循的，比如$TLE$的時候可以考慮吧初溫調小，$WA$

的時候調大之類的）

流程如下：

令 $T=Tb$
進行轉移
$T*=D$
重複步驟$2、3$ 直到 $T<Te$

轉移（重點qwq）

設估價函式為 $f(x)$ （其實就是對你隨機出來的方案計算答案的一個函式），$x$為原先得到的

解，$x1=x+T0*rand()$ 為新得到的隨機出來的解（ $-1<=R<=1$,且$R!=0$ ），$\Delta E =f(x1)- f(x)$

~~雖然我覺得到這裡還是一頭霧水，但是看到程式碼應該就比較清楚~~

放個虛擬碼

if  ( DelE>0 )  更新解；
else if ( exp(-DelE/T)<(double)rand()/RAND_MAX ) 更新解；

然後發現一個奇怪的東西：if ( exp(-DelE/T)<(double)rand()/RAND_MAX ) 這是幹嘛用的

$exp$其實就是一個$math$函式，得到以$e$為底的自然對數，不用管它；

$-DelE/T$ 就是剛剛所說，$\Delta E <0$的情況，得到的一個概率，後面用rand()確保概率的隨機性，

也就是以一定的概率隨機接受這個解。需要注意的是，這裡的$<$很有講究。

首先，顯然的是 $-DelE$ 是一個正數，然後 $/T$保證了這個概率和當前降溫的進度有關：

當降溫剛剛開始的時候，$T$比較大，得到的概率比較小，容易使小於號成立，從而在開始時更有可能跳

出局部最優解也就是貪心的弊端得到全域性最優解；快要結束的時候反之。然後，當 $-\Delta E$比較

大，也就是和當前最有解相差較大時，接受這個解的概率小，反之亦然。

需要注意的是，這裡針對的是求最大值的情況，如果是最小值需要把小於號改成大於號。原因同理。

最後的最後——一些細節

隨機的初始化真的很重要
調參這個東西有時候很考經驗
其實模擬退火可以解決基本所有最優解問題 ~~(所以快去，用SA全做一遍~~
有時候為了求得最優解，需要多次SA，這裡有個不TLE的好辦法：

 ``` while ( (double)clock()/CLOCKS_PER_SEC<=0.8 ) SA();   ```
 
 0.8根據時限而定，一般比時限小一點點

模板題——luogu P5544 炸彈攻擊1

題目連結

/*
P5544 [JSOI2016]炸彈攻擊1（作為模擬退火的板子qwq
題意：有不超過10個己方建築和1e3個敵方（以圓心半徑形式給出
可以任意放置一個半徑不超過r的炸彈，求不傷害己方能消滅敵方建築的最大值（邊界上算敵不算己
思路：模擬退火一個位置。然後列舉計算最大半徑和個數。
*/
#include <bits/stdc++.h>
using namespace std;
const int N=11,M=1010;
const double delta=0.996,Te=1e-10;			//D和最終溫度
int n,m,x[N],y[N],r[N],p[M],q[M],R,ansout=0;
double ansx,ansy;

int read()
{
    int x=0,f=1; char ch=getchar();
    while ( ch<'0' || ch>'9' ) 
    {
        if ( ch=='-' ) f=-1; ch=getchar();
    }
    while ( ch>='0' && ch<='9' ) x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    return x*f;
}

double dis( double ax,double ay,double bx,double by )		//求兩點間距離
{
    double res=(bx-ax)*(bx-ax)+(by-ay)*(by-ay);
    res=sqrt(res);
    return res;
}

int calc( double xx,double yy )			//估價函式
{
    double nowr=R;
    for ( int i=1; i<=n; i++ )
        nowr=min( nowr,dis(xx,yy,x[i],y[i])-(double)r[i] );
    int cnt=0;
    for ( int i=1; i<=m; i++ )
        if ( dis(xx,yy,p[i],q[i])<=nowr ) cnt++;
    return cnt;
}

void SA()
{
    double T=6000;		//初始溫度
    int ans=0;
    while ( T>Te )
    {
        double nowx=ansx+(rand()*2-RAND_MAX)*T;
        double nowy=ansy+(rand()*2-RAND_MAX)*T;
        int nowans=calc( nowx,nowy );
        int DelE=nowans-ans;
        if ( DelE>0 ) ansx=nowx,ansy=nowy,ans=nowans,ansout=max( ansout,ans );
        //更優的解，必然接受
        else if ( exp(-DelE/T)<(double)rand()/RAND_MAX ) ansx=nowx,ansy=nowy,ans=nowans;
        //按概率接受，保證了越到後期，和最優解的差距越大，越難被接受
        T*=delta;
    }
}

int main()
{
    srand(time(0)); ansx=ansy=0;
    n=read(); m=read(); R=read();
    for ( int i=1; i<=n; i++ )
        x[i]=read(),y[i]=read(),r[i]=read();
    for ( int i=1; i<=m; i++ )
        p[i]=read(),q[i]=read(),ansx+=x[i],ansy+=y[i];
    
    ansx/=m; ansy/=m; ansout=calc( ansx,ansy );
    while ( (double)clock()/CLOCKS_PER_SEC<=0.8 ) SA();

    printf( "%d",ansout );
}

習題1——UVA10228 A Star not a Tree?

題目連結

注意多組資料，計時的CLOCK是總時間。

如果用 while ( (double)clock()/CLOCKS_PER_SEC<=0.8 ) SA();

這樣的話第一組資料就跑了很多次SA,後面幾乎沒跑，就會出問題。。。

所以就改成固定次數了。

/*
給定一個N邊形所有頂點座標x,y，求其費馬點到所有頂點距離和
費馬點是指到多邊形所有頂點距離和最小的點
第一行為T， T組資料
第二行正整數N，其後N行，每行兩個整數x,y。
*/
#include <bits/stdc++.h>
using namespace std;
const int N=110;
const double delta=0.998,Te=1e-10;
int n,m,x[N],y[N];
double ansx,ansy;

int read()
{
    int x=0,f=1; char ch=getchar();
    while ( ch<'0' || ch>'9' ) 
    {
        if ( ch=='-' ) f=-1; ch=getchar();
    }
    while ( ch>='0' && ch<='9' ) x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    return x*f;
}

double dis( double ax,double ay,double bx,double by )
{
    double res=(bx-ax)*(bx-ax)+(by-ay)*(by-ay);
    res=sqrt(res);
    return res;
}

double calc( double xx,double yy )
{
    double res=0;
    for ( int i=1; i<=n; i++ )
        res+=dis( xx,yy,x[i],y[i] );
    return res;
}

void SA()
{
    double T=1e7; 
    while ( T>Te )
    {
        double nowx=ansx+(rand()*2-RAND_MAX)*T;
        double nowy=ansy+(rand()*2-RAND_MAX)*T;
        double nowans=calc( nowx,nowy ),ans=calc(ansx,ansy);
        double DelE=nowans-ans;
        if ( DelE<0 ) ansx=nowx,ansy=nowy;
        else if ( exp(-DelE/T)>(double)rand()/RAND_MAX ) ansx=nowx,ansy=nowy;
        T*=delta;
    }
}

int main()
{
    srand(time(0));  int Case=read();
    for ( int cas=0; cas<Case; cas++ )
    {
        if ( cas ) printf( "\n" );
        ansx=ansy=0; n=read(); 
        for ( int i=1; i<=n; i++ )
            x[i]=read(),y[i]=read(),ansx+=x[i],ansy+=y[i];
        
        ansx/=n; ansy/=n;
        for ( int i=1; i<=100; i++ ) SA();
        printf( "%.0lf\n",calc( ansx,ansy ) );
    }
}