Robust Stable Economic MPC with Applications in Engine Control
53rd IEEE Conference on Decision and Control December 15-17, 2014. Los Angeles, California, USA
Robust Stable Economic MPC with Applications in Engine Control Timothy Broomhead1 , Chris Manzie1 , Rohan Shekhar1 , Michael Brear1 and Peter Hield2
Abstract— Economic Model Predictive Controllers have shown to improve a plant’s economic performance using state dependant economic stage costs. Recent extensions have provided continual feasibility guarantees, despite changes in economic parameters, however, perfect plant models have been assumed. This assumption is invalid in practise due to modelling errors or un-modelled disturbances and can therefore lead to infeasibility of the optimisation problem. This paper proposes a robust economic model predictive controller, which takes advantage of constraint tightening techniques to guarantee feasibility despite modelling errors. Inputto-state stability is proven using a Lyapunov function. The advantages of this method are highlighted against alternative control structures in the application of power tracking for diesel engines in series hybrid type applications.
I. I NTRODUCTION To minimise a plant’s economic cost, Model Predictive Control (MPC) is often employed as part of a hierarchical control architecture, where a high level, static Real Time Optimiser (RTO), chooses operating set-points which minimise an economic cost [1]. The lower level Tracking Model Predictive Control (TMPC) then steers the plant to the chosen set-points in an admissible way, taking into account plant dynamics and constraints. Typically, the TMPC problem is formulated using an error system, with convergence of the closed-loop system guaranteeing set point tracking. Demand for higher efficiencies and reduced economic costs is driving development of alternative approaches. The TMPC is designed to reject disturbances and track the system as closely as possible to the given set-points, however it does not necessarily take into account the economic cost of the state trajectories. Therefore, if the set-points change frequently with respect to the plant dynamics, this approach may result in poor economic performance. An alternative is to combine the RTO and TMPC into a single controller, known as Economic Model Predictive Control (EMPC) [2]. In [2] a Lyapunov function is used to show stability of the EMPC formulation under the assumption of strong duality with the steady state problem and using a terminal constraint. The assumption of strong duality has been relaxed to a dissipativity assumption in [3]. The need for terminal constraints can be relaxed for sufficiently long prediction horizons where controllability and turnpike properties hold [4], however in some applications the required horizon length for acceptable performance may be computationally prohibitive.
Economic parameters may change as a result of variation in material prices or requested plant production rates. With a change in economic parameters, the optimal steady state condition may shift and potentially result in a loss of feasibility for controllers with terminal constraints as dependant on this state. Recent work in [5] extends the results of [2] relaxing the terminal constraint to be any point on the steady state manifold and including a terminal weight. Referred to as Economic Tracking Model Predictive Control (ETMPC) in this paper, this approach requires modification of the stage cost to ensure convergence to the optimal steady state, resulting in sub-optimal transient behaviour, though feasibility is guaranteed for changes in economic parameters. Each of the approaches introduced so far assumes a perfect prediction model to ensure continual feasibility. In practice, this assumption is invalid due to external disturbances, measurement error and modelling errors, potentially resulting in infeasibility through constraint violation. Many existing robustness techniques, such as open or closed loop Min-Max MPC, are overly conservative or computationally intractable [6]. Constraint Tightening (CT) [7] is a robustness technique for which robust feasibility can be guaranteed. In the general case, computational requirements for CT do not exceed that of nominal robustness. The basic principle is that constraints are systematically tightened along the horizon to reserve a margin which can then be used to correct for errors resulting from model mismatch at some future time. This paper combines the results of [5] with the constraint tightening principles in [7] to produce a practical control structure capable of directly addressing economic objectives and constraints whilst being robust to modelling errors. In the application of diesel generator control, the economic objective is to minimise the penalty associated with not meeting the requested output power, whilst minimising the cost associated with fuel consumption under emissions and safety constraints. While previous approaches utilising MPC for engine control often represent TMPC formulations [8], [9], [10], an economic formulation can increase overall performance. Representation of diesel engines requires complex models [11], which can be impractical for control design. Model reduction techniques are therefore required, resulting in modelling error, motivating the development of a Constraint Tightened Tracking Economic Model Predictive Control (CT-TEMPC).
The authors would like to thank the Australian Defence Science and Technology Organisation, the Defence Science Institute and the Elizabeth and Vernon Puzey trust for their support of this research. 1 Department of Mechanical Engineering, The University of Melbourne, Australia. [email protected] 2 Defence Science and Technology Organisation, Melbourne, Australia
978-1-4673-6090-6/14/$31.00 ©2014 IEEE
2511
II. P ROBLEM P RELIMINARIES Consider a linear-time invariant system xk+1 = Axk + Buk + wk
(1)
where x ∈ Rn is the state, u ∈ Rm is the control input and w ∈ Rn is an unknown but bounded disturbance, obeying wk ∈ W. The model outputs y ∈ Rl are given by yk = Cxk + Duk
(2)
the optimal trajectory. The following assumption is therefore made in order to demonstrate asymptotic stability. Assumption 5 (Strong Duality of Steady-State problem): Let ˜l(x, u) be the rotated stage cost function given by ˜l(x, u, p) = l(x, u, p) + λT (x − Ax − Bu) − l(x∗ , u∗ , p) s s
and constraints must be satisfied for all k, described by yk ∈ Y
(3)
where Y is a compact and convex set. Disturbances on the output (2) could be incorporated as in [10], however are omitted due to space limitations. Assumption 1: W is compact, convex and contains 0. Assumption 2: There exists a sequence of nilpotent controller gains Kj which drive the system xj+1 = (A + BKj )xj , j ∈ I0:Nnp to the origin in Nnp < N − 1 steps, where Ia:b denotes the set of integers from a to b. The tightened constraint sets for robustness are defined as Y0 = Y, Yj+1 = Yj ∼ [C + DKj ]Lj W, ∀j ∈ I1:N
(4)
where λ is a multiplier ensuring a unique minimum at (x∗s , u∗s ) ∀x, u given Cx + Du ∈ Y s . Then there exists a κ∞ -function α(.) such that ˜l(x, u, p) ≥ α(|x − xs |). III. C ONTROL F ORMULATION Definition 1: For the a given constant γ > |λ| a convex offset function VO (x, u) satisfies VO (x, u) − VO (xs , us ) ≥ γ|x − xs | and for which min VO (x, u)
(x,u)∈Zs
has a unique solution at (x∗s , u∗s ). The economic cost function for changing parameters is [5] lt (x, u, p) = l(x + x∗s , u + u∗s , p)
where ”∼” denotes the Pontryagin difference [12] and Lj is L0 = I, Lj+1 = (A + BKj )Lj , ∀j ∈ I1:N .
VN (x, p; u) =
(5)
(6b)
s = YN ∼ B and B is an arbitrarily small ball of where YN radius > 0, required for stability guarantees. We wish to generate a robustly feasible feedback control law which aims to minimise cost P∞ (7) V∞ = k=0 l(xk , uk , p)
where economic parameters p may vary in time. The best robustly achievable steady state are (x∗s (p), u∗s (p)) = arg min l(x, u, p) s.t (x, u) ∈ Zs . x,u
(10a)
lt (xk+j − xk+N −1 , uk+j − uk+N −1 , p)
+ VO (xk+N −1 , uk+N −1 )
(8)
Henceforth, dependence of x∗s , u∗s on p is removed for brevity. Assumption 4: The stage cost function l(x, u) is Lipschitz continuous in (X × U) with Lipschitz constant Ml , i.e.: |l(x1 , u1 , p) − l(x2 , u2 , p)| ≤ Ml (|x1 − x2 | + |u1 − u2 |) ∀x1 , x2 ∈ X , ∀u1 , u2 ∈ U. In EMPC, standard Lyapunov arguments cannot be used since the economic cost is not necessarily decreasing along
(10b)
The optimisation problem PN (x, p) is given as PN (x, p) : arg min VN (x, p; u) u
s.t xk = x
(6a)
ys = Cxs + Dus } Xs = {xs : ∃u(x, u) ∈ Zs }
N −1 X j=0
Remark 1: Since the controller gains are nilpotent, Lj = 0 ∀j ≥ Nnp . If a disturbance feedback controller is utilised for disturbance rejection, the controller gains may be optimised as part of a convex optimisation [13]. Assumption 3: The constraint set (3) is sufficiently large with respect to the disturbance set W and Kj are chosen such that the resulting tightened sets are non-empty. Define the set of robustly admissible equilibrium conditions and states for the nominal system (where wk = 0) as Zs = {(xs , us ) : ys ∈ Y s , xs = Axs + Bus ,
(9)
(11a) (11b)
xk+N ∈ Xs
(11c)
xk+j+1 = Axk+j + Buk+j
(11d)
yk+j = Cxk+j + Duk+j
(11e)
yk+j ∈ Yj , ∀j ∈ I0:N −1
(11f)
while the robust economic MPC control law κeN (x, p) is implicitly defined as the optimal solution to PN (x, p). IV. S TABILITY A NALYSIS A. Recursive Feasibility Theorem 1 (Feasibility): Consider a system described by (1) and (2), subject to constraints (3) and let Assumptions 1, 2 and 3 hold. Under the control law given by PN (x, p) and given that the system is initially within the feasibility set XN , the system will remain within XN for all time. Proof: Consider any xk0 ∈ XN . Define the optimal solution to PN (xk0 , p) as u0k0 +j , with corresponding states x0k0 +j and outputs yk00 +j , ∀j ∈ I0:N −1 . Consider the candidate solution for k0 + 1 u ˆk0 +j = u0k0 +j + Kj−1 Lj−1 wk0 ∀j ∈ I1:N −1 u ˆk0 +N =
u0k0 +N −1
(12a) (12b)
This solution is created by shifting the previous solution by one step, adding perturbations for disturbance rejection and appending one step of control to the end. Due to the steady state terminal condition, the appended control step can be the final control action of the previous solution. The initial state
2512
at k0 + 1 is given by the system dynamics (1). Substitution into the internal model (11d) gives xk0 +1 = x0k0 +1 + wk0 . Substituting this initial condition and the candidate solution (12) gives the candidate state sequences ∀j ∈ I1:N x ˆk0 +j =
x0k0 +j
+ Lj−1 wk0
x ˆk0 +N +1 = Aˆ xk0 +N + B u ˆk0 +N
Proof: Using the Cauchy-Schwarz’s inequality, λT (x − ≤ |λ||x − x∗s |. Let γˆ = γ − |λ|, then V˜O (x, u) = VO (x, u) + λT (x − x∗s ) − VO (x∗s , u∗s )
x∗s )
≥ γ|x − x∗s | + λT (x − x∗s )
(13a) (13b)
and output sequences ∀j ∈ I1:N −1 yˆk0 +j = yk00 +j + (C + DKj−1 )Lj−1 wk0
(14a)
yˆk0 +N = C x ˆk0 +N + Dˆ uk0 +N
(14b)
Initial condition, state equality and output equality constraints are therefore satisfied by construction. The satisfaction of the terminal constraint (11c) uses that for a nilpotent controller Lj = 0 for j ≥ Nnp . Therefore, since Nnp < N − 1, (13a) implies x ˆk0 +N = x0k0 +N . Feasibility 0 at k0 implies xk0 +N is a steady state with input u0k0 +N −1 , substitution into (13b) gives x ˆk0 +N +1 = x0k0 +N , therefore x ˆk0 +N +1 = x ˆk0 +N and the terminal constraint is satisfied. Feasibility at k0 requires yk00 +j ∈ Yj ∀j ∈ I1:N −1 . Since Yj−1 is related to Yj by (4), (14a) implies yˆk0 +j ∈ Yj−1 ∀j ∈ I1:N −1 . This satisfies the output constraints for all but the final step N . Since x ˆk0 +N = x0k0 +N −1 0 and u ˆk0 +N |k0 +1 = uk0 +N −1|k0 then yˆk0 +N = yk00 +N −1 . Since at time k0 feasibility requires yk00 +N −1 ∈ YN −1 then yk00 +N ∈ YN −1 ⊆ YN . This demonstrates the candidate solution satisfies all constraints, and hence for any x ∈ XN a feasible solution exists. B. Input-to-state stability
= γˆ |x − x∗s | + |λ||x∗s − x| − λT (x∗s − x) ≥ γˆ |x − x∗s |. Lemma 3: Consider a system described by (1) and (2), subject to constraints (3). Let Assumptions 1 - 5 hold. Under the control law given by (11) and given x ∈ XN , there exists a κ∞ -function, α(.), and κ-function, ρ(.), such that V˜ 0 (x+ , p) − V˜ 0 (x, p) ≤ −α(|x − x0 |) + ρ(|w|) (20) N
Proof: Consider at k0 any xk0 ∈ XN and recall the candidate solution (12) for the problem at k0 + 1. Compare the optimal cost V˜N0 (xk0 ) with that of the candidate solution V˜N (ˆ xk0 +1 ) − V˜N0 (xk0 ) = PN −1 ˜ 0 − j=0 lt (xk0 +j − x0k0 +N −1 , u0k0 +j − u0k0 +N −1 , p) PN ˜ + lt (ˆ xk +j − x ˆ0 ,u ˆk +j − u ˆ0 , p) j=1
which has the following properties l˜t (0, 0, p) = ˜l(x∗s , u∗s , p) = 0 l˜t (z, v, p) ≥ α1 (|z|) for the κ∞ -function α1 (.).
(16a) (16b)
A rotated offset function is used, defined as V˜O (x, u) = VO (x, u) + λT (x − x∗s ) − VO (x∗s , u∗s ),
(17)
with λ from Assumption 5. The auxiliary cost function is V˜N (x, p; u) =
N −1 X
l˜t (xk+j − xk+N −1 , uk+j − uk+N −1 , p)
j=0
+ V˜O (xk+N −1 , uk+N −1 ),
(18)
and the auxiliary optimisation problem is P˜N (x, p) : arg min V˜N (x, p; u) s.t (11b) − (11f ). (19a) u
Lemma 1 (Lemma 1 of [5]): The problem P˜N (x, p) gives the same optimal control sequence as PN (x, p). Lemma 2: Due to Definition 1 there exists a positive γˆ such that V˜O (x, u) ≥ γˆ |x − x∗s |
0
k0 +N
Since Nnp < N − 1 and by the definition of u ˆ, we know that x0k0 +N −1 = x ˆk0 +N and u0k0 +N −1 = u ˆk0 +N so the offset functions equate. Furthermore, from Assumptions 4 and 5 and (15) we can establish a bound on the Lipschitz ˜ l ≤ Ml + |λT (I − A)| + |λT B|). constant for l˜t (x, u, p) as M Additionally, taking into account that the disturbance is nilpotent for j ≥ Nnp and property (16a) we find VˆN (xk +1 ) − V 0 (xk ) 0
(15)
k0 +N
0
−VO (x0k0 +N −1 , u0k0 +N −1 ) + VO (ˆ xk0 +N , u ˆk0 +N ).
N
0
≤ − l˜t (x0k0 − x0k0 +N −1 , u0k0 − u0k0 +N −1 , p) ˜ l PNnp (|Lj−1 | + |Kj−1 Lj−1 |)|wk |. +M
For the purposes of showing stability, an auxiliary problem is introduced, with stage cost l˜t (x, u, p) = ˜l(x + x∗s , u + u∗s , p)
N −1
N
j=1
0
xk0 +1 ) and due to property By optimality V˜N0 (xk0 +1 ) ≤ V˜N (ˆ (16b) we can state V˜N (x+ ) − V˜N0 (x) ≤ −α1 (x − xN −1 ) + ˜ l PNnp (|Lj−1 | + |Kj−1 Lj−1 |)ν. ρ(|w|), where ρ(ν) = M j=1 Define the constant PNnp γ0 = Ml j=0 (|(1 + Kj )Lj |)|. (21) Lemma 4: Consider a system described by (1) and (2), subject to constraints (3). Let Assumptions 1 - 5 hold. Let (ˆ xs , u ˆs ) ∈ Zs and u ˆ be a sequence of control actions u ˆj = Kj (xj − x ˆs ) + u ˆs , where xj+1 = Axj + B u ˆj . Then under Assumption 2, (ˆ xs , u ˆs ) can be chosen such that u ˆ is a feasible solution to P˜N (x) and such that |x − x ˆs | ≤ ω(), ω > 0 given |xj − x∗s | ≤ ∀j ∈ I0:N −1 . It then follows that ˆ ) ≤ γ0 |x − x V˜N (x, p; u ˆs | + VO (ˆ xs , u ˆs ) + d, where d = λT (x − xs ) − VO (x∗s , u∗s ), and γ0 is as defined in (21). Proof: Noting from (5) that uj = Kj Lj (x − x∗s ) + u∗s , xj = Lj (x − x∗s ) + x∗s , and Lj = 0, ∀j > Nnp , then PN −1 ∗ ∗ j=0 (l(xj , uj , p) − l(xs , us , p)) PNnp ≤ Ml j=0 (|Lj (x − x∗s )| + |Kj Lj (x − x∗s )|) PNnp ≤ Ml j=0 (|(1 + Kj )Lj |)|x − x∗s | = γ0 |x − x∗s |.
2513
It can then be shown that P −1 ˆ) = N V˜N (x, p; u ˆs + x∗s , uj − u ˆs + u∗s , p) j=0 (l(xj − x − l(x∗s , u∗s , p)) + VO (ˆ xs , u ˆs ) + d ≤ γ0 |x − x ˆs | + VO (ˆ xs , u ˆs ) + d, where d = λ (x − x∗s ) − VO (x∗s , u∗s ). Lemma 5: Consider a system described by (1) and (2), subject to constraints (3). Let Assumptions 1 - 5 hold. Consider the constant γ > γ0 as defined in (21). Assume for an initial state x the optimal solution to P˜N (x, p) is (x0 ,u0 ). If x0N −1 = x and u0N −1 = κeN (x, p), then T
x0N −1 = x∗s , u0N −1 = u∗s
(22)
Proof: Consider that the optimal solution to P˜N (x, p) is (x0 , u0 ), x = x0N −1 and κeN (x, p) = u0N −1 . Then (x, u) ∈ Zs and the optimal cost function is V˜N0 (x, p) = V˜O (x0N −1 , u0N −1 ). In order to prove the lemma by contradiction, first consider that (x0N −1 , u0N −1 ) 6= (x∗s , u∗s ). Define (ˆ xs , u ˆs ) = β(x0N −1 , u0N −1 ) + (1 − β)(x∗s , u∗s ), and note that due to convexity of Zs , (ˆ xs , u ˆs ) ∈ Zs . Following from Assumption 2, there exists a βˆ ∈ (0, 1) such that for any ˆ 1), the control law uj = Kj (xj − x β ∈ [β, ˆs ) + u ˆs drives ˆs in an admissible way in ≤ N the system from x0N −1 to x steps. Define u ˆ as the control law uj = Kj (xj − x ˆs ) + u ˆs , u ˆ is a feasible solution to P˜N (x0N −1 , p). From Lemma 4 we have that V˜N0 (x0N −1 , p) = V˜O (x0N −1 , u0N −1 ) ≤ V˜N (x0N −1 , p; u ˆ) ≤ γ0 (1 − β)|x0N −1 − xs | + VO (ˆ xs , u ˆs ) + d where d = λT (x0N −1 − xs ) − VO (xs , us ). Define W (x0N −1 , β) = γ0 (1 − β)|x0N −1 − x∗s | + VO (ˆ xs , u ˆs ) + d. The partial derivative of W about β is ∂W /∂β =
−γ0 |x0N −1
− xs | + g
0
(x0N −1
−
x∗s , u0N −1
−
u∗s )
Theorem 2 (Input-to-state stability): Consider a system described by (1) and (2), subject to constraints (3). Let Assumptions 1 - 5 hold. Under the control law given by (11) and given x ∈ XN , the closed loop system will be input-to-state stable, with a domain of attraction XN for all γ > γ0 . Proof: Due to property (16b) and Lemma 2 we have that V˜N0 (x, p) ≥ α1 (|x − xN −1 |) + γˆ |xN −1 − x∗s |, for all x ∈ XN . Hence there is a κ∞ -function, αJ (.), such that V˜N0 (x, p) ≥ αJ (|x − xN −1 | + |xN −1 − xs |) ≥ αJ (|x − x∗s |). Since V˜N0 (x∗s ) = 0 and V˜N0 (x) is continuous in XN , there exists a κ∞ -function βJ such that for all x ∈ XN [14, Proposition 2.18] V˜N0 (x, p) ≤ βJ (|x − x∗s |). Let e(x) = x − x0N −1 for all x ∈ XN and x0N −1 ∈ Xs . Then from Lemma 3, V˜N0 (x+ , p) − V˜N0 (x, p) >= −α1 (|e(x)|) + ρ1 (|w|), where ρ1 is positive definite. If e(x) = 0, then x = x0N −1 , and from Lemma 5 this requires that x0N −1 = x∗s ; while if x = x∗s , then by optimality x0N −1 = x∗s and hence e(x) = 0. Hence we can conclude that |e(x)| = 0 iff x = x∗s . Since XN is compact, |e(x)| = 0 when x = x∗s and otherwise |e(x)| > 0, there exists an κ-function, αe (.), such that |e(x)| ≥ αe (|x − x∗s |). It then follows that α1 (|e(x)|) ≥ α1 (αe (|x − x∗s |)) = αV (|x − x∗s |) where αV (.) is a κ∞ -function. It can now be concluded that: 1) αJ (|x − xs |) ≤ V˜N0 (x, p) ≤ βJ (|x − xs |), ∀x ∈ XN 2) V˜N0 (x+ , p) − V˜N0 (x, p) ≤ −αV (|x − xs |) + ρ1 (|w|), ∀x ∈ XN , ∀w ∈ W Therefore V˜N0 (x, p) is an ISS-Lyapunov function in XN , and xs is a local input-to-state (ISS) stable equilibrium point for the closed loop system controlled by P˜N (x) [15]. As a result of Lemma 1, ISS stability also extends to PN (x). V. A PPLIED E XAMPLE A. Problem Formulation
where g 0 is a sub-differential of VO (ˆ xs , u ˆs ). Evaluating this derivative for β = 1 ∂W /∂β|β=1 = −γ0 |x0N −1 −x∗s |+¯ g 0 (x0N −1 −x∗s , u0N −1 −u∗s ) where g¯0 is a sub-differential of VO (x0N −1 , u0N −1 ). Since VO is convex, and due to Definition 1 g¯0 (x0N −1 − x∗s , u0N −1 − u∗s ) ≥ VO (x0N −1 , u0N −1 ) − VO (x∗s , u∗s ) ≥ γ|x0N −1 − x∗s | Therefore ∂W /∂β|β=1 ≥ −γ0 |x0N −1 − x∗s | + γ|x0N −1 − x∗s | = (γ − γ0 )|x0N −1 − x∗s | Since γ > γ0 and x0N −1 6= x∗s , ∂W /∂β|β=1 > 0. This ˆ 1) such that the cost would mean that there exists a β ∈ [β, to move the system from x0N −1 to x ˆs , is smaller than to remain in x0N −1 . This contradicts optimality of the solution to P˜N (x, p) and hence x0N −1 = x∗s .
In this example, the developed robust economic MPC formulation will be utilised in a diesel generator power tracking application. The economic objective is associated with minimising the error between the requested power and the brake power produced by the generator, with a secondary economic objective to minimise fuel consumption. To simulate the plant, a calibrated and validated high order Mean Value Engine Model (MVEM) is utilised, similar to that described in [10]. The model inputs u := [tδ τload uvgt ]T are the injection duration, load torque and Variable Geometry Turbocharger (VGT) demand. The system states x := [ωe pim pem ωt xvgt ]T are engine speed, intake and exhaust manifold pressure, turbocharger speed and VGT position. Since the MVEM is highly nonlinear, and computationally impractical for use in an online control strategy, model reduction techniques described in [16] are utilised. The relatively fast manifold pressure dynamics are eliminated before the reduced model is linearised around a fixed operating point. Since the controller will be tested against the full order MVEM, the model reduction results in modelling error,
2514
l(x, u, p) = (ωe τeng − Pref )2 + (Ψ · m ˙ f (ωe , tδ ))2
10
(23)
where the requested power is the economic parameter Pref and m ˙ f (ωe , tδ ) is the fuel consumption map with weighting Ψ, assumed small in this example. In the sequel, a number of different control strategies will be examined, and their ability to minimise (7) compared. Where relevant, the steady state conditions (xs , us ) used in controller formulations are found for each Pref using (8). The reference Pref is constructed with step changes every 5 seconds of a random size, to represent an aggressive test-case. The power request is limited between 20 kW and 57 kW with a maximum absolute step size of 20 kW. A timestep of 100 ms has been chosen, based on the bandwidth of the engine’s sensors and actuators. Denote the controller as CT-TEMPC, constructed using optimisation problem (11) and the offset function VO = kQV (xk+N − xs )k1 + kRV (uk+N − us )k1 ,
(24)
where QV , RV are chosen for good transient performance. Tightened constraint sets are calculated using (4). A Monte Carlo approach is used to identify the disturbance set W, capturing the error between the reduced and full order models. The controller CT-TEMPC, initially with a conservative estimate for W, is requested to track the randomly generated reference signal, Pref , for a large number of samples. The modelling errors between the reduced and full order models are captured at each time-step. To reduce conservatism, lower and upper error bounds were chosen at the 0.1 and 99.9 percentiles to define the set W. This process was iterated upon several times to generate the final error set, shown
−2
−1 rpm
0
Turbine Speed
1
−4
−2
0 2 3 rpm×10
−0.1
0 %
2 0
25
10
15
20
25
N
B. Comparison of economic control formulations One approach to economic MPC is to use the pure economic cost as the stage cost without terminal constraints, as described in [4]. This eliminates the need to determine optimal steady state points explicitly, potentially leading to a good approximation of (7) assuming an accurate internal model, and sufficiently long horizon. Denote the robust version of this controller as CT-EMPC, constructed using the optimisation PN −1 arg minu j=0 l(xj , uj , p) s.t (11b), (11d) − (11f ). (25) In Fig. 2 this controller is demonstrated with a power request at 20 kW, changing to 57 kW at 10 seconds. In this application, a horizon length less than 16 results in unacceptable steady state tracking error. Furthermore, with a horizon length of 16, the nonlinear optimisation is computationally intractable for online control of diesel engines. A horizon length of around 10 is more realistic in terms of computational load, motivating alternative controller formulations. One option is to use a robustified version of the EMPC in [2], constructed using optimisation (25) with an additional terminal constraint xk+N = x∗s , denoted as CT-EMPC-D. However, it was found that the controller became infeasible after changes in power request with a reasonable horizon length of 10, due to a small region of attraction. As shown in [5], a sufficiently large weighting on the offset function results in the output of (11) being equal to CT-EMPC-D Ref
CT−TEMPC−H
CT−TEMPC
TEMPC
50 40 30 0
10
20
30
40
50
60
10
20
30
40
50
60
10
20
30 Time (s)
40
50
60
10 0 −10 0
0.1
Fig. 1. Prediction errors for each state. Shown are 0.1 and 99.9 percentiles
4
in Fig. 1. Due to the Monte Carlo approach and the error bounds used, absolute guarantees can no longer be made for the controller; however good practical robustness remains, as demonstrated later.
VGT Position
4
15 20 Time (s)
6
Fig. 2. Left: Closed-loop output of CT-EMPC after a step-change in power request. Right: Steady state tracking error.
Engine Speed 3 (rpm) × 10
Engine Speed
N=10 N=13 N=16 Ref
40
20
which will be accounted for through constraint tightening. To represent measurement noise, a small normally distributed disturbance is added to the simulated states, which can be treated as an additive disturbance [7]. The engine speed, turbine speed and each input has lower and upper operating bounds, while VGT actuator also has a maximum slew rate constraint, summarised in Table I. The normalised air-to-fuel ratio of the engine, λ, has a lower limit in order to prevent excessive smoke production, set to 1.4 in this application. This non-linear, non-convex constraint is approximated using a convex polytopic constraint. The power tracking economic objective can be summarised with the following stage cost
Final Tracking Error (kW)
Power Output (kW)
∆uvgt (%) −5 5
Power Output (kW)
Lower Upper
60
Power Error (kW)
Bounds
TABLE I L OWER AND U PPER C ONSTRAINTS ωe ωt uvgt tδ τload (rpm) (rpm) (%) (ms) (Nm) 1500 50 × 103 65 0.35 0 2500 140 × 103 85 1 400
2
1.5 0
Fig. 3. Comparison of Economic MPC formulations. The TEMPC controller fails after 3 seconds, indicated on the Engine Speed plot.
2515
when inside its region of attraction. Denote this controller as CT-TEMPC-H, where H indicates that QV is chosen to be large. In Fig. 3, CT-TEMPC-H is compared against CT-TEMPC and its non-robust equivalent, TEMPC, both with the same offset function (24). All controllers use a horizon length of N = 10. The non-robust controller, TEMPC, fails due to an engine speed constraint violation, indicated on Fig. 3 at around 3 seconds. Due to the finite horizon lengths, offset functions, economic stage cost for tracking and model mismatch, each of the controllers presented are suboptimal; this is the price paid for guaranteed feasibility and convergence. Selection of the offset function directly influences the closed loop performance, as is made clear in the comparison between CT-TEMPC-H and CT-TEMPC. During large step changes, CT-TEMPC-H increases tracking error in order to ensure fast convergence of the terminal state to x∗s . This turns out to be undesirable behaviour, which can be removed with a different offset function, as seen in the performance of CT-TEMPC.
It can be seen that the CT-TEMPC controller out-performs the tracking controller during the transients, reaching the target power output in a shorter period of time, this is particularly noticeable when power demand is reduced. Despite the apparent similarity of the trajectories, the cumulative economic cost (23) of the economic controller over the simulation was 28 % less than that of the tracking controller. Additionally, the time-consuming tuning process for the tracking controllers must be performed for any changes in economic parameters, further motivating the use of the economic controller. VI. C ONCLUSIONS Through the use of constraint tightening, a robust economic MPC has been presented. The controller is guaranteed to remain feasible despite modelling errors and changes in economic parameters. An ISS result is given based on a Lyapunov function, ensuring convergence to a set around the optimal steady state. The efficacy of the developed controller has been demonstrated against alternative formulations in the application of diesel generator power tracking.
C. Tracking controller comparison
R EFERENCES
The developed economic controller CT-TEMPC will now be compared against a robustified TMPC formulation [17], representing the hierarchical set-point tracking approach. The robust tracking controller, denoted CT-TMPC, is described by the optimisation PN −1 arg minu j=0 |xk+j − xk+N |2Q + |uk+j − uk+N |2R + VO s.t (11b) − (11f ),
(26)
where VO is from (24), with the same tuning as CT-TEMPC. The weights Q and R have been tuned using Matlab’s Optimisation Toolbox to minimise the overall economic cost over a 150 second simulation using a randomly generated reference power trajectory. The optimisation used the internal model as the plant model, representing an off-line tuning process, and a Cholesky decomposition of Q and R to enforce positive definiteness. In Fig. 4, CT-TMPC is compared against CT-TEMPC. Both controllers exhibit negligible steady state tracking error, however their closed-loop behaviour differs during transients.
Power Output (kW)
Ref
CT−TMPC
CT−TEMPC
50 40 30 0
10
20
30
40
50
60
10
20
30
40
50
60
10
20
30 Time (s)
40
50
60
Power Error (kW)
10 0
Engine Speed 3 (rpm) × 10
−10 0 2
1.5 0
Fig. 4.
Comparison of tracking and economic MPC formulations
[1] S. J. Qin and T. A. Badgwell, “A survey of industrial model predictive control technology,” Contr. Eng. Prac., vol. 11, pp. 733–764, 2003. [2] M. Diehl, R. Amrit, and J. B. Rawlings, “A lyapunov function for economic optimizing model predictive control,” IEEE Trans. on Auto. Contr., vol. 56, no. 3, pp. 703–707, 2011. [3] D. Angeli, R. Amrit, and J. B. Rawlings, “On average performance and stability of economic model predictive control,” IEEE Trans. on Auto. Contr., vol. 57, no. 7, pp. 1615–1626, 2012. [4] L. Gr¨une, “Economic receding horizon control without terminal constraints,” Automatica, vol. 49, no. 3, pp. 725–734, 2013. [5] A. Ferramosca, D. Limon, and E. Camacho, “Economic MPC for a changing economic criterion for linear systems,” Automatic Control, IEEE Transactions on, vol. In press, 2014. [6] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. Scokaert, “Constrained model predictive control: Stability and optimality,” Automatica, vol. 36, no. 6, pp. 789–814, 2000. [7] A. Richards and J. How, “Robust model predictive control with imperfect information,” in Proceedings of the ACC, 2005, p. 268. [8] S. Di Cairano, W. Liang, I. Kolmanovsky, M. Kuang, and A. Phillips, “Engine power smoothing energy management strategy for a series hybrid electric vehicle,” in Proc. of the ACC, 2011, pp. 2101–2106. [9] P. Ortner and L. del Re, “Predictive control of a diesel engine air path,” IEEE Trans. on CST, vol. 15, no. 3, pp. 449–456, 2007. [10] T. Broomhead, C. Manzie, L. Eriksson, M. Brear, and P. Hield, “A robust model predictive control framework for diesel generators,” in IFAC World Congress, 2014. [11] J. Wahlstr¨om and L. Eriksson, “Modelling diesel engines with a variable-geometry turbocharger and exhaust gas recirculation by optimization of model parameters for capturing non-linear system dynamics,” J. of Auto. Eng., vol. 225, no. 7, pp. 960–986, 2011. [12] I. Kolmanovsky and E. G. Gilbert, “Maximal output admissible sets for discrete-time systems with disturbance inputs,” in Proc. of the ACC, vol. 3, 1995. [13] P. J. Goulart, E. C. Kerrigan, and J. M. Maciejowski, “Optimization over state feedback policies for robust control with constraints,” Automatica, vol. 42, no. 4, pp. 523–533, 2006. [14] D. Q. Mayne and J. B. Rawlings, Model Predictive Control: Theory and Design. Wisconsin, USA: Nob Hill Publishing, LLC, 2009. [15] Z.-P. Jiang and Y. Wang, “Input-to-state stability for discrete-time nonlinear systems,” Automatica, vol. 37, no. 6, pp. 857–869, 2001. [16] R. Sharma, D. Nesic, and C. Manzie, “Model reduction of turbocharged (TC) spark ignition (SI) engines,” IEEE Trans. on CST, vol. 19, no. 2, pp. 297–310, 2011. [17] D. Lim´on, I. Alvarado, T. Alamo, and E. F. Camacho, “MPC for tracking piecewise constant references for constrained linear systems,” Automatica, vol. 44, no. 9, pp. 2382–2387, 2008.
2516
Robust Stable Economic MPC with Applications in Engine Control Timothy Broomhead1 , Chris Manzie1 , Rohan Shekhar1 , Michael Brear1 and Peter Hield2
Abstract— Economic Model Predictive Controllers have shown to improve a plant’s economic performance using state dependant economic stage costs. Recent extensions have provided continual feasibility guarantees, despite changes in economic parameters, however, perfect plant models have been assumed. This assumption is invalid in practise due to modelling errors or un-modelled disturbances and can therefore lead to infeasibility of the optimisation problem. This paper proposes a robust economic model predictive controller, which takes advantage of constraint tightening techniques to guarantee feasibility despite modelling errors. Inputto-state stability is proven using a Lyapunov function. The advantages of this method are highlighted against alternative control structures in the application of power tracking for diesel engines in series hybrid type applications.
I. I NTRODUCTION To minimise a plant’s economic cost, Model Predictive Control (MPC) is often employed as part of a hierarchical control architecture, where a high level, static Real Time Optimiser (RTO), chooses operating set-points which minimise an economic cost [1]. The lower level Tracking Model Predictive Control (TMPC) then steers the plant to the chosen set-points in an admissible way, taking into account plant dynamics and constraints. Typically, the TMPC problem is formulated using an error system, with convergence of the closed-loop system guaranteeing set point tracking. Demand for higher efficiencies and reduced economic costs is driving development of alternative approaches. The TMPC is designed to reject disturbances and track the system as closely as possible to the given set-points, however it does not necessarily take into account the economic cost of the state trajectories. Therefore, if the set-points change frequently with respect to the plant dynamics, this approach may result in poor economic performance. An alternative is to combine the RTO and TMPC into a single controller, known as Economic Model Predictive Control (EMPC) [2]. In [2] a Lyapunov function is used to show stability of the EMPC formulation under the assumption of strong duality with the steady state problem and using a terminal constraint. The assumption of strong duality has been relaxed to a dissipativity assumption in [3]. The need for terminal constraints can be relaxed for sufficiently long prediction horizons where controllability and turnpike properties hold [4], however in some applications the required horizon length for acceptable performance may be computationally prohibitive.
Economic parameters may change as a result of variation in material prices or requested plant production rates. With a change in economic parameters, the optimal steady state condition may shift and potentially result in a loss of feasibility for controllers with terminal constraints as dependant on this state. Recent work in [5] extends the results of [2] relaxing the terminal constraint to be any point on the steady state manifold and including a terminal weight. Referred to as Economic Tracking Model Predictive Control (ETMPC) in this paper, this approach requires modification of the stage cost to ensure convergence to the optimal steady state, resulting in sub-optimal transient behaviour, though feasibility is guaranteed for changes in economic parameters. Each of the approaches introduced so far assumes a perfect prediction model to ensure continual feasibility. In practice, this assumption is invalid due to external disturbances, measurement error and modelling errors, potentially resulting in infeasibility through constraint violation. Many existing robustness techniques, such as open or closed loop Min-Max MPC, are overly conservative or computationally intractable [6]. Constraint Tightening (CT) [7] is a robustness technique for which robust feasibility can be guaranteed. In the general case, computational requirements for CT do not exceed that of nominal robustness. The basic principle is that constraints are systematically tightened along the horizon to reserve a margin which can then be used to correct for errors resulting from model mismatch at some future time. This paper combines the results of [5] with the constraint tightening principles in [7] to produce a practical control structure capable of directly addressing economic objectives and constraints whilst being robust to modelling errors. In the application of diesel generator control, the economic objective is to minimise the penalty associated with not meeting the requested output power, whilst minimising the cost associated with fuel consumption under emissions and safety constraints. While previous approaches utilising MPC for engine control often represent TMPC formulations [8], [9], [10], an economic formulation can increase overall performance. Representation of diesel engines requires complex models [11], which can be impractical for control design. Model reduction techniques are therefore required, resulting in modelling error, motivating the development of a Constraint Tightened Tracking Economic Model Predictive Control (CT-TEMPC).
The authors would like to thank the Australian Defence Science and Technology Organisation, the Defence Science Institute and the Elizabeth and Vernon Puzey trust for their support of this research. 1 Department of Mechanical Engineering, The University of Melbourne, Australia. [email protected] 2 Defence Science and Technology Organisation, Melbourne, Australia
978-1-4673-6090-6/14/$31.00 ©2014 IEEE
2511
II. P ROBLEM P RELIMINARIES Consider a linear-time invariant system xk+1 = Axk + Buk + wk
(1)
where x ∈ Rn is the state, u ∈ Rm is the control input and w ∈ Rn is an unknown but bounded disturbance, obeying wk ∈ W. The model outputs y ∈ Rl are given by yk = Cxk + Duk
(2)
the optimal trajectory. The following assumption is therefore made in order to demonstrate asymptotic stability. Assumption 5 (Strong Duality of Steady-State problem): Let ˜l(x, u) be the rotated stage cost function given by ˜l(x, u, p) = l(x, u, p) + λT (x − Ax − Bu) − l(x∗ , u∗ , p) s s
and constraints must be satisfied for all k, described by yk ∈ Y
(3)
where Y is a compact and convex set. Disturbances on the output (2) could be incorporated as in [10], however are omitted due to space limitations. Assumption 1: W is compact, convex and contains 0. Assumption 2: There exists a sequence of nilpotent controller gains Kj which drive the system xj+1 = (A + BKj )xj , j ∈ I0:Nnp to the origin in Nnp < N − 1 steps, where Ia:b denotes the set of integers from a to b. The tightened constraint sets for robustness are defined as Y0 = Y, Yj+1 = Yj ∼ [C + DKj ]Lj W, ∀j ∈ I1:N
(4)
where λ is a multiplier ensuring a unique minimum at (x∗s , u∗s ) ∀x, u given Cx + Du ∈ Y s . Then there exists a κ∞ -function α(.) such that ˜l(x, u, p) ≥ α(|x − xs |). III. C ONTROL F ORMULATION Definition 1: For the a given constant γ > |λ| a convex offset function VO (x, u) satisfies VO (x, u) − VO (xs , us ) ≥ γ|x − xs | and for which min VO (x, u)
(x,u)∈Zs
has a unique solution at (x∗s , u∗s ). The economic cost function for changing parameters is [5] lt (x, u, p) = l(x + x∗s , u + u∗s , p)
where ”∼” denotes the Pontryagin difference [12] and Lj is L0 = I, Lj+1 = (A + BKj )Lj , ∀j ∈ I1:N .
VN (x, p; u) =
(5)
(6b)
s = YN ∼ B and B is an arbitrarily small ball of where YN radius > 0, required for stability guarantees. We wish to generate a robustly feasible feedback control law which aims to minimise cost P∞ (7) V∞ = k=0 l(xk , uk , p)
where economic parameters p may vary in time. The best robustly achievable steady state are (x∗s (p), u∗s (p)) = arg min l(x, u, p) s.t (x, u) ∈ Zs . x,u
(10a)
lt (xk+j − xk+N −1 , uk+j − uk+N −1 , p)
+ VO (xk+N −1 , uk+N −1 )
(8)
Henceforth, dependence of x∗s , u∗s on p is removed for brevity. Assumption 4: The stage cost function l(x, u) is Lipschitz continuous in (X × U) with Lipschitz constant Ml , i.e.: |l(x1 , u1 , p) − l(x2 , u2 , p)| ≤ Ml (|x1 − x2 | + |u1 − u2 |) ∀x1 , x2 ∈ X , ∀u1 , u2 ∈ U. In EMPC, standard Lyapunov arguments cannot be used since the economic cost is not necessarily decreasing along
(10b)
The optimisation problem PN (x, p) is given as PN (x, p) : arg min VN (x, p; u) u
s.t xk = x
(6a)
ys = Cxs + Dus } Xs = {xs : ∃u(x, u) ∈ Zs }
N −1 X j=0
Remark 1: Since the controller gains are nilpotent, Lj = 0 ∀j ≥ Nnp . If a disturbance feedback controller is utilised for disturbance rejection, the controller gains may be optimised as part of a convex optimisation [13]. Assumption 3: The constraint set (3) is sufficiently large with respect to the disturbance set W and Kj are chosen such that the resulting tightened sets are non-empty. Define the set of robustly admissible equilibrium conditions and states for the nominal system (where wk = 0) as Zs = {(xs , us ) : ys ∈ Y s , xs = Axs + Bus ,
(9)
(11a) (11b)
xk+N ∈ Xs
(11c)
xk+j+1 = Axk+j + Buk+j
(11d)
yk+j = Cxk+j + Duk+j
(11e)
yk+j ∈ Yj , ∀j ∈ I0:N −1
(11f)
while the robust economic MPC control law κeN (x, p) is implicitly defined as the optimal solution to PN (x, p). IV. S TABILITY A NALYSIS A. Recursive Feasibility Theorem 1 (Feasibility): Consider a system described by (1) and (2), subject to constraints (3) and let Assumptions 1, 2 and 3 hold. Under the control law given by PN (x, p) and given that the system is initially within the feasibility set XN , the system will remain within XN for all time. Proof: Consider any xk0 ∈ XN . Define the optimal solution to PN (xk0 , p) as u0k0 +j , with corresponding states x0k0 +j and outputs yk00 +j , ∀j ∈ I0:N −1 . Consider the candidate solution for k0 + 1 u ˆk0 +j = u0k0 +j + Kj−1 Lj−1 wk0 ∀j ∈ I1:N −1 u ˆk0 +N =
u0k0 +N −1
(12a) (12b)
This solution is created by shifting the previous solution by one step, adding perturbations for disturbance rejection and appending one step of control to the end. Due to the steady state terminal condition, the appended control step can be the final control action of the previous solution. The initial state
2512
at k0 + 1 is given by the system dynamics (1). Substitution into the internal model (11d) gives xk0 +1 = x0k0 +1 + wk0 . Substituting this initial condition and the candidate solution (12) gives the candidate state sequences ∀j ∈ I1:N x ˆk0 +j =
x0k0 +j
+ Lj−1 wk0
x ˆk0 +N +1 = Aˆ xk0 +N + B u ˆk0 +N
Proof: Using the Cauchy-Schwarz’s inequality, λT (x − ≤ |λ||x − x∗s |. Let γˆ = γ − |λ|, then V˜O (x, u) = VO (x, u) + λT (x − x∗s ) − VO (x∗s , u∗s )
x∗s )
≥ γ|x − x∗s | + λT (x − x∗s )
(13a) (13b)
and output sequences ∀j ∈ I1:N −1 yˆk0 +j = yk00 +j + (C + DKj−1 )Lj−1 wk0
(14a)
yˆk0 +N = C x ˆk0 +N + Dˆ uk0 +N
(14b)
Initial condition, state equality and output equality constraints are therefore satisfied by construction. The satisfaction of the terminal constraint (11c) uses that for a nilpotent controller Lj = 0 for j ≥ Nnp . Therefore, since Nnp < N − 1, (13a) implies x ˆk0 +N = x0k0 +N . Feasibility 0 at k0 implies xk0 +N is a steady state with input u0k0 +N −1 , substitution into (13b) gives x ˆk0 +N +1 = x0k0 +N , therefore x ˆk0 +N +1 = x ˆk0 +N and the terminal constraint is satisfied. Feasibility at k0 requires yk00 +j ∈ Yj ∀j ∈ I1:N −1 . Since Yj−1 is related to Yj by (4), (14a) implies yˆk0 +j ∈ Yj−1 ∀j ∈ I1:N −1 . This satisfies the output constraints for all but the final step N . Since x ˆk0 +N = x0k0 +N −1 0 and u ˆk0 +N |k0 +1 = uk0 +N −1|k0 then yˆk0 +N = yk00 +N −1 . Since at time k0 feasibility requires yk00 +N −1 ∈ YN −1 then yk00 +N ∈ YN −1 ⊆ YN . This demonstrates the candidate solution satisfies all constraints, and hence for any x ∈ XN a feasible solution exists. B. Input-to-state stability
= γˆ |x − x∗s | + |λ||x∗s − x| − λT (x∗s − x) ≥ γˆ |x − x∗s |. Lemma 3: Consider a system described by (1) and (2), subject to constraints (3). Let Assumptions 1 - 5 hold. Under the control law given by (11) and given x ∈ XN , there exists a κ∞ -function, α(.), and κ-function, ρ(.), such that V˜ 0 (x+ , p) − V˜ 0 (x, p) ≤ −α(|x − x0 |) + ρ(|w|) (20) N
Proof: Consider at k0 any xk0 ∈ XN and recall the candidate solution (12) for the problem at k0 + 1. Compare the optimal cost V˜N0 (xk0 ) with that of the candidate solution V˜N (ˆ xk0 +1 ) − V˜N0 (xk0 ) = PN −1 ˜ 0 − j=0 lt (xk0 +j − x0k0 +N −1 , u0k0 +j − u0k0 +N −1 , p) PN ˜ + lt (ˆ xk +j − x ˆ0 ,u ˆk +j − u ˆ0 , p) j=1
which has the following properties l˜t (0, 0, p) = ˜l(x∗s , u∗s , p) = 0 l˜t (z, v, p) ≥ α1 (|z|) for the κ∞ -function α1 (.).
(16a) (16b)
A rotated offset function is used, defined as V˜O (x, u) = VO (x, u) + λT (x − x∗s ) − VO (x∗s , u∗s ),
(17)
with λ from Assumption 5. The auxiliary cost function is V˜N (x, p; u) =
N −1 X
l˜t (xk+j − xk+N −1 , uk+j − uk+N −1 , p)
j=0
+ V˜O (xk+N −1 , uk+N −1 ),
(18)
and the auxiliary optimisation problem is P˜N (x, p) : arg min V˜N (x, p; u) s.t (11b) − (11f ). (19a) u
Lemma 1 (Lemma 1 of [5]): The problem P˜N (x, p) gives the same optimal control sequence as PN (x, p). Lemma 2: Due to Definition 1 there exists a positive γˆ such that V˜O (x, u) ≥ γˆ |x − x∗s |
0
k0 +N
Since Nnp < N − 1 and by the definition of u ˆ, we know that x0k0 +N −1 = x ˆk0 +N and u0k0 +N −1 = u ˆk0 +N so the offset functions equate. Furthermore, from Assumptions 4 and 5 and (15) we can establish a bound on the Lipschitz ˜ l ≤ Ml + |λT (I − A)| + |λT B|). constant for l˜t (x, u, p) as M Additionally, taking into account that the disturbance is nilpotent for j ≥ Nnp and property (16a) we find VˆN (xk +1 ) − V 0 (xk ) 0
(15)
k0 +N
0
−VO (x0k0 +N −1 , u0k0 +N −1 ) + VO (ˆ xk0 +N , u ˆk0 +N ).
N
0
≤ − l˜t (x0k0 − x0k0 +N −1 , u0k0 − u0k0 +N −1 , p) ˜ l PNnp (|Lj−1 | + |Kj−1 Lj−1 |)|wk |. +M
For the purposes of showing stability, an auxiliary problem is introduced, with stage cost l˜t (x, u, p) = ˜l(x + x∗s , u + u∗s , p)
N −1
N
j=1
0
xk0 +1 ) and due to property By optimality V˜N0 (xk0 +1 ) ≤ V˜N (ˆ (16b) we can state V˜N (x+ ) − V˜N0 (x) ≤ −α1 (x − xN −1 ) + ˜ l PNnp (|Lj−1 | + |Kj−1 Lj−1 |)ν. ρ(|w|), where ρ(ν) = M j=1 Define the constant PNnp γ0 = Ml j=0 (|(1 + Kj )Lj |)|. (21) Lemma 4: Consider a system described by (1) and (2), subject to constraints (3). Let Assumptions 1 - 5 hold. Let (ˆ xs , u ˆs ) ∈ Zs and u ˆ be a sequence of control actions u ˆj = Kj (xj − x ˆs ) + u ˆs , where xj+1 = Axj + B u ˆj . Then under Assumption 2, (ˆ xs , u ˆs ) can be chosen such that u ˆ is a feasible solution to P˜N (x) and such that |x − x ˆs | ≤ ω(), ω > 0 given |xj − x∗s | ≤ ∀j ∈ I0:N −1 . It then follows that ˆ ) ≤ γ0 |x − x V˜N (x, p; u ˆs | + VO (ˆ xs , u ˆs ) + d, where d = λT (x − xs ) − VO (x∗s , u∗s ), and γ0 is as defined in (21). Proof: Noting from (5) that uj = Kj Lj (x − x∗s ) + u∗s , xj = Lj (x − x∗s ) + x∗s , and Lj = 0, ∀j > Nnp , then PN −1 ∗ ∗ j=0 (l(xj , uj , p) − l(xs , us , p)) PNnp ≤ Ml j=0 (|Lj (x − x∗s )| + |Kj Lj (x − x∗s )|) PNnp ≤ Ml j=0 (|(1 + Kj )Lj |)|x − x∗s | = γ0 |x − x∗s |.
2513
It can then be shown that P −1 ˆ) = N V˜N (x, p; u ˆs + x∗s , uj − u ˆs + u∗s , p) j=0 (l(xj − x − l(x∗s , u∗s , p)) + VO (ˆ xs , u ˆs ) + d ≤ γ0 |x − x ˆs | + VO (ˆ xs , u ˆs ) + d, where d = λ (x − x∗s ) − VO (x∗s , u∗s ). Lemma 5: Consider a system described by (1) and (2), subject to constraints (3). Let Assumptions 1 - 5 hold. Consider the constant γ > γ0 as defined in (21). Assume for an initial state x the optimal solution to P˜N (x, p) is (x0 ,u0 ). If x0N −1 = x and u0N −1 = κeN (x, p), then T
x0N −1 = x∗s , u0N −1 = u∗s
(22)
Proof: Consider that the optimal solution to P˜N (x, p) is (x0 , u0 ), x = x0N −1 and κeN (x, p) = u0N −1 . Then (x, u) ∈ Zs and the optimal cost function is V˜N0 (x, p) = V˜O (x0N −1 , u0N −1 ). In order to prove the lemma by contradiction, first consider that (x0N −1 , u0N −1 ) 6= (x∗s , u∗s ). Define (ˆ xs , u ˆs ) = β(x0N −1 , u0N −1 ) + (1 − β)(x∗s , u∗s ), and note that due to convexity of Zs , (ˆ xs , u ˆs ) ∈ Zs . Following from Assumption 2, there exists a βˆ ∈ (0, 1) such that for any ˆ 1), the control law uj = Kj (xj − x β ∈ [β, ˆs ) + u ˆs drives ˆs in an admissible way in ≤ N the system from x0N −1 to x steps. Define u ˆ as the control law uj = Kj (xj − x ˆs ) + u ˆs , u ˆ is a feasible solution to P˜N (x0N −1 , p). From Lemma 4 we have that V˜N0 (x0N −1 , p) = V˜O (x0N −1 , u0N −1 ) ≤ V˜N (x0N −1 , p; u ˆ) ≤ γ0 (1 − β)|x0N −1 − xs | + VO (ˆ xs , u ˆs ) + d where d = λT (x0N −1 − xs ) − VO (xs , us ). Define W (x0N −1 , β) = γ0 (1 − β)|x0N −1 − x∗s | + VO (ˆ xs , u ˆs ) + d. The partial derivative of W about β is ∂W /∂β =
−γ0 |x0N −1
− xs | + g
0
(x0N −1
−
x∗s , u0N −1
−
u∗s )
Theorem 2 (Input-to-state stability): Consider a system described by (1) and (2), subject to constraints (3). Let Assumptions 1 - 5 hold. Under the control law given by (11) and given x ∈ XN , the closed loop system will be input-to-state stable, with a domain of attraction XN for all γ > γ0 . Proof: Due to property (16b) and Lemma 2 we have that V˜N0 (x, p) ≥ α1 (|x − xN −1 |) + γˆ |xN −1 − x∗s |, for all x ∈ XN . Hence there is a κ∞ -function, αJ (.), such that V˜N0 (x, p) ≥ αJ (|x − xN −1 | + |xN −1 − xs |) ≥ αJ (|x − x∗s |). Since V˜N0 (x∗s ) = 0 and V˜N0 (x) is continuous in XN , there exists a κ∞ -function βJ such that for all x ∈ XN [14, Proposition 2.18] V˜N0 (x, p) ≤ βJ (|x − x∗s |). Let e(x) = x − x0N −1 for all x ∈ XN and x0N −1 ∈ Xs . Then from Lemma 3, V˜N0 (x+ , p) − V˜N0 (x, p) >= −α1 (|e(x)|) + ρ1 (|w|), where ρ1 is positive definite. If e(x) = 0, then x = x0N −1 , and from Lemma 5 this requires that x0N −1 = x∗s ; while if x = x∗s , then by optimality x0N −1 = x∗s and hence e(x) = 0. Hence we can conclude that |e(x)| = 0 iff x = x∗s . Since XN is compact, |e(x)| = 0 when x = x∗s and otherwise |e(x)| > 0, there exists an κ-function, αe (.), such that |e(x)| ≥ αe (|x − x∗s |). It then follows that α1 (|e(x)|) ≥ α1 (αe (|x − x∗s |)) = αV (|x − x∗s |) where αV (.) is a κ∞ -function. It can now be concluded that: 1) αJ (|x − xs |) ≤ V˜N0 (x, p) ≤ βJ (|x − xs |), ∀x ∈ XN 2) V˜N0 (x+ , p) − V˜N0 (x, p) ≤ −αV (|x − xs |) + ρ1 (|w|), ∀x ∈ XN , ∀w ∈ W Therefore V˜N0 (x, p) is an ISS-Lyapunov function in XN , and xs is a local input-to-state (ISS) stable equilibrium point for the closed loop system controlled by P˜N (x) [15]. As a result of Lemma 1, ISS stability also extends to PN (x). V. A PPLIED E XAMPLE A. Problem Formulation
where g 0 is a sub-differential of VO (ˆ xs , u ˆs ). Evaluating this derivative for β = 1 ∂W /∂β|β=1 = −γ0 |x0N −1 −x∗s |+¯ g 0 (x0N −1 −x∗s , u0N −1 −u∗s ) where g¯0 is a sub-differential of VO (x0N −1 , u0N −1 ). Since VO is convex, and due to Definition 1 g¯0 (x0N −1 − x∗s , u0N −1 − u∗s ) ≥ VO (x0N −1 , u0N −1 ) − VO (x∗s , u∗s ) ≥ γ|x0N −1 − x∗s | Therefore ∂W /∂β|β=1 ≥ −γ0 |x0N −1 − x∗s | + γ|x0N −1 − x∗s | = (γ − γ0 )|x0N −1 − x∗s | Since γ > γ0 and x0N −1 6= x∗s , ∂W /∂β|β=1 > 0. This ˆ 1) such that the cost would mean that there exists a β ∈ [β, to move the system from x0N −1 to x ˆs , is smaller than to remain in x0N −1 . This contradicts optimality of the solution to P˜N (x, p) and hence x0N −1 = x∗s .
In this example, the developed robust economic MPC formulation will be utilised in a diesel generator power tracking application. The economic objective is associated with minimising the error between the requested power and the brake power produced by the generator, with a secondary economic objective to minimise fuel consumption. To simulate the plant, a calibrated and validated high order Mean Value Engine Model (MVEM) is utilised, similar to that described in [10]. The model inputs u := [tδ τload uvgt ]T are the injection duration, load torque and Variable Geometry Turbocharger (VGT) demand. The system states x := [ωe pim pem ωt xvgt ]T are engine speed, intake and exhaust manifold pressure, turbocharger speed and VGT position. Since the MVEM is highly nonlinear, and computationally impractical for use in an online control strategy, model reduction techniques described in [16] are utilised. The relatively fast manifold pressure dynamics are eliminated before the reduced model is linearised around a fixed operating point. Since the controller will be tested against the full order MVEM, the model reduction results in modelling error,
2514
l(x, u, p) = (ωe τeng − Pref )2 + (Ψ · m ˙ f (ωe , tδ ))2
10
(23)
where the requested power is the economic parameter Pref and m ˙ f (ωe , tδ ) is the fuel consumption map with weighting Ψ, assumed small in this example. In the sequel, a number of different control strategies will be examined, and their ability to minimise (7) compared. Where relevant, the steady state conditions (xs , us ) used in controller formulations are found for each Pref using (8). The reference Pref is constructed with step changes every 5 seconds of a random size, to represent an aggressive test-case. The power request is limited between 20 kW and 57 kW with a maximum absolute step size of 20 kW. A timestep of 100 ms has been chosen, based on the bandwidth of the engine’s sensors and actuators. Denote the controller as CT-TEMPC, constructed using optimisation problem (11) and the offset function VO = kQV (xk+N − xs )k1 + kRV (uk+N − us )k1 ,
(24)
where QV , RV are chosen for good transient performance. Tightened constraint sets are calculated using (4). A Monte Carlo approach is used to identify the disturbance set W, capturing the error between the reduced and full order models. The controller CT-TEMPC, initially with a conservative estimate for W, is requested to track the randomly generated reference signal, Pref , for a large number of samples. The modelling errors between the reduced and full order models are captured at each time-step. To reduce conservatism, lower and upper error bounds were chosen at the 0.1 and 99.9 percentiles to define the set W. This process was iterated upon several times to generate the final error set, shown
−2
−1 rpm
0
Turbine Speed
1
−4
−2
0 2 3 rpm×10
−0.1
0 %
2 0
25
10
15
20
25
N
B. Comparison of economic control formulations One approach to economic MPC is to use the pure economic cost as the stage cost without terminal constraints, as described in [4]. This eliminates the need to determine optimal steady state points explicitly, potentially leading to a good approximation of (7) assuming an accurate internal model, and sufficiently long horizon. Denote the robust version of this controller as CT-EMPC, constructed using the optimisation PN −1 arg minu j=0 l(xj , uj , p) s.t (11b), (11d) − (11f ). (25) In Fig. 2 this controller is demonstrated with a power request at 20 kW, changing to 57 kW at 10 seconds. In this application, a horizon length less than 16 results in unacceptable steady state tracking error. Furthermore, with a horizon length of 16, the nonlinear optimisation is computationally intractable for online control of diesel engines. A horizon length of around 10 is more realistic in terms of computational load, motivating alternative controller formulations. One option is to use a robustified version of the EMPC in [2], constructed using optimisation (25) with an additional terminal constraint xk+N = x∗s , denoted as CT-EMPC-D. However, it was found that the controller became infeasible after changes in power request with a reasonable horizon length of 10, due to a small region of attraction. As shown in [5], a sufficiently large weighting on the offset function results in the output of (11) being equal to CT-EMPC-D Ref
CT−TEMPC−H
CT−TEMPC
TEMPC
50 40 30 0
10
20
30
40
50
60
10
20
30
40
50
60
10
20
30 Time (s)
40
50
60
10 0 −10 0
0.1
Fig. 1. Prediction errors for each state. Shown are 0.1 and 99.9 percentiles
4
in Fig. 1. Due to the Monte Carlo approach and the error bounds used, absolute guarantees can no longer be made for the controller; however good practical robustness remains, as demonstrated later.
VGT Position
4
15 20 Time (s)
6
Fig. 2. Left: Closed-loop output of CT-EMPC after a step-change in power request. Right: Steady state tracking error.
Engine Speed 3 (rpm) × 10
Engine Speed
N=10 N=13 N=16 Ref
40
20
which will be accounted for through constraint tightening. To represent measurement noise, a small normally distributed disturbance is added to the simulated states, which can be treated as an additive disturbance [7]. The engine speed, turbine speed and each input has lower and upper operating bounds, while VGT actuator also has a maximum slew rate constraint, summarised in Table I. The normalised air-to-fuel ratio of the engine, λ, has a lower limit in order to prevent excessive smoke production, set to 1.4 in this application. This non-linear, non-convex constraint is approximated using a convex polytopic constraint. The power tracking economic objective can be summarised with the following stage cost
Final Tracking Error (kW)
Power Output (kW)
∆uvgt (%) −5 5
Power Output (kW)
Lower Upper
60
Power Error (kW)
Bounds
TABLE I L OWER AND U PPER C ONSTRAINTS ωe ωt uvgt tδ τload (rpm) (rpm) (%) (ms) (Nm) 1500 50 × 103 65 0.35 0 2500 140 × 103 85 1 400
2
1.5 0
Fig. 3. Comparison of Economic MPC formulations. The TEMPC controller fails after 3 seconds, indicated on the Engine Speed plot.
2515
when inside its region of attraction. Denote this controller as CT-TEMPC-H, where H indicates that QV is chosen to be large. In Fig. 3, CT-TEMPC-H is compared against CT-TEMPC and its non-robust equivalent, TEMPC, both with the same offset function (24). All controllers use a horizon length of N = 10. The non-robust controller, TEMPC, fails due to an engine speed constraint violation, indicated on Fig. 3 at around 3 seconds. Due to the finite horizon lengths, offset functions, economic stage cost for tracking and model mismatch, each of the controllers presented are suboptimal; this is the price paid for guaranteed feasibility and convergence. Selection of the offset function directly influences the closed loop performance, as is made clear in the comparison between CT-TEMPC-H and CT-TEMPC. During large step changes, CT-TEMPC-H increases tracking error in order to ensure fast convergence of the terminal state to x∗s . This turns out to be undesirable behaviour, which can be removed with a different offset function, as seen in the performance of CT-TEMPC.
It can be seen that the CT-TEMPC controller out-performs the tracking controller during the transients, reaching the target power output in a shorter period of time, this is particularly noticeable when power demand is reduced. Despite the apparent similarity of the trajectories, the cumulative economic cost (23) of the economic controller over the simulation was 28 % less than that of the tracking controller. Additionally, the time-consuming tuning process for the tracking controllers must be performed for any changes in economic parameters, further motivating the use of the economic controller. VI. C ONCLUSIONS Through the use of constraint tightening, a robust economic MPC has been presented. The controller is guaranteed to remain feasible despite modelling errors and changes in economic parameters. An ISS result is given based on a Lyapunov function, ensuring convergence to a set around the optimal steady state. The efficacy of the developed controller has been demonstrated against alternative formulations in the application of diesel generator power tracking.
C. Tracking controller comparison
R EFERENCES
The developed economic controller CT-TEMPC will now be compared against a robustified TMPC formulation [17], representing the hierarchical set-point tracking approach. The robust tracking controller, denoted CT-TMPC, is described by the optimisation PN −1 arg minu j=0 |xk+j − xk+N |2Q + |uk+j − uk+N |2R + VO s.t (11b) − (11f ),
(26)
where VO is from (24), with the same tuning as CT-TEMPC. The weights Q and R have been tuned using Matlab’s Optimisation Toolbox to minimise the overall economic cost over a 150 second simulation using a randomly generated reference power trajectory. The optimisation used the internal model as the plant model, representing an off-line tuning process, and a Cholesky decomposition of Q and R to enforce positive definiteness. In Fig. 4, CT-TMPC is compared against CT-TEMPC. Both controllers exhibit negligible steady state tracking error, however their closed-loop behaviour differs during transients.
Power Output (kW)
Ref
CT−TMPC
CT−TEMPC
50 40 30 0
10
20
30
40
50
60
10
20
30
40
50
60
10
20
30 Time (s)
40
50
60
Power Error (kW)
10 0
Engine Speed 3 (rpm) × 10
−10 0 2
1.5 0
Fig. 4.
Comparison of tracking and economic MPC formulations
[1] S. J. Qin and T. A. Badgwell, “A survey of industrial model predictive control technology,” Contr. Eng. Prac., vol. 11, pp. 733–764, 2003. [2] M. Diehl, R. Amrit, and J. B. Rawlings, “A lyapunov function for economic optimizing model predictive control,” IEEE Trans. on Auto. Contr., vol. 56, no. 3, pp. 703–707, 2011. [3] D. Angeli, R. Amrit, and J. B. Rawlings, “On average performance and stability of economic model predictive control,” IEEE Trans. on Auto. Contr., vol. 57, no. 7, pp. 1615–1626, 2012. [4] L. Gr¨une, “Economic receding horizon control without terminal constraints,” Automatica, vol. 49, no. 3, pp. 725–734, 2013. [5] A. Ferramosca, D. Limon, and E. Camacho, “Economic MPC for a changing economic criterion for linear systems,” Automatic Control, IEEE Transactions on, vol. In press, 2014. [6] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. Scokaert, “Constrained model predictive control: Stability and optimality,” Automatica, vol. 36, no. 6, pp. 789–814, 2000. [7] A. Richards and J. How, “Robust model predictive control with imperfect information,” in Proceedings of the ACC, 2005, p. 268. [8] S. Di Cairano, W. Liang, I. Kolmanovsky, M. Kuang, and A. Phillips, “Engine power smoothing energy management strategy for a series hybrid electric vehicle,” in Proc. of the ACC, 2011, pp. 2101–2106. [9] P. Ortner and L. del Re, “Predictive control of a diesel engine air path,” IEEE Trans. on CST, vol. 15, no. 3, pp. 449–456, 2007. [10] T. Broomhead, C. Manzie, L. Eriksson, M. Brear, and P. Hield, “A robust model predictive control framework for diesel generators,” in IFAC World Congress, 2014. [11] J. Wahlstr¨om and L. Eriksson, “Modelling diesel engines with a variable-geometry turbocharger and exhaust gas recirculation by optimization of model parameters for capturing non-linear system dynamics,” J. of Auto. Eng., vol. 225, no. 7, pp. 960–986, 2011. [12] I. Kolmanovsky and E. G. Gilbert, “Maximal output admissible sets for discrete-time systems with disturbance inputs,” in Proc. of the ACC, vol. 3, 1995. [13] P. J. Goulart, E. C. Kerrigan, and J. M. Maciejowski, “Optimization over state feedback policies for robust control with constraints,” Automatica, vol. 42, no. 4, pp. 523–533, 2006. [14] D. Q. Mayne and J. B. Rawlings, Model Predictive Control: Theory and Design. Wisconsin, USA: Nob Hill Publishing, LLC, 2009. [15] Z.-P. Jiang and Y. Wang, “Input-to-state stability for discrete-time nonlinear systems,” Automatica, vol. 37, no. 6, pp. 857–869, 2001. [16] R. Sharma, D. Nesic, and C. Manzie, “Model reduction of turbocharged (TC) spark ignition (SI) engines,” IEEE Trans. on CST, vol. 19, no. 2, pp. 297–310, 2011. [17] D. Lim´on, I. Alvarado, T. Alamo, and E. F. Camacho, “MPC for tracking piecewise constant references for constrained linear systems,” Automatica, vol. 44, no. 9, pp. 2382–2387, 2008.
2516
