启发式算法 之 模拟退火原理及实践

一、初窥其貌

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 退火效果