# Reduce API Debugging Overhead via Knowledge Prepositioning

Shujun Wang<EMAIL_ADDRESS>0000-0001-8161-1695 Alibaba
GroupBeijingChina100000 ,  Yongqiang Tian<EMAIL_ADDRESS>Alibaba
GroupHangZhouChina100000 and  Dengcheng He<EMAIL_ADDRESS>Alibaba GroupHangZhouChina100000

(2023)

###### Abstract.

OpenAPI indicates a behavior where producers offer Application Programming
Interfaces (APIs) to help end-users access their data, resources, and
services. Generally, API has many parameters that need to be entered. However,
it is challenging for users to understand and document these parameters
correctly. This paper develops an API workbench to help users learn and debug
APIs. Based on this workbench, much exploratory work has been proposed to
reduce the overhead of learning and debugging APIs. We explore the knowledge,
such as parameter characteristics (e.g., enumerability) and constraints (e.g.,
maximum/minimum value), from the massive API call logs to narrow the range of
parameter values. Then, we propose a fine-grained approach to enrich the API
documentation by extracting dependency knowledge between APIs. Finally, we
present a learning-based prediction method to predict API execution results
before the API is called, significantly reducing user debugging cycles. The
experiments evaluated on the online system show that this work’s approach
substantially improves the user experience of debugging OpenAPIs.

API, Debugging, Constraints, Knowledge, Prediction

††journalyear: 2023††copyright: acmlicensed††conference: Companion Proceedings
of the ACM Web Conference 2023; April 30-May 4, 2023; Austin, TX,
USA††booktitle: Companion Proceedings of the ACM Web Conference 2023 (WWW ’23
Companion), April 30-May 4, 2023, Austin, TX, USA††price: 15.00††doi:
10.1145/3543873.3587311††isbn: 978-1-4503-9419-2/23/04††ccs: Networks Cloud
computing

##  1\. Introduction

Application Programming Interfaces (APIs) play an essential role in modern
software development(Huang et al., 2018; Lamothe, 2020; Zhong et al., 2020;
Shen et al., 2020). With the help of APIs, developers can complete their tasks
more efficiently(Gu et al., 2016). OpenAPI indicates a behavior where
producers offer Application Programming Interfaces (APIs) to help end-users
access their data, resources, and services. Increasingly, many companies have
published an ocean of OpenAPIs on the internet(Murphy et al., 2017). For
example, as of October 2022, Alibaba Cloud lists over 15560 OpenAPIs for web
services . Moreover, more than 60,000 OpenAPIs will be launched soon.

However, ”Learning OpenAPIs” is an empirical challenge for end-users of the
knowledge gap between API designers and API users(Lamothe, 2020). To bridge
this gap, we developed an OpenAPI developer workbench equipped with multiple
auxiliary tools to help users learn and employ the OpenAPIs. For instance,API
documenration(Sohan et al., 2017; Zhong and Mei, 2019), API debugging tool and
SDK generation tool111https://github.com/aliyun/darabonba. Unfortunately,
although we provide detailed OpenAPI usage documentation and tools, the
success rate of OpenAPI calls is still unsatisfactory(Meng et al., 2019).

By analyzing the user behavior sequence of the API workbench, we discovered
frequent user behavior (Read Documentation $\rightarrow$ Debug APIs). Inspired
by the above, in this paper, we address the problem of improving the success
rate of API debugging. We conclude that the success rate of API debugging is
an important indicator that intuitively reflects the ease of API learning.
Crucially, the contributions are outlined as follows:

  * •

We design a workbench (https://next.api.alibabacloud.com) to help users learn
and use APIs, and we are willing to share all the data from this workbench.

  * •

We extract parameter constraints and rules from massive API call logs. Based
on this, we reduced the complexity of some parameters from input to selection,
e.g., enumeration.

  * •

We propose a fine-grained approach to extract dependencies between APIs to
ensure that the API is in an executable state.

  * •

We present a classification-based API call result prediction method for saving
debugging time.

##  2\. Parameter Constraints

OpenAPIs may have constraints on parameters, such that not all parameters are
always required or always optional. The presence of one parameter could cause
another parameter to be required. A clear overview of the constraints helps
API consumers integrate without additional support and with fewer integration
faults. We detect two levels of parameter constraints in our work:

  1. (1)

API Level: Parameter combination problem, i.e., for the same OpenAPI, the
sequence of input parameters may be different since there are optional
parameters.

  2. (2)

Parameter Level: Mandatory parameter, Enumeration parameter extraction, and
Numeric parameter range.

###  2.1. Example Extraction

The examples of the parameter are essential. An intuitive way to obtain
example values of parameters is to extract them from the API call log
directly(i.e., Extract parameter values from a correct API call). However,
random extracting examples from massive values is insufficient. We believe
that the example values of parameters should be universal. For instance, 90%
parameter values are composed of pure English letters, and the examples should
be one of the values of these 90% rather than others.

This subsection proposes a MapReduce-based Parameter abstraction approach to
extract the common features of an ocean of parameter values. 1

1

Input : A set of parameter values $V=v_{1}\cup v_{2}\cup\ldots\cup v_{n}$

Output : A common subsequence $s$

A set of parameter patterns $P=p_{1}\cup\ldots\cup p_{n}$

A set of parameter lengths $L=l_{1}\cup l_{2}\cup\ldots\cup l_{n}$

2

3Mapper

4 $s_{i}\leftarrow$ExtractCommonSubsequence(_$v_{i}$_);

5 $p_{i}\leftarrow$TransformAndCompress(_$v_{i}$_);

6 $l_{i}\leftarrow$LengthComputing(_$v_{i}$_);

7

8Reducer

9 $s\leftarrow$MergeSubsequence(_$s_{i}$_);

10 $P\leftarrow$MergePatterns(_$p_{i}$_);

11 $L\leftarrow$MergeLengths(_$l_{i}$_);

return $s$, $P$, $L$

Algorithm 1 Parameter Abstraction

As shown in Algorithm 1, we take a two-stage approach: Mapper and Reducer, to
extract common parameter features from an ocean of examples. We emphasize
three types of parameter features: Component Element Analysis, Element
Arrangement Analysis, and Length Analysis.

There are two main stages of parameter abstraction:

Mapper: In the mapper stage, we mainly focus on extracting local features of
parameter values. Specifically,

  * •

Extracting Common Subsequence will extract the longest common subsequences in
the parameter set $v_{i}$.

  * •

Transformer And Compress will translate a specific parameter value into an
abstract value. The translation rules are as follows:

Table 1. Transform Rules Representation | Description
---|---
z | The Chinese character
x | English lowercase characters
X | English uppercase characters
d | Number
| Other characters reserved

The transformer stage is responsible for abstract parameters. For example,
$123$ will be abstracted as $ddd$, and the compressing stage will further
translate it into a character $d$.

  * •

LengthComputing mainly computes the length information of parameters and their
frequency.

Reducer: In this stage, the local parameter features are combined to extract
the global characteristics of the parameters.

###  2.2. Parameter Combination

In order to successfully employ an OpenAPI, there may be different parameter
sequences. Take Alibaba Cloud’s OpenAPI SendSms (Send Message) as an example,
which has two parameter sequences (See Table 2). Note that Outid is an
optional parameter.

Table 2. Parameter Sequences OpenAPI | Parameter Sequences
---|---
SendSms | PhoneNumbers, SignName, TemplateCode, TemplateParam
SendSms | PhoneNumbers, SignName, TemplateCode, TemplateParam, OutId
Table 3. Parameter Sequence Rate OpenAPI | Parameter Sequences | Rate
---|---|---
SendSms | PhoneNumbers, SignName, TemplateCode, TemplateParam | 0.7
SendSms | PhoneNumbers, SignName, TemplateCode, TemplateParam, OutId | 0.3

Algorithm 2 outlines the parameter sequence extraction process. We present a
MapReduce approach to gather all parameter sequences and compute the
frequency. We first extract all parameters from an API request (Line 2). Then,
we filter the useless parameters, for example, private or system parameters
(Line 3). After this, we can obtain all valuable parameters. Further, we sort
these parameters to merge them as keys (Line 4). Finally, we collect all keys
(sorted parameter sequence) and compute their frequency (Line 5). In the
reduce stage, we combined the intermediate results (Line 7).

1

Input : API Requests $R=r_{1}\cup r_{2}\cup\ldots\cup r_{n}$

Output : An Object $O$, the key is the parameter sequence and value is the
frequency

2

3Mapper

4 $t_{i}\leftarrow$ExtractParameters(_$r_{i}$_);

5 $p_{i}\leftarrow$FilterParameters(_$t_{i}$_);

6 $k_{i}\leftarrow$SortedParametersAsKey(_$p_{i}$_);

7 $m_{i}\leftarrow$MergeKeyCount(_$k_{i}$_);

8 Reducer

9 $O\leftarrow$MergeKeyCount(_$m_{i}$_);

return $O$

Algorithm 2 Parameter Sequences Extraction

###  2.3. Enumeration Parameter

Enumeration value is an essential parameter type, which reduces the complexity
from the input level to the select level. In the case of sufficient computing
resources, a simple SQL can complete the enumeration value extraction.

⬇

1SELECT COUNT(DISTINCT value) as number,

2DISTINCT value as values

3FROM An ocean of parameter values

4WHERE number < 20

Where 20 is a hyperparameter. In the case of insufficient computing resources,
suppose we have obtained a parameter $p_{i}$ of OpenAPI $a_{i}$ is enumerable,
and the enumeration value set is $v_{i}$. After that, we process a small part
of the data every day and get an enumeration value set of $v_{j}$, and we only
need to guarantee — $v_{i}\cup v_{j}$—¡ 20.

##  3\. API Dependencies Computing

We propose a fine-grained approach to compute the parameter-level OpenAPI
relevance in this subsection. We have noticed that the input parameter of one
OpenAPI is often the output parameter of some other OpenAPIs.

Table 4. SendSms Input | Output
---|---
PhoneNumbers | Code
SignName | Message
TemplateCode | BizId
TemplateParam | RequestId
OutId |

Table 5.  AddSmsSign Input | Output
---|---
SignName | Code
Remark | Message
SignFileList | SignName
TemplateParam | RequestId

Table 5, 5 exhibits the input parameters and output parameters of two OpenAPIs
(i.e., SendSms, AddSmsSign). The output parameter SignName of AddSmsSign is
also the input parameter of SendSms. Regarding query semantics, if the user
wants to use a SignName, it should be created first.

We design a generation and ranking approach to compute this co-occurrence.

  * •

Candidate Generation: For any OpenAPI input parameter $p_{i}$, we output a set
of OpenAPIs $C$, whose output parameter contains $p_{i}$.

  * •

Ranking: Given an OpenAPI $a_{i}$ and an input parameter $p_{i}$ of $a_{i}$,
the task of this stage is to rank the OpenAPI in $C$, and then recommend a
subset of $C$ to end-users. Given an OpenAPI $c$ of $C$, the score of
$s(a_{i},p_{i},c_{i})$ as follows:

(1) $s(a_{i},p_{i},c)=\frac{\alpha sim(a_{i},c)\cdot\beta
sim(p_{i},c)\cdot\sigma s(a_{i},c)}{c_{i}+c_{o}}$

Where $\alpha,\beta,\sigma$ is the hyperparameter, $sim(x,y)$ means the string
similarity between $x$ and $y$, and $s(x,y)$ is the OpenAPI relevance between
$x$ and $y$.

##  4\. Char-Level Error Prediction

In this subsection, we mainly discuss how to predict the result of an API
request before the request is executed. We reduce the result prediction task
to a multi-classification job to solve this problem. The categories are
various Error Codes and the $\langle$right$\rangle$ label, representing the
correct execution.

An example of SendMessage as follows:

⬇

1{

2 "PhoneNumbers":"186xxx9602",

3 "SignName":"hanxing",

4 "TemplateCode":"SMS_209470795"

5}

The parameter value does not have actual meaning, such as the value of
TemplateCode, i.e., SMS_209470795. This paper explores treating text as a raw
signal at the character level and applying temporal (one-dimensional) ConvNets
to solve it.

###  4.1. Key Modules

The main component is the temporal convolutional module, which computes a 1-D
convolution. Suppose we have a discrete input function
$g(x)\in[1,l]\rightarrow\mathbb{R}$ and a discrete kernel function
$f(x)\in[1,k]\rightarrow\mathbb{R}$ The convolution
$h(y)\in[1,\lfloor(l-\dot{k})/d\rfloor+1]\rightarrow\mathbb{R}$ between $f(x)$
and $g(x)$ with stride $d$ is defined as

(2) $h(y)=\sum_{x=1}^{k}f(x)\cdot g(y\cdot d-x+c)$

Where $c=k-d+1$ is an offset constant, just as in traditional convolutional
networks in vision, the module is parameterized by a set of such kernel
functions $f_{ij}(x)(i=1,2,\ldots,m$ and $j=$ $1,2,\ldots,n)$ which we call
weights, on a set of inputs $g_{i}(x)$ and outputs $h_{j}(y).$ We call each
$g_{i}$ (or $h_{j}$ ) input (or output) features, and $m$ (or $n$ ) input (or
output) feature size. The outputs $h_{j}(y)$ is obtained by a sum over $i$ of
the convolutions between $g_{i}(x)$ and $f_{ij}(x)$

One key module that helped us to train deeper models is temporal max-pooling.
It is the 1 -D version of the max-pooling module used in computer vision [2].
Given a discrete input function $g(x)\in$ $[1,l]\rightarrow\mathbb{R}$, the
max-pooling function $h(y)\in[1,\lfloor(l-k)/d\rfloor+1]\rightarrow\mathbb{R}$
of $g(x)$ is defined as

$h(y)=\max_{x=1}^{k}g(y\cdot d-x+c)$

Where $c=k-d+1$ is an offset constant, this pooling module enabled us to train
ConvNets deeper than six layers, where all others fail. The non-linearity used
in our model is the rectifier or thresholding function $h(x)=\max\\{0,x\\}$,
which makes our convolutional layers similar to rectified linear units
(ReLUs). The algorithm used is stochastic gradient descent (SGD) with a
minibatch of size 128, using momentum $0.9$ and initial step size $0.01$,
which is halved every three epochs for ten times.

###  4.2. Character quantization

Our models accept a sequence of encoded characters as input. The encoding is
done by prescribing an alphabet of size $m$ for the input language and then
quantizing each character using 1 -of- $m$ encoding (or ”one-hot” encoding).
Then, the sequence of characters is transformed to a sequence of such $m$
sized vectors with fixed length $l_{0}$. Any character exceeding length
$l_{0}$ is ignored, and any characters not in the alphabet, including blank
characters, are quantized as all-zero vectors. The character quantization
order is backward. The latest reading on characters is always placed near the
beginning of the output, making it easy for fully connected layers to
associate weights with the latest reading. The alphabet used in our models
consists of 96 characters, including 52 english letters (Case Sensitive), 10
digits, 33 other characters.

⬇

1 abcdefghijklmnopqrstuvwxyz

2 ABCDEFGHIJKLMNOPQRSTUVWXYZ

3 0123456789

4 -,;/!?:/\|_@#$%^&*~‘+-=<>(){}

##  5\. Evaluation

In this section, we briefly introduce the metrics used in subsequent
experiments. $sr$ represents the success rate of OpenAPI calls.

(3) $sr=\dfrac{|Call~{}Success|}{|Call~{}Number|}$

###  5.1. Overall Evaluation

Table 6. Success Rate Comparision Model | Period 1 | Period 2 | Period 3 | Period 4 | Period 5
---|---|---|---|---|---
Baseline | $23.2\%$ | $24.5\%$ | $25.3\%$ | ${\char 35\relax}$ | ${\char 35\relax}$
Base-DE | ${\char 35\relax}$ | ${\char 35\relax}$ | ${\char 35\relax}$ | $49.7\%$ | $53.2\%$

The core goal of this work is to reduce the cost of users learning API and
improve the API success rate.

Table 6 shows that the proposed optimization strategies (w.r.t Debugging
Experience Optimization and Documentation experience Optimization) achieved
better $sr$ scores than pure baseline. The Baseline model refers to the
workbench that does not add optimization solutions. Base-DE means that it
increased the API debugging strategy. One period means one week. Our
experiments are based on real statistics from online platforms. In Period 4,
we launched the API parameter constraints and rules optimization strategy, and
in Period 5, we launched the prediction optimization strategy.

###  5.2. Prediction Evaluation

We designed 2 ConvNets - one large and one small. They are 9 layers deep with
6 convolutional layers and 3 fully-connected layers. The input has several
features equal to 96 due to our character quantization method, and the input
feature length is 1014 . It seems that 1014 characters could already capture
most of the texts of interest. We also insert 2 dropout modules in between the
3 fully-connected layers to regularize, and the dropout probability is $0.5.$

We initialize the weights using a Gaussian distribution. The mean and standard
deviation for initializing the large model is $(0,0.02)$ and the small model
$(0,0.05)$.

Table 7. Fully-connected layers used in our experiments. The problem determines the number of output units for the last layer. For example, for a 10 -class classification problem, it will be $10.$ Layer | Output Units Large | Output Units Small
---|---|---
7 | 2048 | 1024
8 | 2048 | 1024
9 | Depends on the problem

For different problems, the input lengths may be different (for example, in
our case, $l_{0}=1014$ ), and so are the frame lengths. From our model design,
it is easy to know that given input length $l_{0}$, the output frame length
after the last convolutional layer (but before any of the fully-connected
layers) is $l_{6}=\left(l_{0}-96\right)/27$. This number multiplied with the
frame size at layer 6 will give the input dimension the first fully-connected
layer accepts.

Table 8. Precision of all the models Model | Precision
---|---
Lg. Full Conv | $80.61\%$
Sm. Full Conv | $85.36\%$

Table 8 shows the reliable results of our char-level models.

## CONCLUSIONS

APIs have been widely used but are often poorly learned since many parameters
need to be entered. A direct proof is that the success rate of API debugging
is extremely low (¡37%). This phenomenon causes the loss of a large number of
potential users. To reduce the overhead of learning APIs, we developed an API
workbench to help users understand and debug APIs; based on this workbench, we
collected massive user debug logs and conducted exploratory research. The
practice has shown that parameter limitation, rule mining, and parameter range
limitation are practical solutions to reduce the cost of learning APIs for
users. In addition, helping users know the API’s possible execution results in
advance is a remarkably effective solution for improving the API debugging
experience.

## Acknowledgments

This work was supported by Alibaba Group through Alibaba Innovative Research
Program.

## References

  * (1)
  * Gu et al. (2016) Xiaodong Gu, Hongyu Zhang, Dongmei Zhang, and Sunghun Kim. 2016.  Deep API learning. In _Proceedings of the 24th ACM SIGSOFT International Symposium on Foundations of Software Engineering, FSE 2016, Seattle, WA, USA, November 13-18, 2016_ , Thomas Zimmermann, Jane Cleland-Huang, and Zhendong Su (Eds.). ACM, 631–642.  https://doi.org/10.1145/2950290.2950334
  * Huang et al. (2018) Qiao Huang, Xin Xia, Zhenchang Xing, David Lo, and Xinyu Wang. 2018.  API method recommendation without worrying about the task-API knowledge gap. In _2018 33rd IEEE/ACM International Conference on Automated Software Engineering (ASE)_. IEEE, 293–304.
  * Lamothe (2020) Maxime Lamothe. 2020\.  Bridging the divide between API users and API developers by mining public code repositories. In _ICSE ’20: 42nd International Conference on Software Engineering, Companion Volume, Seoul, South Korea, 27 June - 19 July, 2020_ , Gregg Rothermel and Doo-Hwan Bae (Eds.). ACM, 178–181.  https://doi.org/10.1145/3377812.3382124
  * Meng et al. (2019) Michael Meng, Stephanie Steinhardt, and Andreas Schubert. 2019\.  How developers use API documentation: an observation study.  _Communication Design Quarterly Review_ 7, 2 (2019), 40–49.
  * Murphy et al. (2017) Lauren Murphy, Tosin Alliyu, Andrew Macvean, Mary Beth Kery, and Brad A Myers. 2017.  Preliminary Analysis of REST API Style Guidelines.  _Ann Arbor_ 1001 (2017), 48109.
  * Shen et al. (2020) Qi Shen, Shijun Wu, Yanzhen Zou, Zixiao Zhu, and Bing Xie. 2020.  From API to NLI: A new interface for library reuse.  _J. Syst. Softw._ 169 (2020), 110728.  https://doi.org/10.1016/j.jss.2020.110728
  * Sohan et al. (2017) S. M. Sohan, Craig Anslow, and Frank Maurer. 2017\.  Automated Example Oriented REST API Documentation at Cisco. In _39th IEEE/ACM International Conference on Software Engineering: Software Engineering in Practice Track, ICSE-SEIP 2017, Buenos Aires, Argentina, May 20-28, 2017_. IEEE Computer Society, 213–222.  https://doi.org/10.1109/ICSE-SEIP.2017.11
  * Zhong and Mei (2019) Hao Zhong and Hong Mei. 2019.  An Empirical Study on API Usages.  _IEEE Trans. Software Eng._ 45, 4 (2019), 319–334.  https://doi.org/10.1109/TSE.2017.2782280
  * Zhong et al. (2020) Hao Zhong, Na Meng, Zexuan Li, and Li Jia. 2020\.  An empirical study on API parameter rules. In _ICSE ’20: 42nd International Conference on Software Engineering, Seoul, South Korea, 27 June - 19 July, 2020_ , Gregg Rothermel and Doo-Hwan Bae (Eds.). ACM, 899–911.  https://doi.org/10.1145/3377811.3380922


# Towards a Machine Learning-Based Approach to Predict Space Object Density
Distributions

Victor Rodriguez-Fernandez111Associate Professor, Department of Computer
Systems Engineering<EMAIL_ADDRESS>Universidad Politécnica de
Madrid, Madrid 28038, Spain  Sumiyajav Sarangerel222B.S. in Computer Science
and Engineering<EMAIL_ADDRESS>Massachusetts Institute of Technology,
Cambridge, Massachusetts 02139, USA  Peng Mun Siew333Research Scientist,
Department of Aeronautics and Astronautics<EMAIL_ADDRESS>Massachusetts
Institute of Technology, Cambridge, Massachusetts 02139, USA  Pablo
Machuca444Visiting Assistant Professor, Department of Aerospace Engineering,
<EMAIL_ADDRESS>San Diego State University, San Diego, California 92182, USA
Daniel Jang555Ph.D. Candidate, Department of Aeronautics and Astronautics,
<EMAIL_ADDRESS>Massachusetts Institute of Technology, Cambridge, Massachusetts
02139, USA  and Richard Linares.666Rockwell International Career Development
Professor and Associate Professor, Department of Aeronautics and Astronautics,
<EMAIL_ADDRESS>Senior Member AIAA. Massachusetts Institute of Technology,
Cambridge, Massachusetts 02139, USA

###### Abstract

With the rapid increase in the number of Anthropogenic Space Objects (ASOs),
Low Earth Orbit (LEO) is facing significant congestion, thereby posing
challenges to space operators and risking the viability of the space
environment for varied uses. Current models for examining this evolution,
while detailed, are computationally demanding. To address these issues, we
propose a novel machine learning-based model, as an extension of the MIT
Orbital Capacity Tool (MOCAT). This advanced model is designed to accelerate
the propagation of ASO density distributions, and it is trained on hundreds of
simulations generated by an established and accurate model of the space
environment evolution. We study how different deep learning-based solutions
can potentially be good candidates for ASO propagation and manage the high-
dimensionality of the data. To assess the model’s capabilities, we conduct
experiments in long term forecasting scenarios (around 100 years), analyze how
and why the performance degrades over time, and discuss potential solutions to
make this solution better.

##  1 Introduction

In recent years, the Low Earth Orbit (LEO) has been experiencing significant
congestion due to the rapid increase in the number of Anthropogenic Space
Objects (ASOs) such as satellites and space debris. This escalating trend is
projected to continue as multiple companies, including SpaceX, Amazon, and
Astra Space, plan to launch large constellations of hundreds to thousands of
satellites. The resulting dense and complex operating environment elevates the
risk of collisions and debris generation, posing substantial challenges for
space operators. Not only does this situation threaten the safety of flight
and mission success in the short run, but it also jeopardizes the long-term
viability of the LEO environment for scientific, commercial, and national
security uses. Hence, understanding and modeling the evolution of the space
environment is crucial for ensuring its sustainability and informing
strategies for effective space traffic management.

A variety of models have emerged to examine this evolution and calculate the
orbital capacity, which is referred to as the number of satellites that can
feasibly be situated in LEO [1]. Traditional proprietary models, developed by
organizations like NASA’s LEGEND [2], ESA’s DELTA [3], JAXA’s IMPACT [4] among
others [5, 6], have been complemented by newer open-source initiatives such as
the MIT Orbital Capacity Tool (MOCAT) and its various versions [7, 8]. Most of
these models operate by propagating all the ASOs forward in time, utilizing
physical models of spacecraft dynamics. This methodology incorporates factors
such as atmospheric drag, solar radiation pressure, third-body perturbations,
and space weather, in addition to simulated collisions and explosions.
Although this technique provides detailed information about individual ASOs,
it is computationally demanding and time-intensive. There are also more recent
approaches that propose a Multi-layer Temporal Network Model (MTN) to explore
space environment dynamics through network interactions among space objects,
offering a different perspective on modeling and analyzing the complex
dynamics of space environment evolution [9].

Figure 1: MOCAT-ML is designed to propagate phase space density distributions
up to epoch $t+n+m$, given a preliminary sequence of distributions up to
$t+n$, which are computed by a reliable and accurate method such as MOCAT-MC.

To address the inherent computational efficiency challenge in these methods,
we suggest a machine learning-based strategy to study the evolution of the
space environment. Our framework, named MOCAT-ML, has been created to extend
the capabilities of the MIT Orbital Capacity Tool. In this work, we show the
current development of a predictive model in MOCAT-ML, specifically aimed at
forecasting ASO density distributions over the phase space of orbital elements
and the area-to-mass ratio $A/m$. It does this by utilizing a dataset of
simulations of the evolution of ASOs in different scenarios and initial
conditions, generated by the validated and precise Monte Carlo variant of
MOCAT, MOCAT-MC [7]. Once trained and validated, the primary purpose of the
new model is to accelerate the propagation of the density distribution of ASOs
for new scenarios. Various data representations, neural architectures, and
training strategies are examined and compared with respect to efficiency and
accuracy, spanning from autoencoders to recurrent neural networks and
transformers. As depicted in Fig. 1, this new model isn’t intended to replace
other models, but to be utilized in a cooperative manner, where the initial
epochs of the evolution are provided by a reliable method (e.g. MOCAT-MC), and
then the machine learning model completes the rest of the propagation.

The rest of the paper is structured as follows: Section 2 provides a
comprehensive overview of the existing literature and foundational concepts in
both the MIT Orbital Capacity Tool (MOCAT) and deep learning, the main
technology utilized in the paper. Section 3 details the dataset utilized in
our study, emphasizing its significance and sources. In Section 4, we
introduce our novel deep learning approaches designed to support the evolution
of the space environment through the prediction of phase space density
distributions. Section 5 presents a thorough examination of the experiments
conducted, including setup, implementation details, and preliminary results,
benchmarking the proposed approaches and testing the best of them for long-
term forecasts. Section 6 discusses the implications of our findings and
potential future research directions. Finally, in Section 7 we summarize the
key insights and contributions of our study.

##  2 Backgrounds

###  2.1 MIT Orbital Capacity Tool

The MIT Orbital Capacity Analysis Tool (MOCAT) is a set of modeling tools
created to understand the long-term sustainability of the LEO environment.
There are two distinct methods – the Source-Sink Evolutionary Model (SSEM) and
the Monte-Carlo (MC) methods. MOCAT-SSEM [10] divides the entire LEO
population into a few object types, which interact at the population level
within each altitude shell. The dynamics that govern the population is
represented by a set of differential equations, which makes this method
computationally efficient. MOCAT-MC [7] is a Monte-Carlo tool that propagates
all LEO objects and models interaction at each timestep. As an MC-based tool,
MOCAT-MC simulates the evolution of the LEO environment and models the
interaction of each object with one another to understand the population
evolution. The Cube method is used for stochastically sample collisions, and
the NASA Standard Breakup model is used for collision debris generation, while
a variety of atmospheric model and propagators can be used. For this work,
MOCAT-MC is used to create a time-series set of propagated orbital objects.

###  2.2 Deep learning concepts relevant to this work

Since the inception of neural networks, both network architecture and training
algorithms have undergone numerous iterations, resulting in substantial
advancements in performance, training efficiency, and data utilization.
Meanwhile, some architectures have been crafted to specialize in handling
specific data types or tasks, leveraging spatial and/or temporal relationships
within the data. Examples of such architectures include Convolutional Neural
Networks (CNNs) [11], Recurrent Neural Networks (RNNs) [12], and Transformer
networks [13]. Deep learning is a natural extension of neural networks, where
multiple neural networks are connected in sequence or in parallel to form deep
architectures. These deep architectures can learn hierarchical representations
of data, with each layer capturing abstractions of the data at different
“depth” levels. CNNs are a class of deep neural networks specially designed to
handle grid-structured data. CNNs utilize small filters to slide or convolve
across the input data, capturing local patterns and features efficiently. On
the other hand, RNNs are designed to handle temporal data, where they process
sequential data by maintaining a hidden state that is updated at each time
step. This mechanism allows RNNs to capture temporal dependencies within the
data. Variations of RNN architectures include Long Short-Term Memory (LSTM)
Networks [14] and Gated Recurrent Units (GRUs) [15]. Transformer networks are
also designed to handle temporal data, but in a more efficient fashion
compared to traditional RNNs. Transformer networks utilize self-attention
mechanisms to weigh the importance of different data points in a sequence,
enabling them to capture dependencies regardless of the distance between data
points. Unlike traditional RNNs, Transformer networks process the entire
sequence of data in parallel, further enhancing their computational
efficiency. In this work, we also utilize an autoencoder, a type of
unsupervised machine learning [16]. An autoencoder consists of an encoder and
decoder pair, with the primary objective of learning an efficient encoding of
the input data. The model aims to facilitate the reconstruction of the input
data from compressed encoding.

##  3 Dataset

As aforementioned, MOCAT-MC is used to propagate both an initial population
and a yearly number of newly launched satellites along time, considering
interactions between objects (i.e., collisions) that result in added debris
objects to the Earth-orbit environment, and natural dynamics such as those
induced by Earth’s oblateness ($J_{2}$ effect) and atmosphere, which result in
changes in each object’s orbital elements and potential decay into the
atmosphere. The initial population is based on two-line element sets (TLEs)
corresponding to January 2023, and the yearly launch rate repeats the
frequency and altitude pattern of launches in the 2018–2022 period. Both the
initial population and number of yearly launches can be multiplied by
arbitrary factors to simulate alternative scenarios. Here, both the initial
population and yearly launches are multiplied by a factor of $15$, which
reduces the sparsity of interactions between objects.

MOCAT-MC considers various stochastic processes such as probability of
collision or number and physical parameters (e.g., mass and size) of newly
generated debris objects. Multiple (Monte Carlo) simulations are then
performed to gain insight into the expected time evolution of a given
scenario. Ultimately, for each simulation of a scenario, MOCAT-MC provides a
sequence of each object’s state (e.g., position, velocity, or orbital elements
along time).

Following a similar representation to that in [17], at each time step, we can
process MOCAT-MC results to generate density distributions of objects in
various phase spaces: e.g., area-to-mass ratio of objects, $A/m$, versus their
perigee, $r_{p}$ (as in Fig. 1), inclination versus apogee, etc., providing a
graphical and intuitive representation of number of objects within certain
phase-space regions, along time. This work primarily focuses on predictions of
the density distribution over the area-to-mass ratio and perigee phase space:
$A/m$ is discretized into 99 evenly distributed bins in a log-scale (from
$10^{-5}$ $m^{2}/kg$ to $10^{1}$ $m^{2}/kg$, which is comprehensive of all
objects in the MOCAT-MC simulations), and $r_{p}$ is discretized into 36 bins
of 50 km in height, from 200 km to 2,000 km altitudes (objects with a perigee
below 200 km are considered to have decayed into the atmosphere, and objects
above 2,000 km are considered to be outside the LEO region and are not
simulated).

Conclusively, the dataset employed in this work consists of multiple Monte
Carlo simulations ($100$) of the same scenario (initial population and launch
rate multiplied by a factor of $15$), where, at each time step within each
simulation, the number of objects within each two-dimensional bin of area-to-
mass ratio and perigee is known.

##  4 Proposed approaches

We frame the problem of forecasting phase space density distributions from a
pure data-driven perspective. No expert knowledge is included in any of the
models at any part of the stage. The common factor among the three approaches
proposed in this section is the use of deep learning. Each of the approaches
explore different representation learning techniques and training strategies
that range from classic architectures for dealing with sequential data such as
the LSTM or the Transformer, to a combination of convolutional and recurrent
layers that exploit the spatio-temporal nature of the data.

###  4.1 Decoupled one-dimensional autoencoder and TST forecaster (Decoupled
1D AE + TST)

Based on the expertise of the authors in deep learning applied to time series
in the space domain [18, 19], we decided to start tackling this problem from a
time series forecasting point of view. However, as seen in section 3, the data
used in this problem is not a one-dimensional time series, but a sequence of
two-dimensional density distributions. Therefore, we have to encode first the
density distributions into one-dimensional vectors.

Figure 2: Encoder-decoder architecture to lower the dimensionality of the
input space before applying a Deep Neural Network (DNN) suitable for time
series forecasting.

As a result, we designed an encoder-decoder architecture, which can be seen in
Fig. 2, where we first encode the sequence of input density distributions into
a lower one-dimensional latent space. Next, a time series forecasting deep
neural network is employed to propagate that latent state into the future.
Finally, a decoder reverts the future states into the corresponding
distributions in the original space of density distributions.

We chose a Transformer architecture for time series forecasting, also known as
TST [20] to be used as the forecaster component. This work is important for
time series forecasting as it introduces a novel, efficient design for
Transformer-based models, leveraging techniques like patching and channel
independence, which have successfully been applied in other fields such as
computer vision. These innovations enhance the model’s ability to perform
long-term forecasting with higher accuracy and efficiency compared to state-
of-the-art models.

As for the encoder-decoder structure itself, we designed a simple autoencoder
that processes one density distribution at a time. The encoder compresses each
two-dimensional density distribution using a sequence of convolutional layers,
while the decoder employs transposed convolutional (or devconvolutional)
layers [21] for upscaling, effectively reconstructing the original two-
dimensional distributions from the latent space representation.

The reason why we call this section “Decoupled" is that we did not train the
two main components of the architecture, namely the autoencoder and the
forecaster, at the same time. Instead, we trained each component separately in
its own training process, keeping the validation set in common between them so
that there is no data leakage from one to another. This approach is simpler
and worked effectively in similar works where we tried to forecast the
atmospheric density following the same strategy [19].

###  4.2 End to end one-dimensional autoencoder and LSTM forecaster (e2e 1D AE
+ LSTM)

In this network, a batch of image sequence serves as the input, traversing
through an encoder layer, LSTM, and decoder layer. Each encoder and decoder
layer is composed of both convolutional and linear components. The
convolutional layers incorporate three Con2d layers, facilitating upsampling
of the input. Subsequently, the output undergoes flattening through the
subsequent linear layer, effectively downsampling the input into the encoding
dimension. In this setup, the encoder processes a batch of image sequences,
treating the input as a matrix with dimensions of batch size x sequence
length. It then encodes this input into a matrix with dimensions of batch size
x sequence length x encoding dimension. This process ensures that the encoder
doesn’t engage with time series data directly; instead, its role is to
transform the input into vectors. Next, LSTM with hidden size 500 is used to
capture the time series part of the data. Its output is in the shape of batch
size x horizon x encoding dimension. Then, it runs through linear and
convolutional parts of the decoder layer. Every convolutional layer has
batchnorm layer and elu activation function. Dropout layers with p = 0.5 are
used in the linear layers of encoder and decoder layer of the architecture.
During the training, we observed a lot of instability and used gradient
clipping to ensure a stable training process. Furthermore, we did not use one
cycle policy which resulted in slower convergence. Instead, training with
constant learning rate that was found from lr-finder produced much quicker
convergence.

###  4.3 End to end bi-dimensional autoencoder and forecaster

In our last proposed model, we integrate the strengths of convolutional and
recurrent neural networks to process spatiotemporal data, keeping the latent
space that goes from the encoder to the decoder as bi-dimensional feature
maps. The input for this architecture is a sequence of density distributions
of length $l$, and, unlike the two approaches presented above, the output in
this case will not be a single density distribution but a sequence of the same
length. In forecasting words, we are setting a horizon for the forecasts equal
to the lookback period. We leave for future work the design of a more flexible
architecture like this that allows the use of different values for the
lookback and horizons.

Figure 3: Overview of the encoder-decoder ConvGRU neural network. Encoder
blocks (three in this example) downsample and encode temporal features, which
are then upsampled and decoded by decoder blocks (also three) to predict the
output sequence. Temporal coherence is maintained by ConvGRU cells, while
spatial resolution is restored via PixelShuffle upsampling.

A graphical view of the architecture can be found in Fig. 3. The architecture,
again, adopts an encoder-decoder structure, where both components contain, as
a key layer, the so-called ConvGRU cells. Each ConvGRU cell in the encoder
processes the spatial data across time, maintaining the spatial structure
while capturing temporal dependencies. This is crucial for understanding
changes in the density distributions over time.

In the encoder design, the flow of data through convolutional layers and
ConvGRU cells is carefullyw orchestrated to capture both spatial and temporal
aspects of the input densities. The encoder processes the input data in
several stages, with each stage consisting of a convolutional operation
followed by a recurrent (ConvGRU) processing step. Initially, the input
sequence is optionally passed through a CoordConv [22] layer, in case spatial
context wants to be emphasized. Following this, the sequence enters the first
stage of the encoder, where a convolutional layer processes each density of
the sequence independently. This step ensures consistent extraction of spatial
features across all time steps, maintaining the integrity of spatial
information in each density. Immediately after the convolutional processing,
the sequence is passed to a ConvGRU cell. Unlike the preceding convolutional
operation, the ConvGRU cell does not treat each density in isolation. Instead,
it integrates information across the sequence, capturing how the spatial
features evolve over time. This recurrent processing is crucial for
understanding temporal dynamics, allowing the model to track changes and
patterns that unfold through the sequence. This process of alternating between
convolutional layers and ConvGRU cells continues through the encoder for
several stages (three by default), progressively refining the representation
of the encoded data.

The decoder, on the other hand, reconstructs the forecasted sequence from the
encoded representations. At each stage, a ConvGRU cell is applied first,
focusing on integrating temporal information across the encoded densities.
Then, the sequence is passed through an upsample block. Each upsample block is
applied independently over the time dimension, and utilizes PixelShuffle
layers, introduced in [23], designed to avoid common issues like
“checkerboard" artifacts in upscaled images. This process of alternating
between ConvGRU cells and upsample blocks is repeated across the decoder
stages.

Finally, the last layer of the architecture takes the output of the decoder
and applies, independently over each time step, a $1\times 1$ convolution to
match the desired spatial dimensions (the same as the input), and the channel
depth, which in our case is one, since we only process one phase space at a
time, the combination of perigee and area-to-mass ratio.

Because in this architecture the horizon of the forecast is greater than one,
the loss function (MSE) is aggregated across all time steps, providing a
comprehensive measure of the model’s performance over the entire sequence.
Across the entire architecture, group normalization [24] is utilized as the
normalization layer, enhancing training stability and model performance
without the dependencies on batch size typical of batch normalization. Indeed,
we found that replacing the normalization layers to batch norm completely
prevented the network from learning anything during training process. The
activation function utilized was ReLU in every layer of the network.

##  5 Experimentation

This section presents the experiments done with the proposed approaches in
section 4. The experimentation is divided into two primary experiments. The
first one focuses on benchmarking the effectiveness of the proposed deep
learning architectures in accurately performing a simple one-step ahead
prediction of the phase space density distribution. The second experiment uses
the best model from the aforementioned benchmark and explor its capability in
long-term forecasting.

We ran every experiment of this study in a GPU node from the MIT Supercloud
[25]. Each node is powered by an Intel Xeon Gold 6248 processor and features
40 CPU cores. The node is equipped with 384 GB of RAM, averaging 9 GB per
core, and includes two Nvidia Volta V100 GPUs, each with 32 GB of RAM. We
implemented the code using the libraries fastai [26] and tsai [27], which are
both written on top of Pytorch. Our source code has been made publicly
available on Github 777https://github.com/ARCLab-MIT/mocat-ml

###  5.1 Benchmark of proposed approaches

In this section, we perform a benchmark among the approaches proposed in
section 4, under the same data and conditions. We utilize a dataset comprising
$100$ simulations from MOCAT-MC. For each of the simulations, we take the
sequence of density distributions over the phase space of the perigee
($r_{p}$) and the area-to-mass ratio ($A_{m}$). Therefore, our data can be
represented as a four-dimensional array with dimensions corresponding to the
number of simulations, time steps, and a discretized density distribution grid
defined by area-to-mass ratio (height of the grid) and perigee (width of the
grid). Each time step in our simulations represents a two-week interval. Since
we are only interested in a side-by-side comparison in this experiment, we
only get the first 2 years of data of each simulation (i.e., around 100
distributions per sequence).

To prepare this dataset for input into our approaches, we employ a sliding
window process. This process creates overlapping windows across the time
dimension of each simulation. Each window encompasses a specific number of
consecutive time steps, defined by the sum of the lookback and the horizon.
The lookback parameter determines the number of past time steps included in
each window, capturing the historical data necessary for the model to learn
temporal patterns. The horizon parameter represents the number of future time
steps the model aims to predict. The choice of lookback and horizon values is
critical, as it balances the amount of historical information the model can
access against its predictive capacity for future states. For this benchmark,
we keep the horizon fixed as one (i.e., the task will be a next-step
prediction), while we do multiple training runs with increasing values of the
lookback, ranging from 2 steps to 16. The stride to perform the sliding window
is fixed as 2 steps. Given that the bi-dimensional ConvGRU approach works only
in a setup where the horizon is equal to the lookback (i.e., it forecasts the
same number of densities it gets), for this experiment we train that approach
restricting the loss function to be computed only in the first element of the
output, to compare side by side the performance in just next-step prediction.

Table 1: Configuration for each of the proposed approaches on the benchmark.
AE stands for AutoEncoder, LR stands for Learning Rate. The number of filters
are given as an array with the value per layer.

(a) Decoupled 1D AE + TST Batch size | 512
---|---
Conv. kernel size | (3,3)
Conv. padding | zero padding
# Filters encoder | [8, 16, 32]
# Filters decoder | [32, 16, 8]
LR | 0.01
N. Epochs | 20
# Parameters (AE) | 1,652,769
# Parameters (TST) | 536,896
TST model type | Plus
Optimizer | Adam
Loss function | MSE

(b) e2e 1D AE + LSTM Batch size | 128
---|---
Conv. kernel size | (3,3)
Conv. padding | zero padding
LSTM hidden size | 500
Dropout rate | 0.5
# Filters encoder | [24, 48, 16]
# Filters decoder | [16, 48, 24]
Norm | BatchNorm
LR | 0.01
# Epochs | 20
# Parameters | 709,985
Act. function | eLU
Optimizer | Adam
Loss function | MSE

(c) e2e 2D ConvGRU Batch size | 32
---|---
Conv. kernel size | (3,3)
ConvGRU kernel size | (5,5)
ConvGRU padding | zero padding
Dropout rate | 0.2
N. Filters encoder | [16, 64, 96]
N. Filters decoder | [96, 64, 16]
LR schedule | 1-cycle
Max. LR | 0.001
N. Epochs | 20
N. Parameters | 4,198,401
Act. function | ReLU
Optimizer | Adam
Loss function | MSE

The individual configuration for each of the proposed approaches can be seen
in Table 1. For the optimizer, i.e., the component that updates the parameters
of the neural network, we employed the popular Adam optimizer [28] in all
cases. For the learning rate, we used a dynamic learning rate with a “1cycle”
policy [29], which schedules the learning rate and momentum with cosine
annealing. The maximum learning rate is set as $1e^{-3}$ in all cases too. As
for the loss function, we employ the Mean Squared Error (MSE), that averages
the difference between the predicted number of ASOs and the real value for
every cell of the density distribution.

Table 2: Comparison of the different methods for varying lookback values, in
terms of the average validation loss at the end of training (MSE) and the
inference time per sample. Cells with best values are highlighted.

L | Approach a (d=64) | Approach a (d=128) | Approach b (d=64) | Approach b (d=128) | Approach d
---|---|---|---|---|---
MSE | Time (ms) | MSE | Time (ms) | MSE | Time (ms) | MSE | Time (ms) | MSE | Time (ms)
1 | 1276.4k±308.3k | 3.38±0.17 | 595.1k±124.6k | 3.46±0.08 | 589.27±51.31 | 1.48±0.52 | 583.88±45.86 | 1.79±0.59 | 47.37±22.66 | 24.95±20.35
2 | 1252.5k±474.8k | 3.72±0.18 | 652.8k±115.5k | 3.76±0.12 | 588.48±43.45 | 1.52±0.50 | 581.85±39.60 | 1.55±0.54 | 49.12±19.10 | 17.66±0.25
4 | 488.8k±238.8k | 3.69±0.12 | 795.1k±586.0k | 3.76±0.12 | 597.49±42.94 | 1.58±0.47 | 570.94±31.22 | 1.61±0.48 | 29.53±8.23 | 27.50±0.72
8 | 690.1k±375.0k | 3.72±0.11 | 340.1k±133.1k | 3.74±0.08 | 493.28±32.89 | 1.80±0.51 | 540.87±37.19 | 1.80±0.55 | 40.98±6.54 | 45.56±0.55
16 | 262.1k±169.2k | 3.78±0.19 | 89.3k±25.6k | 3.75±0.38 | 462.98±16.68 | 2.80±0.52 | 453.70±18.00 | 2.80±0.53 | 45.23±8.28 | 84.77±0.19

The comparison among approaches was done in terms of the value of the
validation loss after the last training epoch (the lower the better). The
validation set was obtained by holding out randomly all the samples for 20% of
the simulations (i.e., 20 in this case). Additionally, we also measure the
inference time per sample to have a sense of how fast this solution could be,
which, at the end, it is the main benefit of this whole work over the use of
just a Monte Carlo solution like the one provided by MOCAT-MC. The results are
shown in Table 2. As it can be seen, there are orders of magnitude of
difference in the MSE among the three approaches, where the decoupled approach
performs the worst and the end-to-end 2D ConvGRU performing best. It is
noteworthy that, even when using a sophisticated and state of the art
Transformer architecture in the forecaster module of the decoupled approach,
the fact that the training was not a single back-propagation but two separated
ones is a key element to failure, regardless of the accuracy of each component
and of the dimensionality of the latent space.

We also see that the evolution of the MSE with respect to the lookback does
not follow a clear decreasing trend as one might expect. In fact, if we focus
on the best performing approach (last column of the table), we see that the
benefit of using more than just the last density to predict the next one is
not aligned with a meaningful decrease of the loss. This might be due to the
low pace with which, naturally, the space environment situation changes in
such a small period of time as two weeks (one step in the dataset). On the
other hand, we find a sweet spot at $L=4$ for approach that achieves the
minimum loss of the entire table. We argue that this is just a matter of the
configuration of the benchmark, and that, given enough time to train, the
configurations with longer horizons would lead to the best performances.

(a) Decoupled 1D AE + TST

(b) e2e 1D AE + LSTM

(c) e2e 2D ConvGRU

Figure 4: Comparison of the predictions (along with last input of the lookback
sequence) of the sample with worst validation loss among the different
proposed approaches

To get a visual grasp on the differences in the quality of each of the
proposed approaches, we show in Fig. 4 the samples of the validation set that
got the largest validation loss (i.e., the worst predictions), for the models
trained with lookback 4. As it can be seen, we can confirm here the poor
performance of the decoupled approach, which makes it valueless even for such
a simple benchmark, where, at the end, just using an identity function would
be a strong baseline. On the other hand, both the e2e 1D AE+LSTM and the e2e
2D ConvGRU approaches seem to respect the shape of the input and create, at
least, faithful reproductions of it, even if we saw in Table 2 that they were
far from each other in terms of MSE.

###  5.2 Long term prediction

In this section, we select the best model from the previous benchmark and
train it again on the full dataset, i.e., on a dataset of 100 simulations with
a total of 2435 epochs per simulation, which span a total of around 93.4 years
of data (the time-step between epochs is two weeks). The selected model is the
e2e 2D ConvGRU, trained with a lookback and a horizon of 4 steps (as a
reminder, this approach always generates the same number of steps it is fed
with). The rest of the training configuration is kept the same as in Table
1(c), with the exception of the stride, for which we employ a value of 4 this
time (therefore there are no overlapping windows in the dataset).
Additionally, this time the loss function (MSE) is computed for the whole
forecasted density sequence, and not only for the first element of the
forecast as it was done before.

Figure 5: Comparative evolution of loss functions for overlapping (2 weeks per
step) and non-overlapping (2 months per step) forecasting strategies: full
period and first year analysis.

Long term predictions are achieved by performing iterative forecasting. This
simple idea consists in calling the model many times in a loop, using in each
iteration the outputs that the model produced in the previous one.

In Fig. 5, we show the evolution of the average validation loss of the
forecasts over time, for two different forecasting strategies, applied on two
different trained models. The first one, that we refer to as “non
overlapping", outputs a sequence of the four predicted distributions that come
directly after the input data (another sequence of four densities), i.e., it
creates a two months-ahead forecast. The “overlapping" strategy, on the other
hand, produces a sequence of four density distributions too, but instead of
being subsequent to the input, they actually overlap with it partially, in a
way that effectively, the predicted data only spans a period- of two weeks
ahead ahead of the last element of the input (i.e., it is effectively a one-
step forecast as we did in the benchmark). This strategy is common in RNNs,
since it allows for fine-grained one-step predictions in models that
architecturally output always the same length as the input (a common pattern
in seq2seq models).

On the left side of the figure, we see how both strategies fail in forecasting
the full period of data (around 100 years), reaching values of the loss six
orders of magnitude larger than the ones we saw in the previous benchmark (see
section 5.1). In fact, only in the first year (right panel of the figure) we
can appreciate a huge increase of the loss function, especially in the
overlapping approach. This can be due to, despite the overlapping task being
easier at each iteration (part of the prediction can be learned from the
leaked input itself), it needs four times more iterations to span the same
period of time than the non-overlapping approach, also consuming four times
more time in the process. This leads us to conclude that, in order to explore
long term forecasting strategies for this problem, we might investigate the
creation of models that forecast data with “gaps" between the input and the
output, in order to reduce the number of iterations needed to span long
periods of time, which is what causes a bigger degradation in this type of
architectures.

(a) Predictions

(b) Targets (truth) for years 1 (26 epochs-ahead), 2 (52 epochs-ahead) and 3
(80-steps ahead)

Figure 7: Iterative yearly forecasts (top) vs true densities (bottom),
produced by a ConvGRU model with a lookback and horizon of four time steps
(with a step size of two weeks). Each epoch represents a two-week step

To have a visual grasp of this degradation of performance in the models on
long-term horizons, we show the predictions for one random simulation of the
validation set in Fig. 7, and contrast them with the real densities at those
epochs. As it can be seen, already in the first year, after 26 epochs (1 year,
considering that each epoch spans two weeks), the background of the density
map starts getting “washed", and we can see that areas which should be
constantly zero across the evolution of the space environment (e.g., edges of
the density map) merge themselves, creating a lose of sharpness that is more
accentuated at epoch 52 (year 2). Finally, the moment we analyze the
predictions on the third year (right side of the figure), the characteristic
shape of the scenario that we are using in this experimentation (See the
description in Section 3) is completely lost. This can explain why the
evolution of the loss we saw in Fig. 5 plateaus after the first years. The
prediction becomes then a blurry distribution where the loss function is
already too high to represent anything meaningful. Further study would be
needed to explain the second bump in the loss in the last iterations of the
curve.

##  6 Perspectives

From a machine learning point of view, there are several challenges associated
with the task of propagating density distributions that we faced in this work,
and that we will face in subsequent iterations of it too:

  * •

The long term horizon of the propagation (up to a hundred years ahead).

  * •

The high-dimensionality of the problem.

  * •

The absence of real-world data and benchmarks to assess the performance and
generalization of the method.

With respect to the first challenge, we showed in our experimentation (See
section 5) how, even a relatively sophisticated architecture such as the end-
to-end 2D ConvGRU (See section 4), which achieves decent performance in the
one-step prediction baseline, fails when used iteratively to span forecasting
periods of more than two years. Needless to say, there is much work to do in
terms of hyperparameter tuning, design of new architectures, and
implementation of more efficient forecasting strategies other than a naive
iteration. However, a different perspective that we plan to implement to solve
this problem is to hybridize the Monte Carlo (MOCAT-MC)/Machine learning
(MOCAT-ML) predictions. This would imply that the first $n$ epochs of the
propagation are provided by MOCAT-MC and the remaining $m$ (with $m>n$) are
provided by MOCAT-ML (See Fig. 1) until the uncertainty of the latter is
estimated to be over a tolerable threshold, where the MC approach could be
used again as a way of recalibrating the process. This opens up a new research
line on how to properly quantify uncertainty in this application of deep
learning models.

Figure 8: Perspectives on how to tackle the high-dimensionality of the
problem: represent directly the input space as a graph, and train using a
Graph Neural Network (GNN).

Concerning the high-dimensional data problem, we employ deep learning
techniques, which are well known for managing learning tasks with non-linear
high-dimensional data such as language, images or video. However, as seen in
the experimentation, and more specifically, in Fig. 7, the surroundings of the
densities, which are constantly zero over time, can cause problems to the
important parts of the phase space if they are not kept as zero during the
propagation, creating artifacts in the main part of the distribution and
eventually leading the forecast to an unrecognizable shape. Given that most of
the phase space is in fact populated with zeros, and therefore, it has a
sparse nature, our perspectives for future iterations of this work go through
the use of graph-based representations. By representing each nonzero value of
the phase space as a node in a graph, where adjacent cells of the phase space
are edges, we can then apply a Graph Neural Network (GNN) that predicts how
that graph will look like in future epochs (See Fig. 8). This method learns
directly from the original input space, taking advantage of the well-known
scalability of GNNs for massive spatio-temporal data [30]

Lastly, in relation to the final challenge previously highlighted, we propose
evaluating MOCAT-ML using simulations of unexplored scenarios distinct from
those utilized in the training set. This is done to gauge the generalizability
of the model. All training and test simulations generated in this work by
MOCAT-MC belong to one single scenario, namely one with an initial population
and launch rate multiplied by a factor of $15$, where, at each time step
within each Monte Carlo simulation, the number of objects within each two-
dimensional bin of area-to-mass ratio and perigee is known. In order to proof
the generalization capabilities of the model, this has to be tested on
different scenarios to the ones used for training, to ensure that it learns
the dynamics behind the evolution of the space environment, and not simply
“memorize" an optimal encoding to reproduce Monte Carlo simulations.

In addition, to objectively measure the model’s performance relative to
published results, we plan to reproduce two test cases examined by Giudici et
al. in [17], and pass them through our approach. Both cases involve
fragmentation events of satellites from Starlink and OneWeb respectively,
where the phase-space density distribution is propagated for up to one year
following the fragmentation.

##  7 Conclusions

This study presents an advancement in the field of space environment modeling
through the creation of MOCAT-ML, a machine learning-based framework designed
to extends the capabilities of the MIT Orbital Capacity Tool (MOCAT) from a
machine learning point of view. In this work, we focus on the predictive
modeling of ASO density distributions in the Low Earth Orbit (LEO)
environment.

The innovative methodologies proposed in this research represent a
breakthrough in addressing the computational efficiency challenge inherent in
traditional space environment models. The experimental results demonstrate
that some of the proposed deep learning architectures, particularly the end-
to-end bi-dimensional autoencoder with ConvGRU cells, have a good performance
in short-term forecasts, which they can deliver in the order of milliseconds.
In contrast, they fail in forecasting densities up to more than two years
ahead, for which several insights and ideas for improvement have been
discussed.

The implications of this study are far-reaching. By improving the accuracy and
reducing the computational demands of ASO density distribution modeling,
MOCAT-ML offers a more efficient and effective tool for space traffic
management and the long-term sustainability of the LEO environment. This is
particularly crucial in the face of escalating challenges posed by the
proliferation of satellites and space debris.

Moreover, our work underscores the potential of hybridizing Monte Carlo and
machine learning predictions for space environment evolution. This perspective
allows for initial epochs of propagation provided by traditional methods like
MOCAT-MC, with subsequent predictions enhanced by the ML model, thereby
offering a balanced strategy that maximizes accuracy and computational
efficiency. We plan on exploring more in depth this integration of the two
models in future iterations of the tool.

Future directions for this research also include refining the machine learning
models to handle long-term forecasting horizons more effectively, testing the
generalization capabilities of the model to unseen scenarios, quantifying
prediction uncertainty, and exploring graph-based representations for the
high-dimensional phase space data. By continuing to innovate in these areas,
we aim to further advance our understanding, scalability and robustness of the
applications of machine and deep learning in this field.

## Acknowledgments

This research was conducted during a fellowship at MIT funded by the
Universidad Politécnica de Madrid’s program for promoting international
academic collaboration in the USA, 2022/2023. The authors acknowledge the MIT
SuperCloud and Lincoln Laboratory Supercomputing Center for providing (HPC,
database, consultation) resources that have contributed to the research
results reported within this paper/report. This work was partially supported
by the National Science Foundation under award NSF-PHY-2028125.

## References

  * Arnas and Linares [2023] Arnas, D., and Linares, R., “Non-Self-Intersecting Trajectories and Their Applications to Satellite Constellation Design and Orbital Capacity,” _Journal of Guidance, Control, and Dynamics_ , Vol. 46, No. 3, 2023, pp. 417–430.  10.2514/1.G007208, URL https://doi.org/10.2514/1.G007208.
  * Liou et al. [2004] Liou, J. C., Hall, D. T., Krisko, P. H., and Opiela, J. N., “LEGEND – a Three-dimensional LEOto-GEO Debris Evolutionary Model,” _Advances in Space Research_ , Vol. 34, 2004, pp. 981–986.  10.1016/J.ASR.2003.02.027.
  * Martin et al. [2004] Martin, C. E., Cheese, J. E., and Klinkrad, H., “Space Debris Environment Analysis with DELTA 2.0,” _55th International Astronautical Congress_ , 2004.  10.2514/6.IAC-04-IAA.5.12.5.04.
  * Sorge and Mains [2016] Sorge, M. E., and Mains, D. L., “IMPACT fragmentation model developments,” _Acta Astronautica_ , Vol. 126, 2016, pp. 40–46.  https://doi.org/10.1016/j.actaastro.2016.04.011, URL https://www.sciencedirect.com/science/article/pii/S0094576515303349, space Flight Safety.
  * Dolado-perez et al. [2013] Dolado-perez, J., Costanzo, R. D., and Revelin, B., “Introducing MEDEE - a New Orbital Debris Evolutionary Model,” _Proc. 6th European Conference on Space Debris_ , 2013, pp. 22–25.
  * Wang and Liu [2019] Wang, X. W., and Liu, J., “An Introduction to a New Space Debris Evolution Model: SOLEM,” _Advances in Astronomy_ , Vol. 2019, 2019.  10.1155/2019/2738276.
  * Jang et al. [2023] Jang, D., Gusmini, D., Siew, P. M., D’Ambrosio, A., Servadio, S., Machuca, P., and Linares, R., “Monte Carlo Methods to Model the Evolution of the Low Earth Orbit Population,” _33rd AAS/AIAA Space Flight Mechanics Meeting_ , Austin, Texas, USA, 2023.
  * D’Ambrosio et al. [0] D’Ambrosio, A., Servadio, S., Mun Siew, P., and Linares, R., “Novel Source–Sink Model for Space Environment Evolution with Orbit Capacity Assessment,” _Journal of Spacecraft and Rockets_ , Vol. 0, No. 0, 0, pp. 1–15.  10.2514/1.A35579, URL https://doi.org/10.2514/1.A35579.
  * Wang et al. [2023] Wang, Y., Wilson, C., and Vasile, M., “Multi-Layer Temporal Network Model of the Space Environment,” _2023 AAS/AIAA Astrodynamics Specialist Conference_ , Aerospace Centre of Excellence, Big Sky, Montana, USA, 2023.
  * Lifson et al. [2023] Lifson, M., Jang, D., Pasiecznik, C., and Linares, R., “MOCAT-SSEM: A Source-Sink Evolutionary Model for Space Debris Environment Evolutionary Modeling,” _9th Annual Space Traffic Management Conference_ , 2023.
  * LeCun et al. [1989] LeCun, Y., et al., “Generalization and network design strategies,” _Connectionism in perspective_ , Vol. 19, No. 143-155, 1989, p. 18.
  * Medsker and Jain [1999] Medsker, L., and Jain, L. C., _Recurrent neural networks: design and applications_ , CRC press, 1999.
  * Vaswani et al. [2017] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I., “Attention is all you need,” _Advances in neural information processing systems_ , Vol. 30, 2017.
  * Hochreiter and Schmidhuber [1997] Hochreiter, S., and Schmidhuber, J., “Long short-term memory,” _Neural computation_ , Vol. 9, No. 8, 1997, pp. 1735–1780.
  * Gers et al. [2000] Gers, F. A., Schmidhuber, J., and Cummins, F., “Learning to forget: Continual prediction with LSTM,” _Neural computation_ , Vol. 12, No. 10, 2000, pp. 2451–2471.
  * Hinton and Salakhutdinov [2006] Hinton, G. E., and Salakhutdinov, R. R., “Reducing the dimensionality of data with neural networks,” _science_ , Vol. 313, No. 5786, 2006, pp. 504–507.
  * Giudici and Colombo [2023] Giudici, L., and Colombo, C., “Space debris density propagation through a finite volume method,” _33rd AAS/AIAA Space Flight Mechanics Meeting_ , Austin, Texas, USA, 2023.
  * Stevenson et al. [2023] Stevenson, E., Rodriguez-Fernandez, V., Urrutxua, H., and Camacho, D., “Benchmarking deep learning approaches for all-vs-all conjunction screening,” _Advances in Space Research_ , Vol. 72, No. 7, 2023, pp. 2660–2675.  https://doi.org/10.1016/j.asr.2023.01.036, URL https://www.sciencedirect.com/science/article/pii/S0273117723000601.
  * Briden et al. [2023] Briden, J., Siew, P. M., Rodriguez-Fernandez, V., and Linares, R., “Transformer-based Atmospheric Density Forecasting,” _24th Advanced Maui Optical and Space Surveillance Technologies_ , arXiv, Maui, Hawaii, United States, 2023.  10.48550/arXiv.2310.16912, URL http://arxiv.org/abs/2310.16912, arXiv:2310.16912 [physics].
  * Nie et al. [2023] Nie, Y., H. Nguyen, N., Sinthong, P., and Kalagnanam, J., “A Time Series is Worth 64 Words: Long-term Forecasting with Transformers,” _International Conference on Learning Representations_ , 2023.
  * Long et al. [2015] Long, J., Shelhamer, E., and Darrell, T., “Fully convolutional networks for semantic segmentation,” _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2015, pp. 3431–3440.
  * Liu et al. [2018] Liu, R., Lehman, J., Molino, P., Petroski Such, F., Frank, E., Sergeev, A., and Yosinski, J., “An intriguing failing of convolutional neural networks and the coordconv solution,” _Advances in neural information processing systems_ , Vol. 31, 2018.
  * Shi et al. [2016] Shi, W., Caballero, J., Huszár, F., Totz, J., Aitken, A. P., Bishop, R., Rueckert, D., and Wang, Z., “Real-time single image and video super-resolution using an efficient sub-pixel convolutional neural network,” _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 1874–1883.
  * Wu and He [2018] Wu, Y., and He, K., “Group normalization,” _Proceedings of the European conference on computer vision (ECCV)_ , 2018, pp. 3–19.
  * Reuther et al. [2018] Reuther, A., Kepner, J., Byun, C., Samsi, S., Arcand, W., Bestor, D., Bergeron, B., Gadepally, V., Houle, M., Hubbell, M., Jones, M., Klein, A., Milechin, L., Mullen, J., Prout, A., Rosa, A., Yee, C., and Michaleas, P., “Interactive supercomputing on 40,000 cores for machine learning and data analysis,” _2018 IEEE High Performance extreme Computing Conference (HPEC)_ , IEEE, 2018, pp. 1–6.
  * Howard and Gugger [2020] Howard, J., and Gugger, S., _Deep Learning for Coders with Fastai and Pytorch: AI Applications Without a PhD_ , O’Reilly Media, Incorporated, 2020.  URL https://books.google.no/books?id=xd6LxgEACAAJ.
  * Oguiza [2023] Oguiza, I., “tsai - A state-of-the-art deep learning library for time series and sequential data,” Github, 2023.  URL https://github.com/timeseriesAI/tsai.
  * Kingma and Ba [2014] Kingma, D. P., and Ba, J., “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * Smith [2017] Smith, L. N., “Cyclical learning rates for training neural networks,” _2017 IEEE winter conference on applications of computer vision (WACV)_ , IEEE, 2017, pp. 464–472.
  * Cini et al. [2022] Cini, A., Marisca, I., Bianchi, F. M., and Alippi, C., “Scalable Spatiotemporal Graph Neural Networks,” _NeurIPS 2022 Temporal Graph Learning Workshop_ , 2022.  URL https://openreview.net/forum?id=UEANz_37Vo.