企业官网建设流程全解析

简介：这套基于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) + q
q_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δ δe
C_D = C_D0 + C_Dα α²
C_m = C_m0 + C_mα α + C_mq q c/(2V) + C_mδ δe

q c/(2V) 是把无量纲俯仰角速率项换算成弧度。C_mα 通常是负的，代表静稳定；C_mδ 的符号取决于舵面偏转定义，设计动态逆控制器之前一定要确认这两个符号，符号反了逆出来的舵面方向就反了。

2.2 一套可用于练习的小飞机气动参数

下面给出一组参考参数，量级对应总重约 1200 kg 的小型无人机或通用航空飞机。把这些值存入 MATLAB 结构体，后续所有函数共用这一个结构体。

拿到新机型数据时，第一件事不是仿真，而是检查 C_mα 和 C_mδ 的符号。C_mδ = -0.8 /rad 表示升降舵上偏（负偏角）产生抬头力矩，这是常规飞机的标准约定。一旦符号传递错误，后面所有调参都是白做。

2.3 用 MATLAB 函数实现状态导数

这一步把方程变成可供 ode45 或 RK4 调用的函数。输入 x 为 4×1 状态向量，u 为 2×1 控制向量，输出 xdot。函数命名为small_aircraft_dynamics.m：

function 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.m：

function 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(:,i+1) = 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); end

dt = 0.001 s在这里是仿真步长，也相当于控制器更新周期。实际飞控中控制器按 50~100 Hz 运行，也就是 0.01~0.02 s 的周期。改成离散控制器时，需要把 x 和上一次的舵指令作为状态保存。

4.2 从结果图中评估控制效果：超调量、调节时间和舵面余量

用上面这组参数跑完，你会在第一张图上看到迎角从 2° 到 5° 的平滑过渡。预期响应指标大致如下：

这些数值合理与否，关键看两处：第一，超调量不能太大，小型飞机迎角余量本来就小，动态逆控制下的超调超过 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_true：

de = ndi_lon_controller(x, alpha_cmd, Kp_a, Kp_q, air_nom); u = [de; 0.16]; x_log(:,i+1) = 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 的变化率而不是单纯调增益。这比盯着超调量调参数要快得多。

本文还有配套的精品资源，点击获取