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
:
where
is the state,
is the control input,
is the measured output and
are real matrices of suitable dimensions. We assume that
is controllable and
is observable. By applying a persistently exciting input sequence
of length
T to system (
1) we obtain a corresponding output sequence
.
If one considers an input–output model corresponding to (
1), it is necessary to introduce the parameter
that limits the window of past input–output data necessary to compute the current output, i.e.,
for some real–valued coefficients. For simplicity of exposition we assume the same
for inputs and outputs.
Next, we introduce some instrumental notation. For any finite number
of vectors
we will make use of the operator
. For any
(starting time instant in the data vector) and
(length of the data vector), define
Let
denote the prediction horizon. Then we can define the Hankel data matrices:
According to the DeePC design [
6], one must choose
and
, which implicitly requires an assumption on the system order (number of states). Given the measured output
at time
and
we define the sequences of known input–output data at time
, which are trajectories of system (
1):
Next, we define the sequences of predicted inputs and outputs at time
, which should also be trajectories of system (
1):
where the notation
denotes the predicted value of
using measured output data available up to
.
For a positive definite matrix
L let
denote its Cholesky factorization. At time
, given
, the regularized DeePC algorithm [
12] computes a sequence of predicted inputs and outputs as follows:
Above
for some positive definite
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.,
,
. The references
and
can be constant or time–varying. We assume that the sets
and
contain
and
in their interior, respectively. The cost
is a regularization cost proposed in [
12], where
denotes a generalized pseudo–inverse. Notice that using such a regularization cost requires
in order to ensure that the matrix
has a sufficiently large null–space. If a shorter data length
T is desired, alternatively, the regularization cost
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
and the variables
. We observe that the computational complexity of DeePC is dominated by the vector of variables
, with
. 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
as the sum of two input sequences, i.e.,
where
is a known, base line input sequence typically chosen as the shifted, optimal input sequence from the previous time, i.e.,
where common choices for the last element
are
or
. At time
, when enough input–output data is available to run the GDPC algorithm,
is initialized using a zero input sequence (or an educated guess). The sequence of inputs
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.,
where the length
T of the first part/sequence can be taken as large as desired and the choice of the length
of the second part/sequence is flexible. That is,
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
Online, at time
, given
,
and
, compute
and solve the
GDPC optimization problem:
Above, the cost functions
,
,
and
are defined in the same way as in (
5) for some
.
Since the size of the matrix
does not depend on the data length
T, i.e., the number of columns of the Hankel matrix
, the computation as in (
11) of the predicted output corresponding to
is efficient even for a large
T. Thus, GDPC benefits from an unbiased base line output prediction, which allows choosing
, 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
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 sequencewhere is chosen as the shifted optimal input sequence from the previous time, i.e.,The input data blocks in the Hankel matrices and H must be replaced with incremental input data, i.e., , and , , respectively. The input applied to the system is then . 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
Lemma 1 (GDPC well–posedness).
Consider one (or two) persistently exciting input sequence(s) of sufficient length(s) and construct two Hankel matrices and H as in (9) with T and columns, respectively, and such that has full row rank. Consider also the corresponding output sequence(s) generated using system (1). For any given input sequence , and initial conditions and , let be defined as in (11). Then there exists a real vector such that (14b) holds if and only if and are trajectories of system (1). Proof. Define
. Then it holds that:
Since the matrix
has full row rank, the concatenated matrix
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
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
results in a larger null space of the data matrix
, which confines
in the deterministic case. However, high performance can be achieved in the case of noisy data even for a smaller
, 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
is to use an unconstrained SPC control law [
5]. To this end, notice that the matrix
can be partitioned, see, e.g., [
21], into
such that
Then, by defining
,
,
and
, we obtain:
where
and
. Since the inverse of
G is computed offline, computing
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
can be used in (
12a), without loosing consistency. When the base line input sequence is computed as in (
16), the optimized input sequence
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
.
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., for all , and a sufficiently large , the optimization problem (14) becomes the DeePC problem (4) and yields the corresponding . Then, if in GDPC one selects and in (9), we have that satisfies (15), and the DeePC solution is recovered for any baseline sequence . Alternatively, if the baseline sequence is chosen as in (16) and the regularization cost is added in (14a), then as , where 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 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 . The data length T for the subpace predictor can be determined via standard subspace identification criteria, see, e.g., [3,22]. 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 and . A systematic method for tuning 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
, let
and let
and
denote optimal trajectories at time
. Given an optimal input sequence at time
, i.e.,
, define a suboptimal input sequence at time
as
and let
denote the corresponding suboptimal output trajectory. Note that the last output in the suboptimal output sequence satisfies:
Definition 1 (Class functions). A function belongs to class if it is continuous, strictly increasing and . A function belongs to class if and . id denotes the identity function, i.e., .
Next, as proposed in [
11], we define a non–minimal state:
and the function
. In what follows we assume that
and
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 there exists a function , with , a prediction horizon and such that and Theorem 1 (Stability of GDPC).
Suppose that there exist such thatfor some proper set with the origin in its interior. Furthermore, let Assumption 1 hold and suppose that problem (14) is feasible for all . 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
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
such that
Then, define the supply function
By the principle of optimality, it holds that
Hence, the storage function
satisfies a dissipation inequality along closed–loop trajectories. Since by (
18) the supply function satisfies
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
.
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 , , respectively. This also requires including the required additional data in the corresponding 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.,
, then
is a feasible choice at time
, which yields
and hence,
. Hence, (
18) trivially holds since the stage cost
is positive definite and
. Alternatively, one could employ a data–driven method to compute a suitable invariant terminal set in the space of
, 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
we obtain a discrete–time linear model as in (
1) with:
with two inputs, the throttle
and
, the angle of the elevator, and two outputs, the longitudinal velocity and the climb rate, respectively. The inputs and outputs are constrained as follows:
The cost function of the predictive controllers is defined as in (4) and (
5) using
,
,
,
,
and
. For GDPC, if the suboptimal input sequence is computed as in (
8), using shifted optimal sequence from the previous time with
, the regularization cost
is defined as in (
6). If the analytic SPC solution (
16) is used to compute the suboptimal input sequence, then the regularization cost
is used. For DeePC the regularization cost
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 - 1:
Generate persistently exciting input , see Section 13.3 in [ 22] - 2:
Apply input, measure and collect system response - 3:
Construct Hankel matrices and H as in ( 9) - 4:
Compute as in ( 10) using
|
Algorithm 2. Implementation of GDPC |
Input: , H, and the system constraints for
do 1: Measure 2: Apply to the system and wait for the next sampling interval end for for do 1: Measure 2: Build , from the last data samples 3: Compute using either ( 8) or ( 16) 4: Compute 5: Solve QP (12) with any QP solver to obtain 6: Apply 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 samples have been collected. For the time instants up to , the system is actuated by a small random input. For the sake of a sound comparison, the random input signal used up to is identical for all simulations/predictive controllers. In the data generation experiment, the system inputs are constructed using two PRBS signals, both constrained between .
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
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
as defined in (
8) and (
16), respectively, for low variance noise.
Although the two different methods to calculate the base line input sequence show little difference in the resulting total input , the optimized part of the input, 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
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
(relevant for online complexity),
.
The three tested predictive controllers yield the following average computational time for the high–variance noise simulation: GDPC with
as in (
8) −59 ms; GDPC with
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:
where
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
. 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
(Hankel matrix based predictor) and
(least squares based predictor) can match the performance of DeePC with data length
, while the DeePC performance with
is lower. Also, the input energy is lower for GDPC compared with the the same cost for DeePC with
. Two extreme outliers were removed from the results of DeePC with
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
and
is of the same order as the average CPU time of DeePC with
, while the average CPU time of DeePC with
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:
,
,
and
for GPDC algorithms, and for DeePC we test
and
. 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
outperforms DeePC for the same data size
, 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
outperforms DeePC even with
, while it is computationally faster than DeePC with
. The fact that GDPC with baseline (
16) uses more input energy than DeePC with
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 and tuning the hyper parameters and . Also, we will research methods for providing robust stability guarantees for GDPC using the input-to-state stability framework.