Differentiation.java

  1. /* Copyright 2002-2019 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.utils;

  19. import org.hipparchus.analysis.UnivariateFunction;
  20. import org.hipparchus.analysis.UnivariateVectorFunction;
  21. import org.hipparchus.analysis.differentiation.DSFactory;
  22. import org.hipparchus.analysis.differentiation.DerivativeStructure;
  23. import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
  24. import org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction;
  25. import org.orekit.attitudes.AttitudeProvider;
  26. import org.orekit.orbits.Orbit;
  27. import org.orekit.orbits.OrbitType;
  28. import org.orekit.orbits.PositionAngle;
  29. import org.orekit.propagation.SpacecraftState;
  30. import org.orekit.propagation.numerical.NumericalPropagator;

  32. /** Utility class for differentiating various kinds of functions.
  33.  * @author Luc Maisonobe
  34.  * @since 8.0
  35.  */
  36. public class Differentiation {

  38.     /** Factory for the DerivativeStructure instances. */
  39.     private static final DSFactory FACTORY = new DSFactory(1, 1);

  41.     /** Private constructor for utility class.
  42.      */
  43.     private Differentiation() {
  44.     }

  46.     /** Differentiate a scalar function using finite differences.
  47.      * @param function function to differentiate
  48.      * @param driver driver for the parameter
  49.      * @param nbPoints number of points used for finite differences
  50.      * @param normalizedStep step for finite differences, in <em>normalized</em> units
  51.      * @return scalar function evaluating to the derivative of the original function
  52.      * @deprecated as of 9.3, replaced by {@link #differentiate(ParameterFunction, int, double)}
  53.      */
  54.     @Deprecated
  55.     public static ParameterFunction differentiate(final ParameterFunction function, final ParameterDriver driver,
  56.                                                   final int nbPoints, final double normalizedStep) {
  57.         return differentiate(function, nbPoints, normalizedStep * driver.getScale());
  58.     }

  60.     /** Differentiate a scalar function using finite differences.
  61.      * @param function function to differentiate
  62.      * @param nbPoints number of points used for finite differences
  63.      * @param step step for finite differences, in <em>physical</em> units
  64.      * @return scalar function evaluating to the derivative of the original function
  65.      * @since 9.3
  66.      */
  67.     public static ParameterFunction differentiate(final ParameterFunction function,
  68.                                                   final int nbPoints, final double step) {

  70.         return new ParameterFunction() {

  72.             /** Finite differences differentiator to use. */
  73.             private final FiniteDifferencesDifferentiator differentiator  =
  74.                             new FiniteDifferencesDifferentiator(nbPoints, step);

  76.             /** {@inheritDoc} */
  77.             @Override
  78.             public double value(final ParameterDriver driver) {

  80.                 final UnivariateFunction uf = new UnivariateFunction() {
  81.                     /** {@inheritDoc} */
  82.                     @Override
  83.                     public double value(final double value) {
  84.                         final double saved = driver.getValue();
  85.                         driver.setValue(value);
  86.                         final double functionValue = function.value(driver);
  87.                         driver.setValue(saved);
  88.                         return functionValue;
  89.                     }
  90.                 };

  92.                 final DerivativeStructure dsParam = FACTORY.variable(0, driver.getValue());
  93.                 final DerivativeStructure dsValue = differentiator.differentiate(uf).value(dsParam);
  94.                 return dsValue.getPartialDerivative(1);

  96.             }
  97.         };

  99.     }

  101.     /** Differentiate a vector function using finite differences.
  102.      * @param function function to differentiate
  103.      * @param provider attitude provider to use for modified states
  104.      * @param dimension dimension of the vector value of the function
  105.      * @param orbitType type used to map the orbit to a one dimensional array
  106.      * @param positionAngle type of the position angle used for orbit mapping to array
  107.      * @param dP user specified position error, used for step size computation for finite differences
  108.      * @param nbPoints number of points used for finite differences
  109.      * @return matrix function evaluating to the Jacobian of the original function
  110.      */
  111.     public static StateJacobian differentiate(final StateFunction function, final int dimension,
  112.                                               final AttitudeProvider provider,
  113.                                               final OrbitType orbitType, final PositionAngle positionAngle,
  114.                                               final double dP, final int nbPoints) {
  115.         return new StateJacobian() {

  117.             @Override
  118.             public double[][] value(final SpacecraftState state) {
  119.                 final double[] tolerances =
  120.                         NumericalPropagator.tolerances(dP, state.getOrbit(), orbitType)[0];
  121.                 final double[][] jacobian = new double[dimension][6];
  122.                 for (int j = 0; j < 6; ++j) {

  124.                     // compute partial derivatives with respect to state component j
  125.                     final UnivariateVectorFunction componentJ =
  126.                             new StateComponentFunction(j, function, provider, state,
  127.                                                        orbitType, positionAngle);
  128.                     final FiniteDifferencesDifferentiator differentiator =
  129.                             new FiniteDifferencesDifferentiator(nbPoints, tolerances[j]);
  130.                     final UnivariateDifferentiableVectorFunction differentiatedJ =
  131.                             differentiator.differentiate(componentJ);

  133.                     final DerivativeStructure[] c = differentiatedJ.value(FACTORY.variable(0, 0.0));

  135.                     // populate the j-th column of the Jacobian
  136.                     for (int i = 0; i < dimension; ++i) {
  137.                         jacobian[i][j] = c[i].getPartialDerivative(1);
  138.                     }

  140.                 }

  142.                 return jacobian;

  144.             }

  146.         };
  147.     }

  149.     /** Restriction of a {@link StateFunction} to a function of a single state component.
  150.      */
  151.     private static class StateComponentFunction implements UnivariateVectorFunction {

  153.         /** Component index in the mapped orbit array. */
  154.         private final int             index;

  156.         /** State-dependent function. */
  157.         private final StateFunction   f;

  159.         /** Type used to map the orbit to a one dimensional array. */
  160.         private final OrbitType       orbitType;

  162.         /** Tpe of the position angle used for orbit mapping to array. */
  163.         private final PositionAngle   positionAngle;

  165.         /** Base state, of which only one component will change. */
  166.         private final SpacecraftState baseState;

  168.         /** Attitude provider to use for modified states. */
  169.         private final AttitudeProvider provider;

  171.         /** Simple constructor.
  172.          * @param index component index in the mapped orbit array
  173.          * @param f state-dependent function
  174.          * @param provider attitude provider to use for modified states
  175.          * @param baseState base state, of which only one component will change
  176.          * @param orbitType type used to map the orbit to a one dimensional array
  177.          * @param positionAngle type of the position angle used for orbit mapping to array
  178.          */
  179.         StateComponentFunction(final int index, final StateFunction f,
  180.                                final AttitudeProvider provider, final SpacecraftState baseState,
  181.                                final OrbitType orbitType, final PositionAngle positionAngle) {
  182.             this.index         = index;
  183.             this.f             = f;
  184.             this.provider      = provider;
  185.             this.orbitType     = orbitType;
  186.             this.positionAngle = positionAngle;
  187.             this.baseState     = baseState;
  188.         }

  190.         /** {@inheritDoc} */
  191.         @Override
  192.         public double[] value(final double x) {
  193.             final double[] array = new double[6];
  194.             final double[] arrayDot = new double[6];
  195.             orbitType.mapOrbitToArray(baseState.getOrbit(), positionAngle, array, arrayDot);
  196.             array[index] += x;
  197.             final Orbit orbit = orbitType.mapArrayToOrbit(array, arrayDot,
  198.                                                           positionAngle,
  199.                                                           baseState.getDate(),
  200.                                                           baseState.getMu(),
  201.                                                           baseState.getFrame());
  202.             final SpacecraftState state =
  203.                     new SpacecraftState(orbit,
  204.                                         provider.getAttitude(orbit, orbit.getDate(), orbit.getFrame()),
  205.                                         baseState.getMass());
  206.             return f.value(state);
  207.         }

  209.     }

  211. }

Created with JaCoCo 0.8.2.201808211720