Maximal terminal region approach for MPC using subsets sequence

Abstract

For the sake of enlarging the domain of attraction of model predictive control (MPC), a novel approach of gradually approximating the maximal terminal state region is proposed in this paper. Given the terminal cost, both surrounding set sequence and subsets sequence of the maximal terminal region were constructed by employing one-step set expansion iteratively. It was theoretically proved that when the iteration time goes to infinity, both the surrounding set sequence and the subsets sequence will converge to the maximal terminal region. All surrounding and subsets sets in these sequences are extracted from the state space by exploiting support vector machine (SVM) classifier. Finally, simulations are implemented to validate the effectiveness of this approach.

Keywords

model predictive control (MPC) / terminal region / domain of attraction / support vector machine (SVM)

Introduction

Model predictive control (MPC), which approximates the infinite time-domain optimal control by finite horizon optimization, is well known as a suboptimal control algorithm. Due to the ability to handle control with state constraints, MPC has become quite popular recently. During the approximating process in MPC, the stability problem arises. To guarantee the stability of MPC, a terminal constraint and a terminal cost are added to the online optimization problem such that the terminal region is a positively invariant set for the system and the terminal cost is an associated Lyapunov function [1,2]. It is well accepted that the attraction domain of MPC can be enlarged either by increasing the prediction horizon or by enlarging the terminal region. Whereas, the former will subsequently increase the computational complexity, see Ref. [3] for example. Thus, the latter has attracted much more attentions in previous literatures.

In Refs. [4] and [5], an ellipsoidal set and a polytopic set, which were included in the stabilizable region of using linear feedback controller, were designed as the terminal region, respectively. While in Ref. [6], a saturated local control law was utilized to enlarge the terminal region. In Ref. [7], support vector machine (SVM) was employed to estimate the stabilizable region of linear feedback controller, and the subsequently estimated stabilizable region was utilized as the terminal region; thus, enlarged the terminal region dramatically. Moreover, it was proved in Ref. [8] that without terminal constraints, the terminal region can be enlarged by weighting the terminal cost. In Ref. [9], the enlargement of the domain of attraction was obtained by employing a contractive terminal constraint. Most recently, the domain of attraction was enlarged by the inclusion of an appropriate set of slacked terminal constraints into the control problem, as described in Ref. [2]. These previous literatures are of the same framework that the control law is supposed as known, and the corresponding terminal state region is calculated. Thus, there is no guarantee that the consequent terminal region is the maximal one in different control laws.

According to the given terminal cost-function, this paper proposes a novel approach that gradually approximates the maximal terminal state region based on the stability constraints. The fundamental principle of this method is the iterative procedure: At first, the sufficient conditions to guarantee the stability of MPC are presented and the maximal terminal region satisfying these conditions is defined. Then, given the terminal cost and an initial surrounding set or subset of the maximal terminal region according to the Lyapunov condition, a subsets sequence is obtained by using one-step set expansion iteratively. Moreover, it is proved that when the iteration time goes to infinity, this subsets sequence will converge to the maximal terminal region. Finally, the surrounding sets or subsets sequence are sequentially separated from the state space by exploiting SVM classifier (see Refs. [10,11] for details of SVM), and the maximal terminal state region is approximated gradually.

Model predictive control

Consider the discrete-time system as follows:

(1)x_{k+\mathrm{1}}=f(x_k,u_k),

where x_k\in R^n, u_k\in R^m are the state and the input of the system at the sampling time k, respectively. x_{k+\mathrm{1}}\in R^n is the successor state and the mapping f:R^{n+m}\mapsto R^n satisfying f(\mathbf{0},\mathbf{0})=\mathbf{0} is known. The system is subject to constraints on both state and control action. They are given by x_k\in X, u_k\in U, where X is a closed and bounded set, U is a compact set. Both of them contain the origin.

The online optimization problem of MPC at the sample time k, denoted by P_N(x_k), is stated as

(2)\begin{array}{l}\underset{u(i,x_k)\in U}{\min }{J}_N(u,x_k)={\displaystyle \sum _{i=\mathrm{0}}^{N-\mathrm{1}}q(x(i,x_k),u(i,x_k))}+F(x(N,x_k)),\hfill \\ s.t. x(i+\mathrm{1},x_k)=f(x(i,x_k),u(i,x_k)),\hfill \\ x(i+\mathrm{1},x_k)\in X,u(i,x_k)\in U,x(N,x_k)\in X_f,\hfill \end{array}

where x(\mathrm{0},x_k)=x_k is the state at the sample time k, q(x,u) denotes the stage cost and it is positive definite, N is the prediction horizon, X_f denotes the terminal region and it is closed and satisfies \mathbf{0}\in X_f\subseteq X, F(\cdot ) satisfying F(\mathbf{0})=\mathrm{0} is the terminal cost and it is continuous and positive definite.

Consider the following assumption at first.

Assumption 1　For the terminal region and the terminal cost, the following two conditions are satisfied [1]:

(C1) For any x\in X_f, there exists a Lyapunov function F(\cdot ) that satisfying the inequality

F(x)\ge \underset{u\in U}{\min }\{q(x,u)+F(f(x,u))\}.

(C2) X_f is a positively invariant set. For any x\in X_f, by using the optimal control resulting from the minimization problem showed in condition (C1), denoted by u_{opt}, we have f(x,u_{opt})\in X_f.

Let J_N^{*}(x_k) be the minimum of P_N(x_k) and {\mathit{u}}_N^{*}(x_k)=\{u_N^{*}(\mathrm{0},x_k),u_N^{*}(\mathrm{1},x_k),\dots ,u_N^{*}(N-\mathrm{1},x_k)\} be the optimal control trajectory. The control strategy of MPC is that, at the sample time k, u_N^{*}(\mathrm{0},x_k) is inputted into the real system, and at the sample time k+\mathrm{1}, the control inputted into the system is not u_N^{*}(\mathrm{1},x_k) but the first element of the optimal control trajectory resulting from the similar online optimization problem. At the sample time k+\mathrm{1}, the state is x_{k+\mathrm{1}}=f(x_k,u_N^{*}(\mathrm{0},x_k)) and the online optimization problem, denoted by P_N(x_{k+\mathrm{1}}), is the same as Eq. (2) except that x_k is replaced by x_{k+\mathrm{1}}. Similarly, let J_N^{*}(x_{k+\mathrm{1}}) be the minimum of P_N(x_{k+\mathrm{1}}) and u_N^{*}(x_{k+\mathrm{1}})=\{u_N^{*}(\mathrm{0},x_{k+\mathrm{1}}),u_N^{*}(\mathrm{1},x_{k+\mathrm{1}}),\dots ,u_N^{*}(N-\mathrm{1},x_{k+\mathrm{1}})\} be the optimal control trajectory. The control inputted into the system at the sample time k+\mathrm{1} is u_N^{*}(\mathrm{0},x_{k+\mathrm{1}}). Thus, the control law of MPC can be stated as

u_{RH}(x_k)=u_N^{*}(\mathrm{0},x_k),k=\mathrm{0},\mathrm{1},\mathrm{2},\dots ,\infty .

The closed-loop stability of the controlled system is shown in Lemma 1.

Lemma 1　For any x_{\mathrm{0}}\in X, if x_{\mathrm{0}} satisfies x_N^{*}(N,x_{\mathrm{0}})\in X_f and Assumption 1 is satisfied, it is guaranteed that x_{\mathrm{0}} will be steered to0by using the control law of MPC.

The proof can be found in Ref. [1].

Gradually approximation of the maximal terminal region

Using SVM classifier to estimate the terminal region, as described in Ref. [7], is considered as an effective approach for solving the approximating problem in MPC. However, the method in Ref. [7] is somewhat conservative. The reason is that, the obtained terminal region actually is the stabilizable region of using a predetermined linear feedback controller.

In this section, a novel method of computing a terminal region is proposed. Given the terminal cost function and an initial set of the maximal terminal region, both surrounding sequence and subsets sequence can be constructed by using one-step set expansion iteratively. When some conditions are satisfied, the iteration ends and the consequent set can be seen as the suboptimal estimation of the maximal terminal region.

Construction of surrounding set sequence

The proposed approach constructs the surrounding set sequence directly from conditions (C1) and (C2). The terminal region X_f can be defined as

(3)X_f:=\{x\in X|F(x)\ge F_{X_f}^{*}(x)\},

where F_{X_f}^{*}(x) is the solution of the following optimization problem:

(4)\begin{array}{l}\underset{u\in U}{\min }{F}_{X_f}(x)=q(x,u)+F(f(x,u)),\\ s.t. f(x,u)\in X_f.\end{array}

Remark 1　The definition of X_f has two meanings: (I) the optimization problem (4) has feasible solution, that is to say, \exists u\in U, s.t. f(x,u)\in X_f; (II) the minimum of the optimization problem satisfies that F_{X_f}^{*}(x)\le F(x).

Remark 2　From the definition of X_f, it is obvious that, the terminal region is essentially a positively invariant set of using the optimal control resulting from the optimization problem (4) when F(\cdot ) is given.

Remark 3　In Refs. [4,5,7], the linear feedback control is attached to the construction of X_f, and X_f is the stabilizable region of using the linear feedback controller. In Ref. [6], a saturated local control law was used. However, in this paper, there is no explicit control attached to the definition of X_f. Thus, the requirement on X_f is lower than that in Refs. [4-7] while it guarantees the stability of the controlled system.

From the definition of X_f, it cannot be determined whether a state point belongs to X_f. The difficulty arises that the X_f acts as the constraint in the optimization problem (4). To solve this problem, the method of using one-step set expansion iteratively is adopted.

The iterative approximating process begins with an initial estimation satisfying X_f^{\mathrm{0}}\supseteq X_f, and from which, a more accurate estimation X_f^{\mathrm{1}} is achieved. Repeatedly, the surrounding set sequence X_f^{\mathrm{2}},X_f^{\mathrm{3}},…,X_f^{+\infty }, which satisfies X_f^{\mathrm{0}}\supseteq X_f^{\mathrm{1}}\supseteq \cdot \cdot \cdot \supseteq X_f^j\supseteq X_f, can be estimated sequentially. And the convergence of this sequence is proved in previous section. The gradual approximation process can be depicted in Fig. 1, the specific of this approach can be found in the authors’ previous literature [12].

The initial surrounding set is determined by

(5)X_f^{\mathrm{0}}:=\{x\in X|F(x)\ge F_{X_f^{\mathrm{0}}}^{*}(x)\}.

The iterative construction of surrounding set sequence is enforced based on the following assumptions.

Assumption 2　A surrounding set of X_{f,\max }, denoted by X_f^{\mathrm{0}} and containing the origin, is initialized, where X_{f,\max } is defined as the largest terminal region.

Assumption 3　X_f^{\mathrm{0}} is a positively invariant set, that is to say, for any x\in X_f^{\mathrm{0}}, \exists u\in U, s.t.F(x)\ge q(x,u)+F(f(x,u)) and f(x,u)\in X_f^{\mathrm{0}}.

According to the known X_f^{\mathrm{0}}, another surrounding set of X_{f,\max }, denoted by X_f^{\mathrm{1}}, can be constructed as

(6)X_f^{\mathrm{1}}:=\{x\in X_f^{\mathrm{0}}|F(x)\ge F_{X_f^{\mathrm{1}}}^{*}(x)\},

where F_{X_f^{\mathrm{1}}}^{*}(x) is the minimum of

(7)\begin{array}{l}\underset{u\in U}{\min }{F}_{X_f^{\mathrm{1}}}(x)=q(x,u)+F(f(x,u)),\\ s.t. f(x,u)\in X_f^{\mathrm{0}}.\end{array}

As mentioned in Remark 1, the construction of X_f^{\mathrm{1}} contains two meanings: (I) for any x\in X_f^{\mathrm{1}}, \exists u\in U, s.t. f(x,u)\in X_f^{\mathrm{0}}; (II) the minimum of Eq. (7) satisfies F_{X_f^{\mathrm{0}}}^{*}(x)\le F(x). The constructions of X_f^j in sequel have the similar meanings. And the convergence of this sequence can be proved from the following theorem.

Theorem 1　If Assumption 2 and Assumption 3 are satisfied, the sequence \{X_f^j\}, similarly defined as X_f^{\mathrm{0}} and X_f^{\mathrm{1}}, will converge to X_{f,\max } when j goes to infinity.

Proof　This theorem is proved by contradiction.

Supposes there exists X_{spo}\supset X_{f,\max } to guarantee that X_f^j\to X_{spo}. Then, for \forall x\in X_{spo}, satisfies

F(x)\ge \underset{u\in U}{\min }\{q(x,u)+F(f(x,u))\},

where f(x,u)\in X_{spo}. Moreover, it is contractive to the fact that X_{f,\max } is the maximal set, which satisfies conditions (C1) and (C2) as listed in Assumption 1.

And vice versa, supposes there exists X_{spo}\subset X_{f,\max } to guarantee that X_f^j\to X_{spo}. Then, there exists a certain integer N that X_f^N\supseteq X_{f,\max } and X_f^{N+\mathrm{1}}\subset X_{f,\max }, for \forall x\in X_{f,\max }\backslash X_f^{N+\mathrm{1}}, we have

F(x)\ge \underset{u\in U}{\min }\{q(x,u)+F(f(x,u))\},

where f(x,u)\in X_{f,\max }\subseteq X_f^N. Because x\in X_f^N, then, x\in X_f^{N+\mathrm{1}}, which leads to paradox.

End proof.

From the previous description, we can find that the surrounding approximating approach will probably contain some unstable point by mistake. Then, we naturally proposed another gradually extending-approximation approach, which will be detailedly illustrated in the following subsection.

Construction of subsets sequence

Given F(\cdot ) and from conditions (C1) and (C2), the terminal region X_f can be defined as

(8)X_f:=\{x\in X|F(x)\ge F_{X_f}^{*}(x)\},

where F_{X_f}^{*}(x) is the minimum of the following optimization problem:

(9)\begin{array}{l}\underset{u\in U}{\min }{F}_{X_f}(x)=q(x,u)+F(f(x,u)),\\ s.t. f(x,u)\in X_f.\end{array}

Remind here that Remark 1 and Remark 2 are also consistent in the problem (9). From the definition of X_f, it cannot be determined whether a state point belongs to X_f. The difficulty lies in that, the X_f itself acts as the constraint in the optimization problem (9).

To avoid this problem, the method of using one-step set expansion iteratively is adopted. The gradual approximation process can be depicted in Fig. 2.

Assumption 4　A subset of X_{f,\max }, denoted by X_f^{\mathrm{0}} and containing the origin, is known, and X_{f,\max } is denoted as the largest terminal region.

Assumption 5　X_f^{\mathrm{0}} is a positively invariant set, that is to say, for any x\in X_f^{\mathrm{0}}, \exists u\in U, s.t. F(x)\ge q(x,u)+F(f(x,u))and f(x,u)\in X_f^{\mathrm{0}}.

Given X_f^{\mathrm{0}}, another subset of X_{f,\max }, denoted by X_f^{\mathrm{1}}, can be constructed as

(10)X_f^{\mathrm{1}}:=\{x\in X|F(x)\ge F_{X_f^{\mathrm{0}}}^{*}(x)\},

where F_{X_f^{\mathrm{0}}}^{*}(x) is the minimum of

(11)\begin{array}{l}\underset{u\in U}{\min }{F}_{X_f^{\mathrm{0}}}(x)=q(x,u)+F(f(x,u)),\\ s.t. f(x,u)\in X_f^{\mathrm{0}}.\end{array}

As mentioned in Remark 1, the construction of X_f^{\mathrm{1}} contains two meanings: (I) for any x\in X_f^{\mathrm{1}}, \exists u\in U, s.t. f(x,u)\in X_f^{\mathrm{0}}; (II) the minimum of Eq. (11) satisfies F_{X_f^{\mathrm{0}}}^{*}(x)\le F(x). The constructions of X_f^j in sequel have the similar meanings.

Lemma 2　If Assumption 3 is satisfied, there is X_f^{\mathrm{0}}\subseteq X_f^{\mathrm{1}}.

Proof　If Assumption 5 is satisfied, it is obvious that, for any x\in X_f^{\mathrm{0}}, \exists u\in U, s.t. F(x)\ge q(x,u)+F(f(x,u)) and f(x,u)\in X_f^{\mathrm{0}}. It follows that, F(x)\ge F_{X_f^{\mathrm{0}}}^{*}(x). From the construction of X_f^{\mathrm{1}}, we can know x\in X_f^{\mathrm{1}}, namely, X_f^{\mathrm{0}}\subseteq X_f^{\mathrm{1}}.

End proof.

Remark 4　From the construction of X_f^{\mathrm{1}}, it is obvious that, if Assumption 5 is satisfied, X_f^{\mathrm{1}} is a positively invariant set. We know that, for any x\in X_f^{\mathrm{1}}, \exists u\in U, s.t. F(x)\ge q(x,u)+F(f(x,u)) and f(x,u)\in X_f^{\mathrm{0}}. Because of X_f^{\mathrm{0}}\subseteq X_f^{\mathrm{1}} as shown in Lemma 2, we have f(x,u)\in X_f^{\mathrm{1}}.

Similarly, by replacing X_f^{\mathrm{0}} with X_f^{\mathrm{1}} in the constraint of Eq. (11), another subset, denoted by X_f^{\mathrm{2}}, can be obtained as follows:

(12)X_f^{\mathrm{2}}:=\{x\in X|F(x)\ge F_{X_f^{\mathrm{1}}}^{*}(x)\},

where F_{X_f^{\mathrm{1}}}^{*}(x) is the minimum of

(13)\begin{array}{l}\underset{u\in U}{\min }{F}_{X_f^{\mathrm{1}}}(x)=q(x,u)+F(f(x,u)),\\ s.t. f(x,u)\in X_f^{\mathrm{1}}.\end{array}

Repeatedly, X_f^j, j=\mathrm{3},\mathrm{4},\dots ,\infty, can be constructed as

(14)X_f^j:=\{x\in X|F(x)\ge F_{X_f^{j-\mathrm{1}}}^{*}(x)\},

where F_{X_f^{j-\mathrm{1}}}^{*}(x) is the minimum of

(15)\begin{array}{l}\underset{u\in U}{\min }{F}_{X_f^{j-\mathrm{1}}}(x)=q(x,u)+F(f(x,u)),\\ s.t. f(x,u)\in X_f^{j-\mathrm{1}}.\end{array}

This method of constructing X_f^j given X_f^{j-\mathrm{1}} is defined as one-step set expansion in this paper. By employing it iteratively, a subsets sequence of largest terminal region, denoted by \{X_f^j\}, j=\mathrm{1},\mathrm{2},\dots ,\infty, can be achieved.

Remark 5　Similar with Lemma 2 and Remark 3, any subset in this sequence is positively invariant and any two neighboring subsets satisfy X_f^{j-\mathrm{1}}\subseteq X_f^j.

As j increases, \{X_f^j\} will converge to a set, denoted by X_f^{+\infty }. Theorem 2 will show that, X_f^{+\infty } is equal to the largest terminal region.

Theorem 2　If Assumption 4 and Assumption 5 are satisfied, forX_f^j constructed in Eqs. (14) and (15), whenj goes to infinity, \{X_f^j\}will converge toX_{f,\max }.

Proof　This theorem is proved by contradiction.

(A) Assume that, there exists a set which is denoted by X_{spo} satisfying X_{spo}\subset X_{f,\max } and X_f^j\to X_{spo} when j\to +\infty. From Remark 5, we can know X_f^{\mathrm{0}}\subseteq X_{spo}. It is obvious that \mathbf{0}\in X_{spo} because of \mathbf{0}\in X_f^{\mathrm{0}} as shown in Assumption 4. It follows that, \mathbf{0}\notin X_{f,\max }\backslash X_{spo} and for any x\in X_{f,\max }\backslash X_{spo}, we have F(x)>\mathrm{0} since F(\cdot ) is positive definite. Define \xi as the minimum of \{F(x)|x\in X_{f,\max }\backslash X_{spo}\}, it is satisfied that, \xi >\mathrm{0}.

From the construction of X_f^j, we know that, for any x_{\mathrm{0}}\in X_{f,\max }\backslash X_{spo}, there exists no such a u\in U satisfying F(x_{\mathrm{0}})\ge q(x_{\mathrm{0}},u)+F(f(x_{\mathrm{0}},u)) and f(x_{\mathrm{0}},u)\in X_{spo} because of X_{spo}\subset X_{f,\max }. However, from conditions (C1) and (C2), we know that, \exists u(x_{\mathrm{0}})\in U, s.t. F(x_{\mathrm{0}})\ge q(x_{\mathrm{0}},u(x_{\mathrm{0}}))+F(x_{\mathrm{1}}) and x_{\mathrm{1}}\in X_{f,\max }, where x_{\mathrm{1}}=f(x_{\mathrm{0}},u(x_{\mathrm{0}})). It is obvious that, x_{\mathrm{1}}\notin X_{spo}. Thus, we have x_{\mathrm{1}}\in X_{f,\max }\backslash X_{spo}. Similarly, we can know, \exists u(x_{\mathrm{1}})\in U, s.t. F(x_{\mathrm{1}})\ge q(x_{\mathrm{1}},u(x_{\mathrm{1}}))+F(x_{\mathrm{2}}) and x_{\mathrm{2}}\in X_{f,\max }\backslash X_{spo}, where x_{\mathrm{2}}=f(x_{\mathrm{1}},u(x_{\mathrm{1}})), since x_{\mathrm{1}}\in X_{f,\max }\backslash X_{spo}.

Repeatly, for x_i\in X_{f,\max }\backslash X_{spo},\exists u(x_i)\in U, s.t. F(x_i)\ge q(x_i,u(x_i))+F(x_{i+\mathrm{1}}) and x_{i+\mathrm{1}}\in X_{f,\max }\backslash X_{spo}, where x_{i+\mathrm{1}}=f(x_i,u(x_i)), i=\mathrm{2},\mathrm{3},\dots ,\infty. It is clear that, F(x_i)\to \mathrm{0} when i\to \infty. We know that, for the minimum of \{F(x)|x\in X_{f,\max }\backslash X_{spo}\}, defined as \xi, there is a positive real number \delta satisfying . Since F(x_i)\to \mathrm{0} when i\to \infty, \exists N_{\delta }>\mathrm{0}, s.t. for any i\ge N_{\delta }, we have . Obviously, this is contradicted with that \xi is the minimum of \{F(x)|x\in X_{f,\max }\backslash X_{spo}\}.

(B) Similarly, assume that, there exists an X_{spo} satisfying X_{spo}\supset X_{f,\max } and X_f^j\to X_{spo} when j\to +\infty. For any x\in X_{spo}, we have that F(x)\ge \underset{u\in U}{\min }\{q(x,u)+F(f(x,u))\} and f(x,u)\in X_{spo}. Obviously, this is contradicted with that X_{f,\max } is the largest one satisfying conditions (C1) and (C2).

End proof.

Remark 6　In this paper, the largest terminal region means the positively invariant set satisfying conditions (C1) and (C2). However, conditions (C1) and (C2) are sufficient conditions to guarantee the stability of the controlled system, not the necessary conditions. There may be a set larger than X_{f,\max } and the stability of the controlled system can be guaranteed by using this set as the terminal region.

Remark 7　In the calculation of X_{f,\max }, it is impossible to keep iteration computation until j\to +\infty. When the iteration time goes to j=E (E is a positive integer), if X_f^E is equal to X_f^{E-\mathrm{1}} in principle, it can be deemed that \{X_f^j\} converges to X_f^E in rough. Hence, X_f^E can be taken as the terminal region and it is a good approximation to X_{f,\max }.

Remark 8　If the iteration time does not go to infinity, the obtained set may be just a large positively invariant subset of X_{f,\max }. This has no effect on the stability of the controlled system. The only negative influence is that its corresponding domain of attraction is smaller than that corresponding to X_{f,\max }.

Until now, it seems that we can choose any X_f^j in the subsets sequence as the terminal region. This is infeasible. Since X_f^j is not described in explicit expression, it cannot serve as the terminal constraint in the optimization problem (2) directly. Then, an estimated one described in explicit expression is needed. Due to the strong optimizing ability of SVM, SVM is exploited to separate each X_f^j from the state space.

Support vector machine

SVM is the youngest part in the statistical learning theory. It is an effective approach for pattern recognition. In SVM approach, the main aim is to obtain a function, which determines the decision boundary or hyperplane. This hyperplane optimally separates two classes of input data points.

Take the example of separating X into A and X\backslash A. For each x_i\in A, an additional variable y_i=+\mathrm{1} is introduced. Similarly, for each x_i\in X\backslash A, y_i=-\mathrm{1} is introduced. Define I^{+}:=\{i:y_i=+\mathrm{1}\} and I^{-}:=\{i:y_i=-\mathrm{1}\}, SVM will find a separating hyperplane, denoted by O(x):=w\cdot \varphi (x_i)+b=\mathrm{0}, between A and X\backslash A. Therefore, A can be estimated as \widehat{A}=\{x\in X|O(x)\ge \mathrm{0}\}, where O(x) is determined by solving the following problem:

(16)\begin{array}{l}\underset{\alpha }{\min }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \sum _i^{}{\displaystyle \sum _j^{}{\alpha }_i{\alpha }_j{y}_i{y}_{j}ker(x_i,x_j)}}-{\displaystyle \sum _i^{}{\alpha }_i},\\ s.t. {\displaystyle \sum _i^{}{\alpha }_i{y}_i}=\mathrm{0},\\ \mathrm{0}\le {\alpha }_i\le C,\begin{array}{l}\end{array}\forall i\in I^{+};\begin{array}{l}\end{array}{\alpha }_i\ge \mathrm{0},\begin{array}{l}\end{array}\forall i\in I^{-},\end{array}

where \mathit{ker}(\cdot ,\cdot ) denotes the kernel function and the Gaussian kernel as follows is adopted in this paper:

(17)ker(x,x_i)=\exp \left(-{\displaystyle \frac{{\Vert x-x_i\Vert }^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}}}\right),

with \sigma being the positive Gaussian kernel width.

When \{{\alpha }_i\} are computed out, some support vectors are chosen from \{x_i\} and the optimal hyperplane can be determined with these support vectors and their relevant weights. Denote P_s as the number of support vectors and X_s as the support vectors set, the optimal hyperplane is described as

(18)O(x)={\displaystyle \sum _{i=\mathrm{1}}^{P_s}{w}_i\cdot ker(x_i,x)}+b

where x_i\in X_s is the support vector and w_i={\alpha }_i{y}_i satisfying {\displaystyle {\sum }_{i=\mathrm{1}}^{P_s}{w}_i}=\mathrm{0} is the relevant weight.

There are many software packages of SVM available on internet. They can be downloaded and used directly. To save space, it is not introduced in detail in this paper. For more details, please refer to Refs. [10] and [11].

Estimating the approximating sequence by employing SVM

From Subsection 3.3, we know that, SVM find a separating hyperplane between \{x_i|i\in I^{+}\} and \{x_i|i\in I^{-}\}. This hyperplane is used to separate X into A and X\backslash A. All of \{x_i\} and their relevant \{y_i\} compose a set, named the training points set. This subsection will show how to achieve the training points set when estimating X_f^j and how to determine X_f^j when the separating hyperplane is known.

First, choose arbitrary points x_i\in X, i=\mathrm{1},\mathrm{2},\mathrm{...},P (P is the number of training points); then, assign y_i to each x_i by implementing the following procedure:

IF　(I) The following optimization problem has feasible solution

\begin{array}{l}\underset{u\in U}{\min }{F}_{X_f^j}(x_i)=q(x_i,u)+F(f(x_i,u)),\\ s.t. f(x_i,u)\in {\widehat{X}}_f^{j-\mathrm{1}}.\end{array}

\left(\mathrm{W}\mathrm{h}\mathrm{e}\mathrm{n}j=1,{\widehat{X}}_f^0=X_f^0.\right)

(II) Its minimum satisfies

F(x_i)\ge F_{X_f^j}^{*}(x_i).

THEN　 y_i=+\mathrm{1}

ELSE　 y_i=-\mathrm{1}

END IF

By implementing this procedure for every x_i, each y_i is known. Input \{x_i\} and \{y_i\} into SVM classifier, an optimal hyperplane O^j(x)=\mathrm{0} will be obtained. Therefore, the estimated set of X_f^j can be achieved as {\widehat{X}}_f^j=\{x\in X|O^j(x)\ge \mathrm{0}\}.

When {\widehat{X}}_f^j is known, the training points for separating X_f^{j+\mathrm{1}} from X can be computed by the similar procedure. By inputting them into SVM classifier, a hyperplane O^{j+\mathrm{1}}(x)=\mathrm{0} and an estimated set of X_f^{j+\mathrm{1}}, denoted by {\widehat{X}}_f^{j+\mathrm{1}}=\{x\in X|O^{j+\mathrm{1}}(x)\ge \mathrm{0}\}, will be obtained. Repeatedly, \{O^j(x)\},j=\mathrm{1},\mathrm{2},\dots ,\infty, and \left\{{\widehat{X}}_f^j\right\} can be achieved by the similar technology.

Estimating the terminal region

Section 3 showed how to achieve the subsets sequence by employing SVM. Theoretically, the larger the iteration time j is, the higher the precision of {\widehat{X}}_f^j approaching to X_{f,\max } is. However, it is impossible to keep computation until j\to +\infty. To avoid this problem, the iteration should be ended when some conditions are satisfied.

When j=E, if it is satisfied that, for x_i\in X_{s,E-\mathrm{1}}, i=\mathrm{1},\mathrm{2},\dots ,P_{s,E-\mathrm{1}}, there exists

(19){\displaystyle \sum _{i=\mathrm{1}}^{P_{s,E-\mathrm{1}}}\Vert O^E(x_i)-O^{E-\mathrm{1}}(x_i)\Vert }\le ϵP_{s,E-\mathrm{1}},

and we can find that {\widehat{X}}_f^E is equal to {\widehat{X}}_f^{E-\mathrm{1}} in principle and {\widehat{X}}_f^j converges to {\widehat{X}}_f^E. In Eq. (19), X_{s,E-\mathrm{1}} is the support vectors set at j=E-\mathrm{1}, P_{s,E-\mathrm{1}} is the number of support vectors and ϵ is a tunable threshold. The smaller ϵ is, the higher the precision of {\widehat{X}}_f^E approximating to X_{f,\max } is. Finally, {\widehat{X}}_f^E is used to serve as the terminal region.

Remark 9　Here, we used the information that, in SVM classifier, the hyperplanes are only determined on the support vectors.

Now, the concrete algorithm of estimating the largest terminal region is displayed as follows.

Step 1　Set the number of training points P used in SVM and the tunable threshold ϵ.

Step 2　For j=\mathrm{1},\mathrm{2},\dots ,\infty, use SVM to achieve the optimal hyperplane O^j(x)=\mathrm{0} and the estimated set of X_f^j, denoted by {\widehat{X}}_f^j.

Substep 1　Choose arbitrary points x_i\in X, i=\mathrm{1},\mathrm{2},\mathrm{...},P.

Substep 2　Assign y_i to each x_i by implementing the procedure in Subsection 3.4.

Substep 3　Input \{x_i,y_i\} into the SVM. An optimal hyperplane O^j(x)=\mathrm{0} will be obtained and X_f^j can be approximated by {\widehat{X}}_f^j=\{x\in X|O^j(x)\ge \mathrm{0}\}, where

O^j(x)={\displaystyle \sum _{i=\mathrm{1}}^{P_{s,j}}{w}_i\cdot ker(x_i,x)}+b_j,

with P_{s,j} denoting the number of support vectors, x_i being the support vector, w_i denoting its relevant weight, and b_j denoting the classifier threshold.

Step 3　Check the iteration status. When j=E, if inequality (14) is satisfied, end iteration and take {\widehat{X}}_f^E as the largest terminal region.

Remark 10　It is obvious that, {\widehat{X}}_f^j is achieved one by one. Namely, {\widehat{X}}_f^j can only be achieved when {\widehat{X}}_f^{j-\mathrm{1}} is known.

Simulation experiment

In this section, we implemented both surrounding and subset approximating approach in simulated model, and the comparison to representative alternatives are shown in results. The model is a discrete-time realization of the continuous-time system used in Refs. [4,7]:

\begin{array}{l}\left[\begin{array}{l}{x}_{\mathrm{1}}(k+\mathrm{1})\\ x_{\mathrm{2}}(k+\mathrm{1})\end{array}\right]=\left[\begin{array}{ll}\mathrm{1}& T\\ T& \mathrm{1}\end{array}\right]\left[\begin{array}{l}{x}_{\mathrm{1}}(k)\\ x_{\mathrm{2}}(k)\end{array}\right]+\left[\begin{array}{l}T\mu \\ T\mu \end{array}\right]u(k)\\ +\left[\begin{array}{cc}T(\mathrm{1}-\mu )& \mathrm{0}\\ \mathrm{0}& -\mathrm{4}T(\mathrm{1}-\mu )\end{array}\right]\left[\begin{array}{l}{x}_{\mathrm{1}}(k)\\ x_{\mathrm{2}}(k)\end{array}\right]u(k),\end{array}

where \mu =\mathrm{0.5}, T=\mathrm{0.1}\mathbf{s}, and the state constraint and control constraint are X=\{x|{\Vert x\Vert }_{\mathrm{1}}\le \mathrm{4}\}, U=\{u|\left|u\right|\le \mathrm{2}\}, respectively.

The stage cost is q(x,u)=x^{\mathrm{T}}Qx+u^{\mathrm{T}}Ru, where Q=\mathrm{0.5}I and R=\mathrm{1}. The terminal cost is chosen as F(x)=x^{\mathrm{T}}Gx, where G=[\mathrm{1107.356}\mathrm{857.231};\mathrm{857.231}\mathrm{1107.356}] and X_f^{\mathrm{0}} is given as the terminal region in Ref. [4], which is

X_f^{\mathrm{0}}=\left\{x\in X|x^{\mathrm{T}}\left[\begin{array}{ll}\mathrm{16.5926}& \mathrm{11.5926}\\ \mathrm{11.5926}& \mathrm{16.5926}\end{array}\right]x\le \mathrm{0.7}\right\}

.

To estimate each X_f^j, 4000 training points are generated, and set ϵ=\mathrm{2.5}. From the simulating process we find that in the surrounding approximating approach, when j is iterated to 15, X_f^{\mathrm{15}} and X_f^{\mathrm{14}} are near enough. There only exists minor difference between X_f^{\mathrm{15}} and X_f^j, then,

{\displaystyle \sum _{i=\mathrm{1}}^{P_{s,\mathrm{14}}}\Vert O^{\mathrm{15}}(x_i)-O^{\mathrm{14}}(x_i)\Vert }\le ϵP_{s,\mathrm{14}},

where x_i\in X_{\mathrm{s},\mathrm{14}}, X_{\mathrm{s},\mathrm{14}} is the support vectors set, and P_{s,\mathrm{14}} is the number of support vectors at j=\mathrm{14}. For \forall j>\mathrm{15}, so we can choose X_f^{\mathrm{15}} as the terminal estimation of X_{f,\max }. The approximating process can be well identified in Fig. 3. On the other hand, if the subsets approximating approach is adopted, if j is as large as 22, the following iterations various slowly enough. It can be rearranged as

{\displaystyle \sum _{i=\mathrm{1}}^{P_{s,\mathrm{21}}}\Vert O^{\mathrm{22}}(x_i)-O^{\mathrm{21}}(x_i)\Vert }\le ϵP_{s,\mathrm{21}},

where x_i\in X_{\mathrm{s},\mathrm{21}}, X_{\mathrm{s},\mathrm{21}} is the support vectors set, and P_{s,\mathrm{21}} is the number of support vectors at j=\mathrm{21}. Thus, {\widehat{X}}_f^{\mathrm{22}} can be seen as the final estimation of X_{f,\max }, and the approximating process is well depicted in Fig. 4.

In Figs. 3 and 4, the blue circles, which stem from the method that presented in Ref. [4], are selected as the initial set of surrounding and subset approaches, respectively. We can find that the converged sets in both approaches are nearly the same, as depicted by black solid lines, so in the following illustration, the estimation results of surrounding method is neglected for clarity.

Figure 5 depicts the convergent trajectory of some state points, which are selected to implement the comparison between representative alternatives and the proposed approach. These state points are \left(\begin{array}{ll}-\mathrm{4}& \mathrm{2}\end{array}\right), \left(\begin{array}{ll}-\mathrm{4}& \mathrm{4}\end{array}\right), \left(\begin{array}{ll}-\mathrm{3}& \mathrm{1.5}\end{array}\right), \left(\begin{array}{ll}-\mathrm{2}& \mathrm{3.5}\end{array}\right), \left(\begin{array}{ll}\mathrm{0.15}& -\mathrm{4}\end{array}\right), \left(\begin{array}{ll}\mathrm{0.3}& \mathrm{0}\end{array}\right) and \left(\begin{array}{ll}\mathrm{3.7}& -\mathrm{0.5}\end{array}\right), as denoted by “*” in Fig. 5.

In Fig. 5, the blue ellipsoid is the terminal region in Ref. [4] and the region encompassed by black dash lines is the result in Ref. [7]. The region encompassed by black solid lines is the terminal region in this paper. From this figure, we can find that all the selected states, which obviously lie in the domain of attraction, are out of the estimating regions of both Refs. [4] and [7], but they are all included in the corresponding results of this paper. The terminal region in this paper contains the solution that obtained in Ref. [4], but not entirely contains the result in Ref. [7] although it is much larger than that in Ref. [7]. The reason is that, the terminal region in this paper is the largest one satisfying conditions (C1) and (C2). However, conditions (C1) and (C2) are just the sufficient conditions to guarantee the stability of the controlled system, not the necessary conditions as shown in Remark 9. The red solid lines denote the closed-loop trajectories of the selected points. Note that, with the same sampling interval and prediction horizon as those in this paper, these points are not in the regions of attraction of MPC in Refs. [4] and [7]. However, they can be leaded to the origin by using the control law of MPC in this paper.

Conclusion

Given the terminal cost, a sequence of surrounding sets and subsets of the maximal terminal region is extracted from state space iteratively by employing SVM classifier. The convergence of these approximation processes are proved in analytical way. Simulations also show that the convergent region is better approximation to the maximal terminal region, when compared to some mature representative alternatives.

References

[1] Mayne D Q, Rawlings J B, Rao C V, Scokaert P O M. Constrained model predictive control: Stability and optimality. Automatica, 2000, 36(6): 789-814

[2] González A H, Odloak D. Enlarging the domain of attraction of stable MPC controllers, maintaining the output performance. Automatica, 2009, 45(4): 1080-1085

[3] Magni L, De Nicolao G, Magnani L, Scattolini R. A stabilizing model based predictive control algorithm for nonlinear systems. Automatica, 2001, 37(9): 1351-1362

[4] Chen H, Allgower F. A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability. Automatica, 1998, 34(10): 1205-1217

[5] Cannon M, Deshmukh V, Kouvaritakis B. Nonlinear model predictive control with polytopic invariant sets. Automatica, 2003, 39(8): 1487-1494

[6] De Doná J A, Seron M M, Mayne D Q, Goodwin G C. Enlarged terminal sets guaranteeing stability of receding horizon control. Systems & Control Letters, 2002, 47(1): 57-63

[7] Ong C J, Sui D, Gilbert E G. Enlarging the terminal region of nonlinear model predictive control using the support vector machine method. Automatica, 2006, 42(6): 1011-1016

[8] Limon D, Alamo T, Salas F, Camacho E F. On the stability of constrained MPC without terminal constraint. IEEE Transactions on Automatic Control, 2006, 51(5): 832-836

[9] Limon D, Alamo T, Camacho E F. Enlarging the domain of attraction of MPC controllers. Automatica, 2005, 41(4): 629-635

[10] Burges C C. A tutorial on support vector machines for pattern recognition. Data Mining and Knowledge Discovery, 1998, 2(2): 121-167

[11] Vapnik V. The Nature of Statistical Learning Theory. New York: Springer, 1995

[12] Wang Y F, Sun F C, Zhang Y A, Liu H P, Min H B. Getting a suitable terminal cost and maximizing the terminal region for MPC. Mathematical Problems in Engineering, 2010, 2010: 1-16

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg