Orbital covariance propagation over one orbital revolution

Hello,
I am approaching for the first time the theory of covariance propagation and I have some doubts about the results I am obtaining when propagating numerically both the orbit and its covariance matrix.
Following what suggested in Covariance matrix before last measurements, I implemented this piece of code:

public class test {
	
    public static void main(String[] args) {
        
        // Initial state
        TimeScale utc = TimeScalesFactory.getUTC();
        AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, utc);
        Vector3D position = new Vector3D(-605000.7922166, -5870000.2295111, 3493000.0531990);
        Vector3D velocity = new Vector3D(-1568.254290, -3702.348910, -6479.483950);
        PVCoordinates cartCoord = new PVCoordinates(position, velocity);
        CartesianOrbit initialOrbit = new CartesianOrbit(cartCoord, FramesFactory.getGCRF(), initialDate, Constants.IERS2010_EARTH_MU);
        double period = initialOrbit.getKeplerianPeriod();
        
        // Initial covariance
        double[][] initialCovariance = {{1.000000e+02, 1.000000e-02, 1.000000e-02, 1.000000e-04, 1.000000e-04, 1.000000e-04},
				  						{1.000000e-02, 1.000000e+02, 1.000000e-02, 1.000000e-04, 1.000000e-04, 1.000000e-04},
				  						{1.000000e-02, 1.000000e-02, 1.000000e+02, 1.000000e-04, 1.000000e-04, 1.000000e-04},
				  						{1.000000e-04, 1.000000e-04, 1.000000e-04, 1.000000e-04, 1.000000e-06, 1.000000e-06},
				  						{1.000000e-04, 1.000000e-04, 1.000000e-04, 1.000000e-06, 1.000000e-04, 1.000000e-06}, 
				  						{1.000000e-04, 1.000000e-04, 1.000000e-04, 1.000000e-06, 1.000000e-06, 1.000000e-04}};
        RealMatrix initCovariance = new Array2DRowRealMatrix(initialCovariance);
        
        // Numerical propagator
        NumericalPropagator prop = new NumericalPropagator(new ClassicalRungeKuttaIntegrator(1));
        SpacecraftState initialSpacecraftState = new SpacecraftState(initialOrbit);
        prop.setOrbitType(OrbitType.CARTESIAN);
        prop.setInitialState(initialSpacecraftState);
        
        PartialDerivativesEquations PDE = new PartialDerivativesEquations("partial-derivatives", prop);
        SpacecraftState initialSpacecraftStateWithDerivatives = PDE.setInitialJacobians(initialSpacecraftState);
        prop.resetInitialState(initialSpacecraftStateWithDerivatives);
        JacobiansMapper jacobianMapper = PDE.getMapper();
        
        prop.setStepHandler(new PropagationStepHandler(jacobianMapper, initCovariance));
        prop.propagate(initialDate.shiftedBy(period));
        
    }
    
    // Step handler
    private static class PropagationStepHandler implements OrekitStepHandler {
    	
    	private JacobiansMapper jacobianMapper;
    	private RealMatrix initCovariance;

        public PropagationStepHandler(JacobiansMapper jacobianMapper, RealMatrix initCovariance) {

        	this.jacobianMapper = jacobianMapper;
        	this.initCovariance = initCovariance;
        	
        }

		@Override
		public void handleStep(OrekitStepInterpolator interpolator) {
			
			SpacecraftState currentSpacecraftState = interpolator.getCurrentState();
			
			// Current state transition matrix
			double[][] stm = new double[6][6];
			jacobianMapper.getStateJacobian(currentSpacecraftState, stm);
			RealMatrix dYdY0 = new Array2DRowRealMatrix(stm, false);
			
			// Current covariance
			RealMatrix newCovariance = dYdY0.multiply(initCovariance.multiplyTransposed(dYdY0));
			System.out.println(newCovariance.getEntry(0, 0));
			
		}
        
    }
        
}

As you can see, I am running the propagation over one orbital revolution. I was therefore expecting to obtain the final covariance matrix equal to the initial one since the cartesian elements at the end of the propagation are exactly those of the initial state (I was expecting the dYdY0 matrix to be the identity matrix). Instead I get that, for example, the element [0][0] of the final covariance matrix is equal to 2774.37 instead of 100.
Did I implement something wrong in the code or am I just missing something from the theory of covariance propagation?

Thank you in advance!
Kind regards,

Claudio Toquinho Campana

Hi @toquika

What you obtain is expected. In addition, your code looks good!
The state transition matrix is equal to the identity matrix at the initial epoch. However, because you are propagating your orbit, the state transition matrix is not equal to the identity matrix after one orbital period. Indeed, it takes into account the impact of the Newtonian attraction.

Generally speaking, based on Gaussian assumption, the covariance matrix represents the uncertainty of the orbital state. Propagating an orbit increases the uncertainty. As a result, the covariance matrix increases too.

Best regards,
Bryan

When we deal with covariance, we state that their are some uncertainties in the initial data
and that will affect the orbit. Covariance propagation is a way to quantify how these
uncertainties evolve (usually, they grow up).

Under these assumption, even the notion of “over one orbital revolution” is cautious,
because even the orbital period is uncertain. So if you have an error on velocity for example,
then the partial derivatives of position Y with respect to position Y0 will change because
velocity will change the time it takes to travel “over one orbital revolution”. Everything is
coupled and mixed together.

Thank you both @bcazabonne and @luc for your kind answers. I have an additional question related to this topic. I think I am being a bit confused about the difference between the state transition matrix and the Jacobian matrix. In the link I provided you, it is said that the state transition matrix is retrieved by using the JacobiansMapper.getStateJacobian method. Why is it called “state transition matrix” since it is indeed the Jacobian? What would be the theoretical difference between the state transition matrix and the Jacobian one in this case?
I am asking this because I saw that starting from Orekit11.1 PartialDerivativesEquations class will be deprecated and therefore I should modify my code to make use of the MatricesHarvester interface.
Thank you!

Hi @toquika

The covariance matrix represents the uncertainty of the orbital state.
The state transition matrix represents the partial derivatives of the orbital parameters with respect to the initial ones. This matrix highlights the dynamical evolution of the parameters.

The previous architecture (i.e., based on JacobianMapper and PartialDerivativesEquations classes) was not user-friendly and was ambiguous. Using the JacobianMapper-based architecture, the getStateJacobian() method gives the state transition matrix. It is called jacobian because it is a partial derivatives matrix.
This method has been replace in the 11.1 version by MatricesHarvester.getStateTransitionMatrix(). It gives the same result. Only the signatures change between the two implementations. You can find an example on how to use this new architecture in the CovariancePropagation tutorial.

Best regards,
Bryan

Thank you!

Best regards,
Toquinho