1 问题描述

旅行商问题(Travelling Salesman Problem, TSP)是这样一个问题:给定一系列城市和每对城市之间的距离,求解访问每一座城市一次并回到起始城市的最短回路。它是组合优化中的一个NP困难问题。本文选择公开的a280数据集进行测试,求解访问全部城市的最短路径。数据预处理及数据散点图如下:

import math,random,time,copy

# 城市坐标
cities = []
# 邻接矩阵,distance[i][j]表示cities[i]到cities[j]的距离
distance = []

def init():
with open("a280.tsp",'r') as f:
lines = f.readlines()
i=0
while "NODE_COORD_SECTION" not in lines[i]:
i += 1
for j in range(280):
c = lines[i+1+j].split()
cities.append([int(c[1]),int(c[2])])
for i in range(len(cities)):
distance.append([])
for j in range(len(cities)):
if i<=j:
distance[i].append(math.sqrt((cities[i][0]-cities[j][0])**2+(cities[i][1]-cities[j][1])**2))
else:
distance[i].append(distance[j][i])

# 计算一个解的路径总长度
def totalDistance(arr):
res = 0
for i in range(len(arr)-1):
res += distance[arr[i]][arr[i+1]]
res += distance[arr[-1]][arr[0]]
return res

在最优解搜索问题中,首先需要明确的是解的形式是什么以及评价解好坏的标准是什么。在旅行商问题中,解的形式是遍历城市的顺序,也就是城市数组排列的一种顺序,评价解的标准是路径总长度。

2 精确算法

精确算法得到的结果是全局最优解,因而在求解NP完全问题时必然会遇到时间复杂度爆炸的瓶颈。

2.1 暴力搜索

暴力搜索是采用深度优先遍历,依次遍历所有可能的路径,保证最终能得到全局最优解。算法结构简单,占用空间小,但在问题包含 $N$ 个节点的情况下总路径数高达 $N!$ ,时间复杂度也是 $O(N!)$ 。

2.2 动态规划

首先要证明动态规划的可行性,即问题的最优解所包含的子问题的解也是最优的。在本题中就是要证明,当 $s_0\rightarrow s_1\rightarrow s_2 … \rightarrow s_{n-1} \rightarrow s_0$ 是从 $s_0$ 出发经过 $n$ 个城市回到 $s_0$ 的最短路径时,子问题 $s_1\rightarrow s_2…\rightarrow s_{n-1} \rightarrow s_0$ 就是从 $s_1$ 出发经过 $n-1$ 个城市到达 $s_0$ 的最短路径,使用反证法就能证明。

设 $dp(i,V)$ 表示从顶点 $s_i$ 出发经过 $V$ 中各个顶点一次且仅一次,最后回到出发点 $s_0$ 的最短路径长度。当 $V$ 是空集时,直接回到 $s_0$ 即可,路径长度就是 $s_i$ 到 $s_0$ 的距离,记为 $dp(i,\emptyset)=c_{i0}$ ,当 $V$ 不是空集时,就要遍历当前可选择的所有子问题,即 $dp(i,V)=min{c_{ik}+dp(k,V-{k})}$

时间复杂度 $O(2^nn^2)$ ,空间复杂度 $O(n2^n)$

3 个体启发式算法

3.1 爬山法

爬山法是一种简单的局部择优的算法。算法从一个随机位置开始(也就是先随机生成一个解),在每次迭代中随机访问当前最优解附近的一个解(也就是对当前最优解做一些微小的改变),如果新的解更优就更新为当前最优解,如果新的解更差就丢弃。直到在当前最优解附近找不到更优解,说明算法找到了一个局部最优解,就把这个解当做最终结果返回。在实际应用时,用这种随机方法寻找局部最优的时间可能很长,为了缩短程序的运行时间,需要设定最大迭代次数,当算法已经达到最大迭代次数,即使还未找到局部最优解也要停止迭代,返回当前最优解。Python3实现如下,迭代次数设为十万次:

def hillClimb():
startTime = time.time()
# 初始的随机解
arr = list(range(len(cities)))
random.shuffle(arr)
maxIter = 100000
currIter = 0
minDistance = -1
history = []
while currIter<maxIter:
if currIter==1:
history.append([currIter,minDistance])
if currIter%1000==0 and currIter>0:
history.append([currIter,minDistance])
# 每次对解的改变就是随机交换两个位置的值
change = random.sample(range(len(cities)),2)
arr[change[0]],arr[change[1]] = arr[change[1]],arr[change[0]]
newDistance = totalDistance(arr)
if newDistance<minDistance or minDistance<0:
minDistance = newDistance
else:
# 不是更优解就恢复原样
arr[change[0]],arr[change[1]] = arr[change[1]],arr[change[0]]
currIter += 1
history.append([maxIter,minDistance])
endTime = time.time()
print("最短路程:"+str(minDistance))
print("运行时间:"+str(endTime-startTime))
print("最优路线:"+str([list(cities[i]) for i in arr]))
print("迭代历史:"+str(history))

最优结果对应的路线如下:

最优解搜索过程如下:

3.2 模拟退火

模拟退火算法(Simulated Annealing,SA)借鉴了物理中固体物质的退火过程,将固体加热到足够高的温度,使分子呈随机排列状态,然后逐步降温使之冷却,最后分子以低能状态排列,固体达到某种稳定状态。

模拟退火算法首先初始化一个随机解、一个较高的初始温度和一个降温速率,在每次迭代中对当前最优解做微小改变,如果新的解是更优解,就更新为当前最优解,与爬山法不同的是,如果新的解不是更优解,则以一个与当前温度呈正相关的概率接受新的解,一般根据Metropolis接受准则把概率设为 $p=exp(-\Delta F/T), \Delta F>0$ ,其中 $\Delta F$ 是两个解在评估函数下的差值, $T$ 是当前温度。

以一定概率接受更差的解让算法有机会逃离局部最优,理论上能比爬山法找到更好的解。而随着温度下降,接受概率越来越低,最终趋近于0,表示不再接受差解,算法将收敛到一个更好的局部最优解。Python3实现如下,迭代次数设为十万次:

def SA():
startTime = time.time()
# 初始的随机解
arr = list(range(len(cities)))
random.shuffle(arr)
maxIter = 100000
currIter = 0
minDistance = -1
T = 10000
alpha = 0.99
history = []
while currIter<maxIter:
if currIter==1:
history.append([currIter,minDistance])
if currIter%1000==0 and currIter>0:
history.append([currIter,minDistance])
# 每次对解的改变就是随机交换两个位置的值
change = random.sample(range(len(cities)),2)
arr[change[0]],arr[change[1]] = arr[change[1]],arr[change[0]]
newDistance = totalDistance(arr)
if newDistance<minDistance or minDistance<0:
minDistance = newDistance
# 一定概率接受差解
elif random.random()<math.exp(-(newDistance-minDistance)/T):
minDistance = newDistance
else:
arr[change[0]],arr[change[1]] = arr[change[1]],arr[change[0]]
currIter += 1
T *= alpha
history.append([maxIter,minDistance])
endTime = time.time()
print("最短路程:"+str(minDistance))
print("运行时间:"+str(endTime-startTime))
print("最优路线:"+str([list(cities[i]) for i in arr]))
print("迭代历史:"+str(history))

最优结果对应的路线如下:

最优解搜索过程如下:

3.3 禁忌搜索

禁忌搜索(Tabu Search或Taboo Search,简称TS)是对局部搜索(LS)的一种扩展,是一种全局寻优算法,其特点是采用禁忌技术,即用一个禁忌表记录下已经到达过的局部最优点及求解过程,在下一次搜索中,利用禁忌表中的信息不再或有选择地搜索这些点,以此来跳出局部最优点。该算法可以克服爬山算法全局搜索能力不强的弱点。

算法流程如下,

arg3-0

Python3实现如下,迭代次数设为四万次:

def TS():
startTime = time.time()
maxIter = 40000
currIter = 0
searchNum = 200
tabuLen = 20
tabuTable = []
history = []
# 初始的随机解
arr = list(range(len(cities)))
random.shuffle(arr)
minDistance = totalDistance(arr)
localBestPath = []
while currIter<maxIter:
if currIter==1:
history.append([currIter,minDistance])
if currIter%400==0 and currIter>0:
history.append([currIter,minDistance])
i = 0
# 当前邻域搜索
while i<searchNum:
tmp = copy.deepcopy(arr)
change = random.sample(range(len(cities)),2)
tmp[change[0]],tmp[change[1]] = tmp[change[1]],tmp[change[0]]
# 搜索不在禁忌表中的解
if tmp not in tabuTable:
newDistance = totalDistance(tmp)
if newDistance<minDistance:
localBestPath = copy.deepcopy(tmp)
minDistance = newDistance
i += 1
arr = copy.deepcopy(localBestPath)
# 更新禁忌表
if len(tabuTable)==tabuLen:
tabuTable = tabuTable[1:]
tabuTable.append(localBestPath)
currIter += 1
history.append([maxIter,minDistance])
endTime = time.time()
print("最短路程:"+str(minDistance))
print("运行时间:"+str(endTime-startTime))
print("最优路线:"+str([list(cities[i]) for i in arr]))
print("迭代历史:"+str(history))

最优结果对应的路线如下:

最优解搜索过程如下:

4 群体启发式算法

4.1 遗传算法

遗传算法(GeneticAlgorithm ,GA )借鉴了生物进化过程来优化随机搜索策略,与上述两种算法不同的是,遗传算法维护并更新的不是单个解,而是用一组解模拟种群。在迭代中主要模拟选择、遗传、变异、免疫等过程实现种群中个体信息的交换,提升整个种群的质量。在本题中的主要步骤如下:

  • 1、初始化种群规模、迭代次数、突变概率三个参数
  • 2、随机生成一组解,代表初代种群。根据每个解的路径长度设定它的适应度,适应度表示个体对环境的适应程度,优胜劣汰,越是精英的个体适应度应该越高,所以路径长度越短的解适应度就越高。设 $dis[i]$ 表示种群中第 $i$ 个解的路径长度, $totalDis$ 表示种群中所有解的路径长度总和,那么第 $i$ 个解的适应度就是:
  • 3、父代选择。选择将要生成下一代的两个父代个体,同样为了保证优胜劣汰,要让适应度高的个体容易被选中。一般采用轮盘赌法选择个体,首先产生一个随机概率 $rand$ ,依次累种群个体的适应度,当出现 $rand < \sum_{k=0}^jP(k)$ 时,说明从第0个个体到第 $j$ 个个体的累积适应度超过了随机概率 $rand$ ,此时选择第 $j$ 个个体,完成一次选择过程。之所以叫轮盘赌法,是因为这个选择过程就好像以一个随机的力度转轮盘,结果从0转到 $j$ 停止,非常形象。
  • 4、两个父代个体杂交产生一个子代个体。方法不固定,只要是综合两个解的信息生成新的解就可以,根据方法的不同,既可以生成一个子代,也可以生成两个子代。本文采取的方法是,随机生成一个交叉点 $i$ ,两个父代中位于交叉点之前的城市顺序不变,而交叉点之后的城市顺序改成在对方中的顺序,交叉完毕后两个父代变成了两个新的子代。
  • 5、子代变异。对于杂交产生的每个子代,生成一个随机数,如果随机数小于预设的变异概率,就对该子代做微小改变,视为遗传变异。本文采取的变异方法是随机交换两个城市的位置。
  • 6、种群更新。保留子代中的精英个体,本文采取的方式是通过杂交变异生成两倍种群大小的子代,根据每个子代的路径长度排序,选择路径最短的一半作为新一代种群。

Python3实现如下,迭代次数设为四万次:

def GA():
startTime = time.time()
history = []
groupSize = 1000
# 突变概率
mutation = 0.01
maxIter = 40000
currIter = 0
fatherGroup = []
sonGroup = []
fitness = []
minDistance = -1
minPath = []
# 初始化种群
totalDis = 0
for i in range(groupSize):
arr = list(range(len(cities)))
random.shuffle(arr)
fatherGroup.append(arr)
d = totalDistance(arr)
fitness.append(d)
totalDis += d
# 计算适应度
fitness = [totalDis/i for i in fitness]
totalFit = sum(fitness)
fitness = [i/totalFit for i in fitness]
while currIter<maxIter:
if currIter==1:
history.append([currIter,minDistance])
if currIter%400==0 and currIter>0:
print("ga: "+str(currIter))
history.append([currIter,minDistance])
sonDis = []
# 杂交
for i in range(groupSize):
# 轮盘赌选择父代
f1 = -1
r = random.random()
p = 0
while r>p:
f1 += 1
p += fitness[f1]
f2 = f1
while f1==f2:
f2 = -1
r = random.random()
p = 0
while r>p:
f2 += 1
p += fitness[f2]
f1 = fatherGroup[f1]
f2 = fatherGroup[f2]
# 杂交生成两个子代
r = random.randint(0,len(cities)-2)
s1 = f1[:i]
s2 = f2[:i]
for k in f1:
if k not in s2:
s2.append(k)
for k in f2:
if k not in s1:
s1.append(k)
r = random.random()
if r<mutation:
# 突变就是随机交换两个位置的值
change = random.sample(range(len(cities)),2)
s1[change[0]],s1[change[1]] = s1[change[1]],s1[change[0]]
s2[change[0]],s2[change[1]] = s2[change[1]],s2[change[0]]
sonGroup.append(s1)
sonGroup.append(s2)
sonDis.append(totalDistance(s1))
sonDis.append(totalDistance(s2))
# 子代按路径长度排序
zipped=zip(sonGroup,sonDis)
sort_zipped = sorted(zipped,key=lambda x:(x[1],x[0]))
result = zip(*sort_zipped)
sonGroup, sonDis = [list(x) for x in result]
sonGroup = sonGroup[:groupSize]
sonDis = sonDis[:groupSize]
# 种群更新
minDistance = min(sonDis)
minPath = sonGroup[sonDis.index(minDistance)]
totalDis = sum(sonDis)
fitness = [totalDis/i for i in sonDis]
totalFit = sum(fitness)
fitness = [i/totalFit for i in fitness]
fatherGroup = sonGroup
sonGroup = []
currIter += 1
history.append([maxIter,minDistance])
endTime = time.time()
print("最短路程:"+str(minDistance))
print("运行时间:"+str(endTime-startTime))
print("最优路线:"+str([list(cities[i]) for i in minPath]))
print("迭代历史:"+str(history))

最优结果对应的路线如下:

最优解搜索过程如下:

4.2 蚁群算法

蚁群算法(Ant Colony Optimization ,ACO )的灵感来源于蚂蚁觅食。蚂蚁觅食是一种群体性行为,蚂蚁在寻找食物时,会在其经过的路径上释放一种信息素,并能够感知其他蚂蚁释放的信息素。信息素浓度的大小体现了路径的远近,因为蚂蚁会以较大的概率优先选择信息素浓度高的路径,并释放一定量的信息素,以增强该路径上的信息素浓度,这形成了一种正反馈,最终使得蚂蚁能够找到一条到食物源的最短路径。这种思路很适合解决旅行商问题,只需要把蚂蚁的行走路径表示为问题的解,每次迭代的过程就是把所有蚂蚁随机放置在不同的位置,根据路径的信息素浓度选择下一个访问的城市,当所有蚂蚁访问完所有的城市后,一个迭代结束。所以信息素的更新和路径的选择策略是蚁群算法的两大核心。

路径选择:设蚂蚁的数量为 $m$ ,城市的数量为 $n$ ,城市 $i$ 与城市 $j$ 之间的距离为 $d_{ij}$ , $t$ 时刻城市 $i$ 与城市 $j$ 连接路径上的信息素浓度是 $\tau_{ij}(t)$ ,初始时刻,所有路径上的信息素浓度相同,即 $\tau_{ij}(0)=\tau_0$ 。设 $P_{ij}^k(t)$ 是 $t$ 时刻蚂蚁 $k$ 从城市 $i$ 转移到城市 $j$ 的概率,其计算公式为,

其中, $\eta_{ij}(t)$ 为启发函数, $\eta_{ij}(t)=1/d_{ij}$ ,表示蚂蚁从城市 $i$ 转移到城市 $j$ 的期望程度。 $allow_k$ 是蚂蚁 $k$ 待访问城市的集合,分段公式的第二段就表示蚂蚁转移到已经访问过的城市的概率是0。蚂蚁对下一个目的地的选择取决于信息素浓度和启发函数,而这里用的启发函数是距离的倒数,所以可以说蚂蚁的选择兼顾了群体行为留下的信息和当前城市到其他城市的客观距离,两个启发因子 $\alpha$ 和 $\beta$ 就是用来调节这两方面在决策中的重要程度。

信息素更新:当蚂蚁经过一条路径时,会释放信息素,同时路径上原有的信息素会有一定程度的挥发,设参数 $\rho(0<\rho<1)$ 表示信息素的挥发率,当所有蚂蚁完成一次循环后,各个路径上信息素的更新公式为,

其中, $\tau_{ij}^k$ 表示第 $k$ 只蚂蚁在城市 $i$ 到城市 $j$ 之间的路径上释放的信息素浓度, $\tau_{ij}$ 就表示所有蚂蚁在城市 $i$ 到城市 $j$ 之间的路径上释放的信息素浓度之和。对于 $\tau_{ij}^k$ 的计算方法,有三种模型:

  • 1、ant cycle system 模型其中 $Q$ 为常数,表示蚂蚁循环一次所释放的信息素总量, $L_K$ 表示第 $k$ 只蚂蚁经过路径的长度
  • 2、ant quantity system 模型
  • 3、ant density system 模型在上述三种模型中,一般采用ant cycle system 模型,让信息素的浓度随蚂蚁迄今走过的总长度的增加而下降,使得路径上的信息素浓度更包含一种全局性的信息。

Python3实现如下,迭代次数设为一万次:

def ACO():
startTime = time.time()
history = []
maxIter = 500
currIter = 0
minDistance = -1
minPath = []
# 初始化参数
antNum = len(cities)
alpha = 1
beta = 5
rho = 0.1
Q = 1000
# 初始化信息素矩阵
pheromone = [[1 for i in range(len(cities))] for j in range(len(cities))]
while currIter<maxIter:
if currIter==1:
history.append([currIter,minDistance])
if currIter%5==0 and currIter>0:
print("aco: "+str(currIter)+" "+str(minDistance))
history.append([currIter,minDistance])
# 初始化行走的路径
path = [[] for j in range(antNum)]
visit = [[0 for i in range(len(cities))] for j in range(antNum)]
antDis = []
# 随机放置蚂蚁
arr = list(range(len(cities)))
random.shuffle(arr)
for i in range(antNum):
path[i].append(arr[i])
visit[i][arr[i]] = 1
# 所有蚂蚁都走完一条完整路径时才算一次迭代
for i in range(antNum):
while len(path[i])<len(cities):
# 路径选择概率表
weight = [0 for k in range(len(cities))]
for k in range(len(cities)):
if visit[i][k]==0:
weight[k] = (pheromone[path[i][-1]][k]**alpha)/(distance[path[i][-1]][k]**beta+0.0001)
# 轮盘赌法选下个城市
totalWeight = sum(weight)
p = [w/totalWeight for w in weight]
r = random.random()
sumP = 0
nextCity = -1
while r>sumP:
nextCity += 1
sumP += p[nextCity]
path[i].append(nextCity)
visit[i][nextCity] = 1
for i in path:
antDis.append(totalDistance(i))
minDistance = min(antDis)
minPath = path[antDis.index(minDistance)]
# 更新信息素
delta = [[0 for i in range(len(cities))] for j in range(len(cities))]
for i in range(antNum):
tau = Q/antDis[i]
for j in range(1,len(cities)):
fromCity = path[i][j-1]
toCity = path[i][j]
delta[fromCity][toCity] += tau
delta[toCity][fromCity] += tau
for i in range(len(cities)):
for j in range(len(cities)):
pheromone[i][j] = pheromone[i][j]*(1-rho)+delta[i][j]
currIter += 1
history.append([maxIter,minDistance])
endTime = time.time()
print("最短路程:"+str(minDistance))
print("运行时间:"+str(endTime-startTime))
print("最优路线:"+str([list(cities[i]) for i in minPath]))
print("迭代历史:"+str(history))

最优结果对应的路线如下:

最优解搜索过程如下,在几种算法里是效果最好的。因为蚁群算法的原始思想就是找最短路径,从初始化开始就是找最短路径,而不是像其他算法一样要修改成适合旅行商问题的形式,从一个随机初始化的解开始迭代。

4.3 粒子群算法

粒子群算法( Particle Swarm Optimization, PSO)源于对鸟群觅食行为的研究。一群鸟在随机搜索食物,在这个区域里只有一块食物。所有的鸟都不知道食物在那里。但是他们知道当前的位置离食物还有多远。那么找到食物的最优策略是什么呢。最简单有效的就是搜寻离食物最近的鸟的周围区域。用粒子模拟单个鸟类个体,粒子的位置对应优化问题的一个候选解,粒子的飞行过程即为该个体的解空间搜索过程。粒子的飞行速度可根据粒子历史最优位置和种群历史最优位置进行动态调整。粒子仅具有两个属性:速度和位置,速度代表移动的快慢,位置代表移动的方向。每个粒子单独搜寻的最优解叫做个体极值,粒子群中最优的个体极值作为当前全局极值。不断迭代,更新速度和位置,最终得到满足终止条件的最优解。

更新规则:粒子群算法的核心就是根据两个极值更新粒子的速度和位置。首先速度的更新公式为,

其中, $v_i$ 是粒子 $i$ 的速度, $x_i$ 是粒子 $i$ 的位置, $rand()$ 是介于 $(0,1)$ 的随机数, $pbest_i$ 是粒子 $i$ 的个体极值, $gbest$ 是粒子群的全局极值。组成该式子的三个部分分别叫做记忆项、自身认知项和群体认知项。记忆项代表上次的速度和位置,自身认知项代表了个体自身的历史经验,群体认知项代表了粒子间协同合作和知识共享的结果。 $\omega$ 、 $c1$ 和 $c2$
是三个权重因子,用来调节这三部分在粒子速度更新过程中的重要程度。最后,位置的更新公式相对简单, $x_i=x_i+v_i$

但是粒子群算法的原生公式不适用于旅行商问题,在解决旅行商问题时,要先定义粒子、位置、速度分别代表什么。一般思路是把粒子定义成问题的一个解,也就是城市序列,这个序列同时也代表粒子的位置。粒子移动的本质就是所有粒子向极值靠近,所以序列向序列靠近可以解释成序列片段的复制。这样看来,用粒子群算法解决旅行商问题就变得很像遗传算法的思路,只不过是把父代之间的序列片段交换,变成了粒子和极值之间的交换,三个权重因子也就代表着三次交换的片段的长度。

4.4 免疫算法

和遗传算法大同小异,都是优胜劣汰的种群更新思路。

4.5 鱼群算法

4.6 人工蜂群算法

4.7 蛙跳算法

4.8 烟花算法

4.9 萤火虫算法

4.10 细菌觅食算法