1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.estimation.leastsquares;
18
19 import org.hipparchus.exception.LocalizedCoreFormats;
20 import org.hipparchus.geometry.euclidean.threed.Vector3D;
21 import org.hipparchus.linear.RealMatrix;
22 import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation;
23 import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
24 import org.junit.jupiter.api.Assertions;
25 import org.junit.jupiter.api.Test;
26 import org.orekit.attitudes.LofOffset;
27 import org.orekit.errors.OrekitException;
28 import org.orekit.errors.OrekitMessages;
29 import org.orekit.estimation.DSSTContext;
30 import org.orekit.estimation.DSSTEstimationTestUtils;
31 import org.orekit.estimation.measurements.EstimationsProvider;
32 import org.orekit.estimation.measurements.GroundStation;
33 import org.orekit.estimation.measurements.ObservedMeasurement;
34 import org.orekit.estimation.measurements.PVMeasurementCreator;
35 import org.orekit.estimation.measurements.Range;
36 import org.orekit.estimation.measurements.TwoWayRangeMeasurementCreator;
37 import org.orekit.estimation.measurements.RangeRateMeasurementCreator;
38 import org.orekit.estimation.measurements.modifiers.PhaseCentersRangeModifier;
39 import org.orekit.frames.LOFType;
40 import org.orekit.gnss.antenna.FrequencyPattern;
41 import org.orekit.orbits.Orbit;
42 import org.orekit.propagation.Propagator;
43 import org.orekit.propagation.conversion.DSSTPropagatorBuilder;
44 import org.orekit.propagation.semianalytical.dsst.DSSTPropagator;
45 import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
46 import org.orekit.time.AbsoluteDate;
47 import org.orekit.utils.ParameterDriver;
48 import org.orekit.utils.ParameterDriversList;
49
50 import java.util.ArrayList;
51 import java.util.List;
52
53 public class DSSTBatchLSEstimatorTest {
54
55
56
57
58 @Test
59 public void testKeplerPV() {
60
61 DSSTContext context = DSSTEstimationTestUtils.eccentricContext("regular-data:potential:tides");
62
63 final DSSTPropagatorBuilder propagatorBuilder =
64 context.createBuilder(true, 60.0, 600.0, 1.0);
65
66
67 final Propagator propagator = DSSTEstimationTestUtils.createPropagator(context.initialOrbit,
68 propagatorBuilder);
69 final List<ObservedMeasurement<?>> measurements =
70 DSSTEstimationTestUtils.createMeasurements(propagator,
71 new PVMeasurementCreator(),
72 0.0, 1.0, 300.0);
73
74
75 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
76 propagatorBuilder);
77 for (final ObservedMeasurement<?> measurement : measurements) {
78 estimator.addMeasurement(measurement);
79 }
80 estimator.setParametersConvergenceThreshold(1.0e-2);
81 estimator.setMaxIterations(10);
82 estimator.setMaxEvaluations(20);
83
84 DSSTEstimationTestUtils.checkFit(context, estimator, 1, 2,
85 0.0, 4.8e-9,
86 0.0, 2.6e-8,
87 0.0, 8.9e-9,
88 0.0, 4.4e-12);
89
90 RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
91 RealMatrix physicalCovariances = estimator.getPhysicalCovariances(1.0e-10);
92 Assertions.assertEquals(6, normalizedCovariances.getRowDimension());
93 Assertions.assertEquals(6, normalizedCovariances.getColumnDimension());
94 Assertions.assertEquals(6, physicalCovariances.getRowDimension());
95 Assertions.assertEquals(6, physicalCovariances.getColumnDimension());
96 Assertions.assertEquals(0.00258, physicalCovariances.getEntry(0, 0), 1.0e-5);
97
98 }
99
100
101 @Test
102 public void testKeplerPVBackward() {
103
104 DSSTContext context = DSSTEstimationTestUtils.eccentricContext("regular-data:potential:tides");
105
106 final DSSTPropagatorBuilder propagatorBuilder =
107 context.createBuilder(true, 60.0, 600.0, 1.0);
108
109
110 final Propagator propagator = DSSTEstimationTestUtils.createPropagator(context.initialOrbit,
111 propagatorBuilder);
112 final List<ObservedMeasurement<?>> measurements =
113 DSSTEstimationTestUtils.createMeasurements(propagator,
114 new PVMeasurementCreator(),
115 0.0, -1.0, 300.0);
116
117
118 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
119 propagatorBuilder);
120 for (final ObservedMeasurement<?> measurement : measurements) {
121 estimator.addMeasurement(measurement);
122 }
123 estimator.setParametersConvergenceThreshold(1.0e-2);
124 estimator.setMaxIterations(10);
125 estimator.setMaxEvaluations(20);
126
127 DSSTEstimationTestUtils.checkFit(context, estimator, 1, 3,
128 0.0, 4.8e-9,
129 0.0, 2.7e-8,
130 0.0, 3.9e-9,
131 0.0, 1.9e-12);
132
133 RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
134 RealMatrix physicalCovariances = estimator.getPhysicalCovariances(1.0e-10);
135 Assertions.assertEquals(6, normalizedCovariances.getRowDimension());
136 Assertions.assertEquals(6, normalizedCovariances.getColumnDimension());
137 Assertions.assertEquals(6, physicalCovariances.getRowDimension());
138 Assertions.assertEquals(6, physicalCovariances.getColumnDimension());
139 Assertions.assertEquals(0.00258, physicalCovariances.getEntry(0, 0), 1.0e-5);
140
141 }
142
143
144
145
146 @Test
147 public void testKeplerRange() {
148
149 DSSTContext context = DSSTEstimationTestUtils.eccentricContext("regular-data:potential:tides");
150
151 final DSSTPropagatorBuilder propagatorBuilder =
152 context.createBuilder(true, 60.0, 600.0, 1.0);
153
154
155 final Propagator propagator = DSSTEstimationTestUtils.createPropagator(context.initialOrbit,
156 propagatorBuilder);
157 final List<ObservedMeasurement<?>> measurements =
158 DSSTEstimationTestUtils.createMeasurements(propagator,
159 new TwoWayRangeMeasurementCreator(context),
160 1.0, 3.0, 300.0);
161
162
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(10);
170 estimator.setMaxEvaluations(20);
171 estimator.setObserver(new BatchLSObserver() {
172 int lastIter = 0;
173 int lastEval = 0;
174
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 Assertions.assertEquals(lastEval + 1, evaluationscount);
184 } else {
185 Assertions.assertEquals(lastIter + 1, iterationsCount);
186 }
187 lastIter = iterationsCount;
188 lastEval = evaluationscount;
189 Assertions.assertEquals(measurements.size(), evaluationsProvider.getNumber());
190 try {
191 evaluationsProvider.getEstimatedMeasurement(-1);
192 Assertions.fail("an exception should have been thrown");
193 } catch (OrekitException oe) {
194 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
195 }
196 try {
197 evaluationsProvider.getEstimatedMeasurement(measurements.size());
198 Assertions.fail("an exception should have been thrown");
199 } catch (OrekitException oe) {
200 Assertions.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 Assertions.assertTrue(current.compareTo(previous) >= 0);
206 previous = current;
207 }
208 }
209 });
210
211 ParameterDriver aDriver = estimator.getOrbitalParametersDrivers(true).getDrivers().get(0);
212 Assertions.assertEquals("a", aDriver.getName());
213 aDriver.setValue(aDriver.getValue() + 1.2);
214 aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
215
216 DSSTEstimationTestUtils.checkFit(context, estimator, 2, 3,
217 0.0, 3.1e-6,
218 0.0, 5.7e-6,
219 0.0, 1.5e-6,
220 0.0, 6.1e-10);
221
222
223
224 for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
225 if ("a".equals(driver.getName())) {
226
227 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
228 } else {
229
230 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(propagatorBuilder.getInitialOrbitDate()), 1.0e-15);
231 }
232 }
233
234 }
235
236
237
238
239 @Test
240 public void testKeplerRangeWithOnBoardAntennaOffset() {
241
242 DSSTContext context = DSSTEstimationTestUtils.eccentricContext("regular-data:potential:tides");
243
244 final DSSTPropagatorBuilder propagatorBuilder =
245 context.createBuilder(true, 60.0, 600.0, 1.0);
246 propagatorBuilder.setAttitudeProvider(new LofOffset(propagatorBuilder.getFrame(), LOFType.LVLH));
247 final Vector3D antennaPhaseCenter = new Vector3D(-1.2, 2.3, -0.7);
248
249
250 final Propagator propagator = DSSTEstimationTestUtils.createPropagator(context.initialOrbit,
251 propagatorBuilder);
252 final List<ObservedMeasurement<?>> measurements =
253 DSSTEstimationTestUtils.createMeasurements(propagator,
254 new TwoWayRangeMeasurementCreator(context,
255 Vector3D.ZERO, null,
256 antennaPhaseCenter, null,
257 0),
258 1.0, 3.0, 300.0);
259
260
261 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
262 propagatorBuilder);
263 final PhaseCentersRangeModifier obaModifier = new PhaseCentersRangeModifier(FrequencyPattern.ZERO_CORRECTION,
264 new FrequencyPattern(antennaPhaseCenter,
265 null));
266 for (final ObservedMeasurement<?> range : measurements) {
267 ((Range) range).addModifier(obaModifier);
268 estimator.addMeasurement(range);
269 }
270 estimator.setParametersConvergenceThreshold(1.0e-2);
271 estimator.setMaxIterations(10);
272 estimator.setMaxEvaluations(20);
273 estimator.setObserver(new BatchLSObserver() {
274 int lastIter = 0;
275 int lastEval = 0;
276
277 @Override
278 public void evaluationPerformed(int iterationsCount, int evaluationscount,
279 Orbit[] orbits,
280 ParameterDriversList estimatedOrbitalParameters,
281 ParameterDriversList estimatedPropagatorParameters,
282 ParameterDriversList estimatedMeasurementsParameters,
283 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
284 if (iterationsCount == lastIter) {
285 Assertions.assertEquals(lastEval + 1, evaluationscount);
286 } else {
287 Assertions.assertEquals(lastIter + 1, iterationsCount);
288 }
289 lastIter = iterationsCount;
290 lastEval = evaluationscount;
291 Assertions.assertEquals(measurements.size(), evaluationsProvider.getNumber());
292 try {
293 evaluationsProvider.getEstimatedMeasurement(-1);
294 Assertions.fail("an exception should have been thrown");
295 } catch (OrekitException oe) {
296 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
297 }
298 try {
299 evaluationsProvider.getEstimatedMeasurement(measurements.size());
300 Assertions.fail("an exception should have been thrown");
301 } catch (OrekitException oe) {
302 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
303 }
304 AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
305 for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
306 AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
307 Assertions.assertTrue(current.compareTo(previous) >= 0);
308 previous = current;
309 }
310 }
311 });
312
313 ParameterDriver aDriver = estimator.getOrbitalParametersDrivers(true).getDrivers().get(0);
314 Assertions.assertEquals("a", aDriver.getName());
315 aDriver.setValue(aDriver.getValue() + 1.2);
316 aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
317
318 DSSTEstimationTestUtils.checkFit(context, estimator, 2, 3,
319 0.0, 2.3e-5,
320 0.0, 5.9e-5,
321 0.0, 2.7e-5,
322 0.0, 1.1e-8);
323
324
325
326 for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
327 if ("a".equals(driver.getName())) {
328
329 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
330 } else {
331
332 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(propagatorBuilder.getInitialOrbitDate()), 1.0e-15);
333 }
334 }
335
336 }
337
338 @Test
339 public void testWrappedException() {
340
341 DSSTContext context = DSSTEstimationTestUtils.eccentricContext("regular-data:potential:tides");
342
343 final DSSTPropagatorBuilder propagatorBuilder =
344 context.createBuilder(true, 60.0, 600.0, 1.0);
345
346
347 final Propagator propagator = DSSTEstimationTestUtils.createPropagator(context.initialOrbit,
348 propagatorBuilder);
349 final List<ObservedMeasurement<?>> measurements =
350 DSSTEstimationTestUtils.createMeasurements(propagator,
351 new TwoWayRangeMeasurementCreator(context),
352 1.0, 3.0, 300.0);
353
354
355 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
356 propagatorBuilder);
357 for (final ObservedMeasurement<?> range : measurements) {
358 estimator.addMeasurement(range);
359 }
360 estimator.setParametersConvergenceThreshold(1.0e-2);
361 estimator.setMaxIterations(10);
362 estimator.setMaxEvaluations(20);
363 estimator.setObserver(new BatchLSObserver() {
364
365 @Override
366 public void evaluationPerformed(int iterationsCount, int evaluationscount,
367 Orbit[] orbits,
368 ParameterDriversList estimatedOrbitalParameters,
369 ParameterDriversList estimatedPropagatorParameters,
370 ParameterDriversList estimatedMeasurementsParameters,
371 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) throws DummyException {
372 throw new DummyException();
373 }
374 });
375
376 try {
377 DSSTEstimationTestUtils.checkFit(context, estimator, 3, 4,
378 0.0, 1.5e-6,
379 0.0, 3.2e-6,
380 0.0, 3.8e-7,
381 0.0, 1.5e-10);
382 Assertions.fail("an exception should have been thrown");
383 } catch (DummyException de) {
384
385 }
386
387 }
388
389 private static class DummyException extends OrekitException {
390 private static final long serialVersionUID = 1L;
391 public DummyException() {
392 super(OrekitMessages.INTERNAL_ERROR);
393 }
394 }
395
396
397
398
399 @Test
400 public void testKeplerRangeRate() {
401
402 DSSTContext context = DSSTEstimationTestUtils.eccentricContext("regular-data:potential:tides");
403
404 final DSSTPropagatorBuilder propagatorBuilder =
405 context.createBuilder(true, 60.0, 600.0, 1.0);
406
407
408 final Propagator propagator = DSSTEstimationTestUtils.createPropagator(context.initialOrbit,
409 propagatorBuilder);
410 final double groundClockDrift = 4.8e-9;
411 for (final GroundStation station : context.stations) {
412 station.getClockDriftDriver().setValue(groundClockDrift);
413 }
414 final double satClkDrift = 3.2e-10;
415 final List<ObservedMeasurement<?>> measurements1 =
416 DSSTEstimationTestUtils.createMeasurements(propagator,
417 new RangeRateMeasurementCreator(context, false, satClkDrift),
418 1.0, 3.0, 300.0);
419
420 final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
421 measurements.addAll(measurements1);
422
423
424 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
425 propagatorBuilder);
426 for (final ObservedMeasurement<?> rangerate : measurements) {
427 estimator.addMeasurement(rangerate);
428 }
429 estimator.setParametersConvergenceThreshold(1.0e-3);
430 estimator.setMaxIterations(10);
431 estimator.setMaxEvaluations(20);
432
433 DSSTEstimationTestUtils.checkFit(context, estimator, 1, 2,
434 0.0, 5.4e-7,
435 0.0, 1.2e-6,
436 0.0, 8.3e-4,
437 0.0, 4.5e-7);
438 }
439
440
441
442
443 @Test
444 public void testKeplerRangeAndRangeRate() {
445
446 DSSTContext context = DSSTEstimationTestUtils.eccentricContext("regular-data:potential:tides");
447
448 final DSSTPropagatorBuilder propagatorBuilder =
449 context.createBuilder(true, 60.0, 600.0, 1.0);
450
451
452 final Propagator propagator = DSSTEstimationTestUtils.createPropagator(context.initialOrbit,
453 propagatorBuilder);
454
455 final List<ObservedMeasurement<?>> measurementsRange =
456 DSSTEstimationTestUtils.createMeasurements(propagator,
457 new TwoWayRangeMeasurementCreator(context),
458 1.0, 3.0, 300.0);
459 final double groundClockDrift = 4.8e-9;
460 for (final GroundStation station : context.stations) {
461 station.getClockDriftDriver().setValue(groundClockDrift);
462 }
463 final double satClkDrift = 3.2e-10;
464 final List<ObservedMeasurement<?>> measurementsRangeRate =
465 DSSTEstimationTestUtils.createMeasurements(propagator,
466 new RangeRateMeasurementCreator(context, false, satClkDrift),
467 1.0, 3.0, 300.0);
468
469
470 final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
471 measurements.addAll(measurementsRange);
472 measurements.addAll(measurementsRangeRate);
473
474
475 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
476 propagatorBuilder);
477 for (final ObservedMeasurement<?> meas : measurements) {
478 estimator.addMeasurement(meas);
479 }
480 estimator.setParametersConvergenceThreshold(1.0e-3);
481 estimator.setMaxIterations(10);
482 estimator.setMaxEvaluations(20);
483
484
485 DSSTEstimationTestUtils.checkFit(context, estimator, 1, 3,
486 0.0, 5.1e-7,
487 0.0, 3e-6,
488 0.0, 9.e-8,
489 0.0, 6e-11);
490 }
491
492 @Test
493 public void testIssue359() {
494 DSSTContext context = DSSTEstimationTestUtils.eccentricContext("regular-data:potential:tides");
495
496 final DSSTPropagatorBuilder propagatorBuilder =
497 context.createBuilder(true, 60.0, 600.0, 1.0);
498
499
500
501 propagatorBuilder.getPropagationParametersDrivers().getDrivers().get(0).setSelected(true);
502
503 final DSSTPropagator propagator = (DSSTPropagator) DSSTEstimationTestUtils.createPropagator(context.initialOrbit,
504 propagatorBuilder);
505 final List<ObservedMeasurement<?>> measurements =
506 DSSTEstimationTestUtils.createMeasurements(propagator,
507 new PVMeasurementCreator(),
508 0.0, 1.0, 300.0);
509
510
511 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
512 propagatorBuilder);
513 for (final ObservedMeasurement<?> measurement : measurements) {
514 estimator.addMeasurement(measurement);
515 }
516 ParameterDriversList estimatedParameters = estimator.getPropagatorParametersDrivers(true);
517
518 final String driverName = DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT;
519 Assertions.assertTrue(propagator.getAllForceModels().get(0) instanceof DSSTNewtonianAttraction);
520 Assertions.assertTrue(propagatorBuilder.getAllForceModels().get(0) instanceof DSSTNewtonianAttraction);
521 Assertions.assertNotNull(estimatedParameters.findByName(driverName));
522 Assertions.assertTrue(propagator.getAllForceModels().get(0).getParametersDrivers().get(0).isSelected());
523 Assertions.assertTrue(propagatorBuilder.getAllForceModels().get(0).getParametersDrivers().get(0).isSelected());
524 }
525
526 }