I’m working on STM-based covariance propagation as part of a personal study project, using Orekit as a reference. I understand that Orekit computes the Jacobian matrix (STM) using Hipparchus’ automatic differentiation.
I’ve implemented an autodiff-based Jacobian in my own code and, when tested independently (i.e., outside of propagation), the Jacobian output closely matches the one from Hipparchus. However, when I integrate this Jacobian computation into my full propagator and compare the resulting STM and propagated covariance with Orekit’s, I observe significant discrepancies. In particular, my STM shows oscillatory patterns over a 1-day propagation window, whereas Orekit’s output remains stable and linearly increasing.
I’ve reviewed the documentation and relevant literature, but haven’t been able to determine what’s causing this divergence. I’d appreciate any suggestions or insights on what might be leading to this behavior, especially during the propagation phase.
Thanks in advance !
Which propagator and force models are you using?
If applicable, which numerical integrator are you using (both from Hipparchus for the Orekit reference and your own implementation)?
I haven’t used the STM methods in a long time now, but from what I remember the STM is computed in the same coordinate type used in the integration, which by default is Equinoctial (when building the propagator with NumericalPropagator, if you use the NumericalPropagatorBuilder I think it uses the coordinate type of the orbit passed as input).
Are you sure you’re comparing the matrices in the same coordinates type?
Also, if I’m not mistaken the STM is not computed with automatic differentiation, but rather from the variational equations, but again, I’m not sure.
Yes, I’m aware that Orekit’s STM is integrated via the variational equations, and those are built using AD to get the partial derivatives of acceleration.
In my implementation, I’m also computing the variational equations using AD to generate the Jacobian of acceleration w.r.t. state. I then integrate both the state and STM forward.
Emiliano, good point about the coordinate system. I’m using Cartesian coordinates (J2000) for both integration and STM. I’ll double-check that Orekit’s setup is aligned with mine.
What still puzzles me is the oscillatory behavior in my STM entries over a 1-day propagation, while Orekit’s STM grows more smoothly (almost linearly in some cases). Has anyone seen this sort of thing before?
Happy to share a test case or plots if that helps.
If I have understood well, if we denote with \dot{\boldsymbol{x}} = \boldsymbol{f}(\boldsymbol{x}(t), t) the ODE that describes the dynamics, and with \boldsymbol{\Phi}(t, t_0) = \partial \boldsymbol{x}(t)/\partial \boldsymbol{x}(t_0) the STM, then the STM satisfies
Yes, with the subtility that AD is applied to the Cartesian coordinates and that \Phi is computed with them. A conversion is applied if applicable when querying for the STM with the MatricesHarvester.
Ah interesting, I thought the variational equations were given with the analytical formulation, not via AD. I guess AD makes more sense if you need to include drag, since I don’t think there’s an analytical form for that one.
Thanks for the correction!
Actually this calls for another clarification. With drag, Orekit has two ways to compute the derivatives of the atmospheric density w.r.t. the Cartesian vector (the default one being the first):
Finite differences
AD
This is because the Field version of some models e.g. NRLMSISE00 can be computationnally intensive. In fact it can easily become the bottleneck of propagation.
Last year I proposed a potential improvement but it received mixed opinions so I’ve not implemented it.
Out of curiosity @Serrof , do those bottlenecks come into play only when the NRLMSISE derivatives are used, or would this issue also occur if you used a computation method that did not require derivative calculations, such as unscented propagation?
@Ningguang you might already have picked up on this, but just in case you have not, drag and solar are not included in default covariance propagation. You can manually propagate the covariance by propagating the STM including the drag and solar additional parameters, and then using the STM to manually propagate the covariance. This might or might not be contributing to your issues.
What do you mean by this? Aren’t the 6x6 STM and 6x6 state covariance computed using the same force model setup for the propagator to which the StateCovarianceMatrixProvider is attached to?
I would think that manual setup is only needed if you want to compute the partials of these parameters, e.g. to estimate them in the OD process.
Do the nominal trajectories computed with your implementation and Orekit always agree? Then I would start debugging by evaluating the partials \partial\boldsymbol{f}(\boldsymbol{x}(t),t)/\partial\boldsymbol{x}(t) along the same reference, and compare the outputs between the two implementations.
Even the non-Field version of NRLMSISE00 is computationally expansive, which is why even with finite differences, the STM is still costly. Yeah some estimators like the UKF are derivative free. But propagating N sigma points can still be pretty slow. Especially if measurements are close to each other in time, then starting and stopping propagators a lot is not very efficient.
@afvy Yes you’re right - what I said is only relevant if you want the full covariance; it doesn’t matter for the STM. I am just so used to assuming that if you want the STM that it’s for the covariance that my brain jumped ahead I think.
Thanks all — this has been incredibly informative!
As a follow-up, I’ve run a focused test case with only the J2 perturbation included — no drag, solar pressure, or other force models. The goal is to isolate where the divergence originates.
Here’s a comparison of the propagated covariance matrices from Orekit and my own propagator, using the same initial conditions and force model (J2-only). I’m integrating the STM via variational equations, where A(t) = ∂f/∂x is computed using autodiff.
Would appreciate any feedback — especially on whether you’ve seen similar nonlinear divergence in STM-based covariance propagation with only J2 involved.
