Estimates of the Prediction Horizon Length in MPC: a Numerical Case Study

: In this paper we are concerned with estimates of the prediction horizon length in nonlinear model predictive control (MPC) without terminal constraints or costs for systems governed by ordinary diﬀerential equations. A growth bound — which is known to be the crucial condition in order to determine a horizon length for which asymptotic stability or a desired performance of the MPC closed loop is guaranteed — is numerically deduced for an example of a synchronous generator. Then, the system dynamics are discretized and the computations are repeated for the resulting sampled data system. We investigate how the obtained estimates are related — in particular, for sampling periods tending to zero. Furthermore, it is shown that a suitable design of the running costs in the sampled data setting can lead to improved performance bounds and, thus, can ensure stability for signiﬁcantly shorter prediction horizons.


INTRODUCTION
Model predictive control (MPC) is a control strategy in order to approximately solve an optimal control problem on an infinite time horizon.To this end, a sequence of optimal control problems on a truncated and, thus, finite time horizon is iteratively solved.MPC is particularly attractive due to its ability to take hard constraints directly into account, cf.Maciejowski (2002).However, stability of the original problem may get lost if the prediction horizon is chosen too short, see, e.g., Raff et al. (2006).Such drawbacks can be excluded if (artificial) terminal constraints and/or costs are incorporated as shown in Keerthi and Gilbert (1988); Chen and Allgöwer (1998).In this setting the stability analysis is quite mature, even for nonlinear and infinite dimensional systems, cf.Camacho and Bordons (1999); Rawlings and Mayne (2009); Ito and Kunisch (2002).However, the construction of suitable stabilizing terminal constraints and/or costs remains a challenging task which explains why so called unconstrained MPC is often used in practice, cf.Qin and Badgwell (2003).Here, the term unconstrained emphasizes that no additional terminal constraints and/or costs are used, control and state constraints can, however, be taken into account.
For unconstrained MPC schemes stability can be concluded for sufficiently long prediction horizons, cf.Jadbabaie and Hauser (2005).The first estimates of the required horizon length in a nonlinear setting can be found in Grimm et al. (2005); Tuna et al. (2006).A further improvement was obtained by using the methodology introduced in Grüne (2009); Grüne et al. (2010), cf. Worthmann (2012) for a comparison.All these techniques have in common that a growth condition is assumed in which optimal open loop costs are set in relation to a "reference quantity" depending on the corresponding initial state.In these references a discrete time framework is considered, the approach, however, can also be applied in a continuous time setting, cf.Reble and Allgöwer (2011).The connection between the resulting estimates of the prediction horizon is investigated in Worthmann et al. (2012) for systems satisfying an exponential controllability condition.
In this paper this investigation is continued, however, without supposing exponential controllability.Instead a growth condition is numerically verified.Then, different approximations of the continuous time cost functional are considered.Two main conclusions can be drawn: firstly, convergence of the horizon estimates seems to hold for sampling periods tending to zero which substantiates the conjecture that the results from Worthmann et al. (2012) can be transferred to our setting based on the weaker growth condition.Secondly, for a given sampling period, the use of a suitably chosen discrete time approximation may significantly enhance the corresponding performance estimates -even in comparison to the continuous time setting.This allows to ensure asymptotic stability for much shorter prediction horizons although the MPC cost functional and, thus, the (numerical) effort in order to verify the required growth condition remain unchanged.The reason for this phenomenon is that both the current state and the influence of the control are reflected by the reference quantity in the growth condition if the discrete time running costs are suitably designed.
The paper is organized as follows.In the ensuing Section 2 an MPC scheme and a methodology in order to compute a prediction horizon, for which asymptotic stability of the MPC closed loop is guaranteed, are presented.Then, in Section 3 the results are applied for a nonlinear example.In the following Section 4 the technique introduced in Section 2 is transferred to sampled data systems.To this end, the concept of multistep feedback laws is required.In Section 5 the main results dealing with the relation of the continuous setting and its sampled data implementation are presented.Finally, conclusions are drawn in Section 6.

PROBLEM FORMULATION & RECAP
A continuous function η : R ≥0 → R ≥0 is said to be of class K ∞ if it is strictly increasing, unbounded, and zero at zero.Furthermore, a continuous function β : R ≥0 × R ≥0 → R ≥0 is called a KL-function if β(•, t) is of class K ∞ for each t ≥ 0 and β(r, t) → 0 for t → ∞ holds for all r > 0.
We consider systems governed by nonlinear autonomous (ordinary) differential equations ẋ(t) = f (x(t), u(t)). (1) Here, x(t) ∈ R n and u(t) ∈ R m represent the state and the control at time t, respectively.State and control constraints are modeled by suitably chosen subsets X ⊆ R n and U ⊆ R m .For a given state x and time T a control function u : [0, T ) → R m is said to be admissible -which is denoted by u ∈ U(x, T ) -if the corresponding solution x(t; x, u) exists for each t ∈ [0, T ] and the conditions are satisfied.Furthermore, u : R ≥0 → R m is called admissible on the infinite horizon, written u ∈ U(x, ∞), if the restriction u| [0,T ) is contained in U(x, T ) for all T > 0.

Problem Formulation & Model Predictive Control
Let x ∈ X be a state for which a control input u ∈ U exists such that f (x , u ) = 0 holds.(x , u ) is said to be an equilibrium pair.Furthermore, let running costs : R n ×R m → R ≥0 with (x , u ) = 0 be defined.Then, for given initial state x 0 ∈ X and admissible control function u ∈ U(x 0 , ∞), the performance on the infinite time horizon can be evaluated by using the cost functional (3) Our goal is to find a feedback map µ : X → U such that, for each feasible initial state x 0 ∈ X, the resulting closed loop solution x µ (•; x 0 ) generated by ẋµ (t; x 0 ) = f (x µ (t; x 0 ), µ(x µ (t; x 0 ))), x µ (0; x 0 ) = x 0 , exists, satisfies the constraints (2), minimizes the costs (3), and asymptotically converges to x , i.e. there exists β ∈ KL such that x µ (t; x) − x ≤ β( x − x , t) holds for all t ≥ 0. In order to ensure that this task can be satisfied, the following assumption is necessary.
Since infinite horizon optimal control problems are, in general, hard to solve, we use model predictive control (MPC) as a remedy.For a given initial condition x := x 0 , an MPC algorithm typically consists of three steps: (1) Solve the optimal control problem on a truncated and, thus, finite time horizon, i.e. compute a control function û which minimizes over all admissible control functions, i.e. u ∈ U(x, T ).
Here, in order to keep the presentation technically simple, existence of a minimizer is assumed.(2) For given δ ≤ T , define a feedback law µ of the computed control function û is implemented at the plant by using the MPC feedback law µ T .(3) Then, the current state is updated by setting x := x µ T (δ; x) = x(δ; x; û) and the prediction horizon is shifted forward in time.
Iterative application of this procedure generates a solution x µ T (•; x 0 ) on the infinite time horizon.The corresponding control function is denoted by µ M P C (•; x 0 ).

Estimates of the Prediction Horizon & Stability
In Reble and Allgöwer (2011) a methodology was introduced which allows to determine a prediction horizon length T for which asymptotic stability or a desired performance of the MPC closed loop is guaranteed.To this end, the notation of the optimal value function V for all x ∈ X be given.In addition, let T > δ be chosen such that α T > 0 holds with Then, the relaxed Lyapunov inequality holds for all x ∈ X and, as a consequence, the performance estimate hold for x ∈ X, the MPC closed loop is asymptotically stable with prediction horizon T .
Proof.Theorem 2 mainly summarizes results from Reble and Allgöwer (2011) with two modifications.Firstly, continuity of the function B(•) is not assumed.Secondly, state constraints are included in our setting.Both changes are based on the observation that the respective proofs presented in Reble and Allgöwer (2011) remain valid, see also the corresponding results derived for a discrete time system in Grüne (2009) and Grüne et al. (2010).
The main assumption needed in order to apply Theorem 2 and, thus, to conclude asymptotic stability or a performance guarantee is the growth condition (5).If this inequality is satisfied, α T > 0 always holds for a sufficiently large prediction horizon.Corollary 3. Let B : R ≥0 → R ≥0 be a monotone and bounded function with B(t) > 0 for t > 0.Then, the performance bound α T given by ( 6) converges to unity for T approaching infinity, i.e. lim T →∞ α T = 1.
Proof.Let an arbitrary but fixed δ > 0 be given.Since B(•) is monotone, bounded and satisfies In addition, the inequality e can be concluded.Hence, Corollary 3 can be shown by combining these two assertions.
Corollary 3 generalizes (Reble and Allgöwer, 2011, Section 4) to systems which are not exponentially controllable in terms of the stage costs.Remark 4. We point out that Condition ( 5) is only needed on the interval [0, T ] in order to prove Theorem 2 -a fact which facilitates verifying this assumption numerically.

SYNCHRONOUS GENERATOR
We consider the example of the synchronous generator given by with parameters b 1 = 34.29,b 2 = 0.0, b 3 = 0.149, b 4 = 0.3341, P = 28.22,and E = 0.2405, cf.Galaz et al. (2003).We want to stabilize this system at the equilibrium pair (x , u ) with x ≈ (1.124603730, 0, 0.9122974248) T and u = 0.The running costs are used where • denotes the Euclidean norm.The parameter λ = 0.01 penalizes the control effort used to manipulate the system behavior.Furthermore, the physically motivated state constraints 0 ≤ x 1 < π/2 and x 3 ≥ 0 have to be taken into account.
Our goal is to numerically verify the introduced growth condition (5) for the sub-level set X := {x ∈ R 3 : V 0.6 (x) ≤ 0.09 2 } which is control invariant according to our numerical experiments, see, e.g., Blanchini and Miani (2008) for a definition of control invariance.To this end, we proceed in two steps: (1) For each state x ∈ X a function B x : [0, T ] → R ≥0 is determined such that (5) holds.Here, the subscript x indicates the dependence on the state x.
(2) Then, the desired function B(•) is obtained as pointwise supremum, i.e.B(t The first step is only carried out approximately.The sublevel set X is contained in the interior of the cube with a 1 := 0.4, a 2 := 0.5, and a 3 := 0.9.This cube is discretized with stepsize ∆x i = 0.02, i ∈ {1, 2, 3}, in each coordinate direction which results in a grid G.Then, the set X is defined as the intersection X ∩G, i.e. a sub-level set of the constructed grid consisting of 8309 grid points, see Fig. 1 for an illustration.Hence, the following derivation is rigorous apart from verifying (5) on the set X instead of X and, thus, only for a finite number of points.For a given state x ∈ X\{x } and time t, the condition has to be satisfied.Hence, an upper bound of V t (x) = inf u∈U (x,t) J t (x, u) has to be computed.To this end, we make use of the fact that optimality of B(t) and, thus, B x (t) is not needed.In principal, arbitrary L ∞ ([0, t), R)functions are contained in the set of admissible control functions U(x, t).However, we restrict ourselves to piecewise constant control functions.To be more precise, the time domain is discretized with stepsize τ = 0.0125 and control functions u : [0, T ) → R are confined to be constant on each interval [nτ, (n + 1)τ ), n = 0, 1, 2, . . ., T /τ − 1.On the one hand doing so leads to slightly more conservative bounds B x (t), t ∈ τ N, since the set of admissible control functions becomes smaller.On the other hand, for each x ∈ X and t ∈ τ N, the value V t (x) and, thus, B x (t) can be estimated by solving a nonlinear optimization problem.1 In addition, since neither terminal constraints nor terminal costs were used, V t (x) ≤ V (n+1)τ (x) and, as a consequence, B x (t) ≤ B x ((n + 1)τ ) hold for each t ∈ (nτ, (n + 1)τ ].
Then, a bound B(t) is given by sup x∈ X B x (t) for each t ∈ [0, T ].In conclusion, this approach yields a function B(t) satisfying (5) for each x ∈ X ⊂ X, cf.Fig. 2.
Next, Theorem 2 is applied in order to determine a prediction horizon for which asymptotic stability is guaranteed for the MPC closed loop with δ = 0.1.Here, we want to determine T as a multiple of δ.Then, T = 2.1 ensures the relaxed Lyapunov inequality ( 7) with α = 0.048518.Note that the coarser discretizations depicted in Fig. 2 require a prolongation to T = 2.2 and T = 2.3, respectively, which shows that improving the bound B pays off.

SAMPLED DATA SYSTEMS
Applying Theorem 2 requires that the input signal can be switched at each time instant t ∈ [0, T ) -an assumption which typically cannot be satisfied in practice.Instead, only a limited number of such switching instants is possible, a restriction which fits well to the derivation of the upper bounds with piecewise constant control functions as carried out in the previous Section 3. Here, we consider a sampled data system with zero order hold in order to take this into account, i.e. the successor state is given by the discrete time dynamics f : R n × R m → R n defined as x(0) = x 0 , with ũ(t) = u(n) for all t ∈ [0, τ ).Here, τ > 0 denotes the discretization parameter (sampling period) and Φ(•; x(n), ũ) denotes the solution of the differential equation ( 1) with constant control function ũ emanating from initial state x(n).Again, the state and control are denoted by x and u with a slight abuse of notation.
A sequence u = (u(0), u(1), . . ., u(N − 1)) of N control values is called admissible for x and N , denoted by U N (x), if the state trajectory x u (n; x), n = 1, 2, . . ., N , governed by system dynamics (12) exists and the conditions f (x u (n; x), u(n)) ∈ X and u(n) ∈ U (13) hold for each n ∈ {0, 1, 2, . . ., N − 1}.Furthermore, U ∞ (x) is defined analogously to the continuous time setting.In addition, Assumption 1 is still required for the optimal value function for each x ∈ X with running costs ˜ : R n × R m → R ≥0 , ˜ (x , u ) = 0.A necessary condition for this assumption is that the set X is control invariant, see, e.g., Kerrigan and Maciejowski (2000); Primbs and Nevistić (2000) for methods to ensure control invariance in MPC.
The MPC algorithm described in Section 2 has to be modified in step (2).Since now a sequence of control values is computed in step (1), the first p = δ/τ ∈ N elements of this sequence are implemented.To this end, the definition of a multistep feedback law µ N,p : {0, 1, 2, . . ., p − 1} × X → U is required, cf.Grüne (2009).For the considered MPC algorithm µ N,p (n, x) is defined as the (n + 1)-st element of the control sequence u ∈ U N (x) satisfying Next, we give conditions under which asymptotic stability or a desired performance bound can be guaranteed, cf.Grüne et al. (2010) and Worthmann (2012) for a proof.Theorem 5. Let p ∈ N be given and suppose that a monotone bounded sequence (M i ) i∈N ⊂ R ≥1 exists satisfying for each x ∈ X.In addition, let N > p be chosen such that α N,p > 0 holds with α N,p given by .
Then, the relaxed Lyapunov inequality and, as a consequence, a performance estimate analogously to (8) with α N,p instead of α T hold for all x 0 ∈ X and the MPC feedback law with prediction horizon N .If, additionally, K ∞ -functions , exist such that the inequalities V 1 (x) ≥ (x) and V N (x) ≤ (x) hold for all x ∈ X, the MPC closed loop is asymptotically stable.
The origin of Theorems 5 and 2 is the construction of a linear program which characterizes optimal solutions by Bellman's principle of optimality and a controllability condition instead of the growth conditions ( 16) and ( 5), cf.Grüne (2009).If the involved coefficients in this controllability condition exhibit a "submultiplicativity condition", the expressions α N,p and α T coincide with the corresponding solution.Otherwise solving the respective linear programs may improve the deduced bounds, cf.Grüne et al. (2010).We like to point out that the linear program has infinitely many constraints in the continuous time setting which explains why switching at each time instant is a required condition in order to ensure their satisfaction, cf.Reble and Allgöwer (2011) for further details.

CONNECTION BETWEEN THE CONTINUOUS AND THE DISCRETE TIME SETTING
The discrete time system dynamics ( 12) are a sampled data implementation with zero order hold of their contin-uous time companion piece (1).Furthermore, the integral in the optimal value function V T , T ∈ R >0 ∪ {∞}, is straightforwardly replaced by a sum of appropriate length.Intuitively, the running costs ˜ employed in the discrete time setting should then approximate the integral over on a sampling interval.Hence, one option is to use the running costs with ũ(t) = u(n), t ∈ [0, τ ) which reproduces the continuous time cost functional with the restriction to piecewise constant control functions ũ.However, the right hand side of Assumption ( 5) is based on (x) = x − x 2 for the running costs defined by ( 10) whereas its counterpart in (16), i.e.V 1 (x) := inf u∈U : f (x,u)∈X ˜ 1 (x, u), can be influenced by the control input u.
Another possible choice for the running costs is Since sampled data systems with zero order hold are considered, this definition leads to an approximation error with respect to the state trajectory.On the other hand, assuming admissibility of u = 0 the equality τ (x) = V 1 (x) holds.Hence, the reference quantity used in the right hand side of Assumptions ( 5) and ( 16) remains essentially the same.
For systems exponentially controllable in terms of their stage costs it has been shown that the continuous time estimates can be arbitrarily well approximated for sampling periods tending to zero, cf.Worthmann et al. (2012).
In addition, the corresponding performance index α T has proven to be an upper bound for the ones resulting from the discrete time setting.Do these assertions also hold without such an exponential controllability condition?
Moreover, we want to investigate whether the possibility to influence the reference quantity V 1 in ˜ 1 can be exploited in order to obtain better suboptimality estimates?

Numerical Results for ˜ 2
We begin our investigation with ˜ 2 and δ = 0.1.The control is allowed to change its value p ∈ {1, 2, 4, 8} times on the interval [0, δ), i.e. τ ∈ {0.1, 0.05, 0.025, 0.0125}.Hence, the concept of multistep feedback laws is required in order to take this into account.Then, for each point x ∈ X, a sequence M i is computed analogously to the continuous time setting considered in Section 3. Assuming f (x, 0) ∈ X the equality V 1 (x) = τ (x) holds.Hence, the reference quantity V 1 (x) does not depend on the control.In this setting, one observes that using a finer discretization leads to improved bounds (M i ) i∈N , as shown in Fig. 2 for the continuous time setting, and, thus, to improved suboptimality estimates.Indeed, the continuous time estimates are approximated by using smaller sampling periods in combination with multistep feedback laws, cf.Fig. 3.
We conjecture that carrying out a few additional refinement steps leads to estimates converging to the continuous time ones -analogously to the results rigorously derived for a setting based on systems exponentially controllable in terms of their stage costs, cf.Worthmann et al. (2012).This observation meets our expectations.

Design of suitable running costs
In this subsection we repeat our computations for the discrete time setting based on running costs ˜ 1 instead of ˜ 2 .Then, the optimal value functions from the left hand side of our growth condition ( 16) coincide with their counterparts from (5) apart from the restriction to piecewise control inputs.However, V 1 reflects the influence of the first control input which may lead to using a control input not equal to zero.As a consequence, the reference quantity may not be a scaled version of anymore.
The numerical computations lead to the estimates depicted in Fig. 4. Here, we observe a significant decrease in the required prediction horizon length needed in order to guarantee asymptotic stability or a desired performance specification in comparison to the continuous time setting for τ = δ = 0.1 -despite the restriction to piecewise constant control inputs which implies that changes are only permitted at multiples of the sampling rate τ .Moreover, τ = 0.1 yields the best estimates although these bounds result from the smallest class of control inputs.Hence, at least a slight improvement seems to be attainable by allowing for multirate sampling, cf. Lee et al. (1992).However, using smaller τ and, thus, larger p in order to keep δ constant worsens the performance bounds resulting from Theorem 5. Here, we point out that the reference quantity V 1 (•) depends on τ and, as a consequence, changes.Again, the suboptimality estimates obtained in the sampled data setting seem to converge to their continuous time counterparts for τ tending to zero.
The reference quantity V 1 (x) reflects both the state x as well as the dynamical behavior of the system which is determined by the interplay of the system dynamics and the control.Although the involved bounds (M i ) i∈N are worse due to the coarser discretization, the use of this accumulated information in the reference quantity reduces the required prediction horizon length in order to ensure asymptotic stability compared to their counterparts resulting from the continuous time setting.In conclusion, the discrete time approach allows to derive improved performance estimates by choosing the stage costs appropriately, e.g.˜ 1 instead of ˜ 2 for the considered example, and, thus, using additional information on the system behavior.

CONCLUSIONS AND OUTLOOK
In this paper a methodology in order to determine a prediction horizon length for which asymptotic stability or a desired degree of suboptimality is guaranteed for the MPC closed loop is numerically investigated for the example the synchronous generator.Firstly, we showed that the results obtained in a sampled data setting converge to their counterparts derived in a purely continuous time setting for sampling periods tending to zero.Secondly, we observed that using the discrete time approach can lead to significantly better estimates since its growth condition allows to take more information into account by using a suitably constructed reference quantity.

Fig. 2 .
Fig. 2. Visualization of the bound B(•) in (5) resulting from our numerical computations in dependence of the parameter τ .

Fig. 3 .
Fig. 3. Visualization of α N,p in comparison to α T based on the estimates resulting from Theorem 5 with ˜ 2 and Theorem 2.

Fig. 4 .
Fig. 4. Visualization of α N,p in comparison to α T based on the estimates resulting from Theorem 5 with ˜ 1 and Theorem 2.