use of org.orekit.utils.ParameterDriversList in project Orekit by CS-SI.
the class BatchLSEstimatorTest method testKeplerRange.
/**
* Perfect range measurements with a biased start
* @throws OrekitException
*/
@Test
public void testKeplerRange() throws OrekitException {
Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 1.0);
// create perfect range measurements
final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new RangeMeasurementCreator(context), 1.0, 3.0, 300.0);
// create orbit estimator
final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder);
for (final ObservedMeasurement<?> range : measurements) {
estimator.addMeasurement(range);
}
estimator.setParametersConvergenceThreshold(1.0e-2);
estimator.setMaxIterations(10);
estimator.setMaxEvaluations(20);
estimator.setObserver(new BatchLSObserver() {
int lastIter = 0;
int lastEval = 0;
/**
* {@inheritDoc}
*/
@Override
public void evaluationPerformed(int iterationsCount, int evaluationscount, Orbit[] orbits, ParameterDriversList estimatedOrbitalParameters, ParameterDriversList estimatedPropagatorParameters, ParameterDriversList estimatedMeasurementsParameters, EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) throws OrekitException {
if (iterationsCount == lastIter) {
Assert.assertEquals(lastEval + 1, evaluationscount);
} else {
Assert.assertEquals(lastIter + 1, iterationsCount);
}
lastIter = iterationsCount;
lastEval = evaluationscount;
Assert.assertEquals(measurements.size(), evaluationsProvider.getNumber());
try {
evaluationsProvider.getEstimatedMeasurement(-1);
Assert.fail("an exception should have been thrown");
} catch (OrekitException oe) {
Assert.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
}
try {
evaluationsProvider.getEstimatedMeasurement(measurements.size());
Assert.fail("an exception should have been thrown");
} catch (OrekitException oe) {
Assert.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
}
AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
Assert.assertTrue(current.compareTo(previous) >= 0);
previous = current;
}
}
});
ParameterDriver aDriver = estimator.getOrbitalParametersDrivers(true).getDrivers().get(0);
Assert.assertEquals("a", aDriver.getName());
aDriver.setValue(aDriver.getValue() + 1.2);
aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
EstimationTestUtils.checkFit(context, estimator, 2, 3, 0.0, 1.1e-6, 0.0, 2.8e-6, 0.0, 4.0e-7, 0.0, 2.2e-10);
// got a default one
for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
if ("a".equals(driver.getName())) {
// user-specified reference date
Assert.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
} else {
// default reference date
Assert.assertEquals(0, driver.getReferenceDate().durationFrom(propagatorBuilder.getInitialOrbitDate()), 1.0e-15);
}
}
}
use of org.orekit.utils.ParameterDriversList in project Orekit by CS-SI.
the class ModelTest method testPerfectValue.
@Test
public void testPerfectValue() throws OrekitException {
final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
final NumericalPropagatorBuilder[] builders = { propagatorBuilder };
// create perfect PV measurements
final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new PVMeasurementCreator(), 0.0, 1.0, 300.0);
final ParameterDriversList estimatedMeasurementsParameters = new ParameterDriversList();
for (ObservedMeasurement<?> measurement : measurements) {
for (final ParameterDriver driver : measurement.getParametersDrivers()) {
if (driver.isSelected()) {
estimatedMeasurementsParameters.add(driver);
}
}
}
// create model
final ModelObserver modelObserver = new ModelObserver() {
/**
* {@inheritDoc}
*/
@Override
public void modelCalled(final Orbit[] newOrbits, final Map<ObservedMeasurement<?>, EstimatedMeasurement<?>> newEvaluations) {
Assert.assertEquals(1, newOrbits.length);
Assert.assertEquals(0, context.initialOrbit.getDate().durationFrom(newOrbits[0].getDate()), 1.0e-15);
Assert.assertEquals(0, Vector3D.distance(context.initialOrbit.getPVCoordinates().getPosition(), newOrbits[0].getPVCoordinates().getPosition()), 1.0e-15);
Assert.assertEquals(measurements.size(), newEvaluations.size());
}
};
final Model model = new Model(builders, measurements, estimatedMeasurementsParameters, modelObserver);
model.setIterationsCounter(new Incrementor(100));
model.setEvaluationsCounter(new Incrementor(100));
// Test forward propagation flag to true
assertEquals(true, model.isForwardPropagation());
// evaluate model on perfect start point
final double[] normalizedProp = propagatorBuilder.getSelectedNormalizedParameters();
final double[] normalized = new double[normalizedProp.length + estimatedMeasurementsParameters.getNbParams()];
System.arraycopy(normalizedProp, 0, normalized, 0, normalizedProp.length);
int i = normalizedProp.length;
for (final ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
normalized[i++] = driver.getNormalizedValue();
}
Pair<RealVector, RealMatrix> value = model.value(new ArrayRealVector(normalized));
int index = 0;
for (ObservedMeasurement<?> measurement : measurements) {
for (int k = 0; k < measurement.getDimension(); ++k) {
// the value is already a weighted residual
Assert.assertEquals(0.0, value.getFirst().getEntry(index++), 1.6e-7);
}
}
Assert.assertEquals(index, value.getFirst().getDimension());
}
use of org.orekit.utils.ParameterDriversList in project Orekit by CS-SI.
the class KalmanEstimatorTest method testEquinoctialRightAscensionDeclination.
/**
* Perfect right-ascension/declination measurements with a perfect start
* Equinoctial formalism
* @throws OrekitException
*/
@Test
public void testEquinoctialRightAscensionDeclination() throws OrekitException {
// Create context
Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
// Create initial orbit and propagator builder
final OrbitType orbitType = OrbitType.EQUINOCTIAL;
final PositionAngle positionAngle = PositionAngle.TRUE;
final boolean perfectStart = true;
final double minStep = 1.e-6;
final double maxStep = 60.;
final double dP = 1.;
final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(orbitType, positionAngle, perfectStart, minStep, maxStep, dP);
// Create perfect range measurements
final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new AngularRaDecMeasurementCreator(context), 1.0, 4.0, 60.0);
// Reference propagator for estimation performances
final NumericalPropagator referencePropagator = propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters());
// Reference position/velocity at last measurement date
final Orbit refOrbit = referencePropagator.propagate(measurements.get(measurements.size() - 1).getDate()).getOrbit();
// Cartesian covariance matrix initialization
final RealMatrix cartesianP = MatrixUtils.createRealDiagonalMatrix(new double[] { 1e-4, 1e-4, 1e-4, 1e-10, 1e-10, 1e-10 });
// Jacobian of the orbital parameters w/r to Cartesian
final Orbit initialOrbit = orbitType.convertType(context.initialOrbit);
final double[][] dYdC = new double[6][6];
initialOrbit.getJacobianWrtCartesian(positionAngle, dYdC);
final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
// Keplerian initial covariance matrix
final RealMatrix initialP = Jac.multiply(cartesianP.multiply(Jac.transpose()));
// Process noise matrix
final RealMatrix cartesianQ = MatrixUtils.createRealDiagonalMatrix(new double[] { 1.e-6, 1.e-6, 1.e-6, 1.e-12, 1.e-12, 1.e-12 });
final RealMatrix Q = Jac.multiply(cartesianQ.multiply(Jac.transpose()));
// Build the Kalman filter
final KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
kalmanBuilder.builder(propagatorBuilder);
kalmanBuilder.estimatedMeasurementsParameters(new ParameterDriversList());
kalmanBuilder.initialCovarianceMatrix(initialP);
kalmanBuilder.processNoiseMatrixProvider(new ConstantProcessNoise(Q));
final KalmanEstimator kalman = kalmanBuilder.build();
// Filter the measurements and check the results
final double expectedDeltaPos = 0.;
final double posEps = 1.53e-5;
final double expectedDeltaVel = 0.;
final double velEps = 5.04e-9;
final double[] expectedSigmasPos = { 0.356902, 1.297507, 1.798551 };
final double sigmaPosEps = 1e-6;
final double[] expectedSigmasVel = { 2.468745e-4, 5.810027e-4, 3.887394e-4 };
final double sigmaVelEps = 1e-10;
EstimationTestUtils.checkKalmanFit(context, kalman, measurements, refOrbit, positionAngle, expectedDeltaPos, posEps, expectedDeltaVel, velEps, expectedSigmasPos, sigmaPosEps, expectedSigmasVel, sigmaVelEps);
}
use of org.orekit.utils.ParameterDriversList in project Orekit by CS-SI.
the class PartialDerivativesTest method doTestParametersDerivatives.
private void doTestParametersDerivatives(String parameterName, double tolerance, OrbitType... orbitTypes) throws OrekitException {
OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
ForceModel drag = new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(), earth), new IsotropicDrag(2.5, 1.2));
NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
ForceModel gravityField = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
Orbit baseOrbit = new KeplerianOrbit(7000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.TRUE, FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
double dt = 900;
double dP = 1.0;
for (OrbitType orbitType : orbitTypes) {
final Orbit initialOrbit = orbitType.convertType(baseOrbit);
for (PositionAngle angleType : PositionAngle.values()) {
NumericalPropagator propagator = setUpPropagator(initialOrbit, dP, orbitType, angleType, drag, gravityField);
propagator.setMu(provider.getMu());
for (final ForceModel forceModel : propagator.getAllForceModels()) {
for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
driver.setValue(driver.getReferenceValue());
driver.setSelected(driver.getName().equals(parameterName));
}
}
PartialDerivativesEquations partials = new PartialDerivativesEquations("partials", propagator);
final SpacecraftState initialState = partials.setInitialJacobians(new SpacecraftState(initialOrbit));
propagator.setInitialState(initialState);
final JacobiansMapper mapper = partials.getMapper();
PickUpHandler pickUp = new PickUpHandler(mapper, null);
propagator.setMasterMode(pickUp);
propagator.propagate(initialState.getDate().shiftedBy(dt));
double[][] dYdP = pickUp.getdYdP();
// compute reference Jacobian using finite differences
double[][] dYdPRef = new double[6][1];
NumericalPropagator propagator2 = setUpPropagator(initialOrbit, dP, orbitType, angleType, drag, gravityField);
propagator2.setMu(provider.getMu());
ParameterDriversList bound = new ParameterDriversList();
for (final ForceModel forceModel : propagator2.getAllForceModels()) {
for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
if (driver.getName().equals(parameterName)) {
driver.setSelected(true);
bound.add(driver);
} else {
driver.setSelected(false);
}
}
}
ParameterDriver selected = bound.getDrivers().get(0);
double p0 = selected.getReferenceValue();
double h = selected.getScale();
selected.setValue(p0 - 4 * h);
propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
propagator2.getMu(), initialState.getAttitude()));
SpacecraftState sM4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
selected.setValue(p0 - 3 * h);
propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
propagator2.getMu(), initialState.getAttitude()));
SpacecraftState sM3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
selected.setValue(p0 - 2 * h);
propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
propagator2.getMu(), initialState.getAttitude()));
SpacecraftState sM2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
selected.setValue(p0 - 1 * h);
propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
propagator2.getMu(), initialState.getAttitude()));
SpacecraftState sM1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
selected.setValue(p0 + 1 * h);
propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
propagator2.getMu(), initialState.getAttitude()));
SpacecraftState sP1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
selected.setValue(p0 + 2 * h);
propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
propagator2.getMu(), initialState.getAttitude()));
SpacecraftState sP2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
selected.setValue(p0 + 3 * h);
propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
propagator2.getMu(), initialState.getAttitude()));
SpacecraftState sP3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
selected.setValue(p0 + 4 * h);
propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
propagator2.getMu(), initialState.getAttitude()));
SpacecraftState sP4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
fillJacobianColumn(dYdPRef, 0, orbitType, angleType, h, sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
for (int i = 0; i < 6; ++i) {
Assert.assertEquals(dYdPRef[i][0], dYdP[i][0], FastMath.abs(dYdPRef[i][0] * tolerance));
}
}
}
}
use of org.orekit.utils.ParameterDriversList in project Orekit by CS-SI.
the class Model method getSelectedPropagationDriversForBuilder.
/**
* Get the selected propagation drivers for a propagatorBuilder.
* @param iBuilder index of the builder in the builders' array
* @return the list of selected propagation drivers for propagatorBuilder of index iBuilder
* @exception OrekitException if orbit cannot be created with the current point
*/
public ParameterDriversList getSelectedPropagationDriversForBuilder(final int iBuilder) throws OrekitException {
// Lazy evaluation, create the list only if it hasn't been created yet
if (estimatedPropagationParameters[iBuilder] == null) {
// Gather the drivers
final ParameterDriversList selectedPropagationDrivers = new ParameterDriversList();
for (final DelegatingDriver delegating : builders[iBuilder].getPropagationParametersDrivers().getDrivers()) {
if (delegating.isSelected()) {
for (final ParameterDriver driver : delegating.getRawDrivers()) {
selectedPropagationDrivers.add(driver);
}
}
}
// List of propagation drivers are sorted in the BatchLSEstimator class.
// Hence we need to sort this list so the parameters' indexes match
selectedPropagationDrivers.sort();
// Add the list of selected propagation drivers to the array
estimatedPropagationParameters[iBuilder] = selectedPropagationDrivers;
}
return estimatedPropagationParameters[iBuilder];
}
Aggregations