$(\varepsilon,d)=(+1,0)$ and $r=2$, low dimensionality forces the extra
relation $q(\tilde{\phi}_{4}(\pi_{4}(\mathbf{p})))=0$, which causes that
$G_{2}\\!\setminus\\!\mathcal{P}_{2}^{5}\cong\mathbb{C}^{4}$.

## References

  * [1] S.J. Bloch. _Higher regulators, algebraic K-theory, and zeta functions of elliptic curves._ CRM Monograph Series, 11. American Mathematical Society, Providence, RI, 2000.
  * [2] A. Borel, N. Wallach. _Continuous cohomology, discrete subgroups, and representations of reductive groups_. Second edition. Mathematical Surveys and Monographs, 67. American Mathematical Society, Providence, RI, 2000.
  * [3] M. Bucher, M. Burger, A. Iozzi. _The bounded Borel class and 3-manifold groups_. Duke Math. J. 167, no. 17 (2018), 3129–3169.
  * [4] M. Burger, A. Iozzi, A. Wienhard. _Surface group representations with maximal Toledo invariant_. Ann. Math. 172 (2010), no. 1, 517–566.
  * [5] M. Burger, N. Monod. _Bounded cohomology of lattices in higher rank Lie groups._ J. Eur. Math. Soc. 1 (1999), no. 2, 199–235.
  * [6] M. Burger, N. Monod. _Continuous bounded cohomology and applications to rigidity theory_. Geom. Funct. Anal. 12 (2002), no. 2, 219–280.
  * [7] M. Burger, N. Monod. _On and around the bounded cohomology of $\operatorname{SL}_{2}$_. Rigidity in dynamics and geometry, 19–37, Springer, Berlin, 2002.
  * [8] C. Campagnolo, F. Fournier-Facio, N. Heuer, M. Moraschini. _Bounded Cohomology and Simplicial Volume_. London Mathematical Society Lecture Note Series, Cambridge University Press, Cambridge, 2022.
  * [9] J. L. Dupont. _Bounds for characteristic numbers of flat bundles_. Algebraic topology, Aarhus 1978 (Proc. Sympos., Univ. Aarhus, Aarhus, 1978), pp. 109-119, Lecture Notes in Math., 763, Springer, Berlin, 1979.
  * [10] C. De la Cruz Mengual, _On Bounded-Cohomological Stability for Classical Groups._ ETH Zurich (2019), Zurich.
  * [11] C. De la Cruz Mengual, T. Hartnick. _A Quillen stability criterion for bounded cohomology._ arXiv preprint (2023), arXiv:2307.12808.
  * [12] C. De la Cruz Mengual, T. Hartnick. _Stabilization of Bounded Cohomology for Classical Groups_. arXiv preprint (2022), arXiv:2201.03879.
  * [13] J. Farre. _Relations in bounded cohomology._ J. Topol. 13 (2020), no. 3, 1084–1118.
  * [14] R. Frigerio. _Bounded cohomology of discrete groups_. Mathematical Surveys and Monographs, 227. American Mathematical Society, Providence, RI, 2017.
  * [15] S. Galatius, A. Kupers, O. Randal-Williams. $\mathrm{E}_{2}$-cells and mapping class groups. Publ. Math., Inst. Hautes Étud. Sci. 130 (2019), 1–61.
  * [16] W. M. Goldman. _Complex hyperbolic geometry._ Oxford Mathematical Monographs. Oxford Science Publications. Clarendon Press. xx, 316 p.m Oxford, 1999.
  * [17] A.B. Goncharov. _Explicit construction of characteristic classes_. I. M. Gel’fand Seminar, 169–210, Adv. Soviet Math., 16 (1993), Part 1, Amer. Math. Soc., Providence, RI.
  * [18] W. Greub, S. Halperin, R. Vanstone. _Connections, curvature and cohomology, Volume III: Cohomology of principal bundles and homogeneous spaces_. Pure and Applied Mathematics, Vol. 47-III. Academic Press [Harcourt Brace Jovanovich, Publishers], New York–London, 1976.
  * [19] M. Gromov. _Volume and bounded cohomology_. Inst. Hautes Études Sci. Publ. Math. 56 (1982), 5–99.
  * [20] T. Hartnick, A. Ott. _Bounded cohomology via partial differential equations, I._ Geom. Topol. 19 (2015), no. 6, 3603–3643.
  * [21] T. Hartnick, A. Ott. _Surjectivity of the Comparison Map in Bounded Cohomology for Hermitian Lie Groups_. Int. Math. Res. Not. IMRN 2012 (2012), no. 9, 2068–2093.
  * [22] T. Kastenholz, R. J. Sroka. _Simplicial bounded cohomology and stability_. arXiv preprint (2023), arXiv:2309.05024.
  * [23] A. Korányi, H. M. Reimann. _The complex cross ratio on the Heisenberg group._ Enseign. Math., II. Sér. (1987), no. 33, 291–300.
  * [24] B. Kostant. _The principal three-dimensional subgroup and the Betti numbers of a complex simple Lie group._ Am. J. Math. 81 (1959), 973–1032.
  * [25] J-F. Lafont, B. Schmidt. _Simplicial volume of closed locally symmetric spaces of noncompact type._ Acta Math. 197 (2006), no. 1, 129–143.
  * [26] M. Mimura, H. Toda. _Topology of Lie groups. I, II_. Translated from the 1978 Japanese edition by the authors. Translations of Mathematical Monographs, 91. American Mathematical Society, Providence, RI, 1991.
  * [27] N. Monod. _An invitation to bounded cohomology_. International Congress of Mathematicians. Vol. II, 1183–1211, Eur. Math. Soc., Zürich, 2006.
  * [28] N. Monod. _Continuous bounded cohomology of locally compact groups_. Lecture Notes in Mathematics, 1758. Springer-Verlag, Berlin, 2001.
  * [29] N. Monod. _Equivariant measurable liftings_. Fundam. Math. 230 (2015), no. 2, 149–165.
  * [30] N. Monod. _Stabilization for $\operatorname{SL}_{n}$ in bounded cohomology_. Contemp. Math. 347 (2004), 191–202.
  * [31] A. L. Onishchik. _Topology of Transitive Transformation Groups._ Johann Ambrosius Barth Verlag GmbH, Leipzig (1994).
  * [32] D. I. Panyushev. _On the Dynkin index of a principal $\mathfrak{sl}_{2}$-subalgebra_. Adv. Math. 221 (2009), no. 4, 1115–1121.
  * [33] J. D. Stasheff. _Continuous cohomology of groups and classifying spaces_. Bull. Amer. Math. Soc. 84 (1978), no. 4, 513–530.
  * [34] T. Tao. _Exceptional isogenies between the classical groups_. Blog entry in "What’s new" (11 March, 2011) https://terrytao.wordpress.com/2011/03/11/exceptional-isogenies-between-the-classical-lie-groups/.
  * [35] D.E. Taylor. _The geometry of the classical groups._ Sigma Series in Pure Mathematics, 9. Heldermann Verlag, Berlin, 1992.
  * [36] B. L. van der Waerden. _Gruppen von Linearen Transformationen_. Springer-Verlag, Berlin, 1935.
  * [37] D. Zagier. _The dilogarithm function_. Frontiers in number theory, physics, and geometry. II, 3–65, Springer, Berlin, 2007.
  * [38] R. J. Zimmer. _Ergodic theory and semisimple groups_. Monographs in Mathematics, Vol. 81. Birkhäuser, Boston-Basel-Stuttgart, 1984.


21 D2.3.8 72 Dubai, United Arab Emirates 25-29 October 2021 2021Mr. Boris
Benedikter

Boris Benedikter${}^{\text{a}\ast}$, ${}^{\text{a}}$PhD Student, Department of
Mechanical and Aerospace Engineering, Sapienza University of Rome, Via
Eudossiana 18, 00184, Rome, Italy<EMAIL_ADDRESS>
Alessandro Zavoli${}^{\text{b}}$, ${}^{\text{b}}$Research Assistant,
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>
# Autonomous Upper Stage Guidance with Robust Splash-Down Constraint

###### Abstract

This paper presents a novel algorithm, based on model predictive control
(MPC), for the optimal guidance of a launch vehicle upper stage. The proposed
strategy not only maximizes the performance of the vehicle and its robustness
to external disturbances, but also robustly enforces the splash-down
constraint. Indeed, uncertainty on the engine performance, and in particular
on the burn time, could lead to a large footprint of possible impact points,
which may pose a concern if the reentry points are close to inhabited regions.
Thus, the proposed guidance strategy incorporates a neutral axis maneuver
(NAM) that minimizes the sensitivity of the impact point to uncertain engine
performance. Unlike traditional methods to design a NAM, which are
particularly burdensome and require long validation and verification tasks,
the presented MPC algorithm autonomously determines the neutral axis direction
by repeatedly solving an optimal control problem (OCP) with two return phases,
a nominal and a perturbed one, constrained to the same splash-down point. The
OCP is transcribed as a sequence of convex problems that quickly converges to
the optimal solution, thus allowing for high MPC update frequencies. Numerical
results assess the robustness and performance of the proposed algorithm via
extensive Monte Carlo campaigns.

##  1 Introduction

Launch vehicle dynamics are subject to significant uncertainties, due to hard-
to-predict aerodynamic coefficients, scattering of propulsion system
performance, and sudden variations of the local environment, to name a few.
Therefore, a robust guidance algorithm that ensures the accurate injection of
the payload into the desired orbit while meeting the numerous system
requirements, even in presence of model uncertainties, is an absolute
necessity for the success of a launch vehicle (LV) ascent mission. Also, since
launch vehicles, as most aerospace systems, must operate at the boundaries of
their performance envelope to minimize the overall mission cost, their
guidance should be _optimal_ , meaning that the vehicle should fly
trajectories that are not only safe but also convenient (e.g., in terms of
propellant consumption, to maximize the carrying capacity and the
responsiveness in case of off-nominal conditions).

Model predictive control (MPC) is one of the few control synthesis strategies
that can optimize the system performance while systematically accomodating
mission constraints [1]. Specifically, MPC consists of solving repeatedly an
optimal control problem (OCP), updated with the onboard system measurements,
and implementing the computed optimal control law in the time frames between
the optimization procedures. The closed-loop architecture inherently provides
robustness to model uncertainties and in-flight disturbances, as the
continuous update of the optimal path compensates for deviations from the
nominal one.

Among the aerospace community, the interest toward MPC has grown constantly
over the years [2], featuring successful employment over diverse problems,
including rocket landing [3, 4], spacecraft landing [5, 6], rendezvous [7, 8],
and formation flying [9, 10], to name a few. Indeed, since aerospace vehicles
are typically subject to tight mission requirements, equipped with limited
computational resources, and often operate in uncertain conditions, MPC
represents an uniquely effective tool to control such systems. The earliest
applications of MPC were limited to relatively simple problems, as the
computational cost associated with the solution of the OCP can impair the
effectiveness of the MPC scheme, which can compensate for random disturbances
only if the update frequency is sufficiently high. For this reason, applying
the appealing properties of MPC to nonlinear systems has always been a
challenging task. However, in the last decades, the advances in computing
hardware and optimization theory gave a new drive to _nonlinear MPC_ [11],
which still poses serious challenges in terms of robustness guarantees and
computational efficiency, but ongoing research efforts provide meaningful
theoretical results [12] and numerical evidence over numerous successful
applications [13].

A quite common approach to nonlinear problems is _explicit MPC_ [14], which
relies on a set of precomputed solutions rather than solving onboard the OCP.
However, despite successful applications to aerospace problems are documented
in the present literature [15, 16], explicit MPC is an intrinsically
suboptimal methodology and it may require formidable resources when dealing
with high-dimensional OCPs. Recently, machine learning techniques have been
also proposed to address nonlinear guidance and control problems [17]. For
instance, a deep neural network (DNN) can be trained to imitate an expert
behavior, as in behavioral cloning [18], or learn an optimal control policy
from repeated interactions with the environment, as in reinforcement learning
[19, 20]. However, despite the low evaluation times and the high accuracy in
function approximation of DNNs, the absence of theoretical guarantees and the
extensive computational effort associated with the training, still limit their
application to much easier problems than the one here investigated.

On the other hand, convex optimization techniques are natural candidates for
real-time applications, as convex OCPs can be solved by means of highly
efficient interior point algorithms, which converge to the optimal solution in
polynomial time regardless of the initialization point [21]. The main
challenge consists in formulating a nonlinear (specifically, nonconvex)
problem as an equivalent convex one. Lossless convexification techniques are
extremely powerful in this sense, as they can transform problems without
introducing any approximations thanks to a convenient change of variables and
suitable constraint relaxations that are theoretically guaranteed to retain
the same solution as the original problem [22]. When no lossless
convexification exist for the problem at hand, one can always rely on
successive convexification methods (first and foremost, successive
linearization), which define a sequence of convex problems that eventually
converges to the optimal solution of the original problem [23]. Convergence of
successive convexification methods is not guaranteed in general, except under
appropriate assumptions [24, 25, 26], but extensive numerical evidence
suggests the validity of successive convexification in a broad range of
applications, including rocket ascent [27, 28], spacecraft rendezvous [29],
proximity operations [30], low-thrust transfers [31], rocket landing [32, 33,
34], and atmospheric entry [35, 36, 37].

This paper investigates the use of MPC for the computational guidance of the
upper stages of a four-stage VEGA-like vehicle, which features three solid
rocket motors (SRMs) and a small liquid rocket engine, AVUM, that performs the
final orbit insertion maneuver [38]. Due to this stage configuration, this
case study poses some unique challenges. Indeed, the third stage burns out at
an extremely high velocity, close to the orbital one, and ends up falling very
distant from the launch site. As a result, the return trajectory of the spent
stage must be accurately predicted and constrained to an uninhabited region,
which is usually the sea. To obtain a valuable estimate of the splash-down
location, a complete simulation of the return trajectory must be included in
the OCP, greatly increasing its complexity.

Uncertainties on the solid rocket motor performance, and especially on its
burn time, do not allow to easily pinpoint a splash-down location but rather
define a finite-dimension footprint of possible impact points. Thus, a further
system requirement is bounding the extent of this region. To this aim, robust
optimization and robust MPC can be used to endow the guidance and control with
some robustness guarantees. However, these methods are usually overly
conservative since they are primarily based on min-max OCP formulations [39]
or on constraint tightening [40]. On the other hand, more recent tube-based
MPC [41] or stochastic optimal control methods, such as chance-constrained
optimization [42] or covariance control [43, 44], could also be considered as
viable options, but their application to the case of non-additive
disturbances, such as those that arise from dynamics under uncertain time-
lengths, can be quite challenging.

For all these reasons, VEGA currently relies on a neutral axis maneuver to
robustly minimize the footprint size and make the return as insensitive as
possible to the third stage performance dispersions [45]. This maneuver is
based on the _null miss condition_ developed for ballistic missiles [46]: over
the last few seconds of operation, the stage is constrained to hold on to an
attitude such that the splash-down point is retained regardless of any
additional velocity increments. The maneuver reduces the carrying capacity of
the launch vehicle since part of the propulsive energy is spent in a non-
optimal direction, but it robustly guarantees that the actual splash-down
location is sufficiently close to the predicted one, despite additional
unpredicted seconds of operation of the SRM. However, the neutral axis
maneuver can be quite difficult to design with traditional methods and usually
extensive trajectory validation and verification tasks are required before
each launch. As an alternative to the neutral axis maneuver, the use of retro-
rockets to ignite after the third stage separation has been proposed [47].
However, such a solution would need the integration of additional hardware and
mass into a consolidated architecture.

This paper presents an improved version of the MPC algorithm previously
proposed by the authors in Ref. [48]. The main difference with the past work
concerns the introduction of a constraint on the control direction at the
stage burnout, which must correspond to the neutral-axis attitude. In
practice, two return trajectories (a nominal and a perturbed one) are
considered, whose initial state (i.e., conditions at motor cut-off) differs
only by an excess of velocity provided in the control direction. By
constraining the impact point of the two return phases to the same value, the
neutral axis attitude is indirectly attained. The envisaged approach overcomes
traditional issues in the design of the neutral axis maneuver, which relies on
cumbersome and time-consuming design procedures, and allows its integration in
a modern computational guidance scheme, ensuring the system robustness to
extra burn seconds of the SRM.

The paper is organized as follows. Section 2 describes the guidance and
control problem, detailing the phases, the dynamical model, and mission
requirements of the OCP to solve. Section 3 illustrates the state-of-the-art
convexification techniques and the $hp$ pseudospectral method used to cast the
nonlinear OCP as a sequence of convex problems that rapidly converges to the
optimal solution. In Section 4, the convex optimization algorithm is embedded
into the MPC architecture and a few relaxations are introduced to ensure the
recursive feasibility of the algorithm. Finally, Section 5 analyzes the
effectiveness of the proposed guidance through a Monte Carlo analysis. In
particular, the robustness of the employed algorithm is verified in the
presence of off-nominal initial conditions, generic uncertainties on the
dynamical model, and random SRM performance. A conclusion section ends the
manuscript.

##  2 Problem Description

MPC is based on solving at every control step an OCP. In this section, the
guidance problem of the third stage of the VEGA-like launch vehicle is
described and formulated as a multi-phase OCP. First, the phase sequence of
the OCP is outlined; then, the considered dynamical model is presented; and
finally, the objective function and all the constraints of the OCP are
reported.

###  2.1 Phase Sequence

Figure 1: Phases of the optimal control problem.

Figure 1 illustrates schematically the considered phase sequence. The phases
are numbered from 1 to 7 in a chronological order, while Phases 8 and 9
simulate the nominal and perturbed uncontrolled return of the spent third
stage and start at the burnout of the third stage. The third stage operation
is split into three phases to account for the mass discontinuity associated
with the fairing jettisoning, which takes place few seconds after the stage
ignition and divides Phase 1 from Phase 2, and for the transition to the
neutral axis maneuver, between Phases 2 and 3. In Phase 3, the control rate is
constrained to a small value, to accomodate the neutral axis maneuver. The
introduction of Phase 3 near the end of the stage operation prevents numerical
issues to the discretization of rapidly changing dynamics since, during the
transition to the neutral axis attitude, the rocket rotates by approximately
$40\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$ in a few seconds.

A short coasting (Phase 4) follows the third stage separation, then the fourth
stage performs the orbit insertion maneuver, which is composed of two
propelled arcs (Phases 5 and 7) separated by a long coasting (Phase 6). Two
return phases are included in the OCP to robustly ensure the splash-down
constraint even in the case of additional burn seconds of the SRM. The initial
position and velocity of the nominal return (Phase 8) correspond to the state
values at the third stage burnout, while the perturbed return (Phase 9) is
assumed to start with a larger velocity. The magnitude of the velocity
perturbation is fixed and its value is in the order of the maximum extra
$\Delta V$ that can be provided by the SRM, while the direction is constrained
to be the same as the control direction at the burnout. Note that since both
returns are constrained to fall at the same splash-down point, the optimal
direction of the $\Delta V$ is the neutral axis direction.

Hereinafter, let $t_{0}^{(i)}$ and $t_{f}^{(i)}$ denote the initial and final
time of the $i$-th phase. 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}^{(7{})}}$, respectively. Likewise, let $t_{R}$ denote the
nominal return time of the spent stage $t_{f}^{(8{})}$ and $t_{PR}$ the
perturbed return time $t_{f}^{(9{})}$.

###  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, with the $x$
and $y$ axes in the equatorial plane and the $z$ axis 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 [49], is the fact that any trajectory can be studied, even
missions to polar and high inclination orbits, which are the typical targets
for VEGA, without running into singularities.

The only forces assumed to act on the launch vehicle are gravity, the
aerodynamic drag, and the engine 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}SC_{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 the
system at hand.

The thrust $T$ generated by the SRMs and the associated mass flow rate
$\dot{m}_{e}$ generally vary over time according to the geometry of the solid
propellant and, in general, to the design of the stage. So, every SRM is
characterized by a vacuum thrust law $T_{\text{vac}}(t)$ and a corresponding
mass flow rate law $\dot{m}_{e}(t)$. Although upper stages operate at high
altitudes, the net thrust is slightly different from the vacuum one due to the
the atmospheric pressure $p$ contribution, thus $T=T_{\text{vac}}(t)-pA_{e}$,
where $A_{e}$ is the nozzle exit area.

Instead, the thrust and mass flow rate of the last stage, which is equipped
with a small liquid rocket engine, are assumed to be constant and equal to the
maximum attainable values. Also, differently from the other stages, AVUM can
be cut off and re-ignited. So, its burn time can be actively controlled and is
determined by the optimization process with the aim of minimizing the consumed
propellant mass.

In the propelled phases, the thrust direction $\bm{\hat{T}}$ is a control
variable and its components are expressed in the same frame as the velocity.
Being a unit vector, the components satisfy the following identity at any time

$\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 Optimal Control Problem

The goal of the optimization is to minimize the propellant consumed during the
fourth stage operation. Indeed, if the guidance algorithm minimizes the
propellant mass, not only it will provide better responsiveness in off-nominal
conditions, but it will also increase the launch vehicle nominal payload, as
smaller fuel tanks will lead to an increment of the vehicle carrying capacity.
The objective can be equivalently formulated as maximizing the final mass or
minimizing its opposite. Thus, a convenient cost function is

$J=-m(t_{f})$ (5)

In an MPC framework, the initial condition of the OCP is updated at every
control step with the real-time measurements $\tilde{\bm{x}}_{0}$, therefore
the initial state is completely assigned

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

The final conditions of Phase 7 are the following

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

In the present study, we consider a circular target orbit with semi-major axis
$a_{\text{des}}$, thus Eq. (8) with $v_{\text{des}}=\sqrt{\mu/a_{\text{des}}}$
and Eq. (9) 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}}}$.

To constrain the splash-down of the spent stage, the following conditions are
imposed at the end of Phases 8 and 9

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

with $z_{R,\text{des}}=R_{E}\sin\varphi_{R,\text{des}}$. Equations (11) and
(12) constrain the final radius to be equal to the Earth radius $R_{E}$.
Instead, Eq. (13) 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, this 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
would be very expensive since it would require a maneuver to change the plane
of the ascent trajectory.

The phases of the OCP must be linked by some conditions at the interfaces. In
particular, all state variables are continuous across phases, with the
relevant exception of the mass, which is discontinuous at the fairing
jettisoning and at the third stage separation. Thus,

$\displaystyle m(t_{0}^{(2{})})$
$\displaystyle=m(t_{f}^{(1{})})-m_{\text{fairing}}$ (14) $\displaystyle
m(t_{0}^{(4{})})$ $\displaystyle=m(t_{f}^{(2{})})-m_{\text{dry},3}$ (15)

As for Phases 8 and 9, their initial position and velocity are equal to the
burnout conditions of the third stage, with a velocity increment $\bm{\Delta
V}_{NA}$ for Phase 9, while their mass is equal to the dry mass of the stage:

$\displaystyle\bm{r}(t_{0}^{(8{})})$
$\displaystyle=\bm{r}(t_{0}^{(9{})})=\bm{r}(t_{f}^{(3{})})$ (16)
$\displaystyle\bm{v}(t_{0}^{(8{})})$ $\displaystyle=\bm{v}(t_{f}^{(3{})})$
(17) $\displaystyle\bm{v}(t_{0}^{(9{})})$
$\displaystyle=\bm{v}(t_{f}^{(3{})})+\bm{\Delta V}_{NA}$ (18) $\displaystyle
m(t_{0}^{(8{})})$ $\displaystyle=m(t_{0}^{(9{})})=m_{\text{dry},3}$ (19)

We also impose a continuity constraint between the velocity perturbation
$\bm{\Delta V}_{NA}$ and the final control direction at the nominal burnout of
the third stage. Thus,

$\bm{\Delta V}_{NA}/\Delta V_{NA}=\bm{\hat{T}}(t_{f}^{(3{})})$ (20)

This ensures that in case of an SRM longer burn time, the launch vehicle would
hold on to this direction, thus minimizing the shift of the impact point.

Finally, since the fairing is jettisoned as soon as possible to reduce the dry
mass of the system, from that moment on, the thermal flux on the rocket must
be lower than a given value to preserve the integrity of the payload. Thus, a
constraint on the maximum heat flux is included in the formulation for Phases
2 to 7

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

##  3 Convex Formulation

This section outlines the convexification strategy devised to convert the OCP
into a sequence of convex problems that quickly converges to the optimal
solution. This is a crucial step in the design of the guidance algorithm, as
it greatly reduces the cost of the onboard optimization and enables the MPC
controller to attain high update frequencies.

First, lossless convexification strategies, such as changes of variables and
constraint relaxations, are employed to reduce the nonconvexity of the
original problem. Instead, the remaining nonconvexities are tackled via
successive linearization. Safeguarding modifications, such as virtual controls
and a trust region on the time-lengths of the phases, are also included in the
formulation. The resulting convex OCP is discretized via a $hp$ Radau
pseudospectral method to set up a finite-dimensional optimization problem.
Finally, the reference solution update strategy and the convergence criteria
are outlined.

###  3.1 Change of Variables

In Eqs. (2)–(4), the control variable $\bm{\hat{T}}$ is coupled with the state
variable $m$. This may lead to numerical issues, such as high-frequency
jitters in the control law [35]. Instead, the problem should be formulated in
a way such that the dynamics are linear in the controls, for instance

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

This can be achieved by introducing a new control

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

which includes both the thrust-to-mass ratio $T/m$ and the thrust direction
vector $\bm{\hat{T}}$. By replacing it in Eqs. (2)–(4), control-affine
dynamics are obtained

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

and the $\tilde{B}$ matrix in Eq. (22) 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$.

Also the new control variables must satisfy Eq. (1), which is reformulated as

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

Note that we introduced an additional variable $u_{N}$, constrained as

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

The introduction of this nonlinear path constraint is the price to pay to
obtain control-affine dynamics. Nevertheless, the advantages of this change of
variables outweigh the disadvantages, since, despite the inclusion of the
additional constraint, the complexity of the problem is reduced compared to a
coupling of state and control variables in the dynamics.

###  3.2 Constraint Relaxation

Equation (27) is a nonlinear equality constraint that must be convexified.
This constraint belongs to a class of constraints that is suitable for a
lossless relaxation [32, 22, 50, 28]. So, by replacing the equality with an
inequality, a second-order cone constraint is formulated

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

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. (29) is
satisfied with the equality sign. A theoretical proof of the lossless property
of this relaxation can be found in Proposition 1 of Ref. [28]. It should be
noted that Proposition 1 holds as long as the heat flux path constraint (21)
is inactive almost everywhere (i.e., it is active, at most, at isolated points
in time). However, this is a mild assumption for the problem at hand, as the
heat flux profile is rarely active over arcs of finite duration.

###  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

The equations of motion (24)–(26) are affine in the control but still
nonlinear in the state; thus, they must be linearized. First, we exploit a
common strategy in direct optimization methods to account for free-time phases
[51]. 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 the
following

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

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)}$ (31)

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

The dynamics are thus 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}$ (32)

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)$
(33) $\displaystyle B$
$\displaystyle=\bar{s}\frac{\partial\bm{f}}{\partial\bm{u}}(\bar{\bm{x}},\bar{\bm{u}},\tau)$
(34) $\displaystyle\Sigma$
$\displaystyle=\bm{f}(\bar{\bm{x}},\bar{\bm{u}},\tau)$ (35)
$\displaystyle\bm{c}$ $\displaystyle=-(A\bar{\bm{x}}+B\bar{\bm{u}})$ (36)

Finally, since the linearization may lead to artificial infeasibility [25], a
virtual control signal $\bm{q}$ is added to the dynamics

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

This guarantees the feasibility of the convex problem. Naturally, as the
virtual controls are unphysical variables, their use must be highly penalized.
Thus, the following penalty term is added to the cost function

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

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 (7)–(11) 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}$ (39)
$\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}$ (40)
$\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$ (41)
$\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}}$ (42)

Likewise, the terminal radius constraints at the splash-down, Eqs. (11) and
(12), are formulated as

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

with $i=8{},9{}$.

The linearization of the above constraints may cause artificial infeasibility.
Thus, virtual buffer zones are introduced. In particular, Eqs. (39)–(43) 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}$ (44)

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

####  3.3.3 Path Constraints

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

$\frac{\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}})}{\dot{Q}_{\text{max}}}\leq
1$ (45)

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}}$
(46) $\displaystyle\frac{\partial\dot{\bar{Q}}}{\partial\bm{v}}$
$\displaystyle=\frac{3}{2}\bar{\rho}\bar{v}_{\text{rel}}\bm{\bar{v}}_{\text{rel}}$
(47)

Note that to avoid numerical issues, we normalized the constraint with respect
to the maximum heat flux.

Likewise, the auxiliary variable $u_{N}$ must satisfy Eq. (28) 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}})$
(48)

###  3.4 Trust Region on Time-Lengths

Thanks to the change of variables (23), $\bm{f}$ is linear in the control and
the $A$ and $B$ matrices do not depend on $\bar{\bm{u}}$. However, $\Sigma$
and $\bm{c}$ are functions of the reference controls and the successive
convexification algorithm may suffer from instability issues, as intermediate
controls can change significantly in the first iterations [29]. Nevertheless,
when $s\approx\bar{s}$, Eq. (32) reduces to

$\bm{x}^{\prime}=A\bm{x}+B\bm{u}+\tilde{\bm{c}}$ (49)

with $\tilde{\bm{c}}=\bar{s}\tilde{\bm{f}}(\bar{\bm{x}},\tau)-A\bar{\bm{x}}$,
thus entirely dropping the dependence on $\bar{\bm{u}}$.

To prevent excessive changes from the reference value, a soft trust region
constraint is imposed on the time-lengths $s$

$|s^{(i)}-\bar{s}^{(i)}|\leq\delta^{(i)}\qquad i=6{},7{}$ (50)

This constraint is imposed only on the duration of Phases 6 and 7 since only
these arcs resulted particularly sensitive to instability phenomena [28]. 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=6{},7{}$ (51)

In the authors’ experience, a suitable choice of the upper bound
$\delta_{\text{max}}^{(i)}$ is usually somewhere between 1% and 10% of
$\bar{s}^{(i)}$. Instead, the penalty weights $\lambda_{\delta}$ should be
picked as small as possible in order not to converge toward suboptimal
solutions (e.g., in the range $\smash{10^{-6}}\div\smash{10^{-3}}$).

###  3.5 Discretization

The problem formulated so far is infinite-dimensional since state and controls
are continuous-time functions. However, to solve it numerically, the problem
must be cast as a finite set of variables and constraints. In this paper, we
use a $hp$ pseudospectral method, which is a particularly performing
transcription method since it combines the advantages of $h$ and $p$ schemes.
Indeed, pseudospectral methods exhibit exponential convergence rate only in
regions where the solution is smooth, but combining them with $h$ methods
allows to introduce mesh nodes near potential discontinuities, thus extending
their appealing properties to problems with rapidly changing dynamics [52].
These features are key to reduce the computational cost of the onboard
optimization, as such a discretization method limits the problem dimension and
yet it accurately represents the continuous-time problem.

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) [53] is used as collocation scheme since it is one
of the most accurate and performing pseudospectral methods [54]. Furthermore,
the Legendre-Gauss-Radau (LGR) points avoid redundant control variables at the
segment boundaries and approximate the optimal control at each mesh point,
except for the final node of the final segment that must be extrapolated _a
posteriori_ from the polynomial approximation of the control signal.

Once a grid is defined, the state and control are discretized over it,
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 (37) 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$ (52)

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 [55].

Note that also the virtual control $\bm{q}$ and the linearization matrices
(33)–(36) are discretized over the mesh. We now define also the penalty
function $P(\bm{q})$ from Eq. (38) as the 1-norm of the vector containing all
the discrete virtual control variables $\bm{q}_{i}^{k}$. Likewise, path
constraints (45) and (48) are imposed at every node, generating a finite set
of algebraic constraints. Further details on the implementation of the $hp$
Radau pseudospectral method can be found in Ref. [56].

###  3.6 Update of the Reference Solution

Successive linearization is based on Taylor series expansions around a
reference solution. After each convex problem is solved, the reference
solution must be updated somehow. The most common update method simply
consists of replacing the reference solution with the new one. However, to
provide further algorithmic robustness to the successive convexification
approach, we employ an alternative approach, referred to as _filtering_ , that
was successfully used by the authors in previous works [29, 28, 48].

Filtering consists of updating the reference solution with a weighted sum of
the previous solutions. By averaging out the values over different solutions,
the sequence of reference solutions is smoother and instabilities due to
possible diverging intermediate iterations are filtered out. As a result, the
algorithm is less sensitive to artificial unboundedness and trust regions [25]
or other expedients [37], typically necessary in such a complex problem, can
be avoided. Indeed, the major appeal of filtering is the fact that it does not
introduce any additional constraints or penalty terms to the problem,
retaining its original formulation.

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)\\}}$ (53)

where $K$ is the number of previous solutions to account for, $\alpha_{k}$ are
the corresponding weights, and $x^{(i)}$ denotes the solution to the $i$-th
problem. Note that until $i\geq K$ the initial reference solution $x^{(0)}$
appears multiple times in the sum. Three previous solutions are used for the
reference solution update, weighted in the same way as in the authors’
previous works [29, 28, 48]. The values of the weights are $\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}}$
(54)

  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}$
(55)

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}$ (56)

##  4 Model Predictive Control

###  4.1 Guidance Strategy

In the present work, the MPC controller is applied only to the third stage
operation, thus the control cycle repeats until the stage burnout, which is
the MPC stopping criterion. However, the _prediction_ horizon (i.e., the time
domain of the OCP) must last until the first point in time that allows to
evaluate performance and account for all the mission constraints, that is, the
payload release into orbit and the splash-down of the third stage.

The peculiarity of the present application is a possible mismatch between the
expected and the actual burnout time of the stage. To design a guidance
algorithm that is robust to this uncertainty, we devised a strategy that
consists of a closed-loop guidance until the minimum burn time and an open-
loop neutral axis maneuver in the extra burn seconds. The neutral axis
direction is determined onboard by the OCP, which provides the optimal
direction of an additional $\Delta V$ such that the return point is unchanged.

###  4.2 Noise Model

To simulate the actual stage operation, a numerical integration of the
original dynamical model, Eqs. (2)–(4), is carried out. To account for model
uncertainty and external disturbance, white Gaussian noise is added to the
dynamics, resulting in the following stochastic differential equation (SDE)

$d\bm{x}=\bm{f}(\bm{x},\bm{u},t)dt+G(\bm{x},\bm{u},t)d\bm{w}$ (57)

where $\bm{w}$ is a $n_{w}$-dimensional Wiener process and the matrix $G$
determines how the external disturbances affect the system.

In the present case, we assume that the standard deviation of the random
disturbance is proportional to the thrust level by a factor $\alpha_{G}$.
Also, we model two different noise contributions: the first one is aligned
with the thrust direction while the other one is in a random direction. We
suppose that the former is $\alpha_{G_{\parallel}}$ times larger than the
latter, as most of the noise is usually in the thrust direction. The resulting
$G$ matrix is

$G=\alpha_{G}\frac{T}{m}\begin{bmatrix}\bm{0}_{3\times 1}&\bm{0}_{3\times
3}\\\ \alpha_{G_{\parallel}}\bm{\hat{T}}&\bm{I}_{3\times 3}\\\ \bm{0}_{1\times
1}&\bm{0}_{1\times 3}\end{bmatrix}$ (58)

Note that $n_{w}=4$.

Also the initial conditions at the stage ignition are perturbed to account for
possible performance dispersions of the previous stages. The position error is
modeled as a perturbation of the altitude only, as perturbations in other
directions do not significantly affect performance. Errors on $r$ are
uniformly sampled in the range $[-\Delta_{r},\Delta_{r}]$. Instead, the
scattered initial conditions on the velocity are generated by uniformly
sampling noise from a sphere of radius $\Delta_{v}$.

As for the SRM thrust profile used in the simulation, every simulation
considers a perturbed thrust law that features an additional burn time $\Delta
t_{b}$ of up to $5\text{\,}\mathrm{s}\text{/}$ and a small increase in the
engine total impulse $\Delta I_{\text{tot}}$ up to $0.5$% with respect to the
nominal one used as a reference in the OCP. By randomly scattering these
quantities, a constant, nonzero thrust level is assumed in the extra seconds
of operation

$T_{\text{vac}}(t>t_{b,\text{nom}})=\frac{\Delta I_{\text{tot}}}{\Delta
t_{b}}$ (59)

The nominal thrust profiles and some examples of perturbed laws are shown in
Fig. 2.

Figure 2: Nominal and perturbed thrust profiles

###  4.3 Relaxation

Recursive feasibility of the OCP is key to the design of an MPC controller, as
the optimization procedure must provide a control signal at every control
step. The convex problem formulated in Section 3 already includes virtual
controls and buffer zones to ensure the feasibility of each instance. Under
nominal circumstances, these expedients would guarantee the feasibility of the
OCP. However, when the system operates in off-nominal conditions, satisfying
some constraints may be impossible, and, as a result, the optimization may
provide unphysical solutions that actively exploit virtual variables or no
solution at all. A common strategy to avoid this undesired phenomenon is
relaxing the initial condition constraint (6), e.g., as proposed in Ref. [4].
This is a general approach that can be readily employed regardless of the
application. However, as such, it does not exploit any knowledge on the
problem specificities. Rather, if possible, it is more efficient to identify
the critical mission requirements and formulate a properly relaxed problem.

For the problem under investigation, only the constraint on maximum thermal
flux (45) may cause an infeasible problem instance. Indeed, the nominal
trajectory is such that the heat flux is greater than or equal to the
threshold just until the fairing jettisoning. So, 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 keeping it below
threshold may be impossible. The infeasibility is solved by relaxing the heat
flux constraint (45) as

$\frac{\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}})}{\dot{Q}_{\text{max}}}\leq
1+\delta_{\dot{Q}}$ (60)

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}}$ (61)

In the authors’ experience, high values should not be assigned to the penalty
weight $\lambda_{\dot{Q}}$, since the optimization would then exploit virtual
controls and may never meet the converge criterion (55). 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.

###  4.4 Update

At the end of each control cycle, the OCP must be updated with the
measurements coming from the navigation system. The update concerns the
initial condition of the OCP, Eq. (6), which is updated with the state vector
resulting from the numerical integration of the SDE (57). Also, the 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 is also reduced
linearly at every step from the values in Table 1 until a minimum value of $5$
nodes per segment. Also, Phases 1 and 2 are removed from the OCP as soon as
they are over.

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

As soon as the nominal burn time $t_{b,\text{nom}}$ is passed, the guidance
switches from a closed-loop MPC architecture to an open-loop neutral axis
maneuver for the extra seconds of stage operation. The control direction
during the neutral axis maneuver is constant and corresponds to the terminal
control action computed in the last MPC cycle.

It is worthwhile mentioning that, in the present work, the time to solve the
optimization problem is supposed null, but actually the simulation should
introduce a delay between the measurement of the updated state and the
actuation of the optimal control law. This delay can be quite long and
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 Numerical Results

In this section, numerical results are presented to assess the performance and
robustness of the devised MPC framework. In particular, we compare the devised
guidance strategy with the one previously proposed by the authors [48]. We
refer to the latter as _single-return_ algorithm and to the novel one as
_multi-return_ since it includes an additional perturbed return phase in the
OCP to autonomously determine the neutral axis direction. The algorithms were
implemented in C++ and use Gurobi’s [57] second-order cone programming solver
to solve the OCP. All the simulations were run on a computer equipped with
Intel® Core™ i7-9700K CPU @ $3.60\text{\,}\mathrm{GHz}\text{/}$.

The parameters of the successive convexification algorithm are the same as in
Ref. [28]. Specifically, the penalty weights on the time-lengths,
$\smash{\lambda_{\delta}^{(6{})}}$ and $\smash{\lambda_{\delta}^{(7{})}}$,
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. (54)–(56), were prescribed as
$\epsilon_{\text{tol}}=${10}^{-4}$$ and $\epsilon_{f}=${10}^{-6}$$.

###  5.1 Nominal Trajectories

The investigated scenario is the same as in Refs. [28] and [48], 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 constrained to
$\varphi_{R,\text{des}}=$60\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$$.
Also, the same data are used to model the launch vehicle stages and
aerodynamic coefficients, with the exception of the uncertain additional burn
time of the third stage, which is modeled as detailed in Section 4. The U.S.
Standard Atmosphere 1976 model is used to evaluate the air density and
pressure as functions of the altitude [58].

The third stage nominal burn time is
$t_{b,3}=$110\text{\,}\mathrm{s}\text{/}$$, which is also the minimum burn
time considered in the simulations. The fairing is released $\Delta
t^{(1{})}=$5.4\text{\,}\mathrm{s}\text{/}$$ after ignition, while Phase 3 is
prescribed to last $\Delta t^{(3{})}=$5\text{\,}\mathrm{s}\text{/}$$. Thus,
Phase 2 lasts $\Delta t^{(2{})}=$99.6\text{\,}\mathrm{s}\text{/}$$. The
duration of the coasting arc after the third stage separation (Phase 4) is
fixed to $\Delta t^{(4{})}=$15.4\text{\,}\mathrm{s}\text{/}$$. The time-
lengths of all the following phases are left to be optimized.

In the OCP, the threshold on the heat flux used to compute the solution is set
to $\dot{Q}_{\text{max}}=$900\text{\,}\mathrm{W}\text{/}{\mathrm{m}}^{2}$$.
Since $\dot{Q}_{\text{max}}$ is lower than the maximum bearable heat flux that
must not be exceeded in simulations
($1135\text{\,}\mathrm{W}\text{/}{\mathrm{m}}^{2}$) [38], the relaxation of
the heat flux constraint with a relatively small value of the penalty
coefficient ($\lambda_{\dot{Q}}=${10}^{-2}$$) is justified.

The nominal trajectories of the single and multi-return algorithms were
computed via the convex optimization approach described in Ref. [28], with the
only difference of including Phases 3 and 9 in the OCP for the multi-return
scenario. The considered velocity increment at the start of Phase 9 is $\Delta
V_{NA}=$25\text{\,}\mathrm{m}\text{/}\mathrm{s}$$, which is comparable to the
(expected) maximum velocity increment provided by the SRM compared to the
nominal performance. The nominal trajectories are used as initial reference
solution in the first control cycle and provide the nominal payload masses:
$1400.1\text{\,}\mathrm{kg}\text{/}$ for the single-return problem and
$1328.3\text{\,}\mathrm{kg}\text{/}$ for the multi-return one. The carrying
capacity of the multi-return strategy is smaller due to the fact that the
third stage must be oriented in the neutral axis direction at burnout, while
the single-return problem does not take into account such requirement. Thus,
the difference in the payload mass represents the price to pay for attaining
the neutral axis attitude.

###  5.2 Monte Carlo Campaigns

A Monte Carlo analysis on the combined effect of off-nominal initial
conditions, in-flight disturbance, and uncertain thrust profiles as described
in Section 4.2 was carried out for both the single-return strategy and the
present one. The initial radius error is scattered uniformly in a range of
size $\Delta_{r}=$500\text{\,}\mathrm{m}\text{/}$$ and the initial velocity
error is sampled from a sphere of radius
$\Delta_{v}=$40\text{\,}\mathrm{m}\text{/}\mathrm{s}$$. Three levels of in-
flight disturbance were considered as reported in Table 2 and called Low (L),
Medium (M), and High (H). In each case, the standard deviation of the noise in
the thrust direction is 5 times greater than the one in a random direction,
thus $\alpha_{G_{\parallel}}=5$.

Table 2: Standard deviations of the random-direction Gaussian in-flight disturbance in terms of $T/m$ percentage Parameter | Case L | Case M | Case H
---|---|---|---
$\alpha_{G}$ | $0.25\text{\,}\mathrm{\text{\textperthousand}}\text{/}$ | $0.5\text{\,}\mathrm{\text{\textperthousand}}\text{/}$ | $1\text{\,}\mathrm{\text{\textperthousand}}\text{/}$

In each simulation, the additional burn time $\Delta t_{b}$ is sampled
uniformly in the range $[0,5]$$\mathrm{s}\text{/}$, while the additional total
impulse $\Delta I_{\text{tot}}$ is a fraction of the nominal one uniformly
sampled in the range $[0,5]$$\mathrm{\text{\textperthousand}}\text{/}$. The
update frequency of the MPC is set to $1\text{\,}\mathrm{Hz}\text{/}$, meaning
that the OCP is solved every $T=$1\text{\,}\mathrm{s}\text{/}$$.

Table 3: Results of the Monte Carlo campaigns Algorithm |  Case | Payload mass (kg) |  $\varphi_{R}$ (deg)
---|---|---|---
Min | Mean | Max | $\sigma$ | Min | Mean | Max | $\sigma$
Single return | L | $1359.19$ | $1409.74$ | $1449.70$ | $17.10$ | $59.74$ | $63.60$ | $68.15$ | $1.87$
M | $1335.17$ | $1405.21$ | $1482.21$ | $24.41$ | $59.30$ | $63.33$ | $68.65$ | $2.00$
H | $1254.00$ | $1407.49$ | $1503.37$ | $39.77$ | $57.78$ | $63.63$ | $75.40$ | $2.73$
Multi return | L | $1278.39$ | $1331.44$ | $1375.78$ | $18.16$ | $59.72$ | $60.00$ | $60.21$ | $0.09$
M | $1266.97$ | $1329.81$ | $1389.62$ | $23.56$ | $59.17$ | $60.01$ | $60.62$ | $0.20$
H | $1215.90$ | $1331.64$ | $1456.62$ | $41.06$ | $58.25$ | $60.02$ | $61.41$ | $0.38$

The results of the Monte Carlo campaigns are summarized in Table 3. For each
case, 400 independent MPC simulations were carried out. For each simulation,
the final conditions at the burnout of the SRM are used to estimate the
attainable payload mass by solving the ascent OCP from Phases 4 to 7 and to
evaluate the splash-down point through forward propagation of the equations of
motion, Eqs. (2)–(4). As expected, for both algorithms, as the intensity of
the in-flight disturbance increases, larger dispersions on the payload mass
and on the splash-down point are observed. The average payload mass attained
by the algorithms is slightly greater than the nominal one due to the
additional $I_{\text{tot}}$ that provides more energy to the system. The gain
in payload mass of the multi-return strategy is smaller due to the fact that
the additional impulse is provided in the neutral axis direction, which is not
optimal. Instead, the payload mass ranges of the two algorithms are
comparable, except for a shift of approximately
$75\text{\,}\mathrm{kg}\text{/}$ due to the cost of the neutral axis maneuver.

The major difference in the behavior of the two strategies lies in the splash-
down point. The additional impulse provided by the SRM shifts the mean splash-
down point attained by the single-return strategy farther by approximately
$3.6\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$ in all three cases, meaning
that this deviation depends only on the SRM performance uncertainty. Instead,
the multi-return algorithm retains the mean splash-down location with great
accuracy (the largest error is
$0.02\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$ in case H) thanks to the
fact that the thrust vector is oriented in the neutral axis direction during
the extra seconds of operation. The dispersion of the spent stage impact
latitude grows for both approaches as the disturbance intensity increases,
but, while the lower bounds are comparable, the upper bounds of the return
points corresponding to the single-return strategy are much more distant from
the desired value. This asymmetry is due to the fact that the additional SRM
impulse moves the return point only farther; thus, the two disturbances add up
for the single-return approach. Instead, the neutral axis maneuver compensates
for the SRM performance uncertainty and the range of the attained
$\varphi_{R}$ is symmetric. Also, the dispersion of the return point is
significantly reduced when employing the multi-return algorithm, with standard
deviations $\sigma$ one order of magnitude smaller than in the single-return
approach.

(a) Single-return algorithm

(b) Multi-return algorithm

Figure 3: Altitude profiles for Case H

The trajectories for case H, which is associated with the largest envelopes,
are reported in Fig. 3 and compared with the nominal one, which is denoted by
the black line. The trajectories of both guidance strategies deviate
significantly from the nominal trajectory during the central portion of the
flight, as the MPC approach autonomously recomputes optimal paths to
efficiently compensate for the encountered disturbances. However, while the
terminal conditions at the payload release are met by both algorithms, the
single-return strategy features a much greater dispersion on the return point,
as previously commented.

(a) Single-return algorithm

(b) Multi-return algorithm

Figure 4: Control elevation profiles for Case H

The control histories are shown in Fig. 4 in terms of the elevation angle,
which is defined as the angle between the thrust direction $\bm{\hat{T}}$ and
the local horizontal. The elevation profiles are close to the nominal controls
for most of the third stage operation, with larger deviations toward the end.
Indeed, due to the external disturbance, meeting the splash-down constraint
requires rapid variations of the thrust direction as the burnout time
approaches. Figure 4 shows the different orientation of the neutral axis
direction compared to the optimal one, as the maneuver requires the rocket to
quickly rotate by approximately
$40\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$ downward.

(a) Single-return algorithm

(b) Multi-return algorithm

Figure 5: Heat flux profiles for Case H

The heat flux profiles are illustrated in Fig. 5. Despite the relaxation, the
thermal requirement is satisfied within an acceptable error. In particular,
the heat flux is never greater than
$1000\text{\,}\mathrm{W}\text{/}{\mathrm{m}}^{2}$, except for a few
simulations of the single-return strategy in the proximity of the burnout
where, due to the additional seconds of SRM operations, the relative velocity
increases more than predicted. However, the violation only depends on the
scattering of the initial conditions, as increased initial velocity and
reduced altitude expose the system to more critical thermal conditions.

As a final remark, the average solution time of the single-return OCP was
$1.50\text{\,}\mathrm{s}\text{/}$, while solving the multi-return problem
required $3.20\text{\,}\mathrm{s}\text{/}$ on average. The multi-return OCP is
thus more expensive to solve, due to the additional return phase and the
related constraints that significantly increase the problem dimension and
complexity. Even though both times are incompatible with a
$1\text{\,}\mathrm{Hz}\text{/}$ update frequency, this issue could be easily
addressed by using optimized software, dedicated hardware, and custom convex
solvers, for which the computation times are expected to be at least one order
of magnitude lower. Also, previous investigations on this topic show that good
MPC performance could be obtained even with less frequent updates, in the
range $1$–$3\text{\,}\mathrm{s}\text{/}$ [48].

##  6 Conclusions

This paper presented a model predictive control algorithm for the optimal
guidance of a launch vehicle upper stage that robustly controls the splash-
down location of the spent stage even in presence of uncertainties on the
motor burn time, by incorporating in the OCP a set of boundary conditions that
allows to numerically find the thrust direction corresponding to the null miss
condition. Specifically, two return phases of the spent stage are considered
in every OCP, a nominal return and another perturbed by a small velocity
impulse that is applied in the same direction as the thrust at the end of the
nominal cut-off. This approach avoids the solution of the neutral axis
equation, which enforces the same condition algebraically for a ballistic
return phase, and can be directly introduced in the MPC guidance scheme, thus
greatly reducing the validation and verification burden of such a complex
system.

The OCP was mindfully cast as a sequence of convex problems that converges to
the optimal solution in a short time, allowing the MPC algorithm to attain an
high update frequency and effectively compensate for external random
disturbances. Also, a $hp$-pseudospectral discretization was adopted to
provide an accurate representation of the continuous-time ascent and return
dynamics with a limited computational burden.

A Monte Carlo analysis has been carried out to investigate the effectiveness
of the proposed approach in presence of in-flight disturbances, off-nominal
operating conditions, and uncertain SRM performance. The results of this
analysis provide significant evidence of the robustness of the MPC algorithm.
In particular, the benefits in terms of robustness of the splash-down location
of the present approach compared to the previous one that did not incorporate
a null miss constraint into the OCP are apparent. Indeed, even though the
magnitude of the in-flight disturbances directly affects the dispersion on the
payload mass and on the splash-down location, the deviations are one order of
magnitude smaller.

Even though further tests are still necessary to demonstrate the suitability
of the proposed MPC algorithm for practical onboard implementation, the
preliminary results are very encouraging, as the design of the neutral axis
maneuver is still a burdensome, yet necessary, task for several launch
vehicles, first and foremost, VEGA, and having an autonomous strategy to
include the null miss condition within the optimization loop represents a
significant step forward in the design of launch vehicles guidance and control
algorithms.

## References

  * [1] B. Kouvaritakis and M. Cannon, Model Predictive Control.  Advanced Textbooks in Control and Signal Processing, Springer International Publishing, 2016.
  * [2] 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.
  * [3] C. A. Pascucci, S. Bennani, and A. Bemporad, “Model predictive control for powered descent guidance and control,” in European Control Conference (ECC), pp. 1388–1393, 2015.
  * [4] J. Wang, N. Cui, and C. Wei, “Optimal rocket landing guidance using convex optimization and model predictive control,” Journal of Guidance, Control, and Dynamics, vol. 42, no. 5, pp. 1078–1092, 2019.
  * [5] J. Carson and A. Acikmese, “A model-predictive control technique with guaranteed resolvability and required thruster silent times for small-body proximity operations,” in AIAA Guidance, Navigation, and Control Conference and Exhibit, 2006.
  * [6] T. Reynolds and M. Mesbahi, “Small body precision landing via convex model predictive control,” in AIAA SPACE and Astronautics Forum and Exposition, 2017.
  * [7] E. N. Hartley, “A tutorial on model predictive control for spacecraft rendezvous,” in European Control Conference (ECC), pp. 1355–1361, 2015\.
  * [8] A. Weiss, M. Baldwin, R. S. Erwin, and I. Kolmanovsky, “Model predictive control for spacecraft rendezvous and docking: Strategies for handling constraints and case studies,” IEEE Transactions on Control Systems Technology, vol. 23, no. 4, pp. 1638–1647, 2015.
  * [9] T. Keviczky, F. Borrelli, K. Fregene, D. Godbole, and G. J. Balas, “Decentralized receding horizon control and coordination of autonomous vehicle formations,” IEEE Transactions on Control Systems Technology, vol. 16, no. 1, pp. 19–33, 2008.
  * [10] D. Morgan, S.-J. Chung, and F. Y. Hadaegh, “Model predictive control of swarms of spacecraft using sequential convex programming,” Journal of Guidance, Control, and Dynamics, vol. 37, no. 6, pp. 1725–1740, 2014.
  * [11] D. Mayne, “Nonlinear model predictive control: Challenges and opportunities,” in Nonlinear Model Predictive Control, pp. 23–44, Springer, 2000.
  * [12] G. De Nicolao, L. Magni, and R. Scattolini, “Stability and robustness of nonlinear receding horizon control,” in Nonlinear Model Predictive Control, pp. 3–22, Springer, 2000.
  * [13] S. J. Qin and T. A. Badgwell, “An overview of nonlinear model predictive control applications,” in Nonlinear Model Predictive Control, pp. 369–392, Springer, 2000.
  * [14] A. Alessio and A. Bemporad, “A survey on explicit model predictive control,” in Nonlinear Model Predictive Control, pp. 345–369, Springer, 2009.
  * [15] Øyvind Hegrenæs, J. T. Gravdahl, and P. Tøndel, “Spacecraft attitude control using explicit model predictive control,” Automatica, vol. 41, no. 12, pp. 2107 – 2114, 2005.
  * [16] M. Leomanni, E. Rogers, and S. B. Gabriel, “Explicit model predictive control approach for low-thrust spacecraft proximity operations,” Journal of Guidance, Control, and Dynamics, vol. 37, no. 6, pp. 1780–1790, 2014.
  * [17] 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.
  * [18] 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.
  * [19] 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.
  * [20] 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.
  * [21] S. Boyd, S. P. Boyd, and L. Vandenberghe, Convex Optimization.  Cambridge University Press, 2004.
  * [22] 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.
  * [23] X. Liu, P. Lu, and B. Pan, “Survey of convex optimization for aerospace applications,” Astrodynamics, vol. 1, no. 1, pp. 23–40, 2017.
  * [24] X. Liu and P. Lu, “Solving nonconvex optimal control problems by convex optimization,” Journal of Guidance, Control, and Dynamics, vol. 37, no. 3, pp. 750–765, 2014.
  * [25] Y. Mao, M. Szmuk, and B. Açıkmeşe, “Successive convexification of non-convex optimal control problems and its convergence properties,” in 55th Conference on Decision and Control (CDC), pp. 3636–3641, IEEE, 2016\.
  * [26] R. Bonalli, A. Cauligi, A. Bylard, T. Lew, and M. Pavone, “Trajectory optimization on manifolds: A theoretically-guaranteed embedded sequential convex programming approach,” in Proceedings of Robotics: Science and Systems, 2019.
  * [27] 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), July 2019.
  * [28] 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,” in AAS/AIAA Astrodynamics Specialist Conference, 2020.
  * [29] 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.
  * [30] P. Lu and X. Liu, “Autonomous trajectory planning for rendezvous and proximity operations by conic optimization,” Journal of Guidance, Control, and Dynamics, vol. 36, no. 2, pp. 375–389, 2013.
  * [31] 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.
  * [32] 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.
  * [33] 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\.
  * [34] M. Sagliano, “Generalized hp pseudospectral-convex programming for powered descent and landing,” Journal of Guidance, Control, and Dynamics, vol. 42, no. 7, pp. 1562–1570, 2019.
  * [35] 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.
  * [36] Z. Wang and M. J. Grant, “Constrained trajectory optimization for planetary entry via sequential convex programming,” Journal of Guidance, Control, and Dynamics, vol. 40, no. 10, pp. 2603–2615, 2017.
  * [37] Z. Wang and Y. Lu, “Improved sequential convex programming algorithms for entry trajectory optimization,” Journal of Spacecraft and Rockets, vol. 57, no. 6, pp. 1373–1386, 2020.
  * [38] Arianespace, Vega User’s Manual, 2014.
  * [39] A. Bemporad and M. Morari, “Robust model predictive control: A survey,” in Robustness in identification and control, pp. 207–226, 1999.
  * [40] 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.
  * [41] W. Langson, I. Chryssochoos, S. Raković, and D. Mayne, “Robust model predictive control using tubes,” Automatica, vol. 40, no. 1, pp. 125–133, 2004.
  * [42] 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.
  * [43] 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.
  * [44] 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.
  * [45] M. Giannini and I. Cruciani, “VEGA LV Qualification Process: GNC aspects on HWIL Testing and Analysis,” in 5th European Conference for Aeronautics and Space Sciences, 2013.
  * [46] F. J. Regan and S. M. Anandakrishnan, Dynamics of Atmospheric Re-Entry, ch. 6, pp. 175–177.  American Institute of Aeronautics and Astronautics, 1993.
  * [47] G. Martens, E. Vellutini, and I. Cruciani, “Innovative strategy for z9 reentry,” in 6th International Conference on Astrodynamics Tools and Techniques (ICATT), 2016.
  * [48] 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.
  * [49] 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.
  * [50] X. Liu, Z. Shen, and P. Lu, “Exact convex relaxation for optimal flight of aerodynamically controlled missiles,” IEEE Transactions on Aerospace and Electronic Systems, vol. 52, pp. 1881–1892, August 2016.
  * [51] J. T. Betts, Practical methods for optimal control and estimation using nonlinear programming, vol. 19, ch. 4, pp. 132–134.  Siam, 2010.
  * [52] 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.
  * [53] D. Garg, Advances in Global Pseudospectral Methods for Optimal Control.  PhD thesis, University of Florida, Gainesville, FL, 2011.
  * [54] 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.
  * [55] J.-P. Berrut and L. N. Trefethen, “Barycentric lagrange interpolation,” SIAM Review, vol. 46, no. 3, pp. 501–517, 2004.
  * [56] M. A. Patterson and A. V. Rao, “Gpops-ii: A matlab software for solving multiple-phase optimal control problems using hp-adaptive gaussian quadrature collocation methods and sparse nonlinear programming,” ACM Trans. Math. Softw., vol. 41, no. 1, 2014.
  * [57] Gurobi Optimization, LLC, Gurobi Optimizer Reference Manual, 2020.
  * [58] NOAA, NASA, and USAF, U.S. Standard Atmosphere.  U.S. Government Printing Office, 1976.