# Station-keeping under conical constraint on the control force

Alesia Herasimenka, Ariadna Farrés, Lamberto Dell'Elce

► To cite this version:

Alesia Herasimenka, Ariadna Farrés, Lamberto Dell'Elce. Station-keeping under conical constraint on the control force. 9th International Conference on Astrodynamics Tools and Techniques (ICATT), European Space Agency, Jun 2023, Sopot, Poland. hal-04156856

HAL Id: hal-04156856

https://inria.hal.science/hal-04156856

Submitted on 9 Jul 2023

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
STATION-KEEPING UNDER CONICAL CONSTRAINT ON THE CONTROL FORCE

Alesia Herasimenka<sup>(1)</sup>, Ariadna Farrés<sup>(2)</sup>, and Lamberto Dell'Elce<sup>(3)</sup>

<sup>(1)</sup>PhD researcher student, Université Côte d'Azur, CNRS, Inria, LJAD, 06100 Nice, France; ESA contract no 4000134950/21/NL/GLC/my

<sup>(2)</sup>Science Collaborator, University of Maryland Baltimore County, 1000 Hilltop Circle, Baltimore, Maryland, USA

<sup>(3)</sup>Researcher, Inria center at Université Côte d'Azur, 06902 Valbonne, France

ABSTRACT

Many observation satellites are subject to attitude constrains arising from peculiar mission re-
quirements or environmental conditions. These obstructions often constrain the direction of the
thrust vector to remain within a cone. In this study, we investigate the local controllability of
station-keeping maneuvers of satellites with low thrust capabilities or small chemical impulsions
on nominal periodic orbits subject to such constraints. We offer a numerical methodology based
on convex optimization to identify the minimum cone angle guaranteeing local controllability for
a specific orbit. We compare the results with a dynamical approach using Floquet Modes. An
illustrative example inspired by the James Webb Space Telescope is proposed. Specifically, we
consider a satellite is on a Halo orbit around L₂ in the Sun-Earth circular restricted three-body
problem.

# 1 INTRODUCTION

Due to specific mission goals, many satellites are subject to cone constraints on the thrust direction. For example, James Webb Space Telescope (JWST), launched on December 25, 2021 toward a Halo orbit around the Sun-Earth L₂ libration point, has a thermal shield that is designed to prevent the telescope and other instruments from overheating [1]. Therefore, JWST is constrained to always keep its attitude such that the angle between the normal to the shield and the Sun direction is smaller than 53 deg, which results in conical constraints for the propulsion directions. Using chemical propulsion to perform small impulsive corrections maneuvers or a low-thrust satellite with very specific constraints on the control does not always allow to perform any desirable maneuver, as we showed in [2], where the controllability of non-ideal solar sails in orbit about a planet was investigated.

In [2], we considered elliptic Keplerian orbits, and formulated a convex optimization problem aimed at assessing whether some functions of the integrals of motion could not be decreased after one orbital period. Existence of such functions implies that there is a half-space of the neighborhood orbit's coordinates (orbital elements) where motion is locally forbidden [3]. In that paper, we strongly relied on the super-integrability of the Kepler problem. Here, we extend that methodology to infer local controllability of station-keeping satellites for *any* periodic orbit, regardless the dynamical system at hand. Given the projection of the nominal orbit on a surface of section, the methodology aims at verifying if a half space of such projection exists where the motion is forbidden after one orbital period. Variation of parameters is used to achieve a convex optimization problem that investigates the
existence of obstructions to variations of local integrals of motion. Conical constraints are enforced by leveraging on the formalism of positive polynomials postulated by Nesterov [4], so that a finite-dimensional formulation of the convex program is achieved. As a case study, we consider a Halo orbit around the Sun-Earth $L_2$ libration point, and show that there exist a minimum cone thrust angle ($\alpha_{min}$) necessary for the controllability of the orbit. By this we mean, that for spacecraft with a thrust cone angle limitation $\alpha < \alpha_{min}$, we are not able to move everywhere around the neighborhood of an initial position. Even though Halo orbit in the circular restricted three-body problem (CRTBP) is considered in the case study, we emphasize again that the methodology is developed for a generic locally-integrable system.

In this paper, we will also look at the problem from a dynamical system theory approach, as studied in [5]. The dynamics around a periodic orbit in the CRTBP can be described with the Floquet modes, a periodic reference frame that splits the linear dynamics in 3 reference planes: a saddle, a centre and a neutral plane. The instability of the orbit is due to the hyperbolic nature and is contained in the saddle plane. In [5] we also showed how the projection of a thrust direction on the saddle plane is related to the cost of a station-keeping maneuver, and how restrictions on this direction affected the station-keeping strategy. Here we will describe how to use this formulation to derive the cone angle limitations and how it compares to the previous approach.

## 2 EQUATIONS OF MOTION

Consider the equations of motion of a control-affine dynamical system of dimension $n$ with $m$ controls subject to cone constraint on the control, namely

$$ \frac{d x}{d t} = f(x) + B(x) u, \quad x \in M \subseteq \mathbb{R}^n, \quad u \in K_\alpha \subseteq \mathbb{R}^m, \|u\| \le \varepsilon \quad (1) $$

Here, $K_\alpha$ is a cone of revolution characterized by an opening angle $\alpha$, $\varepsilon$ is thrust magnitude, which is assumed to be small, and $f(x)$ denotes a generic drift, e.g., for the CRTBP we have

$$ f(x) = \begin{pmatrix} v_x \\ v_y \\ v_z \\ \left(-\frac{1-\mu}{r_1^3} - \frac{\mu}{r_2^3} + 1\right) r_x - \left(\frac{1}{r_1^3} - \frac{1}{r_2^3}\right) (1-\mu)\mu + 2v_y \\ \left(-\frac{1-\mu}{r_1^3} - \frac{\mu}{r_2^3} + 1\right) r_y - 2v_x \\ \left(-\frac{1-\mu}{r_1^3} - \frac{\mu}{r_2^3}\right) r_z \end{pmatrix}, \quad B(X) = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \quad (2) $$

where $\mu$ is the mass ratio of the system, $x = (r_x, r_y, r_z, v_x, v_y, v_z)$ are position and velocity coordinates in the classical synodic frame, and $r_1$ and $r_2$ are distances between the satellite and the two main bodies:

$$ r_1 = \sqrt{(r_x + \mu)^2 + r_y^2 + r_z^2}, $$

$$ r_2 = \sqrt{(r_x - 1 + \mu)^2 + r_y^2 + r_z^2}. $$

## 3 NECESSARY CONDITION FOR LOCAL CONTROLLABILITY

Given the conical constraint on the thrust vector, $u \in K_\alpha$, we are interested in determining if System (1) is locally controllable. Specifically, given a periodical (uncontrolled) reference orbit $y(t)$ of
<figure><img src="image_2.png" /><figcaption>Figure 1: Forbidden half-space of $\delta x_0$ generated by $p_{\delta x_0}$</figcaption></figure>

period $T$ and a surface of section $S(x)$, and denoting $x_0$ the coordinates of the orbit at the crossing of $S(x)$, namely

$$ \left\{ \begin{array}{l} \frac{d y}{d t} = f(y) \\ y(0) = y(T) = x_0 \\ S(x_0) = 0 \end{array} \right. \qquad (3) $$

we are interested in determining if controls in $K_\alpha$ are capable of moving the crossing point on $S(x)$ in an open neighborhood of $x_0$ after a period $T$, as shown in Fig. 1. To this purpose, we introduce a necessary condition on $\alpha$ for the given orbit in order to have local controllability under the constraint $u \in K_\alpha$.

Denoting by $\Phi(t, x_0)$ the state transition matrix of the system, and by $\delta x_0 \in T_{x_0} S$ a perturbation of the initial state $x_0$, uncontrolled linearized motion in proximity of the periodic orbit is governed by

$$ \delta x(t) = \Phi(t, x_0) \delta x_0. \qquad (4) $$

Linearization of Eq. (1) gives:

$$ \frac{d \delta x}{d t} = \left. \frac{\partial f}{\partial x} \right|_{y} \delta x + B(y) u. \qquad (5) $$

Recalling that $\frac{d \Phi}{d t} = \frac{\partial f}{\partial x} \Phi$, differentiation of Eq. (4) and substitution in Eq. (5) yields the classical variation of parameters

$$ \frac{d \delta x_0}{d t} = \Phi^{-1}(t, x_0) B(y(t)) u, \quad \delta x_0 \in T_{x_0} S, \quad u \in K_\alpha. \qquad (6) $$

The necessary condition for local controllability of the satellite is written in terms of possible displacements of the system on the Poincaré map, i.e. by verifying if the system can be moved everywhere
in the tangent space $T_{x_0}S$ after one orbital period. For mathematical proof of the necessary condition please refer to [3]. Negation of this condition implies the existence of a not accessible half-space in the neighborhood of $x_0$, as shown in Fig. 1. Since the interior thrust directions of $K_\alpha$ can be approximated by combinations of vectors on the boundary of the cone, $\partial K_\alpha$, we propose to solve the following problem in order to verify the necessary condition:

$$ \text{if } \exists p_{\delta x_0} \in T_{x_0}^*S, p_{\delta x_0} \neq 0 \text{ such that} \\ \left\langle p_{\delta x_0}, \frac{\mathrm{d}\delta x_0}{\mathrm{d}t} \right\rangle \geq 0, \quad \forall u \in \partial K_\alpha, \|u\| = 1, t \in [0, T) \quad (7) $$

then System (1) is not locally controllable in one orbit.

If $p_{\delta x_0}$ solution of Problem (7) exists, then the linear functional

$$ V(t, u) = \langle p_{\delta x_0}, \Phi^{-1}(t, x_0) B(y(t)) u \rangle $$

cannot be decreased for any $u \in K_\alpha$ and $t \in [0, T)$, hence motion is forbidden in the half-space with normal $p_{\delta x_0}$, and the satellite cannot move in any direction pointing inside this half-space after one orbital period. Absence of forbidden directions for control of satellites is crucial for station-keeping.

# 4 CONVEX OPTIMIZATION PROBLEM TO VERIFY THE NECESSARY CONDITION

A practical check of the necessary condition is carried out by solving

$$ \max_{J, \|p_{\delta x_0}\| \le 1} J \quad \text{s.t.} \\ \langle p_{\delta x_0}, \Phi^{-1}(t, x_0) B(y(t)) u \rangle \ge J, \quad \forall u \in \partial K_\alpha, \|u\| = 1, t \in [0, T]. \quad (8) $$

Problem (8) is convex and semi-infinite, because inequality constraints need to be enforced for all $u$ on the surface of the cone and for all time between 0 and the period $T$. Evaluating inequalities in the interior of the cone is not necessary because dynamics is affine in $u$. If $J^*$, solution of Problem (8), is strictly positive, then the necessary condition is not satisfied and the system is not locally controllable for the given $\alpha$ and $x_0$. The constraint $\|p_{\delta x_0}\| \le 1$ is preferred to the equality condition $\|p_{\delta x_0}\| = 1$ to preserve the convexity properties of Problem (8).

For mission design purposes, it is interesting to know which is the minimum $\alpha$ angle of the thrust cone satisfying the necessary condition. This angle can be identified by solving

$$ \min_{\alpha} \alpha \quad \text{s.t.} \\ J^*(\alpha) = 0 \quad (9) $$

where $J^*(\alpha)$ denotes solution of Problem (8) for a given $\alpha$. Problem (9) can be efficiently solved by means of a simple bisection method.

Numerical solution of Problem (8) is achieved by:

1. Parametrizing $K_\alpha$ by means of an angle $\delta$, as shown in Fig. 2, to avoid discretization of the cone by using, for example, a polyhedral cone with a finite number of generators. Thus, vectors of $u$ on the surface of the cone can be expressed as:

$$ u = \begin{bmatrix} \cos \alpha \\ \cos \delta \sin \alpha \\ \sin \delta \sin \alpha \end{bmatrix} \quad (10) $$

with $\alpha \in [-\frac{\pi}{2}, \frac{\pi}{2}]$ and $\delta \in [0, 2\pi]$;
<figure><img src="image_3.png" /><figcaption>Figure 2: Parametrization of the control vector.</figcaption></figure>

2. Given that $u$ is trigonometric in $\delta$, using Fourier transform for Eq. (6):

$$ \Phi^{-1}(t, x_0) B(y(t)) u = \sum_{l=-1}^{1} \sum_{k=-d}^{d} C^{(k,l)} e^{ikt} e^{il\delta} $$

where $C^{(k,l)}$ is the $kl$-th coefficient of the Fourier transform of $\Phi^{-1}(t, x_0) B(y(t)) u$ and $d$ the degree of truncation of the series in $t$. Note that $u$ is already an exact trigonometric polynomial of degree 1 in $\delta$. Thus, the inequality from Eq. (8) becomes:

$$ \langle p_{\delta x_0}, \Phi^{-1}(t, x_0) B(y(t)) u \rangle \geq J \iff p_{\delta x_0}^T \left( \sum_{l=-1}^{1} \sum_{k=-d}^{d} C^{(k,l)} e^{ikt} e^{il\delta} \right) - J \geq 0 \quad (11) $$

In the example of a Halo orbit given in this paper, we decide to truncate the Fourier series at $d = 30$, as the convergence of the coefficients is enough to find the minimum cone angle, as shown in Fig. 3;

3. Using the formalism of positive polynomials [4, 6] to enforce positivity constraints.

Consider the basis of bivariate trigonometric polynomials of degree $d$ in $t$ and 1 in $\delta$: $\mathcal{P}(t, \delta) = [1, e^{i\delta}]^T \otimes [1, e^{it}, e^{2it}, \dots, e^{dit}]^T = [1, e^{it}, e^{2it}, \dots, e^{dit}, e^{i\delta}, e^{it}e^{i\delta}, e^{2it}e^{i\delta}, \dots, e^{dit}e^{i\delta}]^T$ and $C$ vector of coordinates of the polynomial in the basis. Its corresponding squared functional system is $S^2(t, \delta) = \mathcal{P}(t, \delta)\mathcal{P}^H(t, \delta)$, where $\mathcal{P}^H(t, \delta)$ denotes conjugate transpose of $\mathcal{P}(t, \delta)$. Let $N$ be the dimension of $\mathcal{P}(t, \delta)$ ($N = 2 \times (d+1)$ in our application) and $\Lambda_H: \mathbb{C}^N \to \mathbb{C}^{N \times N}$ a linear operator
<figure><img src="image_4.png" /><figcaption>Figure 3: Convergence of Fourier coefficients</figcaption></figure>

mapping coefficients of polynomials in $\mathcal{P}(t, \delta)$ to the squared base, so that application of $\Lambda_H$ on $\mathcal{P}(t, \delta)$ yields

$$ \Lambda_H(\mathcal{P}(t, \delta)) = \mathcal{P}(t, \delta)\mathcal{P}^H(t, \delta) \quad (12) $$

and define its adjoint operator $\Lambda_H^* : \mathbb{C}^{N \times N} \to \mathbb{C}^N$ as

$$ \langle Y, \Lambda_H(C) \rangle_H = \langle \Lambda_H^*(Y), C \rangle_H, \quad Y \in \mathbb{C}^{N \times N}, \quad C \in \mathbb{C}^N. \quad (13) $$

Theory of squared functional systems postulated by Nesterov [4] proves that trigonometric polynomial is non-negative if and only if a Hermitian positive semidefinite matrix $Y$ exists such that $C = \Lambda_H^*(Y)$. Dumitrescu extends this theory for multivariate trigonometric polynomials in [6] and shows that all nonnegative bivariate trigonometric polynomials can be written as sum-of-squares. This equivalence is false for three or more variables.

Thus, $\langle \mathcal{P}(t, \delta), C \rangle_H$ is non-negative for all $t \in [0, T)$ and for all $u \in K_\alpha$ if and only if a Hermitian positive semidefinite matrix $Y$ exists such that $C = \Lambda_H^*(Y)$, namely

$$ \langle \mathcal{P}(t, \delta), C \rangle_H \ge 0, \quad t \in [0, T), \quad u \in K_\alpha \iff \exists Y \succeq 0 : C = \Lambda_H^*(Y). \quad (14) $$

In fact, it holds in this case that

$$ \begin{aligned} \langle \mathcal{P}(f, \delta), C \rangle_H &= \langle \mathcal{P}(f, \delta), \Lambda_H^*(Y) \rangle_H = \langle \Lambda_H(\mathcal{P}(f, \delta)), Y \rangle_H, \\ &= \langle \mathcal{P}(f, \delta)\mathcal{P}^H(f, \delta), Y \rangle_H = \mathcal{P}^H(f, \delta)Y\mathcal{P}(f, \delta) \ge 0. \end{aligned} \quad (15) $$

For trigonometric polynomials $\Lambda^*$ is given by

$$ \Lambda_H^*(Y) = \begin{bmatrix} \operatorname{tr}(\langle Y, T_{00} \rangle) \\ \vdots \\ \operatorname{tr}(\langle Y, T_{kl} \rangle) \\ \vdots \\ \operatorname{tr}(\langle Y, T_{21} \rangle) \end{bmatrix} \quad k = 0, 1, 2, \ l = 0, 1. \quad (16) $$
where $T_j, j = 0, 1, 2$ are the elementary Toeplitz matrices with ones on the j-th diagonal and zeros elsewhere and $T_{kl}$ are obtained from a Kronecker product of such matrices, e.g.,

$$ T_0 = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}, \quad T_1 = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{pmatrix} \\ T_{10} = T_0 \otimes T_1 = \begin{pmatrix} 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} \tag{17} $$

Finally, the inequality in Eq. (8) is rewritten as an linear matrix inequalities (LMI):

$$ \langle p_{\delta x_0}, \Phi^{-1}(t, x_0) Bu \rangle - J \geq 0, \ t \in [0, T), \ u \in \partial K_{\alpha} \\ \iff \exists Y \succeq 0 \text{ such that } C p_{\delta x_0} - e_1 J = \Lambda_H^*(Y) \tag{18} $$

where $Y \in \mathbb{C}^{N \times N}$ is a Hermitian matrix to be determined, with $N = 2 \times (d+1) = 62$, and $e_1$ is a vector of dimension $N$ with 1 in the first position and zeros elsewhere. Hence, the finite-dimensional counterpart of Problem (8) is

$$ \begin{array}{ll} \displaystyle \min_{J, \|p_{\delta x_0}\| \le 1, Y \in \mathbb{C}^{62 \times 62}} & J \quad \text{s.t.:} \\ & Y \succeq 0 \\ & \Lambda_H^*(Y) = C p_{\delta x_0} - e_1 J \end{array} \tag{19} $$

Solution of Problem (9) is carried out by means of a simple bisection algorithm, which does not require the evaluation of derivatives of the non-smooth function $J^*(\alpha)$ (we note that Problem (8) has trivial solution $J = 0$, $p_{\delta x_0} = 0$ for $\alpha > \alpha_{min}$). The CVX software [7, 8] is used to solve the convex Problem (19). Fourier coefficients of $\Phi^{-1}(t, x_0) B(y(t)) u$ are evaluated by means of the fast Fourier transform (FFT) algorithm. The only relaxation of Problem (19) with respect to Problem (8) is truncation of the Fourier series. Remarkably, no discretization was done to approximate $u$ on the surface of a cone.

# 5 DYNAMICS AROUND HALO ORBITS

Understanding the local behavior of the flow associated to a dynamical system near a periodic orbit gives insight on station-keeping strategies, as well as the effects of different error sources. The study of the local behaviour is generally done throughout the first order variational equations. To fix notation, let us denote by $\varphi$ the flow associated to the equations of motion (Eq. 2). Hence, $\varphi_{\tau}(x_0)$ is the image at a time $t = \tau$ of a point $x_0 \in \mathbb{R}^6$ at $t = 0$, and $\Phi(\tau, x_0) = D\varphi_{\tau}(x_0)$ is the first order variational of $\varphi_{\tau}(x_0)$ with respect to the initial condition $(x_0)$. For $h \in \mathbb{R}^6$, we have

$$ \varphi_{\tau}(x_0 + h) = \varphi_{\tau}(x_0) + D\varphi_{\tau}(x_0) \cdot h + O(|h|^2). $$

Provided $h$ small, $\varphi_{\tau}(x_0) + \Phi(\tau, x_0) \cdot h$ gives a good approximation of $\varphi_{\tau}(x_0 + h)$. Finally, the linear dynamics around a periodic orbit can be described through the eigenvalues and eigenvectors of the monodromy matrix of the orbit $\Phi(T)$, where $T$ is the period of the orbit.
It is well known, that for most Halo orbits in the CRTBP, the eigenvalues ($\lambda_1,...,\lambda_6$) of the monodromy matrix, $\Phi(T)$, satisfy: $\lambda_1 > 1$, $\lambda_2 < 1$ are real, $\lambda_3 = \bar{\lambda}_4$ are complex with modulus 1 and $\lambda_5 = \lambda_6 = 1$. Where the geometrical meaning of these three pairs of eigenvalues and their associated eigenvectors is [9]:

* The first pair ($\lambda_1, \lambda_2 \in \mathbb{R}$, verify $\lambda_1 \cdot \lambda_2 = 1$, and are related to the hyperbolic character of the orbit. The value $\lambda_1$ is the largest in absolute value, and is related to the eigenvalue $e_1(0)$, which gives the most expanding direction. After one period, a given distance to the nominal orbit in this direction is amplified by a factor of $\lambda_1$. Moreover, at each point of the orbit, the vector $e_1(\tau) = D\varphi_\tau e_1(0)$ together with the vector tangent to the orbit, span a plane that is tangent to the local unstable manifold. In the same way $\lambda_2$ and its related eigenvector $e_2(0)$ are related to the local stable manifold with $e_2(\tau) = D\varphi_\tau e_2(0)$.

* The second pair ($\lambda_3, \lambda_4 \in \mathbb{C}$, verify $\lambda_3 = \bar{\lambda}_4$ and of modulus 1. Together with the other two eigenvalues ($\lambda_5 = \lambda_6$) describe the central motion around the periodic orbit. The monodromy matrix restricted to the plane spanned by the real $e_3(0)$ and imaginary $e_4(0)$ parts of the eigenvectors associated to $\lambda_3, \lambda_4$ along the orbits is a rotation of of angle $\Gamma = \arctan\left(\frac{Im(\lambda_3)}{Re(\lambda_3)}\right)$.

* The third couple ($\lambda_5 = \lambda_6 = 1$) is associated to the neutral directions (i.e. non-unstable modes). There is only one eigenvector of $\Phi(T)$ with eigenvalue 1, this vector is the tangent vector to the orbit and we call it $e_5(0)$. The other eigenvalue is associated to variations of the period, or any other variable which parameterised the family of periodic orbits.

Finally, the functions $e_i(\tau) = D\varphi_\tau e_i(0)$, $i = 1, ..., 6$ can be used to describe the linear dynamics of the phase space around a periodic orbit, and from them we can define the Floquet Modes around the Halo orbit:

$$ \bar{e}_{1,2}(\tau) = e_{1,2}(\tau) \exp \left( -\tau \frac{\ln \lambda_{1,2}}{T} \right), \qquad (20) $$

$$ \bar{e}_3(\tau) = \cos \left( -t \frac{\Gamma}{T} \right) e_3(\tau) - \sin \left( -\tau \frac{\Gamma}{T} \right) e_4(\tau), \qquad (21) $$

$$ \bar{e}_4(\tau) = \sin \left( -t \frac{\Gamma}{T} \right) e_3(\tau) + \cos \left( -\tau \frac{\Gamma}{T} \right) e_4(\tau), \qquad (22) $$

$$ \bar{e}_5(\tau) = e_5(\tau), \qquad (23) $$

$$ \bar{e}_6(\tau) = e_6(\tau) + \varepsilon(\tau)e_5(\tau), \qquad (24) $$

where T is the period of the orbit (for further details see [9]).

The Floquet Modes are a periodic reference frame that helps us describe the natural dynamics in the vicinity of a periodic orbit. Using the Floquet Modes reference frame, the dynamics in the vicinity of the orbit is simple: on the planes generated by $\{\bar{e}_1(t), \bar{e}_2(t)\}$ the motion is a saddle, having the trajectory escape, with an exponential rate, along the unstable direction; on the planes generated by $\{\bar{e}_3(t), \bar{e}_4(t)\}$ the dynamics consists of a rotation around the periodic orbit; and on the planes generated by $\{\bar{e}_5(t), \bar{e}_6(t)\}$ the dynamics is neutral. Figure 4 shows the evolution of a schematic representation of the linear dynamics around a Halo orbit using this reference frame, where the origin of coordinates (in each projection) corresponds to the periodic orbit.

## 5.1 Station-keeping strategies

The goal of any station-keeping strategy is to keep the spacecraft close to a nominal orbit. Given the instability around a Halo orbit, in order to remain close to it, small correction maneuvers are required
<figure><img src="image_5.png" /><figcaption>Figure 4: Schematic representation of the dynamics around a halo orbit using the Floquet modes $\bar{e}_i(t)$.</figcaption></figure>

from time to time, where the size of these maneuvers depends on their frequency and thrust direction [1, 5]. In the literature we find different approaches to determine the optimal delta-v maneuvers. A classical approach used in many missions like the JWST, is the velocity constraint approach, where the goal is to find a delta-v maneuver that ensures that after 1.5 orbital revolutions, the $v_x$ component of the trajectory when crossing the $y=0$ plane is zero, which can be solved with a classical targeting method. The underlining idea behind this strategy is to keep orbiting about the libration point, since a Halo orbit satisfies $v_x = 0$ at the plane $y=0$.

We use the Floquet Modes reference frame to visualize the behaviour of the trajectory each time a station-keeping maneuver is applied. Note that a delta-v maneuver is an instantaneous change in velocity, and is seen as a jump in the phase space. The first row of Figure 5 shows (in dimensionless units (DU)) the evolution of a trajectory close to a Halo orbit when no station-keeping maneuvers are applied. While the second row of Figure 5 is the trajectory when station-keeping maneuvers are applied. On the left-hand side of the plot we have the evolution of the trajectory on the XYZ plane. The other three plots (from left to right) are the projection of the trajectory on the saddle, centre and neutral plane defined by the Floquet modes. As we can see on the top plot, the trajectory escapes along the unstable direction on the saddle plane projection, while rotation around the origin in the centre plane. Note that once the trajectory is far from the nominal orbit the linear dynamics no longer holds, showing an unexpected behaviour. On the bottom plots, the station-keeping maneuvers are jumps in the trajectory. Notice how on the saddle projection, these maneuvers cancel the instability by bringing the trajectory towards the stable direction ($e_2$). On the centre plane we see a sequence of rotations around the centre, and each time a maneuver is applied the distance to the centre varies.

In [5] we compared this classical approach (velocity constraint) with the Floquet Mode approach, showing that their behaviour in the Floquet modes reference frame was the same: each time a station-keeping maneuver was performed, the trajectory will come close to the stable manifold of the nominal orbit, canceling the unstable component of the trajectory, as we can see in Figure 5. We also saw that the magnitude of the station-keeping maneuver is related to the jump on the saddle plane, from the time of the maneuver towards the stable direction, along a given delta-v direction.

To illustrate the relation between the size of a station-keeping magnitude and the thrust direction we have performed the following simulations. Take an initial condition on the halo orbit displacing $10^{-6}$ DU (≈ 150 km) along the unstable direction (i.e., its position in the Floquet mode basis is $s = (\pm10^{-6}, 0, 0, 0, 0, 0)^T$). Then, for a set of different thrust directions, determine the required station-keeping maneuver along that thrust direction. To simplify the analysis, we have considered only thrust directions in the $xy$-plane $u = (0, 0, 0, \cos\alpha, \sin\alpha, 0)^T$, with $\alpha \in [-90, 90]^{\circ}$ (Note that $\alpha$ corresponds to the cone angle).

Figure 6 shows the relation between the magnitude of delta-v maneuver and the jump on the saddle plane. The two left subplots show the projection on the saddle plane of the delta-v maneuver, while the
<figure><img src="image_6.png" /><figcaption>(a) Evolution of the Floquet components for a trajectory close to halo orbit with no station-keeping maneuvers.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Evolution of the Floquet components for a trajectory close to a halo orbit with station-keeping maneuvers.</figcaption></figure>

Figure 5: Comparison between trajectories with and without station-keeping maneuvers. Black arrows show the sense of orbital motion, and cyan arrows show the jumps associated to maneuvers.

right subplots show the size of the maneuver as a function of the cone angle $α$. Note how on the left plots, depending on the initial location on the saddle plane ($c_1 = +10^{-6}$ or $c_1 = -10^{-6}$), there is only a set of directions that allows one to control the spacecraft, and that in all cases the maneuver cancels the unstable mode (i.e. produces a jump in velocities such that $c_1 = 0$). The two right plots use the same color-code as the corresponding left plot for the size of the delta-v maneuver. Where the vertical red and green lines correspond to the $α$ values of the direction given by the position components of the unstable and stable directions, respectively. Notice that the maneuver directions with smaller magnitudes correspond to those directions whose projection on the saddle plane is close to $(±1, 0)$ and are close to the stable direction (green dotted line). While the largest maneuver magnitudes are close to $(0, ±1)$ on the saddle plane and are close to the unstable direction (red dotted line). Finally, describing the station-keeping strategy throughout its projection on the saddle plane can help us determine for which cone angles the system is not controllable. As we can see in Figure 6, the set of delta-v directions that can bring the trajectory towards the stable direction will depend on where the spacecraft is at the time of the maneuver on the saddle plane projection. However, regardless of where we are on the saddle plane, maneuvers whose projection on the saddle plane are parallel to $(0, ±1)$ will not be able to bring the trajectory towards the stable manifold.

# 6 CASE STUDY

Let us consider a periodical Halo orbit situated around Sun-Earth L2 point, as shown in Fig. 7. It is the same point where James Webb Space Telescope was sent. We suppose that a satellite has to perform station-keeping around this orbit. The given satellite can produce either small impulsions using chemical propulsion or low-thrust engines, and has a conical constraint on the directions of the thrust. Our goal is to determine what is the maximum conical constraint that can be imposed on the
<figure><img src="image_8.png" /><figcaption>Figure 6: Relation between the maneuver magnitude and the jump on the saddle plane.</figcaption></figure>

Left: Projection of delta-v maneuvers on the saddle plane. Right: delta-v magnitude as a function of the $\alpha = \arctan(\Delta v_y / \Delta v_x)$.

propulsion, i.e. what is the minimum cone angle for thrust directions that allows local controllability after one orbital period.

To find out the minimum requirement, we apply the two different approaches described in section 4 and 5 on the given periodical Halo orbit. The initial data of the orbit is $x_0 = (1.0083, 5.15 \times 10^{-19}, 0.0010, 1.3714 \times 10^{-16}, 0.0102, -4.1015 \times 10^{-17})$ in AU according to the Sun-centered reference frame.

<figure><img src="image_9.png" /><figcaption>Figure 7: Halo orbit used for the simulation. Size of the Sun is schematical.</figcaption></figure>

## 6.1 Minimum Cone angle using convex optimization

The results given by Fig. 8 show that the minimum thrust cone angle $\alpha_{min} = 43$ deg exists, and is a necessary requirement for local controllability of a station-keeping satellite using low-thrust or small chemical impulsions. The results mean that a satellite with thrust directions limited by a cone of
<figure><img src="image_10.png" /><figcaption>Figure 8: Solution of Problems (8) (black curve) and (9) (red dot).</figcaption></figure>

less than 43 degrees is not capable of moving anywhere around the neighborhood of an intial position after one orbital period. Moreover, it indicates that there exists a half-space of the initial configuration neighborhood which includes all forbidden directions. For example, the satellite might not be capable of raising its velocity in *y*-direction, or decreasing its *z*-position after one orbital period, therefore it is not locally controllable in one orbital period. It might be possible that global controllability still holds, in other words that the satellite has to move away from the initial orbit to perform the necessary maneuver and then to come back. Nevertheless, it would require an important amount of propellant or it is probably not feasible by the low-thrust engines.

## 6.2 Minimum Cone angle using the Floque Modes

Let us consider that the thrust vector *u* parameterized by the cone angle *α* and the clock angle *δ*, having $u = (0, 0, 0, \cos\theta, \sin\theta \cos\psi, \sin\theta \sin\psi)$. Where we constrain the thrust vector to $u_x > 0$, i.e. $\alpha \in [0, 90]^{\circ}$ and $\delta \in [-180, 180]^{\circ}$. Figure 9 shows the projection of *u* in the RLP reference frame (top) and on the saddle plane (bottom). The color-code in the plots is used to identify vectors with the same cone angle (left) and the same clock angle (right). Notice that most of the thrust directions in the saddle plane have $c_1 > 0$, having a limited set of thrust directions that ensure jumps in the saddle plane in the direction toward the Sun. As we can see, if the cone angle is close to 50 deg (light green arrows), the projection on the saddle plane is close to (0, -1) which would not allow to cancel the unstable component of the trajectory.

Notice that the saddle projection of the thrust vector *u* will vary as the spacecraft moves along the periodic orbit, as the stable and unstable directions vary over one orbital period. In order to determine the minimum cone angle that enables controllability, one must compute for different points along the orbit, the delta-v vector that is projected into (0, ±1) on saddle plane, and its corresponding cone angle. Figure 10 show the variation over one orbital period of this minimum cone angle, where we can see that the minimum value corresponds to 44 deg.
<figure><img src="image_11.png" /><figcaption>Figure 9: Projections of the thrust vectors $u = (0, 0, 0, \cos \alpha, \sin \alpha \cos \delta, \sin \alpha \sin \delta)$ as a function of pitch angle $\alpha$ and clock angle $\delta$. Top: projection of the velocity components on the RLP frame. Bottom: projection of $u$ on the saddle plane $\bar{e}_1, \bar{e}_2$. Left: color code groups same pitch angles. Right color code groups same clock angles.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Variation of the cone angle of the thrust vectors $u = (0, 0, 0, dvx, dvy, dvz)$ whose projection on the saddle plane is (0,1).</figcaption></figure>

# 7 CONCLUSION

In this paper we propose a methodology to find the minimum requirement for station-keeping of the satellites with cone-constrained thrust. Our analysis is inspired by the James Webb Space Telescope, which has to maintain the imposed attitude towards the Sun because of the solar shield protecting
its instruments. We formulate a convex optimization problem giving a solution in terms of a mini-
mum cone angle of the thrust directions allowing local controllability. In other words, the minimum
condition is a necessary condition for the satellite to be capable of moving anywhere to maintain its
position on the orbit. The proposed methodology consists in finding a forbidden half-space in the
neighborhood of the initial configuration of the satellite on the Poincaré map where it cannot move
after one orbital period. The optimization problem is solved using convex programming and theory
of positive bivariate trigonometric polynomials.

We also compare the methodology with Floque modes approach, where the projection of a thrust
direction on the saddle plane allows to evaluate station-keeping maneuvers. The obtained results are
very close, but a small gap exists. More analysis should be done to better understand this difference.
The minimum requirement that we propose can be used for a design of space missions around any
periodic orbit for satellites that have specific constraints on the thrust directions. It can be applied to
a low-thrust satellite or even those with chemical propulsion under condition of using small impulses,
so that the linearization of the dynamics holds.

REFERENCES

[1] J. Petersen, “L2 Station Keeping Maneuver Strategy For The James Webb Space Telescope,” *AIAA/AAS Astrodynamics Specialist Conference*, American Institute of Aeronautics and Astronautics, 2019.

[2] A. Herasimenka, L. Dell'Elce, J.-B. Caillau, and J.-B. Pomet, “Controllability Properties of Solar Sails,” *Journal of Guidance, Control, and Dynamics*, 2022. in press.

[3] J.-B. Caillau, L. Dell'Elce, A. Herasimenka, and J.-B. Pomet, “On the Controllability of Nonlinear Systems with a Periodic Drift,” 2022. HAL preprint no. 03779482.

[4] Y. Nesterov, “Squared Functional Systems and Optimization Problems,” *High Performance Optimization* (P. M. Pardalos, D. Hearn, H. Frenk, K. Roos, T. Terlaky, and S. Zhang, eds.), Vol. 33, pp. 405–440, Boston, MA: Springer US, 2000. Series Title: Applied Optimization, 10.1007/978-1-4757-3216-0_17.

[5] A. Farrés, C. Gao, J. J. Masdemont, G. Gómez, D. C. Folta, and C. Webster, “Geometrical Analysis of Station-Keeping Strategies About Libration Point Orbits,” *Journal of Guidance, Control, and Dynamics*, Vol. 45, No. 6, 2022, pp. 1108–1125, 10.2514/1.G006014.

[6] B. Dumitrescu, *Positive Trigonometric Polynomials and Signal Processing Applications*. Springer Netherlands, 2007, 10.1007/978-1-4020-5125-8.

[7] M. Grant and S. Boyd, “CVX: Matlab Software for Disciplined Convex Programming, version 2.1,” http://cvxr.com/cvx, Mar. 2014.

[8] M. Grant and S. Boyd, “Graph implementations for nonsmooth convex programs,” *Recent Advances in Learning and Control* (V. Blondel, S. Boyd, and H. Kimura, eds.), Vol. 371 of *Lecture Notes in Control and Information Sciences*, pp. 95–110, Springer-Verlag, 2008. http://stanford.edu/~boyd/graph_dcp.html.

[9] G. Gómez, A. Jorba, C. Simó, and J. Masdemont, *Dynamics and Mission Design Near Libration Points*. WORLD SCIENTIFIC, 2001, 10.1142/4337.