Towards explainable data-driven predictive control with regularizations

2.1 Fundamentals of direct data-driven predictions

Instead of utilizing a discrete-time state-space model with input u ∈ R m , state x ∈ R n , and output y ∈ R p as in traditional MPC, predictions in DPC are realized based on previously collected trajectory data u ( 1 ) , y ( 1 ) , … , u ( ℓ ) , y ( ℓ ) ∈ R m L × R p L via linear combinations

Here, the dimensions of the data matrix D ∈ R L ( m + p ) × ℓ and generator vector a ∈ R ℓ are specified by the length L of recorded (and generated) trajectories and the number ℓ of data trajectories used for predictions. Furthermore, we refer to a pair u , y as an input/output-sequence (I/O-sequence), and use roman notation to distinguish individual variables like u , y from their associated sequences u, y. Assuming exact data generated by an LTI system, which can be described, e.g., by the state-space representation

and assuming L is greater than the lag ι of the system, the image R D is equivalent to the set of all possible system trajectories (of length L) if and only if [15]

Note that the lag ι is an integer invariant [5] of the system, i.e., invariant with respect to the considered representation. In context of the state-space representation (2), it is also known as its observability index, which is defined as the smallest integer ι, for which r a n k O ι = r a n k O ι + 1 with

is satisfied. The generalized persistency of excitation condition (3) provides the theoretical foundation for the linear combinations (1), since

Representing system trajectories in this way is also known as an image representation of the system, as opposed to, e.g., a state-space representation (2). A popular sufficient condition for data to satisfy (3) is known as Willems’ Fundamental Lemma [5], which has become synonymous with using image representations. We note that, contrary to the Fundamental Lemma, the generalized persistency of excitation condition (3) is both sufficient and necessary, and neither requires controllability nor a Hankel structure for the data matrix. To include the current initial condition of the system as a starting point for predicted trajectories, the generated I/O-sequence is typically partitioned into a past section (u _p , y _p ) and a future section (u _f , y _f ) with N _p respectively N _f time-steps yielding

The past section of a predicted trajectory is then forced to match the I/O-data ξ recorded in the most recent N _p time-steps during closed-loop operation, i.e., the constraints

force any predicted trajectory to start with the most recently witnessed behavior of the system. Note that ξ is a (non-minimal) state of the LTI system (2), if N _p is chosen greater or equal to its lag ι. From now on, we omit the “future” subscript from u _f , y _f , U _f , and Y _f , since their “past” counterparts are already incorporated in W and ξ , eliminating any risk of confusion. For more concise notation, we also define

Additionally, with a slight abuse of notation, we redefine D ≔ Z ⊤   Y ⊤ ⊤ . The latter is simply a block-row permuted version of the original data matrix introduced in (1), used to facilitate the corresponding partitioning of u_gen, y_gen into ξ , u, y. To summarize these new notations, W contains data columns for the state ξ , U contains data columns for the future input sequence u, and Y contains data columns for the future output sequence y. Furthermore, Z contains data for the state-input sequence z , and D contains data for the state-input-output sequence w .

Notably, in the ideal deterministic LTI setting with condition (3), the data matrix always has the rank deficiency

such that the image representation (4) implies a unique (and exact) linear predictor mapping y = y ̂ ( ξ , u ) . However, this rank deficiency, and with it the unique predictions, are typically lost in the presence of noise or nonlinearities, which is visualized in Figure 2. One remedy for this might be the use of (structured) low-rank approximations [17] for D , which sometimes even implicitly occurs in DPC schemes, e.g., in γ -DDPC [18] with γ ₃ = 0 (see [13]). Specific to the context of DPC, another remedy is given by the addition of a regularization term h( a ) to the objective function (6a) [1], [2], which is in the focus of this work.

Figure 2:

Visualization of the data matrix D for the low-dimensional state-space example in Section 2.3. The green marks show data columns x 0 ( i ) , u ( i ) , x ( i ) for (a) the ideal deterministic LTI setting and (b) the same data with added measurement noise. The gray subspaces depict R D , showing r a n k D = r a n k Z = 2 in (a) and full rank D with r a n k D = 3 > r a n k Z = 2 in (b).

2.2 Regularized DPC

The optimal control problem (OCP) that is solved for DPC in every time step can be stated as

with control objective J( ξ , u, y), regularization h( a ), and input-output constraints U × Y . Conditions for equivalence of DPC and MPC are well established for some special cases. In particular, as mentioned in the introduction, the equivalence with MPC for the LTI system (2) holds with exact data satisfying (3) and setting the regularization to h( a ) = 0 [2]. In the presence of noise or nonlinearities, the unregularized OCP may use the (unrealistic) additional degrees of freedom, available to it by (5) not being satisfied, to greedily minimize the objective function J( ξ , u, y). Given enough data (columns), it is established and realistic (see, e.g. [12], [18], [19], for a discussion) to make the following assumption.

Note that full row-rank of D renders (6b) meaningless without regularization, since there is an a solving (6b) for any arbitrary left-hand side trajectory.

2.3 Running numerical example

We emphasize that this work does not propose a new DPC scheme; rather, it provides tools for analyzing the structure of the OCPs in existing schemes. Therefore, and in light of the following discussion in Section 3 regarding point (i), we do not consider extensive closed-loop simulations to be the best demonstration of our results. Such simulations can be found in the works we reference when analyzing the respective schemes. Instead, we use a low-dimensional system as a running example to visualize the structure of the data matrix D itself in Figure 2, and the structure of implicit predictors in Figures 3 and 4 (see page 373, page 378). Since using a past I/O sequence as a non-minimal state ξ yields a minimum dimension of (p + m)N _p ≥ 2, we instead treat the example system in a state-space setting (see Remark 1). The data-generating system is LTI with dimensions n = m = 1 and parameters ( A , B ) = (2, −0.5). To generate Figure 2(a), we drew 500 i.i.d. samples x 0 ( i ) , u ( i ) ∼ U [ − 1,1 ] and computed the subsequent state x⁽ⁱ⁾ via (2). For Figure 2(b), we added i.i.d. measurement noise N ( 0,0.01 ) to each x 0 ( i ) , x ( i ) . To visualize the implicit predictors and optimal parametric solutions in Figure 3, we chose N _f = 1 and set up the control objective (18) with reference x_ref = 0, weight Q = 1 , and quadratic input costs J _u ( ξ , u) = u². The data matrix D is a block Hankel matrix (see, e.g., [5]), constructed from the input-state sequence

which was generated by drawing an initial state x 0 ∼ U [ − 1,1 ] , applying the i.i.d. input sequence u d ( i ) ∼ U [ − 1,1 ] via (2), and then adding i.i.d. measurement noise N ( 0,0.01 ) to the resulting state sequence. To visualize the implicit predictors in Figure 4, we set N _f = 2, with reference x ref = 0.5 0.5 ⊤ , weight Q = I 2 , and arbitrary input costs, using a regularization weight of λ = 1. The data matrix D is again a block Hankel matrix, constructed from the input-state sequence

where the two additional samples are generated in the same way as before. These extra samples are needed to allow that the data matrix D satisfies Assumption 1 for the increased prediction horizon N _f = 2.

4.1 Trajectory-specific effect of quadratic regularization

It was first shown in [12] with additional details in [13], Prop. 1] that, under Assumption 1, the trajectory-specific effect of quadratic regularization h ( a ) = λ ‖ a ‖ 2 2 is given by

Here, y ̂ LS ( ξ , u ) and u ̂ LS ( ξ ) are the (multistep) predictor/controller “favored” by the regularization, as initially mentioned in Section 3. The regularization pushes predicted y and u towards them, which is visualized in Figure 3 for the low-dimensional state-space example discussed in Section 2.3. By parametrically solving (1) and decomposing block-matrix expressions (see [12] for details), one can show that both are given by linear mappings, which can be represented as

via the least-squares solutions

where ‖.‖_F denotes the Frobenius norm. Note that y ̂ LS ( ξ , u ) is equivalent to the subspace predictor used as an equality constraint in subspace predictive control (SPC, [24]). The involved weighing matrices are given by

where

are the residual matrices associated with the least-square problems (13), and the inverses exist under Assumption 1. Therefore, Q reg and R reg can be interpreted as inverses of scaled (because they are not normalized w.r.t. the amount of data columns ℓ) empirical second moment matrices for the output prediction error E _y and “input prediction error” E _u of the least-squares solutions G _LS, K _LS based on the data D . Hence, quadratic regularization h ( a ) = λ ‖ a ‖ 2 2 pushes y and u towards the least-squares estimates y ̂ LS ( ξ , u ) and u ̂ LS ( ξ ) . It does so more (less) harshly in directions, where the latter are believed to be more (less) accurate, based on the available data. The same second moment (or covariance)-based interpretation appears in [19], with the distinction that the state-input dependent part in [19], Eq. (21)] is not yet decomposed into the two terms (10b) and (10c), similar to [12]. However [19], also shows that the results naturally generalize to rank-deficient data matrices, i.e., D not satisfying Assumption 1, by simply using the Moore–Penrose inverse in (14). In that case, deviations e y = y − y ̂ LS ( ξ , u ) are not penalized in singular directions of E y E y ⊤ , i.e., directions in which no deviation from y ̂ LS ( ξ , u ) has been observed in the data matrix D . However, due to (7), one can show that y may not deviate at all in these directions. That is, in unexplored directions of the state-input-output space, the penalization of deviations in (10) is replaced with hard constraints, since deviations from it are not possible in R D . The same considerations apply for deviations e u = u − u ̂ LS ( ξ ) if Z does not have full row-rank.

When discussing the role of these cost terms, first note that the last cost term (10c) is irrelevant to the OCP, since ξ is a parameter determined in closed-loop and not an optimization variable. Regarding the usefulness of (10a), we believe that y being pushed towards the least-squares (multistep) predictor y ̂ LS ( ξ , u ) in (10a) is quite intuitive. However, while y ̂ LS ( ξ , u ) is conceptually not a bad choice for a predictor, it should be noted that (without further modifications such as in [25]) it is typically neither causal nor time-invariant. Furthermore, although it is conceptually nice to include additional statistical information via Q reg and R reg , this information is utilized very greedily in DPC. That is, since the optimal solution is a trade-off between control objective and regularization cost, the DPC cost objective (6a) promotes deviations from y ̂ LS ( ξ , u ) primarily in uncertain directions in order to achieve lower costs for the control objective J( ξ , u, y). We believe that it remains an interesting avenue for future research to explore under what conditions penalizing predictions y towards a (potentially well- or ill-suited) predictor y ̂ ( ξ , u ) , as in (10a), outperforms enforcing it as an equality constraint y = y ̂ ( ξ , u ) . While the intuition behind penalizing the input u towards u ̂ LS ( ξ ) in the same fashion may not be immediately clear, it was shown in the γ -DDPC setting that a cost term equivalent to (10b) (proven in [12]) can indeed be utilized to increase DPC performance [26]. Furthermore, a functionally similar input penalty appears in another framework focused on minimizing the final control error (FCE) [27]. Intuitively, the performance improvement can be explained by the term (10b) pushing the chosen input u towards the best-explored region of the state-input-space, i.e., where the most confident predictions can be made based on the available data. For a more precise quantification of the relation between prediction error variance and ‖ γ 2 ‖ 2 2 , we refer to [26].

Finally, we want to highlight the discrepancy between (10a) and (10b) in terms of proper tuning for λ. The output-related cost term (10a) needs a large weight λ because, without the previously discussed rank deficiency (5), it is the sole factor keeping output predictions y from being greedily and unrealistically (i.e., without considering the data in D ) pushed toward the minimum of the control objective J( ξ , u, y). This effect is visualized in Figure 3(a). On the other hand, the input-related cost term (10b) may not have a large weight λ, since this would favor always sticking to the “best-explored” choice u ̂ LS ( ξ ) over choosing an input sequence u that benefits the control objective J( ξ , u, y), which is visualized in Figure 3(d) and (e). These considerations highlight the need for isolation of the individual cost terms. In the following section, we present two modified DPC schemes in which, as revealed through an analysis of the trajectory-specific effect of regularization, this isolation of regularization cost terms naturally occurs.

4.2 Isolating trajectory-specific effects via projections or γ-DDPC

In [20], the orthogonal projection matrices

were introduced to yield a regularization h ( a ) = λ ‖ Π ⊥ a ‖ 2 2 for which the closed-loop behavior is consistent with the subspace predictor as λ → ∞. While it is noted that norm-based regularizers such as h ( a ) = λ ‖ a ‖ 2 2 “are not consistent and bias the optimal solution” [20], Section IV.C], the nature of this bias is not further explored. In [12], this nature is uncovered by showing that (under Assumption 1) the trajectory-specific effect of projection-based quadratic regularization h ( a ) = λ ‖ Π ⊥ a ‖ 2 2 is given by

Since ‖ a ‖ 2 2 = ‖ Π a ‖ 2 2 + ‖ Π ⊥ a ‖ 2 2 , we also immediately see that the remaining cost terms (10b) and (10c) are associated with Π. That is, (under Assumption 1) the trajectory-specific effect of projection-based quadratic regularization h ( a ) = λ ‖ Π a ‖ 2 2 is given by

This separation explains the observations in [20], Figure 2]. There, the performance of h ( a ) = λ ‖ a ‖ 2 2 deteriorates for λ → ∞ since the chosen input u is aggressively pushed towards u ̂ LS ( ξ ) , while neglecting the control objective J( ξ , u, y). In contrast, the projection-based regularization h ( a ) = λ ‖ Π ⊥ a ‖ 2 2 avoids this problem by fully dropping the input-related term (10b). However, the “closed-loop view” (discussed in Section 3) in [20], Figure 2] hides the fact that two functionally very different cost terms (10a) and (10b) are responsible for the observed performance. Furthermore [20], only considers the two options of either using both (10a) and (10b) with the same weight λ via h ( a ) = λ ‖ a ‖ 2 2 , or only using (10a) via h ( a ) = λ ‖ Π ⊥ a ‖ 2 2 . However, an individually weighted mix of the two terms may provide the best of both worlds, which is coincidentally explored in the regularization schemes proposed for γ -DDPC [18].

In γ -DDPC, the constraint (6b) and variable a are replaced via LQ decomposition as

where the diagonal blocks L _ii for i ∈ {1, 2, 3} are non-singular (under Assumption 1) and the matrices Q _i have orthonormal rows. Furthermore, Q ₄ and γ ₄ are typically omitted, since they do not affect the generated trajectory. The idea behind γ -DDPC can be summarized as re-parameterizing the OCP with a lower dimensional variable γ , and decoupling the matching of the initial condition ξ , since γ 1 = L 11 − 1 ξ is uniquely determined [18]. The proposed regularization strategies for γ -DDPC are based on a mix of quadratic regularization h ̃ ( γ ) = λ 2 ‖ γ 2 ‖ 2 2 + λ 3 ‖ γ 3 ‖ 2 2 , or constraining γ ₃ = 0 and using only h ̃ ( γ ) = λ 2 ‖ γ 2 ‖ 2 2 [23], [26]. The connection between this regularization of γ -variables and projection-based regularization of a was revealed in [13] by analyzing their trajectory-specific effect.

We note that the difference in notation w.r.t. to Definition 1 only comes from the fact that ( γ ₁, γ ₂, γ ₃) are uniquely determined by ( ξ , u, y), and therefore γ i * ( ξ , u , y ) = γ i . From Corollary 1, we can see that both the projection-based regularization and the regularization of γ _i isolate specific terms of the trajectory-specific effect of standard quadratic regularization. Their relation can be summarized as follows

As a guideline for practitioners, we advocate to use a mixed regularization h ( a ) = λ 2 ‖ Π a ‖ 2 2 + λ 3 ‖ Π ⊥ a ‖ 2 2 (or, equivalently in the γ -DDPC framework h ̃ ( γ ) = λ 2 ‖ γ 2 ‖ 2 2 + λ 3 ‖ γ 3 ‖ 2 2 ). This approach generalizes both the quadratic regularization h ( a ) = λ ‖ a ‖ 2 2 , and the purely projection based h ( a ) = λ ‖ Π ⊥ a ‖ 2 2 by choosing λ = λ₂ = λ₃ or λ = λ₃, λ₂ = 0, respectively. To avoid output predictions blindly following the control objective, λ₃ should be chosen very large, and (except for potential numerical solver instabilities) there are no adverse effects to be expected even for λ₃ → ∞. On the contrary, λ₂ should be chosen very carefully to avoid the inputs blindly following the least-squares controller u ̂ LS ( ξ ) . While λ₂ = 0 (recovering the projecion-based regularization h ( a ) = λ 3 ‖ Π ⊥ a ‖ 2 2 ) can be considered a safe choice in this regard, it also foregoes the potential performance improvements associated with proper tuning of this term. Here, we refer to [26] for more specific tuning guidelines. Finally, as also observed in [19], adding or removing trajectory columns from D after the tuning requires a re-scaling of all regularization weights

to retain a quantitatively similar effect, since Q reg , R reg are not normalized w.r.t. the number of trajectory columns.

An implicit predictor for DPC with quadratic regularization

In the following, we assume that the control objective is a quadratic output-tracking formulation

with reference y_ref, positive semidefinite weighing matrix Q , and arbitrary input control objective J _u ( ξ , u). We further assume that no additional output constraints are present, i.e., Y = R p N f . For some general observations on the effect of output constraints, see [12], Section III.C], and specific results regarding (terminal) equality constraints and general feasibility are given in Sections 6.4 and 6.5, respectively. In [12], Thm 3, 4], results for an implicit predictor with reference y_ref = 0 are given, which were later extended to possibly nonzero reference in [14], Thm. 2]. Under Assumption 1, solving (16) comes down to solving the unconstrained optimization problem

Furthermore, since the trajectory-specific cost of h ( a ) = λ ‖ a ‖ 2 2 and h ( a ) = λ ‖ Π ⊥ a ‖ 2 2 only differ in terms independent of y (see (10) and (15)), it turns out that the predictive behavior of both DPC schemes can be characterized by the same implicit predictor

Note that the involved inverse exists because Q reg is positive definite (under Assumption 1) and Q is positive semidefinite. Structurally, (19) can be seen as a (matrix-)weighed sum, shifting output predictions towards y_ref and y ̂ LS ( ξ , u), depending on the weights Q and λ Q reg , respectively. This structure is visualized in Figure 3 for the low-dimensional state-space example discussed in Section 2.3. The predictions blindly follow y_ref (independent of ξ , u) for λ = 0 and tend towards y ̂ LS ( ξ , u ) for λ → ∞. Furthermore, y ̂ DPC ( ξ , u ) is generally non-causal and time-variant (which cannot be observed in Figure 3 due to the limited prediction horizon N _f = 1). For finite λ and y_ref ≠ 0 , it is also an affine predictor (otherwise linear), regardless of the true system class generating the data. This last point highlights the interpretation of implicit predictors being “the predictive behavior implicitly attributed to the data-generating system by the DPC scheme” [12], which may not necessarily match the true system if the scheme is poorly chosen.

Figure 3:

Implicit predictor, optimal parametric solutions, and least-squares mappings for the DPC problem discussed in Section 2.3 (a–d). The optimal parametric DPC solutions ( x ₀, u*( x ₀), x*( x ₀)) for the different regularizations h ( a ) = λ ‖ a ‖ 2 2 (orange) and h ( a ) = λ ‖ Π ⊥ a ‖ 2 2 (green) naturally evolve on the implicit predictor x ̂ DPC ( x 0 , u ) (gray). (e) Visualization of the least-square solutions x ̂ LS ( x 0 , u ) (gray) and u ̂ LS ( x 0 ) . The latter is shown via the tuple x 0 , u ̂ LS ( x 0 ) , x ̂ LS ( x 0 , u ̂ LS ( x 0 ) ) (orange).

Figure 4:

Implicit predictor x ̂ DPC ( x 0 , u ) for the DPC example discussed in Section 2.3 with prediction horizon N _f = 2 and reference. (a) No terminal constraints. (b) Terminal constraints. x ref = ( 0.5 0.5 ) ⊤ . To deal with the higher dimensionality, the first (green) and second (orange) prediction step of x ̂ DPC ( x 0 , u ) are visualized individually. Furthermore, we constrain the second input to u₂ = 0 and only visualize the first input u₁.

6.1 DPC for affine systems

While standard (linear) DPC can yield exact predictions for deterministic LTI systems (see the discussion in Section 2.1), exact extensions to particular classes of nonlinear systems have been proposed, e.g., in [6], [7], [22]. Among these, we want to briefly discuss the case of affine time-invariant (ATI) systems

proposed in [22], which is also used for other nonlinear systems with continuously updated trajectory data in order to approximate a local (affine) linearization of the nonlinear system for predictions as in [10], [22]. Similarly to how trajectories of LTI systems (2) can be generated by linear combinations, trajectories of ATI systems (20) can be generated by affine combinations of trajectory data. That is, in addition to (1), generated trajectories must also satisfy

Intuitively, this condition can be explained by noting that the effect of e , r is present exactly once in each data trajectory and, accordingly, should be present exactly once in the generated trajectory. Assuming exact data generated by an ATI system, the affine hull a ff ( D ) of data columns in D is equivalent to the set of all possible system trajectories iff [29]

However, in the presence of noise and (other) nonlinearities, the same discussions as in Section 2.1 apply, i.e., the unique and exact predictions are no longer possible. Accordingly, our analysis is not confined to affine DPC applied to data from ATI systems (20), but rather extends to affine DPC with data generated by any system, including the nonlinear tracking case in [10]. To understand the features of affine DPC in the presence of such realistic data, we extend our results from the linear DPC case. As discussed in [19], Section II] and, in particular, [19], Rem. 4], many analysis results for linear DPC also apply to such nonlinear systems, which are linear in known (nonlinear) transformations of the state, input, and output. This also applies to the affine DPC scheme at hand, where we can simply consider ξ ̌ ≔ 1 ξ ⊤ ⊤ to be a nonlinear transformation of ξ , and W ̌ ≔ 1 ℓ W ⊤ ⊤ contains the corresponding data. Similarly, we define

for the transformed state-input data. Instead of the linear least-squares estimates y ̂ LS ( ξ , u ) , u ̂ LS ( ξ ) , consider their affine counterparts

with

and the corresponding residual matrices

The analysis of affine DPC with regularizations in terms of their trajectory-specific effect and implicit predictors then follows accordingly. In the following, we present the case of (projection-based) quadratic regularization.

The interpretation of these cost terms also directly follows from the discussion below (10). Importantly, note that instead of the least-squares estimates for a linear predictor/controller y ̂ LS ( ξ , u ) , u ̂ LS ( ξ ) the regularization now favors the least-squares estimates for the affine predictor/controller y ̂ ALS ( ξ , u ) , u ̂ ALS ( ξ ) and shifts the predicted trajectory towards them. Similarly to (15), the effect of the cost terms (23b) and (23c) can be eliminated by considering a projection-based quadratic regularization h ( a ) = λ ‖ Π ̌ ⊥ a ‖ 2 2 with

Finally, the analysis of predictive behavior via implicit predictors follows accordingly.

Again, the interpretation of this predictor follows from the linear DPC case discussed below (19).

6.2 Regularization with offset

Some DPC schemes, in which the tracking of a non-zero equilibrium is desired, modify the regularization by including an offset, i.e., h ̃ ( a ) = h ( a − a ̄ ) [10], [30]. Although this adjustment is based on the assumption that the original regularization shifts predicted trajectories toward zero (see our discussion in Section 3), it introduces some interesting features, which can be analyzed via its trajectory-specific effect. Although [30] considers various norm choices, we focus on the quadratic regularization h ( a ) = λ ‖ a − a ̄ ‖ 2 2 used in both [10], [30]. Additionally, while the results in [10] are framed for an affine DPC setting, our analysis will focus on linear DPC as in [30]. Note that these results can be readily extended to the affine case, as discussed in Section 6.1. Analogously to how a generates the trajectory w = D a , we denote the trajectory generated by a ̄ as w ̄ ≔ D a ̄ , and similarly with z ̄ , u ̄ , y ̄ , ξ ̄ . Furthermore, we denote the difference between the two as

and similarly with Δ z , Δu, Δy, Δ ξ .

Hence, an offset by a ̄ in the regularization simply translates to the same trajectory-specific effect with an offset w ̄ for the considered trajectory. Note that (25) can also be equivalently expressed as

While the third cost term is always irrelevant, we can see that the first two terms also may have no additional effect (compared to the usual quadratic regularization) if the trajectory w ̄ adheres to

respectively. Similarly to (15), the effect of the cost terms (25b) and (25c) can be eliminated by considering a projection-based quadratic regularization. On that note, we briefly remark that both

yield the same effect, since they only differ in a constant term unrelated to a . Furthermore, using the first term of the alternative cost expression (26), we can also state an implicit predictor as follows.