Robust Model Predictive Control for Takagi-Sugeno Model with Bounded Disturbances – Pólya Approach

This paper proposes a general robust model predictive control (MPC) approach for the constrained Takagi-Sugeno (T-S) fuzzy model with additive bounded disturbances. We adopt the homogeneous polynomially parameter-dependent (HPP) Lyapunov matrix with the arbitrary complexity degree and the corresponding HPP control law for the controller design. By applying the Pólya’s theorem and the extended nonquadratic boundedness property, a systematic approach to construct a set of sufficient conditions for assessing robust stability described by parameter-dependent linear matrix inequalities (LMIs) is established. The proposed approach is an improvement over existing approaches in terms of control performance and stabilizable model range. Numerical examples are provided to show the effectiveness of the proposed robust MPC approach.


I. INTRODUCTION
Takagi-Sugeno (T-S) fuzzy model has been widely used to approximate or even exactly represent nonlinear systems, whose basic idea is transforming the original system into a family of linear submodels [1]- [3]. For stability analysis of T-S model, many efforts are made based on the Lyapunov function, such as the common quadratic Lyapunov function in [4], the parameter-dependent quadratic Lyapunov function in [5], the piecewise-quadratic Lyapunov function in [6], the nonquadratic Lyapunov function in [7], and barrier Lyapunov functions in [8], [9]. To further improve the performance and reduce the conservatism, a general nonquadratic stabilization conditions are presented by the multiple-parameterization approach in [10]. In [11], [12], other general forms of relaxed stabilization conditions are derived by means of affine parameter-dependent Lyapunov functions. More details are available in [2], [3].
Model predictive control (MPC), as a widespread control technique being implemented in a receding horizon fashion, has advantages in constraints handling for multivariable plants, such as distributed MPC [14], industrial hierarchical MPC [15], and stochastic MPC [16]. Usually, at each sampling interval, MPC solves an optimal control problem, with the performance index being associated with the system evolutions over a prediction horizon, subject to physical constraints, where a sequences of control moves are treated as decision variables, but only the first control move among this sequence is implemented on the actuator. These actions are repeated in a receding horizon fashion. Since the future predictions for input/state/output are needed, which are obtained based on the system dynamic model, the accuracy of this model is crucial for the future prediction, which as a result, can influence the control performance [13]. The stability analysis and control synthesis for T-S fuzzy model by MPC approaches have been studied with variety. In [17], an interval type-2 fuzzy MPC approach is proposed for nonlinear networked systems. The model and controller are not required to share the same lower and upper membership functions. In [18], to improve the performance, a local stability approach is applied and an estimation of the domain of attraction is provided. The work of [19] investigates the robust fuzzy MPC, which uses the nonlinear local models. More relaxed results are achieved, and on-line computational cost is significantly reduced. The authors in [20] proposes a fuzzy generalized predictive control for T-S systems based on Kernel Ridge Regression strategy which learns the T-S fuzzy parameters from the input and output data. [21] utilizes the zonotopic set and interval matrices to bound the membership function errors. The controller parameters are stored in an off-line table for searching and robust tubes can be time-varying. In [22], the nonlinear multivariable predictive control is proposed for vehicle systems. For maintaining the robustness and stability, the controller design is based on LMI convex optimization. The cooperative fuzzy MPC is represented in [23] where the overall nonlinear plant consists of a group of parallel input-coupled T-S fuzzy models. For this cooperation, convergence and stability are guaranteed.
In order to deal with unknown disturbances, a series of paradigms have been elaborated (see e.g., [26]- [31], [33]). In [26], the homothetic tube-based approach is proposed, which maintains the state predictions of the linear model in the presence of disturbance within an on-line scaling tube centered at the disturbance-free model trajectory. The work of [27] proposes a tube-based MPC for nonlinear continuous-time model, and the feedback control law is optimized off-line. The authors in [28] utilized an integral non-squared stage cost and a non-squared terminal cost, so that the robustness of the resultant MPC is ensured without additional stability constraints. In [29], the input-to-state stability property is utilized for the quasi-min-max MPC design, and the first control move from the control sequence can be optimized directly. In [30], the notion of quadratic boundedness (QB) is utilized, based on which, the system is guaranteed to be quadratically bounded in the presence of disturbance. [31] proposes the full and partial dynamic output feedback MPC applying full Lyapunov matrix. The elliptical estimation error set is refreshed on-line based on optimized information of the last sampling instant. [32] aims at the norm-bounded model parametric uncertainty. The estimated state feedback gain and state estimator matrix are optimized on-line while the state estimator gain is designed off-line. In [33], sufficient conditions for computing the positively invariant set for T-S fuzzy systems are derived, and the terminal constraint set for 0-step and N -step control strategies are obtained. This paper characterizes MPC synthesis, based on improving the Lyapunov function, for constrained T-S fuzzy model with the bounded disturbance. Some general results for the positiveness of polynomials with matrix-valued coefficients (based on Pólya's theorem) is given in [35], where some complete characterization of the solution of parameter-dependent LMIs, usually arising in the robust stability analysis, is proposed. However, when the model parameters are timevarying uncertain, the results in [35] cannot be directly invoked. We deal with this issue in this paper, and contributions are summarized as follows.
1) The potentiality of Pólya's theorem is exploited. The general homogeneous polynomially parameterdependent (HPP) Lyapunov matrix whose complexity degree is tunable, and corresponding HPP control law, are applied.
2) By a generalization of the methods based-on Pólya's theorem and parameter-dependent LMIs, a series of finite-dimensional LMI relaxations, as sufficient stability conditions, are developed to robustly stabilize the resultant closed-loop system.
In [34], a general robust MPC approach for linear parameter varying (LPV) systems in the absence of bounded disturbance has been proposed, which can include many existing approaches with common quadratic Lyapunov matrices and state feedback laws (e.g., [24], [25]) as special cases. As compared with [12], [34], this paper handles the unknown but bounded disturbance. Since the bounded disturbance is incorporated in the model, the characterization of the closedloop stability is different from that in [34]. There are two other differences as compared with [34].
• While [34] introduces a free control move (i.e., the control move is the immediate decision variable), this paper does not. • While [34] utilizes the dilution parameter G, this paper does not.
The rationale of approximation is the same as in [34]: the complexity degree of HPP solutions is tunable for the proposed approach and, when it increases, the conservatism of the results reduces; the HPP Lyapunov matrix and HPP feedback gain matrix can asymptotically approximate any Lyapunov and feedback gain matrices which are continuous on the combining coefficient functions. Notation: The symbol induces a symmetric structure in the matrix inequalities. A variable with superscript * means the optimal solution to the optimization problem. m×n denotes the m × n-dimensional real matrix set. N + is the set of nonnegative integers. For the vector x and positive-definite matrix P > 0, x 2 P = x T P x. M ! denotes factorial of M . x(i|k) is the value of x at the future interval k + i, predicted at interval k. I is the identity matrix with appropriate dimension. For the column vectors x and y, [x; y] = [x T , y T ] T . ε P = {ξ|ξ T P ξ ≤ 1} denotes the ellipsoid that is associated with the symmetric positive-definite matrix P . The timedependence of the MPC decision variables is often omitted for brevity.

II. PROBLEM STATEMENT
Consider a class of T-S fuzzy systems, with its jth rule represented by where j ∈ {1, . . . , r} with r rules. Let θ(k) = [θ 1 (k); θ 2 (k); · · · ; θ ϑ (k)] be the measurable premise variable. H 2 , · · · , H (j) ϑ are the fuzzy sets. x ∈ R n , u ∈ R m , and w(k) are measurable state, input, and bounded disturbance vectors, respectively. The disturbance is persis-tent and satisfies w(k) ∈ ε Pw . Then the acquisition of T-S model can be described as where H (j) τ (θ τ (k)) denotes the grade of membership of θ τ (k) in H (j) τ , and h j (θ(k)) is the normalized membership function (abbreviated as the membership function), with h j (θ(k)) ≥ 0, r j=1 h j (θ(k)) = 1. For the sake of simplicity, h j is short for h j (θ(k)) and h j+ is short for h j (θ(k + 1)) for the rest of the paper.
The physical constraints to be considered is as follows, Ψ ∈ r×n can be any pre-specified forms of constraints on state x. The corresponding disturbance-free (reference) model is expressed as wherex andũ are the disturbance-free state and input signals. The control objective is to design an MPC controller which steers x(k), with the increase of sampling interval k, to converge to a neighborhood of the origin, while satisfying the constraints (4)- (5). Accordingly, the cost function is given as where Q > 0 and R > 0 are the given symmetric weighting matrices.x(0|k) = x(k). Note that the penalized signals arex andũ from the disturbance-free model (6), any other signals linearly dependent onx andũ are also allowed.

III. MODEL PREDICTIVE CONTROL SYNTHESIS
In this section, we propose the main result. Firstly, the standard MPC synthesis approach with the nonquadratic Lyapunov function is proposed with guaranteed recursive feasibility and closed-loop stability. Then, we apply the Pólya's theorem to extend the result such that the LMI conditions in the MPC optimization problem is relaxed.

A. FORMULATION OF OPTIMIZATION PROBLEM
This paper utilizes a g-degree HPP control law in the form of g is the parameter-dependent feedback gain, which is define as follows.
Firstly, introduce some definitions in order to be consistent with [35]. Let the HPP matrix with tunable complexity degree g be where η ∈ Ω r , p i ∈ N + , i = 1, 2, . . . , r, each η p1 1 η p2 2 · · · η pr r is a monomial. For all p ∈ K(g), G p ∈ n×n are parameterbased matrices. K(g) is the family of r-tuples, which is comprised of all the terms For example, for the case with g = 2 and r = 2, where P 0 (·), P 1 (·), . . . , P r (·) are continuous functions with respect to parameters η = [η 1 , η 2 , . . . , η r ] T ∈ Ω r and unknown µ ∈ R M . If there exists µ(η) ensuring P (µ(η), η) > 0 for all η, then a homogenous polynomial functionμ(η) exists such that P (μ(η), η) > 0 holds. Proof. See [35]. Remark 1: Lemma 1 implies that a homogeneous polynomial form of solutions, whose parameters lie in the unit simplex, is a very general form, i.e., can be readily transformed into any other continuous solutions of parameter-dependent LMIs. By considering the famous Weierstrass approximation theorem and applying Lemma 1, it infers that P g and corresponding F g can represent any Lyapunov matrix and feedback gain matrix that are parameterized by h j , j ∈ {1, 2 . . . r} as g increases.
By applying (8) on (2), the closed-loop system is obtained as In order to guarantee the stability of closed-loop system (10), the technique of quadratic boundedness (QB) ( [37]), which is primarily utilized for the state estimation problem and can be particularly useful for handling the system with bounded disturbance, is utilized to ensure that the state will stay in a quadratically bounded set. The definition and theorem of QB according to [37] are reviewed as follows.
Definition 1: For all allowable w(k) ∈ ε Pw , k ≥ 0, the autonomous linear system The QB condition is applied to many MPC problems with common quadratic Lyapunov functions. However, the resulting control performance can be conservative due to the utilization of the form of Lyapunov function. Since it is known that the extended nonquadratic Lyapunov method outperforms the common quadratic one, the extended nonquadratic boundedness is applied to characterize the closedloop property in this paper.
Let P g be the nonquadratic Lyapunov matrix. At inter- (10) is strictly nonquadratically bounded with an extended nonquadratic Lyapunov matrix By inheriting the results in [37] and extending to the case with the extended nonquadratic boundedness , the following conclusion can be obtained Lemma 2: For all allowable w(k + i), i ≥ 0, the following statements are equivalent. a) System (2) is nonquadratically bounded with an extended nonquadratic Lyapunov matrix P g . b From Lemma 2, we can obtain that the system (2) is nonquadratically bounded if at each sampling interval k, the following condition are satisfied: Since Hence, the condition (11) is equivalent to According to (10), (12) can be expressed in quadratic form as where ∆ : g . By eliminating the variables [x(i|k) T w(k + i) T ] T and invoking the S-procedure, it is shown that (12) is satisfied if and only if there exists a scalar α > 0 such that (see [37]) By substitute P g = γS −1 g , P g+ = γS −1 g+ , pre-and postmultiplying (14) by diag{S g , I} (which leaves the inequality unaffected), and applying the Schur complement, it is shown that (14) is guaranteed by Proposition 1: Reference model (6) is quadratically stable with the Lyapunov matrix P g .
Proof. According to Lemma 2, if (2) is quadratically bounded with the Lyapunov matrix P g , then condition (15) is satisfied. By neglecting the disturbance, (15) is reduced to g . Thus, by defining the HPP quadratic function V (x(i|k)) = x(i|k) T P gx (i|k), and substituting V (x(i + 1|k)) =x(i + 1|k) T P g+x (i + 1|k), we have Thus, the conclusion holds. Remark 2: The HPP quadratic function V (·) is in a very general form, which implies that it covers many existing Lyapunov functions. For example, by taking g = 0, the Lyapunov function in [24] is recovered; by taking g = 1, [25] is recovered. Based on Proposition 1, the state and input prediction of reference model (6) will converge to the origin, i.e., lim i→∞x (i|k) = 0 and lim i→∞ũ (i|k) = 0. Hence, summing (16) Since P g = γS −1 g andx(0|k) = x(k), let γ be the upper bound of (17), i.e., V (x(0|k)) ≤ γ, then the following holds: Similar to the procedure in [24], the input and state constraints in (4) and (5) are guaranteed by where z s,inf = min{u s ,ū s }, ψ s,inf = min{ψ s ,ψ s } and Z ss (Γ ss ) is the sth diagonal element of Z (Γ).
As the usual practice in robust MPC, the optimization problem is formulated as the following min-max form: Noth that by specifying α in the interval (0, 1), problem (21) will become a convex optimization which can be solved efficiently by the interior point method. The MPC optimization (21) is simply an extension of the approach in [24] to the case with nonquadartic Lyapunov function. We will further extend this approach by the utilization of the Pólya's theorem in the next section. Theorem 1: For (2), at any sampling interval, once there exists a feasible solution to optimization problem (21) at any interval k, then it will be feasible at k+1, and x will converge to a neighborhood of the origin. Proof. The proof contains two steps. 1)recursive feasibility. Supposed there exists feasible solution to (21) at sampling interval k. At next interval k + 1, we need to check that whether the constraints in (21) can still be satisfied. It is noted the optimization problems at k and k + 1 are different only in constraint (18) which contains the state x(k). At interval k + 1, let feasible boundγ(k + 1) = γ * (k). Constraints (15) and (18) at k imply V (x(i + 1|k)) ≤ V (x(i|k)) ≤ γ * (k) =γ(k +1) which, by applying the Schur complement, is shown to be equivalent to Thus, the optimal solution at sampling interval k still satisfies constraint (18) at interval k + 1. Other constraints can be naturally satisfied at k+1 by substituting the optimal solution at k. Hence, the constraints of the optimization problem are satisfied at k + 1.
Remark 3: For simplifying the presentation, we only consider the case when switching horizon N = 0 which is consistent with the benchmark work [24]. when N > 0, free control moves are added before the control law (8), this can be achieved easily by generalization, which is omitted here for brevity.
In this paper, we consider the homogeneous polynomial matrix X with degree g × g in the following form where each X p,q is a parameter-based matrix. Also, we have l j=1 η j = 1, l j=1 η j+ = 1. p = p 1 p 2 . . . p r , p 1 + p 2 + · · · + p r = g, and q = q 1 q 2 . . . q r , q 1 + q 2 + · · · + q r = g. Proposition 2: For the condition when (23) is positive, it is guaranteed that a set of sufficiently large d and d + exist which ensures the positiveness of all the coefficients of (η 1 + · · · + η r ) d (η 1+ + · · · + η r+ ) d+ X.
According to Proposition 2, the proof of (25) contains two steps.
Thus, the proof is complete. In summary, optimization problem (21) Similar to (21), problem (30) can be solved via convex optimization if α is pre-specified. However, note that the computational burden of (30) can be much heavier but the result is non-conservative. Corollary 3.1: If optimization problem (30) has feasible solution for a particular set of {g 0 , d 0 , d 0+ }, then it also holds for g > g 0 , d > d 0 , d + > d 0+ .
Proof. see [35]. Remark 4: An alternative methodology to calculate the HPP Lyapunov solution to (30) is by only increasing g and choosing d = 0, d + = 0. However, more decision variables are emerged by increasing g while the increase in {d, d + } brings the larger number of LMIs. If the computational efficiency is a crucial factor, one can simply reduce g, d and d + to a satisfactory level.

Example 1.
Consider the following discrete-time nonlinear system: where x 1 (k) ∈ [−β, β] with β > 0 and disturbance |w(k)| ≤ 0.5. Let membership function h 1 = (β + x 1 (k))/(2β) and h 2 = (β − x 1 (k))/(2β) be the combination coefficients. Then, nonlinear system (31) can be represented by the T-S VOLUME 4, 2016 model in the form of (2) with The T-S model can be different by simply changing β. Larger β implies the model in a larger region, which is more difficult to control. Choose Q = I, R = I, P w = 2, α = 0.3, w(k) = 0.5sin(k). In order to show the effectiveness of the proposed approach, (30)   of T-S model is enlarged by applying the approach in this paper. Moreover, in order to get the exact same β 0 , one can increase g while maintaining d and d + to be zero, or increase the set {g, d, d + } as a whole. To illustrate the effectiveness of different g, three special cases with {d = 0, d + = 0} for the same T-S model are considered, i.e., the case g = 0 that is applied in benchmark work [24], the case g = 1 that is applied in another work [25], x(k 1 ) 2 Q + u(k 1 ) 2 R to assess the performance . Calculate J sum for g = 0, g = 1, g = 2, respectively. As can be seen in Figure 3, the value of J sum tends to be small as g increases. Hence, it is concluded that a larger value of g brings the performance improvement.
To validate the effectiveness of {d, d + }, we consider three different cases with fixed g = 1, i.e., {d = 0, d + = 0} in [25] Figure 6. It can be observed that, the control performance is improved with {d, d + } increasing. Example 2.
Consider another benchmark example, i.e., a continuous stirred tank reactor (CSTR) whose continuous dynamics iṡ where C A is the concentration of irreversible exothermic reaction A in the reactor, T is the measurable reactor temperature, T c is the temperature of the coolant stream. The disturbance |w| ≤ 1. The objective is to regulate C A and T by manipulating T c . Corresponding parameters used are summarized in Table 2.
). By discretizing the continuous system with sampling period ∆t s = 0.2minutes, the nonlinear system (32)    Then choose Q = I, R = I, w(k) = sin(k). We use another way to modify the model through multiplying submodels by σ simultaneously. The maximal values of σ for the stabilizable T-S models are listed in Table 3. From the Table 3, it is shown that with values of {g, d, d + } increasing, σ becomes larger. Thus, the stabilizable model range is enlarged.
To illustrate the effectiveness via different complexity degrees g, we consider three cases for the same T-S model, i.e., {g = 0, d = 0, d + = 0}, {g = 1, d = 0, d + = 0}, {g = 2, d = 0, d + = 0}. Choose σ = 1, x(0) = [0.2, 4] T . The trajectories of states and input are depicted in Figures  7-8. Figures 7-8 show that the state and input evolve to the neighborhood of the origin without the constraints violation over the whole simulation horizon, so the closedloop system is nonquadratically bounded. Choose J g = 30 k=0 x(k) 2 Q + u(k) 2 R as the performance criterion. Calculate J g via g = 0, g = 1, g = 2, respectively, and we obtain that J (g=0) = 154.3521, J (g=1) = 133.0989, J (g=2) = 110.1759 (see Figure 9). As can be shown J g reduces as the value of g increases. It implies that performance is improved with a larger value of g. In this paper, a less conservative MPC approach for T-S fuzzy model with bounded disturbance is proposed. A general form of HPP Lyapunov function and the corresponding HPP control law are adopted to extend the previous approaches which are taken as special cases in this paper. The complexity degree is allowed to be tuned in order to balance the control performance and the computational efficiency. The proposed technique brings less conservatism as well as enlarging the stabilizable model range. Controlled systems subject to measurement noises widely exist in practice, our future attention is therefore paid on extending the proposed method to other systems with this issue, such as Markov system [38], [39], linear parameter-varying system [40], and networked control systems [41].