# Once upon a time step: A closed-loop approach to robust MPC design

Anilkumar Parsi  Marcell Bartos  Amber Srivastava  Sebastien Gros  Roy S.
Smith  Fellow, IEEE This work was supported by the Swiss National Science
Foundation under Grant 200021_178890 and NCCR Automation Grant 180545.
Corresponding author Anilkumar Parsi.Anilkumar Parsi, Amber Srivastava, Roy
Smith and Marcell Bartos are with the Automatic Control Laboratory, Swiss
Federal Institute of Technology (ETH Zürich), 8092 Zurich, Switzerland
(e-mail<EMAIL_ADDRESS>mbartos@student.ethz.ch).Sebastien Gros is with the Department of Engineering
Cybernetics, Norwegian University of Science and Technology (NTNU), Trondheim,
Norway<EMAIL_ADDRESS>
###### Abstract

A novel perspective on the design of robust model predictive control (MPC)
methods is presented, whereby closed-loop constraint satisfaction is ensured
using recursive feasibility of the MPC optimization. Necessary and sufficient
conditions are derived for recursive feasibility, based on the effects of
model perturbations and disturbances occurring at one time step. Using these
conditions and Farkas’ lemma, sufficient conditions suitable for design are
formulated. The proposed method is called a closed-loop design, as only the
existence of feasible inputs at the next time step is enforced by design. This
is in contrast to most existing formulations, which compute control policies
that are feasible under the worst-case realizations of all model perturbations
and exogenous disturbances in the MPC prediction horizon. The proposed method
has an online computational complexity similar to nominal MPC methods while
preserving guarantees of constraint satisfaction, recursive feasibility and
stability. Numerical simulations demonstrate the efficacy of our proposed
approach.

Predictive control for linear systems, Robust control, Uncertain systems,
Constrained control

##  1 Introduction

Model predictive control (MPC) is one of the popular control strategies in use
today. Owing to the modeling and computational flexibility, MPC has been
successfully implemented in a diverse range of applications — such as oil
refineries [1], autonomous vehicles [2], robotics [3], and supply chain
management [4]. A distinctly highlighting feature of MPC is its ability to
efficiently handle constraints on the system states and control inputs, while
also guaranteeing the stability of the system [5]. In particular, MPC poses
decision making as an optimization problem, where one optimizes an objective
function subject to the constraints resulting from the system dynamics, and
admissible system states and control inputs. The optimization problem is
solved at each time step to compute an open-loop input sequence over a finite
number of time steps (called the prediction horizon), of which only the first
input is applied in closed loop. Hence MPC produces an optimization-based
closed-loop policy from the current system state to an input to be applied to
the system. A popular approach in practice is to use an approximate linear
dynamical model of the system, as it results in convex optimization problems
suitable for real-time implementation [6]. Such an approximation can be
compensated for, by introducing structured uncertainty and disturbances into
the model.

Robust MPC algorithms explicitly consider the effects of model uncertainties
and disturbances in the controller design so that closed-loop stability and
constraint satisfaction can be guaranteed [7]. In addition to these
properties, an important consideration for MPC is whether the optimization
problem remains feasible in closed-loop operation, which is known as recursive
feasibility [8, 9]. The design of robust MPC thus has three desired
properties: (i) constraint satisfaction, (ii) recursive feasibility, and (iii)
closed-loop stability. In this paper, we highlight that generally (ii) already
implies (i), and use this observation to propose a novel method for designing
robust MPC schemes.

The design for (i) and (ii) is coupled and follows an open-loop approach in
most existing MPC works [10, 11, 12]. This involves two steps. First,
constraint satisfaction is explicitly enforced on the open-loop trajectories
computed by the MPC optimizer, considering the worst-case model perturbations
and disturbances along the prediction horizon of MPC. Second, the state of the
system at the end of the prediction horizon is driven into a precomputed
terminal set, which is robustly invariant under a known terminal controller
[7]. This strategy ensures that the open-loop trajectories do not violate
constraints for an infinite time, thus achieving both recursive feasibility
and constraint satisfaction.

Although such an open-loop design approach works, a large degree of
conservatism can originate from ensuring feasibility over multiple time steps
in the prediction horizon. This is because recursive feasibility of MPC is a
one-step property. That is, recursive feasibility only requires that under the
action of the first control input computed by MPC at the current time step,
the MPC optimization is also feasible at the next time step. However, the
open-loop design for recursive feasibility operates via ensuring feasibility
throughout the prediction horizon, instead of the next time step alone. The
conservatism of open-loop approaches is arguably stemming from this
discrepancy, as depicted in Figure 1.

In Figure 1, it can be seen that tube MPC approaches (such as [13, 14])
construct a state tube which accounts for all possible disturbances
$w_{k},w_{k+1},w_{k+2}$ and model perturbations
$\Delta_{k},\Delta_{k+1},\Delta_{k+2}$ in the prediction horizon of three time
steps, before the open-loop state trajectories reach the terminal set.

In addition to the conservatism, many open-loop approaches have a large
computational complexity of online optimization. This is because in an open-
loop approach, MPC either has to recursively outer-approximate the reachable
sets under worst-case perturbations (as in tube MPC methods), or use scenario
trees which grow exponentially along the horizon [15].

$\mathbf{k}$-th
step$x_{k}$$w_{k},\Delta_{k}$$w_{k+1},\Delta_{k+1}$$x_{k+1}$Constraint
set$x_{k+2}$$w_{k+2},\Delta_{k+2}$Terminal set$x_{k+3}$ Figure 1: Classical
tube MPC methods, where the predicted open-loop trajectories remain feasible
under _all_ model perturbations $\Delta_{i}$ and disturbances $w_{i}$
occurring in the prediction horizon.

This work is motivated by the insight that under mild conditions, recursive
feasibility implies constraint satisfaction. That is, if an MPC problem
remains feasible in closed loop, then the system constraints are also
satisfied. Therefore, one need not explicitly ensure constraint satisfaction
along the prediction horizon under the open-loop MPC policy. Although this
observation is not particularly difficult to verify, it is not exploited in
the controller design in most existing methods.

To the best of our knowledge, this insight is not formally presented in
previous literature, but a few early robust MPC methods implicitly use it for
systems affected by additive disturbances [16, 17, 8]. All these works study
the recursive feasibility of the MPC optimization and do not explicitly prove
constraint satisfaction, as recursive feasibility already implies constraint
satisfaction. The methods in [16] and Chapter 6.5 of [8] propose ways to
robustify nominal MPC controllers, if a robust control invariant set is known
a priori. However, these sets are very difficult to compute [18, 19]. Whereas
[17] studies the problem of ensuring recursive feasibility using set-
invariance, it only provides tools for analysis of a given robust MPC
algorithm. The recent work [20] extends the ideas from [17] for systems with
multiplicative perturbations, but also provides conditions suitable for
analysis, not design. In contrast to these early methods, most later works
explicitly design for constraint satisfaction following the open-loop approach
described earlier [10]-[14].

In this work, we present a robust MPC design which explicitly ensures
recursive feasibility in closed-loop operation, and hence refer to it as
closed-loop robust MPC. The proposed method in this work uses constraint
tightening to design a recursively feasible MPC. In constraint tightening
methods, the system constraints along the prediction horizon are tightened to
account for model uncertainty, such that if a nominal trajectory satisfies the
modified constraints, the closed-loop trajectories remain feasible under model
perturbations and disturbances [21, 22]. Because the online optimization
computes only the nominal trajectories, the computational complexity is
similar to that of nominal MPC schemes.

The closed-loop robust MPC proposed in this work is depicted in Figure 2. In
this approach, the robust MPC design accounts only for the disturbance $w_{k}$
and model perturbation $\Delta_{k}$ occurring at the current ($k$-th) time
step. In particular, the MPC constraints are tightened offline such that, for
any feasible nominal trajectory at a time step $k$, disturbances $w_{k}$ and
perturbations $\Delta_{k}$ result in a feasible nominal trajectory at time
step $k+1$, guaranteeing recursive feasibility. Note that, although constraint
tightening techniques already exist in literature ([21, 22]), the tightenings
therein are computed by accounting for all perturbations and disturbances in
the MPC prediction horizon, instead of those occurring at the first time step
as done in this work.

$\mathbf{k}$-th stepNominal trajectory$x_{k}$Tightened constraintsTerminal
set$\mathbf{k+1}$-th step$x_{k}$$w_{k},\Delta_{k}$$x_{k+1}$ Perturbed nominal
---
trajectory for $k+1$-th step
Figure 2: Closed-loop robust MPC (our proposed approach). Here, MPC only
computes a nominal trajectory online such that tightened constraints are
satisfied. The tightenings are designed offline such that, under perturbation
$\Delta_{k}$ and disturbance $w_{k}$ occurring at the _current_ ($k$-th) time
step, the MPC problem will be feasible at the next ($k+1$-th) time step.

In this work, a constraint tightening technique is developed for systems whose
dynamics can be described by a nominal linear model combined with a time-
varying linear fractional perturbation [23] and bounded exogenous
disturbances. The proposed algorithm has two phases: offline and online. In
the offline phase, a non-convex optimization is solved to compute constraint
tightenings. In the online phase, a nominal MPC problem with tightened
constraints is solved at each time step to compute the control inputs. The
proposed framework is designed to ensure recursive feasibility, constraint
satisfaction and input-to-state stability. Numerical examples demonstrate that
the proposed design results in an order of magnitude reduction in online
computation time and improved performance compared to state-of-the-art
methods.

Our work has three novel contributions. First, we formally present the
conditions under which recursive feasibility of MPC optimization implies
closed-loop constraint satisfaction. Second, necessary and sufficient
conditions for recursive feasibility are presented for systems affected by
linear fractional perturbations and disturbances. Third, recursive feasibility
is formulated as a design constraint on the robust MPC controller, by
analyzing only the effects of perturbations and disturbances occurring at one
time step.

###  1.1 Notation

The sets of real numbers is denoted by $\mathbb{R}$ and the sequence of
integers from $n_{1}$ to $n_{2}$ by $\mathbb{N}_{n_{1}}^{n_{2}}$. For a vector
$b$, $\left\|b\right\|_{2}$ represents the $2-$norm. The $i^{th}$ row of a
matrix $A$ is denoted by $[A]_{i}$. The notation $a_{l|k}$ denotes the value
of $a$ at time step $k+l$ computed at the time step $k$. The identity matrix
of size $n\times n$ is denoted by $I_{n}$. A block of zeros of size $n\times
m$ inside a larger matrix will be denoted by $0_{n,m}$ when the size cannot be
inferred from the other blocks, and $0$ otherwise. The minimum eigenvalue of a
matrix $A$ is denoted as $\rho_{\min}(A)$. The Kronecker product of two
matrices $A$ and $B$ is denoted by $A\otimes B$, and $\text{co}\\{\cdot\\}$
represents the convex hull operator. A continuous function
$\alpha:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ is a $\mathcal{K}$
function if $\alpha(0)=0$, $\alpha(s)>0$ for all $s>0$ and it is strictly
increasing; it is a $\mathcal{K}_{\infty}$ function if it is a $\mathcal{K}$
function and $\alpha(s)\rightarrow\infty$ as $s\rightarrow\infty$. A
continuous function $\beta:\mathbb{R}_{\geq
0}\times\mathbb{N}_{0}^{\infty}\rightarrow\mathbb{R}_{\geq 0}$ is a
$\mathcal{K}\mathcal{L}$ function if $\beta(s,t)$ is a $\mathcal{K}$ function
in $s$ for every $t\geq 0$, it is strictly decreasing in $t$ for every $s>0$
and $\beta(s,t)\rightarrow 0$ when $t\rightarrow\infty$.

##  2 Problem formulation

We consider uncertain linear, time-invariant systems of the form:

$\displaystyle x_{k+1}$ $\displaystyle=Ax_{k}+Bu_{k}+B_{p}p_{k}+B_{w}w_{k},$
(1a) $\displaystyle q_{k}$ $\displaystyle=D_{x}x_{k}+D_{u}u_{k}+D_{w}w_{k},$
(1b) $\displaystyle p_{k}$ $\displaystyle=\Delta_{k}q_{k},$ (1c)

where $x_{k}\in\mathbb{R}^{n_{x}}$ represents the state of the system,
$u_{k}\in\mathbb{R}^{n_{u}}$ represents the control input, and
$w_{k}\in\mathbb{R}^{n_{w}}$ represents an exogenous disturbance acting on the
system’s state. In addition, the model uncertainty is captured using a linear
fractional transformation (see [23]), described by
$p_{k},q_{k}\in\mathbb{R}^{n_{p}}$ and the matrix
$\Delta_{k}\in\mathbb{R}^{n_{p}\times n_{p}}$. The vectors $p_{k},q_{k}$ and
the matrix $\Delta_{k}$ cannot be measured, but $\Delta_{k}$ is known to lie
inside the set

$\displaystyle\mathcal{P}:=\text{co}\\{\Delta^{1},\Delta^{2},\ldots,\Delta^{n_{\Delta}}\\},$
(2)

where $n_{\Delta}$ represents the number of vertices defining the convex hull
$\mathcal{P}$.

The exogenous disturbance $w_{k}$ is assumed to lie within a compact set
containing the origin, defined as

$\displaystyle\mathcal{W}:=\\{w\in\mathbb{R}^{n_{w}}|H_{w}w\leq h_{w}\\},$ (3)

where $H_{w}\in\mathbb{R}^{n_{H_{w}}\times n_{w}}$ and
$h_{w}\in\mathbb{R}^{n_{H_{w}}}$.

Moreover, the states and inputs of the system are required to lie in a compact
polytopic constraint set containing the origin, defined as

$\displaystyle\mathcal{C}:=\\{(x,u)\in\mathbb{R}^{n_{x}}{\times}\mathbb{R}^{n_{u}}\>\>|\>\>Fx+Gu\leq
b\\},$ (4)

where $F\in\mathbb{R}^{n_{c}\times n_{x}},G\in\mathbb{R}^{n_{c}\times n_{u}}$
and $b\in\mathbb{R}^{n_{c}}$. The control task is subject to a stage cost
function of the form $l(x,u)$, whose cumulative sum has to be minimized over a
(possibly) infinite horizon. Directly minimizing this cost is intractable, as
it requires optimization over infinite number of variables and constraints.

Instead, model predictive control (MPC) is used to find suboptimal input
sequences [24]. In this approach, a receding horizon strategy is used where
the control inputs over the next $N$ timesteps (called the prediction horizon)
are optimized while also ensuring that the state after $N$ timesteps reaches a
terminal set $\mathcal{X}_{N}:=\\{x\>|\>Yx\leq z\\}$ which contains the
origin. In this work, robustness is ensured using constraint tightening [21,
22]. That is, the MPC optimization problem only uses a nominal model of the
system, and tightens the constraints (4) to ensure robustness in closed loop.

Thus, an optimization problem of the following form is solved at each time
step $k$ using the available state measurement $x_{k}$,

$\displaystyle\min_{\hat{\mathbf{x}}_{k},\hat{\mathbf{u}}_{k}}\quad
l_{N}(\hat{x}_{N|k})+\displaystyle\sum_{i=0}^{N-1}$ $\displaystyle
l(\hat{x}_{i|k},\hat{u}_{i|k})$ (5a) $\displaystyle\text{s.t.}\quad
A\hat{x}_{i|k}+B\hat{u}_{i|k}$
$\displaystyle=\hat{x}_{i+1|k},\quad\hat{x}_{0|k}{=}x_{k},$ (5b)
$\displaystyle F\hat{x}_{i|k}+G\hat{u}_{i|k}$ $\displaystyle\leq b-t_{i},\quad
i\in\mathbb{N}_{0}^{N-1},$ (5c) $\displaystyle Y\hat{x}_{N|k}$
$\displaystyle\leq z-t_{N},$ (5d)

where
$\hat{\mathbf{x}}_{k}{:=}[\hat{x}_{0|k}^{\intercal},\ldots,\hat{x}_{N|k}^{\intercal}]^{\intercal}$
and
$\hat{\mathbf{u}}_{k}{:=}[\hat{u}_{0|k}^{\intercal},\ldots,\hat{u}_{N-1|k}^{\intercal}]^{\intercal}$
represent the nominal state and input trajectories, as predicted by the
nominal dynamics in (5b). In addition, (5d) is a terminal constraint imposed
on the last state in the prediction horizon, and $l_{N}(\cdot)$ is a terminal
cost function. The terms $\\{t_{i}\\}_{i=0}^{N}$ are tightenings imposed on
the constraints along the prediction horizon and the terminal set. The MPC
components $l_{N}(\cdot)$, $Y$, $z$ and $\\{t_{i}\\}_{i=0}^{N}$ are design
choices, and must be computed offline, i.e, before (5) is solved.

Let the optimal solution to (5) be defined by the superscript $(^{*})$ on the
optimization variables. The MPC control law is then defined as
$\pi(x_{k})=\hat{u}_{0|k}^{*}$. That is, only the first control input from
$\hat{\mathbf{u}}_{k}^{*}$ is applied to the system, and (5) is solved again
at the following time step.

In order to achieve desired closed-loop properties, relevant design criteria
must be imposed on the MPC components $l_{N}(\cdot)$, $Y$, $z$ and
$\\{t_{i}\\}_{i=0}^{N}$. The following sections present methods formulating
such design criteria.

##  3 Constraint satisfaction using recursive feasibility

Two main properties that are desired to be satisfied by robust MPC controllers
are (i) constraint satisfaction and (ii) recursive feasibility. Existing
robust MPC methods are designed by explicitly accounting for both (i) and
(ii), for example, in [21]. In this work, a novel design strategy is proposed
such that satisfying (ii) implicitly ensures that (i) is satisfied. For this
purpose, the notion of recursive feasibility (RF) is formally presented [9].

######  Definition 1 (Recursive feasibility).

An MPC controller is recursively feasible if, for any state $x_{k}$ such that
the MPC optimization is feasible, the optimal control input computed by MPC
when applied to the system results in a state $x_{k+1}$ which is also feasible
for the MPC optimization problem.

In the context of the system (1)-(3) and the MPC (5), RF means that for every
$x_{k}$ such that (5) is feasible, the MPC (5) is also feasible for the state
$x_{k+1}$ under the control law $\pi(x_{k})=\hat{u}_{0|k}^{*}$ and for any
realization of $w_{k}\in\mathcal{W}$ and $\Delta_{k}\in\mathcal{P}$.

Note that $\pi(x_{k})$ uses the optimal solution of (5), which cannot be
easily expressed as an explicit function of $x_{k}$. Thus, RF is difficult to
use as a design criterion or even verify [9]. In the light of these
difficulties, it is useful to consider a stronger condition, pertaining to
only the feasibility of MPC, which will ensure that RF is satisfied.

######  Definition 2 (Strong recursive feasibility).

An MPC controller is strongly recursively feasible if, for any state $x_{k}$
such that MPC optimization is feasible, any feasible control input computed by
MPC when applied to the system results in a state $x_{k+1}$ which is also
feasible for the MPC optimization problem.

Note that strong recursive feasibility (SRF) is a sufficient condition for RF,
and has been referred to in literature as _strong recursive feasibility_ [9,
Definition 2.2], or _robust strong feasibility_ [17, Definition 12]. Because
SRF does not depend on the optimality of the control input, it can be used in
the design of MPC controllers. Although using SRF for design could lead to
conservatism, it must be noted that most existing robust MPC methods satisfy
SRF [12, 14].

In the context of the system (1)-(3) and the MPC (5), SRF can be formalized as
follows. Let $\mathcal{U}(x_{k})$ be the set of input sequences such that
(5b)-(5d) are satisfied for a given state $x_{k}$. Additionally, define
$\mathcal{U}_{0}(x_{k})$ as the set of feasible inputs at the first predicted
time step, that is,

$\displaystyle\mathcal{U}_{0}(x_{k}):=\left\\{\hat{u}_{0|k}\>\>|\>\>\exists\\{\hat{u}_{l|k}\\}_{l=0}^{N-1}\in\mathcal{U}(x_{k})\right\\}.$
(6)

The MPC (5) is said to be SRF if, for every state $x_{k}$ such that (5) is
feasible, under any control input $u_{k}\in\mathcal{U}_{0}(x_{k})$, and for
any realization of $w_{k}\in\mathcal{W}$ and $\Delta_{k}\in\mathcal{P}$, the
optimization (5) is also feasible for the state $x_{k+1}$ following the
dynamics (1).

In this paper, we will rely on the fact that, under very mild conditions, RF
implies closed-loop constraint satisfaction. That is, if the MPC scheme is
recursively feasible, then the conditions under which constraints (4) are
always satisfied by the closed-loop trajectories are trivial. These statements
are formalized for the closed loop (1)-(5) in the following theorem.

######  Theorem 1.

Let the optimization (5) be recursively feasible and also be feasible at time
$k=0$. In addition, let $t_{0}\geq 0$. Then, the closed-loop trajectories
resulting from the dynamics (1) and MPC control policy
$\pi(x_{k})=\hat{u}_{0|0}^{*}$ never violate the constraints (4).

###### Proof.

Because the optimization (5) is feasible at time step $k{=}0$ and
$t_{0}{\geq}0$, the applied input $u_{0}=\hat{u}_{0|0}^{*}$ satisfies

$\displaystyle Fx_{0}+Gu_{0}\leq b-t_{0}\leq b.$ (7)

Therefore, (4) is satisfied at time $k{=}0$. Moreover, because RF is
satisfied, the optimization (5) always remains feasible under the MPC control
law $\pi(x)$. By using a similar argument as the one presented above for
$k{=}0$, (4) is satisfied for all $k{\geq}0$. ∎

It can be useful to stress here that Theorem 1 – beyond RF – requires only
that problem (5) is initially feasible, and that the initial constraint (5c)
is “properly” tightened, i.e. that $t_{0}{\geq}0$. The latter requirement is
arguably trivial.

######  Remark 1.

Theorem 1 relies only on the feasibility of the solutions to (5) and not on
their optimality. As a consequence, a straightforward corollary is that
Theorem 1 also holds under SRF. Hence, SRF of the closed loop (1)-(5) will
guarantee constraint satisfaction under any control input
$u_{k}\in\mathcal{U}_{0}(x_{k})$. That is, even if problem (5) cannot be
solved to optimality, constraint satisfaction is guaranteed by SRF.

Most robust MPC methods explicitly build outer bounds on worst-case
trajectories of the system to ensure constraint satisfaction [21, 10], which
can be interpreted as open-loop approaches to robust MPC design. This is
because trajectories the predicted by MPC in open loop are designed to lie
inside the system constraints even with perturbations acting along the
prediction horizon. Theorem 1 enables a design which circumvents this
construction altogether, and is the key element of the closed-loop design
presented in this work.

In the following sections, the MPC components $Y,z,\\{t_{i}\\}_{i=0}^{N}$ and
$l_{N}(\cdot)$ will be designed such that SRF is satisfied. The objective is
formulate the design of these components as an optimization problem that needs
to be solved only once, offline. This offline optimization will be presented
in Section 6.

##  4 Necessary and sufficient conditions for SRF

In order to utilize Theorem 1 for design, the SRF of (5) needs to be
guaranteed. In this section, necessary and sufficient conditions for SRF are
derived. The constraints (5b)-(5d) can be compactly written as

$\displaystyle\mathbf{H}_{xu}\mathbf{S}\mathbf{s}_{k}\leq\mathbf{b}-\mathbf{t},$
(8)

where

$\displaystyle\mathbf{s}_{k}=\Big{[}x_{k}^{\intercal}~{}~{}\hat{\mathbf{u}}_{k}^{\intercal}\Big{]}^{\intercal},~{}\mathbf{H}_{xu}=\begin{bmatrix}I_{N}\otimes
F&0&I_{N}\otimes G\\\ 0&Y&0\end{bmatrix},$
$\displaystyle\mathbf{S}=\left[\begin{array}[]{c|cccc}I&0&...&...&0\\\
A&B&0&...&0\\\ \vdots&\vdots&\ddots&\ddots&\vdots\\\
A^{N-1}&A^{N-2}B&...&B&0\\\ A^{N}&A^{N-1}B&...&AB&B\\\ 0&\lx@intercol\hfil
I_{Nn_{u}}\hfil\lx@intercol\end{array}\right]=\left[\begin{array}[]{c|c}\mathbf{S}_{x}&\mathbf{S}_{u}\end{array}\right]$
(16) $\displaystyle\mathbf{b}=\Big{[}\textbf{1}_{N}^{\intercal}\otimes
b^{\intercal}~{}~{}z^{\intercal}\Big{]}^{\intercal},~{}\text{and
}\mathbf{t}=[t_{0}~{}~{}t_{1}~{}~{}\ldots~{}~{}t_{N}]^{\intercal}.$ (17)

One can readily observe that the vector $\mathbf{s}_{k}$ collects the initial
state $x_{k}$ and input trajectory $\hat{\mathbf{u}}_{k}$. The matrix
$\mathbf{S}$ maps the initial state $x_{k}$ and input trajectory
$\hat{\mathbf{u}}_{k}$, i.e., $\mathbf{s}_{k}$, to the predicted state-input
trajectories; more precisely,
$\mathbf{Ss}_{k}=[\hat{\mathbf{x}}_{k}^{\intercal}~{}~{}\hat{\mathbf{u}}_{k}^{\intercal}]^{\intercal}$.
The matrix $\mathbf{H}_{xu}$ maps the predicted state-input trajectories
($\mathbf{Ss}_{k}$) into the left-hand side of constraints (5c)-(5d). Vectors
$\mathbf{b}$ and $\mathbf{t}$ collect the right-hand side of constraints
(5c)-(5d). Let $\mathbb{P}_{0}$ denote the set of all feasible vectors
$\mathbf{s}_{k}$ satisfying (8), that is,

$\displaystyle\mathbb{P}_{0}$
$\displaystyle:=\\{\,\mathbf{s}\in\mathbb{R}^{n_{x}{+}Nn_{u}}\;|\;\mathbf{H}_{xu}\mathbf{S}\mathbf{s}\leq\mathbf{b}-\mathbf{t}\,\\}.$
(18)

Observe that SRF is satisfied if, for every $\mathbf{s}_{k}\in\mathbb{P}_{0}$,
and for all possible realizations of $x_{k+1}$ under the control input
$u_{k}=\hat{u}_{0|k}$, disturbance $w_{k}\in\mathcal{W}$ and perturbation
$\Delta_{k}\in\mathcal{P}$, there exists an input sequence
$\hat{\mathbf{u}}_{k+1}$ satisfying

$\displaystyle\mathbf{H}_{xu}\mathbf{S}\mathbf{s}_{k+1}\leq\mathbf{b}-\mathbf{t},$
(19)

where
$\mathbf{s}_{k+1}=[x_{k+1}^{\intercal}~{}\hat{\mathbf{u}}_{k+1}^{\intercal}]^{\intercal}$.
Note that (19) can be expressed in terms of $\mathbf{s}_{k}$, $w_{k}$, and
$\Delta_{k}$ by re-writing the state-input trajectories $\mathbf{Ss}_{k+1}$ as

$\displaystyle\mathbf{S}\mathbf{s}_{k+1}$
$\displaystyle=\mathbf{S}_{x}x_{k+1}+\mathbf{S}_{u}\hat{\mathbf{u}}_{k+1},$
$\displaystyle=\mathbf{S}_{x}(Ax_{k}+Bu_{k})+\mathbf{S}_{x}\big{(}B_{w}w_{k}+B_{p}\Delta_{k}(D_{x}x_{k}$
$\displaystyle\quad+D_{u}u_{k}+D_{w}w_{k})\big{)}+\mathbf{S}_{u}\hat{\mathbf{u}}_{k+1}$
$\displaystyle=(\mathbf{L}{+}\mathbf{S}_{x}B_{p}\Delta_{k}D_{xu})\mathbf{s}_{k}$
$\displaystyle\quad+\mathbf{S}_{x}(B_{w}{+}B_{p}\Delta_{k}D_{w})w_{k}+\mathbf{S}_{u}\hat{\mathbf{u}}_{k+1}$
$\displaystyle=\mathbf{C_{s}s}_{k}+\mathbf{C_{w}}w_{k}+\mathbf{S}_{u}\hat{\mathbf{u}}_{k+1},$
(20)

where

$\displaystyle\mathbf{C}_{\mathbf{s}}=\mathbf{L}+\mathbf{S}_{x}B_{p}\Delta_{k}D_{xu},\>\>\mathbf{C}_{\mathbf{w}}=\mathbf{S}_{x}(B_{w}+B_{p}\Delta_{k}D_{w}),$
$\displaystyle\mathbf{L}=\begin{bmatrix}A&B&0_{n_{x},(N-1)n_{u}}\\\
A^{2}&AB&0_{n_{x},(N-1)n_{u}}\\\ \vdots&\vdots&\vdots\\\
A^{N+1}&A^{N}B&0_{n_{x},(N-1)n_{u}}\\\ \lx@intercol\hfil
0_{Nn_{u},n_{x}+Nn_{u}}\hfil\lx@intercol\end{bmatrix},D_{xu}=\begin{bmatrix}D_{x}\\\
D_{u}\\\ 0_{n_{p},(N-1)n_{u}}\end{bmatrix}.$ (21)

Note that the matrix $\mathbf{C}_{\mathbf{s}}$ maps $\mathbf{s}_{k}$ into the
predicted state-input trajectories at time $k+1$. Similarly,
$\mathbf{C}_{\mathbf{w}}$ provides the gain from disturbance $w_{k}$ at time
$k$ to the predicted state-input trajectories at time $k+1$.

Using the above definitions, necessary and sufficient conditions for SRF are
presented in the following theorem.

######  Theorem 2.

The closed loop system (1)-(5) satisfies SRF if, and only if, for all
$\mathbf{s}_{k}{\in}\mathbb{P}_{0}$, $w_{k}{\in}\mathcal{W}$, and
$j{\in}\mathbb{N}_{1}^{n_{\Delta}}$, there exist $\hat{\mathbf{u}}_{k+1}^{j}$
such that

$\displaystyle\mathbf{H}_{xu}\left(\mathbf{C}^{j}_{\mathbf{s}}\mathbf{s}_{k}+\mathbf{C}^{j}_{\mathbf{w}}w_{k}+\mathbf{S_{u}}\hat{\mathbf{u}}_{k+1}^{j}\right)\leq\mathbf{b}-\mathbf{t},$
(22)

where
$\mathbf{C}^{j}_{\mathbf{s}}=\mathbf{L}+\mathbf{S}_{x}B_{p}\Delta^{j}D_{xu}$,
and $\mathbf{C}^{j}_{\mathbf{w}}=\mathbf{S}_{x}(B_{w}+B_{p}\Delta^{j}D_{w})$.

###### Proof.

The idea of the proof is briefly described first. By definition, SRF requires
that for every $\mathbf{s}_{k}{\in}\mathbb{P}_{0}$, $w_{k}{\in}\mathcal{W}$
and $\Delta_{k}{\in}\mathcal{P}$, there exists an input sequence
$\hat{\mathbf{u}}_{k+1}$ satisfying (19). Note that the qualifier for all
$\Delta_{k}\in\mathcal{P}$ can be rewritten as for all
$j{\in}\mathbb{N}_{1}^{n_{\Delta}}$, with $\Delta_{k}{=}\Delta^{j}$. Thus,
existence of independent input sequences $\hat{\mathbf{u}}_{k+1}^{j}$, each
corresponding to $\Delta_{k}{=}\Delta^{j}$, ensures SRF. Using (19) and (4),
this condition can be written as (22). The detailed proof follows.
Sufficiency:
Let there exist inputs $\hat{\mathbf{u}}_{k+1}^{j}$ such that (22) is
satisfied for all $\mathbf{s}_{k}{\in}\mathbb{P}_{0}$,
$w_{k}{\in}\mathcal{W}$, and $j{\in}\mathbb{N}_{1}^{n_{\Delta}}$. Moreover,
let the true perturbation $\Delta_{k}$ be given as

$\displaystyle\Delta_{k}=\sum_{j=1}^{n_{\Delta}}\tau_{j}\Delta^{j},$ (23)

where $\sum_{j=1}^{n_{\Delta}}\tau_{j}=1$ and $\tau_{j}\geq 0$. Such a
representation is possible as $\Delta_{k}$ lies inside $\mathcal{P}$ defined
in (2). Consider the input candidate sequence
$\hat{\mathbf{u}}_{k+1}=\sum_{j=1}^{n_{\Delta}}\tau_{j}\hat{\mathbf{u}}_{k+1}^{j}$
for verifying (19). Then, using (4),

$\displaystyle\mathbf{S}\mathbf{s}_{k+1}$
$\displaystyle=(\mathbf{L}{+}\mathbf{S}_{x}B_{p}\Delta_{k}D_{xu})\mathbf{s}_{k}{+}\mathbf{S}_{u}\hat{\mathbf{u}}_{k{+}1}{+}\mathbf{S}_{x}(B_{w}{+}B_{p}\Delta_{k}D_{w})w_{k}$
$\displaystyle=(\mathbf{L}{+}\mathbf{S}_{x}B_{p}\sum_{j=1}^{n_{\Delta}}\tau_{j}\Delta^{j}D_{xu})\mathbf{s}_{k}{+}\mathbf{S}_{u}\sum_{j=1}^{n_{\Delta}}\tau_{j}\hat{\mathbf{u}}_{k{+}1}^{j}{+}$
$\displaystyle\qquad\qquad\mathbf{S}_{x}(B_{w}{+}B_{p}\sum_{j=1}^{n_{\Delta}}\tau_{j}\Delta^{j}D_{w})w_{k}$
$\displaystyle=\sum_{j=1}^{n_{\Delta}}\tau_{j}\big{(}(\mathbf{L}{+}\mathbf{S}_{x}B_{p}\Delta^{j}D_{xu})\mathbf{s}_{k}{+}\mathbf{S}_{u}\hat{\mathbf{u}}_{k{+}1}^{j}{+}$
$\displaystyle\qquad\qquad\mathbf{S}_{x}(B_{w}{+}B_{p}\Delta^{j}D_{w})w_{k}\big{)}$
$\displaystyle=\sum_{j=1}^{n_{\Delta}}\tau_{j}\left(\mathbf{C}^{j}_{\mathbf{s}}\mathbf{s}_{k}+\mathbf{S_{u}}\hat{\mathbf{u}}_{k+1}^{j}+\mathbf{C}^{j}_{\mathbf{w}}w_{k}\right).$
(24)

Therefore, using (22) and (4),

$\displaystyle\mathbf{H}_{xu}\mathbf{S}\mathbf{s}_{k+1}=\mathbf{H}_{xu}\sum_{j=1}^{n_{\Delta}}\tau_{j}\left(\mathbf{C}^{j}_{\mathbf{s}}\mathbf{s}_{k}+\mathbf{S_{u}}\hat{\mathbf{u}}_{k+1}^{j}+\mathbf{C}^{j}_{\mathbf{w}}w_{k}\right)$
$\displaystyle=\sum_{j=1}^{n_{\Delta}}\tau_{j}\mathbf{H}_{xu}\left(\mathbf{C}^{j}_{\mathbf{s}}\mathbf{s}_{k}+\mathbf{S_{u}}\hat{\mathbf{u}}_{k+1}^{j}+\mathbf{C}^{j}_{\mathbf{w}}w_{k}\right)$
$\displaystyle\leq\sum_{j=1}^{n_{\Delta}}\tau_{j}(\mathbf{b}-\mathbf{t})=\mathbf{b}-\mathbf{t},$
(25)

where the inequality is obtained as (22) is satisfied by the inputs
$\hat{\mathbf{u}}_{k+1}^{j}$. Thus, (19) is satisfied for all
$\Delta_{k}{\in}\mathcal{P}$, $\mathbf{s}_{k}{\in}\mathbb{P}_{0}$, and
$w_{k}{\in}\mathcal{W}$, and the closed-loop satisfies SRF.

Necessity:
Let the closed-loop (1)-(5) satisfy SRF. That is, for all
$\mathbf{s}_{k}{\in}\mathbb{P}_{0}$, $\Delta_{k}{\in}\mathcal{P}$ and
$w_{k}{\in}\mathcal{W}$, there exists $\hat{\mathbf{u}}_{k+1}$ satisfying
(19). Then, consider that for each $j\in\mathbb{N}_{1}^{n_{\Delta}}$,
$\hat{\mathbf{u}}_{k+1}^{j}$ represents the feasible input sequences for (19)
when $\Delta_{k}{=}\Delta^{j}$, and for $\mathbf{s}_{k}{\in}\mathbb{P}_{0}$,
$w_{k}{\in}\mathcal{W}$. Using (19) and (4), one obtains

$\displaystyle\mathbf{H}_{xu}(\mathbf{L}{+}\mathbf{S}_{x}B_{p}\Delta^{j}D_{xu})\mathbf{s}_{k}+\mathbf{S}_{u}\hat{\mathbf{u}}_{k{+}1}^{j}+\quad$
$\displaystyle\mathbf{S}_{x}(B_{w}{+}B_{p}\Delta^{j}D_{w})w_{k}$
$\displaystyle\leq\mathbf{b}-\mathbf{t}$
$\displaystyle\implies\mathbf{H}_{xu}\left(\mathbf{C}^{j}_{\mathbf{s}}\mathbf{s}_{k}+\mathbf{C}^{j}_{\mathbf{w}}w_{k}+\mathbf{S_{u}}\hat{\mathbf{u}}_{k+1}^{j}\right)$
$\displaystyle\leq\mathbf{b}-\mathbf{t},$ (26)

proving that (22) is satisfied for all $j{\in}\mathbb{N}_{1}^{n_{\Delta}}$. ∎

Although Theorem 2 provides necessary and sufficient conditions for the
characterization of SRF, (22) is difficult to directly use in the design of
the MPC components $Y$,$z$ and $\mathbf{t}$. The difficulty arises as one must
guarantee the existence of input sequences satisfying (22) for each
$\mathbf{s}_{k}\in\mathbb{P}_{0}$ and $w_{k}\in\mathcal{W}$.

Theorem 2 explicitly captures the conditions under which SRF holds. Whereas
Theorem 2 formulates these conditions as existence of input sequences, the
criteria can be written as set inclusions using ideas from set-invariance
[25]. The early work [17] presents necessary and sufficient conditions to
achieve SRF when uncertainty occurs as additive disturbances. However, the
resulting conditions are not suited for design. Recently, analogous conditions
were formulated for multiplicative perturbations in [20], where polytope
projections were used to verify SRF, given the MPC parameters.

##  5 Sufficient conditions for SRF using feedback

Theorem 2 is difficult to directly use as a design condition because of the
requirement to guarantee the existence of input sequences. In this section,
such a guarantee is achieved by parameterizing the inputs using feedback
gains.

Note that Theorem 2 allows one to design independent input sequences
$\hat{\mathbf{u}}_{k+1}^{j}$ for each $\Delta^{j}$. This is a powerful
feature, and differs from most existing open-loop strategies, where a single
feedback law is designed to compensate for all perturbations and disturbances
[12]. Exploiting the fact that $n_{\Delta}$ independent input sequences are
sufficient to ensure SRF, for each $j\in\mathbb{N}_{1}^{n_{\Delta}}$, the
input sequences $\hat{\mathbf{u}}_{k+1}^{j}$ are parameterized as

$\displaystyle\hat{u}_{i|k+1}^{j}$
$\displaystyle=\hat{u}_{i+1|k}^{j}+M_{i}^{j}B_{w}w_{k}+K^{j}_{i,\Delta}y_{k},\>\forall
i\in\mathbb{N}_{0}^{N-2},$ $\displaystyle\hat{u}_{N-1|k+1}^{j}$
$\displaystyle=K^{j}\hat{x}_{N|k}+M_{N{-}1}^{j}B_{w}w_{k}+K^{j}_{N{-}1,\Delta}y_{k},$
(27)

where
$y_{k}=\begin{bmatrix}x_{k}^{\intercal}&u_{k}^{\intercal}\end{bmatrix}^{\intercal}$,
and for $j\in\mathbb{N}_{1}^{n_{\Delta}}$ and $i\in\mathbb{N}_{0}^{N-1}$,
$K^{j}$ represents a terminal state feedback gain, $M_{i}^{j}$ represents a
disturbance feedback gain and $K_{i,\Delta}^{j}$ represents a perturbation
feedback gain. The use of disturbance feedback gains is common in robust MPC
[26]. Moreover, each $K_{i,\Delta}^{j}$ compensates for the effect of
$\Delta^{j}$ on the nominal trajectory. The term $y_{k}$ is used for feedback
as the effect of the perturbation $\Delta_{k}$ on the nominal trajectories
also depends on $y_{k}$ (see (4)).

Condition (22) can now be rewritten as the existence of suitable feedback
gains, which can be solved for offline. The following notation is introduced
to account for the newly introduced feedback gains. Let for all
$j\in\mathbb{N}_{1}^{n_{\Delta}}$

$\displaystyle A_{K}^{j}$
$\displaystyle:=A+BK^{j},\>\mathbf{M}^{j}:=\begin{bmatrix}{M_{0}^{j}}^{\intercal}&{M_{1}^{j}}^{\intercal}&\ldots&{M_{N{-}1}^{j}}^{\intercal}\end{bmatrix}^{\intercal},$
$\displaystyle\mathbf{K}^{j}_{\Delta}$
$\displaystyle:=\begin{bmatrix}{K^{j}_{0,\Delta}}^{\intercal}&{K^{j}_{1,\Delta}}^{\intercal}&\ldots&{K^{j}_{N-1,\Delta}}^{\intercal}\end{bmatrix}^{\intercal},$
$\displaystyle\mathbf{L}_{K}^{j}$
$\displaystyle:=\begin{bmatrix}A&B&0&...&0\\\ A^{2}&AB&B&...&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ A^{N}&A^{N-1}B&A^{N-2}B&...&B\\\
A_{K}^{j}A^{N}&A_{K}^{j}A^{N-1}B&A_{K}^{j}A^{N-2}B&...&A_{K}^{j}B\\\\[5.69046pt]
0&0&\lx@intercol\hfil\quad\quad <EMAIL_ADDRESS>K^{j}A^{N}&K^{j}A^{N-1}B&K^{j}A^{N-2}B&...&K^{j}B\end{bmatrix},$
$\displaystyle\mathbf{C}_{K}^{j}$
$\displaystyle:=\mathbf{L}_{K}^{j}+\mathbf{S}_{x}B_{p}\Delta^{j}D_{xu}+$
$\displaystyle\qquad\quad\mathbf{S}_{u}\mathbf{K}^{j}_{\Delta}\begin{bmatrix}I_{n_{x}+n_{u}}&0_{n_{x}+n_{u},n_{u}(N-1)}\end{bmatrix},$
$\displaystyle\mathbf{C}_{M}^{j}$
$\displaystyle:=\mathbf{S}_{x}(B_{w}+B_{p}\Delta^{j}D_{w})+\mathbf{S}_{u}\mathbf{M}^{j}B_{w}.$
(28)

Here matrix $\mathbf{L}_{K}^{j}$ maps a pair of initial state and the nominal
input trajectory (that is, $\mathbf{s}_{k}$) into the corresponding nominal
predicted state-input trajectory at the next time step, under the
parameterization (5). Matrix $\mathbf{C}_{K}^{j}$ maps $\mathbf{s}_{k}$ into a
perturbed state-input trajectory at the next time step, when the perturbation
$\Delta_{k}=\Delta^{j}$ and the input is parameterized as (5). Similarly,
$\mathbf{C}_{M}^{j}$ maps the effects of a disturbance $w_{k}$ on the state-
input trajectory at the next time step.

Consider the following polytopes in $\mathbb{R}^{n_{x}+Nn_{u}+n_{w}}$

$\displaystyle\mathbb{P}_{1}$
$\displaystyle:=\left\\{\,\begin{bmatrix}\mathbf{s}\\\
w\end{bmatrix}\>\;\middle|\;\>\begin{bmatrix}\mathbf{H}_{xu}\mathbf{S}&0\\\
0&H_{w}\end{bmatrix}\begin{bmatrix}\mathbf{s}\\\
w\end{bmatrix}{\leq}\begin{bmatrix}\mathbf{b}-\mathbf{t}\\\
h_{w}\end{bmatrix}\,\right\\},$ $\displaystyle\text{and }\mathbb{P}_{2}^{j}$
$\displaystyle:=\left\\{\,\begin{bmatrix}\mathbf{s}\\\
w\end{bmatrix}\;\middle|\;\mathbf{H}_{xu}\begin{bmatrix}\mathbf{C}_{K}^{j}&\mathbf{C}_{M}^{j}\end{bmatrix}\begin{bmatrix}\mathbf{s}\\\
w\end{bmatrix}\leq\mathbf{b}-\mathbf{t}\,\right\\},$ (29)

for $j\in\mathbb{N}_{1}^{n_{\Delta}}$. Polytopes $\mathbb{P}_{1}$ and
$\mathbb{P}_{2}^{j}$ are depicted in Figure 3, and can be interpreted in the
following manner. Recall that $\mathbb{P}_{0}$ is the set of feasible states
and input sequences for MPC, as defined in (18). Then, $\mathbb{P}_{1}$ is the
augmentation of $\mathbb{P}_{0}$ with the set of possible disturbances (3).
Moreover, $\mathbb{P}_{2}^{j}$ is the set of states, inputs and disturbances
acting at time step $k$ such that, if $\Delta_{k}{=}\Delta^{j}$,
$\mathbf{s}_{k+1}$ lies inside $\mathbb{P}_{0}$. That is, $\mathbb{P}_{2}^{j}$
is the pre-image of $\mathbb{P}_{0}$ under the dynamics (1) with
$\Delta_{k}{=}\Delta^{j}$, $w_{k}{\in}\mathcal{W}$ and $u_{k}$ as in (5).

Therefore, it can be seen that, if $\mathbb{P}_{1}\subseteq\mathbb{P}_{2}^{j}$
and $\Delta_{k}{=}\Delta^{j}$, for any $s_{k}\in\mathbb{P}_{0}$ and
$w_{k}\in\mathcal{W}$, $s_{k+1}$ lies inside $\mathbb{P}_{0}$. Using this
approach, polytopic set inclusions can be defined to guarantee SRF as shown in
the following theorem.

$\mathbb{P}_{1}$$s_{k}\in\mathbb{P}_{0}$$w\in\mathcal{W}$
$\pi_{\text{proj}}(\mathbb{P}_{1})=\mathbf{P_{0}}$
---
(a)$\mathbb{P}_{1}$$s_{k+1}\in\mathbb{P}_{0}$$\mathbb{P}_{2}^{j}$ Under
(1),(5),
---
and $\Delta_{k}=\Delta^{j}$
(b) Figure 3: (a) Projection of $\mathbb{P}_{1}$ on the space of $s_{k}$
overlaps with $\mathbb{P}_{0}$ by definition. (b) $\mathbb{P}_{2}^{j}$
contains all $(s_{k},w_{k})$ such that $s_{k+1}\in\mathbb{P}_{0}$ holds under
the dynamics (1), with $\Delta_{k}=\Delta^{j}$, $w_{k}\in\mathcal{W}$ and
input parameterization (5).

######  Theorem 3.

Consider the polytopes $\mathbb{P}_{1}$ and
$\\{\mathbb{P}_{2}^{j}\\}_{j=1}^{n_{\Delta}}$ proposed in (29). Then, the
closed loop (1)-(5) satisfies SRF if there exist feedback gains
$\\{K^{j}\\}_{j=1}^{n_{\Delta}}$,
$\\{\\{M_{i}^{j}\\}_{j=i}^{n_{\Delta}}\\}_{i=0}^{N-1}$ and
$\\{\\{K_{i,\Delta}^{j}\\}_{j=i}^{n_{\Delta}}\\}_{i=0}^{N-1}$ such that

$\displaystyle\mathbb{P}_{1}\subseteq\mathbb{P}_{2}^{j},\quad\forall
j\in\mathbb{N}_{1}^{n_{\Delta}}.$ (30)

###### Proof.

For (5) to be strongly recursively feasible, (19) must be satisfied for every
$\begin{bmatrix}\mathbf{s}_{k}^{\intercal}&w_{k}^{\intercal}\end{bmatrix}^{\intercal}\in\mathbb{P}_{1}$
and $\Delta_{k}\in\mathcal{P}$.

Let $\Delta_{k}$ be defined as in (23) and a candidate control input sequence
be defined as
$\hat{\mathbf{u}}_{k+1}=\sum_{j=1}^{n_{\Delta}}\tau_{j}\hat{\mathbf{u}}_{k+1}^{j}$,
where $\hat{\mathbf{u}}_{k+1}^{j}$ are parameterized according to (5).
Substituting the candidate sequence in (4),

$\displaystyle\mathbf{S}\mathbf{s}_{k+1}=(\mathbf{L}{+}\mathbf{S}_{x}B_{p}\sum_{j=1}^{n_{\Delta}}\tau_{j}\Delta^{j}D_{xu})\mathbf{s}_{k}$
$\displaystyle{+}\mathbf{S}_{x}(B_{w}{+}B_{p}\sum_{j=1}^{n_{\Delta}}\tau_{j}\Delta^{j}D_{w})w_{k}+\mathbf{S}_{u}\sum_{j=1}^{n_{\Delta}}\tau_{j}\hat{\mathbf{u}}_{k+1}^{j}$
(31a)
$\displaystyle=\sum_{j=1}^{n_{\Delta}}\tau_{j}\bigg{(}\bigg{(}\mathbf{L}_{K}^{j}+\mathbf{S}_{x}B_{p}\Delta^{j}D_{xu}$
$\displaystyle\quad+\mathbf{S}_{u}\mathbf{K}^{j}_{\Delta}\begin{bmatrix}I_{n_{x}+n_{u}}&0_{n_{x}+n_{u},n_{u}(N-1)}\end{bmatrix}\bigg{)}\mathbf{s}_{k}$
$\displaystyle\quad+\textbf{S}_{x}\big{(}B_{w}+B_{p}\Delta^{j}D_{w}\big{)}w_{k}+\mathbf{S}_{u}\mathbf{M}^{j}B_{w}w_{k}\bigg{)}$
(31b)
$\displaystyle=\sum_{j=1}^{n_{\Delta}}\tau_{j}\bigg{(}\mathbf{C}_{K}^{j}\mathbf{s}_{k}+\mathbf{C}_{M}^{j}w_{k}\bigg{)}.$
(31c)

In (31), (31b) is obtained by substituting the control law (5) in (31a) and
using the definition of $\mathbf{L}_{K}^{j}$. Using (31), it can be seen that

$\displaystyle
H_{xu}\mathbf{S}\mathbf{s}_{k+1}=H_{xu}\sum_{j=1}^{n_{\Delta}}\tau_{j}\bigg{(}\mathbf{C}_{K}^{j}\mathbf{s}_{k}+\mathbf{C}_{M}^{j}w_{k}\bigg{)}$
(32a)
$\displaystyle=\sum_{j=1}^{n_{\Delta}}\tau_{j}H_{xu}\big{(}\mathbf{C}_{K}^{j}\mathbf{s}_{k}+\mathbf{C}_{M}^{j}w_{k}\big{)}$
(32b)
$\displaystyle\leq\sum_{j=1}^{n_{\Delta}}\tau_{j}(\mathbf{b}-\mathbf{t})\leq\mathbf{b}-\mathbf{t},$
(32c)

where (32c) is obtained because any
$\begin{bmatrix}\mathbf{s}_{k}^{\intercal}&w_{k}^{\intercal}\end{bmatrix}^{\intercal}$
inside $\mathbb{P}_{1}$ also lies inside $\mathbb{P}_{2}^{j}$ from (30).
Therefore, (19) is satisfied for all feasible realizations of
$\mathbf{s}_{k}$, $\Delta_{k}$ and $w_{k}$ under the proposed feedback if (30)
holds, and so the closed loop (1)-(5) satisfies SRF. ∎

Observe that unlike (22), (30) does not depend on the state $x_{k}$. Thus,
(30) can be used to compute the MPC components offline. The following result
based on Farkas’ Lemma will be used to reformulate (30).

######  Lemma 1.

[7, Lemma 5.6] Let $\mathbb{X}_{i}=\\{\,x\;|\;H_{i}x\leq h_{i}\,\\}$ for
$i=1,2$ be two polytopes. Then, $\mathbb{X}_{1}\subseteq\mathbb{X}_{2}$ is
satisfied if, and only if, there exists a matrix $\Lambda$ such that

$\displaystyle\Lambda\geq 0,\quad\Lambda H_{1}=H_{2},\quad\Lambda h_{1}\leq
h_{2}.$ (33)

Note that the first inequality in (33) is element-wise.

######  Proposition 1.

The set inclusions (30) are satisfied if, and only if, there exist matrices
$\\{\Lambda^{j}\\}_{j=1}^{n_{\Delta}}$ such that

$\displaystyle\Lambda^{j}\geq 0,\quad\Lambda^{j}$
$\displaystyle\begin{bmatrix}\mathbf{b}-\mathbf{t}\\\
h_{w}\end{bmatrix}\leq\mathbf{b}-\mathbf{t},$ (34a)
$\displaystyle\Lambda^{j}\begin{bmatrix}\mathbf{H}_{xu}\mathbf{S}&0\\\
0&H_{w}\end{bmatrix}$
$\displaystyle=\mathbf{H}_{xu}\begin{bmatrix}\mathbf{C}_{K}^{j}&\mathbf{C}_{M}^{j}\end{bmatrix},\>\>\forall
j\in\mathbb{N}_{1}^{n_{\Delta}}.$ (34b)

###### Proof.

The proof follows from a direct application of Lemma 1 to each set inclusion
in (30). ∎

Conditions (34) will be used in the following section to build a robust MPC
scheme satisfying SRF.

##  6 Design of robust MPC satisfying SRF

In this section, (34) will be used as constraints in the design of the MPC
components to ensure SRF and therefore constraint satisfaction. We present
viable algorithmic solutions to enforce (34) while minimizing the tightening
of the constraints, so as to promote a larger region of attraction for the
closed-loop problem.

###  6.1 Terminal components

It can already be seen that Theorem 1 imposes the constraint $t_{0}\geq 0$.
Additionally, Proposition 1 explicitly formulates constraints under which the
MPC problem satisfies SRF. The MPC components in (5) which are to be computed
are $Y,z,\mathbf{t}$, and the terminal cost $l_{N}(\cdot)$. The design for SRF
does not depend on the choice of the cost function, and hence, $l_{N}(\cdot)$
will be defined later.

Ideally, all the components $Y,z,\mathbf{t}$ would be set as decision
variables in the offline design problem. However, this would result in large
number of bilinear constraints in (34), as $Y$ is a component of
$\mathbf{H}_{xu}$ defined in (4). Thus, in this work, a popular strategy from
the literature is used to choose $Y$ and $z$ a priori [7]. Specifically, they
are chosen to be of the structure of maximal robust positively invariant sets
under a terminal feedback law $K_{Y}$, given as

$\displaystyle Y:=\begin{bmatrix}F+GK_{Y}\\\ (F+GK_{Y})(A+BK_{Y})\\\ \vdots\\\
(F+GK_{Y})(A+BK_{Y})^{k^{\prime}}\end{bmatrix},z:=\begin{bmatrix}b\\\ b\\\
\vdots\\\ b\end{bmatrix},$ (35)

where $K_{Y}$ is a feedback chosen such that $A+BK_{Y}$ is Schur stable, and
$k^{\prime}$ is a positive scalar chosen by the designer depending on the
desired flexibility. Although the components $Y$ and $z$ are fixed, the
tightenings $t_{N}$ allow the optimizer to modify the terminal set.

###  6.2 Offline optimization for constraint tightening

Given an uncertain model (1)-(3) and constraint set (4), the terminal set is
computed as in (35), and the following optimization problem can be solved to
compute the tightenings, $\mathbf{t}$.

$\displaystyle\begin{split}\min_{\begin{subarray}{c}\mathbf{t},\\\
\\{K^{j},\Lambda^{j},\mathbf{M}^{j},\mathbf{K}^{j}_{\Delta}\\}_{j=1}^{n_{\Delta}}\end{subarray}}&l_{t}(\mathbf{t})\\\
\text{s.t.}\qquad\eqref{eq:LambdaRefo},\quad&t_{0}\geq 0,\quad
f_{t}(\mathbf{t})\leq 0.\end{split}$ (36)

In (36), $l_{t}(\cdot)$ represents an objective function, and $f_{t}(\cdot)$
represents additional constraints which can be imposed on the tightenings
based on the desired goal of the control application. Moreover, the constraint
$t_{0}\geq 0$ needs to be explicitly enforced, so that recursive feasibility
of (5) results in robust constraint satisfaction, as shown in Theorem 1. Note
that the guarantees derived in Theorem 3 and Proposition 1 are independent of
$l_{t}(\cdot)$ and $f_{t}(\cdot)$, as long as the constraints (34) are
satisfied.

Here, we describe one possible strategy of choosing $l_{t}(\cdot)$ and
$f_{t}(\cdot)$, where the design objective is to maximize the region of
attraction (ROA) of the closed-loop system. The ROA is defined by the feasible
region of the state for the MPC optimization (5), and thus is a convex
polytope. Because the volume of a generic polytope is difficult to compute
even when the hyperplanes are known a priori [27], maximizing ROA is a
difficult task. Under the chosen structure of the optimization problem (5),
minimization of $\left\|\mathbf{t}\right\|_{2}^{2}$ can be used as an
approximation for maximizing ROA as proposed in [22].

However, solely approximating ROA maximization with minimization of
$\left\|\mathbf{t}\right\|_{2}^{2}$ can sometimes lead to poor results. In
this work, the approximation is improved by additionally maximizing the size
of an $l_{1}$-norm ball that can fit inside the feasible region of the state
space. Consider the vertices of the unit-$l_{1}$-norm ball in
$\mathbb{R}^{n_{x}}$, represented by $\\{x^{i}\\}_{i=1}^{2n_{x}}$, which lie
on the positive and negative directions of the principal axes of the state
space. For each vertex of the $l_{1}$-norm ball, define the corresponding MPC
optimization variable $\mathbf{s}^{i}_{\alpha}$ as

$\mathbf{s}^{i}_{\alpha}=\begin{bmatrix}\alpha{x^{i}}\\\
\mathbf{u}^{i}\end{bmatrix},\quad\forall i\in\mathbb{N}_{1}^{2n_{x}}$ (37)

where $\alpha$ is a positive scalar defining the size of the $l_{1}$-norm
ball, and $\mathbf{u}^{i}$ represents a feasible input sequence for the
$i^{\text{th}}$ vertex. Then, the offline optimization problem can be written
as

$\displaystyle\min_{\begin{subarray}{c}\mathbf{t},\alpha,\\{\mathbf{u}^{i}\\}_{i=1}^{2n_{x}},\\\
\\{\Lambda^{j},K^{j},\mathbf{M}^{j},\mathbf{K}^{j}_{\Delta}\\}_{j=1}^{n_{\Delta}}\end{subarray}}$
$\displaystyle\left\|\mathbf{t}\right\|^{2}_{2}-\mu\alpha$ (38a) s.t.
$\displaystyle\Lambda^{j}\geq 0,\quad\alpha>0,\quad t_{0}\geq 0,$ (38b)
$\displaystyle\Lambda^{j}\begin{bmatrix}\mathbf{H}_{xu}\mathbf{S}&0\\\
0&H_{w}\end{bmatrix}$
$\displaystyle=\mathbf{H}_{xu}\begin{bmatrix}\mathbf{C}_{K}^{j}&\mathbf{C}_{M}^{j}\end{bmatrix},$
(38c) $\displaystyle\Lambda^{j}\begin{bmatrix}\mathbf{b}-\mathbf{t}\\\
h_{w}\end{bmatrix}$ $\displaystyle\leq\mathbf{b}-\mathbf{t},\quad\forall
j\in\mathbb{N}_{1}^{n_{\Delta}},$ (38d)
$\displaystyle\mathbf{H}_{xu}\mathbf{S}\mathbf{s}^{i}_{\alpha}$
$\displaystyle\leq\mathbf{b-t},\quad\forall i\in\mathbb{N}_{1}^{2n_{x}},$
(38e)

where $\mu$ is a positive weighting factor between minimizing
$\left\|\mathbf{t}\right\|_{2}$ and maximizing $\alpha$. Additionally, (38e)
represent the feasibility constraints for each vertex of the $l_{1}$-norm
ball.

The optimization problem (38) is non-convex, due to the presence of bilinear
constraints (38d). Although this increases the complexity of design, problem
(38) needs to be solved only once, in the offline phase. Note that the number
of bilinearities is greatly reduced because $Y$ and $z$ are chosen beforehand,
and consequently (38c) are affine constraints. Moreover, given a good initial
guess of the tightenings $\mathbf{t}$, feasible solutions to all the other
optimization variables can be computed by solving a convex optimization
problem (by fixing $\mathbf{t}$ in (38)). This provides a simple way to
initialize non-convex optimization techniques which can solve (38).

######  Remark 2 (Initial guess).

It can be difficult to compute a feasible initial guess for the tightenings
$\mathbf{t}$ in (38), as the tightenings are problem dependent. In practice,
one strategy which produced good results was to use an existing constraint
tightening technique, such as [21] or [22] to guess $\mathbf{t}$. Note that
these techniques were designed for systems affected by additive disturbances
alone, and may not always result in feasible initial guesses for $\mathbf{t}$.

###  6.3 Cost functions and online optimization

Whereas the constraints in (5) are designed to ensure recursive feasibility,
the cost function must be chosen to guarantee closed-loop stability. A
quadratic tracking cost formulation is presented here for simplicity, but the
method proposed here can easily be extended to economic cost functions.
Consider the stage cost defined as

$l(x,u)=x^{\intercal}Q_{x}x+u^{\intercal}Q_{u}u,$ (39)

where the matrices $Q_{x},Q_{u}$ are positive definite. The design task is to
choose a terminal cost of the form $l_{N}(x)=x^{\intercal}Q_{N}x$ such that
the closed-loop system is stable. This can be achieved by choosing $Q_{N}$ to
be a positive definite matrix satisfying $\forall
j\in\mathbb{N}_{1}^{n_{\Delta}}$

$(1+\epsilon)\mathbf{C}_{K}^{j\intercal}Q_{s}\mathbf{C}_{K}^{j}-\mathbf{S}^{\intercal}Q_{s}\mathbf{S}\prec-\begin{bmatrix}P_{x}&0\\\
0&P_{u}\end{bmatrix}=-P,$ (40)

where $Q_{s}:=\text{diag}\\{Q_{x}\otimes I_{N},Q_{N},Q_{u}\otimes I_{N}\\}$,
$\epsilon\in\mathbb{R}_{>0}$ is a design parameter and
$P_{x}\in\mathbb{R}^{n_{x}\times n_{x}}$, $P_{u}\in\mathbb{R}^{Nn_{u}\times
Nn_{u}}$ are positive definite matrices. The values of $\epsilon$ and $P$
affect the resulting stability margin as will be shown Section 6.4. For a
given $\epsilon$ and known values of the feedback gains
$\\{K^{j},\mathbf{M}^{j},\mathbf{K}^{j}_{\Delta}\\}_{j=1}^{n_{\Delta}}$,
constraints (40) are linear matrix inequalities in the variables $Q_{N}$ and
$P$. There are multiple strategies which can be used to choose $Q_{N}$ such
that (40) holds, for example, by minimizing the trace of $Q_{N}$ or minimizing
the deviation of $Q_{N}$ from the infinite horizon LQR gain.

The online optimization can then be written as

$\displaystyle\min_{\hat{\mathbf{x}}_{k},\hat{\mathbf{u}}_{k}}\quad\hat{x}_{N|k}^{\intercal}Q_{N}\hat{x}_{N|k}+$
$\displaystyle\displaystyle\sum_{i=0}^{N-1}\left(\hat{x}_{i|k}^{\intercal}Q_{x}\hat{x}_{i|k}+\hat{u}_{i|k}^{\intercal}Q_{u}\hat{u}_{i|k}\right)$
(41a) $\displaystyle\text{s.t.}\quad A\hat{x}_{i|k}+B\hat{u}_{i|k}$
$\displaystyle=\hat{x}_{i|k+1},\quad\hat{x}_{k,0}=x_{k},$ (41b) $\displaystyle
F\hat{x}_{i|k}+G\hat{u}_{i|k}$ $\displaystyle\leq b-t_{i},\quad
i\in\mathbb{N}_{0}^{N-1},$ (41c) $\displaystyle Y\hat{x}_{N|k}$
$\displaystyle\leq z-t_{N}.$ (41d)

The computational complexity of (41) is the same as that of a nominal MPC
problem. This enables a large computational advantage compared to existing
robust MPC methods, such as tube MPC [14].

###  6.4 Algorithm and properties

Using the offline (38) and online (41) optimization problems, a robust MPC
algorithm is formally described in Algorithm 1. It can be seen that the
offline phase involves the computation of $\mathbf{t}$ by solving (38), and
subsequently choosing the terminal cost matrix such that (40) holds. Although
a less conservative approach would be to include (40) as an additional
constraint in (38), this results in a semi-definite program with bilinear
constraints, and is difficult to solve. An alternative approach is to first
compute $Q_{N}$ and
$\\{K^{j},\mathbf{M}^{j},\mathbf{K}^{j}_{\Delta}\\}_{j=1}^{n_{\Delta}}$
satisfying (40), and then solving (38) with the resulting control gains.

Algorithm 1 Closed-loop robust MPC

1:Offline:

2:Choose $N,K_{Y}$, and compute $Y$ and $z$

3:Choose $\mu$ and $\epsilon$

4:Solve (38) to compute $\mathbf{t}$ and control gains

5:Choose $Q_{N}$ satisfying (40)

1:Online: At each time-step $k\geq 0$:

2:Obtain the measurement $x_{k}$

3:Solve (41)

4:Apply $\pi(x_{k})=\hat{u}^{*}_{k,0}$

The properties of the proposed algorithm will now be discussed. For this
purpose, the notion of regional input-to-state stability (ISS) is first
defined.

######  Definition 3 (Regional ISS in $\mathbb{X}$ [28]).

Given a system with dynamics (1)-(3) and a compact set
$\mathbb{X}\subseteq\mathbb{R}^{n_{x}}$ including the origin as an interior
point, the system is said to be regionally ISS in $\mathbb{X}$ with respect to
$w_{k}$ if $\mathbb{X}$ is a robust positively invariant set and if there
exist a $\mathcal{KL}$ function $\beta(\cdot,\cdot)$ and a $\mathcal{K}$
function $\lambda_{1}(\cdot)$ such that, for all $x_{0}\in\mathbb{X}$ and
$k\geq 0$

$\left\|x_{k}\right\|\leq\beta(\left\|x_{0}\right\|,k)+\lambda_{1}(\left\|\mathbf{w}_{[0:k-1]}\right\|),$
(42)

where
$\left\|\mathbf{w}_{[0:k-1]}\right\|=\sup_{i\in\mathbb{N}_{0}^{k-1}}\left\|w_{i}\right\|$.

A sufficient condition to guarantee ISS of a system is the existence of an ISS
Lyapunov function [28], defined below.

######  Definition 4.

Consider the system (1)-(3) and a compact set
$\mathbb{X}\subseteq\mathbb{R}^{n_{x}}$ including the origin as an interior
point. A function $V:\mathbb{R}^{n}\rightarrow\mathbb{R}_{\geq 0}$ is called
an ISS Lyapunov function in $\mathbb{X}$ with respect to $w_{k}$ if there
exist $\mathcal{K}_{\infty}$ functions
$\sigma_{1}(\cdot),\sigma_{2}(\cdot),\sigma_{3}(\cdot)$ and a $\mathcal{K}$
function $\lambda_{2}(\cdot)$ such that $\forall x_{k}\in\mathbb{X}$

$\displaystyle\sigma_{1}(\left\|x_{k}\right\|)\leq V(x_{k})$
$\displaystyle\leq\sigma_{2}(\left\|x_{k}\right\|),$ (43a) $\displaystyle
V(x_{k+1})-V(x_{k})$
$\displaystyle\leq-\sigma_{3}(\left\|x_{k}\right\|)+\lambda_{2}(w_{k}).$ (43b)

######  Theorem 4.

Let the offline optimization problem (38) be feasible, the cost function
satisfy (40), and the online optimization (41) have a feasible solution at
time $k=0$. Moreover, let $\mathbb{P}_{0,x}$ be the projection of
$\mathbb{P}_{0}$ on the first $n_{x}$ elements of $\mathbf{s}_{k}$, and let
the origin be an interior point of $\mathbb{P}_{0,x}$. Then, the closed loop
formed by a system with dynamics (1)-(3) and the robust MPC controller in
Algorithm 1 satisfies the following properties:

  1. (a)

The problem (41) remains feasible for all $k>0$.

  2. (b)

The constraints (4) are satisfied for all $k>0$.

  3. (c)

The closed-loop system is ISS in $\mathbb{P}_{0,x}$ with respect to $w_{k}$.

###### Proof.

In order to prove property (a), note that the constraints in Proposition 1 are
satisfied because the offline optimization (38) is feasible. Therefore, the
set inclusions (30) are satisfied, and Theorem 3 implies that the online MPC
optimization (41) is strongly recursively feasible. Moreover, because (41)
admits a feasible solution at time $k=0$, it remains feasible for all $k>0$.

The property (b) follows from Theorem 1. This is because the offline
optimization (38) explicitly enforces $t_{0}\geq 0$, and SRF implies RF.

The proof of (c) can be seen as follows. Let the optimal MPC cost function be
a candidate ISS Lyapunov function denoted by

$\displaystyle\begin{split}V(x_{k})&=\hat{x}_{N|k}^{*\intercal}Q_{N}\hat{x}_{N|k}^{*}+\displaystyle\sum_{i=0}^{N-1}\left(\hat{x}_{i|k}^{*\intercal}Q_{x}\hat{x}_{i|k}^{*}+\hat{u}_{i|k}^{*\intercal}Q_{u}\hat{u}_{i|k}^{*}\right)\\\
&=\mathbf{s}_{k}^{*\intercal}\mathbf{S}^{\intercal}Q_{s}\mathbf{S}\mathbf{s}_{k}^{*}.\end{split}$
(44)

In the following, it will be shown that $V(x)$ satisfies (43). Consider
$\sigma_{1}(\left\|x\right\|)=\rho_{\min}(Q_{x})\left\|x\right\|^{2}$, which
is guaranteed to be a lower bound of $V(x)$. Moreover, $V(x)$ is a continuous
and piecewise quadratic function in $\mathbb{P}_{0,x}$ as shown in [29].
Therefore, the existence of $\sigma_{2}(\cdot)$ satisfying (43a) is a direct
consequence of [28, Lemma 4].

In order to show (43b), let the perturbation $\Delta_{k}$ be given as (23) and
consider the candidate input sequences in (5). The cost function satisfies

$\displaystyle
V(x_{k+1})-V(x_{k})=\mathbf{s}_{k+1}^{*\intercal}\mathbf{S}^{\intercal}Q_{s}\mathbf{S}\mathbf{s}_{k+1}^{*}-\mathbf{s}_{k}^{*\intercal}\mathbf{S}^{\intercal}Q_{s}\mathbf{S}\mathbf{s}_{k}^{*}$
$\displaystyle\leq\sum_{j=1}^{n_{\Delta}}\tau_{j}\bigg{(}\mathbf{C}_{K}^{j}\mathbf{s}_{k}^{*}+\mathbf{C}_{M}^{j}w_{k}\bigg{)}^{\intercal}Q_{s}\sum_{l=1}^{n_{\Delta}}\tau_{l}\bigg{(}\mathbf{C}_{K}^{l}\mathbf{s}_{k}^{*}+\mathbf{C}_{M}^{l}w_{k}\bigg{)}$
$\displaystyle\qquad\qquad-\mathbf{s}_{k}^{*\intercal}\mathbf{S}^{\intercal}Q_{s}\mathbf{S}\mathbf{s}_{k}^{*}$
(45a)
$\displaystyle\leq\left(1+\frac{1}{\epsilon}\right)\left(w_{k}^{\intercal}\mathbf{C}_{M}(\tau)^{\intercal}Q_{s}\mathbf{C}_{M}(\tau)w_{k}\right)+$
$\displaystyle\quad(1+\epsilon)\left(\mathbf{s}_{k}^{*\intercal}\mathbf{C}_{K}(\tau)^{\intercal}Q_{s}\mathbf{C}_{K}(\tau)\mathbf{s}_{k}^{*}\right)-\mathbf{s}_{k}^{*\intercal}\mathbf{S}^{\intercal}Q_{s}\mathbf{S}\mathbf{s}_{k}^{*}$
(45b)
$\displaystyle\leq\lambda_{2}(w_{k})+\mathbf{s}_{k}^{*\intercal}((1+\epsilon)\mathbf{C}_{K}(\tau)^{\intercal}Q_{s}\mathbf{C}_{K}(\tau)-\mathbf{S}^{\intercal}Q_{s}\mathbf{S})\mathbf{s}_{k}^{*}$
(45c)

where
$\mathbf{C}_{K}(\tau):=\sum_{j=1}^{n_{\Delta}}\tau_{j}\mathbf{C}_{K}^{j}$ and
$\mathbf{C}_{M}(\tau):=\sum_{j=1}^{n_{\Delta}}\tau_{j}\mathbf{C}_{M}^{j}$, and
(45b) is obtained by applying Cauchy-Schwarz theorem and Fenchel’s inequality
[30], by which,
$\left\|x+w\right\|\leq(1+\epsilon)\left\|x\right\|+(1+\frac{1}{\epsilon})\left\|w\right\|$
for any $\epsilon>0$. The function $\lambda_{2}(w_{k})$ in (45c) can chosen as
the first term in (45b). Moreover, using (40), it can be seen that

$\displaystyle\mathbf{S}^{\intercal}Q_{s}\mathbf{S}-P-(1{+}\epsilon)\mathbf{C}_{K}^{j\intercal}Q_{s}\mathbf{C}_{K}^{j}\succ
0,\quad\forall j\in\mathbb{N}_{1}^{n_{\Delta}}$ (46a)
$\displaystyle\implies\begin{bmatrix}(1{+}\epsilon)^{-1}Q_{s}^{-1}&\mathbf{C}_{K}^{j}\\\
\mathbf{C}_{K}^{j\intercal}&\mathbf{S}^{\intercal}Q_{s}\mathbf{S}-P\end{bmatrix}\succ
0,\>\forall j\in\mathbb{N}_{1}^{n_{\Delta}},$ (46b)
$\displaystyle\implies\begin{bmatrix}\displaystyle\sum_{j=1}^{n_{\Delta}}\tau_{j}(1{+}\epsilon)^{-1}Q_{s}^{-1}&\displaystyle\sum_{j=1}^{n_{\Delta}}\tau_{j}\mathbf{C}_{K}^{j}\\\
\displaystyle\sum_{j=1}^{n_{\Delta}}\tau_{j}\mathbf{C}_{K}^{j\intercal}&\displaystyle\sum_{j=1}^{n_{\Delta}}\tau_{j}(\mathbf{S}^{\intercal}Q_{s}\mathbf{S}-P)\end{bmatrix}\succ
0,$ (46c)
$\displaystyle\implies\begin{bmatrix}(1{+}\epsilon)^{-1}Q_{s}^{-1}&\mathbf{C}_{K}(\tau)\\\
\mathbf{C}_{K}(\tau)^{\intercal}&(\mathbf{S}^{\intercal}Q_{s}\mathbf{S}-P)\end{bmatrix}\succ
0,$ (46d)
$\displaystyle\implies\mathbf{S}^{\intercal}Q_{s}\mathbf{S}-P-(1{+}\epsilon)\mathbf{C}_{K}(\tau)^{\intercal}Q_{s}\mathbf{C}_{K}(\tau)\succ
0,$ (46e)

where the Schur complement lemma [31] is used to obtain (46b) and (46e). Using
(45c) and (46e),

$\displaystyle\begin{split}&V(x_{k+1})-V(x_{k})\leq\lambda_{2}(w_{k})-s_{k}^{*\intercal}Ps_{k}^{*}\\\
&\leq\lambda_{2}(w_{k})-x_{k}^{\intercal}P_{x}x_{k}=\lambda_{2}(w_{k})-\sigma_{3}(\left\|x_{k}\right\|).\end{split}$
(47)

Thus, $V(x)$ satisfies (43), and is an ISS Lyapunov function of the system.
Therefore, applying [32, Lemma 3.5], the closed-loop system is ISS. ∎

It can be seen that $\epsilon$ is used to define $\lambda_{2}(\cdot)$, which
affects the function $\lambda_{1}(\cdot)$ in (42) [32]. Thus, the choice of
$\epsilon$ affects the stability margin of the robust MPC controller, and can
be used as a tuning parameter by the designer.

###  6.5 Discussion

The central idea of this work is that RF is a one-step property, and can be
used to guarantee constraint satisfaction in closed loop as shown in Theorem
1. However, barring a few results in early MPC literature such as [21, 17],
this idea has not been fully utilized. Hence, some discussion is warranted on
interesting research directions.

####  6.5.1 Beyond constraint tightening

Theorem 1 does not restrict the robust MPC optimization to be of the structure
(41), which is a characteristic of constraint tightening techniques. The
theorem is also valid for any MPC optimization of the form

$\displaystyle\min_{\hat{\mathbf{x}}_{k},\hat{\mathbf{u}}_{k}}\quad\hat{x}_{N|k}^{\intercal}$
$\displaystyle
Q_{N}\hat{x}_{N|k}+\textstyle\sum_{i=0}^{N-1}\hat{x}_{i|k}^{\intercal}Q_{x}\hat{x}_{i|k}+\hat{u}_{i|k}^{\intercal}Q_{u}\hat{u}_{i|k}$
(48a) s.t. $\displaystyle\hat{x}_{k,0}=x_{k},$ (48b) $\displaystyle A$
$\displaystyle\hat{x}_{i|k}+B\hat{u}_{i|k}=\hat{x}_{i|k+1},\>\forall
i\in\mathbb{N}_{0}^{N-1},$ (48c) $\displaystyle
F\hat{x}_{0|k}+G\hat{u}_{0|k}\leq b,$ (48d)
$\displaystyle\Phi_{x}\hat{x}_{0|k}+\Phi_{u}\hat{\mathbf{u}}_{k}\leq\phi,$
(48e)

where $\Phi_{x},\Phi_{u}$ and $\phi$ can be chosen by the designer. In (48),
the system constraints (4) are only imposed on the first state and input in
the prediction horizon in (48d), which can be seen as selecting $t_{0}=0$ in
(5). Using a similar approach as in Theorem 3 and Proposition 1, sufficient
conditions can also be derived for SRF of (48). However, in practice, solving
the resulting non-convex optimization problem is harder, due to the higher
number of bilinearities.

####  6.5.2 Robust control invariance

A key concept which is known in control literature is that of robust control
invariance (RCI) [25, 8]. As opposed to robust invariance (RI), which
considers invariance of autonomous systems subject to uncertainty, RCI
considers controlled dynamical systems. For this reason, computation of RCI
sets is challenging [18, 19].

RCI provides an alternative interpretation of Theorem 2. A necessary and
sufficient condition for strong recursive feasibility of MPC is that the
feasible region is RCI. Because generic RCI sets are difficult to compute, the
proposed method uses parameterized feedback laws to formulate the problem as
computation of the resulting parameterized RI sets. Note that all RI sets are
also RCI sets, and are easier to compute in practice.

##  7 Illustrative example

In this section, Algorithm 1 will be used to generate a controller for a mass-
spring-damper system, where two masses are connected to each other by a spring
and a damper, and the spring constant and damping coefficient are uncertain.

The example has been adapted from [12], where similar systems with increasing
number of masses are considered, and a tube MPC strategy was used to design
the controller. The state of the system is a vector of the position and
velocity of each mass, and the control input is a vector of forces acting on
each mass. The system dynamics can be represented in discrete time by

$\displaystyle A$ $\displaystyle=\left[\begin{array}[]{cccc}1&T_{s}&0&0\\\
{-}\frac{k_{12}T_{s}}{m_{1}}&{-}\frac{c_{12}T_{s}}{m_{1}}{+}1&\frac{k_{12}T_{s}}{m_{1}}&\frac{c_{12}T_{s}}{m_{1}}\\\
0&0&1&T_{s}\\\
\frac{k_{12}T_{s}}{m_{2}}&\frac{c_{12}T_{s}}{m_{2}}&{-}\frac{k_{12}T_{s}}{m_{2}}&-\frac{c_{12}T_{s}}{m_{2}}{+}1\\\
\end{array}\right],$ (53) $\displaystyle B_{p}$
$\displaystyle{=}\left[\begin{array}[]{cccc}0&0\\\
\frac{k_{u}T_{s}}{m_{1}}&\frac{c_{u}T_{s}}{m_{1}}\\\ 0&0\\\
{-}\frac{k_{u}T_{s}}{m_{2}}&{-}\frac{c_{u}T_{s}}{m_{2}}\\\
\end{array}\right],B{=}\begin{bmatrix}0&0\\\ \frac{T_{s}}{m_{1}}&0\\\ 0&0\\\
0&\frac{T_{s}}{m_{2}}\end{bmatrix},B_{w}{=}w_{b}B,$ (58) $\displaystyle D_{x}$
$\displaystyle=\begin{bmatrix}{-}1&0\\\ 0&1\\\ \end{bmatrix},\quad
D_{u}=0,\quad D_{w}=0,$ (59)

where $m_{1}=m_{2}=0.2\mathrm{kg}$ is the mass, $k_{12}=0.5\mathrm{Nm^{-1}}$
and $c_{12}=0.5\mathrm{Nsm^{-1}}$ are the nominal spring constant and damping
coefficient. In addition, $T_{s}=0.1\mathrm{s}$ represents the sampling time
used for discretization, and $w_{b}=0.2$ represents the bound on a disturbance
acting on the velocity states. The magnitudes of all the states and inputs are
bounded by 2. The terms $k_{u}{=}0.04k_{12}$ and $c_{u}=0.02c_{12}$ model the
structured uncertainty matrix $B_{p}$. The matrices $\Delta^{j}$ defining
$\mathcal{P}$ in (2) were chosen as diagonal matrices of size $2\times 2$ with
$\pm 1$ on each diagonal entry, modeling a $\pm 4\%$ uncertainty for the
spring constant, and $\pm 2\%$ uncertainty for the damping coefficient.

Two controllers are designed for the system. The first uses Algorithm 1, and
is referred to as CLR, for closed-loop robust MPC. The second uses the method
proposed in [12], where an ellipsoidal tube is constructed at each time step
by the controller to bound the states and inputs in the prediction horizon.
This method will be referred to as ET, for ellipsoidal tube MPC. Both the
controllers use $N=5$ time steps. The terminal constraints of CLR are computed
using $K_{Y}$ chosen as the LQR gain, and $k^{\prime}=2$. The constants $\mu$
and $\epsilon$ were chosen as 2 and 0.1 respectively. The terminal cost is
computed such that it satisfies (40) and its deviation from the infinite
horizon LQR gain is minimized.

The initial state of the system is chosen to be
$x_{0}=[1.9,0.5,-1.7,1.7]^{\intercal}$. The disturbance affecting the system
is the same for both the controllers, and is sampled from a uniform
distribution in $\mathcal{W}$. The simulations were performed for 25
realizations of the spring constants, damping coefficients and disturbance
sequences chosen using a uniform distribution within the specified bounds.

It was observed that CLR requires higher offline computation time compared to
ET, and a much lower online computation time. This is because the ET algorithm
solves semi definite programs offline and online, whereas CLR solves a non-
convex optimization offline and a quadratic program online. All the
optimization problems were solved on a laptop with Intel i7-8550U processor,
and formulated using YALMIP [33]. The convex optimization problems were solved
using Gurobi [34] and MOSEK [35], and the non-convex optimization using IPOPT
[36]. The offline computation time for the CLR controller is $45.1\mathrm{s}$,
whereas that for ET is $0.06\mathrm{s}$. For solving (38), an initial guess
for the tightenings was computed as follows. Ignoring the perturbation from
$\Delta_{k}$, the approach from [22] was used to obtain a tightening. This
tightening was then scaled by a factor of 1.7, which was a feasible initial
guess for (38). The average online computation time for CLR is
$2.2\mathrm{ms}$, whereas that for ET is $53.9\mathrm{ms}$.

The closed-loop trajectories of the system are shown in Figure 4, where
highlighted regions show the upper and lower bounds of state and input
variables at each time step in closed loop over all the simulations. It can be
seen that both controllers successfully regulate the state to the origin
without constraint violations. However, the ET controller is aggressive at the
start of the simulation, because the initial position is too close to the
constraint, and the corresponding velocity is positive. Moreover, the ET
controller recursively outer-approximates the effect of worst-case
perturbations on state trajectories over the prediction horizon, resulting in
the aggressive behavior. The average closed-loop cost achieved by CLR is 83.0
and that for ET is 93.7, which is a $11.4\%$ cost reduction.

Figure 4:  Upper and lower bounds on the closed-loop trajectories with CLR and
ET MPC algorithms over 25 random realizations of disturbances and true system
parameters.

##  8 Conclusions

A novel design approach was presented for robust MPC, using the property that
recursive feasibility can imply constraint satisfaction in closed loop.
Recursive feasibility is a one-step property, meaning that it depends on
perturbations occurring at _one_ time step. Moreover, recursive feasibility
only requires the existence of a feasible input sequence at the subsequent
time step, as the MPC optimizer can compute the input sequence in closed loop.
Exploiting these ideas, it was shown that recursive feasibility can be imposed
as a constraint in the robust MPC design. Such a design reduces the
conservatism of most existing MPC methods, which ensure constraint
satisfaction by computing feasible open-loop trajectories under multiple
perturbations and disturbances occurring along the prediction horizon.

## References

  * [1] S. J. Qin and T. A. Badgwell, “A survey of industrial model predictive control technology,” _Control engineering practice_ , vol. 11, no. 7, pp. 733–764, 2003.
  * [2] G. Williams, P. Drews, B. Goldfain, J. M. Rehg, and E. A. Theodorou, “Information-theoretic model predictive control: Theory and applications to autonomous driving,” _IEEE Transactions on Robotics_ , vol. 34, no. 6, pp. 1603–1622, 2018.
  * [3] R. Pieters, S. Lombriser, A. Alvarez-Aguirre, and B. J. Nelson, “Model predictive control of a magnetically guided rolling microrobot,” _IEEE Robotics and Automation Letters_ , vol. 1, no. 1, pp. 455–460, 2016.
  * [4] F. Lejarza and M. Baldea, “Economic model predictive control for robust optimal operation of sparse storage networks,” _Automatica_ , vol. 125, p. 109346, 2021.
  * [5] J. B. Rawlings and D. Q. Mayne, _Model predictive control: theory, computation, and design_. Nob Hill Publishing Madison, WI, 2009.
  * [6] M. L. Darby and M. Nikolaou, “MPC: Current practice and challenges,” _Control Engineering Practice_ , vol. 20, no. 4, pp. 328–342, 2012.
  * [7] B. Kouvaritakis and M. Cannon, _Model predictive control: Classical, robust and stochastic_. New York, NY: Springer International Publishing, 2015.
  * [8] E. C. Kerrigan, “Robust constraint satisfaction: Invariant sets and predictive control,” Ph.D. dissertation, University of Cambridge UK, 2001.
  * [9] J. Löfberg, “Oops! I cannot do it again: Testing for recursive feasibility in MPC,” _Automatica_ , vol. 48, no. 3, pp. 550–555, 2012\.
  * [10] J. Fleming, B. Kouvaritakis, and M. Cannon, “Robust tube MPC for linear systems with multiplicative uncertainty,” _IEEE Transactions on Automatic Control_ , vol. 60, no. 4, pp. 1087–1092, 2014.
  * [11] J. Hu and B. Ding, “Output feedback robust MPC for linear systems with norm-bounded model uncertainty and disturbance,” _Automatica_ , vol. 108, p. 108489, 2019.
  * [12] A. Parsi, A. Iannelli, and R. S. Smith, “Scalable tube model predictive control of uncertain linear systems using ellipsoidal sets,” _International Journal of Robust and Nonlinear Control_ , 2022, doi: 10.1002/rnc.6485.
  * [13] D. Q. Mayne, M. M. Seron, and S. V. Raković, “Robust model predictive control of constrained linear systems with bounded disturbances,” _Automatica_ , vol. 41, no. 2, pp. 219–224, 2005.
  * [14] L. Schwenkel, J. Köhler, M. A. Müller, and F. Allgöwer, “Model predictive control for linear uncertain systems using integral quadratic constraints,” _IEEE Transactions on Automatic Control_ , 2022.
  * [15] M. Maiworm, T. Bäthge, and R. Findeisen, “Scenario-based model predictive control: Recursive feasibility and stability,” _IFAC-PapersOnLine_ , vol. 48, no. 8, pp. 50–56, 2015.
  * [16] L. Chisci and G. Zappa, “Robustifying a predictive controller against persistent disturbances,” in _1999 European Control Conference (ECC)_ , 1999, pp. 2419–2424.
  * [17] E. C. Kerrigan and J. M. Maciejowski, “Robust feasibility in model predictive control: necessary and sufficient conditions,” in _Proceedings of the 40th IEEE Conference on Decision and Control (Cat. No.01CH37228)_ , vol. 1, 2001, pp. 728–733 vol.1.
  * [18] S. V. Rakovic, E. C. Kerrigan, D. Q. Mayne, and K. I. Kouramas, “Optimized robust control invariance for linear discrete-time systems: Theoretical foundations,” _Automatica_ , vol. 43, no. 5, pp. 831–841, 2007.
  * [19] M. Rungger and P. Tabuada, “Computing robust controlled invariant sets of linear systems,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 7, pp. 3665–3670, 2017.
  * [20] Y. Abdelsalam, S. Subramanian, and S. Engell, “Synthesis of computationally efficient recursively feasible multi-stage predictive controllers for uncertain linear systems,” _IEEE Control Systems Letters_ , vol. 6, pp. 1796–1801, 2021.
  * [21] L. Chisci, J. A. Rossiter, and G. Zappa, “Systems with persistent disturbances: Predictive control with restricted constraints,” _Automatica_ , vol. 37, no. 7, pp. 1019–1028, 2001.
  * [22] A. Parsi, P. Anagnostaras, A. Iannelli, and R. S. Smith, “Computationally efficient robust MPC using optimized constraint tightening,” in _61st IEEE Conference on Decision and Control_ , 2022.
  * [23] J. Doyle, A. Packard, and K. Zhou, “Review of LFTs, LMIs, and mu,” in _Proceedings of the 30th IEEE Conference on Decision and Control_ , vol. 2, 1991, pp. 1227–1232.
  * [24] F. Borrelli, A. Bemporad, and M. Morari, _Predictive control for linear and hybrid systems_. Cambridge University Press, 2017.
  * [25] F. Blanchini and S. Miani, _Set-theoretic methods in control_. Boston: Birkhäuser, 2008.
  * [26] 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.
  * [27] J. Lawrence, “Polytope volume computation,” _Mathematics of computation_ , vol. 57, no. 195, pp. 259–271, 1991.
  * [28] D. Limon, T. Alamo, D. M. Raimondo, D. M. D. L. Pena, J. M. Bravo, A. Ferramosca, and E. F. Camacho, “Input-to-state stability: A unifying framework for robust model predictive control,” _Nonlinear model predictive control_ , pp. 1–26, 2009.
  * [29] A. Bemporad, M. Morari, V. Dua, and E. Pistikopoulos, “The explicit linear quadratic regulator for constrained systems,” _Automatica_ , vol. 38, no. 1, pp. 3–20, 2002.
  * [30] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [31] F. Zhang, _The Schur complement and its applications_. Springer Science & Business Media, 2006, vol. 4.
  * [32] Z. Jiang and Y. Wang, “Input-to-state stability for discrete-time nonlinear systems,” _Automatica_ , vol. 37, no. 6, pp. 857–869, 2001.
  * [33] J. Löfberg, “YALMIP : A toolbox for modeling and optimization in MATLAB,” _In Proceedings of the CACSD Conference_ , 2004.
  * [34] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2022. [Online]. Available: https://www.gurobi.com
  * [35] MOSEK ApS, _The MOSEK optimization toolbox for MATLAB manual. Version 9.0._ , 2019. [Online]. Available: http://docs.mosek.com/9.0/toolbox/index.html
  * [36] 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.


# Are all Dirac-harmonic maps uncoupled?

Bernd Ammann  B. Ammann, Fakultät für Mathematik, Universität Regensburg,
93040 Regensburg, Germany<EMAIL_ADDRESS>
###### Abstract.

Dirac-harmonic maps $(f,{\varphi})$ consist of a map $f:M\to N$ and a twisted
spinor ${\varphi}\in\Gamma(\Sigma M\otimes f^{*}TN)$ and they are defined as
critical points of the super-symmetric energy functional. A Dirac-harmonic map
is called _uncoupled_ , if $f$ is a harmonic map. We show that under some
minimality assumption Dirac-harmonic maps defined on a closed domain are
uncoupled.

The author was supported by SPP 2026 (Geometry at infinity) and the SFB 1085
(Higher Invariants), both funded by the DFG (German Science Foundation).

MSC: 58E20 (Primary), 53C43, 53C27 (Secondary)

##  1\. Introduction

Let $M$ and $N$ be Riemannian manifolds, and $f:M\to N$ a $C^{1}$. If $M$ is
compact, we can define the energy of $f$ as

$\mathcal{E}_{1}(f)=\frac{1}{2}\int_{M}|df|^{2}\,\mathrm{dvol}^{g}$

considered as a functional on $C^{1}(M,N)$. Critical points of this functional
are called _harmonic maps_ , and there is an extensive literature about
harmonic maps, with many interesting applications, also weak solutions were
studied.

In the recent years, there is a growing number of publications about a super-
symmetric analogue of harmonic maps, called _Dirac-harmonic maps_. They are
defined as critical points of the sypersymmetric energy functional
$\mathcal{E}$ defined in (1), see Section 2 for more details. In particular, a
Dirac-harmonic map is a pair $(f,{\varphi})$ of a map $f:M\to N$ and a twisted
spinor ${\varphi}\in\Gamma(\Sigma M\otimes f^{*}TN)$. An early article [6]
about such maps was written in 2005 by Chen, Jost, Li and Wang, studying
regularity issues for Dirac-harmonic maps, followed by [7, 9] and stimulating
many associated questions and results. Between 2005 and September 2022,
MathSciNet lists about 45 publications with “Dirac-harmonic” in the title.
Many questions answered for harmonic maps can be discussed in the Dirac-
harmonic context, often it is hard to include the spinorial part into the
estimates, and good progress was achieved.

Obiously, Dirac-harmonic maps with ${\varphi}\equiv 0$ are uninteresting, as
then $(f,{\varphi})$ is Dirac-harmonic if and only if $f$ is harmonic; such
solution are called _spinor-trivial_. Similarly, solutions with $f$ constant,
called _map-trivial_ solutions, are not in the focus of interest within this
subject either; in this case the problem is equivalent to finding harmonic
spinors in the classical sense, see e. g. [10, 4, 1]. We say that a Dirac-
harmonic map is _uncoupled_ if $f$ is harmonic, otherwise we say that this
Dirac-harmonic map is _coupled_. Being uncoupled is equivalent to the
vanishing of the $\mathcal{R}$-term, see Section 2. A major question discussed
in this article is whether coupled Dirac-harmonic maps with compact domain
exist. Let us discuss examples of Dirac-harmonic maps in the literature first.

Progress in constructing Dirac-harmonic maps was achieved in [12] by using
(untwisted) harmonic spinors, twistor spinors and similar solutions of other
spinorial equations, although it remained unclear, how one could get solutions
of these spinorial equations. This was analyzed later by Ginoux and the author
in [3]: we showed – see [3, Theorems 1.1 and 1.3] – that the conditions in
[12, Theorems 1 and 3] can be satisfied only in exceptional cases.111Note that
the proof of [12, Theorems 1] uses that the immersion is isometric, although
this is not stated in this theorem [12, Theorems 1] explicitly. In particular,
for $\dim M\geq 3$ and $M$ complete, one concludes that $M$ has to be simply-
connected of constant negative curvature. Other strong obstructions exist for
$\dim M=2$, in particular $f$ has to be necessarily harmonic, thus any Dirac-
harmonic map with compact domain obtained this way is uncoupled. The solutions
in [12, Theorems 2], reproven as [3, Corollary 2.3] are uncoupled as well.

Many more Dirac-harmonic maps $(f,{\varphi})$ were constructed in [2]. Here,
the domain $M$ of $f$ is closed. The method starts with a harmonic map $f:M\to
N$ and then one uses index theory to get a non-vanishing harmonic spinor
${\varphi}\in\Gamma(\Sigma M\otimes f^{*}T^{*}N)$. In particular, these
solutions are uncoupled.

In the special case $M=N=S^{2}$ with the round metric, it was shown in [18]
that every Dirac-harmonic map is uncoupled, and this was recently reproven in
[13, Proposition 1.1].

Summarizing this, the author is not aware of any publication in which the
existence of any coupled Dirac-harmonic map with compact domain was proven.

In the current article, we will show, that “generically” every Dirac-harmonic
map with closed domain is uncoupled, see Corollary 3 for details. The argument
may be adapted to compact manifolds $M$ with boundary for many suitable
homogeneous boundary conditions, an extension of the results that we will not
work out.

It thus remains questionable whether coupled Dirac-harmonic maps with closed
domain actually exist. Suppose such solutions did not exist, then the main
results in several recent publications would not provide new result, compared
to what is known already for harmonic maps. However, spectral flow methods
might possibly be used to construct coupled Dirac-harmonic maps with closed
domain. Thus it is hard to predict whether an overseen reason implies
$\mathcal{R}=0$ for $M$ closed, or whether we did not yet find the right tools
to find coupled solutions.

What we do not discuss here. In order to avoid missunderstandings, we want to
add here some issues that we _do not_ prove in this article. First, there is a
variant of the super-symmetric energy functional that includes a quartic
spinor term, depending on the curvature of the target manifold $N$, see e. g.
[11, (2.4)]. Orally, I was told, that this variant is even more important from
the perspective of applications to physics. I am unable to adapt the arguments
of the current article to this modified functional. On the other hand I expect
that many publications about Dirac-harmonic maps are sufficiently robust in
order to be applicable also to critical points of this modified functional.
Similar arguments apply to other perturbations of the supersymmetric energy
functional. As a consequence, the techniques in the articles listed above
would provide valuable contributions even when a non-existence result for
coupled solution on closed domains were finally proven.

Second, we did not yet discuss the case of non-compact domains. Our argument
relies essentially on the fact that the Euler-Lagrange equations (2) describe
the stationary points of a well-defined functional, that $\not{D}^{f}$ is
self-adjoint and that we do not get boundary terms by partial integration. In
[12] and [3] an example of a Dirac-harmonic map with non-vanishing
$\mathcal{R}$-term is (implicitly) given: here $M^{m}$ is a simply connected
complete Riemannian manifold of constant negative sectional curvature
$-4/(m+2)$ (i. e. a space form), and $f$ is an isometric embedding into
$(m+1)$-dimensional hyperbolic space $N$, for which $f(M)$ is totally umbillic
with parallel shape tensor in $N$. Due to the conformal covariance of the
Dirac operator and of the Penrose operator, $M$ carries many twistor spinors
and harmonic spinors, thus the Jost-Mo-Zhu method [12] may be applied, and one
easily checks that $\tau(f)$ is a non-zero constant. However, in this case
$|df|^{2}$ is not integrable, and thus the super-symmetric energy functional
$\mathcal{E}$ does not converge. The Dirac-harmonic map $(f,{\varphi})$ solves
the system of partial differential equations (2) without being the
stationarity equation of a functional defined on all pairs of maps with
spinors.

Third, we only considered solutions to the Dirac-harmonic map equation in the
strong sense. We do not know to which extend our methods also generalize to
weak settings. For the harmonic maps, weak solutions gave rise to involved
research, see e.g. Bethuel’s work [5].

Acknowledgements. The result of this note grew out of discussions related to a
pair of conferences in Shanghai and Sanya in spring 2018, organized by Jürgen
Jost, Miaomiao Zhu and others. My thanks goes to these organizers. I also
thank my former PhD student Johannes Wittmann for his contributions. Thanks
also to Christian Bär and Johannes Wittmann for comments on an earlier version
of this preprint. Thanks also to Lei Liu for pointing out that I had overseen
a boundary term in a previous version.

##  2\. Dirac-harmonic maps with closed domain are generically uncoupled

Let $(M,g)$ be a compact Riemannian spin manifold, and $(N,h)$ a Riemannian
manifold. For a map $f:M\to N$ and a twisted spinor ${\varphi}\in\Gamma(\Sigma
M\otimes f^{*}TN)$ we define

$\displaystyle\mathcal{E}_{1}(f)$ $\displaystyle:=$
$\displaystyle\frac{1}{2}\int_{M}|df|^{2}\,\mathrm{dvol}^{g}$
$\displaystyle\mathcal{E}_{2}(f,{\varphi})$ $\displaystyle:=$
$\displaystyle\frac{1}{2}\int_{M}\langle{\varphi},\not{D}^{f}{\varphi}\rangle\,\mathrm{dvol}^{g}$
(1) $\displaystyle\mathcal{E}(f,{\varphi})$ $\displaystyle:=$
$\displaystyle\mathcal{E}_{1}(f)+\mathcal{E}_{2}(f,{\varphi})\,.$

The functional $\mathcal{E}$ is called the _super-symmetric energy
functional_. We define

$\displaystyle\tau(f)$ $\displaystyle:=$
$\displaystyle-\frac{\partial\mathcal{E}_{1}}{\partial f}=\mathrm{tr}\nabla
df\in\Gamma(f^{*}TN)$ $\displaystyle\mathcal{R}(f,{\varphi})$
$\displaystyle:=$ $\displaystyle\frac{\partial\mathcal{E}_{2}}{\partial f}$
$\displaystyle\Psi(f,{\varphi})$ $\displaystyle:=$
$\displaystyle\frac{\partial\mathcal{E}_{2}}{\partial{\varphi}}=\not{D}^{f}{\varphi}\,.$

In these equations, $\not{D}^{f}$ denotes the twisted Dirac operator acting on
$\Gamma(\Sigma M\otimes f^{*}TN)$. Further, we used the canonical
$L^{2}$-scalar product on $\Gamma(f^{*}TN)$ and on $\Gamma(\Sigma\otimes
f^{*}TN)$, in order to identify the partial derivatives of $\mathcal{E}_{1}$
and $\mathcal{E}_{2}$, given above, with the corresponding gradients.

Stationary points of $\mathcal{E}$ are called _Dirac-harmonic maps_. Thus this
condition is equivalent to

(2) $\not{D}^{f}{\varphi}=0\text{ and }\tau(f)=\mathcal{R}(f,{\varphi}).$

We say that $(f,{\varphi})$ is _uncoupled_ , if $\mathcal{R}(f,{\varphi})=0$.
In this publication we only consider _smooth_ Dirac-harmonic maps.

######  Remark.

Note that $\Sigma M\otimes f^{*}TN$ will always denote the real tensor
product, even when $f^{*}TN$ and $\Sigma M$ carry natural complex structures.
Note that $\Sigma M$ may be defined as a real, a complex, a quaternionic or a
$\mathrm{Cl}_{m}$-linear Dirac operator, this does not matter. For simplicity
we only restrict to the complex version in this article, as it is classically
the most studied one.

######  Remark.

The terminology for the $\mathcal{R}$-term is standard, see e.g. [12, Page
1514, after (9)], [8, (1.3)].

######  Theorem 1.

Let $(f_{0},{\varphi}_{0})$ be given. We assume that there is an open
neighborhood $U$ of $f_{0}$, and a $C^{1}$-map $U\ni
f\mapsto\hat{\varphi}(f)$, such that $\not{D}^{f}(\hat{\varphi}(f))=0$ and
such that $\hat{\varphi}(f_{0})={\varphi}_{0}$. Then $(f_{0},{\varphi}_{0})$
is uncoupled.

###### Proof.

Obviously we have $\mathcal{E}_{2}(f,\hat{\varphi}(f))=0$. In the following we
write $d/df$ for the total222In the sense used in Lagrangian mechanics: we
derive the expression $f\mapsto\mathcal{E}(f,\hat{\varphi}(f))$ derivative
with respect to $f$ and – as before – $\partial/\partial f$ for the partial
derivative with respect to $f$. We thus get

$\displaystyle-\tau(f)$ $\displaystyle=$
$\displaystyle\frac{d}{df}\mathcal{E}_{1}(f)$ $\displaystyle=$
$\displaystyle\frac{d}{df}\Bigl{(}\mathcal{E}\bigl{(}f,\hat{\varphi}(f)\bigr{)}\Bigr{)}$
$\displaystyle=$ $\displaystyle\frac{\partial\mathcal{E}}{\partial
f}\mathchoice{\lower 3.9pt\hbox{\Big{|}}\lower 7.40005pt\hbox{$\scriptstyle
f$}}{\lower 3.30002pt\hbox{\big{|}}\lower 6.60004pt\hbox{$\scriptstyle
f$}}{\lower 3.00003pt\hbox{\big{|}}\lower 6.00006pt\hbox{$\scriptscriptstyle
f$}}{\lower 2.5pt\hbox{\big{|}}\lower 5.09995pt\hbox{$\scriptscriptstyle
f$}}+\frac{\partial\mathcal{E}}{\partial{\varphi}}\mathchoice{\lower
3.9pt\hbox{\Big{|}}\lower
7.40005pt\hbox{$\scriptstyle(f,\hat{\varphi}(f))$}}{\lower
3.30002pt\hbox{\big{|}}\lower
6.60004pt\hbox{$\scriptstyle(f,\hat{\varphi}(f))$}}{\lower
3.00003pt\hbox{\big{|}}\lower
6.00006pt\hbox{$\scriptscriptstyle(f,\hat{\varphi}(f))$}}{\lower
2.5pt\hbox{\big{|}}\lower
5.09995pt\hbox{$\scriptscriptstyle(f,\hat{\varphi}(f))$}}\,\frac{\partial\hat{\varphi}}{\partial
f}\mathchoice{\lower 3.9pt\hbox{\Big{|}}\lower 7.40005pt\hbox{$\scriptstyle
f$}}{\lower 3.30002pt\hbox{\big{|}}\lower 6.60004pt\hbox{$\scriptstyle
f$}}{\lower 3.00003pt\hbox{\big{|}}\lower 6.00006pt\hbox{$\scriptscriptstyle
f$}}{\lower 2.5pt\hbox{\big{|}}\lower 5.09995pt\hbox{$\scriptscriptstyle f$}}$
$\displaystyle=$ $\displaystyle\frac{\partial\mathcal{E}_{1}}{\partial
f}\mathchoice{\lower 3.9pt\hbox{\Big{|}}\lower 7.40005pt\hbox{$\scriptstyle
f$}}{\lower 3.30002pt\hbox{\big{|}}\lower 6.60004pt\hbox{$\scriptstyle
f$}}{\lower 3.00003pt\hbox{\big{|}}\lower 6.00006pt\hbox{$\scriptscriptstyle
f$}}{\lower 2.5pt\hbox{\big{|}}\lower 5.09995pt\hbox{$\scriptscriptstyle
f$}}+\frac{\partial\mathcal{E}_{2}}{\partial f}\mathchoice{\lower
3.9pt\hbox{\Big{|}}\lower
7.40005pt\hbox{$\scriptstyle(f,\hat{\varphi}(f))$}}{\lower
3.30002pt\hbox{\big{|}}\lower
6.60004pt\hbox{$\scriptstyle(f,\hat{\varphi}(f))$}}{\lower
3.00003pt\hbox{\big{|}}\lower
6.00006pt\hbox{$\scriptscriptstyle(f,\hat{\varphi}(f))$}}{\lower
2.5pt\hbox{\big{|}}\lower
5.09995pt\hbox{$\scriptscriptstyle(f,\hat{\varphi}(f))$}}+\underbrace{\frac{\partial\mathcal{E}_{2}}{\partial{\varphi}}\mathchoice{\lower
3.9pt\hbox{\Big{|}}\lower
7.40005pt\hbox{$\scriptstyle(f,\hat{\varphi}(f))$}}{\lower
3.30002pt\hbox{\big{|}}\lower
6.60004pt\hbox{$\scriptstyle(f,\hat{\varphi}(f))$}}{\lower
3.00003pt\hbox{\big{|}}\lower
6.00006pt\hbox{$\scriptscriptstyle(f,\hat{\varphi}(f))$}}{\lower
2.5pt\hbox{\big{|}}\lower
5.09995pt\hbox{$\scriptscriptstyle(f,\hat{\varphi}(f))$}}}_{=\Psi(f,\hat{\varphi}(f))=0}\frac{\partial\hat{\varphi}}{\partial
f}\mathchoice{\lower 3.9pt\hbox{\Big{|}}\lower 7.40005pt\hbox{$\scriptstyle
f$}}{\lower 3.30002pt\hbox{\big{|}}\lower 6.60004pt\hbox{$\scriptstyle
f$}}{\lower 3.00003pt\hbox{\big{|}}\lower 6.00006pt\hbox{$\scriptscriptstyle
f$}}{\lower 2.5pt\hbox{\big{|}}\lower 5.09995pt\hbox{$\scriptscriptstyle f$}}$
$\displaystyle=$ $\displaystyle-\tau(f)+\mathcal{R}(f,\hat{\varphi}(f)).$

and this obviously implies $\mathcal{R}(f_{0},{\varphi}_{0})=0$. ∎

######  Corollary 2.

Suppose $f_{0}$ has a neighborhood $U$ such that for all $f\in U$ we have
$\dim\not{D}^{f}\geq\dim\not{D}^{f_{0}}$. Then for every
${\varphi}_{0}\in\mathrm{ker}\not{D}^{f_{0}}$ the pair $(f_{0},{\varphi}_{0})$
is uncoupled.

Or as a special case we get.

######  Corollary 3.

Suppose $(f_{0},{\varphi}_{0})$ is a Dirac-harmonic map with $M$ closed. Then

  1. (1)

$(f_{0},{\varphi}_{0})$ is uncoupled, or

  2. (2)

any neighborhood $U$ of $f_{0}$ contains an $f\in U$ with
$\dim\not{D}^{f}<\dim\not{D}^{f_{0}}$.

This corollary shows that the construction of coupled Dirac-harmonic maps with
compact domain is difficult. This fits to the fact that to the author’s
knowledge no coupled Dirac-harmonic map with closed domain has been found so
far.

Note that these results may be easily generalized to compact manifolds with
boundary with suitable homogeneous boundary conditions. However, our methods
are not robust enough to generalize to more general type of functionals
$\mathcal{E}$, e. g. to the one with quartic spinor term discussed in [11,
(2.4)].

######  Remark 4.

For studying harmonic maps, Sacks and Uhlenbeck [17] applied a perturbation of
the energy functional successfully, namely they defined for $\alpha\in(1,2)$

$\mathcal{E}_{\alpha}(f)=\frac{1}{2}\int_{M}(1+|df|^{2})^{\alpha}\,\mathrm{dvol}^{g}.$

Recently, a similar perturbation was also used for the super-symmetric energy
functional, i. e. one defines
$\mathcal{E}(f,{\varphi})=\mathcal{E}_{\alpha}(f)+\mathcal{E}_{2}(f,{\varphi})$,
see e. g. [14, Eq. (1.2)], [15] and [16]. Obviously Theorem 1 remains true for
this modification. The condition of Theorem 1 are satisfied as index
theoretical methods and minimal kernel methods are used. As a consequence all
$\alpha$-Dirac-harmonic maps in [15] are uncoupled and most of the results for
($\alpha$)-Dirac-harmonic maps in [15] directly follow from the corresponding
statements for ($\alpha$)-harmonic maps of the same article. In the results in
[14, Eq. (1.1)] an additional perturbation $F$ is allowed, which makes our
trick non-applicable, while it applies to stationary points of (1.2) in the
same article. Again, it seems unanswered whether coupled $\alpha$-Dirac-
harmonic maps with closed domain exist.

## References

  * [1] Ammann, B., Dahl, M., and Humbert, E. Surgery and harmonic spinors.  Adv. Math. 220 (2009), 523–539.
  * [2] Ammann, B., and Ginoux, N. Dirac-harmonic maps from index theory.  Calc. Var. Partial Differ. Equ. 47, 3-4 (2013), 739–762.
  * [3] Ammann, B., and Ginoux, N. Some examples of Dirac-harmonic maps.  Lett. Math. Phys. 109, 5 (2019), 1205–1218.
  * [4] Bär, C. Metrics with harmonic spinors.  Geom. Funct. Anal. 6 (1996), 899–942.
  * [5] Bethuel, F. On the singular set of stationary harmonic maps.  Manuscripta Math. 78, 4 (1993), 417–443.
  * [6] Chen, Q., Jost, J., Li, J., and Wang, G. Regularity theorems and energy identities for Dirac-harmonic maps.  Math. Z. 251, 1 (2005), 61–84.
  * [7] Chen, Q., Jost, J., Li, J., and Wang, G. Dirac-harmonic maps.  Math. Z. 254, 2 (2006), 409–432.
  * [8] Chen, Q., Jost, J., Sun, L., and Zhu, M. Estimates for solutions of Dirac equations and an application to a geometric elliptic-parabolic problem.  J. Eur. Math. Soc. (JEMS) 21, 3 (2019), 665–707.
  * [9] Chen, Q., Jost, J., and Wang, G. Liouville theorems for Dirac-harmonic maps.  J. Math. Phys. 48, 11 (2007), 113517, 13.
  * [10] Hitchin, N. Harmonic spinors.  Adv. Math. 14 (1974), 1–55.
  * [11] Jost, J., Liu, L., and Zhu, M. Energy identity and necklessness for a nonlinear supersymmetric sigma model.  Calc. Var. Partial Differential Equations 61, 3 (2022), Paper No. 112, 26.
  * [12] Jost, J., Mo, X., and Zhu, M. Some explicit constructions of Dirac-harmonic maps.  J. Geom. Phys. 59, 11 (2009), 1512–1527.
  * [13] Jost, J., Sun, L., and Zhu, J. Dirac-harmonic maps with trivial index.  arXiv: 2209.04106, 2022.
  * [14] Jost, J., and Zhu, J. $\alpha$-Dirac-harmonic maps from closed surfaces.  Calc. Var. Partial Differential Equations 60, 3 (2021), Paper No. 111, 41.
  * [15] Jost, J., and Zhu, J. Existence of (Dirac-)harmonic maps from degenerating (spin) surfaces.  J. Geom. Anal. 31, 11 (2021), 11165–11189.
  * [16] Li, J., Liu, L., Zhu, C., and Zhu, M. Energy identity and necklessness for $\alpha$-Dirac-harmonic maps into a sphere.  Calc. Var. Partial Differential Equations 60, 4 (2021), Paper No. 146, 19.
  * [17] Sacks, J., and Uhlenbeck, K. The existence of minimal immersions of $2$-spheres.  Ann. of Math. (2) 113, 1 (1981), 1–24.
  * [18] Yang, L. A structure theorem of Dirac-harmonic maps between spheres.  Calc. Var. Partial Differential Equations 35, 4 (2009), 409–420.