1   /* Copyright 2002-2022 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.List;
20  
21  import org.hipparchus.exception.MathRuntimeException;
22  import org.hipparchus.filtering.kalman.extended.ExtendedKalmanFilter;
23  import org.hipparchus.linear.MatrixDecomposer;
24  import org.hipparchus.linear.RealMatrix;
25  import org.hipparchus.linear.RealVector;
26  import org.orekit.errors.OrekitException;
27  import org.orekit.estimation.measurements.ObservedMeasurement;
28  import org.orekit.propagation.conversion.DSSTPropagatorBuilder;
29  import org.orekit.propagation.semianalytical.dsst.DSSTPropagator;
30  import org.orekit.time.AbsoluteDate;
31  import org.orekit.utils.ParameterDriver;
32  import org.orekit.utils.ParameterDriversList;
33  import org.orekit.utils.ParameterDriversList.DelegatingDriver;
34  
35  /**
36   * Implementation of an Extended Semi-analytical Kalman Filter (ESKF) to perform orbit determination.
37   * <p>
38   * The filter uses a {@link DSSTPropagatorBuilder}.
39   * </p>
40   * <p>
41   * The estimated parameters are driven by {@link ParameterDriver} objects. They are of 3 different types:<ol>
42   *   <li><b>Orbital parameters</b>:The position and velocity of the spacecraft, or, more generally, its orbit.<br>
43   *       These parameters are retrieved from the reference trajectory propagator builder when the filter is initialized.</li>
44   *   <li><b>Propagation parameters</b>: Some parameters modelling physical processes (SRP or drag coefficients).<br>
45   *       They are also retrieved from the propagator builder during the initialization phase.</li>
46   *   <li><b>Measurements parameters</b>: Parameters related to measurements (station biases, positions etc...).<br>
47   *       They are passed down to the filter in its constructor.</li>
48   * </ol>
49   * <p>
50   * The Kalman filter implementation used is provided by the underlying mathematical library Hipparchus.
51   * All the variables seen by Hipparchus (states, covariances, measurement matrices...) are normalized
52   * using a specific scale for each estimated parameters or standard deviation noise for each measurement components.
53   * </p>
54   *
55   * @see "Folcik Z., Orbit Determination Using Modern Filters/Smoothers and Continuous Thrust Modeling,
56   *       Master of Science Thesis, Department of Aeronautics and Astronautics, MIT, June, 2008."
57   *
58   * @see "Cazabonne B., Bayard J., Journot M., and Cefola P. J., A Semi-analytical Approach for Orbit
59   *       Determination based on Extended Kalman Filter, AAS Paper 21-614, AAS/AIAA Astrodynamics
60   *       Specialist Conference, Big Sky, August 2021."
61   *
62   * @author Julie Bayard
63   * @author Bryan Cazabonne
64   * @author Maxime Journot
65   * @since 11.1
66   */
67  public class SemiAnalyticalKalmanEstimator {
68  
69      /** Builders for orbit propagators. */
70      private DSSTPropagatorBuilder propagatorBuilder;
71  
72      /** Kalman filter process model. */
73      private final SemiAnalyticalKalmanModel processModel;
74  
75      /** Filter. */
76      private final ExtendedKalmanFilter<MeasurementDecorator> filter;
77  
78      /** Observer to retrieve current estimation info. */
79      private KalmanObserver observer;
80  
81      /** Kalman filter estimator constructor (package private).
82       * @param decomposer decomposer to use for the correction phase
83       * @param propagatorBuilder propagator builder used to evaluate the orbit.
84       * @param covarianceMatrixProvider provider for process noise matrix
85       * @param estimatedMeasurementParameters measurement parameters to estimate
86       * @param measurementProcessNoiseMatrix provider for measurement process noise matrix
87       */
88      public SemiAnalyticalKalmanEstimator(final MatrixDecomposer decomposer,
89                                           final DSSTPropagatorBuilder propagatorBuilder,
90                                           final CovarianceMatrixProvider covarianceMatrixProvider,
91                                           final ParameterDriversList estimatedMeasurementParameters,
92                                           final CovarianceMatrixProvider measurementProcessNoiseMatrix) {
93  
94          this.propagatorBuilder = propagatorBuilder;
95          this.observer          = null;
96  
97          // Build the process model and measurement model
98          this.processModel = new SemiAnalyticalKalmanModel(propagatorBuilder, covarianceMatrixProvider,
99                                                            estimatedMeasurementParameters,  measurementProcessNoiseMatrix);
100 
101         // Extended Kalman Filter of Hipparchus
102         this.filter = new ExtendedKalmanFilter<>(decomposer, processModel, processModel.getEstimate());
103 
104     }
105 
106     /** Set the observer.
107      * @param observer the observer
108      */
109     public void setObserver(final KalmanObserver observer) {
110         this.observer = observer;
111     }
112 
113     /** Get the current measurement number.
114      * @return current measurement number
115      */
116     public int getCurrentMeasurementNumber() {
117         return processModel.getCurrentMeasurementNumber();
118     }
119 
120     /** Get the current date.
121      * @return current date
122      */
123     public AbsoluteDate getCurrentDate() {
124         return processModel.getCurrentDate();
125     }
126 
127     /** Get the "physical" estimated state (i.e. not normalized)
128      * @return the "physical" estimated state
129      */
130     public RealVector getPhysicalEstimatedState() {
131         return processModel.getPhysicalEstimatedState();
132     }
133 
134     /** Get the "physical" estimated covariance matrix (i.e. not normalized)
135      * @return the "physical" estimated covariance matrix
136      */
137     public RealMatrix getPhysicalEstimatedCovarianceMatrix() {
138         return processModel.getPhysicalEstimatedCovarianceMatrix();
139     }
140 
141     /** Get the orbital parameters supported by this estimator.
142      * <p>
143      * If there are more than one propagator builder, then the names
144      * of the drivers have an index marker in square brackets appended
145      * to them in order to distinguish the various orbits. So for example
146      * with one builder generating Keplerian orbits the names would be
147      * simply "a", "e", "i"... but if there are several builders the
148      * names would be "a[0]", "e[0]", "i[0]"..."a[1]", "e[1]", "i[1]"...
149      * </p>
150      * @param estimatedOnly if true, only estimated parameters are returned
151      * @return orbital parameters supported by this estimator
152      */
153     public ParameterDriversList getOrbitalParametersDrivers(final boolean estimatedOnly) {
154 
155         final ParameterDriversList estimated = new ParameterDriversList();
156         for (final ParameterDriver driver : propagatorBuilder.getOrbitalParametersDrivers().getDrivers()) {
157             if (driver.isSelected() || !estimatedOnly) {
158                 driver.setName(driver.getName());
159                 estimated.add(driver);
160             }
161         }
162         return estimated;
163     }
164 
165     /** Get the propagator parameters supported by this estimator.
166      * @param estimatedOnly if true, only estimated parameters are returned
167      * @return propagator parameters supported by this estimator
168      */
169     public ParameterDriversList getPropagationParametersDrivers(final boolean estimatedOnly) {
170 
171         final ParameterDriversList estimated = new ParameterDriversList();
172         for (final DelegatingDriver delegating : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
173             if (delegating.isSelected() || !estimatedOnly) {
174                 for (final ParameterDriver driver : delegating.getRawDrivers()) {
175                     estimated.add(driver);
176                 }
177             }
178         }
179         return estimated;
180     }
181 
182     /** Get the list of estimated measurements parameters.
183      * @return the list of estimated measurements parameters
184      */
185     public ParameterDriversList getEstimatedMeasurementsParameters() {
186         return processModel.getEstimatedMeasurementsParameters();
187     }
188 
189     /** Process a single measurement.
190      * <p>
191      * Update the filter with the new measurement by calling the estimate method.
192      * </p>
193      * @param observedMeasurements the list of measurements to process
194      * @return estimated propagators
195      */
196     public DSSTPropagator processMeasurements(final List<ObservedMeasurement<?>> observedMeasurements) {
197         try {
198             processModel.setObserver(observer);
199             return processModel.processMeasurements(observedMeasurements, filter);
200         } catch (MathRuntimeException mrte) {
201             throw new OrekitException(mrte);
202         }
203     }
204 
205 }
206