AbstractFixedBoundaryCartesianSingleShooting.java

  1. /* Copyright 2022-2024 Romain Serra
  2.  * Licensed to CS GROUP (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.control.indirect.shooting;

  19. import org.hipparchus.Field;
  20. import org.hipparchus.analysis.differentiation.Gradient;
  21. import org.hipparchus.analysis.differentiation.GradientField;
  22. import org.hipparchus.util.FastMath;
  23. import org.hipparchus.util.MathArrays;
  24. import org.orekit.attitudes.Attitude;
  25. import org.orekit.control.indirect.adjoint.CartesianAdjointDerivativesProvider;
  26. import org.orekit.control.indirect.adjoint.FieldCartesianAdjointDerivativesProvider;
  27. import org.orekit.control.indirect.shooting.boundary.CartesianBoundaryConditionChecker;
  28. import org.orekit.control.indirect.shooting.boundary.FixedTimeBoundaryOrbits;
  29. import org.orekit.control.indirect.shooting.boundary.FixedTimeCartesianBoundaryStates;
  30. import org.orekit.control.indirect.shooting.propagation.ShootingPropagationSettings;
  31. import org.orekit.frames.Frame;
  32. import org.orekit.orbits.CartesianOrbit;
  33. import org.orekit.orbits.Orbit;
  34. import org.orekit.propagation.FieldSpacecraftState;
  35. import org.orekit.propagation.SpacecraftState;
  36. import org.orekit.propagation.numerical.FieldNumericalPropagator;
  37. import org.orekit.propagation.numerical.NumericalPropagator;
  38. import org.orekit.time.FieldAbsoluteDate;
  39. import org.orekit.utils.AbsolutePVCoordinates;
  40. import org.orekit.utils.TimeStampedPVCoordinates;

  42. /**
  43.  * Abstract class for indirect single shooting methods with Cartesian coordinates for fixed time fixed boundary.
  44.  * Inheritors must implement the iteration update, assuming derivatives are needed.
  45.  * Terminal mass is assumed to be free, thus corresponding adjoint must vanish at terminal time.
  46.  * On the other hand, other terminal adjoint variables are free because the Cartesian state is fixed.
  47.  *
  48.  * @author Romain Serra
  49.  * @since 12.2
  50.  * @see org.orekit.control.indirect.adjoint.CartesianAdjointDerivativesProvider
  51.  * @see org.orekit.control.indirect.adjoint.FieldCartesianAdjointDerivativesProvider
  52.  */
  53. public abstract class AbstractFixedBoundaryCartesianSingleShooting extends AbstractIndirectShooting {

  55.     /** Default value for defects scaling. */
  56.     private static final double DEFAULT_SCALE = 1.;

  58.     /** Template for initial state. Contains the initial Cartesian coordinates. */
  59.     private final SpacecraftState initialSpacecraftStateTemplate;

  61.     /** Terminal Cartesian coordinates. */
  62.     private final TimeStampedPVCoordinates terminalCartesianState;

  64.     /** Condition checker. */
  65.     private final CartesianBoundaryConditionChecker conditionChecker;

  67.     /** Scale for velocity defects (m). */
  68.     private double scaleVelocityDefects;

  70.     /** Scale for position defects (m/s). */
  71.     private double scalePositionDefects;

  73.     /** Tolerance for convergence on terminal mass adjoint, if applicable to dynamics. */
  74.     private double toleranceMassAdjoint = DEFAULT_TOLERANCE_MASS_ADJOINT;

  76.     /**
  77.      * Constructor with boundary conditions as orbits.
  78.      * @param propagationSettings propagation settings
  79.      * @param boundaryConditions boundary conditions as {@link FixedTimeCartesianBoundaryStates}
  80.      * @param conditionChecker boundary condition checker
  81.      */
  82.     protected AbstractFixedBoundaryCartesianSingleShooting(final ShootingPropagationSettings propagationSettings,
  83.                                                            final FixedTimeCartesianBoundaryStates boundaryConditions,
  84.                                                            final CartesianBoundaryConditionChecker conditionChecker) {
  85.         super(propagationSettings);
  86.         this.conditionChecker = conditionChecker;
  87.         this.initialSpacecraftStateTemplate = buildInitialStateTemplate(boundaryConditions.getInitialCartesianState());
  88.         this.terminalCartesianState = boundaryConditions.getTerminalCartesianState().getPVCoordinates(propagationSettings.getPropagationFrame());
  89.         this.scalePositionDefects = DEFAULT_SCALE;
  90.         this.scaleVelocityDefects = DEFAULT_SCALE;
  91.     }

  93.     /**
  94.      * Constructor with boundary conditions as orbits.
  95.      * @param propagationSettings propagation settings
  96.      * @param boundaryConditions boundary conditions as {@link FixedTimeBoundaryOrbits}
  97.      * @param conditionChecker boundary condition checker
  98.      */
  99.     protected AbstractFixedBoundaryCartesianSingleShooting(final ShootingPropagationSettings propagationSettings,
  100.                                                            final FixedTimeBoundaryOrbits boundaryConditions,
  101.                                                            final CartesianBoundaryConditionChecker conditionChecker) {
  102.         super(propagationSettings);
  103.         this.conditionChecker = conditionChecker;
  104.         this.initialSpacecraftStateTemplate = buildInitialStateTemplate(boundaryConditions.getInitialOrbit());
  105.         this.terminalCartesianState = boundaryConditions.getTerminalOrbit().getPVCoordinates(propagationSettings.getPropagationFrame());
  106.         this.scalePositionDefects = DEFAULT_SCALE;
  107.         this.scaleVelocityDefects = DEFAULT_SCALE;
  108.     }

  110.     /**
  111.      * Setter for scale of position defects.
  112.      * @param scalePositionDefects new scale
  113.      */
  114.     public void setScalePositionDefects(final double scalePositionDefects) {
  115.         this.scalePositionDefects = scalePositionDefects;
  116.     }

  118.     /**
  119.      * Getter for scale of position defects.
  120.      * @return scale
  121.      */
  122.     public double getScalePositionDefects() {
  123.         return scalePositionDefects;
  124.     }

  126.     /**
  127.      * Setter for scale of velocity defects.
  128.      * @param scaleVelocityDefects new scale
  129.      */
  130.     public void setScaleVelocityDefects(final double scaleVelocityDefects) {
  131.         this.scaleVelocityDefects = scaleVelocityDefects;
  132.     }

  134.     /**
  135.      * Getter for scale of velocity defects.
  136.      * @return scale
  137.      */
  138.     public double getScaleVelocityDefects() {
  139.         return scaleVelocityDefects;
  140.     }

  142.     /**
  143.      * Getter for the boundary condition checker.
  144.      * @return checker
  145.      */
  146.     protected CartesianBoundaryConditionChecker getConditionChecker() {
  147.         return conditionChecker;
  148.     }

  150.     /**
  151.      * Getter for the target terminal Cartesian state vector.
  152.      * @return expected terminal state
  153.      */
  154.     protected TimeStampedPVCoordinates getTerminalCartesianState() {
  155.         return terminalCartesianState;
  156.     }

  158.     /**
  159.      * Setter for mass adjoint tolerance.
  160.      * @param toleranceMassAdjoint new tolerance value
  161.      */
  162.     public void setToleranceMassAdjoint(final double toleranceMassAdjoint) {
  163.         this.toleranceMassAdjoint = FastMath.abs(toleranceMassAdjoint);
  164.     }

  166.     /**
  167.      * Create template initial state (without adjoint varialbles) for propagation from orbits.
  168.      * @param initialOrbit initial orbit
  169.      * @return template propagation state
  170.      */
  171.     private SpacecraftState buildInitialStateTemplate(final Orbit initialOrbit) {
  172.         final Frame frame = getPropagationSettings().getPropagationFrame();
  173.         final CartesianOrbit cartesianOrbit = new CartesianOrbit(initialOrbit.getPVCoordinates(frame), frame,
  174.             initialOrbit.getDate(), initialOrbit.getMu());
  175.         final Attitude attitude = getPropagationSettings().getAttitudeProvider()
  176.                 .getAttitude(cartesianOrbit, cartesianOrbit.getDate(), cartesianOrbit.getFrame());
  177.         return new SpacecraftState(cartesianOrbit, attitude);
  178.     }

  180.     /**
  181.      * Create template initial state (without adjoint varialbles) for propagation.
  182.      * @param initialCartesianState initial Cartesian state
  183.      * @return template propagation state
  184.      */
  185.     private SpacecraftState buildInitialStateTemplate(final AbsolutePVCoordinates initialCartesianState) {
  186.         final Frame frame = getPropagationSettings().getPropagationFrame();
  187.         final AbsolutePVCoordinates absolutePVCoordinates = new AbsolutePVCoordinates(frame,
  188.                 initialCartesianState.getPVCoordinates(frame));
  189.         final Attitude attitude = getPropagationSettings().getAttitudeProvider()
  190.                 .getAttitude(absolutePVCoordinates, absolutePVCoordinates.getDate(), absolutePVCoordinates.getFrame());
  191.         return new SpacecraftState(absolutePVCoordinates, attitude);
  192.     }

  194.     /** {@inheritDoc} */
  195.     @Override
  196.     public ShootingBoundaryOutput solve(final double initialMass, final double[] initialGuess) {
  197.         // check initial guess
  198.         final SpacecraftState initialState = createStateWithMassAndAdjoint(initialMass, initialGuess);
  199.         final ShootingBoundaryOutput initialGuessSolution = computeCandidateSolution(initialState, 0);
  200.         if (initialGuessSolution.isConverged()) {
  201.             return initialGuessSolution;
  202.         } else {
  203.             return iterate(initialMass, initialGuess);
  204.         }
  205.     }

  207.     /** {@inheritDoc} */
  208.     @Override
  209.     protected NumericalPropagator buildPropagator(final SpacecraftState initialState) {
  210.         final NumericalPropagator propagator = super.buildPropagator(initialState);
  211.         final CartesianAdjointDerivativesProvider derivativesProvider = (CartesianAdjointDerivativesProvider)
  212.             getPropagationSettings().getAdjointDynamicsProvider().buildAdditionalDerivativesProvider();
  213.         derivativesProvider.getCost().getEventDetectors().forEach(propagator::addEventDetector);
  214.         return propagator;
  215.     }

  217.     /** {@inheritDoc} */
  218.     @Override
  219.     protected FieldNumericalPropagator<Gradient> buildFieldPropagator(final FieldSpacecraftState<Gradient> initialState) {
  220.         final FieldNumericalPropagator<Gradient> fieldPropagator = super.buildFieldPropagator(initialState);
  221.         final Field<Gradient> field = fieldPropagator.getField();
  222.         final FieldCartesianAdjointDerivativesProvider<Gradient> derivativesProvider = (FieldCartesianAdjointDerivativesProvider<Gradient>)
  223.             getPropagationSettings().getAdjointDynamicsProvider().buildFieldAdditionalDerivativesProvider(field);
  224.         derivativesProvider.getCost().getFieldEventDetectors(field).forEach(fieldPropagator::addEventDetector);
  225.         return fieldPropagator;
  226.     }

  228.     /**
  229.      * Form solution container with input initial state.
  230.      * @param initialState initial state
  231.      * @param iterationCount iteration count
  232.      * @return candidate solution
  233.      */
  234.     private ShootingBoundaryOutput computeCandidateSolution(final SpacecraftState initialState,
  235.                                                             final int iterationCount) {
  236.         final NumericalPropagator propagator = buildPropagator(initialState);
  237.         final SpacecraftState actualTerminalState = propagator.propagate(getTerminalCartesianState().getDate());
  238.         final boolean converged = checkConvergence(actualTerminalState);
  239.         return new ShootingBoundaryOutput(converged, iterationCount, initialState, getPropagationSettings(), actualTerminalState);
  240.     }

  242.     /**
  243.      * Iterate on initial guess to solve boundary problem.
  244.      * @param initialMass initial mass
  245.      * @param initialGuess initial guess for initial adjoint variables
  246.      * @return solution (or null pointer if not converged)
  247.      */
  248.     private ShootingBoundaryOutput iterate(final double initialMass, final double[] initialGuess) {
  249.         double[] initialAdjoint = initialGuess.clone();
  250.         int iterationCount = 0;
  251.         int maximumIterationCount = getConditionChecker().getMaximumIterationCount();
  252.         SpacecraftState initialState = createStateWithMassAndAdjoint(initialMass, initialGuess);
  253.         while (iterationCount < maximumIterationCount) {
  254.             final FieldSpacecraftState<Gradient> fieldInitialState = createFieldStateWithMassAndAdjoint(initialMass,
  255.                 initialAdjoint);
  256.             final Field<Gradient> field = fieldInitialState.getDate().getField();
  257.             final FieldAbsoluteDate<Gradient> fieldTerminalDate = new FieldAbsoluteDate<>(field, getTerminalCartesianState().getDate());
  258.             final FieldNumericalPropagator<Gradient> fieldPropagator = buildFieldPropagator(fieldInitialState);
  259.             final FieldSpacecraftState<Gradient> fieldTerminalState = fieldPropagator.propagate(fieldTerminalDate);
  260.             initialState = fieldInitialState.toSpacecraftState();
  261.             if (checkConvergence(fieldTerminalState.toSpacecraftState())) {
  262.                 // make sure non-Field version also satisfies convergence criterion
  263.                 final ShootingBoundaryOutput solution = computeCandidateSolution(initialState, iterationCount);
  264.                 if (solution.isConverged()) {
  265.                     return solution;
  266.                 }
  267.             }
  268.             initialAdjoint = updateAdjoint(initialAdjoint, fieldTerminalState);
  269.             if (Double.isNaN(initialAdjoint[0])) {
  270.                 return computeCandidateSolution(initialState, iterationCount);
  271.             }
  272.             iterationCount++;
  273.             maximumIterationCount = getConditionChecker().getMaximumIterationCount();
  274.         }
  275.         return computeCandidateSolution(initialState, maximumIterationCount);
  276.     }

  278.     /**
  279.      * Checks convergence.
  280.      * @param actualTerminalState achieved terminal state
  281.      * @return convergence flag
  282.      */
  283.     private boolean checkConvergence(final SpacecraftState actualTerminalState) {
  284.         final boolean isCartesianConverged = getConditionChecker().isConverged(getTerminalCartesianState(),
  285.                 actualTerminalState.getPVCoordinates());
  286.         if (isCartesianConverged) {
  287.             final String adjointName = getPropagationSettings().getAdjointDynamicsProvider().getAdjointName();
  288.             final double[] terminalAdjoint = actualTerminalState.getAdditionalState(adjointName);
  289.             if (terminalAdjoint.length == 7) {
  290.                 return FastMath.abs(terminalAdjoint[6]) < toleranceMassAdjoint;
  291.             } else {
  292.                 return true;
  293.             }
  294.         } else {
  295.             return false;
  296.         }
  297.     }

  299.     /**
  300.      * Create initial state with input mass and adjoint vector.
  301.      * @param initialMass initial mass
  302.      * @param initialAdjoint initial adjoint variables
  303.      * @return state
  304.      */
  305.     protected SpacecraftState createStateWithMassAndAdjoint(final double initialMass, final double[] initialAdjoint) {
  306.         final String adjointName = getPropagationSettings().getAdjointDynamicsProvider().getAdjointName();
  307.         return createStateWithMass(initialMass).addAdditionalState(adjointName, initialAdjoint);
  308.     }

  310.     /**
  311.      * Create initial state with input mass.
  312.      * @param initialMass initial mass
  313.      * @return state
  314.      */
  315.     private SpacecraftState createStateWithMass(final double initialMass) {
  316.         if (initialSpacecraftStateTemplate.isOrbitDefined()) {
  317.             return new SpacecraftState(initialSpacecraftStateTemplate.getOrbit(),
  318.                 initialSpacecraftStateTemplate.getAttitude(), initialMass);
  319.         } else {
  320.             return new SpacecraftState(initialSpacecraftStateTemplate.getAbsPVA(),
  321.                 initialSpacecraftStateTemplate.getAttitude(), initialMass);
  322.         }
  323.     }

  325.     /**
  326.      * Create initial state with input mass and adjoint vector.
  327.      * @param initialMass initial mass
  328.      * @param initialAdjoint initial adjoint variables
  329.      * @return state
  330.      */
  331.     protected FieldSpacecraftState<Gradient> createFieldStateWithMassAndAdjoint(final double initialMass,
  332.                                                                                 final double[] initialAdjoint) {
  333.         final int parametersNumber = initialAdjoint.length;
  334.         final GradientField field = GradientField.getField(parametersNumber);
  335.         final FieldSpacecraftState<Gradient> stateWithoutAdjoint = new FieldSpacecraftState<>(field, createStateWithMass(initialMass));
  336.         final Gradient[] fieldInitialAdjoint = MathArrays.buildArray(field, initialAdjoint.length);
  337.         for (int i = 0; i < parametersNumber; i++) {
  338.             fieldInitialAdjoint[i] = Gradient.variable(parametersNumber, i, initialAdjoint[i]);
  339.         }
  340.         final String adjointName = getPropagationSettings().getAdjointDynamicsProvider().getAdjointName();
  341.         return stateWithoutAdjoint.addAdditionalState(adjointName, fieldInitialAdjoint);
  342.     }

  344.     /**
  345.      * Update initial adjoint vector.
  346.      * @param originalInitialAdjoint original initial adjoint (before update)
  347.      * @param fieldTerminalState final state of gradient propagation
  348.      * @return updated initial adjoint vector
  349.      */
  350.     protected abstract double[] updateAdjoint(double[] originalInitialAdjoint,
  351.                                               FieldSpacecraftState<Gradient> fieldTerminalState);
  352. }
Created with JaCoCo 0.8.12.202403310830