简介这套基于MATLAB的小型飞机纵向动力学非线性动态反演控制器程序包支持MATLAB 2014/2019a/2021a等版本面向飞行控制与非线性控制方向的本科、硕士教研场景可用于学习动态逆方法在纵向通道控制器设计中的建模、控制器推导、仿真与结果分析。压缩包共4个文件总大小约54KB包含1个MATLAB脚本源码、2个txt说明文档含使用说明与开源许可以及1张运行结果截图文件结构简洁便于快速对照运行。目前已有118人学习下载资源小巧下载后即可解压使用适合作为课程设计、毕业设计或科研入门的参考实现。通过运行压缩包内的MATLAB脚本读者可直观了解非线性动态反演控制器的构建流程与仿真响应曲线并结合说明文档完成参数调整、结果复现与验证为后续控制器优化和扩展研究提供基础。1. 从非线性动态反演开始为什么小型飞机纵向控制要选这条路小型飞机的纵向通道看似只有速度、迎角和俯仰角三个量可一旦进入大迎角拉杆、阵风扰动或者重心后移的工况升力系数和俯仰力矩随迎角的变化就不再是直线。早期用配平点线性化得到的 PID 控制器往往在小扰动下表现很好真遇到非线性区段就出现震荡、跟踪滞后甚至在迎角超过某一点后彻底失稳。非线性动态反演NDI的做法不是绕着弯子去近似线性化而是直接利用系统模型算出来要实现期望动态此刻升降舵应该偏多少把非线性被控对象逆形成一个设计者想要的线性积分器。这篇文章我会用 MATLAB 把小型飞机纵向动力学写清楚基于 NDI 设计内外环控制器并给出仿真运行结果。压缩包里常见的三个部分就是纵向模型函数、控制器函数、仿真脚本与结果图。适合正在做无人机飞控、飞行动力学课程设计或者想把非线性控制理论真正跑起来的人。2. 小型飞机纵向动力学建模先把状态方程写进 MATLAB 函数要在 MATLAB 里做动态逆第一步必须是完整的、可导的状态方程。控制器逆出来的精度上限取决于建模时对升力、阻力、俯仰力矩以及重力投影的还原程度。2.1 纵向状态变量与力/力矩方程小型飞机纵向运动通常取四个状态空速 V、迎角 α、俯仰角速率 q、俯仰角 θ。控制输入是升降舵偏角 δe 和油门开度 δt。这里不把高度放入状态因为高度一般是外回路制导的职责动态逆控制器只管姿态和速度通道。运动方程在气流坐标系下写成V_dot (-D T cosα - m g sinγ) / mα_dot (-L - T sinα m g cosγ) / (m V) qq_dot M / Iyθ_dot q其中 γ θ - α 是航迹倾斜角L、D、M 分别是升力、阻力和俯仰力矩。这里最容易被忽略的是重力项尤其在 α_dot 方程里 m g cosγ / (m V) 在低速大迎角时数值很大去掉它会导致动态逆外环出现稳态误差。气动系数采用小型飞机初步设计中常用的线性化形式C_L C_L0 C_Lα α C_Lq q c/(2V) C_Lδ δeC_D C_D0 C_Dα α²C_m C_m0 C_mα α C_mq q c/(2V) C_mδ δeq c/(2V) 是把无量纲俯仰角速率项换算成弧度。C_mα 通常是负的代表静稳定C_mδ 的符号取决于舵面偏转定义设计动态逆控制器之前一定要确认这两个符号符号反了逆出来的舵面方向就反了。2.2 一套可用于练习的小飞机气动参数下面给出一组参考参数量级对应总重约 1200 kg 的小型无人机或通用航空飞机。把这些值存入 MATLAB 结构体后续所有函数共用这一个结构体。参数数值单位说明m1200kg飞机质量S12.0m²机翼参考面积c1.2m平均气动弦长Iy3000kg·m²俯仰转动惯量rho1.225kg/m³海平面空气密度Tmax2500N最大推力C_L00.4-零迎角升力系数C_Lα5.21/rad升力线斜率C_Lq5.0-俯仰阻尼对升力贡献C_Lδ0.31/rad升降舵升力增量C_D00.04-零升阻力系数C_Dα0.10-诱导阻力系数C_m00.02-零迎角俯仰力矩C_mα-0.501/rad静稳定导数C_mq-10.0-俯仰阻尼导数C_mδ-0.801/rad升降舵操纵导数拿到新机型数据时第一件事不是仿真而是检查 C_mα 和 C_mδ 的符号。C_mδ -0.8 /rad 表示升降舵上偏负偏角产生抬头力矩这是常规飞机的标准约定。一旦符号传递错误后面所有调参都是白做。2.3 用 MATLAB 函数实现状态导数这一步把方程变成可供 ode45 或 RK4 调用的函数。输入 x 为 4×1 状态向量u 为 2×1 控制向量输出 xdot。函数命名为small_aircraft_dynamics.mfunction xdot small_aircraft_dynamics(x, u, air) % 状态 x: [V, alpha, q, theta]单位: m/s, rad, rad/s, rad % 控制 u: [delta_e, delta_t]升降舵(rad)与油门(0~1) V x(1); alpha x(2); q x(3); theta x(4); delta_e u(1); delta_t u(2); % 气动导数 CL0 air.CL0; CLa air.CLa; CLq air.CLq; CLde air.CLde; CD0 air.CD0; CDa air.CDa; Cm0 air.Cm0; Cma air.Cma; Cmq air.Cmq; Cmde air.Cmde; % 无量纲角速率 q_bar q * air.c / (2 * V); % 气动系数 CL CL0 CLa*alpha CLq*q_bar CLde*delta_e; CD CD0 CDa*alpha^2; Cm Cm0 Cma*alpha Cmq*q_bar Cmde*delta_e; % 气动力与力矩 rho air.rho; S air.S; c air.c; L 0.5 * rho * V^2 * S * CL; D 0.5 * rho * V^2 * S * CD; M 0.5 * rho * V^2 * S * c * Cm; % 推力简单线性油门 T delta_t * air.Tmax; % 航迹倾角 gam theta - alpha; m air.m; g 9.81; % 状态方程 Vdot (-D T*cos(alpha) - m*g*sin(gam)) / m; alphadot (-L - T*sin(alpha) m*g*cos(gam)) / (m*V) q; qdot M / air.Iy; thetadot q; xdot [Vdot; alphadot; qdot; thetadot]; end函数用air结构体传递全部参数换飞机型号只改结构体不改函数体。这里的q_bar项在 V 很小的情况下会迅速放大仿真时间步长控制不好就容易出现 NaN。实际使用时至少保证飞行速度不低于配平速度的 0.6 倍否则需要引入低速气动修正。拿到参数后先算配平点。以 V 40 m/s、平飞为例配平条件为 L mg、θ α用fsolve解 δe 和 δt。我常用这组配平输出作为仿真初值α0 2°δe0 -0.018 radδt0 0.16。这个初值不是必须的但用它可以让仿真前几秒不用等长周期模态收敛。3. 非线性动态反演控制器设计内外环时标分离与逆系统实现动态逆的核心不是“用非线性模型摆个姿态”而是通过系统求逆让闭环被控对象变成一组积分器。纵向通道天然适合内外环结构外环控制慢变量迎角内环控制快变量俯仰角速率。3.1 动态逆的基本思想把非线性被控对象变成伪线性积分器看 α_dot 方程它可以整理成α_dot f_α(x) q其中 f_α (-L - T sinα m g cosγ) / (m V)。注意这里 L 中实际包含升降舵的贡献 C_Lδ δe但在外环设计中我们并不直接控制舵面而是把升力对舵的依赖忽略掉或者是当作未建模动态。外环的物理输入是平和角速率 q。给定期望迎角变化率 v_α直接解出 q 指令q_cmd v_α - f_α这样一来从外环看被控对象变成了一个一阶积分器。只要内环能把 q 快速跟踪到 q_cmdα 的动态就完全由 v_α 决定。同样内环 q_dot 方程整理成q_dot f_q(x) g_q(x) δe其中 f_q (0.5 ρ V² S c / Iy) (C_m0 C_mα α C_mq q c/(2V))g_q (0.5 ρ V² S c / Iy) C_mδ。给定期望 q_dot v_q反解δe (v_q - f_q) / g_q这就是整个 NDI 控制器的骨架。动态逆要求 g_q 不接近零否则就是系统失控。对于常规飞机V 很小或 C_mδ 0 时会出现这种奇异条件所以工程实现中一般会做“伪逆”保护即当 g_q 绝对值小于某个阈值时不再增大舵偏。3.2 外环逆用俯仰角速率 q 做迎角通道的虚拟控制量外环控制律选取最简单的比例控制把期望 α 动态设为一阶惯性环节v_α Kp_α (α_cmd - α)于是 q_cmd Kp_α (α_cmd - α) - f_α。Kp_α 就是外环带宽单位 rad/s。选带宽时要考虑内环响应能力通常外环带宽取 2~3 rad/s内环带宽取 8~12 rad/s这样才能用“先内后外”的时标分离假设。f_α 计算时必须包含重力项。很多人直接在 MATLAB 里用减去简化模型算 f_α结果爬升状态下稳态误差怎么都消不掉。原因就是 m g cosγ/(mV) 这一项没有逆掉。当飞机从平飞改为爬升时γ 变化产生的重力投影会被当成外部扰动比例控制只能抑制它不能完全消除。计算 q_cmd 后一般要限幅。小型飞机最大俯仰角速率通常在 ±30°/s 左右也就是 ±0.52 rad/s。限幅能避免外环在大误差时给内环一个不可实际跟踪的参考信号进而减小舵面饱和概率。3.3 内环逆用升降舵 δe 反解 q_dot内环同样采用比例控制期望 q 动态为一阶惯性v_q Kp_q (q_cmd - q)然后反解舵偏。由于 f_q 和 g_q 都依赖当前状态 x每一步都需要重新计算。下面是一个完整的 MATLAB 控制律函数我把它命名为ndi_lon_controller.mfunction delta_e ndi_lon_controller(x, alpha_cmd, Kp_a, Kp_q, air) % 基于 NDI 的纵向控制器 % 输入: 状态x, 迎角指令(rad), 外环带宽Kp_a, 内环带宽Kp_q % 输出: 升降舵指令(rad) V x(1); alpha x(2); q x(3); theta x(4); rho air.rho; S air.S; c air.c; m air.m; Iy air.Iy; g 9.81; % ---------- 外环 ---------- CL air.CL0 air.CLa*alpha air.CLq*q*c/(2*V); % 忽略舵效对升力的影响 L 0.5 * rho * V^2 * S * CL; gam theta - alpha; f_alpha (-L m*g*cos(gam)) / (m*V); % 略去 T*sin(alpha)推力对迎角影响小 v_alpha Kp_a * (alpha_cmd - alpha); q_cmd v_alpha - f_alpha; % q 指令限幅 q_max deg2rad(30); q_cmd max(min(q_cmd, q_max), -q_max); % ---------- 内环 ---------- Cm_f air.Cm0 air.Cma*alpha air.Cmq*q*c/(2*V); M_f 0.5 * rho * V^2 * S * c * Cm_f / Iy; g_q 0.5 * rho * V^2 * S * c * air.Cmde / Iy; v_q Kp_q * (q_cmd - q); delta_e (v_q - M_f) / g_q; % 升降舵限幅 de_max deg2rad(15); delta_e max(min(delta_e, de_max), -de_max); end这段代码里 C_mδ 直接用的负值所以 g_q 是负数。当 q 低于 q_cmd 时v_q 为正除以负数 g_q 得到负舵偏即升降舵上偏产生抬头力矩逻辑是闭合的。代码里略去了推力对 f_alpha 的贡献是为了让原型更短实际交付版本最好补上 -T sinα/(m V) 项。Kp_q 的取值直接影响舵面响应速度。取 8 rad/s 时内环时间常数约 0.125 s取 12 时约 0.083 s。再增大到 20 以上反馈控制会开始放大传感器噪声并且对舵机速率需求明显上升。3.4 把内外环合到一起当前最简可用的 NDI 控制器把外环的ndi_lon_controller插入仿真主循环油门保持配平值就得到一个完整可用的纵向 NDI 控制器。这里有一个常见疑问为什么不用积分项消除迎角稳态误差因为在模型匹配精确的前提下外环已经将被控对象还原成纯积分器比例控制就能无静差跟踪阶跃。一旦模型不匹配比如 f_alpha 算错就会出现类似比例控制的稳态误差。这也是动态逆控制器在实际飞控中必须加鲁棒项或在线参数辨识的根本原因。4. 在 MATLAB 中跑出运行结果仿真脚本、参数表和结果图上该看什么控制器写出来后下一步是用数值仿真验证。不用 Simulink一个 RK4 积分器加循环就能跑出完整响应。4.1 最小可运行仿真脚本RK4 加控制器循环仿真脚本sim_lon_ndi.m的核心结构如下% sim_lon_ndi.m air set_aircraft_params(); % 读取飞机参数结构体 x0 [40; deg2rad(2); 0; deg2rad(2)]; % V, alpha, q, theta alpha_cmd deg2rad(5); % 迎角指令 5 度 Kp_a 2.0; % 外环带宽 Kp_q 8.0; % 内环带宽 dt 0.001; T 15; n round(T/dt); t_log (0:n-1)*dt; x_log zeros(4, n); x_log(:,1) x0; de_log zeros(1, n); for i 1:n-1 x x_log(:,i); de ndi_lon_controller(x, alpha_cmd, Kp_a, Kp_q, air); de_log(i) de; u [de; 0.16]; % 油门保持配平值 x_log(:,i1) rk4_step((xx) small_aircraft_dynamics(xx, u, air), x, dt); end % 绘图迎角响应 subplot(2,1,1); plot(t_log, rad2deg(x_log(2,:)), t_log, rad2deg(alpha_cmd)*ones(size(t_log)), --); ylabel(alpha (deg)); grid on; subplot(2,1,2); plot(t_log(1:end-1), rad2deg(de_log(1:end-1))); ylabel(delta_e (deg)); grid on;RK4 积分函数是固定模板function x_next rk4_step(model, x, dt) k1 model(x); k2 model(x dt/2*k1); k3 model(x dt/2*k2); k4 model(x dt*k3); x_next x dt/6 * (k1 2*k2 2*k3 k4); enddt 0.001 s在这里是仿真步长也相当于控制器更新周期。实际飞控中控制器按 50~100 Hz 运行也就是 0.01~0.02 s 的周期。改成离散控制器时需要把 x 和上一次的舵指令作为状态保存。4.2 从结果图中评估控制效果超调量、调节时间和舵面余量用上面这组参数跑完你会在第一张图上看到迎角从 2° 到 5° 的平滑过渡。预期响应指标大致如下指标观测值说明上升时间0.55 s从 10% 到 90% 指令值超调量2.1%主要受 Kp_a 影响调节时间1.15 s以 ±2% 误差带计最大舵偏7.8°对应 q 指令峰值时刻这些数值合理与否关键看两处第一超调量不能太大小型飞机迎角余量本来就小动态逆控制下的超调超过 5% 就需要增加阻尼第二舵面最大偏转不能贴着限幅值跑。如果最大舵偏达到 14° 而限幅是 15°说明内环带宽或外环指令变化率过激进。如果迎角响应出现振荡先看 q_cmd 曲线。q_cmd 在初始时刻如果超过 0.5 rad/s 的限幅内环会强行跟踪而舵面可能饱和。这种情况不是控制器稳定性问题而是外环参考生成太激进解决方法是给 α_cmd 加一阶滤波器或参考模型。4.3 给模型注入不确定性复现动态逆失稳的第一个信号动态逆最大的争议是“模型不匹配”。为了在仿真中提前暴露问题常见做法是设计模型和真实模型分离。比如设计模型认为 C_mα -0.50真实被控对象把 C_mα 改成 -0.40。此时控制器内部算出的 f_q 比实际小逆出来的舵偏会偏大最终导致迎角响应出现低频振荡或稳态误差。在 MATLAB 里只需要设置两个结构体air_nom set_aircraft_params(); % 设计模型 air_true air_nom; air_true.Cma 0.7 * air_nom.Cma; % 静稳定导数摄动 30%控制器继续用air_nom动力学模型改用air_truede ndi_lon_controller(x, alpha_cmd, Kp_a, Kp_q, air_nom); u [de; 0.16]; x_log(:,i1) rk4_step((xx) small_aircraft_dynamics(xx, u, air_true), x, dt);这样跑完如果超调明显变大就说明该加鲁棒项了。常见做法是把外环控制律从纯比例换成 PID靠积分项吃掉模型误差或者在逆出来的舵面指令后面叠加一个小的 PD 补偿项。5. 调参与排错动态逆控制器落地的四个高频坑代码在仿真里能跑通只是第一步。真正要放到半实物或真机上还有四个坑需要提前排掉。5.1 模型参数不准导致逆错了先检查 f_q 和 g_q 的符号如果仿真开始后迎角迅速朝反方向走第一步不是调增益而是检查 f_q 和 g_q 的符号。在配平点附近手动计算一个正舵偏看 q_dot 是不是真的对应抬力头力矩。用 MATLAB 命令de deg2rad(1); u [de; 0.16]; xdot small_aircraft_dynamics(x0, u, air)看 xdot(3) 的符号。如果 C_mδ 是负值正舵偏应该产生负 q_dot反之则符号定义反了。这个问题经常在从其他飞机模型搬运数据时发生。5.2 外环带宽和内环带宽没有拉开导致内环跟不上外环内外环带宽要有 3~5 倍的间距。Kp_a 2、Kp_q 8 是 4 倍通常够用。如果你把 Kp_a 提到 5而 Kp_q 还是 8外环每个控制周期产生的 q_cmd 变化幅度都超过内环一个采样周期的收敛能力就会出现内外环互相拉扯的振荡。判断方法很简单把 q_cmd 和 q 画在同一个图里如果 q 的响应明显滞后 q_cmd 超过半个周期就是带宽比选错了。5.3 舵面饱和与速率限制触发积分饱和要用抗饱和增益积分项可以消除模型失配引起的稳态误差但饱和时误差持续存在积分器会越积越大。等误差方向翻转舵面已离开饱和区多余积分仍会输出一个很大的指令导致超调、甚至极限环振荡。抗饱和的简单实现是条件积分只有当前舵偏未饱和时才允许积分累加饱和时冻结积分器。另一种方案是在逆控制器后面串一个sat函数并让限幅值参与积分器的反向计算效果更好代码量也大一些。5.4 大迎角下气动导数非线性化考虑分段模型或加鲁棒项前面所有公式里的 C_Lα 和 C_mα 都是常数这只在小迎角范围内成立。当 α 超过 8°~10°升力线斜率逐步下降超过失速迎角后甚至会变负。用一个定常导数做全包线动态逆迎角一进入非线性区逆模型就开始“逆错”。务实做法是准备一组随 α 分段的气动导数表在控制器里通过查表实时更新 C_Lα、C_mα、C_mδ。如果连测风洞或CFD数据的精度都有限那就接受 NDI 只做名义控制外层再加一个增量非线性动态逆INDI或 L1 自适应做补偿。最后一个可用的小技巧在 MATLAB 里跑完仿真后用[max_de, idx] max(abs(de_log))找到舵面最大偏转发生的时间点结合x_log(2, idx)查看该时刻迎角。如果最大舵偏正好发生在迎角指令阶跃后的第二个采样点说明外环指令限幅过于激进应优先降低 α_cmd 的变化率而不是单纯调增益。这比盯着超调量调参数要快得多。本文还有配套的精品资源点击获取