Analytical Confidence Boundaries for Non-Gaussian Uncertainty in Perturbed Spacecraft Dynamics

Abstract: This work investigates nonlinear uncertainty propagation in perturbed astrodynamics, focusing on the rapid characterization of non-Gaussian distributions and the construction of three-dimensional "banana-shaped" confidence boundaries. To bridge the gap between computationally intensive high-fidelity methods and inaccurate linear approximations, this paper introduces a fully analytical, sample-free framework for higher-order moments extraction. Leveraging Differential Algebra to bypass repeated numerical integration, statistical moments are extracted analytically via Isserlis' theorem and a monomial-to-Hermite basis transformation. A pair-product projection strategy is exploited to overcome the severe computational bottleneck of full fourth-order tensor contractions and compute only relevant terms via efficient polynomial algebra. The extracted skewness and kurtosis components directly parameterize non-elliptical confidence geometries that capture spatial bending and out-of-plane coupling of typical non-Gaussian distributions in astrodynamics. The approach is validated in high-fidelity environments including a cislunar Near-Rectilinear Halo Orbit and close-proximity trajectories around Apophis during Earth's flyby, where the analytical approach achieves geometric accuracy comparable to expensive Monte Carlo simulations while reducing computational runtime by orders of magnitude.

• The paper presents a method for analytically computing higher-order uncertainty boundaries in spacecraft dynamics using sample-free techniques.
• It employs Differential Algebra and Isserlis’ theorem to extract skewness and kurtosis, forming accurate 'banana-shaped' confidence contours.
• Validation in cislunar and Apophis scenarios demonstrates improved sample containment and reduced computational cost compared to Monte Carlo and LinCov methods.

Analytical Confidence Boundaries for Non-Gaussian Uncertainty in Perturbed Spacecraft Dynamics

Problem Definition and Motivation

Quantifying uncertainty in spacecraft state propagation is fundamental for mission robustness, especially when dynamics exhibit strong nonlinearity, weak stability, or proximity to multi-body gravitational resonances. Standard covariance-based geometric bounding is inadequate in these conditions, as initially Gaussian uncertainties rapidly evolve into complex non-Gaussian distributions characterized by significant skewness, kurtosis, and spatial "banana" geometries. While Monte Carlo (MC) methods provide high-fidelity reference statistics, their computational cost prohibits routine large-scale analyses or onboard evaluation, and even advanced sampling-based schemes (e.g., GMMs, PCE, CUT, UT) still entail repeated, expensive integrations. Traditional linear covariance approximations (LinCov) can fail catastrophically in bounding state dispersion and uncertainties in these settings.

This work systematically addresses the construction of fully analytic, sample-free, higher-order uncertainty boundaries that can efficiently characterize generic non-Gaussian features in perturbed orbital dynamics, particularly for cislunar and small-body environments where such phenomena are acute. The principal advance is the derivation and validation of closed-form, 3D confidence boundaries incorporating analytically computed skewness and kurtosis, informed by efficient exploitation of Differential Algebra (DA), Isserlis' theorem, and polynomial chaos expansions.

Analytical Framework for Confidence Boundary Construction

The analytic process synthesizes multiple mathematical innovations:

• Basis Construction: Propagation of initial Gaussian perturbations is carried out with high-order DA, which yields the dynamical flow as a multivariate polynomial, as opposed to sampling the ODE repeatedly. Expansion is handled in whitened coordinates, which diagonalizes initial covariance.
• Moment Extraction: Statistical moments (mean, covariance, skewness, kurtosis) are extracted directly from these polynomials via Isserlis' theorem, which yields expectations without sampling. For higher-order moments, the computational bottleneck in direct tensor contraction is addressed by a Hermite basis transformation, which sparsifies the moment tensor.
• Project-then-Product Strategy: Recognizing that only certain moments (e.g., those aligned with principal directions of stretching and bending) are required for geometric boundary construction, the authors introduce a method in which polynomial maps are first projected onto principal axes, after which moment computation reduces to 1D polynomial algebra (circumventing full tensor formation and contraction; see Section 3).
• 3D Banana Boundary Synthesis: The principal axes (length, width, thickness) of the final uncertainty distribution are identified by eigendecomposition of the covariance matrix. The analytic boundary is constructed by deforming an initial ellipsoid according to the projected higher-order moments: in-plane bending (cross-skewness), tail asymmetry (skewness, Cornish-Fisher expansion), and out-of-plane torsion (out-of-plane skewness), producing the physically observed "banana-shape" uncertainty contours.

Results depend critically on accurate evaluation of terms such as $E_{uuu}$, $E_{uuv}$, $E_{uuw}$, and $E_{uuuu}$, representing higher-order moment projections along and across the major axes.

Figure 1: Cislunar NRHO scenario. Left: Nominal trajectory with uncertainty propagation distribution at different instants. Right: Different viewing angles of the analytical boundary reconstruction informed by higher moments at final time.

Numerical Validation in Perturbed Dynamical Regimes

The method is validated in two difficult problem classes:

• Cislunar NRHO (Earth-Moon $L_2$): The strong nonlinearities and rapid variation around perilune stretch uncertainty contours into highly non-Gaussian distributions. The analytic 3D banana boundary robustly encloses the MC ensemble, effectively capturing spatial bending and out-of-plane deformation.
• Asteroid Apophis Proximity: Both deep-space orbits and Earth-flyby arcs are considered in a regime dominated by rapid, irregular, third-body and higher-order gravitational perturbations. Once again, the analytic boundaries derived from DA-based higher moments tightly fit the MC clouds, outperforming LinCov ellipsoids and Unscented Transform boundaries.

The authors show that ellipsoidal approximations dramatically underestimate the volume required for containment (as little as 69% of MC samples within "3$\sigma$" ellipsoid vs. $>96\%$ for analytic banana). The analytic method reliably achieves containment rates matching the theoretical confidence level of the non-Gaussian boundary.

Figure 2: Comparison of different UQ methods at the final time instant, with respect to the MC solution (gray dots) and the DA analytical approximation (black solid line).

Figure 3: Computational times of different UQ methods for cislunar scenario. Note that DA methods times do not include the DA expansion time, but are related to the single UQ query. DA Full Tensors refers to the faster Hermite-DA approach.

Computational Performance and Scaling

A critical contribution is the drastic reduction in online computational runtime afforded by the analytic, tensor-sparsifying, project-then-product method:

• MC Methods: Require $10^4$--$10^5$ ODE integrations for convergence, resulting in prohibitive timescales.
• GMM, UT, Standard PCE: Reduce the sample count, but still perform large numbers of ODE solves and cannot efficiently extract high-order moments.
• DA Full Tensor (Hermite): Eliminates sampling noise, but the scaling of tensor contractions remains prohibitive in high-dimension, high-order expansions.
• DA Project-Then-Product: Reduces online evaluation to 1D polynomial algebra; online UQ for a new uncertainty distribution is orders of magnitude faster while maintaining accuracy.

As shown in the comparison figures and tables, relative errors in skewness and kurtosis, and thus in the boundary geometry, are reduced to a few percent or less, matching the best possible MC estimates. For example, in the cislunar scenario, the analytic DA methods have covariance, skewness, and kurtosis relative errors of $O(10^{-2})$ without any sampling.

Figure 4: Runtime breakdown for analytical higher-order moments recovery approaches using DA framework.

Figure 5: Statistical moments relative errors with respect to the MC reference for different UQ methods for cislunar scenario. Note that the error of skewness and kurtosis of DA project-then-product approach refers only to the extracted terms.

In computational cost breakdowns, the analytic method's online query time is reduced from $E_{uuv}$0 (for MC) and $E_{uuv}$1 (for full-tensor DA) to sub-second ($E_{uuv}$2) for a 6D, 4th-order expansion, with no sacrifice in boundary fidelity.

Theoretical and Practical Implications

The methodology enables rigorous, real-time, onboard uncertainty quantification suitable for robust GNC, autonomous collision avoidance, and trajectory optimization under realistic, highly perturbed dynamical environments. Analytical confidence boundaries support the rigorous satisfaction of chance constraints in stochastic optimization, even under severe deviation from Gaussian statistics. The analytic nature of the construction allows for gradient-based optimization and formal verification of enclosure.

Furthermore, the framework generalizes to any system in which DA or high-order Taylor expansions are available and is not restricted to astrodynamics. The modular definition of the boundary deformation allows straightforward extension to higher-order corrections, more complex topological uncertainty features, or non-unimodal distributions, if the relevant higher-order moments can be evaluated.

Speculation on Future Developments

Future directions include the extension to:

• Fully non-unimodal (e.g., multi-modal or mixture) distributions using analytic Gaussian mixtures or higher-order chaos expansions.
• Automated error control and adaptive order selection in DA expansion for certified real-time onboard applications.
• Augmentation with state transition tensors, integrating robust DA-based UQ directly into closed-loop stochastic optimization and feedback synthesis.
• Direct integration with stochastic guidance and autonomous onboard decision-making under non-Gaussian bounded uncertainty.

Conclusion

This work establishes a sample-free, analytic, and computationally scalable framework for propagating non-Gaussian uncertainty and for constructing rigorously justified confidence boundaries in highly perturbed nonlinear dynamical systems, with specific application to astrodynamics. By exploiting sparse-tensor moment extraction and axis projection, the method achieves both the geometric accuracy of high-fidelity MC while reducing computational cost by several orders of magnitude. The methodology directly supports high-tempo, onboard, and large-scale robust design and operation activities where conventional methods are computationally prohibitive, and enables new avenues in analytic, real-time, and provable uncertainty management for autonomous spacecraft operations.

Figure 6: Apophis proximity scenarios. Left: During deep-space. Right: During Earth's flyby.

Figure 7: Computational times of different UQ methods for Apophis (Earth's flyby scenario). Note that DA methods times do not include the DA expansion time, but are related to the single UQ query. DA Full Tensors refers to the faster Hermite-DA approach.

Figure 8: Statistical moments relative errors with respect to the MC reference for different UQ methods for Apophis (Earth's flyby scenario). Note that the error of skewness and kurtosis of DA project-then-product approach refers only to the extracted terms.