Closed-Form Solutions to the Fokker-Planck Equation for Orbital Uncertainty Propagation

This paper presents a computationally efficient, closed-form, grid-free solution to the Fokker-Planck equation that propagates orbital uncertainty under stochastic forcing by preserving an exponential-of-quadratic-form ansatz, thereby accurately capturing non-Gaussian features and asymmetric tails without the need for Monte Carlo sampling.

The Gist

Imagine you are trying to predict where a lost hiker will be tomorrow.

In the world of space, instead of a hiker, we have a satellite. Instead of a forest, it's orbiting Earth. The problem is that we never know the satellite's exact position or speed; there's always a little bit of "fuzziness" or uncertainty. Furthermore, the air around Earth (even high up) isn't perfectly smooth; it has tiny, random gusts that push the satellite off course. This is called stochastic forcing (random noise).

The Old Way: The "Gaussian" Guess

Traditionally, scientists have treated this uncertainty like a perfect, round balloon (a Gaussian distribution). They assume the satellite is most likely in the center, and the chance of it being far away drops off symmetrically, like a bell curve.

The Problem: Space isn't a straight line. Gravity bends paths, and satellites speed up and slow down wildly as they get closer to Earth. These forces stretch and twist that "balloon" into weird, banana-shaped clouds.

• If you use the old "round balloon" math, you might think the chance of the satellite hitting space junk is zero because the "balloon" doesn't reach the junk.
• In reality, the "tail" of the cloud (the weird, stretched-out part) might actually be touching the junk.
• To get the right answer using old methods, you'd have to run a supercomputer simulation with millions of fake satellites (Monte Carlo method) just to see if any of them hit the target. This takes forever and uses a lot of power.

The New Way: "Taylor Map Diffusion"

This paper introduces a clever new trick called Taylor Map Diffusion. Instead of guessing or running millions of simulations, the authors found a mathematical "shortcut" that solves the problem exactly (or very close to it) using a single, elegant formula.

Here is how it works, using an analogy:

1. The Stretchy Rubber Sheet (The Map)

Imagine the uncertainty of the satellite's position as a rubber sheet.

• The Old Way: You just stretch the sheet evenly in all directions.
• The New Way: The authors realized that even when the rubber sheet gets twisted, bent, and stretched into a complex shape by gravity and wind, it still follows a specific mathematical pattern. They call this an "exponential-of-quadratic-form."

Think of it like this: If you take a piece of dough and stretch it, it gets thin in some places and thick in others. The authors proved that you don't need to track every single grain of flour in the dough. You only need to track how the dough is being stretched and twisted at a few key points.

2. The "Recipe" (The ODEs)

The paper proves that you can describe this entire stretching process using a small set of rules (equations) that tell you how the "shape" of the uncertainty changes over time.

• Instead of tracking 1,000,000 fake satellites, their method tracks 94 numbers (variables) that describe the shape of the cloud.
• These numbers update themselves automatically as the satellite flies, accounting for gravity's pull and the random wind gusts.

3. The "Tail" Problem

The most important part of this paper is that it captures the tails of the distribution.

• In a normal bell curve, the tails are thin and predictable.
• In space, the tails can get huge and lopsided (asymmetric).
• Because this new method tracks the shape of the rubber sheet, it knows exactly how far the "tail" stretches. It can tell you, "There is a 1 in 100,000 chance the satellite is right here," even if that spot is very far from the center.

Why This Matters

• Speed: The old way (Monte Carlo) takes minutes or hours to run a simulation with enough samples to be accurate. This new way takes less than a second on a standard laptop.
• Accuracy: It doesn't just guess the average; it captures the weird, stretched-out shapes that happen in real life.
• Safety: In space, knowing the "tails" is crucial. If you miss a collision warning because your math was too simple, you could lose a satellite worth millions of dollars. If you overestimate the risk, you waste fuel maneuvering unnecessarily. This method gives the right answer, fast.

The Bottom Line

The authors found a way to turn a messy, chaotic problem (predicting a satellite's path with random wind and gravity) into a clean, solvable puzzle. They proved that the "shape" of the uncertainty follows a specific rule that can be calculated instantly, without needing to simulate millions of possibilities. It's like having a crystal ball that tells you exactly where a lost hiker might be, including the unlikely places, without having to send out a search party of a million people.

Technical Summary

1. Problem Statement

The paper addresses the critical challenge of orbital uncertainty propagation in space operations, specifically for conjunction assessment (collision avoidance) and space domain awareness (SDA).

• The Core Issue: Current operational methods rely on Gaussian assumptions propagated through linearized dynamics (e.g., covariance propagation). However, nonlinear orbital dynamics (such as the $1/r^2$ gravitational force) distort initially Gaussian uncertainty clouds into non-Gaussian shapes with asymmetric tails.
• The Consequence: Gaussian models fail to accurately represent the "tails" of the probability distribution. Since collision probabilities are often extremely low (e.g., $10^{-5}$), errors in modeling these tails lead to either false alarms (wasting fuel on unnecessary maneuvers) or missed warnings for dangerous encounters.
• Limitations of Existing Methods:
  • Monte Carlo (MC): While accurate, MC requires $10^5$ to $10^6$ samples to resolve low-probability tails, making it computationally prohibitive for real-time operations.
  • Gaussian Mixture Models (GMM): Require an increasing number of components to capture non-Gaussian tails, leading to high computational costs.
  • Grid-based Solvers: Suffer from the "curse of dimensionality" (e.g., a 6D state on a coarse grid requires $\approx 10^{10}$ evaluations per step).
  • Probability Transformation Method (PTM): Excellent for deterministic systems but struggles to incorporate process noise (stochastic forcing like atmospheric drag fluctuations) without reintroducing high-dimensional convolution integrals.

2. Methodology: Taylor Map Diffusion

The authors propose a new framework called Taylor Map Diffusion, which provides a closed-form, grid-free solution to the Fokker-Planck Equation (FPE) for systems with additive Gaussian noise.

A. Theoretical Foundation

The method is based on proving that a specific ansatz (structural form) for the probability density function (PDF) is preserved under the combined action of nonlinear advection (deterministic drift) and diffusion (stochastic noise).

• The Ansatz: The PDF $p(x, t)$ is modeled as the exponential of a quadratic form of a nonlinear map:
  $$p(x, t) = \mathcal{N}(t) \mu(x, t) \exp\left( -F(x, t)^\top Q(t) F(x, t) \right)$$
  Where:
  • $F(x, t)$: A smooth, invertible map (represented as a Taylor expansion) centered at the distribution's mode.
  • $Q(t)$: A symmetric positive-definite precision matrix (inverse covariance).
  • $\mu(x, t)$: A scalar field accounting for volume divergence (constant for conservative systems).
  • $\mathcal{N}(t)$: A normalization constant.

B. Key Theoretical Result (Theorem 1)

The authors prove that if the initial density follows this form, it remains in this form over time. The evolution of the parameters ($F, Q, \mu$) is governed by a coupled system of Ordinary Differential Equations (ODEs):

1. Map Evolution ($F$): Describes how the shape of the distribution deforms due to advection and diffusion.
2. Precision Matrix Evolution ($Q$): Describes how the distribution broadens (uncertainty grows) based on the local curvature of the flow and noise intensity.
3. Volume Evolution ($\mu$): Accounts for compressibility (relevant for non-conservative forces like drag).

Crucially, these ODEs are derived by direct substitution into the FPE, yielding zero residual. No approximation is made at the continuous level; the only approximation in practice comes from the truncation order of the Taylor expansion used to represent the map $F$.

C. Implementation

• Representation: The map $F$ is represented as a finite-order Taylor expansion (e.g., up to second order).
• State Vector: The augmented state includes the nominal trajectory, the Jacobian, the Hessian tensor, and the precision matrix. For a 4D planar orbit, this results in a system of 94 scalar ODEs.
• Computation: The system is integrated using standard ODE solvers (e.g., RK4). Derivatives are computed efficiently using automatic differentiation (via JAX), avoiding complex manual derivation of high-dimensional analytical expressions.

3. Key Contributions

1. Structural Preservation Proof: The paper establishes, for the first time, that the exponential-of-quadratic-form solution is structurally preserved under the full FPE dynamics (nonlinear advection + additive diffusion).
2. Grid-Free Non-Gaussian Propagation: The method captures the full non-Gaussian shape of the PDF, including asymmetric tails, without spatial discretization or sampling.
3. Closed-Form Tail Evaluation: Unlike Monte Carlo methods, the solution provides an analytical expression for the PDF. This allows for the direct evaluation of collision probabilities in the deep tails (e.g., $5\sigma$ or $10\sigma$) with negligible computational cost.
4. Unified Framework: It generalizes classical covariance propagation (Kalman filter) to nonlinear dynamics, where the precision matrix $Q$ and the nonlinear map $F$ replace the simple mean and covariance.

4. Results

The method was validated on a planar Keplerian orbit problem with the following setup:

• Scenario: An eccentric orbit ($a \approx 0.74$) propagated for 9.5 revolutions ($t_f \approx 37.7$) with additive stochastic noise on velocity (representing drag fluctuations).
• Benchmark: Compared against a Monte Carlo simulation with 400,000 stochastic trajectories.
• Performance:
  • Accuracy: The Taylor Diffusion solution showed close agreement with the Monte Carlo benchmark across position and velocity subspaces. It successfully captured banana-shaped distributions, asymmetric tails, and stochastic broadening caused by the nonlinear $1/r^2$ dynamics.
  • Efficiency: The Taylor Diffusion solver integrated 94 ODEs in less than one second. In contrast, the Monte Carlo benchmark required over a minute and still suffered from sampling noise.
  • Scalability: The authors note that extending this to a full 6D orbital state would require only 279 ODEs, which remains computationally trivial compared to the exponential growth of Monte Carlo samples needed for tail resolution.

5. Significance and Future Outlook

• Operational Impact: This method offers a computationally efficient alternative to Monte Carlo for conjunction assessment, enabling accurate collision probability estimation in the critical low-probability tails without the prohibitive cost of sampling.
• Orbit Determination: The framework provides a non-Gaussian prior that can be updated analytically, offering a robust solution for sequential orbit determination in scenarios with sparse measurements and long propagation gaps where Gaussian filters (EKF/UKF) typically degrade.
• Future Work: The authors identify extensions to non-conservative forces (requiring the evolution of $\mu$), higher-order Taylor expansions for longer horizons, and the application of Generalized Equinoctial Orbital Elements (GEqOE) to reduce effective nonlinearity in the state representation.

In summary, Taylor Map Diffusion bridges the gap between the accuracy of Monte Carlo methods and the speed of analytical covariance propagation, providing a rigorous, closed-form tool for managing uncertainty in complex, stochastic orbital environments.