[學習筆記] 模擬退火 (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]$ 的隨機值做爲差值)

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

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

還有就是最後燒完以後能夠再在全局最優解的基礎上進行登山.

本題的參考實現:

 1 #include <bits/stdc++.h>
 2  
 3 const int MAXN=1e4+10;
 4  
 5 struct Point{
 6     double x;
 7     double y;
 8     Point(double x=0,double y=0){
 9         this->x=x;
10         this->y=y;
11     }
12 };
13 Point P[MAXN];
14  
15 int n;
16 Point ans;
17 int g[MAXN];
18 double minAns=DBL_MAX;
19  
20 double Rand();
21 double Sqr(double);
22 double Calc(const Point&);
23 bool Accept(double,double);
24 double EucDis(const Point&,const Point&);
25 Point SimulatedAnnealing(Point,double,double,double);
26  
27 int main(){
28     scanf("%d",&n);
29     Point init;
30     for(int i=0;i<n;i++){
31         scanf("%lf%lf%d",&P[i].x,&P[i].y,g+i);
32         init.x+=P[i].x;
33         init.y+=P[i].y;
34     }
35     init.x/=n;
36     init.y/=n;
37     SimulatedAnnealing(init,1e5,1-7e-3,1e-3);
38     printf("%.3f %.3f\n",ans.x,ans.y);
39     return 0;
40 }
41  
42 inline double Calc(const Point& origin){
43     double ans=0;
44     for(int i=0;i<n;i++){
45         ans+=EucDis(origin,P[i])*g[i];
46     }
47     if(ans<minAns){
48         ::ans=origin;
49         minAns=ans;
50     }
51     return ans;
52 }
53  
54 bool Accept(double delta,double tmp){
55     return delta<0||Rand()<exp(-delta/tmp);
56 }
57  
58 Point SimulatedAnnealing(Point initAns,double initT,double dec,double end){
59     double tmp=initT;
60     Point now=initAns;
61     double nowAns=Calc(now);
62     while(tmp>end){
63         Point next=Point(now.x+tmp*(Rand()*2-1),now.y+tmp*(Rand()*2-1));
64         double ans=Calc(next);
65         if(Accept(ans-nowAns,tmp)){
66             now=next;
67             nowAns=ans;
68         }
69         tmp*=dec;
70     }
71     for(int i=0;i<1000;i++){
72         Point rnd=Point(ans.x+tmp*(Rand()*2-1),ans.y+tmp*(Rand()*2-1));
73         Calc(rnd);
74     }
75     return now;
76 }
77  
78 inline double Rand(){
79     return double(rand())/double(RAND_MAX);
80 }
81  
82 inline double Sqr(double x){
83     return x*x;
84 }
85  
86 inline double EucDis(const Point& a,const Point& b){
87     return sqrt(Sqr(a.x-b.x)+Sqr(a.y-b.y));
88 }

View Code

可是調參的過程仍是比較看臉的...究極非洲大酋長慎用(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)