Andrea C Morelli

Carmine Giordano
email: carmine.giordano@polimi.it.aiaa

Riccardo Bonalli
email: riccardo.bonalli@cnrs.fr

Francesco Topputo
email: francesco.topputo@polimi.it.aiaa

Ph D Candidate

Characterization of Singular Arcs in Spacecraft Trajectory Optimization

I. Introduction

The recent advancements in the field of low-thrust engines for interplanetary spacecraft transfers have opened completely new mission scenarios. The high efficiency of such thrusters and their continuous thrust enabled cost-effective space missions with flexible launch and arrival dates [START_REF] Topputo | Envelop of reachable asteroids by M-ARGO CubeSat[END_REF]. Consequently, an increasing number of missions are expected to adopt this technology in the near future.

In space missions, trajectories are usually designed to minimize some key objective, such as the fuel mass consumption [START_REF] Morelli | Robust Low-Thrust Trajectory Optimization Using Convex Programming and a Homotopic Approach[END_REF]. Due to their continuous thrust, low-thrust engines require an Optimal Control Problem (OPC) to be solved to find such trajectories [START_REF] Bryson | Applied optimal control[END_REF]. State-of-the-art methods that are used to solve the problem divide into direct and indirect ones.

The former discretize the continuous-time problem and solve the resulting (non)linear program, whereas the latter solve a two-point boundary value problem formulated using the calculus of variation [START_REF] Betts | Survey of Numerical Methods for Trajectory Optimization[END_REF].

By employing the Pontryagin's Maximum Principle (PMP) [START_REF] Kirk | Optimal Control Theory: An Introduction[END_REF] it can be proved that the optimal control of the low-thrust trajectory optimization (LTO) problem follows a bang-bang structure, i.e., the throttle factor should be either at the maximum or the minimum value [START_REF] Topputo | Survey of Direct Transcription for Low-Thrust Space Trajectory Optimization with Applications[END_REF]. However, there could also be cases where the first-order optimality conditions cannot provide any information on the structure of the optimal control. In that case, the control is usually referred to as singular [START_REF] Kirk | Optimal Control Theory: An Introduction[END_REF]. A thorough analysis of singular control in the powered descent and landing (PDG) problem was recently performed [START_REF] Leparoux | Structure of optimal control for planetary landing with control and state constraints[END_REF]. In that case, due to the constant free dynamics, there could only be one singular arc in a given trajectory, and only for specific initial conditions. However, the results are not directly applicable to the LTO problem given the highly nonlinear free dynamics. In practical applications, if unperturbed two-body dynamics with thrust acceleration are considered, singular thrust arcs are rare and therefore it is not uncommon to solve the LTO problem as if they did not exist [START_REF] Russell | Primer vector theory applied to global low-thrust trade studies[END_REF][START_REF] Taheri | Enhanced smoothing technique for indirect optimization of minimumfuel low-thrust trajectories[END_REF][START_REF] Bertrand | New smoothing techniques for solving bang-bang optimal control problems numerical results and statistical interpretation[END_REF][START_REF] Tang | Fuel-Optimal Low-Thrust Trajectory Optimization Using Indirect Method and Successive Convex Programming[END_REF][START_REF] Nurre | Duty-cycle-aware low-thrust trajectory optimization using embedded homotopy[END_REF], though this is only supported by empirical results rather than theoretical evidence. Previous work has shown that singular arcs can in fact theoretically happen in a two-body dynamical environment and trajectories with several singular arcs can actually be designed [START_REF] Azimov | Extremal analytical solutions for intermediate-thrust arcs in a newtonian field[END_REF]. A further example are the Lawden's spirals for intermediate thrust arcs [START_REF] Lawden | Optimal Trajectories for Space Navigation[END_REF].

Others have formulated necessary conditions for the optimality of singular arcs in the case of multiple gravitational bodies [START_REF] Park | Necessary conditions for the optimality of singular arcs of spacecraft trajectories subject to multiple gravitational bodies[END_REF]. So far, researchers have focused on investigating the optimality of singular arcs rather than quantifying, theoretically and numerically, their occurrence. Moreover, in many cases, necessary conditions and expressions of the singular controls were expressed as a function of both states and costates [START_REF] Park | Necessary conditions for the optimality of singular arcs of spacecraft trajectories subject to multiple gravitational bodies[END_REF], making it impossible to have a physical grasp on the problem. This work proposes for the first time a full characterization of singular arcs in spacecraft LTO in a two-body dynamical environment, but the logic applies to more complex dynamics as well. In reality, more complex dynamics are usually considered when designing a spacecraft trajectory (e.g., 𝑛-body problem). The bang-bang structure of the optimal control does not depend on the specific dynamics, but the same does not apply to the frequency of singular arcs, which could appear more often. It is therefore relevant to investigate the behaviour of singular arcs to understand whether mission analysts can safely assume that a bang-bang control accurately captures the solution of a LTO problem. The contribution of this article is threefold. First, necessary conditions for having singular arcs that solely depend on three physical variables are provided. Moreover, the singular control is defined through an expression that depends on a limited set of physical variables. Finally, based on our theoretical results, we provide insightful information about the reasons why in practical applications singular arcs rarely appear. In particular, we show that the necessary conditions have a limited number of solutions, given the orbit eccentricity and true anomaly. In addition, we show that for typical inner Solar System missions, even in the case the necessary conditions are satisfied, the expression of the singular throttle factor is singular in the minority of the cases.

The remainder of the paper is organized as follows. Section II formulates the considered LTO problem. Section III describes the approach that has been used to characterize the thrust arcs, as well as the theoretical findings. Section IV presents the numerical simulations. Finally, V concludes the work.

II. Problem Statement

The two-body dynamics of a spacecraft around a primary body and equipped with a low-thrust engine can be expressed in Cartesian coordinates as [START_REF] Jiang | Practical Techniques for Low-Thrust Trajectory Optimization with Homotopic Approach[END_REF] out-of-plane angles of the thrust vector. 𝑇 max is the constant maximum thrust of the engine, 𝐼 sp is the specific impulse, 𝑔 0 is the Earth gravity acceleration at sea level, and 𝜇 is the gravitational parameter of the primary body. Equation [START_REF] Topputo | Envelop of reachable asteroids by M-ARGO CubeSat[END_REF] can be written in vectorial form as

                                                     𝑥 = 𝑣 𝑥 𝑦 = 𝑣 𝑦 𝑧 = 𝑣 𝑧 𝑣 𝑥 = -𝜇 𝑥 ∥r∥ 3 + 𝑇 max 𝑚 𝑐 cos 𝛿 cos 𝛾 𝑣 𝑦 = -𝜇 𝑦 ∥r∥ 3 + 𝑇 max 𝑚 𝑐 cos 𝛿 sin 𝛾 𝑣 𝑧 = -𝜇 𝑧 ∥r∥ 3 + 𝑇 max 𝑚 𝑐 sin 𝛿 𝑚 = - 𝑇 max 𝐼 sp 𝑔 0 𝑐 (1
x = f(x) + g(x, u) (2)
where

f(x) =             0 -𝜇 r ∥r∥ 3 0             , g(x, u) =             0 𝑇 max 𝑚 𝑐n -𝑇 max 𝐼 sp 𝑔 0 𝑐             (3)
The vector

n =             cos 𝛿 cos 𝛾 cos 𝛿 sin 𝛾 sin 𝛿             (4)
represents the thrust direction. In all the above equations, the time dependency has been dropped for brevity.

We consider the problem of finding the spacecraft trajectory that minimizes the fuel consumption from a fixed initial boundary condition to a fixed final boundary condition and no path constraints. The objective function can be expressed in Meyer form as [START_REF] Longuski | Optimal control with aerospace applications[END_REF] 𝐽 = -𝑚(𝑡 𝑓 ) [START_REF] Kirk | Optimal Control Theory: An Introduction[END_REF] where 𝑡 𝑓 is the time of flight. Note that minimizing 𝐽 is equivalent to maximizing the final mass 𝑚(𝑡 𝑓 ). The constraints of the problem are the boundary conditions

x(𝑡 0 ) = x 0 , x(𝑡 𝑓 ) = x 𝑓 (6)
and the control bounds, namely

0 ≤ 𝑐 ≤ 1 (7)
The optimization problem is formulated as min u∈U 𝐽 s.t. Eqs. ( 2), [START_REF] Topputo | Survey of Direct Transcription for Low-Thrust Space Trajectory Optimization with Applications[END_REF], and ( 7) [START_REF] Russell | Primer vector theory applied to global low-thrust trade studies[END_REF] where U is the set of admissible controls. The Hamiltonian of the system is [3]

𝐻(x, u, p) = p 𝑟 • v + p 𝑣 • -𝜇 r ∥r∥ 3 + 𝑇 max 𝑚 𝑐n -𝑝 𝑚 𝑇 max 𝐼 sp 𝑔 0 𝑐 (9)
where p = [p 𝑟 , p 𝑣 , 𝑝 𝑚 ] is the vector that collects the position, velocity, and mass costate variables. According to the Pontryagin's Maximum Principle (PMP), the optimal control u * maximizes the Hamiltonian [START_REF] Kirk | Optimal Control Theory: An Introduction[END_REF]. Therefore,

u * = argmax u∈U (p 𝑣 • n) 𝑇 max 𝑚 -𝑝 𝑚 𝑇 max 𝐼 sp 𝑔 0 𝑐 (10)
Let us indicate the optimal thrust direction with n * , and let us define the switching function

𝑆 = (p 𝑣 • n * ) 𝑇 max 𝑚 -𝑝 𝑚 𝑇 max 𝐼 sp 𝑔 0 (11)
In turn,

𝑐 * = argmax 𝑐∈U 𝑆𝑐 (12)
As a consequence,

𝑐 * =                  1 if 𝑆 > 0 0 if 𝑆 < 0 ∈ (0, 1) if 𝑆 = 0 (13)
The third of the cases in Eq. ( 13) represents the singular case. If p 𝑣 ̸ = 0, then n * = p 𝑣 ∥p 𝑣 ∥ and the switching function 𝑆 becomes

𝑆 = ∥p 𝑣 ∥ 𝑚 - 𝑝 𝑚 𝐼 sp 𝑔 0 (14)
Finally, the dynamics of the costates can also be retrieved from the Hamiltonian function and be written as [START_REF] Bryson | Applied optimal control[END_REF]

                   p 𝑟 = - 𝜕𝐻 𝜕r = - 3𝜇 ∥r∥ 5 (r • p 𝑣 )r + 𝜇 ∥r∥ 3 p 𝑣 p 𝑣 = - 𝜕𝐻 𝜕v = -p 𝑟 𝑝 𝑚 = - 𝜕𝐻 𝜕𝑚 = 𝑐 ∥p 𝑣 ∥𝑇 max 𝑚 2 (15) Lemma 1. Let 𝑆 = 0 on 𝐼 𝑠 ⊂ [𝑡 0 , 𝑡 𝑓 ]. Then, p 𝑣 ̸ = 0 and n * = p 𝑣 ∥p 𝑣 ∥ on 𝐼 𝑠 .
Proof. To prove the claim by contradiction, assume that p 𝑣 = 0. From Eq. ( 11), since 𝑆 = 0, then also 𝑝 𝑚 = 0. Moreover, since p 𝑣 = 0, then p 𝑣 = 0. From the second of Eqs. [START_REF] Park | Necessary conditions for the optimality of singular arcs of spacecraft trajectories subject to multiple gravitational bodies[END_REF], p 𝑟 = 0. However, this would violate the nontriviality condition

[p 𝑟 , p v , 𝑝 𝑚 ] ̸ = 0 on 𝐼 𝑠 .
Therefore, the optimal thrust vector n * and the costate p 𝑣 associated with the spacecraft velocity are always parallel on singular arcs. This result will be used throughout the rest of the paper.

III. General Results

In this section, two of our main findings are presented in the form of theorems, along with several lemmas that will be used to proof the theorems.

A. Main Statements

Theorem 1. Let a LTO problem be described by Eq. [START_REF] Russell | Primer vector theory applied to global low-thrust trade studies[END_REF], with the further assumption of planar dynamics. Moreover, let 𝑆 = 0 on 𝐼 𝑠 for that problem. If the angle between the thrust direction n * and the spacecraft position vector r is denoted as 𝛽, then the closed-form surface that relates 𝛽 to the eccentricity 𝑒 and the true anomaly 𝜃 of the spacecraft on 𝐼 𝑠 is given by

Ψ(𝑒, 𝜃, 𝛽) = 2 cos 𝛽 sin 𝛽 -1± √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 ∓ √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 -(1 -3 cos 2 𝛽) 𝑒 sin 𝜃 1 + 𝑒 cos 𝜃 = 0 ( 16
)
Remark 1. As it will be shown later on, the signs ± and ∓ come from different sources and can change independently.

Therefore, the two terms do not cancel out a priori.

Equation ( 16) represents a set of necessary conditions to have singular arcs. This means that the angle between the optimal thrust direction and the position vector must assume precise values depending on the value of the eccentricity and true anomaly of the spacecraft orbit.

2) the term 𝐴(𝛽) can only be zero when either cos 𝛽 = 0 or sin 𝛽 = ± √ 2 5 . The first case corresponds to the conditions 𝑒 = 0 or 𝜃 = 2𝑘 𝜋, 𝑘 ∈ Z. Moreover, for that case, a. if 𝑒 = 0, there might exist an interval Ĩ𝑠 ∈ 𝐼 𝑠 such that 𝐴(𝛽) = 0 on Ĩ𝑠 if cos 𝜃 = 0; b. if 𝜃 = 2𝑘 𝜋, 𝑘 ∈ Z, then 𝐴(𝛽) = 0 only at isolated points, i.e., Ĩ𝑠 ∈ 𝐼 𝑠 such that 𝐴(𝛽) = 0 on Ĩ𝑠 .

The terms 𝐵 and 𝐴 will be given explicitly in the proof of the theorem. They only depend on the physical variables ∥r∥, 𝛽, 𝑒, 𝜃, and 𝑚. Therefore, the computation of the singular control 𝑐 𝑠 can be easily performed in case singular arcs appear while solving a LTO problem. Nonetheless, we will show that the factor 𝑐 𝑠 is actually singular (i.e, 0 < 𝑐 𝑠 < 1) in limited regions of interest of the state space.

B. Main Lemmas and Corollaries

Lemma 2. If 𝑆 = 0 on 𝐼 𝑠 , then the costates p 𝑣 an p 𝑟 are perpendicular on 𝐼 𝑠 .

Proof. The result from Lemma 1 allows to write the switching function as in Eq. ( 14). Moreover, it must also be

𝑆 = - ∥p 𝑣 ∥ 𝑚 2 𝑚 + p 𝑣 • p 𝑣 ∥p 𝑣 ∥𝑚 - 𝑝 𝑚 𝐼 sp 𝑔 0 = 0 (18)
By substituting the expressions for 𝑚 and 𝑝 𝑚 contained in Eqs. ( 1) and ( 15), it is obtained

𝑆 = ∥p 𝑣 ∥ 𝑚 2 𝑇 max 𝐼 sp 𝑔 0 𝑐 + p 𝑣 • p 𝑣 ∥p 𝑣 ∥𝑚 - ∥p 𝑣 ∥ 𝑚 2 𝑇 max 𝐼 sp 𝑔 0 𝑐 = p 𝑣 • p 𝑣 ∥p 𝑣 ∥𝑚 (19)
Using the second costate equation in Eq. ( 15), one gets

𝐷 1 = p 𝑣 • p 𝑣 = -p 𝑣 • p 𝑟 = 0 (20)
Therefore, the costates p 𝑣 and p 𝑟 are perpendicular along singular arcs.

Lemma 3. Let a LTO problem be described by Eq. ( 8), with the further assumption of planar dynamics. Moreover, let 𝛼 indicate the angle between the reference direction and the spacecraft radius r. If 𝑆 = 0 on 𝐼 𝑠 , then on 𝐼 𝑠 it must be

-∥r∥ 3 ( 𝛼 + 𝛽) 2 + 𝜇(1 -3 cos 2 𝛽) = 0 (21)
Proof. Using the result of Lemma 2, from the hypothesis 𝑆 = 0 follows that 𝐷 1 = 0 (see Eq. ( 20)). Consequently,

𝐷 1 = p 𝑣 • p 𝑟 + p 𝑣 • p 𝑟 = -p 𝑟 • p 𝑟 + p 𝑣 • - 3𝜇 ∥r∥ 5 (r • p 𝑣 )r + 𝜇 ∥r∥ 3 p 𝑣 = -∥p 𝑟 ∥ 2 - 3𝜇 ∥r∥ 5 (r • p 𝑣 )(r • p 𝑣 ) + 𝜇 ∥r∥ 3 (p 𝑣 • p 𝑣 ) = 0 (22)
Introducing the hypothesis of planar dynamics, the angles described in Fig. 1 can be defined. Therefore,

-∥p 𝑟 ∥ 2 - 3𝜇 ∥r∥ 3 ∥p 𝑣 ∥ 2 cos 2 𝛽 + 𝜇 ∥r∥ 3 ∥p 𝑣 ∥ 2 = 0 (23)
From which comes

-∥r∥ 3 ∥p 𝑟 ∥ 2 +𝜇∥p 𝑣 ∥ 2 (1 -3 cos 2 𝛽) = 0. ( 24
)
Consider the costate p 𝑣 = ∥p 𝑣 ∥q, where ∥q∥= 1. Deriving the vector with respect to time yields

p 𝑣 = d d𝑡 ∥p 𝑣 ∥q + ∥p 𝑣 ∥ q = d d𝑡 ∥p 𝑣 ∥q + ∥p 𝑣 ∥ 𝛿s, ( 25
)
where s is a unitary vector perpendicular to q. The dynamics of the velocity costate in Eqs. [START_REF] Park | Necessary conditions for the optimality of singular arcs of spacecraft trajectories subject to multiple gravitational bodies[END_REF] and Eq. [START_REF] Junkins | Exploration of alternative state vector choices for low-thrust trajectory optimization[END_REF] show that the component of p 𝑣 along q must be zero. In turn,

p 𝑣 = -p 𝑟 = ∥p 𝑣 ∥ 𝛿s (
From Fig. 1, it can be observed that the angle 𝛿 = 𝛼 + 𝛽, so

∥p 𝑟 ∥= ∥p 𝑣 ∥| 𝛼 + 𝛽| (28)
Finally, Eq. ( 24) becomes

0 = -∥r∥ 3 ∥p 𝑣 ∥ 2 ( 𝛼 + 𝛽) 2 + 𝜇∥p 𝑣 ∥ 2 (1 -3 cos 2 𝛽) = -∥r∥ 3 ( 𝛼 + 𝛽) 2 + 𝜇(1 -3 cos 2 𝛽) = 𝐷 2 (29)
Corollary 1. In order to have singular arcs, it must be verified

1 -3 cos 2 𝛽 ≥ 0 (30)
Proof. The proof directly follows from Lemma 3.

8 and

∥ p 𝑣 ∥= ∥p 𝑟 ∥= ∥p 𝑣 ∥| 𝛿| (27)
From Fig. 1, it can be observed that the angle 𝛿 = 𝛼 + 𝛽, so

∥p 𝑟 ∥= ∥p 𝑣 ∥| 𝛼 + 𝛽| (28)
Finally, Eq. ( 24) becomes

0 = -∥r∥ 3 ∥p 𝑣 ∥ 2 ( 𝛼 + 𝛽) 2 + 𝜇∥p 𝑣 ∥ 2 (1 -3 cos 2 𝛽) = -∥r∥ 3 ( 𝛼 + 𝛽) 2 + 𝜇(1 -3 cos 2 𝛽) = 𝐷 2 (29)
Corollary 1. In order to have singular arcs, it must be verified

1 -3 cos 2 𝛽 ≥ 0 (30)
Proof. The proof directly follows from Lemma 3.

Remark 2. The condition in Corollary 1 corresponds to the following values of 𝛽:

𝐼 𝛽 = 𝐼 1 𝛽 ∪ 𝐼 2 𝛽 (31)
where

𝐼 1 𝛽 = {𝛽 : 𝛽 0 + 2𝑘 𝜋 ≤ 𝛽 ≤ 𝜋 -𝛽 0 + 2𝑘 𝜋, 𝑘 ∈ N} 𝐼 2 𝛽 = {𝛽 0 + (2𝑘 + 1)𝜋 ≤ 𝛽 ≤ 2𝜋(𝑘 + 1) -𝛽 0 , 𝑘 ∈ N} (32)
with 𝛽 0 = arccos √

C. Proof of Theorem 1

Proof. In order to have singular arcs, it must be 𝐷 1 = 0 (see Eq. ( 22)). Considering Eq. ( 24),

-2∥r∥ 3 (p 𝑟 • p 𝑟 ) -3∥r∥ ∥p 𝑟 ∥ 2 (r • r) + 2𝜇(1 -3 cos 2 𝛽)(p 𝑣 • p 𝑣 ) + 6𝜇∥p 𝑣 ∥ 2 cos 𝛽 sin 𝛽 𝛽 = 0 (33)
Using the result from Lemma 2, the equation becomes

𝐷 3 = -2∥r∥ 3 (p 𝑟 • p 𝑟 ) -3∥r∥ ∥p 𝑟 ∥ 2 (r • r) + 6𝜇∥p 𝑣 ∥ 2 cos 𝛽 sin 𝛽 𝛽 = 0 (34)
Making use of the costate equation for p 𝑟 , the term (p 𝑟 • p 𝑟 ) can be expressed as

(p 𝑟 • p 𝑟 ) = p 𝑟 • - 3𝜇 ∥r∥ 5 (r • p 𝑣 )r + 𝜇 ∥r∥ 3 p 𝑣 = - 3𝜇 ∥r∥ 5 (r • p 𝑣 )(r • p 𝑟 ) = - 3𝜇 ∥r∥ 3 ∥p 𝑣 ∥ ∥p 𝑟 ∥cos 𝛽 cos 𝛽± 𝜋 2 (35)
In the above equation, Lemma 2 has been used to write (r • p 𝑟 ) as ∥r∥ ∥p 𝑟 ∥cos 𝛽± 𝜋 2 . Note that as also shown in Fig. 1, the orientation of p 𝑟 is not known, and therefore both the signs ± need to be considered. Since cos 𝛽± 𝜋 2 = ∓ sin 𝛽,

(p 𝑟 • p 𝑟 ) = ± 3𝜇 ∥r∥ 3 ∥p 𝑣 ∥ ∥p 𝑟 ∥cos 𝛽 sin 𝛽 = ± 3𝜇 ∥r∥ 3 ∥p 𝑣 ∥ 2 | 𝛿|cos 𝛽 sin 𝛽 ← from Eq. (27) = ± 3𝜇 ∥r∥ 3 ∥p 𝑣 ∥ 2 cos 𝛽 sin 𝛽| 𝛼 + 𝛽| ← from Fig. 1 (36)
By substituting this expression back in Eq. (34):

∓6𝜇∥p 𝑣 ∥ 2 cos 𝛽 sin 𝛽| 𝛼 + 𝛽|-3∥r∥ ∥p 𝑟 ∥ 2 (r • r) + 6𝜇∥p 𝑣 ∥ 2 cos 𝛽 sin 𝛽 𝛽 = 0 ∓6𝜇∥p 𝑣 ∥ 2 cos 𝛽 sin 𝛽| 𝛼 + 𝛽|-3∥r∥ ∥p 𝑣 ∥ 2 ( 𝛼 + 𝛽) 2 (r • r) + 6𝜇∥p 𝑣 ∥ 2 cos 𝛽 sin 𝛽 𝛽 = 0 ∓2𝜇 cos 𝛽 sin 𝛽| 𝛼 + 𝛽|-∥r∥( 𝛼 + 𝛽) 2 (r • r) + 2𝜇 cos 𝛽 sin 𝛽 𝛽 = 0 (37) From Lemma 3, ( 𝛼 + 𝛽) 2 = 𝜇 ∥r∥ 3 (1 -3 cos 2 𝛽) (38)
| 𝛼 + 𝛽|= √ 𝜇 ∥r∥ 3 (1 -3 cos 2 𝛽) (39) 𝛽 = -𝛼± √ 𝜇 ∥r∥ 3 (1 -3 cos 2 𝛽) (40)
Now, we want to express 𝛼 as a function of physical variables only, such as orbital parameters. In case of planar transfers, the angle 𝛼 can be expressed as [START_REF] Curtis | Orbital mechanics for engineering students[END_REF]

𝛼 = 𝜃 + 𝜔, (41)
where 𝜔 is the argument of periapsis, which is defined as [START_REF] Curtis | Orbital mechanics for engineering students[END_REF] 𝜔 = arccos

N • e ∥N∥ ∥e∥ (42)
N is a vector pointing towards the ascending node of the orbit, and e is the eccentricity vector. In case of a planar orbit, N is undefined. By convention, it is assumed that it coincides with the reference direction 𝑥 and therefore the angle 𝜔 is defined as the angle between the reference direction and the eccentricity vector, hence Eq. (42). By deriving Eq. ( 41):

𝛼 = 𝜃 + 𝜔 = 𝜃 2B + 𝜃 P + 𝜔 2B + 𝜔 P ( 43
)
where the terms related to the two-body motion and to orbital perturbations have been highlighted and identified with (•) 2B and (•) P , respectively. In case of planar orbits, it can be proved that 𝜔 P = -𝜃 P [START_REF] Curtis | Orbital mechanics for engineering students[END_REF]. Moreover, since in two-body motion the argument of periapsis does not change, 𝜔 2B = 0. Therefore,

𝛼 = 𝜃 2B = ∥h∥ ∥r∥ 2 = √ 𝜇 ∥r∥ 3 (1 + 𝑒 cos 𝜃) (44)
where

∥h∥= √ 𝜇∥r∥(1 + 𝑒 cos 𝜃) (45)
is the norm of the orbital specific angular momentum [START_REF] Curtis | Orbital mechanics for engineering students[END_REF]. Therefore, one can write 𝛽 as

𝛽 = √ 𝜇 ∥r∥ 3 (1 + 𝑒 cos 𝜃) -1± √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 (46)
Moreover, we have that

(r • r) = (r • v) = ∥r∥𝑣 𝑟 = ∥r∥ 𝜇 ∥h∥ 𝑒 sin 𝜃 = 𝑒 sin 𝜃 √ 𝜇∥r∥ 1 + 𝑒 cos 𝜃 , (47)
where

𝑣 𝑟 = 𝜇 ∥h∥ 𝑒 sin 𝜃 (48)
is the component of the velocity along the radius direction [START_REF] Curtis | Orbital mechanics for engineering students[END_REF]. By substituting Eqs. ( 38), ( 39), ( 46) and (47) in Eq.

(37) we obtain

Ψ(𝑒, 𝜃, 𝛽) = 2 cos 𝛽 sin 𝛽 -1± √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 ∓ √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 -(1 -3 cos 2 𝛽) 𝑒 sin 𝜃 1 + 𝑒 cos 𝜃 = 0 (49)
Figure 2 shows Ψ as a function of 𝛽 for 𝑒 = 0.2 and 𝜃 = 5 12 𝜋 for the different combinations of the signs.

Remark 3. It can be observed that:

1) the function is only defined when 1 -3 cos 2 𝛽 ≥ 0, and therefore in the two sub-domains of [0, 2𝜋] defined in Remark 2;

2) since cos(𝛽 + 𝜋) = -cos 𝛽 and sin(𝛽 + 𝜋) = -sin 𝛽, the function has the same behaviour for both the sub-domains.

3) when 𝛽 = 𝛽 0 , the function Ψ = 2(cos 2 𝛽 0 -sin 2 𝛽 0 ) < 0 and when 𝛽 = 𝜋 -𝛽 0 then Ψ = 2[cos 2 (𝜋 -𝛽 0 ) -sin 2 (𝜋 -

𝛽 0 )] > 0.
Corollary 2. For a fixed couple ē, θ , the equation Ψ ē, θ, 𝛽 = 0 defined in Theorem 1 has at least 6 zeros and at most 10.

Proof. The function Ψ ē, θ, 𝛽 = 0 should be studied in the interval 𝐼 𝛽 . However, due to the second point of Remark 3, it can be studied inside the interval 𝐼 1 𝛽 and the same results apply for 𝐼 2 𝛽 . Let us compute the derivative of the function Ψ( ē, θ, 𝛽) with respect to 𝛽:

dΨ d𝛽 ( ē, θ, 𝛽) = 2(cos 2 𝛽 -sin 2 𝛽) -1± √ 1 -3 cos 2 𝛽 1 + ē cos θ ∓ √ 1 -3 cos 2 𝛽 1 + ē cos θ + 6 cos 2 𝛽 sin 2 𝛽 √ 1 -3 cos 2 𝛽 ± 1 √ 1 + ē cos θ ∓ 1 √ 1 + ē cos θ -6 cos 𝛽 sin 𝛽 𝑒 sin θ 1 + ē cos θ (50)
Now, depending on the signs inside the function, we have three cases, correspondent to the three cases in Fig. 2.

Opposite Signs

Let us first consider the case in which the two signs are opposite. We get:

dΨ 2 d𝛽 = -2(cos 2 𝛽 -sin 2 𝛽) -6 cos 𝛽 sin 𝛽 𝑒 sin θ 1 + ē cos θ = -2(cos 2 𝛽 -sin 2 𝛽) -6 cos 𝛽 sin 𝛽Γ, Γ = const. ( 51
)
Let us analyze the first term of the function, i.e., 𝑇 1 = -2(cos 2 𝛽 -sin 2 𝛽). Inside 𝐼 1 𝛽 , it is always positive, has a maximum in 𝛽 = 𝜋/2, and it is symmetric with respect to the 𝛽 = 𝜋/2 axis. The term 𝑇 2 = -6 cos 𝛽 sin 𝛽Γ inside the same interval is instead monotonic and 1. negative in 𝐼 (1,1) 𝛽 = {𝛽 0 ≤ 𝛽 < 𝜋/2} and positive in the interval 𝐼 (1,2) Case Γ = 0 In this case, 𝑇 2 ≡ 0 and the derivative of Ψ is always positive. It follows that Ψ has one root only.

𝛽 = {𝜋/2 < 𝛽 ≤ 𝜋 -𝛽 0 } if Γ > 0;

Both Signs Negative

Consider now the case when both signs are negative. The derivative of Ψ takes the form

dΨ 3 d𝛽 = -2(cos 2 𝛽 -sin 2 𝛽) 1 + 2 √ 1 -3 cos 2 𝛽 1 + ē cos θ -12 cos 2 𝛽 sin 2 𝛽 √ 1 -3 cos 2 𝛽 √ 1 + ē cos θ -6 cos 𝛽 sin 𝛽Γ = 𝑇 1 + 𝑇 2 + 𝑇 3 (52)
First note that when 𝛽 → 𝛽 0 and 𝛽 → 𝜋 -𝛽 0 , then dΨ 3 d𝛽 → -∞. Let us consider two values 𝛽 1 , 𝛽 2 such that 1) dΨ 3 d𝛽 < 0 in (𝛽 0 , 𝛽 1 ) and dΨ 3 d𝛽 = 0 at 𝛽 = 𝛽 1 < 𝜋/2;

2) dΨ 3 d𝛽 < 0 in (𝛽 2 , 𝜋 -𝛽 0 ) and dΨ 3 d𝛽 = 0 at 𝛽 = 𝛽 2 > 𝜋/2;

3) dΨ 3 d𝛽 can change sign in (𝛽 1 , 𝛽 2 ).

We want to show that dΨ 3 d𝛽 > 0 always in (𝛽 1 , 𝛽 2 ), from which would follow that Ψ 3 has only one root. The following statements hold in the interval 𝐼 1 𝛽 : 1) 𝑇 1 (𝛽) = 𝑇 1 (𝜋 -𝛽), 𝑇 2 (𝛽) = 𝑇 2 (𝜋 -𝛽), and 𝑇 3 (𝛽) = -𝑇 3 (𝜋 -𝛽);

2) 𝑇 1 is always positive and has a maximum at 𝛽 = 𝜋/2;

3) 𝑇 2 is always negative and it is zero at 𝛽 = 𝜋/2.

Both Signs Positive

Finally consider the case when both of the signs inside Ψ are positive. The derivative becomes:

dΨ 1 d𝛽 = -2(cos 2 𝛽 -sin 2 𝛽) 1 -2 √ 1 -3 cos 2 𝛽 1 + ē cos θ + 12 cos 2 𝛽 sin 2 𝛽 √ 1 -3 cos 2 𝛽 √ 1 + ē cos θ -6 cos 𝛽 sin 𝛽Γ = 𝑇 1 + 𝑇 2 + 𝑇 3 (53)
First note that when 𝛽 → 𝛽 0 and 𝛽 → 𝜋 -𝛽 0 , then Ψ 1 → +∞. Let us define two angles 𝛽 1 , 𝛽 2 such that 1) dΨ 1 d𝛽 > 0 in (𝛽 0 , 𝛽 1 ) and dΨ 3 d𝛽 = 0 at 𝛽 = 𝛽 1 < 𝜋/2;

2) dΨ 1 d𝛽 > 0 in (𝛽 2 , 𝜋 -𝛽 0 ) and dΨ 1 d𝛽 = 0 at 𝛽 = 𝛽 2 > 𝜋/2;

3) dΨ 1 d𝛽 can change sign in (𝛽 1 , 𝛽 2 ).

The idea is to show that dΨ 1 d𝛽 < 0 inside (𝛽 1 , 𝛽 2 ). If this happens, then Ψ has at most three zeros and at least one. The following statements hold in the interval 𝐼 1 𝛽 : 1) 𝑇 1 (𝛽) = 𝑇 1 (𝜋 -𝛽), 𝑇 2 (𝛽) = 𝑇 2 (𝜋 -𝛽), and 𝑇 3 (𝛽) = -𝑇 3 (𝜋 -𝛽);

2) 𝑇 1 has a negative minimum at 𝛽 = 𝜋/2;

3) 𝑇 2 is always positive and it is zero at 𝛽 = 𝜋/2. 4) 𝑇 3 has already been analyzed before.

Regardless of the sign of the term 𝑇 3 , the angle 𝛽 1 < 𝜋/2 due to Properties 2), 3), and 4) of the above item list.

Case Γ > 0 Due to the monotonicity of the functions 𝑇 1 and 𝑇 2 in the interval (𝛽 1 , 𝜋/2) and since 𝑇 3 < 0 inside the same interval, we have that dΨ 1 d𝛽 ≤ 0 in the interval [𝛽 1 , 𝜋/2). In 𝛽 = 𝜋/2, 𝑇 1 < 0, 𝑇 2 = 𝑇 3 = 0, therefore the derivative is negative. Since 𝑇 1 , 𝑇 2 , and 𝑇 3 are monotonic, dΨ 1 d𝛽 is negative and increasing in the interval [𝜋/2, 𝛽 2 ). Consequently, Ψ can have maximum 3 zeros and minimum 1, depending on the values of 𝑇 1 , 𝑇 2 , and 𝑇 3 at 𝛽 = 𝛽 1 and 𝛽 = 𝛽 2 .

Case Γ < 0 Due to the monotonicity of the functions 𝑇 1 , 𝑇 2 , and 𝑇 3 in the interval (𝛽 1 , 𝜋/2), we have that dΨ 1 d𝛽 is negative and increasing in the interval [𝛽 1 , 𝜋/2). In 𝛽 = 𝜋/2, 𝑇 1 < 0, 𝑇 2 = 𝑇 3 = 0, therefore the derivative is negative.

Since 𝑇 1 , 𝑇 2 are monotonic and 𝑇 3 < 0 in the interval [𝜋/2, 𝛽 2 ), then dΨ 1 d𝛽 is negative and increasing. Consequently, Ψ can have maximum 3 zeros and minimum 1, depending on the values of 𝑇 1 , 𝑇 2 , and 𝑇 3 at 𝛽 = 𝛽 1 and 𝛽 = 𝛽 2 .

Case Γ = 0 The same conclusion of the first two cases directly follows.

Therefore, the function Ψ can have at most 5 zeros and at least 3 inside interval 𝐼 1 𝛽 . The same reasoning applies to the interval 𝐼 2 𝛽 , hence the total number of zeros of the function Ψ is at most 10 and at least 6. Corollary 2 provides insightful information about the rareness of singular arcs: for a given couple ( ē, θ), there only is a discrete set of angles 𝛽 that satisfy the necessary condition. Therefore, the thruster must assume a specific direction, otherwise no singular arcs can happen. Not only this must be verified, but the tuple (𝑒, 𝜃, 𝛽) must move on a surface.

D. Proof of Theorem 2.1

Proof. Let us consider Eq. (34). Using the expression for (p 𝑟 • p 𝑟 ) from Eq. ( 36), it becomes

𝐷 3 = ∓2𝜇∥p 𝑣 ∥ ∥p 𝑟 ∥cos 𝛽 sin 𝛽 -∥r∥ ∥p 𝑟 ∥ 2 (r • v) + 2𝜇∥p 𝑣 ∥ 2 cos 𝛽 sin 𝛽 𝛽 = 0 (54)
By deriving the above equation, one obtains

𝐷 3 = ∓2𝜇 (p 𝑣 • p 𝑣 ) ∥p 𝑣 ∥ ∥p 𝑟 ∥cos 𝛽 sin 𝛽 + (p 𝑟 • p 𝑟 ) ∥p 𝑟 ∥ ∥p 𝑣 ∥cos 𝛽 sin 𝛽 + ∥p 𝑣 ∥ ∥p 𝑟 ∥(cos 2 𝛽 -sin 2 𝛽) 𝛽 - (r • v) 2 ∥r∥ ∥p 𝑟 ∥ 2 -2∥r∥(p 𝑟 • p 𝑟 )(r • v) -∥r∥ ∥p 𝑟 ∥ 2 [∥v∥ 2 +(r • v)] + 2𝜇 2(p 𝑣 • p 𝑣 ) cos 𝛽 sin 𝛽 𝛽 + ∥p 𝑣 ∥ 2 (cos 2 𝛽 -sin 2 𝛽) 𝛽 2 + ∥p 𝑣 ∥ 2 cos 𝛽 sin 𝛽 𝛽 = 0 (55) Next consider that r • v = r • -𝜇 r ∥r∥ 3 + 𝑇 max 𝑚 𝑐n - 𝜇 ∥r∥ + 𝑇 max 𝑚 𝑐∥r∥cos 𝛽 (56)
Using Lemma 2 and Eqs. ( 36), ( 39) and ( 56), the equation 𝐷 3 = 0 becomes According to the classical orbital mechanics, the component of the velocity along the radius is expressed in Eq. ( 48), whereas the component of the velocity along the direction perpendicular to the radius is [START_REF] Curtis | Orbital mechanics for engineering students[END_REF]

𝐷 3 = - 6𝜇 2 ∥r∥ 3 cos 2 𝛽 sin 2 𝛽∓ 2𝜇 2 ∥r∥ 3 √ (1 -3 cos 2 𝛽)(1 + 𝑒 cos 𝜃)(cos 2 𝛽 -sin 2 𝛽) -1± √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 - 𝜇 ∥r∥ 2 (1 -3 cos 2 𝛽)∓ 6𝜇 ∥r∥ 2 cos 𝛽 sin 𝛽 √ 𝜇 ∥r∥ 3 (1 -3 cos 2 𝛽)(r • v) - 𝜇 ∥r∥ 2 (1 -3 cos 2 𝛽)[∥v∥ 2 - 𝜇 ∥r∥ + 𝑇 max 𝑚 𝑐∥r∥cos 𝛽] + 2 𝜇 2 ∥r∥ 3 (cos 2 𝛽 -sin 2 𝛽)(1 + 𝑒 cos 𝜃) -1± √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃
𝑣 ⊥ = 𝜇 ∥h∥ (1 + 𝑒 cos 𝜃) = 𝜇 √ 𝜇∥r∥(1 + 𝑒 cos 𝜃) (1 + 𝑒 cos 𝜃) = √ 𝜇(1 + 𝑒 cos 𝜃) ∥r∥ ( 58
)
From the above equations we have that

∥v∥ 2 = 𝜇 ∥r∥ 𝑒 2 sin 2 𝜃 1 + 𝑒 cos 𝜃 + 𝜇(1 + 𝑒 cos 𝜃) ∥r∥ = 𝜇 ∥r∥ (1 + 𝑒 cos 𝜃) 2 + 𝑒 2 sin 2 𝜃 (1 + 𝑒 cos 𝜃) (59)
By substituting Eqs. ( 47) and (59) into Eq. ( 57), we get

𝐷 3 = - 6𝜇 ∥r∥ 2 cos 2 𝛽 sin 2 𝛽∓ 2𝜇 ∥r∥ 2 √ (1 -3 cos 2 𝛽)(1 + 𝑒 cos 𝜃)(cos 2 𝛽 -sin 2 𝛽) -1± √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 - 𝜇 ∥r∥ 2 𝑒 2 sin 2 𝜃 (1 + 𝑒 cos 𝜃) (1 -3 cos 2 𝛽)∓ 6𝜇 ∥r∥ 2 cos 𝛽 sin 𝛽 √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 (𝑒 sin 𝜃) - 𝜇 ∥r∥ 2 (1 -3 cos 2 𝛽) (1 + 𝑒 cos 𝜃) 2 + 𝑒 2 sin 2 𝜃 (1 + 𝑒 cos 𝜃) -1 -(1 -3 cos 2 𝛽) 𝑇 max 𝑚 𝑐 cos 𝛽 + 2 𝜇 ∥r∥ 2 (cos 2 𝛽 -sin 2 𝛽)(1 + 𝑒 cos 𝜃) -1± √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 2 + 2 cos 𝛽 sin 𝛽 𝛽 = 0 (60)
Therefore, we have expressed the term 𝐷 3 only as a function of ∥r∥, 𝑒, 𝛽, 𝜃, and 𝛽. In order to cancel out the dependency on 𝛽, we proceed as follows. Let us consider Eq. ( 29) and take its derivative

-3∥r∥(r • v)( 𝛼 + 𝛽) 2 -2( 𝛼 + 𝛽)( 𝛼 + 𝛽)∥r∥ 3 +6𝜇 cos 𝛽 sin 𝛽 𝛽 = 0 (61)
By substituting the expressions for 𝛽, ( 𝛼 + 𝛽), and (r • v):

-3 √ 𝜇 3 ∥r∥ 3 𝑒 sin 𝜃 √ 1 + 𝑒 cos 𝜃 (1 -3 cos 2 𝛽)∓2 √ 𝜇∥r∥ 3 (1 -3 cos 2 𝛽)( 𝛼 + 𝛽) + 6𝜇 cos 𝛽 sin 𝛽 √ 𝜇 ∥r∥ 3 √ 1 + 𝑒 cos 𝜃 -1± √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 = 0 (62)
From Eq. (44),

𝛼 = d d𝑡 ∥h∥ ∥r∥ 2 = d d𝑡 √ 𝜇∥r∥(1 + 𝑒 cos 𝜃) ∥r∥ 2 ← using Eq. (45) = √ 𝜇 d d𝑡 √ (1 + 𝑒 cos 𝜃) ∥r∥ 3 = √ 𝜇 - 3 2 ∥r∥ -7/2 (r • v) √ 1 + 𝑒 cos 𝜃 + 𝑒 cos 𝜃 -𝑒 𝜃 sin 𝜃 2 √ 1 + 𝑒 cos 𝜃 ∥r∥ -3/2 = - 3 2 𝜇 ∥r∥ 3 𝑒 sin 𝜃 + √ 𝜇 ∥r∥ 3 𝑒 cos 𝜃 -𝑒 𝜃 sin 𝜃 2 √ 1 + 𝑒 cos 𝜃 ← using Eq. (47) = - 3 2 𝜇 ∥r∥ 3 𝑒 sin 𝜃 + 𝐾(∥r∥, 𝑒, 𝜃, 𝑒, 𝜃) (63)
From the classical orbital mechanics [START_REF] Curtis | Orbital mechanics for engineering students[END_REF], we have that the evolution in time of the eccentricity and the true anomaly is described by

𝑒 = √ ∥r∥(1 + 𝑒 cos 𝜃) 𝜇 sin 𝜃 𝑇 𝑟 𝑚 + √ ∥r∥ 𝜇(1 + 𝑒 cos 𝜃) [(2 + 𝑒 cos 𝜃) cos 𝜃 + 𝑒] 𝑇 𝑠 𝑚 𝜃 = √ 𝜇 ∥r∥ 3 (1 + 𝑒 cos 𝜃) + √ ∥r∥ 𝑒𝜇(1 + 𝑒 cos 𝜃) (1 + 𝑒 cos 𝜃) cos 𝜃 𝑇 𝑟 𝑚 -(2 + 𝑒 cos 𝜃) sin 𝜃 𝑇 𝑠 𝑚 (64)
In the above equations, 𝑇 𝑟 and 𝑇 𝑠 are the components of the thrust along the radius and its perpendicular direction, respectively. That is,

𝑇 𝑟 = 𝑐𝑇 max cos 𝛽 𝑇 𝑠 = 𝑐𝑇 max sin 𝛽 (65)
Many works in the literature use the Gauss variational equations or some modified version to represent the equations of motion of a spacecraft equipped with low-thrust engines [START_REF] Hudson | Reduction of low-thrust continuous controls for trajectory dynamics[END_REF][START_REF] Junkins | Exploration of alternative state vector choices for low-thrust trajectory optimization[END_REF][START_REF] Gurfil | Nonlinear feedback control of low-thrust orbital transfer in a central gravitational field[END_REF]. Now, let us analyze the term 𝐾 in Eq. ( 63). Consider

first that 𝐾 = 𝐾 1 + 𝐾 2 ,
where

𝐾 1 = √ 𝜇 ∥r∥ 3 𝑒 cos 𝜃 2 √ 1 + 𝑒 cos 𝜃 𝐾 2 = - √ 𝜇 ∥r∥ 3 𝑒 𝜃 sin 𝜃 2 √ 1 + 𝑒 cos 𝜃 (66)
First consider 𝐾 1 . If we substitute the expression for 𝑒, we get

𝐾 1 = √ 𝜇 ∥r∥ 3 𝑒 cos 𝜃 2 √ 1 + 𝑒 cos 𝜃 = √ 𝜇 ∥r∥ 3 cos 𝜃 2 √ 1 + 𝑒 cos 𝜃 √ ∥r∥(1 + 𝑒 cos 𝜃) 𝜇 sin 𝜃 𝑇 𝑟 𝑚 + √ ∥r∥ 𝜇(1 + 𝑒 cos 𝜃) [(2 + 𝑒 cos 𝜃) cos 𝜃 + 𝑒] 𝑇 𝑠 𝑚 = cos 𝜃 sin 𝜃 2∥r∥ 𝑇 𝑟 𝑚 + cos 2 𝜃(2 + 𝑒 cos 𝜃) 2∥r∥(1 + 𝑒 cos 𝜃) 𝑇 𝑠 𝑚 + 𝑒 cos 𝜃 2∥r∥(1 + 𝑒 cos 𝜃) 𝑇 𝑠 𝑚 = cos 𝜃 sin 𝜃 2∥r∥ 𝑐 𝑇 max 𝑚 cos 𝛽 + cos 2 𝜃(2 + 𝑒 cos 𝜃) 2∥r∥(1 + 𝑒 cos 𝜃) 𝑐 𝑇 max 𝑚 sin 𝛽 + 𝑒 cos 𝜃 2∥r∥(1 + 𝑒 cos 𝜃) 𝑐 𝑇 max 𝑚 sin 𝛽 ← from Eq. (65) = 𝑐 𝑇 max 𝑚 1 2∥r∥ cos 𝜃 sin 𝜃 cos 𝛽 + cos 2 𝜃(2 + 𝑒 cos 𝜃) 1 + 𝑒 cos 𝜃 sin 𝛽 + 𝑒 cos 𝜃 1 + 𝑒 cos 𝜃 sin 𝛽 = 𝑐 𝑇 max 𝑚 𝐻 1 (∥r∥, 𝑒, 𝜃, 𝛽) (67)
We can proceed in the same way for 𝐾 2 :

𝐾 2 = - √ 𝜇 ∥r∥ 3 𝑒 𝜃 sin 𝜃 2 √ 1 + 𝑒 cos 𝜃 = - √ 𝜇 ∥r∥ 3 𝑒 sin 𝜃 2 √ 1 + 𝑒 cos 𝜃 √ 𝜇 ∥r∥ 3 (1 + 𝑒 cos 𝜃) - √ 𝜇 ∥r∥ 3 𝑒 sin 𝜃 2 √ 1 + 𝑒 cos 𝜃 √ ∥r∥ 𝑒𝜇(1 + 𝑒 cos 𝜃) (1 + 𝑒 cos 𝜃) cos 𝜃 𝑇 𝑟 𝑚 -(2 + 𝑒 cos 𝜃) sin 𝜃 𝑇 𝑠 𝑚 = - 𝜇 ∥r∥ 3 𝑒 sin 𝜃 2 - sin 𝜃 cos 𝜃 2∥r∥ 𝑇 𝑟 𝑚 + sin 2 𝜃(2 + 𝑒 cos 𝜃) 2∥r∥(1 + 𝑒 cos 𝜃) 𝑇 𝑠 𝑚 ← from Eq. (65) = - 𝜇 ∥r∥ 3 𝑒 sin 𝜃 2 + 𝑐 𝑇 max 𝑚 1 2∥r∥ -cos 𝜃 sin 𝜃 cos 𝛽 + sin 2 𝜃(2 + 𝑒 cos 𝜃) 1 + 𝑒 cos 𝜃 sin 𝛽 = - 𝜇 ∥r∥ 3 𝑒 sin 𝜃 2 + 𝑐 𝑇 max 𝑚 𝐻 2 (∥r∥, 𝑒, 𝜃, 𝛽) (68)
Therefore, using Eqs. ( 66)-( 68), the expression of 𝛼 in Eq. ( 63) becomes:

𝛼 = - 3
Finally, Eq. ( 60) can be rewritten by substituting the expression of 𝛽 as

𝐷 3 = - 6𝜇 ∥r∥ 3 cos 2 𝛽 sin 2 𝛽∓ 2𝜇 ∥r∥ 3 √ 1 -3 cos 2 𝛽(cos 2 𝛽 -sin 2 𝛽) √ 1 + 𝑒 cos 𝜃 -1± √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 - 𝜇 ∥r∥ 3 𝑒 2 sin 2 𝜃 (1 + 𝑒 cos 𝜃) (1 -3 cos 2 𝛽)∓ 6𝜇 ∥r∥ 3 cos 𝛽 sin 𝛽 √ 1 -3 cos 2 𝛽 𝑒 sin 𝜃 √ 1 + 𝑒 cos 𝜃 - 𝜇 ∥r∥ 3 (1 -3 cos 2 𝛽) (1 + 𝑒 cos 𝜃) 2 + 𝑒 2 sin 2 𝜃 (1 + 𝑒 cos 𝜃) -1 - 1 ∥r∥ (1 -3 cos 2 𝛽) 𝑇 max 𝑚 𝑐 cos 𝛽 + 2 𝜇 ∥r∥ 3 (cos 2 𝛽 -sin 2 𝛽)(1 + 𝑒 cos 𝜃) -1± √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 2 + 2 cos 𝛽 sin 𝛽 𝐷(∥r∥, 𝑒, 𝜃, 𝛽) -𝑐 𝑇 max 𝑚 sin 𝛽 ∥r∥ = 0 (73)
The above equation only depends on ∥r∥, 𝑒, 𝑚, 𝜃, 𝛽, and 𝑐. Therefore, we can write

𝑐 𝑠 = 𝐵(∥r∥, 𝑒, 𝑚, 𝜃, 𝛽) 𝐴(𝛽) (74)
where

𝐴 = (1 -3 cos 2 𝛽) cos 𝛽 + 2 cos 𝛽 sin 2 𝛽 (75)
and

𝐵 = ∥r∥ 𝑚 𝑇 max - 6𝜇 ∥r∥ 3 cos 2 𝛽 sin 2 𝛽∓ 2𝜇 ∥r∥ 3 √ 1 -3 cos 2 𝛽(cos 2 𝛽 -sin 2 𝛽) √ 1 + 𝑒 cos 𝜃 -1± √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃 - 𝜇 ∥r∥ 3 𝑒 2 sin 2 𝜃 (1 + 𝑒 cos 𝜃) (1 -3 cos 2 𝛽)∓ 6𝜇 ∥r∥ 3 cos 𝛽 sin 𝛽 √ 1 -3 cos 2 𝛽 𝑒 sin 𝜃 √ 1 + 𝑒 cos 𝜃 - 𝜇 ∥r∥ 3 (1 -3 cos 2 𝛽) (1 + 𝑒 cos 𝜃) 2 + 𝑒 2 sin 2 𝜃 (1 + 𝑒 cos 𝜃) -1 + 2 𝜇 ∥r∥ 3 (cos 2 𝛽 -sin 2 𝛽)(1 + 𝑒 cos 𝜃) -1± √ 1 -3 cos 2 𝛽 1 + 𝑒 cos 𝜃
Now, if the case 𝑒 = 0 is considered,

𝑒 = ±2 √ ∥r∥ 𝜇 cos 𝜃 𝑐𝑇 max 𝑚 (80)
which can only be zero if either 𝑐 = 0 (and thus in a non-singular arc case) or when cos 𝜃 = 0. Consider now the case in which sin 𝜃 = 0. From Eq. ( 81),

𝜃 = √ 𝜇 ∥r∥ 3 (1 + 𝑒 cos 𝜃) > 0 (81)

F. Summary of the Theoretical Results

The theoretical results obtained in Section III can be summarized as follows.

1) We have found the algebraic necessary condition expressed in Eq. ( 16) to have singular arcs that depends on three physical variables only, namely the eccentricity of the spacecraft, its true anomaly, and the angle 𝛽. Previous works in literature have only found necessary conditions that depend on the state and the costates [START_REF] Park | Necessary conditions for the optimality of singular arcs of spacecraft trajectories subject to multiple gravitational bodies[END_REF]. This represents a major improvement because an easier evaluation of the necessary condition can be performed.

Moreover, the evaluation of the algebraic condition provides a physical grasp on the problem.

2) In Corollary 2, we have shown that the solutions of the algebraic necessary conditions for each fixed couple of eccentricity and true anomaly are, in number, between six and ten. This means that the thruster must assume specific directions, hence suggesting the reasons of the rareness of singular arcs.

3) Leveraging the Gauss variational equations, we have expressed the singular throttle factor as the ratio of two algebraic expressions that solely depend on a limited set of physical variables. The expression of the singular throttle factor is provided in Theorem 2. As per the previous point, past works in literature were only able to express it using the costates.

So far, reference works in the field of low-thrust trajectory optimization have assumed singular arcs are negligible when computing interplanetary trajectories, although a complete theoretical framework to justify this assumption was still missing.

IV. Numerical Simulations

A. Algebraic Necessary Conditions

This section shows the solutions of the algebraic condition in Eq. ( 49). Taking into account the fact that two signs can vary, the equation actually represents four conditions. However, due to the symmetry of the condition itself, they reduce to three. Since the equation is not solvable in closed form, the solutions are obtained numerically with the MATLAB® function fzero and the default solver, which requires an initial guess. We relied on Corollary 2 to provide the correct number of initial guesses for each couple (𝑒, 𝜃) such that all the solutions of the equation were found. Figure 3 shows the angle 𝛽 that respects the condition for the case sin 𝛽 > 0. Figure 4 shows the same as the previous figures, respectively, for the case sin 𝛽 < 0. Note that, due to the symmetry of the necessary conditions, the values 𝛽 2 correspondent to the cases sin 𝛽 < 0 are such that 𝛽 2 = 𝛽 1 + 𝜋, where 𝛽 1 are the angles that correspond to the case sin 𝛽 > 0. In the plots, the eccentricity varies between the values 0 and 0.5, as our simulations on several low-thrust trajectories show that it rarely overcomes the value of 0.5. Note that for the selected intervals of 𝑒 and 𝜃, the algebraic necessary condition has always ten zeros. If cases with 𝑒 > 0.5 are considered, it may happen that it has less, as Fig. 5 shows.

B. Value of the Singular Control

In this section, we use Eq. ( 74) to find the value of the singular control for a grid of parameters ∥r∥, 𝑒, 𝑚, 𝜃. We 2) eccentricity values ranging from 𝑒 min = 10 -3 to 𝑒 max = 0.9; the minimum eccentricity value is not zero to avoid the singularity of the term 𝐴(𝛽) presented in Theorem 2.2.

3) a uniform grid of true anomaly values from 𝜃 min = 10 -2 to 𝜃 max = 1.99𝜋; the minimum and the maximum true Sun and the inclination to the ecliptic of the major Solar System celestial bodies. Although we developed our work using the assumption of planar dynamics, the inclination to the ecliptic of the majorities of the celestial bodies is low, therefore it is likely that the results could still be applied to real missions. Finally, when the radius reaches the value of approximately 15, Fig. 6 shows that if the necessary condition to have singular arcs is satisfied, then the control is indeed singular in 100% of the cases because, for those values, 0 ≤ 𝑐 𝑠 ≤ 1.

V. Conclusions

In this work, algebraic necessary conditions to have singular arcs for the planar two-body low-thrust trajectory optimization problem were presented. The approach, which exploits Gauss variational equations, has allowed to express the necessary conditions as a function of three physical variables only. An analytical expression on the singular control has also been found, which only depends on a limited set of physical variables too. We have shown that the necessary condition is only satisfied if the angle between the thrust direction and the spacecraft radius assumes discrete values, which are in number between six and ten. This suggests the rarity of singular arcs. Moreover, through numerical simulations, we have shown that singular arcs can indeed happen but the cases when the associated throttle factor is singular are relatively rare when trajectories in the inner Solar System need to be designed. Although our approach is specific for planar cases, it is likely that the results can also be extended to three-dimensional cases in practical applications due to the small inclination of most Solar System's celestial bodies. Further work will consist of theoretical investigation of such cases.

  ) r = [𝑥, 𝑦, 𝑧] and v = [𝑣 𝑥 , 𝑣 𝑦 , 𝑣 𝑧 ] are the position and velocity vectors of the spacecraft, respectively, and 𝑚 is the spacecraft mass. u = [𝑐, 𝛿, 𝛾] are the controls, where 𝑐 ∈ [0, 1] is the thrust throttle factor and 𝛿 and 𝛾 are the in-and

Fig. 1

 1 Fig. 1 Definition of the angles.

Fig. 1

 1 Fig. 1 Definition of the angles.

2 .

 2 positive in 𝐼(1,1) 𝛽 = {𝛽 0 ≤ 𝛽 < 𝜋/2} and negative in the interval 𝐼(1,2) 𝛽 = {𝜋/2 < 𝛽 ≤ 𝜋 -𝛽 0 } if Γ < 0; 3. constantly zero if Γ = 0. Case Γ > 0 If |𝑇 2 (𝛽 0 )|> |𝑇 1 (𝛽 0 )|,the derivative dΨ 2 d𝛽 is negative at 𝛽 0 . As 𝑇 2 increases, the derivative crosses zero and becomes positive. Since 𝑇 1 > 0 and 𝑇 2 is monotonic and using the third point of Remark 3, the function has exactly one root. If |𝑇 2 (𝛽 0 )|< |𝑇 1 (𝛽 0 )| the derivative is always positive. Therefore, the function Ψ only has one root. Case Γ < 0 In this case, the derivative of Ψ is always positive in 𝐼 (1,1) 𝛽 . The derivative of Ψ can change sign at most once in 𝐼 (1,2) 𝛽 because 𝑇 2 is monotonic and 𝑇 1 > 0 and, using the third point of Remark 3, Ψ at 𝛽 = 𝜋 -𝛽 0 is positive. It follows that the function has only one root.

4 )

 4 𝑇 3 has already been analyzed before. Case Γ > 0 Let us consider the interval (𝛽 1 , 𝜋/2]. Inside the interval, all the terms 𝑇 𝑖 , 𝑖 = 1, . . . , 3 are growing monotonically. Since dΨ 3 d𝛽 (𝛽 1 ) = 0, then dΨ 3 d𝛽 is always positive inside the considered interval. Due to the symmetries of 𝑇 1 , 𝑇

  Algebraic constraints with one sign positive and one sign negative. Algebraic constraints with both signs negative.

Fig. 2

 2 Fig. 2 Algebraic constraint function for a couple of values of 𝑒 and 𝜃.

2 +

 2 2𝜇 cos 𝛽 sin 𝛽 𝛽 = 0 (57) Note that in the above equation the control throttle factor 𝑐 appears. Now, let us analyze the expression of the term ∥v∥ 2 .

2 + 2 2 Proof.

 222 cos 𝛽 sin 𝛽𝐷(∥r∥, 𝑒, 𝜃, 𝛽) For 𝐴(𝛽) to be zero it must be(1 -3 cos 2 𝛽) cos 𝛽 + 2 cos 𝛽 sin 2 𝛽 = 0 cos 𝛽(1 -3 cos 2 𝛽 + 2 sin 2 𝛽) = 0 cos 𝛽(1 -3 + 3 sin 2 𝛽 + 2 sin 2 𝛽) = 0 cos 𝛽(-2 + 5 sin 2 𝛽) = 0 (77) Which means cos 𝛽 = 0 ∨ sin 𝛽 = ± √ 2 5 (78)However, the second case is not compatible with the condition in Corollary 1. Therefore, the only case for which 𝐴 = 0 happens when 𝛽 = 𝜋 2 + 𝑘 𝜋, 𝑘 ∈ Z. According to Eq. (49), this corresponds to the following cases: 𝑒 = 0 or sin 𝜃 = 0, hence the first point of Theorem 2.1. To demonstrate the second point, consider that if 𝛽 = 𝜋 2 + 𝑘 𝜋, 𝑘 ∈ Z, then 𝑇 𝑟 = 0𝑇 𝑠 = ±𝑐𝑇 max .Consequently, Eqs. (64) become 𝑒 = √ ∥r∥ 𝜇(1 + 𝑒 cos 𝜃) [(2 + 𝑒 cos 𝜃) cos 𝜃 + 𝑒] 𝑇 𝑠 𝑚 𝜃 = √ 𝜇 ∥r∥ 3 (1 + 𝑒 cos 𝜃) + √ ∥r∥ 𝑒𝜇(1 + 𝑒 cos 𝜃) -(2 + 𝑒 cos 𝜃) sin 𝜃 𝑇 𝑠 𝑚

Fig. 3

 3 Fig. 3 Solutions of the equation Ψ(𝑒, 𝜃, 𝛽) = 0 when sin 𝛽 > 0.

Fig. 4

 4 Fig. 4 Solutions of the equation Ψ(𝑒, 𝜃, 𝛽) = 0 when sin 𝛽 < 0.

Theorem 2 .

 2 Let 𝑆 = 0 and (𝑒, 𝜃, 𝛽) satisfy Eq. (16) on 𝐼 𝑠 . Moreover, let the term 𝐴(𝛽) = (1 -3 cos 2 𝛽) cos 𝛽 + 2 cos 𝛽 sin 2 𝛽 ̸ = 0 on 𝐼 𝑠 . Then, 1) the singular thrust control action 𝑐 𝑠 on 𝐼 𝑠 can be expressed as

	𝑐 𝑠 =	𝐵(∥r∥, 𝛽, 𝑒, 𝜃, 𝑚) 𝐴(𝛽)

  Let us analyze the term 𝐻 1 (∥r∥, 𝑒, 𝜃, 𝛽) + 𝐻 2 (∥r∥, 𝑒, 𝜃, 𝛽): 𝜃(2 + 𝑒 cos 𝜃) + 𝑒 cos 𝜃 + sin 2 𝜃(2 + 𝑒 cos 𝜃)

	Using Eq. (71), from Eq. (62) we can find
	-6 cos 𝛽 sin 𝛽 𝜇 ∥r∥ 3 𝑒 sin 𝜃 -𝑐 √ 𝛽 = ∓ + 2 𝑇 max 𝜇 3 ∥r∥ 3 (1 + 𝑒 cos 𝜃) -1± √ 2 √ 𝜇∥r∥ 3 (1 -3 cos 2 𝛽) 1-3 cos 2 𝛽 1+𝑒 cos 𝜃 + 3 sin 𝛽 𝑚 ∥r∥ 𝑇 max sin 𝛽 = 𝐷(∥r∥, 𝑒, 𝜃, 𝛽) -𝑐 𝑚 ∥r∥	√	𝜇 3 ∥r∥ 3	𝑒 sin 𝜃 √ 1+𝑒 cos 𝜃	(1 -3 cos 2 𝛽)
	2 = -2	𝜇 ∥r∥ 3 𝑒 sin 𝜃 + 𝑐 𝜇 ∥r∥ 3 𝑒 sin 𝜃 + 𝑐 𝑇 max 𝑇 max 𝑚 𝑚	𝐻 1 (∥r∥, 𝑒, 𝜃, 𝛽) -[𝐻 1 (∥r∥, 𝑒, 𝜃, 𝛽) + 𝐻 2 (∥r∥, 𝑒, 𝜃, 𝛽)] 𝜇 ∥r∥ 3 𝑒 sin 𝜃 2 + 𝑐 𝑇 max 𝑚	𝐻 2 (∥r∥, 𝑒, 𝜃, 𝛽)	(69)
	𝐻 1 (∥r∥, 𝑒, 𝜃, 𝛽) + 𝐻 2 (∥r∥, 𝑒, 𝜃, 𝛽) =	1 2∥r∥	cos 𝜃 sin 𝜃 cos 𝛽 +	cos 2 𝜃(2 + 𝑒 cos 𝜃) 1 + 𝑒 cos 𝜃	sin 𝛽 +	𝑒 cos 𝜃 1 + 𝑒 cos 𝜃	sin 𝛽
		+	1 2∥r∥	-cos 𝜃 sin 𝜃 cos 𝛽 +	sin 2 𝜃(2 + 𝑒 cos 𝜃) 1 + 𝑒 cos 𝜃	sin 𝛽
		=	1 2∥r∥	cos 2 1 + 𝑒 cos 𝜃	sin 𝛽	(70)
		= =	1 2∥r∥ 1 2∥r∥	cos 2 𝜃(2 + 𝑒 cos 𝜃) + 𝑒 cos 𝜃 + (1 -cos 2 𝜃)(2 + 𝑒 cos 𝜃) 1 + 𝑒 cos 𝜃 sin 𝛽 2 sin 𝛽 = ∥r∥	sin 𝛽
	Eq. (69) reduces therefore to	𝛼 = -2	𝜇 ∥r∥ 3 𝑒 sin 𝜃 + 𝑐	𝑇 max 𝑚	sin 𝛽 ∥r∥	(71)

Table 1 Solar System celestial bodies data.

 1 Celestial body Semi-major axis [AU] Inclination to the ecliptic [deg]

	Mercury	0.39	7.00
	Venus	0.72	3.39
	Earth	1.00	0.00
	NEOs	< 1.30	Variable
	Mars	1.52	1.85
	Main belt	2.20-3.20	Variable
	Ceres	2.77	10.60
	Jupiter	5.20	1.30
	Saturn	9.50	2.49
	Uranus	19.20	0.77
	Neptune	30.10	1.77

Acknowledgment A. C. M., C. G., and F. T. acknowledge EXTREMA, a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant Agreement No. 864697). R.B. acknowledges the French National Research Agency who provided support with funding ANR-22-CE46-0006.

anomaly values are not zero and 2𝜋, respectively, to avoid the singularity of the term 𝐴(𝛽) presented in Theorem 2.2.

4) a mass of value 1;

Note that all the above values are dimensionless. The value of 𝛽 to be used comes from the algebraic necessary condition given a couple of parameters (𝑒, 𝜃). For all combinations of parameters, we evaluate the singular control 𝑐 𝑠 in Eq. ( 74).

If 𝑐 𝑠 > 1 -𝜀 or 𝑐 𝑠 < 𝜀, with 𝜀 = 10 -3 , this means that the controls are in fact non-singular. Figure 6 shows the overall probability of encountering singular controls as a function of the radius ∥r∥ and of eccentricity. The current low-thrust missions are usually sent to Mars, Venus, or the Main Asteroid Belt. Although further destinations are possible, most of low-thrust missions use solar arrays to power engines, and the available power decreases with the inverse of the square of the distance from the Sun. Our simulations show that typical eccentricity values for low-thrust interplanetary missions are in the range of 0.1 -0.5. In these ranges, the overall percentage of possible singular controls is below 10%.

Note that this does not mean that singular arcs happen with this frequency, but that in case the necessary algebraic condition is verified, there is a probability of at most 10% to have singular values. Table 1 shows the distance to the


Paul Laffont-Lozes

Florian Salipante

Geraldine Leguelinel-Blache

Catherine Dunyach-Remy

Jean-Philippe Lavigne

Albert Sotto

Romaric Larcher
email: romaric.larcher@chu-nimes.fr

Effect of antimicrobial consumption on Escherichia coli resistance: assessment and forecasting using Dynamic Regression models in a French university hospital (2014-2019)

Keywords: Dynamic regression Antimicrobial resistance Antimicrobial consumption Escherichia coli Antimicrobial stewardship AAPBI, anti-pseudomonal activity penicillin and beta-lactamase inhibitor, AMC, antimicrobial consumption, AMR, antimicrobial resistance, ARIMA, autoregressive integrated moving average, ASP, antimicrobial stewardship programs, AST, antimicrobial susceptibility testing, DDD, defined daily doses, DR, dynamic regression, DRPC, dynamic regression by principal components, EUCAST, European committee on antimicrobial susceptibility testing, LTF, linear transfer function, MAPE, mean absolute percentage error, PBI, penicillin

Effect of Antimicrobial Consumption on Escherichia coli

resistance: assessment and forecasting using Dynamic Regression models in a French university hospital (2014-2019)

Introduction

Antimicrobial resistance (AMR) is increasing worldwide [START_REF] Murray | Global burden of bacterial antimicrobial resistance in 2019: a systematic analysis[END_REF]. AMR leads to increased complexity of care, mortality and hospital length of stay [START_REF] Barrasa-Villar | Impact on morbidity, mortality, and length of stay of hospital-acquired infections by resistant microorganisms[END_REF]. Human activity, particularly antimicrobial consumption (AMC), is one of the main factors responsible for AMR [START_REF] Laxminarayan | Antibiotic resistance-the need for global solutions[END_REF] . As a result, the World Health Organization (WHO) requested that countries monitor AMC and implement antimicrobial steward-ship programs (ASPs) to reduce AMR [START_REF]Global antimicrobial resistance surveillance system (GLASS) report: early implementation 2020[END_REF]. In Europe, after several years of effort, AMR is reaching steady-state or is even starting to decline [START_REF] Peñalva | Decreasing and stabilising trends of antimicrobial consumption and re-sistance in Escherichia coli and Klebsiella pneumoniae in segmented regression analysis, European Union/European Economic Area, 2001 to 2018[END_REF]. However, AMR remains a concern in hospital settings [START_REF]WHO Regional Office for, Europe/European Centre for Disease Prevention and Control Antimicrobial resistance surveillance in Europe 2022-2020 data[END_REF], prompting healthcare authorities to continue AMR surveillance and encourage novel ASPs.

Digital interventions are promising tools for surveillance and ASP development [START_REF] Van Dort | The impact of digital interventions on antimicrobial stewardship in hospitals: a qualitative synthesis of systematic reviews[END_REF]. At hospital level, dynamic regression (DR) models have been used to assess the impact of AMC on AMR [START_REF] Aldeyab | Modelling the impact of antibiotic use and infection control practices on the incidence of hospital-acquired methicillin-resistant Staphylococcus aureus: a time-series analysis[END_REF][START_REF] Gallini | Influence of fluoroquinolone consumption in inpatients and outpatients on ciprofloxacinresistant Escherichia coli in a university hospital[END_REF][START_REF] Monnet | Antimicrobial drug use and methicillin-resistant Staphylococcus aureus , Aberdeen, 1996-20 0 0[END_REF][START_REF] O'riordan | Antimicrobial use and antimicrobial resistance in Enterobacterales and Enterococcus faecium: a time series analysis[END_REF]. DR models are considered the most relevant methodology [START_REF] López-Lozano | Modelling and forecasting antimicrobial resistance and its dynamic relationship to antimicrobial use: a time series analysis[END_REF] and can be used to predict AMR using AMC data [START_REF] Willmann | Time se-ries analysis as a tool to predict the impact of antimicrobial restriction in an-tibiotic stewardship programs using the example of multidrug-resistant Pseu-domonas aeruginosa[END_REF]. However, the predictive ability of these models remains questionable [START_REF] Colson | Quantifying uncertainty about future antimicrobial resistance: Comparing structured expert judgment and statistical forecasting methods[END_REF].

The aim of this study was to describe trends in AMC and AMR rates in Escherichia coli , the most frequently isolated nosocomial pathogen in France [START_REF]WHO Regional Office for, Europe/European Centre for Disease Prevention and Control Antimicrobial resistance surveillance in Europe 2022-2020 data[END_REF], using DR models in a French teaching hospital, and to investigate the correlation between AMC and AMR, and the burden of each antibiotic class on AMR acquisition in E. coli strains. Also assessed was the ability of DR models to predict AMR, with a view to using them in future ASP implementations.

Methods

Study design and settings

A retrospective ecological study in a French teaching hospital was conducted between January 2014 and December 2019, using the ConsoRes software (a tool for monitoring data on AMC and AMR provided by the French Public Health Agency).

The Nimes University Hospital has 1773 beds, including 46 for critical care, 24 for haematooncology, 235 for surgery and 190 for long-term. During the study period, between 41 300 and 50 100 inpatients per year were admitted and the average days of hospitalisation was 55 200 days/month and increased each year up to 57 500 days/month in 2019.

Bacterial samples

All microbiological samples that were positive for E. coli collected during the study period were analysed. Samples obtained from outpatients, or collected in the emergency department or within 48 h from hospital admission were excluded, unless the patient had been hospitalized within the previous 30 days, according to the definition of nosocomial infection [START_REF] Garner | CDC definitions for nosocomial infections, 1988[END_REF] . Paediatric and ambulatory ward samples were also excluded. All samples were deduplicated based on patient identity and antimicrobial phenotype.

Microbiological species and antimicrobial susceptibility testing (AST) results were obtained from the Department of Microbiology using chromogenic medium or mass spectrometry Vitek®MS (bioMérieux, Marcy-l'Etoile, France) and Vitek 2®(AST-372, bioMérieux) automated system and/or disk diffusion method on Mueller-Hinton agar (Bio-Rad, Hercules, CA, USA), according to the European committee on antimicrobial susceptibility testing (EU-CAST) guidelines [START_REF]EUCAST: Clinical breakpoints and dosing of antibiotics n[END_REF]. Microorganisms were classified as susceptible or resistant using the EUCAST breakpoints [START_REF]EUCAST: Clinical breakpoints and dosing of antibiotics n[END_REF]. All isolates classified as intermediate were reclassified as resistant.

The resistance rate to a given class of antimicrobial was defined as the percentage of isolate that was non-susceptible to ≥1 antimicrobial of this class. The resistance rate to each class of antimicrobial was studied monthly.

Antimicrobial consumption

Data on AMC were extracted monthly from the pharmacy soft-ware. The AMC was determined by counting the number of antimicrobials delivered in care units and expressed in defined daily doses (DDD)/10 0 0 inpatient-days according to the ATC/DDD classi-fication [START_REF]Collaborating Centre for Drug Statistics MethodologyATC Index with DDDs[END_REF].

The focus of the study was on the consumption of penicillin extended-spectrum (J01CA), penicillin and beta-lactamase inhibitor (PBI) (J01CR), anti-pseudomonal activity penicillin and beta-lactamase inhibitor (AAPBI) (J01CR05), cephalosporin (J01D), carbapenem (J01DH), monobactam (J01DF), fluoroquinolone (J01MA), aminoglycoside (J01GB) and sulfamethoxazole-trimethoprim (J01EE01). Consumption of fourth-generation cephalosporin (J01DE01) was not pooled with other cephalosporins because of its broader spectrum of activity and was studied separately. The AMC of the emergency department and paediatric wards was excluded. Antimicrobials dispensed to outpatients were also excluded.

Statistical analysis

All statistical analyses were performed with R software (version 4.2.0). P -values < 0.05 were considered statistically significant. Time series analyses were managed with the Forecast package [START_REF] Hyndman | Automatic time series forecasting: The forecast package for R[END_REF] . AMC trends were evaluated by linear regression and Cochran-Armitage tests were used for AMR. Autoregressive integrated moving average (ARIMA) models were used to describe E. coli resistance to several antimicrobial classes. DR models were used to study the relationship between the use of several antimicrobial classes and resistance to a chosen antimicrobial class. The Linear Transfer Function (LTF) was used to link the output series (antimicrobial resistance) to the input series (the use of several antimicrobials) considering possible time lags (delay of antimicrobial effect on AMR acquisition). A dimensionality reduction strategy was used to avoid multicollinearity. More precisely, DR by principal components (DRPC) models were used as described by del Moral & Valderrama [START_REF] Del Moral | A principal component approach to dynamic regression models[END_REF]. This model enables several antimicrobial classes at different time lags to be considered as input to predict the out-put (resistance to a given antimicrobial class). In general, one or two principal components were used to summarise the information of the use of different antimicrobial classes. Cross-correlations were used to describe the temporal relationship between AMC and AMR. The data were split into a train set (2014-2018 period) and a test set (2019 period). The models were built in the train set and R ²(determination coefficient) was reported to evaluate the percentage of the variance of the observed time-series explained by the model. Forecasts were made with the test set and the accuracy of predictions was evaluated using mean absolute percentage error (MAPE). For each monthly plotted time series, a moving average time series with a window of 7 months (the average of the value observed for one month, the three previous months and the three next months) was added to give a better visual representation.

Ethics approval

The Institutional Review Board of Nimes University Hospital ap-proved the study protocol (No. 210307) and waived the need for patient signed consent. This study was conducted according to the guidelines of the Declaration of Helsinki.

Results

Trends in antimicrobial consumption

Trends in AMC are shown in Table 1. Extended-spectrum penicillins, PBIs and cephalosporins were the most commonly pre-scribed antimicrobials. The AMC of half the antimicrobial classes increased: extended-spectrum penicillins (ß= 2.35, P < 0.001), AAPBI (ß= 1.43, P < 0.001), fourth-generation cephalosporins (ß= 0.127, P < 0.0 01), monobactams (ß= 0.03, P = 0.0 02), and sulfonamides (ß= 0.25, P = 0.006). Interestingly, a significant decrease was observed in fluoroquinolone consumption (ß= -2.6, P < 0.001). The AMC of other antimicrobial classes was unchanged

Trends in antimicrobial resistance

A total of 12 721 samples were screened, 7661 of which were included in the study. Among them, 1241 samples (16%) were collected within the first 48 hours of admission from patients dis-charged within the last 30 days. Included samples were provided mainly from rehabilitation units (N = 1214, 16%), long-term care beds (N = 869, 11%), the department of obstetrics and gynaecology (N = 591, 8%), intensive care units (N = 552, 7%) and the department of gerontology (N = 541, 7%).

Annual rates of AMR are presented in Table 2. The highest average rates of AMR were observed for resistance to extended-spectrum penicillins and PBIs (56% and 37%, respectively). Resistance to other beta-lactams ranged from 16% for AAPBI to 11% for cephalosporins, whereas carbapenem resistance was rare, found in only 35 samples (0.4%). Average rates of resistance to sulfonamides, fluoroquinolones and aminoglycosides were 28%, 25% and 8%, respectively.

Rates of resistance to cephalosporins (from 13% to 9%), fluoroquinolones (from 30% to 19%), and aminoglycosides (from 9% to 6%) decreased, whereas the resistance rates to other classes remained stable. Interestingly, decreases in resistance to cephalosporins and fluoroquinolones were correlated (R ²= 63%).

Correlation between antimicrobial consumption and resistance in E. coli

Cross-correlations are reported in Table 3.

Fluoroquinolone resistance was significantly correlated with the use of fourth-generation cephalosporins (lag 0; r = -0.26), AAPBIs (lag 0 to 3; r = -0.52 to -0.39), PBIs (lag 1; r = 0.24) and fluoroquinolones (lag 0 to 5; r = 0.58 to 0.50). The best multivariate DR model included fluoroquinolone and AAPBI consumption and explained around 54% of the variability in resistance using principal component analysis (Figure 1 and Table 4).

Similarly, cephalosporin resistance was significantly correlated with the consumption of fourth-generation cephalosporins (lag 0; r = -0.26), AAPBIs (lag 0 to 2; r= -0.35 to -0.29) and fluoroquinolones (lag 0 to 4; r = 0.34 to 0.29). DR model analysis included these three antimicrobial classes and explained 15% of variability in resistance (Figure 1 and Table 4). PBI resistance was correlated with the use of AAPBIs (lag 4 to 10; R ²= 21% to 28%), fluoroquinolones (lag 0, 1; r = 0.25 and 0.26) and carbapenems (lag 6; r = -0.25). However, the best DR model integrated only PBI consumption and explained 53% of PBI resistance ( Figure 1 and Table 4 ).

Extended-spectrum penicillin resistance was correlated with the consumption of cephalosporins (lag 5, 8; r = -0.26 and -0.27), PBIs (lag 0, 5; r = 0.26 and -0.24) and carbapenems (lag 0; r = 0.25).

Multivariate DR model explained 36% of resistance to this class of antimicrobial ( Figure 1 and Table 4 ).

No correlation was found between AMC of other class and AMR acquisition.

Assessment of DR models as AMR forecasting tools

The comparison of resistance rates predicted for 2019 using the DR models to those observed in 2019 highlighted an error associated with the forecasted value of 12% for fluoroquinolone resistance, 34% for cephalosporin resistance, 20% for PBI resistance and 8% for extendedspectrum penicillin resistance. Curves of observed versus predicted resistance rates are depicted in Figure 1 .

Discussion

This study investigated the correlation between AMC and acquisition of AMR in E. coli strains in a teaching hospital and reported coefficients of correlation r < 0.6. Over a six-year period, AMC in-creased overall and AMR decreased. The DR models highlighted that the decrease in fluoroquinolone use and the increase in AAPBI consumption explained 54% of the decrease in fluoroquinolone resistance and 15% of the decrease in cephalosporin resistance. In addition, 53% of PBI resistance was explained by PBI consumption and 36% of extended-spectrum penicillin resistance by the use of PBIs, cephalosporins and carbapenems, whereas AMR and AMC remained stable for these classes. Moreover, DR models have also demonstrated variable forecasting capacity that could limit their use for ASP conception.

Certain phenomena may have influenced the variation in AMC during the study period. First, local and national ASPs encouraged physicians to increase the use of extended-spectrum penicillins and sulfonamides and to decrease fluoroquinolone consumption [START_REF]Reported decrease in antibiotic consumption across EU/EEA during COVID-19 pandemic[END_REF]. Second, 20 additional intensive care beds and 7 additional haematology beds were opened in the hospital in 2016. This may explain, at least in part, the increased consumption of PBIs, AAPBIs and fourth-generation cephalosporins, as these care units are among the highest consumers of beta-lactams. Also, some prescriptions for fluoroquinolone may have been replaced by prescriptions for AAPBI and fourth-generation cephalosporin. These changes in AMC have induced changes in the AMR of several bacteria, particularly E. coli , also reported nationally [START_REF]Reported decrease in antibiotic consumption across EU/EEA during COVID-19 pandemic[END_REF].

The current study results are in agreement with those of Mahamat et al., who reported that an increase in 1 DDD/10 0 0 inpatient-days of ofloxacin, ciprofloxacin and norfloxacin increased E. coli resistance to ofloxacin by 0.81%, 0.65% and 0.53%, respectively, with lag times of 4 to 6 months [START_REF] Mahamat | Evolution of fluoroquinolone resistance among Escherichia coli urinary tract isolates from a French university hospital: application of the dynamic regression model[END_REF]. Similar correlations between fluoroquinolone consumption and resistance have been reported both in hospital settings and in the community [START_REF] Guo | Antimicrobial resistance surveillance and prediction of Gram-negative bacteria based on antimicrobial consumption in a hospital setting: A 15-year retrospective study[END_REF][START_REF] Vernaz | Modelling the impact of antibiotic use on antibiotic-resistant Escherichia coli using population-based data from a large hospital and its surrounding community[END_REF][START_REF] Batard | Tetracycline use in the community may promote decreased susceptibility to quinolones in Escherichia coli isolates[END_REF], and appear to be mainly related to the selection for gyr and par mutant strains of E. coli [START_REF] Redgrave | Fluoroquinolone resistance: mechanisms, impact on bacteria, and role in evolutionary success[END_REF]. Thus, the current study confirmed the usefulness of limiting the use of fluoroquinolones to decrease the rate of resistance to this class [START_REF] Sarma | Effects of fluoroquinolone restriction (from 2007 to 2012) on resistance in Enterobacteriaceae: interrupted time-series analysis[END_REF].

In accordance with previous studies [START_REF] Kaier | The impact of antimicrobial drug consumption and alcohol-based hand rub use on the emergence and spread of extendedspectrum -lactamase-producing strains: a time-series analysis[END_REF][START_REF] Kim | Impact of antibiotic usage on extended-spectrum β-lactamase producing Escherichia coli prevalence[END_REF], a correlation was found between fluoroquinolone consumption and resistance to cephalosporins and PBIs. The cross-selection pressure effect that is, at least in part, mediated by fluoroquinolone resistance genes carried on betalactamase plasmids, such as bla CTX-M , supports this finding [START_REF] Canton | The CTX-M β-lactamase pandemic[END_REF] . However, the effect of decreased cephalosporin consumption in limiting fluoroquinolone resistance in E. coli needs to be confirmed [START_REF] O'riordan | Antimicrobial use and antimicrobial resistance in Enterobacterales and Enterococcus faecium: a time series analysis[END_REF][START_REF] Guo | Antimicrobial resistance surveillance and prediction of Gram-negative bacteria based on antimicrobial consumption in a hospital setting: A 15-year retrospective study[END_REF]. On the contrary, the negative correlation of AAPBI consumption with resistance to fluoroquinolones and cephalosporins strengthens the current study result, highlighting the potential activity of tazobactam on extended-spectrum beta-lactamase (ESBL) [START_REF] O'riordan | Antimicrobial use and antimicrobial resistance in Enterobacterales and Enterococcus faecium: a time series analysis[END_REF][START_REF] Bantar | Replacement of broad-spectrum cephalosporins by piperacillin-tazobactam: Impact on sustained high rates of bacterial resistance[END_REF]. Similarly, the negative correlation between cephalosporin resistance and consumption of fourth-generation cephalosporin is supported by its activity on AmpC betalactamase [START_REF] Kaki | Impact of antimicrobial stewardship in critical care: a systematic review[END_REF].

Overall, the impact of beta-lactam consumption on AMR seemed limited in the current study. These results are in agreement with the literature, in which few studies have reported correlations between the consumption of beta-lactams and the acquisition of resistance to this class in E. coli [START_REF] Guo | Antimicrobial resistance surveillance and prediction of Gram-negative bacteria based on antimicrobial consumption in a hospital setting: A 15-year retrospective study[END_REF][START_REF] Hsueh | Relationships between antimicrobial use and antimicrobial resistance in Gram-negative bacteria causing nosocomial infections from 1991-2003 at a university hospital in Taiwan[END_REF]. Thus, beta-lactams seem to generate less resistance in E. coli than fluoroquinolones, which seem to have the greatest influence on the emergence of AMR.

To the best of our knowledge, this study was the first to assess the potential of DR models as AMR forecasting tools. Herein are re-ported the results of an external validation, showing an error percentage for AMR prediction ranging from 8% to 34%. Thus, within the framework of ASP development, the current study results underlined the need to consider the uncertainty of AMR estimates when DR models are used to plan restriction strategies targeting the AMC with the greatest impact.

This study has some limitations. First, its single-centre retrospective design may limit the generalisability of the results. Second, AMC was low in our hospital [11 , 22 , 33], below the reported threshold inducing AMR, which may have led to misestimating some correlations [START_REF] Lopez-Lozano | THRESHOLDS study group A nonlinear time-series analysis approach to identify thresholds in associations between population antibiotic use and rates of resistance[END_REF] . Third, the current study model did not include some mechanisms of AMR acquisition (colonisation pressure and out-hospital AMC), nor individual factors (co-morbidities, antimicrobial dose and duration) [START_REF] Turnidge | Antibiotic use and resistance-proving the obvious[END_REF], and it focused only on E. coli, which could have underestimated the impact of AMC on global AMR [START_REF] Mclaughlin | Correlations of antibiotic use and carbapenem resistance in Enterobacteriaceae[END_REF] . However, DR models are reputed to be more informative than other methodologies [START_REF] López-Lozano | Modelling and forecasting antimicrobial resistance and its dynamic relationship to antimicrobial use: a time series analysis[END_REF] , to assess the link between AMC and AMR, explaining the fraction of AMR related to AMC [START_REF] Turnidge | Antibiotic use and resistance-proving the obvious[END_REF] and accounting for autocorrelation of dependent resistance rates series and lag times. Finally, the COVID-19 pandemic restricted the avail-ability of AMC data for 2020 and 2021, and therefore limited the external validation of the model [START_REF] Lai | Increased antimicrobial resistance during the COVID-19 pandemic[END_REF] . Although the predictive abil-ities were suboptimal, the results highlighted the antibiotics most at risk for AMR, prompting physicians to limit their use in future ASPs.

Conclusions

Over a six-year period in a French teaching hospital, the in-crease in AMC was mainly driven by the increase in consumption of extended-spectrum penicillin (namely amoxicillin) and AAPBI (namely piperacillin-tazobactam), whereas the use of fluoroquinolones dramatically decreased. Among E. coli strains isolated during this period, there was a decrease in resistance to fluoroquinolone and cephalosporin, which was correlated with the de-crease in fluoroquinolone use and the increase in AAPBI use. Resistance to penicillin remained stable and was correlated with penicillin consumption.

The current study results underlined the complex relationship between AMC and AMR, indicating that it involves multiple mechanisms and depends on antimicrobial classes, and highlighted the burden of fluoroquinolone use. The study also showed that DR models are potentially helpful in ASP development. Nevertheless, the margin of error of the models should prompt clinicians to interpret the forecasted results with caution. Further studies focusing on other bacteria, in one or more hospitals, are mandatory to con-firm and extend the current study results.

Acknowledgements

We thank Sarah Kabani for her editing assistance. We also thank the team of the microbiology laboratory for their contribution to this work.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Competing Interests

None.

Ethical Approval

This study was conducted according to the guidelines of the Declaration of Helsinki. The Institutional Review Board of Nimes University Hospital approved the study (No. 210307).

Sequence Information

Not applicable.