predictPositionVelocity returns wrong partials for nearly circular orbits

Orekit development
1

Hello everybody,

I am trying to compute the partials of slightly shifted PVCoordinates w.r.t. to the initial position and velocity. For this, I am using coordinates of type Gradient.

However, I noticed that the implementation in KeplerianMotionCartesianUtility.predictPositionVelocity, which is used by FieldOrbit<T>.shiftedBy(T dt) and in the Keplerian propagator, returns a wrong result if the input FieldPVCoordinates<T> define a circular (or almost circular) orbit.

To investigate this further, I computed the partials of the shifted PVCoordinates using different methods, namely:

  1. FieldAbsolutePVCoordinates<T>.shiftedBy(T dt), which uses a constant acceleration
  2. KeplerianMotionCartesianUtility.predictPositionVelocity directly
  3. FieldCartesianOrbit<T>.shiftedBy(T dt), which uses (2)
  4. FieldKeplerianPropagator<T>, which also uses (2)
  5. FieldNumericalProagator<T>, which integrates the Keplerian dynamics numerically

I am considering a very small time shift, such that (1) is a good approximation of the other approaches if the initial acceleration is set equal to the nominal Keplerian acceleration.

For nearly circular orbits, (1) and (5) produce consistent results, while (2) to (4) fail to compute the correct partials.

For more eccentric orbits, all approaches return approximately the same results (which is the expected behavior).

These are the results I obtained in two different cases (for the second case I just perturb the velocity vector to increase the eccentricity of the first orbit)

Initial Keplerian parameters: {a: 6797121.104302704; e: 2.2204485879284054E-16; i: 100.39951330616847; pa: 89.91335946442716; raan: 26.6094141509654; v: 0.08664053557284496;}

Partials of final position w.r.t. initial position and velocity

constant acceleration:
dx/dX0[1.0000168994862317, 4.5283252297977215E-6, -1.468012823409825E-5, -6.004370269625362E-4, -2.7189281896099027E-9, 8.814343594433244E-9]
dy/dX0[8.466108202097147E-6, 1.0000022685477443, -7.354280026392985E-6, -5.083278934008411E-9, -6.004282421537878E-4, 4.415707411311126E-9]
dz/dX0[-1.2576069630502057E-11, -3.3698381488348877E-12, 1.0000000000109246, 7.55101025161933E-15, 2.0233414060012877E-15, -6.004268800633028E-4]

predicted PV (analytical Keplerian motion):
dx/dX0[37.8013609081357, -1.6448618582351173, -700.6896803759753, 142874.88808466596, -331305.3883924205, 8.81419456915055E-9]
dy/dX0[41.64139772915594, -50.59850875272008, -316.82260566518426, -204016.05429537073, -159168.66324828798, 4.415611375196497E-9]
dz/dX0[4.320477920155774, -8.42424405988989, -10.427034967150316, -35713.15754554635, -18633.058849253375, -6.004268800944921E-4]

shifted orbit (analytical Keplerian motion):
dx/dX0[37.8013609081357, -1.6448618582351173, -700.6896803759753, 142874.88808466596, -331305.3883924205, 8.81419456915055E-9]
dy/dX0[41.64139772915594, -50.59850875272008, -316.82260566518426, -204016.05429537073, -159168.66324828798, 4.415611375196497E-9]
dz/dX0[4.320477920155774, -8.42424405988989, -10.427034967150316, -35713.15754554635, -18633.058849253375, -6.004268800944921E-4]

analytical Keplerian propagator:
dx/dX0[37.8013609081357, -1.6448618582351173, -700.6896803759753, 142874.88808466596, -331305.3883924205, 8.81419456915055E-9]
dy/dX0[41.64139772915594, -50.59850875272008, -316.82260566518426, -204016.05429537073, -159168.66324828798, 4.415611375196497E-9]
dz/dX0[4.320477920155774, -8.42424405988989, -10.427034967150316, -35713.15754554635, -18633.058849253375, -6.004268800944921E-4]

numerical propagator:
dx/dX0[1.000016899486007, 4.528325220798379E-6, -1.4680128179234142E-5, -6.004370272592041E-4, -2.7189912543690298E-9, 8.81443895650591E-9]
dy/dX0[8.466108193105892E-6, 1.000002268547533, -7.354280135291047E-6, -5.083336418465478E-9, -6.004282424783014E-4, 4.415767307364149E-9]
dz/dX0[-1.252131731632744E-11, -3.4786895586336186E-12, 1.0000000000113594, 3.410605131648481E-13, 1.1368683772161603E-13, -6.004268803110335E-4]

Initial Keplerian parameters: {a: 7849161.160914415; e: 0.1354594254432164; i: 100.33817753639391; pa: 81.59899074987786; raan: 32.90381183966837; v: 9.535005927822688;}

Partials of final position w.r.t. initial position and velocity

constant acceleration:
dx/dX0[1.0000168994862317, 4.5283252297977215E-6, -1.468012823409825E-5, -6.004370269625362E-4, -2.7189281896099027E-9, 8.814343594433244E-9]
dy/dX0[1.0934365593905585E-5, 1.0000029299330708, -9.49838870103491E-6, -6.5652870189490684E-9, -6.004286392673159E-4, 5.703087893329081E-9]
dz/dX0[-1.2576069630502057E-11, -3.3698381488348877E-12, 1.0000000000109246, 7.55101025161933E-15, 2.0233414060012877E-15, -6.004268800633028E-4]

predicted PV (analytical Keplerian motion):
dx/dX0[1.0000168994860035, 4.52832521805049E-6, -1.4680128180495782E-5, -6.00437026792361E-4, -2.7187660806215795E-9, 8.817182268901528E-9]
dy/dX0[1.0934365582690384E-5, 1.0000029299328588, -9.498388810428104E-6, -6.564782477943558E-9, -6.004286392050358E-4, 5.704714862079537E-9]
dz/dX0[-1.2521245326361835E-11, -3.4787149276699216E-12, 1.00000000001136, -1.0051532092441608E-13, 5.221940995658252E-14, -6.004268803226246E-4]

shifted orbit (analytical Keplerian motion):
dx/dX0[1.0000168994860035, 4.52832521805049E-6, -1.4680128180495782E-5, -6.00437026792361E-4, -2.7187660806215795E-9, 8.817182268901528E-9]
dy/dX0[1.0934365582690384E-5, 1.0000029299328588, -9.498388810428104E-6, -6.564782477943558E-9, -6.004286392050358E-4, 5.704714862079537E-9]
dz/dX0[-1.2521245326361835E-11, -3.4787149276699216E-12, 1.00000000001136, -1.0051532092441608E-13, 5.221940995658252E-14, -6.004268803226246E-4]

analytical Keplerian propagator:
dx/dX0[1.0000168994860035, 4.52832521805049E-6, -1.4680128180495782E-5, -6.00437026792361E-4, -2.7187660806215795E-9, 8.817182268901528E-9]
dy/dX0[1.0934365582690384E-5, 1.0000029299328588, -9.498388810428104E-6, -6.564782477943558E-9, -6.004286392050358E-4, 5.704714862079537E-9]
dz/dX0[-1.2521245326361835E-11, -3.4787149276699216E-12, 1.00000000001136, -1.0051532092441608E-13, 5.221940995658252E-14, -6.004268803226246E-4]

numerical propagator:
dx/dX0[1.0000168994860075, 4.528325220271023E-6, -1.4680128178568008E-5, -6.004370274865778E-4, -2.7189628326595994E-9, 8.814375007659692E-9]
dy/dX0[1.0934365584802652E-5, 1.0000029299328599, -9.498388809092706E-6, -6.565585408679908E-9, -6.004286401548597E-4, 5.703114425159583E-9]
dz/dX0[-1.252176140553729E-11, -3.478661803058003E-12, 1.000000000011359, 1.1368683772161603E-12, 9.094947017729282E-13, -6.004268807373592E-4]

I would guess the error comes from the solution to the Kepler equation, but I haven’t looked into it yet. Do you have any hint?

This is the code I used to produce the above results:
KeplerianMotionDerivativesTest.java (6.8 KB)

Thank you,
Alberto

2

Hi Alberto,

this might be a direct consequence of this issue.

There is also problems with hyperbolic orbits so I think in general there’s a need for more robustness in Orekit’s Keplerian propagation with Cartesian coordinates for some particular cases.

Cheers,
Romain.

3

Yes it is likely the case. Thank you for pointing this out!

Cheers,
Alberto

4

Contributions welcomed!

Note that for (4), if you define your orbit with equinoctial elements, you shouldn’t have any problems.

Cheers
Romain.

5

Hi Romain,

as you mentioned, it is indeed a limitation of the method used to compute the shifted PVCoordinates.

In particular for #1086, atan2(y, x) is not even defined in (0, 0). Other NaNs also appear in the computation of e few lines below, as sqrt(x) is not differentiable in 0.
DerivativeStructure is indeed also affected.

I did some research on alternative approaches, and I came up with two possible solutions:

  1. The method described in [1] expresses the Lagrange coefficients in terms of change in the eccentric anomaly. This avoids having to compute quantities which are undefined or not differentiable if the orbit is circular. There is also an equivalent method for hyperbolic orbits.
  2. The method described in [2] which uses the universal anomaly. The advantage here is that it does not require branching for different types of orbits.

I implemented prototypes for both, and tested them on my problem and on issues 1086 and 1642. They both succeed in all cases.
Code is available in this repository, in particular in the julia/ folder.

As these are quite substantial changes, I leave open for discussion if and which method should be implemented in Orekit. If there is consensus, I will be happy to contribute!

Important note:

Both methods require to solve a nonlinear equation in the form f(x,p)=0.

  1. to solve for the change in eccentric anomaly given the change in mean anomaly
  2. to solve for the universal anomaly

In my code I use a bare-bone Newton-Raphson algotihm, which does its job perfectly. However, I am aware that Orekit implements very specialized solvers to find the eccentric anomaly given the mean one…

Moreover, extra-care must be taken when using Field elements to ensure that all the partials, and not only the scalar part, are converged. This is easy for Gradient, but less for DerivativeStructure at high orders.

Cheers,
Alberto

[1] R. H. Battin, An Introduction to the Mathematics and Methods of Astrodynamics, Revised Edition, 1999, pp. 162-170.
[2] H. D. Curtis, Orbital Mechanics for Engineering Students, 4th ed., 2021, pp. 167–177.

1 Like
6

Hi Alberto,

thank you for this.
The thing is that this code had already been modified a few years ago to better handle some other special cases, in particular near hyperbolic (see issue 1015).
And it looks as tho a generalized Kepler equation like in Battin was removed.

Given that the circular orbit case problem is blocking for issue #1774, I’ve merged a 5 lines modification to predictPositionVelocity fixing the “symptoms” of 1086. This way we’re sure we’ll have both fixed for Orekit 13.1. The MR for 1774 is here and is gonna need some test tolerance modifications.

Now, in the longer term, it could be beneficial to revamp this Cartesian, analytical propagation, for instance like you suggest with universal variables. We would just need to make sure there is no regression and all special cases work well. This could be done for the next major release, as it is likely to introduce some numerical noise w.r.t. the current results.

Cheers,
Romain.

7

Thank you Romain,

I quickly tested issue 1015, and these are the results I obtain with the two algorithms described above (universal is Curtis, elliptic is Battin):

Test case: Orekit issue 1015

Orbit type: elliptic
Specific energy: -238.82842582464218 m^2/s^2

Solving universal Kepler equation...
Newton-Raphson converged in 11 iterations. Error: 0.0.
...equation solved! Error: 0.0.

Solving elliptic Kepler equation...
Newton-Raphson converged in 6 iterations. Error: 4.36208562054358e-16.
...equation solved! Error: 4.96970112021859e-21.

If I compare the specific energy of the initial state with that of the final one, the absolute errors are 4\times 10^{-4}\, m^2/s^2 for the universal variable approach, and 2.5\times 10^{-5} m^2/s^2 for the other approach.

Great, I’ll test it on my problem.

Agreed. I will try to dedicate some spare time to test these new algorithms against the Orekit test suite, plus all mentioned issues.

I took a quick look at what was implemented before MR 335, and to me it seems yet another method.

Cheers,
Alberto

8

I rebuilt Orekit with the changes you merged in !850, and unfortunately it does not solve the issue at the very beginning of this discussion.

The threshold you set, Precision.SAFE_MIN, is well below the eccentricity of my test orbit, which is about e=2.22e-16. So I still obtain the same results.

Moreover, I tested issue 1086 with \mu=1, initial state equal to \boldsymbol{x}_0 = [1,0,0,0,1,0]^T, and propagation duration \Delta t=\pi/4. These are the results I obtained with Orekit (using Gradient):

KeplerianMotionCartesianUtility.predictPositionVelocity():

x:  0.7071067811865475
y:  0.7071067811865476
z:  0.0
vx: -0.7071067811865476
vy: 0.7071067811865475
vz: 0.0

dx/(dX0, dt): [2.3731878829959347, -0.7071067811865476, 0.0, 0.0, 1.6660811018093873, 0.0, -0.7071067811865476]
dy/(dX0, dt): [-0.9589743206228394, 0.7071067811865475, 0.0, 0.0, -1.666081101809387, 0.0, 0.7071067811865475]
dz/(dX0, dt): [0.0, 0.0, 0.7071067811865475, 0.0, 0.0, 0.7071067811865476, 0.0]
dx/(dX0, dt): [1.666081101809387, 0.0, 0.0, 0.7071067811865475, 0.9589743206228394, 0.0, -0.7071067811865475]
dy/(dX0, dt): [1.6660811018093873, 0.0, 0.0, 0.7071067811865476, 2.3731878829959347, 0.0, -0.7071067811865476]
dz/(dX0, dt): [0.0, 0.0, -0.7071067811865476, 0.0, 0.0, 0.7071067811865475, 0.0]

Numerical propagator (Keplerian dynamics):

x:  0.7071067811865468
y:  0.7071067811865482
z:  0.0
vx: -0.7071067811865482
vy: 0.7071067811865468
vz: 0.0

dx/(dX0, dt): [1.580294664182335, 0.2071067811866306, 0.0, 0.914213562373179, 0.08029466418218623, 0.0, -0.7071067811865491]
dy/(dX0, dt): [0.24813246056385752, 0.7928932188133698, 0.0, 0.08578643762682292, 0.7481324605640052, 0.0, 0.7071067811865477]
dz/(dX0, dt): [0.0, 0.0, 0.7071067811865468, 0.0, 0.0, 0.7071067811865482, 0.0]
dx/(dX0, dt): [1.3731878829957873, 0.7071067811866312, 0.0, 1.414213562373178, 0.37318788299563954, 0.0, -0.7071067811865477]
dy/(dX0, dt): [0.9589743206226936, -0.2928932188133694, 0.0, 0.41421356237317897, 0.9589743206225438, 0.0, -0.7071067811865491]
dz/(dX0, dt): [0.0, 0.0, -0.7071067811865482, 0.0, 0.0, 0.7071067811865468, 0.0]

As you can see, the partials returned by predictPositionVelocity are wrong.

Those obtained with numerical propagation do instead coincide with the analytical ones I obtained with the new methods:

6Ă—7 Matrix{Float64}:
 1.58029    0.207107   0.0       0.914214   0.0802947  0.0       -0.707107
 0.248132   0.792893   0.0       0.0857864  0.748132   0.0        0.707107
 0.0        0.0        0.707107  0.0        0.0        0.707107   0.0
 1.37319    0.707107   0.0       1.41421    0.373188   0.0       -0.707107
 0.958974  -0.292893   0.0       0.414214   0.958974   0.0       -0.707107
 0.0        0.0       -0.707107  0.0        0.0        0.707107   0.0

Cheers,
Alberto

9

Hi Alberto,

thanks for your analysis. I wil have another look.
I must have rushed things too much, my bad.
I think #1086 can be done for a patch, whilst 1774 can wait for Orekit 13.2 or 14.0.

Cheers,
Romain.

10

I have been finally able to dedicate some free time to this.

I have implemented both alternative algorithms, the one that uses the difference in eccentric anomaly, and the one that uses the universal variables, in this branch.

They both seem to solve all mentioned issues.

Cheers,

Alberto