Data-conforming data-driven control: avoiding premature generalizations beyond data

Abstract

Data-driven and adaptive control approaches face the problem of introducing sudden distributional shifts beyond the distribution of data encountered during learning. Therefore, they are prone to invalidating the very assumptions used in their own construction. This is due to the linearity of the underlying system, inherently assumed and formulated in most data-driven control approaches, which may falsely generalize the behavior of the system beyond the behavior experienced in the data. This paper seeks to mitigate these problems by enforcing consistency of the newly designed closed-loop systems with data and slow down any distributional shifts in the joint state-input space. This is achieved through incorporating affine regularization terms and linear matrix inequality constraints to data-driven approaches, resulting in convex semi-definite programs that can be efficiently solved by standard software packages. We discuss the optimality conditions of these programs and then conclude the paper with a numerical example that further highlights the problem of premature generalization beyond data and shows the effectiveness of our proposed approaches in enhancing the safety of data-driven control methods.

Data-driven control, adaptive control, system identification, robust control, offline reinforcement learning.

1 Introduction

The development of adaptive control systems, while marked by significant theoretical advancements, has historically faced skepticism from practitioners, rooted in concerns over the reliability and robustness of such systems in real-world applications. Brian Anderson, in his article [4], explains reasons for this distrust and the inherent dangers of adaptive control algorithms. In this paper we highlight another fundamental vulnerability in adaptive learning and data-driven control algorithms: premature and often false generalization beyond the seen data. That is, hastily generalizing the behavior of the system beyond its behavior that was observed in the data. This is pervasive in modern data-driven control methods and can lead to catastrophic results in real-world applications. We propose practical methods for overcoming this vulnerability in a computationally efficient manner and using standard off-the-shelf software packages.

Data has become a cornerstone across every scientific domain, and it is foundational for the advancements in artificial intelligence [23]. In control theory, the reliance on data in control design is not new; and the subfields of system identification [26], robust control [17], and adaptive control [6], pioneered decades ago, have long been central to the field of control. The historical insights these subfields have generated are still crucial to the development of modern data-driven control and reinforcement learning (RL) algorithms. Take, for example, the “exciting” input for the stability of adaptive control [5, 27], which is analogous to the role of stochastic policies in RL algorithms [35], or the problem of the chaos phenomenon resulting from learning and control loops operating in comparable time scales, discovered in adaptive control [28] and later in RL [39]. These insights have also been crucial to the development of dual control, and its automatic experiment design aspect [32, 21], which have their connections to the exploration vs exploitation trade-off in RL [36].

The aforementioned insights, particularly the difficulties they are coupled with, hindered the spread of adaptive control algorithms in real-world applications, and instead, practitioners and control theorists alike found refuge in the more established and guaranteed robust control algorithms [34]. The importance of adaptive control cannot be excluded, however; and the recent surge in the application of data-driven control and RL algorithms (in power systems, for instance [18, 41] ) and RL’s role in the development of large language models [7] revive this importance. Thus, we are motivated to further understand the problems of adaptive control and expand their remedies.

In this paper we highlight an extra vulnerability of adaptive and data-driven control, not addressed explicitly by Anderson in [4]. This vulnerability is the premature and possibly false generalizations beyond data, resulting from extrapolating the behavior of the system beyond what was experienced in the data. This in turn can lead to compromising safety and performance under the existence of unmodeled nonlinearities in the true underlying system.

Modern data-driven control methods, whether direct (model-free) [13] or indirect (model-based) [40], mostly assume the linearity of the data-generating dynamics (although a few approaches extend direct methods to very limited forms of nonlinearities [15]). The problem with the linearity assumption is that it imposes the universality of the data, in the sense that the underlying dynamics must “behave similarly” (according to the same linear dynamics), beyond data and in any region of the state space. This assumption is often invalid, however, because various systems admit nonlinear effects beyond their engineered operating conditions and typical control design does not account for this situation.

We emphasize that in a system identification and data collection experiment, an identified system (or recorded data in the direct approaches case) is not necessarily valid under different experimental and data collection conditions. For example, a controller designed to be stabilizing/optimal with respect to the identified model (or recorded data) is guaranteed to be stabilizing/optimal only with respect to this specific identified model. However, eventually this controller will be connected to the actual data-generating system. This connection likely will change the operating conditions of the system and therefore may invalidate the exact assumptions this controller was based on. This can possibly activate unmodeled nonlinearities that can result in instability or severely degraded performance. Therefore, even with batch learning and iterative identification in a much smaller time scale than control [2], a change in the controller may still result in a rapid shift in the system’s operating conditions, compromising safety and performance.

The methods we propose in this paper enforce the consistency between the designed closed-loop system and the data encountered during learning. We achieve this consistency by augmenting the standard linear quadratic regulator (LQR) problem with affine regularization terms and linear matrix inequalities (LMIs) in the corresponding decision variables. The result is an affine semi-definite program (SDP) with LMI constraints that can be solved, even for high-dimensional systems, efficiently and using off-the-shelf software packages. Furthermore, we introduce a single hyperparameter that enforces the consistency level; in other words, this hyperparameter is an exploration vs exploitation balance factor. Our method allows for an iterative identification (data collection) and control that prevents rapid distributional shifts and the related sudden violations of data consistency.

In mathematical terms, this augmentation of the standard LQR problem includes working with the parametrization given by the controllability-type Gramian [25], which is also the steady-state covariance matrix of the system state. This makes it possible to relate the closed-loop state-input distribution to that of the data, using LMI constraints and/or regularization terms such as the Frobenius norm on covariance matrices or the Kullback–Liebler divergence between distributions. Throughout the paper we call a control design approach data-conforming if it enforces the consistency with data according to our aforementioned methods.

Inspired by the problems addressed in [4], our problem is closest in spirit to the so-called Rohrs’ counterexample [33], which sought to invalidate adaptive control algorithms in the 1980s. Part of Rohrs’ argument is that every physical system has unmodeled (parasitic) high-frequency dynamics. Because contemporary adaptive control led to highly nonlinear loops even for simple linear systems, these loops could internally generate high-frequency signals, exciting the unmodeled dynamics and possibly leading to instability. Our concern regarding modern data-driven control is analogous to that of Rohrs’ with adaptive control, but we target unmodeled nonlinear dynamics instead of unmodeled high-frequency ones. Through iterative identification and carefully expanding the closed-loop bandwidth, one can avoid exciting high-frequency dynamics [3]. Analogously, through iterative identification and careful data-conforming control design, one can avoid sudden exposures to new areas of the state space, which potentially contain harmful, unmodeled nonlinearities.

Our data-conforming control design approach also bears resemblance to the unfalsified control approach [12] in enforcing consistency with past data. Contrasting with this approach, which categorizes consistency through Boolean values, our approach is built around consistency in distances between distributions, allowing for generalized formulations that can cope with stochastic dynamics. Moreover, the method of [12] relies on model reference adaptive control, whereas our approach augments modern optimal control design via extending the standard LQR problem with affine regularization terms and LMI constraints. Our approach readily integrates with many modern control design formulations, allows for multivariate state-space systems, and can be extended to handle receding-horizon and input-output formulations.

The idea of consistency with data is also fundamental in the field of offline RL [1], but there the main motivation is data efficiency by initially limiting exploration. Offline RL algorithms are generally complex [20] and rely on stochastic nonlinear programming methods. Instead, our approach insists on computationally efficient methods that blend with modern control design techniques.

One can argue that modern data-driven policy gradient methods [42, 19] with small enough step sizes can prevent sudden distributional shifts and therefore enhance consistency with data. This is indeed a plausible argument. However, it is difficult to relate between the control gain step size and the corresponding distributional shift beyond a simple ($\epsilon-\delta$) continuity argument. Instead, our approach deals explicitly with distributions and directly dampens their shifts. Moreover, a very dampened policy gradient procedure may be slow and hence inadequate for many time-varying dynamics, limiting its applicability to constant or pseudo-constant dynamics. On the other hand, our approach decouples distributional shifts from learning system variations and therefore can accomplish both tasks effectively and simultaneously.

We conclude the paper with a simple yet telling numerical example that explains the premature generalization problem inherent in modern data-driven approaches and shows how our proposed solutions can mitigate this problem.

2 Problem Formulation

Consider the dynamic system

$$x_{k+1} =f(x_{k},u_{k},w_{k}),$$ (1)

where $x_{k}\in\mathbb{R}^{r_{x}}$ is the state and $u_{k}\in\mathbb{R}^{r_{u}}$ is the control input. The function $f$ is bounded but unknown and possibly nonlinear. The exogenous disturbance $w_{k}\in\mathbb{R}^{r_{x}}$ is independent and identically distributed. It is also independent from $x_{0}$, the initial condition, which is also random and has finite mean and covariance.

The goal of this paper is to design a control law that minimizes the cost function

$$J=\lim_{k\to\infty}\mathbb{E\,}\left\{x_{k}^{\top}Qx_{k}+u_{k}^{\top}Ru_{k}\right\},$$ (2)

where $Q\succeq 0$ and $R\succ 0$ are positive semi-definite and positive definite, respectively, and the expectation averages over all the possible realizations of $x_{0}$ and $w_{k}$.

Since $f$ is unknown, we assume that we have access to it only through experiments. Hence, we rely on the data-driven paradigm for the control design. During the first data collection experiment, we assume an initial control law of the form

$$u_{k}=K_{0}x_{k}+v_{k},$$ (3)

where this law is locally stabilizing (not necessarily optimal) and $v_{k}$ is a persistently exciting (PE) signal, in the following sense.

Assumption 1.

(The PE assumption [14, condition (6)]). The signal $v_{k}$ in (3) is PE for (1). That is, given a natural number $N\geq(r_{u}+1)r_{x}+r_{u}$, every control realization $\{u_{k}\}_{k=0}^{N}$ and the corresponding state realization $\{x_{k}\}_{k=0}^{N}$ produce the matrix

$${D}:=\begin{bmatrix}X\\ U\end{bmatrix}=\begin{bmatrix}x_{0}&\ldots&x_{N}\\ u_{0}&\ldots&u_{N}\end{bmatrix}$$ (4)

with full row rank $r_{x}+r_{u}$.

The PE assumption is typical in system identification and data-driven control design and equates to an effective experiment design and a minimal descriptive state definition.

Problem statement

Starting from the initial control law of the form (3) and given access to the resulting input and state data (4), with the possibility of iteratively changing the control law and collecting new data, design a linear state feedback gain $K$ that minimizes the cost $J$ in (2).

Assumption 2.

The data (4) is centered around zero, and the state-input joint distribution can be approximated by a multivariate Gaussian with a positive definite covariance.

This assumption is true if the underlying system is linear Gaussian. For nonlinear systems, however, this point can be achieved, for instance, if the data was collected locally within a specified region of the state-input space under some operational limits. In other words, this assumption is satisfied if there is a “good enough” linear approximation of the underlying system within each operating conditions (4).

One may ask about the purpose of considering a general nonlinear system in (1) if the data can be predicted by an approximate linear model. The reason is to remember that once a new controller is designed, the resulting closed-loop system may operate in regions of the state-input space not adequately explored by the learning data and over which this approximate linear model is no longer a valid approximation. This point can be easily forgotten if we start with a linear system in (1).

3 Background

In this section we start by presenting the standard LQR problem. Although the problem formulation targets the nonlinear dynamics (1), our derivations in the subsequent sections assume the existence of a valid linear approximation under each operating condition.

3.1 Standard LQR

Suppose, in this subsection only, that the state-space system (1) is linear,

$$x_{k+1}=f(x_{k},u_{k},w_{k})=Ax_{k}+Bu_{k}+w_{k},$$

where $x_{0},w_{k}$ are random variables and their statistics are as in (1), the disturbance $w_{k}$ has zero mean and a known (in this subsection only, it will be approximated from data in the later subsections) covariance $W\succeq 0$, and the pair $(A,B)$ is known and controllable [25]. The LQR problem is obtaining a controller of the form $u_{k}=Kx_{k}+v_{k}$, where $v_{k}$ is a white noise of zero mean and a fixed user-defined covariance $V\succ 0$, to minimize (2). The signal $v_{k}$ is to preserve Assumption 1 for future or iterative control design.

The cost function (2), up to an additive constant in $K$, can be rewritten as (for a detailed derivation, check Appendix A.1)

$$J=\textit{tr}\left(P\left[W+BVB^{\top}\right]\right),$$ (5)

where $P$ is the observability-type Gramian [25], given by

$$P=\sum_{k=0}^{\infty}[A+BK]^{k\,\top}\left[Q+K^{\top}RK\right][A+BK]^{k},$$ (6)

where $Q\succeq 0$ and $R\succ 0$ are the user-defined weighting matrices, as in (2), and $(Q,A)$ is assumed detectable. The optimal control gain $K$ (which minimizes $J$) is then given by the solution to the following problem.

Problem 1.

Standard LQR (Observability-type Gramian) [8, Sec. 4.1]:

$$\text{Evaluate }K=-\left[R+B^{\top}PB\right]^{-1}B^{\top}PA,$$

when $P$ satisfies the algebraic Riccati equation

$$0=Q+A^{\top}PA-P-A^{\top}PB\left[R+B^{\top}PB\right]^{-1}B^{\top}PA.$$

The problem above results from the square completion in $K$ of the Lyapunov equation

$$0=[A+BK]^{\top}P[A+BK]-P+Q+K^{\top}RK,$$ (7)

then chosing $K=-\left[R+B^{\top}PB\right]^{-1}B^{\top}PA$, which corresponds to the minimum $P$ (in the Loewner ($\succeq$) ordering sense).

There is an equivalent characterization to solving the LQR problem. Using the linearity and the cyclic property of the trace, the cost can be written as (check Appendix A.1 for more details)

$$J =\textit{tr}\left(\left[Q+K^{\top}RK\right]\Sigma\right),$$ (8)

$$=\textit{tr}\left(Q\Sigma\right)+\textit{tr}\left(RK\Sigma K^{\top}\right),$$

where $\Sigma$ is the controllability-type Gramian, given by

$$\Sigma =\lim_{k\to\infty}\mathbb{E\,}\left\{x_{k}x_{k}^{\top}\right\},$$ (9)

$$=\sum_{k=0}^{\infty}\left[A+BK\right]^{k}\left[W+BVB^{\top}\right]\left[A+BK\right]^{k\,\top},$$

and is also the solution of the Lyapunov equation

$$\Sigma=\left[A+BK\right]\Sigma\left[A+BK\right]^{\top}+W+BVB^{\top}.$$ (10)

Since $V\succ 0$ and the pair $(A,B)$ is controllable, the controllability Gramian satisfies $\Sigma\succ 0$.

Using this parametrization, we first define the extra variable $Z_{0}$, such that $Z_{0}\succeq K\Sigma K^{\top}$, which is equivalent to (the Schur complement corresponding to) the LMI

$$\begin{bmatrix}Z_{0}&K\Sigma\\ \Sigma K^{\top}&\Sigma\end{bmatrix}\succeq 0,\text{ or, }\begin{bmatrix}Z_{0}&L\\ L^{\top}&\Sigma\end{bmatrix}\succeq 0,$$

using the change of variables $L=K\Sigma$. The stability-imposing inequality relaxation of (10) (we discuss the consequences of this relaxation in Section 5) can also be described as an LMI, that is,

$$\Sigma\succeq\left[A+BK\right]\Sigma\left[A+BK\right]^{\top}+W+BVB^{\top}$$

$$\iff\Sigma\succeq\left[A\Sigma+BL\right]\Sigma^{-1}\left[A\Sigma+BL\right]^{\top}+W+BVB^{\top}$$

$$\iff\begin{bmatrix}\Sigma-W-BVB^{\top}&A\Sigma+BL\\ \Sigma A^{\top}+L^{\top}B^{\top}&\Sigma\end{bmatrix}\succeq 0.$$

Therefore, the LQR problem can be described as an SDP with affine cost and LMI constraints [9].

Problem 2.

Standard LQR (Controllability-type Gramian):

$$\min_{\Sigma,\,L,\,Z_{0}}\textit{tr}\left(Q\Sigma\right)+\textit{tr}\left(RZ_{0}\right)$$

$$\text{s.t. }\Sigma\succ 0,\quad\begin{bmatrix}Z_{0}&L\\ L^{\top}&\Sigma\end{bmatrix}\succeq 0,$$

$$\begin{bmatrix}\Sigma-W-BVB^{\top}&A\Sigma+BL\\ \Sigma A^{\top}+L^{\top}B^{\top}&\Sigma\end{bmatrix}\succeq 0,$$

where the change of variables $L=K\Sigma$ has been used and the optimal control is then recovered from the optimal values $\Sigma_{\star}$ and $L_{\star}$ through $K_{\star}=L_{\star}\Sigma_{\star}^{-1}$.

The choice of using the controllability-type parametrization of the cost (8) is motivated by (9); the matrix $\Sigma$ is the steady-state state covariance matrix and thus has a statistical interpretation. Since this paper is in the context of data-driven control, the matrix $\Sigma$ is of great importance, and having it appear explicitly is a key. We show this in Section 4.

3.2 Data-driven certainty equivalence

In the case of unknown system matrices $A$ and $B$, conventional (indirect) data-driven approaches identify a model $(\widehat{A},\widehat{B})$ and $\widehat{W}$, then solve the standard LQR problem with the identified model assumed to represent the true system exactly, a procedure termed certainty equivalence LQR.

The estimate model $(\widehat{A},\widehat{B})$ and $\widehat{W}$ can be obtained as a solution to the least-squares problem [42]

$$\begin{aligned} &=\operatorname*{arg\,min}_{A,B}\left\lVert X_{err}\right\rVert_{F},\\ \widehat{W}&=\frac{1}{N}X_{err}X_{err}^{\top},\\ X_{err}&=X_{1}-[A,B]\begin{bmatrix}X_{0}\\ U_{0}\end{bmatrix},\end{aligned}$$ (11)

where $\lVert\cdot\rVert_{F}$ is the Frobenius norm. Using the acquired data (4), we let $X_{0}:=[x_{0},\ldots,x_{N-1}]$, $U_{0}:=[u_{0},\ldots,u_{N-1}]$, and $X_{1}:=[x_{1},\ldots,x_{N}]$.

The certainty equivalence LQR problem can be decomposed into an identification problem and then a control design problem. We describe this two-level approach as the following problem.

Problem 3.

Certainty equivalence LQR:

$$\min_{\Sigma,\,L,\,Z_{0}}\textit{tr}\left(Q\Sigma\right)+\textit{tr}\left(RZ_{0}\right)$$

$$\text{s.t. }\Sigma\succ 0,\quad\begin{bmatrix}Z_{0}&L\\ L^{\top}&\Sigma\end{bmatrix}\succeq 0,$$

$$\begin{bmatrix}\Sigma-\widehat{W}-\widehat{B}V\widehat{B}^{\top}&\widehat{A}\Sigma+\widehat{B}L\\ \Sigma\widehat{A}^{\top}+L^{\top}\widehat{B}^{\top}&\Sigma\end{bmatrix}\succeq 0,$$

$$\widehat{A},\widehat{B},\widehat{W}\text{ as in \eqref{eq:linear state-space ID}}.$$

The optimal control is recovered from the optimal values $\Sigma_{\star}$ and $L_{\star}$ through $K_{\star}=L_{\star}\Sigma_{\star}^{-1}$.

This data-driven approach is a basic certainty equivalence control design approach. For a discussion of the robustness and the sample efficiency, and to see how to enforce robustness under the existence of noise ($w_{k}\neq 0$) when the underlying model (1) is linear, an interested reader can consult [16] and [38]. In this paper, however, we tackle the problem of the false generalization of nonlinear systems by linear ones, beyond data, which is implicit in data-driven control algorithms.

Notice if the underlying system is indeed linear, the region of the state-input space from which the data (4) has been acquired is irrelevant, as long as the PE assumption, Assumption 1, is met. The reason is that if the system is linear, predicting its behavior in any region equates to predicting its behavior universally across the state-input space. This generally does not hold if the underlying system is nonlinear.

One shortcoming of using Problem 3 as the optimal control model is that once a linear model is identified and the resulting gain matrix $K$ is applied in feedback, the region in the state space, where the new closed-loop system is operating, may not be within or close to the region explored by the data (4) in the identification step. Therefore, if the underlying system is nonlinear, the new controller may force the system to go to regions of the domain of $f$ over which the system admits behaviors that are very different from what was experienced in the data. This implicit and inherent false universal generalization may in turn result in serious violations of safety or lead to a degraded performance.

4 Data-conforming data-driven control

In practice, for a controller to conform to the learning data, a measure of similarity between distributions has to be adopted and augmented to the cost/constraints of the LQR problem. For this purpose, several approaches and simplifications can be considered. We explore some of them in this section.

4.1 Conforming to the state data

Notice that Problem 3 is a convex optimization problem and that the decision variable $\Sigma$, the steady-state state covariance matrix, can be enforced to be similar to the state empirical covariance from data, via convex constraints or regularization.

Using (4), let $\Sigma_{data}$ be the state empirical covariance:

$$\Sigma_{data}=\frac{1}{N+1}{X}{X}^{\top}.$$

The data-conforming (in the state) version of Problem 3 can now be stated as follows.

Problem 4.

Data-conforming (in the state) LQR via a hard constraint:

$$\min_{\Sigma,\,L,\,Z_{0}}\textit{tr}\left(Q\Sigma\right)+\textit{tr}\left(RZ_{0}\right)$$

$$\text{s.t. }\Sigma\succ 0,\quad\begin{bmatrix}Z_{0}&L\\ L^{\top}&\Sigma\end{bmatrix}\succeq 0,\quad\Sigma=\Sigma_{data},$$

$$\begin{bmatrix}\Sigma-\widehat{W}-\widehat{B}V\widehat{B}^{\top}&\widehat{A}\Sigma+\widehat{B}L\\ \Sigma\widehat{A}^{\top}+L^{\top}\widehat{B}^{\top}&\Sigma\end{bmatrix}\succeq 0,$$

$$\widehat{A},\widehat{B},\widehat{W}\text{ as in \eqref{eq:linear state-space ID}},$$

and, similarly, the optimal control is recovered from the optimal values $\Sigma_{\star}$ and $L_{\star}$ through $K_{\star}=L_{\star}\Sigma_{\star}^{-1}$.

The new hard linear constraint $\Sigma=\Sigma_{data}$ has two shortcomings: (i) possible numerical instabilities and/or feasibility issues and (ii) limits to exploration beyond learning data. That is, potential future improvements via some exploration is not possible because the new controller will always seek to generate the same state distribution. This can be relaxed by using some margin, say, $\Sigma\succeq\Sigma_{data}-\epsilon\mathbb{I}$ and $\Sigma\preceq\Sigma_{data}+\epsilon\mathbb{I}$, for some $\epsilon>0$, predetermined or possibly included in the cost as a slack variable. This flexibility in the constraint, resembled in $\epsilon$, can determine the exploration vs exploitation balance of the new control design.

One can also design a convex regularization term using some norm on the covariance matrices. For example, the squared Frobenius norm $\lVert\Sigma-\Sigma_{data}\rVert_{F}^{2}=\textit{tr}\left(\left[\Sigma-\Sigma_{data}\right]\left[\Sigma-\Sigma_{data}\right]^{\top}\right)$ can be used by minimizing $\textit{tr}\left(Z^{\prime}\right)$, where $Z^{\prime}\succeq\left[\Sigma-\Sigma_{data}\right]\left[\Sigma-\Sigma_{data}\right]^{\top}$, or equivalently as an LMI [10],

$$\begin{bmatrix}Z^{\prime}&\Sigma-\Sigma_{data}\\ \Sigma-\Sigma_{data}&I\end{bmatrix}\succeq 0.$$

We use this regularization term in the following modified problem.

Problem 5.

Data-conforming (in the state) LQR via regularization:

$$\min_{\Sigma,\,L,\,Z_{0},\,Z^{\prime}}\textit{tr}\left(Q\Sigma\right)+\textit{tr}\left(RZ_{0}\right)+\gamma^{\prime}\textit{tr}\left(Z^{\prime}\right)$$

$$\text{s.t. }\Sigma\succ 0,\quad\begin{bmatrix}Z_{0}&L\\ L^{\top}&\Sigma\end{bmatrix}\succeq 0,$$

$$\begin{bmatrix}\Sigma-\widehat{W}-\widehat{B}V\widehat{B}^{\top}&\widehat{A}\Sigma+\widehat{B}L\\ \Sigma\widehat{A}^{\top}+L^{\top}\widehat{B}^{\top}&\Sigma\end{bmatrix}\succeq 0,$$

$$\begin{bmatrix}Z^{\prime}&\Sigma-\Sigma_{data}\\ \Sigma-\Sigma_{data}&I\end{bmatrix}\succeq 0,$$

$$\widehat{A},\widehat{B},\widehat{W}\text{ as in \eqref{eq:linear state-space ID}},$$

where $\gamma^{\prime}>0$ and the optimal control is recovered from the optimal values $\Sigma_{\star}$ and $L_{\star}$ through $K_{\star}=L_{\star}\Sigma_{\star}^{-1}$.

The user-defined weight $\gamma^{\prime}$ plays an analogous role to $\epsilon$ in the discussion succeeding Problem 4, in determining the exploration vs exploitation balance of the new control design.

Only the state data distribution is used in Problems 4 and 5. This puts no emphasis on the input data distribution and carries the implicit assumption that the underlying true system (1) is linear with constant coefficients in $u_{k}$. That is, the seen inputs have the same effect regardless of the system state at their occurrence. In the following subsection, however, we enforce data conformation in the state-input joint distribution instead, which is more suited to general nonlinear systems.

4.2 Conforming to the joint state-input data

We denote the steady-state state-input joint density of the new design and the one estimated from the data by $\mathcal{N}_{des}=\mathcal{N}(\mu_{des},\Gamma_{des})$ and $\mathcal{N}_{data}=\mathcal{N}(\mu_{data},\Gamma_{data})$, respectively. The Gaussian property of these densities and $\mu_{des}=\mu_{data}=0$ follow from Assumption 2. Therefore, these densities are fully characterized by their covariance matrices.

Since $u_{k}=Kx_{k}+v_{k}$, the design covariance matrix satisfies

$$\Gamma_{des} =\begin{bmatrix}\lim_{k\to\infty}\mathbb{E\,}x_{k}x_{k}^{\top}&\lim_{k\to\infty}\mathbb{E\,}x_{k}u_{k}^{\top}\\ \lim_{k\to\infty}\mathbb{E\,}u_{k}x_{k}^{\top}&\lim_{k\to\infty}\mathbb{E\,}u_{k}u_{k}^{\top}\end{bmatrix}$$

$$=\begin{bmatrix}\Sigma&\Sigma K^{\top}\\ K\Sigma&K\Sigma K^{\top}+V\end{bmatrix},$$ (12)

where $\Sigma$ is as defined in (9). The empirical covariance matrix satisfies

$$\begin{aligned} \Gamma_{data}&\approx\frac{1}{N+1}\left[D-\mu_{data}\right]\left[D-\mu_{data}\right]^{\top}\\ &=\begin{bmatrix}\Sigma_{data}&H_{data}\\ H_{data}^{\top}&M_{data}\end{bmatrix},\end{aligned}$$ (13)

where $H_{data}=(N+1)^{-1}XU^{\top}$ and $M_{data}=(N+1)^{-1}UU^{\top}$.

It is not clear whether one can enforce the closeness of $\Gamma_{data}$ to $\Gamma_{des}$, as was done in Problems 4 and 5, via LMIs in $\Sigma$ and $K$ (or $L$). Instead, we use the Kullback–Leibler divergence, notated $d_{KL}(\cdot\mid\mid\cdot)$, between the densities $\mathcal{N}_{des}$ and $\mathcal{N}_{data}$, which reduces to an expression in terms of their covariances. This expression, under some relaxation, can be posed as an affine regularization term and LMI constraints in $\Sigma$ and $L$.

Although the KL divergence is not a full-fledged metric, it is non-negative, and zero if and only if the two distributions are equal (almost everywhere) [37, Sec. 8.6]. It is a measure of the inefficiency of assuming one distribution when the true distribution is another, and it has been used in the contexts of exploration vs. exploitation in artificial intelligence and RL [35] and in constructing ambiguity sets in distributionally robust optimization techniques [29].

Lemma 1.

[30] The Kullback–Leibler divergence metric between $\mathcal{N}_{des}$ and $\mathcal{N}_{data}$ is

$$d_{KL}(\mathcal{N}_{des}\mid\mid\mathcal{N}_{data})=\int_{\mathbb{R}^{r_{x}+r_{u}}}\mathcal{N}_{des}(x)\log\frac{\mathcal{N}_{des}(x)}{\mathcal{N}_{data}(x)}dx$$

$$=\frac{1}{2}\Bigg{[}\textit{tr}(\Gamma_{data}^{-1}\Gamma_{des})+(\mu_{des}-\mu_{data})^{\top}\Gamma_{data}^{-1}(\mu_{des}-\mu_{data})$$

$$\hskip 42.67912pt-(r_{x}+r_{u})+\log\frac{\text{det }\Gamma_{data}}{\text{det }\Gamma_{des}}\Bigg{]}.$$ (14)

∎

We drop the second term in (14), since $\mu_{des}=\mu_{data}=0$, and the additive constant $-(r_{x}+r_{u})$, since it does not alter the optimal control. The modified KL-based regularization term is then

$$\bar{d}_{KL}(\mathcal{N}_{des}\mid\mid\mathcal{N}_{data})=\frac{1}{2}\Bigg{[}\textit{tr}(\Gamma_{data}^{-1}\Gamma_{des})+\log\frac{\text{det }\Gamma_{data}}{\text{det }\Gamma_{des}}\Bigg{]}.$$

We now employ an approximation to reshape $\bar{d}_{KL}(\mathcal{N}_{des}\mid\mid\mathcal{N}_{data})$ into a form amenable to the linear control design machinery, using the following results.

Lemma 2.

The algebraic equality

$$\log\frac{\text{det }\Gamma_{data}}{\text{det }\Gamma_{des}}=\textit{tr}\log\left(\Gamma_{data}\Gamma_{des}^{-1}\right)$$ (15)

holds.

Proof.

This is immediate from a result in linear algebra for positive definite matrices. We can write

$$\log\frac{\text{det }\Gamma_{data}}{\text{det }\Gamma_{des}} =\log\text{det }\Gamma_{data}-\log\text{det }\Gamma_{des},$$

$$=^{(1)}\log\text{det }\Gamma_{data}+\log\text{det }\Gamma_{des}^{-1},$$

$$=^{(2)}\textit{tr}\log\Gamma_{data}+\textit{tr}\log\Gamma_{des}^{-1},$$

$$=^{(3)}\textit{tr}\log\left(\Gamma_{data}\Gamma_{des}^{-1}\right),$$

where $(1),(2),(3)$ are implied by the positive definiteness of $\Gamma_{des}$ and $\Gamma_{data}$, and the identity $\log\text{det}\,\cdot=\textit{tr}\log\cdot$ for nonsingular matrices. ∎

Toward an approximation of (15) that is convex in $\Gamma_{des}$, suppose $\Gamma_{data}\Gamma_{des}^{-1}\approx I$. That is, suppose $\Gamma_{data}\Gamma_{des}^{-1}=I+\Theta$ for a small perturbation matrix $\Theta$. The matrix logarithm can be written in terms of its Taylor expansion as

$$\log\left(I+\Theta\right)=\Theta-\frac{1}{2}\Theta^{2}+\frac{1}{3}\Theta^{3}+\ldots,$$

which, for a first-order truncation, satisfies

$$\log\left(\Gamma_{data}\Gamma_{des}^{-1}\right)=\log\left(I+\Theta\right) \approx\Theta=\Gamma_{data}\Gamma_{des}^{-1}-I.$$ (16)

Substituting the approximation (16), after dropping the identity (additive constant), in (15), we can use a surrogate version of $\bar{d}_{KL}$, call it $F$, as a regularization term. For a fixed $\Gamma_{data}$,

$$F(\Gamma_{des})=\textit{tr}\left(\Gamma_{data}^{-1}\Gamma_{des}+\Gamma_{data}\Gamma_{des}^{-1}\right),$$ (17)

where the $1/2$ factor in $\bar{d}_{KL}$ is dropped and will be replaced by a hyperparameter $\gamma>0$. The new version of the cost (8), with $F$ as a regularization term and $\gamma>0$ a user-defined hyperparameter, is now designated by

$$J_{F} =\textit{tr}\left(Q\Sigma\right)+\textit{tr}\left(RK\Sigma K^{\top}\right)+\gamma F.$$ (18)

Next we show that the regularization term $F$, although an approximation of the KL divergence, has favorable properties for our purposes.

4.3 Properties of $F$

We show that $\Gamma_{des}=\Gamma_{data}$ is the unique minimizer of the regularization term $F$, or equivalently of $J_{F}$ as $\gamma\to\infty$. That is, the regularization term enforces conforming to the data distribution, even though $F$ is not the KL divergence, after adopting the approximation (16). We take as given the arithmetic and geometric means inequality, which we state as the following lemma.

Lemma 3.

Let $\lambda>0$. Then $\lambda+1/\lambda\geq 2$, with equality if and only if $\lambda=1$. ∎

From this, we conclude the subsequent lemma about positive definite matrices.

Lemma 4.

Let $E$ be an $n\times n$ symmetric positive definite matrix. Then $E=\mathbb{I}$ is the unique minimizer of $\bar{F}(E)=\textit{tr}\left(E+E^{-1}\right)$.

Proof.

Since $E$ is symmetric positive definite, we may write the unitary diagonalization $E=\Psi\Lambda\Psi^{T}$, where $\Lambda\succ 0$ is diagonal and $\Psi$ is unitary. Therefore, by algebraic identities $\bar{F}(E)=\textit{tr}\left(\Psi\left(\Lambda+\Lambda^{-1}\right)\Psi^{T}\right)=\textit{tr}\left(\left(\Lambda+\Lambda^{-1}\right)\Psi^{T}\Psi\right)=\textit{tr}\left(\Lambda+\Lambda^{-1}\right)$. That is,

$$\bar{F}(E)=\textit{tr}\left(\Lambda+\Lambda^{-1}\right)=\sum_{i}^{n}\left(\lambda_{i}+\frac{1}{\lambda_{i}}\right),$$

where $\lambda_{i}>0,\,i=1,\ldots,n,$ are the eigenvalues of $E$. The unique minimum happens when each $\lambda_{i}=1$, per Lemma 3. That is, when $E=\Lambda=\mathbb{I}$ ($\mathbb{I}$ is the only positive definite matrix of eigenvalues of $1$), it is the unique minimizer of $\bar{F}$. ∎

Theorem 1.

$\Gamma_{des}=\Gamma_{data}$ is the global minimizer of

$$F(\Gamma_{des})=\textit{tr}\left(\Gamma_{data}^{-1}\Gamma_{des}+\Gamma_{data}\Gamma_{des}^{-1}\right).$$

Proof.

By the linearity and the cyclic property of the trace, we can rewrite $F$ as follows:

$$F(\Gamma_{des})=\textit{tr}\left(\Gamma_{data}^{-\frac{1}{2}}\Gamma_{des}\Gamma_{data}^{-\frac{1}{2}}+\Gamma_{data}^{\frac{1}{2}}\Gamma_{des}^{-1}\Gamma_{data}^{\frac{1}{2}}\right).$$

If we let $E=\Gamma_{data}^{-\frac{1}{2}}\Sigma_{des}\Sigma_{data}^{-\frac{1}{2}}$, then we have

$$\bar{F}(E)=\textit{tr}\left(E+E^{-1}\right).$$

Showing that $\Gamma_{des}=\Gamma_{data}$ is the unique global minimizer of $F$ equates to showing that $E=\mathbb{I}$ is the unique global minimizer of $\bar{F}$, which is true by Lemma 4. ∎

Next we show that the regularization term $F$ not only has a unique global minimizer that we seek to enforce, that is, $\Gamma_{des}=\Gamma_{data}$, but is also computationally appealing, convex in particular.

Theorem 2.

$F$ is convex in $\Gamma_{des}$.

Proof.

The domain—the set of all positive definite matrices—is convex. The first term of $F(\Gamma_{des})$ is linear in $\Gamma_{des}$, while the second term

$$\textit{tr}\left(\Gamma_{data}\Gamma_{des}^{-1}\right) =\textit{tr}\left(\Gamma_{data}^{\frac{1}{2}}\Gamma_{des}^{-1}\Gamma_{data}^{\frac{1}{2}}\right)=\sum_{i}\xi_{i}^{T}\Gamma_{des}^{-1}\xi_{i},$$

where $\xi_{i}$ is the $i$th column of $\Gamma_{data}^{\frac{1}{2}}$. Each entry in the sum is convex in $\Gamma_{des}$, as shown in [11, p. 76]. ∎

Notice that $\frac{1}{2}F$ is the Jeffreys divergence, a discrepancy metric between distributions [30] (Gaussians in this case). The above derivations highlight its relation to the KL divergence, which is more familiar in the exploration vs exploitation context.

The above results are concerned with the regularization term as a function of $\Gamma_{des}$. However, the original LQR problem, whether Problem 2 or 3, has $\Sigma$ and $L$ ($L=K\Sigma$) as decision variables. We show next that, under some relaxation, $F$ can be related back to the original decision variables affinely in the cost together with LMI constraints.

4.4 Representing $F$ in terms of $\Sigma$ and $L$

The first term of $F$, $\textit{tr}\left(\Gamma_{data}^{-1}\Gamma_{des}\right)$, can be upper-bounded by $\textit{tr}\left(\Gamma_{data}^{-1}Z_{1}\right)$, where $Z_{1}\succeq\Gamma_{des}$, or equivalently, committing to the change of variables $L=K\Sigma$,

$$Z_{1}\succeq\Gamma_{des} =\begin{bmatrix}\Sigma&\Sigma K^{\top}\\ K\Sigma&K\Sigma K^{\top}+V\end{bmatrix},$$

$$=\begin{bmatrix}\Sigma\\ L\end{bmatrix}\Sigma^{-1}\begin{bmatrix}\Sigma&L^{\top}\end{bmatrix}+\mathcal{V},$$

which can be equivalently described by the LMI

$$\begin{bmatrix}Z_{1}-\mathcal{V}&\begin{bmatrix}\Sigma\\ L\end{bmatrix}\\ \begin{bmatrix}\Sigma&L^{\top}\end{bmatrix}&\Sigma\end{bmatrix}\succeq 0,$$ (19)

where $\mathcal{V}=\text{block-diag}(0_{r_{x}\times r_{x}},V)$.

The second term is more involved since it contains the inverse of $\Gamma_{des}$. Applying the inverse of a partitioned matrix [24], we have

$$\Gamma_{des}^{-1}=\begin{bmatrix}\Sigma^{-1}+K^{\top}V^{-1}K&-K^{\top}V^{-1}\\ -V^{-1}K&V^{-1}\end{bmatrix},$$

which exists and is positive definite, since $\Gamma_{des}\succ 0$. Hence, the second term of $F$

$$\textit{tr}\left(\Gamma_{data}\Gamma_{des}^{-1}\right) =\textit{tr}\left(\Sigma_{data}\Sigma^{-1}+\Sigma_{data}K^{\top}V^{-1}K\right)+$$

$$\hskip-28.45274pt\textit{tr}\left(M_{data}V^{-1}-2H^{\top}_{data}K^{\top}V^{-1}\right),$$

$$=\textit{tr}\left(\Sigma_{data}\Sigma^{-1}\right)+\textit{tr}\left(V^{-1}K\Sigma_{data}K^{\top}\right)+$$

$$\hskip-28.45274pt\textit{tr}\left(M_{data}V^{-1}\right)+\textit{tr}\left(-2K^{\top}V^{-1}H^{\top}_{data}\right).$$ (20)

Lemma 5.

We can rewrite the following term as

$$\textit{tr}\left(\Sigma_{data}K^{\top}V^{-1}K-2K^{\top}V^{-1}H^{\top}_{data}\right)=$$ (21)

$$\textit{tr}\Big{(}V^{-1}\left[K\Sigma-H_{data}^{\top}\Sigma_{data}^{-1}\Sigma\right]\Sigma^{-1}\Sigma_{data}\Sigma^{-1}\times$$

$$\hskip-28.45274pt\left[K\Sigma-H_{data}^{\top}\Sigma_{data}^{-1}\Sigma\right]^{\top}-V^{-1}H_{data}^{\top}\Sigma_{data}^{-1}H_{data}\Big{)}.$$

Proof.

Using the linearity and the cyclic property of the trace, we reduce the right-hand side to the left one. ∎

After dropping the additive constants in $\Sigma$ and $K$, from (4.4) and (LABEL:eq:_completing_the_squares), and toward forming LMIs, we relax the right-hand side of (LABEL:eq:_completing_the_squares) by assuming $\Sigma\approx\Sigma_{data}$; thus, $\Sigma^{-1}\Sigma_{data}\Sigma^{-1}\approx\Sigma^{-1}$. We use the extra variable $Z_{2}$ such that

$$Z_{2}\succeq\left[K\Sigma-H_{data}^{\top}\Sigma_{data}^{-1}\Sigma\right]\Sigma^{-1}\left[K\Sigma-H_{data}^{\top}\Sigma_{data}^{-1}\Sigma\right]^{\top},$$

which, equivalently, as an LMI, is

$$\begin{bmatrix}Z_{2}&L-H_{data}^{\top}\Sigma_{data}^{-1}\Sigma\\ \left[L-H_{data}^{\top}\Sigma_{data}^{-1}\Sigma\right]^{\top}&\Sigma\end{bmatrix}\succeq 0.$$ (22)

We show empirically in Section 6 that, even after the above relaxation, the consistency with state-input joint distribution is still effectively enforced.

The term $\textit{tr}\left(\Sigma_{data}\Sigma^{-1}\right)$ can also be described by an LMI by using an extra variable $Z_{3}\succeq\Sigma^{-1}$. Hence,

$$\begin{bmatrix}Z_{3}&I\\ I&\Sigma\end{bmatrix}\succeq 0.$$ (23)

4.5 The joint state-input data-conforming data-driven LQR

Toward a state-input data-conforming LQR control, we adjust Problem 3 by including the regularization term $F$ and the accompanying extra variables and LMIs defined in (19), (22), and (23).

Problem 6.

Data-conforming (jointly in the state-input) data-driven LQR:

$$\min_{\Sigma,\,L,\,Z_{0},\,Z_{1},\,Z_{2},\,Z_{3}}\textit{tr}\left(Q\Sigma\right)+\textit{tr}\left(RZ_{0}\right)+$$

$$\hskip 14.22636pt\gamma\Big{\{}\textit{tr}\left(\Gamma_{data}^{-1}Z_{1}\right)+\textit{tr}\left(V^{-1}Z_{2}\right)+\textit{tr}\left(\Sigma_{data}Z_{3}\right)\Big{\}}$$

$$\text{s.t. }\Sigma\succ 0,\quad\begin{bmatrix}Z_{0}&L\\ L^{\top}&\Sigma\end{bmatrix}\succeq 0,$$

$$\begin{bmatrix}\Sigma-\widehat{W}-\widehat{B}V\widehat{B}^{\top}&\widehat{A}\Sigma+\widehat{B}L\\ \Sigma\widehat{A}^{\top}+L^{\top}\widehat{B}^{\top}&\Sigma\end{bmatrix}\succeq 0,$$

$$\begin{bmatrix}Z_{1}-\mathcal{V}&\begin{bmatrix}\Sigma\\ L\end{bmatrix}\\ \begin{bmatrix}\Sigma&L^{\top}\end{bmatrix}&\Sigma\end{bmatrix}\succeq 0,\quad\begin{bmatrix}Z_{3}&I\\ I&\Sigma\end{bmatrix}\succeq 0,$$

$$\begin{bmatrix}Z_{2}&L-H_{data}^{\top}\Sigma_{data}^{-1}\Sigma\\ \left[L-H_{data}^{\top}\Sigma_{data}^{-1}\Sigma\right]^{\top}&\Sigma\end{bmatrix}\succeq 0,$$

$$\widehat{A},\widehat{B},\widehat{W}\text{ as in \eqref{eq:linear state-space ID}}.$$

The optimal control is recovered from the optimal values $\Sigma_{\star}$ and $L_{\star}$ through $K_{\star}=L_{\star}\Sigma_{\star}^{-1}$.

Remark 1.

Notice that if the zero mean condition is to be relaxed in Assumption 2, the corresponding term with $\mu_{des},\mu_{data}\neq 0$ in (14) and the cost (2) are quadratic in the variables $\bar{x}$ and $\bar{u}$, where $\mu_{des}=\left[\bar{x}^{\top},\,\bar{u}^{\top}\right]^{\top}$, and subject to the linear constraint given by steady-state mean dynamics

$$\bar{x}=(A+BK)\bar{x}+B\bar{u}.$$

One therefore can solve Problem 6 to find the optimal gain $K$ and covariance $\Sigma$ and then, as a second stage, form a quadratic program in $\mu_{des}$.

This is common in scenarios where the current operating point during the data collection $\mu_{data}=\left[\bar{x}^{\top}_{data},\,\bar{u}^{\top}_{data}\right]^{\top}$ is different from the intended one. In such a case, we center the data and the system (1) around the latter and incorporate the quadratic program discussed above.

Remark 2.

As mentioned previously, the identification step, inside Problems 3, 4, 5, and 6, resembles a naive certainty equivalence approach [16]. One can enhance the robustness of these problems by incorporating a convex hull or a ball, according to some systems’ norm, centered around $(\widehat{A},\widehat{B})$, then search for a controller that stabilizes them all, while minimizing the cost and satisfying other constraints. This adjustment, say, to Problem 6, can be done through more LMI constraints [10] in the same decision variables. One also can, again, through LMI constraints, enforce robustness under state measurement noise (state is measured through $\zeta_{k}$, where $\zeta_{k}=x_{k}+\eta_{k}$, for some white noise $\eta_{k}$) as in [14] for when the covariance of the state evolution disturbance $W\neq 0$, as in [38].

Remark 3.

The robustness results referred to in Remark 2 hold when the underlying true system (1) is linear. This motivates a new interpretation of our data-conforming framework.

Suppose data has been collected, a system has been identified, and a convex hull or a ball $\mathbb{B}$ has been constructed around this system’s estimate. When the new robust (e.g., according to [10]) controller is applied, new regions of the joint state-input space might be visited, due to the distributional shift, and new modes of nonlinearity activated. If the new system is to be identified and a new convex hull or ball to be constructed, they might be different from $\mathbb{B}$ used in the design process.

In other words, because of the nonlinearity of the underlying model in our case, the application of the new control law equates to the transformation of $\mathbb{B}$. Hence, this invalidates the premise on which these robust control design approaches were built.

Our data-conforming framework helps in “damping” this transformation in $\mathbb{B}$; and therefore, when incorporated to the data-driven robust design approaches [10, 38], it increases their validity via solidifying their premises.

5 On the optimality conditions of the proposed problems

In the proposed SDPs, we relaxed the Lyapunov equality (10) into an inequality such that we can represent it as an LMI in the decision variables. This relaxation is of no effect to the standard LQR SDP, as we will show shortly, but may have some consequences to the proposed data-conforming approaches.

Notice that Problems 3, 5 and 6, before the change of variables and before applying the equivalent LMI formulations, can be put in the following form (we use $A,B,W$ instead of $\widehat{A},\widehat{B},\widehat{W}$ for ease of notation)

$$\min_{\Sigma,\,K}\textit{tr}\left(Q\Sigma\right)+\textit{tr}\left(RK\Sigma K^{\top}\right)+\Xi\left(\Sigma,K\right)$$

s.t.

$$\Sigma\succeq\left[A+BK\right]\Sigma\left[A+BK\right]^{\top}+W+BVB^{\top},$$ (24)

where $\Xi\left(\Sigma,K\right)=0$ in Problem 3,

$$\Xi\left(\Sigma,K\right)=\gamma^{\prime}\textit{tr}\left(\left[\Sigma-\Sigma_{data}\right]\left[\Sigma-\Sigma_{data}\right]^{\top}\right)$$

in Problem 5, and

$$\Xi\left(\Sigma,K\right)=\gamma\Big{(}\textit{tr}\left(\Sigma_{data}\Sigma^{-1}\right)+\textit{tr}\left(V^{-1}K\Sigma_{data}K^{\top}\right)$$

$$\hskip 28.45274pt+\textit{tr}\left(-2K^{\top}V^{-1}H^{\top}_{data}\right)$$

$$+\textit{tr}\left(\mathcal{A}\Sigma+\mathcal{B}K\Sigma\right)+\textit{tr}\left(\mathcal{B}^{\top}\Sigma K^{\top}+\mathcal{C}K\Sigma K^{\top}\right)\Big{)}$$

(up to an additive constant in $K,\Sigma$) in Problem 6, where the matrices in calligraphic font denote the block matrices of the appropriate size of $\Gamma_{data}^{-1}$, that is,

$$\Gamma_{data}^{-1}=\begin{bmatrix}\mathcal{A}&\mathcal{B}\\ \mathcal{B}^{\top}&\mathcal{C}\end{bmatrix}.$$

In the above problems, the Langrangian is given by

$$\mathcal{L}(\Sigma,K,\Upsilon)=\textit{tr}\left(Q\Sigma\right)+\textit{tr}\left(RK\Sigma K^{\top}\right)+\Xi\left(\Sigma,K\right)+$$

$$\textit{tr}\Bigg{(}\Upsilon\Big{[}\left[A+BK\right]\Sigma\left[A+BK\right]^{\top}-\Sigma+W+BVB^{\top}\Big{]}\Bigg{)},$$

where $\Upsilon$ is the dual $r_{x}\times r_{x}$ matrix variable corresponding to the Lyapunov inequality constraint. From the dual feasibility and the zero gradient Karush-Kuhn-Tucker (KKT) [11] optimality conditions (derivatives of traces of matrices [31]), the optimal solution $(K_{\star},\Sigma_{\star},\Upsilon_{\star})$ has to satisfy [22]:

$$\Upsilon_{\star}\succeq 0,$$

$$\frac{\partial}{\partial\Sigma}:\,Q+K_{\star}^{\top}RK_{\star}+\left[A+BK_{\star}\right]^{\top}\Upsilon_{\star}\left[A+BK_{\star}\right]-\Upsilon_{\star}$$

$$\hskip 56.9055pt+\frac{\partial}{\partial\Sigma}\Xi\left(\Sigma_{\star},K_{\star}\right)=0.$$

Notice that in the standard LQR formulation, Problem 3, when $\Xi\left(\Sigma,K\right)=0$, $\Upsilon_{\star}$ is nothing but the observability-type Gramian $P$ in (7), that is, $\Upsilon_{\star}=P\succ 0$. The positive definiteness of $\Upsilon_{\star}$ then implies that in the complementary slackness condition,

$$\Upsilon_{\star}\Big{[}\left[A+BK_{\star}\right]\Sigma_{\star}\left[A+BK_{\star}\right]^{\top}-\Sigma_{\star}+W+BVB^{\top}\Big{]}=0,$$

the constraint has to be active to satisfy the above condition, that is, $(\Sigma_{\star},K_{\star})$ satisfy equation (10),

$$\left[A+BK_{\star}\right]\Sigma_{\star}\left[A+BK_{\star}\right]^{\top}-\Sigma_{\star}+W+BVB^{\top}=0,$$

or in other words, $\Sigma_{\star}$ is exactly the steady-state covariance of the system state we expect to see in future data when applying the law $u_{k}=K_{\star}x_{k}+v_{k}$.

This, however, is not generally true for Problems 5 and 6, where the Lyapunov inequality may be inactive in certain situations. The zero gradient condition for these two problems, respectively,

$$Q+K_{\star}^{\top}RK_{\star}+\left[A+BK_{\star}\right]^{\top}\Upsilon_{\star}\left[A+BK_{\star}\right]-\Upsilon_{\star}$$

$$\hskip 56.9055pt+2\gamma^{\prime}\left[\Sigma_{\star}-\Sigma_{data}\right]=0,$$

and

$$Q+\gamma\mathcal{A}+K_{\star}^{\top}\left[R+\gamma\mathcal{C}\right]K_{\star}+\gamma K_{\star}^{\top}\mathcal{B}^{\top}+\gamma\mathcal{B}K_{\star}$$

$$\hskip 0.0pt-\gamma\Sigma_{\star}^{-1}\Sigma_{data}\Sigma_{\star}^{-1}+\left[A+BK_{\star}\right]^{\top}\Upsilon_{\star}\left[A+BK_{\star}\right]-\Upsilon_{\star}=0.$$

Using the complementary slackness condition again, to guarantee the satisfaction of the equality (10) in each problem, the terms

$$Q+K_{\star}^{\top}RK_{\star}+2\gamma^{\prime}\left[\Sigma_{\star}-\Sigma_{data}\right],$$ (25)

and

$$Q+\gamma\mathcal{A}+K_{\star}^{\top}\left[R+\gamma\mathcal{C}\right]K_{\star}+\gamma K_{\star}^{\top}\mathcal{B}^{\top}+\gamma\mathcal{B}K_{\star}$$ (26)

$$-\gamma\Sigma_{\star}^{-1}\Sigma_{data}\Sigma_{\star}^{-1},$$

have to be positive definite, which might not be the case if $\Sigma_{data}$ is significantly larger than $\Sigma_{\star}$. However, we show next that even if the terms (25) and (26) are not positive definite, it is still beneficial to have the inequality (24) inactive.

By the contrapositive argument, if the inequality constraint (24) is inactive, the terms (25) and (26) are either indefinite, or negative semi-definite. This can be implied when some eigenvalues of $\Sigma_{data}$ are significantly larger than some of $\Sigma_{\star}$’s. In such case, two important observations can be made: (i) when $\Sigma_{data}$ has larger eigenvalues than $\Sigma_{\star}$, the data show that aggressive exploration has been done along the direction of the corresponding eigenvectors, while less aggressive (conservative) control design is achieved across the same directions, (ii) suppose $\Sigma_{actual}$ is the solution of the Lyapunov equation (10) when $K=K_{\star}$, then $\Sigma_{actual}\preceq\Sigma_{\star}$, that is, the actual achieved design is more conservative in its exploration tendencies than the intended design. The inequality $\Sigma_{actual}\preceq\Sigma_{\star}$ holds because when the Lyapunov inequality (24) is inactive, it is implied that the optimal solution $\Sigma_{\star}$ to the SDP is not the fixed point solution of the Lyapunov equation (10), hence, $\Sigma_{\star}$ is not the steady-state covariance, but rather an upper bound. That is, $\Sigma_{actual}$ is the smallest covariance that satisfies the Lyapunov inequality (24). This can also be seen when writing the optimality conditions of the SDP with the cost $\textit{tr}\left(\Sigma\right)$ and subject to the Lyapunov inequality (24).

The above observations show that the case of inactive constraints corresponds to a more conservative achieved control. We show empirically in the next section that the Lyapunov inequality constraint (24) is mostly active in our simulations, and that our approaches still yield favorable results for nonlinear systems, over the certainty equivalence LQR case, even when the Lyapunov constraint is inactive.

6 Numerical simulations

The scalability of our proposed data-conforming formulations is self-evident; as convex SDPs with affine costs and LMI constraints, they can be scaled to handle problems with high state-input dimensions. Instead, we choose to work with a simple, yet telling problem, crafted to further clarify and illustrate our data-conforming paradigm.

Suppose the data-generating system (1) is of the form

$$\begin{pmatrix}x_{1,k+1}\\ x_{2,k+1}\end{pmatrix}=x_{k+1}=\begin{bmatrix}.98&.1\\ 0&.95\end{bmatrix}x_{k}+\begin{bmatrix}0\\ .1\end{bmatrix}u_{k}+w_{k},$$ (27)

with the noise $w_{k}$ having a covariance $W=\text{diag}(0.2,\,0.1)$. Toward solving Problems 3 and 5, an experimental simulation was conducted using $N=2000$ collected state/input samples. The control $u_{k}=K_{0}x_{k}+v_{k}$, where $K_{0}=\left[-0.2,\,-9.0\right]$ is stabilizing, and $v_{k}$ is a white noise of zero mean and covariance $V=0.5$ (meeting the PE condition, Assumption 1).

Using these simulation samples to form the data matrices (4), we then solve Problems 3 and 5 and show their results in Figure 1. Each control gain resulting from these problems is used in a control law of the form $u_{k}=Kx_{k}+v_{k}$ and is run in a new simulation of the same number of time steps ($2000$). Figure 1 shows the resulting design distributions compared with the learning data. Notice that when $\gamma^{\prime}$ starts having bigger and bigger values, the state design distribution converges to the state data distribution because data conformity is more heavily weighted in the control design.

Notice also in Figure 1 (a), that the (empirical) support of the state distribution of the new designed closed-loop system goes beyond that of the data. Suppose, for instance, that these not well-explored regions contain issues such as nonlinearities, discontinuities, or some safety or design condition violations. These issues, if not known and modeled by the operator in the control design, are also not discovered by the initial experiment. Hence, a non-data-conforming control, as the one in Figure 1 (a), can lead to activating these issues due to its inherent generalization beyond data. On the other hand, in (b), (c) and (d), the new controller conforms to the learning data to different levels dictated by the values of $\gamma^{\prime}$ in each.

Figure 1: The x-axis and y-axis of each figure correspond to the state space, $x_{1,k}$ and $x_{2,k}$, respectively. The orange squares are the state data resulting from the initial simulation experiment, using $K_{0}$. The blue x-crosses are the state data resulting from solving Problem 3 (or, equivalently, Problem 5 with $\gamma^{\prime}=0$) in (a), and Problem 5 with $\gamma^{\prime}=5$ in (b), $\gamma^{\prime}=20$ in (c), and $\gamma^{\prime}=100$ in (d).

To better illustrate the effect of unknown/unmodeled nonlinearities, we adjust the dynamics (27) such that it is now

$$x_{k+1}=\begin{bmatrix}.98&.1\\ 0&.95\end{bmatrix}x_{k}+\begin{pmatrix}\theta x_{2,k}^{2}\\ 0\end{pmatrix}+\begin{bmatrix}0\\ .1\end{bmatrix}u_{k}+w_{k},$$ (28)

where $\theta=1/9$. Notice now that the nonlinearity $\theta x_{2,k}^{2}$ is relatively larger in the not well-explored regions. Naively applying Problem 3 can lead to unexpected behavior and possibly to instability.

Again, we simulate the new nonlinear system with an equivalent control $u_{k}=K_{0}x_{k}+v_{k}$ as in the previous case and for $N=2000$ time steps. Then, at the last time step, we evaluate the new control laws, using the recorded past $2000$ data points, and apply each control law immediately in feedback for a new $2000$ time steps. The results of each simulation are illustrated in Figure 2. The instability resulting from relying on Problem 3 (or, equivalently, Problem 5 with $\gamma^{\prime}=0$) is shown in Figure 2: (1). In contrast, Figure 2: (2), (3), and (4) correspond to Problem 5 with $\gamma^{\prime}=5,\,20,\,100$, respectively. The latter three cases enforce conforming to data and hence avoid unexplored regions, in particular those where the nonlinearity is dominant and dangerous.

Figure 2: Values of $x_{1,k}$ for the $N=2000$ samples of the initial experiment then under the control law resulting from solving (1) Problem 3 (or, equivalently, Problem 5 with $\gamma^{\prime}=0$), and Problem 5 with (2) $\gamma^{\prime}=5$, (3) $\gamma^{\prime}=20$, and (4) $\gamma^{\prime}=100$.

The steps of the above procedure, from (i) an initial experiment of $N=2000$ bounded data, to (ii) the control design using the above four problems (Problems 3, and Problem 5 with $\gamma^{\prime}=5,\,20,\,100$) and then (iii) apply each control in feedback for $2000$ time steps, are repeated ((i)–(iii)) for $1000$ repetitions. Then, we calculate the percentage of the stable simulations resulting from each procedure. The stability of each simulation is decided by the boundedness within some high threshold ($50$) in each state coordinate and pointwise in the $2000$ time steps of step (iii). The Lyapunov inequality constraint (24) was active in almost all of these simulations. The percentages of the stable simulations are shown in Table 1, further highlighting the importance of data-conforming control design.

Table 1: Percentages of stable simulations (out of 1,000 simulations)
Problem 3	Problem 5	Problem 5	Problem 5
(Prob. 5, $\gamma^{\prime}=0$)	($\gamma^{\prime}=5$)	($\gamma^{\prime}=20$)	($\gamma^{\prime}=100$)
21.7%	99.2%	99.9%	100.0%

The examples in (27) and (28) share, in their structure, that the input enters the system linearly and with constant coefficients. Therefore, the effect of the input is the same on the state evolution, regardless of the contemporaneous state of the system. This is not true in general, and the effect of the input may be coupled with the current state—for instance, in bilinear systems. The following is a modification of (28) with a coupled state-input term:

$$x_{k+1} =\begin{bmatrix}.98&.1\\ 0&.95\end{bmatrix}x_{k}+\begin{pmatrix}\theta x_{2,k}^{2}\\ 0\end{pmatrix}+$$ (29)

$$\hskip 42.67912pt\begin{bmatrix}0\\ 0.1+\theta\tanh{x_{1,k}}\end{bmatrix}u_{k}+w_{k}.$$

Similar to the previous two examples, we run an experimental simulation of this system for $N=2000$ time steps and collect the input and state data as in (4). At the final time step $k=N=2000$, we use the collected data in solving Problems 3, 5, and 6, then use the resulting control laws for another $2000$ samples. The results of these simulations are shown in Figure 3, showing the importance of introducing Problem 6 to account for the coupled state-input effect.

Similar to the construction of Table 1, we run the above procedure for $1000$ repetitions of bounded experiments and calculate the percentages of stable simulations resulting from each of the four control design procedures (Problem 3, Problem 5 (with $\gamma^{\prime}=20,\,100$) and Problem 6 (with $\gamma=10$)). The results are recorded in Table 2, which shows the importance of using the joint state-input data distribution when the system’s nonlinearities contain a coupled state-input effect, even though the Lyapunov inequality constraint (24) was mostly inactive.

Table 1: Percentages of stable simulations (out of 1000 simulations)
Problem 3	Problem 5	Problem 5	Problem 6
(Prob. 5, $\gamma^{\prime}=0$)	($\gamma^{\prime}=20$)	($\gamma^{\prime}=100$)	($\gamma=10$)
0.07%	22.9%	22.5%	97.9%

We note that in all of these simulations, none of the tested problems returned any feasibility issues. Moreover, the computation time of each problem (about six milliseconds for Problem 6 on an Apple M1 Max MacBook Pro with 64 GB or RAMs) is comparable to solving a standard LQR problem.

These examples show how the concept of data-conforming control can enhance the safety of the data-driven approaches. In practice, a control engineer can start with a high value of $\gamma$ (or $\gamma^{\prime}$) then gradually reduce it, if needed, to allow some exploration beyond data. One of the important recommendations of the adaptive control literature [5, 4] is to model as much as possible, and leave any model reduction or abstraction for later. Our data-conforming framework can provide the safety and the time necessary for the control engineer to observe, model/interpret, and react to any potential activation of unknown/unmodeled nonlinearities.

Figure 3: Value of $x_{1,k}$ for the $N=2000$ samples of the initial experiment then under the control law resulting from solving Problem 3, Problem 5 with $\gamma^{\prime}=10,\,100$, and Problem 6 ($\gamma=10$). The last results in the only stabilizing control since it accounts for the joint state-input data.

7 Conclusion

This paper addresses a problem inherent in many modern data-driven and adaptive control approaches, namely, the problem of the premature and possibly false generalization beyond data, together with its consequences in terms of sudden distributional shifts in the state-input space. We present methods for mitigating this problem through enforcing consistency with the learning data, and we numerically test them. Because of the formulation of our methods as solutions to SDPs of affine costs and LMI constraints, they are computationally efficient, can scale up to systems with hundreds in dimension, and can be easily integrated with modern control design approaches.

Further work is being pursued to (i) apply the data-conformation concept in the problem of robust data-driven control, in the sense described in Remarks 2 and 3; and (ii) investigate the possibility of developing data-conforming policy gradient surrogates to dampen distributional shifts resulting from standard policy gradient methods.

Government License: The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan. http://energy.gov/downloads/doe-public-access-plan.

Appendix A Appendix

A.1 The cost in the $P-$parametrization

Starting from the cost (2),

$$J=\lim_{k\to\infty}\mathbb{E\,}\left\{x_{k}^{\top}Qx_{k}+u_{k}^{\top}Ru_{k}\right\},$$

substituting $u_{k}=Kx_{k}+v_{k}$, and using the fact that $v_{k}$ is a zero mean white noise,

$$J=\lim_{k\to\infty}\mathbb{E\,}\left\{x_{k}^{\top}\left[Q+K^{\top}RK\right]x_{k}+v_{k}^{\top}Rv_{k}\right\}.$$

Using the cyclic property of the trace and the linearity of $\mathbb{E\,}$, $\mathbb{E\,}\left\{v_{k}^{\top}Rv_{k}\right\}=\textit{tr}\left(RV\right)$, which is an additive constant (for a fixed $V$) and can be omitted from the cost. Now, using the cyclic property of the trace and the linearity of $\mathbb{E\,}$ again, we have

$$J \propto\lim_{k\to\infty}\mathbb{E\,}\left\{x_{k}^{\top}\left[Q+K^{\top}RK\right]x_{k}\right\},$$

$$=\textit{tr}\left(\lim_{k\to\infty}\mathbb{E\,}\left\{\left[Q+K^{\top}RK\right]x_{k}x_{k}^{\top}\right\}\right),$$

$$=\textit{tr}\left(\left[Q+K^{\top}RK\right]\lim_{k\to\infty}\mathbb{E\,}\left\{x_{k}x_{k}^{\top}\right\}\right),$$

$$=\textit{tr}\left(\left[Q+K^{\top}RK\right]\Sigma\right).$$ (30)

The steady-state covariance of the state is designated by $\Sigma=\lim_{k\to\infty}\mathbb{E\,}\left\{x_{k}x_{k}^{\top}\right\}$, which is also known as the controllability-type Gramian. Since

$$x_{k}=\left[A+BK\right]^{k}x_{0}+\sum_{i=0}^{k-1}\left[A+BK\right]^{k-1-i}\left(w_{i}+Bv_{i}\right),$$

and using the fact that $v_{i},\,w_{i}$ are white, of zero mean and mutually independent from each other and from $x_{0}$, this covariance can be described by

$$\Sigma =\lim_{k\to\infty}\Big{\{}\left[A+BK\right]^{k}\Sigma_{x_{0}}\left[A+BK\right]^{k,\,\top}+$$

$$\hskip-14.22636pt\sum_{i=0}^{k-1}\left[A+BK\right]^{k-1-i}\left[W+BVB^{\top}\right]\left[A+BK\right]^{k-1-i,\,\top}\Big{\}},$$

or, if $A+BK$ is hurwitz, the transient term vanishes and we have

$$\Sigma=\sum_{k=0}^{\infty}\left[A+BK\right]^{k}\left[W+BVB^{\top}\right]\left[A+BK\right]^{k,\,\top}.$$

After substituting this description of $\Sigma$ in (30), and using the cyclic property of the trace operator, we have

$$J \propto\textit{tr}\Big{(}\left[Q+K^{\top}RK\right]\times$$

$$\hskip 14.22636pt\sum_{k=0}^{\infty}\left[A+BK\right]^{k}\left[W+BVB^{\top}\right]\left[A+BK\right]^{k,\,\top}\Big{)},$$

$$=\textit{tr}\Big{(}\sum_{k=0}^{\infty}\left[A+BK\right]^{k,\,\top}\left[Q+K^{\top}RK\right]\left[A+BK\right]^{k}\times$$

$$\hskip 14.22636pt\left[W+BVB^{\top}\right]\Big{)},$$

$$=\textit{tr}\left(P\left[W+BVB^{\top}\right]\right),$$

where $P$ is the observability-type Gramian, as in (6).