1   /* Copyright 2002-2025 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
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      /** Orbit type for propagation. */
51      private final OrbitType orbitType = OrbitType.EQUINOCTIAL;
52  
53      /** Position angle for propagation. */
54      private final PositionAngleType positionAngleType = PositionAngleType.MEAN;
55  
56      /** Initial orbit. */
57      private Orbit orbit0;
58  
59      /** Propagator builder. */
60      private DSSTPropagatorBuilder propagatorBuilder;
61  
62      /** Covariance matrix provider. */
63      private CovarianceMatrixProvider covMatrixProvider;
64  
65      /** Estimated measurement parameters list. */
66      private ParameterDriversList estimatedMeasurementsParameters;
67  
68      /** Kalman estimator containing models. */
69      private SemiAnalyticalUnscentedKalmanEstimator kalman;
70  
71      /** Kalman observer. */
72      private ModelLogger modelLogger;
73  
74      /** State size. */
75      private int M;
76  
77      /** Range after t0. */
78      private Range range;
79  
80      /** Driver for satellite range bias. */
81      private ParameterDriver satRangeBiasDriver;
82  
83      /** Driver for SRP coefficient. */
84      private ParameterDriver srpCoefDriver;
85  
86      /** Tolerance for the test. */
87      private final double tol = 5.0e-9;
88  
89      @BeforeEach
90      public void setup() {
91          // Create context
92          final DSSTContext context = DSSTEstimationTestUtils.eccentricContext("regular-data:potential:tides");
93  
94          // Initial orbit and date
95          this.orbit0 = context.initialOrbit;
96          ObservableSatellite sat = new ObservableSatellite(0);
97  
98          // Create propagator builder
99          this.propagatorBuilder = context.createBuilder(PropagationType.MEAN, PropagationType.OSCULATING, true,
100                                                        1.0e-6, 60.0, 10., DSSTForce.SOLAR_RADIATION_PRESSURE);
101 
102         //  t0
103         final AbsoluteDate date0 = context.initialOrbit.getDate();
104 
105         // Create one 0m range measurement at t0 + 10s
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         // Exact range value is 1.8616150246470984E7 m
110 
111         // Add sat range bias to PV and select it
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         // Gather list of meas parameters (only sat range bias here)
126         this.estimatedMeasurementsParameters = new ParameterDriversList();
127         for (final ParameterDriver driver : range.getParametersDrivers()) {
128             if (driver.isSelected()) {
129                 estimatedMeasurementsParameters.add(driver);
130             }
131         }
132         // Select SRP coefficient
133         this.srpCoefDriver = propagatorBuilder.getPropagationParametersDrivers().
134                         findByName(RadiationSensitive.REFLECTION_COEFFICIENT);
135         srpCoefDriver.setReferenceDate(date);
136         srpCoefDriver.setSelected(true);
137 
138         // Create a covariance matrix using the scales of the estimated parameters
139         final double[] scales = getParametersScale(propagatorBuilder, estimatedMeasurementsParameters);
140         this.M = scales.length;
141         this.covMatrixProvider = setInitialCovarianceMatrix(scales);
142 
143         // Initialize Kalman
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         // Check model at t0 before any measurement is added
157         // -------------------------------------------------
158         checkModelAtT0();
159 
160     }
161 
162     /** Check the model physical outputs at t0 before any measurement is added. */
163     private void checkModelAtT0() {
164 
165         // Instantiate a Model from attributes
166         final SemiAnalyticalUnscentedKalmanModel model = new SemiAnalyticalUnscentedKalmanModel(propagatorBuilder,
167                                                                                                 covMatrixProvider,
168                                                                                                 estimatedMeasurementsParameters,
169                                                                                                 null);
170         model.setObserver(modelLogger);
171 
172         // Evaluate at t0
173         // --------------
174 
175         // Observer
176         Assertions.assertNotNull(model.getObserver());
177 
178         // Time
179         Assertions.assertEquals(0., model.getEstimate().getTime(), 0.);
180         Assertions.assertEquals(0., model.getCurrentDate().durationFrom(orbit0.getDate()), 0.);
181 
182         // Measurement number
183         Assertions.assertEquals(0, model.getCurrentMeasurementNumber());
184 
185         // Physical state and predicted filter correction
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         // Normalized covariance - filled with 1
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         // Physical covariance = initialized
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         // State transition matrix (null for unscented kalman filter)
211         Assertions.assertNull(model.getPhysicalStateTransitionMatrix());
212 
213         // Measurement matrix (null for unscented kalman filter)
214         Assertions.assertNull(model.getPhysicalMeasurementJacobian());
215 
216         // Check that other "physical" matrices are null
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     /** Get an array of the scales of the estimated parameters.
228      * @param builder propagator builder
229      * @param estimatedMeasurementsParameters estimated measurements parameters
230      * @return array containing the scales of the estimated parameter
231      */
232     private double[] getParametersScale(final DSSTPropagatorBuilder builder,
233                                        ParameterDriversList estimatedMeasurementsParameters) {
234         final List<Double> scaleList = new ArrayList<>();
235 
236         // Orbital parameters
237         for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
238             if (driver.isSelected()) {
239                 scaleList.add(driver.getScale());
240             }
241         }
242 
243         // Propagation parameters
244         for (ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) {
245             if (driver.isSelected()) {
246                 scaleList.add(driver.getScale());
247             }
248         }
249 
250         // Measurement parameters
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     /** Create a covariance matrix provider with initial and process noise matrix constant and identical.
265      * Each element Pij of the initial covariance matrix equals:
266      * Pij = scales[i]*scales[j]
267      * @param scales scales of the estimated parameters
268      * @return covariance matrix provider
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     /** Observer allowing to get Kalman model after a measurement was processed in the Kalman filter. */
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 }