# Data enabled Predictive Control of LPV systems

Taleb Bou Hamdan, Patrick Coirault, Guillaume Mercère, Thibault Dairay

► To cite this version:

Taleb Bou Hamdan, Patrick Coirault, Guillaume Mercère, Thibault Dairay. Data enabled Predictive Control of LPV systems. Control Engineering Practice, 2024, 24, 149, 10.1016/j.conengprac.2024.105969. hal-04582708v2

HAL Id: hal-04582708

https://hal.science/hal-04582708v2

Submitted on 22 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Data enabled Predictive Control of LPV systems

Taleb Bou Hamdan<sup>a,b</sup>, Patrick Coirault<sup>a,*</sup>, Guillaume Mercère<sup>a</sup>, Thibault Dairay<sup>b,c</sup>

## ARTICLE INFO

### Keywords:

Data enabled predictive control
Linear Parameter Varying Systems
Willem's lemma

## ABSTRACT

This paper introduces a data-enabled predictive control methodology designed for Linear Parameter Varying (LPV) systems. By leveraging a polytopic representation of the LPV system, we formulate Willem's lemma to accommodate parameter variations. The system trajectory is predicted over a finite horizon using specific trajectories generated offline. A notable advantage of this approach is its independence from a priori knowledge of parametric variations. Simulation studies conducted on both a numerical example and a simulator of a calendering process governed by partial differential equations substantiate the effectiveness of this approach.

## 1. Introduction

The effectiveness of Model Predictive Control (MPC) in industrial system control is well-established, attributed to its favorable balance between performance and ease of tuning. Notably, MPC adeptly handles Multi-Input Multi-Output (MIMO) systems, addressing both equality and inequality constraints on control signal and output. Despite its known computational complexity, the diminishing costs of computational resources, especially on-board computers, have further bolstered the adoption of MPC. A significant challenge in deploying MPC lies in obtaining a representative model. This model enables MPC to predict system outputs over a finite horizon based on control variables. Typically, this model is derived through an identification step reliant on input/output measurements from a properly excited system (Van Over-schee & De Moor, 1994, 1996). The execution of this identification step demands proficiency in signal processing, statistics and optimization, making its accurate implementation potentially complex. Consequently, the ability to represent system dynamics without the prerequisite of the model identification step holds substantial value in this context.

Addressing this demand, recent years have witnessed the emergence of Data Driven methods. Utilizing Willem's fundamental lemma, these methods offer a solution to the predictive control problem by relying solely on input/output sequences (Willems, 2007; Willems, Rapisard, Markovsky, & De Moor, 2005), with extensions found in Dinkla, Mulders, van Wingerden, and Oomen (2023), Markovsky and Dörfler (2021) and Van Waarde, De Persis, Camlibel, and Tesi (2020). The key concept revolves around the notion that a trajectory (an input/output sequence) from a Linear Time-Invariant (LTI) system, obtained through sufficiently rich excitation, inherently contains all

relevant information on the system dynamics *a priori*. This trajectory can effectively serve as a parametrization tool for any other trajectory of the system. The extensive literature on this approach, known as Data Driven Predictive Control or Data Enabled Predictive Control (DeePC), has predominantly focused on discrete-time LTI systems, particularly exploring its robustness in the presence of noise (Breschi, Chiuso, & Formentin, 2023a; Markovsky & Dörfler, 2021). This article principal contribution lies in extending this approach to Linear Parameter Varying (LPV) systems.

Linear Parameter Varying (LPV) systems are linear dynamic systems with a state-space representation dependent on nonstationary exogenous parameters (Morato, Normey-Rico, & Sename, 2020; Toth, 2010). Conceptually, an LPV system can be visualized as a collection of Linear Time-Invariant (LTI) systems interconnected by these exogenous parameters, resulting in nonstationary dynamics. This framework enables the modeling of system nonlinearities through a series of LTI models (Gidon et al., 2021; Kapsalis, Sename, Milanes, & Molina, 2022; Li, Nguyen, Guerra, & Kruszewski, 2023). Leveraging this approach, control methods originally designed for linear systems can be applied to manage nonlinear counterparts effectively. Given the close proximity of LPV modeling to LTI modeling, it is logical to explore the conditions under which the DeePC approach, initially developed for LTI systems, can be extended to LPV systems. Recent studies have delved into the application of Willem's fundamental Lemma to LPV systems (Verhoek, Abbas, Toth, & Haesaert, 2021). However, a notable limitation arises in these investigations as the exogenous parameters must be known and measurable at every moment, posing a significant drawback to LPV modeling.

<footnote><sup>a</sup> Laboratoire d'Automatique et d'Informatique pour les Systèmes, Poitiers University, Poitiers, France</footnote>

<footnote><sup>b</sup> Manufacture Française des Pneumatiques Michelin, Clermont-Ferrand, France</footnote>

<footnote><sup>c</sup> Centre Borelli, CNRS, Université Paris Saclay, ENS Paris Saclay, Gif-sur-Yvette, France</footnote>

<footnote>* Correspondence to: Bâtiment B25, 2 rue Pierre Brousse, TSA 41105, 86073 Poitiers cedex 9, France.
E-mail address: patrick.coirault@univ-poitiers.fr (P. Coirault).</footnote>
The main contribution of this paper is the introduction of an extension of DeePC to systems represented by an LPV model, obviating the requirement to measure exogenous parameters. We consider in this paper the noise free case only. In achieving this objective, output prediction over a finite horizon is conducted based on a set of specific trajectories acquired offline. Key distinctions from the LTI case include the number of trajectories needed to capture system dynamics and the methodology used for generating these trajectories. Although the optimization problem remains quadratic, its size expands due to the increased number of decision variables corresponding to the higher number of trajectories.

The subsequent sections of this paper are organized as follows: Section 2 provides an overview of notations and the model structure. Section 3 initially revisits Data Enabled Predictive Control for Linear Time-Invariant systems and then introduces a novel methodology for applying DeePC to Linear Parameter Varying systems without prior knowledge of the varying parameter trajectories. Section 4 evaluates the applicability and efficiency of the introduced methodology through two numerical examples, including an LPV system with two varying parameters and a calendering process. The paper concludes with a summary of remarks in Section 5.

## 2. Preliminaries

This section presents the notations used in the paper, the model structure and useful definitions.

### 2.1. Notations and definitions

Let $\mathbb{N}$ and $\mathbb{R}$ be the sets of positive integers and real numbers, respectively. $\mathbb{N}^*$ denotes the set of positive non-zero integers. Let $I_{k,s} = \{k, k+1, \dots, k+s-1\}$. The set of real column vectors of dimension $n \in \mathbb{N}^*$ is denoted by $\mathbb{R}^n$ and the set of real matrices of $n \in \mathbb{N}^*$ rows and $m \in \mathbb{N}^*$ columns is denoted by $\mathbb{R}^{m \times m}$. Given a rectangular matrix $\mathbf{A} \in \mathbb{R}^{m \times m}$, its transpose is denoted by $\mathbf{A}^\top \in \mathbb{R}^{m \times n}$, $\mathbf{A}^{(i)} \in \mathbb{R}^{1 \times m}$ represents its $i$th row. For a vector $\mathbf{x}(k) \in \mathbb{R}^{n_x}$, $\Delta\mathbf{x}(k) = \mathbf{x}(k) - \mathbf{x}(k-1)$. For any vector $\mathbf{x}(k) \in \mathbb{R}^{n_x}$, with $k \in \mathbb{N}$, the finite vector over a specific window of size $\rho$ steps ($\rho \in \mathbb{N}^*$) starting from a specified instant $k \in \mathbb{N}$ is denoted as

$$ \mathbf{X}_{k,\rho,1} = \begin{bmatrix} \mathbf{x}(k) \\ \mathbf{x}(k+1) \\ \vdots \\ \mathbf{x}(k+\rho-1) \end{bmatrix} \in \mathbb{R}^{n_x \rho}. \quad (1) $$

Accordingly, the block Hankel matrix containing the available data starting from instant $k \in \mathbb{N}$ distributed over $\rho \in \mathbb{N}^*$ rows and $M \in \mathbb{N}^*$ columns is denoted as

$$ \mathbf{X}_{k,\rho,M} = [\mathbf{X}_{k,\rho,1} \quad \mathbf{X}_{k+1,\rho,1} \quad \dots \quad \mathbf{X}_{k+M-1,\rho,1}] \in \mathbb{R}^{n_x \rho \times M}. \quad (2) $$

**Definition 1.** An input sequence $\{\mathbf{u}(k)\}_{k=0}^{N-1}$, $\mathbf{u}(k) \in \mathbb{R}^{n_u}$, is persistently exciting of order $\rho$ if the matrix $U_{0,\rho,N}$ is of rank $n_u \rho \le N \in \mathbb{N}$ (see Verhaegen and Verdult (2007, chapter 10)).

The norm of the vector $\|\mathbf{X}_{k,\rho,1}\|_Q^2$ denotes the quadratic form $\mathbf{X}_{k,\rho,1}^\top \mathbf{Q} \mathbf{X}_{k,\rho,1}$ where $\mathbf{Q} \in \mathbb{R}^{n_x \rho \times n_x \rho}$ is a symmetric and strictly positive matrix. The following matrices are defined by

$$ S_{\ell,n} = \begin{bmatrix} \mathbb{1}_{n \times n} & \mathbb{0}_{n \times n} & \dots & \mathbb{0}_{n \times n} \\ \mathbb{0}_{n \times n} & \mathbb{1}_{n \times n} & \mathbb{0}_{n \times n} & \dots \\ \vdots & \vdots & \ddots & \vdots \\ \mathbb{0}_{n \times n} & \mathbb{0}_{n \times n} & \mathbb{0}_{n \times n} & \mathbb{1}_{n \times n} \end{bmatrix} \in \mathbb{R}^{\ell n \times \ell n}, \quad (3) $$

$$ 1_{\ell,n} = \begin{bmatrix} \mathbb{1}_{n \times n} \\ \vdots \\ \mathbb{1}_{n \times n} \end{bmatrix} \in \mathbb{R}^{\ell n \times n}. \quad (4) $$

### 2.2. Model structure

This section introduces a behavioral framework for LPV systems (see Toth (2010) for more details).

**Definition 2.** A discrete-time parameter varying dynamical system $\mathcal{S}$ is defined as a quadruple

$$ \mathcal{S} = (\mathbb{T}, \mathbb{P}, \mathbb{W}, \mathcal{B}) \quad (5) $$

with $\mathbb{T} = \mathbb{Z}$ the discrete time axis, $\mathbb{P} \subseteq \mathbb{R}^{n_p}$ the scheduling space, $\mathbb{W} = \mathbb{R}^{n_w}$ the signal space and $\mathcal{B} \subseteq (\mathbb{W} \times \mathbb{P})^\top$ the behavior ($\mathbb{X}^\top$ is the standard notation for the collection of all maps from $\mathbb{T}$ to $\mathcal{X}$)

$\mathbb{p} \in \mathbb{P}^\top$ is the scheduling variable, an external free signal of the system that governs the dynamical behavior. In the sequel, we assume that $\mathbb{p}$ evolves in a polytope of dimension $N_p = 2^{n_p}$. Inside the polytope, the admissible trajectories of $\mathbb{p}$ are restricted to a subset of $\mathbb{P}^\top$, the set of admissible scheduling trajectories. This set, denoted by $\mathcal{B}_{\mathbb{P}}$, is described as follows

$$ \mathcal{B}_{\mathbb{P}} = \{ \mathbb{p} \in \mathbb{P}^\top \mid \exists \mathbf{w} \in \mathbb{W}^\top \text{ s.t. } (\mathbf{w}, \mathbb{p}) \in \mathcal{B} \}. \quad (6) $$

Additionally, for a given trajectory $\mathbb{p} \in \mathcal{B}_{\mathbb{P}}$, the projected behavior

$$ \mathcal{B}_{\mathbb{p}} = \{ \mathbf{w} \in \mathbb{W}^{\top} \mid (\mathbf{w}, \mathbb{p}) \in \mathcal{B} \} \quad (7) $$

defines all the signal trajectories compatible with the fixed scheduling trajectory $\mathbb{p}$, In case of a constant scheduling trajectory $\mathbb{p} \in \mathcal{B}_{\mathbb{p}}$, $\mathbb{p}(k) = \mathbb{p}_f \quad \forall k \in \mathbb{Z}$, the projection behavior $\mathcal{B}_{\mathbb{p}_f}$ is called a frozen behavior, i.e.,

$$ \mathcal{B}_{\mathbb{p}_f} = \{ \mathbf{w} \in \mathbb{W}^{\top} \mid (\mathbf{w}, \mathbb{p}) \in \mathcal{B} \text{ with } \mathbb{p}(k) = \mathbb{p}_f \quad \forall k \in \mathbb{Z} \}. \quad (8) $$

In this case, $\mathcal{F}_{\mathbb{p}_f} = (\mathbb{T}, \mathbb{W}, \mathcal{B}_{\mathbb{p}_f})$ defines a frozen system of $\mathcal{S}$ for each $\mathbb{p}_f \in \mathbb{P}$. $\mathcal{S}$ is called Linear Parameter Varying (LPV) if the following conditions are satisfied:

* $\mathbb{W}$ is a vector space and $\mathcal{B}_{\mathbb{p}}$ is linear subspace of $\mathbb{W}^\top \forall \mathbb{p} \in \mathcal{B}_{\mathbb{P}}$,

* $\mathbb{T}$ is closed under addition,

* for any $(\mathbf{w}, \mathbb{p}) \in \mathcal{B}$ and any $k \in \mathbb{Z}$, it holds that $(\mathbf{w}.(k), \mathbb{p}.(k)) \in \mathcal{B}$.

A discrete time state-space representation of an LPV system, denoted as $\mathcal{S} = (\mathbf{p}, \mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})$, is given by

$$ \mathbf{x}(k+1) = \mathbf{A}(\mathbb{p}(k))\mathbf{x}(k) + \mathbf{B}(\mathbb{p}(k))\mathbf{u}(k), \quad (9a) $$

$$ \mathbf{y}(k) = \mathbf{C}(\mathbb{p}(k))\mathbf{x}(k) + \mathbf{D}(\mathbb{p}(k))\mathbf{u}(k), \quad (9b) $$

where $\mathbf{x}(k) \in \mathbb{R}^{n_x}$ is the state vector of the system at each instant $k \in \{0, \dots, n_t\}$ ($n_t$ is the concerned time domain and $k$ stands for the discrete time step), $\mathbf{u}(k) \in \mathbb{R}^{n_u}$ and $\mathbf{y}(k) \in \mathbb{R}^{n_y}$ are the input and output vectors of the system, respectively. $\mathbf{A} \in \mathbb{R}^{n_x \times n_x}$, $\mathbf{B} \in \mathbb{R}^{n_x \times n_u}$, $\mathbf{C} \in \mathbb{R}^{n_y \times n_x}$ and $\mathbf{D} \in \mathbb{R}^{n_y \times n_u}$ are the system, input, output and feedthrough matrices, respectively. $(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})$ are matrix functions with static dependence on $\mathbb{p} \in \mathbb{R}^{n_p}$.

As $\mathbb{p}$ lies in a polytope, each component $\mathbb{p}_h$ of $\mathbb{p}$ is bounded, i.e.,

$$ \underline{\mathbb{p}}_h \le \mathbb{p}_h \le \bar{\mathbb{p}}_h \quad \forall h \in \{1, n_p\}. $$

For LPV systems, unlike Linear Time Varying (LTV) systems, the variation of the system dynamics is not directly linked to time but is instead associated with the variation of the parameter vector $\mathbb{p}$.

The class of LTI systems can be derived from the LPV ones through the following approach.

**Definition 3.** Let $\mathcal{S} = (\mathbb{T}, \mathbb{P}, \mathbb{W}, \mathcal{B})$ be an LPV system. The set of LTI systems

$$ \mathcal{F}_f = \{ \mathcal{F} = (\mathbb{T}, \mathbb{W}, \mathcal{B}') \mid \exists \mathbb{p}_f \in \mathbb{P} \text{ s.t. } \mathcal{B}' = \mathcal{B}_{\mathbb{p}_f} \} $$

is called the frozen set of $\mathcal{S}$.
It follows that the structure of an element of $\mathcal{F}_s$ is given by

$$ \mathbf{x}(k+1) = \mathbf{A}\mathbf{x}(k) + \mathbf{B}\mathbf{u}(k), \quad (10a) $$

$$ \mathbf{y}(k) = \mathbf{C}\mathbf{x}(k) + \mathbf{D}\mathbf{u}(k). \quad (10b) $$

This LTI system is denoted by $\mathcal{S}(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})$. The following assumptions are considered for each element of $\mathcal{F}_s$, i.e. an LTI system $\mathcal{S}(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})$.

**Assumption 1.** The pair $(\mathbf{A}, \mathbf{B})$ is controllable.

**Assumption 2.** The pair $(\mathbf{A}, \mathbf{C})$ is observable.

### 3. Data enabled predictive control

This section first recalls the DeePC problem for deterministic LTI systems. It begins with a review of Willem's fundamental lemma and its application to predictive control. Subsequently, an extension to deterministic LPV systems is introduced. Finally, a systematic approach to embed an integral action into data driven predictive control is presented. These two last aspects are the primary contributions of this paper.

#### 3.1. Trajectory based representation for LTI systems

To implement Predictive Control, it is necessary to be able to calculate the derivatives of a constraint cost function with respect to the optimization variables over a predictive horizon. This involves expressing the output and/or input as a function of these optimization variables (Toth, 2010). The classic approach is to use a model. This model can be deduced from a transfer function, a state-space realization or from a subspace approach (Favoreel, Moor, & Gevers, 1999). Obtaining this model requires an identification step which can be complex to implement. In contrast, data driven approaches do away with the model and the identification step, starting from the idea that all the information on the dynamics of the system is contained in the input-output data. This is the main idea of Willem's fundamental lemma which is reviewed in Lemma 1.

The notion of system trajectory plays a key role in this framework. An output-input sequence $\{\mathbf{y}(k), \mathbf{u}(k)\}_{k=0}^{N-1}$ for $N \in \mathbb{N}^*$ is a trajectory of the LTI system $\mathcal{S}$ if it satisfies Definition 4.

**Definition 4.** An output-input sequence $\{\mathbf{y}(k), \mathbf{u}(k)\}_{k=0}^{N-1}$ for $N \in \mathbb{N}^*$ (represented by $\{\mathbf{Y}_{0,N,1}, \mathbf{U}_{0,N,1}\}$) is a trajectory of the LTI system $\mathcal{S}(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})$ if an initial condition $\mathbf{x}(0) \in \mathbb{R}^{n_x}$ and a state sequence $\{\mathbf{x}(k)\}_{k=0}^N$ exist such that the sequence $\{\mathbf{x}(k)\}_{k=0}^N$ satisfies the state space equations of the system given in Eq. (10).

Starting with the LTI system (10) and from Definition 4, Lemma 1 gives a trajectory based representation which is known as Willem's Lemma (Willems, 2007), or the Fundamental lemma as established in Berberich and Allgöwer (2020).

**Lemma 1.** Consider an LTI system $\mathcal{S}(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})$. Assume that the input signal used to build the Hankel matrix $\mathbf{U}_{0,s,N}$ is persistently exciting of order $s$ with $s \in \mathbb{N}^*$. Then, for $N \in \mathbb{N}^*$, the system output-input trajectory at each instant $k \in \mathbb{N}$ over $s$ steps $[\mathbf{Y}_{k,s,1}^\top \quad \mathbf{U}_{k,s,1}^\top]^\top$ belongs to the column space of the matrix $[\mathbf{Y}_{0,s,N}^\top \quad \mathbf{U}_{0,s,N}^\top]^\top$.

In other words, there exists a vector $\alpha(k) \in \mathbb{R}^N$ such that

$$ \begin{bmatrix} \mathbf{Y}_{k,s,1} \\ \mathbf{U}_{k,s,1} \end{bmatrix} = \begin{bmatrix} \mathbf{Y}_{0,s,N} \\ \mathbf{U}_{0,s,N} \end{bmatrix} \alpha(k). \quad (11) $$

The key result of Lemma 1 is that any output-input trajectory of the LTI system $\mathcal{S}$ can be embedded into a linear combination of sufficiently excited past output-input trajectory.

#### 3.2. Data driven predictive control for LTI systems

The main objective behind the predictive control approach is to find, at each time $k$, the optimal trajectory $[\mathbf{Y}_{k,L,1} \quad \mathbf{U}_{k,L,1}]$ over the prediction horizon $L \in \mathbb{N}^*$ which minimizes the quadratic cost function

$$ J(\mathbf{Y}_{k,L,1}, \mathbf{U}_{k,L,1}) = \|\mathbf{Y}_{k,L,1} - \mathbf{Y}_r\|_Q^2 + \|\mathbf{U}_{k,L,1}\|_R^2, \quad (12a) $$

$$ \text{s.t. } \mathbf{u}(i) \in \mathcal{V}, \quad (12b) $$

$$ \mathbf{y}(i) \in \mathcal{Y}, \quad (12c) $$

$\forall i \in I_{k,L}$, where $\mathbf{Q} \in \mathbb{R}^{n_y L \times n_y L}$ and $\mathbf{R} \in \mathbb{R}^{n_u L \times n_u L}$ are user-defined output and input error penalizing positive definite matrices. They are tuned based on a trade-off between the degree of importance of each of the outputs and inputs terms (Rawlings, Mayne, & Diehl, 2017). $\mathcal{Y}$ and $\mathcal{V}$ are polytopes defining the applicable lower and upper boundaries of the system output and input respectively. $\mathbf{Y}_r \in \mathbb{R}^{n_y L}$ stands for the reference trajectory over the prediction horizon such that

$$ \mathbf{Y}_r = \begin{bmatrix} \mathbf{y}_r(k) \\ \vdots \\ \mathbf{y}_r(k + L - 1). \end{bmatrix} \quad (13) $$

Splitting up the trajectories between past and future trajectories as well as the Hankel matrices as shown in Fig. 1 and assuming that $\text{rank}(\mathbf{U}_{0,\rho+L,N}) = \rho + L$, a vector $\alpha(k) \in \mathbb{R}^N$ exists such that

$$ \begin{bmatrix} \mathbf{Y}_{k-\rho,\rho,1} \\ \mathbf{Y}_{k,L,1} \\ \mathbf{U}_{k-\rho,\rho,1} \\ \mathbf{U}_{k,L,1} \end{bmatrix} = \begin{bmatrix} \mathbf{Y}_{0,\rho,N} \\ \mathbf{Y}_{\rho,L,N} \\ \mathbf{U}_{0,\rho,N} \\ \mathbf{U}_{\rho,L,N} \end{bmatrix} \alpha(k). \quad (14) $$

Referring to Eq. (12), the decision variables are the future output and input vectors that can be written in terms of $\alpha(k)$ based on Eq. (14) such that

$$ \mathbf{Y}_{k,L,1} = \mathbf{Y}_{\rho,L,N} \alpha(k), \quad (15a) $$

$$ \mathbf{U}_{k,L,1} = \mathbf{U}_{\rho,L,N} \alpha(k). \quad (15b) $$

Using the result of Eq. (15), the cost function (12) can be reformulated with $\alpha(k)$ as the sole decision variable. The DeePC problem boils down to finding $\alpha(k)$ such that

$$ \arg \min_{\alpha(k)} (\|\alpha(k)\|_{\mathbf{W}}^2 + 2\mathbf{C}_{\rho}^{\top}\alpha(k)), \quad (16) $$

subject to

$$ \begin{bmatrix} \mathbf{Y}_{k-\rho,\rho,1} \\ \mathbf{U}_{k-\rho,\rho,1} \end{bmatrix} = \begin{bmatrix} \mathbf{Y}_{0,\rho,N} \\ \mathbf{U}_{0,\rho,N} \end{bmatrix} \alpha(k), \quad (17a) $$

$$ \mathbf{Y}_{L,\rho,N} \alpha(k) = \begin{bmatrix} \mathbf{y}_r(k + L - \rho) \\ \vdots \\ \mathbf{y}_r(k + L - 1) \end{bmatrix}, \quad (17b) $$

$$ \mathbf{U}_{\rho,L,N}^{(i)} \alpha(k) \in \mathcal{V}, \quad (17c) $$

$$ \mathbf{Y}_{\rho,L,N}^{(i)} \alpha(k) \in \mathcal{Y}, \quad (17d) $$

where $i \in I_{1,L}$. The terminal equality constraint is added to guarantee closed loop stability (Berberich, Köhler, Müller, & Allgöwer, 2021). $\mathbf{W}$ and $\mathbf{C}_\rho$ are given as

$$ \mathbf{W} = \mathbf{Y}_{\rho,L,N}^{\top} \mathbf{Q} \mathbf{Y}_{\rho,L,N} + \mathbf{U}_{\rho,L,N}^{\top} \mathbf{R} \mathbf{U}_{\rho,L,N}, \quad (18a) $$

$$ \mathbf{C}_{\rho} = -\mathbf{Y}_{\rho,L,N}^{\top} \mathbf{Q} \mathbf{Y}_{r}, \quad (18b) $$

The resulting predictive control algorithm is summed up in Algorithm 1.

**Remark 1.** The former results as well as Algorithm 1 are valid for noise free data only. If the available data are noisy, a regularization term
<figure><img src="image_2.png" /><figcaption>Fig. 1. Timeline showing the collected data for the Hankel matrices $Y_{0,p,N}$ and $Y_{p,L,N}$ (top) and the past and future vectors at each instant $k$, $Y_{k-\rho,\rho,1}$ and $Y_{k,L,1}$, respectively (bottom).</figcaption></figure>

should be added to the optimization problem (16) in order to mitigate the noise effect. For instance, inspired by Verheijen, Breschi, and Lazar (2023), the following problem can be suggested

$$ \arg \min_{\alpha(k)} (\|\alpha(k)\|_{\mathbf{W}}^2 + 2C_{\rho}^{\mathrm{T}}\alpha(k) + \lambda h(\alpha)), \quad (19) $$

where $\lambda \in \mathbb{R}^+$ and $h(.)$ a nonnegative function. The regularization term is often chosen as a squared two-norm $\|\alpha(k)\|_2^2$ but this choice modifies the initial problem and does not leads to the exact solution. To overcome this problem, the following alternative is introduced in Dörfler, Coulson, and Markovsky (2023)

$$ h(\alpha) = \|(\mathbf{I} - \Pi)\alpha(k)\|_2^2, \quad \text{where } \Pi = \begin{bmatrix} U_{0,\rho,N} \\ Y_{0,\rho,N} \\ U_{\rho,L,N} \end{bmatrix}^{\mathrm{T}} \begin{bmatrix} U_{0,\rho,N} \\ Y_{0,\rho,N} \\ U_{\rho,L,N} \end{bmatrix}, \quad (20) $$

but this specific choice of $h(.)$ increases the time load of the corresponding DeePC.

### Algorithm 1 Data Enabled Predictive Control

**Phase 1 – Before closing the control loop:**

1. Choose $Q, R, \rho, L, \mathcal{Y}, \mathcal{V}$, and $Y_r$ suitably

2. Collect the open loop data sets $Y_{0,\rho+L,N}, U_{0,\rho+L,N}$

3. Compute $W$ and $C_p$ from (18)

**Phase 2 – During the control loop**

4. **for** $k=1, 2, \dots$ **do**

5. Measure $y(k)$

6. Minimize the optimization function in Eq. (16) for $\alpha(k)$

7. Compute $U_{k,L,1} = U_{\rho,L,N} \alpha(k)$

8. Apply the input $u(k) = U_{k,1,1}$ to the system

9. Update $C_p$ if $Y_r$ changes

10. Update the past vectors $U_{k-\rho,\rho,1}$ and $Y_{k-\rho,\rho,1}$

11. **end for**

### 3.3. Extension of data enabled predictive control for LPV systems

LPV systems are highly beneficial for modeling parametric uncertainties, nonlinearities, or variations in operating points (Toth, 2010). In state representation, the matrices $(A, B, C, D)$ are contingent upon varying parameters. For the analysis, and more particularly for robustness concerning parametric uncertainties, it is essential to be aware of the variation intervals of these parameters. This requirement holds even when defining a robust control law that encompasses the entire parameter variation domain. Nevertheless, this approach tends to exhibit a conservative nature when applied to the real-time adaptation of control in response to parametric variations. In analysis, particularly in the study of robustness concerning parametric uncertainties, it is imperative to have knowledge of the variation intervals of the parameters. This requirement holds true even when a robust control law is established to span the entire domain of parameter variation.

This section aims at addressing the challenge of real-time awareness of parameter variations while maintaining an adaptive control law aligned with the system dynamics at each moment. To achieve this, we present an extension of the Fundamental Lemma 1 tailored for LPV systems.

The considered model structure is noted as $S = (p, A, B, C, D)$ and defined by Eq. (9). Let matrix $E$ be defined by

$$ E = \begin{bmatrix} \mathbf{A} \\ \mathbf{C} \\ \vdots \\ \mathbf{D} \end{bmatrix}. $$

Without loss of generality, it is assumed that $E$ has affine dependence, i.e.,

$$ E = E_0 + \sum_{h=1}^{n_p} p_h E_h, \quad (21) $$

where $E_0$ is the nominal term. Let us define

$$ m_h = \frac{p_h + \bar{p}_h}{2}, \quad w_h = \frac{\bar{p}_h - p_h}{2}, \quad \forall h \in \{1, \dots, n_p\}. $$

Then Eq. (21) can be written as

$$ \begin{aligned} E &= E_0 + \sum_{h=1}^{n_p} \left( \frac{p_h - m_h}{w_h} w_h E_h \right) + \sum_{h=1}^{n_p} m_h E_h, \\ &= Z_0 + Z, \end{aligned} $$

with

$$ Z = \sum_{h=1}^{n_p} \theta_h Z_h, \qquad \theta_h = \frac{p_h - m_h}{w_h}, $$

$$ Z_h = w_h E_h, \qquad Z_0 = E_0 + \sum_{h=1}^{n_p} m_h E_h. $$

Note that $\theta_h$ are normalized parameter variation such that

$$ -1 \le \theta_h \le 1, \quad \forall h \in \{1, \dots, n_p\}. $$

Consider the $\mathbb{R}^{n_p}$-vectors $\phi_k$, $k \in \mathbb{N}^*$, with components equaling either 0 or 1, associated with the binary numbers, and then define

$$ \Psi_k = \mathbb{1}_{n_p,1} - 2\phi_k, \quad \forall k \in \mathbb{N}^*. $$

Define $N_p = 2^{n_p}$ and note $\Psi = [\psi_1 \quad \dots \quad \psi_{N_p}]$. Let $G$ be defined by $G = [Z_1^T \quad \dots \quad Z_{N_p}^T]^T$. When every $p_h$ varies in its range, $\theta = [\theta_1 \quad \dots \quad \theta_{N_p}]^T$ describes an hypercube. As a consequence, $G$ describes an orthotope of matrices whose vertices, denoted by $M_i$, $i = 1, \dots, N_p$, correspond to extreme values of various $p_h$ and are given by

$$ M_i = \psi_i^T G, \quad \forall i \in \{1, N_p\}. $$

In other words, $Z$ can be written as a convex combination of vertices $M_i$, i.e.,

$$ Z = \sum_{i=1}^{N_p} \beta_i M_i, $$

with

$$ \beta_i \ge 0 \quad \forall i \in \{0, \dots, N_p\}, \quad \sum_{i=1}^{N_p} \beta_i = 1, $$

or equivalently,

$$ E = \sum_{i=1}^{N_p} \beta_i \left( E_0 + \sum_{h=1}^{n_p} m_h E_h + M_i \right) \quad (22a) $$

$$ = \sum_{i=1}^{N_p} \beta_i E^{[i]}, \quad (22b) $$

where $E$ describes an orthotope of matrices with its $N_p$ vertices

$$ E^{[i]} = \begin{bmatrix} \mathbf{A}^{[i]} \\ \mathbf{C}^{[i]} \\ \vdots \\ \mathbf{D}^{[i]} \end{bmatrix}. $$
Let us denote $p^{[i]} \in \mathbb{R}^{n_p}$ as the fixed value representing the scheduling parameter $p$ at the $i$th vertex. Each component of $p^{[i]}$ corresponds to one of the $N_p$ combinations of $\bar{p}_h$ and $p_h$, where $h$ ranges from 1 to $N_p$. $E^{[i]}$ is derived by evaluating $E$ at the $i$th vertex with $p = p^{[i]}$. Note that $E$ is entirely described by a convex combination of its vertices. The coefficients $\beta_i$ are the parameters of the convex combination and are function of $p(k)$. Indeed, it can be shown that

$$ \beta_i = \frac{1}{N_p} \left( \prod_{h=1}^{n_p} \left( \psi_{i_h} \theta_h + 1 \right) \right) \quad \forall i \in \{1, \dots, N_p\}, \qquad (23) $$

where $\psi_{i_h}$ is the $h^{\text{th}}$ component of $\psi_i$. To simplify the notation, this dependence is omitted. In the paper, it is denoted $\beta_i$ instead of $\beta_i(p(k))$.

In prior LPV control methods documented in the literature (refer, for instance, to Verhoek et al. (2021)), the trajectory of $p$ was presumed to be known, or equivalently, the coefficients $\beta_i$ were assumed to be known at each time. In contrast, this work considers the trajectory of the parameter vector $p$ as unknown. The achievement of data driven predictive control for LPV systems with unknown parameter variations involves extending the fundamental lemma to encompass LPV systems, as shown by Proposition 1.

**Proposition 1.** Consider the LPV system (9) where the parameter $p \in \mathbb{R}^{n_p}$. Then, any output-input trajectory $[\mathbf{Y}_{k,s,1}^\top \quad \mathbf{U}_{k,s,1}^\top]^\top$ of (9), for $k \in \mathbb{N}$, belongs to the column space of the matrix

$$ H_{yu} = \begin{bmatrix} \mathbf{Y}_{0,s,N}^{[1]} & \dots & \mathbf{Y}_{0,s,N}^{[N_p]} \\ \mathbf{U}_{0,s,N}^{[1]} & \dots & \mathbf{U}_{0,s,N}^{[N_p]} \end{bmatrix} \in \mathbb{R}^{s(n_u+ny) \times N_* N_p}, \qquad (24) $$

where $\mathbf{Y}_{0,s,N}^{[1]}$ and $\mathbf{U}_{0,s,N}^{[1]}$ are the Hankel matrices of the frozen system $S^{[1]} = (\mathbf{p}^{[1]}, \mathbf{E}^{[1]})$, assuming that $\text{rank}(\mathbf{U}_{0,s,N}^{[1]}) = s n_u$ for $i = 1, \dots, N_p$.

**Proof.** Using Eqs. (22), (9) can be equivalently rewritten in the polytopic form

$$ \mathbf{x}(k+1) = \sum_{i=1}^{N_p} \beta_i \mathbf{A}^{[i]} \mathbf{x}(k) + \sum_{i=1}^{N_p} \beta_i \mathbf{B}^{[i]} \mathbf{u}(k), \qquad (25a) $$

$$ \mathbf{y}(k) = \sum_{i=1}^{N_p} \beta_i \mathbf{C}^{[i]} \mathbf{x}(k) + \sum_{i=1}^{N_p} \beta_i \mathbf{D}^{[i]} \mathbf{u}(k). \qquad (25b) $$

Let us define

$$ \mathbf{x}^{[i]}(k) = \beta_i \mathbf{x}(k), \quad \mathbf{u}^{[i]}(k) = \beta_i \mathbf{u}(k), \quad \mathbf{y}^{[i]}(k) = \beta_i \mathbf{y}(k). \qquad (26) $$

As $\sum_{i=1}^{N_p} \beta_i = 1$, it is straightforward to write

$$ \mathbf{x}(k) = \sum_{i=1}^{N_p} \mathbf{x}^{[i]}(k), \quad \mathbf{u}(k) = \sum_{i=1}^{N_p} \mathbf{u}^{[i]}(k), \quad \mathbf{y}(k) = \sum_{i=1}^{N_p} \mathbf{y}^{[i]}(k). \qquad (27) $$

The vectors $\mathbf{x}^{[i]}(k)$, $\mathbf{u}^{[i]}(k)$ and $\mathbf{y}^{[i]}(k)$ correspond to the state, input and output vectors of a frozen model at instant $k \in \mathbb{N}$ and corresponds to a vertex $S^{[i]} = (\mathbf{A}^{[i]}, \mathbf{B}^{[i]}, \mathbf{C}^{[i]}, \mathbf{D}^{[i]})$ of the orthotope $S$. This frozen model can be written in the state space form

$$ \mathbf{x}^{[i]}(k+1) = \mathbf{A}^{[i]} \mathbf{x}^{[i]}(k) + \mathbf{B}^{[i]} \mathbf{u}^{[i]}(k), \qquad (28a) $$

$$ \mathbf{y}^{[i]}(k) = \mathbf{C}^{[i]} \mathbf{x}^{[i]}(k) + \mathbf{D}^{[i]} \mathbf{u}^{[i]}(k). \qquad (28b) $$

Consider now an excitation signal persistently exciting of order $s \in \mathbb{N}^*$. Then, based on Definition 4, the output-input sequence $\{\mathbf{y}^{[i]}(k), \mathbf{u}^{[i]}(k)\}_{k=0}^{N-1}$ is a trajectory of the frozen model $S^{[i]}$.

From Lemma 1, the system output-input sequence over $s$ steps $[\mathbf{Y}_{k,s,1}^{[i]\top} \quad \mathbf{U}_{k,s,1}^{[i]\top}]^\top$ belongs to the column space of the matrix $[\mathbf{Y}_{0,s,N}^{[i]\top} \quad \mathbf{U}_{0,s,N}^{[i]\top}]^\top$. In other words, each trajectory of the frozen model $S^{[i]}$ is a linear combination of the column space of the Hankel matrices $\mathbf{Y}_{0,s,N}^{[i]}$ and $\mathbf{U}_{0,s,N}^{[i]}$.

As a result of Lemma 1, for each trajectory and each instant $k$, over $s \in \mathbb{N}^*$ steps of $S^{[i]}$ for $i = 1, \dots, 2^{n_p}$, there exists a vector $\alpha^{[i]}(k) \in \mathbb{R}^N$ such that

$$ \begin{bmatrix} \mathbf{Y}_{k,s,1}^{[i]} \\ \mathbf{U}_{k,s,1}^{[i]} \end{bmatrix} = \begin{bmatrix} \mathbf{Y}_{0,s,N}^{[i]} \\ \mathbf{U}_{0,s,N}^{[i]} \end{bmatrix} \alpha^{[i]}(k). \qquad (29) $$

Substituting Eq. (29) into Eq. (27) leads to write the trajectory of the LPV system at instant $k \in \mathbb{N}$ over $s \in \mathbb{N}^*$ steps as

$$ \begin{bmatrix} \mathbf{Y}_{k,s,1} \\ \mathbf{U}_{k,s,1} \end{bmatrix} = \sum_{i=1}^{N_p} \begin{bmatrix} \mathbf{Y}_{k,s,1}^{[i]} \\ \mathbf{U}_{k,s,1}^{[i]} \end{bmatrix} \qquad (30a) $$

$$ = \sum_{i=1}^{N_p} \begin{bmatrix} \mathbf{Y}_{0,s,N}^{[i]} \\ \mathbf{U}_{0,s,N}^{[i]} \end{bmatrix} \alpha^{[i]}(k). \qquad (30b) $$

Eq. (30) can be summed up in a general form as

$$ \begin{bmatrix} \mathbf{Y}_{k,s,1} \\ \mathbf{U}_{k,s,1} \end{bmatrix} = H_{yu} \alpha(k), \qquad (31) $$

where $\alpha(k) \in \mathbb{R}^{N_p N_*}$. This proves that any trajectory of the LPV system $S$ belongs to the column space of $H_{yu}$. $\square$

**Remark 2.** The coefficients $\beta_i$ do not appear explicitly in Eq. (31). However, they are implicitly present in the vectors $\alpha^{[i]}(k)$. As in the DeePC algorithm, $\alpha(k)$ is recalculated at each time $k$. This means that at each time a new convex combination is calculated to position matrices $(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})$ in the orthotope $S$ without an explicit measurement of $p(k)$.

Using this result, Proposition 2 can be established to implement an Data Enabled Predictive Control to the LPV systems with no prior knowledge of the parameter trajectory.

**Proposition 2.** The Data Enabled Predictive Control problem for the LPV system (9) where the parameter vector $p$ is unknown is given by

$$ \arg\min_{\alpha(k)} \|\alpha(k)\|_{W}^{2} + 2C_{p}^{\top}\alpha(k), \qquad (32) $$

subject to

$$ \begin{bmatrix} \mathbf{Y}_{k-\rho,\rho,1} \\ \mathbf{U}_{k-\rho,\rho,1} \end{bmatrix} = \begin{bmatrix} \mathbf{Y}_{0,\rho,N}^{[1]} & \dots & \mathbf{Y}_{0,\rho,N}^{[N_p]} \\ \mathbf{U}_{0,\rho,N}^{[1]} & \dots & \mathbf{U}_{0,\rho,N}^{[N_p]} \end{bmatrix} \alpha(k), \qquad (33a) $$

$$ \begin{bmatrix} \mathbf{Y}_{L,\rho,N}^{[1]} & \dots & \mathbf{Y}_{L,\rho,N}^{[N_p]} \end{bmatrix} \alpha(k) = \begin{bmatrix} \mathbf{y}_r(k+L-\rho) \\ \vdots \\ \mathbf{y}_r(k+L-1) \end{bmatrix}, \qquad (33b) $$

$$ \begin{bmatrix} \mathbf{U}_{\rho,L,N}^{[1]} & \dots & \mathbf{U}_{\rho,L,N}^{[N_p]} \end{bmatrix}^{(i)} \alpha(k) \in \mathcal{V}, \qquad (33c) $$

$$ \begin{bmatrix} \mathbf{Y}_{\rho,L,N}^{[1]} & \dots & \mathbf{Y}_{\rho,L,N}^{[N_p]} \end{bmatrix}^{(i)} \alpha(k) \in \mathcal{Y}, \qquad (33d) $$

where $\mathcal{V}$ and $\mathcal{Y}$ are the polytopes governing the input and output trajectories respectively and $i \in I_{1,L}$. $\mathbf{W}$ and $\mathbf{C}_p$ are now given as

$$ \mathbf{W} = \left\| \begin{bmatrix} \mathbf{Y}_{\rho,L,N}^{[1]} & \dots & \mathbf{Y}_{\rho,L,N}^{[N_p]} \end{bmatrix} \right\|_{Q}^{2} + \left\| \begin{bmatrix} \mathbf{U}_{\rho,L,N}^{[1]} & \dots & \mathbf{U}_{\rho,L,N}^{[N_p]} \end{bmatrix} \right\|_{R}^{2}, \qquad (34a) $$

$$ \mathbf{C}_p = - \begin{bmatrix} \mathbf{Y}_{\rho,L,N}^{[1]} & \dots & \mathbf{Y}_{\rho,L,N}^{[N_p]} \end{bmatrix}^{\top} \mathbf{Q} \mathbf{Y}_r, \qquad (34b) $$

where $\mathbf{Q} \in \mathbb{R}^{n_y L}$ and $\mathbf{R} \in \mathbb{R}^{n_u L}$ are the positive definite diagonal penalizing matrices depending on the trade off between the input and the output of the LPV system.

**Proof.** The predictive control problem is given by minimizing Eq. (12) where the relation between the system output and input trajectory, $[\mathbf{Y}_{k,L,1}^\top, \mathbf{U}_{k,L,1}^\top]^\top$ is dependent on the system to be controlled. For an LPV system $S = (\mathbf{p}, \mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})$, any output-input trajectory of the system can be written as a linear combination of the columns of $H_{yu}$ as illustrated in Proposition 1. As a result, Eq. (31) can be reformulated
with $s = \rho + L$ corresponding to $\rho \in \mathbb{N}^*$ steps before the instant $k \in \mathbb{N}^*$ and $L \in \mathbb{N}^*$ steps ahead of the instant $k$ such that

$$ \begin{bmatrix} \mathbf{Y}_{k-\rho,\rho,1} \\ \mathbf{Y}_{k,L,1} \\ \mathbf{U}_{k-\rho,\rho,1} \\ \mathbf{U}_{k,L,1} \end{bmatrix} = \begin{bmatrix} \mathbf{Y}_{0,\rho,N}^{[1]} & \cdots & \mathbf{Y}_{0,\rho,N}^{[N_p]} \\ \mathbf{Y}_{\rho,L,N}^{[1]} & \cdots & \mathbf{Y}_{\rho,L,N}^{[N_p]} \\ \mathbf{U}_{0,\rho,N}^{[1]} & \cdots & \mathbf{U}_{0,\rho,N}^{[N_p]} \\ \mathbf{U}_{\rho,L,N}^{[1]} & \cdots & \mathbf{U}_{\rho,L,N}^{[N_p]} \end{bmatrix} \alpha(k), \quad (35) $$

where $\alpha(k) \in \mathbb{R}^{N_p}$. This shows that the future output and input trajectories can be written as a function of $\alpha(k)$ such that

$$ \mathbf{Y}_{k,L,1} = \begin{bmatrix} \mathbf{Y}_{\rho,L,N}^{[1]} & \cdots & \mathbf{Y}_{\rho,L,N}^{[N_p]} \end{bmatrix} \alpha(k), \quad (36a) $$

$$ \mathbf{U}_{k,L,1} = \begin{bmatrix} \mathbf{U}_{\rho,L,N}^{[1]} & \cdots & \mathbf{U}_{\rho,L,N}^{[N_p]} \end{bmatrix} \alpha(k). \quad (36b) $$

Based on Eq. (36), the minimization problem in Eq. (12) is reformulated as (32) subject to (33a). ☐

### 3.4. Integral action for DeePC

In practical applications, addressing deterministic disturbances necessitates the incorporation of an integrator for offset-free tracking. The optimal strategy for rejecting a deterministic disturbance involves integrating its model into the design model. Consequently, to eliminate a static error, an integral action must be integrated into the design model. Unfortunately, in DeePC, there is no predefined design model available a priori. To bypass this challenge, let us consider the innovation form of a state space representation of an LTI system

$$ \mathbf{x}(k+1) = \mathbf{A}\mathbf{x}(k) + \mathbf{B}\mathbf{u}(k) + \mathbf{K}\mathbf{a}(k), \quad (37a) $$

$$ \mathbf{y}(k) = \mathbf{C}\mathbf{x}(k) + \mathbf{D}\mathbf{u}(k) + \mathbf{a}(k), \quad (37b) $$

where $\mathbf{K}$ is the Kalman gain and $\mathbf{a}$ is a random walk, i.e., $\mathbf{a}(k+1) = \mathbf{a}(k) + \mathbf{e}(k)$ with $\mathbf{e}$ a zero mean white noise. It follows that

$$ \mathbf{\Delta}\mathbf{x}(k+1) = \mathbf{A}\mathbf{\Delta}\mathbf{x}(k) + \mathbf{B}\mathbf{\Delta}\mathbf{u}(k) + \mathbf{K}\mathbf{e}(k), \quad (38a) $$

$$ \mathbf{\Delta}\mathbf{y}(k) = \mathbf{C}\mathbf{\Delta}\mathbf{x}(k) + \mathbf{D}\mathbf{\Delta}\mathbf{u}(k) + \mathbf{e}(k), \quad (38b) $$

where $\mathbf{\Delta}$ is the differentiation operator, i.e., $\mathbf{\Delta}\mathbf{z}(k+1) = \mathbf{z}(k+1) - \mathbf{z}(k)$. Let us consider an augmented state

$$ \mathbf{x}_a(k) = \begin{bmatrix} \mathbf{\Delta}\mathbf{x}(k) \\ \mathbf{y}(k) \end{bmatrix}. \quad (39) $$

The corresponding augmented model is given by

$$ \mathbf{x}_a(k+1) = \mathbf{A}_a \mathbf{x}_a(k) + \mathbf{B}_a \mathbf{\Delta}\mathbf{u}(k) + \mathbf{K}_a \mathbf{e}(k), \quad (40a) $$

$$ \mathbf{y}(k) = \mathbf{C}_a \mathbf{x}_a(k), \quad (40b) $$

with

$$ \mathbf{A}_a = \begin{bmatrix} \mathbf{A} & \mathbb{0}_{n_y, n_y} \\ \mathbf{C} \mathbf{A} & \mathbb{0}_{n_y, n_y} \end{bmatrix}, \quad \mathbf{B}_a = \begin{bmatrix} \mathbf{B} \\ \mathbf{C} \mathbf{B} \end{bmatrix}, \quad (41a) $$

$$ \mathbf{C}_a = \begin{bmatrix} \mathbb{0}_{n_y, n_x} & \mathbb{0}_{n_y, n_y} \\ \mathbb{0}_{n_x, n_x} & \mathbb{0}_{n_x, n_x} \end{bmatrix}, \quad \mathbf{K}_a = \begin{bmatrix} \mathbf{K} \\ \mathbf{C} \mathbf{K} \end{bmatrix}. \quad (41b) $$

The eigenvalues of the augmented model (40) comprise both the eigenvalues of the plant model (10) and an additional $n_y$ eigenvalues located on the unit circle. This implies the presence of $n_y$ integrators incorporated into the augmented design model. Although the augmented system and the initial system share the same outputs, their inputs differ. To introduce an integral action in DeePC, it suffices to consider Hankel matrices using incremental inputs, as detailed in Lazar and Verheijen (2022). In the LPV case, this involves defining matrices at each vertex

as follows:

$$ \begin{bmatrix} \mathbf{Y}_{0,\rho,N}^{[1]} \\ \mathbf{Y}_{\rho,L,N}^{[1]} \\ \mathbf{\Delta}\mathbf{U}_{0,\rho,N}^{[1]} \\ \mathbf{\Delta}\mathbf{U}_{\rho,L,N}^{[1]} \end{bmatrix}. \quad (42) $$

The DeePC problem with embedded integrator is recast as

$$ \arg \min_{\alpha(k)} (\|\alpha(k)\|_{\mathbf{W}_I}^2 + 2\mathbf{C}_{\rho_I}^{\top} \alpha(k)), \quad (43) $$

subject to

$$ \mathbf{Y}_{k-\rho,\rho,1} = \begin{bmatrix} \mathbf{Y}_{0,\rho,N}^{[1]} & \cdots & \mathbf{Y}_{0,\rho,N}^{[N_p]} \end{bmatrix} \alpha(k), $$

$$ \mathbf{\Delta}\mathbf{U}_{k-\rho,\rho,1} = \begin{bmatrix} \mathbf{\Delta}\mathbf{U}_{0,\rho,N}^{[1]} & \cdots & \mathbf{\Delta}\mathbf{U}_{0,\rho,N}^{[N_p]} \end{bmatrix} \alpha(k), $$

$$ \begin{bmatrix} \mathbf{y}_r(k + L - \rho) \\ \vdots \\ \mathbf{y}_r(k + L - 1) \end{bmatrix} = \begin{bmatrix} \mathbf{Y}_{L,\rho,N}^{[1]} & \cdots & \mathbf{Y}_{L,\rho,N}^{[N_p]} \end{bmatrix} \alpha(k), $$

$$ (\mathbf{S}_{L,n_u} \begin{bmatrix} \mathbf{\Delta}\mathbf{U}_{0,\rho,N}^{[1]} & \cdots & \mathbf{\Delta}\mathbf{U}_{0,\rho,N}^{[N_p]} \end{bmatrix} \alpha(k) + \mathbb{1}_{L,n_u} \mathbf{u}(k-1))^{(i)} \in \mathcal{V}, $$

$$ (\begin{bmatrix} \mathbf{Y}_{0,\rho,N}^{[1]} & \cdots & \mathbf{Y}_{0,\rho,N}^{[N_p]} \end{bmatrix} \alpha(k))^{(i)} \in \mathcal{Y}, $$

where $i \in I_{1,L}$. $\mathbf{W}_I$ and $\mathbf{C}_{\rho_I}$ are given by

$$ \mathbf{W}_I = \begin{bmatrix} \mathbf{Y}_{\rho,L,N}^{[1]} & \cdots & \mathbf{Y}_{\rho,L,N}^{[N_p]} \end{bmatrix}^{\top} \mathbf{Q} \begin{bmatrix} \mathbf{Y}_{0,L,N}^{[1]} & \cdots & \mathbf{Y}_{0,L,N}^{[N_p]} \end{bmatrix}, \\ + \begin{bmatrix} \mathbf{\Delta}\mathbf{U}_{\rho,L,N}^{[1]} & \cdots & \mathbf{\Delta}\mathbf{U}_{\rho,L,N}^{[N_p]} \end{bmatrix}^{\top} \mathbf{R} \begin{bmatrix} \mathbf{\Delta}\mathbf{U}_{0,L,N}^{[1]} & \cdots & \mathbf{\Delta}\mathbf{U}_{0,L,N}^{[N_p]} \end{bmatrix}, \\ \mathbf{C}_{\rho_I} = -\mathbf{Y}_r^{\top} \mathbf{Q} \begin{bmatrix} \mathbf{Y}_{\rho,L,N}^{[1]} & \cdots & \mathbf{Y}_{\rho,L,N}^{[N_p]} \end{bmatrix}. $$

The resulting algorithm is summed up in Algorithm 3.

### Algorithm 2 LPV-Data Enabled Predictive Control

#### Phase 1 – Before closing the control loop:

1. Choose $\mathbf{Q}$, $\mathbf{R}$, $\rho$, $L$, $\mathcal{Y}$, $\mathcal{V}$, and $\mathbf{Y}_r$ suitably

2. **for** $i = 1, \dots, N_p$ **do**

3. Collect the open loop data sets $\mathbf{Y}_{0,\rho+L,N}^{[1]}$, $\mathbf{U}_{0,\rho+L,N}^{[1]}$

4. **end for**

5. Build $\mathbf{H}_{yu}$ as in Eq. (24)

6. Compute $\mathbf{W}_I$ and $\mathbf{C}_{\rho_I}$ from Eq. (18)

#### Phase 2 – During the control loop

7. **for** $k = 1, 2, \dots$ **do**

8. Measure $\mathbf{y}(k)$

9. Minimize the optimization function in Eq. (43) for $\alpha(k)$

10. Compute $\mathbf{\Delta}\mathbf{U}_{k,L,1} = \mathbf{\Delta}\mathbf{U}_{\rho,L,N} \alpha(k)$

11. Apply the input $\mathbf{u}(k) = \mathbf{\Delta}\mathbf{U}_{k,1,1} + \mathbf{u}(k-1)$ to the system

12. Update $\mathbf{C}_{\rho_I}$ if $\mathbf{Y}_r$ changes

13. Update the past vectors $\mathbf{\Delta}\mathbf{U}_{k-\rho,\rho,1}$ and $\mathbf{Y}_{k-\rho,\rho,1}$

14. **end for**

### 3.5. Computational complexity

The formulation of the LPV system in the polytopic form (25) needs $N_p = 2^{n_p}$ vertices. This implies that for trajectory definition, we require $\alpha(k) \in \mathbb{R}^{N_p}$, The exponential increase in the number of vertices $N_p$ as the number of scheduling variables $p$ increase can lead to large optimization problem and computational challenges. To
overcome this issue, one way is to consider the γ-DDPC algorithm first introduced in Breschi, Chiuso, and Formentin (2023b) and Breschi, Fabris, Formentin and Chiuso (2023) to parametrize the solution in terms of a lower dimensional parameter vector α(k). Let us define

$$Z_{ini}(k) = \begin{bmatrix} \mathbf{Y}_{k-\rho,\rho,1} \\ \Delta U_{k-\rho,\rho,1} \end{bmatrix}, \quad Z_p^{[i]} = \begin{bmatrix} \mathbf{Y}_{0,\rho,N}^{[i]} \\ \Delta U_{0,\rho,N}^{[i]} \end{bmatrix},$$

and

$$\begin{align*} Z_p &= [Z^{[1]} \quad \dots \quad Z^{[N_p]}], \\ Y^f &= [Y_{\rho,L,N}^{[1]} \quad \dots \quad Y_{\rho,L,N}^{[N_p]}], \\ \Delta U^f &= [\Delta U_{\rho,L,N}^{[1]} \quad \dots \quad \Delta U_{\rho,L,N}^{[N_p]}]. \end{align*}$$

Consider the LQ decomposition

$$\begin{bmatrix} Z_p \\ Y^f \\ \Delta U^f \end{bmatrix} = \begin{bmatrix} L_{11} & 0 & 0 \\ L_{21} & L_{22} & 0 \\ L_{31} & L_{32} & L_{33} \end{bmatrix} \begin{bmatrix} Q_1 \\ Q_2 \\ Q_3 \end{bmatrix}. \quad (44)$$

By multiplying both sides of Eq. (44) by α(k), we have

$$\begin{bmatrix} Z_{ini}(k) \\ Y_{k,L,1} \\ \Delta U_{k,L,1} \end{bmatrix} = \begin{bmatrix} L_{11} & 0 & 0 \\ L_{21} & L_{22} & 0 \\ L_{31} & L_{32} & L_{33} \end{bmatrix} \begin{bmatrix} \gamma_1(k) \\ \gamma_2(k) \\ \gamma_3(k) \end{bmatrix}, \quad (45)$$

where $\gamma_i(k) = Q_i\alpha(k)$ for $i = 1, 2, 3$. The dimensions of $\gamma_i(k)$ are independent of $N$ and $N_p$. The optimization problem becomes computationally tractable by using $\gamma_i(k)$ as optimization variables instead of $\alpha(k)$, whatever the number of vertices of the polytope. Moreover, since $Z_{ini}(k)$ is known at time $k$, $\gamma_1(k)$ can be obtained from

$$\gamma_1(k) = L_{11}^\top Z_{ini}(k). \quad (46)$$

By denoting

$$\gamma(k) = \begin{bmatrix} \gamma_2(k) \\ \gamma_3(k) \end{bmatrix},$$

problem (43) thus reduces to

$$\arg\min_{\gamma(k)} (\|\gamma(k)\|_{\mathbf{W}_{\gamma}}^2 + 2C_{p_{\gamma}}^{\top}\gamma(k)), \quad (47)$$

subject to

$$\begin{bmatrix} \mathbf{y}_r(k + L - \rho) \\ \vdots \\ \mathbf{y}_r(k + L - 1) \end{bmatrix} = [\mathbf{A}_1 \gamma(k) + \mathbf{L}_{21} \gamma_1(k)]_{(L-\rho+1:L)}, \\ (\mathbf{S}_{L,n_u}(\mathbf{A}_2 \gamma(k) + \mathbf{L}_{31} \gamma_1(k)) + \mathbf{1}_{L,n_u} \mathbf{u}(k-1))^{(i)} \in \mathcal{V}, \\ (\mathbf{A}_1 \gamma(k) + \mathbf{L}_{21} \gamma_1(k))^{(i)} \in \mathcal{Y},$$

where $i \in I_{1,L}$, $[\cdot]_{(n:m)}$ represents the rows $n$ to $m$ of $[\cdot]$, whereas

$$\mathbf{A}_1 = [L_{22} \quad 0], \quad (48)$$

$$\mathbf{A}_2 = [L_{32} \quad L_{33}], \quad (49)$$

$$\mathbf{W}_{\gamma} = \mathbf{A}_1^{\top} \mathbf{Q} \mathbf{A}_1 + \mathbf{A}_2^{\top} \mathbf{R} \mathbf{A}_2, \quad (50)$$

$$C_{p_{\gamma}} = \mathbf{A}_1^{\top} \mathbf{Q} (\mathbf{L}_{21} \gamma_1(k) - \mathbf{Y}_r) + \mathbf{A}_2^{\top} \mathbf{R} \mathbf{L}_{31} \gamma_1(k). \quad (51)$$

The resulting algorithm is given in Algorithm 3.

## 4. Numerical simulations

To demonstrate the effectiveness of this novel method, two examples are presented in this section. Firstly, the proposed control method undergoes testing on a numerical LPV system, as detailed in Section 4.1. Subsequently, the effectiveness of the proposed control strategy is evaluated using simulated data generated from a high-fidelity simulator replicating a calendering process.

### Algorithm 3 γ-LPV-Data Enabled Predictive Control

**Phase 1 – Before closing the control loop:**

1: Choose $\mathbf{Q}, \mathbf{R}, \rho, L, \mathcal{Y}, \mathcal{V}$, and $\mathbf{Y}_r$ suitably

2: **for** $i = 1, \dots, N_p$ **do**

3: Collect the open loop data sets $\mathbf{Y}_{0,\rho+L,N}^{[i]}$, $\mathbf{U}_{0,\rho+L,N}^{[i]}$

4: **end for**

5: Build $\mathbf{H}_{y_u}$ as in Eq. (24)

6: Compute the LQ decomposition as in Eq. (44)

7: Compute $\mathbf{W}_{\gamma}$ and $C_{p_{\gamma}}$ from Eq. (18)

**Phase 2 – During the control loop**

8: **for** $k = 1, 2, \dots$ **do**

9: Measure $\mathbf{y}(k)$

10: Compute $\gamma_1(k)$ from Eq. (46)

11: Minimize the optimization function in Eq. (43) for $\gamma(k)$

12: Compute $\Delta U_{k,L,1} = \mathbf{L}_{31}\gamma_1(k) + \mathbf{A}_2\gamma(k)$

13: Apply the input $\mathbf{u}(k) = \Delta U_{k,1,1} + \mathbf{u}(k-1)$ to the system

14: Update $C_{p_{\gamma}}$ and the past vectors $Z_{ini}(k)$

15: **end for**

### 4.1. Numerical example

Consider the LPV system $S = (\mathbf{p}, \mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})$ where

$$\mathbf{A} = \begin{bmatrix} 1 - 0.001 p_1(k) & 0.01 & 0.02 p_2(k) & 0 \\ 0 & 1 - 0.01025 p_1(k) & 0 & 0.033 \\ 0.01 p_2(k) & 0 & 1 - 0.0095 p_1(k) & 0 \\ 0 & 0.05 p_2(k) & 0.0375 p_1(k) & -0.053 \end{bmatrix}, \\ \mathbf{B} = \begin{bmatrix} 0.017 \\ 0.001 \\ 0.5 \\ 0.072 \end{bmatrix}, \quad \mathbf{C} = \begin{bmatrix} 1 & 0 & 0 & 0 \end{bmatrix}, \quad \mathbf{D} = 0.$$

All the possible trajectories of the scheduling parameter $\mathbf{p}(k) = \begin{bmatrix} p_1(k) \\ p_2(k) \end{bmatrix}$ evolve inside a convex polytope such that

$$\begin{bmatrix} 8 \\ 0.1 \end{bmatrix} \le \mathbf{p}(k) \le \begin{bmatrix} 16 \\ 2 \end{bmatrix}.$$

From Eq. (22), an orthotope $\mathbf{E}$ can be constructed with four vertices. Each vertex $\mathbf{E}^{[i]}$, $i = 1, \dots, 4$ defines an LTI system which is excited with an input signal $\mathbf{u}_i^{\top}$ chosen as a Pseudo Random Binary Signal (PRBS) of length $N + \rho + L$. The data length is $N = 300$ time steps and the prediction horizon $L$ is set as 40 time steps. This choice ensures the persistency of excitation of the input signal. The resulting Hankel matrices of the $N_p$ LTI systems are horizontally stacked as per Eq. (30). The output and input penalty weights are chosen as $Q = 10$ and $R = 1$. The control objective is to track a reference trajectory $\mathbf{y}_r(k)$ while the scheduling parameter $\mathbf{p}(k)$ follows a trajectory demonstrated in Fig. 2. The input $\mathbf{u}$ is bounded between $[-20, 20]$.

The control input and output results from the closed-loop implementation of the control algorithm with the LPV system over 6000 time steps are presented in Figs. 3 and 4 respectively. As evident in Fig. 4, the DeePC algorithm (43)–(41) adeptly guides the system output (depicted by the blue line) toward the reference trajectory (indicated by the red line), effectively handling variations in the scheduling parameters. For comparison, the LPV system is also controlled using a DeePC controller tuned for an LTI system, considering the nominal LTI model of the studied system. The nominal model, representing the barycenter of $\mathbf{E}$, corresponds to the frozen model defined for $\mathbf{p} = [12, 105]$. The resulting output is depicted by the black line in Fig. 4. Notably, incorporating the
<figure><img src="image_3.png" /><figcaption>Fig. 2. Trajectory of the scheduling parameter $p(k)$. The nominal value $p = [12, 1.05]$ corresponds to the starting point ($k = 1$).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Control inputs with DeePC for LPV system (blue line) and DeePC for nominal LTI system (black line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

LTI models on the vertices of the orthotope $E$ surpasses the standard solution with the nominal LTI model. In the latter case, the controller fails to guide the system output to the required trajectory. It is important to note that the algorithm does not always find a solution to the optimization problem and may yield a saturated control signal, as shown in Fig. 3 (black line). In contrast, the DeePC algorithm for LPV systems yields an unsaturated control signal, facilitating the system output to precisely track its reference trajectory, as illustrated by the blue line in Fig. 3.

## 4.2. Calendering process control

### 4.2.1. Problem formulation

Through the tire manufacturing process, a pivotal calendering stage is employed to create thin rubber layers for subsequent processing. This crucial step entails compressing the rubber material through a series of rotating cylinders. Its importance is underscored by its role in preparing

<figure><img src="image_5.png" /><figcaption>Fig. 4. Controlled outputs with DeePC for LPV system (blue line) and DeePC for nominal LTI system (black line). Red line corresponds to reference trajectory. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. A two dimensional scheme showing the calendering process with the spatial variation of the temperature at a given calendering velocity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

thin and seamless layers, ultimately shaping the overall quality of the manufactured tires. A two-dimensional scheme of the calendering rolls is illustrated in Fig. 5 depicting the spatial variation of the temperature within the rubber at a specific calendering velocity. Lower temperatures are visually represented in blue and higher temperatures are depicted in red. Within the calendering process, the viscoelastic properties of rubber can give rise to the phenomenon known as viscous heating. This heat generation is a function of the production rate, namely the rotational speed of the calendering rolls. Notably, substantial temperature fluctuations occur, particularly at the inlet, where the rubber mixture accumulates to form a rubber bed, and at the outlet, where the rubber exits the gap between the two cylindrical wheels. These specific regions experience the most significant temperature changes throughout the calendering process.

Maintaining the rubber temperature within a specific setpoint is crucial to preserve the quality of the rubber, underscoring the critical role of efficient control temperature dynamics in optimizing the tire manufacturing process. Despite the significant impact of this mechanical process on overall production, achieving full automation proves challenging due to the high nonlinearities inherent in the system. In this particular scenario, the objective is to regulate the rubber temperature at the process outlet while adhering to the constraints imposed by
<figure><img src="image_7.png" /><figcaption>Fig. 6. Scheme of studied domain and the boundary conditions used for calculations.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Diagram showing the controller implementation in closed loop with the process simulator.</figcaption></figure>

the maximum temperature the mixture can endure before breaking down. This control is implemented using a high-fidelity simulator of the calendering process. It is important to note that temperature variations at the inlet, influenced by the rubber feeding system, are not addressed in this example as they fall beyond its scope.

### 4.2.2. Solution framework

Within the scope of this example, the governing equations that model the dynamics of the rubber flow through the calendering process are formulated based on the following assumptions: (i) two dimensional flow, (ii) unsteady flow, (iii) the elastic property of the rubber is neglected thus eventually, the rubber is treated as an incompressible non Newtonian viscous fluid flow. Consequently, the equations characterizing the rubber behavior encompass the generalized Stokes equation for momentum conservation, a divergence-free velocity equation for mass conservation, and a convection-diffusion equation that incorporates a term for viscous heating for energy conservation. The constitutive law applied across all scenarios presented is a power-law coupled with an Arrhenius law, effectively accounting for temperature-dependent viscosity. (For the comprehensive set of physical equations, please refer to Appendix). To achieve precise High-Fidelity (HF) solutions for temperature profiles, we utilize an in-house finite element solver called MEF++ (developed collaboratively by Michelin and GIREF Laboratory Laval & GIREF, 2021). This simulator effectively solves the discretized physical equations that govern the calendering process with a sampling time of 1 s. The discretized domain and the associated boundary conditions utilized for the finite elements simulations are detailed in Fig. 6. As previously mentioned, the rotational velocities of the cylinders, denoted by $V_{cyl1}$ for the upper cylinder and $V_{cyl2}$ for the lower cylinder in Fig. 6, serve as the inputs of the system and they are to be adjusted to control the temperature of the rubber. Both of these cylinders in reality rotate together with a constant ratio, expressed as

$$ \frac{V_{cyl1}}{V_{cyl2}} = fr \qquad (52) $$

where $0 < fr \le 1$ represents the frictional constant between the two calendering wheels. In the scope of this example, the frictional constant is chosen unitary. As a result, the rotational speeds are the same, leading to $V_{cyl1} = V_{cyl2}$. This leads to a unified input for the system, denoted as $V_{cyl}$. Based on this, the LPV data driven controller can be implemented as shown in Fig. 7.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Random cylinder rolls velocity signal used to excite the system dynamics.</figcaption></figure>

### 4.2.3. Control implementation and results

The dynamic viscosity of the rubber, denoted as $\eta$, is a function of both the rubber temperature and the calendering velocity in the physical domain, as explained in Appendix. In order to find the bounds of this varying parameter, a multi step signal is used to excite the system dynamics. Specifically, for the considered application, the rotational speeds of the calendering wheels range between 0.01 m/s and 0.7 m/s. Accordingly, the viscosity of the rubber (expressed in Pa.s) is found to fall within the interval of [30, 600]. The established viscosity bounds result in a Reynolds number, denoted as Re, that is less than 1 aligning with the physical assumptions inherent in the equations presented in Appendix. Based on this, the calendering system can be modeled as an LPV system where the scheduling parameter is the dynamic viscosity $\eta$. All the feasible trajectories of $\eta$ evolve in the convex polytopic set, $\eta^v$ defined by $N_p = 2$ vertices such that

$$ \eta^v = \{30, 600\}. \qquad (53) $$

At each bound of the dynamic viscosity range, a constant dynamic viscosity is assumed. From a physical standpoint, this can be achieved by employing two different materials, each characterized by a constant viscosity. Subsequently, the resulting LTI systems are simulated using the signal presented in Fig. 8.

On the one hand, a prediction horizon of 15 time steps is adopted, $L = 15$. On the other hand, the number of past time steps required must be at least equal to the number of states corresponding to the system. The number of actual states is equal to the number of nodes based on the spatial discretization of the physical domain. In this example, the discretized domain in Fig. 6 contains 8387 space nodes. In order to determine the number of states of this system, the simulator is excited over 200 time steps and a snapshot of the complete system states (i.e, temperature field) is collected. A singular value decomposition of the temperature snapshot matrix shows that only 9 singular values are sufficient to capture 99.99% of the temperature variance. Consequently, the system order can be reduced to $n_x = 9$ states, (for further details about model order reduction, the reader is referred to Schilders, van der Vorst, and Rommes (2008)). Based on this, the past data length is set as 10 time steps.

In addressing the minimization problem for the controller, the weighting matrices $Q$ and $R$ are selected as 10 and 0.1 respectively. The upper and lower limits for the calendering velocity are established at 0.01 m/s and 0.7 m/s, aligning with the physical constraints of the calendering process. Our new solution performance is benchmarked against DeePC for LTI systems, considering the system behavior as linear time-invariant (LTI), and DeePC for LTI systems based on a
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Temporal variation of the output temperature (K) (top) and the calendering speed as input (m/s) (bottom) after the implementation of three different controllers for a given mixture of rubber.</figcaption></figure>

nominal system. The latter assumes a fixed viscosity value at $\eta = 450$ Pa.s, selected to ensure that the temperature variation aligns with that of the considered rubber mixture. The outcomes of applying the three controllers are illustrated in Fig. 9 for a given rubber mix. Fig. 9 illustrates the effective application of the proposed method, showcasing its ability to adjust the system input and drive the output to converge accurately toward various reference values with minimal error. Notably, the DeePC based on the nominal system successfully controls the system for the first two reference values but encounters difficulty with the third. Conversely, the controller assuming the system as a Linear Time Invariant (LTI) model fails to regulate the system when the reference is set at $T_r = 410$ K. These outcomes, stemming from the assumption of the system as an LTI model, are expected due to the inherent limitations of LTI models in capturing system nonlinearities, contrasting with the high efficiency of modeling the system using an LPV model.

In real-world applications, the viscosity of rubber material varies across different mixtures. To assess the efficacy of our proposed method, the controller undergoes testing with two distinct rubber mixtures. The viscosity variations of the two considered mixtures is shown in Fig. 10 due to the excitation of the calendering system using the calendering velocity signal depicted in Fig. 8.

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Temporal variation of the viscosity of two different rubber mixtures due to the excitation of the calendering process using the calendering speed in Fig. 8.</figcaption></figure>

It is crucial to emphasize that the temporal variation of viscosity for both materials around the output occurs within the polytope $\eta^p$, staying within the prescribed excitation limits. The outcomes of implementing the introduced control method on the two materials are illustrated in Figs. 11 and 12 accordingly. Figs. 11 and 12 depict the system output converging to the setpoint with slight oscillations around it for the two types of rubber mixtures undergoing calendering. Notably, when the required calendering rolls velocity reaches 0.275 m/s, exceeding the maximum velocity used during excitation (0.25 m/s, as shown in Fig. 8), the controller adeptly manages system dynamics variations, showcasing its robustness. This underscores the effectiveness of the proposed approach across different rubber types without prior knowledge of the corresponding viscosity temporal variation. Additionally, the controller exhibits inherent versatility, enabling its direct application to various materials without requiring an initialization procedure. As anticipated, controllers utilizing Hankel matrices from system excitation with a base material or assuming the system is an LTI system fail to regulate the system effectively, struggling to steer the output toward the predefined setpoint.

## 5. Conclusion

This paper confronts two core challenges in predictive control for LPV systems: the identification of a system model and the real-time measurement of time-varying parameters. Through a polytopic formulation of the model, a novel derivation of Willem's lemma is introduced. The proposed approach involves measuring a system trajectory at each vertex of the polytope, necessitating the freezing of parameters at all combinations of their extreme values. However, the notable advantage lies in eliminating the need to know the scheduling parameter's value at each moment for effective system control. This methodology preserves the advantages of Data Driven Predictive Control approaches, such as the absence of an identification stage, inherent consideration of multi-variable systems, an output feedback control law without the need for designing an observer, and the incorporation of equality and/or inequality constraints on inputs/outputs.

The primary drawback of this method is its dependency on the ability to generate a trajectory on each vertex with persistent excitation. This requires accessibility and tuning of the components of $p$, which may not always be feasible or practical. Another limitation lies in the requirement for $p$ to evolve within a polytope, stemming from the necessity of modeling the LPV system in polytopic form. However, the challenge of exponential computational complexity growth with the increasing number of parameters can be solved by the way of an LQ
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Figure showing the variation of the output temperature and the adjusted input after the implementation of the LPV controller using the first mixture.</figcaption></figure>

transformation conducted offline, i.e. prior to closing the control loop. The resulting computational complexity akin to the initial DeePC for LTI systems, irrespective of the scheduling parameter dimension.

The effectiveness of this approach is substantiated by results obtained from a high-fidelity, realistic process simulator modeling a rubber calendering process employed in tire manufacturing. The subsequent phase of this study involves implementing this novel method in an actual operational process and to evaluate its robustness in presence of noisy data.

### CRediT authorship contribution statement

**Taleb Bou Hamdan:** Writing – original draft. **Patrick Coirault:** Supervision, Methodology, Conceptualization. **Guillaume Mercère:** Supervision, Methodology. **Thibault Dairay:** Supervision, Project administration.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Coirault Patrick reports was provided by University of Poitiers. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. Figure showing the variation of the output temperature and the adjusted input after the implementation of the LPV controller using the second mixture.</figcaption></figure>

## Appendix. Physical equations

Based on the assumptions presented in Section 4.2.2, the conservation equations can be formulated in a simplified manner. The conservation of momentum (Batchelor, 1967) can be written in the form

$$ - \nabla \cdot (2\eta(\mathbf{u})\dot{\mathbf{u}}(\mathbf{u})) - \nabla \cdot \mathbf{p} = 0, \quad (A.1) $$

where $\nabla$ stands for the divergence operator, $\mathbf{u}$ is the two dimensional velocity vector, $\mathbf{p}$ is the pressure applied on the rubber, $\dot{\mathbf{u}}$ is the general deformation rate related to the shear rate, and $\eta$ is the dynamic viscosity of the rubber. This equation is applicable for very low Reynolds number, denoted as Re (Lee, Lee, & Yeo, 2014). The Reynolds number is a dimensionless quantity used in fluid mechanics to characterize the flow of a fluid (liquid or gas). It provides information about the relative importance of inertial forces to viscous forces in the fluid flow. The definition of the Reynolds number is:

$$ \mathrm{Re} = \frac{\rho V L}{\eta} \quad (A.2) $$

where $L$ is a characteristic length. $V$ is a reference velocity. $\rho$ is the density of the fluid and $\eta$ is the fluid dynamic viscosity. The dynamic viscosity of the rubber follows a power law (George & Qureshi, 2013)
with

$$ \eta(\mathbf{u}) = \eta(\dot{\gamma}, T) = K(\tau\dot{\gamma})^{n-1} \exp\left(\alpha\left(\frac{1}{T} - \frac{1}{T_{ref}}\right)\right), \quad (A.3) $$

where $\tau$ is a time characteristic constant, $K$ is the fluid consistency, $n$ is the exponent of the power law such that $n < 1$, $\alpha$ is the ratio of the material activation energy over the perfect gas constant, $\alpha = E_a/R$, $T$ is the temperature of the rubber and $T_{ref}$ is the reference temperature of the Arrhenius law. The mass conservation equation for an incompressible flow is given by

$$ \nabla \cdot \mathbf{u} = 0. \quad (A.4) $$

The energy conservation equation is given by (Fletcher, 1991):

$$ \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T = \frac{\lambda}{\rho C_p} \nabla^2 T + \frac{\eta(\mathbf{u})\dot{\gamma}^2(\mathbf{u})}{\rho C_p} \quad (A.5) $$

where $\lambda$ is the rubber thermal conductivity, $C_p$ is the specific heat capacity of the rubber. The additional term $\frac{\eta(\mathbf{u})\dot{\gamma}^2(\mathbf{u})}{\rho C_p}$ is the viscous heating source term. The dynamic viscosity $\eta(\mathbf{u})$ is calculated using the power law in Eq. (A.3). This induces the coupling between the conservation equations, i.e. Eq. (A.1), Eqs. (A.4) and (A.5) (thermo-mechanical coupling).

## References

Batchelor, G. (1967). *An introduction to fluid dynamic: vol. 52*, Cambridge University Press.

Berberich, J., & Allgöwer, F. (2020). A trajectory-based framework for data-driven system analysis and control. In *European control conference* (pp. 1365–1370). http://dx.doi.org/10.23919/ECC51009.2020.9143608.

Berberich, J., Köhler, J., Müller, M. A., & Allgöwer, F. (2021). On the design of terminal ingredients for data-driven MPC. *IFAC-PapersOnLine*, 54(6), 257–263. http://dx.doi.org/10.1016/j.ifacol.2021.08.554, 7th IFAC Conference on Nonlinear Model Predictive Control NMPC 2021.

Breschi, V., Chiuso, A., & Formentin, S. (2023a). Data-driven predictive control in a stochastic setting: a unified framework. *Automatica*, (152), 1–16. http://dx.doi.org/10.1016/j.automatica.2023.110961.

Breschi, V., Chiuso, A., & Formentin, S. (2023b). Data-driven predictive control in a stochastic setting: a unified framework. *Automatica*, 152, Article 110961. http://dx.doi.org/10.1016/j.automatica.2023.110961.

Breschi, V., Fabris, M., Formentin, S., & Chiuso, A. (2023). Uncertainty-aware data-driven predictive control in a stochastic setting. *IFAC-PapersOnLine*, 56(2), 10083–10088. http://dx.doi.org/10.1016/j.ifacol.2023.10.878, 22nd IFAC World Congress.

Dinkla, R., Mulders, S., van Wingerden, J.-W., & Oomen, T. (2023). Closed-loop aspects of data-enabled predictive control. In *IFAC 22nd triennial world congress, Yokohama, Japan*, 2023 (pp. 1388–1393). http://dx.doi.org/10.1016/j.ifacol.2023.10.1806.

Dörfler, F., Coulson, J., & Markovsky, I. (2023). Bridging direct and indirect data-driven control formulations via regularizations and relaxations. *IEEE Transactions on Automatic Control*, 68(2), 883–897. http://dx.doi.org/10.1109/TAC.2022.3148374.

Favoreel, W., Moor, B. D., & Gevers, M. (1999). SPC: Subspace predictive control. *IFAC Proceedings Volumes*, 32(2), 4004–4009. http://dx.doi.org/10.1016/S1474-6670(17)56683-5, 14th IFAC World Congress 1999, Beijing, China, 5-9 July.

Fletcher, C. A. J. (1991). *Fluid dynamics: The governing equations. In Computational techniques for fluid dynamics 2: specific techniques for different flow categories* (pp. 1–46). Springer Berlin, chapter 1.

George, H. F., & Qureshi, F. (2013). Newton’s law of viscosity, Newtonian and non-Newtonian fluids. In *Encyclopedia of tribology* (pp. 2416–2420). Boston, MA: Springer.

Gidon, D., Abbas, H. S., Bonzanini, A. D., Graves, D. B., Mohammadpour Velni, J., & Mesbah, A. (2021). Data-driven LPV model predictive control of a cold atmospheric plasma jet for biomaterials processing. *Control Engineering Practice*, 109, Article 104725. http://dx.doi.org/10.1016/j.conengprac.2021.104725.

Kapsalis, D., Sename, O., Milanes, V., & Molina, J. J. (2022). A reduced LPV polytopic look-ahead steering controller for autonomous vehicles. *Control Engineering Practice*, 129, Article 105360. http://dx.doi.org/10.1016/j.conengprac.2022.105360.

Laval, U., & GIREF (2021). MEF+, URL: https://fr.wikipedia.org/wiki/MEF%2B%2B.

Lazar, M., & Verheijen, P. (2022). Offset-free data-driven predictive control. In *IEEE 61st conference on decision and control* (pp. 1099–1104). Cancun, Mexico: http://dx.doi.org/10.1109/CDC51059.2022.9992453.

Lee, S. H., Lee, K.-K., & Yeo, I. (2014). Assessment of the validity of Stokes and Reynolds equations for fluid flow through a rough-walled fracture with flow imaging. *Geophysical Research Letters*, 41, http://dx.doi.org/10.1002/2014GL060481.

Li, S., Nguyen, A.-T., Guerra, T.-M., & Kruszewski, A. (2023). Reduced-order model based dynamic tracking for soft manipulators: Data-driven LPV modeling, control design and experimental results. *Control Engineering Practice*, 138, Article 105618. http://dx.doi.org/10.1016/j.conengprac.2023.105618.

Markovsky, I., & Dörfler, F. (2021). Behavioral systems theory in data-driven analysis, signal processing, and control. *Annual Reviews in Control*, 52, 42–64. http://dx.doi.org/10.1016/j.arcontrol.2021.09.005.

Morato, M. M. M., Normey-Rico, J., & Sename, O. (2020). Model predictive control design for linear parameter varying systems: A survey. *Annual Reviews in Control*, (49), 64–80. http://dx.doi.org/10.1016/j.arcontrol.2020.04.016.

Rawlings, J., Mayne, D., & Diehl, M. (2017). *Model predictive control: Theory, computation, and design*. Nob Hill Publishing, URL: https://books.google.fr/books?id=MrJctAEACAAJ.

Schilders, W., van der Vorst, H., & Rommes, J. (2008). *Model order reduction: theory, research aspects and applications*. Springer.

Toth, R. (2010). *Modeling and identification of linear parameter varying systems*. Springer.

Van Overschee, P., & De Moor, B. (1994). N4SID: subspace algorithms for the identification of combined deterministic stochastic systems. *Automatica*, 30, 75–93. http://dx.doi.org/10.1016/0005-1098(94)90230-5.

Van Overschee, P., & De Moor, B. (1996). *Subspace identification for linear systems. Theory, implementation, applications*. Kluwer Academic Publishers.

Van Waarde, H. J., De Persis, C., Camlibel, M. K., & Tesi, P. (2020). Willem's fundamental lemma for state-space systems and its extension to multiple datasets. *IEEE Control Systems Letters*, 4(3), 602–607. http://dx.doi.org/10.1109/LCSYS.2020.2986991.

Verhaegen, M., & Verdult, V. (2007). *Filtering and system identification: a least squares approach*. Cambridge University Press.

Verheijen, P., Breschi, V., & Lazar, M. (2023). Handbook of linear data-driven predictive control: Theory, implementation and design. *Annual Reviews in Control*, 56, Article 100914. http://dx.doi.org/10.1016/j.arcontrol.2023.100914.

Verhoek, C., Abbas, H., Toth, R., & Haesaert, S. (2021). Data-driven predictive control for linear parameter varying systems. *IFAC PapersOnline*, 101–108. http://dx.doi.org/10.1016/j.ifacol.2021.08.588.

Willems, J. (2007). *System theoretic models for the analysis of physical systems. Ricerche di Automatica*, 10, 287–301.

Willems, J., Rapisard, P., Markovsky, I., & De Moor, B. (2005). A note on persistency of excitation. *Systems & Control Letters*, 54, 325–329. http://dx.doi.org/10.1016/j.sysconle.2004.09.003.