Hi,
I have used the Orekit Python Wrapper in order to create a small simulation bench, and then started by implementing basic measurement simulation and orbit estimation with the BatchLeastSquareEstimator.
I wanted to get a little further an try to implement maneuver detection / estimation in the OD process. In fact I wanted to try out a method which differs from what is suggested in the Orekit tutorials on the same topic. The idea is to consider constant by segment empirical forces in the estimation, which seems to be wonderfully modeled by the class TimeSpanParametricAcceleration.java, combined with PolynomialAccelerationModel.java. These empirical forces should be close to zero most of the time, and raise whenever a maneuver has been done.
So I implemented this and tried to apply it on “no maneuver case” scenarios, and get these kind of exceptions :
- org.orekit.errors.OrekitException: altitude (71 662,495 m) au dessous de la limite autorisée de 120 000 m
- org.orekit.errors.OrekitException: impossible de calculer l’anomalie excentrique hyperbolique à partir de l’anomalie moyenne après 1 000 000 itérations
Anyway, it turns out that the empirical forces are reaching tremendous values, thus making completely aberrant dynamics at the cost of the orbit realism. (when looking at it, the values for empirical forces raise up to 1e7 N !!)
Have you got any idea of where the problem could be? I don’t think it comes from the measurements because the estimation is doing very well in the “classic” OD process without empirical forces. In my opinion it is probably a bad understanding/implementation of these empirical forces on my side, or maybe this maneuver estimation method is just not meant to succeed?
Here are extracts of the scripts I wrote :
Method to add empirical forces in by "NumericalPropagatorBuilder wrapper"
def add_empirical_forces(self, time_span: List[datetime], time_step: float):
"""
Add empirical forces to the force model. The corresponding forces are set to 0 in the propagator in itself, but
can be tuned in an orbit determination process to match observation in order to model a maneuver for example.
In its current implementation, the empirical forces will be modelled with a constant acceleration on a segment.
There use is not recommended in a classical restitution without maneuvers, since it allows to adjust a force
model to the not imperfect observation, at the cost of realism.
:param time_span: Time window on which the empirical forces should be added, typically the measurement arc
:param time_step: Time step in seconds, corresponding to constant acceleration segment length
:return:
"""
self.consider_empirical_forces = True
# segment time window
start_date = time_span[0]
stop_date = time_span[1]
dt = timedelta(seconds=time_step)
# list to store parameters identifiers : will be needed for estimator settings
identifiers = []
# build empirical forces for the 3 axes, with a polynomial with degree zero (i.e. constant)
param_prefix_x = f"empirical_segment_1_x"
force_model_x = PolynomialAccelerationModel(param_prefix_x, datetime_to_absolute_date(start_date), 0)
empirical_force_x = TimeSpanParametricAcceleration(Vector3D.PLUS_I, False, force_model_x)
param_prefix_y = f"empirical_segment_1_y"
force_model_y = PolynomialAccelerationModel(param_prefix_y, datetime_to_absolute_date(start_date), 0)
empirical_force_y = TimeSpanParametricAcceleration(Vector3D.PLUS_J, False, force_model_y)
param_prefix_z = f"empirical_segment_1_z"
force_model_z = PolynomialAccelerationModel(param_prefix_z, datetime_to_absolute_date(start_date), 0)
empirical_force_z = TimeSpanParametricAcceleration(Vector3D.PLUS_K, False, force_model_z)
# store parameter identifier, the syntax "[0]" is added upon force model creation
identifiers.extend([f"{param_prefix_x}[0]", f"{param_prefix_y}[0]", f"{param_prefix_z}[0]"])
current_date = start_date + dt
count = 2
while current_date < stop_date:
# add empirical forces segment (segment are [t, inf[ and override old ones on overlapping periods)
param_prefix_x = f"empirical_segment_{count}_x"
force_model_x = PolynomialAccelerationModel(param_prefix_x, datetime_to_absolute_date(current_date), 0)
param_prefix_y = f"empirical_segment_{count}_y"
force_model_y = PolynomialAccelerationModel(param_prefix_y, datetime_to_absolute_date(current_date), 0)
param_prefix_z = f"empirical_segment_{count}_z"
force_model_z = PolynomialAccelerationModel(param_prefix_z, datetime_to_absolute_date(current_date), 0)
empirical_force_x.addAccelerationModelValidAfter(force_model_x, datetime_to_absolute_date(current_date))
empirical_force_y.addAccelerationModelValidAfter(force_model_y, datetime_to_absolute_date(current_date))
empirical_force_z.addAccelerationModelValidAfter(force_model_z, datetime_to_absolute_date(current_date))
# store parameter identifier, the syntax "[0]" is added upon force model creation
identifiers.extend([f"{param_prefix_x}[0]", f"{param_prefix_y}[0]", f"{param_prefix_z}[0]"])
# updates for next iteration
current_date += dt
count += 1
# add empirical forces to the propagator
self.orekit_builder.addForceModel(empirical_force_x)
self.orekit_builder.addForceModel(empirical_force_y)
self.orekit_builder.addForceModel(empirical_force_z)
return identifiers
This method is then called in a method from my estimator :
Add the empirical forces to the estimated parameters
def set_empirical_forces_estimation(self, time_step: float = 3600.):
"""
Activates empirical forces estimation, in order to detect maneuvers for example.
:param time_step: Segment length in seconds on which empirical forces are defined, default is 1 hour.
:return:
"""
if len(self.measure_list) == 0:
print(f"Measures should be added before setting this option.")
else:
# measurement list is chronologically sorted
time_window = [absolute_date_to_datetime(self.measure_list[0].getDate()),
absolute_date_to_datetime(self.measure_list[-1].getDate())]
# add empirical forces to the propagator
param_identifiers = self.propagator_builder.add_empirical_forces(time_window, time_step)
# add empirical forces to the estimated parameters
for param_id in param_identifiers:
param = self.propagator_builder.orekit_builder.getPropagationParametersDrivers().findByName(param_id)
param.setSelected(True)
Any kind of idea or explanation would be of great help !
Thank you very much ! 
