# Getting a Suitable Terminal Cost and Maximizing the Terminal Region for MPC

Authors: Wang Ya-feng; Sun Fu-chun; Zhang You-an; Liu Hua-ping; Min Haibo
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

The model predictive control (MPC) subject to control and state constraint is studied. Given a terminal cost, a terminal region is obtained through iterative estimation by using support vector machine (SVM). It is proved that the obtained terminal region is the largest terminal region when the terminal cost is given. The relationships between terminal cost and terminal region and between terminal cost and total cost are discussed, respectively. Based on these relationships, a simple method to get a suitable terminal cost is proposed and it can be adjusted according to our need. Finally, some experiment results are presented.


---

## Body

## 1. Introduction

Model predictive control also known as receding horizon control has become quite popular recently. The key advantage is its ability to handle control and state constraints. It was pointed out in [1] that what MPC solves is the standard optimal control (SOC) problem except that it uses the finite-horizon optimization (some paper used quasi-infinite horizon optimization, such as [2, 3]) to replace the infinite-horizon optimization of SOC and the control is computed online.Generally, to say that an MPC is good or bad, it is contrasted with SOC from two aspects: the domain of attraction and the total cost (the wasted performance index from initial time to infinity). If an MPC has a larger domain of attraction and for any initial state point the total cost is fewer than those of another MPC, it is considered as a better MPC. In MPC, there are three factors playing important roles in its performances on the two aspects as just mentioned: the prediction horizon, the terminal region, and the terminal cost. As known to all, lengthening the prediction horizon, the domain of attraction will be enlarged and the total cost will decrease, but the computation burden of online optimization will increase. Recently, many attentions have concentrated on the latter two factors:(I) how to get a large terminal region? (II) and how to get a suitable terminal cost? Here, only some typical papers are listed. Chen and Allgöwer [2] presented a terminal cost by using SOC method based on the linearized model of system and took an ellipsoidal set in which the state can be drived to the equilibrium point by linear feedback control as the terminal region. Cannon et al. [4] used a polytopic set to replace the ellipsoidal set. De Doná et al. [5] took the stabilizable region of using saturated control as the terminal region. Ong et al. [6] got a very large stabilizable region and a terminal cost of using linear feedback control via support vector machine. Limon et al. [7] proved that, for MPC without terminal state constraint in its on-line optimization, the terminal region will be enlarged by weighting the terminal cost. Most of these papers have a common shortage; the terminal region is computed under a precondition; some explicit controller was given in advance, like linear feedback controller and saturated controller. So, the computed terminal region is somewhat conservative, in other words, it is not the largest one.In this paper, a novel method is proposed to get a terminal state region. Given a terminal cost, a set sequence is obtained by using one-step set contraction iteratively. It is proved that, when the iteration time goes to infinity, this set sequence will converge to the maximal terminal region. In this sequence, each set is estimated by using support vector machine (SVM, see [7, 8] for details). Next, the relationships between terminal cost and terminal region and between terminal cost and total cost are discussed, respectively. Then, a simple method to get a suitable terminal cost according to our need is given. Finally, some experiment results are presented.


## 2. The Relationship between SOC and MPC

As mentioned, MPC is an approximation to SOC, and SOC is the counterpoint to evaluate MPC. Here, the study on MPC begins with the comprehension of SOC. Consider the discrete-time system as follows:(2.1)xk+1=f(xk,uk),
where xk∈Rn,uk∈Rm are the state and the input of the system at sampling time t=k, respectively.  xk+1∈Rn is the successor state and the mapping f:Rn+m↦Rn with f(0,0)=0 is known. The system is subject to constraints on both state and control action, and they are given by xk∈X,uk∈U, where X is a closed set and U a compact set, both of them containing the origin. The control objective is usually to steer the state to the origin.The optimization problemP∞(x0) of SOC at the initial state x0 can be stated as follows:
(2.2)minu(i,x0)∈UJ∞(u,x0)=∑i=0∞q(x(i,x0),u(i,x0)),s.t.x(i+1,x0)=f(x(i,x0),u(i,x0)),x(i+1,x0)∈X,u(i,x0)∈U,
where x(0,x0)=x0, q(x,u) is the stage cost and its form is chosen as q(x,u)=xTQx+uTRu in which Q,R are positive definites.It is well known that the stability is guaranteed ifP∞(x0) has feasible solution. Here, assume that the solution of P∞(x0) is known if it has the following. J∞*(x0) is defined as the solution, u∞*(i,x0), i=1,2,…,∞, the optimal control trajectory and x∞*(i,x0), the corresponding state trajectory. As is well known, J∞*(x0) is the total cost of using u∞*(i,x0) to drive x0 to 0.But solvingP∞(x0) is not an easy job, especially when f is nonlinear. To avoid this problem, the infinite-horizon optimization of SOC can be approximated by the finite-horizon optimization of MPC (As mentioned, quasi-infinite horizon optimization was used in [2, 3]. For convenience, we consider it to belong to the frame of finite-horizon optimization).Similarly, the optimization problemPN(x0) of MPC at the initial state x0 can be stated as
(2.3)minu(i,x0)∈UJN(u,x0)=∑i=0N-1q(x(i,x0),u(i,x0))+F(x(N,x0)),s.t.x(i+1,x0)=f(x(i,x0),u(i,x0)),x(i+1,x0)∈X,u(i,x0)∈U,x(N,x0)∈Xf,
where N is prediction horizon, Xf is terminal region, and F is terminal cost satisfying F(0)=0 and F(x)≥α(∥x∥), (the mapping α:R+↦R+ satisfying α(0)=0 is continuous and strictly increasing, where R+={x∈R∣x>0}).There exist many optimization algorithms to compute the solution ofPN(x0). Let JN*(x0) be the solution, uN*(x0)={uN*(0,x0),…,uN*(N-1,x0)},xN*(x0)={xN*(0,x0),…,xN*(N-1,x0)} be the optimal control trajectory and corresponding predicted state trajectory of PN(x0), respectively.At sampling timet=0, uN*(0,x0) is inputted to the system. At the next sampling time t=1, x1=f(x0,uN*(0,x0)) is outputted, and the control input uN*(0,x1) can be computed by solving the optimization problem PN(x1). By repeating this procedure, two trajectories can be obtained:  uRH(x0)={uN*(0,x0),uN*(0,x1),…},xRH(x0)={x0,x1,…}. Here, for convenience, they are called as the receding horizon control trajectory and the receding horizon state trajectory of MPC with PN(x0), respectively.The introduction ofXf and F in (2.3) is to guarantee the closed loop stability of system.Lemma 2.1.
DefineΓN:={x0∈X∣xN*(N,x0)∈Xf}. For any x0∈ΓN if Xf and F satisfy two conditions as follows:
(C1)
F being a Lyapunov function, more strictly, for any x∈Xf, there exists
(2.4)F(x)≥minu∈U{q(x,u)+F(f(x,u))},
(C2)
Xf being an invariant set. In other words, one has f(x,u)∈Xf. Here, u is the control in (C1).It is guaranteed that,x will be led to 0 by using uRH(x0).  ΓN is called the domain of attraction.The proof can be found in [1].Obviously, the optimal choice ofF(x(N,x0)) is J∞(x(N,x0)). The total cost with this choice is J∞(x0), namely it is the least one. But, as mentioned, J∞ can not be obtained generally, so other ways should be found to get F. Certainly, the closer it approaches to J∞, the better it is. When F is given, Xf can be found to satisfy the conditions (C1) and (C2). There are many Xf satisfying the conditions and different Xf will be obtained by different methods. But, to make the domain of attraction the largest, it is wished that the largest Xf can be obtained. Here, define Xf,max as the largest terminal region when F is given, and in the next section, a novel method to get Xf,max will be proposed.


## 3. Maximizing the Terminal Region for MPC

Until now, there exist many methods to constructXf. As mentioned, these methods have a common basic idea: some controller like linear feedback controller or saturation controller is given in advance, then a stabilizable domain of using this controller is computed and works as the terminal region of MPC. It is obvious that this kind of construction is somewhat conservative and the Xf computed by using this method does not approximate to Xf,max to the largest extent.In this paper, a novel method is proposed in whichXf is constructed directly from conditions (C1) and (C2).


### 3.1. Approximating the Largest Terminal Region Asymptotically

DefineXf as
(3.1)Xf:={x∈X∣F(x)≥FXf*(x)},
where FXf*(x)=minFXf(x) is the solution of the following optimization problem:
(3.2)minu∈UFXf(x)=q(x,u)+F(f(x,u)),s.t.f(x,u)∈Xf.Obviously, for anx∈X, it cannot be decided whether x belongs to Xf from (3.1) and (3.2) when Xf is unknown. The difficulty is that the state constraint in the optimization (3.2) uses the Xf itself. To avoid it, the method of asymptotic approximation is adopted. Firstly, an initial set Xf0 which can be obtained by the following discriminant is given:
(3.3)Xf0:={x∈X∣F(x)≥FXf0*(x)},
where FXf0*(x)  is the solution of
(3.4)minu∈UFXf0(x)=q(x,u)+F(f(x,u)),s.t.f(x,u)∈X.Then, usingXf0 instead of X in the state constraint in optimization of (3.4), Xf1 can be obtained. One by one, Xf2，Xf3,… will be obtained. The whole procedure can be pictured as in Figure 1.Figure 1
The approximation ofXf,max.In Figure1, Xf1 is defined as
(3.5)Xf1:={x∈Xf0∣F(x)≥FXf1*(x)},
where FXf1*(x) is the solution of
(3.6)minu∈UFXf1(x)=q(x,u)+F(f(x,u)),s.t.f(x,u)∈Xf0.Similarly,Xfj can be defined as：(3.7)Xfj:={x∈Xfj-1∣F(x)≥FXfj*(x)},
and FXfj*(x) is the solution of
(3.8)minu∈UFXfj(x)=q(x,u)+F(f(x,u)),s.t.f(x,u)∈Xfj-1.Obviously, as Figure1 shows, there exists Xf0⊇Xf1⊇⋯⊇Xfj⊇Xf. As j increases, Xfj will converge to a set denoted by Xf+∞. But, whether is Xf+∞ the terminal region we want? Theorem 3.1 provides the answer to this question.Theorem 3.1.
ForXfj defined in (3.7), when j goes to infinity, Xfj will converge to Xf,max, namely, as j→+∞, one has Xfj→Xf,max.This theorem is proved by contradiction.Proof.
(A) Assume that there existsXspo which satisfies Xspo⊃Xf,max and Xfj→Xspo, as j→+∞ then for any x∈Xspo, we have F(x)≥minu∈U{q(x,u)+F(f(x,u))} and f(x,u)∈Xspo. Obviously this is contradicted with Xf,max being the largest one satisfying (C1) and (C2).
(B)
Similarly, assume that there existsXspo which satisfies Xspo⊂Xf,max and Xfj→Xspo, as j→+∞, then there exists N with 0≤N<+∞ satisfying XfN⊇Xf,max and Xf,max∖XfN+1≠Φ, where Φ denotes the empty set. Choose any x∈Xf,max∖XfN+1, it is obviouse that, x satisfies F(x)≥minu∈U{q(x,u)+F(f(x,u))} and f(x,u)∈Xf,max⊆XfN. On the other hand, we know that x∈XfN, so x meets the conditions in the definition of XfN+1 and we have x∈XfN+1. This is contradicted with x∈Xf,max∖XfN+1.Remark 3.2.
Generally, in the computation ofXf,max, it is impossible to keep computation until j→+∞. So, when the iteration goes to j=N, if XfN is equal to XfN-1 in principle, XfN can be taken as the terminal region we want.Remark 3.3.
The terminal region computed through the method in Remark3.2 is not Xf,max itself, but its enclosing set. Then, the corresponding domain of attraction may include some points which should not be in the real domain of attraction. To avoid this problem, the outspread skill from smaller region to larger one can be used to turn the contraction skill from larger region to smaller one in this paper. The concrete algorithm is not presented here, just the general idea is stated; giving a known subset of Xf,max, denoted by Xf,0 in advance and using Xf,0 to serve as the state constraint in (3.4), then a larger region Xf,1 will be computed. By the same procedure as presented in Section 3.1, a terminal region which is a subset of Xf,max will be gotten.To obtainXfj is not an easy job. The only tool is statistical learning method. Here, the SVM is used.


### 3.2. Support Vector Machine

SVM (see [8, 9]) is the youngest part of statistical learning theory. It is an effective approach for pattern recognition. In SVM approach, the main aim of an SVM classifier is obtaining a function, which determines the decision boundary or hyperplane. This hyperplane optimally separates two classes of input data points.Take the example of separatingX into A and X∖A. For each xi∈A, an additional variable yi=+1 is introduced. Similarly, for each xi∈X∖A, yi=-1 is introduced. Define I+:={i:yi=+1}, I-:={i:yi=-1}. SVM is used to find a separating hyperplane O(x):=w·ϕ(xi)+b=0, between A and X∖A. Then, we obtain an estimated set of A, Â={x∈X∣O(x)≥0}.  O(x) can be obtained by solving the following problem:
(3.9)minα12∑i∑jαiαjyiyjker(xi,xj)-∑iαi,s.t.∑iαiyi=0,0≤αi≤C,∀i∈I+,αi≥0,∀i∈I-,
where ker denotes the kernel function.In this paper, the following Gaussian kernel is used:(3.10)ker(x,xi)=exp(-∥x-xi∥22σ2),
with σ being the positive Gaussian kernel width.There are many software packages of SVM available on internet. They can be downloaded and used directly. By using SVM, the support vectors are extracted from{xi} and their relevant weights are exported. Denote Ps as the number of support vectors and Xs as the support vectors set, the optimal hyperplane is described as follows:
(3.11)O(x)=∑i=1Pswi·ker(xi,x)+b,
where xi∈Xs is a support vector, and wi=αiyi satifying ∑i=1Pswi=0 is the relevant weight.


### 3.3. Estimating the Largest Terminal Region

In SVM classifier, the training data is inputted and the hyperplane will be outputted. To us, the training of data is the only job.Take the separation ofXf0 from X as an example. Firstly, choose arbitrary points xi∈X, i=1,2,…,Np (Np is the number of training points) then decide the value of yi corresponding to xi by using the following procedureIfF(xi)≥FXf0*(xi),yi=+1;  else yi=-1endif.When all theyi for all the xi are gotten, they can be packed to constitute the training data. Then, by inputting the training data into SVM classifier, an optimal hyperplane O0(x)=0 and an estimated set of Xf0, X̂f0={x∈X∣O0(x)≥0} will be obtained.WhenX̂f0 is known, the training data for separating Xf1 from Xf0 can be known by the similar procedure. By inputting them into SVM classifier, a hyperplane O1(x)=0 and an estimated set of Xf1, X̂f1={x∈X∣O1(x)≥0} will be gotten.Repeatedly, a hyperplane seriesOj(x)=0,j=1,2,… will be obtained. When j=N if it is satisfied that, for xi∈Xsup,N-1, i=1,2,…,Nsup,N-1, there exists
(3.12)∑i=1Nsup,N-1∥ON(xi)-ON-1(xi)∥≤εNsup,N-1,
it is deemed that X̂fN is equal to X̂fN-1 in principle and X̂fN is taken as the final estimation of Xf,max. Where Xsup,N-1 is the support vectors set at j=N-1, Nsup,N-1 is the number of support vectors and ε is a tunable threshold. The smaller it is, the higher the precision of X̂fN approximating to Xf,max is.Remark 3.4.
Here, we used the information that, in SVM classifier, the hyperplanes are decided just on the support vectors.


## 3.1. Approximating the Largest Terminal Region Asymptotically

DefineXf as
(3.1)Xf:={x∈X∣F(x)≥FXf*(x)},
where FXf*(x)=minFXf(x) is the solution of the following optimization problem:
(3.2)minu∈UFXf(x)=q(x,u)+F(f(x,u)),s.t.f(x,u)∈Xf.Obviously, for anx∈X, it cannot be decided whether x belongs to Xf from (3.1) and (3.2) when Xf is unknown. The difficulty is that the state constraint in the optimization (3.2) uses the Xf itself. To avoid it, the method of asymptotic approximation is adopted. Firstly, an initial set Xf0 which can be obtained by the following discriminant is given:
(3.3)Xf0:={x∈X∣F(x)≥FXf0*(x)},
where FXf0*(x)  is the solution of
(3.4)minu∈UFXf0(x)=q(x,u)+F(f(x,u)),s.t.f(x,u)∈X.Then, usingXf0 instead of X in the state constraint in optimization of (3.4), Xf1 can be obtained. One by one, Xf2，Xf3,… will be obtained. The whole procedure can be pictured as in Figure 1.Figure 1
The approximation ofXf,max.In Figure1, Xf1 is defined as
(3.5)Xf1:={x∈Xf0∣F(x)≥FXf1*(x)},
where FXf1*(x) is the solution of
(3.6)minu∈UFXf1(x)=q(x,u)+F(f(x,u)),s.t.f(x,u)∈Xf0.Similarly,Xfj can be defined as：(3.7)Xfj:={x∈Xfj-1∣F(x)≥FXfj*(x)},
and FXfj*(x) is the solution of
(3.8)minu∈UFXfj(x)=q(x,u)+F(f(x,u)),s.t.f(x,u)∈Xfj-1.Obviously, as Figure1 shows, there exists Xf0⊇Xf1⊇⋯⊇Xfj⊇Xf. As j increases, Xfj will converge to a set denoted by Xf+∞. But, whether is Xf+∞ the terminal region we want? Theorem 3.1 provides the answer to this question.Theorem 3.1.
ForXfj defined in (3.7), when j goes to infinity, Xfj will converge to Xf,max, namely, as j→+∞, one has Xfj→Xf,max.This theorem is proved by contradiction.Proof.
(A) Assume that there existsXspo which satisfies Xspo⊃Xf,max and Xfj→Xspo, as j→+∞ then for any x∈Xspo, we have F(x)≥minu∈U{q(x,u)+F(f(x,u))} and f(x,u)∈Xspo. Obviously this is contradicted with Xf,max being the largest one satisfying (C1) and (C2).
(B)
Similarly, assume that there existsXspo which satisfies Xspo⊂Xf,max and Xfj→Xspo, as j→+∞, then there exists N with 0≤N<+∞ satisfying XfN⊇Xf,max and Xf,max∖XfN+1≠Φ, where Φ denotes the empty set. Choose any x∈Xf,max∖XfN+1, it is obviouse that, x satisfies F(x)≥minu∈U{q(x,u)+F(f(x,u))} and f(x,u)∈Xf,max⊆XfN. On the other hand, we know that x∈XfN, so x meets the conditions in the definition of XfN+1 and we have x∈XfN+1. This is contradicted with x∈Xf,max∖XfN+1.Remark 3.2.
Generally, in the computation ofXf,max, it is impossible to keep computation until j→+∞. So, when the iteration goes to j=N, if XfN is equal to XfN-1 in principle, XfN can be taken as the terminal region we want.Remark 3.3.
The terminal region computed through the method in Remark3.2 is not Xf,max itself, but its enclosing set. Then, the corresponding domain of attraction may include some points which should not be in the real domain of attraction. To avoid this problem, the outspread skill from smaller region to larger one can be used to turn the contraction skill from larger region to smaller one in this paper. The concrete algorithm is not presented here, just the general idea is stated; giving a known subset of Xf,max, denoted by Xf,0 in advance and using Xf,0 to serve as the state constraint in (3.4), then a larger region Xf,1 will be computed. By the same procedure as presented in Section 3.1, a terminal region which is a subset of Xf,max will be gotten.To obtainXfj is not an easy job. The only tool is statistical learning method. Here, the SVM is used.


## 3.2. Support Vector Machine

SVM (see [8, 9]) is the youngest part of statistical learning theory. It is an effective approach for pattern recognition. In SVM approach, the main aim of an SVM classifier is obtaining a function, which determines the decision boundary or hyperplane. This hyperplane optimally separates two classes of input data points.Take the example of separatingX into A and X∖A. For each xi∈A, an additional variable yi=+1 is introduced. Similarly, for each xi∈X∖A, yi=-1 is introduced. Define I+:={i:yi=+1}, I-:={i:yi=-1}. SVM is used to find a separating hyperplane O(x):=w·ϕ(xi)+b=0, between A and X∖A. Then, we obtain an estimated set of A, Â={x∈X∣O(x)≥0}.  O(x) can be obtained by solving the following problem:
(3.9)minα12∑i∑jαiαjyiyjker(xi,xj)-∑iαi,s.t.∑iαiyi=0,0≤αi≤C,∀i∈I+,αi≥0,∀i∈I-,
where ker denotes the kernel function.In this paper, the following Gaussian kernel is used:(3.10)ker(x,xi)=exp(-∥x-xi∥22σ2),
with σ being the positive Gaussian kernel width.There are many software packages of SVM available on internet. They can be downloaded and used directly. By using SVM, the support vectors are extracted from{xi} and their relevant weights are exported. Denote Ps as the number of support vectors and Xs as the support vectors set, the optimal hyperplane is described as follows:
(3.11)O(x)=∑i=1Pswi·ker(xi,x)+b,
where xi∈Xs is a support vector, and wi=αiyi satifying ∑i=1Pswi=0 is the relevant weight.


## 3.3. Estimating the Largest Terminal Region

In SVM classifier, the training data is inputted and the hyperplane will be outputted. To us, the training of data is the only job.Take the separation ofXf0 from X as an example. Firstly, choose arbitrary points xi∈X, i=1,2,…,Np (Np is the number of training points) then decide the value of yi corresponding to xi by using the following procedureIfF(xi)≥FXf0*(xi),yi=+1;  else yi=-1endif.When all theyi for all the xi are gotten, they can be packed to constitute the training data. Then, by inputting the training data into SVM classifier, an optimal hyperplane O0(x)=0 and an estimated set of Xf0, X̂f0={x∈X∣O0(x)≥0} will be obtained.WhenX̂f0 is known, the training data for separating Xf1 from Xf0 can be known by the similar procedure. By inputting them into SVM classifier, a hyperplane O1(x)=0 and an estimated set of Xf1, X̂f1={x∈X∣O1(x)≥0} will be gotten.Repeatedly, a hyperplane seriesOj(x)=0,j=1,2,… will be obtained. When j=N if it is satisfied that, for xi∈Xsup,N-1, i=1,2,…,Nsup,N-1, there exists
(3.12)∑i=1Nsup,N-1∥ON(xi)-ON-1(xi)∥≤εNsup,N-1,
it is deemed that X̂fN is equal to X̂fN-1 in principle and X̂fN is taken as the final estimation of Xf,max. Where Xsup,N-1 is the support vectors set at j=N-1, Nsup,N-1 is the number of support vectors and ε is a tunable threshold. The smaller it is, the higher the precision of X̂fN approximating to Xf,max is.Remark 3.4.
Here, we used the information that, in SVM classifier, the hyperplanes are decided just on the support vectors.


## 4. Choosing an Appropriate Terminal Cost

In the previous chapter, a method to maximize the terminal region was proposed, but the method has a premise: the terminal cost is given in advance. In this chapter, how to get a terminal cost will be shown. Before this, some properties of terminal cost will be analyzed.


### 4.1. Weighting Terminal Cost, the Domain of Attraction Will Be Enlarged

From conditions (C1) and (C2), it is known that the terminal region is based on the choice of the terminal cost. We want to know what the relationship between them is and this relationship will give us what messages when we choose terminal cost. Theorem4.1 will give us the answer.DenoteXf,max,Xλ,max as the terminal regions of PN(x0) and PN,λ(x0) with a weighted terminal cost λF, λ≥1, respectively. Limon et al. [7] proved that the terminal region will be enlarged by weighting the terminal cost for MPC without terminal constraint in its on-line optimization. Here, we will show that, this property will also hold water in our case.Theorem 4.1.
ConsiderF and Xf,max satisfying conditions (C1) and (C2). When a weighted terminal cost λF,λ≥1 is used, the corresponding terminal region Xλ,max is larger than Xf,max, namely,  Xf,max⊆Xλ,max.Proof.
For anyx0∈Xf,max, the conditions of (C1) and (C2) are equivalent to the fact that there exists a control trajectory u={u0,u1,u2,…} to make the following inequalities stand up:
(4.1)F(xi)≥q(xi,ui)+F(xi+1),i=0,1,2,…,
where xi+1=f(xi,ui).
It is obvious that, when the terminal cost isλF, these inequalities also stand up using the same control trajectory
(4.2)λF(xi)≥q(xi,ui)+λF(xi+1).
So, we can see that x0∈Xλ,max, namely, Xf,max⊆Xλ,max.Remark 4.2.
When the prediction horizon is given, the larger the terminal region is, the larger the domain of attraction is. So Theorem4.1 shows that, giving the prediction horizon, by weighting terminal cost, the domain of attraction will be enlarged.


### 4.2. Weighting Terminal Cost, the Total Cost Will Increase

It is known from the last section that, when the terminal cost is weighted, the terminal region is enlarged. So, can we weight the terminal cost arbitrarily? The answer is no. This section will tell us the reason the total cost will be increased by weighting the terminal cost.LetJN,λ*(x0) be the solution of PN,λ(x0), uN,λ*(x0)={uN,λ*(0,x0),…,uN,λ*(N,x0)}, xN,λ*(x0)={xN,λ*(0,x0),…,xN,λ*(N,x0)} be the optimal control trajectory and corresponding state trajectory, respectively, and let uRH,λ(x0)={uN*(0,x0),uN*(0,x1,λ),…},  xRH,λ(x0)={x0,x1,λ,x2,λ,…} be the receding horizon control trajectory and receding horizon state trajectory of using MPC with PN,λ(x0), respectively. Define JRH,λ(x0),  JRH(x0) as the total costs of using MPC with PN,λ(x0) and with PN(x0), respectively.For convenience, consider an assumption.Assumption.
For anyx0∈ΓN, where ΓN is the domain of attraction of MPC with PN(x0), the terminal state by solving PN,λ(x0) belongs to Xf,max, that is to say xN,λ*(N,x0)∈Xf,max.Remark 4.3.
Assumption1 means, for any x0∈ΓN, that the solution of PN,λ(x0) with Xλ,max as its terminal region is equal to that with Xf,max as its terminal region. A few points in ΓN may not satisfy this assumption, for convenience, their influence is neglected. Under this assumption, it is obviouse that JN,λ*(x0)≥JN*(x0) and the following lemma holds water.Lemma 4.4.
For anyx0∈ΓN, there exists
(4.3)(JN-1,λ*-JN-1*)(xN*(1,x0))≥(JN-1,λ*-JN-1*)(xN,λ*(1,x0)).Proof.
From the view of optimality,JN*(x0) can be expressed as
(4.4)JN*(x0)=q(x0,uN*(0,x0))+JN-1*(xN*(1,x0)).
Considering Assumption1 and by optimality, there exists
(4.5)JN*(x0)≤q(x0,uN,λ*(0,x0))+JN-1*(xN,λ*(1,x0)).
Similarly,JN,λ*(x0) can be expressed as
(4.6)JN,λ*(x0)=q(x0,uN,λ*(0,x0))+JN-1,λ*(xN,λ*(1,x0)).
And by optimality, there exists(4.7)JN,λ*(x0)≤q(x0,uN*(0,x0))+JN-1,λ*(xN*(1,x0)).
Obviously, the result of subtractingJN*(x0) from the right hand of (4.7) is bigger than the result of subtracting the right hand of (4.5) from JN,λ*(x0), in other words,
(4.8)(q(x0,uN*(0,x0))+JN-1,λ*(xN*(1,x0)))-JN*(x0)≥JN,λ*(x0)-(q(x0,uN,λ*(0,x0))+JN-1*(xN,λ*(1,x0))).
Finally, the following result can be obtained:(4.9)(JN-1,λ*-JN-1*)(xN*(1,x0))≥(JN-1,λ*-JN-1*)(xN,λ*(1,x0)).DefineF𝒦 as a kind of 𝒦 functions satisfying the following: for any F1,F2∈F𝒦, and x1,x2∈X, there exists the following:  if F1(x1)≥F1(x2), the inequality F2(x1)≥F2(x2) stands up.To continue discussion, another assumption is needed.Assumption.
All of the positive cost functions used in this paper likeJN-1*, JN*, JN-1,λ*, and JN,λ* and the results of the addition or subtraction between them like JN-1,λ*-JN-1*, JN-1,λ*+JN*  belong to F𝒦.Based on Assumption2 and Lemma 4.4, it is known that, for any x0∈ΓN and F𝒦, there exists(4.10)F𝒦(xN*(1,x0))≥F𝒦(xN,λ*(1,x0)).Then, by using (4.10) and Assumptions 1 and 2, another lemma which is a key for our study on this issue can be gotten.Lemma 4.5.
Under Assumptions1 and 2, for any x0∈ΓN and any positive cost function Ftra∈F𝒦 satisfying Ftra≤JN-1*, there exists
(4.11)q(x0,uN*(0,x0))+Ftra(x1)≤q(x0,uN,λ*(0,x0))+Ftra(x1,λ).Proof.
Here,x1 means xN*(1,x0) and x1,λ means xN,λ*(1,x0). From Assumption 2, it is known that (JN-1*-Ftra)∈F𝒦, so there exists
(4.12)q(x0,uN*(0,x0))+Ftra(x1)=q(x0,uN*(0,x0))+JN-1*(x1)-(JN-1*-Ftra)(x1)≤q(x0,uN,λ*(0,x0))+JN-1*(x1,λ)-(JN-1*-Ftra)(x1,λ)=q(x0,uN,λ*(0,x0))+Ftra(x1,λ).
Here, we used the fact thatJN*(x0)≤JN,λ*(x0) and  (JN-1*-Ftra)(x1)≥(JN-1*-Ftra)(x1,λ).Then, the reason why we cannot weight the terminal cost arbitrarily can be presented.Theorem 4.6.
Under Assumptions1 and 2, for any x0∈ΓN, the following inequality stands up:
(4.13)JRH(x0)≤JRH,λ(x0).Proof.
It is obvious thatJN* meets the condition in Lemma 4.5 because of JN*≤JN-1*. Choose Ftra=JN* in Lemma 4.5. There exists
(4.14)q(x0,uN*(0,x0))+JN*(x1)≤q(x0,uN,λ*(0,x0))+JN*(x1,λ).
Similarly, for anyx∈ΓN, the following result can be obtained:
(4.15)q(x,uN*(0,x))+JN*(xN*(1,x))≤q(x,uN*(0,x))+JN-1*(xN*(1,x))=JN*(x)≤JN-1*(x).
So, choose(4.16)Ftra(x)=q(x,uN*(0,x))+JN*(xN*(1,x)).
By usingx1,  x1,λ to replace x, respectively, there exists
(4.17)q(x0,uN*(0,x0))+q(x1,uN*(0,x1))+JN*(x2)=q(x0,uN*(0,x0))+Ftra(x1)≤q(x0,uN,λ*(0,x0))+Ftra(x1,λ)=q(x0,uN,λ*(0,x0))+q(x1,λ,uN*(0,x1,λ))+JN*(xN*(1,x1,λ)).
From Assumption1, it is known that, for x0∈ΓN, there exists x1,λ∈ΓN. Replacing x0 with x1,λ in (4.14), the following inequality can be gotten:
(4.18)q(x1,λ,uN*(0,x1,λ))+JN*(xN*(1,x1,λ))≤q(x1,λ,uN,λ*(0,x1,λ))+JN*(x2,λ).
So, there exists(4.19)q(x0,uN*(0,x0))+q(x1,uN*(0,x1))+JN*(x2)≤q(x0,uN,λ*(0,x0))+q(x1,λ,uN,λ*(0,x1,λ))+JN*(x2,λ)≤q(x0,uN,λ*(0,x0))+q(x1,λ,uN,λ*(0,x1,λ))+JN,λ*((x2,λ)).
Repeating this procedure, there exists(4.20)q(x0,uN*(0,x0))+⋯+q(xj,uN*(0,xj))+JN*(xj+1)≤q(x0,uN,λ*(0,x0))+⋯+q(xj,λ,uN,λ*(0,xj,λ))+JN,λ*((xj+1,λ)).
Letj→+∞, the final result can be obtained as follows:
(4.21)JRH(x0)≤JRH,λ(x0).Theorem4.6 shows that, weighting the terminal cost, the total cost will be increased. So, when choosing a terminal cost, people should not only take into account the need of enlarging the terminal region.


### 4.3. Getting an Appropriate Terminal Cost

It was pointed out from Theorems4.1 and 4.6 that the terminal cost is a double-edged sword. On its choice, two factors, the terminal region and the total cost, must be considered. With different emphasis, different terminal cost should be chosen.Here, a simple method to get a terminal cost is presented, whose basic idea is getting an initial terminal cost in advance then adjusting it according to our need.As mentioned, a good terminal cost should approximate toJ∞* as close as possible. People can only achieve it in a small neighborhood around the origin by using SOC method, see [2] for continuous-time system.Consider the linearization of the system (2.1) at the origin
(4.22)xk+1=Axk+Buk,
with A=(∂f/∂x)(0,0) and B=(∂f/∂u)(0,0).Here, assume that (4.22) is stabilizable, then a terminal cost which serves as an initial candidate can be found through the following procedure.Step 1.
Solving the Riccati equation to get a preparatoryG0,
(4.23)G0=ATG0A-(ATG0B)(BTG0B+R)-1(BTG0A)+Q.Step 2.
Getting a locally stabilizing linear state feedback gainK,
(4.24)K=-(BTG0B+R)-1(BTG0A).Step 3.
ComputingGK by solving the following Riccati equation:
(4.25)(αAK)TGK(αAK)-GK=-QK,
where AK=A+BK, QK=Q+KTRK, and α∈[1,+∞) is an adjustable parameter satisfying α|λmax(AK)|<1.Then,F(x)=xTGKx can serve as an initial terminal cost. According to our need and the properties of terminal cost, the initial one can be adjusted to get an appropriate one, F(x)=xT(λGK)x, λ>0. For example, if a larger terminal region is wanted and the total cost is not cared, λ can be set to be a larger number; otherwise, if a lower total cost is demanded and the domain of attraction already covers the operating region of system, a small one can be used.


## 4.1. Weighting Terminal Cost, the Domain of Attraction Will Be Enlarged

From conditions (C1) and (C2), it is known that the terminal region is based on the choice of the terminal cost. We want to know what the relationship between them is and this relationship will give us what messages when we choose terminal cost. Theorem4.1 will give us the answer.DenoteXf,max,Xλ,max as the terminal regions of PN(x0) and PN,λ(x0) with a weighted terminal cost λF, λ≥1, respectively. Limon et al. [7] proved that the terminal region will be enlarged by weighting the terminal cost for MPC without terminal constraint in its on-line optimization. Here, we will show that, this property will also hold water in our case.Theorem 4.1.
ConsiderF and Xf,max satisfying conditions (C1) and (C2). When a weighted terminal cost λF,λ≥1 is used, the corresponding terminal region Xλ,max is larger than Xf,max, namely,  Xf,max⊆Xλ,max.Proof.
For anyx0∈Xf,max, the conditions of (C1) and (C2) are equivalent to the fact that there exists a control trajectory u={u0,u1,u2,…} to make the following inequalities stand up:
(4.1)F(xi)≥q(xi,ui)+F(xi+1),i=0,1,2,…,
where xi+1=f(xi,ui).
It is obvious that, when the terminal cost isλF, these inequalities also stand up using the same control trajectory
(4.2)λF(xi)≥q(xi,ui)+λF(xi+1).
So, we can see that x0∈Xλ,max, namely, Xf,max⊆Xλ,max.Remark 4.2.
When the prediction horizon is given, the larger the terminal region is, the larger the domain of attraction is. So Theorem4.1 shows that, giving the prediction horizon, by weighting terminal cost, the domain of attraction will be enlarged.


## 4.2. Weighting Terminal Cost, the Total Cost Will Increase

It is known from the last section that, when the terminal cost is weighted, the terminal region is enlarged. So, can we weight the terminal cost arbitrarily? The answer is no. This section will tell us the reason the total cost will be increased by weighting the terminal cost.LetJN,λ*(x0) be the solution of PN,λ(x0), uN,λ*(x0)={uN,λ*(0,x0),…,uN,λ*(N,x0)}, xN,λ*(x0)={xN,λ*(0,x0),…,xN,λ*(N,x0)} be the optimal control trajectory and corresponding state trajectory, respectively, and let uRH,λ(x0)={uN*(0,x0),uN*(0,x1,λ),…},  xRH,λ(x0)={x0,x1,λ,x2,λ,…} be the receding horizon control trajectory and receding horizon state trajectory of using MPC with PN,λ(x0), respectively. Define JRH,λ(x0),  JRH(x0) as the total costs of using MPC with PN,λ(x0) and with PN(x0), respectively.For convenience, consider an assumption.Assumption.
For anyx0∈ΓN, where ΓN is the domain of attraction of MPC with PN(x0), the terminal state by solving PN,λ(x0) belongs to Xf,max, that is to say xN,λ*(N,x0)∈Xf,max.Remark 4.3.
Assumption1 means, for any x0∈ΓN, that the solution of PN,λ(x0) with Xλ,max as its terminal region is equal to that with Xf,max as its terminal region. A few points in ΓN may not satisfy this assumption, for convenience, their influence is neglected. Under this assumption, it is obviouse that JN,λ*(x0)≥JN*(x0) and the following lemma holds water.Lemma 4.4.
For anyx0∈ΓN, there exists
(4.3)(JN-1,λ*-JN-1*)(xN*(1,x0))≥(JN-1,λ*-JN-1*)(xN,λ*(1,x0)).Proof.
From the view of optimality,JN*(x0) can be expressed as
(4.4)JN*(x0)=q(x0,uN*(0,x0))+JN-1*(xN*(1,x0)).
Considering Assumption1 and by optimality, there exists
(4.5)JN*(x0)≤q(x0,uN,λ*(0,x0))+JN-1*(xN,λ*(1,x0)).
Similarly,JN,λ*(x0) can be expressed as
(4.6)JN,λ*(x0)=q(x0,uN,λ*(0,x0))+JN-1,λ*(xN,λ*(1,x0)).
And by optimality, there exists(4.7)JN,λ*(x0)≤q(x0,uN*(0,x0))+JN-1,λ*(xN*(1,x0)).
Obviously, the result of subtractingJN*(x0) from the right hand of (4.7) is bigger than the result of subtracting the right hand of (4.5) from JN,λ*(x0), in other words,
(4.8)(q(x0,uN*(0,x0))+JN-1,λ*(xN*(1,x0)))-JN*(x0)≥JN,λ*(x0)-(q(x0,uN,λ*(0,x0))+JN-1*(xN,λ*(1,x0))).
Finally, the following result can be obtained:(4.9)(JN-1,λ*-JN-1*)(xN*(1,x0))≥(JN-1,λ*-JN-1*)(xN,λ*(1,x0)).DefineF𝒦 as a kind of 𝒦 functions satisfying the following: for any F1,F2∈F𝒦, and x1,x2∈X, there exists the following:  if F1(x1)≥F1(x2), the inequality F2(x1)≥F2(x2) stands up.To continue discussion, another assumption is needed.Assumption.
All of the positive cost functions used in this paper likeJN-1*, JN*, JN-1,λ*, and JN,λ* and the results of the addition or subtraction between them like JN-1,λ*-JN-1*, JN-1,λ*+JN*  belong to F𝒦.Based on Assumption2 and Lemma 4.4, it is known that, for any x0∈ΓN and F𝒦, there exists(4.10)F𝒦(xN*(1,x0))≥F𝒦(xN,λ*(1,x0)).Then, by using (4.10) and Assumptions 1 and 2, another lemma which is a key for our study on this issue can be gotten.Lemma 4.5.
Under Assumptions1 and 2, for any x0∈ΓN and any positive cost function Ftra∈F𝒦 satisfying Ftra≤JN-1*, there exists
(4.11)q(x0,uN*(0,x0))+Ftra(x1)≤q(x0,uN,λ*(0,x0))+Ftra(x1,λ).Proof.
Here,x1 means xN*(1,x0) and x1,λ means xN,λ*(1,x0). From Assumption 2, it is known that (JN-1*-Ftra)∈F𝒦, so there exists
(4.12)q(x0,uN*(0,x0))+Ftra(x1)=q(x0,uN*(0,x0))+JN-1*(x1)-(JN-1*-Ftra)(x1)≤q(x0,uN,λ*(0,x0))+JN-1*(x1,λ)-(JN-1*-Ftra)(x1,λ)=q(x0,uN,λ*(0,x0))+Ftra(x1,λ).
Here, we used the fact thatJN*(x0)≤JN,λ*(x0) and  (JN-1*-Ftra)(x1)≥(JN-1*-Ftra)(x1,λ).Then, the reason why we cannot weight the terminal cost arbitrarily can be presented.Theorem 4.6.
Under Assumptions1 and 2, for any x0∈ΓN, the following inequality stands up:
(4.13)JRH(x0)≤JRH,λ(x0).Proof.
It is obvious thatJN* meets the condition in Lemma 4.5 because of JN*≤JN-1*. Choose Ftra=JN* in Lemma 4.5. There exists
(4.14)q(x0,uN*(0,x0))+JN*(x1)≤q(x0,uN,λ*(0,x0))+JN*(x1,λ).
Similarly, for anyx∈ΓN, the following result can be obtained:
(4.15)q(x,uN*(0,x))+JN*(xN*(1,x))≤q(x,uN*(0,x))+JN-1*(xN*(1,x))=JN*(x)≤JN-1*(x).
So, choose(4.16)Ftra(x)=q(x,uN*(0,x))+JN*(xN*(1,x)).
By usingx1,  x1,λ to replace x, respectively, there exists
(4.17)q(x0,uN*(0,x0))+q(x1,uN*(0,x1))+JN*(x2)=q(x0,uN*(0,x0))+Ftra(x1)≤q(x0,uN,λ*(0,x0))+Ftra(x1,λ)=q(x0,uN,λ*(0,x0))+q(x1,λ,uN*(0,x1,λ))+JN*(xN*(1,x1,λ)).
From Assumption1, it is known that, for x0∈ΓN, there exists x1,λ∈ΓN. Replacing x0 with x1,λ in (4.14), the following inequality can be gotten:
(4.18)q(x1,λ,uN*(0,x1,λ))+JN*(xN*(1,x1,λ))≤q(x1,λ,uN,λ*(0,x1,λ))+JN*(x2,λ).
So, there exists(4.19)q(x0,uN*(0,x0))+q(x1,uN*(0,x1))+JN*(x2)≤q(x0,uN,λ*(0,x0))+q(x1,λ,uN,λ*(0,x1,λ))+JN*(x2,λ)≤q(x0,uN,λ*(0,x0))+q(x1,λ,uN,λ*(0,x1,λ))+JN,λ*((x2,λ)).
Repeating this procedure, there exists(4.20)q(x0,uN*(0,x0))+⋯+q(xj,uN*(0,xj))+JN*(xj+1)≤q(x0,uN,λ*(0,x0))+⋯+q(xj,λ,uN,λ*(0,xj,λ))+JN,λ*((xj+1,λ)).
Letj→+∞, the final result can be obtained as follows:
(4.21)JRH(x0)≤JRH,λ(x0).Theorem4.6 shows that, weighting the terminal cost, the total cost will be increased. So, when choosing a terminal cost, people should not only take into account the need of enlarging the terminal region.


## 4.3. Getting an Appropriate Terminal Cost

It was pointed out from Theorems4.1 and 4.6 that the terminal cost is a double-edged sword. On its choice, two factors, the terminal region and the total cost, must be considered. With different emphasis, different terminal cost should be chosen.Here, a simple method to get a terminal cost is presented, whose basic idea is getting an initial terminal cost in advance then adjusting it according to our need.As mentioned, a good terminal cost should approximate toJ∞* as close as possible. People can only achieve it in a small neighborhood around the origin by using SOC method, see [2] for continuous-time system.Consider the linearization of the system (2.1) at the origin
(4.22)xk+1=Axk+Buk,
with A=(∂f/∂x)(0,0) and B=(∂f/∂u)(0,0).Here, assume that (4.22) is stabilizable, then a terminal cost which serves as an initial candidate can be found through the following procedure.Step 1.
Solving the Riccati equation to get a preparatoryG0,
(4.23)G0=ATG0A-(ATG0B)(BTG0B+R)-1(BTG0A)+Q.Step 2.
Getting a locally stabilizing linear state feedback gainK,
(4.24)K=-(BTG0B+R)-1(BTG0A).Step 3.
ComputingGK by solving the following Riccati equation:
(4.25)(αAK)TGK(αAK)-GK=-QK,
where AK=A+BK, QK=Q+KTRK, and α∈[1,+∞) is an adjustable parameter satisfying α|λmax(AK)|<1.Then,F(x)=xTGKx can serve as an initial terminal cost. According to our need and the properties of terminal cost, the initial one can be adjusted to get an appropriate one, F(x)=xT(λGK)x, λ>0. For example, if a larger terminal region is wanted and the total cost is not cared, λ can be set to be a larger number; otherwise, if a lower total cost is demanded and the domain of attraction already covers the operating region of system, a small one can be used.


## 5. Simulation Experiment

The model is an approximate discrete-time realization from a continuous-time system used in [2] as follows:
(5.1)[x1(k+1)x2(k+1)]=[1TT1][x1(k)x2(k)]+[TμTμ]uk+[T(1-μ)00-4T(1-μ)][x1(k)x2(k)]uk,
where μ=0.5, T=0.1s, and the state constraint and control constraint are X={x∣∥x∥1≤4}, U={u∣|u|≤2}, respectively.The stage cost is chosen asq(x,u)=xTQx+uTRu with Q=0.5I and R=1. By using SOC method, the locally linear feedback gain is adopted as K=-[2.01072.0107]   and |λmax(AK)|=0.9000 is obtained. Then, choose α=1.11  and get the terminal cost as F(x)=xTGx with G=[1107.356857.231;857.2311107.356].To estimate eachXfj, 4000 training points are generated. Set ε=1, when j=15, there exists(5.2)∑i=1Nsup,14∥O15(xi)-O14(xi)∥≤εNsup,14,
where xi∈Xsup,14, Xsup,14 is the support vectors set at j=14, and Nsup,14 is the number of support vectors. Then, it is deemed that X̂f15 is equal to X̂f14 in principle and X̂f15 can be taken as the final estimation of Xf,max. Figure 2 shows the approximation process of Xf,max. The blue line is the hyperplane at j=1, the black dot line is that at j=15, and the red lines between them are those at j=2,3,…,14. Let the prediction horizon be N=3. Figure 3 shows the closed-loop trajectories of some points chosen from the domain of attraction arbitrarily.Figure 2
The approximation process of terminal region.Figure 3
The closed-loop trajectories of points inΓ3.When the terminal cost is enlarged toF(x)=xT(10G)x, a new terminal region larger than the old one can be obtained. Figure 4 shows it. The red line is the new hyperplane and the black dot line is the old one.Figure 4
Comparison of the terminal regions.For convenience, let (A) denote the MPC usingxTGx as its terminal cost and (B) the MPC using xT(10G)x. For some points chosen from Γ3 of (A) arbitrarily, Figure 5 shows their closed-loop trajectories of using (A) and (B), respectively, where red lines denote the results of using (A), and blue dash-dotted lines denote the results of using (B). Table 1 shows the comparison of the total costs. Obviously, for the same point, the total cost of using (A) is smaller than that of using (B).Table 1
Comparison of the total costs.
Initial points[2-1][1-2][-12][-11](A)34.791518.246016.45287.9936(B)37.012020.671423.285423.8467Figure 5
Comparison of the closed-loop trajectories.


## 6. Conclusion

This paper discussed the relationships between terminal cost and terminal region and between terminal cost and total cost, respectively. It showed that, by enlarging the terminal cost, terminal region will be enlarged, but the total cost will be increased too. A simple method to get a suitable terminal cost was proposed, it can be adjusted according to our need. For example, to get a larger terminal region, it can be weighted; to reduce the total cost, it can be unweighted. When a terminal cost was given, a novel method was proposed to receive a maximal terminal region by using SVM. With the same prediction horizon, its corresponding domain of attraction is the largest one.



---

*Source: 853679-2010-10-11.xml*