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