简介一套基于MATLAB的多无人机群飞行仿真代码包面向无人机编队控制研究者与MATLAB学习者可用于多机协同飞行动力学建模、Simulink仿真与三维可视化演示。压缩包共38个文件大小仅417KB涵盖7个m程序脚本、2个mdl模型、2个fig界面、wrl/obj三维模型、c/h源文件、dll动态库以及AIAA论文和使用说明文档等内容完整且便于查阅。已有126人浏览学习适合快速上手复现仿真过程。资源内置主函数main.m与多个辅助函数可直接在MATLAB 2020b中运行得到多无人机编队飞行结果搭配wrl三维场景与fig图形界面可直观观察集群运动轨迹与姿态变化同时使用说明文档逐项讲解文件作用对理解无人机集群协同控制及MATLAB/Simulink建模仿真具有较高参考价值。1. 群飞行仿真在MATLAB里到底仿真什么多无人机群飞行仿真比单机仿真难在“多机同时在场”。单机只需要管一条航迹群飞却要同时维护队形、机间距离和控制律三件事互相干扰。MATLAB适合做这件事因为状态可以组织成矩阵一次矩阵运算更新所有飞机而且自带绘图命令能实时看轨迹。标题里的“使用说明文档.rar”也代表了这类项目的常态一套MATLAB脚本加一份说明文档脚本跑通仿真文档写清参数怎么改。这篇文章按这个思路展开先讲运动学模型再写仿真主循环然后调编队PID最后加避碰并给出验证指标。新手能照做熟手可以直接看第4章的调参表和最后一章的验证代码。2. 多无人机群飞行的运动学模型与编队关系这一章先把理论基础打实用什么模型表示无人机、怎么在MATLAB里组织多机状态、以及“队形”如何用数学关系描述。如果这两步不做好后面调PID时每改一次参数都要跑一遍全仿真效率极低。2.1 用三自由度平面模型描述定高飞行群飞行仿真里很少直接上六自由度全量模型。六自由度要处理姿态角、角速度、控制力矩参数一多队形控制的逻辑就被淹没了。常见做法是先做定高飞行假设用三自由度平面模型近似每架无人机水平位置(x,y)、前向速度v、偏航角ψ控制量是速度指令和转向角速率。运动学更新公式是x(k1) x(k) v·cos(ψ)·dty(k1) y(k) v·sin(ψ)·dtψ(k1) ψ(k) ω·dt这套模型适合验证编队算法、PID调参和路径跟随但如果要仿真无人机的滚转、俯仰或者风速突变下的姿态变化就需要切回六自由度模型。做群飞课题时我一般先用三自由度把编队逻辑跑通再针对单机替换成六自由度模型。在多机仿真里每架无人机的状态可以放在矩阵的一行列变量含义单位1x水平横坐标m2y水平纵坐标m3vx水平速度x分量m/s4vy水平速度y分量m/s5psi偏航角rad6omega偏航角速率rad/s对应到MATLAB初始化时直接构造一个N×6的零矩阵再按列赋值。如果还要保存每个时刻的飞行数据就再加一个时间维变成 (nt, N, 6) 的三维数组后文的traj就是这样组织的。2.2 领航者-跟随者用偏移矩阵定义队形多机群飞最简单的队形表达是领航者-跟随者。指定其中一架为领航者其余每架飞机相对领航者的期望位置是一个固定的二维偏移。设领航者当前位置是[xL, yL]第i架跟随者的偏移是[dx_i, dy_i]那么该机的期望位置就是desired_i [xL dx_i, yL dy_i]这个方法的好处是队形切换非常直接想从V字变成楔形只要替换偏移矩阵。下面的代码用5架飞机生成一个简单的V字队形N 5; offset [0, 0; % 领航者作为跟随者偏移0 -2, 2; % 左翼 2, 2; % 右翼 -4, 4; 4, 4];offset第一行是全零表示领航者跟随自己。后面四行定义了左右两侧各两架飞机的相对位置单位是米。调整第2列可以改变编队前后间距调整第1列可以改变翼展宽度。offset矩阵的编码顺序很重要如果领航者不是第1行跟随控制程序里就要单独记录leader_index否则desired矩阵的行顺序会错位。真正的跟随控制里领航者自己可以沿航路点飞行跟随者用下一章的主循环去跟踪desired点即可。2.3 一致性控制让速度自动对齐位置偏移法在直线编队时表现很好但一旦领航者转弯跟随者会走内切或外切。一个更稳的做法是在期望位置之外再加一个速度一致性约束让每架飞机的期望速度与邻居飞机平均速度一致。用数学语言表达就是v_des_i v_leader K_align * (v_neighbor_avg - v_i)在实际MATLAB代码里这一项通常表现为速度环的前馈。计算期望速度时先把领航者速度前馈加上再叠加队形保持产生的修正量。这样就避免了一致性算法完全接管队形位置。要注意的是一致性算法对通信延迟敏感延迟过大会出现整队振荡这个坑在第5章里再展开。我一般会把一致性项放在PID期望速度里而不是直接改变水平位置这样队形保持和速度同步互不干扰。3. 在MATLAB中搭建群飞行仿真主循环理论模型定下来之后下一步是让多架飞机在MATLAB里动起来。一个可读性好的群飞仿真主循环应该包含三部分时间步进、控制律计算、状态更新。这里给出一个最小可运行的脚本重点看程序骨架和MATLAB的数据组织。3.1 最小可运行的群飞主循环下面是四架飞机直接朝各自目标点移动的版本。虽然还没加编队控制但循环结构和矩阵索引方式已经完整。T 20; % 仿真总时长秒 dt 0.05; % 仿真步长秒 nt round(T / dt); N 4; % 无人机数量 target [0, 10; -2, 8; 2, 8; 0, 6]; pos [0, 0; -1, 1; 1, 1; 0, -1]; vel zeros(N, 2); traj zeros(nt, N, 2); % 保存全部轨迹 traj(1,:,:) pos; for k 1:nt-1 for i 1:N err target(i,:) - pos(i,:); dist norm(err); if dist 0.1 continue; end v_des err / dist * 2.0; % 目标方向乘以巡航速度 vel(i,:) v_des; % 本版本简单设速度 pos(i,:) pos(i,:) vel(i,:) * dt; end traj(k1,:,:) pos; end逻辑说明err是目标点与当前位置的差dist是距离。如果距离小于0.1米就认为到达目标结束这一架飞机的位置更新避免在目标点附近来回抖。v_des用err/dist归一化方向再乘2m/s巡航速度。这个版本没有加速度环节所以速度变化是瞬时的适合先验证循环是否有问题。参数说明dt是离散化步长它同时决定积分精度和动画帧率。0.05秒意味着每秒20次更新对平面运动学模型足够。如果dt取得过大飞机单步位移超过0.1米时dist判断可能失效。巡航速度2.0是仿真设定值改成3.0后要注意转弯半径是否超过安全距离。写完这个脚本建议先运行一下并看pos和target是否在几分钟内收敛。如果轨迹发散优先检查矩阵维度的顺序很多新手会在这里把N和nt搞混。3.2 把领航者航路点和队形偏移接进循环上面的target是固定目标。群飞仿真里目标点要随时间变化。常见的做法是定义一条路径path用当前时刻k查表得到领航者位置然后加上第2章的offset得到所有飞机的期望位置。path zeros(nt, 2); path(:,2) linspace(0, 20, nt); % 领航者直线北飞 offset [0, 0; -2, 2; 2, 2; -4, 4; 4, 4]; N size(offset, 1); desired zeros(N, 2); for k 1:nt leader path(k,:); desired leader offset; end这段代码里path的每一行是领航者在第k个时刻的位置。desired是一个N×2矩阵每一行是第i架飞机在当前时刻应该处于的位置。注意offset矩阵的行数决定了无人机总数N所以增加飞机只需要在offset里加行。实际控制时不能把desired直接当成目标点因为第3.1节的位置更新没有加速度模型速度跳变更明显。在加入PID控制后desired只是位置环的参考值速度环会另外计算。3.3 用animatedline画实时轨迹调试群飞行仿真只看终点坐标不够必须看中间过程。MATLAB里最顺手的是animatedline它支持追加数据点配合drawnow每步刷新。绘图方式特点适用场景plot一次性画出全部轨迹仿真结束后出图animatedline边算边追加轨迹点调试实时轨迹quiver绘制速度向量箭头观察速度方向一致性scatter绘制离散散点显示当前时刻机群位置h animatedline(Color, r, LineWidth, 1.5); for k 2:nt addpoints(h, traj(k, 1, 1), traj(k, 1, 2)); drawnow limitrate % 降低刷新频率提高速度 end上面的代码只画了第1架飞机的轨迹。要画全部机群可以用hold on创建多个animatedline句柄每个句柄对应一架飞机。drawnow limitrate是MATLAB R2018a之后支持的刷新方式它每秒最多刷新20帧可以避免每步drawnow的卡顿。在飞机数量少的时候直接drawnow也行数量超过10架后建议把轨迹先存入traj再离线重绘。4. 编队控制算法与PID参数调节的现场手法主循环跑起来之后飞机的运动就是“瞬移”式的和真实无人机的惯性不符。要让它在朝队形点移动时平滑、不超调、不偏航需要引入PID控制。这里不讨论Simulink里的模块拖拽直接用脚本写离散PID这样参数可以批量修改。4.1 位置环PD与速度环P组成的串级结构群飞编队控制不推荐一个PID直接输出速度指令。更稳的结构是位置外环产生期望速度速度内环产生加速度或油门指令。位置环用PD调节速度环用P就够。离散形式是v_des Kp * e_pos Kd * (e_pos - e_pos_prev) / dtacc Kvp * (v_des - vel)其中e_pos是期望位置减当前位置。速度环误差e_vel v_des - vel然后乘以Kvp得到加速度。注意这里的acc再乘dt加到vel上避免出现速度跳变。Kp 1.2; Kd 0.4; Kvp 2.0; v_max 3.0; e_prev zeros(N, 2); for k 1:nt-1 % 计算每架飞机的期望位置 desired这里省略更新逻辑 e_pos desired - pos; v_des Kp * e_pos Kd * (e_pos - e_prev) / dt; v_des max(min(v_des, v_max), -v_max); e_vel v_des - vel; acc Kvp * e_vel; vel vel acc * dt; pos pos vel * dt; e_prev e_pos; traj(k1,:,:) pos; end注意这里desired需要在每个时刻更新代码里省略了更新逻辑。v_des进行了限幅避免Kp大时输出过快。e_prev必须在更新前保存否则微分项会变成零。acc是加速度速度环只有比例控制所以有稳态误差但在编队控制里位置环的积分效应会消除它。4.2 参数调节表怎么判断Kp、Kd、Kvp该往哪个方向拧手动调参时只看轨迹很难判断问题。下面表整理最常遇到的现象和调节方向现象原因调整方法队形离目标点远但稳定Kp太小增大Kp或加积分项飞机在目标点附近抖动位置环Kd或速度环Kvp过大降低Kvp或调大Kd增加阻尼转弯时外抛或内切Kd不够速度前馈缺失增大Kd并加领航者速度前馈响应很慢、跟不上领航者Kp/Kvp太低同时提高Kp和Kvp一个相对稳妥的初始参数范围是Kp在0.8~2.0Kd在0.2~0.8Kvp在1.0~3.0。如果发现速度曲线像锯齿先检查v_max是否限幅如果位置误差出现正弦形振荡把Kvp降到1.0以下。提示如果编队中某一架飞机总是最后一个到位单独调它的Kp就行不要动不动就动全局参数。4.3 用误差平方和判断仿真参数是否可用调参不能凭肉眼。我会在每次仿真后计算整个队列的队形误差平方和然后绘制随时间变化的曲线。曲线应该迅速下降并保持平稳。err_all desired - pos; err_sq sum(err_all.^2, 2); % 两列求平方和后按行求和 figure; plot((0:nt-2)*dt, err_sq(1:end-1)); xlabel(时间/s); ylabel(队形误差平方和/m^2);注意desired和pos都需要保存成二维数组这里为了简洁只显示了最终时刻的计算。真正使用时要在每个循环里把(e_pos.^2)按行求和存进err_sq向量。err_sq如果持续下降说明稳定如果在某个值附近等幅振荡说明Kp/Kvp比例过高先按表里的方向调低。5. 多机避碰、避障与鲁棒性设计编队控制可以让机群保持队形但它不负责回答“两架飞机太近怎么办”。在实际飞行中队形控制律和避碰逻辑必须同时运行否则一次GPS跳变就可能让两架飞机撞上。这一章加两层保护机间排斥场和固定障碍物避让并讨论通信延迟的影响。5.1 机间排斥场在PID输出上叠加排斥力最简单的机间防撞是在控制加速度上叠加一个随距离减小的排斥项。当两架飞机距离小于安全距离d_safe时产生排斥加速度距离越近排斥力越大。方向是两架飞机连线的反向。代码如下d_safe 1.5; rep_gain 1.0; for i 1:N for j i1:N delta pos(i,:) - pos(j,:); d norm(delta); if d d_safe d 0.01 dir delta / d; factor rep_gain * (d_safe / d - 1); acc(i,:) acc(i,:) factor * dir; acc(j,:) acc(j,:) - factor * dir; end end end逻辑是只有距离小于安全阈值的飞机对参与计算。factor在d等于d_safe时趋近于0在d远小于d_safe时迅速增大。注意要优先累加避碰加速度再叠加PID控制加速度顺序不同会影响最终行为。参数说明d_safe要大于飞机物理半径加上定位误差通常设为编队间距的40%~60%。rep_gain过大会把飞机弹开导致队形扭曲过小则防不住碰撞。一个起点值是队形间距的0.2倍。这个方案在静态环境中很好用但要注意排斥场会造成队形“呼吸”现象飞机靠近安全距离就被弹开随后PID又往回拉形成周期性抖动。遇到这种情况把rep_gain调小或对factor加一个死区。5.2 圆盘障碍物避让固定障碍物可以用圆盘模型描述。对每个障碍物计算无人机到圆心距离若小于碰撞半径就在加速度上叠加一个远离圆心的分量。这里不再给完整循环只给出核心运算obs_center [5, 5]; obs_radius 1.2; delta pos(i,:) - obs_center; d norm(delta); if d obs_radius margin acc(i,:) acc(i,:) rep_gain * delta / d * (obs_radius margin - d); end注意这个避障项需要与机间排斥场共用同一个acc行向量因此在主循环里要先清空acc再依次加PID、机间排斥、避障项。margin是安全余量可以根据无人机尺寸设0.2~0.5米。障碍物圆心的位置可以放在脚本开头的可配置区域方便换场景。5.3 通信延迟与丢包的仿真方式一致性控制很依赖邻居状态是否新鲜。仿真通信延迟的常见做法是维护一个状态缓冲队列控制律读取delay_steps之前的状态而不是当前状态。丢包则是按一定概率跳过状态更新。下表是三个可直接调整的仿真参数参数含义建议范围delay_steps通信延迟对应的仿真步数1~3drop_rate单次状态更新丢失概率0~0.1broadcast_interval状态广播间隔秒0.1~0.5在MATLAB里可以用一个三维缓冲数组保存最近几个时刻的所有飞机状态。控制律使用时计算delay max(1, k-delay_steps)然后从buffer(delay,:,:)读取数据。丢包则用rand判断本次更新是否被丢弃如果丢弃前一刻的状态就会被重复利用一次。延迟增加时位置环Kp要适当降低否则会出现整队振荡。丢包率超过10%时一致性算法的收敛性会明显下降此时可以让速度前馈更依赖自身历史值而不是邻居状态。6. 用误差指标和动画验证群飞行仿真效果仿真结束后需要回答三个问题队形是否保持、机间是否安全、参数改了以后效果如何。这一章给出几个可以直接复制到MATLAB里的验证代码块。6.1 队形保持误差RMSE与最大误差队形保持的量化指标是每架飞机实际位置与期望位置的偏差。计算每一时刻的全部飞机误差再求均方根和最大值err_hist zeros(nt, 1); for k 1:nt err_hist(k) sqrt(mean(sum((desired(k,:,:) - traj(k,:,:)).^2, 3))); end rmse sqrt(mean(err_hist.^2)); max_err max(err_hist); fprintf(RMSE%.3f m, max_err%.3f m\n, rmse, max_err);这里的desired和traj都是三维数组时间×机号×坐标。sum(...,3)是沿着坐标维求和得到每一时刻每架飞机的误差平方再对机数求mean开方后得到该时刻的均方根误差。最终rmse是全程平均。如果max_err总是小于安全距离说明即使极端的飞行时刻也没有碰撞风险。6.2 最小机间距离统计群飞安全性的另一个硬指标是全程任意两架飞机之间的最小距离。用二维扫描即可min_d inf; for k 1:nt for i 1:N for j i1:N d norm(squeeze(traj(k,i,:) - traj(k,j,:))); if d min_d min_d d; end end end end disp([min distance: , num2str(min_d), m]);注意squeeze去掉单例维确保差的结果是二维向量。如果min_d小于你的安全阈值就要去回看k索引附近的时间段通常问题发生在领航者转弯或编队切换的瞬间。6.3 把参数集中在结构体里方便和使用说明文档对应一个群飞仿真脚本跑通后另一个常见需求是让其他人也能接手改参数。把所有仿真参数放进一个struct命名字段和说明文档保持一致比散落的全局变量好维护得多param.dt 0.05; param.N 6; param.Kp 1.2; param.Kd 0.4; param.Kvp 2.0; param.d_safe 1.5; param.v_max 3.0;说明文档可以按字段逐行解释这些值的作用。使用说明文档通常不需要描述算法本身只要写“调大param.d_safe会增大机间排斥距离队形会显得更松散”这类直接因果关系。这样下来脚本和文档就能形成完整交付。如果机群规模扩大到20架以上要先把traj保存到工作区再离线绘制和导出动画实时drawnow很容易成为性能瓶颈。这能帮你快速区分是算法问题还是渲染问题。本文还有配套的精品资源点击获取