Regularized and Distributionally Robust Data-Enabled Predictive Control

TL;DR

This paper introduces a distributionally robust version of data-enabled predictive control (DeePC) for stochastic systems, providing theoretical guarantees and demonstrating improved robustness through numerical case studies.

Contribution

It develops a distributionally robust formulation of DeePC, linking it to regularization techniques, and offers probabilistic performance guarantees for unknown stochastic LTI systems.

Findings

Robust DeePC coincides with a regularized optimization problem.

The approach provides probabilistic guarantees for control performance.

Numerical case studies demonstrate enhanced robustness.

Abstract

In this paper, we study a data-enabled predictive control (DeePC) algorithm applied to unknown stochastic linear time-invariant systems. The algorithm uses noise-corrupted input/output data to predict future trajectories and compute optimal control policies. To robustify against uncertainties in the input/output data, the control policies are computed to minimize a worst-case expectation of a given objective function. Using techniques from distributionally robust stochastic optimization, we prove that for certain objective functions, the worst-case optimization problem coincides with a regularized version of the DeePC algorithm. These results support the previously observed advantages of the regularized algorithm and provide probabilistic guarantees for its performance. We illustrate the robustness of the regularized algorithm through a numerical case study.

Equations

$$\{\begin{array}{l}x(t+1)=Ax(t)+Bu(t)+Ev(t)\\ y(t)=Cx(t)+Du(t)+Fv(t),\end{array}$$

$$\underset{u\in \mathcal{U}}{\inf }{\mathbb{E}}_{{\mathbb{P}}_v^{T_{\text{f}}}}\left\{f(u,y)\right\},$$

$$\{\begin{array}{l}x(t+1)=Ax(t)+Bu(t)\\ y(t)=Cx(t)+Du(t).\end{array}$$

$$\underset{u\in \mathcal{U}}{\inf }f(u,y).$$

$${\mathcal{H}}_L(w):=\left(\begin{array}{cccc}w(1)& w(2)& \dots & w(T-L+1)\\ w(2)& w(3)& \dots & w(T-L+2)\\ ⋮& ⋮& \ddots & ⋮\\ w(L)& w(L+1)& \dots & w(T)\end{array}\right).$$

$$\left(\begin{array}{c}{\hat{U}}_{\mathrm{p}}\\ {\hat{U}}_{\mathrm{f}}\end{array}\right):={\mathcal{H}}_{T_{\text{ini}}+T_{\text{f}}}(\hat{u}),\left(\begin{array}{c}{\hat{Y}}_{\mathrm{p}}\\ {\hat{Y}}_{\mathrm{f}}\end{array}\right):={\mathcal{H}}_{T_{\text{ini}}+T_{\text{f}}}(\hat{y}),$$

$$\left(\begin{array}{c}{\hat{U}}_{\mathrm{p}}\\ {\hat{Y}}_{\mathrm{p}}\\ {\hat{U}}_{\mathrm{f}}\\ {\hat{Y}}_{\mathrm{f}}\end{array}\right)g=\left(\begin{array}{c}{\hat{u}}_{\text{ini}}\\ {\hat{y}}_{\text{ini}}\\ u\\ y\end{array}\right).$$

$$\underset{g}{\mathrm{minimize}}$$

$$\left(\begin{array}{c}{\hat{U}}_{\mathrm{p}}\\ {\hat{Y}}_{\mathrm{p}}\end{array}\right)g=\left(\begin{array}{c}{\hat{u}}_{\text{ini}}\\ {\hat{y}}_{\text{ini}}\end{array}\right)$$

$${\hat{U}}_{\mathrm{f}}g\in \mathcal{U}.$$

$$\underset{g}{\mathrm{minimize}}$$

$${\hat{U}}_{\mathrm{p}}g={\hat{u}}_{\text{ini}}$$

$${\hat{U}}_{\mathrm{f}}g\in \mathcal{U},$$

$$\underset{g}{\mathrm{minimize}}$$

$${\hat{U}}_{\mathrm{p}}g={\hat{u}}_{\text{ini}}$$

$${\hat{U}}_{\mathrm{f}}g\in \mathcal{U},$$

$$\underset{g\in G}{\mathrm{minimize}}f({\hat{U}}_{\mathrm{f}}g,{\hat{Y}}_{\mathrm{f}}g)+{\lambda }_{\text{ini}}\Vert {\hat{Y}}_{\mathrm{p}}g-{\hat{y}}_{\text{ini}}{\Vert }_1.$$

$${\xi }_j:={\left(\begin{array}{cc}{Y}_{\mathrm{p}}& y_{\text{ini}}\\ Y_{\mathrm{f}}& 0\end{array}\right)}_{j,\cdot },$$

$$\underset{h\in H}{\mathrm{minimize}}$$

$$+{\lambda }_{\text{ini}}\Vert ({\hat{\xi }}_{1}h,\dots ,{\hat{\xi }}_{p{T}_{\text{ini}}}h){\Vert }_1\},$$

$$\underset{h\in H}{\mathrm{minimize}}$$

$${\mathbb{E}}_{\mathbb{P}}\{c(\xi ,h^{\star })\},$$

$$\underset{h\in H}{\inf }\underset{Q\in \hat{\mathcal{P}}}{\sup }{\mathbb{E}}_Q\{c(\xi ,h)\},$$

$${\text{d}}_{\text{W}}(Q_1,Q_2):=\underset{\Pi }{\inf }\left\{{\int }_{{\Xi }^2}\Vert {\xi }_1-{\xi }_2{\Vert }_{\text{W}}\Pi (\text{d}{\xi }_1,\text{d}{\xi }_2)\right\},$$

$$\underset{h\in H}{\inf }\underset{Q\in B_{ϵ}(\hat{\mathbb{P}})}{\sup }{\mathbb{E}}_Q\{c(\xi ,h)\}.$$

$$\mathbb{P}\left\{{\mathbb{E}}_{\mathbb{P}}\{c(\xi ,h)\}\le \underset{Q\in B_{ϵ}(\hat{\mathbb{P}})}{\sup }{\mathbb{E}}_Q\{c(\xi ,h)\}\right\}\ge 1-\beta .$$

$$\underset{Q\in B_{ϵ}(\hat{\mathbb{P}})}{\sup }{\mathbb{E}}_Q\left\{c(\xi ,h)\right\}$$

$$\le c(\hat{\xi },h)+ϵ\cdot \max \left\{\underset{\theta \in {\Theta }_2}{\sup }\Vert \theta {\Vert }_{\infty }\Vert \text{col}(g,0){\Vert }_{\ast },{\lambda }_{\text{ini}}\Vert h{\Vert }_{\ast }\right\},$$

$$\underset{h\in H}{\mathrm{argmin}}\{c(\hat{\xi },h)$$

$$+ϵ\cdot \max \left\{\underset{\theta \in {\Theta }_2}{\sup }\Vert \theta {\Vert }_{\infty }\Vert \text{col}(g,0){\Vert }_{\ast },{\lambda }_{\text{ini}}\Vert h{\Vert }_{\ast }\right\}\}.$$

$$\underset{\underset{\forall j\in \{1,\dots ,r\}}{{\zeta }_j\in {\Omega }_j}}{\sup }L({\zeta }_1^{T}b,\dots ,{\zeta }_r^{T}b)-\lambda \sum _{j=1}^r\Vert {\zeta }_j-{\hat{\zeta }}_j\Vert$$

$$\le \{\begin{array}{ll}L({\hat{\zeta }}_1^{T}b,\dots ,{\hat{\zeta }}_r^{T}b)& \text{if}{\sup }_{\theta \in \Theta }\Vert \theta {\Vert }_{\infty }\Vert b{\Vert }_{\ast }\le \lambda \\ \infty & \text{otherwise},\end{array}$$

$$L^{\ast }(z)$$

$$=\underset{s,\zeta \in \Omega }{\sup }\left\{\sum _{j=1}^r{z}_j{\zeta }_j^{T}b-L(s)|\underset{\forall j\in \{1,\dots ,r\}}{s_j={\zeta }_j^{T}b}\right\},$$

$$\mathcal{L}(s,\zeta ,\theta )=\sum _{j=1}^r\left(z_j{\zeta }_j^{T}b+{\theta }_j(s_j-{\zeta }_j^{T}b)\right)-L(s),$$

$$L^{\ast }(z)$$

Full text

I Introduction

Data-driven control approaches have become increasingly popular in recent years due to the growing complexity of systems [1]. The increase in system complexity poses significant challenges to model-based control design, in which acquiring an accurate model for the system is the most crucial step, and often the most time-consuming and expensive step [2, 3]. The goal of data-driven control approaches is to bypass the system identification step, and directly use system measurements for control design.

The focus of this paper is on the problem of data-driven control for unknown systems. We seek to develop a data-driven control algorithm that uses input/output data from the unknown system to compute control inputs which satisfy system constraints. Furthermore, since uncertainties and noise corruption in the data measurements are inevitable in any real-world application, we seek algorithms that are robust to the uncertain and noisy data measurements used.

Several approaches such as safe-learning, learning based Model Predictive Control (MPC), and stochastic MPC aim at solving similar problems (see, e.g., [4, 5, 6, 7] and references therein). However, these methods rely on having an a-priori known model for the system, require full-state measurement, and often assume knowledge about the disturbances that act on the system (e.g., bounded disturbances within a known set). Another large area of connected literature is reinforcement learning (see [8, 9] and references therein). These methods also address similar problems, but usually require a large number of data measurements. Additionally, these approaches sometimes result in unreliable outcomes due to the sensitivity of the algorithms to hyper-parameters [10]. Hence, none of the approaches above are suitable for real-time, constrained, and optimal control based only on input/output samples.

We employ an alternative approach which does not rely on a particular parametric system model or assume any knowledge of the disturbances. The approach uses behavioural system theory to characterize the behaviour of the unknown system (i.e., the possible trajectories of the system) [11]. This behavioural approach was used for control design in [12] where optimal open-loop control policies were computed using input/output data from the unknown system. This was later extended in [13] to a receding horizon set up, which was proved to have equivalent closed-loop performance when compared to standard MPC in the case of deterministic linear time invariant (LTI) systems.

The approach presented in this paper is built upon the Data-enabled Predictive Control (DeePC) algorithm for deterministic linear systems presented in [13]. The DeePC algorithm uses raw input/output data to compute optimal controls for unknown systems using real-time output feedback via a receding horizon implementation, thus allowing for the consideration of input/output constraints. In turn, this approach is much simpler to implement than the safe-learning, MPC, and reinforcement learning approaches above, which in the case of an unknown system with only output measurements available, require system identification and state observer design. On the other hand, the DeePC algorithm may not be amenable to stochastic systems, as the theory is built on the assumption of a deterministic LTI system. However, it was observed in [13] that after adding heuristic regularizations to the algorithm, DeePC still performed well on systems in the presence of stochastic disturbances, yet no robustness guarantees were given. We focus on extending this algorithm for use on stochastic systems. In doing so, we will give rigorous meaning to the heuristic regularizations presented in [13] by means of probabilistic guarantees on robust performance.

Due to the receding horizon implementation of the DeePC algorithm, we are required to repeatedly solve an optimization problem whose objective function naturally depends on data measurements that are affected by unknown disturbances. In order to be robust to these disturbances, we use distributionally robust optimization techniques [14, 15] which give rise to a regularized version of the DeePC algorithm similar to the heuristics presented in [13].

Contributions: Motivated by distributionally robust optimization and its connection to regularized optimization [14, 15], we develop a novel robust DeePC algorithm which uses input/output data to compute optimal and robust control inputs for unknown linear stochastic systems. We give probabilistic guarantees on its performance, giving rigorous justification for the use of the regularized control algorithm given in [13]. As a by-product, we gain new insights into the design of such regularizers.

Organization: Section II contains the problem statement. In Section III, we recall the DeePC algorithm and show how robust stochastic optimization techniques can be applied to improve the algorithm. Section IV contains the main results showing that regularizations on the DeePC algorithm result in robust performance. We illustrate the performance of the regularized algorithm in Section V. The paper is concluded in Section VI.

Notation: Given $x,y\in{}^{n}$, $\langle x,y\rangle\vcentcolon=x^{T}y$ denotes the usual inner product on ${}^{n}\times{}^{n}$. We denote the dual norm of a norm $\|\cdot\|$ on n by $\|x\|_{\ast}\vcentcolon=\sup_{\|y\|\leq 1}\langle x,y\rangle$. The conjugate function of a function $f:{}^{n}\to\real$ is denoted by $f^{\ast}(\theta)\vcentcolon=\sup_{x\in{}^{n}}\langle\theta,x\rangle-f(x)$. We denote the stacked column vector of $x$ and $y$ by ${\textup{col}}(x,y)=(x^{T},y^{T})^{T}$. We denote by $\delta_{x}$ the Dirac distribution at $x$. We use the notation $\widehat{\cdot}$ to denote objects that depend on measured data.

II Problem Statement

Consider the stochastic LTI system

[TABLE]

where $A\in{}^{n\times n}$, $B\in{}^{n\times m}$, $C\in{}^{p\times n}$, $D\in{}^{p\times m}$, $E\in{}^{n\times q}$, $F\in{}^{p\times q}$, and $x(t)\in{}^{n}$, $u(t)\in{}^{m}$, $y(t)\in{}^{p}$, $v(t)\in{}^{q}$ are respectively the state, control input, output, and stochastic disturbance of the system at time $t\in\mathbb{Z}_{\geq 0}$. The disturbance $v(t)$ is drawn from an unknown probability distribution $\mathbb{P}_{v}$ for all $t\in\mathbb{Z}_{\geq 0}$. We assume throughout the paper that the pair $(A,B)$ is controllable, the pair $(A,C)$ is observable, that system (1) itself is unknown (i.e., $A$, $B$, $C$, $D$, $E$, $F$ unknown), and that we have access only to input/output measurements, which we will denote by $\hat{u}(t),\hat{y}(t)$.

We consider a finite-horizon optimal control problem, in which the task is to design control inputs for the unknown system to minimize a given objective function. More specifically, given a time horizon ${T_{\textup{f}}}\in\mathbb{Z}_{>0}$, an objective function $f\colon{}^{m{T_{\textup{f}}}}\times{}^{p{T_{\textup{f}}}}\to\real$, and the current time $t\in\mathbb{Z}_{\geq 0}$, we wish to choose a sequence of control inputs $u={\textup{col}}(u(t),\dots,u(t+{T_{\textup{f}}}-1))\in\mathcal{U}$, where $\mathcal{U}\subset{}^{m{T_{\textup{f}}}}$ is an input constraint set, such that the resulting stochastic output trajectory, $y={\textup{col}}(y(t),\dots,y(t+{T_{\textup{f}}}-1))$, of (1) minimizes the expectation of the objective function, i.e., we wish to solve the following optimization problem:

[TABLE]

where $\mathbb{P}_{v}^{{T_{\textup{f}}}}=\mathbb{P}_{v}\times\cdots\times\mathbb{P}_{v}$ is the ${T_{\textup{f}}}$-fold product distribution. Such optimization problems also appear in stochastic MPC [16].

Solving problem (2) poses two main challenges: the system (1) is unknown, and the disturbance distribution $\mathbb{P}_{v}$ is unknown. Hence, we are unable to predict future output trajectories $y$, and we are missing necessary information to compute the expectation. Note that even in the case when system (1) and $\mathbb{P}_{v}$ are known, solving (2) would require high-dimensional integration and is often computationally intractable [17]. One could also consider including joint output chance constraints in the problem setup (2) which would pose similar challenges as above. We do not consider such constraints as they are beyond the scope of this paper.

To address these challenges, we first simplify the problem by considering the deterministic version of system (1):

[TABLE]

In this case, problem (2) reduces to

[TABLE]

We recall a data-enabled predictive control (DeePC) algorithm first presented in [13] that can be used to solve (4) for deterministic though unknown systems (3) (see Section III-A). We then show how this algorithm can be extended to approach problem (2) for stochastic systems (1) using distributionally robust optimization techniques (see Section III-B). The main results can be found in Section IV.

III Preliminaries

III-A DeePC Algorithm for Deterministic Systems

In this section, we recall the DeePC algorithm which uses raw input/output data to construct a non-parametric predictive model that was first developed using a behavioural system theory approach in [11]. Let $L,q,T\in\mathbb{Z}_{>0}$. We define the Hankel matrix of a signal $w\in{}^{qT}$ with $L$ block rows as the matrix

[TABLE]

We say that signal $w$ is persistently exciting of order $L$ if $\mathscr{H}_{L}(w)$ has full row rank. Note that in order for $w$ to be persistently exciting of order $L$, we must have $T\geq(q+1)L-1$, i.e., the signal $w$ has to be sufficiently rich and sufficiently long. The lag of system (3) is defined as the smallest integer $\ell\in\mathbb{Z}_{>0}$ such that the observability matrix $\mathscr{O}_{\ell}(A,C)\vcentcolon={\textup{col}}\left(C,CA,\dots,CA^{\ell-1}\right)$ has rank $n$ (see [18, Section 7.2] for equivalent state space free definitions of lag).

Like in Section II, we assume that we have access to input/output measurements. Let ${\textup{col}}(\hat{u},\hat{y})={\textup{col}}(\hat{u}(1),\dots,\hat{u}(T),\hat{y}(1),\dots,\hat{y}(T))\in{}^{(m+p)T}$ be a measured trajectory of (3) of length $T\in\mathbb{Z}_{>0}$. Assume that $\hat{u}$ is persistently exciting of order ${T_{\textup{ini}}}+{T_{\textup{f}}}+n$, where ${T_{\textup{ini}}},{T_{\textup{f}}}\in\mathbb{Z}_{>0}$. We organize the data into the Hankel matrices

[TABLE]

where ${\widehat{U}_{\mathrm{p}}}$ consists of the first ${T_{\textup{ini}}}$ block rows of $\mathscr{H}_{{T_{\textup{ini}}}+{T_{\textup{f}}}}(\hat{u})$ and ${\widehat{U}_{\mathrm{f}}}$ consists of the last ${T_{\textup{f}}}$ block rows of $\mathscr{H}_{{T_{\textup{ini}}}+{T_{\textup{f}}}}(\hat{u})$ (similarly for ${\widehat{Y}_{\mathrm{p}}}$ and ${\widehat{Y}_{\mathrm{f}}}$). Given the current time $t\in\mathbb{Z}_{\geq 0}$, let $\hat{u}_{\textup{ini}}={\textup{col}}(\hat{u}(t-{T_{\textup{ini}}}),\dots,\hat{u}(t-1))\in{}^{m{T_{\textup{ini}}}}$ and $\hat{y}_{\textup{ini}}={\textup{col}}(\hat{y}(t-{T_{\textup{ini}}}),\dots,\hat{y}(t-1))\in{}^{p{T_{\textup{ini}}}}$ be the ${T_{\textup{ini}}}$ most recent input and output measurements of (3), respectively. By [11, Theorem 1], ${\textup{col}}(\hat{u}_{\textup{ini}},u,\hat{y}_{\textup{ini}},y)\in{}^{(m+p)({T_{\textup{ini}}}+{T_{\textup{f}}})}$ is a trajectory of system (3) if and only if there exists $g\in{}^{T-{T_{\textup{ini}}}-{T_{\textup{f}}}+1}$ such that

[TABLE]

Furthermore, if ${T_{\textup{ini}}}\geq\ell$, then by [12, Lemma 1], there exists a unique $x(t)\in{}^{n}$ such that the output trajectory $y$ is uniquely determined by the system (3). In other words, the trajectory ${\textup{col}}(\hat{u}_{\textup{ini}},\hat{y}_{\textup{ini}})$ fixes the underlying initial state $x(t)$ from which the trajectory ${\textup{col}}(u,y)$ evolves. Hence, the Hankel matrix in (6) serves as a non-parametric predictive model for system (3). This observation was first exploited for control purposes in [12], where equation (6) was used to construct open-loop control policies for tracking. This was then extended in [13], where equation (6) was used in a receding horizon optimal control algorithm. We recall the latter approach below.

Given a time horizon ${T_{\textup{f}}}\in\mathbb{Z}_{>0}$, past input/output data ${\textup{col}}(\hat{u}_{\textup{ini}},\hat{y}_{\textup{ini}})\in{}^{(m+p){T_{\textup{ini}}}}$, objective function $f\colon{}^{m{T_{\textup{f}}}}\times{}^{p{T_{\textup{f}}}}\to\real$, we formulate the following optimization problem:

[TABLE]

Note that optimization problem (III-A) is equivalent to (4) since by (6), $u={\widehat{U}_{\mathrm{f}}}g$ and $y={\widehat{Y}_{\mathrm{f}}}g$. The data-enabled predictive control (DeePC) algorithm is then given as follows:

This algorithm was shown to be equivalent to the classical Model Predictive Control (MPC) algorithm when considering deterministic LTI systems of the form (3), also in the presence of additional output constraints. When Algorithm 1 was applied to a nonlinear stochastic quadcopter model in [13], additional heuristic regularization terms had to be included in (III-A) to achieve robust performance. We recall the regularizations below.

[TABLE]

where $\lambda_{\textup{ini}},\lambda_{g}\in{\mathbb{R}}_{\geq 0}$. We will show that this regularized optimization problem coincides with solving a distributionally robust variation of optimization problem (2), giving rigorous justification for the regularizations (see Section IV).

III-B Distributionally Robust DeePC

Consider now the stochastic system (1). Let ${\textup{col}}(\hat{u},\hat{y})\in{}^{(m+p)T}$ be a measured trajectory of length $T$ of system (1) such that $\hat{u}$ is persistently exciting of order ${T_{\textup{ini}}}+{T_{\textup{f}}}+n$, where ${T_{\textup{ini}}}\geq\ell$. Furthermore, assume ${\textup{col}}(\hat{u}_{\textup{ini}},\hat{y}_{\textup{ini}})\in{}^{(m+p){T_{\textup{ini}}}}$ is the last measured trajectory of system (1). The output trajectories $\hat{y}$ and $\hat{y}_{\textup{ini}}$ can be viewed as particular realizations of random variables, $y$ and ${y_{\textup{ini}}}$. If we naively use the noise corrupted trajectory $\hat{y}$ to build Hankel matrices ${\widehat{U}_{\mathrm{p}}}$, ${\widehat{Y}_{\mathrm{p}}}$, ${\widehat{U}_{\mathrm{f}}}$, ${\widehat{Y}_{\mathrm{f}}}$ as in (5) and apply Algorithm 1, we will run into difficulties regarding the consistency of the constraint equations in (III-A). Indeed, since $\hat{y}_{\textup{ini}}$ and ${\widehat{Y}_{\mathrm{p}}}$ consist of noise corrupted outputs, there may not exist $g$ that satisfies the equation ${\widehat{Y}_{\mathrm{p}}}g=\hat{y}_{\textup{ini}}$. Hence, we soften the equality constraint and penalize the slack variable with an appropriate cost function; we recall that the use of such “soft constraints” is common in MPC [19]. This results in the optimization problem

[TABLE]

where we have used the 1-norm as a penalty function on the slack variable. It is well known that if $\lambda_{\textup{ini}}$ is chosen large enough then the solution $g^{\star}$ will violate ${\widehat{Y}_{\mathrm{p}}}g^{\star}=\hat{y}_{\textup{ini}}$ only if the equation is infeasible [20]. Hence, the original constraint will be satisfied if it can be satisfied. Since the constraint is now deterministic, we may define $G=\{g\in{}^{T-{T_{\textup{ini}}}-{T_{\textup{f}}}+1}\mid{\widehat{U}_{\mathrm{p}}}g=\hat{u}_{\textup{ini}},\;{\widehat{U}_{\mathrm{f}}}g\in\mathcal{U}\}$ and rewrite the above as

[TABLE]

For ease of notation, we put all random objects into a matrix whose $j$-th row we denote by $\xi_{j}$, i.e.,

[TABLE]

where $\begin{pmatrix}\cdot\end{pmatrix}_{j,\cdot}$ denotes the $j$-th row of a matrix. For all $j\in\{1,\dots,p({T_{\textup{ini}}}+{T_{\textup{f}}})\}$, let $\mathbb{P}_{j}$ denote the probability distribution of $\xi_{j}$ supported on $\Xi_{j}\subseteq{}^{T-{T_{\textup{ini}}}-{T_{\textup{f}}}+2}$. Note that the distributions $\mathbb{P}_{j}$ and their support sets $\Xi_{j}$ are determined by the dynamics of system (1) and the unknown distribution $\mathbb{P}_{v}$. Define $\xi={\textup{col}}(\xi_{1}^{T},\dots,\xi_{p({T_{\textup{ini}}}+{T_{\textup{f}}})}^{T})$ and let $\mathbb{P}\vcentcolon=\mathbb{P}_{1}\times\cdots\times\mathbb{P}_{p({T_{\textup{ini}}}+{T_{\textup{f}}})}$ denote the unknown probability distribution of $\xi$ supported on $\Xi=\prod_{j=1}^{p({T_{\textup{ini}}}+{T_{\textup{f}}})}\Xi_{j}\subseteq{}^{p({T_{\textup{ini}}}+{T_{\textup{f}}})(T-{T_{\textup{ini}}}-{T_{\textup{f}}}+2)}$. Let $H=\{h={\textup{col}}(g,-1)\mid{\widehat{U}_{\mathrm{p}}}g=\hat{u}_{\textup{ini}},\;{\widehat{U}_{\mathrm{f}}}g\in\mathcal{U}\}$. Substituting this notation into (9) yields

[TABLE]

where $\hat{\xi}={\textup{col}}(\hat{\xi}_{1},\dots,\hat{\xi}_{p({T_{\textup{ini}}}+{T_{\textup{f}}})})$ are our measurements. Denoting the empirical distribution of our measurements $\hat{\xi}$ by $\widehat{\mathbb{P}}=\delta_{\hat{\xi}}$, and substituting this notation into the above yields

[TABLE]

where $c(\xi,h)=f(({\widehat{U}_{\mathrm{f}}},0)h,(\xi_{p{T_{\textup{ini}}}+1}h,\dots,\xi_{p({T_{\textup{ini}}}+{T_{\textup{f}}})}h))+\lambda_{\textup{ini}}\|(\xi_{1}h,\dots,\xi_{p{T_{\textup{ini}}}}h)\|_{1}$.

The quantity in (10) is known as the in-sample performance. Unfortunately, the solution $h^{\star}$ to (10) may result in poor out-of-sample performance

[TABLE]

which if we compare to equation (2) is the real quantity of interest (see, e.g., [14, 21] for examples displaying poor out-of-sample performance).

To alleviate this problem, we focus instead on a robust variation of (10) which will serve as an upper bound for out-of-sample performance (11) with high confidence. In particular, we seek solutions of

[TABLE]

where $\widehat{\mathcal{P}}$ is an ambiguity set which depends on the sampled data trajectories $\hat{\xi}$. The ambiguity set will be constructed in such a way that it contains the actual distribution $\mathbb{P}$ with high confidence. Thus, if $h^{\star}$ and $J^{\star}$ are the solution and optimal value of (12), then $J^{\star}$ will upper bound the out-of-sample performance with high confidence.

Following [14], we define the ambiguity set in (12) in terms of the Wasserstein metric defined on the space $\mathcal{M}(\Xi)$ denoting the set of all distributions $Q$ supported on $\Xi$ such that $\mathbb{E}_{Q}[\|\xi\|_{\textup{W}}]<\infty$, where $\|\cdot\|_{\textup{W}}$ is an arbitrary norm.

Definition III.1

The Wasserstein metric $\textup{d}_{\textup{W}}\colon\mathcal{M}(\Xi)\times\mathcal{M}(\Xi)\to{\mathbb{R}}_{\geq 0}$ is defined as

[TABLE]

where $\Pi$ is a joint distribution of $\xi_{1}$ and $\xi_{2}$ with marginal distributions $Q_{1}\in\mathcal{M}(\Xi)$ and $Q_{2}\in\mathcal{M}(\Xi)$ respectively.

Let $\epsilon\geq 0$. We denote the Wasserstein ball of radius $\epsilon$ centred around distribution $Q$ by $B_{\epsilon}(Q)\vcentcolon=\{Q^{\prime}\in\mathcal{M}(\Xi)\mid\textup{d}_{\textup{W}}(Q,Q^{\prime})\leq\epsilon\}$. The Wasserstein metric can be viewed as a distance between probability distributions, where the distance is calculated via an optimal mass transport plan $\Pi$. Note that there are other ways to construct ambiguity sets (see, e.g., [21] where the ambiguity set is constructed as the confidence region of a goodness-of-fit hypothesis test). Replacing the general ambiguity set $\widehat{\mathcal{P}}$ with a Wasserstein ball around the empirical distribution $\widehat{\mathbb{P}}$ results in the problem

[TABLE]

In the next section, we show that if $h^{\star}$ and $J^{\star}$ are the solution and optimal value of (13), then $J^{\star}$ will upper bound the out-of-sample performance (11) with high confidence. Hence, we obtain probabilistic guarantees that applying control inputs $u^{\star}={\widehat{U}_{\mathrm{f}}}h^{\star}$ to system (1) will result in good performance of the resulting stochastic trajectory $y$. Additionally, we show that (13) is computationally tractable and results in a regularized version of the DeePC algorithm similar to (III-A).

IV Main Results

The following result relates the robust problem (13) to the out-of-sample performance (11). In particular, if $h^{\star}$ and $J^{\star}$ are the solution and the optimal value of the robust problem (13), then $J^{\star}$ upper bounds the out-of-sample performance (11) with high confidence.

Theorem IV.1

(Robust Performance Guarantee):* Assume that distribution $\mathbb{P}$ is light-tailed, i.e., there exists $a>1$ such that $\mathbb{E}_{\mathbb{P}}[e^{\|\xi\|_{\textup{W}}^{a}}]<\infty$. Let $\beta\in(0,1)$. Then there exists $\epsilon=\epsilon(\beta)>0$ such that for all $h\in H$,*

[TABLE]

The proof of the above theorem follows directly from [14, Theorem 3.5]. Note that the light-tailed assumption is satisfied automatically when $\Xi$ is compact. Hence, all distributions truncated to a compact support set satisfy the assumptions of Theorem IV.1. Other examples include Gaussian and exponential distributions.

By adapting the proof methods of [14] and [15] to our setting, we show that for certain objective functions, the semi-infinite optimization problem (13) reduces to a tractable convex program that coincides with a regularized version of (10). Let $\|\cdot\|$ be an arbitrary norm on ${}^{T-{T_{\textup{ini}}}-{T_{\textup{f}}}+2}$ and $\|\xi\|_{\textup{W}}=\sum_{j=1}^{p({T_{\textup{ini}}}+{T_{\textup{f}}})}\|\xi_{j}\|$ be the norm used in the Wasserstein metric.

Theorem IV.2

(Tractable Reformulation):* Assume that the objective function $f$ is separable and can be written as $f(u,y)=f_{1}(u)+f_{2}(y)$ for all $u,y$, where $f_{1}$ and $f_{2}$ are convex and continuous. Furthermore, assume $f_{2}$ is such that $\Theta_{2}=\{\theta\mid f_{2}^{\ast}(\theta)<\infty\}$ is a bounded set in ${}^{p{T_{\textup{f}}}}$. Then*

[TABLE]

where $\|\cdot\|_{\ast}$ denotes the dual norm of $\|\cdot\|$. Equality holds when $\Xi_{j}={}^{T-{T_{\textup{ini}}}-{T_{\textup{f}}}+2}$ for all $j\in\{1,\dots,p{T_{\textup{ini}}}\}$ and $\Xi_{j}={}^{T-{T_{\textup{ini}}}-{T_{\textup{f}}}+1}\times\{0\}$ for all $j\in\{p{T_{\textup{ini}}}+1,\dots,p({T_{\textup{ini}}}+{T_{\textup{f}}})\}$.

Note that the condition for equality holds when the disturbance $v$ affecting system (1) is drawn from a distribution $\mathbb{P}_{v}$ defined on an unbounded support set and the matrix $F$ in (1) is full row rank. In other words, equality will hold if the probability distribution $\mathbb{P}_{v}$ has unbounded support and disturbances $v(t)$ affect all entries of the output vector $y(t)$. In this case, each probability distribution $\mathbb{P}_{j}$ would have support $\Xi_{j}={}^{T-{T_{\textup{ini}}}-{T_{\textup{f}}}+2}$ for $j\in\{1,\dots,p{T_{\textup{ini}}}\}$ and $\Xi_{j}={}^{T-{T_{\textup{ini}}}-{T_{\textup{f}}}+1}\times\{0\}$ for $j\in\{p{T_{\textup{ini}}}+1,\dots,p({T_{\textup{ini}}}+{T_{\textup{f}}})\}$. For example, we would obtain equality in the above if $F$ is full row rank and $\mathbb{P}_{v}$ is a Gaussian distribution which is a common assumption in stochastic MPC (see [7] and references therein). We also see an immediate connection between Theorem IV.2 and (III-A). In fact, when $\|\cdot\|_{\ast}=\|\cdot\|_{1}$ then $\epsilon\cdot\max\{\sup_{\theta\in\Theta_{2}}\|\theta\|_{\infty}\|{\textup{col}}(g,0)\|_{\ast},\lambda_{\textup{ini}}\|h\|_{\ast}\}=\epsilon\cdot\max\{\sup_{\theta\in\Theta_{2}}\|\theta\|_{\infty}\|g\|_{1},\lambda_{\textup{ini}}(\|g\|_{1}+1)\}$. Hence, depending on the known cost function $f_{2}$ (i.e., on the set $\Theta_{2}$), $\lambda_{g}$ in (III-A) plays the role of either $\epsilon\cdot\sup_{\theta\in\Theta_{2}}\|\theta\|_{\infty}$ or $\epsilon\lambda_{\textup{ini}}$. Thus, the result of Theorem IV.2 gives insight into the design of the regularizer on the decision variable $g$ to achieve various robustness goals for DeePC. The results above give rise to the robust DeePC algorithm.

Remark IV.1

The Wasserstein ball radius $\epsilon(\beta)$ in Theorem IV.1 is often larger than necessary, i.e., $\mathbb{P}\not\in B_{\epsilon}(\widehat{\mathbb{P}})$ with probability much less than $\beta$. Furthermore, even when $\mathbb{P}\not\in B_{\epsilon}(\widehat{\mathbb{P}})$, the robust quantity $\sup_{Q\in B_{\epsilon}(\widehat{\mathbb{P}})}\mathbb{E}_{Q}\{c(\xi,h)\}$ may still serve as an upper bound for the out-of-sample performance $\mathbb{E}_{\mathbb{P}}\{c(\xi,h)\}$ [14]. Thus, for practical purposes, one should choose the radius $\epsilon$ of the Wasserstein ball in a data-driven fashion (see Section V). $\bullet$

Remark IV.2

The assumption that $\Theta_{2}=\{\theta\mid f_{2}^{\ast}(\theta)<\infty\}$ is a bounded set in ${}^{p{T_{\textup{f}}}}$ holds for many objective functions of practical interest. Indeed, any arbitrary norm satisfies this assumption, and in fact any Lipschitz continuous function. $\bullet$

Remark IV.3

The norm $\sum_{j=1}^{p({T_{\textup{ini}}}+{T_{\textup{f}}})}\|\xi_{j}\|$ used in the Wasserstein metric can be chosen to achieve various robustness goals. Indeed, being robust in the trajectory space in the $\|\cdot\|_{\infty}$ sense requires regularizing with a $\|\cdot\|_{1}$. Likewise, being robust in the trajectory space in the $\|\cdot\|_{2}$ sense requires regularizing with a $\|\cdot\|_{2}$. If no noise is present in the system, we may set $\epsilon=0$ recovering the DeePC algorithm for deterministic LTI systems. $\bullet$

To prove the theorem above, we require the following lemma which is an extension of [15, Lemma A.3] to functions with vector valued inputs:

Lemma IV.1

Let $q,r\in\mathbb{Z}_{>0}$, and $\Omega_{j}\subseteq{}^{q}$ for $j\in\{1,\dots,r\}$. Let $\hat{\zeta}_{j}\in\Omega_{j}$ be given for $j\in\{1,\dots,r\}$, $\lambda\in{\mathbb{R}}_{>0}$ and $L\colon{}^{r}\to\real$ convex and continuous such that $\{\theta\mid L^{\ast}(\theta)<\infty\}$ is a bounded set in r. Then, for fixed $b\in{}^{q}$,

[TABLE]

where $\|\cdot\|$ is an arbitrary norm on q. Furthermore, the above is an equality when $\Omega_{j}={}^{q}$ for all $j\in\{1,\dots,r\}$.

Proof:

By definition of the conjugate function,

[TABLE]

where $\zeta={\textup{col}}(\zeta_{1},\dots,\zeta_{r})$, $z={\textup{col}}(z_{1},\dots,z_{r})$, $s={\textup{col}}(s_{1},\dots,s_{r})$, and $\Omega=\Omega_{1}\times\cdots\times\Omega_{r}$. The Lagrangian of the above is given by

[TABLE]

where $\theta={\textup{col}}(\theta_{1},\dots,\theta_{r})$. By strong duality (see, e.g., [22, Proposition 5.3.1]),

[TABLE]

Hence, by duality

[TABLE]

where $\Theta=\{\theta\mid L^{\ast}(\theta)<\infty\}$ is the effective domain of $L^{\ast}$. Since $L$ is convex and continuous, the biconjugate $L^{\ast\ast}$ coincides with the function $L$ itself. Hence,

[TABLE]

Thus,

[TABLE]

where the last equality comes from the definition of the dual norm and using homogeneity of the norm. Using the minimax theorem (see, e.g., [22, Proposition 5.5.4]) we switch the supremum and infimum in the above and bring it into the sum giving

[TABLE]

where equality holds when $\Omega_{j}={}^{q}$ for all $j\in\{1,\dots,r\}$. Hence,

[TABLE]

This proves the claimed result. ∎

Remark IV.4

Lemma IV.1 may also be used to extend the results of [15] to learning problems in which the input data is matrix-valued. We do not explore this connection as this is beyond the scope of this paper. $\bullet$

Proof:

(Theorem IV.2): Since $c(\cdot,h)$ is a proper, convex, continuous function for all $h\in H$ then

[TABLE]

This can be shown by studying the dual problem, and noticing that the worst case distribution coincides with a Dirac distribution at the point $\xi$ which results in the largest norm $\|\xi-\hat{\xi}\|_{\textup{W}}$ (see [14, Theorem 6.3]). Define $f_{3}(\cdot)\vcentcolon=\lambda_{\textup{ini}}\|\cdot\|_{1}$. By separability of the objective function,

[TABLE]

where we used the epigraph formulation and the definition of $\|\cdot\|_{\textup{W}}$ for the last equality. By definition of the conjugate function, $f_{3}^{\ast}(\theta)=0$ if $\|\theta\|_{\infty}\leq\lambda_{\textup{ini}}$ and is infinite otherwise. Hence, $\Theta_{3}\vcentcolon=\{\theta\mid f_{3}^{\ast}(\theta)<\infty\}$ is a bounded set in ${}^{p{T_{\textup{f}}}}$. We also have that $\Theta_{2}\vcentcolon=\{\theta\mid f_{2}^{\ast}(\theta)<\infty\}$ by assumption. Hence, by Lemma IV.1,

[TABLE]

with equality when for all $j\in\{1,\dots,p{T_{\textup{ini}}}\}$, $\Xi_{j}={}^{T-{T_{\textup{ini}}}-{T_{\textup{f}}}+2}$. Note that the uncertainties $\xi_{p{T_{\textup{ini}}}+1},\dots,\xi_{p({T_{\textup{ini}}}+{T_{\textup{f}}})}$ affecting $f_{2}$ have 0 as their last entry by definition. Hence, carrying through the steps of the proof of Lemma IV.1 for $f_{2}$ yields

[TABLE]

where $\tilde{h}={\textup{col}}(g,0)$ and equality holds when $\Xi_{j}={}^{T-{T_{\textup{ini}}}-{T_{\textup{f}}}+1}\times\{0\}$ for all $j\in\{p{T_{\textup{ini}}}+1,\dots,p({T_{\textup{ini}}}+{T_{\textup{f}}})\}$. By reversing the epigraph formulation, we have

[TABLE]

Note that the infinite case is dropped since we are taking the infimum over $\lambda\geq 0$. Substituting notation yields

[TABLE]

Noting that $\tilde{h}={\textup{col}}(g,0)$ gives the claimed result. ∎

V Simulations

We illustrate the performance of the Robust DeePC Algorithm 2 on a model of a quadcopter linearized around the hover position. The states of the quadcopter model are given by the 3 spatial coordinates ($x$, $y$, $z$) and their velocities, and the 3 angular coordinates $(\alpha,\beta,\gamma)$ and their velocities, i.e., the state is $(x,y,z,\dot{x},\dot{y},\dot{z},\alpha,\beta,\gamma,\dot{\alpha},\dot{\beta},\dot{\gamma})$. The inputs are given by the individual thrusts from the 4 rotors. Full state measurement was assumed. The states are affected by additive zero-mean Gaussian noise. The state-space matrices used in (1) are

[TABLE]

[TABLE]

[TABLE]

During the data collection process, the input was drawn from a uniform random variable to ensure persistency of excitation. The resulting output data was corrupted by additive noise drawn from a zero-mean Gaussian distribution. We used 214 input/output measurements to populate the matrices ${\widehat{U}_{\mathrm{p}}}$, ${\widehat{Y}_{\mathrm{p}}}$, ${\widehat{U}_{\mathrm{f}}}$, and ${\widehat{Y}_{\mathrm{f}}}$, which is the minimum number of measurements needed to ensure persistency of excitation. We commanded the quadcopter to track a parameterized figure-8 reference trajectory, denoted by $r$.

The following parameters were chosen for the optimization problem (13): the infinity-norm $\|\cdot\|_{\infty}$ was used in the definition of the the norm $\|\cdot\|_{\textup{W}}$ used in the Wasserstein metric, ${T_{\textup{f}}}=30$, ${T_{\textup{ini}}}=1$, $f(u,y)=\|u\|_{1}+200\|y-r\|_{1}$, $\lambda_{\textup{ini}}=10^{5}$, $\epsilon=0.001$, $\mathcal{U}=[-0.7007,0.2993]^{{T_{\textup{f}}}}$. The input constraint set $\mathcal{U}$ is chosen this way to mimic the constraints present in a nonlinear quadcopter model where the normalized rotor thrusts can only vary in the set $[0,1]$. Online measurement noise was drawn from same Gaussian distribution. The performance of the algorithm is seen in Figures 1-3. As can be seen, the algorithm exhibits desirable behaviour despite the measurement noise and process noise.

As discussed in Remark IV.1, computing an optimal radius $\epsilon$ for the Wasserstein ball in a data-driven fashion may increase the performance of the algorithm. We present results of the robust algorithm for various Wasserstein radii when the quadcopter was commanded to follow a step trajectory. By observing the performance of the quadcopter, we are able to estimate a range of optimal Wasserstein radii. We again measured the full state. The states were affected by zero-mean truncated Gaussian noise, where the truncation was made 3 standard deviations from both sides of the mean. The output measurements were also corrupted by a zero-mean truncated Gaussian distribution. The infinity-norm was used in the definition of the Wasserstein metric. A Wasserstein radius $\epsilon$ was fixed and the cost $f(u,y)$ accumulated over the 20 second horizon of the simulation was computed. The cost was averaged over 15 random simulations for each $\epsilon$.

The results shown in Figure 4 indicate that $\epsilon$ should be chosen in the interval $[0.001,0.01]$. The results also support the claim that optimizing the in-sample performance (10) (i.e., setting $\epsilon=0$) displays poor out-of-sample performance, hence justifying the robust approach. It is also clear from Figure 4 that choosing $\epsilon$ large (i.e., being over-conservative) results in poor performance.

VI Conclusion

We studied the problem of controlling an unknown stochastic system with the goal of minimizing an objective function dependent on the input/output trajectories of the system. With no knowledge of the nature of the stochasticity, we proposed a robust DeePC algorithm which uses noise corrupted input/output data to compute optimal and robust control inputs. Robustifying the original DeePC algorithm gave rise to principled regularization terms, supporting the observed superior performance of the regularized algorithm on stochastic systems. Future work includes incorporating multiple measured data sets to improve performance, and including output constraints.