Sijia, Maxime, Loris Aurélien

Nils

Aradana

Rémy

Xavier

Julie Paul

Grégoire

Keywords:

Un stage ? Mais vous ne voulez pas un sujet de thèse plutôt ? J'ai un sujet sur l'atterrissage de fusées réutilisables qui pourrait vous plaire." C'est sur cette phrase de Nicolas Petit qu'est née l'idée de cette thèse, en 2016, lors d'une discussion à propos d'un sujet qui n'avait pas grand chose à voir. A l'époque, j'étais un jeune étudiant travaillant sur l'atténuation des vibrations pour les machines à laver, sous la supervision de Florent Di Meglio.

Je tiens à remercier Nicolas de m'avoir proposé ce sujet et encadré depuis tout ce temps. Merci d'avoir toujours répondu présent, même pendant cette longue période de crise sanitaire. Ce n'est pas seulement avec un chercheur que j'ai pu travailler, mais avec un passionné d'aérospatial. Je tiens aussi à remercier Eric Bourgeois de m'avoir encadré avec Nicolas.

Je tiens à remercier toute l'équipe du Centre Automatique et Systèmes pour ces nombreuses années de collaboration. Tous les permanents m'ont aidé à un moment ou à un autre de ma thèse. Merci à Delphine pour m'avoir fait confiance pour les cours d'optimisation, me laissant une grande liberté dans le choix de mes sujets, et d'avoir été toujours disponible quelle que soit la question. Merci à Florent, pour son soutien sans faille depuis toujours. Merci à Philippe, pour toutes ces discussions autour de la machine à café, et les nombreux coups de main techniques. Merci à Laurent, pour ses relectures précises, et ses conseils techniques avancés.

Merci infiniment à Dilshad, d'avoir été un co-doctorant et ami pendant 3 ans, de m'avoir initié au puit sans fond qu'est tikz, de m'avoir fait découvrir des sujets techniques comme littéraires, et d'avoir refait le monde autour d'un nombre peu raisonnable de cafés. Merci à tous les autres doctorants et anciens doctorants du laboratoire,

x = Rocket states x = (h, v h , z, v z , m) ⊤ ∈ R 5 . u = Guidance controls u = (q r , α) ⊤ ∈ R 2 .

Nomenclature for the 3D model α = Incidence (In 3D) α is unsigned: a nor = Normal acceleration a nor is signed.

0 • ≤ α < 90 • . Ψ = Yaw
x = Rocket states x = (z, y, h, v z , v y , v h , m, q r ) ⊤ ∈ R 8 . u = Guidance controls u = (q, α z , α y ) ⊤ ∈ R 3 . Chapter 1

Nomenclature for the optimization methods

ξ = Input variable Size N ξ , s.t. ξ = (∆x 0 ⊤ , ∆η ⊤ , ∆u init ⊤ ) ⊤ . x =

Introduction

Résumé

Ce chapitre introduit les différents enjeux liés au calcul embarqué de trajectoires d'atterrissage d'urgence pour des lanceurs réutilisables. Après avoir présenté le contexte dans lequel se déroule le développement de ces lanceurs, les différentes familles de méthodes de guidage disponibles dans la littérature sont rappelées. Ceci permet de présenter le problème au coeur de cette thèse, et les outils nécessaires pour le résoudre : l'ordre co-lexicographique (ordre d'urgence), les problèmes (linéaires) de négotiation et le problème (quadratique) de rafinement. Enfin, le plan du manuscrit est annoncé, et un résumé schématique est proposé.

Powered Descent Guidance for reusable launchers

Reusable launchers are autonomous vehicles tailored for complex missions consisting in a succession of distinct phases. Their development is driven by the need for a fast and cost-efficient solutions to put payloads into orbit. Reusability imposes that return to Earth of the launcher's first stage must be accomplished safely and reliably. Powered Descent Guidance (PDG) refers here to the action of making such a flying vehicle land autonomously on a horizontal surface using rocket-like engines for maneuvers and deceleration. PDG has emerged as a paramount topic during the 1960's space race, the prime example of PDG being the lunar module (LM) landing of the Apollo missions, which was successfully performed six times in a row on the moon. LM landing used surprisingly straightforward path planning methods [START_REF] Klumpp | Apollo Lunar Descent Guidance[END_REF][START_REF] Steinfeldt | Guidance, Navigation, and Control System Performance Trades for Mars Pinpoint Landing[END_REF], relying on non-optimized analytic calculations in an atmosphere-free environment.

Though the space race slowed down drastically in the 1970's, the need for PDG on other planets re-emerged with automated exploration missions towards Mars, the Moon and even some large asteroids. These missions require a high level of autonomy, for they are conducted extremely far from Earth, and thus cannot rely on near-instantaneous communications used in teleoperation control. Accuracy is also a key factor, especially for Mars missions, which aims at landing in the vicinity of rich geological features where landing areas are scarce and narrow [START_REF] Blackmore | Autonomous Precision Landing of Space Rockets[END_REF]. Several missions performed a series of increasingly accurate autonomous powered landings, peaking with the impressive performance of the Perseverance Rover which landed . Interestingly, all of these missions operated in the absence of atmosphere or under negligible atmospheric density.

For PDG on Earth's surface, where aerodynamic forces are not negligible, progress is more recent. At the end of the 1990's, the DC-X project conducted collaboratively between McDonnell Douglas and NASA started to explore the possibilities for a reusable rocket that could land back on its "feet", using the same engines for take-off and landing. As shown in Figure 1.1-(Left), this atypical tetrahedron shaped vehicle achieved several short flights, but never reached a significant altitude. It was necessary to wait until the early 2010's to witness the first successful powered landings, achieved by Space X rocket Grasshopper. In parallel, some smaller prototypes were developed and successfully tested, such as the Xombie by Masten Space Systems and JPL 1 . Today, several projects are under developments and have reached various Technology Readiness Levels, such as New Shepard and New Glenn (Blue Origins), FROG (CNES, [START_REF] Rmili | FROG, a Rocket for GNC demonstrations: Firsts flights attempts of the FROG turbojet version and preparation of the future mono-propellant rocket engine[END_REF]), Callisto (CNES/DLR/JAXA) and Themis (Ariane-Group/CNES/ONERA) or even Shenlan Nebula (Deep Blue Aerospace).

In this thesis, we are interested in tossback vehicles and, more precisely, in reusable Two-Stage-To-Orbit (TSTO) rockets, equipped with a single gimbaled engine, whose trajectory slenderness ratio2 is medium to high. As shown in Figure 1.2, this type of launcher typically goes through six different flight phases: take-off, boost-back, ballistic, re-entry-burn, re-entry glide and final-burn. Our work focuses on the last flight phase, the final burn, which starts at a few kilometers of altitude, lasts for 10 seconds to a couple of minutes during which the engine is always turned on and (at least partially) controllable.

To land properly at a desired site, the typical Guidance and Control (G&C) strategy considers two objectives: i) the guidance problem i.e. determination of a high-level guidance trajectory , ii) the control problem i.e. tracking of this trajectory using a low-level controller. This thesis addresses the guidance problem only under the assumption of full knowledge of the state of the system, making low-level control and other estimation topics out of scope.

Providing guidance for the final burn is a challenging and necessary task. Indeed, solving the PDG in-flight must be achieved using limited time and computational power. Only single-CPU methods are studied in the thesis, as opposed to methods requiring the use of GPUs. Further, PDG has to deal with the accumulation of the tracking errors -i.e. the distance between the actual rocket states and the expected reference trajectory -during the flight phases prior to the final burn, which have to be shrunk to zero during this last phase [START_REF] Blackmore | Autonomous Precision Landing of Space Rockets[END_REF]. This must be done under strong disturbances (e.g. wind speed or changes of atmosphere density profile among others) and under multiple constraints: the incidence -or angle-of-attack -must be limited to avoid hazardous regions of the aerodynamic flight domain, the engine flow is mechanically bounded and its internal dynamics is not instantaneous, the low-level controllers actuating the rocket have limited capabilities which requires to bound the normal acceleration, the landing site is small, the available mass is limited, etc. Gathering all these detrimental effects, it is likely that the rocket starts its final burn too far from its reference trajectory to get a landing trajectory strictly meeting all these requirements. However, at the expense of sacrificing some requirements and/or loosening well-chosen constraints, it is possible to successfully land and maintain the reusability of the vehicle.

A key concept to preserve the launcher's reusability is to "maximize the launcher's integrity during landing". To illustrate it, consider the following example, pictured in Figure 1.3. The reference trajectory A is shown in plain black. In a typical nominal scenario, the final burn would start reasonably close to the reference trajectory such that a proper guidance algorithm would provide the rocket B with a trajectory that reaches the desired landing at null speed while satisfying all the constraints listed earlier. However, if the rocket starts its descent farther away from the reference trajectory -e.g. with a large lateral displacement -then, at some point, it is impossible to land properly and to satisfy the whole system of constraints, which are in fact incompatible. Allowing sharper turns, which can be achieved by relaxing the incidence bound (or the normal acceleration bound) may make the system of constraints compatible again, at the cost of a trajectory alteration. Indeed, this relaxation enables the rocket C to land on the desired landing site. If we push this reasoning further, the rocket could start its final burn far enough from the reference trajectory so that even relaxing the incidence bound would not be enough. If the area neighboring the desired landing site is uninhabited and flat enough, one can allow the landing location to be relaxed too. This is what happens in the case of rocket D.

This first point highlighted by this example is that there is a list of quantifiable elements that can help to relax the constraints (to some extent). These will be referred to as negotiable parameters in this thesis. Maximizing the launcher's integrity is choosing the values of these parameters in an optimal way. The negotiable parameters can refer to many different factors, such as the incidence bound, the normal acceleration bound, the landing site location or even the final vertical speed. The set of negotiable parameters can differ depending on the mission requirements, the rocket structural capabilities or the landing site surroundings for instance. In the thesis we develop a generic methodology able to handle different sets of landing parameters, and not only a single particular set.

Further, the above-mentioned example illustrates why these parameters are not equally critical. Indeed, relaxing the incidence bound of a couple degrees will have a much lower impact on the rocket structure than relaxing the maximum final vertical speed by a few meters per second. A classic way to handle this relative importance while solving the PDG problem is to penalize the constraint incompatibility [START_REF] Chinneck | Feasibility and Infeasibility in Optimization: Algorithms and Computational Methods[END_REF]. Computationally, it means that: i) the negotiable parameters are added to the decision variables of the optimization problem describing the descent trajectory, and that ii) the use of these parameter is penalized in the cost of the optimization problem. In this case, the relative importance of the parameters is partially ensured by the difference in the weights associated to each parameter. However, this is a heuristic, and it does not guarantee the relative importance of the parameter to be exactly satisfied. Instead, we propose to mathematically define the relative importance of the negotiable parameters via the introduction of a strict hierarchy between these parameters. In the thesis, this hierarchy is represented by a specific order relation over the set of negotiable parameters, which we call the emergency order below. A challenge for the guidance methodology to be developed is to compute the smallest constraint alteration that recovers feasibility in the sense of the latter order.

Following the discussion above, our main objective in this thesis is formulated as follows: design a method to update the guidance trajectory at the beginning of the final-burn, while maximizing the launcher's integrity, in a computationally efficient way. To address this objective, the thesis proposes a guidance method that performs online trajectory optimization for the final burn, while minimizing the use of negotiable parameters according to a certain hierarchy making sure that their relative importance is strictly respected.

Mathematical programming for online PDG 1.2.1 Current state-of-the-art

The Powered Descent Guidance (PDG) problem is a worked applied mathematics problem, most often tackled from the Optimal Control perspective.

The decision variables belong to infinite-dimensional function spaces (the control law u and the state x), and the scalar time-of-flight (t f ) over which these functions are defined. The goal in such Optimal Control Problems (OCPs) is to find the triplet (x, u, t f ) that minimizes a certain criterion, under multiple constraints.

Several performance indexes (or cost functions) have been considered in the literature: minimum-fuel [START_REF] Açikmeşe | G-FOLD: A Real-Time Implementable Fuel Optimal Large Divert Guidance Algorithm for Planetary Pinpoint Landing[END_REF][START_REF] Carson | Lossless convexification of Powered-Descent Guidance with non-convex thrust bound and pointing constraints[END_REF][START_REF] Leparoux | Structure of optimal control for planetary landing with control and state constraints[END_REF][START_REF] Meditch | On the problem of optimal thrust programming for a lunar soft landing[END_REF][START_REF] Szmuk | Successive Convexification for Fuel-Optimal Powered Landing with Aerodynamic Drag and Non-Convex Constraints[END_REF][START_REF] Wang | Optimal Rocket Landing Guidance Using Convex Optimization and Model Predictive Control[END_REF], minimum-acceleration [START_REF] Steinfeldt | Guidance, Navigation, and Control System Performance Trades for Mars Pinpoint Landing[END_REF], minimum time-of-flight [START_REF] Meditch | On the problem of optimal thrust programming for a lunar soft landing[END_REF][START_REF] Steinfeldt | Guidance, Navigation, and Control System Performance Trades for Mars Pinpoint Landing[END_REF], minimum error-to-reference [START_REF] Song | Survey of autonomous guidance methods for powered planetary landing[END_REF], minimum-landing-error [START_REF] Blackmore | Minimum-Landing-Error Powered-Descent Guidance for Mars Landing Using Convex Optimization[END_REF], or any weighted combinations of these [START_REF] Souza | An optimal guidance law for planetary landing[END_REF][START_REF] Sagliano | Guidance and Control Strategy for the CALLISTO Flight Experiment[END_REF].

The constraints depend on the rocket design and the mission requirements. They can take the form of inequalities (e.g. engine flow upper and lower bounds) as well as equalities (e.g. final horizontal and vertical speeds). As discussed earlier, these constraints can be quantified by several parameters. Among these, the ones that can eventually be modified are called the negotiable parameters. For instance, the incidence bound is a safety requirement regarding mechanical and flight quality aspects, whose value may be negotiated if necessary. On the contrary, the engine flow bound is a typical example of non-negotiable parameters, since it is a physical limitation. Mathematically, these cases are treated as follows: the constraints can be tuned using a vector parameter denoted p, which is defined such that p = 0 when the nominal requirements are met.

The PDG problem has several inputs: the initial state of the rocket -i.e. its initial position, speed, orientation and mass -and other external parameters, such as the wind speed for example. These cover all the above-mentioned sources of disturbances. Mathematically, the inputs are conveyed by a parameter ξ.

From a general perspective, the PDG problem, seeking a guidance trajectory for the final burn, is written as an OCP such that OCP (ξ, p) :=

                         min x,u,t f J (x, u, t f , ξ) (Cost), s.t. ẋ = f (x, u) on [0, t f ] (Rocket dynamics),
x(0) = x 0 + ξ (Initial conditions), Ψ(x(t f ), p) = 0 on [0, t f ] (Final state), C(x, u, p) ≤ 0 on [0, t f ] (Other constraints).

where J is the performance index (cost), f is the right-hand side of the ordinary differential equations, Ψ represents the terminal condition and C conveys the var-ious constraints. The special case OCP (ξ, 0) is the nominal PDG problem, and OCP (0, 0) corresponds to the reference trajectory. This problem having an infinite dimensional decision variable and an infinite number of constraints, it is discretized for numerical resolution (for some given value of ξ and p). Two well-known and dual discretization approaches exist (see e.g. [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF][START_REF] Hull | Optimal control theory for applications[END_REF][START_REF] Rao | Trajectory Optimization: A Survey[END_REF][START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF]). On one hand, for indirect methods, stationary conditions are derived from OCP (ξ, p) in the form of Ordinary Differential Equations (ODEs) which are then discretized and solved. The celebrated Pontryagin Maximum Principle [START_REF] Hartl | A Survey of the Maximum Principles for Optimal Control Problems with State Constraints[END_REF] is the fundamental tool to form the stationary conditions. On the other hand, for direct methods, OCP (ξ, p) is first discretized into the sub-problem DSC (ξ, p), using a finite dimensional variable z representing the unknowns (x, u, t f ) of the initial problem, where DSC (ξ, p) is an optimization problem of the form DSC (ξ, p) :=

           min z J(z, ξ)
s.t. h(z, ξ, p) ≤ 0, g(z, ξ, p) = 0, which is then solved using Non-Linear Programming (NLP).

Though the indirect methods are known to be more accurate, they are also very sensitive 3 and usually are not used in real-time. Currently, for online applications, PDG problems are solved almost exclusively using direct methods. We follow this trend and will focus of the parametric problem DSC (ξ, p).

Still, the main challenges for direct methods are the choice of the decision variable z, the design of the associated functions J, h and g, and the selection of the proper numerical method used to solve DSC (ξ, p). These decisions directly determine whether DSC (ξ, p) can be used for online applications with sufficient accuracy.

A broad spectrum of direct methods have been used to solve OCP (ξ, p) for aerospace applications [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF]. Direct trajectory optimization methods relying on NLP have been explored since the 1980's [START_REF] Hargraves | Direct trajectory optimization using nonlinear programming and collocation[END_REF]. Since then, many methods have emerged. Pseudospectral methods -consisting in discretizing the states and the controls of OCP (ξ, p) at well-chosen time points using polynomial approximations [START_REF] Fahroo | Direct Trajectory Optimization by a Chebyshev Pseudospectral Method[END_REF] to transform OCP (ξ, p) into a NLP -have been investigated from various perspectives [START_REF] Ross | A review of pseudospectral optimal control: From theory to flight[END_REF][START_REF] Sagliano | Onboard Guidance for Reusable Rockets: Aerodynamic Descent and Powered Landing[END_REF][START_REF] Sostaric | Powered Descent Guidance Methods For The Moon and Mars[END_REF][START_REF] Sostaric | Powered Descent Trajectory Guidance and Some Considerations for Human Lunar Landing[END_REF][START_REF] Wang | A Pseudospectral-Convex Optimization Algorithm for Rocket Landing Guidance[END_REF]. Sensitivity-based parametric optimization has been used for Mars atmospheric entry problems [START_REF] Seelbinder | On-board Trajectory Computation for Mars Atmospheric Entry Based on Parametric Sensitivity Analysis of Optimal Control Problems[END_REF]. Recent work also focused on applying learning based methods to PDG [START_REF] Furfaro | Waypoint-Based Generalized Zem/Zev Feedback Guidance For Planetary Landing Via A Reinforcement Learning Approach[END_REF]. A noteworthy state-of-the-art direct method aiming at solving OCP (ξ, p) for rocket landing problems is Successive Convexification [START_REF] Açikmeşe | Convex Programming Approach to Powered Descent Guidance for Mars Landing[END_REF]. This method solves OCP (ξ, p) by successively solving linearized and discretized versions of it [START_REF] Mao | Successive convexification of non-convex optimal control problems and its convergence properties[END_REF][START_REF] Reynolds | Optimal Planar Powered Descent with Independent Thrust and Torque[END_REF][START_REF] Szmuk | Successive Convexification for Real-Time 6-DoF Powered Descent Guidance with State-Triggered Constraints[END_REF]. In other words, this method is a variant of Successive Convex Programming (SCP) methods, tailored and optimized for the requirements of PDG. This method is presented in great details in [START_REF] Malyuta | Convex Optimization for Trajectory Generation[END_REF]. Theoretical guarantees were established by B. Açikmeşe and his fellow researchers. Informally speaking, the main statement can be rephrased as follows: "if the method converges and if some intermediate slack variables are sufficiently penalized and tend to zero, then the solution is the optimal landing method in the sense of the Karush-Kuhn-Tucker conditions and the constraints are satisfied" [START_REF] Malyuta | Convex Optimization for Trajectory Generation[END_REF]Thm. 8]. Superlinear convergence rates have been established under mild assumptions [START_REF] Mao | Successive Convexification: A Superlinearly Convergent Algorithm for Non-convex Optimal Control Problems[END_REF]. Specifically designed interior point solvers have also been implemented [START_REF] Dueri | Customized Real-Time Interior-Point Methods for Onboard Powered-Descent Guidance[END_REF]. These results are rather powerful and the method proved to work correctly even on real vehicles [START_REF] Açikmeşe | G-FOLD: A Real-Time Implementable Fuel Optimal Large Divert Guidance Algorithm for Planetary Pinpoint Landing[END_REF]. It is often used along with Lossless Convexification, which aims at performing an exact convex relaxation of the thrust magnitude constraints [START_REF] Blackmore | Lossless convexification of control constraints for a class of nonlinear optimal control problems[END_REF][START_REF] Carson | Lossless convexification of Powered-Descent Guidance with non-convex thrust bound and pointing constraints[END_REF][START_REF] Malyuta | Convex Optimization for Trajectory Generation[END_REF]. The above-mentioned list of direct methods for PDG is not exhaustive. For further details on the available methods, see the recent survey by Song et al. [START_REF] Song | Survey of autonomous guidance methods for powered planetary landing[END_REF].

All the above-mentioned methods tackle different versions of the PDG problem, but face the same bottleneck when it comes to landing on Earth's surface: the presence of the aerodynamic forces. Without atmosphere, i.e. for Moon and Mars landing 4 , the problem is called planetary landing, and one can safely say that the above-mentioned literature is mature enough to provide efficient and proven PDG methods for real-time usage. A noteworthy example is the full characterization of the reachable set for Mars landing [START_REF] Eren | Constrained Reachability and Controllability Sets for Planetary Precision Landing via Convex Optimization[END_REF], obtained by combining convex analysis and Successive Convexification. Though a few papers have tackled PDG in the past decade, there are still many open questions regarding the best way to perform online PDG in the presence of non-negligible aerodynamic effects. In this thesis, this is the case we consider.

Hierarchy and the emergency problem

In this thesis we address cases when ξ is such that there exist no trajectories (x, u, t f ) that satisfy all the constraints of OCP (ξ, p). This situation, referred to as the emergency problem, is also sometimes called the abort planning problem in the air and space field [START_REF] Calise | Generation of Launch Vehicle Abort Trajectories Using a Hybrid Optimization Method[END_REF][START_REF] Hanson | Ascent, transition, entry, and abort guidance algorithm design for the X-33 vehicle[END_REF][START_REF] Lampazzi | Intact Ascent Aborts Workbook 21002[END_REF][START_REF] Lu | Abort Guidance during Powered Descent for Crewed Lunar Missions[END_REF][START_REF] Shapiro | Survivability of Emergency Escape from a Simulated Shuttle Entry Trajectory[END_REF][START_REF] Vana | Any-Time Trajectory Planning for Safe Emergency Landing[END_REF]. To the best of our knowledge, a single paper has tackled this topic for planetary PDG, for Mars landings [START_REF] Blackmore | Minimum-Landing-Error Powered-Descent Guidance for Mars Landing Using Convex Optimization[END_REF]. Using Lossless and Successive Convexification, the latter paper exposes a method that first minimizes the landing error -i.e. the minimum distance between the feasible landing sites and the targeted one -and then computes the fuel-optimal trajectory among the ones having a minimum landing-error. Our goal is to design a more general methodology to perform emergency guidance by dealing with more than a single parameter to relax.

From a modeling perspective, in this thesis, the above-mentioned vector of negotiable parameters p is decomposed into R vector parameters of possibly different dimensions, ranked with respect to their relative importance, such that p = p (1) , . . . , p (R) = Least critical, . . . , Most critical .

For instance, one can consider that p (1) = ∆α max (the incidence bound) is less critical than p (2) = ∆v f h (the final vertical speed). The elements p (j) can be vectors. For instance, p (j) = (∆z f , ∆y f ) is the vector denoting the landing site location on a map.

The hierarchy among the parameters is defined using a variant of the lexicographic order: the negotiable parameters are compared using the 1-norm of their sub-parameters, by comparing their most critical sub-parameters first. In practice, a vector p a is said to be larger -or more negotiated -than another vector p b , which is denoted

p a ⪰ e p b if and only if ∥p (R) a ∥ 1 > ∥p (R) b ∥ 1 or ∥p (R) a ∥ 1 = ∥p (R) b ∥ 1 and ∥p (R-1) a ∥ 1 > ∥p (R-1) b ∥ 1 , or . . . or ∥p (R) a ∥ 1 = ∥p (R)
b ∥ 1 and . . . and ∥p (1) a ∥ 1 > ∥p

(1) b ∥ 1 , or ∥p (R) a ∥ 1 = ∥p (R)
b ∥ 1 and . . . and ∥p (1) a ∥ 1 = ∥p

(1) b ∥ 1 .
(1.1)

In the context of this thesis, the order ⪰ e is called the extended colexicographic order or simply the emergency order.

Maximizing the launcher's integrity by sacrificing the parameter p according to the hierarchy of importance between the sub-parameters translates into finding the smallest p in the sense of ⪰ e such that there is a feasible trajectory (x, u, t f ) satisfying the constraints of OCP (ξ, p). There are no generic methods yet available that recovers feasibility in OCP (ξ, p) while enforcing such a hierarchy. This is the subject and main contribution of this thesis.

Proposed contribution: Hierarchical Emergency Guidance Optimization

The main contribution of this thesis is an online method that computes the optimal trajectory and the optimal values of the negotiable parameters, in the sense of the extended colexicographic order, using only Linear and Quadratic programming techniques. In the manuscript, this method is demonstrated on two rocket models of increasing complexity: from a planar model with non-trivial aerodynamic effects, to a richer three-dimensional model with non-negligible engine transients. This contribution is obtained in two steps: first a nominal guidance method computing the landing trajectory z for a given p is presented, then an emergency guidance method aiming at computing the best value for p is studied.

Fast PDG for fixed value of constraint parameters : a sensitivity-based approach

To perform nominal guidance, i.e. to compute z for given values of ξ and p, we use a sensitivity-based approach. In this approach, OCP (ξ, p) is described as an optimal correction problem PDG (ξ, p), defined w.r.t. a reference trajectory (x, ū, tf ) (e.g. rocket A in Figure 1.3). The latter trajectory is mission-specific, and its design is out-of-scope for this manuscript. Using handy notations 5 , PDG (ξ, p) is sketched below as an infinite dimensional optimization problem of the form

PDG (ξ, p) :=                          min δu,∆t f J (δu, ∆t f ) s.t. ẋ(t) = f (x(t), ū(t) + δu(t)), x(0) = x0 + ξ, Ψ(x( tf + ∆t f ), p) = 0, C(x(t), ū(t) + δu(t), p) ≤ 0.
The decision variable δu, the ODE of the dynamics and the constraints are discretized using a finite-dimensional variable z, which conveys a parametric description µ of the control change w.r.t. the reference control ū and the implicit time-of-flight change ∆t f . Through this discretization procedure, the state x is expressed as a function of ξ and z using the flow of the ODE defined by f . The problem PDG (ξ, p) is approximated by its discrete non-linear version NLP (ξ, p), which writes

NLP (ξ, p) :=            min z J(z, ξ) h(z, ξ) ≤ H p p, g(z, ξ) = B p p.
where it is stressed that the negotiable parameter p has a linear influence on the constraint right-hand side of the guidance problem tackled in this thesis.

Then, sensitivity analysis is used to approximate the solution of NLP (ξ, p). Revisiting results from the literature, we introduce a Quadratic Program (QP) which writes

QP (ξ, p) :=            min z 1 2 z ⊤ P z + q ⊤ z s.t. Gz ≤ h 0 + H ξ ξ + H p p Az = b 0 + B ξ ξ + B p p
and whose solution is shown to be a local approximation to the solution of NLP (ξ, p), under mild assumptions. Interestingly, Strict Complementary Slackness, often assumed in the literature, is not needed here. Approximating NLP (ξ, p) by QP (ξ, p) is shown to be relevant and usable even for large -non-local -values of ξ, when used with the two above-mentioned rocket models. The problem QP (ξ, p) contains a convenient linear description of the constraints which makes the emergency guidance method described thereafter tractable.

Problems such as QP (ξ, p) can be solved in an online/offline fashion, where the defining constant matrices are computed before the flight for the prescribed reference trajectory, and where the QP itself is solved for given ξ and p on-board using offthe-shelf QP solvers. Numerical resolution can be achieved with confidence as QP solvers are now considered mature and reliable technology [START_REF] Boyd | Convex Optimization[END_REF][START_REF] Mattingley | CVXGEN: a code generator for embedded convex optimization[END_REF][START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF].

Computing the best constraint alteration achieving feasibility

When ξ is such that the nominal guidance problem QP (ξ, 0) has no solutions, we propose a method computing a value of p such that QP (ξ, p) is feasible.

Finding the smallest p in the sense of ⪰ e guaranteeing the existence of a trajectory is achieved by solving the following R different negotiation problems. The idea is to iteratively minimize each sub-parameter p (j) under the 1-norm (consistently with the mathematical definition of ⪰ e ), and memorize the associated optimal value P * j . For each negotiation problem, the constraints are those of NLP (ξ, p), plus new constraints that ensure that the optimal values P * j+1 , . . ., P * R of the previous negotiation problems are reached. In details, the j th negotiation problem writes

P * j ←-min z, p ∥p (j) ∥ 1 (1.2a) s.t. Gz ≤ h 0 + H ξ ξ + H p p (1.2b) Az = b 0 + B ξ ξ + B p p (1.2c) ∥p (i) ∥ 1 = P * i , i = j + 1, . . . , R.
(1.2d)

The HEGO algorithm: main contribution of the thesis

Once this problem has been solved for each j, starting from R down to 1, we minimize z under the original performance index J, such that

z * ←-arg min z, p 1 2 z ⊤ P z + q ⊤ z (1.3a) s.t. Gz ≤ h 0 + H ξ ξ + H p p (1.3b) Az = b 0 + B ξ ξ + B p p (1.3c) ∥p (i) ∥ 1 = P * i , i = 1, . . . , R. (1.3d)
This latter problem is the refine problem. It only returns the value of z. The value of p is not necessarily unique. The "negotiation and refine" problems are gathered into the HEGO algorithm, which is a nominal and emergency guidance algorithm, described using pseudo-code below. Overall, HEGO is a numerical method that provides nominal and emergency guidance, relying on Linear and Quadratic programming solvers only.

Manuscript outline

A high-level summary of the approach with references to the associated chapters is presented in Figure 1.4.

Chapter 2 introduces the dynamic model of the class of reusable launchers under study, with two levels of complexity (2D and 3D rocket models). Among others, the non-trivial aerodynamic model is detailed, and special care is taken to define the 3D model. Then, after discussing the various mission constraints, the PDG problem is formulated as an OCP with respect to a reference trajectory.

Chapter 3 takes a side turn and tackles the problem of optimal thrust programming for the special case of a purely vertical atmospheric flight. By applying Pontryagin Maximum Principle, necessary and sufficient conditions are derived, showing that the optimal thrust program is min-max for the class of launchers that we study. The by-product of this result is a full characterization the reachable set of the rocket for the vertical landing problem. For a first reading, this chapter can be skipped without loss of continuity.

Chapter 4 details the method that performs nominal trajectory planning, sketched in Section 1.3.1. Using a finite-dimension decision variable to describe the trajectory of the PDG, a NLP is derived. Its optimal solution is then approximated using parametric sensitivity analysis, and solved using a single QP. The mathematical formulation of this QP is instrumental in the rest of the thesis.

Chapter 5 presents the core topic of this thesis. After discussing the available negotiation parameters, the Algorithm HEGO is introduced in the LP/QP framework. Its behavior is explained on a detailed toy example. Then, theoretical guarantees are proposed and proved. Among others, the Lipschitz-continuity of the optimal solution z * is established, guaranteeing the absence of "jumps" in the solutions, a desirable property in practice. Also, high-level comments on the emergency guidance method are proposed, to distinguish what is generalizable from what comes from the underlying nominal guidance problem. Finally, several examples are presented to illustrate the various modeling possibilities offered by HEGO, and to visualize the quality of the results.

Chapter 6 offers a quantitative assessment of the performances of HEGO. The inputs of the algorithm are dispersed over wide uncertainty intervals, and the results are analyzed pair-wise. Also, a comparison with the vertical landing problem from Chapter 3 is presented.

A concluding chapter discusses a few topics that have not been detailed in the previous chapters, it also presents possible future research directions, and conclusions.

Hierarchical Emergency

Guidance Optimization (Chapter 5)

P * j ←-min z, p ∥p (j) ∥ 1 s.t.: Gz ≤ h 0 + H ξ ξ + H p p, Az = b 0 + B ξ ξ + B p p, ∥p (i) ∥ 1 = P * i , i = j + 1, . . . , R.
Negotiation problem (loop for j = R, . . . , 1)

z * ←-arg min z, p 1 2 z ⊤ P z + ξ ⊤ Qz s.t.: Gz ≤ h 0 + H ξ ξ + H p p, Az = b 0 + B ξ ξ + B p p, ∥p (i) ∥ 1 = P * i , i = 1, . . . , R.

Refine problem

Discrete approximation of PDG (ξ, p) (from Chapter 2).

Negotiable parameters (e.g. incidence bound, landing site, etc.) p = p (1) , . . . , p (R) least -→ most critical Chapter 2

Dynamic models and the PDG problem

Résumé

Ce chapitre contient une description mathématique du problème de guidage pour l'atterrissage (PDG). Il introduit les modèles dynamiques dérivant la fusée et les contraintes. Tout d'abord, les choix généraux de modélisation sont discutés, ce qui permet de souligner le rôle de l'atmosphère, une question rarement abordée dans le cadre du guidage pour l'atterrissage. Deux modèles de fusée, avec différents niveaux de complexité, sont présentés. Un modèle de fusée dans le plan (2D) et un modèle de fusée tridimensionnel (3D) sont construits. Le modèle de fusée 2D suppose que la fusée reste toujours dans un seul plan. Il servira à illustrer les principes de guidage des chapitres suivants d'une manière beaucoup plus accessible que le modèle 3D. En revanche, le modèle 3D est utilisé dans les exemples avancés des Chapitres 4 et 5 et pour l'ensemble du Chapitre 6, ce dernier contenant une évaluation des performances numériques et vérifie l'applicabilité de la méthodologie proposée.

Un problème général de calcul de trajectoire est présenté à la fin de ce chapitre, sous la forme d'un OCP en temps final libre, défini par rapport à une trajectoire de référence. La résolution de ce problème sera l'objet principal du Chapitre 4. L'importance relative des contraintes du PDG sera examinée plus loin dans le Chapitre 5.

This chapter contains a mathematical description of the PDG problem. It introduces the dynamic models governing the rocket and the constraints.

First, general modeling choices are discussed, which serve to stress the role of the atmosphere, a seldomly looked at issue in PDG. Two rocket models, with different levels of complexity are presented. A planar (2D) and a three-dimensional (3D) rocket model are constructed. The 2D rocket model assumes that the rocket always remains in a single plane. It will serve to illustrate the guidance principles of the next chapters in a much more accessible way than the 3D model. On the other hand, the 3D model is used in the advanced examples of Chapter 4 and 5 and for the whole 15 Chapter 6 containing the numerical results assessing the numerical performance and the applicability of the proposed methodology.

A general trajectory design problem is presented at the end of this chapter, as a free-final time OCP defined w.r.t. a reference trajectory. The resolution of this problem will be the main concern of Chapter 4. The relative importance of the PDG constraints will be discussed later in Chapter 5.

Atmospheric flight dynamics

Here are presented the features shared by both 2D and 3D models. Some notions, such as those regarding the environment and the aerodynamic model are also used in Chapter 3.

The typical flight phases of a tossback vehicle are illustrated in Figure 1.2. In this thesis, we are interested in the last part of the flight: the final burn until landing. It starts a few kilometers above the ground [START_REF] Brendel | Optimal guidance for toss back concepts of Reusable Launch Vehicles[END_REF].

In both 2D and 3D rocket models, which will be used in guidance algorithms, the rocket is assumed to be a punctual mass having an orientation (and thus having an aerodynamic incidence). This conceptual approach is presented in details below.

Earth and atmosphere model

Since we are only interested in the last flight phase, the Earth is considered locally flat, non-rotating, with a constant gravity field of magnitude g. The atmosphere is described via the pressure P a , the density ρ, the temperature and the speed of sound S SP are functions of the altitude. They are computed using linear interpolation of data samples.

Denoting by V r the norm of the relative speed of the rocket, we define M a := V r /S SP (h) the Mach number at a given altitude h.

The wind is assumed to be horizontal with a speed depending on the altitude only. It is assumed null at null altitude. In practice, the wind map is described by its value at three reference altitudes1 , for each direction:

• w z,0 , w z,1 and w z,2 for the first horizontal direction (2D and 3D models),

• w y,0 , w y,1 and w y,2 for the second horizontal direction (3D model only).

For these wind parameters, index 0 (resp. 1, 2) corresponds to an altitude of h = 2 km (resp. 5 km, 10 km).

Aerodynamic model

The aerodynamic model of a rocket moving in the direction of its thrust flame is notoriously hard to determine. Early works from 1966 started to describe the aerodynamic effect of an air jet pushing in front of a body in a supersonic flow.

Later works from the early 2000's from JAXA have completed these observations [START_REF] Nonaka | Vertical Landing Aerodynamics of Reusable Rocket Vehicle[END_REF]. They were followed and corroborated by recent experiments of the DLR [START_REF] Marwege | First Wind Tunnel Data of CALLISTO Reusable VTVL Launcher First Stage Demonstrator[END_REF]. The common findings of [START_REF] Nonaka | Vertical Landing Aerodynamics of Reusable Rocket Vehicle[END_REF] and [START_REF] Marwege | First Wind Tunnel Data of CALLISTO Reusable VTVL Launcher First Stage Demonstrator[END_REF] are that for a sufficiently strong jet flow, its wrapping around the rocket body lowers drastically the drag along the body axis, though orthogonal effects remain strong for non-zero incidences. This can be partially explained by the fact that the air jet creates an air cushion in front of the rocket that helps the boundary layer to stick to the fuselage from its lower end. A qualitative explanation of this observation is presented in Figure 2.1.

For the sake of this thesis, the aerodynamic effects of the air flow around the rocket is taken into account via:

1. A non-trivial lift coefficient C Lift , to account for the aerodynamic effect orthogonal to the drag and in the opposite direction to the relative speed, 2. A small-magnitude drag coefficient C Drag , to account for the aerodynamic effect in the direction of the rocket body, 3. Drag and lift coefficients C Drag and C Lift depending on two parameters: the rocket incidence α and the Mach number M a ,

4. An altitude-corrected expression of the thrust. Denoting by T the thrust magnitude along the rocket body, we will consider that

T = g Isp q -S E P (h) (2.1)
where g is the gravity acceleration, Isp the engine specific impulse, q is the engine flow, S E the nozzle section surface and h the altitude.

5. A lift vector s.t. the effective aerodynamic forces on the rocket are contained in the plane defined by the relative speed vector and the rocket body longitudinal axis.

The functions M a , α → C Drag (M a , α) and M a , α → C Lift (M a , α) considered below for the 2D and 3D models convey the same aerodynamic model.

Remark 1.

Few articles have tackled PDG where drag is not negligible. See for instance [START_REF] Szmuk | Successive Convexification for Fuel-Optimal Powered Landing with Aerodynamic Drag and Non-Convex Constraints[END_REF], [START_REF] Leparoux | Structure of optimal control for planetary landing with control and state constraints[END_REF] or [START_REF] Brendel | Optimal guidance for toss back concepts of Reusable Launch Vehicles[END_REF].

Noteworthy particularities

Engine dynamics

The engine is assumed to generate the only control forces available on the rocket (i.e. aerodynamic grid fins or Reaction Control Systems (RCS) are not considered in this thesis). The output flow and the angular position of the engine are actuated. As far as the flow is concerned, its dynamic is not instantaneous and should not be neglected. It is assumed that its real flow is q r , that the controlled flow -or input signal -is q c , and that they are related via a first-order low-pass dynamics with time constant τ q s.t. Qualitative change of the wind flow with (left) and without (right) an air jet pushing in front of a moving rocket. See [START_REF] Nonaka | Vertical Landing Aerodynamics of Reusable Rocket Vehicle[END_REF] and [START_REF] Marwege | First Wind Tunnel Data of CALLISTO Reusable VTVL Launcher First Stage Demonstrator[END_REF] for experimental results.

qr = q c -q r τ q . ( 2
The dynamics of the angular position of the engine2 is not rigorously instantaneous. However, the time constant of its transient is sufficiently small compared to the other time constants of the guidance problem to be neglected. In the following, it will only be required that the attitude of the rocket (that the rocket control system will track) be a continuous function of time.

Remark 2. Consider some lower and upper bounds q -and q + on the engine flow.

If q r (0) ∈ [q -, q + ] and q c (t) ∈ [q -, q + ], then from Equation (2.2) we get that q r (t) remains within these bounds for all t ≥ 0. However, if q c (t) is not guaranteed to be between q -and q + , the proper version of Equation (2.2) becomes qr = Sat q + q -(q c )q r τ q .

Thrust direction

The thrust is assumed colinear to the rocket longitudinal axis. This assumption is an approximation, and corresponds to the equilibrium of the low-level controllers (out-of-the-scope of this thesis). More precisely, non-zero nozzle gimbal angles would slightly deviate the thrust vector, creating a momentum and in the end a rotation of the rocket. However, the rotation time constant is assumed significantly smaller than the translation time constant involved in the guidance problem.

Thrust dominance

The rocket engine is powerful compared to its weight, and the rocket incidence remains always sufficiently low s.t. the vertical speed is always negative and slowing down. Among others, this assumption prevents hovering maneuvers.

Dynamic equation and parameters choice

The dynamic equation of each model will be described by an ODE of the form

ẋ = f (x, u, η)
where x denotes the rocket states, u its controls, and η its dynamics parameters.

As will be detailed next, the rocket states are the position, the speed and the total mass, plus the real engine flow for the 3D rocket model only. The controls are the controlled engine flow, and one or two variables that convey the rocket incidence, depending on which model is considered (2D or 3D).

The states and controls chosen to describe each model will be made explicit. The choice of dynamics parameters presented in this Chapter is taken arbitrarily wide, to illustrate the various modeling possibilities. However, for the numerical examples of Chapters 4, 5 and 6, only relevant sub-sets of these parameters will be analyzed.

C N C A • α < 0 V r e Or e Vr V T D + L L D • θ > 0 mg w e z e h V = v h v z V r = v h v z -w(h)

Planar rocket model

The planar rocket model -a.k.a. 2D model -is described by its altitude h, its vertical speed v h , its horizontal position z, its horizontal speed v z and its total mass m. As mentioned above, v h ≤ 0.

The rocket is equipped with its own orthonormal frame (C A , C N ), where C A is parallel to the rocket body and oriented towards the engine, as pictured in Figure 2.2. To alleviate the writing, the vector pointing in the direction opposite to C A is noted e A := -C A .

The rocket orientation -or attitude -is defined by a single signed angle θ, formed by the angle between the vertical axis e h and the rocket main axis. Thus, when the rocket flies purely vertically, one has θ = 0 • . The incidence is defined as the signed angle α between the relative speed and the vector C A . The unit vector associated to the relative speed is denoted e Vr . The unit vector orthogonal to e Vr is denoted e Or , and is s.t. (e Vr , e Or ) is positively oriented, as shown in Figure 2

.2. A geometric relation gives tan(θ + α) = v z -w(h) |v h | .
As introduced in Section 2.1.1, the wind map is denoted w(h) and parametrized by the three values w z,0 , w z,1 and w z,2 . The equations of motion (EoM) are

ḣ = v h (2.3a) vh = -g + L sin θ + (T + D) cos θ m (2.3b) ż = v z (2.3c) vz = L cos θ -(T + D) sin θ m (2.3d)
where

T := gIspq r -P a (h)S E , (2.4a) V r := v 2 h + (v z -w(h)) 2 , ( 2.4b
)

D := 1 2 ρ(h)V 2 r S ref C Drag (M a , α), ( 2.4c
)

L := 1 2 ρ(h)V 2 r S ref C Lift (M a , α). (2.4d)
Adding the mass dynamics and using the fact that the vertical speed is assumed to be always negative (see Section 2.1.3.3), we can write the EoM

ḣ = v h vh = -g + ((T + D)|v h | + Lv z ) cos α + ((T + D)v z -L|v h |) sin α V r m ż = v z vz = (-(T + D)v z + F L |v h |) cos α + ((T + D)|v h | + Lv z ) sin α V r m ṁ = -q r
The 2D rocket variables are

(States) x := (h, v h , z, v z , m) ⊤ ∈ R 5 , (Controls) u := (q r , α) ⊤ ∈ R 2 , (Parameters) η := (∆Isp, w z,0 , w z,1 , w z,2 ) ⊤ ∈ R 4 ,
where the parameter ∆Isp is incorporated into Equation (2.4a) s.t. it becomes

T = g(Isp + ∆Isp)q r -S E P (h).
This yields the dynamic function f 2d of the 2D rocket model:

ẋ = f 2d (x, u, η).

Remark 3. Since the 2D rocket model serves illustrative purposes, it relies on the

simplifying assumption that the engine flow dynamics (2.2) is instantaneous, and that q c is its control variable. However, this engine dynamic will not be neglected in the 3D model.

Normal acceleration

The normal acceleration a nor is defined as the non-gravitational acceleration orthogonal to the relative speed. Thus, it equals

a nor = e Or •    vh + g vz    (2.6)
where • is the inner product. For the 2D model, a nor is signed.

Three-dimensional rocket model

The 3D rocket model has been designed to match the 2D model as closely as possible when the rocket trajectory remains in a plane. Thus, some concepts easily translate from one model to the other. Transposing the notions of incidence and attitude in 3D is, however, a delicate part. First, we introduce a series of frames that are necessary to define the rocket orientation in 3D. Then, we express the aerodynamic model, and formulate the associated EoM. Finally, some comments on the specifics of the 3D model are provided.

Orientation frames

The rocket is axially symmetric, making the notion of roll irrelevant. Two angles are used to describe its orientation.

First, we introduce the rocket yaw and pitch using Euler angles, which enables us to define a frame attached to the rocket body. Then, this frame is re-written using projected angles, which are more convenient for our applications.

Rocket orientation frame

As shown in Figure 2.3, the Earth's frame is (e z , e y , e h ). The rocket, initially3 positively colinear to e h , is oriented using the yaw Ψ first, and then using the pitch ξ, as explained in Figure 2.3. Mathematically, this translates into

C P := R (e z , Ψ) e y C N := R (C P , ξ) e z e A := R (C P , ξ) R (e z , Ψ) e h C A := -e A
where the vectors (C N , C A , C P ) define a new direct orthonormal frame, attached to the rocket body. Expressed in the frame (e z , e y , e h ), the latter gives

C N =       cos ξ sin ξ sin Ψ -sin ξ cos Ψ       , C A =       -sin ξ cos ξ sin Ψ -cos ξ cos Ψ       , C P =       0 cos Ψ sin Ψ       .
(2.7)

Projected angles

The above-defined angles Ψ and ξ define the orientation of e A . When projected onto the plane (e h , e z ) (respectively (e y , e h )), the vector e A has an angle ζ y (resp.

ζ z ) with the vector e h . Note that, in terms of vector labeling, the angle on the plane (e h , e z ) corresponds to a rotation on the axis e y . The angles are illustrated in Figure 2.4. The expression of e A writes

e A = 1 1 + tan 2 ζ z + tan 2 ζ y       tan ζ y -tan ζ z 1       (2.8)
and is valid for

|ζ z | < 90 • and |ζ y | < 90 • .
The steps used to derive formula (2.8) are as follows. Consider the square pyramid defined by a square on plane (e z , e y ) with its summit being at the top of e A . As shown in Figure 2.4, denote by v its height, and h z (resp. h y ) its base length orthogonal to e z (resp. e y ). Here, v, h z and h y are taken unsigned. Then, using basic geometry, one has

1 = v 2 + h 2 z + h 2 y , tan ζ z = h z v , tan ζ y = h y v which gives v = 1 1 + tan 2 ζ z + tan 2 ζ y , h z = v tan ζ z , h y = v tan ζ y .
Using the proper signs for e A , as shown in Figure 2

C N = 1 t z,y       1 cos ζz t y sin ζ z -t y cos ζ z       , C A = 1 t z,y       -t y t z -1      
and

C P =       0 cos ζ z sin ζ z      
. (2.10)

Dynamics

To express the aerodynamic forces, the orientation of the relative speed vector has to be defined. Then, lift and drag vectors are formulated. Finally, the dynamic equations are expressed.

Relative speed definition

To describe the orientation of the relative speed vector w.r.t. the rocket body, we introduce the incidence variables for the 3D model. The speed vector is

V := ( ż, ẏ, ḣ) ⊤ = (v z , v y , v h ) ⊤ .
The horizontal wind vector is w := (w z (h), w y (h), 0) ⊤ . Thus, the relative speed is

V r := V -w =       v z -w z (h) v y -w y (h) v h       .
The unit vector of the relative speed vector V r is denoted e Vr . We can define the orientation of e Vr using projected angles. The angles α y and α z are introduced s.t. the angles defining the position of e Vr w.r.t. the base frame using projected angles are ζ y + α y and ζ z + α z . They are represented in Figure 2.5. Since e Vr is defined equivalently as e A in Equation (2.8), its expression is

e Vr = 1 1 + tan 2 (α z + ζ z ) + tan 2 (α y + ζ y )       -tan(α y + ζ y ) tan(α z + ζ z ) -1       . (2.11)
Then, knowing the expressions of C A and e Vr , we can define the incidence as the unsigned angle between these vectors, leading to the expression

α = arcsin     (T z -t z ) 2 + (T y -t y ) 2 + (T z t y -T y t z ) 2 1 + t 2 z + t 2 y 1 + T 2 z + T 2 y     . ( 2
C A × e Vr = 1 (1 + t 2 z + t 2 y )(1 + T 2 z + T 2 y )       T z -t z t y -T y T y t z -T z t y       .
where × is the cross product. Then, from the definition of V r and the definition of the projected angles, we have

tan(α z + ζ z ) = v y -w y (h) |v h | and tan(α y + ζ y ) = - v z -w z (h) |v h | . ( 2
ζ z = -α z + arctan v y -w y (h) |v h | and ζ y = -α y -arctan v z -w z (h) |v h | . (2.14)

Aerodynamic effects

Lift and drag can be defined using (C N , C A , C P ). Indeed, the drag D is colinear to C A , in the opposite direction to V r . Only low-incidence flight is considered, and D is positively colinear with -C A . Moreover, in consistency with the aerodynamic model described in Section 2.1.2, the axial symmetry of the 3D rocket model implies that L, C A and V r are linearly dependent, and that L must belong to the plane (C P , C N ). Therefore, as illustrated in Figure 2.6, one can define the lift and drag vectors as

L = L.e L and D = -D.C A ,
where the magnitudes L and D equal

L = 1 2 ρ(h) V 2 r S ref C Lift (M a , α) and D = 1 2 ρ(h) V 2 r S ref C Drag (M a , α).
The direction e L of the lift requires further attention, since the lift orientation is defined only when the incidence is not zero. When it is well defined, the vector e L is a unit vector, positively colinear to the projection of -e Vr lying in the plane (C P , C N ). Thus, it can be expressed as

e L := - (e Vr • C N )C N + (e Vr • C P )C P ∥(e Vr • C N )C N + (e Vr • C P )C P ∥ = - e Vr -(e Vr • C A )C A ∥e Vr -(e Vr • C A )C A ∥ (2.15)
when e Vr is not colinear to C A and e L = 0 otherwise. This apparent discontinuity is actually not troublesome. Indeed, for α > 0, the vector e L can be equivalently defined as

e L = C A × C A × e Vr ∥C A × e Vr ∥ = C A × (C A × e Vr )
sin α which yields the following expression for the lift

L = L.e L = 1 2 ρ(h)V 2 r S ref C Lift (M a , α) sin α C A × (C A × e Vr ).
For any fixed Mach number M a , the map α → C Lift (M a , α) is assumed continuously differentiable and it equals zero at α = 0. Thus, the ratio C Lift (Ma,α) sin α remains bounded when α tends to zero. Also, C A × e Vr tends to zero when e Vr tends to C A . Consequently, the expression of e L does not matter when α = 0, which is a false singularity.

Remark 5.

Here, the coefficient C Lift is positive and only evaluated for positive values of α. However, note that in the 2D model, α → C Lift (M a , α) is taken odd for any fixed Mach number M a and is evaluated on signed values of α.

Rocket dynamics in 3D

With g = (0, 0, -g) ⊤ denoting the gravity vector, the acceleration vector a equals where the thrust vector is T = -T C A and its magnitude T is defined in Equation (2.1).

a := d dt V = g + T + L + D m C N C A C P D V r e Vr
The dynamics parameters conveyed by the variable η are the Isp via ∆Isp, multiplicative factors for the aerodynamic coefficients (m L , m D ) and the wind parameters (w z,0 , w z,1 , w z,2 ) and (w y,0 , w y,1 , w y,2 ) defined in Section 2.1.1. The first three parameters must be incorporated in the equations s.t. T = gIspq r -S E P (h) becomes g(Isp + ∆Isp)q r -S E P (h)

L = 1 2 ρ(h)V 2 r S ref C Lift (M a , α) becomes 1 2 ρ(h)V 2 r S ref (1 + m L )C Lift (M a , α) D = 1 2 ρ(h)V 2 r S ref C Drag (M a , α) becomes 1 2 ρ(h)V 2 r S ref (1 + m D )C Drag (M a , α) It allows us to define the 3D rocket variables as (States) x := (z, y, h, v z , v y , v h , m, q r ) ⊤ ∈ R 8 , (Controls) u := (q r , α z , α y ) ⊤ ∈ R 3 , (Parameters) η := (∆Isp, m L , m D , w z,0 , w z,1 , w z,2 , w y,0 , w y,1 , w y,2 ) ⊤ ∈ R 9
Then, written using blocks, the dynamic equation equals

ẋ = f 3d (x, u, η) =           V g + T+L+D m -q r q c -q r τ q           . ( 2

Normal acceleration, downrange and attitude

As it will be needed later in our work, we now focus on the acceleration component normal to the relative speed vector. The unsigned normal acceleration a u.s. nor is defined as the norm of the part of the non-gravitational acceleration vector normal to the relative speed. Denoting F = T + L + D yields Compared to the 2D model, this expression is naturally unsigned. However, the norm in this expression may bring differentiation issues. Indeed, it will be needed to differentiate this term later when considering optimization problems (see e.g. (4.6g) in Chapter 4). Instead of a u.s. nor , we consider an alternate expression that remains signed (and thus differentiable), see below Equation (2.17).

First remark that, according to the aerodynamic model and as shown in (2.17) s.t. a u.s. nor = |a nor |. Thanks to Equation (2.17), the 3D model also has a term a nor describing the normal acceleration that is signed.

Downrange

The downrange, or distance of the vertical projection of the rocket to the landing site, equals d := z 2 + y 2 .

Attitude

The attitude θ is defined as the angle between e h and e A . It is unsigned (contrary to the planar rocket model). Its expression equals θ := arcsin sin 2 ξ + cos 2 ξ sin 2 Ψ.

(2.18)

PDG as an Optimal Control Problem

First, we present the general PDG objectives and then write the main OCPs.

Mission goals and constraints

In this sub-section, we present all the constraints of our PDG problem.

Landing site target

It is desired to find a trajectory -i.e. states x(t) and controls u(t) -that steers the rocket from a given initial condition x 0 to the landing site in a time t f . Assuming that the landing site is located at the origin of our coordinate system, the following end-point constraints are introduced z(t f ) = 0 y(t f ) = 0 (For the 3D model only)

h(t f ) = 0 v z (t f ) = 0 v y (t f ) = 0 (For the 3D model only) v h (t f ) = -ε f v
Note that a small non-zero vertical speed at landing (ε f v > 0) is desired. Among others, this serves to avoid singularity of incidence at t f and gives a margin regarding the thrust dominance assumption 4 . These conditions can be written as a linear equality constraint

A f x(t f ) = b f
where A f is filled with zeros and ones only, and b f contains zeros and -ε f v only.

Mechanical bounds

As noted in Section 2.1.3, the rocket has several mechanical limitations. Its engine flow magnitude and rate of change are bounded. These impose that the real and controlled flows satisfy q -≤ q r ≤ q + and q -≤ q c ≤ q + , (

where q + > q -> 0.

Remark 7. It should be noted that in the literature, the decision variables often considered (especially for non-atmospheric missions) are the components of the thrust vector T itself. The bounded flow constraints are then taken into account through the constraint T min ≤ ∥T∥ ≤ T max . In the latter, the lower bound defines an artificially non-convex constraint. Lossless convexification [START_REF] Açıkmeşe | Lossless convexification of a class of optimal control problems with non-convex control constraints[END_REF][START_REF] Carson | Lossless convexification of Powered-Descent Guidance with non-convex thrust bound and pointing constraints[END_REF], a method using an intermediate slack variable, is often used to overcome this problem. However, with our modeling choices, picking the engine flow as one of the decision variables makes the constraints described in Equation (2.19) convex, since the constraints of the shape "q -≤ q ≤ q + " are used for scalar values of q.

Bounded mass

The fuel tank being finite, there are upper and more importantly lower bounds on the mass

m dry ≤ m ≤ m wet .
(2.20)

Safety bounds

For safety reasons, it is also desirable to remain within limited incidences. This constraint is one of the most important difference between planetary and atmospheric landing. It writes

|α| ≤ α max .
This constraint has a straightforward interpretation for the 2D model, since it only implies a single control variable, α itself. However, it is more intricate for the 3D model, since α is defined through Equation (2.12) and depends non-linearly on state and control variables. We choose to impose the following constraints for the 3D model

|α z | ≤ α max and |α y | ≤ α max ,
which is not equivalent but conservative. Moreover, since the guidance trajectory must be tracked by the underlying rocket control system, it is necessary that this trajectory does not exceed prescribed thresholds of normal accelerations. Thus, we impose

|a nor | ≤ a max nor
where a nor is defined in Equation (2.6) for the 2D rocket model, and in Equation (2.17) for the 3D one.

Constraints not considered

Other types of constraints can be found in the literature, such as5 :

• Landing cone (a.k.a. glideslope) constraints [START_REF] Açikmeşe | Flight Testing Of Trajectories Computed By G-FOLD: Fuel Optimal Large Divert Guidance Algorithm For Planetary Landing[END_REF][START_REF] Eren | Constrained Reachability and Controllability Sets for Planetary Precision Landing via Convex Optimization[END_REF][START_REF] Szmuk | Successive Convexification for Real-Time 6-DoF Powered Descent Guidance with State-Triggered Constraints[END_REF],

• Pointing (a.k.a. attitude) constraints [START_REF] Eren | Constrained Reachability and Controllability Sets for Planetary Precision Landing via Convex Optimization[END_REF],

• Thermal flux (a.k.a. heating rate) constraints [START_REF] Brendel | Optimal guidance for toss back concepts of Reusable Launch Vehicles[END_REF][START_REF] Wang | A Pseudospectral-Convex Optimization Algorithm for Rocket Landing Guidance[END_REF][START_REF] Wang | Optimal Rocket Landing Guidance Using Convex Optimization and Model Predictive Control[END_REF],

• Dynamic pressure constraints [START_REF] Brendel | Optimal guidance for toss back concepts of Reusable Launch Vehicles[END_REF][START_REF] Wang | Optimal Rocket Landing Guidance Using Convex Optimization and Model Predictive Control[END_REF].

However, as detailed below, these constraints are not considered here, though this would be possible as natural extensions.

Landing cone constraints are not considered due to the thrust dominance assumption. Indeed, the class of rockets studied in this thesis naturally performs landing trajectories with a high slenderness ratio. The same motivation rules out the pointing constraints.

Thermal flux constraints are critical mechanical requirements for re-entry problems at hypersonic speeds, when the vehicle directly relies on the atmosphere to brake [START_REF] Bonnard | Optimal control of the atmospheric arc of a space shuttle and numerical simulations with multiple-shooting method[END_REF]. Since the speeds involved for our landing scenarios are high (approximately between 0 and Mach 2) but not hypersonic6 and considering the unusual air flow around the rocket, as depicted in Figure 2.1, the thermal flux constraints are not needed. The same arguments apply for the dynamic pressure constraint.

Formulation as an optimal correction problem

The PDG problem is formulated as an OCP in free-final time, w.r.t. a reference trajectory.

Regarding the notations, whatever the chosen model is (2D or 3D), the dynamic function is noted f . The lower and upper control bounds are respectively denoted u - and u + . The mixed state-control constraints are conveyed by a function c.

As mentioned in the Introduction, we consider that all of the above-mentioned constraints can be tuned by a parameter p that allows one to adjust their nominal value. For example, let us say that we parametrize the incidence and the normal acceleration bounds in the planar rocket model, i.e. p = (p 1 , p 2 ) ⊤ = (∆α max , ∆a max nor ) ⊤ . Then, the parameterized constraints become

|α| ≤ α max + p 1 and |a nor | ≤ a max nor + p 2 .
From a general point of view, we say that all the constraints are parameterized by p. Therefore, the dependency on p of the right-hand side vector b f , the control bounds u -and u + , and the mixed state-control constraint c will be highlighted whenever necessary. It is assumed that p = 0 for a nominal landing. The exact choice of p varies depending on the mission needs, and will be discussed extensively in Chapter 5.

Let us consider a reference trajectory, which is described as a quadruplet (x, ū, η, tf ). Such a trajectory can be computed offline, for any given mission, using well-known and possibly time-consuming numerical methods [START_REF] Bryson | Applied optimal control: optimization, estimation and control[END_REF]. It is assumed that this trajectory satisfies the constraints, i.e. that it satisfies the control constraints and c(x(t), ū(t), η, p) ≤ 0 for all times, and that it is dynamically feasible, i.e. that it satisfies the Initial Value Problem (IVP)

     x(0) = x0 , ẋ(t) = f (x(t), ū(t), η), ∀t ∈ [0, tf ].
(2.21)

At the beginning of the final burn, the gap between the current state and the current dynamics parameters and their reference values is denoted ∆x 0 and ∆η. They are conveyed by the input variable7 ξ.

Knowing ξ, and the value of p, the mathematical goal of PDG, when formulated w.r.t. this reference trajectory, is to find a control correction δu and a time-offlight correction ∆t f making the rocket land while satisfying all the constraints and minimizing a certain performance index8 J . For completeness, note that δu belongs to a functional space defined over [0, tf + ∆t f ]. We will consider that this space9 equals U(∆t f ) := L ∞ ([0, tf + ∆t f ], R m ) for the definition below, although we will restrict the problem to a much more specific class of control corrections in Chapter 4.

Definition 1 (Infinite dimensional problem, PDG (ξ, p)). Given a reference de-

scribed by its quadruplet (x, ū, η, tf ), where ū is defined over [0, tf ], find the optimal time-of-flight change ∆t f and the optimal control correction δu ∈ U(∆t f ) for the optimization problem PDG (ξ, p) defined by

min δu,∆t f J (δu, ∆t f ) (2.22a) s.t. ẋ(t) = f (x(t), ū(t . tf /( tf + ∆t f )) + δu(t), η + ∆η), (2.22b
)

x(0) = x0 + ∆x 0 (2.22c) A f x( tf + ∆t f ) = b f (p) (2.22d) u -(p) ≤ ū(t . tf /( tf + ∆t f )) + δu(t) ≤ u + (p), (2.22e) c(x(t), ū(t . tf /( tf + ∆t f )) + δu(t), η + ∆η, p) ≤ 0, (2.22f)
where conditions (2.22b), (2.22e) and (2.22f) are meant for all t ∈ [0, tf + ∆t f ].

Remark 8. J is assumed strictly convex. Moreover, for null inputs (i.e. ∆x 0 = 0 and ∆η = 0) it has a null minimum (i.e. δu * = 0 and ∆t * f = 0). Remark 9. Written this way, the formulation of PDG (ξ, p) suggests that the bounded mass constraint is enforced over the whole interval [0, t f ]. However, since the engine flow is positive, the mass is decreasing, and thus it is only necessary to enforce the simpler condition m( tf + ∆t f ) ≥ m dry in practice.

Written in this format, PDG (ξ, p) falls into the field of perturbation methods for OCPs [START_REF] Bensoussan | Perturbation methods in optimal control[END_REF][START_REF] Deshpande | Directional Input Adaptation in Parametric Optimal Control Problems[END_REF], which is the ground base of Chapter 4. However, this problem is studied for any values of ∆x 0 and ∆η, because we are interested in the solutions of this problem even for non small values of these inputs.

Summary

In this chapter, we have defined dynamic models of the rocket, for 2D and 3D motions. The aerodynamic model is the main difference between the planetary landing problems studied in the literature and atmospheric landing problems. The general, infinite-dimensional version, of the PDG problem has been defined as finding the optimal correction (δu, ∆t f ) w.r.t. a reference trajectory, while satisfying all the problem constraints, for given values of ξ and p.

A simpler version of PDG (ξ, p), considering only vertical motion, is considered in Chapter 3 (which can be skipped without loss of continuity), and Chapter 4 presents a method to compute an approximation of the solution of the general problem.

Chapter 3

Mathematical properties of the optimal vertical descent

Résumé

Ce chapitre se concentre sur un cas particulier du problème PDG présenté précédemment : l'optimisation de la consommation de carburant pour l'atterrissage vertical. Comme nous l'avons déjà vu, le rapport entre la translation latérale et la translation verticale est faible. En tant que cas limite, il est tentant d'étudier d'abord le problème purement vertical. Se concentrer sur le problème purement vertical permet plusieurs simplifications : il n'y a plus qu'une seule variable de décision (le débit du moteur) et le modèle aérodynamique est grandement simplifié. Cela rend l'analyse analytique du problème envisageable.

Le problème de l'atterrissage atmosphérique vertical optimal en termes de carburant est ici étudié en tant que problème de Commande Optimale en temps final libre. La principale contribution établit la nature de la loi de poussée optimale en consommation de carburant, en étendant les résultats de la littérature sur les problèmes sans atmosphère. Des conditions suffisantes et nécessaires sont fournies qui garantissent la nature Min-Max des extremums normaux. Il est également démontré que les extremales anormales -au sens du principe du maximum de Pontryagin (PMP) -sont soit Min, soit Max. Un sous-produit utile de cette étude est une caractérisation de l'ensemble atteignable pour les atterrissages verticaux. Cette notion sera réutilisée plus tard dans le chapitre 6 à des fins d'évaluation des performances.

This chapter focuses on a special case of the PDG problem presented previously: the fuel-optimal vertical landing. As already discussed, the ratio of lateral vs vertical translation is small. As a limit case, it is tempting to study the limit case of the purely vertical problem first. Focusing on the purely vertical problem enables several simplifications: there is only one decision variable left (the engine flow) and the aerodynamic model is greatly simplified. This makes the analytic analysis of the problem tractable.

The vertical fuel-optimal vertical atmospheric landing problem is here studied 35 as a free-final time OCP. The main contribution establishes the nature of the fueloptimal thrust program, extending results from the literature on atmosphere-free problems. Sufficient and necessary conditions are provided that guarantee the Min-Max nature of the normal extremals. Abnormal extremals -in the sense of the Pontryagin Maximum Principle (PMP) -are also shown to be either Min or Max. A useful by-product of this study is a characterization of the reachable set for vertical landings. This notion will be re-used latter in Chapter 6 for performance assessment purposes.

If necessary, the reader can skip directly to page 49 for a summary of the important results. Also, note that this chapter is a detailed version of [START_REF] Ménou | Fuel-optimal program for atmospheric vertical powered landing[END_REF].

Vertical descent

Historically, Meditch [START_REF] Meditch | On the problem of optimal thrust programming for a lunar soft landing[END_REF] and then Shi & Eckstein [START_REF] Shi | An exact solution for optimum controlled soft lunar landing[END_REF] have offered analytic solutions for the (atmosphere-free) vertical Moon landing problem. Since then, due to the spectacular development of reusable launcher technologies, powered landing strategies have been successfully addressed using numerical methods [START_REF] Blackmore | Lossless convexification of control constraints for a class of nonlinear optimal control problems[END_REF][START_REF] Brendel | Optimal guidance for toss back concepts of Reusable Launch Vehicles[END_REF][START_REF] Lee | Constrained Autonomous Precision Landing via Dual Quaternions and Model Predictive Control[END_REF][START_REF] Ross | A review of pseudospectral optimal control: From theory to flight[END_REF][START_REF] Szmuk | Successive Convexification for Real-Time 6-DoF Powered Descent Guidance with State-Triggered Constraints[END_REF][START_REF] Ulybyshev | Optimization of Three-Dimensional Lunar Landing Trajectories and Accessible Area Computation[END_REF]. Due to the non-negligible effects of the atmosphere, the analytic results derived for the Moon landing problem can not be directly adapted to the problem of Earth landing. Yet, analytical results on this problem would still represent valuable assets. On the one hand, analytic solutions are very useful to assess the quality of the numerical methods, by providing well-described reference solutions to standardized problems, see e.g. [START_REF] Bonnard | Optimal control of the atmospheric arc of a space shuttle and numerical simulations with multiple-shooting method[END_REF][START_REF] Souza | An optimal guidance law for planetary landing[END_REF][START_REF] Goddard | A Method Of Reaching Extreme Altitudes[END_REF][START_REF] Graichen | Solving the Goddard problem with thrust and dynamic pressure constraints using saturation functions[END_REF][START_REF] Reynolds | Optimal Planar Powered Descent with Independent Thrust and Torque[END_REF]. Further, when analytical investigations establish the switching structure of the solution, very efficient numerical methods can be employed, using a reduced number of unknown variables [START_REF] Souza | An optimal guidance law for planetary landing[END_REF][START_REF] Reynolds | Optimal Planar Powered Descent with Independent Thrust and Torque[END_REF][START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF]. For complex dynamics and high-dimensional systems, obtaining such analytical results is usually considered as out-of-reach [START_REF] Schättler | Geometric Optimal Control: Theory, Methods and Examples[END_REF]. It is thus of importance to select only dominant factors while leaving out unnecessary details in the modeling. Following this modus operandi, we consider a simplified (but not simplistic) representation of the general powered landing problem and establish a non-trivial result.

The analysis presented in this chapter considers one key element: the effects of atmosphere. The model under study builds upon the variable-mass model of a rocket considered in [START_REF] Meditch | On the problem of optimal thrust programming for a lunar soft landing[END_REF] and incorporates atmospheric effects in the form of an altitudedependent bias of the thrust only, as introduced in Equation (2.1) in Chapter 2. In this model of the final phase of the powered landing, the thrust generator is always turned on1 and the thrust is upper and lower-bounded in a way that prevents hovering flight (according to the Thrust dominance assumption already presented).

Intuitively, one could expect that it is more efficient to wait until the last feasible moment to use maximal thrust, as early efforts trying to slow down the rocket are likely to be less effective due to the varying mass scaling of the dynamics. The contribution of this chapter is to establish conditions under which fuel-optimal vertical powered landing through the atmosphere is indeed of this expected Min-Max nature.

The arguments of proof are as follows. Under simple assumptions on the atmosphere pressure model (decreasingness, convexity), the optimal thrust program is first shown to have a Max-Min-Max structure, based on the PMP. Compared to [START_REF] Meditch | On the problem of optimal thrust programming for a lunar soft landing[END_REF], some sharper differential inequalities on the adjoint states are necessary to obtain a conclusion. Also, both normal and abnormal extremals need to be tackled. Then, using additional inequality constraints derived from the Implicit Function Theorem (IFT), Min-Max structures are proven to be more fuel-optimal than Max-Min-Max structures. These conditions can be checked numerically, over a finite domain. It is also shown that these conditions hold for zero atmosphere (and scarce atmosphere, using a continuity argument), which makes a connection with [START_REF] Meditch | On the problem of optimal thrust programming for a lunar soft landing[END_REF].

The chapter is organized as follows. In Section 3.1, the dynamics and the powered landing problem are summarized, in harmony with Chapter 2. In Section 3.2, the flight envelope is described based on flow analysis and differential inequalities. In Section 3.3, the optimal thrust program is shown to be Min-Max using the PMP, the IFT and mild assumptions. Finally, we provide numerical details in Section 3.4, and concluding remarks in Section 3.5.

Single dimensional rocket model

Following Chapter 2, we describe the rocket having a purely vertical motion by its altitude h, speed v and total mass m. The dynamics write

ḣ = v, v = -g + T (h, q) m , ṁ = -q
where q is the engine flow, and the thrust T defined in Equation (2.1). Because the engine is firing, the sole effect of the atmosphere is through the atmospheric pressure in the expression of T . Recall that negative speed conveys descending movement.

During the powered landing, the rocket engine is always firing and q is bounded.

In the problem setup under consideration, the engine flow bounds are s.t. the net thrust is always positive, i.e. q -is s.t.

a cc := g Isp q --S E max h≥0 P a (h) > 0.
(3.1)

Normalized dynamics

The following normalized variables are introduced

u := 2 q -q - q + -q --1, y 1 := h g Isp , y 2 := v g Isp , y 3 := 2m q + -q -
where y := (y 1 , y 2 , y 3 ) ⊤ denotes the normalized2 states and where

κ := 1 Isp , r := q + + q - q + -q -, π(y 1 ) := P a (g Isp y 1 ) 2S E g Isp (q + -q -) .
This yields the control-affine dynamics in

R 3 ẏ = f (y) + ug(y), ( 3.2)
where |u| ≤ 1 and

(Altitude) ẏ1 = y 2 (3.3a) (Speed) ẏ2 = r + u -π(y 1 ) y 3 -κ (3.3b) (Mass) ẏ3 = -(r + u). (3.3c)
Also, recall that the mass is bounded s.t. m -≤ y 3 ≤ m + .

Assumptions specific to the vertical descent

The problem under study is also described by the two following assumptions.

Assumption 1 (Pressure model properties). The normalized pressure function π,

is of class C 2 , and π > 0, π ′ < 0, π ′′ > 0.
This assumption is very general and holds for all reference Earth atmosphere models, such as [START_REF] Leslie | The NASA Marshall Space Flight Center Earth Global Reference Atmospheric Model, 2010 Version[END_REF]. Then, the thrust dominance assumption from Chapter 2 is re-written as follows.

Assumption 2 (Thrust dominance). ẏ2 ≥ a cc > 0.

Assumption 2 implies condition (3.1), shows that the ratio r is greater than 1 and it also prevents hovering. Reaching null speed at a positive altitude is thus an undesired behavior and is not a steady state.

Optimal Control Problems

A natural goal for rocket landing is to maximize the final mass [START_REF] Meditch | On the problem of optimal thrust programming for a lunar soft landing[END_REF], or equivalently to minimize the fuel consumption. Landing is defined as final null altitude and (vertical) velocity. A constrained optimal control problem in free final time depending on an initial state y 0 can then be formulated.

Problem 1 (Fuel optimal landing with state inequality path constraints).

min u(.),t f t f 0 r + u(s)ds (3.4a) s.t. ẏ = f (y) + ug(y), (3.4b) |u| ≤ 1 (3.4c) y(0) = y 0 , y 1 (t f ) = y 2 (t f ) = 0 (3.4d) y 1 (t) ≥ 0, y 2 (t) ≤ 0, y 3 (t) ∈ [m -, m + ] (3.4e)
State constraints (3.4e) are meant for any t in [0, t f ]. Additionally, we will consider another formulation where the state constraints (3.4e) have been removed, as they will be shown to be automatically satisfied.

Problem 2 (Fuel optimal landing).

min u(.),t f t f 0 r + u(s)ds s.t. ẏ = f (y) + ug(y), |u| ≤ 1 y(0) = y 0 , y 1 (t f ) = y 2 (t f ) = 0
Studying Problem (2) will help us describe the solutions of Problem 1.

Remark 10 (Terminology). In the following proofs, maximal solutions of an ordinary differential equation are the solutions that cannot be extended in time.

Premilinaries on the dynamics

This section aims at describing conditions under which state path inequalities (3.4e) can be ignored. A detailed study of the dynamics is conducted. First, the altitude and speed dynamics are studied using surfaces of R3 s.t. any trajectory that lands must start between these surfaces. This region is called the flight envelope. Then, the mass constraint is discussed.

Let us denote the domain

D := R + × R -× (0, m + ].
Below, we say that a trajectory starting at some y 0 ∈ D lands applying the thrust u(.) if it reaches

y 1 (t f ) = y 2 (t f ) = 0 for some t f > 0.
Note that the minimum mass constraint is not included in the first part of this discussion.

Beforehand, remark that there is a unique time T u associated to a control u(.) s.t.

y 0 3 - Tu 0 r + u(s)ds = 0. ( 3.5)
Since r + u ≥ r -1 > 0, the map t → 1/y 3 (t) is not integrable near T u because of (3.3c). Thus, the maximal solution of (3.2) starting at

y 0 ∈ D is defined on the interval [0, T u ). If u ≡ σ is constant, then T σ = y 0 3 /(r + σ). Lemma 1. Let σ ∈ [-1, 1]
a constant parameter. For any y 0 2 and y 0 3 , there is a unique y 0 1 (σ, y 0 2 , y 0 3 ) s.t. a trajectory starting at (y 0 1 (σ, y 0 2 , y 0 3 ), y 0 2 , y 0 3 ) ⊤ ∈ D lands when applying the constant thrust u ≡ σ.

Proof. The maximal solution y of (3.2) with u ≡ σ, starting at y 0 ∈ D, is defined on [0, T σ ). y 2 (.) is continuous, increasing and diverges to +∞ as t tends to T σ . Thus, there is a unique time, denoted t * (y 0

1 ) ∈ [0, T σ ) s.t. y 2 (t * (y 0 1 )) = 0. The IFT applied with Assumption 2, on equation 3 Φ f +σg t * (y 0 1 ), (y 0 1 , y 0 2 , y 0 3 ) ⊤ 2 = 0 (3.6)
shows that the application that maps y 0 1 into t * is actually continuous, and differentiable, for all y 0 1 ≥ 0. Then, define

η : z ∈ R + → Φ f +σg (t * (z), (z, y 0 2 , y 0 3 ) ⊤ ) 1 ∈ R. (3.7)
From the regularity of f + σg, the flow Φ f +σg is continuous and thus η is continuous. Necessarily, η(0) < 0. Moreover, since the acceleration is lower-bounded by a cc , it is possible to find an altitude y crit 1 > 0 large enough s.t. η(y crit 1 ) > 0. Therefore, there is a y * 1 ≥ 0 s.t. η (y * 1 ) = 0. Using t f = t * (y * 1 ), one has y 1 (t f ) = y 2 (t f ) = 0 by construction of η. A comparison argument, as in the proof of Proposition 1, shows that η is actually increasing, proving the uniqueness of y * 1 . It yields y * 1 = y 0 1 (σ, y 0 2 , y 0 3 ) using the above-mentioned variables.

Let us denote Σ max (respectively Σ min ) the set of initial conditions s.t. landing is successful, at mass y f 3 ∈ (0, m + ], when applying a constant maximum (resp. minimum) thrust. Denoting

y max 1 (y 2 , y 3 ) := y 0 1 (1, y 2 , y 3 ), y min 1 (y 2 , y 3 ) := y 0 1 (-1, y 2 , y 3 ), yields Σ max := {(y max 1 (y 2 , y 3 ), y 2 , y 3 ) : y 2 ≤ 0, y 3 ∈ (0, m + ]}, Σ min := {(y min 1 (y 2 , y 3 ), y 2 , y 3 ) : y 2 ≤ 0, y 3 ∈ (0, m + ]}.
Moreover, using flows of the backward-time dynamics, for σ ∈ [-1, 1], define

Σ σ := Φ -(f +σg) t, (0, 0, y f 3 ) ⊤ : 0 ≤ t ≤ m + -y f 3 r + σ , 0 < y f 3 ≤ m +
which provides the relations Σ max = Σ 1 and Σ min = Σ -1 . It is noteworthy that

y max 1
(y 2 , y 3 ) ≤ y min 1 (y 2 , y 3 ), implying that Σ max is always "below" Σ min , as pictured in Figure 3.2.

Note that the applications y min 1 and y max 1 are continuous: this property stresses the continuity of the flows and the formal definition of Σ σ . Continuity can also be proven using the IFT on function η from Equation (3.7), considering y 0 2 and y 0 3 as variables.

Proposition 1. For any y 0 ∈ D, if y min 1 (y 0 2 , y 0 3 ) < y 0 1 then for any control u(.) in [-1, 1] the dynamics reaches null speed at a positive altitude.

Proof. Consider y 0 ∈ D s.t. y min 1 (y 0 2 , y 0 3 ) < y 0 1 and denote ỹ0 := (y min 1 (y 0 2 , y 0 3 ), y 0 2 , y 0 3 ) ⊤ .

Let ỹ be the maximal solution of (3.2) with u ≡ 1 and y be the maximal solution of (3.2) for some measurable function u satisfying |u| ≤ 1 at all times. y starts at y 0 and ỹ at ỹ0 . They are respectively defined on [0, T 1 ) and [0, T u ), where

T u ≤ T 1 .
Using mass as a time-varying scaling, we get

   ẏ1 (t) ẏ2 (t)    ≥ K   t,    y 1 (t) y 2 (t)       :=    y 2 (t) -κ + r-1-π(y 1 (t)) y 0 3 -t(r-1)   
for any t ∈ [0, T u ). By construction, ỹ satisfies the equality version of this equation. Thus, comparison Lemma 13 (in Appendix) yields

ỹ1 (t) ≤ y 1 (t), ỹ2 (t) ≤ y 2 (t), ∀t ∈ [0, T u ).
Since y 2 is continuous, increasing and diverges as t → T u , there is a unique t * ∈ [0, T u ) s.t. y 2 (t * ) = 0. Therefore, y 1 (t * ) ≥ ỹ1 (t * ) ≥ 0. Using a Taylor expansion on (3.3a) with the initial conditions shows that the last inequality is strict, whence the proposition.

Using a very similar proof, one shows the following result.

Proposition 2. For any y 0 ∈ D, if y max 1 (y 0 2 , y 0 3 ) > y 0 1 then for any control u(.) in [-1, 1] the dynamics reaches null altitude at a negative speed.

Proposition 1 defines the notion of being too high, meaning that if the rocket starts its powered descent above Σ min (in terms of altitude), then it will either lack fuel before reaching null speed, or go back up before touching the ground and then lack fuel at a non-zero altitude. In both cases, landing fails. Proposition 2 is the exact equivalent for the notion of being too low, meaning that the rocket will hit the ground at a non-zero speed if it starts below Σ max .

Further, note that if a trajectory lands s.t. the mass remains in [m -, m + ], then the acceleration is upper-bounded by ācc := -κ + r+1 m -for any positive time. Since the fuel flow is lower-bounded, the mass can remain in [m -, m + ] for at most

T max := m + -m - r-1
Therefore, for any positive time, the speeds are lower-bounded by y 2 and the altitudes are upper-bounded by ȳ1 s.t. Consequently, if y 0 ∈ F and Problem (2) has a solution, then altitude and speed constraints are enforced. If y 0 ∈ D\F, then Problems 1 and 2 cannot have solutions.

Leaving out the limit cases of Σ max and Σ min , for which landing can be achieved by applying, respectively, the maximum and the minimum thrust, for the whole duration of the flight, we introduce

F * := F\(Σ max ∪ Σ min ).
(3.10)

The following result discusses feasibility of the landing. Optimality will be studied later on in Section 3.3.

Proposition 3. If y 0 belongs to F * , then there is always a control u of structure Min-Max that lands.

Proof. The Min-Max structure denotes a 2 step sequence starting with minimum value of the control and ending with maximum value. For such y 0 , denote ỹ0 := (y min 1 (y 0 2 , y 0 3 ), y 0 2 , y 0 3 ) ∈ Σ min . Since y 0 ∈ F * , then y max 1 (y 0 2 , y 0 3 ) < y 0 1 < y min 1 (y 0 2 , y 0 3 ). Let us denote y and ỹ the maximal solutions of Equation (3.2) with u ≡ -1, starting respectively at y 0 and ỹ0 . Then, using similar comparisons as in the previous proof, one obtains y 1 (t) < ỹ1 (t) and y 2 (t) < ỹ2 (t) for all positive times. Thus, one deduces that y 1 reaches zero at some time t ′ > 0, before y 2 does. Moreover, the map

ξ : t ∈ [0, t ′ ] → y 1 (t) -y max 1 (y 2 (t), y 3 (t)) ∈ R (3.11)
is continuous, and satisfies ξ(0) > 0 since the trajectory starts strictly above Σ max , and ξ(t ′ ) < 0 since (0, y 2 (t ′ ), y 3 (t ′ )) is necessarily below Σ max in terms of altitude (recall that y 2 (t ′ ) < 0). Thus, there exists a time t ′′ ≤ t ′ s.t. y(t ′′ ) ∈ Σ max . The desired Min-Max control law equals -1 on [0, t ′′ ) and +1 for times t ≥ t ′′ .

As far as the mass is concerned, since it is a continuous decreasing function of time, enforcing the terminal constraint y 3 (t f ) ≥ m -is sufficient to guarantee the mass constraint (3.4e).

Therefore, only the simplified Problem (2) needs to be solved. If there is a solution that satisfies y 3 (t f ) ≥ m -, then Problem (1) shares the same solution. Otherwise, if y 3 (t f ) < m -, then Problem (1) has no solutions. Indeed, since the solution is fuel-optimal, there is no other way to land with a greater final mass.

Optimal thrust programs

This section focuses on Problem (2) exclusively, which, according to the previous discussion gives an answer to Problem (1) or proves its infeasibility. We aim at proving that optimal controls are of Min-Max nature, where one of the min or max arcs may be absent. To establish this result (Theorem 1), we proceed as follows. First, stationary conditions are derived from the PMP. Then, using properties of the second adjoint state variable, the optimal thrust program is shown to be Max-Min-Max. Finally, the first maximum arc is shown to be absent under one (mild) additional assumption on the atmosphere model (Assumption 4).

Fuel Optimal Landing

Consider y 0 ∈ F. Let u be an optimal thrust program for Problem (2) and y be the corresponding trajectory. Let t f be the time-of-flight. It is assumed that y 3 (t f ) ≥ m -. Thus, from the previous section, y lies in F.

The Hamiltonian of Problem ( 2) is defined as

H := λ 0 (r + u) + λ ⊤ (f (y) + ug(y)) (3.12)
where λ 0 ∈ R and λ : [0, t f ] → R 3 denote the adjoint states. To study the controlaffine Hamiltonian, consider the switching function

Γ(t) := λ 0 + λ(t) ⊤ g(y(t)) = λ 0 + λ 2 (t) y 3 (t) -λ 3 (t). (3.13)
The PMP, as stated in [84, Thm. 2.2.1], yields

(λ 0 , λ(t)) ̸ = 0 R 4 , ∀t ∈ [0, t f ] (3.14a) λ1 = λ 2 π ′ (y 1 ) y 3 (3.14b) λ2 = -λ 1 (3.14c) λ3 = λ 2 y 2 3 (r + u -π(y 1 )) (3.14d) u = -Sgn (Γ(t)) , when Γ(t) ̸ = 0 (3.14e) λ(t f ) = ν 1 ν 2 0 ⊤ , (ν 1 , ν 2 ) ∈ R 2 (3.14f)
Equation (3.14a) states the non-triviality of the adjoint states. Following [START_REF] Bryson | Applied optimal control: optimization, estimation and control[END_REF] and [102, Thm. 7.8.1], since the integral cost, the dynamics and the end-point constraints are time-invariant, the Hamiltonian is constant along the extremals and for such a free time, fixed endpoint problem, this constant is zero:

H(t) ≡ 0, ∀t ∈ [0, t f ]. (3.15)
The optimal pairs (y, u) are called abnormal extremals [START_REF] Agrachev | On Abnormal Extremals for Lagrange Variational Problems[END_REF][START_REF] Montgomery | Abnormal minimizers[END_REF] if λ 0 = 0, and normal extremals if λ 0 ̸ = 0. We now proceed to establish some intermediate results on the adjoint states.

Proposition 4. (λ 1 (t), λ 2 (t)) ̸ = (0, 0) for all t ∈ [0, t f ].
Proof. The linear time-varying dynamics of (λ 1 , λ 2 ) is Lipschitz in (λ 1 , λ 2 ) and continuous in time. Therefore, from the Cauchy-Lipschitz theorem, any maximal solution is unique. Thus, if there is a t 0 s.t. (λ 1 , λ 2 )(t 0 ) = 0, then λ 2 ≡ 0 over [0, t f ], implying λ 3 ≡ 0 from (3.14d) and (3.14f) and then λ 0 = 0 from (3.15), violating (3.14a). Now, remark that the sign of λ 2 and λ2 can be extrapolated from the following second-order equation λ2 = a(t)λ 2 , where a(t

) := - π ′ (y 1 (t)) y 3 (t) > 0. (3.16)
Indeed, the cones R + × R + and R -× R -are both invariant through the dynamics (3.16). This shows that if λ 2 or λ2 is null at some t λ ∈ (0, t f ), then they will both remain in one of these cones after t λ . Further, from Proposition 4 and (3.16), they will actually remain in interior subsets of these cones for times t > t λ . Hence, by an exhaustive enumeration of possible cases we can state the following result.

Proposition 5. (λ 2 , λ2 ) necessarily match one of these conditions, as illustrated in 1. λ 2 and λ2 are never zero on (0, t f ): (b) λ 2 > 0 and λ2 < 0 on (0, t f ),

(a) λ 2 > 0 and λ2 > 0 on (0, t f ),
(c) λ 2 < 0 and λ2 > 0 on (0, t f ), (d) λ 2 < 0 and λ2 < 0 on (0, t f ), 2. There is a unique t λ ∈ (0, t f ) s.t. λ 2 (t λ ) = 0 and λ2 ̸ = 0 on [0, t f ]: (a) Sgn (λ 2 (t)) = -Sgn (t -t λ ) and λ2 < 0, (b) Sgn (λ 2 (t)) = Sgn (t -t λ ) and λ2 > 0, 3. There is a unique t λ ∈ (0, t f ) s.t. λ2 (t λ ) = 0 and λ 2 ̸ = 0 on [0, t f ]: (a) Sgn λ2 (t) = -Sgn (t -t λ ) and λ 2 < 0, (b) Sgn λ2 (t) = Sgn (t -t λ ) and λ 2 > 0.
Note that for scenarios 1a and 1d (resp. scenarios 1b and 1c), either λ 2 or λ2 can be zero at t = 0 (resp. at t = t f ). Also, note that for scenario 2, λ 2 is necessarily non-zero at t = 0 and t = t f since λ2 is of constant sign. The same kind of remark applies to scenario 3 as well. The goal is to state whether these scenarios are consistent with conditions (3.14a)-(3.14f), and if so to what control structure they refer to. Proposition 6. Abnormal extremals are optimal programs of constant thrust.

Proof. For abnormal extremals, λ 0 = 0. Using Equation (3.15) at t = t f yields ν 2 = 0 = λ 2 (t f ). Thus, from Proposition 5, λ 2 has a constant non-zero sign over [0, t f ). Moreover, from (3.14d) and (3.14f), one has

Sgn (λ 3 (t)) = -Sgn (λ 2 (t)) , ∀t ∈ [0, t f ). (3.17) Therefore, for any t ∈ [0, t f ): Sgn(Γ(t)) = Sgn(λ 2 ). Hence, u has a constant value in {-1, +1} over [0, t f ).
This latest proposition shows that abnormal extremals require the initial state y 0 to be on constant thrust trajectories achieving landing, i.e y 0 ∈ Σ max or y 0 ∈ Σ min must hold for these extremals.

From now on, we consider normal extremals only, and, without loss of generality 4 , we consider λ

0 = 1. Let us define b(t) := π(y 1 (t)) y 3 (t)
Note that from Assumption 1 and the sign of y 2 , one can show that a(.) and b(.) are increasing from a study of their derivatives. Also, for all times in [0, t f ], a and b are respectively lower and upper-bounded by

a := - π ′ (ȳ 1 ) m + and b := π(0) m -. (3.18) Let us define γ(t) := λ2 (t) + λ 2 (t)b(t), which satisfies dΓ dt (t) = Γ ′ (t) = γ(t) y 3 (t) (3.19)
Since y 3 is positive, γ carries the sign of Γ ′ . From this point, Γ is the subject of our investigations.

Lemma 2. If γ < 0 over (0, t f ), then Γ is null at most on a single t ∈ [0, t f ]. Lemma 3. Γ < 0 in the left-neighborhood of t f . Proof. Equation (3.15) at t f yields ν 2 = -(r + u(t f ))/ ẏ2 (t f ). Thus, one gets Γ(t f ) = - κy 3 (t f ) + π(0) y 3 (t f ) ẏ2 (t f ) < 0.
The conclusion follows from the continuity of Γ(.).

λ 2 must be non-positive in a neighborhood of t f . Indeed, let us assume that there is a time t ′ s.t. λ 2 is positive on [t ′ , t f ). Note that λ 2 (t f ) may be null. Then, using (3.14d) and (3.14f), λ 3 would necessarily be negative on [t ′ , t f ), leading to Γ(t) > 0 for t in [t ′ , t f ], which contradicts Lemma 3. This eliminates scenarios 1a, 1b, 2b and 3b.

Moreover, note that scenario 1d necessarily corresponds to Min-Max programs, where one arc may be absent, for it satisfies Lemma 2.

Then, the three remaining scenarios, namely 1c, 2a and 3a, require a refined sign study of λ 2 and λ2 . Using differential equations bounding λ 2 , we can establish bounds on γ tight enough to derive valuable sign information. Definition 2. For a constant c > 0 and t 0 ∈ (0, t f ), the C 2 scalar function x c is defined over [0, t f ] as the unique solution of the initial value problem

ẍc = cx c with x c (t 0 ) = λ 2 (t 0 ) and ẋc (t 0 ) = λ2 (t 0 ) which yields x c (t) = λ 2 (t 0 ) cosh √ c(t -t 0 ) + λ2 (t 0 ) √ c sinh √ c(t -t 0 ) .
Inspired from the definition of γ (from (3.19) and before), let us denote

γ c (t) := ẋc (t) + x c (t)b(t) (3.20)
and introduce z λ := (λ 2 , λ2 ) ⊤ and z := (x a , ẋa ) ⊤ s.t.

żλ = F (t, z λ ) :=    0 1 a(t) 0    z λ . (3.21)
The next proofs require the following assumption. Proof. Here, λ2 < 0 and Sgn (λ 2 (.)) = -Sgn (.t λ ), where t λ ∈ (0, t f ). In this proof only, we consider the functions from Definition 2 with t 0 = t λ . It leads to

γ a (t) = λ2 (t λ ) cosh ( √ a(t -t λ )) + b(t) √ a sinh ( √ a(t -t λ )) .
For Lemma 3 and Proposition 7 imply that the sign of Γ changes at most once over [0, t f ] for scenario 2a.

Lemma 4. If λ 2 < 0 over [0, t f ], if γ(t γ ) = 0 for some t γ ∈ (0, t f ) and if Assump- tion 3 holds, then: γ(t) < 0, ∀t > t γ . Proof. By construction λ 2 (t γ ) = -λ2 (t γ )/b(t γ ). Necessarily, λ2 (t γ ) > 0.
In this proof only, we consider the functions from Definition 2 with t 0 = t γ . It yields

γ a (t) = λ2 (t γ ) 1 - b(t) b(t γ ) cosh ( √ a(t -t γ )) + b(t)b(t γ ) -a b(t γ ) √ a sinh √ c(t -t γ )
Since b increases, the factor associated to the cosh term is negative. Also, Assumption 3 yields

b(t)b(t γ ) -a ≤ (b(t) + √ a)(b(t) - √ a) < 0 (3.

Optimality of Min-Max Programs

We shall now discuss under which conditions Min-Max trajectories are always more fuel-optimal than Max-Min-Max trajectories, for some y 0 ∈ F * .

Let us consider a trajectory y starting at y 0 , with thrust structure Max-Min-Max. Denote t 1 its first time of switch (from max to min). The last max arc may be of null duration. Then, for every time t ′ 1 ∈ [0, t 1 ], there is a trajectory with thrust structure Max-Min-Max, with first time of switch t ′ 1 , that lands, which is guaranteed by applying Proposition 3 at t ′ 1 . Below, we derive conditions under which the trajectory having the smallest first time of switch has the highest final mass, showing that the Min-Max trajectory starting from y 0 is fuel-optimal.

The second time of switch, denoted t 2 , and the final time t f are implicitly imposed by t 1 so that the rocket lands. This relation will be given later. For the time being, note that the final mass, denoted y f 3 , satisfies

y 0 3 -y f 3 (t 1 ) = (r + 1)t 1 + (r -1)(t 2 (t 1 ) -t 1 ) + (r + 1)(t f (t 1 ) -t 2 (t 1 )). (3.26)
The first two components of y are collected in µ(y), i.e. µ(y) := (y 1 , y 2 ) ⊤ . The landing condition is simply µ(y(t f )) = 0. Define

L(τ 1 , τ 2 , τ f ) := µ Φ f +g (τ f -τ 2 , Φ f -g (τ 2 -τ 1 , Φ f +g (τ 1 , y 0 ))) .
Then, the landing condition boils down to

L(t 1 , t 2 , t f ) = 0. ( 3
dt 2 dt 1 , dt f dt 1 ⊤ = - ∂L ∂[t 2 , t f ] -1 • ∂L ∂t 1 . (3.28)
To express these derivatives w.r.t. t 1 , intermediate quantities are introduced. The transition matrices M (t f ) and N (t 2 ) are respectively defined as the unique solutions to the matrix initial value problems

Ṁ (t) = ∂(f + g) ∂y (y(t)) • M (t) and M (t 2 ) = I 3 , (3.29) Ṅ (t) = ∂(f -g) ∂y (y(t)) • N (t) and N (t 1 ) = I 3 . (3.30) Let us define R 1 , R 2 , S 1 and S 2 by R 1 R 2 ⊤ := µ (M (t f ) • (f -g)(y(t 2 ))) , (3.31) S 1 S 2 ⊤ := µ (M (t f ) • N (t 2 ) • (f + g)(y(t 1 ))) . (3.32)
Since Assumption 2 holds, the invertibility condition of ∂L ∂[t 2 ,t f ] needed to apply the IFT boils down to R 1 ̸ = 0. Then, one can provide a detailed version of (3.28)

dt 2 dt 1 = 1 - S 1 R 1 , ( 3.33
)

dt f dt 1 = 1 - R 1 S 2 -R 2 S 1 + S 1 ẏ2 (t f ) ẏ2 (t f )R 1 . (3.34)
Using the previous terms with Equation (3.26) yields

dy f 3 dt 1 (t 1 ) = r + 1 ẏ2 (t f )R 1 R 1 (S 2 -ẏ2 (t f )) + S 1 ẏ2 (t f ) r -1 r + 1 -R 2 . (3.35)
The conditions that enable us to state that Min-Max thrust programs are always more fuel-optimal than the Max-Min-Max ones, by allowing us to apply the IFT on L, are thus conveyed by the assumption below Assumption 4. The parameter defined in (3.31) is s.t. R 1 ̸ = 0, and one has

dy f 3 dt 1 (0) < 0 for any y 0 ∈ F * s.t. the rocket lands at y 3 (t f ) ≥ m -.
Note that, since it is formulated for any y 0 ∈ F * , it is sufficient to check these conditions for t 1 = 0 only. Moreover, these conditions can be either checked through (3.35), analytically -if the pressure model is known well enough and tractable -or numerically.

Remark 12. For illustration purposes only, let us check the validity of Assumption 4

when there is no atmosphere. When π ≡ 0, every term from (3.33) and (3.34) can be explicitly written using the fact that r+u y 3 = -ẏ3 y 3 for intermediate integrations, which yields

R 1 = 2 r + 1 1 - y 3 (t f ) y 3 (t 2 ) + log y 3 (t f ) y 3 (t 2 ) , R 2 = -κ + r -1 y 3 (t f ) , S 1 = - 2 r -1 log y 3 (t 2 ) y 3 (t 1 ) , S 2 = -κ + r + 1 y 3 (t f ) . R 1 is negative since y 3 (t 2 ) > y 3 (t f ). Thus, (3.35) becomes dy f 3 dt 1 (t 1 ) = - 4κ ẏ2 (t f )(r -1) 1 R 1 log y 3 (t 2 ) y 3 (t 1 ) < 0. (3.36)
The negativity of this quantity gives the desired conclusion. By continuity, the assumption also holds for scarce atmospheres. Further, an example based on a nonscarce tabulated pressure model is treated in Section 3.4.

Main result

Under Assumptions 1, 2, 3, and 4, if the final mass y f 3 of the landing Min-Max trajectory, starting from a y 0 in F, satisfies y f 3 ≥ m -, then the optimal thrust program of Problem ( 1) is Min-Max, where one arc may be absent. Conversely, if 1) has no solution. Henceforth, it is possible to describe the whole set of feasible initial conditions. Define

y f 3 < m -or if y 0 / ∈ F, then Problem (
Ω y f 3 := Φ -(f -g) τ 1 ,Φ -(f +g) τ 2 , (0, 0, y f 3 ) ⊤ : τ 1 ≥ 0, τ 2 ≥ 0, (r -1)τ 1 + (r + 1)τ 2 ≤ m + -y f 3
which denotes the set of states landing at final mass y f 3 ≤ m + applying a Min-Max control. Minimum (resp. maximum) arcs last for τ 1 (resp. τ 2 ). Thus, the solution set F sol of the initial conditions y 0 s.t. Problem (1) has a solution is

F sol := m -≤y f 3 ≤m + Ω y f 3 (3.37)
The following theorem summarizes this discussion.

y 2 < 0 0 M a s s ( y 3 ) y - 3 y + 3 Altitu de (y 1 ) 0 y 1 > 0 Σ max Σ min Ω(y - 3 )
Figure 3.2: Flight envelope. F sol is delimited by Σ max , Σ min , Ω y - 3 and closed by the constraint y 3 ≤ y + 3 on the last side. For Σ max and Σ min , only the trajectories that land with a mass y 3 (t f ) ≥ m -are represented. The vertical axis conveys the altitude to ease the visualization.

Numerical illustrations

Let us consider the following (normalized) parameters κ = 0.00285 s -1 , r = 4.0, m -= 458.3 s, m + = 520.3 s.

In this example, the engine can be used at 60-100% of its maximum flowrate. Also, κ is taken close to the values of actual reusable launcher engines [START_REF] Onel | Liquid rocket engine performance assessment in the context of small launcher optimisation[END_REF], such as the Merlin (Falcon 9) or the BE-4 (New Glenn). We consider a pressure model describing Earth's atmosphere from tabulated values, satisfying Assumption 1, s.t. π(0) = 6.2 × 10 -1 . Assumption 2 and 3 are satisfied since

a cc = 2.90 × 10 -3 > 0 and b/ √ a = 3.37 × 10 -1 < 1.
F sol is pictured in Figure 3.2 for these values. Assumption 4 is then checked numerically, by computing R 1 and

dy f 3 dt 1 (0)

Comments

The results presented above call for a few comments.

Vertical flight envelope applications

u(y) =            1 if y ∈ Σ max , -1 if y ∈ Σ min ∪ Ω(y - 3 ), v(y) otherwise,
guarantees that a rocket starting inside the vertical flight envelope will land at null vertical speed. This could be used to design safe-set control laws [START_REF] Garcia | A Comprehensive Survey on Safe Reinforcement Learning[END_REF][START_REF] Wabersich | Safe exploration of nonlinear dynamical systems: A predictive safety filter for reinforcement learning[END_REF]. Also, in terms of software design, using the structure where a control "u" supervises another control "v" would help increase the Run-Time Assurance of the rocket G&C system [START_REF] Clark | A Study on Run Time Assurance for Complex Cyber Physical Systems[END_REF][START_REF] Sha | Using simplicity to control complexity[END_REF], but it brings us out-of-the-scope of this thesis.

Non robustness of the Min-Max trajectory

The optimal thrust program being of a bang-bang nature, it is non-robust to some sources of uncertainty, due to the proximity of the ground. For instance, if the switch from the Min to the Max arc is delayed, the rocket necessarily comes out of the vertical flight envelope, which guarantees that there are no thrust programs that allow a proper landing. Given the uncertain and complex environment in which the actual rocket final burn occurs, it can be of high interest to consider more robust thrust programs, that are not too close to the actuators limits, or, in other words, not too close to the vertical flight envelope boundary. dt 1 (0) while dispersing most of the involved parameters, make us conjecture that Assumption 4 could be made much weaker: it is sufficient but not necessary.

Conjecture regarding assumption 4

Recent results on the topic

The results presented in this chapter correspond to the paper [START_REF] Ménou | Fuel-optimal program for atmospheric vertical powered landing[END_REF] published in 2021. Since then, Leparoux et al. published an article on a really close topic [START_REF] Leparoux | Structure of optimal control for planetary landing with control and state constraints[END_REF]. They explored the structure of optimal thrust programs of 3D rocket models for planetary landing. Among others, they showed that the optimal thrust programs were generally of the Max-Min-Max nature, and that it was not sensitive to several model changes. Among the model changes that they considered, they used an atmospheric model that resembles the thrust bias law T = gIspq -S E P (h) that we use in this manuscript, but with a constant pressure term instead.

Other approach for optimal thrust programs The optimal control problem studied in this chapter has a state of dimension 3, which makes the theoretical results from H. J. Sussmann and H. Schättler applicable as well [START_REF] Schättler | The Local Structure of Time-Optimal Trajectories in Dimension Three under Generic Conditions[END_REF][START_REF] Schättler | Geometric Optimal Control: Theory, Methods and Examples[END_REF][START_REF] Sussmann | Time-optimal control in the plane[END_REF][START_REF] Sussmann | Lie Brackets and real analiticity in control theory[END_REF]. They provided a variety of results based on the Lie brackets of the functions f and g involved in Equation (3.2), which can help characterize the optimal solution nature of Problem ( 2), but for a minimum-time criterion only.

Chapter 4

Nominal guidance via Quadratic Programming

Résumé

Le problème général de guidage par descente motorisée PDG (ξ, p), défini dans l'équation (2.22), est un OCP, c'est-à-dire un problème d'optimisation de dimension infinie, en temps libre. Dans ce chapitre, ce problème est réécrit sous la forme d'un problème de dimension finie, aussi simple que possible, afin qu'il puisse être résolu rapidement et de manière fiable en vol. Cette approche est la méthode de guidage par descente motorisée que nous proposons.

Tout d'abord, PDG (ξ, p) est réécrit sous la forme d'un problème d'optimisation non-linéaire (NLP) paramétrique de faible dimension. Cette réécriture nécessite, entre autres, une description des variables d'état basée sur le flot d'équations différentielles ordinaires (ODE), une représentation à dimension finie de la variable de contrôle, et une remise à l'échelle de la variable temporelle pour tenir compte des variations du temps final.

Deuxièmement, la solution de ce dernier NLP est approximée par une expansion directionnelle du premier ordre, en utilisant les résultats classiques de l'analyse de sensibilité des NLP. Étant donné que des variations générales et non infiniment petites des paramètres ξ et p doivent être prises en compte dans l'application, il est nécessaire de traiter les changements dans l'ensemble des contraintes actives. Il s'agit là d'une caractéristique essentielle de l'analyse de sensibilité. Une méthode de calcul basée sur la programmation quadratique (QP) est décrite pour traiter cette question.

Enfin, des commentaires importants concernant l'utilisation offline/online de ce QP sont discutés et illustrés par trois exemples numériques. La méthode de guidage nominal décrite ici est générale et peut être appliquée aux problèmes 2D et 3D.

53

The general Powered Descent Guidance problem PDG (ξ, p), defined in Equation (2.22), is an OCP, i.e. an infinite dimensional optimization problem, in free-final time. In this chapter, this problem is re-written as a finite dimensional problem, as simple as possible, so that it can be solved quickly and reliably in-flight. This approach is our proposed nominal Powered Descent Guidance method.

First, PDG (ξ, p) is re-written as a low dimensional parametric Non-Linear Program (NLP). Among others, this rewriting requires a description of the state variables based on the flow of Ordinary Differential Equations (ODEs), a finitedimensional representation of the control variable, and a re-scaling of the time variable to account for the variations of the free-final time.

Second, the solution of the latter NLP is approximated by a directional first-order expansion, using classic sensitivity analysis results of NLPs. Because general and non-infinitely small variations of the parameters ξ and p must be considered in the application, it is necessary to deal with changes in the set of active constraints. This appears as a critical feature of the sensitivity analysis. A computational method based on Quadratic Programming (QP) is described to handle this.

Finally, important comments regarding the offline/online use of this QP are discussed, and illustrated on three numerical examples. The nominal guidance method described here is general, and can be applied to both the 2D and 3D problems.

This chapter is an updated version of [START_REF] Ménou | Sensitivity Analysis for Powered Descent Guidance: Overcoming degeneracy[END_REF], with enhanced examples and a new application to the 3D rocket model.

Non-Linear Programming formulation for PDG

The goal of this section is to explain how PDG (ξ, p) can be approximated using a NLP with few variables.

As recalled in Chapter 1, OCPs are commonly solved using either direct or indirect methods. On one hand, direct methods first discretize the optimization problem and then solve it using NLP techniques. On the other hand, indirect methods consist in formulating infinite dimensional stationary conditions first, and then discretizing them. The former is often more robust but less accurate than the latter. For both approaches, it is highly difficult to guarantee convergence times for complex nonlinear problems.

We propose to approximate PDG (ξ, p) using a method in-between these two classic approaches. As used in certain direct methods [START_REF] Hargraves | Direct trajectory optimization using nonlinear programming and collocation[END_REF][START_REF] Kraft | On Converting Optimal Control Problems into Nonlinear Programming Problems[END_REF][START_REF] Vlassenbroeck | A chebyshev polynomial method for optimal control with state constraints[END_REF], we discretize the control variable using an interpolation method. However, we do not use a coarse1 discretization scheme to convey the dynamic equation. The latter is here expressed exactly, as the flow of a certain ODE. A point worth special care is that the ODE is defined over a time domain whose endpoint is an unknown of the problem. This representation will help us form a finite dimensional problem denoted NLP (ξ) thereafter. The sensitivity analysis based method used to solve this latter problem will be the matter of the next section. Here is represented a scalar Cubic Spline, described by its values µ 0 , . . . , µ 3 at several time instances, and by its slopes µ 4 and µ 5 at the starting and end-points. The inequality constraints are enforced on the subdivision τ ′ 0 , . . . , τ ′ Nc .

τ Parametric control u µ (τ ) • • • • | τ 0 τ 1 . . . . . . τ Nc τ 0 = 0 τ 1 τ 2 τ 3 = 1 µ 0 µ 1 µ 2 µ 3 µ 4 µ 5 • •

Discretization of the decision variable

Free-final time

First, recall that in PDG (ξ, p), the time-of-flight is an optimization variable implicitly defined by the constraints and the cost. Its change w.r.t. the reference time of flight tf is denoted ∆t f . We scale the time variable t by considering

τ := t tf + ∆t f
The final time is considered as an extra state of null dynamics: ṫf = 0. The augmented state equals x := (x ⊤ , t f ) ⊤ , and satisfies the dynamics

ẋ = f (x, u, η) :=    t f f (x, u, η) 0   
where the variables are here defined for times τ ∈ [0, 1]. The unknown ∆t f is now taken into account as an initial condition, s.t.

x(0) =    x0 + ∆x 0 tf + ∆t f   
To alleviate the writing, the constraint

A f x(1) = b f (p) will be written A f x(1) = b f (p)
as well, where the latter matrix A f is simply the former matrix A f with an extra column of zeros on the right-hand side. Here, b f (p) remains unchanged. Likewise, the notation c(x, u, η, p) will be used to refer to the former constraint c(x, u, η, p) ≤ 0.

Parametric description of the control

We choose to describe the infinite dimensional variable δu using a smooth parametric description. Let us describe δu via a function

(µ, τ ) ∈ R Nµ × [0, 1] → u µ (τ ) ∈ R m
where µ → u µ (τ ) is linear for any fixed τ . There is a matrix valued function M (.) s.t. u µ (τ ) = M (τ )µ. This framework encompasses the use of many interpolation methods, from piecewise constant interpolation to Cubic Splines and Hermite polynomials. Our choice is detailed below.

For N ≥ 2, consider a subdivision of the normalized time interval [0, 1] denoted by the N + 1 time instances τ 0 = 0 < τ 1 < . . . < τ N = 1. We chose to describe u µ as a Cubic Spline defined on the subdivision (τ i ) 0≤i≤N , as represented in Figure 4.1. These Splines are constructed using the classic method from [START_REF] Kraft | On Converting Optimal Control Problems into Nonlinear Programming Problems[END_REF]Sec. 3.3]. Thus, the vector µ embeds the values u k of the corrections at τ k and the slopes u0 and uN s.t.

µ := (u 0 ) ⊤ , (u 1 ) ⊤ , . . . , (u N ) ⊤ , ( u0 ) ⊤ , ( uN ) ⊤ ⊤ ∈ R Nµ . (4.1)
where

N µ = m(N + 3).
The bounds on the engine flow must be discussed differently depending on the engine model considered. On the one hand, for the 2D rocket model, it was decided that the engine flow was directly controlled via q r . Even though this simplifying choice was made for illustration purposes only, choosing a smooth parametric description for u µ has a convenient by-product. Since the controlled flow can be expressed from the real flow s.t. q c = τ q qr + q r , then it will be possible to express exactly the controlled flow, as long as • q r remains differentiable,

• τ q qr + q r remains within [q -, q + ],

• the estimated value qr (0) of the initial real flow is imposed, i.e.: qr (0) + (q r ) µ (0) = qr (0).

On the other hand, the 3D model already takes into account the first-order dynamics on the engine flow. If q -≤ q 0 r ≤ q + and that q -≤ q 0 c ≤ q + , then the real flow will remain within [q -, q + ]. Therefore, imposing the flow bounds via the constraint u -(p) ≤ ū(τ ) + u µ (τ ) ≤ u + (p) is sufficient to guarantee that q c and q r lie in the proper interval.

Moreover, it was previously mentioned that the engine flow dynamics was taken into account but not the orientation dynamics, since the time constant of the latter was significantly much smaller than the former. However, to remain as feasible as possible by the real system, it is of high interest to have a continuous and smooth control law description. Picking a parametric description such as the Cubic Splines ensures the smoothness. Moreover, imposing the initial value of the incidence -or the projected incidences for the 3D model -to be the same as the current estimate ensures the continuity. Therefore, we consider a constraint of the shape ū(0) + u µ (0) = û(0), which also writes

u µ (0) = ∆u init (4.2)
where ∆u init is the gap between the reference control at time 0 and the current control.

Remark 14. When converting OCPs to NLPs, Cubic Splines and other similar discretization methods are often used to approximate the control variable and the state itself [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF][START_REF] Fahroo | Direct Trajectory Optimization by a Chebyshev Pseudospectral Method[END_REF][START_REF] Kraft | On Converting Optimal Control Problems into Nonlinear Programming Problems[END_REF][START_REF] Vlassenbroeck | A chebyshev polynomial method for optimal control with state constraints[END_REF]. However, here, only the variable δu is described by finitely many values, and the state will be described exactly, using Equation (4.7) below. The state is not a Cubic Spline. This choice is motivated by the need to evaluate with a very high accuracy the states, especially at touchdown.

Parametric problem

The control variable of the parametric description of our problem is

z := (µ ⊤ , ∆t f ) ⊤ . (4.3) It is of dimension N z = N µ + 1.
From the previous discussion, let us denote the whole input vector by

ξ := (∆x 0 ⊤ , ∆η ⊤ , ∆u init ⊤ ) ⊤ ∈ R N ξ . (4.4)
Describing the control correction by u µ (.) and imposing the control constraints

u -≤ ū + u µ ≤ u + (4.5)
does not necessarily imply that µ is bounded. Indeed, the slopes of u µ at τ = 0 and τ = 1 are not directly bounded by these constraints, especially if these constraints are enforced on a badly chosen subset of time instances. Likewise, neither the bounds from Equation (4.5) nor the other above-mentioned constraints necessarily imply that ∆t f is bounded. To guarantee that both µ and ∆t f remain bounded, additional constraints are imposed on the decision variable z s.t.

z low ≤ z ≤ z up .
By extending the mixed state-control constraints conveyed by the function c in Equation (2.22f), an approximation of PDG (ξ, p) is min

z∈R Nz J(z, ξ) (4.6a) s.t. ẋ(τ ) = f (x(τ ), ū(τ ) + u µ (τ ), η + ∆η), ∀τ ∈ [0, 1] (4.6b) x(0) =    x 0 + ∆x 0 tf + ∆t f    (4.6c) A f x(1) = b f (p) (4.6d) u µ (0) = ∆u init (4.6e) u -(p) ≤ ū(τ ) + u µ (τ ) ≤ u + (p), ∀τ ∈ [0, 1] (4.6f) c(x(τ )), ū(τ ) + u µ (τ ), η + ∆η, p) ≤ 0, ∀τ ∈ [0, 1] (4.6g) z low ≤ z ≤ z up (4.6h)
Problem (4.6) has a finite-dimensional decision variable, but an infinite number of constraints.

Formulation of the finite dimensional guidance problem

Description of the state as the flow of an ODE

Our goal is to remove x from the description of (4.6). To this purpose, let us introduce a classic notation used for the flow of the ODE following e.g. [START_REF] Bonnans | Course on Optimal Control, Part I: the Pontryagin approach[END_REF][START_REF] Sontag | Mathematical Control Theory[END_REF].

Definition 3 (Flow of f ). Consider the subsets X ⊂ R n (assumed open), U ⊂ R m (assumed compact) and Ω ⊂ R nη . Given a differentiable function f : X ×U ×Ω → R n , vectors (x 0 , η) ∈ X × Ω and a control function 2 u ∈ L ∞ ([0, 1], U)
, the flow of the ODE defined by f is defined using the following Initial Value Problem (IVP)

∀t ∈ [0, 1], x(t) = Φ f t, x 0 , η; u ⇔      x(0) = x 0 ẋ(s) = f (x(s), u(s), η), ∀s ∈ [0, t].
Using this notation, let us describe the extended state

x[τ, z, ξ] ∈ R n+1 as x[τ, z, ξ] := Φ f   τ,    x 0 + ∆x 0 t f + ∆t f    , η + ∆η; ū + u µ    , ∀τ ∈ [0, 1]. (4.7)
The latter notation x[τ, z, ξ] will prove to be handy when writing the inequality and equality constraints (4.8) and (4.9) below. Note also that the hypothesis of Definition 3 guarantee that x[τ, z, ξ] is uniquely defined. For formal results on the existence and uniqueness of Φ f (t, x 0 , η; u), see e.g. [START_REF] Sontag | Mathematical Control Theory[END_REF]Appendix C.3].

Provided that f is continuously differentiable w.r.t. all of its inputs, since µ → u µ (t) is also assumed continuously differentiable for all times t, then x[τ, z, ξ] is continuously differentiable w.r.t. all of its inputs. For detailed properties of Φ, see Appendix A.2.2 recalling some useful classic results.

PDG as a NLP

To alleviate the writing, let us first consider that the constraint parameter p equals zero (the case with p ̸ = 0 will be handled afterwards). The last ingredients that must be discretized in (4.6) are the inequality constraints. Indeed, constraints (4.6f) and (4.6g) are defined on an infinite number of points. We decide to enforce these constraints on a number of N c + 1 times instance

τ ′ 0 = 0 < τ ′ 1 < . . . < τ ′ Nc = 1 s.t. the subdivision (τ ′ i ) 0≤i≤Nc
is an uniform oversampled version of (τ i ) 0≤i≤N , as shown in Figure 4.1. In other words, every interval [τ i , τ i+1 ] is split into several sub-intervals, and the constraints (4.6f) and (4.6g) are enforced at their borders.

Thus, the discretized version of the inequality constraints can be re-written as

h(z, ξ) ≤ 0 where h(z, ξ) :=                           ū(τ ′ 0 ) + u µ (τ ′ 0 ) -u + (0) u -(0) -(ū(τ ′ 0 ) + u µ (τ ′ 0 )) c (x[τ ′ 0 , z, ξ], ū(τ ′ 0 ) + u µ (τ ′ 0 ), η + ∆η, 0) . . . ū(τ ′ Nc ) + u µ (τ ′ Nc ) -u + (0) u -(0) -(ū(τ ′ Nc ) + u µ (τ ′ Nc )) c x[τ ′ Nc , z, ξ], ū(τ ′ Nc ) + u µ (τ ′ Nc ), η + ∆η, 0 z -z up z low -z                           . (4.8)
Moreover, using the same notations, the equality constraints (4.6d) and (4.6e) are conveyed by the condition g(z, ξ) = 0 where

g(z, ξ) :=    A f x[1, z, ξ] -b f (0) u µ (0) -∆u init    . ( 4.9)
Thus, we get an expression of the finite-dimensional constraints approximating PDG (ξ, 0), i.e. for the special case where p = 0. The general case where p ̸ = 0 is then straightforward. Indeed, there is no loss of generality in assuming that the constraint parameter p has a linear influence on the constraints of the original optimization problem, as explained below. For instance, for the 2D rocket model, the final horizontal position in the terminal constraint (4.6d) can be parametrized by a variable ∆z f , s.t.

A f x(1) = b f + 0 0 ∆z f 0 ⊤ .
Likewise, to negotiate the incidence limit, the control bounds are changed by a variable ∆α max s.t.

u --(0, ∆α max ) ⊤ ≤ ū(t) + u µ (t) ≤ u + + (0, ∆α max ) ⊤ .
Therefore, using this assumption of linear influence of p, we assume that there are matrices H p and B p s.t. the original constraints h(z, ξ) ≤ 0 and g(z, ξ) = 0 actually write h(z, ξ) ≤ H p p and g(z, ξ) = B p p when p ̸ = 0. The matrices H p and B p are basically filled by zeros and a few ones.

To sum up all the preceding steps, PDG (ξ, p) is converted into a NLP as follows

• change the time variable from t ∈ [0, t f ] to τ ∈ [0, 1],
• consider a parametric description u µ (.) of the infinite dimensional variable δu,

• describe the dynamic equation and the initial condition constraints through Equation (4.7),

• enforce the inequality constraints on the time instances τ ′ i for i = 0, . . . , N c , instead of enforcing them for all τ ∈ [0, 1],

• enforce bounds on the decision variable z. not usually recommended from a numerical point of view (e.g. [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF]). Any iterative method aiming at solving NLP (ξ, p) (e.g. Successive Quadratic Programming [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF]) requires the evaluation of x[τ, z, ξ] and its derivatives at each iteration, which means solving multiple ODEs at each iteration. However, as it will be detailed below, our goal is to provide an approximation of the solutions of NLP (ξ, p) using a (directional) first-order expansion. Therefore, it is only needed to evaluate these computationally expensive terms once and offline, making the use of x[τ, z, ξ] appropriate in this context.

Sensitivity analysis for degenerate parametric NLP

Without loss of generality, to simplify the exposition, we consider in this section the special case where p = 0. The direct extension to the case p ̸ = 0 will be discussed in Section 4.3.

Let us study the standard problem (4.10), but where the constraint right-hand sides have been simplified s.t.

NLP (ξ) :=

           min z J(z, ξ) s.t. h(z, ξ) ≤ 0, g(z, ξ) = 0.
In fact, NLP (ξ) is a parametric NLP formulated in a standard form [25, Eq. ( 1)]. When it is too hard to compute its exact solution, an alternative approach is to provide a reasonable approximation of it w.r.t. the parameter ξ, near a known solution. This section aims at answering the three following questions. First, if z * denotes the optimal solution of NLP Then, how can we compute this (local) expansion? Finally, how can this be applied to solve NLP (ξ) non-locally (at least approximately)?

An introductory toy example

Let us consider a low-dimensional example that will help us illustrate the challenges of NLP sensitivity. This toy problem has only two variables and one parameter ξ which appears non-linearly.

Basic Example 1.

For a scalar parameter ξ > 0, consider the parametric NLP

min x∈R 2 1 2 (x 2 1 + x 2 2 ) (4.11) s.t. c 1 (x, ξ) = -x 1 -x 2 + ξ 1 + ξ 2 ≤ 0 (4.12) c 2 (x, ξ) = x 1 -x 2 1 + ξ 2 -ξ ≤ 0 (4.13)
This problem posseses a unique solution and can be solved analytically. The solution is illustrated in Figure 4.2. For ξ > 0, the constraint c 1 is active (with associated multiplier λ 1 ) and the optimum is

x * 1 (ξ) = ξ 2 1 + ξ 2 , x * 2 (ξ) = ξ 2 1 + ξ 2 , λ 1 (ξ) = ξ 2 1 + ξ 2 .
For ξ < 0, the constraint c 2 is active (with associated multiplier λ 2 ) and the optimum is

x * 1 (ξ) = ξ 2 + ξ 2 , x * 2 (ξ) = -ξ 1 + ξ 2 2 + ξ 2 , λ 2 (ξ) = - ξ 2 + ξ 2 .
Finally, for ξ = 0, since the point (x 1 , x 2 ) = (0, 0) globally minimizes the cost and satisfies the constraints, it is the optimum. Both constraints are active, and the associated multipliers equal 0.

A first remark is that for all values ξ ̸ = 0, the function ξ → x * (ξ) is continuously differentiable. This is pictured in Figure 4.2, and corresponds to the smooth parts of the black curve. Also, note that the multipliers of the active constraints are positive for ξ ̸ = 0. When ξ = 0, then λ 1 = λ 2 = 0 even though both constraints are active, which is what we will call a degenerate scenario following [START_REF] Jittorntrum | Solution point differentiability without strict complementarity in nonlinear programming[END_REF].

Then, the main point of interest is when ξ = 0, i.e. point A in Figure 4.2. The set of active constraints changes when the sign of ξ changes. The optimal solution x * is not differentiable, but has a left-hand side and right-hand side derivatives. Therefore, the value of x * can be inferred in the neighborhood of ξ = 0 using the following (directional) expansions, for ε ≥ 0:

x * (ε) = x * (0) + dx * dξ (0 + )ε + o(ε) and x * (-ε) = x * (0) - dx * dξ (0 -)ε + o(ε).
The goal of the next sub-section is to highlight the conditions under which x * (.) is at least directionally differentiable. The need to be able to handle degenerate scenarios for PDG -i.e. when some multipliers are zero when the input parameter is zerowill be demonstrated in Section 4.3, and is directly related to local changes in the active set of constraints.

Known results in parametric NLP sensitivity

In this section, we consider the problem NLP (ξ) from a general perspective, i.e. not necessarily expressed by the above-mentioned expressions for h and g. The goal is to present sufficient conditions that enable us to compute an expansion of the solutions of NLP (ξ). Without loss of generality, we seek an expansion in the neighborhood of ξ = 0. First, after recalling necessary conditions for the existence of a minimizer, a classic theorem that gives conditions for the differentiability of the solution is recalled. A discussion on the key aspects of its proof stresses the need of a more general theorem, as its main assumption needs to be relaxed in view of application to PDG. With this view in mind, an alternative result from the literature, based on directional derivatives -a.k.a. Dini derivatives -is presented. Finally, computational aspects are discussed for the evaluation of the solution expansions.

Sensitivity analysis with Strict Complementary Slackness

Introduce the multipliers ν ∈ R n in and λ ∈ R neq and the Lagrangian

L(z, ν, λ, ξ) := J(z, ξ) + ν ⊤ h(z, ξ) + λ ⊤ g(z, ξ). (4.14)
Classically, a tuple (z, ν, λ, ξ) is said to satisfy the Karush-Kuhn-Tucker (KKT) conditions 3 if

L z (z, ν, λ, ξ) = 0, (4.15a) g(z, ξ) = 0, (4.15b) ν ⊤ h(z, ξ) = 0 and ν ≥ 0. (4.15c)
At ξ = 0, the multipliers are denoted ν 0 and λ 0 and are assumed to satisfy the KKT conditions. The compact notation L[0] = L(z 0 , ν 0 , λ 0 , 0) is used if necessary to alleviate the writing. Likewise, h[0] = h(z 0 , 0), h z [0] = h z (z 0 , 0), etc. Lemma 5 (Second-Order Sufficient Conditions (SOSC), [37, Lemma 3.2.1]). If the functions defining Problem NLP (0) are twice continuously differentiable in a neighborhood of z 0 , if there exist multipliers ν 0 ∈ R n in and λ 0 ∈ R neq s.t. the KKT conditions (4.15) hold and, further, if

κ ⊤ L zz (z 0 , ν 0 , λ 0 , 0)κ > 0 (4.16)
for any non-zero vector κ ∈ R Nz that satisfies

[h i ] z (z 0 , 0)κ ≤ 0, ∀i : h i (z 0 , 0) = 0, (4.17a) [h i ] z (z 0 , 0)κ = 0, ∀i : (ν 0 ) i > 0, (4.17b) g z (z 0 , 0)κ = 0. (4.17c)
then, z 0 is a strict local minimizing point of P(0).

Definition 5 (SCS). For a pair (z, ν), Strict Complementary Slackness (SCS) holds

when

ν i > 0, ∀i : h i (z, 0) = 0.
The conditions from Lemma 5 and the Strict Complementary Slackness condition allow one to present the following theorem, adapted from [37, Thm. 3.2.2], which states a well-known NLP sensitivity result.

Theorem 2 (Continuous differentiability, with SCS from [START_REF] Fiacco | Introduction to Sensitivity and Stability Analysis in Nonlinear Programming[END_REF]). Assume that J, h and g are twice continuously differentiable in z and that their gradients w.r.t. z and the constraints are continuously differentiable in ξ in a neighborhood of (z 0 , 0). (c) for ξ near 0, the set of binding inequalities is unchanged, SCS holds, and the binding constraint gradients are linearly independent at z * (ξ).

If (i) [h i ] z (z 0 , 0) (for i s.t. h i (z 0 , 0) = 0)
The proof of this theorem conveys several keys features helping to understand what makes the solution z * continuously differentiable. Precisely, Equation (4.19) below is an highly useful by-product of the proof, providing an explicit formula for the derivative of z * and the associated multipliers.

Remark 16. The SCS assumption, combined with the conditions of Lemma 5, implies that L zz [0] is positive definite on the kernels of both matrices h z a [0] and g z [0]. Denoting by N K the dimension of the intersection of these kernels, and by Q ∈ R Nz×N K the matrix that generates this vector space, the conditions (4.16) and (4.17) can be re-written equivalently as

Q ⊤ L zz [0]Q ≻ 0.
Sketch of proof of Theorem 2, adapted from [START_REF] Fiacco | Introduction to Sensitivity and Stability Analysis in Nonlinear Programming[END_REF] and [START_REF] Büskens | Sensitivity Analysis and Real-Time Optimization of Parametric Nonlinear Programming Problems[END_REF]. The idea of the proof is to use the Implicit Function Theorem (IFT) on a subset of the KKT conditions, which defines functions z * (ξ), ν * (ξ), λ * (ξ) that are then shown to be locally unique minimizers of Problem NLP (ξ) using the conditions of Lemma 5.

At z 0 , the active inequality constraints are denoted with an exponent "a". Their cardinal is denoted n a . For instance, h a denotes the rows of h that are active at z 0 , and ν 0 a are the corresponding multipliers. Without loss of generality, we assume that the active inequality constraints are the first n a components of h and ν.

The main function of interest in this proof is4

K(z, ν a , λ, ξ) :=       L z (z, ν a , λ, ξ) ν a • h a (z, ξ) g(z, ξ)       =              J z (z, ξ) + (ν a ) ⊤ h z a (z, ξ) + λ ⊤ g z (z, ξ) ν 1 h 1 (z, ξ) . . . ν na h na (z, ξ) g(z, ξ)             
and is often referred to as the Kuhn-Tucker matrix [START_REF] Büskens | Sensitivity Analysis and Real-Time Optimization of Parametric Nonlinear Programming Problems[END_REF]. The conditions of Lemma 5 imply that the matrix

∂K ∂[z, ν a , λ] (z 0 , ν 0 a , λ 0 , 0) =              L zz [h 1 ] z ⊤ . . . [h na ] z ⊤ g z ⊤ (ν 0 ) 1 [h 1 ] z h 1 0 0 . . . . . . . . . (ν 0 ) na [h na ] z 0 h na g z 0 . . . 0              (4.18)
is invertible. Note that all the elements on the right-hand side of the latter matrix are evaluated on (z 0 , ν 0 a , λ 0 , 0), but this has been omitted to alleviate the writing.

Thus, the IFT applies on the condition K(z, ν a , λ, ξ) = 0 at (z 0 , ν 0 a , λ 0 , 0), which guarantees the existence of a differentiable function

y : ξ → (z * (ξ), ν * (ξ), λ * (ξ)) s.t. z * (0) = z 0 , ν * (0) = ν 0 , λ * (0) = λ 0 and K(z * (ξ), (ν * ) a (ξ), λ * (ξ), ξ) = 0 in the neighborhood of ξ = 0.
Finally, the strict inequalities of SCS enables us to show that the set of active constraints does not change in the neighborhood of ξ = 0, which implies that y(ξ) locally satisfies the sufficient conditions of Lemma 5, hence the theorem.

This proof has a useful by-product. Indeed, the construction of the triplet (z * (ξ), ν * (ξ), λ * (ξ)) via the use of the IFT directly provides the derivatives of these variables at ξ = 0

        dz * dξ (0) d(ν * ) a dξ (0) dλ * dξ (0)         = - ∂K ∂[z, ν a , λ] (z 0 , ν 0 a , λ 0 , 0) -1 ∂K ∂ξ (z 0 , ν 0 a , λ 0 , 0) . (4.19)
However, as we have seen with Toy Example 1, SCS does not necessarily hold everywhere, which becomes a problem in the above-mentioned arguments. Indeed, if one of the multipliers of the active constraints becomes zero, then the corresponding line in Equation (4.18) becomes zero, leading to a singular Kuhn-tucker matrix. In this scenario, the IFT does not apply anymore.

Sensitivity analysis without Strict Complementary Slackness

To overcome the above-mentioned difficulty, another set of assumptions is needed.

Definition 6 (Strong SOSC).

There exists a scalar a > 0 s.t.

κ ⊤ L zz (z 0 , ν 0 , λ 0 , 0)κ ≥ a∥κ∥ 2 (4.20)
for any κ ∈ R Nz s.t.

g z (z 0 , 0)κ = 0, [h j ] z (z 0 , 0)κ = 0, ∀j s.t. h j (z 0 , 0) = 0 and (ν 0 ) j > 0.
The SOSC property is weaker than the SCS property. As shown by Jittorntrum [START_REF] Jittorntrum | Solution point differentiability without strict complementarity in nonlinear programming[END_REF], relaxing the SCS using Strong SOSC instead eventually leads to directional differentiability properties. To formulate this, let us introduce the following. For a function γ : R p → R q , when it exists, the upper Dini derivative of γ at x in the direction d is a vector of R q and is denoted by

D + d γ(x) := lim ε↓0 γ(x + εd) -γ(x) ε ∈ R q .
For any given direction ξ, consider the directional problem NLP (εξ), where ε ≥ 0 is a scalar.

Theorem 3 (Directional differentiability, without SCS). Let us assume that J, h

and g are twice continuously differentiable in a neighborhood of z 0 . If (i) [h i ] z (z 0 , 0) (for i s.t. h i (z 0 , 0) = 0) and [g j ] z (z 0 , 0) (all j) are linearly independent, (ii) Strong SOSC holds at z 0 , then there exists a unique continuous function ε → (z * (εξ), ν * (εξ), λ * (εξ)) that locally minimizes NLP (εξ), for ε ≥ 0, s.t. z * (0) = z 0 , ν * (0) = ν 0 and λ * (0) = λ 0 . Furthermore, its right-hand derivative at ε = 0 exists, i.e. the upper Dini derivatives

D ξ + z * (0), D ξ + ν * (0) and D ξ + λ * (0) exist.
For a proof of this theorem, see the discussions after the Theorems 3 and 4 in [START_REF] Jittorntrum | Solution point differentiability without strict complementarity in nonlinear programming[END_REF]. A direct consequence of Theorem 3 is that for any ξ, for ε ≥ 0, the following expansion holds:

z * (εξ) = z 0 + D + ξ z * (0)ε + o(ε). (4.21)
The same way Theorem 2 provided an expression for the derivatives of the solution tuple, the proof of Theorem 3 by Jittorntrum provides a way to compute the latter expansion, which will be used later in Equation (4.25). While Theorem 2 requires to inverse a linear system as shown in Equation (4.19), Theorem (3) needs to solve a QP, as indicated in the Proposition below.

Proposition 10 (Adapted from [START_REF] Jittorntrum | Solution point differentiability without strict complementarity in nonlinear programming[END_REF]Eq. 24]). The vector D + ξ z * (0) in Equation (4.21) is the optimal solution of

min ∆z∈R Nz 1 2 ∆z ⊤ L zz [0]∆z + ξ ⊤ L ξz [0]∆z s.t. h a z [0]∆z + h a ξ [0]ξ ≤ 0, g z [0]∆z + g ξ [0]ξ = 0,
whose uniqueness is guaranteed by Assumption (ii) in Theorem 3.

Proposition 10 provides a way to compute a local, directional, expansion of the optimal solution z * (ξ) using QP. Qualitatively, it boils down to linearizing the active inequality constraints, forgetting the inactive inequality constraints, linearizing the equality constraints, and taking into account the second-order approximation of the cost. Therefore, aside from the non-linear aspect of NLP (ξ), the only missing features of Problem NLP (ξ) in Proposition 10 are the inactive inequality constraints. This result can be improved and made more useful for practical applications.

To alleviate the writing, and without loss of generality, we assume for the rest of this section that z 0 = 0. If instead of considering only the active inequality constraints we consider all inequality constraints, i.e. that we solve QP (ξ) := min then we have a more useful tool. On one hand, the result of Proposition 10 is preserved. Indeed, denote by i any (strictly) inactive inequality constraint at ξ = 0. Consider any arbitrary value of ξ. Since ε → z * (εξ) is continuous on a certain right-neighborhood of ε = 0 by Theorem 3, then h i (z * (εξ), εξ) < 0 holds in this neighborhood. By reducing this neighborhood if necessary, the condition

z∈R Nz 1 2 z ⊤ L zz [0]z + ξ ⊤ L ξz [0]z (4.22a) s.t. h[0] + h z [0]z + h ξ [0]ξ ≤ 0 (4.22b) g z [0]z + g ξ [0]ξ = 0 (4.22c)
h i [0] + [h i ] z [0]z * (εξ) + [h i ] ξ [0]εξ < 0
holds as well in the vicinity of ε = 0 and will not affect the local property of Proposition 10. On the other hand, the linearized conditions (4.22b) provide a linear approximation of the inequality constraints, giving a non-local estimation of changes in the active set 5 .

As an intermediate summary, if a parametric NLP described by NLP (ξ) satisfies Strong SOSC and some other mild conditions described in Theorem 3, its solution point is locally uniquely defined, and admits a directional first-order expansion, that can be computed by solving QP (ξ) detailed in (4.22). This latter Quadratic Program consists simply in modifying NLP (ξ) by linearizing the constraints in both the parameter and the decision variable, and by taking a second-order expansion of the cost. Locally, the approximations provided by QP (ξ) are very accurate (in the sense of (4.21)), and can handle active set changes that depend on the direction ξ. Globally, QP (ξ) can keep providing an approximation of the solution, that will be as good as the approximations made in (4.22b) and (4.22c) are. From a high-level point of view, this Chapter describes how to convert an infinite dimensional problem into a finite one, which in turns is solved using parametric sensitivity analysis. However, the converse approach which consists in performing parametric sensitivity analysis on the infinite dimensional problem directly is an alternative that has been explored in the literature. See e.g. [START_REF] Deshpande | Directional Input Adaptation in Parametric Optimal Control Problems[END_REF], that describes how directional differentiability can be used on the stationary conditions of OCPs with multiple types of constraints.

Remark 17. A direct application of QP (ξ) to the Basic Example 1 is shown in

Fast nominal guidance method

Up to here, we have presented a method to re-formulate the original problem PDG (ξ, p) into the finite-dimensional problem NLP (ξ, p). For p = 0, the solutions of NLP (ξ, 0) are approximated by solving QP (ξ).

Let us extend this sensitivity-based method to the case where p ̸ = 0, and summarize how the newly formed problem QP (ξ, p) is intended to be used in practice to provide nominal guidance.

An offline/online approach for nominal guidance

First, note that the role of p in the definition of NLP (ξ, p) is equivalent to the one of ξ, in the sense that both are given parameters. Thus, by introducing convenient intermediate matrices and vectors, we can extend the QP from Equation (4.22) to the general case under the form

QP (ξ, p) :=            min z∈R Nz 1 2 z ⊤ P z + ξ ⊤ Qz s.t. Gz ≤ h 0 + H ξ ξ + H p p Az = b 0 + B ξ ξ + B p p (4.23)
where H p and B p are the same as the ones used in the definition of the constraints of NLP (ξ, p) and where P := J zz (0, 0) Q := J ξz (0, 0)

G := h z (0, 0) A := g z (0, 0) h 0 := -h(0, 0) H ξ := -h ξ (0, 0) b 0 := -g(0, 0) B ξ := -g ξ (0, 0)
It is this expression of QP (ξ, p) that will be used from now on.

Remark 20. Note that the reason why p does not appear in the cost of QP (ξ, p) is that the cost J in Equation (4.10a) does not depend on p, but only on z and ξ. Remark 21 (Shortcuts). In the remainder of this thesis, the dependency of the linear part of the cost and the constraint right-hand side on ξ in QP (ξ) may be omitted to alleviate the writing, by using the vectors q, h and b s.t.

q := J ξz (0, 0) ⊤ ξ, h := h 0 + H ξ ξ, and b := b 0 + B ξ ξ.
To maximize computational efficiency, the constant matrices of QP (ξ, p) are computed offline, for a reference trajectory, using the transition matrices formulas from the Appendix. As described in Section 4.1, such a reference trajectory (x, ū, η, tf ) must satisfy the constraints of NLP (ξ, p), and its design depends critically on the mission goals. The reference may be defined as the solution of a more complex optimization problem, out-of-the-scope of the thesis, solved using indirect methods [START_REF] Bryson | Applied optimal control: optimization, estimation and control[END_REF]. Its computation can take anything from several seconds to hours, depending on the desire accuracy, the optimization solver or the computer used for it. On the other hand, computing the transition matrices, even to a good accuracy, only needs a few seconds, even with a non efficiency-optimized code.

Guidance law

As introduced in (4.3), the optimal values µ * and ∆t * f returned by QP (ξ, p) via z * enable us to describe the guidance control law as a continuous time function. Indeed, interpreting ū and u µ as functions defined on [0, 1], the guidance law becomes

u * (t) = ū t tf + ∆t * f + u µ * t tf + ∆t * f , ∀t ∈ [0, tf + ∆t * f ]. (4.24)

Directional first-order estimate of waypoints

Let us denote by z nlp (ξ, p) the hypothetical solution of NLP (ξ, p) and by z qp (ξ, p) the solution of QP (ξ, p). Then, from Section 4.2, the following directional expansion holds

z nlp (εξ, εp) = z qp (εξ, εp) + o(ε).
Moreover, using the definition of the augmented state x[τ, z, ξ] from Equation (4.7), for any τ ∈ [0, 1] we have

x[τ, z, ξ] =    x[τ ] tf    + ∂ x ∂z [τ, 0, 0]z + ∂ x ∂ξ [τ, 0, 0]ξ + o(∥(z, ξ)∥).
Let us introduce xlin as x lin (τ, ξ, p)

xlin (τ, ξ, p) :=    x[τ ] tf    + ∂ x ∂z [τ, 0, 0]z qp (ξ, p) + ∂ x ∂ξ [τ, 0, 0]ξ ∈ R n+1 . ( 4
t lin f (τ, ξ, p)    . Horizontal position z x 0 = x(τ 0 ) • • • • • •
Reference trajectory, x using ū during tf with η

x(τ 1 )

x(τ 2 )

x(τ 3 )

x(τ 4 )

h(0)- • Altitude h Solving QP (ξ, p) at x 0 + ∆x 0 gives µ * and ∆t * f (Here, p = 0) ∆x 0 • • • • • x 1 x 2 x 3
x 4 Way points, expressed using ∆x 0 , ∆η, µ * and ∆t * f . Note that t lin f does not actually depend on τ by construction of the augmented state, and Equation (4.25) applied to t lin f simply yields t lin f (ξ, p) = tf + ∆t * f (ξ, p). In practice, the solution z * of QP (ξ, p) can be used either to form the guidance law from Equation (4.24), or by computing waypoints x k using the state approximation provided by x lin s.t.

x k := x lin (τ k , ξ, p) (4.26)
where (τ k ) k is a prescribed subdivision of [0, 1] of normalized time instance. The values τ k are equivalent to the non-normalized time instances t k s.t.

t k = τ k . tf + ∆t * f (ξ, p) .

Numerical examples

Let us discuss the three following examples, each highlighting a distinct feature of the solutions of QP (ξ, p). Since the case where p ̸ = 0 is the matter of Chapter 5, these examples focus exclusively on nominal descents, i.e. where p = 0. For all the examples, the cost J that we consider is defined directly by giving its matrices P and Q. We consider that Q = 0, and that the matrix P is a diagonal matrix of positive weights, which favors early corrections (i.e. the weight associated to u(τ 0 ) is smaller than the one of u(τ N )).

The reference trajectory chosen for each example is represented by a black line. Moreover, for the examples using the 3D model, the reference trajectory is assumed to lie in the plane (e z , e h ) (though corrections implying out of plane trajectories are discussed in the last example).

The data have been normalized for all the examples.

Effectiveness of calculated guidance

The first example aims at showing that in some neighborhood of the reference parameters, the expansion (4.25) provides accurate corrections enforcing the terminal constraints, even when using the corrections in open-loop only. It is directly illustrated on the 3D rocket model. Due to the terminal constraints, the six states corresponding to the positions and speeds are expected to be null at the final time (except

v h (t f ) = -ε f v h ).
If the reference guidance law (ū, tf ) is applied directly to a scenario where ξ ̸ = 0 -i.e. no corrections are applied -then strong constraint violation errors will appear. However, if the change in parameter is corrected using u * from (4.24) (i.e. using QP (ξ, p)), the terminal constraints are supposed to be approximately satisfied up to the first-order. These two behaviors are well observed in the sub-figure (a) of Figure 4.5, which represents the terminal horizontal position z(t f ). The other terminal constraints (on y, h, v z , v y and v h ) have voluntarily been omitted, as their correction curve is much flatter than for z. In other words, it means that even though this terminal constraint component has the worst open-loop correction curve, it still demonstrates that the first-order corrections brought by QP (ξ, p) work well in a non-trivial neighborhood of ξ = 0 in practice.

Changes in the active set

The second example aims at showing that the optimal solutions of QP (ξ, p) are indeed only Dini-differentiable and not smooth, even in a standard landing scenario. Let us use the 2D rocket model, and consider that ξ only varies in horizontal position ∆z 0 = z 0 -z0 .

The directional-derivatives of α at the middle point t 2 in the two directions ∆z 0 = 1 and ∆z 0 = -1 differ, as shown in Figure 4. 6-(d). This behavior has a real-world interpretation. When ∆z 0 > 0, using more incidence on α(t 2 ) is not a possible option as the constraint is already active and becomes strictly active (the associated multiplier becomes positive when ∆z 0 > 0). However, when ∆z 0 < 0, lowering α(t 2 ) is possible, allowing the presented corrections.

Non-local constraint satisfaction

This third example aims at showing the behavior of QP (ξ, p) for large values of ξ. Let us consider the 3D model, with a reference trajectory contained in the plane (e z , e h ).

As illustrated in Figure 4.7, we are here interested in the behavior of QP (ξ, p) for values of ξ covering a grid in (∆z 0 , ∆y 0 ), the change in initial horizontal positions. Considering a such choice of inputs allows us to stress that QP (ξ, p) is able to compute out-of-plane trajectories.

As pointed out in the sub-Figures 4.7 If one considers even larger values of ξ in this example, it comes a point where the constraints of QP (ξ, p) are infeasible, at least as long as we keep p = 0. Hence the need for an algorithm that computes the proper value for p according to a hierarchy of objectives, which will be the matter of the next chapter.

Chapter 5

Emergency guidance via Linear and Quadratic Programming

Résumé

Le problème QP (ξ, 0) a été le sujet principal jusqu'à présent. À un moment donné, si ξ est trop grand, alors QP (ξ, 0) peut être infaisable. Par exemple, si la position horizontale initiale est trop éloignée du site d'atterrissage, il n'est pas possible de concevoir une trajectoire qui satisfasse toutes les contraintes en même temps : ces dernières sont incompatibles. La question centrale de cette thèse est donc la suivante : que faire dans cette situation ?

Lorsque QP (ξ, p) est infaisable à p = 0 et pour une valeur donnée de ξ, le problème de l'atterrissage doit être modifié pour retrouver la faisabilité, dans une certaine mesure. Cela peut se faire en relâchant les contraintes, c'est-à-dire en modifiant p. Une stratégie de relaxation doit être définie.

Tout d'abord, nous identifions les principaux paramètres négociables qui peuvent être relâcher, dans la Section 5.1, et nous soulignons leur importance relative. Cette liste ordonnée découle des connaissances préalables et de la compréhension commune des ingénieurs chargés de la réussite de la mission.

Ensuite, à partir de cette liste ordonnée, une suite de problèmes d'optimisation visant à minimiser l'amplitude des paramètres négociables est présentée dans la Section 5.2 et nommée HEGO pour Optimisation Hiérarchique pour le Guidage d'Urgence. L'utilisation de cette suite produit une trajectoire de guidage, tout en imposant une hiérarchie prescrite entre les paramètres négociables sélectionnés, au sens de l'ordre extended colexicographic order présenté au Chapitre 1. La formulation des problèmes sous-jacents repose uniquement sur la programmation linéaire et quadratique. De plus, des garanties théoriques sur le caractère bien posé et la régularité de la fonction qui à ξ associe la trajectoire optimale sont fournies. Comme il est d'un grand intérêt pour les applications pratiques, où la régularité des méthodes numériques est toujours souhaitable, il est démontré que cette dernière fonction est Lipschitz-continue, évitant ainsi les sauts dans la trajectoire de guidage pour des entrées arbitrairement proches. D'autres propriétés peuvent être établies, et la monotonie directionnelle de l'amplitude du paramètre de négociation est examinée dans un cas particulier.

L'Algorithme HEGO est généralisé dans la Section 5.6. Cette généralisation permet de distinguer ce qui est lié à la méthode de guidage nominal sousjacente de ce qui est propre aux problèmes de négociation eux-mêmes. On montre notamment comment HEGO pourrait être utilisé avec un autre choix de méthode de guidage nominal (i.e. au lieu de QP (ξ, p)), parmi les méthodes directes pour les OCP.

Enfin, des résultats numériques sont présentés dans la Section 5.7. Ils mettent en avant divers aspects de la méthode de guidage d'urgence, à la fois sur les modèles de fusée 2D et 3D. Plusieurs ensembles de paramètres négociables sont utilisés pour illustrer les différentes stratégies d'urgence accessibles par la méthodologie HEGO ainsi proposée.

The problem QP (ξ, 0) has been the main topic so far. At some point, if ξ is too large, then QP (ξ, 0) may be infeasible. For instance, if the initial horizontal position is too far from the landing site, then it is not possible to design a trajectory that satisfies all the constraints at the same time: they are inconsistent. Hence, the central question of this thesis is: what should one do in this case?

When QP (ξ, p) is infeasible at p = 0 for a given value of ξ, the landing problem has to be modified to recover feasibility to some extent. This can be done by revising the constraints, i.e. by modifying p. Some revision strategy has to be defined.

First, we identify the main negotiable parameters that can be relaxed, in Section 5.1, and emphasize their relative importance. This ordered list stems from prior knowledge and common understanding between engineers in charge of the mission success.

Then, from this ordered list, a sequence of optimization problems aiming at minimizing the magnitude of the negotiable parameters is introduced in Section 5.2 and denoted HEGO for Hierarchical Emergency Guidance Optimization. Using the sequence produces a guidance trajectory, while enforcing a prescribed hierarchy between the selected negotiable parameters, in the sense of the extended colexicographic order introduced in Chapter 1. Formulating the underlying problems relies only on Linear and Quadratic Programming. Moreover, theoretical guarantees on the well-posedness and the smoothness of the mapping from ξ to the optimal trajectory are provided. As it is of high interest for a practical application where smoothness of numerical calculation procedures is always desirable, the latter map is shown to be Lipschitz-continuous, preventing jumps in the guidance trajectory for arbitrarily similar inputs. Some further properties can be established, e.g. the directional monotonicity of the negotiation parameter magnitude is discussed for a special case.

The HEGO algorithm is generalized in Section 5.6. It distinguishes what is linked to the underlying nominal guidance method, from what is fundamental to the negotiation problems themselves. In details, it is shown how HEGO could be used with another choice of nominal guidance method (i.e. instead of QP (ξ, p)), among direct methods for OCPs.

Finally, numerical results are presented in Section 5.7. They highlight various aspects of the emergency guidance method, on both the 2D and the 3D rocket models. Several sets of negotiable parameters are employed, to illustrate the various emergency policies achievable by the proposed HEGO methodology. [START_REF] Ménou | Nominal And Emergency Rocket Landing Guidance Using Quadratic Programming[END_REF]. However, the last sections, including the theoretical proofs and the numerical simulations, are new (unpublished) elements.

The first sections of this chapter are a detailed version of

Negotiable parameter choices

For a given input ξ, when landing has been declared infeasible by QP (ξ, p), it is necessary to loosen some of the constraints. First, the parameters that can be negotiated are listed and it is shown how they modify the constraints of QP (ξ, p). Then, a model describing their relative importance is proposed.

Negotiated constraints

The parameters that can be negotiated are the ones describing the goals of the landing. The physics-based equations of motion are not negotiable. However, the location of the landing site is, at least partially, negotiable. Indeed, if the landing site is located in a wide and flat area, it is of interest to allow touchdowns in a neighborhood of the ideal landing site1 . Some of the other parameters defining the constraints can be partially loosened. For example, the incidence limit should be seen more as a safety constraint -and a way to limit long-term fatigue of the rocket -and could be slightly widened if necessary, whereas the engine flow limitations are non-negotiable mechanical constraints.

The negotiable parameters are already conveniently conveyed by p, the variable introduced in Chapter 2. The purpose of this chapter is to discuss extensively how these negotiable parameters are mathematically modeled, what are the possible choices for its components, and how its value is chosen.

Let us recall that p has a linear influence on the Right-Hand Side (RHS) of the constraints of QP (ξ, p), s.t. the nominal constraints

     Gz ≤ h 0 + H ξ ξ Az = b 0 + B ξ ξ
are transformed into the negotiated constraints by the action of p

     Gz ≤ h 0 + H ξ ξ + H p p Az = b 0 + B ξ ξ + B p p (5.1)
For any such p, the matrices H p and B p are basically filled with zeros and a few ones, which makes them sparse.

These negotiated constraints come with an extra condition: it is assumed that all negotiable parameters are negotiable within prescribed limits, i.e. that p is bounded If necessary, this cube-like constraint could be modeled as a bounded polytope, described by an inequality Dp ≤ d for some pair (D, d), without changing any of the reasoning below.

For example, as it will be discussed in more details in Example 1 at the end of this chapter, a possible choice for the variable p is the pair (∆α max , ∆z f ) ⊤ for the 2D rocket model. From a more exhaustive point of view, the list of all parameters that can or cannot be loosened in practice is presented in Table 5.1.

On the relative importance of the parameters

Adapting Orwell's Animal Farm quote, all parameters are important but some are more important than others (in (5.1)). It is necessary to enforce a hierarchy of importance between the negotiation parameters. Thus, let us separate the negotiation parameters in R different sub-parameters p (j) of ranked (increasing) importance, as

p = (p (1) ) ⊤ . . . (p (R) ) ⊤ ⊤ ∈ R nneg .
The higher the index j, the more critical p (j) is. The dimension of p (j) is noted n j .

Mathematically speaking, the negotiable parameters are compared using the 1norm of their sub-parameters, by comparing their most critical sub-parameters first. As first introduced in (1.1) in Chapter 1, a vector p a is said to be more negotiated than another vector p b , which is denoted p a ⪰ e p b , if and only if

∥p (R) a ∥ 1 > ∥p (R) b ∥ 1 or ∥p (R) a ∥ 1 = ∥p (R) b ∥ 1 and ∥p (R-1) a ∥ 1 > ∥p (R-1) b ∥ 1 , or . . . or ∥p (R) a ∥ 1 = ∥p (R)
b ∥ 1 and . . . and ∥p (1) a ∥ 1 > ∥p

(1) b ∥ 1 , or ∥p (R) a ∥ 1 = ∥p (R)
b ∥ 1 and . . . and ∥p (1) a ∥ 1 = ∥p (5.

3)

The relation ⪰ e is an extended colexicographic order, that we will eventually refer to as the emergency order in this thesis2 . (1) , p (2) ) ⊤ , where p (1) ∈ R and p (2) ∈ R 2 . Then (Least critical) (1, 2, 1) ⪯ e (1, 3, 1 ) ⪯ e (2, -4, 0 ) ⪯ e (-3, 0, 5) (Most critical).

Basic Example 2. Let us consider negotiable parameters

∥p (2) ∥ 1 = 4 ∥p (2) ∥ 1 = 4
For illustration purposes, let us consider again p = (∆α max , ∆z f ) for the 2D rocket model. Considering that it is less critical to sacrifice a few degrees of incidence limit than to land outside the desired landing site yields p (1) = ∆α max and p (2) = ∆z f . Remark 22. To some extent, recovering feasibility has been tackled from different perspectives. In a seminal paper by Blackmore et al. [START_REF] Blackmore | Minimum-Landing-Error Powered-Descent Guidance for Mars Landing Using Convex Optimization[END_REF], the problem of finding the actual landing site that minimizes the distance to the desired landing site is described in details, for Mars landing (i.e. without atmosphere), using Successive Convexification. Translated into the above-mentioned taxonomy, one can write that they have two different negotiable parameters that were negotiated at the same time, which were the two final horizontal positions of their 3D lander model. Their taxonomy also implies that R = 1 with ours, meaning that they do not use any notion of hierarchy.

A hierarchical negotiation

Now that it has been discussed what the levers available to modify the constraints of QP (ξ, p) are -using Equation (5.1) -there comes the follow-up question: how does one negotiate these parameters?

When the input ξ makes the nominal constraints infeasible -i.e. when nominal landing is not feasible anymore -the goal is to find the smallest change in the negotiable parameters that recovers feasibility. There are two salient expectations regarding the method that computes these negotiable parameters and the associated trajectory that will be called emergency trajectory.

On the one hand, this method must pick the smallest negotiable parameter as possible, in the sense of the emergency order. This aims at using as little negotiable parameters as possible while enforcing the relative importance of the parameters.

On the other hand, the map ξ → z * giving the optimal emergency trajectory should be as smooth as possible (continuity is a minimum), in order to prevent jumps in the trajectory between arbitrarily close values of ξ. Such jumps could be very detrimental and cause serious issues to the control algorithms (out-of-the-scope of the thesis).

With these objectives in mind, we introduce the method HEGO (i.e. Algorithm 2 below), composed of a finite sequence of negotiation problems labeled LP j and a refinement problem labeled Refine, which aims at fulfilling the two latter goals. HEGO is then tested on a low-dimensional example, helping to understand some of the methodology design choices. Finally, several aspects of this algorithm are put into perspective. The theoretical proofs showing that this algorithm behaves as expected will be discussed in Section 5.3.

Algorithmic principle of HEGO

From a high-level perspective, finding the smallest negotiation parameters that satisfy the constraints could take the form of a single optimization problem s.t.

z * , p * ←-min z,p
Penalty(p) s.t. Negotiated constraints for (z, p), from Eq. (5.1).

However, there are several limitations to this strategy. First, nothing guarantees that neither z * nor p * are unique when coming out of such a procedure. Moreover, since the cost of the latter problem differs from the cost of the nominal guidance problem QP (ξ), it is highly likely that the map ξ → z * encounters a discontinuity when the emergency is triggered, i.e. when going from the most inner set to the intermediate set of Figure 5.1. Finally, and most importantly, the hierarchy is ignored with this kind of problem description. At best, the parameters can be non-homogeneously weighted to reflect their importance, but not their strict ranking.

Instead, we propose to successively minimize the magnitude of the negotiable parameters while enforcing the existence of at least one feasible trajectory at each step. This is a penalty-free approach. As required by the emergency order, this minimization procedure will focus sequentially on the sub-parameters, starting by the last one (i.e. p (R) ), down to the first one (i.e. p (1) ). This means that the most critical negotiable parameters are minimized first. At each step, only the proper sub-parameter p (j) is minimized, s.t. the result be p (j) = 0 if it is not necessary to use this parameter to recover feasibility. Also, to enforce the desired hierarchy, a kind of memory effect is needed so that each step takes into account the results of the previous steps, by preserving the negotiation levels. This will be the role of condition (5.4e) below. As will appear, Linear Programming (LP) will play a key role to implement these negotiation steps. Quantitatively, let us first introduce the negotiation problems LP j as

LP j := min z,p ∥p (j) ∥ 1 (5.4a) s.t. Gz ≤ h 0 + H ξ ξ + H p p (5.4b) Az = b 0 + B ξ ξ + B p p (5.4c) p low ≤ p ≤ p up (5.4d) ∥p (i) ∥ 1 = P * i , i = j + 1, . . . , R (5.4e)
where P * i denotes the optimal value of LP i . To make this problem definition wellposed, note that the constraint (5.4e) does not exist when j = R. Moreover, note that the inputs of each LP j are ξ and P * i for i = j + 1, . . . , R. To alleviate the writing, these inputs are omitted wherever the context is clear enough. Finally, it is important to highlight the fact that z and p are both optimization variables in LP j , even though only a few coefficients of p are involved in the cost (5.4a).

The role of LP j is to minimize a cost on the j -th negotiable sub-parameters, while making sure that there are still feasible trajectories z, and that the already determined levels of relaxation of the previous negotiation problems are unchanged in the process. The reason why constraint (5.4e) must be satisfied instead of a constraint of the type "p (j) = p (j) * " is that LP j does not necessarily have a unique solution. Imposing ∥p (j) ∥ 1 = P * i encompasses all the solutions of LP j and makes sure that the level of negotiation reached at step i remains satisfied in the follow-up negotiations. Descending from j = R to j = 1 and imposing this latter constraint guarantees that the emergency order is enforced.

Thus, solving successively LP j for j = R, . . . , 1 (decreasing indices) provides a way to recover feasibility, while hierarchically minimizing what needs to be sacrificed from the original guidance problem, by building the sequence P * 1 , . . . , P * R starting from the end. Among others, a noteworthy property of this process is that if landing is feasible without any negotiation, then solving LP j will return ∥p (j) ∥ 1 = 0 at each step, implying that the overall vector p equals zero. Once this negotiation sequence has been computed, there may be many possible values for p, and for z as well. It is thus necessary to pick the best trajectory among these ones, by solving Refine := min z,p

1 2 z ⊤ P z + ξ ⊤ Qz (5.5a) s.t. Gz ≤ h 0 + H ξ ξ + H p p (5.5b) Az = b 0 + B ξ ξ + B p p (5.5c) p low ≤ p ≤ p up (5.5d
)

∥p (i) ∥ 1 = P * i , i = 1, . . . , R (5.5e)
Like LP j , the problem Refine takes as inputs ξ and P * i for i = 1, . . . , R. These two optimization problems are the basic bricks of the central algorithm of this thesis, defined below.

Algorithm 2 Hierarchical Emergency Guidance Optimization (HEGO)

Require: Difference w.r.t. reference values: ξ = (∆x 0 , ∆η, ∆u init ).

for j = R, . . . , 1 (decreasing indices) do P * j ← min LP j ξ, P * j+1 , . . . , P * R // From definition (5.4). end for z * ← argmin Refine (ξ, P * 1 , . . . , P * R ) // From definition (5.5). return z * For the argmin operation, the value of p is voluntarily ignored, since it is not needed nor unique.

It is important to observe that this guidance procedure provides, depending on the situation, either a nominal guidance or an emergency guidance, as stated in the following proposition. Proposition 11. If ξ is s.t. the constraints of the nominal guidance problem QP (ξ, 0) are feasible, then HEGO and QP (ξ, 0) return the same value z * .

Proof. If the constraints of QP (ξ, 0) are feasible, then there exists at least one pair (z, p) with p = 0 that satisfies the negotiated constraints (5.1). Thus, all the problems LP j will necessarily return P * j = 0, and constraint (5.5e) will subsequently impose p = 0, making Refine and QP (ξ, 0) coincide. Uniqueness of the solutions yield the conclusion.

Qualitatively, the nominal guidance method based solely on QP (ξ, 0) can provide landing trajectories on the set A from Figure 5.1 whereas HEGO is able to do it on the sets A and B. HEGO therefore provides landing guidance for a wider range of input values ξ. It also meets the requirements regarding the negotiable parameter minimization and the hierarchy enforcement. These features make it an "universal" guidance algorithm.

The analysis of the map ξ → z * defined by this algorithm is presented later in Section 5.3. Among others, Proposition 12 will show the well posedness of the algorithm, and Theorem 7 proves the Lipschitz-continuity of ξ → z * (ξ) for HEGO. Before getting into these technical details, let us discuss the following example.

An illustrative toy example

Consider a low-dimensional example that resembles the landing problem, illustrating how HEGO works and why the negotiation problem hierarchy matters. Although this example could seem greatly over-simplified at first glance, its similarities with the actual landing problem are noteworthy, especially when comparing the curves from Figure 5.3 (left) and Figure 5.5 (d) reporting the negociation of parameters.

The problem takes as input a scalar ξ ≥ 0, and aims at minimizing the norm of

z = (z 1 , z 2 ) ⊤ ∈ R 2 , under some constraints min z, p z 2 1 + z 2 2 s.t. (Ineq 1 ) z 1 ≥ 0 (Ineq 2 ) z 2 ≥ 0 (Eq) z 2 = 1 -ξ -z 1
This problem is represented in Figure 5.2 (0). By analogy, let us imagine that z 1 conveys the incidence, z 2 the engine flow and ξ the initial horizontal position error. Therefore,

• (Ineq 1 ) conveys the incidence bound, that may be negotiated by a variable p 1 , s.t. z 1 ≥ -p 1 .

• (Ineq 2 ) conveys the mechanical limits of the engine flow, which are nonnegotiable.

• (Eq) represents the terminal condition on the horizontal position and is directly influenced by ξ. It may be negotiated by a parameter p 2 if necessary, s.t.

z 2 = 1 -ξ -z 1 + 2p 2 .
Moreover, note that the negotiable parameters p = (p 1 , p 2 ) ⊤ ∈ R 2 can be negotiated up to some extent, so we decide to bound them between 0 and 1. More precisely, imposing 0 ≤ p 2 ≤ 1 means that we allow the rocket to eventually land farther of the landing site, but on one side only and up to a certain limit. Finally, as in the actual landing problem, we consider that negotiating p 2 is more critical than p 1 (i.e. R = 2). Let us now review the possible scenarios depending on the input value.

Nominal scenario

When ξ remains low, i.e. 0 ≤ ξ ≤ 1, then there is no need for parameter negotiation, because the problem is feasible when p 1 = p 2 = 0. Running Algorithm 2 will give null values for the negotiation penalties. In this case, the optimal solution is

z * 1 (ξ) = z * 2 (ξ) = (1 -ξ)/2.
This scenario is represented in Figure 5.2 (A).

First negotiation scenario

When 1 < ξ ≤ 2, the initial constraints are not compatible anymore and must be negotiated, whence the need for Algorithm 2. The first step of the latter is min

z, p |p 2 | s.t. z 1 ≥ -p 1 z 2 ≥ 0 z 2 = 1 -ξ -z 1 + 2p 2 0 ≤ p 1 ≤ 1 0 ≤ p 2 ≤ 1
and will result in a null optimal negotiation, i.e. P * 2 (ξ) = 0, since it is possible to recover feasibility without using p 2 . The second step will be the problem min

z, p |p 1 | s.t. z 1 ≥ -p 1 z 2 ≥ 0 z 2 = 1 -ξ -z 1 + 2p 2 0 ≤ p 1 ≤ 1 0 ≤ p 2 ≤ 1 0 = |p 2 |
where "0 = |p 2 |" is here to enforce the result from the former negotiation problem. The latter will return the optimal negotiation P * 1 (ξ) = ξ -1. Thus, the optimization variable z can be re-optimized over the newly negotiated set -which is actually a singleton. In the end, it yields

z * 1 (ξ) = 1 -ξ and z * 2 (ξ) = 0.
This scenario is represented in Figure 5.2 (B).

Mild negotiation scenario

Then, let us consider a scenario that requires even more negotiations, i.e. when 2 < ξ ≤ 4. In this case, the first step of Algorithm 2 gives a non-zero optimal negotiation s.t. P * 2 (ξ) = ξ-2 2 , and the second step gives the result P * 1 (ξ) = 1. Qualitatively, this must be interpreted as: "the parameter p 2 must be modified just enough so that there is a value of p 1 that provides a non-empty feasible set for z". In this scenario, z * 1 (ξ) = -1 and z * 2 (ξ) = 0, which are plotted in Figure 5.2 (C).

Infeasible scenario

Finally, for ξ > 4, there are no possible solutions for the first negotiation problem. Therefore, Algorithm 2 does not return anything, because there are no solutions to this over-constrained problem. Note that the limit p 2 ≤ 1 is the ultimate constraint that makes the negotiation problem infeasible.

The importance of the parameters hierarchy

The values of the variables z and p are represented in Figure 5.3 (left), w.r.t. the input ξ.

The chosen parameters hierarchy is crucial. Indeed, if instead of considering that "p 2 is more critical than p 1 " it was considered that the whole vector (p 1 , p 2 ) could be negotiated at once, the results would have been completely different, as shown in Figure 5.3 (right). In the latter case, there would be only one negotiation problem min z, p

|p 1 | + |p 2 | s.t. z 1 ≥ -p 1 , z 2 ≥ 0, z 2 = 1 -ξ -z 1 + 2p 2 , 0 ≤ p 1 ≤ 1, 0 ≤ p 2 ≤ 1.
Doing this would imply that the variable p 2 would be used to recover feasibility before p 1 .

Noteworthy remarks

Remark 23 (Linear Programs). Using standard material from the literature, such as [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF]Example 1.13], decomposing p = ρ + -ρ -with ρ + , ρ -≥ 0 makes it possible to solve LP j using Linear Programming, whence the name LP j . Indeed, this decomposition yields the convenient re-writing ∥p∥ 1 = 1 ⊤ (ρ + + ρ -). See also the re-writing method in the Lipschitz-continuity proofs below.

Remark 24 (In the literature). From a very general mathematical programming point of view, recovering feasibility in Linear Programming has been discussed extensively by Chinneck [START_REF] Chinneck | Feasibility and Infeasibility in Optimization: Algorithms and Computational Methods[END_REF] for instance. Problem (5.4) builds upon right-hand side constraint "alteration" methods, by exploiting the available levers conveyed through the parameter p, the matrices H p and B p , and the need to enforce the parameter hierarchy.

Remark 25 (Hierarchy does not impact feasibility). Though the hierarchy notion is important, it does not change the set of inputs s.t. there exists at least one value of the negotiable parameters that make the constraints (5.1) feasible. Mathematically, it means that the set of ξ ∈ R N ξ s.t.

∃(z, p) ∈ R Nz × R nneg :            Gz ≤ h 0 + H ξ ξ + H p p Az = b 0 + B ξ ξ + B p p p low ≤ p ≤ p up
does not depend on ⪰ e by construction. The latter order will only determine which values of p will be used for a given ξ. This fact can be observed in Figure 5.3, by remarking that the regions (D) and (C ′ ) correspond to the same sets for ξ.

0 z 1 z 2 | 1 | 1 (Ineq 1 ) (Ineq 2 ) (Eq) • • Feasible set for z when p1 = p2 = 0 and ξ = 0. 0 z 1 z 2 | -1 | 1 ξ • • 0 z 1 z 2 | -1 | 1 ξ p 1 • 0 z 1 z 2 | -1 | 1 ξ p 2 p 1 • 0 A B C
| | | | 1 2 3 4 | | 0 1 | -1 z 1 z 2 ξ Optimal penalties | | | | 1 2 3 4 | | 0 1 | -1 P * 1 P * 2 A B C D With the hierarchy. ξ Optimal values of z | | | | 1 2 3 4 | | 0 1 | -1 z 1 z 2 ξ Optimal penalties | | | | 1 2 3 4 | | 0 1 | -1 P * 1 P * 2 A' B' C'
Without the hierarchy. where both variables p 1 and p 2 are negotiated at the same time. The penalties differ from the previous case, and so do the optimal values z 1 and z 2 .

Remark 26 (Limits of the hierarchy). HEGO guarantees that the most critical parameters are minimized first. However, since the negotiable parameters can have a very different influence on the constraints, it is still possible to have scenarios where a critical parameter is non-zero, but a less critical parameter is zero. For instance, let us consider the following trivial example where the negotiated constraints take the form

z 1 = ξ 1 + p 1 , z 2 = ξ 2 + p 2 , -1 ≤ z 1 ≤ 1, -1 ≤ z 2 ≤ 1.
Since there are no links between the variables indexed by 1 and the ones indexed by 2, even imposing that p 2 be more critical than p 1 does not guarantee that p 1 will ever be used when the magnitude of ξ 2 increases.

To summarize, the least critical negotiable parameters will be used in lieu of the most critical ones only if the fundamental nature of the problem makes it possible. 3 after the thesis defense, trying to minimize the value of a parameter with respect to a lexicographic or co-lexicographic order is similar to some methods used in multiobjective optimization. For instance, see [START_REF] Cococcioni | Lexicographic multiobjective linear programming using grossone methodology: Theory and algorithm[END_REF] for more details on the topic.

Remark 27. (Post-defense note) As noted by one of the jury members

Smoothness of the HEGO algorithm

For the reasons mentioned previously, the fact that the outputs of Algorithm 2 -i.e. HEGO-do not change too fast when its inputs vary is of high interest.

The goal of this section is to prove that the map ξ → z * defined by Algorithm 2 is globally Lipschitz on its definition domain. Proving this property relates to QP sensitivity analysis w.r.t. the constraints RHS and the linear part of the cost. When the cost is defined with a positive definite matrix, this property holds and is a wellknown result [START_REF] Mangasarian | Lipschitz Continuity of Solutions of Linear Inequalities, Programs and Complementarity Problems[END_REF]. However, this has to be adapted to our framework, where only a part of the optimization variable is unique (i.e. z is unique, but p is not necessarily unique).

To alleviate the writing, and without loss of generality, we make the assumption that the cost of Refine only has a quadratic term in z, i.e. that Q = 0 in Equation (5.5a). Indeed, extending the constraints RHS sensitivity results to perturbations in the linear part of the cost is obtained using well-known dualization methods of QPs [START_REF] Boyd | Convex Optimization[END_REF][START_REF] Gauvin | Formulae for the Sensitivity Analysis of Linear Programming Problems[END_REF][START_REF] Mangasarian | Lipschitz Continuity of Solutions of Linear Inequalities, Programs and Complementarity Problems[END_REF].

First, a re-writing of the problems LP j and Refine is introduced, to show that the result z * from Algorithm 2 is well-defined and unique. Then, we show that the optimal negotiation maps giving P * i and the optimal solution maps giving z * are Lipschitz continuous functions of the RHS of their constraints. This is proved using a series of well-known results and by adapting a theorem by Mangasarian & Shiau [START_REF] Mangasarian | Lipschitz Continuity of Solutions of Linear Inequalities, Programs and Complementarity Problems[END_REF] to our framework. Finally, we conclude on the properties of Algorithm 2.

Problem re-writing

The constraints of LP j and Refine can be rewritten into a unified, standard and linear framework. Denote by C j these constraints

C j :=                  Gz ≤ h 0 + H ξ ξ + H p p, Az = b 0 + B ξ ξ + B p p, p low ≤ p ≤ p up , ∥p (i) ∥ 1 = P * i , i = j + 1, .
. . , R. For j = 1, . . . , R, the constraints C j convey those of LP j , and C 0 convey those of Refine. Slack variables, denoted s G , s up , s low , can be used to transform the inequalities of C j into the constraints

Gz + s G = h 0 + H ξ ξ + H p p, p + s up = p up , p -s low = p low , s G , s up , s low ≥ 0 Let us define the variable x = (z + , z -, ρ + , ρ -, s G , s up , s low ) where z = z + -z -, (5.6a
)

p = ρ + -ρ -, (5.6b)
x ≥ 0.

(

and introduce the matrices Āj and rj as

Āj :=              G -G -H p H p I dim s G O O A -A -B p B p O O O O O I nneg -I nneg O I nneg O O O I nneg -I nneg O O -I nneg O O I j+1 I j+1 O O O              , rj :=              h 0 + H ξ ξ b 0 + B ξ ξ p up p low P R j+1             
where I j+1 and P R j+1 denote P R j+1 := P * R . . .

P * j+1 ⊤ ∈ R R-j and I j+1 =          O . . . O 1 ⊤ n R 1 ⊤ n R-1 O . . . . . . . . . O 1 ⊤ n j+1 . . . O          ∈ R (R-j)×nneg .
Note that for j = R, the last line of ĀR and rR is absent. Also, for j = 0, the first column of zeros is absent for I 1 . Recall for the rest of the proofs that rj is the vector bearing all the dependency of the problems LP j and Refine on ξ. Finally, note A j and r j the following matrices

A j :=    Āj -Āj    and r j :=    rj -r j    .
Using the new non-negative variable x, and the notations above, the constraints C j can be re-written under the two equivalent forms

     Āj x = rj x ≥ 0 or      A j x ≥ r j x ≥ 0 (5.7)
Further, for j = 1, . . . , R -1, by construction of the constraints C j , the following property holds

{x | Āj x = rj , x ≥ 0} = {x | Āj+1 x = rj+1 , x ≥ 0, c ⊤ j+1 x = P * j+1 } (5.8)
and will be used later to prove Proposition 12. Regarding the cost, define the vector c j of the same dimension as x s.t.

c j :=      (. . . , 0, 1 ⊤ n j , O 1×(nneg-n j ) , 1 ⊤ n j , 0, . . .) ⊤ if j = 1, . . . , R O if j = 0
where the terms 1 ⊤ n j correspond to the position of the j th negotiable variable p (j) , here conveyed by the corresponding parts of ρ + and ρ -from Equation (5.6). Moreover, let us define

D :=       P -P -P P O O O       ∈ R dim x×dim x .
It is straightforward to verify that, for all P positive definite, the matrix D is positive semidefinite. Thus, the optimization problems LP j and Refine can be re-written in a more standard form than their original definition, respectively (5.4) and (5.5). For any j = 1, . . . , R, LP j can be described as a Linear Program of the form

(Primal LP) j            min x c ⊤ j x s.t. Āj x = rj x ≥ 0 (5.9)
and Refine as a Quadratic Program which is

(Primal QP)            min x 1 2 x ⊤ Dx + c ⊤ 0 x s.t. Ā0 x = r0 , x ≥ 0.
           min x 1 2 x ⊤ Dx + c ⊤ 0 x s.t. A 0 x ≥ r 0 , x ≥ 0. (5.11)
Certainly, in view of numerical implementation, there are more memory-efficient ways to translate Refine into these formats, but the formulations above are handy in the theoretical proof below. Also, note that we keep the same definition for P * j as in Problem (5.4), i.e. P * j denotes the optimal value of the cost of Problem (5.9) (its well-posedness will be established in Proposition 12).

Uniqueness of the optimal trajectory

Recalling that c 0 = 0, let us introduce the dual 4 of Problem (5.10) s.t.

(Dual QP)            max x,µ,λ -1 2 x ⊤ Dx -λ ⊤ r0 s.t. Dx -µ + Ā0 ⊤ λ = 0, µ ≥ 0.
(5.12)

Using the Strong Duality Theorem 5 we get the following property. Proof, adapted from Lemma 2.1 in [START_REF] Berkelaar | Sensitivity analysis in (degenerate) quadratic programming[END_REF]. Since the tuples from the statement are optimal and since there is no duality gap, then the primal and dual cost are equal s.t.

1 2 (x * ) ⊤ Dx * = - 1 2 (x) ⊤ Dx -λ ⊤ r0
Since the tuples are optimal solutions, complementary slackness holds, i.e. µ ⊤ x = 0 and µ ⊤ x * = 0 (see e.g. Proposition 5.1.5 in [START_REF] Bertsekas | Nonlinear Programming[END_REF]). Applying the equality constraints from (5.12) at x, and multiplying it by (x * ) ⊤ to the left gives Proof. Let us prove by induction for j = R, . . . , 1 (with decreasing indices), that "Problem (5.9) has a finite optimal value at index j, and constraints (5.7) are feasible at index j -1". Beforehand, note that the cost function of all Problems 5.9 is lowerbounded by 0. For j = R, since χ(ξ) is assumed non-empty, and since the cost is lower-bounded by 0, then it has a finite optimal value, denoted P * R . Moreover, this minimum is attained, as guaranteed by Lemma 12 in the Appendix, at a point denoted x R . Thanks to the Equation (5.8), we get that x R is a feasible point for the constraints (5.7) at index R -1.

(x * ) ⊤ Dx = (x * ) ⊤ µ -(x * ) ⊤ Ā0 ⊤ λ. Thus 1 2 (x * -x) ⊤ D(x * -x) = -λ ⊤ r0 + (x * ) ⊤ Ā0 ⊤ λ = λ ⊤ ( Ā0 x * -r0 ) =
To prove the induction, let us assume that Problem (5.9) has a finite optimal value at index j, and constraints (5.7) are feasible at index j -1, for some j ≥ 2. Using the relation from Equation (5.8), and by induction, the set x | Āj-1 x = rj-1 , x ≥ 0 is non-empty. Thus, since its cost is lower-bounded by 0, Problem (5.9) at index j -1 has a finite optimal value P * j-1 , also attained at some point denoted x j-1 . The latter being a feasible point for constraints (5.7) at index j -2, one concludes the induction proof.

This shows that the optimal penalties P * 1 , . . . , P * R are well defined, and that the constraints (5.7) are feasible at index 0. Consequently, Problem (5.10) also has a minimum. By Lemma 7, and using the expression of D, we get that any solution of Problem (5.10) has a unique value for z + and z -, showing that z * = z +z -exists and is unique. This concludes the proof. Among others, this last proposition points out that as long as the first-to-becomputed negotiation problem is feasible -i.e. LP R -then Algorithm 2 will necessarily terminate. This is summarized as follows, using the set Λ := {ξ ′ | χ(ξ ′ ) ̸ = ∅}. Proof. Using the definition of rR , there is a constant vector b and a constant matrix B s.t. rR = b + Bξ, where b and B directly depend on h 0 , H ξ , b 0 and B ξ appearing in the formulation of QP (ξ, p) defined in (4.23). Let us assume that ξ 1 and ξ 2 are s.t. χ(ξ 1 ) ̸ = ∅ and χ(ξ 2 ) ̸ = ∅. For x 1 ∈ χ(ξ 1 ) and x 2 ∈ χ(ξ 2 ), the condition

(1 -t)x 1 + tx 2 ∈ χ (1 -t)ξ 1 + tξ 2 , ∀t ∈ [0, 1],
holds due to the linearity of (5.7) w.r.t. x and ξ, which shows the desired result.

It is noteworthy that Λ is not necessarily bounded. Indeed, if the intersection of the kernels of matrices B ξ and H ξ (defined in Problem (4.23)) is wider than the singleton {0}, then the set Λ is even guaranteed to be unbounded.

Regularity w.r.t. the right-hand side of the constraints

The proof that z * is a Lipschitz-continuous map of its inputs can be split in two steps. First, we need to show that the optimal penalties P * j are Lipschitz-continuous maps of their inputs, and then that z * is also a Lipschitz-continuous map of ξ and P * i for i = 1, . . . , R. The former result is rather straightforward and will be dealt with in Lemma 8. However, the latter result requires a bit more attention, and will be detailed in Theorem 6.

Remark 28. The results used below are expressed with both the Euclidean and the maximum norms, respecting their original formulation, whereas the main result is given in an homogeneous form -i.e. using only the Euclidean norm -in Theorem 7.

RHS regularity of the negotiation maps

The following theorem is adapted from [START_REF] Mangasarian | Lipschitz Continuity of Solutions of Linear Inequalities, Programs and Complementarity Problems[END_REF], and applies to the standard LP min

∥x 1 * -x 2 * ∥ ∞ ≤ L    b 1 -b 2 d 1 -d 2    2 (5.14)
A direct corollary is the following.

-p * 2 | ≤ K    b 1 -b 2 d 1 -d 2    2
Lemma 8. For any j = 1, . . . , R, the maps ξ ∈ Λ → P * j are Lipschitz continuous.

Proof. In order for these maps to be well-defined, recall that each value P * j depends on ξ and the preceding values P * j+1 , . . . , P * R , (except for P * R that only depends on ξ). Thus, for all j = 1, . . . , R, we are precisely interested in the maps Γ j : Γ j : Λ → R + ξ → P * j (ξ, Γ j+1 (ξ, . . .), . . . , Γ R (ξ)). where Γ R (ξ) = P * R (ξ). Consider j = R. By composition, since the RHS of Problem (5.9) is affine dependent on ξ, and since the optimal value of Problem (5.9) is Lipschitz-continuous w.r.t. its RHS (due to Corollary 2), then Γ R is Lipschitz continuous. Then, by induction, let us assume that at each step j = 1, . . . , R, the previous functions Γ j+1 , . . . , Γ R are Lispchitz continuous. Using the same composition argument, we obtain the desired result.

RHS regularity of the Linear Complementary Problem

To show that the map (ξ, P * 1 , . . . , P * R ) → z * is Lipschitz continuous, we proceed in three steps. First, we recall that polytopes satisfy a Lipschitz continuity-like property (Theorem 5) w.r.t. their RHS. Then, we show how the former map is related to a certain Linear Complementary Problem (LCP, in Lemma 9). Finally, the Lipschitz-continuity of the uniquely defined components of the LCP is established (Theorem 6). respectively. There exists a constant µ, that depends only on A and C, s.t. for each x 1 ∈ F 1 , there exists an x 2 ∈ F 2 closest to x 1 in the ∞-norm s.t.

∥x 1 -x 2 ∥ ∞ ≤ µ    b 1 -b 2 d 1 -d 2    2
Proof of the following lemma can be found in [START_REF] Murty | Linear and Combinatorial Programming[END_REF]Sec. 16.4.4].

Lemma 9 (Linear Complementary Problem, [START_REF] Mangasarian | Lipschitz Continuity of Solutions of Linear Inequalities, Programs and Complementarity Problems[END_REF][START_REF] Murty | Linear and Combinatorial Programming[END_REF]). Assume that D is positive semidefinite. Then, x is a solution of Problem (5.11) if and only if there exists a vector η s.t.

x :=    x η    (5.15)
is a solution to the following LCP M x + q ≥ 0, x ≥ 0, (M x + q) ⊤ x = 0 (5.16) where

M :=    D -A 0 ⊤ A 0 O    and q :=    c 0 -r 0   .
Given a subset J ⊂ {1, . . . , dim x}, any solution of the following linear system6 M j x + q j ≥ 0, xj = 0, j ∈ J, (5.17a)

M j x + q j = 0, xj ≥ 0, j / ∈ J, (5.17b) is a solution to the LCP (5.16) for (M, q). For such sets J, denote Q(J) the set of all q vectors for which (5.17) has a solution 7 .

Lemma 10 (Active set partitions [START_REF] Mangasarian | Lipschitz Continuity of Solutions of Linear Inequalities, Programs and Complementarity Problems[END_REF]Lemma 3.1]). Let q 1 and q 2 be two distinct vectors and let q(t) := (1-t)q 1 +tq 2 for every t ∈ [0, 1]. Assume that the LCP (5.16) for (M, q(t)) is solvable for t ∈ [0, 1]. Then, there exists a partition 0 = t 0 < . . . < t N = 1 s.t. for i = 1, . . . , N

q(t i-1 ) ∈ Q(J i ), q(t i ) ∈ Q(J i ),
for some J i ⊂ {1, . . . , n}.

The constructive proof of this result can be found in [59, p.592]. Its main purpose is to provide a characterization of the active set changes along [0, 1].

When D is positive definite, then Lemma 10 is instrumental to show that the solutions of the LCP (5.16) are Lipschitz w.r.t. q [START_REF] Mangasarian | Lipschitz Continuity of Solutions of Linear Inequalities, Programs and Complementarity Problems[END_REF]. Among others, the latter reasoning relies on the fact that the positive definiteness of D uniquely defines the solution. However, in our framework, D is only positive semidefinite, and the result can not be applied directly. Instead of using much more abstract results of the literature (see e.g. [START_REF] Aubin | Lipschitz Behavior of Solutions to Convex Minimization Problems[END_REF][START_REF] Lee | Continuity of the Solution Map in Quadratic Programs under Linear Perturbations[END_REF]), we decided to adapt the proofs of Mangasarian & Shiau [START_REF] Mangasarian | Lipschitz Continuity of Solutions of Linear Inequalities, Programs and Complementarity Problems[END_REF] and to focus only on the part of the optimization variable that is uniquely defined.

More precisely, the variable x defined in (5.15) equals x = (z + , z -, ρ + , ρ -, s G , s up , s low , η),

where the first part "z + , z -" is uniquely defined, as pointed out in Proposition 12. Therefore, we consider that x can be decomposed in two parts x u and x m , s.t. x u is uniquely defined and x = (x u , x m ). Let us introduce the following theorem, which is a generalized version of [START_REF] Mangasarian | Lipschitz Continuity of Solutions of Linear Inequalities, Programs and Complementarity Problems[END_REF]Thm. 3.2].

Theorem 6 (Lipschitz continuity of the uniquely defined components of the LCP).

Let q 1 and q 2 be points s.t. the LCP (5.16) for (M, q(t)) with q(t) = (1t)q 1 + tq 2 has a solution x(t) = (x u (t), x m (t)) s.t. x u (t) is unique, for every t ∈ [0, 1]. Then, there exists a constant σ > 0, depending only on M , s.t. any solutions x1 = (x 1 u , x 1 m ) and x2 = (x 2 u , x 2 m ) of (5.16) with respective vectors q 1 and q 2 satisfy

∥x 1 u -x 2 u ∥ ∞ ≤ σ∥q 1 -q 2 ∥ 2 (5.

18)

Proof. There exists a subdivision 0 = t 0 < t 1 < . . . < t N = 1 satisfying the properties of Lemma 10. For i = 0, ..., N , let x(t i ) = (x u (t i ), x m (t i )) be a solution of (5.16) for (M, q(t i )). Note that the x u (t i ) are unique, but the x m (t i ) are not necessarily unique.

For any set of indices J ⊂ {1, . . . , dim x} and any matrix A, denote by A J (resp. A J ) the matrix composed of the rows of A whose indices are in J (resp. in {1, . . . , dim x}\J). With this notation, for any 1 ≤ i ≤ N and for any t ∈ [t i-1 , t i ], the LCP (5.16) reduces to the linear problem

   M J i I Ji    x(t) +    q(t) J i O Ji    ≥ 0 and    M Ji I J i    x(t) +    q(t) Ji O J i    = 0 (5.19)
by construction of J i and Q(J i ). Then, according to Theorem 5, stating the Lipschitz-continuity of feasible points of linear constraints, there exists a solution ŷ(t i-1 ) := (y u (t i-1 ), y m (t i-1 )) of (5.16) for (M, q(t i-1 )), i.e. a feasible point of (5.19) at t = t i-1 , s.t.

∥x(t i ) -ŷ(t i-1 )∥ ∞ ≤ µ i ∥q(t i ) -q(t i-1 )∥ 2
for some µ i > 0. Let us define σ := max{µ i | 1 ≤ i ≤ N }. Since the first part of x(t i-1 ) and ŷ(t i-1 ) is uniquely defined, i.e. x u (t i ) = y u (t i ), the following inequality holds

∥x 1 u -x 2 u ∥ ∞ ≤ N i=1 ∥x u (t i ) -x u (t i-1 )∥ ∞ = N i=1 ∥x u (t i ) -y u (t i-1 )∥ ∞ ≤ N i=1 ∥x(t i ) -ŷ(t i-1 )∥ ∞ ≤ N i=1 µ i ∥q(t i ) -q(t i-1 )∥ 2 ≤ σ N i=1 ∥(t i -t i-1 )(q 1 -q 2 )∥ 2 = σ∥q 1 -q 2 ∥ 2 .
Hence the desired result.

Remark 29. The constants of Theorem 4 and Theorem 6 are actually defined via a constructive approach, which is detailed in [START_REF] Mangasarian | Lipschitz Continuity of Solutions of Linear Inequalities, Programs and Complementarity Problems[END_REF].

Conclusion on the Lipschitz-continuity of HEGO

∥(z * ) 1 -(z * ) 2 ∥ 2 ≤ L∥ξ 1 -ξ 2 ∥ 2 .
Proof. This theorem links all the previously stated results. The optimal function z * is well defined on Λ due to Proposition 12. Its value is the minimum of Problem (5.11), whose RHS is affinely dependent on the optimal penalties P * 1 , . . . , P * R . The latter are Lipchitz-continuous maps of ξ, as shown in Lemma 8. Moreover, the solutions of Problem (5.11) are solutions of the LCP (5.16), as recalled in Lemma 9. The uniquely defined components of the solutions of the LCP -i.e. the variables z + and z -in x -are Lipschitz-continuous functions of the vector q from (5.16), in the sense of Theorem 6. Also, the latter vector q is affinely dependent on the vector r 0 (see Lemma 9), which is affinely dependent on ξ and the optimal penalties P * 1 , . . . , P * R . Thus, by composition, z + and z -are Lipschitz-continuous maps of ξ, in the sense of Equation (5.18). The desired result, i.e. with the 2-norm on both sides of the Lipschitz inequality, stems from the equivalence of the norms in finite dimension.

A direct corollary of Theorem 7 is the following. Corollary 3. Let z * (ξ) denote the value returned by Algorithm 2 with input ξ. Then, the optimal value function ξ → J(z * (ξ)) is Lipschitz-continuous on Λ.

Monotonicity of the optimal negotiations

This section aims at giving a mathematical meaning to the sentence "the farther from the reference trajectory, the higher the negotiation", for the special case R = 1, i.e. when there is no hierarchy involved. Getting farther from the reference is modeled directionally by considering the map t → tξ where t ≥ 0 and ξ is an arbitrary input direction. Among others, we aim at showing that the map t → P * 1 (tξ) is non-decreasing, for non-negative values of t, which is what we call the negotiation monotonicity.

Formally, when R = 1, there is only a single negotiation problem, that writes

P * ←-min z,p ∥p∥ 1 (5.20a) Gz ≤ h 0 + H ξ ξ + H p p (5.20b) Az = b 0 + B ξ ξ + B p p (5.20c) p low ≤ p ≤ p up (5.20d)
where P * is the optimal value. Using the same kind of re-writing technique as in Section 5.3.1, one can show that there are matrices M and Q, and a vector r, s.t. Problem (5.20) is equivalent to the following LP in its standard primal form

V * P (r) ←-min x,y 1 ⊤ x (5.21a) s.t. M x + Qy = r, (5.21b)
x, y ≥ 0.

(5.21c)

where V * P (r) denotes its optimal value function, and r = r 0 + Kξ ∈ R nr for some matrix K and some vector r 0 . The problems are equivalent in the sense that P * (ξ) = V * P (r 0 + Kξ). (5.22) Following Proposition 15 in the Appendix, the dual associated to Problem (5.21) is

V * D (r) ←-max µ µ ⊤ r (5.23a) s.t. M ⊤ µ ≤ 1, (5.23b
)

Q ⊤ µ ≤ 0. ( 5
(r + td) = V * D (r + td) (5.24)
which highlights the absence of duality gap. Now that the directional RHS sensitivity function is well-defined, thanks to Theorem 11, we can use standard results from the literature (as recalled by Proposition 16 in Appendix) to formulate the following theorem.

Theorem 8 (Direct application of [START_REF] Adler | A geometric view of parametric linear programming[END_REF][START_REF] Murty | Linear Programming[END_REF]). W is a closed interval (possibly unbounded). v is a continuous convex function, which is piecewise affine on a finite number of sub-intervals of W. Theorem 9. Consider any ξ and define d := Kξ. Assume that (i) V * P (r 0 ) = 0, (ii) ξ is s.t. V * P (r 0 + td) is defined for t > 0, on a (small) non-trivial interval. Then, the negotiation map t → P * (tξ) := V * P (r 0 + td) is non-decreasing on W ∩ R + . Proof. We proceed by combining the local affine description of t → V * P (r 0 + td) with the above mentioned results. Let us introduce the quantity D d (r) s.t.

D d (r) ←-max λ λ ⊤ d s.t. M ⊤ λ ≤ 1, Q ⊤ λ ≤ 0, λ ⊤ r = V * P (r).
8 Recalled as Theorem 11 in Appendix.

Thanks to the absence of duality gap over W, as shown in Equation (5.24), Theorem 12 used with assumption (ii) states that there exists a t ′ > 0 s.t. for any scalar t ∈ [0, t ′ ] we have V * P (r 0 + td) = V * P (r 0 ) + tD d (r 0 ) where D d (r 0 ) is finite. This formula is what we call Gauvin's formula in the Appendix.

Since V * P (r 0 ) = 0 by assumption (i), then 0 is a feasible vector for D d (r 0 ), and consequently D d (r 0 ) ≥ 0. Finally, since t → P * (tξ) is convex (from Theorem 8) and has a non-negative slope at t = 0, it is necessarily non-decreasing on W ∩ R + . Remark 32. For the general emergency problem -i.e. when there are R ≥ 2 negotiable sub-parameters -the negotiation maps t → P * j (tξ) behave differently. Showing that they are all continuous and piecewise affine can be achieved with little effort, but they are generally not convex, as the introductory example of Section 5.2.2 shows. This latter example also shows that even t → R j=1 P * j (tξ) is not necessarily convex. However, it has been conjectured that the negotiation maps t → P * j (tξ) remain non-decreasing functions of t. This conjecture is currently under investigation, at the time of writing this manuscript.

Remark 30. Assumption (i) in

Non-monotonicity of the optimal trajectories

Contrary to the negotiation maps, the directional optimal solution maps t → z * (tξ) defined by the outputs of HEGO are not necessarily component-wise monotonous. In the same fashion as in Section 5.2.2, consider the following example which illustrates this property. It is one of the reason why the heatmaps of the optimal time-of-flight variation, i.e. Figures B.9

Basic Example 4. Consider an optimization problem with

z = (z 1 , z 2 ) ⊤ ∈ R 2 ,

Emergency guidance method generalization

HEGO, as presented in Algorithm 2, is an emergency guidance method built upon a nominal guidance method that relies on Quadratic Programming. This convenient framework, with its linear constraints, allows the use of mature LP and QP solvers for the implementation of the problems LP j and Refine. However, it is possible to take some distance from this presentation, as other kinds of mathematical programming are often used to provide an approximation of PDG (ξ, p), as discussed in Chapter 1. Thus, let us describe the link between nominal and emergency guidance methods from a high-level perspective.

This section can be skipped with no loss of continuity.

Generalized notations

Let us denote by J(z, ξ) the cost that must be minimized, and by F eas (ξ, p) ⊂ R Nz the feasible set, which conveys the various control constraints, the dynamic model, etc. We insist that this set can convey much more general constraints than Equation (5.1). For example, it could stem from any direct collocation method, as presented in well-known references [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF][START_REF] Hargraves | Direct trajectory optimization using nonlinear programming and collocation[END_REF][START_REF] Stryk | Numerical Solution of Optimal Control Problems by Direct Collocation[END_REF], or other state-of-the-art methods [START_REF] Malyuta | Convex Optimization for Trajectory Generation[END_REF]. Also, it it desired that p ∈ Ω for some set Ω ⊂ R nneg .

Generalized emergency order

The goal is still to provide a trajectory, even when F eas (ξ, 0) is empty. To that purpose, instead of the 1-norm, let us introduce more general negotiation functions γ j : R n j → R + , for j = R, . . . , 1, assumed convex. It helps us define a generalized emergency order denoted ⪰ γ . The latter is also co-lexical, as ⪰ e , but in the sense of the negotiation functions γ j . Like in (5.3), a vector p a is said to be more negotiated than another vector p b s.t. p a ⪰ γ p b , if and only if

γ R p (R) a > γ R p (R) b or γ R p (R) a = γ R p (R) b and γ R-1 p (R-1) a > γ R-1 p (R-1) b , or . . . or γ R p (R) a = γ R p (R) b and . . . and γ 1 p (1) a > γ 1 p (1) b , or γ R p (R) a = γ R p (R) b and . . . and γ 1 p (1) a = γ 1 p (1) b .
Therefore, the general meaning of maximizing the launcher's integrity is to find the smallest p in the sense of ⪰ γ .

Generalized sequence of optimization problems

Let us here define a more general version of LP j and Refine. With the above-defined notations, the generalized nominal guidance problem is the following optimization problem

(Nominal) min z J(z, ξ), s.t. z ∈ F eas (ξ, 0) .
As before, what is of interest here is what happens when ξ is s.t. F eas (ξ, 0) is empty. Thus, we will seek to minimize a cost γ i for each parameter p (i) .

Let us introduce the first negotiation problem, dealing with the most critical negotiation parameter p (R) , s.t.

Γ * R ←-min z, p γ R (p (R) ), s.t. z ∈ F eas (ξ, p) , p ∈ Ω,
where Γ * R denotes the optimal value of this problem. Let S R be the set of points p s.t. there is a z ∈ F eas (ξ, p) and s.t. (z, p) minimizes γ R . In other words, S R is the set of minimizers of the latter problem, but projected on the negotiable parameters set. By successively negotiating the other parameters, in a similar fashion as LP j , we define the follow-up negotiation problems as

Γ * j ←-min z, p γ j (p (j) ), s.t. z ∈ F eas (ξ, p) , p ∈ Ω, p (j+1) ∈ S j+1 ,
where j = 1, . . . , R-1, and where S i denotes the set of minimizers of the i th problem, projected onto the negotiable parameters set (which is, by its recursive definition, a subset of the previous sets S i for i = j + 1, . . . , R). Consequently, the generalized version of the Refine problem becomes

z * , ⋆ ←-argmin z, p J(z, ξ), s.t. z ∈ F eas (ξ, p) , p ∈ Ω, p (1) ∈ S 1
where the notation z * , ⋆ means that the value of p is ignored, since it is not necessarily unique. The two latter problems can be re-written in a more convenient form.

Since the conditions z ∈ F eas (ξ, p) and p ∈ Ω do not change between these problems, it is possible to simplify the above-mentioned notations by switching the abstract condition p (j+1) ∈ S j+1 into

γ i (p (i) ) = Γ * i , ∀i = j + 1, . . . , R.
(5.25)

However, even for a convex function γ i , the equality constraint γ i (p (i) ) = Γ * i is numerically ill-posed, since it defines a non-convex level-set in general (e.g. when γ i (.) = ∥.∥ 2 ). Thankfully, it can be relaxed without loosing generality, by using only the inequality

γ i (p (i) ) ≤ Γ * i . ( 5
γ i (p (i) ) ≤ Γ * i , ∀i = 1, . . . , R (5.28d)

High-level description of safety margins

Enforcing the condition z ∈ F eas (ξ, p) for some value ξ may sometimes bring the variables of the problem to the frontier of what is feasible for a given p. Therefore, given a set M that contains 0, we would like to impose that if z is feasible for a value of p, then for every ∆p ∈ M there is another z ′ feasible for p + ∆p. Consequently, the condition z ∈ F eas (ξ, p) from the previous problems (5.27)

g ξ (z) ≤ K k=1 σ k g ξ (z k ) ≤ K k=1 σ k H(ξ).(p + ∆p k ) = H(ξ).(p + ∆p)
and likewise: A(ξ)z = b(ξ)+B(ξ).(p+∆p). Therefore, z belongs to F eas (ξ, p + ∆p), whence (5.29b). Constraint (5.29a) also holds, since it corresponds to the sub-case ∆p = 0, hence the conclusion.

Therefore, incorporating Proposition 14 into Problems 5.27 and 5.28 is a way to enforce safety margins while performing nominal an emergency guidance, using only a finite number of constraints.

Illustrations

In this section, four numerical examples are proposed along with qualitative discussions. Example 1 illustrates the basic principles of HEGO. Example 2 shows that a wide selection of negotiable parameters can be used together. Example 3 demonstrates the modeling capabilities offered by HEGO. Finally, Example 4 shows how emergency guidance scales to the 3D model. A quantitative analysis of this last example is proposed in Section 6.2 of Chapter 6.

Note that all of these examples have the same number of discretization points, i.e. N = 4, where N is defined in p. 56 right before Equation (4.1). For the 2D model (resp. the 3D model), it means that the size of z is N z = 15 (resp. N z = 22).

All the data presented in the examples below is normalized.

With the 2D model

Example 1 (Basic 2D scenario). Let us consider a simple choice of negotiable parameters, consisting in the incidence limit ∆α max and the final horizontal position ∆z f s.t.

p = ∆α max , ∆z f ⊤ .
In terms of hierarchy, we impose • p (1) = ∆α max (least critical),

• p (2) = ∆z f (most critical). This example is illustrated in Figures 5. [START_REF] Adler | A geometric view of parametric linear programming[END_REF] • p (1) = ∆α max (least critical),

• p (2) = ∆a max nor ,

• p (3) = ∆z f , • p (4) = ∆h f (most critical).
This example is illustrated in Figure 5.7.

Recall that ∆h f denotes the final altitude, which can seem surprising at first site. Why do not we use the final vertical speed ∆v f h instead? The reason is linked to the linearization used when we define QP (ξ). Indeed, ∆v f h appears, in practice, to be a

(τ 1 ) = µ 3 , α(τ 1 ) = µ 4 , qr (t f ) = µ 2(N +3)-1 and α(t f ) = µ 2(N +3
) . The Lispschitz-continuity stated in Theorem 7 can be observed on all the charts, except (g) whose Lipschitz-continuity is related to Corollary 3. Note that nominal guidance is performed up to ∆z 0 ≈ 0.3 (see Figure 5.5-(d)), highlighting the fact that significant active set changes can occur even with nominal guidance. Also, as one can observe in (h) between ∆z 0 = 0.0 and 0.2, ∆t * f has a non-zero though very small slope. This is mostly due to the fact that changing ∆t * f has a strong influence on simultaneously the vertical and the horizontal components of the trajectory, meanwhile ∆z 0 influences (almost only) the horizontal one. less useful lever than ∆h f . Mathematically speaking, it means that the image of the matrix 9 B ∆v f h does not describe the same vector space as the one of B ∆h f , and thus will not have the same ability to recover the changes in ξ. This remark also applies to the matrices B ∆v f h and B ∆h f . In practice, negotiating ∆h f is blindly allowed in the optimization problems of HEGO. However, when HEGO says that ∆h f < 0 is necessary, the optimal trajectory z * will reach the ground before reaching the new final altitude ∆h f . The state x at which the altitude of this trajectory reaches null altitude may be defined as the actual negotiated landing state.

Note that in Figure 5.7, dispersing the inputs w.r.t. ∆z 0 does not trigger the use of ∆h f , whose negotiation curve P * 4 remain flat in the sub-Figure (d). However, this example is the ground base of the numerical assessment provided in Chapter 6, where sufficiently rich scenarios are considered, and show how ∆h f is used.

Example 3 (2D scenario with repeated negotiable parameters). To demonstrate the modeling capabilities that Algorithm 2 offers, let us consider an example where several negotiable parameter are "repeated":

p = ∆α 1 max , ∆z f 1 , ∆α 2 max , ∆z f 2 ⊤
In terms of hierarchy, we impose

• p (1) = ∆α 1 max (least critical), • p (2) = ∆z f 1 , • p (3) = ∆α 2 max , • p (4) = ∆z f 2 (most critical).
This example is illustrated in Figure 5.8. This choice of negotiation parameters allows to negotiate the incidence and the final horizontal position alternatively. This may be helpful when the nature of the area neighboring the landing site becomes increasingly worse when moving away. To some extent, it can be applied to the terrain presented in Figure 1.3 from Chapter 1: the order of priority is to first negotiate the incidence, then the horizontal position (as long as it remains within the crops or the beach), then the incidence again, and finally trying to land in the forest or in the ocean. • p (1) = ∆α max (least critical),

With the 3D model

• p (2) = ∆a max nor ,

• p (3) = (∆z f , ∆y f ) ⊤ ,

• p (4) = ∆h f (most critical).

This example is illustrated in Figure 5.9. To illustrate behaviors that are not shared with the 2D rocket model, the input are dispersed w.r.t. ∆y 0 , the initial horizontal position component that is out-of-plane compared to the reference trajectory. These negotiation parameters are used for the quantitative analysis in Chapter 6.

Summary

In this chapter, we exposed a method to provide emergency guidance. It boils down to a sequential minimization of the amplitude of negotiable parameters, enforcing a prescribed hierarchy between these parameters, and is implemented using a finite number of LPs and a single QP. The method HEGO is capable of producing both nominal (when possible) and emergency guidance solutions thanks to a unified formulation.

Theoretical guarantees prove that the outputs of HEGO are Lipschitz-continuous w.r.t. its inputs, preventing the solutions from varying too fast for small changes in the inputs.

Numerical simulations have demonstrated how HEGO behaves on 2D and 3D examples, and that its outputs are consistent with the above-mentioned theoretical guarantees. Here, we propose quantitative assessments for HEGO. First, we provide highlevel comments regarding the implementation of HEGO. Then, we present a quantitative analysis of Example 4 from Chapter 5, by computing bivariate dispersions of the inputs on the 3D rocket model. Also, since HEGO aims at dealing with infeasible landing scenarios, its outcomes are compared with the vertical flight envelopes introduced in Chapter 3.

To improve readability, the figures of this chapter have been moved at its end.

General comments

The time required to run HEGO depends directly on the design choices presented in Chapters 4 and 5.

HEGO has been implemented in python and tested with cvxopt [START_REF] Andersen | CVXOPT: Convex Optimization[END_REF], mosek [START_REF] Andersen | The MOSEK interior point optimizer for linear programming: an implementation of the homogeneous algorithm[END_REF], glpk [START_REF] Fsf | GNU Linear Programming Kit[END_REF] and qpSWIFT [START_REF] Pandala | qpSWIFT: A real-time sparse quadratic program solver for robotic applications[END_REF] as underlying solvers for the LP and QP solvers. The benchmarks have been forced to run on a single CPU (in practice, tested on an Intel® Core™ i9-9900K, at 3.60GHz, and on an Intel® Core™ i7-8550U CPU at 1.80GHz). Its run time typically ranges between • a few milliseconds for the 2D rocket model with 2 negotiable parameters,

• and ≈ 60 ms for the 3D rocket model with 5 negotiable parameters.1

Major influence on run time

The parameters that naturally affect the computation time are the dimensions N z (i.e. main decision variable z) and n neg (i.e. negotiation parameter p). Let us denote by N opt := N z + n neg the dimension of the optimization variable involved in the sub-problems of HEGO. Recall that m is the dimension of the control variable, and that N is the number of discretization points (as shown in Figure 4.1). Since Cubic Splines are used for the description of the control corrections in Chapter 4, we get

N opt = (m + 3)N + 1 + n neg
showing the relative importance of the above-mentioned sizes. Also, the computation time may be influenced by the number of constraints, which is directly proportional to N c , as defined in Equation (4.8), and by the number of partitions of the negotiable parameter (i.e. R defined in Chapter 5).

Minor or no influence on run time

It is noteworthy that some parameters do not influence, or in a negligible way, the computation time. Among others, the dimension of ξ has a negligible impact on the run time of HEGO. Indeed, it only matters when computing the constraints right-hand side Gz ≤ h 0 + H ξ ξ + H p p, Az = b 0 + B ξ ξ + B p p.

A useful application would be to handle thinner descriptions of the atmosphere parameters, such as the wind for instance. Let us assume that one is able to measure the horizontal wind component at a high resolution for the low atmosphere layers (let us say 1 point of measure per 10 m, up to 10 km, for illustration purposes). Then, ξ would resemble: ξ = (∆x 0 ) ⊤ , w 0 m , w 10 m , w 20 m , , . . . , w 10 000 m

⊤

The dimension of ξ now equals 1009 in this case (when ∆x 0 is of of dimension 8 for the 3D rocket model of Chapter 2), though this has a non-significant impact on HEGO run-time.

Input dispersion on 3D rocket model

Let us consider the same reference trajectory as in Example 4, for the 3D rocket model. We pick the same choices of negotiable parameters as in Example 4, i.e. p = ∆α max , ∆a max nor , ∆z f , ∆y f , ∆h f ⊤ with R = 4 s.t.

• p (1) = ∆α max (least critical),

• p (2) = ∆a max nor ,

• p (3) = ∆z f , ∆y f ⊤ ,

• p (4) = ∆h f (most critical).

The purpose of this section is to quantify and assess the quality of the outputs of HEGO, over bivariate dispersions of the inputs. We consider the 15 inputs presented in Table 6.2. From the 105 pairs of inputs that can be formed from this list, we report results obtained with a selection of 20 pairs, enumerated in Table 6.3. These pairs have been selected for their representativeness, and because they demonstrate a wide variety of behaviors.

To categorize the outputs of HEGO, we use the notion of emergency mode. It is defined as the index of the most critical negotiable sub-parameter that has a nonzero 1-norm, thus taking its values within {0, . . . , R + 1}. Mathematically, using handy notations, it is defined s.t. where the values P * i (ξ) come from negotiation problems of HEGO. The color coding associated with the emergency modes is presented in Table 6.1.

Selection of figures

First, two pairs are detailed in Figures 6.1 and 6.2. They present the outputs of HEGO for two pairs, respectively (∆z 0 , ∆y 0 ) and (∆q 0 r , ∆q 0 c ). The first pair yields a typical collection of inputs that require the negotiation capabilities of HEGO, whereas the second pair shows that in some cases no negotiation is needed.

To ease the visualization, the pairs from Table 6.3 have been split in two batches, labeled A and B. For each batch, the following charts are provided:

• Emergency mode map (Figures 6. Several complementary comments are provided along with the Figures below.

Conclusion on the example

This example demonstrated that the negotiable parameters have very different impacts on how they help solve the emergency problem.

• ∆α max and (∆z f , ∆y f ) appear to be extremely relevant to recover feasible trajectories when inputs influencing the horizontal behavior of the landing change. ∆a max nor has a non-negligible though small influence on these changes.

• ∆h f has an influence on the vertical part of the landing, but the trajectories relaxed using ∆h f are not doable in practice, for the margin it offers is too thin.

These remarks are summarized in Table 6.4.

Comparison with vertical flight envelopes

As detailed in Chapter 5, Algorithm HEGO provide landing guidance for nominal and emergency problems. Thus, the set of inputs that do not require emergency guidance define the "nominal reachable set", in the sense of HEGO. Due to the approximations made in Chapter 4 (such as the discretization and the linearization), the set of inputs that are considered feasible without negotiation differ from the actual reachable set. The full characterization of the reachable set of the "complete" rocket model (i.e. the 2D or 3D rocket model) is not computationally tractable, due to the nonlinearities and the number of states involved. However, as detailed at the end of Chapter 3, it is possible to fully describe the reachable set of the rocket for purely Note that all of these trajectories have been computed using the 3D rocket model, even if it has been used with a purely vertical motion. These trajectories are represented in Figure 6.17.

vertical motion. Due to the high slenderness ratio of the rocket final-burn, this latter set is a coarse approximation of the actual reachable set of the complete model. Thus, we propose a comparison between the classification offered by HEGO in terms of emergency modes, and the vertical flight envelope defined in Chapter 3. We consider four different reference trajectories, which are listed in Table 6.5 and displayed in Figure 6.17.

For all four trajectories, a slice of the reachable set is considered, and HEGO is computed over a grid of inputs on each of these slices. In the chart (a) of their respective, the Max, Min and Min-Max surfaces defining the vertical flight envelope are defined at the end of Chapter 3.

The main observation is that HEGO is conservative. Indeed, on one hand, it declares some inputs "non-nominal" when they would be expected to be feasible without negotiating the constraints, as shown analytically in Chapter 3. On the other hand, their are no "false-positive", i.e. all the inputs declared "nominal" by HEGO are indeed inside the vertical flight envelope. Also, it is noteworthy that the reference trajectory significantly influences the outcomes of HEGO, though it does not change the latter conservative property. Contrary to what the chart (e) suggests, there is no discontinuity between the various scenarios. The fact that one can see a "gap" in the curves is only a matter of mesh refinement: a thinner mesh would show that this gap is continuously filled with curves. However, this gap shows that on some input subsets, the algorithm may have stiff output variations, as discussed in Figure 6.15.

Conclusion

Résumé

Ce chapitre propose un résumé succinct de ce manuscrit, quelques commentaires sur les sujets qui ont été volontairement omis, et des pistes de recherche pour d'éventuels travaux futurs.

The main focus of this thesis is on computing a landing trajectory for a reusable tossback vehicle in response to changes in the flight parameters. When this guidance problem is infeasible, one faces an emergency situation. In such cases it is acceptable to sacrifice some of the constraints originally formulated in the nominal (non emergency) situation.

A methodology (Algorithm HEGO) has been developed in the manuscript to solve this emergency guidance problem. It is applicable in the -relative vast-vicinity of a reference trajectory having a reasonably high slenderness ratio. It is capable of handling the sacrifice of constraints according to predefined extended colexicographic order. As illustration, results have been presented and have served to compute performance charts in Chapter 6. For actual missions, aerospace engineers would be interested in the quantitative version of Table 6.4.

A significant advantage of the proposed method, i.e. Algorithm HEGO, is that it is deterministic, and its implementation relies on mature technologies (LP and QP) for which off-the-shelf solvers are available. There are nearly no heuristic used to make the method work, apart from the tuning of the termination condition tolerances of the numerical solvers.

As announced in the introduction, several topics have been voluntarily considered out-of-scope in this thesis. Some of them represent possible future research directions, and are sketched here. The modeling choices of Chapter 4 could be tailored to other frameworks, to be able to handle singular arcs in the parametric control description for instance.

As far as the overall G&C system is concerned, the control part has been considered out-of-scope for this manuscript. The interplay between the control and the guidance algorithms has been swiftly discussed in [START_REF] Ménou | Nominal And Emergency Rocket Landing Guidance Using Quadratic Programming[END_REF]. Simulations gathering both parts of the G&C system would provide a more complete performance assessment.

Negotiating vertical components (∆v f h or ∆h f ) has proved to be ill-posed, as shown in Chapter 6. To obtain robustness regarding the vertical motion, it is needed to go back to the mission design itself, and eventually consider other landing strategies than the classic tossback trajectory structure, pictured in Figure 1.2. More generally, the choice of the reference trajectory has a strong impact on the final performances.

The sensitivity-based PDG method presented in Chapter 4 has a sufficient accuracy for our application (i.e. final burn guidance of a tossback vehicle with high slenderness ratio) but depends significantly on the chosen reference trajectory and its generalization to more complex maneuvers is limited. Its application to large diverts of high agility vehicles would require a dedicated study with extensive numerical benchmarks. However, it may be of interest to apply this sensitivity-based guidance method to multi-phase problems. Indeed, combining the re-entry glide and the final burn phases in a single guidance problem would be a relevant problem, for which the sensitivity based approach would have a great potential.

Our implementation of HEGO conveniently builds upon the linear description of the constraints. However, following the generalization discussed in Section 5.6, it would be interesting to develop a version of HEGO with other underlying guidance methods, such as successive convexification or pseudo-spectral methods, which could be applied to large-divert problems.

Finally, the emergency problem that we formulated is only one way to describe the problem of landing guidance relaxation. Our modeling is not sufficient to handle disjunctive scenarios, where one would need to choose between two separate landing sites for instance. This represents an important research direction, that could be worth exploring.

The proof of this property is standard material that can be found in most optimization text-books (see e.g. [START_REF] Boyd | Convex Optimization[END_REF]Ch. 5]). It is provided here for completeness and for its tutorial aspect. Also, beware of the fact that the role of µ and λ is inverted in the previous proposition and its proof compared to Equation (A.2). Note that if the condition c -A ⊤ µλ = 0 is not feasible -or equivalently that A ⊤ µ ≤ c is not feasible -then g equals -∞ for any value of µ and λ, and therefore the optimal value of the dual is -∞.

A direct consequence of Proposition 15 and the Strong Duality Theorem presented above is that

v(b) = min{c ⊤ x | Ax = b, x ≥ 0} = max{µ ⊤ b | A ⊤ µ ≤ c}.
whenever the primal problem is feasible and finite. By convention, when the primal problem is infeasible, v takes the optimal value of the dual problem (possibly ±∞).

- 1

 1 Saturates x between a and b. Sgn (.) = Sign function For a ∈ R, Sgn (a) := if a < 0. IVP = Initial Value Problem OCP = Optimal Control Problem STM = State Transition Matrix See Appendix A.2.2. LP = Linear Programming QP = Quadratic Programming NLP = Non Linear Programming RHS = Right-hand side OoP = Out of Plane Mathematical nomenclature specific to Chapter 2 e ⋆ = Unit vector In direction of ⋆. X, X = Vector, Norm For X a vector of R 3 , then X = ∥X∥ 2 . R (a, γ) = Rotation 3 × 3 rotation matrix of axis a and angle γ. Nomenclature for the 2D model α = Incidence (In 2D) α is signed: -90 • < α < 90 • . (e z , e h ) = Earth frame xv xvi Nomenclature (C A , C N ) = Rocket body frame θ = Attitude (In 2D) θ is signed: -90 • < θ < 90 • . h = Altitude v h = Vertical speed z = Horizontal position v z = Horizontal speed q r = Engine flow m = Total mass a nor = Normal acceleration a nor is signed.

  ξ = Pitch (e z , e y , e h ) = Earth frame (C N , C A , C P ) = Rocket body frame θ = Attitude (In 3D) θ is unsigned: 0 • ≤ θ < 90 • . ζ z , ζ y = Projected attitudes Defined in Figures 2.4 and 2.5. α z , α y = Projected incidences Defined in Figure 2.5. z = Horizontal position (With respect to e z .) y = Horizontal position (With respect to e y .) h = Altitude v z = Horizontal speed (With respect to e z .) v y = Horizontal speed (With respect to e y .) v h = Vertical speed m = Total mass q r = Engine flow (Real) q c = Engine flow (Controlled)

  Reference state State of size n. ū = Reference control Control of size m. η = Reference parameter Parameter of size n η . tf = Reference time-of-flight µ = Discretized control Size N µ . For Cubic Splines, N µ = m(N + 3).

τ

  i = Control time instance N elements s.t. 0 = τ 0 < τ 1 < . . . < τ N = 1. τ ′ i = Constraint time instance N c elements s.t. 0 = τ ′ 0 < τ ′ 1 < . . . < τ ′ Nc = 1. ∆t f = Final time change z = Decision variable Size N z = N µ + 1, s.t. z = (µ ⊤ , ∆t f ) ⊤ .

Figure 1 . 1 :

 11 Figure 1.1: (Left) DC-X take-off and landing in 1993. (© New Mexico Museum of Space History) (Right) Perseverance landing site, 2021. (© ESA/DLR/FU-Berlin/NASA/JPL-Caltech)

Figure 1 . 2 :

 12 Figure 1.2: Typical flight phases of a tossback vehicle, from take-off to landing (trajectory not to scale). The higher the slenderness ratio, the more vertical the flight trajectories. (© Google Earth V 9.168.0.0, (July 28, 2022), France, Landsat/Copernicus)

Forest

Figure 1 . 3 :

 13 Figure 1.3: Trying to maximize launcher's integrity by relaxing well-chosen constraints.

Figure 1 . 4 :

 14 Figure 1.4: High-level summary. Nominal and emergency guidance methods presented in this thesis, with references to the associated chapters.

Figure 2 . 1 :

 21 Figure 2.1: Air jet influence.Qualitative change of the wind flow with (left) and without (right) an air jet pushing in front of a moving rocket. See[START_REF] Nonaka | Vertical Landing Aerodynamics of Reusable Rocket Vehicle[END_REF] and[START_REF] Marwege | First Wind Tunnel Data of CALLISTO Reusable VTVL Launcher First Stage Demonstrator[END_REF] for experimental results.

Figure 2 . 2 :

 22 Figure 2.2: Planar rocket model. Axis, angles and forces. e Vr (respectively e Or ) denotes the unit vector parallel (resp. orthogonal) to V r .

Figure 2 . 4 :

 24 Figure 2.4: Angles ζ z and ζ y in the frame (e z , e y , e h ), for the proof of Equation (2.8).

(Figure 2 . 6 :

 26 Figure 2.6: Relations between the lift, the drag and the relative speed vectors. Here, the angle ν denotes the oriented angle between C N and L, defined only for α ̸ = 0.

  a u.s. nor = ∥F -(F • e Vr )e Vr ∥ m .

Fig- ure 2 . 6 ,

 26 T, D, L, e Vr and e L all belong to the plane (e Vr , e L ), or equivalently to the plane (e Or , e Vr ), where e Or and e Vr are orthogonal. Consequently, the term F -(F • e Vr )e Vr only has a component along e Or , which yields ∥F -(F • e Vr )e Vr ∥ = |F • e Or | Using the latter expression with the definitions of T, D, L, e Vr , e L and C A yields a nor := L cos α -(T + D) sin α m .

y 2 : 2 max 2 . ( 3 . 8 )

 22238 = -ā cc T max and ȳ1 := ācc T Let us define F ⊂ D the flight envelope, as the set of states y lying between Σ max and Σ min (in terms of altitude), and satisfying y 1 ≤ ȳ1 , y 2 ≥ y 2 , and m -≤ y 3 ≤ m + .(3.9)

Figure 3

 3 Figure 3.1:

Figure 3 . 1 :

 31 Figure 3.1: Possible scenarios for (λ 2 , λ2 ). The origin is prohibited due to Proposition 4.

Figure 3 . 3 :

 33 Figure 3.3: Three Max-Min-Max trajectories, with varying first time of switch t ′ 1 . Maximum final mass is obtained for t ′ 1 = 0, i.e. for a Min-Max thrust program.

  A useful by-product of the proofs from above is the full characterization of the reachable set for the vertical motion (3.37), as shown in Figure 3.2. Using the notations from Figure 3.2, for any map v(y) : R 3 → [-1, 1], every feedback control law u : R 3 → R s.t.

  Theorem 1 relies on Assumption 4 to rule out the Max-Min-Max programs, and keep the Min-Max programs only. A greater number of numerical simulations, aiming at computing dy f 3

Figure 4 . 1 :

 41 Figure 4.1: Control discretization and time instances of the constraints for N = 3.The correction is described by a parametric function τ → u µ (τ ). Here is represented a scalar Cubic Spline, described by its values µ 0 , . . . , µ 3 at several time instances, and by its slopes µ 4 and µ 5 at the starting and end-points. The inequality constraints are enforced on the subdivision τ ′ 0 , . . . , τ ′ Nc .

Figure 4 . 2 :

 42 Figure 4.2: Representation of the optimal points of Basic Example 1 for -1 ≤ ξ ≤ 1, in the plane (x 1 , x 2 ). Point A is sometimes called a cusp in the literature [23].

  and [g j ] z (z 0 , 0) (all j) are linearly independent, (ii) the conditions of Lemma 5 are satisfied at z 0 for the multipliers ν 0 and λ 0 , (iii) SCS holds for (z 0 , ν 0 ), then (a) z 0 is a local isolated minimizing point of problem NLP (0) and the associated Lagrange multipliers ν 0 and λ 0 are unique, (b) for ξ in a neighborhood of 0, there exists a unique, once continuously differentiable vector function y(ξ) = [z * (ξ), ν * (ξ), λ * (ξ)] ⊤ satisfying the second-order sufficient conditions for a local minimum of NLP (ξ) s.t. y(0) = (z 0 , ν 0 , λ 0 ), and hence z * (ξ) is a locally unique local minimum of NLP (ξ) with associated unique Lagrange multipliers ν * (ξ) and λ * (ξ),

Figure 4 . 3 :

 43 Figure 4.3: Comparison between the exact solution of Example 1 and the solution returned by QP (ξ), for -1 ≤ ξ ≤ 1.

Figure 4 . 3 . 18 .

 4318 Figure 4.3. Remark 18. Similar methods aiming at computing an expansion of z * (ξ) can be found in more recent work. See for example the work of Bonnans & Shapiro [19, Sec. 5.2]. Remark 19. From a high-level point of view, this Chapter describes how to con-

. 25 )

 25 Composing the previous expressions yields the directional first-order expansion of the statex[τ, z nlp (εξ, εp), εξ] = xlin (τ, εξ, εp) + o(ε).Since xlin is an approximation of the augmented state, it can be split in half: the first n components form the state approximation x lin , and the last component form the time-of-flight approximation t lin f xlin (τ, ξ, p) =   

Figure 4 . 4 :

 44 Figure 4.4: Summary of the nominal guidance method, as presented in Chapter 4.

Figures 4 .

 4 5 illustrate a dispersion of the input variable ξ along a single component, namely ∆v 0 h , for both positive and negative values. The inputs are voluntarily dispersed over a small range of values.

  -(c) and (d), non-local constraints are activated for sufficiently large values of the input, as demonstrated by the trajectories that reach the incidence bounds.

Figure 4 . 5 :

 45 Figure 4.5: First-order correction, for ∆v 0 h varying in the vicinity of 0. The sub-Figure (a) is the main purpose of this example, showing that the green curve is a second order residual w.r.t. to the input ∆v 0 h . Note that sub-Figure (c) shows very late corrections in the engine flow. The weighting matrix P has been chosen to favorearly corrections. Since QP (ξ, p) returned late corrections anyway, it means that earlier flow corrections would have required even higher incidence corrections (which is partially explained by the high dynamic pressure at the beginning of the descent).

Figure 4 . 6 :

 46 Figure 4.6: Changes in the input variable leading to local changes in the active set.The input variable used in this example is ∆z 0 .

Figure 4 . 7 :

 47 Figure 4.7: Non-local behavior of QP (ξ, p), when ∆z 0 and ∆y 0 change. For large values of ∆z 0 and ∆y 0 , several constraints start to be triggered, such as the incidence bounds (upper and/or lower bounds).

  s.t. n neg = 2 and R = 2, making both sub-parameters are thus scalars. Then (Least critical) (1, 1) ⪯ e (3, 2) ⪯ e (4, 2) ⪯ e (1, 3) ⪯ e (-2, 3) (Most critical). Basic Example 3. Let us consider negotiable parameters s.t. n neg = 3 and R = 2 s.t. p = (p

Figure 5 . 1 :

 51 Figure 5.1: Pictorial representation of the possible values for the input ξ.

Figure 5 . 2 :

 52 Figure 5.2: Illustration of the constraint for the problem of Section 5.2.2 w.r.t. the input values, using HEGO.

Figure 5 . 3 :

 53 Figure 5.3: Curves associated to the example of Section 5.2.2. (Left)The first remark is that all the quantities displayed are indeed continuous w.r.t. the input ξ. On the bottom chart, there are four distinguishable areas. (A) corresponds to the nominal scenario, when no negotiation is needed. (B) and (C) corresponds to scenarios that require respectively 1 and 2 non-zero negotiation parameter values to recover feasibility. Finally, scenario (D) is when there are no options left., and no allowed values of p 1 or p 2 can help recover feasibility for z. (Right) Without enforcing any hierarchy, the previous zones (B), (C) merge into a single zone (B ′ ), where both variables p 1 and p 2 are negotiated at the same time. The penalties differ from the previous case, and so do the optimal values z 1 and z 2 .

(5. 10 )

 10 It is also possible to re-write Problem (5.10) into the following canonical QP (Canonical QP)

Lemma 6 .Lemma 7 .

 67 If {x | Ā0 x = r0 , x ≥ 0} is not empty, then Problem (5.10) has no duality gap, i.e. Problem (5.10) and Problem (5.12) have optimal values and they are equal. Under the assumption of Lemma 6, let (x * , µ * , λ * ) and (x, μ, λ) be optimal solutions of both Problems 5.10 and 5.12. Then, Dx * = Dx.

4Proposition 12 .

 12 0 which gives D(x * -x) = 0 since D ⪰ 0, hence the desired result. See e.g. [21, Section 5.2]. 5 Theorem 11, recalled in the Appendix.Let us now introduce the setχ(ξ) := x | ĀR x = rR , x ≥ 0 . If χ(ξ) is non-empty, then the value of z * = z +z -computedby Algorithm 2, obtained by successively solving Problems 5.9 and 5.10, exists and is unique.

Corollary 1 .

 1 Algorithm 2 returns a solution for the input ξ if and only if ξ ∈ Λ.Proof. If χ(ξ) ̸ = ∅, Proposition 12 shows that Algorithm 2 has a solution. Otherwise, χ(ξ) = ∅, implies that LP R is infeasible, and then Algorithm 2 fails. The set Λ is not empty, since 0 belongs to Λ by definition of QP (ξ, p). Proposition 13. Λ is convex.

ATheorem 4 (

 4 feasible point for Problem (5.13) is a point x that satisfies (5.13b) and (5.13c). A solution point for Problem (5.13) is a feasible point that is minimal for (5.13a). Adapted from [59, Thm. 2.4]). Let the Linear Program (5.13) have non-empty solution sets S 1 and S 2 for right-hand sides (b 1 , d 1 ) and (b 2 , d 2 ), respectively. There exists a constant L > 0 s.t. for each x 1 * ∈ S 1 , there exists an x 2 * ∈ S 2 s.t.

Theorem 5 (

 5 Adapted from [59, Thm. 2.2], Lipschitz-continuity of feasible points of linear inequalities and equalities). Let the conditions Ax = b and Cx ≤ d have non-empty feasible sets F 1 and F 2 for the right-hand sides (b 1 , d 1 ) and (b 2 , d 2 ),

  and B.20, and the curve of ∆t * f , in Figure B.23, exhibit such complex patterns.

Figure 5 . 4 :

 54 Figure 5.4: Illutration of Basic Example 4. This example shows that there can be components z * i among the outputs of HEGO that are non-monotonous.

  and 5.6, where the input ξ is dispersed for positive values of the change in initial horizontal position ∆z 0 . Example 2 (Advanced 2D scenario). Let us consider a more advanced version of Example 1, where the negotiable parameter is now p = ∆α max , ∆a max nor , ∆z f , ∆h f ⊤ In terms of hierarchy, we impose

2 Figure 5 . 5 :

 255 Figure 5.5: Dispersion of ∆z 0 over [0, 1], for the 2D rocket model from Example 1. The curves of sub-Figures (a), (b) and (c) are plotted for 30 values of ∆z 0 . The similarity between the introductory example of Section 5.2.2 and this actual example is clear in the sub-Figure (d).

Figure 5 . 6 :

 56 Figure 5.6: Dispersion of ∆z 0 over [0, 1], for the 2D rocket model from Example 1.All these charts have very different y-scales. The charts (c) and (d) are zoomed views of (a) and (b). Except for the charts (c) and (d), where the blue dots represent the computed values, all the curves have been plotted for 301 values of ∆z 0 . The charts (a) to (f) shows a subset of µ, of size 2(N + 3). Precisely, using the nomenclature from Equation (4.1), q r (τ 1 ) = µ 3 , α(τ 1 ) = µ 4 , qr (t f ) = µ 2(N +3)-1 and α(t f ) =

4 Figure 5 . 7 :

 457 Figure 5.7: Dispersion of ∆z 0 over [-2, 1], for the 2D rocket model from Example 2.In sub-Figure(c), a thinner mesh (not required for the purpose of this example) would reveal the presence of blue-dot trajectories between the orange and the red ones, on the left of the reference trajectory. The asymmetry in the negotiation maps (sub-Figure(d)) comes from the reference trajectory itself, which is a non-trivial curve in the plane (z, h).

H 2 | 3 | h f | = * 4 Figure 5 . 9 :

 23459 Figure 5.9: Dispersion of ∆y 0 over [-1, 1], for the 3D rocket model. Dispersing the inputs along ∆y 0 leads to the computation of out-of-plane trajectories. The dispersion of the inputs of HEGO will be analyzed quantitatively on this exact same scenario in Chapter 6.

EmergencyMode 1

 1 (ξ) = . . . = P * R (ξ) = 0, R + 1 if LP R (ξ) is not feasible, arg max i=1,...,R { i | P * i (ξ) ̸ = 0} otherwise.(6.1)

3

 3

  and 6.4),• Projected incidence α y (Figures 6.5 and 6.6),

Figure 6 . 1 :

 61 Figure 6.1: Dispersion of (∆z 0 , ∆y 0 ). Dispersing this pair leads to in-plane and out-of-plane trajectories. The negotiation maps (charts (f) to (i)) have a structure deeply linked to the partition described in the emergency mode map from chart (b).

Figure 6 . 2 :

 62 Figure 6.2: Dispersion of (∆q 0 c , ∆q 0 r ). The inputs correspond to purely nominal scenarios. The negotiation maps are not displayed, since they are all constant and equal to zero.

Figure 6 . 3 :

 63 Figure 6.3: Emergency modes. Batch A. Each chart represents the emergency modes obtained by dispersing a given pair of inputs. For instance, in (a), ∆v 0 z w.r.t ∆z 0 means that ∆z 0 is in the x-axis, and ∆v 0z is in the y-axis. Since the role of the outof-plane variables (such as ∆y 0 , ∆v 0 y or w y,0 ) is symmetric in this problem having a planar reference trajectory, their associated charts have an axis of symmetry.

Figure 6 . 4 :

 64 Figure 6.4: Emergency modes. Batch B. Finding purely horizontal or purely vertical lines separating different emergency modes for some specific pairs reveals that the associated input variables are uncorrelated. Among others, the variables conveying the vertical motion (h, v h , m, Isp) are mostly uncorrelated with the ones conveying the horizontal motion (z, v z , y, v y ), especially when considering the emergency mode of ∆h f . The pairs (i) and (m) are typical supporting examples. Moreover, contraryto what appears to be, charts (p) and (q) are not feature-less. The specific case of the pair (q) is detailed above in Figure6.2.

Figure 6 . 5 :

 65 Figure 6.5: Projected incidence α y (in-plane) w.r.t. normalized time. Batch A. In most charts, the curves describe a sort of pivot around a the normalized time τ = 0.5. This behavior is directly linked to the Cubic Spline description of the correction, detailed in Chapter 4. Note that for some pairs -e.g. pair (e) -the incidence is negotiated but not always up to its maximum value, even when more critical parameters are negotiated first. The main reason is due to the way the active set changes w.r.t. to the input, which is illustrated in detail in Figure B.23.

Figure 6 . 6 :

 66 Figure 6.6: Projected incidence α y (in-plane) w.r.t. normalized time. Batch B. It is normal that the pair (l), with purely out-of-plane variables, has nearly zero impact on α y , since the latter conveys the in-plane projected incidence.

Figure 6 . 7 :

 67 Figure 6.7: Projected incidence α z (out-of-plane) w.r.t. normalized time. Batch A.Contrary to what the chart (e) suggests, there is no discontinuity between the various scenarios. The fact that one can see a "gap" in the curves is only a matter of mesh refinement: a thinner mesh would show that this gap is continuously filled with curves. However, this gap shows that on some input subsets, the algorithm may have stiff output variations, as discussed in Figure6.15.

Figure 6 . 8 :

 68 Figure 6.8: Projected incidence α z (out-of-plane) w.r.t. normalized time. Batch B.

Figure 6 . 9 :

 69 Figure 6.9: Controlled engine flow q c w.r.t. normalized time. Batch A. The highly different influence of the horizontal variables and the vertical ones on the engine flow is clearly pictured by these charts. A change on the input pair (c) only implies little changes in the controlled engine flow. However, the pair (e) has a high impact on the engine flow, especially due to the presence of ∆v 0 h .

Figure 6 . 10 :

 610 Figure 6.10: Controlled engine flow q c w.r.t. normalized time. Batch B.

Figure 6 . 11 :

 611 Figure 6.11: Heatmap of the first negotiation, i.e. P * 1 = |∆α max |. Batch A.

Figure 6 . 12 :

 612 Figure 6.12: Heatmap of the first negotiation, i.e. P * 1 = |∆α max |. Batch B.

Figure 6 . 13 :

 613 Figure 6.13: Heatmap of the third negotiation, i.e. P * 3 = ∥(∆z f , ∆y f ) ⊤ ∥ 1 . Batch A. These negotiation maps have a structure deeply linked to the partition described in the emergency mode maps, as shown in Figure 6.3.

Figure 6 . 14 :

 614 Figure 6.14: Heatmap of the third negotiation, i.e. P * 3 = ∥(∆z f , ∆y f ) ⊤ ∥ 1 . Batch B.

Figure 6 . 15 :

 615 Figure 6.15: Several zooms on stiff parts of the pair (∆z 0 , ∆h 0 ), illustrating that the negotiation maps are indeed Lipschitz-continuous in practice. All the x-axis convey ∆z 0 , and all the y-axis convey ∆h 0 .

Figure 6 . 16 :

 616 Normal load a nor on two different input pairs. H o ri z o n ta l p o s it io n z -1 -0.5 0 0.5 1 H o r iz o n t a l p o s it io n y time: = t/tf mdry (f) Total mass m

Figure 6 . 17 :

 617 Figure 6.17: Reference trajectories selected for the benchmark against the vertical flight envelope, as listed in Table6.5.

Figure 6 . 18 :

 618 Figure 6.18: Vertical envelope and associated slices (Reference A). As shown in the chart (b), the red slice shows what happens when dispersing the total mass m w.r.t. the vertical speed v h . The light yellow area conveys the area which is out of the flight envelope. On the contrary, the blue area is the area inside it. The dashed line conveys the separation between these areas.

Figure 6 . 19 :

 619 Figure 6.19: Vertical envelope and associated slices (Reference B).

Figure 6 . 20 :

 620 Figure 6.20: Vertical envelope and associated slices (Reference C).

Figure 6 . 21 :

 621 Figure 6.21: Vertical envelope and associated slices (Reference D).

Proof.

  Introduce the LagrangianL(x, µ, λ) := c ⊤ x + µ ⊤ (b -Ax)λ ⊤ x.Let us form the function g s.t.g(µ, λ) := inf x L(x, µ, λ) Since g(µ, λ) = µ ⊤ b + inf x (c -A ⊤ µλ) ⊤ x, then g(µ, λ) =      µ ⊤ b if c -A ⊤ µλ = 0 -∞ otherwise.Therefore, the dual of (A.5) is sup µ,λ g(µ, λ),s.t. λ ≥ 0.Using the Strong Duality Theorem (Theorem 11 above) with the fact that v(b) is finite and x → c ⊤ x is convex, makes the sup of the latter problem a max. Finally, after simplification of λ = c -A ⊤ µ ≥ 0, ones gets (A.6).

Theorem 12 (. 7 ) 33 .•

 12733 From [41, Thm. 1]). Under the assumption of Proposition 15, for any direction d and for any scalar t > 0 sufficiently small, we always havev(b + td) = v(b) + t sup{λ ⊤ d | A ⊤ λ ≤ c, λ ⊤ b = v(b)}.(ARemark Theorem 12 must be understood as follows: If it is known that t → v(b + td) exists on some interval [0, t ′ ] for t ′ > 0, then the "sup" term in Equation (A.7) is a "max", and Equation (A.7) holds (at least) on a non-trivial sub-interval of [0, t ′ ].

Figure B. 10 :

 10 Figure B.10: Heatmap of second negotiation P * 2 = |∆a max nor |. Batch A.

  .2)

	∧	ROCKET BODY	∧	∧	ROCKET BODY	∧
			Boundary layer stays attached to the body
	∧	ENGINE	∧	∧	ENGINE	∧
			Boundary layer detaches from the body
		∧

  , e y , e h ) is rotated by Ψ around e z , giving (e z , C P , e h ′ ). Then, the later is rotated by ξ around C P to give (C N , C P , e A ). Finally C A = -e A , leading to the orthonormal basis (C N , C A , C P ).

		ξ	e h ′	Ψ	e h
	e A = -C A			Yaw: Ψ Pitch: ξ Here, Ψ > 0, ξ > 0.
					ξ
	e z	Ψ			•	e y Ψ C P
		ξ			V r
		C N			C A
	Figure 2.3: Rocket orientation, based on yaw (Ψ) and pitch (ξ) angles. First, (e z

  A = -C A convey enough information to establish the change of variables between (Ψ, ξ) and (ζ z , ζ y ). By taking the ratio of the first two components of e A , the following relations are obtained

	Equation (2.7) and (2.8) and e Ψ = ζ z and	tan ξ = cos Ψ tan ζ y .	(2.9)
	Then, using the shortcuts t z := tan ζ z , t y := tan ζ y and t z,y := 1 + t 2 z + t 2 y , we get
	hence (2.8).	.4, one has e A = (h y , -h z , v) ⊤ ,
	Remark 4. Equation (2.7) is well defined whatever the values of Ψ and ξ. However, Equation (2.8) has a singular definition when ζ z = 90 • or ζ y = 90 • . This raises
	multiple comments:
	1. In all the scenarios studied in the thesis, only close-to-vertical trajectories are
	considered, where the rocket remains far from these singularities.
	2. The singularity in Equation (2.8) is only a problem of definition. Indeed, e A
	can be extended by continuity everywhere, except when ζ z and ζ y equal ±90 • simultaneously. For example, for a given |ζ y | < 90 • , then e A (ζ z ) -→ ζz↑90 • e z .
	3. One could legitimately wonder whether a numerical method could fall into one
	of these singularities during intermediate computations. As far as this thesis
	is concerned, the guidance methods exposed below only require the evaluation
	of the dynamic function f and its derivatives along a prescribed reference tra-
	jectory before the flight (more details on this topic in Section 4.3). Since the
	evaluation of f is not required on-board, this singular definition at ±90 • does not present any risk for our applications.

  .12) Representation of the projected angles: ζ z , ζ y , α z , α y . Note that C A , e Vr and V r are not coplanar to (e z , e h ) nor (e y , e h ).wheret z = tan ζ z , t y = tan ζ y , T z = tan(α z + ζ z ) andT y = tan(α y + ζ y ), to alleviate the writing. This expression stems from Equations (2.10) and (2.11) and ∥C A × e Vr ∥ = ∥C A ∥ ∥e Vr ∥ sin C A , e Vr = sin α

	e A		e h			e A	e h
		ζ y				ζ z
					Here, ζ z > 0 and ζ y > 0.
	e z	e y	•			e z	•		e y
			ζ y	α y	e Vr	|v h |	ζ z	α z	e Vr	|v h |
			α y + ζ y	C A		α z + ζ z	C A
					V r				V r
			-(v z -w z (h))			v y -w y (h)
	Figure 2.5: where

  Assumption 3 depends on b and a, which depend on the bounds on y 2 and y 1 . Though the estimates of y 2 and ȳ1 provided in Equation (3.8) are coarse,

	Assumption 3. The constants in (3.18) are s.t. b <	√ a.
	Remark 11.

they are sufficient for the numerical application discussed below. If needed, analytic bounds sharper than (3.8) could be computed. Proposition 7. For scenario 2a, γ(t) < 0, ∀t ∈ [0, t f ].

  Proof. γ can be zero at most on an isolated point. Indeed, γ is continuous and if there is t 0 s.t. γ(t 0 ) = 0, then, from Lemma 4, it cannot be zero for greater times. Therefore, from (3.19), Γ can be zero at most on two isolated points.Proposition 8 shows that the two remaining scenarios (1c and 3a) correspond to Max-Min-Max structures. It enables us to state the main result below. Under Assumptions 1, 2 and 3, and for y 0 in F, any solution of Problem (2) is necessarily a Max-Min-Max thrust program, where one or two arcs may be absent.

	(3.25)
	The conclusion stems from comparison Lemma 13.
	Proposition 8. Under the assumptions of Lemma 4, the sign of Γ changes at most twice on [0, t f ].
	Proposition 9.

[START_REF] Bryson | Applied optimal control: optimization, estimation and control[END_REF]

Thus, γ a (t) < 0 for t > t γ . Moreover, z λ (t γ ) = z a (t γ ) holds and since λ 2 < 0, for any t ∈ [t γ , t f ] żλ (t) = F (t, z λ (t)) and ża (t) ≤ F (t, z a (t)).

  .27) It describes the above-mentioned implicit dependence of (t 2 , t f ) on t 1 . When applicable, the IFT used on (3.27) provides us with the differentiability and the value of the derivatives of t 2 and t f w.r.t. t 1 , as

Definition 4 .

 4 The finite-dimensional approximation of the problem PDG (ξ, p) is defined by NLP (ξ, p), which is the following non-linear optimization problem

	min z	J(z, ξ)	(4.10a)
	s.t. h(z, ξ) ≤ H p p, g(z, ξ) = B p p.	(4.10b) (4.10c)

Remark 15 (State modeling choices). Describing the state using Equation (4.7) is

Table 5 .

 5

	Parameter name	Variables	Negotiable? Comment
	Final horizontal positions Final horizontal speeds Final altitude Final vertical speed Incidence bound (2D) Projected incidence bound (3D)	∆z f , ∆y f ∆v f z , ∆v f y ∆h f ∆v f h ∆α max	✓ × × ≈ ✓	If landing area is solid and flat Otherwise the rocket would tilt at landing. See Example 2 for a discussion. Tiny margin, imposed by landing gear design Safety and flight quality bound.
	Engine flow bounds Normal acceleration	∆q min , ∆q max ∆a max nor	× ✓	Physical constraint. Safety and flight quality bound.

1: List of negotiable parameters in PDG (ξ, p). under the form p low ≤ p ≤ p up . (5.2)

  Corollary 2 (Lipschitz continuity of the optimal value function of LPs w.r.t. RHS perturbations). Let the Linear Program (5.13) have non-empty solution sets S 1 and S 2 for right-hand sides (b 1 , d 1 ) and (b 2 , d 2 ), respectively, with associated optimal values p * 1 and p * 2 . Then, there exists a constant K > 0 s.t.

	|p * 1

  Theorem 7 (Lipschitz-continuity of Algorithm 2). There exists a constant L > 0 s.t. for any inputs ξ 1 and ξ 2 in Λ, the unique solutions (z * ) 1 and (z * ) 2 returned by Algorithm 2 satisfy

  If V * P (r + td) is feasible and finite, then V * D (r + td) is finite by the Strong Duality Theorem 8 . For the converse case, note that V * D (r + td) is always feasible since µ = 0 satisfies (5.23b) and 5.23c. Also recall that the Dual of V * D (r + td) is V * P (r + td). Therefore, using once more the Strong Duality Theorem, if V * D (r + td) is finite then V * P (r + td) is feasible and finite. This gives the conclusion. For t ∈ W, let us define the function v s.t.

	.23c)
	Also, we will consider a RHS change in an arbitrary direction d. Let us define
	W := { t ∈ R : V * P (r + td) is feasible and finite} = {t ∈ R : V * D (r + td) is finite}. Lemma 11. For any r, V * P (r + td) is feasible and finite if and only if V * D (r + td) is
	finite.
	Proof. v(t) := V * P

  Theorem 9 is guaranteed by the assumption stating that the reference trajectory must satisfy the constraints of NLP (ξ, p): for null inputs -i.e. RHS equals r 0 in Problem (5.21) -there is no need for negotiation (i.e.

	V * P (r 0 ) = 0).
	Remark 31. Assumption (ii) in Theorem 9 is not over restrictive. It means that we
	consider only inputs ξ s.t. the negotiation problem remains feasible when exploring
	the inputs in this direction.

  To alleviate the writing, we denote by a function "c" the latter constraints s.t. these are satisfied if and only if c(ξ, z, p) ≤ 0. In this case, the (single) negotiation problem associated to these constraints is simply

	min z 1 ,z 2 ,p s.t. c(ξ, z, p) ≤ 0, p p ≥ 0. Moreover, given the arbitrary cost J(z) := 1 P * ←-2 (z 1 -2) 2 + 1 (z 2 + 2) 2 , the refine problem 2 writes
	min z 1 ,z 2 ,p	J(z)
	s.t. c(ξ, z, p) ≤ 0, p ≥ 0, p ≤ P
	with a scalar input ξ ≥ 0, and with a scalar negotiable parameter p ≥ 0. Consider four constraints s.t.
	(Const. 1) (Const. 2) (Const. 3) (Const. 4)	z 2 ≥ ξ z 2 ≥ z 1 z 2 ≤ 2 + p -z 1 z 2 ≤ 2 + z 1

* .

Using HEGO with the two latter problems, we observe that the first component of the optimal solution, i.e. z * 1 , has a non-monotonous behavior w.r.t. ξ, as illustrated in

Figure 5.4.

  The latter modeling falls into the context of robust optimization (see e.g.[START_REF] Ben-Tal | Robust Optimization[END_REF] Ch.1]). It cannot be used as-is, since it conveys an infinite number of constraints. Let us make two further assumptions in order to simplify (5.29)1. The set F eas is convex in z, and linearly influenced by p. For example, we assume the existence of a (possibly non-linear) convex function g

	modified into	and (5.28) must be
	z ∈ F eas (ξ, p) , F eas (ξ, p + ∆p) ̸ = ∅, ∀∆p ∈ M	(5.29a) (5.29b)

ξ and (possibly non-linear) matrix-valued maps H(.), A(.), b(.) and B(.) s.t.

F eas (ξ, p) = {z ∈ R Nz : g ξ (z) ≤ H(ξ).p and A(ξ).z = b(ξ) + B(ξ).p} 2.

The set M is convex, and assumed to have a finite number of extreme point.

In other words, there are K vectors ∆p i ∈ R nneg s.t.

M = ConvexHull(∆p 1 , . . . , ∆p K ).

Proposition 14. Under the two latter assumptions, there are elements z k s.t. z k ∈ F eas ξ, p + ∆p k , ∀k = 1, . . . , K if and only if conditions (5.29a) and (5.29b) are satisfied. Proof. The return implication of the equivalence is guaranteed by construction. To prove the direct implication, let us consider any arbitrary vector ∆p ∈ M. By construction of M, there are K non-negative scalars σ k s.t. σ k = 1 and ∆p = K k=1 σ k ∆p k . Let us denote by z k an element of F eas ξ, p + ∆p k and consider the vector z := K k=1 σ k ∆p k . Using the convexity of g ξ , one gets

Table 6 . 3 :

 63 List of pairs of input components dispersed in the charts below. For example, using the data from Table6.2, there are 81 × 81 = 6561 pairs of values of (∆z 0 , ∆v 0 z ) that have been considered.

	Usefulness parameter Overall Horizontal Vertical Negotiable
	∆α max	+++	+++	0
	∆a max nor	++	++	0
	(∆z f , ∆y f )	+++	+++	0
	∆h f	+	0	+

Table 6 . 4 :

 64 Performance summary. The qualitative scale goes from 0 (useless) to + + + (extremely useful).by HEGO correspond to sharper turns, but does not require to negotiate ∆a max nor by a lot. The corresponding heatmaps (i.e. the heatmaps representing P * 2 ) are in the Appendix.

Table 6 . 5 :

 65 Reference trajectories considered for the benchmark w.r.t. the vertical flight envelope. All four trajectories share the same reference time-of-flight tf .

	Label Figure	Trajectory type	Engine flow structure	Comment
	A	6.18	Vertical	Near Min-Max
	B	6.19	Vertical	Max-Min-Middle
	C	6.20	Planar	Max-Min-Middle	Reference trajectory of Section 6.2.
	D	6.21	Fully 3D	Max-Min-Max	Out-of-plane reference.

See[START_REF] Açikmeşe | G-FOLD: A Real-Time Implementable Fuel Optimal Large Divert Guidance Algorithm for Planetary Pinpoint Landing[END_REF] and the associated video: youtube.com/watch?v=BqXFzVVCSCU.

For trajectories, it is the ratio between maximum flight altitude and maximum downrange.

An initial guess is needed for the adjoint state or the input u, t f and numerical methods are sensitive to these guesses.

Mars atmosphere is often treated as a small disturbance, for it is more than 100 times thinner than Earth's.

The exact definition of problem PDG (ξ, p), presented here in a simplified shape, is in Chapter 2.

Picking three values to describe the wind profile is an arbitrary design choice. If it is needed to consider a finer description of the wind profile in practice, significantly increasing the number of variables describing the wind profile does not change the approach presented in the next chapters, as discussed in Chapter 6.

One gimbal angle for the 2D model, two on the

3D model

In the sense of the successive rotations defining the rocket body frame.

If the rocket reaches v h = 0 before h = 0, the engine must be stopped, since hovering is prevented by the thrust dominance assumption. Thus, seeking v h (t f ) = -ε f v leaves a small safety margin.

Both pointing and landing cone constraints are well discussed in[START_REF] Malyuta | Convex Optimization for Trajectory Generation[END_REF] Fig.19] for instance.

Note that subsonic, supersonic and hypersonic speeds are usually defined as below Mach 1, between Mach 1 and Mach 5, and above Mach 5.

The exact definition of ξ will be detailed later, in Equation (4.4).

The exact performance index used in this thesis is detailed in the examples of Chapter 4.

L ∞ denotes the set of essentially bounded measurable functions.

Note that the ignition time optimization is a different topic.

In the numerical examples of Chapter 4, 5 and 6, the model is normalized for numerical stability, whereas, in this chapter, it aims at simplifying the writing of the dynamics.

Here ⋆| i denotes the i th component of ⋆.

Equations being linear in λ, one can consider λ/λ0 instead of λ.

Here, coarse refers to discretization schemes such as Euler methods with few collocation points.

L ∞ ([0, 1], U) denotes the set of essentially bounded measurable functions.

Recall that the derivatives are denoted by putting the variable of differentiation as an index, e.g. Lz = ∂L/∂z.

x • y denotes the component-wise product of vectors x and y, a.k.a. the Hadamard product.

See[START_REF] Büskens | Sensitivity Analysis and Real-Time Optimization of Parametric Nonlinear Programming Problems[END_REF] Sec. 4.3] for a complete discussion on that aspect.

Note that this would not apply to landings on offshore platforms, for obvious reasons.

Naturally, p b ⪯e pa is equivalent to pa ⪰e p b .

I would like to thank Laurent Pfeiffer for its attention to detail and its rich feedback.

Mj denotes the j th row of M .

If necessary, more details on the role of Q(J) are provided in[START_REF] Mangasarian | Lipschitz Continuity of Solutions of Linear Inequalities, Programs and Complementarity Problems[END_REF] Sec. 3].

These performances have been demonstrated live during the thesis defense, with on-demand initial condition changes.

Remerciements

 [-4, 2], for the 2D rocket model from Example 3. Note that the negotiation maps from sub-Figure (d) are exactly the kind of negotiation structure that support the conjecture presented in Remark 32. The Lipschitz-continuity of the optimal cost function, stated in Corollary 3, can be observed on the same sub-Figure . Also, remark that the cost function variations convey several active set changes within the nominal guidance mode, that can be visualized through the slope changes. The reason why the cost function remains constant when the horizontal final positions are being negotiated is that the optimal trajectory is shifted purely horizontally in these cases. • Projected incidence α z (Figures 6.7 and 6.8),

Color Value

• Controlled engine flow q c (Figures 6.9 and 6.10),

• Heatmap of the first negotiation map: P * 1 (Figures 6.11 and 6.12), • Heatmap of the third negotiation map: P * 3 (Figures 6.13 and 6.14). For the figures having the input pairs as their x and y-axis, the black dot • conveys the origin (0, 0) (for example the emergency mode maps). Also, the reference trajectory x and the reference control ū are represented in plain black.

Additional data is provided in Appendix, regarding the various states for these dispersions, the heatmap of the optimal time-of-flight change ∆t f * , and the other negotiation maps.

Observations and comments

Variable combinations. Some variables combinations have a natural constructive or destructive behavior. The pair (∆z 0 , ∆v 0 z ) is the clearest one. Starting farther away from the landing site but with a greater horizontal speed (or closer but with a lower horizontal speed) is a typical constructive behavior that explains the green strip in Figure 6.3-(a).

Continuity. The negotiation maps ξ → P *

i (ξ) discussed in Chapter 5 are Lispchitz-continuous. The reader might think that some of the maps below would suggest the contrary, due to abrupt changes. It is only a matter of zoom and Lipschitz constants. The detailed example pictured in Figure 6.15 provides various zoom levels to see how the negotiation maps can vary.

Normal acceleration negotiation.

As indicated by the black line in Figure 6.16, the ratio between the maximum normal acceleration of the reference trajectory and the normal acceleration bound is close to 1, leaving little room for negotiations. However, it is not always necessary to negotiate ∆a max nor , as pictured in Figure 6.1. For a large number of values of the pair (∆z 0 , ∆y 0 ), the guidance trajectories provided

A.1 Optimization results

A.1.1 Duality gap

The results presented here are taken from D. P. Bertsekas' book [START_REF] Bertsekas | Nonlinear Programming[END_REF], which contains all the proofs of the theorems and lemmas presented below. Let f , g and h be functions defined over R n , where f has scalar values, g has m components, and h has r components. Let us define the primal optimization problem

where f * denotes its optimal value. An optimization problem is said to be feasible if there exists at least one element satisfying its constraints. An optimization problem is said to be finite if its optimal value exists and does not equal ±∞.

The Lagrangian of Problem (A.1) is defined by

where µ ∈ R m and λ ∈ R r .

Definition 7. Vectors µ * and λ * are said to be Lagrange multipliers for the primal Problem

The dual of Problem (A.1) is the following problem sup

Let us denote by q * the optimal value of Problem (A.2). Here, q * ∈ R ∪ {±∞}.

Theorem 10 (Weak Duality Theorem [START_REF] Bertsekas | Nonlinear Programming[END_REF]Prop. 5.1.3]). Assume that the primal Problem (A.1) is feasible, but not necessarily finite. Then:

The duality gap is defined as the following non-negative quantity

We say that there is no duality gap if ∆ gap = 0, and that there is a duality gap if ∆ gap > 0.

Let us introduce the matrices and vectors E, d, A, b of proper dimensions, assume that f is convex, and define the problem Assume that the primal Problem (A.4) is feasible and its optimal value f * is finite. Let also f be convex over R n . Then, there is no duality gap and there exists at least one Lagrange multiplier.

Lemma 12 (Existence of Primal Optimal Solutions of LPs [13, Lemma 5.2.1]).

Assume that Problem (A.4) is feasible and its optimal value f * is finite. Let also f be linear. Then, Problem (A.4) has at least one optimal solution.

A.1.2 Right-hand side sensitivity of Linear Programs

Consider the following primal LP in its standard form and its optimal value function v (viewed as a function of its right-hand side b) s.t.

By convention, if A ⊤ µ ≤ c is infeasible, then the optimal value of the latter problem is -∞.

• On the contrary, if the "sup" term equals +∞ , then the primal problem with right-hand side b + td is infeasible for any t > 0.

Equation (A.7), that we also refer to as Gauvin's formula in this manuscript, is a powerful tool to analyze the RHS sensitivity of LPs. It is important to remark that the maximization problem in (A.7) is deeply linked to the dual LP from (A.6), s.t.

Proposition 16 (Adapted from [5, Prop.

2.3]). There exists a closed interval [α, β]

(possibly empty) s.t.

A.2 Differential Equations

A.2.1 Comparison theorem

For n ≥ 1, a function F : I × X ⊂ R × R n → R n is said to be quasi-monotone increasing if, for every pair (t, x) and (t, v) in I × X and every i = 1, . . . , n, one gets

Lemma 13 (Adapted from [80, IX.2.6]). Let F be a continuously-differentiable, quasi-monotone increasing function and x :

A.2.2 Flow of Ordinary Differential Equations

The sensitivity computations described below are standard material in the literature. The results recalled below focus mainly on the computational aspects of the different derivatives of the flow. For further references, complementary approaches can be found in [18, [START_REF] Oniki | Comparative Dynamics (Sensitivity Analysis) in Optimal Control Theory[END_REF] for more applied methods. Useful material about the flow properties can also be found in [84, Sec.

4.5] (properties of the flow linked to Lie brackets).

Let us consider a dynamic function f : R n × R m × R nη → R n with state x, control u and parameter η, which defines the ODE ẋ = f (x, u, η).

where x is defined over [0, 1] by the following IVP

x(0) = x 0 .

The non-relevant inputs of Φ f may be omitted to alleviate the writing. For instance, if f depends only on its state, its flow will be denoted Φ f (t, x 0 ).

The propositions below present how to compute the derivatives of Φ f w.r.t. each of its variables. Basically, these derivatives are the results of Initial Value Problems (IVPs) involving the derivatives of f . Note, first, that the time derivative of Φ f is a direct consequence of its definition and, for any t ∈ [0, 1], it equals

A.2.2.1 Sensitivity w.r.t. the initial condition

The derivative of Φ f w.r.t. the initial condition is referred to as the state transition matrix (STM), or simply the transition matrix [START_REF] Bryson | Applied optimal control: optimization, estimation and control[END_REF].

Proposition 17. The derivative of Φ f w.r.t. the initial condition x 0 is ∂Φ f ∂x 0 (T, x 0 ) = M (T ) where the transition matrix M satisfies the matrix-IVP

When the context is clear enough, the following notation is sometimes used (e.g. in Chapter 3, or in [START_REF] Bryson | Applied optimal control: optimization, estimation and control[END_REF]):

where

The IVP of Proposition 19 can be derived from Proposition 17, by considering an extended state composde of the original state, and one state having null dynamics for each component of the parameter η. For further details on this, see e.g. [START_REF] Demailly | Analyse numérique et équations différentielles[END_REF]. Using this remark, the smoothness of the flow w.r.t. its parameters can be derived from Proposition 18, as stated in the corollary below.

Corollary 4. Pick any

A.2.2.3 Sensitivity w.r.t. the control Proposition 20 (Adapted from [18, Sec 3.2]). Consider

When the control depends on a vector parameter, the following corollary holds.

where z is defined by the matrix-IVP

where x(t) = Φ f (t, x 0 ; u(θ, .)).

A.2. Differential Equations 153

Since the variable u belongs to a functional space, it is more technical to state how smooth is the flow w.r.t. to the control. Hence the need for the following proposition, which formally describes to what extent the function z u,v from Equation (A.10) defines the first-order expansion of Φ f w.r.t. u. Proposition 21 (Adapted from [START_REF] Pytlak | Numerical methods for optimal control problems with state constraints[END_REF]Prop. 1.3]). Consider a bounded set Ω ⊂ R m and a vector x 0 ∈ R n . Let us assume that f : R n × R m → R n is differentiable, that f , ∂f ∂x and ∂f ∂u are continuous and that there exists K < ∞ s.t.

Then there exists a function ε : R *

where 1

• z u,v is defined in Equation (A.10).

A.2.2.4 Technical summary

Consider a parametric control t → u µ (t), continuously differentiable in t and µ. The following expansion holds

where ε is a function s.t. ε(⋆)/∥ ⋆ ∥ 2 tends to zero when ⋆ = (∆x 0 , ∆η, µ) tends to zero, and where A, B and C are matrix valued functions defined by the following IVPs

where x(t) = Φ f (t, x 0 , η; u 0 ).

1 Here, L ∞ ([0, 1], R m ) denotes the vector space of essentially bounded measurable functions, for the essential supremum norm: ∥u∥L∞ := inf{C ≥ 0 : ∥u(t)∥2 ≤ C a.e. on [0, 1]}. Similarly, L ∞ ([0, 1], R m ) denotes the vector space of essentially bounded measurable functions, for the 1norm:

Additional data

Résumé

Ce chapitre contient des données complémentaires concernant l'exemple détaillé dans la Section 6.2 du Chapitre 6. This chapter contains complementary data regarding the example detailed in Section 6.2 of Chapter 6. Are displayed in these extra charts:

• the height states with respect to the normalized time,

• the heatmap of the optimal time-of-flight change ∆t * f , • and the heatmap of the second and fourth negotiations (P * 2 and P * 4 ).

These charts are presented for each pairs of Table 6.

ABSTRACT

This thesis studies emergency Powered Descent Guidance (PDG) for reusable launchers, as an Optimal Control Problem in free final-time with constraints. For such a launcher, subject to strong aerodynamic effects and having limited maneuverability, we wish to perform « emergency » trajectory planning by relaxing some negotiable parameters, such as the incidence safety bound, the normal acceleration load, or the landing site location. To this end, a hierarchy between the parameters is introduced and an algorithm, Hierarchical Emergency Guidance Optimization (H.E.G.O.), is developed to enforce it. The algorithm consists of a finite sequence of negotiation Linear Programs, followed by a refinement Quadratic Program.

The rocket is modeled by eight states, and three controls. The flight parameters are the initial conditions of the rocket states and other parameters, such as the Engine Specific Impulse and the wind profile. The user-defined hierarchy is conveyed via a co-lexicographic order. The methodology is theoretically studied.

Among others, the Lipschitz-continuity of the guidance trajectory with respect to the input flight parameters is established. Extensive numerical results serve to quantify the performance and relevance of the methodology. KEYWORDS powered descent guidance, optimal control, NLP sensitivity, emergency guidance