Information about AD in Orekit

[nick]

Dear all,

I’m trying to learn more on AD (i.e. automatic differentiation) features in Orekit, following a couple of tutorials showing the use of these functionalities (specifically, FieldPropagation and GradientComputation, the only I managed to find about this topic).

I must admit that I’m quite a newbie on AD in Orekit, so I’m trying to get more proficient by also playing around with it and implementing some test cases.

For this purpose, I implemented a couple of test cases (that I’m attaching below for your reference) on which I have few questions.

Regarding the first test case (TestAD1.java (4.7 KB)), I have only a doubt about the full correctness of its results. Specifically, in this test case I’m trying to compute (by means of AD) the time derivative (first order) of the Keplerian elements associated to the final state of a numerical propagation. Clearly, I had to setup some force models in the numerical propagator to make this test meaningful. At the end of the simulation, I’m printing the results and comparing them against the output of dedicated methods from KeplerianOrbit to get aDot, eDot, etc.
The difference between the AD results and those from the mentioned KeplerianOrbit getters is quite small, but still the results are not really the same.
Is this expected? Am I doing something wrong here, or misunderstanding something?

Regarding the second tests case (TestAD2.java (5.0 KB)), here the problem concerns instead derivatives that should be computed in theory, but apparently they aren’t. Specifically, in this test case the simulation scenario is based on a Hohmann transfer from an initial 500x500 km orbit to a final 1000x1000km orbit, considering a pure Keplerian dynamics. I’m performing only the first (impulsive) maneuver to inject into the transfer orbit, propagating for half period of this orbit and then retrieving the final position vector (note that in this example I’m not performing any final maneuver to inject to the target circular orbit). Once I get the final position, I’m trying to compute the gradient of each position component w.r.t. the free variables, which I have previously set to be the delta-v components of the initial maneuver. In principle, I would expect some numeric values for these gradients, but I get NaN everywhere.
Why is this happening? Is there something wrong I’m doing in this case?

Lastly, still about second test case, is it possible to use something like the “FieldImpulseManeuver” but for the case of finite maneuvers (to set, for example, the maneuver thrust and/or direction as free variables)?

As always, many thanks in advance to anyone can provide some feedback about these points and for your time

[Vincent]

Hello @nick

You are already getting the hang of it !

Regarding your first question :

@nick

Is this expected? Am I doing something wrong here, or misunderstanding something?

I’ve run your example and the differences are quite small. I’m not surprised as both process are using different methods to compute the same things:

1. The getADot method in your case uses : computeJacobianTrueWrtCartesianElliptical
2. Meanwhile, the AD goes through other steps

Regarding your other problem

@nick

Why is this happening? Is there something wrong I’m doing in this case?

Good news, there is something wrong ! Your initial orbit should be initialized otherwise :

        DerivativeStructure a = dsf.constant(7e6);
        DerivativeStructure e = dsf.constant(0.0);
        DerivativeStructure i = dsf.constant(FastMath.PI / 4);
        DerivativeStructure pa = dsf.constant(0.0);
        DerivativeStructure raan = dsf.constant(0.0);
        DerivativeStructure tan = dsf.constant(0.0);

You have e=0 and pa=0 which is singular and cause your gradients to go NaN further in the internal code. A simple fix is to convert this initial orbit to a CartesianOrbit or, even better, use CircularOrbit or EquinoctialOrbit right away instead.

I have also noticed one “minor issue”:

       double dt = 1.0;
        propagator.propagate(initialDateDS.shiftedBy(-dt));

These lines are not useful, the propagator initial state is still the same:

Before :FieldSpacecraftState{orbit=Keplerian parameters: {a: 6878136.6; e: 0.0; i: 0.0; pa: 0.0; raan: 0.0; v: 0.0;}, attitude=org.orekit.attitudes.FieldAttitude@6ea6d14e, mass=org.hipparchus.analysis.differentiation.Gradient@c3e3eef0, additional={}, additionalDot={}}
After:  FieldSpacecraftState{orbit=Keplerian parameters: {a: 6878136.6; e: 0.0; i: 0.0; pa: 0.0; raan: 0.0; v: 0.0;}, attitude=org.orekit.attitudes.FieldAttitude@6ea6d14e, mass=org.hipparchus.analysis.differentiation.Gradient@c3e3eef0, additional={}, additionalDot={}}

Cheers,
Vincent

[nick]

Hello @Vincent,

many thanks for your kind reply and the feedback!

Regarding the second test case, I tried to implement the fix you suggested about the orbit, but it seems to work only when converting to a FieldEquinoctialOrbit; indeed, I get numeric values for the gradients only in this case, whereas when I convert to other orbit parametrizations I still get NaN everywhere. Don’t know if this is an issue or something expected for some reason.

Still on the second test case, I’m considering this logic you commented

double dt = 1.0;
        propagator.propagate(initialDateDS.shiftedBy(-dt));

because I have an impulsive maneuver located at the really beginning of my transfer trajectory and, in order to let the propagator to “notice” it, I need to first step back for a little and then move forward, otherwise the maneuver won’t be triggered (but this is another story…).

Lastly, regarding this point

is it possible to use something like the “FieldImpulseManeuver” but for the case of finite maneuvers (to set, for example, the maneuver thrust and/or direction as free variables)?

do you maybe know if there’s something already there in Orekit for this?
Otherwise, if not, is it something rather complex to implement on my own (consider that I’m just starting to get proficient with the Orekit AD )?

[Vincent]

You are welcome !

@nick

indeed, I get numeric values for the gradients only in this case, whereas when I convert to other orbit parametrizations I still get NaN everywhere. Don’t know if this is an issue or something expected for some reason.

That is strange, i tested a quick fix by creating a FieldCartesianOrbit out of your initial FieldKeplerianOrbit and it worked. This should also work using a FieldCircularOrbit. I’ll have to investigate further .

@nick

because I have an impulsive maneuver located at the really beginning of my transfer trajectory and, in order to let the propagator to “notice” it, I need to first step back for a little and then move forward, otherwise the maneuver won’t be triggered (but this is another story…).

Oh i get what you wanted to do, and i also thought that it would have required further fixing. Hopefully it worked right away but if you want to be sure you must reset the propagator state after “stepping back”. You could also specify the whole propagation interval instead of only specifying the target date.

@nick

do you maybe know if there’s something already there in Orekit for this?

I totally forgot to answer this one . You would need to implement your own PropulsionModel. The existing BasicConstantThrustPropulsionModel already handles thrust and flow rate as free variables (because they are ParameterDriver) but uses a constant direction so you would need to setup a parameter driver for that as well. Once this is done, you will be able to create your own Maneuver and add it to your propagator.

Feel free to post a message on the forum if you are stuck !

Cheers,
Vincent

[Serrof]

Hi,

@nick

Lastly, still about second test case, is it possible to use something like the “FieldImpulseManeuver” but for the case of finite maneuvers (to set, for example, the maneuver thrust and/or direction as free variables)?

are you interested in first-order derivatives only?

If so, partial derivatives w.r.t. propagation model are available via the MatricesHarvester of NumericalPropagator (AD is used under the hood but you don’t need to worry about it). You “simply” need to activate beforehand the so-called ParameterDriver of interest in your ForceModel. The class ConstantThrustManeuver for instance already has some on the dates, the propulsion ones, etc. It’s the system used internally by estimators using differential information, so you chan check out the OD tutorials to see how to deal with the drivers.

Cheers,
Romain.

[nick]

Hello @Vincent,

regarding the issue due to the orbit parametrization, just for your reference I’m attaching below an updated version of the second test case, where I’m considering the three additional parametrizations (see lines 48-50), so you just need to use the desired one when supplying the orbit as input to the propagator (at line 59). As I was saying in my previous message, the Equinoctial orbit is the only one that proved to work, whereas all the other result in NaN everywhere in the evaluated derivatives.

TestAD2mod.java (5.3 KB)

@Serrof, I may need also to evaluate second-order derivatives, so not sure whether your suggestion would also apply for this case.

Anyway, I’ll try to follow your suggestions to eventually get a Maneuver with the desired free parameters in it, so that I can later perform my derivatives w.r.t. those parameters.
In case of issues, I will try to reach you out in this forum for some help

Thanks a lot all and have a nice day!

[Serrof]

Hi,

If you want 2nd order or above, you’ll need to use FieldPropagator and the like. It gets a bit trickier to have custom parameters for the differentiation then but it’s still doable. See for instance this

Cheers,
Romain.