Simulation¶

In this section, the different stages of a typical simulation setup are described. The user is guided along the creation of the environment, the dynamical model, output variables, termination settings, the integrator, and the actual simulation itself. Code examples will illustrate the usage of the application programming interface (API), both in Python and C++. (C++ is WIP)

Environment Setup¶

In Tudat, the physical environment is defined by a system of bodies, each encapsulated in a Body object. Such an object may represent a celestial body, or a manmade vehicle, and Tudat makes no a priori distinction between the two. The distinction is made by the user when creating the bodies The combination of all Body objects is stored in a bodies object.

The typical procedure to create the environment is the following:

Load default settings for celestial bodies in the simulation
Modify default settings as required for the simulation (if needed)
Create celestial bodies based on these settings
Create any additional bodies which have no defaults
Assign properties to these additional bodies

All the options available to you for creating bodies are given in the following pages:

Details on all available environment models, and how to define them in your simulation, is given below:

Specifics on the environment during propagations can be found here:

Propagation Setup¶

In Tudat, the term ‘propagator’ is used to denote the formulation of the differential equations that are to be numerically solved. This includes the type of dynamics (translational, rotational, mass), but also the formulation for a given type of dynamics (e.g. Encke, Cowell, Kepler elements for translational dynamics).

Propagator Setup¶

Tudat currently supports the propagation of three types of dynamics: translational, rotational and mass. Any combination of any number of types of dynamics, for any number of bodies, may be defined.

On the following page, we give a brief description of how propagate the dynamics of multiple bodies concurrently:

Multi-body Dynamics

The following page provides you with the difference between conventional and propagated coordinates in the propagator settings. It is important to keep in mind that this page covers some details which happen ‘under the hood’.

Conventional vs. Propagated Coordinates

Acceleration Model Setup¶

To propagate translational dynamics, you must provide a set of acceleration models. The acceleration model setup is provided here:

Acceleration Model Setup

Below, a comprehensive list of all available acceleration models in Tudat, and the manner in which to define them, is given

Available Acceleration Models

Torque Model Setup¶

Output Variables¶

By default, propagating the dynamics of a body provides only the numerically integrated state history of your model as output. Tudat has the option to provide any number of additional outputs from your simulation, by defining the output_variables input to the propagator settings. A comprehensive list of available outputs is given below

Available Dependent Variables

Termination Settings¶

The typical propagation is terminated when a specific final time is reached. Tudat provides additional possibilities for defining when the simulator terminates the propagation. A full list is given below:

Available Termination Settings

Integrator Setup¶

The environment and formulation of dynamical equations are now in place. In order to solve these equations you still need to define the numerical integrator settings. These settings specify how the equations are solved. In Tudat(Py) you can choose between different types of integrators.

Since the choice of integrator strongly depends on the nature of the dynamical problem and the requirements of the user, there is no ‘best’ integrator that works in all cases. Details about the different types will not be given here; you are referred to existing literature on the topic of numerical integrators. We only show you how to create the integrator settings, given on the following page:

Integrator Settings

Running the Simulation¶

With all the necessary simulation settings in place, it is time to run the simulation. In Tudat(Py), this is done by means of a DynamicsSimulator object, which handles the setup and execution of the simulation. It also contains functions to retrieve the propagated state history and dependent variables for further analysis and plotting.

In its simplest form, the DynamicsSimulator is used as shown in this example:

Required Show/Hide

# Create simulation object and propagate dynamics.
dynamics_simulator = propagation_setup.SingleArcDynamicsSimulator(
        bodies, integrator_settings, propagator_settings, True)

states = dynamics_simulator.state_history

First, a SingleArcDynamicsSimulator is created using the system of bodies, integrator settings, and propagator settings objects. Tudat will then automatically read and setup the simulation accordingly. The True at the end of the line indicates that the equations of motion should be integrated immediately after creating the object, such that the state history can be retrieved afterwards.

The latter is done in the next line. The simulator will return a dictionary (Python) or map (C++) containing the state of the vehicle at each epoch, which can be exported or used for subsequent analysis.

If the user chose to also export dependent variables, they can be extracted from the dynamics simulator as follows:

Required Show/Hide

# Retrieve dependent variables
dependent_variable_history = dynamics_simulator.get_dependent_variable_history()

This concludes the section on simulation of this API guide. For more detailed information, refer to the pages listed in this section or refer to the next few sections for information and examples on e.g. interpolators, coordinate and time conversions, and interfaces to Spice, JSON, and Sofa.

Linear Sensitivity Analysis¶

Up to this point, we have been concerned with propagating states of bodies only. Tudat is also capable of propagating the so-called variational equations associated with the dynamics to produce the state transition matrix \(\Phi(t,t_{0})\) and sensitivity matrix \(S(t)\), which we define here as:

\[\begin{split}\Phi(t,t_{0}) &= \frac{\partial \mathbf{x}(t)}{\partial\mathbf{x}(t_{0})}\\ S &= \frac{\partial \mathbf{x}(t)}{\partial \mathbf{p }}\\\end{split}\]

where \(\mathbf{x}\) is the propagated state, \(\mathbf{p}\) the vector of a parameter vector (e.g. gravity field parameters, rotation model parameters, etc.), and \(t_{0}\) denotes the initial time. These two matrices are based on linearization of the complex dynamics and can be used to quickly determine the influence of a change in initial state (\(\mathbf{x}(t_{0})\)) and/or parameters (\(\mathbf{p}\)) on the state \(\mathbf{x}(t)\) at time \(t\).

Warning

Note that linearization comes with a caveat. These matrices only produce accurate results within a ‘reasonable’ time frame around \(t_{0}\) and with a ‘small’ change in parameters/initial state. The words between quotation marks indicate that the user should always determine whether the linearization is valid for the case at hand; no rules can be given as the answer strongly depends on the kind of dynamics and problem parameters.

Note

In some literature, the sensitivity matrix is not defined separately, but the state transition matrix \(\Phi(t,t_{0})\) is defined as \(\frac{\partial[\mathbf{x}(t);\text{ }\mathbf{p}]}{\partial[\mathbf{x}(t_{0};\text{ }\mathbf{p}])}\)

The propagation of these equations is done similarly as for the dynamics, and is executed by a (derived class of) VariationalEquationsSolver.

At the moment, the following VariationalEquationsSolver options are available or under development in Tudat:

Single-arc (TudatPy + Tudat);
Multi-arc (Tudat only);
Hybrid (under development).

These are implemented in derived classes and are discussed below. Note that these objects of these classes propagate both the variational equations and dynamics (either concurrently or sequentially).

The parameter estimation framework of Tudat allows an ever increasing variety of parameters to be estimated, these parameters may be:

Properties of a body, such as a gravitational parameter \(\mu\)
Properties of a ground station, such as its body-fixed position \(\mathbf{x}_{GS}^{(B)}\)
Global properties of the simulation, such a Parameterize Post_Newtonian (PPN) parameters \(\gamma\) and \(\beta\)
Acceleration model properties, such as empirical acceleration magnitudes
Observation model properties, such as absolute and relative observation biases

In Tudat, these parameters influence the simulation in a variety of manners, and during propagation and/or observation simulation, information of this parameter is transferred in manner different ways. To provide a unified framework for estimating any type of parameter, the EstimatableParameter class has been set up.

class EstimatableParameter¶: This class has interfaces to retrieve and reset parameters, providing a single interface for modifying/obtaining any of the parameters that Tudat supports. For each estimated parameter, there is a dedicated derived class of EstimatableParameter.

class EstimatableParameterSet¶: The full list of estimated parameters is stored in an object of type EstimatableParameterSet. This class is templated by the state scalar type of the estimated initial state parameters.

Note

For the remainder of this page, we will implicitly assume that the template argument of an EstimatableParameterSet object is double, unless explicitly mentioned otherwise.

As is the case for acceleration models, integration models, environment models, etc., the parameter objects are created by calling factory functions in the estimation_setup.parameter.* namespace. Most parameters need one or more arguments which contain settings for setting up the desired parameter. See Available Estimated Parameters for more details on the available parameters in TudatPy. Most factory functions return a single parameter, such that the list can be constructed using the [parameter1, parameter_2] syntax in Python. There are, however, a few that return a list of parameters, such as those that create initial state parameters. Adding them to the list using the comma operator would result in a multi-dimensional list, which will raise an error when creating the variational equations simulator object. To prevent this, use the following syntax:

parameter_settings = [
                      estimation_setup.parameter.gravitational_parameter("Earth"),
                      estimation_setup.parameter.constant_drag_coefficient("Delfi-C3"),
                      estimation_setup.parameter.radiation_pressure_coefficient("Delfi-C3")
                     ]
                     + estimation_setup.parameter.initial_states(propagator_settings, bodies)

In the snippet above, parameters are created to estimate the gravitational parameter of the Earth, the constant drag and the radiation pressure coefficient of Delfi-C3, and the initial states of the propagated bodies.

Now that the parameter settings have been created, the variational equations solver can be set up as well. This solver acts as a replacement of the normal dynamics simulator:

variational_equations_solver = estimation_setup.SingleArcVariationalEquationsSolver(
        bodies, integrator_settings, propagator_settings,
        estimation_setup.create_parameters_to_estimate(parameter_settings, bodies),
        integrate_on_creation=1 )

The state history, state transition matrices, and sensitivity matrices can then be extracted:

states = variational_equations_solver.state_history
state_transition_matrices = variational_equations_solver.state_transition_matrix_history
sensitivity_matrices = variational_equations_solver.sensitivity_matrix_history

For a complete example of propagating the variational equations, please see the tutorial Propagating Variational Equations.

The EstimatableParameterSet object contains three objects that have EstimatableParameter as base class (one for each parameter). We distinguish two types of EstimatableParameter objects:

Those that represent initial conditions for dynamics (denoted as \(\mathbf{x}_{0}\) below)
Those that represent fixed parameters for environment, acceleration or observation models (denoted as \(\mathbf{q}\) below)

Resetting the full parameter vector \(\mathbf{p}(=[\mathbf{x}_{0};\mathbf{q}])\) is done as follows (for double state scalar type):

// Create parameter set
std::shared_ptr< EstimatableParameterSet< double > > parametersToEstimate = ...

Eigen::VectorXd parameterVector =
     parametersToEstimate->getFullParameterValues< double >( );

While resetting the full parameter vector is done as:

// Create parameter set
std::shared_ptr< EstimatableParameterSet< double > > parametersToEstimate = ...

// Define vector of new values of estimated parameters
Eigen::VectorXd newParameterVector = ...

// Reset parameter values
parametersToEstimate->resetParameterValues< double >( );

When resetting the parameter vector, the change in the values in \(\mathbf{q}\) immediately take effect. For the initial state parameters to take effect, however, the dynamics must be re-propagated. This occurs automatically when estimating parameters. It can also be performed manually by calling the resetParameterEstimate member function of the VariationalEquationsSolver class.

Available Estimated Parameters