# Near-Optimal Design of Safe Output Feedback Controllers from Noisy Data

Luca Furieri  Baiwei Guo⋆ Andrea Martin⋆ and Giancarlo Ferrari-Trecate Authors
are with the Institute of Mechanical Engineering, École Polytechnique Fédérale
de Lausanne, Switzerland. E-mails: {luca.furieri, baiwei.guo, andrea.martin,
<EMAIL_ADDRESS>Guo and Andrea Martin contributed
equally to this work.Research supported by the Swiss National Science
Foundation under the NCCR Automation (grant agreement 51NF40_80545).

###### Abstract

Recent work in data-driven control has revived behavioral theory to perform a
variety of complex control tasks, by directly plugging libraries of past
input-output trajectories into optimal control problems. Despite recent
advances, two key aspects remain unclear when the data are corrupted by noise:
how can safety be guaranteed, and to what extent is the control performance
affected? In this work, we provide a quantitative answer to these questions.
In particular, we formulate a robustly safe version of the recently introduced
Behavioral Input-Output Parametrization (BIOP) for the optimal predictive
control of unknown constrained systems. The proposed framework has three main
advantages: 1) it allows one to safely operate the system while explicitly
quantifying, as a function of the noise level corrupting the data, how much
the performance degrades, 2) it can be used in combination with state-of-the-
art impulse response estimators, and finally, being a data-driven approach, 3)
the state-space parameters and the initial state need not be specified for
controller synthesis. We corroborate our results through numerical
experiments.

Constrained control, Data-driven control, Learning-Based Control, Linear
Systems, Optimal Control.

##  1 Introduction

Many safety-critical engineering systems that play a crucial role in our
modern society are becoming too complex to be accurately modeled through
white-box state-space models [1, 2]. As a consequence, most contemporary
control philosophies envision unknown black-box systems for which a safe and
optimal behavior must be attained by solely relying on a collection of
historical system’s output trajectories in response to different inputs.

Widely speaking, controllers for unknown systems can be designed according to
two paradigms. _Model-based_ methods follow a two-step procedure: first, data
are exploited to identify the system parameters, and then a suitable
controller is computed for the estimated model. On the other hand, _model-
free_ methods learn the best control policy iteratively, without explicitly
reconstructing an internal representation of the dynamical system. For a
description of the advantages and limitations of both approaches, we refer to
[3], among recent surveys.

Given the intricacy of deriving rigorous suboptimality and sample-complexity
bounds, most recent model-based and model-free approaches have focused on
basic Linear Quadratic Regulator (LQR) and Linear Quadratic Gaussian (LQG)
control problems as suitable benchmarks to establish how machine learning can
be interfaced to the continuous action spaces typical of control [4, 5, 6, 7,
8, 9, 10, 11]. When it comes to complex control tasks, it is more challenging
to perform a thorough probabilistic analysis. Recent advances include [12, 13]
for constrained and distributed LQR control with direct state measurements,
respectively, and [14] for distributed output-feedback LQG.

Model-based methods may pose a difficulty when it comes to accurately
identifying the state-space model of a system whose state dimension is very
large; this is the case, for instance, for complex networked systems such as
the power grid, brain and traffic networks [2]. A promising data-driven
approach that aims at bypassing a parametric state-space description of the
system dynamics, while still being conceptually simple to implement for the
users, hinges on the _behavioral framework_ [15]. This approach has gained
renewed interest with the introduction of Data-EnablEd Predictive Control
(DeePC) [16, 17, 18, 19], which established that constrained output reference
tracking can be effectively tackled in a Model-Predictive-Control (MPC)
fashion by plugging adequately generated historical data into a convex
optimization problem. In parallel, [20] introduced data-driven formulations
for some controller design tasks. These works inspired several extensions
including closed-loop control with stability guarantees [21], maximum-
likelihood (ML) identification for control [22, 23], and nonlinear variants
[24].

In many scenarios, however, exact data are not available. For instance, data
can be corrupted by measurement noise or even by malicious attacks intended at
fatally compromising the safety [25], the quality, and the reliability of the
synthesized control policies. It is therefore essential that data-driven
controllers are endowed with quantifiable robustness properties. While some
approaches have been suggested, e.g. [17, 26, 23, 27, 28], it remains largely
unexplored how much noise-corrupted data affect the performance and the safety
of data-driven control systems. Recently, [29, 30] have derived suboptimality
[30] and sample-complexity [29] bounds for LQR through direct behavioral
formulations based on 1) Linear Matrix Inequalities (LMI) and 2) the System
Level Synthesis (SLS) approach, respectively. A limitation is that the
internal system states must be measured directly, which is infeasible for
several large-scale systems [2]. Furthermore, while [30] proves that for low-
enough noise a high-performing and robustly stabilizing controller can be
found, the corresponding suboptimality is not quantified explicitly. To
address these open points, in [31] we have formulated a Behavioral Input-
Output Parametrization (BIOP) of linear control policies which makes it
possible to derive noise-dependent suboptimality analysis for output-feedback
LQG, solely based on a non-parametric estimation of the system open-loop
responses. In order to deploy these systems in real-world scenarios, it is
however important to include safety guarantees in the analysis. All the
mentioned works [29, 30, 31, 7] do not include such safety requirements.

A preliminary version of this work has recently been submitted to the 60th
IEEE Conference on Decision and Control [31]. As opposed to [31], this paper
includes safety constraints in the analysis, thus addressing a novel and
independent set of challenges and results. We outline the new contributions
with more detail in the following subsection. Furthermore, this work includes
all the technical proofs. Last, new numerical experiments are developed to
consider safety constraints and to explicitly include the estimation procedure
of [23].

###  1.1 Contributions

Similar to MPC, a key aspect of the DeePC framework is that safety
requirements can be directly included in the problem formulation. However, it
is difficult to analyze how much the performance degrades and to what extent
safety is compromised in the presence of noise. Motivated by this challenge,
we propose a method for designing safe and near-optimal output-feedback
control policies solely based on noise-corrupted data, while explicitly
characterizing the incurred suboptimality and robustness as a function of the
noise level.

First, we show that worst-case safe and optimal output-feedback controllers
can be designed by extending the Behavioral Input-Output Parametrization
(BIOP) introduced in [31]. We showcase that our approach extends beyond
standard least-squares (LS) system identification [32] and may encompass, for
instance, data-enabled Kalman filtering [26], and the recently proposed Signal
Matrix Model (SMM) [23, 22]. Second, we develop a new approximate optimization
problem that guarantees safety while accounting for noise-corrupted data. We
provide original insight on how this program can be solved in practice by
introducing useful hyper-parameters. Third, we show that the incurred level of
suboptimality converges to zero approximately as a linear function of the size
of noise corrupting the data. The corresponding analysis differs from that of
[31], in that a feasible solution to the proposed optimization problem must be
characterized analytically while taking the safety constraints into account.
On a more general level, our analysis has been inspired by [12], which first
combined safe robust control tools with classical system identification for
state-feedback systems. To the best of our knowledge, no other existing
framework enables designing safe control policies with explicit noise-
dependent suboptimality guarantees in the output-feedback case, neither in a
behavioral context or when combining classical system-identification and
control.

###  1.2 Notation

We use $\mathbb{R}$ and $\mathbb{N}$ to denote the sets of real numbers and
non-negative integers, respectively. We use $I_{n}$ to denote the identity
matrix of size $n\times n$, $0_{m\times n}$ to denote the zero matrix of size
$m\times n$ and $1_{n}$ to denote the vector of all ones of length $n$. We
write $M=\text{blkdiag}(M_{1},\dots,M_{N})$ to denote a block-diagonal matrix
with $M_{1},\dots,M_{N}\in\mathbb{R}^{m\times n}$ on its diagonal block
entries, and for
$\mathbf{M}=\begin{bmatrix}M_{1}^{\mathsf{T}}&\dots&M_{N}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$
we define the block-Toeplitz matrix

$\footnotesize\text{Toep}_{m\times
n}\left(\mathbf{M}\right)=\begin{bmatrix}M_{1}&0_{m\times n}&\dots&0_{m\times
n}\\\ M_{2}&M_{1}&\dots&0_{m\times n}\\\ \vdots&\vdots&\ddots&\vdots\\\
M_{N}&M_{N-1}&\dots&M_{1}\end{bmatrix}.$

More concisely, we will write $\text{Toep}(\cdot)$ when the dimensions of the
blocks are clear from the context. The Kronecker product between
$M\in\mathbb{R}^{m\times n}$ and $P\in\mathbb{R}^{p\times q}$ is denoted as
$M\otimes P\in\mathbb{R}^{mp\times nq}$. For a vector $v\in\mathbb{R}^{n}$ and
a matrix $A\in\mathbb{R}^{m\times n}$ we denote as $\left\lVert
v\right\rVert_{p}$, $\left\lVert A\right\rVert_{p}$, their standard $p$-norm
and induced $p$-norms, respectively [33]. For a symmetric matrix $M$, we write
$M\succ 0$ or $M\succeq 0$ if and only if it is positive definite or positive
semidefinite, respectively. We say that $x\sim\mathcal{D}(\mu,\Sigma)$ if the
random variable $x\in\mathbb{R}^{n}$ follows a distribution with mean
$\mu\in\mathbb{R}^{n}$ and covariance matrix $\Sigma\in\mathbb{R}^{n\times
n},\Sigma\succeq 0$.

A finite-horizon trajectory of length $T$ is a sequence
$\omega(0),\omega(1),\dots,\omega(T-1)$ with $\omega(t)\in\mathbb{R}^{n}$ for
every $t=0,1,\dots,T-1$, which can be compactly written as

$\bm{\omega}_{[0,T-1]}=\begin{bmatrix}\omega^{\mathsf{T}}(0)&\omega^{\mathsf{T}}(1)&\dots&\omega^{\mathsf{T}}(T-1)\end{bmatrix}^{\mathsf{T}}\in\mathbb{R}^{nT}\,.$

When the value of $T$ is clear from the context, we will omit the subscript
$[0,T-1]$. For a finite-horizon trajectory $\bm{\omega}_{[0,T-1]}$ we also
define the Hankel matrix of depth $L$ as

$\footnotesize\mathcal{H}_{L}(\bm{\omega}_{[0,T-1]})=\begin{bmatrix}\omega(0)&\omega(1)&\dots&\omega(T-L)\\\
\omega(1)&\omega(2)&\dots&\omega(T-L+1)\\\ \vdots&\vdots&\ddots&\vdots\\\
\omega(L-1)&\omega(L)&\dots&\omega(T-1)\end{bmatrix}\,.$

###  1.3 Paper structure

The safe optimal control problem is introduced in Section 2. In Section 3 we
formulate the problem in a data-driven way, first assuming perfect data and
then moving to noise-corrupted data. Moreover, we present our controller
synthesis methodology and discuss its numerical implementation. Section 4
quantifies the suboptimality incurred by the new controller as a function of
the level of noise corrupting the data. We present numerical experiments in
Section 5 and conclude the paper in Section 6.

##  2 Problem Statement

We consider a linear system with output observations, whose state-space
representation is given by

$\displaystyle x(t+1)=Ax(t)+Bu(t),~{}~{}y(t)=Cx(t)+v(t)\,,$ (1)

where $x(t)\in\mathbb{R}^{n}$ is the state of the system and $x(0)=x_{0}$ for
a predefined $x_{0}\in\mathbb{R}^{n}$, $u(t)\in\mathbb{R}^{m}$ is the control
input, $y(t)\in\mathbb{R}^{p}$ is the observed output, and
$v(t)\in\mathbb{R}^{p}$ denotes measurement noise
$v(t)\sim\mathcal{D}(0,\Sigma_{v})$, with $\Sigma_{v}\succ 0$. The system is
controlled through a time-varying, dynamic affine control policy

$u(t)=\sum_{k=0}^{t}K_{t,k}y(k)+g_{t}+w(t)\,,$ (2)

where $K_{t,k}$ and $g_{t}$ are the linear and affine parts of the policy,
respectively, and $w(t)\in\mathbb{R}^{m}$ denotes noise on the input
$w(t)\sim\mathcal{D}(0,\Sigma_{w})$ with $\Sigma_{w}\succeq 0$, which acts as
a process noise. Furthermore, we assume that the noise is bounded with

$\left\lVert w\right\rVert_{\infty}\leq w_{\infty}\,,\quad\left\lVert
v\right\rVert_{\infty}\leq v_{\infty}\,,$

where $w_{\infty},v_{\infty}>0$. We consider the standard problem of
synthesizing a feedback control policy that minimizes the expected value with
respect to the disturbances of a quadratic objective defined over future
input-output trajectories of length $N\in\mathbb{N}$:

$J^{2}:=\mathbb{E}_{w,v}\left[\sum_{t=0}^{N-1}\left(y(t)^{\mathsf{T}}Q_{t}y(t)+u(t)^{\mathsf{T}}R_{t}u(t)\right)\right]\,,$
(3)

where $Q_{t}\succeq 0$ and $R_{t}\succ 0$ for every $t=0,\dots,N-1$.

The problem is made more challenging by the requirement that the outputs and
the inputs satisfy safety constraints

$\begin{bmatrix}y(t)\\\
u(t)\end{bmatrix}\in\Gamma\subseteq\mathbb{R}^{p+m}\,,\quad\forall
t=0,\dots,N-1\,,$ (4)

where $\Gamma$ is a nonempty polytope defined as

$\Gamma=\left\\{(y,u)\in(\mathbb{R}^{p},\mathbb{R}^{m})|~{}F_{y}y\leq
b_{y},F_{u}u\leq b_{u}\right\\}\,,$ (5)

with $F_{y}\in\mathbb{R}^{s\times p}$, $F_{u}\in\mathbb{R}^{s\times m}$ and
$b_{y},b_{u}\in\mathbb{R}^{s}$. Despite (3) being convex in the input and
output trajectories and $\Gamma$ being polytopic, we highlight that minimizing
(3) subject to (1), (2) and (4) is a non-convex problem in the control policy
parameters $K_{t,k},g_{t}$. We refer the interested reader to [34, 35].

###### Remark 1

In this work, we focus on solving and analyzing a finite-horizon control
problem, which represents one iteration of a receding-horizon MPC
implementation. It is therefore appropriate to compare the proposed approach
with a _single_ iteration of the DeePC [16, 21]. The main difference is that
we perform _closed-loop predictions_ , i.e., we optimize over feedback
policies $\pi(\cdot)$ such that $u(t)=\pi(y(t),\dots,y(0))$, while the DeePC
[16, 21] performs _open-loop predictions_ , i.e., it directly optimizes over
input sequences $u(0),u(1),u(N-1)$. It is well-known that closed-loop
predictions are less conservative than open-loop ones for MPC [36]. Indeed, by
setting $K_{t,k}=0$ in (2) the closed-loop policy reduces to an open-loop one.
Most notably, closed-loop policies may preserve feasibility for significantly
longer prediction horizons [36]. Naturally, the price to pay for such
performance improvement is an increased computational burden due to the larger
dimensionality of the problem.

###  2.1 Convex design through the IOP

By leveraging tools offered by the framework of the Input-Output
Parametrization111The IOP yields a convex representation of input-output
closed-loop responses. It is numerically stable for the case of infinite-
horizon stable plants and for finite-horizon control problems; we refer the
interest reader to [37] for details. (IOP) [38], we formulate a convex program
that computes the optimal safe feedback control policy by searching over the
input-output closed-loop responses. The state-space equations (1) provide the
following relations between trajectories

$\displaystyle\mathbf{x}_{[0,N-1]}=\mathbf{P}_{A}(:,0)x_{0}+\mathbf{P}_{B}\mathbf{u}_{[0,N-1]}\,,$
(6)
$\displaystyle\mathbf{y}_{[0,N-1]}=\mathbf{C}\mathbf{x}_{[0,N-1]}+\mathbf{v}_{[0,N-1]}\,,$
(7)

where $\mathbf{P}_{A}(:,0)$ denotes the first block-column of $\mathbf{P}_{A}$
and

$\displaystyle\mathbf{P}_{A}=(I-\mathbf{Z}\mathbf{A})^{-1}\,,$
$\displaystyle\quad\mathbf{P}_{B}=(I-\mathbf{Z}\mathbf{A})^{-1}\mathbf{Z}\mathbf{B}\,,$
$\displaystyle\mathbf{A}=I_{N}\otimes A\,,$
$\displaystyle\quad\mathbf{B}=I_{N}\otimes B\,,$
$\displaystyle\mathbf{C}=I_{N}\otimes C\,,$
$\displaystyle\quad\mathbf{Z}=\begin{bmatrix}0_{n\times n(N-1)}&0_{n\times
n}\\\ I_{n(N-1)}&0_{n(N-1)\times n}\end{bmatrix}\,.$

A few comments on the used notation are in order. First, note that the matrix
$\mathbf{Z}$ is the block-downshift operator. Second, from now on we
equivalently denote $\mathbf{G}=\mathbf{CP}_{B}$ to highlight that
$\mathbf{G}$ is a block-Toeplitz matrix containing the first $N$ components of
the impulse response of the plant $\mathbf{G}(z)=C(zI-A)^{-1}B$. Last, the
matrix $\mathbf{CP}_{A}(:,0)$ contains the entries of the observability matrix
$CA^{i}$ for $i=0,\dots,N-1$. We denote the model-based free response of the
system as $\mathbf{y}_{x_{0}}=\mathbf{CP}_{A}(:,0)x_{0}$. The control policy
can be rewritten as:

$\mathbf{u}_{[0,N-1]}=\mathbf{K}\mathbf{y}_{[0,N-1]}+\mathbf{g}+\mathbf{w}_{[0,N-1]}\,,$
(8)

where $\mathbf{K}$ and $\mathbf{g}$ are defined as:

${\footnotesize\mathbf{K}=\begin{bmatrix}K_{0,0}&0_{m\times
p}&\dots&0_{m\times p}\\\ K_{1,0}&K_{1,1}&\ddots&0_{m\times p}\\\
\vdots&\vdots&\ddots&\vdots\\\
K_{N-1,0}&K_{N-1,1}&\dots&K_{N-1,N-1}\end{bmatrix},\mathbf{g}=\begin{bmatrix}g_{0}\\\
g_{1}\\\ \vdots\\\ g_{N-1}\end{bmatrix}.}$ (9)

The safety constraints (4)-(5) take the form

$\max_{\left\lVert\mathbf{v}\right\rVert_{\infty}\leq
v_{\infty},~{}\left\lVert\mathbf{w}\right\rVert_{\infty}\leq
w_{\infty}}\mathbf{F}_{y}\mathbf{y}\leq\mathbf{b}_{y}\,,\max_{\left\lVert\mathbf{v}\right\rVert_{\infty}\leq
v_{\infty},\left\lVert\mathbf{w}\right\rVert_{\infty}\leq
w_{\infty}}\mathbf{F}_{u}\mathbf{u}\leq\mathbf{b}_{u}\,,$ (10)

where $\max(\cdot)$ is intended row-wise, $\mathbf{F}_{y}=I_{N}\otimes F_{y}$,
$\mathbf{F}_{u}=I_{N}\otimes F_{u}$, $\mathbf{b}_{y}=1_{N}\otimes b_{y}$,
$\mathbf{b}_{u}=1_{N}\otimes b_{u}$. By plugging the controller (8) into
(6)-(7), it is easy to derive the relationships

$\displaystyle\begin{bmatrix}\mathbf{y}\\\
\mathbf{u}\end{bmatrix}=\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}\end{bmatrix}\begin{bmatrix}\mathbf{v}+\mathbf{y}_{x_{0}}\\\
\mathbf{w}\end{bmatrix}+\begin{bmatrix}\mathbf{Gq}\\\
\mathbf{q}\end{bmatrix}\,,$ (11)

where

$\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}\end{bmatrix}=\begin{bmatrix}(I-\mathbf{GK})^{-1}&(I-\mathbf{GK})^{-1}\mathbf{G}\\\
\mathbf{K}(I-\mathbf{GK})^{-1}&(I-\mathbf{KG})^{-1}\end{bmatrix}\,,$ (12)

and

$\mathbf{q}=(I-\mathbf{KG})^{-1}\mathbf{g}=\bm{\Phi}_{uu}\mathbf{g}\,.$ (13)

The parameters
$\bm{\Phi}=(\bm{\Phi}_{yy},\bm{\Phi}_{yu},\bm{\Phi}_{uy},\bm{\Phi}_{uu})$
represent the four closed-loop responses defining the relationship between
disturbances and input-output signals, while $\mathbf{q}$ represents the
affine part of the disturbance-feedback control policy [34, 35]. To achieve a
convex reformulation of the control problem under consideration, we extend the
IOP from [38] to account for the safety constraints (10) in a convex way. The
result is summarized in the next Proposition, whose proof is reported in
Appendix .1 for completeness.

###### Proposition 1

Consider the LTI system (1) evolving under the control policy (8) within a
horizon of length $N\in\mathbb{N}$. Then:

  1. $i)$

For any control policy $(\mathbf{K},\mathbf{g})$ that complies with the safety
constraints, there exist four matrices
($\bm{\Phi}_{yy},\bm{\Phi}_{yu},\bm{\Phi}_{uy},\bm{\Phi}_{uu}$) and a vector
$\mathbf{q}$ such that $\mathbf{K}=\bm{\Phi}_{uy}\bm{\Phi}_{yy}^{-1}$,
$\mathbf{g}=\bm{\Phi}_{uu}^{-1}\mathbf{q}$. Furthermore, for all
$j=1,\dots,sN$,

$\displaystyle\begin{bmatrix}I&-\mathbf{G}\end{bmatrix}\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}\end{bmatrix}=\begin{bmatrix}I&0\end{bmatrix},$
(14) $\displaystyle\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}\end{bmatrix}\begin{bmatrix}-\mathbf{G}\\\
I\end{bmatrix}=\begin{bmatrix}0\\\ I\end{bmatrix},$ (15)
$\displaystyle\left\lVert\begin{bmatrix}v_{\infty}(F_{y,j}\bm{\Phi}_{yy})^{\mathsf{T}}\\\
w_{\infty}(F_{y,j}\bm{\Phi}_{yu})^{\mathsf{T}}\end{bmatrix}\right\rVert_{1}+F_{y,j}(\mathbf{Gq}+\bm{\Phi}_{yy}\mathbf{y}_{x0})\leq\mathbf{b}_{y,j}\,,$
(16)
$\displaystyle\left\lVert\begin{bmatrix}v_{\infty}(F_{u,j}\bm{\Phi}_{uy})^{\mathsf{T}}\\\
w_{\infty}(F_{u,j}\bm{\Phi}_{uu})^{\mathsf{T}}\end{bmatrix}\right\rVert_{1}+F_{u,j}(\mathbf{q}+\bm{\Phi}_{uy}\mathbf{y}_{x_{0}})\leq\mathbf{b}_{u,j}\,,$
(17)
$\displaystyle\bm{\Phi}_{yy},\bm{\Phi}_{yu},\bm{\Phi}_{uy},\bm{\Phi}_{uu}\text{
with causal sparsities }\,,$ (18)

where $F_{y,j}\in\mathbb{R}^{1\times Np}$, $F_{u,j}\in\mathbb{R}^{1\times Nm}$
and $\mathbf{b}_{u,j},\mathbf{b}_{y,j}\in\mathbb{R}$ are the $j$-th row of
$\mathbf{F}_{y}$, $\mathbf{F}_{u}$ and $\mathbf{b}_{u},\mathbf{b}_{y}$,
respectively.

  2. $ii)$

For any four matrices
($\bm{\Phi}_{yy},\bm{\Phi}_{yu},\bm{\Phi}_{uy},\bm{\Phi}_{uu}$) lying in the
affine subspace (14)-(18) and any vector $\mathbf{q}\in\mathbb{R}^{mN}$, the
matrix $\mathbf{K}=\bm{\Phi}_{uy}\bm{\Phi}_{yy}^{-1}$ is causal as per (9) and
it yields the closed-loop responses
($\bm{\Phi}_{yy},\bm{\Phi}_{yu},\bm{\Phi}_{uy},\bm{\Phi}_{uu}$). Moreover, the
affine policy $(\mathbf{K},\mathbf{g}$) with
$\mathbf{g}=\bm{\Phi}_{uu}^{-1}\mathbf{q}$ complies with the safety
constraints.

We remark that the IOP is well-suited to a data-driven output-feedback setup,
as all affine control policies are directly parametrized through the impulse
response parameters $\mathbf{G}$, without requiring an internal state-space
representation. This is useful for two reasons. First, when dealing with
unknown systems, the state-space parameters $(A,B,C,x_{0})$ can only be
estimated up to an unknown change of variable, which may be problematic for
defining the cost and the noise statistics [39]. Second, several large-scale
systems feature a very large number of states, but a comparably small number
of inputs and outputs, that is $n>>\max(m,p)$. In such applications, it is
advantageous to bypass a state-space representation and directly deal with
$\mathbf{G}$, whose dimensions do not depend on $n$.

From now on, in the interest of simplifying the expressions appearing
throughout the next sections and without any loss of generality333One can
redefine
$\tilde{\mathbf{y}}=\begin{bmatrix}1&\mathbf{y}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$,
$\overline{\mathbf{v}}=\begin{bmatrix}1&\mathbf{v}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$,
$\overline{\mathbf{C}}=\begin{bmatrix}0_{1\times Nn}\\\
\mathbf{C}\end{bmatrix}$,
$\overline{\mathbf{K}}=\begin{bmatrix}\mathbf{g}&\mathbf{K}\end{bmatrix}$ and
$\overline{\bm{\Phi}}$ as per (12) with $\overline{\mathbf{G}}$ and
$\overline{\mathbf{K}}$ in place of $\mathbf{G}$ and $\mathbf{K}$,
respectively. Minor modifications to (16)-(17) are needed as well., we will
assume that $\mathbf{q}=\mathbf{g}=0_{Nm\times 1}$. Namely, we will focus on
_linear_ control policies. We are now ready to establish a convex formulation
of the optimal control problem under study.

###### Proposition 2

Consider the LTI system (1). The linear control policy that achieves the
minimum of the cost functional (3) is given by
$\mathbf{K}=\bm{\Phi}_{uy}\bm{\Phi}_{yy}^{-1}$, where
$\bm{\Phi}_{uy},\bm{\Phi}_{yy}$ are optimal solutions to the following convex
program:

$\displaystyle~{}\min_{\bm{\Phi}}\left\lVert\begin{bmatrix}\mathbf{Q}^{\frac{1}{2}}&0\\\
0&\mathbf{R}^{\frac{1}{2}}\end{bmatrix}\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}\end{bmatrix}\begin{bmatrix}\bm{\Sigma}^{\frac{1}{2}}_{v}&0&\mathbf{y}_{x_{0}}\\\
0&\bm{\Sigma}^{\frac{1}{2}}_{w}&0\end{bmatrix}\right\rVert_{F}^{2}$ (19)
$\displaystyle\operatorname*{subject~{}to}~{}\eqref{eq:ach1}-\eqref{eq:ach3}\,,$

where $\mathbf{Q}=\text{blkdiag}(Q_{0},\dots,Q_{N-1})$,
$\mathbf{R}=\text{blkdiag}(R_{0},\dots,R_{N-1})$,
$\bm{\Sigma}_{v}=I_{N}\otimes\Sigma_{v}$,
$\bm{\Sigma}_{w}=I_{N}\otimes\Sigma_{w}$ and where (16)-(17) are evaluated at
$\mathbf{q}=0$.

######  Proof 2.1.

We refer to Proposition 2 of [31] for a complete derivation of the cost
function, which is based on standard norm arguments. To conclude the proof, it
suffices to notice that the objective function and the safety constraints
(16)-(17) are convex in $\bm{\Phi}$.

When the system parameters $(A,B,C,x_{0})$ are known, a globally optimal
solution
$(\bm{\Phi}^{\star}_{yy},\bm{\Phi}^{\star}_{yu},\bm{\Phi}^{\star}_{uy},\bm{\Phi}^{\star}_{uu}$)
for problem (19) can be efficiently computed with off-the-shelf solvers. The
corresponding globally optimal and safe control policy is then recovered as
$\mathbf{K}^{\star}=\bm{\Phi}_{uy}^{\star}(\bm{\Phi}^{\star}_{yy})^{-1}$.

The rest of the paper addresses the following two questions:

  * Q1)

How can we compute a safe control policy with performance close to that of
$\mathbf{K}^{\star}$, solely based on libraries of _noisy_ input-output
trajectories?

  * Q2)

How can we _quantify_ the performance gap with respect to $\mathbf{K}^{\star}$
as the noise corrupting the data increases?

##  3 Robustly Safe Controller Synthesis From Noisy Data through the BIOP

In this section, we answer question Q1) by developing a method to synthesize
near-optimal safe controllers directly from noisy data. We exploit behavioral
theory [15] to compute approximate system responses and we cast a safe version
of the BIOP that we have introduced for unconstrained systems in [31].

First, we briefly recall _Willems’ Fundamental Lemma_ and useful results from
behavioral linear system theory. Second, we show how to compute the optimal
safe controller $\mathbf{K}^{\star}$ in a data-driven way when the data are
not corrupted by noise. Last, we establish an optimization problem based on
the BIOP that tightly approximates the optimal and safe control policy in the
presence of noisy data. We conclude the section by offering novel insights on
its properties and its numerical implementation. Hence, we significantly
generalize the theoretical validity and the applicability of previous results
on combining system identification and robust control for state-feedback
systems [12].

###  3.1 Willems’ lemma and behavioral theory for synthesizing safe
controllers

We recall the definition of persistency of excitation and the result known as
the _Fundamental Lemma_ for LTI systems [40].

######  Definition 3.2.

We say that $\mathbf{u}^{h}_{[0,T-1]}$ is _persistently exciting_ (PE) of
order $L$ if the Hankel matrix $\mathcal{H}_{L}(\mathbf{u}^{h}_{[0,T-1]})$ is
full row-rank.

A necessary condition for the matrix
$\mathcal{H}_{L}(\mathbf{u}^{h}_{[0,T-1]})$ to be full row-rank is that it has
at least as many columns as rows. It follows that the input trajectory
$\mathbf{u}^{h}_{[0,T-1]}$ must be long enough to satisfy $T\geq(m+1)L-1$.

######  Lemma 3.3 (Theorem 3.7, [40]).

Consider system (1). Assume that $(A,B)$ is controllable and that there is no
noise. Let $\\{\mathbf{y}^{h}_{[0,T-1]},\mathbf{u}^{h}_{[0,T-1]}\\}$ be a
system trajectory of length $T$ that has been recorded during a past
experiment. Then, if $\mathbf{u}^{h}_{[0,T-1]}$ is PE of order $n+L$, the
signals $\mathbf{y}^{\star}_{[0,L-1]}\in\mathbb{R}^{pL}$ and
$\mathbf{u}^{\star}_{[0,L-1]}\in\mathbb{R}^{mL}$ are trajectories of (1) if
and only if there exists $g\in\mathbb{R}^{T-L+1}$ such that

$\begin{bmatrix}\mathcal{H}_{L}(\mathbf{y}^{h}_{[0,T-1]})\\\
\mathcal{H}_{L}(\mathbf{u}^{h}_{[0,T-1]})\end{bmatrix}g=\begin{bmatrix}\mathbf{y}^{\star}_{[0,L-1]}\\\
\mathbf{u}^{\star}_{[0,L-1]}\end{bmatrix}\,.$ (20)

We proceed by showing how Lemma 3.3 allows one to derive a data-driven
formulation of (19) when the data are not noisy. We work under the following
standard assumptions.

###### Assumption 1

The data-generating LTI system (1) is such that $(A,B)$ is controllable and
$(A,C)$ is observable.

###### Assumption 2

The following data are available:

  1. $i)$

a _recent_ system trajectory of length $T_{ini}$:
$\left\\{\mathbf{y}^{r}_{[0,T_{ini}-1]},\mathbf{u}^{r}_{[0,T_{ini}-1]}\right\\}$,
with $\mathbf{y}^{r}_{[0,T_{ini}-1]}=\mathbf{y}_{[-T_{ini},-1]}$ and
$\mathbf{u}^{r}_{[0,T_{ini}-1]}=\mathbf{u}_{[-T_{ini},-1]}$,

  2. $ii)$

a _historical_ system trajectory of length $T$:
$\left\\{\mathbf{y}^{h}_{[0,T-1]},\mathbf{u}^{h}_{[0,T-1]}\right\\}$, with
$\mathbf{y}^{h}_{[0,T-1]}=\mathbf{y}_{[-T_{h},-T_{h}+T-1]}$ and
$\mathbf{u}^{h}_{[0,T-1]}=\mathbf{u}_{[-T_{h},-T_{h}+T-1]}$ for
$T_{h}\in\mathbb{N}$ such that $T_{h}>T+T_{ini}$.

###### Assumption 3

The historical input trajectory $\mathbf{u}^{h}_{[0,T-1]}$ is PE of order
$n+T_{ini}+N$, where $T_{ini}\geq l$ and $l$ is the smallest integer such that
the matrix

$\begin{bmatrix}C^{\mathsf{T}}&(CA)^{\mathsf{T}}&\dots&(CA^{l-1})^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}\,,$

has full row-rank. Note that if Assumption 1 holds, then $l\leq n$.

The _historical_ data are to be used in substitution of the system model,
while the _recent_ data reflect the system initial state
$x_{0}\in\mathbb{R}^{n}$ [41]. Next, we observe that (20) can be exploited to
compute the safe and optimal controller $\mathbf{K}^{\star}$ in a convex way
through a behavioral version of the IOP.

######  Theorem 3.4 (Safe Behavioral IOP).

Consider the LTI system (1), whose parameters $(A,B,C,x_{0})$ are _unknown_ ,
and let Assumptions 1-3 hold. Further assume that the historical and recent
trajectories are _not_ affected by noise. Let $(G,g)$ be any solutions to the
linear system of equations

$\begin{bmatrix}U_{p}\\\ Y_{p}\\\
U_{f}\end{bmatrix}\begin{bmatrix}G&g\end{bmatrix}=\begin{bmatrix}0_{mT_{ini}\times
m}&\mathbf{u}^{r}_{[0,T_{ini}-1]}\\\ 0_{pT_{ini}\times
m}&\mathbf{y}^{r}_{[0,T_{ini}-1]}\\\ \begin{bmatrix}I_{m}&0_{m\times
m(N-1)}\end{bmatrix}^{\mathsf{T}}&0_{mN\times 1}\end{bmatrix},$ (21)

where $\begin{bmatrix}U_{p}\\\
U_{f}\end{bmatrix}=\mathcal{H}_{T_{ini}+N}(\mathbf{u}^{h}_{[0,T-1]})$ and
$\begin{bmatrix}Y_{p}\\\
Y_{f}\end{bmatrix}=\mathcal{H}_{T_{ini}+N}(\mathbf{y}^{h}_{[0,T-1]})$. Then,
the optimization problem (19) is equivalent to

$\displaystyle\min_{\bm{\Phi}}\left\lVert\begin{bmatrix}\mathbf{Q}^{\frac{1}{2}}&0\\\
0&\mathbf{R}^{\frac{1}{2}}\end{bmatrix}\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}\end{bmatrix}\begin{bmatrix}\bm{\Sigma}^{\frac{1}{2}}_{v}&0&Y_{f}g\\\
0&\bm{\Sigma}^{\frac{1}{2}}_{w}&0\end{bmatrix}\right\rVert_{F}^{2}$ (22)
$\displaystyle\operatorname*{subject~{}to}\begin{bmatrix}I&-\text{Toep}(Y_{f}G)\end{bmatrix}\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}\end{bmatrix}=\begin{bmatrix}I&0\end{bmatrix}\,,$
$\displaystyle\qquad\quad\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}\end{bmatrix}\begin{bmatrix}-\text{Toep}(Y_{f}G)\\\
I\end{bmatrix}=\begin{bmatrix}0\\\ I\end{bmatrix}\,,$
$\displaystyle\qquad\quad\left\lVert\begin{bmatrix}v_{\infty}\left(F_{y,j}\bm{\Phi}_{yy}\right)^{\mathsf{T}}\\\
w_{\infty}\left(F_{y,j}\bm{\Phi}_{yu}\right)^{\mathsf{T}}\end{bmatrix}\right\rVert_{1}+\left(F_{y,j}\bm{\Phi}_{yy}\right)Y_{f}g\leq\mathbf{b}_{y,j}\,,$
$\displaystyle\qquad\quad\left\lVert\begin{bmatrix}v_{\infty}\left(F_{u,j}\bm{\Phi}_{uy}\right)^{\mathsf{T}}\\\
w_{\infty}\left(F_{u,j}\bm{\Phi}_{uu}\right)^{\mathsf{T}}\end{bmatrix}\right\rVert_{1}+\left(F_{u,j}\bm{\Phi}_{uy}\right)Y_{f}g\leq\mathbf{b}_{u,j}\,,$
$\displaystyle\qquad\quad\forall j=1,\dots,sN\,,$
$\displaystyle\qquad\quad\bm{\Phi}_{yy},\bm{\Phi}_{yu},\bm{\Phi}_{uy},\bm{\Phi}_{uu}\text{
with causal sparsities}.$

######  Proof 3.5.

In problem (19), the system parameters $(A,B,C,x_{0})$ appear through the
terms $\mathbf{G}=\mathbf{CP}_{B}$ in the constraints and $\mathbf{y}_{x_{0}}$
in the cost. It is therefore sufficient to show that we are able to substitute
both elements with data as per the theorem statement. Let $G$ be any solution
(21). By rearranging the terms, each column of $G$ can be seen as a solution
to (20) associated with a zero initial condition and input
$\mathbf{u}_{[0,N-1]}=\begin{bmatrix}e_{i}^{\mathsf{T}}&0_{1\times
m(N-1)}\end{bmatrix}^{\mathsf{T}}$, where $e_{i}\in\mathbb{R}^{m}$ is the
$i$-th element of the standard orthogonal basis of $\mathbb{R}^{m}$. Since the
hypotheses of Lemma 3.3 are satisfied for $L=T_{ini}+N$, similar to
Proposition 11 of [41] we deduce that $Y_{f}G$ is the first block-column of
the system impulse response matrix, no matter which solution to (21) is
chosen. Therefore, we can equivalently substitute
$\mathbf{G}=\text{Toep}_{p\times m}(Y_{f}G)$ in the constraints (14)-(15) of
problem (19). Finally, note that $Y_{f}g$ corresponds to the trajectory
starting at $x_{0}$ (as implicitly defined by the recent trajectory
$\mathbf{y}_{[-T_{ini},-1]}$ and $\mathbf{u}_{[-T_{ini},-1]}$) when applying a
zero input [41]. Therefore, it corresponds to the true free response starting
from $x_{0}$.

Since the historical and recent data are not noisy, $Y_{f}G$ and $Y_{f}g$
yield the _true_ impulse response matrix $\mathbf{G}$ and free response
$\mathbf{y}_{0}$ and the optimal solution of (22) recovers the optimal safe
controller $\mathbf{K}^{\star}$ for the real system.

###  3.2 A robust BIOP with noise-corrupted data

In practice, exact historical and recent data are not available. As per the
noise model in the dynamics (1)-(2), we assume that historical and recent
trajectories are affected by additive noise
$w^{h}(t),w^{r}(t),v^{h}(t),v^{r}(t)$444where “$w$” and “$v$” denote input and
output noise, respectively, and the apices $r$ and $h$ denote recent data and
historical data, respectively. at all time instants, with zero expected values
and variances
$\bm{\Sigma}^{h}_{w},\bm{\Sigma}^{r}_{w},\bm{\Sigma}^{h}_{v},\bm{\Sigma}^{r}_{v}$
respectively. Hence, the matrix on the left-hand-side of (21) becomes full
row-rank almost surely and (21) will yield imprecise impulse and free
responses. This issue is well-known in the behavioral theory literature, and
several mitigation strategies have recently been proposed [16, 17, 20, 23, 26,
22]. Specifically, these methods have the goal of yielding accurate and
coherent estimators based on behavioral theory. For instance, the standard LS
estimator is given by

$\displaystyle G_{LS}$ $\displaystyle=\begin{bmatrix}\hat{U}_{p}\\\
\hat{Y}_{p}\\\ \hat{U}_{f}\end{bmatrix}^{+}\begin{bmatrix}0_{mT_{ini}\times
m}\\\ 0_{pT_{ini}\times m}\\\ \begin{bmatrix}I_{m}&0_{m\times
m(N-1)}\end{bmatrix}^{\mathsf{T}}\end{bmatrix}\,,$ (23) $\displaystyle g_{LS}$
$\displaystyle=\begin{bmatrix}\hat{U}_{p}\\\ \hat{Y}_{p}\\\
\hat{U}_{f}\end{bmatrix}^{+}\begin{bmatrix}\mathbf{u}^{r}_{[0,T_{ini}-1]}\\\
\mathbf{y}^{r}_{[0,T_{ini}-1]}\\\ 0_{mN\times 1}\end{bmatrix}\,,$ (24)

while the ML estimator [23] is computed as

$\displaystyle G_{ML}=$ $\displaystyle\operatorname*{arg\,min}_{G}$
$\displaystyle-log\left[p\left(\begin{bmatrix}\Xi_{y}\\\
Y_{f}G\end{bmatrix}|~{}G,Y_{f}\right)\right]$ (25)
$\displaystyle\operatorname*{subject~{}to}$
$\displaystyle\begin{bmatrix}\hat{U}_{p}\\\
\hat{U}_{f}\end{bmatrix}G=\begin{bmatrix}0_{mT_{ini}\times m}\\\
\begin{bmatrix}I_{m}&0_{m\times
m(N-1)}\end{bmatrix}^{\mathsf{T}}\end{bmatrix}\,,$ $\displaystyle g_{ML}=$
$\displaystyle\operatorname*{arg\,min}_{g}$ $\displaystyle-
log\left[p\left(\begin{bmatrix}\xi_{y}\\\
Y_{f}g\end{bmatrix}|~{}g,Y_{f}\right)\right]$ (26)
$\displaystyle\operatorname*{subject~{}to}$
$\displaystyle\begin{bmatrix}\hat{U}_{p}\\\
\hat{U}_{f}\end{bmatrix}g=\begin{bmatrix}\mathbf{u}^{r}_{[0,T_{ini}-1]}\\\
0_{mN\times 1}\end{bmatrix}\,,$

where $p(a|b)$ indicates the probability of event $a$ conditioned to $b$. We
will test the performance of (25)-(26) against (23)-(24) in Section 5. It is
effective to work under the assumption that, with high-probability, the errors
$\left\lVert\mathbf{G}-\widehat{\mathbf{G}}\right\rVert$ and
$\left\lVert\mathbf{y}_{0}-\widehat{\mathbf{y}}_{0}\right\rVert$ are small in
an appropriate sense. Specifically, [4, 7, 12] consider the estimation error
resulting from LS system identification. Here, we abstract from the particular
identification scheme and formalize the following assumption.

###### Assumption 4

Let $\bm{\Delta}=\mathbf{G}-\widehat{\mathbf{G}}$ and
$\bm{\delta}_{0}=\mathbf{y}_{0}-\widehat{\mathbf{y}}_{0}$. There exist
$\epsilon_{2,G},\epsilon_{\infty,G},\epsilon_{2,y},\epsilon_{\infty,y}>0$ such
that, for any sequence of noisy historical and recent data, with high
probability

$\displaystyle\left\lVert\bm{\Delta}\right\rVert_{2}\leq\epsilon_{2,G},~{}~{}~{}$
$\displaystyle\left\lVert\bm{\delta}_{0}\right\rVert_{2}\leq\epsilon_{2,y}\,,$
$\displaystyle\left\lVert\bm{\Delta}\right\rVert_{\infty}\leq\epsilon_{\infty,G},$
$\displaystyle\left\lVert\bm{\delta}_{0}\right\rVert_{\infty}\leq\epsilon_{\infty,y}\,.$

We denote $\epsilon_{2}=\max(\epsilon_{G,2},\epsilon_{0,2})$ and
$\epsilon_{\infty}=\max(\epsilon_{G,\infty},\epsilon_{0,\infty})$.

After condensing the effect of noise-corrupted data into error parameters
$\bm{\Delta},\bm{\delta}$, whose norms are bounded by
$\epsilon_{2},\epsilon_{\infty}$, we can formulate a _robustly safe_ data-
driven control problem as follows. To facilitate the derivations in the rest
of the paper, we select the weights $\mathbf{Q}$, $\mathbf{R}$,
$\bm{\Sigma}_{w}$, $\bm{\Sigma}_{v}$ to be identity matrices with appropriate
dimensions. All our results can be easily adapted to non-identity weights.

######  Proposition 3.6.

Assume that the collected historical and recent data are noisy. Let
$\widehat{\mathbf{G}},\widehat{\mathbf{y}}_{0}$ be estimators of
$\mathbf{G},\mathbf{y}_{0}$, respectively, such that Assumption 4 holds with
$\epsilon_{2},\epsilon_{\infty}>0$. Let
$J(\mathbf{G},\mathbf{K})=\sqrt{\mathbb{E}_{\mathbf{w},\mathbf{v}}\left[\mathbf{y}^{\mathsf{T}}\mathbf{y}+\mathbf{u}^{\mathsf{T}}\mathbf{u}\right]}$
denote the square root of the cost in (3). Consider the following model-based
worst-case optimal control problem:

$\displaystyle\min_{\mathbf{K}}~{}\max_{\tiny\begin{matrix}\left\lVert\bm{\Delta}\right\rVert_{2}\leq\epsilon_{2},~{}\left\lVert\bm{\delta}_{0}\right\rVert_{2}\leq\epsilon_{2}\\\
\left\lVert\bm{\Delta}\right\rVert_{\infty}\leq\epsilon_{\infty},~{}\left\lVert\bm{\delta}_{0}\right\rVert_{\infty}\leq\epsilon_{\infty}\end{matrix}}$
$\displaystyle J(\mathbf{G},\mathbf{K})$ (27)
$\displaystyle\operatorname*{subject~{}to}~{}$
$\displaystyle\eqref{eq:state_compact},\eqref{eq:output_compact},\eqref{eq:control_policy},\eqref{eq:K_sparsity},\eqref{eq:compact_constraints}\,,$

and let

$\displaystyle\bm{\Phi}_{yy}$
$\displaystyle=\widehat{\bm{\Phi}}_{yy}(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}\,,$
$\displaystyle\bm{\Phi}_{yu}$
$\displaystyle=\widehat{\bm{\Phi}}_{yy}(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}(\widehat{\mathbf{G}}+\mathbf{\Delta})\,,$
$\displaystyle\bm{\Phi}_{uy}$
$\displaystyle=\widehat{\bm{\Phi}}_{uy}(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}\,,$
$\displaystyle\bm{\Phi}_{uu}$
$\displaystyle=(I-\widehat{\bm{\Phi}}_{uy}\mathbf{\Delta})^{-1}\widehat{\bm{\Phi}}_{uu}\,.$

Then, problem (27) is equivalent to

$\displaystyle\qquad\min_{\widehat{\bm{\Phi}}}~{}~{}\max_{\bm{\Delta},\bm{\delta}_{0}}~{}~{}\left\lVert\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}\end{bmatrix}\Bigg{[}\begin{matrix}I&0&\mathbf{y}_{0}\\\
0&I&0\end{matrix}\Bigg{]}\right\rVert$ (28)
$\displaystyle\operatorname*{subject~{}to}\quad\begin{bmatrix}I&-\widehat{\mathbf{G}}\end{bmatrix}\begin{bmatrix}\widehat{\bm{\Phi}}_{yy}&\widehat{\bm{\Phi}}_{yu}\\\
\widehat{\bm{\Phi}}_{uy}&\widehat{\bm{\Phi}}_{uu}\end{bmatrix}=\begin{bmatrix}I&0\end{bmatrix},$
(29)
$\displaystyle\qquad\begin{bmatrix}\widehat{\bm{\Phi}}_{yy}&\widehat{\bm{\Phi}}_{yu}\\\
\widehat{\bm{\Phi}}_{uy}&\widehat{\bm{\Phi}}_{uu}\end{bmatrix}\begin{bmatrix}-\widehat{\mathbf{G}}\\\
I\end{bmatrix}=\begin{bmatrix}0\\\ I\end{bmatrix},$ (30)
$\displaystyle\qquad\left\lVert\begin{bmatrix}v_{\infty}\left(F_{y,j}\bm{\Phi}_{yy}\right)^{\mathsf{T}}\\\
w_{\infty}\left(F_{y,j}\bm{\Phi}_{yu}\right)^{\mathsf{T}}\end{bmatrix}\right\rVert_{1}+\left(F_{y,j}\bm{\Phi}_{yy}\right)\mathbf{y}_{0}\leq\mathbf{b}_{y,j}\,,$
(31)
$\displaystyle\qquad\left\lVert\begin{bmatrix}v_{\infty}\left(F_{u,j}\bm{\Phi}_{uy}\right)^{\mathsf{T}}\\\
w_{\infty}\left(F_{u,j}\bm{\Phi}_{uu}\right)^{\mathsf{T}}\end{bmatrix}\right\rVert_{1}+\left(F_{u,j}\bm{\Phi}_{uy}\right)\mathbf{y}_{0}\leq\mathbf{b}_{u,j}\,,$
(32) $\displaystyle\qquad\forall j=1,\dots,sN\,,$
$\displaystyle\qquad\widehat{\bm{\Phi}}_{yy},\widehat{\bm{\Phi}}_{yu},\widehat{\bm{\Phi}}_{uy},\widehat{\bm{\Phi}}_{uu}\text{
with causal sparsities.}$ (33)

The proof of Proposition 3.6 can be found in Appendix .2. We remark that the
closed-responses $\bm{\Phi}$ appearing in (28), (31) and (32) are associated
with the true impulse response, whereas the closed-loop responses
$\widehat{\bm{\Phi}}$ appearing in (29), (30) and (33) are associated with the
_estimated_ impulse response. This is because, while we are interested in
minimizing the cost and satisfying the safety constraints for the _real_
system, we can only parametrize the closed-loop responses for the approximate
system.

The robust optimization problem (28) is highly non-convex in the cost and in
the constraints because $\bm{\Phi}$ is a nonlinear function of the matrix
variables $\widehat{\bm{\Phi}}$, $\bm{\Delta}$ and $\bm{\delta}$. Therefore,
it is challenging to find a feasible solution, let alone one with very good
performance. Our first main result is to reformulate (28) into a form that can
be tackled with simple off-the-shelf optimization techniques.

###  3.3 Proposed reformulation for safe controller synthesis

We proceed to reformulate the intractable problem (28) into a more convenient
form. Specifically, our proposed approach is to 1) upper-bound the cost
function, and 2) tighten the safety constraints with more tractable
expressions. In Section 4 we will explicitly quantify the suboptimality
incurred by these approximations. At its core, this methodology is inspired by
that developed in [12] for the state-feedback case without measurement noise.
However, the addition of output-feedback and measurement noise leads to new
terms both in the cost and the safety constraints which are more challenging
to analyze.

The following two lemmas establish the basis for our reformulation. Lemma 3.7
provides the new expression which upper-bounds $J(\mathbf{G},\mathbf{K})$ and
Lemma 3.8 provides a tightened form of the safety constraints. Their rather
lengthy technical proof is reported in the Appendices .3 and .4, respectively.

######  Lemma 3.7.

Let $\widehat{\bm{\Phi}}$ denote the closed-loop responses obtained by
applying $\mathbf{K}$ to $\widehat{\mathbf{G}}$. Further assume that
$\left\lVert\widehat{\bm{\Phi}}_{uy}\right\rVert_{2}\leq\gamma$, where
$\gamma\in[0,\epsilon_{2}^{-1})$. Then, we have

$\displaystyle
J(\mathbf{G},\mathbf{K})\leq\frac{1}{1-\epsilon_{2}\gamma}\times$ (34)
$\displaystyle\times\left\lVert\begin{bmatrix}\sqrt{1+h(\epsilon_{2},\gamma,\widehat{\mathbf{G}})+h(\epsilon_{2},\gamma,\widehat{\mathbf{y}}_{0})}\widehat{\bm{\Phi}}_{yy}&\widehat{\bm{\Phi}}_{yu}&\widehat{\bm{\Phi}}_{yy}\widehat{\mathbf{y}}_{0}\\\
\sqrt{1+h(\epsilon_{2},\gamma,\widehat{\mathbf{y}}_{0})}\widehat{\bm{\Phi}}_{uy}&\widehat{\bm{\Phi}}_{uu}&\widehat{\bm{\Phi}}_{uy}\widehat{\mathbf{y}}_{0}\end{bmatrix}\right\rVert_{F},$

where

$h(\epsilon,\gamma,\mathbf{Y})=\epsilon^{2}(2+\gamma\|\mathbf{Y}\|_{2})^{2}+2\epsilon\left\lVert\mathbf{Y}\right\rVert_{2}(2+\gamma\left\lVert\mathbf{Y}\right\rVert_{2})\,.$

Lemma 3.7 exploits the upperbound
$\left\lVert\widehat{\bm{\Phi}}_{uy}\right\rVert_{2}\leq\gamma$ to establish
an explicit relationship between the cost $J(\mathbf{G},\mathbf{K})$, obtained
by applying a controller $\mathbf{K}$ to the real system $\mathbf{G}$, and the
cost $J(\widehat{\mathbf{G}},\mathbf{K})$, obtained by applying the same
controller to the estimated system $\widehat{\mathbf{G}}$. To see this, notice
that (34) is equivalently rewritten as

$\displaystyle J(\mathbf{G},\mathbf{K})$
$\displaystyle\leq\frac{\Big{(}J(\widehat{\mathbf{G}},\mathbf{K})+\|\widehat{\bm{\Phi}}_{yy}\|_{F}^{2}(h(\epsilon_{2},\gamma,\widehat{\mathbf{G}})+}{1-\epsilon_{2}\gamma}$
$\displaystyle\frac{+h(\epsilon_{2},\gamma,\widehat{\mathbf{y}}_{0}))+\|\widehat{\bm{\Phi}}_{uy}\|_{F}^{2}h(\epsilon_{2},\gamma,\widehat{\mathbf{y}}_{0})\Big{)}^{\frac{1}{2}}}{1-\epsilon_{2}\gamma}\,.$
(35)

The expression (35) upper-bounds the gap between $J(\mathbf{G},\mathbf{K})$
and $J(\widehat{\mathbf{G}},\mathbf{K})$ as a quantity that increases with
$\epsilon_{2}$ and with the norm of
$\widehat{\mathbf{G}},\widehat{\mathbf{y}}_{0},\widehat{\bm{\Phi}}$. We now
derive a tightened - yet more tractable - expression for the safety
constraints (31)-(32).

######  Lemma 3.8.

Assume $\left\lVert\widehat{\bm{\Phi}}_{uy}\right\rVert_{\infty}\leq\tau$,
where $\tau\in[0,\epsilon_{\infty}^{-1})$. Then, if $\widehat{\bm{\Phi}}$
satisfies the tightened safety constraints

$\displaystyle
f_{1,j}(\widehat{\bm{\Phi}})+f_{2,j}(\widehat{\bm{\Phi}})+f_{3,j}(\widehat{\bm{\Phi}})\leq\mathbf{b}_{y,j}\,,$
(36) $\displaystyle
f_{4,j}(\widehat{\bm{\Phi}})+f_{5,j}(\widehat{\bm{\Phi}})+f_{6,j}(\widehat{\bm{\Phi}})\leq\mathbf{b}_{u,j}\,,\,,$
(37)

where

$\displaystyle f_{1,j}(\widehat{\bm{\Phi}})=\frac{v_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}}{1-\epsilon_{\infty}\tau}\,,\quad
f_{4,j}(\widehat{\bm{\Phi}})=\frac{v_{\infty}\left\lVert
F_{u,j}\widehat{\bm{\Phi}}_{uy}\right\rVert_{1}}{1-\epsilon_{\infty}\tau}\,,$
$\displaystyle
f_{2,j}(\widehat{\bm{\Phi}})=w_{\infty}\left\lVert\begin{bmatrix}\left(F_{y,j}\widehat{\bm{\Phi}}_{yu}\right)^{\mathsf{T}}\\\
\epsilon_{\infty}\frac{1+\tau\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}}{1-\epsilon_{\infty}\tau}\left(F_{y,j}\widehat{\bm{\Phi}}_{yy}\right)^{\mathsf{T}}\end{bmatrix}\right\rVert_{1}\,,$
$\displaystyle
f_{5,j}(\widehat{\bm{\Phi}})=w_{\infty}\left\lVert\begin{bmatrix}\left(F_{u,j}\widehat{\bm{\Phi}}_{uu}\right)^{\mathsf{T}}\\\
\epsilon_{\infty}\frac{1+\tau\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}}{1-\epsilon_{\infty}\tau}\left(F_{u,j}\widehat{\bm{\Phi}}_{uy}\right)^{\mathsf{T}}\end{bmatrix}\right\rVert_{1}\,,$
$\displaystyle
f_{3,j}(\widehat{\bm{\Phi}})=F_{y,j}\widehat{\bm{\Phi}}_{yy}\widehat{\bm{y}}_{0}+\epsilon_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}\left(\frac{1+\tau\left\lVert\widehat{\mathbf{y}}_{0}\right\rVert_{\infty}}{1-\epsilon_{\infty}\tau}\right)\,,$
$\displaystyle
f_{6,j}(\widehat{\bm{\Phi}})=F_{u,j}\widehat{\bm{\Phi}}_{uy}\widehat{\bm{y}}_{0}+\epsilon_{\infty}\left\lVert
F_{u,j}\widehat{\bm{\Phi}}_{uy}\right\rVert_{1}\left(\frac{1+\tau\left\lVert\widehat{\mathbf{y}}_{0}\right\rVert_{\infty}}{1-\epsilon_{\infty}\tau}\right)\,,$

then $\widehat{\bm{\Phi}}$ satisfies also the safety constraints (31)-(32) for
the real system.

Lemma 3.8 exploits the upperbound
$\left\lVert\widehat{\bm{\Phi}}_{uy}\right\rVert_{\infty}\leq\tau$ to take
into account the worst-case increase in the values of the input and output
signals due to uncertainty. One can interpret (36)-(37) as a kind of
constraint tightening approach related to the one that is used in the robust
learning-based MPC literature [42]. In our setup, similar to [12], the
feasible set shrinks in the presence of larger impulse and free response
estimation error $\epsilon_{\infty}$. This is because (36)-(37) are more
restrictive, and will eventually become infeasible for sufficiently large
$\epsilon_{\infty}$. Instead, the effect of increasing the value of $\tau$ is
less intuitive. Indeed, as $\tau$ increases, the constraint
$\left\lVert{\widehat{\mathbf{\Phi}}}_{uy}\right\rVert_{\infty}\leq\tau$
softens while (36)-(37) tightens. It is therefore necessary to explicitly
optimize over $\tau$. Further discussion and comparison with the MPC
literature can be found in [12] for the state-feedback case. We are now ready
to establish a reformulation of problem (28).

######  Theorem 3.9.

Assume estimation errors $\epsilon_{2}$ and $\epsilon_{\infty}$. Consider the
following optimization problem:

$\displaystyle\min_{\gamma\in[0,\epsilon_{2}^{-1}),\tau\in[0,\epsilon_{\infty}^{-1})}$
$\displaystyle\frac{1}{1-\epsilon_{2}\gamma}\min_{\widehat{\bm{\Phi}}}\qquad
J_{inner}$ (38) $\displaystyle\operatorname*{subject~{}to}~{}$
$\displaystyle\eqref{eq:constraints_hatPhi1},~{}\eqref{eq:constraints_hatPhi2},~{}\eqref{eq:constraints_hatPhi3}\,,$
$\displaystyle\left\lVert\widehat{\bm{\Phi}}_{uy}\right\rVert_{2}\leq\gamma\,,\quad\left\lVert\widehat{\bm{\Phi}}_{uy}\right\rVert_{\infty}\leq\tau\,,$
(39)
$\displaystyle\eqref{eq:quasiconvex_constraints1}-\eqref{eq:quasiconvex_constraints2},\quad\forall
j=1,\dots,sN\,,$

where $J_{inner}$ is equal to

$\left\lVert\begin{bmatrix}\sqrt{1+h(\epsilon_{2},\gamma,\widehat{\mathbf{G}})+h(\epsilon_{2},\gamma,\widehat{\mathbf{y}}_{0})}\widehat{\bm{\Phi}}_{yy}&\widehat{\bm{\Phi}}_{yu}&\widehat{\bm{\Phi}}_{yy}\widehat{\mathbf{y}}_{0}\\\
\sqrt{1+h(\epsilon_{2},\gamma,\widehat{\mathbf{y}}_{0})}\widehat{\bm{\Phi}}_{uy}&\widehat{\bm{\Phi}}_{uu}&\widehat{\bm{\Phi}}_{uy}\widehat{\mathbf{y}}_{0}\end{bmatrix}\right\rVert_{F}.$

Then, (38) has the following properties:

  1. $i)$

it can be written as

$\min_{\gamma\in[0,\epsilon_{2}^{-1}),\tau\in[0,\epsilon_{\infty}^{-1})}J_{robust}(\gamma,\tau)\,,$
(40)

where the value $J_{robust}(\gamma,\tau)$ is obtained by solving a convex
program in $\widehat{\bm{\Phi}}$,

  2. $ii)$

all of its feasible solutions yield a controller
$\widehat{\mathbf{K}}=\widehat{\bm{\Phi}}_{uy}\widehat{\bm{\Phi}}_{yy}^{-1}$
complying with the safety constraints (16)-(17) for the real system,

  3. $iii)$

its minimal cost upper-bounds that of (27).

######  Proof 3.10.

Lemma 3.7 shows that the the cost of (38) upper-bounds
$J(\mathbf{G},\mathbf{K})=J(\mathbf{G},\widehat{\bm{\Phi}}_{uy}\widehat{\bm{\Phi}}_{yy}^{-1})$
for every feasible $\mathbf{K}$. Lemma 3.8 shows that (36)-(37) imply
(31)-(32), which are equivalent to the safety constraints (16)-(17) for the
real system. By inspection of (38) it is clear that it can be rewritten as per
(40). When $\gamma$ and $\tau$ are fixed, the cost $J_{inner}$ is convex in
$\widehat{\bm{\Phi}}$ and so are the constraints of the inner optimization
problem.

######  Remark 3.11 (Safe exploration).

In many applications it is desirable not only that the control policy
synthesized from data is safe, but also that the system operates safely during
the data-collection phase. In our setup this amounts to requiring that
historical and recent trajectories verify (4). Similar to [12], we now clarify
that safe data can be collected exploiting the feasible space of the
optimization problem (38) as a corollary of Theorem 3.9.

More in details, assume that rough estimates
$\widehat{\mathbf{G}}_{r},\widehat{\mathbf{y}}_{0,r}$ are given with possibly
large errors $\epsilon_{\infty,r}$ and $\epsilon_{2,r}$. Note that it is
inherently impossible to guarantee safe exploration unless some prior
information is available. Consider now an exploration signal $\eta(t)$ such
that $\left\lVert\eta(t)\right\rVert_{\infty}\leq\eta_{\infty}$ for every
$t\in\mathbb{Z}$ and define $\widetilde{w}_{\infty}=w_{\infty}+\eta_{\infty}$.
Let $\mathbf{K}_{r}=\widehat{\bm{\Phi}}_{uy,r}\widehat{\bm{\Phi}}_{yy,r}^{-1}$
be any feasible solution to the instance of the optimization problem (38)
where we use $\widetilde{w}_{\infty}$ in place of $w_{\infty}$. Then, by
Theorem 3.9, the control policy

$\mathbf{u}=\mathbf{K}_{r}\mathbf{y}+\bm{\eta}+\mathbf{w}\,,$

can be applied to the real system during the exploration phase to generate
safe recent and historical trajectories. Last, we note that if one is
interested in generating arbitrarily long safe trajectories, it is possible to
establish an infinite-horizon version of problem (38) by adding a tail
variable and adopting a standard finite-horizon approximation of stable
transfer functions; we refer the interested reader to [12] for details.

###  3.4 Numerical Implementation

Theorem 3.9 shows that problem (28), which is highly non-convex in its matrix
variables, can be approximated as the problem of minimizing
$J_{robust}(\gamma,\tau)$ in (40) over two scalar variables $\gamma$ and
$\tau$. The $\epsilon$-dependent suboptimality introduced by such an
approximation will be quantified in the next section. The global optimum of
(38) can thus be determined by exhaustive search over the box
$(\gamma,\tau)\in[0,\epsilon_{2}^{-1})\times[0,\epsilon_{\infty}^{-1})$, for
instance through gridding, random search [43] or bisection [44]. The value
$J_{robust}(\gamma,\tau)$ needs to be computed at each location by fixing
$\gamma$ and $\tau$ and solving a convex program in $\widehat{\bm{\Phi}}$.
While being enormously more tractable than trying to tackle (28) directly, the
gridding over $(\tau,\gamma)$ suffers from significant computational
complexity if we are interested in determining a near-optimal solution with
very low tolerance.

In this section we offer novel insight on the properties of
$J_{robust}(\gamma,\tau)$ and discuss a more efficient numerical
implementation. First, as observed in [7, 31], the cost function in problem
(38) can be made independent of $\gamma$ by introducing a new hyper-parameter
$\alpha\in\mathbb{R}$ that acts as an upper-bound to $\gamma$. As a result,
the cost $J_{robust}(\gamma,\tau)$ becomes unimodal in $\gamma$, and the
globally optimal $\gamma^{\star}(\tau)$ for each fixed $\tau$ can be found
efficiently through golden-section search [45].

######  Proposition 3.12.

Assume estimation errors $\epsilon_{2}$, $\epsilon_{\infty}$ and fix
$\alpha\in[0,\epsilon_{2}^{-1})$. Consider the following optimization problem

$\displaystyle\min_{\gamma\in[0,\alpha],\tau\in[0,\epsilon_{\infty}^{-1})}$
$\displaystyle\frac{1}{1-\epsilon_{2}\gamma}\min_{\widehat{\bm{\Phi}}}\qquad
J_{inner}^{\alpha}(\widehat{\bm{\Phi}})$ (41)
$\displaystyle\operatorname*{subject~{}to}~{}$
$\displaystyle\eqref{eq:constraints_hatPhi1},~{}\eqref{eq:constraints_hatPhi2},~{}\eqref{eq:constraints_hatPhi3},~{}\eqref{eq:quasiconvex_constraints1},~{}\eqref{eq:quasiconvex_constraints2},~{}\eqref{eq:norm_constraints_quasiconvex}$
$\displaystyle\forall j=1,\dots,sN\,,$

where $J_{inner}^{\alpha}(\widehat{\bm{\Phi}})$ is defined as

$\left\lVert\begin{bmatrix}\sqrt{1+h(\epsilon_{2},\alpha,\widehat{\mathbf{G}})+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})}\widehat{\bm{\Phi}}_{yy}&\widehat{\bm{\Phi}}_{yu}&\widehat{\bm{\Phi}}_{yy}\widehat{\mathbf{y}}_{0}\\\
\sqrt{1+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})}\widehat{\bm{\Phi}}_{uy}&\widehat{\bm{\Phi}}_{uu}&\widehat{\bm{\Phi}}_{uy}\widehat{\mathbf{y}}_{0}\end{bmatrix}\right\rVert_{F}.$

Then, the statements $i)$, $ii)$ and $iii)$ of Theorem 3.9 hold. Furthermore

  * $iv)$

$J_{robust}(\gamma,\tau)$ is unimodal in $\gamma$.

######  Proof 3.13.

Since $\gamma\leq\alpha$, $\alpha<\epsilon_{2}^{-1}$, and
$h(\epsilon,\gamma,\cdot)$ is a monotonously increasing function of $\gamma$,
then the inequality $\eqref{eq:nonconvex_bound}$ in Lemma 3.7 continues to
hold when putting $\alpha$ in place of $\gamma$ inside the $h(\cdot)$
functionals. The constraints of (38) are unaffected. Hence, $i)$, $ii)$ and
$iii)$ of Theorem 3.9 continue to hold. It remains to prove $iv)$. Let us fix
any value for $\tau$. First, notice that
$J_{inner}^{\alpha}(\widehat{\bm{\Phi}})$ is a convex function of
$\widehat{\bm{\Phi}}$ and does not depend on $\gamma$, and that the feasible
set of the inner minimization in problem (41) is convex. Denote as $g(\gamma)$
the optimal value of the inner optimization problem. We are left with
minimizing the functional $\frac{g(\gamma)}{1-\epsilon_{2}\gamma}$ over
$\gamma$. We know that $g(\gamma)$ is convex in $\gamma$ because it is
obtained as the partial minimization of a convex functional over a convex set
[46], and that $(1-\epsilon_{2}\gamma)$ is concave in $\gamma$. Since the
ratio of a non-negative convex function and a positive concave function is
quasiconvex555A function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is
quasiconvex if and only if $f(\theta
x_{1}+(1-\theta)x_{2})\leq\max(f(x_{1}),f(x_{2}))$ for every
$x_{1},x_{2}\in\mathbb{R}^{n}$ and every $\theta\in[0,1]$. We refer to [47]
for a comprehensive discussion., we conclude that $J_{robust}(\gamma,\tau)$ is
unimodal in $\gamma$.

It is then natural to wonder whether the function $J_{robust}(\gamma,\tau)$ is
also unimodal in $\tau$. By inspection of the safety constraints (36)-(37),
which appear in a similar form as the state-feedback ones in [12], this may
not in general be the case. Indeed, the product of a strictly monotonically
increasing function and a convex function is not jointly quasiconvex in
general. For instance, similar to the constraints (36)-(37) and those of [12],
consider the function $s:\mathbb{R}^{2}\rightarrow\mathbb{R}$ defined as
$s(x,y)=\frac{|y|x}{(1-x)}$. Fixing
$x_{1}=1,x_{2}=0,y_{1}=-\frac{1}{2},y_{2}=\frac{1}{2}$, one can verify that

$s(\theta x_{1}+(1-\theta)x_{2},\theta
y_{1}+(1-\theta)y_{2})=|0.5-\theta|\frac{\theta}{1-\theta}\,,$

is not unimodal for $\theta\in[0,1]$. Based on this reasoning,
$J_{robust}(\gamma,\tau)$, and similarly the objective (2.3) in [12], may not
be unimodal in $\tau$. We conclude with two observations on unimodality in
$\tau$.

  1. 1.

Similar to $\alpha$, it is possible to introduce a fixed hyper-parameter
$\beta\in[0,\epsilon_{\infty}^{-1})$ that upperbounds $\tau$ and slightly
modifies the definitions of the functions $f_{i,j}(\cdot)$ for $i=1,\dots,6$
to express (36)-(37) as
$\frac{g_{y,j}(\widehat{\bm{\Phi}},\beta)}{1-\epsilon_{\infty}\tau}\leq\mathbf{b}_{y}$
and
$\frac{g_{u,j}(\widehat{\bm{\Phi}},\beta)}{1-\epsilon_{\infty}\tau}\leq\mathbf{b}_{u,j}$.
Since the constraints in this new form are jointly quasiconvex in
$(\tau,\widehat{\bm{\Phi}})$, $J_{robust}(\tau,\gamma)$ becomes unimodal in
$\tau$.

  2. 2.

Exhaustive search over $\tau$ is desirable even if $J_{robust}(\gamma,\tau)$
is separately unimodal in $\gamma$ and $\tau$ as per the point above. This is
because there may be multiple local minima despite separate unimodality.
Hence, we can enjoy the benefits of golden-section search by only introducing
the hyper-parameter $\alpha$.

Motivated as above, the next section will study the suboptimality incurred by
(41).

##  4 Suboptimality Analysis

In this section, we tackle question Q2) in Section 2 about performance
degradation as a function of the level of noise corrupting the data. We denote
as $\mathbf{K}^{\star},\bm{\Phi}^{\star}$ the optimal controller and
corresponding closed-loop responses for the real constrained problem (19).
Furthermore, we denote as
$\widehat{\mathbf{K}}^{\star},\widehat{\bm{\Phi}}^{\star}$ the optimal
controller and corresponding closed-loop responses for program (41) and we let
$J^{\star}=J(\mathbf{G},\mathbf{K}^{\star})$ and
$\hat{J}=J(\mathbf{G},\widehat{\mathbf{K}}^{\star})$. We aim to characterize
the relative suboptimality gap
$\frac{\hat{J}^{2}-{J^{\star}}^{2}}{{J^{\star}}^{2}}$, and specifically we
will show that

$\frac{\hat{J}^{2}-{J^{\star}}^{2}}{{J^{\star}}^{2}}\leq\mathcal{O}\left(\epsilon_{2}\right)+\tilde{S}(\epsilon_{\infty},\epsilon_{2})\,,$

where
$\tilde{S}(\epsilon_{\infty},\epsilon_{2})=S(\epsilon_{\infty})(1+\mathcal{O}(\epsilon_{2}))$,
$S(\epsilon_{\infty})$ quantifies the suboptimality incurred by tightening the
constraints and is such that $S(0)=0$. We prove that if $\epsilon_{2}$ and
$\epsilon_{\infty}$ are small enough and the optimal controller
$\mathbf{K}^{\star}$ does not activate the safety constraints, the
suboptimality gap decrease _linearly_ with $\epsilon_{2}$ and
$S(\epsilon_{\infty})=0$. Otherwise, the gap may decrease according to
$S(\epsilon_{\infty})$, for which we provide a numerical plot in Section 5. In
other words, for small estimation errors $\epsilon_{2}$ and
$\epsilon_{\infty}$ on the impulse response, applying controller
$\widehat{\mathbf{K}}^{\star}$ (which is solely computed from noisy data) to
the _real_ plant achieves almost optimal closed-loop performance while
guaranteeing compliance with safety constraints. Surprisingly, despite the
additional complexity of output-feedback and output noise, our bound matches
the scaling with respect to $\epsilon$ that has been derived in [12] for the
state-feedback case without measurement noise.

We characterize a feasible solution to problem (41), which we later exploit to
establish our suboptimality bound. The proof of Lemma 4.14 and Theorem 4.15 is
reported in the Appendices .5 and .6, respectively.

######  Lemma 4.14 (Feasible solution).

Let $\eta=\epsilon_{2}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}$ and
$\zeta=\epsilon_{\infty}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{\infty}$.
Assume that the estimation errors are small enough to guarantee
$\eta<\frac{1}{5}$ and $\zeta<\frac{1}{2}$, and assume
$\alpha\in[\sqrt{2}\frac{\eta}{\epsilon_{2}(1-\eta)},\epsilon_{2}^{-1})$.
Consider the following optimization problem and its optimal solutions
$\bm{\Phi}^{c}$:

$\displaystyle\bm{\Phi}^{c}\in$
$\displaystyle\operatorname*{arg\,min}_{\bm{\Phi}}$
$\displaystyle\left\lVert\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}&\bm{\Phi}_{yy}\mathbf{y}_{0}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}&\bm{\Phi}_{uy}\mathbf{y}_{0}\end{bmatrix}\right\rVert_{F}$
(42) $\displaystyle\operatorname*{subject~{}to}$
$\displaystyle\begin{bmatrix}I&-\mathbf{G}\end{bmatrix}\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}\end{bmatrix}=\begin{bmatrix}I&0\end{bmatrix},$
$\displaystyle\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}\end{bmatrix}\begin{bmatrix}-\mathbf{G}\\\
I\end{bmatrix}=\begin{bmatrix}0\\\ I\end{bmatrix},$
$\displaystyle\left\lVert\bm{\Phi}_{uy}\right\rVert_{2}\leq\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}\,,\left\lVert\bm{\Phi}_{uy}\right\rVert_{\infty}\leq\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{\infty}\,,$
$\displaystyle\phi_{1,j}(\bm{\Phi})+\phi_{2,j}(\bm{\Phi})+\phi_{3,j}(\bm{\Phi})\leq\mathbf{b}_{y,j}\,,$
(43)
$\displaystyle\phi_{4,j}(\bm{\Phi})+\phi_{5,j}(\bm{\Phi})+\phi_{6,j}(\bm{\Phi})\leq\mathbf{b}_{u,j}\,,$
(44) $\displaystyle\forall j=1,\dots,sN\,,$
$\displaystyle\bm{\Phi}_{yy},\bm{\Phi}_{yu},\bm{\Phi}_{uy},\bm{\Phi}_{uu}\text{\emph{
with causal sparsities}}\,.$

where

$\displaystyle\phi_{1,j}(\bm{\Phi})=\frac{v_{\infty}\left\lVert
F_{y,j}\bm{\Phi}_{yy}\right\rVert_{1}}{1-2\zeta}\,,\quad\phi_{4,j}(\bm{\Phi})=\frac{v_{\infty}\left\lVert
F_{u,j}\bm{\Phi}_{uy}\right\rVert_{1}}{1-2\zeta}\,,$
$\displaystyle\phi_{2,j}(\bm{\Phi})=w_{\infty}\left\lVert\begin{bmatrix}\left(F_{y,j}\bm{\Phi}_{yu}\right)^{\mathsf{T}}\\\
2\frac{\epsilon_{\infty}+\zeta\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}}{1-2\zeta}\left(F_{y,j}\bm{\Phi}_{yy}\right)^{\mathsf{T}}\end{bmatrix}\right\rVert_{1}\,,$
$\displaystyle\phi_{5,j}(\bm{\Phi})=w_{\infty}\left\lVert\begin{bmatrix}\left(F_{u,j}\bm{\Phi}_{uu}\right)^{\mathsf{T}}\\\
2\frac{\epsilon_{\infty}+\zeta\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}}{1-2\zeta}\left(F_{u,j}\bm{\Phi}_{uy}\right)^{\mathsf{T}}\end{bmatrix}\right\rVert_{1}\,,$
$\displaystyle\phi_{3,j}(\bm{\Phi})=F_{y,j}\bm{\Phi}_{yy}\widehat{\mathbf{y}}_{0}+2\frac{\epsilon_{\infty}+\zeta\left\lVert\widehat{\mathbf{y}}_{0}\right\rVert_{\infty}}{1-2\zeta}\left\lVert
F_{y,j}\bm{\Phi}_{yy}\right\rVert_{1}\,,$
$\displaystyle\phi_{6,j}(\bm{\Phi})=F_{u,j}\bm{\Phi}_{uy}\widehat{\mathbf{y}}_{0}+2\frac{\epsilon_{\infty}+\zeta\left\lVert\widehat{\mathbf{y}}_{0}\right\rVert_{\infty}}{1-2\zeta}\left\lVert
F_{u,j}\bm{\Phi}_{uy}\right\rVert_{1}\,.$

Then, the following expressions

$\displaystyle\widetilde{\bm{\Phi}}_{yy}$
$\displaystyle=\bm{\Phi}^{c}_{yy}(I+\mathbf{\Delta}\bm{\Phi}^{c}_{uy})^{-1},~{}\widetilde{\bm{\Phi}}_{yu}=\bm{\Phi}^{c}_{yy}(I+\mathbf{\Delta}\bm{\Phi}^{c}_{uy})^{-1}(\mathbf{G}-\mathbf{\Delta}),$
$\displaystyle\widetilde{\bm{\Phi}}_{uy}$
$\displaystyle=\bm{\Phi}^{c}_{uy}(I+\mathbf{\Delta}\bm{\Phi}^{c}_{uy})^{-1},~{}\widetilde{\bm{\Phi}}_{uu}=(I+\bm{\Phi}^{c}_{uy}\mathbf{\Delta})^{-1}\bm{\Phi}^{c}_{uu},$
$\displaystyle\widetilde{\gamma}$
$\displaystyle=\frac{\sqrt{2}\eta}{\epsilon_{2}(1-\eta)},\qquad\widetilde{\tau}=\frac{\zeta}{\epsilon_{\infty}(1-\zeta)}\,,$
(45)

provide a feasible solution to problem (41).

Since the constraints of (41) are more stringent than those of (19), one
cannot expect, in general, the optimal solution $\bm{\Phi}^{\star}$ to (19) to
be feasible for (41). Hence, to characterize a feasible solution to (41), in
(42) we consider constraints that may even be more stringent than those of
(41). Clearly, the optimal solution
$\mathbf{K}^{c}=\bm{\Phi}_{uy}^{c}(\bm{\Phi}_{yy}^{c})^{-1}$ to (42) will
yield a suboptimal cost $J(\mathbf{G},\mathbf{K}^{c})\geq
J(\mathbf{G},\mathbf{K}^{\star})$. We denote the corresponding suboptimality
gap as

$S(\epsilon_{\infty})=\frac{J(\mathbf{G},\mathbf{K}^{c})^{2}-J(\mathbf{G,\mathbf{K}^{\star}})^{2}}{J(\mathbf{G,\mathbf{K}^{\star}})^{2}}\,.$
(46)

Note that, if the estimation error $\epsilon_{\infty}$ is too large, the
program (42) may become infeasible. This is expected as the uncertainty level
might be incompatible with the required safety. On the other hand, if the
optimal solution to the non-noisy problem (22) does not activate the safety
constraints, then the constraints of (42) remain inactive for small enough
$\epsilon_{\infty}$. In such case we have that $S(\epsilon_{\infty})=0$.

We are now ready to state the main suboptimality result.

######  Theorem 4.15.

Let $\eta=\epsilon_{2}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}$ and
$\zeta=\epsilon_{\infty}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{\infty}$.
Assume that the estimation errors are small enough to guarantee
$\eta<\frac{1}{5}$ and $\zeta<\frac{1}{2}$, and assume that
$\alpha\in\left[\sqrt{2}\frac{\eta}{\epsilon_{2}(1-\eta)},5\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}\right]$.
Moreover, assume that $\epsilon_{\infty}$ is small enough to make the
optimization problem (42) feasible. Then, when applying the controller
$\widehat{\mathbf{K}}^{\star}$ optimizing (38) to the true plant $\mathbf{G}$,
the relative error with respect to the true optimal cost is upper bounded as

$\displaystyle\frac{\hat{J}^{2}-{J^{\star}}^{2}}{{J^{\star}}^{2}}$
$\displaystyle\leq 20\eta+4(M^{c}+V^{c})+4S(\epsilon_{\infty})(1+M^{c}+V^{c})$
$\displaystyle=\mathcal{O}\left(\epsilon_{2}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}(\left\lVert\mathbf{G}\right\rVert^{2}_{2}+\left\lVert\mathbf{y}_{0}\right\rVert^{2}_{2})\right)+$
$\displaystyle~{}~{}+4S(\epsilon_{\infty})(1+M^{c}+V^{c})\,,$ (47)

where

$\displaystyle
M^{c}=h(\epsilon_{2},\alpha,\widehat{\mathbf{G}})+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})+h(\epsilon_{2},\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2},\mathbf{G})$
$\displaystyle\qquad+h(\epsilon_{2},\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2},\mathbf{y}_{0})\,,$
$\displaystyle
V^{c}=h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})+h(\epsilon_{2},\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2},\mathbf{y}_{0})\,.$

We have expressed the suboptimality gap in the form (47) to highlight the
presence of two main parts; the first addend scales as
$\mathcal{O}\left(\epsilon_{2}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}(\left\lVert\mathbf{G}\right\rVert^{2}_{2}+\left\lVert\mathbf{y}_{0}\right\rVert^{2}_{2})\right)$
and the second addend $S(\epsilon_{\infty})(1+M^{c}+V^{c})$ is linked to the
suboptimality of the tightened program (42). The most important observation is
that the suboptimality may decrease _linearly_ with $\epsilon_{2}$ when
$\max{(\epsilon_{2},\epsilon_{\infty})}$ is small enough. This has also been
observed for the unconstrained case in [7, 31]. Here, while the term
$S(\epsilon_{\infty})$ may affect the linearity of the scaling, we observe
through examples (cfr. Figure 1(b)) that $S(\epsilon_{\infty})$ sharply
transitions from $0$ to $\infty$ as $\epsilon_{\infty}$ increases. Hence, in
practice, $S(\epsilon_{\infty})$ can be interpreted as a boolean value
indicating whether $\epsilon_{\infty}$ is small enough for our bound to hold.

Further, our suboptimality bound (47) indicates qualities of the underlying
unknown system that make it easier to be safely controlled based on noisy
data. Notably, the suboptimality grows quadratically with the norm of the true
impulse and free responses. This fact implies that an unknown unstable system
will be difficult to control for a long horizon. Last, we note that,
surprisingly, our rate in terms of $\epsilon_{2}$ matches that of [12] which
was valid under the assumption of exact state measurements. In other words,
our analysis shows that near-optimality can be ensured in complex data-driven
control scenarios that combine hard safety requirements with noisy output
measurements.

######  Remark 4.16 (Sample complexity).

In related work, e.g. [4, 12, 7, 29], the authors more precisely quantify
$\epsilon_{2}$, $\epsilon_{\infty}$ and the probability of the estimate to be
within the corresponding norm error intervals as a function of the noise
statistics and the real system parameters. As a result, an end-to-end sample
complexity analysis is obtained. These works achieve such result by focusing
on a specific estimation technique (e.g. least squares system identification
in [4, 7] and Hankel matrix column-averaging in [29]) and the corresponding
non-asymptotic norm error bounds available in the literature. We expect that
analogous results can be derived for the LS behavioral choice
$(G,g)=(G_{LS},g_{LS})$. On the contrary, in this work we put the focus on
approximating the impulse and free responses without binding ourselves to a
specific estimation technique. Hence, we will not further analyze
$\epsilon_{2}$, $\epsilon_{\infty}$ and the success probability, as both are
dependent on the chosen estimation technique. Instead, we provide numerical
examples in Section 5 to compare LS and ML estimators.

##  5 Numerical Experiments

(a) Closed-loop trajectories. The grey region indicates unsafe input and
output values. The green and blue regions contain trajectories for $50$ noise
realizations obtained through $\mathbf{K}^{\star}$ and
$\widehat{\mathbf{K}}^{\star}$, respectively. Green and blue lines represent a
specific trajectory in both settings.

(b) Robust suboptimality gap $S(\epsilon_{\infty})$. This quantity can be
interpreted as an indicator as to whether the guarantee (47) holds for a given
$\epsilon_{\infty}$.

Figure 1: Safe controller synthesis for system (48).

In this section, we demonstrate numerically the effectiveness of the proposed
framework in safely controlling unknown systems. In the experiments, we
consider the single-input single-output unknown LTI system characterized by
the matrices

$\displaystyle A=\rho\begin{bmatrix}1&0.25\\\ 0&1\end{bmatrix},\text{
}B=\begin{bmatrix}0\\\ 0.1\end{bmatrix},\text{
}C=\begin{bmatrix}1&-1\end{bmatrix}\,,$ (48)

where $\rho>0$ corresponds to the spectral radius of $A$. When $\rho<1$, (48)
is asymptotically stable, that is, its output converges to the origin at an
exponential rate when the input is equal to $0$. When $\rho=1$, (48) is an
unstable double-integrator system.

In all the following tests, the cost function is given by (3) for appropriate
choices of the weights. The expectation in (3) is taken over future
input/output noise with variances $\Sigma_{w}=I_{m}$ and $\Sigma_{v}=I_{p}$.
We consider bounded disturbances between $-1$ and $1$, that is,
$w_{\infty}=v_{\infty}=1$. For solving optimization problems we use MOSEK
[48], called through MATLAB via YALMIP [49] on a standard laptop
computer666The code is open-source and available at https://gitlab.nccr-
automation.ch/data-driven-control-epfl/constrained-biop..

###  5.1 Example: safe controller synthesis from noisy data

In our first test, we synthesize a safe output-feedback controller for system
(48) with $\rho=1$ from noisy data. We assume that
$x_{0}=x(1)=\begin{bmatrix}6&0\end{bmatrix}^{\mathsf{T}}$. The goal is to
comply with the safety constraints

$\displaystyle-5.5$ $\displaystyle\leq y(t)\leq 5.5,\quad\forall
t=2,\dots,12\,,$ $\displaystyle-100$ $\displaystyle\leq u(t)\leq
100\,,\quad\forall t=1,\dots,11\,,$

for all realizations of noise
$\|\mathbf{w}\|_{\infty},\|\mathbf{v}\|_{\infty}\leq 1$, while minimizing the
cost (3) with the weights $\mathbf{Q}$ and $\mathbf{R}$ in (19) set to the
identity.

We first synthesize the optimal controller assuming that the historical and
recent data are not affected by noise. To this end, we cast and solve the
convex optimization problem (22). We verify that the optimal controller
$\mathbf{K}^{\star}$ yields a cost $J(\mathbf{G},\mathbf{K}^{\star})=69.88$.
The green tubes in Figure 1(a) show the regions containing 50 realizations of
the optimal closed-loop input and output trajectories. Due to the high level
of noise, we can observe a significant variability in the trajectory values
for different noise realizations. Nonetheless, all trajectories are safe.

We then discuss the case where historical and recent data are both affected by
noise. We consider a noisy estimate
$(\widehat{\mathbf{G}},\widehat{\mathbf{y}}_{0})$ such that
$\epsilon=\max(\epsilon_{2},\epsilon_{\infty})=0.01$ and compute a near-
optimal solution to the proposed optimization problem (38). As discussed in
Section 3.4, this can be achieved by 1) extensive search over $\gamma$ and
$\tau$, 2) extensive search over $\tau$ and golden-section search over
$\gamma$, or 3) golden-section search over $\gamma$ and $\tau$. Even if the
first solution comes without strong theoretical guarantees, extensive search
over $\gamma$ and $\tau$ may be simpler to implement as it avoids the delicate
task of tuning the hyper-parameters $\alpha$ and $\beta$. Specifically, for
this example we have searched over $100$ randomly extracted values of $\gamma$
and $\tau$ in the interval $[0,\epsilon^{-1})$. A potential improvement to
this heuristic could be to use a bisection algorithm, as proposed in [44] for
example.

Proceeding as above, we synthesize a robustly safe controller
$\widehat{\mathbf{K}}^{\star}$ yielding a cost of
$J(\mathbf{G},\widehat{\mathbf{K}}^{\star})=140.54$. The corresponding
suboptimality gap is
$\frac{\hat{J}^{2}-{J^{\star}}^{2}}{{J^{\star}}^{2}}=3.049$. In Figure 1(a),
the trajectories and variability levels resulting from
$\widehat{\mathbf{K}}^{\star}$ for $50$ noise realizations are plotted in
blue. We observe that, since $\widehat{\mathbf{K}}^{\star}$ is synthesized
using noise-corrupted data, it leads to safer, but more conservative
trajectories. Indeed, due to uncertainty, higher control effort is spent to
keep the output further from the constraints.

It is informative to inspect the robust suboptimality gap
$S(\epsilon_{\infty})$ incurred by the tightened optimization problem (42)
that we have used in the analysis to characterize a feasible solution to (38).
In Figure 1(b), we plot $S(\epsilon_{\infty})$ assuming
$x_{0}=x(1)=\begin{bmatrix}1&0\end{bmatrix}^{\mathsf{T}}$ and requiring
$-3\leq y(t)\leq 3$ for $t=2,\dots,7$. The example exhibits a fast transition
from infeasibility for $\epsilon_{\infty}>0.118$ to near-optimality for
$\epsilon<0.115$. This fact leads to the following observation: high-
performing safe controllers can be synthesized by solving (38) even for values
of $\epsilon_{\infty}$ such that $S(\epsilon_{\infty})=\infty$. However, in
such cases the theoretical suboptimality guarantee (47) may not hold anymore.
This phenomenon is consistent with the numerical examples of [12] for the
state-feedback case.

(a) Estimation error in function of the corrupting noise. ML estimation
through the SMM yields significantly smaller errors than LS. The green and
blue regions indicate the gap for the $2$-norm and $\infty$-norm,
respectively.

(b) Suboptimality gap as a function of $\epsilon_{2}$ (obtained through SMM
estimation) for increasing values of the spectral radius $\rho$ of matrix $A$
(on the left). Suboptimality gap as a function of $\rho$ for increasing values
of $\sigma$ (on the right).

Figure 2: Examples for estimation and suboptimality scaling.

###  5.2 Example: suboptimality scaling beyond least-squares estimation

The bound (47) in Theorem 4.15 states that a low estimation error level
$\epsilon$ is crucial in ensuring safety and near-optimality when controlling
unknown systems based on noisy data. One advantage of the proposed BIOP
formulation is that it is directly compatible with behavioral estimation
approaches beyond LS identification for the reconstruction of the impulse and
free responses, such as data-enabled Kalman filtering [26] and SMM [23, 22].
In our last test, we 1) validate the scaling of the suboptimality gap (47),
and 2) showcase that, for instance, SMM-based estimation [22] may lead to
significantly lower error levels given the same amount of data. We consider
the system (48) with different values of $\rho\in[0.9,0.93,0.96,0.99,1]$ and
$x_{0}=x(1)=\begin{bmatrix}6&0\end{bmatrix}^{\mathsf{T}}$, over a time-horizon
of length $N=11$. The cost function weights in (19) are selected as
$Q(t)=I_{p}$ for every $t=1,\dots,11$, $Q(12)=20I_{p}$ and $R(t)=0.05$ for
every $t=1,\dots,11$. Since the addend $S(\epsilon_{\infty})$ in (47) may
modify the scaling rates, for this example we do not include safety
requirements.

####  5.2.1 Behavioral estimation: LS vs SMM

For a fixed value of $\rho$, we gather noisy historical and recent
trajectories with $T_{h}=200$ and $T_{ini}=30$ respectively. We assume that
the data are corrupted by Gaussian noise with zero mean and variances
$\bm{\Sigma}^{h}_{w}=\bm{\Sigma}^{r}_{w}=\sigma I_{m}$ and
$\bm{\Sigma}^{h}_{v}=\bm{\Sigma}^{r}_{v}=\sigma I_{p}$, where $\sigma\geq 0$.
For each experiment, we fix the variance $\sigma\geq 0$ and select a random
exploration control input $\mathbf{u}$. We collect $1000$ different historical
and recent trajectories for different realizations of the corrupting noise.
For each realization of the trajectories, we compute 1) the LS solution
$(G_{LS},g_{LS})$ using (23)-(24) and the corresponding impulse and free
responses $\widetilde{\mathbf{G}}_{LS},\widetilde{\mathbf{y}}_{0,LS}$, and 2)
the ML solution $(G_{ML},g_{ML})$ using (25)-(26) and the corresponding
impulse and free responses
$\widetilde{\mathbf{G}}_{ML},\widetilde{\mathbf{y}}_{0,ML}$. For each
estimation, we determine the incurred error levels $\epsilon_{2,G}$,
$\epsilon_{\infty,G}$, $\epsilon_{2,y}$ and $\epsilon_{\infty,y}$777Since the
real system is unavailable, in practice this can be done using a bootstrap
procedure.. Last, we record the $90$-th percentile of these values both for
SMM and LS estimation.

In Figure 2(a) we compare the values of $\epsilon_{2}$ and $\epsilon_{\infty}$
incurred by both estimation techniques. We observe that SMM may yield
significantly smaller estimation errors than LS identification. While a full
sample-complexity analysis is still unavailable beyond least-squares [4, 32,
29], these examples showcase an advantage in using more sophisticated
estimation techniques for safe data-driven control.

####  5.2.2 Suboptimality scaling

Having exploited ML estimation to construct approximate impulse and free
responses and the corresponding error-levels, we are ready to solve the
optimization problem (41). Since constraints are not present in this example,
(41) can be simplified to the quasiconvex formulation we have proposed in
[31], where the optimization variable $\tau$ is not present. The hyper-
parameter $\alpha$ is tuned empirically in the interval
$\alpha\in[\sqrt{2}\frac{\eta}{\epsilon_{2}(1-\eta)},\epsilon_{2}^{-1})$.888The
value $\eta=\epsilon_{2}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}$ is
unknown in practice because $\bm{\Phi}^{\star}_{uy}$ is unavailable. One can
then tune $\alpha$ according to $\alpha<\epsilon_{2}^{-1}$.

Figure 2(b) shows the suboptimality gap one incurs by applying the controller
$\widehat{\mathbf{K}}^{\star}$ obtained through the proposed robust BIOP. On
the left, we consider increasing levels of the estimation error level
$\epsilon_{2}$ for each choice of the spectral radius $\rho=0.9,0.99,1$. On
the right, we conversely consider increasing levels of the spectral radius for
each choice of the estimation error level $\epsilon_{2}$. In both cases, we
plot the suboptimality gap
$\frac{\hat{J}^{2}-{J^{\star}}^{2}}{{J^{\star}}^{2}}$. It can be observed
that, consistently with Theorem 4.15, 1) the gap linearly converges to $0$ as
$\epsilon_{2}$ converges to $0$, and 2) the gap grows faster than linearly
with the spectral radius $\rho$. We finally remark that, in finite-horizon,
our BIOP formulations are valid for unstable systems with $\rho>1$. However,
it is inherently challenging to collect trajectories of an unstable system, as
the values to be plugged into the corresponding numerical programs will become
too large to be handled by numerical solvers. For unstable systems in a data-
driven scenario, it is common to assume knowledge of a pre-stabilizing
controller [8, 7].

##  6 Conclusions

In this paper we have addressed two open questions on data-driven predictive
control with output measurements, when the data are corrupted by noise. First,
we have established a data-driven optimization problem that exploits
behavioral theory to yield a robustly safe and near-optimal control policy.
Second we have characterized how much the level of noise impacts the
performance and the safety of an output-feedback control system in the form of
a noise-dependent suboptimality bound.

Our results are intended as a first step towards quantifying the effect of
noise-corrupted data when using behavioral models for complex data-driven
prediction and control tasks. Future work includes studying the performance of
the suggested method when moving to a receding-horizon implementation with
closed-loop predictions. It is also interesting to investigate whether it is
preferable to avoid the behavioral impulse-response identification step and
instead optimize directly over the data, and to quantify the sample-complexity
when using state-of-the-art estimators.

## Acknowledgments

We thank Sarah Dean for sharing the implementation of the examples in [12] and
for helpful insights on unimodality.

## References

  * [1] F. Lamnabhi-Lagarrigue, A. Annaswamy, S. Engell, A. Isaksson, P. Khargonekar, R. M. Murray, H. Nijmeijer, T. Samad, D. Tilbury, and P. Van den Hof, “Systems & control for the future of humanity, research agenda: Current and future roles, impact and grand challenges,” _Annual Reviews in Control_ , vol. 43, pp. 1–64, 2017.
  * [2] G. Baggio, D. S. Bassett, and F. Pasqualetti, “Data-driven control of complex networks,” _Nature communications_ , vol. 12, no. 1, pp. 1–13, 2021.
  * [3] B. Recht, “A tour of reinforcement learning: The view from continuous control,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 2, pp. 253–279, 2019.
  * [4] S. Dean, H. Mania, N. Matni, B. Recht, and S. Tu, “On the sample complexity of the Linear Quadratic Regulator,” _Foundations of Computational Mathematics_ , pp. 1–47, 2019.
  * [5] M. Fazel, R. Ge, S. Kakade, and M. Mesbahi, “Global convergence of policy gradient methods for the linear quadratic regulator,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 1467–1476.
  * [6] D. Malik, A. Pananjady, K. Bhatia, K. Khamaru, P. Bartlett, and M. Wainwright, “Derivative-free methods for policy optimization: Guarantees for linear quadratic systems,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2019, pp. 2916–2925.
  * [7] Y. Zheng and L. Furieri, M. Kamgarpour, and N. Li, “Sample complexity of linear quadratic gaussian (LQG) control for output feedback systems,” _Learning for Dynamics and Control, PMLR_ , 2021.
  * [8] M. Simchowitz, K. Singh, and E. Hazan, “Improper learning for non-stochastic control,” in _Conference on Learning Theory_. PMLR, 2020, pp. 3320–3436.
  * [9] S. Lale, K. Azizzadenesheli, B. Hassibi, and A. Anandkumar, “Logarithmic regret bound in partially observable linear dynamical systems,” _arXiv preprint arXiv:2003.11227_ , 2020.
  * [10] A. Tsiamis, N. Matni, and G. Pappas, “Sample complexity of kalman filtering for unknown systems,” in _Learning for Dynamics and Control_. PMLR, 2020, pp. 435–444.
  * [11] K. Zhang, B. Hu, and T. Basar, “Policy optimization for $\mathcal{H}_{2}$ linear control with $\mathcal{H}_{\infty}$ robustness guarantee: Implicit regularization and global convergence,” in _Learning for Dynamics and Control_. PMLR, 2020, pp. 179–190.
  * [12] S. Dean, S. Tu, N. Matni, and B. Recht, “Safely learning to control the constrained Linear Quadratic Regulator,” in _2019 American Control Conference (ACC)_. IEEE, 2019, pp. 5582–5588.
  * [13] S. Fattahi, N. Matni, and S. Sojoudi, “Efficient learning of distributed linear-quadratic control policies,” _SIAM Journal on Control and Optimization_ , vol. 58, no. 5, pp. 2927–2951, 2020.
  * [14] L. Furieri, Y. Zheng, and M. Kamgarpour, “Learning the globally optimal distributed LQ regulator,” in _Learning for Dynamics and Control_. PMLR, 2020, pp. 287–297.
  * [15] J. C. Willems and J. W. Polderman, _Introduction to mathematical systems theory: a behavioral approach_. Springer Science & Business Media, 1997, vol. 26.
  * [16] J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control: In the shallows of the DeePC,” in _2019 18th European Control Conference (ECC)_. IEEE, 2019, pp. 307–312.
  * [17] ——, “Distributionally robust chance constrained data-enabled predictive control,” _arXiv preprint arXiv:2006.01702_ , 2020.
  * [18] F. Dörfler, J. Coulson, and I. Markovsky, “Bridging direct & indirect data-driven control formulations via regularizations and relaxations,” _arXiv preprint arXiv:2101.01273_ , 2021.
  * [19] F. Fiedler and S. Lucia, “On the relationship between data-enabled predictive control and subspace predictive control,” _arXiv preprint arXiv:2011.13868_ , 2020.
  * [20] C. De Persis and P. Tesi, “Formulas for data-driven control: Stabilization, optimality, and robustness,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 3, pp. 909–924, 2020.
  * [21] J. Berberich, J. Köhler, M. A. Müller, and F. Allgöwer, “Data-driven model predictive control with stability and robustness guarantees,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 4, pp. 1702–1717, 2021.
  * [22] A. Iannelli, M. Yin, and R. S. Smith, “Experiment design for impulse response identification with signal matrix models,” _arXiv preprint arXiv:2012.08126_ , 2020.
  * [23] M. Yin, A. Iannelli, and R. S. Smith, “Maximum likelihood estimation in data-driven modeling and control,” _arXiv preprint arXiv:2011.00925_ , 2020\.
  * [24] Y. Lian and C. N. Jones, “Nonlinear data-enabled prediction and control,” _arXiv preprint arXiv:2101.03187_ , 2021.
  * [25] A. Russo and A. Proutiere, “Poisoning attacks against data-driven control methods,” _arXiv preprint arXiv:2103.06199_ , 2021.
  * [26] D. Alpago, F. Dörfler, and J. Lygeros, “An extended Kalman filter for data-enabled predictive control,” _IEEE Control Systems Letters_ , vol. 4, no. 4, pp. 994–999, 2020.
  * [27] H. J. Vanwaarde, M. K. Camlibel, and M. Mesbahi, “From noisy data to feedback controllers: non-conservative design via a matrix S-lemma,” _IEEE Transactions on Automatic Control, to appear_ , 2020.
  * [28] L. Huang, J. Zhen, J. Lygeros, and F. Dörfler, “Quadratic regularization of data-enabled predictive control: Theory and application to power converter experiments,” _arXiv preprint arXiv:2012.04434_ , 2020.
  * [29] A. Xue and N. Matni, “Data-driven system level synthesis,” _arXiv preprint arXiv:2011.10674_ , 2020.
  * [30] C. De Persis and P. Tesi, “Low-complexity learning of linear quadratic regulators from noisy data,” _Automatica_ , vol. 128, p. 109548, 2021.
  * [31] L. Furieri, B. Guo, A. Martin, and G. Ferrari-Trecate, “A behavioral input-output parametrization of control policies with suboptimality guarantees,” _submitted to the 60th Conference on Decision and Control (CDC)_ , 2021.
  * [32] S. Oymak and N. Ozay, “Non-asymptotic identification of LTI systems from a single trajectory,” in _2019 American Control Conference (ACC)_. IEEE, 2019, pp. 5655–5661.
  * [33] J. H. Gallier, _CIS 515: Vector Norms and Matrix Norms_. University of Pennsylvania, CIS 515 lecture notes, 2020\.
  * [34] 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.
  * [35] L. Furieri and M. Kamgarpour, “Unified approach to convex robust distributed control given arbitrary information structures,” _IEEE Transactions on Automatic Control_ , 2019.
  * [36] A. Bemporad, “Reducing conservativeness in predictive control of constrained systems with disturbances,” in _Proceedings of the 37th IEEE Conference on Decision and Control (Cat. No. 98CH36171)_ , vol. 2. IEEE, 1998, pp. 1384–1389.
  * [37] Y. Zheng, L. Furieri, M. Kamgarpour, and N. Li, “System-level, input-output and new parameterizations of stabilizing controllers, and their numerical computation,” _arXiv preprint arXiv:1909.12346_ , 2019.
  * [38] L. Furieri, Y. Zheng, A. Papachristodoulou, and M. Kamgarpour, “An Input-Output Parametrization of stabilizing controllers: amidst Youla and System Level Synthesis,” _IEEE Control Systems Letters_ , vol. 3, no. 4, pp. 1014–1019, 2019.
  * [39] H. Mania, S. Tu, and B. Recht, “Certainty equivalent control of LQR is efficient,” _arXiv preprint arXiv:1902.07826_ , 2019.
  * [40] J. C. Willems, P. Rapisarda, I. Markovsky, and B. L. De Moor, “A note on persistency of excitation,” _Systems & Control Letters_, vol. 54, no. 4, pp. 325–329, 2005.
  * [41] I. Markovsky and P. Rapisarda, “Data-driven simulation and control,” _International Journal of Control_ , vol. 81, no. 12, pp. 1946–1959, 2008\.
  * [42] M. Zanon and S. Gros, “Safe reinforcement learning using robust MPC,” _IEEE Transactions on Automatic Control_ , 2020.
  * [43] J. Bergstra and Y. Bengio, “Random search for hyper-parameter optimization.” _Journal of machine learning research_ , vol. 13, no. 2, 2012.
  * [44] S. Chen, H. Wang, M. Morari, V. M. Preciado, and N. Matni, “Robust closed-loop model predictive control via system level synthesis,” in _2020 59th IEEE Conference on Decision and Control (CDC)_. IEEE, 2020, pp. 2152–2159.
  * [45] J. Kiefer, “Sequential minimax search for a maximum,” _Proceedings of the American mathematical society_ , vol. 4, no. 3, pp. 502–506, 1953.
  * [46] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [47] A. Agrawal and S. Boyd, “Disciplined quasiconvex programming,” _Optimization Letters_ , pp. 1–15, 2020.
  * [48] MOSEK Aps, “The MOSEK optimization toolbox for MATLAB manual. Version 8.1.” 2017.
  * [49] J. Löfberg, “YALMIP : A Toolbox for Modeling and Optimization in MATLAB,” in _In Proc. of the CACSD Conf._ , Taipei, Taiwan, 2004.

###  .1 Proof of Proposition 1

For the first statement, notice that the controller $\mathbf{K}$ achieves the
closed-loop responses (12). Now select
$(\bm{\Phi}_{yy},\bm{\Phi}_{yu},\bm{\Phi}_{uy},\bm{\Phi}_{uu})$ as

$\begin{bmatrix}\bm{\Phi}_{yy}&\bm{\Phi}_{yu}\\\
\bm{\Phi}_{uy}&\bm{\Phi}_{uu}\end{bmatrix}=\begin{bmatrix}(I-\mathbf{GK})^{-1}&(I-\mathbf{GK})^{-1}\mathbf{G}\\\
\mathbf{K}(I-\mathbf{GK})^{-1}&(I-\mathbf{KG})^{-1}\end{bmatrix}\,,$ (49)

and $\mathbf{q}=\bm{\Phi}_{uu}\mathbf{g}$. Clearly,
$\mathbf{K}=\bm{\Phi}_{uy}\bm{\Phi}_{yy}^{-1}$ and
$\mathbf{g}=\bm{\Phi}_{uu}^{-1}\mathbf{q}$, and by plugging the corresponding
expressions, we verify that (14), (15) and (18) are satisfied. It remains to
prove that (16)-(17) are satisfied. In (10), substitute $\mathbf{y}$ and
$\mathbf{u}$ with their closed-loop expressions (11). It follows that the
addends separately depend on $\mathbf{w}$ or $\mathbf{v}$. Hence, (10) can be
rewritten as

$\displaystyle\max_{\left\lVert\mathbf{v}\right\rVert_{\infty}\leq
v_{\infty}}\left(\mathbf{F}_{y}\bm{\Phi}_{yy}\right)\mathbf{v}+\max_{\left\lVert\mathbf{w}\right\rVert_{\infty}\leq
w_{\infty}}\left(\mathbf{F}_{y}\bm{\Phi}_{yu}\right)\mathbf{w}+$
$\displaystyle+\mathbf{F}_{y}\mathbf{Gq}+\left(\mathbf{F}_{y}\bm{\Phi}_{yy}\right)\mathbf{CP}_{A}(:,0)x_{0}\leq\mathbf{b}_{y}\,,$
(50) $\displaystyle\max_{\left\lVert\mathbf{v}\right\rVert_{\infty}\leq
v_{\infty}}\left(\mathbf{F}_{u}\bm{\Phi}_{uy}\right)\mathbf{v}+\max_{\left\lVert\mathbf{w}\right\rVert_{\infty}\leq
w_{\infty}}\left(\mathbf{F}_{u}\bm{\Phi}_{uu}\right)\mathbf{w}+$
$\displaystyle+\mathbf{F}_{u}\mathbf{q}+\left(\mathbf{F}_{y}\bm{\Phi}_{uy}\right)\mathbf{CP}_{A}(:,0)x_{0}\leq\mathbf{b}_{u}\,,$
(51)

where the $\max(\cdot)$ is to be intended row-wise. The expressions (50)-(51)
are already convex in $\bm{\Phi},\mathbf{q}$. To have a more explicit
expression, similar to [12] we utilize the well-known property that the
$\left\lVert\cdot\right\rVert_{1}$ and the
$\left\lVert\cdot\right\rVert_{\infty}$ vector norms are dual of each other
[46], that is $k\left\lVert x\right\rVert_{1}=\max_{\left\lVert
w\right\rVert_{\infty}\leq k}x^{\mathsf{T}}w$. The result follows immediately
by inspecting (50)-(51) and letting $x^{\mathsf{T}}$ be equal to either
$F_{y,j}\bm{\Phi}_{yy}$, $F_{y,j}\bm{\Phi}_{uy}$, $F_{y,j}\bm{\Phi}_{yu}$ or
$F_{y,j}\bm{\Phi}_{uu}$, and letting $k$ be equal to either $v_{\infty}$ or
$w_{\infty}$.

For the second statement, it is easy to notice that $\mathbf{K}$ is causal by
construction because $\bm{\Phi}_{uy}$ and $\bm{\Phi}_{yy}$ are block lower-
triangular. By selecting the controller
$\mathbf{K}=\bm{\Phi}_{uy}\bm{\Phi}_{yy}^{-1}$ one has

$\displaystyle(I$
$\displaystyle-\mathbf{G}\bm{\Phi}_{uy}\bm{\Phi}_{yy}^{-1})^{-1}=(I-\mathbf{G}\bm{\Phi}_{uy}(I+\mathbf{G}\bm{\Phi}_{uy})^{-1})^{-1}$
$\displaystyle=((I+\mathbf{G}\bm{\Phi}_{uy}-\mathbf{G}\bm{\Phi}_{uy})(I+\mathbf{G}\bm{\Phi}_{uy})^{-1})^{-1}$
$\displaystyle=I+\mathbf{G}\bm{\Phi}_{uy}=\bm{\Phi}_{yy}\,,$

which shows that $\bm{\Phi}_{yy}$ is the closed-loop response from
$\mathbf{v}_{[0,N-1]}+\mathbf{CP}_{A}(:,0)x_{0}$ to $\mathbf{y}_{[0,N-1]}$ as
per (12). Similar computations hold for the remaining closed-loop responses.
For the safety constraints, select any $\bm{\Phi}$ and $\mathbf{q}$ complying
with (16)-(17). It is easy to verify by direct computation that, for any
$\mathbf{w}$ and $\mathbf{v}$, the same input and output trajectories defined
at (11) are obtained by letting $\mathbf{K}=\bm{\Phi}_{uy}\bm{\Phi}_{yy}^{-1}$
and $\mathbf{g}=\bm{\Phi}_{uu}^{-1}\mathbf{q}$ in (6), (7), (8). Hence, the
safety constraints are satisfied for any disturbance realization.

###  .2 Proof of Proposition 3.6

First, we verify by direct inspection that for any $\mathbf{K}$, the
parameters

$\widehat{\bm{\Phi}}=\begin{bmatrix}(I-\widehat{\mathbf{G}}\mathbf{K})^{-1}&(I-\widehat{\mathbf{G}}\mathbf{K})^{-1}\widehat{\mathbf{G}}\\\
\mathbf{K}(I-\widehat{\mathbf{G}}\mathbf{K})^{-1}&(I-\mathbf{K}\widehat{\mathbf{G}})^{-1}\end{bmatrix}\,.$

satisfy the constraints of (28) and are such that
$\mathbf{K}=\widehat{\bm{\Phi}}_{uy}\widehat{\bm{\Phi}}_{yy}^{-1}$. Therefore,
every controller $\mathbf{K}$ is parametrized in problem (28), irrespective of
$\widehat{\mathbf{G}}$. Next, remember that
$\mathbf{G}=\widehat{\mathbf{G}}+\bm{\Delta}$,
$\mathbf{y}_{0}=\widehat{\mathbf{y}}_{0}+\bm{\delta}_{0}$. The real closed-
loop responses $\bm{\Phi}$ achieved by a controller $\mathbf{K}$ can be
expressed in terms of the estimated closed-loop responses
$\widehat{\bm{\Phi}}$ as

$\begin{split}\bm{\Phi}_{yy}&=\left(I-(\widehat{\mathbf{G}}+\mathbf{\Delta})\widehat{\bm{\Phi}}_{uy}\widehat{\bm{\Phi}}_{yy}^{-1}\right)^{-1}=\widehat{\bm{\Phi}}_{yy}\left(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy}\right)^{\text{-}1}\,,\\\
\bm{\Phi}_{yu}&=(I-\mathbf{GK})^{\text{-}1}\mathbf{G}=\bm{\Phi}_{yy}\mathbf{G}=\widehat{\bm{\Phi}}_{yy}\left(I\text{-}\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy}\right)^{\text{-}1}(\widehat{\mathbf{G}}+\mathbf{\Delta})\,,\\\
\bm{\Phi}_{uy}&=\widehat{\bm{\Phi}}_{uy}\widehat{\bm{\Phi}}_{yy}^{-1}\widehat{\bm{\Phi}}_{yy}\left(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy}\right)^{-1}=\widehat{\bm{\Phi}}_{uy}\left(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy}\right)^{-1}\,,\\\
\bm{\Phi}_{uu}&=\bm{\Phi}_{uy}\mathbf{G}+I=\widehat{\bm{\Phi}}_{uy}\left(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy}\right)^{-1}(\widehat{\mathbf{G}}+\mathbf{\Delta})+I\\\
&=\left(I-\widehat{\bm{\Phi}}_{uy}\mathbf{\Delta}\right)^{-1}\widehat{\bm{\Phi}}_{uy}(\widehat{\mathbf{G}}+\mathbf{\Delta})+I\\\
&=\left(I-\widehat{\bm{\Phi}}_{uy}\mathbf{\Delta}\right)^{-1}(\widehat{\bm{\Phi}}_{uy}\widehat{\mathbf{G}}+I)=\left(I-\widehat{\bm{\Phi}}_{uy}\mathbf{\Delta}\right)^{-1}\widehat{\bm{\Phi}}_{uu}\,.\end{split}$

By substituting these expressions into the cost and the safety constraints, we
conclude the proof.

###  .3 Proof of Lemma 3.7

The objective function in Proposition 3.6 can be written as the square-root of
the sum of the square of the Frobenius norms of each of its six blocks. For
the upper-left block, we have

$\displaystyle\|\widehat{\bm{\Phi}}_{yy}(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}\|_{F}\leq\|\widehat{\bm{\Phi}}_{yy}\|_{F}\left\lVert\sum_{k=0}^{\infty}(\bm{\Delta}\widehat{\bm{\Phi}}_{uy})^{k}\right\rVert_{2}$
$\displaystyle\leq\|\widehat{\bm{\Phi}}_{yy}\|_{F}\sum_{k=0}^{\infty}\left\lVert(\epsilon_{2}\widehat{\bm{\Phi}}_{uy})\right\rVert^{k}_{2}=\frac{\|\widehat{\bm{\Phi}}_{yy}\|_{F}}{1-\epsilon_{2}\|\widehat{\bm{\Phi}}_{uy}\|_{2}}\,,$

where the convergence of the series follows from $\bm{\Delta}$ and
$\widehat{\bm{\Phi}}_{uy}$ having zero-entries diagonal blocks by
construction. Similarly,

$\displaystyle\|\widehat{\bm{\Phi}}_{uy}(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}\|_{F}\leq\frac{\|\widehat{\bm{\Phi}}_{uy}\|_{F}}{1-\epsilon_{2}\|\widehat{\bm{\Phi}}_{uy}\|_{2}}\,,$
$\displaystyle\|(I-\widehat{\bm{\Phi}}_{uy}\mathbf{\Delta})^{-1}\widehat{\bm{\Phi}}_{uu}\|_{F}\leq\frac{\|\widehat{\bm{\Phi}}_{uu}\|_{F}}{1-\epsilon_{2}\|\widehat{\bm{\Phi}}_{uy}\|_{2}}\,.$

Next, we have

$\displaystyle\|\widehat{\bm{\Phi}}_{yy}(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}(\widehat{\mathbf{G}}+\mathbf{\Delta})\|_{F}$
$\displaystyle\leq\|\widehat{\bm{\Phi}}_{yy}\widehat{\mathbf{G}}\|_{F}\text{+}\|\widehat{\bm{\Phi}}_{yy}\mathbf{\Delta}\|_{F}\text{+}\left\|\widehat{\bm{\Phi}}_{yy}\left(\sum_{k=1}^{\infty}(\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy})^{k}\right)(\widehat{\mathbf{G}}+\mathbf{\Delta})\right\|_{F}$
$\displaystyle=\|\widehat{\bm{\Phi}}_{yu}\|_{F}+\epsilon_{2}\|\widehat{\bm{\Phi}}_{yy}\|_{F}+\|\widehat{\bm{\Phi}}_{yy}\|_{F}\frac{\epsilon_{2}\|\widehat{\bm{\Phi}}_{uy}\|_{2}(\|\widehat{\mathbf{G}}\|_{2}+\epsilon_{2})}{1-\epsilon_{2}\|\widehat{\bm{\Phi}}_{uy}\|_{2}}$
$\displaystyle\leq\frac{\|\widehat{\bm{\Phi}}_{yu}\|_{F}+\epsilon_{2}\|\widehat{\bm{\Phi}}_{yy}\|_{F}(2+\|\widehat{\bm{\Phi}}_{uy}\|_{2}\|\widehat{\mathbf{G}}\|_{2})}{1-\epsilon_{2}\|\widehat{\bm{\Phi}}_{uy}\|_{2}}\,,$

and therefore, by developing the squares and using that
$\left\lVert\widehat{\bm{\Phi}}_{yy}\widehat{\mathbf{G}}\right\rVert_{F}\leq\|\widehat{\bm{\Phi}}_{yy}\|_{F}\|\widehat{\mathbf{G}}\|_{2}$
we obtain

$\displaystyle\quad\|\widehat{\bm{\Phi}}_{yy}(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}(\widehat{\mathbf{G}}+\mathbf{\Delta})\|_{F}^{2}$
$\displaystyle\leq\frac{\left(\|\widehat{\bm{\Phi}}_{yu}\|_{F}^{2}\|_{2}+\|\widehat{\bm{\Phi}}_{yy}\|_{F}^{2}h(\epsilon_{2},\gamma,\widehat{\mathbf{G}})\right)}{(1-\epsilon_{2}\|\widehat{\bm{\Phi}}_{uy}\|_{2})^{2}}\,.$

Proceeding analogously, one can also prove that

$\displaystyle\|\widehat{\bm{\Phi}}_{yy}(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}(\widehat{\bm{y}}_{0}+\bm{\delta}_{0})\|_{F}^{2}$
$\displaystyle\leq\frac{1}{(1-\epsilon_{2}\|\widehat{\bm{\Phi}}_{uy}\|_{2})^{2}}\left(\|\widehat{\bm{\Phi}}_{yy}\widehat{\mathbf{y}}_{0}\|_{F}^{2}+\|\widehat{\bm{\Phi}}_{yy}\|_{F}^{2}h(\epsilon_{2},\gamma,\widehat{\mathbf{y}}_{0})\right)\,,$
$\displaystyle\|\widehat{\bm{\Phi}}_{uy}(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}(\widehat{\bm{y}}_{0}+\bm{\delta}_{0})\|_{F}^{2}\leq$
$\displaystyle\frac{1}{(1-\epsilon_{2}\|\widehat{\bm{\Phi}}_{uy}\|_{2})^{2}}\left(\|\widehat{\bm{\Phi}}_{uy}\widehat{\mathbf{y}}_{0}\|_{F}^{2}+\|\widehat{\bm{\Phi}}_{uy}\|_{F}^{2}h(\epsilon_{2},\gamma,\widehat{\mathbf{y}}_{0})\right)\,.$

Therefore, combining the above inequalities we finally conclude that

$\displaystyle J(\mathbf{G},\mathbf{K})$
$\displaystyle\leq\frac{\Big{(}\left\lVert\begin{bmatrix}\widehat{\bm{\Phi}}_{yy}&\widehat{\bm{\Phi}}_{yu}&\widehat{\bm{\Phi}}_{yy}\widehat{\mathbf{y}}_{0}\\\
\widehat{\bm{\Phi}}_{uy}&\widehat{\bm{\Phi}}_{uu}&\widehat{\bm{\Phi}}_{uy}\widehat{\mathbf{y}}_{0}\end{bmatrix}\right\rVert_{F}^{2}+\|\widehat{\bm{\Phi}}_{yy}\|_{F}^{2}(h(\epsilon_{2},\gamma,\widehat{\mathbf{G}})+}{1-\epsilon_{2}\|\widehat{\bm{\Phi}}_{uy}\|_{2}}$
$\displaystyle\qquad\frac{+h(\epsilon_{2},\gamma,\widehat{\mathbf{y}}_{0}))+\|\widehat{\bm{\Phi}}_{uy}\|_{F}^{2}h(\epsilon_{2},\gamma,\widehat{\mathbf{y}}_{0})\Big{)}^{\frac{1}{2}}}{1-\epsilon_{2}\|\widehat{\bm{\Phi}}_{uy}\|_{2}}\,.$

###  .4 Proof of Lemma 3.8

By using the fact that for $x\in\mathbb{R}^{n}$ and $y\in\mathbb{R}^{m}$ we
have that
$\left\lVert\begin{bmatrix}x^{\mathsf{T}}&y^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}\right\rVert_{1}=\left\lVert
x\right\rVert_{1}+\left\lVert y\right\rVert_{1}$, the left-hand-sides of
(31)-(32) are each made of three addends. The proof hinges on upper-bounding
each one of them. We report the full derivations for the most informative of
them. We have

$\displaystyle v_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}\right\rVert_{1}$
$\displaystyle\leq v_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}+\max_{\left\lVert\mathbf{v}\right\rVert_{\infty}\leq
v_{\infty}}\lvert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\bm{\Delta}\widehat{\bm{\Phi}}_{uy}(I-\bm{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}\mathbf{v}\rvert$
$\displaystyle\leq v_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}+v_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}\left\lVert\bm{\Delta}\widehat{\bm{\Phi}}_{uy}(I-\bm{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}\right\rVert_{\infty}$
$\displaystyle\leq v_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}+v_{\infty}\epsilon_{\infty}\frac{\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}\left\lVert\widehat{\bm{\Phi}}_{uy}\right\rVert_{\infty}}{1-\epsilon_{\infty}\left\lVert\widehat{\bm{\Phi}}_{uy}\right\rVert_{\infty}}$
$\displaystyle\leq v_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}\left(1+\frac{\epsilon_{\infty}\tau}{1-\epsilon_{\infty}\tau}\right)$
$\displaystyle=\frac{v_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}}{1-\epsilon_{\infty}\tau}=f_{1,j}(\widehat{\bm{\Phi}})\,.$

Next, recalling
$\widehat{\bm{\Phi}}_{yu}=\widehat{\bm{\Phi}}_{yy}\widehat{\mathbf{G}}$, we
have

$\displaystyle w_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}(I-\mathbf{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}(\widehat{\mathbf{G}}+\mathbf{\Delta})\right\rVert_{1}$
$\displaystyle=\max_{\left\lVert\mathbf{w}\right\rVert_{\infty}\leq
w_{\infty}}F_{y,j}\left(\widehat{\bm{\Phi}}_{yy}(\widehat{\mathbf{G}}+\bm{\Delta})+\widehat{\bm{\Phi}}_{yy}\bm{\Delta}\widehat{\bm{\Phi}}_{uy}(I-\bm{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}(\widehat{\mathbf{G}}+\bm{\Delta})\right)\mathbf{w}$
$\displaystyle\leq w_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yu}\right\rVert_{1}+\max_{\left\lVert\mathbf{w}\right\rVert_{\infty}\leq
w_{\infty}}|F_{y,j}\widehat{\bm{\Phi}}_{yy}\bm{\Delta}\mathbf{w}|+$
$\displaystyle~{}~{}~{}+\max_{\left\lVert\mathbf{w}\right\rVert_{\infty}\leq
w_{\infty}}\lvert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\bm{\Delta}\widehat{\bm{\Phi}}_{uy}(I-\bm{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}(\widehat{\mathbf{G}}+\bm{\Delta})\mathbf{w}\rvert$
$\displaystyle\leq w_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yu}\right\rVert_{1}+w_{\infty}\epsilon_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}+$
$\displaystyle~{}~{}~{}+w_{\infty}\epsilon_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}\left\lVert\widehat{\bm{\Phi}}_{uy}(I-\bm{\Delta}\widehat{\bm{\Phi}}_{uy})^{-1}(\widehat{\mathbf{G}}+\bm{\Delta})\right\rVert_{\infty}$
$\displaystyle\leq w_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yu}\right\rVert_{1}\text{+}w_{\infty}\epsilon_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}\left(1\text{+}\tau\frac{\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}+\epsilon_{\infty}}{1-\epsilon_{\infty}\tau}\right)$
$\displaystyle=w_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yu}\right\rVert_{1}+w_{\infty}\epsilon_{\infty}\left\lVert
F_{y,j}\widehat{\bm{\Phi}}_{yy}\right\rVert_{1}\left(\frac{1+\tau\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}}{1-\epsilon_{\infty}\tau}\right)$
$\displaystyle=f_{2,j}(\widehat{\bm{\Phi}})\,.$

Lastly, remembering that
$\widehat{\bm{\Phi}}_{uu}=I+\widehat{\bm{\Phi}}_{uy}\widehat{\mathbf{G}}$ and
noticing that

$(I-\widehat{\bm{\Phi}}_{uy}\mathbf{\Delta})^{-1}\widehat{\bm{\Phi}}_{uu}=\widehat{\bm{\Phi}}_{uu}+\widehat{\bm{\Phi}}_{uy}\bm{\Delta}(I-\widehat{\bm{\Phi}}_{uy}\bm{\Delta})^{-1}\widehat{\bm{\Phi}}_{uu}\,,$

we have

$\displaystyle w_{\infty}\left\lVert
F_{u,j}(I-\widehat{\bm{\Phi}}_{uy}\mathbf{\Delta})^{-1}\widehat{\bm{\Phi}}_{uu}\right\rVert_{1}$
$\displaystyle\leq w_{\infty}\left\lVert
F_{u,j}\widehat{\bm{\Phi}}_{uu}\right\rVert_{1}+$
$\displaystyle+w_{\infty}\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}\epsilon_{\infty}\frac{\left\lVert
F_{u,j}\widehat{\bm{\Phi}}_{uy}\right\rVert_{1}\left\lVert\widehat{\bm{\Phi}}_{uy}\right\rVert_{\infty}}{1-\epsilon_{\infty}\left\lVert\widehat{\bm{\Phi}}_{uy}\right\rVert_{\infty}}\text{+}w_{\infty}\epsilon_{\infty}\frac{\left\lVert
F_{u,j}\widehat{\bm{\Phi}}_{uy}\right\rVert_{1}}{1-\epsilon_{\infty}\left\lVert\widehat{\bm{\Phi}}_{uy}\right\rVert_{\infty}}$
$\displaystyle\leq w_{\infty}\left\lVert
F_{u,j}\widehat{\bm{\Phi}}_{uu}\right\rVert_{1}+w_{\infty}\epsilon_{\infty}\left\lVert
F_{u,j}\widehat{\bm{\Phi}}_{uy}\right\rVert_{1}\frac{1+\tau\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}}{1-\epsilon_{\infty}\tau}$
$\displaystyle=f_{5,j}(\widehat{\bm{\Phi}})\,.$

Similar computations allows one to derive the upperbounds for the remaining
terms.

###  .5 Proof of Lemma 4.14

First, it is easy to verify that $\widetilde{\bm{\Phi}}$ satisfies the
constraints in (41); indeed, $\widetilde{\bm{\Phi}}$ comprises the closed-loop
responses when we apply $\mathbf{K}^{c}$ to the estimated plant
$\widehat{\mathbf{G}}$. Next, we have

$\displaystyle\left\lVert\widetilde{\bm{\Phi}}_{uy}\right\rVert_{2}$
$\displaystyle=\left\lVert\bm{\Phi}^{c}_{uy}(I+\bm{\Delta}\bm{\Phi}^{c}_{uy})^{-1}\right\rVert_{2}$
$\displaystyle\leq\frac{\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2}}{1-\epsilon_{2}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2}}\leq\sqrt{2}\frac{\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2}}{1-\epsilon_{2}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2}}$
$\displaystyle\leq\sqrt{2}\frac{\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}}{1-\epsilon_{2}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}}$
(52)
$\displaystyle=\sqrt{2}\frac{\eta}{\epsilon_{2}(1-\eta)}=\widetilde{\gamma}\,.$
(53)

Since $\alpha\in[\sqrt{2}\frac{\eta}{\epsilon_{2}(1-\eta)},\epsilon_{2}^{-1})$
and $\eta<\frac{1}{5}$, then $\tilde{\gamma}\leq\alpha<\epsilon_{2}^{-1}$.
Hence $\widetilde{\gamma}$ is feasible. Similarly,

$\displaystyle\left\lVert\widetilde{\bm{\Phi}}_{uy}\right\rVert_{\infty}$
$\displaystyle=\left\lVert\bm{\Phi}^{c}_{uy}(I+\bm{\Delta}\bm{\Phi}^{c}_{uy})^{-1}\right\rVert_{\infty}$
$\displaystyle\leq\frac{\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{\infty}}{1-\epsilon_{\infty}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{\infty}}\leq\frac{\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{\infty}}{1-\epsilon_{\infty}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{\infty}}$
(54)
$\displaystyle=\frac{\zeta}{\epsilon_{\infty}(1-\zeta)}=\widetilde{\tau}\,.$
(55)

Since $\zeta<\frac{1}{2}$, then $\tilde{\tau}<\epsilon_{\infty}^{-1}$ and
hence it is a feasible value for $\tau$. It remains to show that
$\widetilde{\bm{\Phi}}$ satisfies the safety constraints (36)-(37). We know
that $\bm{\Phi}^{c}$ is feasible for (42), and hence
$\phi_{1,j}(\bm{\Phi}^{c})+\phi_{2,j}(\bm{\Phi}^{c})+\phi_{3,j}(\bm{\Phi}^{c})\leq\mathbf{b}_{y,j}$
and
$\phi_{4,j}(\bm{\Phi}^{c})+\phi_{5,j}(\bm{\Phi}^{c})+\phi_{6,j}(\bm{\Phi}^{c})\leq\mathbf{b}_{u,j}$.
We conclude the proof by showing that
$f_{i,j}(\widetilde{\bm{\Phi}})\leq\phi_{i,j}(\bm{\Phi}^{c})$ for every
$i=1,\dots,6$. We report the full derivations for the most informative terms.

$\displaystyle f_{1,j}(\widetilde{\bm{\Phi}})$
$\displaystyle=\frac{v_{\infty}\left\lVert
F_{y,j}\left(\bm{\Phi}^{c}_{yy}-\bm{\Phi}^{c}_{yy}\bm{\Delta}\bm{\Phi}^{c}_{uy}\left(I+\bm{\Delta}\bm{\Phi}^{c}_{uy}\right)^{-1}\right)\right\rVert_{1}}{1-\epsilon_{\infty}\widetilde{\tau}}$
$\displaystyle\leq\frac{v_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}\right\rVert_{1}+\frac{v_{\infty}\epsilon_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}\right\rVert_{1}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{\infty}}{1-\epsilon_{\infty}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{\infty}}}{1-\epsilon_{\infty}\widetilde{\tau}}$
$\displaystyle\leq\frac{v_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}\right\rVert_{1}+\frac{v_{\infty}\epsilon_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}\right\rVert_{1}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{\infty}}{1-\epsilon_{\infty}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{\infty}}}{1-\epsilon_{\infty}\widetilde{\tau}}$
$\displaystyle\leq\frac{v_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}\right\rVert_{1}+\frac{\zeta}{1-\zeta}v_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}\right\rVert_{1}}{1-\frac{\zeta}{1-\zeta}}$
$\displaystyle\leq\frac{v_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}\right\rVert_{1}}{1-2\zeta}=\phi_{1,j}(\bm{\Phi}^{c})\,.$

Similarly, it is easy to show that
$f_{4,j}(\widetilde{\bm{\Phi}})\leq\phi_{4,j}(\bm{\Phi}^{c})$. Next, by
remembering that

$\displaystyle\widetilde{\bm{\Phi}}_{yu}=\bm{\Phi}^{c}_{yu}-\bm{\Phi}^{c}_{yy}\bm{\Delta}-\bm{\Phi}^{c}_{yy}\bm{\Delta}\bm{\Phi}^{c}_{uy}(I+\bm{\Delta}\bm{\Phi}^{c}_{uy})^{-1}(\mathbf{G}-\bm{\Delta})\,,$
$\displaystyle\widetilde{\bm{\Phi}}_{yy}=\bm{\Phi}^{c}_{yy}-\bm{\Phi}^{c}_{yy}\bm{\Delta}\bm{\Phi}^{c}_{uy}\left(I+\bm{\Delta}\bm{\Phi}^{c}_{uy}\right)^{-1}\,,$

we have

$\displaystyle f_{2,j}(\widetilde{\bm{\Phi}})$ $\displaystyle\leq
w_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}(I+\bm{\Delta}\bm{\Phi}^{c}_{uy})^{-1}(\mathbf{G}-\bm{\Delta})\right\rVert_{1}+$
$\displaystyle~{}~{}+w_{\infty}\epsilon_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}(I+\bm{\Delta}\bm{\Phi}^{c}_{uy})^{-1}\right\rVert_{1}\left(\frac{1+\widetilde{\tau}\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}}{1-\epsilon_{\infty}\widetilde{\tau}}\right)$
$\displaystyle\leq w_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yu}\right\rVert_{1}\text{+}w_{\infty}\epsilon_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}\right\rVert_{1}\left(1\text{+}\frac{\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{\infty}\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}}{1-\epsilon_{\infty}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{\infty}}\right)\text{+}$
$\displaystyle~{}~{}+w_{\infty}\epsilon_{\infty}\frac{\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}\right\rVert_{1}\left(\frac{1+\widetilde{\tau}\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}}{1-\epsilon_{\infty}\widetilde{\tau}}\right)}{1-\epsilon_{\infty}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{\infty}}$
$\displaystyle\leq w_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yu}\right\rVert_{1}+$
$\displaystyle~{}~{}+w_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}\right\rVert_{1}\left(\epsilon_{\infty}+\frac{\zeta\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}}{1-\zeta}+\frac{\epsilon_{\infty}+\frac{\zeta\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}}{1-\zeta}}{\left(1-\frac{\zeta}{1-\zeta}\right)(1-\zeta)}\right)$
$\displaystyle=w_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yu}\right\rVert_{1}\text{+}2w_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}\right\rVert_{1}\frac{(1-\zeta)\left(\epsilon_{\infty}+\zeta\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}\right)}{1-2\zeta}$
$\displaystyle\leq w_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yu}\right\rVert_{1}+2w_{\infty}\left\lVert
F_{y,j}\bm{\Phi}^{c}_{yy}\right\rVert_{1}\frac{\left(\epsilon_{\infty}+\zeta\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}\right)}{1-2\zeta}$
$\displaystyle\leq\phi_{2,j}(\bm{\Phi}^{c})\,.$

Similarly, $f_{3,j}(\widetilde{\bm{\Phi}})\leq\phi_{3,j}(\bm{\Phi}^{c})$ and
$f_{6,j}(\widetilde{\bm{\Phi}})\leq\phi_{6,j}(\bm{\Phi}^{c})$. By only
noticing that $\left\lVert\bm{\Phi}^{c}_{uu}\right\rVert_{\infty}\leq
1+\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{\infty}\left(\left\lVert\widehat{\mathbf{G}}\right\rVert_{\infty}+\epsilon_{\infty}\right)$
and that $(1+\zeta)(1-2\zeta)\leq 1-\zeta$ for every $\zeta>0$, analogous
computations lead to
$f_{5,j}(\widetilde{\bm{\Phi}})\leq\phi_{5,j}(\bm{\Phi}^{c})$.

###  .6 Proof of Theorem 4.15

By denoting as $\widehat{\bm{\Phi}}^{\star}$ the closed-loop responses
obtained by applying $\widehat{\mathbf{K}}^{\star}$ to $\widehat{\mathbf{G}}$,
we have by Lemma 3.7 and by $\gamma\leq\alpha$

$\displaystyle
J(\mathbf{G},\widehat{\mathbf{K}}^{\star})\leq\frac{1}{1-\epsilon_{2}\gamma^{\star}}\times$
$\displaystyle\times\left\lVert\begin{bmatrix}\sqrt{1+h(\epsilon_{2},\alpha,\widehat{\mathbf{G}})+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})}\widehat{\bm{\Phi}}^{\star}_{yy}&\widehat{\bm{\Phi}}^{\star}_{yu}&\widehat{\bm{\Phi}}^{\star}_{yy}\widehat{\mathbf{y}}_{0}\\\
\sqrt{1+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})}\widehat{\bm{\Phi}}^{\star}_{uy}&\widehat{\bm{\Phi}}^{\star}_{uu}&\widehat{\bm{\Phi}}^{\star}_{uy}\widehat{\mathbf{y}}_{0}\end{bmatrix}\right\rVert_{F},$

where $\gamma^{\star}$ is optimal for (41). By Lemma 4.14, under the
assumptions on $\eta,\zeta,\alpha$ we have that
$(\widetilde{\gamma},\widetilde{\tau},\widetilde{\bm{\Phi}})$ belongs to the
feasible set of (41). Hence, by suboptimality of any feasible solution, one
has

$\displaystyle
J(\mathbf{G},\widehat{\mathbf{K}}^{\star})\leq\frac{1}{1-\epsilon_{2}\widetilde{\gamma}}\times$
$\displaystyle\times\left\lVert\begin{bmatrix}\sqrt{1+h(\epsilon_{2},\alpha,\widehat{\mathbf{G}})+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})}\widetilde{\bm{\Phi}}_{yy}&\widetilde{\bm{\Phi}}_{yu}&\widetilde{\bm{\Phi}}_{yy}\widehat{\mathbf{y}}_{0}\\\
\sqrt{1+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})}\widetilde{\bm{\Phi}}_{uy}&\widetilde{\bm{\Phi}}_{uu}&\widetilde{\bm{\Phi}}_{uy}\widehat{\mathbf{y}}_{0}\end{bmatrix}\right\rVert_{F}.$

Using the definition of $\widetilde{\bm{\Phi}}$ from Lemma 4.14, we now relate
the term

$\widetilde{C}=\left\lVert\begin{bmatrix}\sqrt{1+h(\epsilon_{2},\alpha,\widehat{\mathbf{G}})+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})}\widetilde{\bm{\Phi}}_{yy}&\widetilde{\bm{\Phi}}_{yu}&\widetilde{\bm{\Phi}}_{yy}\widehat{\mathbf{y}}_{0}\\\
\sqrt{1+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})}\widetilde{\bm{\Phi}}_{uy}&\widetilde{\bm{\Phi}}_{uu}&\widetilde{\bm{\Phi}}_{uy}\widehat{\mathbf{y}}_{0}\end{bmatrix}\right\rVert_{F},$

to the optimal cost of problem (42). Define the quantities

$\displaystyle M^{c}$
$\displaystyle=h(\epsilon_{2},\alpha,\widehat{\mathbf{G}})+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})+$
$\displaystyle~{}~{}+h(\epsilon_{2},\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2},\mathbf{G})+h(\epsilon_{2},\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2},\mathbf{y}_{0})\,,$
$\displaystyle M^{\star}$
$\displaystyle=h(\epsilon_{2},\alpha,\widehat{\mathbf{G}})+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})+$
$\displaystyle~{}~{}+h(\epsilon_{2},\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2},\mathbf{G})+h(\epsilon_{2},\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2},\mathbf{y}_{0})\,,$
$\displaystyle V^{c}$
$\displaystyle=h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})+h(\epsilon_{2},\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2},\mathbf{y}_{0})\,,$
$\displaystyle V^{\star}$
$\displaystyle=h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})+h(\epsilon_{2},\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2},\mathbf{y}_{0})\,.$

By using analogous expressions as those derived in Lemma 4.14,

$\displaystyle\widetilde{C}=\Big{(}\left\lVert\begin{bmatrix}\widetilde{\bm{\Phi}}_{yy}&\widetilde{\bm{\Phi}}_{yu}&\widetilde{\bm{\Phi}}_{yy}\widehat{\mathbf{y}}_{0}\\\
\widetilde{\bm{\Phi}}_{uy}&\widetilde{\bm{\Phi}}_{uu}&\widetilde{\bm{\Phi}}_{uy}\widehat{\mathbf{y}}_{0}\end{bmatrix}\right\rVert_{F}^{2}+$
$\displaystyle+\left(h(\epsilon_{2},\alpha,\widehat{\mathbf{G}})+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})\right)\left\lVert\widetilde{\bm{\Phi}}_{yy}\right\rVert_{F}^{2}+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})\left\lVert\widetilde{\bm{\Phi}}_{uy}\right\rVert_{F}^{2}\Big{)}^{\frac{1}{2}}$
$\displaystyle\leq\frac{\sqrt{J(\mathbf{G},\mathbf{K}^{c})^{2}+M^{c}\left\lVert\bm{\Phi}^{c}_{yy}\right\rVert_{F}^{2}+V^{c}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{F}^{2}}}{1-\epsilon_{2}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2}}\,.$

Thus, we have established the chain of inequalities

$\displaystyle J(\mathbf{G},\widehat{\mathbf{K}}^{\star})^{2}$
$\displaystyle\leq\frac{1}{(1-\epsilon_{2}\widetilde{\gamma})^{2}}\widetilde{C}^{2}$
$\displaystyle\leq\frac{1}{(1-\epsilon_{2}\widetilde{\gamma})^{2}}\frac{1}{(1-\epsilon_{2}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2})^{2}}\Big{(}J(\mathbf{G},\mathbf{K}^{c})^{2}+$
$\displaystyle~{}~{}+M^{c}\left\lVert\bm{\Phi}^{c}_{yy}\right\rVert_{F}^{2}+V^{c}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{F}^{2}\Big{)}\,.$

Next, notice that, by definition, we have

$J(\mathbf{G},\mathbf{K}^{c})^{2}=(S(\epsilon_{\infty})+1)J(\mathbf{G},\mathbf{K}^{\star})^{2}\,.$

Recalling that $\eta<\frac{1}{5}$,
$\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2}\leq\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}$
and
$\left\lVert\bm{\Phi}^{c}_{yy}\right\rVert_{2}\leq\left\lVert\bm{\Phi}^{\star}_{yy}\right\rVert_{2}$,
and noticing that if $M,V>0$, then

$Ma^{2}+Vb^{2}\leq(M+V)(a^{2}+b^{2})\,,$

we can establish the following inequalities.

$\displaystyle\frac{J(\mathbf{G},\widehat{\mathbf{K}}^{\star})^{2}-J(\mathbf{G},\mathbf{K}^{\star})^{2}}{J(\mathbf{G},\mathbf{K}^{\star})^{2}}\leq$
$\displaystyle\leq\left(\frac{1}{(1-\epsilon_{2}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{2})^{2}(1-\epsilon_{2}\widetilde{\gamma})^{2}}\right)+$
$\displaystyle~{}~{}~{}+\left(S(\epsilon_{\infty})+1+\frac{M^{c}\left\lVert\bm{\Phi}^{c}_{yy}\right\rVert_{F}^{2}+V^{c}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{F}^{2}}{J(\mathbf{G},\mathbf{K}^{\star})^{2}}\right)-1$
$\displaystyle\leq\left(\frac{1}{(1-\eta)^{2}(1-\sqrt{2}\frac{\eta}{1-\eta})^{2}}-1+\frac{S(\epsilon_{\infty})}{(1-\eta)^{2}(1-\sqrt{2}\frac{\eta}{1-\eta})^{2}}\right)+$
$\displaystyle~{}~{}~{}+\frac{M^{c}\left\lVert\bm{\Phi}^{c}_{yy}\right\rVert_{F}^{2}+V^{c}\left\lVert\bm{\Phi}^{c}_{uy}\right\rVert_{F}^{2}}{(1-\eta)^{2}(1-\sqrt{2}\frac{\eta}{1-\eta})^{2}J(\mathbf{G},\mathbf{K}^{\star})^{2}}$
$\displaystyle\leq\eta\left(\frac{2(1+\sqrt{2})-(1+\sqrt{2})^{2}\eta}{(1-(1+\sqrt{2})\eta)^{2}}\right)+\frac{S(\epsilon_{\infty})}{(1-(1+\sqrt{2})\eta)^{2}}+$
$\displaystyle~{}~{}~{}+\frac{(M^{c}+V^{c})J(\mathbf{G},\mathbf{K}^{c})^{2}}{(1-(1+\sqrt{2})\eta)^{2}J(\mathbf{G},\mathbf{K}^{\star})^{2}}$
$\displaystyle\leq
20\eta+4S(\epsilon_{\infty})+4(M^{c}+V^{c})(S(\epsilon_{\infty})+1)\,,$
$\displaystyle\leq
20\eta+4(M^{c}+V^{c})+4S(\epsilon_{\infty})(1+M^{c}+V^{c})\,.$

Last, we prove that
$20\eta+4(M^{c}+V^{c})=\mathcal{O}\left(\epsilon_{2}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}(\left\lVert\mathbf{G}\right\rVert^{2}_{2}+\left\lVert\mathbf{y}_{0}\right\rVert^{2}_{2})\right)$.
First, notice that $M^{c}+V^{c}\leq M^{\star}+V^{\star}\,.$ By considering the
expressions of $M^{\star}$ and $V^{\star}$, using $\alpha\leq
5\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}$, $\eta<\frac{1}{5}$,
$\left\lVert\widehat{\mathbf{G}}\right\rVert_{2}\leq\left\lVert\mathbf{G}\right\rVert_{2}+\epsilon_{2}$
and
$\left\lVert\widehat{\mathbf{y}}_{0}\right\rVert_{2}\leq\left\lVert\mathbf{y}_{0}\right\rVert_{2}+\epsilon_{2}$,
we deduce that

$\displaystyle
M^{\star}=h(\epsilon_{2},\alpha\,\widehat{\mathbf{G}})+h(\epsilon_{2},\alpha,\widehat{\mathbf{y}}_{0})+h(\epsilon_{2},\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2},\mathbf{G})+$
$\displaystyle~{}~{}~{}~{}+h(\epsilon_{2},\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2},\mathbf{y}_{0})$
$\displaystyle\leq
2\Big{[}\epsilon_{2}^{2}(2+5\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}\|\mathbf{G}\|_{2})^{2}+2\epsilon_{2}\left\lVert\mathbf{G}\right\rVert_{2}(2+$
$\displaystyle~{}~{}~{}+5\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}\left\lVert\mathbf{G}\right\rVert_{2})+\epsilon_{2}^{2}(2+5\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}\|\mathbf{y}_{0}\|_{2})^{2}+$
$\displaystyle~{}~{}~{}+2\epsilon_{2}\left\lVert\mathbf{y}_{0}\right\rVert_{2}(2+5\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}\left\lVert\mathbf{y}_{0}\right\rVert_{2})\Big{]}+$
$\displaystyle~{}~{}~{}+\mathcal{O}(\epsilon_{2}^{2}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}(\left\lVert\mathbf{G}\right\rVert^{2}_{2}+\left\lVert\mathbf{y}_{0}\right\rVert^{2}_{2}))$
$\displaystyle=\mathcal{O}\left(\epsilon_{2}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}(\left\lVert\mathbf{G}\right\rVert^{2}_{2}+\left\lVert\mathbf{y}_{0}\right\rVert^{2}_{2})\right)\,,$

and, similarly,
$V^{\star}=\mathcal{O}\left(\epsilon_{2}\left\lVert\bm{\Phi}^{\star}_{uy}\right\rVert_{2}\left\lVert\mathbf{y}_{0}\right\rVert^{2}_{2}\right)$.