1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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 import org.junit.jupiter.api.Assertions;
22 import org.junit.jupiter.api.Test;
23 import org.junit.jupiter.params.ParameterizedTest;
24 import org.junit.jupiter.params.provider.ValueSource;
25 import org.orekit.control.indirect.adjoint.CartesianAdjointDerivativesProvider;
26 import org.orekit.frames.FramesFactory;
27 import org.orekit.orbits.CartesianOrbit;
28 import org.orekit.propagation.SpacecraftState;
29 import org.orekit.propagation.events.EventDetector;
30 import org.orekit.propagation.integration.AdditionalDerivativesProvider;
31 import org.orekit.propagation.integration.CombinedDerivatives;
32 import org.orekit.time.AbsoluteDate;
33 import org.orekit.utils.PVCoordinates;
34
35 import java.util.List;
36 import java.util.stream.Collectors;
37 import java.util.stream.Stream;
38
39 class QuadraticPenaltyCartesianFuelTest {
40
41 private static final String ADJOINT_NAME = "adjoint";
42
43 @ParameterizedTest
44 @ValueSource(booleans = {false, true})
45 void testUpdateFieldAdjointDerivatives(final boolean withMass) {
46
47 final double massFlowRateFactor = withMass ? 1 : 0;
48 final QuadraticPenaltyCartesianFuel cost = new QuadraticPenaltyCartesianFuel("adjoint", massFlowRateFactor, 2, 0.5);
49 final double[] adjoint = new double[withMass ? 7 : 6];
50 adjoint[3] = 1;
51 final double[] derivatives = new double[adjoint.length];
52
53 cost.updateAdjointDerivatives(adjoint, 1, derivatives);
54
55 for (int i = 0; i < 6; ++i) {
56 Assertions.assertEquals(0., derivatives[i]);
57 }
58 if (withMass) {
59 Assertions.assertNotEquals(0., derivatives[derivatives.length - 1]);
60 } else {
61 Assertions.assertEquals(0., derivatives[derivatives.length - 1]);
62 }
63 }
64
65 @ParameterizedTest
66 @ValueSource(doubles = {0., 0.1, 0.5, 0.9})
67 void testEvaluatePenaltyFunction(final double norm) {
68
69 final double unitMaximumThrust = 1;
70 final QuadraticPenaltyCartesianFuel penalizedCartesianFuel = new QuadraticPenaltyCartesianFuel(ADJOINT_NAME,
71 2., unitMaximumThrust, 0.5);
72
73 final double actualPenalty = penalizedCartesianFuel.evaluatePenaltyFunction(norm);
74
75 Assertions.assertEquals(norm * norm / 2 - norm, actualPenalty, 1e-15);
76 }
77
78 @Test
79 void testGetEventDetectors() {
80
81 final QuadraticPenaltyCartesianFuel penalizedCartesianFuel = new QuadraticPenaltyCartesianFuel(ADJOINT_NAME,
82 1., 2., 0.5);
83
84 final Stream<EventDetector> actual = penalizedCartesianFuel.getEventDetectors();
85
86 final List<EventDetector> eventDetectorList = actual.collect(Collectors.toList());
87 Assertions.assertEquals(2, eventDetectorList.size());
88 final SpacecraftState state = buildState(new double[] {0, 0, 0, 1, 2, 3, 4}, 10);
89 final double g1 = eventDetectorList.get(0).g(state);
90 final double g2 = eventDetectorList.get(1).g(state);
91 Assertions.assertEquals(penalizedCartesianFuel.getMaximumThrustMagnitude(), FastMath.abs(g2 - g1));
92 }
93
94 @Test
95 void testGetThrustAccelerationVectorEpsilonCloseToZero() {
96
97 final double massFlowRateFactor = 1.;
98 final double maximumThrustMagnitude = 10.;
99 final double epsilon = 1e-6;
100 final QuadraticPenaltyCartesianFuel penalizedCartesianFuel = new QuadraticPenaltyCartesianFuel(ADJOINT_NAME,
101 massFlowRateFactor, maximumThrustMagnitude, epsilon);
102 final double mass = 100;
103 final double[] adjoint = new double[] {1, 2, 3, 4, 5, 6, 7};
104
105 final Vector3D actualThrustVector = penalizedCartesianFuel.getThrustAccelerationVector(adjoint, mass);
106
107 final CartesianFuelCost fuelCost = new CartesianFuelCost(ADJOINT_NAME, massFlowRateFactor, maximumThrustMagnitude);
108 final Vector3D expectedThrustVector = fuelCost.getThrustAccelerationVector(adjoint, mass);
109 Assertions.assertEquals(0., expectedThrustVector.subtract(actualThrustVector).getNorm(), 1e-10);
110 }
111
112 @ParameterizedTest
113 @ValueSource(doubles = {0, 1})
114 void testAgainstBoundedCartesianEnergyVaryingMassFlowRateFactor(final double massFlowRateFactor) {
115 testTemplateAgainstBoundedCartesianEnergy(1, massFlowRateFactor);
116 }
117
118 @ParameterizedTest
119 @ValueSource(doubles = {1, 1e2, 1e4})
120 void testAgainstBoundedCartesianEnergyVaryingMass(final double mass) {
121 testTemplateAgainstBoundedCartesianEnergy(mass, 1e-2);
122 }
123
124 private void testTemplateAgainstBoundedCartesianEnergy(final double mass, final double massFlowRateFactor) {
125
126 final double maximumThrustMagnitude = 1e-1;
127 final double epsilon = 1;
128 final QuadraticPenaltyCartesianFuel penalizedCartesianFuel = new QuadraticPenaltyCartesianFuel(ADJOINT_NAME,
129 massFlowRateFactor, maximumThrustMagnitude, epsilon);
130 final double[] adjoint = new double[] {1, 2, 3, 4, 5, 6, 7};
131 final SpacecraftState state = buildState(adjoint, mass);
132 final AdditionalDerivativesProvider derivativesProvider = new CartesianAdjointDerivativesProvider(penalizedCartesianFuel);
133
134 final CombinedDerivatives actualDerivatives = derivativesProvider.combinedDerivatives(state);
135
136 final BoundedCartesianEnergy energy = new BoundedCartesianEnergy(ADJOINT_NAME, massFlowRateFactor, maximumThrustMagnitude);
137 final CombinedDerivatives expectedDerivatives = new CartesianAdjointDerivativesProvider(energy).combinedDerivatives(state);
138 Assertions.assertArrayEquals(expectedDerivatives.getAdditionalDerivatives(),
139 actualDerivatives.getAdditionalDerivatives(), 1e-15);
140 Assertions.assertArrayEquals(expectedDerivatives.getMainStateDerivativesIncrements(),
141 actualDerivatives.getMainStateDerivativesIncrements(), 1e-18);
142 }
143
144 private static SpacecraftState buildState(final double[] adjoint, final double mass) {
145 final CartesianOrbit orbit = new CartesianOrbit(new PVCoordinates(Vector3D.MINUS_I, Vector3D.MINUS_K),
146 FramesFactory.getEME2000(), AbsoluteDate.ARBITRARY_EPOCH, 1.);
147 return new SpacecraftState(orbit, mass).addAdditionalData(ADJOINT_NAME, adjoint);
148 }
149 }