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

*To*: <orekit-developers@orekit.org>*Subject*: Re: [Orekit Developers] Taylor differential algebra*From*: Evan Ward <evan.ward@nrl.navy.mil>*Date*: Thu, 3 Nov 2016 10:27:57 -0400*In-reply-to*: <20161103121915.Horde.RuvjrSWmMWyGPAAHIu267Q7@messagerie.si.c-s.fr>*References*: <20161103121915.Horde.RuvjrSWmMWyGPAAHIu267Q7@messagerie.si.c-s.fr>*User-agent*: Mozilla/5.0 (X11; Linux x86_64; rv:45.0) Gecko/20100101 Thunderbird/45.4.0

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

Hi all,As mentioned last week, a student worked for me during the last fewmonthson 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. Itwillbe committed to the main git repository soon (don't ask me when ...).Andrea didpresent it on November 1st during the Stardust conference at ESTEC, inthe Orbitand Uncertainty Propagation track.The approach we followed was to create generic Field-basedpropagators, so userscan do things like:FieldPropagator<SomeField> propagator = newFieldNumericalPropagator<SomeField>(...);(or some other supported analytical propagator if desired) and then touse it just likea regular double-based propagator.In order to use Taylor differential algebra, users just have to useDerivativeStructureas the field for propagated elements. It is up to them to decide uponthe differentiationorder and the number of parameters. The most classical use for this isto select anorder of about 3 and a number of parameters of 6, each parametercorresponding toone of the component of initial states. Of course, any type of orbitrepresentation canbe used (Cartesian, Keplerian, Circular, Equinoctial). As the fieldelements (hereDerivativeStructures instances) are propagated, you have at all time,both duringpropagation in step handlers and in events handlers and at the filnaltime the full setof partial derivatives up to the configured order. This mean youpropagate uncertaintiesto a higher order than classical state transition matrices (orcovariance matrices) whichare inherently first order only. So if for example you select toretrive the output stateas Cartesian coordinates, you still see the uncertainty cloud bend tofollow the orbitshape, you do not get only an ellipsoid mainly tangent to the orbit.Another very importantuse is Monte-Carlo simulations. You do the same as before, propagateonce, and then on theresult DerivativeStructure instances, you can call the taylor(dp1,dp2, ...) method toevaluate the Taylor expansion to small offsets in the parameters. Theoffsets dp1, dp2 ...are drawn randomly to correspond to your initial distribution and youget the correspondingoutput without redoing the propagation: you are just evaluating lowdegree polynomials.This speeds up Monte-Carlo simulations a lot and you can doMonte-Carlo analysis withhundred of thousands or millions of points if you want. Of course, theinitial propagationis much slower, especially if you have high order, but the overhead isdone only once, sothe complete cost is easily paid for for any realistic Monte-Carlosize. As an example,the overhead multiplicative factor with respect to a traditionaldouble-based propagationand derivatives computed with respect to 6 parameters is 68.13 fordegree 2 derivativesand 116.2 for degree 3, using a numerical propagator and all availableforce models. One hasto remember that for a DerivativeStructure with p parameters and ordero, we compute(o+p)!/(o! p!) components instead of 1 component in double-basedcomputation. Taking thiscombinatorial factor to normalize the overhead, we find that forparameters and ordersbetween 1 and 6, the normalized overhead is about 2 to 5, which seemsquite reasonable.Beware! Do not try to use too large o or p, as the combinatorialfactors grows up quickly!For example using 6x6, each double number is replaced by a 924components DerivativeStructure,the un-normalized overhead is more than 3700, and the normalizedoverhead is between 4 and 5.The field output is mainly a FieldSpacecraftState, but of course usercan also use derivedparameters (position/velocity/acceleration, in any frame, inertial orrotating, geodeticcoordinates ...). So for example you could propagate an uncertainty ingeodetic coordinateswith a rotating Earth, for example to estimate a ground impact zone,supposing you areable to reliably propagate down to ground, which of course is notreally easy ...Andrea already implemented Keplerian, Eckstein-Hechler and numericalpropagators with a bunchof force models. We may add SGP4/SDP4 also, it is quitestraightforward to do.A side-effect of our approach is that we can use it to implement otherthings apart fromTaylor differential algebra (high accuracy, interval arithmetics,...). This side effectwas in fact used for validating the implementation using Decimal64 asthe field to get backto 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 hadto be almost duplicatedwith a FieldXYZ implementation that mimicked the regular double-basedimplementation. Thoseof you who followed Apache Commons Math and later Hipparchgusdevelopment know what that means,as it is exactly what has been done for the ODE integrators. By theway, this Taylor differentialalgebra idea was a main driver even at that time, the classes inHipparchus were developed sothis work in Orekit can be done later on. This approach leads to ahuge amount of code and amaintenance hell ahead of us. Perhaps we should revive the approach ofautomatic conversion inthe spirit of what I created years ago with the Nabla project atApache. I don't know if weshould put Andrea code directly in the main library or split Orekit inmultiple modules, onebeing Field classes (including the existing ones added in the last fewyears for other purposes).It is however not as easy as it first looks, because in some casesinstead of new independantclasses, we just added Field methods in existing classes andinterfaces. An already existingexample in 8.0 is for example the accelerationDerivatives method inForceModel interface, whichalready needs FieldVector3D and FieldRotation (but specifically withDerivativeStructure asthe field, which by the way was a stupid design decision I made).So what do you think? Is this feature interesting? Should it be mergedin the main library?Should we spawn some modules? Can we do that for 9.0?

Andrea and Luc,

Best Regards, Evan

best regards, Luc

**References**:**[Orekit Developers] Taylor differential algebra***From:*MAISONOBE Luc <luc.maisonobe@c-s.fr>

- Prev by Date:
**Re: [Orekit Developers] Ephemeris Writer Interface** - Next by Date:
**Re: [Orekit Developers] Ephemeris Writer Interface** - Previous by thread:
**[Orekit Developers] Taylor differential algebra** - Next by thread:
**[Orekit Developers] Need a STOA detector** - Index(es):