简介：面向航天器姿态控制设计者与相关专业学习者，这份压缩包提供基于四元数的自适应滑模控制完整实现方案。其中包含Simulink仿真模型（mdl文件）、MATLAB脚本（m文件）以及仿真结果图（fig文件），共7个文件，并附有r2012a与r2014a两个版本兼容说明，包体仅37KB，便于快速部署与验证。自适应滑模控制器结合四元数姿态表示，可有效应对参数不确定性与外部扰动，读者可依此学习滑动变量设计、参数更新法则及稳定性分析思路，并通过示例模型直接运行仿真。目前已有476人学习下载，适合正在研究航天器鲁棒控制或滑模控制算法的工程师与研究生参考。

1. 航天器姿态控制：为什么四元数搭配自适应滑模是这组代码的正确用法

航天器姿态控制里有个反直觉的经验：用欧拉角做状态量做仿真，模型跑得好好的，一到真实任务就翻车。原因不玄学，就是万向锁和转动顺序导致的数值病态。而用四元数描述姿态，可以绕开这个问题，代价是控制器的非线性项更绕。这个压缩包里的东西，就是一套已经按四元数形式写好的自适应滑模控制仿真，核心文件不复杂——plant.m、contr.m、SMC.mdl、两张.fig图形，另带r2012a/r2014a两个版本模型。适合正在做飞行器或航天器姿态控制课题的人拆开改。资源解决的具体问题是：在转动惯量未知、外界力矩干扰不确定的情况下，怎么用一套能实时调整增益的滑模控制器，把航天器姿态稳定跟踪到目标四元数。下面按我的拆包习惯，从模型到参数到踩坑，一层层讲。

2. 先立住理论：四元数姿态误差怎么进状态方程

2.1 四元数表示姿态：一句话说清为什么不用欧拉角

四元数本质是一个实部加三个虚部，通常写作 q = [q0, q1, q2, q3]ᵀ。它的物理含义是：绕着空间中某个单位轴 e = [e1, e2, e3]ᵀ 旋转角度 θ，等价于四元数分量 q0 = cos(θ/2)，q1 = e1·sin(θ/2)，q2 = e2·sin(θ/2)，q3 = e3·sin(θ/2)。所以四元数天然满足 q0² + q1² + q2² + q3² = 1 这个约束，这个约束在仿真里必须强制维持，否则积分会飘。

欧拉角的问题是转动顺序相关，滚动、俯仰、偏航三个角度按不同顺序转出来的姿态不一样，而且当某个角度接近 ±90° 时会出现万向锁。四元数不存在这两个问题，多一个维度换来了数值上的稳定。代价是控制器的设计变量不再是三个角度，而是四个分量加一个单位约束，滑模面、误差定义都要重新写。航天器姿态动力学里，角速度还耦合着陀螺效应，所以状态方程必然是非线性的。

2.2 姿态动力学与误差四元数方程

航天器姿态动力学最基本的形式是 J·ω_dot + ω×(J·ω) = u + d，其中 J 是转动惯量矩阵，ω 是本体角速度，u 是控制力矩，d 是外部干扰力矩。这个方程很多教材上都一样，但实际建模时的门槛在于：J 拿不准，d 上界也拿不准。卫星铰链展开后质心偏移、燃料消耗造成质量分布变化，转动惯量能偏 10% 到 20%。这就是为什么纯滑模控制器在仿真里表现好、上了真实系统就抖或者静差大的原因——滑模的鲁棒性依赖干扰上界已知，而上界本身是估计值。

自适应滑模的做法是把「上界未知」这个矛盾转移到自适应律上：先假设干扰上界是一个未知常数 ρ，设计一个估计值 ρ̂ 在线更新，然后控制器里用 ρ̂ 来代替 ρ。压缩包里的 plant.m 大概率就是按这个路子写航天器动力学模型。我一般会把它拆成两个文件理解：一个负责四元数运动学方程，一个负责角速度动力学方程。

姿态误差四元数是控制器的输入之一。定义期望姿态 q_d，实际姿态 q，误差四元数 q_e = q_d⁻¹ ⊗ q，其中 ⊗ 是四元数乘法，q_d⁻¹ 是共轭四元数。误差四元数的前三维 q_e(1:3) 接近零时，说明实际姿态追上了期望。这里有个关键参数叫矢量部分增益 k，它决定姿态误差在控制器里放多大权重，调小了下垂，调大了超调。

实际写plant.m时，状态量取 x = [q0 q1 q2 q3 ωx ωy ωz]ᵀ，前四个是四元数分量，后三个是本体系角速度。运动学方程是 q_dot = 0.5 * q ⊗ [0 ω]ᵀ，动力学方程是 ω_dot = J⁻¹ * (u + d - ω×(J·ω))。下面是一个参考实现的外层骨架：

function [sys, x0, str, ts] = plant(t, x, u, flag, param) % 基于四元数的航天器姿态动力学模型（S-Function格式） % x = [q0; q1; q2; q3; wx; wy; wz] % u = [Mx; My; Mz] 控制力矩，param包含惯量等参数 switch flag case 0 % 初始化 sizes = simsizes; sizes.NumContStates = 7; sizes.NumDiscStates = 0; sizes.NumOutputs = 7; sizes.NumInputs = 3; sizes.DirFeedthrough = 0; sys = simsizes(sizes); x0 = [1; 0; 0; 0; 0; 0; 0]; % 初始姿态为单位四元数 str = []; ts = [0 0]; % 连续系统 case 1 % 状态导数 J = param.J; % 惯量矩阵，单位 kg·m^2 d = param.disturb(t); % 外部干扰力矩，时变 q = x(1:4); w = x(5:7); q = q / norm(q); % 强制归一化 % 四元数运动学：qdot = 0.5*Omega(w)*q Omega = [0, -w(1), -w(2), -w(3); w(1), 0, -w(3), w(2); w(2), w(3), 0, -w(1); w(3), -w(2), w(1), 0]; qdot = 0.5 * Omega * q; % 角速度动力学：J*wdot + w x (J*w) = u + d cross_val = cross(w, J*w); wdot = J \ (u(:) + d(:) - cross_val); sys = [qdot; wdot]; case 3 % 系统输出 sys = x; case 9 % 终止 sys = []; otherwise sys = []; end

这段代码里有两个关键点是新手容易看漏的。第一，q = q / norm(q)强制归一化放在每个积分步里，防止四元数因数值积分误差脱离单位球面。第二，cross(w, J*w)是陀螺力矩项，它让系统呈现强耦合，控制器必须考虑耦合而不是简单解耦成三个独立的通道。参数param.disturb(t)是外部干扰的时变函数，在自适应滑模里通常写成有界形式，比如常值加正弦的组合，用来验证控制器补偿干扰的能力。

3. 滑模控制器与自适应律：从原理到 contr.m 落地

3.1 滑模面设计：误差四元数怎么凑出一个二阶滑模面

滑模面设计是整个控制器的核心。对航天器姿态控制，常用的非线性滑模面是 s = ω_e + λ·q_e(1:3)，其中 ω_e 是角速度误差，λ 是正定的增益矩阵。这个面的物理含义是：当 s = 0 时，角速度误差与姿态误差成比例收敛，系统沿着这个面滑动，最终姿态误差和角速度误差同时归零。这里 λ 的取值决定了收敛速度，一般取对角线矩阵，比如 diag(0.5, 0.8, 0.6)，单位是 s⁻¹。

为什么不用线性滑模面 s = ω_e + λ·∫q_e(1:3)dt？因为四元数误差本身已经是一阶积分性质的量，再叠加积分项会让系统阶数变高，自适应律也更难稳定。用 s = ω_e + λ·q_e(1:3)，滑模动态几乎是一阶的，Lyapunov 分析证起来顺，控制律也简洁。

滑模到达条件一般选 s·ṡ ≤ -η·|s|，其中 η 是到达速度系数。这个条件保证系统状态不管从哪里出发，都能在有限时间到达滑模面。到达后系统对匹配型干扰完全不敏感，这就是滑模控制鲁棒性的来源。注意「匹配型」三个字，如果干扰不在控制通道上，滑模就无能为力了，航天器里的干扰力矩都直接作用在本体系上，所以属于匹配型。

3.2 自适应控制律：不确定上界的在线估计与更新法则

控制器完整形式分三部分：等效控制项、切换项、自适应补偿项。等效控制项让系统在标称模型下保持在滑模面上，表达式是 u_eq = ω×(Ĵ·ω) + Ĵ·(q_e_dot - λ·ω_e)，其中 Ĵ 是标称惯量矩阵。切换项是 -K·sgn(s)，K 必须大于干扰上界才能保证到达条件成立。自适应补偿项用估计上界 ρ̂ 来调整切换增益，让 K 不是一个定值，而是跟随 ρ̂ 增长或收缩。

自适应律的更新公式是 ρ̂_dot = Γ·|s|，其中 Γ 是自适应增益，ρ̂ 是干扰上界的估计值。这里有个工程直觉：干扰越大，|s| 越大，ρ̂ 增长越快，切换增益跟着变大，就能覆盖掉实际干扰。但 Γ 不能取得太大，否则 ρ̂ 会震荡甚至发散，这一点我在第 5 章会单独展开。

contr.m 里一般的实现顺序是：输入误差四元数和角速度误差 → 计算滑模面 s → 计算等效控制 → 计算自适应更新量 → 叠加切换项 → 输出控制力矩。下面这段代码是自适应律和控制律拼在一起的参考写法：

function u = contr_MIMO(t, x, qe, we, param) % 自适应滑模控制器：输入误差四元数、角速度误差，输出控制力矩 % qe = [qe0; qe1; qe2; qe3]，we = [wex; wey; wez] lambda = param.lambda; % 滑模面增益，diag矩阵 Gamma = param.Gamma; % 自适应增益，决定上界估计速度 K0 = param.K0; % 切换增益初值，兜底用 s = we + lambda * qe(1:3); % 滑模面，核心 % 等效控制：基于标称惯量Jhat，补偿陀螺项 Jhat = param.Jhat; Omega_w = [0, -we(3), we(2); we(3), 0, -we(1); -we(2), we(1), 0]; % 等效控制 = Jhat*(lambda*w_e_dot - Jhat⁻¹*(wx J w)) 近似为： u_eq = cross(we, Jhat*we) + Jhat * (lambda * (param.w_dot_ref - we)); % 参考旋转忽略 % 自适应估计：干扰上界rho_hat的微分 rho_hat = param.rho_hat; % 状态变量 rho_hat_dot = Gamma * norm(s); % 关键自适应律 % 切换项带边界层，用饱和函数避免抖振 phi = param.phi; % 边界层厚度，0.05典型值 u_sw = -K0 * s / (norm(s) + phi) - rho_hat * s / (norm(s) + phi); % 或者直接用tanh做平滑切换 param.rho_hat = rho_hat + rho_hat_dot * param.dt; % 在循环外更新 u = u_eq + u_sw; end

注意这段代码里我用了s / (norm(s) + phi)而不是sign(s)，这是工程上必须的改动：理想滑模切换频率无穷大，仿真里会造成抖振，MATLAB 求解器根本没法处理高频切换，会报「Matrix is singular or badly scaled」。边界层 phi 越小，滑模越接近理想切换，控制精度越高，但抖振越强；phi 越大，抖动越小，但稳态误差会被放大。实际调参时，先给 phi = 0.1 观察抖振，再逐步减小到 0.01，看力矩曲线是否出现高频毛刺。这个 trade-off 是滑模控制最真实的工程手感所在。

3.3 参数没调对会怎样：Γ、K0、λ 的作用边界

参数太多是控制类代码包里最常见的现象。滑模控制器在姿态问题里真正要调的只有三个数：λ、Γ、K0，以及边界层厚度 phi。λ 决定滑模面的斜率，取大会让姿态响应变快，但也放大噪声。Γ 决定自适应估计的速度，取太大导致增益震荡，通常从 Γ = 0.01 起步，跑 20 秒仿真看 rho_hat 曲线是否单调收敛再决定要不要加大。K0 是切换增益的兜底值，当自适应估计还没跟上时保障基本的鲁棒性，不能取零，我一般取 diag(0.5, 0.5, 0.5)。真正的干扰上界不需要知道精确值，这就是自适应滑模相对普通滑模的最大便利——普通滑模的切换增益直接依赖干扰上界，取小了不稳定，取大了抖振，自适应滑模让系统自己找这个界。

4. 拆包与仿真落地：把六个文件串成一条能跑的回路

4.1 压缩包里每个文件干什么：按我的读码顺序拆

拿到压缩包别急着开 Simulink。先按这个顺序读文件，能少踩一半的坑。文件清单如下：

读码顺序建议是：先读 plant.m 搞清楚状态方程里含哪些项，再看 contr.m 看控制器假设的干扰条件，最后开 SMC.mdl 对照回路核对信号维度。很多新手先双击 mdl 报错，然后怀疑资源有问题，其实不过是版本兼容问题。那两个 .fig 是仿真后保存的图形句柄，双击会打开一个空 Figure，需要自己重新绘图或者从 Workspace 里调变量画图。

4.2 把控制律灌进 Simulink：信号维度与循环闭合检查

SMC.mdl 的内部结构一般是五个模块串联：期望姿态输入 → 误差四元数计算 → 控制器（嵌入 contr.m 的 MATLAB Function 块）→ 力矩限幅/分配 → plant.m S-Function → 状态反馈回期望端。其中误差四元数计算这一步最容易出错，因为四元数乘法不是按元素相减，而是 q_e = q_d⁻¹ ⊗ q。很多人在 Simulink 里直接拿 q_d - q 做差，结果系统看起来能收敛，但姿态轨迹完全不对。正确做法是单独写一个 MATLAB Function 块做四元数乘法：

function qe = quat_error(qd, q) % 误差四元数：qd的共轭乘以当前四元数 qd_inv = [qd(1); -qd(2); -qd(3); -qd(4)]; qe = quat_multiply(qd_inv, q); % 标准四元数乘法 end

在 Simulink 里连接 S-Function 时还有一个隐藏要求：plant.m 的输出维度与输入维度必须和 S-Function 配置一致。如果 plant.m 的输出是 [7x1] 而误差计算模块的输出是 [4x1]，直接连会报维度不匹配。我一般会在 plant.m 的 S-Function 面板里手动设置输出维度为 7，把仿真停在一个 fixed-step 求解器上，比如 ode4（四阶龙格库塔），步长固定 0.001 或 0.01 秒。自适应滑模控制里的不连续切换项，如果用变步长求解器，误差容差控制会让求解器在切换瞬间无限加密计算步长，仿真时间拉长几十倍。固定步长 + 边界层处理是防止仿真卡死的常规做法。

4.3 从模型文件恢复旧版本：r2012a 与 r2014a 双版本怎么选

压缩包里有两个 mdl 副本后缀，这是作者在不同 MATLAB 版本下保存的格式标记。Simulink 模型的 .mdl 文件本质是纯文本格式，但不同版本用了不同的 schema 标准。新版 MATLAB（r2015b+ 以后）打开 r2012a 格式的 mdl 通常能自动升级，但反过来不行。如果你本机装的是 r2016a 以下版本，打开 SMC.mdl 应该没问题；如果装的是 r2018 以上，打开 r2014a 版本更稳妥。遇到加载失败报 「Error loading model file」 或 「Model contains blocks which require a newer version」 时，不要删文件，直接在模型文件第一行把版本号字段改掉有时能强制加载，但我不推荐——模块库的结构变化会导致部分模块丢失。最常见的情况是找不到 S-Function 的对应端口类型，这种问题直接换成 slx 格式重新保存即可。旧的 .fig 文件同理，在 r2012a 之前的 Figure 保存格式与新版兼容性差，双击打不开就用hgload('untitled.fig')手动加载。

5. 避坑与常见问题：跑这套仿真时最容易翻车的五个点

5.1 四元数不归一化：仿真越跑越飘

现象：系统稳定之后姿态误差不收敛，反而有缓慢漂移，四元数分量 q0 慢慢从 1 变成 0.98、0.95，最终发散。原因：Simulink 积分器只对状态分量做积分，四元数的单位模约束不在积分器的约束范围内，每个积分步都会引入数值误差，误差累积后四元数偏离单位球面。解决：在 plant.m 的状态导数段加入强制归一化，就是 2.2 节代码里q = q / norm(q)那行。另一种方案是把四元数归一化放到状态输出的 feedback 路径上，但这样进入控制器的四元数虽然是单位模，被积分的原始四元数仍然逸散，不推荐。强制归一化的频率越高，数值越稳，代价是引入非线性，对固定步长求解器没有影响。

5.2 抖振导致力矩信号在切换阈值附近高频翻转

现象：控制力矩输出曲线像锯齿一样，在某个值上下快速振荡，仿真步长被迫缩小。原因：滑模切换项使用了理想 sign(s) 函数，切换点附近的符号变化频率超出求解器处理能力。解决：把 sign(s) 替换为s/(norm(s)+phi)饱和函数，或使用tanh(s/phi)平滑近似。phi 的取值直接决定抖振幅值，从 0.05 开始试，观察姿态角响应曲线是否变顺滑。还有一个次要原因是自适应增益 Γ 过大，导致 rho_hat 高频波动带动切换项产生高频分量。逐个排除的顺序是：先改 phi，再降 Gamma，最后调 K0。

5.3 自适应增益 Gamma 调大后系统直接发散

现象：仿真前几秒还正常，几秒后角速度快速增长到 1e6 量级，矩阵计算报奇异性。原因：自适应律 rho_hat_dot = Gamma*|s| 是无记忆的纯积分，Gamma 越大，rho_hat 更新越猛，切换增益跟着暴涨，把控制力矩推到远超执行机构实际能力的范围。这属于增益过大导致的相互激励：控制力矩变大→姿态变化更快→误差变大→|s|变大→rho_hat 继续涨。解决：给自适应律加死区，比如只有 |s| 大于某个阈值时才更新 rho_hat，或者用 σ-修正法，在更新律里加一项 -σ·rho_hat 防止无限增长。在 contr.m 里实现死区也就几行代码：

if norm(s) > param.s_deadzone rho_hat_dot = Gamma * norm(s); else rho_hat_dot = 0; % 误差太小不更新，防止增益膨胀 end

这个死区阈值 s_deadzone 一般取 0.01 到 0.05，太小起不到作用，太大静差明显。

5.4 Simulink 旧版本模型在新版里打不开

现象：双击 SMC.mdl.r2012a，MATLAB 报错说模型格式来自未来版本或未知版本。原因：mdl 文件头部的版本标识与当前 MATLAB 的 Simulink 版本不兼容。解决：先尝试打开 SMC.mdl 主文件而非带版本后缀的文件；如果主文件也打不开，用文本编辑器打开 mdl 文件头部，把Simulink.Version字段改为当前版本号，绝大多数情况下能打开，但如果用了旧版本不存在的模块库仍然会报缺失。稳妥的路径是使用版本对称的 MATLAB 环境，r2014a 副本在 r2014a 到 r2016a 内打开最顺畅。新版打开后建议直接另存为 slx 格式，后续修改不会再触发兼容问题。

5.5 用欧拉角观察姿态结果：从图上看收敛了，实际问题一大堆

现象：把仿真输出里的欧拉角转出来看，滚转角和俯仰角都收敛到 0 附近，于是判定控制器设计成功。原因：欧拉角在某些角度区间会跳变，比如俯仰角接近 ±90° 时，滚转角和偏航角定义退化，导致图上曲线出现 180° 的突变或锯齿。解决：判断姿态误差是否收敛，直接用四元数本身的误差范数 |q_e(1:3)| 或者演化成轴角误差看大小，而不是转换成欧拉角观察。另一个陷阱是误差四元数的符号二义性，q 和 -q 表示同一个姿态，所以 4.2 节里的误差四元数计算必须固定约定，一般取 qe(1) ≥ 0，避免控制器看到误差分量跳变。

6. 进阶技巧：给这套仿真加一脚稳健性验证

调试到能收敛只是第一步，离真正交付还差一口气。我拿到这套代码后的习惯做法是，先把控制器参数固定，然后用三组额外的仿真来验证鲁棒性和边界行为。

第一组是转动惯量偏差测试。把 plant.m 里的 J 与 contr.m 里的 Jhat 设成不同的值，偏差幅度从 ±5%、±10%、±20% 逐级增加。记录每个偏差下的姿态误差稳态值和调节时间。如果 J 偏差 20% 时误差还能维持在 0.01 数量级，说明自适应项起效了；如果误差已经超过 0.5，说明等效控制依赖模型的部分占比太大，需要增大 K0 或者放宽滑模面的 λ。这一步能把资源从「能跑」推进到「知道自己能抗住什么程度的不确定性」。

第二组是干扰频谱测试。把 param.disturb(t) 从简单的常值干扰换成一组叠加扰动：常值项 + 周期项 + 随机偏置项，周期项的角频率从低到高拉几个量级。自适应滑模对低频正弦的补偿效果好，但对高频扰动的补偿能力会下降，因为 rho_hat 的更新速度有限。这个测试暴露的通常是参数 Γ 的带宽受限问题，如果高频干扰下力矩曲线已经明显滞后，就要把 Γ 调大，或者增加一个前馈补偿项，而不是一味加大切换增益。

第三组是符号边界验证。把初始姿态误差设置成接近 180°，四元数分量接近 [-1, 0, 0, 0]，这个点位是四元数误差计算的数值敏感区。很多控制器在这个姿态附近会出现一瞬间的大信号脉冲，因为误差四元数的矢量部分在 180° 附近对扰动特别敏感。如果仿真里这个点位的控制力矩超限了，处理办法是给控制力矩做限幅，或者在误差四元数进入死区时用角速度项主导收敛，不让姿态误差信号单独驱动过大输出。

这三组验证跑下来，代码包本身的边界基本就摸清了。从那以后我每次拿到别人的控制仿真代码，都强制走一遍同样的流程：先归一化四元数，再看切换项用没用饱和函数，接着做模型偏差扫描，最后才谈参数调优。这套流程过滤掉了大部分看上去漂亮、实际跑不了实物的验证模型。希望帮到你。

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