# Characterization of the minimum warning time for low-thrust collision avoidance maneuvers

Lamberto Dell'Elce, Frank de Veld, Jean-Baptiste Pomet

► To cite this version:

Lamberto Dell'Elce, Frank de Veld, Jean-Baptiste Pomet. Characterization of the minimum warning time for low-thrust collision avoidance maneuvers. 2024 AAS/AIAA Astrodynamics Specialist Conference, AAS; AIAA, Aug 2024, Broomfield (Colorado), United States. hal-04871764

HAL Id: hal-04871764

https://inria.hal.science/hal-04871764v1

Submitted on 7 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# CHARACTERIZATION OF THE MINIMUM WARNING TIME FOR LOW-THRUST COLLISION AVOIDANCE MANEUVERS

Lamberto Dell'Elce*, Frank de Veld,† and Jean-Baptiste Pomet‡

Owing to the exponential increase of tracked objects in the LEO and GEO regions, the number of alerts for potential collisions is continuously increasing. Satellites equipped with low-thrust propulsion need to maneuver for a finite, possibly large, time to reduce the risk of collision to a safe threshold. This study offers a characterization of minimal requirements for such time in the context of the fast encounter between a low-thrust actuated satellite and a non-cooperative object. We first introduce a state vector consisting of the time of closest approach and the orbital elements of the primary at this time, and we deduce its equations of motion. Necessary conditions of optimality are then disclosed, and an approximate solution is proposed, which does not require to solve a shooting problem, but it is obtained by means of a backward integration from the safe set for any initial state in the proximity of the advertised collision. A thorough analysis of the optimal synthesis of this solution is proposed with special focus on the computation of the cut and conjugate loci.

## INTRODUCTION

Near-Earth space is becoming increasingly congested with man-made space debris, so that satellites occasionally need to perform collision avoidance maneuvers (CAM) whenever a space debris object is predicted to come dangerously close to their position. Typically, the Joint Space Operations Center (JSpOC) tracks debris objects and sends alerts whenever a potentially dangerous conjunction occurs.¹ This alert does not have to imply a perfect collision, as initial uncertainty in the state of satellites and space debris leads to an increasingly uncertain propagated state. The term time of closest approach (TCA) is used to describe the moment of minimum separation distance between the controlled satellite (primary) and uncontrolled debris object (secondary). The goal of the collision maneuver is to perform a certain control on the satellite in order to increase the minimum separation distance to a safe threshold. While high-thrust CAM on large, conventional satellites are relatively straightforward to approach, low-thrust CAM are more complicated to design due to the small acceleration compared to the drift (uncontrolled motion). Low-thrust CAM still have to be considered, because the usage of low-thrust propulsion satellites has steadily increased in recent years, especially due to the emergence of mega-constellations.

Typically, low-thrust CAM are designed such that fuel consumption is minimized.² Autonomous collision avoidance is desirable, as well. These maneuvers are triggered few days before the predicted conjunction. Starting the maneuver too early can result in a waste of fuel, since the estimation

<footnote>*Researcher, Inria & Université Côte d'Azur, 06902 Valbonne, France</footnote>

<footnote>†PhD researcher, Inria & Université Côte d'Azur, 06902 Valbonne, France</footnote>

<footnote>‡Senior researcher, Inria & Université Côte d'Azur, 06902 Valbonne, France</footnote>
of the closest approach is possibly unreliable because of the large uncertainty in the propagation of the trajectories. Conversely, maneuvering too late may be critical for obvious reasons. An exhaustive study of the minimum time necessary to drive a satellite to a safe state is not yet available in the literature, and it is the object of this paper. Satellites relying on very low-thrust capabilities (e.g., differential drag) may particularly benefit from such analysis.

For this purpose, we introduce a 7-dimensional state vector consisting of the “expected TCA” and the orbital elements of the primary at this time. Appending TCA to the state vector is convenient, because it guarantees that the closest approach is exactly tracked. A detailed derivation of the equations of motion of these variables is offered, which, to our knowledge, is not available in the literature. The resulting system is a time-varying, control-affine system with no drift (i.e., the propagated state at TCA and TCA itself do not change during coasting arcs). We stress that, differently from most research, we do not assume that the relative trajectory is rectilinear when deducing these equations. Any perturbation can be included in the vector field characterizing the drift of the satellites (namely, uncontrolled motion).

After formalizing the problem of minimum warning time to achieve a prescribed safe miss distance, necessary conditions of optimality, stemming from the Pontryagin maximum principle (PMP), are deduced. An approximate solution is then achieved in the context of fast encounters. In this case, we are able to approximate costate variables, and to compute extremal trajectories by backward integrating the equations from the safe set. This solution can be used as an initial guess for the original shooting problem. A major contribution of the manuscript is a thorough analysis of the optimal synthesis of the approximate solution, with special focus on the computation of the cut and conjugate loci.

## MINIMUM WARNING TIME: PROBLEM DEFINITION

Let $x(t)$ and $x_s(t)$ be two sets of orbital elements with values in $M$, open subset of $\mathbb{R}^6$, representing the state of the primary and secondary satellites at time $t$, respectively. Both satellites are supposed to be on nominal, non-controlled trajectories, which are predicted to experience a closed approach at the nominal TCA that we arbitrarily set equal to 0 (choice of the origin of time). Orbital elements of the primary and secondary at the nominal TCA are named $y_0$ and $y_s$, respectively, and are known constants in this study.

We are interested in a non-cooperative CAM, where thrusting capabilities are only onboard the primary satellite, so that motion of $x(t)$ and $x_s(t)$ is governed by

$$ \frac{d x}{d t} = F(x) + \varepsilon G(x) u \quad \text{and} \quad \frac{d x_s}{d t} = F(x_s), $$

where the vector field $F: M \to TM$ is the drift, i.e., uncontrolled orbital motion, which includes the Keplerian contribution and possibly any orbital perturbation*, $G(x)$ can be deduced from Gauss variational equations (GVE), $u$ is the control variable with values on the unit ball $U := \{v \in \mathbb{R}^3 : \|v\| \le 1\}$, and $\varepsilon$ is the thrust-to-mass ratio.

We normalize quantities by choosing length and time units equal to the semi-major axis and the orbital period of the primary at the nominal closest approach, namely $y_{0,1}$ and $t_{orb} = 2\pi\sqrt{y_{0,1}^3/\mu}$ ($\mu = 3.98600 \cdot 10^{14} \text{ m}^3\text{s}^{-2}$ being the Earth's gravitational constant), respectively. Hence, denoting

<footnote>*Time-varying perturbations may be considered as well, although we neglected time dependency to avoid unnecessary notations and additional arguments of the functions.</footnote>
by $T$ the maximum thrust force and by $m$ the mass of the primary, the (dimension-less) definition of $\varepsilon$ is

$$ \varepsilon = 4 \pi^2 \frac{T}{m} \frac{y_{0,1}^2}{\mu} \qquad (1) $$

Focusing on low-thrust maneuvers, which are characterized by negligible fuel consumption, we assume that

**Assumption 1 (Low-thrust).** The parameter $\varepsilon$ is constant throughout the CAM, and $\varepsilon \ll 1$.

We denote by $\varphi^t(\cdot)$ the flow of $F$ (uncontrolled motion) at time $t$. Specifically, $\varphi : \mathbb{R} \times M \to M$ is a smooth function such that, for a fixed $x_0$, $\varphi^t(x_0)$ is the solution at time $t$ of

$$ \frac{d\varphi^t(x_0)}{dt} = F(\varphi^t(x_0)), \quad \varphi^0(x_0) = x_0. $$

The trajectory of the secondary object is fully characterized by $\varphi$ and $y_s$, since $x_s(t) = \varphi^t(y_s)$ for any $t$. The fundamental matrix solution of the uncontrolled motion, which is also referred to as state transition matrix (STM), is defined as $\Phi(t, x) := \frac{\partial \varphi^t(x)}{\partial x}$.

We note that, when using matrix representations, vectors are column, co-vectors are raw, and partial derivatives of an $m$-dimensional vector-valued function with respect to an $n$-dimensional argument are $m \times n$.

**Example 1.** Equinoctial elements are defined as

$$ \begin{aligned} x_1 &= a, & x_2 &= e \sin(\omega + \Omega), & x_3 &= e \cos(\omega + \Omega), \\ x_4 &= \tan \frac{i}{2} \sin(\omega), & x_5 &= \tan \frac{i}{2} \cos(\omega), & x_6 &= \omega + \Omega + M \end{aligned} $$

where $a, e, i, \omega$, and $\Omega$ denote classical elements, and $M$ is the mean anomaly. For these elements, Keplerian drift, flow and STM in matrix form are:

$$ F(x) = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \sqrt{\frac{\mu}{x_1^3}} \end{pmatrix}, \quad \varphi^t(x) = \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6 + \sqrt{\frac{\mu}{x_1^3}} t \end{pmatrix}, \quad \Phi(t, x) = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ -\frac{3}{2} \sqrt{\frac{\mu}{x_1^3}} t & 0 & 0 & 0 & 0 & 1 \end{pmatrix}. $$

The expression of $G(x)$ is given by GVE of equinoctial elements in$^3$[Chap. 10, Pag. 492-493].

Alert occurs because the miss distance (i.e., relative distance at the TCA) is below a prescribed threshold, $\sigma$. Specifically, let $r: M \to \mathbb{R}^3$ and $v: M \to \mathbb{R}^3$ be the smooth maps such that $r(x)$ and $v(x)$ are the Cartesian position and velocity corresponding to the orbital elements $x \in M$ (detailed derivation of these maps for equinoctial elements is available in$^3$[Chap. 10, Pag. 494]). Recall that $y_s$, as defined at the beginning of this section, is the value at time 0 (nominal TCA) of the orbital elements of the secondary, and define the relative position and velocity vectors as the functions

$$ \Delta r(t, x) := r(x) - r(\varphi^t(y_s)), \quad \Delta v(t, x) := v(x) - v(\varphi^t(y_s)), $$
respectively. It holds that

$$ \|\Delta r(0, y_0)\| < \sigma. $$

This study is aimed at characterizing the latest time, $t_0 \le 0$, such that it is possible to execute a CAM guaranteeing that relative distance is never below the safe threshold $\sigma$. The naive mathematical formulation of this concept is

$\min_{t_0 \le 0} |t_0|$ subject to:

$$ \left\{ \begin{array}{l} \displaystyle \frac{d}{dt} x = F(x) + \varepsilon G(x) u \\ \displaystyle \|u(t)\| \le 1 \\ \displaystyle \|\Delta r(t, x)\| \ge \sigma \\ x(t_0) = \varphi^{t_0}(y_0) \end{array} \right. \quad \forall t \in [t_0, \infty) $$

Such formulation does not lend itself to a practical numerical implementation because final time is not bounded. In practice, most fast encounters are characterized by a unique, well-defined close approach, which is a local minimum of the function $t \mapsto \|\Delta r(t, x(t))\|^2$. Multiple encounters (either with the same or with different secondary objects) are possible, but they are beyond the scope of this manuscript. Let $g: \mathbb{R} \times M \to \mathbb{R}^+$ be the smooth function

$$ g(t, x) := \Delta r(t, x) \cdot \Delta v(t, x). $$

For any smooth curve $t \mapsto x(t)$, one has

$$ \frac{1}{2} \frac{d}{dt} (\|\Delta r(t, x)\|^2) = g(t, x). $$

**Assumption 2 (Uniqueness of the closest approach).** There is a unique time $t_{CA} \in [t_0, \infty)$ such that

$$ g(t_{CA}, \varphi^{t_{CA}}(x_0)) = 0, \quad \left. \frac{d}{dt} \left( g(t_{CA}, \varphi^{t_{CA}}(x_0)) \right) \right|_{t=t_{CA}} > 0. $$

**Remark 1.** Relative position and velocity are orthogonal at the closest approach.

By virtue of Assumption 2, we can replace the path constraint $\|\Delta r(t, x(t))\| \ge \sigma$ for all $t$ with a single endpoint condition $\|\Delta r(t_{CA}, x(t_{CA}))\| \ge \sigma$. This condition is active at the optimum (i.e., equality holds), so that the problem can be reformulated as

$\min_{t_0 \le 0} |t_0|$ subject to:

$$ \left\{ \begin{array}{l} \displaystyle \frac{d}{dt} x = F(x) + \varepsilon G(x) u \\ \displaystyle \|u(t)\| \le 1 \\ x(t_0) = \varphi^{t_0}(y_0) \\ \|\Delta r(t_{CA}, x(t_{CA}))\| = \sigma \\ g(t_{CA}, x(t_{CA})) = 0 \end{array} \right. \quad \forall t \in [t_0, t_{CA}] \quad (\text{CAM-X}) $$

Here, time spans on the finite window $[t_0, t_{CA}]$, and both $t_0$ and $t_{CA}$ are design variable to be determined. The last equation is used to characterize $t_{CA}$ itself. Specifically, we stress that the TCA
<figure><img src="image_2.png" /><figcaption>Figure 1. Qualitative illustration of Problem (CAM-X). Here, $t_f$ is the TCA after execution of the CAM.</figcaption></figure>

may move from its advertised value (i.e., $t_{CA} = 0$ if no CAM is executed) as depicted in Figure 1. There is no restriction on the sign of $t_{CA}$, but it is understood that $t_0 < t_{CA}$ because of causality. Initial conditions depend on $t_0$, since they are obtained by backward transportation of $y_0$ with the flow of the drift (red line). Terminal conditions are such that the separation of the satellites at the closest approach matches the safe threshold $\sigma$ (blue curve). The green curve (which is intentionally drawn non-monotone) shows the evolution of the miss distance throughout the maneuver in the sense that it depicts the relative distance evaluated at the conjunction that would occur if thrust is interrupted at some time. This point of view is at the core of the change of variables carried out in the next section.

Although Problem (CAM-X) exhibits some uncommon features like variable initial and final times, and time-dependent boundary conditions, it can indeed be solved by means of direct techniques. Nonetheless, we aim at the development of a reliable approximation and a through understanding of the solution.

**Remark 2.** When dealing with “fast encounters”, variations of the TCA due to the CAM are extremely small (i.e., $t_{CA} \approx 0$). The approximate solution that we will develop requires an additional assumption regarding the definition of fast encounter. We were cautious to avoid the introduction of such assumption at this stage because the next two sections are of general purpose, and may be used to tackle problems where the slow-fast nature of the encounter is less clear. For example, no assumption on the relative trajectory was introduced (other than the lack of control capabilities onboard the secondary), and the order of magnitude of $\sigma$ was not yet discussed.

## EQUATIONS OF MOTION FOR THE STATE AT THE CLOSEST APPROACH

Motivated by the uniqueness of the closest approach outlined in Assumption 2, we introduce a new set of state variables which aims at tracking the evolution of the closest approach as control is applied. These variables consists of the instantaneous TCA, namely $t_{CA}(t)$, and the transportation of the state of the primary at TCA by the flow of the drift (uncontrolled motion) $\varphi$, namely

$$y(t) := \varphi^{t_{CA}-t}(x(t)). \qquad (2)$$

Intuitively, $y(t)$ and $t_{CA}(t)$ describe the closest approach that will occur if thrust is stopped at time $t$. In particular, because $t_0$ is the moment when the CAM is started, it holds that $t_{CA}(t_0) = 0$ and
$y(t_0) = \varphi^0(y_0) = y_0$.

At the closest approach, the time derivative of the relative distance vanishes (this is how the closest approach is characterized in Assumption 2). Now, to keep track of the closest approach, we impose that this constraint is identically satisfied by the new variables, namely

$$g(t_{CA}(t), y(t)) \equiv 0 \quad \forall t \in \mathbb{R} \qquad (3)$$

In turn, the orthogonality between relative position and velocity is used to characterize the TCA as a function of time. Equations of motion are obtained by differentiating in time $g(t_{CA}, y) \equiv 0$:

$$\frac{\mathrm{d}g}{\mathrm{d}t}(t_{CA}, y) \equiv 0 \quad \Rightarrow \quad \frac{\partial g}{\partial t_{CA}} \frac{\mathrm{d}t_{CA}}{\mathrm{d}t} + \frac{\partial g}{\partial y} \frac{\mathrm{d}y}{\mathrm{d}t} = 0, \qquad (4)$$

and by differentiating in time the definition of $y$ in Eq. (2):

$$\begin{aligned} \frac{\mathrm{d}y}{\mathrm{d}t} = \frac{\mathrm{d}\varphi^{t_{CA}-t}(x)}{\mathrm{d}t} &= \frac{\partial \varphi^{t_{CA}-t}}{\partial t} + \frac{\partial \varphi^{t_{CA}-t}}{\partial t_{CA}} \frac{\mathrm{d}t_{CA}}{\mathrm{d}t} + \frac{\partial \varphi^{t_{CA}-t}}{\partial x} \frac{\mathrm{d}x}{\mathrm{d}t} \\ &= -F(y) + F(y) \frac{\mathrm{d}t_{CA}}{\mathrm{d}t} + \Phi(t_{CA}-t, x) [F(x) + \varepsilon G(x)u] \end{aligned}$$

Defining (in analogy with the D-matrix introduced in⁴)

$$D(\tau, y) := \Phi(\tau, \varphi^{-\tau}(y)) G(\varphi^{-\tau}(y)),$$

and noting that $F(y) = \Phi(t_{CA}-t, x)F(x)$ because of the invariance of a vector field under its own flow, yields

$$\frac{\mathrm{d}y}{\mathrm{d}t} = \varepsilon D(t_{CA}-t, y) u + F(y) \frac{\mathrm{d}t_{CA}}{\mathrm{d}t} \qquad (5)$$

The equations of motion of $y$ and $t_{CA}$ are finally obtained by combining Eq. (4) and (5):

$$\frac{\mathrm{d}y}{\mathrm{d}t} = \varepsilon B_y(t, t_{CA}, y) u, \quad \frac{\mathrm{d}t_{CA}}{\mathrm{d}t} = -\varepsilon B_{t_{CA}}(t, t_{CA}, y) u, \qquad (6)$$

with $B_y : \mathbb{R} \times \mathbb{R} \times M \to \mathbb{R}^{6 \times 3}$ and $B_{t_{CA}} : \mathbb{R} \times \mathbb{R} \times M \to \mathbb{R}^{1 \times 3}$ defined as

$$B_y(t, t_{CA}, y) := \left[ I - \frac{1}{g'(t_{CA}, y)} F(y) \frac{\partial g}{\partial y} \right] D(t_{CA} - t, y),$$

$$B_{t_{CA}}(t, t_{CA}, y) := - \frac{1}{g'(t_{CA}, y)} \frac{\partial g}{\partial y} D(t_{CA} - t, y),$$

where $I$ denotes the (6-dimensional) identity matrix, and*

$$g'(t_{CA}, y) := \frac{\partial g}{\partial y} F(y) + \frac{\partial g}{\partial t_{CA}}. \qquad (7)$$

Note that, as a control system with state $(y, t_{CA})$, Eq. (6) is non controllable, since $g(t_{CA}, y)$ is constant along solutions regardless the control. In fact, we choose to stay on the sub-manifold $\{(t_{CA}, y) \text{ s.t. } g(t_{CA}, y) = 0\}$. These equations are nonlinear, time-varying, driftless (y and $t_{CA}$

<footnote>*We recall that, according to the notation adopted in this paper, $F(y) \frac{\partial g}{\partial y}$ is 6 × 6, whereas $\frac{\partial g}{\partial y} F(y)$ is a scalar function.</footnote>
only change when thrusting, this is true because the definition of y cancels the drift), and control affine. They do not depend on a specific control law, and may include any drift (although it is convenient to have an explicit expression for its flow and STM). In this study, they are paramount to achieve the approximate solution of Problem (CAM-X) and the optimal synthesis disclosed in the reminder of the manuscript. However, we believe that their exploitation may go beyond the scope of this paper, namely this local closest-approach-oriented point of view:

* Paves the way to the development of feedback laws for autonomous CAM that increase instantaneous rate of change of the miss distance;

* Offers a tool to deal with encounters of unclear status (slow-fast), because there is no approximation of the relative trajectory. The natural motion of both the primary and secondary objects is governed by the vector field $F$;

* Facilitates the treatment of multiple encounters. In this case, instead of discretizing the entire trajectory or including path constraints, it is sufficient to consider a state vector obtained by concatenating the state of each close approach.

Using the new coordinates, and Assumption 2, Problem (CAM-X) is recast into

$$ \begin{array}{ll} \min_{t_0 \le 0} |t_0| & \text{subject to:} \\ & \left\{ \begin{array}{l} \frac{dy}{dt} = \varepsilon B_y(t, t_{CA}, y) u \\ \frac{dt_{CA}}{dt} = \varepsilon B_{t_{CA}}(t, t_{CA}, y) u \\ |u(t)| \le 1 \end{array} \right. \quad \forall t \in [t_0, t_f] \\ & y(t_0) = y_0, \qquad \|\Delta r(t_{CA}(t_f), y(t_f))\| = \sigma, \\ & t_{CA}(t_0) = 0, \qquad t_{CA}(t_f) = t_f. \end{array} \quad (\text{CAM-Y}) $$

Final time, $t_f$, is free, and the boundary conditions are intuitively illustrated by the departure and arrival points of the green line in Figure 1 (i.e., red and blue dots, respectively).

## NECESSARY CONDITIONS FOR OPTIMALITY

Consider a generic “Mayer problem”, i.e. the minimization problem

$$ \min \Psi(t_0, z(t_0), t_f, z(t_f)), $$

among real numbers $t_0 < t_f$ and solutions of

$$ \frac{dz}{dt} = f(t, z, u), \quad z \in N, u \in U, $$

where $N$ is an open subset of $\mathbb{R}^n$ and $U$ a compact subset of $\mathbb{R}^m$, satisfying initial/terminal constraints

$$ \psi(t_0, z(t_0), t_f, z(t_f)) = 0, $$

where $f: \mathbb{R} \times N \to TN$, $\Psi: \mathbb{R} \times N \times \mathbb{R} \times N \to \mathbb{R}$, and $\psi: \mathbb{R} \times N \times \mathbb{R} \times N \to \mathbb{R}^\ell$ ($\ell$ is the number of constraints) are sufficiently smooth. We recall Pontryagin's necessary condition for optimality in this problem; the following form may be found, e.g., in [5, Thm 4.2.i and eqn. (4.4.3)] (up to the sign of $\lambda_0$ and max vs. min in (11)), see also.⁶⁻⁹
**Proposition 1** (Pontryagin Maximum Principle). Define the Hamiltonian $\mathcal{H}$ as the following function of time $t$, the state $z$, the control $u$ and the adjoint state $\lambda \in (\mathbb{R}^n)^*$:

$$ (t, \lambda, z, u) \mapsto \mathcal{H}(t, \lambda, z, u) = \langle \lambda, F(t, z, u) \rangle, \qquad (8) $$

where $(\mathbb{R}^n)^*$ stands for the dual of $\mathbb{R}^n$ and $\langle \lambda, v \rangle$, $v \in \mathbb{R}^n$ for the linear form $\lambda$ applied to the vector $v$, and the maximized Hamiltonian

$$ \mathcal{H}^*(t, \lambda, z) = \max_{u \in U} \mathcal{H}(t, \lambda, z, u). \qquad (9) $$

If $(\overline{t_0}, \overline{t_f}, \overline{z}(\cdot), \overline{u}(\cdot))$ is a solution of the above optimal control problem, then there exist a constant $\lambda_0$ and a function $t \mapsto \lambda(t)$, $[t_0, t_f] \to (\mathbb{R}^N)^*$ such that

1. $\lambda_0 \le 0$ and $(\lambda_0, \lambda(t)) \ne (0, 0)$ for a.e. $t$,

2. $\lambda(\cdot)$ is a solution of (a.e. $t$):

$$ \frac{d \lambda}{dt} = - \frac{\partial \mathcal{H}}{\partial z} \qquad (10) $$

3. For almost all $t$, with $\mathcal{H}^*$ given by (9),

$$ \mathcal{H}(t, \lambda(t), \bar{z}(t), \bar{u}(t)) = \mathcal{H}^*(t, \lambda(t), \bar{z}(t)). \qquad (11) $$

4. There are constant numbers $\nu_1, \dots, \nu_\ell$ such that the following initial/terminal condition holds*:

$$ \lambda_0 \left( \begin{array}{c} \frac{\partial \Psi}{\partial t_0} \\ \frac{\partial \Psi^\top}{\partial z_0} \\ \frac{\partial \Psi}{\partial t_f} \\ \frac{\partial \Psi^\top}{\partial z_f} \end{array} \right) + \left( \begin{array}{c} -\mathcal{H}(t_0) \\ \lambda(t_0)^\top \\ \mathcal{H}(t_f) \\ -\lambda(t_f)^\top \end{array} \right) = \sum_{i=1}^\ell \nu_i \left( \begin{array}{c} \frac{\partial \psi_i}{\partial t_0} \\ \frac{\partial \psi_i^\top}{\partial z_0} \\ \frac{\partial \psi_i}{\partial t_f} \\ \frac{\partial \psi_i^\top}{\partial z_f} \end{array} \right), \qquad (12) $$

where the argument of all partial derivatives is $(\overline{t_0}, \overline{z}(\overline{t_0}), \overline{t_f}, \overline{z}(\overline{t_f}))$, the short-hand notation $\mathcal{H}(t)$ stands for $\mathcal{H}(t, \lambda(t), \bar{z}(t), \bar{u}(t))$ (and could be replaced with $\mathcal{H}^*(t)$ or $\mathcal{H}^*(t, \lambda(t), \bar{z}(t))$).

**Remark 3.** These conditions are positively (because of the sign condition on $\lambda_0$) homogeneous with respect to $(\lambda_0, \lambda(\cdot), \nu_1, \nu_\ell)$, i.e., they are unchanged upon multiplying these by a positive constant. This allows a normalization.

We may now apply this proposition to state the first-order optimality necessary condition for the

<footnote>*A more geometric statement, without multipliers, is to require that the transpose of the right-hand side of (12), seen as a linear form $\mathbb{R}^{2n+2} \to \mathbb{R}$ be zero in restriction to the subspace $T(\overline{t_0}, \overline{z}(\overline{t_0}), \overline{t_f}, \overline{z}(\overline{t_f}))\mathcal{B}$, tangent space at point $(\overline{t_0}, \overline{z}(\overline{t_0}), \overline{t_f}, \overline{z}(\overline{t_f}))$ to $\mathcal{B} := \{(t_0, z_0, t_f, z_f) \in \mathbb{R} \times \mathbb{N} \times \mathbb{R} \times \mathbb{N}, \psi(t_0, z_0, t_f, z_f) = 0\}$, which is a sub-manifold of $\mathbb{R} \times \mathbb{N} \times \mathbb{R} \times \mathbb{N}$ with dimension $2n+2-\ell$.</footnote>
optimal control problem posed in Eq. (CAM-Y), that is of the form above with

$$n = 7, \quad m = 3, \quad l = 9, \quad z = \begin{pmatrix} y \\ t_{CA} \end{pmatrix},$$

$$\Psi(t_0, z(t_0), t_f, z(t_f)) = -t_0,$$

$$f(t, z, u) = \varepsilon B(t, z) u, \quad \text{with } B(t, z) := \begin{pmatrix} B_y(t, t_{CA}, y) \\ B_{t_{CA}}(t, t_{CA}, y) \end{pmatrix},$$

$$\psi(t_0, z(t_0), t_f, z(t_f)) = \begin{pmatrix} \|\Delta r(z(t_f))\|^2 - \sigma^2 \\ \frac{2}{t_{CA}(t_f) - t_f} \\ y(t_0) - y_0 \\ t_{CA}(t_0) - 0 \end{pmatrix}.$$

where $y(\cdot)$ and $t_{CA}(\cdot)$ are now understood as short-hand notations for the canonical projections of $z(\cdot)$, into $M$ and $\mathbb{R}$, respectively, and we introduced the innocent abuse of notation $\Delta r(z) = \Delta r(t_{CA}, y)$ (extensively used in the remainder of the manuscript).

Hence, the Hamiltonian of our problem is, according to Eq. (8):

$$\mathcal{H}(t, \lambda, z, u) = \varepsilon \lambda B(t, z) u$$

Maximization yields

$$\bar{u} = \left( \frac{\lambda B(t, z)}{\|\lambda B(t, z)\|} \right)^{\top} \quad \Rightarrow \quad \mathcal{H}^*(t, \lambda, z) = \varepsilon \|\lambda B(t, z)\|. \qquad (13)$$

Dynamics of the adjoint variables are, according to Eq. (10),

$$\frac{d \lambda_y}{d t} = -\varepsilon \left\langle \lambda, \frac{\partial B_y}{\partial y}(t, t_{CA}, y) \right\rangle u, \quad \frac{d \lambda_{t_{CA}}}{d t} = -\varepsilon \left\langle \lambda, \frac{\partial B_{t_{CA}}}{\partial t_{CA}}(t, t_{CA}, y) \right\rangle u,$$

where $\lambda$ is partitioned into $\lambda = (\lambda_y, \lambda_{t_{CA}})$ in accordance with $z = (y^\top, t_{CA})^\top$.

Now, the transversality conditions outlined in Eq. (12) read

$$\begin{pmatrix} -\lambda_0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} -\mathcal{H}^*(t_0) \\ \lambda_y(t_0)^\top \\ \lambda_{t_{CA}}(t_0) \\ \mathcal{H}^*(t_f) \\ -\lambda_y(t_f)^\top \\ -\lambda_{t_{CA}}(t_f) \end{pmatrix} = \begin{pmatrix} 0 \\ (\nu_3, \dots, \nu_8)^\top \\ \nu_9 \\ -\nu_2 \\ \nu_1 \left( \frac{\partial r}{\partial y}(y(t_f)) \right)^\top \Delta r(t_f, y(t_f)) \\ -\nu_1 \Delta r(t_f, y(t_f))^\top v (\varphi^{t_f}(y_s)) + \nu_2 \end{pmatrix}.$$

The first nine lines may be used to eliminate $\lambda_0, \nu_2, \dots, \nu_9$ ($\lambda_0 \le 0$ is implied by $\lambda_0 = -\mathcal{H}(t_0)$ and Eq. (13)). Substituting $\nu_2 = -\mathcal{H}^*(t_f)$, the last seven lines and first two equations of $\psi$ yield the
shooting problem,

$$ (t_0, t_f, \lambda(t_0), \nu_1) \mapsto \left\{ \begin{array}{l} \|\Delta r(z(t_f))\|^2 = \sigma^2 \\ t_{CA}(t_f) = t_f \\ \lambda_y(t_f) = -\nu_1 \Delta r(z(t_f))^{\top} \frac{\partial r}{\partial y}(y(t_f)) \qquad (\text{SHOOT}) \\ \lambda_{t_{CA}}(t_f) - \mathcal{H}^*(t_f) = \nu_1 \Delta r(z(t_f))^{\top} v(\varphi^{t_f}(y_s)) \\ \nu_1 \in \{-\sigma^{-1}, 0, \sigma^{-1}\} \end{array} \right. $$

Here, $(t_0, t_f, \lambda(t_0), \nu_1)$ are referred to as shooting variables, and trajectories of $z(t)$ are emanated from the Hamiltonian in Eq. (13) with initial conditions $(y_0, 0, \lambda(t_0))$. Solving necessary conditions consists in finding a set of shooting variables solution of Problem (SHOOT), which is a zero search of a 10-dimensional function. The last equation $(\nu_1 \in \{-\sigma^{-1}, 0, \sigma^{-1}\})$ follows from homogeneity of the Hamiltonian (see Remark 3). The magnitude $\sigma^{-1}$ was arbitrarily chosen in view of some simplifications (recall that $\|\Delta r(z(t_f))\| = \sigma$). In the next session, we show that the cases $\nu_1 = 0$ and $\nu_1 = +\sigma^{-1}$ do not yield any feasible solution for fast encounters.

## APPROXIMATE SOLUTION FOR FAST ENCOUNTERS

We now develop an approximate solution of Problem (CAM-Y) consisting of an approximation of the shooting variables of Problem (SHOOT). This solution is valuable for fast encounters, namely

**Assumption 3 (Fast encounter).** The miss distance and relative velocity at the nominal TCA are of order

$$ \|\Delta r(0, y_0)\| \le \sigma \ll 1, \quad \|\Delta v(0, y_0)\| = \mathcal{O}(1). $$

This assumption comes intentionally late in the manuscript to emphasize that previous developments are general and may be used to tackle problem where the status (slow/fast) of the encounter is not well defined.

**Remark 4.** There are now two small variables, namely $\varepsilon$ and $\sigma$, independent and (possibly) of different magnitude. In real-life operational scenarios, a miss distance larger than 10 kilometers is often considered safe, which fixes a reasonable order of magnitude for $\sigma$ ranging from $10^{-3}$ in low-Earth orbit (LEO) to $10^{-4}$ in geostationary orbit (GEO). Hence, the orbit of the primary will be modified to accommodate displacements of the Cartesian position of only few kilometers, which is a small distance with respect to the size of Earth orbits. This legitimates the forthcoming Assumption 4.

Inspection of the boundary conditions of Problem (CAM-Y) suggests that $\varepsilon|t_0| = \mathcal{O}(\sigma) \ll 1$, because $\|\Delta r(t_{CA}(t_f), y(t_f))\| - \|\Delta r(0, y_0)\| \le \sigma$ so that $z(t_f) - z(0) = \mathcal{O}(\sigma)$ (given boundedness of partial derivatives of $\|\Delta r(z)\|$). For this reason, even in the case of multi-revolution CAM, variations of the state are bounded by $\sigma \ll 1$, which motivates the following assumption.

**Assumption 4 (Neighboring).** The state remains $\sigma$-close to the nominal conjunction during the CAM, namely

$$ z(t) = z_0 + \sigma\delta z(t) \quad \forall t \in [t_0, t_f], $$

<footnote>*This estimate neglects that $\Delta r$ is a map from $N \subset \mathbb{R}^7$ to $\mathbb{R}^3$ (so that small variations of $\Delta r$ could be generated with large variations of $z$), but this possibility is definitely hindered by Assumption 4.</footnote>
where $z_0 = \begin{pmatrix} y_0 \\ 0 \end{pmatrix}$, $\delta z(t_0) = 0$, and $\delta z(t)$ is at most of order 1 for all $t \in [t_0, t_f]$.

The next lemma formalizes the consequences that the bare introduction of Assumptions 3 and 4 entails on solutions of Problem (SHOOT).

**Remark 5.** Because $0 \le \|\Delta r(z(t))\| \le \sigma$ during the entire maneuver, it is not possible to linearize the endpoint constraint $\|\Delta r(z(t))\|^2 = \sigma^2$, cause the two leading terms of the expansion have the same order of magnitude of the second-order term, namely

$$ \|\Delta r(z_0 + \sigma \delta z)\|^2 = \underbrace{\|\Delta r(z_0)\|^2}_{\mathcal{O}(\sigma^2)} + 2 \underbrace{\sigma \Delta r(z_0)^\top}_{\mathcal{O}(\sigma^2)} \frac{\partial \Delta r}{\partial z} \delta z + \sigma^2 \underbrace{\left\| \frac{\partial \Delta r}{\partial z} \delta z \right\|^2}_{\text{non negligible}} + \mathcal{O}(\sigma^3). $$

As such, truncation of the quadratic term of the expansion is not possible.

**Lemma 2.** Under Assumptions 1-4:

1. The relative position at the end of the maneuver is close to the circle of radius $\sigma$ orthogonal to the B-plane of the nominal closest approach in the sense:

$$ \frac{\Delta r(z(t_f))}{\sigma} = \hat{n}(\theta) + \mathcal{O}(\sigma), \quad \hat{n}(\theta) \in \mathbb{S}^1, \quad \hat{n}(\theta) \cdot \Delta v(z_0) = 0, \quad (14) $$

for some angle $\theta$.

2. No feasible solution of Problem (SHOOT) exists for $\nu_1 = 0$ and $\nu_1 = +\sigma^{-1}$. For $\nu_1 = -\sigma^{-1}$,

$$ \lambda(t_f) = \hat{n}(\theta)^{\top} \frac{\partial \Delta r}{\partial z}(z_0) + \mathcal{O}(\varepsilon + \sigma) \quad (15) $$

for the same $\hat{n}(\theta)$ in Eq. (14).

3. Approximating the final time with $t_f = 0$ (and neglecting the endpoint constraint $t_{CA}(t_f) = t_f$), introduces a negligible error ($\sigma$-small) in the evaluation of endpoint constraints.

*Proof.* 1. Expanding the defining condition of the closest approach (Eq. (3)), and leveraging on Assumption 3, we get that, for all $t$ in $[t_0, t_f]$,

$$ \begin{align*} g(z(t)) \equiv 0 \quad &\Rightarrow \quad \frac{\Delta r(z(t))}{\sigma} \cdot \Delta v(z_0 + \sigma \delta z(t)) = 0 \\ &\Rightarrow \quad \frac{\Delta r(z(t))}{\sigma} \cdot (\Delta v(z_0) + \mathcal{O}(\sigma)) = 0 \\ &\Rightarrow \quad \frac{\Delta r(z(t))}{\sigma} \cdot \Delta v(z_0) = \mathcal{O}(\sigma) \end{align*} $$

which, in particular, is true at time $t_f$, when $\|\Delta r(z(t_f))\| = \sigma$.

2. About the case $\nu_1 = 0$. Substituting $\lambda_y(t_f) = 0$ and $\lambda_{t_{CA}}(t_f) - \mathcal{H}^*(t_f) = 0$ in Eq. (13) yields (recall that, at time $t = t_f$, $t - t_{CA}(t)$ is zero):

$$ \mathcal{H}^*(t_f) = \varepsilon \frac{\mathcal{H}^*(t_f)}{|g'(t_{CA}(t_f), y(t_f))|} \left\| \frac{\partial g}{\partial y} D(0, y(t_f)) \right\|, $$
where partial derivatives are evaluated at $(y(t_f), t_f)$. Since $H^*(t_f) \neq 0$ (if not, the whole adjoint vector $(\lambda_0, \lambda_y, \lambda_{t_{CA}})$ would be zero at time $t_f$), substitution of the definition of $g'$ given in Eq. (7) gives

$$ \varepsilon \left\| \frac{\partial g}{\partial y} D(0, y(t_f)) \right\| = \left| \frac{\partial g}{\partial y} F(y(t_f)) + \frac{\partial g}{\partial t_{CA}} \right| . \quad (16) $$

Substituting the definition of $g$ gives for the right-hand term:

$$ \begin{aligned} \frac{\partial g}{\partial y} F(y) + \frac{\partial g}{\partial t_{CA}} &= \Delta v(z) \cdot \left( \underbrace{\frac{\partial r}{\partial y}(y) F(y)}_{v(y)} - \underbrace{\frac{\partial r}{\partial y}(\varphi^t(y_s)) F(\varphi^t(y_s))}_{v(\varphi^t(y_s))} \right) + \mathcal{O}(\sigma) \\ &= \|\Delta v(z)\|^2 + \mathcal{O}(\sigma) \end{aligned} $$

Because of Assumptions 1 and 3, the left and right-hand terms of Eq. (16) are of different orders of magnitude ($\varepsilon \ll 1$ and 1, respectively), so that the equality cannot be satisfied.

About the case $\nu_1 = \pm \sigma^{-1}$. First, we note that $H^*(t_f)/\|\lambda(t_f)\| = \mathcal{O}(\varepsilon)$. Hence, the term $H^*(t_f)$ in the fourth line of Problem (SHOOT) can be neglected, so that endpoint conditions on $\lambda$ become

$$ \begin{aligned} \lambda(t_f) &= -\nu_1 \Delta r(z(t_f))^{\top} \frac{\partial \Delta r}{\partial z}(z(t_f)) + \mathcal{O}(\sigma \varepsilon) \\ &= -\nu_1 \sigma \hat{n}(\theta)^{\top} \frac{\partial \Delta r}{\partial z}(z_0) + \mathcal{O}(\sigma \varepsilon + \sigma^2), \end{aligned} \quad (17) $$

where the last approximation holds by linearizing the partial derivative of $\Delta r$ and replacing $\Delta r(z(t_f)) = \sigma \hat{n}(\theta) + \mathcal{O}(\sigma^2)$ (one order of sigma will be simplified when replacing $\nu_1$). Second, consider the time derivative of the relative position vector

$$ \frac{d \Delta r}{dt} = \frac{\partial \Delta r}{\partial z} \frac{dz}{dt} \approx -\nu_1 \frac{\partial \Delta r}{\partial z} B(t, z) \frac{B(t, z)^{\top} \frac{\partial \Delta r(z_0)}{\partial z} \hat{n}(\theta)}{\left\| B(t, z)^{\top} \frac{\partial \Delta r(z_0)}{\partial z} \hat{n}(\theta) \right\|}, \quad (18) $$

where we used Eq. (17) in the expression of the optimal control. Note that the first partial of $\Delta z$ (without explicit argument) is evaluated at the current $z(t)$. At the final time, using the expression of $\Delta r(z(t_f)) \approx \hat{n}(\theta)$ proved above, it holds

$$ \Delta r(z(t_f)) \cdot \frac{d \Delta r(z(t_f))}{dt} \approx -\nu_1 \left\| B(t_f, z(t_f))^{\top} \frac{\partial \Delta r(z(t_f))^{\top}}{\partial z} \hat{n}(\theta) \right\|, $$

which is positive when $\nu_1 = -\sigma^{-1}$. Hence, the only feasible solutions are obtained for $\nu_1 = -\sigma^{-1}$, which guarantees that the flow is exiting the unsafe region. The expression of $\lambda(t_f)$ follows from substitution of $\nu = -\sigma^{-1}$ in Eq. (17).

3. Assumption 4 entails that $t_f = \sigma \delta t_{CA}(t_f) = \mathcal{O}(\sigma)$. A $\sigma$-small variation of $t_f$ on the $\varepsilon$-slow dynamics of $z$ yields a variation of order $\sigma \varepsilon$ in the final state, so that $\delta z(t_f) = \delta z(0) + \mathcal{O}(\varepsilon)$. Substituting this estimate in the constraints on $\|\Delta r(z(t_f))\|$ and $\lambda(t_f)$, and expanding with respect to $\sigma \delta z$ gives the result.

□
From now on, we solve an approximate problem, where the first point of Lemma 2 suggests that relative position at the conjunction evolves in proximity of the the disk of radius $\sigma$, orthogonal to $\Delta v(z_0)$. Hence, the endpoint condition $\|\Delta r(z(t_f))\| = \sigma$, which was supposed to constrain the final state on $\mathbb{S}^2$, is (at first order in $\sigma$) a constraint on $\mathbb{S}^1$, instead. The second point of the lemma offers an approximate expression for $\lambda(t_f)$. We will show below that it can be used to approximate the shooting variable $\lambda(t_0)$, as well. In addition, an approximation for $t_f$ is given ($t_f \approx 0$).

The next step of our approximation consists in linearizing the equations of motion (Eq. (6)). This is achieved by neglecting terms of order $\sigma\varepsilon$. The most natural linearization is possibly

$$ \frac{d\delta z}{dt} = \frac{\varepsilon}{\sigma} B(t, z_0) u. $$

However, this approximation does not guarantee that motion strictly evolves in the B-plane defined by $\Delta v(z_0)$. Instead, we consider the following dynamics:

$$ \begin{aligned} \frac{d\delta y}{dt} &= \frac{\varepsilon}{\sigma} \left[ I - F(y_0) \frac{\Delta v(z_0)^\top}{\|\Delta v(z_0)\|^2} \frac{\partial r}{\partial y} \right] D(-t, y_0) u \\ \frac{d\delta t_{CA}}{dt} &= -\frac{\varepsilon}{\sigma} \frac{\Delta v^\top}{\|\Delta v(z_0)\|^2} \frac{\partial r}{\partial y} D(-t, y_0) u \end{aligned} \quad (19) $$

where partials of $r$ are evaluated at $y_0$. These equations were obtained from Eq. (6) by noting that $g(t, y)' = \|\Delta v(z_0)\|^2 + \mathcal{O}(\sigma)$ (see first point of the proof of Lemma 2), and that

$$ \frac{\partial g}{\partial y} = \Delta v(z_0)^\top \frac{\partial r}{\partial y}(y_0) + \mathcal{O}(\sigma). $$

It is straightforward to verify that with these equations

$$ \Delta v(z_0) \cdot \frac{d\Delta r}{dt}(z_0 + \sigma\delta z) \equiv 0, $$

so that, because $g(z_0) = \Delta v(z_0) \cdot \Delta r(z_0) = 0$, they guarantee that $\Delta v(z_0) \cdot \Delta r(t) \equiv 0$ for all $t$.

Equation (19) is a linear, time-varying, drift-less control system. The last property entails major structural differences with respect to its nonlinear counterpart:

1. Because the right-hand side of the equations does not depend on $\delta z$, the adjoint vector $\lambda$ is constant in time. Specifically, $\lambda(t_0) = \lambda(t_f)$, so that Eq. (15) is a valuable approximation not only of $\lambda(t_f)$, but also for the shooting variable $\lambda(t_0)$.

2. The ordinary differential equation (ODE) becomes a simple quadrature

$$ \delta z(t) = \frac{\varepsilon}{\sigma} \int_{t_0}^{t} \tilde{B}(t, z_0) u(t, \lambda(\theta)) dt, $$

where the definition of $\tilde{B}$ follows immediately from Eq. (19), $\lambda(\theta)$ denotes expression of $\lambda$ given in Eq. (15) (we emphasize that it depends on a single parameter related to the final relative position), and $u(t, \lambda) = \|\lambda \tilde{B}(t, z_0)\|^{-1} (\lambda \tilde{B}(t, z_0))^\top$.
3. Dynamics is decoupled. For the evaluation of the constraint $\Delta r(0) = \hat{n}(\theta)$ it is sufficient to integrate the relative position at conjunction (this expression is obtained replacing the new equations of motion in Eq. (18)):

$$ \frac{d \Delta r}{d t} = \varepsilon U(t) \frac{U(t)^{\top} \hat{n}(\theta)}{\|U(t)^{\top} \hat{n}(\theta)\|}, \qquad (20) $$

where

$$ U(t) := \left( I - \frac{\Delta v(z_0) \Delta v(z_0)^{\top}}{\|\Delta v(z_0)\|^2} \right) \frac{\partial r}{\partial y}(y_0) D(-t, y_0). $$

**Remark 6.** The only parameters that are needed to evaluate $U(t)$ (and then to integrate Eq (20)) are $y_0$, the direction of $\Delta v(z_0)$, and $\varepsilon$.

At this point, only $\theta$ and $t_0$ need to be determined to complete the approximation. Specifically, we need to solve the reduced shooting problem (dimension 2):

$$ (t_0, \theta) \mapsto \begin{cases} \Delta r(z_0 + \sigma \delta z(0)) \cdot \hat{b}_1 = \sigma \cos(\theta) \\ \Delta r(z_0 + \sigma \delta z(0)) \cdot \hat{b}_2 = \sigma \sin(\theta) \end{cases} \qquad (\text{RED-SHOOT}) $$

where $\Delta r$ is integrated with Eq. (20) with the initial condition $\Delta r(z_0)$ at time $t_0$. Here, $\hat{b}_1$ and $\hat{b}_2$ denote two orthogonal directions in the B-plane defined by the direction of $\Delta v(z_0)$. We arbitrarily chose

$$ \hat{b}_1 \parallel \left( I - \frac{\Delta v(z_0) \Delta v(z_0)^{\top}}{\|\Delta v(z_0)\|^2} \right) v(y_0), \quad \hat{b}_2 = \hat{b}_3 \times \hat{b}_1, $$

so that $\hat{b}_1$ is in the direction of the projection of $v(y_0)$ in the B-plane at the conjunction (this choice is indeed not suitable to treat the case $\Delta v(z_0) \parallel v(z_0)$).

Next section is devoted to the solution of Problem (RED-SHOOT). With this solution, the approximation of the shooting variables (which may serve as an initial guess for solving Problem (SHOOT)) is:

* $t_0$ obtained by solving Problem (RED-SHOOT),

* $t_f = 0$,

* $\lambda(t_0) = \hat{n}(\theta)^{\top} \frac{\partial \Delta r}{\partial z}(z_0)$, with $\theta$ obtained by solving Problem (RED-SHOOT),

* $\nu_1 = -\sigma^{-1}$.

## OPTIMAL SYNTHESIS FOR THE APPROXIMATE PROBLEM

This section is devoted to the development of the optimal synthesis for Problem (RED-SHOOT), Precisely, we aim at achieving a solution for any initial condition of $\Delta r(z_0)$ in the disk of radius $\sigma$, namely $\Delta r(z_0) \mapsto (t_0, \theta)$ for all $\|\Delta r(z_0)\| \le \sigma$.
To provide geometric insight in the developments, we leverage on a case study that can be fully
reproduced with the following inputs (see Remark 6):

$$
y_0 = \begin{pmatrix} 1 \\ 0 \\ 0.05 \\ 0 \\ 0 \\ 30 \text{ deg} \end{pmatrix}, \quad \frac{\Delta v(z_0)}{\|\Delta v(z_0)\|} = \begin{pmatrix} \sin \alpha \sin \beta \\ \cos \beta \\ \cos \alpha \sin \beta \end{pmatrix} \quad \text{with } (\alpha, \beta) = (40, 140) \text{ deg,}
$$

where $y_0$ are a set of equinoctial elements, Keplerian drift is assumed (see Example 1), and $\Delta v(z_0)$ is expressed in the local-normal, local-tangential frame of the primary*. The $\sigma$-to-$\epsilon$ ratio is made explicit in figures' captions. Because the semi-major axis is used to scale lengths, results are valid for orbits of arbitrary size (e.g., LEO and GEO) with the caveat that $\epsilon$ corresponds to different levels of thrust-to-mass ratio (according to its definition in Eq. (1)).

Noting that the final relative position is fully defined by the scalar parameter $\theta$ (Lemma 2, first point), our approach consists in backward integrating Eq. (20) from a one-parameter family of final positions on the endpoint manifold ($\|\Delta r\| = \sigma$). This is done by integrating $2K + 1$ equally-distributed final positions which are then interpolated with a trigonometric polynomial. Specifically, we get (for negative $t_0$):

$$
\frac{\Delta r(t_0, \theta_k)}{\sigma} = \hat{n}(\theta_k) + \frac{\varepsilon}{\sigma} \int_0^{t_0} U(t) \frac{U(t)^{\top} \hat{n}(\theta_k)}{\|U(t)^{\top} \hat{n}(\theta_k)\|} dt, \quad \theta_k = \frac{2\pi k}{2K+1}, \quad k = -K, \dots, K, \quad (21)
$$

so that the solution $\Delta r(z_0) \mapsto (t_0, \theta)$ is implicitly contained in

$$
\Delta r(t_0, \theta) = \sum_{k=-K}^{K} \rho_k(t_0) e^{ik\theta}, \qquad (22)
$$

where $i$ is the imaginary unit, and coefficients of the interpolant, $\{\rho_k(t_0), k = -K, \dots, K\}$, can be
obtained by applying the fast Fourier transform (FFT) algorithm to $\{\Delta r(t_0, \theta_k), k = -K, \dots, K\}$
at any time $t_0$. Equation (22), for a fixed $t_0$, is referred to as *wavefront*, and it represents the locus
of points reaching the safe set in the same time $t_0$.

**Remark 7.** If an explicit expression of the STM, which is necessary to compute $U(t)$, is not avail-
able, $\Phi(t, x)$ can be numerically integrated backward (together with Eq. (21)) by means of the
variational system:

$$
\frac{d\Phi}{dt}(t, y_0) = \frac{\partial F}{\partial y}(y_0) \Phi(t, y_0), \quad \Phi(0, y_0) = I.
$$

This process is shown in Figure (2). The safe set (from which trajectories are integrated backward
in time) is depicted in black and the wavefront is in red. At each fixed $t_0$, the shaded-gray region
denotes the surface of the disk that was not yet covered by the integration. Curves with varying
color are extremal trajectories, i.e., solutions of Eq. (21) for fixed $\theta$ (and the gradient of the color

<footnote>*Unit vectors of this frame are {$\hat{N}, \hat{T}, \hat{H}$}, where $\hat{T}$ is the direction of $v(y_0)$, $\hat{H}$ is toward the angular momentum, and $\hat{N}$ completes the right-hand frame.</footnote>
<figure><img src="image_3.png" /><figcaption>(a) The wavefront (in red) starts contracting inside the safe (b) The wavefront is about to self-intersect: appearence of set (in black). Shaded region is not yet covered. the cut and conjugate loci.</figcaption></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>(c) Trajectories in the two lobes were previously visited (d) The disc is about to be entirely filled. Integration can with a smaller $|t_0|$. Thus they have to be discarded. be stopped.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Backward integration of extremal trajectories from the safe set (i.e., the unit circle). Here, $\sigma/\epsilon = 0.25$.</figcaption></figure>

represents $t_0$). The wavefront starts by contracting inside the disk (Figure 2(a)). This contraction systematically occurs because (use Eq. (20))

$$ \hat{n}(\theta) \cdot \frac{d\Delta r}{dt} = \varepsilon \|U(t)^{\top}\hat{n}(\theta)\| \geq 0. $$

In particular, for all $t < 0$, the time derivative is constrained to the half-plane defined by $\hat{n}(\theta)$ (i.e., exiting point of a trajectory), so that trajectories integrated backward from a given $\hat{n}(\theta)$ cannot evolve in the $\hat{n}(\theta)$ direction itself (a kind of “falling wall” constraint characterizing each trajectory). At some $t_0$, two infinitely-close trajectories achieve the same point (Figure 2(b), trajectories at top-left and bottom-right). The crossing point is the beginning of the cut and conjugate loci, which are illustrated in Figure 3 (dots are conjugate and cut points for the depicted wavefront). This figure clearly shows that cut and conjugate loci originate from the same point. When $t_0$ is further decreased, the wavefront self-intersects (Figure 2(c)). These intersections are referred to as cut points (blue dot in Figure 3). At these points, two trajectories started from two different $\theta$, say $\theta_1$ and $\theta_2$, arrive at the same point with the same $t_0$, namely $\Delta r(t_0, \theta_1) = \Delta r(t_0, \theta_2)$. Trajectories with $\theta \in (\theta_1, \theta_2)$ land on portions of the wavefront lying on a region of the disk that was already visited with a smaller $|t_0|$ (figure 2(c), top and bottom lobes outside the grey-shaded region) and, as
<figure><img src="image_7.png" /><figcaption>Figure 3. Cut and conjugate loci.</figcaption></figure>

such, these solutions are not optimal and have to be discarded. In particular, trajectories landing on the portion of the wavefront between the cut point and a cusp (green points in Figure 3) are locally optimal, but they are not global minimizers, whereas trajectories with endpoint on the portion of the wavefront between two cusps are neither global nor local minima. Cusps are conjugate points, which are points where the local mapping $(t_0, \theta) \mapsto (\Delta r_1, \Delta r_2)$ fails. Hence, conjugate points can be detected by the rank deficiency condition:

$$ \left| \begin{array}{cc} \frac{\partial \Delta r_1}{\partial t_0} & \frac{\partial \Delta r_1}{\partial \theta} \\ \frac{\partial \Delta r_2}{\partial t_0} & \frac{\partial \Delta r_2}{\partial \theta} \end{array} \right| = \frac{\partial \Delta r_1}{\partial t_0} \frac{\partial \Delta r_2}{\partial \theta} - \frac{\partial \Delta r_1}{\partial \theta} \frac{\partial \Delta r_2}{\partial t_0} = 0, $$

which, using Eq. (22) and noting that time derivatives of Fourier coefficients $\rho_k(t_0)$ are the Fourier coefficients of Eq. (20) (i.e., obtained via FFT of the evaluation of the equations of motion at the knots), consists in finding real roots of a trigonometric polynomial of degree $2K$.

The integration is carried out until the disk is filled (Figure 2(d)). This point characterizes the worst-case conjunction scenario for the given $y_0$ and $\Delta v(z_0)$. We note that this point does not necessarily occur at the origin if non-keplerian motion is considered, and it may even be a sub-manifold of dimension 1. Once the disk is filled, we have the optimal solution for any point in the disk. Optimal trajectories originate from the cut locus (two trajectories per cut point), and terminate on the safe set.

Figure 4 shows how trajectories and the cut locus change by varying $\sigma/\epsilon$. We note that the integral in Eq. (21) is independent of this ratio, so that a single integration is sufficient to generate solutions for any $\sigma/\epsilon$. When this ratio is very small (either large thrust-to-mass ratio or very small safe distance), trajectories are nearly rectilinear and quickly leave the disk in the radial direction. This is because the contribution of the drift to the CAM is insignificant, so that trajectories almost minimize the Euclidean distance from the safe set. In this case, the cut locus collapses toward the origin. Conversely, as the ratio is increased, the cut locus expands and tends to align to the vertical direction (recall that the horizontal direction is the projection of $v(t_0)$ into the b-plane).
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 4. Impact of the thrust magnitude to safe distance ratio. The cut locus is depicted in blue.</figcaption></figure>

Indeed trajectories need more time to leave the unsafe region and become nearly horizontal, which suggests that escape is carried out by thrusting in the longitudinal direction and leveraging on the differential mean motion of the primary with respect to its nominal orbit. In this scenario, the CAM is executed in the order of one orbit even for $\sigma/\epsilon = 2.5$. Assuming a safe distance of 10 km, we get a thrust-to-mass ratio of about $10^{-4}$ in LEO and $5 \cdot 10^{-7}$ N/kg in GEO, which is very low, and suggests that alerting few days before the conjunction is safe in most applications. However, different collision scenarios should be tested (i.e., $y_0$ and $\Delta v(z_0)$), and the methodology should be extended to accommodate operational constraints (e.g., forbidden thrusting directions, no thrust during eclipses) and orbital perturbations, which may strongly increase $|t_0|$.

## CONCLUSION

We studied the minimum warning time necessary to carry out a CAM with low thrust, whose purpose is to increase the miss distance of a dangerous conjunction up to a safe threshold. We adopted a local point of view, which consists of using a 7-dimensional state vector obtained by propagating at each time the uncontrolled trajectory of the primary up to the TCA, and the TCA itself. After deducing the equations of motion of these variables, we detailed the necessary conditions for optimality of the problem and studied the optimal synthesis of an approximate solution of the problem in the context of fast encounters.
We believe that the equations that we developed (with their localized point of view) are suitable to be used in various applications ranging from the development of autonomous feedback control laws to the study of multiple conjunctions or slow encounters. Our future research will be devoted to the inclusion of operational constraints on the thrusting windows (e.g., no thrust during eclipses) and direction (e.g., tangential thrust, differential drag) in the methodology, and to the study of the problem with metrics other than the miss distance for the definition of the safeness condition (ellipsoid, probability of collision).

## REFERENCES

[1] NASA, "NASA Spacecraft Conjunction Assessment and Collision Avoidance Best Practices Handbook," 2020.

[2] A. De Vittori, M. F. Palermo, P. D. Lizia, and R. Armellin, "Low-thrust collision avoidance maneuver optimization," *Journal of Guidance, Control, and Dynamics*, Vol. 45, No. 10, 2022, pp. 1815-1829.

[3] R. H. Battin, *An Introduction to the Mathematics and Methods of Astrodynamics, Revised Edition*. American Institute of Aeronautics and Astronautics (AIAA), Jan. 1999, 10.2514/4.861543.

[4] C. Bombardelli, J. Hernando-Ayuso, and R. Armellin, "D Matrix for Keplerian State Error Propagation," *Journal of Guidance, Control, and Dynamics*, Vol. 47, may 2024, p. 901-913, 10.2514/1.g007500.

[5] L. Cesari, *Optimization — theory and applications*, Vol. 17 of *Applications of Mathematics (New York)*. Springer-Verlag, 1983, 10.1007/978-1-4613-8165-5.

[6] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko, *The mathematical theory of optimal processes*. Translated from the Russian by K. N. Trirogoff; edited by L. W. Neustadt, Interscience Publishers John Wiley & Sons, Inc. New York-London, 1962.

[7] R. Vinter, *Optimal control*. Modern Birkhäuser Classics, Birkhäuser Boston, Inc., 2000, 10.1007/978-0-8176-8086-2.

[8] A. A. Agrachev and Y. L. Sachkov, *Control theory from the geometric viewpoint*, Vol. 87 of *Encyclopaedia of Mathematical Sciences*. Berlin: Springer-Verlag, 2004. Control Theory and Optimization, II, 10.1007/978-3-662-06404-7.

[9] H. Schättler and U. Ledzewicz, *Geometric optimal control*, Vol. 38 of *Interdisciplinary Applied Mathematics*. Springer, New York, 2012. Theory, methods and examples, 10.1007/978-1-4614-3834-2.