Re-computed estimated measurements

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

Thanks Maxime for your answers :grinning:

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.ICGEMFormatReader;
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();
        manager.addProvider(new DirectoryCrawler(orekitData));

        // 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
        GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader(GravityFieldFactory.ICGEM_FILENAME, false));
        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.addForceModel(holmes);
        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)),
                vel.add(new Vector3D(-3e-4, 5e-5, 1e-4)));
        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);
        propBuilder.addForceModel(holmes);
        propBuilder.addForceModel(dragForce);

        // 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);
        estimator.addMeasurement(new PV(initialDate, pv.getPosition(), pv.getVelocity(),
                sigmaPos, sigmaVel, weight, sat));
        estimator.addMeasurement(new PV(date2, pv2.getPosition(), pv2.getVelocity(),
                sigmaPos, sigmaVel, weight, sat));
        estimator.addMeasurement(new PV(date3, pv3.getPosition(), pv3.getVelocity(),
                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));
        propBis.addForceModel(new HolmesFeatherstoneAttractionModel(geoid2.getBodyFrame(), gravityFields));
        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.addForceModel(new DragForce(atmo2, newDrag));
        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:
    second_estimator.addMeasurement(meas)
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 :wink:

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 :grinning:

I’m not sure you need a sequential batch neither

See my first point :sweat_smile:

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) :wink:

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…

Hi @Serrof,

I’m also stuck with this issue. Have you solved it? Especially for the satellite ajisai with TimeSpanDragForce.

Here is my pseudo-code (python):
(1) measurement_list = read_crd(…);
(2) estimatedPropagators = estimator.estimate();
(3) estProp = estimatedPropagators[0];
(4) estimatedMeasurements = estimator.getLastEstimations().values()
(5) evaluation_list = evaluate_measurement(measurement_list, estProp);

The method evaluate_measurement is quite similar to ResidualsHandler.

Then I see, the residuals of re-computed by (5) are quite different with the residuals computed from estimatedMeasurements by (4).

Additionally, it is better to use Ephemeris mode.

    generator = estProp.getEphemerisGenerator()
    estProp.propagate(absdateStart, datePredEnd)
    ephProp = generator.getGeneratedEphemeris()

The two attachment are pngs of residials. ‘resid_’ is the raw one, ‘qc_est_’ is the result of re-compute.

Hi there,

No I never formally manage to retrieve the exact residuals of the last step. However by recomputing some “by hand” I noticed that as expected with sufficiently tight convergence tolerance for the estimator, the state is basically the same for the last two iterations, and so are the residuals consequently.
Do you specify the mass in your propagation builder? I know it can be a source of discrepancy between estimation and independent repropagation as was my case above

Best,
Romain.

Hi @Serrof,

I used the estimated propagator returned by the estimator.estimate(), not re-build a new one. The mass is set in the propagator builder, which is used for the estimation. So I think the mass should be set properly. I’ll check it.
Also I’ll try to re-build a new one with the estimated state and other estimated parameters.
Thanks.

Best,
lirw

Hi,

I was thinking again on this and I can think of two sources of discrepancy:

  1. Precompensation: the dates when the states are interpolated in the batch LS are pre-compensated with an offset coming from previous iteration.
    Meaning:
  • At iteration 0 the date of the state is initially the date of the measurement, the spacecraft state is interpolated at this date and then corrected for light-time delay within a loop using shiftedBys
  • The date found for the corrected state is then registered and used in the next iteration
  • At iteration 1 the interpolator uses the corrected date, so the state position evaluation is closer to the emitting position and the light-time correction is smaller. Thus the evaluation is “better” because it relies more on the interpolator and less on shiftedBys
    → When we try to recompute the residuals we’ve lost this pre-compensation, so we rely on shiftedBys for light time correction and it can introduce differences
  1. Interpolation itself: when recomputing the residuals if one do not use the exact same propagator/integrator/handler then the interpolation of the states may be different.
    For example, if one propagates measurement by measurement to get the estimated measurement value, then the spacecraft state value (and thus the estimated measurement) is different than the one interpolated inside the MeasurementStepHandler during the batch LS OD.

Best,
Maxime

Hi, @MaximeJ

Thank you so much for reply. I’ll check the two sources carefully.

I checked the two I mentioned before.
(1) The mass of the estimated propagator returned by the estimator.estimate() is already set properly, which is the value set for propagator builder used by estimator.
(2) I tried re-build a new propagator with the estimated orbit and other estimated parameters, the residual is same as the estimated propagator.

ps: Both propagator are used in ephemeris mode firstly, rather than propagated measurement by measurement. Formerly, the re-computed residuals are wired by propagating measurement by measurement. If ephemeris mode is used, the residual is much acceptable.

    estimatedPropagators = estimator.estimate()
    estProp = estimatedPropagators[0]

    # New a propagatorBuilder similar to the one used in the estimator,
    # except the orbit is the estimated orbit.
    estOrbit = estProp.getInitialState().getOrbit()
    propagatorBuilderConfig = NumericalPropagatorBuilderConfig(estOrbit)
    propagatorBuilder = propagatorBuilderConfig.configure(cfg)

    # update the estimated parameters. Here are Cds for ajisai.
    propagatorParameterDriversList = propagatorBuilder.getPropagationParametersDrivers()
    for name in ['Cd0', 'Cd1', 'Cd2']:
        value = estPropagatorParameterDriversList.findByName(name).getValue()
        propagatorParameterDriversList.findByName(name).setValue(value)

    normalizedParameters = propagatorBuilder.getSelectedNormalizedParameters()
    print(normalizedParameters)
    prop = propagatorBuilder.buildPropagator(normalizedParameters)

    # Enter ephemeris mode, 6 days after the od epoch.
    # Data between absdateOd and absdateEnd (3 days later) are used to estimated the orbit.
    # And it is predicted for extractor 3 days.
    absdateStart = absdateOd
    absdateEnd = absdateOd.shiftedBy(3 * Constants.JULIAN_DAY)
    datePredEnd = absdateEnd.shiftedBy(3 * Constants.JULIAN_DAY)
    generator = prop.getEphemerisGenerator()
    prop.propagate(absdateStart, datePredEnd)
    ephProp = generator.getGeneratedEphemeris()