Robust and Kernelized Data-Enabled Predictive Control for Nonlinear Systems
Abstract
This paper presents a robust and kernelized data-enabled predictive control (RoKDeePC) algorithm to perform model-free optimal control for nonlinear systems using only input and output data. The algorithm combines robust predictive control and a non-parametric representation of nonlinear systems enabled by regularized kernel methods. The latter is based on implicitly learning the nonlinear behavior of the system via the representer theorem. Instead of seeking a model and then performing control design, our method goes directly from data to control. This allows us to robustify the control inputs against the uncertainties in data by considering a min-max optimization problem to calculate the optimal control sequence. We show that by incorporating a proper uncertainty set, this min-max problem can be reformulated as a nonconvex but structured minimization problem. By exploiting its structure, we present a projected gradient descent algorithm to effectively solve this problem. Finally, we test the RoKDeePC on two nonlinear example systems – one academic case study and a grid-forming converter feeding a nonlinear load – and compare it with some existing nonlinear data-driven predictive control methods.
Index Terms:
Data-driven control, kernel methods, nonlinear control, predictive control, robust optimization.I Introduction
Data-driven control has attracted extensive attention in recent years, as it enables optimal control in scenarios, where data is readily available, but the system models are too complex to obtain or maintain; this is often the case, for example, in large-scale energy systems or robotics [1, 2, 3].
One standard approach to data-driven control is indirect. Here one first uses data to identify a system model and then performs control design based on the identified model. This paradigm has a long history in control theory, leading to the developments of system identification (SysID) [4], model predictive control (MPC) [5], etc. Prediction error, maximum likelihood, and subspace methods are popular SysID approaches, and have been successfully combined with model-based control design in many industrial applications. Broadly recognized challenges of the indirect approach are incompatible uncertainty estimates and customizing the SysID (e.g., in terms of objective and model order selection) for the ultimate control objective. We refer to [6, 7, 8, 9] for different approaches. Moreover, recent advances in systems theory, statistics, machine learning, and deep learning have introduced many promising techniques to learn the behavior of general dynamical systems from data, e.g., by using Koopman operators, Gaussian processes, kernel methods, and neural networks [10, 11, 12, 13, 14, 15]. The learned models can also be combined with model-based control design (e.g., MPC) to perform optimal control [16, 17, 18, 19, 20, 21]. However, these methods still belong to indirect data-driven control, and one may not be able to provide deterministic guarantees on the control performance, especially when complicated structures, e.g., neural networks, are used to learn the system’s behaviors. We note that compared to using neural networks, function estimation using regularized kernel methods has a tractable formulation and solution, thanks to the well-posedness of the function classes in reproducing kernel Hilbert spaces (RKHSs) [12, 22].
An alternative formulation is direct data-driven control that aims to go directly from data to control, bypassing the identification step. In recent years, a result formulated by Willems and his co-authors in the context of behavioral system theory, known as the Fundamental Lemma [23], has been widely used in formulating direct data-driven control methods such as [24, 25, 26, 1]. The Fundamental Lemma shows that the subspace of input/output trajectories of a linear time-invariant (LTI) system can be obtained from the column span of a data Hankel matrix, which acts as a data-centric representation of the system dynamics. Unlike indirect data-driven control methods which usually involve a bi-level optimization problem (an inner problem for the SysID step and an outer problem for the control design step), these direct data-driven control methods formulate one single data-to-control optimization problem. In this paradigm, one can conveniently provide robust control against uncertainties in data via regularization and relaxation [8]. For instance, in the data-enabled predictive control (DeePC) algorithm [24], a proper regularization leads to end-to-end distributional robustness against uncertainties in data [27, 28]. Moreover, performance guarantees on the realized input/output cost can be provided given that the inherent uncertainty set is large enough [29]. Motivated in part by the desire to extend this direct data-driven approach beyond LTI systems, Fundamental Lemma has been extended for certain classes of nonlinear systems, e.g., Hammerstein-Wiener systems [30], second order discrete Volterra systems [31], flat nonlinear systems [32], and polynomial time-invariant systems [33], which enables the control design of these classes of nonlinear systems. However, to the best of our knowledge, there has not been a version of Fundamental Lemma for general nonlinear systems.
In this paper, we develop a direct data-driven method for nonlinear systems, referred to as robust and kernelized DeePC (RoKDeePC), which combines kernel methods with a robustified DeePC framework. We use regularized kernel methods [22] to implicitly learn the future behavior of nonlinear systems, then directly formulate a data-to-control min-max optimization problem to obtain robust and optimal control sequences. We prove that the realized cost of the system is upper bounded by the optimization cost of this optimization problem, leading to deterministic performance guarantees. Furthermore, we demonstrate how this min-max problem can be reformulated as a nonconvex yet structured minimization problem when appropriate uncertainty sets are considered. To enable a real-time implementation, we develop a projected gradient descent algorithm exploiting the problem structure to effectively solve the RoKDeePC problem. Unlike learning-based control using neural networks, our method does not require a time-consuming off-line learning process and has real-time plug-and-play ability to handle nonlinear systems. We test RoKDeePC on two nonlinear example systems – one academic case study and a grid-forming converter feeding a nonlinear load – and compare its performance with the Koopman-based MPC method in [17] and MPC that uses a kernel-based predictor as the system model. These methods are all based on linear predictors in the lifted high-dimensional space and thus fair comparisons can possibly be made.
The rest of this paper is organized as follows: in Section II we review the linear predictor and then introduce the kernel-based nonlinear predictor using input/output data. Section III proposes the RoKDeePC algorithm and demonstrates its performance guarantees. Section IV develops gradient descent algorithms to efficiently solve the RoKDeePC. Section V tests the RoKDeePC algorithm on two example nonlinear systems. We conclude the paper in Section VI.
Notation
Let denote the set of positive integers, the set of real -by- symmetric (positive definite) matrices. We denote the index set with cardinality as , that is, . We use () to denote the (induced) 2-norm of the vector (the matrix ); we use to denote the 1-norm of the vector ; for a vector , we use to denote ; for a matrix , we use to denote ; we use to denote the Frobenius norm of the matrix . We use to denote the element of and to denote the vector ; we use to denote the column of , to denote the row of , and to denote -row -column element of . We use to denote a vector of ones of length , to denote a -by- matrix of ones, and to denote an -by- identity matrix (abbreviated as when the dimensions can be inferred from the context). We use to denote the matrix .
II Kernel-Based Nonlinear Predictors
II-A Explicit Predictors for LTI Behavior
Consider a discrete-time linear time-invariant (LTI) system
(1) |
where , , , , is the state of the system at , is the input vector, and is the output vector. Recall the extended observability matrix
The lag of the system (1) is defined by the smallest integer such that the observability matrix has rank , i.e., the state can be reconstructed from measurements. Consider with and length- input and output trajectories of (1): and . For the inputs , define the Hankel matrix of depth as
(2) |
Accordingly, for the outputs define the Hankel matrix . Consider the stacked matrix . By [34, Corollary 19], the restricted behavior (i.e., the input/output trajectory of length ) equals the image of if and only if rank.
This behavioral result can be leveraged for data-driven prediction and estimation as follows. Consider , as well as an input/output time series so that rank. Here the superscript “d” denotes data collected offline, and the rank condition is met by choosing to be persistently exciting of sufficiently high order and by assuming that the system is controllable and observable. The Hankel matrix can be partitioned into
where , , , , and . In the sequel, the data in the partition with subscript P (for “past”) will be used to implicitly estimate the initial condition of the system, whereas the data with subscript F will be used to predict the “future” trajectories. In this case, is the length of an initial trajectory measured in the immediate past during on-line operation, and is the length of a predicted trajectory starting from the initial trajectory. The Fundamental Lemma shows that the image of spans all length- trajectories [23], that is, is a trajectory of (1) if and only if there exists a vector so that
(3) |
The initial trajectory can be thought of as setting the initial condition for the future (to be predicted) trajectory . In particular, if , for every given future input trajectory , the future output trajectory is uniquely determined through (3) [35].
In this case, one can consider an explicit linear predictor as
(4) |
where is a linear mapping and can be calculated from data as
(5) |
where denotes the pseudoinverse operator. Note that (5) is the solution of the following linear regression problem
(6) |
where can be considered as one sample trajectory of .
II-B Explicit Kernel-Based Nonlinear Predictors
We now consider nonlinear systems, where the future outputs are nonlinear functions of the initial trajectory and the future inputs
(7) |
In what follows, we focus on how to find the nonlinear functions that best fit the observed nonlinear trajectories contained in .
To reconcile flexibility of the functions class of and well-poseness of the function estimation problem, we assume that the functions belong to an RKHS [38].
Definition II.1.
An RKHS over a non-empty set is a Hilbert space of functions such that for each , the evaluation functional is bounded.
An RKHS is associated with a positive semidefinite kernel, called reproducing kernel, which encodes the properties of the functions in this space.
Definition II.2.
A symmetric function is called positive semidefinite kernel if, for any ,
The kernel section of centered at is .
According to the Moore–Aronszajn theorem [39], an RKHS is fully characterized by its reproducing kernel. To be specific, if a function belongs to an RKHS with kernel , then it can be expressed as
(8) |
for some , for all , possibly with infinite . Then, the inner product in the RKHS reduces to
where and , and the induced norm of to .
We consider now the function estimation problem for the nonlinear system in (7). Since every column in the data matrix is a trajectory of the system and thus an input/output sample for (7), for each , we search for the best (in a regularized least-square sense) in an RKHS with kernel
(9) |
where . Note that to simplify the forthcoming developments, the formulation in (9) does not consider temporal correlation between different indices and causality (i.e., is not affected by ). One can further consider causality by restricting the arguments of to be , which will lead to more complicated formulations than those in this paper, as different is defined over sets with different dimensions.
The cost function in (9) includes the prediction errors and a regularization term, which is used to avoid over-fitting and penalize undesired behaviors. We note that since the functions in an RKHS can be parameterized as in (8), the representer theorem [22] ensures that the minimization problem in (9) can be solved in closed form. In our context, this leads to the following.
Lemma II.3.
Eq. (10) serves as a data-driven nonlinear explicit predictor which can be used to predict future behaviors of the system (7) and embedded in predictive control algorithms. We note that the choice of kernel encodes the properties of the function space. Indeed, (10) shows that the synthesized functions are linear combinations of the kernel sections centered at the sampling points. Hence, one can choose the kernel based on a priori knowledge on the system. For example, if one knows that () are polynomials in inputs and initial data, then a polynomial kernel can be used. When the functions are very high-order polynomials or other complex forms, a Gaussian kernel can be used to span the set of continuous functions.
III DeePC with Kernel Methods
III-A Review of the DeePC algorithm
The DeePC algorithm proposed in [24] directly uses input/output data collected from the unknown LTI system to predict the future behaviour, and performs optimal and safe control without identifying a parametric system representation. More specifically, DeePC solves the following optimization problem to obtain the optimal future control inputs
(11) |
where is the cost function of control actions, for instance, or . The sets of input/output constraints are defined as and . The positive definite matrix and positive semi-definite matrix are the cost matrices. The vector is a prescribed reference trajectory for the future outputs. DeePC involves solving the convex quadratic problem (11) in a receding horizon fashion, that is, after calculating the optimal control sequence , we apply to the system for time steps, then, reinitialize the problem (11) by updating to the most recent input and output measurements, and setting to , to calculate the new optimal control for the next time steps. As in MPC, the control horizon is a design parameter.
The standard DeePC algorithm in (11) assumes perfect output data generated from the unknown system (1). However, in practice, perfect data is not accessible to the controller due to measurement noise. We note that regularized DeePC algorithms can be used to provide robustness against disturbances on the data, which often amount to one-norm or two-norm regularization of in the cost function of (11) [27, 29]. In this paper, we assume no process noise, i.e., the input data is perfectly known.
Perfect Data and Noisy Data. Throughout the paper, we use , , and to denote the perfect (noiseless) output data generated from the system, which accurately captures the system dynamics; we use , , and to denote the corresponding noisy output data. We use to denote the Gram matrix calculated from perfect output data, and to denote the Gram matrix calculated from noisy output data.
III-B RoKDeePC algorithm
The DeePC algorithm above implicitly incorporates a linear predictor. In this section, we develop a RoKDeePC algorithm using an implicit and robustified version of the kernel-based nonlinear predictor (10) to perform data-driven optimal control for nonlinear systems. We first consider the following certainty-equivalence MPC problem using the kernel predictor in (10) (referred to as kernel-based MPC in this paper)
(12) |
In what follows, we assume that when perfect data is available, the predictor (10) can faithfully predict the future behavior of the nonlinear dynamics (7).
Assumption 1.
When perfect data is available in the kernel-based nonlinear predictor (10), the prediction error is bounded, i.e.,
(13) |
where is the predicted output trajectory based on (10), is a control sequence applied to the system, is the realized output trajectory of the system in response to , and is the error bound.
Assumption 1 implicitly requires that sufficiently rich data is used to capture the system dynamics [40]. Moreover, one can possibly use more data to obtain a more accurate prediction and thus reduce . Since perfect data is available, this assumption can generally be satisfied if the unknown dynamics are contained in the RKHS of the chosen kernel , e.g., is polynomial and a sufficiently high-order polynomial kernel is used. Otherwise, this assumption can be met by using kernels with the universal approximation property (e.g., a Gaussian kernel) and by implicitly assuming that is continuous and that is contained in a compact set [41].
When noisy data is used, the prediction error may grow unbounded, and the obtained control sequence from (12) may not result in a robust closed loop. As a first step toward addressing this problem, we rewrite (10) implicitly as
(14) |
where is uniquely determined by since has full rank, and predicts the future behavior of the system. We note that, given , there exists a unique that satisfies (14). Unlike linear systems, however, given an arbitrary , there may not exist a vector that satisfies (14); for linear systems such a vector exists, thanks to the Fundamental Lemma. Substituting (14) into (12) leads to
(15) |
The equality constraint (14) assumes perfect data to obtain a satisfactory prediction. When perfect data is not available, one may need to relax this equality constraint and ensure feasibility, for instance, by penalizing its violation in the cost function (i.e., soft constraints)
(16) |
where is a penalty parameter. We note that with a sufficiently large , the solution to (16) is equal to the solution to (15) [42]. However, a robust solution may still not be obtained when noisy output data is used. Hence, we further robustify (16) and propose the following RoKDeePC formulation that leverages noisy output data for nonlinear, robust, and optimal control
(17) |
where and are disturbances added to the data matrices, which respectively reside in prescribed disturbance sets and , and robustifies the output constraints.
The above min-max formulation considers the worst-case scenario of the uncertainties added to the data matrices, thereby leading to a robust solution. We will later show that performance guarantees can also be provided by employing this formulation. Note that min-max optimization problems could be in general difficult to solve, however, (17) admits a tractable reformulation as a minimization problem if appropriate disturbance sets and are considered. The following result shows that the inner max problem of (17) can be resolved analytically for appropriate disturbance sets and , resulting in regularized problem formulations.
Lemma III.1.
Consider , and . Given a vector and a vector , it holds that
(18) |
(19) |
Proof.
The proof is adapted from [43, Theorem 3.1]. Consider a fixed and a fixed in (18). It follows from triangle inequality that
Choose such that
Since is a rank-one matrix, we have , and
where the first equality follows from triangle inequality and the fact that is a nonnegative scalar multiple of . Since the lower bound coincides with the upper bound on the maximization problem (18) for any fixed and , the equality in (18) holds. By following a similar process, one can prove the equality in (19). ∎
Corollary III.2 (A tractable formulation for RoKDeePC).
The above result shows that one can get robust control inputs by simply solving a regularized minimization problem. The proposed RoKDeePC algorithm involves solving (20) in a receding horizon manner, similar to the DeePC algorithm. The initial trajectory should be updated every time before solving (20). We remark that and are unstructured uncertainties sets which may lead to conservativeness. For instance, is a symmetric matrix, but does not restrict to be symmetric; may be a Hankel data matrix, but a Hankel structure cannot be imposed on by considering . As will be discussed in the simulation section, though the considered unstructured sets are conservative, they can generally lead to satisfactory performance.
III-C Performance guarantees of RoKDeePC
In what follows, we show that the proposed RoKDeePC algorithm provides deterministic performance guarantees when applying a feasible optimal control sequence to the system, thanks to the min-max formulation in (17). We begin by making the following assumption.
Assumption 2.
For any (contained in a compact set), there exists such that and ; there exists such that .
Assumption 2 indicates that sufficiently large uncertainty sets and are used to cover the deviations of the data matrices induced by noise. This is necessary for (17) to provide a robust solution against the realized uncertainty in practice. Since is contained in a compact set and is continuous (by choosing a continuous kernel), Assumption 2 can be easily satisfied when the output noise is bounded and sufficiently large disturbance sets and are considered, which corresponds to choosing sufficiently large regularization parameters and in (20). We define the realized cost of the system as (with defined in Assumption 1). The following result shows that the realized cost is upper-bounded by the optimization cost in (21) plus the prediction error in (13).
Theorem III.3.
Proof.
According to the definition of , if Assumption 2 holds, we have
(23) |
where minimizes (17) in which . Given a vector , it follows from (14) that there exists a such that
(24) |
where is the predicted output trajectory from the kernel-based nonlinear predictor. By defining , we have
(25) |
which implies that
(26) |
By substituting (25) and (26) into (23) we obtain
(27) |
where the second and the forth inequalities hold thanks to the reverse triangle inequality, the third inequality is satisfied if we take large enough (by assumption) to ensure
Hence, , and the fifth inequality follows from Assumption 1 and the definition of . This completes the proof. ∎
The above result links the optimization cost to the realized cost, namely, the realized cost is always bounded by the optimization cost plus the prediction error even when noisy data is used. This also justifies the design of the cost function of the RoKDeePC in (20), and shows that one should try to find not only a feasible solution but ideally a minimizer (or better the global minimizer) to reduce the optimization cost and possibly the realized cost.
IV Gradient Descent Methods to Solve RoKDeePC
The RoKDeePC problem in (20) is nonconvex in unless the kernel function is linear in . One may use a generic nonlinear optimization toolbox (e.g., IPOPT) to solve this problem. However, solving a nonconvex problem is in general time-consuming, especially considering that the decision variable in (20) can be high-dimensional. We propose customized gradient descent methods exploiting the structure of (20) to solve the RoKDeePC problem more efficiently. Although a global minimum may still not be reached, the numerical study in Section V suggests that these methods generally lead to satisfactory performance with a reasonable running time.
IV-A Gradient descent algorithm for RoKDeePC
Notice that the optimization problem in (20) is nonconvex in but convex in . Moreover, the dimension of is generally much higher than that of . Due to this favorable structure, we apply a gradient descent method to update , while at each iteration, we solve a convex optimization problem to update . Note that the input/output constraints can be handled using projected gradient algorithms. We summarize the algorithm for solving RoKDeePC below.
Input: cost function:
(28) |
initial vectors , ; constraint sets , ; step size ; maximal iteration number ; convergence threshold: .
Initialization: .
Iterate until or :
-
1.
-
2.
-
3.
Output: (sub)optimal control sequence .
In Step 2, the projection of into the set admits a closed-form solution if upper and lower bounds for the elements in are considered. Step 3 requires solving a second-order cone program (SOCP), which could make the algorithm time-consuming as one may need to solve the SOCP times. To reduce the computational burden, in what follows, we consider an equivalent cost function which is quadratic in and thus admits a closed-form solution for Step 3 when the output constraints are inactive.
IV-B Quadratic reformulation
To enable a fast calculation of Step 3, we consider the following cost function that is quadratic in
(29) |
The following result shows the connection between the cost functions in (28) and (29).
Proposition IV.1.
For any , if is a minimizer of
(30) |
then also minimizes
(31) |
with
(32) |
(33) |
Proof.
See Appendix A. ∎
The above result shows that using the cost function in the RoKDeePC is equivalent to setting different parameters (, , and ) to the problem in (20) every time solving it during the receding horizon implementation, as the equivalent values of , , and are related to . This may lead to conservativeness as we need to choose a robust pair of . However, we generally observe excellent performance when using as the cost function for RoKDeePC. During the online implementation, after solving (30), one can compute , , and (by fixing a desired value of or and then compute the other) according to (32) and (33). If the obtained , , and are not sufficiently large to satisfy Assumption 2 and Theorem III.3, then and should be increased to implicitly increase , , and .
Moreover, when the output constraints are inactive, (30) admits a closed-form solution as
(34) |
where , and . Notice that (34) requires simply a linear mapping to obtain , which can greatly accelerate the gradient descent algorithm when the output constraints are inactive. In summary, we consider the following modified RoKDeePC optimization problem
(35) |
which is related to the original formulation in (20) via Proposition IV.1, and we implement the following modifications to Algorithm 1 to solve (35):
-
•
Change the cost function from to ;
-
•
Change Step 3 in the iteration to:
if , .
Step 3 becomes easier to solve if the output constraint set is a nonempty box constraint (i.e., upper and lower bounds for the outputs are considered), as the constraint set can be formulated as for some and [44, 29]. The computation cost can be further reduced by assuming that is bounded and (conservatively) restrict to be a polyhedron for some . In our applications, we do not find these assumptions to be limiting.
We remark that one can also assume a kernel function that is nonlinear in but linear in such that (29) becomes a convex quadratic program, which can be more efficiently solved using standard solvers. For instance, one may consider
(36) |
where . Note that (36) is a combination of two positive semidefinite kernels and thus also satisfies the positive semidefinite property. If the future outputs of the system is linear in , such a kernel can be used to accurately capture the system’s dynamics. Otherwise, the obtained kernel-based predictor is the linear-input model that best fits the observed data, which may not lead to an accurate prediction.
V Simulation Results
In this section, we test the RoKDeePC algorithm on two example nonlinear systems to illustrate its effectiveness. We will compare the performance when different kernel functions are used, and we will compare the RoKDeePC algorithm with the traditional DeePC algorithm, the kernel-based MPC in (12), and the Koopman-based MPC in [17].
V-A Example 1: a nonlinear second-order model
We start with an academic case study to compare different predictors. Consider the following discrete-time nonlinear model of a polynomial single-input-single-output system
(37) |
To excite the system and collect input and output data, we inject white noise (mean: 0; variance: 0.01) into the input channel and collect an input/output trajectory of length , which will be used in different predictors below. The observed data may also be subject to additive white measurement noise.
We first test the prediction quality of different predictors, including: i) the linear predictor in (4), ii) the kernel predictor in (10), and iii) the Koopman-based predictor in [17]. The parameters in (4) and (10) are: , , and . We consider three popular kernels for the kernel predictor: 1) Polynomial kernel ; 2) Gaussian kernel ; 3) Exponential kernel . In this section, the parameters in these kernel functions are chosen to ensure that the kernel predictors have satisfactory performance when noiseless data is available. For instance, the Gaussian kernel has only one parameter to tune, and thus an appropriate value can be easily found with a line search.
For the Koopman-based predictor, one can consider as the state variables of the above system. Similar to [17], we consider the lifting functions in the form of
(38) |
where is in general nonlinear but is linear. The lifting functions are chosen to be the linear and quadratic terms of the state variables and 10 thin plate spline radial basis functions (defined by whose center is ) with centers selected randomly from . Under this setting, the input/output data can be used in the extended dynamic mode decomposition (EDMD) [10] to solve for the parametric system representation (Koopman-based predictor). The lifting functions in (38) leads to a quadratic program for the MPC formulation, which can be solved efficiently in practice.
Prediction Error | Without measurement noise | With white noise (variance: ) |
Linear predictor | 0.2378 | 0.2384 |
Koopman predictor (linear in ) | 0.1927 | 0.2175 |
Kernel predictor (Polynomial kernel) | 0.0008 | 0.0486 |
Kernel predictor (Gaussian kernel) | 0.0051 | 0.0463 |
Kernel predictor (Exponential kernel) | 0.0030 | 0.0456 |
noise level | realized costs | RoKDeePC Gaussian kernel | Kernel-based MPC Gaussian kernel | RoKDeePC Exponential kernel | Kernel-based MPC Exponential kernel | RoKDeePC Polynomial kernel | Kernel-based MPC Polynomial kernel |
variance: | mean | 8.92 | 8.51 | 6.07 | 6.16 | 7.31 | 7.36 |
standard deviation | 4.17 | 3.72 | 1.97 | 2.01 | 3.51 | 3.25 | |
variance: | mean | 25.06 | 28.19 | 25.11 | 29.37 | 34.94 | 40.77 |
standard deviation | 17.99 | 19.24 | 15.31 | 17.66 | 24.97 | 31.43 | |
variance: | mean | 30.76 | 37.54 (+22%) | 31.89 | 40.71 (+28%) | 46.76 | 59.36 (+27%) |
standard deviation | 24.29 | 28.06 | 22.15 | 27.03 | 36.32 | 50.50 |
Table I compares the open-loop predictions of 50 time steps obtained by employing different predictors. For the linear predictor and the kernel predictor (whose prediction horizon is ), the time steps from to () are predicted based on the previous predictions. It can be seen that the linear predictor cannot make an accurate prediction due to the strong nonlinearity of the system. The Koopman-based predictor also cannot make an accurate prediction because (38) assumes independence between and , which is not the case in (37). Moreover, we observe similar prediction errors even with more data in the EDMD, e.g., (data not shown). The structure in (38) enables fast calculations of the associated MPC problem, but restricts the applicability in some classes of nonlinear systems. One can also assume a more general form of to improve the prediction accuracy, but it is out of the scope of this paper.
By comparison, all kernel predictors can predict the system’s behaviors with satisfactory accuracy. Moreover, when there is no measurement noise, the predictor with polynomial kernel almost perfectly predicts the system’s behaviors, as the future outputs are indeed polynomial functions of due to the bilinear structure in (37).

We next compare the time-domain performance when different control methods are applied: 1) DeePC with quadratic regularization [28]; 2) Koopman-based MPC [17]; 3) nominal (i.e., non-robustified) Kernel-based MPC in (12) (the cost function is modified to be ); 4) RoKDeePC in (35). We consider the following input and output costs in all these methods: , where is the output reference over time. This cost function penalizes the control inputs, the tracking errors, and the rates of changes of control inputs. The other parameters of the RoKDeePC are: , . The control horizon is in all the methods. We assume no input/output constraints and focus on comparing the tracking performance. Hence, no projection is needed in the gradient descent algorithm.
Fig. 1 shows the time-domain responses of the system when different control methods are applied, where the output reference steps from to at , and then to at . The sampling time of the discrete-time system is 2ms. It can be seen that when DeePC or Koopman-based MPC is applied, the output cannot accurately track the reference because the linear predictor and the Koopman predictor (linear in ) cannot make an accurate prediction (Table I). By comparison, the RoKDeePC and the kernel-based MPC both achieve satisfactory performance even under measurement noise.
We remark that with the gradient descent algorithm in Section IV, the RoKDeePC has a reasonable running time. On an Intel Core i7-9750H CPU with 16GB RAM, it requires approximately 0.2s to solve (35) (with any of the considered kernels) every time in the above simulations. The kernel-based MPC requires similar time to solve (12) using a similar gradient descent algorithm. By comparison, the Koopman-based MPC and DeePC only require to solve quadratic programs, which take about 2ms and 10ms respectively using OSQP [45] as the solver. However, their performance is vastly inferior.
To compare the robustness of the RoKDeePC and the kernel-based MPC, the above simulations were repeated 100 times with different datasets to construct the Hankel matrices and different random seeds to generate the measurement noise. Table II displays the mean and standard deviation of realized closed-loop costs, i.e., where and are measured from the system in the simulations. It shows that with the increase of the measurement noise level, the RoKDeePC achieves superior performance than the (certainty-equivalence) kernel-based MPC. For instance, when the variance of the measurement noise is , the averaged closed-loop costs of the RoKDeePC are respectively (Gaussian kernel), (Exponential kernel), and (Polynomial kernel) lower than those of the kernel-based MPC. We attribute this performance gap to the robustness of the RoKDeePC, namely, the control sequence obtained from RoKDeePC is robust against the uncertainties in data; by comparison, the kernel-based MPC implicitly assumes certainty equivalence and may not lead to a robust solution. Moreover, we observe that the Gaussian kernel and Exponential kernel have better performance than the Polynomial kernel under noisy measurements, even though the future behaviors of the systems can be described by polynomial functions.
V-B Example 2: A grid-forming converter with nonlinear load
We now test the performance of the RoKDeePC in a grid-forming converter that is feeding a local nonlinear load, as shown in Fig. 2. The grid-forming converter employs virtual synchronous machine control for synchronization, which emulates the swing equation of synchronous generators to generate the frequency[46, 47, 48]. The converter is connected to a local nonlinear load, whose active and reactive power consumptions are (in per-unit values)
where is the terminal voltage magnitude of the load, and is the voltage deviation from its nominal value . The other system parameters are: , , , , , , and .

We apply the RoKDeePC algorithm to perform optimal voltage (or reactive power) control in the converter. As shown in Fig. 2, the input of the converter system is the internal voltage magnitude of the converter (i.e., ), and the output can either be the reactive power or the terminal voltage deviation of the load . Since the converter system is more complex than (37), we choose , , , and a control horizon of . The other settings of the RoKDeePC are the same as those in the previous subsection (cost function, excitation signals, etc.). We again assume no input/output constraints and focus on comparing the tracking performance. For the Koopman-based MPC, we consider as the state variables of the system, and the choices of the lifting functions are the same as those in the previous subsection.

Fig. 3 (a) shows the reactive power responses of the converter in grid-connected mode (i.e., the breaker in Fig. 2 is closed), and we choose the reactive power to be the output of the system. The sampling time of the system is 2ms. The reference for steps from to at , and then to at . It can be seen that the converter has satisfactory performance when the RoKDeePC is applied (with Gaussian, Exponential, or Polynomial kernel). However, when the Koopman-based MPC or DeePC is applied, there are severe tracking errors; with the former we also observe some oscillations that persist even with more data to solve the EDMD problem and a longer prediction horizon (e.g., and ). Similar to the results in Section V-A, this is because of the strong nonlinearity of the systems and the special structure of the lifting functions in (38) that are linear in . With the kernel-based MPC in (12) (with Gaussian, Exponential, or Polynomial kernel), the system became unstable and thus the time-domain responses are not displayed in Fig. 3.
Although the RoKDeePC (with Gaussian, Exponential, or Polynomial kernel) has superior performance, it requires solving a nonlinear program in real time. In this instance, it took approximately 0.6s to solve (35). By comparison, the Koopman-based MPC and DeePC require about 2ms and 80ms respectively using OSQP as the solver.
To enable a real-time implementation of the RoKDeePC, we further consider a kernel that is nonlinear in but linear in , as shown in (36), where . With this kernel, (35) becomes a quadratic program, which admits a closed-form solution when the input/output constraints are inactive. Fig. 3 (a) also plots the reactive power response of the converter when this kernel is used, which shows better performance than the Koopman-based MPC and DeePC (with , the tracking error is less then ). Moroever, it only requires about 10ms to calculate the closed-form solutions at each time step, and thus enables a possible real-time implementation with a control horizon .
We next test the control performance when the converter is in an islanded mode (i.e., the breaker in Fig. 2 is open). We choose the voltage deviation to be the output of the system, that is, the RoKDeePC aims at regulating the voltage deviation of the load. Since the voltage dynamics are in general fast, we choose the sampling time to be 0.5ms. Fig. 3 (b) shows the time-domain responses of the voltage deviation when its reference steps from to at , and then to at . We observe that the RoKDeePC (with Gaussian, Exponential, or Polynomial kernel) achieves the best performance, whereas the Koopman-based MPC and DeePC lead to substantial tracking error. The tracking error is about when the kernel in (36) is used in the RoKDeePC to enable faster calculations.
VI Conclusions
RoKDeePC, an algorithm to perform data-driven, robust, and optimal control for general nonlinear systems based on implicitly predicting the future behaviors of the system using the representer theorem was presented. RoKDeePC involves solving a data-to-control optimization problem in a receding-horizon manner, which does not require any time-consuming off-line learning step. Moreover, we demonstrate how to provide end-to-end robustification for the control sequence against measurement noise in the output data, leading to strong performance guarantees. To handle the nonconvexity of the optimization problem, we exploited its structure and developed projected gradient descent algorithms to enable fast calculations of the (sub)optimal solutions. The RoKDeePC was tested in two example nonlinear systems (including a high-fidelity converter model feeding a nonlinear load) and showed excellent performance with reasonable running time, compatible with real-time implementations in real-world applications. A comparison with related nonlinear data-driven MPC methods showed the favorable performance of our approach.
Acknowledgment
The authors would like to thank Liviu Aolaritei, Francesco Micheli, and Jianzhe Zhen for fruitful discussions.
Appendix A Proof of Proposition IV.1
Proof.
Since is a polyhedron, we can compactly rewrite as where for some , , and [44].
Consider the Lagrangian of (39)
where is the vector of the dual variables. Since is nonempty, there exists a solution to the Karush–Kuhn–Tucker (KKT) conditions of (39)
(41c) | |||
(41d) | |||
(41e) | |||
(41f) |
Consider the Lagrangian of (40)
where is the vector of the dual variables. By choosing , , and according to (32) and (33), it can be verified that , where
if , and otherwise, and as before, satisfy the KKT conditions of (40)
(42c) | |||
(42f) | |||
(42g) | |||
(42j) | |||
(42k) | |||
(42l) | |||
(42m) |
Next we prove the monotonic relationship between and . Let be the minimizer of (39) with (and the minimizer of (40) with ), and the minimizer of (39) with (and the minimizer of (40) with ). If , we directly obtain according to (32). If , according to the definitions of and , we have
which can be reorganized as
(43) |
Moreover, we have
which can be reorganized as
(44) |
By combining (43) and (44), we obtain
which indicates that as . Then, we have .
According to the definitions of and , we also have
leading to
It can then be deduced that . Hence, is increasing with the increase of . By following a similar process, one can also prove that () is increasing with the increase of . This completes the proof. ∎
References
- [1] I. Markovsky and F. Dörfler, “Behavioral systems theory in data-driven analysis, signal processing, and control,” Annual Reviews in Control, vol. 52, pp. 42–64, 2021.
- [2] L. Huang, J. Coulson, J. Lygeros, and F. Dörfler, “Decentralized data-enabled predictive control for power system oscillation damping,” IEEE Transactions on Control Systems Technology, 2021.
- [3] E. Elokda, J. Coulson, P. N. Beuchat, J. Lygeros, and F. Dörfler, “Data-enabled predictive control for quadcopters,” International Journal of Robust and Nonlinear Control, vol. 31, no. 18, pp. 8916–8936, 2021.
- [4] L. Ljung, “System identification,” Wiley Encyclopedia of Electrical and Electronics Engineering, pp. 1–19, 1999.
- [5] M. Morari and J. H. Lee, “Model predictive control: past, present and future,” Computers & Chemical Engineering, vol. 23, no. 4-5, pp. 667–682, 1999.
- [6] H. Hjalmarsson, M. Gevers, and F. De Bruyne, “For model-based control design, closed-loop identification gives better performance,” Automatica, vol. 32, no. 12, pp. 1659–1673, 1996.
- [7] S. Formentin and A. Chiuso, “Core: Control-oriented regularization for system identification,” in 2018 IEEE Conference on Decision and Control (CDC). IEEE, 2018, pp. 2253–2258.
- [8] F. Dorfler, J. Coulson, and I. Markovsky, “Bridging direct & indirect data-driven control formulations via regularizations and relaxations,” IEEE Transactions on Automatic Control, 2022.
- [9] L. Campestrini, D. Eckhard, A. S. Bazanella, and M. Gevers, “Data-driven model reference control design by prediction error identification,” Journal of the Franklin Institute, vol. 354, no. 6, pp. 2628–2647, 2017.
- [10] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, “A data–driven approximation of the koopman operator: Extending dynamic mode decomposition,” Journal of Nonlinear Science, vol. 25, no. 6, pp. 1307–1346, 2015.
- [11] C. E. Rasmussen, “Gaussian processes in machine learning,” in Summer school on machine learning. Springer, 2003, pp. 63–71.
- [12] G. Pillonetto, F. Dinuzzo, T. Chen, G. De Nicolao, and L. Ljung, “Kernel methods in system identification, machine learning and function estimation: A survey,” Automatica, vol. 50, no. 3, pp. 657–682, 2014.
- [13] R. Pascanu, C. Gulcehre, K. Cho, and Y. Bengio, “How to construct deep recurrent neural networks,” arXiv preprint arXiv:1312.6026, 2013.
- [14] Y. Lian, R. Wang, and C. N. Jones, “Koopman based data-driven predictive control,” arXiv preprint arXiv:2102.05122, 2021.
- [15] H. J. van Waarde and R. Sepulchre, “Kernel-based models for system analysis,” arXiv preprint arXiv:2110.11735, 2021.
- [16] I. Lenz, R. A. Knepper, and A. Saxena, “Deepmpc: Learning deep latent features for model predictive control.” in Robotics: Science and Systems, vol. 10. Rome, Italy, 2015.
- [17] M. Korda and I. Mezić, “Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control,” Automatica, vol. 93, pp. 149–160, 2018.
- [18] J. Kocijan, R. Murray-Smith, C. E. Rasmussen, and A. Girard, “Gaussian process model based predictive control,” in Proceedings of the 2004 American control conference, vol. 3. IEEE, 2004, pp. 2214–2219.
- [19] Y. Chen, Y. Shi, and B. Zhang, “Optimal control via neural networks: A convex approach,” in International Conference on Learning Representations, 2018.
- [20] L. Hewing, J. Kabzan, and M. N. Zeilinger, “Cautious model predictive control using gaussian process regression,” IEEE Transactions on Control Systems Technology, vol. 28, no. 6, pp. 2736–2743, 2019.
- [21] E. T. Maddalena, P. Scharnhorst, Y. Jiang, and C. N. Jones, “Kpc: Learning-based model predictive control with deterministic guarantees,” in Learning for Dynamics and Control. PMLR, 2021, pp. 1015–1026.
- [22] B. Schölkopf, R. Herbrich, and A. J. Smola, “A generalized representer theorem,” in International conference on computational learning theory. Springer, 2001, pp. 416–426.
- [23] J. C. Willems, P. Rapisarda, I. Markovsky, and B. L. De Moor, “A note on persistency of excitation,” Systems & Control Letters, vol. 54, no. 4, pp. 325–329, 2005.
- [24] J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control: In the shallows of the DeePC,” in 2019 18th European Control Conference (ECC). IEEE, 2019, pp. 307–312.
- [25] C. De Persis and P. Tesi, “Formulas for data-driven control: Stabilization, optimality, and robustness,” IEEE Trans. Autom. Control, vol. 65, no. 3, pp. 909–924, 2019.
- [26] J. Berberich, J. Köhler, M. A. Muller, and F. Allgower, “Data-driven model predictive control with stability and robustness guarantees,” IEEE Trans. Autom. Control, vol. 66, no. 4, pp. 1702–1717, 2021.
- [27] J. Coulson, J. Lygeros, and F. Dörfler, “Regularized and distributionally robust data-enabled predictive control,” in 2019 IEEE 58th Conference on Decision and Control (CDC). IEEE, 2019.
- [28] L. Huang, J. Zhen, J. Lygeros, and F. Dörfler, “Quadratic regularization of data-enabled predictive control: Theory and application to power converter experiments,” arXiv preprint arXiv:2012.04434, 2020.
- [29] ——, “Robust data-enabled predictive control: Tractable formulations and performance guarantees,” arXiv preprint arXiv:2105.07199, 2021.
- [30] J. Berberich and F. Allgöwer, “A trajectory-based framework for data-driven system analysis and control,” in 2020 European Control Conference (ECC). IEEE, 2020, pp. 1365–1370.
- [31] J. G. Rueda-Escobedo and J. Schiffer, “Data-driven internal model control of second-order discrete volterra systems,” in 2020 59th IEEE Conference on Decision and Control (CDC). IEEE, 2020, pp. 4572–4579.
- [32] M. Alsalti, J. Berberich, V. G. Lopez, F. Allgöwer, and M. A. Müller, “Data-based system analysis and control of flat nonlinear systems,” arXiv preprint arXiv:2103.02892, 2021.
- [33] I. Markovsky, “Data-driven simulation of nonlinear systems via linear time-invariant embedding,” Vrije Universiteit Brussel, Tech. Rep, 2021.
- [34] I. Markovsky and F. Dörfler, “Identifiability in the behavioral setting,” 2020, available online.
- [35] I. Markovsky and P. Rapisarda, “Data-driven simulation and control,” International Journal of Control, vol. 81, no. 12, pp. 1946–1959, 2008.
- [36] J. Coulson, J. Lygeros, and F. Dorfler, “Distributionally robust chance constrained data-enabled predictive control,” IEEE Transactions on Automatic Control, 2021.
- [37] H. J. van Waarde, C. De Persis, M. K. Camlibel, and P. Tesi, “Willems’ fundamental lemma for state-space systems and its extension to multiple datasets,” IEEE Control Systems Letters, vol. 4, no. 3, pp. 602–607, 2020.
- [38] S. Saitoh and Y. Sawano, Theory of reproducing kernels and applications. Springer, 2016.
- [39] N. Aronszajn, “Theory of reproducing kernels,” Transactions of the American mathematical society, vol. 68, no. 3, pp. 337–404, 1950.
- [40] E. T. Maddalena, P. Scharnhorst, and C. N. Jones, “Deterministic error bounds for kernel-based learning techniques under bounded noise,” Automatica, vol. 134, p. 109896, 2021.
- [41] C. A. Micchelli, Y. Xu, and H. Zhang, “Universal kernels.” Journal of Machine Learning Research, vol. 7, no. 12, 2006.
- [42] E. C. Kerrigan and J. M. Maciejowski, “Soft constraints and exact penalty functions in model predictive control,” in Control 2000 Conference, Cambridge. Citeseer, 2000, pp. 2319–2327.
- [43] L. El Ghaoui and H. Lebret, “Robust solutions to least-squares problems with uncertain data,” SIAM Journal on matrix analysis and applications, vol. 18, no. 4, pp. 1035–1064, 1997.
- [44] A. Ben-Tal, D. den Hertog, and J.-P. Vial, “Deriving robust counterparts of nonlinear uncertain inequalities,” Mathematical Programming, vol. 149, no. 1, pp. 265–299, 2015.
- [45] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd, “OSQP: An operator splitting solver for quadratic programs,” Mathematical Programming Computation, pp. 1–36, 2020.
- [46] F. Dörfler and F. Bullo, “Synchronization and transient stability in power networks and nonuniform kuramoto oscillators,” SIAM Journal on Control and Optimization, vol. 50, no. 3, pp. 1616–1642, 2012.
- [47] S. D’Arco, J. A. Suul, and O. B. Fosso, “A virtual synchronous machine implementation for distributed control of power converters in smartgrids,” Elect. Power Syst. Res., vol. 122, pp. 180–197, 2015.
- [48] C. Yang, L. Huang, H. Xin, and P. Ju, “Placing grid-forming converters to enhance small signal stability of pll-integrated power systems,” IEEE Transactions on Power Systems, vol. 36, no. 4, pp. 3563–3573, 2021.