SAMPLED–DATA MODEL PREDICTIVE CONTROL: ADAPTIVE TIME–MESH REFINEMENT ALGORITHMS AND GUARANTEES OF STABILITY

. This article addresses the problem of controlling a constrained, continuous–time, nonlinear system through Model Predictive Control (MPC). In particular, we focus on methods to eﬃciently and accurately solve the underlying optimal control problem (OCP). In the numerical solution of a non-linear OCP, some form of discretization must be used at some stage. There are, however, beneﬁts in postponing the discretization process and maintain a continuous-time model until a later stage. This is because that way we can exploit additional freedom to select the number and the location of the discretization node points. We propose an adaptive time–mesh reﬁnement (AMR) algorithm that iteratively ﬁnds an adequate time–mesh satisfying a pre–deﬁned bound on the local error estimate of the obtained trajectories. The algorithm provides a time–dependent stopping criterion, enabling us to impose higher accuracy in the initial parts of the receding horizon, which are more relevant to MPC. Additionally, we analyze the conditions to guarantee closed–loop stability of the MPC framework using the AMR algorithm. The numerical results show that the proposed AMR strategy can obtain solutions as fast as methods using a coarse equidistant–spaced mesh and, on the other hand, as accurate as methods using a ﬁne equidistant–spaced mesh. Therefore, the OCP can be solved, and the MPC law obtained, faster and/or more accurately than with discrete-time MPC schemes using equidistant–spaced meshes.


1.
Introduction. In this work, we address the problem of constructing a sampleddata feedback control law to stabilize a nonlinear continuous-time system using a Model Predictive Control (MPC) technique. We propose a method, an adaptive time-mesh refinement algorithm, to solve the optimal control problems involved in computing the sampled-data MPC law. We show that the proposed approach can improve the computational efficiency and accuracy while preserving the guarantees of stability.
A sampled-data feedback system is obtained whenever the variables in a plant evolve in continuous-time and are controlled using a digital device. This is the most frequent situation in applications with some degree of complexity requiring an advanced control, which is not easy or practical to implement in an analogue device. The resulting sampled-data systems are often addressed, in practice and in the literature, using a simplified discrete-time model. Regarding the MPC literature, in the last decade, frameworks using discrete-time models clearly dominate (see e.g. [10,12,31]). However, some of the earlier fundamental theoretical works on nonlinear MPC (e.g. [8,17,34,35]) use continuous-time models. In the implementation of an MPC scheme, there is the need to perform, eventually, some form of discretization, or at least a finite parameterization, to solve the optimal control problems (OCP's) involved. Nevertheless, we argue that there are several advantages in maintaining a continuous-time model until later stages. That is, the continuous-time model is fed into the optimizer, which is responsible for the discretization occurring afterwards. The benefits of using an MPC framework in a sampled-data setting can be of an analytical nature or of a numerical nature and are discussed below.
Regarding the benefits of analytical nature, we highlight the following: • the plant model is usually obtained in continuous-time: when using a discretetime approximation from the beginning, it might be difficult to analyze the errors introduced in the discrete approximation; • discrete-time models generally have equidistant sampling instants: this might not be adequate for some nonlinear systems that have multiple time-scales of the dynamics, systems that can rapidly change behavior in short time intervals, or even when there is a discontinuity at some critical instant (e.g. impulsive control, [23], or discontinuous feedbacks [18]); • in path-following problems, the speed at which a reference trajectory is followed has to be adapted in real-time: this procedure is more complex in discrete-time ( [41,42]) than in continuous-time ( [7,14,44]); • the inter-sampling behavior of a continuous-time plant is not appropriately described by a discrete-time model: as a consequence, we might have hidden oscillations or constraint violations [24] that are not detected by the model.
The benefits of numerical nature are the focus of this article. They are essentially due to the fact that when solving OCP's, it is possible to devise numerical procedures that are more efficient using continuous-time models than using discrete-time models. One reason is that with continuous-time models we can exploit additional freedom in placing node points, manipulating the accuracy of the trajectory on several different time subintervals. Another benefit relates to the flexibility in manipulating the horizon length. It is well-known to the MPC community that increasing the prediction horizon might benefit, or even be necessary to guarantee, stability or recursive feasibility. Nevertheless, it is often reported that the computational time increases with the horizon length. While this is true for discrete-time systems with constant inter-sampling time, it is not necessarily the case when continuous-time is used. It is worth noting that some of the authors of works that vastly contributed to spread the use of discrete-time MPC ( [28,33,43]) have also explored advantages of frameworks with continuous-time models. In [39], it is shown that a continuous-time method is significantly more efficient than standard discretetime methods when solving constrained linear quadratic problems. In [27], Grüne et. al. provide an MPC scheme with a sampled-data implementation of the control law with the aim to minimize sampling errors. In [26], the authors propose a multi-step MPC scheme, which despite of being a discrete-time scheme, allows the optimization and feedback to be updated at different time instants. It is shown in [26] that such technique can maintain stability and quantifiable robust performance estimates while keeping the computational requirements low. When the optimal control problems maintain a continuous-time model, the numerical solution can use Adaptive Time-Mesh Refinement (AMR) algorithms, which have shown significant computational time savings when compared with solving the OCP with similar accuracy on equidistant time meshes (as the ones used in the discrete-time models). (For a description of AMR algorithms see [4,38,40,48] and the references therein.) Preliminary results on the application of the AMR technique to MPC, have been reported in [36,37].
We argue that this is the right time to start using more often MPC in the sampled-data setting, even in applications. Many theoretical questions of using sampled-data systems are well documented in the literature ( [15,19,22,27,32]). Also, there are now several readily available software packages to solve nonlinear OCP's ( [5,13,29]), which greatly facilitate the implementation of sampled-data nonlinear MPC schemes. Moreover, the choice of the discretization inter-sampling intervals can be automated, by pre-defining some error criteria, if adaptive timemesh refinement algorithms such as the one described here are used.
In this paper, we consider an adaptive time-mesh refinement (AMR) algorithm to solve a nonlinear OCP ( [38]) and extend it to include time-dependent levels of refinement, so that it can be applied in an MPC setting. When using direct methods to solve the OCP, an essential first step is to select an appropriate mesh to discretize the time interval and the dynamic equation. This selection -the number and location of the mesh node points -, has implications on the accuracy of the solution and on the efficiency of the method. Adaptive mesh-refinement algorithms increase iteratively the number of node points in the time subintervals for which the accuracy of the solution must be improved to attain a desired accuracy. The developed algorithm iteratively computes an adequate time-mesh that, in the end, satisfies a pre-defined bound on the local discretization error estimate of the obtained trajectories.
In an MPC context, the prediction can be interpreted in the sense of planning. When we make plans for the future, we establish planning strategies depending on the prediction horizon. When we plan our -professional or private -schedule, we plan with more detail and on a finer time scale in the shorter term: we might do a detailed hourly plan for the next day, a daily plan for the following week, and a draft of a monthly plan for the next year. Combining this idea with the refinement strategy, we obtain an adaptive time-mesh refinement algorithm which generates meshes with higher concentration of node points in the beginning of the prediction horizon and less concentration of node points in the end of the same interval. This way, we are bringing about the idea of having more node points where they are needed while keeping the overall number low. This is an important feature, since we want to increase the accuracy of the solution while maintaining a modest computational time.
The main features in the adaptation to an MPC context of the refinement strategy developed for an OCP are (a) the pre-defined levels of refinement are now time-dependent, which requires a generalization of the initially developed procedure; (b) since similar (time-shifted) OCP's are repeatedly solved, the developed procedure projects in the new time-mesh the solution of a previous OCP, to create a warm start; (c) in addition, we show that the MPC with the AMR algorithm preserves stability, provided that the design parameters of the controller satisfy a given stability condition.
This paper is organized as follows. In Section 2, we characterize the system model, the model assumptions, and the sampled-data feedback control setting. Section 3 introduces the optimal control problem used, while Section 4 discusses its numerical solution methods. In Section 5, we describe the adaptive time-mesh refinement (AMR) algorithm to solve OCP's and extend it to time-dependent levels of refinement. In Section 6, we provide the MPC framework and combine it with the extended AMR algorithm. In Section 7, we analyze the stability of the AMR-MPC scheme. In Section 8, we illustrate the application of the AMR-MPC scheme to a car-like system involving parking maneuvers and we report numerical results attesting the performance of the methodology developed. Finally, some conclusions are drawn in Section 9. Nomenclature.
Optimal Control and MPC frameworks P(x 0 ) Optimal Control Problem (OCP) starting at x 0 A 0 ⊂ R n set of admissible initial states J cost functional of the OCP V value function of the OCP (x,ū) optimal trajectory-control pair solving the OCP (x,ũ) admissible trajectory-control pair for the OCP (x * , u * ) closed-loop trajectory-control pair, numerically computed with MPC A (x 0 ) set of admissible trajectory-control pairs satisfying the constraints of P(x 0 ) k : 0, 1, 2, . . . , K sampling instants index t real time t k sampling instant Π = {t k } k≥0 sequence of sampling instants i : 0, 1, 2, . . . mesh node index s simulation time used in the prediction model and OCP s i mesh node instant S i mesh refining interval τ i extreme of the mesh refining interval j : 1, 2, . . . , J refinement level index ε(s) local discretization error estimate at time s ε j level of refinement j ε max stopping criterion for the refinement procedure r : 0, 1, 2, . . . , R refinement iteration index π r = {s i } i=0,1,...,Nr optimization mesh at r th refinement iteration set of control functions from [t 1 , t 2 ] to U , piecewise constant on π 0 .
2. The addressed system. In this section, we describe the system model, the sample-data feedback control and the model assumptions.
We consider a dynamical system for which the evolution of the state after some time t 0 is described by the following model: a.e. t ≥ t 0 , (1a) The state and input constraints are expressed as where X ⊂ R n is the state constraint set, and U ⊂ R m is the control constraint set. The initial state x 0 , measured at the time instant t 0 , is selected from the set of possible initial states X 0 ⊂ X, The trajectory solution to (1) is here denoted by the function t → x(t; x 0 , u), when is convenient to identify explicitly the dependence of the trajectory with respect to a particular initial state or control function.
We assume that the system satisfies the following set of hypotheses. H1: The sets X, X 0 and U are compact, contain the origin and f (0, 0) = 0. H2: The system (1)-(2) can be driven asymptotically to the origin, from X 0 , by some open-loop control. That is, for any x 0 ∈ X 0 , there exists a control function u with u(t) ∈ U , such that x(t; x 0 , u) → 0 as t → ∞. H3: The function f is continuous, and x → f (x, u) is Lipschitz continuous for any u ∈ U . It should be noted that these hypotheses depend only on the data of the given system. Later, we establish additional conditions depending on the design parameters to achieve certain desired properties.
2.1. Sampled-data feedback and stabilization. Our objective is to drive the state of the system to a desired equilibrium point, here considered without loss of generality to be the origin. This goal is to be achieved using sampled-data (SD) feedback control. In sampled-data (SD) feedback control (e.g. [9]), at any given time, the control is not a function of the state at that time, rather it is a function of the state at the last sampling instant. More precisely, consider a sequence of sampling instants Π := {t k } k∈N with a constant sampling period δ > 0 and t i+1 = t i + δ for all k ∈ N. Consider also a feedback law κ : R + × R n → R m and a function t → t Π that extracts from Π the last sampling instant before t, that is The control obtained via sampled-data feedback is defined for any t ≥ t 0 by u(t) := κ(t, x( t Π )). The solution of system (1)-(2) under this sampled-data feedback, starting from some initial state x 0 is denoted by the function t → x(t; x 0 , κ(t, x( t Π ))), or simply by t → x(t; x 0 , κ). The MPC technique, described below, combines naturally with this sampled-data feedback concept (see [16,18]).
To meet the goal of driving the state of the system to the origin through a sampled-data feedback, we must have open-loop controllability as is imposed in H2. The converse, whether asymptotical controllability implies sampled-feedback stabilization, is less obvious and was shown in [9] (in that result, though, somewhat different notions of controllability, stabilizability and sampled-data feedback were used).
3. The optimal control problem. The Model Predictive Control strategy involves solving a sequence of open-loop Optimal Control Problems (OCP) at each sampling instant t k , every time using the currently measured state of the plant x k . We consider the following OCP starting from an initial state x k ∈ X 0 : Here, in addition to the data already specified in (1)-(2), we use a fixed horizon T in R + , the running cost function L : R n × R m → R, the terminal cost function G : R n → R, and the terminal state constraint set X f ⊂ R n . The quantities T, L, G and X f are, in general, not given with the system plant (1) and can therefore be selected by the control engineer to determine the characteristics of the controller and the corresponding behavior of the controlled system. Hence, we call design parameters to such quantities.
We note that since the OCP is time-invariant, it is indifferent to consider and solve the problem in the interval [t k , t k + T ] or [0, T ], as long as the initial state x k is the correct one. For a time-varying formulation of this problem see e.g. [17]. For a general reference on optimal control see e.g. [25,46].
Regarding the domain of the optimization problem, we might conceptually consider that the control functions u are selected from the space of measurable functions. However, practicality imposes that these problems have to be solved numerically and so the domain of the optimization problem has to be finitely parameterizable (i.e. the control functions have to be representable in a digital computer using a finite number of parameters). The parameterization used might be based on piecewise constant control, which is the most common choice (e.g. [32,38]), and also is the one used here. Alternatively, the parameterization can be based on polynomials or splines ( [40]), based on bang-bang controls ( [19]), or any other representation of the control functions by a finite set of parameters.
4. Numerical solution of optimal control problems. The solution of nonlinear optimal control problems, in most cases, cannot be explicitly characterized analytically, but can only be obtained using numerical methods. Numerical methods to solve these problems include dynamic programming-based algorithms (solving the Hamilton-Jacobi-Bellman equation), indirect methods (applying the maximum principle to convert the problem into a boundary value problem), and direct methods (discretizing and solving the resulting static nonlinear programming problem).
To solve the OCP, we employ direct methods [3], in which the model is initially discretized in a given time-mesh, the controls are parameterized by piecewise constant functions, and then the problem is transcribed into a finite-dimensional nonlinear programming problem (NLP). The resulting NLP can be solved using one of the several widely available optimization solvers.
We stress that even when using a piecewise constant control parameterization, the intervals of constant control (defining what we call the time-mesh of the optimization solver) do not have to coincide with the inter-sampling intervals (defining the instants when each MPC iteration starts). Here, the optimization mesh π 0 = {s i } i∈{0,1,...,N0} defines the node points where the control function can change value; it contains the sampling instants in Π in the interval [0, T ] and possibility some other node points, i.e. π 0 ⊃ (Π ∩ [0, T ]).
We consider here that in the OCP above, we are seeking a minimizer by selecting a control function among the set U of functions from [0, T ] to R m taking values on U which are piecewise constant on π 0 and right-continuous functions. The domain of the optimization problem is the set of admissible pairs (x, u) comprising a control function u ∈ U and the corresponding state trajectory satisfying the constraints (3b)-(3f), defined to be If there exists an admissible process, i.e. if the above set is non-empty, we say that the OCP is feasible. It is convenient to define also the set of admissible initial states to be We say that a pair (x,ū) ∈ A(x k ) is a minimizer or an optimal solution to P( We note that the set U depends on the optimization mesh π 0 and so does the set A(x k ), the set A 0 , as well as the minimum value on the left-hand side of (4). We use the notation U(π 0 ) to make this dependence explicit when useful, and the notation U(π 0 , [t 1 , t 2 ]) when the control functions to consider have the domain restricted to some subinterval [t 1 , t 2 ]. The difference between the minimum value of the problem when the controls are selected from U(π 0 ) and when they are selected from the space of measurable functions is here denoted by optimality gap of π 0 .
The conditions to guarantee the existence of an optimal solution to the discretized optimal control problem are the same as the ones to guarantee the existence of an optimal solution to the equivalent NLP which, by the Weierstrass extreme value theorem, are simply the continuity of the cost functions L and G and compactness and non-emptiness of the domain.
The fact that the ordinary differential equation (1a) is discretized and then integrated numerically leads to a discretization error in the trajectory. We note that this is distinct and does not include the optimality gap due to the use of piecewise constant control functions. Given two consecutive mesh nodes s i , s i+1 , the local discretization error, here denoted by ε(s i+1 ), is the difference between the computed solution and the exact solution of the differential equation on [s i , s i+1 ] when the two solutions coincide at s i . The local discretization error is a measure used to guide the refinement in adaptive mesh-refinement (AMR) procedures (cf. [4,40]). An AMR procedure to decrease the local discretization error until reaching a pre-defined threshold is described in the next section.

5.
Adaptive Mesh Refinement (AMR) for OCP. When numerically solving an OCP, most frequently, regular time meshes having equidistant spacing are used in the discretization procedure. However, in some cases, these meshes are not the most adequate to deal with nonlinear behaviors. One way to improve the accuracy of the results, while maintaining a reasonable computational time and memory requirement, is to construct a mesh having different time steps. The best location for the smaller step sizes is, in general, not known a priori, so the mesh is refined iteratively. In a mesh-refinement procedure the problem is solved, typically, in an initial coarse uniform mesh in order to capture the basic structure of the solution and of the local discretization error. Then, this initial mesh is repeatedly refined according to a chosen strategy until some stopping criterion is attained. Several mesh refinement methods employing direct collocation methods have been described in the recent years ( [3,4,40,48]). The refinement method used here (a) permits several levels of refinement, obtaining a multi-level time-mesh in a single iteration; (b) it also permits different refinement criteria -the local discretization error of the primal variables, the discretization error of the dual variables and a combination of both criteria can be used in the refinement procedure; (c) it considers distinct criteria for refining the mesh and for stopping the refinement procedure -the refinement strategy can be driven by the information given by the dual variables and it can be stopped according to the information given by the primal variables.
As described in [38], there are advantages in choosing the error of the adjoint multipliers as a refinement criterion, namely in lowering the computational time, because the adjoint multipliers are solution to a linear differential equation system, which is easily solved with high accuracy. To decrease the overall computational time, the solution computed in the previous iteration is used as a warm start in the next one, which proved to be of major importance to improve computational efficiency. This adaptive strategy leads to results with higher accuracy and yet with lower overall computational time, when compared to results obtained by meshes having equidistant spacing, as is the case when using discrete-time models from the beginning.
5.1. The AMR algorithm. We describe here a specific multi-level AMR algorithm proposed in [38], where we refer to for further details. The adaptive timemesh refinement procedure starts by discretizing the time interval [0, T ] in an initial coarse mesh where the NLP problem associated to the OCP is obtained and solved in order to capture the main structure of the solution and to estimate the discretization error.
To simplify presentation, we consider that this mesh is used both to define the piecewise constant control intervals and to discretize the differential equation (1a).
In the second step, we need to identify the time intervals that verify the refinement criteria. That is, we partition the time domain into a set of intervals for which the discretization error estimate is above or below a certain threshold ε max .
In the third step, the intervals S i are refined taking under consideration different levels of refinement, i.e., they are divided into smaller intervals according to userdefined levels of refinement ε 1 , ε 2 , . . . , ε J with ε 1 = ε max . We consider that an interval S i,j is in the j th level of refinement if where ε(s) is an estimate of the discretization error. The algorithm, while constructing the next mesh π 1 , adds more node-points to intervals in higher refinement levels and it adds less node points to those in lower refinement levels (see Figure 1), following a block-structured adaptive mesh refinement strategy where no nodes are removed, just added, i.e. π 1 ⊃ π 0 . By defining several levels of refinement in a single iteration, we generate a multi-block time-mesh.
In the example of the figure, we are subdividing each interval in the first level into N = 2 subintervals (adding N − 1 additional nodes to each interval). Each time interval in refining level j is subdivided into N j subintervals (by adding N j −1 additional nodes). The procedure is then repeated for iteration r = 1, 2, . . ., until a stopping criterion is achieved (ε(s) ≤ ε max for all s ∈ [0, T ]). Refinement and stopping criteria. In order to proceed with the mesh refinement strategy, we have to define some refinement criteria and a stopping criterion.
For the stopping criterion, we consider a threshold for the discretization error of the trajectory (primal variables) (ε x ), and the algorithm computes its maximum estimate value (ε x ) ∞ . For the refinement criteria, we could also consider thresholds for (ε x ) for different time instants, as is done in e.g. [4]. However, as shown in [38], similar error information can be extracted from the dual variables and such information can be computed much faster. Therefore, we set as refinement criteria levels of discretization error (ε q ) for the adjoint multipliers.
The adjoint multipliers satisfy a differential equation system which is defined by the Maximum Principle. Such differential equation is linear and can therefore be easily solved in a fast way and with high accuracy. In order to select which time intervals should be further refined at each refinement step, the algorithm computes the local error of the multipliers for each time instant s: where q MP is a solution to the differential equation system given by the Maximum Principle [21,46], and q NLP is the multiplier numerically obtained by solving the associated nonlinear programming problem (NLP) using a numeric solver. We may consider that this algorithm is somehow a hybrid between a direct and indirect method, since it involves an initial discretization as in a direct method, but it is driven by information obtained by the maximum principle, as in an indirect method. We note that, currently, with the existing automatic differentiation software, obtaining the adjoint equation is an easy process that can be automated (see e.g. CasADi [2] or ICLOCS2 [1] for optimal control software with algorithmic differentiation).

5.2.
Extended AMR algorithm. We extend the algorithm described in section 5.1, by considering time-dependent refinement thresholds (ε max (s)) for the mesh refinement algorithm. This procedure adds more node points and thereby obtains a more accurate solution in the time intervals that correspond to tighter refinement thresholds, see Fig. 2. The time-dependent refinement thresholds can be defined within subintervals of the interval [t k , t k + T ] as follows where 0 < β 1 < β 2 < . . . < β J −1 < 1.
For example, we can define a time-dependent refinement thresholds with 3 levels: for j = 1, . . . , 3 with ε 1 = 10 −5 . Considering β 1 = 10%, β 2 = 30%, we obtain: 6. Model predictive control. The MPC algorithm constructs a sampled-data feedback control by concatenating the initial parts of the optimal control functions obtained as solution to each problem P(x k ) along each t k ∈ Π. Starting at t k = t 0 , the following four steps are implemented at each MPC iteration k. , t > δ is disregarded); 4. Repeat for the next sampling time instant t k ← t k + δ and the next MPC iteration k ← k + 1. The pair (x * , u * ), defined for t ∈ [t 0 , +∞), denotes the sampled-data control resulting from applying the MPC algorithm and the corresponding state trajectory.
6.1. Model predictive control coupled with the extended AMR algorithm. The algorithm extension described in Section 5.2 is of relevance in the MPC context, since it is desirable to have a solution with higher accuracy in the initial part of the horizon.
When applying the MPC procedure to solve an OCP, at each time instant t k we compute the solution in the interval [t k , t k + T ] but we just implement the openloop control in its initial part, until t k + δ. Therefore, taking under consideration the planning strategy discussed above, it would be an improvement if we have an adaptive time-mesh able to cope with this feature, i.e., a time-mesh that is highly refined in the initial part of the time horizon and it is coarser in the final remaining part of the same horizon. Then, we would implement the control on the time interval [t k , t k + δ] computed with high accuracy in a refined mesh. For the remaining time interval we have a less accurate estimate of the solution. In the example bellow, with T = 10 and δ = 2, we might have 10 −4 , s ∈ ]t k + 0.4T, t k + T ] Following the described strategy, we obtain an adaptive time-mesh refinement algorithm which generates meshes with higher concentration of node points in the beginning of the interval [t k , t k + T ] and less concentration of node points towards the end of the same interval, enforcing the idea of having more nodes point where they are needed and keeping a low overall number of node points. This is an important issue because we want to increase the accuracy of the solution without compromising low CPU times.
We start the MPC procedure in the usual way but considering an extended adaptive mesh refinement strategy to solve the OCP. We discretize the time interval [0, T ] and we solve our OCP in open-loop. Then, we implement the control in the first sampling interval. When starting the next MPC step, we apply the time-mesh refinement strategy in order to find the best mesh suited to the solve the OCP in the second sampling interval. We repeat this technique until the MPC procedure ends. In the MPC algorithm, step 2 is modified as follows: 2. (a) Select the intervals S i,j to be refined according to the time-dependent refinement thresholds ε 1 , ε 2 , . . . , ε J and generate a new time-mesh; (b) Solve the OCP P (x k ) obtaining the optimal control functionū : [0, T ] → R m as well as the corresponding trajectoryx : [0, T ] → R n , in the new time-mesh; This AMR-MPC algorithm is illustrated in Fig. 3.  Warm start. Since we are solving a sequence of open-loop OCP's, to decrease the CPU time when going from the problem in [t k , t k + T ] to the problem in [t k + δ, t k + T + δ], the solution of the previous is used as a warm start for the problem. To create this warm start, the solution obtained in [t k , t k + T ] is projected in the new mesh in [t k + δ, t k + T + δ]. This action proved to be vital in the decreasing of the overall computational time. In particular, we notice that the number of iterations of the NLP solver remains within the same order of magnitude along the procedure, even when the number of nodes is significantly increased.
Another strategy that could be pursued to further accelerate obtaining the next solution would be to use additional pre-computed information from the optimizer (e.g. Hessian matrices, gradients); see [11].
7. Stability of the AMR-MPC. In this section, we provide a stability guarantee for the AMR-MPC strategy. The main result asserts that the system resulting from the application of the MPC algorithm is stable, as long as the design parameters are selected to satisfy the Stability Condition below. This condition should be satisfied for initial mesh π 0 containing the sampling instants in [0, T ]. SC: Sufficient condition for Stability. The design parameters T, L, G and X f satisfy: 1 The set X f is a subset of X, is closed, and contains the origin. The function G is Lipschitz continuous and positive definite. The function L is continuous and there exists a function M : R n → R + which is continuous, positive definite and radially unbounded, such that L(x, u) ≥ M (x) for all u ∈ U . 2 The horizon T is such that X 0 is contained in A 0 , when controls from U(π 0 ) are used. 3 There exists a control law k f : and Remark 1. Condition 1 on the objective functions is nonrestrictive for stabilizing and tracking problems. A typical choice for L is a quadratic cost. For the terminal set X f , some possible choices are the set X, a subspace of X, or an ellipsoid centred at the origin. This last choice is convenient for systems for which the linearization is stabilizable (see [8]). In case the aim is the convergence to a target point x # different from the origin, in which the pair (x # , u # ) is an equilibrium pair (i.e. f (x # , u # ) = 0), then condition 1 should be adapted as follows: The set X f is a subset of X, is closed, and contains x # . The function G is Lipschitz continuous and x → G(x − x # ) is positive definite. The function L is continuous and there exists a function M : Condition 2 guarantees the feasibility of the first OCP which implies the feasibility of all subsequent problems. It is often ensured by selecting a sufficient large horizon T . Condition (SCa) is similar to a decrease condition of a Lyapunov function. (Note that the pointwise version of this condition ∇G(x(s; for a.e. s ∈ [0, δ], immediately implies (SCa) after integration on [0, δ]). However, it can be much easier to verify on a smaller set X f than globally.
Condition (SCb) is an invariance condition. After the end of the horizon, the trajectory under the feedback κ f has to remain within X and return to X f after time δ. Alternatively, if the set X f is invariant under κ f , i.e. x(s; x f , k f ) ∈ X, for all s ∈ [0, δ], then this condition can replace both conditions in (SCb) since X f ⊂ X.
Condition (SCc) specifies the admissible controls. Note that from the definition of U it implies that k f (s, The stability result making use of this condition follows.
Theorem 7.1. Assume the system satisfies hypotheses H1-H3. If the design parameters T, L, G and X f satisfy the stability condition SC, then applying the AMR-MPC strategy starting from any x 0 ∈ X 0 and with some initial mesh π 0 ⊃ Π we have: 1. all optimal control problems involved in the AMR-MPC strategy, P(x k ) for all k ≥ 0, are feasible and have a minimum. 2. the closed-loop trajectory x * is asymptotically attracted to the origin, that is x * (t) → 0 as t → +∞.

Remark 2.
We note that this result is just a nominal stability result and, as such, the trajectory x * is the closed-loop trajectory obtained numerically by applying the AMR-MPC strategy. That is, it is assumed that the evolution of the plant trajectory follows exactly the behavior described by the model and is computed using the described numerical scheme.
7.1. Proof of the stability result. This theorem is based, with changes related to AMR, on arguments to establish stability of sampled-data MPC in [8] and [17,Thm. 3]. In what follows we establish, respectively, recursive feasibility, existence of an optimal solution, uniform decrease of the value function, and asymptotical convergence.
Lemma 7.2 (Recursive feasibility). Let x 0 ∈ X 0 and x k = x(kδ) for k ∈ N when the nominal model is followed and the AMR-MPC strategy is applied starting with initial mesh π 0 ⊃ (Π ∩ [0, T ]). Then, P(x k ) is feasible for all k ∈ N. Moreover, an optimal solution to P(x k ) exists.
Proof. For the initial MPC iteration, k = 0, t = 0, the initial state x(0) = x 0 is selected from X 0 and X 0 ⊂ A 0 when π 0 is used. By definition of A 0 , problem P(x 0 ) is feasible. Let the pair (x,ū), defined for s ∈ [0, T ], be an optimal solution to P(x 0 ). Then, assuming the nominal model is followed, we have that, for iteration k = 1, t = δ, the new initial state is x 1 =x(δ).
We now show that P(x 1 ) is also feasible. This is accomplished by defining a pair (x,ũ), also in the interval s ∈ [0, T ], which is guaranteed to be admissible if constructed as an extension of (x,ū) to the interval t ∈ [t k +T, t k +T +δ] as follows.
where κ f is the control law in SC. See Fig. 4 The AMR procedure is such that the time-mesh nodes in [δ, T ] of a certain iteration k are maintained in the next iteration k + 1 in [0, T − δ] (i.e. no mesh nodes are removed in our AMR and all meshes contain the sampling instants in the interval [0, T ], π r ⊃ . . . ⊃ π 0 ⊃ (Π∩[0, T ])). The controlũ can be constructed in the new mesh. Also, from (SCc) we have thatũ(s) ∈ U for s ∈ [T − δ, T ]. Therefore,ũ Figure 4. Construction of the (extended) admissible controlũ with Π = {t k } k∈N , t k = kδ, and with π r = {s i } i∈0,1,...Nr , s i = iδ/2.
Since the problems are time-invariant, our construction ensures that after a feasible problem P(x k ) at MPC iteration k, it follows another feasible OCP P(x k+1 ) at iteration k + 1. Using an induction argument, we show that all OCP's for k ≥ 0 are feasible.
To show the existence of an optimal solution to the discretized P(x k ), we have to show existence of an optimal solution of the equivalent NLP. That is, by the Weierstrass extreme value theorem, the minimum exists if in the NLP (i) the objective function is continuous; (ii) the domain is nonempty; (iii) the domain is compact. Condition (i) is guaranteed by the continuity of L and G. Condition (ii) holds from the feasibility of the OCP established above. Condition (iii) is implied by the compactness of U , X and that its images by a continuous function f are also compact. Therefore, the NLP and thereby the OCP has a minimum.
A key step to show the convergence of the trajectory to the origin is to establish that the value function of the optimal control problems decreases uniformly along the sequence {x k } k∈N . This, in turn, can be shown using SC, mainly with conditions (SCa),(SCb), and (SCc).
We define the cost functional and value function of our OCP as follows. The cost functional is defined to be Having guaranteed the existence of an optimal solution (of the discretized OCP), the value function is We have the following lemma.
Since the pair (x,ũ) is admissible for P(x 1 ), but not necessarily optimal, we have We also have Computing the difference, we obtain where in the last inequality we have used (SCa). Finally, recalling the condition on M in SC we have Since the OCP's are time-invariant and x * is the concatenation of the initial part of the optimal trajectories of the sequence of OCP's, we have From the previous lemma, we can then write that for any k ≥ 0 Since V (x 0 ) is finite, we conclude that the function k → V (x k ) is bounded and then that k → kδ 0 M (x * (t))dt is also bounded. Therefore k → x * (kδ) is bounded and, since f is continuous and takes values on bounded sets of (x, u), t →ẋ * (t) and t → x * (t) are also bounded. All the conditions to apply Barbalat's lemma ( [45,20]) are met, yielding that the trajectory asymptotically converges to the origin.
Note that we have just confirmed what is usually known in the literature as attractiveness (i.e. x * (t) → 0 as t → +∞). Stability in the sense of Lyapunov is not established here. It would be achieved if V is continuous, for example. However, in the application explored below, a car-like system with minimum turning radios R, it is not possible to attain stability in the sense of Lyapunov. To see this, note that when the system is at the target point, an arbitrarily small perturbation in y would require a significantly long path (with length bounded below) to return to the target. 8. Application. To test and to validate the proposed algorithm, a problem involving parking manoeuvres in a constrained space is solved. Consider a car-like system with a minimum turning radius r min . At a given time instant s, we define the state x(s) = (x(s), y(s), ψ(s), v(s)) and the control u(s) = (a(s), c(s)) where (x(s), y(s)) is the position of the midpoint of the axle connecting the rear wheels, ψ(s) is the heading angle, v(s) is the speed, a(s) is the acceleration and c(s) is the curvature (see Figure 5). The objective is to drive a car-like vehicle from an initial position x 0 to a target point x # in a constrained environment. This constrained environment, implemented as a pathwise state constraint, represents a parking lot where the vehicle is parked and the target area is the exit (see Figure 6).
We have chosen this application to show how MPC can overcome nonholonomy challenges ( [6,30]) since it involves planning, not just reactive control; it can generate required nonlinear, discontinuous feedback; and it illustrates that the sampleddata MPC framework combines well with sampling-feedbacks [18]. MPC can also overcome constraints challenges since it is known to be a (if not the main) technique to deal appropriately with constraints.
The optimal control problem (P CP ) for this system on s ∈ [0, T ] is written as follows: x y l ψ Figure 5. Car-like system geometry (iii) the initial state constraint (iv) the terminal constraint and (v) the path-wise state constraint The state constraint set X is defined to be the set of states x(s) satisfying where is the target point, which is an equilibrium state with u = 0. The design parameters L, G, X f and T are selected in such a way that the stability condition SC is satisfied, guaranteeing convergence to the equilibrium point x # . Note that we are not changing the coordinates to make the target state coincide with the origin (see Remark 1).
The terminal set X f is the set of points (x, y) in the x axis, heading towards x # with positive speed (see Figure 6), defined to be The running cost is simply a quadratic cost with the control effort weight α ≥ 0. The terminal cost is the cost to reach the target point in timet from a point in X f maintaining constant speed (of at least v min > 0) and zero curvature, that is G (x(T )) = t 0 L (x(s), u(s)) ds, for which an explicit formula is This formula is obtained noting that with u(s) = (a(s), c(s)) = (0, 0), the time to reach the origin ist Finally, selecting the horizon T to satisfy it can be verified that the stability condition SC holds (see [17,Section 4.3] for details of the verification for a similar system).

Numerical implementation.
To test the proposed algorithm, we have implemented it in MATLAB R2016a combined with optimal control solver ICLOCSversion 0.1b ( [1]) using IPOPT as the NLP solver ( [47]). The problems were solved in a computer with a Intel™ Core © i7-4770K CPU @3.50 GHz. We consider the target point the end-point constraints the pathwise state constraint where We consider the horizon T = 10, the sampling period δ = 2 s, and v min = 0.1.

8.2.
Numerical results for the AMR-OCP. We start by discussing the results obtained when solving optimal control problem (P CP ) using the AMR procedure. We consider as stopping criterion the bound on the estimate of the trajectory discretization error and as refinement criteria for the different levels, we consider the vector ε x = 1, 5, 10, 10 2 ε max x . By the AMR algorithm, each time interval of refining level j, is subdivided into N j subintervals by adding N j − 1 additional nodes. We consider N = 3. The problem is solved in the following three meshes: a) π 0 the initial coarse mesh with 21 equidistant nodes; b) π AMR obtained by the multi-level time-mesh refinement strategy; c) π F considering equidistant nodes with the lowest time step of π AMR .
The problem is first solved using the initial coarse mesh π 0 . In Fig. 7, we can see the path obtained in this initial solution. The car starts at position (0.5, −3), heading downwards. It moves in reverse until reaching the point of maximum y, then moves forward turning left until attaining the target.
When we analyze the discretization error, see Fig. 8, we verify that it is above ε 1 = ε max x between the 2nd and 8th nodes, and it is slightly above ε 3 at node 6. Therefore, we apply refining level 3 in the two intervals adjacent to node 6, dividing each into N 3 = 27 subintervals and we apply refining level 1 in the remaining 4 intervals between the 2nd and 8th node, diving each one of those into 3 subintervals. In total, the refining procedure adds 60 new nodes, 4 * (N −1)+2 * (N 3 −1), forming a new mesh π 1 with 82 nodes (one additional node is added by the procedure to compensate for truncating each interval size to maximum precision) and with the smallest subinterval of dimension 1/54 seconds. Table 1 shows for each mesh the number of nodes, the smallest time step of each mesh, the number of iterations needed to solve the NLP problem, the maximum absolute local errors of the trajectory, and CPU times for solving the OCP problem and for computing the local error as well.
The OCP is then solved in the new mesh π 1 . Despite the significant increase in the number of mesh nodes, the optimizer takes about the same number iterations (coincidentally exactly the same) and similar time to obtain a solution when compared to the first run using π 0 (see Table 1). This is because the previous solution was used as a warm start of the optimizer.
The discretization error estimate is again analyzed and it is verified that it is already below the threshold ε max x . Therefore, there is no need for further refining iterations and the AMR algorithms terminates with final mesh π AMR equal to π 1 . In Table 1, the times and number of iterations of the optimizer for π AMR are reported as the sum of the values obtained in all previous iterations. Fig. 9 depicts the optimal path computed in the final mesh π AMR and Fig. 10 plots the trajectory and the associated control. We can see that the π AMR mesh has higher resolution precisely in the time subintervals where the discretization error was higher in the first step. These time intervals coincide with the ones where the speed attains a higher absolute value. To compare the performance of AMR strategy for a similar level of accuracy, we tried to solve the OCP in a mesh π F which is an equidistant mesh with nodes spaced with the smallest time step of π AMR , which is 1/54 seconds. The mesh obtained this way has 541 nodes. With this mesh, we could not obtain a solution from the optimizer within the maximum of 5000 iterations of the NLP when cold started (with no meaningful initial guess) with the initial state x(0) = (0.5, −3, − π /2, 0). We considered then the problem with initial state x(0) = (0.5, −3, π /4, 0.6), which is considerably simpler problem since the car starts already moving towards the origin. The results from this run are reported in the last line of Table 1.
Despite the fact that a simpler problem is being solved and as, as expected, the maximum absolute local errors are of the same order of magnitude, it took considerably more time and more NLP iterations to obtain a solution in the equidistant fine mesh π F . The mesh π AMR has about 15% of the nodes of π F , and computing the solution using π AMR takes only about 12% of the time needed to get a solution using π F , causing significant savings in memory and computational cost.
Moreover, the fact that the AMR algorithm was able to solve a problem that could not be solved in an equidistant fine mesh shows evidence of its robustness.

8.3.
Numerical results for the AMR-MPC. To implement the AMR-MPC scheme, we consider a sampling period of δ = 2 seconds and we consider N = 2 subintervals in each refinement. We also set the discretization error estimate threshold to ε max and we consider three levels of refinement (the tightest level on the initial 20% part of the horizon, coinciding with the sampling period, the second level on the following 20%, and the third level on the remaining part of the horizon) with the following limitsε x , t ∈ ]t k + 0.4T, t k + T ] Before entering the MPC loop, the OCP is solved in the initial coarse mesh π 0 . For the solution obtained in such mesh, we plot the discretization error against the new time-varying threshold, Fig. 11, and we can see where the refinement is expected to occur in the first MPC iterations.
We apply the AMR-MPC scheme. Table 2 shows for each MPC iteration and each AMR mesh the information about the number of nodes, the smallest time step, the number of iterations needed to solve the NLP problem, the maximum absolute local errors of the trajectory, and computational times for solving the OCP problem and for computing the local error. We observe that the MPC iterations are considerably faster than the initial cold-started OCP. The refinement procedure is only activated in the first MPC iteration. The second MPC iteration still uses part of a refined mesh, but after that the algorithm detects that there is no need to refine and uses a coarse mesh. All error estimates are well below the pre-defined threshold ε max x . Figure 12 shows that the path resulting from the AMR-MPC scheme is equal or very similar to the path obtained when solving the open-loop OCP.
We also observe in Fig. 13 that the solution to the several OCP's overlap (different colors are used for the solution to each OCP; the slight difference in the interval from 2 to 4 seconds, corresponding to the 2nd MPC step, is due to different refinement thresholds between the 1st and 2nd MPC step). This is expected by the choice of objective functions we have done. Since the terminal cost is the cost to go till the target and after that the cost is zero, we have that W (x(T )) = T +t T L(x(s), u(s)) ds = ∞ T L(x(s), u(s)) ds.
Hence, the solution to each OCP coincides in its horizon with the solution to the OCP with infinite horizon objective function with the terminal constraint set imposed also at t = T ∞ 0 L(x(s), u(s)) ds. 9. Conclusions. In this article, we develop an extended adaptive time-mesh refinement (AMR) algorithm to solve open-loop optimal control problems (OCP). The algorithm provides local mesh resolution refining only where it is required, thereby guaranteeing accuracy of the procedure in an efficient way. In the extension proposed, we consider a time-dependent stopping criterion for the mesh refinement algorithm with different levels. The extended algorithm is particularly useful to solve the OCP's involved in MPC schemes. Such use is explored in this article. We propose an MPC framework incorporating the extended AMR algorithm. We show that stability of the system resulting from applying the AMR-MPC strategy can be guaranteed. The extended AMR and the AMR-MPC algorithms are tested in a car-like system with limited turning radius, performing manoeuvres in a constrained space. In the end, the OCP's are solved within MPC with an adapted mesh with local mesh resolution, which has less nodes in the overall procedure, yet having maximum absolute local error of the same order of magnitude when compared against a refined mesh with equidistant-spacing. The fast response of the algorithm, even with tight accuracy limits, shows its ability to solve complex problems in realtime. The application demonstrates the advantage of the proposed adaptive mesh strategy, which leads to results obtained approximately as fast as the ones given  Figure 11. Discretization error in the coarse mesh and the MPC refining levels  This framework allows to use continuous-time plant models directly: the discretization procedure is automated and there is no need to select a priori an adequate discretization time-step.
Another advantage of the use of the AMR-MPC is that even if the optimization procedure is forced to stop in an early stage, as might be the case in real-time, it can still provide a solution, although it might be a less accurate one.
The AMR-MPC strategy proposed here might also be useful in situations that are not explored here but are worthy of further research. These include systems with