啟發式算法 之 模拟退火原理及實踐

一、初窺其貌

1.1 啟發式算法和元啟發式算法

啟發式算法是求解優化問題的一類方法，因為經典優化方法存在局限性，有時無法得到最優解，隻能得到一個可以接受的近似最優解，啟發式算法就适合求解這類問題。啟發式算法就是專家的推測，借助過去經驗的歸納推理，以期求得問題的近似最優解或以一定機率求得最優解，是以認為啟發式算法是基于經驗或者實驗算法的統稱。

元啟發式算法則是從自然界的一些現象中獲得靈感，通過這些現象的特點解決實際問題，如模拟退火算法、遺傳算法等。

啟發式算法 = 元啟發式算法 + 問題特征

1.2 局部最優與全局最優

下圖是說明局部最優和全局最優的典型例子，綠色點是局部最優，紅色點是全局最優。

如果采用局部最優算法，從左到右尋找最小值點，到達綠色點時，算法認為已經找到最優，因為它無法忍受在繼續嘗試的過程中使值反而變大，綠色點左右兩側的值都比它大，算法無法跳出局部最優的情況。

如果采用全局最優算法，從左到右尋找最小值點，達到綠色點時，算法仍有機會跳出該點，因為算法允許情況有一定機率發生惡化，即允許下一次嘗試的值比綠色點的值大。當嘗試的值翻過綠點右側的山丘（極大值）時，則可以找到全局最優的紅點。

二、模拟退火

2.1 靈感來源

退火是冶金行業的術語，指的是将金屬加熱到一定溫度後，以适宜速度冷卻（通常較為緩慢），若以較快速度冷卻則稱為淬火。自然界中的物質，溫度越高能量越大，溫度越低能量越小，而采取緩慢冷卻的方式通常可以達到最低能量狀态，模拟退火就是從該現象得到的啟發。

2.2 算法思想

模拟退火（Simulated-Annealing, SA）算法的思想概括如下：

• 若目标函數由第 i 步移動到第 i + 1 步後的結果更優，即 f(Y( i + 1 )) <= f(Y(i))，總是接受該移動。
• 若沒有，即 f(Y( i + 1 )) > f(Y(i))，則以一定機率接受該移動，且該機率随着時間推移逐漸變小。

模拟退火期望得到的是全局最優解，是以它可以容忍結果以一定機率發生惡化，這個機率的計算采用的是 Metropolis 準則。

k 為常數，t 表示溫度，delta E 為兩次結果的能量內插補點，p 為接受較差結果的機率。其它參數相同時，溫度較高，即 t 較大，p 的值較大；溫度較低，即 t 較小，p 的值較小。說明在算法開始階段，接受差結果的機率較大，随着時間推移溫度降低，為了穩定于某個值，接受差結果的機率變小。

2.3 算法流程圖

(非本人所做，侵删)

三、代碼實作

3.1 旅行商問題

旅行商問題（Travelling Salesman Problem）是一道經典算法題，題幹不再贅述，它可以通過動态規劃求解，這裡采用模拟退火。

import numpy as np
import matplotlib.pyplot as plt
from numpy import random

# generate nums(configurable) citys data randomly
def gen_city_data(nums):
    citys = []
    N = 0
    while N < nums:
        # randn() generates float data conform to standard normal distribution, i.e. N(0,1).
        x = random.randn()
        y = random.randn()
        if [x, y] not in citys:
            citys.append([x, y])
            N = N + 1
    return np.array(citys)

# city-52 data
def init_city52():
    citys = np.array([[565.0, 575.0], [25.0, 185.0], [345.0, 750.0], [945.0, 685.0], [845.0, 655.0],
    [880.0, 660.0], [25.0, 230.0], [525.0, 1000.0], [580.0, 1175.0], [650.0, 1130.0],[1605.0, 620.0], 
    [1220.0, 580.0], [1465.0, 200.0], [1530.0, 5.0], [845.0, 680.0],[725.0, 370.0], [145.0, 665.0], 
    [415.0, 635.0], [510.0, 875.0], [560.0, 365.0],[300.0, 465.0], [520.0, 585.0], [480.0, 415.0], 
    [835.0, 625.0], [975.0, 580.0],[1215.0, 245.0], [1320.0, 315.0], [1250.0, 400.0], [660.0, 180.0], 
    [410.0, 250.0],[420.0, 555.0], [575.0, 665.0], [1150.0, 1160.0], [700.0, 580.0], [685.0, 595.0],
    [685.0, 610.0], [770.0, 610.0], [795.0, 645.0], [720.0, 635.0], [760.0, 650.0],[475.0, 960.0], 
    [95.0, 260.0], [875.0, 920.0], [700.0, 500.0], [555.0, 815.0],[830.0, 485.0], [1170.0, 65.0], 
    [830.0, 610.0], [605.0, 625.0], [595.0, 360.0],[1340.0, 725.0], [1740.0, 245.0]])
    ''' shape(citys) = (52, 2), 52 rows and 2 cols'''
    return citys

# travelling salesman problem
class TSP:
    def __init__(self):
        self._init_citys()
        self._init_dist()
        self._init_params()

    def _init_params(self):
        self.T = 280 # temprature
        self.Itr = 100 * self.n # amount of iterations at each temprature
        self.Alpha = 0.92 # SA cooling factor
        # initial state (0, 1, 2, ... , n-1), i.e. the sequence of salesman's travel.
        self.S = np.arange(self.n)

    def _init_citys(self):
        self.citys = init_city52()
        # self.citys = city_generator(10)
        self.n = self.citys.shape[0]

    # calculate the distance between citys
    def _init_dist(self):
        self.Dist = np.zeros((self.n, self.n))
        for i in range(self.n):
            for j in range(i, self.n):
                # linalg = linear algebra, norm() is 2-Norm in default.
                self.Dist[i][j] = self.Dist[j][i] = np.linalg.norm(self.citys[i] - self.citys[j])

    # calculate the energy in state S
    def energy(self, S):
        sum = 0
        for i in range(self.n - 1):
            sum = sum + self.Dist[S[i]][S[i+1]]
        sum = sum + self.Dist[S[self.n-1]][S[0]]
        return sum

    # randomly select a state from S's neighbors
    def neighbors2(self, S):
        S_neignbor = np.copy(S)
        u = np.random.randint(0, self.n)
        v = np.random.randint(0, self.n)
        while u == v:
            v = np.random.randint(0, self.n)
        if u > v:
            u, v = v, u
        t = S_neignbor[u:v]
        S_neignbor[u:v] = t[::-1] # ::-1 operate the sequence in reverse order
        return S_neignbor

    # simulated-annealing search
    def search(self):
        Temps = [] # cooling process
        Engys = [] # energy change process
        while self.T >= 0.1:
            print('search on temperature: {}'.format(self.T))
            # iterate at this temperature
            for _ in range(self.Itr):
                E_pre = self.energy(self.S)
                S_now = self.neighbors2(self.S)
                E_now = self.energy(S_now)
                # E_now better or Metropolis
                if (E_now < E_pre) or (np.exp((E_pre - E_now) / self.T) >= np.random.rand()):
                    self.S = S_now

            Temps.append(self.T)
            E_now = self.energy(self.S)
            Engys.append(E_now)
            print(E_now)

            # cooling down
            self.T = self.T * self.Alpha

        print(self.S)
        print('misson accomplished\n')
        return Temps, Engys

if __name__ == '__main__':
    tsp = TSP()
    Temps, Engys = tsp.search()
    plt.plot(Engys)
    plt.show()
           

3.2 退火效果