I am currently trying to propagate a state using the AdamsBashforthIntegrator.
During the propagation when the propagator reaches the last step I obtain the following exception:
multistep integrator starter stopped early, maybe too large step size.
As it turned out this happens only when I add the Solar Radiation Pressure to the force models of the Propagator.
I have tried to propagate using a different integrator (ClassicalRungeKuttaIntegrator) but in this case no exception appears.
The satellite is on a LEO, therefore I thought that the contribution of the SRP could create problems for the propagator to reach convergence.
I tried to propagate a GEO satellite with the same SRP parameters and in this case the exception does not appear.
Still I don´t understand why in this specific case and with this specific integrator the propagator throws this exception.
I am also using another software to compare the results with the same integrator and the same parameters and in this case the integration proceeds well.
I was not able to find any information related to this behavior, if you know the reason I would be really grateful for any help.
The reason for this behavior is that solar radiation pressure generates discontinuous acceleration at eclipses umbra entry/exit and penumbra entry/exit. These discontinuities are handled by automatically setting up internal event detectors that trigger a RESET_DERIVATIVES action when the events occur.
This reset is forwarded from the numerical propagator to the underlying ODE integrator, which enforces a new step to start at this point. This is straightforward for single step integrators like ClassicalRungeKuttaIntegrator or DormandPrince853Integrator that do not preserve memory of previous steps and for which starting a new step doesn’t involve any fancy computation. This is not as easy for multistep integrators like AdamsBashforthIntegrator or AdamsMoultonIntegrator as they need to wipe away information gathered on previous steps and rebuild a new model (this information is stored in a Nordsieck vector). For rebuilding this model, multistep integrators use an underlying starter integrator
which by default is a DormandPrince853Integrator. This restarting procedure requires the underlying integrator to be able to run for (n + 3) / 2 steps where n is the number for steps of the Adams method used. In your case, n = 8 so the starter integrator must compute 5 steps.
What happens in your case is that due to bad luck, your propagation range goes from 00:00:00.000 to 06:00:00.000 (i.e. 0.25 day) and a penumbra exit event occurs at 05:59:56.973, i.e. 3.027 seconds before propagation end. There are several other umbra/penumbra entry/exit events before, but this last one is the culprit. After this event the Adams-Bashforth integrator experience a restart, but the underlying Dormand-Prince integrator cannot perform 5 steps in the short 3.027 seconds range that remains before integration end. This triggers the exception.
Note also that your setup for step size handling could be improved. You use constants 1e-4 for absolute and relative tolerances on the integrator and do not specify the orbit type to use for integration, which leads to OrbitType.EQUINOCTIAL to be selected by default. We recommend to set up the orbit type explicitly and to use it for both configuring the propagator and the tolerances in adaptive stepsize. You can do it using this code snippet in your case:
OrbitType type = OrbitType.EQUINOCTIAL;
double[][] tol = NumericalPropagator.tolerances(0.1, orbit, type);
AdamsBashforthIntegrator integrator = new AdamsBashforthIntegrator(8, 1e-15, timeSpan, tol[0], tol[1]);
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setInitialState(state);
propagator.setOrbitType(type);
Finally, when running your test case, I found there was indeed a bug in the way Hipparchus handles the reset and hence the restarts. The resetOccurred flag is not reset to false after the reset has been handled and spurious restarts occur between events, greatly increasing computation costs with multistep integrators. Could you open a bug report on https://github.com/Hipparchus-Math/hipparchus/issues for this? I will fix this bug soon but as a workaround, you can just add resetOccurred = false; before the boolean doneWithStep = false; statement at line 314 in Hipparchus AbstractIntegrator class.
So as a conclusion, there are three things you could do:
extend your propagation time so it does not end just after a discontinuity-inducing event
configure tolerances in step size management using NumericalPropagator.tolerances
and an explicit OrbitType
reset resetOccurred in AbstractIntegrator to avoid too many restarts
many many thanks for your reply, I have reported that bug in GitHub. However I have noticed that I don’t have the doneWithStep flag in my AbstractIntegrator at line 314.
My dependency is set to Orekit 9.3 v and the Orekit pom points to the 1.4 v of Hipparcus, which does not have yet the changes introduced in December.
Following your suggestions I have amended my class with the NumericalPropagator.tolerances method.
Furthermore I have added two detectors to detect the penumbra and umbra events that cause the discontinuities.
I have also added two “custom” handlers to my detectors with the return value equal to Action.CONTINUE.
However, I don’t detect a penumbra event at 05:59:56.973, but at “04:58:59.317” (maybe my detectors are set differently or in a wrong way).
So I tried to propagate to the same time and in this case no exception appear.
Then I changed the return value of the handlers to Action.RESET_DERIVATIVES and the exception appears again.
It seems strange to me that just by changing the handlers of the detectors I am able to reach or not the end of the propagation.
Still I don’t know which is the influence on the final result when using the handler with Action.CONTINUE.
Do you know why this is happening?
In addition I would like to ask you if the issue with the handling of the discontinuities is an intrinsic property of the multistep integrators or it is related to the way in which it was implemented ?
Because if this is connected to the integrator I can consider to switch to a different one (i.e. as you suggested the DormandPrince853Integrator), in order to protect the user while propagating.
Furthermore I have added two detectors to detect the penumbra and umbra events that cause the discontinuities.
The detectors you added are just new ones the propagator will handle in addition to the internal ones from the force model. You can do whatever you want with your detectors, they will not prevent the internal detectors to return Action.RESET_DERIVATIVES.
However, I don’t detect a penumbra event at 05:59:56.973, but at “04:58:59.317”
The event date moves slightly with tolerances in the integration. While playing with you initial code, I have
seen changes of a handful of seconds, so there is an event near 06:00, either just before or just after. Note
also that there are other eclipse events before the end (four events per orbit, penumbra entry, umbra entry, penumbra exit, umbra exit) and that the time I wrote was in UTC.
Still I don’t know which is the influence on the final result when using the handler with Action.CONTINUE .
If you would use Action.CONTINUE on the internal force model event that really induces a derivative change (i.e. not on an additional event detector), then the propagator would miss the fact that the dynamics changed. It would most probably see something is wrong during step size computation and would reduce step size over and over until the step around the derivative discontinuity is so small that it finally accept this step. It would work but would cost a lot of additional computation. Action.RESET_DERIVATIVES is a way to notify the integrator that the right hand side of the differential equation has a discontinuity and that the integrator must put a step boundary here instead of trying to find a small continuous step that contains the discontinuity.
In addition I would like to ask you if the issue with the handling of the discontinuities is an intrinsic property of the multistep integrators or it is related to the way in which it was implemented ?
I would say both are true. Multistep integrators don’t manage derivatives discontinuities because they basically set up a polynomial model based on previous steps, and polynomials are continuous, so they don’t match models with derivatives discontinuities. Basic multistep integrators can be tweaked to manage discontinuities, but this implies a restart. The way we handle restart, which is related to our implementation choices, introduces the need for a “rest time” after the event and before the end of the integration. I think it could theoretically be leveraged if we could enforce starting steps small enough to fit within the available time, but first we cannot do it now (we cannot reset step size control of the starter integrator) and second this would not work if there are two consecutive events, which is even much more difficult to foresee.
many many thanks for your detailed explanation, yes as you said there is an event slightly after the 6 hours.
I think that I will switch to a single step Integrator to avoid further issues.