Model predictive control for tracking with implicit invariant sets

Automatica 179 (2025) 112436 Contents lists available at ScienceDirect Automatica journal homepage: www.elsevier.com/locate/automatica Brief paper Model predictive control for tracking with implicit invariant setsI Irene Luque a,∗, Paula Chanfreut b, Daniel Limón a, José M. Maestre a aDepartment of Systems and Automation Engineering, University of Seville, Seville, Spain bDepartment of Mechanical Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands a r t i c l e i n f o Article history: Received 19 June 2024 Received in revised form 13 February 2025 Accepted 15 May 2025 Available online 19 June 2025 Keywords: Model predictive control Implicit invariant sets Tracking systems a b s t r a c t This paper presents a model predictive control (MPC) technique for tracking with implicit terminal components. The controller formulation includes an artificial setpoint as decision variable, and the terminal constraint is defined implicitly for an augmented system that depends on the latter. In this respect, instead of computing an invariant terminal set, we consider an extended prediction horizon whose length can be bounded simply by solving LPs. This approach overcomes size-related limitations associated with the operations needed for computing invariant sets, also simplifying the offline MPC design. The proposed controller is able to drive large systems to admissible setpoints while guaranteeing recursive feasibility and convergence. Finally, the method is illustrated by an academic example, a mass–spring–damper system of variable-size and a more realistic case study of a drone. © 2025 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/). 1. Introduction The theoretical properties of MPC, such as recursive feasibility and stability, are typically guaranteed by the convenient design of a terminal cost function and a terminal positively invariant set, which becomes increasingly difficult to compute as the size of the system grows (Gilbert & Tan, 1991; Mayne, 2013, 2014). In particular, the explicit determination of invariant sets is challenging even for linear time-invariant systems due to the required computations of set intersections and pre-image sets. In fact, there exist scalable methods that find approximations of these invariant sets, e.g., by considering pre-defined polyhedron shapes (Trodden, 2016), ellipsoids based on linear matrix inequalities (Alamo, Cepeda, & Limon, 2005), inner-outer approximations (Comelli, Olaru, & Kofman, 2024), zonotopes (Morato, Cunha, Santos, Normey-Rico, & Sename, 2021) or data-driven approaches (Berberich, Köhler, Müller, & Allgöwer, 2021), or that generate implicit representations of them (Raković & Zhang, 2022, 2023; Wang & Jungers, 2020). Among them, we are particularly interested in the latter (Raković & Zhang, 2022, 2023), for it extends the prediction horizon guaranteeing that the final state belongs to an invariant set. IThis work is supported by the Spanish Training Program for Academic Staff under Grant (FPU21/05299), and by Grants PID2022-141159OB-I00 and PID2023152876OB-I00, funded by MCIN/AEI/10.13039/501100011033 and ERDF/EU. The material in this paper was partially presented at the 63rd IEEE Conference on Decision and Control (CDC), December 16–19, 2024, Milan, Italy. This paper was recommended for publication in revised form by Associate Editor Dominic Liao-McPherson under the direction of Editor Florian Dorfler. ∗Corresponding author. E-mail addresses: [email protected] (I. Luque), [email protected] (P. Chanfreut), [email protected] (D. Limón), [email protected] (J.M. Maestre). Moreover, some MPC formulations further complicate this issue, such as that of tracking (Ferramosca, Limón, Alvarado, Alamo, & Camacho, 2009; Limón, Alvarado, Alamo, & Camacho, 2008), which considers an augmented terminal system — enlarging its dimension — to deal with non-fixed setpoints. Specifically, the MPC for tracking formulation presented in Limón et al. (2008) and Ferramosca et al. (2009) adds an artificial steady state and input as decision variables in the optimization problem to relax the terminal constraint, and ensures recursive feasibility and convergence to the real setpoint by using a modified cost function. While the family of MPC controllers to track changing setpoints is wider, e.g., Bemporad, Casavola, and Mosca (1997), Garone, Di Cairano, and Kolmanovsky (2017), Gilbert and Kolmanovsky (2002), we consider the approach in Ferramosca et al. (2009), Limón et al. (2008) to apply the proposed implicit methodology not only due to its theoretical guarantees, but also because it enlarges the domain of attraction of the controller. The main contribution of this article consists in designing an MPC for tracking that incorporates implicit terminal ingredients, extending the preliminary work introduced in Luque, Chanfreut, Limón, and Maestre (2024). Particularly, the proposed approach relies on replacing the terminal constraint set with an extended prediction horizon of a predefined finite length. The presented methodology allows obtaining this length for a tracking setting by solving linear programs (LPs), thus avoiding set calculations and therefore enabling its application to systems of any size, at the expense of a marginal increase in the online computational burden caused by the use of longer horizons. As will be seen, the use of implicit invariant sets is not straightforward in this context and requires a tailor-made adaptation of the results for regulation regarding the existence of artificial and real setpoints. https://doi.org/10.1016/j.automatica.2025.112436 0005-1098/© 2025 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/). I. Luque, P. Chanfreut, D. Limón et al. Automatica 179 (2025) 112436 The proposed methodology provides guarantees of recursive feasibility and convergence to the real setpoints, which are key properties in the use of tracking MPC. Furthermore, this paper discusses an alternative approach for cases where the maximum admissible length of the extended prediction horizon is limited, requiring alternative strategies to avoid the calculation of the terminal region while enabling to select the desired extension for the extended prediction horizon. The outline of the rest of the article is as follows. Section 2 introduces the preliminaries of the problem presented. Section 3 presents the controller design with implicit terminal components, as well as the alternative approach and theoretical proofs, whose performance is illustrated in Section 4 with three case studies. Lastly, Section 5 provides the discussion. Notation. Vector [x⊤,u⊤]⊤ is denoted as (x,u); In and 0m×n represent the identity and zero matrices of dimension n×n and m×n, respectively, whereas 0n and 1n are column vectors of zeros and ones of size n×1. R, N denote the sets of real and natural numbers, respectively. Likewise, given a,b∈N, with a<b, we define N[a,b]:= {a,a+1, . . . , b−1,b} and Nb is given for N[0,b]. Finally, [ut]T t=0 denotes vector [u⊤ 0,u⊤ 1, . . . , u⊤ T]⊤ for any given T∈N. The support function h(X,·) of a closed, nonempty subset X∈Rn is given for all y∈Rn by h(X,y):= supx{y⊤x:x∈X}. 2. Problem setting In this section, the system dynamics and the MPC for the considered tracking formulation are introduced. 2.1. System dynamics Consider a discrete-time linear system given by xk+1=Axk+Buk,(1) where xk∈Rn and uk∈Rm are respectively the state of the system and the input at instant k, and matrices A and B are of compatible dimensions, i.e., A∈Rn×n and B∈Rn×m, with m and n being positive and possibly different integers. Also, consider the following constraints xk∈X,uk∈U,∀k∈N,(2) being X⊆Rn and U⊆Rm the state and input constraint sets, respectively. Let us introduce the following assumption: Assumption 1. For system (1) subject to constraints (2), the following holds: •The matrix pair (A, B) is known and it is strictly stabilizable. •Constraint sets X and U are convex polytopic sets containing the origin in their interior. •There exists a feedback gain K∈Rm×n such that matrix A+BK is Schur, and a positive definite matrix P∈Rn×n such that (A+BK)⊤P(A+BK)−P= −(Q+K⊤RK).(3) The system performance will be evaluated through stage cost function ℓ(xk,uk,xs,us)= ∥xk−xs∥2 Q+ ∥uk−us∥2 R,(4) where Q∈Rn×n and R∈Rm×m are symmetric positive definite matrices, and (xs,us) denotes the setpoint to which we want to drive the system. In this regard, notice that any setpoint of the system must satisfy the following equation [A−InB][xs us]=0n.(5) Therefore, we can parameterize the pair (xs,us) through variable θ∈Rm, i.e., [xs us]=[Mθx Mθu]    Mθ θ, (6) being Mθ a suitable basis for the null space of [A−InB] that aggregates matrices Mθx∈Rn×m and Mθu∈Rm×m (Ferramosca et al., 2009). 2.2. MPC for tracking with explicit terminal components The MPC for tracking formulation considered in this article is characterized by the following (Limón et al., 2008): (i) An artificial setpoint, say (xa s, ua s), is introduced as optimization variables in the MPC problem. This artificial setpoint will be parametrized by variable θa and introduces m new optimization variables. (ii) An offset cost function is also added to penalize the deviation of the artificial setpoint from the real setpoints, (xr s, ur s). (iii) An augmented terminal invariant set Ψtr f is used, which is defined for augmented system [xk+1 θa]=[A+BK BL 0m×nIm]    Aaug [xk θa],(7) where L= [−K Im]Mθ. In what follows, we will differentiate between the real setpoint, denoted as (xr s, ur s), and the artificial setpoint, i.e., (xa s, ua s). Note that (7) considers that system (1) employs control law uk=ua s+K(xk−xa s)=Kxk+Lθa.(8) Considering the above, the MPC for tracking problem to be solved at every time instant k adopts the following form: V∗ N(xk,xr s)=min u,θaVN(xk,xr s,u, θa) s.t. x0|k=xk,(9a) xj+1|k=Axj|k+Buj|k,j∈N[0,N−1],(9b) xj+1|k∈X,j∈N[0,N−1],(9c) uj|k∈U,j∈N[0,N−1],(9d) [xa s ua s]=Mθθa,[xN|k θa]∈Ψtr f,(9e) where N is the length of the prediction horizon and u= [uj|k]N−1 j=0. In this respect, subscript j|k indicates a prediction on the corresponding variable for instant k+j made at time k. Also, the cost function is defined as VN(xk,xr s,u, θa)= N−1 ∑ j=0(∥xj|k−xa s∥2 Q+ ∥uj|k−ua s∥2 R) + ∥xN|k−xa s∥2 P+ ∥xa s−xr s∥2 O, (10) where O∈Rn×n is a positive definite matrix and P satisfies (3). Note that, unlike MPC for regulation, the deviation of the system with respect to the artificial setpoint is weighted during the prediction horizon, and an offset cost, i.e. ∥xa s−xr s∥2 O, is added to penalize the difference between the artificial and real state reference. Likewise, terminal constraint (9e) is defined by invariant set Ψtr f, which is computed for augmented system (7) considering constraints (2), as detailed in Limón et al. (2008). 2 I. Luque, P. Chanfreut, D. Limón et al. Automatica 179 (2025) 112436 In particular, Ψtr f is a polyhedral approximation to the maximal invariant set, satisfying Ψtr f⊆ {(x, θ)∈Rn+m:(x,Kx +Lθ)∈(X,U), θ ∈Θ}, where Θ:= {θ∈Rm:Mθxθ∈X,Mθuθ∈U}.(11) Because of the unitary eigenvalues of Aaug, set Ψtr f might not be finitely determined (Gilbert & Tan, 1991). Nonetheless, it is possible to scale Θ by factor λ∈(0,1) so that the maximal admissible invariant set becomes a finitely determined convex polyhedron, say Ψtr f,λ (Gilbert & Tan, 1991; Limón et al., 2008). Note that Ψtr f,λ ⊆ {(x, θ)∈Rn+m:(x,Kx +Lθ)∈(X,U), θ ∈λΘ},(12) and [Axk+B(Kxk+Lθa) θa]∈Ψtr f,λ for all [xk θa]∈Ψtr f,λ. As detailed in Limón et al. (2008) and Ferramosca et al. (2009), the tracking formulation (9) allows driving the system state to any admissible target setpoint. However, it requires computing invariant set Ψtr f,λ for augmented dynamics (7), whose dimension is n+m. While the terminal cost is simple to calculate, the construction of such a terminal set becomes intractable for large systems. To solve this, this article reformulates MPC problem (9) using implicit terminal components, i.e., avoiding the need of explicitly characterizing Ψtr f,λ. 3. MPC for tracking with implicit terminal components In what follows, we extend results on implicit terminal components derived for regulation in Raković and Zhang (2023) to tracking problems. 3.1. Implicit terminal set for regulation For the regulation problem, that is (xr s,ur s)=0n+m, the terminal control law can be simply defined as uk=Kxk, and, therefore, the terminal dynamics are given by xk+1=(A+BK)xk,(13) which is also obtained if we fix the artificial and real setpoint to the origin in (7) and (8). Also, given (2), the constraints for the terminal stage can be compactly defined as xk∈Xt:= {x∈Rn:x∈X,Kx ∈U}.(14) Considering the above, let us introduce the following theorem, which establishes a sufficient condition for the existence of the maximal positively invariant set. Recall that a set Ω⊂Rn is defined as a positively invariant set for constraints x∈X if and only if (Kerrigan, 2001; Raković & Zhang, 2022): xk∈Ω⇒ ∃uk∈U such that xk+1∈Ω,xk+1∈X. Then, a set is defined as the maximal positively invariant set if it is positively invariant and contains all the positively invariant sets in Ω (Kerrigan, 2001). Theorem 1 (Raković and Zhang (2023, Theorem 1 and 2)). Suppose Assumption 1 holds. Then, the maximal positively invariant set for system (13) and constraints (14) is finitely determined if and only if for some M∈N some of the following holds M ⋂ j=0 (A+BK)−jXt⊆(A+BK)−(M+1)Xt or Xt⊆(A+BK)−(M+1)Xt. (15) Given Theorem 1, the maximal positively invariant set for system (13) and constraints (14), say Ψf, is a nonempty closed polyhedral set containing the origin in its interior, defined as Ψf= M ⋂ j=0 (A+BK)−jXt.(16) An alternative approach to check whether a given state belongs to Ψf can be introduced using a trajectory of length M (Raković & Zhang, 2023). That is, xN∈Ψf if sequence [xj]N+M j=N is such that ∀j∈N[N,N+M],xj∈Xt, ∀j∈N[N,N+M−1],xj+1=(A+BK)xj,(17) where M,N∈N, with M≥1 satisfying (15). Note that if M=0, then Ψf=X. The latter serves as the basis for the implicit reformulation of the terminal components. Let set Xt⊆Rn (recall (14)) be a closed polyhedral set whose irreducible representation is Xt:= {x∈Rn:(C+DK)x≤1p},(18) where matrix pair (C,D)∈Rp×n×Rp×m is defined accordingly. Note that p denotes the number of inequalities defining Xt. Given (18), Xt can be similarly defined as Xt= {x∈Rn: ∀i∈N[1,p],X⊤ ix≤1},(19) where X⊤ i is the ith row of the matrix (C+DK). Following Raković and Zhang (2023), we have that (15) holds true if and only if one of the following holds for all i∈N[1,p]: h( M ⋂ j=0 (A+BK)−jXt,((A+BK)M+1)⊤Xi)≤1 or h(Xt,((A+BK)M+1)⊤Xi)≤1. (20) The left-hand sides of the inequalities in (20) can be calculated for different i through the following computationally simple LPs, respectively: sup (x,[xj]M j=0){X⊤ i(A+BK)M+1x:xj=(A+BK)jx,xj∈Xt,∀j∈NM} or sup x {X⊤ i(A+BK)M+1x:x∈Xt}. Then, given any integer M, verifying (15) reduces to solving p LPs. In addition, the search for such an integer M can be performed by searching through any suitably generated sequence of positive increasing integers. See Raković and Zhang (2022, 2023) for further details. 3.2. Implicit terminal set for tracking Assume that system (1) is controlled by control law (8). Then, considering the augmented dynamics, the following holds for a given θ: [xk+1 θ]=Aaug [xk θ],(21) where the constraints for this terminal augmented dynamics are defined as follows Xaug,t:= {(x, θ)∈Rn+m:x∈X,Kx +Lθ∈U, θ ∈λΘ}.(22) Similarly to (18), let set Xaug,t⊆Rn+m be a closed polyhedron that can be defined irreducibly as Xaug,t:= {(x, θ)∈Rn+m:˜ G[x θ]≤1˜ p},(23) 3 I. Luque, P. Chanfreut, D. Limón et al. Automatica 179 (2025) 112436 where ˜ G=[˜ C+˜ DK ˜ DL 0˜ W]∈R˜ p×(n+m).(24) Also, ˜ C=C,˜ D=D,˜ W=(CMθx+DMθu)/λ and 0 is a matrix of zeros of the appropriate size. Remember that ˜ p is defined as the number of inequalities in Xaug,t. Given (23), set Xaug,t can also be defined as Xaug,t={(x, θ)∈Rn+m: ∀i∈N[1,˜ p],X⊤ aug,i[x θ]≤1},(25) where X⊤ aug,i is the ith row of the matrix ˜ G. In addition, let us consider some ˜ M∈N satisfying a condition similar to (15) but adapted for the augmented system, i.e., ˜ M ⋂ j=0 A−j augXaug,t⊆A−(˜ M+1) aug Xaug,t or Xaug,t⊆A−(˜ M+1) aug Xaug,t, (26) and the following terminal set Ψtr f,λ = ˜ M ⋂ j=0 A−j aug [Xt λΘ].(27) Given that, for some λ∈(0,1), the augmented terminal set for tracking Ψtr f,λ is finitely determined (Gilbert & Tan, 1991; Limón et al., 2008), it is possible to ensure that there exists a finite value of ˜ M that satisfies Eq. (26) (see also Raković & Zhang, 2022, Corollary 4). Theorem 2. Suppose Assumption 1 holds, and consider some ˜ M,N∈N, with ˜ M≥1 satisfying (26). Then, constraint (xN, θ)∈ Ψtr f,λ is satisfied if and only if there exists a sequence [xj]N+˜ M j=N such that ∀j∈N[N,N+˜ M],(xj, θ)∈Xaug,t, ∀j∈N[N,N+˜ M−1],xj+1=(A+BK)xj+BLθ. (28) To fulfill (26), one of the following conditions must hold, respectively, for all i∈N[1,˜ p]: h( ˜ M ⋂ j=0 A−j augXaug,t,(A˜ M+1 aug )⊤·Xaug,i)≤1 or h(Xaug,t,(A˜ M+1 aug )⊤·Xaug,i)≤1. (29) The following computationally simple LPs allow us to check the inequalities in (29) for any i∈ [1,˜ p]: sup (x,θ,[xj]˜ M j=0,[θj]˜ M j=0){X⊤ aug,i·A˜ M+1 aug [x θ]: [xj θj]=Aj aug [x θ]∈Xaug,t∀j∈N[0,˜ M]} or sup (x,θ) {X⊤ aug,i·A˜ M+1 aug [x θ]:(x, θ)∈Xaug,t}. As a result, a new parameter ˜ M satisfying (26) can be found by searching through positive increasing integers and solving ˜ p LPs for each of them. 3.3. MPC for tracking with implicit terminal components The proposed MPC, including tracking and implicit terminal constraints, is detailed in this subsection. Following Theorem 2, the prediction horizon is partitioned in two stages: a first part of length N, and a terminal part of length ˜ M. That is, sequences of length N+˜ M will be computed. Based on (9), the proposed optimal control problem to be solved at every time step is defined as min u,θaVN+˜ M(xk,xr s,u, θa) s.t. x0|k=xk,(30a) xj+1|k=Axj|k+Buj|k,j∈N[0,N+˜ M−1],(30b) xj+1|k∈X,j∈N[0,N+˜ M−1],(30c) uj|k∈U,j∈N[0,N+˜ M−1],(30d) uj|k=Kxj|k+Lθa,j∈N[N,N+˜ M−1],(30e) [xa s ua s]=Mθθa(30f) θa∈λΘ.(30g) Here, constraint (30e) imposes the terminal control law and is only considered during the second part of the prediction horizon, previously called the terminal stage. Also, the objective function takes into account the performance up to prediction time instant N+˜ M plus the terminal and the offset costs, as it was defined in (10). Finally, note that the extended prediction horizon replaces the explicit terminal constraint (9e). Remark 1. The implicit approach avoids computing offline the explicit representation of the invariant set in exchange for extending the prediction horizon, which implies increasing the online computational burden. However, since QPs can be solved in polynomial time, this may not be a limiting factor in most real applications. Also, notice that once ˜ M is known, the maximal invariant set in closed-form can also be found using (27), avoiding the need to check convergence conditions at every iteration. Remark 2. An implicit terminal cost function for tracking can be defined as follows, inferred from Raković and Zhang (2023): Vtr F,imp(xN|k,xa s,ua s)=(1 −ϵ)−1 N+˜ M−1 ∑ j=N (∥xj|k−xa s∥2 Q+ ∥uj|k−ua s∥2 R), where ϵ∈ [0,1) is minimized so that the implicit bound in the terminal cost becomes tight. Although this is a valid option, computing P as considered in (10) is typically not expensive and provides and more accurate estimation of the cost-to-go. For this reason, it is chosen in this paper. 3.4. Theoretical properties Hereafter, we prove that the initial feasibility of optimization problem (30) also implies recursive feasibility. In addition, convergence to the real setpoints, if admissible, is also proven. Theorem 3. Assume that at instant k there exists a solution (u∗ k, θ∗ k) of problem (30). Then, we can find a feasible solution of (30) at all instants t≥k. Proof. Let xk be the system state at instant k, and consider solution (u∗ k, θ∗ k). Notice that sequence u∗ k, and its associated predicted 4 I. Luque, P. Chanfreut, D. Limón et al. Automatica 179 (2025) 112436 state sequence are given by u∗ k=(u∗ 0|k,u∗ 1|k, . . . , u∗ N|k, . . . , u∗ N+˜ M−1|k), x∗ k=(x∗ 0|k,x∗ 1|k, . . . , x∗ N|k, . . . , x∗ N+˜ M|k), where x∗ 0|k=xk∈X. By construction, it follows that ∀j∈N[0,N+˜ M−1],x∗ j+1|k=Ax∗ j|k+Bu∗ j|k,(31a) ∀j∈N[N,N+˜ M−1],u∗ j|k=Kx∗ j|k+Lθ∗ k,(31b) ∀j∈N[0,N+˜ M−1],x∗ j+1|k∈X,u∗ j|k∈U.(31c) Let us define a candidate solution (˜ uk+1,˜ θk+1) for problem (30) at instant k+1. In particular, consider ˜ θk+1=θ∗ k and ˜ uk+1= [˜ uj|k]N+˜ M−1 j=0=[[u∗ j|k]N+˜ M−1 j=1 Kx∗ N+˜ M|k+Lθ∗ k], ˜ xk+1= [˜ xj|k]N+˜ M j=0=[[x∗ j|k]N+˜ M j=1 Ax∗ N+˜ M|k+BKx∗ N+˜ M|k+BLθ∗ k], (32) where ˜ xk+1 is the associated predicted state sequence. Note that, given x∗ 0|k=xk, and considering (31a), we have that xk+1= Axk+Bu∗ 0|k=x∗ 1|k. By construction, we have that ˜ xj|k+1∈X and ˜ uj|k+1∈U for all j∈N[0,N−2]. Likewise, ˜ xN−1|k+1=x∗ N|k∈Xf⊆X, ˜ uN−1|k+1=Kx∗ N|k+Lθ∗ k∈U. It is also known that constraint ˜ xN|k+1∈Xf is then implicitly applied given (16) and (32). Hence, during the terminal stage, where j∈N[N,N+˜ M−1], it holds that ˜ xj+1|k+1∈Xf⊆X, ˜ uj|k+1=K˜ xj|k+1+L˜ θk+1∈U.(33) Therefore, candidate solution (˜ uk+1,˜ θk+1) is a feasible solution of problem (30) at instant k+1. By induction, recursive feasibility is guaranteed for all time instants. ■ Theorem 4. Let x0∈XN+˜ M, with XN+˜ M being the domain of attraction of controller (30). Then, state xk of the system controlled by (30) will converge to xr s, if admissible, as k tends to infinity. Proof. Note that XN+˜ M represents the set of states for which optimization problem (30) is feasible. From the above, it follows that if xk∈XN+˜ M, then state xk+1 satisfies xk+1∈XN+˜ M. Consequently, XN+˜ M is a positively invariant set for the closedloop system. Likewise, given that X is bounded, set XN+˜ M is also bounded, thereby implying stability of the system. Below, we demonstrate convergence by verifying that the optimal cost is a Lyapunov function for the closed-loop system, and that the chosen artificial setpoint converges to the real one. Let VN+˜ M(xk+1,xr s,˜ uk+1,˜ θk+1) be the cost at instant k+1 associated with candidate solution (˜ uk+1,˜ θk+1) (see (32)). Also, notation zj|k=(xj|k,uj|k) is introduced for clarity. Then, we have that: VN+˜ M(xk+1,xr s,˜ uk+1,˜ θk+1)−VN+˜ M(xk,xr s,u∗ k, θ∗ k)= ℓ(˜ zk+N+˜ M|k+1,˜ θk+1)−ℓ(z∗ k|k, θ∗ k)+VF(˜ zk+N+˜ M+1|k+1,˜ θk+1)− VF(z∗ k+˜ M+N|k, θ∗ k), where ℓ(·) denotes the stage cost function as in (4), and VF(·) is the terminal cost in (30) (recall (10)). Then, given (33), it is possible to state that VN+˜ M(xk+1,xr s,˜ uk+1,˜ θk+1)−VN+˜ M(xk,xr s,u∗ k, θ∗ k)≤ −ℓ(z∗ k|k, θ∗ k). By applying the principle of optimality, we have that: VN+˜ M(xk+1,xr s,u∗ k+1, θ∗ k+1)−VN+˜ M(xk,xr s,u∗ k, θ∗ k)≤ −ℓ(z∗ k|k, θ∗ k). From this, it follows that the optimal cost is strictly decreasing and provides a Lyapunov function of the system. Given this, we infer that limk→∞ ∥xk−xa,∗ s,k∥Q=0. Note that xa,∗ s,k=Mθxθ∗ k is the chosen artificial state at instant k. Finally, in Limón et al. (2008, Lemma 3) and Ferramosca et al. (2009, Lemma 2), it is proved by contradiction that if xk=xa,∗ s,k, then ∥xk−xr s∥O=0. Consequently, convergence of the system state to xr s, if admissible, is demonstrated. ■ 3.5. Alternative implicit design for tracking This section introduces an alternative approach for the design of stabilizing predictive controllers without terminal constraint following Limón, Ferramosca, Alvarado, and Alamo (2018). This alternative is particularly relevant for scenarios where extending the prediction horizon by a length ˜ M is not feasible or desired. Additionally, for such systems, the explicit computation of the invariant set may also prove to be computationally intractable. Let us consider a terminal cost function similar to the one in (10), i.e. VF(xN+¯ M|k, θa)= ∥xN+¯ M|k−xa s∥2 P, as well as the terminal control law defined in (8). It should be noted that the new variable ¯ M is different from M and ˜ M, and can be selected for each system based on the desired performance or any other specific requirements. Also, choose some scalar α > 0 and define set: Ψα= {(x, θ)∈Rn+m:VF(x, θ)≤α},(34) such that Ψα is an invariant set for tracking. Finally, define the following MPC problem, which is similar to (30) but considers an user-defined horizon extension and adds a new parameter, γ, in the objective function: min u,θaVγ N+¯ M(xk,xr s,u, θa) s.t. x0|k=xk,(35a) xj+1|k=Axj|k+Buj|k,j∈N[0,N+¯ M−1],(35b) xj+1|k∈X,j∈N[0,N+¯ M−1],(35c) uj|k∈U,j∈N[0,N+¯ M−1],(35d) uj|k=Kxj|k+Lθa,j∈N[N,N+¯ M−1],(35e) [xa s ua s]=Mθθa(35f) θa∈λΘ.(35g) Particularly, Vγ N,¯ M(xk,xr s,u, θa) represents the cost function with the terminal cost scaled by γ, that is, γVF(xN+¯ M|k, θa). For the controller above, with γ≥1, recursive feasibility and convergence to admissible setpoints are ensured for all x∈Υ¯ M,γ (xr s) (Limón et al. (2018, Theorem 3)), with Υ¯ M,γ (xr s) being defined as: Υ¯ M,γ (xr s)= {x∈Rn:Vγ∗ N+¯ M(x,xr s)−V∗ O(x,xr s)≤(N+¯ M)d+γ α}, where Vγ∗ N+¯ M(x,xr s) is the optimal value of the objective function in (35), V∗ O(x,xr s) is the associated offset cost, and d represents a positive scalar such that ℓ(x,u,xs,us)≥d for all (x, θ)/∈Ψα. Notice that region Υ¯ M,γ (xr s) is enlarged as ¯ M or γ increase. Finally, it is worth mentioning that this controller design only requires selecting the values for ¯ M, γ, and α, avoiding any explicit or implicit estimation of the terminal invariant set. 5 I. Luque, P. Chanfreut, D. Limón et al. Automatica 179 (2025) 112436 Table 1 Cumulative costs for different values of λ. λAcademic example Drone ˜ MCumulative cost (×106)˜ MCumulative cost (×106) 0.99 10 2.0067 203 3.1548 0.89 5 2.4813 88 3.1548 0.79 4 3.0960 64 3.1548 0.69 3 3.8597 50 11.5174 0.59 3 4.7774 39 44.8102 0.49 3 5.8491 31 112.0422 0.39 2 7.0748 25 271.2111 Fig. 1. State evolution on plane (x1,x2). 4. Examples The performance of the controller will be illustrated in three examples of different dimensions and dynamics, simulated using YALMIP with solver GUROBI (Löfberg, 2004). 4.1. Academic example The first example is a low-dimensional system, whose dynamics are defined by matrices A and B in Limón et al. (2008). The system is constrained by ∥xk∥∞≤3 and ∥uk∥∞≤2, for all k≥0. For the controller design, weighting matrices Q=10 I2, R=100 I2 and O=10.000 I2 are used; N is set to 10; and gain K is obtained from the discrete LQR solution. Likewise, the value of ˜ M that satisfy the presented conditions for λ=0.9 is ˜ M=5. The state trajectory on the plane (x1,x2) is shown in Fig. 1, with x1 and x2 being the two components of the system state. We have considered different setpoints (represented with green dots), with the last of them being not admissible. The invariant set for regulation is also shown in Fig. 1, where it can be seen that it is contained into the projection of our augmented invariant set Ψtr f,λ onto the plane (x1,x2). Also, Table 1 presents a comparison of the cumulative performance cost for different values of λ. Specifically, the cumulative performance costs are computed using the following performance indicator: Vcc = Tsim ∑ k=0 (∥xk−xa s,k∥2 Q+ ∥uk−ua s,k∥2 R)+ Tsim ∑ k=0 ∥xa s,k−xr s∥2 O, where (xa s,k,ua s,k) denotes the artificial setpoint computed at instant k and Tsim is used to denote the number of simulated time instants (210 in this example). It is clear that the cumulative cost increases as λ decreases. This is expected since λ scales the set Θ, thus restricting the admissible values for θa. Likewise, an interesting property is inferred from Table 1: ˜ M decreases with the reduction of λ, for this causes the reachable setpoints to move away from the constraints. This however may result in certain equilibrium points of the system not being reached. Although, in terms of performance, the most convenient option is to choose λ close to 1, this observation provides a new degree of freedom. Specifically, if we determine the minimum λ required to parameterize the real setpoints of interest for the system, we can potentially reduce the necessary ˜ M. The time required to compute offline the terminal horizon length ˜ M was 0.5874s, whereas the time to explicitly compute the maximal positively invariant set was 0.6315s. As for the online control, the average time to solve problem (30) was 0.0047s, whereas to solve (9) we needed 0.0044s on average. Recall that (9) refers to the control problem without the extended prediction horizon and an explicit terminal set. While the simplicity of this academic example does not allow to show significant computational benefits, as it might reasonably be expected, it demonstrates the proposed design’s suitability for tracking purposes while opening up the possibility for scaling to more complex systems, where computing invariant sets is more challenging. Finally, the alternative design proposed in Section 3.5 has been implemented, and some key results are presented in Fig. 2. To illustrate the benefits of this method, we selected a value of ¯ M smaller than ˜ M=5. Specifically, we consider ¯ M=2, γ=20, and α=25. With these parameters, the resulting value of d was 5.6322. Fig. 2 shows the stability region Υ¯ M,γ (xr s) for the admissible target state xr s= [−0.5688,−0.0523]⊤. This indicates that if the initial state of the system lies within this region, the MPC for tracking without terminal constraint 3.5 will asymptotically stabilize the system. Notably, set Υ¯ M,γ (xr s) is practically as large as the projection of invariant set Ψtr f,λ. 4.2. Mass–spring–damper system of increasing size This subsection applies the proposed MPC to a modified version of the system in Riverso and Ferrari-Trecate (2012), Trodden and Maestre (2017). It consists on several carts connected by a spring–damper structure, as shown in Chanfreut, Maestre, Ferramosca, Muros, and Camacho (2021). The dynamics of each cart i are modeled by: [˙ ri ˙vi]=[0 1 −1 mikij −1 mihij]    Aii [ri vi]+[0 50]ui+wi,(36) wi=∑ j∈Ni[0 0 1 mikij 1 mihij]    Aij [rj vj],(37) 6 I. Luque, P. Chanfreut, D. Limón et al. Automatica 179 (2025) 112436 Table 2 Computation times for different numbers of carts. Offline times refer to the computation of the maximal RPI or the calculation of ˜ M for the explicit and implicit cases, respectively. Online times refer to the average time of each MPC optimization. Number of carts 3 5 10 15 25 30 50 70 100 200 Offline explicit time (s) 7.6293 15.5986 54.4716 153.2055 541.9198 799.4687 3.0095×1038.7399×1036.9640×104– Offline implicit time (s) 2.7063 3.7123 6.8664 21.5659 52.3099 73.6518 201.9882 390.5777 881.0654 4.6101×103 Online explicit time (s) 0.0085 0.0107 0.0761 0.1170 0.2924 0.3138 0.2052 0.2010 0.3200 – Online implicit time (s) 0.0095 0.0189 0.1122 0.1256 0.1582 0.1904 0.3101 0.4831 0.7917 2.0693 Fig. 2. Stability region reached by applying the alternative design method. where the state is defined by the displacement of cart i from an equilibrium position, ri, and its velocity, vi; and input ui represents a force that can be applied on cart i. There exist coupling terms between adjacent carts due to the springs and dampers that binds them together. In this respect, kij =kji and hij =hji represent the spring stiffnesses ([N/m]) and damping factors ([N/(m×s)]), and mi is the mass ([kg]) of cart i. The continuous-time dynamics are discretized using zero-order hold and a sampling time of 0.02s. We simulated different system sizes by progressively increasing the number of carts, say Ncarts, from 3 to 200. In all cases, the objective was to control all carts towards a target setpoint while satisfying the following constraints in state and control input: |ri| ≤ 4, |vi| ≤ 1 and |ui|<1, for all i∈ {1,2, . . . , Ncarts}. The weighting matrices were Qi= [1 0;0 1.5] and Ri=20, for all i∈ {1,2, . . . , Ncarts}, and λ=0.9. As can be seen in Table 2, the time to explicitly compute offline the invariant set for the system increases significantly, eventually becoming non-viable. In contrast, the time to compute ˜ M with the proposed implicit method remains tractable in all cases. It is important to note that the offline computation time was limited to a maximum of 48 h (2 days). Also, the online computation times are calculated as the average time to solve (30) over the entire simulation for each number of carts, capturing the overall increasing trend. The difference in the online times between the two methods increases progressively as the size of the system does, although it never becomes a limiting factor for the proposed approach. For instance, for a system comprising 100 carts, the offline computation time for the explicit method is 79 times higher than that of the analogous measure for the proposed method, whereas the online computation time for our method is only 2.5 times higher than for the explicit approach, as shown in Table 2. Fig. 3. Quadrotor position trajectory along the simulation. 4.3. Drone The model of a drone with 12 state variables described in Beard (2008), Romagnoli, Krogh, de Niz, Hristozov, and Sinopoli (2023) is employed now to illustrate the applicability of the proposed MPC method to real-world systems. Let us briefly indicate that the state of the system aggregates the position in [m] and linear velocities in [m/s] in a three dimensional space, say (px,py,pz) and (vx, vy, vz), respectively; and also the angles roll, pitch and yaw, in [rad], and their corresponding angular velocities, in [rad/s]. Likewise, u contains the thrust, expressed in [N], and the torques τφ,µ,ψ in [N· m], which are associated with the roll, pitch, and yaw. The state constraints are the ones in Romagnoli et al. (2023), while the following input constraints are considered: |F| ≤ 1.5,|τφ,µ,ψ | ≤ 0.043. The cost function is defined by Q=10 I12,R=100 I4, O=106I12, and N=5. The terminal feedback gain, K, is obtained as the solution of the discrete LQR and λ=0.9. The simulation length is of 900 time steps with Ts=0.04s being the sampling period, resulting in a simulation of 36s. The terminal horizon length ˜ M that satisfies condition (26) is ˜ M=92. Tests were carried out considering a linear model of the quadrotor. Fig. 3 shows the resulting position trajectory, together with the corresponding artificial and real setpoints for the position. It can be seen that the quadrotor is able to reach all of them. The real target setpoints have been selected by finding admissible equilibrium points of the system through matrix Mθ. These setpoints only involve non-zero state values in position and in yaw angle due to the dynamics of the system. Again, in Table 1 it can be clearly seen how decreasing λ makes it much more difficult for the system to reach the real setpoints and, therefore, the cumulative cost increases. Regarding computation time, finding ˜ M took 16.9558s, whereas explicitly computing the maximum invariant set required 137.1241s. That is, it was possible to achieve a reduction of 87.63% in the offline computation time, as finding ˜ M reduces to solving simple LPs, while obtaining the explicit invariant set 7 I. Luque, P. Chanfreut, D. Limón et al. Automatica 179 (2025) 112436 with traditional methods requires iterative set operations. The proposed drone system has been chosen as a limit-case example, where the explicit computation is still feasible though computationally expensive. As expected, the online computation time increases slightly with the proposed method due to the extended prediction horizon and the resulting accumulation of constraints. However, in our simulations, this increase was marginal: solving the proposed MPC problem (30) required 0.0164s on average, whereas solving (9) took 0.0114s on average. Finally, notice that, if implemented in a real-world setting, the computation times could be further decreased by employing a faster programming language. 5. Conclusions An MPC formulation for tracking based on Ferramosca et al. (2009), Limón et al. (2008) with implicit terminal components has been presented. It avoids the need for the explicit characterization of the maximal positively invariant set of the system to define the terminal constraint and uses instead an extended prediction horizon. The proposed controller can be efficiently designed and still benefits from the addition of artificial variables that characterize the tracking formulation. In this way, it offers an alternative that can handle larger systems in a tractable manner. It has been shown that the offline cost of the controller design is significantly reduced for large systems, while the online computation times increase slightly. Finally, recursive feasibility and convergence properties have been proved. As a line for future research, the application of the proposed approach to nonlinear systems and more general constraints will be considered, e.g., a natural extension applies to constraints sets that are not only polyhedral but also ellipsoidal, or intersections of both. References Alamo, Teodoro, Cepeda, Alfonso, & Limon, Daniel (2005). Improved computation of ellipsoidal invariant sets for saturated control systems. In IEEE 44th Conference on Decision and Control (pp. 6216–6221). Beard, Randal W. (2008). Quadrotor dynamics and control. Brigham Young University, 19(3), 46–56. Bemporad, Alberto, Casavola, Alessandro, & Mosca, Edoardo (1997). Nonlinear control of constrained linear systems via predictive reference management. IEEE Transactions on Automatic Control, 42(3). Berberich, Julian, Köhler, Johannes, Müller, Matthias A., & Allgöwer, Frank (2021). Data-driven model predictive control with stability and robustness guarantees. IEEE Transactions on Automatic Control, 66(4). Chanfreut, Paula, Maestre, José María, Ferramosca, Antonio, Muros, Francisco Javier, & Camacho, Eduardo F. (2021). Distributed model predictive control for tracking: A coalitional clustering approach. IEEE Transactions on Automatic Control, 67(12), 6873–6880. Comelli, Román, Olaru, Sorin, & Kofman, Ernesto (2024). Inner–outer approximation of robust control invariant sets. Automatica, 159, Article 111350. Ferramosca, Antonio, Limón, Daniel, Alvarado, Ignacio, Alamo, Teodoro, & Camacho, Eduardo F. (2009). MPC for tracking with optimal closed-loop performance. Automatica, 45(8), 1975–1978. Garone, Emanuele, Di Cairano, Stefano, & Kolmanovsky, Ilya (2017). Reference and command governors for systems with constraints: A survey on theory and applications. Automatica, 75, 306–328. Gilbert, Elmer, & Kolmanovsky, Ilya (2002). Nonlinear tracking control in the presence of state and control constraints: a generalized reference governor. Automatica, 38(12), 2063–2073. Gilbert, Elmer G., & Tan, K. Tin (1991). Linear systems with state and control constraints: The theory and application of maximal output admissible sets. IEEE Transactions on Automatic Control. Kerrigan, Eric Colin (2001). Robust constraint satisfaction: invariant sets and predictive control (Ph.D. thesis), University of Cambridge UK. Limón, Daniel, Alvarado, Ignacio, Alamo, Teodoro, & Camacho, Eduardo F. (2008). MPC for tracking piecewise constant references for constrained linear systems. Automatica, 44(9), 2382–2387. Limón, Daniel, Ferramosca, Antonio, Alvarado, Ignacio, & Alamo, Teodoro (2018). Nonlinear MPC for tracking piece-wise constant reference signals. IEEE Transactions on Automatic Control, 63(11), 3735–3750. Löfberg, Johan (2004). YALMIP: A toolbox for modeling and optimization in MATLAB. In IEEE International Conference on Robotics and Automation (pp. 284–289). Luque, I., Chanfreut, P., Limón, D., & Maestre, J. M. (2024). Designing implicit invariant sets for model predictive control for tracking. In IEEE 63rd Conference on Decision and Control (pp. 1795–1800). Mayne, David (2013). An apologia for stabilising terminal conditions in model predictive control. International Journal of Control, 86(11), 2090–2095. Mayne, David Q. (2014). Model predictive control: Recent developments and future promise. Automatica, 50(12), 2967–2986. Morato, Marcelo M., Cunha, Victor M., Santos, Tito L. M., Normey-Rico, Julio E., & Sename, Olivier (2021). Robust nonlinear predictive control through qLPV embedding and zonotope uncertainty propagation. IFAC-PapersOnLine, 54(8), 33–38. Raković, Saša V., & Zhang, Sixing (2022). The implicit maximal positively invariant set. IEEE Transactions on Automatic Control, 68(8), 4738–4753. Raković, Saša V., & Zhang, Sixing (2023). Model predictive control with implicit terminal ingredients. Automatica, 151, Article 110942. Riverso, Stefano, & Ferrari-Trecate, Giancarlo (2012). Tube-based distributed control of linear constrained systems. Automatica, 48(11), 2860–2865. Romagnoli, Raffaele, Krogh, Bruce H., de Niz, Dionisio, Hristozov, Anton D., & Sinopoli, Bruno (2023). Software rejuvenation for safe operation of cyber– physical systems in the presence of run-time cyberattacks. IEEE Transactions on Control Systems Technology, 31(4), 1565–1580. Trodden, Paul (2016). A one-step approach to computing a polytopic robust positively invariant set. IEEE Transactions Automatic Control, 61(12), 4100–4105. Trodden, Paul A., & Maestre, Jose Maria (2017). Distributed predictive control with minimization of mutual disturbances. Automatica, 77, 31–43. Wang, Zheming, & Jungers, Raphaël M. (2020). Scenario-based set invariance verification for black-box nonlinear systems. IEEE Control Systems Letters, 5(1), 193–198. Irene Luque received the B.S. degree in Industrial Engineering in 2020 and the M.S. degree in Robotics and Automation Engineering in 2022, both from the University of Seville, Spain. She is currently pursuing the Ph.D. degree in Automation Engineering at the same institution under the Spanish University Professor Training Program (FPU). She worked for ERC Advanced Grant OCONTSOLAR between 2021 and 2022. Her research focuses on model predictive control and cybersecurity approaches for cyber–physical systems. Paula Chanfreut is an Assistant Professor at the Department of Mechanical Engineering of Eindhoven University of Technology (The Netherlands). She received her Ph.D. degree in Automation Engineering from the University of Seville (Spain) in 2022, where she was a predoctoral fellow under the Spanish University Professor Training Program (FPU). Between 2022 and 2023, she worked for ERC Advanced Grant OCONTSOLAR. Her research is framed within the field of MPC, with emphasis on its noncentralized implementations. Daniel Limón received the M.Eng. and Ph.D. degrees in electrical engineering from the University of Seville, Seville, Spain, in 1996 and 2002, respectively. From 1999 to 2007, he was an Assistant Professor with the Departamento de Ingeniería de Sistemas y Automática, University of Seville, from 2007 to 2017 Associate Professor and since 2017, Full Professor in the same Department. He has been visiting researcher at the University of Cambridge and the Mitsubishi Electric Research Labs in 2016 and 2018 respectively. Dr. Limon has been a Keynote Speaker at the International Workshop on Assessment and Future Directions of Nonlinear Model Predictive Control in 2008 and Semiplenary Lecturer at the IFAC Conference on Nonlinear Model Predictive Control in 2012. He has been the Chair of the fifth IFAC Conference on Nonlinear Model Predictive Control (2015). His current research interests include model predictive control, stability and robustness analysis, tracking control and data-based control with application to efficient operation of buildings and water distribution networks and spacecraft rendevouz strategies. 8 I. Luque, P. Chanfreut, D. Limón et al. Automatica 179 (2025) 112436 José M. Maestre holds a PhD from the University of Seville, where he currently serves as a full professor. He has held positions at TU Delft, the University of Pavia, Kyoto University, and the Tokyo Institute of Technology. He is the author of Service Robotics within the Digital Home (Springer, 2011), A Programar se Aprende Jugando (Paraninfo, 2017), Sistemas de Medida y Regulación (Paraninfo, 2018), and Model Predictive Control (Springer, 2025). He is also the editor of Distributed Model Predictive Control Made Easy (Springer, 2014) and Control Systems Benchmarks (Springer, 2025). His research focuses on the control of distributed cyber-physical systems, with a special emphasis on integrating heterogeneous agents into the control loop. He has published more than 200 journal and conference papers and has led multiple research projects. His achievements have been recognized with several awards and honors, including the Spanish Royal Academy of Engineering’s medal for his contributions to predictive control in large-scale systems and the distinction of becoming the youngest full professor in the Spanish university system in 2020. 9