Extracting partials after a Batch Least squares estimation

Hello Orekit experts,

I’m using BatchLSEstimator for my OD problem and I have an use-case that demands for getting:
a. State transition matrix, Φ(t, t0), of the state
b. Measurement Jacobian, H-matrix, with respect to the last estimated state

It would appreciated if someone could provide tips on how to efficiently extract / compute them after calling the estimator.estimate(); step. Or is it possible to leverage on OrbitDeterminationObserver to achieve this?

FYI I’m using Orekit v12.1 currently.

There’s a similar topic in this thread: Partial derivatives of a range measurement with respect to the initial PV coordinates , but I believe its outdated as there’s significant changes been done in the classes/interfaces of the code that were suggested, which aligns with the statement quoted in Partial derivatives of a range measurement with respect to the initial PV coordinates - #8 by MaximeJ

I also wonder if there’s any easy API already available. I noticed there exists StateTransitionMatrixGenerator but could not find documentation or relevant tests that supports for how/when to use it.

Any help would be really useful. Thank you in advance.

Link to StateTransitionMatrixGenerator : Orekit/src/main/java/org/orekit/propagation/numerical/StateTransitionMatrixGenerator.java at release-12.1 · CS-SI/Orekit · GitHub

Hi @astroPanda,

Welcome to the forum!

  1. In an OrbitDeterminationObserver , in method evaluationPerformed(...) you can use lspEvaluation.getJacobian().

  2. At the end of the estimation process, use estimator.getOptimum().getJacobian(), where estimator is your BatchLSEstimator.

Beware that both methods return a normalized matrix.
To “un-normalize” it, follow the same logic as in BatchLSEstimator.getPhysicalCovariances(…): divide each line by the σ of the measurement, and each column by the scale of the estimated parameter.

  • Once the estimation is finished.
  • Get the estimated propagator:
    final Propagator estimatedPropagator = estimator.estimate()[0];
  • Configure it for state transition matrix computation
    MatricesHarvester harvester = estimatedPropagator.setupMatricesComputation("STM", null, null);
  • Propagate from t0 to t:
    SpacecraftState state = estimatedPropagator.propagate(t0, t)
  • Get the orbital state transition matrix:
    RealMatrix dYdY0 = harvester.getStateTransitionMatrix(state);
  • Or the parameters’ state transition matrix
    RealMatrix dYdPp = harvester.getParametersJacobian(state);

There’s a place in Orekit where this is a bit more commented and explained, see KalmanModel.getErrorStateTransitionMatrix().

Hope I’ve answered your questions.
Cheers,
Maxime

1 Like

Hello @MaximeJ ,

Thank you so much for the suggestions. I will give them a try and revert back if there is any questions.

In the KalmanModel.getErrorStateTransitionMatrix(), there’s normalization applied, I suppose I should be skipping it right?

Yes, you can skip it.
The link to the Kalman example was just to show the overall structure of the state transition matrix (see the comment for example).

Understood! Thanks.

On the measurement jacobian side, I’m a bit confused. The lspEvaluation.getJacobian() appears to be a matrix of n rows where n representing the number of flattened (angular = 2 rows) measurements where as the physicalCovariance is just one matrix at t0 where nxm .

Would you mind sharing a sample snippet on the un-normalize part please?

I suppose this is the right way to do it??

private RealMatrix computeMeasurementJacobianWrtCurrentOrbit(final EstimatedMeasurement<?> estimated,
                                                                 final Orbit currentOrbit) {
        final double[] sigma = estimated.getObservedMeasurement().getTheoreticalStandardDeviation();
        // Partial derivatives of the current Cartesian coordinates with respect to current orbital state
        final double[][] aCY = new double[6][6];
        currentOrbit.getJacobianWrtParameters(builder.getPositionAngleType(), aCY);   //dC/dY
        final RealMatrix dCdY = new Array2DRowRealMatrix(aCY, false);

        // Jacobian of the measurement with respect to current Cartesian coordinates
        final RealMatrix dMdC = new Array2DRowRealMatrix(estimated.getStateDerivatives(0), false);
        // Jacobian of the measurement with respect to current orbital state
        RealMatrix dMdY = dMdC.multiply(dCdY);
        final ParameterDriversList driversList = getOrbitalParametersDrivers(true);
        // Fill the normalized measurement matrix's columns related to estimated orbital parameters
        for (int i = 0; i < dMdY.getRowDimension(); ++i) {
            for (int j = 0; j < dMdY.getColumnDimension(); ++j) {
                final ParameterDriver driver = driversList.getDrivers().get(j);
                final double partial = dMdY.getEntry(i, j);
                final double unNormalized = partial / sigma[i] / driver.getScale();
                dMdY.setEntry(i, j, unNormalized);
            }
        }
        return dMdY;
    }

  private RealMatrix computeMeasurementJacobianWrtPropagationParameters(final SpacecraftState currentState,
                                                                        final MatricesHarvester harvester,
                                                                        final RealMatrix dMdY,
                                                                        final double[] sigma,
                                                                        final double[] weight) {
      // Jacobian of the measurement with respect to propagation parameters
      final RealMatrix dYdPp = harvester.getParametersJacobian(currentState);
      RealMatrix dMdPp = dMdY.multiply(dYdPp);
      final ParameterDriversList driversList = getOrbitalParametersDrivers(true);
      for (int i = 0; i < dMdPp.getRowDimension(); ++i) {
          for (int j = 0; j < dMdPp.getColumnDimension(); ++j) {
              final ParameterDriver driver = driversList.getDrivers().get(j);
              final double partial = dMdPp.getEntry(i, j);
              final double unNormalized = weight[i] * partial / sigma[i] / driver.getScale();
              dMdPp.setEntry(i, j, unNormalized);
          }
      }
      return dMdPp;
  }

Hi @astroPanda,

Sorry there was a misunderstanding…

What you get from currentOrbit.getJacobianWrtParameters and harvester.getParametersJacobian is not normalized. So you don’t need to un-normalize it.

What is normalized is the output of lspEvaluation.getJacobian() in OrbitDeterminationObserver, or the output of BatchLSEstimator.getOptimum().getJacobian().

But if you don’t use them you don’t need to bother with normalization.

apologies for the confusion.

if I have to rely on lspEvaluation.getJacobian() then the method should look like this I suppose (based on your comment:

)

public RealMatrix getMeasurementJacobian(final LeastSquaresProblem.Evaluation lspEvaluation,
                                             final EstimationsProvider evaluationsProvider,
                                             final ParameterDriversList estimatedOrbitalParameters,
                                             final ParameterDriversList estimatedPropagatorParameters,
                                             final ParameterDriversList estimatedMeasurementsParameters) {
        final RealMatrix jacobians;
        // get the normalized matrix
        jacobians = lspEvaluation.getJacobian().copy();

        // retrieve the measurement sigmas
        final double[] sigma = new double[jacobians.getRowDimension()];
        int index = 0;
        for (int number = 0; number < evaluationsProvider.getNumber(); number++) {
            final EstimatedMeasurement<?> estimated = evaluationsProvider.getEstimatedMeasurement(number);
            for (double s : estimated.getObservedMeasurement().getTheoreticalStandardDeviation()) {
                sigma[index++] = s;
            }
        }

        final double[] scale = getScale(
                estimatedOrbitalParameters,
                estimatedPropagatorParameters,
                estimatedMeasurementsParameters,
                jacobians.getColumnDimension());

        // unnormalize the matrix, to retrieve physical jac
        for (int i = 0; i < jacobians.getRowDimension(); ++i) {
            for (int j = 0; j < jacobians.getColumnDimension(); ++j) {
                jacobians.setEntry(i, j,  jacobians.getEntry(i, j) / (sigma[i] * scale[j]));
            }
        }
        return jacobians;
    }

But the values appears to be inconsistent with dMdY and dMdPp from my previous comment. What could I be possibly going wrong?

Hi @astroPanda,

I think this line

Should be replaced with:

jacobians.getEntry(i, j) * sigma[i] / scale[j]

Thanks! @MaximeJ

Just to make sure I understand correctly. Could you please provide some light on what’s difference between from lspEvaluation.getJacobian() -> (un-normalized) vs this

// for first measurement
estimated = evaluationsProvider.getEstimatedMeasurement(0)

// Partial derivatives of the current Cartesian coordinates with respect to current orbital state
final double[][] aCY = new double[6][6];
currentState.getOrbit().getJacobianWrtParameters(builder.getPositionAngleType(), aCY);   //dC/dY
final RealMatrix dCdY = new Array2DRowRealMatrix(aCY, false);
// Jacobian of the measurement with respect to current Cartesian coordinates
final RealMatrix dMdC = new Array2DRowRealMatrix(estimated.getStateDerivatives(0), false);
final RealMatrix dMdY = dMdC.multiply(dCdY);
final RealMatrix dYdY0 = harvester.getStateTransitionMatrix(estimated.getStates()[0]);
// Jacobian of the measurement with respect to initial
final RealMatrix dMdY0 = dMdY.multiply(dYdY0);

I thought they are supposed to be the same, no? Apologies if my question sounds silly. Still learning OREKIT :slight_smile:

Actually I figured out what’s the issue. Apparently when the model is in reverse propagation mode, the Jacobians are saved in reverse order but the evaluationProvider does not know that. and so my workaround to to force the for loop in the backward order.

Ah, good catch, note that the forward/backward propagation choice is automatic in the OD (see here) to optimize performances.

As for your first question, I think I forgot to link this piece of code that calculates the Jacobian in BatchLSModel. You’ll find here how the normalized matrix is computed.

Your question isn’t silly, it’s not obvious how the Jacobians are computed in Orekit!

Cheers,
Maxime

1 Like

Thanks very much @MaximeJ for sharing this. That was very helpful.