Note

Generated by nbsphinx from a Jupyter notebook. All the examples as Jupyter notebooks are available in the tudatpy-examples repo.

# Covariance analysis of estimated parameters

Copyright (c) 2010-2022, Delft University of Technology. All rights reserved. This file is part of the Tudat. Redistribution and use in source and binary forms, with or without modification, are permitted exclusively under the terms of the Modified BSD license. You should have received a copy of the license with this file. If not, please or visit: http://tudat.tudelft.nl/LICENSE.

## Context

This example implements the following aspects of orbit estimation: - Simulation of a spacecraft orbit. - Modelling of a tracking station on Earth. - Simulation of Doppler data at 1 mm/s every 60 seconds, during periods where the spacecraft is at an elevation angle above 15deg, as viewed from the station. - Use of the simulated data to estimate the spacecraft initial state, drag coefficient and radiation pressure coefficient.

## Import statements

The required import statements are made here, at the very beginning.

Some standard modules are first loaded: numpy and matplotlib.pyplot.

Then, the different modules of tudatpy that will be used are imported. Most notably, the estimation, estimation_setup, and observations modules are new to this example.

:

# Load standard modules
import numpy as np
# Uncomment the following to make plots interactive
# %matplotlib widget
from matplotlib import pyplot as plt

from tudatpy.kernel import constants
from tudatpy.kernel.interface import spice
from tudatpy.kernel import numerical_simulation
from tudatpy.kernel.numerical_simulation import environment_setup
from tudatpy.kernel.numerical_simulation import propagation_setup
from tudatpy.kernel.numerical_simulation import estimation, estimation_setup
from tudatpy.kernel.numerical_simulation.estimation_setup import observation
from tudatpy.kernel.astro import element_conversion


## Configuration

NAIF’s SPICE kernels are first loaded, so that the position of various bodies such as the Earth, the Sun, or the Moon, can be make known to tudatpy.

Then, the start and end simulation epochs are setups. In this case, the start epoch is set to 0, corresponding to the 1st of January 2000. The end epoch is set 3 days after. The times should be specified in seconds since J2000. Please refer to the API documentation of the time_conversion module here for more information on this.

:

# Load spice kernels

# Set simulation start and end epochs
simulation_start_epoch = 0.0
simulation_end_epoch = 3 * constants.JULIAN_DAY


## Environment setup

Let’s create the environment for our simulation. This setup covers the creation of (celestial) bodies, vehicle(s), and environment interfaces.

### Create the main bodies

Bodies can be created by making a list of strings with the bodies that is to be included in the simulation.

The default body settings (such as atmosphere, body shape, rotation model) are taken from SPICE.

These settings can be adjusted. Please refere to the Available Environment Models in the user guide for more details.

Finally, the system of bodies is created using the settings. This system of bodies is stored into the variable bodies.

:

# Create default body settings for "Sun", "Earth", "Moon", "Mars", and "Venus"
bodies_to_create = ["Sun", "Earth", "Moon", "Mars", "Venus"]

# Create default body settings for bodies_to_create, with "Earth"/"J2000" as the global frame origin and orientation
global_frame_origin = "Earth"
global_frame_orientation = "J2000"
body_settings = environment_setup.get_default_body_settings(
bodies_to_create, global_frame_origin, global_frame_orientation)

# Create system of bodies
bodies = environment_setup.create_system_of_bodies(body_settings)


### Create the vehicle and its environment interface

Let’s now create the satellite for which an orbit will be simulated.

This satellite is setup to have mass of 400kg, a reference area (used both for aerodynamic and radiation pressure) of 4m$$^2$$, a radiation pressure coefficient of 1.2, and a drag coefficient also of 1.2.

When setting up the radiation pressure interface, the Earth is set as a body that can occult the radiation emitted by the Sun.

:

# Create vehicle objects.
bodies.create_empty_body("Delfi-C3")
bodies.get("Delfi-C3").mass = 400.0

# Create aerodynamic coefficient interface settings
reference_area = 4.0
drag_coefficient = 1.2
aero_coefficient_settings = environment_setup.aerodynamic_coefficients.constant(
reference_area, [drag_coefficient, 0.0, 0.0]
)
# Add the aerodynamic interface to the environment

occulting_bodies = ["Earth"]
)


## Propagation setup

Now that the environment is created, the propagation itself can be setup.

First, the bodies to be propagated and the central bodies will be defined. Central bodies are the bodies with respect to which the state of the respective propagated bodies is defined.

:

# Define bodies that are propagated
bodies_to_propagate = ["Delfi-C3"]

# Define central bodies of propagation
central_bodies = ["Earth"]


### Create the acceleration model

The acceleration settings that act on Delfi-C3 are now to be defined. In this case, these consist in the followings: * Gravitational acceleration using a Point Mass model from: - The Sun - The Moon - Mars * Gravitational acceleration using a Spherical Harmonic approximation up to degree and order 8 from Earth. * Aerodynamic acceleration from Earth. * Acceleration caused by the radiation pressure of the Sun on the vehicle approximated as a cannonball.

The acceleration settings defined are then applied to Delfi-C3 in a dictionary.

This dictionary is finally input to the propagation setup to create the acceleration models.

:

# Define the accelerations acting on Delfi-C3
accelerations_settings_delfi_c3 = dict(
Sun=[
propagation_setup.acceleration.point_mass_gravity()
],
Mars=[
propagation_setup.acceleration.point_mass_gravity()
],
Moon=[
propagation_setup.acceleration.point_mass_gravity()
],
Earth=[
propagation_setup.acceleration.spherical_harmonic_gravity(8, 8),
propagation_setup.acceleration.aerodynamic()
])

# Create global accelerations dictionary
acceleration_settings = {"Delfi-C3": accelerations_settings_delfi_c3}

# Create acceleration models
acceleration_models = propagation_setup.create_acceleration_models(
bodies,
acceleration_settings,
bodies_to_propagate,
central_bodies)


### Define the initial state

The initial state of the vehicle that will be propagated is now defined.

This initial state always has to be provided as a cartesian state, in the form of a list with the first three elements reprensenting the initial position, and the three remaining elements representing the initial velocity.

In this case, let’s make use of the keplerian_to_cartesian_elementwise() function that is included in the element_conversion module, so that the initial state can be input as Keplerian elements, and then converted to Cartesian elements.

:

# Set the initial state of the vehicle
earth_gravitational_parameter = bodies.get("Earth").gravitational_parameter
initial_state = element_conversion.keplerian_to_cartesian_elementwise(
gravitational_parameter=earth_gravitational_parameter,
semi_major_axis=7500.0E3,
eccentricity=0.1,
)


### Create the propagator settings

The propagator is finally setup.

First, a termination condition is defined so that the propagation will stop when the end epochs that was defined is reached.

Then, the translational propagation settings are defined. These are used to simulate the orbit of Delfi-C3 around Earth.

:

# Create termination settings
termination_condition = propagation_setup.propagator.time_termination(simulation_end_epoch)

# Create propagation settings
propagator_settings = propagation_setup.propagator.translational(
central_bodies,
acceleration_models,
bodies_to_propagate,
initial_state,
termination_condition
)


### Create the integrator settings

The last step before starting the simulation is to setup the integrator that will be used.

In this case, a RKF78 integrator is used with a step fixed at 60 seconds (this is done by setting both the minimum and maximum step size to 60 seconds, and the tolerance to 1).

:

# Create numerical integrator settings
integrator_settings = propagation_setup.integrator.runge_kutta_variable_step_size(
simulation_start_epoch, 60.0, propagation_setup.integrator.rkf_78, 60.0, 60.0, 1.0, 1.0
)


## Observation Setup

Now that the dynamical model of the simulated orbit is set, the observational model needs to be defined.

In order to simulate observations, which are required for an estimation of the system parameters, the simulation environment has to be extended by a ground station model. In this case, the ground station is modelled to be located on Earth, in Delft, at an altitude of 0m, 52.00667°N, 4.35556°E.

More information on how to use the add_ground_station() function can be found on tudatpy API website.

:

# Define the position of the ground station on Earth
station_altitude = 0.0

# Add the ground station to the environment
bodies.get_body("Earth"),
"TrackingStation",
[station_altitude, delft_latitude, delft_longitude],
element_conversion.geodetic_position_type)


### Define Observation Simulation Settings

In the following cell the settings for the simulation of the observations from the previously defined observation models are defined. For each observation model, the observation simulation settings set the times at which observations are simulated and defines the viability criteria and noise of the observation.

:

# Define observation simulation times for each link (separated by steps of 1 minute)
observation_times = np.arange(simulation_start_epoch, simulation_end_epoch, 60.0)
observation_simulation_settings = observation.tabulated_simulation_settings(
observation.one_way_doppler_type,
observation_times
)

# Add noise levels of roughly 3.3E-12 [s/m] and add this as Gaussian noise to the observation
noise_level = 1.0E-3 / constants.SPEED_OF_LIGHT
[observation_simulation_settings],
noise_level,
observation.one_way_doppler_type
)

# Create viability settings
[observation_simulation_settings],
[viability_setting]
)


The actual simulation of the observations requires Observation Simulators, which are created automatically by the Estimator object. Therefore, the observations cannot be simulated before the creation of an Estimator object.

## Estimation Setup

After the setup of the observations, let’s now setup the estimation.

### Defining the estimatable parameters

The estimation of system parameters requires the equations of the dynamical and observational models to be differentiated w.r.t. the estimatable system parameters, resulting in the so-called Variational Equations. More on variational equations can be found in the tudatpy user guide.

In the following cell, the set of estimatable parameters is defined. For this examplary estimation, the initial state of the Delfi s/c, its drag coefficient and the gravitational parameter of the Earth are included as estimatable parameters. A list of the available estimatable parameters is given in tudatpy user guide.

During the creation of the Estimator object, the parameter set will be used to automatically set up the Variational Equations of the simulated orbit and the observations thereof. The propagation of the Variational Equations does not have to be triggered manually, but is also handled automatically by the Estimator object.

:

# Setup parameters settings to propagate the state transition matrix
parameter_settings = estimation_setup.parameter.initial_states(propagator_settings, bodies)

# Add estimated parameters to the sensitivity matrix that will be propagated
parameter_settings.append(estimation_setup.parameter.gravitational_parameter("Earth"))
parameter_settings.append(estimation_setup.parameter.constant_drag_coefficient("Delfi-C3"))

# Create the parameters that will be estimated
parameters_to_estimate = estimation_setup.create_parameter_set(parameter_settings, bodies)


### Creating the Estimator object

The Estimator object consolidates all the information that is required for the estimation of system parameters: * the environment (bodies) * the parameter set (parameters_to_estimate) * observation models (observation_settings_list) * integrator setup (integrator_settings) * dynamical and numerical setup (propagator_settings)

The Estimator object automatically creates the relevant Variational Equations and Observation Simulators, which are required for the parameter estimation and observation simulation, respectively.

:

# Create the estimation object
estimator = numerical_simulation.Estimator(
bodies,
parameters_to_estimate,
observation_settings_list,
integrator_settings,
propagator_settings)


### Perform the observations simulation

Now that the Estimator object exists, we can make use of its observation simulator members. The observation simulations are performed by calling the simulation_observations() function from the observation module. Note that this function also makes use of the previously defined observation simulation settings!

:

# Simulate required observation
simulated_observations = estimation.simulate_observations(
[observation_simulation_settings],
estimator.observation_simulators,
bodies)


## Estimate the parameters

With the simulated observations and the Estimator object, containing the Variational Equations for the estimatable parameters, we now have everything we need in order to conduct the estimation of the parameters.

### Setup the inversion

Here we collect all inputs required for the inversion in a Pod (Precise orbit determination) input object and define some basic settings of the inversion.

:

# Save the true parameters to later analyse the error
truth_parameters = parameters_to_estimate.parameter_vector

# Create input object for estimation, adding observations and parameter set information
pod_input = estimation.PodInput(
simulated_observations, parameters_to_estimate.parameter_set_size)

# set methodological options
pod_input.define_estimation_settings(
reintegrate_variational_equations=False)

# define weighting of the observations in the inversion
weights_per_observable = \
{estimation_setup.observation.one_way_doppler_type: noise_level ** -2}
pod_input.set_constant_weight_per_observable(weights_per_observable)


### Perform the estimation

Finally, let’s run the estimation itself, using the inputs that have been defined earlier.

:

# Perform estimation (this also prints the residuals and partials)
pod_output = estimator.perform_estimation(pod_input)

Calculating residuals and partials 238
Parameter update   -0.200661    -0.179491     -0.32771  0.000201621  0.000221888 -0.000136231 -4.08822e+07 -7.08494e-05
Current residual: 3.22598e-12
Calculating residuals and partials 238
Parameter update  7.2185e-07  1.41211e-06  5.01767e-06 -4.44237e-09 -1.47038e-09   2.0183e-10      293.169  2.40208e-09
Current residual: 3.146e-12
Calculating residuals and partials 238
Parameter update-1.72519e-07 -5.55437e-07 -4.04358e-06  3.53916e-09  1.64167e-09  1.96357e-10     -161.685 -4.42769e-09
Current residual: 3.146e-12
Calculating residuals and partials 238
Parameter update-8.15745e-07  1.14748e-08  6.88105e-07 -5.11148e-10 -1.84499e-10 -2.87299e-10      -51.288  4.31811e-09
Current residual: 3.146e-12
Calculating residuals and partials 238
Parameter update 6.19134e-07  4.95267e-07  1.15434e-06 -1.08884e-09  -5.4084e-10  2.44444e-10       117.12  7.86453e-10
Current residual: 3.146e-12
Maximum number of iterations reached
Final residual: 3.146e-12

:

# Print the estimation error
print(pod_output.formal_errors)
print(truth_parameters - parameters_to_estimate.parameter_vector)

[1.60858847e-01 1.70140339e-01 1.67143774e-01 1.34696394e-04
1.24426709e-04 1.41853340e-04 2.20984801e+07 3.74865311e-05]
[ 2.00660857e-01  1.79489939e-01  3.27707272e-01 -2.01618353e-04
-2.21887790e-04  1.36231084e-04  4.08819638e+07  7.08463522e-05]


## Results post-processing

Let’s now process the results of the estimation and observations.

### Correlation

First, let’s plot the correlation between the outputs of the estimation.

:

plt.figure(figsize=(9,5))
plt.imshow(np.abs(pod_output.correlations), aspect='auto', interpolation='none')
plt.colorbar()
plt.tight_layout()
plt.show() ### Range rate over time

Let’s now make a plot of the observations range rate over time.

:

observation_times = np.array(simulated_observations.concatenated_times)
observations_list = np.array(simulated_observations.concatenated_observations)

plt.figure(figsize=(9, 5))
plt.title("Observations as a function of time.")
plt.scatter(observation_times / 3600.0, observations_list * constants.SPEED_OF_LIGHT)
plt.xlabel("Time [hr]")
plt.ylabel("Range rate [m/s]")
plt.grid()
plt.tight_layout()
plt.show() ### Residuals history

The history of residuals from the estimation can also been analysed by plotting their history.

:

residual_history = pod_output.residual_history
fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(9, 11))
subplots_list = [ax1, ax2, ax3, ax4, ax5]

for i in range(5):
subplots_list[i].scatter(observation_times, residual_history[:, i])
plt.tight_layout()
plt.show() ### Final residuals

Finally, we can analyse the final residuals, print them, and plot them.

:

print(pod_output.formal_errors / (truth_parameters - parameters_to_estimate.parameter_vector))
final_residuals = pod_output.final_residuals

plt.figure(figsize=(9,5))
plt.hist(final_residuals, 25)
plt.tight_layout()
plt.show()

[ 0.80164537  0.94791017  0.51003987 -0.66807606 -0.56076411  1.04126999
0.54054351  0.52912437] 