1   /* Copyright 2022-2025 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.adjoint.cost;
18  
19  import org.hipparchus.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.util.FastMath;
21  
22  /**
23   * Fuel cost penalized with a logarithmic term, which is a barrier so is not defined for epsilon equal to 0 or 1.
24   *
25   * @author Romain Serra
26   * @since 13.0
27   */
28  public class LogarithmicBarrierCartesianFuel extends PenalizedCartesianFuelCost {
29  
30      /**
31       * Constructor.
32       *
33       * @param name                   adjoint name
34       * @param massFlowRateFactor     mass flow rate factor
35       * @param maximumThrustMagnitude maximum thrust magnitude
36       * @param epsilon                penalty weight
37       */
38      public LogarithmicBarrierCartesianFuel(final String name, final double massFlowRateFactor,
39                                             final double maximumThrustMagnitude, final double epsilon) {
40          super(name, massFlowRateFactor, maximumThrustMagnitude, epsilon);
41      }
42  
43      /** {@inheritDoc} */
44      @Override
45      public double evaluatePenaltyFunction(final double controlNorm) {
46          return FastMath.log(controlNorm) + FastMath.log(1. - controlNorm);
47      }
48  
49      /** {@inheritDoc} */
50      @Override
51      public Vector3D getThrustAccelerationVector(final double[] adjointVariables, final double mass) {
52          final double thrustForceMagnitude = getThrustForceMagnitude(adjointVariables, mass);
53          return getThrustDirection(adjointVariables).scalarMultiply(thrustForceMagnitude / mass);
54      }
55  
56      /** {@inheritDoc} */
57      @Override
58      public void updateAdjointDerivatives(final double[] adjointVariables, final double mass,
59                                           final double[] adjointDerivatives) {
60          if (getAdjointDimension() > 6) {
61              adjointDerivatives[6] += getAdjointVelocityNorm(adjointVariables) * getThrustForceMagnitude(adjointVariables,
62                      mass) / (mass * mass);
63          }
64      }
65  
66      /**
67       * Computes the Euclidean norm of the thrust force.
68       * @param adjointVariables adjoint variables
69       * @param mass mass
70       * @return thrust force magnitude
71       */
72      private double getThrustForceMagnitude(final double[] adjointVariables, final double mass) {
73          final double twoEpsilon = getEpsilon() * 2;
74          double otherTerm = getAdjointVelocityNorm(adjointVariables) / mass - 1;
75          if (getAdjointDimension() > 6) {
76              otherTerm -= getMassFlowRateFactor() * adjointVariables[6];
77          }
78          return twoEpsilon * getMaximumThrustMagnitude() / (twoEpsilon + otherTerm + FastMath.sqrt(otherTerm * otherTerm + twoEpsilon * twoEpsilon));
79      }
80  }