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.