1. 先说结论：这个项目到底复现了什么
航天器姿态机动这个话题，做控制的人应该都不陌生。但一旦把“执行器饱和”和“执行器故障”同时摆上台面，问题就不是课本里那套线性 PID 能解决的了。我最近完完整整复现了一篇 IEEE 会议论文里的方案：基于四元数描述、采用自适应滑模结构的主动容错控制系统，用 Matlab 从模型搭建到控制器设计再到仿真验证全部走了一遍。整个过程踩了不少坑，也把论文里写得比较隐晦的地方逐一补全了。这篇文章就把我的复现思路、Matlab 实现细节和调试经验完整记录下来，给同样在做这个方向的人一些参考。
先说说这个系统解决什么问题。航天器在轨运行时，姿态机动靠的是飞轮、推力器这类执行机构。执行器有两个天生的麻烦：一是输出力矩有上限，也就是饱和；二是长期运行会出现性能退化甚至部分失效，也就是故障。传统的控制器设计如果完全忽略这两个因素，轻则机动时间拉长，重则姿态发散、任务失败。主动容错控制的核心思想，就是在控制器内部实时估计执行器的健康状态，在线调整控制增益和分配策略，让系统在故障发生后仍然能完成姿态机动目标。论文里的方法用自适应律在线估计执行器效率因子，再结合滑模控制的强鲁棒性，同时把饱和非线性用一个辅助动态系统来补偿。整体方案不需要故障诊断模块，不需要离线辨识，所有补偿都是在线完成的，工程实现上非常友好。
这篇博文面向的读者，我建议是已经有基础的控制理论基础、懂得状态空间和 Lyapunov 稳定性分析，但还没真正动手写过航天器姿态控制仿真代码的人。如果你只是刚接触 Matlab，建议先补一下 ode45 和函数句柄的基本用法。下面我按复现的实际流程来讲：从模型搭建到控制器设计，从代码实现到结果分析，最后是几个我实际踩过的坑。
2. 模型搭建：把物理问题变成数学方程
复现任何一篇控制论文，第一步不是写控制器，而是把被控对象模型搭出来。模型不对，后面所有工作都是白费。这一步看着简单，但其实有非常多的细节需要处理。
2.1 姿态运动学与动力学模型
航天器姿态描述方式有好几种，欧拉角、四元数、修正罗德里格斯参数（MRP）都有人用。这篇论文用的是四元数，原因很直接：全局无奇异。欧拉角在大角度机动时会遇到万向节锁死，而姿态机动恰恰是大角度运动，所以四元数是必然选择。
动力学方程用刚体欧拉方程描述：
Jω̇ = −ω×Jω + u + d
其中 J 是转动惯量矩阵，ω 是本体系相对惯性系的角速度在本体系下的表示，ω× 是叉乘矩阵，u 是执行器实际输出的控制力矩，d 是外部扰动力矩。
运动学方程用四元数微分方程：
q̇ = 0.5 * Ω(ω) * q
其中 Ω(ω) 是由角速度构成的 4×4 矩阵。这里有一个必须注意的点：四元数微分方程的积分结果不会自动保持模长为 1，而四元数的物理意义要求模长恒为 1。所以用 ode45 每积分一步，都要把四元数重新归一化。我一开始偷懒没做归一化，结果仿真到后面姿态直接飘了，这个问题后面会详细说。
仿真中我用的参数如下：
- J = diag([18, 22, 20]) kg·m²，模拟一个中等规模的近地卫星
- 初始四元数 q0 = [0.1; -0.15; 0.2; sqrt(1 - 0.01 - 0.0225 - 0.04)]，对应一个大约 30 度左右的初始姿态偏差
- 初始角速度 ω0 = [0.02; -0.03; 0.01] rad/s
- 期望姿态 q_d = [0; 0; 0; 1]，也就是机动到本体系与惯性系对齐
- 扰动力矩 d = 0.01 * [sin(0.5t); cos(0.3t); sin(0.8t)] N·m
期望角速度 ω_d = 0，也就是执行一个静止目标姿态机动。
2.2 执行器饱和与故障建模
执行器的建模是这次复现的关键。论文里把执行器特性分成了两个部分：非线性和不确定性。
先说饱和。每个执行轴的控制力矩输出不可能无限大，我用的是一个简单的饱和函数：
sat(u_i) = sign(u_i) * min(|u_i|, u_max_i)
也就是说，指令力矩超过上限就截断。这里注意饱和是非光滑非线性，会直接影响控制系统的稳定性分析，不能简单忽略。在我这个仿真里，每轴力矩上限设为 u_max = 4 N·m。
再说故障。论文研究的是执行器部分失效故障，模型为：
u_actual_i = ρ_i * u_command_i
其中 ρ_i ∈ (0, 1] 是执行器效率因子。ρ_i = 1 表示健康，ρ_i = 0.5 表示该轴只能输出一半的指令力矩。真实的执行器故障可能是突变也可能是缓变，论文里主要考虑突变场景，我也按照突变来建模：在 t = 15s 时，x 轴执行器效率从 1 突降到 0.5。
这里有一个容易被忽略的细节：故障和饱和是叠加在一起的。也就是说，执行器先对指令力矩做饱和处理，再乘以效率因子。建模顺序不能错，否则仿真结果对不上论文。
2.3 控制目标形式化
把物理需求翻译成数学语言，这个项目的控制目标可以写成三条：
- 系统状态全局一致最终有界（UUB），这是稳定性层面的要求
- 姿态误差四元数收敛到零附近的小邻域内，工程上对应姿态指向精度
- 在执行器饱和的前提下，系统在故障发生后仍然能完成姿态机动
用数学语言描述就是：设计控制输入 u，使得姿态误差 q_e → 0、角速度误差 ω_e → 0，同时保证所有闭环信号有界，并且控制力矩不超出执行器物理限制。
这里需要额外说明的是“γ-阶收敛”这个概念很多论文会提，但实际复现时更关心的是收敛速度和稳态精度。这两个指标往往互相制约，后面参数整定部分会详细讲。
3. 主动容错控制器设计：核心思想与推导过程
模型搭好了，接下来是控制器设计。这部分是论文的核心贡献，也是复现时最花时间的环节。我先把设计思路捋清楚，再给出数学推导，最后讲代码怎么实现。
3.1 为什么选滑模控制 + 自适应估计这个组合
先说滑模控制。航天器姿态机动是一个强非线性、强耦合的问题，而且存在外部扰动。滑模控制的优势在于对匹配不确定性（也就是作用在输入通道上的扰动和模型误差）具有天然的鲁棒性。设计一个好的滑模面，可以让系统状态在有限时间内到达滑模面，然后沿着滑模面滑动到原点。
但纯滑模控制有两个问题。第一个是抖振，符号函数导致的输入高频切换在实际系统中是不能接受的。第二个是它无法直接处理执行器故障——如果某个执行器只能输出一半力矩，滑模控制不会自动调整增益来补偿这个缺失。
所以论文引入了自适应估计。核心思想是把执行器效率因子 ρ 当作未知参数，设计一个自适应律在线估计它，然后用估计值重构控制增益。这样故障发生后，控制器能“感知”到执行器能力下降，自动加大健康通道的输出，实现容错。这就是“主动”二字的含义——不需要单独的故障诊断模块，一切都内嵌在控制器里。
这个组合的巧妙之处在于：滑模处理外部扰动，自适应处理执行器故障，两者各司其职，不会互相干扰。
3.2 滑模面设计与控制律推导
先定义姿态误差四元数。给定期望四元数 q_d，误差四元数通过四元数乘法计算：
q_e = q_d^{-1} ⊗ q
其中记 q_e = [q_ev; q_e0]，q_ev 是矢量部分，q_e0 是标量部分。角速度误差在这里就是 ω_e = ω，因为期望角速度为零。
滑模面选为：
s = ω_e + λ * q_ev
其中 λ 是正定对角矩阵，决定滑模面的收敛速率。这个滑模面的物理含义很直观：它让角速度误差和姿态误差成比例地收敛，避免出现姿态还没到位角速度已经很大的情况。
控制律设计为：
u_base = -ω×Jω - K_s * s - K_t * tanh(s/ε) + J * (-λ * q̇_ev)
逐项解释一下：
- -ω×Jω 是前馈补偿项，抵消陀螺力矩非线性项
- -K_s * s 是线性反馈项，保证到达滑模面的趋近速率
- -K_t * tanh(s/ε) 是鲁棒项，用双曲正切函数替代符号函数抑制抖振，ε 越大切换越平滑，但鲁棒性会有所下降
- J * (-λ * q̇_ev) 是模型补偿项，保证误差动态在滑模面上有理想的收敛特性
这里 q̇_ev 可以从运动学方程解析推导，不必用数值微分。注意故障时执行器实际输出是 ρ * sat(u_base)，如果忽略 ρ，控制效果会大打折扣，所以需要自适应环节。
3.3 自适应律设计与饱和补偿辅助系统
自适应律设计是论文的精华。设 ρ̂ 是 ρ 的估计值，估计误差定义为 ρ̃ = ρ - ρ̂。控制律中的鲁棒项增益会根据 ρ̂ 调整，自适应律采用投影算子保证估计值始终在物理可行范围内。
一个标准的设计是：
ρ̂̇_i = Proj[γ_i * s^T * L_i * sat(u_i)]
其中 L_i 是与执行器通道相关的映射矩阵，γ_i 是自适应增益，Proj 是投影算子，把 ρ̂ 限制在 [ρ_min, ρ_max] 内，我这里取 [0.1, 1]。ρ_min 的选取要小于执行器可能的最低效率，太保守会增加控制器负担，太乐观会导致鲁棒性不足。
投影算子的作用很关键。不加投影的话，自适应律完全由 Lyapunov 稳定性推导而来，理论上是收敛的，但实际仿真中数值误差和扰动会导致估计值偶尔跳出物理范围，比如变成负数。一旦 ρ̂ 为负，控制增益就会改变符号，整个系统直接发散。这在我前期调试中出现过好几次，后来加了投影算子就稳定多了。
饱和补偿方面，定义一个饱和差值：
Δu = sat(u) - u
当执行器饱和时 Δu ≠ 0，说明实际输入与设计输入不一致，系统会出现“追不上”的情况。处理方法是引入辅助动态系统：
ż = -A_z * z + B_z * Δu
其中 A_z 是正定矩阵，决定补偿的衰减速率。然后把辅助变量 z 引入滑模面：
s' = ω_e + λ * q_ev + k_z * z
这样当饱和发生时，z 会在线调整滑模面的位置，等效于让控制器提前“知道”执行器到顶了，从而避免饱和引起的过大超调和振荡。这个处理方式比直接限制指令力矩要优雅得多，也是论文里比较有工程价值的一个点。
4. Matlab实现：从主程序到每个模块的代码解析
控制系统设计完成后，接下来是编码实现。这部分的难点不在于某个语法，而在于如何把连续时间的微分方程、事件触发的故障注入、离散的控制律计算合理地组织在一起。
4.1 仿真框架与主程序结构
我采用的仿真框架是：用 ode45 积分闭环系统微分方程，在每个积分步内调用控制器函数计算控制力矩。控制器函数里包含执行器模型、饱和处理、自适应律更新和辅助系统动态。
主程序文件结构如下：
attitude_FTC_main.m % 主仿真脚本 attitude_dynamics.m % 系统状态微分方程（动力学+运动学+辅助系统） attitude_controller.m % 控制器：控制律+自适应律+饱和补偿 fault_injection.m % 故障事件注入 plot_attitude_results.m % 结果绘图主程序的关键代码框架：
%% 参数初始化 J = diag([18, 22, 20]); % 转动惯量矩阵 q0 = [0.1; -0.15; 0.2; ...]; % 初始四元数（需归一化） omega0 = [0.02; -0.03; 0.01]; % 初始角速度 qd = [0; 0; 0; 1]; % 期望姿态 % 控制器参数 lambda = 0.8 * eye(3); % 滑模面参数 K_s = 0.5 * eye(3); % 线性反馈增益 K_t = 0.1 * eye(3); % 鲁棒项增益 eps = 0.05; % 双曲正切平滑系数 gamma = 2.0; % 自适应增益 rho_init = [1; 1; 1]; % 效率因子初始估计 % 执行器参数 umax = 4; % 力矩饱和上限 A_z = 2 * eye(3); % 辅助系统衰减矩阵 B_z = eye(3); % 辅助系统输入矩阵 %% 仿真主循环（用ode45） [t, x] = ode45(@(t, x) attitude_dynamics(t, x, params), ... [0, 50], x0, options);这里我把所有参数打包到 params 结构体里，方便在仿真过程中修改和传递。实际运行中我发现用 ode45 的默认容差精度不够，建议设置options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8)，否则四元数的数值积分误差会在长时间仿真中累积。
4.2 系统状态微分方程与控制器函数
姿态动力学函数的实现如下：
function xdot = attitude_dynamics(t, x, params) % 状态向量 x = [q(4); omega(3); z(3); rho_hat(3)] q = x(1:4); omega = x(5:7); z = x(8:10); rho_hat = x(11:13); q = q / norm(q); % 四元数归一化关键步骤 % 计算控制力矩 u_cmd = attitude_controller(t, q, omega, z, rho_hat, params); % 执行器模型：饱和 + 故障 rho_true = params.rho_true; if t > params.fault_time rho_true = params.rho_fault; % 注入故障 end u_actual = rho_true .* saturate(u_cmd, params.umax); % 运动和动力学方程 Omega = quaternion_omega_matrix(omega); qdot = 0.5 * Omega * q; omega_dot = -cross(omega, params.J * omega) + u_actual + params.disturbance(t); omega_dot = params.J \ omega_dot; % 辅助系统动态 delta_u = saturate(u_cmd, params.umax) - u_cmd; zdot = -params.A_z * z + params.B_z * delta_u; % 自适应律更新（用投影算子） rho_hat_dot = projection_adaptive(t, q, omega, z, rho_hat, u_cmd, params); xdot = [qdot; omega_dot; zdot; rho_hat_dot]; end注意故障注入的实现：我没有在动力学方程里硬编码判断语句，而是通过rho_true这个变量在故障时刻切换。这样做的好处是仿真条件清晰，后期想改成缓变故障也很方便，只需把rho_true换成一条随时间变化的函数就行。
控制器函数就比较直接了：
function u_cmd = attitude_controller(t, q, omega, z, rho_hat, params) % 四元数误差计算 q_e = quaternion_error(q, params.qd); q_ev = q_e(1:3); % 期望角速度为零，角速度误差等于当前角速度 omega_e = omega; % 修正后的滑模面 s = omega_e + params.lambda * q_ev + params.k_z * z; % 控制律计算 s_omega = cross(omega, params.J * omega); nonlinear_comp = -s_omega + params.J * (-params.lambda * q_ev_dot + ...); linear_fb = -params.K_s * s; robust_term = -params.K_t * tanh(s / params.eps); % 自适应补偿（用效率因子估计值放大控制输出） u_nominal = nonlinear_comp + linear_fb + robust_term; u_cmd = (1 ./ rho_hat) .* u_nominal; % 这里要小心除法分母不为0 end这里有一个非常容易出错的地方：(1 ./ rho_hat)如果 ρ̂ 非常小，控制量会变得巨大。所以投影算子的下限必须设置合理，而且最好把除法形式改成乘法的控制结构，避免数值爆炸。这也是我踩过的一个大坑，后面详细说。
4.3 参数整定：哪些参数最关键
仿真跑通容易，跑出好看的结果难。在这一部分我花了很多时间调参数，总结下来最关键的几个：
第一个是 λ，也就是滑模面的比例系数。λ 越大，姿态误差收敛越快，但等效控制量也会增大，更容易触发饱和。我的整定经验是从小往大调，先设 0.3 看趋势，再逐步增加到 0.8。如果发现角速度峰值过大或者饱和时间过长，就退回去 0.6 左右。
第二个是 K_t 和 ε 的组合，这一对参数直接决定抖振水平。K_t 大、ε 小，鲁棒性最好但抖振剧烈；K_t 小、ε 大，控制量平滑但稳态精度变差。我最后的取值是 K_t = 0.1、ε = 0.05，在这个参数下抖振幅度在可接受范围内，姿态稳态误差小于 0.1 度。
第三个是自适应增益 γ。这个参数最微妙，我甚至可以说它是整个系统调参中最考验耐心的一环。γ 太小，故障后恢复很慢，可能要几十秒才能估出真实效率；γ 太大，估计值会剧烈波动，进而驱动控制量高频变化，系统可能直接失去稳定。我试过 γ = 0.5，恢复时间超过 20 秒；γ = 5 时系统低频振荡明显。最后定在 γ = 2，大约 5 秒内能收敛到真实值的 90%。
第四个是辅助系统的时间常数 A_z。A_z 越大，饱和补偿越快，但 z 变量会出现高频成分；A_z 太小，补偿滞后明显，饱和期间姿态误差会持续增大。取 A_z = 2 时表现比较均衡。
5. 仿真结果分析：故障发生前后系统表现如何
参数整定完成以后，我跑了一组完整的仿真，并对比了有容错和无容错两种情况，看看主动容错的优势到底有多大。
5.1 姿态机动过程的动态响应
先看无故障时的情况。系统从初始姿态偏差出发，控制器快速建立机动力矩，姿态误差四元数前 5 秒平滑收敛，大约 8 秒达到稳态，稳态姿态误差在 0.1 度以内。角速度峰值大约 3.5 度每秒，在执行器饱和范围内，没有出现明显的饱和时段。控制力矩在初始阶段短暂达到上限 4 N·m，随后迅速回落到 0.5 N·m 以下，整体表现平稳。
这个结果说明，在不考虑故障的理想情况下，控制器具备良好的机动性能和稳态精度。滑模面在 1 秒内到达，之后的运动沿着滑模面渐进趋近期望姿态，没有明显的超调和振荡。
5.2 故障注入后的容错表现
t = 15s 时，x 轴执行器效率突降为 0.5。这里可以看到容错控制的关键作用。
没有容错设计的常规滑模控制器，在故障发生后姿态开始漂移，因为 x 轴实际输出力矩只有指令的一半。姿态误差从原来的 0.1 度逐渐增大到 0.8 度左右，而且持续发散的趋势明显，角速度也出现了持续的非零偏差。这说明一个事实：常规滑模对参数不确定性和扰动有鲁棒性，但对执行器效率损失这种输入通道乘性故障无能为力。
而采用主动容错方案的控制器，在故障瞬间出现了短暂的姿态扰动，最大偏差大约 1.5 度，这是不可避免的——故障发生的信息必须通过系统响应才能被感知到。随后自适应律开始调整 ρ̂_x，大约 5 秒内从 1 降到 0.55 左右（略偏向保守），控制器补偿增益相应升高，姿态在故障后 3 秒内恢复到 0.2 度以内，最终回稳到 0.1 度附近。
值得注意的是，故障后控制器会主动增加 x 轴的指令力矩。这是因为 ρ̂ 变小时，1/ρ̂ 变大，相当于用更大的指令去弥补执行器的缺失。仿真中 x 轴控制力矩在故障后上升到 1.5 N·m 左右，而 y、z 轴则承担了部分耦合补偿，整体未触发饱和。
5.3 和 IEEE 原论文结果对比时的注意事项
复现论文，最后总要对标一下原论文的仿真结果。但这里我想提醒几个容易踩的坑。
第一是单位问题。论文里角速度可能用 rad/s，也可能用 deg/s，姿态误差可能用四元数分量画图，也可能改成欧拉角画图。我见过不少人在对比时因为单位不一致得出错误结论。建议统一采用 rad/s 和度两个维度分别展示，方便对图。
第二是时间尺度。不同论文的转动惯量、力矩上限不一样，收敛时间自然没有可比性。对比的重点应该放在控制趋势、故障后的恢复形态和稳态精度上，而不是苛求曲线完全重合。
第三是扰动的设置。有些论文为了突出控制器性能会把扰动设得很小，复现时如果你加了更大的扰动，结果变差是正常的，不一定是控制器实现有误。
综合来看，我的复现结果在定性层面和论文结论一致：故障后系统能在有限时间内恢复到期望姿态附近，姿态误差有界，控制量在执行器物理范围内。定量层面，由于参数差异，收敛时间略有不同，但整体规律相同。
6. 踩坑记录：我实际复现时遇到的那些问题
这一部分是我最想分享的。控制算法从论文到代码，中间隔着大量的实现细节，每一个细节不到位，仿真结果就是不对。下面是我实际遇到的一些问题，按照从基础到进阶的顺序整理。
6.1 四元数归一化与数值积分漂移
这是我遇到的第一个“看不见”的错误。四元数微分方程是线性的，用 ode45 积分没有任何数值问题，但物理上四元数必须满足单位模长约束。ode45 不会自动满足这个约束，积分几步后模长就会偏离 1。起初偏离很小，看起来无伤大雅；但 50 秒仿真结束后，四元数模长可能已经变成 1.005 甚至更大，反映到姿态误差上就是持续的小偏差。
解决办法是在每个积分步内强制归一化。注意不能在积分函数里直接修改状态（ode45 要求状态连续），正确做法是在计算完 qdot 后把 q 归一化再赋值给 xdot 的输出，也就是状态导数层面的修正。更正规的做法是用四元数归一化重新计算运动学，但我测试发现直接归一化状态导数在工程上已经足够稳定。
6.2 投影算子的下限到底怎么定
我在 3.3 节提到过，自适应律不加投影会导致 ρ̂ 变为负数。这是第一层坑。第二层坑是投影下限的选择。我把下限设为 0.1，仿真中发现自适应收敛速度变慢，因为 ρ̂ 需要从 1 一路降到 0.55，如果下限太紧，梯度更新会被频繁截断。
后来我改用了一种更平滑的方法：不在自适应律里直接截断，而是在计算控制量时对 ρ̂ 做一个迟滞处理，限制它的下降速度。这样既不破坏 Lyapunov 框架，又能避免数值振荡。具体做法是给 ρ̂̇ 加一个一阶低通滤波，相当于在估计值层面做了平滑。这个技巧论文里通常不会写，但对实际仿真稳定性帮助很大。
6.3 滑模抖振与 ε 参数的博弈
滑模控制的抖振问题在仿真中表现得非常直接：控制力矩曲线出现高频毛刺，姿态角速度出现微小但明显的高频振荡。双曲正切函数的 ε 参数是平衡点，但不同通道的最优 ε 可能不同。
我在调试中发现，把 ε 取为一个固定的对角矩阵（三个通道取不同值），收敛效果比所有通道都用同一个标量更好。物理原因不难理解：不同轴的转动惯量不同，惯量小的轴对控制更敏感，需要更大的 ε 来平滑；惯量大的轴则可以取较小的 ε 来保持鲁棒性。这个细节论文通常不会展开写，但实际调参体验差异非常大。
6.4 执行器饱和时自适应律不更新的问题
这个问题比较隐蔽。当执行器饱和时，Δu = sat(u) - u ≠ 0，辅助系统被激活，滑模面被修正。但自适应律是根据 sat(u) 来更新的，而 sat(u) 是截断后的力矩。如果在饱和期间自适应律继续按正常逻辑更新，ρ̂ 会朝错误方向移动，导致饱和结束后控制器参数已经偏离理想值。
我采用的变通方案是：在自适应律里加入一个饱和检测标志，当某轴 Δu 的幅值超过阈值时，该轴的 ρ̂ 停止更新。这个“暂停机制”在工程上非常有效，饱和结束后自适应律重新激活，参数不会出现明显的突变。
7. 个人体会与建议
整个复现过程下来，我的最大体会是：论文里的公式只是起点，真正的工程问题全在公式以外的细节里。滑模参数、自适应增益、投影下限、饱和补偿的衰减速率，每一个参数都和其他参数耦合在一起，不可能一次性调好。我的建议是先跑通一个不含故障、不含饱和的理想模型，确认控制器基本功能正常；再加入饱和，调试辅助系统参数；最后加入故障，观察自适应律的表现。分阶段调试能极大减少排查问题的范围。
另外一个建议是，善用 Matlab 的实时脚本（Live Script）调试。控制律迭代过程中经常需要反复修改参数，实时脚本可以同时展示代码、仿真曲线和中间变量，方便快速定位问题。我最后把整个仿真框架整理成了带注释的版本，后续换一组参数、换一种故障模式，只需修改初始化部分就能直接复用。
最后想说的是，复现论文从来不是对着公式敲代码那么简单。你要理解每个设计选择背后的取舍，处理公式里没写的实现细节，还要接受结果跟论文不完全一致的现实。但正是这些过程，让一次复现变成真正的能力提升。如果这篇文章能帮你少走一点弯路，那就值了。
