1、旅行商問題(Travelling salesman problem, TSP)
旅行商問題是經典的組合優化問題，要求找到遍歷所有城市且每個城市只訪問一次的最短旅行路線，即對給定的正權完全圖求其總權重最小的Hamilton回路：設有 n個城市和距離矩陣 D=[dij]，其中dij表示城市i到城市j的距離（i,j = 1,2 … n），則問題是要找出遍訪每個城市恰好一次的一條回路並使其路徑長度為最短。旅行商問題屬於NP完全問題，其全局優化解的計算量以問題規模的階乘關系增長。旅行商問題不僅作為一類典型的組合優化問題經常被用作算法研究和比較的案例，許多實際應用問題如路徑規划、交通物流、網絡管理也可轉化為旅行商問題。
目前，旅行商問題的研究主要集中於探索和發展各種高效近似最優解的優化方法，包括基於問題特征信息（如城市位置、距離、角度等）構造的各種啟發式搜索算法，以及通過模擬或解釋自然規律而發展的模擬退火算法、遺傳算法、蟻群算法、神經網絡算法等智能優化算法或將二者相結合的混合算法。
模擬退火算法不僅可以解決連續函數優化問題，KIRKPATRICK在1983年成功將其應用於求解組合優化問題。模擬退火算法現已成為求解旅行商問題的常用方法，通常采用反序、移位和交換等操作算子產生新解。
歡迎關注 Youcans 原創系列，每周更新數模筆記
Python數模筆記-PuLP庫
Python數模筆記-StatsModels統計回歸
Python數模筆記-Sklearn
Python數模筆記-NetworkX
Python數模筆記-模擬退火算法
2、模擬退火算法求解旅行商問題
模擬退火算法要從當前解的鄰域中產生新的候選解，解的表達形式和鄰域結構對算法收斂非常重要。組合優化問題不同於函數優化，其自變量不是連續變化的，目標函數不僅依賴於自變量的數值，而且與變量的排列次序有關。極端地，旅行商問題的路徑長度僅取決於排列次序，因此常用城市編號的序列來表示解。
新解的產生機制是在當前解序列的基礎上進行變換操作，隨機改變序列中某些點的排列次序，常見的基本變換操作有交換算子(Swap Operator)、反序算子（Inversion Operator）、移位算子(Move Operator)等。交換算子將當前路徑 S_now 中的第 i 個城市 Ci 與第 j 個城市 Cj 的位置交換；反序算子也稱2-opt，將當前路徑 S_now 中從第 i 個城市 Ci 到第 j 個城市 Cj 之間的城市排列順序反向翻轉；移位算子相當於 Or-opt 操作t，將當前路徑 S_now 中的第 i 個城市 Ci 移動到第 j 個城市 Cj 之后的位置。
3、 程序說明
下段給出了模擬退火算法求解旅行商問題的 Python程序。對於程序中的一些細節處理，說明如下：
- 數據的獲取。為了避免讓讀者下載數據文件，程序中采用直接賦值方式讀取旅行城市位置的坐標。實際上，通常是使用數據文件給出城市位置的參數，程序中已經給出了讀取 TSPLib（旅行商問題國際標准數據集）的子程序和調用方法。
- 城市間距的計算。按照 TSPLib 的處理規范，一是城市間距按歐式距離計算，不要說地球是圓的要算球面間距了；二是對計算的城市間距離取整（四舍五入），而不是用實數表示城市間距離，這會導致一些優化結果的差異。
- 新路徑的產生。為便於理解，本程序只給出了交換算子(Swap Operator)的子程序。交換操作非常容易理解，但實際上優化效率和效果並不好，反序操作的性能顯著優於交換算子。
4、模擬退火算法求解旅行商問題 Python 程序
# 模擬退火算法求解旅行商問題 Python程序
# === 關注 Youcans，分享更多原創系列 https://www.cnblogs.com/youcans/ ===
# Program: SimulatedAnnealing_v6.py
# Purpose: Simulated annealing algorithm for traveling salesman problem
# v1.0:
# 模擬退火求解旅行商問題（TSP）基本算法
# 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.98 # 設定降溫參數，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
# === 關注 Youcans，分享更多原創系列 https://www.cnblogs.com/youcans/ ===
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()
# === 關注 Youcans，分享更多原創系列 https://www.cnblogs.com/youcans/ ===
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
版權說明：
原創作品
Copyright 2021 YouCans, XUPT
Crated：2021-05-04
歡迎關注 Youcans 原創系列，每周更新數模筆記
Python數模筆記-PuLP庫（1）線性規划入門
Python數模筆記-PuLP庫（2）線性規划進階
Python數模筆記-PuLP庫（3）線性規划實例
Python數模筆記-Scipy庫（1）線性規划問題
Python數模筆記-StatsModels 統計回歸（1）簡介
Python數模筆記-StatsModels 統計回歸（2）線性回歸
Python數模筆記-StatsModels 統計回歸（3）模型數據的准備
Python數模筆記-StatsModels 統計回歸（4）可視化
Python數模筆記-Sklearn （1）介紹
Python數模筆記-Sklearn （2）聚類分析
Python數模筆記-Sklearn （3）主成分分析
Python數模筆記-Sklearn （4）線性回歸
Python數模筆記-Sklearn （5）支持向量機
Python數模筆記-NetworkX（1）圖的操作
Python數模筆記-NetworkX（2）最短路徑
Python數模筆記-NetworkX（3）條件最短路徑
Python數模筆記-模擬退火算法（1）多變量函數優化
Python數模筆記-模擬退火算法（2）約束條件的處理
Python數模筆記-模擬退火算法（3）整數規划問題
Python數模筆記-模擬退火算法（4）旅行商問題