简介一套面向固定翼无人机研发的 MATLAB 代码包适用于航空工程、自动化控制等专业的学生、研究者与工程师用于解决无人机动力学建模、飞行控制律设计及仿真验证等问题。代码覆盖固定翼无人机三维动力学模型搭建、PID 控制器与状态反馈控制器设计、Simulink 仿真环境构建、路径规划与避障算法验证等核心模块可直接在 MATLAB 中运行或二次开发。压缩包共 46 个文件以 39 个 .m 源代码文件为主另含 2 张用于展示模型与结果的 png 图片、1 份 rtf 格式的环境配置与运行说明、1 个自动备份文件及少量辅助文件整体仅 123KB下载与解压十分便捷。已有 231 人学习该代码包其结构清晰、注释得当既适合快速复现仿真结果也便于按需修改控制参数是学习固定翼无人机飞控与建模仿真的实用参考资料。1. 固定翼无人机 matlab 代码.rar先分清你要的是模型、控制还是整套仿真解压一个固定翼无人机 matlab 代码.rar 之前先想清楚一个问题你是来学六自由度建模还是来调 PID还是只要一个能出图的动画演示这个选择直接决定你先打开哪个文件。我见过太多人一上来就盯主脚本结果被几百行的 S-function 和回调函数绕晕最后连模型跑没跑起来都说不清。固定翼无人机 matlab 代码的常见构成其实很固定一个气动参数初始化脚本、一个六自由度状态导数函数、一个配平脚本、一个控制律串级 PID 居多、一个仿真主循环外加画图脚本。其中真正决定仿真能不能收敛的不是控制律而是那一组气动导数和单位。这篇文章按我平时搭一套可复现固定翼仿真环境的顺序讲先把 12 维状态的微分方程写清楚再把配平和串级 PID 用 matlab 代码落下来最后给一组能直接跑的初始参数和排错清单。适合做毕设、入门飞控算法、要给固定翼调参但缺个沙盘的工程师。读完你能让一套模型在 MATLAB 里从零跑出平飞和阶跃响应而不是对着别人代码改注释。2. 六自由度固定翼模型怎么落成 matlab 代码状态、气动力与导数的处理顺序2.1 先定状态向量12 维排列决定后面所有索引对不对固定翼六自由度模型的状态量常规取法是从惯性位置、机体速度、欧拉角到角速度排成 12 维列向量。我一般用这个顺序% 状态向量 x [pn, pe, pd, u, v, w, phi, theta, psi, p, q, r] % pn/pe/pd: 北东地坐标pd 向下为正单位 m % u/v/w: 机体轴速度分量m/s % phi/theta/psi: 滚转/俯仰/偏航角rad % p/q/r: 机体轴角速度rad/s这个排列的好处是后面按索引切片时不用来回改配平脚本、控制律、画图脚本都共用同一个列顺序。改成别的排列比如把欧拉角放前面不是不行但所有下游脚本都要跟着改鸡毛蒜皮的错最容易出在这。固定翼无人机 matlab 代码里另一个常见问题是把高度直接当 pd 用。注意 NED 坐标系下 pd 向下为正真高度 h -pd。控制律里写高度误差时是 h_c - h不是 pd_c - pd这个符号错一次你的高度环会以爬升 → 掉高 → 爬升的振荡结束。2.2 机体轴运动方程力、角速度、姿态传播三个子块有了状态向量状态导数的计算分三块。第一块是力方程把重力、推力、气动力都投影到机体轴第二块是角速度方程用惯量矩阵和力矩第三块是用方向余弦矩阵把角速度转成欧拉角导数。核心代码function xdot fixedwing6dof(x, u, p) % 状态解包 phi x(7); theta x(8); psi x(9); uu x(4); vv x(5); ww x(6); pp x(10); qq x(11); rr x(12); Va sqrt(uu^2 vv^2 ww^2); alpha atan2(ww, uu); % 迎角 beta atan2(vv, Va); % 侧滑角 qbar 0.5 * p.rho * Va^2; % 动压 % ---- 气动力系数纵向为主横航向侧力简化 ---- CL p.CL0 p.CLalpha * alpha p.CLde * u(1) ... p.CLq * (p.cbar / (2*Va)) * qq; CD p.CD0 p.CDalpha * alpha ... p.CDq * (p.cbar / (2*Va)) * qq; CY p.CYbeta * beta p.CYda * u(2) p.CYdr * u(3); % 升力垂直速度方向阻力沿速度反方向先转到机体轴 F_aero_b qbar * p.S * [ -CD * cos(alpha) CL * sin(alpha); CY; -CD * sin(alpha) - CL * cos(alpha) ]; % 重力在机体轴的分量NED 地轴系 [0;0;mg] 旋转到机体轴 F_grav_b p.m * p.g * [-sin(theta); sin(phi)*cos(theta); cos(phi)*cos(theta)]; % 推力默认沿机体 x 轴可加安装角偏移 F_thrust_b [p.Tmax * u(4); 0; 0]; % 线加速度 f_b F_aero_b F_grav_b F_thrust_b; udot f_b(1)/p.m - qq*ww rr*vv; vdot f_b(2)/p.m - rr*uu pp*ww; wdot f_b(3)/p.m - pp*vv qq*uu; % ---- 力矩与角加速度对称飞机惯量近似处理 ---- Cl p.Clbeta * beta p.Clp * (p.b/(2*Va))*pp ... p.Clda * u(2) p.Cldr * u(3); Cm p.Cm0 p.Cmalpha * alpha p.Cmq * (p.cbar/(2*Va))*qq ... p.Cmde * u(1); Cn p.Cnbeta * beta p.Cnr * (p.b/(2*Va))*rr ... p.Cnda * u(2) p.Cndr * u(3); pdot (p.G1*pp*qq - p.G2*qq*rr p.G3*Cl*qbar*p.S*p.b ... p.G4*Cn*qbar*p.S*p.b); qdot (p.G5*pp*rr - p.G6*(pp^2 - rr^2) Cm*qbar*p.S*p.cbar)/p.Jy; rdot (p.G7*pp*qq - p.G1*qq*rr p.G4*Cl*qbar*p.S*p.b ... p.G8*Cn*qbar*p.S*p.b); % ---- 姿态与位置传播 ---- pndot cos(theta)*cos(psi)*uu (sin(phi)*sin(theta)*cos(psi)-cos(phi)*sin(psi))*vv ... (cos(phi)*sin(theta)*cos(psi)sin(phi)*sin(psi))*ww; pedot cos(theta)*sin(psi)*uu (sin(phi)*sin(theta)*sin(psi)cos(phi)*cos(psi))*vv ... (cos(phi)*sin(theta)*sin(psi)-sin(phi)*cos(psi))*ww; pddot -sin(theta)*uu sin(phi)*cos(theta)*vv cos(phi)*cos(theta)*ww; phidot pp qq*tan(theta)*sin(phi) rr*tan(theta)*cos(phi); thetadot qq*cos(phi) - rr*sin(phi); psidot (qq*sin(phi) rr*cos(phi))/cos(theta); xdot [pndot; pedot; pddot; udot; vdot; wdot; ... phidot; thetadot; psidot; pdot; qdot; rdot]; end代码逻辑说明气动力先算到气流系再通过迎角旋转到机体轴。这是多数固定翼仿真代码采用的简化解法直接按机身轴分解 CL/CD 会丢前后力耦合速度越高的飞机误差越明显。重力投影是新手最容易写反的地方。NED 系下重力矢量是 [0;0;mg]转到机体轴后 x 分量是 -mg·sinθ平飞小迎角下它和推力的平衡决定配平速度。横航向力矩用了 G1~G8 这组惯量组合参数它们由 Ixx、Iyy、Izz、Ixz 计算得到。对称面上 Ixz0 时 G2、G5 等项会自行退化为简化公式但保留通用形式方便以后改模型。参数说明p.rho 是空气密度一般初始化成 1.225 kg/m³高原场景按高度查表或直接乘密度比例。p.Tmax 是最大推力u(4) 是油门指令 0~1推力按线性近似。升降舵、副翼、方向舵指令 u(1)、u(2)、u(3) 单位是弧度控制律输出给舵面时要注意限幅通常 ±0.5 rad 内。2.3 气动导数表没有风洞数据时的一组可用起点拿不到真实气动数据时用下面这组参数就能让模型在 20 m/s 附近稳定平飞。单位符号一定要盯住CLa 的单位是 1/rad舵面效率导数同理。参数含义典型值1.5m 翼展电动固定翼单位S机翼面积0.36m²b / cbar翼展 / 平均气动弦长1.5 / 0.24mm质量1.8kgIxx / Iyy / Izz三轴转动惯量0.045 / 0.075 / 0.12kg·m²CL0 / CLa / CLde零升升力 / 升力线斜率 / 升降舵效率0.28 / 5.0 / 0.351/radCD0 / CDa零升阻力 / 迎角阻力0.03 / 0.121/radCm0 / Cma / Cmq / Cmde俯仰力矩各项-0.02 / -1.2 / -8.0 / -0.91/rad 等Clbeta / Clp / Clda滚转气动导数-0.12 / -0.9 / 0.151/radCnb / Cnr / Cnda偏航气动导数0.12 / -0.08 / -0.061/rad使用提示Cma 必须为负这是俯仰静稳定的前提。算出来为正说明重心太靠后仿真会直接俯仰发散。Clp、Cnr 都是阻尼导数为负表示滚转和偏航运动自身消耗能量。如果阻尼导数设成正横航向会以螺旋发散收场。这组数据来自我常用的教学级小固定翼模板。拿真实无人机时把 S、b、cbar、质量惯量换成实测值气动导数仍可先用这组做量级验证再按风洞或 CFD 结果替换。3. 固定翼 matlab 仿真闭环配平、主循环与串级 PID 的实现顺序3.1 配平为什么固定翼代码跑起来前必须先找平飞点固定翼模型不像四旋翼给了油门就能悬停。它在任意给定速度下需要一组迎角 舵面 油门让线加速度和角加速度同时为零这就是配平点。没配平就仿真飞机会立刻抬头或俯冲你根本分不清是控制律没调好还是初始状态不对。配平问题的数学形式是求解 f(x_trim, u_trim) 0 的前若干行。纵向配平最常见固定 θ、油门和升降舵求解 α。用 fsolve 可以把残差函数写得很直观function F trim_residual(y, p) % y [alpha; de; dt] alpha y(1); de y(2); dt y(3); theta alpha; % 平飞时航迹角 gamma0thetaalpha Va 20; % 目标空速 uu Va * cos(alpha); ww Va * sin(alpha); x [0; 0; -100; uu; 0; ww; 0; theta; 0; 0; 0; 0]; u [de; 0; 0; dt]; xdot fixedwing6dof(x, u, p); % 只保留纵向力、俯仰力矩和油门相关的残差 F [xdot(4); xdot(6); xdot(11)]; % udot, wdot, qdot end逻辑说明平飞且无侧滑时φ、β、侧向速度全取零横航向方程自动满足只校纵向。残差选 udot、wdot、qdot 三项对应三个未知量 α、de、dt。这样方程组闭合fsolve 才有唯一解。目标空速 Va 可调。想得到不同巡航速度的配平点改成输入参数循环调用即可这其实就是 matlab 优化工具箱在固定翼无人机 matlab 代码里最常见的用法不要自己去写牛顿迭代。调用配平代码p.Tmax 8; % 最大推力 N y0 [5*pi/180; -0.05; 0.4]; % alpha, de, dt 的初值 opt optimoptions(fsolve, Display, off, ... Algorithm, trust-region, ... MaxIterations, 200); y_trim fsolve((y) trim_residual(y, p), y0, opt); fprintf(alpha %.2f deg, de %.2f deg, dt %.2f\n, ... y_trim(1)*180/pi, y_trim(2)*180/pi, y_trim(3));参数说明fsolve 的初值要落在合理物理区间α 取 2°~8°de 取 -0.1~-0.02 rad升降舵常规是负偏配平dt 取 0.3~0.7。初值离谱时 trust-region 算法容易给出负油门或大舵偏的解那个解数学上成立但物理上不可飞。配平输出直接作为后面初始化脚本里状态和指令的初值。控制律的积分项也从配平舵面开始积能显著缩短初始振荡。3.2 主循环用定步长 RK4 而不是 ode45 跑整段仿真很多现成固定翼无人机 matlab 代码用 ode45 集成。我的做法是定步长 RK4原因有三控制律和状态更新必须同频ode45 变步长会把控制量离散时间戳搞乱中途要随时插参数扰动或注入风场变步长下反而难处理最后画图和判稳都依赖固定时间网格。主循环核心dt 0.01; t_end 30; N round(t_end / dt); x x0; t 0; hist zeros(13, N); for k 1:N % 控制律接收当前状态输出舵面和油门 u controller(x, p, cmd); % RK4 积分 k1 fixedwing6dof(x, u, p); k2 fixedwing6dof(x 0.5*dt*k1, u, p); k3 fixedwing6dof(x 0.5*dt*k2, u, p); k4 fixedwing6dof(x dt*k3, u, p); x x (dt/6) * (k1 2*k2 2*k3 k4); hist(:, k) [t; x]; t t dt; end逻辑与参数说明dt 取 0.01 s100 Hz是平衡点。固定翼短周期模态一般在 2~4 rad/s100 Hz 积分足以覆盖一个周期 150 个采样点取 0.001 s 跑 30 s 仿真要多等很久不划算。RK4 每步调用 4 次导数函数控制律只用当前状态计算一次这在教学仿真里是正确的离散化顺序。注意不要在 k2、k3 里重新计算控制量那会导致控制频率隐式变成 400 Hz。控制律输出限幅放在 controller 内部最后一行限幅后再进积分器。限幅放在里面还是外面对积分饱和的影响完全不同。3.3 串级 PID 怎么接进固定翼代码滚转、俯仰、高度三层环固定翼控制律的通用结构是最外层是高度/航向制导中间层是姿态角环最内层是角速度阻尼环。外环输出给内环做期望值内环输出舵面。以纵向为例典型实现function u controller(x, p, cmd) phi x(7); theta x(8); pp x(10); qq x(11); pd x(3); uu x(4); ww x(6); h -pd; % NED 转高度 Va sqrt(uu^2 ww^2); % 高度环输出期望俯仰角限幅 ±25° theta_c cmd.h_c p.kp_h * (cmd.h_c - h); theta_c max(min(theta_c, 25*pi/180), -25*pi/180); % 俯仰角环输出期望俯仰角速度 q_c p.kp_theta * (theta_c - theta); % 角速度环 配平舵面前馈 de p.kp_q * (q_c - qq) p.de_trim; % 滚转通道同理保持机翼水平 p_c p.kp_phi * (-phi); % 期望 phi0 da p.kp_p * (p_c - pp) p.da_trim; u [de; da; 0; cmd.dt]; end参数说明与整定顺序整定顺序必须先内后外先只保留 q 环调 Kp_q让俯仰角速度响应在 0.5 s 内收敛再加 θ 环调 Kp_theta最后加高度环调 Kp_h。外环增益单位分别是不同量纲不能直接和 Kp_q 比大小。Kp_q 从 0.1 起步每次翻倍直到角速度响应开始出现 1~2 次过冲再往回退 30%Kp_theta 从 0.5 起步出现频率约 1 Hz 的振荡就减半Kp_h 从 0.3 起步主要在 0.05~0.5 之间。高度环输出限幅是必须的。不做限幅高度误差一大期望俯仰角直接顶到 90°模型会进入深失速。限幅值取 ±20°~±25° 对大多数小型固定翼安全。前馈项 de_trim 是配平脚本算出来的舵面偏角。没有前馈纯靠积分项顶配平舵阶跃响应的初始下沉量会大很多。4. 固定翼代码跑起来的初始化资产与 4 个典型排错点4.1 初始化脚本把所有参数集中在一处改参不改逻辑固定翼无人机 matlab 代码最容易失控的地方是参数散落各个脚本。我一般用一个 init_fixedwing.m 把所有气动、惯量、控制增益和初始状态集中定义用结构体 p 传参这样配平、仿真、画图三套脚本共用同一份参数改一次全链路生效。% init_fixedwing.m p.g 9.80665; p.rho 1.225; p.S 0.36; p.b 1.5; p.cbar 0.24; p.m 1.8; p.Ixx 0.045; p.Iyy 0.075; p.Izz 0.12; p.Ixz 0; % 惯量组合参数 p.G1 p.Ixz*(p.Ixx - p.Iyy p.Izz) / p.Ixx; % 按标准公式计算 % ... 气动导数从第 2.3 节参数表逐行写入 ... p.Tmax 8; % 先配平得到 trim 状态和舵面再灌进初始状态 [alpha, de_trim, dt_trim] compute_trim(p, 20); p.de_trim de_trim; x0 [0; 0; -100; 20*cos(alpha); 0; 20*sin(alpha); ... 0; alpha; 0; 0; 0; 0]; cmd.h_c 100; % 期望高度这个脚本本身只定义参数和初值不执行积分。仿真脚本运行时第一行直接 init_fixedwing保证 restart 后状态一致不会出现上一次运行残留变量污染本次仿真的问题。需要批量跑不同速度的配平点时把 20 改成参数传入 compute_trim 即可p 结构体里多存一份目标空速字段。4.2 仿真跑飞前的检查顺序先看气动导数符号再看单位最后看环路模型发散时的检查顺序很重要乱试参数浪费时间。我按这个顺序排错命中率很高现象首要怀疑点检查手段一松开舵就俯冲或拉起Cma 符号、配平初值打印配平 α、de确认 de 为负小量滚转直接翻过去Clda 符号、da 指令方向断开控制器给 da0.1观察滚转角是否朝正方向滚侧滑越摆越大无法收敛Cnb、Cnr 符号或数值量级检查 Cnb 是否为正、Cnr 是否为负高度环高频振荡频率接近控制频率Kp_h 太大或角度环限幅饱和把 Kp_h 减半观察振荡频率是否跟着变只有油门变化飞机状态几乎不动推力单位或 Tmax 量级错算配平所需推力 0.5·ρ·V²·S·CD0对比 Tmax单位问题排查提示角度全部用弧度。气动导数表里迎角单位是 rad 的数据一旦程序里传了角度值升力线斜率会被放大 57 倍结果是任何控制器都拉不住。动压 qbar 结果量级可用手算验证20 m/s 时 qbar 0.5 × 1.225 × 400 ≈ 245 PaS0.36 时总气动力 ≈ 88 N除以质量 1.8 kg 是 49 m/s²和重力一个量级符合小型固定翼的预期。算出来差几个数量级就查单位。4.3 开头 0.5 s 的瞬态怎么消除把配平状态和配平舵面一起作为初值很多固定翼无人机 matlab 代码在仿真开始时会先来一轮大振荡常见原因是初始迎角和配平迎角不一致、或者升降舵从 0 开始而配平舵面是 -3°。消除办法是仿真代码里强制指定初值来源于配平% 仿真开始前做一致性检查 assert(abs(x0(8) - alpha_trim) 1e-6, ... 初始俯仰角未使用配平值请重新运行配平脚本); assert(abs(u_trim(1) - p.de_trim) 1e-6, ... 升降舵初值未对齐配平舵面);两个断言能把这类问题一次性堵死仿真从第一帧就是平飞状态。若你拿到的代码没有配平流程最粗暴但有效的替代方案是跑一段 5 s 无控制开环把第 5 s 末的状态当作新的 x0再叠加上控制器。这样虽然没解析配平精确但对初学验证控制器逻辑已经够用。5. 验证固定翼模型代码可信度的三个手段5.1 把量纲和物理常量单测一遍仿真代码全部写完先别急着看曲线。把 p 结构体里的每个量代入动压、重力、推力公式手算一个基准工况对比脚本输出。我用一个极短的对拍代码做这件事Va 20; alpha 5*pi/180; qbar 0.5*p.rho*Va^2; % 245 Pa L qbar*p.S*(p.CL0 p.CLalpha*alpha p.CLde*p.de_trim); assert(abs(L - p.m*p.g) 1.5, 升力与重力失衡);这行断言如果通过说明平飞点在物理上成立。多数模型代码的问题不是算法错而是这个基本守恒都没检查。升力与重力差超过 1.5 N约 8% 重力时先查 CL0 或 CLa 的量级再查动压是否由于空速单位错了而偏大偏小。5.2 开环阶跃响应验证模态合理性控制器关掉给升降舵一个 -0.05 rad 的阶跃观察高度和俯仰角速度响应。合格模型的短周期响应应表现为 1~2 s 内快速收敛的俯仰角速度振荡长周期表现为 10 s 量级的起伏爬升而不是直接发散。若俯仰角速度瞬间冲顶优先检查 Cmq 是否漏了这个阻尼导数缺失会让短周期模态呈等幅振荡甚至发散。5.3 蒙特卡洛扫气动参数确认控制裕度给气动导数加 ±10% 随机扰动重复跑 100 次统计高度保持误差和是否掉高超过 5 m。这是验证串级 PID 鲁棒性最便宜的做法matlab 里用 parfor 可以并行加速parfor i 1:100 p_i p; p_i.CLa p.CLa * (1 0.1*randn); p_i.Cmde p.Cmde * (1 0.1*randn); p_i.CD0 p.CD0 * (1 0.1*randn); res run_sim(p_i, cmd, dt, t_end); max_h_err(i) max(abs(res.h - cmd.h_c)); lost_height(i) max(0, cmd.h_c - min(res.h)); end fprintf(95%% 高度误差 %.2f m\n, quantile(max_h_err, 0.95));这个统计量比单次仿真曲线更能说明控制器有工程价值。若 95 分位高度误差小于 2 m、掉高小于 1 m这套模型和增益可以直接拿去验证路径跟踪逻辑若扰动下频繁掉高就回头把内环 Kp_q 加大内环硬了外环才敢加压。等内环裕度够了再把 CLa 的扰动范围从 ±10% 放宽到 ±20% 复测一遍控制器如果仍能压住高度误差这套固定翼无人机 matlab 代码才算真正到了可以信任的程度。本文还有配套的精品资源点击获取