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

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 = Harmonics::from_stor(
        almanac.frame_info(IAU_EARTH_FRAME)?,
        HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true)?,
    );
    orbital_dyn.accel_models.push(harmonics);

    let srp_dyn = SolarPressure::default(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