Monte Carlo dispersion¶
A 1000-run Monte Carlo around a nominal injection burn. The four perturbed axes are a launch-time-slip proxy (CoastTime.Value — extra time on the parking orbit before the burn fires) and the three components of the injection delta-V in the spacecraft's VNB local frame (Inj.Element1/2/3). Each run propagates from the parking orbit, applies the perturbed burn, propagates a fixed 1200 s, and reports its final inertial position.
The aggregated DataFrame carries one row per (run_id, time) — final-step rows give the per-run arrival point, and the empirical covariance of those arrival points yields the 3-sigma ellipse in the (X, Y) plane.
Prerequisites. A local GMAT install (R2026a is the primary development target; see Supported versions) and pip install gmat-sweep[examples] for the matplotlib dependency.
Set up the run¶
Resolve the GMAT install once and confirm the script that ships next to this notebook is where we expect it. The script declares a parking orbit, an Inj impulsive burn in VNB, and a CoastTime variable that controls the time spent on the parking orbit before the burn fires.
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
from gmat_run import locate_gmat
from gmat_sweep import expand_monte_carlo_to_run_specs, monte_carlo
install = locate_gmat()
script_path = Path("injection_dispersion.script").resolve()
print(f"GMAT version: {install.version}")
print(f"Script: {script_path.name}")
print(f"Exists: {script_path.exists()}")
GMAT version: R2026a Script: injection_dispersion.script Exists: True
Define the perturbations¶
Four independent perturbed axes:
CoastTime.Value— normal, mean 600 s, 1-sigma 30 s. Models a launch-time slip as extra time on the parking orbit before the burn. Operationally a 30-second slip translates to a 30-second extra coast in the same orbital geometry.Inj.Element1/2/3— normal, mean (1.0, 0.0, 0.0) km/s, 1-sigma 5 m/s on each axis. The (V, N, B) components of the injection delta-V; the nominal is a 1 km/s burn along the velocity vector and per-axis dispersions are 0.5 % of the along-track magnitude.
The seed is recorded on the manifest the sweep writes; two calls at the same (n, perturb, seed) produce bit-equal draws, regardless of process or machine.
SEED = 20260505
N = 1000
perturb = {
"CoastTime.Value": ("normal", 600.0, 30.0), # +/- 30 s 1-sigma launch-time slip
"Inj.Element1": ("normal", 1.0, 0.005), # +/- 5 m/s along-track 1-sigma
"Inj.Element2": ("normal", 0.0, 0.005), # +/- 5 m/s normal 1-sigma
"Inj.Element3": ("normal", 0.0, 0.005), # +/- 5 m/s binormal 1-sigma
}
Run the Monte Carlo¶
One monte_carlo() call dispatches all 1000 runs through the default LocalJoblibPool, drains the per-run outcomes in completion order, and returns the aggregated (run_id, time)-MultiIndexed DataFrame. With out=None the per-run Parquet files land in a temporary directory whose lifetime is tied to the returned DataFrame.
Failed runs would surface as one NaN-filled row each with __status="failed". The perturbations here are small enough relative to the parking orbit's energy that every run completes cleanly.
df = monte_carlo(
script_path,
n=N,
perturb=perturb,
seed=SEED,
progress=False,
)
df["__status"].value_counts()
__status ok 11997 Name: count, dtype: int64
Final-step arrival points¶
Group by run_id and take the last row per group. The (Sat.X, Sat.Y, Sat.Z) of the final time-step is the arrival point in the EarthMJ2000Eq inertial frame.
final = df.groupby(level="run_id").last()
arrivals = final[["Sat.X", "Sat.Y", "Sat.Z"]].to_numpy()
print(f"Arrival points shape: {arrivals.shape}")
print(f"Mean (km): {arrivals.mean(axis=0).round(3)}")
print(f"Stddev (km): {arrivals.std(axis=0).round(3)}")
Arrival points shape: (1000, 3) Mean (km): [-3015.051 6826.621 3706.294] Stddev (km): [249.009 85.78 46.776]
Histogram of miss distances¶
The miss distance for each run is the Euclidean distance from its arrival point to the empirical mean arrival point. The histogram of those distances is the canonical 1-D summary of an N-dimensional dispersion.
miss_km = np.linalg.norm(arrivals - arrivals.mean(axis=0), axis=1)
fig, ax = plt.subplots(figsize=(8, 4))
ax.hist(miss_km, bins=40, color="steelblue", edgecolor="white")
ax.axvline(np.percentile(miss_km, 50), color="black", linestyle="--", label="median")
ax.axvline(
np.percentile(miss_km, 99.7),
color="firebrick",
linestyle="--",
label="99.7 percentile",
)
ax.set_xlabel("Miss distance from mean arrival (km)")
ax.set_ylabel("Number of runs")
ax.set_title(f"Arrival miss distance over {len(miss_km)} runs")
ax.legend()
fig.tight_layout()
plt.show()
3-sigma ellipse in the (X, Y) plane¶
Project the 3-D arrival cloud onto the (X, Y) plane and overlay the 3-sigma ellipse derived from the empirical 2x2 covariance. The ellipse semi-axes are 3-sigma_max and 3-sigma_min, where (sigma_max, sigma_min) are the square roots of the covariance matrix's eigenvalues; the rotation angle is the orientation of the major eigenvector.
xy = arrivals[:, :2]
xy_mean = xy.mean(axis=0)
cov = np.cov(xy, rowvar=False)
eigvals, eigvecs = np.linalg.eigh(cov)
order = np.argsort(eigvals)[::-1]
eigvals = eigvals[order]
eigvecs = eigvecs[:, order]
sigma_major, sigma_minor = np.sqrt(eigvals)
theta = np.arctan2(eigvecs[1, 0], eigvecs[0, 0])
phi = np.linspace(0.0, 2 * np.pi, 200)
ellipse_local = 3.0 * np.array([sigma_major * np.cos(phi), sigma_minor * np.sin(phi)])
rot = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
ellipse = rot @ ellipse_local + xy_mean[:, None]
fig, ax = plt.subplots(figsize=(7, 6))
ax.scatter(xy[:, 0], xy[:, 1], s=4, alpha=0.4, color="steelblue", label="arrivals")
ax.scatter([xy_mean[0]], [xy_mean[1]], color="black", marker="x", label="mean")
ax.plot(ellipse[0], ellipse[1], color="firebrick", lw=1.5, label="3-sigma ellipse")
ax.set_xlabel("Sat.X (km)")
ax.set_ylabel("Sat.Y (km)")
ax.set_title("Arrival cloud in the (X, Y) plane")
ax.legend(loc="upper right")
ax.set_aspect("equal", adjustable="datalim")
ax.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()
Determinism contract¶
monte_carlo() is deterministic in (mission, n, perturb, seed). Two calls at the same seed produce bit-equal draws at every run_id; switching to a different seed produces a different draw set.
expand_monte_carlo_to_run_specs() exposes the per-run RunSpec.overrides dict the workers would receive — useful for inspecting a draw set without paying the cost of running it. Two calls at the same (perturb, n, seed) build identical specs.
specs_a = expand_monte_carlo_to_run_specs(
perturb, n=N, seed=SEED, script_path=script_path, output_dir=Path(".")
)
specs_b = expand_monte_carlo_to_run_specs(
perturb, n=N, seed=SEED, script_path=script_path, output_dir=Path(".")
)
specs_c = expand_monte_carlo_to_run_specs(
perturb, n=N, seed=SEED + 1, script_path=script_path, output_dir=Path(".")
)
same_seed_identical = all(
a.overrides == b.overrides for a, b in zip(specs_a, specs_b, strict=False)
)
different_seed_changes = any(
a.overrides != c.overrides for a, c in zip(specs_a, specs_c, strict=False)
)
print(f"Same seed reproduces every per-run override: {same_seed_identical}")
print(f"Different seed yields a different draw set: {different_seed_changes}")
print()
print("First three runs at SEED:")
for spec in specs_a[:3]:
pretty = {k: round(v, 6) for k, v in spec.overrides.items()}
print(f" run_id={spec.run_id}: {pretty}")
Same seed reproduces every per-run override: True
Different seed yields a different draw set: True
First three runs at SEED:
run_id=0: {'CoastTime.Value': 615.622936, 'Inj.Element1': 1.00177, 'Inj.Element2': 0.002277, 'Inj.Element3': -0.002967}
run_id=1: {'CoastTime.Value': 625.65666, 'Inj.Element1': 1.002615, 'Inj.Element2': 0.001518, 'Inj.Element3': 0.00302}
run_id=2: {'CoastTime.Value': 622.309289, 'Inj.Element1': 0.999758, 'Inj.Element2': -0.002852, 'Inj.Element3': -0.007749}
Median band over time with sweep_summary¶
sweep_summary collapses the parent (run_id, time)-frame across runs at each time step, returning a per-time statistics frame with a two-level ("statistic", "field") column index. The default q=(0.05, 0.5, 0.95) matches the standard 5/50/95 dispersion bands.
sweep_band_plot reads that summary, picks the lowest and highest q* for the requested column, and shades the band around the median. The plot below is Sat.X over the propagation — the dispersion of the X component of the inertial position at each time step.
from gmat_sweep import sweep_summary
from gmat_sweep.plotting import sweep_band_plot
summary = sweep_summary(df, by="time")
fig, ax = plt.subplots(figsize=(8, 4))
sweep_band_plot(summary, "Sat.X", ax=ax)
ax.set_title("Sat.X median and 5-95% band over time across 1000 runs")
ax.legend(loc="best")
fig.tight_layout()
plt.show()
Pair plot of perturbations vs. miss distance with sweep_corner¶
sweep_corner draws a lower-triangular pair plot of the perturbed dotted-paths, coloured by a per-run scalar metric. Diagonal panels are histograms of each perturbed axis.
The perturbed values aren't columns of df (sweep outputs carry only the report columns), so we materialise them onto df once via expand_monte_carlo_to_run_specs — the same function the determinism cell below uses. The metric is a callable that derives miss-from-mean per run.
import pandas as pd
from gmat_sweep.plotting import sweep_corner
specs = expand_monte_carlo_to_run_specs(
perturb, n=N, seed=SEED, script_path=script_path, output_dir=Path(".")
)
overrides_per_run = {spec.run_id: spec.overrides for spec in specs}
for param in perturb:
df[param] = df.index.get_level_values("run_id").map(
lambda rid, p=param: overrides_per_run[rid][p]
)
def miss_distance_series(df_in):
final = df_in.groupby(level="run_id").last()
arrivals = final[["Sat.X", "Sat.Y", "Sat.Z"]].to_numpy()
miss = np.linalg.norm(arrivals - arrivals.mean(axis=0), axis=1)
return pd.Series(miss, index=final.index, name="miss_km")
axes = sweep_corner(df, params=list(perturb), metric=miss_distance_series, s=8, alpha=0.5)
axes[0, 0].figure.suptitle("Perturbation pair plot, coloured by arrival miss (km)", y=1.02)
plt.show()
Monte Carlo convergence diagnostic¶
mc_convergence walks the runs in run_id order and reports the running mean, running standard deviation, and standard error of the mean over the first n = 1..N successful runs. A flattening centre line and a shrinking band indicate that the per-run scalar has converged.
The arrival-miss-from-mean is a global reduction (it depends on every run), so it can't be expressed as a per-run callable. Materialise it as a miss_km column on df first — mc_convergence then uses the column-name code path with terminal_only=True to take each run's last value.
mc_convergence_plot renders the running mean as a line and the ±1·SE envelope as the shaded band.
from gmat_sweep import mc_convergence
from gmat_sweep.plotting import mc_convergence_plot
miss_series = miss_distance_series(df)
df["miss_km"] = df.index.get_level_values("run_id").map(miss_series.to_dict())
conv = mc_convergence(df, "miss_km", terminal_only=True)
fig, ax = plt.subplots(figsize=(8, 4))
mc_convergence_plot(conv, ax=ax)
ax.set_title("Running mean miss distance with ±1 SE envelope")
ax.set_xlabel("Number of runs")
ax.set_ylabel("Miss distance (km)")
fig.tight_layout()
plt.show()
Where to next¶
- Variance reduction at small
n. Notebook 05 compares plain Monte Carlo against Latin hypercube sampling at the sameperturband the samen=64, showing the stratification visually in a pair plot of the unit-cube samples. - Manifest-driven replay. Pass
out=Path("./mc-run")to anchor the per-run Parquet files and themanifest.jsonlunder an explicit directory;Manifest.loadreads the per-run overrides back without re-running. - Resume after a failure. Notebook 03 walks through a kill mid-sweep followed by a programmatic
Sweep.from_manifest(...).resume()that re-runs only the missing cells. - Distribution shapes outside the three shorthands. Parameter spec → Stochastic specs documents the pre-frozen
scipy.statspass-through for triangular, beta, truncated normal, and so on.