Struct MvnSpacecraft 

pub struct MvnSpacecraft {
    pub template: Spacecraft,
    pub dispersions: Vec<StateDispersion>,
    pub mean: SVector<f64, 9>,
    pub sqrt_s_v: SMatrix<f64, 9, 9>,
    pub std_norm_distr: Normal<f64>,
}

A multivariate spacecraft state generator for Monte Carlo analyses. Ensures that the covariance is properly applied on all provided state variables.

Algorithm

The MvnSpacecraft allows sampling from a multivariate normal distribution defined by a template state (mean) and a set of state dispersions (uncertainties). The core difficulty is that dispersions are often defined in non-Cartesian spaces (like Keplerian orbital elements or B-Plane parameters), while the spacecraft state is represented in Cartesian coordinates (Position and Velocity).

The algorithm proceeds as follows:

1. Jacobian Computation: It computes the Jacobian matrix J representing the partial derivatives of the provided dispersion parameters with respect to the Cartesian state elements (x, y, z, vx, vy, vz).
2. Covariance Transformation: It constructs a diagonal covariance matrix P in the parameter space (assuming input dispersions are independent in that space). It then transforms this into the Cartesian covariance matrix C using the linear mapping approximation: C = J_inv * P * J_inv^T where J_inv is the Moore-Penrose pseudo-inverse of J.
3. Decomposition: It performs a Singular Value Decomposition (SVD) on C (or the full state covariance including mass, Cr, Cd) to obtain the square root of the covariance matrix, denoted as L. C = U * S * V^T = (V * sqrt(S)) * (V * sqrt(S))^T implies L = V * sqrt(S)
4. Sampling: To generate a sample state X, it draws a vector Z of independent standard normal variables (N(0, 1)) and applies the transformation: X = mu + L * Z

Correctness vs. Independent Sampling

One might ask: “Why not just sample each Cartesian coordinate independently using a normal distribution?”

1. Correlations: Independent sampling of Cartesian coordinates assumes a diagonal covariance matrix, implying no correlation between position and velocity components. In orbital mechanics, states are highly correlated (e.g., velocity magnitude and radial distance are coupled by energy). MvnSpacecraft preserves these physical correlations by mapping the physically meaningful uncertainties (e.g., in SMA or Inclination) into the Cartesian space.
2. Geometry: Uncertainties defined in orbital elements form complex shapes (like “bananas”) in Cartesian space. A multivariate normal approximation in Cartesian space captures the principal axes and orientation of this uncertainty volume, which an axis-aligned bounding box (implied by independent sampling) effectively destroys.
3. Consistency: By using the Jacobian transformation, we ensure that the generated samples, when mapped back to the parameter space (linearized), reproduce the input statistics (mean and standard deviation) provided by the user.

Fields

template: Spacecraft

The template state

dispersions: Vec<StateDispersion>

mean: SVector<f64, 9>

The mean of the multivariate normal distribution

sqrt_s_v: SMatrix<f64, 9, 9>

The dot product \sqrt{\vec s} \cdot \vec v, where S is the singular values and V the V matrix from the SVD decomp of the covariance of multivariate normal distribution

std_norm_distr: Normal<f64>

The standard normal distribution used to seed the multivariate normal distribution

Implementations

impl MvnSpacecraft

pub fn new( template: Spacecraft, dispersions: Vec<StateDispersion>, ) -> Result<Self, Box<dyn Error>>

Creates a new mulivariate state generator from a mean and covariance on the set of state parameters. The covariance must be positive semi definite.

Algorithm

This function will build the rotation matrix to rotate the requested dispersions into the Spacecraft state space using [OrbitDual]. If there are any dispersions on the Cr and Cd, then these are dispersed independently (because they are iid).

Examples found in repository

nyx-core/examples/03_geo_analysis/stationkeeping.rs (lines 102-108)

fn main() -> Result<(), Box<dyn Error>> {
    pel::init();
    // Set up the dynamics like in the orbit raise.
    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);

    // Define the GEO orbit, and we're just going to maintain it very tightly.
    let earth_j2000 = almanac.frame_info(EARTH_J2000)?;
    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
    println!("{orbit:x}");

    let sc = Spacecraft::builder()
        .orbit(orbit)
        .mass(Mass::from_dry_and_prop_masses(1000.0, 1000.0)) // 1000 kg of dry mass and prop, totalling 2.0 tons
        .srp(SRPData::from_area(3.0 * 6.0)) // Assuming 1 kW/m^2 or 18 kW, giving a margin of 4.35 kW for on-propulsion consumption
        .thruster(Thruster {
            // "NEXT-STEP" row in Table 2
            isp_s: 4435.0,
            thrust_N: 0.472,
        })
        .mode(GuidanceMode::Thrust) // Start thrusting immediately.
        .build();

    // Set up the spacecraft dynamics like in the orbit raise example.

    let prop_time = 30.0 * Unit::Day;

    // Define the guidance law -- we're just using a Ruggiero controller as demonstrated in AAS-2004-5089.
    let objectives = &[
        Objective::within_tolerance(
            StateParameter::Element(OrbitalElement::SemiMajorAxis),
            42_165.0,
            20.0,
        ),
        Objective::within_tolerance(
            StateParameter::Element(OrbitalElement::Eccentricity),
            0.001,
            5e-5,
        ),
        Objective::within_tolerance(
            StateParameter::Element(OrbitalElement::Inclination),
            0.05,
            1e-2,
        ),
    ];

    let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2)?;
    println!("{ruggiero_ctrl}");

    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);

    let mut jgm3_meta = MetaFile {
        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
    };
    jgm3_meta.process(true)?;

    let harmonics = GravityField::from_stor(
        almanac.frame_info(IAU_EARTH_FRAME)?,
        GravityFieldData::from_cof(&jgm3_meta.uri, 8, 8, true)?,
    );
    orbital_dyn.accel_models.push(harmonics);

    let srp_dyn = SolarPressure::default_flux(EARTH_J2000, almanac.clone())?;
    let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
        .with_guidance_law(ruggiero_ctrl.clone());

    println!("{sc_dynamics}");

    // Finally, let's use the Monte Carlo framework built into Nyx to propagate spacecraft.

    // Let's start by defining the dispersion.
    // The MultivariateNormal structure allows us to define the dispersions in any of the orbital parameters, but these are applied directly in the Cartesian state space.
    // Note that additional validation on the MVN is in progress -- https://github.com/nyx-space/nyx/issues/339.
    let mc_rv = MvnSpacecraft::new(
        sc,
        vec![StateDispersion::zero_mean(
            StateParameter::Element(OrbitalElement::SemiMajorAxis),
            3.0,
        )],
    )?;

    let my_mc = MonteCarlo::new(
        sc, // Nominal state
        mc_rv,
        "03_geo_sk".to_string(), // Scenario name
        None, // No specific seed specified, so one will be drawn from the computer's entropy.
    );

    // Build the propagator setup.
    let setup = Propagator::rk89(
        sc_dynamics.clone(),
        IntegratorOptions::builder()
            .min_step(10.0_f64.seconds())
            .error_ctrl(ErrorControl::RSSCartesianStep)
            .build(),
    );

    let num_runs = 25;
    let rslts = my_mc.run_until_epoch(setup, almanac.clone(), sc.epoch() + prop_time, num_runs);

    assert_eq!(rslts.runs.len(), num_runs);

    rslts.to_parquet("03_geo_sk.parquet", ExportCfg::default())?;

    Ok(())
}

pub fn zero_mean( template: Spacecraft, dispersions: Vec<StateDispersion>, ) -> Result<Self, Box<dyn Error>>

Same as new but with a zero mean

pub fn from_spacecraft_cov( template: Spacecraft, cov: SMatrix<f64, 9, 9>, mean: SVector<f64, 9>, ) -> Result<Self, Box<dyn Error>>

Initializes a new multivariate distribution using the state data in the spacecraft state space.

Trait Implementations

impl Clone for MvnSpacecraft

fn clone(&self) -> MvnSpacecraft

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Debug for MvnSpacecraft

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl Distribution<DispersedState<Spacecraft>> for MvnSpacecraft

fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> DispersedState<Spacecraft>

Generate a random value of T, using rng as the source of randomness.

fn sample_iter<R>(self, rng: R) -> Iter<Self, R, T>

where R: Rng, Self: Sized,

Create an iterator that generates random values of T, using rng as the source of randomness.

fn map<F, S>(self, func: F) -> Map<Self, F, T, S>

where F: Fn(T) -> S, Self: Sized,

Map sampled values to type S

impl<'a, 'py> FromPyObject<'a, 'py> for MvnSpacecraft

where Self: Clone,

type Error = PyClassGuardError<'a, 'py>

The type returned in the event of a conversion error.

fn extract( obj: Borrowed<'a, 'py, PyAny>, ) -> Result<Self, <Self as FromPyObject<'a, 'py>>::Error>

Extracts Self from the bound smart pointer obj.

impl<'py> IntoPyObject<'py> for MvnSpacecraft

type Target = MvnSpacecraft

The Python output type

type Output = Bound<'py, <MvnSpacecraft as IntoPyObject<'py>>::Target>

The smart pointer type to use.

type Error = PyErr

The type returned in the event of a conversion error.

fn into_pyobject( self, py: Python<'py>, ) -> Result<<Self as IntoPyObject<'_>>::Output, <Self as IntoPyObject<'_>>::Error>

Performs the conversion.

impl PyClass for MvnSpacecraft

const NAME: &str = "MvnSpacecraft"

Name of the class.

type Frozen = False

Whether the pyclass is frozen.

impl PyClassImpl for MvnSpacecraft

const MODULE: Option<&str> = ::core::option::Option::None

Module which the class will be associated with.

const IS_BASETYPE: bool = false

#[pyclass(subclass)]

const IS_SUBCLASS: bool = false

#[pyclass(extends=…)]

const IS_MAPPING: bool = false

#[pyclass(mapping)]

const IS_SEQUENCE: bool = false

#[pyclass(sequence)]

const IS_IMMUTABLE_TYPE: bool = false

#[pyclass(immutable_type)]

const RAW_DOC: &'static CStr = /// A multivariate spacecraft state generator for Monte Carlo analyses. Ensures that the covariance is properly applied on all provided state variables. /// /// # Algorithm /// /// The `MvnSpacecraft` allows sampling from a multivariate normal distribution defined by a template state (mean) and a set of state dispersions (uncertainties). /// The core difficulty is that dispersions are often defined in non-Cartesian spaces (like Keplerian orbital elements or B-Plane parameters), while the spacecraft state is represented in Cartesian coordinates (Position and Velocity). /// /// The algorithm proceeds as follows: /// 1. **Jacobian Computation**: It computes the Jacobian matrix `J` representing the partial derivatives of the provided dispersion parameters with respect to the Cartesian state elements (x, y, z, vx, vy, vz). /// 2. **Covariance Transformation**: It constructs a diagonal covariance matrix `P` in the parameter space (assuming input dispersions are independent in that space). It then transforms this into the Cartesian covariance matrix `C` using the linear mapping approximation: /// `C = J_inv * P * J_inv^T` /// where `J_inv` is the Moore-Penrose pseudo-inverse of `J`. /// 3. **Decomposition**: It performs a Singular Value Decomposition (SVD) on `C` (or the full state covariance including mass, Cr, Cd) to obtain the square root of the covariance matrix, denoted as `L`. /// `C = U * S * V^T = (V * sqrt(S)) * (V * sqrt(S))^T` implies `L = V * sqrt(S)` /// 4. **Sampling**: To generate a sample state `X`, it draws a vector `Z` of independent standard normal variables (N(0, 1)) and applies the transformation: /// `X = mu + L * Z` /// /// # Correctness vs. Independent Sampling /// /// One might ask: "Why not just sample each Cartesian coordinate independently using a normal distribution?" /// /// 1. **Correlations**: Independent sampling of Cartesian coordinates assumes a diagonal covariance matrix, implying no correlation between position and velocity components. In orbital mechanics, states are highly correlated (e.g., velocity magnitude and radial distance are coupled by energy). `MvnSpacecraft` preserves these physical correlations by mapping the physically meaningful uncertainties (e.g., in SMA or Inclination) into the Cartesian space. /// 2. **Geometry**: Uncertainties defined in orbital elements form complex shapes (like "bananas") in Cartesian space. A multivariate normal approximation in Cartesian space captures the principal axes and orientation of this uncertainty volume, which an axis-aligned bounding box (implied by independent sampling) effectively destroys. /// 3. **Consistency**: By using the Jacobian transformation, we ensure that the generated samples, when mapped back to the parameter space (linearized), reproduce the input statistics (mean and standard deviation) provided by the user.

Docstring for the class provided on the struct or enum.

const DOC: &'static CStr

Fully rendered class doc, including the text_signature if a constructor is defined.

type Layout = <<MvnSpacecraft as PyClassImpl>::BaseNativeType as PyClassBaseType>::Layout<MvnSpacecraft>

Description of how this class is laid out in memory

type BaseType = PyAny

Base class

type ThreadChecker = NoopThreadChecker

This handles following two situations:

type Inventory = Pyo3MethodsInventoryForMvnSpacecraft

type PyClassMutability = <<PyAny as PyClassBaseType>::PyClassMutability as PyClassMutability>::MutableChild

Immutable or mutable

type Dict = PyClassDummySlot

Specify this class has #[pyclass(dict)] or not.

type WeakRef = PyClassDummySlot

Specify this class has #[pyclass(weakref)] or not.

type BaseNativeType = PyAny

The closest native ancestor. This is PyAny by default, and when you declare #[pyclass(extends=PyDict)], it’s PyDict.

fn items_iter() -> PyClassItemsIter

fn lazy_type_object() -> &'static LazyTypeObject<Self>

fn dict_offset() -> Option<PyObjectOffset>

Used to provide the dictoffset slot (equivalent to tp_dictoffset)

fn weaklist_offset() -> Option<PyObjectOffset>

Used to provide the weaklistoffset slot (equivalent to tp_weaklistoffset

impl PyClassNewTextSignature for MvnSpacecraft

const TEXT_SIGNATURE: &'static str = "(template, dispersions)"

impl PyTypeInfo for MvnSpacecraft

const NAME: &str = <Self as ::pyo3::PyClass>::NAME

👎Deprecated since 0.28.0:

prefer using ::type_object(py).name() to get the correct runtime value

Class name.

const MODULE: Option<&str> = <Self as ::pyo3::impl_::pyclass::PyClassImpl>::MODULE

👎Deprecated since 0.28.0:

prefer using ::type_object(py).module() to get the correct runtime value

Module name, if any.

fn type_object_raw(py: Python<'_>) -> *mut PyTypeObject

Returns the PyTypeObject instance for this type.

fn type_object(py: Python<'_>) -> Bound<'_, PyType>

Returns the safe abstraction over the type object.

fn is_type_of(object: &Bound<'_, PyAny>) -> bool

Checks if object is an instance of this type or a subclass of this type.

fn is_exact_type_of(object: &Bound<'_, PyAny>) -> bool

Checks if object is an instance of this type.

impl DerefToPyAny for MvnSpacecraft