[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Orekit Developers] Taylor differential algebra





On 11/03/2016 07:19 AM, MAISONOBE Luc wrote:
Hi all,

As mentioned last week, a student worked for me during the last few months
on Orekit. His name is Andrea Antolino and his task wes to implement
Taylor differential algebra for orekit propagators.

This work is almost complete by now. We are fixing some lose ends. It will be committed to the main git repository soon (don't ask me when ...). Andrea did present it on November 1st during the Stardust conference at ESTEC, in the Orbit
and Uncertainty Propagation track.

The approach we followed was to create generic Field-based propagators, so users
can do things like:

FieldPropagator<SomeField> propagator = new FieldNumericalPropagator<SomeField>(...);

(or some other supported analytical propagator if desired) and then to use it just like
a regular double-based propagator.

In order to use Taylor differential algebra, users just have to use DerivativeStructure as the field for propagated elements. It is up to them to decide upon the differentiation order and the number of parameters. The most classical use for this is to select an order of about 3 and a number of parameters of 6, each parameter corresponding to one of the component of initial states. Of course, any type of orbit representation can be used (Cartesian, Keplerian, Circular, Equinoctial). As the field elements (here DerivativeStructures instances) are propagated, you have at all time, both during propagation in step handlers and in events handlers and at the filnal time the full set of partial derivatives up to the configured order. This mean you propagate uncertainties to a higher order than classical state transition matrices (or covariance matrices) which are inherently first order only. So if for example you select to retrive the output state as Cartesian coordinates, you still see the uncertainty cloud bend to follow the orbit shape, you do not get only an ellipsoid mainly tangent to the orbit. Another very important use is Monte-Carlo simulations. You do the same as before, propagate once, and then on the result DerivativeStructure instances, you can call the taylor(dp1, dp2, ...) method to evaluate the Taylor expansion to small offsets in the parameters. The offsets dp1, dp2 ... are drawn randomly to correspond to your initial distribution and you get the corresponding output without redoing the propagation: you are just evaluating low degree polynomials. This speeds up Monte-Carlo simulations a lot and you can do Monte-Carlo analysis with hundred of thousands or millions of points if you want. Of course, the initial propagation is much slower, especially if you have high order, but the overhead is done only once, so the complete cost is easily paid for for any realistic Monte-Carlo size. As an example, the overhead multiplicative factor with respect to a traditional double-based propagation and derivatives computed with respect to 6 parameters is 68.13 for degree 2 derivatives and 116.2 for degree 3, using a numerical propagator and all available force models. One has to remember that for a DerivativeStructure with p parameters and order o, we compute (o+p)!/(o! p!) components instead of 1 component in double-based computation. Taking this combinatorial factor to normalize the overhead, we find that for parameters and orders between 1 and 6, the normalized overhead is about 2 to 5, which seems quite reasonable. Beware! Do not try to use too large o or p, as the combinatorial factors grows up quickly! For example using 6x6, each double number is replaced by a 924 components DerivativeStructure, the un-normalized overhead is more than 3700, and the normalized overhead is between 4 and 5.

The field output is mainly a FieldSpacecraftState, but of course user can also use derived parameters (position/velocity/acceleration, in any frame, inertial or rotating, geodetic coordinates ...). So for example you could propagate an uncertainty in geodetic coordinates with a rotating Earth, for example to estimate a ground impact zone, supposing you are able to reliably propagate down to ground, which of course is not really easy ...

Andrea already implemented Keplerian, Eckstein-Hechler and numerical propagators with a bunch of force models. We may add SGP4/SDP4 also, it is quite straightforward to do.

A side-effect of our approach is that we can use it to implement other things apart from Taylor differential algebra (high accuracy, interval arithmetics, ...). This side effect was in fact used for validating the implementation using Decimal64 as the field to get back
to a (costly) way to redo exactly the same as regular double computation.

Now for the ugly part ... In order to do this work, *many* classes had to be almost duplicated with a FieldXYZ implementation that mimicked the regular double-based implementation. Those of you who followed Apache Commons Math and later Hipparchgus development know what that means, as it is exactly what has been done for the ODE integrators. By the way, this Taylor differential algebra idea was a main driver even at that time, the classes in Hipparchus were developed so this work in Orekit can be done later on. This approach leads to a huge amount of code and a maintenance hell ahead of us. Perhaps we should revive the approach of automatic conversion in the spirit of what I created years ago with the Nabla project at Apache. I don't know if we should put Andrea code directly in the main library or split Orekit in multiple modules, one being Field classes (including the existing ones added in the last few years for other purposes). It is however not as easy as it first looks, because in some cases instead of new independant classes, we just added Field methods in existing classes and interfaces. An already existing example in 8.0 is for example the accelerationDerivatives method in ForceModel interface, which already needs FieldVector3D and FieldRotation (but specifically with DerivativeStructure as
the field, which by the way was a stupid design decision I made).

So what do you think? Is this feature interesting? Should it be merged in the main library?
Should we spawn some modules? Can we do that for 9.0?

Andrea and Luc,

I think this is a good feature to have and makes these higher order derivatives more accessible. I wrote paper a little while back on MC sensitivity analysis and it seems like the method you implemented could be a good middle ground between linear models and full MC evaluation. I would be interested to read the paper you presented, if you can share it. It would also be interesting to see how the higher order derivatives could be used for OD, eg. preventing local minima, radius of convergence, speed of convergence.

On the implementation, several classes in Orekit currently use finite differences to compute derivatives, e.g. DragForce when computing the density gradient for an Atmosphere. It seems that the finite difference code would become complex to account for the variable differentiation order. Other sections such as HolmesFeatherstone use specialized code to compute derivatives that takes advantage of the structure of the equations to improve efficiency. Have these code sections been converted to use DerivativeStructure, or are they accounted for in some other way?

Are you planning to change the generic signature of ForceModel so that it can use any field and not just DerivativeStructure?

As for maintaining the code having several copies of the same class does present a headache. For several classes, e.g. EventState, this means there will be four slightly different copies of them, 2 in Hipparchus and 2 in Orekit. I'm not sure if there is a better way as the issue seems to be that Java generics can't be used for primitive values and that using primitive values provides a significant speed up. I think the new feature should be kept in the same git repository as Orekit so it is easy to see that updates are replicated to all the relevant classes.

Best Regards,
Evan


best regards,
Luc