4. 仿真引擎与数值积分方法
航天器GNC数字仿真平台的核心是仿真引擎，而仿真引擎的“心脏”是常微分方程（ODE）求解器。GNC系统涉及轨道动力学、姿态动力学、控制律、执行机构与敏感器模型，这些模型最终都被写成形如 ẋ = f(t, x, u) 的一阶微分方程组。求解器负责在时间轴上推进状态量 x，其选型、步长策略与误差控制直接决定仿真结果的精度、稳定性和计算效率。本节围绕求解器选型、定步长与变步长积分、刚性系统处理、积分误差与稳定性分析四个要点展开。
一、求解器选型：从问题特性出发
ODE求解器大致分为显式与隐式两大类。显式方法（如RK4、Dormand-Prince RK45、Bogacki-Shampine RK23）每一步只需计算若干次右端函数 f，实现简单、单步代价低，适合非刚性、光滑的动力学问题，例如轨道递推、姿态运动学积分。隐式方法（如后向欧拉、梯形法、BDF、Radau IIA）每一步需要求解代数方程组，单步代价高，但具有更好的稳定性，适合刚性系统。
选型时先判断问题是否刚性。一个实用的判据是估计系统雅可比矩阵特征值的实部量级：若最大与最小特征值实部之比（刚性比）远大于1，且最小特征值对应快速衰减模态，则问题偏刚性。GNC仿真中，轨道动力学通常非刚性；但包含柔性太阳翼、液体晃动、高增益姿态控制回路或精细执行机构模型时，系统可能出现从毫秒到小时量级的多时间尺度，刚性比可达10³～10⁶，此时应优先选择隐式或半隐式求解器。
常用选型建议如下：轨道长期递推可选RK45或DOP853；姿态快速机动仿真可选RK4或RK45；含柔性附件与高带宽控制回路时选Radau IIA或BDF；实时仿真对单步时间有硬约束时，可选固定步长RK4配合足够小的步长，或采用半隐式方法。
二、定步长与变步长积分
定步长积分以固定时间间隔 h 推进，实现简单、时间确定，适合实时仿真和需要与固定采样周期对齐的GNC闭环测试。其缺点是：当动态变化剧烈时，固定 h 可能过大导致误差超限；当动态平缓时，h 又过小造成计算浪费。定步长RK4的局部截断误差为 O(h⁵)，全局误差为 O(h⁴)，步长减半误差约降为原来的1/16。
变步长积分根据局部误差估计自动调整步长。以嵌入式RK45为例，每一步同时计算4阶和5阶解，二者之差作为局部误差估计 err。若 err 小于容差 tol，则接受该步并按公式调整下一步步长：
h_new = h * min(max_factor, max(min_factor,
0.9 * (tol / err)^(1/5)))
其中 max_factor 常取2～5，min_factor 常取0.1～0.5。若 err 超过 tol，则拒绝该步，缩小步长重算。变步长方法在平缓段自动放大步长，在剧烈段自动缩小步长，兼顾精度与效率，是通用仿真的首选。但其步长不固定，若GNC控制器按固定周期采样，需要在求解器输出与控制器采样之间做插值或事件处理。
三、刚性系统处理
刚性系统的困难在于：显式方法的稳定域有限，为保持稳定，步长必须小于最快模态的时间常数，即使该模态早已衰减完毕，仍被迫使用极小步长，导致计算量剧增。例如某姿态控制回路含时间常数为1 ms的快速模态，若用显式RK4，稳定步长约需小于0.3 ms；而仿真总时长可能达数千秒，计算量不可接受。
处理刚性系统有三条主要途径。第一，改用隐式方法，如BDF或Radau IIA，其稳定域覆盖左半平面大部分区域，允许步长由精度而非稳定性决定。第二，对刚性部分做隐式处理、非刚性部分做显式处理，即IMEX方法，兼顾效率与稳定。第三，对系统进行模型降阶或模态截断，剔除对GNC性能影响很小的快速模态，从源头降低刚性。实际工程中，若柔性附件模态频率远高于控制带宽，常保留前若干阶模态并加阻尼，其余高频模态用准静态近似或直接截断。
四、积分误差与稳定性分析
积分误差分为局部截断误差和全局累积误差。局部截断误差是单步内用数值解代替精确解产生的误差；全局误差是各步误差累积的结果。误差控制通常通过相对容差 reltol 和绝对容差 abstol 实现，接受判据为 err ≤ max(reltol·|x|, abstol)。对于状态量量级差异大的系统（如位置以米计、姿态角以弧度计），应分别设置绝对容差，避免小量被大量“淹没”。
稳定性分析常用线性测试方程 ẋ = λx，λ 为复数。显式RK4的稳定域在复平面上是一个有限区域，要求 hλ 落在该区域内。若 λ 实部为负且模很大，h 必须很小。隐式后向欧拉对任意 Re(λ)<0 都稳定，称为A稳定；Radau IIA还具备L稳定，能快速抑制刚性模态。实际仿真中，可通过监测步长变化、误差估计和能量/动量守恒情况判断积分是否健康：若步长持续被压到极小、误差估计频繁超限，往往说明系统刚性或模型存在不连续，需要检查模型或更换求解器。
五、工程实践中的参数配置示例
以下给出一个GNC仿真中常见的求解器配置示例，用于姿态轨道联合仿真，含中等刚度柔性附件模型。参数为教学示例，实际使用需根据模型量级调整。
| 配置项 | 推荐值 | 说明 |
|---|---|---|
| 求解器 | Radau IIA（5阶） | 刚性系统，L稳定 |
| reltol | 1e-8 | 相对容差 |
| abstol（位置） | 1e-3 m | 轨道位置绝对容差 |
| abstol（姿态角） | 1e-8 rad | 姿态角绝对容差 |
| abstol（角速度） | 1e-10 rad/s | 角速度绝对容差 |
| 初始步长 | 1e-3 s | 由求解器自动调整 |
| 最大步长 | 0.1 s | 限制单步跨度 |
| 最小步长 | 1e-8 s | 防止步长过小死循环 |
若改为实时仿真且控制器采样周期为10 ms，可切换为定步长RK4，步长取1 ms，保证每10个积分步对应一个控制周期。此时需检查最快模态时间常数是否大于约3倍步长，否则应改用半隐式方法或对快速模态做准静态处理。
综上，求解器选型不是“越高级越好”，而是与问题刚性、精度要求、实时约束和模型特性匹配。显式方法适合非刚性、光滑问题；隐式方法适合刚性、多时间尺度问题；定步长适合实时与固定采样场景；变步长适合通用高精度仿真。误差容差应按状态分量分别设置，稳定性分析应结合测试方程和实际步长行为判断。只有把求解器、步长策略与GNC模型特性统一考虑，仿真结果才具备可信度。