- Email: [email protected]
Nonlinear MPC with fixed-point iterations: numerical and experimental validation
16th IFAC workshop on Control Applications of Optimization 16th workshop on of 16th IFAC IFAC6-9, workshop on Control Control Applications ApplicationsGermany of Optimization Optimization October 2015. Garmisch-Partenkirchen, 16th IFAC workshop on Control Applications of Optimization October 6-9, 2015. Garmisch-Partenkirchen, Germany Available online at www.sciencedirect.com October 6-9, 2015. Garmisch-Partenkirchen, Germany October 6-9, 2015. Garmisch-Partenkirchen, Germany
ScienceDirect IFAC-PapersOnLine 48-25 (2015) 062–067
Nonlinear MPC with fixed-point iterations: Nonlinear MPC with fixed-point iterations: Nonlinear MPC with fixed-point iterations: numerical and experimental validation experimental validation numerical and numerical and experimental validation Knut Graichen Knut Graichen Knut Knut Graichen Graichen Institute of Measurement, Control and Microtechnology, Institute of Measurement, Control and Microtechnology, Institute of Control and University of Ulm, Germany (e-mail: [email protected]). Institute of Measurement, Measurement, Control and Microtechnology, Microtechnology, University of Ulm, Germany (e-mail: [email protected]). University of Ulm, Germany (e-mail: [email protected]). University of Ulm, Germany (e-mail: [email protected]). Abstract: The paper demonstrates the applicability of a recently presented MPC scheme for Abstract: The paper the applicability of a recently presented MPC scheme for Abstract: The demonstrates the of presented MPC for nonlinear withdemonstrates control constraints. The algorithm is based on fixed-point iterations Abstract:systems The paper paper demonstrates the applicability applicability of aa recently recently presented MPC scheme scheme for nonlinear systems with control constraints. The algorithm is based on fixed-point iterations nonlinear systemsconsist with control control constraints. The algorithm algorithm is based based on on fixed-point iterations that in principal of combined forward/backward integrations over the MPCiterations horizon. nonlinear systems with constraints. The is fixed-point that in principal consist of combined forward/backward integrations the MPC horizon. that consist of forward/backward integrations over the MPC horizon. Numerical and runtime for a nonlinear chemical reactor modelover show efficiency of that in in principal principal consistresults of combined combined forward/backward integrations over thethe MPC horizon. Numerical and runtime results for aa nonlinear chemical reactor model show the efficiency of Numerical and runtime results for nonlinear chemical reactor model show the efficiency of the fixed-point iterative MPC scheme. Experimental results for a magnetic levitation system Numerical and runtime results for a nonlinear chemical reactor model show the efficiency of the fixed-point iterative MPC scheme. Experimental results for aa magnetic levitation system the fixed-point iterative MPC scheme. Experimental results for magnetic levitation system confirm its practical applicability. the fixed-point iterative MPC scheme. Experimental results for a magnetic levitation system confirm confirm its its practical practical applicability. applicability. confirm its practical applicability. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Model predictive control, fixed-point iterations, control constraints, real-time Keywords: Model Model predictive control, control, fixed-point iterations, iterations, control constraints, constraints, real-time Keywords: implementation. Keywords: Model predictive predictive control, fixed-point fixed-point iterations, control control constraints, real-time real-time implementation. implementation. implementation. 1. INTRODUCTION the intention to minimize the algorithmic and computa1. INTRODUCTION INTRODUCTION the intention intention to minimize minimize the the algorithmic andconditions. computa1. the to the and tional complexity for solving optimality 1. INTRODUCTION the intention to minimize the algorithmic algorithmic and computacomputational complexity for solving the optimality conditions. tional complexity for solving the optimality conditions. The is based on Pontryagin’s Maximum Princitionalalgorithm complexity for solving the optimality conditions. Theand algorithm is based based on onforward/backward Pontryagin’s Maximum Maximum PrinciModel predictive control nowadays can be considered The algorithm is Pontryagin’s Principle uses a combined integration of The algorithm is based on Pontryagin’s Maximum PrinciModel predictive control control nowadays can for be constrained considered ple ple and and usesdynamics a combined combined forward/backward integration of Model predictive can as an established control nowadays methodology uses a forward/backward integration of Model predictive control nowadays can be be considered considered the system and the adjoint system. Convergence ple and uses a combined forward/backward integration of as an established control methodology for constrained the system dynamics and the adjoint system. Convergence as an established control methodology for constrained dynamical systems. control Nevertheless, and despite increase the dynamics and adjoint Convergence as an established methodology for the constrained results as well as stability results forsystem. a suboptimal MPC the system system dynamics and the the adjoint system. Convergence dynamical systems. Nevertheless, and despite the increase results as as well as as stability stability results for aa suboptimal suboptimal MPC dynamical systems. Nevertheless, in computational power over the and last despite decades,the theincrease online results results for dynamical systems. Nevertheless, and despite the increase implementation provided in Graichen (2012). MPC results as well well asare stability results for a suboptimal MPC in computational computational power over over the last last decades, the (OCP) online implementation implementation are provided in Graichen (2012). in power the decades, the online solution of the constrained optimal control problem are provided in Graichen (2012). in computational power over the last decades, the online implementation are provided in Graichen (2012). solution of the the constrained optimal control problem (OCP) The intention of this paper is to numerically and expersolution of control problem (OCP) that underlies the MPC optimal formalism is still challenging solution of the constrained constrained optimal control problem (OCP) The The intention intention of this this paper is to to numerically and scheme experof is and that underlies the MPC formalism is still challenging validate thepaper fixed-point iterative MPC that underlies the MPC formalism is still challenging The intention of this paper is to numerically numerically and experexperfor systems or if sampling in the imentally thathigh-dimensional underlies the MPC formalism is still times challenging imentally validate the fixed-point iterative MPC scheme imentally validate the fixed-point iterative MPC scheme for high-dimensional systems or if sampling times in the in addition to the theoretical results in Graichen (2012) for high-dimensional systems or if sampling times in the imentally validate the fixed-point iterative MPC scheme (sub)millisecond range are concerned. for high-dimensional systems or if sampling times in the in in addition addition to the the theoretical results in in Graichen (2012) to results Graichen (sub)millisecond range range are are concerned. concerned. and to demonstrate the computational efficiency and(2012) prac(sub)millisecond in addition to the theoretical theoretical results in Graichen (2012) (sub)millisecond range are concerned. and to demonstrate the computational efficiency and pracSeveral approaches exist in the literature to enhance the and to demonstrate the computational efficiency and practical realizability. A nonlinear chemical reactor is used and to demonstrate the computational efficiency and pracSeveral approaches approaches exist in in the literature to enhance enhance the tical tical realizability. A nonlinear chemical reactor is used Several exist the literature to the computational efficiency in linear and nonlinear MPC. realizability. A nonlinear chemical reactor is used Several approaches exist in the literature to enhance the as simulation example to illustrate the performance and tical realizability. A nonlinear chemical reactor is used computational efficiency in MPC linear (EMPC) and nonlinear nonlinear MPC. as simulation example to illustrate the performance and computational efficiency in linear and MPC. For linear systems, explicit provides the as simulation example to illustrate the performance and computational efficiency in linear and nonlinear MPC. runtime requirements. In illustrate addition, the measurement results as simulation example to performance and For linear linear systems, explicit explicit MPC (EMPC) (EMPC) provides the runtime runtime requirements. Insystem addition, measurement results For systems, MPC provides option to pre-solve the linear-quadratic OCP over the requirements. In addition, measurement results For linear systems, explicit MPC (EMPC) provides the for a magnetic levitation with a sampling frequency runtime requirements. In addition, measurement results option state to pre-solve pre-solve the linear-quadratic linear-quadratic OCPminimizes over the the for aa magnetic magnetic levitation system with aa sampling sampling frequency option the OCP whole space Bemporad et al. (2002). This system with frequency option to to pre-solve the linear-quadratic OCP over over the for of Hz show levitation the practical applicability of the fixed-point for700 a magnetic levitation system with a sampling frequency whole state effort, space Bemporad Bemporad et al. al. (2002). This minimizes of 700 Hz show the practical applicability of the fixed-point whole state space et (2002). This minimizes the online but restricts the use of explicit MPC of 700 Hz show the practical applicability of the fixed-point whole state space Bemporad et al. (2002). This minimizes iteration scheme onpractical a standard dSPACE hardware. The of 700 Hz show the applicability of the fixed-point the low-dimensional online effort, effort, but restricts the use use of of explicit explicit MPC iteration scheme on onasaafollows: standard dSPACE hardware. The the restricts the MPC to (linear) systems. issue can be iteration scheme standard dSPACE hardware. The the online online effort, but but restricts the useThis of explicit MPC paper is organized Section 2 contains the probiteration scheme on a standard dSPACE hardware. The to low-dimensional low-dimensional (linear) systems. This(Ferreau issue can can be paper paper is organized organized as follows: follows: Section 22iteration contains scheme the probprobto (linear) systems. This be relaxed, e.g., by using active set strategies et al., is as Section contains the to low-dimensional (linear) systems. This issue issue can be lem statement, before the fixed-point is paper is organized as follows: Section 2 contains the probrelaxed, e.g., by using using active set set strategies strategies (Ferreau et al., al., lem statement, statement, before the the fixed-point iteration scheme is is relaxed, by active (Ferreau et 2008) or e.g., choosing a compromise online and offline before fixed-point iteration scheme relaxed, e.g., by using active set between strategies (Ferreau et al., lem summarized in Section 3. Section 4 presents the numerical lem statement, before the fixed-point iteration scheme is 2008) or or choosing choosing a compromise compromise between online online and and offline offline summarized summarized in Section 3. Section 4 presents the numerical 2008) a between computations (Zeilinger et al., 2011). in Section 3. Section 4 presents the numerical 2008) or choosing a compromise between online and offline and experimental results for the fixed-point iterative MPC summarized in Section 3. Section 4 presents the numerical computations (Zeilinger (Zeilinger et et al., al., 2011). 2011). and experimental experimental results for for the fixed-point iterative MPC computations results iterative computations (Zeilinger et al., 2011). MPC schemes often and scheme. Finally, conclusions are fixed-point drawn in Section 5. MPC and experimental results for the the fixed-point iterative MPC In case of nonlinear systems, real-time scheme. Finally, Finally, conclusions conclusions are are drawn drawn in in Section Section 5. 5. scheme. In case of nonlinear systems, real-time MPC schemes often scheme. Finally, conclusions are drawn in Section 5. In case of nonlinear systems, real-time MPC schemes often use tailored algorithms and/or suboptimal optimization. In case of nonlinear systems, real-time MPC schemes often 2. PROBLEM STATEMENT use tailored tailored algorithms and/or suboptimal optimization. use suboptimal optimization. For instance,algorithms a real-timeand/or iteration scheme based on muluse tailored algorithms and/or suboptimal optimization. 2. PROBLEM PROBLEM STATEMENT 2. For instance, a real-time iteration scheme based on mul2. PROBLEM STATEMENT STATEMENT For instance, a real-time iteration scheme based on multiple shootingaand a one-step Newton-type algorithm is For instance, real-time iteration scheme based on mulWe consider a control-affine system of the form tiple shooting and a one-step Newton-type algorithm is tiple shooting and a one-step Newton-type algorithm is presented in Diehl et al. (2005); Houska et al. (2011). tiple shooting and a one-step Newton-type algorithm is We We consider consider aa control-affine control-affinemsystem system of of the form form presented in effect Diehlofet et al. al. (2005); (2005); Houska Houska et et al. (2011). system of the the form presented in (2011). The general in MPC is We consider a control-affine m presented in Diehl Diehl etincomplete al. (2005);optimization Houska et al. al. (2011). m x˙ = f (x, u) = f0 (x) + x(t0 ) = x0 (1) The general general effect effect of incomplete incomplete optimization inand MPC is m gi (x)ui , The of is investigated Graichen and Kugioptimization (2010); Gr¨ unein PanThe general in effect of incomplete optimization in MPC MPC is x˙˙ = = ff (x, (x, u) u) = = ff00 (x) + ggii (x)u x(t )= x (1) ii ,, 0 0 x (x) + (x)u x(t = x (1) 0) 0 investigated in Graichen and Kugi (2010); Gr¨ u ne and Pani=1 x ˙ = f (x, u) = f (x) + g (x)u , x(t ) = x (1) investigated in Graichen and Kugi (2010); Gr¨ u ne and Pan0 i i 0 0 nek (2010). Several concepts solvingGr¨ the investigated in Graichen and aim Kugiat(2010); uneoptimality and Pani=1 i=1 n T m nek (2010). Several concepts aim at solving the optimality nek (2010). Several concepts aim at solving the optimality i=1 with state x ∈ R and control u = [u , . . . , u ] ∈ R conditions of the underlying OCP in a real-time feasible nek (2010). Several concepts aim at solving the optimality with state x ∈ Rnn and control u = [u1 , . . . , um ]T T ∈ Rm m conditions of the underlying OCP in a real-time feasible 1 m n T with state x ∈ R and control u = [u , . . . , u ] ∈ R conditions the OCP real-time feasible 1 , . . . , u m ] ∈ Rm subject to the valued) constraints manner. Forof the continuation in Ohtsuka state x ∈(vector R and controlcontrol u = [u conditions ofinstance, the underlying underlying OCP in in aamethod real-time feasible with 1 m subject to the (vector valued) control constraints manner. For instance, the continuation method in Ohtsuka subject to the (vector valued) control constraints − + manner. For instance, the continuation method in Ohtsuka (2004) traces the solution the discretized optimality con- subject to the (vector uvalued) ∈ [u− ,control u ] . constraints (2) manner. For instance, the of continuation method in Ohtsuka − , u+ + (2004) traces traces the single solution of the the discretized optimality con(2) (2004) the of discretized con−n + ]] .. un∈ ∈→[u [uR ,u uare (2) ditions over the MPC steps, whereasoptimality the approach (2004) traces the solution solution of the discretized optimality con- The functions f0 , gi : Ru u ∈ [u , ] . (2) to be continditions over over the(2008) singleisMPC MPC steps, whereas the theMaximum approach The functions f , g : Rnn → Rnn are assumed ditions the single steps, whereas approach Cannon et al. based on Pontryagin’s assumed to be continditions over the single MPC steps, whereas the approach The 0 ii :in ntheir n are assumed functions f , g R → R to be contin0 uously differentiable arguments. A typical control functions f0 , gi : R → R are assumed to be continCannon et et al. al.linear (2008)(discrete-time) is based based on Pontryagin’s Pontryagin’s Maximum Cannon Principle systems. InMaximum Graichen The uously arguments. typical control Cannon etfor al. (2008) (2008) is is based on on Pontryagin’s Maximum uously differentiable in their arguments. A typical task is,differentiable for instance, in thetheir stabilization of A a stationary setuously differentiable in their arguments. A typical control control Principle for linear (discrete-time) systems. In Graichen Principle for linear (discrete-time) systems. In Graichen and K¨ apernick (2012), a gradient algorithm to task is, for instance, the stabilization of a stationary setPrinciple for linear (discrete-time) systems. Intailored Graichen task is, for instance, the stabilization of a stationary setpoint task is, for instance, the stabilization of a stationary setand K¨ K¨ pernick (2012), gradient algorithm algorithm tailored to to point and aaapernick (2012), aaa gradient nonlinear MPC is described. and K¨ pernick (2012), gradient algorithm tailored tailored to point (x , u ) : f (x , u ) = 0 . (3) SP SP SP SP point nonlinear MPC is described. nonlinear MPC is (x = 0 .. (3) nonlinear MPC is described. described. SP ,, u SP ) SP ,, u SP ) (x u ) :: fff (x (x uSP = In a similar spirit as the last references, a fixed-point The MPC scheme in))this relies (3) on (xSP uSP (xSP = 00paper . (3) SP ,that SP )is: considered SP , uSP In aa similar similar spirit aspresented the last last references, references, fixed-point The MPCformulation scheme that that is considered considered in this this paper relies relies on In spirit fixed-point MPC iteration scheme wasas Graichenaaa (2012) with The the OCP In a similar spirit as the the last in references, fixed-point The MPC scheme scheme that is is considered in in this paper paper relies on on iteration scheme scheme was presented presented in Graichen Graichen (2012) with with the OCP formulation iteration iteration scheme was was presented in in Graichen (2012) (2012) with the the OCP OCP formulation formulation Copyright IFAC 2015 62 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2015, IFAC (International Federation of Automatic Control) Copyright IFAC 2015 62 Copyright © IFAC 2015 62 Peer review© of International Federation of Automatic Copyright ©under IFAC responsibility 2015 62 Control. 10.1016/j.ifacol.2015.11.060
62 62 62 62
IFAC CAO 2015 October 6-9, 2015. Garmisch-Partenkirchen, Knut Graichen et al. / IFAC-PapersOnLine 48-25 (2015) 062–067 Germany
min
J(xk , u ¯) = V (¯ x(T )) +
T
l(¯ x(τ ), u ¯(τ )) dτ
Hi
(4a)
0
x ¯˙ (τ ) = f (¯ x(τ ), u ¯(τ )) , x ¯(0) = xk = x(tk ) (4b) − + u ¯(τ ) ∈ [u , u ] , τ ∈ [0, T ] . (4c) with the MPC-internal variables x ¯, u ¯, the prediction time coordinate τ ∈ [0, T ], and the horizon length T . The initial value xk denotes the measured (or estimated) state at time step tk = t0 +k∆t, where 0 < ∆t ≤ T is the sampling time of the process. The integral and terminal cost functions in (4a) are usually designed as positive definite functions + n l : Rn × Rm → R+ 0 and V : R → R0 . Moreover, for the fixed-point iteration scheme in Section 3, the integral cost function is constructed in the separable form m 1 l(x, u) = l0 (x) + ri (ui − uSP,i )2 , ri > 0 . (5) 2 i=1 s.t.
63
control interval
u0i
u+ i
ui
Fig. 1. Illustration of minimizing control function (11). solution (6), Pontryagin’s Maximum Principle states that ¯ ∗ (τ ), τ ∈ there exists a corresponding adjoint trajectory λ k [0, T ] such that the canonical equations (9) hold and the ∗ optimal control u ¯k (τ ) satisfies ∗ ¯ ∗ (τ ), u) x∗k (τ ), λ (10) u ¯k (τ ) = arg minu∈[u− ,u+ ] H(¯ k for all times τ ∈ [0, T ]. In order to use the optimality conditions for computing the optimal solution (6), note that the particular structure of the integral cost (5) allows one to split (10) into the separate minimization problems ri Hi (x, λ, ui ) := (ui −uSP,i )2 +λT gi (x)ui min 2 ui ∈[u− ,u+ ] i i
In MPC, it is often assumed that the optimal solution of the OCP (4) at time step tk u ¯∗k (τ ) := u ¯∗ (τ ; xk ) , (6) ∗ ¯∗ (τ ; xk , u ¯∗k ) , τ ∈ [0, T ] x ¯k (τ ) := x ¯∗k ) is instantawith the optimal cost J ∗ (xk ) := J(xk , u neously available and is used as control update for the system (1) on the time interval [tk , tk+1 ), i.e. u(tk + τ ) = u ¯∗k (τ ) , tk = t0 + k∆t , τ ∈ [0, ∆t) . (7) The OCP (4) is then re-solved in the next time step tk+1 = tk + ∆t with the new system state xk+1 = x(tk+1 ).
for the single controls ui , i = 1, . . . , m. Since Hi is strictly convex in ui for ri > 0, the minimizing control function is given by 0 − + 0 ui if ui ∈ (ui , ui ) − − 0 ui = hi (x, λ) = ui if ui ≤ ui (11) + + 0 ui if ui ≥ ui with the unconstrained control 1 u0i = uSP,i − λT gi (x) , i = 1, . . . , m . ri The case-dependent definition of (11) is illustrated in Fig. 1. Inserting the optimal control function h1 (x, λ) .. u = h(x, λ) = (12) . hm (x, λ)
The endpoint-free formulation of the OCP (4) is chosen on purpose and will be of advantage for the fixedpoint iteration scheme in Section 3. Moreover, stability of an MPC scheme without terminal constraints can still be guaranteed under certain assumptions. In most cases, it is assumed that the terminal cost V (x) represents a control Lyapunov function (CLF), see e.g. Parisini and Zoppoli (1995); Chen and Allg¨ ower (1998). In the absence of the terminal cost V (x), the prediction time T must be sufficiently long in order to guarantee closed-loop stability (Jadbabaie and Hauser, 2005; Grimm et al., 2005; Gr¨ une and Pannek, 2011). 3. FIXED-POINT ITERATION SCHEME The particular structure of the OCP (4) with the integral cost (5) can be exploited by solving the respective optimality conditions in a fixed-point iterative manner. This section briefly summarizes the concept of the fixedpoint iteration scheme and derives the algorithm. Details as well as convergence and stability results can be found in Graichen (2012).
into (9) leads to the new form of the canonical equations x˙ = F (x, λ) , x(0) = xk (13a) (13b) λ˙ = G(x, λ) , λ(T ) = Vx (x(T )) with the new functions F (x, λ) := f (x, h(x, λ)) and G(x, λ) := −Hx (x, λ, h(x, λ)). The boundary value problem (13) is independent of the control u and consists of 2n ordinary differential equations with n initial conditions for the state x and n terminal conditions for the adjoint state λ.
3.1 Optimality conditions The fixed-point iteration scheme is based on the optimality conditions of the OCP (4). With the Hamiltonian H(x, λ, u) = l(x, u) + λT f (x, u) , the canonical equations follow as x˙ = f (x, u) , x(0) = xk ˙λ = −Hx (x, λ, u) , λ(T ) = Vx (x(T ))
u− i
3.2 Algorithm
(8)
The separated boundary conditions for x and λ can be taken advantage of by solving (13a) and (13b) successively forward and backward in time. Table 1 summarizes the algorithm, which basically consists of two numerical integrations. The additional damping steps involve the damping factor ε ∈ (0, 1] and the intermediate trajectories
(9a)
(9b) with the adjoint state λ ∈ Rn and the partial derivative functions Hx and Vx . If the OCP (4) has an optimal 63
IFAC CAO 2015 64 Knut Graichen et al. / IFAC-PapersOnLine 48-25 (2015) 062–067 October 6-9, 2015. Garmisch-Partenkirchen, Germany
Table 1. Damped fixed-point iteration scheme.
iterations, N , in each MPC step k. This leads to the (N ) ¯ (N ) (τ ), τ ∈ [0, T ] as well suboptimal trajectories x ¯k (τ ), λ k as the corresponding suboptimal control trajectory (12) (N ) (N ) ¯ (N ) (τ ) , τ ∈ [0, T ] . u ¯k (τ ) = h x ¯k (τ ), λ (16) k
1) Initialization • Damping factor ε ∈ (0, 1] (0) • Initial guess x ¯k (τ ), τ ∈ [0, T ], e.g. using (14) (0) ¯ • Initial guess λk (τ ), τ ∈ [0, T ], e.g. using (15)
Thus, the optimal feedback in (7) for the system (1) is replaced by the suboptimal one
2) Damped fixed-point iteration (1 ≤ j ≤ N ) • Forward integration of (j) x ˆ˙ k (τ ) = F
(j)
(j−1)
¯ x ˆk (τ ), λ k
• Damping of state trajectory (j)
(j)
(j)
(τ ) , x ˆk (0) = xk (j−1)
xk (τ ) + (1 − ε)¯ xk x ¯k (τ ) = εˆ
• Backward integration of
(N )
¯k (τ ) , τ ∈ [0, ∆t) . (17) u(tk + τ ) = u In the next MPC step k + 1, the damped fixed-point iteration scheme can be initialized with the previous trajectories (0) (N ) ¯ (0) (τ ) = λ ¯ (N ) (τ ), τ ∈ [0, T ] . (18) (τ ) = x ¯ (τ ), λ x ¯
(τ ) ,
τ ∈ [0, T ]
k+1
(j) (j) ˆ˙ (j) (τ ) = G x ˆ (j) (τ ) , λ ˆ (j) (T ) = Vx (¯ λ ¯k (τ ), λ xk (T )) k k k
k+1
k
The stability properties of the suboptimal MPC scheme with a constant number N of fixed-point iterations are investigated in Graichen (2012) for the undamped case (ε = 1). By defining the optimization error (N ) (N ) ∆¯ ¯k (·) − x ¯∗k (·) Ln [0,T ] (19) xk (·) Ln [0,T ] := x
• Damping of adjoint state trajectory
ˆ (j) (τ ) + (1 − ε)λ ¯ (j−1) (τ ) , ¯ (j) (τ ) = ελ λ k k k
k
Hence, the integration of (14) and (15) as stated in Table 1 is only required in the first MPC step k = 0.
τ ∈ [0, T ]
• Stop if suitable convergence criterion is fulfilled or if j = N . Otherwise, set j ← j +1 and return to 2).
∞
(j) ˆ (j) (τ ), τ ∈ [0, T ]. The fixed-point iterations x ˆk (τ ) and λ k are undamped for ε = 1, whereas strong damping is achieved for ε → 0. It is shown in Graichen (2012) that for a given set of initial values (xk ∈ X ⊂ Rn ) an upper horizon time Tmax > 0 exists such that for all T < Tmax the undamped fixed-point algorithm (ε = 1) has a linear rate of convergence, i.e. (j) (j−1) ¯ ∗ (·) n ¯ ∗ (·) n ¯ λk (·) − λ ≤ p ¯ λk (·) − λ k k L∞ [0,T ] L [0,T ] (j−1) ∞ (j) ∗ ∗ x ¯k (·) Ln [0,T ] ≤ p x ¯k (·) − x ¯k (·) Ln [0,T ] ¯k (·) − x ∞
(i) (ii) (iii)
∞
The fixed-point iterative MPC scheme is numerically and experimentally validated in order to show the computational efficiency as well as the practical applicability. 4.1 Simulation example: chemical reactor The following nonlinear model describes the dynamics of a continuous stirred tank reactor (CSTR) c˙A = −k1 (ϑ)cA − k2 (ϑ)c2A + [cin − cA ]u1
with the (constant) control input uSP of the setpoint (3). The initial trajectory of the adjoint state λ then follows from the backward integration of the adjoint system (13b) (0) (0) ¯˙ (0)(τ ) = G x ¯ (0)(τ ) , λ ¯ (0)(T ) = Vx (¯ λ ¯ (τ ), λ xk (T )) . (15) k k k k
(20a)
(20b) c˙B = k1 (ϑ)cA − k1 (ϑ)cB − cB u1 ˙ ϑ = −δ k1 (ϑ)cA ∆HAB + k1 (ϑ)cB ∆HBC (20c) 2 +k2 (ϑ)cA ∆HAD + α[ϑc − ϑ] + [ϑin − ϑ]u1 ˙ (20d) ϑc = β[ϑ − ϑc ] + γu2
3.3 MPC implementation
k
∞
4. NUMERICAL AND EXPERIMENTAL VALIDATION
The start of the fixed-point iteration scheme requires an (0) ¯ (0) (τ ), τ ∈ initial guess of the trajectories x ¯k (τ ) and λ k (0) [0, T ]. One possibility is to compute x ¯k (τ ) by forward integration of the system dynamics (1) (0) (0) (0) x ¯˙ k (τ ) = f x ¯k (0) = xk ¯k (τ ), uSP , x (14)
k
the fixed-point iterations are convergent (T < Tmax ), N satisfies a lower bound, (0) ∆¯ x0 (·) Ln [0,T ] satisfies an upper bound,
it is shown in Graichen (2012) that asymptotic stability as well as incremental reduction of the optimization error (19) is ensured.
for some p ∈ (0, 1). In practice, the maximum horizon length Tmax for which convergence of the algorithm is ensured can be enlarged by introducing damping as it is done in Table 1. The influence of the damping factor ε ∈ (0, 1] on the convergence is illustrated numerically for a simulation example in Section 4.1.
Typically, the algorithm in Table 1 would iterate until sufficient accuracy is achieved by choosing a suitable convergence criterion. 1 However, in order to implement the algorithm within a real-time MPC framework, it is more convenient to use a constant number of fixed-point (j) 1 For instance, λ ¯ (j−1) (·) ¯ (·) − λ ≤ δλ for some δλ > 0. n
∞
as a measure for the suboptimality in MPC step k and assuming that
with the Arrhenius functions −Ei ki (T ) = ki0 exp , i = 1, 2 . T /◦ C + 273.15 Details on the model in (20) as well as the model parameters can e.g. be found in Klatt and Engell (1998); Graichen et al. (2009). The state vector x = [cA , cB , ϑ, ϑc ]T comprises the educt and product concentration as well as the reactor and coolant temperature, respectively. The control
L∞ [0,T ]
64
IFAC CAO 2015 October 6-9, 2015. Garmisch-Partenkirchen, Knut Graichen et al. / IFAC-PapersOnLine 48-25 (2015) 062–067 Germany
vector u = [u1 , u2 ]T consists of the normalized input flow rate and the cooling power, subject to the constraints MJ (21) u1 ∈ [5 h−1 , 35 h−1 ] , u2 ∈ [−8.5 MJ h ,0 h ]. The control goal that is used to evaluate the fixed-point iterative MPC scheme is to maximize the product concentration cB while maintaining a desired setpoint ϑSP for the reactor temperature. To this end, the cost functional in (4a) is chosen as 2 T 2 ¯ J(xk , u ¯) = −qB c¯2B +qϑ (ϑ−ϑ ¯21 +r2 u ¯22 dτ (22) SP ) +r1 u
65
16 14
Tmax [min]
12 10 8 6 4
0
6
m 2◦ 2 C h, with the coefficients qB = 102 kmol 2 h , qϑ = 10 −1 0 h r1 = 10 h, and r2 = 10 MJ2 . The sampling time for the MPC scheme is set to ∆t = 1 s. The initial state of the reactor is the stationary setpoint (Utz et al., 2007) T ◦ ◦ kmol x0 = 3.5175 kmol . (23) m3 , 0.74 m3 , 87 C , 79.8 C The damped fixed-point iteration scheme in Table 1 is implemented as C code with a Cmex interface for Matlab. The canonical equations (13) including the optimal control function (12) are computed under Mathematica and are exported as optimized Cmex functions. The forward and backward integrations in the algorithm are performed using a Heun integration scheme on an uniform discretization of the MPC interval [0, T ] with Ndiscr points.
2 0 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
ε ∈ (0, 1]
Fig. 2. Maximum horizon length Tmax for convergence of the fixed-point iterations w.r.t. the damping factor ε ∈ (0, 1]. CPU (M620, 2.67 GHz) and Matlab R2012a (Windows 7, 64 bit). The damping factor and horizon length are set to ε = 0.1 and T = 10 min according to the simulation results in Fig. 3. The last column in Table 2 shows the overall cost value Tsim Jsim = ¯22 dτ (24) −qB c¯2B + qϑ (ϑ¯ − ϑdes )2 + r1 u ¯21 + r2 u
In the first step, it is investigated how the damping factor ε ∈ (0, 1] of the fixed-point iteration scheme affects the maximum horizon length Tmax for which convergence is achieved. To this end, T is successively increased for a constant ε until the MPC scheme in closed-loop turns unstable. The resulting ε-T -curve is depicted in Fig. 2. In the undamped case (ε = 1), the maximum horizon length is Tmax = 0.7 min, whereas for increased damping of the fixed-point iterations Tmax can be significantly enlarged.
0
that is obtained by integration of the simulated trajectories over the simulation time Tsim = 70 min. The value of Jsim can be seen as an indicator for the control performance that decreases for larger values of Ndiscr and N . The statistics in Table 2 as well as the trajectories in Fig. 3 show that an excellent performance is already obtained for N = 2 fixed-point iterations per MPC step with a computation time of 32 µs.
Figure 3 shows the simulation results for the CSTR with N = 2 fixed-point iterations per MPC step and the horizon length T = 10 min, which is sufficiently short for the chosen damping factor ε = 0.1 (cf. Fig 2). In the first time interval [0, 10) min, the MPC is inactive and the reactor stays at the initial setpoint (23). In the second time interval [10, 40) min, the MPC is activated with the desired temperature ϑSP = 87◦ C corresponding to the initial setpoint x0 . By minimizing the cost functional (22), the MPC increases the product concentration cB while maintaining the reactor temperature ϑ in the vicinity of the desired value ϑSP = 87◦ C. During the third time interval [40, 70) min, the desired reactor temperature is increased to the setpoint ϑSP = 110◦ C which also results in a larger product concentration. During both transitions, the controls act in an aggressive manner owing to the comparatively small control weights in the cost functional (22).
Table 2. MPC results for the CSTR. Discretization points Ndiscr
Fixed-point iterations N
CPU time [µs]
Simulated cost Jsim [-]
15
2 5 10
32 55 96
17.56 17.47 17.44
30
2 5 10
58 98 183
15.93 15.87 15.86
60
2 5 10
83 186 355
15.78 15.55 15.54
4.2 Experiment: magnetic levitation system Figure 4 shows the schematic drawing of a magnetic levitation system with an electromagnet and a levitating ferromagnetic object with mass m = 133 g. The electromagnet consists of two coils (2000 windings each) on a steel cup core to generate the magnetic field. The electromagnetic force is applied to the levitating mass in the opposite direction of the gravitational force. The control input is the duty cycle τ ∈ [0, 1] of the pulse width modulated (PWM) voltages applied to each coil. The measured variables are the currents i1 and i2 of both coils, the temperature inside the coils, and the levitation height z by means of a laser
In order to evaluate the computational efficiency in addition to the simulation results in Fig. 3, the runtime requirements are recorded for different numbers of discretization points Ndiscr and fixed-point iterations N per MPC step. Table 2 summarizes the runtime results for the Cmex compiled MPC scheme on a computer with Intel Core i7 2
Note that the goal to maximize the product concentration leads to an integral cost function that is not positive definite. This setting somewhat differs from the classical stabilization task and belongs to the framework of economic MPC, see e.g. (Ferramosca et al., 2010).
65
IFAC CAO 2015 66 Knut Graichen et al. / IFAC-PapersOnLine 48-25 (2015) 062–067 October 6-9, 2015. Garmisch-Partenkirchen, Germany
110
3 2.5 2
1.2
105
ϑ [◦ C]
cB [kmol/ m3 ]
cA [kmol/ m3 ]
4 3.5
1 0.8
95 90
0.6
1.5 10
20
30
40
50
60
70
10
20
time [min]
30
40
50
60
85
70
95 90
25 20 15
85
10
80
5 30
40
50
60
u 2 [MJ/h]
u 1 [1/ h]
100
70
30
40
50
60
70
50
60
70
0
30
20
20
time [min]
35
105
10
10
time [min]
110
ϑ c [◦ C]
100
−2 −4 −6 −8
0
10
20
time [s]
30
40
50
60
70
0
10
20
time [min]
30
40
time [min]
Fig. 3. MPC simulation results for the CSTR. for different values of zSP . To this end, the cost functional in (4a) is chosen as T J(xk , u ¯) = qz (¯ z − zSP )2 + qz˙ z¯˙ 2 + r (¯ τ − τSP )2 dτ (25)
sensor. The experimental setup is constructed for large levitation heights of up to 7 cm which necessitates an active air cooling inside the coils. Moreover, no cascaded current controllers are used in order to achieve the fastest possible control performance.
0
with the weights qz = 104 m12 s , qz˙ = 25 ms2 , and r = 0.25 1s . The stationary control input (duty cycle) is set to τSP = 0.5 for the sake of simplicity. An additional integral gain is added to the MPC scheme in order to compensate for steady state errors in the position of the levitating mass. 1 The MPC sampling time ∆t = 700 Hz ≈ 1.4 ms corresponds to the PWM frequency. The prediction horizon and the damping factor for the fixed-point iterations are set to T = 70 ms and ε = 0.1, respectively.
In B¨ achle et al. (2013), a nonlinear model of the magnetic levitation systen is derived and an alternative gradientbased MPC scheme is used to control the system. In the following, the nonlinear model of this levitation system with the state x = [i1 , i2 , z, z] ˙ T and control τ ∈ [0, 1] is used to demonstrate the practical applicability of the fixed-point iteration scheme. The control task is to stabilize a desired setpoint of the levitation height zSP and to perform fast setpoint changes
The MPC scheme uses the Matlab/C implementation of the fixed-point iteration scheme described in Section 4.1. Table 3 shows the runtime results for the MPC implementation using Ndiscr = 30 discretization points for the prediction horizon [0, T ]. For the experimental implementation, a dSPACE MicroAutoBox (800 MHz) is used as real-time hardware. The computation time on the dSPACE box amounts to approximately 400 µs including the MPC scheme and state estimation (Unscented Kalman Filter), which is well below the sampling time.
air gaps for active cooling coil 1
Figure 5 shows the measurement results for different setpoint changes of the levitation height z between −70 mm and −45 mm. Overall, the MPC scheme shows a very good control performance. The noise in the control inputs is mainly due to the noisy measurements and filtered signals that leads to persistent perturbations of the MPC scheme through the initial conditions in the OCP (4). This effect also demonstrates the robustness of the fixed-point iteration scheme. In addition, the last plot of Fig. 5 depicts the CPU requirements per MPC step plotted over time
coil 2
z cup core levitating mass
Table 3. MPC results for the levitation system.
Fig. 4. Schematic of the magnetic levitation system. 66
Discretization points Ndiscr
Fixed-point iterations N
CPU time [µs]
dSPACE time [µs]
30
2
60
378
IFAC CAO 2015 October 6-9, 2015. Garmisch-Partenkirchen, Knut Graichen et al. / IFAC-PapersOnLine 48-25 (2015) 062–067 Germany
67
z [cm]
REFERENCES −5
B¨ achle, T., Hentzelt, S., and Graichen, K. (2013). Nonlinear model predictive control of a magnetic levitation system. Contr. Eng. Pract., 21(9), 1250–1258. Bemporad, A., Morari, M., Dua, V., and Pistikopoulos, E. (2002). The explicit linear quadratic regulator for constrained systems. Automatica, 38(1), 3–20. Cannon, M., Liao, W., and Kouvaritakis, B. (2008). Efficient MPC optimization using Pontryagin’s Minimum Principle. Int. J. Robust Nonl. Contr., 18, 831–844. Chen, H. and Allg¨ ower, F. (1998). A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability. Automatica, 34(10), 1205–1217. Diehl, M., Bock, H., and Schl¨ oder, J. (2005). A real-time iteration scheme for nonlinear optimization in optimal feedback control. SIAM J. Contr. Optim., 43(5), 1714–1736. Ferramosca, A., Rawlings, J., and Camacho, D.L.E. (2010). Economic MPC for a changing economic criterion. In Proc. 49th IEEE Conf. Dec. Contr. (CDC), 6131–6136. Atlanta (USA). Ferreau, H., Bock, H., and Diehl, M. (2008). An online active set strategy to overcome the limitations of explicit MPC. Int. J. Robust Nonl. Contr., 18(8), 816–830. Graichen, K. (2012). A fixed-point iteration scheme for real-time model predictive control. Automatica, 48(7), 1300–1305. Graichen, K., Hagenmeyer, V., and Zeitz, M. (2009). Design of adaptive feedforward control under input constraints for a benchmark CSTR based on a BVP solver. Comput. Chem. Eng., 33(2), 473–483. Graichen, K. and K¨ apernick, B. (2012). A real-time gradient method for nonlinear model predictive control. In Frontiers of Model Predictive Control, 9–28. InTech (open access). http://www. intechopen.com/articles/show/title/a-real-time-gradientmethod-for-nonlinear-model-predictive-control. Graichen, K. and Kugi, A. (2010). Stability and incremental improvement of suboptimal MPC without terminal constraints. IEEE Trans. Automat. Contr., 55(11), 2576–2580. Grimm, G., Messina, M., Tuna, S., and Teel, A. (2005). Model predictive control: For want of a local control Lyapunov function, all is not lost. IEEE Trans. Automat. Contr., 50(5), 546–558. Gr¨ une, L. and Pannek, J. (2010). Analysis of unconstrained NMPC schemes with incomplete optimization. In Proc. 8th IFAC Symp. Nonl. Contr. Syst. (NOLCOS), 238–243. Bologna (Italy). Gr¨ une, L. and Pannek, J. (2011). Nonlinear Model Predictive Control: Theory and Algorithms. Springer, London. Houska, B., Ferreau, H., and Diehl, M. (2011). An auto-generated real-time iteration algorithm for nonlinear MPC in the microsecond range. Automatica, 47(10), 2279–2285. Jadbabaie, A. and Hauser, J. (2005). On the stability of receding horizon control with a general terminal cost. IEEE Trans. Automat. Contr., 50(5), 674–678. Klatt, K.U. and Engell, S. (1998). Gain-scheduling trajectory control of a continuous stirred tank reactor. Contr. Eng. Pract., 22(4-5), 491–502. Ohtsuka, T. (2004). A continuation/GMRES method for fast computation of nonlinear receding horizon control. Automatica, 40(4), 563–574. Parisini, T. and Zoppoli, R. (1995). A receding–horizon regulator for nonlinear systems and a neural approximation. Automatica, 31(10), 1443–1451. Utz, T., Hagenmeyer, V., and Mahn, B. (2007). Comparative evaluation of nonlinear model predictive and flatness-based twodegree-of-freedom control design in view of industrial application. J. Process Contr., 17(2), 129–141. Zeilinger, M., Jones, C., and Morari, M. (2011). Real-time suboptimal model predictive control using a combination of explicit MPC and online optimization. IEEE Trans. Automat. Contr., 56(7), 1524–1534.
−6 −7 0
2
4
6 time [s]
8
12
coil 1 coil 2
6 i 1 & i 2 [A]
10
4 2 0
2
4
6 time [s]
8
10
12
0
2
4
6 time [s]
8
10
12
0
2
4
6 time [s]
8
10
12
τ [-]
1 0.5
CPU time [ms]
0
1
0.5
0
Fig. 5. MPC measurement results for the magnetic levitation system. (also cf. Table 3). The computation time is practically constant, which illustrates the real-time feasibility of the MPC scheme based on a constant number of fixed-point iterations per sampling step ∆t. A video of the MPC-controlled magnetic levitation system is available via the link http://youtu.be/01q4TUVrToc. 5. CONCLUSIONS The numerical and experimental results in the paper demonstrate the low complexity and computational cost of the fixed-point iteration scheme recently presented for MPC of nonlinear systems subject to control constraints. In general, the horizon length cannot be chosen arbitrarily large in order to ensure convergence as well as stability in closed-loop. This can be interpreted as a payoff for the low algorithmic complexity of the fixed-point iteration scheme. On the other hand, the introduction of damping in the algorithm usually allows one to significantly enlarge the horizon length, as demonstrated for the chemical reactor example. A key role for tuning the fixed-point iterative MPC scheme plays the choice of the damping factor in relation to the prediction horizon and the number of fixedpoint iterations per MPC step. Numerical studies have shown that two iterations per MPC step are often sufficient to obtain a good control performance. 67
ScienceDirect IFAC-PapersOnLine 48-25 (2015) 062–067
Nonlinear MPC with fixed-point iterations: Nonlinear MPC with fixed-point iterations: Nonlinear MPC with fixed-point iterations: numerical and experimental validation experimental validation numerical and numerical and experimental validation Knut Graichen Knut Graichen Knut Knut Graichen Graichen Institute of Measurement, Control and Microtechnology, Institute of Measurement, Control and Microtechnology, Institute of Control and University of Ulm, Germany (e-mail: [email protected]). Institute of Measurement, Measurement, Control and Microtechnology, Microtechnology, University of Ulm, Germany (e-mail: [email protected]). University of Ulm, Germany (e-mail: [email protected]). University of Ulm, Germany (e-mail: [email protected]). Abstract: The paper demonstrates the applicability of a recently presented MPC scheme for Abstract: The paper the applicability of a recently presented MPC scheme for Abstract: The demonstrates the of presented MPC for nonlinear withdemonstrates control constraints. The algorithm is based on fixed-point iterations Abstract:systems The paper paper demonstrates the applicability applicability of aa recently recently presented MPC scheme scheme for nonlinear systems with control constraints. The algorithm is based on fixed-point iterations nonlinear systemsconsist with control control constraints. The algorithm algorithm is based based on on fixed-point iterations that in principal of combined forward/backward integrations over the MPCiterations horizon. nonlinear systems with constraints. The is fixed-point that in principal consist of combined forward/backward integrations the MPC horizon. that consist of forward/backward integrations over the MPC horizon. Numerical and runtime for a nonlinear chemical reactor modelover show efficiency of that in in principal principal consistresults of combined combined forward/backward integrations over thethe MPC horizon. Numerical and runtime results for aa nonlinear chemical reactor model show the efficiency of Numerical and runtime results for nonlinear chemical reactor model show the efficiency of the fixed-point iterative MPC scheme. Experimental results for a magnetic levitation system Numerical and runtime results for a nonlinear chemical reactor model show the efficiency of the fixed-point iterative MPC scheme. Experimental results for aa magnetic levitation system the fixed-point iterative MPC scheme. Experimental results for magnetic levitation system confirm its practical applicability. the fixed-point iterative MPC scheme. Experimental results for a magnetic levitation system confirm confirm its its practical practical applicability. applicability. confirm its practical applicability. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Model predictive control, fixed-point iterations, control constraints, real-time Keywords: Model Model predictive control, control, fixed-point iterations, iterations, control constraints, constraints, real-time Keywords: implementation. Keywords: Model predictive predictive control, fixed-point fixed-point iterations, control control constraints, real-time real-time implementation. implementation. implementation. 1. INTRODUCTION the intention to minimize the algorithmic and computa1. INTRODUCTION INTRODUCTION the intention intention to minimize minimize the the algorithmic andconditions. computa1. the to the and tional complexity for solving optimality 1. INTRODUCTION the intention to minimize the algorithmic algorithmic and computacomputational complexity for solving the optimality conditions. tional complexity for solving the optimality conditions. The is based on Pontryagin’s Maximum Princitionalalgorithm complexity for solving the optimality conditions. Theand algorithm is based based on onforward/backward Pontryagin’s Maximum Maximum PrinciModel predictive control nowadays can be considered The algorithm is Pontryagin’s Principle uses a combined integration of The algorithm is based on Pontryagin’s Maximum PrinciModel predictive control control nowadays can for be constrained considered ple ple and and usesdynamics a combined combined forward/backward integration of Model predictive can as an established control nowadays methodology uses a forward/backward integration of Model predictive control nowadays can be be considered considered the system and the adjoint system. Convergence ple and uses a combined forward/backward integration of as an established control methodology for constrained the system dynamics and the adjoint system. Convergence as an established control methodology for constrained dynamical systems. control Nevertheless, and despite increase the dynamics and adjoint Convergence as an established methodology for the constrained results as well as stability results forsystem. a suboptimal MPC the system system dynamics and the the adjoint system. Convergence dynamical systems. Nevertheless, and despite the increase results as as well as as stability stability results for aa suboptimal suboptimal MPC dynamical systems. Nevertheless, in computational power over the and last despite decades,the theincrease online results results for dynamical systems. Nevertheless, and despite the increase implementation provided in Graichen (2012). MPC results as well well asare stability results for a suboptimal MPC in computational computational power over over the last last decades, the (OCP) online implementation implementation are provided in Graichen (2012). in power the decades, the online solution of the constrained optimal control problem are provided in Graichen (2012). in computational power over the last decades, the online implementation are provided in Graichen (2012). solution of the the constrained optimal control problem (OCP) The intention of this paper is to numerically and expersolution of control problem (OCP) that underlies the MPC optimal formalism is still challenging solution of the constrained constrained optimal control problem (OCP) The The intention intention of this this paper is to to numerically and scheme experof is and that underlies the MPC formalism is still challenging validate thepaper fixed-point iterative MPC that underlies the MPC formalism is still challenging The intention of this paper is to numerically numerically and experexperfor systems or if sampling in the imentally thathigh-dimensional underlies the MPC formalism is still times challenging imentally validate the fixed-point iterative MPC scheme imentally validate the fixed-point iterative MPC scheme for high-dimensional systems or if sampling times in the in addition to the theoretical results in Graichen (2012) for high-dimensional systems or if sampling times in the imentally validate the fixed-point iterative MPC scheme (sub)millisecond range are concerned. for high-dimensional systems or if sampling times in the in in addition addition to the the theoretical results in in Graichen (2012) to results Graichen (sub)millisecond range range are are concerned. concerned. and to demonstrate the computational efficiency and(2012) prac(sub)millisecond in addition to the theoretical theoretical results in Graichen (2012) (sub)millisecond range are concerned. and to demonstrate the computational efficiency and pracSeveral approaches exist in the literature to enhance the and to demonstrate the computational efficiency and practical realizability. A nonlinear chemical reactor is used and to demonstrate the computational efficiency and pracSeveral approaches approaches exist in in the literature to enhance enhance the tical tical realizability. A nonlinear chemical reactor is used Several exist the literature to the computational efficiency in linear and nonlinear MPC. realizability. A nonlinear chemical reactor is used Several approaches exist in the literature to enhance the as simulation example to illustrate the performance and tical realizability. A nonlinear chemical reactor is used computational efficiency in MPC linear (EMPC) and nonlinear nonlinear MPC. as simulation example to illustrate the performance and computational efficiency in linear and MPC. For linear systems, explicit provides the as simulation example to illustrate the performance and computational efficiency in linear and nonlinear MPC. runtime requirements. In illustrate addition, the measurement results as simulation example to performance and For linear linear systems, explicit explicit MPC (EMPC) (EMPC) provides the runtime runtime requirements. Insystem addition, measurement results For systems, MPC provides option to pre-solve the linear-quadratic OCP over the requirements. In addition, measurement results For linear systems, explicit MPC (EMPC) provides the for a magnetic levitation with a sampling frequency runtime requirements. In addition, measurement results option state to pre-solve pre-solve the linear-quadratic linear-quadratic OCPminimizes over the the for aa magnetic magnetic levitation system with aa sampling sampling frequency option the OCP whole space Bemporad et al. (2002). This system with frequency option to to pre-solve the linear-quadratic OCP over over the for of Hz show levitation the practical applicability of the fixed-point for700 a magnetic levitation system with a sampling frequency whole state effort, space Bemporad Bemporad et al. al. (2002). This minimizes of 700 Hz show the practical applicability of the fixed-point whole state space et (2002). This minimizes the online but restricts the use of explicit MPC of 700 Hz show the practical applicability of the fixed-point whole state space Bemporad et al. (2002). This minimizes iteration scheme onpractical a standard dSPACE hardware. The of 700 Hz show the applicability of the fixed-point the low-dimensional online effort, effort, but restricts the use use of of explicit explicit MPC iteration scheme on onasaafollows: standard dSPACE hardware. The the restricts the MPC to (linear) systems. issue can be iteration scheme standard dSPACE hardware. The the online online effort, but but restricts the useThis of explicit MPC paper is organized Section 2 contains the probiteration scheme on a standard dSPACE hardware. The to low-dimensional low-dimensional (linear) systems. This(Ferreau issue can can be paper paper is organized organized as follows: follows: Section 22iteration contains scheme the probprobto (linear) systems. This be relaxed, e.g., by using active set strategies et al., is as Section contains the to low-dimensional (linear) systems. This issue issue can be lem statement, before the fixed-point is paper is organized as follows: Section 2 contains the probrelaxed, e.g., by using using active set set strategies strategies (Ferreau et al., al., lem statement, statement, before the the fixed-point iteration scheme is is relaxed, by active (Ferreau et 2008) or e.g., choosing a compromise online and offline before fixed-point iteration scheme relaxed, e.g., by using active set between strategies (Ferreau et al., lem summarized in Section 3. Section 4 presents the numerical lem statement, before the fixed-point iteration scheme is 2008) or or choosing choosing a compromise compromise between online online and and offline offline summarized summarized in Section 3. Section 4 presents the numerical 2008) a between computations (Zeilinger et al., 2011). in Section 3. Section 4 presents the numerical 2008) or choosing a compromise between online and offline and experimental results for the fixed-point iterative MPC summarized in Section 3. Section 4 presents the numerical computations (Zeilinger (Zeilinger et et al., al., 2011). 2011). and experimental experimental results for for the fixed-point iterative MPC computations results iterative computations (Zeilinger et al., 2011). MPC schemes often and scheme. Finally, conclusions are fixed-point drawn in Section 5. MPC and experimental results for the the fixed-point iterative MPC In case of nonlinear systems, real-time scheme. Finally, Finally, conclusions conclusions are are drawn drawn in in Section Section 5. 5. scheme. In case of nonlinear systems, real-time MPC schemes often scheme. Finally, conclusions are drawn in Section 5. In case of nonlinear systems, real-time MPC schemes often use tailored algorithms and/or suboptimal optimization. In case of nonlinear systems, real-time MPC schemes often 2. PROBLEM STATEMENT use tailored tailored algorithms and/or suboptimal optimization. use suboptimal optimization. For instance,algorithms a real-timeand/or iteration scheme based on muluse tailored algorithms and/or suboptimal optimization. 2. PROBLEM PROBLEM STATEMENT 2. For instance, a real-time iteration scheme based on mul2. PROBLEM STATEMENT STATEMENT For instance, a real-time iteration scheme based on multiple shootingaand a one-step Newton-type algorithm is For instance, real-time iteration scheme based on mulWe consider a control-affine system of the form tiple shooting and a one-step Newton-type algorithm is tiple shooting and a one-step Newton-type algorithm is presented in Diehl et al. (2005); Houska et al. (2011). tiple shooting and a one-step Newton-type algorithm is We We consider consider aa control-affine control-affinemsystem system of of the form form presented in effect Diehlofet et al. al. (2005); (2005); Houska Houska et et al. (2011). system of the the form presented in (2011). The general in MPC is We consider a control-affine m presented in Diehl Diehl etincomplete al. (2005);optimization Houska et al. al. (2011). m x˙ = f (x, u) = f0 (x) + x(t0 ) = x0 (1) The general general effect effect of incomplete incomplete optimization inand MPC is m gi (x)ui , The of is investigated Graichen and Kugioptimization (2010); Gr¨ unein PanThe general in effect of incomplete optimization in MPC MPC is x˙˙ = = ff (x, (x, u) u) = = ff00 (x) + ggii (x)u x(t )= x (1) ii ,, 0 0 x (x) + (x)u x(t = x (1) 0) 0 investigated in Graichen and Kugi (2010); Gr¨ u ne and Pani=1 x ˙ = f (x, u) = f (x) + g (x)u , x(t ) = x (1) investigated in Graichen and Kugi (2010); Gr¨ u ne and Pan0 i i 0 0 nek (2010). Several concepts solvingGr¨ the investigated in Graichen and aim Kugiat(2010); uneoptimality and Pani=1 i=1 n T m nek (2010). Several concepts aim at solving the optimality nek (2010). Several concepts aim at solving the optimality i=1 with state x ∈ R and control u = [u , . . . , u ] ∈ R conditions of the underlying OCP in a real-time feasible nek (2010). Several concepts aim at solving the optimality with state x ∈ Rnn and control u = [u1 , . . . , um ]T T ∈ Rm m conditions of the underlying OCP in a real-time feasible 1 m n T with state x ∈ R and control u = [u , . . . , u ] ∈ R conditions the OCP real-time feasible 1 , . . . , u m ] ∈ Rm subject to the valued) constraints manner. Forof the continuation in Ohtsuka state x ∈(vector R and controlcontrol u = [u conditions ofinstance, the underlying underlying OCP in in aamethod real-time feasible with 1 m subject to the (vector valued) control constraints manner. For instance, the continuation method in Ohtsuka subject to the (vector valued) control constraints − + manner. For instance, the continuation method in Ohtsuka (2004) traces the solution the discretized optimality con- subject to the (vector uvalued) ∈ [u− ,control u ] . constraints (2) manner. For instance, the of continuation method in Ohtsuka − , u+ + (2004) traces traces the single solution of the the discretized optimality con(2) (2004) the of discretized con−n + ]] .. un∈ ∈→[u [uR ,u uare (2) ditions over the MPC steps, whereasoptimality the approach (2004) traces the solution solution of the discretized optimality con- The functions f0 , gi : Ru u ∈ [u , ] . (2) to be continditions over over the(2008) singleisMPC MPC steps, whereas the theMaximum approach The functions f , g : Rnn → Rnn are assumed ditions the single steps, whereas approach Cannon et al. based on Pontryagin’s assumed to be continditions over the single MPC steps, whereas the approach The 0 ii :in ntheir n are assumed functions f , g R → R to be contin0 uously differentiable arguments. A typical control functions f0 , gi : R → R are assumed to be continCannon et et al. al.linear (2008)(discrete-time) is based based on Pontryagin’s Pontryagin’s Maximum Cannon Principle systems. InMaximum Graichen The uously arguments. typical control Cannon etfor al. (2008) (2008) is is based on on Pontryagin’s Maximum uously differentiable in their arguments. A typical task is,differentiable for instance, in thetheir stabilization of A a stationary setuously differentiable in their arguments. A typical control control Principle for linear (discrete-time) systems. In Graichen Principle for linear (discrete-time) systems. In Graichen and K¨ apernick (2012), a gradient algorithm to task is, for instance, the stabilization of a stationary setPrinciple for linear (discrete-time) systems. Intailored Graichen task is, for instance, the stabilization of a stationary setpoint task is, for instance, the stabilization of a stationary setand K¨ K¨ pernick (2012), gradient algorithm algorithm tailored to to point and aaapernick (2012), aaa gradient nonlinear MPC is described. and K¨ pernick (2012), gradient algorithm tailored tailored to point (x , u ) : f (x , u ) = 0 . (3) SP SP SP SP point nonlinear MPC is described. nonlinear MPC is (x = 0 .. (3) nonlinear MPC is described. described. SP ,, u SP ) SP ,, u SP ) (x u ) :: fff (x (x uSP = In a similar spirit as the last references, a fixed-point The MPC scheme in))this relies (3) on (xSP uSP (xSP = 00paper . (3) SP ,that SP )is: considered SP , uSP In aa similar similar spirit aspresented the last last references, references, fixed-point The MPCformulation scheme that that is considered considered in this this paper relies relies on In spirit fixed-point MPC iteration scheme wasas Graichenaaa (2012) with The the OCP In a similar spirit as the the last in references, fixed-point The MPC scheme scheme that is is considered in in this paper paper relies on on iteration scheme scheme was presented presented in Graichen Graichen (2012) with with the OCP formulation iteration iteration scheme was was presented in in Graichen (2012) (2012) with the the OCP OCP formulation formulation Copyright IFAC 2015 62 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2015, IFAC (International Federation of Automatic Control) Copyright IFAC 2015 62 Copyright © IFAC 2015 62 Peer review© of International Federation of Automatic Copyright ©under IFAC responsibility 2015 62 Control. 10.1016/j.ifacol.2015.11.060
62 62 62 62
IFAC CAO 2015 October 6-9, 2015. Garmisch-Partenkirchen, Knut Graichen et al. / IFAC-PapersOnLine 48-25 (2015) 062–067 Germany
min
J(xk , u ¯) = V (¯ x(T )) +
T
l(¯ x(τ ), u ¯(τ )) dτ
Hi
(4a)
0
x ¯˙ (τ ) = f (¯ x(τ ), u ¯(τ )) , x ¯(0) = xk = x(tk ) (4b) − + u ¯(τ ) ∈ [u , u ] , τ ∈ [0, T ] . (4c) with the MPC-internal variables x ¯, u ¯, the prediction time coordinate τ ∈ [0, T ], and the horizon length T . The initial value xk denotes the measured (or estimated) state at time step tk = t0 +k∆t, where 0 < ∆t ≤ T is the sampling time of the process. The integral and terminal cost functions in (4a) are usually designed as positive definite functions + n l : Rn × Rm → R+ 0 and V : R → R0 . Moreover, for the fixed-point iteration scheme in Section 3, the integral cost function is constructed in the separable form m 1 l(x, u) = l0 (x) + ri (ui − uSP,i )2 , ri > 0 . (5) 2 i=1 s.t.
63
control interval
u0i
u+ i
ui
Fig. 1. Illustration of minimizing control function (11). solution (6), Pontryagin’s Maximum Principle states that ¯ ∗ (τ ), τ ∈ there exists a corresponding adjoint trajectory λ k [0, T ] such that the canonical equations (9) hold and the ∗ optimal control u ¯k (τ ) satisfies ∗ ¯ ∗ (τ ), u) x∗k (τ ), λ (10) u ¯k (τ ) = arg minu∈[u− ,u+ ] H(¯ k for all times τ ∈ [0, T ]. In order to use the optimality conditions for computing the optimal solution (6), note that the particular structure of the integral cost (5) allows one to split (10) into the separate minimization problems ri Hi (x, λ, ui ) := (ui −uSP,i )2 +λT gi (x)ui min 2 ui ∈[u− ,u+ ] i i
In MPC, it is often assumed that the optimal solution of the OCP (4) at time step tk u ¯∗k (τ ) := u ¯∗ (τ ; xk ) , (6) ∗ ¯∗ (τ ; xk , u ¯∗k ) , τ ∈ [0, T ] x ¯k (τ ) := x ¯∗k ) is instantawith the optimal cost J ∗ (xk ) := J(xk , u neously available and is used as control update for the system (1) on the time interval [tk , tk+1 ), i.e. u(tk + τ ) = u ¯∗k (τ ) , tk = t0 + k∆t , τ ∈ [0, ∆t) . (7) The OCP (4) is then re-solved in the next time step tk+1 = tk + ∆t with the new system state xk+1 = x(tk+1 ).
for the single controls ui , i = 1, . . . , m. Since Hi is strictly convex in ui for ri > 0, the minimizing control function is given by 0 − + 0 ui if ui ∈ (ui , ui ) − − 0 ui = hi (x, λ) = ui if ui ≤ ui (11) + + 0 ui if ui ≥ ui with the unconstrained control 1 u0i = uSP,i − λT gi (x) , i = 1, . . . , m . ri The case-dependent definition of (11) is illustrated in Fig. 1. Inserting the optimal control function h1 (x, λ) .. u = h(x, λ) = (12) . hm (x, λ)
The endpoint-free formulation of the OCP (4) is chosen on purpose and will be of advantage for the fixedpoint iteration scheme in Section 3. Moreover, stability of an MPC scheme without terminal constraints can still be guaranteed under certain assumptions. In most cases, it is assumed that the terminal cost V (x) represents a control Lyapunov function (CLF), see e.g. Parisini and Zoppoli (1995); Chen and Allg¨ ower (1998). In the absence of the terminal cost V (x), the prediction time T must be sufficiently long in order to guarantee closed-loop stability (Jadbabaie and Hauser, 2005; Grimm et al., 2005; Gr¨ une and Pannek, 2011). 3. FIXED-POINT ITERATION SCHEME The particular structure of the OCP (4) with the integral cost (5) can be exploited by solving the respective optimality conditions in a fixed-point iterative manner. This section briefly summarizes the concept of the fixedpoint iteration scheme and derives the algorithm. Details as well as convergence and stability results can be found in Graichen (2012).
into (9) leads to the new form of the canonical equations x˙ = F (x, λ) , x(0) = xk (13a) (13b) λ˙ = G(x, λ) , λ(T ) = Vx (x(T )) with the new functions F (x, λ) := f (x, h(x, λ)) and G(x, λ) := −Hx (x, λ, h(x, λ)). The boundary value problem (13) is independent of the control u and consists of 2n ordinary differential equations with n initial conditions for the state x and n terminal conditions for the adjoint state λ.
3.1 Optimality conditions The fixed-point iteration scheme is based on the optimality conditions of the OCP (4). With the Hamiltonian H(x, λ, u) = l(x, u) + λT f (x, u) , the canonical equations follow as x˙ = f (x, u) , x(0) = xk ˙λ = −Hx (x, λ, u) , λ(T ) = Vx (x(T ))
u− i
3.2 Algorithm
(8)
The separated boundary conditions for x and λ can be taken advantage of by solving (13a) and (13b) successively forward and backward in time. Table 1 summarizes the algorithm, which basically consists of two numerical integrations. The additional damping steps involve the damping factor ε ∈ (0, 1] and the intermediate trajectories
(9a)
(9b) with the adjoint state λ ∈ Rn and the partial derivative functions Hx and Vx . If the OCP (4) has an optimal 63
IFAC CAO 2015 64 Knut Graichen et al. / IFAC-PapersOnLine 48-25 (2015) 062–067 October 6-9, 2015. Garmisch-Partenkirchen, Germany
Table 1. Damped fixed-point iteration scheme.
iterations, N , in each MPC step k. This leads to the (N ) ¯ (N ) (τ ), τ ∈ [0, T ] as well suboptimal trajectories x ¯k (τ ), λ k as the corresponding suboptimal control trajectory (12) (N ) (N ) ¯ (N ) (τ ) , τ ∈ [0, T ] . u ¯k (τ ) = h x ¯k (τ ), λ (16) k
1) Initialization • Damping factor ε ∈ (0, 1] (0) • Initial guess x ¯k (τ ), τ ∈ [0, T ], e.g. using (14) (0) ¯ • Initial guess λk (τ ), τ ∈ [0, T ], e.g. using (15)
Thus, the optimal feedback in (7) for the system (1) is replaced by the suboptimal one
2) Damped fixed-point iteration (1 ≤ j ≤ N ) • Forward integration of (j) x ˆ˙ k (τ ) = F
(j)
(j−1)
¯ x ˆk (τ ), λ k
• Damping of state trajectory (j)
(j)
(j)
(τ ) , x ˆk (0) = xk (j−1)
xk (τ ) + (1 − ε)¯ xk x ¯k (τ ) = εˆ
• Backward integration of
(N )
¯k (τ ) , τ ∈ [0, ∆t) . (17) u(tk + τ ) = u In the next MPC step k + 1, the damped fixed-point iteration scheme can be initialized with the previous trajectories (0) (N ) ¯ (0) (τ ) = λ ¯ (N ) (τ ), τ ∈ [0, T ] . (18) (τ ) = x ¯ (τ ), λ x ¯
(τ ) ,
τ ∈ [0, T ]
k+1
(j) (j) ˆ˙ (j) (τ ) = G x ˆ (j) (τ ) , λ ˆ (j) (T ) = Vx (¯ λ ¯k (τ ), λ xk (T )) k k k
k+1
k
The stability properties of the suboptimal MPC scheme with a constant number N of fixed-point iterations are investigated in Graichen (2012) for the undamped case (ε = 1). By defining the optimization error (N ) (N ) ∆¯ ¯k (·) − x ¯∗k (·) Ln [0,T ] (19) xk (·) Ln [0,T ] := x
• Damping of adjoint state trajectory
ˆ (j) (τ ) + (1 − ε)λ ¯ (j−1) (τ ) , ¯ (j) (τ ) = ελ λ k k k
k
Hence, the integration of (14) and (15) as stated in Table 1 is only required in the first MPC step k = 0.
τ ∈ [0, T ]
• Stop if suitable convergence criterion is fulfilled or if j = N . Otherwise, set j ← j +1 and return to 2).
∞
(j) ˆ (j) (τ ), τ ∈ [0, T ]. The fixed-point iterations x ˆk (τ ) and λ k are undamped for ε = 1, whereas strong damping is achieved for ε → 0. It is shown in Graichen (2012) that for a given set of initial values (xk ∈ X ⊂ Rn ) an upper horizon time Tmax > 0 exists such that for all T < Tmax the undamped fixed-point algorithm (ε = 1) has a linear rate of convergence, i.e. (j) (j−1) ¯ ∗ (·) n ¯ ∗ (·) n ¯ λk (·) − λ ≤ p ¯ λk (·) − λ k k L∞ [0,T ] L [0,T ] (j−1) ∞ (j) ∗ ∗ x ¯k (·) Ln [0,T ] ≤ p x ¯k (·) − x ¯k (·) Ln [0,T ] ¯k (·) − x ∞
(i) (ii) (iii)
∞
The fixed-point iterative MPC scheme is numerically and experimentally validated in order to show the computational efficiency as well as the practical applicability. 4.1 Simulation example: chemical reactor The following nonlinear model describes the dynamics of a continuous stirred tank reactor (CSTR) c˙A = −k1 (ϑ)cA − k2 (ϑ)c2A + [cin − cA ]u1
with the (constant) control input uSP of the setpoint (3). The initial trajectory of the adjoint state λ then follows from the backward integration of the adjoint system (13b) (0) (0) ¯˙ (0)(τ ) = G x ¯ (0)(τ ) , λ ¯ (0)(T ) = Vx (¯ λ ¯ (τ ), λ xk (T )) . (15) k k k k
(20a)
(20b) c˙B = k1 (ϑ)cA − k1 (ϑ)cB − cB u1 ˙ ϑ = −δ k1 (ϑ)cA ∆HAB + k1 (ϑ)cB ∆HBC (20c) 2 +k2 (ϑ)cA ∆HAD + α[ϑc − ϑ] + [ϑin − ϑ]u1 ˙ (20d) ϑc = β[ϑ − ϑc ] + γu2
3.3 MPC implementation
k
∞
4. NUMERICAL AND EXPERIMENTAL VALIDATION
The start of the fixed-point iteration scheme requires an (0) ¯ (0) (τ ), τ ∈ initial guess of the trajectories x ¯k (τ ) and λ k (0) [0, T ]. One possibility is to compute x ¯k (τ ) by forward integration of the system dynamics (1) (0) (0) (0) x ¯˙ k (τ ) = f x ¯k (0) = xk ¯k (τ ), uSP , x (14)
k
the fixed-point iterations are convergent (T < Tmax ), N satisfies a lower bound, (0) ∆¯ x0 (·) Ln [0,T ] satisfies an upper bound,
it is shown in Graichen (2012) that asymptotic stability as well as incremental reduction of the optimization error (19) is ensured.
for some p ∈ (0, 1). In practice, the maximum horizon length Tmax for which convergence of the algorithm is ensured can be enlarged by introducing damping as it is done in Table 1. The influence of the damping factor ε ∈ (0, 1] on the convergence is illustrated numerically for a simulation example in Section 4.1.
Typically, the algorithm in Table 1 would iterate until sufficient accuracy is achieved by choosing a suitable convergence criterion. 1 However, in order to implement the algorithm within a real-time MPC framework, it is more convenient to use a constant number of fixed-point (j) 1 For instance, λ ¯ (j−1) (·) ¯ (·) − λ ≤ δλ for some δλ > 0. n
∞
as a measure for the suboptimality in MPC step k and assuming that
with the Arrhenius functions −Ei ki (T ) = ki0 exp , i = 1, 2 . T /◦ C + 273.15 Details on the model in (20) as well as the model parameters can e.g. be found in Klatt and Engell (1998); Graichen et al. (2009). The state vector x = [cA , cB , ϑ, ϑc ]T comprises the educt and product concentration as well as the reactor and coolant temperature, respectively. The control
L∞ [0,T ]
64
IFAC CAO 2015 October 6-9, 2015. Garmisch-Partenkirchen, Knut Graichen et al. / IFAC-PapersOnLine 48-25 (2015) 062–067 Germany
vector u = [u1 , u2 ]T consists of the normalized input flow rate and the cooling power, subject to the constraints MJ (21) u1 ∈ [5 h−1 , 35 h−1 ] , u2 ∈ [−8.5 MJ h ,0 h ]. The control goal that is used to evaluate the fixed-point iterative MPC scheme is to maximize the product concentration cB while maintaining a desired setpoint ϑSP for the reactor temperature. To this end, the cost functional in (4a) is chosen as 2 T 2 ¯ J(xk , u ¯) = −qB c¯2B +qϑ (ϑ−ϑ ¯21 +r2 u ¯22 dτ (22) SP ) +r1 u
65
16 14
Tmax [min]
12 10 8 6 4
0
6
m 2◦ 2 C h, with the coefficients qB = 102 kmol 2 h , qϑ = 10 −1 0 h r1 = 10 h, and r2 = 10 MJ2 . The sampling time for the MPC scheme is set to ∆t = 1 s. The initial state of the reactor is the stationary setpoint (Utz et al., 2007) T ◦ ◦ kmol x0 = 3.5175 kmol . (23) m3 , 0.74 m3 , 87 C , 79.8 C The damped fixed-point iteration scheme in Table 1 is implemented as C code with a Cmex interface for Matlab. The canonical equations (13) including the optimal control function (12) are computed under Mathematica and are exported as optimized Cmex functions. The forward and backward integrations in the algorithm are performed using a Heun integration scheme on an uniform discretization of the MPC interval [0, T ] with Ndiscr points.
2 0 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
ε ∈ (0, 1]
Fig. 2. Maximum horizon length Tmax for convergence of the fixed-point iterations w.r.t. the damping factor ε ∈ (0, 1]. CPU (M620, 2.67 GHz) and Matlab R2012a (Windows 7, 64 bit). The damping factor and horizon length are set to ε = 0.1 and T = 10 min according to the simulation results in Fig. 3. The last column in Table 2 shows the overall cost value Tsim Jsim = ¯22 dτ (24) −qB c¯2B + qϑ (ϑ¯ − ϑdes )2 + r1 u ¯21 + r2 u
In the first step, it is investigated how the damping factor ε ∈ (0, 1] of the fixed-point iteration scheme affects the maximum horizon length Tmax for which convergence is achieved. To this end, T is successively increased for a constant ε until the MPC scheme in closed-loop turns unstable. The resulting ε-T -curve is depicted in Fig. 2. In the undamped case (ε = 1), the maximum horizon length is Tmax = 0.7 min, whereas for increased damping of the fixed-point iterations Tmax can be significantly enlarged.
0
that is obtained by integration of the simulated trajectories over the simulation time Tsim = 70 min. The value of Jsim can be seen as an indicator for the control performance that decreases for larger values of Ndiscr and N . The statistics in Table 2 as well as the trajectories in Fig. 3 show that an excellent performance is already obtained for N = 2 fixed-point iterations per MPC step with a computation time of 32 µs.
Figure 3 shows the simulation results for the CSTR with N = 2 fixed-point iterations per MPC step and the horizon length T = 10 min, which is sufficiently short for the chosen damping factor ε = 0.1 (cf. Fig 2). In the first time interval [0, 10) min, the MPC is inactive and the reactor stays at the initial setpoint (23). In the second time interval [10, 40) min, the MPC is activated with the desired temperature ϑSP = 87◦ C corresponding to the initial setpoint x0 . By minimizing the cost functional (22), the MPC increases the product concentration cB while maintaining the reactor temperature ϑ in the vicinity of the desired value ϑSP = 87◦ C. During the third time interval [40, 70) min, the desired reactor temperature is increased to the setpoint ϑSP = 110◦ C which also results in a larger product concentration. During both transitions, the controls act in an aggressive manner owing to the comparatively small control weights in the cost functional (22).
Table 2. MPC results for the CSTR. Discretization points Ndiscr
Fixed-point iterations N
CPU time [µs]
Simulated cost Jsim [-]
15
2 5 10
32 55 96
17.56 17.47 17.44
30
2 5 10
58 98 183
15.93 15.87 15.86
60
2 5 10
83 186 355
15.78 15.55 15.54
4.2 Experiment: magnetic levitation system Figure 4 shows the schematic drawing of a magnetic levitation system with an electromagnet and a levitating ferromagnetic object with mass m = 133 g. The electromagnet consists of two coils (2000 windings each) on a steel cup core to generate the magnetic field. The electromagnetic force is applied to the levitating mass in the opposite direction of the gravitational force. The control input is the duty cycle τ ∈ [0, 1] of the pulse width modulated (PWM) voltages applied to each coil. The measured variables are the currents i1 and i2 of both coils, the temperature inside the coils, and the levitation height z by means of a laser
In order to evaluate the computational efficiency in addition to the simulation results in Fig. 3, the runtime requirements are recorded for different numbers of discretization points Ndiscr and fixed-point iterations N per MPC step. Table 2 summarizes the runtime results for the Cmex compiled MPC scheme on a computer with Intel Core i7 2
Note that the goal to maximize the product concentration leads to an integral cost function that is not positive definite. This setting somewhat differs from the classical stabilization task and belongs to the framework of economic MPC, see e.g. (Ferramosca et al., 2010).
65
IFAC CAO 2015 66 Knut Graichen et al. / IFAC-PapersOnLine 48-25 (2015) 062–067 October 6-9, 2015. Garmisch-Partenkirchen, Germany
110
3 2.5 2
1.2
105
ϑ [◦ C]
cB [kmol/ m3 ]
cA [kmol/ m3 ]
4 3.5
1 0.8
95 90
0.6
1.5 10
20
30
40
50
60
70
10
20
time [min]
30
40
50
60
85
70
95 90
25 20 15
85
10
80
5 30
40
50
60
u 2 [MJ/h]
u 1 [1/ h]
100
70
30
40
50
60
70
50
60
70
0
30
20
20
time [min]
35
105
10
10
time [min]
110
ϑ c [◦ C]
100
−2 −4 −6 −8
0
10
20
time [s]
30
40
50
60
70
0
10
20
time [min]
30
40
time [min]
Fig. 3. MPC simulation results for the CSTR. for different values of zSP . To this end, the cost functional in (4a) is chosen as T J(xk , u ¯) = qz (¯ z − zSP )2 + qz˙ z¯˙ 2 + r (¯ τ − τSP )2 dτ (25)
sensor. The experimental setup is constructed for large levitation heights of up to 7 cm which necessitates an active air cooling inside the coils. Moreover, no cascaded current controllers are used in order to achieve the fastest possible control performance.
0
with the weights qz = 104 m12 s , qz˙ = 25 ms2 , and r = 0.25 1s . The stationary control input (duty cycle) is set to τSP = 0.5 for the sake of simplicity. An additional integral gain is added to the MPC scheme in order to compensate for steady state errors in the position of the levitating mass. 1 The MPC sampling time ∆t = 700 Hz ≈ 1.4 ms corresponds to the PWM frequency. The prediction horizon and the damping factor for the fixed-point iterations are set to T = 70 ms and ε = 0.1, respectively.
In B¨ achle et al. (2013), a nonlinear model of the magnetic levitation systen is derived and an alternative gradientbased MPC scheme is used to control the system. In the following, the nonlinear model of this levitation system with the state x = [i1 , i2 , z, z] ˙ T and control τ ∈ [0, 1] is used to demonstrate the practical applicability of the fixed-point iteration scheme. The control task is to stabilize a desired setpoint of the levitation height zSP and to perform fast setpoint changes
The MPC scheme uses the Matlab/C implementation of the fixed-point iteration scheme described in Section 4.1. Table 3 shows the runtime results for the MPC implementation using Ndiscr = 30 discretization points for the prediction horizon [0, T ]. For the experimental implementation, a dSPACE MicroAutoBox (800 MHz) is used as real-time hardware. The computation time on the dSPACE box amounts to approximately 400 µs including the MPC scheme and state estimation (Unscented Kalman Filter), which is well below the sampling time.
air gaps for active cooling coil 1
Figure 5 shows the measurement results for different setpoint changes of the levitation height z between −70 mm and −45 mm. Overall, the MPC scheme shows a very good control performance. The noise in the control inputs is mainly due to the noisy measurements and filtered signals that leads to persistent perturbations of the MPC scheme through the initial conditions in the OCP (4). This effect also demonstrates the robustness of the fixed-point iteration scheme. In addition, the last plot of Fig. 5 depicts the CPU requirements per MPC step plotted over time
coil 2
z cup core levitating mass
Table 3. MPC results for the levitation system.
Fig. 4. Schematic of the magnetic levitation system. 66
Discretization points Ndiscr
Fixed-point iterations N
CPU time [µs]
dSPACE time [µs]
30
2
60
378
IFAC CAO 2015 October 6-9, 2015. Garmisch-Partenkirchen, Knut Graichen et al. / IFAC-PapersOnLine 48-25 (2015) 062–067 Germany
67
z [cm]
REFERENCES −5
B¨ achle, T., Hentzelt, S., and Graichen, K. (2013). Nonlinear model predictive control of a magnetic levitation system. Contr. Eng. Pract., 21(9), 1250–1258. Bemporad, A., Morari, M., Dua, V., and Pistikopoulos, E. (2002). The explicit linear quadratic regulator for constrained systems. Automatica, 38(1), 3–20. Cannon, M., Liao, W., and Kouvaritakis, B. (2008). Efficient MPC optimization using Pontryagin’s Minimum Principle. Int. J. Robust Nonl. Contr., 18, 831–844. Chen, H. and Allg¨ ower, F. (1998). A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability. Automatica, 34(10), 1205–1217. Diehl, M., Bock, H., and Schl¨ oder, J. (2005). A real-time iteration scheme for nonlinear optimization in optimal feedback control. SIAM J. Contr. Optim., 43(5), 1714–1736. Ferramosca, A., Rawlings, J., and Camacho, D.L.E. (2010). Economic MPC for a changing economic criterion. In Proc. 49th IEEE Conf. Dec. Contr. (CDC), 6131–6136. Atlanta (USA). Ferreau, H., Bock, H., and Diehl, M. (2008). An online active set strategy to overcome the limitations of explicit MPC. Int. J. Robust Nonl. Contr., 18(8), 816–830. Graichen, K. (2012). A fixed-point iteration scheme for real-time model predictive control. Automatica, 48(7), 1300–1305. Graichen, K., Hagenmeyer, V., and Zeitz, M. (2009). Design of adaptive feedforward control under input constraints for a benchmark CSTR based on a BVP solver. Comput. Chem. Eng., 33(2), 473–483. Graichen, K. and K¨ apernick, B. (2012). A real-time gradient method for nonlinear model predictive control. In Frontiers of Model Predictive Control, 9–28. InTech (open access). http://www. intechopen.com/articles/show/title/a-real-time-gradientmethod-for-nonlinear-model-predictive-control. Graichen, K. and Kugi, A. (2010). Stability and incremental improvement of suboptimal MPC without terminal constraints. IEEE Trans. Automat. Contr., 55(11), 2576–2580. Grimm, G., Messina, M., Tuna, S., and Teel, A. (2005). Model predictive control: For want of a local control Lyapunov function, all is not lost. IEEE Trans. Automat. Contr., 50(5), 546–558. Gr¨ une, L. and Pannek, J. (2010). Analysis of unconstrained NMPC schemes with incomplete optimization. In Proc. 8th IFAC Symp. Nonl. Contr. Syst. (NOLCOS), 238–243. Bologna (Italy). Gr¨ une, L. and Pannek, J. (2011). Nonlinear Model Predictive Control: Theory and Algorithms. Springer, London. Houska, B., Ferreau, H., and Diehl, M. (2011). An auto-generated real-time iteration algorithm for nonlinear MPC in the microsecond range. Automatica, 47(10), 2279–2285. Jadbabaie, A. and Hauser, J. (2005). On the stability of receding horizon control with a general terminal cost. IEEE Trans. Automat. Contr., 50(5), 674–678. Klatt, K.U. and Engell, S. (1998). Gain-scheduling trajectory control of a continuous stirred tank reactor. Contr. Eng. Pract., 22(4-5), 491–502. Ohtsuka, T. (2004). A continuation/GMRES method for fast computation of nonlinear receding horizon control. Automatica, 40(4), 563–574. Parisini, T. and Zoppoli, R. (1995). A receding–horizon regulator for nonlinear systems and a neural approximation. Automatica, 31(10), 1443–1451. Utz, T., Hagenmeyer, V., and Mahn, B. (2007). Comparative evaluation of nonlinear model predictive and flatness-based twodegree-of-freedom control design in view of industrial application. J. Process Contr., 17(2), 129–141. Zeilinger, M., Jones, C., and Morari, M. (2011). Real-time suboptimal model predictive control using a combination of explicit MPC and online optimization. IEEE Trans. Automat. Contr., 56(7), 1524–1534.
−6 −7 0
2
4
6 time [s]
8
12
coil 1 coil 2
6 i 1 & i 2 [A]
10
4 2 0
2
4
6 time [s]
8
10
12
0
2
4
6 time [s]
8
10
12
0
2
4
6 time [s]
8
10
12
τ [-]
1 0.5
CPU time [ms]
0
1
0.5
0
Fig. 5. MPC measurement results for the magnetic levitation system. (also cf. Table 3). The computation time is practically constant, which illustrates the real-time feasibility of the MPC scheme based on a constant number of fixed-point iterations per sampling step ∆t. A video of the MPC-controlled magnetic levitation system is available via the link http://youtu.be/01q4TUVrToc. 5. CONCLUSIONS The numerical and experimental results in the paper demonstrate the low complexity and computational cost of the fixed-point iteration scheme recently presented for MPC of nonlinear systems subject to control constraints. In general, the horizon length cannot be chosen arbitrarily large in order to ensure convergence as well as stability in closed-loop. This can be interpreted as a payoff for the low algorithmic complexity of the fixed-point iteration scheme. On the other hand, the introduction of damping in the algorithm usually allows one to significantly enlarge the horizon length, as demonstrated for the chemical reactor example. A key role for tuning the fixed-point iterative MPC scheme plays the choice of the damping factor in relation to the prediction horizon and the number of fixedpoint iterations per MPC step. Numerical studies have shown that two iterations per MPC step are often sufficient to obtain a good control performance. 67