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的大小是 (N+1) x (N+1)。其中
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(i+1); 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 * N(N为配送点数)。保证在每个温度下,解空间能得到充分搜索。 - 终止温度(T_end):或设置最大迭代次数。当温度低于此阈值,或连续若干个温度下最优解未改进时,算法停止。
参数调优心得:初始温度和衰减系数对结果影响最大。在我的项目中,对于10-15个点的问题,T_init=1000,alpha=0.95,L=2000是一个不错的起点。如果发现算法总是很快收敛到一个不太好的解,可以尝试提高初始温度或增大衰减系数(如0.98),让冷却过程更慢。
4. Matlab代码实现与逐行解析
下面,我将结合完整的代码框架,详细解释每个部分的作用和实现细节。为了清晰,我将代码模块化。
4.1 主函数框架与初始化
function [best_route, best_dist, history] = sa_for_tsp(coords, depot_idx, sa_params) % 基于模拟退火的无人机配送路径规划 % 输入: % coords: (N+1 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(end+1) = best_dist; history.temperature(end+1) = 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(增量计算):模拟退火中,新解通常只由旧解经过微小扰动得到。我们可以只计算路径中发生变化的那部分距离,而不必重新计算整条路径。例如,对于“逆转”操作,只有逆转片段的边界连接发生了变化。这需要更复杂的逻辑,但能极大提升速度。
- 优化方法1(向量化):利用Matlab的向量运算。对于路径
- 减少不必要的计算:在计算接受概率
exp(-DeltaE/T)时,如果DeltaE是负数(解变好),我们直接接受,不需要计算指数。代码中已经做了这个判断。 - 调整退火计划表:不一定需要非常慢的冷却(
alpha=0.99)和非常长的链(L=5000)。对于很多实际问题,一个更“激进”的计划表(alpha=0.9,L=500)可能在更短的时间内得到一个可接受的解。这需要在解质量和时间成本之间做权衡。
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中实现它,不仅帮助我们快速验证了想法的可行性,其清晰的流程和可视化结果,也成为了我们向非技术背景的合作伙伴解释方案价值的绝佳工具。当你看到算法“思考”出的那条蜿蜒但高效的路径在地图上呈现出来时,那种将抽象算法转化为实际生产力的满足感,正是工程实践的乐趣所在。