The BatchLSEstimator doesn't reach convergence

Good evening
I’m a space engineering student and I’m doing a thesis on Orbit Determination and Initial Orbit Determination on Molniya satellites using python. I use a telescope for RA and Dec measurements, with a pixel scale = 4e-5 arcsec.
I’m trying to use the batchLSEstimator starting from an initial orbit determination that I have done with Laplace method, to obtain the inital position and velocity of the object.
I use the LevenbergMarquardtOptimizer with a set of measurements of RA and DEC, given at the Batch in RaDec form as it want. The problem is that it never reach convergence but with a threshold of 1e5, if I use smaller threshold, at the estimator.estimate() point, it give me an error that request a smaller minimum step (“prop_min_step” in the code), so I tried with a smaller step but it give me the same error.

The code is this:

estimator_position_scale = 10.0 # [m]
prop_min_step = 1e-9 # [s]
prop_max_step = 20.0 # [s]
prop_position_error = 100.0 # [m]

integrator_Builder = DormandPrince853IntegratorBuilder(prop_min_step, prop_max_step, prop_position_error)
propagatorBuilder = init_numerical_propagator(initial_guess_orbit, integrator_Builder, estimator_position_scale, sat_list, j2000, Body, obsData_half)

Observatory site

lon = radians(13.37503)
lat = radians(41.76524)
alt = 555.0

Initial covariance measurementes (arcsec)

sigma = 4e-5

Batch Least Square estimator

estimator_convergence_thres = 1e-2
estimator_max_iterations = 500
estimator_max_evaluations = 800

optimizer = LevenbergMarquardtOptimizer()

estimator = BatchLSEstimator(optimizer, propagatorBuilder)

estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
estimator.setMaxIterations(estimator_max_iterations)
estimator.setMaxEvaluations(estimator_max_evaluations)

Add measurements

for radec in raDecs:
estimator.addMeasurement(radec)

Observer

bls_observer = CustomBatchLSObserver()
estimator.setObserver(bls_observer)

estimator.getOrbitalParametersDrivers(True)
estimator.getPropagatorParametersDrivers(True)
estimator.getMeasurementsParametersDrivers(True)

batch_orbit = estimator.estimate()

#############################

These are some of the measurements that I use:

ANGLE_1 = 2023-05-30T19:48:29.084 225.958399
ANGLE_2 = 2023-05-30T19:48:29.084 -26.342211
ANGLE_1 = 2023-05-30T19:48:31.126 225.966465
ANGLE_2 = 2023-05-30T19:48:31.126 -26.362219
ANGLE_1 = 2023-05-30T19:48:33.168 225.974163
ANGLE_2 = 2023-05-30T19:48:33.168 -26.378820
ANGLE_1 = 2023-05-30T19:48:35.210 225.983265
ANGLE_2 = 2023-05-30T19:48:35.210 -26.396772
ANGLE_1 = 2023-05-30T19:48:37.252 225.992470
ANGLE_2 = 2023-05-30T19:48:37.252 -26.416811
ANGLE_1 = 2023-05-30T19:48:39.294 226.003625
ANGLE_2 = 2023-05-30T19:48:39.294 -26.437235
ANGLE_1 = 2023-05-30T19:48:41.336 226.011380
ANGLE_2 = 2023-05-30T19:48:41.336 -26.454963
ANGLE_1 = 2023-05-30T19:48:43.378 226.020752
ANGLE_2 = 2023-05-30T19:48:43.378 -26.473826
ANGLE_1 = 2023-05-30T19:48:45.420 226.030952
ANGLE_2 = 2023-05-30T19:48:45.420 -26.492258
ANGLE_1 = 2023-05-30T19:48:47.462 226.038145
ANGLE_2 = 2023-05-30T19:48:47.462 -26.510381
ANGLE_1 = 2023-05-30T19:48:49.504 226.048209
ANGLE_2 = 2023-05-30T19:48:49.504 -26.529082
ANGLE_1 = 2023-05-30T19:48:51.546 226.055752
ANGLE_2 = 2023-05-30T19:48:51.546 -26.546325
ANGLE_1 = 2023-05-30T19:48:53.588 226.068137
ANGLE_2 = 2023-05-30T19:48:53.588 -26.567976
ANGLE_1 = 2023-05-30T19:48:55.630 226.076027
ANGLE_2 = 2023-05-30T19:48:55.630 -26.584372
ANGLE_1 = 2023-05-30T19:48:57.672 226.085744
ANGLE_2 = 2023-05-30T19:48:57.672 -26.604844
ANGLE_1 = 2023-05-30T19:48:59.714 226.094550
ANGLE_2 = 2023-05-30T19:48:59.714 -26.622727
ANGLE_1 = 2023-05-30T19:49:01.756 226.103223
ANGLE_2 = 2023-05-30T19:49:01.756 -26.641743

With
ANGLE_1: Right Ascension
ANGLE_2: Declination

The angles are in EME2000 frame

TLEs of the object are:

2023-05-30 13:25:18

tleLine1 = ‘1 22671U 93035A 23150.55923760 .00000178 00000-0 39690-3 0 9992’
tleLine2 = ‘2 22671 63.1320 34.9861 4846547 12.5651 356.2237 5.62935234371833’

2023-05-30 21:56:54

l1_next = ‘1 22671U 93035A 23150.91451846 .00000172 00000-0 38251-3 0 9999’
l2_next = ‘2 22671 63.1322 34.7799 4846541 12.5698 356.2239 5.62935391371856’

I know that TLEs are mean parameters and not osculating parameters, but I use the TLEs to compare the orbit obtained with Batch and the orbit of spacetrack. The most different parameters between these two orbit is the argument of perigee (~100 degrees) and the semi-axis major (difference of 2000 km). I think that is threshold fault, 1e5 is too high and the Batch doesn’t work.

Do you have any advice for me to solve the problem? Or do you recommend using other methods?
Thank you in advance for the help

Hi there,

Could you show us how you create the measurements objects for the OD? All your elevations are negative so there’s something wrong here probably.
Also it’s worth saying that for “usual” Earth orbits used by man-made objects, Laplace IOD is not really recommended. Use Gauss instead. It’s available in Orekit 12.

Best,
Romain.