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