# Re-computed estimated measurements

Hey there,

what does `getLastEstimations` of batch least squares return exactly?
Here is what I have done, after an OD:
1/ I iterate over its values (as it’s in Python I cast them in `EstimatedMeasurement`)
2/ for each, I call `getObservedMeasurement` that I cast into the appropriate class e.g. `PV`
3/ on this, I call `getEstimatedValue`
My goal being to retrieve the theoretical measurements according to the fitted orbit.

The problem is that if I try to recompute these theoretical measurements in an independent process (creating a separate propagator with same force model and initial conditions from the OD and calling the observation model), I get non-negligible differences.
I am aware that there are many potential sources of discrepancies in this, so before digging into each of them, I would like to be clear on what steps 1 to 3 are supposed to give me. My fear is that it returns the penultimate iteration before convergence.

Cheers,
Romain.

Hi Romain,

I quickly look at the code (using my mobile phone… not easy ) and I think you are right. It looks like the method returns the estimated measurements based on the state optimized at the penultimate iteration.
I see that the evaluations (i.e., the object returned by the getLastEstimation() method) is updated before the differential correction is applied to the state vector. In other words, the method gives the observed and estimated measurements used to compute the residual vector of the differential correction for the last estimation.
Taking the parallel with the Kalman Filter, they are, for the last estimation, predicated measurements, not corrected.

Best regards,
Bryan

Hi Bryan,

If that is confirmed, what would be the easiest way to retrieve the estimated measurements according to the final, converged (if applicable) iteration of the OD process (at least in the batch least square case )? I’m a bit afraid that trying to do things manually with the estimator might be made quite tedious by the private status of some attributes/methods…

Cheers,
Romain S.

Hi Romain, Bryan,

I’m not sure I understand your debate.

If this is indeed how it works, it looks like a bug to me.

This is what Romain wants I think, the last estimated measurements when the convergence checker says it is ok to stop the iteration process. No ?
If you look up for example in the GaussNewtonOptimizer `optimize` method it does:

• Calls LeastSquareProblem.evaluate function (l. 163): `current = lsp.evaluate(currentPoint);`
In Orekit, this calls AbstractBatchLSModel.value method which runs the propagation, compute the residuals and Jacobian and updates the last estimations (this is done at the end of the method: `observer.modelCalled(orbits, evaluations);`)
• Then checks convergence at line 169 and returns `current` point if the estimator has converged.

So, in my opinion, you should get the last estimations that lead to convergence. Is this what you are looking for Romain ?
And how did you recompute the theoretical measurements Romain ?

In this former post I uploaded some dedicated handlers to get the map of (observed, estimated) measurements from a propagator and a list of measurements in input. Maybe you could use something similar in Python ?

Maxime

If this is indeed how it works, it looks like a bug to me.

I’d say it’s not intuitive for sure, and should at least be made cleared in the documentation (if this is confirmed of course)

This is what Romain wants I think, the last estimated measurements when the convergence checker says it is ok to stop the iteration process. No ?

So, in my opinion, you should get the last estimations that lead to convergence. Is this what you are looking for Romain ?

I want the theoretical measurements (actually, also other results such as location/clock biases, rejection status from outlier detectors, etc.) according to the determined parameters. So if `getLastEstimations` gives theoretical measurements before the last update is performed (no matter how small it is), that doesn’t quite work for me (but would explain my findings, assuming there is no mistake somewhere in my code, which cannot be completely ruled out).

And how did you recompute the theoretical measurements Romain ?

With dedicated approaches depending on the type of data, not directly related to the estimator itself (that’s why I said independently). For example with what I call Cartesian pseudo-measurements (`PV` class), it is pretty simple: I create a propagator with the fitted initial conditions and same settings than the builder (integrator, force model, etc.) and I propagate to the dates of interest. I still get differences that look non-negligible to me (a few centimeters at epoch can lead to significant errors later in time).

Thanks for the link @MaximeJ I will have a closer look.

Cheers,
Romain.

That is indeed what you should get and I think this is what you get.
When calling `BatchLSEstimator.estimate`, the estimator returns propagators that are already set up with estimated parameters (orbit, propagation etc.) and the measurements in last estimation should also have last estimated measurements parameters set up.

I extracted the log of the first OD of tutorial `SequentialBatchLeastSquare` (see below) to have an example.
In my opinion, after `estimate` you should get a propagator that reproduces what you have on the line for iteration 4. @bryan, do you think this is not the case ?

``````First estimation: Batch Least Squares
iteration evaluations      ΔP(m)        ΔV(m/s)           RMS          nb Range    nb Range-rate  nb Angular     nb PV                        Px                      Py                      Pz                      Vx                      Vy                      Vz
1          1                                 7332.378699696918       0/0            0/0          0/0          0/0         -8470598.000000000       -656365.000000000       8683142.000000000          2370.000000000         -4795.000000000          2057.000000000
2          2           9.138601  0.547831386   1.201133120341       0/0            0/0          0/0          0/0         -8470599.384532047       -656363.895219139       8683150.965297380          2369.667811081         -4795.425227808          2057.094610042
3          3           3.659733  0.001021017   0.196154591045       0/0            0/0          0/0          0/0         -8470598.574341394       -656367.094420379       8683152.547177860          2369.668344354         -4795.424849766          2057.093825708
4          4           0.000001  0.000000000   0.196154583718       0/0            0/0          0/0          0/0         -8470598.574340913       -656367.094420622       8683152.547177862          2369.668344354         -4795.424849766          2057.093825708
Estimation time (s): 1.647
``````

Could you maybe provide an example we could work on ?

Cheers,
Maxime

I think the residuals computed at iteration 4 are computed using the state estimated at iteration 3.

Maybe I’m wrong but I think the getLastEstimations() method returns the observed/estimated measurements for iteration 3 because they are the ones used to compute the residuals vector of iteration 4.
However, if I good understand Romain would like the ones of iteration 4. But they are not computed because the optimizer converged and there is no extra iteration. If I understand well, Romain thought that the estimated measurements of iteration 4 are provided by the method getLastEstimaton() (which can be normal according to the method’s name)

The estimated propagator is indeed representative of iteration 4 in terms of estimated parameters.
But the estimated measurements of getLastEstimation() are, I think, reprensentative of the propagator at iteration 3. To get the estimated measurements based on the state estimated at iteration 4, one have to perform an orbit propagation using the estimated propagator.

This is what I think based on my understanding of the code, but maybe I’m wrong.

You are right.
Most of the time spent in optimization corresponds to the evaluation of residuals, because in our case it needs to perform a full propagation. In other non-space flight dynamics uses, this can be different of course, but the evaluation is often considered a costly operation.

Taking this into account, a classical design choice is to have the costly evaluation step at the start of the iteration, followed by problem inversion (i.e. QR decomposition and solving) and parameters updates.

So if one wants to re-evaluate the residuals with latest parameter, one has to run the evaluation one more time (and in this case, it would be tempting to use that for updating the parameters, and we would be back to square one with just one additional iteration performed).

Thanks everyone for the discussion.

I believe I have pinned down at least three sources of discrepancy. As discussed previously, there is the fact that getLastEstimations returns the theoretical values (with or without observation modifiers though?) according to the parameters before the final update.
@MaximeJ can I force the least squares estimator to compute the measurements again and get them instead of doing the independent procedure you proposed in the other thread?

But also, by writing the same procedure directly in Java, in a single script so very linearly (as opposed to a more object-oriented architecture) I realized that my value of Earth gravitational constant was not kept tracked of consistently in Python (optional arguments are to be blamed in part), leading to small differences in some computations.

Last but not least, how the propagation is performed matters, even for a given integrator object. I don’t know exactly what’s done within the estimation, but outside, using or not the so-called former ephemeris mode (interpolation from integration steps with ephemeris generator) for example can give slightly different results. I guess this is just part of the propagation noise, as would arise also with different tolerances or even different integration schemes.

Cheers,
Romain

Hi Romain,

Yes when you call `ObservedMeasurement.estimate` it takes the modifiers into account.

I don’t think so, not easily at least. This is done in `AbstractBatchLSModel.value` method and this model is built on the fly during the call to `BatchLSModel.estimate()` method.
So unless an attribute and a getter are added in the Java class the model is not reachable by the user.

Cheers,
Maxime

Actually I thought of a trick (haven’t tried it yet): creating a new estimator, sequential least squares so that it can be initialized with the previous result including the covariance matrix. And set up the convergence threshold a bit lower so that one iteration is enough to stop. Then getLastEstimations with that one would give me the theoretical measurement I want.

Cheers,
Romain.

Hi everyone,

I am still experiencing discrepancies between the estimator object itself and what I recompute with a separate propagator. I can’t identify the source as one of things I mentioned previously, but it seems linked to drag.
I’ve created a test case in Java that hopefully you guys can run.

Imports:

``````import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.ode.nonstiff.DormandPrince54Integrator;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.estimation.leastsquares.BatchLSEstimator;
import org.orekit.estimation.measurements.EstimatedMeasurement;
import org.orekit.estimation.measurements.ObservableSatellite;
import org.orekit.estimation.measurements.ObservedMeasurement;
import org.orekit.estimation.measurements.PV;
import org.orekit.forces.drag.DragForce;
import org.orekit.forces.drag.IsotropicDrag;
import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.ITRFVersion;
import org.orekit.models.earth.Geoid;
import org.orekit.models.earth.ReferenceEllipsoid;
import org.orekit.models.earth.atmosphere.NRLMSISE00;
import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation;
import org.orekit.orbits.CartesianOrbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.conversion.DormandPrince54IntegratorBuilder;
import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;

import java.io.File;
import java.util.Locale;
import java.util.Map;

``````

Main:

``````public static void main(final String[] args) {
try {

// configure Orekit
final File home       = new File(System.getProperty("user.home"));
final File orekitData = new File(home, "Dev//orekit-data");
if (!orekitData.exists()) {
System.err.format(Locale.US, "Failed to find %s folder%n",
orekitData.getAbsolutePath());
System.err.format(Locale.US, "You need to download %s from %s, unzip it in %s and rename it 'orekit-data' for this tutorial to work%n",
"orekit-data-master.zip", "https://gitlab.orekit.org/orekit/orekit-data/-/archive/master/orekit-data-master.zip",
home.getAbsolutePath());
System.exit(1);
}
final DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();

// Initial date in UTC time scale
final TimeScale utc = TimeScalesFactory.getUTC();
final AbsoluteDate initialDate = new AbsoluteDate(2021, 10, 2, 0, 0, 00.000, utc);

// Frames
final Frame inertialFrame = FramesFactory.getGCRF();
final Frame itrf = FramesFactory.getITRF(ITRFVersion.ITRF_2014, IERSConventions.IERS_2010, false);

// Initial conditions
final double r0 = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 500e3;
final double v0 = FastMath.sqrt(Constants.EGM96_EARTH_MU / r0);
final double mass = 10.;
final double cs = 1.;
final double dragCoef = 2.;
final Vector3D pos = new Vector3D(r0, 0., 0.);
final Vector3D vel = new Vector3D(0., v0, 0.);
final PVCoordinates pv = new PVCoordinates(pos, vel);

// Forces
final int degree = 2;
final int order = 1;
final NormalizedSphericalHarmonicsProvider gravityFields = GravityFieldFactory.getNormalizedProvider(degree, order);
final Geoid geoid = new Geoid(gravityFields, ReferenceEllipsoid.getWgs84(itrf));
final HolmesFeatherstoneAttractionModel holmes = new HolmesFeatherstoneAttractionModel(geoid.getBodyFrame(),
gravityFields);
final MarshallSolarActivityFutureEstimation atmParams = new MarshallSolarActivityFutureEstimation(MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,
MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
final NRLMSISE00 atmo = new NRLMSISE00(atmParams, CelestialBodyFactory.getSun(), geoid);
final double coeffMin = 0.;
final double coeffMax = 1000.;
final IsotropicDrag drag = new IsotropicDrag(cs, dragCoef, coeffMin, coeffMax);
final DragForce dragForce = new DragForce(atmo, drag);

final double mu = holmes.getMu();
final CartesianOrbit orbit = new CartesianOrbit(pv, inertialFrame, initialDate, mu);
final SpacecraftState spacecraftState = new SpacecraftState(orbit, mass);

// Integrators
final OrbitType orbitType = OrbitType.CARTESIAN;
final PositionAngle positionAngle = PositionAngle.TRUE;
final double minStep = 0.1;
final double maxStep = 60.;
final double positionTolerance = 0.1;
final double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, orbit, orbitType);
final DormandPrince54Integrator dp = new DormandPrince54Integrator(minStep, maxStep,
tolerances[0], tolerances[1]);

// Create Cartesian pseudo-measurements by propagating without drag
final double dt = 10000.;
final AbsoluteDate date2 = initialDate.shiftedBy(dt);
final AbsoluteDate date3 = date2.shiftedBy(dt);
final NumericalPropagator propReal = new NumericalPropagator(dp);
propReal.setInitialState(spacecraftState);
propReal.setOrbitType(orbitType);
propReal.setResetAtEnd(false);
final PVCoordinates pv2 = propReal.propagate(date2).getPVCoordinates(inertialFrame);
final PVCoordinates pv3 = propReal.propagate(date3).getPVCoordinates(inertialFrame);

// Guess for estimator (real orbit with error introduced)
final PVCoordinates modifiedPv = new PVCoordinates(pos.add(new Vector3D(1e3, -2e3, 4e2)),
final CartesianOrbit modifiedOrbit = new CartesianOrbit(modifiedPv, inertialFrame, initialDate, mu);

// Propagator builder for estimator
final DormandPrince54IntegratorBuilder dpBuilder = new DormandPrince54IntegratorBuilder(minStep, maxStep,
positionTolerance);
final double positionScale = 1.;
final NumericalPropagatorBuilder propBuilder = new NumericalPropagatorBuilder(modifiedOrbit, dpBuilder,
positionAngle, positionScale);

// Create estimator
final double weight = 1.;
final double sigmaPos = 1e-3;
final double sigmaVel = 1e-6;
final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
final BatchLSEstimator estimator = new BatchLSEstimator(optimizer, propBuilder);
estimator.setMaxEvaluations(60);
estimator.setMaxIterations(40);
estimator.setParametersConvergenceThreshold(1e-4);
final ObservableSatellite sat = new ObservableSatellite(0);
sigmaPos, sigmaVel, weight, sat));
sigmaPos, sigmaVel, weight, sat));
sigmaPos, sigmaVel, weight, sat));
final String name = "drag coefficient";
estimator.getPropagatorParametersDrivers(false).findByName(name).setSelected(true);
estimator.getPropagatorParametersDrivers(false).findByName(name).setReferenceValue(dragCoef);

// perform OD and post process it
final NumericalPropagator estimatedProp = (NumericalPropagator) estimator.estimate()[0];
estimatedProp.setResetAtEnd(false);
final SpacecraftState estimatedState = estimatedProp.getInitialState();
final PVCoordinates estimatedPV = estimatedState.getPVCoordinates(inertialFrame);
final PVCoordinates pv2Fitted = estimatedProp.propagate(date2).getPVCoordinates(inertialFrame);
final PVCoordinates pv3Fitted = estimatedProp.propagate(date3).getPVCoordinates(inertialFrame);

// Create independent propagator and recompute estimated measurements
final double[][] tolerancesBis = NumericalPropagator.tolerances(positionTolerance,
estimatedState.getOrbit(), estimatedState.getOrbit().getType());
final DormandPrince54Integrator dpBis = new DormandPrince54Integrator(minStep, maxStep,
tolerancesBis[0], tolerancesBis[1]);
final NumericalPropagator propBis = new NumericalPropagator(dpBis);
final CartesianOrbit orbitBis = new CartesianOrbit(new PVCoordinates(estimatedPV.getPosition(),
estimatedPV.getVelocity()), inertialFrame, initialDate, estimatedState.getMu());
final SpacecraftState scBis = new SpacecraftState(orbitBis, spacecraftState.getMass());
propBis.setInitialState(scBis);
propBis.setOrbitType(orbitType);
final Geoid geoid2 = new Geoid(gravityFields, ReferenceEllipsoid.getWgs84(itrf));
final double fittedCoef = estimator.getPropagatorParametersDrivers(true).findByName(name).getValue();
final IsotropicDrag newDrag = new IsotropicDrag(cs, fittedCoef, coeffMin, coeffMax);
final NRLMSISE00 atmo2 = new NRLMSISE00(atmParams, CelestialBodyFactory.getSun(), geoid2);
propBis.setResetAtEnd(false);
final PVCoordinates pv2Repropagated = propBis.propagate(date2).getPVCoordinates(inertialFrame);
final PVCoordinates pv3Repropagated = propBis.propagate(date3).getPVCoordinates(inertialFrame);

// assess differences between propagators
final double[] diffPos2 = (pv2Fitted.getPosition().subtract(pv2Repropagated.getPosition())).toArray();
final double[] diffVel2 = (pv2Fitted.getVelocity().subtract(pv2Repropagated.getVelocity())).toArray();
final double[] diffPos3 = (pv3Fitted.getPosition().subtract(pv3Repropagated.getPosition())).toArray();
final double[] diffVel3 = (pv3Fitted.getVelocity().subtract(pv3Repropagated.getVelocity())).toArray();

// look at estimated measurements (for debug)
final Map<ObservedMeasurement<?>, EstimatedMeasurement<?>> estimatedMeasurement = estimator.getLastEstimations();
for (EstimatedMeasurement meas: estimatedMeasurement.values()) {
final double[] estimatedValue = meas.getEstimatedValue();
final double[] observedValue = meas.getObservedValue();
}

} catch (OrekitException oe) {
System.err.println(oe.getLocalizedMessage());
}
}
``````

Specifically, the issue is that `diffPositionX` is non-negligible.

N.B.: if some objects are re-instantiated rather than reused, it is because this is what my original Python code does and I’ve tried to change things as little as possible when transcribing it into Java.

Cheers,
Romain.

Hi @Serrof,

I found the discrepancy: you should call `setMass(mass)` on the propagator builder object otherwise it takes 1000 kg as default mass.

1 Like

Thanks a lot @MaximeJ, dunno how long I would have taken to figure that one out!

Now I can ask another question. I’ve tried to implement my trick to get the “real” estimated measurements, creating a second estimator for a “useless” iteration. But it still iterates twice and I don’t understand why.

Here is how I have set it up:

``````hipparchus_optim = SequentialGaussNewtonOptimizer().withEvaluation(self._java_estimator.getOptimum())
second_estimator = SequentialBatchLSEstimator(hipparchus_optim, [new_prop_builder.propagatorBuilder])
second_estimator.setMaxEvaluations(300)
second_estimator.setMaxIterations(100)
second_estimator.setParametersConvergenceThreshold(self.threshold * 10.)
for meas in self._measurements:
for driver in self._java_estimator.getOrbitalParametersDrivers(True).getDrivers():
name = driver.getName()
second_driver = second_estimator.getOrbitalParametersDrivers(False).findByName(name)
second_driver.setSelected(True)
second_driver.setReferenceValue(driver.getValue())
for driver in self._java_estimator.getPropagatorParametersDrivers(True).getDrivers():
name = driver.getName()
second_driver = second_estimator.getPropagatorParametersDrivers(False).findByName(name)
second_driver.setSelected(True)
second_driver.setReferenceValue(driver.getValue())
for driver in self._java_estimator.getMeasurementsParametersDrivers(True).getDrivers():
name = driver.getName()
second_driver = second_estimator.getMeasurementsParametersDrivers(False).findByName(name)
second_driver.setSelected(True)
second_driver.setReferenceValue(driver.getValue())
second_estimator.estimate()
``````

It’s in Python, sorry for the Java purists.
`new_prop_builder.propagatorBuilder` is the Orekit propagator builder I create with the fitted initial conditions.
`self._measurements` is a list of Orekit object also passed to the original estimator (`self._java_estimator`).

Cheers,
Romain.

You’re welcome Romain, actually you should thank the Java debugger

Here’s a piece of code from Hipparchus `GaussNewtonOptimizer.optimize` method:

``````        // iterate until convergence is reached
Evaluation current = null;
while (true) {
iterationCounter.increment();

// evaluate the objective function and its jacobian
Evaluation previous = current;
// Value of the objective function at "currentPoint".
evaluationCounter.increment();
current = lsp.evaluate(currentPoint);
final RealVector currentResiduals = current.getResiduals();
final RealMatrix weightedJacobian = current.getJacobian();
currentPoint = current.getPoint();

// Check convergence.
if (previous != null &&
checker.converged(iterationCounter.getCount(), previous, current)) {
return Optimum.of(current,
evaluationCounter.getCount(),
iterationCounter.getCount());
}
``````

In Gauss-Newton optimize method iteration 1 is the first run, with `previous = null`, the convergence is not checked.
Then on iteration 2, previous = iteration 1, current = iteration 2 and then the convergence checker stops the process.

Note:

1. You can use the same propagator builder (you don’t have to create a new one), it should have its parameter drivers already updated with the values that lead the first estimator to convergence;
2. I’m not sure you need a sequential batch neither, it seems to me that you’re using measurements’ info twice (one from the covariance, second from the measurements).
But that should not change the results much since you’re already converged.

Hope this helps,
Maxime

Hi @MaximeJ,

In Gauss-Newton optimize method iteration 1 is the first run, with `previous = null`, the convergence is not checked.

This seems a bit counterintuitive to me, especially if I have initialized the optimizer with an evaluation. Shouldn’t it be able to use it for the convergence check as well?

You can use the same propagator builder (you don’t have to create a new one), it should have its parameter drivers already updated with the values that lead the first estimator to convergence;

Thanks for the tip

I’m not sure you need a sequential batch neither

See my first point

Cheers,
Romain.

Hi @Serrof,

True, maybe line 306 of Hipparchus SequentialGaussNewton:

``````// iterate until convergence is reached
Evaluation current = null;
``````

Should be replaced by:
`Evaluation current = oldEvaluation`

To be discussed with @bcazabonne I’d say.

I still think it would be simpler for you to get the last estimations with your own Python version of a ResidualsHandler.java (3.8 KB)

I forgot to ask: if the second estimator has made two iterations, does it mean it has updated at least once the estimated parameters? In that case I am back to square one with theoretical measurements that are slightly off w.r.t. the fitted orbit.

Cheers,
Romain.

I don’t remember why we initialized the value to `null` but yes using `oldEvaluation` looks good. Interesting to see if tests still pass in both Orekit and Hipparchus after the update.

Good question, I don’t think so, at least the number of iterations will be different.

I thought about this yesterday and I’m not sure we should do it, I’ll have a look at the equations again to see if this would be right with respect to a “canonical” sequential BLS.
At least we could add a method “initalizeWithOldEvaluation” that would do the trick.

Yes I think so Romain…