# Robust Nonlinear Optimal Control
via System Level Synthesis

Antoine P. Leeman1, Johannes Köhler1, Andrea Zanelli1, Samir Bennani2, and
Melanie N. Zeilinger1 This work has been supported by the European Space
Agency under OSIP 4000133352, the Swiss Space Center, and the Swiss National
Science Foundation under NCCR Automation (grant agreement 51NF40
180545).1Antoine P. Leeman, Johannes Köhler, Andrea Zanelli, and Melanie N.
Zeilinger are with the Institute for Dynamic Systems and Control, ETH Zürich,
Zürich 8053, Switzerland (email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>mzeilinger@ethz.ch)2Samir Bennani is with the European Space
Agency, Noordwijk 2201AZ, The Netherlands (email<EMAIL_ADDRESS>
###### Abstract

This paper addresses the problem of finite horizon constrained robust optimal
control for nonlinear systems subject to norm-bounded disturbances. To this
end, the underlying uncertain nonlinear system is decomposed based on a first-
order Taylor series expansion into a nominal system and an error (deviation)
described as an uncertain linear time-varying system. This decomposition
allows us to leverage system level synthesis to optimize an affine error
feedback while planning the nominal trajectory and ensuring robust constraint
satisfaction for the nonlinear system. The proposed approach thereby results
in a less conservative planning compared with state-of-the-art techniques. A
tailored sequential quadratic programming strategy is proposed to solve the
resulting nonlinear program efficiently. We demonstrate the benefits of the
proposed approach to control the rotational motion of a rigid body subject to
state and input constraints.

###### Index Terms:

NL predictive control, Nonlinear systems, Optimal control, Robust control,
System level synthesis

##  I Introduction

Robust nonlinear optimal control represents one of the central problems in
many safety-critical applications, involving, e.g., robotic systems, drones,
spacecraft, and many others. While this problem has been extensively studied
in the literature [1], and rigorous constraint satisfaction properties can be
derived in the presence of disturbances (see robust predictive control
formulations, e.g., [2, 3, 4]), this is commonly achieved at the cost of
introducing conservativeness.

The robust control design task is traditionally divided into two main steps:
the optimization of the nominal trajectory [5, 6] and the offline design of a
stabilizing feedback [3] compensating for modeling errors or disturbances. To
ensure robust satisfaction of safety-critical state constraints, the nominal
trajectory optimization is coupled with the over-approximation of the error
reachable set using, e.g., tubes or funnels [7]. There exists a wide range of
techniques to construct corresponding over-approximations of the tubes/funnels
(cf., e.g., [8, 9, 10, 7, 2, 3, 4]), however, these methods can introduce
significant conservatism, especially due to the choice of an offline fixed
error feedback. We address this limitation by proposing a method to solve the
robust trajectory optimization while jointly optimizing over the error
feedback and ensuring guaranteed robust constraint satisfaction.

The conservativeness of an offline-determined error feedback policy can be
addressed for linear systems by directly predicting robust control invariant
polytopes [11]. Compare also [12] for a recent approach for nonlinear systems.
Another systematic approach to jointly optimize a linear feedback while
considering constraints is presented in [13, 14], which extends approximate
ellipsoidal disturbance propagation [10, 15] to include optimized feedback
policies. Other methods to obtain feedback policies and (optimal)
trajectories, which come without principled guarantees for robust constraint
satisfaction (cf. [16, 17, 18]), are used in practice.

However, all the above mentioned methods that jointly optimize over nominal
trajectories and feedback policies result in an over- or under-approximation
of the true reachable set, even for linear systems. We overcome this
limitation by leveraging system level synthesis (SLS) [19, 20, 21], or
equivalently affine disturbance feedback [22, 23, 24]. In particular, for
linear systems, SLS allows to jointly optimize a linear error feedback policy
and nominal trajectory and thereby provide a tight reachable set at least for
linear systems. There exist conceptual extensions of the SLS framework to
nonlinear systems [25, 26, 27], however these existing approaches do not
consider (robust) constraint satisfaction.

Disturbance and linearization error $\displaystyle
d_{k}=w_{k}+r(x_{k},u_{k},z_{k},v_{k})$ Nominal dynamics $\displaystyle
z_{k+1}=f(z_{k},v_{k})$ SLS-based LTV error reachable set $\displaystyle
x_{k}\in z_{k}\oplus\mathbb{D}_{k}$ Theorem III.1
$\displaystyle(x_{k},u_{k})\in\mathbb{P}$ $\forall w_{k}\in\mathbb{W}$ Figure
1: Decomposition of the uncertain nonlinear dynamics into nominal nonlinear
dynamics, an error term made up of the disturbance and linearization error
(Section III-A), and an LTV error system used for the SLS-based error
reachable sets (Section III-B). This decomposition enables optimization over
affine error feedback with robust constraint satisfaction for the nonlinear
uncertain system (Section III-C).

#### Contribution

We propose a novel approach for optimal control of nonlinear systems with
robust constraint satisfaction, using SLS. As shown in Fig. 1, the nonlinear
system is decomposed into a nominal nonlinear system and a linear time-varying
(LTV) error (deviation) system constructed around the online-optimized nominal
trajectory, which includes an online-optimized error term corresponding to the
linearization error (Section III-A). We apply a linear SLS formulation
(Section III-B) to the LTV error system to jointly optimize the affine error
feedback and the nominal trajectory and obtain an over-approximation of the
reachable set (funnel) for the nonlinear system (Section III-C). The presented
method has the following advantages when compared to the literature:

  * •

The nominal trajectory and the affine error feedback are optimized jointly to
decrease conservativeness.

  * •

The reachable set used for robust constraint satisfaction is tight for linear
systems, i.e., the only source of conservativeness stems from the over-
approximation of the linearization error.

In addition, we provide an inexact sequential quadratic programming (SQP)
algorithm to solve the corresponding nonlinear program (NLP), which can
substantially reduce the computation times compared to other SQP methods, or
IPOPT [28], especially when the dynamics of the system are described by an
expensive to integrate ordinary differential equation (Section IV).
Furthermore, this inexact SQP algorithm ensures that the resulting QP sub-
problems are comparable to linear SLS problems [20]. Overall, the proposed
approach can be applied to any three times continuously differentiable
nonlinear system and requires no complicated offline design.

Finally, we demonstrate the benefits of the proposed method using a nonlinear
numerical example and provide a comparison with open-loop (robust) trajectory
optimization and optimal control based on the linearized dynamics to highlight
the reduced conservativeness (Section V).

#### Notation

We define the set $\mathbb{N}_{T}\vcentcolon=\\{0,\dots,T-1\\}$ where $T$ is a
natural number. We denote stacked vectors or matrices by
$\left(a,b\right)=[a^{\top}\leavevmode\nobreak\ b^{\top}]^{\top}$. For a
vector $r\in\mathbb{R}^{n}$, we denote its $i^{\text{th}}$ component by
$r_{i}$. Let $\mathbb{R}$ be the set of real numbers, and
$0_{p,q}\in\mathbb{R}^{p,q}$ be a matrix of zeros. Let $\mathcal{L}^{T,p\times
q}$ denote the set of all block lower-triangular matrices with the following
structure

$\mathcal{M}=\begin{bmatrix}M^{0,0}&0_{p,q}&\dots&0_{p,q}\\\
M^{1,1}&M^{1,0}&\dots&0_{p,q}\\\ \vdots&\vdots&\ddots&\vdots\\\
M^{T-1,T-1}&M^{T-1,T-2}&\dots&M^{T-1,0}\end{bmatrix},$ (1)

where $M^{i,j}\in\mathbb{R}^{p\times q}$. The block diagonal matrix consisting
of matrices $A_{1},\dots,A_{T}$ is denoted by
$\mathrm{blkdiag}(A_{1},\dots,A_{T})$. The matrix $\mathcal{I}$ denotes the
identity with its dimensions either inferred from the context or indicated by
the subscript, i.e.,
$\mathcal{I}_{n_{\mathrm{x}}}\in\mathbb{R}^{n_{\mathrm{x}}\times
n_{\mathrm{x}}}$. Let $\mathbb{B}_{\infty}^{m}$, be the unit ball defined by
$\mathbb{B}_{\infty}^{m}\vcentcolon=\\{d\in\mathbb{R}^{m}|\leavevmode\nobreak\
\|d\|_{\infty}\leq 1\\}$. For a matrix $M\in\mathbb{R}^{n\times m}$, the
$\infty$-norm is given by
$\|M\|_{\infty}\vcentcolon=\max_{d\in\mathbb{B}_{\infty}^{m}}\|Md\|_{\infty}$.
For two sets $\mathbb{W}_{1},\mathbb{W}_{2}\subseteq\mathbb{R}^{n}$, the
Minkowski sum is defined as
$\mathbb{W}_{1}\oplus\mathbb{W}_{2}\vcentcolon=\\{w_{1}+w_{2}|\leavevmode\nobreak\
w_{1}\in\mathbb{W}_{1},w_{2}\in\mathbb{W}_{2}\\}$. We define
$\mathbb{W}^{k}\vcentcolon=\underbrace{\mathbb{W}\times\dots\times\mathbb{W}}_{\text{$k$
times}}$. For a sequence of vectors
$w_{k}\in\mathbb{W}_{k}\subseteq\mathbb{R}^{m}$ and $k\in\mathbb{N}$, we
define
$\mathbf{w}_{0:k}:=(w_{0},\dots,w_{k})\in\mathbb{W}^{0:k}\vcentcolon=\mathbb{W}_{0}\times\dots\times\mathbb{W}_{k}$.

##  II Problem Formulation

We consider the following robust nonlinear optimal control problem:

$\displaystyle\min_{\pi(\cdot)}\quad$ $\displaystyle
J_{T}(\bar{{x}},\pi(\cdot)),$ (2a) s.t. $\displaystyle
x_{k+1}=f(x_{k},u_{k})+w_{k}\leavevmode\nobreak\ \forall k\in\mathbb{N}_{T},$
(2e) $\displaystyle x_{0}=\bar{{x}},$ $\displaystyle
u_{k}=\pi_{k}(\mathbf{x}_{0:k})\leavevmode\nobreak\ \forall
k\in\mathbb{N}_{T},$
$\displaystyle(x_{k},u_{k})\in\mathbb{P}\leavevmode\nobreak\ \forall
k\in\mathbb{N}_{T+1}\leavevmode\nobreak\ \forall w_{k}\in\mathbb{W}.$

The dynamics are given by (2e), with the state
$x_{k}\in\mathbb{R}^{n_{\mathrm{x}}}$, the input
$u_{k}\in\mathbb{R}^{n_{\mathrm{u}}}$ and the disturbance
$w_{k}\in\mathbb{W}\subseteq\mathbb{R}^{n_{\mathrm{x}}}$ at time step $k$. The
initial condition is given by $\bar{{x}}\in\mathbb{R}^{n_{\mathrm{x}}}$ in
(2e). The control input is obtained by optimizing over general causal policies
$\pi_{k}$ (2e), with
$\pi=(\pi_{0},\dots,\pi_{T}):\mathbb{R}^{(T+1)n_{\mathrm{x}}}\mapsto\mathbb{R}^{(T+1)n_{\mathrm{u}}}$,
and the last input $u_{T}$ is kept in the problem formulation for notational
convenience. We primarily focus on the robust constraint satisfaction (2e)
and, for simplicity, consider a general cost (2a) over the prediction horizon
$T$ which does not depend on $w$.

Problem (2a) is not computationally tractable because of the optimization over
the general feedback policy $\pi(\cdot)$ and the robust constraint
satisfaction required in (2e). Consequently, the goal is to find a feasible,
but potentially sub-optimal, solution to this problem. To this end, we define
a nominal trajectory as

$z_{k+1}=f(z_{k},v_{k})\leavevmode\nobreak\ \forall
k\in\mathbb{N}_{T},\leavevmode\nobreak\ z_{0}=\bar{{x}},$ (3)

and restrict the policy to a causal affine time-varying error feedback

$\pi_{k}(\mathbf{x}_{0:k})=v_{k}+\sum_{j=0}^{k-1}K^{k-1,j}\Delta x_{k-j}$ (4)

with
$\pi_{k}(\mathbf{x}_{0:k}):\mathbb{R}^{(k+1)n_{\mathrm{x}}}\mapsto\mathbb{R}^{n_{\mathrm{u}}}$,
$v_{k}\in\mathbb{R}^{n_{\mathrm{u}}}$, $z_{k}\in\mathbb{R}^{n_{\mathrm{x}}}$,
$K^{i,j}\in\mathbb{R}^{n_{\mathrm{u}}\times n_{\mathrm{x}}}$, and the errors
$\Delta x_{k}\vcentcolon=x_{k}-z_{k}$, $\Delta u_{k}\vcentcolon=u_{k}-v_{k}$.
We also consider the following standard assumptions.

######  Assumption II.1.

The nonlinear dynamics (2e)
$f:\mathbb{R}^{n_{\mathrm{x}}}\times\mathbb{R}^{n_{\mathrm{u}}}\mapsto\mathbb{R}^{n_{\mathrm{x}}}$
are three times continuously differentiable.

######  Assumption II.2.

The constraint set $\mathbb{P}$ (2e) is a compact polytopic set with
$\mathbb{P}\vcentcolon=\\{(x,u)|\leavevmode\nobreak\
c_{i}^{\top}(x,u)+b_{i}\leq 0,\forall i\in\mathbb{N}_{I}\\}$, where
$c_{i}\in\mathbb{R}^{n_{\mathrm{x}}+n_{\mathrm{u}}}$ and $b_{i}\in\mathbb{R}$.

######  Assumption II.3.

The disturbance set is given by

$w_{k}\in\mathbb{W}=\\{Ed|\leavevmode\nobreak\
d\in\mathbb{B}_{\infty}^{n_{\mathrm{w}}}\\}=E\mathbb{B}_{\infty}^{n_{\mathrm{w}}}\subseteq\mathbb{R}^{n_{\mathrm{x}}},$
(5)

with $E\in\mathbb{R}^{n_{\mathrm{x}}\times n_{\mathrm{w}}}$.

##  III Robust Nonlinear Optimal Control via SLS

In this section, we derive the main result of the paper using the steps
depicted in Fig. 1, i.e., we propose a formulation to optimize over affine
policies (4) that provide robust constraint satisfaction for the nonlinear
system (2e). We decompose the nonlinear system equivalently as the sum of a
nominal nonlinear system and an LTV error system that accounts both for the
local linearization error (Section III-A) and the additive disturbance. Using
established SLS tools for LTV systems (Section III-B), we parameterize the
closed-loop response for this LTV system. As a result, we obtain an
optimization problem that jointly optimizes the nominal trajectory (3) and the
error feedback (4), and that guarantees robust constraint satisfaction
(Section III-C).

###  III-A Over-approximation of nonlinear reachable set

The goal of this section is to decompose the uncertain nonlinear system into a
nominal nonlinear system and an LTV error system subject to some disturbance.
The linearization of the dynamics (3) around a nominal state and input $(z,v)$
is characterized by the Jacobian matrices:

$\displaystyle A(z,v)\vcentcolon=\left.\frac{\partial f}{\partial
x}\right|_{(x,u)=(z,v)},\leavevmode\nobreak\
B(z,v)\vcentcolon=\left.\frac{\partial f}{\partial u}\right|_{(x,u)=(z,v)}.$
(6)

Using the Lagrange form of the remainder of the Taylor series expansion, we
obtain

$\displaystyle f(x,u)+w$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:nominal_trajectory}}}{{=}}f(z,v)$
(7) $\displaystyle\leavevmode\nobreak\ +A(z,v)(x-z)+B(z,v)(u-v)$
$\displaystyle\leavevmode\nobreak\ +\underbrace{r(x,u,z,v)+w}_{=\vcentcolon
d},$

with the remainder
$r:\mathbb{R}^{n_{\mathrm{x}}}\times\mathbb{R}^{n_{\mathrm{u}}}\times\mathbb{R}^{n_{\mathrm{x}}}\times\mathbb{R}^{n_{\mathrm{u}}}\mapsto\mathbb{R}^{n_{\mathrm{x}}}$
and both the disturbance $w$ and the remainder $r(x,u,z,v)$ are lumped in the
disturbance $d\vcentcolon=r(x,u,z,v)+w\in\mathbb{R}^{n_{\mathrm{x}}}$.

To bound the remainder, we consider the (symmetric) Hessian
$H_{i}:\mathbb{R}^{n_{\mathrm{x}}+n_{\mathrm{u}}}\mapsto\mathbb{R}^{(n_{\mathrm{x}}+n_{\mathrm{u}})\times(n_{\mathrm{x}}+n_{\mathrm{u}})}$
of the $i^{\text{th}}$ component of $f$, i.e.,

$H_{i}(\xi)=\left.\begin{bmatrix}\frac{\partial^{2}f_{i}}{\partial
x^{2}}&\frac{\partial^{2}f_{i}}{\partial x\partial u}\\\
*&\frac{\partial^{2}f_{i}}{\partial
u^{2}}\end{bmatrix}\right|_{\left(x,u\right)=\xi},$ (8)

where $\xi\in\mathbb{R}^{n_{\mathrm{x}}+n_{\mathrm{u}}}$ lies between the
linearization point $(z,v)$ and the evaluation point $(x,u)$. We define the
constant111Note that $\max_{\|h\|_{\infty}\leq
1}|h^{\top}Hh|\leq\sum_{i}\sum_{j}|H_{ij}|$, with $H_{ij}$, the element on the
$i^{\text{th}}$ row and $j^{\text{th}}$ column of the matrix $H$.
$\mu\in\mathbb{R}^{n_{\mathrm{x}}\times n_{\mathrm{x}}}$ as

$\mu\vcentcolon=\text{diag}(\mu_{1},\dots,\mu_{n_{\mathrm{x}}}),\leavevmode\nobreak\
\mu_{i}\vcentcolon=\frac{1}{2}\max_{\xi\in\mathbb{P},\|h\|_{\infty}\leq
1}|h^{\top}H_{i}(\xi)h|,$ (9)

and the error $e_{k}\vcentcolon=\left(\Delta x_{k},\Delta
u_{k}\right)\in\mathbb{R}^{n_{\mathrm{x}}+n_{\mathrm{u}}}$.

######  Proposition III.1.

Given Assumptions II.1 and II.2, the remainder in (7) satisfies

$|r_{i}(x,u,z,v)|\leq\|e\|_{\infty}^{2}\mu_{i},$ (10)

for any $(x,u)\in\mathbb{P},(z,v)\in\mathbb{P}$.

###### Proof.

Applying the definition of the Lagrange form of the remainder, for all
$(x,u,z,v)$ and for all $i$, there exists a $\xi\in\mathbb{P}$ (using
convexity) such that

$\displaystyle|r_{i}(x,u,z,v)|$
$\displaystyle=\frac{1}{2}|e^{\top}H_{i}(\xi)e|$
$\displaystyle\leq\max_{\xi\in\mathbb{P}}\frac{1}{2}|e^{\top}H_{i}(\xi)e|$
(11)
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:def_mu}}}{{\leq}}\|e\|_{\infty}^{2}\mu_{i}.$

The constants $\mu_{i}$ are computed offline, which is the only offline design
required for the proposed method. Intuitively, the magnitude of $\mu_{i}$
quantifies the nonlinearity of the system, which will be incorporated in the
disturbance propagation in the proposed formulation (Sec. III.C).

Due to Proposition III.1 and Assumption II.3, the combined disturbance $d$
from (7) satisfies

$\displaystyle d\vcentcolon=w+r(x,u,z,v)$ $\displaystyle\in
E\mathbb{B}_{\infty}^{n_{\mathrm{w}}}\oplus\|e\|_{\infty}^{2}\mu\mathbb{B}_{\infty}^{n_{\mathrm{x}}}$
(12)
$\displaystyle=[E,\|e\|_{\infty}^{2}\mu]\mathbb{B}_{\infty}^{n_{\mathrm{w}}+n_{\mathrm{x}}}.$

Given a bound on $\|e\|_{\infty}$, we can compute an outer approximation of
the reachable set of the nonlinear system, using the following LTV error
system

$\Delta x_{k+1}=A_{k}\Delta x_{k}+B_{k}\Delta u_{k}+d_{k}\leavevmode\nobreak\
\forall k\in\mathbb{N}_{T+1},\leavevmode\nobreak\ \Delta
x_{0}=0_{n_{\mathrm{x}}},$ (13)

with $A_{k}\vcentcolon=A(z_{k},v_{k})$, $B_{k}\vcentcolon=B(z_{k},v_{k})$.

Similar LTV error dynamics are used in [8, 9, 12, 10, 14] to over-approximate
the reachable set. To optimize affine error feedback (4) while ensuring robust
constraint satisfaction of the nonlinear system based on this LTV error
dynamics, we next study the robust optimal control problem for the special
case of LTV systems.

###  III-B Robust optimal control for LTV systems

In this section, we study the parameterization of affine feedback policies for
LTV systems, which will provide the basis for the employed parameterization of
affine error feedback for nonlinear systems. We show that by using SLS
techniques [19] we can jointly optimize the error feedback and nominal
trajectory of any LTV system. Consider an LTV system of the form

$\tilde{x}_{k+1}=\tilde{A}_{k}\tilde{x}_{k}+\tilde{B}_{k}\tilde{u}_{k}+\tilde{w}_{k}\leavevmode\nobreak\
\forall k\in\mathbb{N}_{T},\leavevmode\nobreak\
\tilde{x}_{0}=0_{n_{\mathrm{x}}},$ (14)

with $\tilde{A}_{k}\in\mathbb{R}^{n_{\mathrm{x}}\times n_{\mathrm{x}}}$,
$\tilde{B}_{k}\in\mathbb{R}^{n_{\mathrm{x}}\times n_{\mathrm{u}}}$, and

$\tilde{w}_{k}\in\tilde{\mathbb{W}}_{k}\vcentcolon=\tilde{E}_{k}\mathbb{B}_{\infty}^{n_{\tilde{\mathrm{w}}}},$
(15)

with some $\tilde{E}_{k}\in\mathbb{R}^{n_{\mathrm{x}}\times
n_{\tilde{\mathrm{w}}}}$. The dynamics (14) are written compactly as

$\tilde{\mathbf{x}}=\mathcal{Z}\tilde{\mathcal{A}}\tilde{\mathbf{x}}+\mathcal{Z}\tilde{\mathcal{B}}\tilde{\mathbf{u}}+\tilde{\mathbf{w}},$
(16)

with
$\tilde{\mathbf{w}}\vcentcolon=\left(\tilde{w}_{0},\dots,\tilde{w}_{T-1}\right)\in\mathbb{R}^{Tn_{\mathrm{x}}}$,
$\tilde{\mathbf{x}}\vcentcolon=\left(\tilde{x}_{1},\dots,\tilde{x}_{T}\right)\in\mathbb{R}^{Tn_{\mathrm{x}}}$,
$\tilde{\mathbf{u}}\vcentcolon=\left(\tilde{u}_{1},\dots,\tilde{u}_{T}\right)\in\mathbb{R}^{Tn_{\mathrm{u}}}$,
$\tilde{\mathcal{A}}\vcentcolon=\text{blkdiag}(\tilde{A}_{1},\dots,\tilde{A}_{T-1},0_{n_{\mathrm{x}},n_{\mathrm{x}}})\in\mathcal{L}^{T,n_{\mathrm{x}}\times
n_{\mathrm{x}}}$,
$\tilde{\mathcal{B}}\vcentcolon=\text{blkdiag}(\tilde{B}_{1},\dots,\tilde{B}_{T-1},0_{n_{\mathrm{x}},n_{\mathrm{u}}})\in\mathcal{L}^{T,n_{\mathrm{x}}\times
n_{\mathrm{u}}}$ and the block-lower shift matrix
$\mathcal{Z}\in\mathcal{L}^{T,n_{\mathrm{x}}\times n_{\mathrm{x}}}$ is given
by

$\mathcal{Z}\vcentcolon=\begin{bmatrix}0_{n_{\mathrm{x}},n_{\mathrm{x}}}&0_{n_{\mathrm{x}},n_{\mathrm{x}}}&\dots&0_{n_{\mathrm{x}},n_{\mathrm{x}}}\\\
\mathcal{I}_{n_{\mathrm{x}}}&0_{n_{\mathrm{x}},n_{\mathrm{x}}}&\dots&0_{n_{\mathrm{x}},n_{\mathrm{x}}}\\\
\vdots&\ddots&\ddots&\vdots\\\
0_{n_{\mathrm{x}},n_{\mathrm{x}}}&\dots&\mathcal{I}_{n_{\mathrm{x}}}&0_{n_{\mathrm{x}},n_{\mathrm{x}}}\end{bmatrix}.$
(17)

We introduce the causal linear feedback
$\tilde{\mathbf{u}}=\tilde{\mathcal{K}}\tilde{\mathbf{x}}$,
$\tilde{\mathcal{K}}\in\mathcal{L}^{T,n_{\mathrm{u}}\times n_{\mathrm{x}}}$
i.e., $\tilde{u}_{k}=\sum_{j=0}^{k-1}\tilde{K}^{k-1,j}\tilde{x}_{k-j}$,
$\tilde{K}^{i,j}\in\mathbb{R}^{n_{\mathrm{u}}\times n_{\mathrm{x}}}$. Using
this feedback, we can write the closed-loop dynamics as

$\tilde{\mathbf{x}}=\mathcal{Z}(\tilde{\mathcal{A}}+\tilde{\mathcal{B}}\tilde{\mathcal{K}})\tilde{\mathbf{x}}+\tilde{\mathbf{w}},\leavevmode\nobreak\
\tilde{\mathbf{u}}=\tilde{\mathcal{K}}\tilde{\mathbf{x}},\leavevmode\nobreak\
\tilde{x}_{0}=0_{n_{\mathrm{x}}},$ (18)

or equivalently as

$\begin{bmatrix}\tilde{\mathbf{x}}\\\
\tilde{\mathbf{u}}\end{bmatrix}=\begin{bmatrix}(\mathcal{I}-\mathcal{Z}\tilde{\mathcal{A}}-\mathcal{Z}\tilde{\mathcal{B}}\tilde{\mathcal{K}})^{-1}\\\
\tilde{\mathcal{K}}(\mathcal{I}-\mathcal{Z}\tilde{\mathcal{A}}-\mathcal{Z}\tilde{\mathcal{B}}\tilde{\mathcal{K}})^{-1}\end{bmatrix}\tilde{\mathbf{w}}=\vcentcolon\begin{bmatrix}\tilde{\Phi}_{\mathrm{x}}\\\
\tilde{\Phi}_{\mathrm{u}}\end{bmatrix}\tilde{\mathbf{w}},\\\ $ (19)

with

$\displaystyle\tilde{\Phi}_{\mathrm{x}}$
$\displaystyle=\begin{bmatrix}\tilde{\Phi}_{\mathrm{x}}^{0,0}&&\\\
\vdots&\ddots&\\\
\tilde{\Phi}_{\mathrm{x}}^{T-1,T-1}&\cdots&\tilde{\Phi}_{\mathrm{x}}^{T-1,0}\\\
\end{bmatrix}\in\mathcal{L}^{T,n_{\mathrm{x}}\times n_{\mathrm{x}}},$ (20)
$\displaystyle\tilde{\Phi}_{\mathrm{u}}$
$\displaystyle=\begin{bmatrix}\tilde{\Phi}_{\mathrm{u}}^{0,0}&&\\\
\vdots&\ddots&\\\
\tilde{\Phi}_{\mathrm{u}}^{T-1,T-1}&\cdots&\tilde{\Phi}_{\mathrm{u}}^{T-1,0}\\\
\end{bmatrix}\in\mathcal{L}^{T,n_{\mathrm{u}}\times n_{\mathrm{x}}}.$

The matrices $\tilde{\Phi}_{\mathrm{x}}$ and $\tilde{\Phi}_{\mathrm{u}}$ are
called the system responses from the disturbance to the closed-loop state and
input, respectively. The following proposition from the literature shows that
the closed-loop response under arbitrary affine feedback is given by all
system responses in a linear subspace.

######  Proposition III.2.

[19, adapted from Theorem 2.1] Let
$\tilde{\mathbf{w}}\in\tilde{\mathbb{W}}^{0:T-1}$ be an arbitrary disturbance
sequence. Any $\tilde{\mathbf{x}}$, $\tilde{\mathbf{u}}$ satisfying (18), also
satisfy (19) with some $\tilde{\Phi}_{x}\in\mathcal{L}^{T,n_{\mathrm{x}}\times
n_{\mathrm{x}}}$, $\tilde{\Phi}_{u}\in\mathcal{L}^{T,n_{\mathrm{u}}\times
n_{\mathrm{x}}}$ lying on the affine subspace

$\left.\begin{aligned}
&\left[\mathcal{I}-\mathcal{Z}\tilde{\mathcal{A}}\quad-\mathcal{Z}\tilde{\mathcal{B}}\right]\left[\begin{array}[]{c}\tilde{\Phi}_{\mathrm{x}}\\\
\tilde{\Phi}_{\mathrm{u}}\end{array}\right]=\mathcal{I}.\end{aligned}\right.$
(21)

Let $\tilde{\Phi}_{\mathrm{x}}$ and $\tilde{\Phi}_{\mathrm{u}}$ be arbitrary
matrices satisfying (21). Then the corresponding $\tilde{\mathbf{x}}$ and
$\tilde{\mathbf{u}}$ computed with (19) also satisfy (18) with
$\tilde{\mathcal{K}}=\tilde{\Phi}_{\mathrm{u}}\tilde{\Phi}_{\mathrm{x}}^{-1}\in\mathcal{L}^{T,n_{\mathrm{u}}\times
n_{\mathrm{x}}}$.

###### Proof.

The proof follows directly from [19, Theorem 2.1] for systems with zero
initial conditions. ∎

Note that this proposition holds for any LTV system and in particular for
(13), where the matrices $\tilde{\mathcal{A}}$ and $\tilde{\mathcal{B}}$
depend on the nominal trajectory $(z,v)$.

Next, we show how the parameterization (21) can be utilized to exactly solve
Problem (2a) with affine feedback in the case of LTV systems. To this end, we
introduce a nominal LTV system

$\tilde{z}_{k+1}=\tilde{A}_{k}\tilde{z}_{k}+\tilde{B}_{k}\tilde{v}_{k}\leavevmode\nobreak\
\forall k\in\mathbb{N}_{T},\leavevmode\nobreak\ \tilde{z}_{0}=\tilde{x}_{0}.$
(22)

The error dynamics are written as

$\tilde{x}_{k+1}-\tilde{z}_{k+1}=\tilde{A}_{k}(\tilde{x}_{k}-\tilde{z}_{k})+\tilde{B}_{k}(\tilde{u}_{k}-\tilde{v}_{k})+\tilde{w}_{k}\leavevmode\nobreak\
\forall k\in\mathbb{N}_{T}.$ (23)

By considering a causal affine error feedback of the form

$\tilde{\mathbf{u}}=\tilde{\mathbf{v}}+\tilde{\mathcal{K}}(\tilde{\mathbf{x}}-\tilde{\mathbf{z}}),\leavevmode\nobreak\
\tilde{u}_{0}=\tilde{v}_{0},\leavevmode\nobreak\
\tilde{\mathcal{K}}\in\mathcal{L}^{T,n_{\mathrm{u}}\times n_{\mathrm{x}}},$
(24)

we can apply Proposition III.2 to characterize the system response of the LTV
error system (23). Note that the linear feedback on the error system
corresponds to the affine feedback in (24). The closed-loop error on the
states and inputs is expressed using the definition of
$\tilde{\Phi}_{\mathrm{x}}$ and $\tilde{\Phi}_{\mathrm{u}}$ in (19) for the
LTV system (23), i.e.,

$\tilde{e}_{k}=\left(\tilde{x}_{k},\tilde{u}_{k}\right)-\left(\tilde{z}_{k},\tilde{v}_{k}\right)=\sum_{j=0}^{k-1}\tilde{\Phi}^{k-1,j}\tilde{w}_{k-1-j}\leavevmode\nobreak\
\forall k\in\mathbb{N}_{T+1},$ (25)

where
$\tilde{\Phi}^{k,j}\vcentcolon=\left(\tilde{\Phi}_{\mathrm{x}}^{k,j},\tilde{\Phi}_{\mathrm{u}}^{k,j}\right)$,
$\tilde{\Phi}_{\mathrm{x}}\in\mathcal{L}^{T,n_{\mathrm{x}}\times
n_{\mathrm{x}}}$ and
$\tilde{\Phi}_{\mathrm{u}}\in\mathcal{L}^{T,n_{\mathrm{u}}\times
n_{\mathrm{x}}}$. Given Proposition III.2, we can provide tight conditions for
robust state and input constraint satisfaction for the LTV system.

######  Proposition III.3.

There exists a causal error feedback of the form (24) such that for any
$\tilde{\mathbf{w}}\in\tilde{\mathbb{W}}^{T}$

$c_{i}^{\top}\left(\tilde{x}_{k},\tilde{u}_{k}\right)+b_{i}\leq
0\leavevmode\nobreak\ \forall k\in\mathbb{N}_{T+1}\leavevmode\nobreak\ \forall
i\in\mathbb{N}_{I},$ (26)

with $\left(\tilde{x}_{k},\tilde{u}_{k}\right)$ according to (14) if and only
if there exist matrices $\tilde{\Phi}_{\mathrm{x}}$,
$\tilde{\Phi}_{\mathrm{u}}$ and a nominal trajectory satisfying (21), (22),
and

$\displaystyle c_{i}^{\top}\left(\tilde{z}_{k},\tilde{v}_{k}\right)+b_{i}$
(27)
$\displaystyle+\sum_{j=0}^{k-1}\|c_{i}^{\top}\tilde{\Phi}^{k-1,j}\tilde{E}_{k-1-j}\|_{1}\leq
0\leavevmode\nobreak\ \forall k\in\mathbb{N}_{T+1}\leavevmode\nobreak\ \forall
i\in\mathbb{N}_{I}.$

###### Proof.

As per Proposition III.2, the LTV error system can be written with Equation
(25), and we can apply directly [22, Example 8]. Namely,$\leavevmode\nobreak\
\forall k\in\mathbb{N}_{T+1}\leavevmode\nobreak\ \forall i\in\mathbb{N}_{I}$,
we have

$\displaystyle\max_{\tilde{\mathbf{w}}_{0:k-1}\in\tilde{\mathbb{W}}^{0:k-1}}c_{i}^{\top}\left(\tilde{x}_{k},\tilde{u}_{k}\right)+b_{i}$
(28)
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:decomp}}}{{=}}c_{i}^{\top}\left(\tilde{z}_{k},\tilde{v}_{k}\right)+\max_{\tilde{\mathbf{w}}_{0:k-1}\in\tilde{\mathbb{W}}^{0:k-1}}\sum_{j=0}^{k-1}c_{i}^{\top}\tilde{\Phi}^{k-1,j}\tilde{w}_{k-1-j}+b_{i}$
$\displaystyle=c_{i}^{\top}\left(\tilde{z}_{k},\tilde{v}_{k}\right)+\sum_{j=0}^{k-1}\max_{\tilde{w}_{j}\in\tilde{\mathbb{W}}_{k-1-j}}c_{i}^{\top}\tilde{\Phi}^{k-1,j}\tilde{w}_{k-1-j}+b_{i}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq_noise_set_lin}}}{{=}}c_{i}^{\top}\left(\tilde{z}_{k},\tilde{v}_{k}\right)+\sum_{j=0}^{k-1}\left\|c_{i}^{\top}\tilde{\Phi}^{k-1,j}\tilde{E}_{k-1-j}\right\|_{1}+b_{i}.$
∎

######  Remark III.1.

A similar reformulation for ellipsoidal disturbances or more general polytopic
disturbances can be found in [22, Example 7-8]. Likewise, the result can be
naturally extended to time-varying constraints.

We have presented conditions for affine error feedback with guaranteed
constraint satisfaction for a given LTV system using the following three
components in the SLS formulation: the nominal trajectory (22), the affine
error feedback parameterization (21), and a description of the disturbance set
$\tilde{\mathbb{W}}$ (15). In combination, Proposition III.2 characterizes the
system response of the error dynamics, while Proposition III.3 provides tight
bounds on the nominal trajectory and system response to ensure robust
constraint satisfaction. By combining these two results, we can solve the
robust optimal control problem (2a) for LTV dynamics (23), affine error
feedback (24), and disturbances of the form (15) using the following
optimization problem:

$\displaystyle\min_{\tilde{\mathbf{z}},\tilde{v}_{0},\tilde{\mathbf{v}},\tilde{\Phi}}\quad$
$\displaystyle
J_{T}(\bar{{x}},\tilde{\mathbf{z}},\tilde{\mathbf{v}},\tilde{\Phi}),$ (29a)
s.t.
$\displaystyle\left[\mathcal{I}-\mathcal{Z}\tilde{\mathcal{A}}\quad-\mathcal{Z}\tilde{\mathcal{B}}\right]\left[\begin{array}[]{c}\tilde{\Phi}_{\mathrm{x}}\\\
\tilde{\Phi}_{\mathrm{u}}\end{array}\right]=\mathcal{I},$
$\displaystyle\tilde{z}_{k+1}=\tilde{A}_{k}\tilde{z}_{k}+\tilde{B}_{k}\tilde{v}_{k}\leavevmode\nobreak\
\forall k\in\mathbb{N}_{T},z_{0}=\bar{{x}},$
$\displaystyle\sum_{j=0}^{k-1}\|c_{i}^{\top}\tilde{\Phi}^{k-1,j}\tilde{E}_{k-1-j}\|_{1}$
$\displaystyle\leavevmode\nobreak\
+c_{i}^{\top}\left(\tilde{z}_{k},\tilde{v}_{k}\right)+b_{i}\leq
0\leavevmode\nobreak\ \forall k\in\mathbb{N}_{T+1}\leavevmode\nobreak\ \forall
i\in\mathbb{N}_{I},$

where we define
$\tilde{\Phi}\vcentcolon=(\tilde{\Phi}_{\mathrm{x}},\tilde{\Phi}_{\mathrm{u}})$.
Specifically, the solution of (29a) provides a jointly optimized nominal
trajectory $\tilde{\mathbf{z}}$, $\tilde{\mathbf{v}}$ and linear error
feedback
$\tilde{\mathcal{K}}=\tilde{\Phi}_{\mathrm{u}}\tilde{\Phi}_{\mathrm{x}}^{-1}$
which guarantee robust satisfaction of the constraints for the closed-loop
state and input trajectories. Problem (29a) is a reformulation of established
results in the literature (see, e.g., [19, 29]). In the following, we address
the nonlinear problem (29a) by merging the robust optimal control for LTV
systems with the linearization error bounds from Section III-B for nonlinear
systems.

###  III-C Robust nonlinear finite-horizon optimal control problem

In this section, given the parameterization of the affine error feedback for
the LTV system (Section III-B) and the linearization error of the nonlinear
system (Section III-A), we are in a position to introduce the main result of
this paper (cf. Fig. 1). In particular, we will show in Theorem III.1 that the
following NLP provides a feasible solution to the robust optimal control
problem in (2a):

$\displaystyle\min_{\mathbf{z},v_{0},\mathbf{v},\Phi,\bm{\tau}}$
$\displaystyle J_{T}(\bar{{x}},\mathbf{z},\mathbf{v},\Phi),$ (30a) s.t.
$\displaystyle\left[\mathcal{I}-\mathcal{Z}\mathcal{A}(\mathbf{z},\mathbf{v})\quad-\mathcal{Z}\mathcal{B}(\mathbf{z},\mathbf{v})\right]\begin{bmatrix}\Phi_{\mathrm{x}}\\\
\Phi_{\mathrm{u}}\\\ \end{bmatrix}=\mathcal{I},$ (30e) $\displaystyle
z_{k+1}=f(z_{k},v_{k})\leavevmode\nobreak\ \forall
k\in\mathbb{N}_{T},z_{0}=\bar{{x}},$
$\displaystyle\sum_{j=0}^{k-1}\|c_{i}^{\top}\Phi^{k-1,j}[E,\tau_{k-1-j}^{2}\mu]\|_{1}$
$\displaystyle\leavevmode\nobreak\
+c_{i}^{\top}\left(z_{k},v_{k}\right)+b_{i}\leq 0\leavevmode\nobreak\ \forall
k\in\mathbb{N}_{T+1}\leavevmode\nobreak\ \forall i\in\mathbb{N}_{I},$
$\displaystyle\sum_{j=0}^{k-1}\|\Phi^{k-1,j}[E,\tau_{k-1-j}^{2}\mu]\|_{\infty}\leq\tau_{k}\leavevmode\nobreak\
\forall k\in\mathbb{N}_{T},$

where we denote a feasible solution as
$\\{\mathbf{z}^{\diamond},v_{0}^{\diamond},\mathbf{v}^{\diamond},\Phi^{\diamond},\bm{\tau}^{\diamond}\\}$,
with
$\mathbf{z}\vcentcolon=\left(z_{1},\dots,z_{T}\right)\in\mathbb{R}^{Tn_{\mathrm{x}}}$,
$\mathbf{v}\vcentcolon=\left(v_{1},\dots,v_{T}\right)\in\mathbb{R}^{Tn_{\mathrm{u}}}$,
$\mathcal{A}(\mathbf{z},\mathbf{v})\vcentcolon=\text{blkdiag}(A_{1},\dots,A_{T-1},0_{n_{\mathrm{x}},n_{\mathrm{x}}})\in\mathcal{L}^{T,n_{\mathrm{x}}\times
n_{\mathrm{x}}}$,
$\mathcal{B}(\mathbf{z},\mathbf{v})\vcentcolon=\text{blkdiag}(B_{1},\dots,B_{T-1},0_{n_{\mathrm{x}},n_{\mathrm{u}}})\in\mathcal{L}^{T,n_{\mathrm{x}}\times
n_{\mathrm{u}}}$, $\bm{\tau}=(\tau_{0},\dots,\tau_{T-1})\in\mathbb{R}^{T}$,
${\Phi}_{\mathrm{x}}\in\mathcal{L}^{T,n_{\mathrm{x}}\times n_{\mathrm{x}}}$,
${\Phi}_{\mathrm{u}}\in\mathcal{L}^{T,n_{\mathrm{u}}\times n_{\mathrm{x}}}$
and $\mu$ according to (9). The nominal prediction is given by (30e). Equation
(30e) computes the system response for the linearization around the nominal
trajectory $\mathbf{z}$, $\mathbf{v}$ (cf. Proposition III.2). The auxiliary
variable $\tau_{k}$ is introduced to upper bound $\|e_{k}\|_{\infty}$, which
is used to obtain a bound on the linearization error (cf. Proposition III.2),
which depends on all previous $\tau_{j}$, $j=0,\dots,k-1$, giving (30e). Using
Proposition III.3, the constraints are tightened with respect to both the
additive disturbance $w_{k}\in\mathbb{W}$ and the linearization error
$\|e_{k}\|_{\infty}\mu$ combined as in (12), i.e., the additive uncertainties
on the LTV error lie in the set
$E\mathbb{B}_{\infty}^{n_{\mathrm{w}}}\oplus\|e\|_{\infty}^{2}\mu\mathbb{B}_{\infty}^{n_{\mathrm{x}}}$.
As a result, the reachable set of the nonlinear system (2e) at time step $k$,
in closed-loop with the affine error feedback computed as in Theorem III.1
satisfies

$x_{k}\in z_{k}^{\diamond}\bigoplus_{j=0}^{k-1}\Phi_{\mathrm{x}}^{\diamond
k-1,j}[E,{\tau^{\diamond
2}_{k-1-j}}\mu]\mathbb{B}_{\infty}^{n_{\mathrm{x}}+n_{\mathrm{w}}}=\vcentcolon\mathbb{D}_{k}\leavevmode\nobreak\
\forall k\in\mathbb{N}_{T+1}.$ (31)

The following theorem summarizes the properties of the proposed NLP (30a).

######  Theorem III.1.

Given Assumptions II.1, II.2 and II.3, suppose the optimization problem (30a)
is feasible. Then, the affine error feedback
$\mathbf{u}=\mathbf{v}^{\diamond}+\mathcal{K}^{\diamond}(\mathbf{x}-\mathbf{z}^{\diamond})$,
$\mathcal{K}^{\diamond}=\Phi_{\mathrm{u}}^{\diamond}{\Phi_{\mathrm{x}}^{\diamond}}^{-1}$,
$u_{0}=v_{0}^{\diamond}$ provides a feasible solution to Problem (2a), i.e.,
the closed-loop trajectories of system (2e) under this error feedback robustly
satisfy the constraints (2e).

###### Proof.

First, the constraints (30e) ensure that the nominal trajectory satisfies the
dynamics (3). Then, we use a Taylor series approximation with respect to the
nominal trajectory (30e), resulting in the LTV error system (13). We apply
Proposition III.2 to the LTV error system (13) and the constraint (30e)
implies that the closed-loop trajectories of the error system satisfy

$\begin{bmatrix}\mathbf{x}-\mathbf{z}\\\
\mathbf{u}-\mathbf{v}\end{bmatrix}=\begin{bmatrix}\Phi_{\mathrm{x}}\\\
\Phi_{\mathrm{u}}\end{bmatrix}\bm{d},$ (32)

with $\bm{d}\vcentcolon=(d_{0},\dots,d_{T-1})\in\mathbb{R}^{Tn_{\mathrm{x}}}$
given by (12), $\mathbf{x}$, $\mathbf{u}$ satisfying (2e) and $\mathbf{z}$,
$\mathbf{v}$ satisfying (3). In the following, we show by induction that the
auxiliary variables $\bm{\tau}\in\mathbb{R}^{T}$ satisfy

$\|e_{j}\|_{\infty}\leq\tau_{j}\leavevmode\nobreak\ \forall
j\in\mathbb{N}_{T}.$ (33)

Inequality (33) holds for $j=0$ since $\|e_{0}\|_{\infty}=0$ (cf. (30e)) and
$\tau_{0}\geq 0$ (cf. (30e)). Note that, as per Proposition III.1, the
disturbance on the LTV error system $d_{k}$ satisfies (12). Then, assuming
Inequality (33) holds $\forall j\in\mathbb{N}_{k-1}$, we have

$d_{j}\in[E,\|e_{j}\|_{\infty}^{2}\mu]\mathbb{B}_{\infty}^{n_{\mathrm{x}}+n_{\mathrm{w}}}\subseteq[E,\tau_{j}^{2}\mu]\mathbb{B}_{\infty}^{n_{\mathrm{x}}+n_{\mathrm{w}}}\leavevmode\nobreak\
\forall j\in\mathbb{N}_{k-1}.$ (34)

Hence, we obtain

$\displaystyle\|e_{k}\|_{\infty}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:cl_nl}}}{{=}}\left\|\sum_{j=0}^{k-1}\Phi^{k-1,j}d_{k-1-j}\right\|_{\infty}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:new_E}}}{{\leq}}\sum_{j=0}^{k-1}\|\Phi^{k-1,j}[E,\tau_{k-1-j}^{2}\mu]\|_{\infty}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:rob_lin_constraints_b}}}{{\leq}}\tau_{k}.$

Therefore, the constraints (30e) ensure that Inequality (33) holds for any
realization of the disturbance. Finally, the constraint (30e) in combination
with Equation (34) ensures that the constraints are robustly satisfied,
analogously to Proposition III.3, which can be seen by substituting
$\tilde{E}_{k}$ for $[E,\tau_{k}^{2}\mu]$ with
$n_{\mathrm{\tilde{w}}}=n_{\mathrm{x}}+n_{\mathrm{w}}$. ∎

######  Remark III.2.

When the disturbance $w$ is set to zero with
$E=0_{n_{\mathrm{x}},n_{{\mathrm{w}}}}$, we recover a nominal trajectory
optimization problem (cf., e.g., [5]) as the constraints (30e) reduce to a
nominal constraint, with $\bm{\tau}=0_{T}$, independent of $\Phi$. In case the
system is linear ($\mu=0_{n_{\mathrm{x}},n_{\mathrm{x}}}$), we recover the
linear SLS formulation (cf. [20, Eq.(26)]) since $\bm{\tau}$ does not enter
the constraints (30e).

######  Remark III.3.

The considered handling of the nonlinear system is based on a linearization
along an online optimized nominal trajectory and is comparable to [12, 8, 13,
14], where the latter results also provide an affine feedback policy. However,
the over-approximation of the reachable set in [13] is based on the
ellipsoidal propagation in [10], where even the linear case is not tight. In
contrast to [9], both the linearized system $\mathcal{A}$, $\mathcal{B}$ and
the bound on the remainder term $\bm{\tau}$ are adjusted online based on the
jointly optimized nominal trajectory $\mathbf{z}$, $\mathbf{v}$, which avoids
conservativeness.

In the next section, we present an inexact SQP variant to solve the NLP (30a),
with a reduced computational footprint with respect to standard SQP methods or
IPOPT and results in QP sub-problems comparable to linear SLS problems [20].

##  IV Inexact SQP for SLS-based robust nonlinear optimal control

Standard methods to solve the NLP (30a) (nonlinear interior point methods
[28], SQP) use the derivatives of the functions used in the constraints. In
general, for a function in $\mathbb{R}^{n}\mapsto\mathbb{R}$, evaluating its
Jacobian and Hessian has a runtime cost up to respectively $2n$ and $8n$ times
the cost of evaluating the function alone [30]. This section discusses how to
obtain a feasible solution to (30a) via a tailored inexact SQP algorithm to
speed up the computations. The proposed SQP avoids the derivatives of the
constraint (30e), and hence the Hessian of the dynamics, and only uses the
Jacobians of the dynamics. As evaluating the derivatives is often very
expensive computationally when solving an NLP, we can expect a speedup of up
to $\frac{1+2n+8n}{1+2n}$ when the evaluation of the dynamics dominates the
computational complexity. To outline the algorithm and its properties, we
utilize a compact formulation of Problem (30a):

$\displaystyle\min_{\mathbf{y}}\quad$ $\displaystyle
J_{T}(\bar{{x}},\mathbf{y}),$ (36a) s.t. $\displaystyle
g(\mathbf{y})=0,\leavevmode\nobreak\ h(\mathbf{y})\leq 0,\leavevmode\nobreak\
s(\mathbf{y})\leq 0,$ (36b)

where $\mathbf{y}\in\mathbb{R}^{n_{\mathrm{y}}}$ collects $\Phi$,
$\mathbf{z}$, $v_{0}$, $\mathbf{v}$, and $\bm{\tau}$ in a vector, as well as
the auxiliary variables needed to encode the 1- and $\infty$\- norms. Here,
$g(\mathbf{y})=0$ and $h(\mathbf{y})\leq 0$ correspond respectively to the
nonconvex equality constraints (30e)-(30e) and the nonconvex inequality
constraints (30e)-(30e), while $s(\mathbf{y})\leq 0$ corresponds to the
additional linear inequality constraints on the auxiliary variables.

A standard SQP implementation iteratively solves the following QP sub-problems

$\displaystyle\min_{\Delta\mathbf{y}}\quad$
$\displaystyle\left.\frac{\partial}{\partial\mathbf{y}}J_{T}(\bar{{x}},\mathbf{y})\right|_{\mathbf{y}=\hat{\mathbf{y}}}\Delta\mathbf{y}+\frac{1}{2}\Delta\mathbf{y}^{\top}H_{J}\Delta\mathbf{y},$
(37a) s.t. $\displaystyle
g(\hat{\mathbf{y}})+\left.\frac{\partial}{\partial\mathbf{y}}g(\mathbf{y})\right|_{\mathbf{y}=\hat{\mathbf{y}}}\Delta\mathbf{y}=0,$
(37d) $\displaystyle
h(\hat{\mathbf{y}})+\left.\frac{\partial}{\partial\mathbf{y}}h(\mathbf{y})\right|_{\mathbf{y}=\hat{\mathbf{y}}}\Delta\mathbf{y}\leq
0,$ $\displaystyle s(\hat{\mathbf{y}}+\Delta\mathbf{y})\leq 0,$

where we denote the solution at the previous iteration with the symbol
$\hat{\leavevmode\nobreak\ }$ and define
$\Delta\mathbf{y}\vcentcolon=\mathbf{y}-\hat{\mathbf{y}}$. The matrix
$H_{J}\in\mathbb{R}^{n_{\mathrm{y}}\times n_{\mathrm{y}}}$ is an approximation
of the Hessian of the Lagrangian (see, e.g., [31] for practical
approximations).

In the following, we describe an inexact SQP algorithm [32] to improve the
computational efficiency and highlight the numerical similarities with linear
SLS [20, 19]. In particular, the linearization of (30e) is given by

$\displaystyle\left[\mathcal{I}-\mathcal{Z}\mathcal{A}(\hat{\mathbf{z}},\hat{\mathbf{v}})\quad\right.\left.-\mathcal{Z}\mathcal{B}(\hat{\mathbf{z}},\hat{\mathbf{v}})\right]\Phi$
(38) $\displaystyle\leavevmode\nobreak\
+\sum_{i=1}^{T(n_{\mathrm{x}}+n_{\mathrm{u}})}\left.\frac{\partial}{\partial\xi_{i}}\left[-\mathcal{Z}\mathcal{A}\leavevmode\nobreak\
-\mathcal{Z}\mathcal{B}\right]\right|_{(\mathbf{z},\mathbf{v})=(\hat{\mathbf{z}},\hat{\mathbf{v}})}\hat{\Phi}(\xi_{i}-\hat{\xi}_{i})=\mathcal{I},$

where $\xi_{i}$ is the $i^{\text{th}}$ element of the vector
$(\mathbf{z},\mathbf{v})$. The proposed inexact SQP algorithm replaces the
equality constraint (38) by

$\displaystyle\left[\mathcal{I}-\mathcal{Z}\mathcal{A}(\hat{\mathbf{z}},\hat{\mathbf{v}})\quad\right.$
$\displaystyle\left.-\mathcal{Z}\mathcal{B}(\hat{\mathbf{z}},\hat{\mathbf{v}})\right]\Phi=\mathcal{I},$
(39)

within the QP (37a), i.e., the Jacobians with respect to $\xi_{i}$ in the
constraint (38) are set to zero. Thus, the proposed inexact SQP avoids
evaluating second-order derivatives of the dynamics (3), leading to improved
computation times as shown in Section V. The numerical complexity of each QP
sub-problem of the inexact SQP scheme is similar to a linear SLS problem
(29a), as each constraint of the resulting QPs, except (30e), has its analog
in (29a).

Under several mild assumptions, inexact SQP using (39) converges to a
feasible, but suboptimal, solution
$\\{\Phi^{\diamond},\mathbf{z}^{\diamond},v_{0}^{\diamond},\mathbf{v}^{\diamond},\bm{\tau}^{\diamond}\\}$
of the original NLP (30a) (see, e.g., [33] for the proof).

######  Remark IV.1.

To decrease the number of decision variables and improve the numerical
efficiency, one could consider $\mathcal{K}$ to be block-banded, as in [24].
It is important to note that the feedback $\mathcal{K}$ in (24) can be
implemented without explicitly computing the inverse of $\Phi_{\mathrm{x}}$
[19].

######  Remark IV.2.

Both SQP and inexact SQP are commonly combined with a globalization strategy
such as trust-region or line search methods to ensure convergence when the
initial guess is relatively far from a local optimum [31].

##  V Case Study: Satellite Attitude Control

The following example demonstrates the benefits of the proposed approach where
we optimize jointly the error feedback and the nominal trajectory for
nonlinear systems with additive disturbances. Moreover, it demonstrates the
numerical properties of the method. In particular, we study a nonlinear
aerospace problem, the constrained attitude control of a satellite [6].

Figure 2: Robust nonlinear optimal control for the attitude control of a
spacecraft computed with inexact SQP. For each state and input, the reachable
set (see Equation (31)) around the nominal trajectory depicts the online-
computed reachable sets (green) and its LTV approximation (blue). The
reachable sets are designed to remain within the constraints (black), and the
sample trajectory (red) is hence guaranteed to stay therein.

###  V-A System and constraints

We consider the following Euler’s equation of the dynamics

$\dot{z}=\begin{bmatrix}\Omega(\omega)q\\\
I_{S}^{-1}(v-\omega\times(I_{S}\omega))\end{bmatrix},$ (40)

with states $z\vcentcolon=(q,\omega)$, attitude quaternion
$q\in\mathbb{R}^{4}$, angular rotation rate $\omega\in\mathbb{R}^{3}$, input
torque $v\in\mathbb{R}^{n_{\mathrm{u}}}$, $n_{\mathrm{x}}=7$,
$n_{\mathrm{u}}=3$ and $\times$ is the cross product. The symmetric inertia
matrix of the satellite is $I_{S}=\text{diag}(5,2,1)\in\mathbb{R}^{3\times 3}$
and

$\Omega(\omega)\vcentcolon=\frac{1}{2}\begin{bmatrix}0&-\omega_{1}&-\omega_{2}&-\omega_{3}\\\
\omega_{1}&0&\omega_{3}&-\omega_{2}\\\ \omega_{2}&-\omega_{3}&0&\omega_{1}\\\
\omega_{3}&\omega_{2}&-\omega_{1}&0\\\ \end{bmatrix},$ (41)

with $\Omega:\mathbb{R}^{3}\mapsto\mathbb{R}^{4\times 4}$. The dynamics are
discretized using the $4^{\text{th}}$ order Runge–Kutta method, using a time
step of one, which results in a negligible discretization error. We consider a
bounded disturbance applied to the system described by (5) with $E=5\cdot
10^{-3}[0_{3,4}\leavevmode\nobreak\
\mathcal{I}_{3}]^{\top}\in\mathbb{R}^{7\times n_{\mathrm{w}}}$ with
$n_{\mathrm{w}}=3$, which could stem, e.g., from the solar radiation pressure,
the flexible modes of the solar panels, the aerodynamic drag, or any
unmodelled dynamics. Additionally, we consider the constraints
$-0.1\leq\omega_{i}\leq 0.1$ and $-0.1\leq v_{i}\leq 0.1$ for $i=1,2,3$. We
use the nominal cost function
$J_{T}(\bar{{x}},\mathbf{z},\mathbf{v})=\sum_{k=0}^{T-1}\ell(z_{k},v_{k})+\ell_{T}(z_{T})$,
with stage cost

$\ell(z,v)=(z-z_{\text{ref}})^{\top}{Q}(z-z_{\text{ref}})+(v-v_{\text{ref}})^{\top}{R}(v-v_{\text{ref}}),$
(42)

and the terminal cost

$\ell_{T}(z)=(z-z_{\text{ref}})^{\top}{Q}(z-z_{\text{ref}}),$ (43)

with ${Q}=0.7\mathcal{I}\in\mathbb{R}^{7\times 7}$,
${R}=\mathcal{I}\in\mathbb{R}^{3\times 3}$, the reference
$z_{\text{ref}}=(1,0,0,0,0,0,0)^{\top}$, $v_{\text{ref}}=(0,0,0)^{\top}$, and
the horizon is $T=10$. As often seen in numerical optimization, we consider an
additional term to the cost function, i.e., we use
$J_{T}+\alpha\mathbf{y}^{\top}\mathbf{y}$, with $\alpha=10^{-2}$. We
approximate the constant
$\mu\approx\text{diag}(3.699,3.703,3.717,3.635,0.649,4.608,5.635)\in\mathbb{R}^{7\times
7}$ from Equation (9) for the nonlinear dynamics (40) by using a Monte-Carlo
method with $10^{4}$ samples for a total offline computation time of 444
seconds. The initial condition is $\bar{{x}}=(\bar{q},\bar{\omega})$, where
$\bar{q}$ is inferred from the Euler angles $(180,45,45)\cdot\frac{\pi}{180}$
and $\bar{\omega}=(-1,-4.5,4.5)\cdot\frac{\pi}{180}$.

###  V-B Results

Figure 3: Comparison between the proposed robust nonlinear optimal control
(30a) (a), the linear optimal control problem, based on the linearization
around a reference point, applied to the nonlinear system (b), the nominal
trajectory optimization without robustness guarantees (c), and the (reduced-
horizon) open-loop robust nonlinear optimal control (d). The reachable set
(31) around the nominal trajectory (blue) is plotted with a different color
for each time step $k$, together with the trajectory with disturbance (red).

The NLP (30a) is solved using the inexact SQP (Sec. IV) formulated with Casadi
[34] with the QPs solved with Gurobi [35]. In addition, a Gauss-Newton
approximation of the Hessian (37a) is used, with a small regularization term,
i.e., we use $H_{J}+\gamma\mathcal{I}_{n_{\mathrm{y}}}$ instead of $H_{J}$ in
(37a), with $\gamma=10^{-2}$. We assume convergence when the condition
$\|(\Delta\mathbf{y},\Delta\bm{\nu})\|_{\infty}\leq 10^{-6}$ is fulfilled,
where $\Delta\bm{\nu}$ is the decrease in the dual solution [31].

####  V-B1 Nonlinear system with affine error feedback

For the problem considered, Fig. 2 shows the solution of the NLP (30a)
computed with inexact SQP, which ensures that the trajectories of the
nonlinear system remain robustly within the constraints. The shaded areas
correspond to the reachable sets of the nonlinear system (Equation (31))
(green) and its LTV approximation (blue), i.e., $\mu=0$. Illustrative
disturbance sequences have been applied and, as ensured by the proposed
design, the resulting trajectory (red) remains within the reachable sets. The
flexible error feedback parameterization allows the tubes to grow in some
directions and shrink in others to meet the constraints, which illustrates the
flexibility of this method. A phase plot of the states $\omega_{2}$ and
$\omega_{3}$ is also depicted in Fig. 3(a) for comparison with other
approaches.

####  V-B2 Comparison with open-loop counterpart

To highlight the benefits of the proposed method, we solve the NLP (30a), with
$\Phi_{\mathrm{u}}=0_{Tn_{\mathrm{u}},Tn_{\mathrm{x}}}$, resulting in an open-
loop robust formulation, i.e.,
$\mathcal{K}=0_{Tn_{\mathrm{u}},Tn_{\mathrm{x}}}$. For the open-loop case
only, we consider a reduced horizon of $T=6$, as a longer horizon leads to
infeasibilities. Indeed, this open-loop robust method does not apply error
feedback, making the reachable set effectively larger. The open-loop robust
formulation maintains the guarantees of robust constraint satisfaction, as the
reachable set or tube always stays within the constraints as depicted in Fig.
3(d). Because of the large size of the tube, the nominal trajectory is forced
to move away from the constraints, leading to poor performance. Indeed, for
the system considered, the tube cross-section keeps increasing, which limits
both the size of the disturbance and the horizon that can be considered.
Therefore, the affine error feedback is instrumental to ensure robust
constraint satisfaction for large disturbances without acting overly
conservatively.

####  V-B3 Comparison with nominal counterpart

The method is also compared with its nominal counterpart, i.e., where we do
not optimize error feedback and neglect the disturbance
($\Phi_{\mathrm{x}}=0_{Tn_{\mathrm{x}},Tn_{\mathrm{x}}},\Phi_{\mathrm{u}}=0_{Tn_{\mathrm{u}},Tn_{\mathrm{x}}}$)
and hence do not guarantee robust constraint satisfaction. The optimal nominal
trajectory satisfies the constraints, but the trajectory with disturbance
results in significant constraint violations (see Fig. 3(c)). Therefore, it is
crucial to consider the disturbance in the optimal control problem.

####  V-B4 Comparison with linear counterpart

We design a controller based on a linearization of the nonlinear dynamics
(40), i.e., we linearize the dynamics around the reference point
$(z_{\text{ref}},v_{\text{ref}})\in\mathbb{R}^{n_{\mathrm{x}}+n_{\mathrm{u}}}$,
ignoring the nonlinearity ($\mu=0_{n_{\mathrm{x}},n_{\mathrm{x}}}$), and solve
the resulting linear SLS problem. Subsequently, we apply the resulting
controller to the nonlinear dynamics. When the linearization point is far from
the operation point, or when the linear model is not a good approximation of
the nonlinear system dynamics, the resulting controller leads to poor
performance and large constraint violations (see Fig. 3(b)).

####  V-B5 Computation times

To assess the performance of the two SQP variants (cf. Section IV), we compare
them with the well-established NLP solver, IPOPT [28] using just-in-time
compilation (cf. jit [34]) and MUMPS as a linear solver. The problem was
solved on an i9-7940X processor with 32GB of RAM memory. Table I shows the
comparison in the number of iterations required until the convergence
condition is satisfied, the solve time, and the optimal cost for the three
methods. The inexact SQP achieves the fastest convergence time with negligible
difference in minimizer, highlighting the benefit of this SQP variant. Fig. 4
illustrates the convergence rate of both SQP variants. The speedup between
inexact-SQP and SQP largely depends on the computational cost of evaluating
the Hessians of the dynamics $f$. Hence, we expect a larger difference for a
stiff system or other high-order integration methods. We observe linear
convergence for both SQP variants, as predicted by the theory in [32, 31].

Figure 4: Numerical convergence of SQP and inexact SQP in terms of primal-dual step size. | Inexact SQP | SQP | IPOPT
---|---|---|---
Iterations | 18 | 18 | 52
Computation time [s] | 4.45 | 5.79 | 70.88
$J_{T}(\bar{{x}},\mathbf{y}^{\diamond})$ | 12.27 | 12.27 | 12.27
Table I: Numerical comparison between SQP, inexact SQP, and IPOPT for solving
Problem (30a) in terms of the number of iterations to reach convergence, the
total computation time, and the optimal cost achieved. A just in time compiler
was used for both SQP variants but not for IPOPT.

##  VI Conclusion

We proposed a novel approach to solve finite horizon constrained robust
optimal control problems for nonlinear systems subject to additive
disturbances, as commonly encountered in trajectory optimization or MPC. One
of the main novelties lies in the joint optimization of a nominal trajectory
and an affine error feedback policy to compensate for disturbances with robust
constraint satisfaction. We also presented an inexact-SQP variant, which
results in sub-problems comparable to linear SLS and reduces the computation
times compared to state-of-the-art methods (SQP, IPOPT). We showcased the
method for the control of a rigid body in rotation with constraints on the
states and inputs. We showed that a robust formulation is needed for the
constraints to be robustly satisfied, that a linearized model may not be
sufficient, and that the optimized affine error feedback improves the overall
performance.

Considering a receding horizon implementation, as is typical in robust MPC,
including a corresponding recursive feasibility and stability analysis is left
for future work.

## References

  * [1] L. Grüne and J. Pannek, _Nonlinear Model Predictive Control: Theory and Algorithms_. Springer, 2017.
  * [2] S. Yu, C. Maier, H. Chen, and F. Allgöwer, “Tube MPC scheme based on robust control invariant set with application to Lipschitz nonlinear systems,” _Systems and Control Letters_ , vol. 62, no. 2, pp. 194–200, 2013\.
  * [3] S. Singh, A. Majumdar, J. J. Slotine, and M. Pavone, “Robust online motion planning via contraction theory and convex optimization,” _Proceedings \- IEEE International Conference on Robotics and Automation_ , pp. 5883–5890, 2017\.
  * [4] J. Köhler, R. Soloperto, M. A. Müller, and F. Allgöwer, “A Computationally Efficient Robust Model Predictive Control Framework for Uncertain Nonlinear Systems,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 2, pp. 794–801, 2021.
  * [5] D. Malyuta, T. P. Reynolds, M. Szmuk, T. Lew, R. Bonalli, M. Pavone, and B. Açıkmeşe, “Convex Optimization for Trajectory Generation: A Tutorial on Generating Dynamically Feasible Trajectories Reliably and Efficiently,” _IEEE Control Systems Magazine_ , vol. 42, no. 5, pp. 40–113, 2022.
  * [6] D. Malyuta, Y. Yu, P. Elango, and B. Açıkmeşe, “Advances in trajectory optimization for space vehicle control,” _Annual Reviews in Control_ , vol. 52, pp. 282–315, 2021.
  * [7] A. Majumdar and R. Tedrake, “Funnel libraries for real-time robust feedback motion planning,” _International Journal of Robotics Research_ , vol. 36, no. 8, pp. 947–982, 2017.
  * [8] M. Althoff, O. Stursberg, and M. Buss, “Reachability analysis of nonlinear systems with uncertain parameters using conservative linearization,” in _Proc. 47th IEEE Conference on Decision and Control_. IEEE, 2008, pp. 4042–4048.
  * [9] M. Cannon, J. Buerger, B. Kouvaritakis, and S. Raković, “Robust tubes in nonlinear model predictive control,” _IEEE Transactions on Automatic Control_ , vol. 56, no. 8, pp. 1942–1947, 2011.
  * [10] B. Houska, “Robust Optimization of Dynamic Systems,” Ph.D. dissertation, KULeuven, 2011.
  * [11] W. Langson, I. Chryssochoos, S. V. Raković, and D. Q. Mayne, “Robust model predictive control using tubes,” _Automatica_ , vol. 40, no. 1, pp. 125–133, 2004.
  * [12] M. E. Villanueva, R. Quirynen, M. Diehl, B. Chachuat, and B. Houska, “Robust MPC via min–max differential inequalities,” _Automatica_ , vol. 77, pp. 311–321, 2017.
  * [13] F. Messerer and M. Diehl, “An Efficient Algorithm for Tube-based Robust Nonlinear Optimal Control with Optimal Linear Feedback,” in _Proc. 60th IEEE Conference on Decision and Control (CDC)_. IEEE, 2021, pp. 6714–6721.
  * [14] T. Kim, P. Elango, and B. Acikmese, “Joint Synthesis of Trajectory and Controlled Invariant Funnel for Discrete-time Systems with Locally Lipschitz Nonlinearities,” _arXiv preprint arXiv:2209.03535_ , 2022.
  * [15] A. Zanelli, J. Frey, F. Messerer, and M. Diehl, “Zero-order robust nonlinear model predictive control with ellipsoidal uncertainty sets,” in _Proc. 7th IFAC Conference on Nonlinear Model Predictive Control (NMPC 2021)_ , vol. 54, no. 6, 2021, pp. 50–57.
  * [16] M. Neunert, C. De Crousaz, F. Furrer, M. Kamel, F. Farshidian, R. Siegwart, and J. Buchli, “Fast nonlinear Model Predictive Control for unified trajectory optimization and tracking,” in _IEEE Proc. International Conference on Robotics and Automation_ , vol. 2016-June, 2016, pp. 1398–1404.
  * [17] D. Gramlich, C. Scherer, and C. Ebenbauer, “Robust Differential Dynamic Programming,” _arXiv preprint arXiv:2205.12632_ , 2022.
  * [18] W. Li and E. Todorov, “Iterative Linear Quadratic Regulator Design for Nonlinear Biological Movement Systems,” in _Proc. 1st International Conference on Informatics in Control, Automation and Robotics_ , 2011, pp. 222–229.
  * [19] J. Anderson, J. C. Doyle, S. H. Low, and N. Matni, “System level synthesis,” _Annual Reviews in Control_ , vol. 47, pp. 364–393, 2019.
  * [20] J. Sieber, A. Zanelli, S. Bennani, and M. N. Zeilinger, “System Level Disturbance Reachable Sets and their Application to Tube-based MPC,” _European Journal of Control_ , 2022.
  * [21] S. Chen, N. Matni, M. Morari, and V. M. Preciado, “System Level Synthesis-based Robust Model Predictive Control through Convex Inner Approximation,” _arXiv preprint arXiv:2111.05509_ , 2021.
  * [22] P. J. Goulart, E. C. Kerrigan, and J. M. Maciejowski, “Optimization over state feedback policies for robust control with constraints,” _Automatica_ , vol. 42, no. 4, pp. 523–533, 2006.
  * [23] P. J. Goulart, “Affine Feedback Policies for Robust Control with Constraints,” Ph.D. dissertation, University of Cambridge, 2006.
  * [24] J. Skaf and S. P. Boyd, “Design of affine controllers via convex optimization,” _IEEE Transactions on Automatic Control_ , vol. 55, no. 11, pp. 2476–2487, 2010.
  * [25] L. Furieri, C. L. Galimberti, and G. Ferrari-Trecate, “Neural System Level Synthesis: Learning over All Stabilizing Policies for Nonlinear Systems,” _arXiv preprint arXiv:2203.11812_ , 2022.
  * [26] L. Conger, J. Shuang, Li, E. Mazumdar, and S. L. Brunton, “Nonlinear System Level Synthesis for Polynomial Dynamical Systems,” _arXiv preprint arXiv:2205.02187_ , 2022.
  * [27] D. Ho, “A System Level Approach to Discrete-Time Nonlinear Systems,” in _Proc. American Control Conference_ , 2020, pp. 1625–1630.
  * [28] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” _Mathematical programming_ , vol. 106, no. 1, pp. 25–57, 2006\.
  * [29] J. Sieber, S. Bennani, and M. N. Zeilinger, “A System Level Approach to Tube-based Model Predictive Control,” _IEEE Control Systems Letters_ , vol. 6, pp. 776–781, 2021.
  * [30] A. Griewank and A. Walther, _Evaluating Derivatives_ , 2nd ed. Society for Industrial and Applied Mathematics, 2008\.
  * [31] S. Wright, J. Nocedal _et al._ , _Numerical optimization_. Springer Science, 1999.
  * [32] A. Zanelli, “Inexact methods for nonlinear model predictive control: stability, applications, and software,” Ph.D. dissertation, Universität Freiburg, 2021.
  * [33] H. G. Bock, M. Diehl, P. Kühl, E. Kostina, J. P. Schlöder, and L. Wirsching, “Numerical methods for efficient and fast nonlinear model predictive control,” _Lecture Notes in Control and Information Sciences_ , vol. 358, no. 1, pp. 163–179, 2007.
  * [34] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “CasADi -A software framework for nonlinear optimization and optimal control,” _Mathematical Programming Computation_ , vol. 11, no. 1, pp. 1–36, 2019.
  * [35] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2022. [Online]. Available: https://www.gurobi.com

|  Antoine P. Leeman received a Master degree in aerospace engineering from
the University of Liège, Belgium and a Diplôme d’Ingénieur (MSc) from ISAE-
SUPAERO, France, both in 2018. After his research exchange in the Korea
Advanced Institute of Science and Technology (KAIST), South Korea, he was
working as a Young Graduate Trainee at the European Space Agency (ESA), the
Netherlands, within the Guidance, Navigation and Control (GNC) section from
2018 to 2020. Since 2020, he is a Ph.D. candidate at the Institute for Dynamic
Systems and Control (IDSC) at ETH Zürich, Switzerland. His research interests
include model predictive control, and control of robotics and aerospace
systems.
---|---
|  Johannes Köhler received his Master degree in Engineering Cybernetics from
the University of Stuttgart, Germany, in 2017. In 2021, he obtained a Ph.D. in
mechanical engineering, also from the University of Stuttgart, Germany, for
which he received the 2021 European Systems & Control Ph.D. award. He is
currently a postdoctoral researcher at the Institute for Dynamic Systems and
Control (IDSC) at ETH Zürich. His research interests are in the area of model
predictive control and control and estimation for nonlinear uncertain systems.
---|---
|  Andrea Zanelli received a BSc degree in Automation Engineering from
Politecnico di Milano and an MSc degree in Robotics, Systems and Control from
ETH Zurich in 2012 and 2015, respectively. He pursued a Ph.D. at the Systems
Control and Optimization Laboratory at the University of Freiburg, Germany.
Since 2021, he is a postdoctoral researcher at the Institute for Dynamic
Systems and Control at ETH Zurich. His research focuses on the development and
software implementation of efficient numerical methods for embedded
optimization and nonlinear model predictive control with numerical and system
theoretic guarantees.
---|---
|  Samir Benanni has a Ph.D in aerospace engineering from Delft University,
Faculty of Aerospace Engineering. He is specialized in the field of dynamical
systems and robust control theory with application to space flight controls
problem. He is at the European Space Agency (ESA) technology center as a
Guidance Navigation and Control (GNC) Fellow and Senior Advisor to the GNC
division. In the last two decades relevant efforts were made in infusing
robust control technologies enabling challenging ESA missions. Among these
were, LISA Pathfinder, Rosetta and VEGA launch vehicle family and many others.
He is conducting various technology maturation programs in the domain of space
transportation, exploration, science and new space missions. He is currently
working on GNC validation and verification technologies, real-time optimized
and data-driven guidance and control technologies for the next generation ESA
technology demonstrator platforms. He is a member of various professional
organisations such as IEEE, AIAA.
---|---
|  Melanie N. Zeilinger is an Associate Professor at ETH Zurich, Switzerland.
She received the Diploma degree in engineering cybernetics from the University
of Stuttgart, Germany, in 2006, and the Ph.D. degree with honors in electrical
engineering from ETH Zurich, Switzerland, in 2011. From 2011 to 2012 she was a
Postdoctoral Fellow with the Ecole Polytechnique Federale de Lausanne (EPFL),
Switzerland. She was a Marie Curie Fellow and Postdoctoral Researcher with the
Max Planck Institute for Intelligent Systems, Tübingen, Germany until 2015 and
with the Department of Electrical Engineering and Computer Sciences at the
University of California at Berkeley, CA, USA, from 2012 to 2014. From 2018 to
2019 she was a professor at the University of Freiburg, Germany. Her current
research interests include safe learning-based control, as well as distributed
control and optimization, with applications to robotics and human-in-the-loop
control.
---|---