Examples of ForwardDiff or JSMDUtils.Autodiff

[travelingspaceman]

Thanks to the devs for writing a very useful package! One of the most attractive features to me is the compatibility with ForwardDiff, but I can't figure out how to get that to work. I don't see any examples in the docs or tests. It looks like automatic differentiation is a work in progress. A couple useful use cases that I would like to be able to use are:

• Jacobian of spacecraft state in frame A, with respect to spacecraft state in frame B
• Gradient of spacecraft state in frame A, with respect to epoch

The issue I have run into is that update_point! requires the number type of the new state to match the number type of the FrameSystem. However, number types with autodiff get complicated with tagging, and I'm not sure how to construct the FrameSystem to get around this.

[andreapasquale94]

Hi @travelingspaceman! Thank you for giving a shot to FramesTransformations.jl and for the feedback.

We are still working on the documentation on the package as well as improving the AD compatibility, at some point we'll also have those cases fully documented and tested :smile:

Let me give you a bit of context here.

We have mostly treated AD compatibility as "compatibility in time", e.g. being able to differentiate the FrameSystem in the independent variable (i.e., time). This time-compatible AD has been mostly handled from two points of view: the possibility to differentiate within ephemerides (achieved back-ending the frame system with Ephemerides) and the possibility to create frames with those points (the computable axes).

Now that the package is consolidated we are starting to work on the design capabilities, e.g. the one associated with the spacecraft trajectory. Also in this case there are two cases, the most complex being the one associated with Trajectory Optimisation\Navigation Analysis.

In this case, the actual trajectory is under design and the trajectory might be not completely available - this is a complex case as the AD might be associated with a partially propagated trajectory, and we might want to have derivatives for different quantities: time, states, and parameters.

Time

These are already handled by the FrameSystem but a proper orbit representation shall be provided, as you said, to avoid perturbation confusion as well as to properly handle the AD. As you might have noticed, 2 point types could be suited for this: updatable and dynamical points. For this case updatable points are not suited as they are practically constants for the AD system for now, but you can properly handle the case with dynamical points. Here an example:

using Ephemerides
using FrameTransformations
using JSMDInterfaces.Math: interpolate
using JSMDUtils.Math: InterpCubicSplines

# Load ephemeris to memory
eph = load(EphemerisProvider, ["/home/andrew/Documents/Kernels/spk/de440s.bsp"])

# Create the graph
G = FrameSystem{4, Float64}(eph)

# Register the ICRF axes 
add_axes_inertial!(G, :ICRF, 1)

# Register the root node 
add_point_root!(G, :SSB, 0, 1)

# Register the other nodes 
add_point_ephemeris!(G, :EarthB, 3)
add_point_ephemeris!(G, :Sun, 10)
add_point_ephemeris!(G, :Earth, 399)

# Create dummy orbits
c = 2π/86400
de = 0:3600:10*86400;
fv(t) = [cos(c*t), sin(c*t), 0, -c*sin(c*t), c*cos(c*t), 0];
fv2(t) = [sin(c*t), cos(c*t), 0];
const ORB1 = InterpCubicSplines(de, hcat([fv(e) for e in de]...));
const ORB2 = InterpCubicSplines(de, hcat([fv2(e) for e in de]...));

# Insert dynamical points
add_point_dynamical!(G, :SC1, -1, 399, 1, t->interpolate(ORB1, t)[1:3], t->interpolate(ORB1, t))
add_point_dynamical!(G, :SC2, -2, 399, 1, t->interpolate(ORB2, t))

# Gradients of SC1 w.r.t. center (Earth)
vector12(G, 399, -1, 1, 12345.0)[7:12]

# Gradient of SC1 w.r.t. SC2 
vector12(G, -1, -2, 1, 12345.0)[7:12]

Which is answering:

Gradient of spacecraft state in frame A, with respect to epoch

State/parameters

Jacobian of spacecraft state in frame A, with respect to spacecraft state in frame B

These instead are not directly associated with the FrameSystem but to the propagator as requires State Transition Matrice (STM) and sensitivities. You can still use the previous approach to register a propagated point within the frame system but then the AD might be not straightforward to set up.

If you are interested, we are planning to develop those features within the ecosystem this year, starting from the automatic insertion of state transformations in the chain with AstroRepresentations.jl, the different astrodynamical force models (we have wip private package for those) and then interfacing with SciML DifferentialEquations ecosystem.

[MicheleCeresoli]

The issue I have run into is that update_point! requires the number type of the new state to match the number type of the FrameSystem. However, number types with autodiff get complicated with tagging, and I'm not sure how to construct the FrameSystem to get around this.

As @andreapasquale94 has said, the complete AD compatibility of the different transformations with respect to time is still a work-in-progress and should hopefully be finished in the next couple of months. In particular, some issues have indeed arisen with the UpdatablePoints that you mention: in all other cases, the state only depends on time and the time derivative is straightforward. However, some ambiguities can arise when considering a point that is arbitrarily updated by the user (which might then not depend on time).

The current idea we had in mind is to force the user to update the state with a vector of dual numbers when trying to use Autodiff with updatable states, although I'm not totally sure whether this is always feasible.

This means that, although the documentation may suggest differently, you cannot yet apply ForwardDiff.derivative to a vector3 call. However, if you are looking to retrieve the acceleration\jerk you can always get those by calling the vector9 or vector12 functions (which internally leverage ForwardDiff to automatically compute the higher order derivatives if they have not been specified).

[andreapasquale94]

@travelingspaceman does this answer your question?

A dedicated tutorial will be inserted in the docs based on this discussion.