1. 项目背景与核心挑战
航天器姿态控制是空间任务中最关键的技术之一。在实际太空环境中，执行机构（如反作用轮、推力器等）常面临两种典型问题：一是物理饱和（执行力矩达到硬件极限），二是突发故障（如部分失效或完全失效）。这两种情况都会导致控制系统性能急剧下降，甚至引发任务失败。传统控制方法往往将饱和与故障分开处理，而实际工程中二者经常同时出现，这就对控制系统的鲁棒性和容错能力提出了更高要求。
我最近复现的这篇TIE（IEEE Transactions on Industrial Electronics）论文，提出了一种融合状态观测器、反步控制（Backstepping）和自适应滑模控制的主动容错方案。这套方法最吸引我的地方在于：它不需要故障检测与诊断（FDD）系统提供精确的故障参数，而是通过在线自适应机制实时补偿故障和饱和的综合影响。下面我将从工程实现角度，拆解这个系统的设计要点和Matlab实现中的关键技巧。
2. 系统建模与问题描述
2.1 航天器姿态动力学模型
采用修正罗德里格斯参数（MRP）描述姿态，其动力学模型可表示为：
matlab复制% 航天器姿态动力学方程
function dsigma = MRP_dynamics(sigma, omega, J, u, d)
S = skew(sigma);
J_omega = J * omega;
dsigma = 0.25 * ((1 - sigma'*sigma)*eye(3) + 2*S + 2*sigma*sigma') * omega;
domega = inv(J) * (-skew(omega) * J_omega + u + d);
end
function S = skew(v)
S = [0 -v(3) v(2); v(3) 0 -v(1); -v(2) v(1) 0];
end
其中sigma为MRP向量，omega为角速度，J为惯量矩阵，u为控制输入，d为外部扰动。执行器饱和用饱和函数sat(u)表示，故障模型采用乘性+加性形式：
code复制u_actual = ρ(t) * sat(u) + φ(t)
ρ(t)为执行器效率矩阵（对角阵，元素∈[0,1]），φ(t)为偏差故障。
2.2 控制目标与技术路线
核心控制目标：
- 在存在执行器饱和和未知故障情况下，实现姿态稳定跟踪
- 不依赖精确的故障诊断信息
- 抑制外部扰动影响
技术路线分三步：
- 设计状态观测器估计系统"总扰动"（含故障、饱和和外部扰动）
- 基于反步法构建标称控制器
- 引入自适应滑模项在线补偿扰动估计误差
3. 核心算法实现
3.1 状态观测器设计
采用扩展状态观测器（ESO）结构：
matlab复制function [x_hat, z_hat] = ESO(y, u, params)
% y: 系统输出 (sigma, omega)
% u: 控制输入
% params: 观测器参数
persistent xi_hat;
if isempty(xi_hat)
xi_hat = zeros(6,1);
end
% 观测器动态
e = y - C * xi_hat;
dxi_hat = A * xi_hat + B * u + L * e;
xi_hat = xi_hat + params.Ts * dxi_hat;
% 总扰动估计
z_hat = xi_hat(7:end);
x_hat = xi_hat(1:6);
end
关键点在于观测器增益L的选择。论文采用带宽参数化法，将极点配置在-ω_o处（ω_o为观测器带宽）。实际调试中发现，ω_o需要比控制系统带宽大3-5倍，但过大会放大噪声。
3.2 反步控制器设计
反步控制分两步进行：
matlab复制function [u_nominal, alpha] = backstepping_control(sigma, omega, sigma_d, omega_d, domega_d, params)
% 第一步：虚拟控制量
z1 = sigma - sigma_d;
alpha = -params.c1 * z1 + 0.25 * ((1-sigma'*sigma)*eye(3) + 2*skew(sigma) + 2*(sigma*sigma'))' * omega_d;
% 第二步：实际控制量
z2 = omega - alpha;
u_nominal = J * (-params.c2 * z2 - z1 + dalpha) + skew(omega)*J*omega;
end
其中c1, c2为正常数，需要满足稳定性条件。实际实现时，dalpha的计算涉及雅可比矩阵，这是容易出错的环节：
注意：计算
dalpha时需要用到sigma和omega的导数，建议采用解析求导而非数值差分，后者会引入额外噪声。
3.3 自适应滑模补偿
设计滑模面：
code复制s = omega - alpha + k * (sigma - sigma_d)
自适应律采用σ修正策略防止参数漂移：
matlab复制function [u_comp, theta_hat] = adaptive_smc(s, params, theta_hat_prev)
% 自适应律
dtheta_hat = params.gamma * (norm(s) - params.kappa * theta_hat_prev);
theta_hat = theta_hat_prev + params.Ts * dtheta_hat;
% 补偿控制
u_comp = -theta_hat * sign(s);
end
实测中发现，直接用sign(s)会导致抖振。改进方案是用饱和函数sat(s/φ)代替，其中φ为边界层厚度。
4. Matlab实现关键技巧
4.1 执行器饱和建模
matlab复制function u_sat = actuator_saturation(u, limits)
% u: 指令力矩
% limits: [u_min, u_max] per axis
u_sat = min(max(u, limits(1)), limits(2));
% 记录饱和情况（用于故障检测）
persistent sat_count;
if isempty(sat_count)
sat_count = zeros(3,1);
end
for i = 1:3
if abs(u_sat(i) - u(i)) > 1e-6
sat_count(i) = sat_count(i) + 1;
end
end
end
4.2 故障注入模块
matlab复制function u_fault = inject_fault(u, t, fault_type)
switch fault_type
case 'partial_loss'
rho = diag([0.8, 0.5, 1.0]); % 执行器效率下降
phi = [0; 0; 0];
case 'bias'
rho = eye(3);
phi = [0.1; -0.2; 0.05] * (t >= 10);
case 'stuck'
rho = diag([1, 0, 1]); % 第二轴卡死
phi = [0; 0.5; 0] * (t >= 15);
end
u_fault = rho * u + phi;
end
4.3 稳定性分析代码
论文中的稳定性证明可通过数值验证：
matlab复制% 构造Lyapunov函数导数
V_dot = -z1'*c1*z1 - z2'*c2*z2 + z2'*(delta_hat - delta);
% 绘制V_dot随时间变化
plot(t, V_dot);
xlabel('Time (s)');
ylabel('$\dot{V}$', 'Interpreter', 'latex');
title('Lyapunov Function Derivative');
5. 仿真结果与参数整定
5.1 典型工况测试
设置三种测试场景：
- 仅饱和（各轴限幅±2Nm）
- 仅故障（t=10s时Y轴效率降为50%）
- 饱和+故障复合情况
性能指标对比：
| 场景 | 稳定时间(s) | 超调量(%) | 稳态误差(deg) |
|---|---|---|---|
| 仅饱和 | 8.2 | 4.5 | 0.03 |
| 仅故障 | 12.7 | 7.8 | 0.12 |
| 复合情况 | 15.3 | 9.2 | 0.18 |
5.2 参数整定经验
-
观测器带宽ω_o：从系统带宽的3倍开始，逐步增加直到估计误差收敛。过大会引入噪声，建议不超过采样频率的1/10。
-
滑模增益θ：初始值设为预期扰动上界的1.2倍，通过以下公式在线调整：
matlab复制if norm(s) > threshold theta = theta * 1.05; else theta = theta * 0.95; end -
边界层厚度φ：通常取控制精度的5-10倍。实测发现φ=0.05~0.1时能有效抑制抖振。
6. 工程实践中的问题与解决
6.1 计算时延问题
原始算法假设控制量可瞬时计算，实际数字实现存在一个采样周期的延迟。解决方案：
matlab复制% 预测下一步状态
x_pred = x + dx * Ts;
u = controller(x_pred); % 使用预测状态计算控制量
6.2 执行器分配问题
当使用推力器作为执行机构时，需要解决控制力矩到推力器开关指令的映射：
matlab复制function thruster_cmd = control_allocation(u, config)
% config: 推力器配置矩阵
A = config.geometry_matrix;
thruster_cmd = pinv(A) * u; % 伪逆求解
% 考虑最小点火时间约束
thruster_cmd(thruster_cmd < 0.01) = 0;
end
6.3 星敏感器噪声处理
实际角速度测量含噪声，需要在观测器前加低通滤波：
matlab复制% 二阶Butterworth滤波器
[b,a] = butter(2, 0.1); % 截止频率=0.1*Nyquist
omega_filt = filtfilt(b, a, omega_raw);
7. 代码架构建议
推荐采用面向对象设计，提高代码复用性：
matlab复制classdef FaultTolerantAttitudeControl < handle
properties
observer;
controller;
adapt_law;
fault_detect;
end
methods
function obj = FaultTolerantAttitudeControl(params)
obj.observer = ESO(params);
obj.controller = BacksteppingController(params);
obj.adapt_law = AdaptiveLaw(params);
end
function u = compute_control(obj, y, yd, t)
x_hat = obj.observer.estimate(y);
u_nominal = obj.controller.compute(x_hat, yd);
u_comp = obj.adapt_law.compute(x_hat);
u = u_nominal + u_comp;
end
end
end
完整实现代码已开源在GitHub（链接见文末），包含以下模块：
dynamics/：航天器动力学模型control/：控制器实现simulation/：测试场景脚本utils/：辅助函数（坐标转换、滤波等）