Merge branch 'issue-1644' into 'develop'

add fight duration cost for indirect control Closes #1644 See merge request orekit/orekit!719
CS-SI · Jan 14, 2025 · 4eb2c4c · 4eb2c4c
2 parents 15e61bf + f2cee3e
commit 4eb2c4c
Show file tree

Hide file tree

Showing 35 changed files with 798 additions and 104 deletions.
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
@@ -21,6 +21,9 @@
   </properties>
   <body>
     <release version="13.0" date="TBD" description="TBD">
+        <action dev="serrof" type="add" issue="1644">
+            Add duration flight cost for indirect shooting.
+        </action>
         <action dev="serrof" type="add" issue="1641">
             Add tuning option for LU decomposition in Newton-based indirect shooting.
         </action>

diff --git a/src/main/java/org/orekit/control/indirect/adjoint/cost/AbstractCartesianCost.java b/src/main/java/org/orekit/control/indirect/adjoint/cost/AbstractCartesianCost.java
@@ -35,6 +35,9 @@ public abstract class AbstractCartesianCost implements CartesianCost {
     /** Mass flow rate factor (always positive). */
     private final double massFlowRateFactor;
 
+    /** Dimension of adjoint vector. */
+    private final int adjointDimension;
+
     /**
      * Constructor.
      * @param name name
@@ -43,6 +46,15 @@ public abstract class AbstractCartesianCost implements CartesianCost {
     protected AbstractCartesianCost(final String name, final double massFlowRateFactor) {
         this.name = name;
         this.massFlowRateFactor = FastMath.abs(massFlowRateFactor);
+        this.adjointDimension = this.massFlowRateFactor == 0. ? 6 : 7;
+    }
+
+    /**
+     * {@inheritDoc}
+     */
+    @Override
+    public int getAdjointDimension() {
+        return adjointDimension;
     }
 
     /**

diff --git a/src/main/java/org/orekit/control/indirect/adjoint/cost/BoundedCartesianEnergy.java b/src/main/java/org/orekit/control/indirect/adjoint/cost/BoundedCartesianEnergy.java
@@ -74,7 +74,10 @@ public double getMaximumThrustMagnitude() {
     @Override
     protected double getThrustForceNorm(final double[] adjointVariables, final double mass) {
         final double adjointVelocityNorm = getAdjointVelocityNorm(adjointVariables);
-        final double factor = adjointVelocityNorm / mass - getMassFlowRateFactor() * adjointVariables[6];
+        double factor = adjointVelocityNorm / mass;
+        if (getAdjointDimension() > 6) {
+            factor -= getMassFlowRateFactor() * adjointVariables[6];
+        }
         if (factor > maximumThrustMagnitude) {
             return maximumThrustMagnitude;
         } else {

diff --git a/src/main/java/org/orekit/control/indirect/adjoint/cost/CartesianCost.java b/src/main/java/org/orekit/control/indirect/adjoint/cost/CartesianCost.java
@@ -38,12 +38,10 @@ public interface CartesianCost {
      */
     String getAdjointName();
 
-    /** Getter for adjoint vector dimension. Default is 7 (six for Cartesian coordinates and one for mass).
+    /** Getter for adjoint vector dimension.
      * @return adjoint dimension
      */
-    default int getAdjointDimension() {
-        return 7;
-    }
+    int getAdjointDimension();
 
     /** Getter for mass flow rate factor. It is negated and multiplied by the thrust force magnitude to obtain the mass time derivative.
      * The fact that it is a constant means that the exhaust speed is assumed to be independent of time.

diff --git a/src/main/java/org/orekit/control/indirect/adjoint/cost/CartesianEnergyConsideringMass.java b/src/main/java/org/orekit/control/indirect/adjoint/cost/CartesianEnergyConsideringMass.java
@@ -22,7 +22,7 @@
 import org.orekit.propagation.events.EventDetectionSettings;
 
 /**
- * Abstract class for energy cost with Cartesian coordinates and non-zero mass flow rate.
+ * Abstract class for energy cost with Cartesian coordinates.
  * An energy cost is proportional to the integral over time of the squared Euclidean norm of the control vector, often scaled with 1/2.
  * This type of cost is not optimal in terms of mass consumption, however its solutions showcase a smoother behavior favorable for convergence in shooting techniques.
  *
@@ -81,7 +81,9 @@ protected Vector3D getThrustDirection(final double[] adjointVariables) {
     @Override
     public void updateAdjointDerivatives(final double[] adjointVariables, final double mass,
                                          final double[] adjointDerivatives) {
-        adjointDerivatives[6] += getThrustForceNorm(adjointVariables, mass) * getAdjointVelocityNorm(adjointVariables) / (mass * mass);
+        if (getAdjointDimension() > 6) {
+            adjointDerivatives[6] += getThrustForceNorm(adjointVariables, mass) * getAdjointVelocityNorm(adjointVariables) / (mass * mass);
+        }
     }
 
     /** {@inheritDoc} */
@@ -124,7 +126,11 @@ public double g(final SpacecraftState state) {
          */
         private double evaluateVariablePart(final double[] adjointVariables, final double mass) {
             final double adjointVelocityNorm = getAdjointVelocityNorm(adjointVariables);
-            return adjointVelocityNorm / mass - getMassFlowRateFactor() * adjointVariables[6];
+            double variablePart = adjointVelocityNorm / mass;
+            if (getAdjointDimension() > 6) {
+                variablePart -= getMassFlowRateFactor() * adjointVariables[6];
+            }
+            return variablePart;
         }
 
     }

diff --git a/src/main/java/org/orekit/control/indirect/adjoint/cost/CartesianFlightDurationCost.java b/src/main/java/org/orekit/control/indirect/adjoint/cost/CartesianFlightDurationCost.java
@@ -0,0 +1,88 @@
+/* Copyright 2022-2025 Romain Serra
+ * Licensed to CS GROUP (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.control.indirect.adjoint.cost;
+
+
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+
+/**
+ * Class for minimizing the flight duration (a.k.a. time of flight) with Cartesian coordinates.
+ * It is the integral over time of the constant one. The control is assumed to be bounded.
+ * It also assumes that no external acceleration depends on mass.
+ * If the mass flow rate factor is zero, then there is no adjoint for the mass.
+ *
+ * @author Romain Serra
+ * @see CartesianCost
+ * @since 13.0
+ */
+public class CartesianFlightDurationCost extends AbstractCartesianCost {
+
+    /**
+     * Maximum value of thrust force Euclidean norm.
+     */
+    private final double maximumThrustMagnitude;
+
+    /**
+     * Constructor.
+     *
+     * @param name                   name
+     * @param massFlowRateFactor     mass flow rate factor
+     * @param maximumThrustMagnitude maximum thrust magnitude
+     */
+    public CartesianFlightDurationCost(final String name, final double massFlowRateFactor,
+                                       final double maximumThrustMagnitude) {
+        super(name, massFlowRateFactor);
+        this.maximumThrustMagnitude = maximumThrustMagnitude;
+    }
+
+    /**
+     * Getter for maximum thrust magnitude.
+     *
+     * @return maximum thrust
+     */
+    public double getMaximumThrustMagnitude() {
+        return maximumThrustMagnitude;
+    }
+
+    /**
+     * {@inheritDoc}
+     */
+    @Override
+    public Vector3D getThrustAccelerationVector(final double[] adjointVariables, final double mass) {
+        return new Vector3D(adjointVariables[3], adjointVariables[4], adjointVariables[5]).normalize()
+                .scalarMultiply(maximumThrustMagnitude);
+    }
+
+    /**
+     * {@inheritDoc}
+     */
+    @Override
+    public void updateAdjointDerivatives(final double[] adjointVariables, final double mass,
+                                         final double[] adjointDerivatives) {
+        if (getAdjointDimension() > 6) {
+            adjointDerivatives[6] += getAdjointVelocityNorm(adjointVariables) * maximumThrustMagnitude / (mass * mass);
+        }
+    }
+
+    /**
+     * {@inheritDoc}
+     */
+    @Override
+    public double getHamiltonianContribution(final double[] adjointVariables, final double mass) {
+        return -1.;
+    }
+}
diff --git a/src/main/java/org/orekit/control/indirect/adjoint/cost/CartesianFuelCost.java b/src/main/java/org/orekit/control/indirect/adjoint/cost/CartesianFuelCost.java
@@ -86,7 +86,11 @@ public EventDetectionSettings getEventDetectionSettings() {
      * @return value of switch function
      */
     private double evaluateSwitchFunction(final double[] adjointVariables, final double mass) {
-        return getAdjointVelocityNorm(adjointVariables) / mass - adjointVariables[6] * getMassFlowRateFactor() - 1.;
+        double switchFunction = getAdjointVelocityNorm(adjointVariables) / mass - 1.;
+        if (getAdjointDimension() > 6) {
+            switchFunction -= adjointVariables[6] * getMassFlowRateFactor();
+        }
+        return switchFunction;
     }
 
     /**
@@ -113,9 +117,11 @@ public Vector3D getThrustAccelerationVector(final double[] adjointVariables, fin
     @Override
     public void updateAdjointDerivatives(final double[] adjointVariables, final double mass,
                                          final double[] adjointDerivatives) {
-        final double switchFunction = evaluateSwitchFunction(adjointVariables, mass);
-        if (switchFunction > 0.) {
-            adjointDerivatives[6] += getAdjointVelocityNorm(adjointVariables) * maximumThrustMagnitude / (mass * mass);
+        if (getAdjointDimension() > 6) {
+            final double switchFunction = evaluateSwitchFunction(adjointVariables, mass);
+            if (switchFunction > 0.) {
+                adjointDerivatives[6] += getAdjointVelocityNorm(adjointVariables) * maximumThrustMagnitude / (mass * mass);
+            }
         }
     }
 

diff --git a/src/main/java/org/orekit/control/indirect/adjoint/cost/FieldAbstractCartesianCost.java b/src/main/java/org/orekit/control/indirect/adjoint/cost/FieldAbstractCartesianCost.java
@@ -36,6 +36,9 @@ public abstract class FieldAbstractCartesianCost<T extends CalculusFieldElement<
     /** Mass flow rate factor (always positive). */
     private final T massFlowRateFactor;
 
+    /** Dimension of adjoint vector. */
+    private final int adjointDimension;
+
     /**
      * Constructor.
      * @param name name
@@ -44,6 +47,13 @@ public abstract class FieldAbstractCartesianCost<T extends CalculusFieldElement<
     protected FieldAbstractCartesianCost(final String name, final T massFlowRateFactor) {
         this.name = name;
         this.massFlowRateFactor = FastMath.abs(massFlowRateFactor);
+        this.adjointDimension = this.massFlowRateFactor.isZero() ? 6 : 7;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public int getAdjointDimension() {
+        return adjointDimension;
     }
 
     /**

diff --git a/src/main/java/org/orekit/control/indirect/adjoint/cost/FieldBoundedCartesianEnergy.java b/src/main/java/org/orekit/control/indirect/adjoint/cost/FieldBoundedCartesianEnergy.java
@@ -79,7 +79,10 @@ public T getMaximumThrustMagnitude() {
     @Override
     protected T getFieldThrustForceNorm(final T[] adjointVariables, final T mass) {
         final T adjointVelocityNorm = getFieldAdjointVelocityNorm(adjointVariables);
-        final T factor = adjointVelocityNorm.divide(mass).subtract(adjointVariables[6].multiply(getMassFlowRateFactor()));
+        T factor = adjointVelocityNorm.divide(mass);
+        if (getAdjointDimension() > 6) {
+            factor = factor.subtract(adjointVariables[6].multiply(getMassFlowRateFactor()));
+        }
         final double factorReal = factor.getReal();
         final T zero = mass.getField().getZero();
         if (factorReal > maximumThrustMagnitude.getReal()) {

diff --git a/src/main/java/org/orekit/control/indirect/adjoint/cost/FieldCartesianCost.java b/src/main/java/org/orekit/control/indirect/adjoint/cost/FieldCartesianCost.java
@@ -39,12 +39,10 @@ public interface FieldCartesianCost<T extends CalculusFieldElement<T>> {
      */
     String getAdjointName();
 
-    /** Getter for adjoint vector dimension. Default is 7 (six for Cartesian coordinates and one for mass).
+    /** Getter for adjoint vector dimension.
      * @return adjoint dimension
      */
-    default int getAdjointDimension() {
-        return 7;
-    }
+    int getAdjointDimension();
 
     /** Getter for mass flow rate factor. It is negated and multiplied by the thrust force magnitude to obtain the mass time derivative.
      * The fact that it is a constant means that the exhaust speed is assumed to be independent of time.

diff --git a/...in/java/org/orekit/control/indirect/adjoint/cost/FieldCartesianEnergyConsideringMass.java b/...in/java/org/orekit/control/indirect/adjoint/cost/FieldCartesianEnergyConsideringMass.java
@@ -19,7 +19,6 @@
 
 import org.hipparchus.CalculusFieldElement;
 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
-import org.hipparchus.util.FastMath;
 import org.orekit.propagation.FieldSpacecraftState;
 import org.orekit.propagation.events.FieldEventDetectionSettings;
 
@@ -34,13 +33,7 @@
  * @see CartesianEnergyConsideringMass
  * @since 13.0
  */
-abstract class FieldCartesianEnergyConsideringMass<T extends CalculusFieldElement<T>> implements FieldCartesianCost<T> {
-
-    /** Name of adjoint vector. */
-    private final String name;
-
-    /** Mass flow rate factor (always positive). */
-    private final T massFlowRateFactor;
+abstract class FieldCartesianEnergyConsideringMass<T extends CalculusFieldElement<T>> extends FieldAbstractCartesianCost<T> {
 
     /** Detection settings for singularity detection. */
     private final FieldEventDetectionSettings<T> eventDetectionSettings;
@@ -53,26 +46,10 @@ abstract class FieldCartesianEnergyConsideringMass<T extends CalculusFieldElemen
      */
     protected FieldCartesianEnergyConsideringMass(final String name, final T massFlowRateFactor,
                                                   final FieldEventDetectionSettings<T> eventDetectionSettings) {
-        this.name = name;
-        this.massFlowRateFactor = massFlowRateFactor;
+        super(name, massFlowRateFactor);
         this.eventDetectionSettings = eventDetectionSettings;
     }
 
-    /**
-     * Getter for adjoint vector name.
-     * @return name
-     */
-    @Override
-    public String getAdjointName() {
-        return name;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public T getMassFlowRateFactor() {
-        return massFlowRateFactor;
-    }
-
     /**
      * Getter for event detection settings.
      * @return detection settings.
@@ -106,19 +83,11 @@ protected FieldVector3D<T> getFieldThrustDirection(final T[] adjointVariables) {
 
     /** {@inheritDoc} */
     @Override
-    public void updateFieldAdjointDerivatives(final T[] adjointVariables, final T mass,
-                                                                                  final T[] adjointDerivatives) {
-        adjointDerivatives[6] = adjointDerivatives[6].add(getFieldThrustForceNorm(adjointVariables, mass)
-            .multiply(getFieldAdjointVelocityNorm(adjointVariables)).divide(mass.square()));
-    }
-
-    /**
-     * Computes the Euclidean norm of the adjoint velocity vector.
-     * @param adjointVariables adjoint vector
-     * @return norm of adjoint velocity
-     */
-    protected T getFieldAdjointVelocityNorm(final T[] adjointVariables) {
-        return FastMath.sqrt(adjointVariables[3].square().add(adjointVariables[4].square()).add(adjointVariables[5].square()));
+    public void updateFieldAdjointDerivatives(final T[] adjointVariables, final T mass, final T[] adjointDerivatives) {
+        if (getAdjointDimension() > 6) {
+            adjointDerivatives[6] = adjointDerivatives[6].add(getFieldThrustForceNorm(adjointVariables, mass)
+                    .multiply(getFieldAdjointVelocityNorm(adjointVariables)).divide(mass.square()));
+        }
     }
 
     /** {@inheritDoc} */
@@ -161,7 +130,11 @@ public T g(final FieldSpacecraftState<T> state) {
          */
         private T evaluateVariablePart(final T[] adjointVariables, final T mass) {
             final T adjointVelocityNorm = getFieldAdjointVelocityNorm(adjointVariables);
-            return adjointVelocityNorm.divide(mass).subtract(adjointVariables[6].multiply(getMassFlowRateFactor()));
+            T variablePart = adjointVelocityNorm.divide(mass);
+            if (getAdjointDimension() > 6) {
+                variablePart = variablePart.subtract(adjointVariables[6].multiply(getMassFlowRateFactor()));
+            }
+            return variablePart;
         }
 
     }