Shape-Independent Model Predictive Control for Takagi-Sugeno Fuzzy Systems

Predictive control of TS fuzzy systems has been addressed in prior literature with some simplifying assumptions or heuristic approaches. This paper presents a rigorous formulation of the model predictive control of TS systems, so that results are valid for any membership value (shape-independent) with a suitable account of causality (control can depend on current and past memberships and state). As in most fuzzy control results, a family of progressively better controllers can be obtained by increasing Polya-related complexity parameters, generalising over prior proposals.


Introduction
Takagi-Sugeno (TS) systems are a widely-used tool to exactly model nonlinear systems via the so-called sector-nonlinearity approach [1], on a compact modelling region Ω.Subsequently, stability analysis and control design tasks can be carried out on the TS models.Control techniques for TS fuzzy systems based on LMIs [1] have been deeply developed in recent years, see the review paper [2] and references therein.Of course, the TS+LMI approach is conservative with respect to an "ideal" nonlinear control approach [3], but it allows solving the problems via convex optimisation tools, derived from linear systems and related to the linear parameter varying (LPV) approach [4].As the model is only valid on a compact set, the results are also limited to this modelling region, usually the largest Lyapunov level set in Ω (slightly larger sets can be actually proven, see [5]) or an inescapable set in disturbance-rejection problems [6,7].A recent alternative to the LMI approach in TS systems is the asymptotically exact geometric polytope manipulation approach in [8,9], extending well-known results in the robust-polytopic control field [10].
In realistic applications, there is always control saturation; most LMI conditions conservatively require the control action to avoid saturation in the outermost Lyapunov level set [1], others allow saturation [11] but they cannot prove improvement with respect to non-saturating laws or require iteration/Bilinear Matrix Inequalities, see [12].Note also that if the operation point is not at the center of actuator range, constraints are non-symmetric; however, in an LMI setup constraints must be forced to be symmetric in the vast majority of cases.The above-referred geometric approach can naturally handle such nonsymmetric constraints.
Apart from sheer stability, optimality of a quadratic index for TS systems is a problem that has been also addressed in many works in literature in the socalled guaranteed cost setup [1,12,13,14,15]; the "guaranteed cost" terminology stems from the fact that the controllers in such references only prove an "upper bound" to the actual cost.
The above-referred works use a 1-step version of the Belmann equation, i.e., V (x k ) ≥ L(x k , u k ) + V (x k+1 ), L being some quadratic "step cost", to prove that V is a cost bound via some LMIs.Obviously, the natural multi-step extension to the Bellmann equation is V (x k ) ≥ N −1 j=1 L(x k , u k ) + V (x k+N ); such optimisations are the key point behind model predictive control (MPC) approaches, very popular in applications nowadays [16].Linear MPC under quadratic stage and terminal cost are routinely solved via very efficient quadratic programming (QP, convex) optimisation and, also, there are well-known stability and optimality guarantees [17] emanating from classical LQR theory.In the uncertain case, the so-called minimax predictive control is addressed in the well-known references [18,19,10]; state an input constraints are naturally handled in the resulting constrained optimisation setups (QP, LMI).These "uncertain" MPC approaches can, of course, apply to TS systems but, they are conservative in the same way as "robust linear" controllers are, as they do not exploit the fact that membership functions are measurable in on-line operation so control can depend on them.On the other hand, the "ideally" least-conservative framework would be directly applying nonlinear MPC [20,21] to the original nonlinear system from which the TS model was originated; the main drawback is the fact that computational cost of nonlinear MPC is high (convexity cannot be guaranteed) and, in many cases, convergence of the optimisation iterations is, if they actually converge, to a possibly local minimum.
So, fuzzy MPC stands in the "middle ground" between the well-developed linear/robust cases and the elusive nonlinear one.However, quite a few of the fuzzy MPC earlier literature works have significant shortcomings.For instance, the widely-cited work [22] resorts to clustering but, actually, it carries out nonlinear MPC on the resulting identified fuzzy model, via nonconvex branch-and-bound optimisation.The work [23] computes a linear MPC by "freezing" the memberships at a particular instant and assuming they will be constant in the future; this might work in practice, but it lacks theoretical justification in fast transients.The work [24] presents an interesting approach in which a sequence of quadratic cost bounds and state-feedback gains solves (suboptimally) the MPC problem.The great advantage is its computational tractability; however, it is well-known that even for the linear case, under constraints, the optimal value function is not quadratic in the state, so the approach is conservative.Recent works such as [25] discusses networked interval type-2 systems, but their results are still based on the 1-step equation discussed above so they are not solving a "proper" multi-step problem as most MPC literature understands.There are other related works on LMI-based suboptimal MPC for Wiener and Hammerstein models [26] or for input-output LPV systems (uncertain impulse-response coefficients) [27] anyway, these non-TS representations are intentionally out of the scope of the present manuscript, as well as other multi-agent/cooperative/stochastic versions of MPC [28,29].
Actually, the papers whose approach is more similar to the one pursued here are, for instance, [30,31,32,33].Nevertheless, for brevity in this introduction, further discussion and comparative analysis with these related results are presented in later sections of this manuscript, once notation and problem statement have been established.
The objective of this paper is adapting the linear MPC for TS fuzzy models in a rigorous way, suitably posing the problem and proposing linear matrix inequality conditions for the resulting optimisation.The key idea is introducing the shape-independent concept to the model predictive control problem formulation on TS fuzzy systems.This concept implies that the MPC is valid not only for the actual non-linear system originating a TS model, but also for any other possible realization of the membership functions.This is a source of conservatism with respect ot pure nonlinear-MPC but, in exchange, systematic LMI conditions can be posed.Importantly, the proposed controllers make use of the fact that future memberships will be measurable at the moment of computing the future control action, absent in other literature proposals.
The structure of the paper is as follows: Section 2 presents preliminary results, notation and problem statement; Section 3 presents the main contribution on shape-independent fuzzy-MPC, prediction model and constraints; this section is followed by a stability analysis in Section 5, where terminal cost and controller are discussed in depth.Two brief sections present a shape-dependent MPC variation and a comparative analysis and discussion.Section 7 develops an example, and a conclusion section closes the paper.

Preliminaries
Considerer a nonlinear discrete-time system, given by: where x k ∈ R n represents the state vector, u k ∈ R m the control actions and f x (0) = 0.The above system can be expressed in compact regions X, U con-taining the origin as a Takagi-Sugeno fuzzy model with r rules in the form In this representation µ i (x) are the membership functions which, for later manipulations, they will be grouped in the vector of memberships which belongs to the standard simplex: The regions where system (1) is modelled will be assumed to be polyhedra where inequalities in vectors are understood to hold element-wise.
As discussed in the introduction, well-known LMIs have been developed to synthesise state-feedback controllers for systems in the TS form (2). In particular, this paper will root on the so-called guaranteed-cost results [1,12] and will propose predictive controllers which improve on the obtained cost bounds.

Fuzzy predictive control: problem statement
Generic Predictive control problems [16,17] are based on solving finitehorizon constrained optimisation problems.Let us discuss how the predictive control problem can be formally stated in a TS framework.
First, in a TS fuzzy system, let us define a membership-dependent cost index: In this way, some non-quadratic costs can be embedded into the framework of this paper.If a controller u k = g(µ(x k ), x k ) is under operation, the infinite cost is just a function of the initial state, to be denoted as J ∞ (x 0 ).Abusing the notation, if arguments are clear the above index will be shorthanded to J ∞ .Then, the fuzzy version of the finite-time problems in predictive control requires defining a N -step cost as: where, for convenience, memberships are arranged as a matrix: and µ ∈ ∆ N +1 will indicate that each column belongs to ∆.The term V (µ(x N ), x N ) will be denoted as terminal cost and the term L(µ(x k ), x k , u k ) as stage cost; u denotes the set of control actions {u 0 , . . ., u N −1 }.Under mild assumptions, if there exists a controller such that V (µ(x), x) ≥ J ∞ (x), then stability can be proved [17] for initial states such that J N (µ, u, x 0 ) is finite.Generic nonlinear predictive control (NMPC, [16]) would, subsequently, try to obtain a sequence of future inputs u * j , j = 0, . . ., N − 1 so that J N (µ, u, x 0 ) is minimised (understanding µ to be obtained from the simulated states under u * from an initial condition x 0 ).The basic problem of the above generic NMPC approach is the fact that the optimisation problem is, in general, non-convex and without guaranteed convergence or execution time bound in real-time implementation, and also stability guarantees are missing.Disregarding µ, system (2) can be considered as an uncertain polytopic one, so such early linear minimax results cited in the introduction can be applied to them to design robust controllers, but in such a case the knowledge of the membership functions in on-line operation is not exploited.
The objective of this paper is providing a Polya-based asymptotically exact solution to a shape-independent version of the above predictive control problem: under suitable constraints, where u := g(µ, U ) will denote a causal controller parametrisation, detailed in Section 3. Causality refers to the fact that the actual shape of the membership functions is unknown at design time but it will be known when implemented, and future memberships will, too, be eventually known.The proposed solution will be in LMI form.

Homogeneous polynomial notation
In order to solve fuzzy control problems, homogeneous polynomials in memberships are widely used (for instance, in the asymptotically exact solutions in [34,35,9], etc.) Considering a degree d ∈ N, d ≥ 1, and state x, all monomials of degree d in the memberships can be expressed as where α ∈ N r is a vector of natural numbers (understood including zero) such that |α| := r i=1 α i = d.In the sequel, the following notation will be used to denote a degree-d homogeneous polynomial in memberships, say Ξ d µ , which will be expressed as: where n α Ξ α are the polynomial coefficients, factorised as Ξ α and a combinatorial number n α , defined to be the multiset permutation of α in the set |α| = d, i.e.: It can be proved that: In Appendix A, some well-knwon results on products of homogeneous polynomials and Polya-theorem, later used, are recalled.State dependence µ α (x) will be shorthanded to µ α if clear from the context.

Remark 1
In order to use this notation in TS models, we will understand where ξ ∈ N r is, forcedly, a vector with a single element equal to 1, being the rest zero, and we understand A ξ to be equal to A i , with i being the unique index such that ξ i = 1.As n ξ = 1, we omit writing it at the left-hand side of (14).
Homogeneous polynomials in delayed instants.The appearence of memberships of states from x 0 to x N in the problem statement motivates extending the fuzzy summation notation to encompass different instants.This paper will discuss homogeneous polynomials in membership functions, evaluated at several instants of time, arranged as (8), so µ ik := µ i (x k ).The degree of the homogeneous polynomial at different time instants may differ: the homogeneous polynomial in memberships at time k will, by assumption, have degree d k .In order to compactly handle such situation, we will introduce a degree vector d ∈ N 1×(N +1) , conformed with said elements d k .
Considering now (11), as a different vector α will be needed for different instants, the definition of α will be generalised to being a matrix, using its k-th column to index a particular monomial at instant k: considering a matrix of natural numbers α ∈ N r×(N +1) , notation α k for k = 0, . . ., N , will denote the k-th column of a matrix α, considering the first one to be indexed by zero, in order to be consistent with the idea that the degrees correspond to µ(x 0 ).Also, in this matrix case, |α| will denote the vector of dimension 1 × (N + 1) formed by the column-wise sums, i.e., whose element at position j is |α k | = r i=1 α ik .Now, notation µ α will represent a monomial in the membership functions (at different instants of time), given by: and n α is It can be proved in a straightforward way that n α µ α = 1.The generalisation of (11) to the multiple-instant setting will be denoted by: where d will now be a degree vector, d ∈ N N +1 indicating the degree at each instant of time of the overall fuzzy summation.
Example 1 Consider a prediction horizon N = 3 and two rules, r = 2. Example membership monomials for some degree matrices α are, say:

Shape-independent fuzzy predictive control
Let us now completely state the shape-independent problem (9) rooting from ( 7) and its associated constraints.
At the moment of computing u k in ( 7), memberships µ(x 0 ), . . .µ(x k ) will be known; hence, we will search over controller depending on these memberships, i.e., u k := g k (µ(x 0 ), . . ., µ(x k ), U k ), where U k is a vector of numeric parameters for g k ().Juxtaposing all controls, defining the set of control laws u := {u 0 , u 1 , . . ., u N −1 }, and the set of numeric decision variables U := {U 0 , . . ., U N −1 }, where U k are those associated to control u k at instant k, we will express them with the shorthand u := g(µ, U ).A homogeneous polynomial structure for g and U will be detailed shortly below.
Under the above causality constraints, then, the shape-independent fuzzy MPC problem would amount to computing: subject to, for k = 0, . . .N − 1, x k ∈ X (21) where T is a so-called terminal set (see Section 5 for details on how to obtain it to guarantee stability).Note that, if memberships are considered as "uncertainty" and the controller cannot depend on them, the above problem statement would be that of the so-called minimax predictive control [18,19].Intermediate cases with partiallyknown membership components [36, Eq. ( 27)], may also be conceived, but they will be left out of the scope of the discussion for brevity.

Relevant prior results in literature
As discussed in the introduction, some prior results on the topic have already been dealt with in literature.The approaches may be classified as follows: Infinite-horizon (guaranteed cost) fuzzy control.Sets N = 0, and conditions such that V (µ(x), x) ≥ J ∞ (x 0 ) are set up (usually in LMI form) to find both V and an associated controller g(µ, x); if feasible, finiteness of the cost index is guaranteed (and stability under mild well-known conditions).Conditions are valid for any x 0 (usually, in a certain ellipsolidal region), unknown at design time.The reader is referred to, for instance, [1,12,13,14].Also, an outputfeedback setup can be found in [15].
On-line infinite-horizon fuzzy control.Adapting the above to a particular measured x 0 , the resulting LMIs yield Lyapunov functions and controller gains depending on the measured state.The reader is referred to, for instance, [32,37,38,33].These results are considered by their authors, too, to be model predictive control; we are not proposing any notation/terminology debate here, as it is just a matter of taste.The important fact is that, as the cost bound is decreasing in time, it can be used to prove feasibility in future time and stability, see the cited references for details.Note that dependence on the first membership µ(x 0 ) can be removed, as it can be directly measured.This implies a minor modification to the problem statement; for clarity in exposition, it will be deferred, and discussed later on in Section 4.4.
One-step predictive control.We will classify with this label the literature works where, given a measured x 0 , the control action u 0 is itself a decision variable, instead of resulting from a feedback gain as in the above-cited setups.Later on, a terminal state-feedback setup u k = g(µ(x k ), x k ) is assumed from k = 1 onwards.See, for instance [30] for details on the approach.In our setup, this proposal will amount to setting N = 1 (later on discussed in more detail).
Multi-step predictive control.This approach is the closest to the original MPC idea: control variables u 0 , . . ., u N −1 are computed once x(0) is measured, assuming that a terminal state-feedback controller is applied from instant N onwards [16,17].In [31], a membership-dependent future control action is considered.
Next section will present the formulation of a more general proposal of the above multi-step predictive control.In Section 6, a discussion and comparison of our proposal with some of the above-referred approaches will be presented.

Main Result
Causal fuzzy controller parametrisation.As future values of membership functions are unknown, the control action cannot depend on them, as discussed in earlier sections.So, u 0 will be chosen to be an homogeneous polynomial in µ(x 0 ) with, say, degree c 00 , u 1 will be a polynomial in µ(x 0 ) and µ(x 1 ), with degrees c 10 and c 11 in µ(x 0 ) and µ(x 1 ), respectively, and so on with u 2 , . . ., u N −1 .
We will introduce notacion c [0] = (c 00 , 0, 0, 0), c [1] = (c 10 , c 11 , 0, 0), . . .so c kl indicates the degree of u k in µ(x l ), for l ≤ k, 0 ≤ k ≤ N .Based on the above discussion, we will consider a control action to be applied to the TS system (2) to be given by: where u α,k are control decision variables conforming U k , and c [k] ∈ N N is a user-defined degree vector, with the above-discussed causality constraints (i.e., elements k + 1 to N equal to zero).This parametrisation generalises that in [31], being the one in that work recovered by setting c [k] = (1, 1, . . ., 1, 0, 0, . . . ) with a total of k − 1 ones, see Section 6 for further comparison.
4.1.Prediction model 1-step closed-loop model.With the above control law (23), at time k, the closedloop successor state would be: We will vectorise all decision variables in u α,k , i.e., stacking them in a column vector U k defined as: using any arbitrary enumeration ordering for α, for instance lexicographic.Conversely, we will invert the vectorisation using notation u α,k = E α,k U k where E α,k is a matrix that selects the suitable vector elements, according to the chosen order in the vectorisation operation.With this notation, the input at instant k can be written as: Hence, With suitable manipulations, using Corollary 2 and, if so wished Polya expansion (Corollary 1), we can express the closed-loop model (details omitted for brevity) as: where, with S c [k] e k γ defined according to (A.9) as and e k denotes the vector whose elements are all zero except the k-th one, equal to 1. Abusing the notation introduced in Remark 1, A ξ , when indexed by a matrix ξ such that |ξ| = e k , should be understood as the consequent matrix A i , where i is the row number at k-th column of the single element of ξ equal to 1.Note also that, if we choose q [k] = c [k] +e k then no Polya expansion is carried out; for any other larger choice of elements of q [k] the expression of the Polya expansion in Corollary 1 is implicitly considered in (28).
Prediction along the full horizon.In MPC, from an initial state x 0 , predictions of x 1 , x 2 , . . .x N must be made.So, the above 1-step ahead prediction must be nested, and expression (28) generalised.For instance, given x 0 , we could predict x 2 with an homogeneous controller of degree vector c [0] := (c 00 0 . . .0) at instant 0, and a controller which, at instant 1, can depend on memberships at instants 0 and 1, with degrees c [1] := (c 10 c 11 0 . . .0), the prediction of x 2 , without any Polya expansion, would be: Following similar steps, predictions of x 3 , . . .x N can be crafted as follows: In expression (31) and in those of larger horizon, products of model matrices at different instants appear in convolution-like formulae.In order to set a compact notation for such products, let us consider the product of matrices from instant j to k.The auxiliary degree vector 1 jk , used in (34) and ( 35) below, will be defined as the vector whose components j to k are equal to 1, being the rest equal to zero.Thus, matrices ξ such that |ξ| = 1 jk are those in which a single element is equal to one in each of the columns from j to k.Now, for such ξ, we will define: Example 3 Considering an index matrix |ξ| = 1 23 , in a 2-rule model, we would have, say: Note that in the model predictions, the above matrix would be multiplied by With the above notations, the classical convolution formula in linear MPC Bu j is extended to the TS case as follows: Theorem 1 For any degree vector q ∈ N N +1 , with q j ≥ max i (c jj + 1, c ij ), j = 0, . . ., N − 1, the prediction of x k+1 as a function of x 0 and future controls in the form (26) can be expressed as: being where S 1 0k γ is built as defined in (A.4), and: Proof: Carried out by exhaustively repeating the analogous operations to (31) for x 3 , etc., omitted for brevity because of the tedious nature of the developments.
Finally, the k-step prediction, for degree vector q, can be written as been U = (U T 0 , . . ., U T N −1 ) T and G γ,k is formulated by extracting out x 0 and U from Ξ γ in (34), in order to express (33) in matrix form (36).
Hence, juxtaposing G γ,k in column form, resulting in a matrix to be denoted as G γ , the full prediction from k = 1 up to k = N can be expressed as follows:

Constraints on decision variables
This section discusses how to enforce the constraints on states and inputs.Indeed, from prediction model (37), future states depend on x 0 and U ; thus, constraints ( 20)-( 22) must be formulated as constraints on decision variables U .
Input constraints.Carrying out a Polya expansion, if so wished of (26), the input at instant k can expressed, taking any degree vector h ∈ R N , so that h j ≥ max i (c jj , c ij ), j = 0, . . ., N − 1, as: where W γ,k , from (A.6), is defined as: At this point, as (38) is an homogeneous polynomial, in order to enforce u k ∈ U k = 0 . . .N , as required by constraint (20), we formulate the constraint on the polynomial coefficients, in terms of the decision variables U , using matrix S and vector s in (5) as: indeed, if the above holds, as n α µ α = 1, we have Finally, juxtaposing W γ,k in column form, as a matrix to be denoted as W γ , the full prediction inputs from k = 1 to k = N − 1 can be expressed as follows: Future state constraints.In a similar way (details omitted for brevity), the condition x k ∈ X, for k = 1 . . .N − 1, required in (21), is expressed in term of the decision variables U and the current state values x 0 , for any Polya expansion of degree h such that h j ≥ max i (c jj + 1, c ij ), j = 0, . . ., N − 1, as the sufficient condition Finally, for stability reasons to be later discussed, the terminal set must be defined as: where P c µ is a homogeneous polynomial of degree c, following notation (17).Taking the prediction model ( 28) for x N : and applying well-kwown Schur-complement manipulations, we have that x N ∈ T, required in (22), is equivalent to the matrix inequality: If a Polya expansion of the polynomials above is done up to degree vector h ≥ q, h ≥ l, the following (asymptotically exact) sufficient conditions for x N ∈ T in LMI form are obtained, requiring the coefficients of the referred expansion to be positive semidefinite:

LMI formulation of the model predictive control problem
Considering the cost index in (7), in order to cast the MPC problem as LMI, the k-step cost L(µ(x k ), x k , u k ) is defined as the following expression, quadratic in the state and inputs, but involving homogeneous polynomials of degree d in the involved matrices: and the terminal cost V (µ(x N ), x N ) is defined as: We are now in conditions of presenting the main result of this paper: Theorem 2 Given x 0 and the degree-q prediction model (37), the cost index J * in (9) fulfills J * (x 0 ) ≤ δ * if, for a degree vector h ∈ N N +1 , h ≥ q, h ≥ f , h ≥ l (see (54) below for definitions of f and l), the following LMI optimisation is feasible: subject to: Proof: Considering the index J N in ( 7) and ( 9), let us prove that it the LMI conditions in the theorem hold, once u is parametrised as (26).Indeed, replacing in J N the prediction model (linear in x 0 and U ), a quadratic expression in control decision variables U arises.In order to remove the quadratic dependence on the decision variables U in (55), the Schur complement is applied, resulting in the condition: According to Definition 1 and Proposition 2, the degree of the polynomial (56) can be extended to any complexity parameter h, being h ≥ f and h ≥ q.Thus, the Polya expansion of (56) yields: As all µ k ar positive, the inequality (57) will hold if the inequalities (50) hold; henceforth, δ will be an upper bound of the cost index J N .Note that it is also needed that future states belong to X, future control actions must lie in U and the state at instant N must be driven inside the terminal set.In order to ensure that, the constraints (51)-( 53) have been added to the problem as sufficient conditions to constraints (20)-( 22) in Section 4.2.Note that (51) and (53) are elementwise (scalar) constraints instead of full matrix inequalities.
Feasible set of initial conditions.For fixed horizon, terminal cost and terminal sets, the set of x 0 yielding feasible LMIs (50)-(54) in Theorem 2 (if x 0 were now considered to be a decision variable) is a convex LMI set.Such a set will be denoted as shape-independent feasible set.

Shape-dependent solution (known µ(x 0 )).
Note that the MPC problem stated in Section 3 is fully shape-independent.Thus, any x 0 in the shape-independent feasible set discussed on page 15 would be guaranteed feasible for any TS system, whatever the value of µ(x 0 ) happened to be.However, in actual implementation this would be suboptimal, given the fact that µ(x 0 ) would be known so the solution needs not to be valid for "all" possible µ(x 0 ) but only for the currently measured value.So, if matrix µ contained only the unknown memberships from instants 1 to N , the shape-dependent MPC problem would be stated as: subject to the same constraints as in the original problem.Note that the first control action does not need to depend on the membership function µ(x 0 ) so the set of decision variables U can be considered to be an arbitrary scalar u 0 for the first step (instead of using U 0 defined in ( 25)), and the unaltered vectors This entails some minor modifications to the previously-presented setting in Theorem 2: the part of the prediction model G α computing x 1 from x 0 is a single linear model x 1 = A(µ(x 0 ))x 0 + B(µ(x 0 ))u 0 , thus, the degree of all polynomials involving µ(x 0 ) can be set to zero, because all "vertices" are the same.For brevity, details on these modifications are left to the reader.Note that the feasible set of x 0 in problem (58) (and the associated constraints) would not be a convex LMI set as µ(x 0 ) may be any arbitrarily complex nonlinearity, see the example on Section 7.This set will be denoted as shape-dependent feasible set.

Stability
As in classical predictive control, in order to prove stability, the terminal set and terminal costs must be computed assuming there is a so-called terminal controller which guarantees that the infinite cost ( 6) is bounded by the terminal cost.Let us discuss the details.

Terminal Controller Theorem 3
The system (2) is stabilizable and the cost index J ∞ in (6), with step cost (47), is bounded by x 0 , with c ∈ N being a degree parameter, if there exist matrices P α , and K α , for all |α| = c, such that the following LMIs are feasible: for all γ ∈ N r such that |γ| = q, with any arbitrarily chosen q ≥ max(c+1, d), and all θ ∈ N r such that |θ| = c, with S c γ , S d γ and S c1 γ defined in (A.4) and (A.9), and the control action is: Proof: Omitted for brevity, as it closely follows well-known guaranteedcost results, see [1,12,35], plus some standard Polya-expansion argumentations incorporated here.
Optimal controller.It is straightforward to check that, for instance, adding the condition and minimising δ would guarantee that J ∞ ≤ δ * x T 0 x 0 , δ * being the optimal solution.

Computation of terminal sets.
Consider the terminal set defined in (43).The developments below will discuss how to obtain a bound on λ such that T ⊂ X and the control action (60) is admissible, i.e., u ∈ U, for all x ∈ T. Note that, as T is a subset of every Lyapunov level set1 , once T is entered, the state trajectory will never leave T in the future under the terminal control law computed with LMIs in Theorem 3.
Theorem 4 For a TS fuzzy system (2), with input (60), the terminal set T verifies T ⊂ X and u(x) ∈ U for all x ∈ T, if λ is the minimum positive scalar such that, for all |γ| = q, being q ≥ c, where R i denotes the i-th row in matrix R and S j the j-th row of S; and i ranges from one up to the number of rows of R, j from one to the number of rows of S.
Proof: As the shape-independent levels sets of the Lyapunov function from Theorem 3 are symmetric, because V (µ, x) = V (µ, −x), in order for such a level set to lie inside the (possibly non-symmetric X), the referred level set must belong to the symmetric set According to [39], it is well known that the set Substituting P c µ = |γ|=c µ γ P γ by an arbitrary Polya expansion So, from Theorem 8, a sufficient condition for the above is that each of the homogeneous polynomial coefficients fulfills the inequality, which is what (62) states.
On the other hand, the control action must belong to set U for all states inside the sought level set.Again, as done with the state constraints, symmetry of the level set and linearity in the state of the terminal controller needs to enforce u ∈ U sym , where , and the condition below holds, then the control action will be admissible for all x in T. Concentrating on the inequality over the braces, if congruence with matrix P c µ is applied and, later on, a Schur complement, we get the equivalent condition: Finally, an arbitrary Polya expansion can be carried out, yielding This inequality will hold if all the coefficients of the polynomial are positive, which amounts to (63).

Stability: main result
The first result to be stated considers the standard open-loop application of the computed control actions until the terminal set is reached: µ ≥ 0 and H d µ > 0 for any µ ∈ ∆, and control actions (26) are applied once the solution to the optimisation in Theorem 2 is obtained, from instant k = 0 to k = N − 1, and the terminal controller is applied from instants N onwards, then, the system will reach the origin asymptotically and the cost (6) will be bounded by δ * .
Proof: Note that, the theorem 2 includes the condition (52), which forces that control actions (26) steer the system to terminal set (43) within N steps.If once the system is inside the terminal set the terminal controller is applied, the system will reach the origin asymptotically with it.With regard to a bound on cost (6), the terminal cost bounds the infinite cost (Theorem 3) once x N ∈ T and the total cost (including the transient until T is reached) is bounded by δ * discussed in the statement and proof of Theorem 2.
In predictive control setups, receding horizon implementation is a common approach [17].In our context, what we will understand as receding horizon implementation is solving problem (58) at each step, after current state is measured.Under such setting, we can state the following results.
First, inspired in [30], we will consider the particular case N = 1 which renders the setup in Algorithm 2 of the cited work (except for minor changes in the terminal LMIs and the use of a simpler non-fuzzy cost index).
Theorem 6 In the case N = 1, receding horizon implementation of the solution to (58) is recursively feasible and stable.
Proof: With N = 1, application of the first control action u 0 will drive the state x 1 to the terminal set.Following standard MPC argumentations analogous to [30,Theorem 1] or [17], the terminal control action will be feasible at next step, hence the solution to (58) will produce an equal or better cost bound, which can be proven to be a Lyapunov function.For brevity, details are left to the reader.
For receding-horizon implementations with longer horizons, the following theorem can be stated: Theorem 7 Receding horizon implementation of the solution to (58) is recursively feasible and stable, too, in the case the terminal controller (60) is polynomial in the membership functions (quadratic Lyapunov function P α ≡ P ), and controller degree vectors c [k] in (23) fulfill: and c (N −1)(N −1) ≥ c, with c being the complexity of the terminal controller gain K c µ(x) in (60), and the Polya expansion degree vector h in Theorem 2 verifies h i−1 ≥ h i , i = 2, . . ., N .
Proof: Under the conditions in the theorem statement, the terminal controller is polynomial in the memberships.With such controller, a parametrisation of high enough degree on the control actions u k in (25), k = 1, . . ., N − 1 can allow an argumentation similar to the previous theorem, as follows.
Let us consider, for illustration a case with {0, 0, 0, 0, c} Terminal control The symbol indicates that, as the membership µ(x 0 ) is known when solving (58), fuzzy dependence on such variable is not necessary in the control parametrisation.
One step later, the open-loop application of the solution from the previous step would have complexity: {c 21 , c 22 , 0, 0, 0} u 3 : {c 31 , c 32 , c 33 , 0, 0} u 4 : {0, 0, 0, c, 0} Terminal control u 5 : {0, 0, 0, 0, c} Terminal control Hence, if we apply the optimisation problem arising from the first table in receding horizon, in order to ensure feasibility and non-increasing cost bound, computations with the complexities on the first table must be able to replicate the second one (i.e., the second table must produce a feasible sub-optimal result on the problem associated to the first one).To achieve the latter condition, we need to take two facts into account: 1. measuring µ(x 1 ) at next instant will ensure that fuzziness on it is not needed (so can indeed replicate the results with complexities h 1 , c 11 , c 21 and c 31 for any value of these degrees).2. Except the ones, discussed above, each element of the first table must be larger or equal than the corresponding element on the second one, i.e., c 11 ≥ c 22 , c 21 ≥ c 32 , c 22 ≥ c 33 , and for the last action to be more powerful than the terminal one, we need c 33 ≥ c.Additionally, we need the Polya expansion degrees to fulfill The above discussion justifies the inequalities in the complexity parameters in the theorem statement.
Once the complexity issues are cleared out so recursive feasibility and a monotonically decreasing cost bound are proven, the standard MPC stability argumentations can be applied, in the same way as commented in the proof of Theorem 6, omitted, again, for brevity.
In the above theorem, the motivation of considering a quadratic Lyapunov function is to avoid the terminal controller being rational in the memberships.Indeed, as the control actions from k = 1 to N − 1 are polynomial in the memberships we cannot ensure that they can replicate the rational terminal control law (this is why only N = 1 allowed to prove Theorem 6).In this way, the stability results of, for instance [31], can be considered a particular case of Theorem 7.

Discussion and comparative analysis
As discussed in the introduction, other references have dealt with predictive control in a fuzzy context.Some comparative discussion with a few references appears in the introduction and in Section 3.1, focusing the similarities and differences with other guaranteed-cost and infinite-horizon proposals.We deferred to this section a comparison with other more recent LMI-based results on similar problems.
Regarding comparison with [40], their approach pursues similar goals to ours.However, their determination of the terminal set and controller is conservative, in the sense that, first, their proposed controller is a PDC one (see their equation 4) and, second, they pose conditions for all i, j, l in their Theorems 1 and 2, where i affects the current process and Lyapunov function vertex model, j affects the controller vertex and l is the next-instant Lyapunov function vertex.In that way, say, controller vertex 1 proves stability and cost bounds for every i and k so their results cannot be better than those from a robust linear controller.Our terminal controller even for the same level of Polya complexity parameter will, hence, prove a better cost bound.Note, however, that they consider uncertainty in the TS vertex models as well as persistent disturbances so such conservative steps are clearly justified and needed to synthesise the terminal controller; on the other hand, our approach purposely does not consider uncertain models, although incorporating polytopic uncertainty such as that in [36] is indeed possible, by omitting dependence of the controller parametrisation (23) on unknown components; if not component were known, our control actions would be the plain scalars in minimax predictive control [19].
In [32], a state-dependent guaranteed-cost control is proposed so that, once x k is measured, a cost bound γ can be guaranteed solving some LMIs ensuring that control action does not saturate.As states get closer to the origin, better cost bounds can be found by increasing the controller gain.In fact, when our terminal set is eventually reached, the same bound would be obtained as LMIs are equivalent.Our improvement lies in the fact that we allow full saturation until the terminal set is reached.
One of the most relevant recent results related to our proposal is [31].Indeed, the problem statement is similar.The controller parametrisation in the cited work, [31, eq. ( 10)], proposes a control action with multilinear dependence on past memberships (u k depends on memberships up to k − 1).However, it does not incorporate the membership at instant k, as we do.In that way, not being "fuzzy" in the last computation has the advantage of avoiding multiple-sum formulations so an elegant single-sum prediction model ensues, see [31, eq. (11)], but loses expressive power.Actually, as already commented, using our notation, the case in the above work is a particular case of our proposal setting c [k] in (23) to be (1, 1, 1, 1, . . ., 1, 0, 0, . . . ) containing k − 1 ones, the rest zeros.Indeed, the resulting prediction model with q [k] in (28) to be (1,1,1,1,...,1,0,0,...) containing "k" ones would yield [31, eq. ( 11)] as a particular case.The LMI problem [31, eq. ( 27)] is a particular case of our approach (Theorem 2).In summary, in our proposal, we generalise the control action parametrisation from multilinear to polynomials in past and current memberships (and Polya expansions), with arbitrary degree, as well as membership-dependent cost-index weightings.last, in the discussion of Theorem 7, we also argued that the stability results in Ding's work can be considered a particular case of ours.
As a last comparison, the work [30] is also relevant.Indeed, the proposed terminal set in [30, eq. (35)] can be proved to be the same as (43).

Examples
Consider a 2-rule TS system (2) with the local models and membership functions defined as: Input and state constraints are defined as: For simplicity, non-fuzzy weighing matrices H and F will be employed, with d = 0: Terminal set and controller.The terminal controller has been computed with degree c = 1, with LMIs (59), ( 62) and (63) expanded to degree q = 20.The terminal set T (shape-independent) is drawn with black line2 on Figure 1.
Simulation.For illustration, trajectories (using the less-conservative shape-dependent solution) from x 0 = (0, 4) (solid blue) and x 0 = (0.8, −4.9) (dashed blue) are shown in the phase plane in Figure 1, and in the time-domain, jointly with the control action, in Figure 2. A third simulation with x 0 = (4.8,−4.8) from the shape-dependent feasible set also appears in Figure 1.
Comparative analysis.Let us compare our feasible set with that obtained in [40, Theorem 3] with the same horizon N = 4, depicted in orange, in the above figure.Importantly, the cited work tries to maximise the volume of an ellipsoidal feasible set via logdet convex optimisation subject to LMI constraints.Conservatism of the proposed conditions there makes our approach to yield significantly larger feasible sets both in the shape-independent case (purple) and in the shape-dependent one (green).Considering x 0 = (0, 4), the comparison of the guaranteed cost bounds with different approaches is summarised on Table 1 below: Note that the main goal of the paper [40] is not to minimize the cost, but to maximise the feasible area; hence, a large penalty with respect to [32] is incurred.Nevertheless, our approach beats both [32] in cost and [40] in size of the feasible set 4 .
Application example: fuzzy predictive control of a CSTR.Consider the continuous stirred tank reactor (CSTR) in [15, Eq. ( 52)].From a certain operating point, state x 1 is the reactive concentration increment, and state x 2 is its temperature increment.Control action u is a coolant temperature increment.The state constraints (modelling region), as well as the input ones are set as: The reader is referred to [15] for full modelling details.The shape-dependent fuzzy predictive control problem (58) is proposed with N = 3, c [0] = ( , 0, 0, 0), c [1] = ( , 1, 0, 0), c [2] = ( , 1, 1, 0), h = ( , 2, 2, 1), the cost function is quadratic with values Q = 1 and R −1 = 0.25.The trajectories of output and control action are the blue lines displayed in Figure 3.The proposed fuzzy MPC problem had a computational time cost, on an Intel i7-4790 @3.6GHz, 8GB RAM, of 0.56 seconds.
For comparison, we also programmed the controller proposed in [33] whose output and control action are represented by the red dashed lines in Figure 3. Setting an initial state x = (0.35, 5.6) T , the achieved cost value J REAL = 35 0 (y 2 k + 0.25u 2 k ) of our proposal is 345 whereas that of [33] is 379.The outputfeedback version in [15] was also tested, achieving a larger cost of 439, as intuitively expected (simulation omitted as results are not comparable, out of the scope of this work).

Conclusions
This paper has presented a generalisation of the predictive control approach to Takagi-Sugeno systems.The approach obtains better results than prior literature both in achieved cost and in the feasible region, due to a suitable control action parametrisation as a function of future memberships and the use of Polya relaxations.Of course, many guaranteed-cost results in prior literature are particular cases considering an horizon of a single sample; in fact, such results are used to build the terminal controller and terminal cost needed in the developments.
Theorem 8 (Polya theorem [34,41]) The summation M d µ is positive for all µ ∈ ∆ if there exists q such that all coefficients of the expanded degree-q polynomial E(q, M d µ ) are positive.Furthermore, if M d µ > > 0 for all µ ∈ ∆ there exists a finite q such that E(q, M d µ ) has all its coefficients positive.where the set S dl γ is defined as: Obviously, if any of d, c, l were equal to 1, the corresponding combinatorial number n α , n β , n ξ would be omitted, as in Remark 1.

Proposition 2
The summation over a single instant k, H q µ(x k ) in (11), can be represented as a multi-instant sum (i.e., with α being a matrix): where degree vector d is such that d k = q and the other values of the vector d are arbitrary non-negative integers.
Proof: From the fact that |α|=di n α µ(x i ) α = 1 for all i and non-negative natural d i , we have

Figure 1 :
Figure1: Terminal and feasible sets, plus three simulated trajectories of the MPC controller (in blue).Shaded green regions are unfeasible in both shape-dependent and shape-independent solutions.

Figure 2 :
Figure 2: Time response (state and control) of the trajectories on figure 1 for two initial conditions.