💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
🎁完整资源、论文复现、期刊合作、论文辅导及科研仿真定制事宜点击：
👉👉👉本文完整资源下载
⛳️座右铭：行百里者，半于九十。
⛳️赠与读者
👨💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能解答你胸中升起的一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。
或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎
💥第一部分——内容介绍
摘要
本文系统阐述了一种基于灰狼优化算法（Grey Wolf Optimizer, GWO）及其多种群改进版本（MP-GWO）的多无人机协同航迹规划方法。针对复杂战场环境中多架无人机的时间协同、威胁规避、运动学约束及防碰撞等核心问题，构建了包含航程代价、高度惩罚、雷达与通用威胁、时间协同偏差及碰撞次数在内的五维加权目标函数。详细阐述了从决策变量编码到约束处理的完整数学逻辑，分析了标准GWO基于等级制度的位置更新机制与MP-GWO依托子目标聚类实现的多样性保持策略。该实现具备二维与三维空间的自适应能力，为航迹规划算法的工程应用与对比研究提供了可复现的基准平台。
1. 引言
多无人机协同航迹规划是任务规划系统的核心环节，其目标是在满足各型无人机平台固有飞行性能约束的前提下，为编队规划出能够有效规避敌方雷达探测区域与火力杀伤区、并在严格的时间窗口内同步抵达各自目标点的空间飞行路径。该问题本质上是一个高维度、强约束、多目标耦合的非线性优化难题，传统数值优化方法往往难以在可接受的时间内求得满意解。群体智能算法因其对问题数学性质要求宽松、天然具备并行搜索能力等优势，近年来在该领域获得了广泛关注。本文以灰狼优化算法作为核心求解器，并引入多种群协同进化策略，构建了一套层次清晰、接口完备的算法实现框架。
2. 协同航迹规划数学模型
2.1 决策空间与航迹表示
设无人机编队包含若干架异构平台，其具体数量由任务场景动态决定。仿真维度依据任务需求设定为二维平面或三维空间。对于编队中的任意一架无人机，系统均给定其明确的起点坐标与目标点坐标，并在两点之间设置若干个待优化的中间导航节点。这些节点在空间中的具体位置并非固定，而是作为决策变量交由优化算法进行搜索确定。
完整航迹由起点、所有中间导航点以及终点依次连接而成的折线段构成。除了空间坐标外，每架无人机还需分配一个代表其巡航速度的标量值。将所有无人机的全部导航点坐标与速度值按照既定的顺序拼接，即可构成一个长度可观的全局决策向量。该向量的总维度由无人机架数、每架飞机的导航点数量以及空间维度的乘积决定，体现了问题的复杂程度。
2.2 约束条件体系
在实际突防任务中，优化得到的决策向量必须严格满足一系列物理与战术约束条件，任何一项指标的违背都意味着航迹不可行。这些约束涵盖了多个层面：
在运动性能方面，每架无人机的飞行速度必须处于其动力系统所能提供的区间之内；水平面内的转弯角度以及垂直面内的俯仰角度（针对三维场景）均不能超过机体结构所允许的最大机动过载，以防止失速或结构损伤。
在空间使用方面，所有导航点的高度（三维场景下）必须限定在低空突防的安全高度区间内，同时平面位置也需处于任务空域的边界范围以内，避免飞出作战区域。
在航迹几何质量方面，任意两个相邻导航点之间的距离不得小于设定的最小航段长度，这是为了确保航迹具有一定的宏观指向性，避免出现无意义的局部锯齿抖动；同时，每架无人机的总飞行里程不能超出由起点到终点直线距离所决定的最大容许航程，这一限制源于燃油载量的客观约束。
在战场生存方面，任何一段连接相邻导航点的航迹线段都不允许穿越事先标定的各类威胁区域，包括雷达探测锥区和防空武器杀伤区。此外，对于编队内的任意两架无人机，在飞行过程中的任一时刻，它们之间的空间距离都必须维持在安全间隔阈值以上，以防止空中相撞事故的发生。
2.3 约束的柔性修复策略
若在进化过程中对每一个个体都强制其必须同时满足上述所有硬性约束，将会导致可行域的极度碎片化，严重制约算法的搜索效率。为此，本系统采用了一种后验修复策略。当系统检测到某个导航点导致了转角过大、穿越威胁或航段过短等违规情况时，该点会被判定为“问题点”。此时，算法不会简单地丢弃该个体，而是取该点前后相邻两个节点的空间位置进行算术平均，用计算得到的均值点替换原有的问题点。对于靠近航迹首端或末端的问题点，则分别使用固定的起点或终点作为参考依据进行修复。完成上述几何修正后，系统还会对全局决策向量的所有分量执行一次边界裁剪操作，将所有越界的数值强制拉回到预设的上下限边界上。尽管这种启发式的修复手段带有一定的经验色彩，但在实际计算中能够有效维持种群的可行性比率，并引导航迹曲线自然地趋向平滑。
3. 综合评价目标函数设计
航迹规划的优劣难以用单一指标衡量。本系统将飞行效能分解为五个各具物理意义的子目标，并采用经典的线性加权求和法将多目标问题转化为单目标适应度值，以便于群体智能算法的迭代寻优。
第一个子目标是航程代价，它代表了燃油消耗的等效度量，通过将所有无人机的实际飞行距离累加，并除以预先计算的最大容许总航程进行归一化处理，使得该指标的量纲不依赖于具体场景的尺度。
第二个子目标是高度偏离惩罚，仅在三维空间规划中生效。系统对每个导航点的实际高度进行检查，若其超出安全高度带的上限或下限，则根据超出的幅度施加相应的线性惩罚值。
第三个子目标是综合威胁暴露代价，这是战生存性的核心量化指标。系统在模型中区分了两种不同物理特性的威胁源：对于雷达探测威胁，其信号强度随距离衰减极快，因此代价与距离的四次方成反比；而对于高炮火力或恶劣气象等威胁，代价则与距离的一次方成反比。优化算法将迫使航迹节点尽可能地远离这些威胁区域的中心。
第四个子目标是时间协同偏差，这也是多机协同区别于单机规划的关键所在。系统会预先给定一个期望的协同到达时刻，每架无人机根据各自的总航程和分配到的速度计算出实际飞行用时。若该用时落在由速度上下限决定的可行时间窗口之外，则计算其与协同时刻之间的绝对偏差并施加惩罚，从而迫使各无人机的预计到达时间趋向一致。
第五个子目标是碰撞风险代价，即通过前述的多机时空碰撞检测流程统计得到的全局碰撞总次数。该指标直接以次数计量，用以引导算法在空间上错开各机的航迹。
最终，将上述五个经过合理缩放和加权的子目标值进行求和，所得标量值即为综合评价适应度，该值越小代表航迹方案的综合性能越优越。
4. 标准灰狼优化算法核心机制
灰狼优化算法通过模仿自然界狼群严格的社会支配等级与集体狩猎行为来实现对最优解的搜索。在每一代进化过程中，算法将当前种群中适应度最优的个体尊为头狼，适应度次优和第三优的个体分别尊为二狼和三狼，其余个体统称为普通狼。
在捕食行为的数学模拟中，算法首先定义了一个随进化代数从初始值线性递减至零的关键收敛参数。该参数直接调控着搜索步长与勘探开发之间的平衡。对于种群中的每一个普通个体，算法会分别计算它相对于上述三只最优秀头狼的位置矢量差，并通过引入两组由均匀分布随机数生成的系数向量来赋予这些差值以随机扰动。随后，每个个体都会获得三个指向不同领袖方向的引导位置。最终，该个体的下一代位置被确定为这三个引导位置的加权算术平均中心点。当收敛参数的绝对值较大时，个体的移动步幅会显著增加，促使算法在全局范围内搜索未知区域；而当该参数收缩至较小值时，个体的运动将被约束在领袖个体附近的局部空间内，实现了对优质解邻域的精细开采。
5. 多种群灰狼优化算法的改进策略
为应对标准灰狼算法在高维复杂航迹规划场景中容易因种群同质化而过早陷入局部最优的缺陷，多种群版本引入了一种基于子目标性能差异的聚类分裂机制。
传统的单种群算法仅依赖总适应度值进行排序，这将不可避免地丢失关于个体在不同子目标上表现优劣的宝贵信息。多种群策略将每个个体的五维子目标向量作为分类的特征标签。具体操作时，系统首先计算当前全体种群中每一个体的子目标向量，形成一个二维性能矩阵。随后，针对矩阵的每一行，即每一个独立的子目标维度，分别进行由优到劣的排序。最后，依据这些排序序号将种群在行方向上进行分段切割，等量地划分为多个子种群。
经过这种划分方式，第一个子种群将在航程指标上具有相对优势，第二个子种群将在高度控制上表现突出，以此类推。这种设计保证了各个子种群能够沿着不同的性能偏好方向进行专业化进化，从而在整体上维持了种群在多重相互冲突的目标之间的多样性覆盖。
在每一代的进化过程中，各个子种群完全独立地维护属于本群的专属头狼、二狼和三狼，并行执行标准灰狼算法的位置更新公式。各子种群之间互不干扰，直到当前迭代步结束再统一进行合并。最终呈现的全局最优解是所有这些子种群头狼中总适应度最优的个体，而输出的收敛曲线则是对各子种群头狼适应度取算术平均值，此举有效平滑了随机震荡，更能反映出算法整体的宏观收敛态势。
6. 关键辅助模块的工程实现
为了保障优化算法输出的航迹具有实际工程意义，系统配备了两大精密的后端检测模块。
第一个是威胁穿越检测模块。该模块负责判断任意一段连接两个相邻导航点的直线段是否与三维空间中的球体威胁区（或二维平面中的圆形威胁区）发生干涉。检测逻辑结合了点与球心的距离判断以及线段到球心的最小距离几何判据。若线段的某个端点已落入威胁区内则直接判定为穿越；否则计算球心到该线段的垂直距离，若该距离小于威胁半径且垂足落在线段的有效范围之内，则同样判定为穿越。这一双重判据确保了威胁规避检测的数学严谨性。
第二个是多机时空碰撞检测模块。由于各无人机的速度设定不同，它们经过同一地理坐标点的时间各异，因此静态航迹相交并不等同于实际碰撞。该模块采用了基于时间轴插值的精细检测方法。对于目标无人机的某一个导航点及其对应的到达时刻，系统会在另一架无人机的航迹上通过线性插值算法找到该架飞机在同一时刻所处的空间位置。若两者之间的欧氏距离小于预设的安全间隔，则计数器累计一次碰撞事件。这种方法充分考虑了速度差异带来的时变空间位置关系，能够真实反映编队飞行的安全裕度。
7. 实验配置框架与运行指引
本系统在主入口脚本中通过一个简洁的枚举变量供用户切换两种优化算法。默认推荐的标准参数配置兼顾了计算耗时与解的质量，种群规模设定为数十个个体，最大进化代数设定为百余代，足以应对常规规模的突防场景。
系统内置了多种预设的典型作战场景供研究人员直接调用。这些场景覆盖了二维大范围平面机动突防、三维低空地形跟随突防以及高密度威胁区域的饱和突防等多种战术想定，各自配置了不同数量的无人机编队、差异化的导航点个数以及空间分布各异的威胁群组。
此外，系统具备完善的参数自检功能。在正式进入迭代寻优之前，入口程序会自动校验仿真维度的一致性、威胁坐标与空间维度的匹配性、各无人机约束矩阵的行数对齐情况，尤其是严格检查每个无人机设置的导航点数量是否超出了由最小航段长度和最大航程所决定的物理容纳上限。一旦发现参数配置存在逻辑冲突或数值越界，程序将立即抛出明确的错误提示并终止运行，从而有效避免了无效计算对时间的浪费。
8. 总结与展望
本文深入剖析了一套融合了标准灰狼优化算法与多种群灰狼优化算法的多无人机协同航迹规划代码框架。在模型层面，通过将复杂的飞行动力学与战场威胁约束转化为可后验修复的柔性编码策略，并将多目标权衡问题汇聚为线性加权的标量适应度函数，该实现能够在合理的时间开销内生成几何平滑、规避有效且时间同步的可行空间路径。在算法层面，多种群策略通过基于子目标偏好的聚类划分，显著增强了个体在竞争目标之间的差异化探索能力，有效抑制了标准算法在迭代后期易出现的早熟收敛现象。
整个软件框架具有高度模块化和接口灵活的特点，既可作为研究群体智能算法性能比较的验证基底，也可通过简单的场景配置文件快速适配不同战术需求的工程任务。未来的迭代改进方向可以考虑在该框架内引入严格意义上的帕累托非支配排序与拥挤度距离机制，彻底摒弃预先设定的加权因子，从而实现真正的无权重多目标协同优化。同时，针对超大规模无人机集群场景，可研究将聚类进化策略与分布式计算架构进一步融合，以突破单机计算资源的性能瓶颈。
📚第二部分——运行结果
部分代码：
function solution = MP_GWO(UAV, SearchAgents, Max_iter)
%MP_GWO 多种群灰狼优化算法
%Multi Population Gray Wolf Optimization
% 超参数
g = 50; % 动态更新加权系数
% 算法初始化
[WolfPops, ~] = PopsInit(UAV, SearchAgents, false); % 随机生成 初始狼群
ClassPops = PopsCluster(WolfPops, UAV); % 进行初始聚类（用来获得 k 值）
dim = WolfPops.PosDim; % 状态变量维度
cSearchAgents = ClassPops.SearchAgents; % 搜索智能体个数（子种群数量）
SearchAgents = cSearchAgents * ClassPops.k; % 对所有智能体个数进行修正（k的整数倍）
WolfPops.Pos = WolfPops.Pos(1:SearchAgents, :); % 对种群进行修正
% 报错
if cSearchAgents < 4
error('搜索智能体个数过少')
end
% 初始化解
Alpha_pos = zeros(ClassPops.k, dim); % α解
Alpha_score = 1 ./ zeros(ClassPops.k, 1); % α解适应度
Beta_pos = zeros(ClassPops.k, dim); % β解
Beta_score = 1 ./ zeros(ClassPops.k, 1); % β解适应度
Delta_pos = zeros(ClassPops.k, dim); % δ解
Delta_score = 1 ./ zeros(ClassPops.k, 1); % δ解适应度
Fitness_list = zeros(ClassPops.k, Max_iter); % 适应度曲线
Pops.PosDim = WolfPops.PosDim; % 子种群
Pops.lb = WolfPops.lb;
Pops.ub = WolfPops.ub;
% 迭代求解
tic
fprintf('>>MP-GWO 优化中 00.00%%')
for iter = 1 : Max_iter
% ① 更新参数a
a = 2 - iter * 2 / Max_iter; % 线性递减
%a = 2 * cos((iter/Max_iter)*pi/2); % 非线性递减
% ② 聚类
if iter > 1 %初始化时聚完一次了，节省计算量
ClassPops = PopsCluster(WolfPops, UAV);
end
for k = 1 : ClassPops.k
Positions = ClassPops.Pos{k};
% ③ 寻找 α、β、δ 狼
for i = 1 : cSearchAgents
% 读取目标函数
fitness = ClassPops.Fitness(k, i);
% 更新 Alpha、Beta 和 Delta 解
if fitness <= Alpha_score(k) % 适应能力最强（因为性能指标越小越好，因此为小于号）
Alpha_score(k) = fitness;
Alpha_pos(k, :) = Positions(i, :);
end
if fitness > Alpha_score(k) && fitness <= Beta_score(k)
Beta_score(k) = fitness;
Beta_pos(k, :) = Positions(i, :);
end
if fitness > Alpha_score(k) && fitness > Beta_score(k) && fitness <= Delta_score(k)
Delta_score(k) = fitness;
Delta_pos(k, :) = Positions(i, :);
end
end
% ④ 更新位置（朝着前三只狼位置前进）
for i = 1 : cSearchAgents
for j = 1 : dim
r1 = rand();
r2 = rand();
A1 = 2*a*r1 - a;
C1 = 2*r2;
D_alpha = abs(C1*Alpha_pos(k, j) - Positions(i, j));
X1 = Alpha_pos(k, j) - A1*D_alpha;
r1 = rand();
r2 = rand();
A2 = 2*a*r1 - a;
C2 = 2*r2;
D_beta = abs(C2*Beta_pos(k, j) - Positions(i, j));
X2 = Beta_pos(k, j) - A2*D_beta;
r1 = rand();
r2 = rand();
A3 = 2*a*r1 - a;
C3 = 2*r2;
D_delta = abs(C3*Delta_pos(k, j) - Positions(i, j));
X3 = Delta_pos(k, j) - A3*D_delta;
% 静态更新
Positions(i, j) = (X1 + X2 + X3) / 3;
% 动态更新
% q = g * a; %阈值
% if abs(Alpha_score(k)-Delta_score(k)) > q
% Sum_score = Alpha_score(k) + Beta_score(k) + Delta_score(k);
% Positions(i, j) = (Alpha_score(k)*X1 + Beta_score(k)*X2 + Delta_score(k)*X3) / Sum_score;
% else
% Positions(i, j) = (X1 + X2 + X3) / 3;
% end
end
end
% ⑤ 调整不符合要求的状态变量
Pops.Pos = Positions;
ProbPoints = ClassPops.ProbPoints{k};
[Pops, ~] = BoundAdjust(Pops, ProbPoints, UAV);
% ⑥ 存储适应度
Fitness_list(k, iter) = Alpha_score(k);
% ⑦ 合并种群
WolfPops.Pos(cSearchAgents*(k-1)+1:cSearchAgents*k, :) = Pops.Pos;
end
if iter/Max_iter*100 < 10
fprintf('\b\b\b\b\b%.2f%%', iter/Max_iter*100)
else
fprintf('\b\b\b\b\b\b%.2f%%', iter/Max_iter*100)
end
end
fprintf('\n\n>>计算完成！\n\n')
toc
% 寻找 α β δ 位置
n = 3;
A = ClassPops.Fitness;
t = findmin(A, n);
index = cSearchAgents * (t(:, 1) - 1) + t(:, 2);
real_Alpha_no = index(1);
real_Beta_no = index(2);
real_Delta_no = index(3);
Alpha_Data = ClassPops.Data{t(1, 1)}{t(1, 2)} ;
% 输出值
solution.method = 'MP-GWO'; % 算法
% solution.ClassPops = ClassPops; % 分类信息
solution.WolfPops = WolfPops; % 所有解种群信息
solution.Tracks = Pops2Tracks(WolfPops, UAV); % 所有解航迹信息
solution.Fitness_list = mean(Fitness_list, 1); % 所有α解的平均适应度曲线
solution.Alpha_Data = Alpha_Data; % 真 · α 的威胁信息
solution.Alpha_no = real_Alpha_no; % 真 · α 的位置
solution.Beta_no = real_Beta_no; % 真 · β 的位置
solution.Delta_no = real_Delta_no; % 真 · δ 的位置
end
%% 寻找A矩阵中最小n个数的位置
function t = findmin(A, n)
t = sort(A(:));
[x, y] = find(A <= t(n), n);
t = [x, y]; % 前n个最小项在矩阵A中的位置[行,列]
B = zeros(n, 1);
for i = 1 : n
B(i) = A(t(i, 1), t(i, 2));
end
[~, index] = sort(B);
t = t(index, :); % 前n个从小到大排序的位置
end
%% 对种群进行聚类
function [ClassPops] = PopsCluster(WolfPops, UAV)
SearchAgents = size(WolfPops.Pos, 1); % 智能体个数
Dim = WolfPops.PosDim; % 智能体维度
Tracks = Pops2Tracks(WolfPops, UAV); % 智能体转换成航迹信息
% 计算适应度
o_Fitness = zeros(SearchAgents, 1); % 60*1
o_subF = []; % 5*60
o_ProbPoints = cell(SearchAgents, 1); % 60*1
o_Data = cell(SearchAgents, 1); % 60*1
% parfor 并行计算适应度
for i = 1:SearchAgents
[fitness, subF, Data] = ObjFun(Tracks{i}, UAV);
o_ProbPoints{i} = Data.ProbPoint; %cell-cell
o_Data(i) = {Data}; %cell-struct
o_Fitness(i) = fitness; %vector-var
o_subF = [o_subF, subF];
end
% 分类
k = size(subF, 1); % 分 k 类（由objfun决定）
cSearchAgents = floor(SearchAgents / k); % 并行智能体个数
cFitness = zeros(k, cSearchAgents); % 保存每类的适应度
cPositions = cell(k, 1); % 保存每类的位置信息
cTracks = cell(k, 1); % 保存每类的航迹信息
cProbPoints = cell(k, 1); % 保存每类的有问题航迹点
cData = cell(k, 1); % 存储每类的检测报告
% 排序
[~, Index] = sort(o_subF, 2, "ascend") ; % 沿维度2升序排序，返回新矩阵和序号
% fitness越小越好
% 聚类
for i = 1:k
Positions = zeros(cSearchAgents, Dim);
batchTrack = cell(cSearchAgents, 1);
batchProbPoints = cell(cSearchAgents, 1);
batchData = cell(cSearchAgents, 1);
for j = 1:cSearchAgents
idx = Index(i, j);
cFitness(i, j) = o_Fitness(idx); %mat-vector
Positions(j, :) = WolfPops.Pos(idx, :); %mat-mat
batchTrack{j} = Tracks{idx}; %cell-cell
batchProbPoints{j} = o_ProbPoints{idx}; %cell-cell
batchData{j} = o_Data{idx}; %cell-cell
end
cPositions(i) = {Positions}; %cell-mat
cTracks{i} = batchTrack; %cell-cell
cProbPoints{i} = batchProbPoints; %cell-cell
cData{i} = batchData; %cell-cell
end
% 输出
ClassPops.Pos = cPositions;
ClassPops.Tracks = cTracks;
ClassPops.ProbPoints = cProbPoints;
ClassPops.Data = cData;
ClassPops.Fitness = cFitness;
ClassPops.SearchAgents = cSearchAgents;
ClassPops.k = k;
end
🎉第三部分——参考文献
文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。(文章内容仅供参考，具体效果以运行结果为准)
[1]周瑞,黄长强,魏政磊,赵克新.MP-GWO 算法在多 UCAV 协同航迹规划
中的应用[J].空军工程大学学报(自然科学版),2017,18(05):24-29.
[2]胡中华,赵敏,姚敏,李可现,吴蕊.一种改进蚂蚁算法的无人机多目标三
维航迹规划[J].沈阳工业大学学报,2011,33(05):570-575.
[3]柳长安,王晓鹏,刘春阳,吴华.基于改进灰狼优化算法的无人机三维航迹
规 划 [J]. 华 中 科 技 大 学 学 报 ( 自 然 科 学 版 ),2017,45(10):38-
42.DOI:10.13245/j.hust.171007.
🌈第四部分——本文完整资源下载
资料获取，更多粉丝福利，MATLAB|Simulink|Python|数据|文档等完整资源获取
本文完整资源下载
