Robust data-based predictive control of systems with parametric uncertainties: Paving the way for cooperative learning

Journal of Process Control 132 (2023) 103109 Available online 20 October 2023 0959-1524/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Contents lists available at ScienceDirect Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont Robust data-based predictive control of systems with parametric uncertainties: Paving the way for cooperative learning Eva Masero a,b,∗, José M. Maestrea, José R. Salvadorc, Daniel R. Ramireza, Quanyan Zhu d aDepartment of Ingeniería de Sistemas y Automática, Universidad de Sevilla, Seville 41092, Spain bDepartment of Electronics, Information and Bioengineering, Politecnico di Milano, Milan 20133, Italy cDepartment of Ingeniería, Universidad Loyola Andalucía, Dos Hermanas, Seville, 41704, Spain dDepartment of Electrical and Computer Engineering, New York University, Brooklyn, NY 11201, United States ARTICLE INFO Keywords: Predictive control Data-driven control Tube-based control Robustness Cooperative learning ABSTRACT This article combines data and tube-based predictive control to deal with systems with bounded parametric uncertainty. This integration generates robustly feasible control sequences that can also be exploited in cooperative scenarios where controllers learn from each other’s data. In particular, the approach is based on a database that contains information from previous executions of the same and other controllers handling similar systems. By the combination of feasible histories plus an auxiliary control law that deals with bounded uncertainties, which only needs to be stabilizing for at least one of the system realizations within the uncertainty set, this scheme provides a finite-horizon predictive controller that guarantees exponential stability and robust constraint satisfaction. The validity and benefits of the proposed scheme are shown in case studies with linear and non-linear dynamics. 1. Introduction Multiple applications of predictive model control (MPC) in industry, e.g., [1–3], have stimulated renewed interest in approaches combining data-based and learning methods with predictive control schemes [4]. For example, those based on neural networks [5], adaptive models [6], scenarios [7], and the extension of terminal regions [8], to name a few examples. In the recent literature, we can find MPC schemes for repetitive tasks that learn from previous executions to improve its performance. For instance, in [9], the theoretical properties are derived from the use of safe sets, which can be defined as regions of the state space where there exists a control law that guarantees constraint satisfaction for all successive time steps. This topic is also addressed in [10,11], where methods to expand regions of safe states are proposed. Learning is also used to infer the plant model, as in [12], where experimental data on system inputs and outputs are used to feed a non-parametric machine learning method. Multi-step predictors learned from the data for robust MPC are used in [13]. Another strategy is to use a model-free approach such as [14], where the control signal is calculated using combinations of closed-loop trajectories stored from the past in such a way that the combined performance cost is minimized. This approach is related to direct weight optimization (DWO) techniques [15] and also to kriging methods [16, 17]. In general, both techniques rely on expressing a given datum as a ∗Corresponding author at: Department of Ingeniería de Sistemas y Automática, Universidad de Sevilla, Seville 41092, Spain. E-mail address: [email protected] (E. Masero). linear combination of data stored in the database, while some measurement of the combined weights is minimized. The DWO approach has been used in the context of identification in [18], and for data-based control purposes in [19]. On the other hand, kriging methods have been mainly used as an interpolation technique, first in geostatistics [20], but also in some control-related applications, e.g., [21] uses a basic kriged black box prediction model to design a non-linear MPC controller, and [22] introduces a new type of state space modeling technique based on kriging weights. The strategy presented in [14] was extended in [23] to obtain offset-free control under some limiting conditions. However, these strategies are highly dependent on the quality of the data and lack the robustness to track setpoints not included in the database. To address this challenge, several robust data-driven MPC techniques have emerged. For example, for the control of Type 1 diabetes, [24] combines data-driven learning with min–max MPC to obtain patterns of disturbance sets from historical data and calculate insulin infusions in the worst-case disturbances. Moreover, [25] proposes a robust data-driven MPC method that uses Hankel matrices to predict the system behavior by its input–output trajectories and includes a slack variable to relax the constraints and account for noise, i.e., the differences between the online measurements and the data used for prediction. This method is extended in [26] using constraint tightening to enforce closed-loop output constraint satisfaction. Another https://doi.org/10.1016/j.jprocont.2023.103109 Received 21 November 2022; Received in revised form 21 September 2023; Accepted 10 October 2023 Journal of Process Control 132 (2023) 103109 2 E. Masero et al. way of handling disturbances is to rely on stochastic approaches. For instance, [27] proposes a stochastic MPC based on data driven for linear systems subject to unknown additive disturbances. It is also considered a distributionally robust optimization with a Wasserstein ambiguity set based on a dataset consisting of disturbance trajectories. Then, a scenario-based tube MPC is employed to reduce the online complexity of the classical scenario-based stochastic MPC. In this article, we deal with the reutilization of past historical information following the method presented in [14], where data-based predictive control is proposed to manage linear systems through a convex combination of previously applied trajectories assuming that successive executions of the system remain constant and are disturbance-free. The key contribution of our work is to extend [14] to overcome its lack of robustness regarding constraint satisfaction, e.g., due to process noise or parametric uncertainties. To this end, we provide a set of conditions that must be met by the trajectories employed by the databased controller to ensure robust constraint satisfaction. Here we apply a strategy analogous to tube-based MPC [28], but it is a data-based predictive controller that provides us with the nominal disturbancefree trajectory for the current system. Therefore, the control problem is twofold: (i) to calculate the nominal control inputs by a convex combination of previously trajectories collected in the database, and (ii) to apply a stabilizing feedback controller to reject the differences between the real and nominal dynamics. Regarding the former, it is employed a direct weight optimization based on [14], an idea aligned with the rationale of kriging methods [21,22] and recent data-driven MPC approaches [26,29], which rely on the fundamental lemma by Willems et al. [29, Lemma 2]. The lemma establishes that a discretetime LTI system can be fully characterized by its trajectories as long as the input is persistently exciting [30]. The main differences with regard to the previously mentioned articles are that the newly proposed method guarantees recursive feasibility and robust exponential stability, and can even manage one system with trajectories generated from different instances of the system being controlled (e.g., a set of controllers sharing a cloud-based database), leading to a cooperative learning approach. To learn from each other’s trajectories, controllers need trajectories generated in a wide range of conditions and system realizations so that they are representative of the behavior of the system; also, trajectories should be stored in a format that is easily accessible and can be processed efficiently by the cooperative learning algorithm. The idea of using data from one or more systems to design a stabilizing controller for another system is also addressed in [31] in a data-based manner for the particular case of transferring a stabilizing control policy from a source system to a target system. The commonalities of our method and that of [31] are a prior knowledge of the maximum distance between two linear systems, significant data (with noise) from the source system, and the design of a stabilizing controller based on linear matrix inequalities (LMI) for the target system. However, our twist is the integration of data-driven control and predictive control to inherit the optimal and robustness properties of tube-based MPC. Our method also provides robust exponential stability and robust constraint satisfaction, which is a major issue in reusing past historical data. Another approach related to our proposal is the interpolationbased MPC [32], which relies on the convex combination of feedback control laws for a set of uncertain systems. Although this idea has some connections with the current work, note that our trajectories may have been generated by different types of controllers. Furthermore, robust constraint satisfaction is attained in [32] leveraging the wellknown approach of Kothare et al., which is known to be conservative. The referred article also points out some significant differences with standard predictive control approaches. To sum up, the novelty of our work is the combination of a datadriven MPC with a tube-based approach to address systems with polytopic uncertainty. This integration presents a novel contribution that can be leveraged to have multiple controllers learning in a cooperative fashion. It also allows for the design of controllers that are both data-driven and robust, filling a gap in the existing literature. Finally, a preliminary version of the proposed strategy was accepted for presentation [33]. There are significant differences between the current article and the conference version: (i) the problem formulation is now less conservative, with fewer constraints required; (ii) the optimization problem has now more degrees of freedom, improving performance; (iii) the current article does include theoretical analysis of the recursive feasibility and stability; (vi) new simulations and an additional example are also included in the current version. Index of contents: Section 2presents the problem setting and the proposed control strategy. Sections 3and 4detail the design and implementation of the nominal and auxiliary controllers, respectively. Section 5analyzes the recursive feasibility and stability guarantees. Sections 6and 7present the case studies. Section 8ends the paper with concluding remarks. Notation: Sets Rand Nrefer to real numbers and integers, respectively. The distance between vectors 𝑥, 𝑦 ∈R𝑛is d(𝑥, 𝑦) = ‖𝑥−𝑦‖= √(𝑥1−𝑦1)2+⋯+ (𝑥𝑛−𝑦𝑛)2, while the square of the distance is denoted as ‖𝑥−𝑦‖2. The image of a set ⊆R𝑛under a linear mapping 𝐴∶R𝑛↦R𝑚is given by 𝐴 ≜ {𝐴𝑥 ∶𝑥∈}. For sets ,⊆R𝑛, the Minkowski sum is ⊕ ≜ {𝑥+𝑦∶𝑥∈, 𝑦 ∈}, and the Pontryagin difference is ⊖ ≜ {𝑧∈R𝑛∶⊕{𝑧}⊆}for ⊆. The matrix 𝐴 ≻ 0(𝐴⪰ 0) is positive definite (semi-definite). The identity matrix is denoted as I. The set 𝛺∈R𝑛is a robust positively invariant (RPI) set for system 𝑥+=𝑓(𝑥, 𝑤)with constraints and if ∀𝑥∈𝛺 ⊂ and ∀𝑤∈, its evolution satisfies 𝑥+∈𝛺. 2. Problem formulation We consider that the system dynamics are described by the following discrete-time linear time-invariant (LTI) model with an unknown vector of parameters 𝜃∈𝛩 ⊆ R𝑛𝜃: 𝑥+=𝐴(𝜃)𝑥+𝐵(𝜃)𝑢, (1) where 𝑥, 𝑥+∈R𝑛𝑥and 𝑢∈R𝑛𝑢are, respectively, the state, the successor state, and the input of the system. The matrices 𝐴(𝜃) ∈ R𝑛𝑥×𝑛𝑥and 𝐵(𝜃) ∈ R𝑛𝑥×𝑛𝑢are the state-transition and input-to-state matrices, which depend on the realization of parametric uncertainty 𝜃. Assumption 1. All systems (𝐴(𝜃), 𝐵(𝜃)) are controllable. Assumption 2. The differences between the system realizations are bounded by polytopic sets 𝛥𝐴𝜃and 𝛥𝐵𝜃that account for possible parameter variations. Two different realizations of 𝜃, say 𝜃𝑖, 𝜃𝑗∈𝛩, lead to two distinct system realizations (𝐴𝑖, 𝐵𝑖)and (𝐴𝑗, 𝐵𝑗), with 𝐴𝑖−𝐴𝑗∈𝛥𝐴𝜃, 𝐵𝑖−𝐵𝑗∈𝛥𝐵𝜃.(2) The parametric uncertainty of the system can be reinterpreted as an unknown disturbance 𝑤. In particular, the dynamic evolution of the system realizations 𝑖and 𝑗can be related as follows by (2) 𝑥+=𝐴𝑖𝑥+𝐵𝑖𝑢 ∈ (𝐴𝑗+𝛥𝐴𝜃)+ (𝐵𝑗+𝛥𝐵𝜃) =𝐴𝑗𝑥+𝐵𝑗𝑢+𝑤, (3) where the state and inputs are subject to constraints 𝑥∈, 𝑢 ∈,(4) with the disturbance term 𝑤bounded by 𝑤∈ ≜ 𝛥𝐴𝜃⊕ 𝛥𝐵𝜃.(5) Note that, in case of need, can be enlarged to deal with additional sources of uncertainty. Assumption 3. Sets ⊆R𝑛𝑥,⊆R𝑛𝑢, and ⊆R𝑛𝑤are polytopic convex sets that contain the origin in their interiors. Journal of Process Control 132 (2023) 103109 3 E. Masero et al. Table 1 Database structure. Traj. ID 𝑘 𝑥𝑡(𝑘)𝑢𝑡(𝑘) 1 1 𝑥𝑡(1) 𝑢𝑡(1) ⋮ ⋮ ⋮ ⋮ 1𝑘1 end 𝑥𝑡(𝑘1 end)𝑢𝑡(𝑘1 end) 2 1 𝑥𝑡(1) 𝑢𝑡(1) ⋮ ⋮ ⋮ ⋮ 2𝑘2 end 𝑥𝑡(𝑘2 end)𝑢𝑡(𝑘2 end) ⋮ ⋮ ⋮ ⋮ 2.1. Database structure and control goal Any controller has access to a database that contains a set of 𝑇trajectories generated by different instances of system (1). Each trajectory is identified using indices within set  ≜ {1,…, 𝑇 }and is assumed to be feasible for the system realization that generated it, but does not need to be feasible for other system realizations. In particular, a trajectory 𝑡∈is composed of a sequence of samples (𝑥𝑡(𝑘), 𝑢𝑡(𝑘)), for 𝑘∈𝑡≜{1,…, 𝑘𝑡 end}, that satisfy (1) and (4) for the particular system realization. See Table 1 for an example regarding the structure of the considered database. Remark 1. Although we consider the general case of the database being cooperatively updated by multiple realizations of system (1), it is worth noticing that there is no difference if the database is generated by collecting data from a single plant that uses different parametric configurations. Assumption 4. At some time instant 𝑘𝑡 end, the sequences of all the trajectories 𝑡∈finish at a bounded distance from the origin smaller than a threshold 𝜀∈R+. Remark 2. The threshold 𝜀defines a ball around the origin, and must be small enough to satisfy basic conditions such as constraints. It can be empirically chosen as a function of the amplitude of noise and random disturbances that affect the tracking errors. The larger 𝜀, the greater the uncertainty of the system’s state relative to the origin. The control objective is to regulate the system to the origin while minimizing a global infinite-horizon cost: 𝐽∞= ∞ ∑ 𝑛=0 𝓁(𝑥(𝑛), 𝑢(𝑛)),(6) where 𝓁(⋅,⋅)is the stage cost, which is defined by positive definite weighting matrices 𝑄∈R𝑛𝑥×𝑛𝑥and 𝑅∈R𝑛𝑢×𝑛𝑢: 𝓁(𝑥(𝑛), 𝑢(𝑛))=𝑥(𝑛+ 1)⊤𝑄 𝑥(𝑛+ 1) + 𝑢(𝑛)⊤𝑅 𝑢(𝑛).(7) For the sake of implementability, it is desirable to consider a finite horizon 𝑁, so the cost function (6) is replaced by: 𝐽𝑁= 𝑁−1 ∑ 𝑛=0 𝓁(𝑥(𝑛), 𝑢(𝑛))+𝑓(𝑥(𝑁)),(8) where 𝑓(⋅) = 𝑥(𝑁)⊤𝑃 𝑥(𝑁)with 𝑃 ≻ 0is a terminal cost function. In this regard, let us define the 𝑁-length partial sequences of 𝐱and 𝐮, which are obtained by taking 𝑁consecutive states and inputs from trajectory 𝑡∈: 𝑝≜{𝐱= {𝑥𝑡(𝑘), 𝑥𝑡(𝑘+ 1),…, 𝑥𝑡(𝑘+𝑁)} 𝐮= {𝑢𝑡(𝑘), 𝑢𝑡(𝑘+ 1),…, 𝑢𝑡(𝑘+𝑁− 1)} , and the performance of this sequence measured by (8) can be denoted as 𝐽𝑁(𝐱,𝐮). 2.2. Dual control law We control the system realization (𝐴𝑖, 𝐵𝑖)using a data-based controller that combines the trajectories contained in . A robust feedback controller is also employed to deal with disturbance differences 𝑤due to uncertain system realizations. Therefore, the control law becomes 𝑢(𝑥) = 𝑣(𝑥) + 𝑣𝑒(𝑥, 𝑧),(9) where 𝑣(𝑥)corresponds to the first element of a sequence of control actions generated by the data-based law, and 𝑣𝑒(𝑥, 𝑧)rejects the differences between the system state 𝑥and the nominal state 𝑧, which is obtained from the database and follows disturbance-free dynamics: 𝑧+=𝐴𝑗𝑧+𝐵𝑗𝑣(recall (3)). A block diagram of the proposed dual control law is shown in Fig. 1. Sections 3and 4detail, respectively, how 𝑣(𝑥)and 𝑣𝑒(𝑥, 𝑧)are obtained. 3. Data-based control law 𝒗(𝒙) Let us consider a particular system realization (𝐴𝑖, 𝐵𝑖). In the most restrictive case, there is no previous information in the database about previous executions of this system, i.e., we need to rely on information from other system realizations. 3.1. Trajectory robustness check Since the sequences of any other system realization (𝐴𝑗, 𝐵𝑗), with 𝑗≠𝑖might not be feasible for (𝐴𝑖, 𝐵𝑖), we need to guarantee their robust feasibility. Therefore, any candidate trajectory {𝑥𝑡(𝑘), 𝑢𝑡(𝑘)}, with 𝑘∈𝑡= {1,…, 𝑘𝑡 end}, must have enough margin respect to the constraints to accommodate uncertainties given by the closed-loop dynamics of the errors due to system discrepancies. One way to do this is to check the following conditions: 𝑥𝑡(𝑘) ∈ ⊖ 𝛺𝑥, 𝑢𝑡(𝑘) ∈ ⊖ 𝛺𝑢,(10) for 𝑘∈𝑡, where 𝛺𝑥and 𝛺𝑢are given sets that define the desired margin with respect to the state and input constraints, respectively. Details on the calculation of 𝛺𝑥and 𝛺𝑢are given in Section 4. By checking (10), we can find the set of robustly feasible trajectories in the database.1Let ⊆contain the IDs of these feasible trajectories. Assumption 5. The set of robustly feasible trajectories is assumed to be nonempty. 3.2. Optimization problem Conventional tube-based MPC employs a nominal model with tighter constraints. If we considered a database that only contains information about the nominal system, we could define, for nominal dynamics 𝑧+=𝐴𝑗𝑧+𝐵𝑗𝑣, an optimal control problem P𝑁(𝑧)as [14]: 𝐽o 𝑁(𝑧) = min 𝛬∈R𝑀∑ 𝑚∈  𝜆𝑚𝐽𝑁(𝐱𝑚,𝐮𝑚), s.t.∑ 𝑚∈  𝜆𝑚𝑥𝑚(0) = 𝑧, ∑ 𝑚∈  𝜆𝑚= 1, 𝜆𝑚≥0,∀𝑚∈ , (11) where 𝑧is the current state measurement, and 𝛬∈R𝑀is the set of weights {𝜆1,…, 𝜆𝑀}for the convex combination of 𝑀partial control 1Feasible sequences can be extended with zeros so that their length becomes greater or at least equal to 𝑁,i.e.,𝑘𝑡 end ≥𝑁. Journal of Process Control 132 (2023) 103109 4 E. Masero et al. Fig. 1. Block diagram of the combined data-based and feedback controller. sequences listed in set  , which contains the indices of 𝑀≤|𝑡|⋅|| partial trajectories. Here, however, we consider a database that lacks previous information on the uncertain system 𝑥+=𝐴𝑗𝑥+𝐵𝑗𝑢+𝑤being controlled. Therefore, the previous problem P𝑁(𝑧)is transformed into P∗ 𝑁(𝑥): 𝐽∗ 𝑁(𝑥) = min 𝑧, 𝛬 ∑ 𝑚∈  𝜆𝑚𝐽𝑁(𝐱𝑚,𝐮𝑚),(12a) s.t.∑ 𝑚∈  𝜆𝑚𝑥𝑚(0) = 𝑧, (12b) ∑ 𝑚∈  𝜆𝑚= 1,(12c) 𝑥∈𝑧⊕𝛺𝑥,(12d) 𝑧∈⊖ 𝛺𝑥,(12e) 𝜆𝑚≥0,∀𝑚∈ ,(12f) where 𝑥is the current state measurement, 𝑧∈R𝑛𝑥is the associated nominal state, which is now a decision variable, and 𝛬∈R𝑀is the set of weights for the combination of 𝑀partial control sequences contained in the set of robustly feasible trajectories  , which gathers the indices of 𝑀≤|𝑡|⋅||partial trajectories. Here, we bestow problem P∗ 𝑁(𝑥)on an additional degree of freedom by optimizing the nominal state, similar to [34,35]. Problem P∗ 𝑁(𝑥)can be solved very efficiently, although in the case of stringent timing requirements and large datasets, it may be necessary to perform parallel searches in the database using GPU computing or specific hardware, e.g., FPGA. Note that the feasible nominal sets of any given instance of Problem (12) are, respectively, a set of states and inputs that are obtained, respectively, by convex combinations of database state and input trajectories using coefficients of 𝛬that satisfy the constraints. In particular, let 𝑧∗and 𝛬∗= {𝜆∗ 1, 𝜆∗ 2,…, 𝜆∗ 𝑀}be the minimizer of (12), provided that it exists,2then the nominal input partial trajectory for (3) is given by: 𝐯∗=∑ 𝑚∈  𝜆∗ 𝑚𝐮𝑚,(13) with the first element of 𝐯∗= {𝑣∗(0),…, 𝑣∗(𝑁−1)} defining the nominal control law: 𝑣(𝑥) = 𝑣∗(0),(14) and the optimized current nominal state 𝑧∗(0) is defined as: 𝑧∗(0) = ∑ 𝑚∈  𝜆∗ 𝑚𝑥𝑚(0). At the next time instant, it is produced a new combination of partial sequences based on the most recent information available in the database following a receding-horizon strategy. 2The existence of the minimizer requires the database to be large enough and have data representative of the operating conditions of the system. This issue can also be mitigated by increasing the number of samples considered for the combination (i.e.,𝑀) and using soft constraints to solve Problem P∗ 𝑁(𝑥). 3.3. Control algorithm We summarize the robust predictive controller proposed based on historical data in Algorithm 1, which is executed at each time instant. Note that, as a selection criterion to limit the computational burden, we consider that only the samples listed in  can be combined to compute the control inputs. Let d ∶×→Rbe a function that computes the distance between a system’s state 𝑥𝑖and a trajectory’s partial state 𝑥𝑡. The set  contains the indices of the 𝑀feasible partial trajectories with the smallest d(𝑥𝑖, 𝑥𝑡(𝑘))for all 𝑘∈𝑡≜{1,…, 𝑘𝑡 end}and ∀𝑡∈, and | |=𝑀with 𝑀≤∑𝑡∈|𝑡|. Algorithm 1 Robust data-based predictive controller Inputs: 𝑥𝑖,,,𝑄, 𝑅, 𝑃 ,𝑁,,𝑀,𝐾, 𝛺.Output: 𝑢. 1: For all 𝑡∈, compute the distance d(𝑥𝑖, 𝑥𝑡(𝑘))for 𝑘∈𝑡. 2: Gather in  the indices of the 𝑀feasible partial sequences (x𝑚,u𝑚) with the lowest distance to 𝑥𝑖. 3: Evaluate the cost 𝐽𝑁(x𝑚,u𝑚)of each candidate 𝑚∈ . 4: Solve (12) to optimize the weight vector 𝛬∗for the convex combination of the 𝑀control sequences, and the nominal state 𝑧∗. 5: Obtain v∗by (13), and set the first element as 𝑣=𝑣∗(0). 6: Calculate the auxiliary control input term 𝑣𝑒that penalizes the deviation between 𝑥𝑖and 𝑧∗(see Section 4for more details). 7: Obtain (9) and apply 𝑢to the system. 8: Wait for the next sampling time and repeat from Step 1. 4. Auxiliary control law 𝒗𝒆(𝒙, 𝒛) The tube-based strategy requires the auxiliary control law to keep the evolution of the error confined within set 𝛺𝑥using control actions in set 𝛺𝑢despite the role played by parametric uncertainties, which lie within set . Typically, the design of sets 𝛺𝑥and 𝛺𝑢is achieved using a feedback controller 𝐾and a robust positively invariant (RPI) set, which, together with system (3), guarantee (𝐴𝑗+𝐵𝑗𝐾)𝛺𝑥⊕⊆ 𝛺𝑥,(15) where 𝛺𝑥⊂and 𝛺𝑢≜𝐾 𝛺𝑥⊂. That is, considering 𝑣𝑒(𝑥, 𝑧) = 𝐾(𝑥−𝑧)as the auxiliary control law, the system dynamics (𝐴𝑖, 𝐵𝑖)from the viewpoint of 𝑗th system realization (3) and the double control law (9) are 𝑥+=𝐴𝑗𝑥+𝐵𝑗(𝑣+𝐾 𝑒) + 𝑤, (16) where 𝑒= (𝑥−𝑧) ∈ 𝛺𝑥captures the difference between the real and nominal values of the state, and 𝑤∈accounts for the model mismatch between systems 𝑗and 𝑖. As a consequence, the trajectory of system realization 𝑖will lie inside a tube around the nominal trajectory generated by the data-based controller. The system will be recursively feasible if there is a feasible solution in the database for its trajectory toward the origin, thanks to the combination of the database sequences. Finally, the system state will lie within a set 𝛺𝑥around the origin. The cornerstone of this approach is feedback gain 𝐾. The methods for designing 𝐾vary based on the availability and extent of prior information about the system. In the following subsections, we will discuss several strategies that can be effectively employed in this context. Journal of Process Control 132 (2023) 103109 5 E. Masero et al. 4.1. Redundancy of feedback gain The role of 𝐾is to guarantee that the closed-loop dynamics (𝐴𝑗+ 𝐵𝑗𝐾)are stable. However, some systems are inherently stable (e.g., passive systems) and therefore, Eq. (15) holds even if no feedback gain is applied, i.e.,𝐾= 0. This means that robustness can be attained if the set of system realizations (𝐴(𝜃), 𝐵(𝜃))is known to be stable. That being said, even in this case, feedback is often used in practice to improve the system’s performance and robustness features. 4.2. Empirical approach For some systems, it may be possible to find a feedback controller 𝐾 for all system realizations empirically. In such a case, 𝛺𝑥and 𝛺𝑢could be designed as balls of ratio 𝑟𝑥, 𝑟𝑢>0, respectively, i.e.,𝛺𝑥=B𝑥≜{𝑥∈ ∶dist(𝑥, 0) ≤𝑟𝑥}and 𝛺𝑢=B𝑢≜{𝑢∈∶dist(𝑢, 0) ≤𝑟𝑢}. Here, 𝑟𝑥and 𝑟𝑢become tuning parameters that might be adjusted considering the database trajectories to maintain system constraints in a robust manner. 4.3. Data-based design The design of the auxiliary control law and the invariant set can also be derived from data. Nevertheless, typical approaches are not modelfree because they often involve the implicit or explicit identification of a model based on the available data, which is then used to compute the invariant set and the control law. For example, in [36] a semidefinite programming problem (SDP) is formulated to simultaneously find the uncertain model, the feedback gain, and the corresponding invariant sets. A concurrent approach is also proposed in [37,38] to compute an approximation of a minimal robust control invariant set (mRCI) from experimental data via linear programming (LP) and select a suitable model given a model structure. 4.4. Model-based design There are multiple well-known model-based design methods to generate a feedback gain that yields stable closed-loop dynamics for a particular realization, e.g., pole placement. In this regard, note that this is all that is required by the proposed method. Likewise, given the robust nature of the problem, it may be desirable to employ a robust feedback controller 𝐾valid for all system realizations, assuming that information regarding the vertices of the uncertainty set is available. This can be performed using linear matrix inequalities (LMIs) [39] leveraging Assumption 2, which makes the disturbance set of system realizations (𝐴(𝜃), 𝐵(𝜃))polytopic. In particular, let set 𝛶denote the vertices of this set, so that the dynamics of the extreme realizations can be expressed as: 𝑥+=𝐴𝜈𝑥+𝐵𝜈𝑢, 𝜈 ∈𝛶 . (17) Lemma 1 (Optimal Feedback Controller [40]).Given system (𝐴𝜈, 𝐵𝜈)and stage cost matrices 𝑄and 𝑅. If there are matrices 𝐻=𝐻⊤∈R𝑛𝑥×𝑛𝑥and 𝑌∈R𝑛𝑢×𝑛𝑥, such that the following constraints are satisfied ⎡⎢⎢⎢⎢⎣ 𝐻 𝐻𝐴⊤ 𝜈+𝑌⊤𝐵⊤ 𝜈𝐻𝑄1∕2 𝑌⊤𝑅1∕2 𝐴𝜈𝐻+𝐵𝜈𝑌 𝐻 0 0 𝑄1∕2𝐻0 I 0 𝑅1∕2𝑌0 0 I ⎤⎥⎥⎥⎥⎦ ⪰ 0,∀𝜈∈𝛶 , (18) then there exist an optimal feedback control 𝐾=𝑌 𝐻−1, maximizing the trace of 𝐻, that stabilizes the closed-loop system, and a Lyapunov function 𝑓(𝑥)=𝑥⊤𝑃 𝑥, with 𝑃=𝐻−1. See [40, Appendix A] for the proof of Lemma 1. Given the feedback gain 𝐾and the bounded disturbances leveraging Assumption 2, which makes the disturbance set of system realizations (𝐴(𝜃), 𝐵(𝜃))polytopic, we seek to find now an RPI set (𝛺𝑥) as small as possible to reduce conservatism. For polytopic disturbances, methods such as the linear programming method proposed in [41] or the outer approximation presented in [42] can be employed to compute a common minimal RPI set for the system realizations. Similarly, one can find an RPI set for each vertex and then use the union of the RPI sets in the checking steps (10). 5. Recursive feasibility and stability Regardless of the method employed to find feedback gain 𝐾and sets 𝛺𝑥and 𝛺𝑢, it is possible to prove recursive feasibility and stability. Assumption 6. There exists, at least, a feedback controller 𝐾that can stabilize system (3) and a robust positively invariant (RPI) set 𝛺𝑥for constraint sets (,,) under the control law 𝑢=𝐾𝑥 that satisfies (15). Without lost of generality, we here consider that the database trajectories end in the origin although, as state in Assumption 4, it is not necessary. The only difference is the nominal state will be confined into a ball centered at the origin, therefore the uncertain system will be ultimately bounded in that ball plus the set 𝛺𝑥. The robust feasibility and stability properties of the proposed algorithm are inherited from tube-based MPC controllers [28], which impose the following constraints in the computation of the nominal trajectory: (i) the linear model, (ii) the set of tighter constraints (10) and (12e), and (iii) a terminal region to guarantee stability and recursive feasibility. Here, the data-based controller provides a nominal trajectory, which accounts for condition (i). The set of Eqs. (10) guarantees that only robustly feasible sequences are considered, and (12e) ensures robustly feasible nominal states (condition ii). Likewise, there is no need for a terminal region because the generated trajectories reach the origin, so that recursive feasibility is guaranteed (condition iii). Finally, the auxiliary control law is identical to that of tube-based MPC and deals with the disturbances arising from the use of a nominal trajectory that belongs to a combination of trajectories from possibly different realizations of system (1). 5.1. Recursive feasibility For the feasibility analysis of the closed-loop system, let us first define the domain of the optimization problem P∗ 𝑁(𝑥). The feasible input set of problem (12), denoted as V𝑁(𝑧), is a set of nominal input trajectories {𝑣∗(0), 𝑣∗(1),…, 𝑣∗(𝑁− 1)}, which are obtained from a convex combination of the database input trajectories using coefficients of 𝛬∗(recall (13)). Likewise, the feasible state set Z𝑁= {𝑧∈⊖ 𝛺𝑥∶V𝑁(𝑧)≠∅} is composed of nominal states that are obtained from a convex combination of the database state trajectories using coefficients of 𝛬∗. The feasible state set for the real dynamics is X𝑁=Z𝑁⊕ 𝛺𝑥. A state 𝑥= (𝑧+𝑒) ∈ X𝑁is said to be recursively feasible if V𝑁(𝑧+)≠∅ and 𝑒+∈𝛺𝑥. Theorem 1 (Recursive Feasibility).Consider that Assumptions 1–6hold. Suppose also 𝑥= (𝑧+𝑒) ∈ X𝑁at instant 𝑘so that (𝑧∗(0),𝐯∗), where 𝐯∗= {𝑣∗(0), 𝑣∗(1),…, 𝑣∗(𝑁−1)} and its associated nominal state trajectory 𝐳∗= {𝑧∗(0), 𝑧∗(1),…, 𝑧∗(𝑁)} are feasible for problem P∗ 𝑁(𝑥). Then, for all 𝑥+=𝐴𝑗𝑥+𝐵𝑗𝑢+𝑤with 𝑤∈, problem P∗ 𝑁(𝑥+)is also feasible for the next instants. Proof. Given a feasible state 𝑥∈𝑧∗(0) ⊕ 𝛺𝑥, the successor state 𝑥+∈𝑧∗(1) ⊕ 𝛺𝑥remains within the tube with the applied control 𝑢=𝑣∗(0) + 𝐾(𝑥−𝑧∗(0)) (for a detailed proof, see [43, Proposition 1]). Since there is a feasible solution (𝑧∗(0),𝐯∗) at instant 𝑘, the system trajectory toward the origin is recursively feasible for the next instants as it is the combination of database trajectories. Journal of Process Control 132 (2023) 103109 6 E. Masero et al. 5.2. Robust exponential stability Here, we analyze the robust stability of the controlled uncertain system 𝑥+=𝐴𝑗𝑥+𝐵𝑗𝑢+𝑤, ∀𝑤∈, which will inherit the properties of the database trajectories, so let us consider the next assumption to prove stability. Assumption 7. The trajectories in the database are generated using stabilizing controllers, and for all 𝑧∈Z𝑁, the value function 𝐽o 𝑁(𝑧) (recall Eq. (11)) is Lyapunov, i.e., 𝐽o 𝑁(𝑧+)≤𝐽o 𝑁(𝑧) − 𝓁(𝑧, 𝑣),(19) so there exist positive constants 𝑏 > 𝑎 > 0such that 𝐽o 𝑁(𝑧)≥𝑎‖𝑧‖2, 𝐽o 𝑁(𝑧)≤𝑏‖𝑧‖2, 𝐽o 𝑁(𝑧+)≤𝐽o 𝑁(𝑧) − 𝑎‖𝑧‖2. (20) In this way, the origin is exponentially stable for the controlled nominal dynamics 𝑧+=𝐴𝑗𝑧+𝐵𝑗𝑣for all 𝑧∈Z𝑁. Should there be trajectories that do not satisfy Assumption 7, they can be removed from the database. Although Assumption 7 is restrictive, it is necessary to inherit stability from the nominal trajectories generated by the data-based controller. Proposition 1. (i) For all 𝑥∈X𝑁,𝐽∗ 𝑁(𝑥) = 𝐽o 𝑁(𝑧∗(0)) and their corresponding control solutions 𝐯∗(𝑥) = 𝐯o(𝑧∗(0)). (ii) 𝐽∗ 𝑁(⋅)is a Lyapunov function for dynamics 𝑥+=𝐴𝑗𝑥+𝐵𝑗𝑢+𝑤, i.e., 𝐽∗ 𝑁(𝑥+)≤𝐽∗ 𝑁(𝑥) − 𝓁(𝑧∗(0), 𝑣∗(0)).(21) Proof. The result (i) follows from a direct consequence of the definition of problems P𝑁(𝑧)and P∗ 𝑁(𝑥). To prove (ii), consider Theorem 1 and note that 𝑥+=𝑧∗(1) ⊕ 𝛺𝑥such that (𝑧∗(1),𝐯o(𝑧∗(1))) is also feasible for problem P∗ 𝑁(𝑥+), but 𝐽o 𝑁(𝑧∗(1)) ≥𝐽∗ 𝑁(𝑥+). From (19), we have 𝐽o 𝑁(𝑧∗(1)) ≤𝐽o 𝑁(𝑧∗(0)) − 𝓁(𝑧∗(0), 𝑣∗(0)), which is lower bounded by 𝐽o 𝑁(𝑥+)≤𝐽o 𝑁(𝑧∗(1)) ≤𝐽o 𝑁(𝑧∗(0))−𝓁(𝑧∗(0), 𝑣∗(0)). Since 𝐽∗ 𝑁(𝑥) = 𝐽o 𝑁(𝑧∗(0)) by (i), we conclude that (21) holds. Theorem 2 (Robust Exponential Stability).Consider that Assumption 7 holds. For all 𝑥∈X𝑁and 𝑤∈, the control law 𝑢=𝑣∗(0) + 𝐾(𝑥−𝑧∗(0)) stabilizes the uncertain system 𝑥+=𝐴𝑗𝑥+𝐵𝑗𝑢+𝑤exponentially in the set 𝛺𝑥. Proof. As proved in [34], since 𝐽∗ 𝑁(𝑥) = 𝐽o 𝑁(𝑧∗(0)) (let us denote 𝑧∗(0) = 𝑧∗ 0(𝑥)to make explicit the dependence of the state 𝑥) for all 𝑥∈X𝑁, we obtain the following set of equations using (20) and (21): 𝐽∗ 𝑁(𝑥) = 𝐽o 𝑁(𝑧∗ 0(𝑥)) ≥𝑎‖𝑧∗ 0(𝑥)‖2,(22a) 𝐽∗ 𝑁(𝑥) = 𝐽o 𝑁(𝑧∗ 0(𝑥)) ≤𝑏‖𝑧∗ 0(𝑥)‖2,(22b) 𝐽∗ 𝑁(𝑥+) − 𝐽∗ 𝑁(𝑥)≤−𝓁(𝑧∗ 0(𝑥), 𝑣∗ 0(𝑥)) ≤−𝑎‖𝑧∗ 0(𝑥)‖2,(22c) where 𝑏 > 𝑎 > 0. Let 𝑥(𝑘) = 𝑥+denote the state system at instant 𝑘 given an initial state 𝑥(0) = 𝑥and 𝑤∈; from Theorem 1,𝑥(𝑘) ∈ X𝑁 for all 𝑘. Let us define 𝑆𝛼≜{𝑥∶𝐽∗ 𝑁(𝑥)≤𝛼}for all 𝛼≥0. Then, 𝑆0∈𝛺𝑥 for 𝛼= 0, and there is an 𝛼 > 0such that 𝑆𝛼⊂. From (22), we have 𝐽∗ 𝑁(𝑥(𝑘)) ≤𝜌𝑘𝐽∗ 𝑁(𝑥(0)), where 𝜌𝑘≜(1 − (𝑎∕𝑏)) ∈ (0,1). Iterating from 𝑘= 0 to 𝑘and following 𝐽∗ 𝑁(𝑥(𝑘)) ≤𝜌𝑘𝐽∗ 𝑁(𝑥(0)), we have that 𝑎‖𝑧∗ 0(𝑥(𝑘))‖2≤𝜌𝑘𝑏‖𝑧∗ 0(𝑥(0))‖2,e.g.,‖𝑧∗ 0(𝑥(𝑘))‖≤𝛿𝑘𝑐‖𝑧∗ 0(𝑥(0))‖, where 𝛿𝑘=√𝜌and 𝑐=√𝑏∕𝑎 < ∞for all 𝑥(0) = 𝑥∈𝛼. Since X𝑁 is bounded, it yields from (22c) the existence of a finite integer 𝐷 such that 𝑥(𝑘) ∈ 𝛼for all 𝑘∈𝐷,𝑥∈X𝑁and 𝑤∈. Due to 𝑥(𝑘) ∈ 𝑧∗ 0(𝑥(𝑘)) ⊕ 𝛺𝑥, the set 𝛺𝑥is robustly exponentially stable for the uncertain system with all 𝑥∈X𝑁. 6. Cooperative learning example A potential application, for which the proposed approach can be used as-is, is that of cooperative learning, where multiple controllers share their previous experiences to attain better performance. Let us consider a set of robots moving in a plane, so that their position and velocities are, respectively, controlled and manipulated variables. Moreover, the set of robots share a cloud-based database that collects information from different robots’ controllers. We may model such robots as a system with parametric uncertainties that is formed nominally by two discrete-time linear integrators as follows: [𝑥1 𝑥2]+ =𝐴(𝜃)[𝑥1 𝑥2]+𝐵(𝜃)[𝑢1 𝑢2],(23) with state 𝑥= [𝑥1, 𝑥2]⊤and control input 𝑢= [𝑢1, 𝑢2]⊤both subject to constraints: [−2,−2]⊤≤{𝑥, 𝑢}≤[2,2]⊤. Matrices 𝐴(𝜃)and 𝐵(𝜃)are, respectively, composed of the nominal matrices 𝐴= I and 𝐵= I plus parametric uncertainties 𝛥𝐴𝜃=[𝛥𝐴𝜃10 0𝛥𝐴𝜃2], 𝛥𝐵𝜃=[𝛥𝐵𝜃10 0𝛥𝐵𝜃2], where 𝛥𝐴𝜃1, 𝛥𝐴𝜃2, 𝛥𝐵𝜃1, 𝛥𝐵𝜃2∈ [0,0.4]. For simulations, we consider the realizations of the nominal system detailed in Table 2. The simulation length is 𝑁end = 30, the prediction horizon is 𝑁= 5, and the cost function is defined by (8) with 𝑄= 0.5⋅Iand 𝑅= 5 ⋅I. The optimal feedback controller 𝐾 for all system realizations (𝐴(𝜃), 𝐵(𝜃), 𝑄, 𝑅)is calculated by solving the set of LMIs (18) obtaining 𝐾= −0.7257 ⋅I. The RPI set 𝛺𝑥is calculated with the MPT toolbox of MATLAB taking into account the set defined by (5). Furthermore, we have created a database of 𝑇= 100 trajectories from the nominal system (𝐴, 𝐵)by using different proportional controllers. Note that we consider the most restrictive case in which there are no data from other realizations. These trajectories start from different feasible initial states and reach the origin, i.e.,𝑥obj = [0,0]⊤. Each trajectory has 𝑘𝑡 end–step information, where the value of each 𝑘𝑡 end depends on when its state error was below the threshold 𝜀= 10−3. Hence, there is a total of ∑𝑡∈𝑘𝑡 end partial trajectories, but the algorithm only considers 𝑀= 1500 to combine them. We test our robust data-based controller (Robust DBC) with the data-based controller (DBC) previously proposed in [14]. For initial state 𝑥(0) = [−1.8,1.8]⊤,Fig. 2 shows a comparison of state trajectories {𝑥(𝑘)} obtained by implementing both controllers for realizations {(𝐴, 𝐵),(𝐴1, 𝐵1),…,(𝐴4, 𝐵4)}. We also show all trajectories in the database (), the sequence of the optimized nominal state {𝑧∗ 0(𝑥(𝑘))} centered in the RPI set 𝛺𝑥, and the constraint set . Furthermore, for all considered system realizations, Table 2 displays a numerical comparison in terms of the maximum computing time to solve the control problem (denoted as 𝐶𝑇max [s]), the absolute average error (𝐸abs): 𝐸abs = 𝑁end ∑ 𝑘=1 ‖𝑥(𝑘) − 𝑥obj‖ 𝑁end , and the cumulative performance cost (𝐽acc): 𝐽acc = 𝑁end ∑ 𝑘=1 𝑥(𝑘)⊤𝑄 𝑥(𝑘) + 𝑢(𝑘)⊤𝑅𝑢(𝑘). As shown in Table 2, the implementation of the Robust DBC results in lower absolute average error and cumulative performance cost than the DBC method for all system realizations. Moreover, as depicted in Fig. 2a, the DBC method fails to control the system realizations (𝐴2, 𝐵2)and (𝐴4, 𝐵4)driving the states far from the origin. In contrast, the Robust DBC control method successfully steers the state toward the goal for all system realizations. These results Journal of Process Control 132 (2023) 103109 7 E. Masero et al. Fig. 2. State trajectories with initial state 𝑥(0) = [−1.8,1.8]⊤for realizations {(𝐴, 𝐵),(𝐴1, 𝐵1),…,(𝐴4, 𝐵4)} by implementing our Robust DBC and the DBC previously proposed in [14]. Table 2 R2.4 Numerical results comparison between the DBC method previously proposed in [14], and our Robust DBC strategy. System realizations System parametric uncertainties DBC Robust DBC Δ𝐴𝜃1Δ𝐴𝜃2Δ𝐵𝜃1Δ𝐵𝜃2𝐶𝑇max [s] 𝐸abs 𝐽acc 𝐶𝑇max [s] 𝐸abs 𝐽acc (𝐴, 𝐵)0 0 0 0 0.0036 0.13 11.17 0.0256 0.11 10.30 (𝐴1, 𝐵1)0.4 0 0.4 0 0.0022 0.21 20.83 0.0217 0.14 13.88 (𝐴2, 𝐵2)0.4 0 0 0.4 0.0015 2.35 613.49 0.0221 0.20 22.77 (𝐴3, 𝐵3)0.4 0 0.4 0.4 0.0022 0.18 18.66 0.0203 0.12 12.67 (𝐴4, 𝐵4)0 0.4 0.4 0 0.0013 2.69 870.78 0.0192 0.29 24.58 (𝐴5, 𝐵5)0 0.4 0 0.4 0.0014 0.17 14.94 0.0187 0.12 11.70 (𝐴6, 𝐵6)0 0.4 0.4 0.4 0.0009 0.17 13.85 0.0192 0.12 10.13 (𝐴7, 𝐵7)0.4 0.4 0.4 0 0.0017 0.60 55.45 0.0242 0.24 23.56 (𝐴8, 𝐵8)0.4 0.4 0 0.4 0.0013 2.04 369.57 0.0195 0.27 28.92 (𝐴9, 𝐵9)0.4 0.4 0.4 0.4 0.0013 0.27 25.82 0.0186 0.16 15.57 Fig. 3. State trajectory with initial state 𝑥(0) = [−1.8,1.8]⊤for realization (𝐴2, 𝐵2)by implementing our Robust DBC with a feedback gain 𝐾designed solving the LMIs [39] for a particular different realization. demonstrate that the original method can be very sensitive to problems such as uncertainties or noise in the data, where the proposed method is inherently more robust. Finally, Fig. 3 illustrates how the controller is capable to stabilize the system even if the feedback gain is poorly design for a different particular realization. Specifically, we show how the proposed method is applied to system realization: (𝐴2, 𝐵2). Here, it has been used as an auxiliary controller the feedback gain 𝐾= −0.2475 ⋅I, which was designed specifically for a different system realization within the uncertain set, e.g.,(𝐴10, 𝐵10): 𝐴10 =[1 0 0 1], 𝐵10 =[1.4 0 0 1.4], solving the LMIs [39]. This case is interesting because this feedback gain would yield unstable closed-loop dynamics with eigenvalues 0.6536 and 1.1525 if it was used directly as main controller. Therefore, this example illustrates one of the remarkable features of the proposed method, which is to require only a feedback gain that stabilizes at least one of the system realizations. Note also that the change of gain in this example motivates the different shape of the invariant set. 7. Nonlinear-dynamics example We consider the quadruple-tank plant presented in [44] as the case study with non-linear dynamics. The plant is formed by four interconnected tanks, where the upper tanks (#3 and #4) discharge flow rates to the lower tanks (#1 and #2), Journal of Process Control 132 (2023) 103109 8 E. Masero et al. Fig. 4. Scheme of the quadruple-tank plant. and these, in turn, into a sinking tank. The plant is controlled by two pumps that keep the water circulation between the tanks. There are also two three-way valves that divide the flow into two ways (see Fig. 4). Applying the mass balance and Bernoulli’s law to the plant, the non-linear dynamics are: 𝑆𝑑ℎ1 𝑑𝑡 = −𝑎1√2𝑔ℎ1+𝑎3√2𝑔ℎ3+𝛾a 𝑞a 3600 , 𝑆𝑑ℎ2 𝑑𝑡 = −𝑎2√2𝑔ℎ2+𝑎4√2𝑔ℎ4+𝛾b 𝑞b 3600 , 𝑆𝑑ℎ3 𝑑𝑡 = −𝑎3√2𝑔ℎ3+ (1 − 𝛾b)𝑞b 3600 , 𝑆𝑑ℎ4 𝑑𝑡 = −𝑎4√2𝑔ℎ4+ (1 − 𝛾a)𝑞a 3600 , (24) where 𝑆= 0.03 m2is the cross section of the tanks, {ℎ1, ℎ2, ℎ3, ℎ4}are the water levels of the tanks, and {𝑎1, 𝑎2, 𝑎3, 𝑎4}are the cross section of the outlet pipes. The parameters {𝛾a, 𝛾b} ∈ [0,1] refer to the opening of the three-way valves (𝛾a= 0.3,𝛾b= 0.4),𝑔= 9.81 m∕s2is the gravitational constant, and {𝑞a, 𝑞b}are the pump flow rate. The dynamics of the non-linear system (24) can be linearized by its Taylor expansion as: 𝑥+=𝐴 𝑥 +𝐵 𝑢, (25) where 𝑥 = [ℎ1−ℎ0 1,…, ℎ4−ℎ0 4]⊤and 𝑢 = [𝑞a−𝑞0 a, 𝑞b−𝑞0 b]⊤around the following reference point: ⎡⎢⎢⎢⎢⎣ ℎ0 1 ℎ0 2 ℎ0 3 ℎ0 4 ⎤⎥⎥⎥⎥⎦ =⎡⎢⎢⎢⎢⎣ 0.5005 0.4704 0.5206 0.4319 ⎤⎥⎥⎥⎥⎦ m,[𝑞0 a 𝑞0 b]=[1.5355 1.6794]m3∕h. Therefore, the discrete-time linear nominal system (𝐴, 𝐵) with a sampling time 𝑇s= 30 s is characterized by: 𝐴=⎡⎢⎢⎢⎢⎣ 0.6654 0 0.1919 0 0 0.5971 0 0.2250 0 0 0.7643 0 0 0 0 0.7077 ⎤⎥⎥⎥⎥⎦ , 𝐵 =⎡⎢⎢⎢⎢⎣ 0.0684 0.0179 0.0254 0.0868 0 0.1461 0.1644 0 ⎤⎥⎥⎥⎥⎦ . We consider the following state and control input constraints: 0 [m]≤{ℎ1, ℎ2, ℎ3, ℎ4}≤1.2 [m], 0 [m3∕h]≤{𝑞a, 𝑞b}≤3 [m3∕h].(26) Fig. 5. Set of 100 trajectories of the nominal plant (𝐴, 𝐵)that compose the database. 7.1. Control parameters and database The control objective is to drive the states of tanks #1 and #2 to the reference point: ℎ0≜[ℎ0 1 ℎ0 2]=[0.5005 0.4704]m,(27) while minimizing the control cost (8) subject to constraints (26). The prediction horizon is 𝑁= 5, matrices 𝑄and 𝑅are: 𝑄=⎡⎢⎢⎢⎢⎣ 1000 0100 0000 0000 ⎤⎥⎥⎥⎥⎦ , 𝑅 =[0.01 0 0 0.01], the length of simulations is 𝑁end = 2010 s, and the sampling time of the controller is 𝑇s= 30. The database is composed of 𝑇= 100 trajectories, which have been obtained from the nominal plant (𝐴, 𝐵)using PI controllers. These trajectories start at different points and end at the reference point (27), as illustrated in Fig. 5. Since the plant is a multiple-input multipleoutput (MIMO) system, we employ two discrete PI controllers, one to manage tank level ℎ1via pump 𝑞aand another to manage ℎ2through pump 𝑞b, whose parameters are randomized in the intervals: 𝐾p∈ [2.25,3.25] and 𝑇i∈ [275,475]. The procedure carried out to obtain the trajectories can be divided into two steps. First, the system is brought to a point that meets the constraints (26). Once in steady state, we apply the PI controller to bring the system to the reference point (27) until the sum of 100 consecutive absolute errors of ℎ1(𝑘)plus ℎ2(𝑘)is below the threshold 𝜀= 0.35. Note that despite using PIs to generate the database, we may encounter minor offset in some of the trajectories. Each trajectory has 7000 seconds of information on the plant operation. Therefore, there is a total of 𝑇⋅7000∕𝑇scandidate partial trajectories at each time step, but only 𝑀= 1500 partial trajectories are considered by the algorithm to compute the final control sequence. Furthermore, we consider five different system realizations, {(𝐴1, 𝐵1),…,(𝐴5, 𝐵5)}, to test the proposed control method, which are obtained by slightly changing the cross section of the outlet pipes of nominal system (𝐴, 𝐵), as detailed in Table 3. Note that, for a practical application of the control method, the difference between systems should be small enough to apply robustness arguments. The optimal feedback controller 𝐾for all system realizations obtained by solving LMIs (18) is 𝐾=[−2.9758 −0.7966 −0.1686 −1.5448 −1.7242 −1.7282 −2.3000 0.3660 ], and 𝛺𝑥is calculated with MPT toolbox taking into account 𝐾and the polytopic set  ≜ {𝑤∈R𝑛𝑤∶𝑤 < 𝟎.𝟏}. Journal of Process Control 132 (2023) 103109 9 E. Masero et al. Fig. 6. Results for system realization (𝐴2, 𝐵2)by implementing our Robust DBC and the previously proposed DBC. Table 3 R2.4 Numerical results comparison between the DBC method proposed in [14] and our Robust DBC strategy. System realizations Cross section of outlet pipes [m] DBC Robust DBC 𝑎1⋅10−4 𝑎2⋅10−4 𝑎3⋅10−4 𝑎4⋅10−4 𝐶𝑇max [s] 𝐸abs (%) 𝐽acc 𝐶𝑇max [s] 𝐸abs (%) 𝐽acc (𝐴, 𝐵)1.301 1.597 0.876 1.026 0.0078 2.162 0.131 0.0236 0.751 0.039 (𝐴1, 𝐵1)1.341 1.627 0.896 1.108 0.0021 2.162 0.116 0.0238 0.706 0.038 (𝐴2, 𝐵2)1.236 1.488 0.807 0.997 0.0023 2.534 0.137 0.0320 0.751 0.041 (𝐴3, 𝐵3)1.236 1.627 0.807 1.108 0.0023 1.341 0.089 0.0267 0.714 0.038 (𝐴4, 𝐵4)1.683 1.691 0.964 0.940 0.0156 2.460 0.175 0.0219 0.869 0.053 (𝐴5, 𝐵5)1.475 1.845 0.896 1.208 0.0022 2.197 0.116 0.0243 0.668 0.037 7.2. Simulation results We performed simulations for different system realizations to provide a numerical comparison between the data-based controller (DBC) presented in [14], and our robust data-based controller (Robust DBC). Fig. 6 shows the trajectories of tank levels and flow rates {𝑞a, 𝑞b}, considering the initial states (ℎ1(0) = 0.7005 m and ℎ2(0) = 0.6704 m, with their corresponding ℎ3and ℎ4levels) for the system realization (𝐴2, 𝐵2).Table 3 displays the mean absolute error (%) of the tank levels 𝑛= {1,2} calculated as: 𝐸𝑛 abs =∑𝑁end∕𝑇s 𝑘=1 |||ℎ𝑛(𝑘) − ℎ0 𝑛||| 𝑁end∕𝑇s ⋅100, 𝐸abs(%) = (𝐸1 abs +𝐸2 abs)∕2, and also the accumulated cost: 𝐽acc = 𝑁end∕𝑇s ∑ 𝑘=1  ℎ(𝑘)⊤𝑄 ℎ(𝑘) + 𝑞(𝑘)⊤𝑅 𝑞(𝑘), where  ℎ(𝑘)=[ℎ1(𝑘), ℎ2(𝑘)] − [ℎ0 1, ℎ0 2], 𝑞(𝑘)=[𝑞a(𝑘), 𝑞b(𝑘)] − [𝑞0 a, 𝑞0 b]. As displayed in Table 3, the mean absolute errors obtained with Robust DBC are lower than those obtained with the original algorithm. In terms of the accumulated cost, Robust DBC also outperforms the original algorithm for all system realizations. Finally, Figs. 6(a) and 6(b) show the simulated trajectories of state and input for system realization (𝐴2, 𝐵2). Our robust data-based controller is able to manage the plant with the information of the database (see Fig. 6(b)), while the original controller offers a worse control, as shown in Fig. 6(a). These results demonstrate the effectiveness of the proposed robust data-based approach in dealing with challenging aspects of systems, such as noise and parametric uncertainties. 8. Conclusions The proposed approach extends a previously presented historianbased method with tube-based MPC ideas to obtain new properties such as robust constraint satisfaction, also paving the way for extensions to linear parameter-varying systems and applications such as cooperative learning. The main advantage of the new method is that it is straightforward to implement, as it does not require additional experiments, since historical data is typically available in most control systems. Finally, our outcomes in the two case studies show promising results and highlight the potential for data exchange between different realizations of the system in this context. In particular, our proposal uses closed-loop trajectories and an auxiliary control law; two ingredients analogous to the ones of classical MPC methods (closed-loop predictions and semi-feedback strategies [45]) to reduce conservatism, exploiting the duality between parametric uncertainties and additive uncertainties in the proposed setup. The adopted approach somehow navigates between the modelbased and model-free paradigms. In particular, it is based on a dual control approach where a nominal and an auxiliary controller are