Generalized Stabilizing Conditions for Model Predictive Control

This note addresses the tracking problem for model predictive control. It presents simple procedures for both linear and nonlinear constrained model predictive control when the desired equilibrium state is any point in a specified set. The resultant region of attraction is the union of the regions of attraction for each equilibrium state in the specified set and is therefore larger than that obtained when conventional model predictive control is employed.


Introduction
Model predictive control of a dynamic system determines the current control by solving a finite horizon optimal control problem; the optimization yields the optimal control sequence that satisfies the state and control constraints as well as a terminal state constraint. The control applied to the system is the first element of the optimal control sequence. Because the optimal control problem has a finite horizon, the implicit model predictive control law is not necessarily stabilizing. To overcome this obstacle, one of two approaches is usually adopted; either the optimal control problem is modified by the addition of a terminal cost and constraint [1] or a sufficiently large horizon is chosen [2]. With the first approach, for each horizon, the feasible set for the optimal Communicated by Lars Grüne. B David Mayne d.mayne@imperial.ac.uk 1 Imperial College London, London, UK control problem (the set of states for which there exists a control sequence satisfying all constraints) is control invariant (if the terminal constraint set is control invariant) and is a subset of the feasible set for a larger horizon, a feature that ensures recursive feasibility. The size of the feasible set, which is also the region of attraction, depends on the terminal constraint set and can be smaller than desired if the terminal constraint set is small.
In the literature, a variety of approaches to the problem of enlarging the region of attraction, including increasing the horizon length [3], has therefore been proposed. In [4], a version of model predictive control that can cope with a family of set-points is developed; pseudo-linearization is employed to yield a transformed system that has the same linearization at each set-point in the family, permitting use of the same terminal set and terminal penalty for each set-point. In [5], the single terminal constraint set is replaced by a nested set of three terminal sets in which the first set is the terminal constraint set and each set is a subset set of its successor; each point in any set in the sequence can be steered into the preceding set by an admissible control. In the optimal control problem solved online, the last terminal set is employed at the initial time, the last but one set at the next time, followed by the second last and, thereafter, the conventional terminal set. The implicit model predictive control law is therefore time-varying in the first few iterations during which the optimal cost can increase. The domain of attraction is now the set of states that can be steered by an admissible control into the third set in the sequence; this is clearly larger than the region of attraction for conventional model predictive control. In [6], an enlargement of the region of attraction is obtained by constraining the unstable states to a set of states that can be controlled to the terminal set in a fixed number of steps. In [7,8], the terminal set and, hence, the region of attraction is enlarged by using a local saturating control law in place of the conventional linear control law in the terminal constraint set. A larger terminal region can also be obtained by employing a detuned local control law in the terminal set at the price of losing local optimality as shown in [9]. An alternative approach for linear systems is interpolation-based MPC, discussed in [9], in which a set of local stabilizing controllers is employed together with their associated invariant sets in which all constraints are satisfied. Each of these sets could serve as a terminal constraint set. A larger terminal set that is the convex hull of these sets is obtained by expressing the current state as a convex combination of states in these sets and employing the local control law that is the same convex combination of the associated controls. The interpolation procedure is extended in [10] in which a set of state trajectories, each terminating in a standard terminal constraint set, is computed. For each trajectory, the terminal constraint set is translated and scaled; the scaling is chosen to ensure satisfaction of the state and control constraints in each set. The enlarged terminal set is the convex hull of these translated and scaled sets. A similar but simplified approach is adopted in [11] where a single trajectory is employed. For each state in this single trajectory, a set is chosen to ensure satisfaction of the state and control constraints. The terminal set is then defined by the convex hull of these sets. Attention has also been given [12][13][14], to the problem of increasing the region of attraction in the context of the tracking problem. Very recently, a method for increasing the region of attraction for conventional and economic MPC [15] has been proposed and analysed.
The main purpose of this paper is to present simple methods for increasing the region of attraction to any desired equilibrium state in a given set since simplicity is an important requirement in many process applications. The method bears some resemblance to phase 1-phase 2 algorithms [16, §2.6] for constrained optimization in which, in the first phase, the algorithm satisfies the constraints and, in the second phase, minimizes the cost function while maintaining feasibility. It therefore also has some resemblance to dual-mode model predictive control as proposed in [17] for tracking in that, in the initial stages of the algorithm when the terminal constraint is not satisfied, the controller acts to satisfy this constraint even if, as a consequence, the optimal cost increases. Although developed independently, the proposed methods for linear systems are similar in some respects to the procedure described in [18,Chapter 5] and [19] that addresses a different problem, distributed economic model predictive control. Whereas our algorithm generates a sequence of set-points converging to a neighbourhood of a known set-point, the algorithm in [19] does not assume the economic set-point is known; instead, a sequence of set-points converging to the economic set-point is computed by a distributed optimization algorithm. Both methods therefore require generation of a sequence of terminal constraint sets; the centre of each terminal constraint set is the 'virtual' set-point that converges to, or approaches, the desired setpoint ( [19] does not employ the term 'virtual set-point'). In [19], to ensure recursive feasibility, the size of the terminal constraint set is reduced at each time. In this paper, both the centre (the virtual set-point) and size of the terminal constraint set is altered (kept constant in Linear Control Algorithm 1 and increased or decreased in Linear Control Algorithm 2) at each iteration to ensure recursive feasibility; the virtual setpoint reaches the desired set-point in finite time. A similar procedure is also employed in [14]. A more detailed comparison of the two methods and their similarities and differences is given in Sect. 2.1 after relevant notation has been defined.

Increasing the Region of Attraction for Constrained Systems
In the sequel, we use f (·) to denote a function irrespective of the number of its arguments. B(x, c) denotes {z : |z − x| ≤ c}. For any two integers i, j, I i: j denotes {i, i + 1, i + 2, . . . , j − 1, j}.
The system to be controlled is described by Given a reference r * ∈ R, R a compact subset of R p , the control objective is to stabilize the equilibrium state x r * . For each r ∈ R, the corresponding equilibrium pair (x r , u r ) is defined as the solution, if it exists, of For simplicity, we assume that y and u have the same dimension and that Eq. (2) has a unique solution u r = φ u (r ), x r = φ x (r ) for all r ∈ R; in the linear case (i.e. x + = Ax + Bu, y = C x), φ u (·) and φ x (·) are linear (φ u (r ) = M u r, φ x (r ) = M x r ) and our assumption is equivalent to the assumption that the zero frequency gain G(0) = C(I − A) −1 B is non-singular. If a solution to (2) does not exist, there does not exist a constant u such that y(t) = r for all t. The analysis can be extended to the case when there are more outputs than inputs [20, §1.5.1]. The system is subject to state and control constraints x ∈ X and u ∈ U. It is assumed that R satisfies (x r , u r ) ∈ X × U for all r ∈ R and that X is polyhedral and that U is a polytope. For each (x, r ), the optimal control problem P N (x, r ) solved on line is in which is the set of control sequences u that, given (x, r ), satisfy the state and control constraints and the terminal constraint The control applied to the system is u 0 (0; x, r ), the first element of u 0 (x, r ); the implicit model predictive control law is, therefore, κ N (·) defined by κ N (x, r ) := u 0 (0; x, r ). Let x κ N (i; x, r ) denote the solution at time i of the closed loop system ) if the initial state at time 0 is x. The equilibrium state is asymptotically stable with a region of attraction X N (r ) := {x : U N (x, r ) = ∅} for the closed loop system if the terminal 'ingredients' X f and V f (·) satisfy the usual stability condition: The stability conditions are satisfied if V f (·) is the value function for the unconstrained infinite horizon problem or if V f (·) is the infinite horizon cost for the unconstrained system with any locally stabilizing controller We choose R and c > 0 such that the terminal constraint set X f (r ) satisfies ). If f (·) is nonlinear, it is possible, under mild conditions, to determine a terminal cost function V f (·) and associated constraint set X f (r ) satisfying the stability conditions.
We show below that, for any r * ∈ R, it is possible to design a model predictive controller such that x r * is asymptotically stable with a region of attraction X N (R) := ∪ r ∈R X N (r ) [rather than X N (r * )].

Model Predictive Controller with Extended Region of Attraction: Linear Systems
Some simplifications can be made if the system being controlled is linear ( f (x, u) := Ax + Bu, y = C x). We consider first the simple case when the terminal constraint sets are ellipsoidal and differ only by translation since, for each r, ; P and c do not vary with r (P r = P and c r = c for all r ). A consequence is that, for all r ∈ R, the local stabilizing control law is affine ( In conventional model predictive control, the value function of the online optimal control problem decreases monotonically with time. In the modified version described below, the value function, as in phase 1-phase 2 optimization, may increase until feasibility of the terminal inequality constraint with r = r * is achieved; thereafter, the value function decreases monotonically. The control algorithm can be simply stated: Linear Control Algorithm 1: Step 0: (Initialization:) Determine, at the initial state x, a r ∈ R such that x ∈ X N (r ) (i.e. such that P N (x, r ) is feasible).
To establish properties of the algorithm, two preliminary results are required.
Proof Because the state and control constraints for ; it then follows from the inequality above that x + N ∈ X f (r ) ∩ X f (s) so that the initial state Our main result follows.  , (A, B) is stabilizable, and the stability condition is satisfied with X f (r ) ⊂ X and κ f (X f (r ), r ) ⊂ U for all r ∈ R. Then, for all r * ∈ R, x r * is exponentially stable with a region of attraction X N (R).
Proof Let (x i , r i , λ i ) denote the value of (x, r, λ) at the ith iteration of the algorithm. The sequence {x r i : Step 2 of the algorithm and Proposition 2.2 that P N (x i+1 , r i+1 ) is feasible and |x r * − x r i+1 | P ≤ |x r * − x r i | P − δ. Hence, there exists a finite time j such that x r j ∈ X f (r * ); using the definition of λ 0 , x r j+1 = x r * . For all i ≥ j + 1, the algorithm is conventional MPC so that, from known results, x r * is exponentially stable with a region of attraction X N (R).
To compare the results above with those of [19], we note that [19] defines (using our notation) the terminal cost function V f (·) as we do but defines the terminal constraint set by X f (r ) := {x : E(x − x r ) ≤ α r } (since r in our analysis replaces time t in [19]), whereas our terminal constraint set is defined in Sect. 2.1 by X f (r ) := {x : V f (x, r ) ≤ c} (c independent of r ). Lemma 1 in [19] then shows how to choose the virtual reference x r + and α r + so that, for all x ∈ X f (r ), x + := f (x, κ f (x, r )) satisfies x + ∈ X f (r + ); crucially (Assumption 2(ii) in [19]) X f (r ) has the contraction property Using this, [19,Theorem 1] shows x r i converges to x r * as time i → ∞. In Linear Control Algorithm 1, X f (r ) has a different property; Proposition 2.1 shows that if This enables us to show that x r i enters X f (x r * ) in finite time. In [19], the terminal constraint set X f (r i ) decreases in size (α r i decreases) as i increases, whereas in Linear Control Algorithm 1, the terminal constraint set has constant size.
The fact that X f (r ) in Linear Control Algorithm 1 has constant size (since c is constant) constrains each X f (r ) to be small if, for some r ∈ R, x r is close to the boundary of X and/or κ f (x, r ) is close to the boundary of U. To overcome this disadvantage, we replace Step 2 in Control Algorithm 1 by an optimization algorithm that provides a terminal set X f (r ) that increases in size (c r larger) where it is possible to do so without transgressing the state or control constraints. The terminal set is now defined by X f (r ) := {x : V f (x, r ) ≤ c r } in which c r varies with r . The optimization algorithm,P N (x, r ), x ∈ X f (r ), employed in Step 2, is defined for all x by min λ,c in which X := {x : x ∈ X, κ f (x, s) ∈ U} and s := r + λ(r * − r ). Let (λ 0 (x, r ), c 0 (x, r )) denote the solution ofP N (x, r ). It is shown in the 'Appendix' thatP N (x, r ) may be reformulated as an LMI problem for which standard algorithms are available [21]. The modified control algorithm is: Linear Control Algorithm 2: Data: r * and c r * .
Step 0: (Initialization:) Determine, at the initial state x, a r ∈ R and a c = c r > 0 such that x ∈ X N (r ) (i.e such that P N (x, r ) is feasible).
Because c r varies with r , we have to modify the analysis given for the Linear Control Algorithm 1. We assume that R is chosen so that there exists a c > 0 such that c r ≥ c for all r ∈ R. Next, let e be computed as in the proof of Proposition 2.1, and, for all r ∈ R, let d r := c r − e. It can be shown that Proposition 2.1 still holds for all r ∈ R with c replaced by c r and d replaced by d r , i.e. e is independent of r . However, because the size of X f (r ) now varies with r , modification of Proposition 2.2 is required; we need to know how c s varies with s = s(λ) := r + λ(r * − r ) ∈ R.
Consider the function λ → φ(λ) that maps λ to the minimum of |x s(λ) − z| P with respect to z lying in a hyperplane H; this minimum is the distance, in norm |·| P , of x s(λ) from the hyperplane H. Since x s = M x s is linear, the function φ(·) is affine. Hence, the function that maps λ to the distance (in norm | · | P ) of x s(λ) from the boundary of the polyhedron X , which has a finite number of faces, is continuous and piecewise affine with a Lipschitz constant L < ∞. For each x ∈ R n , let c(λ) denote the largest value of c > 0 such that the set {x : V f (x, s(λ)) ≤ c} = {x : |x − x s(λ) | P ≤ √ 2c} is a subset of X . From the discussion above, the function λ → √ 2c(λ) is identical to the function λ → φ(λ) and is, therefore, piecewise affine with Lipschitz constant L x 1 Thus, there exists δ ∈]0, δ] such that, if |x r * − x r | P > δ , |x r + − x r | P ≥ δ . Theorem 2.2 then holds for the modified algorithm.
The fact that the terminal constraint set X f (r ) can increase in size (see Fig. 1) when this is possible means that the region of attraction X N (R) is larger than that for Linear Control Algorithm 1 (size fixed) and for [19,Algorithm 1] in which the size of X f (r ) decreases (in order to attain a more complex objective) as a consequence of the different contraction property for X f (r ).

Model Predictive Controller with Extended Region of Attraction: Nonlinear Systems
With some further assumptions, the above procedure may be modified to be applicable to nonlinear systems. The control problem is described, as above, by Eqs. (1)- (5). Because of the nonlinearity, the control algorithm has to be significantly altered. To see this, consider the obvious extension of the above procedure. For each r ∈ R, the terminal constraint set X f (r ) is defined by X f (r ) The function V f ( · , r ) and the set X f (r ) can be determined, for each r ∈ R, as in [22] using linearization of f (·) at each (x r , u r ) provided that the linearization of f (·) is stabilizable at each x r . Unlike the linear case, P r and c r now vary with x r making simple extensions of Propositions 2.1 and 2.2 impossible. Because of the difficulty in determining X f (r ) for each r ∈ R and devising a suitable algorithm for adjustment of the virtual reference, we do not pursue this approach.
The impasse can be resolved by using a simpler stabilizing condition, namely a terminal equality constraint that is satisfied for each r ∈ R admittedly at some cost since a terminal equality constraint has a smaller region of attraction than a terminal constraint set; X N (r ), the region of attraction, is now the set of states that can be steered to x r in time N along a trajectory that satisfies the state and control constraints. It is desirable that X N (r ) is a large subset of R n ; this is a strong assumption.
The optimal control problem P 1 N (x, r ) solved at each state (x, r ) now has the form: where V N (·) is defined by As before, the set of feasible points for P 1 , r ), and the model predictive control that is applied to the system is κ N (x, r ) := u 0 (0; x, r ). The implicit control law is, therefore, κ N ( · , r ), and the successor state is x + = f (x, κ N (x, r )).
The control algorithm requires, in addition, the solution of another optimization problem P 2 N (x, r ) that yieldsŝ =ŝ(x, r ) ∈ R and, hence, (xŝ, uŝ), a potentially improved equilibrium pair; x is not altered. The optimal control problem P 2 N (x, r ) is defined by is an optimality function [16, §1.1] in that θ(x, r ) < 0 implies |x r * − x r | can be reduced and θ(x, r ) = 0 implies |x r * − x r | cannot be reduced. We assume: Assumption 1 There exist an ε † > 0 and a c 1 ∈]0, 1/2[ such that (i): for all r ∈ R, the ball B(x r , ε † ) ⊂ X N (r ) and (ii): Assumption 1(i) excludes the degenerate case in which the feasible set X N (r ) tends to a set with an empty interior as r tends to some r ∈ R.
Assumption 1(ii) ensures |x r * − x r | can be reduced for all r ∈ R, r = r * , and by a finite amount for all r such that |x r * − x r | ≥ ε † /2 provided |x − x r | ≤ ε † /2. It is easily shown to be satisfied with c 1 = 1/2 if f (·) is linear in which case the set {x r : r ∈ R} is convex. Assumption 1(ii) ensures that Step 3 is executed only a finite number of times and permits the use of standard descent algorithms for solving P 2 N (x, r ); this is a strong assumption since it excludes the case when the problem P 2 N (x, r ) has local minima. Let B(z, ε) := {x : |x − z| ≤ ε}.
The control algorithm for the nonlinear case may now be defined: Nonlinear Control Algorithm: Data: x, r * , c 1 ∈]0, 1/2[.
Step 0: Determine, at the initial state x ∈ X N (R), a r ∈ R such that x ∈ X N (r ).
Step 2: r ) = (x + , r ) and go to Step 1. Else go to Step 3.
We show below that the algorithm repeatedly executes Steps 2 and 3, reducing |x r * − x r | each time until the reference r remains constant with value r * . Execution of Step 1 then causes the state x i , where x i is the value of x at the ith iteration of the algorithm, to converge to x r * .
We can now state our main result.

Theorem 2.2
Suppose that Assumption 1 is satisfied, f (·) is continuous, that the parameters Q and R of the stage cost (·) are positive definite, that (2) has a unique solution for each r ∈ R, that (x r , u r ) ∈ X × U for all r ∈ R, that there exists a K ∞ function α(·) such that, for all r ∈ R, all x ∈ X N (r ), V 0 N (x, r ) ≤ α(|x − x r |) and that the stability condition is satisfied. Then, for all r * ∈ R, x r * is asymptotically stable with a region of attraction X N (R) := ∪ r ∈R X N (r ).
Proof Let x i , r i and ε i denote the values of the state x, the virtual reference r and ε at the ith iteration of the algorithm so that their initial values are x 0 = x and r 0 = r . The reference r i remains constant while the condition in Step 2 is satisfied. During these iterations, the solution of P 1 N (x i , r i ), with r i constant, is standard model predictive control so that x i converges to x r i . After a finite number of iterations, the state x i+1 = x + i is sufficiently close to x r i (x + i ∈ B(x r i , ε † /2)) for Proposition 2.3 to hold so that θ(x + i , r i ) ≤ −c 1 ε † ≤ −c 1 ε i . The condition in Step 3 is therefore satisfied. It follows that r i is updated to r i+1 = r + i =ŝ(x + i , r i ) and |x r * −x r i+1 | ≤ |x r * −x r i |−c 1 ε † . Steps 2 and 3 are then entered iteratively reducing |x i − x r i |, |x r * − x r i |, and ε i until Since both x + i and x r * lie in X N (r * ), r + i =ŝ(x + i , r i ) = r * . Thereafter, the algorithm iteratively solves P 1 N (x i , r * ). This is a standard MPC algorithm with a constant reference r * . Asymptotic stability of x r * with a region of attraction X N (R) follows.
As formulated, problem P 2 N (x, r ) is complex largely because x s is a nonlinear function φ x (s) of the reference s ∈ R. The algorithm can be simplified by merely requiring an approximate solution of P 1 N and P 2 N ; techniques for doing this are given in [16]. A possible alternative is to restrict the choice of s by minimizing with respect to λ ∈ [0, 1] and setting s = s λ = r + λ(r * − r ).

Linear System
Our first example is the second-order linear system x + = Ax + bu, y = C x in which The actual reference is r * = 0. The stage cost is ( P is the solution of the algebraic Riccati equation. The initial state is x(0) = (0.077, 0.385). The initial reference is chosen to be that point in R closest to y(0). The control horizon is N = 7. Feasibility of standard MPC can be recovered using a horizon N = 11. Figure 2a shows the trajectory of the virtual reference as generated by Linear Control Algorithm 2. Figure 2b shows the corresponding state trajectory. Figure 1 shows the terminal sets generated by the subalgorithmP N (x, r ) for sampled values of x and r .

Nonlinear System
The nonlinear system is a continuous stirred tank reactor based on a model described by in which x 1 is the product concentration, in which x j denotes the jth coordinate of the vector x. The controller horizon is N = 10 corresponding to 10 seconds since the sampling period is 1 second. The stage cost is (x, r, u) = |x − x r | 2 Q + |u − u r | 2 with Q = diag{1, 100}. The algorithm parameter c 1 is 0.001. The control is required to steer the output y to the reference r * = 0.4. The initial value of a feasible virtual reference was obtained by solving an optimization problem P 2 N (x(0), r ). The optimization packages described in [23,24] were employed for all our simulations. Figure 3 shows how the plant output y = x 2 converges to r * and also shows the convergence of the virtual reference r to r * . Figure  4 shows the locus of r ∈ R, the convergence of x from x 0 to x r * for the Nonlinear Control Algorithm with N = 10 and the convergence of x to x r * for standard model predictive control with N = 110; the desired state x r * is infeasible for standard model predictive control if the horizon N = 10 is employed; standard MPC requires, in this example, a value of the horizon N that is larger than 100.

Conclusions
The main objective of this note is to present a simple procedure for increasing the region of attraction of an equilibrium state when model predictive control is employed. Simplicity is important in the process industries, the main application area for model predictive control, mainly to increase ease of commissioning, operation and maintenance of controlled processes. It is for this reason that commercial software for model predictive control usually eschews provision for handling hard state constraints. In many applications, the output of a controlled process is required to converge to any reference r lying in a given set R; for such applications, a reasonable objective is the determination of a suitably large region of attraction X N (R) such that every state in X N (R) can be steered to any equilibrium state x r * , r * ∈ R. This note addresses the problem of regulation to any state in {x r : r ∈ R} and, hence, differs from several contributions in the literature where the concern is to increase the region of attraction for one equilibrium state, usually the origin. Because the initial state x does not necessarily lie in X N (r * ), it is necessary for the control algorithm to include a mechanism for adjusting both the current state x and a virtual reference r ; the reference r is adjusted from an initial value that is feasible for the initial state until the desired reference r * is feasible for the current state. This adjustment can cause the cost V N (·) to increase, a feature that causes difficulty in single-mode algorithms. Hence, we employ, as in [17], a dual-mode algorithm in which the initial objective is adjustment of a virtual reference r while maintaining feasibility of the current state until r * is reached. However, this feature of the algorithm depends on the ability of the algorithm to decrease |x r * − x r | at each iteration and this in turn requires the absence of local minima in the function r → |x r * − x r |, necessitating the rather strong Assumption 1 that is difficult to verify.