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

*To*: orekit-developers@orekit.org*Subject*: Re: [Orekit Developers] status on second order derivatives*From*: MAISONOBE Luc <luc.maisonobe@c-s.fr>*Date*: Wed, 19 Nov 2014 14:55:58 +0100*In-reply-to*: <546C9D28.4060309@nrl.navy.mil>*References*: <20141029110217.580826apojmibi95@messagerie.si.c-s.fr> <84b16bf7980e9abdb95a6b5c517df16f@buffalo.edu> <20141029115718.31806x7bx5h2dfym@messagerie.si.c-s.fr> <5451296E.4010106@nrl.navy.mil> <20141029201403.840318u7665khugr@messagerie.si.c-s.fr> <54514F4C.3090906@nrl.navy.mil> <54521CAE.5020401@c-s.fr> <5453D26A.4000604@nrl.navy.mil> <20141101092339.12633h7r1d1bw6mj@messagerie.si.c-s.fr> <545CE6B9.6060603@nrl.navy.mil> <20141110175103.60916i4uq7a3xgjb@messagerie.si.c-s.fr> <5461358B.5020509@nrl.navy.mil> <20141111160332.127758prq3tsats4@messagerie.si.c-s.fr> <546B5275.8070203@c-s.fr> <546C9D28.4060309@nrl.navy.mil>*User-agent*: Internet Messaging Program (IMP) H3 (4.3.9)

Hi Evan, Evan Ward <evan.ward@nrl.navy.mil> a écrit :

On 11/18/2014 09:06 AM, Luc Maisonobe wrote:Hi all, Here are some final updates about the derivatives and the velocity problem in Eckstein-Hechler.Evan Ward <evan.ward@nrl.navy.mil> a écrit : To check that the initial orbit matching is causing the errors, I found that if the initial orbit is shifted by a quarter period (or -0.75 period) that the EH ephemeris stays within 70 m of the numerical propagator.For better fitting, we can always use the PropagatorConverter implementations : they adjust the initial parameters on a complete sample instead of only one point.I have further studied this option and reached a state I consider now acceptable. The Eckstein-Hechler-Cartesian branch includes the latest fixes, as well as some attempt to document the behaviour.So IMHO we can either fit the velocity exactly at the epoch or fit the orbit parameters and tolerate a small velocity discrepancy at the epoch in return for an orbit closer to the numerical propagator. Both can work with option 2.I definitely think option 2 is much simpler, so if you can find a way to meke it work, I would be very happy.The solution implemented in this branch is the one proposed by Evan: directly build a CartesianOrbit, and forget about any new specialized orbit type. Looking at several ways to fit the initial orbits, it appeared that was has been done from the beginning (i.e. fit the circular parameters) works pretty well, so I ave kept it. Furthermore, this means the new Eckstein-Hechler propagator just propagates the exact same positions has it has always done, it only output them as Cartesian. So the circular orbit that are used in the equations of the model and that are also used for the initialization are now only internal ones, they are not output anymore. I compared this initialization with the PropagatorConverter Pascal created a few years ago, and which relies on a full sample and not a single point, and indeed the result are pretty close. So our initialization was good, and could be kept. So if users rebuild circular parameters from the generated Cartesian parameters, they will see a difference with respect to the initial orbit they used for building the propagator, but I now consider these differences to be normal and related to the velocity discrepancy we observed. If we recompute just the position from the circular parameters, they all match to sub-micrometer level (i.e. the position computed from the initial circular parameters, the position extracted from the propagated orbit at t = 0 and the position computed from converting the propagated orbit at t = 0 to circular). I have created an issue <https://www.orekit.org/forge/issues/180> and solved it. I have also updated the wiki page on propagation <https://www.orekit.org/forge/projects/orekit/wiki/Propagation#Eckstein-Hechler-propagation>. I think we could now merge branch Eckstein-Hechler-Cartesian (which includes the branch position-velocity-acceleration) back to master. What do you think ?Looks good to me. Thanks for getting to the bottom of this issue. One small thing, is EHPropagator.toCartesian(...) supposed to be public?

No, I will fix this and then merge the branch to master. thanks for your help on this issue, Luc

Best Regards, EvanFor now, I have pushed the two branches I use for experiments. One is called EHOrbit and implements option 3, the other is called Eckstein-Hechler-Cartesian and implements option 2. In both cases, I have updated the tests as much as possible to match the accuracy achieved. For option 3, now all tests without exception do pass. Interpolation is very good but shiftedBy behaves strangely (I feel the error increases too fast for my taste). For option 2, four tests in EcksteinHechlerPropagatorTest fail, as explained in my previous post. If we could get them working, I think this branch would be a much better solution to merge in master.Do you have access the reference publication for the EH propagator? I think it may help to see how the propagator is initialized there.Unfortunately, no. The implementation was created from an appendix in the CNES 1995 book, not from the original paper. This book is not available anymore as far as I know. I will look again at it but don't think there is anything about initialization. The initialization process we used is a very simple fixed-point iteration as the model creates osculating parameters from mean parameters : so we iterate until the produced osculating match the initial orbit, currently in terms of circular parameters.I have looked in the book. There was nothing about the initialization of the model, only a few pages of raw equations without any explanation. If you have access to this book (Spaceflight Dynamics - part I, 1995, collaborative work from CNES, editor Jean-Pierre Carrou, ISBN 2.85428.376.7), the model is in pp 474-481. best regards, Lucbest regards, LucBest Regards, EvanSo up to now, I was not able to have both initial orbit, propagated orbit and velocity correct. In all my tests, I had only two correct among the three.3. Create a combined orbit that has both the EH elements and the position and velocity. This may be a bit confusing for the users since the EH orbital elements are only define w.r.t. the EH propagation model.In fact, as seen in the previous case, the orbital elements are also visible to the user as they are used to initialize the propagator from an osculating orbit (notwithstanding velocity). The hidden parameters are the mean parameters from wich EH computes the osculating ones.In other words, if a user writes a method to use a Propagator (the base interface, not one of its sub classes) then the user wouldn't be able to use the elements (a, e, i, ...) in the returned orbit since they wouldn't know how the elements are defined. On the other hand, if the user needs the EH parameters then I think this is the definition we should choose. (IIUC this is what Luc is describing.)I tried it for the last few days (before switching to your suggestion number 2 above). It almost works but is quite complicated and I did not yet manage to solve a last problem. The way I've done it was to create an internal EHOrbit class extending CircularOrbit and overriding a few methods. At construction time it preserve the value of the circular osculating parameters (so initial orbit is preserved when propagating to same date), and it does not modify the mean parameters computation loop (so propagated orbit is consistent with other numerical propagators), and it additionally computes Cartesian parameters knowing about the derivatives (so it solves the initial problem). It's a bunch of code, but is quite straightforward thanks to the DerivativeStructure class which greatly simplifies derivatives computation. This however is not enough, and a few methods should be overriden too: shitfetBy and interpolate. Failing to overriding thme creates another type of inconsistency: the velocity of shifted or interpolated instances computed directly by the base CircularOrbit class is wrong. Overriding shiftedBy is trivial and works well. Overriding interpolate is difficult as it could be fed from several types of orbit in the sample, so I tried to rebuild the derivatives if the sample did not include them. This is a nightmare and I utterly failed. The problem is interpolate is automatically used for example when one creates an ephemeris from an EcksteinHechler propagator. It is a feature shared among all propagators and widely used. So I am confused now. I could attempt to reset the specialized EHOrbit class that knows about non-Keplerian effect when converting back and forth between circular and Cartesian, but I have yet to find how interpolate should work. Perhaps we could simply say that the specific implementation of interpolate in EHOrbit works *only* when the sample contains only EHOrbit instances and throws an exception otherwise. It seems fair to me as we cannot rebuild missing derivatives, and it would still work correctly in the regular case where people use consistently EHOrbit for a sample. In fact, it could probably be considered a general contract of the interpolate method to use only one type of orbits at a time and not to convert anything under the hood. Any thought ? LucThe problem may not hold for NumericalPropagator (I still have to check), because basically we do the computation the other way round. We start from x, y, z, vx, vy, vz and deduce a, e, ... from the mapping. So as long as our vx=dx/dt, vy=dy/dt, vz=dz/dt, the initial mapping will still hold.It seems that if we add rates to all the elements then the resulting class could be classified as a Propagator.It is only a way to get a consistent (P, V), and as a propagator it would be really limited as it is a Taylor expansion. I would simply qualify it as a perturbed orbit allowing local expansion.Extending the Taylor series expansion to higher order derivatives makes sense if we need to return the propagator's internal state, and the satellite's PV in one consistent object. Is there a use case for providing the propagator's internal state? I don't know, but I can't think of any. If we need to provide the user with acceleration as well as PV then I agree that we need a bigger container class to hold the extra information. Maybe the new PVACoordinates would be the right choiceI agree this is a good place. However, I did not create an extra class for this either, I simply added a third vector to PVCoordinates. We can discuss this later, when the merge of the branch will be attempted.instead of a PerturbedOrbit class since I don't think we would get the same extrapolation from a PerturbedOrbit stored in Keplerian vs Cartesian elements. To put it another way, If I converted a PerturbedKeplerianOrbit to a PerturbedCartesianOrbit would they both follow the same path?At second order level, in the neighborhood of the matching date, yes they should. Differences should appear at the first ignored derivative level and build up from there.I'm not against the PerturbedOrbit classes; I just want to make sure we don't choose a complex solution to a simple problem.You are right. I have a clear tendency to over-engineer things, this is why I ask for advice from time to time. This is the great force of collaborative development : people don't get stuck in their own errors too long.I think these constructive discussions combine the best ideas from everyone.So to be clear, here is the consensus : - I will simplify what I have done to limit as much as possible the changes to CircularOrbit. - In this setting, orbits will not know about non-Keplerian effects. - As this knowledge is needed for accurate and more importantly consistent velocity computation, the mapping between orbital parameters and Cartesian parameters will be precomputed by the propagators - If the mapping is not precomputed and orbits are built from one set of parameters only, then the mapping will be Keplerian only and will be performed byt orbit classes just as it is nowIf we decide to define the orbit returned from a Propagator according to option 3 above then I think this is a good plan to implement it. Best Regards, Evanbest regards, LucThere might be a case for this intermediate level of fidelity/speed, but is it worth the added complexity? Especially since we already have one that is fast (the Orbit class) and one that is precise (the Propagator classes). Maybe I'm still not understanding your proposal.I detected the problem when checking some pointing attitude modes that are related to spacecraft velocity (alignment of the spacecraft axis with ground drift for Earth observation). This mode transforms something that is a derivative (velocity) into a regular non-differentiated value (an angle). So we go up one order of derivation. For such modes, a wrong velocity induces a wrong angle and a wrong acceleration induces a wrong angular rate. As the angular rate is used for time shift, I needed it to be accurate. You also mentioned issues with Doppler which would also occur. We do need accurate velocity, and we do need a velocity that is consistent with the derivative of the position.Definitely agree that getting the velocity correct is important. Thanks again for all the work you've put into solving this bug. Best Regards, Evanbest regards, LucBest Regards, EvanThen we wouldn't have to modify the existing set of Orbit classes, and the user would see the correct osculating P/V. (This might be equivalent to your second approach.) As far as where to put the code, it seems like the conversion code would be specific to the internal state representation used by the propagator, so I think it makes sense to keep the code as private to the propagator. Though if you think there would be other uses for the conversion, I think a public static factory method would work well.Yes, the conversion code is propagator dependent. The first implementation I played with is Eckstein-Hechler propagator, and it definitely is Eckstein-Hechler specific (it's a simple derivation of the original equations, so each time we did compute something like e = a * c + b * d, now we also have another statement to compute eDot = aDot * c + a * dCot + bDot * d + b * dDot). For numerical propagator we already have the derivatives since we start from the derivatives and integrate them, so its even simpler. For DSST, this will be a mix as the mean elements are integrated and the short periodics terms are computed from Fourier coefficients which are straightforward to differentiate. For ephemeris-based propagator, we will need to compute the derivatives of the underlying polynomials, which is also straightforward.Thanks Luc for finding this issue and doing the analysis. I can see how this would be an issue when computing the Doppler as well as time shifting.your welcome best regards, LucBest Regards, Evan On 10/29/2014 06:57 AM, MAISONOBE Luc wrote:Hi Paul, paulcefo <paulcefo@buffalo.edu> a écrit :Luc, Do I correctly understand that your concern is that Keplerian transformations do apply outside the osculating space?The problem I had was that we did use Keplerian-only expression to set up local Taylor expansions around the current point (a few seconds away). This was slightly wrong when all the parameters were time-dependent and not only the anomaly was time-dependent. Of course, the error increasing with the time offset with respect to the central date at which the Taylor expansion is built. The fix was simply to not forget the derivatives of these other parameters. This Taylor expansion feature is a built-in feature available in all Orekit orbits, it typically allows to do computation in the vicinity of an already computed point without needed to trigger a complete propagator. It can even be used for some computation inside the run of a propagator, as for example when the higher level propagator takes care of the long term propagation and at each step we need some additional points surrounding the current step to compute attitude evolution in some specific modes. My concern was how to implement this fix in our current architecture, and more precisely were to put the code: in an existing class or in a dedicated class which would be used only by propagators. best regards, LucPaul -- Dr. Paul J. Cefola Consultant in Aerospace Systems, Spaceflight Mechanics, & Astrodynamics Adjunct Faculty, Dept. of Mechanical and Aerospace Engineering, University at Buffalo (SUNY) 4 Moonstone Way Vineyard Haven, MA 02568 USA 508-696-1884 (phone on Martha's Vineyard) 978-201-1393 (cell) paulcefo@buffalo.edu paul.cefola@gmail.com On 10/29/2014 6:02 am, MAISONOBE Luc wrote:Hello, As some of you may be aware, I have been working for a few months on second order derivatives in the git branch position-velocity-acceleration. This work is still ongoing but I hope to finish it for 7.0 and merge the branch back to master soon. For now, there are still failing tests so I can't do it. This change should allow us to reach several goals : - improved accuracy in shiftedBy methods - improved accuracy in interpolators (with user-defined choices to use or not first and second derivatives from the sample) - improved accuracy in attitude - removal of ugly hidden finite differences in some classes (most notably attitude modes) with hard-coded steps - hopefully faster Earth transforms, by replacing Hermite interpolation with single point extrapolation - availability of non-Keplerian acceleration everywhere - availability of angular acceleration in attitude and frames - proper composition of dynamics in frames - possibility to propagate orbits in non-inertial frames - possibility to propagate orbits without a central body (interplanetary missions, Lagrange point missions, ...) There is one point that bothers me right now. As I removed some of the ugly finite differences, some non-regression tests started to fail. I finally found the raw cause of these failures and was surprised to discover an old bug in the way we use the osculating orbits produced by the Eckstein-Hechler analytical propagator. This propagator takes zonal terms into account, and produces directly circular parameters a, ex, ey, ... When we compute anything related to geometry, we compute Cartesian coordinates using the Orbit getPVCoordinates method. As the Orbit classes do not know anything about the perturbation, the (P, V) pair does in fact implicitly relies on Keplerian-only expressions. So the velocity part is *not* consistent with the derivative of the position. The real derivative of the position takes the non-Keplerian effects into account which are ignored by getPVCoordinates. The difference is small, but as the tests threshold were deliberatly very tight, the tests started to fail when the various pointing directions were not computed anymore from finite differences mainly involving position and when they relied on the computed velocity. So the problem already happens in the master branch, it is not specific to the introduction of acceleration (it was just detected here during testing). The solution is in fact quite simple. If an orbit has been produced by a non-Keplerian propagator, the propagator already knows about the derivatives of the orbital elements (which are circular in the Eckstein-Hechler model case but can be any kind of parameters for other propagators). The propagator should therefore provide these derivatives to the orbit so they can be used in the PVCoordinates conversion. The code is very simple and straightforward. I have checked this and got very interesting results with Eckstein-Hechler/Circular, as for example a simple interpolation over a 900s arc with proper velocity/acceleration has a 88m error with two base points now whereas it was 5162 m before (and 0.02m vs 650m for 3 points, 1.0e-5m vs 259m for 4 points). Here is what bothers me: Should we create specialized classes for perturbed orbits or should we simply add a constructor to the existing orbits with the parameters derivatives and set them to 0 when they are not known? For my tests, I created PerturbedCircularOrbit which extends CircularOrbit and override the protected initPVCoordinates method and the public shiftedBy and interpolate methods. I could also have simply moved everything into CircularOrbit with a new constructor. I do not like much the PerturbedXxxxOrbit approach, as it forces to create also additional entries in the OrbitType enum with additional converters and it becomes awkward if for example a user configures a NumericalPropagator to generate XxxxOrbit, despite this propagator will in fact really generate PerturbedXxxOrbit because it is what a Numerical propagator is for. So there should be either an internal modification of the user setting from OrbitType.XXXX to OrbitType.PERTURBED_XXXX or an error triggered which would invalidate *all* current user code as it would become forbiddent to generate XXXX orbits now. On the other hand, the drawback of modifying the existing classes to hold the non-Keplerian derivatives is that they will consume more memory. I don't think it is a problem with current computers. In any case, initial orbits created directly from user code or by reading files would not include the derivatives and therefore will be built as usual (by calling the unmodified classes in the first approach, or by using the already existing constructors in the second approach, assuming these constructors will automatically set the derivatives to Keplerian-only values). In any case, full-blown perturbed orbits will be created internally by Orekit propagators, which can easily be modified to provide the derivatives they know (by creating instances of the new derived classes in the first approach, or by using new constructors with additional parameters in the second approach). My humble opinion would be to use the second approach to solve this bug. I will probably do this in the position-velocity-acceleration branch so it will include accelerations right from the start and will be merged to master at the same time as the rest of the branch. Of course, this will be a dedicated commits (Git branches are great!). What do you think ? best regards, Luc ---------------------------------------------------------------- This message was sent using IMP, the Internet Messaging Program.---------------------------------------------------------------- This message was sent using IMP, the Internet Messaging Program.---------------------------------------------------------------- This message was sent using IMP, the Internet Messaging Program.---------------------------------------------------------------- This message was sent using IMP, the Internet Messaging Program.

**References**:**[Orekit Developers] status on second order derivatives***From:*MAISONOBE Luc <luc.maisonobe@c-s.fr>

**Re: [Orekit Developers] status on second order derivatives***From:*paulcefo <paulcefo@buffalo.edu>

**Re: [Orekit Developers] status on second order derivatives***From:*MAISONOBE Luc <luc.maisonobe@c-s.fr>

**Re: [Orekit Developers] status on second order derivatives***From:*Evan Ward <evan.ward@nrl.navy.mil>

**Re: [Orekit Developers] status on second order derivatives***From:*MAISONOBE Luc <luc.maisonobe@c-s.fr>

**Re: [Orekit Developers] status on second order derivatives***From:*Evan Ward <evan.ward@nrl.navy.mil>

**Re: [Orekit Developers] status on second order derivatives***From:*Luc Maisonobe <Luc.Maisonobe@c-s.fr>

**Re: [Orekit Developers] status on second order derivatives***From:*Evan Ward <evan.ward@nrl.navy.mil>

**Re: [Orekit Developers] status on second order derivatives***From:*MAISONOBE Luc <luc.maisonobe@c-s.fr>

**Re: [Orekit Developers] status on second order derivatives***From:*Evan Ward <evan.ward@nrl.navy.mil>

**Re: [Orekit Developers] status on second order derivatives***From:*MAISONOBE Luc <luc.maisonobe@c-s.fr>

**Re: [Orekit Developers] status on second order derivatives***From:*Evan Ward <evan.ward@nrl.navy.mil>

**Re: [Orekit Developers] status on second order derivatives***From:*MAISONOBE Luc <luc.maisonobe@c-s.fr>

**Re: [Orekit Developers] status on second order derivatives***From:*Luc Maisonobe <Luc.Maisonobe@c-s.fr>

**Re: [Orekit Developers] status on second order derivatives***From:*Evan Ward <evan.ward@nrl.navy.mil>

- Prev by Date:
**Re: [Orekit Developers] html entities versus UTF-8** - Next by Date:
**[Orekit Developers] Preparation for next release** - Previous by thread:
**Re: [Orekit Developers] status on second order derivatives** - Next by thread:
**[Orekit Developers] html entities versus UTF-8** - Index(es):