【智能算法】模擬退火算法

1.前言

      隨着機器學習的發展，在求解一些問題的過程當中，問題的規模逐漸加大，對這些問題的求解到很精確的解，須要花費大量的資源，還不必定能獲得，當碰到這些問題時，能夠經過必定的優化算法，獲得問題的最優解或者近似值。例如：

      旅行商問題（traveling salesman proplem TSP）：給定一系列城市和每對城市之間的距離，求解訪問每一座城市一次並回到起始城市的最短迴路。

      超大規模集成電路(Very Large Scale Integration VLSI)：一般須要在幾毫米見方大小的硅片上集成上萬至百萬晶體管、線寬在1微米如下的集成電路。因爲晶體管與連線一次完成，故製做幾個至上百萬晶體管的工時和費用是等同的。大量生產時，硬件費用幾乎可不計，而取決於設計費用。用模擬退火算法幾乎能夠很好地完成全部優化的VLSI設計工做。如全局佈線、布板、佈局和邏輯最小化等等。

      圖像識別與處理問題：可用來進行圖像恢復等工做，即把一幅被污染的圖像從新恢復成清晰的原圖，濾掉其中被畸變的部分。

如此以上這些領域都是用模擬退火算法在求解，能產生使人滿意的近似最優解，並且所用的時間也不很長。既然模擬退火算法有這麼強大的效果，各位看官是否是火燒眉毛的想進入正題了，我跟大家說，先別急，在正式讓模擬退火算法隆重登場以前，咱們先介紹一種簡單的貪心搜索算法----登山算法。

2.登山算法

      前面咱們已經說了，登山算法是一種簡單的貪心搜索算法，該算法每次從當前解的臨近解空間中選擇一個最優解做爲當前解，直到達到一個局部最優解。

      登山算法實現很簡單，其主要缺點是會陷入局部最優解，而不必定能搜索到全局最優解。如圖1所示：假設C點爲當前解，登山算法搜索到A點這個局部最優解就會中止搜索，由於在A點不管向那個方向小幅度移動都不能獲得更優的解。

圖 1

下面咱們具體看一個登山算法的例子

求解下面的式子

                                                     （公式-1）

下面經過Python實如今x在區間[-2,4]和y在區間[-2,4]的圖像

# -*- coding: utf-8 -*-
from mpl_toolkits.mplot3d  import Axes3D
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = Axes3D(fig)
X = np.arange(-2, 4, 0.05)
Y = np.arange(-2, 4, 0.05)
X, Y = np.meshgrid(X, Y)
Z = np.exp(-(np.power(X,2)+np.power(Y,2)))+2*np.exp(-(np.power(X-1.7,2)+np.power(Y-1.7,2)))

ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='hot')
plt.show()

生成的圖像以下：

圖 2

從圖像中發現，這個函數具備兩座山峯，一高一矮，而後咱們寫一個登山算法，基本思路爲：先隨機選擇一個起點，今後點開始爬山，依次尋找該點四周的點，判斷是否比此點高，不斷循環，直到四周沒有比其更高的點爲止，此時找到的點即爲最優解（最大值）。

使用Python寫一個登山函數

def argmax(stx, sty, alisx=0, alisy=0):
    cur = func(stx, sty)
    next = func(alisx, alisy)

    return cur < next and True or False

def climb(X,Y):
    global_X = []
    global_Y = []

    len_x = len(X)
    len_y = len(Y)
    # 隨機爬山點
    st_x = random.randint(0, len_x-1)
    st_y = random.randint(0, len_y-1)

    while (len_x > st_x >= 0) or (len_y > st_y >= 0):
        if st_x + 1 < len_x and argmax(X[0][st_x], Y[st_y][0], X[0][st_x +1],Y[st_y][0]):
            st_x += 1
        elif st_y + 1 < len_x and argmax(X[0][st_x], Y[st_y][0], X[0][st_x +0], Y[st_y + 1][0]):
            st_y += 1
        elif st_x >= 1 and argmax(X[0][st_x], Y[st_y][0], X[0][st_x -1],Y[st_y][0]):
            st_x -= 1
        elif st_y >= 1 and argmax(X[0][st_x], Y[st_y][0], X[0][st_x +0], Y[st_y - 1][0]):
            st_y -= 1
        else:
            break
        global_X.append(X[0][st_x])
        global_Y.append(Y[st_y][0])
    return global_X, global_Y, func(X[0][st_x], Y[st_y][0])

調用此登山函數，並把登山路徑描繪出來

px, py, maxhill = climb(X, Y)
print maxhill
ax.plot(px,py,func(np.array(px),np.array(py)),'b')

經過屢次運行腳本，獲得的如下幾種可能的圖像（由於是隨機值，你在本地運行的登山路徑具體可能不同）

序號	圖像	最大值
1		2.00308871541
2		2.00308871541
3		1.00617743082

從上面表格中序號爲3的圖中，可以清楚的看到，登山的路徑不必定會爬到最高的山頂（全局最優解），有可能會爬到局部的山頂（局部最優解），因此這種算法存在這種缺陷。聰明的科研人員提出了一種改進這種算法的方法，噹噹噹……，下面就有請模擬退火算法登場。

3.模擬退火算法基本思想

模擬退火（simulated annealing SA）算法的思想最先是由Metropolis等提出的，其出發點是基於物理中固體物質的退火過程與通常問題的組合優化問題的類似性。其物理過程由三部分組成：

（1）加溫過程。使粒子增強運動，使其偏離平衡位置。當溫度足夠高時，固體變爲液體，從而消除系統原先存在的非均勻狀態。

（2）等溫過程。溫度不變，系統狀態的自發變化朝自由能減小的方向進行，當自由能達到最小時，保持平衡。

（3）冷卻過程。使粒子熱運動減弱，系統能量下降，獲得晶體結構。

通俗的講就是：溫度越高，出現一次能量差的降溫的機率就越大；溫度越低，則出現降溫的機率就越小。模擬退火算法以必定的機率來接受一個比當前解要差的解，所以有可能會跳出這個局部的最優解，達到全局的最優解。以圖1爲例，模擬退火算法在搜索到局部最優解A後，會以必定的機率接受到E的移動。也許通過幾回這樣的不是局部最優的移動後會到達D點，因而就跳出了局部最大值A。

關於登山算法與模擬退火，有一個有趣的比喻：

    登山算法：兔子朝着比如今高的地方跳去。它找到了不遠處的最高山峯。可是這座山不必定是珠穆朗瑪峯。這就是登山算法，它不能保證局部最優值就是全局最優值。

    模擬退火：兔子喝醉了。它隨機地跳了很長時間。這期間，它可能走向高處，也可能踏入平地。可是，它漸漸清醒了並朝最高方向跳去。這就是模擬退火。

SA算法實現過程：

（1）初始化：取初始溫度T0足夠大，令T=T0.隨機取一個值，得初始解S1，肯定每一個T時的迭代數，即Metropolis鏈長L。

（2）對當前溫度T和k=1,2……，L，重複步驟（3）-（6）

（3）對當前解S1隨機產生一個新解S2

（4）計算S2的增量df=f(S2)-f(S1)，其中f(S1)爲S1的代價函數

（5）若df<0，則接受S2做爲新的當前解，即S1=S2；不然計算S2的接受機率exp(-df/T)，即隨機產生（0,1）區間上均勻分佈的隨機數rand，若exp（-df/T)>rand，也接受S2做爲新的當前解，S1=S2；不然保留當前解S1

（6）若是知足終止條件Stop，則輸出當前解S1做爲最優解，結束程序。終止條件Stop一般爲：在連續若干個Metropolis鏈中新解S2都沒有被接受時，終止算法；或者設定的溫度結束。不然繼續按衰減函數衰減T後返回步驟（2）。

4.模擬退火算法案例與實現

4.1 算法流程

流程圖中Metropolis準則爲：

若是df<0，以1的機率接受值，不然，以機率exp(-df/T)的機率接受新的值。

4.2 算法Python實現

def climb_SA(X,Y):
    global_X = []
    global_Y = []

    # 初始溫度
    temperature = 105.5
    # 溫度降低的比率
    delta = 0.98
    # 溫度精確度
    tmin = 1e-10

    len_x = len(X)
    len_y = len(Y)
    # 隨機爬山點
    st_x = X[0][random.randint(0, len_x - 1)]
    st_y = Y[random.randint(0, len_y-1)][0]
    st_z = func(st_x, st_y)

    while (temperature > tmin):
        # 隨機產生一個新的鄰近點
        # 說明: 溫度越高，鄰近點跳躍的幅度越大
        tmp_x = st_x + (random.random() * 2 - 1) * temperature
        tmp_y = st_y + + (random.random() * 2 - 1) * temperature
        if 4 > tmp_x >= -2 and 4 > tmp_y >= -2:
            if argmax(st_x, st_y, tmp_x, tmp_y):
                st_x = tmp_x
                st_y = tmp_y
            else:  # 有機會跳出局域最優解
                pp = 1.0 / (1.0 + np.exp(-(func(tmp_x, tmp_y) - func(st_x, st_y)) / temperature))
                if random.random() < pp:
                    st_x = tmp_x
                    st_y = tmp_y
        temperature *= delta  # 以必定的速率降低
        global_X.append(st_x)
        global_Y.append(st_y)
    return global_X, global_Y, func(st_x, st_y)

調用此模擬退火算法函數，並把點描繪出來：

px,py,maxhill = climb_SA(X,Y)
print maxhill
ax.plot(px,py,func(np.array(px),np.array(py)),'b.')

經過運行腳本：

打印的結果爲：2.00264256484

再次運行

打印的結果爲：2.00063593674

這樣屢次運行，觀察發現，山頂處匯聚了大量的點，山腳下散落着少許的點，最終基本都能獲得問題的近似最優解。

5.算法評價

模擬退火算法是一種隨機算法，並不必定能找到全局的最優解，能夠比較快的找到問題的近似最優解。 若是參數設置得當，模擬退火算法搜索效率比窮舉法要高。可是，模擬退火算法其實就是靠大數定律才成立的一個方法，初始解和產生新解，判斷新解是否可以接受都是靠隨機數，都是靠機率，因此期間就會產生大量的解，當問題的規模不斷變大的時候，要想獲得質量高的解，搜索的解的數目就會指數級的增長，無限膨脹算法的計算時間。後續可考慮將算法的尋優過程由串行改成並行。