[Python Wrapper] BatchLSEstimator not initializing

Orekit Python Wrapper usage
1

Hello All,

I am a new user of the Orekit library, and doing my best to figure it out from the Java documentation and the few example scripts out there. I’m using the Batch Least Squares Estimator with Azimuth and Elevation angles. I’m working off of this script as an example for how to set up a BatchLSEstimator:


The script uses range values, and the data that I have available is in Azimuth and Elevation. It seems to me that the setup would be the exact same, and instead of using the Range() measurement object, use the AngularAzEl(). I’m able to initialize the addmeasurement() for the estimator, but the error comes when I get to the line estimator.estimate(). I know Python well enough to diagnose errors there, but I don’t know Java at all, so anyone who recognizes the error that could help is greatly appreciated!

Here’s my code that deals with the propagator and estimator set up.

#Orbit Propagator Parameters
prop_min_step = 0.001 
prop_max_step = 300.0 
prop_position_error = 10.0 

#Estimator Parameters
estimator_position_scale = 1.0 
estimator_convergence_thres = 1e-3
estimator_max_iterations = 25
estimator_max_evaluations = 35

#Measurement Parameters
azError = 0.001
elError = 0.001
azBaseWeight = 1.0
elBaseWeight = 1.0

#Defining Integrator
integratorBuilder = DormandPrince853IntegratorBuilder(prop_min_step, 
prop_max_step, prop_position_error)
propagatorBuilder = NumericalPropagatorBuilder(keplerianOrbit,
                                           integratorBuilder, PositionAngle.MEAN, 
estimator_position_scale)
nadirPointing = NadirPointing(eci, wgs84Ellipsoid)
propagatorBuilder.setMass(400.0)
propagatorBuilder.setAttitudeProvider(nadirPointing)

########## ADDING PERTURBATIONS

#Earth Gravity Field
gravityProvider = GravityFieldFactory.getConstantNormalizedProvider(64, 64)
gravityAttractionModel = HolmesFeatherstoneAttractionModel(ecef, gravityProvider)
propagatorBuilder.addForceModel(gravityAttractionModel)

#3rd Body
moon_3dbodyattraction = ThirdBodyAttraction(moon)
propagatorBuilder.addForceModel(moon_3dbodyattraction)
sun_3dbodyattraction = ThirdBodyAttraction(sun)
propagatorBuilder.addForceModel(sun_3dbodyattraction)

########## BATCH LEAST SQUARES ESTIMATOR SETUP

#Defining Optimizer
matrixDecomposer = QRDecomposer(1e-11)
optimizer = GaussNewtonOptimizer(matrixDecomposer, False)

#Defining Estimator
estimator = BatchLSEstimator(optimizer, propagatorBuilder)
estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
estimator.setMaxIterations(estimator_max_iterations)
estimator.setMaxEvaluations(estimator_max_evaluations)

########## ADDING MEASUREMENTS

#Add Measurements to Estimator
for j in range(0,len(az)):
    orekitAzEl = AngularAzEl(GroundStation(stationFrame),
					datetime_to_absolutedate(datetime[j]),
					JArray('double')([radians(az[j]),radians(el[j])]),
                    JArray('double')([radians(azError),radians(elError)]),
                    JArray('double')([azBaseWeight,elBaseWeight]),
                    ObservableSatellite(1))
    estimator.addMeasurement(orekitAzEl)

########## STATISTICAL ORBIT DETERMINATION

estimatedPropagatorArray = estimator.estimate()

Here’s the error I’m getting:

Traceback (most recent call last):
File "EstimationLockheed.py", line 279, in <module>
estimatedPropagatorArray = estimator.estimate()
orekit.JavaError: <super: <class 'JavaError'>, <JavaError object>>
Java stacktrace:
java.lang.ArrayIndexOutOfBoundsException: 1
at org.orekit.estimation.measurements.AngularAzEl.theoreticalEvaluation(AngularAzEl.java:88)
at org.orekit.estimation.measurements.AbstractMeasurement.estimate(AbstractMeasurement.java:204)
at org.orekit.estimation.leastsquares.MeasurementHandler.handleStep(MeasurementHandler.java:95)
at org.orekit.propagation.PropagatorsParallelizer$SinglePropagatorHandler.handleStep(PropagatorsParallelizer.java:322)
at org.orekit.propagation.integration.AbstractIntegratedPropagator$AdaptedStepHandler.handleStep(AbstractIntegratedPropagator.java:856)
at org.hipparchus.ode.AbstractIntegrator.acceptStep(AbstractIntegrator.java:429)
at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:290)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:466)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:414)
at org.orekit.propagation.PropagatorsParallelizer.propagate(PropagatorsParallelizer.java:141)
at org.orekit.estimation.leastsquares.BatchLSModel.value(BatchLSModel.java:245)
at org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresFactory$LocalLeastSquaresProblem.evaluate(LeastSquaresFactory.java:440)
at org.orekit.estimation.leastsquares.BatchLSEstimator$TappedLSProblem.evaluate(BatchLSEstimator.java:602)
at org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.optimize(GaussNewtonOptimizer.java:399)
at org.orekit.estimation.leastsquares.BatchLSEstimator.estimate(BatchLSEstimator.java:422)

Any help diagnosing the problem is greatly appreciated! If you need more explanation I’m glad to help!

2

Hi @JPJohn_2015,

Welcome to the Orekit forum !!

I haven’t tried your code so I cannot be 100% sure but I think your error comes from the way you create your measurements.

You should set ObservableSatellite(0) instead of ObservableSatellite(1) when calling the constructor of AngularAzEl.

That should fix your ArrayIndexOutOfBoundsException bug since satellites are ordered from 0 to n in Orekit’s multi-satellite orbit determination algorithm, not from 1 to n.

Hope this helps,
Maxime

3

Hi @JPJohn_2015,

Like @MaximeJ said, the ObservableSatellite starts from 0. I upgraded the orbit determination example to Orekit 10 (and fixed several issues) a few days ago by the way, you can see the ObservableSatellite(0) parameter here: https://nbviewer.jupyter.org/github/GorgiAstro/laser-orbit-determination/blob/6f0ac2376d31a52bd28acf37a56d96cafd3d3764/02-orbit-determination-example.ipynb

Other than that, your code with AngularAzEl observations looks fine. My example originally included angular measurements from laser ground stations too, but because they actually degraded the quality of the orbit determination (I did not weigh them properly), I removed them. I will add them again to the example.

Cheers
Clément

4

@MaximeJ and @yzokras,

Thank you for the quick response! Yes that was the issue. A followup if you are willing to help out, I was able to get the estimator running, but it quickly errored out. I am adding measurements from a pass which is about a day later, so that was my thought from the True Anomaly issue, so is there anything I need to do differently for these types of problems compared to semi-continuous data? I’m observing LEO objects, so the passes are few and short.

org.orekit.errors.OrekitIllegalArgumentException: true anomaly 3.142 out of hyperbolic range (e = 1, -3.142 < v < 3.142)
5

Hi @JPJohn_2015,

Did you solve your issue ?

If not, what type of orbit are you using in input ? Is it hyperbolic to start with ?

I think we’ll need more data of you want us to help you out on this one.

Regards,
Maxime

6

Hi @MaximeJ

Still stuck on this. I am able to process data that is all from the same pass fine, however, I cannot process data from a second pass that is right afterwards. These objects are LEO, specifically NORAD 877, so I know it is basically circular, inclination at 65ish degrees, I’m not sure how it is going hyperbolic. I am using the Gooding IOD module for the initial guess of the state required by the propagator. My thought is that since it is from a different pass, then the object has completed a full revolution, and the true anomaly might be tracked incorrectly or something.

My first set of measurements start at JD 2458413.60081039
The next pass of measurements is at JD 2458414.03186308

The observations are 10 hours apart, so the object would have made about 6 or 7 revolutions before the next set of observations.

My question I guess is the BatchLSEstimator able to process measurements that are that far apart, or should I switch to a Kalman Filter type system?

7

Hi @JPJohn_2015,

I think the BatchLSEstimator should be fine.
Can you share the code that is bugging so we may look into it ?

I don’t fully understand your setup, could you please elaborate a little on this ?

You’re using Gooding to estimate the initial guess at the beginning of the first pass, the second one or both ?

If your first pass estimates the orbit well have you tried propagating from this estimated orbit to the beginning of your second pass and use the propagated orbit as an input for your second batch LS estimation ?

Maxime

8

Hi @MaximeJ,

Here is my code. I cannot give an of the data that I’m using as there is an NDA associated with that, but I can walk you through my code and what I’m doing.

First step is to define all the coordinate frames, celestial bodies for perturbations and reference frames, then the time coordinate.

########## COORDINATE FRAMES, CELESTIAL BODIES AND TIME

#Ground Station Position
latitude = XX.XXXX
longitude = XXX.XXXX
altitude = XXX.XX

#Defining ECI AND ECEF Coordinate Frames
gcrf = FramesFactory.getGCRF()
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
eci = gcrf
ecef = itrf

#Defining Topocentric Coordinate Frame
tod = FramesFactory.getTOD(IERSConventions.IERS_2010, False)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 
                     Constants.WGS84_EARTH_FLATTENING, 
                     itrf)
station = GeodeticPoint(radians(latitude), 
                    radians(longitude), 
                    altitude)
stationFrame = TopocentricFrame(earth, station, "Esrange")

#Defining Celestial Bodies
wgs84Ellipsoid = ReferenceEllipsoid.getWgs84(ecef)
moon = CelestialBodyFactory.getMoon()
sun = CelestialBodyFactory.getSun()

#Defining Time System
utc = TimeScalesFactory.getUTC()

Next I go through the Initial orbit determination process. I am using the gooding method that takes in three measurements and outputs an Keplerian Orbit object. Gooding requires a range guess, so I have a function that estimates the range based on orbital altitude. The IOD results are matching up very well with published results. Inclination is within 0.5 Degrees, Right Ascension is pretty spot on, and eccentricity and semi-major axis are within range such that statistical orbit determination would improve the estimate.

########## INITIAL ORBIT DETERMINATION

#Best Range Guess
rho1, rho3 = RangeGuess(el[0],el[2],(orbitHeight[0]+idx))

#Gooding Initial Orbit Determination
iod = IodGooding(eci, Constants.WGS84_EARTH_MU)
keplerianOrbit = iod.estimate(rSiteList[0],rSiteList[1],rSiteList[2],
        slantRange[0],datetime_to_absolutedate(datetime[0]),
        slantRange[1],datetime_to_absolutedate(datetime[1]),
        slantRange[2],datetime_to_absolutedate(datetime[2]),
        rho1,rho3)

Next is the propagator setup. Here are the constant values that are used to initialize the propagator. I’m using the same values from the example code earlier for propagator settings, and the azimuth and elevation angle errors and base errors are rough estimates.

########## PROPAGATOR SETUP

#Orbit Propagator Parameters
prop_min_step = 0.001 
prop_max_step = 300.0 
prop_position_error = 10.0 

#Estimator Parameters
estimator_position_scale = 1.0 
estimator_convergence_thres = 1e-3
estimator_max_iterations = 5
estimator_max_evaluations = 35

#Measurement Parameters
azError = 0.001
elError = 0.001
azBaseWeight = 1/azError
elBaseWeight = 1/elError

#Defining Integrator
integratorBuilder = DormandPrince853IntegratorBuilder(prop_min_step, 
prop_max_step, prop_position_error)
propagatorBuilder = NumericalPropagatorBuilder(keplerianOrbit,
                                           integratorBuilder, PositionAngle.MEAN, 
estimator_position_scale)
nadirPointing = NadirPointing(eci, wgs84Ellipsoid)
propagatorBuilder.setMass(400.0)
propagatorBuilder.setAttitudeProvider(nadirPointing)

Next is adding perturbations to the propagator. Currently only have Earth gravity field and n-body perturbations active because they do not require information about spacecraft size, mass and attitude. I understand drag is a large factor for low earth orbiting satellites, but I don’t think that in the span of 10 hours it would be large enough to cause errors. Over longer periods it needs to be accounted for, but for less than a day, I don’t think it should be a problem.

########## ADDING PERTURBATIONS

#Earth Gravity Field
gravityProvider = GravityFieldFactory.getConstantNormalizedProvider(64, 64)
gravityAttractionModel = HolmesFeatherstoneAttractionModel(ecef, gravityProvider)
propagatorBuilder.addForceModel(gravityAttractionModel)

#3rd Body
moon_3dbodyattraction = ThirdBodyAttraction(moon)
propagatorBuilder.addForceModel(moon_3dbodyattraction)
sun_3dbodyattraction = ThirdBodyAttraction(sun)
propagatorBuilder.addForceModel(sun_3dbodyattraction)

'''
#Atmospheric Drag
msafe = MarshallSolarActivityFutureEstimation(    '(?:Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)\p{Digit}\p{Digit}\p{Digit}\p{Digit}F10\.(?:txt|TXT)',
MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)
DM.feed(msafe.getSupportedNames(), msafe) # Feeding the F10.7 bulletins to 
Orekit's data manager
atmosphere = NRLMSISE00(msafe, sun, wgs84Ellipsoid)
isotropicDrag = IsotropicDrag(sat_list[sc_name]['cross_section'], sat_list[sc_name] 
['cd'])
dragForce = DragForce(atmosphere, isotropicDrag)
propagatorBuilder.addForceModel(dragForce)

# Solar radiation pressure
isotropicRadiationSingleCoeff = IsotropicRadiationSingleCoefficient(sat_list[sc_name]['cross_section'], sat_list[sc_name]['cr']);
solarRadiationPressure = SolarRadiationPressure(sun, wgs84Ellipsoid.getEquatorialRadius(),
                                            isotropicRadiationSingleCoeff)
propagatorBuilder.addForceModel(solarRadiationPressure)
'''

Next I set up the Batch Least Squares Estimator. Again this is pretty much copy pasted from the example from my original post.

########## BATCH LEAST SQUARES ESTIMATOR SETUP

#Defining Optimizer
matrixDecomposer = QRDecomposer(1e-11)
optimizer = GaussNewtonOptimizer(matrixDecomposer, False)

#Defining Estimator
estimator = BatchLSEstimator(optimizer, propagatorBuilder)
estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
estimator.setMaxIterations(estimator_max_iterations)
estimator.setMaxEvaluations(estimator_max_evaluations)

Next is adding the measurements to the estimator object. This is the same from the example post, but instead of Range measurements, there is AngularAzEl measurements. No errors here after the initial fix of defining the ObservableSatellite(0).

########## ADDING MEASUREMENTS

#Process Dataframe into lists for initial orbit determination
az, el, ra, dec, datetime, slantRange = orbitalPass.getFirstAndLast()

#Add Measurements to Estimator
for j in range(1,len(az)):
    orekitAzEl = AngularAzEl(GroundStation(stationFrame),
					datetime_to_absolutedate(datetime[j]),
					JArray('double')([radians(az[j]),radians(el[j])]),
                    JArray('double')([radians(azError),radians(elError)]),
                    JArray('double')([azBaseWeight,elBaseWeight]),
                    ObservableSatellite(0))
    estimator.addMeasurement(orekitAzEl)

And then I get here where the estimator actually runs.

########## STATISTICAL ORBIT DETERMINATION

estimatedPropagatorArray = estimator.estimate()
9

Hi @JPJohn_2015,

I understand perfectly. However without a running example of your bug it will be harder to help you.
Maybe you could generate some measurements on a slightly different orbit, with your two passes separated from 10hours, and send us an example that we could actually run.

Or maybe you did set a BatchLSOrbserver to print the outputs of each iteration of the batch least-square algorithm ? If yes, the printed outputs could help understand what is happening.

Just two things though.

  1. I would set the azBaseWeight and elBaseWeight back to 1.0. By doing what you do in your code I think you are cancelling out the normalization in the inner matrix and residuals’ computation of the BLS algorithm.
    These weights should ideally be set between 0 and 1.
    We usually always set them to 1.0.
  2. I would add the drag, it is one of the major perturbation in LEO. After 10hours your along-track position could be off by hundreds of kms. And this doesn’t help the optimizer converging.

Regards,
Maxime