【模擬退火】旅行商問題

1.問題描述

旅行商問題（Travelling Salesman Problem, 簡記TSP，亦稱貨郎擔問題)：設有n個城市和距離矩陣D=[d_ij]，其中d_ij表示城市i到城市j的距離，i，j=1，2 … n，則問題是要找出遍訪每一個城市剛好一次的一條迴路並使其路徑長度爲最短。

利用模擬退火算法（Simulate Anneal，SA）解決上述旅行商問題。

2.算法設計

1)  問題分析

在TSP問題中，n個城市的任何一種排列均是問題的一個可能解。其中表示第i個到達的城市是j，而且默認。當規定了出發城市時，解空間的規模是(n-1)!。所以難以利用深度搜索或廣度搜索等窮舉方法進行求解。

模擬退火算法來源於固體退火原理，是一種基於機率的算法。從某一較高初溫出發，伴隨溫度參數的不斷降低,結合機率突跳特性在解空間中隨機尋找目標函數的全局最優解，即在局部最優解能機率性地跳出並最終趨於全局最優。所以可以有效得解決TSP問題。

2)  參數設置

模擬退火算法的表現和參數設置有着很大關係。在通過實驗分析以後，實驗肯定了一組較爲合適的參數，該組參數權衡了運行時間和解空間的質量等因素。參數以下：

• 初始溫度t0 = 280；
• 每一個溫度下采用固定的迭代次數Lk = 100n，，n爲城市數量；
• 溫度衰減係數α = 0.92，即tk+1 = 0.92 * tk；
• 算法中止準則：溫度t < 0.1。

3)  新解產生

實驗中採用2變換法來產生新解：任選序號u和v (u<v)，將u和v及其之間的順序逆轉，即變爲。

4)  指標函數

實驗中的指標函數設置爲城市序列的總距離。指標函數差爲新解與原解之間的指標函數差。即有：

 

5)  新解接受準則

3.程序流程

 

4.核心僞代碼

Procedure TSP_USE_SIMULATED_ANNEALING;

begin

LOAD_DATA;   { 加載城市數據[(x1, y1), (x2, y2), … , (xn, yn)] }

INITIALIZE;    { 初始化 }

Repeat

    for l := 1 to L do

    begin

       GENERATE;          { 採用2變換法產生新的路徑 }

       CALCULATE(ΔD);     { 計算新路徑下的總距離差}

       if ACCEPT(t, ΔD) then  { 根據Metropolis準則判斷是否接受新的路徑 }

       begin

           S_new := S_now;  { 更新可能解爲當前城市序列 }

       end;

    end;

    t := t * α;    { 衰減函數 }

until t < 0.1      { 中止準則爲溫度小於0.1 }

End;

5.代碼運行及測試

    實驗中採用兩個數據集來測試模擬退火算法在TSP問題上的表現。第一個數據集CITY-52來源於網上，包含52個城市，每一個城市i均由整數型座標來表示。第二個數據集爲CITY-10,是使用0-1均勻分佈隨機函數生成的10個浮點型城市數據。相關實驗參數設置如2部分所述。

    實驗中利用模擬退火算法分別計算在CITY-52和CITY-10上的TSP求解，獲得城市序列總距離於溫度迭代次數的曲線關係分別以下所示。

 

圖 1 城市總距離d與溫度迭代次數t的曲線（左圖爲CITY-52，右圖爲CITY-10）

    從圖中能夠看出，在CITY-52上，模擬退火算法僅通過25次左右溫度迭代便收斂到最優解，而在CITY-10上，因爲城市數據爲浮點型數據，且知足0-1均勻分佈，城市之間距離差別小，所以在搜索過程當中波動較大，可是在通過100屢次溫度迭代後，算法也收斂到了最優解。

實驗代表不管在大值整數型數據仍是在0-1浮點型數據下，模擬退火算法都可以在TSP問題上收斂，證實了算法的有效性和先進性。並且通過屢次測試，算法在CITY-52都可以收斂到全局最優解7544.3659,對應最優序列爲：[10 50 32 42 9 8 7 40 18 44 31 48 0 21 30 17 2 16 20 41 6 1 29 22 19 49 28 15 45 43 33 34 35 38 39 36 37 47 23 4 14 5 3 24 11 27 26 25 46 12 13 51]。

6.結論

    旅行商問題(TSP)可利用模擬退火算法、遺傳算法、蟻羣算法等多種算法求解，在本實驗中，咱們利用模擬退火算法來求解TSP問題。咱們經過實驗分析以及權衡時間和質量等因素肯定了模擬退火算法的初始算法，並採用2變換法來產生新解。實驗測試了算法在CITY-52以及CITY-10兩個數據集上的表現狀況，在少數次溫度迭代以後，算法便收斂到了全局最優解，證實了模擬退火算法在TSP問題上的有效性和先進性。

    實驗時爲了便於實現和觀察，將指標函數差設爲城市序列的總距離差，可是因爲2變換法的對稱性，指標函數差能夠簡化爲變換點距離之差，從而提升算法運行時的計算速度。此外，實驗的初始參數還能夠根據城市數據類型進行相應的量級變換，這些都可以進一步得提升模擬退火算法表現。

7.源碼

import numpy as np
import matplotlib.pyplot as plt


# 隨機生成nums個城市數據
def city_generator(nums):
    from numpy import random
    
    citys = []
    N = 0
    while N < nums:
        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,包含52個城市 
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]])
    return citys
    

class TSP:
    
    def __init__(self):
        self._init_citys()
        self._init_dist()
        self._init_params()

    ## 初始化參數
    def _init_params(self):
        # 初始溫度
        self.T = 280
        # 每一個溫度下的迭代次數
        self.L = 100 * self.n
        # 溫度降低縮減因子
        self.alpha = 0.92
        # 初始狀態
        self.S = np.arange(self.n)

    # 初始化城市
    def _init_citys(self):
        self.citys = init_city52()
        # self.citys = city_generator(30)
        self.n = self.citys.shape[0]

    # 獲得距離矩陣的函數
    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):
                self.Dist[i][j] = self.Dist[j][i] = np.linalg.norm(self.citys[i] - self.citys[j])


    # 計算在狀態下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

    # 從狀態S的鄰域中隨機選擇
    def neighbors(self, S):
        S_neibor = 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)
        S_neibor[u], S_neibor[v] = S_neibor[v], S_neibor[u]
        return S_neibor
    
    # 從狀態S的鄰域中隨機選擇
    def neighbors2(self, S):
        S_neibor = np.copy(S)
        u = np.random.randint(0, self.n)
        v = np.random.randint(0, self.n)
        if u > v:
            u, v = v, u
        while u == v:
            v = np.random.randint(0, self.n)
        temp = S_neibor[u:v]
        S_neibor[u:v] = temp[::-1]
        return S_neibor

    # 模擬退火搜索過程
    def search(self):
        Ts = []
        Es = []
        while self.T >= 0.1:
            print('search on T:{}'.format(self.T))
            for i in range(self.L):
                E_pre = self.energy(self.S)
                S_now = self.neighbors2(self.S)
                E_now = self.energy(S_now)
                if (E_now < E_pre) or (np.exp((E_pre - E_now) / self.T) >= np.random.rand()):
                    self.S = S_now
                    
            Ts.append(self.T)
            E_now = self.energy(self.S)
            Es.append(E_now)
            print(E_now)
            
            # 判斷是否達到終止狀態
            self.T = self.T * self.alpha
        print(self.S)
        print('finished\n')
                
        return Ts, Es


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