Commit 43c50468 authored by Luc Maisonobe's avatar Luc Maisonobe

Merge branch 'parameters-drivers' into orbit-determination

parents 142e188f 93ff64c5
......@@ -284,6 +284,13 @@
<version>${orekit.hipparchus.version}</version>
<type>jar</type>
<optional>false</optional>
</dependency>
<dependency>
<groupId>org.hipparchus</groupId>
<artifactId>hipparchus-stat</artifactId>
<version>${orekit.hipparchus.version}</version>
<type>jar</type>
<optional>false</optional>
</dependency>
<dependency>
<groupId>junit</groupId>
......@@ -551,6 +558,35 @@
</configuration>
</plugin>
</plugins>
<pluginManagement>
<plugins>
<!--This plugin's configuration is used to store Eclipse m2e settings only. It has no influence on the Maven build itself.-->
<plugin>
<groupId>org.eclipse.m2e</groupId>
<artifactId>lifecycle-mapping</artifactId>
<version>1.0.0</version>
<configuration>
<lifecycleMappingMetadata>
<pluginExecutions>
<pluginExecution>
<pluginExecutionFilter>
<groupId>org.apache.felix</groupId>
<artifactId>maven-bundle-plugin</artifactId>
<versionRange>[${orekit.maven-bundle-plugin.version},)</versionRange>
<goals>
<goal>manifest</goal>
</goals>
</pluginExecutionFilter>
<action>
<ignore></ignore>
</action>
</pluginExecution>
</pluginExecutions>
</lifecycleMappingMetadata>
</configuration>
</plugin>
</plugins>
</pluginManagement>
</build>
<reporting>
......
......@@ -132,14 +132,12 @@ public enum OrekitMessages implements Localizable {
OUT_OF_RANGE_EPHEMERIDES_DATE("out of range date for ephemerides: {0}, [{1}, {2}]"),
UNEXPECTED_TWO_ELEVATION_VALUES_FOR_ONE_AZIMUTH("unexpected two elevation values: {0} and {1}, for one azimuth: {2}"),
UNSUPPORTED_PARAMETER_NAME("unsupported parameter name {0}, supported names: {1}"),
DUPLICATED_PARAMETER_NAME("duplicated parameter name {0}"),
UNKNOWN_ADDITIONAL_STATE("unknown additional state \"{0}\""),
UNKNOWN_MONTH("unknown month \"{0}\""),
SINGULAR_JACOBIAN_FOR_ORBIT_TYPE("Jacobian matrix for type {0} is singular with current orbit"),
STATE_JACOBIAN_NOT_INITIALIZED("state Jacobian has not been initialized yet"),
STATE_JACOBIAN_NEITHER_6X6_NOR_7X7("state Jacobian is a {0}x{1} matrix, it should be either a 6x6 or a 7x7 matrix"),
STATE_AND_PARAMETERS_JACOBIANS_ROWS_MISMATCH("state Jacobian has {0} rows but parameters Jacobian has {1} rows"),
STEPS_NOT_INITIALIZED_FOR_FINITE_DIFFERENCES("finite differences steps for force models derivatives not initialized"),
INITIAL_MATRIX_AND_PARAMETERS_NUMBER_MISMATCH("initial Jacobian matrix has {0} columns, but {1} parameters have been selected"),
ORBIT_A_E_MISMATCH_WITH_CONIC_TYPE("orbit should be either elliptic with a > 0 and e < 1 or hyperbolic with a < 0 and e > 1, a = {0}, e = {1}"),
ORBIT_ANOMALY_OUT_OF_HYPERBOLIC_RANGE("true anomaly {0} out of hyperbolic range (e = {1}, {2} < v < {3})"),
......
......@@ -18,16 +18,20 @@ package org.orekit.estimation;
import org.hipparchus.analysis.MultivariateMatrixFunction;
import org.hipparchus.analysis.MultivariateVectorFunction;
import org.hipparchus.analysis.UnivariateFunction;
import org.hipparchus.analysis.UnivariateVectorFunction;
import org.hipparchus.analysis.differentiation.DerivativeStructure;
import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
import org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitExceptionWrapper;
import org.orekit.errors.OrekitMessages;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.utils.ParameterDriver;
/** Utility class for orbit determination.
* @author Luc Maisonobe
......@@ -40,6 +44,60 @@ public class EstimationUtils {
private EstimationUtils() {
}
/** Differentiate a scalar function using finite differences.
* @param function function to differentiate
* @param driver driver for the parameter
* @param nbPoints number of points used for finite differences
* @param step step for finite differences
* @return scalar function evaluating to the derivative of the original function
*/
public static ParameterFunction differentiate(final ParameterFunction function,
final ParameterDriver driver,
final int nbPoints, final double step) {
final UnivariateFunction uf = new UnivariateFunction() {
/** {@inheritDoc} */
@Override
public double value(final double normalizedValue)
throws OrekitExceptionWrapper {
try {
final double saved = driver.getNormalizedValue();
driver.setNormalizedValue(normalizedValue);
final double functionValue = function.value(driver);
driver.setNormalizedValue(saved);
return functionValue;
} catch (OrekitException oe) {
throw new OrekitExceptionWrapper(oe);
}
}
};
final FiniteDifferencesDifferentiator differentiator =
new FiniteDifferencesDifferentiator(nbPoints, step);
final UnivariateDifferentiableFunction differentiated =
differentiator.differentiate(uf);
return new ParameterFunction() {
/** {@inheritDoc} */
@Override
public double value(final ParameterDriver parameterDriver)
throws OrekitException {
if (!parameterDriver.getName().equals(driver.getName())) {
throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME,
parameterDriver.getName(), driver.getName());
}
try {
final DerivativeStructure dsParam = new DerivativeStructure(1, 1, 0, parameterDriver.getNormalizedValue());
final DerivativeStructure dsValue = differentiated.value(dsParam);
return dsValue.getPartialDerivative(1);
} catch (OrekitExceptionWrapper oew) {
throw oew.getException();
}
}
};
}
/** Differentiate a vector function using finite differences.
* @param function function to differentiate
* @param dimension dimension of the vector value of the function
......
/* Copyright 2002-2016 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.estimation;
import org.hipparchus.linear.RealVector;
import org.hipparchus.optim.nonlinear.vector.leastsquares.ParameterValidator;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.Precision;
import org.orekit.errors.OrekitInternalError;
import org.orekit.orbits.OrbitType;
/** Enumerate for validating orbital parameters.
* <p>
* This class prevents the orbit determination engine to use
* inconsistent orbital parameters by trimming them back
* to fit domain boundary. As an example, if an optimizer
* attempts to evaluate a point containing an orbit with
* negative inclination, the inclination will be reset to
* 0.0 by the validator before it is really used in the
* space flight dynamics code.
* </p>
* @author Luc Maisonobe
* @since 8.0
*/
public enum OrbitValidator implements ParameterValidator {
/** Validator for {@link OrbitType#CARTESIAN Cartesian} orbits. */
CARTESIAN() {
/** {@inheritDoc} */
@Override
public RealVector validate(final RealVector params) {
// nothing to validate for Cartesian parameters
return params;
}
},
/** Validator for {@link OrbitType#CIRCULAR circular} orbits. */
CIRCULAR() {
/** {@inheritDoc} */
@Override
public RealVector validate(final RealVector params) {
// ensure semi-major axis is positive
if (params.getEntry(0) <= 0) {
params.setEntry(0, Precision.SAFE_MIN);
}
// ensure eccentricity is less than 1
final double e = FastMath.hypot(params.getEntry(1), params.getEntry(2));
if (e >= 1.0) {
params.setEntry(1, params.getEntry(1) / FastMath.nextUp(e));
params.setEntry(2, params.getEntry(2) / FastMath.nextUp(e));
}
// ensure inclination is non-negative
if (params.getEntry(3) < 0.0) {
params.setEntry(3, 0.0);
}
return params;
}
},
/** Validator for {@link OrbitType#KEPLERIAN Keplerian} orbits. */
KEPLERIAN() {
/** {@inheritDoc} */
@Override
public RealVector validate(final RealVector params) {
// ensure inclination is non-negative
if (params.getEntry(2) < 0.0) {
params.setEntry(2, 0.0);
}
final double a = params.getEntry(0);
final double e = params.getEntry(1);
if (a > 0 && e > 0 && e < 1) {
// regular elliptic orbit, nothing to do
return params;
} else if (a < 0 && e > 1) {
// regular hyperbolic orbit, nothing to do
return params;
} else {
// inconsistent parameters, force orbit to elliptic
params.setEntry(0, FastMath.max(a, Precision.SAFE_MIN));
params.setEntry(1, FastMath.max(0.0, FastMath.min(e, FastMath.nextDown(1.0))));
return params;
}
}
},
/** Validator for {@link OrbitType#EQUINOCTIAL equinoctial} orbits. */
EQUINOCTIAL() {
/** {@inheritDoc} */
@Override
public RealVector validate(final RealVector params) {
// ensure semi-major axis is positive
if (params.getEntry(0) <= 0) {
params.setEntry(0, Precision.SAFE_MIN);
}
// ensure eccentricity is less than 1
final double e = FastMath.hypot(params.getEntry(1), params.getEntry(2));
if (e >= 1.0) {
params.setEntry(1, params.getEntry(1) / FastMath.nextUp(e));
params.setEntry(2, params.getEntry(2) / FastMath.nextUp(e));
}
return params;
}
};
/** Get the validator corresponding to an orbit type.
* @param type orbite type
* @return validator for the orbit type
*/
public static OrbitValidator getValidator(final OrbitType type) {
switch (type) {
case CARTESIAN:
return CARTESIAN;
case CIRCULAR:
return CIRCULAR;
case KEPLERIAN:
return KEPLERIAN;
case EQUINOCTIAL:
return EQUINOCTIAL;
default :
// this should never happen
throw new OrekitInternalError(null);
}
}
}
/* Copyright 2010-2011 Centre National d'Études Spatiales
/* Copyright 2002-2016 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
......@@ -14,61 +14,23 @@
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.propagation.numerical;
package org.orekit.estimation;
import org.orekit.forces.ForceModel;
import org.orekit.errors.OrekitException;
import org.orekit.utils.ParameterDriver;
/** Simple container associating a parameter name with a step to compute its jacobian
* and the provider that manages it.
* @author V&eacute;ronique Pommier-Maurussane
/** Interface representing a scalar function depending on a {@link ParameterDriver}.
* @see EstimationUtils#differentiate(ParameterFunction ParameterDriver, int, double)
* @author Luc Maisonobe
* @since 8.0
*/
public class ParameterConfiguration {
/** Parameter name. */
private String parameterName;
/** Parameter step for finite difference computation of partial derivative with respect to that parameter. */
private double hP;
/** Provider handling this parameter. */
private ForceModel provider;
/** Parameter name and step pair constructor.
* @param parameterName parameter name
* @param hP parameter step */
public ParameterConfiguration(final String parameterName, final double hP) {
this.parameterName = parameterName;
this.hP = hP;
this.provider = null;
}
/** Get parameter name.
* @return parameterName parameter name
*/
public String getParameterName() {
return parameterName;
}
/** Get parameter step.
* @return hP parameter step
*/
public double getHP() {
return hP;
}
/** Set the povider handling this parameter.
* @param provider provider handling this parameter
*/
public void setProvider(final ForceModel provider) {
this.provider = provider;
}
public interface ParameterFunction {
/** Get the povider handling this parameter.
* @return provider handling this parameter
/** Evaluate the function.
* @param parameterDriver friver for the parameter.
* @return scalar value of the function
* @throws OrekitException if evaluation cannot be performed
*/
public ForceModel getProvider() {
return provider;
}
double value(ParameterDriver parameterDriver) throws OrekitException;
}
......@@ -16,13 +16,11 @@
*/
package org.orekit.estimation.leastsquares;
import java.util.List;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem;
import org.orekit.errors.OrekitException;
import org.orekit.estimation.measurements.EvaluationsProvider;
import org.orekit.orbits.Orbit;
import org.orekit.utils.ParameterDriver;
import org.orekit.utils.ParameterDriversList;
/** Observer for {@link BatchLSEstimator batch least squares estimator} iterations.
* <p>
......@@ -49,8 +47,8 @@ public interface BatchLSObserver {
* being able to porvide an evaluation)
*/
void iterationPerformed(int iterationsCount, int evaluationscount, Orbit orbit,
List<ParameterDriver> estimatedPropagatorParameters,
List<ParameterDriver> estimatedMeasurementsParameters,
ParameterDriversList estimatedPropagatorParameters,
ParameterDriversList estimatedMeasurementsParameters,
EvaluationsProvider evaluationsProvider,
LeastSquaresProblem.Evaluation lspEvaluation)
throws OrekitException;
......
......@@ -203,9 +203,7 @@ public abstract class AbstractMeasurement<T extends Measurement<T>> implements M
throws OrekitException {
// combine the measurement parameters and the modifier parameters
for (final ParameterDriver parameterDriver : modifier.getParametersDrivers()) {
parameterDriver.checkAndAddSelf(supportedParameters);
}
supportedParameters.addAll(modifier.getParametersDrivers());
modifiers.add(modifier);
......
......@@ -26,7 +26,6 @@ import org.orekit.frames.Transform;
import org.orekit.propagation.SpacecraftState;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.AngularCoordinates;
import org.orekit.utils.ParameterDriver;
/** Class modeling an Azimuth-Elevation measurement from a ground station.
* The motion of the spacecraft during the signal flight time is taken into
......@@ -53,7 +52,10 @@ public class Angular extends AbstractMeasurement<Angular> {
public Angular(final GroundStation station, final AbsoluteDate date,
final double[] angular, final double[] sigma, final double[] baseWeight)
throws OrekitException {
super(date, angular, sigma, baseWeight, station.getPositionOffsetDriver());
super(date, angular, sigma, baseWeight,
station.getEastOffsetDriver(),
station.getNorthOffsetDriver(),
station.getZenithOffsetDriver());
this.station = station;
}
......@@ -150,28 +152,30 @@ public class Angular extends AbstractMeasurement<Angular> {
evaluation.setStateDerivatives(dAzOndP, dElOndP);
final ParameterDriver positionOffsetDriver = station.getPositionOffsetDriver();
if (positionOffsetDriver.isEstimated()) {
if (station.getEastOffsetDriver().isSelected() ||
station.getNorthOffsetDriver().isSelected() ||
station.getZenithOffsetDriver().isSelected()) {
// partial derivatives with respect to parameters
// Be aware: east; north and zenith are expressed in station parent frame but the derivatives are expressed
// with respect to reference station topocentric frame
// partial derivatives of azimuth with respect to parameters in body frame
final double[] dAzOndQ = new double[] {
if (station.getEastOffsetDriver().isSelected()) {
evaluation.setParameterDerivatives(station.getEastOffsetDriver(),
azimuth.getPartialDerivative(0, 0, 0, 1, 0, 0),
elevation.getPartialDerivative(0, 0, 0, 1, 0, 0));
}
if (station.getNorthOffsetDriver().isSelected()) {
evaluation.setParameterDerivatives(station.getNorthOffsetDriver(),
azimuth.getPartialDerivative(0, 0, 0, 0, 1, 0),
azimuth.getPartialDerivative(0, 0, 0, 0, 0, 1)
};
elevation.getPartialDerivative(0, 0, 0, 0, 1, 0));
}
if (station.getZenithOffsetDriver().isSelected()) {
evaluation.setParameterDerivatives(station.getZenithOffsetDriver(),
azimuth.getPartialDerivative(0, 0, 0, 0, 0, 1),
elevation.getPartialDerivative(0, 0, 0, 0, 0, 1));
}
// partial derivatives of elevation with respect to parameters in body frame
final double[] dElOndQ = new double[] {
elevation.getPartialDerivative(0, 0, 0, 1, 0, 0),
elevation.getPartialDerivative(0, 0, 0, 0, 1, 0),
elevation.getPartialDerivative(0, 0, 0, 0, 0, 1)
};
evaluation.setParameterDerivatives(positionOffsetDriver, dAzOndQ, dElOndQ);
}
return evaluation;
......
......@@ -16,7 +16,8 @@
*/
package org.orekit.estimation.measurements;
import java.util.Arrays;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import org.orekit.errors.OrekitException;
......@@ -29,46 +30,48 @@ import org.orekit.utils.ParameterDriver;
*/
public class Bias<T extends Measurement<T>> implements EvaluationModifier<T> {
/** Parameter holding the bias value. */
private final ParameterDriver driver;
/** Parameters holding the bias value components. */
private final List<ParameterDriver> drivers;
/** Identity matrix, for partial derivatives. */
/** IPartial derivatives. */
private final double[][] derivatives;
/** Simple constructor.
* @param name name of the bias
* @param bias initial value of the bias
* @param scale scale of the bias, for normalization
* @param min minimum value of the bias
* @param max maximum value of the bias
* @exception OrekitException if initial value cannot be set
*/
public Bias(final String name, final double ... bias)
public Bias(final String[] name, final double[] bias, final double[] scale,
final double[] min, final double[] max)
throws OrekitException {
driver = new ParameterDriver(name, bias) {
/** {@inheritDoc} */
@Override
public void valueChanged(final double[] newValue) {
}
};
drivers = new ArrayList<>(bias.length);
for (int i = 0; i < bias.length; ++i) {
drivers.add(new ParameterDriver(name[i], bias[i], scale[i], min[i], max[i]));
}
derivatives = new double[bias.length][bias.length];
for (int i = 0; i < bias.length; ++i) {
// derivatives are computed with respect to the physical parameters,
// not with respect to the normalized parameters (normalization is
// performed later on), so the derivative is really 1.0 and not scale[i]
derivatives[i][i] = 1.0;
}
}
/** Get the bias driver.
* @return driver for the bias
/** {@inheritDoc}
* <p>
* For a bias, there are {@link Measurement#getDimension()} parameter drivers,
* sorted in components order.
* </p>
*/
public ParameterDriver getDriver() {
return driver;
}
/** {@inheritDoc} */
@Override
public List<ParameterDriver> getParametersDrivers() {
return Arrays.asList(driver);
return Collections.unmodifiableList(drivers);
}
/** {@inheritDoc} */
......@@ -77,16 +80,16 @@ public class Bias<T extends Measurement<T>> implements EvaluationModifier<T> {
// apply the bias to the measurement value
final double[] measurementValue = evaluation.getValue();
final double[] biasValue = driver.getValue();
for (int i = 0; i < driver.getDimension(); ++i) {
measurementValue[i] += biasValue[i];
for (int i = 0; i < drivers.size(); ++i) {
final ParameterDriver driver = drivers.get(i);
measurementValue[i] += driver.getValue();
if (driver.isSelected()) {
// add the partial derivatives
evaluation.setParameterDerivatives(driver, derivatives[i]);
}
}
evaluation.setValue(measurementValue);
if (driver.isEstimated()) {
// add the partial derivatives
evaluation.setParameterDerivatives(driver, derivatives);
}
}
......
......@@ -55,7 +55,7 @@ public class Evaluation<T extends Measurement<T>> implements TimeStamped {
private double[][] stateDerivatives;
/** Partial derivatives with respect to parameters. */
private final Map<ParameterDriver, double[][]> parametersDerivatives;
private final Map<ParameterDriver, double[]> parametersDerivatives;
/** Simple constructor.
* @param measurement associated measurement
......@@ -70,7 +70,7 @@ public class Evaluation<T extends Measurement<T>> implements TimeStamped {
this.iteration = iteration;
this.count = count;
this.state = state;
this.parametersDerivatives = new IdentityHashMap<ParameterDriver, double[][]>();
this.parametersDerivatives = new IdentityHashMap<ParameterDriver, double[]>();
}
/** Get the associated measurement.
......@@ -176,12 +176,12 @@ public class Evaluation<T extends Measurement<T>> implements TimeStamped {
* @return partial derivatives of the simulated value
* @exception OrekitIllegalArgumentException if parameter is unknown
*/
public double[][] getParameterDerivatives(final ParameterDriver driver)
public double[] getParameterDerivatives(final ParameterDriver driver)
throws OrekitIllegalArgumentException {
final double[][] p = parametersDerivatives.get(driver);
final double[] p = parametersDerivatives.get(driver);
if (p == null) {
final StringBuilder builder = new StringBuilder();
for (final Map.Entry<ParameterDriver, double[][]> entry : parametersDerivatives.entrySet()) {
for (final Map.Entry<ParameterDriver, double[]> entry : parametersDerivatives.entrySet()) {
if (builder.length() > 0) {
builder.append(", ");
}
......@@ -191,24 +191,16 @@ public class Evaluation<T extends Measurement<T>> implements TimeStamped {
driver.getName(),
builder.length() > 0 ? builder.toString() : "<none>");
}
final double[][] sd = new double[measurement.getDimension()][];
for (int i = 0; i < measurement.getDimension(); ++i) {
sd[i] = p[i].clone();
}
return sd;
return p;
}
/** Set the partial derivatives of the {@link #getValue()
* simulated measurement} with respect to state Cartesian coordinates.
* simulated measurement} with respect to parameter.
* @param driver driver for the parameter
* @param parameterDerivatives partial derivatives with respect to parameter
*/
public void setParameterDerivatives(final ParameterDriver driver, final double[] ... parameterDerivatives) {
final double[][] p = new double[measurement.getDimension()][];
for (int i = 0; i < measurement.getDimension(); ++i) {
p[i] = parameterDerivatives[i].clone();
}
parametersDerivatives.put(driver, p);