Examples with RealMatrix org.hipparchus.linear.RealMatrix used on opensource projects

Example 26 with RealMatrix

use of org.hipparchus.linear.RealMatrix 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);
}

Example 27 with RealMatrix

use of org.hipparchus.linear.RealMatrix in project Orekit by CS-SI.

the class KalmanOrbitDeterminationTest method displayFinalCovariances.

/**
 * Display covariances and sigmas as predicted by a Kalman filter at date t.
 * @throws OrekitException
 */
private void displayFinalCovariances(final PrintStream logStream, final KalmanEstimator kalman) throws OrekitException {
    // // Get kalman estimated propagator
    // final NumericalPropagator kalmanProp = kalman.getProcessModel().getEstimatedPropagator();
    // 
    // // Link the partial derivatives to this propagator
    // final String equationName = "kalman-derivatives";
    // PartialDerivativesEquations kalmanDerivatives = new PartialDerivativesEquations(equationName, kalmanProp);
    // 
    // // Initialize the derivatives
    // final SpacecraftState rawState = kalmanProp.getInitialState();
    // final SpacecraftState stateWithDerivatives =
    // kalmanDerivatives.setInitialJacobians(rawState);
    // kalmanProp.resetInitialState(stateWithDerivatives);
    // 
    // // Propagate to target date
    // final SpacecraftState kalmanState = kalmanProp.propagate(targetDate);
    // 
    // // Compute STM
    // RealMatrix STM = kalman.getProcessModel().getErrorStateTransitionMatrix(kalmanState, kalmanDerivatives);
    // 
    // // Compute covariance matrix
    // RealMatrix P = kalman.getProcessModel().unNormalizeCovarianceMatrix(kalman.predictCovariance(STM,
    // kalman.getProcessModel().getProcessNoiseMatrix()));
    final RealMatrix P = kalman.getPhysicalEstimatedCovarianceMatrix();
    final String[] paramNames = new String[P.getRowDimension()];
    int index = 0;
    int paramSize = 0;
    for (final ParameterDriver driver : kalman.getOrbitalParametersDrivers(true).getDrivers()) {
        paramNames[index++] = driver.getName();
        paramSize = FastMath.max(paramSize, driver.getName().length());
    }
    for (final ParameterDriver driver : kalman.getPropagationParametersDrivers(true).getDrivers()) {
        paramNames[index++] = driver.getName();
        paramSize = FastMath.max(paramSize, driver.getName().length());
    }
    for (final ParameterDriver driver : kalman.getEstimatedMeasurementsParameters().getDrivers()) {
        paramNames[index++] = driver.getName();
        paramSize = FastMath.max(paramSize, driver.getName().length());
    }
    if (paramSize < 20) {
        paramSize = 20;
    }
    // Header
    logStream.format("\n%s\n", "Kalman Final Covariances:");
    // logStream.format(Locale.US, "\tDate: %-23s UTC\n",
    // targetDate.toString(TimeScalesFactory.getUTC()));
    logStream.format(Locale.US, "\tDate: %-23s UTC\n", kalman.getCurrentDate().toString(TimeScalesFactory.getUTC()));
    // Covariances
    String strFormat = String.format("%%%2ds  ", paramSize);
    logStream.format(strFormat, "Covariances:");
    for (int i = 0; i < P.getRowDimension(); i++) {
        logStream.format(Locale.US, strFormat, paramNames[i]);
    }
    logStream.println("");
    String numFormat = String.format("%%%2d.6f  ", paramSize);
    for (int i = 0; i < P.getRowDimension(); i++) {
        logStream.format(Locale.US, strFormat, paramNames[i]);
        for (int j = 0; j <= i; j++) {
            logStream.format(Locale.US, numFormat, P.getEntry(i, j));
        }
        logStream.println("");
    }
    // Correlation coeff
    final double[] sigmas = new double[P.getRowDimension()];
    for (int i = 0; i < P.getRowDimension(); i++) {
        sigmas[i] = FastMath.sqrt(P.getEntry(i, i));
    }
    logStream.format("\n" + strFormat, "Corr coef:");
    for (int i = 0; i < P.getRowDimension(); i++) {
        logStream.format(Locale.US, strFormat, paramNames[i]);
    }
    logStream.println("");
    for (int i = 0; i < P.getRowDimension(); i++) {
        logStream.format(Locale.US, strFormat, paramNames[i]);
        for (int j = 0; j <= i; j++) {
            logStream.format(Locale.US, numFormat, P.getEntry(i, j) / (sigmas[i] * sigmas[j]));
        }
        logStream.println("");
    }
    // Sigmas
    logStream.format("\n" + strFormat + "\n", "Sigmas: ");
    for (int i = 0; i < P.getRowDimension(); i++) {
        logStream.format(Locale.US, strFormat + numFormat + "\n", paramNames[i], sigmas[i]);
    }
    logStream.println("");
}

Example 28 with RealMatrix

use of org.hipparchus.linear.RealMatrix in project Orekit by CS-SI.

the class KalmanOrbitDeterminationTest method testW3B.

@Test
public // Orbit determination for range, azimuth elevation measurements
void testW3B() throws URISyntaxException, IllegalArgumentException, IOException, OrekitException, ParseException {
    // Print results on console
    final boolean print = false;
    // Input in tutorial resources directory/output
    final String inputPath = KalmanOrbitDeterminationTest.class.getClassLoader().getResource("orbit-determination/W3B/od_test_W3.in").toURI().getPath();
    final File input = new File(inputPath);
    // Configure Orekit data access
    Utils.setDataRoot("orbit-determination/W3B:potential/icgem-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
    // Choice of an orbit type to use
    // Default for test is Cartesian
    final OrbitType orbitType = OrbitType.CARTESIAN;
    // Initial orbital Cartesian covariance matrix
    // These covariances are derived from the deltas between initial and reference orbits
    // So in a way they are "perfect"...
    // Cartesian covariance matrix initialization
    final RealMatrix cartesianOrbitalP = MatrixUtils.createRealDiagonalMatrix(new double[] { FastMath.pow(2.4e4, 2), FastMath.pow(1.e5, 2), FastMath.pow(4.e4, 2), FastMath.pow(3.5, 2), FastMath.pow(2., 2), FastMath.pow(0.6, 2) });
    // Orbital Cartesian process noise matrix (Q)
    final RealMatrix cartesianOrbitalQ = MatrixUtils.createRealDiagonalMatrix(new double[] { 1.e-4, 1.e-4, 1.e-4, 1.e-10, 1.e-10, 1.e-10 });
    // Propagation covariance and process noise matrices
    final RealMatrix propagationP = MatrixUtils.createRealDiagonalMatrix(new double[] { // Cd
    FastMath.pow(2., 2), // leak-X
    FastMath.pow(5.7e-6, 2), // leak-X
    FastMath.pow(1.1e-11, 2), // leak-Y
    FastMath.pow(7.68e-7, 2), // leak-Y
    FastMath.pow(1.26e-10, 2), // leak-Z
    FastMath.pow(5.56e-6, 2), // leak-Z
    FastMath.pow(2.79e-10, 2) });
    final RealMatrix propagationQ = MatrixUtils.createRealDiagonalMatrix(new double[] { // Cd
    FastMath.pow(1e-3, 2), // Leaks
    0., // Leaks
    0., // Leaks
    0., // Leaks
    0., // Leaks
    0., // Leaks
    0. });
    // Measurement covariance and process noise matrices
    // az/el bias sigma = 0.06deg
    // range bias sigma = 100m
    final double angularVariance = FastMath.pow(FastMath.toRadians(0.06), 2);
    final double rangeVariance = FastMath.pow(500., 2);
    final RealMatrix measurementP = MatrixUtils.createRealDiagonalMatrix(new double[] { angularVariance, angularVariance, rangeVariance, angularVariance, angularVariance, rangeVariance, angularVariance, angularVariance, rangeVariance, angularVariance, angularVariance, rangeVariance, angularVariance, angularVariance, rangeVariance });
    // Process noise sigma: 1e-6 for all
    final double measQ = FastMath.pow(1e-6, 2);
    final RealMatrix measurementQ = MatrixUtils.createRealIdentityMatrix(measurementP.getRowDimension()).scalarMultiply(measQ);
    // Kalman orbit determination run.
    ResultKalman kalmanW3B = run(input, orbitType, print, cartesianOrbitalP, cartesianOrbitalQ, propagationP, propagationQ, measurementP, measurementQ);
    // Tests
    // -----
    // Definition of the accuracy for the test
    final double distanceAccuracy = 0.1;
    // degrees
    final double angleAccuracy = 1e-5;
    // Number of measurements processed
    final int numberOfMeas = 521;
    Assert.assertEquals(numberOfMeas, kalmanW3B.getNumberOfMeasurements());
    // Test on propagator parameters
    // -----------------------------
    // Batch LS result
    // final double dragCoef  = -0.2154;
    final double dragCoef = 0.1931;
    Assert.assertEquals(dragCoef, kalmanW3B.propagatorParameters.getDrivers().get(0).getValue(), 1e-3);
    final Vector3D leakAcceleration0 = new Vector3D(kalmanW3B.propagatorParameters.getDrivers().get(1).getValue(), kalmanW3B.propagatorParameters.getDrivers().get(3).getValue(), kalmanW3B.propagatorParameters.getDrivers().get(5).getValue());
    // Batch LS results
    // Assert.assertEquals(8.002e-6, leakAcceleration0.getNorm(), 1.0e-8);
    Assert.assertEquals(5.994e-6, leakAcceleration0.getNorm(), 1.0e-8);
    final Vector3D leakAcceleration1 = new Vector3D(kalmanW3B.propagatorParameters.getDrivers().get(2).getValue(), kalmanW3B.propagatorParameters.getDrivers().get(4).getValue(), kalmanW3B.propagatorParameters.getDrivers().get(6).getValue());
    // Batch LS results
    // Assert.assertEquals(3.058e-10, leakAcceleration1.getNorm(), 1.0e-12);
    Assert.assertEquals(1.831e-10, leakAcceleration1.getNorm(), 1.0e-12);
    // Test on measurements parameters
    // -------------------------------
    final List<DelegatingDriver> list = new ArrayList<DelegatingDriver>();
    list.addAll(kalmanW3B.measurementsParameters.getDrivers());
    sortParametersChanges(list);
    // Station CastleRock
    // Batch LS results
    // final double[] CastleAzElBias  = { 0.062701342, -0.003613508 };
    // final double   CastleRangeBias = 11274.4677;
    final double[] CastleAzElBias = { 0.062635, -0.003672 };
    final double CastleRangeBias = 11289.3678;
    Assert.assertEquals(CastleAzElBias[0], FastMath.toDegrees(list.get(0).getValue()), angleAccuracy);
    Assert.assertEquals(CastleAzElBias[1], FastMath.toDegrees(list.get(1).getValue()), angleAccuracy);
    Assert.assertEquals(CastleRangeBias, list.get(2).getValue(), distanceAccuracy);
    // Station Fucino
    // Batch LS results
    // final double[] FucAzElBias  = { -0.053526137, 0.075483886 };
    // final double   FucRangeBias = 13467.8256;
    final double[] FucAzElBias = { -0.053298, 0.075589 };
    final double FucRangeBias = 13482.0715;
    Assert.assertEquals(FucAzElBias[0], FastMath.toDegrees(list.get(3).getValue()), angleAccuracy);
    Assert.assertEquals(FucAzElBias[1], FastMath.toDegrees(list.get(4).getValue()), angleAccuracy);
    Assert.assertEquals(FucRangeBias, list.get(5).getValue(), distanceAccuracy);
    // Station Kumsan
    // Batch LS results
    // final double[] KumAzElBias  = { -0.023574208, -0.054520756 };
    // final double   KumRangeBias = 13512.57594;
    final double[] KumAzElBias = { -0.022805, -0.055057 };
    final double KumRangeBias = 13502.7459;
    Assert.assertEquals(KumAzElBias[0], FastMath.toDegrees(list.get(6).getValue()), angleAccuracy);
    Assert.assertEquals(KumAzElBias[1], FastMath.toDegrees(list.get(7).getValue()), angleAccuracy);
    Assert.assertEquals(KumRangeBias, list.get(8).getValue(), distanceAccuracy);
    // Station Pretoria
    // Batch LS results
    // final double[] PreAzElBias = { 0.030201539, 0.009747877 };
    // final double PreRangeBias = 13594.11889;
    final double[] PreAzElBias = { 0.030353, 0.009658 };
    final double PreRangeBias = 13609.2516;
    Assert.assertEquals(PreAzElBias[0], FastMath.toDegrees(list.get(9).getValue()), angleAccuracy);
    Assert.assertEquals(PreAzElBias[1], FastMath.toDegrees(list.get(10).getValue()), angleAccuracy);
    Assert.assertEquals(PreRangeBias, list.get(11).getValue(), distanceAccuracy);
    // Station Uralla
    // Batch LS results
    // final double[] UraAzElBias = { 0.167814449, -0.12305252 };
    // final double UraRangeBias = 13450.26738;
    final double[] UraAzElBias = { 0.167519, -0.122842 };
    final double UraRangeBias = 13441.7019;
    Assert.assertEquals(UraAzElBias[0], FastMath.toDegrees(list.get(12).getValue()), angleAccuracy);
    Assert.assertEquals(UraAzElBias[1], FastMath.toDegrees(list.get(13).getValue()), angleAccuracy);
    Assert.assertEquals(UraRangeBias, list.get(14).getValue(), distanceAccuracy);
    // Test on statistic for the range residuals
    final long nbRange = 182;
    // statistics for the range residual (min, max, mean, std)
    final double[] RefStatRange = { -12.981, 18.046, -1.133, 5.312 };
    Assert.assertEquals(nbRange, kalmanW3B.getRangeStat().getN());
    Assert.assertEquals(RefStatRange[0], kalmanW3B.getRangeStat().getMin(), distanceAccuracy);
    Assert.assertEquals(RefStatRange[1], kalmanW3B.getRangeStat().getMax(), distanceAccuracy);
    Assert.assertEquals(RefStatRange[2], kalmanW3B.getRangeStat().getMean(), distanceAccuracy);
    Assert.assertEquals(RefStatRange[3], kalmanW3B.getRangeStat().getStandardDeviation(), distanceAccuracy);
    // test on statistic for the azimuth residuals
    final long nbAzi = 339;
    // statistics for the azimuth residual (min, max, mean, std)
    final double[] RefStatAzi = { -0.041441, 0.023473, -0.004426, 0.009911 };
    Assert.assertEquals(nbAzi, kalmanW3B.getAzimStat().getN());
    Assert.assertEquals(RefStatAzi[0], kalmanW3B.getAzimStat().getMin(), angleAccuracy);
    Assert.assertEquals(RefStatAzi[1], kalmanW3B.getAzimStat().getMax(), angleAccuracy);
    Assert.assertEquals(RefStatAzi[2], kalmanW3B.getAzimStat().getMean(), angleAccuracy);
    Assert.assertEquals(RefStatAzi[3], kalmanW3B.getAzimStat().getStandardDeviation(), angleAccuracy);
    // test on statistic for the elevation residuals
    final long nbEle = 339;
    final double[] RefStatEle = { -0.025399, 0.043345, 0.001011, 0.010636 };
    Assert.assertEquals(nbEle, kalmanW3B.getElevStat().getN());
    Assert.assertEquals(RefStatEle[0], kalmanW3B.getElevStat().getMin(), angleAccuracy);
    Assert.assertEquals(RefStatEle[1], kalmanW3B.getElevStat().getMax(), angleAccuracy);
    Assert.assertEquals(RefStatEle[2], kalmanW3B.getElevStat().getMean(), angleAccuracy);
    Assert.assertEquals(RefStatEle[3], kalmanW3B.getElevStat().getStandardDeviation(), angleAccuracy);
    RealMatrix covariances = kalmanW3B.getCovariances();
    Assert.assertEquals(28, covariances.getRowDimension());
    Assert.assertEquals(28, covariances.getColumnDimension());
    // drag coefficient variance
    Assert.assertEquals(0.016349, covariances.getEntry(6, 6), 1.0e-5);
    // leak-X constant term variance
    Assert.assertEquals(2.047303E-13, covariances.getEntry(7, 7), 1.0e-16);
    // leak-Y constant term variance
    Assert.assertEquals(5.462497E-13, covariances.getEntry(9, 9), 1.0e-15);
    // leak-Z constant term variance
    Assert.assertEquals(1.717781E-11, covariances.getEntry(11, 11), 1.0e-15);
    // Test on orbital parameters
    // Done at the end to avoid changing the estimated propagation parameters
    // ----------------------------------------------------------------------
    // Estimated position and velocity
    final Vector3D estimatedPos = kalmanW3B.getEstimatedPV().getPosition();
    final Vector3D estimatedVel = kalmanW3B.getEstimatedPV().getVelocity();
    // Reference position and velocity at initial date (same as in batch LS test)
    final Vector3D refPos0 = new Vector3D(-40541446.255, -9905357.41, 206777.413);
    final Vector3D refVel0 = new Vector3D(759.0685, -1476.5156, 54.793);
    // Gather the selected propagation parameters and initialize them to the values found
    // with the batch LS method
    final ParameterDriversList refPropagationParameters = kalmanW3B.propagatorParameters;
    final double dragCoefRef = -0.215433133145843;
    final double[] leakXRef = { +5.69040439901955E-06, 1.09710906802403E-11 };
    final double[] leakYRef = { -7.66440256777678E-07, 1.25467464335066E-10 };
    final double[] leakZRef = { -5.574055079952E-06, 2.78703463746911E-10 };
    for (DelegatingDriver driver : refPropagationParameters.getDrivers()) {
        switch(driver.getName()) {
            case "drag coefficient":
                driver.setValue(dragCoefRef);
                break;
            case "leak-X[0]":
                driver.setValue(leakXRef[0]);
                break;
            case "leak-X[1]":
                driver.setValue(leakXRef[1]);
                break;
            case "leak-Y[0]":
                driver.setValue(leakYRef[0]);
                break;
            case "leak-Y[1]":
                driver.setValue(leakYRef[1]);
                break;
            case "leak-Z[0]":
                driver.setValue(leakZRef[0]);
                break;
            case "leak-Z[1]":
                driver.setValue(leakZRef[1]);
                break;
        }
    }
    // Run the reference until Kalman last date
    final Orbit refOrbit = runReference(input, orbitType, refPos0, refVel0, refPropagationParameters, kalmanW3B.getEstimatedPV().getDate());
    // Test on last orbit
    final Vector3D refPos = refOrbit.getPVCoordinates().getPosition();
    final Vector3D refVel = refOrbit.getPVCoordinates().getVelocity();
    // Check distances
    final double dP = Vector3D.distance(refPos, estimatedPos);
    final double dV = Vector3D.distance(refVel, estimatedVel);
    // FIXME: debug - Comparison with batch LS is bad
    final double debugDistanceAccuracy = 234.73;
    final double debugVelocityAccuracy = 0.086;
    Assert.assertEquals(0.0, Vector3D.distance(refPos, estimatedPos), debugDistanceAccuracy);
    Assert.assertEquals(0.0, Vector3D.distance(refVel, estimatedVel), debugVelocityAccuracy);
    // Print orbit deltas
    if (print) {
        System.out.println("Test performances:");
        System.out.format("\t%-30s\n", "ΔEstimated / Reference");
        System.out.format(Locale.US, "\t%-10s %20.6f\n", "ΔP [m]", dP);
        System.out.format(Locale.US, "\t%-10s %20.6f\n", "ΔV [m/s]", dV);
    }
}

Example 29 with RealMatrix

use of org.hipparchus.linear.RealMatrix in project Orekit by CS-SI.

the class KalmanOrbitDeterminationTest method run.

/**
 * Function running the Kalman filter estimation.
 * @param input Input configuration file
 * @param orbitType Orbit type to use (calculation and display)
 * @param print Choose whether the results are printed on console or not
 * @param cartesianOrbitalP Orbital part of the initial covariance matrix in Cartesian formalism
 * @param cartesianOrbitalQ Orbital part of the process noise matrix in Cartesian formalism
 * @param propagationP Propagation part of the initial covariance matrix
 * @param propagationQ Propagation part of the process noise matrix
 * @param measurementP Measurement part of the initial covariance matrix
 * @param measurementQ Measurement part of the process noise matrix
 */
private ResultKalman run(final File input, final OrbitType orbitType, final boolean print, final RealMatrix cartesianOrbitalP, final RealMatrix cartesianOrbitalQ, final RealMatrix propagationP, final RealMatrix propagationQ, final RealMatrix measurementP, final RealMatrix measurementQ) throws IOException, IllegalArgumentException, OrekitException, ParseException {
    // Read input parameters
    KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class);
    parser.parseInput(input.getAbsolutePath(), new FileInputStream(input));
    // Log files
    final RangeLog rangeLog = new RangeLog();
    final RangeRateLog rangeRateLog = new RangeRateLog();
    final AzimuthLog azimuthLog = new AzimuthLog();
    final ElevationLog elevationLog = new ElevationLog();
    final PositionLog positionLog = new PositionLog();
    final VelocityLog velocityLog = new VelocityLog();
    // Gravity field
    GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-5c.gfc", true));
    final NormalizedSphericalHarmonicsProvider gravityField = createGravityField(parser);
    // Orbit initial guess
    Orbit initialGuess = createOrbit(parser, gravityField.getMu());
    // Convert to desired orbit type
    initialGuess = orbitType.convertType(initialGuess);
    // 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);
    // 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)));
    }
    // Building the Kalman filter:
    // - Gather the estimated measurement parameters in a list
    // - Prepare the initial covariance matrix and the process noise matrix
    // - Build the Kalman filter
    // --------------------------------------------------------------------
    // Build the list of estimated measurements
    final ParameterDriversList estimatedMeasurementsParameters = new ParameterDriversList();
    for (ObservedMeasurement<?> measurement : measurements) {
        final List<ParameterDriver> drivers = measurement.getParametersDrivers();
        for (ParameterDriver driver : drivers) {
            if (driver.isSelected()) {
                // Add the driver
                estimatedMeasurementsParameters.add(driver);
            }
        }
    }
    // Sort the list lexicographically
    estimatedMeasurementsParameters.sort();
    // Orbital covariance matrix initialization
    // Jacobian of the orbital parameters w/r to Cartesian
    final double[][] dYdC = new double[6][6];
    initialGuess.getJacobianWrtCartesian(propagatorBuilder.getPositionAngle(), dYdC);
    final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
    RealMatrix orbitalP = Jac.multiply(cartesianOrbitalP.multiply(Jac.transpose()));
    // Orbital process noise matrix
    RealMatrix orbitalQ = Jac.multiply(cartesianOrbitalQ.multiply(Jac.transpose()));
    // Build the full covariance matrix and process noise matrix
    final int nbPropag = (propagationP != null) ? propagationP.getRowDimension() : 0;
    final int nbMeas = (measurementP != null) ? measurementP.getRowDimension() : 0;
    final RealMatrix initialP = MatrixUtils.createRealMatrix(6 + nbPropag + nbMeas, 6 + nbPropag + nbMeas);
    final RealMatrix Q = MatrixUtils.createRealMatrix(6 + nbPropag + nbMeas, 6 + nbPropag + nbMeas);
    // Orbital part
    initialP.setSubMatrix(orbitalP.getData(), 0, 0);
    Q.setSubMatrix(orbitalQ.getData(), 0, 0);
    // Propagation part
    if (propagationP != null) {
        initialP.setSubMatrix(propagationP.getData(), 6, 6);
        Q.setSubMatrix(propagationQ.getData(), 6, 6);
    }
    // Measurement part
    if (measurementP != null) {
        initialP.setSubMatrix(measurementP.getData(), 6 + nbPropag, 6 + nbPropag);
        Q.setSubMatrix(measurementQ.getData(), 6 + nbPropag, 6 + nbPropag);
    }
    // Build the Kalman
    KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
    kalmanBuilder.builder(propagatorBuilder);
    kalmanBuilder.estimatedMeasurementsParameters(estimatedMeasurementsParameters);
    kalmanBuilder.initialCovarianceMatrix(initialP);
    kalmanBuilder.processNoiseMatrixProvider(new ConstantProcessNoise(Q));
    final KalmanEstimator kalman = kalmanBuilder.build();
    // Add an observer
    kalman.setObserver(new KalmanObserver() {

        /**
         * Date of the first measurement.
         */
        private AbsoluteDate t0;

        /**
         * {@inheritDoc}
         * @throws OrekitException
         */
        @Override
        @SuppressWarnings("unchecked")
        public void evaluationPerformed(final KalmanEstimation estimation) throws OrekitException {
            // Current measurement number, date and status
            final EstimatedMeasurement<?> estimatedMeasurement = estimation.getCorrectedMeasurement();
            final int currentNumber = estimation.getCurrentMeasurementNumber();
            final AbsoluteDate currentDate = estimatedMeasurement.getDate();
            final EstimatedMeasurement.Status currentStatus = estimatedMeasurement.getStatus();
            // Current estimated measurement
            final ObservedMeasurement<?> observedMeasurement = estimatedMeasurement.getObservedMeasurement();
            // Measurement type & Station name
            String measType = "";
            String stationName = "";
            // Register the measurement in the proper measurement logger
            if (observedMeasurement instanceof Range) {
                // Add the tuple (estimation, prediction) to the log
                rangeLog.add(currentNumber, (EstimatedMeasurement<Range>) estimatedMeasurement);
                // Measurement type & Station name
                measType = "RANGE";
                stationName = ((EstimatedMeasurement<Range>) estimatedMeasurement).getObservedMeasurement().getStation().getBaseFrame().getName();
            } else if (observedMeasurement instanceof RangeRate) {
                rangeRateLog.add(currentNumber, (EstimatedMeasurement<RangeRate>) estimatedMeasurement);
                measType = "RANGE_RATE";
                stationName = ((EstimatedMeasurement<RangeRate>) estimatedMeasurement).getObservedMeasurement().getStation().getBaseFrame().getName();
            } else if (observedMeasurement instanceof AngularAzEl) {
                azimuthLog.add(currentNumber, (EstimatedMeasurement<AngularAzEl>) estimatedMeasurement);
                elevationLog.add(currentNumber, (EstimatedMeasurement<AngularAzEl>) estimatedMeasurement);
                measType = "AZ_EL";
                stationName = ((EstimatedMeasurement<AngularAzEl>) estimatedMeasurement).getObservedMeasurement().getStation().getBaseFrame().getName();
            } else if (observedMeasurement instanceof PV) {
                positionLog.add(currentNumber, (EstimatedMeasurement<PV>) estimatedMeasurement);
                velocityLog.add(currentNumber, (EstimatedMeasurement<PV>) estimatedMeasurement);
                measType = "PV";
            }
            // Header
            if (print) {
                if (currentNumber == 1) {
                    // Set t0 to first measurement date
                    t0 = currentDate;
                    // Print header
                    final String formatHeader = "%-4s\t%-25s\t%15s\t%-10s\t%-10s\t%-20s\t%20s\t%20s";
                    String header = String.format(Locale.US, formatHeader, "Nb", "Epoch", "Dt[s]", "Status", "Type", "Station", "DP Corr", "DV Corr");
                    // Orbital drivers
                    for (DelegatingDriver driver : estimation.getEstimatedOrbitalParameters().getDrivers()) {
                        header += String.format(Locale.US, "\t%20s", driver.getName());
                        header += String.format(Locale.US, "\t%20s", "D" + driver.getName());
                    }
                    // Propagation drivers
                    for (DelegatingDriver driver : estimation.getEstimatedPropagationParameters().getDrivers()) {
                        header += String.format(Locale.US, "\t%20s", driver.getName());
                        header += String.format(Locale.US, "\t%20s", "D" + driver.getName());
                    }
                    // Measurements drivers
                    for (DelegatingDriver driver : estimation.getEstimatedMeasurementsParameters().getDrivers()) {
                        header += String.format(Locale.US, "\t%20s", driver.getName());
                        header += String.format(Locale.US, "\t%20s", "D" + driver.getName());
                    }
                    // Print header
                    System.out.println(header);
                }
                // Print current measurement info in terminal
                String line = "";
                // Line format
                final String lineFormat = "%4d\t%-25s\t%15.3f\t%-10s\t%-10s\t%-20s\t%20.9e\t%20.9e";
                // Orbital correction = DP & DV between predicted orbit and estimated orbit
                final Vector3D predictedP = estimation.getPredictedSpacecraftStates()[0].getPVCoordinates().getPosition();
                final Vector3D predictedV = estimation.getPredictedSpacecraftStates()[0].getPVCoordinates().getVelocity();
                final Vector3D estimatedP = estimation.getCorrectedSpacecraftStates()[0].getPVCoordinates().getPosition();
                final Vector3D estimatedV = estimation.getCorrectedSpacecraftStates()[0].getPVCoordinates().getVelocity();
                final double DPcorr = Vector3D.distance(predictedP, estimatedP);
                final double DVcorr = Vector3D.distance(predictedV, estimatedV);
                line = String.format(Locale.US, lineFormat, currentNumber, currentDate.toString(), currentDate.durationFrom(t0), currentStatus.toString(), measType, stationName, DPcorr, DVcorr);
                // Handle parameters printing (value and error)
                int jPar = 0;
                final RealMatrix Pest = estimation.getPhysicalEstimatedCovarianceMatrix();
                // Orbital drivers
                for (DelegatingDriver driver : estimation.getEstimatedOrbitalParameters().getDrivers()) {
                    line += String.format(Locale.US, "\t%20.9f", driver.getValue());
                    line += String.format(Locale.US, "\t%20.9e", FastMath.sqrt(Pest.getEntry(jPar, jPar)));
                    jPar++;
                }
                // Propagation drivers
                for (DelegatingDriver driver : estimation.getEstimatedPropagationParameters().getDrivers()) {
                    line += String.format(Locale.US, "\t%20.9f", driver.getValue());
                    line += String.format(Locale.US, "\t%20.9e", FastMath.sqrt(Pest.getEntry(jPar, jPar)));
                    jPar++;
                }
                // Measurements drivers
                for (DelegatingDriver driver : estimatedMeasurementsParameters.getDrivers()) {
                    line += String.format(Locale.US, "\t%20.9f", driver.getValue());
                    line += String.format(Locale.US, "\t%20.9e", FastMath.sqrt(Pest.getEntry(jPar, jPar)));
                    jPar++;
                }
                // Print the line
                System.out.println(line);
            }
        }
    });
    // Process the list measurements
    final Orbit estimated = kalman.processMeasurements(measurements).getInitialState().getOrbit();
    // Get the last estimated physical covariances
    final RealMatrix covarianceMatrix = kalman.getPhysicalEstimatedCovarianceMatrix();
    // Parameters and measurements.
    final ParameterDriversList propagationParameters = kalman.getPropagationParametersDrivers(true);
    final ParameterDriversList measurementsParameters = kalman.getEstimatedMeasurementsParameters();
    // Eventually, print parameter changes, statistics and covariances
    if (print) {
        // Display parameter change for non orbital drivers
        int length = 0;
        for (final ParameterDriver parameterDriver : propagationParameters.getDrivers()) {
            length = FastMath.max(length, parameterDriver.getName().length());
        }
        for (final ParameterDriver parameterDriver : measurementsParameters.getDrivers()) {
            length = FastMath.max(length, parameterDriver.getName().length());
        }
        if (propagationParameters.getNbParams() > 0) {
            displayParametersChanges(System.out, "Estimated propagator parameters changes: ", true, length, propagationParameters);
        }
        if (measurementsParameters.getNbParams() > 0) {
            displayParametersChanges(System.out, "Estimated measurements parameters changes: ", true, length, measurementsParameters);
        }
        // Measurements statistics summary
        System.out.println("");
        rangeLog.displaySummary(System.out);
        rangeRateLog.displaySummary(System.out);
        azimuthLog.displaySummary(System.out);
        elevationLog.displaySummary(System.out);
        positionLog.displaySummary(System.out);
        velocityLog.displaySummary(System.out);
        // Covariances and sigmas
        displayFinalCovariances(System.out, kalman);
    }
    // Instantiation of the results
    return new ResultKalman(propagationParameters, measurementsParameters, kalman.getCurrentMeasurementNumber(), estimated.getPVCoordinates(), rangeLog.createStatisticsSummary(), rangeRateLog.createStatisticsSummary(), azimuthLog.createStatisticsSummary(), elevationLog.createStatisticsSummary(), positionLog.createStatisticsSummary(), velocityLog.createStatisticsSummary(), covarianceMatrix);
}

Example 30 with RealMatrix

use of org.hipparchus.linear.RealMatrix 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());
}