[學習筆記] 模擬退火 (Simulated Annealing)

真沒想到這東西真的在考場上用到了...順便水篇blog以示詐屍好了(逃

模擬退火算法

模擬退火是一種隨機化算法, 用於求函數的極值qwq

比如給出一個問題, 我們要求最優解的值, 但是可能的方案數量極大, 直接搜索會T飛(或者方案是連續的總數無窮根本沒法搜), 這種時候我們一般會有兩種選擇:

爬山算法

爬山算法每次在當前找到的方案附近尋找一個新的方案(常見方式是隨機一個差值), 然后如果這個解更優那么直接轉移.

對於單峰函數來說這顯然可以直接找到最優解(不過你都知道它是單峰函數了為啥不三分呢?)

但是對於多數我們求解的函數來說, 它並不一定會長成這個樣子...於是就極其有可能鑽進一個局部的最優解出不來了

算法得出的最優解與初始解的位置以及搜尋的附近解的區域大小有關.

當然如果你尋找新方案的區間很大的話有概率跳出去, 但是太大的話又可能跳來跳去跳亂了從而找不到最優解...

歐皇專用最優化求解方式(@liu_runda)

 

然而並不是所有人都是歐皇, 像博主這樣的非酋要怎么辦捏?

當然是求助於自然規律(大霧)

退火的理論部分

退火其實本來是冶金工業里的術語...大概過程是先把晶體加熱到極高的溫度再緩慢降溫, 在這個過程中減少晶體中的缺陷(達到能量最低的最穩定狀態)

具體理論是這樣的:

設 $E[\{x_i\}]$ 表示某一物質體系在微觀狀態 $\{x_i\}$ 下的內能, 對於給定溫度 $T$, 若體系處於熱平衡態時, $E[\{x_i\}]$ 服從 Boltzmann 分布:

$$ f=c(T)\exp\left(-\frac{E[\{x_i\}]}{kT}\right) $$

其中 $k$ 為 Boltzmann 常數.

$T$ 下降, $E$ 隨之下降. 若 $T$ 下降的速度足夠慢, 則體系總可以保持熱平衡態. 當 $T\to 0$ 時 $E$ 趨近於最小. 這樣的物質降溫過程被稱為退火過程.

然后機智的我們發現這個過程最終和我們的最優化過程類似!

於是我們去模擬這個過程, 按照退火的規律引入更多隨機因素, 這樣我們得到最優解的概率就會增加辣233

可以證明, 模擬退火得到的最終解依概率收斂於最優解.

emmmm...

等等模擬這個過程? 這是計算機又不是實驗室你怎么模擬啊(╯°□°)╯︵┻━┻

拿出物理化學(假裝自己還是個ChO黨)...

根據熱力學規律並結合計算機對離散數據的處理, 我們定義: 如果當前溫度為 $T$ , 當前狀態與新狀態之間的能量差為 $\Delta E$ , 則發生狀態轉移的概率為:

$$ P(\Delta E) = e^{\frac { \Delta E } { kT } } $$

顯然如果 $ \Delta E$ 為正的話轉移是一定會成功的, 但是對於 $\Delta E < 0$ 我們則以上式中計算得到的概率接受這個新解.

然后我們維護溫度 $T$ 即可. 這里我們有三個參數: 初溫 $T_0$ , 降溫系數 $d$ , 終溫 $T_k$

一般 $T_0$ 是個比較大的數, $d$ 是個接近 $1$ 但是小於 $1$ 的值, $T_k$ 是個接近 $0$ 的正值.

首先讓溫度 $T=T_0$ , 然后進行一次轉移嘗試, 然后讓 $T=dT$.

當 $T<T_k$ 時模擬退火過程結束, 當前解作為最優解.

看起來好像並不是很難理解?

Wikipedia上的動圖:

一般來講模擬退火在參數合適的情況下效果拔群, TSP隨便跑(x

模擬退火的實際使用

實際使用里這函數可不一定是個單元函數...而且尋找新解好像是個很模糊的東西, 畢竟很多時候我們會發現我們要求解的問題的所有可能解並不是離散的...

先拿道題說說...

3680: 吊打XXX

Time Limit: 10 Sec  Memory Limit: 128 MBSec  Special Judge
Submit: 4247  Solved: 1566
[Submit][Status][Discuss]

Description

gty又虐了一場比賽，被虐的蒟蒻們決定吊打gty。gty見大勢不好機智的分出了n個分身，但還是被人多勢眾的蒟蒻抓住了。蒟蒻們將
n個gty吊在n根繩子上，每根繩子穿過天台的一個洞。這n根繩子有一個公共的繩結x。吊好gty后蒟蒻們發現由於每個gty重力不同，繩
結x在移動。蒟蒻wangxz腦洞大開的決定計算出x最后停留處的坐標，由於他太弱了決定向你求助。
不計摩擦，不計能量損失，由於gty足夠矮所以不會掉到地上。

Input

輸入第一行為一個正整數n(1<=n<=10000)，表示gty的數目。
接下來n行,每行三個整數xi，yi，wi，表示第i個gty的橫坐標，縱坐標和重力。
對於20%的數據，gty排列成一條直線。
對於50%的數據，1<=n<=1000。
對於100%的數據，1<=n<=10000,-100000<=xi,yi<=100000

Output

輸出1行兩個浮點數（保留到小數點后3位），表示最終x的橫、縱坐標。

Sample Input

3
0 0 1
0 2 1
1 1 1

Sample Output

0.577 1.000

HINT

 

Source

By wangxz

在這個題中我們為了到達穩定狀態要讓整個體系的總重力勢能最低.

重力勢能怎么求呢? 別忘了這繩子的總長度是不會變的...於是某個質點的重力勢能和到繩結的水平距離成一次函數關系.

我們為了簡化問題, 可以將某個質點對最終答案產生的貢獻計算為 $ dis(o,x)\times m_x $ . 然后我們要讓這個值最小化.

這個時候我們可以考慮模擬退火. 首先隨機一個點作為初始解(為了加速收斂, 我們可以直接取各個點坐標的平均值所在的店). 然后隨機兩個值作為差值加到這個點的坐標上作為下一個解.

然后模擬退火直接往上套就可以了233

具體實現就是一個  while  循環, 循環內有4步:

1. 根據當前解找到下一個解
2. 計算下一個解的 "能量" (也就是價值)
3. 決定是否要接受這個新解
4. 降溫

找下一個解的時候有一個提高精度的小技巧: 根據當前溫度決定差值的范圍. 這樣在降溫即將結束接近最優解的時候可以有更大的概率更精確地命中最優解.

當然如果是解是離散的就不能這樣搞了. 以及生成下一個解的時候萬萬不能全部重新隨機生成, 那就和瞎隨沒區別了...要隨機作出一些相對小的修改.

具體做法就是使用一個產生 $[0,1]$ 隨機實數的函數, 將隨機區間轉為 $[-1,1]$ 后乘上 $T$ 作為差值. (也就是生成一個 $[-T,T]$ 的隨機值作為差值)

不過實際操作的時候我們較少直接輸出最終解, 而是選擇在模擬退火的過程中單獨維護一個解, 只在遇到更優解的時候將其更新, 增加正確率.

另一個提高正確率的方法是: 多次進行模擬退火過程(或者說"重新燒熱再退火一遍"), 每次取最優解.

還有就是最后燒完之后可以再在全局最優解的基礎上進行爬山.

本題的參考實現:

#include <bits/stdc++.h>
 
const int MAXN=1e4+10;
 
struct Point{
    double x;
    double y;
    Point(double x=0,double y=0){
        this->x=x;
        this->y=y;
    }
};
Point P[MAXN];
 
int n;
Point ans;
int g[MAXN];
double minAns=DBL_MAX;
 
double Rand();
double Sqr(double);
double Calc(const Point&);
bool Accept(double,double);
double EucDis(const Point&,const Point&);
Point SimulatedAnnealing(Point,double,double,double);
 
int main(){
    scanf("%d",&n);
    Point init;
    for(int i=0;i<n;i++){
        scanf("%lf%lf%d",&P[i].x,&P[i].y,g+i);
        init.x+=P[i].x;
        init.y+=P[i].y;
    }
    init.x/=n;
    init.y/=n;
    SimulatedAnnealing(init,1e5,1-7e-3,1e-3);
    printf("%.3f %.3f\n",ans.x,ans.y);
    return 0;
}
 
inline double Calc(const Point& origin){
    double ans=0;
    for(int i=0;i<n;i++){
        ans+=EucDis(origin,P[i])*g[i];
    }
    if(ans<minAns){
        ::ans=origin;
        minAns=ans;
    }
    return ans;
}
 
bool Accept(double delta,double tmp){
    return delta<0||Rand()<exp(-delta/tmp);
}
 
Point SimulatedAnnealing(Point initAns,double initT,double dec,double end){
    double tmp=initT;
    Point now=initAns;
    double nowAns=Calc(now);
    while(tmp>end){
        Point next=Point(now.x+tmp*(Rand()*2-1),now.y+tmp*(Rand()*2-1));
        double ans=Calc(next);
        if(Accept(ans-nowAns,tmp)){
            now=next;
            nowAns=ans;
        }
        tmp*=dec;
    }
    for(int i=0;i<1000;i++){
        Point rnd=Point(ans.x+tmp*(Rand()*2-1),ans.y+tmp*(Rand()*2-1));
        Calc(rnd);
    }
    return now;
}
 
inline double Rand(){
    return double(rand())/double(RAND_MAX);
}
 
inline double Sqr(double x){
    return x*x;
}
 
inline double EucDis(const Point& a,const Point& b){
    return sqrt(Sqr(a.x-b.x)+Sqr(a.y-b.y));
}

但是調參的過程還是比較看臉的...究極非洲大酋長慎用(x

一般情況下我們會在時間允許的情況下盡量多地嘗試新的解. 一般降溫系數 $d$ 與 $1$ 的差減少一個數量級, 時間消耗大約多 $10$ 倍, $T_0$ 和 $T_k$ 變化一個數量級, 時間消耗不會變化很大.

這種時候我們可以試着先本機跑跑自造數據看看精度怎么樣. 如果發現經常陷入局部最優解的話考慮增大 $T_0$ 和 $d$ , 如果發現最終精度不夠的話考慮減小 $T_k$.

至於模擬退火的正確率計算么...好像只有實驗是最方便的了吧(x

今天上午考試的時候手調一波參數然后極限數據下測試 $100$ 次發現精度達標率有 $60\%$ 就交了...然后A了...

於是借此在高二那邊水了個 $\texttt{rk4}$ (逃

 

模擬退火解旅行商問題

剛剛說模擬退火TSP隨便跑...那么我們就說說TSP怎么跑...

有人可能會問了這個東西怎么求下一個解?

其實還是隨機...

對於TSP, 我們的一個方案其實就是一個遍歷順序(也就是一個排列)

這時我們在生產新解的時候可以選擇隨機選取兩個結點, 然后將它們在排列中的位置交換一下(好暴力啊)

然而事實證明效果拔群...

 

總結

模擬退火在OI中是一種在最優化問題中騙分的好方法

對於一些奇奇怪怪的多元函數也可以用這個方法來求解

其實在上面的例子中也可以體現出來, 這個算法的要點在於新解的選取以及參數的調整...

實際上利用退火過程的性質大膽隨機再配合調參經驗一般效果拔群OwO

但是作為一個隨機化算法並不一定能找到最優解qwq...IOI賽制/ACM賽制的話可能騙分更容易些?(畢竟可以多次提交233)