Stabilizing non-linear model predictive control using linear parameter-varying embeddings and tubes

This paper proposes a model predictive control (MPC) approach for non-linear systems where the non-linear dynamics are embedded inside a linear parameter-varying (LPV) representation. The non-linear MPC problem is therefore replaced by an LPV MPC problem, without using linearization. Compared to general non-linear MPC, advantages of this approach are that it allows for the tractable construction of a terminal set and cost, and that only a single convex program must be solved online. The key idea that enables proving recursive feasibility and stability, is to restrict the state evolution of the non-linear system to a time-varying sequence of state constraint sets. Because in LPV embeddings, there exists a relationship between the scheduling and state variables, these state constraints are used to construct a corresponding future scheduling tube. Compared to non-time-varying state constraints, tighter bounds on the future scheduling trajectories are obtained. Computing a scheduling tube in this setting requires applying a non-linear function to the sequence of constraint sets. Outer approximations of this non-linear projection-based scheduling tube can be found, e.g., via interval analysis. The computational properties of the approach are demonstrated on numerical examples

the terminal set and cost to guarantee recursive feasibility and closed-loop stability are available in the LPV setting.However, because the scheduling variable is considered a free signal that takes its values independently of the other signals in the embedding, some conservatism is introduced.
To give a brief overview of the state-of-the-art based on this concept, the most relevant NMPC approaches based on LPV embeddings are reviewed next.For a comprehensive survey of MPC methods applicable to LPV systems in general, we refer to [4].
In the scheduling quasi-min-max approach of [5], it is assumed that the non-linear system can be represented by an LPV model obtained as a family of linearized models around various operating points.It is shown that if the optimization problem remains feasible, the closed-loop system is asymptotically stable.However, because the family of linearized models is not guaranteed to embed the true non-linear dynamics along closed-loop trajectories, recursive feasibility cannot be established a priori.In the 'gain-scheduling' method of [6], the nonlinear dynamics are embedded in an LPV representation with a discrete scheduling set, effectively representing a collection of uncertain linear models.Rate-of-variation (ROV) bounds on the state and scheduling variable are imposed.
Based on the feedback min-max LPV MPC algorithm of [7] that assumes complete independence of the scheduling variable, the work [8] presents an MPC design applicable to LPV embeddings.The approach relies on the online solutions of linear matrix inequalities (LMIs), and known and bounded ROV on the scheduling variable is assumed.As pointed out in [9], the approach [8] is not guaranteed to satisfy input constraints.
Recently, in [10,11], an 'iterative' MPC scheme was presented to control non-linear systems embedded in an LPV representation.Therein, knowledge of a scheduling map that describes the non-linear relation between the scheduling variable and the state and input signals is exploited in the prediction stage.Based on an initial guess of the future scheduling trajectory, a simple linear time-varying (LTV) MPC problem is solved.The resulting optimal state and input trajectories are used to generate a new predicted scheduling trajectory, and the procedure is iterated until the predicted trajectories converge.This is computationally highly efficient because at each time instant, only the solution of a sequence of LTV MPC sub-problems is required.If the sub-problems converge and the cost decreases, the closed-loop system can be shown to be asymptotically stable.The concept was applied in [12] to the control of a control moment gyroscope using LPV models in both state-space and input-output representations based on velocity-based linearization.A related approach is [13] where an ARX-based least-squares estimate of the future scheduling trajectory is exploited to construct a single LTV MPC problem to be solved online.
In the current paper, an alternative non-linear MPC approach is developed by utilizing LPV embeddings of the non-linear dynamics in a tube-based LPV-MPC framework, as a further development of the ideas first described in [14].The controller requires at each sampling instant the solution of a single linear program (LP) or quadratic program (QP).As a key contribution of the paper, it is shown that by bounding the trajectories of the scheduling variable of the embedding, strong guarantees on recursive feasibility and stability are inherited from the tubebased MPC framework for LPV systems [15].The principal idea that enables bounding the evolution of the scheduling variable, is to first restrict the state evolution of the non-linear system to an initial time-varying sequence of state constraint sets.By applying the known scheduling map to these constraint sets, a corresponding set sequence of possible scheduling trajectories is generated.During the MPC optimization, the scheduling variable is considered as a free signal that can vary inside the constructed scheduling sets.
The approach is not restricted to a specific method of constructing the initial state constraint sequence.In practice, it is proposed to construct it either as a sequence centred around an initially feasible trajectory, or based on a maximum allowed ROV of the state variable.Different ways to compute the corresponding scheduling sets from the selected constraint sets are also proposed: these include straightforward computation of bounds using global optimization or gridding, and a more advanced approach based on interval analysis previously employed in non-linear reachability analysis.
In summary, the usage of an LPV model to embed the nonlinear dynamics allows us to formulate the control problem in terms of a single LP or QP to be solved online with an a priori guarantee on recursive feasibility and closed-loop stability.In return, an additional mechanism to bound the future evolution of the introduced scheduling variable is now required.
Both the controllers presented here and in [10,11] make use of a scheduling map to take into account the relationship between the state and scheduling variables.Where [10] solves a sequence of LTV MPC problems at each time instant, the tube-based method that this paper proposes solves a single LP or QP which is, however, more complex than the problems solved in [10].Therefore, both approaches provide a different trade-off between computational efficiency and a priori stability guarantees.The stability guarantees of the method presented herein follow from a suitable construction of the LPV embedding combined with the properties of the tube-based LPV MPC framework.The guarantees of [10,11] rely on convergence of a sequence of QPs to a local minimum that results in a cost decrease, which is not a priori verifiable in general.
As the final part of this introduction, we elaborate on the differences between the current paper and our earlier work on tube-based MPC approaches for LPV systems [14][15][16].
The paper [15] presents the construction of a terminal set and cost, based on the concept of finite-step contraction, to guarantee stability of a tube model predictive control (TMPC) scheme for LPV systems.Therein, the scheduling variable is considered as an exogenous signal that is not dependent on any state or input variables.As a result, the bounds on the future evolution of the scheduling variable can be arbitrarily fixed.In the current work, the scheduling variable is a signal that depends on the internal state variables of an embedded non-linear system.This means that the bounds on the scheduling variable have to be computed dynamically in order to be consistent with the relationship that exists between the state and scheduling variables.Therefore, the conditions under which recursive feasibility and closed-loop stability can be proven are new and different from [15].
Next, the paper [16] presents alternative tube parameterizations that can be used to enhance the trade-off between computational complexity and control performance.The scheduling variable is considered to be an exogenous signal in the same way as in [15].The development of the parameterizations in [16] was driven by optimizing the trade-off between computational complexity and the volume of the closed-loop region of attraction.In the current paper, the main challenge for the tube construction lies in 'over-approximating' embedded non-linear dynamics by variations of the scheduling variable without introducing too much conservatism.As such, the conditions that the tube must satisfy are different from [16].
The work presented in [14] was an initial development leading up to the work presented in the current paper.The main differences and improvements over [14] are (i) strengthening of all theoretical results through, for example the introduction of the notion of 'convex outer approximator' and a new definition of -contractivity for LPV embeddings, (ii) addition of a section on how to compute the bounds on the scheduling variable, which is critical to practical applicability of the proposed approach, and (iii) new numerical simulations on more challenging examples and a comparison with existing controllers from the literature.

Outline
The remainder of this paper is structured as follows.Preliminaries are discussed in Section 2, which includes notation, an overview of LPV embeddings, and the problem setting considered herein.Section 3 presents the developed algorithm, Section 4 discusses possible initialization procedures to find an initial state constraint sequence, and Section 5 proposes methods to compute the required bounds on the scheduling sets based on the selected state constraints.The approach is demonstrated on two numerical examples (a pendulum driven by a DC motor and a two-tank system) in Section 6.The paper ends with concluding remarks given in Section 7.

PRELIMINARIES
This section introduces notation and basic definitions, gives an overview of LPV embeddings, and introduces the problem setting considered by the authors.

Notation and basic definitions
A set that is convex, compact, and has a non-empty interior that contains the origin is called a PC-set.A polyhedron is a convex set that can be represented as the intersection of finitely many half-spaces.A polytope is a compact polyhedron and can equivalently be described as the convex hull of finitely many vertices.The power set of  ⊆ ℝ n is the set of all subsets of  (including the empty set and  itself), and is denoted by 2  .Sequences are denoted compactly as The Hausdorff distance between a set  ⊆ ℝ n and the origin is where ‖ ⋅ ‖ can be any vector norm on ℝ n .
A function f ∶ ℝ + → ℝ + is of class  if it is continuous, strictly increasing, and is continuous, strictly decreasing, and lim →0 g() = 0. Lastly, a function h ∶ ℝ + × ℝ + → ℝ + is said to be in class  if it is class- in its first argument and class- in its second argument.

LPV embeddings of non-linear systems
This section describes the general concept of an LPV embedding, which will be used throughout this paper.Consider a general constrained system represented by a state-space (SS) representation and an LPV state-space (LPV-SS) representation of the form where k ∈ ℕ.In the above equations u ∶ ℕ →  ⊂ ℝ n u is the control input, x ∶ ℕ →  ⊂ ℝ n x is the state variable, and  ∶ ℕ → Θ ⊂ ℝ n  is called the scheduling signal.Constraints on the input and state variables are represented by the compact sets (, ).
The representation (2) can be constructed such that it is an embedding of the non-linear system (1).The general process of constructing an LPV embedding is depicted in Figure 1.Starting with the original non-linear system, scheduling variables  ∈ Θ ⊂ ℝ n  are introduced such that the non-linear system (1) is equivalently represented by an LPV model as where the dependence of the matrix functions (A, B) ∶ Θ × Θ → ℝ n x ×n x × ℝ n x ×n u can be non-linear.Typically, the introduced scheduling variables are related to the state and input signals, that is there exists a scheduling map  ∶  ×  → Θ such that  = (x, u) (Figure 1.(b)).When an embedding is constructed on a bounded operating region, it is possible to bound the values of  inside of a bounded scheduling set Θ ⊂ ℝ n  .Then, to construct the embedding, the map (x, u) ↦  is discarded and only the inclusion (k) ∈ Θ is assumed (Figure 1.(c)).Following the above high-level explanation, the concept of an LPV embedding that establishes the link between (1) and ( 2) can be formally defined as follows: Definition 1.Let  ⊆  be an open, bounded, and invariant set corresponding to an operating region of the non-linear system (1).Define  = {u ∈  | ∃x ∈  ∶ f (x, u) ∈ }.Let  ∶  ×  → Θ be a scheduling map with Θ ⊂ ℝ n  a bounded set.An LPV-SS representation (2) is an embedding of (1) on , if for all (x, u) ∈  ×  there holds f (x, u) = A((x, u))x + B((x, u))u.From the above definition, it can be concluded that any state trajectory-within the operating domain -of the non-linear system (1) is equal to a state trajectory of the embedding for a particular trajectory of the scheduling signal defined by the scheduling map (⋅, ⋅).However, when a controller is designed for an LPV embedding, it has to guarantee closed-loop stability with respect to all possible scheduling trajectories inside the scheduling set Θ.This leads to the following important result, which the authors will exploit in this article.
Lemma 1 (See [17]).Suppose that there exists a controller that asymptotically stabilizes the origin of an LPV embedding of a non-linear system, for all possible scheduling trajectories inside the bounded scheduling set Θ.Then, the same controller asymptotically stabilizes the origin of the original non-linear system.
The idea of an LPV embedding is reminiscent of that of a linear difference inclusion (LDI).Such inclusions (as well as their continuous-time counterpart, the linear differential inclusion) are well-known tools in the stability analysis of non-linear systems or systems with time delay [18,19].The principal difference between an LPV embedding and a difference inclusion, is the fact that in an embedding the scheduling variable is assumed to be measurable.This measurement can be exploited by a controller, as was done, for example in [20] for flight control design.
In the use of LPV embeddings, there is a trade-off between (computational) simplicity and achievable control performance.On the one hand, using an embedding enables the use of welldeveloped efficient linear design methods to obtain controllers for non-linear systems.On the other hand, the requirement that this linear controller has to be stabilizing with respect to all possible scheduling trajectories can introduce conservatism.
The MPC approach presented herein does not assume one particular method for constructing the LPV embedding of the non-linear system (1), as long as the resulting embedding (2) satisfies Definition 1.However, note that velocity-based embedding approaches cannot be applied directly, as these yield an integral form of the embedding defined in terms of x(k) and u(k).

Problem setting
This paper considers constrained non-linear systems (1) that can be embedded in a constrained LPV-SS representation of the form where the matrix A(⋅) is an affine function given as and where {A i } n  i=1 are real matrices of compatible dimensions.The affine dependence of A(⋅) on  is necessary to enable the derivation of a computationally tractable MPC algorithm later.Note that an embedding with affine scheduling dependence can always be constructed by appropriate choice of the scheduling map (⋅) [26].The following assumptions on the non-linear system (1) are made.
Assumption 1.The system (1) satisfies the following properties: (i) The function f ∶  ×  → ℝ n x is continuously differentiable, and it is bounded on  × .(ii) The origin is an equilibrium of (1), that is f (0, 0) = 0. (iii) The value of the state variable x(k) can be measured for all times k ∈ ℕ. (iv) The state and input constraint sets (, ) are PC-sets.
Knowledge of the scheduling map (⋅) can be exploited to obtain sets of all possible scheduling trajectories based on the possible future state evolution; this is what will be done in the MPC developed herein.Note that by knowing (⋅), Assumption 1.(iii) implies that the value (k) in ( 4) is known for all time k ∈ ℕ as well.
The class of LPV embeddings described by ( 4) is somewhat restricted, because the map (⋅) does not depend on u and because the input matrix B does not depend on .Some remarks on possible future extensions of the results of this paper to a more general setting are provided in the conclusion of Section 7. A class of systems that can always be embedded in the required form of (4) are, for instance, Lur'e systems of the form where (⋅) is Lipschitz continuous and bounded on a compact set  ⊂ ℝ n x .Sufficient conditions under which an embedding of the form (4) can be constructed for an arbitrary continuous-time non-linear system can be found in [3,32].The discrete-time case is still a topic of ongoing research.
The control problem considered here is to design a controller K mpc (⋅, ⋅, ⋅) that renders the origin of the closed-loop representation regionally asymptotically stable in the following sense: Definition 2. Let x(k|x 0 ) denote the solution x(k) of ( 5) for the initial state x(0) = x 0 .The origin is said to be a regionally asymptotically stable equilibrium of (5), if there exists a function  and a PC-set  ⊆  ⊂ ℝ n x such that ‖x(k|x 0 )‖ ≤ (‖x 0 ‖, k) for all x 0 ∈ , and for all k ∈ ℕ.
From the basic properties of embeddings, it follows that if the developed MPC stabilizes the LPV representation (4), then it stabilizes the original non-linear system (1) as well (Lemma 1).It is emphasized that the output of K mpc (⋅, ⋅, ⋅) is the result of an optimization problem that is solved online at every sampling instant.

THE MPC ALGORITHM
In this section, the proposed MPC algorithm for LPV embeddings is presented.

Fundamentals: tubes for embeddings
Before introducing the concept of a tube, the so-called  1property and the controlled image must be defined as mathematical preliminaries.

Definition 4.
The controlled image of a set for a constrained LPV embedding represented by the LPV-SS representation with a given controller The type of tubes we employed can now be defined: where X i ⊆ ℝ n x are sets and where With respect to the tubes that were used in previous work [15,16,33], our approach has the following two important differences: • The state constraints are allowed to be time-dependent In the setting of this paper, this stricter condition is necessary to guarantee consistency between the transition constraints (X i , Θ i |K i ) ⊆ X i+1 and the dynamics (4).
Note that if all involved sets are polytopic, all the conditions appearing in Definition 5 can be expressed in terms of linear inequalities.The role of the scheduling tube Θ ⊆ Θ N is to bound all possible behaviours of  in the embedding (4).To this end, the sequence of time-varying state constraints allowed by Definition 5 is constructed at each time instant k, and an associated scheduling sequence Θ is generated by applying the scheduling map (⋅).So, the signal  is still free, but all its possible future trajectories are known to belong to the scheduling tube T heta in a way that is consistent with all possible future state trajectories allowed by the constraints (6).This setup, that allows constraint sequences consisting of arbitrary sets, is more general than the bounds on the state ROV assumed in [6].
For notational and computational simplicity, this paper considers homothetic tubes [33][34][35]: Definition 6.A tube T satisfying to Definition 5 is called a homothetic tube if (i) For all i ∈ [0..N ] the cross sections are parameterized as X i = z i ⊕  i S with (z i ,  i ) being the centre and scaling, respectively, with S being a fixed polytopic shape set; (ii) For all i ∈ [0..N − 1] the controllers K i ∶ X i × Θ i →  are parameterized as gain-scheduled vertex controllers defined on the corresponding sets X i .
It is, however, possible to use the approach from this paper with the heterogeneously parameterized tubes of [16] or the periodically parameterized tubes of [15] as well.
To obtain a scheduling tube consistent with the state constraints ( 6), the map (⋅) has to be applied to set-valued arguments X ⊆  according to Because (⋅) is a non-linear function in general, Equation 7is not straightforward to compute directly.If X is a polytope, then there is no guarantee that (X ) is also a polytope, or even convex.Therefore, in the sequel, it is assumed that there exists a so-called convex outer-approximator μ(⋅) of (⋅).
In this section, the existence of a convex outer-approximator is simply assumed.Possible constructions of such an approximator are presented later in Section 5.
To derive an MPC approach for LPV embeddings, for a measured initial state x and a given constraint sequence ⇀ , a set of feasible tubes which depends on the time-varying state constraints can be defined as where X f is a terminal set, whose construction is discussed later in Section 3.3.For control, it is useful to select from the set (8) a tube that optimizes a specified performance criterion.This is done by solving the tube synthesis problem where J N (⋅, ⋅) is the cost function, which will be specified later.The function V (⋅|⋅) is called the value function, and an opti-mizer of ( 9) is denoted by Observe that ( 8)-( 9) are not dependent on T heta, as one may expect from Definition 5, but on ⇀ .This is because in the proposed approach, the scheduling sets in the sequence T heta ⊆ Θ N are generated by applying the outer-approximator μ(⋅) to the given sequence This construction is also the reason why the state constraints in Definition 5 have to be of the form X i ⊆  i instead of the more relaxed form (X i , Θ i |K i ) ⊆ .The tube transition constraints are of the form (X i , Θ i |K i ) ⊆ X i+1 where, according to (8), Θ i = μ( i ).If X i would not be a subset of  i , there could exist x ∈ X i such that (x) ∉ Θ i , and therefore the transition constraint (X i , Θ i |K i ) ⊆ X i+1 would not be consistent with the dynamics (4).Therefore the condition X i ⊆  i is necessary.
To realize (9), the cost function is chosen in the form with (⋅, ⋅, ⋅) being the stage cost and F (⋅) being the terminal cost.The design of the stage cost, which specifies a performance target, is discussed next in Section 3.2.The design of the terminal set and cost, essential for proving recursive feasibility and closed-loop stability, will be addressed in Section 3.3.Details on how the tube synthesis problem ( 9) can be formulated as a single LP or QP are given in the Appendix.The exact implementation is dependent on the selected stage cost function, which is addressed next in Section 3.2.

Stage cost design
This section presents two possible designs for the stage cost (⋅, ⋅, ⋅) in ( 10): a worst-case cost used in previous work [15,16,33], and a centre-based cost that is introduced here.Whether to use the worst-case or the centre-based cost is a choice that can be made by the control designer according to the needs of a particular application.The worst-case cost is defined as where ‖ ⋅ ‖ can be any vector norm, c ≥ 1, and (Q, R) ∈ ℝ n q ×n x × ℝ n r ×n x are full column rank matrices corresponding to tuning parameters used to specify performance objectives.Some properties of (11) will be important later to guarantee stability of the MPC algorithm we presented.These are summarized in the following proposition.
The stage cost is homogeneous of degree c in the sense that for all  ∈ ℝ + ,  wc (X , K , Θ) =  c  wc (X , K , Θ).
It can also be of interest to optimize an average-or expected cost instead of the worst-case one.Hence we can also introduce an alternative centre-based stage cost that has properties similar to those of (11).Define r (2) ∶= 2 and r (∞) ∶= 1.This alternative centre-based stage cost function, for sets of the form X = z ⊕ S with S = convh{s 1 , … , sq s }, is then formulated as where p ∈ {2, ∞}, and where the matrices (Q, R) ∈ ℝ n q ×n x × ℝ n r ×n u and the scalar P > 0 are tuning parameters of full column rank.In (12), the vectors xi and θ j represent the vertices of the sets X and Θ that are being passed as arguments of  cb (⋅, ⋅, ⋅).As the centre-based cost ( 12) is defined for sets represented as X = z ⊕ S , it is less general than the worst-case cost (11) which can be applied to arbitrary sets.Nonetheless it has some useful properties: In what follows, let K ∶ X × Θ →  be a  1 -controller.The cost function (12) satisfies the following properties: (i) For p ∈ {2, ∞}, there exists a  ∞ -function  such that Remark 1.It is important to note that the centre-based cost (12) does not satisfy a subset inclusion property like  wc in Lemma 2.(i).This has a consequence for the proof of closedloop stability: in particular, if the cost ( 12) is used, it is necessary to relax the initial condition constraint in the set of feasible tubes (8) from X 0 = {x} to x ∈ X 0 .The full proof of stability, where this point is taken into account, will be presented later for Theorem 1.

Terminal set and cost design
To guarantee recursive feasibility of (9), a terminal set constraint X N ⊆ X f is included in (8).The set X f must be designed to be controlled -contractive for (4).In the case of an LPV embedding, this can be defined in the following sense.
Definition 8. Let X and X both be PC-sets such that X ⊆ X ⊆ .Then, the set X is called controlled contractive for an LPV-SS embedding (4), if there exists a local In light of the above definition, the procedure for computing the terminal set X f consists of first choosing a constraint set X, after which X f can be computed using any of the standard approaches available for LPV models (e.g. the algorithm of [37]).If no X f can be found or if its volume is smaller than desired, the constraint set X can be changed and the procedure is repeated.Observe that X f is always bounded as X f ⊆ X.However, a large X also results in a large scheduling set μ( X) and therefore can also lead to a smaller set X f .
Suppose that X f is -contractive according to Definition 8. Then a suitable terminal cost F (⋅) to be used in (10) can be computed with the same basic structure as the terminal cost used in [15,16].As a preliminary, define a so-called set-gauge function as follows: Next, define the terminal cost as where the maximization is done with respect to sets X = z ⊕ S to be in agreement with Definition 6.The function →  is a local set-induced controller that renders X f -contractive in the sense of Definition 8.Because (⋅, ⋅, ⋅) is convex, (13) reduces to a finite-dimensional optimization problem provided that μ( X) is a polytope.Both possible choices of stage cost (11) or (12) and the terminal cost (13) satisfy conditions that guarantee the existence of  ∞ -bounds on the value function of (9): Lemma 4 ([16, Proposition 19] or [36,Proposition 3.2]).There exist  ∞ -functions v, v such that for all x ∈  and

Main result: The algorithm
The proposed MPC approach is summarized in Algorithm 1, and Theorem 1 gives the main result concerning its properties.

ALGORITHM 1
The MPC algorithm for LPV embeddings Theorem 1. Suppose that μ(⋅) is a convex outer-approximator of the scheduling map (⋅) according to Definition 7. Let X f satisfy Definition 8, let F (⋅) be as in (13), and in Definition 6 set S = X f .Furthermore let Then, the algorithm satisfies the following properties: where x + = A()x + BK 0 (x, ) = A((x))x + BK 0 (x, (x)).
(ii) The state of the controlled system reaches X in N steps or less, that is The origin of ( 1) is regionally asymptotically stabilized in the sense of Definition 2.

Proof of (i).
In what follows, let be the sequence of scheduling sets generated by applying the convex outer-approximator μ(⋅) to the sequence , at the next time instant, the scheduling set sequence becomes where Θ = μ( X) according to Definition 8. (Note that when actually executing Algorithm 1 we additionally have Θ + 0 = { + } ⊆ Θ 1 , but for all following proofs the definition of Θ + given above is sufficient and more convenient notationally.) Next, define  = Ψ X f (X N ) according to (13).It is easily observed that the relationship holds.Therefore, if the worst-case stage cost of ( 11) is employed, a feasible candidate tube can be constructed in the same way as in [15,Proposition 8] as Observe that if the tube T satisfies the state constraints of Definition 5, then T + satisfies the state constraints If the centre-based stage cost ( 12) is used, the set of feasible tubes has a different initial condition constraint (see Remark 1), and it is given by The candidate tube Proof of (ii).Given Property (i), the satisfaction of Property (ii) follows by construction of the constraint set sequences Proof of (iii).This proof is essentially identical to the proof of [15, Theorem 13], but slightly adapted to connect more directly to the setting of the current paper.Substitute the solutions T and T + (see the proof of (i)) in the cost function of the tube synthesis problem and compute the difference between the value functions at time k and time k + 1 to obtain where either X + 0 = {x + } or X + 0 = X 1 depending on the applicable definition of T + (see the proof of (i), also for the definition of , Θ i , and Θ).In both cases, observe that (X + 0 , K 1 , Θ 1 ) ≤ (X 1 , K 1 , Θ 1 ) and therefore where the last inequality follows from the definition of l in (13). Since and therefore where the last inequality follows from Lemma

INITIALIZATION APPROACHES
This section presents two possible methods for determining the initial state constraint sequence ⇀ (0) that is required in Algorithm 1.The first method of Section 4.2 constructs the sets around a given initially feasible 'seed' trajectory, whereas the second method of Section 4.1 is based on restricting the ROV of the state variable.
Note that although this section presents two particular constructions, the feasibility and stability properties established in the previous section are valid for any arbitrary construction of the sequence ⇀ (0).

Initial feasible trajectory
Suppose that a feasible initial state and input trajectory can be computed according to the relationship where In comparison to the bounded ROV initialization discussed next in Section 4.2, this approach is generally more likely to yield a feasible solution to (9), but the assumed availability of the initially feasible trajectory T0 can be considered a disadvantage (note, however, that linearization-based approaches [39,40] also presume the knowledge of such a trajectory).If ⇀ (0) is constructed according to (14) with  = 0, then (9) is always feasible at k = 0. Due to the constraint update Step 6 of Algorithm 1, a value of  = 0 however means that the controller can never deviate from the initial trajectory, thus blocking the feedback action of MPC.On the other hand, a value  that is too large can lead to infeasibility as the corresponding scheduling sets become too large.Therefore, a trade-off for the value of  has to be found.
The trajectory T0 can be obtained by solving a non-linear MPC problem.This could be computationally expensive, but this problem only needs to be solved once at the first sample.An efficient approach that can be used to generate an initial trajectory is [10].The method [10] is not guaranteed to converge to a solution, but if it converges, the result can be used to initialize the MPC we proposed which subsequently guarantees recursive feasibility and stability.

Bounded rate of variation
Let  ∈ ℝ n x + denote a bound on the ROV of the state variable.Then, given a measurement of the initial system state x(0) = x 0 , the initial constraint sequence ⇀ (0) can be computed as with Δ() as in (15).The bounded ROV-construction is similar to what was used in [6].Similarly to the situation in Section 4.1 a value of  needs to be found which is not too restrictive in terms of the allowed state variation, but which also does not lead to too large scheduling sets.Note that if there exists a partitioning of the state vector x = (x 1 , x 2 ) such that (x) = (x 1 ), then it is sufficient to impose a bounded ROV on the x 1 -component only.

BOUNDING THE SCHEDULING TUBE
In this section, methods for constructing outer-approximations μ(X ) are proposed.The easy special case of an affine scheduling map is presented first in Section 5.1.For general nonlinear scheduling maps, some methods to construct the approximations are discussed in Section 5.2.These methods include straightforward computation of bounds using global optimization or gridding, and a more advanced approach based on Taylor expansion and interval analysis.The latter approach was used previously-in a slightly different form-for the purpose of reachability analysis of non-linear systems [41].
Throughout this section, for notational simplicity, all results are given with respect to a general set X ⊆  that is assumed to be convex and compact.When applying these results, use should be made of the state constraint sequence ⇀ (k) in order to obtain tight outer-approximations of the scheduling sets that are valid for all possible state trajectories of the closed-loop system.

Affine scheduling map
Suppose that  ∶  → Θ is affine, that is of the form where  ∈ ℝ n  ×n x and b ∈ ℝ n  .An affine scheduling map will typically result when constructing an embedding of a bilinear system (1).Furthermore suppose that the vertex representation of X is available as Then, the outer-approximation μ(X ) is exact and is given by The assumption that the vertex representation ( 17) is available is not restrictive, because hyperplane and vertex representations of the sets in the sequence ⇀  can be constructed simultaneously when using one of the initialization approaches from Section 4. Hence, no computationally expensive conversion between the two representations has to be performed.

General scheduling map
In the general case, the scheduling map  ∶  → Θ is a nonlinear function of the form where for each i ∈ [1..n  ],  i ∶  → ℝ.For a possibly nonlinear and non-convex function (⋅), obtaining a convex outerapproximator μ(⋅) is not always straightforward.Perhaps the simplest way to get an outer approximation μ(X ) of (X ) is to use gridding or sampling of X .Suppose that a number of M points xi ∈ X , i ∈ [1..M ], has been selected.Then, where the computation of the convex hull only requires the elimination of redundant vertices [42].In low dimensions, gridding can work efficiently.The main caveat is that the grid has to be dense enough in order to guarantee that μ(⋅) is indeed a valid outer-approximation in the sense of Definition 7.
Another conceptually simple way of constructing the convex outer-approximation μ(X ) is to compute the interval hypercube Because the functions  i (⋅) are non-convex in general, computing the minima and maxima in (21) requires the use of an optimization algorithm that is guaranteed to find global optima, for example [43].This optimization has to be done online, but only once when the controller is initialized with the constructed initial constraint sequence.Besides the potentially heavy computational burden of global optimization, a possible downside of this approach is that the restriction to an interval hypercube may lead to a conservative over-approximation.This conservatism may be reduced by also considering rotations, by compressing the scheduling space [44], or via the method of [26].An alternative approach that avoids both the need for global optimization and gridding, can be derived based on the method for non-linear reachability analysis presented in [41].Therein, the idea is to construct outer approximations of non-linear functions applied to sets using a combination of Taylor expansion and interval analysis [45].To apply this approach to the present setting, the following assumption is necessary.

Assumption 2. The function 𝜇(⋅) ∶ 𝕏 → Θ in (19) is twice continuously differentiable.
Let D f (z ) denote the gradient of a function f ∶ ℝ n → ℝ evaluated at the point z with z ∈ ℝ n , and let D 2 f (z ) denote its Hessian matrix evaluated at z.The functions  i (⋅) can be approximated by a first-order Taylor series plus a Lagrange remainder term as where the linearization point x ⋆ ∈ X is arbitrary and where the set of possible values of the Lagrange remainder is given by see [46].The idea is to construct a 'global' bound i for the remainder, such that and, evidently, An expression for the set i can be derived following the same basic approach as in [41, Propositions 1-2].Due to the different problem setting we adopted, the result is summarized for convenience in the following proposition.
Lemma 5 ([36,Proposition 8.2]).An overbound i for the Lagrange remainder ] where e Furthermore, the choice of linearization point is 'optimal' in the sense that it gives the smallest possible interval i .
In Lemma 5, evaluation of the intervals i still requires computing the maximum max x∈X |D 2  i (x)| where D 2  i (⋅) is the second derivative of  i (⋅).In [41], it is proposed to compute this bound using interval analysis [45].The principal idea of interval analysis is to determine, given a (multi-dimensional) interval ] in a computationally efficient manner, that is without using global optimization (the max-and min-operations in the above equation are taken element-wise).A list of how interval approximation for a number of common functions can be computed is given in [47], and a commercially available toolbox implementing computationally efficient interval methods is [48].
It can be expected that this approximation works well for weakly non-linear functions with small second derivatives, such that the linearization error bound from Lemma 5 is small.In [41] it is shown that, if necessary, the linearization error bound can be reduced by splitting the sets X .This approach is not pursued further here.
Based on the discussion above, some other approaches that could be considered in a practical scenario are the following: • Apply interval analysis directly on the non-linear function (19), using, for example the implementation of [48].This yields an interval hypercube of the form ( 21), but depending on the function, it may result in a coarse over-approximation.• Consider the first-order approximation (22).Instead of the interval result in Lemma 5, another method (such as gridding or sampling) to bound the Lagrange remainder could be used.

NUMERICAL EXAMPLES
In this section, the capabilities of the MPC approach we proposed are demonstrated on two numerical examples.

Electrically driven inverted pendulum
In this numerical example, the purpose is to control the angle q of the electrically driven inverted pendulum shown in Figure 2.
Standard gravity g 9.81 ms −2 Define the state vector where q is the angle of the pendulum, ̇q is the angular velocity, and i is the motor current.Then, the system is described by the non-linear differential equation , where u is the input voltage that is applied to the motor.The model parameters are displayed in Table 1.Introduce a scheduling variable The dynamics (25) can now be written in the equivalent LPV embedding form as For control design, ( 26) is discretized using a second-order polynomial approximation [2,Chapter 6].This polynomial approximation is more accurate than the first-order (Euler) approximation that is frequently used.Due to the particular structure of ( 26), the second-order discretization still has an affine dependency on  and a constant B-matrix.The sampling time that is used in this experiment is  = 0.04 s, which leads to a discrete-time embedding with the following numerical values: 1.00 0.04 0.07 0 0.90 2.64 0 −0.02 0.64 The system is subject to the constraints For the simulation, the following settings are used: • The prediction horizon and tuning parameters are set to N = 6, Q = diag{100, 1, 1}, and R = 10.An infinity-norm based worst-case cost ( 11) is used.• A 0.98-contractive terminal set X f is computed using the algorithm of [37], with The set X f is contractive under a robustly stabilizing linear state feedback u = Kx, and has 24 vertices.

𝜋 ∞ ∞
] ⊤ .Because in (27) the scheduling map (⋅) is only dependent on x 1 , it is only necessary to bound the ROV of x 1 .
• The non-linear continuous-time model (25) is used as the real controlled system.
The following controllers are compared: TMPC: The controller from this paper with the infinitynorm cost function (11).qLPV-MPC: The iterative approach of [10].NMPC: A non-linear MPC optimizing a quadratic stage cost, but using the same terminal cost and terminal constraint as the 'TMPC' controller.TMPC Heterogeneous: The tube-based controller from [16] with the same infinity-norm cost function (11).
In [16], there is no mechanism to take into account the relation between the scheduling and state variables.For that reason the static scheduling set Θ full = [min x (x), max x (x)] = [−0.22,1] was used to guar- Simulation result for the electrically driven inverted pendulum antee recursive feasibility.With this scheduling set, no feasible solution was found for the initial state x 0 , so a different initial state x0 = 0.75x 0 was used.Constructing the sequence of scheduling sets based on a bounded ROV was not possible with [16], because this would lead to infeasibility after a few samples due to the neglected relation between scheduling and state variables.The framework of [16] allows considerable freedom in choosing the parameterization of the synthesized state tube.For this example, a 'scenario tube' parameterization was employed for the first 3 prediction time instances, with the subsequent prediction instances using the homothetic tube parameterization that is also presented in the current paper.
The obtained simulation result is shown in Figure 3.All controllers drive the pendulum to its upright unstable equilibrium position x 1 = q = 0 while respecting the constraints.
In terms of the obtained state response, the controllers behave similarly, taking note that the 'TMPC heterogeneous'controller starts from a different initial state.The control inputs applied by both TMPC-based controllers are noticeably more aggressive than with the qLPV and NMPC.This difference is mainly because the tube-based controllers optimize an infinity-norm based cost, whereas the other controllers use a quadratic cost.An example comparing the difference between the infinity-norm and quadratic cost in the TMPC is given in Section 6.2.
The RMS value of the closed-loop state trajectories and the total variance 1 of the control input are summarized in Table 2 to give an idea of the relative closed-loop performance of 1 The RMS value of a signal with a length of n samples is defined as RMS(x) = √  the considered controllers.Note that a direct comparison of these numbers is not necessarily meaningful because all controllers aim to optimize a different cost function under different assumptions.
For the TMPC controller, it was found that the tube synthesis problem is infeasible for the considered initial state if all state variables are allowed to vary arbitrarily fast.Thus, the more complex initialization procedure of Sections 4-5 in this case is really necessary.
The 'TMPC heterogeneous' controller from [16] also reported infeasibility for the considered initial state x 0 .When using a sequence of scheduling sets based on the assumption of a bounded ROV, this controller can give an initially feasible solution only to report infeasibility later when the state has evolved in a way that is no longer consistent with the assumed sequence of scheduling sets.This is the reason that with the 'TMPC heterogeneous' controller from [16], the system could only be stabilized for a smaller set of initial states.The undesirable situation where feasibility is lost after a few samples cannot occur with the TMPC controller developed in the current paper, because the relationship between scheduling and state variables is taken into account explicitly.
Furthermore, it was found that for some combinations of tuning parameters, for example Q = diag{400, 4, 4} and R = 1, the qLPV-MPC diverges after a few samples and fails to stabilize the system.In contrast, the TMPC we proposed provides an a priori stability guarantee and is therefore not susceptible to this type of issue.
To get an idea of the relative computational complexity of the different approaches, the solver computation times per sample are summarized in Table 3.The results were obtained with a 2009 Intel Core i7 CPU at 2.8 GHz, with Gurobi 8.1 being the QP and LP solver, and fmincon from MATLAB 2018b the non-linear solver.

Two-tank system
The two-tank system from [8,49] is depicted in Figure 4.A flow of liquid u with density  is being pumped into the upper tank.This tank has a cross-sectional area (CSA) S 1 .Through a pipe with CSA A 1 , liquid flows out into the lower tank which has a CSA of S 2 .From the lower tank, in turn, liquid flows out through a pipe with CSA A 2 .It is desired to regulate the liquid levels h 1 and h 2 at a given setpoint.For this purpose, the input flow u is available as a control input.The dynamics of the system are described by the non-linear differential equations with the parameters shown in Table 4.
The input flow is subject to a constraint  = {u|0 ≤ u ≤ u} and the liquid levels have to satisfy the bounds described by the state constraint set where the values of the constraints are shown in Table 5.In [8], it is shown that the non-linear dynamics (28) can be embedded in a (continuous-time) LPV-SS representation by introducing the scheduling variables where the scheduling set corresponding to ( 29) is (30) After discretizing the resulting continuous-time embedding using the Euler approach with a sampling time of  = 0.9 s, the discrete-time LPV embedding is obtained where (⋅) is the scheduling map corresponding to (29).The goal in the simulation is to bring the level h 2 of the second tank to a reference value of 115 cm.By solving (28), it is found that this corresponds to the equilibrium state and input values ] , u ss = 1.90.
As was done in [8], the problem of regulating the state to a nonzero equilibrium is converted into a 'stabilization' problem by introducing the translated state and input variables x = x − x ss , ũ = u − u ss , so that regulating x to zero becomes equivalent to regulating x to the desired setpoint.In all of the simulations, the following settings were used: • The tuning parameters are set to Q = diag{0, 2}, R = 0.1, and P = 5 for the centre-based stage cost (12); • For the tube-based MPC, a controlled 0.995-contractive terminal set X f with a fixed representation complexity 6 vertices is computed using the method in [36, Appendix A] with X = ; • The tube-based MPC is initialized according to the bounded ROV approach of Section 4.2, with  = [ ] ⊤ ; • With the scheduling map (⋅) from ( 31), even though it is not an affine function, the scheduling set sequence can be computed by evaluating (⋅) at the vertices of the constraints sets as done in Section 5.1.
For the two-tank system, four different simulations are executed with the following corresponding MPC designs: Casavola: The controller from [8], with N = 2 and ] ⊤ .TMPC 1: The tube-based MPC from this paper with infinity-norm cost function (11), and with N = 4 and ] ⊤ .This simulation allows to compare the behaviour of the TMPC with an infinity-norm cost function to [8], which uses a quadratic cost.TMPC 2: The tube-based MPC from this paper with cost function (11), and with N = 4 and ] ⊤ .This represents an initial state for which [8] is infeasible.For this initial state, the TMPC was also infeasible when an arbitrarily fast scheduling ROV was assumed: thus, the initialization of the state constraint sets with a bounded ROV  was necessary in this case.] ⊤ .This compares the behaviour of the TMPC with a quadratic cost function and a longer horizon to [8].
The suggested correction from [9] to the algorithm [8] was applied, requiring that the initial state is contained inside of an ellipsoidal invariant set.This means that the length of the prediction horizon and the assumed scheduling ROV do not influence the domain of attraction (DOA) of [8].For this reason, and because it does not provide a mechanism to ensure that the assumed scheduling ROV is actually satisfied along realized closed-loop trajectories, the controller [8] is run with an assumed arbitrarily fast scheduling ROV.
The simulation results are depicted in Figure 5.All controllers eventually bring the system to the desired set point, but in slightly different ways.In particular, it can be observed that the TMPC with the infinity-norm based cost function pushes the system closer to its constraints than the controllers with quadratic cost functions.The TMPC with quadratic cost and FIGURE 5 Closed-loop state and input trajectories in the two-tank example FIGURE 6 Closed-loop state trajectories for the two-tank example in a phase plot longer prediction horizon achieves convergence a bit faster than [8], which uses a short horizon.Note that the computational complexity of [8] grows exponentially in the horizon length: therefore, its performance for longer horizons was not investigated.
For the purpose of further illustration, the realized closedloop trajectories are shown in a phase plot in Figure 6.The initial tube for 'TMPC 2' is shown in Figure 7, together with the generated state constraint sets: as expected, every tube cross section X i is contained inside of the corresponding constraint set  i .Similar to Section 6.1, a short summary of generic performance indicators of the closed-loop trajectories is given in Table 6.

CONCLUDING REMARKS
The authors presented a tube-based MPC strategy applicable to the control of non-linear systems embedded in an LPV-SS representation.The principal feature of the approach is the bounding of predicted closed-loop state trajectories by a sequence of state constraints that has to be chosen when initializing the controller.The bounding of state trajectories in this way enables the computation of a corresponding sequence of scheduling sets that are valid along all possible realized closed-loop trajectories, thus allowing the derivation of strong feasibility and stability guarantees.
The determination of a sequence of state constraints and corresponding scheduling sets adds a degree of complexity in the initialization stage of the controller that is not present in approaches such as [5,8], where it is simply assumed that the scheduling sets will be valid along realized closed-loop trajectories.This increased complexity appears, however, to be a price that has to be paid for guaranteed recursive feasibility.
To enlarge the applicability of the approach, and to fully exploit the advantages of the proposed setting, the development of improved methods of selecting the initial state constraint sequences should be considered.Related to this, improved and computationally efficient ways for generating the associated scheduling sets should also be investigated.Extension of the method to LPV embeddings with input non-linearities is also of interest.
Lastly, an extension of the method to the tracking of piecewise-constant reference signals can be considered.This would require the computation of a new state constraint sequence and associated scheduling set sequence every time that the reference changes.Further research is necessary to investigate how this can be done efficiently and how convergence of the state of the controlled non-linear system to the reference can be guaranteed.
where u ( j ,l ) i ∈ ℝ n u are decision variables according to the vertex control parameterization of Definition 6. Vertex representations of the scheduling sets Θ i = convh{ θ1 , … , θq i } are assumed to be available.For the set S from Definition 6, both vertex and hyperplane representations S = convh{s 1  Terminal constraint.The terminal constraint X N ⊆ X f from ( 8) is equivalent to where  ∈ ℝ + is a decision variable and a hyperplane representation X f = {x | H f ≤ 1} is used.
Terminal cost.The decision variable  ∈ ℝ + introduced before corresponds to Ψ X f (X N ).Therefore, F (X N ) =  1− l with l the value corresponding to the maximization in (13).
Full LP with worst-case stage cost.We give the construction for c = ∞ in (11), as used in the numerical examples.The cost (11) can be implemented by introducing non-negative slack variables (, , ).Then,  wc (X i , K i , Θ i ) ≤  i is equivalent to ∀i ∈ [0..N − 1]: .q s ] × [1..q i ] ∶  with d being the collection of all decision variables introduced above.By introducing appropriate slack variables and constraints, a similar construction can be made for the case c = 1.
Full QP with centre-based stage cost.We give the construction for p = 2 in (12), as used in the numerical examples.The total QP realizing (9) then becomes

FIGURE 1
FIGURE 1The construction of an LPV embedding.(a) The original non-linear system.(b) The non-linear system equivalently represented.(c) The scheduling variable is 'disconnected' and becomes a free signal: the LPV embedding is obtained

FIGURE 4
FIGURE 4The two-tank system

TMPC 3 :
The tube-based MPC from this paper with quadratic centre-based cost function, and with N = 10 and x 0 = [ 15 40

FIGURE 7
FIGURE 7Tube T 0 computed at k = 0 in the two-tank simulation 'TMPC 2'

(
j ,l ) i ) subject to (A.1)-(A.6)with d being the collection of all decision variables introduced above.

TABLE 1
Model parameters for Example 2

TABLE 2
Performance summary for the pendulum example.Note that 'TMPC Heterogeneous' used a different initial state

TABLE 3
Summary of solver computation times per sample

TABLE 4
[8]ameters of the two-tank system, the same as used in[8]

TABLE 5
[8]straints of the two-tank system, the same as used in[8]