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