从TSP到VRP:电工杯B题车辆无人机协同与多中心选址求解

从TSP到VRP:电工杯B题车辆无人机协同与多中心选址求解 简介本资源为2022年电工杯数学建模竞赛B题「5G网络环境下应急物资配送问题」的完整论文资料面向备战数学建模竞赛的高校学生及对路径规划算法感兴趣的读者可用于赛题复盘、算法选型参考与论文写作学习。压缩包内为1个PDF文件约1.01MB即该赛题的一等奖级别全文论文含摘要、问题重述、模型假设、符号说明及各问模型建立与求解等完整章节可直接对照题目研读建模思路。论文将类旅行商问题转化为车辆路径规划模型问题一采用模拟退火与深度优先搜索求解得配送里程582km问题二用粒子群优化结合广度优先搜索得配送时间380分钟问题三、四以K-means划分区域后构建遗传算法车辆路径模型并讨论两个物资集中地点的选址。读者可从中获取启发式算法的建模框架、目标函数与约束写法、流程图示及结果验证方法如穷举法与SADFS的耗时对比适合作为赛前模板与算法积累。目前已有4346人学习下载。1. 从类 TSP 到 VRP2022 电工杯 B 题的建模边界在哪2022 年电工杯 B 题给的场景很具体14 个地点、物资集中在第 9 个点、车辆载重 1000kg、总需求 782kg第一问只让车辆跑。很多人拿到邻接矩阵就条件反射地套 TSP结果卡在 2、3、4 这几个节点上——图是非完全图从 3 号点出来只能回 5 号点绕不开重复经过。782kg 没超过 1000kg车不必回配送中心可路径里又必须出现重复节点这已经不是标准 TSP 了。到了第三问载重砍到 500kg必须至少回中心一次第四问 30 个点、1568kg、要选两个集散点问题彻底滑向 VRP 加设施选址。真正决定模型成败的不是算法名字而是先判断清楚图的连通结构、载重与需求量的比例关系再决定要不要拆路径、拆几段。这套判断链条对做物流调度、无人机协同配送的人同样是通用技能。2. 邻接矩阵预处理与 SADFS 求解类 TSP 路径第一问的核心矛盾是图不连通成完全图但需求又允许一趟跑完于是路径必然带重复节点。直接上模拟退火而不做预处理搜索空间里全是非法解迭代效率会被拖死。所以顺序必须是先修矩阵、再搜路径、最后退火微调。2.1 邻接矩阵的两种填充值不能混从附件读进来的距离表空缺项和对角线要分开处理否则算法会把不可达当成一条超长边反复走位置填充值理由自身到自身0不产生额外代价避免自环被计费不可达节点对9999足够大但不溢出贪心搜索能主动避开真实可达边原值单位 km最大载重与里程约束都用它import numpy as np INF 9999 # 14 个地点的邻接矩阵行列为 0 基索引 dist np.array([ [0, INF, INF, INF, 54, INF, 55, INF, INF, INF, 26, INF, INF, INF], [INF, 0, 56, INF, 18, INF, INF, INF, INF, INF, INF, INF, INF, INF], # ...省略中间行按附件 1 逐行填入 ], dtypefloat) # 对角线强制归零防止数据源里残留非零值 np.fill_diagonal(dist, 0) # 不可达统一为 INF别用 np.inf否则累加后出现 nan dist[dist 0] 0这段的关键是INF取 9999 而不是np.inf。路径代价做累加和比较时np.inf一旦参与减法或概率计算会传染成 nan后面 Metropolis 准则里的指数运算直接崩掉。2.2 DFS 先给出一条合法初始路径模拟退火需要一个可行的起点否则前期大量迭代浪费在修复非法路径上。DFS 负责从配送中心出发尽量遍历所有节点def dfs_path(start, dist, n): best [start] visited {start} def dfs(cur, path): nonlocal best if len(path) len(best): best path[:] # 按距离升序试邻居优先走短边 for nxt in sorted(range(n), keylambda x: dist[cur][x]): if nxt in visited or dist[cur][nxt] 9999: continue visited.add(nxt) path.append(nxt) dfs(nxt, path) path.pop() visited.remove(nxt) dfs(start, [start]) return bestvisited用集合而不是布尔数组是因为这个图存在必须重复经过的节点DFS 阶段的访问过只约束本次尝试不能全局锁死。sorted(..., keydist[cur][x])让搜索优先走短边得到的初始路径质量更高退火收敛更快。2.3 扰动加 Metropolis 准则才是退火的本质初始路径定下来后对路径施加扰动生成新路径按代价差决定是否接受。核心是允许以一定概率接受更差的解跳出局部最优import math, random def sa_solve(dist, init, T01000.0, alpha0.95, iters2000): cur init[:] cur_cost path_cost(cur, dist) best, best_cost cur[:], cur_cost T T0 for _ in range(iters): new perturb(cur, dist) # 交换两个节点或插入新节点 new_cost path_cost(new, dist) delta new_cost - cur_cost # delta 0 直接接受否则按 exp(-delta/T) 概率接受 if delta 0 or random.random() math.exp(-delta / T): cur, cur_cost new, new_cost if cur_cost best_cost: best, best_cost cur[:], cur_cost T * alpha return best, best_costT01000给的是初始温度量级参考单段最长边 61km 的十几倍保证早期接受率足够高alpha0.95是降温系数2000 次迭代后温度约剩 1000×0.95²⁰⁰⁰早已趋近 0末段基本只接受更优解。算例里迭代约 200 次结果就稳定不再变化说明iters2000留了足够冗余。2.4 用穷举反向验证而不是自我感觉这类题最怕算法跑出来就说最优。第一问规模小可以穷举对照总共 147258 种可行路径穷举耗时 139.816sSADFS 只要 3.79s拿到的最优路径是 9→13→14→10→6→4→6→5→3→2→5→7→1→11→12→8→9总里程 582km。两套结果一致说明加速没有牺牲精度。规模一上来穷举必然失效这个对照只做一次用来标定模型可信度就够。3. 车辆无人机协同的 PSOBFS 求解与 epoch 调参第二问引入无人机后车辆能走实线无人机还能走虚线。两条路径同时可变直接联合优化维度太高论文选择先定车、再定机这个解耦思路在车机协同调度里是常规操作。3.1 先车后机解耦降低搜索维度物资总量仍是 782kg车辆一趟能装完所以车不用回第 9 点。需求量大、或者只走实线的点交给车辆其余点用无人机覆盖。判断顺序是先把每个地点的物资需求量排序需求量大的地方车辆必要性高排出优先级序列再据此确定车辆主干路径。这样做的好处是把车辆路径 × 无人机路径的组合爆炸压成了两级搜索——外层改车内层重算机。3.2 BFS 枚举无人机的可达投放点车辆主干确定后每个车辆停靠点都可能放无人机。可用节点集合要动态维护from collections import deque def bfs_drone(root, avail, drone_adj): 从车辆停靠点 root 出发BFS 找出无人机可达且尚未被占用的节点 q deque([root]) seen {root} reachable [] while q: cur q.popleft() for nxt in drone_adj[cur]: if nxt in seen or nxt not in avail: continue seen.add(nxt) reachable.append(nxt) q.append(nxt) return reachableavail是全部节点减去车辆路径已用和即将使用的节点每趟循环结束后要把本次无人机用过的点从avail里减掉否则同一个点会被车和无人机重复服务。drone_adj包含虚线边和车辆邻接矩阵不是同一张表必须单独构建这点在代码里最容易出错。3.3 PSO 的粒子就是一组无人机指派方案把 BFS 搜出的候选指派方案作为初始粒子群适应度同时看时间和里程符号含义调参影响适应度时间与里程的加权归一化指标权重偏时间则倾向多派无人机epoch迭代轮数论文取 10000粒子位置各车辆段是否派无人机及其目标点决定解空间覆盖度def fitness(plan, t_norm, s_norm, w10.6, w20.4): # t_norm、s_norm 分别把时间和里程归一化到同一量纲 return w1 * t_norm(plan) w2 * s_norm(plan) def pso_run(particles, epochs10000): gbest min(particles, keyfitness) for _ in range(epochs): for p in particles: p.update(gbest) # 向全局最优靠拢 if fitness(p) fitness(gbest): gbest p return gbestw1取 0.6 是把时间作为主目标因为题面问的是配送时间epochs10000不是拍脑袋论文专门跑了不同轮数对比10000 轮时拿到最优解的概率超过 98%。3.4 结果复盘40 次只中 23 次的真实原因最终车辆路径为 9→8→7→5→2→5→6→10→9各段对应的无人机任务如下车辆段无人机路径9→89→13→88→78→12→77→57→11→1→26→106→3→4→1010→910→14→9总时间 380 分钟。看起来漂亮但实际跑 40 次只有 23 次收敛到该结果。原因是 PSO 的修改是自上而下、只动后方不动前方初始车辆路径一旦选歪后续任何无人机指派都救不回来直接锁死在局部最优。这个坑的通用解法是加多次随机重启或者对车辆路径本身也做小扰动。另外该问题可行解超过 2000000 种穷举验证不现实只能横向比较多次运行结果用 50% 以上的命中率说明模型可用。4. 载重 500kg 下的 K-meansGA 车辆路径规划第三问把载重压到 500kg需求 782kg一趟装不下车辆至少要回配送中心一次。这时候再当成整条 TSP 处理就错了得换成 VRP每一次出中心—服务若干点—回中心算一条独立路径多条路径组成总方案。4.1 染色体编码全排列加分隔符GA 里最讲究的是编码方式这里用非完全排列 若干 0来表示方案import random def encode(nodes, num_routes2): 全排列后插入 num_routes-1 个 00 把序列切成多段每段是一条车路径 perm nodes[:] random.shuffle(perm) for _ in range(num_routes - 1): perm.insert(random.randint(1, len(perm) - 1), 0) return perm def decode(chrom): routes, cur [], [] for g in chrom: if g 0: if cur: routes.append(cur) cur [] else: cur.append(g) if cur: routes.append(cur) return routes比如6 8 7 9 0 6 5 4 1 2 3解码成两条路径第一条服务 8、7、9第二条服务 6、5、4、1、2、3起点和返回中心在解码时自动补齐不用写进染色体能省一大截长度。4.2 合法性校验与适应度函数随机生成的个体经常超载必须在分割阶段按载重和里程重新划分def is_valid(routes, demand, cap500): return all(sum(demand[n] for n in r) cap for r in routes) def fitness(chrom, dist, demand, w(0.5, 0.3, 0.2)): routes decode(chrom) if not is_valid(routes, demand): return -1e9 # 非法个体直接给极低适应度 t total_time(routes, dist) s total_distance(routes, dist) n len(routes) wt, ws, wn w return -(wt * t ws * s wn * n) # 取负越大越好适应度里同时压时间、里程和车辆数三个权重都取正数方向是越小越好所以整体取负号转成最大化问题。is_valid里的 500 就是本题载重上限换成第四问的数值即可复用。4.3 只做变异、不做交叉的取舍标准 GA 有交叉和变异两步这里刻意去掉了交叉理由是交叉容易产生非法个体还得额外写修复逻辑而变异本身就带交换语义反复交换已经足够维持种群多样性def mutate(chrom, dist): new chrom[:] i, j random.sample(range(len(new)), 2) new[i], new[j] new[j], new[i] if not legal_by_graph(decode(new), dist): return chrom # 变异失败则原样返回 return newlegal_by_graph检查相邻两节点在邻接矩阵里是否可达。写法上先交换再校验非法就回滚好处是变异后个体必然合法不需要额外的修复算子。代价是部分变异被浪费但换来的是实现简单和收敛稳定。4.4 聚类预分区与收敛判定问题三先用 K-means 把节点硬聚类成两片8、7、5、2、1、11、12、13 一片3、4、6、10、14 一片两片分别跑 GA。最终两条路径为路径 1 车走 9→6→10→9路径 2 车走 9→8→7→5→2→5→9第一趟 238kg、第二趟 494kg总时间 408 分钟。from sklearn.cluster import KMeans coords load_coords() # 各地点坐标 km KMeans(n_clusters2, n_init100, random_state0) labels km.fit_predict(coords)n_init100是为了降低随机初始化带来的波动问题四里索性把整个聚类加 GA 的过程重复跑 100 次只有 100 次结果的极差落在阈值内才认定收敛。这种内层 GA 外层聚类重复的双层验证代价是时间换来的是结果可信。5. 多中心选址与 OR-Tools 对照校验的实操技巧第四问 30 个点、总需求 1568kg、要选两个集散中心思路是先把 30 个点用 K-means 分成两片再对每片当成独立的 VRP 跑 GA。聚类结果落在第 5 点和第 20 点第一类覆盖 1、2、3、4、5、6、7、8、9、10、11、12、13、14、18第二类覆盖 15、16、17、19、20、21、22、23、24、25、26、27、28、29、30。第二类明显更分散这也是最终两条路径长度差异大的直接原因。实际方案里第一类车辆分两次出发第二类同样分两次无人机按段补盲具体配对如下表类别车辆路径无人机路径第一类第一次5→2→5→7→52→1→2第一类第二次5→9→6→4→6→10→9→59→8→9、9→13→9、9→12→9、4→10→14→10第二类第一次20→16→15→19→24→25→2016→21→16、25→29→25第二类第二次20→26→28→26→30→26→27→26→2020→25→20、27→22→27想让上面这套手写启发式的结果更硬可以用 OR-Tools 的约束求解器做一次独立对照附录里的脚本就是这么用的from ortools.constraint_solver import routing_enums_pb2 from ortools.constraint_solver import pywrapcp def create_data_model(): data {} data[distance_matrix] DIST # 从附件读入不可达填 999999 data[num_vehicles] 1 data[depot] 3 # 注意这里的索引口径 return data def main(): data create_data_model() manager pywrapcp.RoutingIndexManager( len(data[distance_matrix]), data[num_vehicles], data[depot]) routing pywrapcp.RoutingModel(manager) def distance_callback(from_index, to_index): f manager.IndexToNode(from_index) t manager.IndexToNode(to_index) return data[distance_matrix][f][t] cb routing.RegisterTransitCallback(distance_callback) routing.SetArcCostEvaluatorOfAllVehicles(cb) params pywrapcp.DefaultRoutingSearchParameters() params.first_solution_strategy ( routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC) solution routing.SolveWithParameters(params) print(solution.ObjectiveValue() if solution else no solution) if __name__ __main__: main()几个实操要点。第一depot的索引口径必须和主程序统一如果 SA/DFS 用的是从 1 开始的地点编号接入 OR-Tools 前要整体减 1附录脚本里depot 3与手写程序的编号方式并不一致直接照抄会导致起点错位、结果对不上。第二PATH_CHEAPEST_ARC只负责生成第一个可行解属于构造式启发式对多车辆、带载重约束的场景必须再加Dimension约束单靠它不会考虑容量。第三distance_matrix里不可达边填 999999 这类大整数求解器会主动绕开但填得太大叠加后会溢出取值控制在一百万量级比较稳。第四验证时不要只对比总代价把路径序列打出来逐点核对OR-Tools 对等价路径的节点顺序和手写算法经常不同只比数字会误判成结果不一致。重复跑聚类加 GA 100 次、观察极差是否收敛比盯单次运行结果更能说明模型到底稳不稳。本文还有配套的精品资源点击获取