Tube‐enhanced multi‐stage model predictive control for flexible robust control of constrained linear systems with additive and parametric uncertainties

The trade‐off between optimality and complexity has been one of the most important challenges in the field of robust model predictive control (MPC). To address the challenge, we propose a flexible robust MPC scheme by synergizing the multi‐stage and tube‐based MPC approaches. The key idea is to exploit the nonconservatism of the multi‐stage MPC and the simplicity of the tube‐based MPC. The proposed scheme provides two options for the user to determine the trade‐off depending on the application: the choice of the robust horizon and the classification of the uncertainties. Beyond the robust horizon, the branching of the scenario‐tree employed in multi‐stage MPC is avoided with the help of tubes. The growth of the problem size with respect to the number of uncertainties is reduced by handling small uncertainties via an invariant tube that can be computed offline. This results in linear growth of the problem size beyond the robust horizon and no growth of the problem size concerning small magnitude uncertainties. The proposed approach helps to achieve a desired trade‐off between optimality and complexity compared to existing robust MPC approaches. We show that the proposed approach is robustly asymptotically stable. Its advantages are demonstrated for a CSTR example.

and predicts a single control input at every stage. Because of the lack of recourse, the robustness comes at the cost of a significant loss of performance.
Feedback min-max MPC models the presence of feedback information explicitly in the predictions and thus reduces the conservatism of the open-loop schemes. 2 A general feedback min-max MPC optimizes the worst-case value of the cost function over a sequence of control policies, leading to infinite-dimensional optimization problems. One possibility to formulate a feedback MPC method with a finite-dimensional optimization problem is to consider a tree-structure to represent the evolution of the uncertainty 3 because, for each predicted state at every stage, the possibility is considered to adapt the inputs in the predictions. The tree structure grows exponentially with respect to the length of the prediction horizon, making the approach inapplicable in practice for long prediction horizons.
Other related robust MPC approaches optimize the expected value of the cost function 4 or a weighted sum of all the predicted scenarios, as done in multi-stage MPC. 5 The weights of the multi-stage MPC are tuning parameters that provide additional degrees of freedom to improve the closed-loop performance compared to a feedback min-max MPC scheme.
An alternative to the representation of feedback via a scenario tree consists of restricting the optimization to control policies with a fixed structure, linear policies, 6 or affine policies. [7][8][9][10] Tube-based MPC is one of the most discussed robust MPC approaches in the literature that usually considers an affine parameterization of the feedback policies. 10,11 It was shown in Reference 11 that the problem size can be kept the same as that of nominal MPC if the feedback gain is chosen offline and kept constant in the predictions. However, this comes at the cost of performance loss. Tube-based MPC approaches that relax the structure of the control policy or that predict the tube online (as opposed to an invariant tube) can improve the performance as shown in . The performance advantages come at the cost of an increase in computational complexity with respect to the length of the prediction horizon. To handle parametric uncertainties, tube-based MPC based on Farkas' lemma was proposed. 16,17 The complexity of the approach grows linearly with respect to the length of the prediction horizon. Advanced tube-based schemes such as in References 13,14,16, and 17 use contractive sets for the prediction of tubes online. The number of inequalities and the number of vertices that characterize the tube can increase rapidly with respect to the dimension of the states and this makes the approach difficult to implement for high dimensional systems. If low complexity tubes are employed as proposed in References 18 and 19, the schemes can be highly conservative.
The aim of this article is to propose a novel scheme to achieve the dual goal of low computational cost and low conservatism. To achieve this goal, we combined the multi-stage and the tube-based MPC approaches by the classification of uncertainties in Reference 20. The multi-stage MPC is employed to handle significant uncertainties, and the tube-based MPC is used to handle small magnitude disturbances. We extend the scheme proposed in Reference 20 in this contribution in such a way that the rapid increase in problem complexity in multi-stage MPC is addressed both with respect to the number of uncertainties and with respect to the length of the prediction horizon. The proposed scheme gives the user two options that determine the trade-off depending on the requirements of an application. The two options are the choice of the robust horizon in multi-stage MPC and the classification of the uncertainties. The key aspects of the proposed approach are as follows: 1. The branching of the scenario tree is stopped beyond a certain prediction step called robust horizon. An affine feedback policy is employed beyond the robust horizon with the help of tubes to achieve robust constraint satisfaction and recursive feasibility guarantees. 2. Different scenarios are predicted in the proposed framework for any choice of robust horizon greater than or equal to 1. The resulting increase in the number of degrees of freedom enables the employment of low complexity tubes for high dimensional systems without a significant loss of performance when compared to standard tube-based schemes. 3. In addition, the growth of the problem size with respect to the number of uncertainties is reduced by formulating an invariant tube for small disturbances by making use of the ideas proposed in Reference 20.
We investigate in detail the theoretical properties of the proposed approach and demonstrate that the proposed approach is robustly asymptotically stable. We present the resulting tube-enhanced multi-stage MPC scheme as a convex optimization problem that is solved at every time step. This is achieved by employing the tube-based formulations from References 13 and 16 and the multi-stage formulation from Reference 20. The advantages of the scheme are demonstrated for a CSTR example.

PRELIMINARIES
We study discrete-time linear dynamical systems of the form: where x ∈ R n x represents the state, u ∈ R n u represents the input, w ∈ W ⊂ R n x denotes additive disturbances, the matrix A ∈ R n x ×n x represents the uncertain system matrix, and B ∈ R n x ×n u denotes the uncertain input matrix of the controlled system. The system matrix A and the input matrix B are contained in a convex polytope and can be represented as (A, B) ∈ conv({(A i , B i ), ∀i ∈ Γ p }), where conv() denotes the convex-hull operator and Γ p ∶= {1, … , n p }. We assume that there exists a feedback gain K that is stabilizing for all (A, B) ∈ conv({(A i , B i ), ∀i ∈ Γ p }). W is assumed to be a convex polytope with the origin in its interior and is characterized by n w vertices. The bounds of the additive disturbances can be defined in terms of vertices of the set as W ∶= {w|w ∈ conv({w l , ∀l ∈ Γ w })}, where Γ w = {1, … , n w }. The following definitions of invariant sets adapted from References 21 and 22 will be used throughout this article.

Definition 1.
A set S is said to be robust positively invariant (RPI) for the system x + = (A i + B i K)x + w, if ∀x ∈ S, x + ∈ S, ∀w ∈ W, ∀i ∈ Γ p .

Definition 2.
A set S min is said to be the minimal robust positively invariant set (mRPI) if S min is contained in every closed robust positively invariant set.

Definition 3.
A set S max is said to be the maximal robust positively invariant set (MRPI) if S max contains every closed robust positively invariant set.

Multi-stage MPC
The robustness of multi-stage MPC is achieved by modeling the future evolution of the system by a scenario tree as shown in Figure 1. Each branch of the tree denotes a realization of the uncertainties. Each node denotes a predicted state at the corresponding point in time. If all the extreme values of the uncertainties are realized in the predictions, the predicted states form the convex hull of all the possible trajectories in the future until the end of the prediction horizon. Realizations of the uncertainties that are not extreme can also be included to improve the resulting closed-loop performance. The tree branches until the end of the prediction horizon for each predicted node. The availability of feedback information in the predictions is explicitly modeled in the tree structure without restricting the structure of the feedback policy. This makes F I G U R E 1 Scenario tree representation of the uncertainty evolution for multi-stage MPC for the system with additive and parametric uncertainties. The vertex matrices are given by {(A 1 , B 1 ), (A 2 , B 2 )} and the vertices of the additive disturbances are given by w 1 and w 2 the approach less conservative than those which assume linear or affine feedback policies, but exponentially complex. 23 The optimization problem that is solved at every time step is given as follows: subject to: where the set of all indices (j, k) in the scenario tree is denoted as I and the set of indices occurring from a stage k 1 until a certain stage k 2 is denoted by I ⟦k 1 , k 2 ⟧ , where 0 ≤ k 1 ≤ N p and k 1 ≤ k 2 ≤ N p . Also, the set of indices occurring at a stage at k is denoted by The number of branches at every predicted node is given by n d = n p × n w . The uncertainty realization corresponding to each branch of the scenario tree is denoted by d r , r ∈ {1, … , n d }, where d r represents a particular combined realization of the parametric uncertainties (A i , B i ), and of the additive disturbances w l . That is, Each state x c k+1 predicted at the time step k + 1 is the child node in the scenario tree obtained from the parent node x j k , the input u j k , and the realization of the uncertainty d r . The mapping of the index of a child node from the parent node across each prediction step is given by the relationship c = n d (j − 1) + r. Each node in the scenario tree has a nonnegative weight associated with it j k , ∀(j, k) ∈ I. The weighted sum of the stage costs (x, u) along the prediction horizon and the terminal penalty function V f (x) constitute the overall objective function.
The state and the input bounds, and the bounds on the terminal state are enforced via (2c) using polytopic sets X and U, and X f , respectively.
The control input for a particular node must be the same for all the branches to enforce the causality of the control policy, that is, u j k = u l k if x j k = x l k for all (j, k), (l, k) ∈ I ⟦0, N p −1⟧ . However, the future inputs at different nodes can be different as measurement information will be available at the next stages, that is, u j k can be different from u l k if x j k ≠ x l k for all (j, k), (l, k) ∈ I ⟦0, N p −1⟧ . The optimal input at the fist prediction step u 1 * 0 (x) obtained by solving the optimization problem (2) is applied to the plant. The terminal region X f is chosen as the maximal RPI set for a stabilizing control law K f x. The problem size grows rapidly with respect to the number of uncertainties and the length of the prediction horizon. Therefore we investigate solutions of reduced complexity that approximately realize the performance of the multi-stage scheme. This will be discussed in detail in the rest of this article.
In the linear case considered here, since all the extreme realizations of the uncertainties are used in the predictions, the scenario tree predicts the reachable set of state trajectories. Every node of the scenario tree denotes the vertices of the polytope that forms the reachable set of the system for the predicted control policy.

TUBE-ENHANCED MULTI-STAGE MPC
The problem defined in (2) suffers from rapid growth with respect to the number of realizations of the uncertainties and the length of the prediction horizon. We propose to employ two kinds of tubes to deal with the growth in problem complexity as described below: 1. An invariant tube using an affine feedback policy is employed to handle small-magnitude disturbances.
• The invariant tube is obtained offline and hence, the complexity of the optimization problem does not grow with respect to the number of small disturbances.
• The invariant tube is employed only for small disturbances and hence the method does not introduce a large conservatism.
F I G U R E 2 Scenario tree representation of the evolution of the uncertainties for tube-enhanced multi-stage MPC with robust horizon N r = 2 2. Different tubes for each scenario are introduced to handle the significant uncertainties after a predefined horizon (robust horizon) in the prediction, instead of further branching of the scenario tree (see Figure 2).
• The problem complexity grows only linearly with respect to the prediction horizon beyond the robust horizon.
• The formulation optimizes for different feed-forward terms of the predicted tubes at every stage that belong to different scenarios beyond the robust horizon (in addition to modeling full recourse until the robust horizon) and hence the approach is less conservative when compared to a pure tube-based scheme. Low complexity tubes can also be employed for less conservatism making the approach applicable to high dimensional systems.
The following subsections will elaborate the key points discussed above to obtain an improved trade-off between optimality and complexity.

Handling small disturbances using an invariant tube
Because of the multiplicative nature of the parametric uncertainties, they will have a large influence on the state trajectory far away from the origin. Therefore, we classify all parametric uncertainties as large magnitude uncertainties. The additive disturbances can be large or small depending on the application. To reduce the computational complexity, the set of additive disturbances W is decomposed into two polytopes W and W that contain the origin in their interiors such that The set W denotes large magnitude disturbances and the set W denotes small magnitude uncertainties. It is recommended that the uncertain set W is decomposed such that the number of vertices of the large uncertainties W is small while W is of smaller volume. The model that accounts for large uncertainties is defined as follows: for all i ∈ Γ p and for all w ∈ W, where z ∈ R n x is the state and v ∈ R n u is the input of the model (3). The number of vertices of the additive uncertainty set W is denoted by n w and the set W is defined as W ∶= {w|w ∈ conv({w l , ∀l ∈ Γ w )}, where Γ w = {1, … , n w }. The large uncertainties are considered in the predictions using the multi-stage approach. To handle small disturbance set W, an affine feedback policy is employed as follows: where K inv denotes the feedback gain associated with the invariant tube and is chosen such that the parameter-varying closed-loop system {(A i + B i K inv ), ∀i ∈ Γ p } is asymptotically stable. For the system given in (1), the state x and the control law (4), a set S is defined as small disturbance invariant if (A i + B i K inv )S ⊕ W ⊆ S for all i ∈ Γ p . This disturbance invariant set S can be chosen as the convex RPI over-approximation of the minimal RPI set of the model (3). A convex outer-approximation (S) of the minimal RPI set S can be obtained from the algorithm given in Reference 24. This problem is solved offline and hence does not affect the online computation time of the proposed algorithm. In the implementation section, we propose a novel convex optimization problem to over-approximate the minimal RPI set that is disturbance invariant using a linear programming problem. Since the small magnitude uncertainties are handled using an affine feedback policy, they will not be considered in the online optimization problem and hence do not affect the problem complexity.

Handling long prediction horizons using predicted tubes
For the large uncertainties considered in the scenario tree, the problem complexity increases exponentially with respect to the prediction horizon N p . To reduce the problem complexity, the branching of the tree can be stopped beyond a certain prediction step called the robust horizon N r . Beyond the robust horizon, the affine policies (v + K pred z) are employed to handle all large uncertainties considered in the scenario tree. Here K pred denotes the feedback gain associated with the tubes predicted online. The dynamics of the system beyond N r can be described using the following set recursion: for all i ∈ Γ p and Z denotes the tube of states at the current time step and Z + denotes the tube at the successor time step. The idea is illustrated in Figure 2 for a robust horizon N r = 2. It can be seen that the tubes replace branches beyond N r leading to a linear growth of the complexity of the scenario tree with respect to the prediction horizon. The feed-forward terms associated with the affine control law can be different for different scenarios. This can improve the performance of the controller when compared to the tube-based scheme, 16 where the problem formulation considers one feed-forward term per prediction step. The invariant tube is shown on top of the scenario tree for representational purposes. The invariant tube S is used to obtain a suitable back-off that guarantees satisfaction of the original constraints for all possible values of the small magnitude uncertainties W and does not contribute additional costs to the online optimization problem.
Remark 1. The idea of robust horizon has been already proposed for nonlinear systems in Reference 5 by assuming that the uncertainty remains constant beyond the robust horizon. However, no rigorous study on the recursive feasibility and the stability of the closed-loop system was performed. In Reference 25, it was proposed that a terminal set must be employed at the end of robust horizon to achieve stability but this is clearly restrictive. In this work, we propose a rigorous solution using a robust horizon in the multi-stage MPC framework enhanced by a tube-based formulation. We make use of affine control policies beyond the robust horizon and formulate the problem such that the closed-loop system is stable. The proposed formulation is not as restrictive as in Reference 25, the terminal constraints are enforced at the end of the prediction horizon as in conventional MPC schemes and the proposed scheme does not assume that the uncertainties remain constant beyond N r .

Problem formulation
The optimization problem P N p (x) that is solved at every time step is given by subject to: where There are several differences in both the objective and the constraints in (6) compared to that of the standard multi-stage MPC (2). In the optimization problem (6), the objective can be divided into three parts: a multi-stage part J MS k , a tube-based part J tube k and the terminal penalty part J term k . The multi-stage part of the objective function J MS k is applied until N r − 1, and beyond N r − 1, the tube-based part of the objective function J tube k is applied. As always, the terminal penalty J term k is applied at the last prediction step N p . The multi-stage part of the objective function J MS k is the same as in (2). The tube-based part of the objective function J tube k and the terminal part of the objective function J term k have maximization terms associated with them. For the tube-based part of the optimization problem, the worst-case cost associated with the predicted tubes Z j k are obtained at every prediction step. This helps in establishing the optimal value function as a Lyapunov function for the proposed scheme.
The state inside the tube that maximizes the objective function is defined asz j k and is constrained byz j k ∈ Z j k as given in (6f). The variablesz j k are formulated as decision variables for all (j, k) ∈ I ⟦N r , N p −1⟧ . Equations (6e)-(6g) represent the constraints beyond the robust horizon N r . Equation (6e) guarantees the recursive bounding of the state trajectories for the chosen control law v + K pred z. It can also be seen that the affine term v j k can be freely chosen for all (j, k) This can improve the resulting solution even for a robust horizon of N r = 1 when compared to a pure tube-based scheme. Equations (6f)-(6g) denote the state, input, and the terminal constraints. The set of all x ∈ X for which there exists a feasible feedback policy is denoted as X N p . Equation (6d) is formulated at the robust horizon to establish continuity between the predicted scenarios until the robust horizon and the tubes predicted beyond the robust horizon.
The initial state of the scenario tree z 1 0 is a function of the current state x as defined in (6h) and is a decision variable of the optimization problem (6). The optimization problem is solved with tightened constraints Z = X ⊖ S and V = U ⊖ K inv S in (6c) and (6f). The number of branches at every node is given by n d = n p × n w until the robust horizon N r . The uncertainty realization corresponding to each branch of the scenario tree for the proposed formulation is given by Note that the number of branches can be reduced dramatically if n w ≪ n w when compared to the consideration of all uncertainties in a scenario tree. The control input u applied to the system is given by where v 1 * 0 is the first element of the optimal control input sequence obtained by solving (6). Since the polytope W has smaller number of vertices compared to W, the problem size solved using (6) is reduced. Despite the reduced complexity when compared to a standard multi-stage problem, the proposed controller can often achieve a performance comparable to that of multi-stage MPC for a fraction of the computational complexity by the choice of N r and W.
Remark 2. The proposed scheme is flexible and includes options to further improve it which we do not analyze in this work. Two of the possible modifications are listed below: 1. The feedback gain is denoted as a constant K pred for the predicted tubes. This is done only to simplify the presentation. The feedback gains can be different for different scenarios. The only necessary condition is that the gain must be stabilizing and can be chosen freely for different scenarios to improve the performance of the closed-loop. 2. The number of tubes at the robust horizon is formulated equal to the number of nodes predicted until that stage using the scenario tree in (6). The constraint (6d) represents the continuity equation. However, different nodes can be bundled together in one tube and the number of tubes can be smaller than the number of nodes predicted until N r by modifying the continuity constraint (6d). This can help to reduce problem complexity further.

Stabilizing objective function and choice of weights
Since a persisting disturbance w ∈ W is assumed, convergence to the origin cannot be established. Instead, as described in Reference 26, a robust positively invariant set T will be shown to be asymptotically stable using the proposed robust MPC scheme. To achieve stability, we assume that the terminal set, the proposed stage cost and the terminal penalty function satisfy the following properties: A stage cost with these properties as proposed in Reference 26 is given as follows: where Q and R are positive semi-definite matrices. The terminal penalty function can be simply set to zero, that is, V f (z) = 0, ∀z ∈ R n x . The choice of the stage cost is different from the nominal cost nom (z, v) = ||Qz|| p + ||Rv|| p which is generally used in the case of nominal MPC. The nominal stage cost nom penalizes the distance to the origin and the control effort depending on the choice of tuning matrices Q and R. The stage cost (z, v) penalizes the distance to the set Z f and the deviations from the control law K f z. The terminal gain K f is chosen equal to the gain of the predicted tubes K pred . The gain K f can be chosen freely when there is no tube-based part of the scheme. This is further discussed in Section 4.4 and in Section 5. Each uncertainty realization (branch) has a fixed and bounded positive weight r > 0 associated with it for all r ∈ {1, … , n d }. Appropriate weights can be chosen depending on the applications. However, to establish stability the weights must follow certain rules which are formalized in the following assumption. The weights associated with the particular realization of the uncertainty will be assigned to the nodes that result from them in the predictions. For example, if z 1 1 is realized because of (A 1 , B 1 ) and w 1 , the weight associated with the node z 1 1 will be defined as 1 1 = 1 (the corresponding weight of {(A 1 , B 1 ), w 1 }). The requirement for the choice of weights associated with each node in the scenario tree is given in the following assumption. Assumption 1. The weight of the root node 1 0 is smaller than or equal to the minimum of all the weights (i.e., 1 0 ≤ min{ 0 , 1 , … , n d }) and it must be positive 1 0 > 0. The weights j k associated with the other nodes z j k are equal to the weights associated with the realization of the uncertainty from which they are obtained for all (j, k) ∈ I ⟦1, N r −1⟧ . The weights associated with tubes are chosen as Remark 3. The choice of the weights is important in establishing the stability properties of the proposed approach. The weights affect the objective function and thus the value function. If the weights are chosen as per Assumption 1 stability can be proven, see Lemma 4. In the multi-stage part of the scheme, the weights are chosen the same in each stage for the same realization of the uncertainty. Then feasible values of the stage costs for the succeeding step can be obtained by the convex combination of the stage costs that are realized at the current time step. Since the branching stops beyond the robust horizon, the weights must be updated to account for the receding horizon implementation of the MPC scheme. Hence the weights of the tube-based part of the scheme are employed as Assumption 2. The stage cost (z, v) and the terminal penalty V f (z) are convex and positive definite functions and satisfy the following relationships:

IMPLEMENTATION DETAILS
We employ the tube-based MPC approach based on Farkas' Lemma as proposed in Reference 16 to account for the disturbances after the robust horizon. The proposed approach works for any choice of robust horizon and in principle can be chosen equal to N r = 0 in which case, the scheme reduces to the approach in Reference 16 if where F is a n c × n x dimensional matrix, n c denotes the number of state constraints, and 1 denotes a vector with all elements 1 of an appropriate dimension. The tightened input constraint set V is defined as V ∶= {v|Gv ≤ 1}, where G is a n m × n u matrix and n m denotes the number of input constraints. In the following, we discuss three types of tubes: general complexity tube, homothetic tube and low complexity tube. We show that how the proposed formulation (6) can be implemented as convex optimization problem. In addition, we also discuss the subtleties in the reformulation with respect to the proposed approach.

General complexity tube
The complexity of the tubes that are employed beyond the robust horizon can be fixed and defined as where T ∈ R n r ×n x is fixed for all prediction steps and ∈ R n r is a decision variable that is chosen online. Here n r denotes the number of inequalities that describes the tube Z and it is typically larger than 2n x . The matrix T is chosen such that the set Λ ∶= {z|Tz ≤ 1} is -contractive for ∈ [ min , 1), where min is the joint spectral radius of the closed-loop uncertain system matrices for the chosen feedback gain K pred . The details on algorithm for obtaining the largest -contractive set can be found in Reference 22 (pp. 171-184). One can also use off-the-shelf toolboxes such as Multi-Parametric Toolbox 27 to obtain -contractive sets. Now, let us have a look at the problem formulation (6), where the Equations (6e)-(6f) are in the form of set operations. Farkas' lemma can be employed to convert them from set operations to a set of linear equalities and inequalities. This reformulation will help us formulate the optimization problem as a convex optimization problem that does not require any set operations online. Equation (6e) bounds the error dynamics recursively as where the variables j k , ∀(j, k) ∈ I ⟦N r , N p ⟧ can be formulated as decision variables in the optimization problem. Since an affine control law v + K pred z is employed beyond the robust horizon N r , z + ∈ {(A i + B i K pred )z + B i v} ⊕ W for a given i ∈ Γ p and the tube at the next prediction step that bounds all the trajectories for the predicted input and the arbitrary realizations of the uncertainties for all (j, k) ∈ I ⟦N r , N p −1⟧ is given by To achieve the set recursion, Farkas' lemma is employed and it is given below: 16,22 Lemma 1. Given two nonempty sets there exists a nonnegative matrix P that satisfies the equality PT 1 = T 2 and the inequality P 1 ≤ 2 .
Using Lemma 1, we can employ nonnegative matrices P i , ∀i ∈ Γ p , and holds for all i ∈ Γ p iff: The set (9) is reformulated using linear equalities and inequalities as shown in (10). Similarly, the set operations (6f)-(6g) can be reformulated using Lemma 1. The state and input constraints can be formulated using nonnegative matrices P x , P u as follows: The choice of nonnegative matrices P i , ∀i ∈ Γ p , P x , P u if obtained online, results in a nonconvex optimization problem. Hence as proposed in Reference 16, the nonnegative matrices can be obtained offline such that the equality constraints are satisfied. This reduces the computational load of the online problem and convexifies it. The matrices P i , ∀i ∈ Γ p can be obtained by solving the following problem for a given i: Similar to (12), the matrices P x , P u can be obtained by formulating linear programming problems that satisfy the equality constraints (11a), and (11c). The terminal constraints are implemented as follows: where the terminal set is defined as and The terminal gain K f is chosen as the gain of the predicted gain to simplify the implementation and the discussion that follows. Further discussions on the computation of the terminal set can be found in Section 4.4. The minimal RPI set can be obtained using the methods discussed in Reference 24. Another possible over-approximation of the mRPI set can be obtained as S ∶= {z|Tz ≤ }. Here the matrixT is chosen such that {z|Tz ≤ 1} is -contractive for the chosen feedback gain K inv . The right hand side of the inequality that defines S, can be obtained as follows: The objective can be chosen as 1-norm or ∞-norm so that the optimization problem is an LP. The LP (15) guarantees that the set S ∶= {z|Tz ≤ } is robust positively invariant. The nonnegative matricesP i , for all i ∈ Γ p can be chosen as discussed in (11). SinceT andP i , ∀i ∈ Γ p are fixed, the set S can be a conservative over-approximation. The advantage however is that the LP (15) can be solved much faster compared to the algorithm given in Reference 24.
In addition to the reformulation of the set operations in the constraints, the objective function can be simplified by removing the maximization part with the help of slack variables as shown in References 3 and 26. If the stage cost (7) is used with p = 1, min v∈V maxz (z, v + K predz ) can be obtained as follows: subject to: Note thatz is not known a priori and has been added as a decision variable in (16). The constraint (16c) simplifies to − ≤ Rv ≤ because the terminal gain K f is chosen equal to the gain employed in the predicted tubes (i.e., K f = K pred ). The constraints of the state are infinite dimensional. However, a simplification is possible by reformulating the bounds on the state objective (16b) with the help of a nonnegative matrix P Q as discussed earlier and the constraint (16b) can be satisfied for any givenz ∈ Z if all z ∈ Z satisfy (16b) as proposed in Reference 28. Using Lemma 1, we can define a nonnegative matrix P Q such that the following equations hold to satisfy the constraints (16b) and the constraints (16b)-(16d) can be rewritten as follows for all (j, k) ∈ I ⟦N r , N p −1⟧ .: The equality constraint (17a) can be obtained offline and the inequality constraints (17b) and (17c) can be included as the constraints in the optimization problem online. However, this reformulation of the inequality constraints does not provide a tight upper bound which means that the obtained cost is not the same as the one obtained using the original inner maximization in (36a). The formulation, however, retains the theoretical properties of recursive feasibility and stability of the original formulation (6). Combining all the reformulations, the resulting optimization problem P G N p (x) can be formulated as subject to: (6b), (6c), (6h), (10b), (11b), (11d), (13), (17b), (17c), (17d), Tz , ∀(j, N r ) ∈ I N r . In (18), it can be seen that the set operations are replaced by inequality constraints with the help of the nonnegative matrices P i , ∀i ∈ Γ p , P x , P u , P Q that are obtained offline. The constraints (17b)-(17d) bound the objective from above using the slack variables

Homothetic tube
Since the tubes are characterized using inequalities in the general complexity tube, explicit values of the vertices are not available online. Hence, the tight upper bound for the stage costs of the tubes are difficult to obtain. By formulating the predicted tubes as homothetic tubes (i.e., by fixing the shape of the sets and only varying the scaling variable), the characterizations of the vertices can be obtained as follows: where j k ≥ 0 is a nonnegative scalar and the -contractive set Λ is characterized by the vertices ∨ 1 , ∨ 2 , … , ∨ n v , where n v denotes the number of vertices. The constraints for the homothetic tube can be rewritten as follows using Lemma 1.
Note that the terminal feedback gain is chosen as the gain of the predicted tubes K pred as in the general complexity tube case. The resulting terminal set is defined as The homothetic tube enables us to formulate the tube and obtain a tight upper bound on the extreme stage costs of the tubes as follows: where (j, r) denotes the indices associated with the vertices of the tube Z j k for all (j, k) ∈ I ⟦N r , N p −1⟧ . The optimization formulation in the case of homothetic tubes P H N p (x) results as follows: subject to: (6b), (6c), (6h), (20a)-(20d), (22b)-(22e), Tz j N r ≤ j N r 1, ∀(j, N r ) ∈ I N r . As in the general complexity tube case, the nonnegative matrices P i , ∀i ∈ Γ p , P x , P u are obtained offline as defined in (10) Remark 4. The homothetic tubes can also be implemented using the vertices of the tube as proposed in Reference 13. However, the formulation proposed here using the approach in Reference 16 can be advantageous in terms of reduced computational complexity for a small conservatism. For example, a low complexity tube for a system of n x dimensions require 2n x inequalities only, while 2 n x vertices are required to represent the same tube (an exponential increase against a linear increase). By restricting the vertices to get a tight upper bound on the stage cost and not computing the reachable sets, applications to high dimensional systems can be achieved. This is the motivation for the proposed formulation of the tube-enhanced multi-stage MPC.

Low complexity tube
To reduce the computational complexity of the scheme, low complexity tubes can be employed. A polytopic tube for the low complexity tube for a given (j, k) ∈ I ⟦N r , N p −1⟧ can be defined as follows: where T ∈ R n x ×n x and j k , j k ∈ R n x . Notice that the matrix T is a square matrix with n x rows and columns. Because of the reduced number of inequalities describing the tube, the complexity can be reduced significantly. In addition, there exists efficient implementations of the tubes (See 7,29 ) to improve the online computational time. The low complexity tubes can, however, lead to conservatism of the resulting scheme. In the proposed formulation, because of the presence of more degrees of freedom, the conservatism of low complexity tubes can be mitigated to a certain extent while reducing the computational time. A more detailed discussion is presented in Section 6.

Computation of feedback gains and invariant sets
The feedback gains K inv and K pred should be chosen offline such that the joint spectral radii of the systems {(A i + B i K inv ), i ∈ Γ p } and {(A i + B i K pred ), i ∈ Γ p } are in the unit circle of the complex plane. They can be chosen such that they optimize a performance measure: for example, the feedback gain K inv can be chosen such that the volume of the invariant tube S is as small as possible and the prediction gain K pred can be chosen as the optimal gain for the linear quadratic regulator either in the nominal or in the worst-case. As we point out in Remark 2, instead of choosing a prediction gain, one can choose different gains for different scenarios. However, a systematic way to obtain multiple prediction gains is an issue for further investigations. If N r = N p , the terminal gain K f can be chosen to maximize the volume of the terminal set to improve the volume of the feasible domain for a chosen prediction horizon N p . The implementation of the terminal ingredients is simplified and the complexity is reduced by using (14) and (21) because these formulations can be implemented directly in the optimization problem without computing additional sets. As pointed in References 16 and 29, this simplifies the implementation but can lead to a decrease in the feasible domain (for a fixed prediction horizon). In References 16 and 29, mode 2 dynamics are introduced to find a trade-off between improved feasible domain and computational complexity. Another option is to choose the terminal set as large as possible as described in Reference 29 (pp. 216-219). Since the proposed scheme offers additional feed-forward terms, the proposed implementation of the terminal set is not as restrictive as in the standard tube-based MPC schemes. The disturbance invariant set S can be obtained by solving the optimization problem defined in (15).

RECURSIVE FEASIBILITY AND STABILITY PROPERTIES
We formulate the fundamental assumptions to establish the theoretical properties of the proposed scheme as follows.

Assumption 3. A convex compact disturbance invariant polytopic set S is available for the system (3) if
Assumption 3 is required to obtain a nonempty feasible domain. The proposed formulation (6) offers flexibility to decompose the given uncertainty set W. Only the disturbances considered in W is used to build the set S. Hence, satisfying the assumptions is always possible because W = {0} is always a possible choice. If W = {0}, there is no tightening of the constraints required, but this can result in increased computational complexity. Assumption 4. An RPI polytopic set Z f ⊆ X ⊖ S that contains the origin is available for the system (3) for the feedback gain where Z and Z + represent the employed tubes (that over-approximates the reachable set) in the optimization problem (6).
If N r = N p , there is no tube-based part in the predictions. Hence it is sufficient that the terminal set is robustly invariant with respect to the system (3) for a terminal feedback gain K f . In this case, the tube Z in Assumption 4 is defined as a singleton Z = {z} and Z + = {(A i + B i Kf )z + w l , i ∈ Γ p , l ∈ Γ w }. However , if N r < N p , the tube-based part of the scheme requires that the set recursion (5) employed by the tube is also robustly invariant, that is, it is also necessary that the terminal set in robust positively invariant with the employed tubes that over-approximate the reachable sets of the system at every time step for the feedback gain K pred .
Remark 5. We propose to keep the terminal feedback gain K f the same as that of gain of the predicted tubes. If K f is chosen different from K pred , then the requirement in Assumption 4 should be modified as follows to guarantee recursive feasibility: An RPI polytopic set Z f ⊆ X ⊖ S that contains the origin is available for the system (3) for the control law This leads to additional complexity in the terminal set and it is not clear if it leads to advantages in terms of performance. Keeping K f = K pred simplifies the requirement and is consistent with the tube-based schemes proposed in References 16 and 29.

Lemma 2. Suppose Assumptions 3 and 4 hold and
x ∈ X N p such that P N p (x) (6) has a feasible solution, then P N p (x + ) is feasible for all x + ∈ conv({A i x + B i u}) ⊕ W, ∀i ∈ Γ p if the control input applied to the system (1) follows the control policy u = v + K inv (x − z).
Proof. Let the sequence of optimal control inputs obtained by solving the problem (6) be defined as The root node of the scenario tree is a decision variable as defined in (6g). Let the optimal value be denoted as z = z 1 * 0 . The first element in the input sequence sequence v = v 1 * 0 and the optimal root node z = z 1 * 0 are used in the control law , where x is the current state of the plant. The input u is then applied to the plant. The plant evolves from the current state x to the state x + for the applied input u and the realizations of the uncertainties w ∈ W and the system matrices (A, B) ∈ conv{({A i , B i , ∀i ∈ Γ p })}. x + satisfies the constraints because, the additive disturbances w ∈ W and the vertex matrices {(A i , B i ), ∀i ∈ Γ p } are explicitly considered in the scenario tree in the predictions and the invariant set S accounts for the disturbances w ∈ W. At the next time step, the optimization problem (6) is solved again for the realized state x + . Since S is invariant with respect to the additive disturbances w ∈ W, there exists a z + ∈ conv({z j * 1 , ∀(j, 1) ∈ I 1 }) ⊆ Z = X ⊖ S from Assumption 3. There exists a feasible input sequence for the optimization problem P N p (x + ) that is in the convex hull of inputs predicted in the previous time step for the optimization problem P N p (x) for all prediction steps until N p − 1. For the last prediction step, there exists a control law K f z for all z ∈ Z f from Assumption 4. A feasible input sequence for the next time step can be obtained as the convex combination of the predicted inputs as follows: where j k denote the associated convex weights for all (j, k) ∈ I ( j k ≥ 0, and ∑ n d j=1 j k = 1). Since there exists a feasible root node z + and a feasible input sequence v(x + ), problem P N p (x + ) is feasible for all We now show that the recursve feasibility property is retained if the convex optimization problems (18) and (23) related to the formulation (6) are solved.

Corollary 1. Suppose Assumptions 3 and 4 hold and x
Proof. The constraints of the optimization problems (6) and (18) can be compared one to one. The constraints (6b), (6c), (6h) are retained in optimization problem (18). The remaining constraints are direct results of Lemma 1 which establishes sufficient conditions for the set recursion and guaranteeing that a set is a subset of another. Hence, the feasibility arguments discussed in Lemma 2 directly applies to the formulation (18). Hence, P G the control input applied to the system (1) follows the control policy u = v + K inv (x − z). ▪

N p (x) (23) has a feasible solution for the tube
Proof. The same arguments in Corollary 1 directly apply here as well. ▪

Lemma 3. If Assumptions 1-4 hold, then V
is the optimal value function of (6) with the length of the prediction horizon N p .
Proof. Since the stage costs (z, K f z) = 0 and maxz ∈Z (z, K f z) = 0 in the terminal set and the control law K f z is feasible, the additional prediction step does not add any cost to the optimal value function V * N p +1 (x) and hence V *

Lemma 4. If Assumptions 1-4 hold and X N p is compact, then the optimal value function fulfills the following properties
where c 1 is a positive constant.
Proof. Since x ∈ Z f ⊕ S implies that z 1 0 (x) ∈ Z f is a feasible point, the stage cost (z, K f z) = 0 from Assumption 2 and the control law K f z keeps the state z in the terminal set Z f from Assumption 4. Since V f (z) = 0 and all the stage costs are zero V * N p (x) = 0 for all x ∈ Z f ⊕ S. This proves (27).
From Assumption 2, (z, v) ≥ k|z| Z f and the optimal value function satisfies the property V * In the following, the descent property of the optimal value function will be proven. This property establishes the optimal value function as a Lyapunov function and is an important contribution of this article.
As shown in Lemma 1, the initial state z 1 * 0 (x) and the control policy given in (26) are feasible for all x ∈ X N p . Since z + ∈ conv({z j * 1 , ∀(j, 1) ∈ I}), the optimal value function V * N p (x) and a feasible value function at the next time step (V N p (x + )) F I G U R E 3 Illustration of the descent property of the proposed scheme can be written as: The proof will be accomplished in two parts. Initially the multi-stage part of the value function V N p (x + ) will be compared to the optimal value function V * N p (x) obtained one step before. Because of the receding horizon nature of the MPC, the comparison of the multi-stage part of V N p (x + ) will be performed with the multi-stage part and the one-step tube-based part (J MS * k (x) + J tube * N r (x)) of the optimal value function V * N p (x). It will be shown that The multi-stage part of the value function ∑ N r −1 k=0 J MS k (x + ) can be written as (refer to Figure 3): where z Since the proposed stage cost is convex from Assumption 2, the following inequality holds for the stage cost of the root node of the tree at x + .
Multiplying both sides with 1 0 , we get 1 0 From Assumption 1, 1 0 ≤ min{ 1 , … , p }, so 1 ∑ n d j=1 j 1 = 1, which results in: Similarly the stage costs for the prediction steps until N r − 2 of the value function V N p (x + ) can be compared for all (j, k) ∈ I ⟦0, N r −2⟧ by grouping the nodes that are established for a particular realization of the uncertainty i ∈ Γ p and l ∈ Γ w together. Note that, if N r = 0, there is no multi-stage part in the objective. If N r = 1 only the root node belongs to the multi-stage part and it will be compared with the tube-based part (first-stage) of the optimal value function in the previous time-step V * N p (x) as will be discussed later. If N r = 2, the root node will be compared to the multi-stage part as in (34) and for the comparison of the further stages onwards one can skip to (36). The following relationship for the case is discussed whenever N r > 2 until N r − 2 of the value function V N p (x + ). Since the weights j k associated with each node z j k are the same as the weights associated with a realized uncertainty until N r − 1 of the problem P N p (x + ), we have The stage cost at the prediction step N r − 1 of the problem P N p (x + ) can be compared with the tube-based part of the optimal value function V * N p (x). Again using convexity, the following relationship can be established.
where z j+ N r belongs to the convex hull of the tubes Z j * N r predicted in the problem P N p (x). Note that the stage cost of the tube-based part of the optimal value function V * N p (x) contains the maximization objective, the optimal statesz j * N r give the worst-case cost for all (j, N r ) ∈ I N r . Since the weight associated with J tube * N r (x) is greater than the weights associated with all the realizations of the uncertainty from Assumption 1, we have Substituting (34), (35), (37) in (31), we get The tube-based part of the optimal value function V N p (x + ) can be compared with the tube-based part of V * N p (x). Since the predicted tubes of the problem P N p (x + ) from the prediction step N r until N p − 2 belong to the convex hull of P N p (x), we have The stage costs at the prediction step N p − 1 are ( = 0 by definition. Hence the following relationship holds: From (39) and (41), the value function V N p (x + ) can then be written in terms of V N p (x) as follows: This proves (29) and with this Lemma 3 is established. ▪ Remark 6. The proof of descent (29) is valid for all values of robust horizon in the range N r ∈ [0, N p ]. If N r = 0, then the scheme is simplified into a tube-based MPC scheme enhanced by an invariant tube. If N r = N p , the scheme simplifies into a multi-stage MPC scheme enhanced by an invariant tube. The proof was established for a generic case where both the multi-stage and tube components are present. If one of the components is absent, it can be shown that the proof still holds by removing the corresponding elements in the proof.

Lemma 5. If Assumptions 1-4 hold and X G N p
is compact, the optimal value function V G * N p (x) of the optimization problem (18) fulfills the following properties for the system (1) where c 2 is a positive constant.
Proof. The proof can be found in Appendix A.
where c 3 is a positive constant.
Proof. The proof is given in Appendix B. ▪ The terminal set Z f ⊕ S is asymptotically stable for the proposed scheme as it is shown in the following theorem.

Theorem 1.
Suppose the Assumptions 1-4 are satisfied, then the set Z f ⊕ S is robustly asymptotically stable for the controlled uncertain system given in (1) using the proposed scheme (6) with the implementation (18) or (23).
Proof. Since the optimal value function is established as a Lyapunov function in Lemmas 4-6 with respect to the terminal set Z f , the state z of (3) converges to the terminal set asymptotically. The state of the system (1) satisfies the property x ∈ {z} ⊕ S and converges robustly asymptotically to the set Z f ⊕ S. ▪ Finite time reachability of the terminal set and robust asymptotic stability of the minimal RPI set can be proven for the proposed scheme using a dual mode control policy as proposed in Reference 26. The required conditions are formalized in the following assumption.

Assumption 5.
There exists a dual mode control policy that is employed as follows: where X max is an RPI set for the asymptotically stabilizing control law u = K f x for the system (1) that satisfies the conditions Z f ⊕ S ⊂ X max and X max ⊆ X and v and z are optimal solutions of the proposed scheme (6) with an asymptotically stabilizing feedback gain K inv .

Theorem 2.
Suppose the Assumptions 1-5 are satisfied, the minimal RPI set S min of the system (1) is robustly asymptotically stable for the controlled uncertain system defined in (1) employed using the dual mode control policy (51).
Proof. As Z f ⊕ S ⊂ X max and Z f ⊕ S is robustly asymptotically stable from Theorem 1, the state enters X max in finitely many time steps. Since the control policy is switched to K f x when x ∈ X max and that the control law K f x is asymptotically stabilizing for the system (1), the state x converges to the minimal RPI set asymptotically. Hence for the uncertain system defined in (1)  Proof. This follows directly from the proof of Theorem 2. If the minimal RPI set S min is contained in Z f ⊕ S and S min is robustly asymptotically stable, the state reaches Z f ⊕ S in finite time steps. ▪ In the proposed approach, a multi-stage MPC solution is computed on the scenario tree for the large uncertainties with recourse, that is, a tree of future inputs depends on the realization of the uncertainty. The affine feedback is added "on top" to robustify the solution against the small disturbances. The feedback gain K inv is fixed only for small disturbances and the degrees of freedom are increased using the multi-stage approach for large uncertainties resulting in an improved trade-off between optimality and complexity. Also, if the robust horizon is chosen as N r ≥ 1, we have different feed-forward terms at each stage in the predictions. This results in a scheme with the following advantages when compared to multi-stage MPC and tube-based MPC independently: 1. The growth in problem complexity is reduced when compared to a pure multi-stage approach because the small uncertainties are not considered in the scenario tree. 2. The structurally relaxed recourse which is modeled in the prediction for the realizations of the large uncertainties until robust horizon reduces the conservatism compared to pure tube-based MPC. 3. The choice of a robust horizon on the one hand limits the rapid growth of the scenario tree and on the other hand, provides increased degrees of freedom when compared to a standard tube-based MPC resulting in an improved trade-off. 4. The use of low complexity tubes is possible for less conservatism because of increased degrees of freedom in the form of feed-forward terms beyond the robust horizon. This enables the application of the approach to high dimensional systems.

CASE STUDY
The example considered in this article is a continuous stirred-tank reactor (CSTR) with a reaction scheme that is adapted from Reference 30. Two chemical reactions take place in the reactor: The linearized discrete-time model has the form given in (1), where x = [ΔC a , ΔC b , ΔT R , ΔT J ] T and u = ΔF. ΔC a and ΔC b denotes the deviations of concentration of component A and B in mol/l , ΔT R and ΔT J denote the change in the temperature of the reactor and in the jacket temperature in • C with respect to the equilibrium point. The input is the deviation of the feed from the equilibrium in l/h. The uncertain model is given as:

Details of the simulation study
The tube-based part of the proposed scheme was implemented with three different types of tubes and different simulation studies were performed. The types of the tubes studied are given in Table 1.
The additive disturbances are considered as small disturbances and an offline invariant tube S was obtained. The set S is obtained by using (15) Table 2, the volumes and the computation times of the proposed approach for different robust horizons and the tube-based MPC approach are given. In the following, we investigate all aspects of the proposed approach.

The effect of the invariant tube
If the pure multi-stage MPC is applied, it gives rise to 64 branches per node in the scenario tree resulting in more than 100 million scenarios. The additive disturbances therefore are removed from the multi-stage part and are formulated in the (invariant) tube-based part of the scheme. Hence, W = {0} and W = W for the studied example. This results in four branches at every node which is a dramatic reduction when compared to the 64 branches required in the case of a full scenario tree for all the uncertainties. For comparison purposes, tube-based MPC was implemented for the three different types of tubes considered in the proposed framework. The -contractive set resulted in 32 inequalities for the system (1) compared to 18 inequalities for the system (3). The contractive sets were computed using the multiparametric toolbox. 27 Since the additive disturbances are not considered in (3), the resulting complexity of the -contractive set is different. The 3D projections of the -contractive sets and the small disturbances invariant set S are shown in Figure 4.
If the tube-based MPC 16 is employed without the classification of uncertainties proposed in this article, the tube at every time step is characterized by 32 inequalities in the case of homothetic and general complexity tube tubes and the number of constraints increases with the number of additive and parametric uncertainties. In the studied example, we must consider four vertex matrices and 16 vertices of the additive disturbance set. Hence, to formulate the propagation of tubes, 32 × 4 × 16 = 2048 constraints are required per prediction step. In contrast, for the proposed scheme with robust horizon N r = 0, only 18 × 4 = 72 constraints are required to formulate the propagation of the tubes. This is because the small uncertainties are not considered both in the computation of contractive sets and in the online problem. Instead, a suitable back-off is obtained by making use of the disturbance invariant set S. If the low complexity tube is employed, the complexity of the tube-based MPC scheme without the invariant tube is 16 times larger than the proposed scheme with N r = 0 because of the 16 vertices of the additive disturbance bounds. The feasible domains of the tube-based MPC scheme (without the invariant tube) and the proposed tube-enhanced multi-stage MPC scheme with N r = 0 obtained using the general complexity tube predicted tubes are shown in Figure 5. It can be seen that the feasible region of the tube-based MPC scheme that does not employ the invariant tube is larger than the feasible region of the proposed scheme with N r = 0. The volume of the feasible domain of tube-based MPC scheme (without the invariant tube) is 1197.1 whereas the proposed scheme with N r = 0 results in a volume of 1110.7 (approx. 7% smaller). However, the computation time of the proposed approach with N r = 0 is only 0.15 seconds compared to the other scheme which has a computation time of 6.55 seconds. If the homothetic tubes are employed, the feasible domain of the proposed scheme is 6% smaller while the computation time of the proposed approach is approximately 96% smaller (for exact values refer to Table 2). In the case of low complexity tube, there is no difference in the volumes of feasible domains observed.
With this example, it can be seen that the invariant set can introduce a certain conservatism in the closed-loop. However, it reduces significantly the computational complexity and improves the computation time of the approach largely (up to 98% reduction). Thanks to the use an invariant tube for small uncertainties, an important reduction in computation time can be expected at the cost of only minor additional conservativeness.

The effect of the robust horizon
If the proposed tube-enhanced multi-stage approach is employed to handle large uncertainties, the number of scenarios considered by the optimization problem increases with the length of robust horizon. Total number of scenarios in the problem is determined by N r . For N r = 1, the problem has only four scenarios and for N r = 5, the problem has 4 5 = 1024 scenarios. First, we discuss the effect of the types of the predicted tubes that are employed and then summarize the observations.

General complexity tube
The scheme with a robust horizon N r = 4 has the same feasible domain as the scheme with full robust horizon. When the robust horizon decreases further, the volume of the feasible domain decreases monotonically. The comparison of the feasible domains of the proposed scheme with different robust horizons that uses general complexity tubes are shown in Figure 6. The proposed scheme with N r = 1 has a volume that is approximately 12% smaller than the volume of the full horizon case. The volumes of the feasible domains of the tube-based MPC implemented without the invariant tube and the proposed scheme with N r = 0 are significantly smaller (approx. 74% smaller and 76% smaller). In the tube-based MPC scheme, the optimization problem has one feed-forward term per prediction step as degrees of freedom. Whereas, in the cased of N r = 1, there are four feed-forward terms optimized at every stage. This improves the degrees of freedom of the controller and results in an improved performance. The polytopes are plotted with the help of the Multi-Parametric Toolbox. 27 The computation times of the scheme with different robust horizons however does not show a uniform trend. The computation times of robust horizons N r = 0 and N r = 1 are smaller than the full robust horizon. However, the schemes with the robust horizons 2 to 4 have computation times larger than the scheme with full robust horizon. This is because of the difference in complexity associated with the tree and the tube. The scheme with N r = 2 has 16 scenarios. However, from the second prediction step, the propagation of tubes is characterized by 72 inequalities and this is formulated for all 16 scenarios requiring 72 × 16 = 1152 inequalities after the second prediction step. This leads to an increased computational effort when compared to a full tree with 4 × 4 = 16 equality constraints per node (though exponentially increasing every stage).

Homothetic tube
The trends in the volumes of the feasible domain and the computation times are similar to the general complexity tube case, if the homothetic tubes are employed in the predictions. However, the scheme results in an increased conservatism compared to the general complexity tube case. This is expected because the shape of the predicted sets is restricted. The proposed scheme with N r = 1 has a volume of the feasible domain that is smaller by approximately 17% in this case compared to the full robust horizon case. The feasible domains of the proposed scheme implemented with homothetic tubes are given in Figure 7. There is also an increase in computation times when compared to the general complexity tube case for a fixed horizon. This is because of the increase in the number of constraints due to 44 vertices of the tube considered to obtain the tight upper bound of the stage cost. This leads to a proportional increase in computational cost.

Low complexity tube
The computation times are smaller in this case as expected because the tube is represented using the minimal number of inequalities. There is however large conservatism as a result. When N r = 0 is employed, the volume of the feasible region is only 96.02 which is more than 10 times smaller than for the homothetic tube case. Similar reductions in the feasible F I G U R E 6 Comparison of the feasible domains obtained using the proposed scheme that uses general complexity tubes for varying robust horizons [Colour figure can be viewed at wileyonlinelibrary.com] domains can be observed across all robust horizons. The feasible domains of the proposed scheme implemented with the Low complexity tube are given in Figure 8. An interesting point to note here is that when N r = 1 is applied, the volume of the feasible domain is much smaller for the same robust horizon than when general complexity/homothetic tubes are employed. Despite this reduction, the volume is larger than the best feasible domain observed for the tube-based MPC. This clearly demonstrates the advantages of the proposed approach. For robust horizon N r = 2, the feasible domain is larger and the computational cost is still lower than tube-based MPC. Hence the scheme can be used to control high dimensional systems, where the pure tube-based scheme is either highly conservative (in the case of using a Low complexity tube) or intractable (in the case of a higher complexity tube). The proposed scheme offers an alternative with an improved trade-off between optimality and complexity. The computation time and the feasible domain for the case N r = N p = 5 does not depend on the type of tube employed because there is no tube-based part in the predictions of the optimization problems (18) and (23) for this case.

Summary
Below are the summary of the observations: 1. The volumes of the feasible domains of the proposed scheme increases monotonically with the increase in the robust horizon. For a given robust horizon, the volume of the feasible domain of the proposed scheme that employs the Low complexity tube is lower and that of the general complexity tube is higher. The volume of the scheme that employs a homothetic tube is in between but it is closer to the general complexity tube than to the Low complexity tube. 2. The proposed scheme with N r = 1 has a larger volume of the feasible domain and a smaller computational effort than the tube-based MPC that is employed without the invariant tube. Even the proposed scheme with N r = 1 that employs the low complexity tube has a larger volume than the tube-based MPC that employs a general complexity tube. 3. Though the proposed scheme with Low complexity tube shows conservatism with respect to the scheme with full robust horizon, it gives the best computation times. In addition, the proposed scheme employed with Low complexity tubes with N r ≥ 1 shows better performance than the tube-based MPC scheme with complex/general complexity tubes. 4. The scheme offers flexibility with respect to the choice of tube and robust horizon. For example, the proposed scheme that employs a homothetic tube with N r = 1 has a feasible volume and computation times comparable to that of the proposed scheme with Low complexity tube with N r = 2. The scheme offers a wide variety of options to achieve a desired trade-off between optimality and complexity than the existing robust MPC schemes.
For this example, N r = 1 using a homothetic/general complexity tube, N r = 2 using Low complexity tube and N r = 5 are possible choices to obtain good trade-offs overall in optimality and computational complexity. The scheme with N r = 5 gives the best performance for a reasonable computational time. Though the volume of the feasible domain of N r = 1 using general complexity tube is approximately 12% smaller, the computation time is reduced by 62%. Hence N r = 1 using general complexity/homothetic tube is a good choice for this example. The closed-loop state and input trajectories of the proposed scheme employed using a general complexity tube with N r = 1 for random initial conditions and random realizations of the uncertainties for 100 simulation runs is shown in Figure 9.

The effect of the prediction horizon
The growth of the problem complexity with respect to the prediction horizon is analyzed by comparing the average computation times of the schemes with different robust horizons. The results are plotted in Figure 10. It can be seen that the computational cost increases exponentially if a full robust horizon is used and is significantly larger than for the other robust horizons considered for N p = 8. The computation times of the schemes with N r < N p grows linearly in complexity with respect to the prediction horizon. However, the slope is seen increasing when the robust horizon increases. The proposed scheme with robust horizon N r = 1 has a computation time of less than one second and the scheme with N r = 2 has a computation time of less than 3 seconds. The variations in the times are much smaller than in the full robust horizon case. As the horizon grows larger, the computational advantages of the proposed scheme increase.

CONCLUSION
In this article, we have shown that the combination of multi-stage and tube-based MPC schemes offers a flexible framework to manage the trade-off between the performance and computational complexity of the robust scheme. The proposed method uses a tube-based method to handle uncertainties that are small or occur far in the prediction (after the robust horizon), while the multi-stage approach handles the significant and immediate uncertainties to increase the performance. The stability of the scheme for linear systems with parametric and additive disturbances was demonstrated for any choice of robust horizon, including the pure multi-stage case. Simulation results show that the proposed method provides flexibility to obtain a good trade-off between complexity and performance.

ACKNOWLEDGMENTS
The research leading to these results has received funding from the European Commission under grant agreement number 291458 (MOBOCON). RP acknowledges the contribution of the Slovak Research and Development Agency under the project APVV 15-0007.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.