# A novel constraint tightening approach for robust data-driven predictive
control†††thanks: †This work was supported by Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under grant MU 3929/1-2 and AL 316/12-2 -
279734922, under Germany’s Excellence Strategy - EXC 2075 - 390740016, and
under grant 468094890. We acknowledge the support by the Stuttgart Center for
Simulation Science (SimTech). The authors thank the International Max Planck
Research School for Intelligent Systems (IMPRS-IS) for supporting Julian
Berberich.

Christian Klöppelt∗, Julian Berberich∗∗, Frank Allgöwer∗∗, Matthias A. Müller∗
∗Christian Klöppelt and Matthias A. Müller are with the Institute of Automatic
Control, Leibniz University Hannover, Germany.∗∗Julian Berberich and Frank
Allgöwer are with the Institute for Systems Theory and Automatic Control,
University of Stuttgart, Germany

###### Abstract

In this paper, we present a data-driven model predictive control (MPC) scheme
that is capable of stabilizing unknown linear time-invariant systems under the
influence of process disturbances. To this end, Willems’ lemma is used to
predict the future behavior of the system. This allows the entire scheme to be
set up using only a priori measured data and knowledge of an upper bound on
the system order. First, we develop a state-feedback MPC scheme, based on
input-state data, which guarantees closed-loop practical exponential stability
and recursive feasibility as well as closed-loop constraint satisfaction. The
scheme is extended by a suitable constraint tightening, which can also be
constructed using only data. In order to control a priori unstable systems,
the presented scheme contains a pre-stabilizing controller and an associated
input constraint tightening. We first present the proposed data-driven MPC
scheme for the case of full state measurements, and also provide extensions
for obtaining similar closed-loop guarantees in case of output feedback. The
presented scheme is applied to a numerical example.

Keywords— data-driven MPC, robust MPC

##  I Introduction

In recent years, there has been significant interest in designing data-driven
model predictive control (MPC) schemes, in which predictions are not based on
a parametric model of the system, but rather directly on a priori collected
input/output data from the system, thus circumventing the challenging
intermediate step of finding an accurate model. This is done by employing the
so-called Willems’ fundamental lemma[1], which states that for a controllable
linear system, all possible system trajectories can be parametrized in terms
of linear combinations of time-shifts of one single, persistently exciting,
trajectory.

A direct data-based MPC scheme based on Willems’ lemma was first considered by
Yang et al.[2] and Coulson et al.[3]. Guarantees for recursive feasibility,
stability, and robustness (in the presence of measurement noise) of the closed
loop were first proven by Berberich et al.[4]. In recent years, various
further properties and extensions of this data-driven MPC framework have been
studied, compare, e.g., the works by Coulson et al.[5], Huang et al.[6], Yin
et al.[7, 8], Xue and Matni[9], Furieri et al.[10], Berberich et al.[11] and
the overview paper by Markovsky and Dörfler[12].

One of the major strengths of MPC is its ability to take constraints in the
optimization problem into account, and therefore, guarantee their satisfaction
in closed-loop operation. For data-driven MPC schemes, achieving constraint
satisfaction is similarly important. However, in practice we typically only
have access to noisy data. Thus, to achieve closed-loop constraint
satisfaction, a suitable constraint tightening is required, similar to model-
based MPC[13, 14, 15]. Such a constraint tightening takes into account
possible (worst-case) disturbances as well as their influence on the system
dynamics in order to ensure that no disturbance that may occur in the future
can result in constraint violation.

For data-driven MPC schemes relying on an a priori identification of the
system model, there already exist schemes that provide a proper constraint
tightening even in the case of additive process noise[16, 17]. However, in the
direct data-driven setting –based on Willems’ lemma– this problem has not been
conclusively solved so far. In Berberich et al.[18], closed-loop constraint
satisfaction is shown in case of measurement noise; however, no process noise
and no input constraint tightening is considered, resulting in the fact that
the proposed scheme can only be applied (without being overly conservative) to
open-loop stable systems. Process disturbances acting additively on the
dynamics have been considered by Huang et al. [19] and Umenberger et al.[20].
However, both schemes lack the aforementioned closed-loop guarantees, and
moreover, rely on the knowledge of a priori measured disturbances. Recently,
Liu et al.[21] proposed a scheme guaranteeing closed-loop stability and
recursive feasibility in the presence of process disturbances, which, however,
lacks of guarantees for closed-loop constraint satisfaction.

In this paper, we propose a data-driven MPC scheme for which recursive
feasibility, robust stability, and closed-loop constraint satisfaction can be
ensured in the presence of process disturbances. In particular, we adapt the
state/output constraint tightening originally proposed by Berberich et al.[18]
to the case of process disturbances and propose a suitable additional input
constraint tightening to allow for a pre-stabilizing feedback, such that the
scheme can also be applied to open-loop unstable systems. We first consider
the case where full state measurement is available, before presenting
extensions to the output feedback case. We discuss how certain constants
required for the proposed constraint tightening can be computed from a priori
collected data. For both schemes, guarantees for practical exponential
stability, recursive feasibility and closed-loop constraint satisfaction are
proven.

The remainder of the paper is structured as follows. In Section II the problem
setup and preliminaries, such as the concept of persistency of excitation and
Willems’ lemma, are introduced. Next, in Section III, we set up the first
proposed data-driven MPC scheme based on state measurements. To this end, we
elaborate the data-driven parametrization of the constraint tightening,
explain the MPC scheme and prove the aforementioned closed-loop guarantees.
Thereafter, in Section IV we present data-driven output-feedback MPC based on
the case where only output measurements are available. We apply the proposed
scheme to a numerical example in Section V, and end with some concluding
remarks in Section VI.

_Notation:_ For a sequence $\\{z_{k}\\}_{k=0}^{N-1}$, we define the Hankel
matrix of depth $L$ as

$H_{L}(z)=\begin{bmatrix}z_{0}&z_{1}&\dots&z_{N-L}\\\
z_{1}&z_{2}&\dots&z_{N-L+1}\\\ \vdots&\vdots&\ddots&\vdots\\\
z_{L-1}&z_{L}&\dots&z_{N-1}\end{bmatrix},$

and the stacked window from time instant $a$ to $b$ as

$z_{\left[a,b\right]}=\begin{bmatrix}z_{a}\\\ \vdots\\\ z_{b}\end{bmatrix}.$

##  II Problem setup and preliminaries

In this paper, we consider the discrete-time multi-input multi-output LTI
system

$\begin{split}x_{k+1}&=Ax_{k}+Bu_{k}+w_{k},\\\
y_{k}&=Cx_{k}+Du_{k},\end{split}$ (1)

with the state $x_{k}\in\mathbb{R}^{n}$, the input $u_{k}\in\mathbb{R}^{m}$,
the output $y_{k}\in\mathbb{R}^{p}$ and the process disturbance
$w_{k}\in\mathbb{R}^{n}$. The setup can be extended to include measurement
noise as well, see Remark 1 below. Throughout the paper, we assume that (1) is
a minimal realization, i.e., that the pair $\left(A,B\right)$ is controllable
and the pair $\left(A,C\right)$ is observable. Moreover, we consider the
matrices $A,B,C,D$ as being unknown and the only knowledge about the system
available being its order $n$.

Moreover, we assume that the process disturbances belongs to a hypercube
(precise definitions will be given in Sections III and IV). The goal of this
paper is to construct a data-driven MPC scheme that stabilizes the origin and
ensures input and state constraint satisfaction (cf. Section III) or output
constraint satisfaction (cf. Section IV), where the respective constraint sets
are given by hypercubes.

To this end, we apply a persistently exciting (p.e.) input sequence to the
system, and measure the resulting state/output sequence, where a persistently
exciting sequence is defined as follows.

######  Definition 1.

A sequence $\\{u_{k}\\}_{k=0}^{N-1}$, with $u_{k}\in\mathbb{R}^{m}$, is
persistently exciting of order $L$ if $\mathrm{rank}\left(H_{L}(u)\right)=mL$.

We want to make use of Willems’ fundamental lemma for the prediction in an MPC
problem.

######  Lemma 1 (Willems’ lemma [1]).

Suppose $\\{u_{k},\hat{y}_{k}\\}_{k=0}^{N-1}$ is a trajectory of the
controllable system

$\begin{split}\hat{x}_{k+1}&=A\hat{x}_{k}+Bu_{k},\\\
\hat{y}_{k}&=C\hat{x}_{k}+Du_{k},\end{split}$ (2)

and $u$ is persistently exciting of order $L+n$. Then,
$\\{\bar{u}_{k},\bar{y}_{k}\\}_{k=0}^{L-1}$ is a trajectory of System (2) if
and only if there exists $\alpha\in\mathbb{R}^{N-L+1}$ such that

$\begin{bmatrix}H_{L}(u)\\\
H_{L}(\hat{y})\end{bmatrix}\alpha=\begin{bmatrix}\bar{u}\\\
\bar{y}\end{bmatrix}.$ (3)

This lemma states that in the absence of disturbances, i.e., if $w_{k}=0$ for
all $k\geq 0$, all trajectories of system (1) can be parametrized by linear
combinations of time shifts of a priori measured, sufficiently exciting
input/output trajectories. In the following two sections, we set up MPC
schemes that use these trajectories for the prediction of the systems
behavior.

##  III Data-driven state-feedback predictive control

In this section, we present a robust data-driven state-feedback MPC scheme
with closed-loop guarantees on stability and constraint satisfaction in the
presence of process noise. In Subsection III-A, we introduce the data-driven
MPC scheme for the case of available state measurements. Thereafter, in
Subsection III-B we prove the closed-loop guarantees of the introduced control
scheme. Finally, in Subsection III-C we show how the system constants, which
are used to set up the constraint tightening of the MPC scheme, can be
approximated purely from data.

###  III-A Proposed MPC scheme

For the first data-driven predictive control scheme, we consider the
availability of full state measurement, i.e., $C=I$, $D=0$ in (1). Moreover,
we assume that the process disturbance belongs to the hypercube
$w_{t}\in\mathbb{W}=\left\\{w\in\mathbb{R}^{n}\mid\left\lVert
w\right\rVert_{\infty}\leq w_{\mathrm{max}}\right\\}$ for all $t\geq 0$, where
$w_{\mathrm{max}}\geq 0$ is known, and the input and state constraint sets are
given by the hypercubes
$u_{t}\in\mathbb{U}=\left\\{u\in\mathbb{R}^{m}\mid\left\lVert
u\right\rVert_{\infty}\leq u_{\mathrm{max}}\right\\}$ and
$x_{t}\in\mathbb{X}=\left\\{x\in\mathbb{R}^{n}\mid\left\lVert
x\right\rVert_{\infty}\leq x_{\mathrm{max}}\right\\}$ for some
$u_{\mathrm{max}}>0$, $x_{\mathrm{max}}>0$. As will become clear later in this
section, it is crucial for the construction of a proper constraint tightening
that the prediction model is stable. If this is not the case a priori (i.e.,
$A$ is unstable), then this can be enforced via a pre-stabilizing input
parametrization

$u_{k}=Kx_{k}+\nu_{k},$ (4)

as it is common, for example, in tube-based MPC[13, 14]. The state feedback
matrix $K$ is chosen such that all eigenvalues of $A_{K}=A+BK$ strictly lie
inside the unit disc. Such a pre-stabilizing controller can be computed purely
from data, e.g., following the approaches by Berberich et al.[22] or van
Waarde et al.[23]. Throughout this paper, we assume that such a controller is
known a priori. In case of a stable system, the following scheme can be
applied with $K=0$.

To make use of Lemma 1 for the prediction of open-loop state sequences, we
consider the input $\nu_{k}$ of the pre-stabilized system

$x_{k+1}=A_{K}x_{k}+B\nu_{k}+w_{k}.$ (5)

We now apply a p.e. input sequence $\\{\nu_{k}^{d}\\}_{k=0}^{N-1}$ of length
$N$ to System (5), and measure the associated disturbed state sequence
$\\{x^{d}_{k}\\}_{k=0}^{N}$, where the superscript ”$d$” denotes a priori
collected data.

######  Assumption 1.

The input sequence $\\{\nu_{k}^{d}\\}_{k=0}^{N-1}$ is persistently exciting of
order $L+n+1$.

We denote the cumulated disturbance influencing the collected data at time $k$
as

$d^{d}_{k}=\sum_{i=0}^{k-1}A^{k-1-i}_{K}w_{i}^{d}$ (6)

and the undisturbed state at as

$\hat{x}_{k}^{d}=x_{k}^{d}-d^{d}_{k}.$ (7)

######  Remark 1.

Note that the following results can be easily extended to the case where,
apart from the process disturbance $w_{k}$, also additive measurement noise on
the state measurements is present. In this case, the disturbance sequence in
(6) has to be extended by the actual measurement noise instant occurring at
time $k$ and, in case of $K\neq 0$, the past measurement noise which are fed
back into and propagated through the system dynamics. However, for the sake of
simplicity, we consider only process disturbances throughout the paper.

With these a priori generated data sequences, we are now able to set up the
following optimal control problem (OCP), given the measured state $x_{t}$ at
time $t$ and with the prediction horizon $L$

$\displaystyle J_{L}^{\ast}(x_{t})=$
$\displaystyle\min_{\begin{subarray}{c}\alpha(t),\sigma(t),\\\
\bar{\nu}(t),\bar{x}(t)\end{subarray}}$
$\displaystyle\sum_{k=0}^{L-1}\left(\left\lVert\bar{\nu}_{k}(t)\right\rVert_{R}^{2}+\left\lVert\bar{x}_{k}(t)\right\rVert_{Q}^{2}\right)+\lambda_{\alpha}w_{\mathrm{max}}\left\lVert\alpha(t)\right\rVert_{2}^{2}+\frac{\lambda_{\sigma}}{w_{\mathrm{max}}}\left\lVert\sigma(t)\right\rVert_{2}^{2}$
(8a) s.t. $\displaystyle\begin{bmatrix}\bar{\nu}(t)\\\
\bar{x}(t)+\sigma(t)\end{bmatrix}=\begin{bmatrix}H_{L}(\nu^{d})\\\
H_{L+1}(x^{d})\end{bmatrix}\alpha(t),$ (8b)
$\displaystyle\bar{x}_{0}(t)=x_{t},$ (8c) $\displaystyle\bar{x}_{L}(t)=0,$
(8d)
$\displaystyle\left\lVert\bar{x}_{k}(t)\right\rVert_{\infty}+a_{u,k}\left\lVert\bar{\nu}(t)\right\rVert_{1}+a_{\alpha,k}\left\lVert\alpha(t)\right\rVert_{1}+a_{\sigma,k}\left\lVert\sigma_{k}(t)\right\rVert_{\infty}+a_{c,k}\leq
x_{\mathrm{max}},$ (8e)
$\displaystyle\left\lVert\bar{\nu}_{k}(t)\right\rVert_{\infty}+b_{u,k}\left\lVert\bar{\nu}(t)\right\rVert_{1}+b_{\alpha,k}\left\lVert\alpha(t)\right\rVert_{1}+b_{\sigma,k}\left\lVert\sigma_{k}(t)\right\rVert_{\infty}+b_{c,k}+\left\lVert
K\bar{x}_{k}(t)\right\rVert_{\infty}\leq u_{\mathrm{max}},$ (8f)
$\displaystyle\forall k=0,\dots,L-1.$ (8g)

We denote the optimal solution of (8) at time $t$ by $\bar{\nu}^{\ast}(t)$,
$\bar{x}^{\ast}(t)$, $\alpha^{\ast}(t)$, $\sigma^{\ast}(t)$. In (8b), we make
use of Lemma 1 for the prediction of future state sequences of the system.
Note that, a Hankel matrix of depth $L+1$ is used in the second block row of
(8b), since the predicted state sequence contains $L+1$ elements (from $k=0$
to $k=L$), whereas the predicted input sequence only contains $L$ elements
(from $k=0$ to k=$L-1$). Moreover, note that, as it is common in predictive
control based on Willems’ lemma, we make use of a slack variable $\sigma(t)$
(first introduced by Coulson et al.[3]) that renders (8b) feasible, even in
the presence of disturbances. The slack variable $\sigma$ as well as the
variable $\alpha$ are regularized in (8a). This leads to smaller values of
$\sigma$ and $\alpha$, improving the prediction accuracy and reducing the
influence of disturbances in the Hankel matrices. For further discussion on
these issues, see also Section IV.A by Berberich et al.[4] and Section IV by
Dörfler et al. [24]. Problem (8) contains the tightened state and input
constraints (8e) and (8f). Note that both constraints depend on $\bar{u}(t)$,
$\alpha(t)$ and $\sigma(t)$. Together with suitably defined coefficients
$a_{u,k}$, $a_{\alpha,k}$, $a_{\sigma,k}$, $a_{c,k}$, and $b_{u,k}$,
$b_{\alpha,k}$, $b_{\sigma,k}$, $b_{c,k}$, which will be defined later on,
this constrained tightening ensures recursive feasibility and closed-loop
constraint satisfaction (cf. Theorem 1), i.e., $\left\lVert
x_{t}\right\rVert_{\infty}\leq x_{\mathrm{max}}$ and $\left\lVert
u_{t}\right\rVert_{\infty}=\left\lVert Kx_{t}+\nu_{t}\right\rVert_{\infty}\leq
u_{\mathrm{max}}$ for all $t\geq 0$. Constraint (8f) can be dropped if no pre-
stabilizing controller is used, i.e., $K=0$. Problem (8) is similar to the one
by Berberich et al.[18], which, however, does not consider process
disturbances and the input constraint tightening (8f). The predictive control
scheme is used in an $n$-step receding horizon manner, i.e., at time $t$ we
solve (8) and choose $\nu_{t+k}=\bar{\nu}_{k}^{\ast}(t)$ in (4) for
$k=0,\dots,n-1$.

In the following, we introduce the coefficients used to set up the tightened
state and input constraints (8e) and (8f). To this end, we first introduce
some system constants. We denote the disturbance sequence
$\\{w_{i}\\}_{i=0}^{k}$ propagated $k$ steps through the system dynamics by

$d_{k}\coloneqq\sum_{i=0}^{k-1}A^{k-1-i}_{K}w_{i},$ (9)

and an upper bound on its $\infty$-norm as

$\bar{d}_{k}\geq\sum_{i=0}^{k-1}\left\lVert
A^{k-1-i}_{K}\right\rVert_{\infty}w_{\mathrm{max}}\geq\left\lVert
d_{k}\right\rVert_{\infty}.$ (10)

Moreover, we define the constant $c_{pe}=\left\lVert
H_{u\hat{x}}^{\dagger}\right\rVert_{1}$, with

$H_{u\hat{x}}=\begin{bmatrix}H_{L}\left(\nu^{d}\right)\\\
H_{1}\left(\hat{x}^{d}_{[0,N-L-1]}\right)\end{bmatrix},$ (11)

where $H_{u\hat{x}}^{\dagger}$ is the Moore-Penrose inverse of $H_{u\hat{x}}$.
Using these system constants, as well as $\rho_{A,k}\geq\left\lVert
A_{K}^{k}\right\rVert_{\infty}$, we define

$c_{\alpha,k}=\rho_{A,k}\bar{d}_{N-L}+\bar{d}_{N-L+k},\quad
c_{\sigma,k}=\rho_{A,k}+1,$ (12)

for $k=0,\dots,L$. Moreover, we introduce the controllability constant
$\Gamma>0$, which is chosen such that, starting at any $x_{0}$, we can find an
input sequence $\nu_{[0,n-1]}$ steering the state of the pre-stabilized system
(5) to the origin in $n$ steps and satisfying

$\left\lVert\nu_{[0,n-1]}\right\rVert_{1}\leq\Gamma\left\lVert
x_{0}\right\rVert_{\infty}.$ (13)

Note that such a constant exists as the pair $\left(A,B\right)$ is
controllable.

We are now ready to define the coefficients of the state and input constraint
tightening as

$\displaystyle\begin{split}a_{u,k}&=0,\ a_{\alpha,k}=c_{\alpha,k},\
a_{\sigma,k}=c_{\sigma,k},\ a_{c,k}=\bar{d}_{k},\\\ b_{u,k}&=0,\
b_{\alpha,k}=\bar{K}c_{\alpha,k},\ b_{\sigma,k}=\bar{K}c_{\sigma,k},\
b_{c,k}=\bar{K}\bar{d}_{k},\end{split}$ (14)

for $k=0,\dots,n-1$, and

$\displaystyle\begin{split}a_{u,k+n}&=a_{u,k}+a_{\alpha,k}c_{pe}+a_{\sigma,k}c_{pe}\bar{d}_{N-1},\\\
a_{\alpha,k+n}&=a_{u,k+n}\Gamma c_{\alpha,L-1}+c_{\alpha,k+n},\\\
a_{\sigma,k+n}&=a_{u,k+n}\Gamma c_{\sigma,L-1}+c_{\sigma,k+n},\\\
a_{c,k+n}&=a_{c,k}+a_{\alpha,k}c_{pe}\left(nx_{\mathrm{max}}+n\bar{d}_{n}\right)+a_{\sigma,k}\left(\bar{d}_{N-1}c_{pe}\left(nx_{\mathrm{max}}+n\bar{d}_{n}\right)+\bar{d}_{n}\right)+\bar{d}_{n},\\\
b_{u,k+n}&=b_{u,k}+b_{\alpha,k}c_{pe}+b_{\sigma,k}c_{pe}\bar{d}_{N-1},\\\
b_{\alpha,k+n}&=b_{u,k+n}\Gamma c_{\alpha,L-1}+\bar{K}c_{\alpha,k+n},\\\
b_{\sigma,k+n}&=b_{u,k+n}\Gamma c_{\sigma,L-1}+\bar{K}c_{\sigma,k+n},\\\
b_{c,k+n}&=b_{c,k}+b_{\alpha,k}c_{pe}\left(nx_{\mathrm{max}}+n\bar{d}_{n}\right)+b_{\sigma,k}\left(\bar{d}_{N-1}c_{pe}\left(nx_{\mathrm{max}}+n\bar{d}_{n}\right)+\bar{d}_{n}\right)+\bar{K}\bar{d}_{n},\end{split}$
(15)

for $k=0,\dots,L-n-1$, where $\bar{K}=\left\lVert K\right\rVert_{\infty}$.
Note that $\bar{d}_{k}$ and $\rho_{A,k}$ grow exponentially if $A_{K}$ has
eigenvalues outside the unit disc. This is the main motivation for the usage
of the pre-stabilizing controller (4), as diverging $\bar{d}_{k}$ and
$\rho_{A,k}$ would also lead to diverging $c_{\alpha,k}$ and $c_{\sigma,k}$
and, therefore, to large coefficients (14) and (15). This would in general
yield an infeasible OCP (8) even for small prediction horizons $L$. In order
to set up the coefficients above, the system constants $\Gamma$, $c_{pe}$,
$\rho_{A,k}$ for $k=0,\dots,L$, and $d_{k}$ for $k=0,\dots,N-1$ have to be
known. All of these constants can be approximated from data as will be shown
in Subsection III-C.

###  III-B Theoretical guarantees

Firstly, we denote the undisturbed state at time $t+k$ resulting from an open-
loop application of $\bar{\nu}^{\ast}(t)$ as

$\hat{x}_{t+k}^{\ast}\coloneqq
A_{K}^{k}x_{t}+\sum_{i=0}^{k-1}A_{K}^{k-1-i}B\bar{\nu}_{i}^{\ast}(t).$ (16)

An upper bound for the prediction error between this undisturbed open-loop
state trajectory $\hat{x}^{\ast}$ and the predicted optimal state sequence
$\bar{x}^{\ast}(t)$ at time $t$ can be derived by the following lemma.

######  Lemma 2.

If (8) is feasible at time t, then

$\left\lVert\hat{x}_{t+k}^{\ast}-\bar{x}^{\ast}_{k}(t)\right\rVert_{\infty}\leq
c_{\alpha,k}\left\lVert\alpha^{\ast}(t)\right\rVert_{1}+c_{\sigma,k}\left\lVert\sigma^{\ast}(t)\right\rVert_{\infty}$
(17)

holds for all $k=0,\dots,L$.

###### Proof.

Similar to the proof of Lemma 2 in the work of Berberich et al.[4], we start
by considering the error between the undisturbed open-loop state and the state
prediction resulting from undisturbed data in the Hankel matrix

$x^{-}_{[t,t+L]}\coloneqq\hat{x}_{[t,t+L]}^{\ast}-H_{L+1}\left(\hat{x}^{d}\right)\alpha^{\ast}(t)=\hat{x}_{[t,t+L]}^{\ast}-\left(H_{L+1}\left(x^{d}\right)-H_{L+1}\left(d^{d}\right)\right)\alpha^{\ast}(t).$
(18)

Note that $\hat{x}_{[t,t+L-1]}^{\ast}$ and
$H_{L}\left(\hat{x}^{d}\right)\alpha^{\ast}(t)$ are trajectories of the
undisturbed system (5) resulting from the application of $\bar{\nu}^{\ast}(t)$
with different initial states. Due to (8c) and (8b), the initial condition of
$H_{L+1}\left(\hat{x}^{d}\right)\alpha^{\ast}(t)$ is

$H_{1}\left(\hat{x}^{d}_{[0,N-L]}\right)\alpha^{\ast}(t)=x_{t}+\sigma_{0}(t)-H_{1}\left(\hat{d}^{d}_{[0,N-L]}\right)\alpha^{\ast}(t).$
(19)

Thus, for the difference between both trajectories, it holds that

$x^{-}_{t+k}=A_{K}^{k}\left(H_{1}\left(d^{d}_{[0,N-L]}\right)\alpha^{\ast}(t)-\sigma_{0}(t)\right)$
(20)

for $k=0,\dots,L$. To show (17), note that

$\displaystyle\left\lVert\hat{x}_{t+k}^{\ast}-\bar{x}^{\ast}_{k}(t)\right\rVert_{\infty}$
$\displaystyle\leq\left\lVert
x^{-}_{t+k}\right\rVert_{\infty}+\left\lVert\sigma_{k}^{\ast}(t)-H_{1}\left(d^{d}_{[k,N-L+k]}\right)\alpha^{\ast}(t)\right\rVert_{\infty}$
(21) $\displaystyle\leq\left(\left\lVert
A_{K}^{k}\right\rVert_{\infty}\bar{d}_{N-L}+\bar{d}_{N-L-1}\right)\left\lVert\alpha^{\ast}(t)\right\rVert_{1}+\left\lVert
A_{K}^{k}\right\rVert_{\infty}\left\lVert\sigma_{0}^{\ast}(t)\right\rVert_{\infty}+\left\lVert\sigma_{k}^{\ast}(t)\right\rVert_{\infty},$
(22)

where the second inequality holds due to

$\displaystyle\left\lVert
H_{1}\left(d^{d}_{[k,N-L+k]}\right)\alpha^{\ast}(t)\right\rVert_{\infty}$
$\displaystyle\leq\left\lVert
d^{d}_{[k,N-L+k]}\right\rVert_{\infty}\left\lVert\alpha^{\ast}(t)\right\rVert_{1}$
(23)
$\displaystyle\leq\bar{d}_{N-L+k}\left\lVert\alpha^{\ast}(t)\right\rVert_{1}$
(24)

for $k=0,\dots,L$. Therefore, we obtain (17) with
$c_{\alpha,k}=\rho_{A,k}\bar{d}_{N-L}+\bar{d}_{N-L+k}$,
$c_{\sigma,k}=\rho_{A,k}+1$, and $\rho_{A,k}\geq\left\lVert
A_{K}^{k}\right\rVert_{\infty}$. ∎

Using the result of Lemma 2, we can now state our main result, which
establishes recursive feasibility, practical exponential stability, and input
and state constraint satisfaction of the closed-loop system, assuming that the
initial state is feasible for Problem (8) and the disturbance bound is
sufficiently small.

######  Theorem 1.

Suppose that Assumption 1 holds. Then, for any $V_{\mathrm{ROA}}>0$, there
exist $\underline{\lambda}_{\alpha}$, $\overline{\lambda}_{\alpha}$,
$\underline{\lambda}_{\sigma}$, $\overline{\lambda}_{\sigma}$ such that for
all $\lambda_{\alpha}$, $\lambda_{\sigma}$ satisfying

$\underline{\lambda}_{\alpha}\leq\lambda_{\alpha}\leq\overline{\lambda}_{\alpha},\quad\underline{\lambda}_{\sigma}\leq\lambda_{\sigma}\leq\overline{\lambda}_{\sigma},$
(25)

there exist $\bar{w}_{\mathrm{max}}$, $\bar{c}_{pe}>0$ as well as a
continuous, strictly increasing function
$\beta:[0,\bar{w}_{\mathrm{max}}]\to[0,V_{\mathrm{ROA}}]$ with $\beta(0)=0$,
such that for all $w_{\mathrm{max}}$ and $c_{pe}$ satisfying

$w_{\mathrm{max}}\leq\min{\left\\{\bar{w}_{\mathrm{max}},\;\frac{\bar{c}_{pe}}{c_{pe}}\right\\}}$
(26)

the following holds for the closed loop resulting from an application of the
$n$-step MPC scheme:

  1. (i)

If $J_{L}^{\ast}(x_{0})\leq V_{\mathrm{ROA}}$, then OCP (8) is feasible at any
time $t\geq 0$.

  2. (ii)

For any initial condition satisfying $J_{T}^{\ast}(x_{0})\leq
V_{\mathrm{ROA}}$ it holds that $x_{t}\in\mathbb{X}$ and $u_{t}\in\mathbb{U}$
for all $t\geq 0$, and $J_{L}^{\ast}(x_{t})$ converges exponentially to
$J_{L}^{\ast}(x_{t})\leq\beta(\bar{w}_{\mathrm{max}})$.

The proof of this result is similar to the one of Theorem 10 by Berberich et
al.[18]. The main differences lie in the fact that we consider process instead
of measurement noise; moreover, we use a pre-stabilizing controller, and thus,
include the input constraint tightening (8f), which also has to be satisfied
by the considered candidate solution, into the OCP. Furthermore, note that
(ii) only shows exponential convergence of $x_{t}$ to a neighborhood of $x=0$,
however, it is possible to establish a suitable lower as well as an upper
bound on $J_{L}^{\ast}(x_{t})$ analogous to Lemma 1 by Berberich et al.[4],
thus, resulting in practical exponential stability. For a detailed discussion
on the influence of the parameters $\lambda_{\alpha}$, $\lambda_{\sigma}$, and
$c_{pe}$ on the stability properties, we refer to the work of Berberich et
al.[4]. In short, the region of attraction increases for smaller disturbance
bounds or better persistency of excitation of the input signal, the latter
being expressed by a decrease of $c_{pe}$.

###### Proof.

To show (i), we construct a candidate solution for (8) at time $t+n$ and show
that (8b)-(8f) hold for this candidate. Therefore, we define the candidate
solution over the first $L-n$ steps via the previously optimal input shifted
by $n$ steps, i.e., $\bar{\nu}^{\prime}_{k}(t+n)=\bar{\nu}^{\ast}_{k+n}(t)$
for $k=0,\dots,L-n-1$. Moreover, the state candidate is chosen as

$\bar{x}^{\prime}_{[0,L-n]}(t+n)=\begin{bmatrix}x_{t+n}\\\
\hat{x}^{\ast}_{[t+n+1,t+L]}\end{bmatrix}$ (27)

for the first $L-n+1$ steps, with $\hat{x}^{\ast}$ defined as in (16). Note
that, $c_{\alpha,k}$ scales linearly with $w_{\mathrm{max}}$ and

$\left\lVert\alpha^{\ast}(t)\right\rVert_{1}\leq\sqrt{N-L+1}\left\lVert\alpha^{\ast}(t)\right\rVert_{2}\leq\sqrt{\left(N-L+1\right)\frac{V_{\mathrm{ROA}}}{\lambda_{\alpha}w_{\mathrm{max}}}},$
(28)

which implies that the first term on the right hand side of (17) becomes
arbitrarily small for sufficiently small $w_{\mathrm{max}}$. Furthermore, the
same holds true for the second term on the right hand side of (17), since
$c_{\sigma,k}$ is uniformly upper bounded for all $k$ and

$\left\lVert\sigma^{\ast}(t)\right\rVert_{\infty}\leq\left\lVert\sigma(t)\right\rVert_{2}\leq\sqrt{\frac{V_{\mathrm{ROA}}w_{\mathrm{max}}}{\lambda_{\sigma}}}.$
(29)

Due to this, Lemma 2, and $\bar{x}^{\ast}_{L}(t)=0$, the state
$\bar{x}^{\prime}_{L-n}(t+n)=\hat{x}^{\ast}_{t+L}$ becomes arbitrarily small
for sufficiently small $w_{\mathrm{max}}$. Thus, for a sufficiently small
$w_{\mathrm{max}}$, by controllability, there exists an input sequence
$\bar{\nu}^{\prime}_{[L-n,L-1]}(t+n)$ that steers
$\bar{x}^{\prime}_{[L-n,L]}(t+n)$ from $\hat{x}^{\ast}_{t+L}$ to $0$ in $n$
steps. Furthermore, we choose the candidate solution for
$\alpha^{\prime}(t+n)$ and $\sigma^{\prime}(t+n)$ as

$\alpha^{\prime}(t+n)=H_{u\hat{x}}^{\dagger}\begin{bmatrix}\bar{\nu}^{\prime}(t+n)\\\
\hat{x}_{t+n}^{\ast}\end{bmatrix},$ (30)

and

$\begin{split}\sigma^{\prime}(t+n)&=H_{L+1}\left(x^{d}\right)\alpha^{\prime}(t+n)-\bar{x}^{\prime}(t+n),\\\
&=\begin{bmatrix}\hat{d}_{t+n}\\\
H_{L}\left(d_{[1,N]}^{d}\right)\alpha^{\prime}(t+n)\end{bmatrix},\end{split}$
(31)

where $\hat{d}_{t+k}$ satisfies $x_{t+k}=\hat{x}_{t+k}^{\ast}+\hat{d}_{t+k}$.
Thus, the candidate solution satisfies (8b)-(8d).

It remains to be shown that also (8e) and (8f) hold for the candidate
solution. In order to show the satisfaction of (8f), we note that

$\begin{split}\left\lVert\alpha^{\prime}(t+n)\right\rVert_{1}&\stackrel{{\scriptstyle\eqref{eq:alpha_candidate}}}{{\leq}}c_{pe}\left(\left\lVert\bar{\nu}^{\prime}(t+n)\right\rVert_{1}+\left\lVert\hat{x}_{t+n}^{\ast}\right\rVert_{1}\right)\\\
&\leq
c_{pe}\left(\left\lVert\bar{\nu}^{\prime}(t+n)\right\rVert_{1}+nx_{\mathrm{max}}+n\bar{d}_{n}\right)\end{split}$
(32)

holds. Furthermore, it holds that

$\begin{split}\left\lVert\sigma^{\prime}(t+n)\right\rVert_{\infty}&\stackrel{{\scriptstyle\eqref{eq:sigma_candidate}}}{{\leq}}\bar{d}_{N-1}\left\lVert\alpha^{\prime}(t+n)\right\rVert_{1}+\bar{d}_{n},\\\
&\stackrel{{\scriptstyle\eqref{eq:alpha_bound}}}{{\leq}}\bar{d}_{N-1}c_{pe}\left(\left\lVert\bar{\nu}^{\prime}(t+n)\right\rVert_{1}+nx_{\mathrm{max}}+n\bar{d}_{n}\right)+\bar{d}_{n}.\end{split}$
(33)

Using the same arguments as in Proposition 8 (Inequality (18)) by Berberich et
al.[18], we can bound the norm of the candidate input by

$\left\lVert\bar{\nu}^{\prime}(t+n)\right\rVert_{1}\leq\left\lVert\bar{\nu}^{\ast}(t)\right\rVert_{1}+\Gamma
c_{\alpha,L-1}\left\lVert\alpha^{\ast}(t)\right\rVert_{1}+\Gamma
c_{\sigma,L-1}\left\lVert\sigma^{\ast}(t)\right\rVert_{\infty}.$ (34)

Since (8f) holds for the optimal solution at time $t$ and due to
$\left\lVert\bar{\nu}_{k}^{\prime}(t+n)\right\rVert_{\infty}=\left\lVert\bar{\nu}^{\ast}_{k+n}\right\rVert_{\infty}$
for $k=0,\dots,L-n-1$, we obtain

$\left\lVert\bar{\nu}_{k}^{\prime}(t+n)\right\rVert_{\infty}\leq
u_{\mathrm{max}}-\left(b_{u,k+n}\left\lVert\bar{\nu}^{\ast}(t)\right\rVert_{1}+b_{\alpha,k+n}\left\lVert\alpha^{\ast}(t)\right\rVert_{1}+b_{\sigma,k+n}\left\lVert\sigma^{\ast}(t)\right\rVert_{\infty}+b_{c,k+n}+\left\lVert
K\bar{x}^{\ast}_{k+n}(t)\right\rVert_{\infty}\right).$ (35)

Plugging (34) and the coefficients in (15) into (35) yields for
$k=0,\dots,L-n-1$

$\displaystyle\begin{split}\left\lVert\bar{\nu}_{k}^{\prime}(t+n)\right\rVert_{\infty}&\stackrel{{\scriptstyle\eqref{eq:input_bound}}}{{\leq}}u_{\mathrm{max}}-b_{u,k}\left\lVert\bar{\nu}^{\prime}(t+n)\right\rVert_{1}-b_{\alpha,k}\left(c_{pe}\left\lVert\bar{\nu}^{\prime}(t+n)\right\rVert_{1}+c_{pe}(nx_{\mathrm{max}}+n\bar{d}_{n})\right)\\\
&\quad-
b_{\sigma,k}\left(\bar{d}_{N-1}c_{pe}\left\lVert\bar{\nu}^{\prime}(t+n)\right\rVert_{1}+\bar{d}_{N-1}c_{pe}(nx_{\mathrm{max}}+n\bar{d}_{n})+\bar{d}_{n}\right)-b_{c,k}\\\
&\quad-\bar{K}c_{\alpha,k+n}\left\lVert\alpha^{\ast}(t)\right\rVert_{1}-\bar{K}c_{\sigma,k+n}\left\lVert\sigma^{\ast}(t)\right\rVert_{\infty}-\left\lVert
K\bar{x}^{\ast}_{k+n}(t)\right\rVert_{\infty}-\bar{K}\bar{d}_{n},\end{split}$
(36)
$\displaystyle\begin{split}&\stackrel{{\scriptstyle\eqref{eq:alpha_bound},\eqref{eq:sigma_bound}}}{{\leq}}u_{\mathrm{max}}-b_{u,k}\left\lVert\bar{\nu}^{\prime}(t+n)\right\rVert_{1}-b_{\alpha,k}\left\lVert\alpha^{\prime}(t+n)\right\rVert_{1}-b_{\sigma,k}\left\lVert\sigma^{\prime}(t+n)\right\rVert_{\infty}-b_{c,k}\\\
&\quad-\bar{K}c_{\alpha,k+n}\left\lVert\alpha^{\ast}(t)\right\rVert_{1}-\bar{K}c_{\sigma,k+n}\left\lVert\sigma^{\ast}(t)\right\rVert_{\infty}-\left\lVert
K\bar{x}^{\ast}_{k+n}(t)\right\rVert_{\infty}-\bar{K}\bar{d}_{n},\end{split}$
(37)
$\displaystyle\begin{split}&\stackrel{{\scriptstyle\eqref{eq:prediction_error_sf}}}{{\leq}}u_{\mathrm{max}}-b_{u,k}\left\lVert\bar{\nu}^{\prime}(t+n)\right\rVert_{1}-b_{\alpha,k}\left\lVert\alpha^{\prime}(t+n)\right\rVert_{1}-b_{\sigma,k}\left\lVert\sigma^{\prime}(t+n)\right\rVert_{\infty}-b_{c,k}\\\
&\quad-\bar{K}\left\lVert\hat{x}_{t+k+n}^{\ast}-\bar{x}^{\ast}_{k+n}(t)\right\rVert_{\infty}-\left\lVert
K\bar{x}^{\ast}_{k+n}(t)\right\rVert_{\infty}-\bar{K}\bar{d}_{n},\end{split}$
(38) $\displaystyle\begin{split}&\leq
u_{\mathrm{max}}-b_{u,k}\left\lVert\bar{\nu}^{\prime}(t+n)\right\rVert_{1}-b_{\alpha,k}\left\lVert\alpha^{\prime}(t+n)\right\rVert_{1}-b_{\sigma,k}\left\lVert\sigma^{\prime}(t+n)\right\rVert_{\infty}-b_{c,k}-\left\lVert
K\bar{x}^{\prime}_{k}(t+n)\right\rVert_{\infty},\end{split}$ (39)

where the last inequality holds due to
$\bar{x}_{k}^{\prime}(t+n)=\hat{x}_{t+n+k}^{\ast}$ for $k=1,\dots,L-n-1$, and

$\left\lVert\bar{x}^{\prime}_{0}(t+n)-\hat{x}_{t+n}^{\ast}\right\rVert_{\infty}=\left\lVert
x_{t+n}-\hat{x}_{t+n}^{\ast}\right\rVert_{\infty}\leq\bar{d}_{n}.$ (40)

Therefore, the candidate solution satisfies (8f) for $k=0,\dots,L-n-1$.
Showing that also (8e) holds for $k=0,\dots,L-n-1$ can be done following the
analogous steps as above.

To show that (8e) and (8f) are also satisfied for $k=L-n,\dots,L-1$, we recall
from above that $\bar{x}^{\prime}_{L-n}(t+n)$ becomes arbitrarily small for
sufficiently small $w_{\mathrm{max}}$. Thus, due to controllability and (13)
also $\bar{\nu}_{[L-n,L-1]}^{\prime}(t+n)$ and, therefore,
$\bar{x}_{[L-n+1,L-1]}^{\prime}(t+n)$ become arbitrarily small. Moreover, due
to (12) and (15) the coefficients $a_{\alpha,k}$ and $b_{\alpha,k}$ depend
linearly on $w_{\mathrm{max}}$, and thus, they also become arbitrarily small
for sufficiently small disturbance bounds $w_{\mathrm{max}}$. Moreover, the
same holds for $a_{u,k}$, $a_{c,k}$, $b_{u,k}$, and $b_{c,k}$. The
coefficients $a_{\sigma,k}$, $b_{\sigma,k}$ converge to constant values for
$w_{\mathrm{max}}\to 0$.Finally, due to (33), $\sigma^{\prime}(t+n)$ becomes
arbitrarily small if $w_{\mathrm{max}}$ is sufficiently small. Hence, all
terms on the left-hand side of (8e) and (8f) (except for $\bar{x}_{k}(t)$ and
$\bar{v}_{k}(t)$, respectively) become arbitrary small. Therefore, since
$x_{\mathrm{max}},\ u_{\mathrm{max}}>0$, there exists a sufficiently small
bound $\bar{w}_{\mathrm{max}}>0$ such that (8e) and (8f) are also satisfied
for $k=L-n,\dots,L-1$.

To show (ii), we can follow the same arguments as in Theorem 3 by Berberich et
al.[4] to conclude invariance of the sublevel set $J_{L}^{\ast}(x_{t})\leq
V_{\mathrm{ROA}}$ and exponential convergence of $J_{L}^{\ast}(x_{t})$ to
$J_{L}^{\ast}(x_{t})\leq\beta(\bar{w}_{\mathrm{max}})$. This is possible since
the candidate solution used in the first part of our proof is constructed
analogous to the one used in the proof of the above reference. Thus, the
closed-loop scheme is recursively feasible. Closed-loop state constraint
satisfaction follows from the same arguments used in the proof of Theorem 10
by Berberich et al.[18]. In order to show closed-loop input constraint
satisfaction, we note that the optimal solution at time $t$ satisfies (8f) for
the first $n$ steps, i.e.,

$\displaystyle u_{\mathrm{max}}$
$\displaystyle\geq\left\lVert\bar{\nu}^{\ast}_{k}(t)\right\rVert_{\infty}+b_{u,k}\left\lVert\bar{\nu}^{\ast}(t)\right\rVert_{1}+b_{\alpha,k}\left\lVert\alpha^{\ast}(t)\right\rVert_{1}+b_{\sigma,k}\left\lVert\sigma^{\ast}(t)\right\rVert_{\infty}+b_{c,k}+\left\lVert
K\bar{x}^{\ast}_{k}(t)\right\rVert_{\infty}$ (41)
$\displaystyle=\left\lVert\bar{\nu}^{\ast}_{k}(t)\right\rVert_{\infty}+\bar{K}c_{\alpha,k}\left\lVert\alpha^{\ast}(t)\right\rVert_{1}+\bar{K}c_{\sigma,k}\left\lVert\sigma^{\ast}(t)\right\rVert_{\infty}+\bar{K}\bar{d}_{k}+\left\lVert
K\bar{x}_{k}^{\ast}(t)\right\rVert_{\infty},$ (42)

for $k=0,\dots,n-1$, where the second inequality follows from plugging in the
coefficients from (14). With (17), we finally obtain for $k=0,\dots,n-1$

$\displaystyle u_{\mathrm{max}}$
$\displaystyle\geq\left\lVert\bar{\nu}^{\ast}_{k}(t)\right\rVert_{\infty}+\bar{K}\left\lVert\hat{x}_{t+k}^{\ast}-\bar{x}^{\ast}_{k}(t)\right\rVert_{\infty}+\bar{K}\bar{d}_{k}+\left\lVert
K\bar{x}_{k}^{\ast}(t)\right\rVert_{\infty}$
$\displaystyle\geq\left\lVert\bar{\nu}^{\ast}_{k}(t)\right\rVert_{\infty}+\bar{K}\left\lVert\hat{x}_{t+k}^{\ast}-\bar{x}^{\ast}_{k}(t)\right\rVert_{\infty}+\bar{K}\left\lVert
x_{t+k}-\hat{x}_{t+k}^{\ast}\right\rVert_{\infty}+\left\lVert
K\bar{x}_{k}^{\ast}(t)\right\rVert_{\infty}$
$\displaystyle\geq\left\lVert\bar{\nu}^{\ast}_{k}(t)+Kx_{t+k}\right\rVert_{\infty},$

where the second inequality holds due to $\left\lVert
x_{t+k}-\hat{x}^{\ast}_{t+k}\right\rVert_{\infty}\leq\bar{d}_{k}$ and the
third inequality due to the triangle inequality. Thus, the input constraints
are satisfied in closed loop. ∎

###  III-C Data-driven Estimation of System Constants

In order to set up the constraint tightening, we need to compute the
coefficients (14) and (15). These depend on several system constants. While
$\bar{K}$, $x_{\mathrm{max}}$, and $n$ are known a priori, the system
constants $\rho_{A,k}$, $\bar{d}_{k}$, $c_{pe}$ and $\Gamma$ have to be
estimated from data. In the following, we provide corresponding estimation
procedures.

An approach for the computation of $\Gamma$ was shown by Berberich et al.[18],
where, however, the availability of the undisturbed data is assumed. To extend
this approach for the approximation of $\Gamma$ in the presence of
disturbances in the data, we adapt the respective optimization problem by
including a slack variable as well as a regularization, similar to (8). Thus,
we set up the optimization problem

$\displaystyle\max_{x_{0}}\min_{\begin{subarray}{c}\alpha,\sigma,\\\
\bar{\nu},\bar{x}\end{subarray}}$
$\displaystyle\quad\left\lVert\bar{\nu}_{[0,n-1]}\right\rVert_{1}+\lambda_{\alpha}^{\prime}w_{\mathrm{max}}\left\lVert\alpha\right\rVert_{2}^{2}+\frac{\lambda_{\sigma}^{\prime}}{w_{\mathrm{max}}}\left\lVert\sigma\right\rVert_{2}^{2},$
(43a) s.t. $\displaystyle\quad\left\lVert\bar{x}_{0}\right\rVert_{\infty}\leq
x_{\mathrm{max}},$ (43b) $\displaystyle\quad\bar{x}_{n}=0,$ (43c)
$\displaystyle\quad\begin{bmatrix}\bar{\nu}\\\
\bar{x}+\sigma\end{bmatrix}=\begin{bmatrix}H_{n}(\nu^{d})\\\
H_{n+1}(x^{d})\end{bmatrix}\alpha,$ (43d)

which can be solved by solving the inner minimization problem for all vertices
of $\mathbb{X}$. Then, we choose
$\Gamma\approx\frac{\left\lVert\bar{\nu}^{\ast}_{[0,n-1]}\right\rVert_{1}}{x_{\mathrm{max}}}$,
where $\bar{\nu}^{\ast}_{[0,n-1]}$ is the optimal solution of (43). Moreover,
we approximate $c_{pe}\approx\left\lVert H_{ux}^{\dagger}\right\rVert_{1}$,
with

$H_{ux}=\begin{bmatrix}H_{L}\left(\nu^{d}\right)\\\
H_{1}\left(x^{d}_{[0,N-L-1]}\right)\end{bmatrix}.$ (44)

This approximation is possible, as for small disturbances also the error
between $H_{ux}$ and $H_{u\hat{x}}$ is small. So far, both procedures
mentioned above only yield approximations of the real constants $\Gamma$ and
$c_{pe}$, without being guaranteed overapproximations of these constants.
However, as was also confirmed in our numerical examples, the error between
the real values and our estimates remains small for the considered disturbance
level. Obtaining guaranteed overapproximations of the constants $\Gamma$ and
$c_{pe}$ even in the presence of noise is an interesting subject for future
research.

In order to estimate the overapproximations $\rho_{A,k}\geq\left\lVert
A_{K}^{k}\right\rVert_{\infty}$ and
$\bar{d}_{k}\geq\sum_{i=0}^{k-1}\left\lVert
A_{K}^{k-1-i}\right\rVert_{\infty}w_{\mathrm{max}}$ for $i=0,\dots,L-1$ we use
Algorithm 1 by Wildhagen et al.[25], which makes use of the S-lemma[23]. Note
that the setup in the aforementioned reference considers a bound on the
$2$-norm of the disturbance. However, by noting that $\left\lVert
w\right\rVert_{2}\leq\sqrt{n}\left\lVert w\right\rVert_{\infty}$, we can
easily adapt the algorithm to our setting.

##  IV Data-driven output-feedback predictive control

In this section, we construct a robust data-driven predictive control scheme,
in case only output measurements are available. First, in Subsection IV-A we
set up the MPC scheme and prove similar theoretical properties to the ones
shown for the state-feedback case. Thereafter, in Subsection IV-B we show how
the coefficients used for the constraint tightening can be computed only from
data.

###  IV-A Proposed MPC scheme and theoretical guarantees

In contrast to the previous section, we consider the case where no full state
measurements are available, but only output measurements. To this end, we
consider the difference operator form

$y_{k}=-A_{n}y_{k-1}-\ldots-
A_{2}y_{k-n+1}-A_{1}y_{k-n}+Du_{k}+B_{n}u_{k-1}+\ldots+B_{2}u_{k-n+1}+B_{1}u_{k-n}+\tilde{w}_{k},$
(45)

with the process disturbance
$\tilde{w}_{k}\in\mathbb{\tilde{W}}=\left\\{\tilde{w}\in\mathbb{R}^{p}\mid\left\lVert\tilde{w}\right\rVert_{\infty}\leq\tilde{w}_{\mathrm{max}}\right\\}$.
The input and output constraints are given by
$u_{t}\in\mathbb{U}=\left\\{u\in\mathbb{R}^{m}\mid\left\lVert
u\right\rVert_{\infty}\leq u_{\mathrm{max}}\right\\}$ and
$y_{t}\in\mathbb{Y}=\left\\{y\in\mathbb{R}^{p}\mid\left\lVert
y\right\rVert_{\infty}\leq y_{\mathrm{max}}\right\\}$ for some
$u_{\mathrm{max}}>0$, $y_{\mathrm{max}}>0$, similar to the setup in Section
III. Note that the following results also hold if only an upper bound on the
system order is known, in which case $n$ needs to be replaced by this upper
bound. Furthermore, note that (45) is an equivalent characterization of the
input-output behavior of (1), with $C\neq I$ or $D\neq 0$. Moreover, we can
transform (45) into the non-minimal realization

$\begin{split}\xi_{k+1}&=\tilde{A}\xi_{k}+\tilde{B}u_{k}+\tilde{E}\tilde{w}_{k},\\\
\tilde{y}_{k}&=\tilde{C}\xi_{k}+Du_{k}+\tilde{w}_{k},\end{split}$ (46)

with the extended state
$\xi_{k}=\begin{bmatrix}u_{k-n}^{\top}&\dots&u_{k-1}^{\top}&y_{k-n}^{\top}&\dots&y_{k-1}^{\top}\end{bmatrix}^{\top}$,
cf. [26]. Similar to the state feedback case in Section III, we want to make
use of a pre-stabilizing controller in case of an unstable system. Therefore,
we introduce the control law

$u_{k}=\tilde{K}\xi_{k}+\nu_{k},$ (47)

where the stabilizing feedback matrix $\tilde{K}$ can be computed purely from
data, e.g., following the approach by Berberich et al.[27]. However, for
simplicity it is assumed, that such a pre-stabilizing controller is known a
priori. Thus, $\nu_{k}$ is the input to the stabilized system

$\begin{split}\xi_{k+1}&=\tilde{A}_{K}\xi_{k}+\tilde{B}\nu_{k}+\tilde{E}\tilde{w}_{k},\\\
y_{k}&=\tilde{C}_{K}\xi_{k}+D\nu_{k}+\tilde{w}_{k}.\end{split}$ (48)

Similar to the previous section, we define the disturbance propagated $k$
steps through the system dynamics as

$\delta_{k}\coloneqq\sum_{i=0}^{k-1}\tilde{A}^{k-1-i}_{K}\tilde{w}_{i},$ (49)

and an upper bound on its $\infty$-norm as

$\bar{\delta}_{k}=\sum_{i=0}^{k-1}\eta_{A}^{k-1-i}\tilde{w}_{\mathrm{max}}\geq\sum_{i=0}^{k-1}\left\lVert\tilde{A}^{k-1-i}_{K}\right\rVert_{\infty}\tilde{w}_{\mathrm{max}}\geq\left\lVert\delta_{k}\right\rVert_{\infty},$
(50)

with $\eta_{A}\geq\left\lVert\tilde{A}_{K}\right\rVert_{\infty}$.

Again, we apply a p.e. input sequence $\\{\nu_{k}^{d}\\}_{k=0}^{N-1}$ of
length $N$ to System (48), and measure the associated disturbed output
sequence $\\{y^{d}_{k}\\}_{k=0}^{N-1}$. As the OCP, introduced in the
following, now contains $n$ additional steps to fix the initial state, the
following assumption is needed.

######  Assumption 2.

The input sequence $\\{\nu_{k}^{d}\\}_{k=0}^{N-1}$ is persistently exciting of
order $L+2n$.

Using these a priori collected data sequences as well as the constants
$\eta_{A}$ from above,
$\eta_{B}\geq\left\lVert\tilde{B}\right\rVert_{\infty}$,
$\eta_{C}\geq\left\lVert\tilde{C}_{K}\right\rVert_{\infty}$,
$\eta_{D}\geq\left\lVert D\right\rVert_{\infty}$, we set up the OCP for the
output-feedback predictive control problem as

$\displaystyle
J_{L}^{\ast}(\bar{\nu}_{[0,n-1]}(t-n),y_{\left[t-n,t-1\right]})=$ (51a)
$\displaystyle\min_{\begin{subarray}{c}\alpha(t),\sigma(t),\\\
\bar{\nu}(t),\bar{y}(t)\end{subarray}}$
$\displaystyle\sum_{k=0}^{L-1}\left(\left\lVert\bar{\nu}_{k}(t)\right\rVert_{R}^{2}+\left\lVert\bar{y}_{k}(t)\right\rVert_{Q}^{2}\right)+\lambda_{\alpha}\tilde{w}_{\mathrm{max}}\left\lVert\alpha(t)\right\rVert_{2}^{2}+\frac{\lambda_{\sigma}}{\tilde{w}_{\mathrm{max}}}\left\lVert\sigma(t)\right\rVert_{2}^{2}$
(51b) s.t. $\displaystyle\begin{bmatrix}\bar{\nu}(t)\\\
\bar{y}(t)+\sigma(t)\end{bmatrix}=\begin{bmatrix}H_{L+n}(\nu^{d})\\\
H_{L+n}(y^{d})\end{bmatrix}\alpha(t),$ (51c)
$\displaystyle\begin{bmatrix}\bar{\nu}_{\left[-n,-1\right]}(t)\\\
\bar{y}_{\left[-n,-1\right]}(t)\end{bmatrix}=\begin{bmatrix}\bar{\nu}^{\ast}_{\left[0,n-1\right]}(t-n)\\\
y_{\left[t-n,t-1\right]}\end{bmatrix},$ (51d)
$\displaystyle\begin{bmatrix}\bar{\nu}_{\left[L-n,L-1\right]}(t)\\\
\bar{y}_{\left[L-n,L-1\right]}(t)\end{bmatrix}=\begin{bmatrix}0\\\
0\end{bmatrix},$ (51e)
$\displaystyle\left\lVert\tilde{K}\right\rVert_{\infty}\eta_{A}^{k}\left\lVert\xi_{t}\right\rVert_{\infty}+\left\lVert\tilde{K}\right\rVert_{\infty}\sum_{i=0}^{k-1}\eta_{A}^{k-1-i}\eta_{B}\left\lVert\bar{\nu}_{i}(t)\right\rVert_{\infty}+\left\lVert\bar{\nu}_{k}(t)\right\rVert_{\infty}+\left\lVert\tilde{K}\right\rVert_{\infty}\bar{\delta}_{k}\leq
u_{\mathrm{max}},$ (51f)
$\displaystyle\eta_{C}\eta_{A}^{k}\left\lVert\xi_{t}\right\rVert_{\infty}+\eta_{C}\sum_{i=0}^{k-1}\eta_{A}^{k-1-i}\eta_{B}\left\lVert\bar{\nu}_{i}(t)\right\rVert_{\infty}+\eta_{D}\left\lVert\bar{\nu}_{k}(t)\right\rVert_{\infty}+\eta_{C}\bar{\delta}_{k}\leq
y_{\mathrm{max}},$ (51g) $\displaystyle\forall k=0,\dots,L-n-1.$ (51h)

Problem (51) is similar to (8) in the state-feedback case. The key difference
is the new input and output constraint tightening in (51f) and (51g). These
constraints are now independent of $\alpha$ and $\sigma$. Instead, the
tightening only involves the extended state at time $t$ and the input
variables as well as the above defined constants. However, this comes at the
price of potential conservatism, since (51f) and (51g) involve terms
$\eta_{A}^{k}$, which is in general larger than terms of the form
$\left\lVert\tilde{A}_{K}^{k}\right\rVert_{\infty}$ that have been used in the
constraint tightening of Section III. Setting up a constraint tightening
similar to the one of (8) remains an interesting issue for future research.
Moreover, the initial constraint (51d) and the terminal constraint (51e) hold
over $n$ steps. This implies that the internal state of the underlying minimal
realization corresponding to the prediction coincides with the initial state
and the terminal state, respectively (compare the work of Markovsky and
Rapisarda[28]). Note that the constants in (51f) and (51g) can be computed
purely from data following the approach that will be discussed in the next
subsection.

We, again, close the loop by applying the optimal solution of (51) in an
$n$-step manner, i.e., $u_{t+k}=\tilde{K}\xi_{t+k}+\bar{\nu}_{k}^{\ast}(t)$
for $k=0,\dots,n-1$, where $\bar{\nu}_{k}^{\ast}(t)$ is the optimal solution
of (51) for the prediction step $k$. Note that $\xi_{t+k}$ contains the inputs
$u_{[t+k-n,t+k-1]}$ and the measured outputs $y_{[t+k-n,t+k-1]}$. We are now
ready to state practical exponential stability, and input as well as output
constraint satisfaction of the closed loop. To this end, following the
approach by Berberich et al.[4], we now consider the Lyapunov function

$V_{t}\coloneqq
J_{L}^{\ast}(u_{\left[t-n,t-1\right]},y_{\left[t-n,t-1\right]})+\gamma
W(\xi_{t}),$ (52)

for some $\gamma>0$, where $W(\xi)$ is an IOSS Lyapunov function, which exists
due to detectability of $(A,C)$[29].

######  Theorem 2.

Suppose that Assumption 2 holds. Then, for any $V_{\mathrm{ROA}}>0$, there
exist $\underline{\lambda}_{\alpha}$, $\overline{\lambda}_{\alpha}$,
$\underline{\lambda}_{\sigma}$, $\overline{\lambda}_{\sigma}$ such that for
all $\lambda_{\alpha}$, $\lambda_{\sigma}$ satisfying

$\underline{\lambda}_{\alpha}\leq\lambda_{\alpha}\leq\overline{\lambda}_{\alpha},\quad\underline{\lambda}_{\sigma}\leq\lambda_{\sigma}\leq\overline{\lambda}_{\sigma},$
(53)

there exist $\bar{\tilde{w}}_{\mathrm{max}}$, $\bar{c}_{pe}>0$ as well as a
continuous, strictly increasing function
$\beta:[0,\bar{\tilde{w}}_{\mathrm{max}}]\to[0,V_{\mathrm{ROA}}]$ with
$\beta(0)=0$, such that for all $\tilde{w}_{\mathrm{max}}$ and $c_{pe}$
satisfying

$\tilde{w}_{\mathrm{max}}\leq\min{\left\\{\bar{\tilde{w}}_{\mathrm{max}},\;\frac{\bar{c}_{pe}}{c_{pe}}\right\\}}$
(54)

the following holds for the closed loop resulting from the application of the
$n$-step MPC scheme:

  1. (i)

If $V_{0}\leq V_{\mathrm{ROA}}$, then OCP (51) is feasible at any time $t\geq
0$.

  2. (ii)

For any initial condition satisfying $V_{0}\leq V_{\mathrm{ROA}}$ it holds
that $y_{t}\in\mathbb{Y}$ and $u_{t}\in\mathbb{U}$ for all $t\geq 0$, and
$V_{t}$ converges exponentially to
$V_{t}\leq\beta(\bar{\tilde{w}}_{\mathrm{max}})$.

This result is similar to the state-feedback case (Theorem 1). Also the proof
works along the lines of the proof of Theorem 1, where the candidate solution
can be chosen analogously. The main difference lies in the modified constraint
tightening. For a discussion on the role of the parameters in the above
statement, we refer to the discussion below Theorem 1.

###### Proof.

The proof is analogous to the proof of Theorem 1. The only difference lies in
the constraint tightening (51f), (51g). To show (i), we note that

$\bar{\nu}^{\prime}(t+n)=\begin{bmatrix}\bar{\nu}^{\ast}_{[0,L-n-1]}(t)\\\
\bar{\nu}^{\prime}_{[L-2n,L-n-1]}(t+n)\\\ 0_{n}\end{bmatrix}$ (55)

and

$\bar{y}^{\prime}(t+n)=\begin{bmatrix}y_{[t-n,t-1]}\\\
\hat{y}_{[t+n,t+L-n-1]}^{\ast}\\\ \bar{y}^{\prime}_{[L-2n,L-n-1]}(t+n)\\\
0_{n}\end{bmatrix}$ (56)

are an input and output candidate solution. Here, $\hat{y}^{\ast}$ is the
undisturbed output starting at $x_{t}$, resulting from an open-loop
application of $\bar{\nu}^{\ast}(t)$; furthermore,
$\bar{\nu}^{\prime}_{[L-2n,L-n-1]}(t+n)$ is the input steering the system to
the origin in $n$ steps and $\bar{y}^{\prime}_{[L-2n,L-n-1]}(t+n)$ is the
associated output. Note that, analogous to the proof of Theorem 1, such an
input exists due to controllability of the system. Moreover, we choose the
candidate

$\alpha^{\prime}(t+n)=H_{u\hat{x}}^{y\dagger}\begin{bmatrix}\bar{\nu}^{\prime}(t+n)\\\
x_{t}\end{bmatrix},$ (57)

with

$H_{u\hat{x}}^{y}=\begin{bmatrix}H_{L+n}\left(\nu^{d}\right)\\\
H_{1}\left(\hat{x}^{d}_{[0,N-L-n]}\right)\end{bmatrix}.$ (58)

Furthermore, as a candidate for the slack variable, we choose

$\sigma^{\prime}(t+n)=H_{L+n}\left(y^{d}\right)\alpha^{\prime}(t+n)-\bar{y}^{\prime}(t+n).$
(59)

Thus, the candidate solution satisfies (51c)-(51e). To show that (51f) is also
satisfied for the candidate solution, we note that
$\bar{\nu}^{\prime}_{k}(t+n)=\bar{\nu}^{\ast}_{k+n}(t)$ holds for
$k=0,\dots,L-2n-1$. Thus, it holds for these $k$ that

$\displaystyle u_{\mathrm{max}}$
$\displaystyle\geq\left\lVert\tilde{K}\right\rVert_{\infty}\eta_{A}^{k+n}\left\lVert\xi_{t}\right\rVert_{\infty}+\left\lVert\tilde{K}\right\rVert_{\infty}\sum_{i=0}^{k+n-1}\eta_{A}^{k+n-1-i}\eta_{B}\left\lVert\bar{\nu}_{i}^{\ast}(t)\right\rVert_{\infty}+\left\lVert\bar{\nu}_{k+n}^{\ast}(t)\right\rVert_{\infty}+\left\lVert\tilde{K}\right\rVert_{\infty}\bar{\delta}_{k+n},$
(60)
$\displaystyle\begin{split}&\geq\left\lVert\tilde{K}\right\rVert_{\infty}\eta_{A}^{k}\left(\eta_{A}^{n}\left\lVert\xi_{t}\right\rVert_{\infty}+\sum_{i=0}^{n-1}\eta_{A}^{n-1-i}\eta_{B}\left\lVert\bar{\nu}_{i}^{\ast}(t)\right\rVert_{\infty}+\bar{\delta}_{n}\right)+\left\lVert\tilde{K}\right\rVert_{\infty}\sum_{i=0}^{k-1}\eta_{A}^{k-1-i}\eta_{B}\left\lVert\bar{\nu}_{i}^{\prime}(t+n)\right\rVert_{\infty}\\\
&\quad+\left\lVert\bar{\nu}_{k}^{\prime}(t+n)\right\rVert_{\infty}+\left\lVert\tilde{K}\right\rVert_{\infty}\bar{\delta}_{k},\end{split}$
(61)
$\displaystyle\geq\left\lVert\tilde{K}\right\rVert_{\infty}\eta_{A}^{k}\left\lVert\xi_{t+n}\right\rVert_{\infty}+\left\lVert\tilde{K}\right\rVert_{\infty}\sum_{i=0}^{k-1}\eta_{A}^{k-1-i}\eta_{B}\left\lVert\bar{\nu}_{i}^{\prime}(t+n)\right\rVert_{\infty}+\left\lVert\tilde{K}\right\rVert_{\infty}\bar{\delta}_{k}+\left\lVert\bar{\nu}_{k}^{\prime}(t+n)\right\rVert_{\infty},$
(62)

where the last inequality holds due to
$\xi_{t+n}=\tilde{A}_{K}^{n}+\sum_{i=0}^{n-1}\tilde{A}_{K}^{n-1-i}\tilde{B}\nu_{i}^{\ast}(t)+\sum_{i=0}^{n-1}\tilde{A}_{K}^{n-1-i}\tilde{w}_{t+i}$.
Therefore, the candidate input satisfies (51f) for $k=0,\dots,L-2n-1$. Showing
that for these prediction steps (51g) is satisfied by the candidate output
works analogously. Moreover, showing that (51f) and (51g) are also satisfied
for $k=L-2n,\dots,L-n-1$ can be done following analogous steps to the ones in
the proof of Theorem 1, by noting that due to the terminal condition and
controllability there exists a sufficiently small $\tilde{w}_{\mathrm{max}}$
such that $\bar{\nu}^{\prime}_{k}(t+n)$ and $\bar{y}^{\prime}_{k}(t+n)$ become
arbitrarily small, thus satisfying the tightened input and output constraint.

To show (ii), we again follow the same arguments as in Theorem 3 by Berberich
et al.[4] to conclude practical exponential stability and recursive
feasibility. It remains to show closed-loop input and output constraint
satisfaction. Therefore, we note that

$\displaystyle u_{\mathrm{max}}$
$\displaystyle\geq\left\lVert\tilde{K}\right\rVert_{\infty}\eta_{A}^{k}\left\lVert\xi_{t}\right\rVert_{\infty}+\left\lVert\tilde{K}\right\rVert_{\infty}\sum_{i=0}^{k-1}\eta_{A}^{k-1-i}\eta_{B}\left\lVert\bar{\nu}_{i}^{\ast}(t)\right\rVert_{\infty}+\left\lVert\tilde{K}\right\rVert_{\infty}\bar{\delta}_{k}+\left\lVert\bar{\nu}_{k}^{\ast}(t)\right\rVert_{\infty}$
(63)
$\displaystyle\geq\left\lVert\tilde{K}\left(\tilde{A}_{K}^{k}\xi_{t}+\sum_{i=0}^{k-1}\tilde{A}_{K}^{k-1-i}\tilde{B}\bar{\nu}_{i}^{\ast}(t)+\delta_{t+k}\right)+\bar{\nu}_{k}^{\ast}(t)\right\rVert_{\infty}$
(64)
$\displaystyle=\left\lVert\tilde{K}\xi_{t+k}+\nu_{t+k}\right\rVert_{\infty}$
(65)

holds for $k=0,\dots,n-1$, which proves closed-loop input constraint
satisfaction. Showing output constraint satisfaction works analogously. ∎

###  IV-B Data-driven Estimation of System Constants

In the following, we provide data-based procedures to compute
(overapproximations of) the coefficients $\eta_{A}$, $\eta_{B}$, $\eta_{C}$,
and $\eta_{D}$. First, we note that $\eta_{A}$ can be computed based on[25],
analogously to $\eta_{A,k}$ in Subsection III-C. The coefficients $\eta_{B}$,
$\eta_{C}$, and $\eta_{D}$ can be computed similarly by modifying the approach
from[25].

Let us consider the data matrices

$\begin{split}X_{+}&\coloneqq\begin{bmatrix}\xi^{d}_{n+1}&\xi^{d}_{n+2}&\dots&\xi^{d}_{N}\end{bmatrix},\\\
X&\coloneqq\begin{bmatrix}\xi^{d}_{n}&\xi^{d}_{n+1}&\dots&\xi^{d}_{N-1}\end{bmatrix},\\\
Y&\coloneqq\begin{bmatrix}y^{d}_{n}y^{d}_{n+1}&\dots&y^{d}_{N-1}\end{bmatrix},\\\
U&\coloneqq\begin{bmatrix}\nu^{d}_{n}\nu^{d}_{n+1}&\dots&\nu^{d}_{N-1}\end{bmatrix},\end{split}$
(66)

where
$\xi^{d}_{k}=\begin{bmatrix}u_{k-n}^{d\top}&\dots&u_{k-1}^{d\top}&y_{k-n}^{d\top}&\dots&y_{k-1}^{d\top}\end{bmatrix}^{\top}$
for $k=n,\dots,N$. We write $F\succeq 0$ if $F$ is a symmetric and positive
semi-definite matrix. In order to compute overapproximations of the respective
system constants, we solve the optimization problem

$\displaystyle\min_{\tau,\bar{\sigma}^{2}}$
$\displaystyle\quad\bar{\sigma}^{2},$ (67a) s.t. $\displaystyle\quad
P_{1}(\bar{\sigma}^{2})-\tau P_{2}\succeq 0,$ (67b)
$\displaystyle\quad\bar{\sigma}^{2}\geq 0,$ (67c) $\displaystyle\quad\tau\geq
0,$ (67d)

where $P_{1}$ and $P_{2}$ are placeholders, which have to be defined for the
different coefficients as follows

$\displaystyle\eta_{A}:\quad P_{1}(\bar{\sigma}^{2})$
$\displaystyle=\begin{bmatrix}-I&0&0\\\ 0&0&0\\\
0&0&\bar{\sigma}^{2}I\end{bmatrix},$ $\displaystyle P_{2}$
$\displaystyle=\begin{bmatrix}-\begin{bmatrix}X\\\
U\end{bmatrix}^{\top}&X_{+}^{\top}\\\
0&I\end{bmatrix}^{\top}\cdot\begin{bmatrix}-\tilde{E}\tilde{E}^{\top}&0\\\
0&n\tilde{w}_{\mathrm{max}}^{2}N\end{bmatrix}\cdot\begin{bmatrix}-\begin{bmatrix}X\\\
U\end{bmatrix}^{\top}&X_{+}^{\top}\\\ 0&I\end{bmatrix},$ (68)
$\displaystyle\eta_{B}:\quad P_{1}(\bar{\sigma}^{2})$
$\displaystyle=\begin{bmatrix}0&0&0\\\ 0&-I&0\\\
0&0&\bar{\sigma}^{2}I\end{bmatrix},$ $\displaystyle P_{2}$
$\displaystyle=\begin{bmatrix}-\begin{bmatrix}X\\\
U\end{bmatrix}^{\top}&X_{+}^{\top}\\\
0&I\end{bmatrix}^{\top}\cdot\begin{bmatrix}-\tilde{E}\tilde{E}^{\top}&0\\\
0&n\tilde{w}_{\mathrm{max}}^{2}N\end{bmatrix}\cdot\begin{bmatrix}-\begin{bmatrix}X\\\
U\end{bmatrix}^{\top}&X_{+}^{\top}\\\ 0&I\end{bmatrix},$ (69)
$\displaystyle\eta_{C}:\quad P_{1}(\bar{\sigma}^{2})$
$\displaystyle=\begin{bmatrix}-I&0&0\\\ 0&0&0\\\
0&0&\bar{\sigma}^{2}I\end{bmatrix},$ $\displaystyle P_{2}$
$\displaystyle=\begin{bmatrix}-\begin{bmatrix}X\\\
U\end{bmatrix}^{\top}&Y^{\top}\\\
0&I\end{bmatrix}^{\top}\cdot\begin{bmatrix}-I^{\top}&0\\\
0&n\tilde{w}_{\mathrm{max}}^{2}N\end{bmatrix}\cdot\begin{bmatrix}-\begin{bmatrix}X\\\
U\end{bmatrix}^{\top}&Y^{\top}\\\ 0&I\end{bmatrix},$ (70)
$\displaystyle\eta_{D}:\quad P_{1}(\bar{\sigma}^{2})$
$\displaystyle=\begin{bmatrix}0&0&0\\\ 0&-I&0\\\
0&0&\bar{\sigma}^{2}I\end{bmatrix},$ $\displaystyle P_{2}$
$\displaystyle=\begin{bmatrix}-\begin{bmatrix}X\\\
U\end{bmatrix}^{\top}&Y^{\top}\\\
0&I\end{bmatrix}^{\top}\cdot\begin{bmatrix}-I^{\top}&0\\\
0&n\tilde{w}_{\mathrm{max}}^{2}N\end{bmatrix}\cdot\begin{bmatrix}-\begin{bmatrix}X\\\
U\end{bmatrix}^{\top}&Y^{\top}\\\ 0&I\end{bmatrix}.$ (71)

We denote the solutions of (67) as $\bar{\sigma}_{A}^{2}$,
$\bar{\sigma}_{B}^{2}$, $\bar{\sigma}_{C}^{2}$, $\bar{\sigma}_{D}^{2}$ for the
respective configuration of $P_{1}$ and $P_{2}$. It is straight forward to
show that we obtain $\eta_{A}$, $\eta_{B}$, $\eta_{C}$, and $\eta_{D}$ via
$\eta\leq\sqrt{n}\bar{\sigma}$. For a detailed discussion of the approach we
refer to[25].

##  V Numerical example

Figure 1: Coefficients for the state and input constraint tightening (8e),
(8f). The red lines correspond to the ideal coefficients that can be computed
if perfect model knowledge is available. The blue lines correspond to the
coefficients computed purely from data. Figure 2: Ideal coefficients for the
state constraint tightening (8e) without pre-stabilizing controller Figure 3:
Simulation

As an example, we consider the two mass-spring-system suggested by Wie et
al.[30], with the masses $m_{1}=0.5\,\mathrm{kg}$, $m_{2}=1\,\mathrm{kg}$ and
the spring constant $k=2\,\frac{\mathrm{kg}}{\mathrm{s}^{2}}$. Discretizing
the system with a sampling time of $1\,\mathrm{s}$ yields the matrices

$A=\begin{bmatrix}-0.1799&1.1799&0.507&0.493\\\ 0.59&0.41&0.2465&0.7535\\\
-1.0421&1.0421&-0.1799&1.1799\\\ 0.5211&-0.5211&0.59&0.41\\\
\end{bmatrix},\quad B=\begin{bmatrix}0.7266\\\ 0.1367\\\ 1.014\\\
0.493\end{bmatrix}.$ (72)

We assume that full state measurements are available and that
$w_{\mathrm{max}}=10^{-3}$, $u_{\mathrm{max}}=10$, $x_{\mathrm{max}}=10$ hold
for the constraint sets, where the process disturbance $w_{k}$ acting on the
system, during the data generation and in closed-loop operation, at time $k$
is sampled uniformly from $\mathbb{W}$. As the matrix $A$ has two eigenvalues
on the unit disc, the usage of a pre-stabilizing controller is expected to be
be advantageous in order to set up the MPC scheme introduced in Section III.
First, we collect apply a PE input sequence of length $50$ to the open-loop
system and measure the corresponding state sequence. Thereafter, we employ
Theorem 1 by Berberich et al.[27] to compute a robust linear quadratic
regulator for the system based on the available noisy data (using diagonal
multipliers to describe the disturbance bound, compare Equation (21) in[27]).
This yields the state-feedback gain

$K=\begin{bmatrix}0.4345&-0.8439&-0.3665&-0.6986\end{bmatrix}$

which serves as the pre-stabilizing controller, leading to all eigenvalues of
$A_{K}=A+BK$ being located strictly inside the unit disc.

The goal is to set up the OCP (8) with prediction horizon $L=12$, a total
amount of $N=50$ data points (of the pre-stabilized system) in the Hankel
matrices, and the parametrization $Q=I_{n}$, $R=1$,
$\lambda_{\alpha}=\lambda_{\sigma}=100$. To this end, we have to compute the
constants (14) and (15) from data. In order to do so, we apply an input
sequence of length $N^{\prime}=5000$, which is uniformly sampled from
$\mathbb{U}$, to the pre-stabilized system and collect the corresponding
$N^{\prime}$ state measurements. Using these data sequences, we compute
$\rho_{A,k}$ for $k=0,\dots,L-1$ and $\bar{d}_{k}$ for $k=0,\dots,N-1$
following the approach mentioned in Subsection III-C. Note that the estimation
of the system constants also works for a smaller amount of data, at the price
of a more conservative overapproximation of these constants. However, for a
good approximation of these constants, we need a much larger amount of data
than for the prediction via the Hankel matrices (i.e., $N^{\prime}\gg N$). We
now choose an input-state-sequence of length $N$ from the collected data (of
total length $N^{\prime}$), for which the input sequence is persistently
exciting of order $L+n+1$, to construct the Hankel matrices and approximate
the constant $c_{pe}$ as described in Subsection III-C. Moreover, we apply the
method from Subsection III-C, with
$\lambda_{\alpha}^{\prime}=\lambda_{\sigma}^{\prime}=1$, in order to compute
an approximation of the controllability constant $\Gamma$. With these
approximations of the system constants, we compute the coefficients in (14)
and (15). The resulting coefficients as well as the ideal coefficients that
would be computable if perfect model knowledge was available, can be found in
Figure 1, where the red lines correspond to the ideal coefficients, and the
blue lines to the coefficients computed from data. It can be seen that the
coefficients computed from data yield over approximations of the ”real”
coefficients. This is mainly due to the fact that the procedure in Subsection
III-C only yields overapproximations of the constants $\rho_{A,k}$ and
$\bar{d}_{k}$. The approximation of $c_{pe}$ by its disturbed counterpart is
accurate for the present example,compare the discussion in Subsection III-C.

Considering the parameter $b_{c,k}$, it can be seen that for larger $k$, this
coefficient already is close to $u_{\mathrm{max}}$. Even though there is a
non-monotonicity in $k$, which results from the usage of an $n$-step MPC
scheme and the corresponding recursive definition of the constants in (14) and
(15), it is clearly visible that $b_{c,k}$ tends to increase for larger $k$.
Thus, for larger prediction horizons or larger disturbance bounds
$w_{\mathrm{max}}$, this parameter would lead to $b_{c,k}>u_{\mathrm{max}}$,
which would render (8) infeasible due to (8f). The reason for this
conservatism in the constraint tightening lies in the fact that
submultiplicativity and the triangle inequality were exploited multiple times
in the proof of recursive feasibility and constraint satisfaction. Moreover,
the estimates for $\rho_{A,k}$ and $\bar{d}_{k}$ only yield overapproximations
of the parameters and for the sake of recursive feasibility,
$x_{\mathrm{max}}$ instead of $\left\lVert x_{t}\right\rVert_{\infty}$ is used
to define the coefficients (15).

As a motivation for the usage of a pre-stabilizing controller, the
coefficients for $K=0$ (i.e., without input constraint tightening) are plotted
in Figure 2. Note that even for $k=4$ it holds that $a_{c,k}\approx 233$ which
already exceeds $x_{\mathrm{max}}$ and would thus lead to an infeasible OCP
even for the smallest possible prediction horizon of the $n$-step scheme,
$L=4$.

The states and inputs resulting from the $n$-step scheme in closed loop
starting at $x_{0}=\begin{bmatrix}4&-4&0&0\end{bmatrix}^{\top}$ can be found
in Figure 3. It can be seen that the predictive control scheme works as
desired, meaning it stabilizes the states at the origin, while satisfying the
state and input constraints.

##  VI Conclusion

In this paper, we introduced a data-driven predictive control scheme relying
on predictions based on a priori measured data, structured in Hankel matrices.
The scheme is capable of stabilizing the origin of an LTI system, even in the
presence of process disturbances acting on the system state. To this end, a
constraint tightening is proposed, which can be set up using only a priori
measured data. The presented MPC scheme allows for the usage of a pre-
stabilizing controller and an associated input constraint tightening, which
enables the use also for a priori unstable systems. Closed-loop recursive
feasibility, practical exponential stability, and constraint satisfaction of
the control scheme is shown. Moreover, the MPC, initially introduced for
available state measurements, is extended to the case that only output
measurements are available. The numerical experiments illustrated the
applicability of the proposed approach and underlined the necessity to include
a pre-stabilization and corresponding input constraint tightening in order to
design a feasible controller. Interesting issues for future research include
the development of less conservative constraint tightenings, in particular in
the output-feedback case, as well as a data-based techniques for obtaining
estimates of $c_{pe}$ and the controllability constant $\Gamma$ from noisy
data, which are guaranteed overapproximations of the real system constants.

## References

  * [1] Jan C Willems, Paolo Rapisarda, Ivan Markovsky, and Bart LM De Moor.  A note on persistency of excitation.  Systems & Control Letters, 54(4):325–329, 2005.
  * [2] Hua Yang and Shaoyuan Li.  A data-driven predictive controller design based on reduced hankel matrix.  In 2015 10th Asian Control Conference (ASCC), pages 1–7. IEEE, 2015\.
  * [3] Jeremy Coulson, John Lygeros, and Florian Dörfler.  Data-enabled predictive control: In the shallows of the DeePC.  In 2019 18th European Control Conference (ECC), pages 307–312. IEEE, 2019.
  * [4] Julian Berberich, Johannes Köhler, Matthias A Müller, and Frank Allgöwer.  Data-driven model predictive control with stability and robustness guarantees.  IEEE Transactions on Automatic Control, 66(4):1702–1717, 2020.
  * [5] J. Coulson, J. Lygeros, and F. Dörfler.  Distributionally robust chance constrained data-enabled predictive control.  IEEE Trans. Automat. Control, 2021.  doi: 10.1109/TAC.2021.3097706.
  * [6] Linbin Huang, Jianzhe Zhen, John Lygeros, and Florian Dörfler.  Robust data-enabled predictive control: Tractable formulations and performance guarantees.  arXiv preprint arXiv:2105.07199, 2021.
  * [7] Mingzhou Yin, Andrea Iannelli, and Roy S Smith.  Maximum likelihood estimation in data-driven modeling and control.  IEEE Transactions on Automatic Control, 2021.
  * [8] Mingzhou Yin, Andrea Iannelli, and Roy S Smith.  Maximum likelihood signal matrix model for data-driven predictive control.  In Learning for Dynamics and Control, pages 1004–1014. PMLR, 2021\.
  * [9] Anton Xue and Nikolai Matni.  Data-driven system level synthesis.  In Learning for Dynamics and Control, pages 189–200. PMLR, 2021\.
  * [10] Luca Furieri, Baiwei Guo, Andrea Martin, and Giancarlo Ferrari-Trecate.  Near-optimal design of safe output feedback controllers from noisy data.  arXiv preprint arXiv:2105.10280, 2021.
  * [11] Julian Berberich, Johannes Köhler, Matthias A Müller, and Frank Allgöwer.  Linear tracking mpc for nonlinear systems part ii: The data-driven case.  arXiv preprint arXiv:2105.08567, 2021.
  * [12] Ivan Markovsky and Florian Dörfler.  Behavioral systems theory in data-driven analysis, signal processing, and control.  Annual Reviews in Control, 52:42–64, 2021.
  * [13] Basil Kouvaritakis and Mark Cannon.  Model Predictive Control.  Springer International Publishing, 2016.
  * [14] Luigi Chisci, John Anthony Rossiter, and Giovanni Zappa.  Systems with persistent disturbances: predictive control with restricted constraints.  Automatica, 37(7):1019–1028, 2001.
  * [15] David Q Mayne, María M Seron, and SV Raković.  Robust model predictive control of constrained linear systems with bounded disturbances.  Automatica, 41(2):219–224, 2005.
  * [16] Anil Aswani, Humberto Gonzalez, S Shankar Sastry, and Claire Tomlin.  Provably safe and robust learning-based model predictive control.  Automatica, 49(5):1216–1226, 2013.
  * [17] Enrico Terzi, Lorenzo Fagiano, Marcello Farina, and Riccardo Scattolini.  Learning-based predictive control for linear systems: A unitary approach.  Automatica, 108:108473, 2019.
  * [18] J. Berberich, J. Köhler, M. A. Müller, and F. Allgöwer.  Robust constraint satisfaction in data-driven MPC.  In Proc. 59th IEEE Conf. Decision and Control (CDC), pages 1260–1267, 2020.
  * [19] L. Huang, J. Coulson, J. Lygeros, and F. Dörfler.  Decentralized data-enabled predictive control for power system oscillation damping.  IEEE Trans. Control Systems Technology, 2021.  doi: 10.1109/TCST.2021.3088638.
  * [20] J. Umenberger.  Closed-loop data-enabled predictive control.  In Proc. American Control Conf. (ACC), pages 3349–3356. IEEE, 2021\.
  * [21] Wenjie Liu, Jian Sun, Gang Wang, Francesco Bullo, and Jie Chen.  Data-driven resilient predictive control under denial-of-service.  arXiv preprint arXiv:2110.12766, 2021.
  * [22] Julian Berberich, Anne Koch, Carsten W Scherer, and Frank Allgöwer.  Robust data-driven state-feedback design.  In 2020 American Control Conference (ACC), pages 1532–1538. IEEE, 2020.
  * [23] H. J. Van Waarde, M. K. Camlibel, and M. Mesbahi.  From noisy data to feedback controllers: non-conservative design via a matrix S-lemma.  IEEE Trans. Automat. Control, 2020.  doi: 10.1109/TAC.2020.3047577.
  * [24] Florian Dörfler, Jeremy Coulson, and Ivan Markovsky.  Bridging direct & indirect data-driven control formulations via regularizations and relaxations.  IEEE Transactions on Automatic Control, 2022.
  * [25] Stefan Wildhagen, Julian Berberich, Michael Hertneck, and Frank Allgöwer.  Data-driven estimation of the maximum sampling interval: analysis and controller design for discrete-time systems.  arXiv preprint arXiv:2101.00903, 2021.
  * [26] G. C. Goodwin and K. S. Sin.  Adaptive filtering prediction and control.  Courier Corporation, 2014.
  * [27] J. Berberich, C. W. Scherer, and F. Allgöwer.  Combining prior knowledge and data for robust controller design.  IEEE Trans. Automat. Control, 2020.  submitted, preprint: arXiv:2009.05253.
  * [28] I. Markovsky and P. Rapisarda.  Data-driven simulation and control.  Int. J. Control, 81(12):1946–1959, 2008.
  * [29] C. Cai and A. R. Teel.  Input–output-to-state stability for discrete-time systems.  Automatica, 44(2):326–336, 2008.
  * [30] Bong Wie and Dennis S Bernstein.  Benchmark problems for robust control design.  Journal of Guidance, Control, and Dynamics, 15(5):1057–1059, 1992\.


#  The Diophantine Equation

$\mathbf{\arctan\left(\dfrac{1}{x}\right)+\arctan\left(\dfrac{\ell}{y}\right)=\arctan\left(\dfrac{1}{k}\right)}$

Konstantine Zelator
Mathematics, Statistics, and Computer Science
212 Ben Franklin Hall
Bloomsburg University of Pennsylvania
400 East 2nd Street
Bloomsburg, PA 17815
USA
and
P.O. Box 4280
Pittsburgh, PA 15203
<EMAIL_ADDRESS>
e-mails: konstantine<EMAIL_ADDRESS>
##  1 Introduction

The subject matter of this work is the two-variable diophantine equation
$\arctan\left(\dfrac{1}{x}\right)+\arctan\left(\dfrac{\ell}{y}\right)=\arctan\left(\dfrac{1}{k}\right)$
for given positive integers $k$ and $\ell$, such that gcd $(\ell,k^{2}+1)=1$
(i.e., $\ell$ and $k^{2}+1$ are relatively prime). The main objective is to
determine all positive integer pairs $(x,y)$ which satisfy

$\left\\{\begin{array}[]{l}\arctan\left(\dfrac{1}{x}\right)+\arctan\left(\dfrac{\ell}{y}\right)=\arctan\left(\dfrac{1}{k}\right)\\\
x,y\in\mathbb{Z}^{+},\ {\rm gcd}(\ell,k^{2}+1)=1\ {\rm and}\\\ {\rm with\
gcd}\ (\ell,y)=1\ \left({\rm i.e.,}\ \ell\ {\rm and}\ y\ {\rm are}\right.\\\
\left.{\rm relatively\ prime}\right)\end{array}\right\\}$ (1)

This is done in Theorem 1, Section 4. As we will see, there are exacgtly $N$
distinct solutions to (1) where $N$ is the number of positive divisors of the
integer $k^{2}+1$. The $N$ pairs $(x,y)$, which are solutions to (1), are
expressed parametrically in terms of the positive divisors of $k^{2}+1$. Also,
note that when $\ell=1$, equation (1) is symmetric with respect to the two
variables $x$ and $y$. If $(a,b)$ is a solution, then so is $(b,a)$. The
motivating force behind this work is a recent article published in the journal
Mathematics and Computer Education (see [1]). The article, authored by Hasan
Unal, is entitled “Proof without words: an arctangent equality”. It consists
of four illustrations, a purely geometric proof of the equality,

$\arctan\left(\dfrac{1}{3}\right)+\arctan\left(\dfrac{1}{7}\right)=\arctan\left(\dfrac{1}{2}\right).$

From the point of view of (1), the last equality says that the pair $(3,7)$ is
a solution of (1), in the case $\ell=1$ and $k=2$.

According to Theorem 1, $(3,7)$ and $(7,3)$ are the only solutions to (1) for
$\ell=1$ and $k=2$.

This, then, is the other objective of this article. To generate more
arctangent type of equalities. This is done in Section 5, where a listing of
such equalities is offered; an immediate consequence of Theorem 1.

In Section 2, we list two trigonometric preliminaries: the well known identity
for the tangent of the sum of two angles and a couple of basic facts regarding
arctangent function.

In Section 3, we state two well known results from number theory: Euclid’s
lemma; and the formula that gives the number of positive divisors of a
positive integer. We use these in the proof of Theorem 1.

##  2 Trigonometric preliminaries

  1. (a)

If $\theta_{1}$ and $\theta_{2}$ are two angles measured in radians, such that
neither $\theta_{1}$ nor $\theta_{2}$, nor their sum $\theta_{1}+\theta_{2}$
is of the form $k\pi+\dfrac{\pi}{2},\ k$ and integer. Then,

$\tan(\theta_{1}+\theta_{2})=\dfrac{\tan\theta_{1}+\tan\theta_{2}}{1-\tan\theta_{1}\theta_{2}}$

  2. (b)

Let $f$ be the arctangent function, $f(x)=\arctan x$. Then,

    1. (i)

$\arctan 1=\dfrac{\pi}{4}$

    2. (ii)

$\left\\{\begin{array}[]{lr}0<\theta<\dfrac{\pi}{r}&\\\ {\rm and}&\\\
&\theta=\arctan
c\end{array}\right\\}\Leftrightarrow\left\\{\begin{array}[]{l}0<\theta<\dfrac{\pi}{4}\\\
\\\ 0<c=\tan\theta<1\end{array}\right\\}$

    3. (iii)

$\left\\{\begin{array}[]{ll}0<\theta<\dfrac{\pi}{2}\\\ {\rm and}&\\\
&\theta=\arctan
c\end{array}\right\\}\Leftrightarrow\left\\{\begin{array}[]{l}0<\theta<\dfrac{\pi}{2}\\\
\\\ 0<c=\tan\theta\end{array}\right\\}$

##  3 Number theory preliminaries

The following result is commonly known as Euclid’s lemma, and is of great
significance in number theory.

Result 1 (Euclid’s lemma): Let $a,b,c$ be positive integers such that $a$ is a
divisor of the product $bc$; and with $a$ also being relatively prime to $b$.
Then, $a$ is a divisor of $c$.

The next result provides a formula that gives the exact number of positive
divisors of a positive integer.

Result 2 (number of divisors formula) Let $n\geq 2$ be a positive integer, and
let $p_{1},\ldots,p_{t}$ in increasing order, be the distinct prime bases that
appear in the prime factorization of $n$, so that $n=p^{e_{1}}_{1},\ldots
p^{e_{t}}_{t}$, with the exponents $e_{1},\ldots,e_{t}$ being positive
integers. Also, let $N$ be the number of positive divisors of $n$. Then,

  1. (i)

$N=\overset{t}{\underset{i=1}{\Pi}}(e_{i}+1)\ldots(e_{1}+1)\ldots(e_{t}+1)$.

  2. (ii)

In particular, when $e_{1}=\ldots=e_{t}=1$ (i.e., when $n$ is squarefree)

$N=2^{t}$

Both of these two results can be easily found in number theory books and
texts. For example, see reference [2].

##  4 Theorem 1 and its proof

######  Theorem 1.

Let $k$ and $\ell$ be fixed or given positive integers such that
gcd$(\ell,k^{2}+1)=1$. Consider the diophantine equation (1).

  1. (a)

There are exactly $N$ distinct positive integer pairs $(x,y)$ which are
solutions to equation (1) where $N$ is the number of positive integer divisors
of the integer $k^{2}+1$. Specifically, if $(x,y)$ is a positive integer
solution of (1), then

$x=k+\ell\cdot\left(\dfrac{k^{2}+1}{d}\right)\ {\rm and}\ y=k\ell+d$ where $d$
is a positive integer divisor of $k^{2}+1$.

  2. (b)

If $k^{2}+1=p$, a prime number, then equation (1) has exactly two distinct
positive integer solutions. These are

$(x,y)=(k+\ell(k^{2}+1),\ k\ell+1),\ \ (k+\ell,\ k\ell+k^{2}+1).$

  3. (c)

If $k^{2}+1=p_{1}p_{2}$, a product of two distinct primes $p_{1}$ and $p_{2}$,
equation (1) has exactly four distinct positive integer solutions. These are,

$\begin{array}[]{ll}(x,y)=(k+\ell(k^{2}+1),k\ell+1),&(k+\ell,\
k\ell+k^{2}+1),\\\ \\\ (k+\ell p_{2}.\ k\ell+p_{1}),\ {\rm and}&{}\ (k+\ell
p_{1},\ k\ell+p_{2})\end{array}$

###### Proof.

First note that parts (b) and (c) are immediate consequences of part (a) and
Result 2. We omit the details. We prove part (a)

  1. (a)

Let $d$ be a positive integer divisor of $k^{2}+1$. We will show that the
positive integer pair,
$(x_{d},y_{d})=\left(k+\ell\cdot\left(\dfrac{k^{2}+1}{d}\right),\
k\ell+d\right)$ is a solution to (1). First note that $y_{d}=k\ell+d$, is
relatively prime to $\ell$. Indeed, if $y_{d}$ and $\ell$ had a prime factor
$q$ in common then $q$ would divide $y_{d}-k\ell=d$; and thus (since $d$ is a
divisor of $k^{2}+1$) $y_{d}-k\ell=d$, then $q$ would divide $k^{2}+1$
contrary to the hypothesis that gcd$(\ell,k^{2}+1)=1$. Thus,
gcd$(\ell,y_{d})=1$.

It is clear that since $k,\ell$ and $d$ are positive integers, we have
$x_{d}>1,y_{d}>1$ and $k\geq 1$. So,

$\left(0<\dfrac{1}{x_{d}}<1,\ 0<\dfrac{\ell}{y_{d}}<1,\ 0<\dfrac{1}{k}\leq
1\right).$ (2)

Let

$\theta_{1}=\arctan\left(\dfrac{1}{x_{d}}\right),\
\theta_{2}=\arctan\left(\dfrac{1}{y_{d}}\right),\
\theta=\arctan\left(\dfrac{1}{k}\right).$ (3)

Then, by (2), (3) and part (b) of the trigonometric preliminaries, we have

$\left\\{\begin{array}[]{l}0<\theta_{1}<\dfrac{\pi}{4},\
0<\theta_{2}<\dfrac{\pi}{4},\ 0<\theta\leq\dfrac{\pi}{4}\\\ \\\ {\rm and}\
0<\theta_{1}+\theta_{2}<\dfrac{\pi}{2},\ \tan\theta_{1}=\dfrac{1}{x_{d}},\
\tan\theta_{2}=\dfrac{\ell}{y_{d}},\
\tan\theta=\dfrac{1}{k}\end{array}\right\\}$ (4)

From (4) and part (a) of trigonometric preliminaries, it follows that

$\begin{array}[]{rcl}\tan(\theta_{1}+\theta_{2})&=&\dfrac{\frac{1}{x_{d}}+\frac{\ell}{y_{d}}}{1-\frac{1}{x_{d}}\cdot\frac{\ell}{y_{d}}};\\\
\\\ \tan(\theta_{1}+\theta_{2})&=&\dfrac{y_{d}+\ell
x_{d}}{x_{d}y_{d}-\ell};\\\ \\\
\tan(\theta_{1}+\theta_{2})&=&\dfrac{d\cdot(y_{d}+\ell
x_{d})}{dx_{d}y_{d}-d\ell}.\end{array}$ (5)

By (5) and the expressions for $x_{d}$ and $y_{d}$ (see beginning of the
proof) we get

$\begin{array}[]{rcl}\tan(\theta_{1}+\theta_{2})&=&\dfrac{d^{2}+k\ell d+k\ell
d+\ell^{2}\cdot(k^{2}+1)}{[dk+\ell(k^{2}+1)](k\ell+d)-d\ell};\\\ \\\
\tan(\theta_{1}+\theta_{2})&=&\dfrac{d^{2}+2k\ell
d+\ell^{2}\cdot(k^{2}+1)}{d\ell k^{2}+k\ell^{2}(k^{2}+1)+kd^{2}+\ell
dk^{2}+d\ell-d\ell};\\\ \\\ \tan(\theta_{1}+\theta_{2})&=&\dfrac{d^{2}+2k\ell
d+\ell^{2}\cdot(k^{2}+1)}{k\cdot[2dk\ell+d^{2}+\ell^{2}(k^{2}+1)]}=\dfrac{1}{k}=\tan\theta;\\\
\\\ \tan(\theta_{1}+\theta_{2})&=&\tan\theta\end{array}$ (6)

By (6) and part (b) of the trigonometric preliminaries, it follows that
$\theta_{1}+\theta_{2}=\theta$, which combined with (3), clearly establishes
that the pair $(x_{d},y_{d})$ is a solution to (1).

Now, the converse. Suppose that $(x,y)$ is a positive integer solution to (1).

Then

$\left(0<\dfrac{1}{x}\leq 1,\ \ 0<\dfrac{\ell}{y}\leq\ell,\ \
0<\dfrac{1}{k}\leq 1\right)$ (7)

Using (7), the trigonometric preliminaries, parts (a) and (b) and by taking
tangent of both sides of (1), we obtain,

$\dfrac{\frac{1}{x}+\frac{\ell}{y}}{1-\frac{1}{x}\frac{\ell}{y}}=\dfrac{1}{k}$

or equivalently

(Note that since $0<\dfrac{1}{k}\leq 1$. The equal sides of (1) can be utmost
equal to $\dfrac{\pi}{4}$ )

$\begin{array}[]{rcl}xy-\ell&=&k(y+x\ell);\\\ \\\
y\cdot(x-k)&=&\ell\cdot(1+kx)\end{array}$ (8)

Equation (8) shows that $y$ is a divisor of the product $\ell(1+kx)$. But, by
(1), we know that gcd$(\ell,y)=1$. Thus, by Result 1 (Euclid’s lemma), it
follows that $y$ must divide $1+kx$; and so,

$\left\\{\begin{array}[]{l}1+kx=y\cdot v\\\ \\\ v\ {\rm a\ positive\
integer}\end{array}\right\\}$ (9)

By (9) and (8) we have that,

$x=\ell\cdot v+k$ (10)

From (9) and (10) we further get

$1+k(\ell v+k)=yv;$

or equivalently

$k^{2}+1=(y-\ell k)\cdot v$ (11)

Since $v$ is a positive integer, equation (11) shows that $(y-\ell k)$ is a
positive integer divisor of $k^{2}+1$. Let $y-\ell k=d$, $d$ a positive
divisor of $k^{2}+1$. Then $y=\ell k+d$ and by (11) and (10) we also get

$x=k+\ell\cdot\left(\dfrac{k^{2}+1}{d}\right),$

which proves that the solution $(x,y)$ has the required form.

Finally, we see by inspection that the $N$ (number of positive divisors of
$k^{2}+1$) positive integer solutions to (1) are distinct since, obviously,
all the $Ny$-coordinates are distinct. The proof is complete.

∎

##  5  A listing of nine equalities

Let $k$ and $\ell$ be positive integers such that gcd$(\ell,k^{2}+1)=1$.
Applying Theorem 1 with $d=1$ and $d=k^{2}+1$ produces two inequalities.

  1. 1.

$\arctan\left(\dfrac{1}{k+\ell(k^{2}+1)}\right)+\arctan\left(\dfrac{\ell}{k\ell+1}\right)=\arctan\left(\dfrac{1}{k}\right)$

  2. 2.

$\arctan\left(\dfrac{1}{k+\ell}\right)+\arctan\left(\dfrac{\ell}{k\ell+k^{2}+1}\right)=\arctan\left(\dfrac{1}{k}\right)$

Next, applying Theorem 1 with $k=\ell=1$, produces the equality:

  1. 3.

$\arctan\left(\dfrac{1}{3}\right)+\arctan\left(\dfrac{1}{2}\right)=\dfrac{\pi}{4}$

For $\ell=1$ and $k=2$:

  1. 4.

$\arctan\left(\dfrac{1}{3}\right)+\arctan\left(\dfrac{1}{7}\right)=\arctan\left(\dfrac{1}{2}\right)$.

For $\ell=1$ and $k=3$

  1. 5.

$\arctan\left(\dfrac{1}{11}\right)+\arctan\left(\dfrac{1}{4}\right)=\arctan\left(\dfrac{1}{3}\right)$

  2. 6.

$\arctan\left(\dfrac{1}{8}\right)+\arctan\dfrac{1}{5}=\arctan\left(\dfrac{1}{3}\right)$

For $\ell=2$ and $k=4$:

  1. 7.

$\arctan\left(\dfrac{1}{38}\right)+\arctan\left(\dfrac{2}{9}\right)=\arctan\left(\dfrac{1}{4}\right)$

  2. 8.

$\arctan\left(\dfrac{1}{6}\right)+\arctan\left(\dfrac{2}{25}\right)=\arctan\left(\dfrac{1}{4}\right)$

For $\ell=1$ and $k=6$:

  1. 9.

$\arctan\left(\dfrac{1}{43}\right)+\arctan\left(\dfrac{1}{7}\right)=\arctan\left(\dfrac{1}{6}\right)$

## References

  * [1] Unal, Hasan, Proof without words: an arctangent equality, Mathematics and Computer Education, Fall 2011, Vol. 45, No. 3, p 197.
  * [2] Rosen, Kenneth H., Elementary Number Theory and Its Applications, 5th edition, Pearson, Addison Wesley, 2005. For Result 1 (Lemma 3.4 in the above book), see page 109. For Result 2 (Theorem 7.9 in the above book), see page 252.