State Transition Matrix Across Matching Points¶
Sometimes we have a trajectory which was set up for a mission design or manoeuvre optimisation which uses multiple shooting and matching points.
For these trajectories, we also need to sometimes perform other analyses, e.g. Covariance or Guidance, for which we need the state transition matrix across matching points.
He we show how this is done.
First perform all necessary imports
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from godot.core import tempo, util, constants
from godot.core import autodif as ad
from godot.core.autodif import bridge as br
from godot.model import interface, common, geometry
from godot import cosmos
# avoid verbose output
util.suppressLogger()
next load the universe configuration:
uni_config = cosmos.util.load_yaml('universe.yml')
uni = cosmos.Universe(uni_config)
Now load the trajectory configuration and create the trajectory object using the universe object:
This trajectory has a matching point which splits the trajectory into two parts.
Note that when we create the trajectory we set the
provideStateflag in the constructor toTrue. This is so that the gradients are resolved using trajectory, rather than the propagator.
tra_config = cosmos.util.load_yaml('trajectory.yml')
tra = cosmos.Trajectory(uni, tra_config, True) # <<< True here means, state leaves refer to trajectory, not propagator
tra.compute(True)
The following function is taken from MIDAS, and performs the linear mapping from Cartesian virtual states at the user provided epoch to a model_value. The state_model is the trajectory, the model_value is a `ad.Vector which is mapped to.
from typing import Iterable
def get_linear_mapping(
state_model: interface.State,
control_epoch: tempo.Epoch,
target: interface.VectorTimeEvaluable,
target_epoch: tempo.Epoch,
parameters: Iterable[str],
) -> np.ndarray:
"""Computes the linear mapping matrix from a virtual cartesian state at control epoch to a target model value
Parameters
----------
state_model : interface.State
Mapping state model
control_epoch : tempo.Epoch
Control epoch of the mapping
target : VectorTimeEvaluable
Target model VectorTimeEvaluable
target_epoch : tempo.Epoch
Target epoch of the mapping
parameters : Iterable[str]
All other parameters than state to include in the mapping generation
Returns
-------
np.ndarray
Mapping matrix from control_epoch cartesian state to target_epoch model
Raises
------
RuntimeError
If state mapping model does not contain virtual_state_leaf
"""
# evaluate the target at the target epoch
target_value = target.eval(tempo.XEpoch(target_epoch))
state_dim = np.min(np.shape(state_model.transitionMatrix(control_epoch, control_epoch).state))
# create a virtual gradient identity (cartesian)
virtual_state_leaf = ad.BasicLeaf(state_dim)
virtual_gradient = ad.VectorGradient(virtual_state_leaf, np.identity(state_dim))
# create a resolve target
target = common.ResolveTargetXVector(control_epoch, virtual_gradient)
data = {state_model: target}
# resolve to the target epoch and to the virtual cartesian state parameters
grad = common.resolve(target_value, data).gradient()
# get the matrix, raising an error if the resolve was not done (target not dependent on state at control epoch)
if grad.contains(virtual_state_leaf):
mat = grad[virtual_state_leaf]
else:
raise RuntimeError("State mapping model does not have access to required gradients")
# add the parameter parts
for p in parameters:
if p not in grad:
mat = np.column_stack((mat, np.zeros((len(target_value), 1))))
else:
mat = np.column_stack((mat, grad[p]))
return mat
To use it, we first define and evaluate the target TimeEvaluable at the target epoch:
target_epoch = tempo.Epoch('10150 TDB')
target = geometry.Vector6(uni.frames, "Earth", "SC_center", "ICRF")
Calling the linear mapping function:
Note the states vector are not used, other than for their dimension. If present, they simply create virtual states with the same dimension and identity partials.
control_epoch = tempo.Epoch('9750 TDB')
map_control_to_target = get_linear_mapping(tra, control_epoch, target, target_epoch, [])
print(map_control_to_target)
Now, we repeat, this time constructing the trajectory with the False flag. This demonstrates the error which will be raised
uni_config = cosmos.util.load_yaml('universe.yml')
uni = cosmos.Universe(uni_config)
tra_config = cosmos.util.load_yaml('trajectory.yml')
tra2 = cosmos.Trajectory(uni, tra_config, False) # <<< False here means, state leaves refer to propagator, not trajectory
tra2.compute(True)
control_epoch2 = tempo.Epoch('9750 TDB')
target_epoch2 = tempo.Epoch('10150 TDB')
target2 = geometry.Vector6(uni.frames, "Earth", "SC_center", "ICRF")
# This will result in a failure because the propagator will not be able to
try:
map_control_to_target = get_linear_mapping(tra2, control_epoch2, target2, target_epoch2, [])
except Exception as e:
print(f"Exception raised: {e}")