AbstractPropagatorConverter.java

  1. /* Copyright 2002-2013 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.propagation.conversion;

  19. import java.util.ArrayList;
  20. import java.util.Collection;
  21. import java.util.List;

  23. import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
  24. import org.apache.commons.math3.analysis.MultivariateVectorFunction;
  25. import org.apache.commons.math3.exception.MaxCountExceededException;
  26. import org.apache.commons.math3.optim.InitialGuess;
  27. import org.apache.commons.math3.optim.MaxEval;
  28. import org.apache.commons.math3.optim.PointVectorValuePair;
  29. import org.apache.commons.math3.optim.SimpleVectorValueChecker;
  30. import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction;
  31. import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian;
  32. import org.apache.commons.math3.optim.nonlinear.vector.Target;
  33. import org.apache.commons.math3.optim.nonlinear.vector.Weight;
  34. import org.apache.commons.math3.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer;
  35. import org.apache.commons.math3.util.FastMath;
  36. import org.orekit.errors.OrekitException;
  37. import org.orekit.errors.OrekitExceptionWrapper;
  38. import org.orekit.errors.OrekitMessages;
  39. import org.orekit.frames.Frame;
  40. import org.orekit.propagation.Propagator;
  41. import org.orekit.propagation.SpacecraftState;
  42. import org.orekit.time.AbsoluteDate;
  43. import org.orekit.utils.PVCoordinates;

  45. /** Common handling of {@link PropagatorConverter} methods for propagators conversions.
  46.  * <p>
  47.  * This abstract class factorizes the common code for propagators conversion.
  48.  * Only one method must be implemented by derived classes: {@link #getObjectiveFunction()}.
  49.  * </p>
  50.  * <p>
  51.  * The converter uses the LevenbergMarquardtOptimizer from the <a
  52.  * href="http://commons.apache.org/math/">commons math</a> library.
  53.  * Different implementations correspond to different methods for computing the jacobian.
  54.  * </p>
  55.  * @author Pascal Parraud
  56.  * @since 6.0
  57.  */
  58. public abstract class AbstractPropagatorConverter implements PropagatorConverter {

  60.     /** Spacecraft states sample. */
  61.     private List<SpacecraftState> sample;

  63.     /** Reference date for conversion (1st date of the states sample). */
  64.     private AbsoluteDate date;

  66.     /** Target position and velocities at sample points. */
  67.     private double[] target;

  69.     /** Weight for residuals. */
  70.     private double[] weight;

  72.     /** Auxiliary outputData: RMS of solution. */
  73.     private double rms;

  75.     /** Position use indicator. */
  76.     private boolean onlyPosition;

  78.     /** Adapted propagator. */
  79.     private Propagator adapted;

  81.     /** List of the desired free parameters names. */
  82.     private Collection<String> parameters;

  84.     /** List of the available free parameters names. */
  85.     private final Collection<String> availableParameters;

  87.     /** Propagator builder. */
  88.     private final PropagatorBuilder builder;

  90.     /** Frame. */
  91.     private final Frame frame;

  93.     /** Optimizer for fitting. */
  94.     private final LevenbergMarquardtOptimizer optimizer;

  96.     /** Maximum number of iterations for optimization. */
  97.     private final int maxIterations;

  99.     /** Build a new instance.
  100.      * @param builder propagator builder
  101.      * @param threshold absolute convergence threshold for optimization algorithm
  102.      * @param maxIterations maximum number of iterations for fitting
  103.      */
  104.     protected AbstractPropagatorConverter(final PropagatorBuilder builder,
  105.                                           final double threshold,
  106.                                           final int maxIterations) {
  107.         this.builder             = builder;
  108.         this.frame               = builder.getFrame();
  109.         this.availableParameters = builder.getParametersNames();
  110.         this.optimizer           = new LevenbergMarquardtOptimizer(new SimpleVectorValueChecker(-1.0, threshold));
  111.         this.maxIterations       = maxIterations;
  112.         this.sample              = new ArrayList<SpacecraftState>();
  113.     }

  115.     /** Convert a propagator to another.
  116.      * @param source initial propagator
  117.      * @param timeSpan time span for fitting
  118.      * @param nbPoints number of fitting points over time span
  119.      * @param freeParameters names of the free parameters
  120.      * @return adapted propagator
  121.      * @exception OrekitException if propagator cannot be adapted
  122.      * @exception IllegalArgumentException if one of the parameters cannot be free
  123.      */
  124.     public Propagator convert(final Propagator source,
  125.                               final double timeSpan,
  126.                               final int nbPoints,
  127.                               final Collection<String> freeParameters) throws OrekitException {

  129.         checkParameters(freeParameters);
  130.         final List<SpacecraftState> states = createSample(source, timeSpan, nbPoints);
  131.         return convert(states, false, freeParameters);
  132.     }

  134.     /** Convert a propagator to another.
  135.      * @param source initial propagator
  136.      * @param timeSpan time span for fitting
  137.      * @param nbPoints number of fitting points over time span
  138.      * @param freeParameters names of the free parameters
  139.      * @return adapted propagator
  140.      * @exception OrekitException if propagator cannot be adapted
  141.      * @exception IllegalArgumentException if one of the parameters cannot be free
  142.      */
  143.     public Propagator convert(final Propagator source,
  144.                               final double timeSpan,
  145.                               final int nbPoints,
  146.                               final String ... freeParameters) throws OrekitException {

  148.         checkParameters(freeParameters);
  149.         final List<SpacecraftState> states = createSample(source, timeSpan, nbPoints);
  150.         return convert(states, false, freeParameters);
  151.     }

  153.     /** Find the propagator that minimize the mean square error for a sample of {@link SpacecraftState states}.
  154.      * @param states spacecraft states sample to fit
  155.      * @param positionOnly if true, consider only position data otherwise both position and velocity are used
  156.      * @param freeParameters names of the free parameters
  157.      * @return adapted propagator
  158.      * @exception OrekitException if propagator cannot be adapted
  159.      * @exception IllegalArgumentException if one of the parameters cannot be free
  160.      */
  161.     public Propagator convert(final List<SpacecraftState> states,
  162.                               final boolean positionOnly,
  163.                               final Collection<String> freeParameters) throws OrekitException {

  165.         checkParameters(freeParameters);

  167.         parameters = new ArrayList<String>();
  168.         parameters.addAll(freeParameters);

  170.         builder.setFreeParameters(parameters);

  172.         return adapt(states, positionOnly);
  173.     }

  175.     /** Find the propagator that minimize the mean square error for a sample of {@link SpacecraftState states}.
  176.      * @param states spacecraft states sample to fit
  177.      * @param positionOnly if true, consider only position data otherwise both position and velocity are used
  178.      * @param freeParameters names of the free parameters
  179.      * @return adapted propagator
  180.      * @exception OrekitException if propagator cannot be adapted
  181.      * @exception IllegalArgumentException if one of the parameters cannot be free
  182.      */
  183.     public Propagator convert(final List<SpacecraftState> states,
  184.                               final boolean positionOnly,
  185.                               final String ... freeParameters) throws OrekitException {

  187.         checkParameters(freeParameters);

  189.         parameters = new ArrayList<String>();
  190.         for (final String name : freeParameters) {
  191.             parameters.add(name);
  192.         }

  194.         builder.setFreeParameters(parameters);

  196.         return adapt(states, positionOnly);
  197.     }

  199.     /** Get the available free parameters.
  200.      * @return available free parameters
  201.      */
  202.     public Collection<String> getAvailableParameters() {
  203.         return availableParameters;
  204.     }

  206.     /** Check if a parameter can be free.
  207.      * @param name parameter name to check
  208.      * @return true if the parameter can be free
  209.      */
  210.     public boolean isAvailable(final String name) {
  211.         for (final String supportedName : availableParameters) {
  212.             if (supportedName.equals(name)) {
  213.                 return true;
  214.             }
  215.         }
  216.         return false;
  217.     }

  219.     /** Get the adapted propagator.
  220.      * @return adapted propagator
  221.      */
  222.     public Propagator getAdaptedPropagator() {
  223.         return adapted;
  224.     }

  226.     /** Get the Root Mean Square Deviation of the fitting.
  227.      * @return RMSD
  228.      */
  229.     public double getRMS() {
  230.         return rms;
  231.     }

  233.     /** Get the number of objective function evaluations.
  234.      *  @return the number of objective function evaluations.
  235.      */
  236.     public int getEvaluations() {
  237.         return optimizer.getEvaluations();
  238.     }

  240.     /** Get the function computing position/velocity at sample points.
  241.      * @return function computing position/velocity at sample points
  242.      */
  243.     protected abstract MultivariateVectorFunction getObjectiveFunction();

  245.     /** Get the Jacobian of the function computing position/velocity at sample points.
  246.      * @return Jacobian of the function computing position/velocity at sample points
  247.      */
  248.     protected abstract MultivariateMatrixFunction getObjectiveFunctionJacobian();

  250.     /** Check if fitting uses only sample positions.
  251.      * @return true if fitting uses only sample positions
  252.      */
  253.     protected boolean isOnlyPosition() {
  254.         return onlyPosition;
  255.     }

  257.     /** Get the size of the target.
  258.      * @return target size
  259.      */
  260.     protected int getTargetSize() {
  261.         return target.length;
  262.     }

  264.     /** Get the date of the initial state.
  265.      * @return the date
  266.      */
  267.     protected AbsoluteDate getDate() {
  268.         return date;
  269.     }

  271.     /** Get the frame of the initial state.
  272.      * @return the orbit frame
  273.      */
  274.     protected Frame getFrame() {
  275.         return frame;
  276.     }

  278.     /** Get the states sample.
  279.      * @return the states sample
  280.      */
  281.     protected List<SpacecraftState> getSample() {
  282.         return sample;
  283.     }

  285.     /** Get the free parameters.
  286.      * @return the free parameters
  287.      */
  288.     protected Collection<String>  getFreeParameters() {
  289.         return parameters;
  290.     }

  292.     /** Create a sample of {@link SpacecraftState}.
  293.      * @param source initial propagator
  294.      * @param timeSpan time span for the sample
  295.      * @param nbPoints number of points for the sample over the time span
  296.      * @return a sample of {@link SpacecraftState}
  297.      * @exception OrekitException if one of the sample point cannot be computed
  298.      */
  299.     private List<SpacecraftState> createSample(final Propagator source,
  300.                                                final double timeSpan,
  301.                                                final int nbPoints) throws OrekitException {

  303.         final SpacecraftState initialState = source.getInitialState();

  305.         final List<SpacecraftState> states = new ArrayList<SpacecraftState>();

  307.         final double stepSize = timeSpan / (nbPoints - 1);
  308.         final AbsoluteDate iniDate = source.getInitialState().getDate();
  309.         for (double dt = 0; dt < timeSpan; dt += stepSize) {
  310.             states.add(source.propagate(iniDate.shiftedBy(dt)));
  311.         }

  313.         source.resetInitialState(initialState);

  315.         return states;
  316.     }

  318.     /** Check if parameters can be free.
  319.      * @param freeParameters names of the free parameters
  320.      * @exception OrekitException if one of the parameters cannot be free
  321.      */
  322.     private void checkParameters(final Collection<String> freeParameters) throws OrekitException {
  323.         for (String parameter : freeParameters) {
  324.             checkParameter(parameter);
  325.         }
  326.     }

  328.     /** Check if parameters can be free.
  329.      * @param freeParameters names of the free parameters
  330.      * @exception OrekitException if one of the parameters cannot be free
  331.      */
  332.     private void checkParameters(final String ... freeParameters) throws OrekitException {
  333.         for (String parameter : freeParameters) {
  334.             checkParameter(parameter);
  335.         }
  336.     }

  338.     /** Check if parameter can be free.
  339.      * @param parameter name of the free parameter
  340.      * @exception OrekitException if the parameter cannot be free
  341.      */
  342.     private void checkParameter(final String parameter) throws OrekitException {

  344.         if ( !availableParameters.contains(parameter) ) {

  346.             // build the list of supported parameters
  347.             final StringBuilder sBuilder = new StringBuilder();
  348.             for (final String available : availableParameters) {
  349.                 if (sBuilder.length() > 0) {
  350.                     sBuilder.append(", ");
  351.                 }
  352.                 sBuilder.append(available);
  353.             }

  355.             throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME,
  356.                                       parameter, sBuilder.toString());

  358.         }

  360.     }

  362.     /** Adapt a propagator to minimize the mean square error for a set of {@link SpacecraftState states}.
  363.      * @param states set of spacecraft states to fit
  364.      * @param positionOnly if true, consider only position data otherwise both position and velocity are used
  365.      * @return adapted propagator
  366.      * @exception OrekitException if propagator cannot be adapted
  367.      */
  368.     private Propagator adapt(final List<SpacecraftState> states,
  369.                              final boolean positionOnly) throws OrekitException {

  371.         this.date = states.get(0).getDate();
  372.         this.onlyPosition = positionOnly;

  374.         // very rough first guess using osculating parameters of first sample point
  375.         final double[] initial = new double[6 + parameters.size()];
  376.         final PVCoordinates pv = states.get(0).getPVCoordinates(frame);
  377.         initial[0] = pv.getPosition().getX();
  378.         initial[1] = pv.getPosition().getY();
  379.         initial[2] = pv.getPosition().getZ();
  380.         initial[3] = pv.getVelocity().getX();
  381.         initial[4] = pv.getVelocity().getY();
  382.         initial[5] = pv.getVelocity().getZ();
  383.         int i = 6;
  384.         for (String name : parameters) {
  385.             initial[i++] = builder.getParameter(name);
  386.         }

  388.         // warm-up iterations, using only a few points
  389.         setSample(states.subList(0, onlyPosition ? 2 : 1));
  390.         final double[] intermediate = fit(initial);

  392.         // final search using all points
  393.         setSample(states);
  394.         final double[] result = fit(intermediate);

  396.         rms = getRMS(result);
  397.         adapted = buildAdaptedPropagator(result);

  399.         return adapted;
  400.     }

  402.     /** Find the propagator that minimize the mean square error for a sample of {@link SpacecraftState states}.
  403.      * @param initial initial estimation parameters (position, velocity, free parameters)
  404.      * @return fitted parameters
  405.      * @exception OrekitException if propagator cannot be adapted
  406.      * @exception MaxCountExceededException if maximal number of iterations is exceeded
  407.      */
  408.     private double[] fit(final double[] initial)
  409.         throws OrekitException, MaxCountExceededException {

  411.         final PointVectorValuePair optimum = optimizer.optimize(new MaxEval(maxIterations),
  412.                                                                 new ModelFunction(getObjectiveFunction()),
  413.                                                                 new ModelFunctionJacobian(getObjectiveFunctionJacobian()),
  414.                                                                 new Target(target),
  415.                                                                 new Weight(weight),
  416.                                                                 new InitialGuess(initial));
  417.         return optimum.getPoint();
  418.     }

  420.     /** Get the Root Mean Square Deviation for a given parameters set.
  421.      * @param parameterSet position/velocity parameters set
  422.      * @return RMSD
  423.      * @exception OrekitException if position/velocity cannot be computed at some date
  424.      */
  425.     private double getRMS(final double[] parameterSet) throws OrekitException {

  427.         final double[] residuals = getResiduals(parameterSet);
  428.         double sum2 = 0;
  429.         for (final double residual : residuals) {
  430.             sum2 += residual * residual;
  431.         }
  432.         return FastMath.sqrt(sum2 / residuals.length);

  434.     }

  436.     /** Get the residuals for a given position/velocity parameters set.
  437.      * @param parameterSet position/velocity parameters set
  438.      * @return residuals
  439.      * @exception OrekitException if position/velocity cannot be computed at some date
  440.      */
  441.     private double[] getResiduals(final double[] parameterSet) throws OrekitException {
  442.         try {
  443.             final double[] residuals = getObjectiveFunction().value(parameterSet);
  444.             for (int i = 0; i < residuals.length; ++i) {
  445.                 residuals[i] = target[i] - residuals[i];
  446.             }
  447.             return residuals;
  448.         } catch (OrekitExceptionWrapper oew) {
  449.             throw oew.getException();
  450.         }
  451.     }

  453.     /** Build the adpated propagator for a given position/velocity(/free) parameters set.
  454.      * @param parameterSet position/velocity(/free) parameters set
  455.      * @return adapted propagator
  456.      * @exception OrekitException if propagator cannot be build
  457.      */
  458.     private Propagator buildAdaptedPropagator(final double[] parameterSet)
  459.         throws OrekitException {
  460.         return builder.buildPropagator(date, parameterSet);
  461.     }

  463.     /** Set the states sample.
  464.      * @param states spacecraft states sample
  465.      * @exception OrekitException if position/velocity cannot be extracted from sample
  466.      */
  467.     private void setSample(final List<SpacecraftState> states) throws OrekitException {

  469.         this.sample = states;

  471.         if (onlyPosition) {
  472.             target = new double[states.size() * 3];
  473.             weight = new double[states.size() * 3];
  474.         } else {
  475.             target = new double[states.size() * 6];
  476.             weight = new double[states.size() * 6];
  477.         }

  479.         int k = 0;
  480.         for (int i = 0; i < states.size(); i++) {

  482.             final PVCoordinates pv = states.get(i).getPVCoordinates(frame);

  484.             // position
  485.             target[k]   = pv.getPosition().getX();
  486.             weight[k++] = 1;
  487.             target[k]   = pv.getPosition().getY();
  488.             weight[k++] = 1;
  489.             target[k]   = pv.getPosition().getZ();
  490.             weight[k++] = 1;

  492.             // velocity
  493.             if (!onlyPosition) {
  494.                 // velocity weight relative to position
  495.                 final double r2 = pv.getPosition().getNormSq();
  496.                 final double v  = pv.getVelocity().getNorm();
  497.                 final double vWeight = v * r2 / states.get(i).getMu();

  499.                 target[k]   = pv.getVelocity().getX();
  500.                 weight[k++] = vWeight;
  501.                 target[k]   = pv.getVelocity().getY();
  502.                 weight[k++] = vWeight;
  503.                 target[k]   = pv.getVelocity().getZ();
  504.                 weight[k++] = vWeight;
  505.             }

  507.         }

  509.     }

  511. }
Created with JaCoCo 0.6.3.201306030806