22 D2.IPB.12.x71527 73 Paris, France 18–22 September 2022 2022Dr. Boris
Benedikter

Boris Benedikter${}^{\text{a}\ast}$, ${}^{\text{a}}$Postdoctoral Researcher,
Department of Mechanical and Aerospace Engineering, Sapienza University of
Rome, Via Eudossiana 18, 00184, Rome, Italy<EMAIL_ADDRESS>
Alessandro Zavoli${}^{\text{b}}$, ${}^{\text{b}}$Assistant Professor,
Department of Mechanical and Aerospace Engineering, Sapienza University of
Rome, Via Eudossiana 18, 00184, Rome, Italy<EMAIL_ADDRESS>
Guido Colasurdo${}^{\text{c}}$, ${}^{\text{c}}$Full Professor, Department of
Mechanical and Aerospace Engineering, Sapienza University of Rome, Via
Eudossiana 18, 00184, Rome, Italy<EMAIL_ADDRESS>
Simone Pizzurro${}^{\text{d}}$, ${}^{\text{d}}$Research Fellow, Launchers and
Space Transportation Department, Italian Space Agency, Via del Politecnico
snc, 00133, Rome, Italy<EMAIL_ADDRESS>
Enrico Cavallini${}^{\text{e}}$${}^{\text{e}}$Head of Space Transportation
Programs Office, Space Transportation, Space Infrastructures, and In-Orbit
Servicing Department, Italian Space Agency, Via del Politecnico snc, 00133,
Rome, Italy<EMAIL_ADDRESS>
# Stochastic Control of Launch Vehicle Upper Stage with Minimum-Variance
Splash-Down

###### Abstract

This paper presents a novel synthesis method for designing an optimal and
robust guidance law for a non-throttleable upper stage of a launch vehicle,
using a convex approach. In the unperturbed scenario, a combination of
lossless and successive convexification techniques is employed to formulate
the guidance problem as a sequence of convex problems that yields the optimal
trajectory, to be used as a reference for the design of a feedback controller,
with little computational effort. Then, based on the reference state and
control, a stochastic optimal control problem is defined to find a closed-loop
control law that rejects random in-flight disturbance. The control is
parameterized as a multiplicative feedback law; thus, only the control
direction is regulated, while the magnitude corresponds to the nominal one,
enabling its use for solid rocket motors. The objective of the optimization is
to minimize the splash-down dispersion to ensure that the spent stage falls as
close as possible to the nominal point. Thanks to an original convexification
strategy, the stochastic optimal control problem can be solved in polynomial
time since it reduces to a semidefinite programming problem. Numerical results
assess the robustness of the stochastic controller and compare its performance
with a model predictive control algorithm via extensive Monte Carlo campaigns.

##  1 Introduction

Launch vehicle dynamics are subject to significant uncertainties due to
dispersions of the propulsion system performance, random variations of the
local environment, and hard-to-model external perturbations. Therefore, the
design of a guidance and control algorithm that can robustly ensure operation
under uncertainty is essential for the system’s in-flight autonomy. Besides
safety requirements, the algorithm should also be optimal, e.g., in terms of
propellant consumption, to maximize the mass-to-orbit, enhance the vehicle
responsiveness in off-nominal conditions, and reduce the mission cost.

The optimal control of uncertain systems is traditionally an uneasy task, as a
general approach to an uncertain optimal control problem (OCP) would involve
solving a dynamic programming (DP) problem over arbitrary feedback control
laws [1]. However, DP is computationally impractical for real-world
applications due to the curse of dimensionality. Differential dynamic
programming (DDP) can yield an optimal control policy with lower computational
complexity than DP by expanding the OCP around a reference solution [2], but
it requires additional effort to retrieve the analytical formulation of the
second-order expansion of the objective function, large disturbances may
deviate excessively the system from the reference solution, and the inclusion
of nonlinear constraints is generally nontrivial [3].

The interest in machine learning (ML) techniques grew significantly in the
last few years for real-time guidance applications thanks to the low
evaluation times and the high accuracy in function approximation of deep
neural networks (DNNs) [4]. In particular, a DNN can be trained to imitate an
expert behavior (hence, a robust control policy) through behavioral cloning
[5, 6]. Also, reinforcement learning was proposed to train DNNs via repeated
interactions with the environment (i.e., simulations of the stochastic
dynamics) to progressively learn a robust optimal control policy [7, 8, 9].
However, the lack of any theoretical guarantees on the robustness of the
policy and the extensive computational effort associated with the DNN training
process, still limit the application of ML methods to relatively simple
problems.

The potentially prohibitive computational cost of DP or DNN training can be
bypassed if leveraging a model predictive control (MPC) scheme. Indeed, MPC,
rather than relying on a feedback policy computed pre-flight, recursively
solves an OCP onboard to update the control signal based on real-time
measurements. The recursive update provides the algorithm with inherent
robustness to in-flight disturbances and model errors even without necessarily
accounting for uncertainty in the OCP. The ability to simultaneously maximize
system performance and enforce mission requirements made MPC extremely
appealing for aerospace guidance and control applications [10, 11, 12, 13].
However, the degree of robustness of an MPC algorithm depends on the update
frequency, which is limited by the efficiency of the onboard hardware in
solving the OCP.

In this respect, convexification techniques gained great popularity in the
last few years, as they allow for solving originally nonconvex problems with
highly efficient convex programming algorithms [14]. For instance, _lossless_
convexification can be used to formulate problems with a quite common class of
nonconvex control constraints as equivalent convex programs solvable in
polynomial time [15]. More general sources of nonconvexity can be tackled via
_successive_ convexification, which comes with less sound theoretical
guarantees, but proved to be successful in numerous aerospace problems,
including spacecraft rendezvous [16], rocket ascent [17, 18, 19] and landing
[20, 21], low-thrust transfers [22, 23], and atmospheric entry [24, 25].
However, to ensure convergence, the parameters of the successive
convexification algorithm must be carefully picked and some modifications to
the formulation are often necessary to prevent undesired phenomena, such as
artificial infeasibility or unboundedness, resulting in longer computational
times and possibly limiting the update frequency (hence, the robustness) of an
MPC architecture.

To endow the optimal control policy with sound robustness guarantees,
uncertainty must be included in the problem formulation. In this respect,
robust optimal control (ROC) methods account for uncertainty by modeling them
as deterministic variables bounded within finite ranges. In this way, the
associated optimization problem is solvable with traditional programming
algorithms. However, this modeling strategy usually leverages either min-max
OCP formulations [26] or constraint tightening [27], which often lead to
overly conservative policies. More recently, tube-based approaches have been
proposed to design more performing policies than traditional ROC methods.
Specifically, by accounting for (bounded) uncertainties directly in the
control algorithm, tube-based approaches can find a control policy that
encloses all future state realizations in a bounded region (referred to as a
_tube_) [28].

On the other hand, stochastic optimal control (SOC) deals with unbounded
uncertainties, which are the class of uncertainty that affect most real-world
systems, by incorporating their probabilistic descriptions in the OCP [29].
The unboundedness of the random variables requires formulating hard
requirements on the state and controls as relaxed probabilistic constraints,
called _chance constraints_ , which allow for a small probability of
violation. Chance-constrained optimization has been applied to several
aerospace problems, such as UAV path planning with obstacle avoidance
constraints [30], model rocket ascent [31], station-keeping in unstable orbits
[32], and impulsive [33] or finite-thrust [34] orbit transfers. However, most
works in the literature consider a simple open-loop control law that cannot
control the time-evolution of the covariance, which increases uncontrolled
over time due to diffusive stochastic perturbations.

To design performing solutions, the stochastic OCP should account for the
presence of a feedback controller. A simple approach consists in designing the
feedback controller _a priori_ via traditional control theory approaches
(e.g., Linear Quadratic Regulator, LQR) and then optimizing only the
feedforward control signal in the stochastic framework [35]. However, the
quality of the solution depends on the considered controller, which, if not
tuned properly, may still lead to over-conservative policies. _Covariance
control_ , originally proposed in the late 1980s for infinite-horizon problems
[36], gained renewed interest among the control community after being applied
to finite-horizon problems as well [37, 38]. The major appeal of this control
synthesis method lies in the ability to systematically design, with low
computational complexity, a closed-loop optimal control policy that drives a
stochastically perturbed system from an initial probability distribution (in
terms of mean state and covariance) to a desired one at a final time. Indeed,
for a linear system, the finite-horizon stochastic OCP can be cast as a
deterministic convex optimization problem, thus solvable in polynomial time
[39, 40]. Recently, thanks to convexification methods, covariance control has
been applied to nonlinear systems as well, including low-thrust interplanetary
trajectory design problem [41, 42], powered descent and landing [43, 44], and
planetary entry [45].

In most works on covariance control, the control law is parameterized as the
sum of a feedforward term and a linear feedback one. Thus, when operating in
uncertain environments, the system adjusts the control magnitude and direction
to compensate for deviations from the nominal path. However, such a control
policy can only be implemented by a system that can regulate the magnitude of
the control action. Instead, non-throttleable systems must rely on policies
that regulate exclusively the direction of the control vector, without
straying from the nominal control magnitude profile.

This paper presents an algorithm to design a nominal optimal trajectory for a
non-throttleable upper stage of a launch vehicle and a robust controller that
compensates for exogenous disturbances. The considered case study involves the
flight of VEGA’s third stage, which features a solid rocket motor (SRM), thus,
its thrust level cannot be regulated in flight. The system must be actively
controlled to ensure that the spent stage falls in an uninhabited area; this
requirement is generally represented as a _splash-down constraint_.

A multiplicative feedback control law that governs only the thrust vector
direction is considered so that the closed-loop control norm is independent of
the disturbances. Due to the nonlinearity arising from the coupling of the
feedback gains and the covariance matrix in the dynamics, the associated
optimization problem is not convex. Thus, to reduce the computational
complexity of the solution process, the problem is solved in two steps. First,
the optimal trajectory and the associated feedforward control signal are
computed in the deterministic scenario, that is, neglecting uncertainties, and
serve as a reference for the second step. Then, a covariance control problem,
based on the reference trajectory and feedforward control, is solved to find
the multiplicative feedback gains that minimize the dispersion of the spent
stage splash-down point. In either step, the OCP can be solved with low
computational complexity thanks to the introduction of convenient changes of
variables, lossless constraint relaxations, and successive linearization.

The paper is organized as follows. Section 2 describes the upper stage
guidance and control problem, detailing the phases, the dynamical model, and
the mission requirements of the corresponding OCP. Section 3 outlines the
convexification strategy used to cast the nonlinear OCP as a sequence of
convex problems that rapidly yields the nominal trajectory and control law in
the unperturbed scenario. Sections 4 and 5 outline two different strategies to
design a robust control framework for the upper stage flight. In Section 4, an
MPC algorithm that recursively solves the sequential convex optimization
problem is proposed as a means to provide robustness to model mismatches and
random disturbances. Instead, Section 5 outlines an optimal covariance control
problem that yields a multiplicative feedback controller robust to the
considered characterization of the in-flight perturbations. Section 6 presents
numerical results to assess the robustness and performance of the proposed
algorithms through extensive Monte Carlo campaigns and compares the MPC
algorithm with the stochastic control approach. A conclusion section ends the
manuscript.

##  2 Problem Statement

In this section, the upper stage optimal control problem is presented. First,
the phase sequence of the OCP is outlined; then, the considered dynamical
model is presented; and finally, the objective and all mission requirements of
the OCP are discussed.

###  2.1 Phase Sequence

The problem under investigation corresponds to the final portion of the VEGA
launch vehicle ascent flight. In particular, the operation of the first two
stages is neglected and the goal of the optimization is to design an optimal
control policy for the third stage flight. To evaluate and maximize the
vehicle performance in terms of carrying capacity the OCP horizon must extend
until payload release.

Figure 1: Phases of the optimal control problem.

Figure 1 illustrates schematically the flight phases of the VEGA launch
vehicle from the third stage ignition until payload injection into the target
orbit. The phases are numbered from 1 to 5 in chronological order, and Phase 6
simulates the uncontrolled return of the third stage after burnout. The third
stage flight corresponds to Phase 1 and terminates at the burnout of the SRM.
A short coasting (Phase 2) of fixed duration follows the third stage
separation. Then, the fourth and final stage must complete the orbit insertion
maneuver, whose optimal phase sequence is similar to a Hohmann maneuver, i.e.,
composed of two propelled arcs (Phases 3 and 5) separated by a quite long
coasting (Phase 4). The time-lengths of Phases 3–5 must be optimized. The
return (Phase 6) is included in the OCP to predict and constrain the splash-
down point of the spent third stage. Also the duration of Phase 6 must be
determined by the optimization process.

Hereinafter, let $t_{0}^{(i)}$ and $t_{f}^{(i)}$ denote the initial and final
time of the $i$-th phase and $\smash{\Delta t^{(i)}}$ its time-length. For the
sake of simplicity, if no phase superscript is specified, then $t_{0}$ and
$t_{f}$ denote the ignition time of the third stage $\smash{t_{0}^{(1{})}}$
and the fourth stage burnout $\smash{t_{f}^{(5{})}}$, respectively. Likewise,
let $t_{R}$ denote the splash-down time of the spent stage $t_{f}^{(6{})}$.

###  2.2 System Dynamics

The vehicle is modeled as a point mass subject to a 3-DoF translational
motion. The rotational dynamics are neglected and the rocket axis is assumed
to be aligned with the thrust direction at any time. The state $\bm{x}$ of the
rocket is described by its position $\bm{r}$, velocity $\bm{v}$, and mass $m$,
as $\bm{x}=[x\;y\;z\;v_{x}\;v_{y}\;v_{z}\;m]$. The position and velocity are
expressed in Cartesian Earth-Centered Inertial (ECI) coordinates. The $x$ and
$y$ axes lie in the equatorial plane, and the $z$ axis is aligned with the
Earth’s angular velocity $\bm{\omega}_{E}$, forming a right-hand frame. The
main advantage of using this set of state variables over, for example, a
spherical coordinate system [18], is the fact that any trajectory can be
studied, even missions to polar or high inclination orbits, which are the
typical targets for VEGA, without running into singularities.

The only external forces assumed to act on the launch vehicle are gravity, the
aerodynamic drag, and the motor thrust. A simple inverse-square gravitational
model is assumed, thus the gravity acceleration is $\bm{g}=-\mu\bm{r}/r^{3}$,
with $\mu$ denoting the Earth’s gravitational parameter. The drag force is
modeled as
$\bm{D}=-0.5\rho{v}_{\text{rel}}^{2}S_{\text{ref}}C_{D}\bm{\hat{v}}_{\text{rel}}$,
where $\rho$ is the atmospheric density, $S$ is the reference surface, $C_{D}$
is the drag coefficient, assumed to be constant, and
$\bm{v}_{\text{rel}}=\bm{v}-\bm{\omega}_{E}\times\bm{r}$ is the relative-to-
atmosphere velocity. Lift is neglected, as it is much smaller than drag for
launch vehicles.

VEGA’s third stage is an SRM, so the thrust history $T(t)$ and the associated
mass flow rate $\dot{m}_{e}(t)$ depend on the geometry of the solid
propellant. Although upper stages operate at high altitudes, the net thrust is
slightly different from the vacuum one due to the atmospheric pressure $p$
contribution, thus we consider a vacuum thrust law $T_{\text{vac}}(t)$ and a
corresponding mass flow rate law $\dot{m}_{e}(t)$. The net thrust is
$T=T_{\text{vac}}(t)-pA_{e}$, where $A_{e}$ is the nozzle exit area.

Instead, the last stage is equipped with a small liquid rocket engine, AVUM;
hence, its thrust and mass flow rate are assumed to be constant and equal to
the maximum attainable values. Also, AVUM can be cut off and re-ignited; so,
its overall burn time can be minimized by the optimization process to save as
much propellant as possible.

While the thrust magnitude is prescribed, the thrust direction $\bm{\hat{T}}$
is the control variable. Its components are expressed in ECI coordinates and,
being a unit vector, must satisfy the following identity

$\hat{T}_{x}^{2}+\hat{T}_{y}^{2}+\hat{T}_{z}^{2}=1$ (1)

The resulting equations of motion $\bm{\dot{x}}=\bm{f}(\bm{x},\bm{u},t)$ are

$\displaystyle\dot{\bm{r}}$ $\displaystyle=\bm{v}$ (2)
$\displaystyle\dot{\bm{v}}$
$\displaystyle=-\frac{\mu}{r^{3}}\bm{r}+\frac{T}{m}\bm{\hat{T}}-\frac{D}{m}\bm{\hat{v}}_{\text{rel}}$
(3) $\displaystyle\dot{m}$ $\displaystyle=-\dot{m}_{e}$ (4)

###  2.3 Mission Requirements and Objective

The state at the third stage ignition is supposed to be completely assigned,
that is

$\bm{x}(t_{0})=\tilde{\bm{x}}_{0}$ (5)

At the burnout of the fourth stage, the payload must be released into the
target orbit. In this work, we suppose that the target orbit is a circular
orbit of semi-major axis $a_{\text{des}}$ and inclination $i_{\text{des}}$.
Thus, the final conditions to impose at the end of Phase 5 are

$\displaystyle x(t_{f})^{2}+y(t_{f})^{2}+z(t_{f})^{2}$
$\displaystyle=a_{\text{des}}^{2}$ (6) $\displaystyle
v_{x}(t_{f})^{2}+v_{y}(t_{f})^{2}+v_{z}(t_{f})^{2}$
$\displaystyle=v_{\text{des}}^{2}$ (7)
$\displaystyle\bm{r}(t_{f})\cdot\bm{v}(t_{f})$ $\displaystyle=0$ (8)
$\displaystyle x(t_{f})v_{y}(t_{f})-y(t_{f})v_{x}(t_{f})$
$\displaystyle=h_{z,\text{des}}$ (9)

Equation (6) prescribes the final radius, while Eq. (7) with
$v_{\text{des}}=\sqrt{\mu/a_{\text{des}}}$ and Eq. (8) are sufficient to
impose the null eccentricity constraint. As for the inclination,
$i_{\text{des}}$ is imposed by prescribing the $z$-component of the angular
momentum vector $h_{z,\text{des}}=\cos i_{\text{des}}\sqrt{\mu
a_{\text{des}}}$.

After burnout, the third stage features a high, but suborbital, velocity, so
it falls back to Earth. To prevent it from crashing into an inhabited area,
Phase 6 simulates its return and the splash-down point can be constrained by
imposing the following conditions

$\displaystyle x(t_{R})^{2}+y(t_{R})^{2}+z(t_{R})^{2}$
$\displaystyle=R_{E}^{2}$ (10) $\displaystyle z(t_{R})$
$\displaystyle=z_{R,\text{des}}$ (11)

with $z_{R,\text{des}}=R_{E}\sin\varphi_{R,\text{des}}$. Equation (10)
constrains the final radius to be equal to the Earth radius $R_{E}$, while Eq.
(11) constrains the latitude of the splash-down point to a desired value
$\varphi_{R,\text{des}}$. For missions toward polar or quasi-polar orbits
(e.g., Sun-synchronous orbits), which are the scope of the present paper,
constraining the latitude is equivalent to constraining the distance from the
launch site, which is the only quantity that can be actively controlled during
the third stage flight. Indeed, constraining the longitude to an arbitrary
value is impractical, as it would be very expensive because it requires a
maneuver to change the plane of the ascent trajectory.

The multi-phase structure of the problem can account for mass discontinuities,
which take place at the third stage separation. In particular, all state
variables are continuous across phases, while the following linkage conditions
must be imposed

$m(t_{0}^{(2{})})=m(t_{f}^{(1{})})-m_{\text{dry},3}$ (12)

As for the return phase, the initial position and velocity are equal to the
burnout conditions of the third stage, while its mass corresponds to the dry
mass of the stage. Thus,

$\displaystyle\bm{r}(t_{0}^{(6{})})$ $\displaystyle=\bm{r}(t_{f}^{(1{})})$
(13) $\displaystyle\bm{v}(t_{0}^{(6{})})$
$\displaystyle=\bm{v}(t_{f}^{(1{})})$ (14) $\displaystyle m(t_{0}^{(6{})})$
$\displaystyle=m_{\text{dry},3}$ (15)

Assuming that the fairing is jettisoned during the previous phases of the
flight, a constraint on the maximum heat flux that the payload undergoes must
be included in the formulation from Phase 1 to Phase 5. Specifically, modeling
the heat flux as a free molecular flow acting on a plane surface perpendicular
to the relative velocity [46], the constraint is expressed as

$\dot{Q}=\frac{1}{2}\rho v_{\text{rel}}^{3}\leq\dot{Q}_{\text{max}}$ (16)

The goal of the optimization is to minimize the propellant consumed during the
fourth stage flight. Indeed, if the guidance algorithm minimizes the
propellant mass consumption, not only the system will feature better
responsiveness in off-nominal conditions, but it will also be able to carry
additional payload mass, as it requires smaller fuel tanks. The objective can
be equivalently formulated as maximizing the final mass. Thus, the OCP that
must be solved to retrieve the optimal trajectory and control in the
unperturbed scenario is

$\displaystyle\mathcal{P}_{A}:\;\;\min_{\bm{x},\bm{\hat{T}},\Delta t}$
$\displaystyle\quad-m(t_{f})$ (17a) s.t. (1)–(16) (17b)

##  3 Convex Formulation of the Deterministic OCP

This section outlines the convexification strategy devised to convert the
deterministic OCP in Eq. (17) into a sequence of convex problems that quickly
converges to the optimal solution.

###  3.1 Change of Variables

The coupling of the control variable $\bm{\hat{T}}$ with the state variable
$m$ in Eqs. (2)–(4) may lead to numerical issues, such as high-frequency
jitters in the control law [47]. Therefore by introducing a new control
variable

$\bm{u}=\frac{T}{m}\bm{\hat{T}}$ (18)

the dynamics are control-affine, that is, they can be written as

$\bm{f}=\tilde{\bm{f}}(\bm{x},t)+\tilde{B}\bm{u}$ (19)

By replacing Eq. (18) in Eqs. (2)–(4), the equations of motion become

$\displaystyle\dot{\bm{r}}$ $\displaystyle=\bm{v}$ (20)
$\displaystyle\dot{\bm{v}}$
$\displaystyle=-\frac{\mu}{r^{3}}\bm{r}+\bm{u}-\frac{D}{m}\bm{\hat{v}}_{\text{rel}}$
(21) $\displaystyle\dot{m}$ $\displaystyle=-\dot{m}_{e}$ (22)

and the $\tilde{B}$ matrix in Eq. (19) is $\tilde{B}=[\bm{0}_{3\times
3}\;\bm{I}_{3\times 3}\;\bm{0}_{1\times 3}]^{T}$ with $\bm{0}_{m\times n}$ and
$\bm{I}_{m\times n}$ denoting the null and identity matrix of size $m\times
n$.

The new control variable must satisfy Eq. (1), which is now reformulated as

$u_{x}^{2}+u_{y}^{2}+u_{z}^{2}=u_{N}^{2}$ (23)

where $u_{N}$ is an additional variable defined as

$u_{N}=\frac{T}{m}$ (24)

In practice, Eq. (24) must be enforced as a nonlinear path constraint. This is
the price to pay to obtain control-affine dynamics, but the advantages of this
change of variables outweigh the disadvantages, as the complexity of the
problem is reduced compared to a coupling of state and control variables in
the dynamics.

###  3.2 Constraint Relaxation

Equation (23) is a nonlinear equality constraint that must be convexified
somehow. This constraint belongs to a class of constraints that is suitable
for a lossless relaxation [48, 15]. Specifically, by replacing the equality
sign with an inequality sign, a second-order cone constraint is obtained

$u_{x}^{2}+u_{y}^{2}+u_{z}^{2}\leq u_{N}^{2}$ (25)

This relaxation is lossless because, despite it defines a larger feasible set,
the solution of the convex problem is the same as the original and Eq. (25) is
satisfied with the equality sign. A theoretical proof of the lossless property
of this relaxation can be found in Proposition 1 of Ref. [19].

###  3.3 Successive Linearization

The remaining nonlinear expressions are replaced with first-order Taylor
series expansions around a reference solution that is updated at every
iteration.

####  3.3.1 Equations of Motion

To account for free-time phases in the optimization, we introduce an
independent variable transformation from physical time $t$ to a new variable
$\tau$ defined, for each phase, over a fixed unitary domain $[0,1]$. The
relationship between the two independent variables is

$t^{(i)}=t_{0}^{(i)}+(t_{f}^{(i)}-t_{0}^{(i)})\,\tau$ (26)

Note that the time dilation $s$ between $t$ and $\tau$ corresponds to the
time-length of the phase

$s^{(i)}=\frac{dt}{d\tau}=t_{f}^{(i)}-t_{0}^{(i)}$ (27)

The time-length $s^{(i)}$ is then introduced as an additional optimization
variable for every free-time phase.

The equations of motion (20)–(22) are control-affine but still nonlinear.
Thus, they are expressed in terms of $\tau$ and linearized around a reference
solution $\\{\bar{\bm{x}},\bar{\bm{u}},\bar{s}\\}$

$\bm{x}^{\prime}\vcentcolon=\frac{d\bm{x}}{d\tau}=s\bm{f}\left(\bm{x},\bm{u},\tau\right)\approx
A\bm{x}+B\bm{u}+\Sigma s+\bm{c}$ (28)

where the following matrices and vectors were introduced

$\displaystyle A$
$\displaystyle=\bar{s}\frac{\partial\bm{f}}{\partial\bm{x}}(\bar{\bm{x}},\bar{\bm{u}},\tau)$
(29) $\displaystyle B$
$\displaystyle=\bar{s}\frac{\partial\bm{f}}{\partial\bm{u}}(\bar{\bm{x}},\bar{\bm{u}},\tau)$
(30) $\displaystyle\Sigma$
$\displaystyle=\bm{f}(\bar{\bm{x}},\bar{\bm{u}},\tau)$ (31)
$\displaystyle\bm{c}$ $\displaystyle=-(A\bar{\bm{x}}+B\bar{\bm{u}})$ (32)

An undesired phenomenon associated with linearization is artificial
infeasibility, that is, the linearized problem may be infeasible. Therefore, a
virtual control signal $\bm{q}$ is added to the dynamics to guarantee the
feasibility of the convex problem, that is

$\bm{x}^{\prime}=A\bm{x}+B\bm{u}+\Sigma s+\bm{c}+\bm{q}$ (33)

Since virtual controls are unphysical variables, their use must be highly
penalized by adding the following penalty term to the cost function

$J_{q}=\lambda_{q}P(\bm{q})$ (34)

where $\lambda_{q}$ is a sufficiently large penalty weight and $P(\bm{q})$ a
suitable penalty function, which will be defined after discretization in
Section 3.5.

####  3.3.2 Boundary Constraints

The final conditions at the payload release (6)–(10) are linearized as

$\displaystyle\bar{\bm{r}}(t_{f})\cdot\bar{\bm{r}}(t_{f})+2\bar{\bm{r}}(t_{f})\cdot(\bm{r}(t_{f})-\bar{\bm{r}}(t_{f}))$
$\displaystyle=a_{\text{des}}^{2}$ (35)
$\displaystyle\bar{\bm{v}}(t_{f})\cdot\bar{\bm{v}}(t_{f})+2\bar{\bm{v}}(t_{f})\cdot(\bm{v}(t_{f})-\bar{\bm{v}}(t_{f}))$
$\displaystyle=v_{\text{des}}^{2}$ (36)
$\displaystyle\bar{\bm{r}}(t_{f})\cdot\bar{\bm{v}}(t_{f})+\bar{\bm{v}}(t_{f})\cdot(\bm{r}(t_{f})-\bar{\bm{r}}(t_{f}))$
$\displaystyle+\bar{\bm{r}}(t_{f})\cdot(\bm{v}(t_{f})-\bar{\bm{v}}(t_{f}))$
$\displaystyle=0$ (37)
$\displaystyle\bar{v}_{y}(t_{f})(x(t_{f})-\bar{x}(t_{f}))-\bar{v}_{x}(t_{f})(y(t_{f})-$
$\displaystyle\bar{y}(t_{f}))$
$\displaystyle-\bar{y}(t_{f})v_{x}(t_{f})+\bar{x}(t_{f})v_{y}(t_{f})$
$\displaystyle=h_{z,\text{des}}$ (38)

Likewise, the terminal radius constraint at the splash-down, Eq. (10), is
formulated as

$\bar{\bm{r}}(t_{f}^{(6{})})\cdot\bar{\bm{r}}(t_{f}^{(6{})})+2\bar{\bm{r}}(t_{f}^{(6{})})\cdot(\bm{r}(t_{f}^{(6{})})-\bar{\bm{r}}(t_{f}^{(6{})}))=R_{E}^{2}$
(39)

Analogously to the dynamics, also the linearization of these constraints may
cause artificial infeasibility. Thus, we introduce virtual buffer zones. In
particular, Eqs. (35)–(39) are grouped into a vector $\bm{\chi}=\bm{0}$ that
is then relaxed as $\bm{\chi}=\bm{w}$. The vector $\bm{w}$ holds all the
virtual buffers, which are highly penalized by adding the following penalty
term to the cost function

$J_{w}=\lambda_{w}\left\lVert\bm{w}\right\rVert_{1}$ (40)

with $\lambda_{w}$ denoting the penalty weight of the virtual buffers.

####  3.3.3 Path Constraints

The heat flux constraint (16) is another nonlinear expression that must be
linearized. By considering its Taylor series expansion, we obtain

$\dot{\bar{Q}}+\frac{\partial\dot{\bar{Q}}}{\partial\bm{r}}\cdot(\bm{r}-\bar{\bm{r}})+\frac{\partial\dot{\bar{Q}}}{\partial\bm{v}}\cdot(\bm{v}-\bar{\bm{v}})\leq\dot{Q}_{\text{max}}$
(41)

where

$\displaystyle\frac{\partial\dot{\bar{Q}}}{\partial\bm{r}}$
$\displaystyle=\frac{1}{2}\frac{d\bar{\rho}}{d\bm{r}}\bar{v}_{\text{rel}}^{3}+\frac{3}{2}\bar{\rho}\bar{v}_{\text{rel}}\bm{\omega}_{E}\times\bm{\bar{v}}_{\text{rel}}$
(42) $\displaystyle\frac{\partial\dot{\bar{Q}}}{\partial\bm{v}}$
$\displaystyle=\frac{3}{2}\bar{\rho}\bar{v}_{\text{rel}}\bm{\bar{v}}_{\text{rel}}$
(43)

Likewise, the auxiliary variable $u_{N}$ must satisfy Eq. (24) at every time.
This is a nonlinear path constraint and it is linearized as

$u_{N}=\frac{T_{vac}-\bar{p}A_{e}}{\bar{m}}\left(2-\frac{m}{\bar{m}}\right)-\frac{A_{e}}{\bar{m}}\frac{d\bar{p}}{d\bm{r}}\cdot(\bm{r}-\bar{\bm{r}})$
(44)

###  3.4 Trust Region on Time-Lengths

To enhance convergence, a soft trust region constraint is imposed on the time-
lengths $s$ of Phases 4 and 5, which resulted particularly sensitive to
diverging phenomena [19]. The trust region constraint is

$|s^{(i)}-\bar{s}^{(i)}|\leq\delta^{(i)}\qquad i=4{},5{}$ (45)

The trust radii $\delta^{(i)}$ are additional optimization variables bounded
in a fixed interval $[0,\smash{\delta_{\text{max}}^{(i)}}]$ and slightly
penalized in the cost function via the following penalty terms

$J_{\delta}^{(i)}=\lambda_{\delta}^{(i)}\delta^{(i)}\qquad i=4{},5{}$ (46)

A reasonable choice of the upper bound $\delta_{\text{max}}^{(i)}$ is between
1% and 10% of $\bar{s}^{(i)}$. Instead, the penalty weights $\lambda_{\delta}$
should be picked as small as possible (e.g., in the range
$\smash{10^{-6}}\div\smash{10^{-3}}$) to prevent converging to a suboptimal
solution.

The resulting second-order cone programming problem is

$\displaystyle\mathcal{P}_{B}:\;\;\min_{\begin{subarray}{c}\bm{x},\bm{u},s,\\\
\bm{q},\delta^{(4{})},\delta^{(5{})}\end{subarray}}$
$\displaystyle\quad-m(t_{f})+J_{q}+J_{\delta}^{(4{})}+J_{\delta}^{(5{})}$
(47a) s.t. (5), (11)–(15), (25), (33), (35)–(39), (41), (44), (45) (47b)

###  3.5 Discretization

Problem $\mathcal{P}_{B}$ in Eq. (47) is infinite-dimensional, as state and
controls are continuous-time functions. Therefore, to solve numerically the
OCP, we employ a $hp$ pseudospectral discretization method to convert it into
a finite set of variables and constraints. A $hp$ transcription combines the
advantages of $h$ and $p$ schemes, as it can introduce mesh nodes near
potential discontinuities and leverage exponential convergence rate in each
subinterval, resulting in a computationally efficient and accurate
discretization [49].

Specifically, in a $hp$ scheme, the time domain of each phase is split into
$h$ subintervals, or _segments_ , and the differential constraints are
enforced in each segment via local orthogonal collocation. The order $p$ of
the collocation can vary among the segments, but, for the sake of simplicity,
we adopt the same order over all the segments of a phase. The Radau
pseudospectral method (RPM) [50] is used as the collocation scheme since it is
one of the most accurate pseudospectral methods and avoids redundant control
variables at the segment boundaries [51].

The state and control are discretized over the grid, generating a finite set
of variables $(\bm{x}_{j}^{k},\bm{u}_{j}^{k})$. The superscript $k$ denotes
the $k$-th segment, while the subscript $j$ refers to the $j$-th node of the
segment. In each segment $[\tau_{k},\tau_{k+1}]$, the state and control
signals are approximated via Lagrange polynomial interpolation and the
derivative of the state approximation is imposed to be equal to the right-hand
side of the equations of motion (33) at the LGR collocation points

$\sum_{j=1}^{p+1}D_{ij}^{k}\bm{x}_{j}^{k}=\frac{\tau_{k+1}-\tau_{k}}{2}\bm{f}_{i}^{k}\qquad
i=1,\dots,p$ (48)

where
$\bm{f}_{i}^{k}=A_{i}^{k}\bm{x}_{i}^{k}+B_{i}^{k}\bm{u}_{i}^{k}+\Sigma_{i}^{k}s+\bm{c}_{i}^{k}+\bm{q}_{i}^{k}$.
The LGR differentiation matrix $D^{k}$ can be efficiently computed via
barycentric Lagrange interpolation [52].

Note that also the virtual control $\bm{q}$ and the linearization matrices
(29)–(32) are discretized over the mesh. The penalty function $P(\bm{q})$ in
Eq. (34) is chosen as the 1-norm of the vector containing all the discrete
virtual control variables $\bm{q}_{i}^{k}$.

###  3.6 Convergence Criteria

Problem $\mathcal{P}_{B}$ in Eq. (47) must be solved sequentially, i.e.,
recursively updating the reference values
$\\{\bar{\bm{x}},\bar{\bm{u}},\bar{s}\\}$ until convergence.

The update of the reference solution is based on _filtering_ , which consists
of updating the reference solution with a weighted sum of $K$ previous
solutions [16, 19]. Filtering provides further algorithmic robustness to the
successive convexification algorithm since it filters out diverging
intermediate iterations by averaging out the reference values over multiple
solutions. In practice, the reference solution for the $i$-th problem is
computed as

$\bar{x}^{(i)}=\sum_{k=1}^{K}\alpha_{k}x^{\max\\{0,(i-k)\\}}$ (49)

where $x^{(i)}$ denotes the solution to the $i$-th problem and $\alpha_{k}$
the corresponding weight. As in the authors’ previous works [16, 11, 19],
three previous solutions are used for the reference solution update, weighted
as $\alpha_{1}=6/11$, $\alpha_{2}=3/11$, and $\alpha_{3}=2/11$.

Eventually, the sequential algorithm terminates when all the following
criteria are met:

  1. (i)

the difference between the computed solution and the reference one converges
below an assigned tolerance $\epsilon_{\text{tol}}$

$\left\lVert\bm{x}-\bar{\bm{x}}\right\rVert_{\infty}<\epsilon_{\text{tol}}$
(50)

  2. (ii)

the computed solution adheres to the nonlinear dynamics within a tolerance
$\epsilon_{f}$

$\Big{\lVert}\sum_{j=1}^{p+1}D_{ij}^{k}\bm{x}_{j}^{k}-\frac{\tau_{k+1}-\tau_{k}}{2}\bm{f}_{i}^{k}\Big{\rVert}_{\infty}<\epsilon_{f}$
(51)

in each phase, for $i=1,\dots,p$, and $k=1,\dots,h$;

  3. (iii)

the virtual buffers of the computed solution are below the dynamics tolerance
$\epsilon_{f}$

$\left\lVert\bm{w}\right\rVert_{\infty}<\epsilon_{f}$ (52)

##  4 Model Predictive Control

Sections 2 and 3 outlined the optimization problem to solve to find the
optimal trajectory and control law in the unperturbed scenario. However,
because of modeling errors or external disturbances, the actual flight path of
the launch vehicle inevitably deviates from the nominal solution. Thus, a
closed-loop controller must be implemented to compensate for deviations and
ensure meeting all mission requirements even in presence of uncertainties. MPC
represents an appealing solution since it provides robustness to uncertainties
by repeatedly solving onboard problem $\mathcal{P}_{B}$ to update the
trajectory and control law based on the encountered conditions.

The considered MPC controller is implemented only during the third stage
flight, i.e., until stage burnout, which is the MPC stopping criterion.
Nevertheless, the _prediction_ horizon (i.e., the time domain of the OCP)
extends until the payload release into orbit and the splash-down of the third
stage. At the end of each control cycle, the OCP is updated with the
measurements coming from the navigation system. The update concerns the
initial condition of the OCP, i.e., Eq. (5), which is updated with the
measured state. The starting reference solution
$\\{\bar{\bm{x}},\bar{\bm{u}},\bar{s}\\}$ is replaced with the optimal
solution computed at the previous time, removing the portion of flight elapsed
between the two control steps.

Since the time-length of the prediction horizon reduces at each update due to
the receding-horizon implementation of the MPC algorithm, to save
computational resources, the size of the discretization grid of Phase 1 is
reduced linearly at every step from the value in Table 1 until a minimum value
of $p=$5$$.

Table 1: Discretization segments and order in each phase | Phase
---|---
| 1 | 2 | 3 | 4 | 5 | 6
$h$ | 1 | 1 | 1 | 1 | 1 | 5
$p$ | 19 | 9 | 19 | 19 | 19 | 20
Nodes | 20 | 10 | 20 | 20 | 20 | 101

To guarantee recursive feasibility, i.e., the formulation of a feasible OCP at
every time step, the constraint on maximum thermal flux (41) has to be
relaxed. Indeed, under nominal operating conditions, the virtual controls and
buffers would prevent artificial infeasibility. Yet, due to off-nominal
initial conditions or external disturbances in the dynamics, the heat flux may
be greater than expected at some points in the simulation and, since the
nominal heat flux constraint is active at some points in the trajectory, there
is little or no margin to keep it below the threshold value. Thus, the
infeasibility is solved by relaxing the heat flux constraint (41) as

$\dot{\bar{Q}}+\frac{\partial\dot{\bar{Q}}}{\partial\bm{r}}\cdot(\bm{r}-\bar{\bm{r}})+\frac{\partial\dot{\bar{Q}}}{\partial\bm{v}}\cdot(\bm{v}-\bar{\bm{v}})\leq\dot{Q}_{\text{max}}(1+\delta_{\dot{Q}})$
(53)

where $\delta_{\dot{Q}}$ is a non-negative optimization variable that is
penalized in the cost function via the term

$J_{\dot{Q}}=\lambda_{\dot{Q}}\delta_{\dot{Q}}$ (54)

In the authors’ experience, high values should not be assigned to the penalty
weight $\lambda_{\dot{Q}}$ since the optimization would exploit virtual
controls and may never meet the converge criterion (51). Rather, considering a
conservative nominal threshold $\dot{Q}_{\text{max}}$ and assigning a small
penalty to its violation appears as a much more effective strategy.

Finally, it is worthwhile mentioning that, in the present work, the time to
solve the optimization problem is supposed null, but the simulation should
introduce a delay between the measurement of the updated state and the
actuation of the optimal control law. This delay can introduce significant
deviations from the predicted trajectory if the computation time is long.
Nevertheless, thanks to the convexification strategy, the time required to
solve the OCP is assumed to be sufficiently brief. Preliminary numerical
results confirm the validity of this hypothesis.

##  5 Optimal Covariance Control

MPC indirectly provides robustness to the system thanks to the continuous
update of the trajectory and control that compensates for the measured
deviations from the nominal path. However, a sounder approach to designing a
robust controller consists in explicitly accounting for uncertainties in the
optimization process.

Generic in-flight disturbances can be modeled as additive white Gaussian noise
by considering the following (linearized) stochastic dynamics

$d\bm{x}=(A(t)\bm{x}+B(t)\bm{u}+\bm{c}(t))dt+G(t)d\bm{w}$ (55)

where $\bm{w}\in^{n_{w}}$ is a multidimensional Wiener process, i.e.,
$\mathbb{E}[d\bm{w}]=0$ and $\mathbb{E}[d\bm{w}d\bm{w}^{T}]=dt$, and
$G\in^{n_{x}\times n_{w}}$ is a matrix that determines how the external
disturbances affect the system. If the initial state distribution is Gaussian,
then the solution to Eq. (55) is a Gaussian process, that is,
$\bm{x}(t)\sim\mathcal{N}(\bm{\mu}(t),P(t))\;\forall t\geq t_{0}$, where
$\bm{\mu}=\mathbb{E}[\bm{x}]$ is the mean state and
$P=\mathbb{E}[(\bm{x}-\bm{\mu})(\bm{x}-\bm{\mu})^{T}]$ is the covariance
matrix.

As for the upper stage control problem, the objective of covariance control
is, starting from an initial state covariance, that is,

$P(t_{0})=P_{0}$ (56)

to minimize the dispersion of the splash-down location of the spent third
stage, i.e., the covariance at splash-down $P(t_{R})$.

###  5.1 Control Parameterization

Let us consider a zero-order hold (ZOH) control law, that is,
$\bm{u}(t)=\bm{u}(t_{j})\;\forall t\in[t_{j},t_{j+1}]$. Under this assumption,
we can discretize the state over a grid of $N$ nodes,

$t_{0}<t_{1}<\cdots<t_{N-1}$ (57)

The discrete-time stochastic system dynamics are

$\bm{x}_{j+1}=A_{j}\bm{x}_{j}+B_{j}\bm{u}_{j}+\bm{c}_{j}+G_{j}\bm{w}_{j}$ (58)

for $j=0,\dots,N-2$. The terms $A_{j}$, $B_{j}$, $\bm{c}_{j}$ are

$\displaystyle A_{j}$ $\displaystyle=\Phi(t_{j},t_{j+1})$ (59) $\displaystyle
B_{j}$ $\displaystyle=\int_{t_{j}}^{t_{j+1}}\Phi(s,t_{j+1})B(s)ds$ (60)
$\displaystyle\bm{c}_{j}$
$\displaystyle=\int_{t_{j}}^{t_{j+1}}\Phi(s,t_{j+1})\bm{c}(s)ds$ (61)

where $\Phi(t_{j},t_{j+1})$ denotes the state transition matrix from $t_{j}$
to $t_{j+1}$, and $G_{j}$ is a matrix such that $G_{j}\bm{w}_{j}$ is an
$n_{x}$-dimensional random Gaussian vector with covariance

$Q_{j}=\int_{t_{j}}^{t_{j+1}}\Phi(s,t_{j+1})G(s)G^{T}(s)\Phi^{T}(s,t_{j+1})ds$
(62)

To actively control the covariance of a linear system, a feedback control law
must be implemented. Specifically, since the SRM is a non-throttleable engine,
a multiplicative feedback law is considered, that is

$\bm{u}=\left(I+\left[{\bm{\alpha}}_{\times}\right]\right)\bm{\bar{\nu}}$ (63)

where the skew matrix $\left[{\bm{\alpha}}_{\times}\right]$ is

$\left[{\bm{\alpha}}_{\times}\right]=\begin{bmatrix}0&-\alpha_{3}&\alpha_{2}\\\
\alpha_{3}&0&-\alpha_{1}\\\ -\alpha_{2}&\alpha_{1}&0\end{bmatrix}$ (64)

and the vector $\bm{\alpha}$ is a linear state feedback

$\bm{\alpha}=K(\bm{x}-\bm{\bar{\mu}})$ (65)

The vectors $\bm{\bar{\mu}}$ and $\bm{\bar{\nu}}$ denote the nominal state and
control, which are the solution of the unperturbed problem $\mathcal{P}_{B}$.
Note that, if the $k$-th component of the $\bm{\alpha}$ vector, denoted by
$\alpha^{(k)}$, satisfies $\lvert\alpha^{(k)}\rvert\ll 1$ for $k=1,2,3$, then
the matrix $\left[I+\left[{\bm{\alpha}}_{\times}\right]\right]$ approximates a
rotation matrix and
$\left\lVert\bm{u}\right\rVert\approx\left\lVert\bm{\bar{\nu}}\right\rVert$.
Thus, the feedback action only rotates the nominal control $\bm{\bar{\nu}}$
without altering its magnitude.

###  5.2 Covariance Matrix Propagation

After replacing Eq. (63) in Eq. (58), the time-evolution of the state
covariance matrix $P$ is described by the following equation

$\displaystyle P_{j+1}$
$\displaystyle=(A_{j}-B_{j}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}K_{j})P_{j}(A_{j}-B_{j}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}K_{j})^{T}$
$\displaystyle+Q_{j}$ (66)

where the matrix $\left[{\bm{\bar{\nu}}}_{\times}\right]$ is

$\left[{\bm{\bar{\nu}}}_{\times}\right]=\begin{bmatrix}0&-\bar{\nu}_{z}&\bar{\nu}_{y}\\\
\bar{\nu}_{z}&0&-\bar{\nu}_{x}\\\ -\bar{\nu}_{y}&\bar{\nu}_{x}&0\end{bmatrix}$
(67)

Since both $K$ and $P$ are optimization variables of the covariance control
problem, Eq. (66) is a nonlinear constraint. Thus, we introduce a new control
variable

$Y_{j}=K_{j}P_{j}$ (68)

that replaces $K_{j}$ in Eq. (66) to obtain

$\displaystyle P_{j+1}$
$\displaystyle=A_{j}P_{j}A_{j}^{T}+A_{j}Y_{j}^{T}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}B_{j}^{T}-B_{j}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}Y_{j}A_{j}^{T}$
$\displaystyle+B_{j}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}Y_{j}P_{j}^{-1}Y_{j}^{T}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}^{T}B_{j}^{T}+Q_{j}$
(69)

Equation (69) is still nonlinear in $P_{j}$ and $Y_{j}$. Nevertheless, by
replacing the equality sign with the inequality sign, the following expression

$\displaystyle P_{j+1}$ $\displaystyle\succeq
A_{j}P_{j}A_{j}^{T}+A_{j}Y_{j}^{T}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}B_{j}^{T}-B_{j}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}Y_{j}A_{j}^{T}$
$\displaystyle+B_{j}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}Y_{j}P_{j}^{-1}Y_{j}^{T}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}^{T}B_{j}^{T}+Q_{j}$
(70)

can be written as a semidefinite cone constraint using Schur complement’s
lemma [53]

$\begin{bmatrix}P_{j}&Y_{j}^{T}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}^{T}B_{j}^{T}\\\
B_{j}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}Y_{j}&\Pi_{j}\end{bmatrix}\succeq
0$ (71)

where

$\displaystyle\Pi_{j}$
$\displaystyle=P_{j+1}-(A_{j}P_{j}A_{j}^{T}+A_{j}Y_{j}^{T}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}B_{j}^{T}$
$\displaystyle-
B_{j}\left[{\bm{\bar{\nu}}}_{\times}\right]_{j}Y_{j}A_{j}^{T}+Q_{j})$ (72)

The inequality sign in Eqs. (70) and (71) implies that the difference between
the left-hand side and the right-hand side is a positive semidefinite matrix.

Replacing Eq. (69) with Eq. (70) is a relaxation, as the covariance matrices
propagated with the inequality sign may have greater eigenvalues (hence,
greater state uncertainties) than the ones that would have been obtained
through forward propagation with the equality sign. Nevertheless, since the
objective of the optimization is to minimize the final covariance, we expect
that the optimal solution satisfies Eq. (70) with the equality sign (i.e., the
relaxation is lossless) like in the authors’ previous work [42]. Yet, this
statement should be verified by numerical experiments or by applying
Pontryagin’s maximum principle and inspecting the first-order necessary
conditions.

No control can be actuated on the system during the return phase. Thus, the
covariance matrix time-evolution reduces to

$P_{R,j+1}=A_{R,j}P_{R,j}A_{R,j}^{T}$ (73)

where the subscript $R$ indicates that these terms are associated with the
return (Phase 6). The return covariance is propagated over an evenly spaced
grid of $N_{R}$ points. Note that no random disturbance is added to the return
dynamics. Also, we assume that the covariance matrix is continuous across the
boundary between Phases 1 and 6, that is,

$P_{R,0}=P_{N-1}$ (74)

###  5.3 Control Magnitude Constraint

To ensure that
$\left\lVert\bm{u}\right\rVert\approx\left\lVert\bm{\bar{\nu}}\right\rVert$,
we must impose $\lvert\alpha^{(k)}\rvert\ll 1$ for $k=1,2,3$. In practice, one
can set a maximum angle $\alpha_{\text{max}}$ that cannot be exceeded, that
is,

$\left\lVert\bm{\alpha}_{j}\right\rVert_{\infty}\leq\alpha_{\text{max}}$ (75)

However, Eq. (75) cannot be enforced as a hard constraint since $\bm{\alpha}$
is a function of the state, which is normally distributed. Indeed,
$\bm{\alpha}$ is a normally distributed vector itself,
$\bm{\alpha}_{j}\sim\mathcal{N}(\bm{0},K_{j}P_{j}K_{j}^{T})$. Therefore, Eq.
(75) must be posed as a chance constraint, that is

$\text{Pr}\\{\left\lVert\bm{\alpha}_{j}\right\rVert_{\infty}\leq\alpha_{\text{max}}\\}\geq
p$ (76)

Equation (76) is numerically intractable, but it can be replaced with a
tractable expression if we recall that
$\left\lVert\bm{\alpha}_{j}\right\rVert_{\infty}\leq\left\lVert\bm{\alpha}_{j}\right\rVert_{2}$
and consider the following sufficient condition of the $L_{2}$ norm of a
normally distributed vector [41]

$\gamma^{2}(p)\rho(K_{j}P_{j}K_{j}^{T})\leq\alpha_{\text{max}}^{2}\Rightarrow\text{Pr}\\{\left\lVert\bm{\alpha}_{j}\right\rVert_{2}\leq\alpha_{\text{max}}\\}\geq
p$ (77)

where $\rho(M)$ denotes the spectral radius of a matrix $M$, $\gamma(p)$ is a
coefficient that depends on the desired probability level $p$ as

$\gamma(p)=\begin{cases}\sqrt{2\ln\left(1/(1-p)\right)}&\text{if
}n_{\alpha}\leq 2\\\
\sqrt{2\ln\left(1/(1-p)\right)}+\sqrt{n_{\alpha}}&\text{if
}n_{\alpha}>2\end{cases}$ (78)

and $n_{\alpha}$ is the number of components of $\bm{\alpha}$.

The spectral radius of the positive semidefinite matrix
$K_{j}P_{j}K_{j}^{T}=Y_{j}P_{j}^{-1}Y_{j}^{T}$ corresponds to its maximum
eigenvalue, which can be characterized in epigraph form
$\lambda_{\text{max}}(Y_{j}P_{j}^{-1}Y_{j}^{T})\leq\tau_{j}$, where $\tau_{j}$
is an auxiliary variable that satisfies

$\tau_{j}I-Y_{j}P_{j}^{-1}Y_{j}^{T}\succeq 0$ (79)

Equation (79) is particularly convenient since it can be written as a
semidefinite cone constraint by leveraging Schur complement’s lemma, that is,

$\begin{bmatrix}\tau_{j}I&Y_{j}\\\ Y_{j}^{T}&P_{j}\end{bmatrix}\succeq 0$ (80)

Replacing the spectral radius with $\tau_{j}$, the sufficient condition in Eq.
(77) can be written as a linear constraint

$\tau_{j}\leq\alpha_{\text{max}}^{2}/\gamma^{2}(p)$ (81)

###  5.4 Objective

The goal of the controller is to minimize the dispersion of the splash-down
point. Since only the position of the splash-down point is relevant, while
velocity and mass at impact are not a concern, the objective function is the
maximum eigenvalue of the position elements of the covariance matrix at
$t_{R}$, that is, the first 3 rows and columns of $P_{R,N_{R}-1}$. The maximum
eigenvalue is characterized in epigraph form
$\lambda_{\text{max}}(P_{R,N_{R}-1}^{(1:3\times 1:3)})\leq\tau_{R}$ as

$\tau_{R}I-P_{R,N_{R}-1}^{(1:3\times 1:3)}\succeq 0$ (82)

To sum up, given a nominal trajectory $\bm{\bar{\mu}}$ and control
$\bm{\bar{\nu}}$, the optimal covariance control problem to solve is

$\displaystyle\mathcal{P}_{C}:\;\;\min_{\begin{subarray}{c}P_{j},P_{R,j},\\\
Y_{j},\tau_{j},\tau_{R}\end{subarray}}$ $\displaystyle\quad\tau_{R}$ (83a)
s.t. (56), (71), (73), (74), (80)–(82) (83b)

$\mathcal{P}_{C}$ is a semidefinite programming problem that can be solved in
polynomial time with interior-point algorithms.

##  6 Numerical Results

The solution algorithms were implemented in C++. Gurobi’s [54] second-order
cone programming solver was used to solve $\mathcal{P}_{B}$, while MOSEK’s
[55] semidefinite programming solver was used for $\mathcal{P}_{C}$. All the
computations were carried out on a computer equipped with Intel® Core™
i7-9700K CPU @ $3.60\text{\,}\mathrm{GHz}\text{/}$.

###  6.1 Case Study

The investigated scenario is the same as in Ref. [19], that is, an ascent
trajectory from an equatorial launch base to a
$700\text{\,}\mathrm{km}\text{/}$ circular orbit with
$i_{\text{des}}=$90\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$$ and a
splash-down of the third stage targeted at a latitude of
$\varphi_{R,\text{des}}=$60\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$$.
Also, the same data are used to model the launch vehicle stages, which are
summarized in Table 2. In addition to these data, the dry mass of the third
stage, which is the spent stage mass, is $m_{\text{dry,
3}}=$1326.5\text{\,}\mathrm{kg}\text{/}$$. The burn time is denoted by
$t_{b}$, and the vacuum thrust and mass flow rate profiles are modeled as
linear functions of time. As for the aerodynamics, the $C_{D}=$0.38$$ and
$S_{\text{ref}}=$9.08\text{\,}{\mathrm{m}}^{2}\text{/}$$ in all phases. The
U.S. Standard Atmosphere 1976 model is used to evaluate the air density and
pressure as functions of the altitude [56].

Table 2: Engine data Quantity | Stage 3 | Stage 4 | Unit
---|---|---|---
$t_{b}$ | $104.6$ | Free | $\mathrm{s}\text{/}$
$T_{\text{vac, 0}}$ | $295.97$ | $2.45$ | $\mathrm{kN}\text{/}$
$T_{\text{vac, f}}$ | $221.60$ | $2.45$ | $\mathrm{kN}\text{/}$
$\dot{m}_{e,0}$ | $103.27$ | $0.79$ | $\mathrm{kg}\text{/}\mathrm{s}$
$\dot{m}_{e,f}$ | $77.32$ | $0.79$ | $\mathrm{kg}\text{/}\mathrm{s}$
$A_{e}$ | $1.18$ | $0.07$ | ${\mathrm{m}}^{2}\text{/}$

Phase 1 lasts for the entire burn duration of Stage 3, i.e., $\Delta
t^{(1{})}=t_{b,3}$. The duration of the coasting arc after the third stage
separation (Phase 2) is fixed to $\Delta
t^{(2{})}=$15.4\text{\,}\mathrm{s}\text{/}$$. The time-lengths of all the
following phases are left to be optimized. The threshold on the bearable heat
flux is set to
$\dot{Q}_{\text{max}}=$900\text{\,}\mathrm{W}\text{/}{\mathrm{m}}^{2}$$. This
value of $\dot{Q}_{\text{max}}$ is lower than the value
($1135\text{\,}\mathrm{W}\text{/}{\mathrm{m}}^{2}$) reported in the official
VEGA manual [46], to allow for the relaxation of the heat flux constraint as
in Eq. (53) with a relatively small value of the penalty coefficient
($\lambda_{\dot{Q}}=${10}^{-2}$$) in the MPC simulations.

###  6.2 Nominal Trajectory

The nominal trajectory was computed by sequentially solving the convex
optimization problem $\mathcal{P}_{B}$ outlined in Eq. (47). The same
parameters as in Ref. [19] were used for the successive convexification
algorithm. Specifically, the penalty weights on the time-lengths,
$\smash{\lambda_{\delta}^{(4{})}}$ and $\smash{\lambda_{\delta}^{(5{})}}$,
were set to ${10}^{-4}$, while, to highly penalize solutions that actively
exploit virtual variables, we set $\lambda_{q}=\lambda_{w}=${10}^{4}$$. The
tolerances on the convergence criteria, Eqs. (50)–(52), were prescribed as
$\epsilon_{\text{tol}}=${10}^{-4}$$ and $\epsilon_{f}=${10}^{-6}$$

First, the entire ascent, from Stage 1 ignition until payload release into
orbit was optimized neglecting the return of the spent third stage. Thanks to
the robustness of the convex optimization algorithm, even a rough initial
guess leads to the optimal solution in a short time ($25\text{\,}\mathrm{s}$
on average)111A deeper analysis of the performance and efficiency of the
convex optimization of the ascent trajectory is available in Ref. [19]. Then,
the return trajectory was simulated and the splash-down point was gradually
moved from the unconstrained latitude to $\varphi_{R,\text{des}}$. The
resulting solution features a nominal payload mass of
$1400.1\text{\,}\mathrm{kg}\text{/}$ and the corresponding trajectory is
illustrated in Fig. 2.

Figure 2: Visualization of the nominal trajectory, with the third stage return
flight colored in red.

###  6.3 Stochastic Control

In the stochastic optimal control problem, two sources of uncertainty are
considered: dispersion of the initial conditions and in-flight random
disturbances. The initial state is normally distributed and features a
diagonal covariance matrix. Specifically, at the initial time, the position
variance is ${10}^{4}\text{\,}{\mathrm{m}}^{2}\text{/}$, the velocity variance
is ${10}^{2}\text{\,}{\mathrm{m}}^{2}\text{/}{\mathrm{s}}^{2}$, and the mass
variance is $1\text{\,}{\mathrm{kg}}^{2}\text{/}$. On the other hand, the
external disturbance is a random perturbing acceleration (hence, the equations
of position and mass are unperturbed). Thus, the dimension of the Wiener
process is $n_{w}=3$ and the $G$ matrix structure is

$G=\begin{bmatrix}0_{3\times 3}&g_{v}I_{3}&0_{1\times 3}^{T}\end{bmatrix}^{T}$
(84)

where $g_{v}$ determines the intensity of the perturbation. In particular,
three levels of in-flight disturbance intensity were considered as reported in
Table 3 and called Low (L), Medium (M), and High (H).

Table 3: In-flight disturbance intensities Parameter | Case L | Case M | Case H | Unit
---|---|---|---|---
$g_{v}$ | $0.1$ | $0.5$ | $1$ | $\mathrm{m}\mathrm{/}\mathrm{s}^{3/2}$

The stochastic control problem was discretized over an evenly spaced grid of
$N=161$ nodes in Phase 1 and $N_{R}=1001$ nodes in Phase 6. The maximum
control angle was set to
$\alpha_{\text{max}}=$5\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$$ and the
corresponding chance constraint in Eq. (76) was imposed with a probability
level $p=0.95$. Problem $\mathcal{P}_{C}$ was solved in approximately
$40\text{\,}\mathrm{s}\text{/}$.

###  6.4 Monte Carlo Campaigns

A Monte Carlo analysis of the combined effect of off-nominal initial
conditions and in-flight disturbance was carried out for both the MPC
algorithm and the stochastic controller resulting from problem
$\mathcal{P}_{C}$. For each disturbance intensity level, 400 independent
simulations were carried out. Each simulation consisted of a numerical
integration of the nonlinear stochastic differential equations, i.e., the
original Eqs. (2)–(4) perturbed by the Wiener process multiplied by the $G$
matrix in Eq. (84), using the Euler-Maruyama method.

The update frequency of the MPC was set to $1\text{\,}\mathrm{Hz}\text{/}$,
meaning that the OCP is solved every $T=$1\text{\,}\mathrm{s}\text{/}$$. This
is the expected highest achievable rate on a standard flight hardware. Note
that the solution of the OCP requires approximately
$1\text{\,}\mathrm{s}\text{/}$ on the employed working station, but
significant speed-ups may be attained through custom code optimization or
dedicated hardware.

(a) Open-loop

(b) Closed-loop (MPC)

(c) Closed-loop (Stochastic control)

Figure 3: Altitude profiles for Case H.

Figure 3 reports the altitude profile of the third stage from ignition until
splash-down versus the range, i.e., the distance from the launch site measured
in degrees, in case H. The black line denotes the nominal flight path, i.e.,
the solution of Problem $\mathcal{P}_{B}$, while the gray lines are the Monte
Carlo simulations. In the Monte Carlo simulations, the Gaussian noise is added
to the equations of motion only during the propelled flight of the third
stage; after burnout, which is marked by a black dot in the plots, no
perturbation is included and the ballistic return is simulated via fourth-
order Runge-Kutta integration of the ordinary differential equations. It is
apparent from Fig. 3 that the dispersion of the splash-down point is huge if
no feedback controller is implemented. Instead, if an MPC algorithm or the
stochastic closed-loop controller is implemented, the splash-down dispersion
reduces to acceptable values.

Figure 4: Splash-down footprints for Case H.

Figure 4 shows the splash-down footprints attained via MPC and stochastic
control. The footprints are defined as the 95%-confidence ellipses of the
splash-down point. Stochastic control achieves smaller dispersion than MPC, as
the latitude error standard deviation is almost double in the latter case.
Instead, the longitude dispersion is comparable for both control strategies.

Table 4: Results of the Monte Carlo campaigns Method |  Case | Payload, kg |  $\varphi_{R}$, deg | Footprint axes, km
---|---|---|---|---
Min | Mean | Max | $\sigma$ | Min | Mean | Max | $\sigma$ | Major | Minor
MPC | L | $1397.7$ | $1400.0$ | $1402.2$ | $0.71$ | $59.93$ | $60.00$ | $60.10$ | $0.03$ | $15.09$ | $2.24$
M | $1386.4$ | $1399.9$ | $1408.5$ | $3.48$ | $59.49$ | $60.00$ | $60.44$ | $0.14$ | $76.96$ | $10.64$
H | $1375.6$ | $1399.2$ | $1418.5$ | $6.81$ | $59.34$ | $59.99$ | $60.91$ | $0.29$ | $157.57$ | $21.13$
Stochastic Control | L | $1377.8$ | $1399.8$ | $1418.2$ | $7.14$ | $59.91$ | $60.00$ | $60.07$ | $0.03$ | $15.27$ | $6.62$
M | $1375.1$ | $1400.2$ | $1423.0$ | $7.46$ | $59.81$ | $60.00$ | $60.20$ | $0.07$ | $39.51$ | $11.56$
H | $1370.4$ | $1400.3$ | $1427.2$ | $9.98$ | $59.59$ | $59.98$ | $60.47$ | $0.13$ | $71.63$ | $27.73$

Detailed results of the Monte Carlo campaigns are reported in Table 4. As the
intensity of the in-flight disturbance increases, larger dispersions on the
splash-down point are observed for both control strategies.

##  7 Conclusions

This paper presented a guidance and control strategy for the upper stage of a
launch vehicle. Specifically, a mindful convexification process was outlined
to find an optimal nominal trajectory with low computational complexity. Then,
two control synthesis methods were analyzed. The first one consists in
embedding the convex optimization algorithm in an MPC framework, which, thanks
to the onboard recomputation of the flight path and control signal, is
inherently robust to model mismatches and random inflight disturbances. The
second approach explicitly accounts for uncertainties in the optimization
process by considering stochastic system dynamics, perturbed by additive white
Gaussian noise, to find a multiplicative feedback controller that guarantees
robustness to the considered characterization of the disturbances. Thanks to
an original convexification process, the associated stochastic optimal control
problem could be solved in polynomial time.

The performance of the MPC algorithm and the stochastic controller was
compared through a Monte Carlo analysis, with a particular focus on the
ability to minimize the splash-down footprint of the spent stage. While MPC
proved to be quite robust thanks to the high update frequencies achievable
because of the short solution time of the sequential convex optimization
algorithm, the stochastic controller turned out to better reject disturbances
and attain a smaller dispersion of the splash-down point. Also, MPC
necessitates quite performing onboard computation hardware, while the
stochastic controller can be designed preflight and does not require any
onboard computation.

Further tests are necessary to prove the suitability of the proposed
stochastic control synthesis method for onboard implementation, but the
preliminary results are very encouraging. Indeed, a computationally efficient
and systematic design of a nominal trajectory and a robust controller is a
crucial technology for launch vehicle guidance and control algorithms.

## Acknowledgments

This work was supported by the Agreement n. 2019-4–HH.0 CUP F86C17000080005
“Technical Assistance on Launch Vehicles and Propulsion” between the Italian
Space Agency and the Department of Mechanical and Aerospace Engineering of
Sapienza University of Rome.

## References

  * [1] R. Bellman, “Dynamic programming,” Science, vol. 153, no. 3731, pp. 34–37, 1966.
  * [2] D. Jacobson and D. Mayne, Differential Dynamic Programming.  Modern analytic and computational methods in science and mathematics, American Elsevier Publishing Company, 1970.
  * [3] N. Ozaki, S. Campagnola, and R. Funase, “Tube stochastic optimal control for nonlinear constrained trajectory optimization problems,” Journal of Guidance, Control, and Dynamics, vol. 43, no. 4, pp. 645–655, 2020.
  * [4] D. Izzo, M. Märtens, and B. Pan, “A survey on artificial intelligence trends in spacecraft guidance dynamics and control,” Astrodynamics, pp. 1–13, 2018.
  * [5] L. Cheng, Z. Wang, Y. Song, and F. Jiang, “Real-time optimal control for irregular asteroid landings using deep neural networks,” Acta Astronautica, vol. 170, pp. 66 – 79, 2020.
  * [6] L. Federici, B. Benedikter, and A. Zavoli, “Deep learning techniques for autonomous spacecraft guidance during proximity operations,” Journal of Spacecraft and Rockets, pp. 1–12, 2021.
  * [7] A. Zavoli and L. Federici, “Reinforcement learning for robust trajectory design of interplanetary missions,” Journal of Guidance, Control, and Dynamics, vol. 44, no. 8, pp. 1440–1453, 2021.
  * [8] L. Federici, A. Scorsoglio, A. Zavoli, and R. Furfaro, “Autonomous guidance for cislunar orbit transfers via reinforcement learning,” in AAS/AIAA Astrodynamics Specialist Conference, 2021.
  * [9] L. Federici, A. Scorsoglio, A. Zavoli, and R. Furfaro, “Reinforcement meta-learning for adaptive spacecraft guidance during missions to near-earth objects,” in 72nd International Astronautical Congress (IAC), 2021.
  * [10] U. Eren, A. Prach, B. B. Koçer, S. V. Raković, E. Kayacan, and B. Açıkmeşe, “Model predictive control in aerospace systems: Current state and opportunities,” Journal of Guidance, Control, and Dynamics, vol. 40, no. 7, pp. 1541–1566, 2017.
  * [11] B. Benedikter, A. Zavoli, G. Colasurdo, S. Pizzurro, and E. Cavallini, “Autonomous upper stage guidance using convex optimization and model predictive control,” in AIAA ASCEND, 2020.
  * [12] B. Benedikter, A. Zavoli, G. Colasurdo, S. Pizzurro, and E. Cavallini, “Autonomous upper stage guidance with robust splash-down constraint,” in 72nd International Astronautical Congress (IAC), 2021.
  * [13] P. Iannelli, F. Angeletti, and P. Gasbarri, “An attitude/spin control architecture for a spacecraft equipped with a flexible rotating payload based on model predictive control,” in 72nd International Astronautical Congress (IAC), 2021.
  * [14] X. Liu, P. Lu, and B. Pan, “Survey of convex optimization for aerospace applications,” Astrodynamics, vol. 1, no. 1, pp. 23–40, 2017.
  * [15] B. Açıkmeşe and L. Blackmore, “Lossless convexification of a class of optimal control problems with non-convex control constraints,” Automatica, vol. 47, no. 2, pp. 341–347, 2011.
  * [16] B. Benedikter, A. Zavoli, and G. Colasurdo, “A convex optimization approach for finite-thrust time-constrained cooperative rendezvous,” Advances in the Astronautical Sciences, vol. 171, pp. 1483–1498, 2019.
  * [17] B. Benedikter, A. Zavoli, and G. Colasurdo, “A convex approach to rocket ascent trajectory optimization,” in 8th European Conference for Aeronautics and Space Sciences (EUCASS), 2019.
  * [18] B. Benedikter, A. Zavoli, G. Colasurdo, S. Pizzurro, and E. Cavallini, “Convex approach to three-dimensional launch vehicle ascent trajectory optimization,” Journal of Guidance, Control, and Dynamics, vol. 44, no. 6, pp. 1116–1131, 2021.
  * [19] B. Benedikter, A. Zavoli, G. Colasurdo, S. Pizzurro, and E. Cavallini, “Convex optimization of launch vehicle ascent trajectory with heat-flux and splash-down constraints,” Journal of Spacecraft and Rockets, vol. 59, no. 3, pp. 900–915, 2022.
  * [20] M. Szmuk, B. Acikmese, and A. W. Berning, “Successive convexification for fuel-optimal powered landing with aerodynamic drag and non-convex constraints,” in AIAA Guidance, Navigation, and Control Conference, 2016\.
  * [21] M. Sagliano, “Pseudospectral convex optimization for powered descent and landing,” Journal of Guidance, Control, and Dynamics, vol. 41, no. 2, pp. 320–334, 2018.
  * [22] Z. Wang and M. J. Grant, “Optimization of minimum-time low-thrust transfers using convex programming,” Journal of Spacecraft and Rockets, vol. 55, no. 3, pp. 586–598, 2018.
  * [23] Z. Wang and M. J. Grant, “Minimum-fuel low-thrust transfers for spacecraft: A convex approach,” IEEE Transactions on Aerospace and Electronic Systems, vol. 54, no. 5, pp. 2274–2290, 2018.
  * [24] Z. Wang and M. J. Grant, “Autonomous entry guidance for hypersonic vehicles by convex optimization,” Journal of Spacecraft and Rockets, vol. 55, no. 4, pp. 993–1006, 2018.
  * [25] M. Sagliano and E. Mooij, “Optimal drag-energy entry guidance via pseudospectral convex optimization,” in AIAA Guidance, Navigation, and Control Conference, 2018.
  * [26] A. Bemporad and M. Morari, “Robust model predictive control: A survey,” in Robustness in identification and control, pp. 207–226, 1999.
  * [27] L. Chisci, J. Rossiter, and G. Zappa, “Systems with persistent disturbances: predictive control with restricted constraints,” Automatica, vol. 37, no. 7, pp. 1019–1028, 2001.
  * [28] W. Langson, I. Chryssochoos, S. Raković, and D. Mayne, “Robust model predictive control using tubes,” Automatica, vol. 40, no. 1, pp. 125–133, 2004.
  * [29] A. Mesbah, “Stochastic model predictive control: An overview and perspectives for future research,” IEEE Control Systems Magazine, vol. 36, no. 6, pp. 30–44, 2016.
  * [30] L. Blackmore, M. Ono, and B. C. Williams, “Chance-constrained optimal path planning with obstacles,” IEEE Transactions on Robotics, vol. 27, no. 6, pp. 1080–1094, 2011.
  * [31] T. Lew, F. Lyck, and G. Müller, “Chance-constrained optimal altitude control of a rocket,” in European Conference for Aeronautics and Space Sciences (EUCASS), 2019.
  * [32] K. Oguri, M. Ono, and J. W. McMahon, “Convex optimization over sequential linear feedback policies with continuous-time chance constraints,” in IEEE 58th Conference on Decision and Control (CDC), pp. 6325–6331, 2019.
  * [33] K. Oguri and J. W. McMahon, “Risk-aware trajectory design with impulsive maneuvers: Convex optimization approach,” Advances in the Astronautical Sciences, vol. 171, pp. 1985–2004, 2020.
  * [34] K. Oguri and J. W. McMahon, “Risk-aware trajectory design with continuous thrust: Primer vector theory approach,” Advances in the Astronautical Sciences, vol. 171, pp. 2049–2067, 2020.
  * [35] T. Lew, R. Bonalli, and M. Pavone, “Chance-constrained sequential convex programming for robust trajectory optimization,” in European Control Conference (ECC), pp. 1871–1878, 2020.
  * [36] A. Hotz and R. E. Skelton, “Covariance control theory,” International Journal of Control, vol. 46, no. 1, pp. 13–32, 1987.
  * [37] Y. Chen, T. T. Georgiou, and M. Pavon, “Optimal steering of a linear stochastic system to a final probability distribution, part i,” IEEE Transactions on Automatic Control, vol. 61, no. 5, pp. 1158–1169, 2016.
  * [38] Y. Chen, T. T. Georgiou, and M. Pavon, “Optimal steering of a linear stochastic system to a final probability distribution, part ii,” IEEE Transactions on Automatic Control, vol. 61, no. 5, pp. 1170–1180, 2016.
  * [39] E. Bakolas, “Optimal covariance control for stochastic linear systems subject to integral quadratic state constraints,” in American Control Conference (ACC), pp. 7231–7236, 2016.
  * [40] E. Bakolas, “Optimal covariance control for discrete-time stochastic linear systems subject to constraints,” in IEEE 55th Conference on Decision and Control (CDC), pp. 1153–1158, 2016.
  * [41] J. Ridderhof, J. Pilipovsky, and P. Tsiotras, “Chance-constrained covariance control for low-thrust minimum-fuel trajectory optimization,” in AAS/AIAA Astrodynamics Specialist Conference, 2020.
  * [42] B. Benedikter, A. Zavoli, Z. Wang, S. Pizzurro, and E. Cavallini, “Covariance control for stochastic low-thrust trajectory optimization,” in AIAA SciTech Forum, 2022.
  * [43] J. Ridderhof and P. Tsiotras, “Minimum-fuel closed-loop powered descent guidance with stochastically derived throttle margins,” Journal of Guidance, Control, and Dynamics, vol. 44, no. 3, pp. 537–547, 2021.
  * [44] B. Benedikter, A. Zavoli, Z. Wang, S. Pizzurro, and E. Cavallini, “Convex approach to stochastic control for autonomous rocket pinpoint landing,” in AAS/AIAA Astrodynamics Specialist Conference, 2022.
  * [45] J. Ridderhof, P. Tsiotras, and B. J. Johnson, “Stochastic entry guidance,” Journal of Guidance, Control, and Dynamics, vol. 45, no. 2, pp. 320–334, 2022.
  * [46] Arianespace, Vega User’s Manual, 2014.
  * [47] X. Liu, Z. Shen, and P. Lu, “Entry trajectory optimization by second-order cone programming,” Journal of Guidance, Control, and Dynamics, vol. 39, no. 2, pp. 227–241, 2016.
  * [48] B. Acikmese and S. R. Ploen, “Convex programming approach to powered descent guidance for mars landing,” Journal of Guidance, Control, and Dynamics, vol. 30, no. 5, pp. 1353–1366, 2007.
  * [49] C. L. Darby, W. W. Hager, and A. V. Rao, “An hp-adaptive pseudospectral method for solving optimal control problems,” Optimal Control Applications and Methods, vol. 32, no. 4, pp. 476–502, 2011.
  * [50] D. Garg, Advances in Global Pseudospectral Methods for Optimal Control.  PhD thesis, University of Florida, Gainesville, FL, 2011.
  * [51] D. Garg, M. Patterson, W. W. Hager, A. V. Rao, D. A. Benson, and G. T. Huntington, “A unified framework for the numerical solution of optimal control problems using pseudospectral methods,” Automatica, vol. 46, no. 11, pp. 1843 – 1851, 2010.
  * [52] J.-P. Berrut and L. N. Trefethen, “Barycentric lagrange interpolation,” SIAM Review, vol. 46, no. 3, pp. 501–517, 2004.
  * [53] S. Boyd, S. P. Boyd, and L. Vandenberghe, Convex Optimization.  Cambridge University Press, 2004.
  * [54] Gurobi Optimization, LLC, Gurobi Optimizer Reference Manual, 2020.
  * [55] MOSEK ApS, MOSEK Fusion API for C++, 2021.
  * [56] NOAA, NASA, and USAF, U.S. Standard Atmosphere.  U.S. Government Printing Office, 1976.