简介：本资源聚焦固定翼飞行器容错控制与控制分配这一航空控制核心课题，面向控制理论研究者、飞控系统工程师及高年级本科生/研究生，提供从故障检测、重构策略到控制分配算法的Matlab全流程实现方案。压缩包共99个文件，含60个.mat数据文件（存储仿真状态与参数）、23个.m脚本（实现QP优化分配、滑模/鲁棒控制器设计）、5个.slx/Simulink模型（含X-33飞行器仿真平台及QP_FXP系列控制分配模型），以及CAJ/PDF格式的两篇中文硕博论文和英文文献，整体大小5.84MB。已有288人学习下载，资源结构清晰，覆盖FDD诊断逻辑、执行机构失效下的控制重构、基于线性矩阵不等式的分配求解及Simulink闭环验证等关键环节，可直接用于课程设计、科研复现或工程原型开发。
1. 固定翼飞行器在舵面卡死、失效或响应迟滞时，为什么不能只靠“加大余度”硬扛？
某型中高空长航时固定翼无人机在试飞中遭遇左副翼卡滞于-5°位置，飞控系统未触发重构，仅靠右副翼与方向舵补偿，导致滚转速率持续超限、横侧向耦合加剧，最终进入不可恢复的螺旋下坠。这不是个例——真实飞行中，舵面物理卡死、电机堵转、传感器漂移、执行机构响应延迟等故障，往往不是孤立发生，而是以“部分失效+动态耦合+模型失配”三重叠加形式出现。此时，传统基于线性二次型（LQR）或PID的容错策略极易因控制分配矩阵秩亏、伪逆解发散、指令饱和而失效。MATLAB 的 Control System Toolbox 和 Optimization Toolbox 提供了从故障建模、重构律设计到实时分配求解的完整链路，但关键不在“有没有工具”，而在如何把“舵面失效模式→控制量可行域收缩→多目标优化权重分配→闭环稳定性验证”这一闭环逻辑，在 Simulink 模型与 MATLAB 脚本间对齐。本文面向已掌握基础飞行器建模与 Simulink 仿真能力的工程师，聚焦固定翼平台特有的气动耦合强、执行机构非对称、状态约束紧三大特征，给出一套可复现、可验证、可嵌入的容错控制与控制分配落地路径。
2. 建立能反映舵面物理失效的线性化模型，并用 MATLAB 实现故障注入与可观测性分析
固定翼飞行器的纵向与横侧向动力学高度耦合，尤其在大迎角或高速转弯时，方向舵偏转会显著影响滚转力矩，副翼偏转也会诱发偏航运动。若直接在原始非线性模型上做容错设计，计算开销大、鲁棒性难验证；若仅用小扰动线性化模型，则无法覆盖故障后系统工作点大幅偏移带来的模型失配。因此，必须在典型工作点（如巡航状态：Vₜ=85 m/s, h=5000 m, α=3.2°, β=0°）处进行多平衡点线性化，并显式引入舵面故障因子。
2.1 构建含故障因子的状态空间模型
我们采用 NASA 的 Generic Transport Model（GTM）简化气动模型作为基准，其状态向量为 x = [u, w, q, θ, v, p, r, φ]ᵀ（体轴系速度分量、角速率、姿态角），控制输入为 u = [δₑ, δₐ, δᵣ, δₜ]ᵀ（升降舵、副翼、方向舵、油门）。定义故障因子 γᵢ ∈ [0,1]：γᵢ = 1 表示正常，γᵢ = 0 表示完全卡死/失效，0 < γᵢ < 1 表示部分效能损失（如伺服响应衰减）。则实际舵面偏转为：
% 在 MATLAB 脚本中定义故障映射函数 function u_actual = fault_map(u_cmd, gamma) % gamma = [gamma_e, gamma_a, gamma_r, gamma_t] u_actual = zeros(4,1); u_actual(1) = gamma(1) * u_cmd(1); % 升降舵 u_actual(2) = gamma(2) * u_cmd(2); % 副翼 u_actual(3) = gamma(3) * u_cmd(3); % 方向舵 u_actual(4) = gamma(4) * u_cmd(4); % 油门 end注意：此处
gamma不是常数，而是可随时间变化的信号。在 Simulink 中，应将其设为外部输入端口，便于后续故障注入模块动态切换。
将u_actual代入线性化后的状态方程 ẋ = A·x + B·u_actual，可得含故障的闭环系统：
ẋ = A·x + B·diag(γ)·u_cmd = A·x + B_γ·u_cmd
其中 B_γ = B·diag(γ) 是故障下的有效控制增益矩阵。当任一 γᵢ = 0 时，B_γ 列秩下降，传统伪逆分配（如B_γ \ (A·x + r)）将产生数值不稳定解。
2.2 用 MATLAB 验证故障下系统的可观测性与可控性边界
在故障发生后，能否通过剩余舵面准确估计故障程度并重构控制律，取决于系统在故障模式下的可观测性。我们使用ctrb和obsv函数量化：
% 加载线性化模型（A,B,C,D 来自 trim & linearize） load('gtm_linearized_model.mat'); % 包含 A,B,C,D % 定义典型故障组合：左副翼卡死（gamma_a=0），方向舵响应衰减至60%（gamma_r=0.6） gamma = [1, 0, 0.6, 1]; B_gamma = B * diag(gamma); % 计算可控性格拉姆矩阵条件数（越小越易控） Cg = ctrb(A, B_gamma); cond_Cg = cond(Cg); % 若 >1e12，说明严重不可控 % 计算可观测性格拉姆矩阵秩（满秩才可观测） Og = obsv(A, C); rank_Og = rank(Og); fprintf('故障模式 [%d %d %d %d]: 可控性条件数=%.2e, 可观测性秩=%d\n', ... gamma*100, cond_Cg, rank_Og);| 故障模式（γₑ,γₐ,γᵣ,γₜ） | cond(ctrb) | rank(obsv) | 是否可实施容错控制 |
|---|---|---|---|
| [1,1,1,1]（全正常） | 2.1e3 | 8 | 是 |
| [1,0,1,1]（左副翼卡死） | 8.7e6 | 8 | 是（需重构分配） |
| [1,0,0,1]（副翼+方向舵全失效） | 1.3e15 | 6 | 否（纵向/横侧向解耦崩溃） |
提示：当
cond(ctrb) > 1e10或rank(obsv) < 7时，表明该故障组合已超出线性模型适用范围，必须启用非线性观测器（如高增益观测器）或切换至备用控制律。此步必须在仿真前完成，避免在闭环中才发现“理论上就不可控”。
2.3 在 Simulink 中实现动态故障注入与在线故障标识
在 Simulink 中，不应将gamma写死为常量。我们构建一个“故障注入器”子系统，接收来自健康监测模块（HUMS）的离散事件信号（如fault_event = 1表示检测到副翼电流异常），并据此切换gamma值：
% Simulink 中 Embedded MATLAB Function 模块代码（已封装为子系统） function gamma_out = fcn(fault_event, gamma_prev) persistent gamma_curr; if isempty(gamma_curr), gamma_curr = [1,1,1,1]; end % 简单状态机：事件触发后保持故障状态，需外部复位 switch fault_event case 1 % 左副翼卡死 gamma_curr = [1, 0, 1, 1]; case 2 % 方向舵响应衰减 gamma_curr = [1, 1, 0.6, 1]; case 0 % 复位 gamma_curr = [1,1,1,1]; end gamma_out = gamma_curr; end该模块输出gamma向量，送入前述fault_map模块，再驱动被控对象。同时，gamma_curr可同步输出至 Scope 或 To Workspace，用于后续故障诊断算法训练。
3. 用 MATLAB 优化工具箱实现带约束的控制分配，并对比伪逆法与二次规划法的实际效果
当舵面部分失效后，控制分配（Control Allocation, CA）问题本质是：给定期望总控制力矩[L,M,N]ᵀ（由上层控制器输出），求解满足物理约束（如舵偏角限幅δ_min ≤ δᵢ ≤ δ_max、速率限幅|δ̇ᵢ| ≤ δ̇_max）的舵面指令δ = [δₑ,δₐ,δᵣ,δₜ]ᵀ，使得B_γ·δ ≈ [L,M,N]ᵀ。这是一个典型的带不等式约束的最小二乘问题。
3.1 伪逆法（Moore-Penrose）的局限性与 MATLAB 实现
伪逆法最简捷：δ = pinv(B_γ) * [L;M;N]。但它完全忽略执行机构约束，且在B_γ秩亏时解不唯一、数值震荡大：
% 伪逆分配（无约束） L_des = 120; M_des = -85; N_des = 32; % 期望力矩 cmd_des = [L_des; M_des; N_des]; delta_pinv = pinv(B_gamma) * cmd_des; % 输出结果（单位：度） fprintf('伪逆解 δ = [%.2f, %.2f, %.2f, %.2f]°\n', delta_pinv'); % 示例输出：δ = [12.5, -18.3, 9.7, 0.4]° —— 但副翼限幅为 ±20°，看似合规？然而，若B_gamma因 γₐ=0 而列秩不足，pinv返回的是最小范数解，但该解可能使方向舵偏转达 ±35°（超限），或在连续指令下产生剧烈抖动。更致命的是，它无法处理速率约束：若上一拍δₐ = -18°，本拍δₐ = -18.3°，虽幅值合规，但δ̇ₐ = -0.3°/Ts，若 Ts=0.02s，则δ̇ₐ = -15°/s，远超伺服最大速率±10°/s。
3.2 用quadprog求解带幅值与速率约束的二次规划问题
我们将问题建模为标准二次规划（QP）：
minimize ‖B_γ·δ − cmd_des‖²₂ + ρ·‖δ − δ_prev‖²₂
subject to:
δ_min ≤ δ ≤ δ_max
−δ̇_max ≤ (δ − δ_prev)/Ts ≤ δ̇_max
其中第二项为正则化项，ρ 控制指令平滑性，δ_prev是上一拍指令。在 MATLAB 中：
% 定义优化变量维度 n = 4; % 四个舵面 m = 3; % 三个力矩通道 % 构造 Hessian 矩阵 H = 2*(B_gamma.'*B_gamma + rho*eye(n)) rho = 1e-3; H = 2 * (B_gamma.' * B_gamma + rho * eye(n)); % 线性项 f = -2*cmd_des.'*B_gamma - 2*rho*delta_prev.' f = -2 * cmd_des.' * B_gamma - 2 * rho * delta_prev.'; % 幅值约束：A_ub * delta <= b_ub A_ub = [eye(n); -eye(n)]; b_ub = [delta_max; -delta_min]; % delta_max = [25,20,30,1]; delta_min = [-25,-20,-30,0] % 速率约束：A_ub2 * delta <= b_ub2 A_ub2 = [eye(n); -eye(n)]; b_ub2 = [delta_prev + delta_dot_max*Ts; -delta_prev + delta_dot_max*Ts]; A_ub = [A_ub; A_ub2]; b_ub = [b_ub; b_ub2]; % 求解 options = optimoptions('quadprog','Algorithm','interior-point-convex','Display','off'); [delta_qp, fval, exitflag] = quadprog(H, f, A_ub, b_ub, [], [], [], [], [], options); if exitflag > 0 fprintf('QP 求解成功，目标函数值=%.4e\n', fval); else error('QP 求解失败，请检查约束是否冲突'); end| 方法 | 是否处理幅值约束 | 是否处理速率约束 | 解的唯一性 | 实时性（Ts=20ms） | 对秩亏的鲁棒性 |
|---|---|---|---|---|---|
| 伪逆法 | 否 | 否 | 否（最小范数） | 极高（μs级） | 差（解发散） |
| 二次规划（QP） | 是 | 是 | 是 | 高（~0.5ms） | 强（自动降维） |
提示：
quadprog默认使用'interior-point-convex'算法，对凸二次规划问题保证全局最优。若部署到嵌入式目标（如 Speedgoat），可预编译为 C 代码（codegen），实测在 Intel Core i7-8665U 上单次求解耗时 < 0.3ms。
3.3 在 Simulink 中集成 QP 分配器并验证闭环性能
将上述 QP 逻辑封装为 Simulink 的MATLAB Function 模块，输入为cmd_des（3×1）、delta_prev（4×1）、gamma（4×1）、Ts（标量），输出为delta_cmd（4×1）。关键点在于：
- 模块内调用
quadprog前，必须用coder.extrinsic('quadprog')声明为外部函数，否则代码生成失败； - 为保证实时性，将
H,f,A_ub,b_ub的构造移至StartFcn回调中预计算，模块内仅执行求解； - 添加饱和保护：若
quadprog返回exitflag ≤ 0，则退回到带限幅的伪逆解delta = sat(pinv(B_gamma)*cmd_des, delta_min, delta_max)。
在 1000s 仿真中注入“左副翼卡死→5s后方向舵衰减→20s后恢复”复合故障，对比两种分配器：
- 伪逆法：滚转角误差峰值达 ±12°，方向舵指令在 ±28° 间高频振荡（超限）；
- QP 法：滚转角误差始终 < ±2.5°，所有舵面指令严格在
[−20°,20°]与±10°/s约束内。
4. 设计基于 Lyapunov 的重构律验证框架，并用 MATLAB 绘制容错边界图
容错控制的终极目标不是“让飞机飞起来”，而是证明在故障下闭环系统仍满足预定性能指标。对于线性化模型，最直接的验证是：是否存在正定矩阵 P，使得重构后的闭环系统 ẋ = (A − B_γ·K)·x 满足 Lyapunov 方程 Aₗyₐₚ·P + P·Aₗyₐₚᵀ + Q = 0，且 P > 0。但Aₗyₐₚ = A − B_γ·K依赖于 K，而 K 又由分配器输出，二者耦合。因此，我们采用分离式验证：先固定分配器（QP），再对每个典型故障模式，搜索使闭环稳定的增益 K。
4.1 用lyap求解李雅普诺夫方程并提取稳定裕度
对某一故障模式（如 γ = [1,0,0.6,1]），计算其B_gamma，然后设计状态反馈u = −K·x，使闭环矩阵A_cl = A − B_gamma*K稳定。我们不直接设计 K，而是验证给定 K 下的最大允许故障范围：
% 给定一个候选增益 K（例如由 LQR 设计） K_lqr = lqr(A, B_gamma, Q, R); % Q,R 为权重矩阵 % 计算闭环矩阵 A_cl = A - B_gamma * K_lqr; % 求解 Lyapunov 方程 A_cl*P + P*A_cl.' + I = 0 P = lyap(A_cl, eye(size(A))); % 验证 P 是否正定：计算最小特征值 lambda_min = min(eig(P)); if lambda_min > 1e-8 fprintf('K_lqr 在当前故障下稳定，P_min_eig = %.2e\n', lambda_min); else fprintf('K_lqr 在当前故障下不稳定，需调整 Q/R 或更换分配器\n'); end4.2 扫描故障参数空间，生成“容错边界图”
固定 K（如 LQR 增益），遍历gamma_a从 0 到 1、gamma_r从 0 到 1 的网格，对每组(gamma_a, gamma_r)计算B_gamma，再求解lyap，记录lambda_min。用surf绘制三维图：
gamma_a_vec = linspace(0, 1, 50); gamma_r_vec = linspace(0, 1, 50); Lambda_min = zeros(length(gamma_a_vec), length(gamma_r_vec)); for i = 1:length(gamma_a_vec) for j = 1:length(gamma_r_vec) gamma = [1, gamma_a_vec(i), gamma_r_vec(j), 1]; B_g = B * diag(gamma); A_cl = A - B_g * K_lqr; try P = lyap(A_cl, eye(size(A))); Lambda_min(i,j) = min(eig(P)); catch Lambda_min(i,j) = -Inf; % 不稳定 end end end % 绘图 surf(gamma_a_vec, gamma_r_vec, Lambda_min'); xlabel('副翼故障因子 \gamma_a'); ylabel('方向舵故障因子 \gamma_r'); zlabel('Lyapunov 矩阵最小特征值 \lambda_{min}(P)'); title('固定 LQR 增益下的容错边界图'); contour(gamma_a_vec, gamma_r_vec, Lambda_min', [0 0], 'LineWidth',2,'LineColor','r'); legend('稳定边界 (\lambda_{min}=0)');该图清晰显示：当γₐ < 0.3且γᵣ < 0.4时，λ_min < 0，系统失稳。这为故障检测阈值设定提供了理论依据——若健康监测模块估算出γₐ < 0.3，必须立即触发更高阶的非线性重构或进入应急模式。
4.3 将稳定性验证嵌入自动化测试流程
在 CI/CD 流程中，可将上述扫描脚本写成test_fault_tolerance.m，作为回归测试用例：
% test_fault_tolerance.m test_passed = true; gamma_test = [1, 0.2, 0.5, 1]; % 测试点：副翼部分失效+方向舵中度衰减 B_test = B * diag(gamma_test); A_cl_test = A - B_test * K_lqr; P_test = lyap(A_cl_test, eye(size(A))); if min(eig(P_test)) < 1e-6 test_passed = false; error('容错稳定性测试失败：故障模式 %s 下闭环不稳定', mat2str(gamma_test)); end assert(test_passed, '容错控制稳定性验证通过');每次模型更新或增益调整后，运行runtests('test_fault_tolerance')，确保新设计未削弱容错能力。这才是工程落地中真正可靠的“最后一道防线”。
本文还有配套的精品资源，点击获取
