State transition matrix parametrization

MarcoTosato · July 21, 2022, 9:26am

Hello,

I am propagating the covariance matrix making use of the state transition matrix derived from the MatricesHarvester.
I just would like to make sure that the following is true:

If KeplerianPropagator / EcksteinHechlerPropagator → the state transition matrix is in cartesian parametrization;
if NumericalPropagator → the state transition matrix parametrization is function of the setOrbitType and setPositionAngleType methods;
if DSSTPropagator → the state transition matrix is in equinoctial / mean anomaly parametrization.

Since I need to implement a general purpose code where the input covariance may be expressed in any parametrization, do I need therefore to transform the input covariance in the correct parametrization (depending on the propagator, through the getJacobianWrt___ methods), propagate it, and then transform it back? Is there any other simpler way to do this?

Thank you!

bcazabonne · July 21, 2022, 9:45am

Hi @MarcoTosato

We recently implemented a class for that. It is currently available in the development version of Orekit and will be officially released in the 11.3 version.
This class is StateCovarianceMatrixProvider.java. It can be used with different orbit propagators (DSST, numerical, Eckstein-Hechler, Brouwer-Lyddane, SGP4, and Keplerian) and it can provide the covariance in different frames and orbit types.

You can find examples in the test class of this function: StateCovarianceMatrixProviderTest.java

Best regards,
Bryan

MarcoTosato · July 21, 2022, 10:02am

Hi @bcazabonne and thank you!

I’ll make use of that class then, as soon as it is available. Just for my understanding, is anyway true what I said?

Thank you.

bcazabonne · July 21, 2022, 12:03pm

True

By default it is equinoctial / true anomaly but yes you can update it using setOrbitType(OrbitType) and setPositionAngleType(PositionAngle) methods.

True

When computing the propagated covariance, the initial one must be indeed in the same orbit type than the state transition matrix. For the orbit type transformation, you can look the implementations of the changeCovarianceType method in StateCovarianceMatrixProvider.
So, if your state transition matrix is computed in cartesian elements but your input covariance uses keplerian elements, you need:

Transform the input covariance Cov_t0 in cartesian
Compute the propagated covariance using Cov_tk = STM * Cov_t0 * STM^T
Transform the propagated covariance Cov_tk in keplerian

Bryan

MarcoTosato · July 21, 2022, 12:09pm

Thank you!

MarcoTosato · July 25, 2022, 11:29am

Hi,
I just have a last question…
I am trying to propagate using a DSST propagator. As soon as I invoke the method

MatricesHarvester matricesHarvester = dsstPropagator.setupMatricesComputation("stm", null, null);

the propagation does not work anymore.
Am I not supposed to extract the MatricesHarvester (needed for stm retrieval) in this way when propagating with a DSST propagator?

Thank you

bcazabonne · July 25, 2022, 12:10pm

Hi,

Could you try this:

    final DSSTHarvester harvester = (DSSTHarvester) propagator.setupMatricesComputation(stmName, null, null);
    harvester.initializeFieldShortPeriodTerms(initialState);

Best regards,
Bryan

MarcoTosato · July 25, 2022, 1:21pm

Thank you @bcazabonne,

this actually works, but it really takes too long if setting some dsst force models besides the newtonian attraction. Isn’t there any other alternative to this initializeFieldShortPeriodTerms?

Thank you!

bcazabonne · July 25, 2022, 1:40pm

What is your integration step?

Unfortunately no… I’m thinking about a solution.

Bryan

MarcoTosato · July 25, 2022, 2:05pm

I am using a classical Runge Kutta with step 0.5 and I am propagating the orbit for more or less one hour.
If there is no alternative, I’ll just use a numerical propagator then. Thank you!

Marco

bcazabonne · July 25, 2022, 2:07pm

Ok, using DSST you can increase the integration step because in DSST theory, the integrated elements are the mean orbital elements. Typically, the integration step can be equal to several orbital periods.

Bryan

Serrof · July 27, 2022, 3:17pm

Hi,

I am confused by your answer here @bcazabonne. In my understanding, the variables propagated by an instance of KeplerianPropagator are the ones of the underlying initial state (via its Orbit and the corresponding shiftedBy method). So I assumed that the STM would be using the same coordinates. But you’re saying that it is always Cartesian?

Cheers,
Romain.

bcazabonne · July 27, 2022, 4:50pm

Hi Romain,

For analytical orbit propagators the type of the STM is defined in the AbstractAnalyticalMatricesHarvester class. If you look at the implementation of the setReferenceState method you will see that we use the cartesian derivatives.

We decided to use cartesian parameters for two reasons : (1) consistency between propagators (Keplerian propagator uses keplerian elements, Eckstein-Hechler uses the cartesian ones and Brouwer-Lyddane uses keplerian) and (2) to avoid possible singularities.

Now, based on your comment, what could be interesting is to add a method getOrbitType() and getPositionAngle() in the MatricesHarvester interface to help users knowing the type used.

Best regards,
Bryan

Serrof · July 28, 2022, 8:13am

Hi Bryan,

thank you for the clarification. I would indeed appreciate some getters as I find this behaviour a bit counter-intuitive.
However, you confused me even more by stating that a Keplerian propagator uses Keplerian elements!!!
Again, I thought it used the shiftedBy method of the underlying Orbit, so if it is a CartesianOrbit, it uses shiftPVElliptic for non-hyperbolic cases.

Best,
Romain.

bcazabonne · July 28, 2022, 8:25am

The main problem is that we don’t have getOrbitType() and/or getPositionAngle() methods for analytical orbit propagators. So we had to make a choice and we opted for the Cartesian for the two previous reasons.

I’ll create an issue for adding the getter. It will only be applicable for the next major version because MatricesHarvester is an interface and having default implementations is not a good idea for these methods.

Right! I thought the return type of the propgateOrbit() method was a KeplerianOrbit. But it is indeed an Orbit, so there is no imposed type as for the Eckstein-Hechler or Brouwer-Lyddane propagators which return a CartesianOrbit or a KeplerianOrbit respectively.
My apologies.

Best regards,
Bryan

Serrof · July 28, 2022, 8:31am

No worries. And to be clear, I appreciate that the architecture of Orekit is not a trivial problem and that you guys are tackling it as best as you can. This is a great library! Thanks for the amazing work.

Romain.

MarcoTosato · September 9, 2022, 9:29am

Hi @bcazabonne, I exploit this post to add a small point.
I tried to replicate the test “testWithAnalyticalPropagator()” in src/test/java/org/orekit/propagation/StateCovarianceMatrixProviderTest.java · develop · Orekit / Orekit · GitLab. The code I wrote is the following:

        // Initial state
        AbsoluteDate initEpoch = new AbsoluteDate("2016-02-13T16:00:00.000", TimeScalesFactory.getUTC());
        Orbit initialOrbit = new CartesianOrbit(new PVCoordinates(new Vector3D(7526993.581890527, -9646310.10026971, 1464110.4928112086),
                                                                  new Vector3D(3033.79456099698, 1715.265069098717, -4447.658745923895)),
                                                FramesFactory.getEME2000(),
                                                initEpoch,
                                                Constants.WGS84_EARTH_MU);

        double[][] initCov = new double[][] {{8.651816029e+01,   5.689987127e+01, -2.763870764e+01, -2.435617201e-02,  2.058274137e-02, -5.872883051e-03},
                                             {5.689987127e+01,   7.070624321e+01,  1.367120909e+01, -6.112622013e-03,  7.623626008e-03, -1.239413190e-02},
                                             {-2.763870764e+01,  1.367120909e+01,  1.811858898e+02,  3.143798992e-02, -4.963106559e-02, -7.420114385e-04},
                                             {-2.435617201e-02, -6.112622013e-03,  3.143798992e-02,  4.657077389e-05,  1.469943634e-05,  3.328475593e-05},
                                             {2.058274137e-02,   7.623626008e-03, -4.963106559e-02,  1.469943634e-05,  3.950715934e-05,  2.516044258e-05},
                                             {-5.872883051e-03, -1.239413190e-02, -7.420114385e-04,  3.328475593e-05,  2.516044258e-05,  3.547466120e-05}};
        RealMatrix initCovReal = MatrixUtils.createRealMatrix(initCov);

        // Propagator
        EcksteinHechlerPropagator prop = new EcksteinHechlerPropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(6, 0));
        MatricesHarvester harvester = prop.setupMatricesComputation("stmName", null, null);

        // Propagate and compute covariance
        SpacecraftState propagatedState = prop.propagate(initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY));

        RealMatrix dYdY0 = harvester.getStateTransitionMatrix(propagatedState);
        RealMatrix propCov = dYdY0.multiply(initCovReal.multiplyTransposed(dYdY0));

        System.out.println(propCov.getEntry(0, 0));

Running it, I get that the first element of the propagated covariance matrix is 577.3144349040394, which is not within the threshold of 5.0e-4 to your reference 5.770543135e+02. Is this just a numerical error or did I make a mistake in my code?
Thank you.

MaximeJ · September 9, 2022, 12:57pm

Hi @MarcoTosato,

I take the liberty to reply in place of @bcazabonne
If you take a close look at method StateCovarianceMatrixProviderTest#compareCovariance you’ll see that the threshold is relative.

Assert.assertEquals(reference.getEntry(row, column), computed.getEntry(row, column), FastMath.abs(threshold * reference.getEntry(row, column)));

So in the case of the first element of the matrix the absolute threshold is 5.10^{-4} \times 577.3144349040394 = 0.28866

Hope this helps,
Maxime

bcazabonne · September 9, 2022, 4:03pm

Thank you @MaximeJ