Question on State Transition Matrices & Covariance

I’ve tried what I described above and i’m still experiencing mahalanobis distances which is enormous. Here is the new calculation area:

NumericalPropagatorBuilder numPropBuild = getPropagatorEq(orb1.getOrbit());
double[] normalizedParameters = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
NumericalPropagator numProp = numPropBuild.buildPropagator(normalizedParameters);
PartialDerivativesEquations PDE = new PartialDerivativesEquations(“dYdY0”, numProp);
orb1 = PDE.setInitialJacobians(orb1);
numProp.setInitialState(orb1);

  	SpacecraftState orbf = (numProp.propagate(dt2));
  	
      Array2DRowRealMatrix dYdY0 = new Array2DRowRealMatrix(6, 6);
      PDE.getMapper().getStateJacobian(orbf, dYdY0.getDataRef());
      //PDE.getMapper().getParametersJacobian(orbf, dYdY0.getDataRef());
  	System.out.println("\n\n Final Jacobian: " + dYdY0);
  	System.out.println("Please be equinoctial:"+orbf.getOrbit().getType());
  	
      double[][] deqdcart = new Array2DRowRealMatrix(6, 6).getData();
  	orbf.getOrbit().getJacobianWrtCartesian(PositionAngle.MEAN, deqdcart);
  	RealMatrix deqdcartRM = new Array2DRowRealMatrix(deqdcart);
  	RealMatrix idedc = MatrixUtils.inverse(deqdcartRM); 
      RealMatrix stmf = deqdcartRM.multiply(dYdY0).multiply(idedc);
      RealMatrix stmi = MatrixUtils.inverse(stmf);

With this regime I get some NaN outputs and outrageously large kmn number.

My secondary hunch was that .getPhysicalCovariance is not giving me the correct number. When I use the optimum.getCovariance the output is much closer, but I’m not just going to change it because it fits expected behavior rather than fits mathematics.

Thanks again for all your help. If you know any additional info let me know!

Sincerely,
Paul

The initial and final Jacobians are different because they depend on the orbital elements. You need to compute:

d(eq_0)/d(cart_0) * d(cart_0)/d(cart_1) * d(cart_1)/d(eq_1)

The subscripts indicate different times. As you can see the first and last terms are separate quantities. So you need to use orb1 to compute the first term instead of orbf.

Ah thanks so much I’ll try this out!

Sincerely,
Paul

Sorry one more question is the d(cart_0) in reference to the initial jacobian dYdy0 and would be identity, and d(cart_1) the dYdY0 after the propegation?

Then this computation is equal to that state transition matrix?

Hi Paul,

Perhaps a text book will explain it better than I can in a message board. Do you have Montenbruck & Gill? It’s a good book that explains the concepts in chapter 7. It’s essentially just applying the chain rule a few times. The middle term is the Cartesian STM. The first and last terms are the Jacobians of the coordinate transformations you do to convert to Cartesian and from Cartesian, respectively.

Regards,
Evan

I just found a copy of that book! I think you also did a good job of explaining it, especially after this last message.

I’m still getting vastly different answers as to which covariance I use between the optimum.getCovariance and the estimator.getPhysicalCovariance. Looking through the declaration I see that the physicalCovariance uses the get******ParametersDrivers.getDrivers to get the scales which it uses to “unormalize” the matrix.

Why is the matrix normalized in the first place?
What are those drivers?

Thanks again for all your help thus far it has been extraordinarily useful and helpful!

Your welcome. Normalization is designed to avoid numerical issues when computing the state updates. It is easier for the algorithm if all quantities are approximately the same magnitude, but in physical units position in meters is many orders of magnitude larger than a drag coefficient, for example. So the the drivers translate back and forth between the normalized values that the numeric algorithms like to see and the physical values that astrodynamicists like to see. To see the difference try setting the referenceValue to zero, scale to 1 and min/max value to -/+ infinity before running the orbit determination.

1 Like

Hey Evan,

I am attempting to produce a more accurate state transition matrix to propagate the covariance. I think using the fix discussed above is not accurate enough at the level required for my application. I am attempting to use a field propagator, and ds factory to create a derivative structure representative of the state transition field but I haven’t yet had success.

Could you point me to an example of how to set up a state transition field or how the field propagator works? I haven’t been able to figure it out

Hi Paul,

About the use of the field propagator, did you look at the Orekit tutorial about it ? It provides elements on how to set up the DS Factory depending the derivatives you want and how to set up the field propagator.

Regards,
Bryan

1 Like

I think it makes sense how I could use that tutorial to propagate the standard deviation of each orbital element (expanding to basically a dAdEdIdPAdRAANdv version with 6 params), but I’m not sure how I would do so for specifically for finding the state transition matrix, or how to propagate the covariance instead.

Sorry if this is a dumb question, but what derivative structures would be applicable for finding the state transition matrix or alternatively how to propegate the 6 x 6 covariance matrix? I understand the stm would be ∂x(t)/∂x(t0) but it is just unclear to me how to code that up with the derivative structure.

Sincerely,
Paul

I am not familiar with covariance matrix propagation. However, it exists a simple way in Orekit to obtain the state transition matrix without having to use the DerivativeStructrure instance. You can use the PartialDerivativesEquations class. Here this is an example on how to use this class to obtain the state transition matrix.

public static void main(String[] args) {

    // Configure Orekit
    File home       = new File(System.getProperty("user.home"));
    File orekitData = new File(home, "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);
    }
    DataProvidersManager manager = DataProvidersManager.getInstance();
    manager.addProvider(new DirectoryCrawler(orekitData));

    // Force model 5*5 gravity field force model
    NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
    ForceModel gravityField =
        new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
    
    // INITIAL ORBIT
    // a = 8000000 m
    // e = 0.01
    // i = 0.1 rad
    // ω = 0.7 rad
    // Ω, = 0 rad
    // anomaly = 1.2 rad
    Orbit initialOrbit =
            new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.TRUE,
                               FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
                               provider.getMu());

    // Set up propagator
    double dP = 0.001; // position error (user specification)
    final double minStep = 0.001;
    final double maxStep = 1000;
    double[][] tol = NumericalPropagator.tolerances(dP, initialOrbit, OrbitType.KEPLERIAN);
    NumericalPropagator propagator = new NumericalPropagator(new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]));
    propagator.setOrbitType(OrbitType.CARTESIAN); // Set orbit type
    propagator.setPositionAngleType(PositionAngle.MEAN); // Set position angle type
    propagator.addForceModel(gravityField); // Add the force model used
    
    // Add partial derivative equations as additional equations
    PartialDerivativesEquations partials = new PartialDerivativesEquations("partials", propagator);
    final SpacecraftState initialState =
                    partials.setInitialJacobians(new SpacecraftState(initialOrbit));
    propagator.setInitialState(initialState);
    final JacobiansMapper mapper = partials.getMapper();
    PickUpHandler pickUp = new PickUpHandler(mapper, null);
    propagator.setMasterMode(pickUp);
    propagator.propagate(initialState.getDate().shiftedBy(900.));
    
    // Partial derivatives of the state with respect to initial state
    // In this case ∂x(t)/∂x(t0)
    double[][] dYdY0 = pickUp.getdYdY0();

}

private static class PickUpHandler implements OrekitStepHandler {

    private final JacobiansMapper mapper;
    private final AbsoluteDate pickUpDate;
    private final double[][] dYdY0;

    public PickUpHandler(JacobiansMapper mapper, AbsoluteDate pickUpDate) {
        this.mapper = mapper;
        this.pickUpDate = pickUpDate;
        dYdY0 = new double[6][6];
    }

    public double[][] getdYdY0() {
        return dYdY0;
    }

    public void init(SpacecraftState s0, AbsoluteDate t) {
        // nothing
    }

    public void handleStep(OrekitStepInterpolator interpolator, boolean isLast)
        {
        final SpacecraftState interpolated;
        if (pickUpDate == null) {
            // we want to pick up the Jacobians at the end of last step
            if (isLast) {
                interpolated = interpolator.getCurrentState();
            } else {
                return;
            }
        } else {
            // we want to pick up some intermediate Jacobians
            double dt0 = pickUpDate.durationFrom(interpolator.getPreviousState().getDate());
            double dt1 = pickUpDate.durationFrom(interpolator.getCurrentState().getDate());
            if (dt0 * dt1 > 0) {
                // the current step does not cover the pickup date
                return;
            } else {
                interpolated = interpolator.getInterpolatedState(pickUpDate);
            }
        }
        mapper.getStateJacobian(interpolated, dYdY0);
    }

}

This example refers to the validation class of the PartialDerivativesEquations class available here.

I hope that this will help you,
Bryan

1 Like

Hey thanks for the for the suggestion. Unfortunately like I was discussing above I need the state transition matrix in equinoctial and Keplarian elements. Using the methodology you put out only.produces the Cartesian elements. Further in the discussion above we discuss using partials to change the elements but there are small non linearities, which are a problem for my application.

Supposedly using the partial derivative structure doesn’t have those issues which is why Im trying to figure out how to do so

Hi @Paul,

I wrote you up a little code comparing state transition matrix computation with Orekit partial derivatives class on one side and derivative structures on the other side.
TestSTMWithDS.java (10.9 KB)
I’ll let you try it and see if you can use it in your application.
Hope it will help.

@evan.ward I’ve been thinking for a while about this bug in the PDE providing only Cartesian derivatives. And I don’t think that this is the source of the problem.
If you look into the results of the code I sent to Paul you’ll see that the PDE and JacobianMapper classes do give the state transition matrix with the proper orbit type but that the computation is quite bad. There is a lot of numerical “noise”, especially for the derivatives of the semi major axis.
My guess is that it comes from numerical errors in getJacobianWrtParameters which uses a matrix inversion on an un-normalized matrix (see method Orbit.createInverseJacobian)
But I’m not 100% sure and it would require more testing.
What do you think ?

Regards,
Maxime

Thanks so much for the help! I adapted and tried the program you sent me. The methodology looks so much better than I was attempting. Finding the STM this way works great and has eliminated some of my error sources, but not all! Thanks again for your help

Hi Paul,

I’m facing similar issues and am not sure about the way to compute the transition matrix to convert the cartesian matrix into an equinoctial matrix. I was wondering if it would be possible for you to make this piece of code available? This would be very helpful.

Thank you very much!

B.

Hi @benoist,

Do you want to compute:

  1. A state transition matrix, that contains the partial derivatives of the final state with respect to the initial state following a propagation from initial to final ?
    If so, the code sent in TestSTMWithDS.java sent 2 posts above yours should help you. We can talk more if you need more insight.

  2. The Jacobian of the transform from Cartesian to Equinoctial ?
    If so you can use the method getJacobianWrtCartesian from the Orbit class.

Tell us more, we’re here to help!
Maxime

Hi Maxime,

Thank you very much for your answer!

Well to be completely honest I’m far from being an expert on the topic and I’m a bit lost in all these transformations. I’m basically trying to tune a Kalman filter I’m using for orbit determination. It used to handle a state vector in PV (x,y,z,Vx,Vy,Vz) but I would like to change to equinoctial parameters to hopefully have a better propagation.

The Kalman filter has already been developed and validated for Cartesian parameters but I’m truly struggling on the migration to equinoctial.

I’ll try to give as much details as I can, please see below.

I’m first defining the state vector.

double[] stateVectorKalman = {
nextStepState.getA(), 
nextStepState.getEquinoctialEx(), 
nextStepState.getEquinoctialEy(), 
nextStepState.getHx(), 
nextStepState.getHy(), 
nextStepState.getLv(),
clockBiasEst,
clockDriftEst,
clockLongTermBiasEst,
thrustEstX,
thrustEstY,
thrustEstZ};

Then I want to create the transition matrix F to propagate the covariance matrix, following the equation below:

// Predicted covariance estimate */
kalmanP = (kalmanF.multiply(kalmanP).multiply(kalmanF.transpose())).add(kalmanQ); // P = F*P*Ft + Q

Here is what I used to do to create the matrix F when I used to handle Cartesian parameters. I was hoping introducing equinoctial parameters would not change anything, but there seems to be issues and I’m not sure what’s the best way forward…

// Definition of partials
PartialDerivativesEquations partials = new PartialDerivativesEquations("partials", propagatorReference);
previousStepState = partials.setInitialJacobians(previousStepState);		
propagator.setInitialState(previousStepState);

// State transition (Jacobian) matrix computation */
JacobiansMapper mapper = partials.getMapper();
double [][] jacobian = new double[JacobiansMapper.STATE_DIMENSION][JacobiansMapper.STATE_DIMENSION];
mapper.getStateJacobian(nextStepState, jacobian);

		
kalmanF = new Array2DRowRealMatrix(nbOfState, nbOfState);
kalmanF.setSubMatrix(jacobian, 0, 0);
kalmanF.setEntry(6, 6, 1);
kalmanF.setEntry(6, 7, dt);
kalmanF.setEntry(6, 8, 0.5*FastMath.pow(dt, 2));
kalmanF.setEntry(7, 7, 1);
kalmanF.setEntry(7, 8, dt); 
kalmanF.setEntry(8, 8, 1);

for (int j = 0; j < 5; j++){
   for(int i = 0; i < 6; i++){
	kalmanF.setEntry(i, 9 + j, jacobianParam[i][j]); // da/dTx, dex/dTx, dey/dTx ...
   }
}
kalmanF.setEntry(9, 9, 1); 	// dTx/dTx
kalmanF.setEntry(10, 10, 1); 	// dTy/dTy
kalmanF.setEntry(11, 11, 1); 	// dTz/dTz

If you have any suggestion they are more than welcome!

Thank you very much for your help.

B.

Have you tried using Orekit’s KalmanEstimator ? :wink:

Or at least have you looked how it’s done there ? The computation of the state transition matrix is done in the method KalmanModel#getErrorStateTransitionMatrix.

If you’ve set the OrbitType of your propagator to equinoctial then your derivatives in mapper.getStateJacobian(nextStepState, jacobian); and mapper.getParametersJacobian(nextStepState, jacobianParam); should be with respect to equinoctial parameters.

So your code looks ok to me. Where exactly do you get an error ?

I don’t get an error per se but I get weird values… see the first line below. I would expect quite low values and potentially higher values on the last line (higher dependency of longitude wrt the other parameters):

jacobian(1,:): 0.9999929232378609 | -11179.473292758054 | 188284.5748993603 | -4.410319010145031 | -78.73060850174079 | 379122.96532085486 
jacobian(2,:): -5.037507530591889E-12 | 1.0000093190774522 | 0.002245084641970569 | -9.850898632262402E-8 | -1.3923259590635902E-6 | 0.0067217660739083765 
jacobian(3,:): 8.001539429059155E-11 | -0.0022902433206608034 | 1.0000051170604998 | 1.018971333109224E-9 | -1.28362688457843E-7 | 4.230109765430311E-4 
jacobian(4,:): -1.210830110639109E-16 | 1.0153031218483577E-8 | 6.209851278636877E-10 | 0.9999976080898223 | -0.0022226301236119364 | 6.447792758188713E-10 
jacobian(5,:): -7.349109716878387E-18 | 6.275543681919286E-10 | 3.840289417930982E-11 | 0.002245104871814292 | 0.9999974329999803 | -1.01699326272539E-8 
jacobian(6,:): -5.3448826054231894E-11 | 0.002263795142480988 | 1.3450309191657896E-4 | 4.614463776866222E-7 | 1.5249828449193317E-8 | 0.999994594142847``

Hi @benoist,

What kind of force models did you add to your propagator ?

I agree with you, but if you have only the central body gravity as a force model (pure Keplerian context).

Example: Equinoctial (mean anomaly) STM (obtained with derivative structures) for an elliptical orbit (a = 15000km, e = 0.125) after an hour of propagation.

  1. With only central body gravity:
	 1.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00  
	 0.000000000000e+00   1.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00  
	 0.000000000000e+00   0.000000000000e+00   1.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00  
	 0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   1.000000000000e+00   0.000000000000e+00   0.000000000000e+00  
	 0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   1.000000000000e+00   0.000000000000e+00  
	-1.237184578301e-07  -1.325762416515e-14   2.442490654175e-15   0.000000000000e+00   0.000000000000e+00   1.000000000000e+00
  1. With a 20x20 Earth gravity potential added:
	 9.998719140385e-01   3.326050735226e+04  -1.993246493362e+03   1.309131780669e+04  -8.681147887369e+03   1.446851409025e+04  
	-3.153414290762e-11   9.997243930087e-01   3.895823310557e-04   1.585024642661e-04   6.013024534925e-04   8.130568846682e-05  
	-5.161113467635e-11   6.209358038102e-04   1.000557228680e+00   7.967632300848e-04   7.677996657246e-05   1.070777006354e-03  
	-3.455127225711e-11  -9.842738774145e-05   3.946960610649e-04   1.000000564658e+00  -3.391442457066e-04   1.314364659125e-04  
	 1.226147172524e-11   1.758417653152e-04  -9.046224373569e-05  -6.844086571906e-05   1.000254686174e+00   6.758629548939e-05  
	-1.236655265848e-07  -2.389731214895e-03  -7.607708945321e-04  -2.535056208532e-03  -1.371055164247e-04   9.979464548385e-01

So your STM looks good to me.
If you’re still unconvinced and think there is a bug in your code we’ll need more context (orbit, force model, propagation configuration etc.) to try and help you.

Cheers,
Maxime