Tube Model Predictive Control with an Auxiliary Sliding Mode Controller

This paper studies Tube Model Predictive Control (MPC) with a Sliding Mode Controller (SMC) as an auxiliary controller. It is shown how to calculate the tube widths under SMC control, and thus how much the constraints of the nominal MPC have to be tightened in order to achieve robust stability and constraint fulfillment. The analysis avoids the assumption of infinitely fast switching in the SMC controller.


Introduction
Model Predictive Control (MPC) has been a great industrial success, particularly in the process industries (Qin and Badgwell, 2003). Still, robustness of MPC controllers continues to be an active research issue. There are many approaches to robust MPC, including dynamic programming, optimization over feedback policies, min-max MPC, etc. Providing an overview over these various approaches to robust MPC is beyond the scope of this paper, the interested reader is instead referred to the bibliographic notes in Chapter 3 of Rawlings and Mayne (2009). As is noted in (Rawlings and Mayne, 2009), optimal or close-tooptimal approaches to robust MPC generally have prohibitive online calculation requirements. More practical approaches are therefore a compromise between performance and online calculation requirements.
A popular approach among such MPC formulations that trade off performance against online calculation requirements is the so-called Tube MPC (Rawlings and Mayne, 2009). In the 'basic' Tube MPC, the MPC controller essentially controls the nominal plant, while there is an auxiliary controller that keeps 'all possible' plants inside a 'tube' close to the nominal plant. To ensure robust constraint satisfaction, the state/output constraints of the nominal MPC have to be modified to account for the width of the tube, and the magnitudes of control inputs that are available to the nominal MPC also have to be reduced to account for the additional input component coming from the auxiliary controller.
It is clearly desirable that the auxiliary controller should add little to the overall computational requirements of the control. If this is the case, the auxiliary controller may execute at a significantly higher sampling frequency than the nominal MPC, thereby taking advantage of feedback also between the samples of the nominal MPC. There is substantial ongoing research effort aiming to reduce the online calculation requirements of MPC. However, approaches like explicit MPC (Bemporad et al., 2002) appear to be limited to a fairly modest number of states, and approaches focusing on online solution of an optimization problem will inevitably have a significant calculation requirement. The higher calculation requirements for MPC will typically result in modest sampling frequencies compared to what is achievable with linear state feedback or controllers of similar complexity. The computational requirements of the Sliding Mode Control (SMC) are trivial, and combined with the good robustness properties of SMC this should make SMC an attractive candidate for the auxiliary controller in Tube MPC.
SMC is known for its robustness to parameter variations and external disturbances. It belongs to a special class of nonlinear discontinuous control algorithms, known as variable structure control (Utkin, 1978;Young et al., 1999;Yu and Kaynak, 2009). In its basic form, the SMC input is simply a relay output, depending on the location of the present state relative to a switching surface, which forces a system state to move along this surface, known also as a sliding surface. The modern realization of SMC, by using microcontrollers or digital signal processors, causes a quasi-sliding motion (Milosavljević, 1985) usually in O(T) vicinity of sliding surface, with T denoting a sampling time, which could induce a chattering, manifested by high-frequency control signal exciting unmodelled system dynamics and reducing plant lifecycle. An overview of existing digital SMC algorithms is given in (Milosavljević, 2004).
In recent years, researchers have developed several control methodologies based on the combination of SMC and MPC. The combination of SMC and generalized predictive control (GPC), as a subclass of MPC, is discussed in (Corradini and Orlando, 1997;Mitić et al., 2013) for systems with discrete-time transfer function models. In (Garcia-Gabin et al., 2009) the cost function is partially optimized with respect to only the predictive part of controller, while sliding mode control is not involved in the optimization problem. Moreover, the reaching and existence conditions of sliding mode are not derived and the stability issues are not discussed. Unfortunately, all these approaches cannot deal with MIMO systems.
In digital SMC based on state-space models, and hence applicable to MIMO systems, one approach to control law design is to force the system to reach the sliding surface at the very next sampling instant (Su et al., 2000). This digital SMC method provides a O(T 2 ) sliding mode accuracy when the system disturbances are known at each time instant or a unit-step delayed disturbance estimate is used in the presence of unmodeled disturbances. However, when the disturbance depends on the control input, the system becomes unstable or causes chattering. That is why two different approaches in integrating digital SMC and MPC are proposed in (Neelakantan, 2005). The first one applies direct optimization of a cost function criterion with respect to the equivalent control. The second control method splits the controller into the equivalent control part, ensuring the system to stay on sliding sur-face once reached, and the reaching control part that guides the system towards the sliding surface. The cost function is optimized with respect to the latter control term.
In (Incremona et al., 2015;Benattia et al., 2015) hierarchical control schemes, consisting of a high level MPC and a low level SMC (Rubagotti et al., 2011), are considered. The role of the SMC component is to reject the matched disturbances acting the plant, and to reduce uncertainty for the MPC design in that way.
The accurate calculation of the tube widths is complicated even when using linear state feedback. A set of states needs to be calculated, inside which the auxiliary controller is able to keep the states of the real system. Once this set is calculated, the corresponding set inside which the input from the auxiliary control will remain has to be calculated. To avoid restricting the nominal MPC unnecessarily, the calculated set should be as small as possible, i.e, we wish to calculate the socalled minimal robustly positive invariant set (mRPI) for the system under the auxiliary control. Unfortunately, the mRPI will typically be excessively complex and demanding to calculate, but good outer approximations can usually be found when linear state feedback is used as an auxiliary controller (Raković et al., 2005).
To the authors' knowledge, there is no previous work on calculating robustly positively invariant sets for systems under SMC control with a finite sampling frequency for the SMC. However, we will also take an alternative approach: following an idea in Rawlings and Mayne (2009), the robustly positively invariant set will not be calculated directly. Instead, we will show how to calculate, for each constraint, how far in the direction of the constraint the true system can be driven by the model uncertainty. This will give a direct measure of how much each constraint will need to be tightened.
Tube MPC was originally proposed with the MPC handling the nominal system only, and the auxiliary controller handling the deviations from nominal behavior. Alternative formulations have been proposed later, where feedback is introduced also into the MPC part of Tube MPC. In this paper, including the example, we have chosen to use a nominal Tube MPC as in the original Tube MPC formulation, in order to highlight the robustness improvement from the auxiliary controller. However, the SMC-based auxiliary controller proposed in the paper would be equally applicable to a Tube MPC formulation where feedback is used also in the MPC part of the Tube MPC.
The paper is organized as follows. In Section 2, the control problem is introduced. Section 3 briefly describes two proposed digital sliding mode control approaches. The algorithms for calculating the con-straints tightening is presented in Section 4. The proposed Tube MPC with an auxiliary SMC has been applied to the real DC servo system (Inteco, 2011), and the digital simulation and experimental results are given in Section 5. Section 6 contains some concluding remarks.

Problem description
Consider the discrete time system described by the model with the system state x ∈ R nx , the input u ∈ R nu , and the disturbance w ∈ R nw . There are also constraints on the allowable state constraints on the allowable input and constraints on the possible range of disturbances It follows from the description above that the sets X,U, and W are polyhedral; we will also assume that they are bounded (and thus that the sets are polytopic), of full dimension, and contain the origin in their interior. In the following, the system state is split into two components, a nominal component z and a deviation from nominal Similarly, the input is split into the input from the nominal MPC v, and the input ν from the auxiliary controller The dynamics, described by eq. (1), may therefore be split into the nominal dynamics and the deviation from nominal Clearly, eqs. (7) and (8) add to eq. (1). The control scheme is illustrated in Figure 1.
At each timestep, the nominal MPC solves the problem subject to The vector δ z i quantifies how much the nominal state constraints have to be restricted at time k + i in order to ensure that the true state adheres to the original constraints, while the vector δ v i quantifies how much the constraints on the nominal MPC input have to be restricted at time k + i in order to ensure that the total input adheres to the original constraints. The simplest Tube MPC formulations treat δ z i and δ v i as constants over the prediction horizon, whereas other formulations allow these to vary to account for the fact that the disturbance will typically be able to drive the true state further from the nominal state (also under the action of the auxiliary control) over a time period of several timesteps than over a single timestep.
Remark. A special terminal set for the state is a common ingredient in MPC formulations guaranteeing closed loop stability. Such a terminal set is ignored in eq. (9) for reasons of notational simplicity, but it adding such a terminal set would be straight forward.
From eq. (9), it is clear that δ z i and δ v i have to be found in order to be able to formulate the nominal MPC. This will be addressed in Section 4.

Digital Sliding Mode Control
To design the auxiliary digital SMC, we consider the deviated system dynamics described by eq. (8). Two control algorithms are implemented herein. The first one is a traditional relay based sliding mode control defined by denotes switching function and is the equation for the sliding surface or the intersection of sliding surfaces if n u > 1 . Notice that K k is usually selected as an auxiliary control in Tube MPC and a matrix K has dimension n u × n x . Here sign(g k ) is understood to be a vector with elements ±1, and ∆ u is a diagonal matrix with constants representing the relay outputs. The sliding surface, i.e. K, should be selected so that the system (Furuta, 1990) is stable. Eqs. (13) and (14), describing system dynamics in sliding mode, are obtained by implementing the well-known equivalent control in eq. (8). Substituting eq. (10) in eqs. (8) and (11) yields defining the switching function dynamics at time instant k, whereas in the prediction horizon it is determined by In order to provide stable switching function dynamics, ∆ u should be calculated according to the following theorem.
Theorem 3.1. If ∆ u is chosen to satisfy the following inequality where Ω is a positive real vector and 1 is a vector of 1's, then, for every initial state g k , there exists a positive integer number k 0 = k 0 (g k ) < N , such that the system phase trajectory, described by eqs. (17) and (18), enters the domain defined by after k 0 time steps and remains in this domain for all i > k 0 .
Proof. See Appendix A.
The second auxiliary digital SMC, used in this paper, is so-called robust discrete-time chattering free sliding mode control (Golo and Milosavljević, 2000) Implementing eq. (20) in eq. (8), the switching function dynamics becomes (21) and, inside the prediction horizon, we have The next theorem gives sufficient conditions for stable sliding motion in prediction horizon.
Theorem 3.2. The system phase trajectory, described by eqs. (22) and (18), reaches the domain G defined by eq. (19) in k 0 = k 0 (g k ) < N time steps for every initial g k , and remains in it for all i > k 0 .
Proof. See Appendix B.

Calculating the required constraint tightening with SMC-based auxiliary controllers
We will first describe how to calculate the required constraint tightening with control defined by eq. (10). However, relay feedback is known to often result in very fast switching, which for some applications will not be desirable. A common remedy is then to replace the 'infinite gain' at the switching surface with a steep linear function, leading to a chattering free SMC described by eq. (20). The second subsection will address the constraint tightening for this type of auxiliary controller. It is noted above that an auxiliary controller with low calculation requirements may operate at a higher sampling rate (shorter timestep) than the MPC. However, we will use the same sampling frequency both for MPC and SMC.
With the simple SMC-inspired auxiliary controllers considered here, the determination of δ v i is trivial, which will become apparent below. However, the calculation of δ z i is more challenging.

Constraint tightening for traditional SMC
Denote the relay term in eq. (10) as To proceed, we define for each element ϑ j of ϑ a binary variable s j , such that ϑ j > 0 ⇒ s j = 1, and s j = 0 otherwise. Also needed are upper and lower bounds on each component of the vector K . From eqs.
(2) and (5), it is clearly safe to assume F ≤ f , and the lower bounds m j and upper bounds µ j on element j of K can be found from the LPs and where K j is row j of K. Let 1 denote a vector of 1's. Equation (23) is then implied by Mignone (2001) ϑ where the value of the binary variables s follow from the constraints Note that numerical optimization solvers cannot distinguish between strict and non-strict inequalities. The formulation above will leave the value of s j undecided if K j = 0. It will then be left for the optimization routine to choose the optimal value.
The furthest from the origin the disturbance sequence {w k } may drive the deviation state k in the direction of the state constraint F l x k ≤ f l over a horizon of N timesteps can then be found by solving subject to For each state constraint j, this optimization should be solved for a number of horizon lengths N . For each timestep i, the elements of δ z i are given by δf j,i . If the system under the relay feedback is stable, δf j,N will approach an upper bound as N grows large.

Constraint tightening for chattering free SMC
Instead of eq. (23) we now use where ∆ uj denotes the j'th element on the main diagonal of the diagonal matrix ∆ u , and K j as before refers to row j of K. We rewrite eq. (30) as To capture this behavior, we need two binary variables, s j and t j for each auxiliary input ϑ j , such that Define where m j is calculated as in eq. (24). We note that q j is non-negative in the domain of interest. The actual input from the auxiliary controller may then be calculated from the expression (37) where we note that s j ≥ t j . The difficulty in the above equation lies in the bilinear terms s j q j and t j q j , both being the product of a binary variable and a nonnegative real. To proceed, we introduce the auxiliary variables σ j = s j q j and τ j = t j q j . From (Bemporad and Morari, 1999), we have that the set can equivalently be expressed as and similarly for (q j , t j , τ j ). Recognizing that in this case a j = µ j − m j , and introducing m 1j = m j + ∆ uj (40) Defining the diagonal matrices Λ = diag(a j ), M 1 = diag(m 1j ),M 1 = diag(µ 1j ), M 2 = diag(m 2j ), and M 2 = diag(µ 2j ), and forming the column vectors m = vec(m j ), s k = vec(s kj ), t k = vec(t kj ), σ k = vec(σ kj ), and τ k = vec(τ kj ), we obtain the optimization formulation subject to (45) K k + Λs k − Λ1 − m < σ k ; k = 0, . . . , N − 1 (60) Clearly, with ν k as in eqs. (10) and (20) taking into account F ≤ f as specified above, we have

Digital simulation and experimental results
The validation of the proposed control methods is performed by using the modular servo system (Inteco, 2011) shown in Figure 2. The objective is to control the angular position of the DC motor shaft. The system consists of the following components: a tachogenerator, a DC motor, an encoder and an inertia load. This modular experimental setup supports real-time design and implementation of advanced control algorithms, and is interfaced with the MATLAB/Simulink using specific RT-DAC4/USB board for transferring the measured signals from the tachogenerator and encoder, and the control signals to the power interface unit. The angular position θ of the DC motor shaft is measured by the incremental encoder, and the angular velocity ω is where |U (t)| ≤ 1 and V max = 12 [V ].
In order to identify the model of the system, the identification tool within Modular Servo Toolbox, which operates directly in the MATLAB/Simulink environment, is used. The identification procedure is also given in (Inteco, 2011). The following transfer function is obtained where K s = 184.73 and T s = 1.3s.
By denoting x 1 = θ and x 2 = ω, the state space model of servo system iṡ x 2 = ax 2 + bu + w where a = −1/T s and b = K s /T s , and w represents the Coulomb friction defined by treated as the unmodeled disturbance. The sampling period is set to T = 0.01s, and the discrete-time state space model is given by In all three sets, the nominal MPC, v, is calculated by the nominal model only, and SMC, ν, is used as the auxiliary controller to eliminate the disturbance.

A. Nominal MPC
In order to show the system response, when only the nominal MPC is applied, the first set of the digital simulation and real-time experiment is conducted. The following control It is shown that both nominal and real states respect the constraints defined by eq. (74), but there is discrepancy between the responses of the nominal model and real plant. This demonstrates the lack of robustness of the nominal MPC when it is applied to the real-time DC servo system in the presence of disturbance.

B. Tube MPC with traditional SMC
The traditional SMC defined by eq. (10), as an auxiliary controller of Tube MPC, is applied to cope with the disturbance. The two control components are now constrained separately.
The constraints for the nominal MPC and SMC are defined by which satisfy eq. (73), i.e. −1 ≤ v + ν ≤ 1. The new state constraints are calculated by using the tightening procedure described in Section 4. The tightened state constraints used for the nominal system are now defined by

C. Tube MPC with chattering free SMC
The same control and state constraints, defined by eqs. (75), (76) and (77), respectively, are used herein. The nominal MPC is applied to the nominal model first. The obtained nominal and real-time system responses are illustrated in Figures 10 and 11. After that, the Tube based MPC with the chattering free SMC is applied to the real DC servo system. The nate from noise existing in angular velocity signal taken from tachogenerator ( Figure 11 ). Therefore, they are not caused by chattering phenomenon. As in the previous experiments, all states and control signals respect the defined constraints.

Conclusion
In this paper, the Tube MPC with a SMC as an auxiliary controller is studied in order to improve the robustness of the overall system. Due to the presence of the SMC component, it is necessary to tighten the constraints of the nominal MPC part. The online calculations for SMC are not time consuming and this is also true for the nominal MPC, handling only the nominal model, which results in lower online calculation requirements. The traditional and chattering-free SMC algorithms are introduced in Tube MPC in order to reject disturbances and to achieve better perfor-mances of the real system. Therefore, the procedures for calculating the required constraints tightening are derived for the both cases. The good characteristics of the proposed control algorithms are demonstrated by conducting several digital simulations and real-time experiments on the DC servo system. Proof of Theorem 3.1 The vector sequence (g k , g k+1 , . . . , g k+i , . . .) , denoted by (g k+i ), converges point-wise to the limit g ∈ R nu if each element of g k+i converges to the corresponding element in g. In other words, (g k+i ) is convergent if lim g k+i = g, i.e. if for every real vector > 0 there exists the natural number N u ( ), such that is the positive (negative) vector sequence if g k+i ≥ 0 (g k+i ≤ 0) for i = 0, 1, 2, . . .. For multipleinput systems, it is probable that the elements of vector g k have different signs, as they represent the switching functions of SMC inputs. After splitting the vector g k onto two sub vectors g + k and g − k with separated positive and negative elements of g k , respectively, eq. (17) can be rewritten as where ∆ + u , ∆ − u ,(KE) + ,(KE) − , w + k+i and w − k+i are diagonal matrices and sub vectors obtained from ∆ u , KE and w k+i by extraction. Similarly, the theorem's condition given in eq. (18) can be expressed by .5) and the domain G defined by eq. (19) as It is obvious that (g + k+i ) and (g − k+i ), defined by eq. (A.2), are positive and negative sequences, respectively.
Let us prove now that (g + k+i ) enters the domain G + in finite time for k 0 ≤ i ≤ N and remains in that area. The proof is similar in the case of (g − k+i ) with respect to G − . If eq. (A.4) is true, then and g + k+i+1 < g + k+i so there exists a positive diagonal matrix Q k+i = diag{q 1 k+i q 2 k+i . . . q n + u k+i }, (0 < q j k+i < 1, j = 1, 2, . . . , n + u , n + u + n − u = n u ) such that where g + k+i and g − k+i+p (p ∈ N ) can be written as giving the following inequality ( > 0) According to Cauchy's theorem, the convergence of vector sequence (g k+i ), satisfying eq. (A.13), is proved. Its convergence domain is directly satisfying eq. (A.10).
Let us now show that system trajectory enters the domain G + in finite time. The sequence (g + k+i ) converges inside domain G + , so it is limited and lim i→∞ g + k+i = g + ∞ . Assume that g + k > ∆ + u 1 + Ω + . According to eq. (A.2) Suppose that g + k+i never enters the domain G + . For i → ∞ , we obtain

Equation (A.16) implies that the vector series
is convergent, and its general element ∆ + u 1 − (KE) + w + k+j converges to zero as j → ∞ , i.e.
that contradicts eq. (A.4), and the initial assumption that g + k+i never enters the domain G + is false. Moreover, g + k+i enters the domain G + at time instant k 0 which is bounded by the maximal element of vector It is obvious that the length of the prediction horizon N should be greater than k 0 and selected in accordance with eq. (A.18).
Therefore, we have proven that g + k+k0+1 ∈ G + and, by induction, the latter can be generalized to for every m > 0. The sign of g k may change at each time step, causing the chattering in that way, but g k will stay in G + . Having demonstrated that eq. (A.23) is satisfied if eq. (A.4) is valid, the proof ends.

Appendix B Proof of Theorem 3.2
Assume that g k / ∈ G. Then, eq. (22) becomes eq. (17) and the proof is similar to the one discussed in Appendix A. This means that g k+k0 ∈ G where k 0 is determined by eq. (A.18). Let g j k be the j th element of g k and assume that corresponding element (KEw k+k0 ) j < 0 . Then g j k+k0+1 = g j k+k0 − δ j u − |(KEw k+k0 ) j | (B.1) < Ω j − |(KEw k+k0 ) j | < δ j u − |(KEw k+k0 ) j | < δ j u where δ j u and Ω j are the j th elements in the diagonal of ∆ u and in vector Ω, respectively . Then, from eqs. (B.1) and (22) we have If (KEw k+k0 ) j > 0, g j k+i will continue to decrease and, after k 1 time instants g j k+k0+k1 ∈ {g j k+i : −δ j u − Ω j < g j k+i < 0} and 2) and (B.5) we have that once g k enters G, it will stay in it, i.e. g k+i+1 = KEw k+i ∈ G (B.6) and, therefore, there is no chattering in sliding mode.