Detector failing to find root

Emiliano · September 28, 2023, 12:12pm

Hello everyone,

I recently encountered this error:

JavaError: <super: <class 'JavaError'>, <JavaError object>>
    Java stacktrace:
org.orekit.errors.OrekitException: org.orekit.propagation.events.EventSlopeFilter@4d94fc95 failed to find root between 2024-03-17T05:55:05.5923675099728Z (g=6.51925802230835E-9) and 2024-03-17T06:26:43.72339071493213Z (g=-5.943172694678251E6)
Last iteration at 2024-03-17T05:55:05.59236751003797Z (g=6.51925802230835E-9)
	at org.orekit.propagation.events.EventState.findRoot(EventState.java:359)
	at org.orekit.propagation.events.EventState.evaluateStep(EventState.java:222)
	at org.orekit.propagation.analytical.AbstractAnalyticalPropagator.acceptStep(AbstractAnalyticalPropagator.java:297)
	at org.orekit.propagation.analytical.AbstractAnalyticalPropagator.propagate(AbstractAnalyticalPropagator.java:151)
	at org.orekit.propagation.AbstractPropagator.propagate(AbstractPropagator.java:276)
Caused by: org.hipparchus.exception.MathIllegalStateException: maximal count (100) exceeded
	at org.hipparchus.util.Incrementor.lambda$static$0(Incrementor.java:41)
	at org.hipparchus.util.Incrementor.increment(Incrementor.java:237)
	at org.hipparchus.analysis.solvers.BaseAbstractUnivariateSolver.incrementEvaluationCount(BaseAbstractUnivariateSolver.java:318)
	at org.hipparchus.analysis.solvers.BaseAbstractUnivariateSolver.computeObjectiveValue(BaseAbstractUnivariateSolver.java:165)
	at org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver.doSolveInterval(BracketingNthOrderBrentSolver.java:297)
	at org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver.solveInterval(BracketingNthOrderBrentSolver.java:427)
	at org.hipparchus.analysis.solvers.BracketedUnivariateSolver.solveInterval(BracketedUnivariateSolver.java:126)
	at org.orekit.propagation.events.EventState.findRoot(EventState.java:350)
4 more

For context: I’m using a NodeDetector on an AggregateBoundedPropagator (obtained via EphemerisGenerator of a NumericalPropagator). This is how I create the detector:

... other stuff ...

ephemeris = self.ephemeris
initState = ephemeris.getInitialState()
ephemeris.clearEventsDetectors()
logger = EventsLogger()
nodeCrossing = NodeDetector(initState.getOrbit(), initState.getFrame()).withHandler(ContinueOnEvent())
nodeCrossing = EventSlopeFilter(nodeCrossing, FilterType.TRIGGER_ONLY_DECREASING_EVENTS)
nodeCrossing = logger.monitorDetector(nodeCrossing)
ephemeris.addEventDetector(nodeCrossing)
ephemeris.propagate(ephemeris.getMaxDate())
nodes = list(logger.getLoggedEvents())

... other stuff ...

I don’t understand what’s happening. I also don’t understand why the two dates reported are so distant from each other. Are those the dates of the two closest ephemeris saved by the generator? Also, I feel like the first value of g reported is close enough to the real root. Is there a way to change the tolerance on the root-search?

As a side note, I’ve encountered similar issues in the past but I don’t really know how to reproduce them. Sometimes just removing the filter on the detector and filtering afterwards based on the isIncreasing method is enough to sidestep the problem.

Best regards,
Emiliano

Emiliano · September 28, 2023, 2:17pm

Update:

Sometimes the detector doesn’t log anything on the logger. I checked the initial state of the BoundedPropagator, just to make sure that I was actually propagating to a different date, and I am. Again, I have no idea what might be causing this

Serrof · September 29, 2023, 7:27am

Hi Emiliano,

You can change the detection tolerance via the threshold parameter. Use the withThreshold constructor to set it on top of existing configuration.

Edit: my bad, that will update the tolerance in seconds on time, not on g. Have you tried increasing the maximum number of iterations via withMaxIter?

Best,
Romain

MaximeJ · September 29, 2023, 8:24am

Hi @Emiliano,

Not sure what’s happening here either.

Could you please send us a small self-contained example reproducing the behaviour?
That way we can investigate if it’s a bug or a configuration issue.

Thanks,
Maxime

Emiliano · September 29, 2023, 11:28am

Hello @Serrof,
I tried your suggestion with a maxCheckIter of 1000 but it still failed. Interestingly enough the error log still reported an iteration counter of 100

Emiliano · September 29, 2023, 11:41am

Hello @MaximeJ,

I noticed that this weird behaviour happens when the propagation of the ephemeris is done in a class instead of the main code.

This is a snippet of code that should reproduce the second behaviour I mentioned (the logger is empty):

import orekit;
orekit.initVM();

from orekit.pyhelpers import setup_orekit_curdir;
setup_orekit_curdir();

from org.hipparchus.ode.nonstiff import DormandPrince853Integrator;

from org.orekit.attitudes import BodyCenterPointing;
from org.orekit.bodies import OneAxisEllipsoid;
from org.orekit.frames import FramesFactory;
from org.orekit.orbits import CircularOrbit, PositionAngle;
from org.orekit.propagation import SpacecraftState;
from org.orekit.propagation.numerical import NumericalPropagator;
from org.orekit.propagation.events import EventSlopeFilter, EventsLogger, FilterType, NodeDetector;
from org.orekit.propagation.events.handlers import ContinueOnEvent;
from org.orekit.time import AbsoluteDate, TimeScalesFactory;
from org.orekit.utils import Constants, IERSConventions;

from orekit import JArray_double;
from math import acos, pi, sqrt;

# %% INITIALIZATION
# =============================================================================
# =============================================================================
initDate = "2024-01-01T12:00:00.0Z";
initDate = AbsoluteDate(initDate, TimeScalesFactory.getUTC());
day = Constants.JULIAN_DAY;

# Data
Re = Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
mu = Constants.WGS84_EARTH_MU;
J2 = -Constants.WGS84_EARTH_C20;
eme2000 = FramesFactory.getEME2000();
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True);
posAngle = PositionAngle.MEAN;

def computeSSOinclination(a):
    cosI = -4*pi/Constants.JULIAN_YEAR*a**(7/2)/(3*J2*Re**2*sqrt(mu));
    return float(acos(cosI));

# %% INITIAL STATE
# =============================================================================
# =============================================================================
earth = OneAxisEllipsoid(Re, Constants.WGS84_EARTH_FLATTENING, itrf);
earthPointing = BodyCenterPointing(eme2000, earth);

h0 = 600e3;
a0 = Re + h0;
ex0 = 1e-4;
ey0 = 1e-4;
i0 = computeSSOinclination(a0);
raan0 = 0.1;
u0 = 0.1;

initOrbit = CircularOrbit(a0, ex0, ey0, i0, raan0, u0, posAngle, eme2000, initDate, mu);
initState = SpacecraftState(initOrbit);

tol = NumericalPropagator.tolerances(1e-3, initOrbit, initOrbit.getType());
integrator = DormandPrince853Integrator(1e-3, 1e+3, JArray_double.cast_(tol[0]), JArray_double.cast_(tol[1]));
propagator = NumericalPropagator(integrator, earthPointing);
propagator.setResetAtEnd(False);
propagator.setInitialState(initState);

logger0 = EventsLogger();
nodeCrossing = NodeDetector(initOrbit, eme2000).withHandler(ContinueOnEvent());
nodeCrossing = EventSlopeFilter(nodeCrossing, FilterType.TRIGGER_ONLY_DECREASING_EVENTS);
nodeCrossing = logger0.monitorDetector(nodeCrossing);

timeSpan = 10*day;
finalDate = initDate.shiftedBy(timeSpan);
ephGen = propagator.getEphemerisGenerator();
propagator.addEventDetector(nodeCrossing);
finalState = propagator.propagate(finalDate);
ephemeris = ephGen.getGeneratedEphemeris();
nodes0 = list(logger0.getLoggedEvents());

logger1 = EventsLogger();
nodeCrossing = NodeDetector(initOrbit, eme2000).withHandler(ContinueOnEvent());
nodeCrossing = EventSlopeFilter(nodeCrossing, FilterType.TRIGGER_ONLY_DECREASING_EVENTS);
nodeCrossing = logger1.monitorDetector(nodeCrossing);

ephemeris.addEventDetector(nodeCrossing);
ephemeris.propagate(finalDate);
ephemeris.clearEventsDetectors();
nodes1 = list(logger1.getLoggedEvents());

class dummyClass():
    def getNodes(ephemeris):
        logger2 = EventsLogger();
        nodeCrossing = NodeDetector(initOrbit, eme2000).withHandler(ContinueOnEvent());
        nodeCrossing = EventSlopeFilter(nodeCrossing, FilterType.TRIGGER_ONLY_DECREASING_EVENTS);
        nodeCrossing = logger2.monitorDetector(nodeCrossing);
        
        ephemeris.addEventDetector(nodeCrossing);
        ephemeris.propagate(ephemeris.getMaxDate());
        ephemeris.clearEventsDetectors();
        
        nodes = list(logger2.getLoggedEvents());
        
        return nodes;

nodes2 = dummyClass.getNodes(ephemeris);

When I run this example, even with the simplest case possible (no perturbations, just Keplerian motion) nodes2 is empty, whereas nodes0 and nodes1 are not.

In my original post I am indeed retrieving the nodes in a method of a custom class. The ephemeris are stored internally and accessed by the method (that’s why they were defined as ephemeris = self.ephemeris). However, this doesn’t always happen though. I’ve been able to do things not so different from this case and they worked just fine.
Additionally, if I remove the first propagation of the ephemeris (the ones obtained via the EphemerisGenerator), then nodes2 is not empty.

Finally, I still haven’t been able to create a simple example to recreate the first issue (detector failing on root finding)

MaximeJ · September 29, 2023, 1:16pm

Hi again @Emiliano,

I think I solved it by changing this

ephemeris.propagate(ephemeris.getMaxDate());

By

ephemeris.propagate(initDate, ephemeris.getMaxDate());

in dummyclass().

Now nodes2 has the same events as the other objects.

While constructing nodes1 you propagate until the final date and thus the initial state in the ephemeris is the state at the final date.
Then when you propagate again, nothing happens since you’re already at the final date.

Emiliano · September 29, 2023, 1:44pm

That’s really weird. Yesterday I did check for this. I compared the date of the initial state in the ephemeris with the target date just to be sure that it was actually propagating. I don’t know, I must have done something wrong
Interestingly enough, specifying both the initial and final date when propagating seems to solve also the other issue. It doens’t fail anymore, even though I don’t understand why that would influence the convergence of the root-finding algorithm

On a side note, why doesn’t BoundedPropagator have a setResetAtEnd method? If it’s not possible to change it, I would expect it to maintain the same flag given to the propagator used to generate the ephemeris.

MaximeJ · September 29, 2023, 2:16pm

Interesting indeed. If you manage to come up with a simple example we could probably find the cause.

I think it could be added directly in the Propagator interface, along with a getter.
Then added as an attribute in AbstractPropagator and implemented in AbstractAnalyticalPropagator.
That way it will be available to all propagators in Orekit.
I think @Serrof already asked for a similar feature in an old thread.

Emiliano · September 29, 2023, 2:56pm

For sure, if I find a simple example to give to you I’d love to get an answer on what was going on.

Thanks again a lot for the help!

Serrof · September 29, 2023, 6:16pm

I would appreciate a getter yes. Adding a setter for everyone would require more though work as I believe that all non integration-based propagators automatically reset in the current implementation, with no other available mechanism.

Cheers,
Romain.

MaximeJ · September 29, 2023, 6:44pm

Right, a setter in an interface is probably not a good idea indeed

the idea was specifically to provide that mechanism or that choice for all propagators, including non-integration-based ones.
Maybe I don’t understand your point

Serrof · September 29, 2023, 8:31pm

I just meant that adding the getter is very easy, while the setter needs some actual work.

Romain.

Serrof · October 3, 2023, 8:05pm

I’ve created two separate issues about resetAtEnd :

for upcoming major release, add getter for integrator-based propagators. This is easy to fix, first contributors welcomed
for next minor, extend mechanism to other propagators

Cheers,
Romain