Action.RESET_DERIVATIVES occurring at the end of integration can cause interpolation in DenseOutputModel to produce NaN

[andrewsgoetz]

When running the following snippet:

// ------------------------------------------------------
// Integrate the IVP ẏ = 1, y(0) = 0 from t = 0 to t = 1.
// ------------------------------------------------------

// Define the initial value problem.
final double ti = 0.;
final double tf = 1.;
final double yi = 0.;
final ODEState initialODEState = new ODEState(ti, new double[] { yi });
final OrdinaryDifferentialEquation ode = new OrdinaryDifferentialEquation() {

    @Override
    public int getDimension() {
        return 1;
    }

    @Override
    public double[] computeDerivatives(final double t, final double[] y) {
        return new double[] { 1. };
    }

};

// Construct an integrator. Use a DenseOutputModel as a step handler.
final ODEIntegrator integrator = new ClassicalRungeKuttaIntegrator(0.25);
final DenseOutputModel denseOutputModel = new DenseOutputModel();
integrator.addStepHandler(denseOutputModel);

// Add an event handler which triggers at the final time and which does Action.RESET_STATE.
final ODEEventHandler odeEventHandler = new ODEEventHandler() {

    @Override
    public double g(final ODEStateAndDerivative state) {
        return state.getTime() - tf;
    }

    @Override
    public Action eventOccurred(final ODEStateAndDerivative state, final boolean increasing) {
        return Action.RESET_DERIVATIVES;
    }

};
integrator.addEventHandler(odeEventHandler, 0.5, 1e-10, 100);

// Perform the integration.
final ODEStateAndDerivative finalODEStateAndDerivative = integrator.integrate(ode, initialODEState, tf);

// The value of y(tf) obtained from the integrator is (close to) expected.
System.out.println("y(tf) from integrator: " + finalODEStateAndDerivative.getPrimaryState()[0]);

// The value of y(tf) obtained from the DenseOutputModel is NaN.
System.out.println("y(tf) from DenseOutputModel: " + denseOutputModel.getInterpolatedState(denseOutputModel.getFinalTime()).getPrimaryState()[0]);

the output is

y(tf) from integrator: 0.9999999999999999
y(tf) from DenseOutputModel: NaN

Although this example is rather artificial (why reset derivatives at the end of integration?) it was extracted from this more realistic example (using Orekit):

// Define a circular equatorial LEO orbit.
final AbsoluteDate epoch = AbsoluteDate.ARBITRARY_EPOCH;
final KeplerianOrbit orbit = new KeplerianOrbit(7000000., 0., 0., 0., 0., 0., PositionAngle.TRUE, FramesFactory.getGCRF(), epoch, Constants.IERS2010_EARTH_MU);
final SpacecraftState initialState = new SpacecraftState(orbit);

// Define a simple fixed-step propagator.
final ODEIntegrator integrator = new ClassicalRungeKuttaIntegrator(10.);
final NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setInitialState(initialState);
propagator.setOrbitType(OrbitType.CARTESIAN);

// Define a 60-second maneuver.
final AbsoluteDate maneuverStartDate = epoch;
final double maneuverDuration = 60.;
final double thrust = 1.; // Newtons
final double isp = 500.; // seconds
final AbsoluteDate maneuverStopDate = epoch.shiftedBy(maneuverDuration);
final ConstantThrustManeuver maneuver = new ConstantThrustManeuver(epoch, maneuverDuration, thrust, isp, Vector3D.PLUS_I);
propagator.addForceModel(maneuver);

// Propagate using ephemeris mode and generate ephemeris.
propagator.setEphemerisMode();
propagator.propagate(maneuverStartDate, maneuverStopDate);
final BoundedPropagator ephemeris = propagator.getGeneratedEphemeris();

// Using the generated ephemeris, print out the final spacecraft state.
System.out.println(ephemeris.propagate(maneuverStopDate));

which prints out

SpacecraftState{orbit=Cartesian parameters: {P(NaN, NaN, NaN), V(NaN, NaN, NaN)}, attitude=org.orekit.attitudes.Attitude@53ca01a2, mass=NaN, additional={}}

The first snippet is intended to capture the underlying basics of what is going on in the Orekit example which cause the NaNs. In the Orekit example, the event occurring at the end of integration is the end of the maneuver, which causes a reset of derivatives. From now on I will discuss the first Hipparchus example.

Digging deeper, what happens is that the event at the end of the integration causes a degenerate ODEStateInterpolator to be added to the end of the internal list maintained in the DenseOutputModel whose "previous" and "current" states are the same and both at the end time of the integration. The call to DenseOutputModel.getInterpolatedState(final double time) delegates to the method AbstractODEStateInterpolator.getInterpolatedState(final double time) of this final degenerate interpolator. Here is the current implementation of AbstractODEStateInterpolator.getInterpolatedState(final double time):

@Override
public ODEStateAndDerivative getInterpolatedState(final double time) {
    final double thetaH         = time - globalPreviousState.getTime();
    final double oneMinusThetaH = globalCurrentState.getTime() - time;
    final double theta          = thetaH / (globalCurrentState.getTime() - globalPreviousState.getTime());
    return computeInterpolatedStateAndDerivatives(mapper, time, theta, thetaH, oneMinusThetaH);
}

The NaN appears at line 166 as theta is computed as 0/0.

As for the fix: I was looking into ways of stopping the degenerate ODEStateInterpolator from being added to the DenseOutputModel (for example, in AbstractIntegrator.acceptStep(final AbstractODEStateInterpolator interpolator, final double tEnd)) copying the check for isLastStep from line 425 to occur also at line 359) but my attempts caused various unit tests to fail. The easiest fix I have found is to allow the degenerate ODEStateInterpolator to be added to the DenseOutputModel and to just add simple check to the implementation of AbstractODEStateInterpolator.getInterpolatedState(final double time):

@Override
public ODEStateAndDerivative getInterpolatedState(final double time) {
    if (globalPreviousState.getTime() == globalCurrentState.getTime()) {
        return globalCurrentState;
    }
    final double thetaH         = time - globalPreviousState.getTime();
    final double oneMinusThetaH = globalCurrentState.getTime() - time;
    final double theta          = thetaH / (globalCurrentState.getTime() - globalPreviousState.getTime());
    return computeInterpolatedStateAndDerivatives(mapper, time, theta, thetaH, oneMinusThetaH);
}

I can submit this change as a pull request if it is agreed that this is a bug and the fix is good.

[maisonobe]

Hi @andrewsgoetz

Your analysis seems fine to me. I wonder what made units test fail with the check. Could you set up a first merge request with this first fix attempt, it would be interesting to understand what happens there.

Also for the final merge request, beware that it would be better to check times within one ULP, and that the same fix must be done in the field version of the integrator.

Anyway, its great to have identified this bug!

[andrewsgoetz]

I have made a merge request for my first fix attempt which did not work: https://github.com/Hipparchus-Math/hipparchus/pull/113.

The failing tests are:

DormandPrince54IntegratorTest.testUnstableDerivative DormandPrince853IntegratorTest.testUnstableDerivative HighamHall54IntegratorTest.testUnstableDerivative EventFilterTest.testTwoOppositeFilters EventFilterTest.testIncreasingOnly EventFilterTest.testDecreasingOnly CloseEventsTest.testSimultaneousDiscontinuousEventsAfterReset CloseEventsTest.testSimultaneousDiscontinuousEventsAfterResetReverse FieldCloseEventsTest.testSimultaneousDiscontinuousEventsAfterReset FieldCloseEventsTest.testSimultaneousDiscontinuousEventsAfterResetReverse

If you don't see any obvious way to patch this up I will make a second merge request with my fix that changes AbstractODEStateInterpolator and AbstractFieldODEStateInterpolator.