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调试。控制律迭代过程中经常需要反复修改参数实时脚本可以同时展示代码、仿真曲线和中间变量方便快速定位问题。我最后把整个仿真框架整理成了带注释的版本后续换一组参数、换一种故障模式只需修改初始化部分就能直接复用。最后想说的是复现论文从来不是对着公式敲代码那么简单。你要理解每个设计选择背后的取舍处理公式里没写的实现细节还要接受结果跟论文不完全一致的现实。但正是这些过程让一次复现变成真正的能力提升。如果这篇文章能帮你少走一点弯路那就值了。
