Consecutive propagation error

[metamon]

Hello.

I am performing a simple propagation for two days using orekit propagator.
I have encountered an interesting problem where the following two scenarios output different propagation results.

① propagate for 2 day

This is done using:
myPropagator.propagate(epoch0+dt.timedelta(hours = 2.)) where epoch0 is the initial epoch
myPropagator has Force model which only include NewtonianAttraction.

② propagate for 1 day → stop → propagate for another 1 day

This is done using:
myPropagator.propagate(epoch0+dt.timedelta(hours = 1.))
last_propagation_epoch= after_prop_state.getDate()
dV_execution_state_REF = myPropagator.propagate(absolutedate_to_datetime(last_propagation_epoch)+dt.timedelta(hours=1.))

I have compared the pv coordinate in the J2K frame for the two scenarios (Euclidian norm of pv coordinate of scenario ② ー Euclidian norm of pv coordinate of scenario ①).
As shown below, after the transition after 1 hour of propagation, the ② scenario show different propagation compared to ① scenario.

Why does this difference occur?

In the future, I would like to use a detector to stop propagation and then resume propagation using a new detector.
For example:

propagate until altitude = 580km → stop using detector (plan manoeuvre) → propagate until mean argument of latitude of 4rad → stop and execute manoeuvre

If this sequence of stopping the propagation and then resuming cannot be used, how can two detectors be connected (as in after one detector is detected, do some calculation, then transition to another detector and propagate)?

[Vincent]

Hello @metamon,

This is most probably due to your integrator configuration, could you send us a snippet of you integrator configuration ?

I think it happens because you stop propagating between two integration steps.

Cheers.
Vincent

[metamon]

Thank you for the comment!

This is my integrator’s configuration.

self._integrator = DormandPrince853Integrator(
            self._min_integ_step,
            self._max_integ_step,
            JArray_double._cast(self._tol[0]),
            JArray_double._cast(self._tol[1]),
        )

where

min_integ_step: float = 0.0001,
max_integ_step: float = 1000.0,
init_integ_step: float = 1.0,
pos_tolerance: float = 0.1,

I also set fixed step handler of 60seconds to record the pv of the sattelite for every seconds (which should have no effect to this).

Sincerely,

[Vincent]

You are welcome

Try with the following configuration:

min_integ_step: float = 0.0001,
max_integ_step: float = 100.0,
init_integ_step: float = 1.0,
pos_tolerance: float = 0.001,

This should do the trick.

Note 1: Also, you can play with pos_tolerance between [0.001;1] in order to reduce computing time if you do not propagate for too long/do not need very high accuracy (e.g studies)

Note 2: You can most probably set min_integ_step to 1e-3 as well.

Cheers,
Vincent

[metamon]

Thank you for the suggestion.

Although the error is very small, there still exist some gap between 2hrs propagation and 1hr+1hr propagation.

Is this inevitable until max_integ_step and pos_tolerence are extremely small and close to 0?

Sincerely,

[Vincent]

I would have to dig deeper into the code, especially regarding integrator with variable step size to give you a complete explaination.

However, i guess you could avoid such issue using a fixed step integrator and making sure that you would stop on an integration step before restarting propagation.

Hope this helps,
Vincent

[metamon]

Variable step is used to optimize computational cost and to ensure accuracy in the propagation. Setting the integration step to a very small fixed value would increase computational cost abruptly and fixing the integration step too large would compromise accuracy. To balance this, I would like to use variable integration step.

I think if there is a way of adjusting the last integration step to end at the epoch I want, this error would be solved. Do you think this is possible or does this need research?

Sincerely,

[Serrof]

Hi there,

try using a very small value for the threshold of your detector, like 1e-10 or even less. The default one is 1e-6 for most detectors. Use withThreshold or at construction with the EventDetectionSettings.

As for your use case, you should be able to implement your workflow for impulsive maneuvers without manually stopping the propagation, by using filters on your triggering event and/or with the ImpulseProvider interface (>=13.0).

Cheers,
Romain.

[luc]

@metamon

I think if there is a way of adjusting the last integration step to end at the epoch I want, this error would be solved. Do you think this is possible or does this need research?

If the propagation stops due to the target time and not due to an event detector, then the last step should automatically be adjusted to match target time. This is done between lines 231 and 239 in the case of embedded Runge-Kutta Integrators like Dormand-Prince 8(5,3).

If the propagation stops due to an event and not due to target time, then it is possible to stop in the middle of a step and then the interpolator for this step is used. One surprising fact about interpolators in ODE algorithms is that their error may be much larger in the middle of the step than at the end (meaning error grows up from step start and then decreases again when reaching step end).

One thing I could suggest (but it is a long shot) would be to have a step handler (beware in this case it must be a variable step handler, not a fixed one) and to store the step lengths as new steps are produced. Then when you restart the propagation, use the size of the last step as the initial step size of the integrator. This would ensure the new integrator really starts in the same conditions as the last one.

[metamon]

Thank you for the discussion—it’s always a great help!

In my simulator, I have a propagation setup that includes multiple detectors (e.g., adL detectors → mean argument of latitude detectors → finite burn). I also have a separate propagation mode that stops at a specific epoch to execute a burn. The idea is to first propagate using detectors and record the epoch at which the mean argument of latitude detector triggers. Then, I use that epoch to stop propagation and perform a manoeuvre.

While reducing the maximum integration step size resolves the issue we encountered, it significantly increases computational cost, which is not ideal for our use case.

@Serrof
Would it be possible to provide an example of how to use detectors with filters to execute multiple manoeuvres without stopping the propagation? I’ve looked through the Orekit website and forums and attempted to build it myself, but I couldn’t get the code working.

Below is the sequence I used before attempting to implement filters:

# adL detector + uM detector + two manoeuvres
after_detector = orekitDetectors.adLDetector(5000, Client_TRUTH_file_path).withHandler(handlers.StopOnEvent())
myPropagator_Servicer_TRUE.propagator.addEventDetector(after_detector)
after_state = myPropagator_Servicer_TRUE.propagator.propagate(epoch0+dt.timedelta(hours=48))

orbital_period = after_state.getKeplerianPeriod()
epoch_after = after_state.getDate()
last_propagation_epoch = epoch_after

myPropagator_Servicer_TRUE.propagator.clearEventsDetectors()
dV_execution_detector = orekitDetectors.uM_Detector(4).withHandler(handlers.StopOnIncreasing())
myPropagator_Servicer_TRUE.propagator.addEventDetector(dV_execution_detector)
dV_execution_state_REF = myPropagator_Servicer_TRUE.propagator.propagate(absolutedate_to_datetime(last_propagation_epoch)+dt.timedelta(seconds=orbital_period))
epoch_dV_start_REF = dV_execution_state_REF.getDate()

manoeuvre_class.epoch_ignition = absolutedate_to_datetime(epoch_dV_start_REF)
manoeuvre_class.duration = burn_duration_dV_REF
manoeuvre_class.delta_v = dV_REF
manoeuvre_class.thrust_dir = control.getThrustDirection(dV_REF)
myPropagator_Servicer_TRUE.add_manoeuvre(manoeuvre_class)
end_state_REF = myPropagator_Servicer_TRUE.propagator.propagate(epoch_dV_start_REF.shiftedBy(burn_duration_dV_REF))
last_propagation_epoch = end_state_REF.getDate()

myPropagator_Servicer_TRUE.propagator.clearEventsDetectors()
dV_execution_detector = orekitDetectors.uM_Detector(5).withHandler(handlers.StopOnIncreasing())
myPropagator_Servicer_TRUE.propagator.addEventDetector(dV_execution_detector)
dV_execution_state_REF = myPropagator_Servicer_TRUE.propagator.propagate(absolutedate_to_datetime(last_propagation_epoch)+dt.timedelta(seconds=orbital_period))
epoch_dV_start_REF = dV_execution_state_REF.getDate

manoeuvre_class.epoch_ignition = absolutedate_to_datetime(epoch_dV_start_REF)
manoeuvre_class.duration = burn_duration_dV_REF_2
manoeuvre_class.delta_v = dV_REF_2
manoeuvre_class.thrust_dir = control.getThrustDirection(dV_REF_2)
myPropagator_Servicer_TRUE.add_manoeuvre(manoeuvre_class)
end_state_REF = myPropagator_Servicer_TRUE.propagator.propagate(epoch_dV_start_REF.shiftedBy(burn_duration_dV_REF_2))

@luc
Being able to store the integration step would be extremely useful when I need to stop propagation. While I can use DateDetector and filters to avoid stopping in most cases, there are scenarios in my simulator where stopping propagation is unavoidable—especially when transitioning between mission phases and planning manoeuvres accordingly.

Although this is in Java (and I prefer Python), I found an example on the forum that illustrates how to retrieve integration steps:

private static class TestStepHandler implements OrekitStepHandler {

        @Override
        public void handleStep(OrekitStepInterpolator interpolator) {
            AbsoluteDate t0 = interpolator.getPreviousState().getDate();
            AbsoluteDate tf = interpolator.getCurrentState().getDate();
            
            System.out.println("\n\tt0 = " + t0);
            System.out.println("\ttf = " + tf);
            System.out.println("\tdt = " + tf.durationFrom(t0));
            
            
        }

Is this the only way to retrieve integration step sizes—by comparing the dates of the previous and current states? Or is there a more direct method?

Sincerely,

[Serrof]

@metamon

While reducing the maximum integration step size resolves the issue we encountered, it significantly increases computational cost, which is not ideal for our use case.

Have you tried reducing the detector’s threshold first?

For filters on events used in impulsive maneuvers, check out this Java tutorial.

In Orekit 13.0, there is a native PropagationStepRecorder implementing OrekitStepHandler and storing the propagation steps. And yes you do need to compare the dates.

Cheers,
Romain.

[luc]

@metamon

Is this the only way to retrieve integration step sizes—by comparing the dates of the previous and current states? Or is there a more direct method?

Yes, even at Hipparchus level the step size of the integrator is available only in the integrators class hierarchy, at users level you just get an ODEStateAndDerivative that contains the previous state, current state, and allows you to interpolate an intermediate state, but there are no method to retrieve only the state.

The getPreviousState, getCurrentState, getDate, and durationFrom methods are all very fast, so you should not worry about this.

[metamon]

Have you tried reducing the detector’s threshold first?

Yes, I have changed the detector setting to set the threshold to 1e-10 and the descrptancy still exists.

I have used PropagationStepRecorder to record the date of each integration step and deduced the integration step just before the first detector’s condition is satisfied.

For withThreshold(1e-10), last integ step 3.637968e-12 sec
For default, last integ step 7.29503881296e-07 sec

After this detector’s condition is satisfied and its switched to the second detector and the propagation is resumed with initial integration step of 1 second (which is set by default).

# adL detector + uM detector + two manoeuvres
after_detector = orekitDetectors.adLDetector(5000, Client_TRUTH_file_path).withHandler(handlers.StopOnEvent()).withThreshold(1e-10)
print(after_detector.getThreshold())
myPropagator_Servicer_TRUE.propagator.addEventDetector(after_detector)
after_state = myPropagator_Servicer_TRUE.propagator.propagate(epoch0+dt.timedelta(hours=48))

orbital_period = after_state.getKeplerianPeriod()
epoch_after = after_state.getDate()
last_propagation_epoch = epoch_after

myPropagator_Servicer_TRUE.propagator.clearEventsDetectors()
dV_execution_detector = orekitDetectors.uM_Detector(4).withHandler(handlers.StopOnIncreasing()).withThreshold(1e-10)
print(dV_execution_detector.getThreshold())
myPropagator_Servicer_TRUE.propagator.addEventDetector(dV_execution_detector)
dV_execution_state_REF = myPropagator_Servicer_TRUE.propagator.propagate(absolutedate_to_datetime(last_propagation_epoch)+dt.timedelta(seconds=orbital_period))
epoch_dV_start_REF = dV_execution_state_REF.getDate()

@luc @Serrof @Vincent
To solve the issue, I have then set the initial integration step of the second propagation to match the last integration step used for the first propagation.

The descreptancy between the two scenarios (one without stopping the propagation and another stopping the propagation when one detector’s condition is satisfied) have been RESOLVED!

#adL detector + uM detector
after_detector = orekitDetectors.adLDetector(5000, Client_TRUTH_file_path).withHandler(handlers.StopOnEvent()) #.withThreshold(1e-10)
print(after_detector.getThreshold())
myPropagator_Servicer_TRUE.propagator.addEventDetector(after_detector)
after_state = myPropagator_Servicer_TRUE.propagator.propagate(epoch0+dt.timedelta(hours=48))
steps = integ_step_handler.copyStates()
last_integ_step = steps[-1].getDate().durationFrom(steps[-2].getDate())
print('last integ step',last_integ_step)

epoch_after = after_state.getDate()

# myPropagator_Servicer_TRUE.min_integ_step = 1e-12
myPropagator_Servicer_TRUE.init_integ_step = last_integ_step


orbital_period = after_state.getOrbit().getKeplerianPeriod()
last_propagation_epoch = epoch_after

myPropagator_Servicer_TRUE.propagator.clearEventsDetectors()
dV_execution_detector = orekitDetectors.uM_Detector(4).withHandler(handlers.StopOnIncreasing()) #.withThreshold(1e-10)
print(dV_execution_detector.getThreshold())
myPropagator_Servicer_TRUE.propagator.addEventDetector(dV_execution_detector)
dV_execution_state_REF = myPropagator_Servicer_TRUE.propagator.propagate(absolutedate_to_datetime(last_propagation_epoch)+dt.timedelta(seconds=orbital_period))
epoch_dV_start_REF = dV_execution_state_REF.getDate()

The list below shows the absolute date of each integration until the first detector’s condition (adL = 5000m) is satisfied.

2027-08-17T00:00:01.000Z
2027-08-17T00:00:11.000Z
2027-08-17T00:01:51.000Z
2027-08-17T00:04:42.057234974229800176Z
2027-08-17T00:04:54.662971882642352768Z
2027-08-17T00:12:38.987424993614013146Z
2027-08-17T00:20:29.3271816078242864Z
2027-08-17T00:28:24.822525075061776064Z
2027-08-17T00:35:58.036396197687736272Z
2027-08-17T00:43:27.432449206285127744Z
2027-08-17T00:50:08.535096437840365952Z
2027-08-17T00:56:42.699102822489749056Z
2027-08-17T01:04:05.703248036535114816Z
2027-08-17T01:10:28.540165458343835776Z
2027-08-17T01:10:41.315173362538189376Z
2027-08-17T01:18:16.450299052178706944Z
2027-08-17T01:26:04.660265914857518528Z
2027-08-17T01:33:39.582116291603597376Z
2027-08-17T01:35:56.920731615611657616Z
2027-08-17T01:35:56.920732345115538912Z

I have two questions:

1. By default, the minimum integration step for my propagator/integrator is set to be 0.0001. I have not touched this setting, but the last integration step observed is in the order of 1e-7. Why does this happen? If the threshold of the detector is set to default of 1e-6, does the integration step required to achieve this threshold override the default minimum integration step set? This causes an error during propagation as some integration step used is below the minimum integration set. This is not a big problem as reducing the minimum integration step from 1e-4 to 1e-6 would not change the computational cost as small integration step is already used for computation. However, I would like to know why.

2. I use two step handler simultaneously, one to record the spacecraft state every 60 seconds, and one to record every integration step.

myStepHandler = orekit_interfaces.OrekitStepHandler()
integ_step_handler = PropagationStepRecorder()
# myPropagator_Servicer_TRUE.propagator.setStepHandler(60., myStepHandler, ) # Stores the state every 60 s eventhough propagator uses variable stepsize.

Step_handler_multiplexer = StepHandlerMultiplexer()

Step_handler_multiplexer.add(integ_step_handler)
Step_handler_multiplexer.add(60., myStepHandler)

myPropagator_Servicer_TRUE.propagator.setStepHandler(Step_handler_multiplexer)

When a propagation is stopped, the PropagationStepRecorder stops recording and do not resume automatically.
This is a problem when I have 3 detectors for 3 consecutive propagations where I have to set the initial integration step of the third propagation to equal the last integration step for the second propagator.

Do I have to add a new integ_step_handler like below to resume recording?

Step_handler_multiplexer.add(integ_step_handler)
myPropagator_Servicer_TRUE.propagator.setStepHandler(Step_handler_multiplexer)

This causes a problem where the fixed step handler starts recording from the epoch where the second detector’s condition is satisfied. This is because the code below also resets the fixed step handler from the epoch at which the propagation has been stopped by a detector, which is not what I want. Ideally, I just want to resume or reconfigure this integ_step_handler WITHOUT effecting the fixed_step_handler. Is there any way around this problem?

[Serrof]

Hi,

In its current form, the PropagationStepRecorder resets everything at each call to init, which happens when a propagation starts. I will open an issue to be able to start from a know result.
In the mean time you can implement your own step handler.

Cheers,
Romain.

[metamon]

@luc @Serrof @Vincent

This is an unfortunate news.

I have realised when the code below is used, it not only change the initial integration step, but also initializes the propagator to the initial epoch.

myPropagator_Servicer_TRUE.init_integ_step = last_integ_step

The reason why we didnt see the rapid jump in the aROE plot was because the propagator initializes, and re-runs the propagation from initial condition to the next epoch.

For example, lets say first propagation starts from epoch A and end at epoch B, and the second propagation is to propagate until epoch C. What we want is second propagation starting from epoch B and end in epoch C. However, due to .init_integ_step, it initializes the propagator and starts propagation from epoch A to epoch C.

So it does not resume the propagation from where it ended but re-runs the propagation from scratch.

BEFORE

AFTER

If we cannot use .init_integ_step to change the first integration step used for second propagation, how can we adjust the integration step so it uses the last integration step used in first propagation?

Do you have any suggestion how to change the following code?

thank you always for your support!

#adL detector + uM detector
after_detector = orekitDetectors.adLDetector(5000, Client_TRUTH_file_path).withHandler(handlers.StopOnEvent()) #.withThreshold(1e-10)
print(after_detector.getThreshold())
myPropagator_Servicer_TRUE.propagator.addEventDetector(after_detector)
after_state = myPropagator_Servicer_TRUE.propagator.propagate(epoch0+dt.timedelta(hours=48))
steps = integ_step_handler.copyStates()
last_integ_step = steps[-1].getDate().durationFrom(steps[-2].getDate())
print('last integ step',last_integ_step)

epoch_after = after_state.getDate()

# myPropagator_Servicer_TRUE.min_integ_step = 1e-12
myPropagator_Servicer_TRUE.init_integ_step = last_integ_step


orbital_period = after_state.getOrbit().getKeplerianPeriod()
last_propagation_epoch = epoch_after

myPropagator_Servicer_TRUE.propagator.clearEventsDetectors()
dV_execution_detector = orekitDetectors.uM_Detector(4).withHandler(handlers.StopOnIncreasing()) #.withThreshold(1e-10)
print(dV_execution_detector.getThreshold())
myPropagator_Servicer_TRUE.propagator.addEventDetector(dV_execution_detector)
dV_execution_state_REF = myPropagator_Servicer_TRUE.propagator.propagate(absolutedate_to_datetime(last_propagation_epoch)+dt.timedelta(seconds=orbital_period))
epoch_dV_start_REF = dV_execution_state_REF.getDate()

[Serrof]

Hi there,

the method you’re talking about on the initial step size is your own wrapping of Orekit, we don’t know what’s inside.
In NumericalPropagator though, you can enable or disable the option to erase the initial time by the last propagation time with setResetAtEnd

Cheers,
Romain