University of California San Diego Research
Stationary/Non-stationary Probabilistic Search
Coming soon...
State Estimation of Chaotic Trajectories:
A Higher-Dimensional, Grid-Based, Bayesian Approach to Uncertainty Propagation
Abstract
The current landscape of orbital uncertainty propagation methods inadequately considers the non-Gaussianity of the state estimation problem for nonlinear systems. In relatively low-perturbed regimes, or when state observations are frequent, the assumption of Gaussianity holds, but as mission design begins to exploit the chaoticity of N-body dynamics to efficiently explore new regimes of space, this assumption is often violated. The uncertainty propagation methods that do not assume Gaussianity are computationally expensive when observations are frequent and often disregard epistemic uncertainty, or the uncertainty of the model. To address the current limitations of orbital uncertainty propagation, we introduce higher-dimensional extension to an existing Bayesian estimation algorithm that efficiently propagates the probability distribution function of a nonlinear system. By adjusting the computational architecture of the algorithm and considering the dynamics of the system, we scale the existing, three-dimensional technique with poor time complexity to an efficient, four-dimensional one. The result is a robust, second-order accurate, time-adaptive, explicit time-marching scheme with the capability of propagating uncertainty governed by chaotic, nonlinear dynamics.
Introduction
Since the dawn of space exploration, orbital uncertainty propagation has seen relatively little change. Much of the success of the Apollo lunar missions can be attributed to the work of Kalman, specifically the Extended Kalman Filter [1] (EKF), the primary algorithm used by on-board systems that determined accurate state estimates of the Apollo spacecraft [2]. This foundational algorithm ensured the feasibility of the complex trajectories essential to the Apollo lunar missions. Since the Space Race, spaceflight has undergone a complete upheaval in almost all facets. Advances in additive manufacturing have resulted in fully 3D-printed rocket engines and propellant tanks that are more weight- and cost-efficient than their legacy counterparts [3,4]. Low-thrust relative transfers have motivated the creation of novel low-thrust propulsion methods, such as ionic propulsion [5] and solar-sail systems [6]. There are numerous domains of spaceflight that are almost unrecognizable when compared with their ancestor systems. Conversely, one of the most recent NASA missions, Artemis-1, employed the Orion Absolute Navigation System which utilized four navigation EKFs for estimating the vehicle's state and associated uncertainty [7], the same framework used in the Apollo missions over 60 years ago. In a field where rapid growth is the norm, this sort of stagnation is unprecedented.
By no means is the prior paragraph meant to discredit the EKF. On the contrary, its persistent utilization over the last several decades in the fields of orbit determination, conjunction analysis, GN&C, SSA, etc. serves as a testament to its robustness. However, there are limitations to this robustness. For instance, the EKF approximates the dynamics of a nonlinear system by linearizing about an estimate. This linearization can result in the EKF underestimating uncertainty and depending heavily on precise initial estimates. Additionally, the EKF assumes the uncertainty of an object may be described by a multivariate Gaussian distribution. In the absence of measurement updates, this holds for a period of time, but the length of this period is dependent on the dynamics of the system. Both of these limitations are predicated on the assumption that the measurement update frequency will be high enough that the errors associated with linearizing are negligible and the state uncertainty stays Gaussian. However, the validity of this assumption is brought into question in the presence of chaotic dynamics, as is the case when the motion of an object with negligible mass is governed by two massive bodies, known as the restricted three-body problem [8] (R3BP). Novel space mission design techniques aim to exploit the chaoticity of three- and four-body dynamics via new families of low-energy trajectories that are more fuel efficient, longer in duration, and further in reach than their classical counterparts [9]. CAPSTONE, a CubeSat pathfinder for the Lunar Gateway and a cornerstone of NASA's Artemis Program, utilized a low-energy, ballistic lunar transfer [10] (BLT), a family of these low-energy transfers, to insert into its operational near-rectilinear halo orbit (NRHO) about the L2 Earth-Moon libration point in November, 2022 [11]. Future proposed missions to the Jovian and Saturnian moons [9,12] are expected to utilize the N-body dynamics of the respective systems for low-energy trajectory design. As the standard for mission complexity increases to incorporate these novel trajectory designs, so too does the necessity for effective state estimation and uncertainty propagation techniques that address the limitations of the EKF and provide more precise solutions in these chaotic regimes.
There are other uncertainty propagation methods that attempt to address the shortcomings of the EKF, each with their own limitations. Monte Carlo (MC) simulations, for instance, circumvent the assumption of Gaussianity by generating N random samples from the initial uncertainty and time-marching them, governed by the dynamics of the system [13]. As N goes to infinity, the final distribution of the samples approaches the true probability distribution. However, obtaining this true distribution requires an unknown-but-significant number of samples, and, in especially chaotic regimes that are extremely sensitive to initial conditions, achieving this number may be computationally expensive. Discrete measurement updates require a complete resampling of the MC simulation, resulting in an even larger computational burden when measurements are frequent. Additionally, MC simulations currently do not have a viable way of handling epistemic uncertainty [14], the uncertainty of the model, an important consideration for unexplored, chaotic regimes. Another family of techniques that attempts to rectify the limitations of the EKF are Gaussian Mixture Models [13] (GMMs). GMMs approximate an arbitrary PDF as a collection of weighted Gaussian density functions, then propagate the mean and covariance of each of the distributions in the collection. However, when the true probability distribution ``bananas'' about a nominal trajectory, as has been demonstrated for nonlinear systems [15], many Gaussian distributions are necessary for representing the true distribution, leading to a computational bottleneck. Additionally, determining and updating the weights of the Gaussian mixtures, as well the splitting procedure, can be both expensive and ad hoc.
In summary, having performed a preliminary review of the landscape of methods that represent and propagate orbital uncertainty, we have determined that there exists a need for a new class of methods that do not assume Gaussian uncertainty, are accurate for long periods of time in the absence of observations, are computationally efficient, consider epistemic uncertainty, and are scalable to high-dimensional problems. To address this need, we propose an astrodynamical extension to an existing nonlinear state estimation method that efficiently solves the Fokker-Planck equation (FPE), which describes the time evolution of a PDF governed by arbitrary dynamics, known as Grid-based Bayesian Estimation Exploiting Sparsity [16] (GBEES). Bayesian estimation strategies are foundational in the way they address the state estimation problem, but have often been disregarded due to their computational inefficiency [17]. GBEES addresses this issue by considering the sparsity of a PDF over the majority of phase space, thereby cutting down on the computational expense of numerically solving the FPE over a large grid. However, the legacy implementation is susceptible to computational burden at higher dimensions given its O(n^2) time complexity (where n is the number of grid cells that make up the discretized PDF), making it sub-optimal for orbital uncertainty propagation. The extension we propose reduces this time complexity via efficient data structures and a consideration of the dynamics of the system. We demonstrate the method on a 4D astrodynamical system, and layout the possibility for expanding to higher-fidelity models.
Grid-Based Bayesian Estimation Exploiting Sparsity
To ensure full understanding of the proposed extension, we provide a succinct explanation of the parent algorithm, GBEES (a complete explanation is included in Ref. 16). Consider the state estimation of a nonlinear system
(1)
where x is the state of the system, f(·) is the nonlinear system function, w is the state disturbance, y is the measurement of the system, h(·) is the measurement function, and v is the measurement noise. The evolution of the PDF p_x(x′, t) of the state x at time t is performed via mixed continuous/discrete time-marching [18] and can be described in two steps:
1. Between discrete measurements, p_x(x′, t) is marched via discretization of the Fokker-Planck equation [19,20] (FPE), in Einstein summation notation:
(2a)
where q_ij is the (i, j)th element of the spectral density of the state disturbances, Q. When Q = 0, the PDE is hyperbolic.
2. At measurement interval t_k, p_x(x′, t) is updated via Bayes’ Theorem [21]:
(2b)
where p_x(x′, t_k+) is the a posteriori PDF, p_y(y_k|x′) is the distribution associated with the measurement, p_x(x′, t_k−) is the a priori PDF, and C is a normalization constant.
To discretize the hyperbolic form of Eq. (2a) (i.e. when Q = 0), a Godunov-type, finite-volume method [22] is utilized. In 2D (with higher-dimensional cases following as obvious extensions), this takes the form
(3)
where n represents the time step, ∆x and ∆y are the grid widths in the x- and y-directions respectively, p_ij is the probability at grid cell (i, j), and fluxes F and G are defined at the interfaces of the grid, thus the half-step indexing. The resulting grid of Riemann problems are solved numerically at the interfaces of each grid cell for every time step. Accounting for corner-transport upwind and flux-limiting corrections, the result is an explicit, 2nd-order accurate time-marching scheme that propagates the PDF governed by the dynamics of the system f(x). To save on computational cost, the algorithm only accounts for grid cells that have a probability above some threshold pthresh, thoroughly lessening the number of operations necessary for each time step. As the probability flows throughout phase space, new cells are inserted about those who have surpassed the probability threshold, and cells below pthresh without such neighbors are deleted. Thus, most of the PDF is represented with a small percentage of the domain of phase space.
For validation, the GBEES algorithm is demonstrated on a 3D Lorenz attractor [23], a highly chaotic solution set to the Lorenz system where uncertainty quickly transforms from Gaussian to non-Gaussian. For the 3D Lorenz system
(4)
with parameter values σ = 4, b = 1, and r = 48 resulting in a chaotic system. Figure 1 depicts the chaotic dynamics governing the system, as well as the continuous time-evolution of an initially Gaussian PDF in the absence of discrete measurements, with Q = 0 (assuming no uncertainty in the system). As can be seen, the Gaussian distribution quickly becomes non-Gaussian, and even becomes bimodal near the final epoch, a characteristic that is indescribable by the EKF.
Figure 1. The continuous time-evolution of an initially Gaussian PDF in with no discrete measurement updates. (right) The PDF is described by the blue isosurfaces of probability p = 5e−3, p = 5e−4, and p = 5e−5 at times t = 0, t = 0.2, t = 0.4, t = 0.6, t = 0.8, and t = 1. The green line represents the Lorenz attractor, the dynamics of which govern the flow of probability. The black line represents the nominal trajectory, about which the initial uncertainty is estimated. (left) The outermost isosurface (p = 5e−5) compared to a 1000 particle Monte Carlo simulation with the same initial conditions, represented by the gray trajectories and black dots. A complete list of the simulation parameters can be found in Table A1
The right subfigure of Figure 1 demonstrates that the outermost isosurface (p = 5e−5) of the PDF propagated by GBEES captures the results of the Monte Carlo simulation starting with the same initial conditions and covariance. Due to its extreme nonlinearity, describing the final PDF using a GMM would require a computationally expensive number of Gaussian distributions.
Addressing the O(n^2) time complexity
The primary issue with the legacy GBEES implementation arises as we push towards higher dimensions. Scaling to higher-dimensional systems requires discretizing a larger domain of phase space, leading to increases in the number of cells needed to represent a PDF. As cells are created and deleted based on the probability threshold constraint, neighboring cells must be adjusted to consider the information provided by the new discretization. In the legacy implementation, said adjustments are performed through nested loops, leading to a time complexity of O(n^2) for the entire algorithm, resulting in the inefficiency at higher dimensions. The aim of the proposed extension is to eliminate this computational bottleneck, as well as look for subprocedures that can be made more efficient. As this is a computational limitation, the following discussion is over the modifications and adjustments implemented to optimize the algorithm’s architecture.
Binary Search Trees
Bayesian estimation methods are predicated on the ability to efficiently perform the grid updates associated with time-marching the numerical solution to Eq. (2a). As such, utilizing efficient data structures within the chosen finite volume method is of the utmost importance. The legacy implementation of GBEES stores the discretized Cartesian grid in a list, which has an O(n) time complexity for creation, deletion, and searching. As the algorithm countlessly uses all three of these procedures throughout, improving their individual time complexities is greatly beneficial. One such data structure that has a better time complexity (O(log n) for creation, deletion, and searching) are binary search trees (BSTs).
Introduction to BSTs. BSTs [24] sort data using a positive integer key z. Each datum, known as a node, stores its own key value as well as pointers to the left and right child nodes. The left child has a value z less than the parent node, and the right child has a value z greater than the parent node. Searching a BST to see if a key exists involves first comparing the key being searched for to the root node, or the first node in the tree, then traversing right or left down the tree depending on if the key being searched for is less than or greater than the root node. This process is repeated until the key being search for is reached, or a leaf node is reached, a node with no left or right child. In this case, the key being searched for does not exist in the BST. Creation and deletion of nodes follows in a similar fashion.
Storing grid cell coordinates in a BST requires a conversion from the real coordinate set to the positive key value z, an unexpectedly nontrivial process; consider that the conversion function r(·) must be able to discern between commutative coordinates (i.e. in 2D, r(a, b) does not equal r(b, a) as coordinate (a, b) does not equal (b, a) when a does not equal b). Additionally, the conversion function must be bijective, or establish a one-to-one mapping from the coordinate set to the key value and cover the entire domain of possible key values, to make for an efficient, compact conversion. Of these conversion functions, known as pairing functions, the most apt at handling high-dimensional coordinates is the Rosenberg-Strong pairing function
Figure 2. Two equivalent BSTs, where the left can be transformed into the right via a right rotation, and the right can be transformed into the left via a left rotation. Each letter represents the key value, thus A < b < C < d < E. Triangle nodes represent nodes that are themselves BSTs.
Rosenberg-Strong pairing/unpairing. The Rosenberg-Strong (RS) pairing function [25] compactly converts n-dimensional coordinates to positive, unique keys by setting the shell number equal to the L∞ norm, or the maximum, of the coordinate set. This way, coordinates on the same shell have similar magnitudes.
Figure 3. RS pairing in two dimensions
In 2D, the RS pairing function r_2(x, y) is
(5a)
and the RS unpairing function in 2D r_2^−1(z) is
(5b)
For n-dimensional coordinates, the RS pairing function r_n(x_1, . . . , x_n) is performed recursively
(6a)
and similarly, the RS unpairing function in n-dimensions r_n^-1(z) is
(6b)
Table 1 demonstrates the importance of considering the L∞ norm for high dimensional coordinate sets, as other pairing functions, like the Cantor or Szduzik functions [25], do not. This can quickly lead to computational bit overflow as the discretized grid advects away from the origin. We utilize the RS pairing function to tractably and conveniently store the information at each grid cell in a BST, with assurance that the high-dimensional coordinates far from the origin will not result in bit overflow
Table 1. Pairing a 6-dimensional coordinate point with various pairing functions
Consideration of the dynamics of the system
In addition to changing data structures, modifications to the creation and deletion subprocedures can be made to improve the efficiency of the algorithm. The creation and deletion of cells that neighbor those with probabilities above the threshold value is essential to exploiting sparsity. Decreasing the number of cells created and deleted at each time step while preserving the true distribution of the PDF results in a speedup at no cost. By considering the dynamics of the system at each grid cell, the number of redundant cells created that are then immediately deleted in the subsequent time
step, and vice versa, are reduced
Creation of cells. In the legacy algorithm, if a cell is above the probability threshold, at the following time step, the algorithm creates all the neighboring cells in all grid directions (assuming they do not already exist) and inserts them into the BST. This ensures that, should the velocity of the system advect the probability in any direction, the change is captured by a neighboring cell. However, this method ignores the fact that the direction of the velocity at any given point is known, thus the algorithm need only create cells in the known downwind direction. Figure 4 schematically demonstrates how this consideration improves the creation procedure.
Deletion of cells. For the deletion process, the legacy implementation checks all directions to ensure that a considered cell neighbors no cells with probability above threshold before deleting it from the tree. Checking every neighbor is an oversight that results in saving useless cells, as only the upwind cells can advect probability into the cell of interest. Therefore, the current implementation cuts down on the number of cells checked along the entire grid for each time step, as is demonstrated schematically in Figure 5.
Additionally, the current algorithm considers the frequency of the deletion procedure. The purpose of pruning low probability cells from the BST is to decrease the total number of Riemann problems that need to be solved. However, the deletion procedure requires an exhaustive search of the entire BST, and in certain cases, the time saved on time-marching fewer cells is not worth the process of deletion. Therefore, the current implementation makes the deletion frequency a parameter, and tunes it empirically.
Figure 4. Two-dimensional schematic demonstrating the difference in the creation procedure of the (left) legacy implementation vs. the (right) current implementation. The green cell represents one with a probability above the threshold, and the downwind fluxes F_i+1/2,j and G_i,j+1/2 represent the probability flow at the half-step forward x- and y-interfaces, respectively.
Figure 5. Two-dimensional schematic demonstrating the difference in the deletion procedure of the (left) legacy implementation vs. the (right) current implementation. The white cell represents one with a probability below the threshold, and the yellow cells are the neighbors checked for probability above the threshold value by the algorithm. The upwind fluxes F_i−1/2,j and G_i,j−1/2 represent the probability flow at the half-step back x- and y-interfaces, respectively.
Time step adaptivity
Stability is an important consideration when designing numerical methods. To ensure stability when time-marching discretized functions through phase space, the size of the time step must adhere to the magnitude of the rate of change, so information does not “skip” grid cells, which would cause the scheme to accumulate errors. The Courant–Friedrichs–Lewy (CFL) convergence condition [26] ensures stability for time integration schemes by defining the time step ∆t as a function of C_max:
(7)
where C_max = 1 for explicit time schemes, and Ω represents the entire domain of phase space where p_x(x′, t) ≥ pthresh. When propagating a numerical scheme governed by chaotic dynamics, there are quiescent, slow-changing periods (|f (x)| is small) and there are chaotic, fast-changing periods (|f (x)| is large). Fixed time step schemes must adhere to the most chaotic period of the entire trajectory, even through the more stable quiescent periods, where larger time steps could achieve stability. The legacy implementation of GBEES uses this inefficient fixed time step method tuned empirically to ensure stability of the scheme, a subprocedure that is also computationally wasteful. Therefore, to address these numerical limitations, we utilize an adaptive time step approach, where ∆t is determined via Eq. (7) at each step
Current vs. legacy implementation
Altogether, the changes proposed to the legacy implementation of GBEES improve the overall time complexity of the algorithm, speeding up the process of propagating a three-dimensional PDF in the Lorenz attractor system as demonstrated by Figure 6. Note that although the number of cells required to describe the PDF is lower for the current implementation due to the creation and deletion improvements, the overall shape of growth is the same. However, the shapes of the program runtime are different, demonstrating that the current implementation has an improved time complexity compared to the legacy counterpart. This efficiency improvement allows the algorithm to feasibly be applied to higher-dimensional problems, specifically propagating orbital uncertainty for nonlinear, chaotic systems.
Figure 6. Comparison of legacy implementation with the current implementation of GBEES. This comparison was performed on the 3D Lorenz attractor from Figure 1.
Planar, circular, restricted three-body problem
The restricted three-body problem considers the motion of an object with negligible mass in a system of two massive bodies that orbit about their respective center of mass, unaffected by the object of negligible mass. Assuming the orbit of the two massive bodies about the center of mass is circular, and the motion of the object is constrained to the plane of orbit, the result is the planar circular restricted three-body problem [8] (PCR3BP). In the barycentric, dimensionless, rotating
coordinate frame, the equations of motion of the object are
(8)
μ is the gravitational parameter of the two-body system, r_1 is the distance to the primary body, and r_2 is the distance to the secondary body
In the PCR3BP, there exist five equilibrium points L_1-L_5, known as Lagrange or libration points [9], that lie in the plane of orbit and are at rest in the synodic frame. L_1, L_2, and L_3 are collinear,
while L_4 and L_5 form equilateral triangles with the two massive bodies. In the vicinity of these equilibrium points, there exist initial conditions that result in quasi-periodic orbits (QPOs) about the libration points. Multiple spacecraft have operated in these orbits [27, 28] as they require minimal propulsion to maintain. The class of quasi-periodic orbits that lie entirely in the plane of the two massive bodies are known as Lyapunov orbits [29]. Extensive research has been conducted to accurately calculate the initial conditions of these orbits [30], to such an extent that precomputed catalogs,
like the JPL Three-Body Periodic Orbit Catalog [31], store the initial conditions for QPOs about all three collinear libration points for various planetary systems. As the initial conditions for these periodic orbits are readily available, and the orbits themselves are low-fidelity representations of true spacecraft motion, we view this as a prudent starting point for the application of GBEES to higher dimensions.
To govern the advection term of Eq. 2a by the equations of motion of the PCR3BP, Eq. 9 must be converted from a 2nd-order ODE to a 1st-order ODE. In 1st-order form, the state x and the equations of motion dx/dt of the PCR3BP become
(9)
Figure 7. Schematic of the Lagrange points in the synodic, dimensionless coordinate frame. Sourced from Ref. 9.
The result is a four-dimensional nonlinear system, a useful and relevant test case with increased complexity for the new and improved GBEES algorithm. To visualize the 4D PDFs generated and propagated, we create two 2D PDFs by integrating p_x(x′, t) over the velocity space for the position PDF p_(x,y)(x′, y′, t), and by integrating p_x(x′, t) over the position space for the velocity PDF p_(v_x,v_y)(v_x′, v_y′, t):
(10a)
(10b)
where Ω(x,v) and Ω(v_x,v_y) represent the entire domain of the PDF over position and velocity, respectively. Utilizing the initial conditions precomputed in the JPL Three-Body Periodic Orbit catalog, we possess a complete framework for propagating a four-dimensional PDF representing the uncertainty of an object in a Lyapunov orbit. We now demonstrate the flexibility of this extension of GBEES given various physical scenarios where we speculate the EKF may perform poorly.
Infrequently observed Jupiter-Europa trajectory
As space missions venture further into the depths of the solar system, novel methods for guidance and navigation become imperative due to the limitations of the legacy standard, wherein commands are uplinked based on downlinked measurements received prior, communicated via the Deep Space Network (DSN). One such limitation is that, at such extreme distances from Earth, data communication between spacecraft and ground control has significant light time delays (i.e., approximately an hour and a half at Jupiter) that complicate navigation in environments where uncertain dynamics require trajectory control on timescales under the light time delay (i.e., repeated low altitude flybys of a Jovian or Saturnian moon). Second, a non-autonomous approach requires a reserved DSN antenna for uplink and downlink, which may not always be feasible as the number of deep space missions grows [32, 33]. NASA’s Deep Space 1 [34, 35] (DS1) avoided these limitations by employing an autonomous navigation system, AutoNav, that utilized on-board imagery of near asteroids to autonomously determine its trajectory in space and compute necessary maneuvers (a.k.a. OpNav). One drawback of utilizing OpNav using main belt asteroids is that cataloged asteroids beyond the main belt are fewer and far between, making measurements updates less frequent (DS1 planned to utilize AutoNav at roughly weekly intervals36). In chaotic regimes, this will result in the uncertainty of a spacecraft spreading widely over phase space. To maintain custody of spacecraft operating in QPOs in these regimes using autonomous navigation, operational orbits must be relatively stable
such that infrequent measurements will not result in accidental reentry or escape.
Given the physical scenario described, the nominal trajectory chosen for uncertainty propagation will be relatively stable and discrete measurement updates will be infrequent. The EKF assumes that the errors in estimation due to the utilization of a linearized propagation model are corrected via frequent measurement updates [37]. However, in the trajectory we have outlined, measurement updates partially rely upon the presence of nearby cataloged asteroids, thereby diminishing our confidence in the assumption of frequency. Conversely, conducting uncertainty propagation via a MC simulation requires an unknown-but-substantial number of particles to ensure a sufficient level of confidence in the approximated distribution. Without knowledge of the posterior distribution prior to the following measurement update, defining what qualifies as “substantial” may be fairly unpredictable. For normal distributions, there are formulae following from the Central Limit Theorem that provide the number of particles required to achieve some confidence level of the final,
approximated distribution [38]. However, for non-Gaussian distributions, determining this optimal number is less straightforward [39]. This often leads to an excessive number of initialized particles, ultimately squandering computational resources.
We utilize a set of initial conditions that result in a Lyapunov orbit about the L_3 Jupiter-Europa libration point, sourced from the JPL Three-Body Periodic Orbit catalog (initial conditions can be found in Table 2 and are relative to the Jupiter-Europa barycenter). We assume an initial uncertainty in position of 10^3 km and velocity of 10^−1 km/s, as is the assumed a priori uncertainty for spacecraft near Jupiter utilizing autonomous navigation [40]. To test both the accuracy and efficiency of GBEES, we initialize a MC simulation with the same initial conditions and uncertainty. We propagate the uncertainty using both GBEES and the MC simulation through one full period of the Lyapunov orbit (∼3.5 days), with a discrete measurement update every 1/3 of the orbit (∼1.17 days). Measurements are assumed to be the state of the nominal trajectory at the given epoch, with uncertainty in position and velocity of the measurement being equal to the initial uncertainty.
Table 2. Initial Conditions of L_3 Jupiter-Europa Lyapunov Trajectory
Figure 8 depicts the results of the described L_3 Jupiter-Europa Lyapunov trajectory simulation. The initially Gaussian uncertainty beginning from the M1 epoch rapidly spreads over phase space the uncertainty in y goes from being on the order of 10^3 km to 10^5 km after only about 1.17 days). Via comparison with the MC simulation, we validate that GBEES has accurately propagated the uncertainty between measurement updates
Figure 8. Initially Gaussian uncertainty time-marched with GBEES governed by the PCR3BP in the Jupiter-Europa System. (top) The contour plots represent the PDF propagated by GBEES about the nominal trajectory. M1, M2, and M3 denote the epochs where discrete measurement updates take place, with M1 being the initialization of the simulation. Note that the colors of the contours representing the magnitude of the probability of a PDF are not relative to other PDFs, but instead are representative of probability differences throughout the individual distributions. (bottom) GBEES compared with a 500 particle MC simulation, to confirm the accuracy of the method. A complete list of the simulation parameters can be found in Table A1.
Figure 9 provides the computational efficiency comparison of a RK4 time-marching scheme MC simulation with GBEES over the trajectory of interest. To ensure the MC simulation provides the same “quantity of information” at each measurement epoch M_i for i = 1, 2, 3, we set the number of particles in the MC simulation equal to number
of cells with probability above threshold in the final PDF at M_i for i = 1, 2, 3. Unlike GBEES, which propagates a discretized grid that grows as the uncertainty grows, MC simulations have no mechanism for refining the distribution it is propagating as it spreads over phase space. Therefore, the MC simulation has the disadvantage of having to time-march a large quantity of particles even when the uncertainty is relatively small so that it may represent the uncertainty when it is relatively large. The discrete changes in the number of cells for both methods in Figure 9 is representative of discrete Bayesian updates in the case of GBEES and resamplings from the new measurements in the case of the MC simulation
Figure 9. Comparing the computational time of GBEES vs. MC when propagating uncertainty for the trajectory from Figure 8. Discrete jumps in the number of cells tracked by each uncertainty propagation method are representative of discrete Bayesian updates in GBEES’s case and resamplings in the MC simulation’s case
Consideration of epistemic uncertainty for solar probes
Of the limitations of existing uncertainty propagation methods, the most concerning for space missions delving into unexplored domains may be their failure to consider epistemic uncertainty. Epistemic uncertainty, or uncertainty of a model, may be negligible in regimes where perturbations from solar radiation pressure, atmospheric drag, zonal harmonics, etc. are either well-defined or relatively small. However, this may not be true for all realms of space. Consider the Parker Solar Probe,41 a NASA spacecraft designated with the task of flying through the Sun’s upper atmosphere to take measurements and expand our understanding of solar wind. In this unexplored regime (the Parker Solar probe flies seven times closer to the Sun than any spacecraft has ever flown), the effect of solar radiation pressure on the motion of the spacecraft is not well-defined, thus the model representing the motion of the probe some uncertainty. To account for this, Eq. (2a) includes the diffusion term (i.e. Q does not equal 0) which is representative of random motion. In 2D (with higher- dimensional cases following as obvious extensions), the diffusion term is added to Eq. (3) by updating the fluxes such that, for all (i, j)
(11a)
(11b)
where μ = [μ_x μ_y] is the coefficient of diffusion and can be tuned depending on how uncertain the model definition is.
Again we choose a set of initial conditions that result in a Lyapunov orbit about the L_3 libration point, but of the Sun-Earth system, (initial conditions can be found in Table 3 and are relative to the Sun-Earth barycenter). This specific set of initial conditions from Table 3 takes the spacecraft to a distance of nearly 7 million km from the Sun at perihelion (the Parker Solar Space probe’s closest approach is 6.16 million km from the Sun). Considering epistemic uncertainty is non-negligibile, we include the diffusion term in the numerical solution of Eq. 2a and set μ = [1.49598E4 1.49598E4 1.489 1.489] (this choice of diffusion coefficient is proportional to the grid width and was tuned empirically). To test the limitations of GBEES, we assume an unrealistic, worst-case scenario initial uncertainty in position of 1.5 × 10^4 km and velocity of 1.5 km/s. This measurement uncertainty is much larger than is estimated for the Parker Solar Probe a priori uncertainties [42]. We initialize a MC simulation with the same initial conditions and uncertainty. From the entire orbit of the chosen initial conditions, we propagate uncertainty through a portion of the close approach, as this section of the orbit will be most affected by epistemic uncertainty. No discrete measurement updates are taken in this simulation to demonstrate how the PDF spreads when model uncertainty is non-negligible.
Table 3. Initial Conditions of L_3 Sun-Earth Lyapunov Trajectory
Figure 10. (left) Position and (right) velocity in barycentric, synodic coordinate frame of L_3 Sun-Earth Lyapunov trajectory with initial conditions from Table 3. Period of full Lyapunov orbit is 365.2554 days, and total time of close approach trajectory is 1.8263 days.
As is demonstrated in Figure 11, the PDF propagated by GBEES spreads further in phase space than the distribution approximated by the MC simulation, especially in the velocity space. This behavior is representative of the diffusion term in Eq. (2a) being nonzero. Figure 12 compares the simulation times of GBEES vs. a MC simulation of similar resolution (note that epistemic uncertainty is not considered in the MC simulation).
Figure 11. Initially Gaussian uncertainty time-marched through a close approach of the Sun. To model epistemic uncertainty, the coefficient of diffusion μ = [1.49598E4 1.49598E4 1.489 1.489]. (top) The contour plots represent the PDF propagated by GBEES about the nominal trajectory. Note that the colors of the contours representing the magnitude of the probability of a PDF are not relative to other PDFs, but instead are representative of probability differences throughout the individual distributions. (bottom) GBEES compared with a 500 particle MC simulation, to confirm the accuracy of the method. A complete list of the simulation parameters can be found in Table A1.
Figure 12. Comparing the computational time of GBEES vs. MC when propagating uncertainty for the trajectory from Figure 11.
Conclusion
We have presented a high-dimensional extension to a grid-based, Bayesian estimation algorithm that accurately and efficiently propagates uncertainty for nonlinear systems. The proposed method attempts to address the limitations of the EKF, as well as other commonly-used state estimation methods that propagate uncertainty. To improve the poor time complexity of the legacy implementation, the proposed extension employs a binary search tree that efficiently stores the discretized PDF. Additionally, the algorithm considers the dynamics of the system when handling the grid creation and deletion procedures. The theorized speedup of the extension was validated via application to a chaotic, three-dimensional system, the Lorenz attractor. Having confirmed the improved efficiency of the new algorithm, we developed two realistic scenarios where the dynamics of the PCR3BP (a four-dimensional, nonlinear system where the motion of a negligible mass is governed by two massive bodies) would be utilized for the state estimation of various spacecraft in different regimes of space. These scenarios represent areas where existing uncertainty propagation methods may be lacking, either due to the infrequency of measurement updates or non-negligible epistemic uncertainty. The results of applying the proposed algorithm to these two scenarios are compared with high-resolution MC simulations to validate their accuracy. As the primary limitation of Bayesian estimation is the computational burden, we also compared the program run time of GBEES with the high-resolution MC simulations. To ensure that both techniques were on equal footing for a valid computational efficiency comparison, we set the number of particles in each MC simulation equal to the number of grid cells above the probability threshold value in the final discretized PDF, prior to obtaining a measurement update at Mi.
We acknowledge that there is still work to be done if we are to reasonably argue that the standard for orbit uncertainty propagation be replaced with our proposed Bayesian estimation technique. Said technique was applied to two scenarios where we speculated that the EKF would perform poorly, either due to the infrequency of measurement updates or non-negligible epistemic uncertainty. Having said this, we have yet to conduct a full investigation of the conditions under which the EKF fails. Contemporary mission design is predicated on the measurement duty cycle being frequent enough that the errors of the EKF from linearizing the dynamics are negligible. However, we anticipate that future missions to the Jovian and Saturnian moons will operate spacecraft in chaotic trajectories that require more frequent measurement duty cycles than currently feasible. Confirming this speculation will require a comprehensive review of the EKF applied to the R3BP compared with other estimation techniques (Bayesian estimation, particle filters, GMMs, MC, etc.). This will establish the limitations of the EKF when estimating chaotic trajectories. Additionally, as demonstrated by the computational efficiency comparisons, GBEES fails to reach the speed of the MC simulations for both scenarios. However, we emphasize that for a true MC simulation, the spread of the final distribution will be unknown; in the case stated here, the MC simulations utilized information from the GBEES distribution to provide a reasonable approximated distribution. If the final distribution is unknown, the number of particles may be overestimated to ensure enough are created to form an accurate approximation. That being said, computational efficiency is still the primary limitation of Bayesian estimation.
To address the computational limitation, future work will focus on improving the efficiency of the algorithm in all facets. Concerning computational efficiency, our intention is to parallelize the currently serial algorithm. The probability update step is embarrassingly parallelizable and could easily be sped up by distributing the procedure over a GPU. Concerning the analytical formulation of the problem, time-marching in the Cartesian system is susceptible to inefficiencies, due to the number of fast-changing variables. Rapid, nonlinear changes in state require extremely fine time steps to preserve stability of the scheme. Conversely, orbit elements do not exhibit the normal variability of anomalistic motion as do the coordinates; and these parameters possess a geometric significance clearer than that which can be deduced from the coordinates. By performing calculations in an orbital-element space, uncertainties remain close to linear for longer periods, improving efficiency. Recent work has derived a local action-angle orbit element set for the CR3BP; employing this coordinate system should result in a drastic improvement of computation speed.
Appendix
Listed in Table A1 are the parameters for each GBEES simulation
Table A1. GBEES Simulation Parameters
The code to recreate all of the results provided in this paper can be found at https://github.com/bhanson10/GBEES.
Acknowledgments
This investigation was supported by the NASA Space Technology Graduate Research Opportunity (Grant Number 80NSSC23K1219).
References
[1] R. E. Kalman, “A New Approach to Linear Filtering and Prediction Problems,” Journal of Basic Engineering, Vol. 82, No. 1, 1960, pp. 35–45.
[2] S. F. Schmidt, “The Kalman filter - Its recognition and development for aerospace applications,” Journal of Guidance and Control, Vol. 4, No. 1, 1981, pp. 4–7.
[3] E. Seedhouse, The Engines, pp. 40–57. Cham: Springer International Publishing, 2022.
[4] R. G. Clinton, “Additive Manufacturing and 3D Printing in NASA: An Overview of Current Projects and Future Initiatives for Space Exploration,” Tech. Rep. M15-4218, NASA Marshall Space Flight Center, 2014.
[5] J. Brophy, C. Garner, B. Nakazono, M. Marcucci, M. Henry, and D. Noon, “The ion propulsion system for Dawn,” 39th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, 2003, p. 4542.
[6] Y. Tsuda, O. Mori, R. Funase, H. Sawada, T. Yamamoto, T. Saiki, T. Endo, K. Yonekura, H. Hoshino, and J. Kawaguchi, “Achievement of IKAROS—Japanese deep space solar sail demonstration mission,” Acta Astronautica, Vol. 82, No. 2, 2013, pp. 183–188.
[7] J. Sullivan and C. D’Souza, “Extended Kalman Filter Performance on the Artemis-1 Mission,” 45th Rocky Mountain AAS GN&C Conference, No. AAS 23-053, 2023.
[8] J. Barrow-Green, Poincar ́e and the Three Body Problem. American Mathematical Society/London Mathematical Society, 1997.
[9] S. Ross, W. Koon, M. Lo, and J. Marsden, Dynamical Systems, the Three-Body Problem, and Space Mission Design. 2022.
[10] J. S. Parker, Low-energy ballistic lunar transfers. PhD thesis, University of Colorado at Boulder, 2007.
[11] M. R. Thompson, E. Kayser, J. S. Parker, C. Ott, T. Gardner, and B. Cheetham, “Navigation design of the CAPSTONE mission near NRHO insertion,” AAS/AIAA Astrodynamics Specialist Conference, 2021.
[12] S. Campagnola, A. Boutonnet, W. Martens, and A. Masters, “Mission Design for the Exploration of Neptune and Triton,” Vol. 30, 01 2013.
[13] Y. z. Luo and Z. Yang, “A review of uncertainty propagation in orbital mechanics,” Progress in Aerospace Sciences, Vol. 89, 2017, pp. 23–39.
[14] “NASA Technology Roadmaps TA 11: Modeling, Simulation, Information Technology, and Processing,” tech. rep., National Aeronautics and Space Administration, 2015.
[15] K. Fujimoto, D. J. Scheeres, and K. T. Alfriend, “Analytical nonlinear propagation of uncertainty in the two-body problem,” Journal of Guidance, Control, and Dynamics, Vol. 35, No. 2, 2012, pp. 497–509.
[16] T. R. Bewley and A. S. Sharma, “Efficient grid-based Bayesian estimation of nonlinear low-dimensional systems with sparse non-Gaussian PDFs,” Automatica, Vol. 48, No. 7, 2012, pp. 1286–1290.
[17] M. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking,” IEEE Transactions on Signal Processing, Vol. 50, No. 2, 2002, pp. 174–188.
[18] A. Jazwinski, Stochastic processes and filtering theory. No. 64 in Mathematics in science and engineering, New York, NY [u.a.]: Acad. Press, 1970.
[19] A. D. Fokker, “Die mittlere Energie rotierender elektrischer Dipole im Strahlungsfeld,” Ann. Phys., Vol. 348, 1914, pp. 810–820.
[20] M. Planck, ̈Uber einen Satz der statistischen Dynamik und seine Erweiterung in der Quantentheorie. Sitzungsberichte der K ̈oniglich-Preussischen Akademie der Wissenschaften zu Berlin, Reimer, 1917.
[21] T. Bayes, “An essay towards solving a problem in the doctrine of chances,” Phil. Trans. of the Royal Soc. of London, Vol. 53, 1763, pp. 370–418.
[22] R. J. Leveque, Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems. SIAM, 2007.
[23] E. N. Lorenz, “Deterministic nonperiodic flow,” Journal of atmospheric sciences, Vol. 20, No. 2, 1963, pp. 130–141.
[24] P. F. Windley, “Trees, Forests and Rearranging,” The Computer Journal, Vol. 3, No. 2, 1960, pp. 84–88.
[25] M. P. Szudzik, “The Rosenberg-Strong Pairing Function,” arXiv:1706.04129, 2019.
[26] R. Courant, K. Friedrichs, and H. Lewy, “On the partial difference equations of mathematical physics,” IBM journal of Research and Development, Vol. 11, No. 2, 1967, pp. 215–234.
[27] D. L. Richardson, “Halo Orbit Formulation for the ISEE-3 Mission,” Journal of Guidance and Control, Vol. 3, No. 6, 1980, pp. 543–548.
[28] V. Domingo, B. Fleck, and A. I. Poland, “SOHO: The Solar and Heliospheric Observatory,” Space Science Reviews, Vol. 72, 1995.
[29] V. Szebehely, Theory of Orbits: The Restricted Problem of Three Bodies. New York, NY: Academic Press, 1967.
[30] R. W. Farquhar and A. A. Kamel, “Quasi-periodic orbits about the translunar libration point,” Celestial mechanics, Vol. 7, No. 4, 1973, pp. 458–473.
[31] R. Park, “Jet Propulsion Laboratory Three-Body Periodic Orbit Catalog,” 2023.
[32] S. B. Broschart, N. Bradley, and S. Bhaskaran, “Optical-Based Kinematic Positioning for Deep-Space Navigation,” 2017.
[33] S. B. Broschart, N. Bradley, and S. Bhaskaran, “Kinematic approximation of position accuracy achieved using optical observations of distant asteroids,” Journal of Spacecraft and Rockets, Vol. 56, No. 5, 2019, pp. 1383–1392.
[34] S. Bhaskaran, S. Desai, P. Dumont, B. Kennedy, G. Null, W. Owen Jr, J. Riedel, S. Synnott, and R. Werner, “Orbit determination performance evaluation of the Deep Space 1 autonomous navigation system,” 1998.
[35] J. Riedel, S. Bhaskaran, S. Desai, D. Han, B. Kennedy, T. McElrath, G. Null, M. Ryne, S. Synnott, T. Wang, et al., “Using autonomous navigation for interplanetary missions: The validation of Deep Space 1 AutoNav,” 2000.
[36] S. Bhaskaran, “Autonomous navigation for deep space missions,” SpaceOps 2012, p. 1267135, 2012.
[37] P. Grover and Y. Sato, “Efficient estimation and uncertainty quantification in space mission dynamics,” AIAA/AAS Astrodynamics Specialist Conference, 2012, p. 5062.
[38] W. F. Oberle, Monte Carlo simulations: number of iterations and accuracy. US Army Research Laboratory Aberdeen Proving Ground, MD, 2015.
[39] V. Elvira, J. M ́ıguez, and P. M. Djuri ́c, “Adapting the number of particles in sequential Monte Carlo methods through an online scheme for convergence assessment,” IEEE Transactions on Signal Processing, Vol. 65, No. 7, 2016, pp. 1781–1794.
[40] N. B. Stastny and D. K. Geller, “Autonomous optical navigation at Jupiter: a linear covariance analysis,” Journal of Spacecraft and Rockets, Vol. 45, No. 2, 2008, pp. 290–298.
[41] J. C. Kasper, K. G. Klein, E. Lichko, J. Huang, C. H. K. Chen, S. T. Badman, J. Bonnell, P. L. Whittlesey, R. Livi, D. Larson, M. Pulupa, A. Rahmati, D. Stansby, K. E. Korreck, M. Stevens, A. W. Case, S. D. Bale, M. Maksimovic, M. Moncuquet, K. Goetz, J. S. Halekas, D. Malaspina, N. E. Raouafi, A. Szabo, R. MacDowall, M. Velli, T. D. d. Wit, and G. P. Zank, “Parker Solar Probe Enters the Magnetically Dominated Solar Corona,” Phys. Rev. Lett., Vol. 127, 2021.
[42] D. R. Jones, P. Thompson, P. Valerino, E. Lau, T. Goodson, M.-K. Chung, and N. Mottinger, “Orbit
Determination Covariance Analyses for the Parker Solar Probe Mission,” 2017.