Time adaptivity in model predictive control

In this article we consider a feedback control approach for a linear parabolic equation using a Model Predictive Control (MPC) framework. Each iteration of the MPC algorithm considers an open loop problem. Here, we propose the use of time adaptive strategies based on a-posteriori error estimations for the residual of an elliptic equation, which is equivalent to the optimality system. The aim is to: (i) improve the accuracy and efficiency of the method, (ii) compute the prediction horizon and (iii) determine an adaptive temporal discretization which is related to the optimal solution. Numerical tests are carried out and show the advantages of our approach.


Introduction
Optimization with PDE constraints is a fundamental topic due to its large requests in many applications like e.g. engineering, medical science, economical and industrial issues etc. Many techniques have been developed in the last decades using, mainly, two different approaches such as open and closed loop to compute the optimal control. The first one is based on the Pontryagin's maximum principle and/or Lagrange techniques and provides a control as function of time for unsteady problems. This has been successfully applied to many problem settings for ordinary and partial differential equations and we refer to the monograph [18] for a general overview.
The second approach is based on the dynamic programming principle and provides a control in closed form as a function of the current state. Its aim is the possibility to realize a control which depends on the position of the system and the ability to react to perturbations due to e.g. external influences or errors in the measurements. Clearly, it has a major advantage in applications. This idea goes back to Bellman in late '50 through the solution of an Hamilton-Jacobi-Bellman (HJB) equation. However, this equation does not admit analytical solutions in many cases and its numerical approximation is rather complicated since it suffers from the curse of dimensionality. This approach requires the solution of a nonlinear hyperbolic PDE with a high-dimensional spatial variable when dealing with large systems such as e.g. PDEs. Although a complete list of works is out of the scopes of this work, we refer the interested reader to [6] for the detailed mathematical analysis of the problem and to [11] for the numerical approximation and the references therein. For large-scale problems with a HJB approach we refer e.g. to [2,20].
An alternative method to compute feedback control is given by the Model Predictive Control (MPC) approach. It is based on repeated solutions of open loop optimal control problems. The first part of the resulting open loop input signal is implemented and the whole process is repeated successively. This approach is also known as moving horizon control or receding horizon control. We refer to the monographs [13,23] for more details on this topic. The prediction horizon plays a crucial role in MPC algorithms. For instance, the quasi infinite horizon nonlinear MPC (NMPC) allows an efficient formulation of NMPC while guaranteeing stability and the performances of the closed loop as shown in [8,19] under appropriate assumptions. For the purpose of this paper, we will use a different approach since we will not deal with terminal constraints.
In the present paper, we study the optimal control of a linear time varying parabolic equation. In linear MPC, linear models are used to predict the system dynamics and consider linear constraints on the states and inputs. Note that even if the system is linear, the closed loop dynamics are nonlinear due to the presence of constraints. The novelty of the paper lies in the inclusion of adaptive concepts for the prediction horizon, which are based on a-posteriori error estimates for the optimal state solution. The idea of adaptivity in MPC has been investigated in e.g. [15,21] where the authors took advantage of the structure of the problem using Lyapunov functions and/or the turnpike property. The turnpike property (see e.g. [24]) is often a key tool to prove asymptotic stability of the MPC method and to find the minimal prediction horizon (see e.g. [7,14,17,22]).
Our method aims to present a general framework for linear quadratic control problems where the prediction horizon and the time discretization is computed automatically using an error estimator. Our approach works as follows: in a first step, we rewrite the optimality conditions as a second order in time and fourth order in space elliptic equation for the state variable. Introducing an auxiliary variable, we can split this equation into a coupled system of two second order in space elliptic equations. The advantage of considering a reformulation of the optimality conditions is, in particular, that we can apply classical concepts from residual based a-posteriori error control. This allows to construct a suitable time grid for the state which is related to the optimal state solution. The idea is based on [12], and now transferred for a mixed formulation where the a-posteriori error estimate is obtained from a semi-time discrete mixed form. The method assumes that the structure of the temporal grid is not sensitive against changes in the spatial resolution and verified heuristically by numerical examples (see [3,4]). The use of an a-posteriori error estimator leads us to two different approaches.
1. An adaptive horizon is carried out using the a-posteriori error estimates before starting the MPC algorithm. We prescribe the number of time instances we want to deal with and the algorithm determines where to collocate them. Then, choosing the number of instances we get an adaptive prediction horizon.
2. An adaptive time grid in each subinterval of the MPC method. This way allows to compute the solution in the most important time instances and to obtain a grid which is not equidistant.
Those two methods are fully general and applicable without any particular requirements. We have presented preliminary work towards time adaptivity in MPC together with model reduction utilizing proper orthogonal decomposition (POD) in talks and posters in various workshops and conferences, e.g. ICCOPT 2019, FGS 2019. While editing the final version of our manuscript the work [16] appeared, where goal-oriented adaptation concepts are proposed to adaptively solve the optimization problems appearing in every step of the MPC algorithm. However, the a-posteriori concepts proposed there differ from our approach which relies on residual based a-posteriori error analysis for the elliptic space-time reformulation of the optimality systems appearing in every step of the MPC algorithm. The outline of this paper is as follows. In Section 2, we present the optimal control problem together with the optimality conditions. In Section 3 we describe the reformulation to a second order in time and fourth order in space elliptic equation as well as a mixed variational form. Further, we derive an a-posteriori error estimate for a semi-time discrete form. In Section 4, we recall the basic idea of the MPC method. In Section 5, we propose the novel time adaptive schemes in MPC which are an offline and an online approach. Finally, numerical tests are discussed in Section 6 and conclusions are made in Section 7.
We use the notation H −1 (Ω) for the dual space of H 1 0 (Ω) and denote ·, · H −1 (Ω),H 1 0 (Ω) as the duality pairing of H −1 (Ω) with H 1 0 (Ω). By | · | H 1 (Ω) we denote the H 1 -seminorm given by |u| H 1 (Ω) = ∇u L 2 (Ω) for u ∈ H 1 0 (Ω). We recall that the Poincaré constant is given by the smallest number c p > 0 such that the Poincaré inequality . For a given Banach space X, we denote by L 2 (0, T ; X) the space of square integrable functions from (0, T ) to X with norm u L 2 (0,T ;X) := ( Note that for a given function g in spacetime, we use the short hand notation g(t) to indicate the time dependency and drop the space argument.

Optimal control problem
In this section, we describe the distributed optimal control problem which we consider throughout this work. The governing dynamics for the state y are given by a linear parabolic partial differential equation of the form where ν > 0 is a given constant, f is a given source term and y 0 is a given initial state. The function u will in the following act as the control. The weak form of (1) reads as follows: for given f ∈ L 2 (0, T ; Ω), y 0 ∈ L 2 (Ω) and u ∈ L 2 (0, T ; Ω), find a state y ∈ W (0, T ; H 1 0 (Ω)) satisfying y(0) = y 0 such that holds for all v ∈ H 1 0 (Ω) and almost everywhere in (0, T ]. For existence of a unique solution to (2), we refer to e.g. [10, §7.1.2, Theorems 3 and 4].
We indicate through the index of the weak solution y [u,t0,y0] to (1) its dependency on a given control u, initial time t 0 and initial value y 0 . This allows to introduce a reduced cost functionalĴ defined byĴ Thus, the optimal control problem (4) can be formulated as min u∈L 2 (0,T ;Ω)Ĵ (u, t 0 , y 0 ).
3 Reformulation of the optimality system and time adaptivity 3.1 Reformulation of the optimality system Following along the lines of [12], we can reformulate the optimality system (1)-(7)-(8) as an elliptic equation of fourth order in space and second order in time involving only the state variable y. The adjoint state p as well as the control u are not present in this equation. In particular, it is a two-point boundary value problem in space-time given by in Ω, y(0) = y 0 in Ω.
We note that for ν = 1 and f ≡ 0 this setting coincides with the setting considered in [12]. Under higher regularity assumptions on the data, the following theorem shows that the optimal state y of (1)-(7)-(8) fulfills the elliptic equation (9) a.e. in space-time.
Proof. The proof follows along the lines of the proof of [12,Theorem 2.7] and uses differentiation and insertion of the equations (1)-(7)- (8).
Let us homogenize (9). For this, let g be a function which fulfills the boundary conditions as well as initial and end time conditions of (9) and is sufficiently smooth. Let y satisfy (9). We defineỹ := y − g and arrive at in Ω, whereỹ Now, let us derive a weak formulation of (10). For this purpose we introduce the function space We introduce the following symmetric bilinear form Existence of a solution to (12) and its relation to a solution to (9) is shown in the following theorem.
Theorem 3.2. Let y denote a solution to (9) and let g be a function which fulfills the boundary, initial and end time conditions in (9) and is sufficiently smooth. Then,ỹ = y − g is a solution to (12). On the other hand, ifỹ is a solution to (12) and the assumptions of Lemma 2.1(ii) are fulfilled, then y =ỹ + g satisfies (9) a.e. in space-time.
Proof. Assume y is a solution to (9). By Green's formula and integration by parts it is straight forward to prove thatỹ = y − g satisfies (12). The other direction follows vice versa.
In order to show equivalence of the optimal control problem (4) to the weak formulation of (9) it remains to prove uniqueness of a solution.

Mixed formulation
In order to use piecewise linear, continuous finite elements for discretization and avoid the construction of finite element subspaces in H 2 (Ω), we introduce an auxiliary variablẽ w := −ν∆ỹ. This allows to write (10) as a coupled system inỹ andw as We introduce the function spaces Y := {v ∈ H 1 (0, T ; H 1 0 (Ω)) : v(0) = 0 in Ω} and W := L 2 (0, T ; H 1 0 (Ω)) and the product space X := Y × W . Let us define the following bilinear form and linear form Definition 3.2. The weak formulation of the mixed formulation (13) is given by: By analogy with Theorem 3.2 and Theorem 3.3 it can be shown that the mixed variational form (14) admits at most one solution and that the pair (ỹ,w) withỹ denoting the unique solution to (10) andw := −ν∆ỹ is a solution to the mixed variational form (14). This means that the unique solution to (10) defines the solution to the mixed variational form (14).

Note that it is
For this reason, we define an energy norm associated with the bilinear form A M by

A-posteriori error estimate for the semi-time discrete mixed form
Let us now consider a semi-time discretization of (14) with respect toỹ while the variablẽ w is kept continuous. We introduce a time grid where P p denotes the space of polynomials of degree up to p. We set With arguments similar to those used for (14) we may show that problem (15) admits a unique solution.
Let us now derive a residual based error estimate for the semi-time discrete mixed form (15). We associate with (ỹ k ,w k ) the residuals R k 1 ∈ Y * and R k 2 ∈ W * by (16) and Next, we derive L 2 -representations of R k 1 and R k 2 by elementwise integration by parts The residual R k 1 fulfills the Galerkin orthogonality and it further holds true Further, the residual equation holds true for all (v 1 , v 2 ) ∈ Y × W : where the last step follows from (19). Now, we are in the position to derive a temporal residual based a-posteriori error estimate for the semi-time discrete mixed variational formulation (15).
Theorem 3.4. Let (ỹ,w) ∈ X denote the solution to (14) and let (ỹ k ,w k ) ∈ Y k × W denote the solution to (15). Then, the following residual based a-posteriori error estimate holds true: with a constant C > 0 and Proof. We combine (18) together with (20).
where we use the notation r k 1,int :=ỹ d + (ỹ k ) tt + ν∆w k − 1 αỹ k and r k 1,end := −(ỹ k ) t (T ) + ν∆ỹ k (T ). Note that the last summand vanishes at the interpolation point t = T . We can estimate using Cauchy-Schwarz Next, using standard interpolation properties (see e.g. [1, Theorem 1.7]), we arrive at whereĨ i denotes the set of intervals which share a vertex with I i . We recall that |.| H 1 denotes the H 1 -seminorm. Together with the Cauchy-Schwarz inequality for sums, we arrive at where we use Hölder's inequality in the last step. We note that In (23) we choose v 1 :=ỹ −ỹ k and v 2 :=w −w k and denote e := (ỹ −ỹ k ,w −w k ), which leads to By the definition of the energy norm ||| · |||, it follows A M (e, e) ≥ |||e||| 2 which yields the a-posteriori error estimate Remark 3.1 (Adaptive cycle). In order to construct an adaptive time grid, we follow the standard solve → estimate → mark → refine cycle. In practice, we solve (15) using space-time finite elements. Then, the error in each time interval is estimated using (21). The intervals with the largest errors are marked using the Dörfler marking strategy [9]. For refinement, we perform a bisection of the marked intervals. We iterate this loop until the time grid has a prescribed number of time instances.
Remark 3.2 (Heuristic assumption). Note that we derived an error estimate (21) for a time discrete formulation in y whereas w is kept continuous. In practice, we solve a fully space-time discrete mixed variational formulation, but still use the error estimate for the semi-time discrete form to construct an adaptive time grid. For this, we assume that the temporal discretization of y k is insensitive with respect to the spatial discretization. In fact, numerical studies in [3,4] show that temporal and spatial discretization decouple for the considered problem settings. In addition, we also assume that a temporal discretization of w k does not strongly influence the error estimate. Of course, these heuristic assumptions might not hold in general. For this reason, we will in future research derive a-posteriori error estimates for a fully space-time discrete mixed variational form.
With the help of (21), we are able to refine the time grid by means of the residual of the system (13). This property will constitute a building block in the Model Predictive Control (MPC) framework as discussed in the next Sections 4,5 and 6.

State equation with depletion term
Let us now consider an optimal control problem of the form (4), where an additional depletion term in the state equation (1) appears as with µ > 0. The reformulation of the associated optimality system into an elliptic equation and an associated mixed formulation, respectively, follows along the lines of Sections 3.1 and 3.2. In particular, the mixed formulation reads as in Ω, y(0) = 0 in Ω.
Let us define the bilinear form The semi-discrete mixed variational formulation then reads as With similar arguments as in the previous sections, one can show existence of a unique solution of the involved equations provided sufficient regularity of the data.
In analogy to Theorem 3.4 we can derive a temporal residual based a-posteriori error estimate for (27). Theorem 3.5. Let (ỹ,w) ∈ X denote the solution to (26) and let (ỹ k ,w k ) ∈ Y k × W denote the solution to (27). Further, let µ ≤ ν/c 2 p , where c p denotes the Poincaré constant. Then, the following residual based a-posteriori error estimate holds true: with a constant C > 0 and Proof. The proof follows along the lines of the proof of Theorem 3.4. Note that it holds Using Green's formula, the definition ofw and Young's inequality, we can estimate the second summand in (30) by With the choice δ := 1 + 2αµ 2 8αµ 2 , it holds that −4δµ 2 + 1 α + µ 2 ≥ 0 and − 1 4δ + 1 ≥ 0.

Model predictive control
In this section, we recall the MPC method to solve the optimal control problem (4). For a comprehensive study we refer to e.g. [13,23]. The basic idea of the MPC approach is to split the (large) time horizon [0, T ] of the optimal control problem (4) into a sequence of smaller time intervals and solve corresponding open loop control problems which allows to compute a state feedback law.
Let us first focus on a uniform time grid 0 = t 0 < t 1 < · · · < t m := T with fixed time step size ∆t = T /m. In Section 5, we will propose a time adaptive method. Let us define a time domain with small time horizon as [t i , t The MPC method works as follows. We start by solving the optimal control problem (32) in [t 0 , t N 0 ] × Ω and we store the optimal control u N on the first subinterval [t 0 , t 0 + ∆t] together with the associated optimal state trajectory. Then, we start a new finite horizon optimal control problem on the shifted time horizon [t 1 , t N 1 ] with t 1 = t 0 + ∆t where the initial condition y 1 is given by the optimal trajectory y [u N ,t0,y0] (t) at t = t 0 + ∆t using the optimal control u N (t) for t ∈ (t 0 , t 0 +∆t]. We store the feedback map φ N which is defined as φ N (y [u N ,t0,y0] (t)) := u N (t). We iterate this process. Note that (32) Compute the associated state y N = y [φ N ,ti,yi] (t) by solving (24) on (t i , t i + ∆t]. Set y i+1 = y N (t i + ∆t) and t i+1 = t i + ∆t. end for In general, it is well-known that the larger the prediction horizon, the better the feedback law one can obtain, see e.g. [13]. However, one is interested in short prediction horizons to minimize the computational cost of the method (or even horizon of minimal length) while guaranteeing certain properties of the MPC scheme such as e.g. stability of the method. The computation of this minimal horizon is related to a relaxed dynamic programmic principle (see [13,Chapter 6]).

Time adaptivity in MPC
In this section, we propose time adaptive techniques within MPC. In Algorithm 1 we have to choose a-priori the prediction horizon and the number of degrees of freedom in each temporal subinterval. Here, we would like to reply to the following question: How to choose a suitable time discretization for the prediction horizon in each level of the MPC?
We aim at computing the temporal discretization to identify the important dynamical structures according to the optimization goal. We will propose adaptive strategies which avoid very small uniform temporal discretizations and realize an efficient implementation. The proposed approaches will lead to adaptive time discretizations which are related to the optimal state. Therefore, we propose different approaches to combine the computation of an adaptive time grid using the a-posteriori estimate provided in Section 3.3 within a MPC framework. The idea of adaptivity leads to different combinations using the error estimate (28). Here, we will deal with (i) an adaptive length of the horizon where the time discretization is computed a-priori (offline approach) and with (ii) an adaptive grid in each subinterval for a fixed prediction horizon where the time discretization is computed on the fly (online approach). For a different adaptive concept based on goal-oriented adaptivity, see the recent work [16].
(1) Offline approach: Compute an adaptive time grid a-priori before starting the MPC We compute an adaptive time grid {t i } m+N −1 i=0 according to (28), which provides a suitable time grid tailored for the optimal state. Note that the final time is t m = T , whereas t m+N −1 is the last time instance necessary to compute the feedback control at time t m . Under the heuristic assumptions of Remark 3.2 we will use a very coarse spatial mesh to solve (27) adaptively in time for the whole time domain before we start the MPC procedure. The advantage of using a coarse spatial resolution is that the construction of an adaptive time grid is computationally very cheap. Moreover, the resulting adaptive time grid is related to the optimal state of the original control problem which in general is our target in MPC. This grid is then used as the time grid in a MPC framework as shown in Figure 1, where we use a chosen number of time instances N in the subintervals [t i , t i+N −1 ] for i = 0, . . . , m − 1. Note that this procedure allows prediction horizons with different lengths t i+N −1 − t i =T i . Note that in tracking MPC the desired state y d frequently stems from an open loop optimal control problem. This further justifies the offline approach, since the grid well represents y d . The approach is summarized below in Algorithm 2.   . In order to make computations even more efficient, the information of the previous MPC iteration can be used as a warm start for the next MPC iteration. In particular, after a coarsening step of the previous adaptive time grid, this grid can be used as an initial adaptive time grid for the next prediction horizon. Furthermore, to improve the inner open-loop solver in each iteration one can use as initial control the one computed at the previous step.

· · ·
Remark 5.2 (Comparison of the two approaches). The two proposed approaches have different features. The offline approach computes the temporal grid before the MPC procedure begins. For this reason, the error estimate (28) is used in an adaptive cycle for the original optimal control problem. The online approach, however, computes an adaptive time grid at each MPC iteration for, eventually, small temporal domains. This means that the error estimate is used in an adaptive cycle for each open loop MPC subproblem. Those approaches, in general, lead to different grids. The offline method provides an adaptive subinterval; the size changes according to the error estimator, whereas the online approach has a fixed size for the subinterval but it provides different time instances.

Numerical example
In the following numerical tests, we compare the different adaptive time schemes within MPC proposed in Section 5. We will investigate numerically the influence of the time grid on the approximation quality of the MPC methods. In particular, we compare Algorithms 1, 2 and 3.
In all numerical examples, the considered spatial domain is the open interval Ω = (0, 1) and the temporal domain is [0, T ] = [0, 1]. In order to solve the mixed form (25), we introduce a partitioning of the space-time domain into regular orthotopes and use Q 1 space-time finite elements for discretization, where Q 1 is the space of polynomials of separate degree up to 1. We solve the equation with a direct solver using a coarse spatial resolution of ∆x = 1/5. For the solution of the MPC open loop subproblems, we use an implicit Euler scheme for the temporal discretization and use P 1 finite elements for the spatial discretization. The optimal control problem is solved with a direct solver. We take as fine spatial resolution an equidistant discretization with ∆x = 1/100. All coding is done in Matlab R2019a.
For small values of ε (we use ε = 10 −3 ), the state y develops a very steep gradient at t = 0.5, which can be seen in Figure 3.  Figure 4. We can see that the standard MPC Algorithm 1 with equidistant time grid fails whereas using either of the Algorithms 2 or 3, it is possible to capture the layer at t = 0.5 and the solutions comply much better with the true state solution Let us now provide more details about the grids we obtained with the methods proposed. In the offline approach (Algorithm 2), the temporal grid is computed before starting the MPC algorithm. The adaptive grids using a coarse and a fine spatial resolution for solving (13) are shown, respectively, in the middle and right panel of Figure 5. As already discussed in [4,12] and recalled in Remark 3.2, we observe numerically that the temporal and spatial  Figure 5: Test 1: Uniform space-time grid with fine spatial resolution (left), offline adaptive grid with coarse (middle) and fine (right) spatial resolution resolution decouple and we are able to compute cheaply the temporal grid by means of a coarse spatial resolution for the solution of (13). Examples of the corresponding equidistant and adaptive temporal intervals are compared in the top and bottom panels of Figure 6. It is clear that the adaptive method refines where the solution exhibits a steep gradient and, therefore, the solution of the problem will be better approximated. Furthermore, from the top pictures in Figure 6 it is clear how the length of the prediction horizonT j is different for each iteration of the MPC method.
The online adaptive grid with a coarse and a fine spatial resolution is shown in the middle and right panel of Figure 7. Again, we observe that the time adaptivity is very insensitive with respect to the spatial resolution.
The corresponding adaptive time horizons are shown in the top panels of Figure 8. As a comparison, the uniform time horizons of the same lengths are shown in the bottom panels of Figure 8 using the same number of degrees of freedom in each interval. It is clear that the a-posteriori error estimate (21) leads to a time grid associated with the open loop optimal state which benefits the accuracy of the control problem.
Finally, we provide an error analysis for the computation of the approximate solutions using Algorithms 1, 2 and 3 for different choices of degrees of freedom in time and prediction horizons. For this, we compute the error between the analytical optimal state solution to (6) and its numerical approximation using the different MPC approaches measured in the L 2 (0, T ; Ω)−norm.
We would like to mention that it is not easy to make a completely fair comparison due to the differences in the strategies and the influence of the time grid (and prediction horizon) on the MPC approach.
In Figure 9 we show the error for Algorithms 1 and 2. We fixed an equidistant time grid in the first case and an adaptive one with the same number of degrees of freedom and we compare the results for different choices of N . Depending on whether the layer at t = 0.5 is a time discretization point or not, the approximation quality can differ strongly leading to the illustrated zig-zag behavior. Since the exact location of the layer is usually not known a-priorily, an equidistant time grid approach (Algorithm 1) is easy to fail. Moreover, we never reach the same level of accuracy of the adaptive grid unless we further increase the number of degrees of freedom in time. On the other hand, on the right panel of Figure 9 we In Figure 10 we compare the computational time (measured in seconds) for Algorithms 1 and 2. Note that the elapsed time for the offline adaptive MPC approach also includes the computational time to generate the adaptive grid. As one can see the times are comparable. Without any doubts we can conclude that the offline approach with an adaptive prediction horizon leads to faster and more accurate results than using an equidistant grid. To conclude our comments on the offline approach we show the results of the terms y − y d 2 L 2 (0,T ;Ω) in the top panels and u 2 L 2 (0,T ;Ω) in the bottom panels of Figure 11 comparing the MPC state and control solution for different horizon lengths and degrees of freedom. As a reference, we provide the values for the terms using the true state and control solution to the original optimal control problem evaluated on the same meshes ("True") and the analytical value which is y true − y d 2 L 2 (0,T ;Ω) = 26.8197 and u true 2 L 2 (0,T ;Ω) = 0.25 ("Analytical"). Moreover, Figure 12 compares the evolution of y(t) − y d 2 L 2 (Ω) and u(t) 2 L 2 (Ω) over time comparing the MPC state and control solution with the true solution evaluated at the given time grid.
In Figure 13 we compare the L 2 -error for Algorithms 1 and 3. We fixed the prediction horizonT and modified the choice of the instances in each sub interval using the equidistant and adaptive method. As one can see, with this approach we need a small prediction horizon and a large number of time instances to obtain an error of order 10 −1 with an equidistant grid whereas the adaptive method provides a more flexible approach for a small T = {0.1, 0.2, 0.3, 0.4}. Again, we can see a zig-zag phenomena in the case of an equidistant choice of the grid which further emphasizes the importance of an adaptive method.
In Figure 14 we compare the computational time in seconds of Algorithm 3 using an equidistant or adaptive time grid including the computational time needed to create the adaptive time discretization within each MPC iteration.
Clearly, to obtain a more accurate solution it is computationally more expensive but we also want to remark that the minimum error with the equidistant grid is 0.0872 computed in 25.87s whereas, with the adaptive approach, to get an error of 0.0216 we needed 16.06s. This shows that our method is more accurate and also more efficient computationally without

Test 2: State equation with depletion term
In this numerical test, we consider the optimal control of (24) with µ > 0. The regularization parameter in the cost is chosen as α = 1 and the desired state is given by The initial condition for the state is y 0 (x) ≡ 10 sin(πx)exp(−1/(4ε 2 )) and the source term in the state equation is chosen in order to get as exact optimal solution (y, u) for the original control problem: Let us note that the Poincaré constant c p and the first eigenvalue λ 1 of the Laplace-Dirichlet operator are related by λ 1 = 1/c 2 p (see e.g [5,Proposition 8.4.3]). For Ω = (0, 1), the first eigenvalue λ 1 of the Laplace-Dirichlet operator is given by λ 1 = π 2 (see e.g. [5,Proposition 8.5.2]). Then, since Theorem 3.5 is applicable if µ ≤ ν/c 2 p , for this setting it requires µ ≤ 0.1 · π 2 . In our numerical simulation we set ν = 0.1, µ = 3, ε = 10 −2 . Thus, we show the unstable case for µ = 3 which should not hold according to our theoretical findings. Nevertheless, the numerical results under this configuration provide accurate results, very similar to a stable case with µ ≤ νπ 2 as required in Theorem 3.5.
The optimal solution is shown on the left panel of Figure 6.2. As one can see we have a layer at t = 0.5. Then we show the results of Algorithm 2 choosing N = 5, m = 45 with an equidistant grid (middle panel of Figure 6.2) and an adaptive grid (right panel of Figure  6.2). It is easy to see that the adaptive grid catches better the behavior of the optimal solution. The grid for this setting is shown in Figure 6.2 (right) and, as expected, our adaptive approach refines consistently in a neighborhood of the layer. The quality of our results are also confirmed by the L 2 −error plots in Figure 17. Under the settings considered we can see that the error with an equidistant grid is of order, around, O(10 0 ) whereas with an adaptive grid we can reach O(10 −1 ). Clearly decreasing ∆x we can improve the accuracy of our approximation.

Conclusions and Outlook
In this work we have proposed two approaches to generate adaptive time discretization in the MPC framework. Our approach is fully flexible and relies on a reformulation of the optimal control problem into a second order in time and fourth order in space equation. Our approach does not require further assumptions on the control problem. The use of a-posteriori error estimates to generate the time grid in the MPC method is the important novelty of our work. Numerical tests have shown the efficiency of the method for both accuracy and computational time. We also want to remark that our approach is particularly suitable when a layer is shown in the solution. Other experiments with mild temporal variations did not always show a clear difference between equidistant and adaptive grid. We can say that the offline approach with an adaptive grid did provide slightly better results. The a-posteriori error indicator delivers an appropriate adaptive time grid even providing a coarse spatial resolution. In the future we plan to derive an a-posteriori error estimator for a fully space-time discrete form and to use that indicator for a fully adaptive and automatic MPC scheme, where the idea is to avoid an a-priori choice of the prediction horizon and/or the number of degrees of freedom in each sub-iteration. Another goal is to extend these results to nonlinear control problem and as soon as we increase the dimension of the problem make use of efficient model reduction techniques, such as POD, to decrease the computational time.