Question on State Transition Matrices & Covariance

Orekit usage
1

Hello everyone!

I have some questions on consistency. I am producing a program in orekit which simulates the process of producing and correlating uncorrelated tracks.

I simulate a satellite orbit and take data when it passes though a radar fov. These “tracks” are on the order of 2-5 minutes long and the radar records the range and angular position every ten seconds while in the fov. The data points have random error associated with them ( 0.015 degrees, 40 m of std dev).

The first, last and middle radar data for each track is then plugged into a herrick-gibbs script I wrote to produce an initial orbital guess. That orbit is then cast into the EquinoctialOrbit type, as that suits well for bls.

To refine the guess I use the orekit batch-leasts-quares estimation package. I take that initial herrick-gibbs equinoctial orbit and use it as the initial guess. The rest of the tracks range and angular measurements are then fed into the bls estimator and a refined orbital guess is produced.

For each track from the same satellite, one would expect the bls estimator to produce relatively similar orbits. To check, I am using the mahalanobis distance like in this paper (Comparison of Covariance Based Track Association Approaches Using Simulated Radar Data | The Journal of the Astronautical Sciences)

To compare two tracks and see if they correlate, I take the covariance and state of the first occurring track and propagate it to the epoch of the second track. This is where my questions lie and where I believe I am having issues.

The propagation of the the epoch 1 state to the second is not an issue. The state vector I use is (a, ex,ey,hx,hy,lv) since those are the orekit elements. To propagate the covariance, however, I need the state transition matrix which I get from the propagation of epoch 1 to epoch 2 shown in the following code block:

  	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);

Then I get the covariance using estimator.getPhysicalCovariance

Then the mahalanobis distance is calculated in the following manor:

Pm(tn) = stm(tn,tm) *Pm( tm) *stm(tn,tm)^T
Pmn(tn) = Pn(tn) + Pm(tn)
z(tn) = xn(tn) - xm(tn), tn>tm
(km,n)^2 = z^T *Pmn(tn)^(-1) *z(tn)

where km,n is the mahalanobis distance. The more related the track the lower this value. tn is the later epoch and tm the earlier. x is the state and n or m mean which track they belong to, and finally stm is the state transition matrix.

Unfortunately, the output I am getting is wrong by orders of magnitude, and points towards an error. I believe that the error arises from the state transition matrix, the covariance or both. From browsing the archived mailing list it seems to me that the stm I’m producing likely might be the cartesian state transition matrix. Additionally, I’ve found both the optimum covariance and the physicalcovariance. I get the feeling the physical covariance may also be for cartesian elements.

Are my hunches correct or is something else happening?

Here I’ve posted the rest of my code, if anyone is interested or thinks something deeper is happening: https://pastebin.com/YfC2v0Cd

Thanks anyone for their help and I hope everyone has a great week!

Sincerely,
Paul

2

Hi Paul,

Last time I looked in to this only Cartesian elements are supported when computing the STM.[1] Bug report is at [2]. Contributions are welcome.

The linked thread suggests a couple different methods, either using the d(Cartesian)/d(Equinoctial) Jacobian on both ends of the propagation, or using a FieldNumericalPropagator with DerivativeStructure.

[1] https://www.orekit.org/mailing-list-archives/orekit-developers/msg00114.html
[2] https://gitlab.orekit.org/orekit/orekit/issues/351

3

@evan.ward Ah I think I understand from your links, but I still think I need some help. From the earlier orbit propagated to the later orbit I should use something like:
> double deqdcart = new Array2DRowRealMatrix(6, 6).getData();

  	orbf.getOrbit().getJacobianWrtCartesian(PositionAngle.MEAN, deqdcart);

where orbf is the earlier orbit propegated to the second epoch.
This gives me the d(Eqiu)/d(cart) at later epoch. Then I take the dYd0 after the

line and multiply to get
stm = (d(eq)/d(cart))dYdY0(d(eq)/d(cart))^-1

correct? That’s what I gather from the links but you say

What do you mean by both ends?

Do you mean I should also get d(eq)/d(cart) before propagation like the following:

        double[][] deqdcart = new Array2DRowRealMatrix(6, 6).getData();
		orbi.getOrbit().getJacobianWrtCartesian(PositionAngle.MEAN, deqdcart);

where orbi is the orbit before being propagated to epoch 2?

Sincerely,
Paul

4

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

5

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.

6

Ah thanks so much I’ll try this out!

Sincerely,
Paul

7

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?

8

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

9

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!

10

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
11

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

12

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
13

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

14

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
15

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

16

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

17

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

18

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.

19

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

20

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.