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.filtering.kalman.extended.NonLinearEvolution;
28 import org.hipparchus.filtering.kalman.extended.NonLinearProcess;
29 import org.hipparchus.linear.Array2DRowRealMatrix;
30 import org.hipparchus.linear.ArrayRealVector;
31 import org.hipparchus.linear.MatrixUtils;
32 import org.hipparchus.linear.RealMatrix;
33 import org.hipparchus.linear.RealVector;
34 import org.hipparchus.util.FastMath;
35 import org.orekit.errors.OrekitException;
36 import org.orekit.errors.OrekitMessages;
37 import org.orekit.estimation.measurements.EstimatedMeasurement;
38 import org.orekit.estimation.measurements.EstimationModifier;
39 import org.orekit.estimation.measurements.ObservableSatellite;
40 import org.orekit.estimation.measurements.ObservedMeasurement;
41 import org.orekit.estimation.measurements.modifiers.DynamicOutlierFilter;
42 import org.orekit.orbits.Orbit;
43 import org.orekit.propagation.MatricesHarvester;
44 import org.orekit.propagation.PropagationType;
45 import org.orekit.propagation.Propagator;
46 import org.orekit.propagation.SpacecraftState;
47 import org.orekit.propagation.conversion.OrbitDeterminationPropagatorBuilder;
48 import org.orekit.propagation.integration.AbstractJacobiansMapper;
49 import org.orekit.time.AbsoluteDate;
50 import org.orekit.utils.ParameterDriver;
51 import org.orekit.utils.ParameterDriversList;
52 import org.orekit.utils.ParameterDriversList.DelegatingDriver;
53
54
55
56
57
58
59
60
61 public abstract class AbstractKalmanModel implements KalmanEstimation, NonLinearProcess<MeasurementDecorator> {
62
63
64 private final List<OrbitDeterminationPropagatorBuilder> builders;
65
66
67 private final ParameterDriversList allEstimatedOrbitalParameters;
68
69
70 private final ParameterDriversList allEstimatedPropagationParameters;
71
72
73
74
75 private final ParameterDriversList[] estimatedOrbitalParameters;
76
77
78 private final ParameterDriversList[] estimatedPropagationParameters;
79
80
81 private final ParameterDriversList estimatedMeasurementsParameters;
82
83
84 private final int[] orbitsStartColumns;
85
86
87 private final int[] orbitsEndColumns;
88
89
90 private final Map<String, Integer> propagationParameterColumns;
91
92
93 private final Map<String, Integer> measurementParameterColumns;
94
95
96 private final List<CovarianceMatrixProvider> covarianceMatricesProviders;
97
98
99 private final CovarianceMatrixProvider measurementProcessNoiseMatrix;
100
101
102 private final int[][] covarianceIndirection;
103
104
105 private final double[] scale;
106
107
108 private MatricesHarvester[] harvesters;
109
110
111 private Propagator[] referenceTrajectories;
112
113
114 private ProcessEstimate correctedEstimate;
115
116
117 private int currentMeasurementNumber;
118
119
120 private AbsoluteDate referenceDate;
121
122
123 private AbsoluteDate currentDate;
124
125
126 private SpacecraftState[] predictedSpacecraftStates;
127
128
129 private SpacecraftState[] correctedSpacecraftStates;
130
131
132 private EstimatedMeasurement<?> predictedMeasurement;
133
134
135 private EstimatedMeasurement<?> correctedMeasurement;
136
137
138 private PropagationType propagationType;
139
140
141 private PropagationType stateType;
142
143
144
145
146
147
148
149
150
151
152 protected AbstractKalmanModel(final List<OrbitDeterminationPropagatorBuilder> propagatorBuilders,
153 final List<CovarianceMatrixProvider> covarianceMatricesProviders,
154 final ParameterDriversList estimatedMeasurementParameters,
155 final CovarianceMatrixProvider measurementProcessNoiseMatrix,
156 final MatricesHarvester[] harvesters) {
157 this(propagatorBuilders, covarianceMatricesProviders, estimatedMeasurementParameters,
158 measurementProcessNoiseMatrix, harvesters, PropagationType.MEAN, PropagationType.MEAN);
159 }
160
161
162
163
164
165
166
167
168
169
170
171
172 protected AbstractKalmanModel(final List<OrbitDeterminationPropagatorBuilder> propagatorBuilders,
173 final List<CovarianceMatrixProvider> covarianceMatricesProviders,
174 final ParameterDriversList estimatedMeasurementParameters,
175 final CovarianceMatrixProvider measurementProcessNoiseMatrix,
176 final MatricesHarvester[] harvesters,
177 final PropagationType propagationType,
178 final PropagationType stateType) {
179
180 this.builders = propagatorBuilders;
181 this.estimatedMeasurementsParameters = estimatedMeasurementParameters;
182 this.measurementParameterColumns = new HashMap<>(estimatedMeasurementsParameters.getDrivers().size());
183 this.currentMeasurementNumber = 0;
184 this.referenceDate = propagatorBuilders.get(0).getInitialOrbitDate();
185 this.currentDate = referenceDate;
186 this.propagationType = propagationType;
187 this.stateType = stateType;
188
189 final Map<String, Integer> orbitalParameterColumns = new HashMap<>(6 * builders.size());
190 orbitsStartColumns = new int[builders.size()];
191 orbitsEndColumns = new int[builders.size()];
192 int columns = 0;
193 allEstimatedOrbitalParameters = new ParameterDriversList();
194 estimatedOrbitalParameters = new ParameterDriversList[builders.size()];
195 for (int k = 0; k < builders.size(); ++k) {
196 estimatedOrbitalParameters[k] = new ParameterDriversList();
197 orbitsStartColumns[k] = columns;
198 final String suffix = propagatorBuilders.size() > 1 ? "[" + k + "]" : null;
199 for (final ParameterDriver driver : builders.get(k).getOrbitalParametersDrivers().getDrivers()) {
200 if (driver.getReferenceDate() == null) {
201 driver.setReferenceDate(currentDate);
202 }
203 if (suffix != null && !driver.getName().endsWith(suffix)) {
204
205
206 driver.setName(driver.getName() + suffix);
207 }
208 if (driver.isSelected()) {
209 allEstimatedOrbitalParameters.add(driver);
210 estimatedOrbitalParameters[k].add(driver);
211 orbitalParameterColumns.put(driver.getName(), columns++);
212 }
213 }
214 orbitsEndColumns[k] = columns;
215 }
216
217
218 allEstimatedPropagationParameters = new ParameterDriversList();
219 estimatedPropagationParameters = new ParameterDriversList[builders.size()];
220 final List<String> estimatedPropagationParametersNames = new ArrayList<>();
221 for (int k = 0; k < builders.size(); ++k) {
222 estimatedPropagationParameters[k] = new ParameterDriversList();
223 for (final ParameterDriver driver : builders.get(k).getPropagationParametersDrivers().getDrivers()) {
224 if (driver.getReferenceDate() == null) {
225 driver.setReferenceDate(currentDate);
226 }
227 if (driver.isSelected()) {
228 allEstimatedPropagationParameters.add(driver);
229 estimatedPropagationParameters[k].add(driver);
230 final String driverName = driver.getName();
231
232 if (!estimatedPropagationParametersNames.contains(driverName)) {
233 estimatedPropagationParametersNames.add(driverName);
234 }
235 }
236 }
237 }
238 estimatedPropagationParametersNames.sort(Comparator.naturalOrder());
239
240
241 propagationParameterColumns = new HashMap<>(estimatedPropagationParametersNames.size());
242 for (final String driverName : estimatedPropagationParametersNames) {
243 propagationParameterColumns.put(driverName, columns);
244 ++columns;
245 }
246
247
248 for (final ParameterDriver parameter : estimatedMeasurementsParameters.getDrivers()) {
249 if (parameter.getReferenceDate() == null) {
250 parameter.setReferenceDate(currentDate);
251 }
252 measurementParameterColumns.put(parameter.getName(), columns);
253 ++columns;
254 }
255
256
257 this.covarianceMatricesProviders = covarianceMatricesProviders;
258 this.measurementProcessNoiseMatrix = measurementProcessNoiseMatrix;
259 this.covarianceIndirection = new int[covarianceMatricesProviders.size()][columns];
260 for (int k = 0; k < covarianceIndirection.length; ++k) {
261 final ParameterDriversList orbitDrivers = builders.get(k).getOrbitalParametersDrivers();
262 final ParameterDriversList parametersDrivers = builders.get(k).getPropagationParametersDrivers();
263 Arrays.fill(covarianceIndirection[k], -1);
264 int i = 0;
265 for (final ParameterDriver driver : orbitDrivers.getDrivers()) {
266 final Integer c = orbitalParameterColumns.get(driver.getName());
267 covarianceIndirection[k][i++] = (c == null) ? -1 : c.intValue();
268 }
269 for (final ParameterDriver driver : parametersDrivers.getDrivers()) {
270 final Integer c = propagationParameterColumns.get(driver.getName());
271 if (c != null) {
272 covarianceIndirection[k][i++] = c.intValue();
273 }
274 }
275 for (final ParameterDriver driver : estimatedMeasurementParameters.getDrivers()) {
276 final Integer c = measurementParameterColumns.get(driver.getName());
277 if (c != null) {
278 covarianceIndirection[k][i++] = c.intValue();
279 }
280 }
281 }
282
283
284 this.scale = new double[columns];
285 int index = 0;
286 for (final ParameterDriver driver : allEstimatedOrbitalParameters.getDrivers()) {
287 scale[index++] = driver.getScale();
288 }
289 for (final ParameterDriver driver : allEstimatedPropagationParameters.getDrivers()) {
290 scale[index++] = driver.getScale();
291 }
292 for (final ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
293 scale[index++] = driver.getScale();
294 }
295
296
297 this.harvesters = harvesters.clone();
298 updateReferenceTrajectories(getEstimatedPropagators(), propagationType, stateType);
299 this.predictedSpacecraftStates = new SpacecraftState[referenceTrajectories.length];
300 for (int i = 0; i < predictedSpacecraftStates.length; ++i) {
301 predictedSpacecraftStates[i] = referenceTrajectories[i].getInitialState();
302 };
303 this.correctedSpacecraftStates = predictedSpacecraftStates.clone();
304
305
306 final RealVector correctedState = MatrixUtils.createRealVector(columns);
307
308 int p = 0;
309 for (final ParameterDriver driver : allEstimatedOrbitalParameters.getDrivers()) {
310 correctedState.setEntry(p++, driver.getNormalizedValue());
311 }
312 for (final ParameterDriver driver : allEstimatedPropagationParameters.getDrivers()) {
313 correctedState.setEntry(p++, driver.getNormalizedValue());
314 }
315 for (final ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
316 correctedState.setEntry(p++, driver.getNormalizedValue());
317 }
318
319
320 final RealMatrix physicalProcessNoise = MatrixUtils.createRealMatrix(columns, columns);
321 for (int k = 0; k < covarianceMatricesProviders.size(); ++k) {
322
323
324 final int nbMeas = estimatedMeasurementParameters.getNbParams();
325
326
327 final int nbDyn = orbitsEndColumns[k] - orbitsStartColumns[k] +
328 estimatedPropagationParameters[k].getNbParams();
329
330
331 final RealMatrix noiseK = MatrixUtils.createRealMatrix(nbDyn + nbMeas, nbDyn + nbMeas);
332 final RealMatrix noiseP = covarianceMatricesProviders.get(k).
333 getInitialCovarianceMatrix(correctedSpacecraftStates[k]);
334 noiseK.setSubMatrix(noiseP.getData(), 0, 0);
335 if (measurementProcessNoiseMatrix != null) {
336 final RealMatrix noiseM = measurementProcessNoiseMatrix.
337 getInitialCovarianceMatrix(correctedSpacecraftStates[k]);
338 noiseK.setSubMatrix(noiseM.getData(), nbDyn, nbDyn);
339 }
340
341 checkDimension(noiseK.getRowDimension(),
342 builders.get(k).getOrbitalParametersDrivers(),
343 builders.get(k).getPropagationParametersDrivers(),
344 estimatedMeasurementsParameters);
345
346 final int[] indK = covarianceIndirection[k];
347 for (int i = 0; i < indK.length; ++i) {
348 if (indK[i] >= 0) {
349 for (int j = 0; j < indK.length; ++j) {
350 if (indK[j] >= 0) {
351 physicalProcessNoise.setEntry(indK[i], indK[j], noiseK.getEntry(i, j));
352 }
353 }
354 }
355 }
356
357 }
358 final RealMatrix correctedCovariance = normalizeCovarianceMatrix(physicalProcessNoise);
359
360 correctedEstimate = new ProcessEstimate(0.0, correctedState, correctedCovariance);
361
362 }
363
364
365
366
367
368
369 protected abstract void updateReferenceTrajectories(Propagator[] propagators,
370 PropagationType pType,
371 PropagationType sType);
372
373
374
375
376
377
378 protected void analyticalDerivativeComputations(final AbstractJacobiansMapper mapper, final SpacecraftState state) {
379
380 }
381
382
383
384
385
386
387
388 private void checkDimension(final int dimension,
389 final ParameterDriversList orbitalParameters,
390 final ParameterDriversList propagationParameters,
391 final ParameterDriversList measurementParameters) {
392
393
394
395 int requiredDimension = orbitalParameters.getNbParams();
396 for (final ParameterDriver driver : propagationParameters.getDrivers()) {
397 if (driver.isSelected()) {
398 ++requiredDimension;
399 }
400 }
401 for (final ParameterDriver driver : measurementParameters.getDrivers()) {
402 if (driver.isSelected()) {
403 ++requiredDimension;
404 }
405 }
406
407 if (dimension != requiredDimension) {
408
409 final StringBuilder builder = new StringBuilder();
410 for (final ParameterDriver driver : orbitalParameters.getDrivers()) {
411 if (builder.length() > 0) {
412 builder.append(", ");
413 }
414 builder.append(driver.getName());
415 }
416 for (final ParameterDriver driver : propagationParameters.getDrivers()) {
417 if (driver.isSelected()) {
418 builder.append(driver.getName());
419 }
420 }
421 for (final ParameterDriver driver : measurementParameters.getDrivers()) {
422 if (driver.isSelected()) {
423 builder.append(driver.getName());
424 }
425 }
426 throw new OrekitException(OrekitMessages.DIMENSION_INCONSISTENT_WITH_PARAMETERS,
427 dimension, builder.toString());
428 }
429
430 }
431
432
433 @Override
434 public RealMatrix getPhysicalStateTransitionMatrix() {
435
436
437
438
439
440
441 final RealMatrix normalizedSTM = correctedEstimate.getStateTransitionMatrix();
442
443 if (normalizedSTM == null) {
444 return null;
445 } else {
446
447 final int nbParams = normalizedSTM.getRowDimension();
448 final RealMatrix physicalSTM = MatrixUtils.createRealMatrix(nbParams, nbParams);
449
450
451 for (int i = 0; i < nbParams; ++i) {
452 for (int j = 0; j < nbParams; ++j) {
453 physicalSTM.setEntry(i, j,
454 normalizedSTM.getEntry(i, j) * scale[i] / scale[j]);
455 }
456 }
457 return physicalSTM;
458 }
459 }
460
461
462 @Override
463 public RealMatrix getPhysicalMeasurementJacobian() {
464
465
466
467
468
469
470
471
472 final RealMatrix normalizedH = correctedEstimate.getMeasurementJacobian();
473
474 if (normalizedH == null) {
475 return null;
476 } else {
477
478 final double[] sigmas = correctedMeasurement.getObservedMeasurement().getTheoreticalStandardDeviation();
479
480
481 final int nbLine = normalizedH.getRowDimension();
482 final int nbCol = normalizedH.getColumnDimension();
483 final RealMatrix physicalH = MatrixUtils.createRealMatrix(nbLine, nbCol);
484
485
486 for (int i = 0; i < nbLine; ++i) {
487 for (int j = 0; j < nbCol; ++j) {
488 physicalH.setEntry(i, j, normalizedH.getEntry(i, j) * sigmas[i] / scale[j]);
489 }
490 }
491 return physicalH;
492 }
493 }
494
495
496 @Override
497 public RealMatrix getPhysicalInnovationCovarianceMatrix() {
498
499
500
501
502
503
504 final RealMatrix normalizedS = correctedEstimate.getInnovationCovariance();
505
506 if (normalizedS == null) {
507 return null;
508 } else {
509
510 final double[] sigmas = correctedMeasurement.getObservedMeasurement().getTheoreticalStandardDeviation();
511
512
513 final int nbMeas = sigmas.length;
514 final RealMatrix physicalS = MatrixUtils.createRealMatrix(nbMeas, nbMeas);
515
516
517 for (int i = 0; i < nbMeas; ++i) {
518 for (int j = 0; j < nbMeas; ++j) {
519 physicalS.setEntry(i, j, normalizedS.getEntry(i, j) * sigmas[i] * sigmas[j]);
520 }
521 }
522 return physicalS;
523 }
524 }
525
526
527 @Override
528 public RealMatrix getPhysicalKalmanGain() {
529
530
531
532
533
534
535
536
537 final RealMatrix normalizedK = correctedEstimate.getKalmanGain();
538
539 if (normalizedK == null) {
540 return null;
541 } else {
542
543 final double[] sigmas = correctedMeasurement.getObservedMeasurement().getTheoreticalStandardDeviation();
544
545
546 final int nbLine = normalizedK.getRowDimension();
547 final int nbCol = normalizedK.getColumnDimension();
548 final RealMatrix physicalK = MatrixUtils.createRealMatrix(nbLine, nbCol);
549
550
551 for (int i = 0; i < nbLine; ++i) {
552 for (int j = 0; j < nbCol; ++j) {
553 physicalK.setEntry(i, j, normalizedK.getEntry(i, j) * scale[i] / sigmas[j]);
554 }
555 }
556 return physicalK;
557 }
558 }
559
560
561 @Override
562 public SpacecraftState[] getPredictedSpacecraftStates() {
563 return predictedSpacecraftStates.clone();
564 }
565
566
567 @Override
568 public SpacecraftState[] getCorrectedSpacecraftStates() {
569 return correctedSpacecraftStates.clone();
570 }
571
572
573 @Override
574 public int getCurrentMeasurementNumber() {
575 return currentMeasurementNumber;
576 }
577
578
579 @Override
580 public AbsoluteDate getCurrentDate() {
581 return currentDate;
582 }
583
584
585 @Override
586 public EstimatedMeasurement<?> getPredictedMeasurement() {
587 return predictedMeasurement;
588 }
589
590
591 @Override
592 public EstimatedMeasurement<?> getCorrectedMeasurement() {
593 return correctedMeasurement;
594 }
595
596
597 @Override
598 public RealVector getPhysicalEstimatedState() {
599
600
601
602 final RealVector physicalEstimatedState = new ArrayRealVector(scale.length);
603 int i = 0;
604 for (final DelegatingDriver driver : getEstimatedOrbitalParameters().getDrivers()) {
605 physicalEstimatedState.setEntry(i++, driver.getValue());
606 }
607 for (final DelegatingDriver driver : getEstimatedPropagationParameters().getDrivers()) {
608 physicalEstimatedState.setEntry(i++, driver.getValue());
609 }
610 for (final DelegatingDriver driver : getEstimatedMeasurementsParameters().getDrivers()) {
611 physicalEstimatedState.setEntry(i++, driver.getValue());
612 }
613
614 return physicalEstimatedState;
615 }
616
617
618 @Override
619 public RealMatrix getPhysicalEstimatedCovarianceMatrix() {
620
621
622
623
624
625
626
627 final RealMatrix normalizedP = correctedEstimate.getCovariance();
628
629
630 final int nbParams = normalizedP.getRowDimension();
631 final RealMatrix physicalP = MatrixUtils.createRealMatrix(nbParams, nbParams);
632
633
634 for (int i = 0; i < nbParams; ++i) {
635 for (int j = 0; j < nbParams; ++j) {
636 physicalP.setEntry(i, j, normalizedP.getEntry(i, j) * scale[i] * scale[j]);
637 }
638 }
639 return physicalP;
640 }
641
642
643 @Override
644 public ParameterDriversList getEstimatedOrbitalParameters() {
645 return allEstimatedOrbitalParameters;
646 }
647
648
649 @Override
650 public ParameterDriversList getEstimatedPropagationParameters() {
651 return allEstimatedPropagationParameters;
652 }
653
654
655 @Override
656 public ParameterDriversList getEstimatedMeasurementsParameters() {
657 return estimatedMeasurementsParameters;
658 }
659
660
661
662
663 public ProcessEstimate getEstimate() {
664 return correctedEstimate;
665 }
666
667
668
669
670
671
672
673 private RealMatrix getErrorStateTransitionMatrix() {
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697 final RealMatrix stm = MatrixUtils.createRealIdentityMatrix(correctedEstimate.getState().getDimension());
698
699
700 for (int k = 0; k < predictedSpacecraftStates.length; ++k) {
701
702
703 final int[] indK = covarianceIndirection[k];
704
705
706 harvesters[k].setReferenceState(predictedSpacecraftStates[k]);
707
708
709 final int nbOrbParams = estimatedOrbitalParameters[k].getNbParams();
710 if (nbOrbParams > 0) {
711 final RealMatrix dYdY0 = harvesters[k].getStateTransitionMatrix(predictedSpacecraftStates[k]);
712
713
714 for (int i = 0; i < dYdY0.getRowDimension(); ++i) {
715 for (int j = 0; j < nbOrbParams; ++j) {
716 stm.setEntry(indK[i], indK[j], dYdY0.getEntry(i, j));
717 }
718 }
719 }
720
721
722 final int nbParams = estimatedPropagationParameters[k].getNbParams();
723 if (nbParams > 0) {
724 final RealMatrix dYdPp = harvesters[k].getParametersJacobian(predictedSpacecraftStates[k]);
725
726
727 for (int i = 0; i < dYdPp.getRowDimension(); ++i) {
728 for (int j = 0; j < nbParams; ++j) {
729 stm.setEntry(indK[i], indK[j + 6], dYdPp.getEntry(i, j));
730 }
731 }
732
733 }
734
735 }
736
737
738
739 for (int i = 0; i < scale.length; i++) {
740 for (int j = 0; j < scale.length; j++ ) {
741 stm.setEntry(i, j, stm.getEntry(i, j) * scale[j] / scale[i]);
742 }
743 }
744
745
746 return stm;
747
748 }
749
750
751
752
753
754
755 private RealMatrix getMeasurementMatrix() {
756
757
758 final SpacecraftState[] evaluationStates = predictedMeasurement.getStates();
759 final ObservedMeasurement<?> observedMeasurement = predictedMeasurement.getObservedMeasurement();
760 final double[] sigma = observedMeasurement.getTheoreticalStandardDeviation();
761
762
763
764
765 final RealMatrix measurementMatrix = MatrixUtils.
766 createRealMatrix(observedMeasurement.getDimension(),
767 correctedEstimate.getState().getDimension());
768
769
770 for (int k = 0; k < evaluationStates.length; ++k) {
771 final int p = observedMeasurement.getSatellites().get(k).getPropagatorIndex();
772
773
774 final Orbit predictedOrbit = evaluationStates[k].getOrbit();
775
776
777
778
779
780 final double[][] aCY = new double[6][6];
781 predictedOrbit.getJacobianWrtParameters(builders.get(p).getPositionAngle(), aCY);
782 final RealMatrix dCdY = new Array2DRowRealMatrix(aCY, false);
783
784
785 final RealMatrix dMdC = new Array2DRowRealMatrix(predictedMeasurement.getStateDerivatives(k), false);
786
787
788 final RealMatrix dMdY = dMdC.multiply(dCdY);
789
790
791 for (int i = 0; i < dMdY.getRowDimension(); ++i) {
792 int jOrb = orbitsStartColumns[p];
793 for (int j = 0; j < dMdY.getColumnDimension(); ++j) {
794 final ParameterDriver driver = builders.get(p).getOrbitalParametersDrivers().getDrivers().get(j);
795 if (driver.isSelected()) {
796 measurementMatrix.setEntry(i, jOrb++,
797 dMdY.getEntry(i, j) / sigma[i] * driver.getScale());
798 }
799 }
800 }
801
802
803
804
805
806 final int nbParams = estimatedPropagationParameters[p].getNbParams();
807 if (nbParams > 0) {
808 final RealMatrix dYdPp = harvesters[p].getParametersJacobian(evaluationStates[k]);
809 final RealMatrix dMdPp = dMdY.multiply(dYdPp);
810 for (int i = 0; i < dMdPp.getRowDimension(); ++i) {
811 for (int j = 0; j < nbParams; ++j) {
812 final ParameterDriver delegating = estimatedPropagationParameters[p].getDrivers().get(j);
813 measurementMatrix.setEntry(i, propagationParameterColumns.get(delegating.getName()),
814 dMdPp.getEntry(i, j) / sigma[i] * delegating.getScale());
815 }
816 }
817 }
818
819
820
821
822
823
824 for (final ParameterDriver driver : observedMeasurement.getParametersDrivers()) {
825 if (driver.isSelected()) {
826
827 final double[] aMPm = predictedMeasurement.getParameterDerivatives(driver);
828
829
830 if (measurementParameterColumns.get(driver.getName()) != null) {
831
832 final int driverColumn = measurementParameterColumns.get(driver.getName());
833
834
835 for (int i = 0; i < aMPm.length; ++i) {
836 measurementMatrix.setEntry(i, driverColumn,
837 aMPm[i] / sigma[i] * driver.getScale());
838 }
839 }
840 }
841 }
842 }
843
844
845 return measurementMatrix;
846
847 }
848
849
850
851
852
853
854
855
856 private RealMatrix normalizeCovarianceMatrix(final RealMatrix physicalCovarianceMatrix) {
857
858
859 final int nbParams = physicalCovarianceMatrix.getRowDimension();
860 final RealMatrix normalizedCovarianceMatrix = MatrixUtils.createRealMatrix(nbParams, nbParams);
861
862
863 for (int i = 0; i < nbParams; ++i) {
864 for (int j = 0; j < nbParams; ++j) {
865 normalizedCovarianceMatrix.setEntry(i, j,
866 physicalCovarianceMatrix.getEntry(i, j) /
867 (scale[i] * scale[j]));
868 }
869 }
870 return normalizedCovarianceMatrix;
871 }
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887 private <T extends ObservedMeasurement<T>> void applyDynamicOutlierFilter(final EstimatedMeasurement<T> measurement,
888 final RealMatrix innovationCovarianceMatrix) {
889
890
891 final ObservedMeasurement<T> observedMeasurement = measurement.getObservedMeasurement();
892
893
894
895 for (EstimationModifier<T> modifier : observedMeasurement.getModifiers()) {
896 if (modifier instanceof DynamicOutlierFilter<?>) {
897 final DynamicOutlierFilter<T> dynamicOutlierFilter = (DynamicOutlierFilter<T>) modifier;
898
899
900 final double[] sigmaDynamic = new double[innovationCovarianceMatrix.getColumnDimension()];
901 final double[] sigmaMeasurement = observedMeasurement.getTheoreticalStandardDeviation();
902
903
904
905
906
907
908
909
910
911
912
913
914 for (int i = 0; i < sigmaDynamic.length; i++) {
915 sigmaDynamic[i] = FastMath.sqrt(innovationCovarianceMatrix.getEntry(i, i)) * sigmaMeasurement[i];
916 }
917 dynamicOutlierFilter.setSigma(sigmaDynamic);
918
919
920 modifier.modify(measurement);
921
922
923 dynamicOutlierFilter.setSigma(null);
924 }
925 }
926 }
927
928
929 @Override
930 public NonLinearEvolution getEvolution(final double previousTime, final RealVector previousState,
931 final MeasurementDecorator measurement) {
932
933
934 final ObservedMeasurement<?> observedMeasurement = measurement.getObservedMeasurement();
935 for (final ParameterDriver driver : observedMeasurement.getParametersDrivers()) {
936 if (driver.getReferenceDate() == null) {
937 driver.setReferenceDate(builders.get(0).getInitialOrbitDate());
938 }
939 }
940
941 ++currentMeasurementNumber;
942 currentDate = measurement.getObservedMeasurement().getDate();
943
944
945
946
947
948
949
950
951
952
953 final RealVector predictedState = predictState(observedMeasurement.getDate());
954
955
956 final RealMatrix stateTransitionMatrix = getErrorStateTransitionMatrix();
957
958
959
960
961
962
963
964
965
966
967 predictedMeasurement = observedMeasurement.estimate(currentMeasurementNumber,
968 currentMeasurementNumber,
969 filterRelevant(observedMeasurement, predictedSpacecraftStates));
970
971
972 final RealMatrix measurementMatrix = getMeasurementMatrix();
973
974
975 final RealMatrix physicalProcessNoise = MatrixUtils.createRealMatrix(previousState.getDimension(),
976 previousState.getDimension());
977 for (int k = 0; k < covarianceMatricesProviders.size(); ++k) {
978
979
980 final int nbMeas = estimatedMeasurementsParameters.getNbParams();
981
982
983 final int nbDyn = orbitsEndColumns[k] - orbitsStartColumns[k] +
984 estimatedPropagationParameters[k].getNbParams();
985
986
987 final RealMatrix noiseK = MatrixUtils.createRealMatrix(nbDyn + nbMeas, nbDyn + nbMeas);
988 final RealMatrix noiseP = covarianceMatricesProviders.get(k).
989 getProcessNoiseMatrix(correctedSpacecraftStates[k],
990 predictedSpacecraftStates[k]);
991 noiseK.setSubMatrix(noiseP.getData(), 0, 0);
992 if (measurementProcessNoiseMatrix != null) {
993 final RealMatrix noiseM = measurementProcessNoiseMatrix.
994 getProcessNoiseMatrix(correctedSpacecraftStates[k],
995 predictedSpacecraftStates[k]);
996 noiseK.setSubMatrix(noiseM.getData(), nbDyn, nbDyn);
997 }
998
999 checkDimension(noiseK.getRowDimension(),
1000 builders.get(k).getOrbitalParametersDrivers(),
1001 builders.get(k).getPropagationParametersDrivers(),
1002 estimatedMeasurementsParameters);
1003
1004 final int[] indK = covarianceIndirection[k];
1005 for (int i = 0; i < indK.length; ++i) {
1006 if (indK[i] >= 0) {
1007 for (int j = 0; j < indK.length; ++j) {
1008 if (indK[j] >= 0) {
1009 physicalProcessNoise.setEntry(indK[i], indK[j], noiseK.getEntry(i, j));
1010 }
1011 }
1012 }
1013 }
1014
1015 }
1016 final RealMatrix normalizedProcessNoise = normalizeCovarianceMatrix(physicalProcessNoise);
1017
1018 return new NonLinearEvolution(measurement.getTime(), predictedState,
1019 stateTransitionMatrix, normalizedProcessNoise, measurementMatrix);
1020
1021 }
1022
1023
1024
1025 @Override
1026 public RealVector getInnovation(final MeasurementDecorator measurement, final NonLinearEvolution evolution,
1027 final RealMatrix innovationCovarianceMatrix) {
1028
1029
1030 applyDynamicOutlierFilter(predictedMeasurement, innovationCovarianceMatrix);
1031 if (predictedMeasurement.getStatus() == EstimatedMeasurement.Status.REJECTED) {
1032
1033 return null;
1034 } else {
1035
1036 final double[] observed = predictedMeasurement.getObservedMeasurement().getObservedValue();
1037 final double[] estimated = predictedMeasurement.getEstimatedValue();
1038 final double[] sigma = predictedMeasurement.getObservedMeasurement().getTheoreticalStandardDeviation();
1039 final double[] residuals = new double[observed.length];
1040
1041 for (int i = 0; i < observed.length; i++) {
1042 residuals[i] = (observed[i] - estimated[i]) / sigma[i];
1043 }
1044 return MatrixUtils.createRealVector(residuals);
1045 }
1046 }
1047
1048
1049
1050
1051
1052 public void finalizeEstimation(final ObservedMeasurement<?> observedMeasurement,
1053 final ProcessEstimate estimate) {
1054
1055
1056 correctedEstimate = estimate;
1057 updateParameters();
1058
1059
1060
1061 final Propagator[] estimatedPropagators = getEstimatedPropagators();
1062 for (int k = 0; k < estimatedPropagators.length; ++k) {
1063 correctedSpacecraftStates[k] = estimatedPropagators[k].getInitialState();
1064 }
1065
1066
1067 correctedMeasurement = observedMeasurement.estimate(currentMeasurementNumber,
1068 currentMeasurementNumber,
1069 filterRelevant(observedMeasurement, correctedSpacecraftStates));
1070
1071
1072 updateReferenceTrajectories(estimatedPropagators, propagationType, stateType);
1073
1074 }
1075
1076
1077
1078
1079
1080
1081
1082 private SpacecraftState[] filterRelevant(final ObservedMeasurement<?> observedMeasurement, final SpacecraftState[] allStates) {
1083 final List<ObservableSatellite> satellites = observedMeasurement.getSatellites();
1084 final SpacecraftState[] relevantStates = new SpacecraftState[satellites.size()];
1085 for (int i = 0; i < relevantStates.length; ++i) {
1086 relevantStates[i] = allStates[satellites.get(i).getPropagatorIndex()];
1087 }
1088 return relevantStates;
1089 }
1090
1091
1092
1093
1094
1095
1096 private RealVector predictState(final AbsoluteDate date) {
1097
1098
1099 final RealVector predictedState = correctedEstimate.getState().copy();
1100
1101
1102 int jOrb = 0;
1103
1104 for (int k = 0; k < predictedSpacecraftStates.length; ++k) {
1105
1106
1107 predictedSpacecraftStates[k] = referenceTrajectories[k].propagate(date);
1108
1109
1110
1111 builders.get(k).resetOrbit(predictedSpacecraftStates[k].getOrbit());
1112
1113
1114
1115
1116
1117 for (DelegatingDriver orbitalDriver : builders.get(k).getOrbitalParametersDrivers().getDrivers()) {
1118 if (orbitalDriver.isSelected()) {
1119 predictedState.setEntry(jOrb++, orbitalDriver.getNormalizedValue());
1120 }
1121 }
1122
1123 }
1124
1125 return predictedState;
1126
1127 }
1128
1129
1130
1131
1132 private void updateParameters() {
1133 final RealVector correctedState = correctedEstimate.getState();
1134 int i = 0;
1135 for (final DelegatingDriver driver : getEstimatedOrbitalParameters().getDrivers()) {
1136
1137 driver.setNormalizedValue(correctedState.getEntry(i));
1138 correctedState.setEntry(i++, driver.getNormalizedValue());
1139 }
1140 for (final DelegatingDriver driver : getEstimatedPropagationParameters().getDrivers()) {
1141
1142 driver.setNormalizedValue(correctedState.getEntry(i));
1143 correctedState.setEntry(i++, driver.getNormalizedValue());
1144 }
1145 for (final DelegatingDriver driver : getEstimatedMeasurementsParameters().getDrivers()) {
1146
1147 driver.setNormalizedValue(correctedState.getEntry(i));
1148 correctedState.setEntry(i++, driver.getNormalizedValue());
1149 }
1150 }
1151
1152
1153
1154
1155 public List<OrbitDeterminationPropagatorBuilder> getBuilders() {
1156 return builders;
1157 }
1158
1159
1160
1161
1162 public Propagator[] getReferenceTrajectories() {
1163 return referenceTrajectories.clone();
1164 }
1165
1166
1167
1168
1169 public void setReferenceTrajectories(final Propagator[] referenceTrajectories) {
1170 this.referenceTrajectories = referenceTrajectories.clone();
1171 }
1172
1173
1174
1175
1176
1177 @Deprecated
1178 public AbstractJacobiansMapper[] getMappers() {
1179 return null;
1180 }
1181
1182
1183
1184
1185
1186 @Deprecated
1187 public void setMappers(final AbstractJacobiansMapper[] mappers) {
1188 setHarvesters(mappers);
1189 }
1190
1191
1192
1193
1194
1195 public void setHarvesters(final MatricesHarvester[] harvesters) {
1196 this.harvesters = harvesters.clone();
1197 }
1198
1199
1200
1201
1202 public Propagator[] getEstimatedPropagators() {
1203
1204 final Propagator[] propagators = new Propagator[getBuilders().size()];
1205 for (int k = 0; k < getBuilders().size(); ++k) {
1206 propagators[k] = getBuilders().get(k).buildPropagator(getBuilders().get(k).getSelectedNormalizedParameters());
1207 }
1208 return propagators;
1209 }
1210
1211 }