I am trying to perform an orbit propagation using the built-in PartialDerivativesEquations class to compute the State Transition Matrix wrt some initial conditions. My intention is to propagate a given covariance ellipsoid and see how it evolves after a certain period of time. This is fairly straigh-forward when working in Cartesian parameters, but would it be possible to do using Keplerian parameters? My goal is to check if the results using Keplerian would be in any way better (given certain statistical figures of merit) than the results using Cartesian, due to a (possible) better linearization of the dynamics.

I have seen an old (circa Aug. 2017) post regarding this subject, but no clear solution was given. For now, what I have tried is to set the OrbitType to â€śKeplerianâ€ť for the NumericalPropagator, without much success (after closer inspection, the â€ścomputeDerivativesâ€ť method always uses Cartesian anyway, just like the link above suggests).

There has been, indeed, several discussions around this subject.
The latest one on the forum gives a solution using Taylor algebra that gives good results.

PartialDerivativesEquations and JacobianMapper classes should work with Keplerian parameters (or any orbit type) but the results do not seem to be numerically stable.
There is an issue on the Orekit forge about that (see #351).

The derivatives are inded computed in Cartesian elements in PartialDerivativesEquations class.
But if you look closely in the JacobiansMapper class here and here youâ€™ll see that the conversion from userâ€™s orbit type to Cartesian (and back) are performed before and after propagation. Thatâ€™s why Iâ€™m saying that it should work with Keplerian parameters.

However, until a proper fix or an explanation is given for the numerical discrepancies we observe using PartialDerivativesEquationâ€¦ I suggest you use the Taylor algebra example given above.