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 java.util.ArrayList;
20 import java.util.Arrays;
21 import java.util.Comparator;
22 import java.util.HashMap;
23 import java.util.List;
24 import java.util.Map;
25
26 import org.hipparchus.filtering.kalman.ProcessEstimate;
27 import org.hipparchus.linear.ArrayRealVector;
28 import org.hipparchus.linear.MatrixUtils;
29 import org.hipparchus.linear.RealMatrix;
30 import org.hipparchus.linear.RealVector;
31 import org.orekit.errors.OrekitException;
32 import org.orekit.errors.OrekitMessages;
33 import org.orekit.estimation.measurements.EstimatedMeasurement;
34 import org.orekit.propagation.Propagator;
35 import org.orekit.propagation.SpacecraftState;
36 import org.orekit.propagation.conversion.PropagatorBuilder;
37 import org.orekit.time.AbsoluteDate;
38 import org.orekit.utils.ParameterDriver;
39 import org.orekit.utils.ParameterDriversList;
40 import org.orekit.utils.ParameterDriversList.DelegatingDriver;
41
42
43
44
45
46
47 abstract class AbstractKalmanEstimationCommon implements KalmanEstimation {
48
49
50 private final List<PropagatorBuilder> builders;
51
52
53 private final ParameterDriversList allEstimatedOrbitalParameters;
54
55
56 private final ParameterDriversList allEstimatedPropagationParameters;
57
58
59
60
61 private final ParameterDriversList[] estimatedOrbitalParameters;
62
63
64 private final ParameterDriversList[] estimatedPropagationParameters;
65
66
67 private final ParameterDriversList estimatedMeasurementsParameters;
68
69
70 private final int[] orbitsStartColumns;
71
72
73 private final int[] orbitsEndColumns;
74
75
76 private final Map<String, Integer> propagationParameterColumns;
77
78
79 private final Map<String, Integer> measurementParameterColumns;
80
81
82 private final List<CovarianceMatrixProvider> covarianceMatricesProviders;
83
84
85 private final CovarianceMatrixProvider measurementProcessNoiseMatrix;
86
87
88 private final int[][] covarianceIndirection;
89
90
91 private final double[] scale;
92
93
94 private ProcessEstimate correctedEstimate;
95
96
97 private int currentMeasurementNumber;
98
99
100 private final AbsoluteDate referenceDate;
101
102
103 private AbsoluteDate currentDate;
104
105
106 private final SpacecraftState[] predictedSpacecraftStates;
107
108
109 private final SpacecraftState[] correctedSpacecraftStates;
110
111
112 private EstimatedMeasurement<?> predictedMeasurement;
113
114
115 private EstimatedMeasurement<?> correctedMeasurement;
116
117
118
119
120
121
122
123 protected AbstractKalmanEstimationCommon(final List<PropagatorBuilder> propagatorBuilders,
124 final List<CovarianceMatrixProvider> covarianceMatricesProviders,
125 final ParameterDriversList estimatedMeasurementParameters,
126 final CovarianceMatrixProvider measurementProcessNoiseMatrix) {
127
128 this.builders = propagatorBuilders;
129 this.estimatedMeasurementsParameters = estimatedMeasurementParameters;
130 this.measurementParameterColumns = new HashMap<>(estimatedMeasurementsParameters.getDrivers().size());
131 this.currentMeasurementNumber = 0;
132 this.referenceDate = propagatorBuilders.getFirst().getOrbitalParameterFactory().getDate();
133 this.currentDate = referenceDate;
134
135 final Map<String, Integer> orbitalParameterColumns = new HashMap<>(6 * builders.size());
136 orbitsStartColumns = new int[builders.size()];
137 orbitsEndColumns = new int[builders.size()];
138 int columns = 0;
139 allEstimatedOrbitalParameters = new ParameterDriversList();
140 estimatedOrbitalParameters = new ParameterDriversList[builders.size()];
141 for (int k = 0; k < builders.size(); ++k) {
142 estimatedOrbitalParameters[k] = new ParameterDriversList();
143 orbitsStartColumns[k] = columns;
144 final String suffix = propagatorBuilders.size() > 1 ? "[" + k + "]" : null;
145 final ParameterDriversList drivers = builders.get(k).
146 getOrbitalParameterFactory().
147 getOrbitalParametersDrivers();
148 for (final ParameterDriver driver : drivers.getDrivers()) {
149 if (driver.getReferenceDate() == null) {
150 driver.setReferenceDate(currentDate);
151 }
152 if (suffix != null && !driver.getName().endsWith(suffix)) {
153
154
155 driver.setName(driver.getName() + suffix);
156 }
157 if (driver.isSelected()) {
158 allEstimatedOrbitalParameters.add(driver);
159 estimatedOrbitalParameters[k].add(driver);
160 orbitalParameterColumns.put(driver.getName(), columns++);
161 }
162 }
163 orbitsEndColumns[k] = columns;
164 }
165
166
167 allEstimatedPropagationParameters = new ParameterDriversList();
168 estimatedPropagationParameters = new ParameterDriversList[builders.size()];
169 final List<String> estimatedPropagationParametersNames = new ArrayList<>();
170 for (int k = 0; k < builders.size(); ++k) {
171 estimatedPropagationParameters[k] = new ParameterDriversList();
172 for (final ParameterDriver driver : builders.get(k).getPropagationParametersDrivers().getDrivers()) {
173 if (driver.getReferenceDate() == null) {
174 driver.setReferenceDate(currentDate);
175 }
176 if (driver.isSelected()) {
177 allEstimatedPropagationParameters.add(driver);
178 estimatedPropagationParameters[k].add(driver);
179 final String driverName = driver.getName();
180
181 if (!estimatedPropagationParametersNames.contains(driverName)) {
182 estimatedPropagationParametersNames.add(driverName);
183 }
184 }
185 }
186 }
187 estimatedPropagationParametersNames.sort(Comparator.naturalOrder());
188
189
190 propagationParameterColumns = new HashMap<>(estimatedPropagationParametersNames.size());
191 for (final String driverName : estimatedPropagationParametersNames) {
192 propagationParameterColumns.put(driverName, columns);
193 ++columns;
194 }
195
196
197 for (final ParameterDriver parameter : estimatedMeasurementsParameters.getDrivers()) {
198 if (parameter.getReferenceDate() == null) {
199 parameter.setReferenceDate(currentDate);
200 }
201 measurementParameterColumns.put(parameter.getName(), columns);
202 ++columns;
203 }
204
205
206 this.covarianceMatricesProviders = covarianceMatricesProviders;
207 this.measurementProcessNoiseMatrix = measurementProcessNoiseMatrix;
208 this.covarianceIndirection = new int[builders.size()][columns];
209 for (int k = 0; k < covarianceIndirection.length; ++k) {
210 final ParameterDriversList orbitDrivers = builders.get(k).
211 getOrbitalParameterFactory().
212 getOrbitalParametersDrivers();
213 final ParameterDriversList parametersDrivers = builders.get(k).
214 getPropagationParametersDrivers();
215 Arrays.fill(covarianceIndirection[k], -1);
216 int i = 0;
217 for (final ParameterDriver driver : orbitDrivers.getDrivers()) {
218 final Integer c = orbitalParameterColumns.get(driver.getName());
219 if (c != null) {
220 covarianceIndirection[k][i++] = c;
221 }
222 }
223 for (final ParameterDriver driver : parametersDrivers.getDrivers()) {
224 final Integer c = propagationParameterColumns.get(driver.getName());
225 if (c != null) {
226 covarianceIndirection[k][i++] = c;
227 }
228 }
229 for (final ParameterDriver driver : estimatedMeasurementParameters.getDrivers()) {
230 final Integer c = measurementParameterColumns.get(driver.getName());
231 if (c != null) {
232 covarianceIndirection[k][i++] = c;
233 }
234 }
235 }
236
237
238 this.scale = new double[columns];
239 int index = 0;
240 for (final ParameterDriver driver : allEstimatedOrbitalParameters.getDrivers()) {
241 scale[index++] = driver.getScale();
242 }
243 for (final ParameterDriver driver : allEstimatedPropagationParameters.getDrivers()) {
244 scale[index++] = driver.getScale();
245 }
246 for (final ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
247 scale[index++] = driver.getScale();
248 }
249
250
251 this.predictedSpacecraftStates = new SpacecraftState[builders.size()];
252 for (int i = 0; i < builders.size(); ++i) {
253 predictedSpacecraftStates[i] = builders.get(i).buildPropagator().getInitialState();
254 }
255 this.correctedSpacecraftStates = predictedSpacecraftStates.clone();
256
257
258 final RealVector correctedState = MatrixUtils.createRealVector(columns);
259
260 int p = 0;
261 for (final ParameterDriver driver : allEstimatedOrbitalParameters.getDrivers()) {
262 correctedState.setEntry(p++, driver.getNormalizedValue());
263 }
264 for (final ParameterDriver driver : allEstimatedPropagationParameters.getDrivers()) {
265 correctedState.setEntry(p++, driver.getNormalizedValue());
266 }
267 for (final ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
268 correctedState.setEntry(p++, driver.getNormalizedValue());
269 }
270
271
272 final RealMatrix physicalProcessNoise = MatrixUtils.createRealMatrix(columns, columns);
273 for (int k = 0; k < covarianceMatricesProviders.size(); ++k) {
274
275
276 final int nbMeas = estimatedMeasurementParameters.getNbParams();
277
278
279 final int nbDyn = orbitsEndColumns[k] - orbitsStartColumns[k] +
280 estimatedPropagationParameters[k].getNbParams();
281
282
283 final RealMatrix noiseK = MatrixUtils.createRealMatrix(nbDyn + nbMeas, nbDyn + nbMeas);
284 if (nbDyn > 0) {
285 final RealMatrix noiseP = covarianceMatricesProviders.get(k).
286 getInitialCovarianceMatrix(correctedSpacecraftStates[k]);
287 if (measurementProcessNoiseMatrix == null && noiseP.getRowDimension() != nbDyn + nbMeas) {
288 throw new OrekitException(OrekitMessages.WRONG_PROCESS_COVARIANCE_DIMENSION,
289 nbDyn + nbMeas, noiseP.getRowDimension());
290 } else if (measurementProcessNoiseMatrix != null && noiseP.getRowDimension() != nbDyn) {
291 throw new OrekitException(OrekitMessages.WRONG_PROCESS_COVARIANCE_DIMENSION,
292 nbDyn, noiseP.getRowDimension());
293 }
294 noiseK.setSubMatrix(noiseP.getData(), 0, 0);
295 }
296 if (measurementProcessNoiseMatrix != null) {
297 final RealMatrix noiseM = measurementProcessNoiseMatrix.
298 getInitialCovarianceMatrix(correctedSpacecraftStates[k]);
299 if (noiseM.getRowDimension() != nbMeas) {
300 throw new OrekitException(OrekitMessages.WRONG_MEASUREMENT_COVARIANCE_DIMENSION,
301 nbMeas, noiseM.getRowDimension());
302 }
303 noiseK.setSubMatrix(noiseM.getData(), nbDyn, nbDyn);
304 }
305
306 KalmanEstimatorUtil.checkDimension(noiseK.getRowDimension(),
307 builders.get(k).getOrbitalParameterFactory().getOrbitalParametersDrivers(),
308 builders.get(k).getPropagationParametersDrivers(),
309 estimatedMeasurementsParameters);
310
311 final int[] indK = covarianceIndirection[k];
312 for (int i = 0; i < indK.length; ++i) {
313 if (indK[i] >= 0) {
314 for (int j = 0; j < indK.length; ++j) {
315 if (indK[j] >= 0) {
316 physicalProcessNoise.setEntry(indK[i], indK[j], noiseK.getEntry(i, j));
317 }
318 }
319 }
320 }
321
322 }
323 final RealMatrix correctedCovariance = KalmanEstimatorUtil.normalizeCovarianceMatrix(physicalProcessNoise, scale);
324
325 correctedEstimate = new ProcessEstimate(0.0, correctedState, correctedCovariance);
326 }
327
328
329
330 @Override
331 public RealMatrix getPhysicalStateTransitionMatrix() {
332
333
334
335
336 return correctedEstimate.getStateTransitionMatrix() == null ?
337 null : KalmanEstimatorUtil.unnormalizeStateTransitionMatrix(correctedEstimate.getStateTransitionMatrix(), scale);
338 }
339
340
341 @Override
342 public RealMatrix getPhysicalMeasurementJacobian() {
343
344
345
346
347
348
349 return correctedEstimate.getMeasurementJacobian() == null ?
350 null : KalmanEstimatorUtil.unnormalizeMeasurementJacobian(correctedEstimate.getMeasurementJacobian(),
351 scale,
352 correctedMeasurement.getObservedMeasurement().getTheoreticalStandardDeviation());
353 }
354
355
356 @Override
357 public RealMatrix getPhysicalInnovationCovarianceMatrix() {
358
359
360
361
362 return correctedEstimate.getInnovationCovariance() == null ?
363 null : KalmanEstimatorUtil.unnormalizeInnovationCovarianceMatrix(correctedEstimate.getInnovationCovariance(),
364 predictedMeasurement.getObservedMeasurement().getTheoreticalStandardDeviation());
365 }
366
367
368 @Override
369 public RealMatrix getPhysicalKalmanGain() {
370
371
372
373
374
375
376 return correctedEstimate.getKalmanGain() == null ?
377 null : KalmanEstimatorUtil.unnormalizeKalmanGainMatrix(correctedEstimate.getKalmanGain(),
378 scale,
379 correctedMeasurement.getObservedMeasurement().getTheoreticalStandardDeviation());
380 }
381
382
383 @Override
384 public SpacecraftState[] getPredictedSpacecraftStates() {
385 return predictedSpacecraftStates.clone();
386 }
387
388
389 @Override
390 public SpacecraftState[] getCorrectedSpacecraftStates() {
391 return correctedSpacecraftStates.clone();
392 }
393
394
395 @Override
396 public int getCurrentMeasurementNumber() {
397 return currentMeasurementNumber;
398 }
399
400
401 @Override
402 public AbsoluteDate getCurrentDate() {
403 return currentDate;
404 }
405
406
407 @Override
408 public EstimatedMeasurement<?> getPredictedMeasurement() {
409 return predictedMeasurement;
410 }
411
412
413 @Override
414 public EstimatedMeasurement<?> getCorrectedMeasurement() {
415 return correctedMeasurement;
416 }
417
418
419 @Override
420 public RealVector getPhysicalEstimatedState() {
421
422
423
424 final RealVector physicalEstimatedState = new ArrayRealVector(scale.length);
425 int i = 0;
426 for (final DelegatingDriver driver : getEstimatedOrbitalParameters().getDrivers()) {
427 physicalEstimatedState.setEntry(i++, driver.getValue());
428 }
429 for (final DelegatingDriver driver : getEstimatedPropagationParameters().getDrivers()) {
430 physicalEstimatedState.setEntry(i++, driver.getValue());
431 }
432 for (final DelegatingDriver driver : getEstimatedMeasurementsParameters().getDrivers()) {
433 physicalEstimatedState.setEntry(i++, driver.getValue());
434 }
435
436 return physicalEstimatedState;
437 }
438
439
440 @Override
441 public RealMatrix getPhysicalEstimatedCovarianceMatrix() {
442
443
444
445
446
447 return KalmanEstimatorUtil.unnormalizeCovarianceMatrix(correctedEstimate.getCovariance(), scale);
448 }
449
450
451 @Override
452 public ParameterDriversList getEstimatedOrbitalParameters() {
453 return allEstimatedOrbitalParameters;
454 }
455
456
457 @Override
458 public ParameterDriversList getEstimatedPropagationParameters() {
459 return allEstimatedPropagationParameters;
460 }
461
462
463 @Override
464 public ParameterDriversList getEstimatedMeasurementsParameters() {
465 return estimatedMeasurementsParameters;
466 }
467
468
469
470
471 public ProcessEstimate getEstimate() {
472 return correctedEstimate;
473 }
474
475
476
477
478 public List<PropagatorBuilder> getBuilders() {
479 return builders;
480 }
481
482
483
484
485 public Propagator[] getEstimatedPropagators() {
486
487 final Propagator[] propagators = new Propagator[getBuilders().size()];
488 for (int k = 0; k < getBuilders().size(); ++k) {
489 propagators[k] = getBuilders().get(k).buildPropagator();
490 }
491 return propagators;
492 }
493
494
495
496
497
498
499 protected RealMatrix getNormalizedProcessNoise(final int stateDimension) {
500 final RealMatrix physicalProcessNoise = MatrixUtils.createRealMatrix(stateDimension, stateDimension);
501 for (int k = 0; k < covarianceMatricesProviders.size(); ++k) {
502
503
504 final int nbMeas = estimatedMeasurementsParameters.getNbParams();
505
506
507 final int nbDyn = orbitsEndColumns[k] - orbitsStartColumns[k] +
508 estimatedPropagationParameters[k].getNbParams();
509
510
511 final RealMatrix noiseK = MatrixUtils.createRealMatrix(nbDyn + nbMeas, nbDyn + nbMeas);
512 if (nbDyn > 0) {
513 final RealMatrix noiseP = covarianceMatricesProviders.get(k).
514 getProcessNoiseMatrix(correctedSpacecraftStates[k],
515 predictedSpacecraftStates[k]);
516 if (measurementProcessNoiseMatrix == null && noiseP.getRowDimension() != nbDyn + nbMeas) {
517 throw new OrekitException(OrekitMessages.WRONG_PROCESS_COVARIANCE_DIMENSION,
518 nbDyn + nbMeas, noiseP.getRowDimension());
519 } else if (measurementProcessNoiseMatrix != null && noiseP.getRowDimension() != nbDyn) {
520 throw new OrekitException(OrekitMessages.WRONG_PROCESS_COVARIANCE_DIMENSION,
521 nbDyn, noiseP.getRowDimension());
522 }
523 noiseK.setSubMatrix(noiseP.getData(), 0, 0);
524 }
525 if (measurementProcessNoiseMatrix != null) {
526 final RealMatrix noiseM = measurementProcessNoiseMatrix.
527 getProcessNoiseMatrix(correctedSpacecraftStates[k],
528 predictedSpacecraftStates[k]);
529 if (noiseM.getRowDimension() != nbMeas) {
530 throw new OrekitException(OrekitMessages.WRONG_MEASUREMENT_COVARIANCE_DIMENSION,
531 nbMeas, noiseM.getRowDimension());
532 }
533 noiseK.setSubMatrix(noiseM.getData(), nbDyn, nbDyn);
534 }
535
536 KalmanEstimatorUtil.checkDimension(noiseK.getRowDimension(),
537 builders.get(k).getOrbitalParameterFactory().getOrbitalParametersDrivers(),
538 builders.get(k).getPropagationParametersDrivers(),
539 estimatedMeasurementsParameters);
540
541 final int[] indK = covarianceIndirection[k];
542 for (int i = 0; i < indK.length; ++i) {
543 if (indK[i] >= 0) {
544 for (int j = 0; j < indK.length; ++j) {
545 if (indK[j] >= 0) {
546 physicalProcessNoise.setEntry(indK[i], indK[j], noiseK.getEntry(i, j));
547 }
548 }
549 }
550 }
551
552 }
553 return KalmanEstimatorUtil.normalizeCovarianceMatrix(physicalProcessNoise, scale);
554 }
555
556
557
558
559 protected int[] getOrbitsStartColumns() {
560 return orbitsStartColumns;
561 }
562
563
564
565
566 protected Map<String, Integer> getPropagationParameterColumns() {
567 return propagationParameterColumns;
568 }
569
570
571
572
573 protected Map<String, Integer> getMeasurementParameterColumns() {
574 return measurementParameterColumns;
575 }
576
577
578
579
580 protected ParameterDriversList[] getEstimatedPropagationParametersArray() {
581 return estimatedPropagationParameters;
582 }
583
584
585
586
587 protected ParameterDriversList[] getEstimatedOrbitalParametersArray() {
588 return estimatedOrbitalParameters;
589 }
590
591
592
593
594 protected int[][] getCovarianceIndirection() {
595 return covarianceIndirection;
596 }
597
598
599
600
601 protected double[] getScale() {
602 return scale;
603 }
604
605
606
607
608 protected ProcessEstimate getCorrectedEstimate() {
609 return correctedEstimate;
610 }
611
612
613
614
615 protected void setCorrectedEstimate(final ProcessEstimate correctedEstimate) {
616 this.correctedEstimate = correctedEstimate;
617 }
618
619
620
621
622 protected AbsoluteDate getReferenceDate() {
623 return referenceDate;
624 }
625
626
627 protected void incrementCurrentMeasurementNumber() {
628 currentMeasurementNumber += 1;
629 }
630
631
632
633
634 protected void setCurrentDate(final AbsoluteDate currentDate) {
635 this.currentDate = currentDate;
636 }
637
638
639
640
641
642
643 protected void setCorrectedSpacecraftState(final SpacecraftState correctedSpacecraftState, final int index) {
644 this.correctedSpacecraftStates[index] = correctedSpacecraftState;
645 }
646
647
648
649
650
651
652 protected void setPredictedSpacecraftState(final SpacecraftState predictedSpacecraftState, final int index) {
653 this.predictedSpacecraftStates[index] = predictedSpacecraftState;
654 }
655
656
657
658
659 protected void setPredictedMeasurement(final EstimatedMeasurement<?> predictedMeasurement) {
660 this.predictedMeasurement = predictedMeasurement;
661 }
662
663
664
665
666 protected void setCorrectedMeasurement(final EstimatedMeasurement<?> correctedMeasurement) {
667 this.correctedMeasurement = correctedMeasurement;
668 }
669 }