1   /* Copyright 2002-2025 CS GROUP
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.estimation.leastsquares;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  
22  import org.hipparchus.exception.LocalizedCoreFormats;
23  import org.hipparchus.linear.RealMatrix;
24  import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation;
25  import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
26  import org.junit.jupiter.api.Assertions;
27  import org.junit.jupiter.api.Test;
28  import org.orekit.attitudes.FrameAlignedProvider;
29  import org.orekit.errors.OrekitException;
30  import org.orekit.errors.OrekitMessages;
31  import org.orekit.estimation.TLEContext;
32  import org.orekit.estimation.TLEEstimationTestUtils;
33  import org.orekit.estimation.measurements.EstimationsProvider;
34  import org.orekit.estimation.measurements.GroundStation;
35  import org.orekit.estimation.measurements.ObservedMeasurement;
36  import org.orekit.estimation.measurements.PVMeasurementCreator;
37  import org.orekit.estimation.measurements.RangeRateMeasurementCreator;
38  import org.orekit.estimation.measurements.TwoWayRangeMeasurementCreator;
39  import org.orekit.frames.FramesFactory;
40  import org.orekit.orbits.Orbit;
41  import org.orekit.orbits.OrbitType;
42  import org.orekit.propagation.Propagator;
43  import org.orekit.propagation.analytical.tle.TLEPropagator;
44  import org.orekit.propagation.conversion.TLEPropagatorBuilder;
45  import org.orekit.time.AbsoluteDate;
46  import org.orekit.utils.ParameterDriver;
47  import org.orekit.utils.ParameterDriversList;
48  
49  public class TLEBatchLSEstimatorTest {
50  
51      /**
52       * Perfect PV measurements with a perfect start
53       */
54      @Test
55      void testPV() {
56  
57          TLEContext context = TLEEstimationTestUtils.eccentricContext("regular-data:potential:tides");
58  
59          final TLEPropagatorBuilder propagatorBuilder =
60                          context.createBuilder(1.0);
61  
62          // create perfect PV measurements
63          final Orbit initialOrbit = TLEPropagator.selectExtrapolator(context.initialTLE).getInitialState().getOrbit();
64          final Propagator propagator = TLEEstimationTestUtils.createPropagator(initialOrbit,
65                                                                             propagatorBuilder);
66          final List<ObservedMeasurement<?>> measurements =
67                          TLEEstimationTestUtils.createMeasurements(propagator,
68                                                                 new PVMeasurementCreator(),
69                                                                 0.0, 1.0, 300.0);
70  
71          // create orbit estimator
72          final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
73                                                                  propagatorBuilder);
74          for (final ObservedMeasurement<?> measurement : measurements) {
75              estimator.addMeasurement(measurement);
76          }
77          estimator.setParametersConvergenceThreshold(1.0e-2);
78          estimator.setMaxIterations(10);
79          estimator.setMaxEvaluations(20);
80  
81          TLEEstimationTestUtils.checkFit(context, estimator, 1, 1,
82                                          0.0, 1.0e-15,
83                                          0.0, 1.0e-15,
84                                          0.0, 4.97e-6,
85                                          0.0, 2.32e-9);
86  
87          RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
88          RealMatrix physicalCovariances   = estimator.getPhysicalCovariances(1.0e-10);
89          Assertions.assertEquals(6,       normalizedCovariances.getRowDimension());
90          Assertions.assertEquals(6,       normalizedCovariances.getColumnDimension());
91          Assertions.assertEquals(6,       physicalCovariances.getRowDimension());
92          Assertions.assertEquals(6,       physicalCovariances.getColumnDimension());
93          Assertions.assertEquals(0.03071, physicalCovariances.getEntry(0, 0), 1.0e-5);
94  
95      }
96  
97      /** Test PV measurements generation and backward propagation in least-square orbit determination. */
98      @Test
99      void testPVBackward() {
100 
101         TLEContext context = TLEEstimationTestUtils.eccentricContext("regular-data:potential:tides");
102 
103         final TLEPropagatorBuilder propagatorBuilder =
104                         context.createBuilder(1.0);
105 
106         // create perfect PV measurements
107         final Orbit initialOrbit = TLEPropagator.selectExtrapolator(context.initialTLE).getInitialState().getOrbit();
108         final Propagator propagator = TLEEstimationTestUtils.createPropagator(initialOrbit,
109                                                                            propagatorBuilder);
110         final List<ObservedMeasurement<?>> measurements =
111                         TLEEstimationTestUtils.createMeasurements(propagator,
112                                                                new PVMeasurementCreator(),
113                                                                0.0, -1.0, 300.0);
114 
115         // create orbit estimator
116         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
117                                                                 propagatorBuilder);
118         for (final ObservedMeasurement<?> measurement : measurements) {
119             estimator.addMeasurement(measurement);
120         }
121         estimator.setParametersConvergenceThreshold(1.0e-2);
122         estimator.setMaxIterations(10);
123         estimator.setMaxEvaluations(20);
124 
125         TLEEstimationTestUtils.checkFit(context, estimator, 1, 1,
126                                      0.0, 1.0e-15,
127                                      0.0, 1.0e-15,
128                                      0.0, 4.97e-6,
129                                      0.0, 2.32e-9);
130 
131         RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
132         RealMatrix physicalCovariances   = estimator.getPhysicalCovariances(1.0e-10);
133         Assertions.assertEquals(6,       normalizedCovariances.getRowDimension());
134         Assertions.assertEquals(6,       normalizedCovariances.getColumnDimension());
135         Assertions.assertEquals(6,       physicalCovariances.getRowDimension());
136         Assertions.assertEquals(6,       physicalCovariances.getColumnDimension());
137         Assertions.assertEquals(0.03420, physicalCovariances.getEntry(0, 0), 1.0e-5);
138 
139     }
140 
141     /**
142      * Perfect range measurements with a biased start
143      */
144     @Test
145     void testRange() {
146 
147         TLEContext context = TLEEstimationTestUtils.eccentricContext("regular-data:potential:tides");
148 
149         final TLEPropagatorBuilder propagatorBuilder =
150                         context.createBuilder(1.0);
151         // this test based on range measurements seems to have an attitude dependence?
152         propagatorBuilder.setAttitudeProvider(new FrameAlignedProvider(FramesFactory.getEME2000()));
153 
154         // create perfect range measurements
155         final Orbit initialOrbit = TLEPropagator.selectExtrapolator(context.initialTLE).getInitialState().getOrbit();
156         final Propagator propagator = TLEEstimationTestUtils.createPropagator(initialOrbit,
157                                                                            propagatorBuilder);
158         final List<ObservedMeasurement<?>> measurements =
159                         TLEEstimationTestUtils.createMeasurements(propagator,
160                                                                new TwoWayRangeMeasurementCreator(context),
161                                                                1.0, 3.0, 300.0);
162 
163         // create orbit estimator
164         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
165                                                                 propagatorBuilder);
166         for (final ObservedMeasurement<?> range : measurements) {
167             estimator.addMeasurement(range);
168         }
169         estimator.setParametersConvergenceThreshold(1.0e-2);
170         estimator.setMaxIterations(30);
171         estimator.setMaxEvaluations(30);
172         estimator.setObserver(new BatchLSObserver() {
173             int lastIter = 0;
174             int lastEval = 0;
175             /** {@inheritDoc} */
176             @Override
177             public void evaluationPerformed(int iterationsCount, int evaluationscount,
178                                             Orbit[] orbits,
179                                             ParameterDriversList estimatedOrbitalParameters,
180                                             ParameterDriversList estimatedPropagatorParameters,
181                                             ParameterDriversList estimatedMeasurementsParameters,
182                                             EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
183                 if (iterationsCount == lastIter) {
184                     Assertions.assertEquals(lastEval + 1, evaluationscount);
185                 } else {
186                     Assertions.assertEquals(lastIter + 1, iterationsCount);
187                 }
188                 lastIter = iterationsCount;
189                 lastEval = evaluationscount;
190                 Assertions.assertEquals(measurements.size(), evaluationsProvider.getNumber());
191                 try {
192                     evaluationsProvider.getEstimatedMeasurement(-1);
193                     Assertions.fail("an exception should have been thrown");
194                 } catch (OrekitException oe) {
195                     Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
196                 }
197                 try {
198                     evaluationsProvider.getEstimatedMeasurement(measurements.size());
199                     Assertions.fail("an exception should have been thrown");
200                 } catch (OrekitException oe) {
201                     Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
202                 }
203                 AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
204                 for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
205                     AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
206                     Assertions.assertTrue(current.compareTo(previous) >= 0);
207                     previous = current;
208                 }
209             }
210         });
211 
212         ParameterDriver xDriver = estimator.getOrbitalParametersDrivers(true).getDrivers().get(0);
213         Assertions.assertEquals(OrbitType.POS_X, xDriver.getName());
214         xDriver.setValue(xDriver.getValue() + 10.0);
215         xDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
216 
217         TLEEstimationTestUtils.checkFit(context, estimator, 2, 3,
218                                         0.0, 1.67e-5,
219                                         0.0, 5.11e-5,
220                                         0.0, 1.77e-5,
221                                         0.0, 2.96e-9);
222 
223         // after the call to estimate, the parameters lacking a user-specified reference date
224         // got a default one
225         for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
226             if (OrbitType.POS_X.equals(driver.getName())) {
227                 // user-specified reference date
228                 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
229             } else {
230                 // default reference date
231                 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(propagatorBuilder.getInitialOrbitDate()), 1.0e-15);
232             }
233         }
234     }
235 
236     @Test
237     void testWrappedException() {
238 
239         TLEContext context = TLEEstimationTestUtils.eccentricContext("regular-data:potential:tides");
240 
241         final TLEPropagatorBuilder propagatorBuilder =
242                         context.createBuilder(1.0);
243 
244         // create perfect range measurements
245         final Orbit initialOrbit = TLEPropagator.selectExtrapolator(context.initialTLE).getInitialState().getOrbit();
246         final Propagator propagator = TLEEstimationTestUtils.createPropagator(initialOrbit,
247                                                                                propagatorBuilder);
248         final List<ObservedMeasurement<?>> measurements =
249                         TLEEstimationTestUtils.createMeasurements(propagator,
250                                                                new TwoWayRangeMeasurementCreator(context),
251                                                                1.0, 3.0, 300.0);
252 
253         // create orbit estimator
254         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
255                                                                 propagatorBuilder);
256         for (final ObservedMeasurement<?> range : measurements) {
257             estimator.addMeasurement(range);
258         }
259         estimator.setParametersConvergenceThreshold(1.0e-2);
260         estimator.setMaxIterations(10);
261         estimator.setMaxEvaluations(50);
262         estimator.setObserver(new BatchLSObserver() {
263             /** {@inheritDoc} */
264             @Override
265             public void evaluationPerformed(int iterationsCount, int evaluationscount,
266                                            Orbit[] orbits,
267                                            ParameterDriversList estimatedOrbitalParameters,
268                                            ParameterDriversList estimatedPropagatorParameters,
269                                            ParameterDriversList estimatedMeasurementsParameters,
270                                            EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) throws DummyException {
271                 throw new DummyException();
272             }
273         });
274 
275         try {
276             TLEEstimationTestUtils.checkFit(context, estimator, 3, 4,
277                                          0.0, 1.5e-6,
278                                          0.0, 3.2e-6,
279                                          0.0, 3.8e-7,
280                                          0.0, 1.5e-10);
281             Assertions.fail("an exception should have been thrown");
282         } catch (DummyException de) {
283             // expected
284         }
285 
286     }
287 
288     private static class DummyException extends OrekitException {
289         private static final long serialVersionUID = 1L;
290         public DummyException() {
291             super(OrekitMessages.INTERNAL_ERROR);
292         }
293     }
294 
295     /**
296      * Perfect range and range rate measurements with a perfect start
297      */
298     @Test
299     void testRangeAndRangeRate() {
300 
301         TLEContext context = TLEEstimationTestUtils.eccentricContext("regular-data:potential:tides");
302 
303         final TLEPropagatorBuilder propagatorBuilder =
304                         context.createBuilder(1.0);
305 
306         // create perfect range measurements
307         final Orbit initialOrbit = TLEPropagator.selectExtrapolator(context.initialTLE).getInitialState().getOrbit();
308         final Propagator propagator = TLEEstimationTestUtils.createPropagator(initialOrbit,
309                                                                            propagatorBuilder);
310 
311         final List<ObservedMeasurement<?>> measurementsRange =
312                         TLEEstimationTestUtils.createMeasurements(propagator,
313                                                                new TwoWayRangeMeasurementCreator(context),
314                                                                1.0, 3.0, 300.0);
315         final double groundClockDrift =  4.8e-9;
316         for (final GroundStation station : context.stations) {
317             station.getClockDriftDriver().setValue(groundClockDrift);
318         }
319         final double satClkDrift = 3.2e-10;
320         final List<ObservedMeasurement<?>> measurementsRangeRate =
321                         TLEEstimationTestUtils.createMeasurements(propagator,
322                                                                new RangeRateMeasurementCreator(context, false, satClkDrift),
323                                                                1.0, 3.0, 300.0);
324 
325         // concat measurements
326         final List<ObservedMeasurement<?>> measurements = new ArrayList<>();
327         measurements.addAll(measurementsRange);
328         measurements.addAll(measurementsRangeRate);
329 
330         // create orbit estimator
331         final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
332                                                                 propagatorBuilder);
333         for (final ObservedMeasurement<?> meas : measurements) {
334             estimator.addMeasurement(meas);
335         }
336         estimator.setParametersConvergenceThreshold(1.0e-3);
337         estimator.setMaxIterations(10);
338         estimator.setMaxEvaluations(20);
339 
340         // we have low correlation between the two types of measurement. We can expect a good estimate.
341         TLEEstimationTestUtils.checkFit(context, estimator, 4, 5,
342                                         0.0, 6.0e-6,
343                                         0.0, 4.6e-5,
344                                         0.0, 6.1e-6,
345                                         0.0, 4.1e-9);
346     }
347 
348 }