Numerical STM speed up options

[Serrof]

Hi folks,

Following discussions with @dgondelach, I’ve opened a couple of issues (1462 and 1463) regarding options for state transition matrix performance, starting with drag and geopotential. The idea is to use a simpler model for the derivatives, e.g. a less complex atmospheric model as some are very expensive with automatic differentiation.

Cheers,
Romain.

[evan.ward]

Hi Romain,

Could you provide a discussion of the limitations of that approach? Perhaps a reference that discusses the impacts on accuracy and convergence?

Rice tried the “hybrid” approach (different models for state and STM) in 1967 in [1] and concluded that it causes “rapid divergence”. It’s a good read. Based on his work, I wouldn’t recommend including the hybrid approach in a library like Orekit. But perhaps in the last 50+ years someone else has figured out how to address the divergence issues.

Regards,
Evan

[1] https://arc.aiaa.org/doi/10.2514/6.1967-123

[bcazabonne]

Hi!

It’s a good idea, but we need to be very careful with current orekit users.
If it’s just an option and the default behavior (i.e. calculation) remains as it is now, I think it’s acceptable to introduce this feature for users who favor calculation time over accuracy.
On the other hand, current users shouldn’t have to do anything - nothing should change for them. If an option is available, the action to be taken (i.e. a method to be called or other) must be taken by the users wishing to limit performance and increase computation speed.

Best regards,
Bryan

[luc]

I agree with Evan here.
From my experience, simplifying the force models when computing derivatives has a bad influence on convergence of orbit determination. It may help at start, while we are still far from the solution and just want to come closer. But when we are near convergence and globally derivatives are close to zero (because we are looking for the minimum of a cost function), it seems to me that having simplified derivatives just either makes the algorithm lose its path as the Jacobian is inconsistent with the evolution of the value or it makes the algorithm converge to the wrong solution.

[Serrof]

Hi all,

Thanks for your comments.

First of all, as Bryan said, it would be an option, and certainly not the default one.

As for accuracy, well it’s always a trade-off right? I’m pretty sure I’ve seen the option to cut off the geopotential order for STMs in GMAT, although I can’t get my hand on a proper link at the moment. As for drag, the proposed approximation is only on the density, which at the moment is with finite differences by default so I wouldn’t say it’s much more accurate. Besides, orbit determination is not the only place where one needs covariance matrices and/or STMs, and it’s not always in an iterative process requiring convergence. I believe users are not all doing operational flight dynamics. Anyhow if it doesn’t make consensus I’ll ditch the issues from the 12.2 milestone for now.

Cheers,
Romain.

[Serrof]

Here’s the option I was talking about in GMAT, field StmLimit:
https://documentation.help/GMAT/Propagator.html#Propagator_ForceModel

[afvy]

Hello,

I recently did some testing on the runtime performance of numerical STM propagation, and I wanted to revamp this discussion.

What I did was:

• tweak HolmesFeatherstoneAttractionModel and AbstractDragForceModel such that the speed-optimized implementations of acceleration with FieldSpacecraftState are called every time the field type is Gradient, and not only if the derivatives represent first-order independent variations in the components of the position vector (code is here Files · faster-partials · Alberto Fossà / Orekit · GitLab )
• Measured the runtime required to propagate state and STM for an SSO-like orbit at 500km altitude for 24h, using both the variational equations approach used by the OD and FieldNumericalPropagator<Gradient> (so “full” automatic differentiation)
• The force models I used were:
  • 48 x 48 non-spherical gravity potential
  • isotropic drag with NRLMSISE00 model (same density model for state and partials, both using AD and FD)
  • cannonball SRP

These are the runtimes in milliseconds, and the runtime ratio between variational and full AD approaches

1. gravitational attraction only:
   evaluating acceleration(double) is as expensive as evaluating accelerationWrtState(Gradient)
   as the variational approach evaluates both (primary and secondary ODEs), the field propagator is almost as twice as fast

   

2. drag only:
   evaluating acceleration(double) is way faster than evaluating the partials, so the advantage of not evaluating the main equations (in the variational sense) with the field propagator is less pronounced
   using FD does not make a huge difference compared to AD (i.e. Gradient)

   

3. SRP only:
   the results are reversed because of the eclipses event detection with FieldAbsoluteDate in the second case
   the partials of the event date w.r.t. the state are non-null with a field propagator, but these partials are neglected with the variational approach

   

4. gravitational attraction + drag + SRP:
   the full field approach is penalized by the fact that state epoch have now non-zero derivatives due to the eclipse event detection, so the speed-optimized implementation of the gravitational attraction cannot be used (as the gravitational attraction depends on the epoch via the GCRF → ITRF transformation)

So my thoughts are:

1. in the variational approach, as the secondary ODEs necessarily compute the primary acceleration, would it be possible to avoid evaluating the latter twice?
2. specifically for the gravitational potential, the similar runtime of primary and variational ODEs smells to me as a direct consequence of the recursive implementation of the Legendre polynomials used to model the spherical harmonics
   if this is indeed the case, I would think that it should be possible to implement custom differentiation rules for these polynomials such that the analytical partials become obsolete
3. would it be possible to avoid the use of event detectors for the eclipses? See for example the OPT_SRP_CORRECT flag in neptune GitHub - Space-Systems/neptune: NPI Ephemeris Propagation Tool with Uncertainty Extrapolation

[Serrof]

Hi Alberto,

thanks for the quantitative analysis!
Some comments about your points:

1. I guess, with some caching. Need to check the code in-depth tho, it might be tricky or at least very dirty
2. No idea
3. There are other ForceModel for SRP. The RadiationPressureModel class allows any LightFluxModel, including the cylindrical one that has only one detector instead of 2*N. Alternatively, with ForceModelModifier, you can easily override get(Field)EventDetectors and return an empty Stream.

Edit: for drag you can now have a different (cheaper) Atmosphere for the derivatives used by the STM.

Cheers,
Romain.

[evan.ward]

Re 1., I like the idea of evaluating value and partials at the same time. That is almost always faster. And I thought it was what happened now. That sounds like NumericalPropagator could be updated to avoid double evaluating the force model.

I have no doubt that the variational equations will be faster, if correctly implemented. The data show that the gravity model can be evaluated twice using only 1.7x the time as the Field version.