Examples with BatchLSEstimator org.orekit.estimation.leastsquares.BatchLSEstimator used on opensource projects

Example 1 with BatchLSEstimator

use of org.orekit.estimation.leastsquares.BatchLSEstimator in project Orekit by CS-SI.

the class HarmonicParametricAccelerationTest method testCoefficientsDetermination.

@Test
public void testCoefficientsDetermination() throws OrekitException {
    final double mass = 2500;
    final Orbit orbit = new CircularOrbit(7500000.0, 1.0e-4, 1.0e-3, 1.7, 0.3, 0.5, PositionAngle.TRUE, FramesFactory.getEME2000(), new AbsoluteDate(new DateComponents(2004, 2, 3), TimeComponents.H00, TimeScalesFactory.getUTC()), Constants.EIGEN5C_EARTH_MU);
    final double period = orbit.getKeplerianPeriod();
    AttitudeProvider maneuverLaw = new LofOffset(orbit.getFrame(), LOFType.VNC);
    SpacecraftState initialState = new SpacecraftState(orbit, maneuverLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
    double dP = 10.0;
    double minStep = 0.001;
    double maxStep = 100;
    double[][] tolerance = NumericalPropagator.tolerances(dP, orbit, orbit.getType());
    // generate PV measurements corresponding to a tangential maneuver
    AdaptiveStepsizeIntegrator integrator0 = new DormandPrince853Integrator(minStep, maxStep, tolerance[0], tolerance[1]);
    integrator0.setInitialStepSize(60);
    final NumericalPropagator propagator0 = new NumericalPropagator(integrator0);
    propagator0.setInitialState(initialState);
    propagator0.setAttitudeProvider(maneuverLaw);
    ForceModel hpaRefX1 = new HarmonicParametricAcceleration(Vector3D.PLUS_I, true, "refX1", null, period, 1);
    ForceModel hpaRefY1 = new HarmonicParametricAcceleration(Vector3D.PLUS_J, true, "refY1", null, period, 1);
    ForceModel hpaRefZ2 = new HarmonicParametricAcceleration(Vector3D.PLUS_K, true, "refZ2", null, period, 2);
    hpaRefX1.getParametersDrivers()[0].setValue(2.4e-2);
    hpaRefX1.getParametersDrivers()[1].setValue(3.1);
    hpaRefY1.getParametersDrivers()[0].setValue(4.0e-2);
    hpaRefY1.getParametersDrivers()[1].setValue(0.3);
    hpaRefZ2.getParametersDrivers()[0].setValue(1.0e-2);
    hpaRefZ2.getParametersDrivers()[1].setValue(1.8);
    propagator0.addForceModel(hpaRefX1);
    propagator0.addForceModel(hpaRefY1);
    propagator0.addForceModel(hpaRefZ2);
    final List<ObservedMeasurement<?>> measurements = new ArrayList<>();
    propagator0.setMasterMode(10.0, (state, isLast) -> measurements.add(new PV(state.getDate(), state.getPVCoordinates().getPosition(), state.getPVCoordinates().getVelocity(), 1.0e-3, 1.0e-6, 1.0)));
    propagator0.propagate(orbit.getDate().shiftedBy(900));
    // set up an estimator to retrieve the maneuver as several harmonic accelerations in inertial frame
    final NumericalPropagatorBuilder propagatorBuilder = new NumericalPropagatorBuilder(orbit, new DormandPrince853IntegratorBuilder(minStep, maxStep, dP), PositionAngle.TRUE, dP);
    propagatorBuilder.addForceModel(new HarmonicParametricAcceleration(Vector3D.PLUS_I, true, "X1", null, period, 1));
    propagatorBuilder.addForceModel(new HarmonicParametricAcceleration(Vector3D.PLUS_J, true, "Y1", null, period, 1));
    propagatorBuilder.addForceModel(new HarmonicParametricAcceleration(Vector3D.PLUS_K, true, "Z2", null, period, 2));
    final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder);
    estimator.setParametersConvergenceThreshold(1.0e-2);
    estimator.setMaxIterations(20);
    estimator.setMaxEvaluations(100);
    for (final ObservedMeasurement<?> measurement : measurements) {
        estimator.addMeasurement(measurement);
    }
    // we will estimate only the force model parameters, not the orbit
    for (final ParameterDriver d : estimator.getOrbitalParametersDrivers(false).getDrivers()) {
        d.setSelected(false);
    }
    setParameter(estimator, "X1 γ", 1.0e-2);
    setParameter(estimator, "X1 φ", 4.0);
    setParameter(estimator, "Y1 γ", 1.0e-2);
    setParameter(estimator, "Y1 φ", 0.0);
    setParameter(estimator, "Z2 γ", 1.0e-2);
    setParameter(estimator, "Z2 φ", 1.0);
    estimator.estimate();
    Assert.assertTrue(estimator.getIterationsCount() < 15);
    Assert.assertTrue(estimator.getEvaluationsCount() < 15);
    Assert.assertEquals(0.0, estimator.getOptimum().getRMS(), 1.0e-5);
    Assert.assertEquals(hpaRefX1.getParametersDrivers()[0].getValue(), getParameter(estimator, "X1 γ"), 1.e-12);
    Assert.assertEquals(hpaRefX1.getParametersDrivers()[1].getValue(), getParameter(estimator, "X1 φ"), 1.e-12);
    Assert.assertEquals(hpaRefY1.getParametersDrivers()[0].getValue(), getParameter(estimator, "Y1 γ"), 1.e-12);
    Assert.assertEquals(hpaRefY1.getParametersDrivers()[1].getValue(), getParameter(estimator, "Y1 φ"), 1.e-12);
    Assert.assertEquals(hpaRefZ2.getParametersDrivers()[0].getValue(), getParameter(estimator, "Z2 γ"), 1.e-12);
    Assert.assertEquals(hpaRefZ2.getParametersDrivers()[1].getValue(), getParameter(estimator, "Z2 φ"), 1.e-12);
}

Example 2 with BatchLSEstimator

use of org.orekit.estimation.leastsquares.BatchLSEstimator in project Orekit by CS-SI.

the class GroundStationTest method testEstimateStationPosition.

@Test
public void testEstimateStationPosition() throws OrekitException, IOException, ClassNotFoundException {
    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);
    // 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);
    // move one station
    final RandomGenerator random = new Well19937a(0x4adbecfc743bda60l);
    final TopocentricFrame base = context.stations.get(0).getBaseFrame();
    final BodyShape parent = base.getParentShape();
    final Vector3D baseOrigin = parent.transform(base.getPoint());
    final Vector3D deltaTopo = new Vector3D(2 * random.nextDouble() - 1, 2 * random.nextDouble() - 1, 2 * random.nextDouble() - 1);
    final Transform topoToParent = base.getTransformTo(parent.getBodyFrame(), (AbsoluteDate) null);
    final Vector3D deltaParent = topoToParent.transformVector(deltaTopo);
    final String movedSuffix = "-moved";
    final GroundStation moved = new GroundStation(new TopocentricFrame(parent, parent.transform(baseOrigin.subtract(deltaParent), parent.getBodyFrame(), null), base.getName() + movedSuffix), context.ut1.getEOPHistory(), context.stations.get(0).getDisplacements());
    // create orbit estimator
    final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder);
    for (final ObservedMeasurement<?> measurement : measurements) {
        final Range range = (Range) measurement;
        final String name = range.getStation().getBaseFrame().getName() + movedSuffix;
        if (moved.getBaseFrame().getName().equals(name)) {
            estimator.addMeasurement(new Range(moved, range.getDate(), range.getObservedValue()[0], range.getTheoreticalStandardDeviation()[0], range.getBaseWeight()[0]));
        } else {
            estimator.addMeasurement(range);
        }
    }
    estimator.setParametersConvergenceThreshold(1.0e-3);
    estimator.setMaxIterations(100);
    estimator.setMaxEvaluations(200);
    // we want to estimate station offsets
    moved.getEastOffsetDriver().setSelected(true);
    moved.getNorthOffsetDriver().setSelected(true);
    moved.getZenithOffsetDriver().setSelected(true);
    EstimationTestUtils.checkFit(context, estimator, 2, 3, 0.0, 5.6e-7, 0.0, 1.4e-6, 0.0, 4.8e-7, 0.0, 2.6e-10);
    Assert.assertEquals(deltaTopo.getX(), moved.getEastOffsetDriver().getValue(), 4.5e-7);
    Assert.assertEquals(deltaTopo.getY(), moved.getNorthOffsetDriver().getValue(), 6.2e-7);
    Assert.assertEquals(deltaTopo.getZ(), moved.getZenithOffsetDriver().getValue(), 2.6e-7);
    GeodeticPoint result = moved.getOffsetGeodeticPoint(null);
    GeodeticPoint reference = context.stations.get(0).getBaseFrame().getPoint();
    Assert.assertEquals(reference.getLatitude(), result.getLatitude(), 1.4e-14);
    Assert.assertEquals(reference.getLongitude(), result.getLongitude(), 2.9e-14);
    Assert.assertEquals(reference.getAltitude(), result.getAltitude(), 2.6e-7);
    RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
    RealMatrix physicalCovariances = estimator.getPhysicalCovariances(1.0e-10);
    Assert.assertEquals(9, normalizedCovariances.getRowDimension());
    Assert.assertEquals(9, normalizedCovariances.getColumnDimension());
    Assert.assertEquals(9, physicalCovariances.getRowDimension());
    Assert.assertEquals(9, physicalCovariances.getColumnDimension());
    Assert.assertEquals(0.55431, physicalCovariances.getEntry(6, 6), 1.0e-5);
    Assert.assertEquals(0.22694, physicalCovariances.getEntry(7, 7), 1.0e-5);
    Assert.assertEquals(0.13106, physicalCovariances.getEntry(8, 8), 1.0e-5);
    ByteArrayOutputStream bos = new ByteArrayOutputStream();
    ObjectOutputStream oos = new ObjectOutputStream(bos);
    oos.writeObject(moved.getEstimatedEarthFrame().getTransformProvider());
    Assert.assertTrue(bos.size() > 155000);
    Assert.assertTrue(bos.size() < 160000);
    ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
    ObjectInputStream ois = new ObjectInputStream(bis);
    EstimatedEarthFrameProvider deserialized = (EstimatedEarthFrameProvider) ois.readObject();
    Assert.assertEquals(moved.getPrimeMeridianOffsetDriver().getValue(), deserialized.getPrimeMeridianOffsetDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPrimeMeridianDriftDriver().getValue(), deserialized.getPrimeMeridianDriftDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPolarOffsetXDriver().getValue(), deserialized.getPolarOffsetXDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPolarDriftXDriver().getValue(), deserialized.getPolarDriftXDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPolarOffsetYDriver().getValue(), deserialized.getPolarOffsetYDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPolarDriftYDriver().getValue(), deserialized.getPolarDriftYDriver().getValue(), 1.0e-15);
}

Example 3 with BatchLSEstimator

use of org.orekit.estimation.leastsquares.BatchLSEstimator in project Orekit by CS-SI.

the class OrbitDetermination method createEstimator.

/**
 * Set up estimator.
 * @param parser input file parser
 * @param propagatorBuilder propagator builder
 * @return estimator
 * @throws NoSuchElementException if input parameters are missing
 * @throws OrekitException if some propagator parameters cannot be retrieved
 */
private BatchLSEstimator createEstimator(final KeyValueFileParser<ParameterKey> parser, final NumericalPropagatorBuilder propagatorBuilder) throws NoSuchElementException, OrekitException {
    final boolean optimizerIsLevenbergMarquardt;
    if (!parser.containsKey(ParameterKey.ESTIMATOR_OPTIMIZATION_ENGINE)) {
        optimizerIsLevenbergMarquardt = true;
    } else {
        final String engine = parser.getString(ParameterKey.ESTIMATOR_OPTIMIZATION_ENGINE);
        optimizerIsLevenbergMarquardt = engine.toLowerCase().contains("levenberg");
    }
    final LeastSquaresOptimizer optimizer;
    if (optimizerIsLevenbergMarquardt) {
        // we want to use a Levenberg-Marquardt optimization engine
        final double initialStepBoundFactor;
        if (!parser.containsKey(ParameterKey.ESTIMATOR_LEVENBERG_MARQUARDT_INITIAL_STEP_BOUND_FACTOR)) {
            initialStepBoundFactor = 100.0;
        } else {
            initialStepBoundFactor = parser.getDouble(ParameterKey.ESTIMATOR_LEVENBERG_MARQUARDT_INITIAL_STEP_BOUND_FACTOR);
        }
        optimizer = new LevenbergMarquardtOptimizer().withInitialStepBoundFactor(initialStepBoundFactor);
    } else {
        // we want to use a Gauss-Newton optimization engine
        optimizer = new GaussNewtonOptimizer(Decomposition.QR);
    }
    final double convergenceThreshold;
    if (!parser.containsKey(ParameterKey.ESTIMATOR_NORMALIZED_PARAMETERS_CONVERGENCE_THRESHOLD)) {
        convergenceThreshold = 1.0e-3;
    } else {
        convergenceThreshold = parser.getDouble(ParameterKey.ESTIMATOR_NORMALIZED_PARAMETERS_CONVERGENCE_THRESHOLD);
    }
    final int maxIterations;
    if (!parser.containsKey(ParameterKey.ESTIMATOR_MAX_ITERATIONS)) {
        maxIterations = 10;
    } else {
        maxIterations = parser.getInt(ParameterKey.ESTIMATOR_MAX_ITERATIONS);
    }
    final int maxEvaluations;
    if (!parser.containsKey(ParameterKey.ESTIMATOR_MAX_EVALUATIONS)) {
        maxEvaluations = 20;
    } else {
        maxEvaluations = parser.getInt(ParameterKey.ESTIMATOR_MAX_EVALUATIONS);
    }
    final BatchLSEstimator estimator = new BatchLSEstimator(optimizer, propagatorBuilder);
    estimator.setParametersConvergenceThreshold(convergenceThreshold);
    estimator.setMaxIterations(maxIterations);
    estimator.setMaxEvaluations(maxEvaluations);
    return estimator;
}

Example 4 with BatchLSEstimator

use of org.orekit.estimation.leastsquares.BatchLSEstimator 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();
    }
}

Example 5 with BatchLSEstimator

use of org.orekit.estimation.leastsquares.BatchLSEstimator in project Orekit by CS-SI.

the class GroundStationTest method testEstimateEOP.

@Test
public void testEstimateEOP() throws OrekitException {
    Context linearEOPContext = EstimationTestUtils.eccentricContext("linear-EOP:regular-data/de431-ephemerides:potential:tides");
    final AbsoluteDate refDate = new AbsoluteDate(2000, 2, 24, linearEOPContext.utc);
    final double dut10 = 0.3079738;
    final double lod = 0.0011000;
    final double xp0 = 68450.0e-6;
    final double xpDot = -50.0e-6;
    final double yp0 = 60.0e-6;
    final double ypDot = 2.0e-6;
    for (double dt = -2 * Constants.JULIAN_DAY; dt < 2 * Constants.JULIAN_DAY; dt += 300.0) {
        AbsoluteDate date = refDate.shiftedBy(dt);
        Assert.assertEquals(dut10 - dt * lod / Constants.JULIAN_DAY, linearEOPContext.ut1.getEOPHistory().getUT1MinusUTC(date), 1.0e-15);
        Assert.assertEquals(lod, linearEOPContext.ut1.getEOPHistory().getLOD(date), 1.0e-15);
        Assert.assertEquals((xp0 + xpDot * dt / Constants.JULIAN_DAY) * Constants.ARC_SECONDS_TO_RADIANS, linearEOPContext.ut1.getEOPHistory().getPoleCorrection(date).getXp(), 1.0e-15);
        Assert.assertEquals((yp0 + ypDot * dt / Constants.JULIAN_DAY) * Constants.ARC_SECONDS_TO_RADIANS, linearEOPContext.ut1.getEOPHistory().getPoleCorrection(date).getYp(), 1.0e-15);
    }
    final NumericalPropagatorBuilder linearPropagatorBuilder = linearEOPContext.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
    // create perfect range measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(linearEOPContext.initialOrbit, linearPropagatorBuilder);
    final List<ObservedMeasurement<?>> linearMeasurements = EstimationTestUtils.createMeasurements(propagator, new RangeMeasurementCreator(linearEOPContext), 1.0, 5.0, 60.0);
    Utils.clearFactories();
    Context zeroEOPContext = EstimationTestUtils.eccentricContext("zero-EOP:regular-data/de431-ephemerides:potential:potential:tides");
    for (double dt = -2 * Constants.JULIAN_DAY; dt < 2 * Constants.JULIAN_DAY; dt += 300.0) {
        AbsoluteDate date = refDate.shiftedBy(dt);
        Assert.assertEquals(0.0, zeroEOPContext.ut1.getEOPHistory().getUT1MinusUTC(date), 1.0e-15);
        Assert.assertEquals(0.0, zeroEOPContext.ut1.getEOPHistory().getLOD(date), 1.0e-15);
        Assert.assertEquals(0.0, zeroEOPContext.ut1.getEOPHistory().getPoleCorrection(date).getXp(), 1.0e-15);
        Assert.assertEquals(0.0, zeroEOPContext.ut1.getEOPHistory().getPoleCorrection(date).getYp(), 1.0e-15);
    }
    // create orbit estimator
    final NumericalPropagatorBuilder zeroPropagatorBuilder = linearEOPContext.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
    final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), zeroPropagatorBuilder);
    for (final ObservedMeasurement<?> linearMeasurement : linearMeasurements) {
        Range linearRange = (Range) linearMeasurement;
        for (final GroundStation station : zeroEOPContext.stations) {
            if (station.getBaseFrame().getName().equals(linearRange.getStation().getBaseFrame().getName())) {
                Range zeroRange = new Range(station, linearRange.getDate(), linearRange.getObservedValue()[0], linearRange.getTheoreticalStandardDeviation()[0], linearRange.getBaseWeight()[0]);
                estimator.addMeasurement(zeroRange);
            }
        }
    }
    estimator.setParametersConvergenceThreshold(1.0e-3);
    estimator.setMaxIterations(100);
    estimator.setMaxEvaluations(200);
    // we want to estimate pole and prime meridian
    GroundStation station = zeroEOPContext.stations.get(0);
    station.getPrimeMeridianOffsetDriver().setReferenceDate(refDate);
    station.getPrimeMeridianOffsetDriver().setSelected(true);
    station.getPrimeMeridianDriftDriver().setSelected(true);
    station.getPolarOffsetXDriver().setReferenceDate(refDate);
    station.getPolarOffsetXDriver().setSelected(true);
    station.getPolarDriftXDriver().setSelected(true);
    station.getPolarOffsetYDriver().setReferenceDate(refDate);
    station.getPolarOffsetYDriver().setSelected(true);
    station.getPolarDriftYDriver().setSelected(true);
    // just for the fun and to speed up test, we will use orbit determination, *without* estimating orbit
    for (final ParameterDriver driver : zeroPropagatorBuilder.getOrbitalParametersDrivers().getDrivers()) {
        driver.setSelected(false);
    }
    estimator.estimate();
    final double computedDut1 = station.getPrimeMeridianOffsetDriver().getValue() / EstimatedEarthFrameProvider.EARTH_ANGULAR_VELOCITY;
    final double computedLOD = station.getPrimeMeridianDriftDriver().getValue() * (-Constants.JULIAN_DAY / EstimatedEarthFrameProvider.EARTH_ANGULAR_VELOCITY);
    final double computedXp = station.getPolarOffsetXDriver().getValue() / Constants.ARC_SECONDS_TO_RADIANS;
    final double computedXpDot = station.getPolarDriftXDriver().getValue() / Constants.ARC_SECONDS_TO_RADIANS * Constants.JULIAN_DAY;
    final double computedYp = station.getPolarOffsetYDriver().getValue() / Constants.ARC_SECONDS_TO_RADIANS;
    final double computedYpDot = station.getPolarDriftYDriver().getValue() / Constants.ARC_SECONDS_TO_RADIANS * Constants.JULIAN_DAY;
    Assert.assertEquals(dut10, computedDut1, 4.3e-10);
    Assert.assertEquals(lod, computedLOD, 4.9e-10);
    Assert.assertEquals(xp0, computedXp, 5.6e-9);
    Assert.assertEquals(xpDot, computedXpDot, 7.2e-9);
    Assert.assertEquals(yp0, computedYp, 1.1e-9);
    Assert.assertEquals(ypDot, computedYpDot, 2.8e-11);
// thresholds to use if orbit is estimated
// (i.e. when commenting out the loop above that sets orbital parameters drivers to "not selected")
// Assert.assertEquals(dut10, computedDut1,  6.6e-3);
// Assert.assertEquals(lod,   computedLOD,   1.1e-9);
// Assert.assertEquals(xp0,   computedXp,    3.3e-8);
// Assert.assertEquals(xpDot, computedXpDot, 2.2e-8);
// Assert.assertEquals(yp0,   computedYp,    3.3e-8);
// Assert.assertEquals(ypDot, computedYpDot, 3.8e-8);
}