# Theory of Functional Connections and Nelder-Mead optimization methods
applied in satellite characterization

Allan K. de Almeida Jr<EMAIL_ADDRESS>Safwan Aljbaae
<EMAIL_ADDRESS>Timothée Vaillant<EMAIL_ADDRESS>Jhonathan M.
Piñeros<EMAIL_ADDRESS>Bruno Coelho<EMAIL_ADDRESS>Domingos Barbosa<EMAIL_ADDRESS>Miguel Bergano<EMAIL_ADDRESS>João
Pandeirada<EMAIL_ADDRESS>Francisco C. Carvalho
<EMAIL_ADDRESS>Leonardo B. T. Santos<EMAIL_ADDRESS>Antonio F. B. A. Prado<EMAIL_ADDRESS>Anna Guerman<EMAIL_ADDRESS>Alexandre C. M. Correia<EMAIL_ADDRESS>Instituto de Telecomunicações,
Universidade de Aveiro, 3810-193 Aveiro, Portugal  Postgraduate Division -
National Institute for Space Research (INPE), São Paulo, Brazil  CFisUC,
Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra, Portugal
ICT-UNIFESP, Institute of Science and Technology, São José dos Campos, Brazil
Instituto Superior Técnico, Avenida Rovisco Pais 1, 1049-001 Lisboa Portugal
Polytechnic School, University of Pernambuco, 50720-001 Recife, PE, Brazil
University of Beira Interior, Covilhã, Portugal  IMCCE, UMR8028 CNRS,
Observatoire de Paris, PSL Université, 77 Avenue Denfert-Rochereau, 75014
Paris, France

###### Abstract

The growing population of man-made objects with the build up of mega-
constellations not only increases the potential danger to all space vehicles
and in-space infrastructures (including space observatories), but above all
poses a serious threat to astronomy and dark skies. Monitoring of this
population requires precise satellite characterization, which is is a
challenging task that involves analyzing observational data such as position,
velocity, and light curves using optimization methods. In this study, we
propose and analyze the application of two optimization procedures to
determine the parameters associated with the dynamics of a satellite: one
based on the Theory of Functional Connections (TFC) and another one based on
the Nelder-Mead heuristic optimization algorithm. The TFC performs linear
functional interpolation to embed the constraints of the problem into a
functional. In this paper, we propose to use this functional to analytically
embed the observational data of a satellite into its equations of dynamics.
After that, any solution will always satisfy the observational data. The
second procedure proposed in this research takes advantage of the Nealder-Mead
algorithm, that does not require the gradient of the objective function, as
alternative solution. The accuracy, efficiency, and dependency on the initial
guess of each method is investigated, analyzed, and compared for several
dynamical models. These methods can be used to obtain the physical parameters
of a satellite from available observational data and for space debris
characterization contributing to follow-up monitoring activities in space and
astronomical observatories.

###### keywords:

astrodynamics, numerical methods, satellite characterization

††journal: Acta Astronautica

##  1 Introduction

The number of space activities, including the launch of constellations of
satellites and the associated space debris population, has been steadily
increasing, raising concerns over the impact on the sky quality for
astronomical purposes. This growing population of man-made objects not only
increases the potential danger to all space vehicles and in-space
infrastructures, from expensive astronomical space observatories to
communications satellites, Earth Observation satellite constellations, Space
Stations, and other spacecraft with humans aboard, but also severely impacts
the quality of Earth’s night sky and radio sky. This threat to the orbital
space environment affects those using space assets and those who look through
it as well as work within it [1, 2]. In fact, since 2009, the International
Astronomical Union (IAU)’s Resolution B5, “In Defence of the Night Sky and the
Right to Starlight” (2009), asserts that “an unpolluted night sky that allows
the enjoyment and contemplation of the firmament should be considered a
fundamental socio-cultural and environmental right”. The IAU Centre for the
Protection of Dark and Quiet Sky from Satellite Interference (IAU CPS), co-
hosted by the NSF’s NOIRLab (the US center for ground-based optical
astronomy), and the SKA Observatory (SKAO) was created to respond to the
concerns acutely raised by astronomers after the launch of the first 60
Starlink satellites in May 2019. The IAU CPS assessed the impact of satellite
constellations on astronomy out from the Dark & Quiet Skies for Science and
Society reports [3, 4]. As man-made objects like satellites and space debris
move across the field of view of any astronomical instrument, the astronomical
exposure causes streaks across the images and interferes with astronomical
data. This is a problem that will be exacerbated for astronomical
observatories and many science projects based on wide field instruments and
requiring larger time exposures [1, 3, 5, 6]. For instance, the Zwicky
Transient Facility has already seen an increase in affected images from 0.5%
in late 2019 to 18% in August 2021 [7]. The 3.5 degrees wide-field imager of
the near-future Vera C. Rubin Observatory in Chile will contain at least one
streak in the majority of exposures, further complicated by cascading
electronics effects [8]. SKAO Simulations antenna taking into account main
beams and sidelobes pick up from undesired satellite beacons and radio relay
data emissions back and forth with ground stations suggest that, once the
mature Starlink population is in orbit, every observation in the relevant
bands will, on average, take 70% longer111https://www.nsf.gov/news/special
reports/jasonreportconstellations [1, 9].

Several studies have been used to model the impact of satellite constellations
in wide-field instruments, not only for Optical and Infrared Astronomy [10,
11, 12], but also in Radio Astronomy HI Intensity mapping [13]. Projects and
tools such as the Astrianet telescope sensor
network222https://doi.org/10.18738/T8/GV0ASD provide important information on
LEO objects properties and enable further modeling of satellite or debris
attitude patterns. This is particularly important to assess impacting brief
bright flashes that can occur when a facet or particularly reflective surface
of a satellite briefly reflects more sunlight to an observer on the ground
[14]. Near-future data fusion concepts using space radar and very wide-field
tracking telescopes for Space Surveillance and Tracking (SST), such as the
PASO telescopes [15], will provide additional information that is much needed
to accurately model the orbital parameters and attitudes of satellites and
debris.

Information on the characteristics of a satellite, including space debris, can
be obtained from its dynamics, such as its motion around Earth and rotation,
or from fusion data with light curve [16, 17]. Accurate knowledge of the
specific parameters of a satellite in mathematical models is crucial for
predicting its motion. Numerical optimization can be used to determine the
values of these parameters that best fit both the mathematical models and
external data obtained from observations. To achieve this, an algorithm is
created to find the best solution while adhering to the model and given
constraints by iteratively adjusting the decision variables until the
objective function is minimized or maximized. In general, optimization
processes may be derived from the gradient method, both for linear systems
like the conjugate gradient [e.g. 18] or for nonlinear systems like the
descent gradient [e.g. 19]. The nonlinear least squares (NLLS) is a well-known
gradient based optimization procedure [e.g. 20], which suitable addresses a
large range of astrodynamics problems [21]. Automatic differentiation [22, 23,
24] facilitates the evaluations of the derivatives inherent to the NLLS,
increasing its computational efficiency to solve astrodynamics problems [25].
On the other hand, non-analytic methods can also be used to address
astrodynamics problems. Such kind of optimization procedures do not require
the evaluations of derivatives, which may be advantageous to lower
computational consumption and / or to be applied to non derivable dynamical
models. For instance, the heuristic optimization procedures Simulated
Annealing [e.g. 26] and Genetic Algorithm [e.g. 27] have also been proven
useful to solve astrodynamics problems [28, 29]. The success of the
optimization process depends on the choice of the optimization algorithm, the
analytical model, and the quality of the input data, among others factors.

In this paper, two new applications of two different optimization methods are
proposed to obtain information on the characteristics of a satellite (or space
debris) from observational data. The first method takes advantage of the
Theory of Functional Connections (TFC) [30], and the second method uses a
numerical optimization process through the Nelder-Mead (NM) algorithm [31].

The first method uses the recently developed mathematical framework named the
Theory of Functional Connections to solve constrained problems with high
efficiency [32]. The TFC derives constrained expressions, which are a class of
functionals that satisfy linear constraints [33]. These constrained
expressions are used as functional interpolation to transform a set of
nonlinear differential equations subject to constraints into a new set of
nonlinear differential equations with no constraints [34]. The resulting set
of unconstrained nonlinear differential equations can be solved using
traditional optimization methods, such as the nonlinear least squares method
adopted in this paper. In the application of the TFC presented here, it is
proposed to embed the observational data into the dynamical model. Thus, any
mathematical solution to the dynamical model analytically satisfies the
observational data. Note that the solution includes not only the position as a
function of time, but also parameters of the dynamics directly related to the
properties of the satellite. This is an important state-of-the-art proposed
application of the theory, since it allows for a numerical search (based in
the dynamical models) of the parameters in a subspace that satisfies the
observational data. Although TFC has been applied in astrodynamics to compute
periodic orbits [35, 36] and to evaluate the costs of transfers between the
Earth and Moon [37] and other systems [38], it has not yet been used to solve
the satellite characterization problem.

The second method employed in this study involves using a numerical
optimization process based on the Nelder-Mead algorithm. This algorithm was
first developed in the 1960s [31] and is still being widely used today [e.g.
39]. The NM algorithm is a heuristic and empirical procedure that operates on
a structure called a simplex using various operations such as ordering,
reflection, expansion, contraction, and shrinkage in order to minimize a given
function [40]. In the context of this study, the NM algorithm is used to
minimize the differences between the observational data and the theoretical
predictions obtained from the area-to-mass ratio model. This is achieved
through iterative operations of the NM algorithm, which result in the values
of the model parameters that best fit the observational data. The optimization
process is, for instance, carried out using the optimization functions
provided by the SciPy library [39], a popular tool for scientific computing in
Python, which includes a range of optimization algorithms, including the NM
algorithm.

The methodology employed in this research aims to demonstrate and analyze the
feasibility of applications of the TFC and NM optimization procedures for
satellite characterization. The approach involves simulating observational
data and extracting satellite dynamics parameters from these simulated
observations. To achieve this, a mathematical model describing the satellite
dynamics is utilized within the optimization procedure. Firstly, the
observational data is simulated using both the aforementioned mathematical
model (employed in the TFC and NM optimization procedures) and an external,
well-known, independent source of data. The methodology and mathematical tools
are elaborated in Section 2, which includes a brief explanation of the TFC and
NM methods. Furthermore, the mathematical model utilized to describe the
satellite dynamics around the Earth is presented. In this study, only the two
most significant specific forces, namely the effects resulting from the
geopotential and drag, are considered. This simplification helps in
visualizing the results and analyzing the primary source of errors. It is
important to note that this approach does not compromise generality, as
similar behavior is expected for the efficiency of both the TFC and NM methods
when accounting for other perturbations such as third-body effects and solar
radiation pressure, among others. Section 3 examines the codes based on the
presented TFC and NM methods, evaluating their ability to determine the
satellite’s area-to-mass parameter using a given dataset. A comparison between
the TFC and NM methods is conducted in terms of efficiency, encompassing
accuracy and resource consumption. Lastly, the conclusions derived from the
research are presented in Section 4.

##  2 Mathematical Tools and Methods

In this section, the procedures for obtaining parameters related to the
dynamics of a satellite are presented. These procedures are based on the
application of the TFC and NM optimization methods to satellite
characterization. The methods rely on a mathematical model for the dynamics
and utilize observational data. The purpose of the application of the TFC is
to embed the observational data into a new dynamical model, whose mathematical
form is to be derived by the theory. One key characteristic of this method is
that any potential solution satisfies the observational data analytically. On
the other hand, the proposition based on the NM method offers the advantage of
not requiring the evaluation of the gradient. This means that abrupt changes
in the trajectory, resulting from one or several impulses, or abrupt
variations in the satellite’s parameters can be handled.

###  2.1 The Theory of Functional Connections

Optimization methods, in general, use optimization procedures to satisfy both
the differential equations and the constraints. In contrast, optimization
methods combined with TFC search for solutions in a subspace that always
satisfies constraints analytically [30]. In this paper, it is proposed that
the constraints of the boundary value problem are given by the observational
data, which are then embed into the dynamical model (equations of motion) by
following the outlined procedures. Any subsequent numerical optimization
procedure would then search for solutions in a subspace that satisfies the
observational data analytically.

The system described by the independent variable $t$ and dependent variable
$x(t)$ is subject to a set of $n$ linear constraints on $x$ and its derivative
$\dot{x}$. The general equation to obtain the constrained expression involves
a free function $g(t)$, a set of $n$ linearly independent support functions
$s_{k}(t)$, and unknown functional coefficients $\eta_{k}(t,g(t))$, according
to [32]

$\displaystyle
x(t,g(t))=g(t)+\displaystyle\sum_{k=1}^{n}\eta_{k}(t,g(t))\,s_{k}(t).$ (1)

Using a given set of support functions, the functional coefficients are
determined from the constraints, and the constrained expression is obtained.
For example, when there is only one constraint for the initial time $t_{0}$
given by $x(t_{0})=x_{0}$, the functional coefficient can be obtained from
Eq.(1) as $\eta_{1}(t,g(t))=(x_{0}-g(t_{0}))/s_{k}(t_{0})$. Using the
functional coefficient into Eq.(1), the constrained expression is derived as
$x(t,g(t))=g(t)+(x_{0}-g(t_{0}))s_{k}(t)/s_{k}(t_{0})$.

If the constrained expression is used in an ordinary differential equation
(ODE) given by $\ddot{x}=f(\lambda,\dot{x},x,t)$, where $\lambda$ is a set of
$i+1$ unknown coefficients
($\lambda=\lambda_{0},\lambda_{1},...,\lambda_{i}$), then it becomes
$\ddot{g}=f^{\prime}(\lambda,\dot{g},g,t)$, where
$f^{\prime}(\lambda,\dot{g},g,t)=f(\lambda,\dot{x(t,g(t))},x(t,g(t)),t)$, and
$x(t,g(t))$ is given by the constrained expression. It is important to note
that there is no constraints in the resulting ODE.

In order to numerically solve the ODE, the free function is expressed as a
linear combination of given orthogonal polynomials according to
$g(t)=\displaystyle\Sigma_{j=0}^{m}\xi_{j}\,h_{j}(t)$ truncated at the $m$-th
term, where $\xi_{j}$ are unknown coefficients. After discretizing the
resulting ODE for a set of $N$ values of time $t$, the nonlinear least squares
optimization method is used to obtain the $[(m+1)+(i+1)]$ unknown coefficients
from the $N$ equations generated by the procedure, with $N\geq[(m+1)+(i+1)]$.

###  2.2 The Nelder-Mead optimization procedure

The Nelder-Mead optimization procedure is a popular iterative algorithm used
to minimize an objective function of several variables [31] \- see also [41,
39]. It is considered a direct search algorithm because it does not require
the computation of gradient information, making it useful in cases where the
objective function is non-differentiable or when the gradient is too difficult
to compute.

The algorithm begins by selecting an initial set of points, called a simplex,
to cover a wide range of the parameter space. For example, consider an
objective function with two variables. The simplex would consist of three
points, forming a triangle in the parameter space. The algorithm then
evaluates the objective function at each of the simplex points and identifies
the best point with the smallest function value.

The point with the highest objective function value is then reflected through
the centroid of the remaining points in the simplex. This reflected point is
compared with the other points in the simplex. If the reflected point has a
lower objective function value than the lowest point in the simplex, the
reflected point is expanded further in the same direction to create a new
point. This new point is then evaluated to determine if it is better than the
previous simplex points. If the reflected point does not improve the objective
function value, the simplex is contracted towards the lowest point. This
contraction reduces the size of the simplex in the direction of the highest
point, creating a smaller simplex that is closer to the minimum point.

The algorithm continues to iterate through these steps until a stopping
criterion is met. This criterion may be a maximum number of iterations or a
change small enough in the objective function value, which is defined by a
value for the tolerance to satisfy the minimization procedure. The Nelder-Mead
algorithm is relatively simple to implement and has been shown to be effective
in a wide range of optimization problems. However, it may not converge to the
global minimum of the objective function and may get stuck in local minima. It
is also sensitive to the choice of the initial simplex and the parameters that
control the size of the operations on the simplex. Therefore, careful tuning
of these parameters may be necessary to achieve good performance.

The Nelder-Mead method is a gradient-free optimization algorithm suitable for
problems with a moderate number of parameters and noisy objective functions.
It is effective for non-linear, non-convex, and complex objective functions.
The method is simple to implement and computationally efficient as it does not
require gradient or Hessian calculations. Other optimization algorithms such
as Powell’s method, COBYLA, and Basin Hopping can also provide good results,
but comparing these methods is not the focus of this work.

###  2.3 Dynamical model for satellite characterization

The dynamics of a satellite in orbit around the Earth can be described in the
International Celestial Reference Frame (ICRF) as

$\displaystyle\ddot{{\bm{r}}}={\bm{g}}_{10}-\frac{1}{2}C_{D}\frac{A}{m}\rho
v_{r}{\bm{v}}_{r}+{\bm{P}}.$ (2)

In the above equation, ${\bm{r}}$ denotes the position of the satellite in the
ICRF, and the dots indicate derivatives with respect to time. The term
${\bm{g}}_{10}$ represents the gravitational interaction with the Earth, which
is modeled using the EGM96 model, truncated in the degree and order 10. Its
calculation is performed using the algorithm outlined in [42], along with
another algorithm proposed in this paper. The purpose of incorporating the
additional algorithm is to reduce the computational resources required for
implementing the TFC method code, as presented in A. The area-to-mass ratio
$A/m$ is the transverse area of the satellite divided by its mass. The
coefficient of drag $C_{D}$ is a measurement of the resistance that the
satellite experiences due to the atmosphere (the usual value $C_{D}=2.2$ is
used [43, 44]) and $\rho({\bm{r}},t)$ is the density of the atmosphere, which
is a function of both position and time according to the Jacchia-Roberts
atmospheric dynamic model [45, 46, 47]. The variable ${\bm{v}}_{r}$, with
magnitude $v_{r}$, is the velocity of the satellite relative to the
International Terrestrial Reference Frame (ITRF), which rotates with the Earth
and is connected to the ICRF through the Earth Orientation Parameters [48].
The term ${\bm{P}}$ represents the other perturbative specific forces acting
on the satellite, such as third-body perturbations, solar radiation pressure,
etc. In order to simplify and enhance comprehension of the outcomes obtained
from the application of the two approaches proposed in this study, the term
represented by ${\bm{P}}$ is not considered when generating the results of
this paper.

###  2.4 Source of the simulated observational data

In this paper, the TFC and NM methods are used to estimate the area-to-mass
parameter using simulated observational data. It is important to emphasize
that the application of the TFC and NM optimization procedures is based on Eq.
(2). Therefore, if the simulated observational data is generated also using
the mathematical model described by Eq. (2), it is expected that the methods
will retrieve the area-to-mass parameter with high precision. This scenario
represents an ideal situation where the observational data perfectly aligns
analytically with the model. However, it is worth noting that this ideal
situation does not represent real scenarios, where perturbation forces act on
the satellite and the observational data cannot be analytically described by
the available mathematical model. Consequently, the TFC and NM methods are
also applied in situations where the simulated observational data is generated
independently, rather than using the model described by Eq. (2). In this case,
the General Mission Analysis Tool333The GMAT software is developed and
maintained by NASA’s Goddard Space Flight Center to model, simulate, and
optimize spacecraft trajectories from low Earth orbits to interplanetary
missions. https://software.nasa.gov/software/GSC-17177-1 (GMAT) is used to
simulate the observational data. Although efforts are made to adjust the GMAT
software to generate results as close as possible to the mathematical model
shown in Eq. (2), the results obtained do not perfectly match due to the
complexity of the algorithms involved and the inherent differences between the
mathematical models. Thus, in this paper, the observational data is simulated
using both our data generator \- based on the dynamical model given by Eq. (2)
- and GMAT as reference solutions for comparison. The optimization process
with the generation of simulated observational data using GMAT is described
below:

  1. 1.

An initial (original) value for the area-to-mass ratio is given.

  2. 2.

GMAT is employed with the initial area-to-mass ratio and a set of initial
conditions to generate a set of boundary conditions.

  3. 3.

The set of boundary conditions generated in the previous step is used as
simulated observational data in the TFC and NM optimization processes to
estimate (retrieve) the area-to-mass ratio.

  4. 4.

The estimated (retrieved) area-to-mass ratio is compared with the initial
value given in step 1.

In the case where the observational data is simulated using our data
generator, the process is similar to the one described above, except that, in
step 2, our data generator (based on Eq. 2) is used to generate the set of
boundary conditions, instead of GMAT.

##  3 Results

In this paper, the TFC and Nelder-Mead optimization procedures are used to
determine the parameters of a spacecraft using available position and velocity
data. Specifically, the area-to-mass ratio is estimated using position and
velocity data at two different times. These required data are generated by
integrating initial conditions with a specific value of the area-to-mass ratio
for a given time period, T. Hence, the inputs for the implemented methods are
the initial position and velocity (at $t_{0}=$ “03 Feb 2023 00:00:00 (UTC)”)
and the final position and velocity (at $t_{f}=t_{0}+T$, with $T=3$ minutes),
while the output is the area-to-mass ratio.

The original parameter refers to the value used to simulate the observational
data (using GMAT software or our data generator), while the retrieved
parameter denotes the value obtained by the TFC and NM methods from the
simulated observational data. The disparities between the retrieved area-to-
mass ratio and the original value are analyzed in this section with respect to
altitude, latitude, longitude, and the area-to-mass ratio for both TFC and NM
methods. This difference is better visualized using a relative error index,
defined as

$\displaystyle\text{relative
error}=\bigg{|}\frac{{(A/m)}_{r}-{(A/m)}_{o}}{{(A/m)}_{o}}\bigg{|},$ (3)

where ${(A/m)}_{r}$ represents the retrieved area-to-mass ratio and
${(A/m)}_{o}$ denotes the original area-to-mass ratio used to generate the
data.

It is important to note that altitudes ranging from 150 km to 800 km are
considered a critical range where the atmosphere’s effect can significantly
influence the orbit within a short time [49]. The magnitude of this effect
varies depending on the altitude; satellites below 300 km experience a 12 km
decrease in altitude per day, while those around 800 km may encounter a
similar decrease in altitude for over 20 days [44]. To evaluate the methods’
performance, orbits with altitudes ranging from 150 km to 800 km are chosen,
with a step size of 50 km, for various latitude and longitude values. No
problem of convergence for the initial guess has been observed for both
methods.

###  3.1 Accuracy of the retrieved parameter as function of the altitude

The relative errors (Eq. (3)) are shown in Table 1 for the original area-to-
mass parameters $A/m=1$ m2/kg and $A/m=0.01$ m2/kg. The TFC and NM methods are
compared against the data generated using GMAT and our data generator, with
latitude and longitude equal to zero. As altitude increases, the density of
the atmosphere decreases, and then, also the drag perturbation, which makes
the optimization process of a parameter proportional to the drag more
difficult, and leads to an increase in the relative error shown in Table 1.
Despite using a simple integration algorithm based on a simplified Bulirsch-
Stoer method in the Nelder-Mead method or minimizing residuals through the
nonlinear least squares method in the TFC method, the retrieved parameters are
reasonably close to the original values used by GMAT, showing a maximum
relative error of about 5% and 2% for extreme cases of high altitudes (close
to 800 km) and a small area-to-mass ratio (0.01 m2/kg) for the TFC and NM
methods, respectively. Table 1 shows that the relative error is up to the
order of $10^{-4}$ for altitudes up to 450 km and $10^{-6}$ for a higher area-
to-mass (1 m2/kg) for both the TFC and NM methods. However, when the
observational data is simulated using our data generator, the relative error
is much smaller. Its order of magnitude varies from $10^{-12}$ to $10^{-7}$
for the TFC method and for the original area-to-mass parameter ($A/m=0.01$
m2/kg). The relative error decreases by two orders of magnitude when the
original parameter is $A/m=1$ m2/kg. The relative error for the NM method is
about two orders of magnitude higher compared to the TFC method when the
parameter is retrieved from observational data simulated using our data
generator.

Table 1: Relative error between the retrieved and original values of the area-
to-mass parameter (Eq. 3) obtained using the TFC and NM methods from data of
position and velocity evaluated at two times with difference of 3 minutes
simulated for $A/m$=0.01 m2/kg and $A/m$=1 m2/kg.

| Relative error
---|---
| TFC method | NM method
Altitude | Parameter retrieved from observational simulated data generated by | Parameter retrieved from observational simulated data generated by
(km) | GMAT | our data generator | GMAT |  our data generator
|  $A/m$=0.01 m2/kg |  $A/m$=1 m2/kg |  $A/m$=0.01 m2/kg |  $A/m$=1 m2/kg |  $A/m$=0.01 m2/kg |  $A/m$=1 m2/kg |  $A/m$=0.01 m2/kg |  $A/m$=1 m2/kg
150 | $3.39\times 10^{-07}$ | $3.28\times 10^{-07}$ | $7.41\times 10^{-12}$ | $6.91\times 10^{-14}$ | 5.96$\times 10^{-08}$ | 3.22$\times 10^{-07}$ | 1.51$\times 10^{-09}$ | 4.66$\times 10^{-11}$
200 | $4.80\times 10^{-05}$ | $5.28\times 10^{-05}$ | $5.59\times 10^{-12}$ | $1.08\times 10^{-13}$ | 1.08$\times 10^{-04}$ | 1.09$\times 10^{-04}$ | 1.07$\times 10^{-08}$ | 7.57$\times 10^{-11}$
250 | $1.73\times 10^{-05}$ | $1.64\times 10^{-07}$ | $1.27\times 10^{-10}$ | $2.18\times 10^{-12}$ | 6.20$\times 10^{-06}$ | 2.79$\times 10^{-07}$ | 4.74$\times 10^{-08}$ | 4.66$\times 10^{-11}$
300 | $4.87\times 10^{-05}$ | $1.44\times 10^{-07}$ | $5.94\times 10^{-10}$ | $6.94\times 10^{-12}$ | 2.18$\times 10^{-05}$ | 1.46$\times 10^{-07}$ | 8.41$\times 10^{-08}$ | 2.25$\times 10^{-09}$
350 | $1.19\times 10^{-04}$ | $8.50\times 10^{-07}$ | $1.23\times 10^{-09}$ | $6.42\times 10^{-12}$ | 4.69$\times 10^{-05}$ | 1.66$\times 10^{-07}$ | 1.57$\times 10^{-07}$ | 1.02$\times 10^{-09}$
400 | $2.73\times 10^{-04}$ | $2.38\times 10^{-06}$ | $1.69\times 10^{-09}$ | $2.47\times 10^{-12}$ | 1.17$\times 10^{-04}$ | 6.91$\times 10^{-07}$ | 3.32$\times 10^{-07}$ | 1.03$\times 10^{-08}$
450 | $5.97\times 10^{-04}$ | $5.61\times 10^{-06}$ | $6.83\times 10^{-09}$ | $6.09\times 10^{-11}$ | 2.57$\times 10^{-04}$ | 2.19$\times 10^{-06}$ | 7.23$\times 10^{-07}$ | 1.03$\times 10^{-08}$
500 | $1.26\times 10^{-03}$ | $1.22\times 10^{-05}$ | $2.65\times 10^{-08}$ | $1.14\times 10^{-10}$ | 5.04$\times 10^{-04}$ | 5.13$\times 10^{-06}$ | 7.23$\times 10^{-07}$ | 2.98$\times 10^{-09}$
550 | $2.60\times 10^{-03}$ | $2.56\times 10^{-05}$ | $1.69\times 10^{-10}$ | $6.48\times 10^{-10}$ | 1.01$\times 10^{-03}$ | 1.22$\times 10^{-05}$ | 1.48$\times 10^{-07}$ | 8.31$\times 10^{-08}$
600 | $5.23\times 10^{-03}$ | $5.19\times 10^{-05}$ | $6.11\times 10^{-09}$ | $3.30\times 10^{-10}$ | 2.46$\times 10^{-03}$ | 2.80$\times 10^{-05}$ | 4.74$\times 10^{-08}$ | 7.62$\times 10^{-08}$
650 | $1.01\times 10^{-02}$ | $1.01\times 10^{-04}$ | $1.47\times 10^{-08}$ | $7.89\times 10^{-10}$ | 4.75$\times 10^{-03}$ | 4.43$\times 10^{-05}$ | 1.53$\times 10^{-05}$ | 4.41$\times 10^{-07}$
700 | $1.90\times 10^{-02}$ | $1.90\times 10^{-04}$ | $4.66\times 10^{-08}$ | $1.02\times 10^{-10}$ | 8.17$\times 10^{-03}$ | 8.84$\times 10^{-05}$ | 5.62$\times 10^{-05}$ | 1.07$\times 10^{-06}$
750 | $3.39\times 10^{-02}$ | $3.38\times 10^{-04}$ | $3.31\times 10^{-07}$ | $2.36\times 10^{-09}$ | 1.32$\times 10^{-02}$ | 1.15$\times 10^{-04}$ | 1.07$\times 10^{-08}$ | 5.66$\times 10^{-07}$
800 | $5.65\times 10^{-02}$ | $5.65\times 10^{-04}$ | $1.13\times 10^{-07}$ | $5.04\times 10^{-09}$ | 1.98$\times 10^{-02}$ | 2.36$\times 10^{-04}$ | 1.31$\times 10^{-04}$ | 4.07$\times 10^{-06}$

###  3.2 Accuracy for several values of latitude and longitude

Figure 1 displays the altitude-dependent relative error obtained using the TFC
method (top panel) and NM method (bottom panel) with various initial latitudes
and longitudes. The data was generated through the GMAT software under the
condition of $A/m=1$ m2/kg, and the velocity was selected such that the
initial orbit’s eccentricity is zero, while its direction depends on the
longitude. The TFC and NM methods yield a relative error below 0.01 for
altitudes up to roughly 700 km, with the relative error order of magnitude
ranging from $10^{-7}$ to $10^{-3}$ in this range. However, for altitudes
exceeding 500 km, the relative error increases.

Figure 1: Relative error of the parameter retrieved using the TFC (above) and
the NM (below) methods as function of the altitude for the area-to-mass ratio
$A/m=1$ m2/kg. The observational data are simulated using GMAT.

Figure 2 illustrates the argument for an area-to-mass ratio ($A/m$) of 0.01
m2/kg. The TFC and Nelder-Mead methods produce comparable outcomes. However,
the precision of the recovered parameter’s relative error increases at
altitudes greater than 250 km for both methods. Again, the relative error
rises at 200 km altitude. Notably, the relative error exhibits a two-order-of-
magnitude increase compared to the results showed in Figure 1 for $A/m$ = 1
m2/kg.

Figure 2: Relative error of the parameter retrieved using the TFC (above) and
the NM (below) methods as function of the altitude for the area-to-mass ratio
$A/m=0.01$ m2/kg. The observational data are simulated using GMAT.

Overall, the relative error increases significantly for altitudes exceeding
500 km when $A/m$=1 m2/kg or 250 km when $A/m$=0.01 m2/kg. Trajectories with
low area-to-mass ratios resemble those with zero drag, particularly at high
altitudes. Although differences in relative error exist depending on latitude
or longitude, the trend of increasing error with altitude remains consistent
globally.

###  3.3 Accuracy for observational data simulated using our data generator

The results for observational data simulated using our data generator are
presented in Figure 3 for an area-to-mass ratio of
$A/m=1~{}\text{m}^{2}/\text{kg}$. For the TFC method, the order of magnitude
of the relative error varies from $10^{-14}$ for lower altitudes to $10^{-9}$
for orbits with altitudes close to 800 km. In comparison, the relative error
of the retrieved parameter using the NM method is a few orders of magnitude
higher, as can be seen in Figure 3. Note that the relative errors are much
lower compared to the data generated by GMAT, as expected. This figure
illustrates the accuracy limits of the methods in the ideal scenario where the
mathematical models describe the data with the highest possible accuracy
(machine-level accuracy). It is important to note that even in this case, the
models are not simple, and there are multiple steps of complex numerical
evaluations required to obtain the optimized parameter’s final value. The
relative error as a function of the area-to-mass ratio is shown in Figure 4,
for the case where the initial position has an altitude and longitude equal to
zero and the altitude is 400 km for both TFC (blue) and NM (red) methods using
observational data simulated by our data generator, i.e., a similar model.

Figure 3: Relative error as function of the altitude for an area-to-mass ratio
$A/m=1$ m2/kg evaluated through the TFC (above) and NM (below) methods. The
data is simulated using our data generator. Figure 4: Relative error as
function of the area-to-mass ratio using the TFC (blue) and the NM (red)
methods for an orbit with altitude of 400 km. The data is generated using our
data generator.

###  3.4 Influence of the quality of the observational data over the retrieved
parameter

This section aims to assess the impact of data quality on the performance of
the proposed methods by incorporating simulated noises. To account for various
observation techniques, a random error, denoted as $x$, is introduced to both
the position and velocity data. This error follows a Gaussian distribution
with a probability function given by

$p(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}},$
(4)

where $\mu$ represents the mean and $\sigma$ denotes the standard deviation of
the distribution. The specific standard deviation value depends on the level
of accuracy and the characteristics of the observation system or technique
being simulated. It is important to note that in practice, the accuracy of Low
Earth Orbit (LEO) objects can vary significantly based on factors such as the
instrumentation, tracking method, atmospheric conditions, altitude, hardware
limitations, and other factors.

Satellite Laser Ranging (SLR) utilizes laser pulses to measure the range
between a ground station and a satellite, providing highly accurate
measurements typically within a range of a few centimeters to a few
millimeters in terms of range accuracy [50]. Therefore, for simulating SLR
observations, a possible standard deviation value could be in the range of
0.01 m to 0.1 m for positions. Considering the uncertainties for positions, it
is assumed that the uncertainties for velocity are in the range of 0.001 m/s
to 0.01 m/s.

Global Navigation Satellite Systems (GNSS) such as GPS, GLONASS, and Galileo
offer precise positioning capabilities. The accuracy of GNSS-based position
measurements for LEO satellites can range from one decimeter using GNSS
broadcast ephemerides [51] to a few meters or up to tens of meters, depending
on factors such as the number and quality of received satellite signals,
atmospheric conditions, and receiver hardware [52]. Therefore, for simulating
GNSS-based observations, a standard deviation value in the range of 0.1 to 10
meters could be considered for positions. For the velocity, an uncertainty in
the range of 0.01 m/s to 1 m/s is assumed.

Furthermore, differences between the positions and velocities obtained with
Two-Line Elements (TLEs) when compared with information obtained from GPS is
up to 140 meters in position [53, 54] and from 0.1 m/s [53, pg. 38] up to 140
m/s [54] in velocity.

By incorporating simulated noises into the positional and velocity data, we
can explore the impact of data quality on the performance of the proposed
methods and gain insights into their robustness and effectiveness under
realistic conditions. To simulate real-world scenarios and observe the
influence of the standard deviation $\sigma$ on the retrieved parameter, an
initial position with latitude, longitude, and altitude set to 0°, 0°, and 200
km, respectively, is considered.

Table 2 shows the median and standard deviation of the retrieved parameter
with respect to the original values $A/m=1,\text{m}^{2}/\text{kg}$ and
$A/m=0.01,\text{m}^{2}/\text{kg}$ for several combination values of the
standard deviation $\sigma$ in the range of $10^{-2}$ m to $10^{2}$ m and
$10^{-3}$ m/s to $10^{0}$ m/s for position and velocity, respectively.

For SLR observations, with uncertainties in the range of $10^{-2}$ m or
$10^{-1}$ m for position and $10^{-2}$ m/s for velocity, Table 2 indicates
that the standard deviation for the retrieved parameter is on the order of
$10^{-3}$ m2/kg.

For GNSS-based observations, with uncertainties in the range of $10^{-1}$ m or
$10^{0}$ m for position and $10^{-2}$ m/s or $10^{-1}$ m/s for velocity, the
standard deviation of the retrieved parameter is also on the order of
$10^{-3}$ m2/kg when the uncertainty in the initial velocity is $10^{-2}$ m/s.
This is sufficient to retrieve the area-to-mass ratio from simulated
observational data with a standard deviation lower than the value of the
parameter. However, for larger uncertainties in the velocity, on the order of
$10^{-1}$ m/s, the standard deviation of the retrieved parameter is on the
order of $10^{-2}$ m2/kg, which is lower than the original value of the area-
to-mass ratio for $A/m=1,\text{m}^{2}/\text{kg}$ but higher if the original
parameter is $A/m=0.01,\text{m}^{2}/\text{kg}$. Thus, the data obtained from
GNSS-based observations could still be used to retrieve the parameter for
larger area-to-mass ratios using the information given by the assumed
constraints (two positions and two velocities).

In comparison with SLR and GNSS-based observations, the uncertainties obtained
from TLEs are much higher, typically on the order of $10^{1}$ m for position
and from $10^{-1}$ to $10^{0}$ m/s for velocity. In this case, the standard
deviation of the retrieved area-to-mass parameter is on the order of $10^{-1}$
m2/kg.

As shown in this subsection, the influence of errors due to the real
observations is much larger than the influence of the errors inherent to the
procedures of the numerical applications of the proposed methods (both TFC and
NM). This makes the methods suitable for the purpose of satellite
characterization.

Table 2: The median and standard deviation of the retrieved parameter with
respect to the original area-to-mass $A/m$=1 m2/kg and $A/m$=0.01 m2/kg
obtained with the TFC and NM methods from data with noise. The initial
position is such that its latitude, longitude, and altitude are 0º, 0º, and
200 km, respectively.

Original | Standard deviation of | Parameter retrieved from TFC method | Parameter retrieved from NM method
---|---|---|---
value | the initial distribution |  |
$A/m$ ($\text{m}^{2}/\text{kg}$) | position (m) | velocity (m/s) | Median $\pm$ Standard deviation ($\text{m}^{2}/\text{kg}$) | Median $\pm$ Standard deviation ($\text{m}^{2}/\text{kg}$)
| $1\times 10^{2}$ | $1\times 10^{0}$ | $0.94\pm 0.52$ | $0.98\pm 0.75$
| $1\times 10^{1}$ | $1\times 10^{0}$ | $0.94\pm 0.41$ | $0.99\pm 0.56$
| $1\times 10^{0}$ | $1\times 10^{0}$ | $1.10\pm 0.44$ | $1.01\pm 0.52$
1 | $1\times 10^{0}$ | $1\times 10^{-1}$ | $1.002\pm 0.043$ | $1.000\pm 0.062$
| $1\times 10^{-1}$ | $1\times 10^{-1}$ | $1.005\pm 0.046$ | $1.013\pm 0.059$
| $1\times 10^{-1}$ | $1\times 10^{-2}$ | $1.0007\pm 0.0045$ | $1.0005\pm 0.0064$
| $1\times 10^{-2}$ | $1\times 10^{-2}$ | $0.9994\pm 0.0041$ | $1.0000\pm 0.0062$
| $1\times 10^{-2}$ | $1\times 10^{-3}$ | $1.00002\pm 0.00045$ | $1.00004\pm 0.00061$
| $1\times 10^{0}$ | $1\times 10^{-2}$ | $(1.06\pm 0.42)\times 10^{-2}$ | $(1.16\pm 0.74)\times 10^{-2}$
0.01 | $1\times 10^{-1}$ | $1\times 10^{-2}$ | $(0.90\pm 0.44)\times 10^{-2}$ | $(0.96\pm 0.56)\times 10^{-2}$
| $1\times 10^{-2}$ | $1\times 10^{-2}$ | $(1.00\pm 0.47)\times 10^{-2}$ | $(1.10\pm 0.57)\times 10^{-2}$
| $1\times 10^{-2}$ | $1\times 10^{-3}$ | $(1.002\pm 0.049)\times 10^{-2}$ | $(0.995\pm 0.059)\times 10^{-2}$

###  3.5 Simulated observational data given only by positions

The results of applying the optimization procedures developed in this paper
are shown above for the case where there is information on the initial
position and velocity, as well as the position and velocity after 3 minutes.

In this subsection, the constraints are limited to four positions evaluated at
four different times: 0, T/3, 2T/3, and T. It should be noted that in this
case, the codes retrieve the parameter assuming that no information on the
velocity is contained in the simulated observational data. The optimization
codes estimate the three components of the velocity together with the area-to-
mass parameter. This consideration increases the execution time of the Nelder-
Mead algorithm, but it does not affect the execution time of the TFC-based
method.

Table 3 demonstrates that the standard deviation values of the retrieved
parameters are similar in magnitude to those shown in Table 2 for the case
where the standard deviation of velocity (in m/s) is one order of magnitude
lower than the standard deviation of position (in m). The retrieved values are
one order of magnitude lower than those shown in Table 2 when the standard
deviation of the initial velocity distribution (in m/s) is of the same
magnitude as the standard deviation of the initial position distribution (in
m). Hence, if the standard deviation of velocity (in m/s) is of the same
magnitude or higher than the standard deviation of position (in m), it is more
accurate to retrieve the parameter using only position information.

Table 3: The median and standard deviation of the retrieved parameter with
respect to the original area-to-mass $A/m$=1 m2/kg and $A/m$=0.01 m2/kg
obtained with the TFC and NM methods from data with noise. The simulated
observational data are given by 4 values of position (for
$t=t_{0},t_{0}+T/3,t_{0}+2T/3,t_{0}+T$, with $T=3$ minutes). The initial
position is such that its latitude, longitude, and altitude are 0º, 0º, and
200 km, respectively.

Original | Standard deviation of | Parameter retrieved from TFC method | Parameter retrieved from NM method
---|---|---|---
value | the initial distribution |  |
$A/m$ ($\text{m}^{2}/\text{kg}$) | position (m) | Median $\pm$ Standard deviation ($\text{m}^{2}/\text{kg}$) | Median $\pm$ Standard deviation ($\text{m}^{2}/\text{kg}$)
| $1\times 10^{1}$ | 1.01$\pm 0.17$ | 1.00$\pm 0.13$
1 | $1\times 10^{0}$ | 0.999$\pm 0.014$ | 0.999$\pm 0.015$
| $1\times 10^{-1}$ | 0.9999$\pm 0.0016$ | 1.0002$\pm 0.0016$
| $1\times 10^{-2}$ | 0.99999$\pm 0.00018$ | 1.00000$\pm 0.00071$
| $1\times 10^{-3}$ | 1.000000$\pm 0.000015$ | 1.000000$\pm 0.000124$
0.01 | $1\times 10^{-1}$ | $(1.00\pm 0.13)\times 10^{-2}$ | $(1.01\pm 0.17)\times 10^{-2}$
| $1\times 10^{-2}$ | $(0.999\pm 0.014)\times 10^{-2}$ | $(1.004\pm 0.015)\times 10^{-2}$
| $1\times 10^{-3}$ | $(1.0000\pm 0.0014)\times 10^{-2}$ | $(1.0001\pm 0.0015)\times 10^{-2}$

###  3.6 Execution time

The results for the TFC method were obtained using the Python programming
language with the assistance of the TFC module [55] that uses automatic
differentiation and a just-in-time (JIT) compiler [56]. On the other hand, the
results for the NM method were obtained using the Python programming language
with the help of the SciPy library [39]. Since every iteration in the NM
procedure depends on the results of the integration of a set of initial
conditions, its performance is highly dependent on the numerical integration
method adopted into the procedure. Several integration methods were tested,
like the Radau or Bulirsch–Stoer, and the best efficiency was obtained with a
Runge-Kutta method of order 8 (DOP853) [e.g. 57]. The computational speed of
the TFC and NM procedures were tested on a machine equipped with an i7-4790
processor. To analyze the speed of the methods, the computational time
required to retrieve the 324 area-to-mass ratios from a set of simulated
observational data (used to generate Fig. 1) was measured. The code based in
the TFC retrieved the 324 values of the area-to-mass parameter in 1.836634
seconds. In contrast, the code based on the NM method retrieved the 324 values
of the area-to-mass parameter values from the same set of simulated
observational data (used for Fig. 1) in 621.743514 seconds.

##  4 Conclusions

This paper presents techniques for finding the parameters of the dynamical
model of a satellite from observational data. The study adopts the TFC (Theory
of Functional Connections) and NM (Nelder-Mead) simplex algorithm as
optimization methods to determine the optimal area-to-mass ratio of a
satellite using an inverse problem approach. The accuracy of these methods is
evaluated for different altitudes and initial conditions, and the results are
compared with those generated using GMAT (General Mission Analysis Tool). The
Nelder-Mead method is chosen as one of the optimization methods due to its
suitability for solving the orbit determination problem with a moderate number
of parameters and its ability to handle noisy objective functions, which is
crucial when dealing with real-world observations and measurements. On the
other hand, TFC is selected for its efficiency, potential in solving complex
optimization problems with a large number of variables, and its ability to
incorporate observational data into the dynamical model.

Comparisons with observational data simulated using our data generator
demonstrate the highest level of accuracy achievable by these methods when the
constraints align with the dynamical model. In this scenario, the TFC method
exhibits higher accuracy than the NM method, with the retrieved parameter
being orders of magnitude closer to the original parameter. However, this high
level of accuracy may not be significant when the observational data does not
perfectly fit the dynamical model, as is the case with real satellites subject
to imperfectly known perturbations. To address this problem, the paper
analyzes the case where the constraints are not perfect solutions of the
dynamical model using observational data simulated with GMAT. In this case,
both the TFC and NM methods demonstrate a similar and reasonably good
accuracy, albeit inferior to the observational data generated by our data
generator. This implies that both the TFC and NM methods proposed in this
paper have the potential to accurately retrieve the parameters of a real
satellite, primarily depending on external factors such as the quality of the
observational data and the dynamical model, rather than the method itself. The
paper also examines the computational resource consumption of the methods. The
TFC procedure has shown a few orders of magnitude faster convergence than the
NM method to retrieve the parameter.

## Acknowledgements

We acknowledge support from the National Institute for Space Research - INPE.
This work is supported by the European Regional Development Fund (FEDER),
through the Competitiveness and Internationalization Operational Programme
(COMPETE 2020) of the Portugal 2020 framework [Project SmartGlow with Nr.
069733 (POCI-01-0247-FEDER-069733)]; The team acknowledges further support
from ENGAGE-SKA Research Infrastructure, ref. POCI-01-0145-FEDER-022217,
funded by COMPETE 2020 and FCT, Portugal; IT team members acknowledge support
from Projecto Lab. Associado UID/EEA/50008/2019. We acknowledge support by the
European Commission H2020 Programme under the grant agreement 2-3SST2018-20.
Support from Center for Mechanical and Aerospace Science and Technologies -
C-MAST, funded by FCT-Fundação para a Ciência e a Tecnologia through project
UIDB/00151/2020. This work is also supported by CFisUC (UIDB/04564/2020 and
UIDP/04564/2020). This work is supported by the Fundação para a Ciência e
Tecnologia (FCT), Ph.D. grant No.2022.12341.BDANA.

## References

  * [1] A. Lawrence, M. L. Rawls, M. Jah, A. Boley, F. Di Vruno, S. Garrington, M. Kramer, S. Lawler, J. Lowenthal, J. McDowell, M. McCaughrean, The case for space environmentalism, Nature Astronomy 6 (2022) 428–435.  arXiv:2204.10025, doi:10.1038/s41550-022-01655-6.
  * [2] D. Barbosa, B. Coelho, M. Bergano, C. Magalhães, D. Mendonça, D. Silva, A. C. M. Correia, J. Pandeirada, V. Ribeiro, T. Esposito, F. Marchis, Cyber-Cosmos: A new citizen science concept in a dark sky destination, Acta Astronautica 200 (2022) 612–619.  arXiv:2111.03458, doi:10.1016/j.actaastro.2022.09.008.
  * [3] C. Walker, S. Di Pippo, M. Aubé, J. Barentine, Z. Benkhaldoun, P. Benvenuti, C. Bouroussis, R. Green, J. Hearnshaw, A. Mackenzie, H. Liszt, J. D. Lowenthal, C. Muñoz-Tuñón, K. Nield, N. Ricard, J. M. Rodriguez Espinosa, P. Sanhueza, A. M. Varela Pérez, A. Williams, S. Lau, C. Muñoz Tuñón, J. M. Rodríguez Espinosa, C. Muñoz-Tuñon, J. P. Armas, C. Benn, P. Blattner, J. Diaz-Castro, M. Donners, B. Douglas, D. Gašparovský, J. Hall, Z. Kolláth, T. Novak, S. Ortolani, J. M. Rodríguez Espinosa, J. Lowenthal, S. Bara, A. Jägerbrand, A. Jechow, T. Longcore, M. Motta, L. Schlangen, S. Schroer, S. Lucatello, D. Galadí-Enríquez, G. Micheva, A. Otarola, L. Allen, M. Bannister, M. Barstow, C. Bassa, R. Davies, D. Devost, O. Hainaut, M. Jah, K. Krafton, H. Krantz, D. Kucharski, T. Maclay, R. Massey, J. Mudd, Charles, J. McDowell, J. Parriott, E. Pearce, M. Rawls, G. Rotola, R. Seaman, P. Seitzer, J. Siminski, L. Storrie-Lombardi, R. Street, J. Tregloan-Reed, J. A. Tyson, E. Unda-Sanzana, R. Wainscoat, O. Zamora, Dark & Quiet Skies I (2020), Dark & Quiet Skies I (2020), Report of the conference held 5-9 October, 2020. (Oct. 2020).  doi:10.5281/zenodo.5898785.
  * [4] C. Walker, S. Di Pippo, M. Aubé, J. Barentine, Z. Benkhaldoun, P. Benvenuti, C. Bouroussis, F. Di Vruno, R. Green, J. Hearnshaw, M. Aoraki, H. Liszt, J. Lowenthal, C. Muñoz-Tuñón, K. Nield, M. Ohishi, N. Ricard, J. M. Rodriguez Espinosa, G. Rotola, P. Sanhueza, A. Varéla Pérez, A. Williams, P. Cameron, A. Hänel, R. Hartley, A. Krop-Benesch, A. Lien, M. Motta, S. J. Ribas, M. Rol, Y. Yakushina, K. M. Zielinska-Dabkowska, G. Zhao, T. Aganaba, A. Boley, E. Cirkovic, A. Falle, S. Freeland, L. L. Hallet, M. Hanlon, R. Jakhu, L. Jamschon MacGarry, C. Johnson, M. Kazlouskaya, A.-S. Martin, T. Masson-Zwaan, N. Mishra, R. Moro-Aguilar, C. Mudd, A. Simon-Butler, C. Van Eijk, R. Venn, Z. Wager, S. Alis, S. Berry, T. Bharadwaj, J. Davis, E. Lagadec, L. Le Breuilly, S. Lucatello, R. Massey, J. Parriott, L. Valdivia, C. Weeden, C. Hofer, L. Andrade, J. Askari, S. Baruch, J. Bohnlein, M. Dunn, J. Grunsfeld, G. Jansson, D. Knox, V. Raval, R. Sedwick, J. Smith, M. Vanotti, C. Bassa, C. Cui, E. D. Galadi, O. Hainaut, P. Hasan, F. Kébé, S. Lawler, J. McDowell, A. Otarola, M. Rawls, P. Seitzer, R. Street, C. Pree, T. Gergely, G. Hellbourg, B. Indermuehle, M. Lindqvist, A. Tzioumis, B. A. Vanderley, B. Winkel, H. Zhang, Dark & Quiet Skies II (2021), Dark & Quiet Skies II (2021), Report of the conference held 3-7 October, 2021. (Oct. 2021).  doi:10.5281/zenodo.5874725.
  * [5] J. Hall, C. Walker, M. Rawls, J. McDowell, R. Seaman, A. Venkatesan, J. Lowenthal, R. Green, K. Krafton, J. Parriott, SATCON2: Executive Summary, in: Bulletin of the American Astronomical Society, 2021, p. 2.205.  doi:10.3847/25c2cfeb.4554c01f.
  * [6] M. L. Rawls, D. Barron, I. Birdwell, E. Cirkovic, T. Deck, F. Di Vruno, V. Gokhale, M. Goodman, S. Kafka, F. Kebe, D. Knox, H. Krantz, S. Kruk, S. Lawler, D. Monet, M. Peel, J. Tregloan-Reed, O. Zamora, L. Allen, C. Walker, J. Hall, SATCON2: Observations Working Group Report, in: Zenodo report, Vol. 56, 2021, p. 5608826.  doi:10.5281/zenodo.5608826.
  * [7] P. Mróz, A. Otarola, T. A. Prince, R. Dekany, D. A. Duev, M. J. Graham, S. L. Groom, F. J. Masci, M. S. Medford, Impact of the SpaceX Starlink Satellites on the Zwicky Transient Facility Survey Observations, The Astrophysical Journal Letters 924 (2) (2022) L30.  arXiv:2201.05343, doi:10.3847/2041-8213/ac470a.
  * [8] J. A. Tyson, Ž. Ivezić, A. Bradshaw, M. L. Rawls, B. Xin, P. Yoachim, J. Parejko, J. Greene, M. Sholl, T. M. C. Abbott, D. Polin, Mitigation of LEO Satellite Brightness and Trail Effects on the Rubin Observatory LSST, The Astronomical Journal 160 (5) (2020) 226.  arXiv:2006.12417, doi:10.3847/1538-3881/abba3e.
  * [9] A. Lawrence, Losing the Sky, Edinburgh: Photon Productions, 2021.
  * [10] C. G. Bassa, O. R. Hainaut, D. Galadí-Enríquez, Analytical simulations of the effect of satellite constellations on optical and near-infrared observations, Astronomy & Astrophysics 657 (2022) A75.  arXiv:2108.12335, doi:10.1051/0004-6361/202142101.
  * [11] A. Mallama, Starlink Satellite Brightness – Characterized From 100,000 Visible Light Magnitudes, arXiv e-prints (2021) arXiv:2111.09735arXiv:2111.09735, doi:10.48550/arXiv.2111.09735.
  * [12] A. Mallama, The Sky Distribution and Magnitudes of Starlink Satellites by the Year 2027, arXiv e-prints (2022) arXiv:2209.12060arXiv:2209.12060, doi:10.48550/arXiv.2209.12060.
  * [13] S. E. Harper, C. Dickinson, Potential impact of global navigation satellite services on total power H I intensity mapping surveys, Monthly Notices of the Royal Astronomical Society 479 (2) (2018) 2024–2036.  arXiv:1803.06314, doi:10.1093/mnras/sty1495.
  * [14] D. Kucharski, G. Kirchner, M. K. Jah, J. C. Bennett, F. Koidl, M. A. Steindorfer, P. Wang, Full attitude state reconstruction of tumbling space debris TOPEX/Poseidon via light-curve inversion with Quanta Photogrammetry, Acta Astronautica 187 (2021) 115–122.  doi:10.1016/j.actaastro.2021.06.032.
  * [15] B. Coelho, D. Barbosa, M. Bergano, J. Pandeirada, P. Marques, A. C. M. Correia, J. Matias de Freitas, Developing a data fusion concept for radar and optical ground based SST station, arXiv e-prints (2022) arXiv:2211.04443arXiv:2211.04443, doi:10.48550/arXiv.2211.04443.
  * [16] M. Jah, R. Madler, Satellite Characterization: Angles and Light Curve Data Fusion for Spacecraft State and Parameter Estimation, in: S. Ryan (Ed.), Advanced Maui Optical and Space Surveillance Technologies Conference, 2007, p. E49.
  * [17] R. Linares, M. K. Jah, J. L. Crassidis, C. K. Nebelecky, Space object shape characterization and tracking using light curve and angles data, Journal of Guidance, Control, and Dynamics 37 (1) (2014) 13–25.  arXiv:https://doi.org/10.2514/1.62986, doi:10.2514/1.62986.
URL https://doi.org/10.2514/1.62986

  * [18] M. R. Hestenes, E. Stiefel, et al., Methods of conjugate gradients for solving linear systems, Journal of research of the National Bureau of Standards 49 (6) (1952) 409–436.
  * [19] H. B. CURRY, The method of steepest descent for non-linear minimization problems, Quarterly of Applied Mathematics 2 (3) (1944) 258–261.
URL http://www.jstor.org/stable/43633461

  * [20] C. L. Lawson, R. J. Hanson, Solving Least Squares Problems, Society for Industrial and Applied Mathematics, 1995.  arXiv:https://epubs.siam.org/doi/pdf/10.1137/1.9781611971217, doi:10.1137/1.9781611971217.
URL https://epubs.siam.org/doi/abs/10.1137/1.9781611971217

  * [21] B. Tapley, B. Schutz, G. Born, in: B. D. Tapley, B. E. Schutz, G. H. Born (Eds.), Statistical Orbit Determination, Academic Press, Burlington, 2004, pp. xi–xv.  doi:https://doi.org/10.1016/B978-012683630-1/50019-9, [link].
URL https://www.sciencedirect.com/science/article/pii/B9780126836301500199

  * [22] R. D. Wilkins, Investigation of a new analytical method for numerical derivative evaluation, Commun. ACM 7 (8) (1964) 465–471.  doi:10.1145/355586.364792.
URL https://doi.org/10.1145/355586.364792

  * [23] R. E. Wengert, A simple automatic derivative evaluation program, Commun. ACM 7 (8) (1964) 463–464.  doi:10.1145/355586.364791.
URL https://doi.org/10.1145/355586.364791

  * [24] A. Griewank, et al., On automatic differentiation, Mathematical Programming: recent developments and applications 6 (6) (1989) 83–107.
  * [25] A. Bani Younes, M. Alhulayil, J. Turner, A. M. Atallah, High order nonlinear least-squares for satellite pose estimation, The Journal of the Astronautical Sciences 70 (15) (2023).  doi:10.1007/s40295-023-00378-8.
  * [26] H. Szu, R. Hartley, Fast simulated annealing, Physics Letters A 122 (3) (1987) 157–162.  doi:https://doi.org/10.1016/0375-9601(87)90796-1.
URL https://www.sciencedirect.com/science/article/pii/0375960187907961

  * [27] M. Mitchell, An Introduction to Genetic Algorithms, The MIT Press, 1998.  doi:10.7551/mitpress/3927.001.0001.
URL https://doi.org/10.7551/mitpress/3927.001.0001

  * [28] S. W. Paek, S. Kim, O. de Weck, Optimization of reconfigurable satellite constellations using simulated annealing and genetic algorithm, Sensors 19 (4) (2019).  doi:10.3390/s19040765.
URL https://www.mdpi.com/1424-8220/19/4/765

  * [29] T. Savitri, Y. Kim, S. Jo, H. Bang, Satellite constellation orbit design optimization with combined genetic algorithm and semianalytical approach, International Journal of Aerospace Engineering 2017 (17) (2017) Article ID 1235692\.  doi:10.1155/2017/1235692.
  * [30] D. Mortari, The Theory of Connections: Connecting Points, Mathematics 5 (4), https://www.mdpi.com/2227-7390/5/4/57 (2017).  doi:10.3390/math5040057.
  * [31] J. A. Nelder, R. Mead, A Simplex Method for Function Minimization, The Computer Journal 7 (4) (1965) 308–313.  arXiv:https://academic.oup.com/comjnl/article-pdf/7/4/308/1013182/7-4-308.pdf, doi:10.1093/comjnl/7.4.308.
URL https://doi.org/10.1093/comjnl/7.4.308

  * [32] D. Mortari, C. Leake, The Multivariate Theory of Connections, Mathematics 7 (3) (2019) 296/1–22, https://www.mdpi.com/2227-7390/7/3/296.  doi:10.3390/math7030296.
  * [33] C. Leake, H. Johnston, D. Mortari, The multivariate theory of functional connections: Theory, proofs, and application in partial differential equations, Mathematics 8 (8) (2020) 1303, http://dx.doi.org/10.3390/math8081303.  doi:10.3390/math8081303.
  * [34] D. Mortari, H. Johnston, L. Smith, High Accuracy Least-squares Solutions of Nonlinear Differential Equations, Journal of Computational and Applied Mathematics 352 (2019) 293 – 307, http://www.sciencedirect.com/science/article/pii/S0377042718307325.  doi:https://doi.org/10.1016/j.cam.2018.12.007.
  * [35] H. Johnston, M. W. Lo, D. Mortari, A functional interpolation approach to compute periodic orbits in the circular-restricted three-body problem, Mathematics 9 (11) (2021).  doi:10.3390/math9111210.
URL https://www.mdpi.com/2227-7390/9/11/1210

  * [36] A. K. de Almeida, A. F. Prado, D. Mortari, Using the theory of functional connections to create periodic orbits with a linear variable thrust, New Astronomy 104 (2023) 102068.  doi:https://doi.org/10.1016/j.newast.2023.102068.
URL https://www.sciencedirect.com/science/article/pii/S1384107623000696

  * [37] A. K. de Almeida Junior, H. Johnston, C. Leake, D. Mortari, Fast 2-impulse non-keplerian orbit transfer using the theory of functional connections, European Physical Journal Plus 136 (2021).  doi:10.1140/epjp/s13360-021-01151-2.
  * [38] A. K. de Almeida Junior, A. F. B. de Almeida Prado, Comparisons between the circular restricted three-body and bi-circular four body problems for transfers between the two smaller primaries, Scientific Reports 12 (4148) (2022) 1–19, https://doi.org/10.1038/s41598-022-08046-x.  doi:10.1038/s41598-022-08046-x.
  * [39] P. Virtanen, R. Gommers, T. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. van der Walt, M. Brett, J. Wilson, J. Millman, N. Mayorov, J. Jokela, A. Virtanen, R. R. Gommers, T. Oliphant, Scipy, https://www.scipy.org/ (2020).
  * [40] W. Spendley, G. R. Hext, F. R. Himsworth, Sequential application of simplex designs in optimisation and evolutionary operation, Technometrics 4 (4) (1962) 441–461.  arXiv:https://www.tandfonline.com/doi/pdf/10.1080/00401706.1962.10490033, doi:10.1080/00401706.1962.10490033.
URL https://www.tandfonline.com/doi/abs/10.1080/00401706.1962.10490033

  * [41] W. Press, S. Teukolsky, W. Vetterling, B. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd Edition, Cambridge University Press, 2007.
  * [42] H. Kuga, V. Carrara, Fortran and C-codes for higher order and degree geopotential and derivatives computation, in: Proceedings of XVI SBSR-Brazilian Symposium on Remote Sensing - Iguacu Falls, Brazil, 2013.
  * [43] G. Cook, Satellite drag coefficients, Planetary and Space Science 13 (10) (1965) 929–946.  doi:https://doi.org/10.1016/0032-0633(65)90150-9.
URL https://www.sciencedirect.com/science/article/pii/0032063365901509

  * [44] E. M. Gaposchkin, A. J. Coster, Analysis of satellite drag, Lincoln Laboratory Journal 1 (1988) 203–224.
  * [45] L. G. Jacchia, Static Diffusion Models of the Upper Atmosphere with Empirical Temperature Profiles, SAO Special Report 170 (Dec. 1964).
  * [46] L. G. Jacchia, Static Diffusion Models of the Upper Atmosphere with Empirical Temperature Profiles, Smithsonian Contributions to Astrophysics 8 (1965) 215.
  * [47] C. E. Roberts, An analytic model for upper atmosphere densities based upon jacchia’s 1970 models, Celestial Mechanics 4 (3-4) (1971) 368–377.  doi:10.1007/bf01231398.
  * [48] P. K. Seidelmann, 1980 iau theory of nutation: The final report of the iau working group on nutation, Celestial Mechanics 27 (1982) 79–106.  doi:10.1007/BF01228952.
  * [49] J. Ashenberg, On the effects of time-varying aerodynamical coefficients on satellite orbits, Acta Astronautica 38 (2) (1996) 75–86.  doi:https://doi.org/10.1016/0094-5765(96)00002-1.
URL https://www.sciencedirect.com/science/article/pii/0094576596000021

  * [50] M. R. Pearlman, J. J. Degnan, J. M. Bosworth, The International Laser Ranging Service, Advances in Space Research 30 (2) (2002) 135–143.  doi:10.1016/S0273-1177(02)00277-6.
  * [51] A. Hauschild, O. Montenbruck, Precise real-time navigation of leo satellites using gnss broadcast ephemerides, NAVIGATION 68 (2) (2021) 419–432.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/navi.416, doi:https://doi.org/10.1002/navi.416.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/navi.416

  * [52] H. Hofmann-Wellenhof, E. Wasle, GNSS — Global Navigation Satellite Systems: GPS, GLONASS, Galileo, and more, Springer, 2008.
  * [53] L. Chen, X.-Z. Bai, Y.-G. Liang, K.-B. Li, Orbital Data Applications for Space Objects, 2017.  doi:10.1007/978-981-10-2963-9.
  * [54] P. Kovář, P. Puričer, K. Kovářová, Study of the two-line element accuracy by 1u cubesat with a gps receiver, Sensors 22 (8) (2022).  doi:10.3390/s22082902.
URL https://www.mdpi.com/1424-8220/22/8/2902

  * [55] C. Leake, H. Johnston, TFC: A Functional Interpolation Framework, https://github.com/leakec/tfc (2021).
  * [56] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, S. Wanderman-Milne, JAX: composable transformations of Python+NumPy programs (2018).
URL http://github.com/google/jax

  * [57] E. Hairer, S. Norsett, G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, Vol. 8, 1993.  doi:10.1007/978-3-540-78862-1.
  * [58] M. Abramowitz, I. A. Stegun, Handbook of Mathematical Functions, Dover Publications, 1972.

##  Appendix A Evaluation of the acceleration $\textit{g}_{10}$ for the point
$(x,y,z)$

The acceleration ${\bm{g}}_{N}$ up to the order $N$ for the point of
coordinates $(x,y,z)$ in the ITRF is given by

${\bm{g}}_{N}=\mathbf{\nabla}V=\frac{\partial V}{\partial
r}\mathbf{u}_{r}+\frac{1}{r}\frac{\partial
V}{\partial\theta}\mathbf{u}_{\theta}+\frac{1}{r\sin\theta}\frac{\partial
V}{\partial\lambda}\mathbf{u}_{\lambda},$ (5)

with the spherical coordinates $(r,\theta,\lambda)$ where $\theta$ is the
colatitude, the orthonormal spherical basis
$(\mathbf{u}_{r},\mathbf{u}_{\theta},\mathbf{u}_{\lambda})$, and the potential
$V$ given by

$V\left(r,\theta,\lambda\right)=\frac{GM}{r}\left(1+\sum^{N}_{n=2}\left(\frac{a}{r}\right)^{n}\sum^{n}_{m=0}\left(\overline{C}_{nm}\cos
m\lambda+\overline{S}_{nm}\sin
m\lambda\right)\overline{P}_{nm}\left(\tau\right)\right),$ (6)

for the order $n$, the degree $m$, the normalized spherical harmonic
coefficients of the gravitational potential
$(\overline{C}_{nm},\overline{S}_{nm})$, the normalized associated Legendre
polynomials $\overline{P}_{nm}$, and $\tau=\cos\theta$

By definition, the Legendre polynomials $P_{n}$ are given by [e.g. 58]

$P_{n}\left(\tau\right)=\frac{1}{2^{n}}\sum^{\left\lfloor\frac{n}{2}\right\rfloor}_{l=0}\left(-1\right)^{l}\left(\begin{matrix}n\\\
l\end{matrix}\right)\left(\begin{matrix}2n-2l\\\
n\end{matrix}\right)\tau^{n-2l},$ (7)

where $\lfloor\rfloor$ is the floor function, and the normalized associated
Legendre polynomials $\overline{P}_{nm}$ by

$\overline{P}_{nm}\left(\tau\right)=\frac{1}{\mathcal{N}_{nm}}\left(1-\tau^{2}\right)^{\frac{m}{2}}\frac{d^{m}}{d\tau^{m}}P_{n}\left(\tau\right),$
(8)

with

$\mathcal{N}_{nm}=\sqrt{\frac{\left(n+m\right)!}{\left(n-m\right)!\left(2n+1\right)\left(2-\delta_{0m}\right)}},$
(9)

and the Kronecher delta $\delta$. The normalized associated Legendre
polynomials $\overline{P}_{nm}$ can be then expressed as

$\overline{P}_{nm}\left(\tau\right)=\left(1-\tau^{2}\right)^{\frac{m}{2}}\sum^{\left\lfloor\frac{n-m}{2}\right\rfloor}_{j=0}\gamma_{nmf\left(j\right)}\tau^{f\left(j\right)},$
(10)

where

$\gamma_{nml}=\frac{\prod^{n}_{p=1}\left(n+m+l+1-2p\right)}{l!\left(n-l-m\right)!\mathcal{N}_{nm}},$
(11)

and

$f\left(j\right)=2j+\left\lfloor\frac{n+m+1}{2}\right\rfloor-\left\lfloor\frac{n+m}{2}\right\rfloor.$
(12)

The derivatives of the potential with the respect to the spherical coordinates
are then given with $\rho=r/a$ by

$\frac{\partial V}{\partial
r}=-\frac{GM}{a^{2}\rho^{2}}\Bigg{(}1+\sum^{N}_{n=2}\frac{1}{\rho^{n}}\sum^{n}_{m=0}\left(1-\tau^{2}\right)^{\frac{m}{2}}\sum^{\left\lfloor\frac{n-m}{2}\right\rfloor}_{j=0}\left(\alpha_{1nmf\left(j\right)}\cos
m\lambda+\beta_{1nmf\left(j\right)}\sin
m\lambda\right)\tau^{f\left(j\right)}\Bigg{)},$ (13) $\frac{\partial
V}{\partial\theta}=\frac{GM}{a\rho}\Bigg{(}\sum^{N}_{n=2}\frac{1}{\rho^{n}}\sum^{n}_{m=0}\left(1-\tau^{2}\right)^{\frac{m-1}{2}}\sum^{\left\lfloor\frac{n+1-m}{2}\right\rfloor}_{j=0}\left(\alpha_{2nmg\left(j\right)}\cos
m\lambda+\beta_{2nmg\left(j\right)}\sin
m\lambda\right)\tau^{g\left(j\right)}\Bigg{)},$ (14) $\frac{\partial
V}{\partial\lambda}=\frac{GM}{a\rho}\Bigg{(}\sum^{N}_{n=2}\frac{1}{\rho^{n}}\sum^{n}_{m=1}\left(1-\tau^{2}\right)^{\frac{m}{2}}\sum^{\left\lfloor\frac{n-m}{2}\right\rfloor}_{j=0}\left(\alpha_{3nmf\left(j\right)}\cos
m\lambda+\beta_{3nmf\left(j\right)}\sin
m\lambda\right)\tau^{f\left(j\right)}\Bigg{)},$ (15)

with

$g\left(j\right)=2j+\left\lfloor\frac{n+m}{2}\right\rfloor-\left\lfloor\frac{n+m-1}{2}\right\rfloor,$
(16)

and the coefficients

$\alpha_{1nml}=\gamma_{nml}\left(n+1\right)\overline{C}_{nm},\
\qquad\beta_{1nml}=\gamma_{nml}\left(n+1\right)\overline{S}_{nm},$ (17)
$\alpha_{2nml}=\epsilon_{nml}\overline{C}_{nm},\
\qquad\beta_{2nml}=\epsilon_{nml}\overline{S}_{nm},$ (18)
$\alpha_{3nml}=\gamma_{nml}m\overline{S}_{nm},\
\qquad\beta_{3nml}=-\gamma_{nml}m\overline{C}_{nm},$ (19)

with

$\epsilon_{nml}=\frac{\left(n+1-m\right)\left(n+m+l\right)-l\left(n+1\right)}{l!\left(n+1-l-m\right)!}\left(\prod^{n}_{p=1}\left(n+m+l-2p\right)\right).$
(20)