# Maneuver estimation with empirical forces

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)

# 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 += dt
count += 1

# add empirical forces to the propagator

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

# 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 !

Hi @garyy,

Welcome to the forum!

Maybe your model has too many degrees of freedom thatâ€™s why itâ€™s reaching aberrant dynamics.
Other than that I would not know.
What type of measurements are you using? Are you estimating anything else apart from the orbit and the maneuvers?
I havenâ€™t seen anything wrong in your snippets of code from a quick reading.
Could you share with us a runnable program (with simulated PV measurements for example)?

Cheers,
Maxime

Hi there,

Iâ€™m not familiar with those classes but I assume you can give bounds to the estimator. Canâ€™t you use that to avoid non physical accelerations?

Cheers,
Romain.

Hi @MaximeJ,

Thank you for your answer. I used Range measurements, and wasnâ€™t estimating anything but the orbital elements and the empirical forces.
In fact I guess youâ€™re right, the model has just way too many degrees of freedom. For instance with 1 hour segments over a day, the additional number of parameters reaches 24 x 3 = 72. Considering the fact that the time windows will be probably way larger in real life applications, this method may be doomed to fail.
Anyway if youâ€™re still interested in the topic Iâ€™m joining a main script on which I gathered all the components for such an estimation, using PVT measurements.

Cheers,
Yvan

Main script for empirical forces estimation
``````import orekit

vm = orekit.initVM()

from orekit.pyhelpers import setup_orekit_curdir

setup_orekit_curdir("../orekit-data/orekit-data.zip")

from org.orekit.orbits import KeplerianOrbit, PositionAngleType
from org.orekit.estimation.measurements import PV, ObservableSatellite, EstimatedMeasurement
from org.orekit.frames import FramesFactory
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import Constants, IERSConventions
from org.orekit.bodies import OneAxisEllipsoid, CelestialBodyFactory
from org.orekit.propagation.conversion import NumericalPropagatorBuilder, DormandPrince853IntegratorBuilder
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel, ThirdBodyAttraction
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.drag import DragForce, IsotropicDrag
from org.orekit.forces.empirical import TimeSpanParametricAcceleration, PolynomialAccelerationModel
from org.orekit.models.earth.atmosphere import DTM2000
from org.orekit.models.earth.atmosphere.data import MarshallSolarActivityFutureEstimation
from org.hipparchus.linear import QRDecomposer
from org.hipparchus.optim.nonlinear.vector.leastsquares import GaussNewtonOptimizer
from org.orekit.estimation.leastsquares import BatchLSEstimator

from org.hipparchus.geometry.euclidean.threed import Vector3D

# propagation parameters for LEO objects
PROPAGATION_PARAMS_LEO = {
"name": "LEO",
"force_model": {
"earth_potential": 40,
"sun_attraction": True,
"moon_attraction": True,
"srp": False,
"drag": True
},
"satellite_model": {
"cross_section": 0.10,
"cr": 3e-2,
"cd": 0.025
}
}

# default estimation settings
DEFAULT_ESTIMATION_PARAMS = {
"convergence_threshold": 1e-3,
"max_iterations": 25,
"max_evaluations": 35
}

# start and stop date for the study
T_SPAN = 3.0 * 86400.
START_DATE = AbsoluteDate(2024, 1, 1, 0, 0, 0., TimeScalesFactory.getUTC())
STOP_DATE = START_DATE.shiftedBy(T_SPAN)

def create_propagator_builder(orbit, propagation_params):
"""
Create a propagator builder with the given propagator parameters.
:param orbit: orbit from which the propagator initial state will be computed
:param propagation_params: Dictionary containing the necessary elements to build a propagator builder.
:return: A NumericalPropagatorBuilder from orekit
"""
# integrator settings
min_step = 0.001
max_step = 100.0
position_error = 1.0

# integrator builder
integrator_builder = DormandPrince853IntegratorBuilder(min_step, max_step, position_error)

# propagator initialization
prop_builder = NumericalPropagatorBuilder(orbit, integrator_builder, PositionAngleType.TRUE, 1.0)

# force model
# earth gravity
Constants.WGS84_EARTH_FLATTENING,
FramesFactory.getITRF(IERSConventions.IERS_2010, True))
gravity_provider = GravityFieldFactory.getNormalizedProvider(propagation_params["force_model"]["earth_potential"],
propagation_params["force_model"]["earth_potential"])

# third body attraction : moon
if propagation_params["force_model"]["moon_attraction"]:

# third body attraction : sun
if propagation_params["force_model"]["sun_attraction"]:

if propagation_params["force_model"]["srp"]:
propagation_params["satellite_model"]["cr"])

# atmospheric drag
if propagation_params["force_model"]["drag"]:
msafe = MarshallSolarActivityFutureEstimation(MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,
MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)
atmosphere = DTM2000(msafe, CelestialBodyFactory.getSun(), earth)
isotropic_drag = IsotropicDrag(propagation_params["satellite_model"]["cross_section"],
propagation_params["satellite_model"]["cd"])
drag_model = DragForce(atmosphere, isotropic_drag)

return prop_builder

"""
Add empirical forces to an existing propagator builder.
:param prop_builder: Propagator builder on which empirical forces should be added
:param time_window: Time window on which empirical forces should be considered
:return:
"""
# segment time window
start_date = time_window[0]
stop_date = time_window[1]
t_step = 3600.

# 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, 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, 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, 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.shiftedBy(t_step)
count = 2
while current_date.compareTo(stop_date) <= 0:
# 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, current_date, 0)

param_prefix_y = f"empirical_segment_{count}_y"
force_model_y = PolynomialAccelerationModel(param_prefix_y, current_date, 0)

param_prefix_z = f"empirical_segment_{count}_z"
force_model_z = PolynomialAccelerationModel(param_prefix_z, current_date, 0)

# 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 = current_date.shiftedBy(t_step)
count += 1

# add empirical forces to the propagator

return identifiers

def create_ls_estimator(initial_orbit, prop_params, estimation_params, estimate_empirical_forces=False):
"""
Create an estimator with the desired settings.
:param initial_orbit: Initial orbit for the estimation
:param prop_params: Propagation params used for propagator builder generation
:param estimation_params: Estimation parameters
:param estimate_empirical_forces : If True, empirical forces will be estimated as well as orbital elements
:return: A BatchLSEstimator from orekit
"""
# create a specific propagator builder
prop_builder = create_propagator_builder(initial_orbit, prop_params)

# account for empirical forces if asked
if estimate_empirical_forces:

for param_id in ids:
param = prop_builder.getPropagationParametersDrivers().findByName(param_id)
param.setSelected(True)
param.setScale(1.0)

# create estimator
optimizer = GaussNewtonOptimizer(QRDecomposer(1e-20), False)
estimator = BatchLSEstimator(optimizer, prop_builder)

# estimator settings
estimator.setParametersConvergenceThreshold(estimation_params["convergence_threshold"])
estimator.setMaxIterations(estimation_params["max_iterations"])
estimator.setMaxEvaluations(estimation_params["max_evaluations"])

return estimator

if __name__ == '__main__':
# build reference orbit
reference_date = START_DATE
reference_orbit = KeplerianOrbit(8000. * 1e3, 0.01, radians(30.), 0., 0., 0.,
PositionAngleType.TRUE,
FramesFactory.getGCRF(),
reference_date,
Constants.WGS84_EARTH_MU)

# create a propagator
propagator_builder = create_propagator_builder(reference_orbit, PROPAGATION_PARAMS_LEO)
leo_propagator = propagator_builder.buildPropagator(propagator_builder.getSelectedNormalizedParameters())

# build pvt measurement from the reference orbit (every hour over 3 days)
print("Generating PVT measurements...")
std_pos = 1.0   # m
std_vel = 1e-3  # m/s

end_date = STOP_DATE
dt = 3600.

generator = leo_propagator.getEphemerisGenerator()
leo_propagator.propagate(end_date)
ephemeris = generator.getGeneratedEphemeris()

measurement_list = []
current_date = reference_date

while current_date.compareTo(end_date) <= 0.:
pv_coordinates = ephemeris.getPVCoordinates(current_date, FramesFactory.getGCRF())
pv_measurement = PV(current_date, pv_coordinates.getPosition(), pv_coordinates.getVelocity(), std_pos, std_vel,
1., ObservableSatellite(0))

measurement_list.append(pv_measurement)
current_date = current_date.shiftedBy(dt)

print(f"Done. Number of measures generated : {len(measurement_list)}")

# build a classic estimator and add the pv measurements
PositionAngleType.TRUE,
FramesFactory.getGCRF(),
reference_date,
Constants.WGS84_EARTH_MU)

classic_ls_estimator = create_ls_estimator(noisy_orbit, PROPAGATION_PARAMS_LEO, DEFAULT_ESTIMATION_PARAMS)
for measurement in measurement_list:

# run a classic orbit determination to assess convergence in this case
print("Launching Classical Orbit Determination...")
classic_ls_estimator.estimate()
print("Done.")

# build an estimator accounting for empirical forces and add the pv measurements
ls_estimator = create_ls_estimator(noisy_orbit, PROPAGATION_PARAMS_LEO, DEFAULT_ESTIMATION_PARAMS, True)
for measurement in measurement_list:

print("Launching Orbit Determination with Empirical Forces...")
ls_estimator.estimate()
print(ls_estimator.getPropagatorParametersDrivers(True).getDrivers())

print("Done.")

``````
1 Like

Hi @Serrof,

Thank you for your answer as well. Indeed I tried to set bounds to the parameters corresponding to the empirical forces with small values, such as [-1e8, 1e8]. I used the setMinValue() and the setMaxValues() functions alongside the setScale() function in order to do so. Whatâ€™s really curious is that the values seems to remain in the specified bounds (tho reaching them each time), but errors such as org.orekit.errors.OrekitException: impossible de calculer lâ€™anomalie excentrique hyperbolique Ă  partir de lâ€™anomalie moyenne aprĂ¨s 1 000 000 itĂ©rations are raised anyway.

Iâ€™m not sure such whether small perturbations like this can have such influence on the orbit, so maybe something strange is happening here.

Cheers,
Yvan

Hi @garyy,

Thanks for the code.
It works if the number of empirical forces is reduced. I tried with `time_window = [START_DATE.shiftedBy(0.), START_DATE.shiftedBy(24. * 3600.)]` in `add_empirical_forces` and it converged with maximum empirical forces lower than 1e-13 N.
So that rules out the fact that there is a mistake somewhere in the code.

Iâ€™ve tried on 2 or 3 days with more measurements (one every 10 minutes) to try to overdetermine the problem but I ended up with a failure (negative eccentricity).

I suggest you set up an observer to follow the convergence during the batch LS estimator (here is a Java example) because the OD is very long I find it fishy. I found the OD is very long and the observer will show you what happens after each iteration, itâ€™s very handy.

Maybe you could activate this on a selected number of segments when an OD fails or when the orbit from one day/pass to another has changed so much that you suspect a maneuver occurred.

Hope this helps,
Maxime

PS:

I guess you meant 1e-8 here because 10^{8} is a very large value for empirical forces!

Hi @MaximeJ ,

Indeed convergence is reached with realistic values for empirical forces when the number of segments is decreasing. Thank you for the tip of the observer, Iâ€™ll try it out.

Maybe you could activate this on a selected number of segments when an OD fails or when the orbit from one day/pass to another has changed so much that you suspect a maneuver occurred.

That could be a solution, Iâ€™ll think about it !

I guess you meant 1e-8 here because 10^8 is a very large value for empirical forces!

Yep sorry for the typo I used the bounds [-1e-8, 1e-8]
By the way this may highlight another issue, because even with the lowest possible bounds, the dynamics seems to be broken as well in the end when the number of segments is too large.