Learning-Based Efficient Approximation of Data-enabled Predictive Control
Yihan Zhou Yiwen Lu Zishuo Li Jiaqi Yan Yilin Mo




High-fidelity and Fault-tolerant Teleportation of a Logical Qubit using Transversal Gates and Lattice Surgery on a Trapped-ion Quantum Computer
    D. Hayes
    April 29, 2024
===============================================================================================================================================








Data-Enabled Predictive Control (DeePC) bypasses the need for system identification by directly leveraging raw data to formulate optimal control policies. However, the size of the optimization problem in DeePC grows linearly with respect to the data size, which prohibits its application due to high computational costs.
In this paper, we propose an efficient approximation of DeePC, whose size is invariant with respect to the amount of data collected, via differentiable convex programming. Specifically, the optimization problem in DeePC is decomposed into two parts: a control objective and a scoring function that evaluates the likelihood of a guessed I/O sequence, the latter of which is approximated with a size-invariant learned optimization problem. The proposed method is validated through numerical simulations on a quadruple tank system, illustrating that the learned controller can reduce the computational time of DeePC by 5x while maintaining its control performance.








§ INTRODUCTION


Optimal trajectory tracking has long been a fundamental problem in control systems. Classical model-based control approaches, such as Model Predictive Control (MPC), typically rely on a state-space model of the system to predict future trajectories and determine the optimal control decision. However, defining the state space model can be a laborious part of the control design <cit.>, and even after it is defined, identifying its parameters from noisy observations is known to be difficult <cit.>.

In recent years, direct data-driven control <cit.> has emerged as a promising solution to the aforementioned challenges, with Data-Enabled Predictive Control (DeePC) <cit.> being a representative algorithm encompassing this concept. A key idea of direct data-driven control is to construct a non-parametric model of the system using a data matrix formed by Input/Output (I/O) trajectories of the system collected a priori. According to Willem's fundamental lemma <cit.>, any feasible I/O trajectory of a noise-free linear system must be spanned by these offline trajectories, as long as the persistency of excitation condition <cit.> is satisfied, which justifies a non-parametric representation. Based on this non-parametric model, DeePC tackles the tracking problem by solving an optimization problem in a receding-horizon manner similarly to MPC, but with the non-parametric model replacing the state-space model for prediction.

DeePC has received significant attention from both theoretical <cit.> and empirical <cit.> perspectives. However, a notable challenge associated with the practical deployment of DeePC lies in its computation demand <cit.>. Since the I/O data collected from real-world systems are typically contaminated with noise, a large amount of data is required to generate an accurate non-parametric model of the system. However, as the amount of data used for building the nonparametric model increases, width of the data matrix grows accordingly, resulting in an increase in the dimension of a decision variable in DeePC.
Therefore, the computational burden of DeePC scales with the amount of data collected, which can be prohibitive for real-time applications such as robotics and power electronics <cit.>.

Previous research efforts have been devoted to mitigating the computational challenge of DeePC by compressing the data matrix. On the one hand, for lossless compression, a representation of predicted trajectories based on the null-space of the data matrix is proposed <cit.>. On the other hand, approximation methods have been developed to perform lossy compression of the data matrix, including those based on Singular Value Decomposition (SVD) <cit.>, Proper Orthogonal Decomposition (POD) <cit.>, and LQ decomposition <cit.>. Overall, these methods amount to diminishing the number of effective trajectory segments in the nonparametric model, i.e., reducing the width of the data matrix.

In this paper, we propose an alternative view of DeePC: we show that the DeePC problem is equivalent to minimizing the sum of a control objective ℓ(τ) and an anomaly metric S(τ), both as functions of the predicted I/O sequence τ. The former encompasses the tracking error, the penalty on control effort, and the constraints on I/O signals, which is independent from the data collected. The latter maps an I/O sequence to a scalar representing how unlikely τ is a trajectory generated by the system, whose evaluation is based on the data. We dub S(τ) as a scoring function or scoring model, since it can serve as an oracle whose queries return a score indicating the fitness of a guessed I/O sequence to the system dynamics. Although this reformulation is equivalent to the original DeePC, it motivates us to accelerate computation of the controller via a learning-based approximation Ŝ(τ) of the scoring function whose computational complexity is fixed even with more data collected. Specifically, we propose to parameterize the approximate scoring function as a differentiable convex program, and train its parameters via supervised learning against the true scoring function S.
With our proposed learning objective based on proximal operators, the minimizer of ℓ + Ŝ is close to the minimizer of ℓ + S, indicating that the approximate scoring function can be used to solve the overall control problem effectively.

An advantage of our approach is its computational merit in decoupling the scale of the control problem from the amount of data used. That is, once the form of the approximate scoring function is determined, one can train its parameters offline with a large dataset, without affecting the number of variables or constraints in the control problem to be solved online. This stands in contrast with the original DeePC, which, as mentioned, has a set of variables whose dimension scales with the amount of data. Hence, the learned approximate scoring function Ŝ can be viewed as a condensed representation of the system dynamics, which, combined with the control objective, enables computationally efficient control. The above intuitions are supported by numerical simulations, which show that our proposed controller can achieve a similar control performance to DeePC in less computational time.

Our contributions can be summarized as follows: i) We present a reformulation of DeePC based on the notion of scoring model, contributing to a new perspective on the data-enabled control problem. ii) We propose an approximate DeePC control law derived from a learning-based approximation of the scoring function. iii) We demonstrate through numerical simulations that the proposed method can significantly reduce the computational time required for DeePC, with only a minor compromise in control performance.

The rest of the paper is organized as follows. Section II provides a brief overview of the DeePC algorithm and its computational restriction. Section III presents our learning-based approximation of DeePC. Section IV presents the simulation results. Section V concludes the paper.

Notations:
The set of all integers is denoted by ℤ, with subsets of positive integers and non-negative integers represented as ℤ_>0 and ℤ_≥0, respectively. In the context of matrices and vectors, X^⊤ denotes the transpose of matrix X, and X^† represents the Moore-Penrose pseudo-inverse of matrix X. The operator diag(x_1, …, x_n) creates a diagonal matrix with elements x_1, …, x_n on its diagonal. The indicator function 𝕀_𝒳(x) equals 0 if x ∈𝒳 and +∞ otherwise. The gradient, subdifferential, and proximal operator of a function f are denoted by ∇ f, ∂ f, and Prox_f, respectively. The operator col is used to denote the column stack of vectors. I denotes the identity matrix of appropriate dimension or the identity operator according to the context. The p-norm of a vector or matrix or operator is denoted by ·_p.
The quadratic form x_Q^2 = x^⊤Qx signifies the weighted norm squared of vector x.
The symbol ∘ denotes the composition of two operators, ⊗ represents the Kronecker product, and ⟨·,·⟩ denotes the inner product. The symbol 𝒩(μ, Σ) denotes a Gaussian distribution with mean μ and covariance Σ.


§ PRELIMINARIES




 §.§ Problem Formulation


In this paper, we consider a discrete-time Linear Time-Invariant (LTI) system with state-space representation


    x(t+1)=Ax(t)+Bu(t)+w(t), 
       y(t)=Cx(t)+v(t),


where A ∈ℝ^n × n, B ∈ℝ^n × m, C ∈ℝ^p × n are state-space matrices, x(t)∈ℝ^n, u(t)∈ℝ^m, y(t)∈ℝ^p are state, input, and output vectors at time t ∈ℤ_≥ 0, respectively, and w(t) i.i.d.∼𝒩(0, Σ_w) and v(t) i.i.d.∼𝒩(0, Σ_v) are process and measurement noises. We assume that (A, B) is controllable and (A, C) is observable.

We consider an optimal tracking problem of the system (<ref>), whose objective is to minimize a quadratic function

    lim sup_T→∞1/T∑_t = 0^T-1(  y(t) - r _Q^2 +  u(t) _R^2 ),

subject to the constraints

    y(t) ∈𝒴, u(t) ∈𝒰, ∀ t ∈ℤ_≥ 0,

where r is a reference signal to be tracked, 𝒴⊆ℝ^p and 𝒰⊆ℝ^m are convex sets representing the output and input constraints respectively, and Q ,R are symmetric positive definite weighting matrices.


   The reference signal can be a time-varying signal r(t) for all of the methods in this paper, but we assume that r is a constant for notational simplicity.


We consider a setting where the state-space representation (<ref>) is unavailable to the system operator, i.e., none of the state-space matrices A, B, C, the state vectors x(t), or the state dimension n is known. Instead, only the input and output vectors u(t) and y(t) are observed.



 §.§ Data-Enabled Predictive Control (DeePC)


DeePC <cit.> tackles the aforementioned optimal tracking problem in a receding horizon manner, where the optimization problem described below is solved at every time step t. We use the convention that the index in the parentheses, e.g., y(t), denotes the actual value of the variable at time t., while the index in subscript, e.g., y_k, denotes the predicted value of the variable at the k-th step in the optimization horizon.

[DeePC Problem at Time Step t]  
   
    u_0:N-1, y_0:N-1, g, σ_ymin∑_k=0^N-1(y_k-r_Q^2+u_k_R^2) +
           λ_g1g_1 + λ_g2g_2^2  + λ_y1σ_y_1 + λ_y2σ_y_2^2 
        s.t. [ mT_ini{;     mN{; pT_ini{;     pN{ ][[ U_p; U_f; Y_p; Y_f ]]
             _ M
             g=[[   u_ini; u_0:N-1;   y_ini; y_0:N-1 ]]+[[   0;   0; σ_y;   0 ]], 
       90mu u_k∈𝒰, ∀ k ∈{0, …, N-1}, 
       90mu y_k∈𝒴, ∀ k ∈{0, …, N-1},

   where the new notations are defined and explained as follows: 
      
  * N ∈ℤ_>0 is the prediction horizon, M ∈ℤ_>0 is the number of trajectory segments collected offline, and T_ini≥ l is the length of the initialization sequence, where l is the lag <cit.> of the system.
      
  * The matrix ℋ≜[ U_p^⊤ U_f^⊤ Y_p^⊤ Y_f^⊤ ]^⊤ is the data matrix, each column of which is formed by concatenating the Input/Output (I/O) signals collected offline at T_ini + N consecutive time steps. We call such a column a trajectory segment. The data matrix has M columns, and is partitioned into row block matrices U_p, U_f, Y_p, Y_f, with number of rows mT_ini, mN, pT_ini, pN respectively, where the subscripts p and f stand for past and future, respectively. It is assumed that the matrix [ U_p^⊤ U_f^⊤ ]^⊤ is of full row rank.
      
  * u_ini∈ℝ^mT_ini and y_ini∈ℝ^pT_ini are column vectors formed by concatenating the I/O signals of the past T_ini time steps, which are inputs to the DeePC problem.
      
  * u_0:N-1∈ℝ^mN, y_0:N-1∈ℝ^pN stand for the predicted I/O sequence over the horizon N, which are optimization variables of the DeePC problem. After solving the DeePC problem, u_0^⋆ from the optimal solution is applied to the system, i.e., u(t) = u_0^⋆.
      
  * g ∈ℝ^M and σ_y ∈ℝ^pN are auxiliary optimization variables of the DeePC problem, which represent the coefficients for spanning the predicted I/O trajectory sequence from the data matrix and the slack variable for the output constraints, respectively. The scalar parameters λ_g1, λ_g2, λ_y1, λ_y2∈ℝ_≥ 0 are predefined coefficients. When brought together, λ_g1g_1 + λ_g2g_2^2  + λ_y1σ_y_1 + λ_y2σ_y_2^2 is viewed as a regularizer essential for the robustness of the controller in the presence of noises <cit.>.
      In different versions of DeePC, the regularizer form has been chosen to be l_1-norm <cit.>, squared l_2-norm <cit.>, or the hybrid of both <cit.>. Therefore, we adopt a general formulation that can cover the above cases.
   


Since the collected I/O data are contaminated with noise, a large amount of data is required for an accurate non-parametric representation of the system dynamics. However, the dimension M of the auxiliary optimization variable g grows linearly with respect to the amount of data, which poses a computational challenge to the online optimization of DeePC.



§ COMPUTATIONALLY EFFICIENT DEEPC VIA LEARNING-BASED APPROXIMATION


We propose a computationally efficient approximate DeePC by introducing the notion of the scoring function.
We first show that the DeePC objective can be viewed as minimizing the sum of the control cost and the score, and the evaluation of the latter is costly due to the scale of data matrix.
Therefore, we propose to approximate the scoring function with a reduced-order one via differentiable convex programming.
The parameters of the reduced scoring function are learned offline.
Finally, an approximate control law is formulated as minimizing the sum of the control cost and the learned approximate scoring function.
The overall framework is illustrated in Fig <ref>.





 §.§ DeePC through the Lens of Scoring Function


Let L=T_ini+N be the length of the I/O trajectory sequence in the data matrix. we denote 
    τ≜col(𝐮, 𝐲)
 as the predicted I/O trajectory sequence, where 𝐮∈ℝ^mL and 𝐲∈ℝ^pL are the predicted input and output sequences, respectively.

We propose to reformulate Problem <ref> and separate the objective into two parts:


    min _τℓ(τ)+S(τ),

where the objective functions are defined as

    ℓ(τ)  = 𝐲 - 𝐫_𝐐^2+𝐮_𝐑^2+ _𝒰^L ×𝒴^L(τ) + _{=τ_ini}(E_iniτ),


    S(τ)  = min _g, σ_yλ_g1g_1+λ_g2g_2^2+λ_y1σ_y_1+λ_y2σ_y_2^2 ,
           s.t. ℋ̃[[   g; σ_y ]]
       =τ,

and 𝐫≜ r⊗1_L, 𝐐≜ Q⊗ I_L, 𝐑≜ R ⊗ I_L, τ_ini≜col(u_ini, y_ini),

E_ini≜[ I 0 0 0; 0 0 I 0 ] is a selection matrix that extracts the initial part of I/O sequence from τ, ℋ̃≜[ U_p^⊤ Y_p^⊤ U_f^⊤ Y_f^⊤;     0    -I     0     0 ]^⊤.
In the notation of _^L×^L(τ), “×" denotes Cartesian power and 𝒰^L denotes L-ary Cartesian power of a set 𝒰.

Here, ℓ(·) and S(·) fulfill two orthogonal roles:



  * ℓ(·) represents the cost function of the control problem, which is distinct from and unaffected by the system dynamics.

  * S(·) is referred to as the . For any I/O trajectory τ, it evaluates the fitness of τ to the trajectory segments collected in the data matrix. It depends only on the system and is not influenced by any specific costs or constraints.


Note that g, σ_y are auxiliary variables required for the evaluation of the scoring function, not part of the control policy,
and the dimension of g grows with the number of trajectory segments in the data matrix.
Therefore, we desire to accelerate the evaluation of S(τ) using a learning-based approach which approximates the  with fewer auxiliary variables.
We denote the approximate scoring function as Ŝ(τ), then the overall control problem becomes:

    min_τℓ(τ)+Ŝ(τ).




 §.§ Learning Objective


Since the approximation goal is to make control input, i.e., the optimal solution, of the learned problem min l(τ) + Ŝ(τ) close to that of min l(τ) + S(τ), we set the learning target to be minimizing the error between the proximal operators of the approximate and true scoring functions Ŝ and S. 

The proximal operator <cit.> of Ŝ is defined as:

    Prox_Ŝ(τ)=min _τ̂Ŝ(τ̂)+1/2τ̂-τ^2.

Consequently, the optimality condition of (<ref>), i.e, 0 ∈ (∂ℓ + ∂Ŝ)(τ̂^⋆), can be equivalently expressed with the proximal operators of ℓ and Ŝ using operator splitting methods. Therefore, the optimal solution τ̂^⋆ to (<ref>) can be determined once Prox_Ŝ is given.



 §.§ Approximator Form: Differentiable Convex Program


Given the learning target Prox_S, the subsequent step is to select a parameterized family of Ŝ that satisfies the following conditions:
i) Ŝ(τ) is capable of approximately representing the true  S(τ);
ii) Ŝ is a convex function of τ, and Prox_Ŝ can be easily parameterized; iii) the control problem min_τℓ(τ)+Ŝ(τ) should be computationally tractable.
Considering the above requirements, we adopt a differentiable convex programming approach to represent Ŝ.
Specially, we adopt the following form:

   
    Ŝ(τ)=
       min_z ∈ℝ^n_z   diag(d_1)z_1 + diag(d_2)z_2^2 
    s.t.     Gz + Wτ = 0.



Here, Ŝ is a mapping from an I/O sequence τ to a scalar standing for the optimal value of a convex optimization problem with n_z variables and m_z constraints, whose decision variable is z and whose parameters include τ as well as learnable parameters d_1 ∈ℝ^n_z, d_2∈ℝ^n_z, G ∈ℝ^m_z × n_z, W ∈ℝ^m_z × L(m+p). To train the differentiable convex program means to update the parameters d_1, d_2, G, W via gradient-based methods, such that the learning target (Prox_Ŝ≈Prox_S in our case) is achieved. The size of the convex program, i.e., n_z and m_z, can be arbitrarily chosen according to a trade-off between the approximation accuracy and the computational cost, but once fixed, it is invariant against the amount of data used for training.

We first show that the chosen form of Ŝ has adequate representational power to approximate the true  S:


   With sufficiently large n_z and m_z, the true S(τ) is representable by (<ref>).



   Let d_1 = [λ_g1, …, λ_g1_M times, λ_σ_y1, …, λ_σ_y1_pT_ini times]^⊤,
   d_2 = [√(λ_g2), …, √(λ_g2)_M times, √(λ_σ_y2), …, √(λ_σ_y2)_pT_ini  times]^⊤, G = ℋ̃, W = -I, then Ŝ(τ) defined in (<ref>) is equivalent to S(τ).


We next show how to evaluate and learn the proximal operator of the function Ŝ given in (<ref>). By definition (<ref>) of proximal operator, we have:
   
    (z^⋆, τ̂^⋆) = min_z, τ̂diag(d_1)z_1 + diag(d_2)z_2^2+ 
       100mu  τ̂- τ^2/2,  s.t.  Gz + Wτ̂= 0. 
       Prox_Ŝ(τ) = τ̂^⋆.


To evaluate (<ref>) and find the gradient of τ^⋆ with respect to the learnable parameters d_1, d_2, G, W, we adopt an algorithm unrolling <cit.> approach, which has been proved efficient for controller learning <cit.>.
In particular, we unroll the Douglas-Rachford Splitting (DRS) <cit.>, which, when applied to the problem (<ref>), yields the following iterations:

   
    [[  z^k+1/2; τ̂^k+1/2 ]] = [[ sh_d_1, d_2(ξ^k);        τ + η^k/2 ]], 
       [[  z^k+1; τ̂^k+1 ]] = (I - G̃^†G̃)[[  2z^k+1/2-ξ^k; 2τ̂^k+1/2-η^k ]], 
       [[ ξ^k+1; η^k+1 ]] = [[     ξ^k+z^k+1-z^k+1/2; η^k + τ̂^k+1-τ̂^k+1/2 ]],

   

where ξ, η are auxiliary variables, G̃ = [ G W ], and the soft-thresholding operator sh_d_1, d_2(x) is defined as:

    sh_d_1, d_2(x)_i= x_i-|d_1 i|/1+2 (d_2 i)^2    if  x_i-|d_1 i|>0, 
    x_i+|d_1 i|/1+2 (d_2 i)^2    if  x_i+|d_1 i|<0, 
    
       0     otherwise .

Here, d_1 i and d_2 i are the i-th elements of d_1 and d_2, respectively.

The unrolled iterations resemble a deep neural network that interleaves affine transformations with soft-thresholding nonlinearities, whose weights are functions of the learnable parameters d_1, d_2, G, W. This resemblance, visualized in Fig <ref>, facilitates the training of the learnable parameters using established deep learning frameworks.





 §.§ Learning Algorithm


Given a data matrix ℋ consisting of I/O trajectory sequences, we first collect a dataset 𝒟 of I/O trajectory sequences by linear combination of the trajectories in ℋ and adding random noise to the trajectories to increase the diversity of the dataset. Then, we compute the ground truth Prox_S(τ) for each trajectory τ in 𝒟. Finally, we learn the parameters θ = (d_1, d_2, G, W) of the approximate scoring function Ŝ by minimizing the mean squared error between Prox_Ŝ and the ground truth Prox_S across all entries in the dataset 𝒟.

The learning procedure of the  is summarized in Algorithm <ref>.

 

Note that learning the approximate scoring function is an offline process, and the learned parameters θ are used to substitute the original scoring function S(τ) to solve the overall control problem (<ref>) online.



§ SIMULATION


In this section, we present the simulation results to demonstrate the effectiveness of our approach.
All computations are performed on a computer with an AMD Ryzen 9 5900X Processor clocked at 3.7GHz.
All controllers are implemented using the MOSEK optimizer.
Our code is available at <https://github.com/zhou-yh19/redpc>.

Consider a quadruple-tank system <cit.> whose simulation employs the same linearized system equations and system parameters as the LTI system used in <cit.>:

    A = [ 0.921     0 0.041     0;     0 0.918     0 0.033;     0     0 0.924     0;     0     0     0 0.937 ], 
       B = [ 0.017 0.001; 0.001 0.023;     0 0.061; 0.072     0 ],
    C = [ 1 0 0 0; 0 1 0 0 ].


The system is subject to process noise w(t) ∼𝒩(0, 0.01I_4) and measurement noise v(t) ∼𝒩(0, 0.1I_2).

We follow the same setting as in <cit.>, such as the control cost matrices Q = 35· I_2 and R = 10^-4I_2, the control input and output constraints 𝒰 = 𝒴 = [-2, 2], the setpoint r(t) = [ 0.65,  0.77 ]^T, the prediction horizon N = 20, and the initial sequence length T_ini = 10.


Data Collection: The dataset is collected by applying random input sequences to the system. The input sequences are generated by sampling from a uniform distribution over 𝒰. The corresponding output sequences are obtained by simulating the system with the input sequences. In total, 1500 I/O data points are collected to construct the data matrix ℋ in the form of a Hankel matrix.

We choose the regularization parameters of DeePC as λ_g1 = 1, λ_g2 = 100, λ_y1 = 100, λ_y2 = 100000, and incorporates all 1500 I/O data points into the data matrix. And the approximate scoring function is learned with 55 equations, 110 variables and 200 DRS iterations. The training procedure is completed within 260s on a single NVIDIA GeForce RTX 4090 GPU.
Fig <ref> illustrates the first two states of the system under the control of both the DeePC controller and our proposed approach. The results demonstrate that both controllers track the setpoint effectively.







In Table <ref>, we benchmark the control performance of our method with DeePC and MPC, detailing the sample average accumulative cost 1/K∑_k∑_t(y(t)-r(t)_Q^2+u(t)_R^2), the average and worst computation time per step. Our approach achieves control performance comparable to MPC and DeePC, while significantly enhancing computational efficiency.



Furthermore, we vary the number of I/O sequences in the data matrix of DeePC and the size of the learned approximate scoring function in our method to examine their effects on control performance.
The results are shown in Fig. <ref>.
For DeePC, an increase in the number of I/O sequences in the data matrix results in better control performance but also longer computation time.
For our approach, the computation time is mainly determined by the size of the learned approximate scoring function, and there is also trade-off between control performance and computation time.
Nevertheless, it can be observed that under the same computational time constraints, our method achieves superior control performance compared to DeePC, and when aiming for comparable control quality, our approach significantly reduces the computation time. This demonstrates the effectiveness of our approach in enhancing the computational efficiency of DeePC while maintaining control performance. 







§ CONCLUSIONS


In this paper, we propose a computationally efficient approximate DeePC by introducing the notion of the scoring function and replacing it with a learned approximate one via differentiable convex programming. The parameters of the reduced scoring function are learned offline, and the control law is formulated as minimizing the sum of the control cost and the learned approximate scoring function. Simulation results demonstrate our approach can achieve a comparable tracking performance to DeePC while significantly reducing the computation time. In the future, we will apply our approach to the real-world control problems and investigate the generalization ability of the learned approximate scoring function.


IEEEtran


