模擬退火演算法Python程式設計（4）旅行商問題

youcans發表於2021-05-04

1、旅行商問題(Travelling salesman problem, TSP)

　　旅行商問題是經典的組合優化問題，要求找到遍歷所有城市且每個城市只訪問一次的最短旅行路線，即對給定的正權完全圖求其總權重最小的Hamilton迴路：設有 n個城市和距離矩陣 D=[dij]，其中dij表示城市i到城市j的距離（i,j = 1,2 … n），則問題是要找出遍訪每個城市恰好一次的一條迴路並使其路徑長度為最短。旅行商問題屬於NP完全問題，其全域性優化解的計算量以問題規模的階乘關係增長。旅行商問題不僅作為一類典型的組合優化問題經常被用作演算法研究和比較的案例，許多實際應用問題如路徑規劃、交通物流、網路管理也可轉化為旅行商問題。
　　目前，旅行商問題的研究主要集中於探索和發展各種高效近似最優解的優化方法，包括基於問題特徵資訊（如城市位置、距離、角度等）構造的各種啟發式搜尋演算法，以及通過模擬或解釋自然規律而發展的模擬退火演算法、遺傳演算法、蟻群演算法、神經網路演算法等智慧優化演算法或將二者相結合的混合演算法。
　　模擬退火演算法不僅可以解決連續函式優化問題，KIRKPATRICK在1983年成功將其應用於求解組合優化問題。模擬退火演算法現已成為求解旅行商問題的常用方法，通常採用反序、移位和交換等操作運算元產生新解。

2、模擬退火演算法求解旅行商問題

　　模擬退火演算法要從當前解的鄰域中產生新的候選解，解的表達形式和鄰域結構對演算法收斂非常重要。組合優化問題不同於函式優化，其自變數不是連續變化的，目標函式不僅依賴於自變數的數值，而且與變數的排列次序有關。極端地，旅行商問題的路徑長度僅取決於排列次序，因此常用城市編號的序列來表示解。
　　新解的產生機制是在當前解序列的基礎上進行變換操作，隨機改變序列中某些點的排列次序，常見的基本變換操作有交換運算元(Swap Operator)、反序運算元（Inversion Operator）、移位運算元(Move Operator)等。
　　交換運算元將當前路徑 S_now 中的第 i 個城市 C_i 與第 j 個城市 C_j 的位置交換，得到新路徑 S_swap ：
　　S_now = C_1⋯C_(i-1)∙(C_i)∙C_(i+1)⋯C_(j-1)∙(C_j)∙C_(j+1)⋯C_n
　　S_swap= C_1⋯C_(i-1)∙(C_j)∙C_(i+1)⋯C_(j-1)∙(C_i)∙C_(j+1)⋯C_n
　　反序運算元也稱2-opt，將當前路徑 S_now 中從第 i 個城市 C_i 到第 j 個城市 C_j 之間的城市排列順序反向翻轉，得到新路徑 S_inv ：
　　S_now = C_1⋯C_(i-1)∙(C_i∙C_(i+1)⋯C_(j-1)∙C_j)∙C_(j+1)⋯C_n
　　S_inv = C_1⋯C_(i-1)∙(C_j∙C_(j-1)⋯C_(i+1)∙C_i)∙C_(j+1)⋯C_n
　　移位運算元相當於 Or-opt 操作t，將當前路徑 S_now 中的第 i 個城市 C_i 移動到第 j 個城市 C_j 之後的位置，得到新路徑 S_move ：
　　S_now = C_1⋯C_(i-1)∙(C_i)∙C_(i+1)⋯C_(j-1)∙C_j∙C_(j+1)⋯C_n
　　S_move= C_1⋯C_(i-1)∙C_(i+1)⋯C_(j-1)∙C_j∙(C_i)∙C_(j+1)⋯C_n

3、程式說明

　　下段給出了模擬退火演算法求解旅行商問題的 Python程式。對於程式中的一些細節處理，說明如下：

資料的獲取。為了避免讓讀者下載資料檔案，程式中採用直接賦值方式讀取旅行城市位置的座標。實際上，通常是使用資料檔案給出城市位置的引數，程式中已經給出了讀取 TSPLib（旅行商問題國際標準資料集）的子程式和呼叫方法。
城市間距的計算。按照 TSPLib 的處理規範，一是城市間距按歐式距離計算，不要說地球是圓的要算球面間距了；二是對計算的城市間距離取整（四捨五入），而不是用實數表示城市間距離，這會導致一些優化結果的差異。
新路徑的產生。為便於理解，本程式只給出了交換運算元(Swap Operator)的子程式。交換操作非常容易理解，但實際上優化效率和效果並不好，反序操作的效能顯著優於交換運算元。
路徑長度的計算。路徑長度計算，當然就是把每一段路徑相加，包括從最後一個城市回到起點的這段路程。這樣處理當然是正確的，但並不是最好的。對於 N 城市的TSP問題，這樣計算路徑長度每次都要將 N 段路程相加。而在模擬退火演算法中，每次通過變換操作產生新的路徑後，其實並不需要計算新路徑的長度，只需要計算新路徑與原有路徑的差 deltaE。因此我們可以直接計算變換操作導致的路徑差，對於交換運算元只要計算 8段距離的加減，對於反序操作只要計算 4段距離的加減，當 N很大的時候，這種處理方法可以極大地減小程式的計算量。本程式已經給出了該子程式 mutateSwapE 和呼叫方法，只要使用程式中相應的註釋語句即可。

4、模擬退火演算法求解旅行商問題 Python 程式

# 模擬退火演算法求解旅行商問題 Python程式
# Program: SimulatedAnnealing_v6.py
# Purpose: Simulated annealing algorithm for traveling salesman problem
# v1.0:
#   模擬退火求解旅行商問題（TSP）基本演算法
#   (1) 直接讀取城市座標資料（CTSP31）
#   (2) 僅採用交換(Swap)方式產生新解
#   (3) 圖形輸出：最優路徑圖，優化過程圖
#   (4) 城市間距離取整
# Copyright 2021 YouCans, XUPT
# Crated：2021-05-01

#  -*- coding: utf-8 -*-
import math                         # 匯入模組 math
import random                       # 匯入模組 random
import pandas as pd                 # 匯入模組 pandas 並簡寫成 pd
import numpy as np                  # 匯入模組 numpy 並簡寫成 np YouCans
import matplotlib.pyplot as plt     # 匯入模組 matplotlib.pyplot 並簡寫成 plt

np.set_printoptions(precision=4)
pd.set_option('display.max_rows', 20)
pd.set_option('expand_frame_repr', False)
pd.options.display.float_format = '{:,.2f}'.format


# 子程式：初始化模擬退火演算法的控制引數
def initParameter():
    # custom function initParameter():
    # Initial parameter for simulated annealing algorithm
    tInitial = 100.0                # 設定初始退火溫度(initial temperature)
    tFinal  = 1                     # 設定終止退火溫度(stop temperature)
    nMarkov = 1000                # Markov鏈長度，也即內迴圈執行次數
    alfa    = 0.99                 # 設定降溫引數，T(k)=alfa*T(k-1)

    return tInitial,tFinal,alfa,nMarkov


# 子程式：讀取TSPLib資料
def read_TSPLib(fileName):
    # custom function read_TSPLib(fileName)
    # Read datafile *.dat from TSPlib
    # return coordinates of each city by YouCans, XUPT

    res = []
    with open(fileName, 'r') as fid:
        for item in fid:
            if len(item.strip())!=0:
                res.append(item.split())

    loadData = np.array(res).astype('int')      # 資料格式：i Xi Yi
    coordinates = loadData[:,1::]
    return coordinates


# 子程式：計算各城市間的距離，得到距離矩陣
def getDistMat(nCities, coordinates):
    # custom function getDistMat(nCities, coordinates):
    # computer distance between each 2 Cities
    distMat = np.zeros((nCities,nCities))       # 初始化距離矩陣
    for i in range(nCities):
        for j in range(i,nCities):
            # np.linalg.norm 求向量的範數（預設求 二範數），得到 i、j 間的距離
            distMat[i][j] = distMat[j][i] = round(np.linalg.norm(coordinates[i]-coordinates[j]))
    return distMat                              # 城市間距離取整（四捨五入）


# 子程式：計算 TSP 路徑長度
def calTourMileage(tourGiven, nCities, distMat):
    # custom function caltourMileage(nCities, tour, distMat):
    # to compute mileage of the given tour
    mileageTour = distMat[tourGiven[nCities-1], tourGiven[0]]   # dist((n-1),0)
    for i in range(nCities-1):                                  # dist(0,1),...dist((n-2)(n-1))
        mileageTour += distMat[tourGiven[i], tourGiven[i+1]]
    return round(mileageTour)                     # 路徑總長度取整（四捨五入）


# 子程式：繪製 TSP 路徑圖
def plot_tour(tour, value, coordinates):
    # custom function plot_tour(tour, nCities, coordinates)

    num = len(tour)
    x0, y0 = coordinates[tour[num - 1]]
    x1, y1 = coordinates[tour[0]]
    plt.scatter(int(x0), int(y0), s=15, c='r')      # 繪製城市座標點 C(n-1)
    plt.plot([x1, x0], [y1, y0], c='b')             # 繪製旅行路徑 C(n-1)~C(0)
    for i in range(num - 1):
        x0, y0 = coordinates[tour[i]]
        x1, y1 = coordinates[tour[i + 1]]
        plt.scatter(int(x0), int(y0), s=15, c='r')  # 繪製城市座標點 C(i)
        plt.plot([x1, x0], [y1, y0], c='b')         # 繪製旅行路徑 C(i)~C(i+1)

    plt.xlabel("Total mileage of the tour:{:.1f}".format(value))
    plt.title("Optimization tour of TSP{:d}".format(num))  # 設定圖形標題
    plt.show()


# 子程式：交換操作運算元
def mutateSwap(tourGiven, nCities):
    # custom function mutateSwap(nCities, tourNow)
    # produce a mutation tour with 2-Swap operator
    # swap the position of two Cities in the given tour

    # 在 [0,n) 產生 2個不相等的隨機整數 i,j
    i = np.random.randint(nCities)          # 產生第一個 [0,n) 區間內的隨機整數
    while True:
        j = np.random.randint(nCities)      # 產生一個 [0,n) 區間內的隨機整數
        if i!=j: break                      # 保證 i, j 不相等

    tourSwap = tourGiven.copy()             # 將給定路徑複製給新路徑 tourSwap
    tourSwap[i],tourSwap[j] = tourGiven[j],tourGiven[i] # 交換 城市 i 和 j 的位置————簡潔的實現方法

    return tourSwap


# 子程式：交換操作運算元 —— 計算 deltaE
def mutateSwapE(tour, n, dist):
    # custom function mutateSwapE(tour, n, dist) by YouCans, XUPT
    # produce a mutation tour with 2-Swap operator
    # swap the position of two Cities in the given tour

    # 在 [0,n) 產生 2個不相等的隨機整數 i,j
    i = np.random.randint(1,n-1)            # 產生第一個 [1,n-1) 區間內的隨機整數
    while True:                             # [1,n-1) 是為了避免 0,(n-1)的特殊情況
        j = np.random.randint(1,n-1)        # 產生一個 [1,n-1) 區間內的隨機整數
        if i != j: break                    # 保證 i, j 不相等
    if i > j: i, j = j, i                   # 整理使 i < j (便於後續處理)

    tourSwap = tour.copy()                  # 將給定路徑複製給新路徑 tourSwap
    tourSwap[i],tourSwap[j] = tour[j],tour[i] # 交換 城市 i 和 j 的位置————簡潔的實現方法
    # 特別注意：深入理解深拷貝和淺拷貝的區別，注意記憶體中的變化，謹慎使用

    dESwap = dist[tour[i-1],tour[j]] - dist[tour[i-1],tour[i]]\
           + dist[tour[i+1],tour[j]] - dist[tour[i+1],tour[i]]\
           + dist[tour[j-1],tour[i]] - dist[tour[j-1],tour[j]]\
           + dist[tour[j+1],tour[i]] - dist[tour[j+1],tour[j]]
    # 特別注意： j=i+1 是特殊情況，不適用以上公式計算 dESwap
    if i+1==j:                            # 特殊處理：i,j相鄰時相當於INV
        dESwap = dist[tour[i-1],tour[j]] - dist[tour[i-1],tour[i]]\
               + dist[tour[j+1],tour[i]] - dist[tour[j+1], tour[j]]

    return tourSwap, dESwap


def main():
    # 主程式

    # # 讀取旅行城市位置的座標
    coordinates = 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]])
    # fileName = "../data/eil76.dat"                      # 資料檔案的地址和檔名
    # coordinates = read_TSPLib(fileName)                 # 呼叫子程式，讀取城市座標資料檔案

    # 模擬退火演算法引數設定
    tInitial,tFinal,alfa,nMarkov = initParameter()      # 呼叫子程式，獲得設定引數

    nCities = coordinates.shape[0]              # 根據輸入的城市座標 獲得城市數量 nCities
    distMat = getDistMat(nCities, coordinates)  # 呼叫子程式，計算城市間距離矩陣
    nMarkov = nCities                           # Markov鏈 的初值設定
    tNow    = tInitial                          # 初始化 當前溫度(current temperature)

    # 初始化準備
    tourNow   = np.arange(nCities)   # 產生初始路徑，返回一個初值為0、步長為1、長度為n 的排列
    valueNow  = calTourMileage(tourNow,nCities,distMat) # 計算當前路徑的總長度 valueNow
    tourBest  = tourNow.copy()                          # 初始化最優路徑，複製 tourNow
    valueBest = valueNow                                # 初始化最優路徑的總長度，複製 valueNow
    recordBest = []                                     # 初始化 最優路徑記錄表
    recordNow  = []                                     # 初始化 最優路徑記錄表

    # 開始模擬退火優化過程
    iter = 0                        # 外迴圈迭代次數計數器
    while tNow >= tFinal:           # 外迴圈，直到當前溫度達到終止溫度時結束
        # 在當前溫度下，進行充分次數(nMarkov)的狀態轉移以達到熱平衡

        for k in range(nMarkov):    # 內迴圈，迴圈次數為Markov鏈長度
            # 產生新解：
            tourNew = mutateSwap(tourNow, nCities)      # 通過 交換操作 產生新路徑
            # tourNew,deltaE = mutateSwapE(tourNow,nCities,distMat)   # 通過 交換操作 產生新路徑（計算 deltaE）
            valueNew = calTourMileage(tourNew,nCities,distMat) # 計算當前路徑的總長度
            deltaE = valueNew - valueNow

            # 接受判別：按照 Metropolis 準則決定是否接受新解
            if deltaE < 0:                          # 更優解：如果新解的目標函式好於當前解，則接受新解
                accept = True
                if valueNew < valueBest:            # 如果新解的目標函式好於最優解，則將新解儲存為最優解
                    tourBest[:] = tourNew[:]
                    valueBest = valueNew
            else:                                   # 容忍解：如果新解的目標函式比當前解差，則以一定概率接受新解
                pAccept = math.exp(-deltaE/tNow)    # 計算容忍解的狀態遷移概率
                if pAccept > random.random():
                    accept = True
                else:
                    accept = False

            # 儲存新解
            if accept == True:                      # 如果接受新解，則將新解儲存為當前解
                tourNow[:] = tourNew[:]
                valueNow = valueNew

        # 平移當前路徑，以解決變換操作避開 0,（n-1）所帶來的問題，並未實質性改變當前路徑。
        tourNow = np.roll(tourNow,2)                # 迴圈移位函式，沿指定軸滾動陣列元素

        # 完成當前溫度的搜尋，儲存資料和輸出
        recordBest.append(valueBest)                # 將本次溫度下的最優路徑長度追加到 最優路徑記錄表
        recordNow.append(valueNow)                  # 將當前路徑長度追加到 當前路徑記錄表
        print('i:{}, t(i):{:.2f}, valueNow:{:.1f}, valueBest:{:.1f}'.format(iter,tNow,valueNow,valueBest))

        # 緩慢降溫至新的溫度，
        iter = iter + 1
        tNow = tNow * alfa                              # 指數降溫曲線：T(k)=alfa*T(k-1)

    # 結束模擬退火過程

    # 圖形化顯示優化結果
    figure1 = plt.figure()     # 建立圖形視窗 1
    plot_tour(tourBest, valueBest, coordinates)
    figure2 = plt.figure()     # 建立圖形視窗 2
    plt.title("Optimization result of TSP{:d}".format(nCities)) # 設定圖形標題
    plt.plot(np.array(recordBest),'b-', label='Best')           # 繪製 recordBest曲線
    plt.plot(np.array(recordNow),'g-', label='Now')             # 繪製 recordNow曲線
    plt.xlabel("iter")                                          # 設定 x軸標註
    plt.ylabel("mileage of tour")                               # 設定 y軸標註
    plt.legend()                                                # 顯示圖例
    plt.show()

    print("Tour verification successful!")
    print("Best tour: \n", tourBest)
    print("Best value: {:.1f}".format(valueBest))

    exit()

if __name__ == '__main__':
    main()

5、執行結果

程式的執行結果只供參考，顯然這並不是最優結果。

Tour verification successful!
Best tour: 
 [18  7 40  2 16 17 31 38 39 36 24 27 26 11 50  3  5  4 23 47 37 14 42  9
  8 32 10 51 13 12 25 46 28 29  1  6 41 20 30 21  0 48 35 34 33 43 45 15
 49 19 22 44]
Best value: 9544.0

模擬退火演算法Python程式設計（3）整數規劃問題
2021-05-02
演算法Python程式設計
模擬退火演算法Python程式設計（2）約束條件的處理
2021-05-02
演算法Python程式設計
模擬退火演算法（1）Python 實現
2021-05-01
演算法Python
模擬退火演算法解析
2021-01-01
演算法
c++實現的模擬退火演算法
2018-08-27
C++演算法
模擬退火與爬山法
2024-10-12
解析·玄學模擬退火
2018-10-23
帶約束條件的運籌規劃問題求解（模擬退火演算法實現）
2023-04-18
演算法
使用模擬退火演算法優化 Hash 函式
2020-10-03
演算法優化函式
模擬退火學習筆記
2024-07-28
筆記
模擬退火學習筆記
2023-01-04
筆記
自己用C語言寫的一個模擬退火演算法
2018-08-27
C語言演算法
Python程式設計常見問題與解答
2019-02-03
Python程式設計
程式碼設計問題
2020-07-07
BZOJ2428 [HAOI2006]均分資料模擬退火
2019-01-20
Python 計算機二級模擬題賞析
2019-04-03
Python計算機
基於免疫演算法的TSP問題求解matlab模擬
2024-11-10
演算法Matlab
遺傳演算法解決旅行商問題(TSP)
2023-01-02
演算法
好程式設計師Python教程分享Python常見面試問題
2019-12-04
程式設計師Python面試
3n+1問題的考察及python程式設計
2020-11-12
Python程式設計
Python爬蟲程式設計常見問題解決方法
2018-09-07
Python爬蟲程式設計
一個51程式設計和模擬——流水燈
2020-10-03
程式設計
COST231-WI模型通道模擬,原始碼模擬matlab程式設計原始碼
2020-12-24
模型原始碼Matlab程式設計
Python退火演算法在高次方程的應用
2020-05-07
Python演算法
[20181031]模擬網路問題.txt
2018-10-31
BNDS 2024/4/6模擬賽題解
2024-04-06
Python教程之網路程式設計及前端常見問題!
2021-05-19
Python程式設計前端
藍橋杯模擬題——長草問題
2020-10-04
Android studio虛擬模擬器安裝問題
2021-01-03
Android
旅行商問題（TSP）概述
2024-10-12
Python中常用網路程式設計模組
2024-03-19
Python程式設計
Python基礎程式設計(十六)——函式4
2020-12-17
Python程式設計函式
noip模擬4
2024-11-03
Python程式設計基礎練習——撲克牌發牌問題
2020-11-11
Python程式設計
Proteus實現簡單51程式的設計與模擬
2020-10-08
用 150 行 Python 程式碼寫的量子計算模擬器
2019-02-16
Python
圖形程式設計問題記錄
2024-03-19
程式設計
SciTech-EECS-電設計- PCB設計-電路設計與模擬系統 + SPICE 模擬描述與模型
2024-07-17
模型