Estimation convergence differences with Cartesian and Equinoctial orbits

1

I am using BatchLSEstimator to refine an initial guess at an orbit and I’m seeing different convergence results with CartesianOrbit vs EquinoctialOrbit. I’m sure that I’m doing something wrong, but can’t spot the error?

I have a very short arc (less than 2s) of az/el/range measurements. I make an initial guess at position/velocity through simple geometry, then refine that with BLS. I just use an analytical Keplerian propagator due to the short time-span. If I initialise with a CartesianOrbit, then the estimator converges in 4 iterations and the residuals are close to zero. If I use an EquinoctialOrbit, then the estimator converges in 36 iterations and the residuals in range are larger than I’d expect.

This is what I mean in numbers (residuals in az - radians, el - radians, range - m) for Cartesian

4
[-1.2240578425970483E-5, 7.092818138071388E-5, -0.2394066087435931]
[1.7907838861574987E-5, -5.546324877797604E-5, -0.551989545347169]
[-2.7295146693973038E-5, -1.9620036056877677E-5, 1.1203486216254532]
[-1.7782575352676133E-5, 3.693302303431345E-7, 0.30834507080726326]
[1.5954533016948602E-5, 1.2710312892127984E-5, 0.40259317052550614]
[4.103404743904804E-5, -5.351441522638556E-6, -0.9828750300221145]
[3.476988051254182E-5, -2.0326263433345648E-5, -0.6837378430645913]
[-2.527240480620918E-5, 1.4317318524359468E-5, 0.15829963283613324]
[-3.342359139768547E-5, -6.86127512401491E-5, 0.5306026425678283]
[6.347995177069521E-6, 7.104873892149755E-5, -0.062180155189707875]

and for equinoctial

36
[-9.26115061883337E-6, 7.553046300529109E-5, 3.252938335062936]
[2.0872866197230167E-5, -5.145506439457881E-5, 2.216291958699003]
[-2.4344761155425942E-5, -1.620611649721937E-5, 3.1640153906773776]
[-1.4847073456447646E-5, 3.1888177171035537E-6, 1.6268466557376087]
[1.887490891183674E-5, 1.493520138723392E-5, 0.9953799655195326]
[4.393905445398971E-5, -3.7213186055429226E-6, -1.11635178537108]
[3.765927525334334E-5, -1.929107235226768E-5, -1.5440260577015579]
[-2.2398866253769967E-5, 1.4757411848032032E-5, -1.4293470934499055]
[-3.0566153465194645E-5, -6.876792126614406E-5, -1.7849487927742302]
[9.189087538885587E-6, 7.029814028652881E-5, -3.106181636918336]

I’ve attached the code that I used. The only difference that I can see is that I create either a CartesianOrbit or EquinoctialOrbit from my initial position/velocity guess. I’ve played around with the positionAngle and positionScale options to the KeplerianPropagatorBuilder, but with very little difference in the resulting residuals.

StandaloneTest.java (7.9 KB)

2

Hi @markrutten,

I’ve run your code and reproduced the behavior you observed.
And I think you found a bug!
Could you please open an issue for it?

Note that if you use a NumericalPropagator instead of a KeplerianPropagator in the equinoctial estimation it will converge to the same residuals and orbit as with the Cartesian estimation (less than 3e-15 rad for angle residuals and 4nm on the range residuals).
But it will take 19 iterations and 22 evaluations instead of 3 iterations and 4 evaluations.
I tried with a circular orbit, achieving convergence in 27 iterations and 37 evaluations.
With a keplerian orbit it took 105 iterations and 117 evaluations…

So I think the problem is two-fold:

  1. A (big) bug in KeplerianPropagator, probably linked to the state transition matrix computation, when the orbit type is not Cartesian.
    Which strongly looks like this old issue (#351: PartialDerivativesEquations does not work with non-Cartesian elements) although the way the STM is computed now is different than it was by that time.
    Unfortunately, the link to the old mailing list in the issue is dead so I couldn’t read the original poster problem.
  2. A numerical issue, probably related to conversion from Cartesian to Keplerian-like parameters (or the reversed conversion)

Both points definitely need more investigation.

Thanks,
Maxime

3

Thanks @MaximeJ. I have created an issue on gitlab (issue 1073).

4

@MaximeJ, I’ve had a bit of a closer look and your option 1 for the problem looks exactly right. Both AbstractAnalyticalGradientConverter and AbstractAnalyticalMatricesHarvester implicitly assume a Cartesian orbit.

If I hack (really, really hack) the setReferenceState method in AbstractAnalyticalMatricesHarvester to pre/post multiply the Cartesian STM by the appropriate Jacobians, then the estimator converges with an Equinoctial orbit as expected. The code is below, but please don’t use it … it’s a hack because we don’t have access to the orbit’s PositionAngle parameter in that method (I’ve hard-coded it) and I’m not modifying the parameters Jacobians at all.

        ...
        addToRow(derivativesVx, 3);
        addToRow(derivativesVy, 4);
        addToRow(derivativesVz, 5);

        // Jacobians wrt Cartesian
        double[][] stateJacobianStartArr = new double[STATE_DIMENSION][STATE_DIMENSION];
        double[][] stateJacobianTargetArr = new double[STATE_DIMENSION][STATE_DIMENSION];
        reference.getOrbit().getType().convertType(gState.getOrbit().toOrbit()).getJacobianWrtParameters(PositionAngle.MEAN, stateJacobianStartArr);
        reference.getOrbit().getJacobianWrtCartesian(PositionAngle.MEAN, stateJacobianTargetArr);
        
        // Map STM to orbital parameters
        RealMatrix stateJacobianStart = MatrixUtils.createRealMatrix(stateJacobianStartArr);
        RealMatrix stateJacobianTarget = MatrixUtils.createRealMatrix(stateJacobianTargetArr);
        RealMatrix stm = stateJacobianTarget
                .multiply(MatrixUtils.createRealMatrix(analyticalDerivativesStm))
                .multiply(stateJacobianStart);
        for (int row = 0; row < STATE_DIMENSION; ++row) {
            for (int col = 0; col < STATE_DIMENSION; ++col) {
                analyticalDerivativesStm[row][col] = stm.getEntry(row, col);
            }
        }

        // Partial derivatives of the state with respect to propagation parameters
        int paramsIndex = converter.getFreeStateParameters();
        ...
5

Thanks @markrutten for the investigations and for opening the issue.

I’ve tried with other analytical propagators (Eckstein-Hechler and Brouwer-Lyddane) and the results between Cartesian on one side and Keplerian-like on the other side are also very different.
I’ve updated your issue with this information.