简介:本资源面向物流优化、智能算法与无人机系统研究领域的高校师生及工程技术人员,聚焦多无人机协同包裹递送中的区域划分与路径规划两大核心问题。通过融合K-means聚类实现中心辐射式配送区域自动划分,结合遗传算法求解全局最优配送路径,在载荷、续航与地理约束下提升整体配送效率并降低成本。压缩包仅含1个HTML文件(10KB),完整嵌入算法原理说明、流程图解、Matlab代码实现及仿真结果可视化,代码支持地理坐标输入、动态任务调整与多场景路径重优化,可直接运行验证算法有效性。目前已有23人学习下载,适用于算法教学演示、课程设计实践及科研原型快速验证,为无人配送系统建模与仿真提供即开即用的技术支撑。 最近带学生做无人机物流仿真项目时,把以前那套“中心辐射模型 + K-means聚类 + 遗传算法优化”的Matlab代码翻出来重写了一遍。这三个词组合在一起,几乎是无人机配送路径规划领域的“老三样”,论文里一抓一大把,但真到了自己上手从零写到能跑、能调、能可视化,我才发现里面全是细节:K-means的K到底怎么定?遗传算法的染色体怎么编码才能表达“多架无人机同时配送”?适应度函数里的续航惩罚项该给多大权重?这些问题在大多数参考代码里都不会直说,只有踩过坑才明白。
这篇文章就把我重构三版之后稳定运行的完整思路整理出来。从需求建模、聚类选址细节、GA路径规划的编码与算子设计,到Matlab调试过程里的高频报错和避坑方法,一次性讲清楚。代码结构也会拆开说明,适合正在做物流无人机、智能优化算法课程设计,或者想用K-means + 遗传算法解决组合优化问题的读者对照参考。
1. 中心辐射:为什么无人机送货要先“聚”再“送”
中心辐射模型不是无人机行业发明的,它最早大规模用在航空客运网络里:几个大型枢纽机场之间用主干航线连接,乘客再从枢纽坐支线飞机去中小城市。这套思路放在快递物流里就是经典的hub-and-spoke结构:包裹先集中到区域分拣中心,再分拨到末端网点,最后送到收件人手里。无人机配送之所以天然适配这种结构,核心原因只有一个:电池续航扛不住点对点长距离飞行。
1.1 无人机配送的续航悖论
假设一个城市里有100个客户点,仓库在城郊,如果让一架无人机从仓库出发挨个配送,路径总长度很可能超过100公里,而市面上多数多旋翼物流无人机的有效航程也就15到30公里。哪怕任务分配给多架无人机,每架都从仓库直飞客户点再返回,仍然会出现大量长距离冗余飞行。实际项目中我们通常会在配送区域内设置若干个起降站点,无人机从站点出发,只负责服务站点周边半径3到5公里内的客户。这样虽然多了一次“中转”,但每架无人机的飞行距离大幅缩短,换电、装载、充电都可以在站点完成,整体效率反而更高。
1.2 中心辐射网络的分层逻辑
这套模型在我们的Matlab代码里是这样分层的:
- 第一层:大型分拨中心。相当于“根节点”,负责接收所有包裹,并向各站点补货。
- 第二层:配送站点。通过K-means聚类生成,每个站点服务一片客户簇,选址目标是最小化簇内客户到站点的平均距离。
- 第三层:末端客户点。即需要交付包裹的具体坐标,它们被聚到某个站点名下,再由站点派出无人机逐一访问。
那么整个路径规划问题就拆解为两个子问题:一是站点选址与客户分配(K-means负责),二是每个站点内多无人机的访问顺序(遗传算法负责)。实际代码中我们把这两个问题串成一个Pipeline:先用kmeans函数把客户点聚成K类,类中心作为初始站点;然后把每类客户看作一个TSP(旅行商问题)实例,用GA规划无人机访问顺序。
1.3 模型假设与输入数据设计
为了让问题可解,我采用了以下假设条件,这些也是大多数同类仿真论文的默认设置:
- 所有客户点坐标已知,无人机飞行路径为直线(等效为欧氏距离)。
- 每架无人机运载能力相同,单次最多携带一定数量的包裹后必须返回站点补货。
- 无人机飞行速度恒定,续航里程有上限,超过上限的路线会被惩罚。
- 站点位置可自由选择,暂不考虑空域限制、楼顶平台、信号覆盖等真实条件。
这些假设在第一章先定下来,后面写遗传算法适应度函数时就不至于反复推翻重来。实际代码中我定义了一个全局结构体param,把无人机数量、最大航程、载重上限、客户坐标等放进去,后续所有脚本都从param读取,避免硬编码。
2. K-means选址:聚类中心不等于起降点
很多初学K-means的人以为,把客户坐标丢进kmeans函数,输出聚类中心就可以直接当无人机起降点了。实际使用中完全不是这么回事,这里有两个很容易被忽略的点:一是K值怎么确定,二是聚类中心是几何中心,不代表现实中可以落地的位置。
2.1 K-means在这一步到底在算什么
我们的目标是:把N个客户点划分到K个站点,使得每个客户到它所属站点的距离之和最小。用数学语言描述就是:
[ \min \sum_{i=1}^{K} \sum_{x \in C_i} |x - \mu_i|^2 ]
其中(C_i)是第i个簇,(\mu_i)是第i个簇的中心。这个目标函数对应的是“簇内平方和”,也就是Matlab里kmeans默认优化的目标。在无人机场景里,站点是终点也是起点,这一步的核心价值是把“服务范围”圈定出来,让每个站点的客户数量大致均衡,从而避免某架无人机累死、另一架闲死。
2.2 K值选取:肘部法则里的数学
K值不是越大越好。K值太小,站点覆盖范围太大,无人机单程飞不到;K值太大,站点数量多,基础设施成本过高,而且每个簇里客户太少,遗传算法容易退化。项目里我一般用肘部法则辅助判断。
Matlab中可以用以下代码绘制簇内平方和随K值变化的曲线:
rng(42); X = [customers.x, customers.y]; K_list = 1:10; sumd_list = zeros(size(K_list)); for i = 1:length(K_list) [~, ~, sumd] = kmeans(X, K_list(i), 'Replicates', 10, 'MaxIter', 300); sumd_list(i) = sum(sumd); end plot(K_list, sumd_list, 'o-'); xlabel('K'); ylabel('Sum of within-cluster distances');这里我刻意使用了Replicates参数。kmeans默认只跑一次,初始中心完全随机,结果方差很大。设置Replicates=10表示从10组不同的初始中心出发,返回目标函数最小的一次结果,这样结果更稳定。
肘部法则的“肘”出现在曲线下降趋势明显变缓的位置。比如我测试某个城市数据时,K从1到4下降非常陡峭,从4到10逐渐平缓,那个拐点K=4就是合理的选择。当然,把K取自变量还要结合实际业务约束,比如手里只有5架无人机,那K就不能超过5,否则站点多但无人机调度不过来。
2.3 Matlab聚类实现与初始中心敏感性
K-means对初始中心敏感,这是老话题了。Matlab的kmeans虽然内置了k-means++策略(Start参数默认为'plus'),但依然可能陷入局部最优。常见的优化方式是多组初始中心并行跑,然后取目标函数最小的解,这正是前面代码里Replicates=10的作用。
另一个容易被忽略的细节是坐标尺度对聚类结果的影响。如果客户点分布在经度、纬度这种量纲不一致的坐标上,直接聚会产生严重偏差。我的做法是把经纬度投影成平面直角坐标,常见选择是Web Mercator或者高斯-克吕格投影,然后再做标准化处理:
X = [customers.x, customers.y]; X = zscore(X); % 标准化注意:zscore会把坐标中心化、缩放成单位方差,这样聚出来的簇结构更稳健,但聚类中心返回的是标准化后的坐标,使用时必须反向还原才能得到真实坐标。
2.4 从簇心到真实起降点的修正
这一步是我在实际项目中最想强调的经验:K-means给出的聚类中心是一个“几何质心”,但在真实城市场景里,这个位置可能位于马路中央、湖面上或者建筑群正中间,无人机根本无法起降。所以我会在聚类之后增加一个“最近可用点匹配”步骤:提前准备一份“候选起降点”列表,例如城市中允许无人机起降的屋顶平台、空地和站场,然后对每个聚类中心,在候选点列表里找距离最近且未被占用的点作为实际站点。
这个逻辑在Matlab里用dsearchn就能实现:
candidate_sites = [sites.x, sites.y]; centroids = C; selected_sites = zeros(K,2); for i = 1:K nearest_idx = dsearchn(candidate_sites, centroids(i,:)); selected_sites(i,:) = candidate_sites(nearest_idx,:); end这一步虽然属于工程细节,但在论文里写“站点选址采用K-means聚类与最近可用点匹配相结合”会让评审觉得你考虑到了真实约束,而不只是跑了个数学模型。
3. 遗传算法跑路线:编码、适应度、算子三件套
站点分好了,接下来每个站点内的客户点顺序就交给遗传算法。为什么不用精确算法?因为这是个NP难问题,站点内客户一旦超过15个,暴力枚举就变得不现实。遗传算法虽然没有理论上的最优解保证,但胜在通用、容易实现,而且和Matlab的矩阵运算天然契合。
3.1 为什么不用精确算法而选GA
对于单个站点10个客户,TSP的路径数量是10! = 3,628,800条,计算机还可以枚举;但站点内客户到20个时,20!已经超出任何民用机器的计算能力。精确算法里像分支定界法,能在小规模问题下找到最优解,但问题规模稍微变大就会指数爆炸。遗传算法不一样,它不追求枚举全部路径,而是从一组随机初始解出发,通过选择、交叉、变异一步步逼近较优解,虽然无法保证找到全局最优,但在工程上“够用”且计算时间可接受。
3.2 多无人机任务分配的双染色体编码
单站点只有一架无人机时,染色体用一条整数排列表示客户访问顺序即可。但实际场景往往是多个站点、每个站点有多架无人机,怎么设计编码就各有门道了。我采用了一种双染色体编码方案:
- 染色体1(任务分配染色体):长度等于客户数,每个基因位存储该客户被分配给哪个站点。
- 染色体2(访问顺序染色体):长度等于客户数,存储一个全局客户序号的排列。
解码时,先根据染色体1把客户按站点分组,再按染色体2中各客户的相对顺序进行排序,从而得到每个站点内各无人机的访问序列。这个方案的好处是:交叉和变异可以分别作用在两条染色体上,既保证了站点分配的变化,又保证了访问顺序的多样性。
举个具体例子:假设有6个客户,2个站点。染色体1为[1,2,1,2,2,1],染色体2为[4,6,2,1,5,3]。那么客户1、3、6属于站点1,客户2、4、5属于站点2。按染色体2的顺序排序,站点1的访问顺序是6→1→3,站点2的访问顺序是2→4→5。
3.3 适应度函数:距离只是底线,约束才是关键
适应度函数直接决定GA搜索方向,我只把“总飞行距离短”作为主要目标还不够,必须把续航约束、载重约束一起纳入,否则算法会一厢情愿地规划出“理论距离最短但无人机根本飞不到”的路线。
设某条完整路径的总距离为(D),超出最大续航的距离为(E_{exceed}),超出载重量的包裹数为(P_{exceed}),那么适应度函数定义为:
[ fitness = D + \lambda_1 \cdot E_{exceed} + \lambda_2 \cdot P_{exceed} ]
其中(\lambda_1)和(\lambda_2)是惩罚系数。代码中我会这样处理:
function fit = compute_fitness(chrom1, chrom2, param) % 解码 routes = decode(chrom1, chrom2, param); total_distance = 0; penalty = 0; for i = 1:length(routes) route = routes{i}; if isempty(route) continue; end dist = route_distance(route, param.distance_matrix); total_distance = total_distance + dist; if dist > param.max_range penalty = penalty + (dist - param.max_range); % 超出部分线性惩罚 end % 载重约束 package_num = param.package_demand(route); if sum(package_num) > param.max_payload penalty = penalty + (sum(package_num) - param.max_payload) * 10; end end fit = total_distance + penalty; end这里有个关键点:惩罚系数的量纲要和距离保持兼容。param.max_payload的单位是“个包裹”,如果罚款系数设得太大,比如1000,那么算法会完全无视距离而只想着别超载;如果设得太小,比如0.1,那么超载的路线会被当作好路线保留下来。经验上,我会先统计初始种群中的距离均值,再把超载惩罚系数设置成距离均值的1/10左右,比如距离均值是1000,那超载1个包裹罚100。
实际上在测试中,我还会对超出续航的部分采取“硬惩罚”和“软惩罚”结合:超出1%以内只加一个固定惩罚值,超出超过1%则返回一个极大值,直接淘汰该个体。这样处理的好处是避免算法在续航边界反复试探。具体阈值要根据无人机型号参数来设,我这里用的是10%。
3.4 交叉与变异:针对顺序编码的算子设计
遗传算法的算子必须和编码方式匹配。对顺序编码的染色体2,不能直接使用单点交叉,否则很容易产生重复客户点、丢失部分客户点的非法解。我常用的顺序交叉(Ordered Crossover,OX)思路是:
- 随机选择两个交叉点;
- 将父代1两交叉点之间的基因段保留到子代1对应位置;
- 子代1剩余位置按照父代2中非重复基因的相对顺序依次填充;
- 父代2同理生成子代2。
这种操作保证了每一个客户点只出现一次,适合路径顺序类编码。MATLAB里我写了一个简单的OX交叉函数,核心代码如下:
function child = ox_crossover(p1, p2) n = length(p1); point1 = randi(n - 1); point2 = point1 + randi(n - point1); child = zeros(1, n); child(point1:point2) = p1(point1:point2); % 从p2中提取非重复基因 remaining = p2(~ismember(p2, child)); idx = 1; for i = 1:n if child(i) == 0 child(i) = remaining(idx); idx = idx + 1; end end end变异算子我用的是“两点交换+片段倒置”混合变异:以一定概率随机交换两个基因位,或者随机选一段基因倒置。两类变异以概率各半触发,这样既能带来扰动,又能保持路径的局部结构。
对染色体1,也就是站点分配基因,我用的是多点变异:随机挑若干基因位,把其值改成另一个站点的编号。这个变异力度不能太大,否则每个站点客户量严重失衡,导致优化过程剧烈震荡。
3.5 参数表与调参经验
遗传算法参数我常用以下初始值,实际调优建议在此基础上微调:
| 参数 | 取值 | 说明 |
|---|---|---|
| 种群规模 | 100 | 客户多时可增加到200 |
| 最大迭代次数 | 200 | 看收敛曲线决定 |
| 交叉概率 | 0.85 | 不宜超过0.9,避免随机搜索 |
| 变异概率 | 0.15 | 太大容易破坏优秀解 |
| 精英保留数 | 2 | 每代直接复制到下一代 |
| 惩罚系数λ1 | 10 | 按续航超出距离比例调 |
| 惩罚系数λ2 | 100 | 按平均距离的1/10设定 |
调参的经验是:先固定变异概率和交叉概率,只调惩罚系数;等路线不再出现超续航的个体之后,再调种群规模和迭代次数。不要一上来就同时改所有参数,否则你根本不知道是哪个参数让结果变好的。
4. 调试与避坑:Matlab实现中的高频问题
这部分是我实际调试代码过程中遇到的最多、最典型的坑。每一个我都在代码注释里标了醒目的警示符号,这里专门整理出来。
4.1 坐标格式带来的距离失真
第一版代码里我直接用经纬度做欧氏距离运算,导致结果完全失真。经度1度在赤道附近大约是111公里,但在高纬度地区可能只有几十公里;纬度1度则基本恒定在111公里左右。如果直接把经纬度丢进sqrt((x1-x2)^2+(y1-y2)^2),在哈尔滨和在三亚算出来的同度数距离会差很多。
解决办法是先用deg2km函数或者等距投影转成平面坐标。Matlab中可以用lldistkm函数(需要Mapping Toolbox)或者自定义投影函数:
% 简单等距投影,适用于城市尺度 lat0 = mean(customers.lat); lon0 = mean(customers.lon); x = (customers.lon - lon0) * 111.32 * cosd(lat0); y = (customers.lat - lat0) * 110.94;如果你的客户坐标只是仿真生成的二维坐标,这一步可以省略,但代码里我还是保留了投影函数入口,为的是以后接入真实试验数据时不用改大框架。
4.2 K-means空簇问题
kmeans在某些随机初始中心下可能产生空簇,也就是某个聚类中心没有分配到任何客户。这种情况下遗传算法里对应站点的客户列表为空,计算适应度时会出现NaN或索引越界报错。
我在实际代码中的处理方式是:在kmeans之后加一个检查,如果任何簇为空,重新运行聚类并递增Replicates。还有一个更实用的思路:先给每个簇强制分配距离该簇中心最近的客户,然后再跑一轮迭代修正簇心。这有点类似带种子约束的K-means变体,在Matlab中实现也不复杂。
4.3 GA早熟:如何判断和干预
遗传算法最常见的现象是迭代到50代左右,最优适应度就不再下降了,但解明显还不是最优解。这就是早熟收敛。我通常用两个指标判断:
- 最优个体适应度连续50代没有变化;
- 种群平均适应度和最优适应度差值小于某个阈值,说明种群多样性丧失。
遇到早熟,我不会盲目加大变异概率,而是采用“自适应变异”:当检测到最优值停滞时,把变异概率临时提升到0.3,维持若干代后再降回来。另外,还可以引入“移民策略”,每20代随机生成一批新个体加入种群,替代最差的一些个体。代码里我用一个if mod(gen, 20) == 0的触发器来做,简单有效。
4.4 惩罚项权重,或无穷大警告
适应度值出现Inf或者NaN时,大概率是惩罚系数设置不合理。在较早期版本中,我直接用fit = 1/(total_distance + penalty)作为适应度,当total_distance + penalty为0时就会爆炸。后来我改成直接使用total_distance + penalty作为目标函数,遗传算法内以“最小化”为方向,只在选择操作时把目标值越小视为适应度越高。
如果你坚持用“适应度越大越好”的形式,一个更稳妥的转换是:
fit = 1 / (1 + total_distance + penalty);这样能保证分母恒大于0,不会因为惩罚过大出现除零问题。
4.5 结果可复现性
Matlab的rand和randi每次运行结果不同,导致两次跑出的路径差异很大,写报告时很尴尬。我在所有随机操作前统一加了一行:
rng(2025);这样每次跑出的初始种群、聚类中心、交叉变异位置都完全一致,方便对比不同参数对结果的影响。注意rng(2025)的位置要放在所有随机操作之前,最好放在脚本第一行,因为kmeans内部也用了随机数生成器,如果你的rng放在kmeans之后才设置,控制不了聚类的随机性。
可视化部分我也做了个简单处理:用不同颜色区分不同站点服务的客户簇,用带箭头的线条表示无人机飞行轨迹。Matlab里可以用quiver画箭头,也可以用plot加annotation实现。实跑效果来看,一张图能把站点、客户点、航线全部展示出来,比单纯输出数字直观得多。
5. 从代码到天空:这套方案的边界与延伸方向
虽然这套仿真跑起来效果不错,但我心里清楚它离真实可部署的无人机配送系统还有一段距离。最明显的问题是第一章节我就提到的“直线距离假设”——城市里无人机虽然不受地面道路限制,但受空域管制、高楼遮挡、信号干扰、气象条件等多重因素影响,实际飞行轨迹永远不可能是标准直线。
5.1 当前模型的简化假设
除了直线距离,还有几个容易忽视的简化:
- 电池能耗没有和载重、风速关联,真实无人机带2kg包裹和带5kg包裹的续航完全不同;
- 客户时间窗没有考虑,所有订单都要求“尽快送达”,没有优先级;
- 站点容量没有约束,假设所有站点都能无限停放无人机、无限补货;
- 无人机之间的冲突规避没有考虑,多机同时飞行时可能存在航线交叉。
这些假设决定了这套代码适合做算法验证和技术预研,但不适合直接拿去做真实空域规划。如果真想落地,至少要把电池能耗模型和空域栅格化建模加进去。
5.2 从静态到动态
当前流程是“先聚类,再规划”,属于静态一次性优化。实际配送场景是订单不断到达的,所以更合理的方式是多阶段滚动优化:每当一批新订单到达,重新执行聚类和路径规划,保留已经起飞的无人机任务不变,只重规划未分配订单。这个方向上可以用Matlab的Timer或实时脚本实现,但计算压力会明显增加。
5.3 可能的扩展
这套代码框架稍微改一改就能扩展到很多方向。比如把目标函数从“总距离最短”改为“总能耗最低”,就需要给每个客户点增加“重量需求”属性,并在适应度函数中引入载荷相关的能耗模型;又比如加入时间窗约束,就需要在适应度函数里加入时间惩罚项,同时把遗传算法的解码过程改成带时间轴的任务插空排序。我后续打算把多站点之间的补货路径也纳入优化,让无人机既送末端客户,也承担站点间调拨任务,这样整个中心辐射模型就更完整了。
如果只是想应付课程设计或论文仿真,这套方案已经完全够用;如果想进一步挑战自己,可以顺着“多目标优化(NSGA-II)+ 无人机调度”的方向继续深挖。毕竟K-means和遗传算法都只是工具箱,解决实际问题时真正值钱的是对约束的建模能力和对算法行为的理解。
本文还有配套的精品资源,点击获取