Volunteer Computing for Global Stability Census of One Million Cislunar Orbits and a Deep Learning Surrogate for Rapid Stability Prediction

1. Introduction

The cislunar realm, the gravitational corridor encompassing Earth and the Moon, has re-emerged as a critical domain for future space exploration, satellite deployment, and deep space gateway operations. Its dynamics are primarily governed by the Circular Restricted Three-Body Problem (CR3BP), which gives rise to a complex tapestry of orbital families, interwoven with chaotic manifolds and instability regions. Navigating this environment requires a precise, global understanding of orbital stability. [1,2]

Traditional approaches rely on numerical integration of equations of motion, which, while accurate, are computationally prohibitive when conducting the vast parameter studies necessary for a comprehensive stability census. This creates a critical bottleneck: due to the sheer volume of calculations required, high-fidelity mapping of the cislunar space at the scale needed for systematic mission design remains an elusive goal, even with centralized supercomputing resources.

Meanwhile, the volunteer computing paradigm has matured over the past two decades, proving its capability to tackle problems of unprecedented scale. Pioneered by projects such as SETI@home, which harnessed idle personal computers to analyze radio telescope data in the search for extraterrestrial intelligence, the model demonstrated the potential of aggregating global, heterogeneous computing power. [3] This was extended to astrophysics by projects such as Einstein@Home, which uses volunteer computing to search for gravitational waves and pulsars, leading to numerous peer-reviewed discoveries. [4]

These successes demonstrate the capability of volunteer computing for high-throughput computing across diverse domains. However, within the field of classical orbital mechanics—a discipline with inherently parallelizable, high-throughput computational needs perfectly suited to this distributed model—volunteer computing campaigns focused on cislunar orbital stability analysis have been notably absent.

This work addresses this gap by applying volunteer computing to a foundational problem in astrodynamics. We use the fully public “Million Orbit” data set—a product of a massive supercomputing campaign by Lawrence Livermore National Laboratory [5,6]—as the input for a novel, distributed stability analysis. While the data set provides the raw trajectories, our project repurposes it for global stability analysis. We developed a containerized application running on the BOINC platform to calculate the Jacobi constant time series for all sampled points along each trajectory in the data set, performing a stability census at a scale and cost-effectiveness impractical for traditional centralized resources.

Beyond producing a static stability map, we recognize that real-time mission planning and large-scale orbit screening require even faster inference than recomputing the full Jacobi time series. Therefore, we further develop a machine learning surrogate that predicts whether an orbit belongs to the stable Region I based solely on the first K Jacobi constant values (prefix). Our model combines an LSTM encoder with an attention mechanism to extract compact features, followed by an XGBoost classifier. We evaluate prefix lengths $K=10,20,50,100,200,500$ and compare against a pure XGBoost baseline using the same train/validation/test splits. The results show that the LSTM+Attention+XGBoost model achieves near-perfect AUC with longer prefixes and remains highly informative even with very short prefixes, demonstrating its utility as a rapid screening tool.

2. Methods

2.1. The “Million Orbit” Reference Data Set

The computational foundation of this census is the fully public “One Million Open-source Cislunar Orbits” data set, generated through a massive supercomputing campaign by Lawrence Livermore National Laboratory (LLNL). [5,6] This data set provides a basis for the study of cislunar dynamics within the Circular Restricted Three-Body Problem (CR3BP) framework.

The data set comprises one million unique spacecraft trajectories, each numerically integrated over a simulated period of six years. The orbits are initialized from geosynchronous Earth orbit to regions beyond the Moon, encompassing the dynamically complex regions between the Earth and the Moon, including Lagrange points and their associated manifolds. Each orbit is stored in a standardized Hierarchical Data Format version 5 (HDF5) file, [7] containing time-series arrays of the spacecraft’s three-dimensional position and velocity vectors in the Earth-Moon rotating frame. [5,6]

For our distributed computing paradigm, this data set presents ideal characteristics:

1.

Web access: All data files are hosted on a public LLNL repository with unrestricted access.

2.

Inherent modularity: The million orbits are stored in independently accessible files and logical groups, allowing for trivial decomposition into discrete computational units.

3.

High scientific value: As a product of a verified high-performance computing (HPC) integration, it provides a benchmark-quality source for secondary analysis, distinguishing our stability metrics from those derived from ad-hoc simulations.

In this study, we treat this data set not as the subject of new numerical integration, but as the primary input for a subsequent, large-scale analytical computation—the derivation of Jacobi constant time series and stability indices for every data point across all trajectories.

2.2. Task Definition: Aligning BOINC Work Units with Existing Data Structure

The architecture of the source “Million Orbit” data set is modular, with trajectories divided into 20,000 HDF5 files, each containing 50 independent orbital groups. This structure enabled a direct and efficient mapping to the BOINC model [8]: each job (or `work unit’) handled one complete HDF5 file, thereby eliminating the need for further data segmentation or pre-processing.

However, because the total size of all HDF5 files amounts to terabytes, we did not embed them in work units. Instead, each work unit contains only a JavaScript Object Notation (JSON) data interchange format file [9] (input.json), which contains two key parameters:

• Data Source Parameters, including the URL of the target .h5 file on the LLNL server [10] and the specific file identifier.
• Task ID, which is used to identify tasks and to enable debugging.

The task of each work unit was to: (i) fetch the assigned HDF5 file, (ii) calculate the time series of the Jacobi constant for all 50 trajectories within it, and (iii) return the derived data set (also in JSON) along with integrity-validation metadata in a result file.

This “one-file-per-task” granularity—corresponding to 50 orbits per computation—was chosen to balance multiple factors. It ensured that individual tasks required several hours of computation on typical volunteer hardware, optimizing the use of donated resources while remaining within practical completion time windows. This granularity also made the management overhead of 20,000 total tasks tractable and provided fault tolerance. The failure of any single task affected only 0.005% of the total data set, allowing for swift re-computation with negligible impact on the overall campaign. To ensure result accuracy, we used redundant computing: we randomly run some of the jobs twice and check if the results agree.

2.3. Design and Implementation of the BOINC Application

To execute the stability census on volunteer resources, we developed a custom application using BOINC’s Universal Docker Application (BUDA) framework. The application was designed with three primary goals: (1) integration with the existing BOINC infrastructure, (2) fault tolerance to handle the heterogeneous and intermittent nature of volunteer nodes, and (3) reproducibility of the computational environment.

The code for fetching the HDF5 files, parsing orbit data, and computing the Jacobi constant time series (as described in the jacobi_calculator.py module), was written in Python. To ensure a consistent runtime environment across the diverse range of volunteer operating systems (Windows, macOS, Linux), we encapsulated the application, along with all its dependencies (NumPy, SciPy, h5py, etc.), in a Docker container. The Dockerfile defines an image based on python:3.9-bullseye, which includes all the required libraries. This containerization strategy guarantees that each work unit executes in an identical software environment, eliminating the “it works on my machine” problem inherent in distributed computing. [11]

Each BOINC work unit corresponds to a single JSON task file (input.json), as detailed in the Task Definition section. The application workflow, orchestrated by buda_main.py, is as follows:

• Initialization: The application sets up the configurations and reads input.json to obtain the target HDF5 file’s URL and other parameters.
• Remote Data Fetching: It downloads the required HDF5 file from the public LLNL repository.
• Checkpointed Computation: It iterates through each of the 50 orbit groups within the downloaded file. A useful feature for volunteer computing is a checkpointing mechanism. Every 300 seconds, the application saves its intermediate state—including the computed results for that orbit and the index of the last completed group—to a checkpoint file. If a volunteer’s computer is interrupted and later resumes the same work unit, the application detects the existing checkpoint file and resumes processing from where it left off, rather than starting over. This prevents the loss of computational progress and ensures efficient use of donated cycles.
• Output and Cleanup: Upon successful processing of all 50 orbit groups, the application generates the final output.json file containing all results and metadata. If the process encounters a non-recoverable error, it logs the error to the output file and exits gracefully.

2.4. Stability Classification from Jacobi Constants

From the computed Jacobi constant time series for each orbit, we classified the orbit into one of five classical stability regions based on the energy thresholds $C_1,C_2,C_3,C_4$ for the Earth-Moon CR3BP with mass ratio $\mu =0.01215$. These thresholds are:

$$C_1\approx 3.18834,C_2\approx 3.17216,C_3\approx 3.02413,C_4\approx 2.98796.$$

For each time step, the instantaneous Jacobi constant C is compared to these thresholds. An orbit is assigned to the region that appears most frequently over its entire time series. Additionally, the maximum Jacobi constant was recorded to identify high-energy events.

2.5. Machine Learning Surrogate: LSTM+Attention+XGBoost

To enable rapid prediction of Region I stability from only the first K Jacobi constant values, we developed a hybrid model combining a sequential encoder and a gradient-boosted classifier.

2.5.1. Data Preparation

From the full set of one million orbits, we extracted the first K Jacobi constants as the feature vector $X\in {\mathbb{R}}^K$, and the binary label $y=1$ if the orbit belongs to Region I, else 0. We used the same train/validation/test split strategy across all K to ensure comparability:

• Test set: 15% of orbits, stratified by label.
• Validation set: 15% of the remaining 85% (i.e., 12.75% of total), also stratified.
• Training set: the rest (72.25% of total).

The splits are deterministic (random seed 42) and are shared between the LSTM+XGBoost model and the pure XGBoost baseline.

2.5.2. LSTM Encoder with Attention

The LSTM encoder transforms the input sequence (length K, 1 feature per time step) into a fixed-dimensional context vector. Our architecture:

• Input shape: $(batch,K,1)$.
• LSTM: 2 layers, hidden size 64, batch-first.
• Attention: A linear layer maps each LSTM hidden state to a scalar score; softmax over time produces attention weights ${\alpha }_t$. The context vector is the weighted sum of hidden states: $\mathbf{c}={\sum }_t{\alpha }_t{\mathbf{h}}_t$.

The attention mechanism allows the model to focus on the most informative time steps for stability classification.

The LSTM is trained jointly with a linear classification head (cross-entropy loss) on the validation set, with early stopping (patience=5) and class weights to handle imbalance. The encoder weights are saved at the best validation loss.

2.5.3. XGBoost on LSTM Features

After training the LSTM encoder, we discard the classification head and use the encoder as a feature extractor. For each orbit, we compute the 64-dimensional context vector $\mathbf{c}$. We then train an XGBoost classifier on these features using the training set, with hyperparameters:

• max_depth = 6, learning_rate = 0.1, subsample = 0.8, colsample_bytree = 0.8.
• Early stopping on validation AUC (5-fold cross-validation on training set) to select the optimal number of boosting rounds.
• scale_pos_weight automatically set to balance positive/negative class.

The final model is evaluated on the held-out test set.

2.5.4. Baseline: Raw Prefix + XGBoost

For comparison, we also train a pure XGBoost model directly on the raw prefix vectors $X\in {\mathbb{R}}^K$ using the same train/validation/test splits and identical XGBoost hyperparameters (except no feature extraction). This baseline quantifies the benefit of the LSTM+Attention encoding.

3. Results

3.1. Volunteer Computing Campaign Overview

The volunteer computing campaign successfully processed all 20,000 work units, covering the complete set of one million orbits. A total of approximately 16 billion individual Jacobi constant values were computed, with an average of ∼16,000 points per orbit.

The software and result dataset has been published to Zenodo. [11,12]

3.2. Stability Region Distribution

Based on the classical CR3BP energy thresholds, the distribution of the one million orbits across the five stability regions is shown in Table 1.

Notably, Region II is completely absent, and Region IV is represented by a single orbit (ID 754482). The Jacobi constant time series for this unique orbit is shown in Figure 1. Its behavior exhibits a distinctive pattern that may indicate a transitional state between stability and escape.

A total of 173 orbits were found to have maximum Jacobi constant exceeding 100. All of these belong to Region V and correspond to high-energy (possibly escape) trajectories. These values have been verified independently.

3.3. Predicting Region I from Prefix Sequences [16,17]

We evaluated both the proposed LSTM+Attention+XGBoost model and the raw XGBoost baseline for prefix lengths $K=10,20,50,100,200,500$. Table 2 reports test set metrics for the LSTM+XGBoost model, and Table 3 for the raw XGBoost baseline.

The LSTM+Attention+XGBoost model consistently outperforms the raw XGBoost baseline for all $K<500$, especially for small K (e.g., at $K=10$: AUC 0.9286 vs. 0.9109). At $K=500$ both models achieve nearly identical AUC (0.9840 vs. 0.9849), but the LSTM-based model uses a much lower dimensional feature space (64 instead of 500) and provides interpretable attention weights.

Figure 2 plots AUC as a function of K for LSTM+Attention+XGBoost models while Figure 3 plots AUC as a function of K for XGBoost models, showing that the LSTM+Attention model reaches high performance more rapidly with increasing prefix length.

3.4. Attention Interpretation

The attention mechanism assigns weights to each time step in the input prefix. For the test set, we visualized attention weights for randomly selected stable and unstable orbits. Figure 4 shows example attention patterns for $K=500$. Stable orbits often exhibit attention concentrated on early or mid-sequence segments, while unstable orbits show more diffuse or late-focused attention. This suggests that the model learns to identify dynamical signatures indicative of long-term stability from only the initial Jacobi constant evolution.

3.5. Feature Importance of LSTM Encoded Features

The XGBoost model trained on LSTM features (64 dimensions) produces feature importance scores (by weight). The top 20 most important features are shown in Figure 5. No single feature dominates, indicating that the LSTM encoder distributes information across multiple dimensions.

4. Discussion

4.1. Stability Census Findings

The overwhelming majority of orbits (91.68%) fall into Region V, the highest energy region, indicating that most trajectories in the LLNL dataset are either distant retrograde orbits or escape orbits. This is consistent with the dataset’s design, which spans from GEO to beyond the Moon.

Region I, corresponding to the most stable low-energy orbits (e.g., halo, Lyapunov families), comprises 8.07% of the dataset. These orbits are of particular interest for long-term missions such as lunar gateways.

The complete absence of Region II orbits suggests a dynamic gap: no trajectories occupy the energy range between the stable low-energy region and the intermediate Region III. This may reflect a natural separation in the phase space of the Earth-Moon system.

The single Region IV orbit (ID 754482) is a rare find. Its Jacobi constant hovers just above $C_4$ for most of the six-year integration, occasionally dipping below. This behavior could indicate a near-escape trajectory that remains bound for a long duration, possibly due to weak chaos or resonance effects. Further analysis of its Lyapunov exponent would be valuable.

The 173 high-energy orbits ($C>100$) likely represent numerical artifacts or extremely fast escapes. Manual verification confirmed they are not the result of software errors.

4.2. Machine Learning Surrogate Performance

The LSTM+Attention+XGBoost model demonstrates that stable Region I orbits can be reliably identified from surprisingly short prefixes. Even with only the first 10 Jacobi constants, the AUC reaches 0.9286, and with $K=20$ it exceeds 0.968. This means that, in practice, one only needs to integrate a trajectory for a few hundred time steps (or even fewer) to decide whether it will remain in the stable region for the entire six-year simulation.

The raw XGBoost baseline, while strong at large K, suffers more from the curse of dimensionality when K is small. The LSTM encoder effectively compresses the prefix into a fixed-size representation, regularizing the model and improving generalization. The attention mechanism further provides interpretability, revealing which parts of the early trajectory are most predictive of long-term stability.

One limitation of the current model is that it only predicts membership in Region I (the stable class). In practice, mission designers might also want to distinguish between, e.g., Region III (intermediate) and Region V (escape). Extending the model to multi-class prediction is straightforward and left for future work.

4.3. Volunteer Computing as an Enabler

The volunteer computing approach proved highly effective: 20,000 tasks were completed in approximately two months, at zero financial cost. This demonstrates the potential of BOINC for large-scale astrodynamics studies. The resulting dataset (16 billion Jacobi constants) is freely available and can serve as a benchmark for machine learning models in celestial mechanics.

4.4. Future Directions

Several extensions are promising:

• Multi-class prediction for all five stability regions.
• Using raw position/velocity time series instead of Jacobi constants.
• Recurrent or transformer-based models that can handle variable-length inputs.
• Deploying the trained model as a lightweight service for real-time orbit screening.

5. Conclusion

We have presented the first volunteer computing campaign dedicated to cislunar orbital stability analysis, processing one million orbits and producing over 16 billion Jacobi constant values. The resulting stability map quantifies the prevalence of Region V (high-energy) and Region I (stable) orbits, and identifies a single rare Region IV orbit. Furthermore, we developed a deep learning surrogate (LSTM+Attention+XGBoost) that predicts Region I membership from only the first K Jacobi constants, achieving high accuracy even with $K=10$. This model is orders of magnitude faster than full numerical integration and can be used for rapid screening in mission design. Together, our work establishes volunteer computing as a transformative tool for astrodynamics and provides a practical machine learning solution for real-time stability assessment.