Next Article in Journal
Studying the Harmonic Functions Associated with Quantum Calculus
Previous Article in Journal
A Data-Driven Decision-Making Model for Configuring Surgical Trays Based on the Likelihood of Instrument Usages
Previous Article in Special Issue
Resources Relocation Support Strategy Based on a Modified Genetic Algorithm for Bike-Sharing Systems
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Generalized Data–Driven Predictive Control: Merging Subspace and Hankel Predictors

Control Systems Group, Eindhoven University of Technology, 5612 AZ Eindhoven, The Netherlands
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(9), 2216; https://doi.org/10.3390/math11092216
Submission received: 31 March 2023 / Revised: 28 April 2023 / Accepted: 30 April 2023 / Published: 8 May 2023
(This article belongs to the Special Issue Information Theory Applied in Scientific Computing)

Abstract

:
Data–driven predictive control (DPC) is becoming an attractive alternative to model predictive control as it requires less system knowledge for implementation and reliable data is increasingly available in smart engineering systems. Two main approaches exist within DPC: the subspace approach, which estimates prediction matrices (unbiased for large data) and the behavioral, data-enabled approach, which uses Hankel data matrices for prediction (allows for optimizing the bias/variance trade–off). In this paper we develop a novel, generalized DPC (GDPC) algorithm by merging subspace and Hankel predictors. The predicted input sequence is defined as the sum of a known, baseline input sequence, and an optimized input sequence. The corresponding baseline output sequence is computed using an unbiased, subspace predictor, while the optimized predicted output sequence is computed using a Hankel matrix predictor. By combining these two types of predictors, GDPC can achieve high performance for noisy data even when using a small Hankel matrix, which is computationally more efficient. Simulation results for a benchmark example from the literature show that GDPC with a reduced size Hankel matrix can match the performance of data–enabled predictive control with a larger Hankel matrix in the presence of noisy data.

1. Introduction

Reliable data is becoming increasingly available in modern, smart engineering systems, including mechatronics, robotics, power electronics, automotive systems, and smart infrastructures, see, e.g., [1,2] and the references therein. For these application domains, model predictive control (MPC) [3,4] has become the preferred advanced control method for several reasons, including constraints handling, anticipating control actions, and optimal performance. Since obtaining and maintaining accurate models requires effort and reliable data becomes readily available in engineering systems, it is of interest to develop data–driven predictive control (DPC) algorithms that can be implemented in practice. An indirect data–driven approach to predictive control design was already developed in [5] more than 20 years ago, i.e., subspace predictive control (SPC). The SPC approach skips the identification of the prediction model and identifies the complete prediction matrices from input–output data using least squares. This provides an unbiased predictor for sufficiently large data.
More recently, a direct data–driven approach to predictive control design was developed in [6] based on behavioral systems theory and Willems’ fundamental lemma [7], i.e., data–enabled predictive control (DeePC). The idea to use (reduced order) Hankel matrices as predictors has been put forward earlier in [8], but the first well–posed constrained data–enabled predictive control algorithm was formulated in [6], to the best of the authors’ knowledge. The DeePC approach skips the identification of prediction models or matrices all together and utilizes Hankel matrices built from input–output data to parameterize predicted future inputs and outputs. In the deterministic, noise free case, equivalence of MPC and DeePC was established in [6,9], while equivalence of SPC and DeePC was shown in [10]. Stability guarantees for DeePC were first obtained in [9] by means of terminal equality constraints and input–output–to–state stability Lyapunov functions. Alternatively, stability guarantees for DeePC were provided in [11] using terminal inequality constraints and dissipation inequalities involving storage and supply functions. An important contribution to DeePC is the consistent regularization cost introduced in [12], which enables reliable performance in the presence of noisy data. Indeed, since the DeePC algorithm jointly solves estimation and controller synthesis problems, the regularization derived in [12] allows one to optimize the bias/variance trade–off if data is corrupted by noise. A systematic method for tuning the regularization cost weighting parameter was recently presented in [13].
Computationally, SPC has the same number of optimization variables as MPC, which is equal to the number of control inputs times the prediction horizon. In DeePC, the number of optimization variables is also driven by the data length, which must be in general much larger than the prediction horizon. Especially in the case of noisy data, a large data size is required to attain reliable predictions, see, e.g., [9,12,13]. As this hampers real–time implementation, it is of interest to improve computational efficiency of DeePC. In [14], a computationally efficient formulation of DeePC was provided via LQ factorization of the Hankel data matrix, which yields the same online computational complexity as SPC/MPC. In this approach, DeePC yields an unbiased predictor, similar to SPC. In [15], a singular value decomposition is performed on the original Hankel data matrix and a DeePC algorithm is designed based on the resulting reduced Hankel matrix. Therein, it was shown that this approach can significantly reduce the computational complexity of DeePC, while improving the accuracy of predictions for noisy data. In [16], an efficient numerical method that exploits the structure of Hankel matrices was developed for solving quadratic programs (QPs) specific to DeePC. Regarding real–life applications of DeePC, the minimal data size required for persistence of excitation is typically used, see, e.g., [17,18], or an unconstrained solution of DeePC is used instead of solving a QP, see, e.g., [19]. These approaches however limit the achievable performance in the presence of noisy data and hard constraints, respectively.
In this paper we develop a novel, generalized DPC (GDPC) algorithm by merging subspace and Hankel predictors with the goal of reducing online computational complexity without sacrificing control performance in the presence of noisy data. The predicted input sequence is defined as the sum of a known, baseline input sequence and an optimized input sequence. The corresponding baseline output sequence is computed using an unbiased, subspace predictor based on a large data set, while the optimized predicted output sequence is computed using a Hankel matrix predictor based on a reduced data set. Via the extension of Willems’ fundamental lemma to multiple data sets [20], the sum of the two trajectories spanned by two (possibly different) data sets will remain a valid system trajectory, as long as the combined data matrices are collectively persistently exciting and the system is linear. By combining these two types of predictors, GDPC can achieve high performance in the presence of noisy data even when using Hankel matrices of smaller size, which is computationally efficient and preserves the bias–variance tradeoff benefits of DeePC. The performance and computational complexity of GDPC with a minimal (according to DeePC design criteria) size Hankel matrix is evaluated for a benchmark example from the MPC literature and compared to DeePC with a Hankel matrix of varying size.
The remainder of this paper is structured as follows. The necessary notation and the DeePC approach to data–driven predictive control are introduced in Section 2. The GDPC algorithm is presented in Section 3, along with design guidelines, stability analysis and other relevant remarks. Simulation results and a comparison with DeePC are provided in Section 4 for a benchmark example from the literature. Conclusions are summarized in Section 5. Frequently used abbreviations are summarized in Table 1.

2. Preliminaries

Consider a discrete–time linear dynamical system subject to zero–mean Gaussian noise w ( k ) N ( 0 , σ w 2 I ) :
x ( k + 1 ) = A x ( k ) + B u ( k ) , k N , y ( k ) = C x ( k ) + w ( k ) ,
where x R n is the state, u R n u is the control input, y R n y is the measured output and ( A , B , C ) are real matrices of suitable dimensions. We assume that ( A , B ) is controllable and ( A , C ) is observable. By applying a persistently exciting input sequence { u ( k ) } k N [ 0 , T ] of length T to system (1) we obtain a corresponding output sequence { y ( k ) } k N [ 0 , T ] .
If one considers an input–output model corresponding to (1), it is necessary to introduce the parameter T ini that limits the window of past input–output data necessary to compute the current output, i.e.,
y ( k ) = i = 1 T ini a i y ( k i ) + i = 1 T ini b i u ( k i ) ,
for some real–valued coefficients. For simplicity of exposition we assume the same T ini for inputs and outputs.
Next, we introduce some instrumental notation. For any finite number q N 1 of vectors { ξ 1 , , ξ q } R n 1 × × R n q we will make use of the operator col ( ξ 1 , , ξ q ) : = [ ξ 1 , , ξ q ] . For any k 0 (starting time instant in the data vector) and j 1 (length of the data vector), define
u ¯ ( k , j ) : = col ( u ( k ) , , u ( k + j 1 ) ) , y ¯ ( k , j ) : = col ( y ( k ) , , y ( k + j 1 ) ) .
Let N T ini denote the prediction horizon. Then we can define the Hankel data matrices:
U p : = u ¯ ( 0 , T ini ) u ¯ ( T 1 , T ini ) , Y p : = y ¯ ( 1 , T ini ) y ¯ ( T , T ini ) , U f : = u ¯ ( T ini , N ) u ¯ ( T ini + T 1 , N ) , Y f : = y ¯ ( T ini + 1 , N ) y ¯ ( T ini + T , N ) .
According to the DeePC design [6], one must choose T ini n and T ( n u + 1 ) ( T ini + N + n ) 1 , which implicitly requires an assumption on the system order (number of states). Given the measured output y ( k ) at time k N and T ini N 1 we define the sequences of known input–output data at time k T ini , which are trajectories of system (1):
u ini ( k ) : = col ( u ( k T ini ) , , u ( k 1 ) ) , y ini ( k ) : = col ( y ( k T ini + 1 ) , , y ( k ) ) .
Next, we define the sequences of predicted inputs and outputs at time k T ini , which should also be trajectories of system (1):
u ( k ) : = col ( u ( 0 | k ) , , u ( N 1 | k ) ) , y ( k ) : = col ( y ( 1 | k ) , , y ( N | k ) ) ,
where the notation y ( i | k ) denotes the predicted value of y ( k + i ) using measured output data available up to y ( k ) .
For a positive definite matrix L let L 1 2 denote its Cholesky factorization. At time k T ini , given u ini ( k ) , y ini ( k ) , the regularized DeePC algorithm [12] computes a sequence of predicted inputs and outputs as follows:
min g ( k ) , u ( k ) , y ( k ) , σ ( k ) l N ( y ( N | k ) ) + i = 0 N 1 l ( y ( i | k ) , u ( i | k ) ) + λ g l g ( g ( k ) ) + λ σ l σ ( σ ( k ) ) subject to constraints :
U p Y p U f Y f g ( k ) = u ini ( k ) y ini ( k ) + σ ( k ) u ( k ) y ( k ) ,
( y ( k ) , u ( k ) ) Y N × U N .
Above
l ( y , u ) : = Q 1 2 ( y r y ) 2 2 + R 1 2 ( u r u ) 2 2 , l σ ( σ ) : = σ 2 2
for some positive definite Q , R matrices. The terminal cost is typically chosen larger than the output stage cost, to enforce convergence to the reference; a common choice is a scaled version of the output stage cost, i.e., l N ( y ) : = α l ( y , 0 ) , α 1 . The references r y R n y and r u R n u can be constant or time–varying. We assume that the sets Y and U contain r y and r u in their interior, respectively. The cost
l g ( g ) : = ( I Π ) g 2 2 , Π : = U p Y p U f U p Y p U f ,
is a regularization cost proposed in [12], where [ · ] denotes a generalized pseudo–inverse. Notice that using such a regularization cost requires T > T ini ( n u + n y ) + N n u in order to ensure that the matrix I Π has a sufficiently large null–space. If a shorter data length T is desired, alternatively, the regularization cost l g ( g ) : = g 2 2 can be used. However, this regularization is not consistent, as shown in [12].
In the deterministic, noise–free case, the DeePC algorithm [6] does not require the costs l g , l σ and the variables σ . We observe that the computational complexity of DeePC is dominated by the vector of variables g R T , with T ( n u + 1 ) ( T ini + N + n ) 1 . Hence, ideally one would prefer to work with the minimal value of data length T, but in the presence of noise, typically, a rather large data length T is required for accurate predictions [9,12,13].

3. Generalized Data–Driven Predictive Control

In this section we develop a novel, generalized DPC algorithm by constructing the predicted input sequence u ( k ) as the sum of two input sequences, i.e.,
u ( k ) : = u ¯ ( k ) + u g ( k ) , k N ,
where u ¯ is a known, base line input sequence typically chosen as the shifted, optimal input sequence from the previous time, i.e.,
u ¯ ( k ) : = { u * ( 1 | k 1 ) , , u * ( N 1 | k 1 ) , u ¯ ( N 1 | k ) } ,
where common choices for the last element u ¯ ( N 1 | k ) are r u or u * ( N 1 | k 1 ) . At time k = T ini , when enough input–output data is available to run the GDPC algorithm, u ¯ ( k ) is initialized using a zero input sequence (or an educated guess). The sequence of inputs u g can be freely optimized online by solving a QP, as explained next.
Using a single persistently exciting input sequence split into two parts, or two different persistently exciting input sequences, we can define two Hankel data matrices as in (3), i.e.,
H ¯ : = U ¯ p Y ¯ p U ¯ f Y ¯ f R ( n u + n y ) ( T ini + N ) × T , H : = U p Y p U f Y f R ( n u + n y ) ( T ini + N ) × T g ,
where the length T of the first part/sequence can be taken as large as desired and the choice of the length T g of the second part/sequence is flexible. That is, T g should be small enough to meet computational requirements, but it should provide enough degrees of freedom to optimize the bias/variance trade off. Offline, compute the matrix
Θ : = Y ¯ f U ¯ p Y ¯ p U ¯ f R ( n y N ) × ( T ini ( n u + n y ) + n u N ) .
Online, at time k T ini , given u ini ( k ) , y ini ( k ) and u ¯ ( k ) , compute
y ¯ ( k ) = Θ u ini ( k ) y ini ( k ) u ¯ ( k )
and solve the GDPC optimization problem:
min g ( k ) , u ( k ) , y ( k ) , σ ( k ) l N ( y ( N | k ) ) + i = 0 N 1 l ( y ( i | k ) , u ( i | k ) ) + λ g l g ( g ( k ) ) + λ σ l σ ( σ ( k ) ) subject to constraints :
U p Y p U f Y f g ( k ) = 0 σ ( k ) u ( k ) u ¯ ( k ) y ( k ) y ¯ ( k ) ,
( y ( k ) , u ( k ) ) Y N × U N .
Above, the cost functions l N ( y ) , l ( y , u ) , l g ( g ) and l σ ( σ ) are defined in the same way as in (5) for some Q , R 0 .
Since the size of the matrix Θ does not depend on the data length T, i.e., the number of columns of the Hankel matrix H ¯ , the computation as in (11) of the predicted output corresponding to u ¯ ( k ) is efficient even for a large T. Thus, GDPC benefits from an unbiased base line output prediction, which allows choosing T g , i.e., the number of columns of the Hankel matrix H, much smaller than T. In turn, this reduces the online computational complexity of GDPC, without sacrificing performance in the presence of noisy data. It can be argued that the selection of T g provides a trade–off between computational complexity and available degrees of freedom to optimize the bias/variance trade off.
Remark 1
(Offset–free GDPC design). In practice it is of interest to achieve offset–free tracking. Following the offset–free design for SPC developed in [21], which was further applied to DeePC in [13], it is possible to design an offset–free GDPC algorithm by defining an incremental input sequence
Δ u ( k ) : = Δ u ¯ ( k ) + Δ u g ( k ) , k N ,
where Δ u ¯ is chosen as the shifted optimal input sequence from the previous time, i.e.,
Δ u ¯ ( k ) : = { Δ u * ( 1 | k 1 ) , , Δ u * ( N 1 | k 1 ) , Δ u ¯ ( N 1 | k ) } .
The input data blocks in the Hankel matrices H ¯ and H must be replaced with incremental input data, i.e., Δ U ¯ p , Δ U ¯ f and Δ U p , Δ U f , respectively. The input applied to the system is then u ( k ) : = Δ u * ( 0 | k ) + u ( k 1 ) .
In what follows we provide a formal analysis of the GDPC algorithm.

3.1. Well-Posedness and Design of GDPC

In this subsection we show that in the deterministic case GDPC predicted trajectories are trajectories of system (1). In this case, the GDPC optimization problem can be simplified as
min g ( k ) , u ( k ) , y ( k ) l N ( y ( N | k ) ) + i = 0 N 1 l ( y ( i | k ) , u ( i | k ) ) subject to constraints :
U p Y p U f Y f g ( k ) = 0 0 u ( k ) u ¯ ( k ) y ( k ) y ¯ ( k ) ,
( y ( k ) , u ( k ) ) Y N × U N .
Lemma 1
(GDPC well–posedness). Consider one (or two) persistently exciting input sequence(s) of sufficient length(s) and construct two Hankel matrices H ¯ and H as in (9) with T and T g columns, respectively, and such that H ¯ has full row rank. Consider also the corresponding output sequence(s) generated using system (1). For any given input sequence u ¯ ( k ) , and initial conditions u ini ( k ) and y ini ( k ) , let y ¯ ( k ) be defined as in (11). Then there exists a real vector g ( k ) R T g such that (14b) holds if and only if u ( k ) and y ( k ) are trajectories of system (1).
Proof. 
Define g ¯ ( k ) : = U ¯ p Y ¯ p U ¯ f u ini ( k ) y ini ( k ) u ¯ ( k ) . Then it holds that:
H ¯ H g ¯ ( k ) g ( k ) = U ¯ p Y ¯ p U ¯ f Y ¯ f g ¯ ( k ) + U p Y p U f Y f g ( k ) = u ini ( k ) y ini ( k ) u ¯ ( k ) y ¯ ( k ) + 0 0 u ( k ) u ¯ ( k ) y ( k ) y ¯ ( k ) = u ini ( k ) y ini ( k ) u ( k ) y ( k ) .
Since the matrix H ¯ has full row rank, the concatenated matrix H ¯ H has full row rank and as such, the claim follows from [6] if one input sequence is used and from [20] if two different input sequences are used to build the Hankel matrices. □
The selection of T g enables a trade off between computational complexity and available degrees of freedom to improve the output sequence generated by the known, base line input sequence. Indeed, a larger T g results in a larger null space of the data matrix U p Y p , which confines g ( k ) in the deterministic case. However, high performance can be achieved in the case of noisy data even for a smaller T g , because the base line predicted output is calculated using an unbiased least squares predictor, i.e., as defined in (11).
An alternative way to define the known input sequence u ¯ ( k ) is to use an unconstrained SPC control law [5]. To this end, notice that the matrix Θ can be partitioned, see, e.g., [21], into P 1 P 2 Γ such that
y ¯ ( k ) = P 1 P 2 u ini ( k ) y ini ( k ) + Γ u ¯ ( k ) .
Then, by defining Ψ : = diag { R , , R } , Ω : = diag { Q , , Q , α Q } , G : = 2 Ψ + Γ T Ω Γ and F : = 2 Γ T Ω , we obtain:
u ¯ spc ( k ) : = G 1 F P 1 P 2 u ini ( k ) y ini ( k ) r y 2 Ψ r u ,
where r y : = col ( r y , , r y ) and r u : = col ( r u , , r u ) . Since the inverse of G is computed offline, computing u ¯ spc ( k ) online is numerically efficient even for a large prediction horizon N. In this case, since the corresponding base line predicted output trajectory is unbiased, the simpler regularization cost l g ( g ) : = g 2 2 can be used in (12a), without loosing consistency. When the base line input sequence is computed as in (16), the optimized input sequence u g ( k ) acts to enforce constraints, when the analytic SPC control law (16) results in violation of input or state constraints, and it can also optimize the bias/variance trade off under appropriate tuning of  λ g .
Remark 2
(Relation with SPC and DeePC). The main novelty of GDPC with respect to existing data–driven predictive controllers SPC [5] and DeePC [6] is the construction of a hybrid data-driven predictor, i.e., a combination of subspace and Hankel predictors. This comes with the benefit that online computational complexity can be reduced, by reducing the size of the Hankel matrix, while robust control performance is still achieved due to the unbiased subspace predictor. Moreover, GDPC can recover SPC or DeePC under specific settings as follows. First, notice that for a zero baseline input sequence, i.e., u ¯ ( k ) = 0 for all k N , and a sufficiently large T g , the optimization problem (14) becomes the DeePC problem (4) and yields the corresponding g DeePC * ( k ) . Then, if in GDPC one selects T g = T and H = H ¯ in (9), we have that g ( k ) : = g DeePC * ( k ) g ¯ ( k ) satisfies (15), and the DeePC solution is recovered for any baseline sequence u ¯ ( k ) . Alternatively, if the baseline sequence is chosen as in (16) and the regularization cost λ g l g ( g ( k ) ) = λ g g 2 2 is added in (14a), then u * ( k ) u SPC * ( k ) as λ g , where u SPC * ( k ) denotes optimal solution of the corresponding SPC algorithm, see, e.g., [13].
Remark 3
(GDPC open challenges). A systematic method for choosing the data length T for the subspace predictor versus the data length T g for the Hankel predictor is required. This choice represents an additional degree of freedom and depends on the noise level, system order and prediction horizon, while it affects the null space of U p Y p . The data length T for the subpace predictor can be determined via standard subspace identification criteria, see, e.g., [3,22]. T g can be chosen via the singular value decomposition method proposed in [15] or as the minimal data length required by DeePC [6]. Another open problem is the tuning of the hyper-parameters λ g and λ σ . A systematic method for tuning λ g was presented for DeePC in [13]. However, the same method does not apply directly to GDPC, due to the presence of the time-varying baseline sequence.

3.2. Stability of GDPC

In this section we will provide sufficient conditions under which GDPC is asymptotically stabilizing. To this end define J ( y ( k ) , u ( k ) ) : = l N ( y ( N | k ) ) + i = 0 N 1 l ( y ( i | k ) , u ( i | k ) ) , let l N ( y ) : = l ( y , 0 ) and let y * ( k ) and u * ( k ) denote optimal trajectories at time k T ini . Given an optimal input sequence at time k T ini , i.e., u * ( k ) = u ¯ ( k ) + u g * ( k ) , define a suboptimal input sequence at time k + 1 as
u s ( k + 1 ) = u ¯ ( k + 1 ) + u g ( k + 1 ) = col ( u * ( 1 | k ) , , u * ( N 1 | k ) , u ¯ ( N | k ) ) + col ( 0 , , 0 , 0 ) ,
and let
y s ( k + 1 ) = y ¯ ( k + 1 ) = col ( y * ( 2 | k ) , , y * ( N | k ) , y ¯ ( N + 1 | k ) )
denote the corresponding suboptimal output trajectory. Note that the last output in the suboptimal output sequence satisfies:
y ¯ ( N + 1 | k ) = i = 1 T ini a i y * ( N + 1 i | k ) + i = 1 T ini b i u * ( N + 1 i | k ) .
Definition 1
(Class K functions). A function φ : R + R + belongs to class K if it is continuous, strictly increasing and φ ( 0 ) = 0 . A function φ : R + R + belongs to class K if φ K and lim s φ ( s ) = . id denotes the identity K function, i.e., id ( s ) = s .
Next, as proposed in [11], we define a non–minimal state:
x ini ( k ) : = col ( y ( k T ini ) , , y ( k 1 ) , u ( k T ini ) , , u ( k 1 ) ) ,
and the function W ( x ini ( k ) ) : = i = 1 T ini l ( y ( k i ) , u ( k i ) ) . In what follows we assume that r y = 0 and r u = 0 for simplicity of exposition. However, the same proof applies for any constant references that are compatible with an admissible steady–state.
Assumption 1
(Terminal stabilizing condition). For any admissible initial state x ini ( k ) there exists a function ρ K , with ρ < id , a prediction horizon N T ini and u ¯ ( N | k ) U such that y ¯ ( N + 1 | k ) Y and
l ( y ¯ ( N + 1 | k ) , u ¯ ( N | k ) ) ( id ρ ) l ( y ( k T ini ) , u ( k T ini ) ) 0 .
Theorem 1
(Stability of GDPC). Suppose that there exist α 1 , l , α 2 , l , α 2 , J K such that
α 1 , l ( col ( y , u ) ) l ( y , u ) α 2 , l ( col ( y , u ) ) , ( y , u ) Y × U
J ( y * ( k ) , u * ( k ) ) α 2 , J ( x ini ( k ) ) , x ini ( k ) N ,
for some proper set N with the origin in its interior. Furthermore, let Assumption 1 hold and suppose that problem (14) is feasible for all k T ini . Then system (1) in closed–loop with the GDPC algorithm that solves problem (14) is asymptotically stable.
Proof. 
As done in [11] for the DeePC algorithm, we consider the following storage function
V ( x ini ( k ) ) : = J ( u * ( k ) , y * ( k ) ) + W ( x ini ( k ) ) ,
and we will prove that it is positive definite and it satisfies a dissipation inequality. First, from (19a) and by Lemma 14 in [11] we obtain that there exist α 1 , V , α 2 , V K such that
α 1 , V ( x ini ( k ) ) V ( x ini ( k ) ) α 2 , V ( x ini ( k ) ) .
Then, define the supply function
s ( y ( k T ini ) , u ( k T ini ) ) : = l ( y ¯ ( N + 1 | k ) , u ¯ ( N | k ) ) l ( y ( k T ini ) , u ( k T ini ) ) .
By the principle of optimality, it holds that
V ( x ini ( k + 1 ) ) V ( x ini ( k ) ) = J ( y * ( k + 1 ) , u * ( k + 1 ) ) + W ( x ini ( k + 1 ) ) J ( y * ( k ) , u * ( k ) ) W ( x ini ( k ) ) J ( y s ( k + 1 ) , u s ( k + 1 ) ) J ( y * ( k ) , u * ( k ) ) + W ( x ini ( k + 1 ) ) W ( x ini ( k ) ) = l ( y ¯ ( N + 1 | k ) , u ¯ ( N | k ) ) l ( y ( 0 | k ) , u ( 0 | k ) ) + l ( y ( 0 | k ) , u ( 0 | k ) ) l ( y ( k T ini ) , u ( k T ini ) ) = s ( y ( k T ini ) , u ( k T ini ) ) .
Hence, the storage function V ( x ini ( k ) ) satisfies a dissipation inequality along closed–loop trajectories. Since by (18) the supply function satisfies
s ( y ( k T ini ) , u ( k T ini ) ) ρ l ( y ( k T ini ) , u ( k T ini ) ) ,
the claim then follows from Corollary 17 in [11].    □
Notice that condition (18) corresponds to a particular case of condition (25) employed in Corollary 17 in [11], i.e., for M = 1 .
Remark 4
(Terminal stabilizing condition). The terminal stabilizing condition (18) can be regarded as an implicit condition, i.e., by choosing the prediction horizon N sufficiently large, this condition is more likely to hold. Alternatively, it could be implemented as an explicit constraint in problem (14) by adding one more input and output at the end of the predicted input and output sequences u ( k ) , y ( k ) , respectively. This also requires including the required additional data in the corresponding H ¯ and H Hankel matrices. This yields a convex quadratically constrained QP, which can still be solved efficiently. The stabilizing condition (18) can also be used in the regularized GDPC problem (12) to enforce convergence, but in this case a soft constraint implementation is recommended to prevent infeasibility due to noisy data.
It is worth to point out that the conditions invoked in [9] imply that condition (18) holds, i.e., Assumption 1 is less conservative. Indeed, if a terminal equality constraint is imposed in problem (14), i.e., y ( N | k ) = r y , then u ¯ ( N | k ) = r u is a feasible choice at time k + 1 , which yields y ¯ ( N + 1 | k ) = r y and hence, l ( y ¯ ( N + 1 | k ) , u ¯ ( N | k ) ) = 0 . Hence, (18) trivially holds since the stage cost l ( y , u ) is positive definite and ρ < id . Alternatively, one could employ a data–driven method to compute a suitable invariant terminal set in the space of x ini , as proposed in [23]. Tractable data–driven computation of invariant sets is currently possible for ellipsoidal sets, via linear matrix inequalities, which also yields a convex quadratically constrained QP that has to be solved online.

4. Simulations Results

In this section we consider a benchmark MPC illustrative example from Section 6.4 in [4] based on the flight control of the longitudinal motion of a Boeing 747. After discretization with zero–order–hold for T s = 0.1 [ s ] we obtain a discrete–time linear model as in (1) with:
A = 0.9997 0.0038 0.0001 0.0322 0.0056 0.9648 0.7446 0.0001 0.0020 0.0097 0.9543 0.0000 0.0001 0.0005 0.0978 1.0000 , B = 0.0010 0.1000 0.0615 0.0183 0.1133 0.0586 0.0057 0.0029 , C = 1.0000 0 0 0 0 1.0000 0 7.7400 ,
with two inputs, the throttle u 1 and u 2 , the angle of the elevator, and two outputs, the longitudinal velocity and the climb rate, respectively. The inputs and outputs are constrained as follows:
U : = u R 2 : 20 20 u 20 20 Y : = y R 2 : 25 15 Y 25 15 .
The cost function of the predictive controllers is defined as in (4) and (5) using λ g = 10 5 , λ σ = 10 7 , N = 20 , T i n i = 20 , Q = 10 · I n y and R = 0.01· I n u . For GDPC, if the suboptimal input sequence is computed as in (8), using shifted optimal sequence from the previous time with u ¯ ( k + N 1 ) = u * ( k 2 + N ) , the regularization cost l g ( g ( k ) ) is defined as in (6). If the analytic SPC solution (16) is used to compute the suboptimal input sequence, then the regularization cost l g ( g ( k ) ) = λ g g ( k ) 2 2 is used. For DeePC the regularization cost l g ( g ( k ) ) is defined as in (6). In this way, all 3 compared predictive control algorithms utilize consistent output predictors.
The QP problems corresponding to the data–driven predictive controllers are solved using the quadprog Matlab QP solver on a laptop with an Intel i7-9750H CPU and 16GB of RAM, with MATLAB version R2021a. The quadprog solver is not the fastest QP solver in general, but it has a high accuracy, which is useful to analyze with precision the control performance of different algorithms. The OSQP solver [24] can provide a faster alternative, and it was used in this work only to solve quadratically constrained QPs. Such optimization problems arise if the stabilizing constraint (18) is included in the GDPC problem.
The GDPC algorithm is implemented in Matlab based on Algorithms 1 and 2.
Algorithm 1. Data generation method, prior to controlling the system
Input: 
Measurable system, Desired sampling period T s
    1:
Generate persistently exciting input U = { u ( k ) } k N [ 0 , T ] , see Section 13.3 in [22]
    2:
Apply input, measure and collect system response Y = { y ( k ) } k N [ 0 , T ]
    3:
Construct Hankel matrices H ¯ and H as in (9)
    4:
Compute Θ as in (10) using H ¯
Algorithm 2. Implementation of GDPC
  • Input:  H ¯ , H, Θ and the system constraints
  •     for  k = 0 : T i n i do
  •       1: Measure y ( k )
  •       2: Apply u ( k ) = γ · rand ( k ) to the system and wait for the next sampling interval
  •     end for
  • for  k = T i n i + 1 : end do
  •       1: Measure y ( k )
  •       2: Build y i n i ( k ) , u i n i ( k ) from the last T i n i data samples
  •       3: Compute u ¯ ( k ) using either (8) or (16)
  •       4: Compute y ¯ ( k ) = Θ u i n i ( k ) T y i n i ( k ) T u ¯ ( k ) T T
  •       5: Solve QP (12) with any QP solver to obtain u ( k )
  •       6: Apply u ( k ) = u ( k ) [ 0 : n u ] to the system and wait for the next sampling interval
  •     end for
In what follows, the simulation results are structured into 3 subsections, focusing on nominal, noise–free data performance, noisy data performance for low and high noise variance and comparison with DeePC. The controllers will only start when T ini samples have been collected. For the time instants up to k = T ini , the system is actuated by a small random input. For the sake of a sound comparison, the random input signal used up to T ini is identical for all simulations/predictive controllers. In the data generation experiment, the system inputs are constructed using two PRBS signals, both constrained between [ 3 , 3 ] .

4.1. Noise–Free Data GDPC Performance

In this simulation we implement the GDPC algorithm with the suboptimal input sequence computed as in (8). Figure 1 shows the outputs, inputs and optimized inputs u g ( k ) over time.
We observe that the GDPC closed–loop trajectories converge to the reference values and that the optimized inputs are active only at the start, after which the suboptimal shifted sequence becomes optimal. The stabilizing condition (18) is implicitly satisfied along trajectories; when imposed online, the GDPC problem is recursively feasible and yields the same trajectories.

4.2. Noisy Data GDPC Performance

Next we illustrate the performance of GDPC for noisy data with a low and high variance. Figure 2 and Figure 3 show the GDPC response using u ¯ ( k ) as defined in (8) and (16), respectively, for low variance noise.
Although the two different methods to calculate the base line input sequence u ¯ ( k ) show little difference in the resulting total input u ( k ) , the optimized part of the input, u g ( k ) shows a notable difference. For the GDPC algorithm that uses the unconstrained SPC we see that the optimized input only acts to enforce constraints, while around steady state the unconstrained SPC becomes optimal.
Next, we show the performance of GDPC for high–variance noise, which also requires a suitable increase of the data size. Figure 4 and Figure 5 show the GDPC response using u ¯ ( k ) as defined in (8) and (16), respectively, for high variance noise. We see that both GDPC algorithms achieve robust control performance, while in the case of the SPC baseline input sequence, the optimized input is active also around steady state. This shows that GDPC indeed optimizes the bias variance trade off with respect to the SPC solution.
In Figure 6 it can be observed that DeePC requires a data sequence of length 750 to achieve similar performance with GDPC with data lengths T g = 250 (relevant for online complexity), T = 5000 .
The three tested predictive controllers yield the following average computational time for the high–variance noise simulation: GDPC with u ¯ as in (8) −59 ms; GDPC with u ¯ as in (16) −30 ms; DeePC −371 ms. This suggests that for high noise and large scale systems GDPC provides a computationally efficient alternative to DeePC.

4.3. Comparison with DeePC over Multiple Runs

In this section, we compare the performance of GDPC with DeePC for different sizes of T over multiple runs. The performance can be expressed using the integral squared error (ISE) or the integral absolute error (IAE) in relation with the input energy (InEn), formally defined as:
ISE = k = T ini t max y ( k ) r ( k ) 2 2 , IAE = k = T ini t max y ( k ) r ( k ) 1 , InEn = k = T ini t max u ( k ) 2 2 ,
where t max is the simulation time and note that the performance scores are computed after the simulation ends, thus using simulated data (not predicted data). Furthermore, both the data-collecting experiment and the simulation are influenced by noise with variance σ w 2 I = 0.05I . All algorithms are designed based on the same offline collected data set and tested online using the same noise realization.
As shown in Figure 7, we notice that the GDPC predictor with data length T g = 150 (Hankel matrix based predictor) and T = 1000 (least squares based predictor) can match the performance of DeePC with data length T = 500 , while the DeePC performance with T = 150 is lower. Also, the input energy is lower for GDPC compared with the the same cost for DeePC with T = 500 . Two extreme outliers were removed from the results of DeePC with T = 150 to make the plots more clear.
From a computational complexity point of view, as shown in Table 2, the average CPU time of GDPC with T g = 150 and T = 500 is of the same order as the average CPU time of DeePC with T = 150 , while the average CPU time of DeePC with T = 500 is about 8 times higher.
The obtained results validate the fact that GDPC offers more flexibility to optimize the trade off between control performance in the presence of noisy data and online computational complexity. As such, GDPC provides engineers with a practical and robust data–driven predictive controller suitable for real–time implementation.

5. Conclusions

In this paper we developed a generalized data–driven predictive controller by merging a subspace unbiased predictor with a Hankel matrix based predictor. We have shown that this formulation results in a well–posed data–driven predictive controller with similar properties as DeePC. Sufficient conditions for closed–loop asymptotic stability of GDPC have been establsihed. Also, we have shown that compared to DeePC, the developed GDPC algorithm provides a trade–off between computational complexity and robust control performance.
To demonstrate this, Table 3 summarizes the performance versus computational complexity trade-off of the developed GDPC algorithms and DeePC, respectively, for system (20) and the following settings: N = 50 , σ w 2 I = 0.2I , T = 5000 and T g = 250 for GPDC algorithms, and for DeePC we test T = 250 and T = 500 . All algorithms have been designed using the same offline collected data set and tested online using the same noise realization.
We observe that the performance of DeePC increases / improves with the increase of the data length, but its computational complexity increases as well. For GDPC with a baseline input sequence as in (8), i.e., the shifted sequence from previous time, and with T g = 250 outperforms DeePC for the same data size T = 250 , with a slight computational time increase. Moreover, GDPC with a baseline input sequence as in (16), i.e., the unconstrained SPC control law, and with T g = 250 outperforms DeePC even with T = 500 , while it is computationally faster than DeePC with T = 250 . The fact that GDPC with baseline (16) uses more input energy than DeePC with T = 500 to achieve improved control performance, suggests that the merged predictor in this case is superior, as it is able to correctly exploit the input energy for controlling the system, despite the large noise variance in both offline and online data.
In future work we will research systematic methods for choosing the data lengths T , T g and tuning the hyper parameters λ g and λ σ . Also, we will research methods for providing robust stability guarantees for GDPC using the input-to-state stability framework.

Author Contributions

Methodology, M.L. and P.C.N.V.; Software, M.L. and P.C.N.V.; Formal analysis, M.L.; Writing original draft, M.L.; Review and editing, P.C.N.V. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hou, Z.S.; Wang, Z. From model-based control to data-driven control: Survey, classification and perspective. Inf. Sci. 2013, 235, 3–35. [Google Scholar] [CrossRef]
  2. Lamnabhi-Lagarrigue, F.; Annaswamy, A.; Engell, S.; Isaksson, A.; Khargonekar, P.; Murray, R.M.; Nijmeijer, H.; Samad, T.; Tilbury, D.; Van den Hof, P. Systems & Control for the future of humanity, research agenda: Current and future roles, impact and grand challenges. Annu. Rev. Control 2017, 43, 1–64. [Google Scholar]
  3. Maciejowski, J.M. Predictive Control with Constraints; Prentice Hall: Harlow, UK, 2002. [Google Scholar]
  4. Camacho, E.F.; Bordons, C. Model Predictive Control; Springer: Berlin/Heidelberg, Germany, 2007. [Google Scholar]
  5. Favoreel, W.; De Moor, B.; Gevers, M. SPC: Subspace predictive control. In Proceedings of the 14th IFAC World Congress, Beijing, China, 5–9 July 1999; Volume 32, pp. 4004–4009. [Google Scholar]
  6. Coulson, J.; Lygeros, J.; Dörfler, F. Data-Enabled Predictive Control: In the Shallows of the DeePC. In Proceedings of the 18th European Control Conference, Naples, Italy, 25–28 June 2019; pp. 307–312. [Google Scholar]
  7. Willems, J.C.; Rapisarda, P.; Markovsky, I.; De Moor, B.L. A note on persistency of excitation. Syst. Control. Lett. 2005, 54, 325–329. [Google Scholar] [CrossRef]
  8. Yang, H.; Li, S. A data–driven predictive controller design based on reduced Hankel matrix. In Proceedings of the 10th Asian Control Conference (ASCC), Kota Kinabalu, Malaysia, 31 May–3 June 2015; pp. 1–7. [Google Scholar]
  9. Berberich, J.; Köhler, J.; Müller, M.A.; Allgöwer, F. Data-Driven Model Predictive Control With Stability and Robustness Guarantees. IEEE Trans. Autom. Control 2021, 66, 1702–1717. [Google Scholar] [CrossRef]
  10. Fiedler, F.; Lucia, S. On the relationship between data–enabled predictive control and subspace predictive control. In Proceedings of the 2021 European Control Conference (ECC), Delft, The Netherlands, 29 June–2 July 2021; pp. 222–229. [Google Scholar]
  11. Lazar, M. A Dissipativity-Based Framework for Analyzing Stability of Predictive Controllers. IFAC PapersOnLine 2021, 54, 159–165. [Google Scholar] [CrossRef]
  12. Dorfler, F.; Coulson, J.; Markovsky, I. Bridging direct & indirect data-driven control formulations via regularizations and relaxations. IEEE Trans. Autom. Control 2022, 68, 883–897. [Google Scholar] [CrossRef]
  13. Lazar, M.; Verheijen, P.C.N. Offset–free data–driven predictive control. In Proceedings of the 2022 IEEE 61st Conference on Decision and Control (CDC), Cancun, Mexico, 6–9 December 2022; pp. 1099–1104. [Google Scholar] [CrossRef]
  14. Breschi, V.; Chiuso, A.; Formentin, S. Data-driven predictive control in a stochastic setting: A unified framework. Automatica 2023, 152, 110961. [Google Scholar] [CrossRef]
  15. Zhang, K.; Zheng, Y.; Li, Z. Dimension Reduction for Efficient Data-Enabled Predictive Control. arXiv 2022. [Google Scholar] [CrossRef]
  16. Baros, S.; Chang, C.Y.; Colón-Reyes, G.E.; Bernstein, A. Online data-enabled predictive control. Automatica 2022, 138, 109926. [Google Scholar] [CrossRef]
  17. Elokda, E.; Coulson, J.; Beuchat, P.N.; Lygeros, J.; Dörfler, F. Data-enabled predictive control for quadcopters. Int. J. Robust Nonlinear Control 2021, 31, 8916–8936. [Google Scholar] [CrossRef] [PubMed]
  18. Berberich, J.; Köhler, J.; Müller, M.A.; Allgöwer, F. Data-driven model predictive control: Closed-loop guarantees and experimental results. at—Automatisierungstechnik 2021, 69, 608–618. [Google Scholar] [CrossRef]
  19. Carlet, P.G.; Favato, A.; Bolognani, S.; Dörfler, F. Data-Driven Continuous-Set Predictive Current Control for Synchronous Motor Drives. IEEE Trans. Power Electron. 2022, 37, 6637–6646. [Google Scholar] [CrossRef]
  20. van Waarde, H.J.; De Persis, C.; Camlibel, M.K.; Tesi, P. Willems’ Fundamental Lemma for State-Space Systems and Its Extension to Multiple Datasets. IEEE Control Syst. Lett. 2020, 4, 602–607. [Google Scholar] [CrossRef]
  21. Verheijen, P.C.N.; Gonçalves da Silva, G.R.; Lazar, M. Data–driven rate–based integral predictive control with estimated prediction matrices. In Proceedings of the 25th International Conference on System Theory, Control and Computing (ICSTCC), Iasi, Romania, 20–23 October 2021; pp. 630–636. [Google Scholar]
  22. Ljung, L. System Identification: Theory for the User, 2nd ed.; Prentice Hall PTR: Upper Saddle River, NJ, USA, 1999. [Google Scholar]
  23. Berberich, J.; Köhler, J.; Müller, M.A.; Allgöwer, F. On the design of terminal ingredients for data-driven MPC. IFAC–PapersOnLine 2021, 54, 257–263. [Google Scholar] [CrossRef]
  24. Stellato, B.; Banjac, G.; Goulart, P.; Bemporad, A.; Boyd, S. OSQP: An operator splitting solver for quadratic programs. Math. Program. Comput. 2020, 12, 637–672. [Google Scholar] [CrossRef]
Figure 1. GDPC tracking performance: T g = 150 , T = 1000 , σ w 2 I = 0 I .
Figure 1. GDPC tracking performance: T g = 150 , T = 1000 , σ w 2 I = 0 I .
Mathematics 11 02216 g001
Figure 2. GDPC tracking performance: T g = 150 , T = 1000 , σ w 2 I = 0.05I , u ¯ as in (8).
Figure 2. GDPC tracking performance: T g = 150 , T = 1000 , σ w 2 I = 0.05I , u ¯ as in (8).
Mathematics 11 02216 g002
Figure 3. GDPC tracking performance: T g = 150 , T = 1000 , σ w 2 I = 0.05I , u ¯ as in (16).
Figure 3. GDPC tracking performance: T g = 150 , T = 1000 , σ w 2 I = 0.05I , u ¯ as in (16).
Mathematics 11 02216 g003
Figure 4. GDPC tracking performance: T g = 250 , T = 5000 , σ w 2 I = 0.5I , u ¯ as in (8).
Figure 4. GDPC tracking performance: T g = 250 , T = 5000 , σ w 2 I = 0.5I , u ¯ as in (8).
Mathematics 11 02216 g004
Figure 5. GDPC tracking performance: T g = 250 , T = 5000 , σ w 2 I = 0.5I , u ¯ as in (16).
Figure 5. GDPC tracking performance: T g = 250 , T = 5000 , σ w 2 I = 0.5I , u ¯ as in (16).
Mathematics 11 02216 g005
Figure 6. DeePC tracking performance: T = 750 , σ w 2 I = 0.5I .
Figure 6. DeePC tracking performance: T = 750 , σ w 2 I = 0.5I .
Mathematics 11 02216 g006
Figure 7. 1: GDPC with u ¯ as in (8), T = 1000 , T g = 150 ; 2: DeePC with T = 150 ; 3: DeePC with T = 500 . σ w 2 I = 0.05I , using 50 Monte-Carlo runs.
Figure 7. 1: GDPC with u ¯ as in (8), T = 1000 , T g = 150 ; 2: DeePC with T = 150 ; 3: DeePC with T = 500 . σ w 2 I = 0.05I , using 50 Monte-Carlo runs.
Mathematics 11 02216 g007
Table 1. Table of abbreviations.
Table 1. Table of abbreviations.
MPCModel Predictive Control
DPCData–driven Predictive Control
SPCSubspace Predictive Control
DeePCData–enabled Predictive Control
GDPCGeneralized Data–driven Predictive Control
QPQuadratic Programming
Table 2. Average CPU time (Quadprog) over all runs for various data–driven predictive controllers.
Table 2. Average CPU time (Quadprog) over all runs for various data–driven predictive controllers.
GDPC (8), T g = 150 DeePC T = 150 DeePC T = 500
CPU33 ms32 ms260 ms
Table 3. Average CPU time (Quadprog) and performance indicators for various DPC algorithms.
Table 3. Average CPU time (Quadprog) and performance indicators for various DPC algorithms.
GDPC with (8)GDPC with (16)DeePC T = 250 DeePC T = 500
CPU131 ms68 ms106 ms386 ms
ISE290258397265
IAE867415975
InEn3400440039003400
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Lazar, M.; Verheijen, P.C.N. Generalized Data–Driven Predictive Control: Merging Subspace and Hankel Predictors. Mathematics 2023, 11, 2216. https://doi.org/10.3390/math11092216

AMA Style

Lazar M, Verheijen PCN. Generalized Data–Driven Predictive Control: Merging Subspace and Hankel Predictors. Mathematics. 2023; 11(9):2216. https://doi.org/10.3390/math11092216

Chicago/Turabian Style

Lazar, M., and P. C. N. Verheijen. 2023. "Generalized Data–Driven Predictive Control: Merging Subspace and Hankel Predictors" Mathematics 11, no. 9: 2216. https://doi.org/10.3390/math11092216

APA Style

Lazar, M., & Verheijen, P. C. N. (2023). Generalized Data–Driven Predictive Control: Merging Subspace and Hankel Predictors. Mathematics, 11(9), 2216. https://doi.org/10.3390/math11092216

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop