I see your points regarding a generic solution that is not using
derivatives. Obviously I have really underestimated the power of
Apache commons-math's DerivativeStructure: Using it and with the help
of your ElevationExtremumDetector example, I was now able to implement
another use case (minimum range) successfully and very elegantly in
only a few lines of code. Taken together the DerivativeStructure
option and the wrapper you described as a kind of backup solution, I
agree there shouldn't be much need for a more generic extremum event
search.
Thanks also for directly adding the ElevationExtremumDetector to the
repository! I have just checked it and its results agree very well
with our reference data.
Best regards,
Thomas
PS (just out of curiosity): The new ElevationExtremumDetector requires
an OneAxisEllipsoid and TopocentricFrame on initialization, whereas
ElevationDetector is fine with the TopocentricFrame only.
Since the OneAxisEllipsoid is not really used within
ElevationExtremumDetector either, I just wondered why it is needed?
-----Ursprüngliche Nachricht-----
Von: orekit-users-request@orekit.org
[mailto:orekit-users-request@orekit.org] Im Auftrag von Luc Maisonobe
Gesendet: Mittwoch, 9. September 2015 09:39
An: orekit-users@orekit.org
Betreff: Re: [Orekit Users] Generic extremum events
Hi Thomas,
Le 08/09/2015 20:43, MAISONOBE Luc a écrit :
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.
Another problem with the extremum solver is that we would have problem
with the pre-filtering of the search interval.
Currently, the root solver is not triggered at each step, bit only
when a root is known to lie in the step (according to the
maxCheckInterval setting). This two stages search saves a *lot* of
computation.
If we want to search for extremum, we would have to do the same and
identifying an extremum occurs would imply computing the derivative
anyway, so if we have it, we can just use it all the way through,
first to bracket the search interval, then to find the root.
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.
I have just added the ExtremumElevationDetector in the development
version in the git repository, you can give it a try now.
best regards,
Luc
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/>