Python数模笔记—求解旅行商问题的联合算子模拟退火算法（完整例程）
文章目录
0 摘要
为了提高模拟退火算法求解旅行商问题的搜索空间和优化效率，分析现有反序、移位和交换等操作算子的特征与相互关系，发现交换操作等价于两个嵌套的反序操作的叠加复合。
基于这种等价关系提出一种新的交换-反序联合算子，可以不增大计算量而同时获得由交换操作和反序操作所产生的 3条新路径，择优作为新解，从而提高算法的优化性能。
以Eil51、Eil76、Eil101、Ch150等不同规模的TSP问题进行测试，仿真结果表明联合算子的性能优于现有的移位、交换、反序算子及组合方案。
1 引言
旅行商问题(Travelling salesman problem, TSP) 是经典的组合优化问题，要求找到遍历所有城市且每个城市只访问一次的最短旅行路线，即对给定的正权完全图求其总权重最小的Hamilton回路[1]。旅行商问题属于NP完全问题，其全局优化解的计算量以问题规模的阶乘关系增长。旅行商问题不仅作为一类典型的组合优化问题经常被用作算法研究和比较的案例，很多实际应用问题如路径规划[2-3]、交通物流[4]、网络管理[5]也可以转化为旅行商问题来解决。
旅行商问题目前的研究集中于探索和发展各种高效求近似最优解的优化方法，主要包括基于问题特征信息（如城市位置、距离、角度等）构造的各种启发式搜索算法[6-7]，以及通过模拟或解释自然规律而发展出的模拟退火算法、遗传算法、蚁群算法、神经网络算法等独立于问题的智能优化算法[8-110]或将二者相结合的混合算法[121]。
模拟退火算法不仅可以解决连续函数优化问题，文[1213]在1983年成功将其应用于求解组合优化问题。模拟退火算法现已成为求解旅行商问题的常用方法，通常采用反序、移位和交换等操作算子产生新解。本文通过分析这些变换操作的性能和关系，发现交换操作与反序操作在机理和算法实现中的内在联系，提出了一种新的交换/反序联合算子，并通过不同规模的TSP问题对新的联合算子进行实验研究。
2 模拟退火算法求解旅行商问题
2.1 模拟退火算法
模拟退火算法结构简单，由温度更新函数、状态产生函数、状态接受函数和内循环、外循环终止准则构成。温度更新函数是指退火温度缓慢降低的实现方案，也称冷却进度表；状态产生函数是指由当前解随机产生新的候选解的方法；状态接受函数是指接受候选解的机制，通常采用Metropolis准则。外循环是由冷却进度表控制的温度循环；内循环是在每一温度下循环迭代产生新解的次数，也称Markov链长度。
模拟退火算法的基本步骤如下：
（1） 随机产生初始解 X_0；
（2）在温度 T_k 下进行 L_k 次循环迭代：由当前解 X_old产生新的候选解X_new，并按 Metropolis 接受准则以一定的概率接受候选解 X_new作为当前解：
式中 E(X_old )、E(X_new )分别是当前解 X_old、新解 X_new的目标函数值。
（3）按照冷却进度表控制退火温度，从温度初值缓慢降低，直到达到温度终值结束。
2.2 多个新解的竞争机制
模拟退火算法不断在当前解的邻域随机产生一个新解，并以一定概率接受较差的新解。虽然理论上模拟退火算法在降温足够慢、终温足够低、Markov链足够大的条件下能实现全局最优，但实际运行时无法满足这些条件，因此常常也会陷入局部极值，还可能丢失优化过程中的最优解。
为了提高优化效率和性能，文[14]提出了多粒子模拟退火算法，每次产生新解时通过多次对当前解扰动得到多个新解，选取其中最优者作为候选解。文[15]在求解 TSP 问题时，每次通过对当前路径进行断裂、倒置、重组得到 6 个新路径，从中选择最优的一条作为新解。类似地，文[16]提出多条马尔可夫链并行的模拟退火算法，通过多线程产生多个解，然后从中找出最优解。
这类同时产生多个新解的竞争机制，可以大幅降低丢失最优解的可能性，也有利于提高优化性能，但需要付出增大计算量的代价。由于计算量大致按同时产生的新解个数的倍数增长，这与直接将 Markov链长度增大数倍相比并不一定更有利，因此并未得到广泛的应用。
2.3 求解旅行商问题的操作算子
模拟退火算法要从当前解的邻域中产生新的候选解，解的表达形式和邻域结构对算法收敛非常重要。组合优化问题不同于函数优化，其自变量不是连续变化的，目标函数不仅依赖于自变量的数值，而且与变量的排列次序有关。极端地，旅行商问题的路径长度仅取决于排列次序，因此常用城市编号的序列来表示解。
新解的产生机制是在当前解序列的基础上进行变换操作，随机改变序列中某些点的排列次序，常见的基本变换操作有交换算子(Swap Operator)、反序算子（Inversion Operator）、移位算子(Move Operator)等[1][17]。
交换算子将当前路径 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 联合操作算子
3.1 交换、反序和移位操作的关系
模拟退火算法由新解与当前解的目标函数差 ∆E 计算新解的接受概率。在旅行商问题中，使用交换算子、反序算子和移位算子产生新路径时，并不需要直接计算新路径的总长度 E(X_new )，可以通过不同变换操作的特征直接计算 ∆E，从而极大地减少计算量。
∆ E s w a p = [ d ( C ( i − 1 ) , C j ) + d ( C j , C ( i + 1 ) ) + d ( C ( j − 1 ) , C i ) + d ( C i , C ( j + 1 ) ) ] − [ d ( C ( i − 1 ) , C i ) + d ( C i , C ( i + 1 ) ) + d ( C ( j − 1 ) , C j ) + d ( C j , C ( j + 1 ) ) ] ) （ 3 ） ∆ E i n v = [ d ( C ( i − 1 ) , C j ) + d ( C i , C ( j + 1 ) ) ] − [ d ( C ( i − 1 ) , C i ) + d ( C j , C ( j + 1 ) ) ] （ 4 ） ∆ E m o v e = [ d ( C ( i − 1 ) , C ( i + 1 ) ) + d ( C j , C i ) + d ( C i , C ( j + 1 ) ) ] − [ d ( C ( i − 1 ) , C i ) + d ( C i , C ( i + 1 ) ) + d ( C j , C ( j + 1 ) ) ] ) （ 5 ） ∆E_{swap}=[d(C_{(i-1)},C_j )+d(C_j,C_{(i+1)} )+d(C_{(j-1)},C_i )+d(C_i,C_{(j+1)} )] - [d(C_{(i-1)},C_i )+d(C_i,C_{(i+1)} )+d(C_{(j-1)},C_j )+d(C_j,C_{(j+1)} )] ) （3）\\ ∆E_{inv} =[d(C_{(i-1)},C_j )+d(C_i,C_{(j+1)} )]-[d(C_{(i-1)},C_i )+d(C_j,C_{(j+1)} )] （4）\\ ∆E_{move}=[d(C_{(i-1)},C_{(i+1)} )+d(C_j,C_i )+d(C_i,C_{(j+1)} )] - [d(C_{(i-1)},C_i )+d(C_i,C_{(i+1)} )+d(C_j,C_{(j+1)} )] ) （5） ∆Eswap=[d(C(i−1),Cj)+d(Cj,C(i+1))+d(C(j−1),Ci)+d(Ci,C(j+1))]−[d(C(i−1),Ci)+d(Ci,C(i+1))+d(C(j−1),Cj)+d(Cj,C(j+1))])（3）∆Einv=[d(C(i−1),Cj)+d(Ci,C(j+1))]−[d(C(i−1),Ci)+d(Cj,C(j+1))]（4）∆Emove=[d(C(i−1),C(i+1))+d(Cj,Ci)+d(Ci,C(j+1))]−[d(C(i−1),Ci)+d(Ci,C(i+1))+d(Cj,C(j+1))])（5）
式（3）～（5）中的 d(C_i,C_j ) 表示城市 C_i 与 C_j之间的距离。
对比这三种操作所产生的路径差 ∆E，可以发现其中有很多相同的距离项。特别地，反序操作的路径差 ∆E_inv 中的全部 4 段距离，都出现在交换操作的路径差 ∆E_swap 中，并且都具有相同的正负号。因此，在计算交换操作的路径差 ∆E_swap 时，可以先计算 ∆E_inv 中的 4 段距离，再计算剩下的 4 段距离。这并未增大计算量，就在计算 ∆E_swap 的同时得到了 ∆E_inv 的值，相当于同时获得了 2 条新路径的 ∆E。
进一步地，考察交换操作的路径差 ∆E_swap 中除 ∆E_inv 之外的 4 段距离，可以发现这正好是另一个反序操作的路径差 ∆E_inv^‘：
∆E_inv^’ =[d(C_j,C_(i+1) )+d(C_(j-1),C_i )]-[d(C_i,C_(i+1) )+d(C_(j-1),C_j )] （6）
∆E_swap= ∆E_inv+ ∆E_inv^’ （7）
不难看出，交换操作可以分解为两步，先将当前路径 S_now 中从城市 C_i 到城市 C_j 之间的城市排列顺序反向翻转得到新路径 S_inv，再将路径 S_inv 中从城市 C_(j-1) 到城市 C_(i+1) 之间的城市排列顺序再次翻转，所得到新路径就是通过交换算子所得到的新路径 S_swap ：
S_inv = C_1⋯C_(i-1)∙C_j∙(C_(j-1)⋯C_(i+1))∙C_i∙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
也就是说，交换操作等价于两个嵌套的反序操作的叠加复合。
3.2 交换-反序联合算子
更有意义的是，式（6）路径差也可以解释理解为对于当前路径 S_now 的反序操作：
将当前路径 S_now 中从城市 C_(i+1) 到城市 C_(j-1) 之间的城市排列顺序反向翻转，得到新路径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
因此，可以先分别计算两个反序操作所产生的新路径 S_inv 和 S_inv^’ 的路径差 ∆E_inv、∆E_inv^’ ，再由（7）计算交换操作的路径差 ∆E_swap，其计算量与按式（3）直接计算 ∆E_swap 是相同的。这样，就可以同时获得 3 条新路径 S_inv、S_inv^’ 和 S_swap 的 ∆E，可以选择最优的一条作为新解。
将这种同时获得 3 条新路径 S_inv、S_inv^’ 和 S_swap 并择优作为新解的操作，称为交换-反序联合算子。根据之前的分析，这种产生多个新解的竞争机制，可以提高模拟退火算法的优化性能。而联合算子的优点还在于，是不需要增大计算量就可以即可获得多个新解。
4 实验及结果分析
4.1 实验方案
为检验交换-反序联合算子的优化性能，从 TSPLib 案例库中选取 4 个不同规模的旅行商问题 Eil51、Eil76、Eil101、Ch150 作为测试案例，其城市数量分别为 51、76、101、150，已知的最优路径长度分别为 426、538、629和6528。
本文算法实验环境为：CPU：Intel®Core™ i5-4460S@2.90GHz，内存：8 GB/ 1600 MHz，操作系统：Windows7 (64位/Service Pack 1)，采用 Python 3.8编程。
对比测试方案为：将交换-反序联合算子与通过其它操作算子产生新解的方案进行对比测试：（1）移位算子（Move），（2）交换算子（Swap），（3）反序算子（Inv），（4）组合移位算子/交换算子/反序算子（M/S/I）：每次随机使用一种算子产生一个新解，（5）交换-反序联合算子（Joint）：每次同时产生 3个新路径 S_inv、S_inv^’ 和 S_swap 并择优作为新解。
使用模拟退火算法的基本方案：控制温度按照 T_k=α*T_(k-1) 指数衰减，衰减系数取 α∈(0.9,1.0)；如式（2）按照 Metropolis 准则接受新解。
模拟退火算法的性能受到优化参数的影响，对每个测试案例分别采用不同的参数各进行 2 组实验。第一组参数为：初温 100，终温 1，衰减系数 0.95，共 90个温度状态，Markov链长度 1000；第一组参数为：初温 100，终温 1，衰减系数 0.98，共 228个温度状态，Markov链线性增大，平均链长 1000。模拟退火算法的优化结果有一定的随机性，对每组参数各做 10次重复实验，重复实验优化结果的最小值、平均值、最大值和标准差如表1～表4 所示。
4.2 实验结果的分析
4 个不同规模的旅行商问题 Eil51、Eil76、Eil101、Ch150、在不同的测试参数下的结果，使用交换-反序联合算子的最小值、平均值、最大值都是最好的，表明联合算子的性能显著优于现有的移位、交换、反序算子及其组合的方案。这与前文的分析一致，得益于联合算子与其它方案相比搜索了 3倍数量的路径。
所有方案在第二组的性能比第一组都有提高，但计算时间也都更长，且时间增大与迭代次数大致成正比。这也意味着如果通过增大迭代次数来改善性能，将会直接付出计算量增大的代价。各种方案在相同参数下运行时间随问题规模有所增大，但差异并不显著，这是因为算法中仅采用路径差 ∆E，不需要计算完整的路径长度，因此计算量与问题规模的关系不大。
在不同问题和不同参数条件下，联合算子的计算时间比其它方案略有增大，但并不显著。这与前文的分析一致，表明联合算子虽然增大了搜索数量、提高了优化性能，但基本上没有增加计算量。
在不同规模问题的测试中，联合算子提高优化性能的幅度不同。较小规模的问题，其它算子也已经获得较好的结果，性能提升的幅度就小；较大规模的问题复杂程度高，其它算子的优化结果较差，联合算子提高性能的幅度更大。
联合算子在Eil51问题找到了已知的最优路径，在Eil76问题很接近已知的最优路径，但对 Eil101、Ch150 问题的优化结果与已知的最优路径还有 2～3%的误差。这一方面说明大规模 TSP问题的复杂程度更高，需要更大的搜索量；另一方面是因为测试目的是对比不同操作算子的性能，而不是测试算子的最优能力。在调整联合算子的参数后，可以获得更好的优化路径。
4 结论
本文研究模拟退火算法求解旅行商问题产生新解的操作算子，分析这些变换操作之间的关系，发现交换操作等价于两个嵌套的反序操作的叠加复合。基于这种关系，本文提出了一种新的交换-反序联合算子，可以不增大计算量而同时获得由交换操作和反序操作所产生的 3条新路径，择优作为新解，从而提高算法的优化性能。通过4个不同规模的TSP问题的实验，验证了联合算子的性能优于现有的移位、交换、反序算子，也优于这些算子的组合方案。
交换-反序联合算子不仅可以用于模拟退火算法，也可以应用于其它进化算法，如遗传算法、粒子群算法中的变异操作以扩展搜索空间。
附录1：参考文献
附录2：源程序
# SATSP_demo_V1C_20201215.py
# 模拟退火算法求解旅行商问题 Python程序
# v2.0:
# (1) 引入新的 反序&交换联合算子(Inv&Swap)，同时计算invA、invB和Swap=invA+invB
# v1.0:
# 模拟退火求解旅行商问题（TSP）基本算法
# v1.0B:
# (1) 增加产生新解的方式：反序(Inv)，移位(Move)
# (2) 可以独立采用 反序/移动/交换操作，也可以按一定概率联合使用
# (3) 通过读取数据文件获得输入数据，如城市坐标
# (4) 计算运行时间
# v1.0C:
# (1) 计算新解与原解的路径差（变换操作的能量差）作为接受判据
# 不再计算新路径的里程（目标函数），极大地减小了计算量
# (2) 可选的不同降温曲线
# (4) 优化结果校验，检验最优路径是否准确
# Programmed by youcans@qq.com
# 日期：2020-12-16
# -*- coding: utf-8 -*-
import math # 导入模块 math
import random # 导入模块 random
import sys # 导入模块 sys
import time # 导入模块 time
import pandas as pd # 导入模块 pandas 并简写成 pd
import numpy as np # 导入模块 numpy 并简写成 np
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 initParameter2():
# custom function initParameter():
# Initial parameter for simulated annealing algorithm
tInitial = 100.0 # 设定初始退火温度(initial temperature)
tFinal = 1 # 设定终止退火温度(stop temperature)
maxMarkov = 2000 # Markov链长度，也即内循环运行次数
alfa = 0.985 # 设定降温参数，T(k)=alfa*T(k-1)
beta = 1.015 # Markov链增长曲线的参数
return tInitial,tFinal,alfa,beta,maxMarkov
# 子程序：读取TSPLib数据
def read_TSPLib(fileName):
# custom function read_TSPLib(fileName)
# Read datafile *.dat from TSPlib
# return coordinates of each city
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) # 路径总长度取整（四舍五入）
# 子程序：校验路径及路径长度的准确性
def verifyTour(tourGiven, valueGiven, nCities, distMat):
# custom function verifyTour(tourGiven, mileageGiven, nCities, distMat):
# to verify the given tour and its mileage
errorFlag = False # 初始化错误标识
sumTour = nCities*(nCities-1)/2 # 0,...n-1
if sum(tourGiven) != sumTour:
print("Error 1: Illegal tour!")
print(sumTour,sum(tourGiven))
errorFlag = True
valueTour = calTourMileage(tourGiven, nCities, distMat) # 计算给定路径的总长度
if abs(valueGiven - valueTour) > 1e-3: # 检验路径长度
print("Error 2: Wrong total millage!")
errorFlag = True
return errorFlag # 返回校验结果
# 子程序：绘制 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] = tourGiven[j] # 交换 城市 i 和 j 的位置
#tourSwap[j] = tourGiven[i]
# 交换 城市 i 和 j 的位置————简洁的实现方法
# 特别注意：深入理解深拷贝和浅拷贝的区别，注意内存中的变化，谨慎使用
tourSwap[i],tourSwap[j] = tourGiven[j],tourGiven[i]
return tourSwap
# 子程序：移动操作算子
def mutateMove(tourGiven, nCities):
# custom function mutateMove(nCities, tourNow)
# produce a mutation tour with Move operator
# Move a random tour fragment i~j and insert it behind the ity k
# 在 [0,n) 产生 3个不相等的随机整数 i,j,k
i = np.random.randint(nCities) # 产生第一个 [0,n) 区间内的随机整数
while True:
j = np.random.randint(nCities) # 产生一个 [0,n) 区间内的随机整数
k = np.random.randint(nCities) # 产生一个 [0,n) 区间内的随机整数
if (i!=j)&(i!=k)&(j!=k): # i,j,k 互不相等时跳出 while循环
break
if i>j: i,j = j,i # 整理使 i < j < k (便于后续处理)
if i>k: i,k = k,i
if j>k: j,k = k,j
# 将 [i, j) 区间的子路径片段插入到 k 之后
#tourMove = tourGiven.copy()
#fragment = tourGiven[i:j].copy() # 复制 [i,j) 子路径片段
#tourMove[i:i+k-j+1] = tourGiven[j:k+1].copy() # 平移 [j,k] 子路径片段
#tourMove[k-j+1+i:k+1] = fragment.copy() # 插入 [i,j) 子路径片段
# 将 [i, j) 区间的子路径片段插入到 k 之后————简洁的实现方法
frag1, frag2, frag3, frag4 = tourGiven[0:i], tourGiven[i:j], tourGiven[j:k+1], tourGiven[k+1:]
tourMove = np.concatenate([frag1, frag3, frag2, frag4])
return tourMove
# 子程序：反序操作算子
def mutateInv(tourGiven, nCities):
# custom function mutateInv(nCities, tourGiven)
# produce a mutation tour with Inverse operator(2-Opt)
# Reversing a sub-permutation between 2 random Cities
i = np.random.randint(nCities) # 产生第一个 [0,n) 区间内的随机整数
while True:
j = np.random.randint(nCities) # 产生一个 [0,n) 区间内的随机整数
if i!=j: break # 保证 i, j 不相等
if i>j: i,j = j,i # 整理使 i < j (便于后续处理)
tourInv = tourGiven.copy() # 将给定路径复制给新路径 tourInv
# 直接用 for 循环对 array 的位 进行赋值，是最初级但却是最安全的处理
for v in range(i,j+1):
tourInv[v] = tourGiven[j+i-v] # 将 [i:j] 逆序后复制给 新路径
# 将 [i:j] 逆序后复制给 新路径————简洁的实现方法
# 特别注意：深入理解深拷贝和浅拷贝的区别，注意内存中的变化，谨慎使用
#tourInv[i:j] = tourGiven[i:j][::-1]
return tourInv
# 子程序：反序操作算子 —— 计算 deltaE
def mutateInvE(tour, n, dist):
# custom function mutateInvE(tour, n, dist)
# produce a mutation tour with Inverse operator(2-Opt)
# Reversing a sub-permutation between 2 random Cities
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 (便于后续处理)
tourInv = tour.copy() # 将给定路径复制给新路径 tourInv
# 直接用 for 循环对 array 的位 进行赋值，是最初级但却是最安全的处理
for v in range(i, j+1): # 共 j-i+1 项
tourInv[v] = tour[j+i-v] # 将 [i:j] 逆序后复制给 新路径
dEInv = 0.0
if i==0: # i-1: n-1
dEInv += dist[tour[n-1],tour[j]] - dist[tour[n-1],tour[i]]
else:
dEInv += dist[tour[i-1],tour[j]] - dist[tour[i-1],tour[i]]
if j==n-1: # j+1: 0
dEInv += dist[tour[0],tour[i]] - dist[tour[0],tour[j]]
else:
dEInv += dist[tour[j+1],tour[i]] - dist[tour[j+1],tour[j]]
# 特别注意： i=0,j=n-1 是特殊情况，不适用以上公式计算 dEInv
return tourInv, dEInv
# 子程序：交换操作算子 —— 计算 deltaE
def mutateSwapE(tour, n, dist):
# custom function mutateSwapE(tour, n, dist)
# 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] = tourGiven[j] # 交换 城市 i 和 j 的位置
#tourSwap[j] = tourGiven[i]
# 交换 城市 i 和 j 的位置————简洁的实现方法
tourSwap[i],tourSwap[j] = tour[j],tour[i]
# 特别注意：深入理解深拷贝和浅拷贝的区别，注意内存中的变化，谨慎使用
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
# 子程序：移动操作算子 —— 计算 deltaE
def mutateMoveE(tour, n, dist):
# custom function mutateMove(nCities, tourNow)
# produce a mutation tour with Move operator
# Move a random tour fragment i~j and insert it behind the ity k
# 在 [0,n) 产生 3个不相等的随机整数 i,j,k
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) 区间内的随机整数
k = np.random.randint(1,n-1) # 产生一个 [1,n-1) 区间内的随机整数
if (i!=j)&(i!=k)&(j!=k): # i,j,k 互不相等时跳出 while循环
break
if i>j: i,j = j,i # 整理使 i < j < k (便于后续处理)
if i>k: i,k = k,i
if j>k: j,k = k,j
# 将 [i, j) 区间的子路径片段插入到 k 之后————简洁的实现方法
frag1, frag2, frag3, frag4 = tour[0:i], tour[i:j], tour[j:k+1], tour[k+1:]
tourMove = np.concatenate([frag1, frag3, frag2, frag4])
dEMove = dist[tour[i-1],tour[j]] - dist[tour[i-1],tour[i]]\
+ dist[tour[i], tour[k]] - dist[tour[j-1],tour[j]]\
+ dist[tour[k+1],tour[j-1]] - dist[tour[k+1],tour[k]]
return tourMove, dEMove
# 子程序：交换操作算子 —— 计算 deltaE
def mutateInvSwapE(tour, n, dist):
# custom function mutateInvSwapE(tour, n, dist):
# produce mutation routes with Inv operator
# Reversing a sub-permutation between 2 random cities
# InvA: Reverse from city i to city j (i,j included)
# InvA: Reverse from city i to city j (i,j included)
# InvA + InvB: Equivalent to Swap operation
# Output the best route among InvA, InvB and Swap.
# 随机整数 i,j 应满足：0<i<j<n-1, j-i>2
i = np.random.randint(3,n-4) # 产生 [3,n-4) 区间内的随机整数
gap = np.random.randint(2,n-1) # 产生 [2,n-1) 区间内的随机步长（间隔）
j = i + gap # 保证 j-i>2
if j>=(n-1): # 保证 j<n-1
j = i # 调换 i,j 保证 i>j
i = max(i+gap-n,1) # 保证 i>0
# 分段计算新旧路径差 S2,S3,S6,S7
S3 = dist[tour[j+1],tour[i]] - dist[tour[j+1],tour[j]]
S6 = dist[tour[i-1],tour[j]] - dist[tour[i-1],tour[i]]
S2 = dist[tour[j-1],tour[i]] - dist[tour[j-1],tour[j]]
S7 = dist[tour[i+1],tour[j]] - dist[tour[i+1],tour[i]]
# 计算 InvA, InvB 和 Swap 操作与给定路径的差 dE
dEInvA = S3 + S6 # dE_InvA (i,j Included)
dEInvB = S2 + S7 # dE_InvB (i,j Excluded)
dESwap = dEInvA + dEInvB # Swap is equivalent to InvA & InvB
tourNew = tour.copy() # 将给定路径复制给新路径 tourNew
dEmin = min(dEInvA,dEInvB,dESwap)
if dEInvA==dEmin: # --> InvA 最优
for v in range(i,j+1): # 共 j-i+1 项
tourNew[v] = tour[j+i-v] # 将 [i:j] 逆序后复制给 新路径
elif dEInvB==dEmin: # --> InvB 最优
for v in range(i+1,j): # 共 j-i-1 项
tourNew[v] = tour[j+i-v] # 将 [i+1:j-1] 逆序后复制给 新路径
else: # --> Swap 最优
tourNew[i], tourNew[j] = tour[j], tour[i]
return tourNew, dEmin
def main():
# 主程序
# 读取旅行城市位置的坐标
"""
coordinates = np.array([[1164.6,399.2],[1172.0,391.3],[1214.8,312.2],[1065.4,295.9],[911.1,299.7],
[876.8,437.7],[1062.7,384.7],[1116.5,408.2],[1083.3,228.4],[1266.3,457.5],
[1253.5,438.8],[1233.8,418.0],[1144.8,380.3],[1125.3,378.7],[1017.4,365.6],
[1170.0,366.5],[1136.0,347.6],[1187.8,320.4],[1172.7,318.6],[1201.9,302.6],
[1193.0,260.8],[1158.9,286.8],[1130.0,282.1],[1143.1,305.2],[1132.3,231.6],
[1103.5,200.2],[1037.3,360.3],[1089.5,342.7],[1040.6,306.7],[1067.1,265.7],
[1027.3,250.4]]) # CTSP31
"""
fileName = "../data/eil51.dat" # 数据文件的地址和文件名
coordinates = read_TSPLib(fileName) # 调用子程序，读取城市坐标数据文件
# 模拟退火算法参数设置
tInitial,tFinal,alfa,beta,maxMarkov = initParameter2() # 调用子程序，获得设置参数
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 # 外循环迭代次数计数器
timeStart = time.process_time() # 启动计时器，用于计算算法运行时间
while tNow >= tFinal: # 外循环，直到当前温度达到终止温度时结束
# 在当前温度下，进行充分次数(nMarkov)的状态转移以达到热平衡
for k in range(nMarkov): # 内循环，循环次数为Markov链长度
# 产生新解：通过在当前解附近随机扰动而产生新解，新解必须在 [min,max] 范围内
"""
turbRand = np.random.rand()
if turbRand<0.4:
tourNew,dE = mutateInvE(tourNow,nCities,distMat) # 通过 反序操作 产生新路径
elif turbRand<0.7:
tourNew,dE = mutateMoveE(tourNow,nCities,distMat) # 通过 移动操作 产生新路径
else:
tourNew,dE = mutateSwapE(tourNow,nCities,distMat) # 通过 交换操作 产生新路径
"""
#tourNew,dE = mutateInvE(tourNow,nCities,distMat) # 通过 反序操作 产生新路径
#tourNew,dE = mutateMoveE(tourNow,nCities,distMat) # 通过 移动操作 产生新路径
#tourNew,dE = mutateSwapE(tourNow,nCities,distMat) # 通过 交换操作 产生新路径
tourNew,dE = mutateInvSwapE(tourNow,nCities,distMat) # 通过 交换操作 产生新路径
valueNew = calTourMileage(tourNew,nCities,distMat) # 计算当前路径的总长度
deltaE = valueNew - valueNow
if abs(deltaE - dE) > 1e-3: # 检验路径长度
print("Error in deltaE verify!")
print(deltaE,dE)
sys.exit(0)
#deltaE = dE
#valueNew = valueNow + deltaE
# 接受判别：按照 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))
# 缓慢降温至新的温度，
#tNow = tNow * alfa # 指数降温曲线：T(k)=alfa*T(k-1)
#tNow = tInitial / (iter+1.0) # Cauchy 降温曲线：T(k)=T0/(k+1)
#tNow = tInitial / np.log(iter+1.0) # Boltzmann 降温曲线：T(k)=T0/log(k+1)
if iter<= 100: alfa = 0.97 # 按一定概率进行回火操作（重升温：升高当前温度）
elif iter<= 200: alfa = 0.98
elif iter<= 300: alfa = 0.985
else: alfa = 0.99
tNow = tInitial * (alfa**iter)
iter = iter + 1
#tNow = tNow * alfa # 指数降温曲线：T(k)=alfa*T(k-1)
nMarkov = min(round(nMarkov*beta), maxMarkov) # Markov链增长曲线：M(k)=beta*M(k-1)
# 结束模拟退火过程
# 模拟退火优化结果校验
errorFlag = verifyTour(tourBest,valueBest,nCities,distMat) # 校验 tourBest 路径及路径长度
if errorFlag:
print("Error in tour verify!" )
sys.exit(0)
else:
print("Tour verification successful!")
timeEnd = time.process_time()
timeRun = timeEnd - timeStart
print("CPU time: {:.2f} s".format(timeRun)) # 不计休眠时间
print("Best tour: \n", tourBest)
print("Best value: {:.1f}".format(valueBest))
# 图形化显示优化结果
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()
exit()
if __name__ == '__main__':
main()
附录3：运行结果
程序的运行结果只供参考，显然这并不是最优结果。
Tour verification successful!
CPU time: 32.11 s
Best tour:
[12 40 39 18 41 43 16 36 14 44 32 38 9 29 33 8 48 4 37 10 31 0 21 1
15 49 20 28 19 34 35 2 27 30 7 25 6 42 23 22 47 5 26 50 45 11 46 3
17 13 24]
Best value: 428.0
版权说明：
本文为论文投稿版本及附件资料，刊印全文请参见：知网论文：求解旅行商问题的联合算子模拟退火算法
Copyright 2021 YouCans, XUPT
Crated：2023-01-20
关注 Youcans，分享原创系列 https://blog.csdn.net/youcans
Python数模笔记-模拟退火算法（1）多变量函数优化
Python数模笔记-模拟退火算法（2）约束条件的处理
Python数模笔记-模拟退火算法（3）整数规划问题
Python数模笔记-模拟退火算法（4）旅行商问题