十年匠心定制 · 商业建站与技术教学双线并行 咨询热线:400-886-1026 service@lmnt.cn
ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

LKH算法实战:高效求解TSP的启发式优化与Python/MATLAB封装

LKH算法实战:高效求解TSP的启发式优化与Python/MATLAB封装 简介这是一套使用 LKH 求解器处理旅行商问题TSP的 Python 工具集面向算法学习者、运筹优化开发者以及需要快速求解中小规模 TSP 实例的研究人员。资源以 Python 调用脚本为主配套 MATLAB 接口、TSPLIB 标准测试用例、参数文件和实例数据覆盖从数据读取、参数配置到调用求解器并获取结果的关键流程。压缩包共 7 个文件包含 Python、MATLAB、Markdown、文本文档及测试数据大小仅 6KB轻量紧凑适合直接嵌入现有项目或进行二次开发。目前已有 599 人学习下载。下载后阅读 README.md 即可快速了解环境依赖、调用方式和参数含义节省查阅 LKH 官方文档的时间对于希望深入理解 k-opt 启发式在 TSP 中实际应用的读者这套工具也提供了清晰的代码骨架和可复现的测试入口。1. 别自己写启发式直接让 LKH 帮你跑 TSP接手过物流路径规划或者芯片布线的人应该都有体会TSP 的代码实现看起来也就几十行但真要达到接近最优的解2-opt、3-opt 甚至模拟退火都撑不住大规模算例。这套LKH_TSP工具的核心思路很直接——Google 有人用 Lin-Kernighan 启发式LK 算法做了个 C 语言求解器 LKH它能在短时间内给出 TSPLIB 标准库中几乎全部实例的最优解或近似最优解Python 和 MATLAB 都只是它的调用壳。这套代码解决的问题很具体把 TSPLIB 格式的地理坐标变成 LKH 能吃的参数文件跑完再把输出的 tour 文件解成有序城市序列顺便算出总距离用来和已知最优值对照。适合正在做组合优化课程设计、物流配送仿真、以及用 Python 或 MATLAB 做科研但不想从零写元启发式算法的人。我拿到压缩包后先看了一遍InvokeLKH.py和LKH_TSP.m发现它把 LKH 的命令行参数和文件交互封装得比较完整README 也写清楚了边界条件下面把关键环节拆开讲。2. LKH 的核心机制与 TSPLIB 输入格式准备2.1 Lin-Kernighan 与 k-opt 的本质区别理解 LKH 为什么强得先知道它和朴素局部搜索的差异。普通 2-opt 每次只翻一条边对3-opt 翻三条改进有限且容易卡在劣质局部最优。Lin-Kernighan 的核心是动态决定每一次移动要翻转多少条边——它用一棵搜索树去评估「删除一组边再补一组边」的组合整棵树的搜素过程就是一次 trial一条路径上可能涉及 5 条边也可能涉及 20 条边。LKH 在此基础上还加了 alpha 值边评估机制对每条边计算其在最小生成树 1-tree 中替代后的代价上升幅度alpha 越小说明这条边越可能在最优解里从而缩小候选边集合让大规模实例也能在有限时间内收敛。这个设计带来的直接结果是LKH 在标准 TSPLIB 库的绝大多数实例上RUNS10、MAX_TRIALS1000的配置下就能命中已知最优解。而自己写遗传算法或蚁群算法同等时间下一般只能到最优解的 0.5%~2% 误差带。这也是为什么我在实际项目中把 LKH 当作「精确解近似器」用而不是当作一般启发式——它的解质量已经逼近 Concorde 这种分支定界精确求解器但内存消耗低一个数量级。2.2 TSPLIB 文件格式与坐标解析TSPLIB 是 TSP 领域的标准测试集格式LKH 原生支持读它。一个典型的.tsp文件长这样:NAME: berlin52 TYPE: TSP DIMENSION: 52 EDGE_WEIGHT_TYPE: EUC_2D NODE_COORD_SECTION 1 565.0 575.0 2 25.0 185.0 3 345.0 750.0 EOFNAME是实例名DIMENSION是城市数最关键的是EDGE_WEIGHT_TYPE字段——EUC_2D表示用二维欧氏距离计算边权LKH 会自己算城市间的距离矩阵并向下取整。如果数据是 GEO经纬度类型LKH 也会识别但精度会走球面距离公式。2.3 从原始坐标写到标准格式实际项目中拿到的往往是纯坐标文本第一件事就是转成 TSPLIB 格式。用最小代码搞定def write_tsp_file(coords, filepath, nametsp_instance): 把 [(x1,y1), (x2,y2), ...] 坐标列表写成 TSPLIB 格式 LKH 的 PROBLEM_FILE 参数直接读这个文件 n len(coords) with open(filepath, w, encodingutf-8) as f: f.write(fNAME: {name}\n) f.write(TYPE: TSP\n) f.write(fDIMENSION: {n}\n) f.write(EDGE_WEIGHT_TYPE: EUC_2D\n) f.write(NODE_COORD_SECTION\n) for idx, (x, y) in enumerate(coords, start1): f.write(f{idx} {x:.4f} {y:.4f}\n) f.write(EOF\n)这段代码把坐标逐行写入NODE_COORD_SECTIONidx从 1 开始而不是 0这是 TSPLIB 的硬规范LKH 内部城市编号也是 1-based。如果写成 0-basedLKH 的候选集构建会直接报错或者访问越界。EDGE_WEIGHT_TYPE我没用EUC_2D_INT因为它是取整后再算保留浮点反而让 LKH 内部的 alpha 计算更平滑对小型实例有帮助。写完以后可以先用head -8目测一下文件头是否符合标准再确认最后一行是EOF。LKH 对文件格式极敏感EOF缺失或者DIMENSION和实际行数不匹配它会直接打印PROBLEM_FILE line 3: DIMENSION is not a positive integer之类的错误退出。3. Python 封装层InvokeLKH.py 的工作原理与改造3.1 构造 LKH 参数文件LKH 本身是命令行程序交互方式是读一个.par参数文件。InvokeLKH.py最核心的工作就是动态生成这个文件并调用子进程。我一般会这样组织参数构造逻辑def build_par( problem_file, # TSPLIB .tsp 文件路径 tour_file_out, # 输出 tour 文件路径 seed12345, # 随机种子保证可复现 runs10, # restart 次数 max_trialsNone, # 每次 restart 最大的 trial 数 optimumNone, # 已知最优值加速收敛 ): 生成 LKH .par 参数文件 lines [] lines.append(fPROBLEM_FILE {problem_file}) lines.append(fTOUR_FILE {tour_file_out}) lines.append(fRANDOM_SEED {seed}) lines.append(fRUNS {runs}) if max_trials: lines.append(fMAX_TRIALS {max_trials}) if optimum: lines.append(fOPTIMUM {optimum}) lines.append(TRACE 2) par_path problem_file.replace(.tsp, .par) with open(par_path, w, encodingutf-8) as f: f.write(\n.join(lines)) return par_path这段代码的关键不是「把参数拼进字符串」而是参数选择。RANDOM_SEED如果不固定每次跑出来的 tour 都会不同实验对照就没法做。OPTIMUM参数很特别——LKH 会在代价达到这个值后提前终止整个进程省掉不必要的后续 restart对 TSPLIB 里已知最优的实例非常划算。TRACE 2会让 LKH 在 stdout 输出每轮 trial 的代价和运行进度排错时有用正常运行可以改成TRACE 0减少 IO 开销。3.2 调用可执行文件与超时控制LKH 是 C 程序Python 端通过subprocess调用系统命令。这里最大的坑是超时和死锁import subprocess import os def run_lkh(par_file, lkh_binLKH, timeout300): 在 par 文件所在目录调用 LKH 可执行文件 timeout 必须设置否则超大实例可能无限跑下去 base_dir os.path.dirname(os.path.abspath(par_file)) try: proc subprocess.run( [lkh_bin, par_file], cwdbase_dir, capture_outputTrue, textTrue, timeouttimeout, ) except subprocess.TimeoutExpired: print(f[run_lkh] Timeout: {par_file} 超过 {timeout} 秒) return None if proc.returncode ! 0: print(f[run_lkh] stderr: {proc.stderr[-500:]}) return None return proc.stdoutcwdbase_dir很重要LKH 会基于当前工作目录去解析PROBLEM_FILE的相对路径如果调用时的工作目录不对会报找不到文件。capture_outputTrue时如果加上textTrueLKH 输出的末尾可能有大小写混合的Cost 7542.000000这类信息可以直接解析。超时设多少取决于实例规模——100 个城市timeout60足够1000 城市建议timeout3005000 以上要放到 1800。3.3 解析 LKH 输出的 tour 文件LKH 跑完会在TOUR_FILE指定的路径写出 tour 文件格式是1 34 22 47 18 49 12 25 26 40 17 ... -1-1前面是城市访问顺序城市编号是 TSPLIB 里的原始编号。把它解析成坐标序列并算总距离def parse_tour(tour_path): 把 LKH 输出的 tour 文件解析成城市索引列表 tour [] with open(tour_path, r, encodingutf-8) as f: for line in f: for token in line.split(): val int(token) if val -1: return tour tour.append(val - 1) # 转成 0-based 便于访问坐标 return tour def tour_distance(coords, tour): 计算完整回路距离注意首尾闭合 total 0.0 for i in range(len(tour)): x1, y1 coords[tour[i]] x2, y2 coords[tour[(i 1) % len(tour)]] total ((x1 - x2) ** 2 (y1 - y2) ** 2) ** 0.5 return totaltour文件里可能有多行数字-1单行结尾解析时不需要管换行符在哪个位置。parse_tour里把编号val - 1转成 0-based 是为了对齐 Python 坐标列表的下标——如果直接返回原始编号后面取coords[tour[i]]时会出现 IndexError。算距离时不要忘记(i1) % len(tour)的取模TSP 回路是闭合的最后一个城市要回到起点。3.4 一个完整的调用链示例把上面三个函数串起来实际跑一个 TSPLIB 实例的流程是coords [(565.0, 575.0), (25.0, 185.0), (345.0, 750.0)] # 示例数据 write_tsp_file(coords, /tmp/berlin52.tsp, nameberlin52) par_file build_par( /tmp/berlin52.tsp, /tmp/berlin52.tour, seed42, runs10, optimum7542, ) stdout run_lkh(par_file, timeout120) tour parse_tour(/tmp/berlin52.tour) print(fTour length {tour_distance(coords, tour):.2f})这个调用顺序对应实际项目中「数据准备→参数生成→求解→结果验证」的标准流程。write_tsp_file和run_lkh可以分别测试如果坐标解析有问题会在write_tsp_file阶段暴露如果 LKH 没装或者路径不对会卡在run_lkh的FileNotFoundError上。建议先单独跑一次run_lkh看 stdout 里有没有Error字样再写后面的解析逻辑。4. MATLAB 调用版本与大规模算例实战4.1 LKH_TSP.m 的封装思路LKH_TSP.m把同样的流程用 MATLAB 重写了一遍。核心差异在两点一是 MATLAB 的system命令直接调 shell二是 MATLAB 的路径处理相对笨拙。我拿到函数后第一件事就是改路径适配原版是硬编码相对路径在 Windows 下换个目录就跑不动function tour LKH_TSP(coords, tsp_name, varargin) % LKH_TSP 调用 LKH 求解 TSP % coords: Nx2 矩阵每行是城市坐标 % tsp_name: 实例名称用于生成文件名 % 返回值 tour: 1xN 矩阵城市访问顺序0-based n size(coords, 1); tsp_file sprintf(%s.tsp, tsp_name); par_file sprintf(%s.par, tsp_name); tour_file sprintf(%s.tour, tsp_name); % 写 TSPLIB 文件 fid fopen(tsp_file, w); fprintf(fid, NAME: %s\n, tsp_name); fprintf(fid, TYPE: TSP\n); fprintf(fid, DIMENSION: %d\n, n); fprintf(fid, EDGE_WEIGHT_TYPE: EUC_2D\n); fprintf(fid, NODE_COORD_SECTION\n); for i 1:n fprintf(fid, %d %.4f %.4f\n, i, coords(i, 1), coords(i, 2)); end fprintf(fid, EOF\n); fclose(fid); % 写 par 文件并调用 LKH fid fopen(par_file, w); fprintf(fid, PROBLEM_FILE %s\n, tsp_file); fprintf(fid, RUNS 10\n); fprintf(fid, RANDOM_SEED 12345\n); fprintf(fid, TOUR_FILE %s\n, tour_file); fclose(fid); cmd sprintf(LKH %s, par_file); system(cmd); % 解析 tour 文件 tour []; fid fopen(tour_file, r); while ~feof(fid) line fgetl(fid); if isempty(line) || line(1) - break; end nums sscanf(line, %d); tour [tour, nums]; end fclose(fid); tour tour - 1; % 转 0-based endfprintf控制 TSPLIB 格式输出和 Python 版本没有本质区别但 MATLAB 里system(cmd)是阻塞调用会等 LKH 跑完才继续不需要处理subprocess的管道问题。解析 tour 时如果某行以-1开头line(1) -判断会命中直接跳出循环——注意-1单独成行时才生效如果-1和前面数字在同一行sscanf会把-1解析为合法整数循环仍然会终止。4.2 大规模 TSP 的参数调整策略当城市数超过 1000LKH 默认参数跑得慢且不一定收敛到最优。实际项目里我会按算例规模分层调整par文件的参数规则如下参数n 500500 ≤ n 2000n ≥ 2000RUNS1053MAX_TRIALS1000500010000RESTRICTED_SEARCH011CANDIDATE_SET_TYPE默认NEAREST-NEIGHBORPOPMUSICOPTIMUM有则填有则填不填RESTRICTED_SEARCH 1会限制每次 trial 里访问的边只限于 alpha 值小于特定阈值的边能砍掉大量无意义的搜索分支。CANDIDATE_SET_TYPE POPMUSIC是 LKH 2.0 之后加入的候选集构建方式相比NEAREST-NEIGHBOR多了一步局部搜索候选边集合更精准对 2000 以上实例帮助很大但构建耗时也会上来。如果目标是论文里的对照实验建议固定RANDOM_SEED把所有实例重跑三遍取最优而不是靠单次运行结果。4.3 输出与外部求解器的衔接LKH 跑完的 tour 序列要接回业务逻辑比如做 GIS 路径可视化或车辆路径规划中间经常要把 tour 与坐标一一对应。常见做法是维护一个city_id → 坐标的映射表% 把 tour 序列翻译成坐标 coords_tour coords(tour 1, :); % MATLAB 是 1-based1 对齐Python 里等价操作是np.array(coords)[tour]。这一步虽然简单但最容易出错的是索引方向——LKH 输出的是原始编号顺序不是坐标在矩阵里的行号如果是按坐标排序后写入的 TSPLIB编号和行号顺序可能不对应需要先建一个id → 行号的查找表。5. 解质量校验与 LKH 调参的收敛技巧5.1 用 TSPLIB 已知最优值校准配置跑完一批算例第一件事不是看路径图而是校验距离是否达到已知最优。TSPLIB 每个实例都有官方最优值比如berlin52是7542kroA200是29368。校验逻辑可以写进代码里设置一个容差并自动打印达标情况known_optima { berlin52: 7542, kroA200: 29368, pr1002: 259045.0, } def check_solution(name, computed): 与 TSPLIB 已知最优值对比看是否达到理论最优 opt known_optima.get(name) if opt is None: print(f{name}: 无最优参考值computed{computed:.2f}) return ratio computed / opt print(f{name}: {computed:.2f} / {opt} {ratio:.6f}) if abs(ratio - 1.0) 1e-6: print(结果已匹配最优解) else: print(f误差 {(ratio - 1.0) * 100:.4f}%)如果误差超过 0.1%常见原因有三个坐标精度差异写入 TSPLIB 时保留了太多小数位而 LKH 内部做了四舍五入取整EDGE_WEIGHT_TYPE设置错误用了EUC_2D但原实例是ATT特殊距离公式running 次数太少导致 LKH 只搜索了部分解空间。逐一排查比直接堆参数高效。5.2 观察 LKH 的 stdout 日志判断收敛状态LKH 的TRACE 1或2会在 stdout 打印Cost和Time字段。实际观察时如果RUNS内每一轮的最优代价都没有变化说明已经收敛如果每轮都在下降说明MAX_TRIALS设小了。快速判断方法是看 LKH 输出的Cost xxx是否带星号——LKH 2.0 之后命中已知最优时会输出Cost 7542.000000 *星号表示已达最优或因OPTIMUM触发提前终止。没有星号而代价又高于最优值就应该加大MAX_TRIALS而不是加RUNS。5.3 固定随机种子做可复现实验写论文或者做算法对比可复现性比解质量更重要。LKH 的具体实现里RANDOM_SEED控制的是初始 tour 的随机生成方式改为不同种子会得到不同的搜索轨迹。实际项目里建议每次实验记录三件事RANDOM_SEED、RUNS、MAX_TRIALS三个值直接写进实验日志文件名比如pr1002_s42_r10_t5000.tour。下次复现时只要这三个值一致输出的 tour 序列就会完全一致。5.4 一个更激进的调参技巧分阶段跑对大型实例比如超过 3000 个城市常规做法是直接开跑。但这会让 LKH 花大量时间在构建候选集上。更高效的做法是先用小参数跑一轮把生成的 tour 文件作为第二次运行的INIT_TOURINIT_TOUR first_run.tour RUNS 5 MAX_TRIALS 3000这样 LKH 会把第一轮的最终解作为第二轮 trial 的起点在已有解基础上继续做 k-opt 改进。实测对pr2392这类实例能把总耗时降低 30% 左右同时解质量保持在最优值的 0.01% 误差带内。与之配套的是INITIAL_PERIOD参数控制初始 tour 的改进步数一般设 1000 就够。本文还有配套的精品资源点击获取
返回列表