1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.estimation.sequential;
18
19 import org.hipparchus.geometry.euclidean.threed.Vector3D;
20 import org.hipparchus.linear.LUDecomposition;
21 import org.hipparchus.linear.MatrixUtils;
22 import org.hipparchus.linear.RealMatrix;
23 import org.hipparchus.linear.RealVector;
24 import org.junit.jupiter.api.Assertions;
25 import org.junit.jupiter.api.BeforeEach;
26 import org.junit.jupiter.api.Test;
27 import org.orekit.estimation.Context;
28 import org.orekit.estimation.EstimationTestUtils;
29 import org.orekit.estimation.Force;
30 import org.orekit.estimation.measurements.EstimatedMeasurement;
31 import org.orekit.estimation.measurements.GroundStation;
32 import org.orekit.estimation.measurements.ObservableSatellite;
33 import org.orekit.estimation.measurements.ObservedMeasurement;
34 import org.orekit.estimation.measurements.PV;
35 import org.orekit.estimation.measurements.Range;
36 import org.orekit.estimation.measurements.modifiers.Bias;
37 import org.orekit.forces.maneuvers.ConstantThrustManeuver;
38 import org.orekit.forces.radiation.RadiationSensitive;
39 import org.orekit.orbits.Orbit;
40 import org.orekit.orbits.OrbitType;
41 import org.orekit.orbits.PositionAngleType;
42 import org.orekit.propagation.MatricesHarvester;
43 import org.orekit.propagation.Propagator;
44 import org.orekit.propagation.SpacecraftState;
45 import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
46 import org.orekit.time.AbsoluteDate;
47 import org.orekit.utils.PVCoordinates;
48 import org.orekit.utils.ParameterDriver;
49 import org.orekit.utils.ParameterDriversList;
50 import org.orekit.utils.TimeSpanMap.Span;
51
52 import java.util.ArrayList;
53 import java.util.Collections;
54 import java.util.List;
55
56
57
58
59
60
61
62
63 public class KalmanModelTest {
64
65
66 private final OrbitType orbitType = OrbitType.CARTESIAN;
67
68
69 private final PositionAngleType positionAngleType = PositionAngleType.TRUE;
70
71
72 private Orbit orbit0;
73
74
75 private NumericalPropagatorBuilder propagatorBuilder;
76
77
78 private CovarianceMatrixProvider covMatrixProvider;
79
80
81 private ParameterDriversList estimatedMeasurementsParameters;
82
83
84 private KalmanEstimator kalman;
85
86
87 private ModelLogger modelLogger;
88
89
90 private int M;
91
92
93 private PV pv;
94
95
96 private Range range;
97
98
99 private ParameterDriver satRangeBiasDriver;
100
101
102 private ParameterDriver srpCoefDriver;
103
104
105 private final double tol = 2e-16;
106
107
108
109
110
111
112
113
114
115
116 @BeforeEach
117 public void setup() {
118
119 final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
120
121
122 this.orbit0 = context.initialOrbit;
123 ObservableSatellite sat = new ObservableSatellite(0);
124
125
126 this.propagatorBuilder = context.createBuilder(orbitType, positionAngleType, true,
127 1.0e-6, 60.0, 10., Force.SOLAR_RADIATION_PRESSURE);
128
129
130 final AbsoluteDate date0 = context.initialOrbit.getDate();
131 this.pv = new PV(date0,
132 context.initialOrbit.getPosition(),
133 context.initialOrbit.getPVCoordinates().getVelocity(),
134 new double[] {1., 2., 3., 1e-3, 2e-3, 3e-3}, 1.,
135 sat);
136
137
138 final AbsoluteDate date = date0.shiftedBy(10.);
139 final GroundStation station = context.stations.get(0);
140 this.range = new Range(station, true, date, 18616150., 10., 1., sat);
141
142
143
144 final Bias<Range> satRangeBias = new Bias<Range>(new String[] {"sat range bias"},
145 new double[] {100.},
146 new double[] {10.},
147 new double[] {0.},
148 new double[] {100.});
149 this.satRangeBiasDriver = satRangeBias.getParametersDrivers().get(0);
150 satRangeBiasDriver.setSelected(true);
151 satRangeBiasDriver.setReferenceDate(date);
152 range.addModifier(satRangeBias);
153 for (ParameterDriver driver : range.getParametersDrivers()) {
154 driver.setReferenceDate(date);
155 }
156
157
158 this.estimatedMeasurementsParameters = new ParameterDriversList();
159 for (final ParameterDriver driver : range.getParametersDrivers()) {
160 if (driver.isSelected()) {
161 estimatedMeasurementsParameters.add(driver);
162 }
163 }
164
165 this.srpCoefDriver = propagatorBuilder.getPropagationParametersDrivers().
166 findByName(RadiationSensitive.REFLECTION_COEFFICIENT);
167 srpCoefDriver.setReferenceDate(date);
168 srpCoefDriver.setSelected(true);
169
170
171 final double[] scales = getParametersScale(propagatorBuilder, estimatedMeasurementsParameters);
172 this.M = scales.length;
173 this.covMatrixProvider = setInitialCovarianceMatrix(scales);
174
175
176 final KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
177 kalmanBuilder.addPropagationConfiguration(propagatorBuilder, covMatrixProvider);
178 kalmanBuilder.estimatedMeasurementsParameters(estimatedMeasurementsParameters, null);
179 this.kalman = kalmanBuilder.build();
180 this.modelLogger = new ModelLogger();
181 kalman.setObserver(modelLogger);
182 }
183
184
185
186
187
188
189 @Test
190 public void ModelPhysicalOutputsTest() {
191
192
193
194 checkModelAtT0();
195
196
197
198
199
200
201 final RealMatrix Q = covMatrixProvider.getProcessNoiseMatrix(new SpacecraftState(orbit0),
202 new SpacecraftState(orbit0));
203
204
205 final RealMatrix P0 = covMatrixProvider.getInitialCovarianceMatrix(new SpacecraftState(orbit0));
206
207
208
209 RealMatrix Ppred = P0.add(Q);
210
211
212 Orbit orbitPred = orbit0;
213
214
215
216 RealMatrix expH = MatrixUtils.createRealMatrix(6, M);
217 for (int i = 0; i < 6; i++) {
218 expH.setEntry(i, i, 1.);
219 }
220
221
222
223 RealMatrix expPhi = MatrixUtils.createRealIdentityMatrix(M);
224
225 checkModelAfterMeasurementAdded(1, pv, Ppred, orbitPred, expPhi, expH);
226
227
228
229
230
231
232
233 final Propagator propagator = propagatorBuilder.buildPropagator();
234
235
236 final String equationName = KalmanEstimator.class.getName() + "-derivatives-";
237 final MatricesHarvester harvester = propagator.setupMatricesComputation(equationName, null, null);
238
239
240 final SpacecraftState scPred = propagator.propagate(range.getDate());
241 orbitPred = scPred.getOrbit();
242
243
244 expPhi = MatrixUtils.createRealIdentityMatrix(M);
245
246
247 final double[][] dYdY0 = harvester.getStateTransitionMatrix(scPred).getData();
248 expPhi.setSubMatrix(dYdY0, 0, 0);
249
250
251 final double[][] dYdPp = harvester.getParametersJacobian(scPred).getData();
252 expPhi.setSubMatrix(dYdPp, 0, 6);
253
254
255 RealMatrix Pest = kalman.getPhysicalEstimatedCovarianceMatrix();
256
257
258 Ppred = expPhi.multiply(Pest.multiplyTransposed(expPhi)).add(Q);
259
260
261 expH = MatrixUtils.createRealMatrix(1, M);
262
263 EstimatedMeasurement<Range> rangeEstimated = range.estimate(0, 0, new SpacecraftState[] {scPred});
264 final RealMatrix dMdY = MatrixUtils.createRealMatrix(rangeEstimated.getStateDerivatives(0));
265 expH.setSubMatrix(dMdY.getData(), 0, 0);
266
267 final SpacecraftState scTransition = scPred.shiftedBy(-rangeEstimated.getTimeOffset());
268 final double[][] dYdPpTransition = harvester.getParametersJacobian(scTransition).getData();
269 final RealMatrix dMdCr = dMdY.multiply(MatrixUtils.createRealMatrix(dYdPpTransition));
270 expH.setEntry(0, 6, dMdCr.getEntry(0, 0));
271
272 expH.setEntry(0, 7, rangeEstimated.getParameterDerivatives(satRangeBiasDriver, new AbsoluteDate())[0]);
273
274
275 checkModelAfterMeasurementAdded(2, range, Ppred, orbitPred, expPhi, expH);
276 }
277
278
279 private void checkModelAtT0() {
280
281
282 final KalmanModel model = new KalmanModel(Collections.singletonList(propagatorBuilder),
283 Collections.singletonList(covMatrixProvider),
284 estimatedMeasurementsParameters,
285 null);
286
287
288
289
290
291 Assertions.assertEquals(0., model.getEstimate().getTime(), 0.);
292 Assertions.assertEquals(0., model.getCurrentDate().durationFrom(orbit0.getDate()), 0.);
293
294
295 Assertions.assertEquals(0, model.getCurrentMeasurementNumber());
296
297
298 final RealVector stateN = model.getEstimate().getState();
299 Assertions.assertArrayEquals(new double[M], stateN.toArray(), tol);
300
301
302 final RealVector x = model.getPhysicalEstimatedState();
303 final RealVector expX = MatrixUtils.createRealVector(M);
304 final double[] orbitState0 = new double[6];
305 orbitType.mapOrbitToArray(orbit0, positionAngleType, orbitState0, null);
306 expX.setSubVector(0, MatrixUtils.createRealVector(orbitState0));
307 expX.setEntry(6, srpCoefDriver.getReferenceValue());
308 expX.setEntry(7, satRangeBiasDriver.getReferenceValue());
309 final double[] dX = x.subtract(expX).toArray();
310 Assertions.assertArrayEquals(new double[M], dX, tol);
311
312
313 final double[][] Pn = model.getEstimate().getCovariance().getData();
314 final double[][] expPn = new double[M][M];
315 for (int i = 0; i < M; i++) {
316 for (int j = 0; j < M; j++) {
317 expPn[i][j] = 1.;
318 }
319 Assertions.assertArrayEquals(expPn[i], Pn[i], tol, "Failed on line " + i);
320 }
321
322
323 final RealMatrix P = model.getPhysicalEstimatedCovarianceMatrix();
324 final RealMatrix expP = covMatrixProvider.getInitialCovarianceMatrix(new SpacecraftState(orbit0));
325 final double[][] dP = P.subtract(expP).getData();
326 for (int i = 0; i < M; i++) {
327 Assertions.assertArrayEquals(new double[M], dP[i], tol, "Failed on line " + i);
328 }
329
330
331 Assertions.assertNull(model.getEstimate().getInnovationCovariance());
332 Assertions.assertNull(model.getPhysicalInnovationCovarianceMatrix());
333 Assertions.assertNull(model.getEstimate().getKalmanGain());
334 Assertions.assertNull(model.getPhysicalKalmanGain());
335 Assertions.assertNull(model.getEstimate().getMeasurementJacobian());
336 Assertions.assertNull(model.getPhysicalMeasurementJacobian());
337 Assertions.assertNull(model.getEstimate().getStateTransitionMatrix());
338 Assertions.assertNull(model.getPhysicalStateTransitionMatrix());
339 }
340
341
342
343
344 private void checkModelAfterMeasurementAdded(final int expMeasurementNumber,
345 final ObservedMeasurement<?> meas,
346 final RealMatrix expPpred,
347 final Orbit expOrbitPred,
348 final RealMatrix expPhi,
349 final RealMatrix expH) {
350
351
352
353 final double srpCoefPred = srpCoefDriver.getValue();
354 final double satRangeBiasPred = satRangeBiasDriver.getValue();
355
356
357 final double[] expMeasPred =
358 meas.estimate(0, 0,
359 new SpacecraftState[] {new SpacecraftState(expOrbitPred)}).getEstimatedValue();
360
361
362 kalman.processMeasurements(Collections.singletonList(meas));
363 KalmanEstimation model = modelLogger.estimation;
364
365
366 final int N = meas.getDimension();
367
368
369 Assertions.assertEquals(0., model.getCurrentDate().durationFrom(expOrbitPred.getDate()), 0.);
370
371
372 Assertions.assertEquals(expMeasurementNumber, model.getCurrentMeasurementNumber());
373
374
375 final RealMatrix phi = model.getPhysicalStateTransitionMatrix();
376 final double[][] dPhi = phi.subtract(expPhi).getData();
377 for (int i = 0; i < M; i++) {
378 Assertions.assertArrayEquals(new double[M], dPhi[i], tol*100, "Failed on line " + i);
379 }
380
381
382 final RealMatrix H = model.getPhysicalMeasurementJacobian();
383 final double[][] dH = H.subtract(expH).getData();
384 for (int i = 0; i < N; i++) {
385 Assertions.assertArrayEquals(new double[M], dH[i], tol, "Failed on line " + i);
386 }
387
388
389 final double[] measSigmas = meas.getTheoreticalStandardDeviation();
390 final RealMatrix R = MatrixUtils.createRealMatrix(N, N);
391 for (int i = 0; i < N; i++) {
392 R.setEntry(i, i, measSigmas[i] * measSigmas[i]);
393 }
394
395
396 final RealMatrix expS = expH.multiply(expPpred.multiplyTransposed(expH)).add(R);
397 final RealMatrix S = model.getPhysicalInnovationCovarianceMatrix();
398 final double[][] dS = S.subtract(expS).getData();
399 for (int i = 0; i < N; i++) {
400 Assertions.assertArrayEquals(new double[N], dS[i], tol*1e4, "Failed on line \" + i");
401 }
402
403
404 final RealMatrix expK = expPpred.multiplyTransposed(expH).multiply(new LUDecomposition(expS).getSolver().getInverse());
405 final RealMatrix K = model.getPhysicalKalmanGain();
406 final double[][] dK = K.subtract(expK).getData();
407 for (int i = 0; i < M; i++) {
408 Assertions.assertArrayEquals(new double[N], dK[i], tol*1e5, "Failed on line " + i);
409 }
410
411
412 final Orbit orbitPred = model.getPredictedSpacecraftStates()[0].getOrbit();
413 final PVCoordinates pvOrbitPred = orbitPred.getPVCoordinates();
414 final PVCoordinates expPVOrbitPred = expOrbitPred.getPVCoordinates();
415 final double dpOrbitPred = Vector3D.distance(expPVOrbitPred.getPosition(), pvOrbitPred.getPosition());
416 final double dvOrbitPred = Vector3D.distance(expPVOrbitPred.getVelocity(), pvOrbitPred.getVelocity());
417 Assertions.assertEquals(0., dpOrbitPred, tol);
418 Assertions.assertEquals(0., dvOrbitPred, tol);
419
420
421 final double[] measPred = model.getPredictedMeasurement().getEstimatedValue();
422 Assertions.assertArrayEquals(expMeasPred, measPred, tol);
423
424
425 final double[] orbitPredState = new double[6];
426 orbitPred.getType().mapOrbitToArray(orbitPred, PositionAngleType.TRUE, orbitPredState, null);
427 final RealVector expXpred = MatrixUtils.createRealVector(M);
428 for (int i = 0; i < 6; i++) {
429 expXpred.setEntry(i, orbitPredState[i]);
430 }
431 expXpred.setEntry(6, srpCoefPred);
432 expXpred.setEntry(7, satRangeBiasPred);
433
434
435 final RealVector observedMeas = MatrixUtils.createRealVector(model.getPredictedMeasurement().getObservedValue());
436 final RealVector predictedMeas = MatrixUtils.createRealVector(model.getPredictedMeasurement().getEstimatedValue());
437 final RealVector innovation = observedMeas.subtract(predictedMeas);
438
439
440 final RealVector expectedXcor = expXpred.add(expK.operate(innovation));
441 final RealVector Xcor = model.getPhysicalEstimatedState();
442 final double[] dXcor = Xcor.subtract(expectedXcor).toArray();
443 Assertions.assertArrayEquals(new double[M], dXcor, tol);
444
445
446 final RealMatrix expectedPcor =
447 (MatrixUtils.createRealIdentityMatrix(M).
448 subtract(expK.multiply(expH))).multiply(expPpred);
449 final RealMatrix Pcor = model.getPhysicalEstimatedCovarianceMatrix();
450 final double[][] dPcor = Pcor.subtract(expectedPcor).getData();
451 for (int i = 0; i < M; i++) {
452 Assertions.assertArrayEquals(new double[M], dPcor[i], tol*1e5, "Failed on line " + i);
453 }
454 }
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498 private double[] getParametersScale(final NumericalPropagatorBuilder builder,
499 ParameterDriversList estimatedMeasurementsParameters) {
500 final List<Double> scaleList = new ArrayList<>();
501
502
503 for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
504 if (driver.isSelected()) {
505 for (Span<Double> span = driver.getValueSpanMap().getFirstSpan(); span != null; span = span.next()) {
506 scaleList.add(driver.getScale());
507 }
508 }
509 }
510
511
512 for (ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) {
513 if (driver.isSelected()) {
514 for (Span<Double> span = driver.getValueSpanMap().getFirstSpan(); span != null; span = span.next()) {
515 scaleList.add(driver.getScale());
516 }
517 }
518 }
519
520
521 for (ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
522 if (driver.isSelected()) {
523 for (Span<Double> span = driver.getValueSpanMap().getFirstSpan(); span != null; span = span.next()) {
524 scaleList.add(driver.getScale());
525 }
526 }
527 }
528
529 final double[] scales = new double[scaleList.size()];
530 for (int i = 0; i < scaleList.size(); i++) {
531 scales[i] = scaleList.get(i);
532 }
533 return scales;
534 }
535
536
537
538
539
540
541
542 private CovarianceMatrixProvider setInitialCovarianceMatrix(final double[] scales) {
543
544 final int n = scales.length;
545 final RealMatrix cov = MatrixUtils.createRealMatrix(n, n);
546 for (int i = 0; i < n; i++) {
547 for (int j = 0; j < n; j++) {
548 cov.setEntry(i, j, scales[i] * scales[j]);
549 }
550 }
551 return new ConstantProcessNoise(cov);
552 }
553
554
555
556 public class ModelLogger implements KalmanObserver {
557 KalmanEstimation estimation;
558
559 @Override
560 public void evaluationPerformed(KalmanEstimation estimation) {
561 this.estimation = estimation;
562 }
563 }
564
565
566
567
568
569 @Test
570 public void MassDepletionTest() {
571
572 AbsoluteDate initialDate = this.propagatorBuilder.getInitialOrbitDate();
573 double initialMass = this.propagatorBuilder.getMass();
574 AbsoluteDate maneuverDate = initialDate.shiftedBy(5.0);
575
576 ConstantThrustManeuver maneuver = new ConstantThrustManeuver(maneuverDate, 10.0, 10.0, 100.0, new Vector3D(1.0, 0.0, 0.0));
577 this.propagatorBuilder.addForceModel(maneuver);
578
579 kalman.processMeasurements(Collections.singletonList(pv));
580 kalman.processMeasurements(Collections.singletonList(range));
581
582 double postManeuverPredictedStateMass = modelLogger.estimation.getPredictedSpacecraftStates()[0].getMass();
583 double postManeuverCorrectedStateMass = modelLogger.estimation.getCorrectedSpacecraftStates()[0].getMass();
584 double postManeuverPropagatorMass = propagatorBuilder.getMass();
585
586
587 Assertions.assertTrue(postManeuverPredictedStateMass < initialMass);
588
589
590
591 Assertions.assertEquals(postManeuverPredictedStateMass, postManeuverPropagatorMass);
592
593
594 Assertions.assertEquals(postManeuverCorrectedStateMass, postManeuverPropagatorMass);
595 }
596 }