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