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.GaussNewtonOptimizer;
23 import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer;
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.TestUtils;
29 import org.orekit.attitudes.LofOffset;
30 import org.orekit.errors.OrekitException;
31 import org.orekit.errors.OrekitMessages;
32 import org.orekit.estimation.Context;
33 import org.orekit.estimation.EstimationTestUtils;
34 import org.orekit.estimation.Force;
35 import org.orekit.estimation.measurements.AngularAzElMeasurementCreator;
36 import org.orekit.estimation.measurements.EstimationsProvider;
37 import org.orekit.estimation.measurements.GroundStation;
38 import org.orekit.estimation.measurements.InterSatellitesRangeMeasurementCreator;
39 import org.orekit.estimation.measurements.MultiplexedMeasurement;
40 import org.orekit.estimation.measurements.ObservableSatellite;
41 import org.orekit.estimation.measurements.ObservedMeasurement;
42 import org.orekit.estimation.measurements.PV;
43 import org.orekit.estimation.measurements.PVMeasurementCreator;
44 import org.orekit.estimation.measurements.Range;
45 import org.orekit.estimation.measurements.TwoWayRangeMeasurementCreator;
46 import org.orekit.estimation.measurements.RangeRateMeasurementCreator;
47 import org.orekit.estimation.measurements.modifiers.OutlierFilter;
48 import org.orekit.estimation.measurements.modifiers.PhaseCentersRangeModifier;
49 import org.orekit.forces.gravity.NewtonianAttraction;
50 import org.orekit.forces.radiation.RadiationSensitive;
51 import org.orekit.frames.LOFType;
52 import org.orekit.gnss.antenna.FrequencyPattern;
53 import org.orekit.orbits.CartesianOrbit;
54 import org.orekit.orbits.KeplerianOrbit;
55 import org.orekit.orbits.Orbit;
56 import org.orekit.orbits.OrbitType;
57 import org.orekit.orbits.PositionAngleType;
58 import org.orekit.propagation.BoundedPropagator;
59 import org.orekit.propagation.EphemerisGenerator;
60 import org.orekit.propagation.Propagator;
61 import org.orekit.propagation.conversion.KeplerianPropagatorBuilder;
62 import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
63 import org.orekit.propagation.numerical.NumericalPropagator;
64 import org.orekit.time.AbsoluteDate;
65 import org.orekit.time.ChronologicalComparator;
66 import org.orekit.utils.PVCoordinates;
67 import org.orekit.utils.ParameterDriver;
68 import org.orekit.utils.ParameterDriversList;
69 import org.orekit.utils.ParameterDriversList.DelegatingDriver;
70 import org.orekit.utils.TimeStampedPVCoordinates;
71
72 import java.util.ArrayList;
73 import java.util.List;
74 import java.util.stream.IntStream;
75
76 class BatchLSEstimatorTest {
77
78 @Test
79 void testIssue1864() {
80
81 EstimationTestUtils.eccentricContext("regular-data:potential:tides");
82 final Orbit orbit = TestUtils.getDefaultOrbit(AbsoluteDate.ARBITRARY_EPOCH);
83 final KeplerianPropagatorBuilder propagatorBuilder = new KeplerianPropagatorBuilder(OrbitType.CARTESIAN.convertType(orbit),
84 PositionAngleType.TRUE, 1.);
85 final BatchLSEstimator estimator = new BatchLSEstimator(new GaussNewtonOptimizer(), propagatorBuilder);
86 estimator.setParametersConvergenceThreshold(1e-2);
87 estimator.setMaxIterations(100);
88 estimator.setMaxEvaluations(100);
89 final ObservableSatellite satellite = new ObservableSatellite(0);
90 for (int i = 0; i < 10; i++) {
91 final TimeStampedPVCoordinates shifted = orbit.shiftedBy(i * 100).getPVCoordinates();
92 estimator.addMeasurement(new PV(shifted.getDate(), shifted.getPosition(), shifted.getVelocity(),
93 1e1, 1.e-1, 1., satellite));
94 }
95 final PVCoordinates pvCoordinates = orbit.shiftedBy(1e4).getPVCoordinates();
96 final PV outlier = new PV(orbit.getDate().shiftedBy(-1), pvCoordinates.getPosition(), pvCoordinates.getVelocity(),
97 1e1, 1.e-1, 1., satellite);
98 final OutlierFilter<PV> outlierFilter = new OutlierFilter<>(1, 1.);
99 outlier.addModifier(outlierFilter);
100 estimator.addMeasurement(outlier);
101
102 estimator.estimate();
103 final LeastSquaresOptimizer.Optimum optimum = estimator.getOptimum();
104
105 Assertions.assertEquals(6, optimum.getIterations());
106 Assertions.assertEquals(9.3177e-19, optimum.getChiSquare(), 1e-22);
107 final double actualReduced = optimum.getReducedChiSquare(0);
108 Assertions.assertEquals(1.5275e-20, actualReduced, 1e-22);
109 }
110
111
112
113
114 @Test
115 void testKeplerPVMultipleDrag() {
116
117 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
118
119 final NumericalPropagatorBuilder propagatorBuilder =
120 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
121 1.0e-6, 60.0, 1.0, Force.DRAG);
122
123 for (ParameterDriver driver:propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
124 if (driver.getName().equals("drag coefficient")) {
125 driver.setSelected(true);
126 driver.addSpanAtDate(context.initialOrbit.getDate());
127 }
128 }
129
130
131 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
132 propagatorBuilder);
133 final List<ObservedMeasurement<?>> measurements =
134 EstimationTestUtils.createMeasurements(propagator,
135 new PVMeasurementCreator(),
136 -3.0, 3.0, 300.0);
137
138
139
140
141 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
142 propagatorBuilder);
143 measurements.forEach(estimator::addMeasurement);
144 estimator.setParametersConvergenceThreshold(1.0e-2);
145 estimator.setMaxIterations(10);
146 estimator.setMaxEvaluations(20);
147
148 EstimationTestUtils.checkFit(context, estimator, 1, 2,
149 0.0, 7.8e-8,
150 0.0, 6.0e-7,
151 0.0, 3.2e-7,
152 0.0, 1.3e-10);
153
154
155 List<DelegatingDriver> Orbparameters = estimator.getOrbitalParametersDrivers(true).getDrivers();
156 Assertions.assertEquals(context.initialOrbit.getA(), Orbparameters.get(0).getValue() , 1.0e-8);
157 Assertions.assertEquals(context.initialOrbit.getE(), Orbparameters.get(1).getValue() , 1.0e-12);
158 Assertions.assertEquals(context.initialOrbit.getI(), Orbparameters.get(2).getValue() , 1.0e-12);
159
160 RealMatrix jacobian = estimator.getOptimum().getJacobian();
161 Assertions.assertEquals(8, jacobian.getColumnDimension());
162
163 }
164
165
166
167
168 @Test
169 void testKeplerPV() {
170
171 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
172
173 final NumericalPropagatorBuilder propagatorBuilder =
174 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
175 1.0e-6, 60.0, 1.0);
176
177
178 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
179 propagatorBuilder);
180 final List<ObservedMeasurement<?>> measurements =
181 EstimationTestUtils.createMeasurements(propagator,
182 new PVMeasurementCreator(),
183 0.0, 1.0, 300.0);
184
185
186 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
187 propagatorBuilder);
188 measurements.forEach(estimator::addMeasurement);
189 estimator.setParametersConvergenceThreshold(1.0e-2);
190 estimator.setMaxIterations(10);
191 estimator.setMaxEvaluations(20);
192
193 EstimationTestUtils.checkFit(context, estimator, 1, 4,
194 0.0, 2.2e-8,
195 0.0, 1.1e-7,
196 0.0, 1.4e-8,
197 0.0, 6.3e-12);
198
199 RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
200 RealMatrix physicalCovariances = estimator.getPhysicalCovariances(1.0e-10);
201 Assertions.assertEquals(6, normalizedCovariances.getRowDimension());
202 Assertions.assertEquals(6, normalizedCovariances.getColumnDimension());
203 Assertions.assertEquals(6, physicalCovariances.getRowDimension());
204 Assertions.assertEquals(6, physicalCovariances.getColumnDimension());
205 Assertions.assertEquals(0.00258, physicalCovariances.getEntry(0, 0), 1.0e-5);
206
207 }
208
209
210 @Test
211 void testKeplerPVBackward() {
212
213 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
214
215 final NumericalPropagatorBuilder propagatorBuilder =
216 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
217 1.0e-6, 60.0, 1.0);
218
219
220 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
221 propagatorBuilder);
222 final List<ObservedMeasurement<?>> measurements =
223 EstimationTestUtils.createMeasurements(propagator,
224 new PVMeasurementCreator(),
225 0.0, -1.0, 300.0);
226
227
228 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
229 propagatorBuilder);
230 measurements.forEach(estimator::addMeasurement);
231 estimator.setParametersConvergenceThreshold(1.0e-2);
232 estimator.setMaxIterations(10);
233 estimator.setMaxEvaluations(20);
234
235 EstimationTestUtils.checkFit(context, estimator, 1, 2,
236 0.0, 8.3e-9,
237 0.0, 5.3e-8,
238 0.0, 5.6e-9,
239 0.0, 1.6e-12);
240
241 RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
242 RealMatrix physicalCovariances = estimator.getPhysicalCovariances(1.0e-10);
243 Assertions.assertEquals(6, normalizedCovariances.getRowDimension());
244 Assertions.assertEquals(6, normalizedCovariances.getColumnDimension());
245 Assertions.assertEquals(6, physicalCovariances.getRowDimension());
246 Assertions.assertEquals(6, physicalCovariances.getColumnDimension());
247 Assertions.assertEquals(0.00258, physicalCovariances.getEntry(0, 0), 1.0e-5);
248
249 }
250
251
252
253
254 @Test
255 void testKeplerRange() {
256
257 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
258
259 final NumericalPropagatorBuilder propagatorBuilder =
260 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
261 1.0e-6, 60.0, 1.0);
262
263
264 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
265 propagatorBuilder);
266 final List<ObservedMeasurement<?>> measurements =
267 EstimationTestUtils.createMeasurements(propagator,
268 new TwoWayRangeMeasurementCreator(context),
269 1.0, 3.0, 300.0);
270
271
272 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
273 propagatorBuilder);
274 measurements.forEach(estimator::addMeasurement);
275 estimator.setParametersConvergenceThreshold(1.0e-2);
276 estimator.setMaxIterations(10);
277 estimator.setMaxEvaluations(20);
278 estimator.setObserver(new BatchLSObserver() {
279 int lastIter = 0;
280 int lastEval = 0;
281
282 @Override
283 public void evaluationPerformed(int iterationsCount, int evaluationscount,
284 Orbit[] orbits,
285 ParameterDriversList estimatedOrbitalParameters,
286 ParameterDriversList estimatedPropagatorParameters,
287 ParameterDriversList estimatedMeasurementsParameters,
288 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
289 if (iterationsCount == lastIter) {
290 Assertions.assertEquals(lastEval + 1, evaluationscount);
291 } else {
292 Assertions.assertEquals(lastIter + 1, iterationsCount);
293 }
294 lastIter = iterationsCount;
295 lastEval = evaluationscount;
296 Assertions.assertEquals(measurements.size(), evaluationsProvider.getNumber());
297 try {
298 evaluationsProvider.getEstimatedMeasurement(-1);
299 Assertions.fail("an exception should have been thrown");
300 } catch (OrekitException oe) {
301 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
302 }
303 try {
304 evaluationsProvider.getEstimatedMeasurement(measurements.size());
305 Assertions.fail("an exception should have been thrown");
306 } catch (OrekitException oe) {
307 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
308 }
309 AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
310 for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
311 AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
312 Assertions.assertTrue(current.compareTo(previous) >= 0);
313 previous = current;
314 }
315 }
316 });
317
318 ParameterDriver aDriver = estimator.getOrbitalParametersDrivers(true).getDrivers().get(0);
319 Assertions.assertEquals("a", aDriver.getName());
320 aDriver.setValue(aDriver.getValue() + 1.2);
321 aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
322
323 EstimationTestUtils.checkFit(context, estimator, 2, 3,
324 0.0, 1.2e-6,
325 0.0, 3.0e-6,
326 0.0, 7.0e-7,
327 0.0, 3e-10);
328
329
330
331 for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
332 if ("a".equals(driver.getName())) {
333
334 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
335 } else {
336
337 Assertions.assertEquals(0,
338 driver.getReferenceDate().durationFrom(propagatorBuilder.getInitialOrbitDate()), 1.0e-15);
339 }
340 }
341
342 }
343
344
345
346
347 @Test
348 void testKeplerRangeWithOnBoardAntennaOffset() {
349
350 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
351
352 final NumericalPropagatorBuilder propagatorBuilder =
353 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
354 1.0e-6, 60.0, 1.0);
355 propagatorBuilder.setAttitudeProvider(new LofOffset(propagatorBuilder.getFrame(), LOFType.LVLH));
356 final Vector3D antennaPhaseCenter = new Vector3D(-1.2, 2.3, -0.7);
357
358
359 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
360 propagatorBuilder);
361 final List<ObservedMeasurement<?>> measurements =
362 EstimationTestUtils.createMeasurements(propagator,
363 new TwoWayRangeMeasurementCreator(context,
364 Vector3D.ZERO, null,
365 antennaPhaseCenter, null,
366 0),
367 1.0, 3.0, 300.0);
368
369
370 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
371 propagatorBuilder);
372 final PhaseCentersRangeModifier obaModifier = new PhaseCentersRangeModifier(FrequencyPattern.ZERO_CORRECTION,
373 new FrequencyPattern(antennaPhaseCenter,
374 null));
375 measurements.forEach(m -> {
376 ((Range) m).addModifier(obaModifier);
377 estimator.addMeasurement(m);
378 });
379 estimator.setParametersConvergenceThreshold(1.0e-2);
380 estimator.setMaxIterations(10);
381 estimator.setMaxEvaluations(20);
382 estimator.setObserver(new BatchLSObserver() {
383 int lastIter = 0;
384 int lastEval = 0;
385
386 @Override
387 public void evaluationPerformed(int iterationsCount, int evaluationscount,
388 Orbit[] orbits,
389 ParameterDriversList estimatedOrbitalParameters,
390 ParameterDriversList estimatedPropagatorParameters,
391 ParameterDriversList estimatedMeasurementsParameters,
392 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
393 if (iterationsCount == lastIter) {
394 Assertions.assertEquals(lastEval + 1, evaluationscount);
395 } else {
396 Assertions.assertEquals(lastIter + 1, iterationsCount);
397 }
398 lastIter = iterationsCount;
399 lastEval = evaluationscount;
400 Assertions.assertEquals(measurements.size(), evaluationsProvider.getNumber());
401 try {
402 evaluationsProvider.getEstimatedMeasurement(-1);
403 Assertions.fail("an exception should have been thrown");
404 } catch (OrekitException oe) {
405 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
406 }
407 try {
408 evaluationsProvider.getEstimatedMeasurement(measurements.size());
409 Assertions.fail("an exception should have been thrown");
410 } catch (OrekitException oe) {
411 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
412 }
413 AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
414 for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
415 AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
416 Assertions.assertTrue(current.compareTo(previous) >= 0);
417 previous = current;
418 }
419 }
420 });
421
422 ParameterDriver aDriver = estimator.getOrbitalParametersDrivers(true).getDrivers().get(0);
423 Assertions.assertEquals("a", aDriver.getName());
424 aDriver.setValue(aDriver.getValue() + 1.2);
425 aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
426
427 EstimationTestUtils.checkFit(context, estimator, 2, 3,
428 0.0, 2.0e-5,
429 0.0, 5.3e-5,
430 0.0, 2.7e-5,
431 0.0, 1.1e-8);
432
433
434
435 for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
436 if ("a".equals(driver.getName())) {
437
438 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
439 } else {
440
441 Assertions.assertEquals(0,
442 driver.getReferenceDate().durationFrom(propagatorBuilder.getInitialOrbitDate()), 1.0e-15);
443 }
444 }
445
446 }
447
448 @Test
449 void testMultiSat() {
450
451 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
452
453 final NumericalPropagatorBuilder propagatorBuilder1 =
454 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
455 1.0e-6, 60.0, 1.0e-3);
456 final NumericalPropagatorBuilder propagatorBuilder2 =
457 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
458 1.0e-6, 60.0, 1.0e-3);
459
460
461 final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
462 final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
463 original.getPosition().add(new Vector3D(1000, 2000, 3000)),
464 original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
465 context.initialOrbit.getFrame(),
466 context.initialOrbit.getMu());
467 final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit,
468 propagatorBuilder2);
469 final EphemerisGenerator generator = closePropagator.getEphemerisGenerator();
470 closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
471 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
472 Propagator propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
473 propagatorBuilder1);
474
475 final double localClockOffset = 0.137e-6;
476 final double remoteClockOffset = 469.0e-6;
477 final List<ObservedMeasurement<?>> r12 =
478 EstimationTestUtils.createMeasurements(propagator1,
479 new InterSatellitesRangeMeasurementCreator(ephemeris,
480 localClockOffset,
481 remoteClockOffset),
482 1.0, 3.0, 300.0);
483
484
485 propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
486 propagatorBuilder1);
487 final List<ObservedMeasurement<?>> r1 =
488 EstimationTestUtils.createMeasurements(propagator1,
489 new TwoWayRangeMeasurementCreator(context),
490 1.0, 3.0, 300.0);
491
492
493 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
494 propagatorBuilder1,
495 propagatorBuilder2);
496 r1.forEach(estimator::addMeasurement);
497 r12.forEach(estimator::addMeasurement);
498 estimator.setParametersConvergenceThreshold(1.0e-2);
499 estimator.setMaxIterations(10);
500 estimator.setMaxEvaluations(20);
501 estimator.setObserver(new BatchLSObserver() {
502 int lastIter = 0;
503 int lastEval = 0;
504
505 @Override
506 public void evaluationPerformed(int iterationsCount, int evaluationscount,
507 Orbit[] orbits,
508 ParameterDriversList estimatedOrbitalParameters,
509 ParameterDriversList estimatedPropagatorParameters,
510 ParameterDriversList estimatedMeasurementsParameters,
511 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
512 if (iterationsCount == lastIter) {
513 Assertions.assertEquals(lastEval + 1, evaluationscount);
514 } else {
515 Assertions.assertEquals(lastIter + 1, iterationsCount);
516 }
517 lastIter = iterationsCount;
518 lastEval = evaluationscount;
519 Assertions.assertEquals(r12.size() + r1.size(), evaluationsProvider.getNumber());
520 try {
521 evaluationsProvider.getEstimatedMeasurement(-1);
522 Assertions.fail("an exception should have been thrown");
523 } catch (OrekitException oe) {
524 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
525 }
526 try {
527 evaluationsProvider.getEstimatedMeasurement(r12.size() + r1.size());
528 Assertions.fail("an exception should have been thrown");
529 } catch (OrekitException oe) {
530 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
531 }
532 AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
533 for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
534 AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
535 Assertions.assertTrue(current.compareTo(previous) >= 0);
536 previous = current;
537 }
538 }
539 });
540
541 List<DelegatingDriver> parameters = estimator.getOrbitalParametersDrivers(true).getDrivers();
542 ParameterDriver a0Driver = parameters.get(0);
543 Assertions.assertEquals("a[0]", a0Driver.getName());
544 a0Driver.setValue(a0Driver.getValue() + 1.2);
545 a0Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
546
547 ParameterDriver a1Driver = parameters.get(6);
548 Assertions.assertEquals("a[1]", a1Driver.getName());
549 a1Driver.setValue(a1Driver.getValue() - 5.4);
550 a1Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
551
552 final Orbit before = new KeplerianOrbit(parameters.get( 6).getValue(),
553 parameters.get( 7).getValue(),
554 parameters.get( 8).getValue(),
555 parameters.get( 9).getValue(),
556 parameters.get(10).getValue(),
557 parameters.get(11).getValue(),
558 PositionAngleType.TRUE,
559 closeOrbit.getFrame(),
560 closeOrbit.getDate(),
561 closeOrbit.getMu());
562 Assertions.assertEquals(4.7246, Vector3D.distance(closeOrbit.getPosition(),
563 before.getPosition()), 1.0e-3);
564 Assertions.assertEquals(0.0010514, Vector3D.distance(closeOrbit.getVelocity(),
565 before.getVelocity()), 1.0e-6);
566 EstimationTestUtils.checkFit(context, estimator, 3, 4,
567 0.0, 5e-06,
568 0.0, 1.3e-05,
569 0.0, 8.3e-07,
570 0.0, 3.7e-10);
571
572 final Orbit determined = new KeplerianOrbit(parameters.get( 6).getValue(),
573 parameters.get( 7).getValue(),
574 parameters.get( 8).getValue(),
575 parameters.get( 9).getValue(),
576 parameters.get(10).getValue(),
577 parameters.get(11).getValue(),
578 PositionAngleType.TRUE,
579 closeOrbit.getFrame(),
580 closeOrbit.getDate(),
581 closeOrbit.getMu());
582 Assertions.assertEquals(0.0, Vector3D.distance(closeOrbit.getPosition(),
583 determined.getPosition()), 6.2e-6);
584 Assertions.assertEquals(0.0, Vector3D.distance(closeOrbit.getVelocity(),
585 determined.getVelocity()), 1.6e-9);
586
587
588
589 for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
590 if (driver.getName().startsWith("a[")) {
591
592 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
593 } else {
594
595 Assertions.assertEquals(0,
596 driver.getReferenceDate().durationFrom(propagatorBuilder1.getInitialOrbitDate()), 1.0e-15);
597 }
598 }
599
600 }
601
602
603
604
605
606
607 @Test
608 void testMultiSatWithParameters() {
609
610
611 final boolean muEstimated = true;
612 final boolean crEstimated1 = false;
613 final boolean caEstimated1 = true;
614 final boolean crEstimated2 = true;
615 final boolean caEstimated2 = false;
616
617
618
619 final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
620 final NumericalPropagatorBuilder propagatorBuilder1 =
621 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
622 1.0e-6, 60.0, 1.0e-3, Force.POTENTIAL, Force.SOLAR_RADIATION_PRESSURE);
623
624
625 String satName = "sat 1";
626 for (DelegatingDriver driver:propagatorBuilder1.getPropagationParametersDrivers().getDrivers()) {
627 if (driver.getName().equals("central attraction coefficient")) {
628 driver.setSelected(muEstimated);
629 }
630 if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
631 driver.setName(driver.getName() + " " + satName);
632 driver.setSelected(crEstimated1);
633 }
634 if (driver.getName().equals(RadiationSensitive.ABSORPTION_COEFFICIENT)) {
635 driver.setName(driver.getName() + " " + satName);
636 driver.setSelected(caEstimated1);
637 }
638 }
639
640
641 final Context context2 = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
642 final NumericalPropagatorBuilder propagatorBuilder2 =
643 context2.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
644 1.0e-6, 60.0, 1.0e-3, Force.POTENTIAL, Force.SOLAR_RADIATION_PRESSURE);
645
646
647 satName = "sat 2";
648 for (ParameterDriver driver:propagatorBuilder2.getPropagationParametersDrivers().getDrivers()) {
649 if (driver.getName().equals("central attraction coefficient")) {
650 driver.setSelected(muEstimated);
651 }
652 if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
653 driver.setName(driver.getName() + " " + satName);
654 driver.setSelected(crEstimated2);
655 }
656 if (driver.getName().equals(RadiationSensitive.ABSORPTION_COEFFICIENT)) {
657 driver.setName(driver.getName() + " " + satName);
658 driver.setSelected(caEstimated2);
659 }
660 }
661
662
663 final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
664 final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
665 original.getPosition().add(new Vector3D(1000, 2000, 3000)),
666 original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
667 context.initialOrbit.getFrame(),
668 context.initialOrbit.getMu());
669 final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit,
670 propagatorBuilder2);
671 final EphemerisGenerator generator = closePropagator.getEphemerisGenerator();
672 closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
673 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
674 Propagator propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
675 propagatorBuilder1);
676 final List<ObservedMeasurement<?>> r12 =
677 EstimationTestUtils.createMeasurements(propagator1,
678 new InterSatellitesRangeMeasurementCreator(ephemeris, 0., 0.),
679 1.0, 3.0, 120.0);
680
681
682 propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
683 propagatorBuilder1);
684 final List<ObservedMeasurement<?>> r1 =
685 EstimationTestUtils.createMeasurements(propagator1,
686 new TwoWayRangeMeasurementCreator(context),
687 1.0, 3.0, 120.0);
688
689
690 propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
691 propagatorBuilder1);
692 final List<ObservedMeasurement<?>> a1 =
693 EstimationTestUtils.createMeasurements(propagator1,
694 new AngularAzElMeasurementCreator(context),
695 1.0, 3.0, 120.0);
696
697
698 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
699 propagatorBuilder1,
700 propagatorBuilder2);
701 r12.forEach(estimator::addMeasurement);
702 r1.forEach(estimator::addMeasurement);
703 a1.forEach(estimator::addMeasurement);
704 estimator.setParametersConvergenceThreshold(1.0e-2);
705 estimator.setMaxIterations(10);
706 estimator.setMaxEvaluations(20);
707 estimator.setObserver(new BatchLSObserver() {
708 int lastIter = 0;
709 int lastEval = 0;
710
711 @Override
712 public void evaluationPerformed(int iterationsCount, int evaluationscount,
713 Orbit[] orbits,
714 ParameterDriversList estimatedOrbitalParameters,
715 ParameterDriversList estimatedPropagatorParameters,
716 ParameterDriversList estimatedMeasurementsParameters,
717 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
718 if (iterationsCount == lastIter) {
719 Assertions.assertEquals(lastEval + 1, evaluationscount);
720 } else {
721 Assertions.assertEquals(lastIter + 1, iterationsCount);
722 }
723 lastIter = iterationsCount;
724 lastEval = evaluationscount;
725
726 AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
727 for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
728 AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
729 Assertions.assertTrue(current.compareTo(previous) >= 0);
730 previous = current;
731 }
732 }
733 });
734
735 List<DelegatingDriver> parameters = estimator.getOrbitalParametersDrivers(true).getDrivers();
736 ParameterDriver a0Driver = parameters.get(0);
737 Assertions.assertEquals("a[0]", a0Driver.getName());
738 a0Driver.setValue(a0Driver.getValue() + 1.2);
739 a0Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
740
741 ParameterDriver a1Driver = parameters.get(6);
742 Assertions.assertEquals("a[1]", a1Driver.getName());
743 a1Driver.setValue(a1Driver.getValue() - 5.4, null);
744 a1Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
745
746 final Orbit before = new KeplerianOrbit(parameters.get( 6).getValue(),
747 parameters.get( 7).getValue(),
748 parameters.get( 8).getValue(),
749 parameters.get( 9).getValue(),
750 parameters.get(10).getValue(),
751 parameters.get(11).getValue(),
752 PositionAngleType.TRUE,
753 closeOrbit.getFrame(),
754 closeOrbit.getDate(),
755 closeOrbit.getMu());
756 Assertions.assertEquals(4.7246, Vector3D.distance(closeOrbit.getPosition(),
757 before.getPosition()), 1.0e-3);
758 Assertions.assertEquals(0.0010514, Vector3D.distance(closeOrbit.getVelocity(),
759 before.getVelocity()), 1.0e-6);
760 EstimationTestUtils.checkFit(context, estimator, 4, 5,
761 0.0, 5.3e-06,
762 0.0, 1.4e-05,
763 0.0, 8.8e-07,
764 0.0, 3.6e-10);
765
766 final Orbit determined = new KeplerianOrbit(parameters.get( 6).getValue(),
767 parameters.get( 7).getValue(),
768 parameters.get( 8).getValue(),
769 parameters.get( 9).getValue(),
770 parameters.get(10).getValue(),
771 parameters.get(11).getValue(),
772 PositionAngleType.TRUE,
773 closeOrbit.getFrame(),
774 closeOrbit.getDate(),
775 closeOrbit.getMu());
776 Assertions.assertEquals(0.0, Vector3D.distance(closeOrbit.getPosition(),
777 determined.getPosition()), 5.3e-6);
778 Assertions.assertEquals(0.0, Vector3D.distance(closeOrbit.getVelocity(),
779 determined.getVelocity()), 3.2e-9);
780
781
782
783 for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
784 if (driver.getName().startsWith("a[")) {
785
786 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
787 } else {
788
789 Assertions.assertEquals(0,
790 driver.getReferenceDate().durationFrom(propagatorBuilder1.getInitialOrbitDate()), 1.0e-15);
791 }
792 }
793
794 }
795
796
797
798
799
800 @Test
801 void testIssue617() {
802
803 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
804
805 final NumericalPropagatorBuilder propagatorBuilder1 =
806 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
807 1.0e-6, 60.0, 1.0);
808 final NumericalPropagatorBuilder propagatorBuilder2 =
809 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
810 1.0e-6, 60.0, 1.0);
811
812
813 final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
814 final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
815 original.getPosition().add(new Vector3D(1000, 2000, 3000)),
816 original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
817 context.initialOrbit.getFrame(),
818 context.initialOrbit.getMu());
819 final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit,
820 propagatorBuilder2);
821 final EphemerisGenerator generator = closePropagator.getEphemerisGenerator();
822 closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
823 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
824 Propagator propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
825 propagatorBuilder1);
826
827 final List<ObservedMeasurement<?>> r12 =
828 EstimationTestUtils.createMeasurements(propagator1,
829 new InterSatellitesRangeMeasurementCreator(ephemeris, 0., 0.),
830 1.0, 3.0, 300.0);
831
832
833 propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
834 propagatorBuilder1);
835 final List<ObservedMeasurement<?>> r1 =
836 EstimationTestUtils.createMeasurements(propagator1,
837 new TwoWayRangeMeasurementCreator(context),
838 1.0, 3.0, 300.0);
839
840
841 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
842 propagatorBuilder1,
843 propagatorBuilder2);
844
845 final List<ObservedMeasurement<?>> independentMeasurements = new ArrayList<>();
846 independentMeasurements.addAll(r1);
847 independentMeasurements.addAll(r12);
848
849
850
851
852
853 final List<ObservedMeasurement<?>> multiplexed = multiplexMeasurements(independentMeasurements, 60.0);
854
855 for (final ObservedMeasurement<?> measurement : multiplexed) {
856 estimator.addMeasurement(measurement);
857 }
858
859 estimator.setParametersConvergenceThreshold(1.0e-2);
860 estimator.setMaxIterations(10);
861 estimator.setMaxEvaluations(20);
862 estimator.setObserver(new BatchLSObserver() {
863 int lastIter = 0;
864 int lastEval = 0;
865
866 @Override
867 public void evaluationPerformed(int iterationsCount, int evaluationscount,
868 Orbit[] orbits,
869 ParameterDriversList estimatedOrbitalParameters,
870 ParameterDriversList estimatedPropagatorParameters,
871 ParameterDriversList estimatedMeasurementsParameters,
872 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
873 if (iterationsCount == lastIter) {
874 Assertions.assertEquals(lastEval + 1, evaluationscount);
875 } else {
876 Assertions.assertEquals(lastIter + 1, iterationsCount);
877 }
878 lastIter = iterationsCount;
879 lastEval = evaluationscount;
880 Assertions.assertEquals(multiplexed.size(), evaluationsProvider.getNumber());
881 try {
882 evaluationsProvider.getEstimatedMeasurement(-1);
883 Assertions.fail("an exception should have been thrown");
884 } catch (OrekitException oe) {
885 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
886 }
887 try {
888 evaluationsProvider.getEstimatedMeasurement(r12.size() + r1.size());
889 Assertions.fail("an exception should have been thrown");
890 } catch (OrekitException oe) {
891 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
892 }
893 AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
894 for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
895 AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
896 Assertions.assertTrue(current.compareTo(previous) >= 0);
897 previous = current;
898 }
899 }
900 });
901
902 List<DelegatingDriver> parameters = estimator.getOrbitalParametersDrivers(true).getDrivers();
903 ParameterDriver a0Driver = parameters.get(0);
904 Assertions.assertEquals("a[0]", a0Driver.getName());
905 a0Driver.setValue(a0Driver.getValue() + 1.2);
906 a0Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
907
908 ParameterDriver a1Driver = parameters.get(6);
909 Assertions.assertEquals("a[1]", a1Driver.getName());
910 a1Driver.setValue(a1Driver.getValue() - 5.4);
911 a1Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
912
913 final Orbit before = new KeplerianOrbit(parameters.get( 6).getValue(),
914 parameters.get( 7).getValue(),
915 parameters.get( 8).getValue(),
916 parameters.get( 9).getValue(),
917 parameters.get(10).getValue(),
918 parameters.get(11).getValue(),
919 PositionAngleType.TRUE,
920 closeOrbit.getFrame(),
921 closeOrbit.getDate(),
922 closeOrbit.getMu());
923 Assertions.assertEquals(4.7246, Vector3D.distance(closeOrbit.getPosition(),
924 before.getPosition()), 1.0e-3);
925 Assertions.assertEquals(0.0010514, Vector3D.distance(closeOrbit.getVelocity(),
926 before.getVelocity()), 1.0e-6);
927 EstimationTestUtils.checkFit(context, estimator, 2, 3,
928 0.0, 2.9e-06,
929 0.0, 8.1e-06,
930 0.0, 7.1e-07,
931 0.0, 3.2e-10);
932
933 final Orbit determined = new KeplerianOrbit(parameters.get( 6).getValue(),
934 parameters.get( 7).getValue(),
935 parameters.get( 8).getValue(),
936 parameters.get( 9).getValue(),
937 parameters.get(10).getValue(),
938 parameters.get(11).getValue(),
939 PositionAngleType.TRUE,
940 closeOrbit.getFrame(),
941 closeOrbit.getDate(),
942 closeOrbit.getMu());
943 Assertions.assertEquals(0.0, Vector3D.distance(closeOrbit.getPosition(),
944 determined.getPosition()), 4.6e-6);
945 Assertions.assertEquals(0.0, Vector3D.distance(closeOrbit.getVelocity(),
946 determined.getVelocity()), 1.6e-9);
947
948
949
950 for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
951 if (driver.getName().startsWith("a[")) {
952
953 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
954 } else {
955
956 Assertions.assertEquals(0,
957 driver.getReferenceDate().durationFrom(propagatorBuilder1.getInitialOrbitDate()), 1.0e-15);
958 }
959 }
960
961 }
962
963 @Test
964 void testWrappedException() {
965
966 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
967
968 final NumericalPropagatorBuilder propagatorBuilder =
969 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
970 1.0e-6, 60.0, 1.0);
971
972
973 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
974 propagatorBuilder);
975 final List<ObservedMeasurement<?>> measurements =
976 EstimationTestUtils.createMeasurements(propagator,
977 new TwoWayRangeMeasurementCreator(context),
978 1.0, 3.0, 300.0);
979
980
981 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
982 propagatorBuilder);
983 measurements.forEach(estimator::addMeasurement);
984 estimator.setParametersConvergenceThreshold(1.0e-2);
985 estimator.setMaxIterations(10);
986 estimator.setMaxEvaluations(20);
987 estimator.setObserver(new BatchLSObserver() {
988
989 @Override
990 public void evaluationPerformed(int iterationsCount, int evaluationscount,
991 Orbit[] orbits,
992 ParameterDriversList estimatedOrbitalParameters,
993 ParameterDriversList estimatedPropagatorParameters,
994 ParameterDriversList estimatedMeasurementsParameters,
995 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
996 throw new DummyException();
997 }
998 });
999
1000 try {
1001 EstimationTestUtils.checkFit(context, estimator, 3, 4,
1002 0.0, 1.5e-6,
1003 0.0, 3.2e-6,
1004 0.0, 3.8e-7,
1005 0.0, 1.5e-10);
1006 Assertions.fail("an exception should have been thrown");
1007 } catch (DummyException de) {
1008
1009 }
1010
1011 }
1012
1013 private static class DummyException extends OrekitException {
1014 private static final long serialVersionUID = 1L;
1015 public DummyException() {
1016 super(OrekitMessages.INTERNAL_ERROR);
1017 }
1018 }
1019
1020
1021
1022
1023 @Test
1024 void testKeplerRangeRate() {
1025
1026 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
1027
1028 final NumericalPropagatorBuilder propagatorBuilder =
1029 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1030 1.0e-6, 60.0, 1.0);
1031
1032
1033 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
1034 propagatorBuilder);
1035 final double groundClockDrift = 4.8e-9;
1036 for (final GroundStation station : context.stations) {
1037 station.getClockDriftDriver().setValue(groundClockDrift);
1038 }
1039 final double satClkDrift = 3.2e-10;
1040 final List<ObservedMeasurement<?>> measurements =
1041 EstimationTestUtils.createMeasurements(propagator,
1042 new RangeRateMeasurementCreator(context, false, satClkDrift),
1043 1.0, 3.0, 300.0);
1044
1045
1046
1047 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
1048 propagatorBuilder);
1049 for (final ObservedMeasurement<?> rangerate : measurements) {
1050 estimator.addMeasurement(rangerate);
1051 }
1052 estimator.setParametersConvergenceThreshold(1.0e-3);
1053 estimator.setMaxIterations(10);
1054 estimator.setMaxEvaluations(20);
1055
1056 EstimationTestUtils.checkFit(context, estimator, 1, 2,
1057 0.0, 5.3e-7,
1058 0.0, 1.3e-6,
1059 0.0, 8.4e-4,
1060 0.0, 5.1e-7);
1061 }
1062
1063
1064
1065
1066 @Test
1067 void testKeplerRangeAndRangeRate() {
1068
1069 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
1070
1071 final NumericalPropagatorBuilder propagatorBuilder =
1072 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1073 1.0e-6, 60.0, 1.0);
1074
1075
1076 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
1077 propagatorBuilder);
1078
1079 final List<ObservedMeasurement<?>> measurementsRange =
1080 EstimationTestUtils.createMeasurements(propagator,
1081 new TwoWayRangeMeasurementCreator(context),
1082 1.0, 3.0, 300.0);
1083 final double groundClockDrift = 4.8e-9;
1084 for (final GroundStation station : context.stations) {
1085 station.getClockDriftDriver().setValue(groundClockDrift);
1086 }
1087 final double satClkDrift = 3.2e-10;
1088 final List<ObservedMeasurement<?>> measurementsRangeRate =
1089 EstimationTestUtils.createMeasurements(propagator,
1090 new RangeRateMeasurementCreator(context, false, satClkDrift),
1091 1.0, 3.0, 300.0);
1092
1093
1094 final List<ObservedMeasurement<?>> measurements = new ArrayList<>();
1095 measurements.addAll(measurementsRange);
1096 measurements.addAll(measurementsRangeRate);
1097
1098
1099 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
1100 propagatorBuilder);
1101 for (final ObservedMeasurement<?> meas : measurements) {
1102 estimator.addMeasurement(meas);
1103 }
1104 estimator.setParametersConvergenceThreshold(1.0e-3);
1105 estimator.setMaxIterations(10);
1106 estimator.setMaxEvaluations(20);
1107
1108
1109 EstimationTestUtils.checkFit(context, estimator, 1, 2,
1110 0.0, 4.6e7,
1111 0.0, 1.8e-6,
1112 0.0, 5.8e-7,
1113 0.0, 2.7e-10);
1114 }
1115
1116
1117
1118
1119 @Test
1120 void testIssue359() {
1121
1122 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
1123
1124 final NumericalPropagatorBuilder propagatorBuilder =
1125 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1126 1.0e-6, 60.0, 1.0);
1127
1128
1129
1130 propagatorBuilder.getPropagationParametersDrivers().getDrivers().get(0).setSelected(true);
1131
1132 final NumericalPropagator propagator = (NumericalPropagator) EstimationTestUtils.createPropagator(context.initialOrbit,
1133 propagatorBuilder);
1134 final List<ObservedMeasurement<?>> measurements =
1135 EstimationTestUtils.createMeasurements(propagator,
1136 new PVMeasurementCreator(),
1137 0.0, 1.0, 300.0);
1138
1139
1140 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
1141 propagatorBuilder);
1142
1143 measurements.forEach(estimator::addMeasurement);
1144
1145 ParameterDriversList estimatedParameters = estimator.getPropagatorParametersDrivers(true);
1146
1147 final String driverName = NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT;
1148 Assertions.assertInstanceOf(NewtonianAttraction.class, propagator.getAllForceModels().get(0));
1149 Assertions.assertInstanceOf(NewtonianAttraction.class, propagatorBuilder.getAllForceModels().get(0));
1150 Assertions.assertNotNull(estimatedParameters.findByName(driverName));
1151 Assertions.assertTrue(propagator.getAllForceModels().get(0).getParameterDriver(driverName).isSelected());
1152 Assertions.assertTrue(propagatorBuilder.getAllForceModels().get(0).getParameterDriver(driverName).isSelected());
1153 }
1154
1155 @Test
1156 void testEstimateOnlyOneOrbitalParameter() {
1157 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ true, false, false, false, false, false });
1158 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ false, true, false, false, false, false });
1159 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ false, false, true, false, false, false });
1160 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ false, false, false, true, false, false });
1161 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ false, false, false, false, true, false });
1162 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ false, false, false, false, false, true });
1163 }
1164
1165 @Test
1166 void testEstimateOnlyFewOrbitalParameters() {
1167 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ false, true, false, true, false, false });
1168 }
1169
1170
1171
1172
1173
1174
1175 @Test
1176 void testErrorThrownWhenNoEnabledMeasurements() {
1177
1178
1179 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
1180 final NumericalPropagatorBuilder propagatorBuilder =
1181 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1182 1.0e-6, 60.0, 1.0e-3);
1183
1184
1185 final List<DelegatingDriver> drivers = propagatorBuilder.getOrbitalParametersDrivers().getDrivers();
1186 drivers.forEach(driver -> driver.setSelected(true));
1187 drivers.forEach(driver -> driver.setValue(1.0001 * driver.getValue()));
1188
1189
1190 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
1191 propagatorBuilder);
1192 estimator.setParametersConvergenceThreshold(1.0e-2);
1193 estimator.setMaxIterations(10);
1194 estimator.setMaxEvaluations(20);
1195
1196
1197
1198 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
1199 propagatorBuilder);
1200 final List<ObservedMeasurement<?>> measurements =
1201 EstimationTestUtils.createMeasurements(propagator,
1202 new PVMeasurementCreator(),
1203 0.0, 1.0, 300.0);
1204
1205
1206 measurements.forEach(measurement -> measurement.setEnabled(false));
1207
1208
1209 measurements.forEach(estimator::addMeasurement);
1210
1211
1212 Assertions.assertThrows(IndexOutOfBoundsException.class, estimator::estimate);
1213 }
1214
1215 private void doTestEstimateOnlySomeOrbitalParameters(boolean[] orbitalParametersEstimated) {
1216
1217 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
1218
1219 final NumericalPropagatorBuilder propagatorBuilder =
1220 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1221 1.0e-6, 60.0, 1.0e-3);
1222
1223
1224 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
1225 propagatorBuilder);
1226 final List<ObservedMeasurement<?>> measurements =
1227 EstimationTestUtils.createMeasurements(propagator,
1228 new PVMeasurementCreator(),
1229 0.0, 1.0, 300.0);
1230
1231 List<ParameterDriversList.DelegatingDriver> orbitalElementsDrivers = propagatorBuilder.getOrbitalParametersDrivers().getDrivers();
1232 IntStream.range(0, orbitalParametersEstimated.length).forEach(i -> {
1233 final ParameterDriver driver = orbitalElementsDrivers.get(i);
1234 if (orbitalParametersEstimated[i]) {
1235 driver.setSelected(true);
1236 driver.setValue(1.0001 * driver.getReferenceValue());
1237 } else {
1238 driver.setSelected(false);
1239 }
1240 });
1241
1242
1243 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
1244 propagatorBuilder);
1245 measurements.forEach(estimator::addMeasurement);
1246 estimator.setParametersConvergenceThreshold(1.0e-2);
1247 estimator.setMaxIterations(10);
1248 estimator.setMaxEvaluations(20);
1249
1250
1251 estimator.estimate();
1252
1253
1254 IntStream.range(0, orbitalParametersEstimated.length).forEach(i -> {
1255 if (orbitalParametersEstimated[i]) {
1256 final ParameterDriver driver = orbitalElementsDrivers.get(i);
1257 Assertions.assertEquals(driver.getReferenceValue(), driver.getValue(),
1258 driver.getReferenceValue() * 1.0e-3);
1259 }
1260 });
1261
1262 }
1263
1264
1265
1266
1267
1268
1269 private List<ObservedMeasurement<?>> multiplexMeasurements(final List<ObservedMeasurement<?>> independentMeasurements,
1270 final double tol) {
1271 final List<ObservedMeasurement<?>> multiplexed = new ArrayList<>();
1272 independentMeasurements.sort(new ChronologicalComparator());
1273 List<ObservedMeasurement<?>> clump = new ArrayList<>();
1274 for (final ObservedMeasurement<?> measurement : independentMeasurements) {
1275 if (!clump.isEmpty() && measurement.getDate().durationFrom(clump.get(0).getDate()) > tol) {
1276
1277
1278 if (clump.size() == 1) {
1279 multiplexed.add(clump.get(0));
1280 } else {
1281 multiplexed.add(new MultiplexedMeasurement(clump));
1282 }
1283
1284
1285 clump = new ArrayList<>();
1286
1287 }
1288 clump.add(measurement);
1289 }
1290
1291 if (clump.size() == 1) {
1292 multiplexed.add(clump.get(0));
1293 } else {
1294 multiplexed.add(new MultiplexedMeasurement(clump));
1295 }
1296 return multiplexed;
1297 }
1298
1299 }
1300
1301