# Collision Avoidance Maneuvers Optimization in the Presence of Multiple
Encounters

Zeno Pavanello111PhD candidate, Te Pūnaha Ātea - Space Institute, 20 Symonds
Street, Auckland Central, Auckland 1010, New Zealand;
<EMAIL_ADDRESS>(Corresponding Author). Te Pūnaha Ātea - The Space
Institue, The University of Auckland, Auckland 1010, New Zealand  Laura
Pirovano222Research Fellow, Te Pūnaha Ātea - Space Institute, 20 Symonds
Street, Auckland Central, Auckland 1010, New Zealand;
<EMAIL_ADDRESS>and Roberto Armellin333Professor, Te Pūnaha
Ātea - Space Institute, 20 Symonds Street, Auckland Central, Auckland 1010,
New Zealand<EMAIL_ADDRESS>Member AIAA. Te Pūnaha Ātea -
The Space Institue, The University of Auckland, Auckland 1010, New Zealand
Andrea De Vittori 444Ph.D. candidate, Department of Aerospace Science and
Technology, Via La Masa 34, Milan 20156. Italy. and Pierluigi di Lizia
555Associate Professor, Department of Aerospace Science and Technology, Via La
Masa 34, Milan 20156. Italy. Polytechnic University of Milan, Milan 20156,
Italy

###### Abstract

The optimization of fuel-optimal low-thrust collision avoidance maneuvers
(CAMs) in scenarios involving multiple encounters between spacecraft is
addressed. The optimization’s objective is the minimization of the total fuel
consumption while respecting constraints on the total probability of
collision. The solution methodology combines sequential convex programming,
second-order cone programming, and differential algebra to approximate the
non-convex optimal control problem progressively. A Gaussian mixture model
method is used to propagate the initial covariance matrix of the secondary
spacecraft, allowing us to split it into multiple mixands that can be treated
as different objects. This leads to an accurate propagation of the
uncertainties. No theoretical guarantee is given for the convergence of the
method to the global optimum of the original optimal control problem.
Nonetheless, good performance is demonstrated through case studies involving
multiple short- and long-term encounters, showcasing the generation of fuel-
efficient CAMs while respecting operational constraints.††A part of this work
was presented at the 2024 AIAA Science and Technology Forum and Exposition,
held in Orlando, Florida, January 8-12 2024. Paper AIAA-2024-0845 was entitled
”A Convex Optimization Method for Multiple Encounters Collision Avoidance
Maneuvers”.

##  1 Introduction

Spacecraft operations are encountering growing challenges due to the
proliferation of space debris and the saturation of valuable orbital regimes,
notably the geostationary Earth orbit (GEO) and the low Earth orbit (LEO)
ones, where large constellations reside [1]. Ensuring the safety of spacecraft
has become a paramount concern. In this regard, the collision avoidance
maneuver (CAM)has emerged as a critical requirement for spacecraft operators,
aiming to autonomously and efficiently navigate through potentially hazardous
encounters. While the existing body of research extensively addresses single
short-term encounters [2, 3, 4, 5], recent attention has shifted toward
scenarios involving multiple close encounters [6, 7, 8, 9, 10]. This paper
presents an innovative approach to address the complex problem of fuel-
efficient CAMs in the context of multiple short- and long-term encounters. In
short-term encounters, the problem can be simplified, and each conjunction can
be framed in its B-plane, whereas this is not possible in the long-term case
because the path of relative motion inside the combined covariance ellipsoid
is far from rectilinear, and the encounter is not instantaneous [11].

Conventional strategies for CAM optimization primarily focus on mitigating a
single short-term encounter through the reduction of the probability of
collision (PoC) or the increase of the miss distance [2, 4, 12, 13, 14, 15,
5]. Very few strategies have been developed for long-term encounters [16, 17,
18]. However, the increasing density of space debris and the prospect of
multiple close encounters demand a more comprehensive and sophisticated
approach. Lately, NASA is taking steps in this direction, with the
introduction of the concept of total probability of collision (TPoC) [8] to
quantify the contribution of multiple conjunctions to the overall collision
probability: they recognize that in an environment that is densely populated,
treating each conjunction as a separate event will become unsustainable for
operators. The challenge lies in developing methodologies that can efficiently
compute optimal trajectories while considering both onboard propulsion’s
practical constraints and the encounters’ combined probabilistic nature.
Regarding long-term encounters, few publications have addressed the problem of
computing a CAM [16, 17, 18, 19]. In this circumstance, instead of minimizing
PoC, the alternative instantaneous probability of collision (IPoC) is
controlled and kept below a certain threshold for the whole window of
interest.

The multiple short-term problem has been faced in previous literature. Kim et
al. [6] developed a genetic algorithm approach that deals with up to four
conjunctions but must perform a separate $\Delta v$ impulse for each
encounter. Their method does not include a check on TPoC; instead, they
minimize PoC of each conjunction separately. This results in a higher final
TPoC, and thus the method underestimates the risk of a collision. The authors
do not provide a run-time for their method, but genetic algorithms are
notoriously slower to converge with respect to other methods. Sánchez and
Vasile [9] developed a comprehensive tool that is based on a linear
approximation of the dynamics with impulsive maneuvers and finds the optimal
CAM using a min-max optimization. Their method combines machine learning and
multi-criteria decision-making techniques to assess the collision risk and
rank possible optimal CAMs. They propose two main strategies: either computing
a single maneuver or treating each conjunction as a separate event. In the
first case, only the worst-case scenario is considered; therefore, the
computed impulse may not be optimal. In the second case, a separate maneuver
is computed for each one of the encounters: this means that after the maneuver
for the first conjunction has been computed, PoC must be re-calculated for
each of the following conjunctions. This process must be repeated as many
times as the number of the conjunctions. Sánchez and Vasile’s approach can be
very beneficial for operators on the ground to assist in the selection of the
preferable CAM for a specific scenario. Still, it cannot be used to compute
autonomous maneuvers due to the high computational load and the necessary
human interaction that it requires. Massòn et al. [10] cast the problem as a
linear program (LP) with non-convex quadratic constraints, and they solve it
by preallocating the maneuvering opportunities. They linearize the dynamics
using the Yamanaka-Ankersen state transition matrix (STM), and they impose
linear station keeping (SK) constraints. Since they employ the $l^{1}$ norm of
the control to keep the problem linear, they could recover a solution with
higher propellant consumption if compared with the $l^{2}$ norm. To summarize,
the state-of-the-art still lacks a method that can deal with an arbitrary
number of conjunctions comprehensively (using TPoC) and does not need a-priori
information on the maneuvering time. Moreover, all of the above-referenced
works considered Keplerian dynamics, disregarding the contributions of orbital
perturbations. Finally, to the best of the authors’ knowledge, the multiple
long-term conjunction CAM optimization problem has never been addressed.

An aspect that is closely related to this topic is the nonlinear propagation
of the uncertainty in the state of the objects involved in the conjunctions.
When the time span of the scenario involves multiple orbits and the
uncertainties of the secondary spacecraft are discretely large, the initial
Gaussian approximation can lose validity [20, 21]. We assume that accurate
orbit determination is performed on the primary so that its uncertainty can be
linearly propagated and remain Gaussian in the considered time window.
Different methods can be used to recover an accurate description of the
uncertainty of the secondary after a long propagation. Among these, the Monte
Carlo analysis is the most accurate but most computationally expensive;
unscented transforms reduce the computational load but assume Gaussian
distribution, Gaussian mixture model (GMM) provide a good compromise between
accuracy and computational load [22]. The research papers previously mentioned
either do not consider the evolution of the uncertainty (because they use the
instantaneous value of covariance given by a Conjunction Data Message (CDM)),
or they only consider a linear propagation. In this regard, a first attempt at
including the nonlinear propagation of the covariance matrix in the CAM
optimization was taken by Dutta and Misra, who recently developed a pseudo-
squared Mahalanobis distance (SMD) approach to treat a single conjunction with
GMM [23]. This work will show how the approach developed for multiple
conjunctions defined by CDMs can be adapted to treat multiple conjunctions
where the initial covariance is propagated nonlinearly using a GMM. This
allows for a more accurate depiction of the keep-out-zone (KOZ) when defining
the convexified PoC constraint.

The solution methodology is based on convex optimization methods and is
tailored to address the intricacies of multiple encounters and nonlinear
uncertainty propagation. This direct optimization approach is chosen for its
proven capability in solving complex aerospace engineering problems
efficiently [24, 25, 26]. These findings highlight the convenience of
employing a second-order cone program (SOCP) to formulate the convexified
problem, as it allows for the minimization of the $l^{2}$ norm of the control
history, which is a crucial requirement in spacecraft trajectory optimization
problems. Turning the problem into an SOCP involves convexifying all non-
convex aspects of the original optimal control problem (OCP), notably the
fuel-optimal objective function, the nonlinear dynamics constraint, and the
PoC constraint [27]. In the proposed approach, lossless relaxation is used to
linearize the objective function. Furthermore, the nonlinear dynamics are
automatically linearized using differential algebra (DA), which gives absolute
freedom on the dynamics model considered. The convexification approach
proposed in reference [5] is used to linearize the PoC constraint in the
short-term problem, and the one proposed in reference [18] is employed for the
IPoC constraint in the long-term problem.

To enhance convergence, a trust region constraint, drawing from references
[28, 29, 30], is employed to confine the solution space. However, due to the
linearization of the dynamics, the optimizer cannot attain the optimal
solution in a single SOCP run. Hence, a sequential convex program (SCP) is
constructed as a sequence of SOCPs, where the solution from the previous
iteration serves as the linearization point for the new iteration. This
iterative process may necessitate a discrete number of iterations before
converging to the optimal solution. Moreover, to account for the combined
nature of TPoC between the multiple encounters, a third outer iteration
process is used to find the optimal maneuver. This process automatically
adapts the single PoC limits associated with the conjunctions. Alternatively,
a single SCP run is performed where, in place of the method proposed in [5]
and [18], a direct linearization of the PoC or IPoC constraint is used.

In the subsequent sections, the paper will present the proposed methodology’s
mathematical framework and elaborate on the SCP approach for approximating the
optimal control problem. Simulated case studies will demonstrate the
approach’s effectiveness in generating fuel-efficient CAMs for multiple short-
and long-term encounters. By addressing this increasingly pertinent challenge,
this paper aims to offer a valuable contribution to the scope of spacecraft
collision avoidance (CA) in cluttered orbital environments.

##  2 Formulation of the Collision Avoidance Problem

We consider multiple conjunctions between a primary spacecraft and
$n_{conj}\in\mathbb{N}$ secondary objects. The primary is subject to a close
encounter with each secondary at a certain time of closest approach (TCA),
indicated with $t_{CA,s}$, where $s\in\\{1,...\hskip 2.0pt,n_{conj}\\}$. The
primary and the secondary objects are characterized by a mean state and a
covariance matrix at the starting time of the simulation $t_{0}$ or the
conjunction’s TCA, $t_{CA,s}$. The states are represented using any arbitrary
set of elements, e.g. Cartesian, Keplerian, or generalized equinoctial
elements

$\displaystyle\bm{x}_{p}(t_{0})\sim\mathcal{N}(\bm{\zeta}_{p}(t_{0}),\bm{C}_{p}(t_{0})),\quad$
$\displaystyle\bm{x}_{s}(t_{CA,s})\sim\mathcal{N}(\bm{\zeta}_{s}(t_{CA,s}),\bm{C}_{s}(t_{CA,s}))$
$\displaystyle\text{for }s\in\\{1,...\hskip 2.0pt,n_{conj}\\}$ (1)

where $\bm{x}_{p}(t_{0})$ and $\bm{x}_{s}(t_{CA,s})\in\mathbb{R}^{6}$ are
Gaussian multivariate random variable (MRV) variables, $\bm{\zeta}_{p}(t_{0})$
and $\bm{\zeta}_{s}(t_{CA,s})\in\mathbb{R}^{6}$ are the mean values,
$\bm{C}_{p}(t_{0})$ and $\bm{C}_{s}(t_{CA,s})\in\mathbb{R}^{6\times 6}$ are
the covariance matrices.

The propagation of the state of the primary and the secondaries can be
performed under a generic dynamical model

$\displaystyle\dot{\bm{x}}_{p}(t)=\bm{f}_{p}(t,\bm{x}_{p}(t),\bm{u}(t),\bm{p}_{p}),\quad$
$\displaystyle\dot{\bm{x}}_{s}(t)=\bm{f}_{s}(t,\bm{x}_{s}(t),\bm{p}_{s})$
$\displaystyle\text{for }s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (2)

where $t\in\mathbb{R}_{[t_{0},t_{f}]}$ is the independent time variable,
$\bm{x}_{p}(t),\bm{x}_{s}(t),\dot{\bm{x}}_{p}(t)$ and
$\dot{\bm{x}}_{s}(t)\in\mathbb{R}^{6}$ are the time-dependent states of the
primary and the secondaries and their derivatives,
$\bm{u}(t)\in\mathbb{R}^{3}$ is the control action,
$\bm{p}_{p}\in\mathbb{R}^{m_{p}}$ and $\bm{p}_{s}\in\mathbb{R}^{m_{s}}$ are
vectors of parameters,
$\bm{f}_{p}(\cdot):\mathbb{R}_{[t_{0},t_{f}]}\times\mathbb{R}^{6}\times\mathbb{R}^{3}\times\mathbb{R}^{m_{p}}\rightarrow\mathbb{R}^{6}$
and
$\bm{f}_{s}(\cdot):\mathbb{R}_{[t_{0},t_{f}]}\times\mathbb{R}^{6}\times\mathbb{R}^{m_{s}}\rightarrow\mathbb{R}^{6}$
are continuous functions..

To precisely establish the PoC constraint, it is crucial to explicitly express
the relative position of the primary object w.r.t. to each secondary. If the
state of the two objects is expressed in Cartesian Earth Centered Inertial
(ECI) coordinates ($\vec{x}_{p}(t)=[\bm{r}_{p}(t);\hskip
4.0pt\vec{v}_{p}(t)]$), as it is the case in Section 4 the relative position
is the subtraction of the first three elements of the two normally distributed
MRV

$\displaystyle\Delta\bm{r}_{s}(t)=\bm{r}_{p}(t)-\bm{r}_{s}(t)\quad$
$\displaystyle s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (3)

Given that the subtraction is a linear transformation, the relative position
is also normally distributed

$\displaystyle\Delta\bm{r}_{s}(t)\sim\mathcal{N}\big{(}\bm{\mu}_{p}(t)-\bm{\mu}_{s}(t),\bm{P}_{p}(t)+\bm{P}_{s}(t)\big{)}\quad$
$\displaystyle s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (4)

The standard OCP in the continuous domain for the short-term conjunction CAM
problem is stated as follows

$\displaystyle\min_{u}\quad$ $\displaystyle
J=\int_{t_{0}}^{t_{f}}u(t)\mathrm{d}t,$ (5a) s.t.
$\displaystyle\dot{\bm{x}}_{p}=\bm{f}(t,\bm{x}_{p}(t),\bm{u}(t)),$ (5b)
$\displaystyle P_{TC}(t)\leq\bar{P}_{TC},$ (5c)
$\displaystyle\bm{x}(t_{0})=\bm{x}_{0},$ (5d) $\displaystyle
u(t)=\sqrt{u_{1}(t)^{2}+u_{2}(t)^{2}+u_{3}(t)^{2}},$ (5e) $\displaystyle
u(t)\leq u_{max},$ (5f)

where $P_{TC}\in\mathbb{R}(\cdot):\mathbb{R}\rightarrow\mathbb{R}$ is the
TPoC, which is monotonically increasing with time, and
$\bar{P}_{TC}\in\mathbb{R}$ is the limit value imposed to the TPoC. In Problem
(5), Eq. 5a is the fuel minimization objective function, Eq. 5b is the
dynamics constraint, Eq. 5c is the PoC constraint, Eq. 5d is the initial state
bound, Eq. 5e is an auxiliary constraint to define the norm of the control,
and Eq. 5f is the bound on the maximum value of the control action. The
collision probability can include $n_{conj}$ conjunctions. Since CAMs
typically involve a $\Delta v$ on the order of magnitude of
[$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$] or
[$\mathrm{c}\mathrm{m}\mathrm{/}\mathrm{s}$], the mass loss due to the
maneuver is not considered in the equations of motion and the optimization
problem [5].

###  2.1 Multiple Short-Term Conjunctions

Conjunctions can be categorized into two types: short-term and long-term.
Short-term conjunctions, being more frequent, have received a more
comprehensive investigation. In such instances, the relative velocity is
notably high, resulting in near-instantaneous collisions. Consequently, the
dynamics of the conjunction can be effectively approximated as linear without
sacrificing accuracy, and the event is typically studied on the two-
dimensional B-plane [31, 5, 4, 32].

In a scenario where $n_{conj}$ consecutive conjunctions take place, the
B-plane of a conjunction $s$ is centered on the secondary object, the $\eta$
axis is aligned with the direction of the relative velocity of the primary
w.r.t. the secondary, and the $\xi\zeta$ plane is perpendicular to the
relative velocity axis. Since TCA is the time when the miss distance is
lowest, it follows that
$\Delta\bm{r}_{s}(t_{CA,s})\cdot\Delta\bm{v}_{s}(t_{CA,s})=0$ and thus the
relative position lies on the $\xi\zeta$ plane. The following discussion
assumes that all the uncertainty is concentrated around the secondary object,
and all the mass is concentrated around the primary [31, 33]. In this way, the
primary’s state, which has to be optimized, can be treated as a deterministic
variable.

The PoC of a single encounter is given by the integration of the Probability
Density Function (PDF) of the relative position over the circle
$\mathbb{C}_{\mathrm{HBR}}$ defined by the combined hard body radius (HBR) of
the two spacecraft and centered in the primary

$P_{C,s}=\frac{1}{(2\pi)\sqrt{\mathrm{det}\big{(}\bm{P}_{\mathcal{B}}(t_{CA,s})\big{)}}}\iint_{\mathbb{C}_{\mathrm{HBR}}}\mathrm{exp}\left(-\frac{\Delta\bm{r}_{\mathcal{B}}^{\mathrm{T}}(t_{CA,s})\bm{P}_{\mathcal{B}}^{-1}(t_{CA,s})\Delta\bm{r}_{\mathcal{B}}(t_{CA,s})}{2}\right)\mathrm{d}A,$
(6)

where $\Delta\bm{r}_{\mathcal{B}}(t_{CA,s})\in\mathbb{R}^{2}$
$\bm{P}_{\mathcal{B}}(t_{CA,s})\in\mathbb{R}^{2\times 2}$ are the projection
of the relative position mean value and covariance on the encounter B-plane at
TCA. The argument of the exponential function is the SMD at TCA,
$d_{m,s}^{2}$, divided by $2$. Various approaches have been proposed to
numerically solve this integral over the years; in this work, we will use
Chan’s method [34]. By employing a numerical inversion of Chan’s formula [32],
it is possible to get the SMD limit from a PoC limit for the conjunction
($\bar{d}_{m}^{2}=f(\bar{P}_{C})$). Chan’s PoC formula is generally preferred
over others because of its accuracy and computational efficienty [35].

Assuming that the conjunctions are uncorrelated, the TPoC of the multiple
encounter scenario is computed as the complement of the probability of not
colliding with any object

$P_{TC}=1-\prod_{s=1}^{n_{conj}}\big{(}1-P_{C,s}\big{)},$ (7a)
$P_{TC}\approx\sum_{s=1}^{n_{conj}}P_{C,s}.$ (7b)

In the hypothesis of small probabilities, it can be shown that the product in
Eq. 7a can be replaced by the summation in Eq. 7b. We will always use Eq. 7a
when computing the TPoC, but we will leverage the approximated formula Eq. 7b
to compute first guesses for the limits values of the PoC of the single
conjunctions. One might assume that the hypothesis of uncorrelation is a
strong one because the trajectory of the primary would be affected by a
collision, thus changing the PoC of the subsequent conjunctions. Nonetheless,
as Frigm et al. highlight [8], from an operation point of view, "there is no
reason to suppose that a close approach with a particular secondary is likely
to promote a close approach with some other secondary."

###  2.2 Multiple Long-Term Conjunctions

Patera [11] demonstrates that when the relative trajectory of the primary
spacecraft diverges significantly from a straight line within the combined
covariance ellipsoid, the encounter qualifies as long-term, requiring the use
of a distinct model for the computation of PoC. Despite numerous recent
efforts to develop efficient methods for computing PoC in long-term encounters
(references [36, 37, 38, 39]), the presence of nonlinear relative trajectories
within the combined uncertainty ellipsoid makes the risk quantification and,
consequently, the design of long-term CAMs very complex. Quantitative
thresholds to distinguish between the two types of conjunctions have been
proposed. Among these, Chan’s [40] states that conjunctions are short-term if
the path of relative motion may be considered a straight line over a distance
of $8-25$ $\mathrm{k}\mathrm{m}$ with a deviation lower than $2$ $\mathrm{m}$;
Dolado et al.’s [41], instead, identify conjunctions as short-term if the
relative velocity at TCA is higher than $10$ $\mathrm{m}\mathrm{/}\mathrm{s}$.
In recent literature, it has been shown that IPoC is a valuable alternative to
PoC when addressing long term encounters [42, 19, 43, 44, 33]: in general, an
increase in PoC is accompanied by a high value of IPoC, so controlling the
evolution of IPoC allows to indirectly control the growth of PoC. For this
reason, in agreement with a good portion of the state-of-the-art, we use IPoC
as a risk metric rather than the more complex PoC when dealing with long-term
encounters.

The integral of the IPoC at any time $t\in\mathbb{R}_{[t_{0},t_{f}]}$ must be
performed over the the three spatial dimensions

$P_{IC,s}(t)=\frac{1}{(2\pi)^{3/2}\mathrm{det}(\bm{P}(t))^{1/2}}\iiint_{\mathbb{S}_{\mathrm{HBR}}}\mathrm{exp}\left(-\frac{\Delta\bm{r}_{s}^{\mathrm{T}}(t)\bm{P}_{rel,s}^{-1}(t)\Delta\bm{r}_{s}(t)}{2}\right)\mathrm{d}V,$
(8)

where both the relative position mean value
$\Delta\bm{r}_{s}(t)\in\mathbb{R}^{3}$ and its covariance
$\bm{P}_{rel,s}(t)\in\mathbb{R}^{3\times 3}$ evolve in time. Multiple
approaches have been proposed to solve this integral. We will consider one
analogous to Alfriend’s [31] for long-term encounters (its derivation can be
found in [18]). Under the assumption of this method, i.e., the PDF of the
relative position is considered constant inside the hard body sphere, IPoC has
a direct dependence on SMD

$P_{IC,s}(t)=\sqrt{\frac{2}{\pi\mathrm{det}\big{(}\bm{P}(t)\big{)}}}\frac{R^{3}}{3}\mathrm{exp}\left(-\frac{d_{m,s}^{2}(t)}{2}\right).$
(9)

In this case, the SMD is computed using the three-dimensional relative
position variable. Analogously to Eq. 7a, the total instantaneous probability
of collision (TIPoC) in a multiple long-term encounters scenario is given by
the combination of the single uncorrelated conjunctions

$P_{TIC}(t)=1-\prod_{s=1}^{n_{conj}}\big{(}1-P_{IC,s}(t)\big{)}\approx\sum_{s=1}^{n_{conj}}P_{IC,s}(t),$
(10)

Where the approximated equation is obtained using the same approximation of
Eq. 7b.

###  2.3 Gaussian Mixture Model for Uncertainty Propagation

When dealing with multiple encounters with the same secondary or with a long-
term conjunction, the time frame of interest can be very extended, spanning in
the order of magnitude of the orbital periods. In the former case, the B-plane
configuration depends on the propagated covariance, whereas in the latter, it
is necessary to know the values of the elements of the covariance matrix of
the relative position at every time to compute the IPoC from Eq. 9. If the
propagation period is sufficiently long, the nonlinearities of the dynamics
propagation become predominant, and the Gaussianity of the state is lost.

An effective way to recover a precise propagation of the uncertainty, which
avails in the computation of the collision metrics, is the use of a GMM.
Through this technique, the PDF of the secondary spacecraft is split into
$n_{mix}\in 2\mathbb{N}+1$ mixands. Each mixand is characterized by a mean
value and a covariance matrix derived from the original distribution. The
propagations of the mixands happen independently, and they are handled as
different secondaries. The initial PDF split of the secondary’s state yields

$p[\bm{x}(t)]=\sum_{c=1}^{n_{mix}}\gamma_{c}p_{g}\left[\bm{x}(t);\bm{\mu}_{c}(t),\bm{P}_{c}(t)\right],$
(11)

where $c\in\\{1,...\hskip 2.0pt,n_{mix}\\}$ indicates the single mixand, $p$
is the PDF of the secondary object, $p_{g}$ indicates the Gaussian PDF, and
$\gamma_{c}\in\mathbb{R}$ are a set of weights such that
$\sum_{c=1}^{n_{mix}}\gamma_{c}=1$. It is assumed that most of the uncertainty
is concentrated in the state of the secondary, so there is no need to split
the state of the primary.

Vittaldev’s uni-variate library method is used to perform the split [22].
Given the original PDF $p[\bm{x}(t)]$, the uni-variate splitting library is
applied along the direction $\hat{\bm{a}}$

$\bm{\mu}_{c}=\bm{\mu}+\mu_{c}\bm{S}\hat{\bm{a}}^{\star},$ (12a)
$\bm{P}_{c}=\bm{S}\left[\bm{I}+\left(\sigma_{c}^{2}-1\right)\hat{\bm{a}}^{\star}\hat{\bm{a}}^{\star
T}\right]\bm{S}^{T},$ (12b)

where the matrix $\bm{S}$ is the Cholesky decomposition of $\bm{P}$, $\mu_{c}$
and $\sigma_{c}$ are from the uni-variate splitting library and
$\hat{\bm{a}}^{\star}$ is

$\hat{\bm{a}}^{\star}=\frac{S^{-1}\hat{\bm{a}}}{\left\|S^{-1}\hat{\bm{a}}\right\|_{2}}.$
(13)

Once the split has been performed, $n_{mix}$ fictitious secondary objects are
defined at TCA, which can be treated independently in the CA framework. Each
of these objects encounters the primary spacecraft each of the conjunction
times $n_{conj}$. The TPoC or TIPoC, then, is derived from Eq. 7a or Eq. 10
with the introduction of a nested product function and the mixand weights

$P_{TC}=1-\prod_{s=1}^{n_{conj}}\prod_{c=1}^{n_{mix}}\left(1-\gamma_{c}P_{C,cs}\right))\approx\sum_{s=1}^{n_{conj}}\sum_{c=1}^{n_{mix}}\left(\gamma_{c}P_{C,cs}\right),$
(14a)
$P_{TIC}(t)=1-\prod_{s=1}^{n_{conj}}\prod_{c=1}^{n_{mix}}\left(1-\gamma_{c}P_{IC,cs}(t)\right)\approx\sum_{s=1}^{n_{conj}}\sum_{c=1}^{n_{mix}}\left(\gamma_{c}P_{IC,cs}(t)\right),$
(14b)

where $P_{C,cs}$ and $P_{IC,cs}\in\mathbb{R}$ are the PoC and IPoC of the
$c$-th mixand in the $s$-th conjunction. The approximated equations are
obtained using the same approximation of Eq. 7b. The splitting direction
$\hat{\bm{a}}$ is arbitrary, but it determines how good an approximation the
GMM provides w.r.t. the real distribution. Following works like Losacco’s [29]
and Vittaldev’s [45], a good selection of the splitting direction depends both
on the nonlinearity index (NLI) of the propagation and on the components of
the initial covariance matrix. The NLI, indicated with $\bm{\nu}$, gives a
quantitative indication of the effect of the nonlinearities in a certain
direction; in reference [29], how to compute it is effectively shown. A vector
$\bm{\phi}\in\mathbb{R}^{6}$ made of the 2-norm of each column of the initial
covariance matrix is used to quantify the entity of the uncertainty in each
direction

$\displaystyle\phi_{k}=\sqrt{\sum_{l=1}^{6}P_{lk}^{2}}$ $\displaystyle
k\in[1,6],$ (15)

The splitting direction, then, considers the information from both the NLI and
the covariance

$\hat{\bm{a}}=\frac{\vec{\nu}\circ\bm{\phi}}{||\vec{\nu}||\cdot||\bm{\phi}||},$
(16)

where $\circ$ indicates the Hadamard product.

##  3 Convex Formulation

The OCP Problem 5 is transformed into a SCP to be solved iteratively. The
methods applied are described in detail in [19]: the SCP iteratively
approximates the global solution of the OCP via a sequence of SOCPs. Here, we
recall the constraints and the objective function and adapt it to the multiple
encounters problem. Three sources of non-convexities are present in the
original OCP as formulated in Problem 5:

  1. 1.

orbital dynamics in Eq. 5b are highly nonlinear, thus non-convex.

  2. 2.

Eq. 5c is a nonlinear constraint because it involves the use of Eq. 7a and the
inversion of Chan’s formula for short-term encounters or Eq. 9 for long-term
ones.

  3. 3.

Eq. 5a and Eq. 5e yield a nonlinear objective function.

These issues are addressed in the following sections, where the dynamics
constraint, the PoC constraint, and the objective function are convexified
using different techniques. A more detailed explanation of the convexification
approach can be found in [5] and [19].

###  3.1 Convexification of the dynamics

The dynamics of the problem are first discretized and then automatically
linearized using DA.

####  3.1.1 Discretization of the Dynamics

The continuous time variable $t\in\mathbb{R}_{[t_{0},t_{f}]}$ is substituted
by the discrete time variable $t_{i}\in\\{t_{0},t_{1},...,t_{N}\\}$, where
$N+1$ is the number of nodes of the discretization. Following Eq. 2 and via
the use of an integration scheme, like Runge-Kutta 7-8, one obtains the states
of the two spacecraft at node $i+1$, which depend on the state and the control
at node $i$:

$\displaystyle\bm{x}_{p,i+1}=\bm{f}_{p,i}(t_{i},\bm{x}_{p,i},\bm{u}_{i},\bm{p}_{p})\quad$
$\displaystyle i\in\\{0,...\hskip 2.0pt,N-1\\},\quad$ (17a)
$\displaystyle\bm{x}_{s,i+1}=\bm{f}_{s,i}(t_{i},\bm{x}_{s,i},\bm{p}_{s})\quad$
$\displaystyle i\in\\{0,...\hskip 2.0pt,N-1\\}\quad$ $\displaystyle
s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (17b)

where
$\bm{f}_{p,i}(\cdot):\mathbb{R}_{[t_{0},t_{f}]}\times\mathbb{R}^{6}\times\mathbb{R}^{3}\times\mathbb{R}^{m_{p}}\rightarrow\mathbb{R}^{6}$
and
$\bm{f}_{s,i}(\cdot):\mathbb{R}_{[t_{0},t_{f}]}\times\mathbb{R}^{6}\times\mathbb{R}^{m_{s}}\rightarrow\mathbb{R}^{6}$
are the functions that describe the dynamics at node $i$,
$\bm{x}_{p,i}=\bm{x}_{p}(t_{i})$, $\bm{x}_{s,i}=\bm{x}_{s}(t_{i})$ and
$\bm{u}_{i}=\bm{u}(t_{i})$. The DA tool is used to introduce perturbations on
the primary state ($\bm{x}_{p,i}+\delta\bm{x}_{p,i}$) and acceleration
($\bm{u}_{i}+\delta\bm{u}_{i}$) at each node, effectively expressing Eq. 17a
through Taylor polynomials:

$\displaystyle\bm{x}_{p,i+1}=\mathcal{T}^{q}_{\bm{x}_{p,i+1}}(\bm{x}_{p,i},\bm{u}_{i})\quad$
$\displaystyle i\in\\{0,...\hskip 2.0pt,N-1\\},$ (18)

where in general the expression $\mathcal{T}^{q}_{y}(x)$ indicates the
$q^{\text{th}}$-order Taylor expansion of the variable $y$ as a function of
$x$, around the expansion point $\tilde{x}$ in which the polynomial is
computed. The reader can find a detailed explanation of the use of DA in [46].

####  3.1.2 Linearization of the Dynamics

The discretized dynamics are linearized around a reference trajectory, so that
a linear constraint can be set in the SOCP. This is an iterative process which
can require numerous runs, referred to as major iterations and denoted by
index $j$. Once the optimization problem of a major iteration $j-1$ is solved,
the solution ($\mathbb{x}^{j-1}$, $\mathbb{u}^{j-1}$) becomes available, a
column vector comprising the state and control at each node. Employing DA,
this solution serves as an expansion point for constructing linear dynamics
maps for iteration $j$. The continuity condition is enforced by requiring that
the state after the propagation of node $i$ is equal to the state before the
propagation of node $i+1$. The state and acceleration of each node are
expanded around the reference points, which are the output of the previous
major iteration $\tilde{\vec{x}}_{i}^{j}=\bm{x}_{i}^{j-1}$ and
$\tilde{\vec{u}}_{i}^{j}=\bm{u}_{i}^{j-1}$

$\displaystyle\bm{x}_{i}^{j}=\tilde{\vec{x}}_{i}^{j}+\delta\bm{x}_{i}^{j},\quad$
$\displaystyle\bm{u}_{i}^{j}=\tilde{\vec{u}}_{i}^{j}+\delta\bm{u}_{i}^{j}$
$\displaystyle i\in\\{0,...\hskip 2.0pt,N-1\\},$ (19)

The state at the subsequent node is obtained through the propagation of the
first-order dynamics. The linear continuity constraint imposes that the
linearly propagated state be equal to the state at the subsequent node

$\displaystyle\bm{x}_{i+1}^{j}=\bm{A}_{i+1}^{j}\bm{x}_{i}^{j}+\bm{B}_{i+1}^{j}\bm{u}_{i}^{j}+\bm{c}_{i}^{j}\quad$
$\displaystyle i\in\\{0,...\hskip 2.0pt,N-1\\},$ (20)

where $\bm{A}_{i+1}^{j}\in\mathbb{R}^{6\times 6}$ is the state transition
matrix, $\bm{B}_{i+1}^{j}\in\mathbb{R}^{6\times 3}$ is the control-state
transition matrix,
$\bm{c}_{i}^{j}=\bar{\bm{x}}_{i+1}-\bm{A}_{i+1}^{j}\tilde{\bm{x}}_{i}^{j}-\bm{B}_{i+1}^{j}\tilde{\bm{u}}_{i}^{j}$
is the residual of the linearization, and
$\bar{\bm{x}}_{i+1}=\bm{f}_{p,i}(t_{i},\tilde{\bm{x}}_{i},\tilde{\bm{u}}_{i})$
is the constant part of the propagation of the state. The initial condition is
fixed because the maneuver cannot alter it

$\bm{x}_{0}^{j}=\bm{x}_{0}^{0}.$ (21)

####  3.1.3 Lossless Relaxation of the Control Magnitude Constraint

Equation Eq. 5e is a non-convex equality constraint. Following the work from
[47], we introduce a lossless relaxation to convexify the constraint: Eq. 5e
is transformed into an inequality constraint, and the control magnitude is
added to the optimization vector. In this way, Eq. 5e becomes a second-order
cone constraint. The variable $u_{i}$ is now allowed to take values higher
than the norm of the control that acts on the dynamics

$u_{i}\geq\sqrt{u_{i,1}^{2}+u_{i,2}^{2}+u_{i,3}^{2}}\quad i\in\\{0,...\hskip
2.0pt,N-1\\}.$ (22)

The discretized forms of Eq. 5a and Eq. 5f become respectively

$J=\sum_{i=0}^{N-1}u_{i}\cdot\Delta t_{i},$ (23)

$\displaystyle 0\leq u_{i}\leq 1$ $\displaystyle i\in\\{0,...\hskip
2.0pt,N-1\\},$ (24)

where $\Delta t_{i}=t_{i+1}-t_{i}$ is the time for which the control is active
from node $i$ to $i+1$. This linearization of the objective function is
lossless, meaning that the optimal solution for the convexified problem is
also optimal for the original problem.

####  3.1.4 Trust region approach

Using a trust region constraint in formulating the SOCP is greatly beneficial
for the convergence of the algorithm [28]. In this work, the approach
originally introduced in [30] and adapted in [19] is used, which is reliant on
the definition of the NLI from [29]. The interested reader is invited to refer
to these works for details on the derivation of the constraint. The gist of
the approach is to consider the radius of the trust region of an optimization
variable as inversely proportional to its nonlinearity. The resulting
constraints are

$\displaystyle\bm{\xi}_{i}\circ[\bm{x}_{i}^{\mathrm{T}};\hskip
2.0pt\bm{u}_{i}^{\mathrm{T}}]^{\mathrm{T}}\preceq\bm{\xi}_{i}\circ[\tilde{\bm{x}}_{i}^{\mathrm{T}};\hskip
2.0pt\tilde{\bm{u}}_{i}^{\mathrm{T}}]^{\mathrm{T}}+\bar{\nu}\cdot\mathbf{1}\quad$
$\displaystyle i\in\\{1,...\hskip 2.0pt,N\\},$ (25a)
$\displaystyle\bm{\xi}_{i}\circ[\bm{x}_{i}^{\mathrm{T}};\hskip
2.0pt\bm{u}_{i}^{\mathrm{T}}]^{\mathrm{T}}\succeq\bm{\xi}_{i}\circ[\tilde{\bm{x}}_{i}^{\mathrm{T}};\hskip
2.0pt\tilde{\bm{u}}_{i}^{\mathrm{T}}]^{\mathrm{T}}-\bar{\nu}\cdot\mathbf{1}\quad$
$\displaystyle i\in\\{1,...\hskip 2.0pt,N\\},$ (25b)

where $\bm{\xi}_{i}$ is the measure of the nonlinearity of the associated
state vector in the second-order DA propagation of the dynamics [30];
$\bar{\nu}$ is a user-defined value that imposes a limit to the maximum NLI
available to the solution.

The introduction of the trust region constraint can cause artificial
infeasibility, so virtual controls are added to the dynamics constraints, Eq.
20. The new constraints with virtual controls become

$\displaystyle\bm{x}_{i+1}^{j}-\bm{A}_{i+1}\bm{x}_{i}^{j}-\bm{B}_{i+1}^{j}\bm{u}_{i}^{j}+\vec{\upsilon}_{i+1}^{j}-\bm{c}_{i}^{j}=0\quad$
$\displaystyle i\in\\{0,...\hskip 2.0pt,N-1\\},$ (26)

where $\vec{\upsilon}_{i}^{j}\in\mathbb{R}^{6}$ is the virtual control vector
for the node $i$ at major iteration $j$.

In the objective function, a term is added which is proportional to the entity
of the virtual controls in order to minimize them. This term is multiplied by
a weight $\kappa_{vc}\in\mathbb{R}$, which must be sufficiently high.

$\displaystyle\upsilon_{i}\geq||(\upsilon_{ca,i}^{j})||\quad$ $\displaystyle
i\in\\{1,...\hskip 2.0pt,N\\},$ (27a) $\displaystyle
J_{vc}=\kappa_{vc}\sum_{i=0}^{N}\upsilon_{i},$ (27b)

###  3.2 Convexifcation of the PoC Constraint for Short-Term Encounters

Two alternative approaches are proposed to convexify the PoC constraint. The
first one is based on the linearization of the SMD, and it consists of setting
a separate constraint for each conjunction; the second is based on a direct
linearization of the TPoC formula of Eq. 7a.

####  3.2.1 Linearized SMD constraint

The PoC constraint for a single conjunction can be formulated by leveraging
the methodology proposed in references [48] and [5]. The reader is invited to
refer to this last reference for an in-depth explanation of the process. An
iterative projection and linearization algorithm is utilized to convexify the
nonlinear PoC constraint. The iterations are nested inside the major
iterations used to linearize the dynamics and take the name of minor
iterations, denoted by the symbol $k$. For each conjunction node, the
projection convex sub-problem aims to find the point on the surface of an
ellipse that is closest to the relative position $\Delta\bm{r}_{s}^{j,k-1}$
from the previous minor iteration (where $s\in\\{1,...\hskip
2.0pt,n_{conj}\\}$). In the formulation of the convex sub-problem, the indices
$s$ and $j$ are dropped since one of these problems is solved multiple times
inside the same major iteration and for each conjunction.

Let the covariance and relative position be expressed in the B-plane reference
frame $\mathcal{B}$. It is possible to define a transformation matrix
$\bm{V}\in\mathbb{R}^{2}$ that diagonalizes the covariance matrix. We call the
reference frame in which the covariance matrix is diagonal $\mathcal{C}$. A
simple quadratic optimization problem is used to find the point on the ellipse
that is closest to the reference trajectory point

$\displaystyle\min_{\hat{\bm{z}}_{\mathcal{C}}}\quad$
$\displaystyle||\bm{z}_{\mathcal{C}}^{k}-\Delta\bm{r}_{q,\mathcal{C}}^{k-1}||$
(28a) s.t.
$\displaystyle(\bm{z}_{\mathcal{C}}^{k})^{\mathrm{T}}(\bm{P}_{\mathcal{C}})^{-1}\bm{z}_{\mathcal{C}}^{k}\leq\bar{d}_{m}^{2},$
(28b)

The objective Eq. 28a imposes the minimization of the Euclidean distance
between the relative position of the previous iteration and the optimization
variable $\bm{z}_{\mathcal{C}}$. Eq. 28b is a relaxed condition on the
optimization variable to be inside the ellipsoid. These relaxed conditions are
lossless since the minimization of the objective guarantees that the optimized
variable is always positioned on the ellipse’s surface, maximizing the
distance from its center.

Once $\bm{z}_{\mathcal{C}}$ is determined, the solution is transformed back
into the original reference frame using the equation
$\bm{z}_{\mathcal{B}}^{k}=\bm{V}^{\mathrm{T}}\bm{z}_{\mathcal{C}}^{k}$. After
obtaining $\bm{z}_{s}^{j,k}=\bm{z}_{\mathcal{B}}^{k}$, a linearization of the
SMD constraint is applied, effectively turning it into a KOZ constraint.
Specifically, we ensure that the new optimized relative position lies within
the semi-plane defined by the line tangent to the ellipse on
$\bm{z}_{s}^{j,k}$. The equation of the constraint is, then

$\displaystyle\nabla(d_{m}^{2})_{s}^{j,k}\Big{|}_{\bm{z}_{s}^{j,k}}\cdot(\Delta\bm{r}_{i}^{j,k}-\bm{z}_{s}^{j,k})\geq
0$ $\displaystyle i\in\\{1,...\hskip 2.0pt,N\\}.$ (29)

####  3.2.2 Linearized TPoC constraint

An alternative to the use of the linearized SMD constraint is the linearized
TPoC constraint. The formula of the TPoC from Eq. 7a is linearized using DA as
a function of the relative position in each of the encounters

$\delta P_{TC}=\bm{\nabla}\delta\mathbb{r},$ (30)

where $\bm{\nabla}$ is the gradient of the expression that gives the total PoC
as a function of the relative position at each TCA and $\mathbb{r}$ is the
stacked vector of the positions of the primary in the nodes in which the
conjunctions happen, i.e., the TCAs

$\mathbb{r}=\begin{bmatrix}\bm{r}_{1}^{\mathrm{T}}&...&\bm{r}_{n_{conj}}^{\mathrm{T}}.\end{bmatrix}^{\mathrm{T}}$
(31)

The PoC constraint, in this case, is a scalar one

$P_{TC}=\bm{\nabla}\mathbb{r}+d\leq\bar{P}_{TC},$ (32)

where $d=\tilde{P}_{TC}-\bm{\nabla}\tilde{\mathbb{r}}$ is the residual of the
linearization (the tilde indicates the values computed at the expansion
point). Note that using the absolute or the relative position of the primary
is equivalent because in Eq. 32, the position of the secondary would be
eliminated. A trust region constraint is used to favor the convergence of the
optimization with this constraint.

The trust region method is analogous to the one used to impose a trust region
constraint based on the dynamics in [19]. Eq. 7a is computed using a second-
order DA expansion, which yields the gradient of the function with a first
order-expansion

$\bm{\nabla}=\bar{\bm{\nabla}}+\delta\bm{\nabla},$ (33)

where $\bar{\bm{\nabla}}\in\mathbb{R}^{3n_{conj}}$ is the constant part of the
gradient and $\delta\bm{\nabla}\in\mathbb{R}^{3n_{conj}}$ is its first-order
expansion w.r.t. the position vector.

$\displaystyle\nabla_{i}=\bar{\nabla}_{i}+\delta\nabla_{i}=\bar{\nabla}_{i}+\sum_{j=1}^{3n_{conj}}a_{ij}\delta
r_{j}\quad$ $\displaystyle i\in\\{1,...\hskip 2.0pt,3n_{conj}\\}$ (34)

The non-linearity vector for the single relative position variable $r_{j}$ is
computed as the ratio between the norm of the second-order terms and that of
the first-order ones

$\nu_{j}=\frac{\sqrt{\sum_{i=1}^{3n_{conj}}a_{ij}^{2}}}{||\bar{\bm{\nabla}}||}|\delta
r_{j}|=\xi_{j}|\delta r_{j}|.$ (35)

It is possible to write two vector equations for each conjunction, which bound
the position of the given node

$\displaystyle\bm{\xi}_{k}\circ\Delta\bm{r}_{k}\preceq\bm{\xi}_{k}\circ\Delta\tilde{\bm{r}}_{k}+\bm{\nu}\quad$
$\displaystyle k\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (36a)
$\displaystyle\bm{\xi}_{k}\circ\Delta\bm{r}_{k}\succeq\bm{\xi}_{k}\circ\Delta\tilde{\bm{r}}_{k}-\bm{\nu}\quad$
$\displaystyle k\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (36b)

where $\bm{\xi}_{k}$ and $\bm{\nu}_{k}\in\mathbb{R}^{3}$ include the three
components of the single encounter.

####  3.2.3 Refinement of the solution

The PoC constraint of the original OCP Eq. 5c bound the TPoC to below a limit
value $\bar{P}_{TC}$. This limit can be split, and a portion of it can be
assigned to the single conjunctions

$\bar{P}_{TC}=1-\prod_{s=1}^{n_{conj}}(1-\bar{P}_{C,s}).$ (37)

Theoretically, in Eq. 37, the single limits $\bar{P}_{C,s}$ can be optimized
to build constraints that allow for the minimization of the $\Delta v$. Under
the assumption that a higher $\Delta v$ leads to a lower TPoC, we need to find
a method to recover a solution for which the final TPoC gets as close as
possible to the total limit.

A first SCP is solved using a first guess for the limits of the single
conjunctions. They are all set to the same value by leveraging the
approximated Eq. 7b

$\displaystyle\bar{P}_{C,s}=\frac{\bar{P}_{TC}}{n_{conj}}\quad$ $\displaystyle
s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (38)

so that Eq. 37 is automatically satisfied. This first convex run’s total PoC
is typically much lower than $\bar{P}_{TC}$, so a refinement is executed. We
have two alternatives for the PoC constraint of this second SCP, which are
reported schematically in Fig. 1. The first option is to include the
linearized TPoC constraint on the right side of Fig. 1. In fact, since it is
based on a heavy linearization of the PoC formula, it is unlikely to converge
when a good reference solution is not provided, and it can only be used to
refine the solution obtained with the first convex run. An alternative
approach, reported in the left side of Fig. 1, consists of adapting the limit
PoC values $\bar{P}_{C,s}$ and running a new SCP with the new optimized
values.

Solve first SCP run with first guess for
$\displaystyle\overline{P}_{C,s}$Select constraintAdapt limits with Problem
41StopUpdate dynamics with Eq. 17SCP run - Major iterationsSCP
runConvergencewith Eq. 57Convergencewith Eq. 58Linearized SMDLinearized
TPoCMinor iterationsCheck TPoCcomplianceLimits refinementSolve SOCP Eq. 56with
TPoC constraintSolve SOCP Eq. 56with SMD constraintUpdate dynamics with Eq.
17Convergencewith Eq. 57 Figure 1: High-level flow of the algorithm.

To tune the limit of each conjunction to optimal values, a non-linear program
(NLP) is solved. In Fig. 2, the purpose of the NLP is shown. The PoC limits
imposed in the first SCP yield SMD limits that identify the dashed ellipses.
Given these two limits, the convex solver may find a solution like the blue
dots: the first conjunction dominates the scenario because $\bm{r}_{1}^{0}$ is
on the border of the ellipse, while $\bm{r}_{2}^{0}$ is far from it. After the
limits adaptation, more importance is given to the first conjunction. The
limit PoC of the first conjunction is increased so that its SMD limit
decreases, and the solution to the new SCP might be found closer to the
ballistic miss distance. The minimization of the miss distance of the second
conjunction is of no concern because, in the first run, the optimized point
was found far from the limit. So, the SMD threshold of the second is allowed
to grow (or the PoC threshold to shrink) until it reaches a value for which
the optimized point is on the limit.

$\displaystyle\Delta\bm{r}_{1}^{0}$$\displaystyle\Delta\bm{r}_{1}^{ball}$$\displaystyle\Delta\bm{r}_{1}$$\displaystyle\Delta\bm{r}_{2}^{ball}$$\displaystyle\Delta\bm{r}_{2}$$\displaystyle\Delta\bm{r}_{2}^{0}$$\displaystyle{\color[rgb]{0.55,0,0.07}\definecolor[named]{pgfstrokecolor}{rgb}{0.55,0,0.07}\bm{r}_{1}}$$\displaystyle{\color[rgb]{0.55,0,0.07}\definecolor[named]{pgfstrokecolor}{rgb}{0.55,0,0.07}\bm{r}}{\color[rgb]{0.55,0,0.07}{}_{2}}$
Figure 2: Scheme of the adaptation of the PoC limits for two conjunctions.

The optimization variables of the NLP are collected in two vectors of
parameters, $\vec{\alpha}$ and $\vec{\rho}\in\mathbb{R}^{n_{conj}}$.
$\vec{\alpha}$ is such that for the $s$-th conjunction
$\bar{P}_{C,s}=\bar{P}_{TC}\cdot 10^{-\alpha_{s}}$. We choose to introduce
this parameter as an exponent to favor the algorithm’s convergence since the
limit PoC values can span a wide range of orders of magnitude; the minus sign
of the exponent is useful to keep the optimization variables always positive.
Our working hypothesis affirms that getting closer to the ballistic relative
position would result in a lower $\Delta v$. So, the optimized relative
position is bound to lie in the same direction as the point resulted from the
previous SCP w.r.t. the ballistic relative position.

$\Delta\hat{\bm{r}}_{s}=\frac{\Delta\vec{r}^{0}_{s}-\Delta\bm{r}_{s}^{ball}}{||\Delta\vec{r}^{0}_{s}-\Delta\bm{r}_{s}^{ball}||},$
(39)

where $\Delta\bm{r}_{s}^{0}$ and $\Delta\bm{r}_{s}^{ball}\in\mathbb{R}^{3}$
are the relative position of the previous SCP run and of the ballistic
trajectory w.r.t. the secondary $s$. In this way, we only need to optimize the
magnitude of the difference between the optimized relative position and the
ballistic one ($\rho_{s}$) for each conjunction $s$

$\displaystyle\rho_{s}=||\Delta\bm{r}_{s}-\Delta\bm{r}_{s}^{ball}||\quad$
$\displaystyle s\in\\{1,...\hskip 2.0pt,n_{conj}\\}.$ (40)

$\vec{\rho}$ collects the $\rho_{s}$ for every conjunction. Now we can write
the NLP that is used to optimize the limit probabilities. In Problem 41, the
variables with superscript $0$ indicate the output of the previous SCP

$\displaystyle\min_{\bm{\rho},\bm{\alpha}}\quad$
$\displaystyle\bm{\rho}\cdot\bm{\beta}\quad$ (41a) s.t.
$\displaystyle\Delta\bm{r}_{s}=\Delta\bm{r}_{s}^{ball}+\rho_{s}\Delta\hat{\bm{r}}_{s}\quad$
$\displaystyle s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (41b)
$\displaystyle\bar{P}_{C,s}=\bar{P}_{TC}\cdot 10^{-\alpha_{s}}\quad$
$\displaystyle s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (41c) $\displaystyle
P_{C,s}(\Delta\bm{r}_{s},\bm{P}_{s},HBR_{s})\leq\bar{P}_{C,s}\quad$
$\displaystyle s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (41d) $\displaystyle
1-\bar{P}_{TC}-\prod_{s=1}^{n_{conj}}(1-\bar{P}_{C,s})=0,\quad$ (41e)
$\displaystyle\beta_{s}=\begin{cases}1&\mbox{if
}|P^{0}_{C,s}-\bar{P}^{0}_{C,s}|/\bar{P}^{0}_{C,s}\leq\varepsilon\\\
0&\mbox{otherwise}\end{cases}\quad$ $\displaystyle s\in\\{1,...\hskip
2.0pt,n_{conj}\\},$ (41f) $\displaystyle 0\leq\rho_{s}\leq\rho_{s}^{0}$
$\displaystyle s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (41g) $\displaystyle
0\leq\alpha_{s}\leq 10\quad$ $\displaystyle s\in\\{1,...\hskip
2.0pt,n_{conj}\\}.$ (41h)

The objective function in Eq. 41a is the norm of the scaled distance of the
optimized point from the ballistic one. With Eq. 41b, the optimized point is
bound to lie in the same direction as the point resulted from the previous SCP
w.r.t. the ballistic relative position. Eq. 41c defines the PoC limits as
functions of the optimization variable $\bm{\alpha}$. Eq. 41d sets the
collision avoidance constraint for each secondary. Eq. 41e imposes that the
combination of the PoC limits must give the TPoC limit. In Eq. 41f, the vector
$\bm{\beta}\in\mathbb{R}^{n_{conj}}$ is used to ensure that the conjunctions
that affect the objective function are the ones for which the PoC contributes
in a non-negligible manner to the TPoC; $\varepsilon\in\mathbb{R}$ is a
tolerance value, e.g., $10^{-2}$. The last two constraints, Eq. 41g and Eq.
41h, are bound constraints on the optimization variables.

As a solution of Problem 41, one obtains the weights in $\bm{\alpha}$ that are
used to find the limit PoC for each conjunction. Using this objective function
requires the underlying hypothesis that a solution to the convex problem that
is closer to the original ballistic trajectory would require a lower $\Delta
v$, which might not always be the case. Lastly, it is important to feed the
NLP a suitable first guess in terms of $\bm{\alpha}_{0}$ and $\vec{\rho}_{0}$:
the value of $\vec{\rho}_{0}$ is taken from the result of the previous SCP,
and $\bm{\alpha}_{0}$ is derived from the values of the PoC associated with
the single conjunctions

$\displaystyle\alpha_{0,s}=-\log\left(\frac{P_{C,s}^{0}}{\bar{P}_{TC}}\right)\quad$
$\displaystyle s\in\\{1,...\hskip 2.0pt,n_{conj}\\}.$ (42)

###  3.3 Conjunction representation with GMM

The GMM representation of the uncertainty can be used in the framework of the
multiple encounters with just a few adjustments. It is easy to visualize every
mixand as a separate secondary, thus treating the conjunctions separately,
like in the scenario where multiple encounters happen with different objects.
As shown in Fig. 3, depending on the direction of the split and on the
direction of relative motion between the two bodies, the TCAs of the different
mixands can be very close, if not identical. This is equivalent to the primary
being required to avoid the collision with more than two very close objects at
times that are close to one another. The issue arises of the previously
established time grid (on which the dynamics have been discretized) not being
able to catch the exact TCAs. This calls for refining the grid near the
nominal TCA. Fig. 4 depicts the idea of the grid refinement for a case with
three mixands: the original equally spaced time grid is represented in black.
While for the first conjunction, in this particular case, the three mixands
all arrive at the close approach at the same time, for the second conjunction,
they have had the time to drift from each other; thus, the close encounter
with the primary happens at the nominal time for the central mixand (dark
green line) and at slightly different times for the other two (light green
lines). A grid thickening algorithm is used to be able to catch the offset
TCAs, yielding new nodes like the ones represented in red.

$t_{CA,12}$$t_{CA,22}$$t_{CA,32}$$c=1$$c=2$$c=3$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}v_{12}}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}v}{\color[rgb]{0.82,0.01,0.11}{}_{22}}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}v}{\color[rgb]{0.82,0.01,0.11}{}_{32}}$
((a))

$t_{CA,12}\approx t_{CA,22}\approx
t_{CA,32}$$c=1$$c=2$$c=3$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}v_{12}}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}v}{\color[rgb]{0.82,0.01,0.11}{}_{22}}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}v}{\color[rgb]{0.82,0.01,0.11}{}_{32}}$
((b))

Figure 3: Different geometric configurations of the conjunctions depending on
the direction of relative motion.

First of all, the real TCA must be computed. One could do this by propagating
the trajectory from node $t_{CA,s}$ previous to the TCA one and the one
subsequent to it and finding the time for which the minimum separation is
reached. This method, though, is highly inefficient. We leverage DA to obtain
the same result. For a conjunction $s\in\\{1,...\hskip 2.0pt,n_{conj}\\}$, let
us take the nominal states of the two spacecraft at TCA and expand the
secondary state

$\bm{x}_{s,CA}^{j}=\tilde{\bm{x}}_{s,CA}^{j}+\delta\bm{x}_{s,CA}^{j},$ (43)

We propagate the states of the two satellites by a DA variable $\delta t_{s}$,
obtaining, with a Taylor polynomial notation

$\displaystyle\bm{x}_{s,CA}^{j*}=\mathcal{T}^{1}_{\bm{x}_{s,CA}^{j}}(\delta
t_{s},\delta\bm{x}_{s,CA}^{j}),\quad$
$\displaystyle\bm{x}_{p,CA}^{j*}=\mathcal{T}^{1}_{\bm{x}_{p,CA}^{j*}}(\delta
t_{s}),$ (44)

where the stars indicate the variables at the real TCA. The relative state is
a new Taylor polynomial dependent on $\bm{x}_{s,CA}^{j}$,
$\tilde{\bm{x}}_{p,CA}^{j}$ and $\delta t_{s}$:
$\Delta\bm{x}_{s,CA}^{j*}=\bm{x}_{p,CA}^{j*}-\bm{x}_{s,CA}^{j*}$. It can be
split into the position and velocity parts to write the equation of the
minimum distance. As mentioned in Section 3.2, the relative position is
minimum when the relative velocity is orthogonal to it

$\Delta\bm{r}_{CA}^{*}\cdot\Delta\bm{v}_{CA}^{*}=\mathcal{T}^{1}_{\Delta\bm{r}_{CA}^{*}\cdot\Delta\bm{v}_{CA}^{*}}(\delta
t_{s},\delta\bm{x}_{s,CA}^{j})=0.$ (45)

This parametric implicit equation can be solved for $\delta t_{s}$ using a
polynomial partial inversion technique, as the one used in [46] and [49]. The
equation for $\delta t_{s}$ reads

$\delta t_{s}=\mathcal{T}^{1}_{\delta
t_{s}}(\Delta\bm{r}_{CA}^{*}\cdot\Delta\bm{v}_{CA}^{*},\delta\bm{x}_{s,CA}^{j}).$
(46)

The exact TCA of the single mixand, then, is

$t_{CA,s}^{*}=t_{CA,s}+\delta t_{s}.$ (47)

$t_{0}$
---
$t_{CA,1}$
---
$t_{CA,2}$
---
Figure 4: Grid adjustment to account for different TCAs of the mixands.

According to the new TCAs, nodes are added to the time grid, and the states of
the primary and the relevant mixands are propagated to such times.

If a GMM model is used to propagate the uncertainty, we can treat each mixand
as an independent secondary object. This way, we can navigate back to the
algorithm described in the previous section. If we want to adapt the limits of
the mixands to refine the solution, we can use Problem 41, with two slight
differences. Firstly, in the first guess, we use the approximated version of
Eq. 14a and the limits are set proportional to the weights

$\displaystyle\bar{P}_{C,cs}=\frac{\gamma_{c}}{n_{conj}}\bar{P}_{TC}\quad$
$\displaystyle c\in\\{1,...\hskip 2.0pt,n_{mix}\\}$ $\displaystyle
s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (48)

where $n_{mix}$ is the number of GMM mixands, $s$ indicates the conjunction
and $c$ indicates the mixand; this is not a necessary condition since any
combination of thresholds is valid as long as it respects the total
probability threshold, but it is found to favor convergence in practice.666For
example, the equation $\bar{P}_{C,cs}=\bar{P}_{TC}/(n_{mix}\cdot n_{conj})$
would be equally valid. Secondly, in Eq. 41d, the PoC of the conjunction
obtained with Chan’s formulas must be multiplied by its weight

$\displaystyle\gamma_{c}P_{C,cs}(\bm{r}_{cs},\bm{C}_{cs},HBR_{s})\leq\bar{P}_{C,cs}\quad$
$\displaystyle c\in\\{1,...\hskip 2.0pt,n_{mix}\\}$ $\displaystyle
s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (49)

and in Eq. 41e, the nested product must be included to account for the GMM
components in each conjunction

$1-\bar{P}_{TC}-\prod_{s=1}^{n_{conj}}\prod_{c=1}^{n_{mix}}(1-\bar{P}_{C,cs})=0.$
(50)

If Eq. 48 is substituted into Eq. 50 we see that the total limit is respected
by the initial arbitrary split.

One last remark about the optimization with the GMM method regards the
computation of the SMD limit. To obtain a realistic SMD limit from the PoC
limit with the inversion of Chan’s formula, the mixand weight must be taken
into account

$\displaystyle\bar{d}_{m,cs}^{2}=f\left(\frac{\bar{P}_{C,cs}}{\gamma_{c}}\right)\quad$
$\displaystyle c\in\\{1,...\hskip 2.0pt,n_{mix}\\}$ $\displaystyle
s\in\\{1,...\hskip 2.0pt,n_{conj}\\}.$ (51)

###  3.4 Convexifcation of the PoC Constraint for Long-Term Encounters

When dealing with the long-term encounter problem, the adopted approach is
derived from the one described in the previous Section 3.2, with some
substantial differences presented in this section. A detailed description of
the projection and linearization algorithm for the single conjunction in long-
term encounters is given in reference [19]. Here, it suffices to say that it
is analogous to the one presented in Section 3.2.1, but instead of working in
the B-plane, we set the problem in the 3D ECI reference frame. There is no
particular need to use an ECI frame, but it is convenient because it does not
require transformations with respect to the frame used for the propagation of
the dynamics.

####  3.4.1 Linearized TIPoC Constraint

A similar approach to Section 3.2.2 is used for the TIPoC constraint in the
case of long-term encounters. In this case, the constraint has to be applied
indistinctly at every node of the optimization, so the stack of positions
changes node-wise $\mathbb{r}\in\mathbb{R}^{3n_{conj}\times N}$ and so does
the Jacobian $\bm{J}\in\mathbb{R}^{3n_{conj}\times N}$. The node-wise
linearized TIPoC constraint reads

$\displaystyle
P_{TIC,i}=\bm{\nabla}_{i}\mathbb{r}_{i}+d_{i}\leq\bar{P}_{TIC},\quad$
$\displaystyle i=\\{1,...\hskip 2.0pt,N\\}$ (52)

where $\bm{\nabla}_{i}\in\mathbb{R}^{3n_{conj}}$ is the gradient of the
expression that gives the TIPoC as a function of the relative position w.r.t.
each secondary and $d=\tilde{P}_{TIC,i}-\bm{\nabla}_{i}\tilde{\mathbb{r}}_{i}$
is the residual of the linearization. The same trust region method from
Section 3.2.2 is employed, bearing in mind that, in this case, the trust
region must be applied to the relative position variables at each node of the
optimization.

####  3.4.2 Refinement of the solution

When dealing with multiple secondaries or propagating GMM components, also in
long-term encounters, it is important to refine the solution with one of the
two methods introduced in Section 3.2.3. The first one, based on the
linearization of the TIPoC, was shown in Section 3.4.1. The second one
consists of adapting the IPoC threshold of each component to obtain the
minimum $\Delta v$ maneuver. An analogous NLP to the one presented in Problem
41 is used to tune the weights of the limits after the problem is solved with
a first guess for the limits. The first guess, as in the short-term scenario,
is for the IPoC limits to be proportional to the mixand weights (using the
approximated version of Eq. 14b)

$\displaystyle\bar{P}_{IC,cs}=\frac{\gamma_{c}}{n_{conj}}\bar{P}_{TIC}\quad$
$\displaystyle c\in\\{1,...\hskip 2.0pt,n_{mix}\\}$ $\displaystyle
s\in\\{1,...\hskip 2.0pt,n_{conj}\\}.$ (53)

Only the optimization node where the TIPoC is maximum is used in the NLP
because the maximum point in the profile is expected to stay the same or be
shifted by a short time. This single value of TIPoC is indicated as
$\hat{P}_{TIC}=\max_{i}(P_{TIC_{i}})$, where $i\in[0,N]$ and the correspondent
values of IPoC for the single secondaries are $\hat{P}_{IC_{s}}$

$\displaystyle\min_{\rho,\bm{\alpha}}\quad$ $\displaystyle|\rho|\quad$ (54a)
s.t. $\displaystyle\bm{r}=\bm{r}^{ball}+\rho\hat{\bm{r}},\quad$ (54b)
$\displaystyle\bar{P}_{IC,cs}=\bar{P}_{TIC}\cdot 10^{-\alpha_{cs}}\quad$
$\displaystyle c\in\\{1,...\hskip 2.0pt,n_{mix}\\}\quad$ $\displaystyle
s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (54c)
$\displaystyle\gamma_{c}\hat{P}_{IC,cs}(\Delta\bm{r}_{cs},\bm{P}_{cs},HBR_{s})\leq\bar{P}_{IC,cs}\quad$
$\displaystyle c\in\\{1,...\hskip 2.0pt,n_{mix}\\}\quad$ $\displaystyle
s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$ (54d) $\displaystyle
1-\bar{P}_{TIC}-\prod_{s=1}^{n_{conj}}\prod_{c=1}^{n_{mix}}(1-\bar{P}_{IC,cs})=0,\quad$
(54e) $\displaystyle 0\leq\rho\leq\rho^{0},$ (54f) $\displaystyle
0\leq\alpha_{cs}\leq 10\quad$ $\displaystyle c\in\\{1,...\hskip
2.0pt,n_{mix}\\}\quad$ $\displaystyle s\in\\{1,...\hskip 2.0pt,n_{conj}\\},$
(54g)

where the position variables
$\bm{r},\vec{r}^{ball},\Delta\bm{r}_{cs}\in\mathbb{R}^{3}$ are all computed at
the time step of maximum TIPoC. In this case, the objective function in Eq.
54a is the minimization of the scalar variable $\rho$. The remaining equations
from Eq. 54b to Eq. 54g can be interpreted in the same way as those in Problem
41. The same considerations from Section 3.3 are also valid for the long-term
encounters when propagating with a GMM.

###  3.5 Final Form of the Problem

This section summarizes the problem we face and provides some insights into
the optimization process. First, we finalize the low-level optimization
problem, the SOCP. The complete objective function is given by the sum of Eq.
23 and Eq. 27b:

$J=\kappa_{vc}\sum_{i=0}^{N}\upsilon_{i}+\sum_{i=0}^{N-1}u_{i}\cdot\Delta
t_{i},$ (55)

where the node-wise control action is weighted to account for variable time
steps. The final SOCP is reported in Problem (56).

$\displaystyle\min_{\mathbb{x},\mathbb{u}}$ Eq. 55 (56) s.t. Eq. 17, Eq. 29 or
Eq. 32 with Eq. 36, Eqs. 21, 22, 24, 25 and 27a

The major iterations’ convergence criterion is computed as the difference
between the maximum control acceleration of two consecutive iterations

$e_{M}=\max_{i}(||\bm{u}_{i}^{j}-\bm{u}_{i}^{j-1}||).$ (57)

Likewise, the minor iterations’ criterion (for the linearized SMD method) is
the difference between the relative position vectors of two consecutive minor
iterations

$e_{m}=\max_{i}(||\bm{r}_{i}^{j}-\bm{r}_{i}^{j-1}||).$ (58)

It is relevant to remember that the successive linearization process is
performed on the output of the previous optimization problem, so there is no
online check on the accuracy of the forward-propagated solution. Such a check
is performed after the optimization has reached convergence in all its levels.
A validation error is computed as the maximum difference between the states at
the last output of the optimization and the states obtained by forward-
propagating the trajectory using the converged control profile

$e_{validation}=\max_{i}(||\bm{x}_{i}^{j=end}-\bm{x}_{i}^{val}||).$ (59)

To conclude the discussion, the high-level work-flow of the algorithm, as it
has been presented in the previous sections, is shown in Fig. 1.

##  4 Results

Realistic test cases for multiple encounters in LEO are analyzed. Three
scenarios are considered. In the first one, the primary must deal with
multiple short-term conjunctions with different secondaries; in the second
one, the primary encounters the secondary more than once, and the
uncertainties are propagated using GMM; in the last one, a long term encounter
with a single secondary is studied using GMMs. Results from both the
refinement methods from Section 3.2 and Section 3.4 are compared. In all the
scenarios, the discretization time step is such that the orbit of the primary
is divided into $60$ nodes. The threshold TPoC and TIPoC are both set to
$10^{-6}$. The tolerance on the major iterations is set to $10^{-3}$, and the
one on the minor iterations for the linearized SMD method is set to $10^{-6}$.
The simulations are run with MATLAB r2022b on AMD Ryzen 9 6900HS @ 3.3GHz. The
optimization is performed using MOSEK 10.0.24, which implements a state-of-
the-art primal-dual interior-point solver. Scaled variables are used both
during the integration of the equations of motion and in the optimization to
favor numerical stability and convergence: the scaling constants are
$L_{sc}=a_{p}$, $V_{sc}=\sqrt{\mu_{E}/a_{p}}$,
$T_{sc}=\sqrt{a_{P}^{3}/\mu_{E}}$, and $A_{sc}=\mu_{E}/a_{P}^{2}$, where
$a_{P}$ and $\mu_{E}$ are the semi-major axis of the primary’s orbit and
Earth’s gravitational constant. In all the considered scenarios, the
trajectories are discretized using a constant time step into 60 nodes per
orbit.

Table 1: Case 1: conjunctions’ parameters. | $t_{CA}$ |  $\Delta\bm{r}$ [$\mathrm{m}$] |  $\Delta\bm{v}$ [$\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$] |  $d_{miss}$ [$\mathrm{m}$] |  $P_{C}$ [-] |  $\gamma_{v}$ [deg]
---|---|---|---|---|---|---
CDM 1 | 1 $\mathrm{h}$ 34 $\mathrm{m}\mathrm{i}\mathrm{n}$ 3 $\mathrm{s}$ | [-11.45; -19.44 7.44] | [ 0.37; -2.85; 4.32] | 23.8 | 0.0017 | 40
CDM 2 | 2 $\mathrm{h}$ 4 $\mathrm{m}\mathrm{i}\mathrm{n}$ 20 $\mathrm{s}$ | [ 99.90; 14.60 2.78] | [-3.94; -5.85; 3.36] | 101.0 | 0.0018 | 62
CDM 3 | 3 $\mathrm{h}$ 43 $\mathrm{m}\mathrm{i}\mathrm{n}$ 11 $\mathrm{s}$ | [ 45.70; -3.73 28.29] | [-0.15; -1.30; 0.90] | 53.9 | 0.0017 | 12
CDM 4 | 4 $\mathrm{h}$ 46 $\mathrm{m}\mathrm{i}\mathrm{n}$ 55 $\mathrm{s}$ | [ 8.90 49.27 -50.88] | [ 0.93; -0.63; 11.56] | 71.4 | 0.0025 | 100
CDM 5 | 5 $\mathrm{h}$ 53 $\mathrm{m}\mathrm{i}\mathrm{n}$ 54 $\mathrm{s}$ | [-8.10 -3.02 -20.62] | [ 13.86; -3.71; -4.92] | 22.4 | 0.0138 | 180
CDM 6 | 9 $\mathrm{h}$ 22 $\mathrm{m}\mathrm{i}\mathrm{n}$ 43 $\mathrm{s}$ | [-12.64 -14.81 13.76] | [ 5.64; -3.00 8.62] | 23.8 | 0.0023 | 90
CDM 7 | 11 $\mathrm{h}$ 47 $\mathrm{m}\mathrm{i}\mathrm{n}$ 47 $\mathrm{s}$ | [ 2.09 13.74 1.89] | [ 7.93; -7.06; 3.21] | 14.0 | 0.0022 | 94
CDM 8 | 13 $\mathrm{h}$ 17 $\mathrm{m}\mathrm{i}\mathrm{n}$ 3 $\mathrm{s}$ | [-4.18 -14.19 14.49] | [ 0.37; -2.23; 1.35] | 20.7 | 0.0020 | 20
CDM 9 | 14 $\mathrm{h}$ 46 $\mathrm{m}\mathrm{i}\mathrm{n}$ 20 $\mathrm{s}$ | [ 8.92 -1.02 -12.78] | [ 1.43; 4.52; 13.58] | 15.6 | 0.0050 | 143
CDM 10 | 15 $\mathrm{h}$ 56 $\mathrm{m}\mathrm{i}\mathrm{n}$ 28 $\mathrm{s}$ | [-19.88 9.10 -4.18] | [ 0.63; -9.52; -0.31] | 22.3 | 0.0012 | 78
|  |  |  |  |  |
Table 2: Case 1: elements of the relative position covariance matrix of each conjunction. |  $P_{rr}$ [$\mathrm{k}\mathrm{m}^{2}$] |  $P_{tt}$ [$\mathrm{k}\mathrm{m}^{2}$] |  $P_{nn}$ [$\mathrm{k}\mathrm{m}^{2}$] |  $P_{rt}$ [$\mathrm{k}\mathrm{m}^{2}$] |  $P_{tn}$ [$\mathrm{k}\mathrm{m}^{2}$] |  $P_{nr}$ [$\mathrm{k}\mathrm{m}^{2}$]
---|---|---|---|---|---|---
CDM 1 | 0.0078 | 0.2887 | 0.0145 | 0.0466 | 0.0103 | 0.0644
CDM 2 | 0.0662 | 0.1279 | 0.1170 | -0.0918 | 0.0879 | -0.1220
CDM 3 | 0.2473 | 0.0022 | 0.0616 | 0.0215 | 0.1232 | 0.0109
CDM 4 | 0.0004 | 0.1457 | 0.1649 | -0.0055 | 0.0057 | -0.1548
CDM 5 | 0.2579 | 0.0199 | 0.0333 | -0.0710 | -0.0924 | 0.0256
CDM 6 | 0.0001 | 0.2000 | 0.1110 | -0.0007 | 0.0005 | -0.1486
CDM 7 | 0.0015 | 0.2027 | 0.1069 | -0.0162 | 0.0120 | -0.1469
CDM 8 | 0.1814 | 0.0012 | 0.1284 | -0.0142 | -0.1523 | 0.0119
CDM 9 | 0.0028 | 0.0002 | 0.3081 | -0.0003 | 0.0278 | -0.0017
CDM 10 | 0.0004 | 0.1327 | 0.1779 | -0.0036 | 0.0040 | -0.1535
|  |  |  |  |  |

###  4.1 Case 1: Multiple Short-Term Encounters in LEO

Most impacts in LEO are head-to-head [1], so they involve very high relative
velocities. The test case presented simulates a series of endogenous
conjunctions in the Starlink constellation under $J2$ dynamics. A total of ten
conjunctions happening in the course of ten orbits of the primary is
synthetically generated. The primary is on a circular orbit with an
inclination of $53$ $\mathrm{d}\mathrm{e}\mathrm{g}$, a radius of $6928$
$\mathrm{k}\mathrm{m}$, a RAAN, and an argument of periapsis equal to $0$
$\mathrm{d}\mathrm{e}\mathrm{g}$; at the starting time, $t_{0}$, its true
anomaly is $0$ $\mathrm{d}\mathrm{e}\mathrm{g}$. The parameters that define
each conjunction are reported in Table 1 and Table 2: the relative states are
expressed in the ECI reference frame, whereas the covariances are in the
radial, along-track, cross-track (RTN) of the secondary. The covariance of the
primary is assumed to be already absorbed by that of the secondary, so the
number of parameters needed to set the scenario is reduced.
$d_{miss}=||\Delta\bm{r}_{s,CA}||$ is the miss distance of the un-maneuvered
case, and
$\gamma_{v}=\mathrm{acos}\left(\vec{v}_{p,CA}\cdot\bm{v}_{s,CA}\right)$ is the
angle between the velocity vectors of the two spacecraft; since the encounters
happen between satellites on circular orbits at the same altitude, the
magnitude of the velocity of the two spacecraft is equal. The HBR is the same
for each Starlink satellite and is conservatively assumed to be $3$
$\mathrm{m}$. The propagation is performed using the harmonics up to the
second degree of the gravitational field. The spacecraft is supposed to mount
a low-thruster with a maximum thrust of $5.2$ $\mathrm{m}\mathrm{N}$, which,
considering the Starlink’s mass of 230 kg, corresponds to a maximum
acceleration of $0.02$ $\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}^{2}$.

((a))

((b))

((c))

((d))

Figure 5: Case 1: equivalent B-plane representation of the conjunctions.

The algorithm is run with both the refinements methods introduced in Section
3.2, and the results are compared.

In Fig. 5, the equivalent B-planes of the conjunctions for each case
considered are shown. To represent all the conjunctions in the same plot, a
coordinate transformation is performed on the KOZ ellipse and on the relative
position points. The transformation comprises a rotation to diagonalize the
covariance matrix and a stretch to circularize it. Considering a generic
vector $\bm{w}=[\xi\hskip 4.0pt\zeta]^{\mathrm{T}}$ expressed in the B-plane
reference frame of a conjunction $s$, the transformed vector
$\tilde{\bm{w}}=[\tilde{\xi}\hskip 4.0pt\tilde{\zeta}]^{\mathrm{T}}$ is

$\tilde{\bm{w}}=\mathrm{diag}\left([1/a;\hskip 4.0pt1/b]\right)\vec{V}\bm{w},$
(60)

where $a,b\in\mathbb{R}_{\mathrm{+}}$ are the semiaxes of the ellipse
associated with the optimized probability threshold, and
$\vec{V}\in\mathbb{R}^{2\times 2}$ is the rotation matrix that diagonalizes
the covariance. The ballistic trajectory points are shown as diamonds, and the
optimized points obtained using the linearized SMD constraint with limit
adaptation are shown as circles. The conjunction number is indicated for each
marker, and the color of the marker is associated with the value of the limit
PoC of the corresponding conjunction after the limit adaptation. The result of
the same limit adaptation is reported in Fig. 6: after the first SCP is
solved, in all the considered cases, a single SCP is needed to find a solution
that grants a TPoC value sufficiently close to the threshold. We notice from
Fig. 6 that depending on the conjunctions configuration and their number, the
algorithm automatically finds the most important ones and assigns them a
higher weight. For example, the second conjunction is equally important to the
first one in the first scenario, but it becomes irrelevant in the next ones.
As a second example, the last three conjunctions that differentiate the last
scenario from the third have no effect on the final result since their
$\bar{P}_{C}$ value drops significantly after the adaptation. This is why the
computed maneuvers for the last two scenarios are almost identical (see Fig.
7(c) and Fig. 7(d)).

Through an analysis of this type, it is also possible to identify clusters of
conjunctions that cause a maneuver. Indeed, the time, direction, and magnitude
of the first two burns remain almost unchanged across the four cases. In
general, the maneuvers are all mainly in the tangential direction. This
confirms the idea proposed by other researchers [2, 4] that raising or
lowering the orbit is the most effective way of performing the CAM if enough
warning time is given.

((a))

((b))

((c))

((d))

Figure 6: Case 1: evolution of the $\bar{P}_{C,s}$ of the single conjunctions.

In Table 3, the results in terms of final TPoC and $\Delta v$ are shown for
the four combinations of conjunctions considered. Convergence results are also
compared in the number of iterations, convergence time, and maximum positional
validation error. The two refinement methods do not differ greatly in the
resulting $\Delta v$ profile or in the final equivalent B-plane configuration,
as is shown in Table 3. Nonetheless, the method that relies on the linearized
SMD constraint with limit adaptation yields more accurate results and can
target the desired TPoC threshold with higher precision. The PoC value at TCA
for the first method is always kept below the threshold value with a margin of
$.002\%$, while the linearized TPoC method reaches TPoC values up to $42\%$
below the threshold. This means that the CAM optimization does not find the
optimal result with the second method, leading to higher $\Delta v$ impulses
than the first method. The run times of the two methods are similar when a low
number of conjunctions is considered, but the linearized SMD method is
definitely faster when we start to consider a high $n_{conj}$. In general, we
can infer that the linearized SMD constraint is more appropriate when dealing
with more than two conjunctions because each PoC constraint is linearized
separately, so the coarseness of approximation of the nonlinearity does not
scale with the complexity of the problem.

Table 3: Case 1: compared results for the SMD and TPoC linearized constraints. Constraint | Linearized SMD | Linearized TPoC
---|---|---
$n_{conj}$ | 2 | 5 | 7 | 10 | 2 | 5 | 7 | 10
$\Delta v$ [$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$] | $21.16$ | $60.05$ | $102.04$ | $102.14$ | $21.16$ | $62.32$ | $103.18$ | $105.58$
$P_{TC}\times 10^{6}$ [-] | $1.001$ | $1.001$ | $1.001$ | $0.998$ | $1.002$ | $0.801$ | $0.824$ | $0.583$
$n_{major}$ [-] | 4 | 6 | 6 | 8 | 5 | 7 | 10 | 20
$n_{minor}$ [-] | 5 | 9 | 10 | 13 | - | - | - | -
$e_{validation}$ [$\mathrm{m}\mathrm{m}$] | $0.92$ | $1.88$ | $0.72$ | $5.12$ | $37.61$ | $5.48$ | $43.40$ | $1940.80$
$T_{sim}$ [$\mathrm{s}$] | $4.53$ | $8.34$ | $14.52$ | $28.03$ | $4.47$ | $8.82$ | $20.66$ | $62.15$
|  |  |  |  |  |  |  |

((a))

((b))

((c))

((d))

Figure 7: Case 1: $\Delta v$ to perform the CAM.

###  4.2 Case 2: Repeating Short-Term Encounter in LEO with GMM

Let us consider a scenario in which a secondary spacecraft’s orbit intersects
the primary’s twice, thus generating a repeating encounter. The first
encounter is described by a CDM, whereas the second one is detected by
propagating the dynamics. The covariance of the secondary at the first TCA is
split into a number of GMM mixands, which are propagated separately until the
second TCA. Here, we assume that the covariance of the primary is negligible
w.r.t. to the secondary’s. The state of the primary is described at TCA by
$a=6800$ $\mathrm{k}\mathrm{m}$, $e=0$, $\omega=\Omega=\theta=0$
$\mathrm{d}\mathrm{e}\mathrm{g}$, and $i=91.67$
$\mathrm{d}\mathrm{e}\mathrm{g}$; the orbit of the secondary is elliptical,
has opposite inclination w.r.t. the primary and its period is a multiple of
the period of the primary. The full secondary’s covariance matrix at the first
TCA is defined in RTN coordinates in Table 4. The primary covariance matrix is
assumed to be null.

Table 4: Case 2: elements of the relative position covariance matrix at the first TCA. Spacecraft |  $P_{rr}$ [$\mathrm{m}^{2}$] |  $P_{tt}$ [$\mathrm{m}^{2}$] |  $P_{nn}$ [$\mathrm{m}^{2}$] |  $P_{\dot{r}\dot{r}}$ [$\mathrm{m}\mathrm{m}^{2}\mathrm{/}\mathrm{s}^{2}$] |  $P_{\dot{t}\dot{t}}$ [$\mathrm{m}\mathrm{m}^{2}\mathrm{/}\mathrm{s}^{2}$] |  $P_{\dot{n}\dot{n}}$ [$\mathrm{m}\mathrm{m}^{2}\mathrm{/}\mathrm{s}^{2}$]
---|---|---|---|---|---|---
Primary | 0 | 0 | 0 | 0 | 0 | 0
Secondary | 2.025 | 10.000 | 2.500 | 0.225 | 5.625 | 0.625
|  |  |  |  |  |

According to reference [21], when propagating Cartesian coordinates in LEO,
for typical initial uncertainty, after five orbits ca., the normality of the
distribution is lost. So, we construct the test such that after six orbits of
the primary, the head-to-head encounter happens again. The eccentricity of the
secondary orbit is $0.6971$; its semi-major axis is $22453.1$
$\mathrm{k}\mathrm{m}$, the inclination is $181.67$
$\mathrm{d}\mathrm{e}\mathrm{g}$, the argument of perigee is $0$
$\mathrm{d}\mathrm{e}\mathrm{g}$, the RAAN is $1.72\cdot 10^{-3}$
$\mathrm{d}\mathrm{e}\mathrm{g}$ and the true anomaly at the first conjunction
is $0.11$ $\mathrm{d}\mathrm{e}\mathrm{g}$. The orbits are propagated backward
for one orbital period of the primary and forward up to the second
conjunction. In this way, the first conjunction is treated as a standard one,
i.e., without splitting, and the second one is treated with the propagated GMM
components. In this case, we are considering a chemical propulsion system with
a maximum instantaneous acceleration of $1$
$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}^{2}$. Different simulations with
different $n_{mix}$ are compared to one where no split is performed, i.e., the
second conjunction is treated with a linearly propagated covariance. Moreover,
the solutions obtained with the two constraints are compared in Table 5.

((a))

((b))

((c))

((d))

Figure 8: Case 2: equivalent B-plane representation of the conjunctions.

In Table 5, we notice that the use of the GMM permits better to represent the
KOZ and reduce the required $\Delta v$ by at least $4$
$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$. Fig. 8 shows the equivalent
B-planes of the two conjunctions for different values of $n_{mix}$. The first
conjunction is always dominant: since it appears first avoiding it means that
the satellite will most likely avoid the subsequent one as well. Nonetheless,
in the other plots, it is useful to see which mixands play a role in defining
the limits: in particular, in Fig. 8(b) mixands 3 and 2, in Fig. 8(c) mixands
4 and 3 and in Fig. 8(d) mixands 5 and 4 are the ones that are not led to the
lowest threshold. This means that in the second (extended) conjunction, the
primary passes closest to the central mixand and the one right after it
compared to all the others. For all the other mixands, the threshold PoC is
brought down to the minimum value of $10^{-9}$, meaning they do not play a
role in defining the maneuver.

((a))

((b))

Figure 9: Case 2: evolution of the $\bar{P}_{C,s}$ of the single conjunctions
and mixands.

Depending on the number of mixands used, the TCAs of the single mixands are
between $1.32$ and $2.91$ $\mathrm{s}$ apart from each other, which is
accounted for in the selection of the time grid for the discretization. This
happens because the encounter is head-to-head, and the split is performed very
close to the direction of the velocity of the secondary so that the primary
runs into each mixand consecutively, as Fig. 3(a) illustrates. In Fig. 9, the
adaptation of the limits for two representative numbers of mixand is shown:
with no GMM split, the first conjunction dominates the final TPoC computation,
while in the case with $n_{mix}=7$, the second conjunction is considered
slightly more important. In Table 6, one sees that, in accordance with most
published research on single encounters [4], in all the considered cases, the
$\Delta v$ is almost purely tangential. Nonetheless, in the case where no
split is performed, the direction of thrust is opposite to that of the other
cases. In Fig. 8(a), the reader can see that the tangential maneuver moves the
relative position along the $\tilde{\xi}$ axis, so there is little to no
difference in performing a maneuver on the positive or the negative tangential
direction and the optimizer prefers one over the other. In all cases with a
split, the preferred maneuver is in the positive tangential direction. The
number of mixands considered does not affect the maneuvering time. The
maneuver is always performed half an orbit before the first TCA, and an almost
negligible second maneuver is performed half an orbit before the second TCA.

Regarding the convergence properties, the two methods behave similarly to the
case in Section 4.1, with the first method being slightly faster when more
mixands are considered. In this case, though, it is important to point out
that because of the coarseness of the linearization, the linearized TPoC
method always reaches a solution that is around $5\%$ over the TPoC threshold.
This allows it to find a lower $\Delta v$ at the cost of precision. The
validation error for both methods is always acceptable.

Table 5: Case 2: compared results for the SMD and TPoC linearized constraints. Constraint | Linearized SMD | Linearized TPoC
---|---|---
$n_{mix}$ | 1 | 3 | 5 | 7 | 1 | 3 | 5 | 7
$\Delta v$ [$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$] | $57.06$ | $52.64$ | $52.98$ | $53.37$ | $56.52$ | $52.23$ | $52.55$ | $52.15$
$P_{TC}\times 10^{6}$ [-] | $1.002$ | $0.968$ | $0.936$ | $0.910$ | $1.053$ | $1.039$ | $1.038$ | $1.048$
$n_{major}$ [-] | 5 | 8 | 8 | 8 | 7 | 7 | 8 | 9
$n_{minor}$ [-] | 12 | 16 | 16 | 19 | - | - | - | -
$e_{validation}$ [$\mathrm{m}\mathrm{m}$] | $9.92$ | $3.64$ | $3.21$ | $2.98$ | $0.46$ | $3.81$ | $4.81$ | $10.61$
$T_{sim}$ [$\mathrm{s}$] | $6.57$ | $19.76$ | $23.01$ | $27.68$ | $9.15$ | $16.98$ | $24.58$ | $30.17$
|  |  |  |  |  |  |  |
Table 6: Case 2: $\Delta v$ to perform the CAM. $n_{mix}$ |  $t_{1}$ [s] |  $\mathrm{R}_{1}$ [$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$] |  $\mathrm{T}_{1}$ [$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$] |  $\mathrm{N}_{1}$ [$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$] |  $t_{2}$ [s] |  $\mathrm{R}_{2}$ [$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$] |  $\mathrm{T}_{2}$ [$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$] |  $\mathrm{N}_{2}$ [$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$]
---|---|---|---|---|---|---|---|---
$1$ | $2790.26$ | $-0.69$ | $-56.63$ | $-0.04$ | $25112.31$ | $0$ | $-0.04$ | $0$
$3$ | $2790.26$ | $-0.68$ | $52.60$ | $0.01$ | $25112.31$ | $0$ | $0.03$ | $0$
$5$ | $2790.26$ | $-0.67$ | $53.04$ | $0.01$ | $25112.31$ |  $0$ $0.08$ |  | $0$
$7$ | $2790.26$ | $-0.01$ | $53.37$ | $0.25$ | - | - | - | -
|  |  |  |  |  |  |  |

###  4.3 Case 3: Single Long-Term Encounter in LEO with GMM

In this scenario, a single long-term conjunction in LEO is studied with GMM.
The problem becomes equivalent to having a long-term conjunction with multiple
secondaries. The scenario is taken from the LEO test case of reference [19]:
the primary spacecraft is on a circular LEO with $a=6800$
$\mathrm{k}\mathrm{m}$, and $\omega=\Omega=i=\theta=0$
$\mathrm{d}\mathrm{e}\mathrm{g}$; the relative state of the secondary as
expressed in the primary’s RTN at TCA, i.e., at the starting time is

$\bm{x}_{rel}=[-1.004814\hskip 5.0pt0\hskip 5.0pt0\hskip 5.0pt0.112790\hskip
5.0pt3.388089\hskip 5.0pt0.011228]^{\mathrm{T}},$

where the first three elements are in [$\mathrm{k}\mathrm{m}$] and the last
three in [$\mathrm{m}\mathrm{/}\mathrm{s}$]. The covariance matrices of the
two spacecraft are reported in Table 7. The combined HBR is $32$ $\mathrm{m}$,
and the maximum instantaneous acceleration achievable by the primary is 0.2
$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}^{2}$, corresponding to a maximum
thrust of $40$ $\mathrm{m}\mathrm{N}$ in a $200$ $\mathrm{k}\mathrm{g}$
satellite. The dynamics include J2 perturbation, lunisolar attraction,
atmospheric drag, and solar radiation pressure (SRP). The models used to
compute the perturbation forces are computed with Accurate Integrator for
Debris Analysis (AIDA) [50]. The physical parameters of the spacecraft are
reported in Table 8.

Table 7: Case 3: elements of the relative position covariance matrix of the two spacecraft. Spacecraft |  $C_{rr}$ [$\mathrm{m}^{2}$] |  $C_{tt}$ [$\mathrm{m}^{2}$] |  $C_{nn}$ [$\mathrm{m}^{2}$] |  $C_{\dot{r}\dot{r}}$ [$\mathrm{m}^{2}\mathrm{/}\mathrm{s}^{2}$] |  $C_{\dot{t}\dot{t}}$ [$\mathrm{m}^{2}\mathrm{/}\mathrm{s}^{2}$] |  $C_{\dot{n}\dot{n}}$ [$\mathrm{m}^{2}\mathrm{/}\mathrm{s}^{2}$]
---|---|---|---|---|---|---
Primary | $0.625$ | $10$ | $3.025$ | $0.00625$ | $0.05625$ | $0.00225$
Secondary | $5.625$ | $90$ | $27.225$ | $0.05625$ | $0.50625$ | $0.02025$
|  |  |  |  |  |
Table 8: LEO scenario: physical properties of the spacecraft. Spacecraft |  $m$ [$\mathrm{k}\mathrm{g}$] |  $A_{drag}$ [$\mathrm{m}^{2}$] |  $C_{D}$ [-] |  $A_{SRP}$ [$\mathrm{m}^{2}$] |  $C_{r}$ [-] | HBR$\mathrm{m}$
---|---|---|---|---|---|---
Primary | $200$ | $1$ | $2.2$ | $1$ | $1.31$ | $25$
Secondary | $50$ | $0.05$ | $2$ | $0.05$ | $1.31$ | $7$
|  |  |  |  |  |

The propagation spans six orbits of the primary from TCA. Different orders of
the GMM are compared in terms of $\Delta v$, TIPoC profile, and convergence
properties in Table 9 and Fig. 10. It is interesting to notice from Fig. 10(a)
that, when a low $n_{mix}$ is used, after four orbits, the profile tends to
lose its periodicity property, which is instead preserved for a high
$n_{mix}$. This suggests that the 3D KOZ starts to become too wide when no GMM
is used, while it is represented more accurately when $n_{mix}$ is high
enough. In Fig. 11, the contribution of each mixand to the ballistic TIPoC is
shown for three (Fig. 11(a)) and seven (Fig. 11(b)) mixands. It is clear that
having more mixands allows us to refine where the majority of the probability
resides. In Fig. 11(a), the total probability is mostly given by the first
mixand. Instead, in Fig. 11(b), we see that most of the mixands do not
contribute at all, being below $10^{-9}$, and the highest contributions are
given by mixands 2, 3, and 1. Moreover, using seven mixands allows us to
preserve the profile’s periodicity, which starts to fade when three mixands
are used because $P_{IC,3}$ has a non-negligible effect given by the inflation
of the corresponding covariance.

((a))

((b))

((c))

Figure 10: Case 3: TIPoC profile for different $n_{mix}$ in the long-term
encounter scenario.

In Fig. 12, the $\Delta v$ profiles for the four considered cases are
reported. The computed maneuver is very similar, and the normal component is
always predominant. From Table 9, we see that the linearized SMD approach is
much more efficient in this test case: the solution requires fewer major
iterations, and the $\Delta v$ is always lower. Since there is no guarantee
for either method to find the global optimum, we observe that the solution
found in the two cases falls into two different local minimum wells. In fact,
the TIPoC profile in the case with $n_{mix}=1$ had three equivalent points of
global maxima. The linearized SMD method’s solution for any $n_{m}ix$
preserves the second maximum (around $1.3$ orbital periods), while the
linearized TIPoC one preserves the first maximum (around $0.73$ orbital
periods). The solution found by the first method is closer to the global
optimum. Moreover, most importantly, the second method with $n_{mix}=7$ finds
a solution that is far from the local optimum achieved by $n_{mix}=3$ and
$n_{mix}=5$, so the linearization starts to be too coarse when more than five
mixands are considered in this test case.

((a))

((b))

Figure 11: Case 3: Breakdown of the contribution of each mixand to the TIPoC.
Zoom in a central orbit.

((a))

((b))

((c))

((d))

Figure 12: Case 3: $\Delta v$ to perform the CAM. Table 9: Case 3: compared results for the SMD and TIPoC linearized constraints. Constraint | Linearized SMD | Linearized TIPoC
---|---|---
$n_{mix}$ [-] | $1$ | $3$ | $5$ | $7$ | $3$ | $5$ | $7$
$\Delta v$ [$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$] | $339.42$ | $342.30$ | $344.12$ | $367.77$ | $401.76$ | $420.78$ | $469.08$
$\max_{i}(P_{TIC,i})\times 10^{-6}$ [-] | $1.006$ | $9.93$ | $0.989$ | $0.829$ | $1.000$ | $1.000$ | $0.196$
$n_{major}$ [-] | 4 | 9 | 6 | 7 | 7 | 7 | 20
$n_{minor}$ [-] | 15 | 21 | 30 | 15 | - | - | -
$e_{validation}$ [$\mathrm{m}\mathrm{m}$] | $12.73$ | $47.06$ | $1.63$ | $30.79$ | $53.43$ | $44.39$ | $20.32$
$T_{sim}$ [$\mathrm{s}$] | $52.85$ | $134.75$ | $164.58$ | $142.02$ | $100.63$ | $114.00$ | $272.64$
|  |  |  |  |  |  |

##  5 Conclusions

This work presented a sequential convex program (SCP)method to compute
collision avoidance maneuvers (CAMs) in scenarios involving multiple
encounters between satellites. Both short-term and long-term encounters were
considered. In the first case, the total probability of collision (TPoC)was
used as a collision metric, whereas in the second, we used the total
instantaneous probability of collision (TIPoC). Moreover, for the scenarios in
which uncertainty propagation is needed for long periods, we employed a
Gaussian mixture model (GMM)method to propagate the uncertainty of the
secondary spacecraft and define multiple keep-out zones.

The fuel-optimal CAM optimization problem was cast as a second-order cone
program (SOCP), where different techniques were used to convexify the non-
convex constraints. Two methods were proposed to refine the initial SCP
solution. The first one relies on a DA-based linearization of the TPoC or
TIPoC constraint; the second one updates the probability of collision (PoC)or
instantaneous probability of collision limits associated with the single
conjunctions by means of a non-linear program.

The proposed methods were tested on three scenarios. The first one comprised
multiple consecutive short-term conjunctions: Chan’s PoC method was used to
assess the risk, and the results for the two refinements algorithms were
compared. The second scenario was a double encounter with the same secondary
where the covariance was propagated with a GMM. The last one was a single
long-term encounter studied with the aid of GMMs with different numbers of
mixands. The algorithms can always find an optimal solution for the fuel-
optimal problem, but global optimality is not guaranteed. The SMD
linearization refinement method was deemed more suitable because it behaved
much better in the presence of significant non-linearities of the original
problem.

The proposed method is general, being able to deal with very diverse settings,
including (i) both short- and long-term encounters, (ii) arbitrary dynamics
models, (iii) arbitrary number of encounters, (iv) nonlinear propagation of
the uncertainty of the secondary. The primary’s uncertainty is propagated
linearly when necessary. Since we assume that accurate orbit determination is
performed on the maneuverable satellite, the initial uncertainty of the
primary is almost negligible with respect to the secondary. If the uncertainty
of the primary is larger, a nonlinear propagation model must also be used for
the uncertainty of the primary.

In an ever more cluttered orbital environment, we believe that the complexity
of the CA scenarios for spacecraft will increase dramatically. Given its light
computational burden and proven efficacy and reliability in dealing with
multiple conjunctions, the proposed algorithm can be an enabling tool for the
safe use of Earth’s space in the future.

## Funding Sources

This material is based upon work supported by the Air Force Office of
Scientific Research under award number FA2386-21-1-4115.

## References

  * ESOC [2023] ESOC, “ESA’s Annual Space Environment Report,” Tech. rep., ESA, 2023.
  * Bombardelli and Hernando-Ayuso [2015] Bombardelli, C., and Hernando-Ayuso, J., “Optimal impulsive collision avoidance in low earth orbit,” _Journal of Guidance, Control, and Dynamics_ , Vol. 38, AIAA International, 2015, pp. 217–225.  10.2514/1.G000742.
  * Gonzalo et al. [2019] Gonzalo, J. L., Colombo, C., and Di Lizia, P., “A semi-analytical approach to low-thrust collision avoidance manoeuvre design,” _70th International Astronautical Congress (IAC)_ , Vol. 70, 70th International Astronautical Congress, Washington D.C., 2019, pp. 21–25.
  * Hernando-Ayuso and Bombardelli [2021] Hernando-Ayuso, J., and Bombardelli, C., “Low-thrust collision avoidance in circular orbits,” _Journal of Guidance, Control, and Dynamics_ , Vol. 44, No. 5, 2021, pp. 983–995.  10.2514/1.G005547.
  * Armellin [2021] Armellin, R., “Collision avoidance maneuver optimization with a multiple-impulse convex formulation,” _Acta Astronautica_ , Vol. 186, 2021, pp. 347–362.  10.1016/j.actaastro.2021.05.046.
  * Kim et al. [2012] Kim, E. H., Kim, H. D., and Kim, H. J., “A Study on the collision avoidance maneuver optimization with multiple space debris,” _Journal of Astronomy and Space Sciences_ , Vol. 29, No. 1, 2012, pp. 11–21.  10.5140/JASS.2012.29.1.011.
  * Lee et al. [2019] Lee, K., Park, H., Park, C., and Park, S. Y., “Sub-optimal cooperative collision avoidance maneuvers of multiple active spacecraft via discrete-time generating functions,” _Aerospace Science and Technology_ , Vol. 93, 2019, p. 105298.  10.1016/J.AST.2019.07.031.
  * Frigm et al. [2015] Frigm, R. C., Hejduk, M. D., Johnson, L. C. ., and Plakalovic, D., “Total Probability of Collision as a Metric for Finite Conjunction Assessment and Collision Risk Management,” _Advanced Maui Optical and Space Surveillance Technologies (AMOS) Conference_ , Maui, Hawaii, 2015.
  * Sánchez and Vasile [2021] Sánchez, L., and Vasile, M., “Constrained Optimal Collision Avoidance Manoeuvre Allocation under Uncertainty for Subsequent Conjunction Events,” _72nd International Astronautical Congress, IAC_ , Vol. A6, Dubai, United Arab Emirates, 2021.
  * Masson et al. [2023] Masson, M., Arzelier, D., Joldes, M., Revelin, B., and Thomassin, J., “Multi-maneuver algorithms for multi-risk collision avoidance via nonconvex quadratic optimization,” _IFAC-PapersOnLine_ , Vol. 56, No. 2, 2023, pp. 1995–2000.  10.1016/j.ifacol.2023.10.1094, URL https://linkinghub.elsevier.com/retrieve/pii/S2405896323014970.
  * Patera [2002] Patera, R., “Satellite Collision Probability for Non-Linear Relative Motion,” _AIAA/AAS Astrodynamics Specialist Conference and Exhibit_ , Vol. 26, American Institute of Aeronautics and Astronautics, Reston, Virigina, 2002.  10.2514/6.2002-4632.
  * De Vittori et al. [2023a] De Vittori, A., Dani, G., Di Lizia, P., and Armellin, R., “Numerically Efficient Low-Thrust Collision Avoidance Maneuver Design in GEO Regime with Equinoctial Orbital Elements,” _2nd ESA NEO and Debris Detection Conference_ , January, ESA/ESOC, Darmstadt, Germany, 2023a.
  * De Vittori et al. [2023b] De Vittori, A., Dani, G., Di Lizia, P., and Armellin, R., “Low-Thrust Collision Avoidance Design for Leo Missions With Return To Nominal Orbit,” _33rd AAS/AIAA Space Flight Mechanics Meeting_ , Austin, TX, 2023b, pp. 1–18.
  * De Vittori et al. [2023c] De Vittori, A., Omodei, M., Di Lizia, P., Armellin, R., Gago, P., Torras Ribell, M., Antón, J., and Antón, D., “Numerically Efficient Low-Thrust Fuel-Optimal Collision Avoidance Maneuvers With Tangential Firings,” _AAS/AIAA Astrodynamics Specialist Conference_ , 2023c.
  * Dutta and Misra [2022] Dutta, S., and Misra, A. K., “Convex optimization of collision avoidance maneuvers in the presence of uncertainty,” _Acta Astronautica_ , Vol. 197, 2022, pp. 257–268.  10.1016/j.actaastro.2022.05.038.
  * Mueller [2009] Mueller, J., “Onboard Planning of Collision Avoidance Maneuvers Using Robust Optimization,” _AIAA Infotech@Aerospace Conference_ , American Institute of Aeronautics and Astronautics, Reston, Virigina, 2009.  10.2514/6.2009-2051.
  * Serra et al. [2015] Serra, R., Arzelier, D., Joldes, M., and Rondepierre, A., “Probabilistic Collision Avoidance for Long-term Space Encounters via Risk Selection,” _Advances in Aerospace Guidance, Navigation and Control_ , Springer International Publishing, 2015, pp. 679–698.  10.1007/978-3-319-17518-8_39.
  * Pavanello et al. [2023a] Pavanello, Z., Pirovano, L., and Armellin, R., “Long-Term Encounters Collision Avoidance Maneuver Optimization with a Convex Formulation,” _33rd AAS/AIAA Space Flight Mechanics Meeting_ , Austin, TX, 2023a, pp. 1–20.
  * Pavanello et al. [2023b] Pavanello, Z., Pirovano, L., and Armellin, R., “Fuel-Optimal Collision Avoidance Maneuvers in Long-Term Encounters with Station-Keeping Constraints,” , 2023b.  URL http://arxiv.org/abs/2307.06004.
  * Luo and Yang [2017] Luo, Y.-z., and Yang, Z., “A review of uncertainty propagation in orbital mechanics,” _Progress in Aerospace Sciences_ , Vol. 89, 2017, pp. 23–39.  10.1016/j.paerosci.2016.12.002.
  * Yanez et al. [2019] Yanez, C., Gupta, M., Morand, V., and Dolado, J. C., “On the Gaussianity validity time for orbital uncertainty propagation,” _1st NEO and Debris Detection Conference_ , Vol. 2019-Octob, 2019, pp. 22–24.
  * Vittaldev [2015] Vittaldev, V., “Uncertainty propagation and conjunction assessment for resident space objects,” Ph.D. thesis, The University of Texas at Austin, 2015\.
  * Dutta and Misra [2023] Dutta, S., and Misra, A. K., “Comparison between convex and non-convex optimization methods for collision avoidance maneuvers by a spacecraft,” _Acta Astronautica_ , Vol. 202, No. May 2022, 2023, pp. 900–908.  10.1016/j.actaastro.2022.06.020.
  * Malyuta et al. [2021] Malyuta, D., Yu, Y., Elango, P., and Açıkmeşe, B., “Advances in trajectory optimization for space vehicle control,” _Annual Reviews in Control_ , Vol. 52, 2021, pp. 282–315.  10.1016/j.arcontrol.2021.04.013, URL https://linkinghub.elsevier.com/retrieve/pii/S1367578821000377.
  * Li et al. [2022] Li, J. S., Yang, Z., and Luo, Y. Z., “A review of space-object collision probability computation methods,” _Astrodynamics_ , Vol. 6, No. 2, 2022, pp. 95–120.  10.1007/s42064-021-0125-x.
  * Pirovano and Armellin [2024] Pirovano, L., and Armellin, R., “Detection and estimation of spacecraft maneuvers for catalog maintenance,” _Acta Astronautica_ , Vol. 215, No. December 2023, 2024, pp. 387–397.  10.1016/j.actaastro.2023.12.016.
  * Boyd and Vandenberghe [2004] Boyd, S., and Vandenberghe, L., _Convex Optimization_ , Cambridge University Press, Cambridge, NY, 2004.
  * Malyuta et al. [2022] Malyuta, D., Reynolds, T. P., Szmuk, M., Lew, T., Bonalli, R., Pavone, M., and Açıkmeşe, B., “Convex Optimization for Trajectory Generation: A Tutorial on Generating Dynamically Feasible Trajectories Reliably and Efficiently,” _IEEE Control Systems_ , Vol. 42, No. 5, 2022, pp. 40–113.  10.1109/MCS.2022.3187542.
  * Losacco et al. [2024] Losacco, M., Fossà, A., and Armellin, R., “Low-Order Automatic Domain Splitting Approach for Nonlinear Uncertainty Mapping,” _Journal of Guidance, Control, and Dynamics_ , 2024, pp. 1–20.  10.2514/1.G007271.
  * Bernardini et al. [2023] Bernardini, N., Wijayatunga, M. C., Baresi, N., and Armellin, R., “State-Dependent Trust Region for Successive Convex Optimization of Spacecraft Trajectories,” _33rd AAS/AIAA Space Flight Mechanics Meeting_ , Austin, TX, 2023, pp. 1–20.
  * Alfriend et al. [1999] Alfriend, K. T., Akella, M. R., Frisbee, J., Foster, J. L., Lee, D.-J., and Wilkins, M., “Probability of collision error analysis,” _Space Debris_ , 1999.  https://doi.org/10.1023/A:1010056509803.
  * De Vittori et al. [2022] De Vittori, A., Palermo, M. F., Di Lizia, P., and Armellin, R., “Low-Thrust Collision Avoidance Maneuver Optimization,” _Journal of Guidance, Control, and Dynamics_ , Vol. 45, No. 10, 2022, pp. 1815–1829.  10.2514/1.G006630.
  * Zhang et al. [2020] Zhang, S., Fu, T., Chen, D., and Cao, H., “Satellite instantaneous collision probability computation using equivalent volume cuboids,” _Journal of Guidance, Control, and Dynamics_ , Vol. 43, No. 9, 2020, pp. 1757–1763.  10.2514/1.G004711.
  * Chan [2009] Chan, K., “International Space Station Collision Probability,” _The Aerospace Corporation, Chantilly, VA, USA_ , Vol. 6, No. 2, 2009, pp. 307–314.
  * Alfano [2007] Alfano, S., “Review of conjunction probability methods for short-term encounters,” _Advances in the Astronautical Sciences_ , 2007.
  * Coppola [2012] Coppola, V. T., “Including velocity uncertainty in the probability of collision between space objects,” _AIAA/AAS Astrodynamics Specialist Conference 2014_ , Vol. 143, San Diego, California, 2012, pp. 2159–2178.
  * Alfano [2014] Alfano, S., “Eliminating Assumptions Regarding Satellite Conjunction Analysis,” _The Journal of the Astronautical Sciences_ , Vol. 59, 2014, pp. 676–705.  10.1007/s40295-014-0002-4.
  * Xu and Xiong [2011] Xu, X.-l., and Xiong, Y.-q., “A Research on the Collision Probability Calculation of Space Debris for Nonlinear Relative Motions†,” _Chinese Astronomy and Astrophysics_ , Vol. 35, No. 3, 2011, pp. 304–317.  10.1016/j.chinastron.2011.07.008.
  * Wen and Qiao [2022] Wen, C., and Qiao, D., “Calculating collision probability for long-term satellite encounters through the reachable domain method,” _Astrodynamics_ , Vol. 6, No. 2, 2022, pp. 141–159.  10.1007/s42064-021-0119-8.
  * Chan [2004] Chan, K., “Short-term vs long-term spacecraft encounters,” _Collection of Technical Papers - AIAA/AAS Astrodynamics Specialist Conference_ , Vol. 3, 2004, pp. 1732–1752.
  * Dolado et al. [2011] Dolado, J. C., Legendre, P., Garmier, R., Revelin, B., and Pena, X., “Satellite collision probability computation for long term encounters,” _22nd AAS/AIAA Spaceflight Mechanics Meeting_ , 2011.
  * Núñez Garzón and Lightsey [2022] Núñez Garzón, U. E., and Lightsey, E. G., “Relating Collision Probability and Separation Indicators in Spacecraft Formation Collision Risk Analysis,” _Journal of Guidance, Control, and Dynamics_ , Vol. 45, No. 3, 2022, pp. 517–532.  10.2514/1.G005744.
  * Jones and Doostan [2013] Jones, B. A., and Doostan, A., “Satellite collision probability estimation using polynomial chaos expansions,” _Advances in Space Research_ , Vol. 52, No. 11, 2013, pp. 1860–1875.  10.1016/j.asr.2013.08.027.
  * Adurthi and Singla [2015] Adurthi, N., and Singla, P., “Conjugate unscented transformation-based approach for accurate conjunction analysis,” _Journal of Guidance, Control, and Dynamics_ , Vol. 38, No. 9, 2015, pp. 1642–1658.  10.2514/1.G001027.
  * Vittaldev and Russell [2016] Vittaldev, V., and Russell, R. P., “Space object collision probability using multidirectional Gaussian mixture models,” _Journal of Guidance, Control, and Dynamics_ , Vol. 39, No. 9, 2016, pp. 2161–2167.  10.2514/1.G001610.
  * Armellin et al. [2010] Armellin, R., Di Lizia, P., Bernelli-Zazzera, F., and Berz, M., “Asteroid close encounters characterization using differential algebra: The case of Apophis,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 107, No. 4, 2010, pp. 451–470.  10.1007/s10569-010-9283-5.
  * Wang and Grant [2018] Wang, Z., and Grant, M. J., “Minimum-Fuel Low-Thrust Transfers for Spacecraft: A Convex Approach,” _IEEE Transactions on Aerospace and Electronic Systems_ , Vol. 54, No. 5, 2018, pp. 2274–2290.  10.1109/TAES.2018.2812558.
  * Mao and Acikmese [2021] Mao, Y., and Acikmese, B., “SCvx-fast: A Superlinearly Convergent Algorithm for A Class of Non-Convex Optimal Control Problems,” _ArXiv_ , Vol. 2112, 2021.  URL http://arxiv.org/abs/2112.00108.
  * He and Xu [2018] He, Y., and Xu, M., “Bounded relative orbits in the zonal gravitational field for formation flying,” _Journal of Guidance, Control, and Dynamics_ , Vol. 41, No. 7, 2018, pp. 1629–1641.  10.2514/1.G003107.
  * Morselli et al. [2014] Morselli, A., Armellin, R., Di Lizia, P., and Bernelli Zazzera, F., “A high order method for orbital conjunctions analysis: Sensitivity to initial uncertainties,” _Advances in Space Research_ , Vol. 53, No. 3, 2014, pp. 490–508.  10.1016/j.asr.2013.11.038.