Re: [Orekit Users] Generic extremum events

[MAISONOBE Luc <luc.maisonobe@c-s.fr>]

Thomas.Fruth@dlr.de a écrit :

> Hello,

Hi Thomas,

> at the DLR GSOC we are presently evaluating whether Orekit could be  
> integrated into our Mission Planning tool suite.
>
> One important functionality that we need is the calculation of  
> extremum events, such as for example the time of maximum elevation  
> in a given topocentric frame or the time of minimum range with  
> respect to a given Earth observation target.
>
> From investigating Orekit we identified two different approaches to  
> achieve this task:
>
>
> 1)      To return the derivative of the value that is to be  
> minimized or maximized in the function g in the specific  
> implementation of AbstractDetector. This approach is for instance  
> chosen by the new Latitude/LongitudeExtremumDetector classes;  
> however, it requires some knowledge about the derivative (as in this  
> example handled by OneAxisEllipsoid.transform(PVCoordinates, Frame,  
> AbsoluteDate)), which might not always be readily available.
>
> 2)      To extend the present event handling mechanism such that it  
> allows to configure whether the given function g should be searched  
> for its root, or be minimized or maximized. My naïve guess is that  
> an extremum search without passing an explicit derivative might be  
> achieved relatively easily by switching between the solver presently  
> used in EventState.evaluateStep(...) and an optimizer (e.g. the  
> BrentOptimizer from Apache commons-math) as needed.
>
> I would be glad if some experienced Orekit users could comment  
> whether my assessment on using Orekit for extremum events is correct,

Yes, your assumptions are correct.

> or if I might have missed some other possibility to achieve this  
> task with Orekit's present functionality. If there yet doesn't exist  
> a generic solution to find extremum events (such as described in  
> option 2 above), are there any plans to include such functionality?

In fact, I think the first approach is the more straightforward one.
When the g function of your ExtremumElevationDetector will be called,
it will get a SpacecraftState which does contain the information about
position, velocity and acceleration which are sufficient to compute
simply all derivatives. You don't even need to care by yourself about
the exact expression for the derivatives since DerivativeStructure from
Apache Commons Math can do it for you. Here is how it can be done:

  public double g(final SpacecraftState s) throws OrekitException {

    // get position, velocity acceleration of spacecraft in topocentric frame
    final Transform inertToTopo = s.getFrame().getTransformTo(topo,  
s.getDate());
    final TimeStampedPVCoordinates pvTopo =  
inertToTopo.transformPVCoordinates(s.getPVCoordinates());

    // convert the coordinates to DerivativeStructure based vector
    // instead of having vector position, then vector velocity then  
vector acceleration
    // we get one vector and each coordinate is a DerivativeStructure  
containing
    // value, first time derivative, second time derivative
    final FieldVector3D<DerivativeStructure> pvDS =  
pvTopo.toDerivativeStructureVector(2);

    // compute elevation and its first and second derivatives
    final DerivativeStructure elevation =
                    pvDS.getZ().divide(pvDS.getNorm()).asin();

    // return elevation first time derivative
    return elevation.getPartialDerivative(1);

  }

Note that in the method above, we never write explicitely any  
derivative, they are
computed analytically using chain rule thanks to Apache Commons Math.


The second approach as you describe it would lead to some problems as when the
propagator used is a numerical propagator, it is not the Orekit EventState
class that is used but the one from Apache Commons Math in the ode  
package. This
means that the same change would have to be done in both libraries. Its not
impossible, but cumbersome. There may be a better approach in similar  
cases when
direct differentiation is not possible, it is to use a wrapper  
function that does
perform the differentiation for you and use the regular event detection on the
wrapped function.

In other words, you could implement an f function that computes the  
function you
want and then wrap it using a finite differences differentiator that  
would create
function g by wrapping function f so that each time g(state) is  
called, in fact
you would get two calls to f. Basically you would end up with somethig like:

  double g(state) {
    double fPlus  = wrapped.f(state.shiftedBy(+step));
    double fMinus = wrapped.f(state.shiftedBy(-step));
    return (fPlus - fMinus) / (2 * step);
  }

With such a wrapper, implementors can easily get extremums without  
worrying about
difficult derivatives, and no changes to any library is needed.
The shiftedBy method is well suited for small offsets like the ones  
used in finite
differences.

As the detector you mention is really simple, I think we may add it in  
the next
few days, using the code snippet above for the first approach. It is  
clearly in
line with the new LatitudeExtremum and LongitudeExtremum added recently.

Hope this helps,
Luc

> Best regards
>
> Thomas
>
> --------------------------
> Deutsches Zentrum für Luft- und Raumfahrt e.V. (DLR)
> German Aerospace Center
> Space Operations and Astronaut Training | Mission Operations |  
> Oberpfaffenhofen | 82234 Wessling | Germany
>
> Dr. Thomas Fruth
> Telephone +49 8153 28-2432 | Telefax +49 8153 28-1456 |  
> thomas.fruth@dlr.de<mailto:thomas.fruth@dlr.de>
> DLR.de<http://www.dlr.de/>