# Two-phase Trajectory Planning Method for Robust Planetary Landing in a Sensor-equipped Area

Clara Leparoux, Bruno Hérissé, Frédéric Jean

► To cite this version:

Clara Leparoux, Bruno Hérissé, Frédéric Jean. Two-phase Trajectory Planning Method for Robust Planetary Landing in a Sensor-equipped Area. 2025 European Control Conference (ECC), Jun 2025, Thessaloniki, Greece. pp.1296-1301, 10.23919/ECC65951.2025.11186925. hal-04752772

HAL Id: hal-04752772

https://ensta.hal.science/hal-04752772v1

Submitted on 24 Oct 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Two-phase Trajectory Planning Method for Robust Planetary Landing in a Sensor-equipped Area

Clara Leparoux<fnref n="1" />,<fnref n="2" />, Bruno Hérissé<fnref n="1" /> and Frédéric Jean<fnref n="2" />

*Abstract*— This article addresses the planetary landing problem by considering uncertainties and leveraging the presence of a detection area where precise measurements are available. The flight consists of two distinct phases: the first phase, subject to a high level of uncertainties, and the second phase, during which the vehicle is feedback controlled to ensure precise landing. We propose a method to compute the optimal control for the initial phase, aiming to minimize fuel consumption for the entire trajectory while satisfying a probabilistic constraint that ensures the vehicle reaches the detection zone with a specified threshold.

## I. INTRODUCTION

Planetary landing problems have gained significant attention recently, driven by the development of reusable launch vehicles and new space exploration missions that demand extremely high landing accuracy. Achieving this level of precision requires careful trajectory planning that accounts for various factors such as parameter uncertainties, external perturbations, and measurement errors, all of which can deviate the vehicle from its intended path. Additionally, trajectory planning must often incorporate system constraints or mission-specific limitations. Fuel consumption, another critical consideration, necessitates to be optimized, so that landing trajectories are typically computed within an optimal control framework.

To address the uncertainties inherent in landing dynamics, several robust guidance methods have emerged. Some approaches specifically target uncertainties in initial conditions, such as the method proposed by [1], which minimizes the solution's sensitivity to these uncertainties. Additionnaly, [2] proposes a replanification method to account for innacurate modelling and initial state errors. Other approaches model uncertainties using stochastic differential equations, like the algorithm presented by [3], which computes control strategies for landing under constraints on both mean and covariance. Similarly, [4] introduces a method for drone trajectory planning that includes probabilistic constraints to avoid collisions.

However, many of these methods assume that uncertainties affect the entire trajectory. Under this assumption, landing accuracy is constrained by the magnitude of these uncertainties. In practice, precise measurements, such as those from differential GPS, become available when the

vehicle is near the target, allowing for greater accuracy. This lead us to the problem treated in this paper: two-phase trajectory planning for a vehicle system being subject to uncertainty until its reaches a designated "detection zone". Several multi-phase trajectory planning methods already exist, particularly for spacecraft landings. For instance, [5] introduces a method to plan consecutive phases of a return mission using a Toss-back strategy. Similarly, [6] optimizes the three phases of a planetary landing by reformulating the problem as a nonconvex quadratically constrained quadratic program. Then, [7] also tackles lunar landing trajectory planning, using Model Predictive Static Programming to optimize the different phases. Another approach, presented by [8], focuses on minimizing fuel consumption up to a certain altitude, then applying a suboptimal state feedback controller to achieve a zero-attitude landing. Additionally, [9] addresses a multi-phase consensus problem, applied to interception missions. Our approach differs by considering uncertainties during a portion of the trajectory.

This paper introduces a novel approach for solving the planetary landing problem by incorporating a detection zone and dividing the vehicle's dynamics into two distinct phases. The detection zone is defined to guarantee that, once the vehicle enters it, successful arrival at the target is ensured. In the first phase, the vehicle operates under uncertainties, modeled using stochastic differential equations. To address the high computational complexity of this model, we apply a statistical linearization method, enabling us to operate within a deterministic optimal control framework. Once the vehicle enters the detection zone, it gains access to precise measurements, marking the transition to the second phase. During this phase, the dynamics are governed by feedback control, so that the impact of uncertainties is considered negligible compared with the first phase. Our focus is on optimizing the control strategy for the initial phase while meeting a probabilistic constraint of the vehicle reaching the detection zone. Simultaneously, we aim to minimize the cost of the entire trajectory. Consequently, the selection of the second-phase controller plays a crucial role in determining the overall success and cost-efficiency of the landing. This aspect is thoroughly discussed in a dedicated section.

The remainder of this paper is structured as follows. Section II introduces the general problem formulation. Section III applies this framework to the landing of reusable launch vehicles. In Section IV, we focus on the cost computation for the second phase, followed by numerical results in Section V.

<footnote n="1">¹Clara Leparoux (clara.leparoux@ensta-paris.fr) and Bruno Hérissé (bruno.herisse@onera.fr) are with DTIS, ONERA, Université Paris-Saclay, F-91123 Palaiseau, France.</footnote>

<footnote n="2">² Clara Leparoux and Frédéric Jean (frederic.jean@ensta-paris.fr) are with UMA, ENSTA Paris, Institut Polytechnique de Paris, 91120 Palaiseau, France.</footnote>
## II. PROBLEM FORMULATION

We investigate the trajectory planning problem for a landing vehicle with the requirement that it reaches a target S with a probability exceeding a specified threshold p. We assume that there exists a zone Z surrounding the target, within which the vehicle has access to precise sensor measurements. These measurements enable a predefined feedback control law to be applied from a time $t_1$, which ensures that the vehicle enters zone Z. Thus, the vehicle's dynamics consist of two phases:

* Phase 1 (from 0 to $t_1$): the vehicle operates under uncertainty, and its dynamics are governed by a controlled stochastic differential equation:

$$dx_t = f(x_t, u(t))dt + g(x_t)dW_t \quad (1)$$

where $x_0$ is a fixed random variable representing the initial state, $u(t)$ is the control input, and $dW_t$ is the increment of a Wiener process modeling the uncertainties;

* Phase 2 (from $t_1$ to $t_f$): the vehicle has access to precise measurements, allowing it to operate in a closed-loop configuration. As a result, the impact of uncertainties on the dynamics becomes negligible. The vehicle is controlled by a predefined feedback law, providing the control input $u$, which may depend on the state at the start of the second phase, $x(t_1)$. The dynamics during this phase are described by the following deterministic differential equation:

$$\dot{x}(t) = f(x(t), u_2(x(t))), \quad (2)$$

where $u_2(x) = u(x, x(t_1))$ denotes the controller used in Phase 2.

We formulate the robust trajectory planning as an optimal control problem, where the objective is to minimize a total cost that includes contributions from both phases with a probability constraint that the vehicle enters the detection zone. Before formally defining this optimal control problem, several key elements need to be introduced.

a) **Definition of $t_1$**: A first approach would be to define $t_1$ as the first time the vehicle's state enters the zone Z, triggering the transition from the dynamics in (1) to those in (2). In this scenario, $t_1$ would inherently depend on the randomness in the system, since it would represent a stopping time for the stochastic process $x_t$ governed by (1). However, in this work, we take a different approach: we treat $t_1$ as a pre-planned, deterministic time, independent of the randomness in the system. This allows us to formulate the guidance problem within a deterministic framework, as explained later in this section of the article, and to consider $t_1$ as a free optimization parameter. However, this requires certain assumptions about the uncertainties to ensure that the second phase does not start too late. To ensure the vehicle reaches the precise measurement zone Z with sufficient likelihood, we impose the constraint that the probability of being in Z at $t_1$ exceeds a given threshold p.

This constraint is expressed as

$$Pr[x_{t_1} \in Z] > p, \quad (3)$$

where $x_t$ represents the trajectory of the system governed by (1) under the control law $u(\cdot)$.

b) **Definition of Z**: The zone where precise measurements are available is denoted by the set Z. This set must necessarily include the target S. Additionally, to guarantee that once the vehicle enters Z it remains there, and that the deterministic dynamics given by (2) govern the system for all $t \ge t_1$, we impose the following assumption.

**Assumption 1:** The set Z is invariant under the dynamics described by (2).

By the definition of $t_1$, we ensure that the vehicle enters Z with a probability greater than p. The invariance of Z then guarantees that the vehicle ultimately reaches the target with a probability exceeding p.

c) **Choice of both the second-phase controller and $t_f$**: The choice of the second-phase controller plays a key role, as it directly influences the properties of the zone Z through Assumption 1. In this section, we do not impose a specific form for the controller, as our goal is to provide a general framework. However, a detailed comparison of various controller options is presented in Section III, specifically within the context of the landing problem. Moreover, the definition of the final time $t_f$ is intrinsically linked to the controller design. Since the initial state of the second phase is a realization of the random variable $x_{t_1}$, the final time $t_f$ can take several forms:

* it could be infinite, as in the case of an infinite-horizon LQ controller;

* it may be a function of the random variable $x_{t_1}$, typically dependent on its mean and covariance;

* it could be left free, thus depending on the realization of $x_{t_1}$;

* or it could be fixed, meaning that $t_f$ is the same for all realizations of $x_{t_1}$. In that case, the cost of the second-phase controller depends on the choice of $t_f$.

d) **Cost expression**: We seek to minimize the total cost of the trajectory, expressed as the expectation

$$J = \min_{u, t_1} \mathbb{E}[C(u, t_1)]. \quad (4)$$

Typically, the cost C represents factors such as fuel consumption along the trajectory. It is often practical to break this cost into two components corresponding to the two phases, $C_1$ for the first phase and $C_2$ for the second:

$$C = C_1(u, t_1) + C_2(x_{t_1}). \quad (5)$$

Since the second-phase control is fixed after $t_1$, $C_2$ depends only on the final state $x_{t_1}$ of the first phase. The form (5) allows for different cost criteria in each phase, such as fuel consumption for the first phase and final error for the second. In such cases, weighting or normalization terms may be needed to make the quantities comparable.

The cost (4) is computed as an expectation over all possible trajectories under control u, including those that
might not reach the target. Note that in theory, only trajectories that reach the target should be considered, especially since the second-phase control $u_2$ and cost $C_2$ may not be defined for those that miss the zone at $t_1$. To account for this, the cost should be reformulated as a conditional expectation. That said, this refinement is not necessary here, as we extend the definition of the controller $u_2$ outside of $\mathcal{Z}$, assigning a higher cost to realizations that goes outside of $\mathcal{Z}$. Thus, trajectories failing to enter $\mathcal{Z}$ have minimal influence on the optimization, in particular when the probability threshold $p$ is close to 1. The choice of a common $t_1$ for all trajectories is crucial to justify this approximation.

Finally, we formulate the planning problem as follows.

### Problem 1 (Stochastic 2-Phase Planning Problem):

$$ \min_{u,t_1} \mathbb{E}[C_1(u, t_1)] + \mathbb{E}[C_2(x_{t_1})] $$

subject to the constraints:

$$ \left\{ \begin{array}{l} x_t, \ t \in [0, t_1], \text{ trajectory of (1),} \\ u(t) \in \mathcal{U} \text{ a.e. on } [0, t_1], \\ \Pr[x_{t_1} \in \mathcal{Z}] > p. \end{array} \right. $$

This is a stochastic open-loop optimal control problem on the first phase, with free time $t_1$, a final probability constraint, and a cost that includes a term distributed along the trajectory $(\mathbb{E}[C_1(u, t_1)])$ and a final term $(\mathbb{E}[C_2(x_{t_1})])$. This type of problem is generally very difficult to solve. That's why we use here an approach similar to the one presented in [10], which is based on statistical linearization. Thus, we approximate the distribution of $x_t$ by a Gaussian distribution with mean $m$ and covariance $P$, solutions of the deterministic system

$$ \left\{ \begin{array}{l} \dot{m} = f(m, u), \\ \dot{P} = D_x f(m, u) P + P D_x f(m, u)^\top + g(m) g(m)^\top. \end{array} \right. \qquad (6) $$

Let $\phi(t)$ be the Gaussian distribution associated with $(m(t), P(t))$. In this way, we replace the cost (4) by

$$ \bar{C} = \mathbb{E}^{\phi}[C_1(u, t_1)] + \mathbb{E}^{\phi}[C_2(x_{t_1})], $$

where the expectations of $C_1$ and $C_2$ are calculated with the Gaussian density associated with $(m, P)$. Similarly, we approximate the constraint (3) by $\Pr^{\phi}[x_{t_1} \in \mathcal{Z}] > p$, expressed in terms of $m(t_1)$ and $P(t_1)$. Assuming that the zone $\mathcal{Z}$ is defined by affine constraints, i.e., that

$$ \mathcal{Z} = \{x \mid a_i x \le c_i \text{ for } i = 1, \dots, n_c\}, $$

the constraint (3) becomes

$$ a_i^\top m(t_1) + \Psi^{-1}(p) \sqrt{a_i^\top P(t_1) a_i} \le c_i, \qquad (7) $$

for $i = 1, \dots, n_c$ and where $\Psi^{-1}$ denotes the inverse cumulative function associated to the Gaussian distribution.

Thus, we approximate the solution of Problem 1 by the solution of the following deterministic optimal control problem (see [10] for more details on error estimates in this approximation).

### Problem 2 (Deterministic 2-Phase Planning Problem):

$$ \min_{u,t_1} \mathbb{E}^{\phi}[C_1(u, t_1)] + \mathbb{E}^{\phi}[C_2(x_{t_1})] $$

subject to the constraints:

$$ \left\{ \begin{array}{l} (m, P)(\cdot) \text{ follows (6) } \forall t \in [0, t_1], \\ u(t) \in \mathcal{U} \text{ almost everywhere on } [0, t_1], \\ a_i^\top m(t_1) + \Psi^{-1}(p) \sqrt{a_i^\top P(t_1) a_i} \le c_i \text{ for } i = 1, \dots, n_c. \end{array} \right. $$

This is a standard optimal control problem with free time $t_1$, a constraint on the final state, and a cost that includes an integral term and a final term.

## III. APPLICATION TO THE LAUNCHER

We now apply the general approach from the previous section to the problem of landing a reusable launcher. For simplicity, we focus on a two-dimensional scenario, whose state vector is defined as $x = (r, v, \mu)$, where $r = (r_y, r_z) \in \mathbb{R}^2$ is the position, $v = (v_y, v_z) \in \mathbb{R}^2$ is the velocity, and $\mu \in \mathbb{R}$ is the mass. The deterministic dynamics are defined by

$$ f(x, u) = \begin{pmatrix} v \\ \frac{T}{\mu}u - g_a \\ -q\|u\| \end{pmatrix}, \qquad (8) $$

where $T$ represents the maximum magnitude of the thruster, $q$ the mass flow rate, $g_a = (0, g_0)$, $g_0$ the gravitational acceleration, and $u \in \mathbb{R}^2$ the control variable. The actuators being constrained in magnitude and direction, the set of admissible controls is defined by

$$ \mathcal{U} = \left\{ \begin{array}{l} u = (u_y, u_z) \in \mathbb{R}^2 : \begin{array}{l} u_{min} \le \|u\| \le u_{max} \\ \text{and } u_z \ge \|u\| \cos(\theta_{max}) \end{array} \end{array} \right\}, $$

where $u_{min}$ and $u_{max}$ represent the minimum and the maximum allowable control magnitude, and $\theta_{max}$ defines the limit on the direction of the actuator output.

One of the objective while planning the trajectory is to minimize fuel consumption. Accordingly, the costs for the first and second phases are defined as follows:

$$ C_1 = \int_0^{t_1} \|u(t)\| dt, $$

and $C_2$ represents fuel consumption or a comparable quantity that needs to be specified. Finally the target $\mathcal{S}$ is defined as the set of states where $r = v = 0$.

### A. Modeling of the First Phase

From 0 to $t_1$, the vehicle is subject to uncertainties. It can measure its position and velocity, allowing for partial state feedback; however, these measurements are subject to noise, which must be incorporated into the modeling using the stochastic differential equation (1). We consider as control variable $\omega = (\rho, \theta, K_n, K_d)$ and express the partial state
feedback control as

$$u_{FB}(x, \omega) = (\rho + K_n \bar{x}) \begin{pmatrix} \cos(\theta + K_d \bar{x}) \\ \sin(\theta + K_d \bar{x}) \end{pmatrix},$$

where $\rho, \theta \in \mathbb{R}$, $K_n, K_d \in \mathbb{R}^4$, and $\bar{x} = (r, v)$. In the following, we write $u_\rho = \rho + K_n \bar{x}$ and $u_\theta = \theta + K_d \bar{x}$. Furthermore, we account for actuator limitations through saturation in the dynamics, following the method described in [10]. Thus we express the unperturbed dynamics of the launcher as $f_{FB}^{\text{sat}}(x, \omega) = f^{\text{sat}}(x, u_{FB}(x, \omega))$, where the saturated dynamics are given by

$$f^{\text{sat}}(x, u) = f \left( x, \text{sat}_{u_{\text{min}}}^{u_{\text{max}}}(u_{\rho}) \begin{pmatrix} \cos(\text{sat}\frac{\pi}{2} + \theta_{\text{max}}(u_{\theta})) \\ \sin(\text{sat}\frac{\pi}{2} + \theta_{\text{max}}(u_{\theta})) \end{pmatrix} \right). \quad (9)$$

## B. Choice of the Controller for the Second Phase

During the second phase, the vehicle benefits from accurate measurements. We model its closed-loop dynamics using a deterministic differential equation, which is derived by substituting $u$ in (8) with a closed-loop control $u_2(x)$,

$$\dot{x} = \begin{pmatrix} \frac{v}{\mu} & u_2(x) - g_a \\ -q \|u_2(x)\| & \end{pmatrix}. \quad (10)$$

Additionally, $u_2$ must make it possible to reach the target $S$ at time $t_f$, starting from the initial conditions $(r_1, v_1, \mu_1)$ at $t_1$. There are several approaches to design a closed-loop controller to solve this problem. For instance, an LQ controller or a ZEM/ZEV controller [11] could be used, or an MPC method. To ensure Assumption 1, it would be relevant to use state feedback associated with a Lyapunov function. We emphasize that an explicit formulation of $u_2$ is not required at this stage, our primary goal is to compute the cost of the trajectory.

## C. Definition of Zone Z

At $t_1$, the vehicle's state must reach an invariant zone $Z$ containing $S$, as per Assumption 1. Thus, we impose the following constraints to $Z$.

1) The position $r$ must belong to a zone $Z_r$, containing 0. We define this zone as limited by a maximum altitude $h_{max}$ and a glide-slope constraint of angle $\gamma$, i.e.,

$$Z_r \subseteq \{(r_y, r_z) \mid 0 \le r_z \le h_{max} \text{ and } r_z \ge \tan(\gamma) |r_y|\}.$$

2) The velocity must be directed towards the ground and inside the zone, and each component of the velocity vector should have a norm lower than a maximum speed $v_{max}$. The purpose of this constraint is to empirically ensure that the zone $Z$ remains invariant under the dynamics of the second phase. This will be verified numerically during simulations.

The zone $Z$ is then defined as follows:

$$Z = \{(r, v) \mid r \in Z_r, \begin{array}{l} -v_{max} \le v_y \le v_{max}, \\ -v_{max} \le v_z \le 0 \text{ and } r_y v_y \le 0 \end{array}\}.$$

The next step is to approximate $\text{Pr}[x \in Z] > p$ by a set of affine constraints. To simplify, we approximate the constraints on $r$ with the constraint $\text{Pr}^{\phi}[r \in Z_r] > p$. The latter is equivalent to the constraints on $m$ and $P$, derived from expression (7). This allows us to replace the constraint $\text{Pr}^{\phi}[x_{t_1} \in Z] > p$ by constraints on $(m, P)(t_1)$.

Ideally, we aim to ensure that there is enough fuel at time $t_1$ to carry out the second phase. This could be done by adding a constraint on the mass at $t_1$, or by verifying afterward that it is satisfied.

## D. Constraint on $t_1$

The constraint $\text{Pr}[x_{t_1} \in Z] > p$ implies that the confidence ellipse defined by $m$ and $P$ for a given confidence level $p$ is entirely contained within $Z$ at $t_1$. We tighten this constraint here in order to ensure that this condition holds at $t_1$, but also that $t_1$ represents the moment when we enter the zone $Z$ with a probability $p$. This amounts to asking that the ellipse has at least one point of contact with one of the edges of $Z$. In practice, for implementation purposes, we will simplify this constraint by requiring that only the top of the ellipse, i.e. $P_{zz}$, touches the upper boundary of $Z$ at time $t_1$. This is expressed mathematically as

$$m_{r_z}(t_1) = h_{max} - \Psi^{-1}(p) \sqrt{P_{zz}(t_1)}.$$

In practice, this simplification is justified if $P_{zz}(t_1)$ is small, i.e., if all trajectories are close at $t_1$.

## IV. ESTIMATION OF THE COST OF THE SECOND PHASE

To solve Problem 2, we need to compute $\mathbb{E}^{\phi}[C_2]$, which depends on how we design the controller $u_2$. We propose two methods for doing this. The first consists in choosing $u_2$ and $C_2$ from an LQ problem, and the second in determining them numerically.

### A. Exact expression of the value function of an LQ problem

Since the second phase is short, the vehicle's mass varies little. Assuming the mass remains constant, the dynamics (10) become control-affine. Consider the control variable $a \in \mathbb{R}^2$ defined by $a = \frac{T}{\mu}u - (0, g_0)$. The vehicle's dynamics is then linear and can be expressed as

$$\begin{cases} \dot{r} = v, \\ \dot{v} = a. \end{cases} \quad (11)$$

Choosing a control that solves an LQ problem is efficient and provides a state feedback. In this case, in Problem 2, we choose the cost of the second phase as the corresponding quadratic cost, with a normalization coefficient:

$$C_2 = \bar{x}(t_f)^{\top} Q_f \bar{x}(t_f) + \int_{t_1}^{t_f} a(t)^{\top} R a(t) dt, \quad (12)$$

where $R \in \mathbb{R}^2$ and $Q_f \in \mathbb{R}^4$ are weighting matrices. There exists an explicit expression for the finite-time value function associated with the dynamics (11) and the cost (12) for deterministic initial conditions at $t_1$, which is

$$C_{2, \det}(\bar{x}(t_1)) = \bar{x}(t_1)^{\top} S(0) \bar{x}(t_1)$$
where $S$ is the solution of the Riccati differential equation

$$ egin{cases} \frac{dS}{dt}(t) + A(t)^{\top}S(t) + S(t)A(t) \\ -S(t)B(t)R^{-1}B(t)^{\top}S(t) = 0, \ t \in [0, t_f - t_1], \\ S(t_f - t_1) = Q_f. \end{cases} $$

It follows that for $\bar{x}(t_1)$ stochastic with mean $\bar{m}(t_1)$ and covariance $\bar{P}(t_1)$, the expression for the expected value of the value function is

$$ \mathbb{E}^{\phi}[C_2(\bar{x}(t_1))] = \bar{m}(t_1)^{\top} S(0) \bar{m}(t_1) + \text{Tr}(S(0) \bar{P}(t_1)). $$

This explicit expression allows for very efficient computation of the expected cost of the second phase. However, it relies on the assumption that the mass is constant or varies little and does not account for constraints on control or state in the calculations. Moreover, it requires imposing $t_f$ and adjusting the normalization of $C_2$ with respect to mass consumption. A discussion on how to choose $t_f$ is provided in Section IV-C.

## B. Calculation by Numerical Solution of an Optimal Control Problem

We now consider a case where $u_2$ is only given numerically as a function of the initial condition $x(t_1)$. Consequently, we must numerically compute the expectation $\mathbb{E}^{\phi}[C_2]$. This can be done, for example, using a Monte Carlo method by calculating $C_2(x(t_1))$ for initial conditions $x(t_1)$ randomly generated according to a Gaussian distribution with mean $m(t_1)$ and covariance $P(t_1)$. The expectation of $C_2$ is then estimated by the average of the obtained costs. This method is flexible as it does not require assumptions about the dynamics. However, it is computationally intensive as it requires a large number of solutions to obtain accurate estimates.

We propose another approach using an "unscented" filter, originally introduced in [12]. The principle is similar to a Monte Carlo method, with the difference that the initial conditions are not randomly generated but chosen based on $m(t_1)$ and $P(t_1)$ to accurately represent the distribution of $x(t_1)$. Using this method, it is only necessary to generate $2n + 1$ solutions to estimate $\mathbb{E}^{\phi}[C_2]$, where $n$ is the dimension of the state. The initial conditions are called sigma points and are chosen as follows.

First, a sigma point is placed at the mean: $x_1^0 = m(t_1)$. Then, for $k$ ranging from 1 to $2n$,

$$ x_1^k = x_1^0 \pm \left( \sqrt{\frac{n}{1 - W^0} P(t_1)} \right)_k, $$

where the notation $(\cdot)_k$ denotes here the $k$-th column of the matrix. Finally, an estimate of the expectation of the cost is obtained by

$$ \mathbb{E}^{\phi}[C_2] \approx \sum_{k=0}^{2n} W^k C_{2,x_1^k}, $$

where $C_{2,x_1^k}$ is the cost of the second phase for the initial condition $x_1^k$, $W^0$ is fixed, and $W^k = \frac{1-W^0}{2n}$ for $k$ ranging from 1 to $2n$. We will use this method in a case where we

<table><caption>TABLE I: Parameters used for numerical results</caption><thead><tr><th>Property</th><th>Value</th><th>Property</th><th>Value</th></tr></thead><tbody><tr><td>T</td><td>1e6N</td><td>r<sup>0</sup></td><td>(1000, 4000)m</td></tr><tr><td>q</td><td>300kgs<sup>-1</sup></td><td>v<sup>0</sup></td><td>(-75, -200)ms<sup>-1</sup></td></tr><tr><td>u<sub>min</sub></td><td>0.2</td><td>μ<sup>0</sup></td><td>40000kg</td></tr><tr><td>u<sub>max</sub></td><td>0.8</td><td>g</td><td>diag(0, 0, 100, 10, 0)N</td></tr><tr><td>θ<sub>max</sub></td><td>65°</td><td>P<sup>0</sup></td><td>diag(100m<sup>2</sup>, 100m<sup>2</sup>, 1m<sup>2</sup>s<sup>-2</sup>, 1m<sup>2</sup>s<sup>-2</sup>, 1600kg<sup>2</sup>)</td></tr><tr><td>g<sub>0</sub></td><td>9.81ms<sup>-2</sup></td><td>γ</td><td>45°</td></tr><tr><td>h<sub>max</sub></td><td>600m</td><td>p</td><td>0.95</td></tr><tr><td>v<sub>max</sub></td><td>30m.s<sup>-1</sup></td><td></td><td></td></tr></tbody></table>

assume that $u_2$ is obtained by solving an LQ problem with additional constraints (or another optimal control problem).

### C. Choice of $t_f$

In the two methods presented to construct $u_2$, it is necessary to set the value of $t_f$ based on $m(t_1)$ and $P(t_1)$. There is no universal method to do this, but heuristics exist. For example, [8] proposes a calculation approximating the duration of the second phase to that of a gravity turn, whose expression is explicit:

$$ t_f = k \frac{\|m_v(t_1)\|}{2} \left( \frac{1 + \sin(\Theta(t_1))}{\alpha + g_0} + \frac{1 - \sin(\Theta(t_1))}{\alpha - g_0} \right), \quad (13) $$

where $\Theta(t_1)$ is the initial angle of the trajectory. In this expression, $\alpha$ corresponds to the acceleration value allowing for a gravity turn with final velocity and altitude being zero, and it is the solution of the following equation:

$$ \alpha^2 + \left( \frac{\|m_v(t_1)\|^2 \sin(\Theta(t_1))}{2m_{rz}(t_1)} \right) \alpha - \left( \frac{\|m_v(t_1)\|^2 g_0 (1 + \sin(\Theta(t_1)))^2}{4m_{rz}(t_1)} + g_0^2 \right) = 0. $$

The value of $t_f$ can also be chosen based on known numerical values.

## V. NUMERICAL RESULTS

In this section, we present numerical results illustrating the two-phase method and using the sigma-point method to calculate $\mathbb{E}^{\phi}[C_2]$ as presented. The results were obtained with the parameters and initial conditions detailed in Table I. The zone $\mathcal{Z}$ is defined with a maximum altitude $h_{max}$ of $600m$, a maximum speed $v_{max}$ of $30m.s^{-1}$. It is represented in gray in the figures. We first computed the solution of Problem 2 using the IPOPT library in Python. To regularize the problem, we penalize the state feedback gains in the cost of Problem 2, so that it is written as:

$$ \mathbb{E}^{\phi}[C_1(u, t_1)] + \mathbb{E}^{\phi}[C_2(x_{t_1})] + \int_{0}^{t_1} 2\|K_n(t)\|^2 + \|K_d(t)\|^2 dt. $$

The solution is used to generate trajectories following the dynamics (9) from 0 to $t_1$, with random initial conditions at 0 with mean $(r^0, v^0, \mu^0)$ and covariance $P^0$. The control and trajectory of the first phase are plotted in blue. From $t_1$
<figure><img src="image_2.png" /><figcaption>Fig. 1: Solution and trajectory samples</figcaption></figure>

to $t_f$, we calculate trajectories as solutions to the following problem.

Problem 3 (Simplified computation of phase 2):

$$ \min_{a} \int_{t_1}^{t_f} a(t)^{\top} a(t) dt + 200 \bar{x}(t_f)^{\top} \bar{x}(t_f) $$

subject to the constraints:

$$ \begin{cases} (r, v, \mu)(t_1) = (r_1, v_1, \mu_1), \\ \bar{x}(\cdot) \text{ follows (11)} \quad \forall t \in [t_1, t_f], \\ u_2(t) = (a(t) - g_a) \frac{\mu}{T} \in \mathcal{U} \text{ a.e. on } [t_1, t_f]. \end{cases} $$

We recall that $\bar{x} = (r, v)$. Additionally, $t_f$ is fixed to a value calculated by (13) from the initial conditions $(m, P)(t_1)$. The control and trajectory samples of the second phase are plotted in red.

We present results for which we numerically compute the term $E^{\phi}[C_2]$ by taking $m(t_1)$ as the initial condition at $t_1$. To do this, we computed $u_2$ using CasADi [13] with the IPOPT solver, and then we compute the contribution $C_2$ as the fuel consumption for the obtained $u_2$ controller. The trajectories generated by the obtained solution are plotted in Figure 1. We observe that trajectories arrive near the center of the zone at $t_1$, i.e., near $r_y = 0$, and trajectories obtained in the second phase are more vertical than those that would be obtained by a one-phase method [10]. Moreover, during the first phase, the control in Figure

1 has a Min-Max-like shape, which is fuel optimal for the unperturbed problem [14], but it is smoother and slightly away from the bounds $u_{min}$ and $u_{max}$. This comes from the saturation modeling in the dynamics. The control of the second phase is far from the bounds $u_{min}$ and $u_{max}$ and does not correspond to a Max-Min-Max shape, providing greater maneuvering margin near the landing point to correct the effects of disturbances.

## VI. CONCLUSIONS

We proposed a bi-level method for computing robust trajectories in a landing problem that involves a sensor-equipped zone. After detailing the process to obtain a general formulation of the problem, we applied it to the landing of a reusable launcher. Two distinct approaches were developed to compute the cost of the second phase of the trajectory. The first approach utilizes analytical calculations, which makes it numerically efficient. However, this method relies on some assumptions and does not focus on minimizing fuel consumption. The second approach is numerical, allowing for a more realistic modeling of the landing problem. The numerical results obtained with the second approach demonstrate the relevance of the method.

## REFERENCES

[1] H. Shen, H. Seywald, and R. W. Powell, "Desensitizing the minimum-fuel powered descent for mars pinpoint landing," *Journal of Guidance, Control, and Dynamics*, vol. 33, pp. 108-115, Jan. 2010.

[2] H. Menou, E. Bourgeois, and N. Petit, "Sensitivity analysis for powered descent guidance: Overcoming degeneracy," in *2022 European Control Conference (ECC)*, vol. 5, pp. 1218-1223, July 2022.

[3] J. Ridderhof and P. Tsiotras, "Minimum-fuel powered descent in the presence of random disturbances," Jan. 2019. AIAA Scitech 2019 Forum.

[4] T. Lew, F. Lyck, and G. Muller, "Chance-constrained optimal altitude control of a rocket," in *EUCASS 2019*, July 2019.

[5] E. Brendel, B. Hérissé, and E. Bourgeois, "Optimal guidance for toss back concepts of Reusable Launch Vehicles," in *EUCASS 2019*, July 2019.

[6] C. Wan, G. Jing, R. Dai, and J. R. Rea, "Fuel-optimal guidance for end-to-end human-mars entry, powered-descent, and landing mission," in *AIAA Scitech 2020 Forum*, American Institute of Aeronautics and Astronautics, Jan. 2020.

[7] A. Banerjee and R. Padhi, "Multi-phase MPSP guidance for lunar soft landing," *Transactions of the Indian National Academy of Engineering*, vol. 5, pp. 61-74, Mar. 2020.

[8] P. Lu, "Propellant-optimal powered descent guidance," *Journal of Guidance, Control, and Dynamics*, vol. 41, no. 4, pp. 813-826, 2018.

[9] X. Ai, L. Wang, J. Yu, and Y. Shen, "Field-of-view constrained two-stage guidance law design for three-dimensional salvo attack of multiple missiles via an optimal control approach," *Aerospace Science and Technology*, vol. 85, pp. 334-346, Feb. 2019.

[10] C. Leparoux, R. Bonalli, B. Hérissé, and F. Jean, "Statistical Linearization for Robust Motion Planning," *Systems & Control Letters*, vol. 189, p. 105825, 2024.

[11] H. Xiangyu, C. Xu, J. Hu, I. Maodeng, M. Guo, X. Wang, Y. Zhao, B. Hua, and Y. Wang, "Powered-descent landing gnc system design and flight results for tianwen-1 mission," *Astrodynamics*, vol. 6, pp. 3-16, Mar. 2022.

[12] S. Julier, J. Uhlmann, and H. Durrant-Whyte, "A new approach for filtering nonlinear systems," in *ACC'95*, 1995.

[13] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, "CasADi – A software framework for nonlinear optimization and optimal control," *Mathematical Programming Computation*, vol. 11, no. 1, pp. 1-36, 2019.

[14] C. Leparoux, B. Hérissé, and F. Jean, "Structure of optimal control for planetary landing with control and state constraints," *ESAIM: COCV*, vol. 28, p. 67, 2022.