Data–driven predictive control with estimated prediction matrices and integral action

P. C. N. Verheijen, G. R. Gonçalves da Silva and M. Lazar The authors are with the Department of Electrical Engineering, Eindhoven University of Technology, The Netherlands: p.c.n.verheijen@student.tue.nl, g.goncalves.da.silva@tue.nl, m.lazar@tue.nl

Abstract

This paper presents a data–driven approach to the design of predictive controllers. The prediction matrices utilized in standard model predictive control (MPC) algorithms are typically constructed using knowledge of a system model such as, state–space or input–output models. Instead, we directly estimate the prediction matrices relating future outputs with current and future inputs from measured data, off–line. On–line, the developed data–driven predictive controller reduces to solving a quadratic program with a similar structure and complexity as linear MPC. Additionally, we develop a new procedure for estimating prediction matrices from data for predictive controllers with integral action, corresponding to the rate–based formulation of linear MPC. The effectiveness of the developed data–driven predictive controller is illustrated on position control of a linear motor model.

Index Terms:

Data–driven control, Model predictive control, Markov parameters estimation, Integral control

I Introduction

Data–driven (DD) control design methods have recently gained an increased attention from the control systems community due to the digital technology trends involving big–data and artificial intelligence systems, see, for example, [1, 2] and the references therein. The main idea of data–driven control is to eliminate the standard separation in model–based controller design, i.e., first obtain a system model, by identification or first principles, and then design a controller. Instead, in DD controller design, more freedom is allowed in mixing identification (or estimation) with controller design, and even direct controller synthesis from data.

Much of the existing data–driven methods for controller design, see, e.g., [3, 4, 5, 6] and the references therein, have originated within the field of adaptive control and make use of the so–called model reference control framework. In these approaches, instead of first deriving a system model, measured input–output (I/O) data is directly used to identify a parameterized controller chosen a priori, given a performance requirement, i.e., the reference model. Since it is hard to guarantee closed–loop stability while tuning the controller with these methods without largely compromising performance, an a posteriori data–driven certification procedure has also been reported in [7].

Recently, these ideas have been extended to other types of robust control problems, such as state–feedback stabilization [8], linear quadratic regulator (LQR) design [9, 10], and output feedback with an embedded prediction step [11]. Combining ideas of state–feedback with prediction theory fits well in the framework of model predictive control (MPC) [12], which can also be described by a given performance (whether a reference model or weighting matrices) with a parameterized structure (the prediction matrices, typically constructed using a system model).

MPC approaches based on input–output step–response system models [13] or finite impulse response (FIR) system models [14], have been at the core of the MPC research for a long time. These methods require an identification step to obtain the FIR model (usually by applying some impulses to the system or by differentiating its step response) and an estimation of the system order, and then plugging in these parameters into the prediction matrices. Such an approach, which still resorts to a system identification procedure, is referred to as an indirect data–driven approach to controller design. Another possibility is to combine model–based and data–driven approaches within predictive control algorithms, see, for example, [15] and the references therein, where a linear state–space model is used for describing known system dynamics and a learning data–driven approach is used to model unknown disturbances.

Direct approaches to data–driven predictive control, inspired by [16], were proposed recently in [17] and [18]. These approaches parameterize future predicted outputs as a linear combination of past inputs and outputs, and future inputs, via Hankel matrices of I/O data and a set of coefficients (parameters). Such methods do not require any off–line identification, at the cost of computing the set of coefficients on–line, simultaneously with computing future control inputs and predicted outputs. This makes the number of optimization variables dependent of both the prediction horizon and the size of the Hankel data matrices. It is worth to mention that [17] computes on–line an additive control input that adjusts an unconstrained predictive control law, while [18] computes on–line the constrained predictive control input directly, as in standard MPC. Also, [17] proposes a method to include integral action, while [18] develops solutions for dealing with noisy data. Stability and robustness of Hankel matrices based predictive controllers have been recently studied in [19].

Overall, direct data–driven predictive controllers are sensitive to noisy data, since they do not inherit any of the standard properties of an intermediate step of estimation (and also due to multiplicative noise terms in the prediction step). Another relevant remark to MPC design [12] versus data–driven predictive control (DPC) design is that stability and optimal performance is not anymore easily determined by the cost function tuning. Indeed, while in linear MPC the LQR control law can be recovered if constraints are not active via a suitable terminal penalty, this is no longer the case in DPC, even under the assumption that the full state is measurable.

An unconstrained solution to an input-output data-driven Subspace Predictive Control algorithm has been derived in [20]. However, the authors suggest to compute the prediction matrices by using a pre-step of QR decomposition of the data matrices, and it is not clear how this computation is affected by noise. They then propose a second step to reduce the state prediction matrix to an estimate of the system order. [21] builds on the input-output GPC approach, including only upper and lower bound type of constraints and proposing an integral action assuming a particular ARIMAX model structure, so they can take the discrete output difference as feedback. Therein, integral prediction matrices are not directly identified from data, instead, they identify the non-integral prediction matrices first, and then they use summation of the specific terms. The method developed in this work do not require any type of system order estimate and we formulate it such that inclusion of constraints follow the model-based approach. We also derive the identification procedure for identifying prediction matrices for rate-based integral action state-space MPC with formal guarantees (Theorem 1) without assuming a specific noise model.

Motivated by the current status in DPC design, in this paper we develop an approach to data–driven predictive control that offers an attractive compromise between indirect and direct approaches to data–driven control. The approach is inspired by the recent data–driven LQR design proposed in [10]. Firstly, instead of identifying a state–space or FIR system model and then building the MPC prediction matrices, the developed approach directly estimates the full prediction matrices from measured data. Estimation of the prediction matrices makes the developed approach less sensitive to noisy data. We also show that these matrices are comprised of the system Markov parameters, so they relate to the MPC approaches using FIR models, but without requiring a priori knowledge of the system order. Moreover, we show that if the state sequence can be measured or the prediction horizon is long enough, then one can recover the corresponding LQR controller by using the derived DPC algorithm. Motivated by the need of off-set free control in practice, another contribution of this paper is the derivation of a rate–based DPC algorithm with integral action, similar to linear MPC with integral action [22]. The performance and robustness of the developed data–driven predictive controller is illustrated for position control of a linear motor model.

II Preliminaries

In this section we recall the standard linear MPC formulation [12] using prediction matrices based on linear discrete–time state–space models, i.e.,

\begin{aligned} x(k+1)&=Ax(k)+Bu(k)\\ y(k)&=Cx(k)\end{aligned},\quad k\in\mathbb{N},

(1)

where $x(k)\in\mathbb{R}^{n}$ is the state, $u(k)\in\mathbb{R}^{m}$ is the input, $y(k)\in\mathbb{R}^{q}$ is the output and $k$ denotes the discrete–time index. Assuming for simplicity that the regulation objective is to control the output to the origin, in MPC one computes a sequence of optimal control inputs $U_{k}^{\ast}:=U^{\ast}(x(k))$ every time instant $k$ by minimizing a cost function, e.g.,

$\displaystyle J(x(k),U_{k})$	$\displaystyle:=y_{N\|k}^{T}Py_{N\|k}+\sum_{i=0}^{N-1}(y_{i\|k}^{T}Qy_{i\|k}+u_{i\|k}^{T}Ru_{i\|k})$	(2)
	$\displaystyle=U_{k}^{T}(\Psi+\Gamma^{T}\Omega\Gamma)U_{k}+2U_{k}^{T}\Gamma^{T}\Omega\Phi x(k)+$
	$\displaystyle~{}~{}~{}~{}x(k)^{T}(C^{T}QC+\Phi^{T}\Omega\Phi)x(k),$

where $U_{k}:=\{u_{0|k},\ldots,u_{N-1|k}\}$ , $x_{0|k}=x(k)$ , $u(k)=u_{0|k}$ , $x_{i+1|k}=Ax_{i|k}+Bu_{i|k}$ and $y_{i|k}=Cx_{i|k}$ . Above $Q$ and $R$ are positive definite symmetric weighting matrices of appropriate dimensions, $N$ is the prediction horizon (we assume that a control horizon equal to the prediction horizon is used for brevity) and $P$ is a terminal weight matrix, which is typically taken equal to the solution of a corresponding discrete–time algebraic Riccati equation (DARE). Moreover, the matrices utilized in (2) are defined as follows:

	$\displaystyle\Phi$	$\displaystyle=\begin{bmatrix}CA\\ CA^{2}\\ \vdots\\ CA^{N}\end{bmatrix},~{}~{}\Gamma=\begin{bmatrix}CB&0&\dddot{}&0\\ CAB&CB&\dddot{}&0\\ \vdots&\vdots&\ddots&\vdots\\ CA^{N-1}B&CA^{N-2}B&\dddot{}&CB\end{bmatrix},$		(3)
	$\displaystyle\Omega$	$\displaystyle=\text{diag}(Q,\dddot{},Q,P),~{}~{}~{}~{}\Psi=\text{diag}(R,R,\dddot{},R).$		(3)

The unconstrained MPC control law (i.e., if output and input constraints are neglected) can be computed by taking the gradient of $J$ with respect to $U_{k}$ and set it equal to $0$ , i.e.:

$\displaystyle\nabla_{U_{k}}J(x(k),U_{k})$	$\displaystyle=2(\Psi+\Gamma^{T}\Omega\Gamma)U_{k}+2\Gamma^{T}\Omega\Phi x(k)=0,$	(4)
$\displaystyle U_{k}^{*}(x(k))$	$\displaystyle=-G^{-1}F\Phi x(k),\quad\text{where}$
$\displaystyle G$	$\displaystyle=2(\Psi+\Gamma^{T}\Omega\Gamma),\quad F=2\Gamma^{T}\Omega.$

When output and input constraints are added to the minimization of $J$ , the constrained MPC control law is computed on–line, by solving a quadratic program, i.e.,

	$\displaystyle\min_{U_{k}}$	$\displaystyle\frac{1}{2}U_{k}^{T}GU_{k}+U_{k}^{T}F\Phi x(k)$		(5)
	subject to:	$\displaystyle\mathcal{L}U_{k}\leq c-\mathcal{D}y(k)-\mathcal{M}\Phi x(k)$		(5)

where the derivation of the matrices $\mathcal{L}$ , $\mathcal{D}$ , $\mathcal{M}$ and the vector $c$ is illustrated next. It is worth mentioning that the terms in (2) that do not depend on $U_{k}$ are omitted in (5), as they do not influence the corresponding optimum.

Consider linear constraints in the outputs and inputs, i.e.

		$\displaystyle M_{i}y_{i\|k}+E_{i}u_{i\|k}\leq b_{i},~{}~{}\forall i=0,1,\dddot{},N-1,\quad\text{and}$		(6)
		$\displaystyle M_{N}y_{N\|k}\leq b_{N},$		(6)

for suitable matrices $M_{i},E_{i}$ and vectors $b_{i}$ . These constraints can be aggregated as follows

$\displaystyle\begin{bmatrix}M_{0}\\ 0\\ \vdots\\ 0\end{bmatrix}y_{0\|k}+\begin{bmatrix}0&\dddot{}&0\\ M_{1}&\dddot{}&0\\ \vdots&\ddots&\vdots\\ 0&\dddot{}&M_{N}\end{bmatrix}\begin{bmatrix}y_{1\|k}\\ \vdots\\ y_{N\|k}\end{bmatrix}$	$\displaystyle+$	(7)
$\displaystyle\begin{bmatrix}E_{0}&\dddot{}&0\\ \vdots&\ddots&\vdots\\ 0&\dddot{}&E_{N-1}\\ 0&\dddot{}&0\end{bmatrix}\begin{bmatrix}u_{0\|k}\\ \vdots\\ u_{N-1\|k}\end{bmatrix}$	$\displaystyle\leq\begin{bmatrix}b_{0}\\ \vdots\\ b_{N}\end{bmatrix},\quad\text{or}$
$\displaystyle\mathcal{D}y(k)+\mathcal{M}Y_{k}+\mathcal{E}U_{k}$	$\displaystyle\leq c.$

Substituting the future predicted outputs $Y_{k}=\Phi x(k)+\Gamma U_{k}$ into equation (7) results in

\displaystyle\mathcal{L}U_{k}

\displaystyle\leq c-\mathcal{M}\Phi x(k)-\mathcal{D}y(k),~{}~{}\text{with}~{}~{}\mathcal{L}=(\mathcal{M}\Gamma+\mathcal{E}).

(8)

The MPC control law is extracted from $U^{\ast}(x(k))$ as follows:

\displaystyle u^{*}(x(k))

\displaystyle=[I_{m}\quad\mathbf{0}_{(N-m)\times m}]U_{k}^{*}(x(k)).

(9)

Next, we will present a data–driven approach to predictive control design, which is based on estimating the prediction matrix $\Gamma$ and the predicted “free response” $\Phi x(k)$ directly from measured input–output (or state) data.

III Data–driven predictive control with estimated prediction matrices

Instead of using an identified state–space model to build the prediction matrices $\Phi$ and $\Gamma$ (see (4) and (8)), we aim to directly compute these matrices from I/O data. The data can either be collected whilst controlling the system on–line or obtained off–line from an experiment on the system. Even though collecting data on–line is preferable, as no prior experiment is required, there is no guarantee of persistence of excitation for the inputs generated by a controller. This can yield not only ill–conditioned prediction matrices, but it can strongly degrade initial performance due to the short measurement horizon. Therefore, in this work we assume the data is measured off–line.

III-A Estimating the prediction matrices

We start by showing how to obtain $\Gamma$ from measured data, using an ARMarkov model as reported in [10, 9, 23]. Suppose the input and output vectors are structured as (for any $k\geq 0,N\geq 1$ ):

\mathbf{y}(k)=\begin{bmatrix}y(k)\\ y(k+1)\\ \vdots\\ y(k+N-1)\end{bmatrix},\quad\mathbf{u}(k)=\begin{bmatrix}u(k)\\ u(k+1)\\ \vdots\\ u(k+N-1)\end{bmatrix}.

These vectors can be stacked in the following Hankel matrices:

$\displaystyle\mathbf{Y}_{p}$	$\displaystyle=\begin{bmatrix}\mathbf{y}(1)&\dddot{}&\mathbf{y}(L+1)\end{bmatrix},$	(10)
$\displaystyle\mathbf{U}_{p}$	$\displaystyle=\begin{bmatrix}\mathbf{u}(0)&\dddot{}&\mathbf{u}(L)\end{bmatrix},$
$\displaystyle\mathbf{Y}_{f}$	$\displaystyle=\begin{bmatrix}\mathbf{y}(N+1)&\dddot{}&\mathbf{y}(N+L+1)\end{bmatrix},$
$\displaystyle\mathbf{U}_{f}$	$\displaystyle=\begin{bmatrix}\mathbf{u}(N)&\dddot{}&\mathbf{u}(N+L)\end{bmatrix},$

representing the so–called “past” and “future” I/O data. Above $L$ represents some chosen measurement horizon, which we will define later. According to [10, 24], the relation between these matrices can be described as follows:

$\displaystyle\mathbf{Y}_{p}$	$\displaystyle=\Phi X_{p}+\Gamma\mathbf{U}_{p},$	(11)
$\displaystyle\mathbf{Y}_{f}$	$\displaystyle=\Phi X_{f}+\Gamma\mathbf{U}_{f},$
$\displaystyle X_{f}$	$\displaystyle=A^{N}X_{p}+D\mathbf{U}_{p},$

where $X_{p}=\begin{bmatrix}x(0)&x(1)&\dddot{}&x(L)\end{bmatrix},~{}X_{f}=\begin{bmatrix}x(N)&x(N+1)&\dddot{}&x(N+L)\end{bmatrix}$ , and $D$ is the extended controllability matrix $D=\begin{bmatrix}A^{N-1}B&\dddot{}&AB&B\end{bmatrix}$ . According to [23], as long as $(N+1)q\geq n$ , it is guaranteed for an observable system that there exists a matrix $M$ such that $A^{N}+M\Phi=0.$ This allows us to rewrite $\mathbf{Y}_{f}$ without state information:

	$\displaystyle\mathbf{Y}_{f}$	$\displaystyle=\Phi A^{N}X_{p}+\Phi D\mathbf{U}_{p}+\Gamma\mathbf{U}_{f}$
		$\displaystyle=-\Phi M\Phi X_{p}+\Phi D\mathbf{U}_{p}+\Gamma\mathbf{U}_{f}$
		$\displaystyle=-\Phi M(\mathbf{Y}_{p}-\Gamma\mathbf{U}_{p})+\Phi D\mathbf{U}_{p}+\Gamma\mathbf{U}_{f}$
		$\displaystyle=\begin{bmatrix}\Phi(D+M\Gamma)&-\Phi M&\Gamma\end{bmatrix}\underbrace{\begin{bmatrix}\mathbf{U}_{p}^{T}&\mathbf{Y}_{p}^{T}&\mathbf{U}_{f}^{T}\end{bmatrix}^{T}}_{W}.$

For given $W$ and $\mathbf{Y}_{f}$ we can obtain the least–squares estimated prediction matrix $\Gamma$ :

\begin{bmatrix}P_{1}&P_{2}&\Gamma\end{bmatrix}=\mathbf{Y}_{f}W^{\dagger},

(12)

where $[.]^{\dagger}$ represents the Moore-Penrose pseudo-inverse, $P_{1}=\Phi(D+M\Gamma)$ and $P_{2}=-\Phi M$ . Hereby it should be noted that $L\geq(N)(2m+q)$ to ensure (12) has a solution.

III-A1 Estimating $\Phi$ : measurable state

If the system state is measurable, an estimate of $\Phi$ can be accurately obtained from data. The following ways to estimate it are derived from [10]. The first one uses the estimated $\hat{\Gamma}$ and (11), as follows:

\mathbf{Y}_{p}=\Phi X_{p}+\Gamma\mathbf{U}_{p},~{}~{}\text{which yields}~{}~{}\Phi=(\mathbf{Y}_{p}-\hat{\Gamma}\mathbf{U}_{p})X_{p}^{\dagger}.

(13)

Another approach is using the orthogonal complement of the row space of the matrix $\mathbf{U}_{p}$ :

$\displaystyle\mathbf{U}_{po}$	$\displaystyle:=I-\mathbf{U}_{p}^{T}(\mathbf{U}_{p}\mathbf{U}_{p}^{T})^{-1}\mathbf{U}_{p}$	(14)
$\displaystyle\mathbf{Y}_{p}\mathbf{U}_{po}$	$\displaystyle=\Phi X_{p}\mathbf{U}_{po}+\Gamma\mathbf{U}_{p}\mathbf{U}_{po}=\Phi X_{p}\mathbf{U}_{po},$
$\displaystyle\Phi$	$\displaystyle=(\mathbf{Y}_{p}\mathbf{U}_{po})(X_{p}\mathbf{U}_{po})^{\dagger}.$

Remark 1

(The LQR equivalence) When the state is measurable, one can also retrieve the stability and (sub–)optimal performance properties of the LQR controller even for a small horizon $N$ . Consider the following LQR and MPC cost functions, respectively:

	$\displaystyle J(y(k),U_{k})$	$\displaystyle=\sum_{i=0}^{\infty}(y_{i\|k}^{T}Qy_{i\|k}+u_{i\|k}^{T}Ru_{i\|k})$
	$\displaystyle J(y(k),U_{k})$	$\displaystyle=x_{N\|k}^{T}Q_{p}x_{N\|k}+\sum_{i=0}^{N-1}(y_{i\|k}^{T}Qy_{i\|k}+u_{i\|k}^{T}Ru_{i\|k})$

Notice that here $Q_{p}$ is a terminal penalty on the state vector, not on the output. If $Q_{p}$ equals the solution of the associated DARE for the LQR problem, then we obtain the infinite horizon LQR control law as the unconstrained MPC control law¹¹1For sufficient large $N$ there exists a closed-form solution for $Q_{p}$ (see, e.g., [10]), but then the MPC cost function also converges to the LQR case..

In order to introduce the terminal state penalty in the developed data–driven predictive controller, one can create an unobservable output $y_{o}=Vx$ , where $V$ is the Cholesky decomposition of $Q_{p}$ , such that $y_{o}^{T}y_{o}=x^{T}Q_{p}x$ . Thus, we have an augmented output vector $y_{a}=[y^{T}~{}~{}y_{o}^{T}]^{T}$ which needs only to be re-substituted in (10). This also yields the following performance matrix changes: $\tilde{Q}=\text{diag}(Q,\mathbf{0}_{n}),$ $\tilde{P}=\text{diag}(\mathbf{0}_{q},\mathbf{1}_{n}).$ For sub-optimality it suffices to choose a $\tilde{Q}_{p}>Q_{p}$ .

III-A2 Estimating $\Phi x(k)$ : only I/O data available

If the state is not measurable, the method of [9] can be used to estimate $\Phi x(k)$ on–line, using previous inputs and measured outputs, as follows:

$\displaystyle\Phi x(k)$	$\displaystyle=\Phi A^{N}x(k-N)+\Phi D\mathbf{u}(k-N)$	(15)
	$\displaystyle=-\Phi M\Phi x(k-N)+\Phi D\mathbf{u}(k-N)$
	$\displaystyle=-\Phi M\mathbf{y}(k-N+1)+\Phi(D+M\Gamma)\mathbf{u}(k-N).$

Since both $\Phi(D+M\Gamma)$ and $-\Phi M$ can be directly extracted from the matrix obtained in equation (12), then one can compute $\Phi x(k)$ in real–time as

\Phi x(k)=\begin{bmatrix}P_{1}&P_{2}\end{bmatrix}\begin{bmatrix}\mathbf{u}(k-N)\\ \mathbf{y}(k-N+1)\end{bmatrix}.

(16)

On–line estimation of $\Phi x(k)$ depends on the unknown $M$ and the estimate accuracy is not guaranteed to improve for any sequence length shorter than $N$ .

Remark 2

(Initial inputs and feasibility) If $\Phi x(k)$ is computed using previous I/O data, it should be noted that due to this poor initial estimation, the predictive control quadratic program may not be feasible at start. This can be circumvented via soft–constraints.

III-B Integral action

Embedding integral action within a model–based predictive controller is known to be beneficial both with respect to tracking time–varying references and removing off–sets. To this end, the rate–based MPC formulation was developed (see, e.g., [22]) by defining an augmented state $x_{I}(k):=\begin{bmatrix}\Delta x(k)^{T}&y(k)^{T}\end{bmatrix}^{T}$ , where $\Delta x(k):=x(k)-x(k-1)$ and an incremental input $\Delta u(k):=u(k)-u(k-1)$ . Then we can construct an augmented system model as follows:

	$\displaystyle x_{I}(k+1)$	$\displaystyle=\underbrace{\begin{bmatrix}A&0\\ CA&I\end{bmatrix}}_{A_{I}}x_{I}(k)+\underbrace{\begin{bmatrix}B\\ CB\end{bmatrix}}_{B_{I}}\Delta u(k)$		(17)
	$\displaystyle y(k)$	$\displaystyle=\underbrace{\begin{bmatrix}0&I\end{bmatrix}}_{C_{I}}x_{I}(k).$		(17)

The cost function changes to:

	$\displaystyle J(x(k),\Delta U_{k})$	$\displaystyle:=(y_{N\|k}-r_{N\|k})^{T}P(y_{N\|k}-r_{N\|k})$
	$\displaystyle+\sum_{i=0}^{N-1}$	$\displaystyle\left((y_{i\|k}-r_{i\|k})^{T}Q(y_{i\|k}-r_{i\|k})+\Delta u_{i\|k}^{T}R\Delta u_{i\|k}\right)$
		$\displaystyle=\Delta U_{k}^{T}(\Psi+\Gamma_{I}^{T}\Omega\Gamma_{I})\Delta U_{k}$
		$\displaystyle~{}~{}+2\Delta U_{k}^{T}\Gamma_{I}^{T}\Omega(\Phi_{I}x_{I}(k)-\mathcal{R}_{k})$
		$\displaystyle~{}~{}+(y(k)-r(k))^{T}Q(y(k)-r(k))$
		$\displaystyle~{}~{}+(\Phi_{I}x_{I}(k)-\mathcal{R}_{k})^{T}\Omega(\Phi_{I}x_{I}(k)-\mathcal{R}_{k}).$

where $\Delta U_{k}=\{\Delta u_{0|k},\ldots,\Delta u_{N-1|k}\}$ , $r_{0|k}=r(k)$ and $\mathcal{R}_{k}=\begin{bmatrix}r_{1|k}^{T}&r_{2|k}^{T}&\dddot{}&r_{N|k}^{T}\end{bmatrix}^{T}$ is a known future reference. The actual control input is obtained as $u(k)=u(k-1)+\Delta u^{*}(k)$ . By substituting the integral model (17) into the definitions in equation (3), we obtain

	$\displaystyle\Phi_{I}$	$\displaystyle=\begin{bmatrix}CA&I\\ CA^{2}+CA&I\\ \vdots\\ \sum_{i=1}^{N}CA^{i}&I\end{bmatrix},$		(18)
	$\displaystyle\Gamma_{I}$	$\displaystyle=\begin{bmatrix}CB&0&\dddot{}&0\\ CAB+CB&CB&\dddot{}&0\\ \vdots&\vdots&\ddots&\vdots\\ \sum_{i=0}^{N-1}CA^{i}B&\sum_{i=0}^{N-2}CA^{i}B&\dddot{}&CB\end{bmatrix}.$		(18)

We now formulate a method that allows us to derive these integral action prediction matrices directly from data by manipulating the input used in the experiment.

Theorem 1

Let $[U_{I}]_{i|k}=\begin{bmatrix}u_{0|k}\\ u_{0|k}+u_{1|k}\\ \vdots\\ \sum_{j=0}^{i}u_{j|k}\end{bmatrix}+\begin{bmatrix}1\\ 1\\ \vdots\\ 1\end{bmatrix}\otimes\underbrace{\left(\sum_{j=0}^{k-1}u(j)\right)}_{u_{I}(k-1)},$ and $[\Phi]_{i}$ , $[\Gamma]_{i}$ (respectively $[\Phi_{I}]_{i}$ , $[\Gamma_{I}]_{i}$ ) represent the corresponding row-block w.r.t. output $y_{i|k}$ . Then it holds that

	$\displaystyle y_{i\|k}$	$\displaystyle=[\Phi_{I}]_{i}[x_{I}]_{0\|k}+[\Gamma_{I}]_{i}U_{i\|k}$		(19)
		$\displaystyle=[\Phi]_{i}x_{0\|k}+[\Gamma]_{i}[U_{I}]_{i\|k},~{}~{}\forall i=1,\dddot{},N.$		(19)

Proof:

Consider first the trivial case, $i=1$ , and let $y_{0|k}=CAx(k-1)+CBu_{I}(k-1)$ :

$\displaystyle y_{1\|k}$	$\displaystyle=\begin{bmatrix}CA&I\end{bmatrix}[x_{I}]_{0\|k}+CBu_{0\|k}$	(20)
	$\displaystyle=CAx_{0\|k}-CAx(k-1)+y_{0\|k}+CBu_{0\|k}$
	$\displaystyle=CAx_{0\|k}-CAx(k-1)+CAx(k-1)$
	$\displaystyle~{}~{}~{}~{}+CBu_{I}(k-1)+CBu_{0\|k}$
	$\displaystyle=CAx_{0\|k}+CB(u_{I}(k-1)+u_{0\|k})$
	$\displaystyle=[\Phi]_{1}x_{0\|k}+[\Gamma]_{1}[U_{I}]_{1\|k}.$

Next, consider the case $i=2$ , and the previous relation for $y_{0|k}$ :

$\displaystyle y_{2\|k}$	$\displaystyle=\begin{bmatrix}CA^{2}+CA\ &I\end{bmatrix}[x_{I}]_{0\|k}$	(21)
	$\displaystyle~{}~{}~{}~{}+\begin{bmatrix}CAB+CB&CB\end{bmatrix}\begin{bmatrix}u_{0\|k}\\ u_{1\|k}\end{bmatrix}$
	$\displaystyle=(CA^{2}+CA)(Ax(k-1)+Bu_{I}(k-1))$
	$\displaystyle~{}~{}~{}~{}-(CA^{2}+CA)x(k-1)+CAx(k-1)$
	$\displaystyle~{}~{}~{}~{}+CBu_{I}(k-1)+(CAB+CB)u_{0\|k}+CBu_{1\|k}$
	$\displaystyle=CA^{2}(Ax(k-1)+Bu_{I}(k-1))$
	$\displaystyle~{}~{}~{}~{}+(CAB+CB)(u_{I}(k-1)+u_{0\|k})+CBu_{1\|k}$
	$\displaystyle=CA^{2}{x}_{0\|k}+\begin{bmatrix}CAB&CB\end{bmatrix}\begin{bmatrix}u_{0\|k}+u_{I}(k-1)\\ u_{1\|k}+u_{0\|k}+u_{I}(k-1)\end{bmatrix}$
	$\displaystyle=[\Phi]_{2}x_{0\|k}+[\Gamma]_{2}[U_{I}]_{2\|k}.$

Thus, we obtained the prediction of output as a function of the original prediction matrices, but for the integrated input. The proof for any other $i\geq 3$ follows using the same derivations and is omitted due to space limitation. ∎

The result of Theorem 1 suggests that it is possible to estimate the prediction matrices corresponding to the rate–based predictive control algorithm by manipulating the data fed to the system in the experiment. More specifically, the following procedure can be applied. First, create an input set $U$ and form a second input set $U_{I}$ from $U$ such that for any $k\geq 0$ , $U_{I}(k)=\sum_{i=0}^{k}U(i)$ . Next, perform the experiment on the system using $U_{I}$ as input and collect the output, which we will denote $Y_{I}$ . Finally, compute $\Gamma_{I}$ as in (12) using $Y_{I}$ for $\mathbf{Y}_{p}$ and $\mathbf{Y}_{f}$ , but using $U$ for the input matrices $\mathbf{U}_{p}$ and $\mathbf{U}_{f}$ , instead of the applied $U_{I}$ . Using $U_{I}$ in the experiment, however, could lead to a high magnitude input. An alternative procedure can be used to circumvent this problem, as follows: apply $U$ in the experiment, collect $Y$ and use $U_{\Delta}(k)=U(k)-U(k-1)$ to compose $\mathbf{U}_{p}$ and $\mathbf{U}_{f}$ . The proof that the developed estimation procedure can still be applied follows the same derivations as above, but we do not include it here due to space limitations.

Remark 3

The implementation of the integral action DPC still requires estimation of $\Phi_{I}x_{I}(k)$ . If the states are measurable (see Section III-A1) the unconstrained integral DPC law becomes:

\displaystyle\Delta U_{k}^{*}(x(k))

\displaystyle=-G^{-1}F\left(\Phi_{I}\begin{bmatrix}x(k)-x(k-1)\\ y(k)\end{bmatrix}-\mathcal{R}_{k}\right).

(22)

To obtain an estimate of $\Phi_{I}$ , equations (13) or (14) can be used in combination with the state vector corresponding to the integral action augmented state, i.e., $x_{I}(k)=[\Delta x(k)^{T}y(k)^{T}]^{T}$ . If the state is not measurable, the unconstrained integral DPC law changes to

	$\displaystyle\Delta U_{k}^{*}(x(k))$	$\displaystyle=-G^{-1}F\left(\begin{bmatrix}P_{1}&P_{2}\end{bmatrix}\begin{bmatrix}\Delta\mathbf{u}(k-N)\\ \mathbf{y}(k-N+1)\end{bmatrix}-\mathcal{R}_{k}\right)$
		$\displaystyle=-G^{-1}F\left(\Phi_{I}x_{I}(k)-\mathcal{R}_{k}\right),$

which is similar to the method used in Section III-A2. If constraints are added to the optimization problem, the same approaches to estimate $\Phi_{I}$ or $\Phi_{I}x_{I}(k)$ , respectively, apply, as the constraints matrices (8) (corresponding to the integral action formulation) only depend on the estimated prediction matrices.

IV Illustrative examples

IV-A Position control for a linear motor

In this section we provide an illustrative example of the proposed methodology for the position control of a linear motor, see, e.g., [25]. The simplified motion dynamics of the mechanical part of the actuator can be represented by the continuous–time model

\dot{x}(t)=\begin{bmatrix}0&1\\ 0&-50\end{bmatrix}x(t)+\begin{bmatrix}0\\ 0.5\end{bmatrix}u(t)+\begin{bmatrix}0\\ 0.5\end{bmatrix}d(t),

y(t)=\begin{bmatrix}1&0\end{bmatrix}x(t),

where the first state (the output) is the position and the second state is the velocity. The system is discretized with a sampling time of $T_{s}=0.02$ s. In order collect data from the system we set an experiment where the input is a PRBS signal with an amplitude of $\pm 20$ N at $50$ Hz and length 6022 samples. The output is corrupted by a white noise sequence such that the Signal-to-Noise Ratio (SNR) is $26$ dB.

The control performance parameters are given by $N=10$ , $Q=6\text{e}5$ , and $R=0.005$ . Also the input force is constrained between $-500~{}\text{N}\leq u(k)\leq 500~{}\text{N}$ , and displacement is limited between $-0.165~{}\text{m}\leq y(k)\leq 0.165~{}\text{m}$ . We then perform the estimation of the prediction matrices in (12) and (14) with $L=3000$ , and similarly for the integral DPC case.

Fig. 1 portrays the closed–loop results for a sequence of steps using 3 predictive controllers: standard model–based PC, the developed DPC with measured state and measured output, respectively. Also, between $t=3$ s and $t=5$ s we apply a constant disturbance force of amplitude $d(k)=-100$ N and we also consider the output corrupted by white noise with SNR = $25$ dB. It can be observed that the DPC with measured state converges faster to the reference compared to the DPC with measured output, and it matches the standard MPC result. Due to the poor initial estimation of the output DPC, the initial response is a bit slower. When the disturbance becomes active all 3 predictive controllers result in an off–set, which is considerably larger for the output DPC.

Refer to caption — Figure 1: Linear motor comparison of 3 predictive controllers: standard MPC; state DPC (with measured state); output DPC (with measured output).

Fig. 2 shows the closed–loop simulation results for the corresponding 3 integral predictive controllers. All controllers successfully remove the off-set caused by the disturbance, while the output integral DPC exhibits longer transients due to the on–line estimation part. This really demonstrates the need of incorporating integral action in data–driven predictive controllers and also, the effectiveness of the developed DPC algorithm.

IV-B UPS system with periodic disturbance

To provide another practical example, the data–driven control algorithm with integral action is also applied to an uninterruptible power supply (UPS). This plant has been studied before in [10].

Fig. 3 shows a simplified electrical diagram of the output stage of a single-phase UPS system. The load effect on the system output is modeled by a parallel connection of an uncertain admittance $Y_{0}(t)$ and an unknown periodic disturbance given by the current source $i_{d}(t)$ . Again for simplicity, the uncertain admittance $Y_{0}(t)$ is presumed to be constant at $0.1519\Omega$ . The input is a PWM control voltage and the output is the voltage over the capacitor. The states are the inductor current and the capacitor voltage, $x(t)=\begin{bmatrix}i(t)&v(t)\end{bmatrix}^{T}$ . The system can be represented as the following state-space model:

\dot{x}(t)=\begin{bmatrix}-\frac{R_{Lf}}{L_{f}}&-\frac{1}{L_{f}}\\ \frac{1}{C_{f}}&-\frac{Y_{0}}{C_{f}}\end{bmatrix}x(t)+\begin{bmatrix}\frac{K_{PWM}}{L_{f}}\\ 0\end{bmatrix}u(t)+\begin{bmatrix}0\\ \frac{-1}{C_{f}}\end{bmatrix}i_{d}(t)

	$\displaystyle\dot{x}(t)$	$\displaystyle=\begin{bmatrix}-15&-1000\\ \frac{1}{3e-4}&-506.46\end{bmatrix}x(t)+\begin{bmatrix}1000\\ 0\end{bmatrix}u(t)+\begin{bmatrix}0\\ \frac{-1}{3e-4}\end{bmatrix}i_{d}(t)$		(23)
	$\displaystyle y(t)$	$\displaystyle=\begin{bmatrix}1&0\end{bmatrix}x(t)$		(23)

The system is discretized with a sampling time of $T_{s}=1/15000~{}s$ using a zero-order hold. The input voltage is constrained between $-260\leq u(t)\leq 260$ and there are no output constraints. The closed–loop is set to follow a sinusoidal reference $r(t)=127\sqrt{2}\text{cos}(120\pi t)$ . Additionally, the following control parameters are used:

Q=200,~{}~{}R=50,~{}~{}N=15,~{}~{}L=3000

The system is measured off–line using a PRBS with an amplitude of $\pm 104V$ with a length $7500$ samples at $15000$ Hz. The experiment is corrupted with output noise such that the SNR is $28$ dB, the disturbance current was not present in the experiment. During the simulation, the system is also subjected to an output noise with a SNR of $36$ dB. The periodic disturbance is defined as multi–sine with 12 components with random frequencies between $2000$ and $12000$ Hz and has an amplitude of $\pm 20A$ . See Fig. 4 for the result of the simulation. The various controllers all achieve reference tracking even for a periodic reference. The total harmonic distortion (THD) caused by the periodic disturbance on the system is between 3.8% and 4%.

V Conclusion

In this work we proposed a new data–driven predictive control approach with estimated prediction matrices and integral action. Our approach lies between the direct data–enabled predictive control and the identified model–based predictive control. In this way, we are able to avoid identification of a complete state–space model and the respective system order, while still enjoying standard properties of least–squares estimation. Also, comparatively to the data–enabled approach, our method deals well with noisy data and requires less in general computations for the on–line implementation. We also presented a formulation that allows direct incorporation of integral action in the DPC algorithm by only manipulating data that is fed to algorithm, without altering the single experiment required for estimation. Simulation results for position control of a linear motor illustrated the effectiveness of our methodology in terms of tracking performance, robustness to noisy measurements and disturbance rejection.

References

[1] F. Lamnabhi-Lagarrigue, A. Annaswamy, S. Engell, A. Isaksson, P. Khargonekar, R. M. Murray, H. Nijmeijer, T. Samad, D. Tilbury, and P. Van den Hof, “Systems & control for the future of humanity, research agenda: Current and future roles, impact and grand challenges,” Annual Reviews in Control, vol. 43, pp. 1–64, 2017.
[2] I. Markovsky and P. Rapisarda, “Data-driven simulation and control,” International Journal of Control, vol. 81, no. 12, pp. 1946–1959, 2008.
[3] A. S. Bazanella, L. Campestrini, and D. Eckhard, Data-driven controller design: the H2 approach. Netherlands: Springer Science & Business Media, 2011.
[4] Z.-S. Hou and Z. Wang, “From model-based control to data-driven control: Survey, classification and perspective,” Information Sciences, vol. 235, pp. 3–35, 2013.
[5] M. C. Campi, A. Lecchini, and S. M. Savaresi, “Virtual reference feedback tuning: a direct method for the design of feedback controllers,” Automatica, vol. 38, no. 8, pp. 1337 – 1346, 2002.
[6] L. Campestrini, D. Eckhard, A. S. Bazanella, and M. Gevers, “Data-driven model reference control design by prediction error identification,” Journal of the Franklin Institute, vol. 354, no. 6, pp. 2828–2647, 2017.
[7] G. R. Gonçalves da Silva, A. S. Bazanella, and L. Campestrini, “One-shot data-driven controller certification,” ISA Transactions, vol. 99, pp. 361–373, 2020.
[8] H. J. Van Waarde, J. Eising, H. L. Trentelman, and M. K. Camlibel, “Data informativity: a new perspective on data-driven analysis and control,” IEEE Transactions on Automatic Control, 2020.
[9] W. Aangenent, D. Kostic, B. de Jager, R. van de Molengraft, and M. Steinbuch, “Data-based optimal control,” in Proceedings of the American Control Conference, vol. 2, 2005, pp. 1460–1465.
[10] G. R. Gonçalves da Silva, A. S. Bazanella, C. Lorenzini, and L. Campestrini, “Data-driven LQR control design,” IEEE Control Systems Letters, vol. 3, no. 1, pp. 180–185, 2019.
[11] C. d. Persis and P. Tesi, “Formulas for data-driven control: Stabilization, optimality, and robustness,” IEEE Transactions on Automatic Control, vol. 65, no. 3, pp. 909–924, 2019.
[12] J. B. Rawlings, D. Q. Mayne, and M. M. Diehl, Model Predictive Control: Theory, Computation, and Design, 2nd ed. Nob Hill Publishing, 2017.
[13] D. W. Clarke, C. Mohtadi, and P. S. Tuffs, “Generalized predictive control–Part I. the basic algorithm,” Automatica, vol. 23, no. 2, pp. 137–148, 1987.
[14] G. Prasath and J. B. Jørgensen, “Model predictive control based on finite impulse response models,” in 2008 American Control Conference. IEEE, 2008, pp. 441–446.
[15] A. Carron, E. Arcari, M. Wermelinger, L. Hewing, M. Hutter, and M. N. Zeilinger, “Data-driven model predictive control for trajectory tracking with a robotic arm,” IEEE Robotics and Automation Letters, vol. 4, no. 4, pp. 3758–3765, 2019.
[16] I. Markovsky, J. C. Willems, S. Van Huffel, and B. De Moor, Exact and approximate modeling of linear systems: A behavioral approach. SIAM, 2006.
[17] H. Yang and S. Li, “A data-driven predictive controller design based on reduced hankel matrix,” in 2015 10th Asian Control Conference (ASCC). IEEE, 2015, pp. 1–7.
[18] J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control: In the shallows of the DeePC,” in 2019 18th European Control Conference (ECC). IEEE, 2019, pp. 307–312.
[19] J. Berberich, J. Köhler, M. A. Muller, and F. Allgower, “Data-driven model predictive control with stability and robustness guarantees,” IEEE Transactions on Automatic Control, 2020.
[20] W. Favoreel, B. De Moor, and M. Gevers, “SPC: Subspace predictive control,” IFAC Proceedings Volumes, vol. 32, no. 2, pp. 4004–4009, 1999.
[21] R. Kadali, B. Huang, and A. Rossiter, “A data driven subspace approach to predictive controller design,” Control engineering practice, vol. 11, no. 3, pp. 261–278, 2003.
[22] L. Wang, Model predictive control system design and implementation using MATLAB®. Springer Science & Business Media, 2009.
[23] R. K. Lim, M. Q. Phan, and R. W. Longman, “State estimation with ARMarkov models,” Department of mechanical and aerospace engineering, Princeton University, Princeton, NJ, Tech. Rep. 3046, October 1998.
[24] P. Van Overschee and B. De Moor, Subspace identification for linear systems. Theory, implementation, applications. Springer, Boston, MA, 01 1996, vol. xiv, pp. xiv + 254.
[25] T. T. Nguyen, M. Lazar, and H. Butler, “Nonlinear model predictive control for ironless linear motors,” in IEEE Conference on Control Technology and Applications. IEEE, 2018, pp. 927–932.

	$\displaystyle y_{i\|k}$	$\displaystyle=[\Phi_{I}]_{i}[x_{I}]_{0\|k}+[\Gamma_{I}]_{i}U_{i\|k}$		(19)
		$\displaystyle=[\Phi]_{i}x_{0\|k}+[\Gamma]_{i}[U_{I}]_{i\|k},~{}~{}\forall i=1,\dddot{},N.$		(19)

$\displaystyle y_{1\|k}$	$\displaystyle=\begin{bmatrix}CA&I\end{bmatrix}[x_{I}]_{0\|k}+CBu_{0\|k}$	(20)
	$\displaystyle=CAx_{0\|k}-CAx(k-1)+y_{0\|k}+CBu_{0\|k}$
	$\displaystyle=CAx_{0\|k}-CAx(k-1)+CAx(k-1)$
	$\displaystyle~{}~{}~{}~{}+CBu_{I}(k-1)+CBu_{0\|k}$
	$\displaystyle=CAx_{0\|k}+CB(u_{I}(k-1)+u_{0\|k})$
	$\displaystyle=[\Phi]_{1}x_{0\|k}+[\Gamma]_{1}[U_{I}]_{1\|k}.$

$\displaystyle y_{2\|k}$	$\displaystyle=\begin{bmatrix}CA^{2}+CA\ &I\end{bmatrix}[x_{I}]_{0\|k}$	(21)
	$\displaystyle~{}~{}~{}~{}+\begin{bmatrix}CAB+CB&CB\end{bmatrix}\begin{bmatrix}u_{0\|k}\\ u_{1\|k}\end{bmatrix}$
	$\displaystyle=(CA^{2}+CA)(Ax(k-1)+Bu_{I}(k-1))$
	$\displaystyle~{}~{}~{}~{}-(CA^{2}+CA)x(k-1)+CAx(k-1)$
	$\displaystyle~{}~{}~{}~{}+CBu_{I}(k-1)+(CAB+CB)u_{0\|k}+CBu_{1\|k}$
	$\displaystyle=CA^{2}(Ax(k-1)+Bu_{I}(k-1))$
	$\displaystyle~{}~{}~{}~{}+(CAB+CB)(u_{I}(k-1)+u_{0\|k})+CBu_{1\|k}$
	$\displaystyle=CA^{2}{x}_{0\|k}+\begin{bmatrix}CAB&CB\end{bmatrix}\begin{bmatrix}u_{0\|k}+u_{I}(k-1)\\ u_{1\|k}+u_{0\|k}+u_{I}(k-1)\end{bmatrix}$
	$\displaystyle=[\Phi]_{2}x_{0\|k}+[\Gamma]_{2}[U_{I}]_{2\|k}.$