
简介这套MATLAB教学资源面向卫星姿态控制初学者以飞轮为执行机构实现三轴稳定控制核心采用经典PID算法动力学建模参考《航天器姿态动力学与控制》标准模型适合辅助课程设计或入门实践。包内提供两套完整实现路径纯M文件脚本便于理解状态方程推导与控制律编程Simulink配合S函数可开展图形化实时闭环仿真覆盖姿态角、角速度、飞轮转速与控制力矩等动态响应并能验证典型扰动下的系统稳定性。资源共14个文件以M文件、SLX模型、Python脚本和效果图为主压缩包仅594KB轻量易用模块命名清晰、变量注释规范便于调试、参数调整也为后续引入非线性补偿或切换至LQR等进阶策略打下基础。已有32人学习下载适合需要从原理到仿真快速上手的读者。 搞卫星姿控仿真的人应该都有过这种经历模型在M文件里跑得好好的换到Simulink重新搭一遍结果对不上或者反过来Simulink里看着直观想批量扫参数时效率又低得让人抓狂。我这个基于MATLAB的飞轮式卫星姿态控制仿真包就是冲着这个矛盾来的。它用同一套被控对象模型做了两条路线——纯M文件闭环仿真和SimulinkS函数模块化仿真两条路都能完整跑完一次从初始姿态到目标姿态的捕获过程结果可以互相校验。正在做姿控课程设计、毕业设计或者刚开始接触小卫星姿态控制仿真的人可以直接拿这份代码做底子按自己的需求改控制律、换执行机构参数。1. 飞轮姿态控制的底层逻辑角动量交换与模型选型1.1 反应轮是怎么把卫星“转”过来的飞轮不是火箭发动机那种“往外扔东西”的执行机构它靠电机驱动一个质量盘旋转在卫星本体和飞轮之间交换角动量。用大白话说卫星想往左转飞轮就往相反方向加速转把本体的角动量“吸”到飞轮盘上整个过程系统总角动量守恒。只要飞轮转速没碰到上限就能一直靠电力维持姿态。相比喷气执行机构飞轮不消耗推进剂适合长寿命任务相比磁力矩器它响应快、控制精度高所以绝大多数三轴稳定卫星都会用飞轮做主力执行机构。工程上最常见的布局是三只反作用飞轮正交安装分别对应卫星本体的三个轴。如果某个轮子转速达到饱和就需要靠磁力矩器或喷气卸载。这个仿真包里我做的是理想飞轮加一阶惯性环节飞轮转速上限和反作用力矩限幅都留了参数接口你可以直接改。执行机构力矩来源优势局限反应轮角动量交换精度高、响应快、不耗燃料会饱和、需要卸载手段喷气质量排出力矩大、无饱和问题消耗燃料、有羽流污染磁力矩器地磁场作用重量轻、可靠性高力矩小、受轨道位置限制1.2 仿真包的整体架构两条路线怎么分工设计这个仿真包的时候我给自己定的要求是M文件和Simulink模型必须共享同一套物理方程和同一套控制逻辑而不是各写各的。M文件路线用来做参数扫描、批处理、离线分析因为脚本天然适合循环跑上百组工况、自动生成图表Simulink路线用来做模块化展示方便后面接更复杂的执行机构模型、敏感器模型也方便在工程里跟其他子系统联调。整体架构分三层最底层是姿态动力学与运动学模型中间是飞轮执行机构模型上层是控制器。控制器输出力矩指令飞轮模型根据指令输出实际反作用力矩力矩反馈到动力学方程中姿态角速度变化又更新四元数形成一个闭环。这个闭环在两条路线里是一模一样的只是在代码组织方式上分道扬镳。这样做最大的好处是两边结果对不上时能迅速锁定问题出在数值求解还是模块连接上不用怀疑模型本身写错了。2. 纯M文件实现四元数、欧拉方程与控制律的一条龙打通2.1 状态向量怎么选四元数加角速度再加飞轮转速M文件实现的第一步是定义状态向量。我用的状态是x [q0; q1; q2; q3; ωx; ωy; ωz; Ω1; Ω2; Ω3]前四个是卫星本体的姿态四元数中间三个是本体系角速度最后三个是三只飞轮的转速。为什么不用欧拉角因为欧拉角在俯仰角接近±90°时会出现万向锁姿态机动仿真里经常要跨越大角度一锁就会出数值问题。四元数虽然多一个维度但没有奇异性而且四元数微分方程是线性的qdot 0.5 * [-q1 -q2 -q3; q0 -q3 q2; q3 q0 -q1; -q2 q1 q0] * [ωx; ωy; ωz]我在代码里直接把它写成矩阵乘法避免用四元数乘法符号视觉上更直观。每次积分之后还要把四元数归一化一次否则数值误差积累会让四元数模长偏离1最终导致姿态估计偏差。这个归一化步骤看起来不起眼但所有大角度机动仿真都离不开它。2.2 动力学方程与控制律怎么落到代码里刚体动力学用欧拉方程I * ωdot ω × (I * ω J_w * Ω_w) -J_w * Ωdot_w T_dist其中 I 是卫星主惯量矩阵J_w 是飞轮转动惯量Ω_w 是飞轮转速向量T_dist 是外部干扰力矩。等号左边是本体的角动量变化率右边是飞轮施加的反作用力矩和外部干扰。我在模型里把飞轮等效成一阶惯性环节τ * Ωdot_w Ω_w Ω_cmdΩ_cmd 是控制器输出的期望飞轮转速τ 是电机响应时间常数。代码实现是把这些方程写进一个ODE函数控制律放在同一个函数里function xdot sat_ode(t, x, param) q x(1:4); q q / norm(q); w x(5:7); rw_w x(8:10); % 控制律PD反馈 前馈解耦 T_cmd controller(q, w, param); % 飞轮一阶惯性模型 hw param.J_w * rw_w; w_dot_wheel (T_cmd - hw) / param.tau; % 刚体动力学 wdot param.I \ (-cross(w, param.I * w hw) T_cmd); % 四元数运动学 qdot 0.5 * omega_matrix(q) * w; xdot [qdot; wdot; w_dot_wheel]; end控制器我选的是经典PD加前馈解耦没有一上来就上滑模、自适应这些高级算法。原因很简单先确保基础物理是对的再谈优化。控制律表达式T_cmd -kp * qe_vec - kd * ωe ω × (I * ω J_w * Ω_w)其中 qe_vec 是误差四元数 qe 的矢量部分ωe 是角速度误差。前馈项 ω × (I * ω J_w * Ω_w) 用来抵消陀螺耦合力矩如果没有这一项卫星做大幅度姿态机动时会被耦合项带着跑偏PD反馈要花很大力气才能拉回来。误差四元数怎么求如果当前姿态四元数是 q目标四元数是 q_ref那么误差四元数qe q_ref^{-1} ⊗ q四元数求逆就是共轭加归一化代码里两行就能写完。注意四元数乘法的顺序不能反这在M文件里写起来就是一组嵌套的线性运算建议写成独立函数别图省事直接展开在主函数里不然读代码的人会被弄疯。2.3 用ode45跑闭环的固定步长细节M文件方案的核心是那个ODE函数主程序用ode45求解设置相对误差和绝对误差来控制精度。但这里有一个很容易踩的坑如果控制频率设计要求是固定步长比如100Hz而ode45是变步长就不能让控制器代码混在连续微分方程里否则控制器的采样节奏会被求解器打乱。我的做法是把控制律单独拆出来在主脚本里用固定步长循环调用。每个控制周期内先根据当前状态计算一次控制力矩再把飞轮指令转换成一阶惯性模型的输入然后用ode45推进一个控制周期内的小步长。这样模拟了星载计算机每个控制周期采样一次、计算一次、输出一次的节奏更接近真实飞行软件的行为。循环结构大概是for k 1:N t_now (k-1) * dt_ctrl; T_cmd controller(x_now, param); [t_out, x_out] ode45((t,x) sat_ode(t, x, T_cmd, param), ... [t_now, t_now dt_ctrl], x_now, opts); x_now x_out(end, :); record(k, :) [t_now dt_ctrl, x_now(1:7), param.J_w * x_now(8:10)]; end这样控制器是离散的被控对象是连续的两者之间的时序关系非常清楚。很多人直接在ODE函数里放控制器跑出来的效果也能看但仔细分析会发现控制力矩每个积分步都在变等效于一个极高频率的连续控制器这在实际工程里是做不到的。3. SimulinkS函数实现把同一套方程装进模块里3.1 Level-2 S函数的基本骨架Simulink路线里我坚持把被控对象模型写成一个S函数而不是用一堆积分器和乘法器搭。原因有三一是能和M文件共享方程代码不容易引入搭模块时的手误二是S函数有明确的连续状态定义积分交给求解器处理三是S函数是文本代码比几百根连线的子系统好做版本管理。Level-2 MATLAB S函数的核心是setup方法注册输入输出端口、连续状态数量和采样时间function att_sfun(block) setup(block); block.NumInputPorts 1; block.NumOutputPorts 1; block.NumContStates 10; block.NumSampleTimes 1; block.InputPort(1).Dimensions 3; block.InputPort(1).DirectFeedthrough true; block.OutputPort(1).Dimensions 7; block.RegBlockMethod(InitializeConditions, InitConditions); block.RegBlockMethod(Derivatives, Derivatives); block.RegBlockMethod(Outputs, Outputs); end连续状态导数写在Derivatives回调里输出写在Outputs回调里。我在仿真包里采用的是“被控对象在S函数里、控制器在外部子系统”的方案S函数输入是飞轮指令力矩输出是当前四元数和角速度。控制器子系统用什么算法、写成什么样都和前向通道隔离方便你替换成别的控制方法。3.2 四元数归一化和连续状态的那些坑S函数里连续状态是10维和M文件一样前四维是四元数。四元数归一化不能直接在Derivatives回调里做因为导数函数是给求解器用的如果改变了状态值会影响求解器的一致性判断。正确的做法是在Outputs回调里对输出的四元数做归一化或者把原始四元数作为状态所有用到姿态的地方都用归一化后的值。我实测过不归一化会怎样仿真开始阶段差别很小但跑了大角度机动之后四元数模长会漂移到1.002甚至更大控制误差随之变大最终姿态稳定精度变差。这类问题特别隐蔽因为它不会报错只会让你觉得“控制器参数哪里没调好”。我还踩过另一个更隐晦的坑S函数里连续状态的顺序必须和初始条件向量严格对齐一旦状态索引写错仿真结果会出现莫名其妙的发散而且发散点通常在好几秒之后排查起来非常费劲。3.3 求解器与采样时间设置的工程习惯Simulink环境里求解器我建议先用变步长、ode45把相对误差设到1e-6。等模型验证通过后再切到定步长模拟真实星载计算机的离散控制周期。这里有一个容易被忽略的点S函数如果不显式声明采样时间可能会被当作连续采样处理而控制器子系统如果用离散模块两者之间会出现离散/连续混合问题。我的处理是在控制器出口加一个零阶保持器让控制量在每个控制周期内保持不变同时把S函数的输入端口采样时间设成继承。这样整个环路的时序逻辑清晰控制器按固定周期更新被控对象在周期之间连续积分完全符合物理过程。算例验证下来这个结构和纯M文件的固定步长循环结果几乎完全一致差异只在求解器内部误差容限层面。4. 双实现结果的互验与调试心得4.1 同一组初值下两边的曲线对得上吗我拿同一组参数跑了两条路线初始姿态四元数 q0[0.8832; 0.3; -0.2; 0.3]目标姿态是单位四元数[1; 0; 0; 0]初始角速度[0; 0; 0] rad/s飞轮初始转速为零惯量矩阵取标准小卫星参数控制周期0.01秒。M文件跑完把姿态角和角速度曲线打出来Simulink模型用完全相同的控制器参数和初值条件跑一遍。两者在姿态机动到稳定的过程中曲线几乎完全重合稳态误差都在0.01度以内。有一两处中间时刻的瞬时差异主要来自ode45变步长和Simulink内部积分器的绝对误差容限策略不同属于正常数值差异不影响工程结论。如果你的两边结果差异很大不要第一时间怀疑控制律先检查初值尤其是Simulink里积分器初始条件是否设置成10维向量、S函数状态索引有没有错位。这类“两边都对却对不上”的问题九成是把四元数顺序填错了。4.2 我踩过的几个坑和排查方法第一个大坑是飞轮饱和没建模。刚开始为了省事飞轮模型是理想积分器仿真5秒后转速持续上升控制器还在不断输出结果姿态虽然稳了飞轮转速却高得离谱。加了角动量上限之后饱和工况下控制精度立刻恶化你才真正意识到卫星为什么需要卸载手段。仿真模型加不加饱和结论完全两样这个一定要在模型里体现。第二个坑是控制器增益的量纲。姿态控制里kp和kd带的是不同量纲很多人直接把kp取成10、kd取成20仿出来振荡得厉害。我习惯先把问题归一化用二阶系统自然频率和阻尼比来换算kp I * ω_n²kd 2 * I * ζ * ω_n。比如取ω_n0.1 rad/s、ζ0.9反推出kp和kd调参就有依据而不是瞎试。三轴惯量不同时按每个主轴惯量分别算。第三个坑是Simulink里信号维度错误。S函数输入如果是力矩指令前面控制器输出维度如果设置不对Simulink经常会报维度不匹配。排查方法是在S函数setup里显式声明端口维度block.InputPort(1).Dimensions 3;这行代码能省掉大量让人抓狂的连线错误。最后再分享一个小技巧双实现模型里我特意把参数集中放在一个MATLAB脚本里M文件和Simulink模型都从这个脚本加载参数改一个参数、两条路线同步生效。这样调参时就不会面临“Simulink里改了、M文件忘了改”的尴尬。做姿控仿真最怕的就是模型可信度打折扣两条路线互相校验至少能帮你把大部分低级错误挡在门外。本文还有配套的精品资源点击获取