An improved design of aggregation-based model predictive control
Systems & Control Letters 62 (2013) 1082–1089
Contents lists available at ScienceDirect
Systems & Control Letters journal homepage: www.elsevier.com/locate/sysconle
An improved design of aggregation-based model predictive control Dewei Li a,∗ , Yugeng Xi a , Zongli Lin b,a a
Department of Automation, Shanghai Jiao Tong University, Key Laboratory of System Control and Information Processing, Ministry of Education, Shanghai, 200240, China b
Charles L. Brown Department of Electrical and Computer Engineering, University of Virginia, Charlottesville, VA 22904-4743, USA
article
info
Article history: Received 28 January 2013 Received in revised form 12 August 2013 Accepted 27 August 2013 Available online 30 September 2013 Keywords: Model predictive control Aggregation Equivalent aggregation Quasi-equivalent aggregation Blocking
abstract Linear aggregation in the input is an effective method to reduce the online computational burden of model predictive control (MPC) but at the cost of degradations in the closed-loop performance. In this paper, an improved aggregation-based MPC algorithm is developed to reduce these degradations. In this algorithm, a time-varying base vector is utilized in conjunction with the quasi-equivalent aggregation strategy. Furthermore, by relaxing the constraints with a sequence of reachable sets, a switching strategy is adopted to enlarge the attractive region of the resulting aggregation-based MPC. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Model predictive control (MPC), also referred to as receding horizon control, emerged in 1970s. Because of its ability to tackle the constraints explicitly and to achieve good control performance, MPC has been attracting more and more attention from both the industrial and academic communities (see, e.g., [1,2]). MPC is characterized by rolling horizon optimization. An optimal control problem over a finite horizon is solved online at each sampling time and the first optimal control input is executed on the controlled plant. This feature leads to several advantages of MPC over other control algorithms but at the cost of heavy online computational burden. Increasing the number of optimization variables, i.e., the control inputs to be optimized, may further improve the control performance due to the increased degree of freedom, but the demand on online computation will also further intensify. As a result, the online computational burden has always been an important issue in the practical application of MPC. Input linear aggregation [3], also known as input parametrization [2], has proven to be an effective strategy to reduce the online computational burden of an MPC algorithm. The idea can be traced back to the blocking technology in [4], where the control horizon is divided into several blocks and the control inputs within each block
∗
Corresponding author. Tel.: +86 21 34204550. E-mail addresses: [email protected] (D. Li), [email protected] (Y. Xi), [email protected] (Z. Lin). 0167-6911/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.sysconle.2013.08.007
are assumed to be the same to reduce the number of online optimization variables. Similar strategies include the Predictive Functional Control (PFC) in [5]. However, the aggregation strategy also leads to performance degradations, including poor control performance, loss of guaranteed closed-loop stability and reduction in size of the attractive region. In order to overcome the difficulty caused by the structural restriction imposed by an aggregation strategy, which makes the usual methods (e.g., [1]) of analyzing and guaranteeing the closedloop stability inapplicable, Refs. [6,7] introduce a time-varying structure named the move-blocking technique. Meanwhile, Refs. [8,9] suggest to relax some constraints to enlarge the feasible region but at the cost of losing guaranteed closed-loop stability. On the other hand, Refs. [10,11] propose the equivalent/quasiequivalent aggregation concept to reduce the degradation in the control performance by utilizing the characteristic of MPC that only the first control input is actually implemented. However, a design of an aggregation-based MPC algorithm that reduces all these degradations remains unavailable. In this paper, an aggregation-based MPC algorithm is developed, starting from a given MPC algorithm, that can provide the required performance and attractive region and guarantee closedloop stability. Meanwhile, the quasi-equivalent aggregation strategy is developed to be embedded into the stable algorithm to reduce the degradation in control performance with guaranteed closed-loop stability. Furthermore, by relaxing the state constraints with a sequence of reachable sets, a switching strategy is adopted to enlarge the attractive region of the resulting aggregation-based MPC.
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089
The remainder of this paper is organized as follows. Section 2 contains the problem statement. The aggregation-based MPC algorithm with guaranteed stability and embedded quasi-equivalent aggregation strategy is developed in Section 3. Section 4 proposes a switching strategy to achieve an enlarged attractive region. A numerical example in Section 5 demonstrates the effectiveness of the proposed aggregation-based MPC algorithms. Notation. In the paper, we will use standard notion. Also, x(k + i|k) denotes the value of x at time k + i, predicted at time k. For a given Q ≻ 0, ∥x∥2Q := xT Qx. For brevity, x(k|k) = x(k). In addition, 0 and I are respectively the zero and identity matrices of appropriate dimensions.
1083
F as the optimal feedback control gain of an LQR problem and calculate its corresponding polyhedral invariant set as Xf . For (5), J (k) can be equivalently written as J (k) = X T (k)QX (k) + U T (k)RU (k), where Q = diag{Q , . . . , Q , QN }, R = diag{R , . . . , R }, X (k) = [xT(k + 1|k), xT(k + 2|k), . . ., xT(k + N |k)]T = Sx(k) + GU (k), and
where A ∈ R ,B ∈ R and (A, B) is stabilizable. The control inputs and the states are constrained by
matrices S and G can be obtained according to the system model (1). For a MPC controller with optimization problem (5), its control performance J∞ (0) is determined by chosen QN , F , and N. For the given weighted matrices (Q , R ) and designed terminal items (Xf , QN ), increasing the control horizon is a familiar method to improve the control performance and attractive region of MPC. But the online computational burden will be also increased as it is mainly determined by the number of optimization variables Nm. Ref. [3] suggests to adopt an aggregation strategy to reduce the number of optimization variables so as to reduce the online computational burden. This aggregation strategy is characterized by the following linear mapping from U (k) to V (k),
u(k) ∈ U := {u|Gu u ≤ 1},
(2)
U (k) = HV (k),
x(k) ∈ X := {x|Gx x ≤ 1},
(3)
where H ∈ R , V (k) = [v1 (k), v2 (k), . . . , vs (k)] ∈ R (m ≤ s ≪ Nm) is the new optimization vector with a much lower dimension than that of U (k). Each vi (k) is a scalar and is called an aggregation variable. Matrix H is called the aggregation matrix and is of full column rank to make full use of all aggregation variables. Here, we make the assumption of s ≥ m to design the equivalent/quasiequivalent aggregation matrix in the following part. As mentioned earlier, while the aggregation strategy greatly reduces the online computational burden, it may also lead to significant performance degradations. Our goal is to propose an MPC algorithm with carefully chosen aggregation matrix H to reduce such degradations.
2. Problem statement Consider the following system x(k + 1) = Ax(k) + Bu(k), n×n
(1)
n×m
where Gx and Gu are constant matrices, and ‘‘≤’’ defines an elementwise inequality. For system (1) subject to (2) and (3), the control goal is to steer the system state to the origin. Accordingly, the optimal control problem can be formulated as min
u(0),u(1),...
J∞ (0),
s.t. (1)–(3), k = 0, 1, . . . , ∞
(4)
2 2 where J∞ (0) = k=0 ∥x(k)∥Q + ∥u(k)∥R with given weighted matrices Q ≻ 0, R ≻ 0. Solving the above problem can give an LQR ∗ control law and achieve the optimal control performance J∞ (0). As pointed out by [1], the essence of MPC is similar to that of optimal control except for the implementation strategy, i.e., the receding horizon optimization. Generally speaking, the MPC technique substitutes the following finite horizon problem for the infinite horizon optimization problem (4) by dual-mode method [1] and then obtains the control law by online solving it at each time k, that is, by solving
∞
min J (k), U (k)
s.t. (1)–(3), x(k + N |k) ∈ Xf ,
(5)
where U (k) = [uT (k|k), uT (k + 1|k), . . . , uT (k + N − 1|k)]T , and J (k) =
N −1
∥x(k + i|k)∥2Q + ∥u(k + i|k)∥2R + ∥x(k + N |k)∥2QN
Nm×s
(9) T
s
3. The stable MPC algorithm with a quasi-equivalent aggregation strategy Although the move-blocking technique in [6] guarantees stability of the aggregation-based MPC, it does not take control performance into consideration. Since the concept of equivalent/ quasi-equivalent aggregation strategy in [10] complies with the goal of reducing the degradation on control performance, it is a natural idea to combine these two techniques. However, the timevarying structure does not combine with the equivalent aggregation strategy. Therefore, in this section, we will resort to the technique of including the tail [7] to develop a stable aggregationbased MPC algorithm. On top of this MPC algorithm, the quasiequivalent aggregation strategy is developed and embedded. 3.1. The stable aggregation-based MPC
i =0
is the objective function, with N being the length of the optimization horizon, QN being the terminal weighting matrix, and Xf (0 ∈ Xf ) being the terminal set with the local control law u(k) = Fx(k). Without loss of generality, we will assume Xf to be a polyhedral set {x|As x ≤ 1} such that Gx x ≤ 1
and Gu Fx ≤ 1,
(A + BF )x ∈ Xf
if x ∈ Xf ,
if x ∈ Xf ,
(7)
QN − (A + BF ) QN (A + BF ) ≽ Q + F R F . T
(6)
T
(8)
According to [1], under the above conditions, the closed-loop system with the MPC based on (5) is stable. The set Xf and matrices QN and F can be designed by Algorithm 1 in [12] or we can choose
The algorithm of stable aggregation-based MPC is presented as follows. Algorithm 1. Step 1. Let U0 (0) = 0. Step 2. At time k, obtain V ∗ (k) and a∗ by solving the optimization problem (5) with U (k) = aU0 (k) + HV (k), where a is an additional optimization variable. Step 3. Act [I , 0, . . . , 0](a∗ U0 (k) + HV ∗ (k)) on the plant. Step 4. Let U0 (k + 1) = [(u∗ (k + 1|k))T , (u∗ (k + 2|k))T , . . . , (u∗ (k + N − 1|k))T , (Fx∗ (k + N |k))T ]T , where u∗ (k + 1|k), u∗ (k + 2|k), . . . , u∗ (k+N −1|k) are elements of U ∗ (k) = a∗ U0 (k)+HV ∗ (k) and x∗ (k + N |k) is the corresponding terminal state. Return to Step 2.
1084
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089
The main idea of Algorithm 1 is to replace U (k) with U (k) = aU0 (k)+ HV (k), where U0 (k) is a time-varying vector referred to as the base vector. When a = 0, Algorithm 1 reduces to a traditional aggregation-based MPC. Theorem 1. Consider Algorithm 1. If the optimization problem in Step 2 is feasible at time k with x(k), then the resulting closed-loop system is asymptotically stable.
By comparing Eq. (11) with (12), we get H T GT QG + R (HV (k) − Uˆ (k)) = 0.
(13) T
T
Since each of the first m columns of matrix H G QG + R is linearly independent of the other Nm − 1 columns, the first m elements of solution HV (k)−Uˆ (k) to (13) must be zero. This concludes the proof.
Proof. The optimization problem in Step 2 of Algorithm 1 is feasible at time k. Let V ∗ (k) and a∗ be the optimal solution. Then, the control input sequence is given by
Following a similar proof of Theorem 2, we can arrive at the following corollary.
U ∗ (k) = a∗ U0 (k) + HV ∗ (k)
Corollary 1. Consider Algorithm 1 with no constraints. If each of the first s columns of matrix H T GT QG + R is linearly independent of the other Nm − 1 columns, then the first s elements of HV (k) must be equal to those of the original MPC.
T := u∗ T (k|k), u∗ T (k + 1|k), . . . , u∗ T (k + N − 1|k) .
Also, let the corresponding sequence of system states be [x∗ (k + 1|k), x∗ (k + 2|k), . . . , x∗ (k + N |k)]T and the cost value be J ∗ (k). Step 4 of Algorithm 1 yields U0 (k + 1) = (u∗ (k + 1|k))T , (u∗ (k + 2|k))T , . . . ,
The following algorithm can be given to construct matrix H ∈ RNm×s that satisfies conditions in Theorem 2 and Corollary 1 by referring to [10].
(u∗ (k + N − 1|k))T , (Fx∗ (k + N |k))T
T
.
Since x∗ (k + N |k) ∈ Xf , according to (6)–(8), (a, V (k + 1)) = (1, 0) is a feasible solution at time k + 1. The result of the theorem then follows from [1]. 3.2. The embedded quasi-equivalent aggregation strategy
Lemma 1. Consider the optimization problem (5) with x(k). Let U (k) = U¯ (k) + Uˆ (k), where U¯ (k) is any known vector and Uˆ (k) is the optimization variable, the optimal solution U ∗ (k) = U¯ (k) + Uˆ ∗ (k) equals to the optimal solution of the original optimization problem. Definition 1 ([10]). If an aggregation-based MPC gives the same control input for state x(k) at time instant k as the original MPC, the aggregation strategy is called an equivalent aggregation strategy and the aggregation matrix is called the equivalent aggregation matrix. Theorem 2. ConsiderAlgorithm 1 with no constraints. If each of the first m columns of H T GT QG + R is linearly independent of the other Nm − 1 columns, then H is an equivalent aggregation matrix. Proof. Since Lemma 1 implies that the time-varying base vector does not influence the result of the original MPC, to prove the theorem, we need only to show that, if the condition of the theorem holds, Algorithm 1 with the given H results in the same control inputs as the original one with the time-varying base vector. In the absence of constraints, the optimization problem of the original MPC with U (k) = aU0 (k) + Uˆ (k) is Uˆ (k),a
(10)
Denote aU0 (k) as U¯ 0 (k). For the original MPC, its optimal solution ∂ J (k) can be obtained from ˆ = 0 as ∂ U (k)
(GT QG + R)Uˆ (k) = −GT QSx(k) − GT QGU¯ 0 (k) − RU¯ 0 (k).
(11)
Similarly, the solution of the aggregation-based MPC is
(GH )T Q (GH ) + H T RH V (k) = −(GH )T QSx(k) − (GH )T QGU¯ 0 (k) − H T RU¯ 0 (k).
Step 3. Divide W2 into W2 = (Nm−s)×(Nm−s)
W20 W21
, where W20 ∈ Rs×(Nm−s) and
W21 ∈ R . Represent the s row vectors of W20 , with the row number ij , j = 1, 2, . . . , s, as the linear combinations of
(12)
T
the rows in W21 with coefficients λj = λ1j , λ2j , . . . , λ(Nm−s)j , i.e., [λ1 λ2 · · · λs ]T W21 = W20 . Step 4. Construct matrix Λ = [−λ1 ,−λ 2 , . . . , −λs ], where λj ’s
In this subsection we will develop the quasi-equivalent aggregation strategy and embed it in Algorithm 1. We first state a lemma, whose proof is straightforward and is thus omitted, and recall the definition of the equivalent aggregation strategy from [10].
min J (k) = X T (k)QX (k) + U T (k)RU (k).
Algorithm 2. Step 1. Choose s and let W = GT QG + R. Step 2. Divide W into W = [W1 W2 ], where W1 ∈ RNm×s and W2 ∈ RNm×(Nm−s) .
I
are obtained in Step 3. Obtain H =
Λ
, where I is an identity
matrix of dimension s. Remark 1. In Algorithm 2, since W is positive definite,W21 must be of full rank. Then, since H T W = H T [W1 W2 ] = H T W1 0 according to the structure of obtained H, rank(H ) = s ≥ m, and W is of full rank, Algorithm 2 yields an equivalent aggregation matrix H for the unconstrained case by Theorem 2. The above algorithm includes the result in [10] with s = m as a special case. systems, if an aggregation strategy with H = For constrained H1 H12 ( H ∈ Rm×m ) is applied, we can find that u(k) = 1 H21 H2 H1 v1 (k)+ H12 v2 (k) and [uT (k + 1), uT (k + 2), . . . , uT (k + N − 1)]T = H21 v1 (k)+ H2 v2 (k), where V (k) = [v1T (k) v2T (k)]T . It is obvious that the current control input is determined by both v1 (k) and v2 (k). Since the constraints are about the whole control input sequence, the constraints on [uT (k + 1), uT (k + 1), . . . , uT (k + N − 1)]T are also imposed on the input u(k) by the structure of the aggregation matrix through v2 (k). In comparison with the original MPC, if the additional relationship on the first control input is removed, the restrictions on the first control input will be closer to the original MPC. Hence, we introduce the following transformation on the matrix H resulting from Algorithm 2
H0 H = H1
∈R
Nm×s
I H⇒ H = 0
0 H1
∈ RNm×(s+l) ,
(14)
where H0 ∈ Rl×s (l ≥ m) and I is an l × l identity matrix. It can be easily proven that if the resulting matrix H from the transformation (14) is of full-rank, it also satisfies the conditions in Theorem 2. Since the transformation in (14) increases the degree of freedom and guarantees the equivalent feature for the unconstrained case, we call the resulting matrix from (14) the quasi-equivalent aggregation matrix. Applying this quasi-equivalent aggregation matrix to Algorithm 1 results in the following algorithm.
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089
Algorithm 3. Off-line part: Choose s and l, use Algorithm 2 to design the equivalent aggregation matrix, and then obtain the quasiequivalent aggregation matrix H by (14). Online part: Carry out Steps 2–5 in Algorithm 1 with an H obtained in the off-line part. 4. Enlargement of the attractive region of aggregation-based MPC To enlarge the attractive region of the aggregation-based MPC, Ref. [9] suggests to relax the constraints of the optimization problem at the cost of loss of guaranteed closed-loop stability. In this section, we will develop an aggregation-based MPC algorithm based on Algorithm 3 and a switching strategy to achieve both an enlarged attractive region and guaranteed closed-loop stability. According to Theorem 1, the attractive region resulting from Algorithm 3 is the feasible region of the optimization problem in the online part with U0 (k) = 0, which can be calculated by linear difference inclusions [13] and is denoted as C0 . For C0 , its one-step set C1 is defined as a set such that for x(k) ∈ C1 , there is an admissible control input u(k) such that x(k + 1) ∈ C0 . Similarly, we can define the one-step set of Cj as Cj+1 . According to [14,13], we can compute off-line a sequence of sets {C0 , C1 , . . . , Ch } one by one from C0 . Here, h is the number of one-step reachable sets, which can be chosen according to applications and a larger h may correspond to a larger attractive region at the cost of a large amount of required memory. Since U0 (k) in Algorithm 3 can be not equal to zero, C0 may not be an invariant set for the algorithm in the online part of Algorithm 3. Hence, Cj cannot be guaranteed to belong to Cj+1 , which makes the switching design in [14] inapplicable. Here, we give the following MPC algorithm based on the aggregation matrix H obtained by the offline part of Algorithm 3 and the sets {C0 , C1 , . . . , Ch } calculated off-line. Algorithm 4. Step 1. At k = 0, U0 (0) = 0 and p(0) = 0. Step 2. When k > 0, check {x(k), x(k + N − 1|k − 1), p(k − 1)}. For case with x(k) ∈ C0 (k ≥ 0) or case with x(k + N − 1|k − 1) ∈ Xf and p(k − 1) = 1, set p(k) = 1. Otherwise, set p(k) = 0 and choose the minimum j such that x(k) ∈ Cj (h ≥ j > 0). Step 3. If p(k) = 1, calculate the optimal solution U ∗ (k) = a∗ U0 (k) + HV ∗ (k) by Step 2 of Algorithm 1 and then go to Step 5. Step 4. If p(k) = 0, solve the following optimization problem to obtain the optimal solution U ∗ (k) = a∗ U0 (k) + HV ∗ (k), min J (k),
V (k),a
(15)
s.t. u(k + i|k) ∈ U , i = 0, 1, . . . , N − 1, x(k + 1|k) ∈ Cj−1 , U (k) = aU0 (k) + HV (k). Step 5. Execute [I , 0, . . . , 0]U ∗ (k) on the controlled plant. Step 6. Calculate x(k + N |k) according to U ∗ (k), let U0 (k + 1) = [(u∗ (k + 1|k))T , (u∗ (k + 2|k))T , . . . , (u∗ (k + N − 1|k))T , (Fx∗ (k + N |k))T ]T . Set k = k + 1 and return to Step 2. Remark 2. In Algorithm 4, the constraints on system states other than x(k + 1|k) and the terminal constraints are relaxed in Step 4. Hence, the feasible region of the above algorithm, i.e., C0 ∪ C1 ∪ · · ·∪ Ch , must be larger than that of Algorithm 3. Since checking the switching conditions in Step 2 involves only algebraic calculation without iteration, the online computational burden of Algorithm 4 is much lighter than that of the original MPC. Theorem 3. Consider the aggregation-based MPC Algorithm 4. If the algorithm is feasible at time k for x(k), then the closed-loop system is asymptotically stable.
1085
Proof. At time k, Algorithm 4 is feasible. There are three cases to be considered. Case 1: When p(k) = 0, x(k) ∈ Cj (h ≥ j > 0) and optimization problem (15) is implemented. Since the algorithm is feasible, there exists an admissible control input u(k) to steer x(k) to Cj−1 . Since Cj−1 is also the one-step set of Cj−2 , if j − 2 > 0, the optimization problem (15) is feasible at the next time. This procedure will continue until x(k) ∈ C0 and the controller switches to the case with p(k) = 1. Case 2: When x(k) ∈ C0 , p(k) = 1 and Step 3 is implemented. Since the algorithm is feasible, the optimal control inputs will steer x(k + N |k) into Xf , i.e., p(k) will be kept as 1. Thus, at time k + 1, Step 3 is still implemented. Then, according to Theorem 1, the closedloop system is asymptotically stable. Case 3: If x(k + N − 1|k − 1) ∈ Xf and p(k − 1) = 1, Step 3 of Algorithm 4 is implemented. Since p(k − 1) = 1, the feasibility of Algorithm 4 at time k means that U0 (k + 1) obtained in Step 5 satisfies all the constraints. So it is a feasible solution at time k + 1 with V (k) = 0 and a = 1 and p(k) will be kept as 1. Then, according to Theorem 1, the closed-loop system is asymptotically stable. The above analysis concludes the proof. Remark 3. In the proof of Theorem 3, we see that, although C0 may not be an invariant set and Cj cannot be guaranteed to belong to Cj+1 , the condition x(k + N − 1|k − 1) ∈ Xf (when k > 0) and p(k − 1) = 1 in Step 2 of Algorithm 4 avoid the possible cyclic switching among several sets of {C0 , C1 , . . . , Ch }. This guarantees the closed-loop stability of Algorithm 4 and is the main difference from the previous design in [14]. 5. Numerical simulation This section will verify the results through two examples. The first one illustrates the effectiveness by a mathematical case. The second one considers a real world system (i.e., a two open channel network system) to verify the practical value of the proposed method. 5.1. Case 1 We first consider the following unstable system with two inputs, x( k + 1 ) =
0.8 −0.2
1 0.1 x(k) + 1.1 −0.2
0.2 u(k). 0.5
The constraints are |ui (k)| ≤ 1, i = 1, 2, and |x2 (k)| ≤ 4. Let Q and R be an identity matrix and N = 10. In order to achieve a well-designed MPC, the terminal set and QN are designed in accordance with the feedback law resulting from an LQR design. Set s = m. Then, by the off-line design of Algorithm 3 with l = 2 and l = 4, respectively, we obtain two quasi-equivalent aggregation matrices. That is, the numbers of aggregation variables are respectively 4 and 6. Shown in Fig. 1 and enclosed by the solid curve is the attractive region resulting from the MPC Algorithm 4 with h = 7 and l = 2, i.e., C7 . Meanwhile, Corg is the attractive region of the original MPC. It is obvious that C0 , which is the feasible region of Algorithm 3, is much smaller than Corg due to the reduction of freedom caused by the aggregation strategy. Comparatively, the attractive region is enlarged by Algorithm 4. For clarity, we have not plotted all the sets C0 , C1 , . . . , C7 in Fig. 1. Meanwhile, in Fig. 1, we see that C0 does not strictly belong to C2 , which verifies the discussion in Section 4. Under Algorithm 4, the state trajectories starting from several initial conditions (i.c.) near the boundary of C7 are plotted in dotted lines in Fig. 2, which shows that the closed-loop system is stable and thus verifies the conclusion of Theorem 3.
1086
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089 Table 1 Comparisons of control performance (percentage difference from that of original MPC).
Avg. cost (%) max cost (%) min cost (%) # i.c. (< 0.1%)
l=2
l=4
[9] with H1
[9] with H2
2.13 7.47 0.05 1
0.70 5.48 0 4
7.49 22.64 1.76 0
3.37 13.72 0.14 0
Table 2 Comparisons of the number of online optimization variables. Original MPC
l=2
l=4
[9] with H1
[9] with H2
20
5
7
6
8
respectively, the algorithm in [9] with I 0
Fig. 1. The attractive region.
.
H1 = .. 0 0
Fig. 2. The states and trajectories.
0 I
.. .
I 0
0 0
.. 20×6 , .∈R 0 I
I 0 0
H2 = ..
. 0 0
0 I 0
0 0 I
0 0
I 0
.. .
.. .
0 0 0
20×8 , .. ∈R . 0
I
respectively. In addition, the original MPC, which achieves the same control performance as the MPC with the infinite control horizon, is used here as a benchmark. The comparison of the accumulated cost indices is shown in Fig. 3 and Table 1. The comparison of the online optimization variables, which largely determine the online computational burden, is given in Table 2. In Fig. 3 and Table 1, results are shown for 17 randomly selected initial conditions. These conditions are marked in Fig. 2 by circle and cross (state [5, −3]T , the sixth point in Fig. 3), respectively. The result obtained from the original MPC is denoted as 100%. In Fig. 3, we can see that the proposed aggregation-based MPC with l = 4 achieves almost the same control performance as the original MPC. The control performance of the proposed aggregation-based MPC with l = 2 is also better than that of the algorithm in [9] with H1 and H2 for most initial conditions. Table 1 also reflects the same conclusion as Fig. 3. In summary, the proposed aggregation-based MPCs achieve better performance with lower online computational burden. In order to illustrate the comparison in more detail, the state responses and control inputs for initial state [5, 0]T are shown in Figs. 4 and 5. From Fig. 5, it is obvious that the control inputs resulting from the proposed MPC are similar as those resulting from the original MPC, especially for the case with l = 4, which leads to the similar state responses in Fig. 4. This owes to the similar restrictions on the control input sequence of Algorithm 4 and equivalent feature of the adopted aggregation strategy. Moreover, the degradation on control performance is reduced as the values of l increases. This comparison reinforces the comparison made in Fig. 3 and Table 1, and again validates the effectiveness of the proposed aggregation-based MPC algorithm. 5.2. Case 2
Fig. 3. Comparison of the accumulated index.
To compare the control performance, we consider two aggregation-based MPCs proposed in this paper with l = 2 and l = 4,
Consider a canal network system with two reaches in [15–17]. The structure of the canal system is shown in Fig. 6, where Hup = 2.5 m, Hdo = 1.14 m, the bed slope is 0.0003, manning roughness coefficient is taken as n = 0.015, and the third gate is fixed at U3 = 0.85 m2 . The parameters of each reach and gate are given in Tables 3 and 4, respectively. According to [16], for the canal system in Fig. 6, its mathematical model around a reference condition of uniform flow can been derived by first linearizing the Saint-Venant equations for the
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089
Fig. 4. The state responses.
Fig. 5. The control inputs.
Fig. 6. The structure of canal network system.
1087
1088
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089
Table 3 Canal parameters. Reach No.
K1i
K3i
W1i
W2i
1 2
3.76×10−5 2.85×10−5
6.38×10−5 7.51×10−5
6.5386 × 10−2 7.0538 × 10−2
−5.8505×10−2 −10.018×10−2
Table 4 Gate parameters. Gate No.
ai
1 2 3
2.41 1.91
bi
ci
−3.76 −2.99
2.99 5.08
unsteady flow of water in open-channels. With consideration of the variations of the user flow rates p, the model can be written as
v˙ (t ) = A1 v(t ) + B1 u(t ) − Ip(t )
Fig. 7. The user flow rate.
where v(t ) = [v1 (t ), v2 (t )]T is the vector of storage volume variations, u(t ) = [u1 (t ), u2 (t )]T is the vector of gate opening size variations, and p(t ) = [p1 (t ), p2 (t )]T is the vector of user flow rates. For details about the procedure to achieve (A1 , B1 ) and the physical meaning of parameters in Tables 3 and 4, refer to [15–17]. According to the continuous-time model of the canal system, by choosing the sampling time as T = 5 s, we can obtain the discretetime model as
v(k + 1) =
0.9977 0.0019
+
0.0012 v(k) 0.9966
15.1211 −3.7230
−13.9767 5 u(k) − 16.9188 0
0 p(k). 5
The input constraints are −0.3 ≤ u1 ≤ 0.3 and −0.15 ≤ u2 ≤ 0.15. Here, we choose N = 200 to approximate the MPC with infinite horizon and do not add the terminal set constraints to avoid the infeasibility of the optimization problem in the real world practice. Another reason for choosing a long horizon of MPC is that the future user flow rate can be known or scheduled in real world applications and MPC with a long horizon can make full use of the known user flow rate to improve the performance. For the original MPC, the number of online optimization variables is 400. Meanwhile, according to [16], the weighting matrices are chosen as Q = diag{1, 1.2233} and R = diag{3 × 105 , 3 × 105 }. In addition, the terminal weighted matrix is adopted to guarantee the closed-loop stability. By choosing s = 2 and l = 0, we can get an MPC controller based on Algorithm 3. Note that, since only system input variables are subject to the constraints, the controlled system is stable, and no terminal constraints are considered, the feasible region of Algorithm 3 based MPC is as large as the original one with infinite horizon. The aggregation strategy reduces the number of online optimization variables from 400 of the original MPC to 3 of Algorithm 3 based MPC. The control results and control inputs from the initial state [10, 22]T are shown in Figs. 8 and 9. Fig. 7 shows the variation of user flow rate. From these figures, we can easily see that, although the number of online optimization variables of aggregation based MPC is much smaller than the original one, the control performance is almost the same as the original MPC. This means that, by the proposed approach, we can achieve a similar control performance as the original MPC but at a much less online computation cost.
Fig. 8. The control inputs of the original MPC and aggregation based MPC.
6. Conclusions Online computational burden has always been a critical issue for the practical application of MPC algorithms. Although the aggregation strategy can reduce the online computational burden,
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089
1089
Acknowledgments The work was supported in part by the National Natural Science Foundation of China (Grant Nos. 61333009, 61104078, 61074060, 60934007, and 61221003) and the Specialized Research Fund for the Doctoral Program of Higher Education (Grant No. 20120073110017). References
Fig. 9. The control results of the original MPC and aggregation based MPC.
it may lead to loss of guaranteed closed-loop stability, poor control performance and a smaller attractive region. This paper proposed an aggregation-based MPC algorithm that addresses these issues associated with the aggregation strategy by combining the timevarying base vector, the quasi-equivalent aggregation matrix, and a switching strategy with relaxed constraints.
[1] D.Q. Mayne, J.B. Rawlings, C.V. Rao, P.O. Scokaert, Constrained model predictive control: stability and optimality, Automatica 36 (6) (2000) 789–814. [2] S.J. Qin, T.A. Badgwell, A survey of industrial model predictive control technology, Control Engineering Practice 11 (7) (2003) 733–764. [3] D. Li, Y. Xi, The general framework of aggregation strategy of model predictive control and stability analysis, in: Proceedings of the 11th IFAC Symposium on Large Scale Systems Theory and Applications, 2007, pp. 192–197. [4] N.L. Ricker, Use of quadratic programming for constrained internal model control, Industrial & Engineering Chemistry Process Design and Development 24 (4) (1985) 925–936. [5] J. Richalet, S. Abu el Ata-Doss, Ch. Arber, H.B. Kuntze, A. Jacubasch, W. Schill, Predictive functional control: application to fast and accurate robots, in: Proceedings of the 10th IFAC Congress, 1987, pp. 251–258. [6] R. Cagienard, P. Grieder, E.C. Kerrigan, M. Morari, Move blocking strategies in receding horizon control, Journal of Process Control 17 (6) (2007) 563–570. [7] G. Valencia-Palomo, M. Pelegrinis, J.A. Rossiter, R. Gondhalekar, A moveblocking strategy to improve tracking in predictive control, in: Proceeding of American Control Conference, ACC 2010, 2010, pp. 6293–6298. [8] F. Oldewurtel, R. Gondhalekar, C.N. Jones, M. Morari, Blocking parameterizations for improving the computational tractability of affine disturbance feedback MPC problems, in: Proceeding of the 48th IEEE Conference on Decision and Control, 2009, pp. 7381–7386. [9] R. Gondhalekar, J. Imura, Least-restrictive move-blocking model predictive control, Automatica 46 (7) (2010) 1234–1240. [10] D. Li, Y. Xi, Quality guaranteed aggregation based model predictive control and stability analysis, Science in China Series F: Information Sciences 52 (7) (2009) 1145–1156. [11] D. Li, Y. Xi, Aggregation based closed-loop MPC with guaranteed performance, in: The 48th IEEE Conference on Decision and Control, 2009, pp. 7400–7405. [12] D. Li, Y. Xi, The synthesis of robust model predictive control with qp formulation, International Journal of Modelling, Identification and Control 13 (1) (2011) 1–8. [13] F. Blanchini, Set invariance in control, Automatica 35 (11) (1999) 1747–1767. [14] D. Limon, T. Alamo, E.F. Camacho, Enlarging the domain of attraction of MPC controllers, Automatica 41 (4) (2005) 629–635. [15] C. Seatzu, Decentralized controllers design for open-channel hydraulic systems via eigenstructure assignment, Applied Mathematical Modelling 24 (12) (2000) 915–930. [16] G. Corriga, S. Sanna, G. Usai, Sub-optimal constant-volume control for open channel networks, Applied Mathematical Modelling 7 (4) (1983) 262–267. [17] G. Corriga, S. Sanna, G. Usai, Estimation of uncertainty in an open-channel network mathematical model, Applied Mathematical Modelling 13 (11) (1989) 651–657.
Contents lists available at ScienceDirect
Systems & Control Letters journal homepage: www.elsevier.com/locate/sysconle
An improved design of aggregation-based model predictive control Dewei Li a,∗ , Yugeng Xi a , Zongli Lin b,a a
Department of Automation, Shanghai Jiao Tong University, Key Laboratory of System Control and Information Processing, Ministry of Education, Shanghai, 200240, China b
Charles L. Brown Department of Electrical and Computer Engineering, University of Virginia, Charlottesville, VA 22904-4743, USA
article
info
Article history: Received 28 January 2013 Received in revised form 12 August 2013 Accepted 27 August 2013 Available online 30 September 2013 Keywords: Model predictive control Aggregation Equivalent aggregation Quasi-equivalent aggregation Blocking
abstract Linear aggregation in the input is an effective method to reduce the online computational burden of model predictive control (MPC) but at the cost of degradations in the closed-loop performance. In this paper, an improved aggregation-based MPC algorithm is developed to reduce these degradations. In this algorithm, a time-varying base vector is utilized in conjunction with the quasi-equivalent aggregation strategy. Furthermore, by relaxing the constraints with a sequence of reachable sets, a switching strategy is adopted to enlarge the attractive region of the resulting aggregation-based MPC. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Model predictive control (MPC), also referred to as receding horizon control, emerged in 1970s. Because of its ability to tackle the constraints explicitly and to achieve good control performance, MPC has been attracting more and more attention from both the industrial and academic communities (see, e.g., [1,2]). MPC is characterized by rolling horizon optimization. An optimal control problem over a finite horizon is solved online at each sampling time and the first optimal control input is executed on the controlled plant. This feature leads to several advantages of MPC over other control algorithms but at the cost of heavy online computational burden. Increasing the number of optimization variables, i.e., the control inputs to be optimized, may further improve the control performance due to the increased degree of freedom, but the demand on online computation will also further intensify. As a result, the online computational burden has always been an important issue in the practical application of MPC. Input linear aggregation [3], also known as input parametrization [2], has proven to be an effective strategy to reduce the online computational burden of an MPC algorithm. The idea can be traced back to the blocking technology in [4], where the control horizon is divided into several blocks and the control inputs within each block
∗
Corresponding author. Tel.: +86 21 34204550. E-mail addresses: [email protected] (D. Li), [email protected] (Y. Xi), [email protected] (Z. Lin). 0167-6911/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.sysconle.2013.08.007
are assumed to be the same to reduce the number of online optimization variables. Similar strategies include the Predictive Functional Control (PFC) in [5]. However, the aggregation strategy also leads to performance degradations, including poor control performance, loss of guaranteed closed-loop stability and reduction in size of the attractive region. In order to overcome the difficulty caused by the structural restriction imposed by an aggregation strategy, which makes the usual methods (e.g., [1]) of analyzing and guaranteeing the closedloop stability inapplicable, Refs. [6,7] introduce a time-varying structure named the move-blocking technique. Meanwhile, Refs. [8,9] suggest to relax some constraints to enlarge the feasible region but at the cost of losing guaranteed closed-loop stability. On the other hand, Refs. [10,11] propose the equivalent/quasiequivalent aggregation concept to reduce the degradation in the control performance by utilizing the characteristic of MPC that only the first control input is actually implemented. However, a design of an aggregation-based MPC algorithm that reduces all these degradations remains unavailable. In this paper, an aggregation-based MPC algorithm is developed, starting from a given MPC algorithm, that can provide the required performance and attractive region and guarantee closedloop stability. Meanwhile, the quasi-equivalent aggregation strategy is developed to be embedded into the stable algorithm to reduce the degradation in control performance with guaranteed closed-loop stability. Furthermore, by relaxing the state constraints with a sequence of reachable sets, a switching strategy is adopted to enlarge the attractive region of the resulting aggregation-based MPC.
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089
The remainder of this paper is organized as follows. Section 2 contains the problem statement. The aggregation-based MPC algorithm with guaranteed stability and embedded quasi-equivalent aggregation strategy is developed in Section 3. Section 4 proposes a switching strategy to achieve an enlarged attractive region. A numerical example in Section 5 demonstrates the effectiveness of the proposed aggregation-based MPC algorithms. Notation. In the paper, we will use standard notion. Also, x(k + i|k) denotes the value of x at time k + i, predicted at time k. For a given Q ≻ 0, ∥x∥2Q := xT Qx. For brevity, x(k|k) = x(k). In addition, 0 and I are respectively the zero and identity matrices of appropriate dimensions.
1083
F as the optimal feedback control gain of an LQR problem and calculate its corresponding polyhedral invariant set as Xf . For (5), J (k) can be equivalently written as J (k) = X T (k)QX (k) + U T (k)RU (k), where Q = diag{Q , . . . , Q , QN }, R = diag{R , . . . , R }, X (k) = [xT(k + 1|k), xT(k + 2|k), . . ., xT(k + N |k)]T = Sx(k) + GU (k), and
where A ∈ R ,B ∈ R and (A, B) is stabilizable. The control inputs and the states are constrained by
matrices S and G can be obtained according to the system model (1). For a MPC controller with optimization problem (5), its control performance J∞ (0) is determined by chosen QN , F , and N. For the given weighted matrices (Q , R ) and designed terminal items (Xf , QN ), increasing the control horizon is a familiar method to improve the control performance and attractive region of MPC. But the online computational burden will be also increased as it is mainly determined by the number of optimization variables Nm. Ref. [3] suggests to adopt an aggregation strategy to reduce the number of optimization variables so as to reduce the online computational burden. This aggregation strategy is characterized by the following linear mapping from U (k) to V (k),
u(k) ∈ U := {u|Gu u ≤ 1},
(2)
U (k) = HV (k),
x(k) ∈ X := {x|Gx x ≤ 1},
(3)
where H ∈ R , V (k) = [v1 (k), v2 (k), . . . , vs (k)] ∈ R (m ≤ s ≪ Nm) is the new optimization vector with a much lower dimension than that of U (k). Each vi (k) is a scalar and is called an aggregation variable. Matrix H is called the aggregation matrix and is of full column rank to make full use of all aggregation variables. Here, we make the assumption of s ≥ m to design the equivalent/quasiequivalent aggregation matrix in the following part. As mentioned earlier, while the aggregation strategy greatly reduces the online computational burden, it may also lead to significant performance degradations. Our goal is to propose an MPC algorithm with carefully chosen aggregation matrix H to reduce such degradations.
2. Problem statement Consider the following system x(k + 1) = Ax(k) + Bu(k), n×n
(1)
n×m
where Gx and Gu are constant matrices, and ‘‘≤’’ defines an elementwise inequality. For system (1) subject to (2) and (3), the control goal is to steer the system state to the origin. Accordingly, the optimal control problem can be formulated as min
u(0),u(1),...
J∞ (0),
s.t. (1)–(3), k = 0, 1, . . . , ∞
(4)
2 2 where J∞ (0) = k=0 ∥x(k)∥Q + ∥u(k)∥R with given weighted matrices Q ≻ 0, R ≻ 0. Solving the above problem can give an LQR ∗ control law and achieve the optimal control performance J∞ (0). As pointed out by [1], the essence of MPC is similar to that of optimal control except for the implementation strategy, i.e., the receding horizon optimization. Generally speaking, the MPC technique substitutes the following finite horizon problem for the infinite horizon optimization problem (4) by dual-mode method [1] and then obtains the control law by online solving it at each time k, that is, by solving
∞
min J (k), U (k)
s.t. (1)–(3), x(k + N |k) ∈ Xf ,
(5)
where U (k) = [uT (k|k), uT (k + 1|k), . . . , uT (k + N − 1|k)]T , and J (k) =
N −1
∥x(k + i|k)∥2Q + ∥u(k + i|k)∥2R + ∥x(k + N |k)∥2QN
Nm×s
(9) T
s
3. The stable MPC algorithm with a quasi-equivalent aggregation strategy Although the move-blocking technique in [6] guarantees stability of the aggregation-based MPC, it does not take control performance into consideration. Since the concept of equivalent/ quasi-equivalent aggregation strategy in [10] complies with the goal of reducing the degradation on control performance, it is a natural idea to combine these two techniques. However, the timevarying structure does not combine with the equivalent aggregation strategy. Therefore, in this section, we will resort to the technique of including the tail [7] to develop a stable aggregationbased MPC algorithm. On top of this MPC algorithm, the quasiequivalent aggregation strategy is developed and embedded. 3.1. The stable aggregation-based MPC
i =0
is the objective function, with N being the length of the optimization horizon, QN being the terminal weighting matrix, and Xf (0 ∈ Xf ) being the terminal set with the local control law u(k) = Fx(k). Without loss of generality, we will assume Xf to be a polyhedral set {x|As x ≤ 1} such that Gx x ≤ 1
and Gu Fx ≤ 1,
(A + BF )x ∈ Xf
if x ∈ Xf ,
if x ∈ Xf ,
(7)
QN − (A + BF ) QN (A + BF ) ≽ Q + F R F . T
(6)
T
(8)
According to [1], under the above conditions, the closed-loop system with the MPC based on (5) is stable. The set Xf and matrices QN and F can be designed by Algorithm 1 in [12] or we can choose
The algorithm of stable aggregation-based MPC is presented as follows. Algorithm 1. Step 1. Let U0 (0) = 0. Step 2. At time k, obtain V ∗ (k) and a∗ by solving the optimization problem (5) with U (k) = aU0 (k) + HV (k), where a is an additional optimization variable. Step 3. Act [I , 0, . . . , 0](a∗ U0 (k) + HV ∗ (k)) on the plant. Step 4. Let U0 (k + 1) = [(u∗ (k + 1|k))T , (u∗ (k + 2|k))T , . . . , (u∗ (k + N − 1|k))T , (Fx∗ (k + N |k))T ]T , where u∗ (k + 1|k), u∗ (k + 2|k), . . . , u∗ (k+N −1|k) are elements of U ∗ (k) = a∗ U0 (k)+HV ∗ (k) and x∗ (k + N |k) is the corresponding terminal state. Return to Step 2.
1084
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089
The main idea of Algorithm 1 is to replace U (k) with U (k) = aU0 (k)+ HV (k), where U0 (k) is a time-varying vector referred to as the base vector. When a = 0, Algorithm 1 reduces to a traditional aggregation-based MPC. Theorem 1. Consider Algorithm 1. If the optimization problem in Step 2 is feasible at time k with x(k), then the resulting closed-loop system is asymptotically stable.
By comparing Eq. (11) with (12), we get H T GT QG + R (HV (k) − Uˆ (k)) = 0.
(13) T
T
Since each of the first m columns of matrix H G QG + R is linearly independent of the other Nm − 1 columns, the first m elements of solution HV (k)−Uˆ (k) to (13) must be zero. This concludes the proof.
Proof. The optimization problem in Step 2 of Algorithm 1 is feasible at time k. Let V ∗ (k) and a∗ be the optimal solution. Then, the control input sequence is given by
Following a similar proof of Theorem 2, we can arrive at the following corollary.
U ∗ (k) = a∗ U0 (k) + HV ∗ (k)
Corollary 1. Consider Algorithm 1 with no constraints. If each of the first s columns of matrix H T GT QG + R is linearly independent of the other Nm − 1 columns, then the first s elements of HV (k) must be equal to those of the original MPC.
T := u∗ T (k|k), u∗ T (k + 1|k), . . . , u∗ T (k + N − 1|k) .
Also, let the corresponding sequence of system states be [x∗ (k + 1|k), x∗ (k + 2|k), . . . , x∗ (k + N |k)]T and the cost value be J ∗ (k). Step 4 of Algorithm 1 yields U0 (k + 1) = (u∗ (k + 1|k))T , (u∗ (k + 2|k))T , . . . ,
The following algorithm can be given to construct matrix H ∈ RNm×s that satisfies conditions in Theorem 2 and Corollary 1 by referring to [10].
(u∗ (k + N − 1|k))T , (Fx∗ (k + N |k))T
T
.
Since x∗ (k + N |k) ∈ Xf , according to (6)–(8), (a, V (k + 1)) = (1, 0) is a feasible solution at time k + 1. The result of the theorem then follows from [1]. 3.2. The embedded quasi-equivalent aggregation strategy
Lemma 1. Consider the optimization problem (5) with x(k). Let U (k) = U¯ (k) + Uˆ (k), where U¯ (k) is any known vector and Uˆ (k) is the optimization variable, the optimal solution U ∗ (k) = U¯ (k) + Uˆ ∗ (k) equals to the optimal solution of the original optimization problem. Definition 1 ([10]). If an aggregation-based MPC gives the same control input for state x(k) at time instant k as the original MPC, the aggregation strategy is called an equivalent aggregation strategy and the aggregation matrix is called the equivalent aggregation matrix. Theorem 2. ConsiderAlgorithm 1 with no constraints. If each of the first m columns of H T GT QG + R is linearly independent of the other Nm − 1 columns, then H is an equivalent aggregation matrix. Proof. Since Lemma 1 implies that the time-varying base vector does not influence the result of the original MPC, to prove the theorem, we need only to show that, if the condition of the theorem holds, Algorithm 1 with the given H results in the same control inputs as the original one with the time-varying base vector. In the absence of constraints, the optimization problem of the original MPC with U (k) = aU0 (k) + Uˆ (k) is Uˆ (k),a
(10)
Denote aU0 (k) as U¯ 0 (k). For the original MPC, its optimal solution ∂ J (k) can be obtained from ˆ = 0 as ∂ U (k)
(GT QG + R)Uˆ (k) = −GT QSx(k) − GT QGU¯ 0 (k) − RU¯ 0 (k).
(11)
Similarly, the solution of the aggregation-based MPC is
(GH )T Q (GH ) + H T RH V (k) = −(GH )T QSx(k) − (GH )T QGU¯ 0 (k) − H T RU¯ 0 (k).
Step 3. Divide W2 into W2 = (Nm−s)×(Nm−s)
W20 W21
, where W20 ∈ Rs×(Nm−s) and
W21 ∈ R . Represent the s row vectors of W20 , with the row number ij , j = 1, 2, . . . , s, as the linear combinations of
(12)
T
the rows in W21 with coefficients λj = λ1j , λ2j , . . . , λ(Nm−s)j , i.e., [λ1 λ2 · · · λs ]T W21 = W20 . Step 4. Construct matrix Λ = [−λ1 ,−λ 2 , . . . , −λs ], where λj ’s
In this subsection we will develop the quasi-equivalent aggregation strategy and embed it in Algorithm 1. We first state a lemma, whose proof is straightforward and is thus omitted, and recall the definition of the equivalent aggregation strategy from [10].
min J (k) = X T (k)QX (k) + U T (k)RU (k).
Algorithm 2. Step 1. Choose s and let W = GT QG + R. Step 2. Divide W into W = [W1 W2 ], where W1 ∈ RNm×s and W2 ∈ RNm×(Nm−s) .
I
are obtained in Step 3. Obtain H =
Λ
, where I is an identity
matrix of dimension s. Remark 1. In Algorithm 2, since W is positive definite,W21 must be of full rank. Then, since H T W = H T [W1 W2 ] = H T W1 0 according to the structure of obtained H, rank(H ) = s ≥ m, and W is of full rank, Algorithm 2 yields an equivalent aggregation matrix H for the unconstrained case by Theorem 2. The above algorithm includes the result in [10] with s = m as a special case. systems, if an aggregation strategy with H = For constrained H1 H12 ( H ∈ Rm×m ) is applied, we can find that u(k) = 1 H21 H2 H1 v1 (k)+ H12 v2 (k) and [uT (k + 1), uT (k + 2), . . . , uT (k + N − 1)]T = H21 v1 (k)+ H2 v2 (k), where V (k) = [v1T (k) v2T (k)]T . It is obvious that the current control input is determined by both v1 (k) and v2 (k). Since the constraints are about the whole control input sequence, the constraints on [uT (k + 1), uT (k + 1), . . . , uT (k + N − 1)]T are also imposed on the input u(k) by the structure of the aggregation matrix through v2 (k). In comparison with the original MPC, if the additional relationship on the first control input is removed, the restrictions on the first control input will be closer to the original MPC. Hence, we introduce the following transformation on the matrix H resulting from Algorithm 2
H0 H = H1
∈R
Nm×s
I H⇒ H = 0
0 H1
∈ RNm×(s+l) ,
(14)
where H0 ∈ Rl×s (l ≥ m) and I is an l × l identity matrix. It can be easily proven that if the resulting matrix H from the transformation (14) is of full-rank, it also satisfies the conditions in Theorem 2. Since the transformation in (14) increases the degree of freedom and guarantees the equivalent feature for the unconstrained case, we call the resulting matrix from (14) the quasi-equivalent aggregation matrix. Applying this quasi-equivalent aggregation matrix to Algorithm 1 results in the following algorithm.
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089
Algorithm 3. Off-line part: Choose s and l, use Algorithm 2 to design the equivalent aggregation matrix, and then obtain the quasiequivalent aggregation matrix H by (14). Online part: Carry out Steps 2–5 in Algorithm 1 with an H obtained in the off-line part. 4. Enlargement of the attractive region of aggregation-based MPC To enlarge the attractive region of the aggregation-based MPC, Ref. [9] suggests to relax the constraints of the optimization problem at the cost of loss of guaranteed closed-loop stability. In this section, we will develop an aggregation-based MPC algorithm based on Algorithm 3 and a switching strategy to achieve both an enlarged attractive region and guaranteed closed-loop stability. According to Theorem 1, the attractive region resulting from Algorithm 3 is the feasible region of the optimization problem in the online part with U0 (k) = 0, which can be calculated by linear difference inclusions [13] and is denoted as C0 . For C0 , its one-step set C1 is defined as a set such that for x(k) ∈ C1 , there is an admissible control input u(k) such that x(k + 1) ∈ C0 . Similarly, we can define the one-step set of Cj as Cj+1 . According to [14,13], we can compute off-line a sequence of sets {C0 , C1 , . . . , Ch } one by one from C0 . Here, h is the number of one-step reachable sets, which can be chosen according to applications and a larger h may correspond to a larger attractive region at the cost of a large amount of required memory. Since U0 (k) in Algorithm 3 can be not equal to zero, C0 may not be an invariant set for the algorithm in the online part of Algorithm 3. Hence, Cj cannot be guaranteed to belong to Cj+1 , which makes the switching design in [14] inapplicable. Here, we give the following MPC algorithm based on the aggregation matrix H obtained by the offline part of Algorithm 3 and the sets {C0 , C1 , . . . , Ch } calculated off-line. Algorithm 4. Step 1. At k = 0, U0 (0) = 0 and p(0) = 0. Step 2. When k > 0, check {x(k), x(k + N − 1|k − 1), p(k − 1)}. For case with x(k) ∈ C0 (k ≥ 0) or case with x(k + N − 1|k − 1) ∈ Xf and p(k − 1) = 1, set p(k) = 1. Otherwise, set p(k) = 0 and choose the minimum j such that x(k) ∈ Cj (h ≥ j > 0). Step 3. If p(k) = 1, calculate the optimal solution U ∗ (k) = a∗ U0 (k) + HV ∗ (k) by Step 2 of Algorithm 1 and then go to Step 5. Step 4. If p(k) = 0, solve the following optimization problem to obtain the optimal solution U ∗ (k) = a∗ U0 (k) + HV ∗ (k), min J (k),
V (k),a
(15)
s.t. u(k + i|k) ∈ U , i = 0, 1, . . . , N − 1, x(k + 1|k) ∈ Cj−1 , U (k) = aU0 (k) + HV (k). Step 5. Execute [I , 0, . . . , 0]U ∗ (k) on the controlled plant. Step 6. Calculate x(k + N |k) according to U ∗ (k), let U0 (k + 1) = [(u∗ (k + 1|k))T , (u∗ (k + 2|k))T , . . . , (u∗ (k + N − 1|k))T , (Fx∗ (k + N |k))T ]T . Set k = k + 1 and return to Step 2. Remark 2. In Algorithm 4, the constraints on system states other than x(k + 1|k) and the terminal constraints are relaxed in Step 4. Hence, the feasible region of the above algorithm, i.e., C0 ∪ C1 ∪ · · ·∪ Ch , must be larger than that of Algorithm 3. Since checking the switching conditions in Step 2 involves only algebraic calculation without iteration, the online computational burden of Algorithm 4 is much lighter than that of the original MPC. Theorem 3. Consider the aggregation-based MPC Algorithm 4. If the algorithm is feasible at time k for x(k), then the closed-loop system is asymptotically stable.
1085
Proof. At time k, Algorithm 4 is feasible. There are three cases to be considered. Case 1: When p(k) = 0, x(k) ∈ Cj (h ≥ j > 0) and optimization problem (15) is implemented. Since the algorithm is feasible, there exists an admissible control input u(k) to steer x(k) to Cj−1 . Since Cj−1 is also the one-step set of Cj−2 , if j − 2 > 0, the optimization problem (15) is feasible at the next time. This procedure will continue until x(k) ∈ C0 and the controller switches to the case with p(k) = 1. Case 2: When x(k) ∈ C0 , p(k) = 1 and Step 3 is implemented. Since the algorithm is feasible, the optimal control inputs will steer x(k + N |k) into Xf , i.e., p(k) will be kept as 1. Thus, at time k + 1, Step 3 is still implemented. Then, according to Theorem 1, the closedloop system is asymptotically stable. Case 3: If x(k + N − 1|k − 1) ∈ Xf and p(k − 1) = 1, Step 3 of Algorithm 4 is implemented. Since p(k − 1) = 1, the feasibility of Algorithm 4 at time k means that U0 (k + 1) obtained in Step 5 satisfies all the constraints. So it is a feasible solution at time k + 1 with V (k) = 0 and a = 1 and p(k) will be kept as 1. Then, according to Theorem 1, the closed-loop system is asymptotically stable. The above analysis concludes the proof. Remark 3. In the proof of Theorem 3, we see that, although C0 may not be an invariant set and Cj cannot be guaranteed to belong to Cj+1 , the condition x(k + N − 1|k − 1) ∈ Xf (when k > 0) and p(k − 1) = 1 in Step 2 of Algorithm 4 avoid the possible cyclic switching among several sets of {C0 , C1 , . . . , Ch }. This guarantees the closed-loop stability of Algorithm 4 and is the main difference from the previous design in [14]. 5. Numerical simulation This section will verify the results through two examples. The first one illustrates the effectiveness by a mathematical case. The second one considers a real world system (i.e., a two open channel network system) to verify the practical value of the proposed method. 5.1. Case 1 We first consider the following unstable system with two inputs, x( k + 1 ) =
0.8 −0.2
1 0.1 x(k) + 1.1 −0.2
0.2 u(k). 0.5
The constraints are |ui (k)| ≤ 1, i = 1, 2, and |x2 (k)| ≤ 4. Let Q and R be an identity matrix and N = 10. In order to achieve a well-designed MPC, the terminal set and QN are designed in accordance with the feedback law resulting from an LQR design. Set s = m. Then, by the off-line design of Algorithm 3 with l = 2 and l = 4, respectively, we obtain two quasi-equivalent aggregation matrices. That is, the numbers of aggregation variables are respectively 4 and 6. Shown in Fig. 1 and enclosed by the solid curve is the attractive region resulting from the MPC Algorithm 4 with h = 7 and l = 2, i.e., C7 . Meanwhile, Corg is the attractive region of the original MPC. It is obvious that C0 , which is the feasible region of Algorithm 3, is much smaller than Corg due to the reduction of freedom caused by the aggregation strategy. Comparatively, the attractive region is enlarged by Algorithm 4. For clarity, we have not plotted all the sets C0 , C1 , . . . , C7 in Fig. 1. Meanwhile, in Fig. 1, we see that C0 does not strictly belong to C2 , which verifies the discussion in Section 4. Under Algorithm 4, the state trajectories starting from several initial conditions (i.c.) near the boundary of C7 are plotted in dotted lines in Fig. 2, which shows that the closed-loop system is stable and thus verifies the conclusion of Theorem 3.
1086
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089 Table 1 Comparisons of control performance (percentage difference from that of original MPC).
Avg. cost (%) max cost (%) min cost (%) # i.c. (< 0.1%)
l=2
l=4
[9] with H1
[9] with H2
2.13 7.47 0.05 1
0.70 5.48 0 4
7.49 22.64 1.76 0
3.37 13.72 0.14 0
Table 2 Comparisons of the number of online optimization variables. Original MPC
l=2
l=4
[9] with H1
[9] with H2
20
5
7
6
8
respectively, the algorithm in [9] with I 0
Fig. 1. The attractive region.
.
H1 = .. 0 0
Fig. 2. The states and trajectories.
0 I
.. .
I 0
0 0
.. 20×6 , .∈R 0 I
I 0 0
H2 = ..
. 0 0
0 I 0
0 0 I
0 0
I 0
.. .
.. .
0 0 0
20×8 , .. ∈R . 0
I
respectively. In addition, the original MPC, which achieves the same control performance as the MPC with the infinite control horizon, is used here as a benchmark. The comparison of the accumulated cost indices is shown in Fig. 3 and Table 1. The comparison of the online optimization variables, which largely determine the online computational burden, is given in Table 2. In Fig. 3 and Table 1, results are shown for 17 randomly selected initial conditions. These conditions are marked in Fig. 2 by circle and cross (state [5, −3]T , the sixth point in Fig. 3), respectively. The result obtained from the original MPC is denoted as 100%. In Fig. 3, we can see that the proposed aggregation-based MPC with l = 4 achieves almost the same control performance as the original MPC. The control performance of the proposed aggregation-based MPC with l = 2 is also better than that of the algorithm in [9] with H1 and H2 for most initial conditions. Table 1 also reflects the same conclusion as Fig. 3. In summary, the proposed aggregation-based MPCs achieve better performance with lower online computational burden. In order to illustrate the comparison in more detail, the state responses and control inputs for initial state [5, 0]T are shown in Figs. 4 and 5. From Fig. 5, it is obvious that the control inputs resulting from the proposed MPC are similar as those resulting from the original MPC, especially for the case with l = 4, which leads to the similar state responses in Fig. 4. This owes to the similar restrictions on the control input sequence of Algorithm 4 and equivalent feature of the adopted aggregation strategy. Moreover, the degradation on control performance is reduced as the values of l increases. This comparison reinforces the comparison made in Fig. 3 and Table 1, and again validates the effectiveness of the proposed aggregation-based MPC algorithm. 5.2. Case 2
Fig. 3. Comparison of the accumulated index.
To compare the control performance, we consider two aggregation-based MPCs proposed in this paper with l = 2 and l = 4,
Consider a canal network system with two reaches in [15–17]. The structure of the canal system is shown in Fig. 6, where Hup = 2.5 m, Hdo = 1.14 m, the bed slope is 0.0003, manning roughness coefficient is taken as n = 0.015, and the third gate is fixed at U3 = 0.85 m2 . The parameters of each reach and gate are given in Tables 3 and 4, respectively. According to [16], for the canal system in Fig. 6, its mathematical model around a reference condition of uniform flow can been derived by first linearizing the Saint-Venant equations for the
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089
Fig. 4. The state responses.
Fig. 5. The control inputs.
Fig. 6. The structure of canal network system.
1087
1088
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089
Table 3 Canal parameters. Reach No.
K1i
K3i
W1i
W2i
1 2
3.76×10−5 2.85×10−5
6.38×10−5 7.51×10−5
6.5386 × 10−2 7.0538 × 10−2
−5.8505×10−2 −10.018×10−2
Table 4 Gate parameters. Gate No.
ai
1 2 3
2.41 1.91
bi
ci
−3.76 −2.99
2.99 5.08
unsteady flow of water in open-channels. With consideration of the variations of the user flow rates p, the model can be written as
v˙ (t ) = A1 v(t ) + B1 u(t ) − Ip(t )
Fig. 7. The user flow rate.
where v(t ) = [v1 (t ), v2 (t )]T is the vector of storage volume variations, u(t ) = [u1 (t ), u2 (t )]T is the vector of gate opening size variations, and p(t ) = [p1 (t ), p2 (t )]T is the vector of user flow rates. For details about the procedure to achieve (A1 , B1 ) and the physical meaning of parameters in Tables 3 and 4, refer to [15–17]. According to the continuous-time model of the canal system, by choosing the sampling time as T = 5 s, we can obtain the discretetime model as
v(k + 1) =
0.9977 0.0019
+
0.0012 v(k) 0.9966
15.1211 −3.7230
−13.9767 5 u(k) − 16.9188 0
0 p(k). 5
The input constraints are −0.3 ≤ u1 ≤ 0.3 and −0.15 ≤ u2 ≤ 0.15. Here, we choose N = 200 to approximate the MPC with infinite horizon and do not add the terminal set constraints to avoid the infeasibility of the optimization problem in the real world practice. Another reason for choosing a long horizon of MPC is that the future user flow rate can be known or scheduled in real world applications and MPC with a long horizon can make full use of the known user flow rate to improve the performance. For the original MPC, the number of online optimization variables is 400. Meanwhile, according to [16], the weighting matrices are chosen as Q = diag{1, 1.2233} and R = diag{3 × 105 , 3 × 105 }. In addition, the terminal weighted matrix is adopted to guarantee the closed-loop stability. By choosing s = 2 and l = 0, we can get an MPC controller based on Algorithm 3. Note that, since only system input variables are subject to the constraints, the controlled system is stable, and no terminal constraints are considered, the feasible region of Algorithm 3 based MPC is as large as the original one with infinite horizon. The aggregation strategy reduces the number of online optimization variables from 400 of the original MPC to 3 of Algorithm 3 based MPC. The control results and control inputs from the initial state [10, 22]T are shown in Figs. 8 and 9. Fig. 7 shows the variation of user flow rate. From these figures, we can easily see that, although the number of online optimization variables of aggregation based MPC is much smaller than the original one, the control performance is almost the same as the original MPC. This means that, by the proposed approach, we can achieve a similar control performance as the original MPC but at a much less online computation cost.
Fig. 8. The control inputs of the original MPC and aggregation based MPC.
6. Conclusions Online computational burden has always been a critical issue for the practical application of MPC algorithms. Although the aggregation strategy can reduce the online computational burden,
D. Li et al. / Systems & Control Letters 62 (2013) 1082–1089
1089
Acknowledgments The work was supported in part by the National Natural Science Foundation of China (Grant Nos. 61333009, 61104078, 61074060, 60934007, and 61221003) and the Specialized Research Fund for the Doctoral Program of Higher Education (Grant No. 20120073110017). References
Fig. 9. The control results of the original MPC and aggregation based MPC.
it may lead to loss of guaranteed closed-loop stability, poor control performance and a smaller attractive region. This paper proposed an aggregation-based MPC algorithm that addresses these issues associated with the aggregation strategy by combining the timevarying base vector, the quasi-equivalent aggregation matrix, and a switching strategy with relaxed constraints.
[1] D.Q. Mayne, J.B. Rawlings, C.V. Rao, P.O. Scokaert, Constrained model predictive control: stability and optimality, Automatica 36 (6) (2000) 789–814. [2] S.J. Qin, T.A. Badgwell, A survey of industrial model predictive control technology, Control Engineering Practice 11 (7) (2003) 733–764. [3] D. Li, Y. Xi, The general framework of aggregation strategy of model predictive control and stability analysis, in: Proceedings of the 11th IFAC Symposium on Large Scale Systems Theory and Applications, 2007, pp. 192–197. [4] N.L. Ricker, Use of quadratic programming for constrained internal model control, Industrial & Engineering Chemistry Process Design and Development 24 (4) (1985) 925–936. [5] J. Richalet, S. Abu el Ata-Doss, Ch. Arber, H.B. Kuntze, A. Jacubasch, W. Schill, Predictive functional control: application to fast and accurate robots, in: Proceedings of the 10th IFAC Congress, 1987, pp. 251–258. [6] R. Cagienard, P. Grieder, E.C. Kerrigan, M. Morari, Move blocking strategies in receding horizon control, Journal of Process Control 17 (6) (2007) 563–570. [7] G. Valencia-Palomo, M. Pelegrinis, J.A. Rossiter, R. Gondhalekar, A moveblocking strategy to improve tracking in predictive control, in: Proceeding of American Control Conference, ACC 2010, 2010, pp. 6293–6298. [8] F. Oldewurtel, R. Gondhalekar, C.N. Jones, M. Morari, Blocking parameterizations for improving the computational tractability of affine disturbance feedback MPC problems, in: Proceeding of the 48th IEEE Conference on Decision and Control, 2009, pp. 7381–7386. [9] R. Gondhalekar, J. Imura, Least-restrictive move-blocking model predictive control, Automatica 46 (7) (2010) 1234–1240. [10] D. Li, Y. Xi, Quality guaranteed aggregation based model predictive control and stability analysis, Science in China Series F: Information Sciences 52 (7) (2009) 1145–1156. [11] D. Li, Y. Xi, Aggregation based closed-loop MPC with guaranteed performance, in: The 48th IEEE Conference on Decision and Control, 2009, pp. 7400–7405. [12] D. Li, Y. Xi, The synthesis of robust model predictive control with qp formulation, International Journal of Modelling, Identification and Control 13 (1) (2011) 1–8. [13] F. Blanchini, Set invariance in control, Automatica 35 (11) (1999) 1747–1767. [14] D. Limon, T. Alamo, E.F. Camacho, Enlarging the domain of attraction of MPC controllers, Automatica 41 (4) (2005) 629–635. [15] C. Seatzu, Decentralized controllers design for open-channel hydraulic systems via eigenstructure assignment, Applied Mathematical Modelling 24 (12) (2000) 915–930. [16] G. Corriga, S. Sanna, G. Usai, Sub-optimal constant-volume control for open channel networks, Applied Mathematical Modelling 7 (4) (1983) 262–267. [17] G. Corriga, S. Sanna, G. Usai, Estimation of uncertainty in an open-channel network mathematical model, Applied Mathematical Modelling 13 (11) (1989) 651–657.
