1. 项目缘起当紧急药品配送遇上城市“毛细血管”去年我参与了一个社区医疗应急保障系统的前期调研。当时一个很具体的问题摆在我们面前如何在老旧城区或突发交通管制的情况下将急救药品从社区卫生中心快速送到多个分散的居民点传统的车辆配送受限于道路拥堵和“最后一公里”的通行难题而单纯的人工徒步又效率太低。这时无人机配送就成了一个非常值得探讨的技术选项。但问题没那么简单。我们手头有几架续航约25分钟、载重2公斤的多旋翼无人机需要服务的点可能有8到15个分布在大约3平方公里范围内。目标很明确在电池耗尽前让无人机访问所有配送点并且总飞行距离尽可能短。这听起来就像经典的“旅行商问题”TSP——找到访问所有城市一次并回到起点的最短路径。然而TSP是出了名的计算难题当配送点超过10个时精确求解的计算量会爆炸式增长对于需要快速响应的应急场景来说这显然不现实。于是我们转向了启发式算法也就是不求最优解但求在可接受时间内找到一个“足够好”的方案。在尝试了贪婪算法、遗传算法后模拟退火算法因其原理清晰、实现相对简单、且容易跳出局部最优解的特性成为了我们的重点试验对象。这个项目就是基于那次实践用Matlab将“距离优先”的无人机药品配送路线规划过程完整复现和解析一遍。它不是天马行空的设想而是针对一个具体约束条件距离近优先下的工程化求解思路。2. 核心问题建模把现实世界抽象成算法能理解的语言在写任何代码之前我们必须把实际的配送问题翻译成数学模型和算法可以处理的数据。这一步是后续所有工作的基石如果模型建偏了代码再漂亮也是南辕北辙。2.1 场景定义与输入数据首先我们需要定义场景的基本要素配送中心Depot药品的出发点和最终返回点通常只有一个。我们将其坐标设为 (0, 0)。配送点Customers需要接收药品的位置。我们随机生成N个点的二维坐标 (x, y)模拟它们在区域内的分布。在实际项目中这些坐标来自GPS或地理信息系统。距离矩阵Distance Matrix这是整个问题的核心数据。我们计算所有点两两之间的欧几里得距离。对于一个有N个配送点1个配送中心的问题距离矩阵D的大小是 (N1) x (N1)。其中D(i, j)表示从点i到点j的直线飞行距离。这里假设无人机可以在点与点之间直线飞行这是与车辆路径规划最大的不同之一也是无人机配送的优势所在。% 示例生成10个随机配送点范围在[0, 100]的平面内 num_points 10; points 100 * rand(num_points, 2); % 生成10行2列的坐标矩阵 depot [0, 0]; % 配送中心坐标 all_points [depot; points]; % 将所有点合并第1行是配送中心 % 计算距离矩阵 num_all size(all_points, 1); dist_matrix zeros(num_all); for i 1:num_all for j 1:num_all dist_matrix(i, j) norm(all_points(i, :) - all_points(j, :)); end end2.2 目标函数如何评价一条路线的好坏我们的目标是“距离最短”。因此对于任意一条配送路线一个访问所有点的顺序我们称之为一个“解”或“路径”我们需要一个函数来量化它的好坏这个函数就是目标函数。一条路线可以表示为一个序列例如[0, 3, 1, 4, 2, 0]其中0代表配送中心其他数字代表配送点编号。这条路径表示从中心出发 - 访问点3 - 点1 - 点4 - 点2 - 返回中心。那么这条路径的总距离就是总距离 D(0,3) D(3,1) D(1,4) D(4,2) D(2,0)在Matlab中我们可以这样实现目标函数的计算function total_dist calculate_total_distance(route, dist_matrix) % route: 路径序列如 [1, 4, 2, 3, 1]注意这里1代表配送中心索引为1的点 % dist_matrix: 预先计算好的距离矩阵 total_dist 0; for i 1:(length(route)-1) from_node route(i); to_node route(i1); total_dist total_dist dist_matrix(from_node, to_node); end end为什么选择欧几里得距离在城区低空无人机配送的初步规划中直线距离是一个合理且高效的近似。它忽略了起飞/降落、爬升/下降的能耗差异以及风、禁飞区等复杂因素但作为路线结构的初步优化目标是完全可行的。在实际部署前还需要用更精细的能耗模型或考虑三维地形的路径进行二次优化。2.3 问题约束算法必须遵守的规则我们的问题有一个隐含但至关重要的约束每个配送点必须被访问且仅被访问一次。这保证了药品能送到每个需求点且没有重复配送的浪费。在路径序列的表示上这就要求除了代表配送中心的节点通常是序列首尾外其他所有节点在序列中只出现一次。另一个现实约束是无人机的续航距离。在我们的目标函数追求总距离最小化的过程中实际上也间接地在优化续航利用率。但更严谨的做法是在目标函数中增加一个惩罚项如果计算出的某条路径总距离超过了无人机的最大航程就给目标函数值加上一个巨大的惩罚数这样模拟退火算法就会自动淘汰这种不可行的解。3. 模拟退火算法精解从冶金原理到寻优策略模拟退火算法的灵感来源于金属冶炼中的退火过程将材料加热到高温然后缓慢冷却以消除内部应力获得更稳定的晶体结构。算法将组合优化问题中的“解”类比为材料的“状态”将“目标函数值”类比为系统的“能量”通过引入一个不断下降的“温度”参数来控制搜索过程。3.1 算法核心流程与比喻我们可以把寻找最短路径的过程想象成在一个多山的地形解空间里寻找最低的谷底最优解。初始解你随机站在一个地方比如随机生成一条路径。当前能量你所在位置的海拔当前路径的总距离。新解你随机决定往某个方向迈一步对当前路径做一个小的改动产生一条新路径。新能量新位置的海拔新路径的总距离。关键就在于你如何决定是否要迈出这一步如果新位置更低新距离更短那当然要过去这对应算法中的DeltaE 0直接接受新解。如果新位置更高新距离更长在传统的“贪心”算法里这一步会被拒绝你永远只走下坡路。但这很容易导致你困在一个小坑里局部最优解而看不到远处更深的峡谷全局最优解。模拟退火的智慧它说即使新位置更高我也以一定的概率接受它。这个概率取决于两个因素一是“高度差”有多大DeltaE二是当前的“温度”T有多高。高温时即使爬很高的山接受的概率也很大。这相当于在退火初期算法有很强的“跳跃”能力可以在解空间里大范围勘探避免过早陷入局部最优。低温时接受爬山的概率变得极小。这相当于退火后期算法主要在当前位置附近精细搜索趋近于一个稳定的“低点”。这个概率由Metropolis准则决定P exp(-DeltaE / T)。DeltaE是能量差新距离-旧距离T是当前温度。3.2 路径“扰动”策略如何生成新解在TSP问题中如何从当前路径“迈出一步”生成新路径是算法效率的关键。常用的扰动策略有交换Swap随机选择路径中两个非中心点的位置交换它们。例如路径[0, A, B, C, D, 0]交换B和D的位置变成[0, A, D, C, B, 0]。逆转Reverse随机选择路径中一段子序列不包括首尾的中心点将其顺序完全颠倒。例如对[0, A, B, C, D, 0]的B到C段进行逆转得到[0, A, C, B, D, 0]。这种操作在TSP中往往效果很好因为它能较大程度地改变路径结构。插入Insert随机选择一个点将其插入到另一个随机位置。在我的实践中对于中小规模点位数20的配送问题采用两种扰动策略混合的方式效果更佳以一定概率如70%执行“逆转”操作以剩余概率30%执行“交换”操作。这样既保证了搜索的广度又能进行有效的局部调整。function new_route generate_new_route(old_route) % 复制旧路径 new_route old_route; % 注意首尾是配送中心索引1不能动 inner_indices 2:(length(old_route)-1); if rand() 0.7 % 70%概率使用逆转操作 % 随机选择逆转片段的起止索引 idx sort(randperm(length(inner_indices), 2)); start_idx inner_indices(idx(1)); end_idx inner_indices(idx(2)); % 逆转片段 new_route(start_idx:end_idx) fliplr(old_route(start_idx:end_idx)); else % 30%概率使用交换操作 % 随机选择两个不同的内部点索引 swap_idx randperm(length(inner_indices), 2); idx1 inner_indices(swap_idx(1)); idx2 inner_indices(swap_idx(2)); % 交换位置 new_route([idx1, idx2]) new_route([idx2, idx1]); end end3.3 退火计划表控制算法的“火候”退火计划表是模拟退火算法的调度器它决定了温度如何下降以及每个温度下要进行多少次搜索尝试。一个典型的计划表包括初始温度T_init需要足够高使得算法初期几乎能接受任何恶化解。一个经验法则是让初始接受概率在80%以上。可以通过随机采样一些扰动计算平均的DeltaE然后根据T_init -avg_DeltaE / ln(0.8)来估算。温度衰减系数alpha通常取0.8到0.99之间。每次外循环结束时温度更新为T alpha * T。系数越大冷却越慢搜索越细致但耗时也越长。每个温度的迭代次数L也称为Markov链长度。通常与问题规模相关例如设为L 100 * NN为配送点数。保证在每个温度下解空间能得到充分搜索。终止温度T_end或设置最大迭代次数。当温度低于此阈值或连续若干个温度下最优解未改进时算法停止。参数调优心得初始温度和衰减系数对结果影响最大。在我的项目中对于10-15个点的问题T_init1000,alpha0.95,L2000是一个不错的起点。如果发现算法总是很快收敛到一个不太好的解可以尝试提高初始温度或增大衰减系数如0.98让冷却过程更慢。4. Matlab代码实现与逐行解析下面我将结合完整的代码框架详细解释每个部分的作用和实现细节。为了清晰我将代码模块化。4.1 主函数框架与初始化function [best_route, best_dist, history] sa_for_tsp(coords, depot_idx, sa_params) % 基于模拟退火的无人机配送路径规划 % 输入 % coords: (N1 x 2)矩阵所有点的坐标第一行是配送中心 % depot_idx: 配送中心在coords中的索引默认为1 % sa_params: 结构体包含算法参数初始温度、衰减系数等 % 输出 % best_route: 最优路径序列 % best_dist: 最优路径总距离 % history: 记录迭代过程中最优距离的变化用于绘图分析 % 参数设置与初始化 if nargin 2 depot_idx 1; end if nargin 3 sa_params struct(); sa_params.T_init 1000; % 初始温度 sa_params.alpha 0.95; % 温度衰减系数 sa_params.L 2000; % 每个温度迭代次数 sa_params.T_end 1e-8; % 终止温度 sa_params.max_stagnation 50; % 最大停滞迭代次数 end % 计算距离矩阵 num_points size(coords, 1); dist_mat pdist2(coords, coords); % 使用统计工具箱函数计算欧氏距离矩阵 % 生成初始解一个随机的路径排列 % 注意配送中心depot_idx固定在路径首尾 inner_points 1:num_points; inner_points(depot_idx) []; % 移除配送中心 random_route inner_points(randperm(length(inner_points))); % 内部点随机排列 current_route [depot_idx, random_route, depot_idx]; % 构成完整回路 % 计算初始路径距离 current_dist calculate_total_distance(current_route, dist_mat); best_route current_route; best_dist current_dist; % 初始化记录器 T sa_params.T_init; iter 0; stagnation_count 0; history.best_dist [best_dist]; history.temperature [T]; % 退火过程主循环 while (T sa_params.T_end) (stagnation_count sa_params.max_stagnation) for i 1:sa_params.L iter iter 1; % 1. 产生新解 new_route generate_new_route(current_route); new_dist calculate_total_distance(new_route, dist_mat); % 2. 计算能量差 delta_dist new_dist - current_dist; % 3. Metropolis准则判断是否接受新解 if delta_dist 0 % 新解更优直接接受 accept true; else % 新解更差以一定概率接受 accept_prob exp(-delta_dist / T); if rand() accept_prob accept true; else accept false; end end % 4. 更新当前解 if accept current_route new_route; current_dist new_dist; % 5. 更新历史最优解 if current_dist best_dist best_route current_route; best_dist current_dist; stagnation_count 0; % 找到更优解重置停滞计数器 end end end % 记录当前温度下的最优结果 history.best_dist(end1) best_dist; history.temperature(end1) T; % 温度衰减 T sa_params.alpha * T; % 停滞检查如果当前温度下最优解没有提升计数器加1 if history.best_dist(end) history.best_dist(end-1) stagnation_count stagnation_count 1; end end fprintf(算法结束。迭代次数%d 最终温度%.6f 找到最优距离%.4f\n, ... iter, T, best_dist); end关键点解析pdist2函数来自Statistics and Machine Learning Toolbox能高效计算两组点之间的成对距离。如果未安装此工具箱可以用前面展示的双重循环代替。初始解生成通过randperm对内部点进行随机排列这是一个简单有效的策略。更复杂的策略如“最近邻法”可以生成更好的起点但随机起点更能体现模拟退火“不依赖初始解”的优势。停滞计数器这是一个实用的改进。如果连续多个温度周期最优解都没有改善可以提前终止算法节省计算时间。4.2 可视化与结果分析模块算法跑完了我们得看看结果怎么样。可视化不仅能验证结果更是分析和展示的利器。function plot_route_and_history(coords, best_route, history) % 绘制最优路径和算法收敛过程 figure(Position, [100, 100, 1200, 500]); % 子图1路径可视化 subplot(1, 2, 1); hold on; grid on; box on; % 绘制所有点 scatter(coords(:,1), coords(:,2), 70, b, filled); % 高亮配送中心 scatter(coords(1,1), coords(1,2), 120, r, ^, filled); % 绘制路径连线 route_coords coords(best_route, :); plot(route_coords(:,1), route_coords(:,2), k-o, ... LineWidth, 1.5, MarkerSize, 8, MarkerFaceColor, g); % 添加标签 for i 1:size(coords, 1) text(coords(i,1)1, coords(i,2)1, sprintf(%d, i), ... FontSize, 10, FontWeight, bold); end title(sprintf(最优配送路径 (总距离: %.2f), ... calculate_total_distance(best_route, pdist2(coords, coords)))); xlabel(X坐标); ylabel(Y坐标); legend(配送点, 配送中心, 飞行路径, Location, best); axis equal; % 子图2收敛过程可视化 subplot(1, 2, 2); yyaxis left; plot(history.best_dist, b-, LineWidth, 1.5); ylabel(最优距离, Color, b); xlabel(温度下降阶段); title(模拟退火算法收敛过程); grid on; yyaxis right; semilogy(history.temperature, r--, LineWidth, 1.0); ylabel(温度 (对数尺度), Color, r); legend(最优距离, 温度, Location, best); end可视化解读左图路径图直观展示了无人机飞行的顺序。一个“好”的路径应该看起来交叉很少线条相对规整没有明显的长途折返。这是快速评估结果合理性的第一印象。右图收敛图这是诊断算法运行状态的“心电图”。蓝色实线最优距离应该呈现一个总体下降并在后期趋于平稳的趋势。如果曲线早期下降非常陡峭说明初始解很差算法在快速改进如果后期还有频繁的剧烈波动可能意味着终止温度设得过高或衰减过快。红色虚线温度呈指数下降。在温度高的初期最优距离曲线允许有向上的“跳动”接受恶化解随着温度降低曲线逐渐稳定。5. 实战调优与避坑指南理论很美好但把代码跑起来总会遇到各种预期之外的情况。下面分享几个我在项目中实际遇到的坑和解决方法。5.1 算法不收敛或收敛至糟糕解现象算法运行后最优距离曲线几乎是一条平线或者收敛到一个明显不合理的路径交叉严重距离很长。排查与解决检查初始温度温度太低是首要嫌疑。如果初始温度T_init设置过低比如10那么exp(-DeltaE/T)在DeltaE稍大时就会变得极小算法在初期就失去了“爬山”能力迅速陷入最近的局部最优。解决方法按照3.3节提到的方法动态估算一个合适的初始温度。一个简单的测试是手动设置一个较大的T_init如10000再跑一次看曲线前期是否出现波动。检查扰动策略你的generate_new_route函数可能产生的“新解”与“旧解”差异太小或者变化方式无效。例如如果只交换相邻的两个点搜索空间可能受限。解决方法确保扰动能产生足够大的变化混合使用“逆转”和“交换”策略并确保操作的随机索引范围覆盖整个路径除首尾中心点。检查距离矩阵确保距离矩阵计算正确。特别是当坐标值很大时比如经纬度欧氏距离的数值也会很大这会影响DeltaE的量级进而影响接受概率。解决方法可以考虑将坐标归一化到 [0,1] 或 [0,100] 区间或者根据距离矩阵的尺度来调整初始温度。增加迭代次数每个温度下的迭代次数L可能不足。算法在一个温度下还没充分搜索就降温了。解决方法适当增加L例如设为100 * N或200 * N。5.2 运行速度太慢现象当配送点超过30个时程序运行时间显著变长。瓶颈分析与优化目标函数计算是热点在calculate_total_distance函数中我们使用了一个for循环。每次产生新解都要调用它而每次迭代都会产生新解。这是最耗时的部分。优化方法1向量化利用Matlab的向量运算。对于路径route我们可以一次性计算所有相邻点对的距离。function total_dist calculate_total_distance_fast(route, dist_matrix) idx_from route(1:end-1); idx_to route(2:end); % 使用线性索引从距离矩阵中快速提取 linear_indices sub2ind(size(dist_matrix), idx_from, idx_to); total_dist sum(dist_matrix(linear_indices)); end优化方法2增量计算模拟退火中新解通常只由旧解经过微小扰动得到。我们可以只计算路径中发生变化的那部分距离而不必重新计算整条路径。例如对于“逆转”操作只有逆转片段的边界连接发生了变化。这需要更复杂的逻辑但能极大提升速度。减少不必要的计算在计算接受概率exp(-DeltaE/T)时如果DeltaE是负数解变好我们直接接受不需要计算指数。代码中已经做了这个判断。调整退火计划表不一定需要非常慢的冷却alpha0.99和非常长的链L5000。对于很多实际问题一个更“激进”的计划表alpha0.9,L500可能在更短的时间内得到一个可接受的解。这需要在解质量和时间成本之间做权衡。5.3 如何融入实际约束我们目前只考虑了总距离最短。真实的无人机配送还有更多约束载重约束每个配送点的药品重量不同无人机有最大载重限制。时间窗约束某些药品需要在特定时间窗口内送达。续航约束路径总距离必须小于无人机单次飞行的最大航程。融入方法修改目标函数。不再是单纯的最小化距离而是最小化一个“代价”函数这个函数包含距离成本和对违反约束的惩罚。例如处理续航约束function cost calculate_cost(route, dist_matrix, max_range) total_dist calculate_total_distance(route, dist_matrix); penalty 0; if total_dist max_range % 如果超出最大航程施加一个巨大的惩罚 penalty 1e10 * (total_dist - max_range); end cost total_dist penalty; end在模拟退火的主循环中用calculate_cost代替calculate_total_distance来计算能量。这样算法在搜索时会自动避开那些不可行的、超航程的路径。5.4 结果的可重复性与随机性模拟退火算法具有随机性每次运行的结果可能略有不同。这是正常的因为它从随机初始解开始并且随机接受恶化解。如何应对多次运行取最优对于关键任务可以独立运行算法多次如10次然后选择所有运行中找到的最优解。设置随机数种子在开发调试阶段使用rng(123)固定随机数种子可以确保每次运行结果相同便于调试和比较不同参数的效果。关注平均性能在评估算法或参数时不应只看单次运行的最好结果而应看多次运行的平均结果和稳定性。6. 超越基础从单机到多机与动态场景的思考我们的模型是单无人机、静态点、一次性配送。现实场景往往更复杂。多无人机车队路径规划当配送点很多或区域很大时需要多架无人机协同。问题就变成了车辆路径问题VRP。一种常见的思路是“先聚类后路径”聚类阶段根据点的地理分布将它们划分成若干组每组由一架无人机负责。划分的原则可以是组内点距离近且各组的总任务量如点数或总重量均衡。路径规划阶段对每个组分别运行单机TSP算法如我们实现的模拟退火来规划路径。 这仍然可以使用模拟退火但“解”的定义和“扰动”策略会更复杂。解需要包含分组信息和组内路径。扰动策略可能包括将一个点从一个组移到另一个组交换两个点所属的组以及组内路径的优化。动态实时路径规划在实际配送中可能有新的订单随时加入。这就需要算法能动态调整。一种方法是滚动时域优化无人机按照当前计划飞行。当新订单到达时算法立即以“当前无人机位置”为新的起点将“未完成的配送点”“新订单点”作为新的点集重新快速规划一条最优路径。由于对实时性要求高这时可能需要更快的启发式算法或者使用并行计算来加速模拟退火过程。从直线距离到实际航路最终基于直线距离的规划路径需要转化为无人机可执行的实际航路。这需要考虑空域限制避开禁飞区、高楼、高压线。飞行高度不同阶段爬升、巡航、降落的能耗不同。天气因素风会影响飞行速度和能耗。 这通常需要接入地理信息系统GIS和更专业的飞行管理软件。我们的模拟退火算法可以作为上层“任务规划器”输出一个理想的访问顺序然后由下层的“路径规划器”去生成具体的安全航点。最后我想说的是模拟退火算法解决无人机路径规划其魅力不在于它能保证找到数学上的最优解而在于它在计算复杂度和解的质量之间提供了一个优雅且可控的权衡。对于像我们遇到的社区应急配送这类问题它能在几秒到几分钟内给出一个远超人工经验的、切实可用的方案。在Matlab中实现它不仅帮助我们快速验证了想法的可行性其清晰的流程和可视化结果也成为了我们向非技术背景的合作伙伴解释方案价值的绝佳工具。当你看到算法“思考”出的那条蜿蜒但高效的路径在地图上呈现出来时那种将抽象算法转化为实际生产力的满足感正是工程实践的乐趣所在。