We discuss connections between sequential system identification and control for linear time-invariant systems,
often termed
indirect data-driven control, as well as a  contemporary direct data-driven control approach seeking an optimal decision compatible with recorded data assembled in a Hankel matrix and robustified through suitable regularizations. We formulate these two problems in the language of behavioral systems theory and parametric mathematical programs, and we bridge them through a multi-criteria formulation trading off system identification and control objectives. We illustrate our results with two methods from subspace identification and control: namely, subspace predictive control and low-rank approximation which constrain trajectories to be consistent with a non-parametric predictor derived from (respectively, the column span of) a data Hankel matrix. In both cases we conclude that direct and regularized data-driven control can be derived as convex relaxation of the indirect approach, and the regularizations account for an implicit identification step. Our analysis further reveals a novel regularizer and a plausible   hypothesis explaining the remarkable empirical performance of direct methods on nonlinear systems.

§ INTRODUCTION

The vast realm of data-driven control methods can be classified into indirect data-driven control approaches consisting of sequential system identification and model-based control as well as direct data-driven control approaches seeking an optimal decision compatible with recorded data. Both approaches have a rich history, and they have received renewed interest cross-fertilized by novel methods and widespread interest in machine learning. Representative recent surveys are [1, 2, 3, 4, 5, 6, 7, 8].

The pros and cons of both paradigms have often been elaborated on. The indirect approach is modular with well understood subtasks, though modeling and identification are cumbersome, their results are often not useful for control (due to, e.g., incompatible uncertainty quantifications), and practitioners  often prefer end-to-end methods.  Direct approaches promise to resolve these problems by learning control policies directly from data. However, they are often analytically and computationally less tractable and rarely apply to real-time and safety-critical control systems. Selected direct methods that proved themselves in theory and practice are iterative feedback tuning and virtual reference feedback tuning [9, 10, 11].

Quite a few approaches have bridged the direct and indirect data-driven control paradigms. Of relevance to this article, we note the literature on identification for control [7, 12, 13, 14] and control-oriented regularized identification [15], which propose that the control objective should bias the identification task. Likewise, dual control dating to [16] addresses the  exploration vs. exploitation trade-offs in simultaneous identification and optimal control; see [17, 18, 19] for  recent contributions. Furthermore, [20] formulates data-driven model reference control as an identification problem, where various degrees of prior information can be incorporated so that the method can range between the direct and the indirect approach.

We take a similar perspective here:  the sequential identification and control tasks can be abstracted as nested bi-level optimization problem: find the best control subject to a model, where the model is the best fit to a data set within some hypothesis class. This approach is modular and both steps admit tractable formulations, but generally it is also suboptimal: there is no separation principle – aside from special cases, see <cit.> – for these two nested optimization problems. An end-to-end direct algorithmic approach may thus outperform indirect methods if a tractable formulation were available.
For the latter we resort to a paradigm square in between behavioral system theory and subspace system identification methods.

Behavioral system theory [21, 22, 23] takes an abstract view on dynamical systems as sets of trajectories, and it does not require parametric representations which makes it appealing from a data-centric perspective. For example, linear time-invariant (LTI) systems are characterized as shift-invariant subspaces within an ambient  space of time series. The role of identification is to find such a low-dimensional feature from data. Subspace methods take a similar (albeit more algorithmic) viewpoint [24, 25, 26] and extract parametric models from the range and null spaces of a low-rank data Hankel matrix.

Both lines of work come together in a result known as the Fundamental Lemma [27]; see also [28, 29, 6] for recent extensions. It states that, under some assumptions, the set of all finite-length trajectories (the restricted behavior) of an LTI system equals the range space of a data Hankel matrix. This result serves as the theoretic underpinning for work in subspace identification [29, 30, 31]
and data-driven control, in particular subspace predictive control based on non-parametric models [32, 33, 34],  explicit feedback policies parametrized by data matrices [35, 36, 37],
and data-enabled predictive control (DeePC) seeking compatibility of predicted trajectories with the range space of a data Hankel matrix. The latter methods have first been established for deterministic LTI systems in [38, 39] and have recently been extended by suitably regularizing the optimal control problems. Closed-loop stability was certified in [40]. The regularizations were first mere heuristics [41] but have later been constructively derived by robust control and optimization [42, 43, 44, 45, 46].  These approaches, albeit recent, have proved themselves in practical nonlinear problems in multiple domains [46, 47, 48, 49, 45]. We also note the recent maximum-likelihood perspective [50].

We refer to [6] surveying results surrounding the fundamental lemma.

In this paper, we explore the following questions: how does data-enabled predictive control relate to a prior system identification? What are principled regularizations? And why does it work so well in the nonlinear case? We start our investigations from indirect data-driven control formulated as a bi-level optimization problem  in the general output feedback setting.

As a vehicle to transition between indirect and direct approaches, we consider a multi-criteria problem trading off identification and control objectives reminiscent of similar approaches [7, 12, 13, 14, 15, 16, 17, 18, 19, 20] blending the two. We formally show that one tail of its Pareto front corresponds to the bi-level problem, and a convex relaxation results in the regularized data-enabled predictive control formulations used in [40, 41, 42, 43, 44, 45, 47, 48, 49, 46].

Most of our results are formulated in the abstract language of behavioral systems theory and parametric mathematical programs, but we also specialize our treatment to two concrete methods: subspace predictive control (SPC) [32, 33, 34] and low-rank approximation [39]. In both cases we conclude that the direct regularized data-driven control can be derived as a convex relaxation of the indirect approach, where $(i)$ LTI complexity specifications (selecting the model class) are dropped, and $(ii)$ the projection of the data on the set of LTI systems is replaced by regularizations accounting for implicit identification.
In particular, starting from indirect data-driven control based on low-rank approximation of a Hankel matrix, we arrive at a DeePC formulation with an $\ell_{1}$-regularizer (Theorem <ref>). When formulating indirect data-driven control via the SPC framework, our analysis reveals a novel regularizer for DeePC promoting a least-square data fit  by projecting on the null space of the Hankel matrix (Theorem <ref>).

We illustrate our results with numerical studies illustrating the role of regularization, superiority of the new regularizer, and comparisons. Informed by our analysis, we hypothesize and numerically confirm that the indirect approach is superior in case of “variance” error, e.g., for LTI stochastic systems, and the direct approach wins in terms of “bias” error, e.g., for nonlinear systems supporting the empirical observations in [46, 47, 48, 49, 45].
Similar bias-variance trade-offs can also be found in the recent pre-print [51] discussing sub-optimality of direct and indirect methods as  function of the data size. These findings also resonate with those of data-driven model reference control [20] concluding that the direct approach is superior in reducing the bias whereas the indirect one gives better variance – especially if an erroneous model class is selected.

The remainder of this paper is organized as follows: Section <ref> reviews representations of LTI systems. Section <ref> formulates the direct and indirect data-driven control problems, and Section <ref> bridges them. Section <ref> contains our numerical studies. Finally, Section <ref> concludes the paper. Readers familiar with the behavioral approach may skip Section <ref>.

§ LTI SYSTEMS AND THEIR REPRESENTATIONS

We adopt a behavioral perspective which allows for system theory independent of parametric representations.
We aim at a concise exposition and refer to [21, 22, 23, 6] for details.

§.§ Behavioral Perspective on Discrete-Time LTI systems

Consider the discrete time axis $\mathbb Z$, the signal space $\real^{q}$, and the associated space of trajectories  $\real^{q\mathbb Z}$ consisting of all $q$-variate sequences $(\dots,w(-1),w(0), w(1),\dots)$ with $w(i) \in \real^{q}$.
Consider a permutation matrix $P$ partitioning each
= P
\left[\begin{smallmatrix}
u(i) \\ y(i)
\end{smallmatrix}\right]

where $u(i) \in \real^{m}$ and  $y(i) \in \real^{q-m}$ are free and dependent variables that will later serve as inputs and outputs. The behavior $\bv$ is defined as a subset of the space of trajectories, $\bv \subset \real^{q\mathbb Z}$, and a system as the triple $(\mathbb Z,\real^{q},\bv)$.

In what follows, we denote a system merely by its behavior $\bv$, keeping the  signal space  $ \real^{q\mathbb Z}$ fixed throughout. A system is linear if $\bv$ is a subspace of $ \real^{q\mathbb Z}$. Let $\sigma$ denote the shift operator with action $\sigma w({t}) = w({t+1})$. A system is time-invariant if $\bv$ is shift-invariant: $\sigma \bv = \bv$.

Finally, $\bv_{L}$ is the restriction of $\bv$ to $ \real^{qL}$, i.e., to trajectories of length $L\in \mathbb Z_{>0}$.

§.§ Kernel Representations and Parametric Models

Rather than a mere set-theoretic descriptions, one typically works with explicit parametric representations (colloquially termed models) of LTI systems.
For instance, a kernel representation with lag $\ell$ specifies an LTI behavior  as

\begin{equation*}
\bv = \text{kernel}(R(\sigma)) = \bigl\{ w \in \real^{q\mathbb Z}\,:\; R(\sigma) w = 0 \bigr\}\,,
\end{equation*}

where $R(\sigma) = R_{0}+R_{1} \sigma + \dots + R_{\ell} \sigma^{\ell}$ is a polynomial matrix of degree $\ell$, and the matrices $R_{0}, R_{1}, \dots, R_{\ell}$ take values in $\real^{(q-m) \times q}$.
Alternatively, one can unfold the kernel representation  by revealing a latent variable: the state $x(t) \in \real^{n}$. The input/state/output (or state-space) representation is

\begin{align*}
\bv = & \bigl\{  w = P \left[\begin{smallmatrix} u \\ y\end{smallmatrix}\right]  \in \real^{q\mathbb Z}\,:\; \exists x \in \real^{n\mathbb Z} \mbox{ such that }
\\
&  \quad \sigma x = Ax + Bu\,,\, y = Cx+Du \bigr\}\,,
\end{align*}

where $A \in \real^{n \times n}$,  $B \in \real^{n \times m}$, $C \in \real^{q-m \times n}$, and  $D \in \real^{q-m \times m}$.
We assume that the lag $\ell$ (resp., the state dimension $n$) is minimal, i.e., there is no other kernel (resp., state-space) representation with smaller lag (resp., state dimension).

The dimension $n$ of a minimal state-space representation manifests itself in a minimal kernel representation as $n=\sum_{i=1}^{q-m} \ell_{i}$, where $\ell_{i}$ is the degree
of the $i$th row of $R(\sigma)$.

§.§ Representation-Free Estimation and Behavior Dimension

Given a state-space representation with $m$ inputs, order $n$, and lag $\ell$, the extended observability and convolution matrices

\begin{equation*}
\obs_{L}
\left[\begin{smallmatrix}
C \\ CA \\ \vdots \\ C A^{L-1}
\end{smallmatrix}\right]
\quad\mbox{and}\quad
\conv_{L}
\left[\begin{smallmatrix}
D & 0 & \cdots & & 0
\\
CB & D & 0 & \cdots & 0
\\
CAB & CB & D & \ddots & \vdots
\\
\vdots & \ddots & \ddots &\ddots & 0
\\
CA^{L-2}B & \cdots & CAB & CB & D
\end{smallmatrix}\right]
\end{equation*}

parametrize all length-$L$ trajectories in $\bv_{L}$ as

\begin{equation}
%w = \Pi \left[\begin{smallmatrix} u \\ y\end{smallmatrix}\right]
\begin{bmatrix}
u \\ y
\end{bmatrix}
\begin{bmatrix}
I & 0 \\ \conv_{L} & \obs_{L}
\end{bmatrix}
\begin{bmatrix}
u \\ \alphani
\end{bmatrix}
\,,
\label{eq: IOS representation}
\end{equation}

where $\alphani \in \real^{n} $ is the initial state.

Recall the observability problem: given length-$L$ time series of inputs and outputs, can $\alphani$ be reconstructed?  Equation (<ref>) gives a succinct answer: namely, $\alphani$ can be reconstructed if and only if $\obs_{L}$ has full column-rank. The minimum $L$ so that $\obs_{L}$ has full rank $n$ equals the lag $\ell$ of a minimal kernel representation.

As readily deducible from (<ref>) and formalized in <cit.>, in a representation-free setting, the  initial condition $\alphani$ for a trajectory $w \in \bv_{L}$, can be estimated via a  prefix trajectory $\wini = \bigl(w(-\Tini+1),  \dots, w(-1), w(0) \bigr)$  of length $\Tini \geq \ell$ so that the concatenation $\wini \wedge w \in \bv_{\Tini+L}$ is a valid trajectory.

Hence, an LTI system is characterized by the complexity parameters $(	q,m,n,\ell)$, and we denote the corresponding  class of LTI systems by $\lti_{m,\ell}^{q,n}$: namely, LTI systems with $m$ inputs, $q-m$ outputs, minimal state dimension $n$, and minimal lag $\ell$.

The following lemma characterizes  the dimension of $\bv_{L} \in \lti_{m,\ell}^{q,n}$ in terms of the complexity parameters $(	q,m,n,\ell)$.

Let $\bv \in \lti_{m,\ell}^{q,n}$. Then $\bv_{L}$ is a subspace of $ \real^{qL}$, and for $L \geq \ell$ its  dimension is $mL + n$.

Due to linearity of $\bv$, $\bv_{L} \subset  \real^{qL}$ is a subspace. To show that the dimension of $\bv_{L}$ equals $mL + n$ for $L \geq \ell$, we appeal to a minimal state-space representation of $\bv$ — a state-space-independent proof is in <cit.>. We have  $w = P \left[\begin{smallmatrix} u \\ y\end{smallmatrix}\right] \in \bv_{L}$ if and only if (<ref>)  holds

for some $\alphani \in \real^{n}$. Since the representation is minimal,  $\obs_{L} \in \real^{(q-m)L \times n}$ is of full column-rank for $L \geq \ell$. Therefore, the matrix $
\left[\begin{smallmatrix}
I & 0 \\ \conv_{L} & \obs_{L}
\end{smallmatrix}\right]
\in \real^{qL \times (mL + n)}$
is of full rank $mL+n$ for  $L \geq \ell$ and forms a basis for $\bv_{L}$. Thus, $\bv_{L}$ has dimension $mL+n$.

All forthcoming results assume known complexity $(q,m,n,\ell)$.
When only data and no prior information is available, it is reasonable to assume upper bounds on $(q,m,n,\ell)$. In this case, the anticipated dimension of $\bv_{L}$ is at most $mL + n$, and the forthcoming rank equalities the behavior dimension should be replaced by  inequalities.

§.§ Image Representation of Restricted Behavior

The restricted behavior $\bv_{L}$, the set of all trajectories of length $L$, can be described by a kernel or state-space representation. As an interesting alternative, we recall  the  image representation of $\bv_{L}$ by a data matrix of a time series.

Consider the sequence $w = \bigl(w(1),w(2),\dots,w(T)\bigr)$ with elements  $w(i) \in \real^{q}$, and define the (block) Hankel matrix $\Han_{L}(w)
\in \real^{qL \times (T-L+1)}$ of depth $L$, for some $L \leq T$, as

\begin{equation*}
\Han_{L}(w) =
\begin{bmatrix}
w(1) & w(2) & \dots & w(T-L+1)
\\
w(2)  & w(3) & \dots & w(T-L+2)
\\
\vdots & \vdots  & \ddots & \vdots
\\
w(L) & w(L+1) & \dots & w(T)
\end{bmatrix}
%\in \real^{Lq \times (T-L+1)}
\,.
\end{equation*}

A result due to [27] that became known as the Fundamental Lemma  offers an image representation of the restricted behavior in terms of the column span of a data Hankel matrix. We present a necessary and sufficient version here assuming:

*  rank$\left(\Han_L(w)\right)=mL + n$.

Consider an LTI system $\bv \in \lti_{m,\ell}^{q,n}$ and an associated trajectory  $w = \bigl(w(1), w(2),$ $ \dots, w(T)\bigr) \in \real^{qT}$. The following are equivalent for $L > \ell$:

\begin{equation*}
\text{colspan} \left(\Han_{L}(w) \right) = \bv_{L}
\quad\Longleftrightarrow\quad
\text{Assumption \ref{ass:L+n pe}}
\end{equation*}

In words, the Hankel matrix $\Han_L(w)$ composed of a single $T$-length  trajectory parametrizes all $L$-length trajectories if and only if rank$\left(\Han_L(w)\right)=mL + n$. A plausible reasoning leading up to Lemma <ref> is that every column of $\Han_L(w)$ is a trajectory of length $L$, and the set of all such trajectories has at most dimension $mL+n$; see Lemma <ref>.
Lemma <ref> extends the original Fundamental Lemma  <cit.> which requires input/output partitioning, controllability, and persistency of excitation of order $L+n$ (i.e., $\Han_{L+n}(u)$ must have full row rank) as sufficient conditions. Lemma <ref> also extends to mosaic Hankel, Page, and trajectory matrices <cit.>.

It is debatable whether the image representation via the Hankel matrix $\Han_{L}(w)$ should be called a “model”, as it is readily available from raw data. Hence, we call $\text{colspan} \left(\Han_{L}(w) \right)$ a data-driven representation of $\bv_{L}$ and reserve the term “model” for parametric (kernel or state-space) representations. Models are useful  for many reasons: first and foremost the availability of powerful analysis and design methods. Another readily discernible advantage is that models are vastly compressed compared to the image representation, and the latter holds only on finite horizons unless trajectories are weaved together [31]; see also Remark <ref>.

§ DIRECT AND INDIRECT DATA-DRIVEN CONTROL

We present different data-driven control formulations along with assumptions under which the formulations are consistent. These assumptions are used only for consistency statements and not for our main results, but they will prove insightful.

§.§ Optimal Control Problem

Given a plant with plant behavior $\bv^P \in \lti_{m,\ell}^{q,n}$, a $\Tini$-length prefix trajectory $\wini = \bigl(w(-\Tini+1),  \dots, w(0) \bigr) \in \bv_{\Tini}$, a $\Tr$-length reference trajectory $w_r\in\real^{q\Tr}$ in a reference behavior $\bv^R$,
and a set of admissible trajectories $\mathcal W \subset \real^{q\Tr}$,
consider the finite-time optimal control problem

w ∈𝒲(w-w_r) C:

∧w ∈^P_+

For  $\Tini \geq \ell$ the prefix trajectory $\wini$ implicitly sets the initial condition for the optimal control problem (<ref>); see Section <ref>. In case of uncertain initial condition, the prefix $\wini$ can be made a decision variable and included via a penalty term in the cost; c.f., [40, 39, 41, 43, 42]. We refrain from such extensions here.

Typically, the cost $\ctr:\, \mathbb R^{q\Tr} \to \real_{\geq 0}$ includes a running and a terminal cost. The set $\mathcal W \subset \real^{q\Tr}$ captures constraints  on admissible trajectories (e.g., capturing input saturation). We denote a minimizer (if it exists) of the optimization problem $\boldsymbol{C}$ in (<ref>) by $w^{\star}_{C}$.
We make the following regularity assumptions:

*  $\ctr:\, \mathbb R^{q\Tr} \to \real_{\geq 0}$ is a convex function that achieves its minimum when $w = w_{r}$;  $\mathcal W \subset \real^{q\Tr}$ is closed, convex, and non-empty; and $(\real^{q\Tini} \oplus \mathcal W) \cap \bv^{P}_{\Tini+\Tr}$ is non-empty.

The last assumption ensures that $\mathcal W$ is viable, i.e., a trajectory of $\bv^{P}$ originating anywhere can be contained within $\mathcal W$ for $L$ steps.
Problem (<ref>) is thus convex with closed, convex, and non-empty feasible set due to Assumption <ref> and because ${\bv^{P}_{\Tini+\Tr}}$ is a subspace; see Lemma <ref>. Under further standard assumptions existence and uniqueness of a (global) minimum can be assured, but we do not impose further structure.

For problem (<ref>), we do not necessarily assume $\bv^P = \bv^R$, since we often ask systems to track non-plant behavior (e.g., steps). Likewise, we  generally do not assume feasibility: $w_{r} \in \mathcal W$.  However, such assumptions connect to model reference control and allow to state  consistency results as presented next.

*  $\wini \wedge w_{r} \in {(\real^{q\Tini} \oplus \mathcal W)} \cap \bv^{P}_{\Tini+\Tr}$, i.e., the reference $w_{r} \in \bv^R_{\Tr}$ is compatible with  the prefix trajectory $\wini$, the plant $\bv^{P}$, and the constraints $\mathcal W$.

Under Assumptions <ref> and <ref>, the minimum of the  control problem $\boldsymbol{C}$ in (<ref>) is   achieved for $w^{\star}_{C}=w_{r}$.

Fact <ref> (and similar consistency results later) follows since $w^{\star}_{C}=w_{r}$ is feasible and achieves the minimum of the cost.
Fact <ref> (and consistency Assumption <ref>) serve to establish ground-truth for comparing different problem formulations.

Problem (<ref>) becomes a “classical” control problem if  a parametric model for the plant $\bv^P$ is available. The latter is usually obtained from data through system identification.

§.§ Indirect Data-Driven Control via System Identification

Given a $\Td$-length trajectory $w_{d} \in\real^{q\Td}$ as  identification data,  conventional system identification and control consists of three steps. The first step, model class selection, amounts to choosing the set of candidate models, e.g., $\lti_{m,\ell}^{q,n}$ specified by the complexity $(q,n,m,\ell)$.

The second step, model fitting, chooses an element from the model class that fits the data best in some specified sense, e.g., distance between data $w_{d}$ and model $\bv$.

This step is often synonymous to learning a parametric model (e.g., PEM), though some classic (e.g., ETFE) and modern (e.g., kernel-based) methods are non-parametric and by-pass the model order selection; see [2] for a review (and the acronyms). However, for control design the non-parametric models again have to be projected on a behavior in $\lti_{m,\ell}^{q,n}$.

Both approaches can be abstracted as

ŵ_d,(ŵ_d-w_d) ID:

ŵ_d ∈_ ,  ∈_m,ℓ^q,n .

It is useful to think of the identification loss $\cid: \mathbb R^{q\Td} \!\to\! \real_{\geq 0}$ as a distance. Given the data $w_{d}$, problem (<ref>) seeks the closest
LTI behavior  within the class $\lti_{m,\ell}^{q,n}$, i.e., the closest subspace with dimensions as in Lemma <ref>.
We denote a minimizer of (<ref>) by $\bigl(\hat w^{\star}_{d,ID},\widehat\bv_{ID}^{\star}\bigr)$ and assume the following about

*  $\cid(\cdot)$ achieves its minimum when $\hat w_{d} \!=\! w_{d}$.

Note that existence and uniqueness of minimizers of (<ref>) does not only hinge upon the regularity of cost and constraint functions, but also on the data. In general, identification problems are non-convex. For now we keep problem (<ref>) abstract and general and resort to more specific formulations in Section <ref>.

Exact identification of the true system requires exact data $w_{d} \in \bv^P_{\Td}$ and an identifiability assumption <cit.> which assures that $\bv^P$ can be  recovered from $w_{d}$:

*  $w_{d}   \in \bv^P_{\Td}$, i.e., $w_{d}$ is a valid trajectory of $\bv^P_{\Td}$	; and

*  rank$\left(\Han_{\ell+1}(w_{d})\right)=m(\ell+1) + n$.

Under assumptions <ref>–<ref>, the minimum value of the system identification problem $\boldsymbol{ID}$ in (<ref>) is achieved for $\hat w_{d,ID}^{\star}=w_{d}$ and $\widehat\bv^{\star}_{ID} = \bv^P$.

We again note that the (arguably strong) Assumptions  <ref>, <ref>, and <ref> are used only for consistency statements (such as Fact <ref>) and not for our later main results and simulations.

Finally, equipped with an identified behavior $\widehat\bv^{\star} \in \lti_{m,\ell}^{q,n}$, the third step is certainty-equivalence control: solve the optimal control problem (<ref>) subject to the identified model:

w ∈𝒲( w-w_r)

∧w  ∈^⋆_+

In (<ref>), $\ctr( w-w_{r})$ is merely a surrogate (predicted) control error since  $ w \in \widehat\bv^{\star}$, the identified model, rather than $ w \in \bv^{P}$.

Putting both the system identification (<ref>) and  certainty-equivalence control (<ref>) together, we arrive at  indirect data-driven control formulated as the bi-level problem

w ∈𝒲(w-w_r) BL:

∧w  ∈^⋆_+

^⋆ ∈_ŵ_d,      (ŵ_d-w_d)

ŵ_d ∈ _ ,

The bi-level problem structure in (<ref>) reflects the sequential system identification and control tasks, that is, first a model is fitted to the data in the inner  identification problem before the model is used for control in the outer problem.

We denote a minimizer for  the inner problem of (<ref>) by $\bigl(\hat w^{\star}_{d,BL},\widehat \bv^{\star}_{BL}\bigr)$ and a minimizer for the outer problem of (<ref>) by $w^{\star}_{BL}$.

The bi-level formulation (<ref>) is only the tip of the iceberg, and the overall design may feature further nested levels, e.g., optimization of the model selection hyper-parameters $(n,\ell)$, uncertainty quantification, etc. We deliberately neglect these levels here and focus on identification and control.

Since our ultimate interest is control, we treat models in a disregarding manner, i.e., they serve merely an auxiliary purpose. Of course, models are desired for other reasons: system design, analysis, the reasons in Remark <ref>, etc.

Under suitable consistency assumptions, the sequential system identification and control approach in (<ref>) is optimal.

Consider the optimal control problem $\boldsymbol{C}$ in (<ref>) and the bi-level problem $\boldsymbol{BL}$ in (<ref>). Then

* under Assumptions <ref>–<ref>, the bi-level problem $\boldsymbol{BL}$ reduces to the optimal control $\boldsymbol{C}$; and

* under the additional Assumptions <ref> and <ref>, the minimum value of the  bi-level problem $\boldsymbol{BL}$ is achieved for $\hat w_{d,BL}^{\star}=w_{d}$, and $\widehat\bv^{\star}_{BL} = \bv^P$, $w^{\star}_{BL}= w_{r}$.

The first statement echos the “model as well as possible” paradigm and a separation of control and identification, albeit in a simple setting; see <cit.> for further reading.

§.§ Direct Data-Driven Control via the Image Representation

The direct data-driven control approach pursued here hinges upon the Fundamental Lemma <ref>. A direct corollary of the latter is that the prediction and estimation trajectories have to be within the column span of the data Hankel matrix.

Assume that Assumptions <ref> and <ref> hold with $L$ replaced by $\Tini+L$, then the optimal control problem $\boldsymbol{C}$ in (<ref>) is equivalent to

w ∈𝒲(w-w_r) D:

∈ colspan (_+(w_d) ),

i.e., the minimizers and minima of (<ref>) and (<ref>) coincide.

Under Assumptions <ref> with $L$ replaced by $\Tini+L$, <ref>, <ref>, and <ref>, the minimum value of (<ref>) is  achieved for $w^{\star}_{D} = w_{r}$.

It is instructive to compare the sample complexity of direct and indirect approaches (<ref>) and (<ref>).

Due to Assumption <ref>,  (<ref>) requires more data than the identification Assumption <ref>. This discrepancy is due to (<ref>) seeking a multi-step predictor, whereas identification (<ref>) seeks a single-step predictor to be applied recursively.
By weaving multiple trajectories of length $\ell+1$, Assumption <ref> can be eased so that the data lengths coincide;   see <cit.>.

In comparison, with system identification, the model order selection is implicit in Assumption <ref> and encoded in the rank of the Hankel matrix $\Han_{\Tini+\Tr}(w_{d})$ – at least, for exact data $w_{d} \in \bv^{P}_{\Td}$. If the data $w_{d}$ is noisy, then $\Han_{\Tini+\Tr}(w_{d})$ likely has full rank, and the constraint  of (<ref>) is vacuous. Thus, $w=w_{r}$ uniquely minimizes the surrogate control error, but the realized control error may be arbitrarily different. In short, certainty-equivalence can fail arbitrarily poorly in direct data-driven control, and the direct approach has to be robustified. This is a major difference with the indirect (first identify, then control) approach (<ref>): one purpose of identification is to filter noisy data by projecting on a deterministic behavior.

To go beyond certainty equivalence, the DeePC approaches [40, 42, 43, 44, 41, 45, 46] reformulate the constraint in (<ref>) as $\text{col}(\wini,w) = \Han_{\Tini+\Tr}(w_{d})  g$ for some $g$ and add a robustifying regularizer.

w ∈𝒲,g(w-w_r) + λ·h(g)D_λ:

w = _+(w_d)  g

To provide an intuition, every column of $\Han_{\Tini+\Tr}(w_{d})$ is a trajectory of $\bv^{P}_{\Tini+\Tr}$, and the decision variable $g$ linearly combines these columns for the optimal trajectory $w$ – consistent with the prefix trajectory $\wini$ and  regularized by $h(g)$.
The regularization function $h(\cdot)$ and parameter $\lambda$ are nonnegative. Choices for $h(\cdot)$ are one-norms [41], two-norms [44], squared two-norms [40, 45], or arbitrary $p$-norms [42, 43, 46].

The regularizers can be related to robust optimization formulations in  deterministic [44, 45, 46] or stochastic settings [42, 43], where $\lambda$ is a design parameter specifying  the size of the assumed uncertainty set. The regularized formulation (<ref>) has proved itself in  practical (nonlinear) control systems [46, 47, 48, 49, 45].

§ BRIDGING DIRECT & INDIRECT APPROACHES

§.§ Multi-Objective Data-Driven Control

From an optimization perspective it is natural to lift the bi-level problem (<ref>) to a multi-criteria problem simultaneously optimizing for identification and control objectives. Using weighted sum scalarization, the multi-criteria problem is

w ∈𝒲,ŵ_d,γ·(ŵ_d-w_d)  +  (w-w_r) MC_γ:

∧w  ∈_+ ,  ŵ_d ∈_ ,

where the trade-off parameter $\gamma \geq 0$ traces the Pareto front between the identification and optimal control objectives.

The multi-criteria problem (<ref>) can be interpreted as fitting a model $\widehat\bv$ simultaneously to two data sets: the identification data $w_{d}$ and the  reference $w_{r}$. From a control perspective, the identification criterion biases the solution $w \in \widehat\bv$ to adhere to the observed data $w_{d}$ rather than merely matching the to be tracked reference $w_{r}$. Likewise, from the other side, the identification criterion is biased by the control objective.

In short, control and identification regularize each other, in the spirit of identification for control [7, 12, 13, 14]. A  similar  formulation has been proposed in [15] interpolating between PEM identification and a model-reference control objective. Likewise, the data-driven model reference control formulation in [20] interpolates between a direct and an indirect approach. Finally, dual control approaches consider similar multi-criteria formulations balancing exploration (for identification) and exploitation (i.e., optimal control) [16, 17, 18, 19].

We denote a minimizer of  (<ref>) by $\bigl(w^{\star}_{MC},\hat w^{\star}_{d,MC},\widehat \bv^{\star}_{MC}\bigr)$.

Under Assumptions <ref>–<ref>, for any $\gamma \geq 0$ the minimum of the parametric multi-criteria problem $\boldsymbol{MC}_{\gamma}$ is  achieved for $\hat w_{d,MC}^{\star}=w_{d}$, $\widehat\bv^{\star}_{MC} = \bv^P$, and $w^{\star}_{MC} = w_{r}$.

Different points on the Pareto front of (<ref>) have different emphasis regarding the control and identification objectives. Below we formalize that for $\gamma$ sufficiently large, the multi-criteria problem (<ref>) recovers the bi-level problem (<ref>) corresponding to sequential system identification and control.

We follow standard penalty arguments from bi-level optimization [52, 53], which are particularly tractable here since  (<ref>) is only weakly coupled: the inner problem does not depend on the decision variable $w$ of the outer problem. Assume there is a minimum (termed value function) of the inner problem:

ŵ_d,   (ŵ_d-w_d) φ  =

ŵ_d ∈   _ ,    ∈_m,ℓ^q,n .

The bi-level problem (<ref>) reads then equivalently as

∧w  ∈_+ ,  ŵ_d ∈_ ,

∈_m,ℓ^q,n ,  (ŵ_d-w_d) - φ= 0  .

At this point the reader is encouraged to review the definition and salient properties of a constraint qualification termed partial calmness [52, 53]; see the appendix. If problem (<ref>) is partially calm at a local minimizer and $\ctr(\cdot)$ is continuous, then there is $\gamma^{\star}>0$ so that, for all $\gamma > \gamma^{\star}$, then (<ref>) equals

w∈𝒲,ŵ_d,γ·| (ŵ_d-w_d) - φ|  +  (w-w_r)

∧w  ∈_+ ,  ŵ_d ∈_ ,   ∈_m,ℓ^q,n  ,

that is, the local minimizers of (<ref>) and (<ref>) coincide; see Proposition <ref>.

We now drop the absolute value (since $ \cid(\hat w_{d}-w_{d}) - \varphi \geq 0$) and the constant $\varphi$ (which in our case does not depend on the variable $w$ of the outer problem) from the objective of (<ref>) to recover problem (<ref>).

We have thus established a chain of equivalences relating the bi-level and multi-criteria problems.
We summarize our discussion below.

Consider the parametric multi-criteria problem $\boldsymbol{MC}_{\gamma}$ in (<ref>) and  the bi-level problem $\boldsymbol{BL}$ in (<ref>). Assume that the inner identification problem admits a minimum as in (<ref>), (<ref>) is partially calm at any local minimizer, and $\ctr(\cdot)$ is continuous. Then there is $\gamma^{\star}>0$ so that for $\gamma > \gamma^{\star}$ the problem $\boldsymbol{MC}_{\gamma}$ is equivalent to $\boldsymbol{BL}$, i.e., $ w^{\star}_{MC} =  w^{\star}_{BL}$ , $\hat w^{\star}_{d,MC} = \hat w^{\star}_{d,BL}$, and $\widehat \bv^{\star}_{MC} = \widehat \bv^{\star}_{BL}$.
Moreover, the optimal values of $\boldsymbol{MC}_{\gamma}$ and $\boldsymbol{BL}$ coincide up to the constant $\gamma \cdot \varphi$ with $\varphi$ defined in (<ref>).

The following comments are in order regarding partial calmness. As discussed in Proposition <ref>, partial calmness is equivalent to the constraint $ \cid(\hat w_{d}-w_{d}) - \varphi \geq 0$ serving as an exact penalty. Partial calmness is satisfied, for instance, appealing to Proposition <ref>, if the identification cost $\cid(\cdot)$ can be phrased as a distance (see the discussion following the identification problem (<ref>)) and $\ctr(\cdot)$ is Lipschitz continuous over the feasible set, e.g., the feasible set is either compact (due to constraints) or the control performance is measured by a norm or Huber loss. The Lipschitz constant then serves as a lower estimate for $\gamma^{\star}$.

A non-Lipschitz cost requires $\gamma \to \infty$ as a sufficient condition. Note that for $\gamma \to \infty$ Proposition <ref> holds without assumptions, since (<ref>) is merely an indicator function reformulation of  (<ref>). Our relaxations in the next sections will, among others, drop the requirement on $\gamma$ sufficiently large  as well as the LTI complexity specification $\widehat\bv \in \lti_{m,\ell}^{q,n}$.

Even if the identification (<ref>) is convex the multi-criteria problem (<ref>) is not, since it simultaneously optimizes over the to-be-identified model $\bv$ and  the to-be-designed trajectory $w$. This can be spotted in a kernel representation: the constraint $ w \in \widehat\bv_{\Tr}$ takes the form $\widehat R(\sigma) w = 0$, where both $\widehat R$ and $ w$ are variables. Other representations lead to the same conclusions.

Consider the multi-criteria problem (<ref>) and a kernel representation of the to-be-identified behavior: $\widehat\bv = \text{kernel}(\widehat R(\sigma))$. Then the feasible set of (<ref>)  is not convex.

We believe that the multi-criteria problem is interesting in its own right: studying its Pareto front and choosing an optimal trade-off parameter may possibly yield superior  performance.

Our problem setup thus far was conceptual rather than practically useful. Below, we consider concrete problem formulations and turn our conceptual insights into concise results.

§.§ Bridging Towards Subspace Predictive Control (SPC)

We explain SPC  from the perspective of the Fundamental Lemma <ref> stating that any trajectory $\wini \wedge w \in \bv_{\Tini+\Tr}^{P}$ lies in $\text{colspan} \left(\Han_{\Tini+\Tr}(w_{d}) \right)$. Recall that $\wini$ is a prefix trajectory of length $\Tini \geq \ell$ setting the initial condition, and $w$ is a future trajectory  of length $\Tr>1$ to be designed via optimal control. Accordingly, permute and partition $w$ and the Hankel matrix

\begin{equation*}
\begin{bmatrix} \wini \\ w \end{bmatrix} \sim \begin{bmatrix}
\uini \\ u \\ \hline \yini  \\ y
\end{bmatrix}
%w_{d} = \begin{bmatrix}
%u_\text{ini}^{d} \\ y_\text{ini}^{d} \\ u^{d} \\ y^{d}
\Han_{\Tini+L}(w_{d})
\sim
\begin{bmatrix}
\Up \\ \Uf \\\hline \Yp  \\ \Yf
\end{bmatrix}
\begin{bmatrix}
\Han_{\Tini + \Tr}(u_{d}) \\\hline \Han_{\Tini + \Tr}(y_{d})
\end{bmatrix}
\,,
\end{equation*}

where $\uini \in \real^{m\Tini}$, $\yini \in \real^{(q-m)\Tini}$, and $\sim$ denotes similarity under a coordinate  permutation.

The subscripts “p” and “f” are synonymous to “past” and “future”.
We seek a linear model, i.e., a matrix $K$, relating past and future as

\begin{equation}
y =
\underbrace{
\begin{bmatrix}
K_{p} & \vline& K_{f}
\end{bmatrix}}_{=K}
\cdot \begin{bmatrix}
\uini \\ \yini \\\hline u
\end{bmatrix}
\label{eq: ARX transition model}
\,.
\end{equation}

The multi-step predictor $K$ is found from  Hankel matrix data by means of the least-square criterion  <cit.>

K	-  K ·

where $\|\cdot\|_{F}$ is the Frobenius norm.
Via the Moore-Penrose inverse, the solution of (<ref>) is the classic SPC predictor [32]

\begin{equation}
K =
\Yf
\begin{bmatrix}
\Up \\ \Yp \\\Uf
\end{bmatrix}^{\dagger}
\,.
\label{eq: SPC predictor}
\end{equation}

It is insightful to compare equation (<ref>) and the matrices $K_{p},K_{f}$ to equation (<ref>) and the extended observability and impulse response matrices $\obs_{L}$ and $\conv_{L}$, respectively. One realizes that for exact data, (<ref>) is an ARX model with rank$(K_{p})=n$ assuring LTI behavior of desired complexity and a lower block-triangular zero pattern of $K_{f}$ assuring causality. For inexact data,  LTI behavior of desired complexity is promoted by low-rank approximation (typically via singular-value thresholding of $K_{p}$) [32]; and one aims to gain causality by heuristically thresholding $K_{f}$ towards a desired zero pattern
<cit.>, <cit.>. The causality requirement can also be omitted for offline or receding horizon control, but it is useful to condition the data on the set of causal models. These steps bring the linear relation (<ref>) half-way towards an LTI model. Though a model has further structure, e.g., $K_{f}$ is Toeplitz, and the entries of $K_{p}$ and $K_{f}$ are coupled; see (<ref>).

Hence, in this case the identification problem (<ref>) is relaxed to the single, monolithic, and non-convex program

K	-  K ·

K = 	 K_p          K_f
K_f  lower-block triangular
rank( K_p) = n

where the lower-block triangular specification means that all entries above the diagonal $(q-m)\times m$ blocks equal zero.

We obtain a parametric version of the indirect data-driven approach (<ref>), where $\wini \wedge w  \in \widehat\bv^{\star}_{\Tini+\Tr}$ and $w \in \mathcal W = \mathcal U \times \mathcal Y$ are replaced by (<ref>) and $(u,y) \in \mathcal U \times \mathcal Y$, respectively:

u ∈𝒰,y ∈𝒴( [ y - y_r
u - u_r    ] )

y = K^⋆ ·

K^⋆ ∈_K
-   K ·

K = 	 K_p          K_f
K_f  lower-block triangular
rank( K_p) = n

We stress that (<ref>) is generally not an equivalent reformulation of (<ref>) since the inner identification does not necessarily lead to an LTI model; see the comments following equation (<ref>).

For comparison, consider  also an instance of the direct regularized problem (<ref>) with regularizer $h(g) = \|(I-\Pi)g\|_{p}$:

u ∈𝒰,y ∈𝒴,g( [ y - y_r
u - u_r    ] )  +  λ·(I-Π)g_p

g =

Here, $\|\cdot\|_{p}$ is any $p$-norm, $\Pi =
\left[\begin{smallmatrix}
\Up \\ \Yp \\\Uf
\end{smallmatrix}\right]^{\dagger}\!
\left[\begin{smallmatrix}
\Up \\ \Yp \\\Uf
\end{smallmatrix}\right]$, and $(I-\Pi)$ is an orthogonal projector on the kernel of the first three block-constraint equations.

The proof of Theorem <ref> will later show that this regularizer is in fact  induced by the least-square identification (<ref>), i.e., $\|(I-\Pi)g\|_{p}=0$ if and only if the least-square criterion is minimized. Hence, it robustifies the problem akin to least squares. We state the following consistency result.

Under Assumptions <ref> with $L$ replaced by $\Tini+L$, <ref>, <ref>, and <ref>, for any $\lambda \geq 0$ the minimum of the regularized  problem (<ref>) is achieved for $y^{\star} = Y_{f}g^{\star} = y_{r}$ and $u^{\star} = U_{f}g^{\star} = u_{r}$, where $\|(I-\Pi)g^{\star}\|_{p}=0$.

Fact <ref> may not appear insightful at first glance, but it highlights an important fact. The projection-based regularizer $h(g)= \|(I-\Pi)g\|_{p}$ is consistent since it penalizes only the homogenous solution to the constraint equations (<ref>) and does not affect the variables $(u,y)$. In comparison,
the conventional norm-based regularizer $h(g) = \|g\|_{p}$ is not consistent: it penalizes the heterogeneous solution of the constraint equations in (<ref>) and thus also $(u,y)$. Hence, even with ideal consistency Assumptions <ref>, <ref>, <ref>, and <ref> in place, the norm-based regularizer $h(g) = \|g\|_{p}$ with $\lambda \neq 0$ does not lead to the ground-truth solution $y^{\star} = Y_{f}g^{\star} = y_{r}$, $u^{\star} = U_{f}g^{\star} = u_{r}$; see also Remark <ref>.

The following is the main result of this subsection.

Consider the indirect data-driven control problem (<ref>) and the direct data-driven control problem (<ref>) parameterized by $\lambda \geq 0$. Let Assumption <ref> hold and assume that $\ctr(\cdot)$ is Lipschitz continuous. For $\lambda$ sufficiently small,  (<ref>) is a convex relaxation of  (<ref>), that is,

$(i)$ (<ref>) is convex,

$(ii)$  any feasible $(u,y)$ in (<ref>)  is feasible for (<ref>), and

$(iii)$  the optimal value of (<ref>) lower-bounds that of (<ref>).

First, we perform a convex relaxation by dropping the rank and block-triangularity constraints in (<ref>). Second, observe that the explicit solution of the inner problem, the predictor (<ref>), is equivalently derived as least-norm solution

\begin{align*}
y = \Yf g^{\star}
\;\text{ where}\;
&g^{\star} = \argmin_{g} \| g\|_{2}
\nonumber
\\&\text{subject to}
\begin{bmatrix}
\Up \\ \Yp \\\Uf
\end{bmatrix}
g = \begin{bmatrix}
\uini \\ \yini \\ u
\end{bmatrix}
\,.
%\label{eq:least-norm form}
\end{align*}

u ∈𝒰,y ∈𝒴( [ y - y_r
u - u_r    ] )
y = g^⋆

g^⋆ ∈_ g    g_2

g =

We now follow the arguments from Section <ref> to reduce the bi-level problem (<ref>) to a single-level multi-criteria problem.

As in (<ref>), the inner problem can be replaced by a constraint assuring that it achieves its minimum.

Here, we add an orthogonality constraint to the constraints of the inner problem:

\begin{equation*}
%0 = g-g^{\star} = g - \Pi g
\begin{bmatrix}
\Up \\ \Yp \\\Uf
\end{bmatrix}
g = \begin{bmatrix}
\uini \\ \yini \\ u
\end{bmatrix}
\quad\text{and}\quad
0 = \| (I - \Pi) g \|_{p}
\end{equation*}

The orthogonality constraint $0 = \| (I - \Pi) g \|_{p}$ poses the inner optimality constraint as the distance to the subspace containing the minimizers of the inner problem.

Retaining all constraints, (<ref>) can then be formulated as the single-level problem

u ∈𝒰,y ∈𝒴,g( [ y - y_r
u - u_r    ] )

g =

(I - Π) g _p = 0

We now apply Proposition <ref>, lift the distance constraint $\| (I - \Pi) g \|_{p} = 0$ to the objective, and recover  problem (<ref>) with $\lambda$ larger than the Lipschitz constant of $\ctr(\cdot)$.

Hence, (<ref>) is equivalent to (<ref>) for $\lambda$ sufficiently large. Our final convex relaxation is to choose $\lambda$ small rather than large. Namely, from the view-point of the objective: it lowers the cost; or from the bi-level viewpoint: it turns the inner optimality constraint into a weaker sub-optimality constraint, i.e., we allow for solutions satisfying $\| (I - \Pi) g \|_{p} \geq  0$.

Conclusion $(i)$ now follows since (<ref>) is convex; $(ii)$ follows since we have only enlarged the feasible set when passing from (<ref>) to (<ref>); and $(iii)$ follows due to the enlarged feasible set, since the costs of (<ref>) and (<ref>)  coincide, and since (<ref>) is a relaxation of (<ref>) if $\lambda$ is not sufficiently large.

First, we summarize the salient arguments to pass from indirect to direct data-driven control: we
relaxed problem (<ref>) by dropping causality (block-triangularity) and LTI complexity (rank) specifications, replaced the least-square criterion (<ref>) by the equivalent least-norm formulation (<ref>),
and lifted the problem from bi-level to multi-criteria, where the least-square objective induces the regularization $\|(I-\Pi)g\|_{p}$.
For equivalence to the least-square objective, the proof requires $\lambda$ larger than the (global) Lipschitz constant of $\ctr(\cdot)$, similar to  robustification-induced regularizations [42, 43]. If  $\ctr(\cdot)$ is only locally Lipschitz, e.g., in case of a quadratic cost, then choosing a finite (small) $\lambda$ is a relaxation that allows the predicted trajectory to not adhere to the least-square fit of the data. Though as we will see in Section <ref>, its effect is minor for $\lambda$ not overly small.

Second, continuing on the magnitude of $\lambda$: For exact data and under consistency assumptions, (<ref>) achieves the exact minimizer for  any $\lambda \geq 0$; see Fact <ref>. When departing from these ideal assumptions,
the least-square fit of the data is enforced only for $\lambda$ sufficiently large. Generally, $\lambda$ should be regarded as a tune-able hyper-parameter chosen by the designer to control how much the predicted trajectory should adhere to the data (versus the control objective) and to ultimately improve the realized performance. The proof of Theorem <ref> suggests a sufficiently large value, which is also confirmed by our later empirical findings (see e.g.  Figure <ref>).

Third, the regularization based on the projector $\|(I-\Pi)g\|_{p}$ differs from the standard $p$-norm regularizers $h(g) = \|g\|_{p}$ [44, 42, 43] (or squared 2-norms $\|g\|_{2}^{2}$ [45, 40]). Actually, it is this projection  which recovers
the least-square criterion (<ref>). In contrast, norm-based regularizers $\|g\|_{p}$ are not consistent and bias the optimal solution $(u^{\star},y^{\star})$; see Remark <ref>.

This is undesirable from an identification perspective: the regularizer should induce a least-square fit of the data. While for small values of $\lambda$  both regularizers have a similar  effect, for sufficiently large $\lambda$  the identification-induced regularizer $\|(I-\Pi)g\|_{p}$ demonstrates a superior performance; see Figure <ref> later.

Fourth, our proof strategy reveals an entire class of regularizers. In fact, we can choose any $p$-norm $\|(I-\Pi)g\|_{p}$, use more general penalty functions such as the (squared) merit functions in [53], or attack problem (<ref>) with other penalty or augmented Lagrangian methods. These degrees of freedom reflect the intuition that the Pareto-front of (<ref>) is invariant under certain (e.g., monotone) transformations of objectives such as taking squares; see <cit.> for a formal reasoning. For our later simulations in Section <ref>, we choose the computationally attractive regularization $\|(I-\Pi)g\|_{2}^{2}$.

Fifth and finally, our proof arguments are obviously “qualitative” crossing out rank and causality constraints similar to most SPC implementations and using non-quantifiable “sufficiently large” reasoning. Hence, the convex relaxation (<ref>) of  (<ref>) should not be expected to be tight.
Nevertheless, the formulation (<ref>) (without projector) has proved itself in many case studies and often outperforms (<ref>), as testified in [46, 47, 48, 49, 45]. Section <ref> will compare the different formulations.

§.§ Bridging Towards Structured Low-Rank  Approximation

We now present an entirely non-parametric problem formulation, namely a version of subspace identification based on structured low-rank approximation [39], and we relate the resulting bi-level problem to direct data-driven control (<ref>).

Given the model class $\lti_{m,\ell}^{q,n}$, we project the identification data $w_{d} \in \real^{qT}$ on $ \widehat\bv_{\Tini+\Tr} \in \lti_{m,\ell}^{q,n}$.  By Lemma <ref>, the latter set is characterized by all trajectories $\hat w \in \real^{q(\Tini+\Tr)}$ so that the associated Hankel matrix satisfies $\text{rank} \left(\Han_{\Tini+\Tr}(\hat w) \right) \leq m(\Tini+\Tr)+n$ for $(\Tini+\Tr) > \ell$. An implicit assumption is, of course, $\Td \gg \Tini+\Tr$: the identification data is much longer than the estimation plus control prediction horizons.

In presence of noise,  $\Han_{\Tini+\Tr}( w_{d}) $ will not have low rank and has to be approximated by a low-rank matrix in an identification step.
Thus, the identification problem (<ref>) reads as

rank (_+(ŵ_d) ) ≤ m(+)+n.

Problem (<ref>) is to be read as low-rank approximation problem: given the identification data assorted in a Hankel matrix $\Han_{\Tini+\Tr}(w_{d})$, we seek the closest sequence $\hat w_{d}$ so that the Hankel matrix $\Han_{\Tini+\Tr}(\hat w_{d})$ has rank no more than $m(\Tini+\Tr)+n$.

Since $\hat w_{d} \in \widehat \bv_{T}$, we have  $\text{rank} \left(\Han_{\Tini+\Tr}(\hat w_{d}) \right) {\leq}  m(\Tini+\Tr)+n$. Since also  $\wini \in \widehat\bv_{\Tini}$ and $w \in \widehat\bv_{\Tr}$, we conclude

\begin{equation*}
\text{rank} \left(\left[\Han_{\Tini+\Tr}(\hat w_{d}) \,~\, \text{col}(\wini,w) \right]\right) {\leq} m(\Tini+\Tr)+n \,.
\end{equation*}

Assuming that $\text{rank} \left(\Han_{\Tini+\Tr}(\hat w_{d})\right) = m(\Tini+\Tr)+n$, which is generically the case, $\Han_{\Tini+\Tr}(\hat w_{d}) g = \text{col}(\wini,w) $ for some vector $g$. Hence, the bi-level problem (<ref>) takes the form

w ∈𝒲,g(w-w_r)

= _+(ŵ_d^⋆)g

ŵ_d^⋆ ∈_ŵ_d   ( ŵ_d-w_d)

(_+( ŵ_d)) = m(+)+n

Consider the indirect data-driven control problem (<ref>) and the direct data-driven control problem (<ref>) for $h(g) = \|g\|_{1}$ and parameterized by $\lambda\geq0$. Let Assumptions <ref> and <ref> hold. For $\lambda$ sufficiently small, (<ref>) is a convex relaxation of  (<ref>), that is,

$(i)$ (<ref>) is convex,

$(ii)$ any feasible $(w,g)$ in (<ref>) is also feasible for (<ref>), and

$(iii)$ the optimal value of (<ref>) lower-bounds that of (<ref>).

To prove the claim, one can resort to a proof strategy via the multi-criteria problem (<ref>), as in the previous section. Instead, we present a more direct approach here.

We start by massaging the rank constraint in (<ref>). First, since $ \text{rank} \left(\Han_{\Tini+\Tr}(\hat w_{d})\right) = m(\Tini+\Tr)+n$, we may without loss of generality add the constraint $\|g\|_0 \leq n+ m(\Tini+\Tr)$ to the outer problem, where $\|g\|_0$ denotes the cardinality (number of nonzero entries) of $g$. Second, we perform a convex relaxation and drop the rank constraint.
Third, another convex relaxation (popular in LASSO problems [55]) is to replace $\|g\|_0 \leq n+ m(\Tini+\Tr)$ by $\|g\|_{1} \leq \alpha$ for $\alpha>0$  sufficiently large. As a result of these three steps, (<ref>) is relaxed to

= _+(ŵ_d^⋆)g
,  g_1 ≤α

ŵ_d^⋆ ∈_ŵ_d   ( ŵ_d-w_d)

Observe that under Assumption <ref> the inner problem admits a trivial solution:  $\hat w_{d}^{\star}=w_{d}$. Thus, (<ref>) reduces to

w∈𝒲,g( w-w_r)

= _+( w_d)g
,  g_1 ≤α.

Next, we lift the 1-norm constraint to the objective

w∈𝒲,g( w-w_r) + λ·g_1

= _+( w_d)g

where $\lambda\geq0$ is a scalar weight.
In particular, for each value of $\alpha$  in (<ref>), there is $\lambda\geq0$ so that the solution of (<ref>) coincides with (<ref>), and vice versa. These equivalences are standard in $\ell_{1}$-regularized problems and follow from strong duality (applicable since $\ctr(\cdot)$ is convex and Slater's condition holds)  [55].

The precise value of $\lambda$ depends on the Lagrange multiplier of the constraint $\|g\|_1 \leq \alpha$ and thus  on the data. In either case, there is a selection of parameters so that both problems are equivalent, and choosing $\lambda$ sufficiently small is a relaxation.

Thus, we arrived at the direct data-driven control (<ref>) for $\lambda$ sufficiently small and $h(g) = \|g\|_{1}$.

Conclusion $(i)$  follows due to convexity (<ref>); $(ii)$ follows since we have enlarged the feasible set passing from (<ref>) to (<ref>); and $(iii)$ follows due to the enlarged feasible set, since the costs of (<ref>) and (<ref>)  coincide, and since (<ref>) is a relaxation of (<ref>) for $\lambda$ small.

In summary, to pass from indirect data-driven control (<ref>) to direct data-driven control (<ref>), we performed a sequence of convex relaxations effectively replacing the rank constraint of the system identification by a $\ell_{1}$-norm regularizer. Hence, the 1-norm regularizer accounts for selecting the model complexity. Similar remarks as those following Theorem <ref> on tightness of the relaxation apply to Theorem <ref>, too; see Remark <ref>.

§.§ Hybrid relaxations

Theorems <ref> and <ref> reveal the roles of the two regularizers: $\|g\|_{1}$ controls the model
complexity, whereas $\|(I-\Pi)g\|_{2}$ accounts for least-square fitting the data.

To blend the two, consider a hybrid formulation of (<ref>) and (<ref>)

w∈𝒲,g(w-w_r)   +  λ_1 ·(I-Π)g^2_2

= _+(ŵ_d^⋆)g

ŵ_d^⋆ ∈_ŵ_d   ( ŵ_d-w_d)

(_+( ŵ_d)) = m(+)+n

where $\lambda_{1} \geq 0$. Observe that this formulation is consistent:

Under Assumptions <ref> with $L$ replaced by $\Tini+L$, <ref>, <ref>, and <ref>, for any $\lambda \geq 0$ the minimum of (<ref>) and achieved for  $w^{\star} = w_{r}$ and $\|(I-\Pi)g^{\star}\|^{2}_{2}=0$.

The arguments in the previous section then lead us to

w ∈𝒲,g( w-w_r)   +  λ_1 ·(I-Π)g^2_2   +  λ_2 ·g_1

= _+( w_d)g

where $\lambda_{2} \geq 0$. We will validate the performance of the hybrid regularizer in Section <ref> below; see specifically Figure <ref>.

§.§ Possible pitfalls of relaxations

Note that the two convex relaxation results in Theorems <ref> and <ref> are trivially true in the limit when $\lambda = 0$. In fact, even the abstract multi-criteria formulation (<ref>) can be related to a  relaxation of the abstract bi-level problem (<ref>)  in the limit $\gamma = 0$. Namely, for $\gamma = 0$,  (<ref>)   reduces to

w,ŵ_d,( w-w_r)

∧w ∈_+ ,  ŵ_d ∈_ ,  ∈_m,ℓ^q,n

The variable $\hat w_{d}$ and the constraint $\hat w_{d} \in \widehat\bv_{\Td}$ can be removed, and (<ref>) amounts to matching the model $\widehat\bv$ to the reference $w_{r}$. The next result is followed by a discussion on regularizers:

Consider the indirect data-driven control (<ref>) and multi-criteria problem (<ref>) in the limit $\gamma = 0$, and let Assumption <ref> hold. Then problem (<ref>) is a  relaxation of problem (<ref>), that is,

$(i)$ any feasible $( w,\hat w_{d},\hat \bv)$ in (<ref>) is also feasible for (<ref>),

$(ii)$ and the optimal value of (<ref>) lower-bounds that of (<ref>).

Consider the equivalent formulation (<ref>) of (<ref>), and note that (<ref>) equals (<ref>) when the inner optimality constraint $\cid(\hat w_{d}-w_{d}) - \varphi = 0$ is dropped. The conclusions now follow analogously as in Theorems <ref> and <ref>.

Analogous corollaries can be stated for Theorems <ref> and <ref> for $\lambda = 0$.
Given such results, one may wonder whether Theorems <ref> and <ref> are vacuous since they are trivially true for $\lambda = 0$. We offer several answers.

First, the limit $\lambda =0$ clearly leads to a better solution $w^{\star}$ (i.e., a lower surrogate tracking error) for the open-loop optimal control problem. However, this solution merely matches the reference $w_{r}$ and does not adhere to the identification data $w_{d}$ in the sense of meeting any fitting criterion. Hence, the optimal solution $w^{\star}$ may not be a trajectory of the true system behavior, and the actual realized control performance can be arbitrarily poor. Obviously, such a situation is not desirable, and one may want to regularize with a small but non-zero $\lambda$ – an observation consistent with [41, 42, 43, 44, 45, 46] albeit derived from a different perspective.

Second, Theorems <ref> and <ref> require $\lambda$ to be sufficiently small, but not  zero. According to the proofs, depending on Lipschitz constants and multipliers of the respective  problems, there is a smallest value for $\lambda$ so that the behavior $\widehat\bv$ matches (in the $\cid(\cdot)$ fitting criterion) the plant behavior $\bv^{P}$.
In [41, 42, 43, 44, 45, 46] the coefficient $\lambda$ relates to a desired robustness level. In either case, $\lambda$ can hardly be quantified a priori and without cross-validation; see also Remark <ref>.

We follow up on this set of questions in the next section.

§ NUMERICAL ANALYSIS AND COMPARISONS

We now numerically investigate the effect of the hyper-parameter $\lambda$, confirm the superiority of the regularizer $h(g) = \|(I-\Pi)g\|_{2}^{2}$, and compare direct and indirect approaches.

§.§ Choice of Regularization Parameter

We first study the parameter $\lambda$ regularizing direct data-driven control (<ref>).
Consider the benchmark single-input, single-output, 5th order, linear time-invariant system [56].
Denoting the $t$-th element of the concatenated input and output by ${w}(t)=({u}(t),{y}(t))$, the control cost was chosen as $c_\textup{ctrl}({w}-w_r) = ({w}-w_r)^{\top}W({w}-w_r)$ with reference $w_r(t) = (u_r(t),y_r(t)) = (0,\sin(2\pi t/(L-1)))$ for $t\in\{0,1,\dots,L-1\}$, prediction horizon $L=20$, $W=I_L\otimes\textup{diag}(0.01,2000)$, where $I_L$ is the $L \times L$ identity, and $\otimes$ denotes the Kronecker product. In this entire section, we disregard constraints, i.e., $\mathcal W \equiv \real^{q\Tr}$. We used a 1-norm regularizer $h(g) = \|g\|_{1}$ in (<ref>) and a prefix-trajectory of length $\Tini=5$ (see Section <ref>).

We collected one noise-free input/output time series of length $T=250$ by applying a random Gaussian input. From this noise-free data set, 100 independent noisy data sets were constructed by adding Gaussian noise with a noise-to-signal ratio of 5%. For each data set and each value of $\lambda\in(0,10^3)$, optimal control inputs were computed from (<ref>). We define the predicted error as $c_{\textup{ctrl}}(w^{\star}-w_r)$, where $w^{\star}$ is an optimizer of (<ref>). We define the realized error as $c_{\textup{ctrl}}(w_{\textup{true}}-w_r)$, where $w_{\textup{true}}$ is the realized trajectory of the system after applying the computed optimal inputs. The predicted and realized errors were converted to a percentage increase in error with respect to the ground-truth optimal performance (i.e., if the deterministic system was exactly known), and were averaged over the 100 independent data sets. The results are plotted in Figure <ref>.

It is apparent that choosing $\lambda$ too small leads to an optimistic predicted error but very poor realized performance. Furthermore, the performance is poor for  large values of $\lambda$ indicating that the regularization parameter should be chosen carefully (though a wide range delivers equally good results). These observations are consistent with those in [41, 44, 42, 43, 45, 46] and the hypotheses discussed at the end of Section <ref>.

Predicted and realized errors (relative to the ground-truth optimal performance and averaged over 100 data sets) with 1-norm regularizer $\lambda \|g\|_{1}$.

§.§ Role of Projection in Two-Norm Regularization

Theorem <ref> suggests that the identification-induced regularizer $h(g) = \|(I-\Pi)g\|_{2}^{2}$ is superior to a two-norm regularizer $h(g) = \|g\|_{2}^{2}$ if one is interested in consistency and the predicted trajectory adhering to a least-square fit of the data.

To test this hypothesis, we consider the same case study from Section <ref>
and report the averaged cost in Figure <ref>.

Comparison of the realized performance (relative to the ground-truth optimal performance and averaged over 100 data sets) for the two-norm $\|g\|_{2}^{2}$ and identification-induced regularization $\|(I-\Pi)g\|_{2}^{2}$ as function of $\lambda$.

Both regularizers perform similarly for small $\lambda$, but the identification-induced regularizer shows a superior and surprisingly constant performance for sufficiently large  $\lambda$. By the proof of Theorem <ref>, for $\lambda$ sufficiently large, the direct and indirect problems (<ref>) and (<ref>) are equivalent – up to causality and complexity constraints. Thus, a sufficiently large $\lambda$ forces the least-square fit (<ref>) and results in excellent performance independent of the specific value of $\lambda$. While there is a small window where the two-norm excels, the identification-induced regularizer shows overall much more robust performance.

Realized error (relative to the ground-truth optimal performance and averaged over 100 data sets) for a hybrid regularizer $\lambda_{1}   \|(I-\Pi)g\|_{2}^{2}  + \lambda_{2}   \|g\|_{1}$

Next we study the merits of hybrid regularization (<ref>). For the same case study Figure <ref> shows the averaged realized performance plotted over the regularization parameters. The $\{\lambda_{1}=0\}$ and $\{\lambda_{2}=0\}$ slices recover Figures <ref> and <ref>. As before, the regularizer $ \|(I-\Pi)g\|_{2}^{2}$ is more robust though a hybrid regularizer yields a minor albeit robust improvement. A closer examination of the data underlying Figure <ref> reveals that a hybrid regularization can improve up to 15% over the best results achievable with the regularizer  $ \|(I-\Pi)g\|_{2}^{2}$  only.

§.§ Effect of data length

We continue with the same case study and discuss the effect of data-length on direct and indirect methods.
For the direct method, we used the identification-induced regularizer $h(g) = \|(I-\Pi)g\|_2^2$ with sufficiently large weight $\lambda = 10000$, as indicated in Figure <ref>. For the indirect method, the inner system identification problem (<ref>) is solved using the subspace approach N4SID [26] with prefix horizon $\Tini=5$, prediction horizon $L=20$, and (correct) model-order selection $n=5$.

For our case study Lemma <ref> demands at least $T = 59$ data points. Figure <ref> below shows the beneficial effects of including more data on the realized median performance of the direct and indirect methods. The main findings are as follows: First,  both methods are asymptotically consistent. Second, the indirect method is superior in the low data regime echoing that models are compressed and de-noised  representations, see Remark <ref>. Third and finally, when an incorrect model-order $n=6$ is selected for the indirect method (resulting in an over-parameterization and thus a bias), then consistency is lost, and the direct method is superior. This effect is even more pronounced when studying the average (as opposed to the median) error due to several outliers of the indirect method.

This third point hints at a bias-variance trade-off between the direct and indirect methods, which will be studied below.

Realized median error (over 100 data sets) for the direct and indirect (with different model order selections) methods for varying amount of data.

§.§ Comparison and Bias-Variance Hypotheses

We now compare the direct and indirect approaches through two case studies. The first study evaluates the performance of both methods on the basis of “variance” error, i.e., on a linear system with noisy measurements. The second study evaluates the performance on the basis of “bias” error, i.e., on a nonlinear system with noise-free measurements.

We expect the direct method to perform better on the nonlinear system since the indirect method erroneously selects a linear model class thus leading to a larger “bias” error. On the other hand, we expect the indirect method to perform better on the linear system with noisy outputs since the identification step filters  noise thus  leading to a lower “variance” error.

§.§ Comparison: Stochastic Linear System

Consider the same case study as in the  Section <ref>, i.e., same LTI system, cost, and reference.

We collected data for varying levels of noise-to-signal ratio, i.e., we considered measurements that were affected by Gaussian noise with noise-to-signal ratio in the set $\{0\%,1\%,\dots,15\%\}$. For each noise-to-signal ratio, $T=250$ input/output data samples were collected by applying a random Gaussian input. This data was then used for both the direct and indirect methods.

For the indirect method, the inner system identification problem (<ref>) is again solved using
N4SID [26] with prefix horizon $\Tini=5$ and prediction horizon $L=20$. Equipped with a  (correct) 5th-order identified model, optimal control inputs are computed by solving (<ref>). The indirect method was compared to the direct method (<ref>), with $h(g) = \|g\|_{1}$, $\Tini=5$, and $\lambda=27$. The hyper-parameters of both methods were kept constant for all simulations below and chosen to give good realized control  performance for all noise-to-signal ratios.

For both methods we recorded the realized performance after applying the open-loop inputs and converted it to a percentage error with respect to the best possible performance (i.e., if the deterministic model was exactly known). For each noise-to-signal ratio, 100 simulations were conducted with different random data sets. The results are displayed in the box plot in Figure <ref> and show that both methods perform well for low levels of noise (up to approximately $2\%$ noise-to-signal ratio). As the data becomes noisier, the performance of the direct method degrades significantly, while the performance of the indirect method remains relatively constant.  We remark that a slightly better albeit qualitatively similar result is obtained with the regularizer  $ \|(I-\Pi)g\|_{2}^{2}$.

We attribute these observations to the fact that identification de-noises the data. These results confirm our hypothesis that the indirect method is superior in terms of “variance” error.

Comparison of direct and indirect methods for varying noise.

§.§ Comparison: Deterministic Nonlinear System

We now consider the scenario where the direct and indirect methods are subject to a “bias” error, but not a “variance” error. Consider the discrete-time nonlinear Lotka-Volterra dynamics  considered for direct data-driven control in [57]

\[
\begin{aligned}
&= f_{\textup{nonlinear}}(x(t_k),u(t_k))
\\
x_1(t_k) + \Delta t(ax_1(t_k) - bx_1(t_k)x_2(t_k)) \\
x_2(t_k) + \Delta t(dx_1(t_k)x_2(t_k) -cx_2(t_k) + u(t_k))
\end{smallmatrix}\right]\,,
\end{aligned}
\]

where $t_{k+1} - t_k = \Delta t = 0.01$, $a=c=0.5$, $b=0.025$, $d=0.005$,  and $x(t_k) = \begin{bmatrix} x_1(t_k) & x_2(t_k)\end{bmatrix}^{\top}$. Here, $(x_1(t_k),x_2(t_k))$ denote prey and predator populations, and $u(t_k)$ is the input. A linearization about the equilibrium  $(\bar{u},\bar{x}_1,\bar{x}_2)=(0,c/d,a/b)$
yields the affine linear system

\[
\begin{aligned}
= f_{\textup{linear}}(x(t_k),u(t_k),\bar{x}_1,\bar{x}_2)
\\
&= \left[\begin{smallmatrix}
x_1(t_k) + \Delta t\left((a-b\bar{x}_2)(x_1(t_k)-\bar{x}_1) - b\bar{x}_1(x_2(t_k)-\bar{x}_2)\right) \\
x_2(t_k) + \Delta t\left(d\bar{x}_2(x_1(t_k)-\bar{x}_1) +(d\bar{x}_1 -c)(x_2(t_k)-\bar{x}_2) + u(t_k)\right)
\end{smallmatrix}\right]\,.
\end{aligned}
\]

We expect direct data-driven control (<ref>) to perform well on such a nonlinear system for two reasons: $(i)$ nonlinear systems can be well approximated by LTI systems of sufficiently high complexity; and $(ii)$ the direct method (<ref>) does not specify the LTI system complexity (e.g., by enforcing rank constraints).

We compare the direct and indirect methods for varying degree of nonlinearity by interpolating between $f_{\textup{nonlinear}}$ and $f_{\textup{linear}}$, i.e., we study the interpolated system

\begin{equation}
\begin{aligned}
x(t_{k+1}) &= \epsilon \cdot f_{\textup{linear}}(x(t_k),u(t_k),\bar{x}_1,\bar{x}_2)\\
&\quad + (1-\epsilon) \cdot f_{\textup{nonlinear}}(x(t_k),u(t_k))
\end{aligned}
\label{eq:interpolated_dynamics}
\end{equation}

for $\epsilon\in [0,1]$. For $\epsilon=1$ (resp. $\epsilon=0$), the dynamics are purely affine (resp. nonlinear). For each $\epsilon\in\{0,0.1,\dots,1\}$, $T=2415$ data points were collected by applying a noisy sinusoidal input $u(t_k) = 2(\sin(t_k)+\sin(0.1t_k)))^2 + v(t_k)$ with $v(t_k)$ sampled from a Gaussian random variable. Full state measurement was assumed. The data collection was repeated for 100 different initial conditions. For each degree of nonlinearity $\epsilon\in\{0,0.1,\dots,1\}$ and each initial condition, the data was used to compute optimal open-loop control inputs using direct and indirect methods. The control cost was chosen as $c_\textup{ctrl}({w}-w_r) = \| {w}-w_r\|_{2}^{2}$ with equilibrium reference $w_r = (0,100,20,\dots,0,100,20)$, $L=600$, and $w=(u,x)$.

For the indirect method, the inner system identification optimization problem given by (<ref>) is solved using the subspace approach N4SID [26] with initial condition horizon $\Tini=4$, and prediction horizon $L=600$. A model order of 4 was chosen, as it produced the best performance as measured by the realized control cost.
Optimal control inputs were then computed by solving (<ref>).
For comparison, we chose the direct method (<ref>) with $h(g) = \|g\|_{1}$, $\Tini=4$, and $\lambda=8000$.
The performance was measured with the realized control cost after applying the open-loop inputs to system (<ref>). As before, the hyper-parameters of both direct/indirect methods were judicially chosen and kept constant for all simulations.

Comparison of direct and indirect methods for varying nonlinearity.

The results displayed in Figure <ref> show that both methods perform well for low levels of nonlinearity: $\epsilon \in [0.7, 1]$. As the system becomes increasingly nonlinear, the performance of the indirect method degrades significantly, while the performance of the direct method remains relatively constant. We attribute this observation to the fact that the indirect method incurs a “bias” error from selecting a linear model class and applying certainty-equivalence control, while the direct method uses data from the nonlinear system without bias.

These findings confirm our earlier bias-variance observations from Figure <ref>.

§ DISCUSSION AND CONCLUSIONS

We studied the relationship between indirect and direct data-driven control formulated as bi-level (first-identify, then control) and single-level regularized (based on the Fundamental Lemma) optimization problems, respectively.  An intermediate multi-criteria problem allowed us to efficiently transition between both formulations. We concluded that the regularized direct approach can be viewed as a convex relaxation of the indirect approach, where the choice regularizer depended on the problem formulation and accounted for an implicit identification step. We also discovered a novel regularizer that is consistent and accounts for least-square identification.

Our results suggested the use of the indirect method in case of “variance” errors and the use of the direct method in presence of “bias” errors (e.g., a nonlinear system or when selecting a wrong model order). These insights echo the bias-variance trade-offs previously encountered for direct and indirect methods in [20, 51], and they
shed some partial light on the remarkable empirical performance of (direct) data-enabled predictive control applied to nonlinear systems.

As a  limitation, our results  concern only the open-loop predictive control problem, though we ultimately care about the realized performance, especially in a receding horizon closed-loop implementation. Some preliminary results on the realized performance of regularized control formulations were obtained in [46] through the lens of robust optimization, but the topic remains largely open. Moreover, we believe that the proposed multi-criteria data-driven control formulation is important in its own right and may deliver excellent performance if one were to find a convex formulation and appropriate trade-off parameter. Both of these are formidable tasks for future work.

Finally, we believe that our approach is also applicable to other identification and control formulations and may deliver interesting and novel direct data-driven control formulations.

§ ACKNOWLEDGEMENTS

The authors acknowledge their colleagues at ETH Zürich, in particular Miguel Picallo Cruz, for fruitful discussions.

Consider the mathematical program (MP)

x  ∈Cf(x)MP:

g(x) ≤0
h(x) = 0

where $C \subset \real^{n}$ is closed, and $f,g,h$ are lower semicontinuous maps from $\real^{n}$ to $\real$, $\real^{m}$, and $\real$. Consider the perturbation

x  ∈Cf(x)MP_ϵ:

g(x) ≤0
h(x) = ϵ

for $\epsilon \in \real$. We recall the definition of partial calmness [53]:

Let $x^{\star}$ solve $\boldsymbol{MP}$, and let $\mathbb B_{n}$ denote the open unit ball in $\mathbb R^{n}$. Then $\boldsymbol{MP}$ is said to be partially calm  at $x^{\star}$ provided that there are $\mu>0$ and $\delta>0$ such that, for all $\epsilon \in \delta \mathbb B_{1}$ and all $x \in x^{\star}+\delta \mathbb B_{n}$ feasible for $\boldsymbol{MP}_{\epsilon}$, one has

\begin{equation}
f(x) + \mu | h(x) | \geq f(x^{\star})
\label{eq: calmness}
\,.
\end{equation}

Partial calmness is equivalent to exact penalization. In particular, consider for $\mu \geq 0$ the penalized mathematical program

x  ∈Cf(x)  +  μ·|h(x)|PMP_μ:

g(x) ≤0

We summarize <cit.>,<cit.>:

Assume that $f$ is continuous, and let $x^{\star}$ be a local minimizer of $\boldsymbol{MP}$.
Then $\boldsymbol{MP}$ is partially calm at $x^{\star}$ if and only if there is $\mu^{\star} >0$ so that $x^{\star}$ is a local minimizer of $\boldsymbol{PMP}_{\mu}$ for all $\mu \geq \mu^{\star}$.

Moreover, any local minima of $\boldsymbol{PMP}_{\mu}$ with $\mu > \mu^{\star}$ are also local minima of $\boldsymbol{MP}$.

Partial calmness has been studied for a range of problems, particularly bi-level problems [52, 53, 58]. Of specific importance to us is a related result due to Clarke <cit.> which allows for exact penalization (or equivalently partial calmness) and reads in our notation as follows.

Consider the mathematical program $\boldsymbol{MP}$ and its penalized version $\boldsymbol{PMP}_{\mu}$. Assume that $f$ is Lipschitz continuous with Lipschitz constant $L$, the equality constraint takes the form of a distance to a closed set $S \subset C$,

\begin{equation*}
0 = h(x) = \text{distance}(x,S) = \text{inf}_{y \in S} \|x-y\| \,,
\end{equation*}

and $\boldsymbol{MP}$ attains a minimum. Then for $\mu > \mu^{\star} = L$, any local minimum of $\boldsymbol{PMP}_{\mu}$ is also a local minimum of $\boldsymbol{MP}$.

We note that $\|\cdot\|$ in Proposition <ref> can be an arbitrary norm. For a more general problem setup with a parametric set $S$ depicting the value function of an (inner) optimization problem, the reader is referred to <cit.>. Additionally, the setup can be extended to (squared) merit functions as penalty functions <cit.>. These generalize the notion of distance but are easier to formulate and compute.

[1]
L. Hewing, K. P. Wabersich, M. Menner, and M. N. Zeilinger, “Learning-based
model predictive control: Toward safe learning in control,” Annual
Review of Control, Robotics, and Autonomous Systems, vol. 3, pp. 269–296,

[2]
G. Pillonetto, F. Dinuzzo, T. Chen, G. De Nicolao, and L. Ljung, “Kernel
methods in system identification, machine learning and function estimation: A
survey,” Automatica, vol. 50, no. 3, pp. 657–682, 2014.

[3]
A. Chiuso and G. Pillonetto, “System identification: A machine learning
perspective,” Annual Review of Control, Robotics, and Autonomous
Systems, vol. 2, pp. 281–304, 2019.

[4]
Z.-S. Hou and Z. Wang, “From model-based control to data-driven control:
Survey, classification and perspective,” Information Sciences, vol.
235, pp. 3–35, 2013.

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

[6]
I. Markovsky and F. Dörfler, “Behavioral systems theory in data-driven
analysis, signal processing, and control,” July 2021. [Online]. Available:
<http://homepages.vub.ac.be/ imarkovs/publications/overview-ddctr.pdf>

[7]
H. Hjalmarsson, “From experiment design to closed-loop control,”
Automatica, vol. 41, no. 3, pp. 393–438, 2005.

[8]
S. Meyn, Control Systems and Reinforcement Learning.1em plus
0.5em minus 0.4emCambridge University Press, 2022.

[9]
H. Hjalmarsson, M. Gevers, S. Gunnarsson, and O. Lequin, “Iterative feedback
tuning: Theory and applications,” IEEE Control Systems Magazine,
vol. 18, no. 4, pp. 26–41, 1998.

[10]
M. C. Campi, A. Lecchini, and S. M. Savaresi, “Virtual reference feedback
tuning: A direct method for the design of feedback controllers,”
Automatica, vol. 38, no. 8, pp. 1337–1346, 2002.

[11]
A. S. Bazanella, L. Campestrini, and D. Eckhard, Data-driven controller
design: the H2 approach.1em plus 0.5em minus 0.4emSpringer
Science & Business Media, 2011.

[12]
H. Hjalmarsson, M. Gevers, and F. De Bruyne, “For model-based control design,
closed-loop identification gives better performance,” Automatica,
vol. 32, no. 12, pp. 1659–1673, 1996.

[13]
M. Gevers, “Identification for control: From the early achievements to the
revival of experiment design,” European Journal of Control, vol. 11,
pp. 1–18, 2005.

[14]
R. J. Schrama, “Accurate identification for control: The necessity of an
iterative scheme,” IEEE Transactions on Automatic Control, vol. 37,
no. 7, pp. 991–994, 1992.

[15]
S. Formentin and A. Chiuso, “CoRe: control-oriented regularization for
system identification,” in 2018 IEEE Conference on Decision and
Control (CDC).1em plus 0.5em minus 0.4emIEEE, 2018, pp.

[16]
A. Feldbaum, “Dual control theory problems,” IFAC Proceedings Volumes,
vol. 1, no. 2, pp. 541–550, 1963.

[17]
M. Ferizbegovic, J. Umenberger, H. Hjalmarsson, and T. B. Schön, “Learning
robust lq-controllers using application oriented exploration,” IEEE
Control Systems Letters, vol. 4, no. 1, pp. 19–24, 2019.

[18]
C. A. Larsson, A. Ebadat, C. R. Rojas, X. Bombois, and H. Hjalmarsson, “An
application-oriented approach to dual control with excitation for closed-loop
identification,” European Journal of Control, vol. 29, pp. 1–16,

[19]
A. Iannelli, M. Khosravi, and R. S. Smith, “Structured exploration in the
finite horizon linear quadratic dual control problem,”
IFAC-PapersOnLine, vol. 53, no. 2, pp. 959–964, 2020.

[20]
L. Campestrini, D. Eckhard, A. S. Bazanella, and M. Gevers, “Data-driven model
reference control design by prediction error identification,” Journal
of the Franklin Institute, vol. 354, no. 6, pp. 2628–2647, 2017.

[21]
J. C. Willems, “The behavioral approach to open and interconnected systems,”
IEEE Control Systems Magazine, vol. 27, no. 6, pp. 46–99, 2007.

[22]
——, “Paradigms and puzzles in the theory of dynamical systems,”
IEEE Transactions on Automatic Control, vol. 36, no. 3, pp. 259–294,

[23]
J. C. Willems and J. W. Polderman, Introduction to mathematical systems
theory: A behavioral approach.1em plus 0.5em minus 0.4emSpringer, 1997, vol. 26.

[24]
P. Van Overschee and B. De Moor, Subspace identification for linear
systems: Theory, Implementation, Applications.1em plus 0.5em minus
0.4emSpringer Science & Business Media, 2012.

[25]
T. Katayama, Subspace methods for system identification.1em plus
0.5em minus 0.4emSpringer Science & Business Media, 2006.

[26]
P. Van Overschee and B. De Moor, “N4SID: Subspace algorithms for the
identification of combined deterministic-stochastic systems,”
Automatica, vol. 30, no. 1, pp. 75–93, 1994.

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

[28]
H. J. van Waarde, C. De Persis, M. K. Camlibel, and P. Tesi, “Willems'
fundamental lemma for state-space systems and its extension to multiple
datasets,” IEEE Control Systems Letters, vol. 4, no. 3, pp. 602–607,

[29]
I. Markovsky and F. Dörfler, “Identifiability in the behavioral setting,”
September 2020, Submitted. Available at
http://homepages.vub.ac.be/ imarkovs/publications/identifiability.pdf.

[30]
I. Markovsky, J. C. Willems, S. Van Huffel, and B. De Moor, Exact and
approximate modeling of linear systems: A behavioral approach.1em
plus 0.5em minus 0.4emSIAM, 2006, vol. 11.

[31]
I. Markovsky, J. C. Willems, P. Rapisarda, and B. L. De Moor, “Algorithms for
deterministic balanced subspace identification,” Automatica, vol. 41,
no. 5, pp. 755–766, 2005.

[32]
W. Favoreel, B. De Moor, and M. Gevers, “SPC: subspace predictive control,”
IFAC Proceedings Volumes, vol. 32, no. 2, pp. 4004–4009, 1999.

[33]
S. J. Qin, W. Lin, and L. Ljung, “A novel subspace identification approach
with enforced causal models,” Automatica, vol. 41, no. 12, pp.
2043–2053, 2005.

[34]
B. Huang and R. Kadali, Dynamic modeling, predictive control and
performance monitoring: A data-driven subspace approach.1em plus
0.5em minus 0.4emSpringer, 2008.

[35]
J. Berberich, C. W. Scherer, and F. Allgöwer, “Combining prior knowledge
and data for robust controller design,” arXiv preprint
arXiv:2009.05253, 2020.

[36]
H. J. van Waarde, M. K. Camlibel, and M. Mesbahi, “From noisy data to
feedback controllers: non-conservative design via a matrix S-lemma,”
arXiv preprint arXiv:2006.00870, 2020.

[37]
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, 2019.

[38]
I. Markovsky and P. Rapisarda, “Data-driven simulation and control,”
International Journal of Control, vol. 81, no. 12, pp. 1946–1959,

[39]
I. Markovsky, “A missing data approach to data-driven filtering and control,”
IEEE Transactions on Automatic Control, vol. 62, no. 4, pp.
1972–1978, 2016.

[40]
J. Berberich, J. Köhler, M. A. Muller, and F. Allgöwer, “Data-driven
model predictive control with stability and robustness guarantees,”
IEEE Transactions on Automatic Control, 2020.

[41]
J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control:
In the shallows of the DeePC,” in European Control Conference,
2019, pp. 307–312.

[42]
——, “Regularized and distributionally robust data-enabled predictive
control,” in IEEE Conference on Decision and Control, 2019, pp.

[43]
——, “Distributionally robust chance constrained data-enabled predictive
control,” 2020, In press. DOI 10.1109/TAC.2021.3097706.

[44]
A. Xue and N. Matni, “Data-driven system level synthesis,” arXiv
preprint arXiv:2011.10674, 2020.

[45]
L. Huang, J. Zhen, J. Lygeros, and F. Dörfler, “Quadratic regularization
of data-enabled predictive control: Theory and application to power converter
experiments,” 2020, In press. Available at

[46]
L. Huang, Z. Jianzhe, J. Lygeros, and F. Dörfler, “Robust data-enabled
predictive control: Tractable formulations and performance guarantees,”
2021, Submitted. Available at <https://arxiv.org/abs/2105.07199>.

[47]
L. Huang, J. Coulson, J. Lygeros, and F. Dörfler, “Decentralized
data-enabled predictive control for power system oscillation damping,”
IEEE Transactions on Control Systems Technology, 2021, In press. DOI

[48]
P. G. Carlet, A. Favato, S. Bolognani, and F. Dörfler, “Data-driven
predictive current control for synchronous motor drives,” in 2020 IEEE
Energy Conversion Congress and Exposition (ECCE), 2020, pp. 5148–5154.

[49]
E. Elokda, J. Coulson, P. Beuchat, J. Lygeros, and F. Dörfler,
“Data-enabled predictive control for quadcopters,” International
Journal of Robust and Nonlinear Control, 2019, In press. DOI

[50]
M. Yin, A. Iannelli, and R. S. Smith, “Maximum likelihood estimation in
data-driven modeling and control,” arXiv preprint arXiv:2011.00925,

[51]
V. Krishnan and F. Pasqualetti, “On direct vs indirect data-driven predictive
control,” arXiv preprint arXiv:2103.14936, 2021.

[52]
J. Ye and D. Zhu, “Optimality conditions for bilevel programming problems,”
Optimization, vol. 33, no. 1, pp. 9–27, 1995.

[53]
J. Ye, D. Zhu, and Q. J. Zhu, “Exact penalization and necessary optimality
conditions for generalized bilevel programming problems,” SIAM Journal
on optimization, vol. 7, no. 2, pp. 481–507, 1997.

[54]
H. Xu, C. Caramanis, and S. Mannor, “Robust regression and lasso,” IEEE
Transactions on Information Theory, vol. 56, no. 7, pp. 3561–3574, 2010.

[55]
T. Hastie, R. Tibshirani, and M. Wainwright, Statistical learning with
sparsity: the lasso and generalizations.1em plus 0.5em minus
0.4emCRC press, 2015.

[56]
I. Landau, D. Rey, A. Karimi, A. Voda, and A. Franco, “A flexible transmission
system as a benchmark for robust digital control,” European Journal of
Control, vol. 1, no. 2, pp. 77–96, 1995.

[57]
E. Kaiser, J. N. Kutz, and S. L. Brunton, “Sparse identification of nonlinear
dynamics for model predictive control in the low-data limit,”
Proceedings of the Royal Society A, vol. 474, no. 2219, p. 20180335,

[58]
P. Mehlitz, L. I. Minchenko, and A. B. Zemkoho, “A note on partial calmness
for bilevel optimization problems with linearly structured lower level,”
Optimization Letters, pp. 1–15, 2020.

[59]
F. H. Clarke, Optimization and nonsmooth analysis.1em plus 0.5em
minus 0.4emSIAM, 1990.

[]Florian Dörfler (S'09-M'13-S'21) is an Associate Professor at the Automatic Control Laboratory at ETH Zürich and the Associate Head of the Department of Information Technology and Electrical Engineering. He received his Ph.D. degree in Mechanical Engineering from the University of California at Santa Barbara in 2013, and a Diplom degree in Engineering Cybernetics from the University of Stuttgart in 2008. From 2013 to 2014 he was an Assistant Professor at the University of California Los Angeles. His primary research interests are centered around control, optimization, and system theory with applications in network systems, in particular electric power grids. He is a recipient of the distinguished young research awards by IFAC (Manfred Thoma Medal 2020) and EUCA (European Control Award 2020). His students were winners or finalists for Best Student Paper awards at the European Control Conference (2013, 2019), the American Control Conference (2016), the Conference on Decision and Control (2020), the PES General Meeting (2020), and the PES PowerTech Conference (2017). He is furthermore a recipient of the 2010 ACC Student Best Paper Award, the 2011 O. Hugo Schuck Best Paper Award, the 2012-2014 Automatica Best Paper Award, the 2016 IEEE Circuits and Systems Guillemin-Cauer Best Paper Award, and the 2015 UCSB ME Best PhD award.

[]Jeremy Coulson (S'09-M'13)  is a PhD student with the Automatic Control Laboratory at ETH Zürich. He received his Master of Applied Science in Mathematics & Engineering from Queen's University, Canada in August 2017. He received his B.Sc.Eng degree in Mechanical Engineering & Applied Mathematics from Queen's University in 2015. His research interests include data-driven control methods, and stochastic optimization.

[]Ivan Markovsky
is an Associate Professor at the department ELEC of the Vrije Universiteit Brussel. He received his Ph.D. degree in Electrical Engineering from the Katholieke Universiteit Leuven in February 2005. From 2006 to 2012 he was an Assistant Professor at the School of Electronics and Computer Science of the University of Southampton. He is a recipient of an ERC starting grant "Structured low-rank approximation: Theory, algorithms, and applications" 2010–2015, Householder Prize honorable mention 2008, and research mandate by the Vrije Universiteit Brussel research council 2012–2022. His main research interests are computational methods for system theory, identification, and data-driven control in the behavioral setting.