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.filtering.kalman.ProcessEstimate;
20 import org.hipparchus.filtering.kalman.unscented.UnscentedEvolution;
21 import org.hipparchus.filtering.kalman.unscented.UnscentedKalmanFilter;
22 import org.hipparchus.filtering.kalman.unscented.UnscentedProcess;
23 import org.hipparchus.linear.ArrayRealVector;
24 import org.hipparchus.linear.MatrixUtils;
25 import org.hipparchus.linear.RealMatrix;
26 import org.hipparchus.linear.RealVector;
27 import org.hipparchus.util.FastMath;
28 import org.orekit.estimation.measurements.EstimatedMeasurement;
29 import org.orekit.estimation.measurements.EstimatedMeasurementBase;
30 import org.orekit.estimation.measurements.ObservedMeasurement;
31 import org.orekit.orbits.Orbit;
32 import org.orekit.orbits.OrbitType;
33 import org.orekit.orbits.PositionAngleType;
34 import org.orekit.propagation.PropagationType;
35 import org.orekit.propagation.SpacecraftState;
36 import org.orekit.propagation.conversion.DSSTPropagatorBuilder;
37 import org.orekit.propagation.semianalytical.dsst.DSSTPropagator;
38 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
39 import org.orekit.propagation.semianalytical.dsst.forces.ShortPeriodTerms;
40 import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
41 import org.orekit.time.AbsoluteDate;
42 import org.orekit.time.ChronologicalComparator;
43 import org.orekit.utils.ParameterDriver;
44 import org.orekit.utils.ParameterDriversList;
45 import org.orekit.utils.ParameterDriversList.DelegatingDriver;
46
47 import java.util.ArrayList;
48 import java.util.Comparator;
49 import java.util.List;
50
51
52
53
54
55
56 public class SemiAnalyticalUnscentedKalmanModel implements KalmanEstimation, UnscentedProcess<MeasurementDecorator>, SemiAnalyticalProcess {
57
58
59 private final DSSTPropagatorBuilder builder;
60
61
62 private final ParameterDriversList estimatedOrbitalParameters;
63
64
65 private final ParameterDriversList estimatedPropagationParameters;
66
67
68 private final ParameterDriversList estimatedMeasurementsParameters;
69
70
71 private final CovarianceMatrixProvider covarianceMatrixProvider;
72
73
74 private final CovarianceMatrixProvider measurementProcessNoiseMatrix;
75
76
77 private final PositionAngleType angleType;
78
79
80 private final OrbitType orbitType;
81
82
83 private ProcessEstimate correctedEstimate;
84
85
86 private KalmanObserver observer;
87
88
89 private int currentMeasurementNumber;
90
91
92 private AbsoluteDate currentDate;
93
94
95 private SpacecraftState nominalMeanSpacecraftState;
96
97
98 private SpacecraftState previousNominalMeanSpacecraftState;
99
100
101 private SpacecraftState predictedSpacecraftState;
102
103
104 private SpacecraftState correctedSpacecraftState;
105
106
107 private EstimatedMeasurement<?> predictedMeasurement;
108
109
110 private EstimatedMeasurement<?> correctedMeasurement;
111
112
113 private RealVector predictedFilterCorrection;
114
115
116 private RealVector correctedFilterCorrection;
117
118
119 private final DSSTPropagator dsstPropagator;
120
121
122 private RealVector shortPeriodicTerms;
123
124
125
126
127
128
129
130 protected SemiAnalyticalUnscentedKalmanModel(final DSSTPropagatorBuilder propagatorBuilder,
131 final CovarianceMatrixProvider covarianceMatrixProvider,
132 final ParameterDriversList estimatedMeasurementParameters,
133 final CovarianceMatrixProvider measurementProcessNoiseMatrix) {
134
135 this.builder = propagatorBuilder;
136 this.angleType = propagatorBuilder.getPositionAngleType();
137 this.orbitType = propagatorBuilder.getOrbitType();
138 this.estimatedMeasurementsParameters = estimatedMeasurementParameters;
139 this.currentMeasurementNumber = 0;
140 this.currentDate = propagatorBuilder.getInitialOrbitDate();
141 this.covarianceMatrixProvider = covarianceMatrixProvider;
142 this.measurementProcessNoiseMatrix = measurementProcessNoiseMatrix;
143
144
145 int columns = 0;
146
147
148 this.estimatedOrbitalParameters = new ParameterDriversList();
149 for (final ParameterDriver driver : propagatorBuilder.getOrbitalParametersDrivers().getDrivers()) {
150
151
152 if (driver.getReferenceDate() == null) {
153 driver.setReferenceDate(currentDate);
154 }
155
156
157 if (driver.isSelected()) {
158 estimatedOrbitalParameters.add(driver);
159 columns++;
160 }
161
162 }
163
164
165 this.estimatedPropagationParameters = new ParameterDriversList();
166 final List<String> estimatedPropagationParametersNames = new ArrayList<>();
167 for (final ParameterDriver driver : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
168
169
170 if (driver.getReferenceDate() == null) {
171 driver.setReferenceDate(currentDate);
172 }
173
174
175 if (driver.isSelected()) {
176 estimatedPropagationParameters.add(driver);
177 final String driverName = driver.getName();
178
179 if (!estimatedPropagationParametersNames.contains(driverName)) {
180 estimatedPropagationParametersNames.add(driverName);
181 ++columns;
182 }
183 }
184
185 }
186 estimatedPropagationParametersNames.sort(Comparator.naturalOrder());
187
188
189 for (final ParameterDriver parameter : estimatedMeasurementsParameters.getDrivers()) {
190 if (parameter.getReferenceDate() == null) {
191 parameter.setReferenceDate(currentDate);
192 }
193 ++columns;
194 }
195
196
197 final int nbMeas = estimatedMeasurementParameters.getNbParams();
198
199
200 final int nbDyn = estimatedOrbitalParameters.getNbParams() +
201 estimatedPropagationParameters.getNbParams();
202
203
204 this.dsstPropagator = getEstimatedPropagator();
205 final SpacecraftState meanState = dsstPropagator.initialIsOsculating() ?
206 DSSTPropagator.computeMeanState(dsstPropagator.getInitialState(), dsstPropagator.getAttitudeProvider(), dsstPropagator.getAllForceModels()) :
207 dsstPropagator.getInitialState();
208 this.nominalMeanSpacecraftState = meanState;
209 this.predictedSpacecraftState = meanState;
210 this.correctedSpacecraftState = predictedSpacecraftState;
211 this.previousNominalMeanSpacecraftState = nominalMeanSpacecraftState;
212
213
214 this.predictedFilterCorrection = MatrixUtils.createRealVector(columns);
215 this.correctedFilterCorrection = predictedFilterCorrection;
216
217
218 final RealMatrix noiseK = MatrixUtils.createRealMatrix(nbDyn + nbMeas, nbDyn + nbMeas);
219 final RealMatrix noiseP = covarianceMatrixProvider.getInitialCovarianceMatrix(nominalMeanSpacecraftState);
220 noiseK.setSubMatrix(noiseP.getData(), 0, 0);
221 if (measurementProcessNoiseMatrix != null) {
222 final RealMatrix noiseM = measurementProcessNoiseMatrix.getInitialCovarianceMatrix(nominalMeanSpacecraftState);
223 noiseK.setSubMatrix(noiseM.getData(), nbDyn, nbDyn);
224 }
225
226 KalmanEstimatorUtil.checkDimension(noiseK.getRowDimension(),
227 propagatorBuilder.getOrbitalParametersDrivers(),
228 propagatorBuilder.getPropagationParametersDrivers(),
229 estimatedMeasurementsParameters);
230
231
232 this.correctedEstimate = new ProcessEstimate(0.0, correctedFilterCorrection, noiseK);
233
234 }
235
236
237 @Override
238 public KalmanObserver getObserver() {
239 return observer;
240 }
241
242
243
244
245 public void setObserver(final KalmanObserver observer) {
246 this.observer = observer;
247 }
248
249
250
251
252
253
254
255
256
257 public ProcessEstimate getEstimate() {
258 return correctedEstimate;
259 }
260
261
262
263
264
265
266 public DSSTPropagator processMeasurements(final List<ObservedMeasurement<?>> observedMeasurements,
267 final UnscentedKalmanFilter<MeasurementDecorator> filter) {
268
269
270 observedMeasurements.sort(new ChronologicalComparator());
271 final AbsoluteDate tStart = observedMeasurements.get(0).getDate();
272 final AbsoluteDate tEnd = observedMeasurements.get(observedMeasurements.size() - 1).getDate();
273 final double overshootTimeRange = FastMath.nextAfter(tEnd.durationFrom(tStart),
274 Double.POSITIVE_INFINITY);
275
276
277 final SemiAnalyticalMeasurementHandler stepHandler = new SemiAnalyticalMeasurementHandler(this, filter, observedMeasurements, builder.getInitialOrbitDate(), true);
278 dsstPropagator.getMultiplexer().add(stepHandler);
279 dsstPropagator.propagate(tStart, tStart.shiftedBy(overshootTimeRange));
280
281
282 return getEstimatedPropagator();
283
284 }
285
286
287
288
289 public DSSTPropagator getEstimatedPropagator() {
290
291 return (DSSTPropagator) builder.buildPropagator();
292 }
293
294
295 @Override
296 public UnscentedEvolution getEvolution(final double previousTime, final RealVector[] sigmaPoints,
297 final MeasurementDecorator measurement) {
298
299
300 final ObservedMeasurement<?> observedMeasurement = measurement.getObservedMeasurement();
301 for (final ParameterDriver driver : observedMeasurement.getParametersDrivers()) {
302 if (driver.getReferenceDate() == null) {
303 driver.setReferenceDate(builder.getInitialOrbitDate());
304 }
305 }
306
307
308 ++currentMeasurementNumber;
309
310
311 currentDate = measurement.getObservedMeasurement().getDate();
312
313
314 final RealMatrix stm = getStm();
315
316
317 final RealVector[] predictedStates = new RealVector[sigmaPoints.length];
318 for (int k = 0; k < sigmaPoints.length; ++k) {
319
320 final RealVector predicted = stm.operate(sigmaPoints[k]);
321 predictedStates[k] = predicted;
322 }
323
324
325 return new UnscentedEvolution(measurement.getTime(), predictedStates);
326
327 }
328
329
330 @Override
331 public RealMatrix getProcessNoiseMatrix(final double previousTime, final RealVector predictedState,
332 final MeasurementDecorator measurement) {
333
334
335 final int nbMeas = getNumberSelectedMeasurementDrivers();
336
337
338 final int nbDyn = getNumberSelectedOrbitalDrivers() + getNumberSelectedPropagationDrivers();
339
340
341 final RealMatrix noiseK = MatrixUtils.createRealMatrix(nbDyn + nbMeas, nbDyn + nbMeas);
342 final RealMatrix noiseP = covarianceMatrixProvider.getProcessNoiseMatrix(previousNominalMeanSpacecraftState, nominalMeanSpacecraftState);
343 noiseK.setSubMatrix(noiseP.getData(), 0, 0);
344 if (measurementProcessNoiseMatrix != null) {
345 final RealMatrix noiseM = measurementProcessNoiseMatrix.getProcessNoiseMatrix(previousNominalMeanSpacecraftState, nominalMeanSpacecraftState);
346 noiseK.setSubMatrix(noiseM.getData(), nbDyn, nbDyn);
347 }
348
349
350 KalmanEstimatorUtil.checkDimension(noiseK.getRowDimension(),
351 builder.getOrbitalParametersDrivers(),
352 builder.getPropagationParametersDrivers(),
353 estimatedMeasurementsParameters);
354
355 return noiseK;
356 }
357
358
359 @Override
360 public RealVector[] getPredictedMeasurements(final RealVector[] predictedSigmaPoints, final MeasurementDecorator measurement) {
361
362
363 final ObservedMeasurement<?> observedMeasurement = measurement.getObservedMeasurement();
364
365
366 final RealVector[] predictedMeasurements = new RealVector[predictedSigmaPoints.length];
367
368
369 for (int k = 0; k < predictedSigmaPoints.length; ++k) {
370
371
372 final RealVector osculating = computeOsculatingElements(predictedSigmaPoints[k], nominalMeanSpacecraftState, shortPeriodicTerms);
373 final Orbit osculatingOrbit = orbitType.mapArrayToOrbit(osculating.toArray(), null, angleType,
374 currentDate, builder.getMu(), builder.getFrame());
375
376
377 final EstimatedMeasurement<?> estimated = estimateMeasurement(observedMeasurement, currentMeasurementNumber,
378 new SpacecraftState[] { new SpacecraftState(osculatingOrbit) });
379 predictedMeasurements[k] = new ArrayRealVector(estimated.getEstimatedValue());
380
381 }
382
383
384 return predictedMeasurements;
385
386 }
387
388
389 @Override
390 public RealVector getInnovation(final MeasurementDecorator measurement, final RealVector predictedMeas,
391 final RealVector predictedState, final RealMatrix innovationCovarianceMatrix) {
392
393
394 predictedFilterCorrection = predictedState;
395
396
397 final RealVector osculating = computeOsculatingElements(predictedFilterCorrection, nominalMeanSpacecraftState, shortPeriodicTerms);
398 final Orbit osculatingOrbit = orbitType.mapArrayToOrbit(osculating.toArray(), null, angleType,
399 currentDate, builder.getMu(), builder.getFrame());
400 predictedSpacecraftState = new SpacecraftState(osculatingOrbit);
401 predictedMeasurement = estimateMeasurement(measurement.getObservedMeasurement(), currentMeasurementNumber,
402 getPredictedSpacecraftStates());
403 predictedMeasurement.setEstimatedValue(predictedMeas.toArray());
404
405
406 KalmanEstimatorUtil.applyDynamicOutlierFilter(predictedMeasurement, innovationCovarianceMatrix);
407
408
409 return KalmanEstimatorUtil.computeInnovationVector(predictedMeasurement);
410
411 }
412
413
414
415 @Override
416 public void finalizeEstimation(final ObservedMeasurement<?> observedMeasurement,
417 final ProcessEstimate estimate) {
418
419 correctedEstimate = estimate;
420
421 correctedFilterCorrection = estimate.getState();
422
423
424 previousNominalMeanSpacecraftState = nominalMeanSpacecraftState;
425
426
427
428 final RealVector osculating = computeOsculatingElements(correctedFilterCorrection, nominalMeanSpacecraftState, shortPeriodicTerms);
429 final Orbit osculatingOrbit = orbitType.mapArrayToOrbit(osculating.toArray(), null, builder.getPositionAngleType(),
430 currentDate, builder.getMu(), builder.getFrame());
431
432
433 correctedSpacecraftState = new SpacecraftState(osculatingOrbit);
434 correctedMeasurement = estimateMeasurement(observedMeasurement, currentMeasurementNumber,
435 getCorrectedSpacecraftStates());
436
437
438 if (observer != null) {
439 observer.evaluationPerformed(this);
440 }
441 }
442
443
444
445
446
447
448
449
450
451
452 private static <T extends ObservedMeasurement<T>> EstimatedMeasurement<?> estimateMeasurement(final ObservedMeasurement<T> observedMeasurement,
453 final int measurementNumber,
454 final SpacecraftState[] spacecraftStates) {
455 final EstimatedMeasurementBase<T> estimatedMeasurementBase = observedMeasurement.
456 estimateWithoutDerivatives(measurementNumber, measurementNumber,
457 KalmanEstimatorUtil.filterRelevant(observedMeasurement, spacecraftStates));
458 final EstimatedMeasurement<T> estimatedMeasurement = new EstimatedMeasurement<>(estimatedMeasurementBase.getObservedMeasurement(),
459 estimatedMeasurementBase.getIteration(), estimatedMeasurementBase.getCount(),
460 estimatedMeasurementBase.getStates(), estimatedMeasurementBase.getParticipants());
461 estimatedMeasurement.setEstimatedValue(estimatedMeasurementBase.getEstimatedValue());
462 return estimatedMeasurement;
463 }
464
465
466
467
468
469
470
471
472 private RealMatrix getStm() {
473
474
475 final int nbDym = getNumberSelectedOrbitalDrivers() + getNumberSelectedPropagationDrivers();
476 final int nbMeas = getNumberSelectedMeasurementDrivers();
477 final RealMatrix stm = MatrixUtils.createRealIdentityMatrix(nbDym + nbMeas);
478
479
480 final double mu = builder.getMu();
481 final double sma = previousNominalMeanSpacecraftState.getOrbit().getA();
482 final double dt = currentDate.durationFrom(previousNominalMeanSpacecraftState.getDate());
483 final double contribution = -1.5 * dt * FastMath.sqrt(mu / FastMath.pow(sma, 5));
484 stm.setEntry(5, 0, contribution);
485
486
487 return stm;
488
489 }
490
491
492 @Override
493 public void finalizeOperationsObservationGrid() {
494
495 updateParameters();
496 }
497
498
499 @Override
500 public ParameterDriversList getEstimatedOrbitalParameters() {
501 return estimatedOrbitalParameters;
502 }
503
504
505 @Override
506 public ParameterDriversList getEstimatedPropagationParameters() {
507 return estimatedPropagationParameters;
508 }
509
510
511 @Override
512 public ParameterDriversList getEstimatedMeasurementsParameters() {
513 return estimatedMeasurementsParameters;
514 }
515
516
517
518
519
520
521 @Override
522 public SpacecraftState[] getPredictedSpacecraftStates() {
523 return new SpacecraftState[] {predictedSpacecraftState};
524 }
525
526
527
528
529
530
531 @Override
532 public SpacecraftState[] getCorrectedSpacecraftStates() {
533 return new SpacecraftState[] {correctedSpacecraftState};
534 }
535
536
537 @Override
538 public RealVector getPhysicalEstimatedState() {
539
540
541
542 final RealVector physicalEstimatedState = new ArrayRealVector(getEstimate().getState().getDimension());
543 int i = 0;
544 for (final DelegatingDriver driver : getEstimatedOrbitalParameters().getDrivers()) {
545 physicalEstimatedState.setEntry(i++, driver.getValue());
546 }
547 for (final DelegatingDriver driver : getEstimatedPropagationParameters().getDrivers()) {
548 physicalEstimatedState.setEntry(i++, driver.getValue());
549 }
550 for (final DelegatingDriver driver : getEstimatedMeasurementsParameters().getDrivers()) {
551 physicalEstimatedState.setEntry(i++, driver.getValue());
552 }
553
554 return physicalEstimatedState;
555 }
556
557
558 @Override
559 public RealMatrix getPhysicalEstimatedCovarianceMatrix() {
560 return correctedEstimate.getCovariance();
561 }
562
563
564 @Override
565 public RealMatrix getPhysicalStateTransitionMatrix() {
566 return null;
567 }
568
569
570 @Override
571 public RealMatrix getPhysicalMeasurementJacobian() {
572 return null;
573 }
574
575
576 @Override
577 public RealMatrix getPhysicalInnovationCovarianceMatrix() {
578 return correctedEstimate.getInnovationCovariance();
579 }
580
581
582 @Override
583 public RealMatrix getPhysicalKalmanGain() {
584 return correctedEstimate.getKalmanGain();
585 }
586
587
588 @Override
589 public int getCurrentMeasurementNumber() {
590 return currentMeasurementNumber;
591 }
592
593
594 @Override
595 public AbsoluteDate getCurrentDate() {
596 return currentDate;
597 }
598
599
600 @Override
601 public EstimatedMeasurement<?> getPredictedMeasurement() {
602 return predictedMeasurement;
603 }
604
605
606 @Override
607 public EstimatedMeasurement<?> getCorrectedMeasurement() {
608 return correctedMeasurement;
609 }
610
611
612 @Override
613 public void updateNominalSpacecraftState(final SpacecraftState nominal) {
614 this.nominalMeanSpacecraftState = nominal;
615
616 shortPeriodicTerms = new ArrayRealVector(dsstPropagator.getShortPeriodTermsValue(nominalMeanSpacecraftState));
617
618 builder.resetOrbit(nominal.getOrbit(), PropagationType.MEAN);
619 }
620
621
622 @Override
623 public void updateShortPeriods(final SpacecraftState state) {
624
625 for (final DSSTForceModel model : dsstPropagator.getAllForceModels()) {
626 model.updateShortPeriodTerms(model.getParameters(), state);
627 }
628 }
629
630
631 @Override
632 public void initializeShortPeriodicTerms(final SpacecraftState meanState) {
633 final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<>();
634 for (final DSSTForceModel force : builder.getAllForceModels()) {
635 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(new AuxiliaryElements(meanState.getOrbit(), 1), PropagationType.OSCULATING, force.getParameters()));
636 }
637 dsstPropagator.setShortPeriodTerms(shortPeriodTerms);
638 }
639
640
641
642
643
644
645
646 private RealVector computeOsculatingElements(final RealVector filterCorrection,
647 final SpacecraftState meanState,
648 final RealVector shortPeriodTerms) {
649
650
651 final RealVector stateVector = toRealVector(meanState);
652
653
654 return stateVector.add(filterCorrection).add(shortPeriodTerms);
655
656 }
657
658
659
660
661
662 private RealVector toRealVector(final SpacecraftState state) {
663
664
665 final double[] stateArray = new double[6];
666 orbitType.mapOrbitToArray(state.getOrbit(), angleType, stateArray, null);
667
668
669 return new ArrayRealVector(stateArray);
670
671 }
672
673
674
675
676 public int getNumberSelectedOrbitalDrivers() {
677 return estimatedOrbitalParameters.getNbParams();
678 }
679
680
681
682
683 public int getNumberSelectedPropagationDrivers() {
684 return estimatedPropagationParameters.getNbParams();
685 }
686
687
688
689
690 public int getNumberSelectedMeasurementDrivers() {
691 return estimatedMeasurementsParameters.getNbParams();
692 }
693
694
695
696
697 private void updateParameters() {
698 final RealVector correctedState = correctedEstimate.getState();
699 int i = 0;
700 for (final DelegatingDriver driver : getEstimatedOrbitalParameters().getDrivers()) {
701
702 driver.setValue(driver.getValue() + correctedState.getEntry(i++));
703 }
704 for (final DelegatingDriver driver : getEstimatedPropagationParameters().getDrivers()) {
705
706 driver.setValue(driver.getValue() + correctedState.getEntry(i++));
707 }
708 for (final DelegatingDriver driver : getEstimatedMeasurementsParameters().getDrivers()) {
709
710 driver.setValue(driver.getValue() + correctedState.getEntry(i++));
711 }
712 }
713
714 }