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

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.

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: 𝑥_𝑘+1=𝑓(𝑥_𝑘,𝑢_𝑘), (2.1) where 𝑥_𝑘∈𝑅^𝑛,𝑢_𝑘∈𝑅^𝑚 are the state and the input of the system at sampling time 𝑡=𝑘, respectively. 𝑥_𝑘+1∈𝑅^𝑛 is the successor state and the mapping 𝑓∶𝑅^𝑛+𝑚↦𝑅^𝑛 with 𝑓(𝟎,𝟎)=𝟎 is known. The system is subject to constraints on both state and control action, and they are given by 𝑥_𝑘∈𝑋,𝑢_𝑘∈𝑈, where 𝑋 is a closed set and 𝑈 a compact set, both of them containing the origin. The control objective is usually to steer the state to the origin.

The optimization problem 𝑃_∞(𝑥₀) of SOC at the initial state 𝑥₀ can be stated as follows: min_{𝑢(𝑖,𝑥0)∈𝑈}𝐽_∞(𝐮,𝑥₀)=^∞_𝑖=0𝑞(𝑥(𝑖,𝑥₀),𝑢(𝑖,𝑥₀)),s.t.𝑥(𝑖+1,𝑥₀)=𝑓(𝑥(𝑖,𝑥₀),𝑢(𝑖,𝑥₀)),𝑥(𝑖+1,𝑥₀)∈𝑋,𝑢(𝑖,𝑥₀)∈𝑈, (2.2) where 𝑥(0,𝑥₀)=𝑥₀, 𝑞(𝑥,𝑢) is the stage cost and its form is chosen as 𝑞(𝑥,𝑢)=𝑥^𝑇𝑄𝑥+𝑢^𝑇𝑅𝑢 in which 𝑄,𝑅 are positive definites.

It is well known that the stability is guaranteed if 𝑃_∞(𝑥₀) has feasible solution. Here, assume that the solution of 𝑃_∞(𝑥₀) is known if it has the following. 𝐽^∗_∞(𝑥₀) is defined as the solution, 𝑢^∗_∞(𝑖,𝑥₀), 𝑖=1,2,…,∞, the optimal control trajectory and 𝑥^∗_∞(𝑖,𝑥₀), the corresponding state trajectory. As is well known, 𝐽^∗_∞(𝑥₀) is the total cost of using 𝑢^∗_∞(𝑖,𝑥₀) to drive 𝑥₀ to 0.

But solving 𝑃_∞(𝑥₀) is not an easy job, especially when 𝑓 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 problem 𝑃_𝑁(𝑥₀) of MPC at the initial state 𝑥₀ can be stated as min_{𝑢(𝑖,𝑥0)∈𝑈}𝐽_𝑁(𝐮,𝑥₀)=^𝑁−1_𝑖=0𝑞(𝑥(𝑖,𝑥₀),𝑢(𝑖,𝑥₀))+𝐹(𝑥(𝑁,𝑥₀)),s.t.𝑥(𝑖+1,𝑥₀)=𝑓(𝑥(𝑖,𝑥₀),𝑢(𝑖,𝑥₀)),𝑥(𝑖+1,𝑥₀)∈𝑋,𝑢(𝑖,𝑥₀)∈𝑈,𝑥(𝑁,𝑥₀)∈𝑋_𝑓, (2.3) where 𝑁 is prediction horizon, 𝑋_𝑓 is terminal region, and 𝐹 is terminal cost satisfying 𝐹(𝟎)=0 and 𝐹(𝑥)≥𝛼(‖𝑥‖), (the mapping 𝛼∶𝑅₊↦𝑅₊ satisfying 𝛼(0)=0 is continuous and strictly increasing, where 𝑅₊={𝑥∈𝑅∣𝑥>0}).

There exist many optimization algorithms to compute the solution of 𝑃_𝑁(𝑥₀). Let 𝐽^∗_𝑁(𝑥₀) be the solution, 𝐮^∗_𝑁(𝑥₀)={𝑢^∗_𝑁(0,𝑥₀),…,𝑢^∗_𝑁(𝑁−1,𝑥₀)},𝐱^∗_𝑁(𝑥₀)={𝑥^∗_𝑁(0,𝑥₀),…,𝑥^∗_𝑁(𝑁−1,𝑥₀)} be the optimal control trajectory and corresponding predicted state trajectory of 𝑃_𝑁(𝑥₀), respectively.

At sampling time 𝑡=0, 𝑢^∗_𝑁(0,𝑥₀) is inputted to the system. At the next sampling time 𝑡=1, 𝑥₁=𝑓(𝑥₀,𝑢^∗_𝑁(0,𝑥₀)) is outputted, and the control input 𝑢^∗_𝑁(0,𝑥₁) can be computed by solving the optimization problem 𝑃_𝑁(𝑥₁). By repeating this procedure, two trajectories can be obtained: 𝐮_RH(𝑥₀)={𝑢^∗_𝑁(0,𝑥₀),𝑢^∗_𝑁(0,𝑥₁),…},𝐱_RH(𝑥₀)={𝑥₀,𝑥₁,…}. Here, for convenience, they are called as the receding horizon control trajectory and the receding horizon state trajectory of MPC with 𝑃_𝑁(𝑥₀), respectively.

The introduction of 𝑋_𝑓 and 𝐹 in (2.3) is to guarantee the closed loop stability of system.

Lemma 2.1. Define Γ_𝑁∶={𝑥₀∈𝑋∣𝑥^∗_𝑁(𝑁,𝑥₀)∈𝑋_𝑓}. For any 𝑥₀∈Γ_𝑁 if 𝑋_𝑓 and 𝐹 satisfy two conditions as follows:
(C1)   𝐹 being a Lyapunov function, more strictly, for any 𝑥∈𝑋_𝑓, there exists 𝐹(𝑥)≥min_𝑢∈𝑈{𝑞(𝑥,𝑢)+𝐹(𝑓(𝑥,𝑢))}, (2.4)  
(C2)   𝑋_𝑓 being an invariant set. In other words, one has 𝑓(𝑥,𝑢)∈𝑋_𝑓. Here, 𝑢 is the control in (C1).

It is guaranteed that, 𝑥 will be led to 0 by using 𝐮_RH(𝑥₀).  Γ_𝑁 is called the domain of attraction.

The proof can be found in [1].

Obviously, the optimal choice of 𝐹(𝑥(𝑁,𝑥₀)) is 𝐽_∞(𝑥(𝑁,𝑥₀)). The total cost with this choice is 𝐽_∞(𝑥₀), namely it is the least one. But, as mentioned, 𝐽_∞ can not be obtained generally, so other ways should be found to get 𝐹. Certainly, the closer it approaches to 𝐽_∞, the better it is. When 𝐹 is given, 𝑋_𝑓 can be found to satisfy the conditions (C1) and (C2). There are many 𝑋_𝑓 satisfying the conditions and different 𝑋_𝑓 will be obtained by different methods. But, to make the domain of attraction the largest, it is wished that the largest 𝑋_𝑓 can be obtained. Here, define 𝑋_𝑓,max as the largest terminal region when 𝐹 is given, and in the next section, a novel method to get 𝑋_𝑓,max will be proposed.

3. Maximizing the Terminal Region for MPC

Until now, there exist many methods to construct 𝑋_𝑓. 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 𝑋_𝑓 computed by using this method does not approximate to 𝑋_𝑓,max to the largest extent.

In this paper, a novel method is proposed in which 𝑋_𝑓 is constructed directly from conditions (C1) and (C2).

3.1. Approximating the Largest Terminal Region Asymptotically

Define 𝑋_𝑓 as 𝑋_𝑓∶={𝑥∈𝑋∣𝐹(𝑥)≥𝐹^∗_𝑋𝑓(𝑥)}, (3.1) where 𝐹^∗_𝑋𝑓(𝑥)=min𝐹_𝑋𝑓(𝑥) is the solution of the following optimization problem: min_𝑢∈𝑈𝐹_𝑋𝑓(𝑥)=𝑞(𝑥,𝑢)+𝐹(𝑓(𝑥,𝑢)),s.t.𝑓(𝑥,𝑢)∈𝑋_𝑓. (3.2)

Obviously, for an 𝑥∈𝑋, it cannot be decided whether 𝑥 belongs to 𝑋_𝑓 from (3.1) and (3.2) when 𝑋_𝑓 is unknown. The difficulty is that the state constraint in the optimization (3.2) uses the 𝑋_𝑓 itself. To avoid it, the method of asymptotic approximation is adopted. Firstly, an initial set 𝑋⁰_𝑓 which can be obtained by the following discriminant is given: 𝑋⁰_𝑓∶={𝑥∈𝑋∣𝐹(𝑥)≥𝐹^∗_𝑋0𝑓(𝑥)}, (3.3) where 𝐹^∗_𝑋0𝑓(𝑥)  is the solution of min_𝑢∈𝑈𝐹_𝑋⁰_𝑓(𝑥)=𝑞(𝑥,𝑢)+𝐹(𝑓(𝑥,𝑢)),s.t.𝑓(𝑥,𝑢)∈𝑋. (3.4)

Then, using 𝑋⁰_𝑓 instead of 𝑋 in the state constraint in optimization of (3.4), 𝑋¹_𝑓 can be obtained. One by one, 𝑋²_𝑓，𝑋³_𝑓,… will be obtained. The whole procedure can be pictured as in Figure 1.

Figure 1 

The approximation of 𝑋𝑓,max.

In Figure 1, 𝑋¹_𝑓 is defined as 𝑋¹_𝑓∶={𝑥∈𝑋⁰_𝑓∣𝐹(𝑥)≥𝐹^∗_𝑋1𝑓(𝑥)}, (3.5) where 𝐹^∗_𝑋1𝑓(𝑥) is the solution of min_𝑢∈𝑈𝐹_𝑋¹_𝑓(𝑥)=𝑞(𝑥,𝑢)+𝐹(𝑓(𝑥,𝑢)),s.t.𝑓(𝑥,𝑢)∈𝑋⁰_𝑓. (3.6)

Similarly, 𝑋^𝑗_𝑓 can be defined as：𝑋^𝑗_𝑓∶={𝑥∈𝑋_𝑓^𝑗−1∣𝐹(𝑥)≥𝐹^∗_𝑋𝑗𝑓(𝑥)}, (3.7) and 𝐹^∗_𝑋𝑗𝑓(𝑥) is the solution of min_𝑢∈𝑈𝐹_𝑋^𝑗_𝑓(𝑥)=𝑞(𝑥,𝑢)+𝐹(𝑓(𝑥,𝑢)),s.t.𝑓(𝑥,𝑢)∈𝑋_𝑓^𝑗−1. (3.8)

Obviously, as Figure 1 shows, there exists 𝑋⁰_𝑓⊇𝑋¹_𝑓⊇⋯⊇𝑋^𝑗_𝑓⊇𝑋_𝑓. As 𝑗 increases, 𝑋^𝑗_𝑓 will converge to a set denoted by 𝑋_𝑓^+∞. But, whether is 𝑋_𝑓^+∞ the terminal region we want? Theorem 3.1 provides the answer to this question.

Theorem 3.1. For 𝑋^𝑗_𝑓 defined in (3.7), when 𝑗 goes to infinity, 𝑋^𝑗_𝑓 will converge to 𝑋_𝑓,max, namely, as 𝑗→+∞, one has 𝑋^𝑗_𝑓→𝑋_𝑓,max.

This theorem is proved by contradiction.

Proof. (A) Assume that there exists 𝑋_spo which satisfies 𝑋_spo⊃𝑋_𝑓,max and 𝑋^𝑗_𝑓→𝑋_spo, as 𝑗→+∞ then for any 𝑥∈𝑋_spo, we have 𝐹(𝑥)≥min_𝑢∈𝑈{𝑞(𝑥,𝑢)+𝐹(𝑓(𝑥,𝑢))} and 𝑓(𝑥,𝑢)∈𝑋_spo. Obviously this is contradicted with 𝑋_𝑓,max being the largest one satisfying (C1) and (C2).
(B) Similarly, assume that there exists 𝑋_spo which satisfies 𝑋_spo⊂𝑋_𝑓,max and 𝑋^𝑗_𝑓→𝑋_spo, as 𝑗→+∞, then there exists 𝑁 with 0≤𝑁<+∞ satisfying 𝑋^𝑁_𝑓⊇𝑋_𝑓,max and 𝑋_𝑓,max⧵𝑋_𝑓^𝑁+1≠Φ, where Φ denotes the empty set. Choose any 𝑥∈𝑋_𝑓,max⧵𝑋_𝑓^𝑁+1, it is obviouse that, 𝑥 satisfies 𝐹(𝑥)≥min_𝑢∈𝑈{𝑞(𝑥,𝑢)+𝐹(𝑓(𝑥,𝑢))} and 𝑓(𝑥,𝑢)∈𝑋_𝑓,max⊆𝑋^𝑁_𝑓. On the other hand, we know that 𝑥∈𝑋^𝑁_𝑓, so 𝑥 meets the conditions in the definition of 𝑋_𝑓^𝑁+1 and we have 𝑥∈𝑋_𝑓^𝑁+1. This is contradicted with 𝑥∈𝑋_𝑓,max⧵𝑋_𝑓^𝑁+1.

Remark 3.2. Generally, in the computation of 𝑋_𝑓,max, it is impossible to keep computation until 𝑗→+∞. So, when the iteration goes to 𝑗=𝑁, if 𝑋^𝑁_𝑓 is equal to 𝑋_𝑓^𝑁−1 in principle, 𝑋^𝑁_𝑓 can be taken as the terminal region we want.

Remark 3.3. The terminal region computed through the method in Remark 3.2 is not 𝑋_𝑓,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 𝑋_𝑓,max, denoted by 𝑋_𝑓,0 in advance and using 𝑋_𝑓,0 to serve as the state constraint in (3.4), then a larger region 𝑋_𝑓,1 will be computed. By the same procedure as presented in Section 3.1, a terminal region which is a subset of 𝑋_𝑓,max will be gotten.

To obtain 𝑋^𝑗_𝑓 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 separating 𝑋 into 𝐴 and 𝑋⧵𝐴. For each 𝑥_𝑖∈𝐴, an additional variable 𝑦_𝑖=+1 is introduced. Similarly, for each 𝑥_𝑖∈𝑋⧵𝐴, 𝑦_𝑖=−1 is introduced. Define 𝐼⁺∶={𝑖∶𝑦_𝑖=+1}, 𝐼⁻∶={𝑖∶𝑦_𝑖=−1}. SVM is used to find a separating hyperplane 𝑂(𝑥)∶=𝑤⋅𝜙(𝑥_𝑖)+𝑏=0, between 𝐴 and 𝑋⧵𝐴. Then, we obtain an estimated set of 𝐴, 𝐴={𝑥∈𝑋∣𝑂(𝑥)≥0}.  𝑂(𝑥) can be obtained by solving the following problem: min_𝛼12_𝑖𝑗𝛼_𝑖𝛼_𝑗𝑦_𝑖𝑦_𝑗ker(𝑥_𝑖,𝑥_𝑗)−_𝑖𝛼_𝑖,s.t._𝑖𝛼_𝑖𝑦_𝑖=0,0≤𝛼_𝑖≤𝐶,∀𝑖∈𝐼⁺,𝛼_𝑖≥0,∀𝑖∈𝐼⁻, (3.9) where ker denotes the kernel function.

In this paper, the following Gaussian kernel is used: ker(𝑥,𝑥_𝑖)=exp(−‖𝑥−𝑥_𝑖‖²2𝜎²), (3.10) 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 {𝑥_𝑖} and their relevant weights are exported. Denote 𝑃_𝑠 as the number of support vectors and 𝑋_𝑠 as the support vectors set, the optimal hyperplane is described as follows: 𝑂(𝑥)=^𝑃𝑠_𝑖=1𝑤_𝑖⋅ker(𝑥_𝑖,𝑥)+𝑏, (3.11) where 𝑥_𝑖∈𝑋_𝑠 is a support vector, and 𝑤_𝑖=𝛼_𝑖𝑦_𝑖 satifying ∑^𝑃𝑠_𝑖=1𝑤_𝑖=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 of 𝑋⁰_𝑓 from 𝑋 as an example. Firstly, choose arbitrary points 𝑥_𝑖∈𝑋, 𝑖=1,2,…,𝑁_𝑝 (𝑁_𝑝 is the number of training points) then decide the value of 𝑦_𝑖 corresponding to 𝑥_𝑖 by using the following procedure

If 𝐹(𝑥_𝑖)≥𝐹^∗_𝑋0𝑓(𝑥_𝑖), 

𝑦_𝑖=+1; else 𝑦_𝑖=−1

endif.

When all the 𝑦_𝑖 for all the 𝑥_𝑖 are gotten, they can be packed to constitute the training data. Then, by inputting the training data into SVM classifier, an optimal hyperplane 𝑂⁰(𝑥)=0 and an estimated set of 𝑋⁰_𝑓, 𝑋⁰_𝑓={𝑥∈𝑋∣𝑂⁰(𝑥)≥0} will be obtained.

When 𝑋⁰_𝑓 is known, the training data for separating 𝑋¹_𝑓 from 𝑋⁰_𝑓 can be known by the similar procedure. By inputting them into SVM classifier, a hyperplane 𝑂¹(𝑥)=0 and an estimated set of 𝑋¹_𝑓, 𝑋¹_𝑓={𝑥∈𝑋∣𝑂¹(𝑥)≥0} will be gotten.

Repeatedly, a hyperplane series 𝑂^𝑗(𝑥)=0,𝑗=1,2,… will be obtained. When 𝑗=𝑁 if it is satisfied that, for 𝑥_𝑖∈𝑋_sup,𝑁−1, 𝑖=1,2,…,𝑁_sup,𝑁−1, there exists ^{𝑁sup,𝑁−1}_𝑖=1‖𝑂^𝑁(𝑥_𝑖)−𝑂^𝑁−1(𝑥_𝑖)‖≤𝜀𝑁_sup,𝑁−1, (3.12) it is deemed that 𝑋^𝑁_𝑓 is equal to 𝑋_𝑓^𝑁−1 in principle and 𝑋^𝑁_𝑓 is taken as the final estimation of 𝑋_𝑓,max. Where 𝑋_sup,𝑁−1 is the support vectors set at 𝑗=𝑁−1, 𝑁_sup,𝑁−1 is the number of support vectors and 𝜀 is a tunable threshold. The smaller it is, the higher the precision of 𝑋^𝑁_𝑓 approximating to 𝑋_𝑓,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. Theorem 4.1 will give us the answer.

Denote 𝑋_𝑓,max,𝑋_𝜆,max as the terminal regions of 𝑃_𝑁(𝑥₀) and 𝑃_𝑁,𝜆(𝑥₀) with a weighted terminal cost 𝜆𝐹, 𝜆≥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. Consider 𝐹 and 𝑋_𝑓,max satisfying conditions (C1) and (C2). When a weighted terminal cost 𝜆𝐹,𝜆≥1 is used, the corresponding terminal region 𝑋_𝜆,max is larger than 𝑋_𝑓,max, namely, 𝑋_𝑓,max⊆𝑋_𝜆,max.

Proof. For any 𝑥₀∈𝑋_𝑓,max, the conditions of (C1) and (C2) are equivalent to the fact that there exists a control trajectory 𝐮={𝑢₀,𝑢₁,𝑢₂,…} to make the following inequalities stand up: 𝐹(𝑥_𝑖)≥𝑞(𝑥_𝑖,𝑢_𝑖)+𝐹(𝑥_𝑖+1),𝑖=0,1,2,…, (4.1) where 𝑥_𝑖+1=𝑓(𝑥_𝑖,𝑢_𝑖).
It is obvious that, when the terminal cost is 𝜆𝐹, these inequalities also stand up using the same control trajectory 𝜆𝐹(𝑥_𝑖)≥𝑞(𝑥_𝑖,𝑢_𝑖)+𝜆𝐹(𝑥_𝑖+1). (4.2) So, we can see that 𝑥₀∈𝑋_𝜆,max, namely, 𝑋_𝑓,max⊆𝑋_𝜆,max.

Remark 4.2. When the prediction horizon is given, the larger the terminal region is, the larger the domain of attraction is. So Theorem 4.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.

Let 𝐽^∗_𝑁,𝜆(𝑥₀) be the solution of 𝑃_𝑁,𝜆(𝑥₀), 𝐮^∗_𝑁,𝜆(𝑥₀)={𝑢^∗_𝑁,𝜆(0,𝑥₀),…,𝑢^∗_𝑁,𝜆(𝑁,𝑥₀)}, 𝐱^∗_𝑁,𝜆(𝑥₀)={𝑥^∗_𝑁,𝜆(0,𝑥₀),…,𝑥^∗_𝑁,𝜆(𝑁,𝑥₀)} be the optimal control trajectory and corresponding state trajectory, respectively, and let 𝐮_RH,𝜆(𝑥₀)={𝑢^∗_𝑁(0,𝑥₀),𝑢^∗_𝑁(0,𝑥_1,𝜆),…},  𝐱_RH,𝜆(𝑥₀)={𝑥₀,𝑥_1,𝜆,𝑥_2,𝜆,…} be the receding horizon control trajectory and receding horizon state trajectory of using MPC with 𝑃_𝑁,𝜆(𝑥₀), respectively. Define 𝐽_RH,𝜆(𝑥₀),  𝐽_RH(𝑥₀) as the total costs of using MPC with 𝑃_𝑁,𝜆(𝑥₀) and with 𝑃_𝑁(𝑥₀), respectively.

For convenience, consider an assumption.

Assumption. For any 𝑥₀∈Γ_𝑁, where Γ_𝑁 is the domain of attraction of MPC with 𝑃_𝑁(𝑥₀), the terminal state by solving 𝑃_𝑁,𝜆(𝑥₀) belongs to 𝑋_𝑓,max, that is to say 𝑥^∗_𝑁,𝜆(𝑁,𝑥₀)∈𝑋_𝑓,max.

Remark 4.3. Assumption 1 means, for any 𝑥₀∈Γ_𝑁, that the solution of 𝑃_𝑁,𝜆(𝑥₀) with 𝑋_𝜆,max as its terminal region is equal to that with 𝑋_𝑓,max as its terminal region. A few points in Γ_𝑁 may not satisfy this assumption, for convenience, their influence is neglected. Under this assumption, it is obviouse that 𝐽^∗_𝑁,𝜆(𝑥₀)≥𝐽^∗_𝑁(𝑥₀) and the following lemma holds water.

Lemma 4.4. For any 𝑥₀∈Γ_𝑁, there exists (𝐽^∗_{𝑁−1,𝜆}−𝐽^∗_𝑁−1)(𝑥^∗_N(1,𝑥₀))≥(𝐽^∗_{𝑁−1,𝜆}−𝐽^∗_𝑁−1)(𝑥^∗_N,𝜆(1,𝑥₀)). (4.3)

Proof. From the view of optimality, 𝐽^∗_𝑁(𝑥₀) can be expressed as 𝐽^∗_𝑁(𝑥₀)=𝑞(𝑥₀,𝑢^∗_𝑁(0,𝑥₀))+𝐽^∗_𝑁−1(𝑥^∗_𝑁(1,𝑥₀)). (4.4)
Considering Assumption 1 and by optimality, there exists 𝐽^∗_𝑁(𝑥₀)≤𝑞(𝑥₀,𝑢^∗_𝑁,𝜆(0,𝑥₀))+𝐽^∗_𝑁−1(𝑥^∗_𝑁,𝜆(1,𝑥₀)). (4.5)
Similarly, 𝐽^∗_𝑁,𝜆(𝑥₀) can be expressed as 𝐽^∗_𝑁,𝜆(𝑥₀)=𝑞(𝑥₀,𝑢^∗_𝑁,𝜆(0,𝑥₀))+𝐽^∗_{𝑁−1,𝜆}(𝑥^∗_𝑁,𝜆(1,𝑥₀)). (4.6)
And by optimality, there exists 𝐽^∗_𝑁,𝜆(𝑥₀)≤𝑞(𝑥₀,𝑢^∗_𝑁(0,𝑥₀))+𝐽^∗_{𝑁−1,𝜆}(𝑥^∗_𝑁(1,𝑥₀)). (4.7)
Obviously, the result of subtracting 𝐽^∗_𝑁(𝑥₀) from the right hand of (4.7) is bigger than the result of subtracting the right hand of (4.5) from 𝐽^∗_𝑁,𝜆(𝑥₀), in other words, (𝑞(𝑥₀,𝑢^∗_𝑁(0,𝑥₀))+𝐽^∗_{𝑁−1,𝜆}(𝑥^∗_𝑁(1,𝑥₀)))−𝐽^∗_𝑁(𝑥₀)≥𝐽^∗_𝑁,𝜆(𝑥₀)−(𝑞(𝑥₀,𝑢^∗_𝑁,𝜆(0,𝑥₀))+𝐽^∗_𝑁−1(𝑥^∗_𝑁,𝜆(1,𝑥₀))). (4.8)
Finally, the following result can be obtained: (𝐽^∗_{𝑁−1,𝜆}−𝐽^∗_𝑁−1)(𝑥^∗_𝑁(1,𝑥₀))≥(𝐽^∗_{𝑁−1,𝜆}−𝐽^∗_𝑁−1)(𝑥^∗_𝑁,𝜆(1,𝑥₀)). (4.9)

Define 𝐹_𝒦 as a kind of 𝒦 functions satisfying the following: for any 𝐹₁,𝐹₂∈𝐹_𝒦, and 𝑥₁,𝑥₂∈𝑋, there exists the following: if 𝐹₁(𝑥₁)≥𝐹₁(𝑥₂), the inequality 𝐹₂(𝑥₁)≥𝐹₂(𝑥₂) stands up.

To continue discussion, another assumption is needed.

Assumption. All of the positive cost functions used in this paper like 𝐽^∗_𝑁−1, 𝐽^∗_𝑁, 𝐽^∗_{𝑁−1,𝜆}, and 𝐽^∗_𝑁,𝜆 and the results of the addition or subtraction between them like 𝐽^∗_{𝑁−1,𝜆}−𝐽^∗_𝑁−1, 𝐽^∗_{𝑁−1,𝜆}+𝐽^∗_𝑁 belong to 𝐹_𝒦.

Based on Assumption 2 and Lemma 4.4, it is known that, for any 𝑥₀∈Γ_𝑁 and 𝐹_𝒦, there exists𝐹_𝒦(𝑥^∗_𝑁(1,𝑥₀))≥𝐹_𝒦(𝑥^∗_𝑁,𝜆(1,𝑥₀)). (4.10)

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 Assumptions 1 and 2, for any 𝑥₀∈Γ_𝑁 and any positive cost function 𝐹_tra∈𝐹_𝒦 satisfying 𝐹_tra≤𝐽^∗_𝑁−1, there exists 𝑞(𝑥₀,𝑢^∗_N(0,𝑥₀))+𝐹_tra(𝑥₁)≤𝑞(𝑥₀,𝑢^∗_N,𝜆(0,𝑥₀))+𝐹_tra(𝑥_1,𝜆). (4.11)

Proof. Here, 𝑥₁ means 𝑥^∗_𝑁(1,𝑥₀) and 𝑥_1,𝜆 means 𝑥^∗_𝑁,𝜆(1,𝑥₀). From Assumption 2, it is known that (𝐽^∗_𝑁−1−𝐹_tra)∈𝐹_𝒦, so there exists 𝑞(𝑥₀,𝑢^∗_𝑁(0,𝑥₀))+𝐹_tra(𝑥₁)=𝑞(𝑥₀,𝑢^∗_𝑁(0,𝑥₀))+𝐽^∗_𝑁−1(𝑥₁)−(𝐽^∗_𝑁−1−𝐹_tra)(𝑥₁)≤𝑞(𝑥₀,𝑢^∗_𝑁,𝜆(0,𝑥₀))+𝐽^∗_𝑁−1(𝑥_1,𝜆)−(𝐽^∗_𝑁−1−𝐹_tra)(𝑥_1,𝜆)=𝑞(𝑥₀,𝑢^∗_𝑁,𝜆(0,𝑥₀))+𝐹_tra(𝑥_1,𝜆). (4.12)
Here, we used the fact that 𝐽^∗_𝑁(𝑥₀)≤𝐽^∗_𝑁,𝜆(𝑥₀) and (𝐽^∗_𝑁−1−𝐹_tra)(𝑥₁)≥(𝐽^∗_𝑁−1−𝐹_tra)(𝑥_1,𝜆).

Then, the reason why we cannot weight the terminal cost arbitrarily can be presented.

Theorem 4.6. Under Assumptions 1 and 2, for any 𝑥₀∈Γ_𝑁, the following inequality stands up: 𝐽_RH(𝑥₀)≤𝐽_RH,𝜆(𝑥₀). (4.13)

Proof. It is obvious that 𝐽^∗_𝑁 meets the condition in Lemma 4.5 because of 𝐽^∗_𝑁≤𝐽^∗_𝑁−1. Choose 𝐹_tra=𝐽^∗_𝑁 in Lemma 4.5. There exists 𝑞(𝑥₀,𝑢^∗_𝑁(0,𝑥₀))+𝐽^∗_𝑁(𝑥₁)≤𝑞(𝑥₀,𝑢^∗_𝑁,𝜆(0,𝑥₀))+𝐽^∗_𝑁(𝑥_1,𝜆). (4.14)
Similarly, for any 𝑥∈Γ_𝑁, the following result can be obtained: 𝑞(𝑥,𝑢^∗_𝑁(0,𝑥))+𝐽^∗_𝑁(𝑥^∗_𝑁(1,𝑥))≤𝑞(𝑥,𝑢^∗_𝑁(0,𝑥))+𝐽^∗_𝑁−1(𝑥^∗_𝑁(1,𝑥))=𝐽^∗_𝑁(𝑥)≤𝐽^∗_𝑁−1(𝑥). (4.15)
So, choose 𝐹_tra(𝑥)=𝑞(𝑥,𝑢^∗_𝑁(0,𝑥))+𝐽^∗_𝑁(𝑥^∗_𝑁(1,𝑥)). (4.16)
By using 𝑥₁,  𝑥_1,𝜆 to replace 𝑥, respectively, there exists 𝑞(𝑥₀,𝑢^∗_𝑁(0,𝑥₀))+𝑞(𝑥₁,𝑢^∗_𝑁(0,𝑥₁))+𝐽^∗_𝑁(𝑥₂)=𝑞(𝑥₀,𝑢^∗_𝑁(0,𝑥₀))+𝐹_tra(𝑥₁)≤𝑞(𝑥₀,𝑢^∗_𝑁,𝜆(0,𝑥₀))+𝐹_tra(𝑥_1,𝜆)=𝑞(𝑥₀,𝑢^∗_𝑁,𝜆(0,𝑥₀))+𝑞(𝑥_1,𝜆,𝑢^∗_𝑁(0,𝑥_1,𝜆))+𝐽^∗_𝑁(𝑥^∗_𝑁(1,𝑥_1,𝜆)). (4.17)
From Assumption 1, it is known that, for 𝑥₀∈Γ_𝑁, there exists 𝑥_1,𝜆∈Γ_𝑁. Replacing 𝑥₀ with 𝑥_1,𝜆 in (4.14), the following inequality can be gotten: 𝑞(𝑥_1,𝜆,𝑢^∗_𝑁(0,𝑥_1,𝜆))+𝐽^∗_𝑁(𝑥^∗_𝑁(1,𝑥_1,𝜆))≤𝑞(𝑥_1,𝜆,𝑢^∗_𝑁,𝜆(0,𝑥_1,𝜆))+𝐽^∗_𝑁(𝑥_2,𝜆). (4.18)
So, there exists 𝑞(𝑥₀,𝑢^∗_𝑁(0,𝑥₀))+𝑞(𝑥₁,𝑢^∗_𝑁(0,𝑥₁))+𝐽^∗_𝑁(𝑥₂)≤𝑞(𝑥₀,𝑢^∗_𝑁,𝜆(0,𝑥₀))+𝑞(𝑥_1,𝜆,𝑢^∗_𝑁,𝜆(0,𝑥_1,𝜆))+𝐽^∗_𝑁(𝑥_2,𝜆)≤𝑞(𝑥₀,𝑢^∗_𝑁,𝜆(0,𝑥₀))+𝑞(𝑥_1,𝜆,𝑢^∗_𝑁,𝜆(0,𝑥_1,𝜆))+𝐽^∗_𝑁,𝜆((𝑥_2,𝜆)). (4.19)
Repeating this procedure, there exists 𝑞(𝑥₀,𝑢^∗_𝑁(0,𝑥₀))+⋯+𝑞(𝑥_𝑗,𝑢^∗_𝑁(0,𝑥_𝑗))+𝐽^∗_𝑁(𝑥_𝑗+1)≤𝑞(𝑥₀,𝑢^∗_𝑁,𝜆(0,𝑥₀))+⋯+𝑞(𝑥_𝑗,𝜆,𝑢^∗_𝑁,𝜆(0,𝑥_𝑗,𝜆))+𝐽^∗_𝑁,𝜆((𝑥_𝑗+1,𝜆)). (4.20)
Let 𝑗→+∞, the final result can be obtained as follows: 𝐽_RH(𝑥₀)≤𝐽_RH,𝜆(𝑥₀). (4.21)

Theorem 4.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 Theorems 4.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 to 𝐽^∗_∞ 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 𝑥_𝑘+1=𝐴𝑥_𝑘+𝐵𝑢_𝑘, (4.22) with 𝐴=(𝜕𝑓/𝜕𝑥)(0,0) and 𝐵=(𝜕𝑓/𝜕𝑢)(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 preparatory 𝐺₀, 𝐺₀=𝐴^𝑇𝐺₀𝐴−(𝐴^𝑇𝐺₀𝐵)(𝐵^𝑇𝐺₀𝐵+𝑅)⁻¹(𝐵^𝑇𝐺₀𝐴)+𝑄. (4.23)

Step 2. Getting a locally stabilizing linear state feedback gain 𝐾, 𝐾=−(𝐵^𝑇𝐺₀𝐵+𝑅)⁻¹(𝐵^𝑇𝐺₀𝐴). (4.24)

Step 3. Computing 𝐺_𝐾 by solving the following Riccati equation: (𝛼𝐴_𝐾)^𝑇𝐺_𝐾(𝛼𝐴_𝐾)−𝐺_𝐾=−𝑄_𝐾, (4.25) where 𝐴_𝐾=𝐴+𝐵𝐾, 𝑄_𝐾=𝑄+𝐾^𝑇𝑅𝐾, and 𝛼∈[1,+∞) is an adjustable parameter satisfying 𝛼|𝜆_max(𝐴_𝐾)|<1.

Then, 𝐹(𝑥)=𝑥^𝑇𝐺_𝐾𝑥 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, 𝐹(𝑥)=𝑥^𝑇(𝜆𝐺_𝐾)𝑥, 𝜆>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: 𝑥₁𝑥(𝑘+1)₂=𝑥(𝑘+1)1𝑇𝑇1𝑥¹(𝑘)₂+𝑢(𝑘)𝑇𝜇𝑇𝜇^𝑘+𝑇(1−𝜇)0𝑥0−4𝑇(1−𝜇)𝑥¹(𝑘)₂𝑢(𝑘)^𝑘, (5.1) where 𝜇=0.5, 𝑇=0.1s, and the state constraint and control constraint are 𝑋={𝑥∣‖𝑥‖₁≤4}, 𝑈={𝑢∣|𝑢|≤2}, respectively.

The stage cost is chosen as 𝑞(𝑥,𝑢)=𝑥^𝑇𝑄𝑥+𝑢^𝑇𝑅𝑢 with 𝑄=0.5𝐼 and 𝑅=1. By using SOC method, the locally linear feedback gain is adopted as 𝐾=−[2.01072.0107]   and |𝜆_max(𝐴_𝐾)|=0.9000 is obtained. Then, choose 𝛼=1.11  and get the terminal cost as 𝐹(𝑥)=𝑥^𝑇𝐺𝑥 with 𝐺=[1107.356857.231;857.2311107.356].

To estimate each 𝑋^𝑗_𝑓, 4000 training points are generated. Set 𝜀=1, when 𝑗=15, there exists^𝑁sup,14_𝑖=1‖𝑂¹⁵(𝑥_𝑖)−𝑂¹⁴(𝑥_𝑖)‖≤𝜀𝑁_sup,14, (5.2) where 𝑥_𝑖∈𝑋_sup,14, 𝑋_sup,14 is the support vectors set at 𝑗=14, and 𝑁_sup,14 is the number of support vectors. Then, it is deemed that 𝑋_𝑓¹⁵ is equal to 𝑋_𝑓¹⁴ in principle and 𝑋_𝑓¹⁵ can be taken as the final estimation of 𝑋_𝑓,max. Figure 2 shows the approximation process of 𝑋_𝑓,max. The blue line is the hyperplane at 𝑗=1, the black dot line is that at 𝑗=15, and the red lines between them are those at 𝑗=2,3,…,14. Let the prediction horizon be 𝑁=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 to 𝐹(𝑥)=𝑥^𝑇(10𝐺)𝑥, 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 using 𝑥^𝑇𝐺𝑥 as its terminal cost and (B) the MPC using 𝑥^𝑇(10𝐺)𝑥. For some points chosen from Γ₃ 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.

Figure 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.