Convergence Properties of Constrained Linear ... - Semantic Scholar
Convergence Properties of Constrained Linear System under MPC Control Law using Affine Disturbance Feedback Chen Wang, Chong-Jin Ong, Melvyn Sim January 2009
Technical Report
National University of Singapore Department of Mechanical Engineering C09-001
National University of Singapore
M.E. Dept.
Abstract This paper shows new convergence properties of constrained linear discrete time system with bounded disturbances under Model Predictive Control (MPC) law. The MPC control law is obtained using an affine disturbance feedback parametrization with an additional linear state feedback term. This parametrization has the same representative ability as some recent disturbance feedback parametrization, but its choice together with an appropriate cost function results in a different closed-loop convergence property. More exactly, the state of the closed-loop system converges to a minimal invariant set with probability one. Deterministic convergence to the same minimal invariant set is also possible if a less intuitive cost function is used. Numerical experiments are provided that validate the results.
2
National University of Singapore
1
M.E. Dept.
Introduction
This paper considers the system: xt+1 = Axt + But + wt ,
(1)
(xt , ut ) ∈ Y, wt ∈ W, ∀ t ≥ 0
(2)
where xt ∈ Rn , ut ∈ Rm and wt ∈ W ⊂ Rn are the state, control and disturbance of the system at time t respectively, and Y represents the joint constraint set on xt and ut . The study of such a system under the Model Predictive Control (MPC) framework has been an active area of research. See, for example, [3, 8, 23, 15, 18, 1, 19] and the references cited therein. An important research issue remains the choice of control parametrization within the control horizon. Several choices have been proposed in the literature [23, 8, 17, 18] and a popular choice is ut = Kxt + ct where K is a fixed feedback gain and ct is the new optimization variable. However, such a choice is known to be conservative [5, 17, 12, 24, 25] and its use will result in a relatively small domain of attraction. A natural extension of the fixed-gain parametrization is the time-varying affine state feedback ut = Kt xt + ct where both Kt and ct are the optimization variables within the horizon. However, such a parametrization is not computationally amiable as the resulting optimization problem is not convex. More recently, control parametrization based on affine function of disturbances have appeared [17, 4, 12, 24, 25]. This parametrization is appealing as the resulting problem is convex and solvable via standard numerical routines. Specifically, [17] proposes the control parametrization uL i
=
i X
Mij wi−j + vi ,
i = 0, · · · , N − 1
(3)
j=1
where Mij and vi are the optimization variables and N is the length of the horizon. [12] show that (3) is equivalent, in terms of the set of states reachable within the horizon, to that of timevarying affine state feedback. They also show that, under mild assumptions, the origin of the closed-loop system is input-to-state stable (ISS) with respect to the disturbance input under the MPC control law derived using (3) and a cost function that corresponds to the Linear Quadratic (LQ) cost for system (1) with wt ≡ 0. Recently, [24] 2008 propose an extended disturbance feedback parametrization uW i = Kf xi + ci +
N −1 X
Cij wi−j , i = 0, · · · , N − 1
(4)
j=1
where Kf is a fixed feedback gain, ci and Cij are the optimization variables. They show that parametrization (4) under the MPC framework has the same domain of attraction as that of using (3). Using an appropriate cost function, they also show a stronger stability result: state of the closed-loop system converges to the minimal disturbance invariance set, F∞ , of the system xt+1 = (A + BKf )xt + wt . Unlike (3), index j runs from 1 to N − 1 in (4) and therefore i − j can be negative in wi−j . When this happens, wi−j refers to past realized disturbances. This also means that the resulting MPC control law derived from (4) requires the values of xt and wt−1 , · · · , wt−N +1 for its evaluation at time t. 3
National University of Singapore
M.E. Dept.
This work proposes a new control parametrization based on (4) and a new cost function. The use of which results in an MPC control law requiring only the measurement of xt for its evaluation at time t. The resulting closed-loop system has the same domain of attraction as [17] and [12] but with a stronger convergence result: the closed-loop system state converges to F∞ with probability one; and deterministic convergence to the same set if a less intuitive cost function is used. The rest of this paper is organized as follows. This section ends with notations used, assumptions needed and a brief review of standard results. Section 2 states the control parametrization, the finite horizon (FH) optimization problem and the cost function used. Section 3 discusses the computation of the MPC problem. The probabilistic convergence of the state of the closed-loop system is given in section 4. Section 5 shows a formulation that strengthens the convergence result under a weaker set of assumptions. This, however, requires the use of a somewhat less intuitive cost function. Numerical examples are the contents of section 6 and they are followed by the conclusions. The following notations are used. Z+ k := {1, · · · , k} and Zk := {0, 1, · · · , k} are the respective sets of positive and non-negative integers up to k. k · k is the standard 2-norm for matrices and vectors. Given matrices A ∈ Rn×m and B ∈ Rp×q : kAkF is the Frobenious £ ¤T norm, vec(A) = AT1 · · · ATm ∈ Rnm where Ai is the ith column of A is the stacked vector of columns of A and A ⊗ B ∈ Rnp×qm is the Kronecker product of A and B. A Â (º)0 means that square matrix A is positive definite (semi-definite). For any A Â 0, kxk2A = xT Ax. 1r is a r-element column vector with all elements being 1 and In is the n × n identity matrix. For any set X, Y ⊂ Rn , X + Y := {x + y : x ∈ X, y ∈ Y } is the Minkowski sum of X and Y . The system (1)-(2) is assumed to satisfy the following assumptions: (A1) system (A, B) is stabilizable; (A2) the set Y is a polytope having a characterization Y := {(x, u)| Yx x + Yu u ≤ 1q } ⊂ Rn+m
(5)
for some Yx ∈ Rq×n and Yu ∈ Rq×m ; (A3) the disturbance wt , t ≥ 0 are independent and identically distributed (i.i.d.) with zero mean and W is a polytope characterized by W := {w| Hw ≤ 1r } ⊂ Rn ,
(6)
for some H ∈ Rr×n ; (A4) a compact and non-empty constraint-admissible disturbance invariant set exists for system (1)-(2) under the feedback law u = Kf x and takes the form Xf := {x| Gx ≤ 1g } ⊂ Rn
(7)
for some G ∈ Rg×n and that Xf contains the minimal disturbance invariant set of xt+1 = (A + BKf )xt + wt in its interior. 4
National University of Singapore
M.E. Dept.
The above assumptions can be rationalized in the following ways. (A1) is standard. The characterization of Y in (A2) is made out of the need for a concrete computational representation. (A3) is a typical assumption on the disturbance and has been used in several prior works [11, 22]. That W is a polytope and contains the origin in its interior is an assumption made out of convenience for the presentation of the main result. Relaxation of (A3) is possible and the details are given in section 5 and illustrated via a numerical example in section 6. The existence of Xf in (A4) has been shown by [13] 1998 provided that W is sufficiently small. More exactly, for any feedback gain Kf ∈ Rm×n such that Φ := A + BKf is strictly stable and sufficiently small W , Xf is the maximal, constraint admissible and disturbance invariant in the sense that Φx + w ∈ Xf , (x, Kf x) ∈ Y for all x ∈ Xf and for all w ∈ W . It is also known [14, 9] that the state of the system xt+1 = Φxt + wt converges to the minimal disturbance invariant set, F∞ , given by F∞ (Kf ) := W + (A + BKf )W + (A + BKf )2 W + · · ·
(8)
and that F∞ is compact. The assumption that F∞ (Kf ) ⊂ Xf is also a natural consequence when W is not too large.
2
Control parametrization
MPC formulation solves an N -stage finite horizon (FH) optimization problem. Let xi and ui , i ∈ ZN −1 denote the predicted state and predicted control at the ith stage respectively within the horizon. The proposed control parametrization within the FH optimization problem takes the form i X ui = Kf xi + di + Dij wi−j for all i ∈ ZN −1 (9) j=1
where di ∈ Rm , Dij ∈ Rm×n , j ∈ Z+ i , i ∈ ZN −1 are the variables of the FH problem and Kf is the given feedback gain in (A4). Since i − j ≥ 0, wi−j is the (i − j)th predicted disturbance at each stage i. In this regard, (9) is similar to (3) in that only predicted disturbances are used W in the parametrization. In addition, ui is equivalent to uL i and ui in terms of the family of functions that they represent. To state this precisely, let 0 ··· 0 0 D1 ··· 0 0 £ T T ¤T 1 T Nm N m×N n d := d0 d1 · · · dN −1 ∈ R , D := (10) .. .. .. ∈ R .. . . . . N −1 DN −1 · · ·
1 DN −1 0
and v and M to be similarly defined as d and D in structure but having entries vi instead W of di and Mij instead of Dij respectively. The equivalence of ui , uL i and ui are stated in the following theorem. N −1 Theorem 1. Suppose x0 and Kf are given. Then (i) the set {ui }i=0 of (9) is equivalent N −1 L to {ui }i=0 of (3) in the sense that for any (v, M), there exists a unique (d, D) such that −1 N −1 L N −1 W N −1 {ui }N i=0 = {ui }i=0 and vice versa. (ii) {ui }i=0 is equivalent to {ui }i=0 of (4).
5
National University of Singapore
M.E. Dept.
Proof. See Appendix A. The FH optimization problem under parametrization (9), referred hereafter as PN (x), is min J(d, D)
(11)
s.t. x0 = x,
(12)
d,D
xi+1 = Axi + Bui + wi , i ∈ ZN −1 ui = Kf xi + di +
i X
Dij wi−j , i ∈ ZN −1
(13) (14)
j=1
(xi , ui ) ∈ Y, ∀wi ∈ W, i ∈ ZN −1
(15)
x N ∈ Xf ,
(16)
∀wi ∈ W, i ∈ ZN −1
where Y and Xf are the corresponding sets given by (5) and (7) respectively. The cost function J(d, D) takes the form N −1 i X X kdi k2Ψ + J(d, D) := (17) kvec(Dj )k2Λ i
i=0
j=1
for any choice of Ψ ∈ Rm×m and Λ ∈ Rmn×mn that satisfy Ψ = ΨT Â 0, Λ º Σw ⊗ Ψ
(18)
where Σw is the covariance matrix of wt and vec(·) is stacking operator defined in Section 1. Clearly, J(d, D) is a measure of the deviation of ui of (14) from the linear control law Kf xi and the motivation for it as the objective function is clear: penalizing the use of nonzero (d, D) forces the asymptotic behavior of the closed-loop system to approach that of xt+1 = (A + BKf )xt . The technical condition (18) is to ensure convergence of the closed-loop states and its exact role will become clear in the proof of Theorem 3. Several comments on J(d, D) are in order. Remark 1. A connection between J(d, D) and the standard LQ cost can be established. Specifically, suppose Q º 0 and R Â 0 are given and let P = AT P A−AT P B(R+B T P B)−1 B T P A+Q, the solution of the algebraic Riccati equation, Ψ = R + B T P B, Λ = Σw ⊗ Ψ and Kf = −(R + B T P B)−1 B T P A. It is shown [10, 24, 25] that "N −1 # X Ew (kxi k2Q + kui k2R ) + kxN k2P = xT0 P x0 + N trace(Σw P ) + J(d, D) (19) i=0
where Ew is the expectation taken over {w0 , · · · , wN −1 } within the horizon. Since the first two terms on the right hand side of (19) are independent of (d, D), minimizing J(d, D) is equivalent to minimizing the expected infinite horizon LQ cost over the disturbance input. Remark 2. From (18) and (19), it may appear that Σw is needed for the determination of Λ. This is not true. The choice of Λ can be made to satisfy (18) even when Σw is not known accurately. One simple choice is to let Λ = α2 In ⊗ Ψ where α := maxw∈W kwk. Then it follows that Λ º Σw ⊗ Ψ because α2 In º wwT for all w ∈ W and α2 In ⊗ Ψ º E[wwT ] ⊗ Ψ. Consequently, (A3) provides for conditions that guarantee the computability of maxw∈W kwk. 6
National University of Singapore
M.E. Dept.
Further discussion on the choice of Ψ and Λ and their influence on closed-loop system trajectories are discussed in section 4. Several associated sets, needed to facilitate the discussions in the sequel, are introduced. Let the feasible set of optimization problem PN (x) be TN := {(x, d, D)|(d, D) is feasible to PN (x)}
(20)
and the set of admissible initial states, or equivalently, the domain of attraction of the MPC controller is XN := {x|∃ (d, D) such that (x, d, D) ∈ TN }. (21) The rest of the MPC formulation is standard: PN (xt ) is solved at each time t to obtain the optimizer (d∗t , D∗t ) := (d∗ (xt ), D∗ (xt )) and the first control, u∗0 , is applied to (1) at time t resulting in the MPC control law, ut = u∗0 = Kf xt + d∗0 .
(22)
Remark 3. From (20) and (21), it is easy to see that TN and XN depend only on the conW straints (12)-(16). If uL i or ui replaces ui in (14), it follows from Theorem 1 that the corresponding domain of attraction, XNL and XNW is the same as XN . However, the stability of the corresponding closed-loop systems can differ. See Remark 6.
3
The Computation of PN (x)
Following (5), (7), (12)-(14), inequalities (15) and (16) can be collectively restated as £ ¤ ¯ + Bd ¯ + max BD ¯ + G¯ w ≤ 1s Ax w∈WN
(23)
£ ¤T T N ¯ ¯ where s = N q + g, w := w0T w1T · · · wN −1 , W := {w|Hw ≤ 1` } with H = IN ⊗ H, ` = ¯ B¯ and G¯ are appropriate matrices given in Appendix B and the max N r following (6), A, operator is meant to be taken row-wise. Correspondingly, the ith row of (23) can be rewritten ¯ ≤ 1` } ≤ bi for some ei ∈ RN m , bi ∈ R that depend on x, d and D. Let as max{eTi w|Hw ¯ Then, it follows by duality zi ∈ R` be the Lagrange multiplier corresponding to the rows of H. T T T ¯ ¯ that max{ei w|Hw ≤ 1` } = min{zi 1` |H zi = ei , zi ≥ 0}. Collecting over all the rows of (23), the PN (x) can be equivalently stated as min J(d, D)
d,D,Z
(24)
¯ + Bd ¯ + ZT 1` ≤ 1s s.t. Ax ¯ = BD ¯ + G¯ ZT H
(25)
zi ≥ 0, i = 1, . . . , s
(27)
(26)
where Z = [z1 · · · zs ] ∈ R`×s and the minimization of Z is relaxed since the existence of any one feasible Z is enough to guarantee that (x, d, D) ∈ TN . Remark 4. The above duality results can be extended to W sets that are non-polyhedral. See, for example, treatments of such sets in [6, 20]. If W is a second-order cone [16, 2] representable 7
National University of Singapore
M.E. Dept.
bounded set with non-empty interior such that WN = {w| kLi w − li k ≤ λTi w − θi , i ∈ Z+ k} T w| w ∈ for some matrices Li , li , λi and θi , i ∈ Z+ , then it follows from duality that max{e P Pk k T + WN } = min(µi ,ηi ) { ki=1 (µTi li −ηi θi )| i=1 (Li µi −ηi λi ) = e, kµi k ≤ ηi , i ∈ Zk }. Similarly, if W is a bounded set with non-empty interior such that WN = PN nsemi-definite cone representable + N n {Ω ∈ R | i=1 Ωi Ci − F < 0, i ∈ ZN n } where Ci and F are symmetrical matrices of appropriate dimension, then max{eT w| w ∈ WN } = minΥ {Trace(F Υ)| Trace(Ci Υ) = ei , i ∈ Z+ N n , Υ 4 0}. Remark 5. While the duality result is available for W being a second-order or semi-definite cone representable set, the availability of Xf satisfying (A4) deserves some clarifications. When W is non-polyhedral, computation of a constraint-admissible disturbance invariant set Xf may not be easy. A simple approach is to construct a polytope Wp such that Wp ⊃ W and Wp ≈ W . In that case, a Xf satisfying (A4) can be constructed using Wp following existing computational methods [14]. Using this Xf in (16) and Remark 4, PN (x) becomes either a second-order cone or a semi-definite cone programming problem. It is worthy to note that the use of such an Xf in (16) and with wi ∈ W for all i ∈ ZN −1 in both (15) and (16) is less conservative than replacing W by Wp throughout (12)-(16). An example using such an approach is illustrated in Section 6.
4
Feasibility and Probabilistic Convergence
The feasibility of PN (xt ) at different time instants and stability of the closed-loop system under the feedback law (22) are addressed in this section. Theorem 2. Suppose (A1)-(A4) are satisfied, the FH optimization problem PN (x) has the following properties (i) TN is convex and compact. (ii) If x ∈ XN , the optimal solution of PN (x) exists. (iii) If PN (xt ) admits an optimal solution, so does PN (xt+1 ) under the feedback law (22) for all wt ∈ W . Proof. See Appendix C. The main result of probabilistic convergence of the state of the closed-loop system is stated in the next theorem. Theorem 3. Suppose x0 ∈ XN and (A1)-(A4) are satisfied. System (1) under MPC control law (22) obtained from the solution of PN (x) using cost function (17) with condition (18) satisfied has the following properties: (i) (xt , ut ) ∈ Y for all t ≥ 0, (ii) xt → F∞ (Kf ) with probability one as t → ∞ (iii) xt enters Xf in finite time with probability one. Proof. See Appendix D. Remark 6. It is of interest to know if the results of Theorem 3 can be extended to the case where uL i is used in (14) in view of Theorem 1. As seen in the proof of Theorem 3, the convergence property depends on the choices of the cost function and the control parametrization. Since (9) and (3) are equivalent when condition (36) is satisfied, the results of Theorem 3 is applicable L under the following conditions: uL i is used in (14); a new cost function J (M, v) :=J((I + 8
National University of Singapore
M.E. Dept.
KϕB)−1 (v − KϕAx), (I + KϕB)−1 (M − KϕG)) is used in (11) where A, B, K, G and ϕ are those given in (32) and (33); and (M, v) becomes the optimization variables for PN (x). The approach by [12] uses the nominal LQ cost as the cost function and it is not clear if the results of Theorem 3 remains true under that situation. One associated issue in the formulation of PN (x) is the choices of Ψ and Λ in J(d, D). How should Ψ and Λ be chosen and how do these choices affect the closed-loop system trajectories? As xt → F∞ with probability one from result (ii) of Theorem 3, it implies that xt enters Xf with probability one and stays within thereafter since F∞ ⊂ Xf . When this happens, the optimal (d, D) are zero in PN (x) and the MPC control law becomes ut = Kf xt for all t thereafter. The closed-loop system behavior then corresponds to that of the system xt+1 = (A + BKf )xt . Clearly, the choices of Λ and Ψ does not affect the asymptotic behavior of the system but only the transient when xt ∈ / Xf . Suppose Λ = Σw ⊗ Ψ. Then admissible changes in Ψ will not result in changes in the system behavior since Σw ⊗ Ψ is linear in Ψ. On the other hand, if Ψ is fixed, Λ can be chosen to be increasingly ”larger” than Σw ⊗ Ψ. In loose terms, a ”larger” choice of Λ penalizes the use of D versus the use of d in J(d, D). Such a preference would mitigate the effect of the disturbance feedback component in the control parametrization, resulting in a parametrization that is closer in spirit to ut = Kf xt +dt of [8]. When this happens, the transient response for the system may become slower even though the domain of attraction XN remains unaffected. This observation together with the associate details used in the experimental study are discussed in section 6.
5
Deterministic Convergence
While the assumption of W being a compact set is reasonable, the assumption of wt being zero mean and i.i.d. is harder to verify in practice. This section is concerned with the relaxation of assumption (A3) while achieving a stronger convergence result than that of Theorem 3. Consider (A3a) wt ∈ W and W is a polytope characterized by W := {w| Hw ≤ 1r } ⊂ Rn for some H ∈ Rr×n . and define the cost function V (d, D) :=
N −1 X
i X kdi k2Ψ + (γ1 kvec(Dij )k2 + γ2 kvec(Dij )k)
i=0
(28)
j=1
for some constants γ1 and γ2 satisfying γ1 ≥ α2 kΨk, γ2 ≥ 2αβkΨk
(29)
where α := maxw∈W kwk and β := max(x,d,D)∈TN ,i∈ZN −1 kdi k. The existence of α and β are guaranteed by compactness of the W and TN sets, provided for in (A3a) and part (i) of Theorem 2 respectively. 9
National University of Singapore
M.E. Dept.
Theorem 4. Suppose x0 ∈ XN and (A1-A2),(A3a) and (A4) are satisfied and J(d, D) is replaced by V (d, D) in PN (x) satisfying condition (29), then system (1) under the MPC control law (22) satisfies (i) (xt , ut ) ∈ Y for all t ≥ 0, (ii) xt → F∞ (Kf ) as t → ∞ (iii) xt enters Xf in finite time. Proof. See Appendix E Remark 7. Several choices of the cost function of (28) are possible. For example, the results of Theorem 4 remain true if kvec(Dij )k is replaced by kDij k. This may be more appealing as less conservative bounds on γ1 and γ2 can be found to ensure the non-negativity of p(wt ). However, its use will result in a semi-definite programming problem for PN (x) and is less desirable computationally. The use of kvec(Dij )k or kDij kF results in a second-order cone programming for PN (x) and is computationally more amiable. Remark 8. The computation of β can be simplified to β = max(x,d,D)∈TN kd0 k, see Appendix F for details. Note that any upper bound of β can be used to guarantee the results of Theorem 4. One such upper bound is β¯ := kσk where σi := max(x,d,D)∈TN |d0 (i)| and d0 (i) is the ith element of d0 .
6
Numerical Examples
Four experiments are conducted on a system to validate the results of the previous sections. The parameters and constraints of the system are: A = [1.1 1; 0 1.3], B = [1 1]T , Kf = [−0.7434 − 1.0922], Y = {(x, u)| |u| ≤ 1, kxk∞ ≤ 8} and W = {H w| ˜ kwk ˜ ∞ ≤ 0.2} where H = [1 − 0.2; 0 1] and w ˜ ∈ R2 is a random vector uniformly distributed over [−0.2, 0.2] × [−0.2, 0.2] with covariance matrix Σw˜ = 0.0133I2 . Expressed in the form of (5), the set Y is Yx = [0 0; 0 0; 0 1/8; 0 − 1/8; 1/8 0; − 1/8 0] and Yu = [1; −1; 0; 0; 0; 0]. The set Xf is the corresponding maximal constraint-admissible disturbance invariant set of (1) under ut = Kf xt given by Xf = {x| Gx ≤ 14 }, where G = [−0.7434 −1.0922; 0.7434 1.0922; 0.8252 − 0.2391; −0.8252 0.2391]. The first experiment, Experiment I, uses the cost function of (17) with Ψ = 1, Λ = Λop := Σw ⊗ Ψ = [0.0139 − 0.0027; −0.0027 0.0133], N = 8 and x0 = [−4 2]T . The simulation results over 15 different disturbance realizations shown in Fig. 1 to 4 by solid lines. It is clear from Fig. 1 and 2 that the constraints are satisfied by all trajectories, in accordance to property (i) of Theorem 3. In addition, Figure 1 shows the convergence of xt into Fˆ∞ (Kf ), a tight outer bound of F∞ (Kf ) obtained using procedures given in [21]. This convergence is further verified in Fig. 3 and 4 where the plots of dis(xt , Fˆ∞ ) := minx∈Fˆ∞ kx − xt k, the minimum distance to Fˆ∞ , and dt := d∗0|t against increasing t are shown respectively. The case where W is non-polyhedral is shown in Experiment II, in connection to Remarks 4 and 5. A different disturbance characteristic is used here: w is uniformly distributed over ¯ := {w| W ¯ kS1 wk ¯ ≤ 1, kS2 wk ¯ ≤ 1} where S1 = [5 1; 0 2.5] and S1 = [2.5 0.5; 0 5]. Note ¯ , is needed for the computation of Xf that a tight bounding polytope, Wp , such that Wp ⊃ W satisfying (A4) and it corresponds to the W set of Experiment I (see Fig. 5). Also, Ψ and Λ of the first experiment are used and it is easy to verify that condition (18) remains true 10
National University of Singapore
M.E. Dept.
¯ is a second-order cone representable set and the conversion because Σw  Σw¯ . In this case, W ¯ follows the expression in Remarks 4, resulting in PN (x) being a of (15) and (16) for all wi ∈ W second-order cone programming problem. The simulation results with N = 8 and x0 = [2 −1]T for 15 different realizations of {wt } are plotted in Fig. 1 to 4 using dash-dot lines. Experiment III is designed to understand the influence of Λ and Ψ of (17) on the performance of the closed-loop system. As stated in section 4, choices of these matrices affect only the transient behavior when xt ∈ / Xf and not the asymptotic behavior of the closed-loop system. To quantify the transient, the average number of time step, tf (x0 ), taken to enter Xf from a given x0 is reported. Here, the average is taken over different realizations of the disturbances. Without loss of generality, values of Λ is increased from Λop (see discussion in section 4). Table 1 shows the values of tf (x0 ) and the associated standard deviations over 20 disturbance realizations for several choices of x0 , N and Λ. For each x0 , the same 20 disturbance realizations are used for the different Λ in computing tf (x0 ) and the standard deviations. From the table, tf (x0 ) generally increases when Λ increases. For comparison purpose, the corresponding trajectories of the system under same settings as the Experiment I except for Λ = 104 Λop are plotted in Fig. 1 to 4 using dash lines. From Fig. 3 and 4, the slower convergence of the state and control trajectories are clearly evident. 3
x(2)
2
Xf
1
0
Fˆ∞
−1
−2 −4
−3
−2
−1
0
1
2
x(1)
Figure 1: State trajectories of the first three experiments: solid line for first experiment, dashdot line for the second and dash line for the third. The last experiment, Experiment IV, considers the case discussed in Section 5. The system parameters are the same as those in the first experiment except that the distribution of w˜ is assumed to be unknown. The parameters of (29) are: α = 0.3124 and β = 2.7307 (when N = 8) and β = 3.5425 (when N = 10). Correspondingly, the weight matrices of (28) are Ψ = 1, γ1op := α2 kΨk = 0.0976, γ2op := 2αβkΨk = 1.7059 (2.213 when N = 10). The values of γ1 and γ2 are increased separately and jointly to assess their influence on the system behavior. The general effect of increasing values of γ1 and γ2 appears to have similar trend on the system as the increase in Λ. The time taken to reach Xf from any given x0 increases, although to a lesser percentage than that by Λ, with increasing values of γ1 and γ2 with γ2 having a heavier influence.
11
National University of Singapore
M.E. Dept.
1
ut
0.5
0
−0.5
−1 1
2
3
4
5
6
7
8
9
10
t
Figure 2: Control trajectories of the first three experiments: solid line for first experiment, dashdot line for the second and dash line for the third.
4
dis(xt , Fˆ∞ )
3
2
1
0
−1
0
2
4
6
8
10
t
Figure 3: Distance between states and F∞ (Kf ) of the first three experiments: solid line for first experiment, dashdot line for the second and dash line for the third.
12
National University of Singapore
M.E. Dept.
0.5
dt = d∗0|t
0
−0.5
−1
−1.5
−2
1
2
3
4
5
6
7
8
9
10
t
Figure 4: Values of dt of the first three experiments: solid line for first experiment, dashdot line for the second and dash line for the third.
Wp = W
0.2 0.15
w(2)
0.1 0.05 0 −0.05 −0.1 −0.15 −0.2 −0.25
¯ W −0.2
−0.1
0
0.1
w(1)
¯ set. Figure 5: Wp set and W
13
0.2
0.3
N 8 8 10 10
x0
[−4 2]T
[−2.5 − 1.2]T
[−4 − 1]T
[−6 2]T
Initial Condition 102 Λop 4.9 (0.4472) 5.15 (0.4894) 6.75 (0.6387) 7.15 (0.6708)
Λop 4.2 (0.4104) 4.9 (0.5525) 6.3 (0.7327) 6.45 (0.6863)
Λ
8.55 (0.6863)
7.45 (0.8256)
5.15 (0.4894)
5.65 (0.4894)
104 Λop
7.9 (0.4472)
7.2 (0.6959)
5.15 (0.4894)
5.2 (0.4104)
(γ1op , γ2op )
8.05 (0.6863)
7.25 (0.7164)
5.15 (0.4894)
5.4 (0.5026)
(10γ1op , γ2op )
8.6 (0.5982)
7.45 (0.8256)
5.15 (0.4894)
5.65 (0.4894)
(γ1op , 10γ2op )
(γ1 , γ2 )
Table 1: Average time step, tf (x0 ), and its standard deviation
8.6 (0.5982)
7.45 (0.8256)
5.15 (0.4894)
5.65 (0.4894)
(10γ1op , 10γ2op )
National University of Singapore M.E. Dept.
14
National University of Singapore
7
M.E. Dept.
Conclusions
Convergence results for constrained linear system under MPC control law using a new control parametrization and a new cost function are presented. The parametrization uses affine disturbance feedback together with a linear state feedback term, Kf x, and is a modification of the parametrization by [25]. Such a parametrization is similar to existing disturbance parameterizations in the literature in terms of the set of state reachable. Using the proposed cost function and the new parametrization, the closed-loop system state converges to the minimal robust invariant set F∞ (Kf ) with probability one. Deterministic convergence to F∞ (Kf ) is also possible using a less intuitive cost function. The asymptotic behavior of the closed-loop system is determined by the choice of Kf so long as the weight matrices of the cost function satisfy some mild conditions.
References [1] Alamo, T., Mu˜ noz de la Pe˜ na, D., Limon, D. and Camacho, E. F. [2005]. Constrained min-max predictive control: modifications of the objective function leading to polynomial complexity, IEEE Transaction on Automatic Control 50(5): 710–714. [2] Alizadeh, F. and Goldfarb, D. [2003]. Second-order cone programming, Methimetical Programming 95(1): 3–51. [3] Bemporad, A. [1998]. Reducing conservativeness in predictive control of constrained systems with disturbances, Proceedings of 37th Conference on Decision and Control, Tampa, Florida, pp. 1384–1389. [4] Ben-Tal, A., Goryashko, A., Guslitzer, E. and Nemirovski, A. [2004]. Ajustable robust solutions of uncertaint linear programs, Mathematical Programming 99(2): 351–376. [5] Ben-Tal, A. and Nemirovski, A. [1999]. Robust solutions to uncertain linear programs, OR Letters 25. [6] Ben-Tal, A. and Nemirovski, A. [2001]. Lectures on Modern Convex Optimization: Analysis, Algorithms, Engineering Applications, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA. [7] Bertsekas, D. P. [2003]. Nonlinear Programming, 2nd edn, Athena Scientific. [8] Chisci, L., Rossiter, J. A. and Zappa, G. [2001]. Systems with persistent disturbances: predictive control with restricted constraints, Automatica 37(7): 1019–1028. [9] Gilbert, E. and Ong, C. J. [2008]. Linear systems with hard constraints and variable set points: Their robustly invariant sets, Technical report, Department of Aerospace Engineering, University of Michigan and Department of Mechanical Engineering, National University of Singapore, available at http://guppy.mpe.nus.edu.sg/ mpeongcj/ongcj.html. C08-002. 15
National University of Singapore
M.E. Dept.
[10] Goulart, P. J. and Kerrigan, E. C. [2005]. On the stability of a class of robust receding horizon control laws for constrained systems, Technical report, Department of Engineering, University of Cambridge. CUED/F-INFENG/TR.532. [11] Goulart, P. J. and Kerrigan, E. C. [2006]. Robust receding horizon control with an expected value cost, Proc. UKACC International Conference (Control 2006), Glasgow, Scotland. [12] Goulart, P. J., Kerrigan, E. C. and Maciejowski, J. M. [2006]. Optimization over state feedback policies for robust control with constraints, Automatica 42(4): 523–533. [13] Kolmanovsky, I. and Gilbert, E. G. [1995]. Maximal output admissible sets for discretetime systems with disturbance inputs, Proceedings of the 1995 American Control Conference, Seattle, pp. 1995–1999. [14] Kolmanovsky, I. and Gilbert, E. G. [1998]. Theory and computation of disturbance invariant sets for discrete-time linear systems, Mathematical Problems in Engineering 4(4): 317– 367. [15] Lee, Y. I. and Kouvaritakis, B. [1999]. Constrained receding horizon predictive control for systems with disturbances, International Journal of Control 72(11): 1027–1032. [16] Lobo, M. S., Vandenberghe, L., Boyd, S. and Lebret, H. [1998]. Applications of secondorder cone programming, Linear Algebra and its Applications 284: 193–228. [17] L¨ofberg, J. [2003]. Approximations of closed-loop minimax mpc, Proceedings of the 42nd IEEE Conference on Decision and Control, Maui, Hawaii, USA, pp. 1438–1442. [18] Mayne, D. Q., Seron, M. M. and Rakovi´c, S. V. [2005]. Robust model predictive control of constrained linear systems with bounded disturbances, Automatica 41(2): 219–224. [19] Mu˜ noz de la Pe˜ na, D., Alamo, T., Bemporad, A. and Camacho, E. F. [2006]. A decomposition algorigthm for feedback min-max model predictive control, IEEE Transaction on Automatic Control 51(10): 1688–1692. [20] Nemirovski, A. [2006]. Advances in convex optimization: conic programming, Proceedings of the International Congress of Mathematicians, Madrid, pp. 413–444. [21] Ong, C. J. and Gilbert, E. G. [2006]. The minimal disturbance invariant set: Outer approximations via its partial sums, Automatica 42(9): 1563–1568. [22] Primbs, J. A. [2007]. A soft constraint approach to stochastic receding horizon control, Proceedings of the 46th IEEE Conference on Decision and Control, New Orleans, Louisiana, USA, pp. 4797 – 4802. [23] Rossiter, J. A., Kouvaritakis, B. and Rice, M. J. [1998]. A numerically robust state-space approach to stable predictive control strategies, Automatica 34(1): 65–73. [24] Wang, C., Ong, C. J. and Sim, M. [2007]. Model predictive control using affine disturbance feedback for constrained linear system, Proceedings of the 46th IEEE Conference on Decision and Control, New Orleans, Louisiana, USA, pp. 1275–1280. 16
National University of Singapore
M.E. Dept.
[25] Wang, C., Ong, C. J. and Sim, M. [2008]. Constrained linear system with disturbance: stability under disturbance feedback, To appear in Automatica . [26] Williams, D. [1991]. Probability with Martingales, Cambridge University Press.
A
Proof of Theorem 1
Proof. (i) Suppose the predicted states, predicted controls and predicted disturbances within £ ¤T £ ¤T the horizon are x := xT0 xT1 · · · xTN ∈ R(N +1)n , u := uT0 uT1 · · · uTN −1 ∈ RN m , and £ ¤T N −1 T w := w0T w1T · · · wN ∈ RN n . The state x and the and the control sequence, {ui }i=0 , −1 defined by (9) can be equivalently stated as
where A :=
In A A2 .. .
, B :=
AN
0 B AB .. .
x = Ax0 + Bu + Gw
(30)
u = Kx + d + Dw
(31)
0 0 B .. .
··· ··· ··· .. .
0 0 0 .. .
AN −1 B AN −2 B · · ·
B
, G :=
0 I A .. .
0 0 I .. .
··· ··· ··· .. .
0 0 0 .. .
AN −1 AN −2 · · ·
I
, K = [IN ⊗ Kf 0] (32)
and d and D are those given by (10). Using (31) in (30), we get x = ϕAx0 +ϕBd+(ϕBD+ϕG)w where ϕ = (I − BK)−1
(33)
u = KϕAx0 + (I + KϕB)d + [(I + KϕB)D + KϕG] w
(34)
and u becomes N −1 Consider the parametrization (3) and the sequence {uL i }i=0 expressed by the variables (M, v) of (10). It follows that uL = v + Mw (35)
Comparing (34) and (35), u = uL if and only if ( KϕAx0 + (I + KϕB)d = v (I + KϕB)D + KϕG = M.
(36)
Note that (I + KϕB) is a lower triangular matrix with all diagonal elements being 1 and is always invertible. In addition, KϕG is a strict lower triangular block matrix like M and its multiplication by (I + KϕB)−1 on the left also results in a strict lower triangular block matrix. Hence, the mapping between (M, v) and (D, d) by (36) is unique or one-to-one for all choices −1 L N −1 of K and x0 . This establishes the equivalence of {ui }N i=0 and {ui }i=0 . 17
National University of Singapore
M.E. Dept.
(ii) Set Cij = 0 for all j > i in (4) and it follows that ui is a special case of uW i . To show the converse, let ( P −1 j di = ci + N j=i+1 Ci wi−j , i ∈ ZN −1 (37) Dij = Cij j ≤ i, i ∈ ZN −1 N −1 W N −1 for any ci , Cij that defines uW i . This establishes the equivalence of {ui }i=0 and {ui }i=0 .
B
Expressions of Matrices in (23)
¯ B¯ and G¯ in (23) are The A, · A¯ = Y
ϕA KϕA
¸
· , B¯ = Y
ϕB I + KϕB
¸
· , G¯ = Y
ϕG KϕG
¸
· , Y=
IN ⊗ Yx 0 IN ⊗ Yu 0 G 0
¸
where A, B, G, K and ϕ are defined in (32) and (33).
C
Proof of Theorem 2
Proof. (i) Since Y is compact from (A2), the projection of Y onto x and u space, denoted by XY and UY respectively, are bounded. From (9) and the fact that W is a polytope in Rn , Dij must be bounded in order for ui ∈ UY . Since the origin is inside W , Kf xi + di must be inside UY . Therefore, di is bounded as xi and UY are bounded. This, together with TN being closed and convex from (25)-(27) leads to the desired result.(ii) Since x ∈ XN , PN (x) is feasible. From (i), this means that ΠN (x) := {(d, D)|(x, d, D) ∈ TN } is compact. This, together with the fact that J(d, D) is continuous with respect to (d, D) means that the optimal solution exists by Weierstrass’ Theorem [7]. (iii) The proof follows standard arguments but the details are given for their relevance to Theorem 3. For clarity, additional subscripts “|t” and “|t + 1” are used to denote the variables at the different times. Let (d∗t , D∗t ) denote the optimal solution ˆ t+1 , D ˆ t+1 ) is chosen as of PN (xt ). At time t + 1, wt is realized and (d dˆi|t+1 = ˆj D i|t+1 =
( i+1 ∗ d∗i+1|t + (Di+1|t ) wt 0 ( j (Di+1|t )∗ 0
i ∈ ZN −2 i=N −1
+ j ∈ Z+ i , i ∈ ZN −2
j ∈ Z+ N −1 , i = N − 1.
(38) (39)
ˆ t+1 , D ˆ t+1 ) is feasible to PN (xt+1 ) for all possible wt ∈ W due to the disturbance This choice of (d invariance of Xf for system (1) under control law ut = Kf xt and that (d∗t , D∗t ) is the optimal solution at time t. That the the optimum of PN (xt+1 ) exists follows from the compactness of ΠN (xt+1 ) and the Weierstrass’ theorem [7]. 18
National University of Singapore
D
M.E. Dept.
Proof of Theorem 3
Proof. (i) The stated result follows directly from Theorem 2. (ii) Let Jt∗ := J(d∗t , D∗t ) and ˆ t+1 (wt ), D ˆ t+1 (wt ), D ˆ t+1 ) where (d ˆ t+1 ) are given by (38)-(39). Then it follows Jˆt+1 (wt ) := J(d that Jt∗ − Jˆt+1 (wt ) =
N −1 X
(kd∗i|t k2Ψ
− kdˆi|t+1 k2Ψ ) +
i=0
N −1 X
i ∗ 2 kvec(Di|t ) kΛ
i=1
=kd∗0|t k2Ψ + =kd∗0|t k2Ψ
+
N −1 X
N −1 X
i=1 N −1 X
i=1
(kd∗i|t k2Ψ − kdˆi−1|t+1 k2Ψ ) + (kd∗i|t k2Ψ
−
kd∗i|t
+
i ∗ (Di|t ) wt k2Ψ )
i=1
=kd∗0|t k2Ψ where g(wt ) =
i ∗ 2 kvec(Di|t ) kΛ
+
N −1 X
i ∗ 2 kvec(Di|t ) kΛ
i=1
+ g(wt ) N −1 X
(40)
i ∗ 2 i ∗ i ∗ (kvec(Di|t ) kΛ − 2(d∗i|t )T Ψ(Di|t ) wt − k(Di|t ) wt k2Ψ ).
(41)
i=1
Taking the expectation of (40) over wt , it follows that h i Jt∗ − kd∗0|t k2Ψ = Ewt Jˆt+1 (wt ) + Ewt [g(wt )] h i ≥ Ewt Jˆt+1 (wt ) £ ∗ ¤ £ ∗ ¤ (wt ) . (wt ) = Et Jt+1 ≥ Ewt Jt+1
(42) (43)
where Et in (43) is the expectation taken over wi , i ≥ t. Inequality (42) follows from the fact that Ewt [g(wt )] ≥ 0. This is so because by taking the expectation of (41), one gets Ewt [g(wt )] =
N −1 X
i ∗ 2 i ∗ 2 i ∗ (kvec(Di|t ) kΛ − kvec(Di|t ) kΣw ⊗Ψ − 2(d∗i|t )T Ψ(Di|t ) E[wt ])
(44)
i=1
where the last term is zero due to (A3) and the rest is non-negative due to (18). The inequality ∗ (w ) for every w ∈ W which implies that in (43) follows from the fact that Jˆt+1 (wt ) ≥ Jt+1 t t ∗ ∗ (w ) depends on ˆ Ewt [Jt+1 (wt )] ≥ Ewt [Jt+1 (wt )]. Equality (43) follows from the fact that Jt+1 t wt only and not on any wi , i > t. Repeating the inequality of (43) for increasing t, one gets, £ ∗ ¤ ∗ Jt+1 (xt+1 ) − kd∗0|t+1 (xt+1 )k2Ψ ≥ Ewt+1 Jt+2 (xt+1 , wt+1 ) where the dependence of the various quantities on xt+1 are added for clarity. Since xt+1 depends on xt and wt , the above can be equivalently written as £ ∗ ¤ ∗ Jt+1 (wt ) − kd∗0|t+1 (wt )k2Ψ ≥ Ewt+1 Jt+2 (wt , wt+1 ) . (45) 19
National University of Singapore
M.E. Dept.
The above inequality holds true for all possible wt , hence ∗ Ewt [Jt+1 (wt )] − Ewt [kd∗0|t+1 (wt )k2Ψ ] £ ∗ ¤ ∗ ≥Ewt [Ewt+1 Jt+2 (wt , wt+1 ) ] = Et [Jt+2 (wt , wt+1 )]
(46)
∗ ∗ Et [Jt+1 (wt )] − Et [kd∗0|t+1 (wt )k2Ψ ] ≥ Et [Jt+2 (wt , wt+1 )]
(47)
or The equality in (46) follows from the i.i.d. assumption in (A3), particularly, £ ∗ ¤ Ewt [Ewt+1 Jt+2 (wt , wt+1 ) ] ·Z ¸ ∗ =Ewt Jt+2 (wt , wt+1 )fwt+1 (wt+1 )dwt+1 Z Z ∗ = Jt+2 (wt , wt+1 )fwt+1 (wt+1 )dwt+1 fwt (wt )dwt Z Z ∗ = Jt+2 (wt , wt+1 )fwt ,wt+1 (wt , wt+1 )dwt+1 dwt ∗ ∗ =Ewt ,wt+1 [Jt+2 (wt , wt+1 )] = Et [Jt+2 (wt , wt+1 )]
where fwt (·), fwt+1 (·) and fwt ,wt+1 (·, ·) are density functions of wt , wt+1 and their joint density function respectively and fwt ,wt+1 (·, ·) = fwt (·)fwt+1 (·) follows from assumption (A3). Summing (43) and (47) leads to ∗ Jt∗ ≥ kd∗0|t k2Ψ + Et [kd∗0|t+1 (wt )k2Ψ ] + Et [Jt+2 (wt , wt+1 )]
Repeating the above procedure infinite times leads to ∞ > Jt∗ ≥ kd∗0|t k2Ψ +
∞ X
h i ∗ Et kd∗0|i (wt , · · · , wi−1 )k2Ψ + lim Et [Jt+r (wt , · · · , wt+r−1 )] r→∞
i=t+1
∗ (w , · · · , w where the left inequality follows from Theorem 7(ii). Using the fact that limr→∞ Et [Jt+r t t+r−1 )] > ∗ 2 0 and kd0|t kΨ is finite, we have
∞>
∞ X
h i Et kd∗0|i (wt , · · · , wi−1 )k2Ψ
i=t+1
By applying Markov bound (given non-negative random variable η and any ² ≥ 0, E[η] ≥ ²Pr{η ≥ ²}) and considering kd∗0|i k2Ψ as a random number, we have ∞>²
∞ X
Pr(kd∗0|i k2Ψ ≥ ²)
(48)
i=t+1
for any arbitrary small ² > 0. From the First Borel-Cantelli Lemma [26], this implies that limi→∞ Pr(kd∗0|i k2Ψ ≥ ²) = 0. Hence d0|i → 0 with probability one as t increases. Consequently, the MPC control law (22) converges to Kf xt with probability one. When this happens, the closed-loop system converges to xt+1 = Φxt +wt and, hence, xt converges to F∞ (Kf ) with probability one. (iii) follows directly from (ii) and assumption (A4) that F∞ (Kf ) ⊂ int(Xf ). 20
National University of Singapore
E
M.E. Dept.
Proof of Theorem 4
Proof. (i) The replacement of cost function J(d, D) by V (d, D) does not affect the feasibility of problem PN (x). This means that part (i) of Theorem 3 remains valid. (ii) Let Vt∗ and Vˆt+1 be defined in the same manner as Jt∗ and Jˆt+1 in the statement of proofs of Theorem 3. Following the same reasoning as in (40), it can be shown that Vt∗ − Vˆt+1 (wt ) = kd∗0|t k2Ψ + p(wt )
(49)
where p(wt ) =
N −1 X
i ∗ 2 i ∗ i ∗ i ∗ (γ1 kvec(Di|t ) k + γ2 kvec(Di|t ) k − 2(d∗i|t )T Ψ(Di|t ) wt − k(Di|t ) wt k2Ψ ).
(50)
i=1
Hence p(wt ) ≥ ≥ = =
N −1 X
i ∗ 2 i ∗ i ∗ i ∗ 2 (γ1 kvec(Di|t ) k + γ2 kvec(Di|t ) k − 2kd∗i|t kkΨkkwt kk(Di|t ) k − kΨkkwt k2 k(Di|t ) k )
i=1 N −1 X
i ∗ 2 i ∗ i ∗ i ∗ 2 (γ1 kvec(Di|t ) k + γ2 kvec(Di|t ) k − 2αβkΨkk(Di|t ) kF − α2 kΨkk(Di|t ) kF ) (51)
i=1 N −1 X
i ∗ 2 i ∗ i ∗ i ∗ 2 (γ1 kvec(Di|t ) k + γ2 kvec(Di|t ) k − 2αβkΨkkvec(Di|t ) k − α2 kΨkkvec(Di|t ) k )
i=1 N −1 X
i ∗ 2 i ∗ ((γ1 − α2 kΨk)kvec(Di|t ) k + (γ2 − 2αβkΨk)kvec(Di|t ) k) ≥ 0
i=1 i )∗ k ≥ k(D i )∗ k and k(D i )∗ k = kvec(D i )∗ k are used. Hence, p(w ) ≥ where the facts k(Di|t t F F i|t i|t i|t 0 under (29). As a consequence, equation (49) implies ∗ Vt∗ − kd∗0|t k2Ψ ≥ Vt+1 ≥0
(52)
Hence, {Vt∗ } is a monotonic non-increasing sequence and is bounded from below by zero. This means that V∞ := limt→∞ Vt∗ ≥ 0 exists. Repeating (52) for t from 0 to ∞ and summing them up, it follows that ∞ X ∗ ∞ > V0 − V∞ ≥ kd∗0|t k2Ψ (53) t=0 limt→∞ d∗0|t
Since Ψ is positive definite, this implies that = 0 and limt→∞ ut = Kf xt . Therefore, the stated result follows. (iii) follows (ii) and assumption (A4) that F∞ (Kf ) ⊂ int(Xf ).
F
Computation of β
β := max(x,d,D)∈TN ,i∈ZN −1 kdi k = max(x,d,D)∈TN kd0 k is due to the fact that for any (x, d, D) ∈ ¯ D) ¯ ∈ TN can be found such that d¯0 = di . SpecifiTN and integer i ∈ Z+ x, d, N −1 , a set of (¯ cally, given (x, d, D) ∈ TN and let the correspondingly defined state and control sequence be 21
National University of Singapore
M.E. Dept.
{x0 , . . . , xN } and {u0 , . . . , uN −1 }. According to (16) xN ∈ Xf for all possible disturbances. ¯ D) ¯ can be defined by Then for any i ∈ Z+ x, d, N −1 , (¯ i
x ¯ = Φ x+
i−1 X j=0
i−1−j
Φ
( dj+i Bdj , d¯j = 0
j ∈ ZN −1−i ¯ jk = ,D N −i≤j ≤N −1
( k Dj+i 0
j ∈ Z+ N −1−i k ∈ Z+ j . N −i≤j ≤N −1
¯ D) ¯ define the control sequence where x ¯ is the nominal state of xi defined by (x, d, D) and (d, {ui , . . . , uN −1 , Kf xN , . . . , Kf xN −1+i }. According to (A4) under controller ut = Kf xt all the constraints are ¯ D) ¯ satisfies (12)-(16), namely satisfied and xt ∈ Xf for t ≥ N since xN ∈ Xf . Therefore, (¯ x, d, ¯ ¯ (¯ x, d, D) ∈ TN . As a result, max(x,d,D)∈TN kd0 k ≥ max(x,d,D)∈TN kdi k, for any i ∈ ZN −1 and β = max(x,d,D)∈TN kd0 k.
22
Technical Report
National University of Singapore Department of Mechanical Engineering C09-001
National University of Singapore
M.E. Dept.
Abstract This paper shows new convergence properties of constrained linear discrete time system with bounded disturbances under Model Predictive Control (MPC) law. The MPC control law is obtained using an affine disturbance feedback parametrization with an additional linear state feedback term. This parametrization has the same representative ability as some recent disturbance feedback parametrization, but its choice together with an appropriate cost function results in a different closed-loop convergence property. More exactly, the state of the closed-loop system converges to a minimal invariant set with probability one. Deterministic convergence to the same minimal invariant set is also possible if a less intuitive cost function is used. Numerical experiments are provided that validate the results.
2
National University of Singapore
1
M.E. Dept.
Introduction
This paper considers the system: xt+1 = Axt + But + wt ,
(1)
(xt , ut ) ∈ Y, wt ∈ W, ∀ t ≥ 0
(2)
where xt ∈ Rn , ut ∈ Rm and wt ∈ W ⊂ Rn are the state, control and disturbance of the system at time t respectively, and Y represents the joint constraint set on xt and ut . The study of such a system under the Model Predictive Control (MPC) framework has been an active area of research. See, for example, [3, 8, 23, 15, 18, 1, 19] and the references cited therein. An important research issue remains the choice of control parametrization within the control horizon. Several choices have been proposed in the literature [23, 8, 17, 18] and a popular choice is ut = Kxt + ct where K is a fixed feedback gain and ct is the new optimization variable. However, such a choice is known to be conservative [5, 17, 12, 24, 25] and its use will result in a relatively small domain of attraction. A natural extension of the fixed-gain parametrization is the time-varying affine state feedback ut = Kt xt + ct where both Kt and ct are the optimization variables within the horizon. However, such a parametrization is not computationally amiable as the resulting optimization problem is not convex. More recently, control parametrization based on affine function of disturbances have appeared [17, 4, 12, 24, 25]. This parametrization is appealing as the resulting problem is convex and solvable via standard numerical routines. Specifically, [17] proposes the control parametrization uL i
=
i X
Mij wi−j + vi ,
i = 0, · · · , N − 1
(3)
j=1
where Mij and vi are the optimization variables and N is the length of the horizon. [12] show that (3) is equivalent, in terms of the set of states reachable within the horizon, to that of timevarying affine state feedback. They also show that, under mild assumptions, the origin of the closed-loop system is input-to-state stable (ISS) with respect to the disturbance input under the MPC control law derived using (3) and a cost function that corresponds to the Linear Quadratic (LQ) cost for system (1) with wt ≡ 0. Recently, [24] 2008 propose an extended disturbance feedback parametrization uW i = Kf xi + ci +
N −1 X
Cij wi−j , i = 0, · · · , N − 1
(4)
j=1
where Kf is a fixed feedback gain, ci and Cij are the optimization variables. They show that parametrization (4) under the MPC framework has the same domain of attraction as that of using (3). Using an appropriate cost function, they also show a stronger stability result: state of the closed-loop system converges to the minimal disturbance invariance set, F∞ , of the system xt+1 = (A + BKf )xt + wt . Unlike (3), index j runs from 1 to N − 1 in (4) and therefore i − j can be negative in wi−j . When this happens, wi−j refers to past realized disturbances. This also means that the resulting MPC control law derived from (4) requires the values of xt and wt−1 , · · · , wt−N +1 for its evaluation at time t. 3
National University of Singapore
M.E. Dept.
This work proposes a new control parametrization based on (4) and a new cost function. The use of which results in an MPC control law requiring only the measurement of xt for its evaluation at time t. The resulting closed-loop system has the same domain of attraction as [17] and [12] but with a stronger convergence result: the closed-loop system state converges to F∞ with probability one; and deterministic convergence to the same set if a less intuitive cost function is used. The rest of this paper is organized as follows. This section ends with notations used, assumptions needed and a brief review of standard results. Section 2 states the control parametrization, the finite horizon (FH) optimization problem and the cost function used. Section 3 discusses the computation of the MPC problem. The probabilistic convergence of the state of the closed-loop system is given in section 4. Section 5 shows a formulation that strengthens the convergence result under a weaker set of assumptions. This, however, requires the use of a somewhat less intuitive cost function. Numerical examples are the contents of section 6 and they are followed by the conclusions. The following notations are used. Z+ k := {1, · · · , k} and Zk := {0, 1, · · · , k} are the respective sets of positive and non-negative integers up to k. k · k is the standard 2-norm for matrices and vectors. Given matrices A ∈ Rn×m and B ∈ Rp×q : kAkF is the Frobenious £ ¤T norm, vec(A) = AT1 · · · ATm ∈ Rnm where Ai is the ith column of A is the stacked vector of columns of A and A ⊗ B ∈ Rnp×qm is the Kronecker product of A and B. A Â (º)0 means that square matrix A is positive definite (semi-definite). For any A Â 0, kxk2A = xT Ax. 1r is a r-element column vector with all elements being 1 and In is the n × n identity matrix. For any set X, Y ⊂ Rn , X + Y := {x + y : x ∈ X, y ∈ Y } is the Minkowski sum of X and Y . The system (1)-(2) is assumed to satisfy the following assumptions: (A1) system (A, B) is stabilizable; (A2) the set Y is a polytope having a characterization Y := {(x, u)| Yx x + Yu u ≤ 1q } ⊂ Rn+m
(5)
for some Yx ∈ Rq×n and Yu ∈ Rq×m ; (A3) the disturbance wt , t ≥ 0 are independent and identically distributed (i.i.d.) with zero mean and W is a polytope characterized by W := {w| Hw ≤ 1r } ⊂ Rn ,
(6)
for some H ∈ Rr×n ; (A4) a compact and non-empty constraint-admissible disturbance invariant set exists for system (1)-(2) under the feedback law u = Kf x and takes the form Xf := {x| Gx ≤ 1g } ⊂ Rn
(7)
for some G ∈ Rg×n and that Xf contains the minimal disturbance invariant set of xt+1 = (A + BKf )xt + wt in its interior. 4
National University of Singapore
M.E. Dept.
The above assumptions can be rationalized in the following ways. (A1) is standard. The characterization of Y in (A2) is made out of the need for a concrete computational representation. (A3) is a typical assumption on the disturbance and has been used in several prior works [11, 22]. That W is a polytope and contains the origin in its interior is an assumption made out of convenience for the presentation of the main result. Relaxation of (A3) is possible and the details are given in section 5 and illustrated via a numerical example in section 6. The existence of Xf in (A4) has been shown by [13] 1998 provided that W is sufficiently small. More exactly, for any feedback gain Kf ∈ Rm×n such that Φ := A + BKf is strictly stable and sufficiently small W , Xf is the maximal, constraint admissible and disturbance invariant in the sense that Φx + w ∈ Xf , (x, Kf x) ∈ Y for all x ∈ Xf and for all w ∈ W . It is also known [14, 9] that the state of the system xt+1 = Φxt + wt converges to the minimal disturbance invariant set, F∞ , given by F∞ (Kf ) := W + (A + BKf )W + (A + BKf )2 W + · · ·
(8)
and that F∞ is compact. The assumption that F∞ (Kf ) ⊂ Xf is also a natural consequence when W is not too large.
2
Control parametrization
MPC formulation solves an N -stage finite horizon (FH) optimization problem. Let xi and ui , i ∈ ZN −1 denote the predicted state and predicted control at the ith stage respectively within the horizon. The proposed control parametrization within the FH optimization problem takes the form i X ui = Kf xi + di + Dij wi−j for all i ∈ ZN −1 (9) j=1
where di ∈ Rm , Dij ∈ Rm×n , j ∈ Z+ i , i ∈ ZN −1 are the variables of the FH problem and Kf is the given feedback gain in (A4). Since i − j ≥ 0, wi−j is the (i − j)th predicted disturbance at each stage i. In this regard, (9) is similar to (3) in that only predicted disturbances are used W in the parametrization. In addition, ui is equivalent to uL i and ui in terms of the family of functions that they represent. To state this precisely, let 0 ··· 0 0 D1 ··· 0 0 £ T T ¤T 1 T Nm N m×N n d := d0 d1 · · · dN −1 ∈ R , D := (10) .. .. .. ∈ R .. . . . . N −1 DN −1 · · ·
1 DN −1 0
and v and M to be similarly defined as d and D in structure but having entries vi instead W of di and Mij instead of Dij respectively. The equivalence of ui , uL i and ui are stated in the following theorem. N −1 Theorem 1. Suppose x0 and Kf are given. Then (i) the set {ui }i=0 of (9) is equivalent N −1 L to {ui }i=0 of (3) in the sense that for any (v, M), there exists a unique (d, D) such that −1 N −1 L N −1 W N −1 {ui }N i=0 = {ui }i=0 and vice versa. (ii) {ui }i=0 is equivalent to {ui }i=0 of (4).
5
National University of Singapore
M.E. Dept.
Proof. See Appendix A. The FH optimization problem under parametrization (9), referred hereafter as PN (x), is min J(d, D)
(11)
s.t. x0 = x,
(12)
d,D
xi+1 = Axi + Bui + wi , i ∈ ZN −1 ui = Kf xi + di +
i X
Dij wi−j , i ∈ ZN −1
(13) (14)
j=1
(xi , ui ) ∈ Y, ∀wi ∈ W, i ∈ ZN −1
(15)
x N ∈ Xf ,
(16)
∀wi ∈ W, i ∈ ZN −1
where Y and Xf are the corresponding sets given by (5) and (7) respectively. The cost function J(d, D) takes the form N −1 i X X kdi k2Ψ + J(d, D) := (17) kvec(Dj )k2Λ i
i=0
j=1
for any choice of Ψ ∈ Rm×m and Λ ∈ Rmn×mn that satisfy Ψ = ΨT Â 0, Λ º Σw ⊗ Ψ
(18)
where Σw is the covariance matrix of wt and vec(·) is stacking operator defined in Section 1. Clearly, J(d, D) is a measure of the deviation of ui of (14) from the linear control law Kf xi and the motivation for it as the objective function is clear: penalizing the use of nonzero (d, D) forces the asymptotic behavior of the closed-loop system to approach that of xt+1 = (A + BKf )xt . The technical condition (18) is to ensure convergence of the closed-loop states and its exact role will become clear in the proof of Theorem 3. Several comments on J(d, D) are in order. Remark 1. A connection between J(d, D) and the standard LQ cost can be established. Specifically, suppose Q º 0 and R Â 0 are given and let P = AT P A−AT P B(R+B T P B)−1 B T P A+Q, the solution of the algebraic Riccati equation, Ψ = R + B T P B, Λ = Σw ⊗ Ψ and Kf = −(R + B T P B)−1 B T P A. It is shown [10, 24, 25] that "N −1 # X Ew (kxi k2Q + kui k2R ) + kxN k2P = xT0 P x0 + N trace(Σw P ) + J(d, D) (19) i=0
where Ew is the expectation taken over {w0 , · · · , wN −1 } within the horizon. Since the first two terms on the right hand side of (19) are independent of (d, D), minimizing J(d, D) is equivalent to minimizing the expected infinite horizon LQ cost over the disturbance input. Remark 2. From (18) and (19), it may appear that Σw is needed for the determination of Λ. This is not true. The choice of Λ can be made to satisfy (18) even when Σw is not known accurately. One simple choice is to let Λ = α2 In ⊗ Ψ where α := maxw∈W kwk. Then it follows that Λ º Σw ⊗ Ψ because α2 In º wwT for all w ∈ W and α2 In ⊗ Ψ º E[wwT ] ⊗ Ψ. Consequently, (A3) provides for conditions that guarantee the computability of maxw∈W kwk. 6
National University of Singapore
M.E. Dept.
Further discussion on the choice of Ψ and Λ and their influence on closed-loop system trajectories are discussed in section 4. Several associated sets, needed to facilitate the discussions in the sequel, are introduced. Let the feasible set of optimization problem PN (x) be TN := {(x, d, D)|(d, D) is feasible to PN (x)}
(20)
and the set of admissible initial states, or equivalently, the domain of attraction of the MPC controller is XN := {x|∃ (d, D) such that (x, d, D) ∈ TN }. (21) The rest of the MPC formulation is standard: PN (xt ) is solved at each time t to obtain the optimizer (d∗t , D∗t ) := (d∗ (xt ), D∗ (xt )) and the first control, u∗0 , is applied to (1) at time t resulting in the MPC control law, ut = u∗0 = Kf xt + d∗0 .
(22)
Remark 3. From (20) and (21), it is easy to see that TN and XN depend only on the conW straints (12)-(16). If uL i or ui replaces ui in (14), it follows from Theorem 1 that the corresponding domain of attraction, XNL and XNW is the same as XN . However, the stability of the corresponding closed-loop systems can differ. See Remark 6.
3
The Computation of PN (x)
Following (5), (7), (12)-(14), inequalities (15) and (16) can be collectively restated as £ ¤ ¯ + Bd ¯ + max BD ¯ + G¯ w ≤ 1s Ax w∈WN
(23)
£ ¤T T N ¯ ¯ where s = N q + g, w := w0T w1T · · · wN −1 , W := {w|Hw ≤ 1` } with H = IN ⊗ H, ` = ¯ B¯ and G¯ are appropriate matrices given in Appendix B and the max N r following (6), A, operator is meant to be taken row-wise. Correspondingly, the ith row of (23) can be rewritten ¯ ≤ 1` } ≤ bi for some ei ∈ RN m , bi ∈ R that depend on x, d and D. Let as max{eTi w|Hw ¯ Then, it follows by duality zi ∈ R` be the Lagrange multiplier corresponding to the rows of H. T T T ¯ ¯ that max{ei w|Hw ≤ 1` } = min{zi 1` |H zi = ei , zi ≥ 0}. Collecting over all the rows of (23), the PN (x) can be equivalently stated as min J(d, D)
d,D,Z
(24)
¯ + Bd ¯ + ZT 1` ≤ 1s s.t. Ax ¯ = BD ¯ + G¯ ZT H
(25)
zi ≥ 0, i = 1, . . . , s
(27)
(26)
where Z = [z1 · · · zs ] ∈ R`×s and the minimization of Z is relaxed since the existence of any one feasible Z is enough to guarantee that (x, d, D) ∈ TN . Remark 4. The above duality results can be extended to W sets that are non-polyhedral. See, for example, treatments of such sets in [6, 20]. If W is a second-order cone [16, 2] representable 7
National University of Singapore
M.E. Dept.
bounded set with non-empty interior such that WN = {w| kLi w − li k ≤ λTi w − θi , i ∈ Z+ k} T w| w ∈ for some matrices Li , li , λi and θi , i ∈ Z+ , then it follows from duality that max{e P Pk k T + WN } = min(µi ,ηi ) { ki=1 (µTi li −ηi θi )| i=1 (Li µi −ηi λi ) = e, kµi k ≤ ηi , i ∈ Zk }. Similarly, if W is a bounded set with non-empty interior such that WN = PN nsemi-definite cone representable + N n {Ω ∈ R | i=1 Ωi Ci − F < 0, i ∈ ZN n } where Ci and F are symmetrical matrices of appropriate dimension, then max{eT w| w ∈ WN } = minΥ {Trace(F Υ)| Trace(Ci Υ) = ei , i ∈ Z+ N n , Υ 4 0}. Remark 5. While the duality result is available for W being a second-order or semi-definite cone representable set, the availability of Xf satisfying (A4) deserves some clarifications. When W is non-polyhedral, computation of a constraint-admissible disturbance invariant set Xf may not be easy. A simple approach is to construct a polytope Wp such that Wp ⊃ W and Wp ≈ W . In that case, a Xf satisfying (A4) can be constructed using Wp following existing computational methods [14]. Using this Xf in (16) and Remark 4, PN (x) becomes either a second-order cone or a semi-definite cone programming problem. It is worthy to note that the use of such an Xf in (16) and with wi ∈ W for all i ∈ ZN −1 in both (15) and (16) is less conservative than replacing W by Wp throughout (12)-(16). An example using such an approach is illustrated in Section 6.
4
Feasibility and Probabilistic Convergence
The feasibility of PN (xt ) at different time instants and stability of the closed-loop system under the feedback law (22) are addressed in this section. Theorem 2. Suppose (A1)-(A4) are satisfied, the FH optimization problem PN (x) has the following properties (i) TN is convex and compact. (ii) If x ∈ XN , the optimal solution of PN (x) exists. (iii) If PN (xt ) admits an optimal solution, so does PN (xt+1 ) under the feedback law (22) for all wt ∈ W . Proof. See Appendix C. The main result of probabilistic convergence of the state of the closed-loop system is stated in the next theorem. Theorem 3. Suppose x0 ∈ XN and (A1)-(A4) are satisfied. System (1) under MPC control law (22) obtained from the solution of PN (x) using cost function (17) with condition (18) satisfied has the following properties: (i) (xt , ut ) ∈ Y for all t ≥ 0, (ii) xt → F∞ (Kf ) with probability one as t → ∞ (iii) xt enters Xf in finite time with probability one. Proof. See Appendix D. Remark 6. It is of interest to know if the results of Theorem 3 can be extended to the case where uL i is used in (14) in view of Theorem 1. As seen in the proof of Theorem 3, the convergence property depends on the choices of the cost function and the control parametrization. Since (9) and (3) are equivalent when condition (36) is satisfied, the results of Theorem 3 is applicable L under the following conditions: uL i is used in (14); a new cost function J (M, v) :=J((I + 8
National University of Singapore
M.E. Dept.
KϕB)−1 (v − KϕAx), (I + KϕB)−1 (M − KϕG)) is used in (11) where A, B, K, G and ϕ are those given in (32) and (33); and (M, v) becomes the optimization variables for PN (x). The approach by [12] uses the nominal LQ cost as the cost function and it is not clear if the results of Theorem 3 remains true under that situation. One associated issue in the formulation of PN (x) is the choices of Ψ and Λ in J(d, D). How should Ψ and Λ be chosen and how do these choices affect the closed-loop system trajectories? As xt → F∞ with probability one from result (ii) of Theorem 3, it implies that xt enters Xf with probability one and stays within thereafter since F∞ ⊂ Xf . When this happens, the optimal (d, D) are zero in PN (x) and the MPC control law becomes ut = Kf xt for all t thereafter. The closed-loop system behavior then corresponds to that of the system xt+1 = (A + BKf )xt . Clearly, the choices of Λ and Ψ does not affect the asymptotic behavior of the system but only the transient when xt ∈ / Xf . Suppose Λ = Σw ⊗ Ψ. Then admissible changes in Ψ will not result in changes in the system behavior since Σw ⊗ Ψ is linear in Ψ. On the other hand, if Ψ is fixed, Λ can be chosen to be increasingly ”larger” than Σw ⊗ Ψ. In loose terms, a ”larger” choice of Λ penalizes the use of D versus the use of d in J(d, D). Such a preference would mitigate the effect of the disturbance feedback component in the control parametrization, resulting in a parametrization that is closer in spirit to ut = Kf xt +dt of [8]. When this happens, the transient response for the system may become slower even though the domain of attraction XN remains unaffected. This observation together with the associate details used in the experimental study are discussed in section 6.
5
Deterministic Convergence
While the assumption of W being a compact set is reasonable, the assumption of wt being zero mean and i.i.d. is harder to verify in practice. This section is concerned with the relaxation of assumption (A3) while achieving a stronger convergence result than that of Theorem 3. Consider (A3a) wt ∈ W and W is a polytope characterized by W := {w| Hw ≤ 1r } ⊂ Rn for some H ∈ Rr×n . and define the cost function V (d, D) :=
N −1 X
i X kdi k2Ψ + (γ1 kvec(Dij )k2 + γ2 kvec(Dij )k)
i=0
(28)
j=1
for some constants γ1 and γ2 satisfying γ1 ≥ α2 kΨk, γ2 ≥ 2αβkΨk
(29)
where α := maxw∈W kwk and β := max(x,d,D)∈TN ,i∈ZN −1 kdi k. The existence of α and β are guaranteed by compactness of the W and TN sets, provided for in (A3a) and part (i) of Theorem 2 respectively. 9
National University of Singapore
M.E. Dept.
Theorem 4. Suppose x0 ∈ XN and (A1-A2),(A3a) and (A4) are satisfied and J(d, D) is replaced by V (d, D) in PN (x) satisfying condition (29), then system (1) under the MPC control law (22) satisfies (i) (xt , ut ) ∈ Y for all t ≥ 0, (ii) xt → F∞ (Kf ) as t → ∞ (iii) xt enters Xf in finite time. Proof. See Appendix E Remark 7. Several choices of the cost function of (28) are possible. For example, the results of Theorem 4 remain true if kvec(Dij )k is replaced by kDij k. This may be more appealing as less conservative bounds on γ1 and γ2 can be found to ensure the non-negativity of p(wt ). However, its use will result in a semi-definite programming problem for PN (x) and is less desirable computationally. The use of kvec(Dij )k or kDij kF results in a second-order cone programming for PN (x) and is computationally more amiable. Remark 8. The computation of β can be simplified to β = max(x,d,D)∈TN kd0 k, see Appendix F for details. Note that any upper bound of β can be used to guarantee the results of Theorem 4. One such upper bound is β¯ := kσk where σi := max(x,d,D)∈TN |d0 (i)| and d0 (i) is the ith element of d0 .
6
Numerical Examples
Four experiments are conducted on a system to validate the results of the previous sections. The parameters and constraints of the system are: A = [1.1 1; 0 1.3], B = [1 1]T , Kf = [−0.7434 − 1.0922], Y = {(x, u)| |u| ≤ 1, kxk∞ ≤ 8} and W = {H w| ˜ kwk ˜ ∞ ≤ 0.2} where H = [1 − 0.2; 0 1] and w ˜ ∈ R2 is a random vector uniformly distributed over [−0.2, 0.2] × [−0.2, 0.2] with covariance matrix Σw˜ = 0.0133I2 . Expressed in the form of (5), the set Y is Yx = [0 0; 0 0; 0 1/8; 0 − 1/8; 1/8 0; − 1/8 0] and Yu = [1; −1; 0; 0; 0; 0]. The set Xf is the corresponding maximal constraint-admissible disturbance invariant set of (1) under ut = Kf xt given by Xf = {x| Gx ≤ 14 }, where G = [−0.7434 −1.0922; 0.7434 1.0922; 0.8252 − 0.2391; −0.8252 0.2391]. The first experiment, Experiment I, uses the cost function of (17) with Ψ = 1, Λ = Λop := Σw ⊗ Ψ = [0.0139 − 0.0027; −0.0027 0.0133], N = 8 and x0 = [−4 2]T . The simulation results over 15 different disturbance realizations shown in Fig. 1 to 4 by solid lines. It is clear from Fig. 1 and 2 that the constraints are satisfied by all trajectories, in accordance to property (i) of Theorem 3. In addition, Figure 1 shows the convergence of xt into Fˆ∞ (Kf ), a tight outer bound of F∞ (Kf ) obtained using procedures given in [21]. This convergence is further verified in Fig. 3 and 4 where the plots of dis(xt , Fˆ∞ ) := minx∈Fˆ∞ kx − xt k, the minimum distance to Fˆ∞ , and dt := d∗0|t against increasing t are shown respectively. The case where W is non-polyhedral is shown in Experiment II, in connection to Remarks 4 and 5. A different disturbance characteristic is used here: w is uniformly distributed over ¯ := {w| W ¯ kS1 wk ¯ ≤ 1, kS2 wk ¯ ≤ 1} where S1 = [5 1; 0 2.5] and S1 = [2.5 0.5; 0 5]. Note ¯ , is needed for the computation of Xf that a tight bounding polytope, Wp , such that Wp ⊃ W satisfying (A4) and it corresponds to the W set of Experiment I (see Fig. 5). Also, Ψ and Λ of the first experiment are used and it is easy to verify that condition (18) remains true 10
National University of Singapore
M.E. Dept.
¯ is a second-order cone representable set and the conversion because Σw  Σw¯ . In this case, W ¯ follows the expression in Remarks 4, resulting in PN (x) being a of (15) and (16) for all wi ∈ W second-order cone programming problem. The simulation results with N = 8 and x0 = [2 −1]T for 15 different realizations of {wt } are plotted in Fig. 1 to 4 using dash-dot lines. Experiment III is designed to understand the influence of Λ and Ψ of (17) on the performance of the closed-loop system. As stated in section 4, choices of these matrices affect only the transient behavior when xt ∈ / Xf and not the asymptotic behavior of the closed-loop system. To quantify the transient, the average number of time step, tf (x0 ), taken to enter Xf from a given x0 is reported. Here, the average is taken over different realizations of the disturbances. Without loss of generality, values of Λ is increased from Λop (see discussion in section 4). Table 1 shows the values of tf (x0 ) and the associated standard deviations over 20 disturbance realizations for several choices of x0 , N and Λ. For each x0 , the same 20 disturbance realizations are used for the different Λ in computing tf (x0 ) and the standard deviations. From the table, tf (x0 ) generally increases when Λ increases. For comparison purpose, the corresponding trajectories of the system under same settings as the Experiment I except for Λ = 104 Λop are plotted in Fig. 1 to 4 using dash lines. From Fig. 3 and 4, the slower convergence of the state and control trajectories are clearly evident. 3
x(2)
2
Xf
1
0
Fˆ∞
−1
−2 −4
−3
−2
−1
0
1
2
x(1)
Figure 1: State trajectories of the first three experiments: solid line for first experiment, dashdot line for the second and dash line for the third. The last experiment, Experiment IV, considers the case discussed in Section 5. The system parameters are the same as those in the first experiment except that the distribution of w˜ is assumed to be unknown. The parameters of (29) are: α = 0.3124 and β = 2.7307 (when N = 8) and β = 3.5425 (when N = 10). Correspondingly, the weight matrices of (28) are Ψ = 1, γ1op := α2 kΨk = 0.0976, γ2op := 2αβkΨk = 1.7059 (2.213 when N = 10). The values of γ1 and γ2 are increased separately and jointly to assess their influence on the system behavior. The general effect of increasing values of γ1 and γ2 appears to have similar trend on the system as the increase in Λ. The time taken to reach Xf from any given x0 increases, although to a lesser percentage than that by Λ, with increasing values of γ1 and γ2 with γ2 having a heavier influence.
11
National University of Singapore
M.E. Dept.
1
ut
0.5
0
−0.5
−1 1
2
3
4
5
6
7
8
9
10
t
Figure 2: Control trajectories of the first three experiments: solid line for first experiment, dashdot line for the second and dash line for the third.
4
dis(xt , Fˆ∞ )
3
2
1
0
−1
0
2
4
6
8
10
t
Figure 3: Distance between states and F∞ (Kf ) of the first three experiments: solid line for first experiment, dashdot line for the second and dash line for the third.
12
National University of Singapore
M.E. Dept.
0.5
dt = d∗0|t
0
−0.5
−1
−1.5
−2
1
2
3
4
5
6
7
8
9
10
t
Figure 4: Values of dt of the first three experiments: solid line for first experiment, dashdot line for the second and dash line for the third.
Wp = W
0.2 0.15
w(2)
0.1 0.05 0 −0.05 −0.1 −0.15 −0.2 −0.25
¯ W −0.2
−0.1
0
0.1
w(1)
¯ set. Figure 5: Wp set and W
13
0.2
0.3
N 8 8 10 10
x0
[−4 2]T
[−2.5 − 1.2]T
[−4 − 1]T
[−6 2]T
Initial Condition 102 Λop 4.9 (0.4472) 5.15 (0.4894) 6.75 (0.6387) 7.15 (0.6708)
Λop 4.2 (0.4104) 4.9 (0.5525) 6.3 (0.7327) 6.45 (0.6863)
Λ
8.55 (0.6863)
7.45 (0.8256)
5.15 (0.4894)
5.65 (0.4894)
104 Λop
7.9 (0.4472)
7.2 (0.6959)
5.15 (0.4894)
5.2 (0.4104)
(γ1op , γ2op )
8.05 (0.6863)
7.25 (0.7164)
5.15 (0.4894)
5.4 (0.5026)
(10γ1op , γ2op )
8.6 (0.5982)
7.45 (0.8256)
5.15 (0.4894)
5.65 (0.4894)
(γ1op , 10γ2op )
(γ1 , γ2 )
Table 1: Average time step, tf (x0 ), and its standard deviation
8.6 (0.5982)
7.45 (0.8256)
5.15 (0.4894)
5.65 (0.4894)
(10γ1op , 10γ2op )
National University of Singapore M.E. Dept.
14
National University of Singapore
7
M.E. Dept.
Conclusions
Convergence results for constrained linear system under MPC control law using a new control parametrization and a new cost function are presented. The parametrization uses affine disturbance feedback together with a linear state feedback term, Kf x, and is a modification of the parametrization by [25]. Such a parametrization is similar to existing disturbance parameterizations in the literature in terms of the set of state reachable. Using the proposed cost function and the new parametrization, the closed-loop system state converges to the minimal robust invariant set F∞ (Kf ) with probability one. Deterministic convergence to F∞ (Kf ) is also possible using a less intuitive cost function. The asymptotic behavior of the closed-loop system is determined by the choice of Kf so long as the weight matrices of the cost function satisfy some mild conditions.
References [1] Alamo, T., Mu˜ noz de la Pe˜ na, D., Limon, D. and Camacho, E. F. [2005]. Constrained min-max predictive control: modifications of the objective function leading to polynomial complexity, IEEE Transaction on Automatic Control 50(5): 710–714. [2] Alizadeh, F. and Goldfarb, D. [2003]. Second-order cone programming, Methimetical Programming 95(1): 3–51. [3] Bemporad, A. [1998]. Reducing conservativeness in predictive control of constrained systems with disturbances, Proceedings of 37th Conference on Decision and Control, Tampa, Florida, pp. 1384–1389. [4] Ben-Tal, A., Goryashko, A., Guslitzer, E. and Nemirovski, A. [2004]. Ajustable robust solutions of uncertaint linear programs, Mathematical Programming 99(2): 351–376. [5] Ben-Tal, A. and Nemirovski, A. [1999]. Robust solutions to uncertain linear programs, OR Letters 25. [6] Ben-Tal, A. and Nemirovski, A. [2001]. Lectures on Modern Convex Optimization: Analysis, Algorithms, Engineering Applications, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA. [7] Bertsekas, D. P. [2003]. Nonlinear Programming, 2nd edn, Athena Scientific. [8] Chisci, L., Rossiter, J. A. and Zappa, G. [2001]. Systems with persistent disturbances: predictive control with restricted constraints, Automatica 37(7): 1019–1028. [9] Gilbert, E. and Ong, C. J. [2008]. Linear systems with hard constraints and variable set points: Their robustly invariant sets, Technical report, Department of Aerospace Engineering, University of Michigan and Department of Mechanical Engineering, National University of Singapore, available at http://guppy.mpe.nus.edu.sg/ mpeongcj/ongcj.html. C08-002. 15
National University of Singapore
M.E. Dept.
[10] Goulart, P. J. and Kerrigan, E. C. [2005]. On the stability of a class of robust receding horizon control laws for constrained systems, Technical report, Department of Engineering, University of Cambridge. CUED/F-INFENG/TR.532. [11] Goulart, P. J. and Kerrigan, E. C. [2006]. Robust receding horizon control with an expected value cost, Proc. UKACC International Conference (Control 2006), Glasgow, Scotland. [12] Goulart, P. J., Kerrigan, E. C. and Maciejowski, J. M. [2006]. Optimization over state feedback policies for robust control with constraints, Automatica 42(4): 523–533. [13] Kolmanovsky, I. and Gilbert, E. G. [1995]. Maximal output admissible sets for discretetime systems with disturbance inputs, Proceedings of the 1995 American Control Conference, Seattle, pp. 1995–1999. [14] Kolmanovsky, I. and Gilbert, E. G. [1998]. Theory and computation of disturbance invariant sets for discrete-time linear systems, Mathematical Problems in Engineering 4(4): 317– 367. [15] Lee, Y. I. and Kouvaritakis, B. [1999]. Constrained receding horizon predictive control for systems with disturbances, International Journal of Control 72(11): 1027–1032. [16] Lobo, M. S., Vandenberghe, L., Boyd, S. and Lebret, H. [1998]. Applications of secondorder cone programming, Linear Algebra and its Applications 284: 193–228. [17] L¨ofberg, J. [2003]. Approximations of closed-loop minimax mpc, Proceedings of the 42nd IEEE Conference on Decision and Control, Maui, Hawaii, USA, pp. 1438–1442. [18] Mayne, D. Q., Seron, M. M. and Rakovi´c, S. V. [2005]. Robust model predictive control of constrained linear systems with bounded disturbances, Automatica 41(2): 219–224. [19] Mu˜ noz de la Pe˜ na, D., Alamo, T., Bemporad, A. and Camacho, E. F. [2006]. A decomposition algorigthm for feedback min-max model predictive control, IEEE Transaction on Automatic Control 51(10): 1688–1692. [20] Nemirovski, A. [2006]. Advances in convex optimization: conic programming, Proceedings of the International Congress of Mathematicians, Madrid, pp. 413–444. [21] Ong, C. J. and Gilbert, E. G. [2006]. The minimal disturbance invariant set: Outer approximations via its partial sums, Automatica 42(9): 1563–1568. [22] Primbs, J. A. [2007]. A soft constraint approach to stochastic receding horizon control, Proceedings of the 46th IEEE Conference on Decision and Control, New Orleans, Louisiana, USA, pp. 4797 – 4802. [23] Rossiter, J. A., Kouvaritakis, B. and Rice, M. J. [1998]. A numerically robust state-space approach to stable predictive control strategies, Automatica 34(1): 65–73. [24] Wang, C., Ong, C. J. and Sim, M. [2007]. Model predictive control using affine disturbance feedback for constrained linear system, Proceedings of the 46th IEEE Conference on Decision and Control, New Orleans, Louisiana, USA, pp. 1275–1280. 16
National University of Singapore
M.E. Dept.
[25] Wang, C., Ong, C. J. and Sim, M. [2008]. Constrained linear system with disturbance: stability under disturbance feedback, To appear in Automatica . [26] Williams, D. [1991]. Probability with Martingales, Cambridge University Press.
A
Proof of Theorem 1
Proof. (i) Suppose the predicted states, predicted controls and predicted disturbances within £ ¤T £ ¤T the horizon are x := xT0 xT1 · · · xTN ∈ R(N +1)n , u := uT0 uT1 · · · uTN −1 ∈ RN m , and £ ¤T N −1 T w := w0T w1T · · · wN ∈ RN n . The state x and the and the control sequence, {ui }i=0 , −1 defined by (9) can be equivalently stated as
where A :=
In A A2 .. .
, B :=
AN
0 B AB .. .
x = Ax0 + Bu + Gw
(30)
u = Kx + d + Dw
(31)
0 0 B .. .
··· ··· ··· .. .
0 0 0 .. .
AN −1 B AN −2 B · · ·
B
, G :=
0 I A .. .
0 0 I .. .
··· ··· ··· .. .
0 0 0 .. .
AN −1 AN −2 · · ·
I
, K = [IN ⊗ Kf 0] (32)
and d and D are those given by (10). Using (31) in (30), we get x = ϕAx0 +ϕBd+(ϕBD+ϕG)w where ϕ = (I − BK)−1
(33)
u = KϕAx0 + (I + KϕB)d + [(I + KϕB)D + KϕG] w
(34)
and u becomes N −1 Consider the parametrization (3) and the sequence {uL i }i=0 expressed by the variables (M, v) of (10). It follows that uL = v + Mw (35)
Comparing (34) and (35), u = uL if and only if ( KϕAx0 + (I + KϕB)d = v (I + KϕB)D + KϕG = M.
(36)
Note that (I + KϕB) is a lower triangular matrix with all diagonal elements being 1 and is always invertible. In addition, KϕG is a strict lower triangular block matrix like M and its multiplication by (I + KϕB)−1 on the left also results in a strict lower triangular block matrix. Hence, the mapping between (M, v) and (D, d) by (36) is unique or one-to-one for all choices −1 L N −1 of K and x0 . This establishes the equivalence of {ui }N i=0 and {ui }i=0 . 17
National University of Singapore
M.E. Dept.
(ii) Set Cij = 0 for all j > i in (4) and it follows that ui is a special case of uW i . To show the converse, let ( P −1 j di = ci + N j=i+1 Ci wi−j , i ∈ ZN −1 (37) Dij = Cij j ≤ i, i ∈ ZN −1 N −1 W N −1 for any ci , Cij that defines uW i . This establishes the equivalence of {ui }i=0 and {ui }i=0 .
B
Expressions of Matrices in (23)
¯ B¯ and G¯ in (23) are The A, · A¯ = Y
ϕA KϕA
¸
· , B¯ = Y
ϕB I + KϕB
¸
· , G¯ = Y
ϕG KϕG
¸
· , Y=
IN ⊗ Yx 0 IN ⊗ Yu 0 G 0
¸
where A, B, G, K and ϕ are defined in (32) and (33).
C
Proof of Theorem 2
Proof. (i) Since Y is compact from (A2), the projection of Y onto x and u space, denoted by XY and UY respectively, are bounded. From (9) and the fact that W is a polytope in Rn , Dij must be bounded in order for ui ∈ UY . Since the origin is inside W , Kf xi + di must be inside UY . Therefore, di is bounded as xi and UY are bounded. This, together with TN being closed and convex from (25)-(27) leads to the desired result.(ii) Since x ∈ XN , PN (x) is feasible. From (i), this means that ΠN (x) := {(d, D)|(x, d, D) ∈ TN } is compact. This, together with the fact that J(d, D) is continuous with respect to (d, D) means that the optimal solution exists by Weierstrass’ Theorem [7]. (iii) The proof follows standard arguments but the details are given for their relevance to Theorem 3. For clarity, additional subscripts “|t” and “|t + 1” are used to denote the variables at the different times. Let (d∗t , D∗t ) denote the optimal solution ˆ t+1 , D ˆ t+1 ) is chosen as of PN (xt ). At time t + 1, wt is realized and (d dˆi|t+1 = ˆj D i|t+1 =
( i+1 ∗ d∗i+1|t + (Di+1|t ) wt 0 ( j (Di+1|t )∗ 0
i ∈ ZN −2 i=N −1
+ j ∈ Z+ i , i ∈ ZN −2
j ∈ Z+ N −1 , i = N − 1.
(38) (39)
ˆ t+1 , D ˆ t+1 ) is feasible to PN (xt+1 ) for all possible wt ∈ W due to the disturbance This choice of (d invariance of Xf for system (1) under control law ut = Kf xt and that (d∗t , D∗t ) is the optimal solution at time t. That the the optimum of PN (xt+1 ) exists follows from the compactness of ΠN (xt+1 ) and the Weierstrass’ theorem [7]. 18
National University of Singapore
D
M.E. Dept.
Proof of Theorem 3
Proof. (i) The stated result follows directly from Theorem 2. (ii) Let Jt∗ := J(d∗t , D∗t ) and ˆ t+1 (wt ), D ˆ t+1 (wt ), D ˆ t+1 ) where (d ˆ t+1 ) are given by (38)-(39). Then it follows Jˆt+1 (wt ) := J(d that Jt∗ − Jˆt+1 (wt ) =
N −1 X
(kd∗i|t k2Ψ
− kdˆi|t+1 k2Ψ ) +
i=0
N −1 X
i ∗ 2 kvec(Di|t ) kΛ
i=1
=kd∗0|t k2Ψ + =kd∗0|t k2Ψ
+
N −1 X
N −1 X
i=1 N −1 X
i=1
(kd∗i|t k2Ψ − kdˆi−1|t+1 k2Ψ ) + (kd∗i|t k2Ψ
−
kd∗i|t
+
i ∗ (Di|t ) wt k2Ψ )
i=1
=kd∗0|t k2Ψ where g(wt ) =
i ∗ 2 kvec(Di|t ) kΛ
+
N −1 X
i ∗ 2 kvec(Di|t ) kΛ
i=1
+ g(wt ) N −1 X
(40)
i ∗ 2 i ∗ i ∗ (kvec(Di|t ) kΛ − 2(d∗i|t )T Ψ(Di|t ) wt − k(Di|t ) wt k2Ψ ).
(41)
i=1
Taking the expectation of (40) over wt , it follows that h i Jt∗ − kd∗0|t k2Ψ = Ewt Jˆt+1 (wt ) + Ewt [g(wt )] h i ≥ Ewt Jˆt+1 (wt ) £ ∗ ¤ £ ∗ ¤ (wt ) . (wt ) = Et Jt+1 ≥ Ewt Jt+1
(42) (43)
where Et in (43) is the expectation taken over wi , i ≥ t. Inequality (42) follows from the fact that Ewt [g(wt )] ≥ 0. This is so because by taking the expectation of (41), one gets Ewt [g(wt )] =
N −1 X
i ∗ 2 i ∗ 2 i ∗ (kvec(Di|t ) kΛ − kvec(Di|t ) kΣw ⊗Ψ − 2(d∗i|t )T Ψ(Di|t ) E[wt ])
(44)
i=1
where the last term is zero due to (A3) and the rest is non-negative due to (18). The inequality ∗ (w ) for every w ∈ W which implies that in (43) follows from the fact that Jˆt+1 (wt ) ≥ Jt+1 t t ∗ ∗ (w ) depends on ˆ Ewt [Jt+1 (wt )] ≥ Ewt [Jt+1 (wt )]. Equality (43) follows from the fact that Jt+1 t wt only and not on any wi , i > t. Repeating the inequality of (43) for increasing t, one gets, £ ∗ ¤ ∗ Jt+1 (xt+1 ) − kd∗0|t+1 (xt+1 )k2Ψ ≥ Ewt+1 Jt+2 (xt+1 , wt+1 ) where the dependence of the various quantities on xt+1 are added for clarity. Since xt+1 depends on xt and wt , the above can be equivalently written as £ ∗ ¤ ∗ Jt+1 (wt ) − kd∗0|t+1 (wt )k2Ψ ≥ Ewt+1 Jt+2 (wt , wt+1 ) . (45) 19
National University of Singapore
M.E. Dept.
The above inequality holds true for all possible wt , hence ∗ Ewt [Jt+1 (wt )] − Ewt [kd∗0|t+1 (wt )k2Ψ ] £ ∗ ¤ ∗ ≥Ewt [Ewt+1 Jt+2 (wt , wt+1 ) ] = Et [Jt+2 (wt , wt+1 )]
(46)
∗ ∗ Et [Jt+1 (wt )] − Et [kd∗0|t+1 (wt )k2Ψ ] ≥ Et [Jt+2 (wt , wt+1 )]
(47)
or The equality in (46) follows from the i.i.d. assumption in (A3), particularly, £ ∗ ¤ Ewt [Ewt+1 Jt+2 (wt , wt+1 ) ] ·Z ¸ ∗ =Ewt Jt+2 (wt , wt+1 )fwt+1 (wt+1 )dwt+1 Z Z ∗ = Jt+2 (wt , wt+1 )fwt+1 (wt+1 )dwt+1 fwt (wt )dwt Z Z ∗ = Jt+2 (wt , wt+1 )fwt ,wt+1 (wt , wt+1 )dwt+1 dwt ∗ ∗ =Ewt ,wt+1 [Jt+2 (wt , wt+1 )] = Et [Jt+2 (wt , wt+1 )]
where fwt (·), fwt+1 (·) and fwt ,wt+1 (·, ·) are density functions of wt , wt+1 and their joint density function respectively and fwt ,wt+1 (·, ·) = fwt (·)fwt+1 (·) follows from assumption (A3). Summing (43) and (47) leads to ∗ Jt∗ ≥ kd∗0|t k2Ψ + Et [kd∗0|t+1 (wt )k2Ψ ] + Et [Jt+2 (wt , wt+1 )]
Repeating the above procedure infinite times leads to ∞ > Jt∗ ≥ kd∗0|t k2Ψ +
∞ X
h i ∗ Et kd∗0|i (wt , · · · , wi−1 )k2Ψ + lim Et [Jt+r (wt , · · · , wt+r−1 )] r→∞
i=t+1
∗ (w , · · · , w where the left inequality follows from Theorem 7(ii). Using the fact that limr→∞ Et [Jt+r t t+r−1 )] > ∗ 2 0 and kd0|t kΨ is finite, we have
∞>
∞ X
h i Et kd∗0|i (wt , · · · , wi−1 )k2Ψ
i=t+1
By applying Markov bound (given non-negative random variable η and any ² ≥ 0, E[η] ≥ ²Pr{η ≥ ²}) and considering kd∗0|i k2Ψ as a random number, we have ∞>²
∞ X
Pr(kd∗0|i k2Ψ ≥ ²)
(48)
i=t+1
for any arbitrary small ² > 0. From the First Borel-Cantelli Lemma [26], this implies that limi→∞ Pr(kd∗0|i k2Ψ ≥ ²) = 0. Hence d0|i → 0 with probability one as t increases. Consequently, the MPC control law (22) converges to Kf xt with probability one. When this happens, the closed-loop system converges to xt+1 = Φxt +wt and, hence, xt converges to F∞ (Kf ) with probability one. (iii) follows directly from (ii) and assumption (A4) that F∞ (Kf ) ⊂ int(Xf ). 20
National University of Singapore
E
M.E. Dept.
Proof of Theorem 4
Proof. (i) The replacement of cost function J(d, D) by V (d, D) does not affect the feasibility of problem PN (x). This means that part (i) of Theorem 3 remains valid. (ii) Let Vt∗ and Vˆt+1 be defined in the same manner as Jt∗ and Jˆt+1 in the statement of proofs of Theorem 3. Following the same reasoning as in (40), it can be shown that Vt∗ − Vˆt+1 (wt ) = kd∗0|t k2Ψ + p(wt )
(49)
where p(wt ) =
N −1 X
i ∗ 2 i ∗ i ∗ i ∗ (γ1 kvec(Di|t ) k + γ2 kvec(Di|t ) k − 2(d∗i|t )T Ψ(Di|t ) wt − k(Di|t ) wt k2Ψ ).
(50)
i=1
Hence p(wt ) ≥ ≥ = =
N −1 X
i ∗ 2 i ∗ i ∗ i ∗ 2 (γ1 kvec(Di|t ) k + γ2 kvec(Di|t ) k − 2kd∗i|t kkΨkkwt kk(Di|t ) k − kΨkkwt k2 k(Di|t ) k )
i=1 N −1 X
i ∗ 2 i ∗ i ∗ i ∗ 2 (γ1 kvec(Di|t ) k + γ2 kvec(Di|t ) k − 2αβkΨkk(Di|t ) kF − α2 kΨkk(Di|t ) kF ) (51)
i=1 N −1 X
i ∗ 2 i ∗ i ∗ i ∗ 2 (γ1 kvec(Di|t ) k + γ2 kvec(Di|t ) k − 2αβkΨkkvec(Di|t ) k − α2 kΨkkvec(Di|t ) k )
i=1 N −1 X
i ∗ 2 i ∗ ((γ1 − α2 kΨk)kvec(Di|t ) k + (γ2 − 2αβkΨk)kvec(Di|t ) k) ≥ 0
i=1 i )∗ k ≥ k(D i )∗ k and k(D i )∗ k = kvec(D i )∗ k are used. Hence, p(w ) ≥ where the facts k(Di|t t F F i|t i|t i|t 0 under (29). As a consequence, equation (49) implies ∗ Vt∗ − kd∗0|t k2Ψ ≥ Vt+1 ≥0
(52)
Hence, {Vt∗ } is a monotonic non-increasing sequence and is bounded from below by zero. This means that V∞ := limt→∞ Vt∗ ≥ 0 exists. Repeating (52) for t from 0 to ∞ and summing them up, it follows that ∞ X ∗ ∞ > V0 − V∞ ≥ kd∗0|t k2Ψ (53) t=0 limt→∞ d∗0|t
Since Ψ is positive definite, this implies that = 0 and limt→∞ ut = Kf xt . Therefore, the stated result follows. (iii) follows (ii) and assumption (A4) that F∞ (Kf ) ⊂ int(Xf ).
F
Computation of β
β := max(x,d,D)∈TN ,i∈ZN −1 kdi k = max(x,d,D)∈TN kd0 k is due to the fact that for any (x, d, D) ∈ ¯ D) ¯ ∈ TN can be found such that d¯0 = di . SpecifiTN and integer i ∈ Z+ x, d, N −1 , a set of (¯ cally, given (x, d, D) ∈ TN and let the correspondingly defined state and control sequence be 21
National University of Singapore
M.E. Dept.
{x0 , . . . , xN } and {u0 , . . . , uN −1 }. According to (16) xN ∈ Xf for all possible disturbances. ¯ D) ¯ can be defined by Then for any i ∈ Z+ x, d, N −1 , (¯ i
x ¯ = Φ x+
i−1 X j=0
i−1−j
Φ
( dj+i Bdj , d¯j = 0
j ∈ ZN −1−i ¯ jk = ,D N −i≤j ≤N −1
( k Dj+i 0
j ∈ Z+ N −1−i k ∈ Z+ j . N −i≤j ≤N −1
¯ D) ¯ define the control sequence where x ¯ is the nominal state of xi defined by (x, d, D) and (d, {ui , . . . , uN −1 , Kf xN , . . . , Kf xN −1+i }. According to (A4) under controller ut = Kf xt all the constraints are ¯ D) ¯ satisfies (12)-(16), namely satisfied and xt ∈ Xf for t ≥ N since xN ∈ Xf . Therefore, (¯ x, d, ¯ ¯ (¯ x, d, D) ∈ TN . As a result, max(x,d,D)∈TN kd0 k ≥ max(x,d,D)∈TN kdi k, for any i ∈ ZN −1 and β = max(x,d,D)∈TN kd0 k.
22
