模擬退火略解

模擬退火

模擬退火 算法（Simulate Anneal，SA）是一種通用機率演算法，用來在一個大的搜尋空間內找尋命題的最優解。
模擬退火的出發點是基於物理中 固體物質的退火過程 與通常組合優化問題之間的類似性。

固體物質的退火過程

lz前天的市物理競賽爆炸（被防AK了）……心情極度不爽
據說oi同志們物理都很好 （翹物理競賽輔導課） ，仍是提一下這東西吧。

根據牛頓冷卻定律 (Newton's law of cooling)，當物體表面與周圍存在溫度差時，單位時間從單位面積散失的熱量與溫度差成正比，即
\[\Delta\tau=-k(\tau-\tau_0)\]其中 \(k\) 爲恆量。
式中 \(\tau\) 表示物體溫度， \(\tau_0\) 表示環境溫度（不變量）。
換句話說，單位時間散失的熱量\[\Delta\tau\propto(\tau-\tau_0)\]



下面咱們以 \(\text{luogu P1337}\) 爲例，講解模擬退火算法的原理和步驟。

題目描述 \(\text{luoguP1337}\)

如圖：有 \(n\) 個重物，每一個重物系在一條足夠長的繩子上。每條繩子自上而下穿過桌面上的洞，而後系在一塊兒。圖中 \(X\) 處就是公共的繩結。假設繩子是徹底彈性的（不會形成能量損失），桌子足夠高（於是重物不會垂到地上），且忽略全部的摩擦。

問繩結 \(X\) 最終平衡於何處。

注意：桌面上的洞都比繩結 \(X\) 小得多，因此即便某個重物特別重，繩結 \(X\) 也不可能穿過桌面上的洞掉下來，最可能是卡在某個洞口處。

輸入格式

文件的第一行爲一個正整數 \(n\)（\(1≤n≤1000\)），表示重物和洞的數目。

接下來的 \(n\) 行，每行是 \(3\) 個整數：\(X_i.Y_i.W_i\)，分別表示第 \(i\) 個洞的座標以及第 \(i\) 個重物的重量。\((-10000≤x,y≤10000, 0<w≤1000 )\)

輸出格式

你的程序必須輸出兩個浮點數（保留小數點後三位），分別表示處於最終平衡狀態時繩結 \(X\) 的橫座標和縱座標。兩個數以一個空格隔開。

輸入樣例

3
0 0 1
0 2 1
1 1 1

輸出樣例

0.577 1.000

\(\text{Solution 1337}\)

定義 \(f(x,y)\) 表示 \(X(x,y)\) 時的答案爲多少。
容易獲得空間能量和\[f(x,y)=\sum_{k=1}^{n}{\sqrt{((X_k-x)^2+(Y_k-y)^2)}\times W_k}\]
由於在一個空間中的能量和越小，該空間越穩定。因此 \(f_{\min}\) 時的 \(X\) 即爲答案。


觀察模擬退火時當前最優解與溫度的關係。

能夠看出，當 \(\tau→0\) 時，咱們找到最優解的概率也愈來愈大。
對於每一個溫度 \(\tau\)，嘗試找一個新解。
若新解更優，則接受；若新解次，則以必定機率接受，這個機率爲

\[e^{\frac{\Delta ans}{k\tau}}\]

其中 \(k\) 是 \(0\) 到 \(1\) 之間的隨機數。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<ctime>

#define reg register

int n;
struct node{
    int x,y,w;
}e[1010];
int cx=0,cy=0;
double ansx,ansy,ans=1e18;

double calc(double x,double y){ //計算f
    double cnt=0;
    for(reg int i=1;i<=n;++i)
        cnt+=sqrt((x-e[i].x)*(x-e[i].x)+(y-e[i].y)*(y-e[i].y))*e[i].w;
    return cnt;
}
void SA(){
    double x=ansx,y=ansy;
    double t=2000.0;    //初溫
    while(t>1e-14){     //末溫
        //獲得新解座標
        double nx=x+((rand()*2)-32767)*t;
        double ny=y+((rand()*2)-32767)*t;
        double nw=calc(nx,ny);
        double delta=nw-ans;    //delta ans
        if(delta<0){    //新解更優，接受
            x=nx;y=ny;ans=nw;
            ansx=nx;ansy=ny;
        }
        //新解更次，以必定機率接受
        else if(exp(-delta/t)*32767>rand()) x=nx,y=ny;
        t*=0.993;   //根據牛頓冷卻定律降溫
    }
}
void work(){
    double px=(double)cx/n,py=(double)cy/n;
    ansx=px;ansy=py;
    //屢次執行SA，獲得正確答案的概率更大
    for(reg int i=1;i<=10;++i) SA();
}
int main(){
    srand(100007);
    scanf("%d",&n);
    for(reg int i=1;i<=n;++i){
        scanf("%d%d%d",&e[i].x,&e[i].y,&e[i].w);
        cx+=e[i].x;cy+=e[i].y;
    }
    //從平均值開始的話，與答案更接近
    work();
    printf("%.3lf %.3lf",ansx,ansy);
}

算法改進

（1）設計合適的狀態產生函數，使其根據搜索進程的須要表現出狀態的全空間分散性或局部區域性；

（2）設計高效的退火策略；

（3）避免狀態的迂迴搜索；

（4）採用並行搜索結構；

（5）爲避免陷入局部極小，改進對溫度的控制方式；

（6）選擇合適的初始狀態；

（7）設計合適的算法終止準則。

題目描述 \(\text{luogu P4360}\)

從山頂上到山底下沿着一條直線種植了 \(n\) 棵老樹。當地的政府決定把他們砍下來。爲了避免浪費任何一棵木材，樹被砍倒後要運送到鋸木廠。

木材只能朝山下運。山腳下有一個鋸木廠。另外兩個鋸木廠將新修建在山路上。你必須決定在哪裏修建這兩個鋸木廠，使得運輸的費用總和最小。假定運輸每公斤木材每米須要一分錢。

你的任務是編寫一個程序，從輸入文件中讀入樹的個數和他們的重量與位置，計算最小運輸費用。

輸入格式

輸入的第一行爲一個正整數 \(n\) ——樹的個數 \((2≤n≤20000)\)。樹從山頂到山腳按照 \(1,2,...,n\) 標號。

接下來 \(n\) 行，每行有兩個正整數（用空格分開）。
第 \(i+1\) 行含有：\(w_i\) ——第 \(i\) 棵樹的重量（公斤爲單位）和 \(d_i\)——第 \(i\) 棵樹和第 \(i+1\) 棵樹之間的距離，\(1≤w_i≤10000,0≤d_i≤10000\)。

最後一顆樹的 \(d_n\) 表示第 \(n\) 棵樹到山腳的鋸木廠的距離。保證全部樹運到山腳的鋸木廠所須要的費用小於 \(2×10^9\) 分。

輸出格式

輸出最小的運輸費用。

輸入樣例

9 
1 2 
2 1 
3 3 
1 1 
3 2 
1 6 
2 1 
1 2 
1 1

輸出樣例

26

說明

樣例圖示

黑點爲鋸木廠。

\(\text{Solution P4360}\)

定義 \(f(a,b)\) 表示鋸木廠在 \(a,b\ (a<b)\) 時的答案。
設第 \(i\) 棵樹到山腳的距離爲 \(D_i\)，則\[D_i=\sum_{j=1}^{i-1}{d_j}\]
則\[\begin{aligned}f(a,b)&=\sum_{i=1}^{a}{w_i·(D_a-D_i)}+\sum_{i=a+1}^{b}{w_i·(D_b-D_i)}+\sum_{i=b+1}^{n}{D_{n+1}-D_i}\\ &=D_a·\sum_{i=1}^{a}{w_i}+D_b·\sum_{i=a+1}^{b}{w_i}+D_{n+1}·\sum_{i=b+1}^{n}{D_i}-\sum_{i=1}^{n}{w_i·D_i} \end{aligned}\]
設\[W_k=\sum_{i=1}^{k}{w_i}\]\[s_k=\sum_{i=1}^{k}{w_i·D_i}\]咱們發現它們可使用前綴和優化，因此咱們能夠用 \(\text{O(1)}\) 的時間複雜度求出 \(f\)。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<ctime>

#define reg register
typedef long long ll;

ll w[200010],d[200010],s;
int n;
struct node{
    int a,b;
    double f(){
        if(a>b) a^=b^=a^=b;
        return w[a]*d[a]+(w[b]-w[a])*d[b]+(w[n]-w[b])*d[n+1]-s;
    }
}x;
int cx=0,cy=0;
int ansx,ansy;
double ans=1e17;
int s1,s2;

void SA(){
    double a=1,b=2;
    double t=(double)n;
    while(t>5e-1){
        int nx=round(a+((rand()*2)-32767)*t);
        int ny=round(b+((rand()*2)-32767)*t);
        x.a=(nx+n)%n;x.b=(ny+n)%n;
        double nw=x.f();
        double delta=nw-ans;
        if(delta<0){
            a=nx;b=ny;ans=nw;
            ansx=(nx+n)%n;ansy=(ny+n)%n;
        }
        else if(exp(-delta/t)*32767>rand()) a=(nx+n)%n,b=(ny+n)%n;
        t*=0.99;
    }
}
void work(){
    for(reg int i=1;i<=10;++i) SA();
}
int main(){
    srand(100007);
    scanf("%d",&n);
    w[0]=d[0]=s=0;
    for(reg int i=1;i<=n;++i){
        scanf("%d%d",&s1,&s2);
        w[i]=w[i-1]+s1;d[i+1]=d[i+1]+s2;
        s=s+s1*d[i];
    }
    work();
    x.a=ansx;x.b=ansy;
    printf("%lld",(ll)x.f());
}

參考資料

序號	標題
\(1\)	淺談玄學算法——模擬退火
\(2\)	題解 P4360 [CEOI2004]鋸木廠選址

題表

序號	標題	題解
\(1\)	P2503 [HAOI2006]均分數據	未解決
\(2\)	P4035 [JSOI2008]球形空間產生器	已解決
\(3\)	P3878 [TJOI2010]分金幣	已解決
\(4\)	SP4587 FENCE3 - Electric Fences	未解決
\(5\)	CF2C Commentator problem	未解決
\(6\)	UVA10228 A Star not a Tree?	已解決
\(7\)	P3936 Coloring	已解決
\(8\)	P2210 Haywire	已解決