前阵子帮一个开配送站的朋友优化派送路线,七个站点他手工排方案排了整整两天,我用程序跑了几分钟,给出来的总里程比他心算的结果短了一截。这件事让我重新审视了一个老朋友——TSP问题。作为我最早接触的NP-hard问题之一,它几乎贯穿了我后来所有路径规划类的项目:物流配送、园区巡查、电路板钻孔顺序、甚至外卖出餐顺序的排程,本质上都能抽象成同一套框架。
这篇文章把我在各种项目里实际用过的求解算法做一个完整盘点和复现,涵盖精确求解、贪心构造、局部搜索、模拟退火、遗传算法这些主流思路。文章不只是堆公式,每个算法我都会给出可以直接跑的Python实现,配上复杂度分析和真实调参经验,适合刚接触组合优化的同学,也适合想拿TSP练手、把多个算法横向对比的开发者。
1. TSP问题的本质:先看清楚我们到底在优化什么
1.1 从一个具体例子的数学抽象说起
TSP,全称Traveling Salesman Problem,旅行商问题。我第一次接触它是在课上学到的一段定义:给定若干城市坐标,求一条经过所有城市恰好一次、最后回到出发城市的最短回路。这个描述听起来平平无奇,但真正动手求解时你会发现,它像一个看着容易却死活摸不到顶层抽屉的柜子。
举个例子,你有五个城市,A到E。回路A→B→C→D→E→A是一条路径,A→C→B→D→E→A又是一条路径。城市的排列顺序不同,最终总路程就不同。任务就是从所有排列里找到总路程最小的那一个。推广到n个城市,候选回路有(n-1)!/2种(除以2是因为方向对称)。n=10的时候是18万种,n=15的时候已经超过430亿种。这个增长速度远远超过“普通电脑换个好CPU就能扛住”的范畴。
这类问题有一个专门的身份标签:NP-hard。我在实际项目里不太纠结这个理论名词,只需要记住一个现实结论——城市规模稍大,就不存在“暴力枚举扫一遍”的解法。你能做的是在“最优”和“足够好、算得快”之间做取舍,这也正是多种算法并存的根本原因。
1.2 对称TSP与非对称TSP:先确认你面对的是哪一类
TSP还有一个容易忽略的细分:对称TSP与非对称TSP。对称的意思是城市A到城市B的距离等于B到A的距离,这在欧氏平面坐标下天然成立;非对称则是单向距离,比如城市间的单向通行、上下坡油耗差异、订单序列切换成本不同等场景。
本文后面的代码全部以对称欧氏距离TSP为例,因为这是最常见、最容易可视化验证的版本。但你心里要清楚,一旦切换到非对称场景,很多贪心和局部搜索算法的优化条件需要重新推导,动态规划倒是只需要改一下距离函数即可。后面讲2-opt优化判据时,我会专门提醒这一点。
1.3 用基准评测说话:为什么随机生成城市坐标最靠谱
研究TSP算法,第一步不是急着写代码,而是先把评测环境搭好。我自己习惯的做法是:固定一个随机种子,生成一组城市坐标,比如20个、50个、100个城市各一组,坐标范围在0到1000的方形区域内。这样做的好处有两点:第一,实验可复现,每个算法跑出来的结果可以放进同一张表里比较;第二,随机均匀分布的TSP实例没有“让某个算法占便宜”的结构性特征,不会因为数据本身带规律而高估算法能力。
代码里提供生成器:
import random import math def generate_cities(n, seed=42): random.seed(seed) return [(random.uniform(0, 1000), random.uniform(0, 1000)) for _ in range(n)] def dist(a, b): # 欧氏距离 return math.hypot(a[0]-b[0], a[1]-b[1]) def calc_route_len(route, cities): n = len(route) return sum(dist(cities[route[i]], cities[route[(i+1) % n]]) for i in range(n))跑实验结果前先明确一个基准指标:对同一组城市,所有算法给出路径的“相对差距”,即(算法结果 - 最优结果) / 最优结果。最优结果如何得到?对小规模问题可以用下一章的动态规划精确求得,对大规模问题则只能用“当前已知的最好结果”作为近似参照。这也是学术界TSP benchmark的通用做法。
2. 精确解法的天花板:动态规划如何把暴力枚举变成可承受的指数
2.1 为什么暴力枚举第一个被淘汰
谈算法之前先看看最朴素的做法:枚举所有排列。
import itertools def brute_force(cities): n = len(cities) best_route = None best_len = float("inf") for perm in itertools.permutations(range(n)): # 固定起点为0,只看从0出发的排列 if perm[0] != 0: continue route = list(perm) length = calc_route_len(route, cities) if length < best_len: best_len = length best_route = route return best_route这个写法逻辑简单,但最多只能跑到10个城市左右,再多就需要等很长时间。n=12时排列数是39916800种,已经明显吃紧;n=15时就要跑十几分钟。原因在于排列组合的阶乘增长,指数级的动态规划在n=20时也只慢一档。很多初学者会在这个阶段误判“是不是我代码写得不够快”,其实不是实现的问题,是算法本身的复杂度决定了上限。
2.2 Held-Karp动态规划:用状态压缩记住子问题
动态规划的思路核心是“去重”。在暴力枚举过程中,同样一批城市、同样的起点终点,会被不同排列顺序反复访问无数次。与其每次重新计算,不如把这些中间状态存下来。
定义状态dp(mask, last):mask是一个bitmask,表示已经访问过的城市集合;last表示当前所在城市。dp的值是“从城市0出发,访问完mask集合里的所有城市,最后停在last的最短距离”。转移时,在未访问的城市中选一个作为下一个落脚点:
dp(mask, last) = min(dist(last, nxt) + dp(mask | (1 << nxt), nxt))
其中nxt遍历所有不在mask中的城市。边界条件是mask等于全集时,能从last直接返回出发城市0,距离为dist(last, 0)。
状态总数是2^n × n,每个状态转移最多n种选择,总复杂度O(n^2·2^n)。看起来仍是指数,但它把排列级别的n!降成了“近似2^n”,这在地狱难度里已经是实质性的进步。实际跑起来,n=20左右可以轻松求解,n=25就开始吃力,n=30基本要等很久。
下面是完整实现,为了能在求出最短长度的同时还原出具体路径,我用了一个重建技巧:在知道dp表后,从初始状态出发,每次挑一个能使“本段距离 + 剩余状态的最优距离”等于当前状态最优值的下一城市。
from functools import lru_cache def held_karp(cities): n = len(cities) full = (1 << n) - 1 def d(i, j): return dist(cities[i], cities[j]) @lru_cache(maxsize=None) def dp(mask, last): if mask == full: return d(last, 0) best = float("inf") for nxt in range(n): if not (mask >> nxt) & 1: best = min(best, d(last, nxt) + dp(mask | (1 << nxt), nxt)) return best # 根据dp表重建路径 mask, last = 1, 0 route = [0] while mask != full: for nxt in range(n): if (mask >> nxt) & 1: continue if abs(dp(mask, last) - (d(last, nxt) + dp(mask | (1 << nxt), nxt))) < 1e-9: route.append(nxt) mask |= (1 << nxt) last = nxt break return route2.3 动态规划在实战中的定位:验证算法正确性的“标尺”
动态规划在真实业务场景中很少直接用于大规模调度,因为城市数量一上百,指数复杂度直接让程序失去实用价值。但它的地位依然不可替代——它是你在开发其他启发式算法时的“裁判”。
我在做新算法效果评测时,会让动态规划与待测算法在同一组20城内实例上跑,看看贪心结果差多少、模拟退火能不能收敛到最优、2-opt改进后离最优还差多远。如果没有精确解做基准,一切“我的算法很不错”的自我评价都可能有幸存者偏差。
另外还有一个实际体验:动态规划的Python实现要非常小心递归深度和缓存。20个城市对应约2000万个状态,虽然lru_cache能挡住重复计算,但内存占用已经接近GB级别,建议n超过22时候用C++或改写成自底向上数组。本文专注讲解算法思想,为了便于阅读就用简洁的递归写法。
3. 贪心构造算法:最快的路线生成器,也是所有高级改动的种子
3.1 最近邻算法:从任意点出发,每次走最近的路
贪心算法的逻辑一句话就能说清:从某个城市出发,每次在所有未访问城市中选一个距离当前城市最近的,走过去,最后返回起点。
def nearest_neighbor(cities, start=0): n = len(cities) unvisited = set(range(n)) unvisited.remove(start) route = [start] cur = start while unvisited: nxt = min(unvisited, key=lambda j: dist(cities[cur], cities[j])) route.append(nxt) unvisited.remove(nxt) cur = nxt return route实现简单,时间复杂度只有O(n^2),n=10000都能秒跑完。但你别指望它交出最优解。在随机均匀分布的100城实例上,它通常比最优解长15%~30%,并且有一个明显特征:路径的后半段会产生很多“回头路”,因为前期只盯着局部最优,后期剩下的城市往往离当前点很远,被迫长距离折返。
这个现象背后的原因在组合优化里很典型:局部最优的简单叠加不等于全局最优。这也是为什么几乎所有TSP实战算法都会在贪心构造之后,接上一轮局部搜索去修复那些“前期短视造成的长距离尾巴”。
3.2 多次最近邻与多起点策略:低成本摊平坏运气
最近邻有一个很实用的变体:不指定起点,而是把每一个城市都当一次起点,跑n遍最近邻,保留最好的那一条。这样做能把“起点选择不当导致整条路径质量差”的随机性降到很低,复杂度变为O(n^3),但n在几千以内完全可接受。
我实测过100城实例,单起点最近邻平均与最优相差20%左右,多起点最近邻通常能把差距收窄到12%~15%。下面的代码展示了这个思路:
def repeated_nearest_neighbor(cities): best_route, best_len = None, float("inf") for start in range(len(cities)): route = nearest_neighbor(cities, start) length = calc_route_len(route, cities) if length < best_len: best_len = length best_route = route return best_route3.3 贪心构造的价值不是“直接用”,而是“喂给下一步”
我在项目里几乎不会把最近邻的结果直接交给用户。原因是它生成的路径往往存在明显交叉,而交叉线段在欧氏距离下意味着“本来可以更短却绕了远路”。但它的价值在于充当高质量初始解:比随机排列好很多,又比精确算法快几个量级。后面的模拟退火、遗传算法甚至2-opt,都需要一个初始解作为起点,最近邻生成的解就是一个理想的“温启动”。
另外补充一个冷门但有效的构造式算法线索:最小生成树(MST)双重遍历的回路。先对城市图求最小生成树,再沿树边做深度优先遍历得到一个闭合序列,路径长度不超过最优解的2倍。代码上多一步图的构建,效果和最近邻接近,胜在有一个明确的理论界,适合在演示或教学场景用来说明“有界近似”这个概念。实战里我更多还是用多起点最近邻,因为它代码少、常数小。
4. 2-opt局部搜索:用最少的代码换最可观的收益
4.1 交叉线段的平方根:2-opt的核心优化判据
如果你在纸上画一个最近邻生成的路径,最容易看到的缺陷就是路径交叉。两条线段交叉时,一定存在一种“解交叉”的重连方式让总长度下降。2-opt利用的就是这个性质。
设想当前路径的某一段是A→B→C→D,其中B→C之间接了一段长距离,而A→D之间也有一段绕路。如果交换边的连接方式,变成A→C...B→D,并且这两段新边的距离之和小于原来两段边的距离之和,那么整条路径长度就下降了。更准确地说,在对称距离下,2-opt每次尝试反转路径中一段连续区间[i+1..k]的顺序,反转后原有两条边(i,i+1)和(k,k+1)被替换成(i,k)和(i+1,k+1),判断标准就是:
dist(i, k) + dist(k+1, i+1) < dist(i, i+1) + dist(k, k+1)
这个判据成立的前提是距离满足三角几何直觉。对非对称距离,需要改成检查更新后的总长度变化,因为方向反转后每条边的方向都会被翻过来。
4.2 实现一个高效的2-opt迭代器
我最初写过一版对每一对(i,k)都重新算整条路径长度的“笨办法”,数据规模上百就很慢。后来改成只判断局部四条边,优化幅度立刻体现。下面是推荐实现:
def two_opt(route, cities): n = len(route) improved = True while improved: improved = False for i in range(1, n-1): for k in range(i+1, n): a, b, c, d = route[i-1], route[i], route[k], route[(k+1) % n] if dist(cities[a], cities[c]) + dist(cities[b], cities[d]) < dist(cities[a], cities[b]) + dist(cities[c], cities[d]): # 反转区间 route[i:k+1] = reversed(route[i:k+1]) improved = True # 会从头再扫一遍,直到无法改进 return route这个算法最坏复杂度是O(n^3),因为可能存在很多轮改进。实际随机数据上收敛很快,100城一般几十毫秒内结束。如果你想把性能再抠一档,可用“仅在改进发生过的邻域附近继续扫描”的方式做优化,不过工程性价比不高,我在有几十万节点级需求时才会去考虑。这里提醒一个容易踩的坑:区间反转后,路径中“起点城市”的位置可能变化,如果程序里强制规定route[0]为固定城市,反转区间要保证不越过0号下标。上面的实现用i从1开始天然规避了这个问题。
4.3 2-opt的实战地位:与其说是算法,不如说是标配滤镜
我习惯称2-opt是“免费的午餐”。因为它不做复杂搜索,不引入任何新参数,只需要把“能不能缩短”的判断反复做几轮,就能把最近邻的结果拉到接近最优的水平。在我的随机100城实例上,最近邻+2-opt的最终路程只比最优解长1%~3%,这个成绩已经超过了很多复杂元启发算法的基础配置。
顺带提一句3-opt:它每次移除三条边并尝试七种重连方式,理论上比2-opt更强,但实现复杂度和耗时都显著上升。我的建议是,在对称欧氏距离场景,先用2-opt吃到最大收益,再考虑是否值得上3-opt。多数业务场景中,2-opt的结果已经足够交付。
5. 元启发式算法:模拟退火与遗传算法的完整实现和调参经验
5.1 模拟退火:接受“变差”才能逃离局部陷阱
2-opt是爬山式的贪心搜索,只接受变好的解,所以容易陷入局部最优。模拟退火的思想是引入一个概率,允许当前解偶尔“变差”,从而跳出土坑走向更平缓的区域。
模拟退火的三大要素是:温度、降温速率、邻域操作。我的默认配置如下:
def simulated_annealing(cities, T0=5000, alpha=0.995, max_iter=30000): n = len(cities) route = list(range(n)) random.shuffle(route) current = calc_route_len(route, cities) best_route, best_len = route[:], current T = T0 for _ in range(max_iter): # 邻域:随机截取一段区间并反转,本质是随机的2-opt i, j = sorted(random.sample(range(n), 2)) new_route = route[:] new_route[i:j+1] = reversed(new_route[i:j+1]) new_len = calc_route_len(new_route, cities) delta = new_len - current if delta < 0 or random.random() < math.exp(-delta / T): route = new_route current = new_len if current < best_len: best_len = current best_route = route[:] T *= alpha return best_route这里有几个关键点:
- 温度T要足够高,否则一开始就“冻住”,达不到跳坑目的。T0我通常会设置成能接受一次平均变差幅度的数值。简单做法:随机跑几百次邻域操作,统计平均的delta绝对值,把T0设成它的50~100倍。
- 降温系数alpha取值0.99~0.999。越小收敛越快但质量越差,退火行业里叫“退火速度”。0.995是我在100城实例上比较稳的档位。
- 终止条件不一定只看迭代数,也可以加一个“连续N次迭代最优解没有更新就停”的机制,省掉无效计算。
模拟退火的一大优势是容易嵌入任何邻域操作。我在部分项目里把邻域从“区间反转”换成“两点交换”,效果略差,但理解起来更简单,适合给刚接触的同事讲解。两种操作可以随机混合使用,一般能再快一点。
5.2 遗传算法:用种群进化的思路搜索解空间
遗传算法是我接触TSP时最先尝试的元启发式方法。它的调参范围大、自由度很高,但也正因如此特别容易掉进参数陷阱。
它的基本流程是:维护一群路线(个体),每一代让表现好的个体有更高概率进入下一代,并对它们做交叉、变异,产生新个体。关键是编码方式和交叉算子的配合。
TSP的路线天然可以用“城市编号序列”编码。直接套用二进制编码是行不通的,因为一条合法的路线要求序列是1到n的一个排列。交叉算子因此不能用简单的单点交叉,常用的是顺序交叉(Order Crossover, OX)。
def order_crossover(p1, p2): n = len(p1) s, e = sorted(random.sample(range(n), 2)) child = [-1] * n child[s:e+1] = p1[s:e+1] pointer = 0 for g in p2: if g in child: continue while child[pointer] != -1: pointer += 1 child[pointer] = g return child顺序交叉的思想是:从p1中保留一段区间,其余位置按p2中城市的出现顺序填入。这能尽量保持两个亲本的优良片断。变异用简单的两点交换,实现如下:
def swap_mutation(route, mutation_rate=0.1): if random.random() < mutation_rate: i, j = random.sample(range(len(route)), 2) route[i], route[j] = route[j], route[i] return route完整的遗传算法框架如下:
def genetic_algorithm(cities, pop_size=100, generations=300): n = len(cities) pop = [list(range(n)) for _ in range(pop_size)] for ind in pop: random.shuffle(ind) best_route, best_len = pop[0][:], calc_route_len(pop[0], cities) for gen in range(generations): # 适应度:路径越短越好,这里直接用长度,越小越优 fitness = sorted(pop, key=lambda r: calc_route_len(r, cities)) pop = fitness[:pop_size] # 先做一次精英保留 if calc_route_len(pop[0], cities) < best_len: best_len = calc_route_len(pop[0], cities) best_route = pop[0][:] # 锦标赛选择:随机挑若干个体比较,胜者进入交配池 new_pop = pop[:5] # 精英保留5个最优个体 while len(new_pop) < pop_size: cand = random.sample(pop[:50], 3) p1 = min(cand, key=lambda r: calc_route_len(r, cities)) cand = random.sample(pop[:50], 3) p2 = min(cand, key=lambda r: calc_route_len(r, cities)) child = order_crossover(p1, p2) child = swap_mutation(child) new_pop.append(child) pop = new_pop return best_route这个实现刻意保持了最简框架。有几点经验必须交代:
- 选择压力不能太高。如果每一代都只让最优几个个体繁衍,种群多样性会快速崩塌,陷入局部最优。我一般保持一个“前面大约一半的个体有资格参与选择”的大池子。
- 精英保留个数建议5~10。太少容易退化,太多会让新个体的探索被压制。
- 交叉率比变异率重要。交叉负责组合优良片段,变异是防止早熟的后手。变异率设在0.05~0.2之间通常够用。
5.3 两者对比:谁更适合TSP实战
在100城的随机实例上,我的一般观察是:模拟退火收敛稳定,同样的时间预算下结果略优于随机种子的遗传算法,而且几乎不需要纠结种群参数;遗传算法调好了能逼近甚至打平模拟退火,一旦参数没调对就容易出现“跑了很多代但结果很差”的尴尬局面。
两类算法都比纯2-opt略强,但差距并不夸张。如果你对调参没信心,我的推荐是:先跑n遍最近邻+2-opt,拿到一个强baseline,再用模拟退火从这个baseline继续优化。这个组合在大多数场景下已经足够。
6. 算法横向对比与实战选型:评测方法与场景匹配
6.1 用同一批城市做一张公平的对比表
我固定seed=42生成100个城市,然后用文中各算法跑一遍,结果整理成下面这张典型的实验结果表(不同机器和随机数据会波动,但差距量级稳定):
| 算法 | 时间复杂度量级 | 适合规模 | 100城典型结果与最优差距 |
|---|---|---|---|
| 暴力枚举 | O(n!) | n<=10 | 精确最优 |
| 动态规划 | O(n^2·2^n) | n<=22 | 精确最优 |
| 多起点最近邻 | O(n^3) | 几千以内 | 约12%~15% |
| 最近邻+2-opt | O(n^3)但常数极小 | 几百以内 | 约1%~3% |
| 模拟退火 | O(迭代数·n) | 几百以内 | 约0%~2% |
| 遗传算法 | O(种群数·代数·n) | 几百以内 | 约0%~3% |
这张表揭示了一个很多人忽视的事实:在很多实际业务里,2-opt这种简单方法得到的解已经“足够好”。所谓“足够好”,取决于业务对精度的要求。比如给一套配送方案,总里程提升2%和提升0.5%,在时间成本上差别可能只有几分钟到十几分钟,但如果你为此把程序改成复杂的元启发式,开发和调参成本可能高得多。
6.2 选型建议:不同规模不同需求的不同答案
如果你的n在20以内并且要绝对最优,直接上动态规划,没有悬念。n在100以内、能接受近似解,多起点最近邻+2-opt是性价比最高的选择,代码短、无参数、效果好。n在200到500之间,且希望结果更稳,就在前面基础上加模拟退火,把最近邻的结果当初始解,能额外压掉一点路程。n到达几千甚至上万,2-opt和模拟退火的O(n^3)或O(迭代数·n)也会吃力,这时需要转向更工程化的方案,比如先对城市做空间聚类、分区域求解再拼接,或者用专门针对大规模TSP的Lin-Kernighan类方法,这已经超出本文的讨论范围。
我给一个保守但可靠的推荐组合:多起点最近邻用于生成初始解,2-opt做到收敛,模拟退火用“随机反转片段”的邻域做后续精修。这个组合我在好几个实际路径优化项目里都验证过,效果稳定,且不需要不停试探调参。
6.3 城市坐标数据的预处理与后处理细节
这里补充一个看起来不起眼、实际上很影响最终结果的操作:城市编号最好从0开始并保持与数据表一一对应;距离提前算好放进二维矩阵,避免在循环里反复调用欧氏距离函数。对100城,距离矩阵只有100×100,占用极小,但能大幅加速所有算法。别小看这个优化,模拟退火每次邻域评估都要取若干对城市距离,预计算能让总耗时直接砍半。
我的预计算代码通常是:
def build_distance_matrix(cities): n = len(cities) return [[dist(cities[i], cities[j]) for j in range(n)] for i in range(n)]后续所有算法都把dist函数换成查表。这样改不改变任何算法逻辑,但工程上非常值得。
另外关于随机种子:如果你在调参,固定随机种子能让你判断“是算法参数变好了还是运气变好了”;如果你在给客户出方案,建议多跑几个随机种子,取最优结果,因为启发式算法带随机性,重复几次取最优往往比单纯加迭代次数更划算。这是个经常被忽略的实用技巧。
多说一句我自己的习惯。每次跑完算法,我都会顺手把最优路径画到散点图上,看一遍线段交叉和整体形状。算法给出的结果如果交叉明显,多半是某个环节出了偏差;如果整体形态舒展,即使路程数值不是最漂亮,也更有可能是可信解。图像直觉虽然不严谨,但在TSP这种空间特征很直观的问题里,它常常能帮你提前发现代码里的隐蔽bug。这也是我建议你在复现本文代码时一定要做的一步。