DerivativeStructures with AdaptiveStepsizeIntegrator for reentry

Orekit usage
1

Hi everyone!

I am currently trying to implement Orekit’s Taylor Algebra functionality into a reentry propagating tool. This tool, in it’s current (non-Taylor Algebra) version, uses a DormandPrince853Integrator, which is generally a good choice for a variable stepsize method. I have successfully used this integrator (or rather its “field” version) for Taylor Algebra-based propagations in the past, but never for reentry.

The problem comes when propagating until an AltitudeDetector (Alt = 0.) event, where the derivatives in the DerivativeStructure variables become null (real order-zero value is still kept). I have made several observations about this:

  • If the AltitudeDetector is set to a considerably higher altitude (e.g. 10 km), derivatives are still there, suggesting that numerics are not stable when the orbit becomes too vertical.
  • Every other variable stepsize method integrator available in Hipparchus suffers the same problem
  • If, however, a fixed stepsize method is chosen instead (RungeKutta in any of its variations), propagation can take place until Alt = 0 and derivatives are kept. This obviously lacks all the benefits of a variable stepsize integrator.

For context, I use a FieldNumericalPropagator with a Holmes-Featherstone attraction model, Sun and Moon as 3rd bodies, solar radiation pressure and an implementation of IsotropicDrag that let’s me input a DerivativeStructure as drag coefficient. Minimum and maximum stepsize for the DormandPrince integrator are respectively [1e-20 s, 3e02 s] and position tolerance is set to 1e-8 m.

Is there any chance my config is wrong here? If not, I would really appreciate suggestions to avoid or improve the fixed-step method usage.

Thank you and cheers,

Mviturro

2

I don’t know for sure, but my wild guess would be that the behaviour shown by variable stepsize algorithms is normal and behaviour shown by fixed step size is just good luck.

From my experience, computing reentry trajectories often needs special handling for two force models: gravity field and atmospheric drag. Spherical harmonics expansion are singular when crossing the Brillouin sphere. Indeed, I would even have guessed you would trigger an exception as soon as the distance to Earth center comes below equatorial radius, even if you are still above the ellipsoid if reentry is not at equator. Atmospheric drag relies on models that often have a lower altitude limitation in the range of 90km or 120km. The NRLMSISE2000 model can go down to ground as can the simple exponential model, though.

As you put the position tolerance to a really small (I would say unrealistic) value of 1.0e-8m, there are chances that when approaching the singularity, the step sizes decreases a lot and you end up evaluating points too close to the singularity and see the derivatives go weird. On the other hand, a fixed step size integrator would most probably not evaluate anything exactly at singularity, but only points slightly away from the singularity, hence not noticing it.

3

As Luc said. I found exactly the same thing happening. A custom ForceModel below the reentry interface altitude could help.

4

Thank you both for the answers and sorry for the delay. Indeed, the behaviour that you describe seems to match the different results I’m getting from fixed and variable step size. The atmospheric model was NRLMSISE2000 all along, which explains why it works even at low altitudes. I’ve also observed that lowering the position tolerance to 1e-3 or 1e-2 allows the variable stepsize solver to (sometimes) “jump” over the singularity. As for the inclusion of a custom ForceModel, I think a solution as simple as changing the gravitational force model to a simple Newton force once it crosses a certain altitude might be good enough for me.

5

Hi all again.

After several days of testing, I come back to you with what I think is the same problem, but this time I’ve investigated it a bit more. I think it’s worth it to re-formulate the problem since last time I hadn’t checked against other derivative computation methods. I have also written a small piece of code that clearly shows the difference between the two and is easily modifiable to see the problem appear and disappear as the configuration is changed (propagation time, initial altitude, etc).

So again, I’m propagating from a low altitude (initial semi-major axis 6550 km and nearly zero eccentricity) for 6000 seconds with a FieldNumericalPropagator. As in the previous case, the integrator I’m using is a DormandPrince853 (field version), but the problem also happens with other variable stepsize integrators available in Hipparchus. Other orbit parameters are basically left at random. As for propagation parameters, tolerances and stepsizes limits are as described in my first post, but changing them to more “reasonable” values doesn’t seem to fix the problem (this is easy to change in the example provided down below). Force models applied are the following:

  • Newtonian grav. attraction
  • HolmesFeatherstone Attraction Model (18,18).
  • Sun and Moon as third bodies.
  • Isotropic Drag configured with a HarrisPriester atmosphere (S = 10m2, Cd = 2.0).
  • SolarRadiationPressure configured with IsotropicRadiationCNES95Convention.

To monitor the value of the derivatives, the ones from the first component of the PV array (x position) are extracted from the DerivativeStructure object (in the Field version). For comparison, a PartialDerivativeEquations object is used in the Non-Field propagation and the first row of the Jacobian obtained after propagation is checked against the extracted derivatives in the Field version.

The basic problem is that after propagating for a certain period of time with the Field integrator with an adaptive stepsize, the derivatives of the PV variables with respect to the initial PV diverge. Trying to isolate the problem has been unsuccessful for me, since it seems to depend on a lot of different parameters. I’ll try to list them as thoroughly as I can:

  • First of all, the problem dissapears when the Drag and SolarRadiationPressure are not added. This is already a big hint that the problem has to do with at least one of these ForceModels.
  • Most of the times, leaving the SolarRadiationPressure on and taking Drag off the propagator fixes the problem (but it’s obsviously not well suited to re-entry cases).
  • Derivatives diverge when the Drag is active, under the cited conditions, but if SolarRadiationPressure is also active the resulting derivatives are even bigger (by some orders of magnitude). This is why I don’t completely rule out SRP as a source of problems.
  • As said in my first post, either making the tolerance big enough in the variable stepsize integrator or integrating with the fixed stepsize one (with no regards to tolerances and a fixed step of around 1 second) makes the derivatives “not diverge”, but the results are also not the same as the ones obtained though the PDE object (which seem to be the best guess of the three).
  • Altitude is a key factor here. Contrary to what I thought when I started this post, it’s not Altitude = 0 that creates the problem. This happens even before the HarrisPriester model runs out of altitude (100 km ). At that altitude, all other models should be comfortable enough to not arrive at any singularities. However, at a higher altitude (initial semi-major axis = 6850 km) the problem goes away.
  • I have also tested this with a SimpleExponentialAtmosphere, and the problem is still there, it only happens at a different altitude.

This quickly-made piece of code should illustrate the issue well enough. What it does is create both propagators (field and non-field) basically side-by-side, with the same input and observing the “same” output. Orbital parameters are set at the beggining and the propagation time (which can be lowered to see the effect on the derivatives) is the same for both cases.

FieldBugTest.java (9.7 KB)
(It will look for the orekit-data folder in the Home directory, but appart from that it should work just fine.)

The output of this code are two printed vectors, showing the values of the previously mentioned derivatives side-by-side. If run as it is, the code will produce vastly different results from both methods. If the initial semi-major axis is increased 100 km, the problem virtually dissapears and both methods produce very similar results.

I’m still not sure this is not something I’ve done rather than a bug, but in that case I guess you can tell me what’s wrong with my whole setup.

Thank you very much in advance. Best regards,

MViturro

6

I was able to reproduce the issue and confirm that using a pair of LutherIntegrator/LutherFieldIntegrator
with a step size of 60.0 seconds, the derivatives are similar but with DormandPrince853Integrator/DormandPrince853FieldIntegrator the derivatives computed with fieds are unrealistic (by a factor 1.0e18…).

Could you open a bug report on the forge at https://gitlab.orekit.org/orekit/orekit/issues?