use of org.orekit.estimation.measurements.ObservedMeasurement in project Orekit by CS-SI.
the class OrbitDetermination method run.
private void run(final File input, final File home) throws IOException, IllegalArgumentException, OrekitException, ParseException {
// read input parameters
KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class);
try (final FileInputStream fis = new FileInputStream(input)) {
parser.parseInput(input.getAbsolutePath(), fis);
}
// log file
final String baseName;
final PrintStream logStream;
if (parser.containsKey(ParameterKey.OUTPUT_BASE_NAME) && parser.getString(ParameterKey.OUTPUT_BASE_NAME).length() > 0) {
baseName = parser.getString(ParameterKey.OUTPUT_BASE_NAME);
logStream = new PrintStream(new File(home, baseName + "-log.out"), "UTF-8");
} else {
baseName = null;
logStream = null;
}
final RangeLog rangeLog = new RangeLog(home, baseName);
final RangeRateLog rangeRateLog = new RangeRateLog(home, baseName);
final AzimuthLog azimuthLog = new AzimuthLog(home, baseName);
final ElevationLog elevationLog = new ElevationLog(home, baseName);
final PositionLog positionLog = new PositionLog(home, baseName);
final VelocityLog velocityLog = new VelocityLog(home, baseName);
try {
// gravity field
final NormalizedSphericalHarmonicsProvider gravityField = createGravityField(parser);
// Orbit initial guess
final Orbit initialGuess = createOrbit(parser, gravityField.getMu());
// IERS conventions
final IERSConventions conventions;
if (!parser.containsKey(ParameterKey.IERS_CONVENTIONS)) {
conventions = IERSConventions.IERS_2010;
} else {
conventions = IERSConventions.valueOf("IERS_" + parser.getInt(ParameterKey.IERS_CONVENTIONS));
}
// central body
final OneAxisEllipsoid body = createBody(parser);
// propagator builder
final NumericalPropagatorBuilder propagatorBuilder = createPropagatorBuilder(parser, conventions, gravityField, body, initialGuess);
// estimator
final BatchLSEstimator estimator = createEstimator(parser, propagatorBuilder);
// measurements
final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
for (final String fileName : parser.getStringsList(ParameterKey.MEASUREMENTS_FILES, ',')) {
measurements.addAll(readMeasurements(new File(input.getParentFile(), fileName), createStationsData(parser, body), createPVData(parser), createSatRangeBias(parser), createWeights(parser), createRangeOutliersManager(parser), createRangeRateOutliersManager(parser), createAzElOutliersManager(parser), createPVOutliersManager(parser)));
}
for (ObservedMeasurement<?> measurement : measurements) {
estimator.addMeasurement(measurement);
}
// estimate orbit
estimator.setObserver(new BatchLSObserver() {
private PVCoordinates previousPV;
{
previousPV = initialGuess.getPVCoordinates();
final String header = "iteration evaluations ΔP(m) ΔV(m/s) RMS nb Range nb Range-rate nb Angular nb PV%n";
System.out.format(Locale.US, header);
if (logStream != null) {
logStream.format(Locale.US, header);
}
}
/**
* {@inheritDoc}
*/
@Override
public void evaluationPerformed(final int iterationsCount, final int evaluationsCount, final Orbit[] orbits, final ParameterDriversList estimatedOrbitalParameters, final ParameterDriversList estimatedPropagatorParameters, final ParameterDriversList estimatedMeasurementsParameters, final EstimationsProvider evaluationsProvider, final LeastSquaresProblem.Evaluation lspEvaluation) {
PVCoordinates currentPV = orbits[0].getPVCoordinates();
final String format0 = " %2d %2d %16.12f %s %s %s %s%n";
final String format = " %2d %2d %13.6f %12.9f %16.12f %s %s %s %s%n";
final EvaluationCounter<Range> rangeCounter = new EvaluationCounter<Range>();
final EvaluationCounter<RangeRate> rangeRateCounter = new EvaluationCounter<RangeRate>();
final EvaluationCounter<AngularAzEl> angularCounter = new EvaluationCounter<AngularAzEl>();
final EvaluationCounter<PV> pvCounter = new EvaluationCounter<PV>();
for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry : estimator.getLastEstimations().entrySet()) {
if (entry.getKey() instanceof Range) {
@SuppressWarnings("unchecked") EstimatedMeasurement<Range> evaluation = (EstimatedMeasurement<Range>) entry.getValue();
rangeCounter.add(evaluation);
} else if (entry.getKey() instanceof RangeRate) {
@SuppressWarnings("unchecked") EstimatedMeasurement<RangeRate> evaluation = (EstimatedMeasurement<RangeRate>) entry.getValue();
rangeRateCounter.add(evaluation);
} else if (entry.getKey() instanceof AngularAzEl) {
@SuppressWarnings("unchecked") EstimatedMeasurement<AngularAzEl> evaluation = (EstimatedMeasurement<AngularAzEl>) entry.getValue();
angularCounter.add(evaluation);
} else if (entry.getKey() instanceof PV) {
@SuppressWarnings("unchecked") EstimatedMeasurement<PV> evaluation = (EstimatedMeasurement<PV>) entry.getValue();
pvCounter.add(evaluation);
}
}
if (evaluationsCount == 1) {
System.out.format(Locale.US, format0, iterationsCount, evaluationsCount, lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
if (logStream != null) {
logStream.format(Locale.US, format0, iterationsCount, evaluationsCount, lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
}
} else {
System.out.format(Locale.US, format, iterationsCount, evaluationsCount, Vector3D.distance(previousPV.getPosition(), currentPV.getPosition()), Vector3D.distance(previousPV.getVelocity(), currentPV.getVelocity()), lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
if (logStream != null) {
logStream.format(Locale.US, format, iterationsCount, evaluationsCount, Vector3D.distance(previousPV.getPosition(), currentPV.getPosition()), Vector3D.distance(previousPV.getVelocity(), currentPV.getVelocity()), lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
}
}
previousPV = currentPV;
}
});
Orbit estimated = estimator.estimate()[0].getInitialState().getOrbit();
// compute some statistics
for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry : estimator.getLastEstimations().entrySet()) {
if (entry.getKey() instanceof Range) {
@SuppressWarnings("unchecked") EstimatedMeasurement<Range> evaluation = (EstimatedMeasurement<Range>) entry.getValue();
rangeLog.add(evaluation);
} else if (entry.getKey() instanceof RangeRate) {
@SuppressWarnings("unchecked") EstimatedMeasurement<RangeRate> evaluation = (EstimatedMeasurement<RangeRate>) entry.getValue();
rangeRateLog.add(evaluation);
} else if (entry.getKey() instanceof AngularAzEl) {
@SuppressWarnings("unchecked") EstimatedMeasurement<AngularAzEl> evaluation = (EstimatedMeasurement<AngularAzEl>) entry.getValue();
azimuthLog.add(evaluation);
elevationLog.add(evaluation);
} else if (entry.getKey() instanceof PV) {
@SuppressWarnings("unchecked") EstimatedMeasurement<PV> evaluation = (EstimatedMeasurement<PV>) entry.getValue();
positionLog.add(evaluation);
velocityLog.add(evaluation);
}
}
System.out.println("Estimated orbit: " + estimated);
if (logStream != null) {
logStream.println("Estimated orbit: " + estimated);
}
final ParameterDriversList orbitalParameters = estimator.getOrbitalParametersDrivers(true);
final ParameterDriversList propagatorParameters = estimator.getPropagatorParametersDrivers(true);
final ParameterDriversList measurementsParameters = estimator.getMeasurementsParametersDrivers(true);
int length = 0;
for (final ParameterDriver parameterDriver : orbitalParameters.getDrivers()) {
length = FastMath.max(length, parameterDriver.getName().length());
}
for (final ParameterDriver parameterDriver : propagatorParameters.getDrivers()) {
length = FastMath.max(length, parameterDriver.getName().length());
}
for (final ParameterDriver parameterDriver : measurementsParameters.getDrivers()) {
length = FastMath.max(length, parameterDriver.getName().length());
}
displayParametersChanges(System.out, "Estimated orbital parameters changes: ", false, length, orbitalParameters);
if (logStream != null) {
displayParametersChanges(logStream, "Estimated orbital parameters changes: ", false, length, orbitalParameters);
}
displayParametersChanges(System.out, "Estimated propagator parameters changes: ", true, length, propagatorParameters);
if (logStream != null) {
displayParametersChanges(logStream, "Estimated propagator parameters changes: ", true, length, propagatorParameters);
}
displayParametersChanges(System.out, "Estimated measurements parameters changes: ", true, length, measurementsParameters);
if (logStream != null) {
displayParametersChanges(logStream, "Estimated measurements parameters changes: ", true, length, measurementsParameters);
}
System.out.println("Number of iterations: " + estimator.getIterationsCount());
System.out.println("Number of evaluations: " + estimator.getEvaluationsCount());
rangeLog.displaySummary(System.out);
rangeRateLog.displaySummary(System.out);
azimuthLog.displaySummary(System.out);
elevationLog.displaySummary(System.out);
positionLog.displaySummary(System.out);
velocityLog.displaySummary(System.out);
if (logStream != null) {
logStream.println("Number of iterations: " + estimator.getIterationsCount());
logStream.println("Number of evaluations: " + estimator.getEvaluationsCount());
rangeLog.displaySummary(logStream);
rangeRateLog.displaySummary(logStream);
azimuthLog.displaySummary(logStream);
elevationLog.displaySummary(logStream);
positionLog.displaySummary(logStream);
velocityLog.displaySummary(logStream);
}
rangeLog.displayResiduals();
rangeRateLog.displayResiduals();
azimuthLog.displayResiduals();
elevationLog.displayResiduals();
positionLog.displayResiduals();
velocityLog.displayResiduals();
} finally {
if (logStream != null) {
logStream.close();
}
rangeLog.close();
rangeRateLog.close();
azimuthLog.close();
elevationLog.close();
positionLog.close();
velocityLog.close();
}
}
use of org.orekit.estimation.measurements.ObservedMeasurement in project Orekit by CS-SI.
the class KalmanEstimatorTest method testCartesianRangeRate.
/**
* Perfect range rate measurements with a perfect start
* Cartesian formalism
* @throws OrekitException
*/
@Test
public void testCartesianRangeRate() throws OrekitException {
// Create context
Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
// Create initial orbit and propagator builder
final OrbitType orbitType = OrbitType.CARTESIAN;
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 RangeRateMeasurementCreator(context, false), 1.0, 3.0, 300.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
// 100m on position / 1e-2m/s on velocity
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.TRUE, dYdC);
final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
// 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 = 9.50e-4;
final double expectedDeltaVel = 0.;
final double velEps = 3.49e-7;
final double[] expectedSigmasPos = { 0.324398, 1.347031, 1.743310 };
final double sigmaPosEps = 1e-6;
final double[] expectedSigmasVel = { 2.856883e-4, 5.765844e-4, 5.056186e-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.estimation.measurements.ObservedMeasurement in project Orekit by CS-SI.
the class KalmanEstimatorTest method testCircularAzimuthElevation.
/**
* Perfect azimuth/elevation measurements with a perfect start
* Circular formalism
* @throws OrekitException
*/
@Test
public void testCircularAzimuthElevation() throws OrekitException {
// Create context
Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
// Create initial orbit and propagator builder
final OrbitType orbitType = OrbitType.CIRCULAR;
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 AngularAzElMeasurementCreator(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.TRUE, dYdC);
final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
// 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 = 4.78e-7;
final double expectedDeltaVel = 0.;
final double velEps = 1.54e-10;
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.estimation.measurements.ObservedMeasurement in project Orekit by CS-SI.
the class KalmanEstimatorTest method testKeplerianRange.
/**
* Perfect range measurements with a biased start
* Keplerian formalism
* @throws OrekitException
*/
@Test
public void testKeplerianRange() throws OrekitException {
// Create context
Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
// Create initial orbit and propagator builder
final OrbitType orbitType = OrbitType.KEPLERIAN;
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 RangeMeasurementCreator(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();
// Change semi-major axis of 1.2m as in the batch test
ParameterDriver aDriver = propagatorBuilder.getOrbitalParametersDrivers().getDrivers().get(0);
aDriver.setValue(aDriver.getValue() + 1.2);
aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
// Cartesian covariance matrix initialization
// 100m on position / 1e-2m/s on velocity
final RealMatrix cartesianP = MatrixUtils.createRealDiagonalMatrix(new double[] { 100., 100., 100., 1e-2, 1e-2, 1e-2 });
// 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.TRUE, dYdC);
final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
// Keplerian initial covariance matrix
final RealMatrix initialP = Jac.multiply(cartesianP.multiply(Jac.transpose()));
// Process noise matrix is set to 0 here
RealMatrix Q = MatrixUtils.createRealMatrix(6, 6);
// 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.77e-4;
final double expectedDeltaVel = 0.;
final double velEps = 7.93e-8;
final double[] expectedSigmasPos = { 0.742488, 0.281910, 0.563217 };
final double sigmaPosEps = 1e-6;
final double[] expectedSigmasVel = { 2.206622e-4, 1.306669e-4, 1.293996e-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.estimation.measurements.ObservedMeasurement in project Orekit by CS-SI.
the class KalmanEstimatorTest method testWrappedException.
/**
* Test of a wrapped exception in a Kalman observer
* @throws OrekitException
*/
@Test
public void testWrappedException() throws OrekitException {
// Create context
Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
// Create initial orbit and propagator builder
final OrbitType orbitType = OrbitType.KEPLERIAN;
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 RangeMeasurementCreator(context), 1.0, 3.0, 300.0);
// Build the Kalman filter
final KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
kalmanBuilder.builder(propagatorBuilder);
kalmanBuilder.estimatedMeasurementsParameters(new ParameterDriversList());
kalmanBuilder.initialCovarianceMatrix(MatrixUtils.createRealMatrix(6, 6));
kalmanBuilder.processNoiseMatrixProvider(new ConstantProcessNoise(MatrixUtils.createRealMatrix(6, 6)));
final KalmanEstimator kalman = kalmanBuilder.build();
kalman.setObserver(estimation -> {
throw new DummyException();
});
try {
// Filter the measurements and expect an exception to occur
EstimationTestUtils.checkKalmanFit(context, kalman, measurements, context.initialOrbit, positionAngle, 0., 0., 0., 0., new double[3], 0., new double[3], 0.);
} catch (DummyException de) {
// expected
}
}
Aggregations