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.List;
21
22 import org.hipparchus.linear.MatrixUtils;
23 import org.hipparchus.linear.RealMatrix;
24 import org.hipparchus.linear.RealVector;
25 import org.hipparchus.util.MerweUnscentedTransform;
26 import org.junit.jupiter.api.Assertions;
27 import org.junit.jupiter.api.BeforeEach;
28 import org.junit.jupiter.api.Test;
29 import org.orekit.estimation.DSSTContext;
30 import org.orekit.estimation.DSSTEstimationTestUtils;
31 import org.orekit.estimation.DSSTForce;
32 import org.orekit.estimation.measurements.GroundStation;
33 import org.orekit.estimation.measurements.ObservableSatellite;
34 import org.orekit.estimation.measurements.Range;
35 import org.orekit.estimation.measurements.modifiers.Bias;
36 import org.orekit.forces.radiation.RadiationSensitive;
37 import org.orekit.orbits.Orbit;
38 import org.orekit.orbits.OrbitType;
39 import org.orekit.orbits.PositionAngleType;
40 import org.orekit.propagation.PropagationType;
41 import org.orekit.propagation.SpacecraftState;
42 import org.orekit.propagation.conversion.DSSTPropagatorBuilder;
43 import org.orekit.time.AbsoluteDate;
44 import org.orekit.utils.ParameterDriver;
45 import org.orekit.utils.ParameterDriversList;
46
47 public class SemiAnalyticalUnscentedKalmanModelTest {
48
49
50
51 private final OrbitType orbitType = OrbitType.EQUINOCTIAL;
52
53
54 private final PositionAngleType positionAngleType = PositionAngleType.MEAN;
55
56
57 private Orbit orbit0;
58
59
60 private DSSTPropagatorBuilder propagatorBuilder;
61
62
63 private CovarianceMatrixProvider covMatrixProvider;
64
65
66 private ParameterDriversList estimatedMeasurementsParameters;
67
68
69 private SemiAnalyticalUnscentedKalmanEstimator kalman;
70
71
72 private ModelLogger modelLogger;
73
74
75 private int M;
76
77
78 private Range range;
79
80
81 private ParameterDriver satRangeBiasDriver;
82
83
84 private ParameterDriver srpCoefDriver;
85
86
87 private final double tol = 5.0e-9;
88
89 @BeforeEach
90 public void setup() {
91
92 final DSSTContext context = DSSTEstimationTestUtils.eccentricContext("regular-data:potential:tides");
93
94
95 this.orbit0 = context.initialOrbit;
96 ObservableSatellite sat = new ObservableSatellite(0);
97
98
99 this.propagatorBuilder = context.createBuilder(PropagationType.MEAN, PropagationType.OSCULATING, true,
100 1.0e-6, 60.0, 10., DSSTForce.SOLAR_RADIATION_PRESSURE);
101
102
103 final AbsoluteDate date0 = context.initialOrbit.getDate();
104
105
106 final AbsoluteDate date = date0.shiftedBy(10.);
107 final GroundStation station = context.stations.get(0);
108 this.range = new Range(station, true, date, 18616150., 10., 1., sat);
109
110
111
112 final Bias<Range> satRangeBias = new Bias<Range>(new String[] {"sat range bias"},
113 new double[] {100.},
114 new double[] {10.},
115 new double[] {0.},
116 new double[] {100.});
117 this.satRangeBiasDriver = satRangeBias.getParametersDrivers().get(0);
118 satRangeBiasDriver.setSelected(true);
119 satRangeBiasDriver.setReferenceDate(date);
120 range.addModifier(satRangeBias);
121 for (ParameterDriver driver : range.getParametersDrivers()) {
122 driver.setReferenceDate(date);
123 }
124
125
126 this.estimatedMeasurementsParameters = new ParameterDriversList();
127 for (final ParameterDriver driver : range.getParametersDrivers()) {
128 if (driver.isSelected()) {
129 estimatedMeasurementsParameters.add(driver);
130 }
131 }
132
133 this.srpCoefDriver = propagatorBuilder.getPropagationParametersDrivers().
134 findByName(RadiationSensitive.REFLECTION_COEFFICIENT);
135 srpCoefDriver.setReferenceDate(date);
136 srpCoefDriver.setSelected(true);
137
138
139 final double[] scales = getParametersScale(propagatorBuilder, estimatedMeasurementsParameters);
140 this.M = scales.length;
141 this.covMatrixProvider = setInitialCovarianceMatrix(scales);
142
143
144 final SemiAnalyticalUnscentedKalmanEstimatorBuilder kalmanBuilder = new SemiAnalyticalUnscentedKalmanEstimatorBuilder();
145 kalmanBuilder.unscentedTransformProvider(new MerweUnscentedTransform(8));
146 kalmanBuilder.addPropagationConfiguration(propagatorBuilder, covMatrixProvider);
147 kalmanBuilder.estimatedMeasurementsParameters(estimatedMeasurementsParameters, null);
148 this.kalman = kalmanBuilder.build();
149 this.modelLogger = new ModelLogger();
150 kalman.setObserver(modelLogger);
151 }
152
153 @Test
154 public void ModelPhysicalOutputsTest() {
155
156
157
158 checkModelAtT0();
159
160 }
161
162
163 private void checkModelAtT0() {
164
165
166 final SemiAnalyticalUnscentedKalmanModel model = new SemiAnalyticalUnscentedKalmanModel(propagatorBuilder,
167 covMatrixProvider,
168 estimatedMeasurementsParameters,
169 null);
170 model.setObserver(modelLogger);
171
172
173
174
175
176 Assertions.assertNotNull(model.getObserver());
177
178
179 Assertions.assertEquals(0., model.getEstimate().getTime(), 0.);
180 Assertions.assertEquals(0., model.getCurrentDate().durationFrom(orbit0.getDate()), 0.);
181
182
183 Assertions.assertEquals(0, model.getCurrentMeasurementNumber());
184
185
186 final RealVector expX = MatrixUtils.createRealVector(M);
187 final double[] orbitState0 = new double[6];
188 orbitType.mapOrbitToArray(orbit0, positionAngleType, orbitState0, null);
189 expX.setSubVector(0, MatrixUtils.createRealVector(orbitState0));
190 expX.setEntry(6, srpCoefDriver.getReferenceValue());
191 expX.setEntry(7, satRangeBiasDriver.getReferenceValue());
192 Assertions.assertArrayEquals(model.getPhysicalEstimatedState().toArray(), expX.toArray(), tol);
193 Assertions.assertArrayEquals(model.getEstimate().getState().toArray(), new double[8], tol);
194
195
196 final double[][] Pn = model.getEstimate().getCovariance().getData();
197 final double[][] expPn = covMatrixProvider.getInitialCovarianceMatrix(null).getData();
198 for (int i = 0; i < M; i++) {
199 Assertions.assertArrayEquals(expPn[i], Pn[i], tol, "Failed on line " + i);
200 }
201
202
203 final RealMatrix P = model.getPhysicalEstimatedCovarianceMatrix();
204 final RealMatrix expP = covMatrixProvider.getInitialCovarianceMatrix(new SpacecraftState(orbit0));
205 final double[][] dP = P.subtract(expP).getData();
206 for (int i = 0; i < M; i++) {
207 Assertions.assertArrayEquals(new double[M], dP[i], tol, "Failed on line " + i);
208 }
209
210
211 Assertions.assertNull(model.getPhysicalStateTransitionMatrix());
212
213
214 Assertions.assertNull(model.getPhysicalMeasurementJacobian());
215
216
217 Assertions.assertNull(model.getEstimate().getInnovationCovariance());
218 Assertions.assertNull(model.getPhysicalInnovationCovarianceMatrix());
219 Assertions.assertNull(model.getEstimate().getKalmanGain());
220 Assertions.assertNull(model.getPhysicalKalmanGain());
221 Assertions.assertNull(model.getEstimate().getMeasurementJacobian());
222 Assertions.assertNull(model.getPhysicalMeasurementJacobian());
223 Assertions.assertNull(model.getEstimate().getStateTransitionMatrix());
224 Assertions.assertNull(model.getPhysicalStateTransitionMatrix());
225 }
226
227
228
229
230
231
232 private double[] getParametersScale(final DSSTPropagatorBuilder builder,
233 ParameterDriversList estimatedMeasurementsParameters) {
234 final List<Double> scaleList = new ArrayList<>();
235
236
237 for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
238 if (driver.isSelected()) {
239 scaleList.add(driver.getScale());
240 }
241 }
242
243
244 for (ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) {
245 if (driver.isSelected()) {
246 scaleList.add(driver.getScale());
247 }
248 }
249
250
251 for (ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
252 if (driver.isSelected()) {
253 scaleList.add(driver.getScale());
254 }
255 }
256
257 final double[] scales = new double[scaleList.size()];
258 for (int i = 0; i < scaleList.size(); i++) {
259 scales[i] = scaleList.get(i);
260 }
261 return scales;
262 }
263
264
265
266
267
268
269
270 private CovarianceMatrixProvider setInitialCovarianceMatrix(final double[] scales) {
271
272 final int n = scales.length;
273 final RealMatrix cov = MatrixUtils.createRealMatrix(n, n);
274 for (int i = 0; i < n; i++) {
275 for (int j = 0; j < n; j++) {
276 cov.setEntry(i, j, scales[i] * scales[j]);
277 }
278 }
279 return new ConstantProcessNoise(cov);
280 }
281
282
283
284 public class ModelLogger implements KalmanObserver {
285 KalmanEstimation estimation;
286
287 @Override
288 public void evaluationPerformed(KalmanEstimation estimation) {
289 this.estimation = estimation;
290 }
291 }
292
293 }