Root finding fails on event detector when event is very close to the current integration step

Orekit usage
1

Hi all! We came accross a failing of the root finding associated with event detection, when the following conditions are met:

  • The force field is non keplerian

  • The event we want to catch is very close to the current integration step: if the g function is going from positive to negative, then g is marginally negative at the current integration step, we are “right after the event”.

  • The detector is going through an event shifter with quite large time shift.

In the hipparchus routine DetectorBasedEventState#evaluateStep, a change of sign is detected and a root finding is triggered. Inside the root finder, the g function is again evaluated at the current integration step. However, the state used is not exactly the same as in when the change of sign is detected: the later is the soft current state sb, while the former is interpolator.getInterpolatedState(sb.getTime()). For some reason, these two states are not exactly equivalent: the primary state are equals to the last decimals, but not the primary derivatives. When these two states go through the event shifter with a quite large time shift, this non-keplerian difference yields variations on the 6-dimensional state, big enough to loose the change of sign, and an exception is raised.

Hereafter is a unit test to reproduce the bug:

@Test
    public void reproduceBracketErrorBug() {
        ClassicalRungeKuttaIntegratorBuilder integratorBuilder = new ClassicalRungeKuttaIntegratorBuilder(1.);
        double positionScale = 10e3;
        double a = 6879290;
        double e = 0.0017;
        double i = 1.7;
        double pa = -5.3;
        double raan = 0.42;
        double anomaly = 3.04;
        AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
        Frame frame = DataContext.getDefault().getFrames().getTOD(false);
        KeplerianOrbit orbit = new KeplerianOrbit(a, e, i, pa, raan, anomaly, PositionAngleType.TRUE, frame, date, Constants.EGM96_EARTH_MU);
        NumericalPropagatorBuilder numericalPropagatorBuilder = new NumericalPropagatorBuilder(orbit, integratorBuilder, PositionAngleType.TRUE, positionScale);
        NumericalPropagator numericalPropagator = numericalPropagatorBuilder.buildPropagator();
        HolmesFeatherstoneAttractionModel geoPotentialModel = new HolmesFeatherstoneAttractionModel(DataContext.getDefault().getFrames().getITRF(IERSConventions.IERS_2010, false),
                                                                                                    DataContext.getDefault().getGravityFields().getNormalizedProvider(2, 0));
        numericalPropagator.addForceModel(geoPotentialModel);
        EventDetectionSettings detectionSettings = new EventDetectionSettings(AbstractDetector.DEFAULT_MAX_CHECK, AbstractDetector.DEFAULT_THRESHOLD,
                                                                              AbstractDetector.DEFAULT_MAX_ITER);
        EventOnIntegratorTimeStepDetector detector = new EventOnIntegratorTimeStepDetector(detectionSettings);
        EventShifter eventShifter = new EventShifter(detector, false, 500, 500);
        numericalPropagator.addEventDetector(eventShifter);
        numericalPropagator.propagate(date.shiftedBy(10));
    }


private class EventOnIntegratorTimeStepDetector extends AbstractDetector {

        protected EventOnIntegratorTimeStepDetector(EventDetectionSettings detectionSettings) {
            super(detectionSettings,  new ContinueOnEvent());
        }

        @Override
        public double g(SpacecraftState s) {
            double exactY = -2391379.2800254025;
            return s.getOrbit().getPosition().getY() - exactY + 1e-8;
        }

        @Override
        protected AbstractDetector create(EventDetectionSettings detectionSettings, EventHandler newHandler) {
            return new EventOnIntegratorTimeStepDetector(detectionSettings);
        }
    }

The expected stacktrace is

org.orekit.errors.OrekitException: org.orekit.propagation.integration.AbstractIntegratedPropagator$AdaptedEventDetector@563f38c4 failed to find root between 1 (g=-1.854245212844352E3) and 2 (g=1.0E-8)
Last iteration at 2 (g=-1.441338860988617E-7)

	at org.orekit.errors.OrekitException.unwrap(OrekitException.java:154)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:572)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:475)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:435)
	at org.orekit.propagation.numerical.NumericalPropagatorTest.reproduceBracketErrorBug(NumericalPropagatorTest.java:132)
	at java.base/java.lang.reflect.Method.invoke(Method.java:569)
	at java.base/java.util.ArrayList.forEach(ArrayList.java:1511)
	at java.base/java.util.ArrayList.forEach(ArrayList.java:1511)
Caused by: org.hipparchus.exception.MathIllegalStateException: org.orekit.propagation.integration.AbstractIntegratedPropagator$AdaptedEventDetector@563f38c4 failed to find root between 1 (g=-1.854245212844352E3) and 2 (g=1.0E-8)
Last iteration at 2 (g=-1.441338860988617E-7)
	at org.hipparchus.ode.events.DetectorBasedEventState.findRoot(DetectorBasedEventState.java:373)
	at org.hipparchus.ode.events.DetectorBasedEventState.evaluateStep(DetectorBasedEventState.java:231)
	at org.hipparchus.ode.AbstractIntegrator.lambda$acceptStep$4(AbstractIntegrator.java:328)
	at java.base/java.util.ArrayList$ArrayListSpliterator.forEachRemaining(ArrayList.java:1625)
	at java.base/java.util.stream.Streams$ConcatSpliterator.forEachRemaining(Streams.java:734)
	at java.base/java.util.stream.ReferencePipeline$Head.forEach(ReferencePipeline.java:762)
	at org.hipparchus.ode.AbstractIntegrator.acceptStep(AbstractIntegrator.java:328)
	at org.hipparchus.ode.nonstiff.FixedStepRungeKuttaIntegrator.integrate(FixedStepRungeKuttaIntegrator.java:168)
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:550)
	... 6 more
Caused by: org.hipparchus.exception.MathIllegalArgumentException: interval does not bracket a root: f(1E0) = -1.854245212844352E3, f(2E0) = -144.13388609886168E-9
	at org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver.doSolveInterval(BracketingNthOrderBrentSolver.java:209)
	at org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver.solveInterval(BracketingNthOrderBrentSolver.java:427)
	at org.hipparchus.analysis.solvers.BracketedUnivariateSolver.solveInterval(BracketedUnivariateSolver.java:126)
	at org.hipparchus.ode.events.DetectorBasedEventState.findRoot(DetectorBasedEventState.java:364)
	... 14 more
2

I ran your test and it works on my machine. Maybe due to differences in EOP or J2 value used? I used the values in src/test/resources/{regular-data,potential}.

Perhaps provide the value of C_2,0 you’re using? Do you still see the problem with zero EOP, or does it require a certain EOP data set to reproduce?

Easiest on the developers would be to integrate your test into the Orekit source tree, make sure it fails when run using mvn clean test and provide the patch.

3

Hello Evan,

Thanks for your reply! I launched the test in NumericalPropagatorTest.java, in which “regular-data:potential/shm-format” is used. Does this answer you question?

4

Hi,

Could this be related to this issue?

Cheers,
Romain.

5

Hello Romain,

Thanks for your answer!
I don’t think so, the issue discussed here is not associated with a large sum of the shifts (increasing - decreasing).

6

Hi Bastien,

like Evan, I’ve run your code (with Orekit 13.1.2) with a different test env and there’s no Exception raised. I confirm that it fails with the one from NumericalPropagatorTest though. Strange…

Cheers,
Romain.

7

Quick update:

  • the failure is very sensitive to the EOP. Swtiching flags in TOD and ITRF can make the test pass. I guess it’s because then the integration step is not close anymore to the event
  • the detection still crashes without an Eventshifter, but a normal detector using this as g:
public double g(SpacecraftState s) {
            double exactY = -2391379.2800254025;
            return s.shiftedBy(-500.).getPosition().getY() - exactY + 1e-8;
        }
8

Would it be possible to put the unit test for this bug in CloseEventsAbstractTest? That way we can test the bug with many propagators and integrators, not just the defaults in NumericalPropagator.

9

We could but I’m not sure it’d be useful because I believe that Bastien pretty much engineered this test around the issue. Basically the integration step must land right next to the shifted event.
Giving it some thought, I think that shifted events are bound to run into difficulties with their detection. Indeed shiftedBy for Orbit relies on Keplerian motion (+ some corrections derived from the acceleration), so if the propagation has perturbations, the two predictions will differ for large timespans, possibly with different signs. At this point I don’t know how to prevent that, but I’m actually wondering if we should provide a fix at all. Maybe there should just be a warning for users about using shiftedBy in the event function…

10

I’m still trying to reproduce it to see what is really going on. The test passes when I run it locally in the NumericalPropagatorTest.

From the exception, unless there is a bug, it appears the g function changed its definition when it was not allowed to according to [1]. The part saying, “2 (g=1.0E-8) Last iteration at 2 (g=-1.441338860988617E-7)” indicates that two different values with different signs were returned for g(2).

Shifting should be fine because it is a function of the state passed to the g function. So passing the same state to g should result in the same value of the g function. But for some reason here it seems to be a function of something else as well.

There are some warnings in [1], but it could be helpful to highlight them or reference them from other parts of the documentation.

[1] org.hipparchus.ode.events (Hipparchus 3.1 API)