\section{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 \cite{HinPUU09} 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 \cite{BarC97} for the detailed mathematical analysis of the problem and to \cite{FalF13} for the numerical approximation and the references therein. For large-scale problems with a HJB approach we refer e.g. to \cite{AllFS19,KalK18}.

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 {\em moving horizon control} or {\em receding horizon control}. We refer to the monographs \cite {GruP17,RawMD09} 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 \cite{Che98, ItoK02} 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.\ \cite{GruSS19,Kre18} 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.\ \cite{Zas06}) is often a key tool to prove asymptotic stability of the MPC method and to find the minimal prediction horizon (see e.g. \cite{BreP20,GruPSetal09,GruSS20,KunP19}).

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 \textit{optimal} state solution. The idea is based on \cite{GonHZ12}, 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 \cite{AllGH16,AllGH18}). The use of an a-posteriori error estimator leads us to two different approaches.
\begin{enumerate}
\item 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.
\item 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.
 \end{enumerate}
 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 \cite{GSSnew} 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~\ref{sec1}, we present the optimal control problem together with the optimality conditions. In Section~\ref{sec2} 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~\ref{secmpc}, we recall the basic idea of the MPC method. In Section~\ref{secadaptivity}, we propose the novel time adaptive schemes in MPC which are an offline and an online approach. Finally, numerical tests are discussed in Section~\ref{secnumerics} and conclusions are made in Section~\ref{secconcl}.\\

\section{Preliminaries and optimal control setting}\label{sec1}

\subsection{Preliminaries}
Let $\Omega \subset \mathbb{R}^n, n \in \{1,2,3\}$ be an open and bounded domain with Lipschitz boundary $\partial \Omega$ and let $T>0$ be a given end time. The Lebesgue space of square integrable functions is denoted by $L^2(\Omega)$  with inner product $(u,v)_{L^2(\Omega)}:= \int_\Omega uv dx$ and norm $\|u\|_{L^2(\Omega)}:=(\int_\Omega |u(x)|^2 dx)^{1/2}$ for $u,v \in L^2(\Omega)$. Further, let $H^k(\Omega)$ defined by
$$H^k(\Omega) := \{u\in L^2(\Omega): u \text{ has }\text{weak }\text{derivatives } D^{\beta}u \in L^2(\Omega) \text{ for all } |\beta | \leq k\}$$
with $k \in \mathbb{N}_0$ and equipped with the norm $\|u\|_{H^{k}(\Omega)}  := (\sum_{|\beta| \leq k} \|D^\beta u \|_{L^2(\Omega)}^2)^{1/2}$  and
$$ H_0^k(\Omega):= \{ u \in H^k(\Omega): D^\beta u = 0 \text{ on } \partial \Omega \text{ in the sense of traces } (|\beta| \leq k-1)\}. $$
\noindent We use the notation $H^{-1}(\Omega)$ for the dual space of $H_0^1(\Omega)$ and denote $\langle \cdot , \cdot \rangle_{H^{-1}(\Omega),H_0^1(\Omega)}$ as the duality pairing of $H^{-1}(\Omega)$ with $H_0^1(\Omega)$. By $| \cdot |_{H^1(\Omega)}$ we denote the $H^1$-seminorm given by $|u|_{H^1(\Omega)} = \|\nabla u \|_{L^2(\Omega)}$ for $u\in H_0^1(\Omega)$. We recall that the Poincar\'{e} constant is given by the smallest number $c_p>0$ such that the Poincar\'{e} inequality
$$\|u\|_{L^2(\Omega)} \leq c_p \|\nabla u \|_{L^2(\Omega)}, \; \forall u \in H_0^1(\Omega)$$
is fulfilled. Thus, $|.|_{H^1(\Omega)}$ is a norm on $H_0^1(\Omega)$ equivalent to the norm $\|.\|_{H^1(\Omega)}$. 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)} := (\int_{0 }^T \|u(t)\|_X^2 dt )^{1/2}$. We define $W(0 ,T;H_0^1(\Omega)):=\{v\in L^2(0 ,T;H_0^1(\Omega)), v_t \in L^2(0 ,T;H^{-1}(\Omega))\}$. Note that for a given function $g$ in space-time, we use the short hand notation $g(t)$ to indicate the time dependency and drop the space argument.

\subsection{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
\begin{equation}\label{heat}
\left\{
\begin{array}{rcll}
y_t-\nu\Delta y & = & f+u &\text{ in } (0 ,T]\times\Omega,\\
y &= & 0 &\text{ on } [0 ,T] \times \partial \Omega,\\
y(0 ) & = & y_0 & \text{ in } \Omega,
\end{array}
\right.
\end{equation}
where $\nu > 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.\\
\noindent The weak form of \eqref{heat} reads as follows: for given $f \in L^2(0,T;\Omega)$, $y_0 \in L^2(\Omega)$ and $u\in L^2(0,T;\Omega)$, find a state $y \in W(0 ,T;H_0^1(\Omega))$ satisfying $y(0 ) = y_0$ such that
\begin{equation}\label{heat_weak}
 \langle y_t(t),v \rangle_{H^{-1}(\Omega),H_0^1(\Omega)} + \nu   \int_\Omega \nabla y(t) \cdot \nabla v dx  = \int_\Omega ( f(t)+u(t))v dx
\end{equation}
holds for all $v \in H_0^1(\Omega)$ and almost everywhere in $(0 ,T]$. For existence of a unique solution to \eqref{heat_weak}, we refer to e.g.\ \cite[\textsection 7.1.2, Theorems 3 and 4]{Eva10}.\\

\noindent We consider the minimization of the following cost functional
\begin{equation}\label{cost}
J(y,u):=\int_{0 }^T \ell(y(t),u(t)) dt,
\end{equation}
where $\ell$ has a tracking-type structure of the form
$$\ell(y(t),u(t)):= \frac{1}{2} \|y(t) - y_d(t)\|_{L^2(\Omega)}^2 + \frac{\alpha}{2}\|u(t)\|_{L^2(\Omega)}^2$$
with given desired state $y_d \in L^2(0 ,T;\Omega)$ and given regularization parameter $\alpha > 0$. Then, the optimal control problem can be formulated as
\begin{equation}\label{ocp}
\min_{(y,u)\in W(0 ,T;H_0^1(\Omega))\times L^2(0 ,T;\Omega)} J(y,u) \quad \text{s.t.} \quad y \text{ satisfies } \eqref{heat_weak} \text{ for } u \in L^2(0,T;\Omega).
\end{equation}
Under these settings the optimal control problem \eqref{ocp} admits a unique solution $u\in L^2(0 ,T;\Omega)$ with associated state $y \in W(0 ,T;H_0^1(\Omega))$, see e.g.\ \cite{HinPUU09}.\\

\noindent We indicate through the index of the weak solution $y_{[u,t_0,y_0]}$ to \eqref{heat} its dependency on a given control $u$, initial time $t_0$ and initial value $y_0$. This allows to introduce a reduced cost functional $\hat{J}$ defined by
\begin{equation}\label{reducedcost}
 \hat{J}(u,t_0,y_0) := \int_{t_0}^T \ell(y_{[u,t_0,y_0]}(t),u(t)) dt.
\end{equation}
Thus, the optimal control problem \eqref{ocp} can be formulated as
\begin{equation}\label{ocpreduced}
\min_{u \in L^2(0,T;\Omega)} \hat{J}(u,t_0,y_0).
\end{equation}
with $t_0=0$. Following the standard Lagrange technique, the optimality system associated with the optimal control problem \eqref{ocpreduced} is given by the following equations understood in a weak sense: the state equation \eqref{heat} together with the adjoint equation
\begin{equation}\label{adj}
\left\{
\begin{array}{rcll}
-p_t-\nu\Delta p  & = & y-y_d &\text{ in } [0 ,T)\times\Omega,\\
p & = & 0 &\text{ on } [0 ,T] \times \partial \Omega,\\
p(T) & = & 0 &\text{ in } \Omega,
\end{array}
\right.
\end{equation}
and the optimality condition
\begin{equation}\label{opt_con}
\alpha u + p=0\quad \mbox{ in } [0 ,T]\times\Omega.
\end{equation}

For the next chapter, we will need the following higher regularity results for the weak solution of \eqref{heat} and \eqref{adj}, respectively.

 \begin{lmm}[Higher regularity \cite{Eva10}]\label{lemma:regularity}
(i) Let $y_0 \in H_0^1(\Omega)$ and let $f$, $u$, $y_d \in L^2(0 ,T;\Omega)$. Then, according to \cite[\textsection 7.1.3. Theorem 5]{Eva10} the weak solution $y$ of \eqref{heat} and the weak solution $p$ of \eqref{adj} fulfill $y,p \in L^2(0 ,T;H^2(\Omega)) \cap L^\infty(0 ,T;H_0^1(\Omega))\cap  H^1(0 ,T;L^2(\Omega))$.\\
 (ii) Let $y_0 \in H_0^1(\Omega) \cap H^3(\Omega)$ and $f,u, y_d \in L^2(0 ,T;H^2(\Omega))\cap H^1(0 ,T;L^2(\Omega))$. Further, let the compatibility assumption $(u+f)(0)+\nu \Delta y_0  \in H_0^1(\Omega)$ hold true. Then according to \cite[\textsection 7.1.3. Theorem 6]{Eva10} the weak solution $y$ of \eqref{heat} and the weak solution $p$ of \eqref{adj} fulfill $y,p \in L^2(0 ,T;H^4(\Omega))\cap H^1(0 ,T;H^2(\Omega))\cap H^2(0 ,T;L^2(\Omega))$.
 \end{lmm}

\section{Reformulation of the optimality system and time adaptivity}\label{sec2}
\subsection{Reformulation of the optimality system}\label{sec2:subreform}
Following along the lines of \cite{GonHZ12}, we can reformulate the optimality system \eqref{heat}-\eqref{adj}-\eqref{opt_con} 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
\begin{equation}\label{2ord}
\left\{
\begin{array}{rcll}
-y_{tt}+\nu^2\Delta^2 y  +\frac{1}{\alpha} y & = &\frac{1}{\alpha}y_d -f_t-\nu \Delta f & \text{ in } (0 ,T)\times\Omega,\\
y & = & 0 &\text{ on } [0 ,T]\times\partial\Omega,\\
\nu \Delta y & = & -f &\text{ on } [0 ,T]\times\partial\Omega,\\
\left(y_t-\nu \Delta y \right)(T) & = & f(T) &\text{ in }\Omega,\\
y(0 ) & = & y_0 &\text{ in } \Omega.
\end{array}
\right.
\end{equation}
We note that for $\nu=1$ and $f\equiv0$ this setting coincides with the setting considered in \cite{GonHZ12}. Under higher regularity assumptions on the data, the following theorem shows that the optimal state $y$ of \eqref{heat}-\eqref{adj}-\eqref{opt_con} fulfills the elliptic equation \eqref{2ord} a.e.\ in space-time.

\begin{thrm}
 Let $(y,u) \in W(0 ,T;H_0^1(\Omega))\times L^2(0 ,T;\Omega)$ with associated adjoint $p\in W(0 ,T;H_0^1(\Omega))$ denote the unique weak solution to \eqref{heat}-\eqref{adj}-\eqref{opt_con}. Further, let the assumptions of Lemma~\ref{lemma:regularity}(ii) be fulfilled. Then, $y$ satisfies \eqref{2ord} a.e.\ in space-time.
\end{thrm}
\begin{proof}
The proof follows along the lines of the proof of \cite[Theorem 2.7]{GonHZ12} and uses differentiation and insertion of the equations \eqref{heat}-\eqref{adj}-\eqref{opt_con}.
\end{proof}

\noindent Let us homogenize \eqref{2ord}. For this, let $g$ be a function which fulfills the boundary conditions as well as initial and end time conditions of \eqref{2ord} and is sufficiently smooth. Let $y$ satisfy \eqref{2ord}. We define $\tilde{y}:=y-g$ and arrive at
\begin{equation}\label{2ord_hom}
\left\{
\begin{array}{rcll}
-\tilde{y}_{tt}+\nu^2\Delta^2 \tilde{y} + \frac{1}{\alpha} \tilde{y} & = &\tilde{y}_d & \text{ in } (0 ,T)\times\Omega,\\
\tilde{y} & = & 0 &\text{ on } [0 ,T]\times\partial\Omega,\\
\nu \Delta \tilde{y} & = & 0 &\text{ on } [0 ,T]\times\partial\Omega,\\
\left(\tilde{y}_t-\nu \Delta \tilde{y} \right)(T) & = & 0 &\text{ in }\Omega,\\
\tilde{y}(0 ) & = & 0 &\text{ in } \Omega,
\end{array}
\right.
\end{equation}
where
\begin{equation}\label{yd}
\tilde{y}_d:= \frac{1}{\alpha} y_d - f_t - \nu \Delta f  + g_{tt} - \nu^2 \Delta^2 g - \frac{1}{\alpha} g.
\end{equation}

\noindent Now, let us derive a weak formulation of \eqref{2ord_hom}. For this purpose we introduce the function space
$$H^{2,1}_0(0 ,T;\Omega):=\left\{v\in H^{2,1}(0 ,T;\Omega): v(0 )=0 \mbox { in }\Omega\right\},$$
where
$$H^{2,1}(0 ,T;\Omega):=L^2(0 ,T;H^2(\Omega)\cap H^1_0(\Omega))\cap H^1(0 ,T; L^2(\Omega)).$$ It is equipped with the norm $$\|v\|_{H^{2,1}(0 ,T;\Omega)}:= \left( \|v\|_{L^2(0 ,T;H^2(\Omega))}^2 + \|v\|_{H^1(0 ,T;L^2(\Omega))}^2 \right)^{1/2}.$$ We introduce the following symmetric bilinear form
$$A:H^{2,1}_0(0 ,T;\Omega)\times H^{2,1}_0(0 ,T;\Omega)\rightarrow\mathbb R,$$
$$A(v_1,v_2):=\displaystyle\int_{0 }^T \int_\Omega \left( (v_1)_t (v_2)_t + \nu^2 \Delta v_1 \Delta v_2  + \frac{1}{\alpha}  v_1 v_2 \right)dxdt + \displaystyle\int_\Omega   \nu\nabla v_1(T) \nabla v_2(T) dx,$$\\
and linear form
$$L:H^{2,1}_0(0 ,T;\Omega)\rightarrow\mathbb R, \quad L(v) := \int_0^T \int_\Omega \tilde{y}_d v \; dxdt$$
where $\tilde{y}_d$ is defined in \eqref{yd}.

\begin{dfntn}(Weak formulation)
 The weak formulation of equation \eqref{2ord_hom} is given by: find $\tilde{y} \in H_0^{2,1}(0 ,T;\Omega)$ which satisfies
\begin{equation}\label{weak_2ord}
A(\tilde{y},v)=L(v)\quad \forall v\in H^{2,1}_0(0 ,T;\Omega).
\end{equation}
\end{dfntn}

Existence of a solution to \eqref{weak_2ord} and its relation to a solution to \eqref{2ord} is shown in the following theorem.

\begin{thrm}\label{thm:existence}
 Let $y$ denote a solution to \eqref{2ord} and let $g$ be a function which fulfills the boundary, initial and end time conditions in \eqref{2ord} and is sufficiently smooth. Then, $\tilde{y}=y-g$ is a solution to \eqref{weak_2ord}. On the other hand, if $\tilde{y}$ is a solution to \eqref{weak_2ord} and the assumptions of Lemma~\ref{lemma:regularity}(ii) are fulfilled, then $y=\tilde{y}+g$ satisfies \eqref{2ord} a.e.\ in space-time.
\end{thrm}
\begin{proof}
Assume $y$ is a solution to \eqref{2ord}. By Green's formula and integration by parts it is straight forward to prove that $\tilde{y}=y-g$ satisfies \eqref{weak_2ord}. The other direction follows vice versa.
\end{proof}

In order to show equivalence of the optimal control problem \eqref{ocp} to the weak formulation of \eqref{2ord} it remains to prove uniqueness of a solution.

\begin{thrm}\label{thm:unique}
 The solution $y$ to \eqref{weak_2ord} is unique.
\end{thrm}
\begin{proof}

The proof follows along the lines of the proof of \cite[Theorem 2.6]{GonHZ12} and uses Lax-Milgram Lemma (see e.g. \cite[\textsection 6.2.1, Theorem 1]{Eva10}).
\end{proof}

\subsection{Mixed formulation}\label{secmixedform}
 In order to use piecewise linear, continuous finite elements for discretization and avoid the construction of finite element subspaces in $H^2(\Omega)$, we introduce an auxiliary variable $\tilde{w}:=-\nu \Delta \tilde{y}$. This allows to write \eqref{2ord_hom} as a coupled system in $\tilde{y}$ and $\tilde{w}$ as
\begin{equation}\label{2ord_decoupled}
 \left\{
 \begin{array}{rcll}
 -\tilde{y}_{tt}-\nu \Delta \tilde{w}  +\frac{1}{\alpha} \tilde{y} & = &\tilde{y}_d & \text{ in } (0 ,T)\times\Omega,\\
 \nu \Delta \tilde{y} + \tilde{w} & = & 0 &  \text{ in } (0 ,T)\times\Omega,\\
 \tilde{y} & = & 0 &\text{ on } [0 ,T]\times\partial\Omega,\\
  \tilde{w} & = & 0 &\text{ on } [0 ,T]\times\partial\Omega,\\
 \left(\tilde{y}_t-\nu \Delta \tilde{y} \right)(T) & = & 0 &\text{ in }\Omega,\\
 \tilde{y}(0 ) & = & 0 &\text{ in } \Omega.
 \end{array}
 \right.
 \end{equation}
 We introduce the function spaces $Y:=\{v \in H^1(0 ,T;H_0^1(\Omega)): v(0 ) = 0 \text{ in } \Omega \}$ and $W:=L^2(0 ,T;H_0^1(\Omega))$ and the product space $X:=Y\times W$. Let us define the following bilinear form
 $$ A_M:X \times X \to \mathbb{R},$$
\begin{align*}
A_M((\tilde{y},\tilde{w}),(v_1,v_2)) &= \displaystyle\int_{0 }^T \int_\Omega
 \tilde{y}_{t}(v_1)_t+\nu \nabla \tilde{w} \nabla v_1  + \frac{1}{\alpha} \tilde{y} v_1  -\nu \nabla \tilde{y} \nabla v_2 + \tilde{w} v_2 dxdt  \\
 & \qquad + \displaystyle\int_\Omega  \nu \nabla \tilde{y}(T)\nabla v_1(T)  dx
\end{align*}
and linear form
$$ L_M: X \to \mathbb{R}, \quad L_M(v_1,v_2) = \displaystyle\int_{0 }^T \int_\Omega \tilde{y}_d v_1 dxdt.$$

\begin{dfntn}
 The weak formulation of the mixed formulation \eqref{2ord_decoupled} is given by: find $(\tilde{y},\tilde{w}) \in X$ which satisfies
 \begin{equation}\label{decoupled_weak}
  A_M((\tilde{y},\tilde{w}),(v_1,v_2)) = L_M(v_1,v_2) \quad \forall (v_1,v_2) \in X.
 \end{equation}

 \end{dfntn}

 By analogy with Theorem~\ref{thm:existence} and Theorem~\ref{thm:unique} it can be shown that the mixed variational form \eqref{decoupled_weak} admits at most one solution and that the pair $(\tilde{y},\tilde{w})$ with $\tilde{y}$ denoting the unique solution to \eqref{2ord_hom} and $\tilde{w}:=-\nu \Delta \tilde{y}$ is a solution to the mixed variational form \eqref{decoupled_weak}. This means that the unique solution to \eqref{2ord_hom} defines the solution to the mixed variational form \eqref{decoupled_weak}.\\

Note that it is
 $$ A_M((y,w),(y,w)) =  \displaystyle\int_{0 }^T \int_\Omega
 y_t^2   + \frac{1}{\alpha}y^2  + w^2 dxdt + \displaystyle\int_\Omega  \nu |\nabla y(T) |^2  dx.
$$
For this reason, we define an energy norm associated with the bilinear form $A_M$ by
  $$ |||(y,w)|||:= \left( \int_{0 }^T \int_\Omega y_t^2 + \frac{1}{\alpha}  y^2 + w^2 dxdt \right)^{1/2}.$$

\subsection{A-posteriori error estimate for the semi-time discrete mixed form}\label{sec2:timedisc}

Let us now consider a semi-time discretization of \eqref{decoupled_weak} with respect to $\tilde{y}$ while the variable $\tilde{w}$ is kept continuous. We introduce a time grid $0=t_0 < t_1 < \dots < t_m = T$  with $m \in \mathbb{N}$, time step sizes $\Delta t_i = t_i - t_{i-1}$ and time intervals $I_i = (t_{i-1},t_i]$ for $i=1,\dots,m$. The time discrete space $V^k$ is defined by
$$ V^{k} = \{v \in C^0(0 ,T;H^1(\Omega)): v|_{I_i} \in \mathbb{P}_p(I_i)\},$$
where $\mathbb{P}_p$ denotes the space of polynomials of degree up to $p$. We set $Y^{k}:=V^k \cap Y$.

\begin{dfntn}(Semi-time discrete mixed form)
The semi-time discrete mixed variational form reads as: find $(\tilde{y}^k,\tilde{w}^k) \in Y^k \times W$ such that
 \begin{equation}\label{weak_dis}
A_M((\tilde{y}^k,\tilde{w}^k),(v_1,v_2)) = L_M(v_1,v_2) \quad \forall (v_1,v_2) \in Y^k\times W.
 \end{equation}
 \end{dfntn}

  With arguments similar to those used for \eqref{decoupled_weak} we may show that problem \eqref{weak_dis} admits a unique solution.

 Let us now derive a residual based error estimate for the semi-time discrete mixed form \eqref{weak_dis}. We associate with $(\tilde{y}^k,\tilde{w}^k)$ the residuals $R_1^k \in Y^*$ and $R_2^k \in W^*$ by
\begin{equation}\label{r1k}
 R_1^k(v_1) = \displaystyle\int_{0 }^T \int_\Omega \tilde{y}_d v_1  -
 (\tilde{y}^k)_t (v_1)_t - \nu \nabla \tilde{w}^k \nabla v_1  -\frac{1}{\alpha} \tilde{y}^k v_1 dxdt
 - \displaystyle\int_\Omega  \nu \nabla \tilde{y}^k(T)\nabla v_1(T)   dx
\end{equation}
and
\begin{equation}\label{r2k}
 R_2^k(v_2) = \displaystyle\int_{0 }^T \int_\Omega  \nu \nabla \tilde{y}^k \nabla v_2 - \tilde{w}^k v_2  dxdt.
\end{equation}
Next, we derive $L^2$-representations of $R_1^k$ and $R_2^k$ by elementwise integration by parts
\begin{align*}
  R_1^k(v_1) &= \displaystyle\sum_{i=1}^m \int_{I_i} \int_\Omega \left\{ \tilde{y}_d + (\tilde{y}^k)_{tt} + \nu \Delta \tilde{w}^k   - \frac{1}{\alpha} \tilde{y}^k \right\}  v_1  dxdt\\
  & \qquad + \int_\Omega \left\{ - (\tilde{y}^k)_t(T)+ \nu \Delta \tilde{y}^k(T)   \right\} v_1(T)  dx
 \end{align*}
and
$$R_2^k(v_2) = \sum_{i=1}^m \int_{I_i} \int_\Omega \left\{ -\nu \Delta \tilde{y}^k - \tilde{w}^k \right\} v_2 dxdt. $$
\noindent The residual $R_1^k$ fulfills the Galerkin orthogonality
\begin{equation}\label{R1k_zero}
R_1^k(v_1) = 0 \quad \forall v_1 \in Y^k.
\end{equation}
and it further holds true
\begin{equation}\label{R2k_zero}
 R_2^k(v_2) = 0 \quad \forall v_2 \in W.
\end{equation}

\noindent Moreover, for $(\tilde{y},\tilde{w})\in Y \times W$ and $(\tilde{y}^k,\tilde{w}^k)\in Y^k \times W$ it holds for all $(v_1,v_2) \in Y^k \times W$:
\begin{equation*}
 \begin{array}{r c l}
  A_M((\tilde{y}-\tilde{y}^k,\tilde{w}-\tilde{w}^k),(v_1,v_2)) & = & A_M((\tilde{y},\tilde{w}),(v_1,v_2)) - A_M((\tilde{y}^k,\tilde{w}^k),(v_1,v_2))\\[0.5em]
  & = & L_M(v_1,v_2) - A_M((\tilde{y}^k,\tilde{w}^k),(v_1,v_2)) \quad = 0.\\[0.5em]
 \end{array}
\end{equation*}
Further, the residual equation holds true for all $(v_1,v_2) \in Y \times W$:
\begin{equation}\label{reseq}
 \begin{array}{r c l}
  A_M((\tilde{y}-\tilde{y}^k,\tilde{w}-\tilde{w}^k),(v_1,v_2))
  & = & R_1^k(v_1) + R_2^k(v_2)  \quad = R_1^k(v_1),\\[0.5em]
 \end{array}
\end{equation}
where the last step follows from \eqref{R2k_zero}. Now, we are in the position to derive a temporal residual based a-posteriori error estimate for the semi-time discrete mixed variational formulation \eqref{weak_dis}.

\begin{thrm}\label{thm:errest}
Let $(\tilde{y},\tilde{w}) \in X$ denote the solution to \eqref{decoupled_weak} and let $(\tilde{y}^k,\tilde{w}^k) \in Y^k\times W$ denote the solution to \eqref{weak_dis}. Then, the following residual based a-posteriori error estimate holds true:
\begin{equation}\label{est-thm31}
 ||| (\tilde{y}-\tilde{y}^k,\tilde{w}-\tilde{w}^k)|||^2 \leq C \eta^2,
\end{equation}
with a constant $C>0$ and
\begin{equation}\label{res:2ord}
 \eta^2 = \sum_{i=1}^m \int_{I_i} \int_\Omega (\Delta t_i)^2 \left| \tilde{y}_d + (\tilde{y}^k)_{tt} + \nu \Delta \tilde{w}^k -  \frac{1}{\alpha} \tilde{y}^k \right|^2 dxdt.
\end{equation}
\end{thrm}

\textit{Proof.}  We combine \eqref{R1k_zero} together with \eqref{reseq}. Let for $v_1 \in Y$ be $I_Y^k v_1$ the approximation to $v_1$ from $Y^k$. Then, it is
\begin{equation*}
 \begin{array}{r c l}
  A_M((\tilde{y}-\tilde{y}^k,\tilde{w}-\tilde{w}^k),(v_1,v_2))  =  R_1^k(v_1 - I_Y^k v_1) =\\
   \displaystyle\sum_{i=1}^m \int_{I_i} \int_\Omega r_{1,int}^k(v_1 - I_Y^k v_1) dxdt +  \int_\Omega r_{1,end}^k (v_1 - I_Y^k v_1)(T) dx,
 \end{array}
\end{equation*}
where we use the notation $r_{1,int}^k:= \tilde{y}_d + (\tilde{y}^k)_{tt} + \nu \Delta \tilde{w}^k- \frac{1}{\alpha} \tilde{y}^k$ and $r_{1,end}^k := -(\tilde{y}^k)_t(T) + \nu \Delta \tilde{y}^k(T) $. Note that the last summand vanishes at the interpolation point $t=T$. We can estimate using Cauchy-Schwarz
$$ |A_M((\tilde{y}-\tilde{y}^k,\tilde{w}-\tilde{w}^k),(v_1,v_2))| \leq \int_\Omega \left( \sum_{i=1}^m \| r_{1,int}^k\|_{L^2(I_i)} \| v_1 - I_Y^k v_1 \|_{L^2(I_i)} \right) dx.$$
Next, using standard interpolation properties (see e.g.\ \cite[Theorem 1.7]{AinO00}), we arrive at
\begin{equation*}
 |A_M((\tilde{y}-\tilde{y}^k,\tilde{w}-\tilde{w}^k),(v_1,v_2))| \leq \int_\Omega \left( \sum_{i=1}^m \| r_{1,int}^k \|_{L^2(I_i)} \; c_1 \; \Delta t_i | v_1 |_{H^1(\tilde{I}_i)} \right) dx,
\end{equation*}
where $\tilde{I}_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
\begin{equation}\label{est:AM}
\begin{array}{l l l}
  |A_M((\tilde{y}-\tilde{y}^k,\tilde{w}-\tilde{w}^k),(v_1,v_2))| \\
 \qquad\qquad\qquad\qquad \leq  c_1 \displaystyle\int_\Omega \left( \sum_{i=1}^m \| r_{1,int}^k \|_{L^2(I_i)}^2 (\Delta t_i)^2 \right)^{1/2} \left( \sum_{i=1}^m | v_1 |_{H^1(\tilde{I}_i)}^2 \right)^{1/2} dx   \\
\qquad\qquad\qquad\qquad   \leq  c_2 \displaystyle\int_\Omega \left( \sum_{i=1}^m \| r_{1,int}^k \|_{L^2(I_i)}^2 (\Delta t_i)^2 \right)^{1/2}  | v_1 |_{H^1(0,T)} dx \\
\qquad\qquad\qquad\qquad   \leq  c_2  \left(\displaystyle\int_\Omega \sum_{i=1}^m \| r_{1,int}^k \|_{L^2(I_i)}^2 (\Delta t_i)^2 dx \right)^{1/2} \left( \displaystyle\int_\Omega | v_1 |_{H^1(0,T)}^2 dx \right)^{1/2},
\end{array}
\end{equation}
where we use H\"older's inequality in the last step. We note that
$$\left( \displaystyle\int_\Omega | v_1 |_{H^1(0,T)}^2 dx \right)^{1/2} \leq \left( \int_0^T \int_\Omega (v_1)_t^2 + \frac{1}{\alpha} v_1^2 + v_2^2 dxdt \right)^{1/2} = ||| (v_1,v_2)|||.$$
In \eqref{est:AM} we choose $v_1:= \tilde{y}-\tilde{y}^k$ and $v_2 := \tilde{w}-\tilde{w}^k$ and denote $e:=(\tilde{y}-\tilde{y}^k,\tilde{w}-\tilde{w}^k)$, which leads to
$$ |A_M(e,e)| \leq c_2 \left( \int_\Omega \sum_{i=1}^m \|r_{1,int}^k \|_{L^2(I_i)}^2 (\Delta t_i)^2 dx \right)^{1/2} \cdot ||| e |||.$$
By the definition of the energy norm $||| \cdot |||$, it follows $A_M(e,e) \geq |||e|||^2$ which yields the a-posteriori error estimate
$$ |||e |||^2 \leq  C \left( \int_\Omega \sum_{i=1}^m \|r_{1,int}^k \|_{L^2(I_i)}^2 (\Delta t_i)^2 dx \right). \hspace{1cm}\square$$

\begin{rmk}[Adaptive cycle]
 In order to construct an adaptive time grid, we follow the standard
 \begin{center}
solve $\to$ estimate $\to$ mark $\to$ refine
 \end{center}
 cycle. In practice, we solve \eqref{weak_dis} using space-time finite elements. Then, the error in each time interval is estimated using \eqref{est-thm31}. The intervals with the largest errors are marked using the D\"orfler marking strategy \cite{Doe96}. 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.
\end{rmk}

\begin{rmk}[Heuristic assumption]\label{rem:heuristic}
Note that we derived an error estimate \eqref{est-thm31} 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 \cite{AllGH16,AllGH18} 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.
\end{rmk}

With the help of \eqref{est-thm31}, we are able to refine the time grid by means of the residual of the system \eqref{2ord_decoupled}. This property will constitute a building block in the Model Predictive Control (MPC) framework as discussed in the next Sections~\ref{secmpc},\ref{secadaptivity}~and~\ref{secnumerics}.

\subsection{State equation with depletion term}\label{secdepletion}
Let us now consider an optimal control problem of the form \eqref{ocp}, where an additional depletion term in the state equation \eqref{heat} appears as
\begin{equation}\label{heatdepletion}
\left\{
\begin{array}{rcll}
y_t-\nu\Delta y - \mu y& = & f+u &\text{ in } (0 ,T]\times\Omega,\\
y &= & 0 &\text{ on } [0 ,T] \times \partial \Omega,\\
y(0 ) & = & y_0 & \text{ in } \Omega,
\end{array}
\right.
\end{equation}
with $\mu > 0$.
The reformulation of the associated optimality system into an elliptic equation and an associated mixed formulation, respectively, follows along the lines of Sections~\ref{sec2:subreform} and \ref{secmixedform}. In particular, the mixed formulation reads as
\begin{equation}\label{2ord_decoupled_depletion}
 \left\{
 \begin{array}{rcll}
 -\tilde{y}_{tt}-\nu \Delta \tilde{w} + 2 \nu \mu \Delta \tilde{y} +\left(\frac{1}{\alpha}+\mu^2\right)\tilde{y} & = &\tilde{y}_d & \text{ in } (0 ,T)\times\Omega,\\
 \nu \Delta \tilde{y} + \tilde{w} & = & 0 &  \text{ in } (0 ,T)\times\Omega,\\
 \tilde{y} & = & 0 &\text{ on } [0 ,T]\times\partial\Omega,\\
  \tilde{w} & = & 0 &\text{ on } [0 ,T]\times\partial\Omega,\\
 \left(\tilde{y}_t-\nu \Delta \tilde{y} - \mu \tilde{y}\right)(T) & = & 0 &\text{ in }\Omega,\\
 \tilde{y}(0 ) & = & 0 &\text{ in } \Omega.
 \end{array}
 \right.
 \end{equation}
Let us define the bilinear form
 $$ A_M^\mu:X \times X \to \mathbb{R},$$
$$A_M^\mu((\tilde{y},\tilde{w}),(v_1,v_2)) = \displaystyle\int_{0 }^T \int_\Omega
\left( \tilde{y}_{t}(v_1)_t+\nu \nabla \tilde{w} \nabla v_1 - 2 \nu \mu \nabla \tilde{y} \nabla v_1 +\left(\frac{1}{\alpha}+\mu^2\right)\tilde{y} v_1 \right. $$
 $$\left. -\nu \nabla \tilde{y} \nabla v_2 + \tilde{w} v_2 dxdt\right)  + \displaystyle\int_\Omega  \nu \nabla \tilde{y}(T)\nabla v_1(T) - \mu \tilde{y}(T) v_1 (T) dx
$$
and linear form
$$ L_M^\mu: X \to \mathbb{R}, \quad L_M^\mu(v_1,v_2) = \displaystyle\int_{0 }^T \int_\Omega \tilde{y}_d v_1 dxdt$$
where $\tilde{y}_d:= \frac{1}{\alpha} y_d - f_t - \nu \Delta f - \mu f + g_{tt} - \nu^2 \Delta^2 g - 2 \nu \mu \Delta g - (\frac{1}{\alpha} + \mu^2) g$.
\begin{dfntn}
 The weak formulation of the mixed formulation \eqref{2ord_decoupled_depletion} is given by: find $(\tilde{y},\tilde{w}) \in X$ which satisfies
 \begin{equation}\label{decoupled_weakdepletion}
  A_M^\mu((\tilde{y},\tilde{w}),(v_1,v_2)) = L_M^\mu(v_1,v_2) \quad \forall (v_1,v_2) \in X.
 \end{equation}

 \end{dfntn}

The semi-discrete mixed variational formulation then reads as
\begin{equation}\label{semitimediscrete_depletion}
 A_M^\mu((\tilde{y}^k,\tilde{w}^k),(v_1,v_2)) = L_M^\mu(v_1,v_2) \quad \forall (v_1,v_2) \in Y^k \times W.
\end{equation}

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~\ref{thm:errest} we can derive a temporal residual based a-posteriori error estimate for \eqref{semitimediscrete_depletion}.

\begin{thrm}\label{thm:errest_depl}
Let $(\tilde{y},\tilde{w}) \in X$ denote the solution to \eqref{decoupled_weakdepletion} and let $(\tilde{y}^k,\tilde{w}^k) \in Y^k\times W$ denote the solution to \eqref{semitimediscrete_depletion}. Further, let $\mu \leq \nu/c_p^2$, where $c_p$ denotes the Poincar\'{e} constant. Then, the following residual based a-posteriori error estimate holds true:
\begin{equation}\label{est-thm31mu}
 ||| (\tilde{y}-\tilde{y}^k,\tilde{w}-\tilde{w}^k)|||^2 \leq C \eta^2,
\end{equation}
with a constant $C>0$ and
\begin{equation}\label{res:2ordmu}
 \eta^2 = \sum_{i=1}^m \int_{I_i} \int_\Omega (\Delta t_i)^2 \left| \tilde{y}_d + (\tilde{y}^k)_{tt} + \nu \Delta \tilde{w}^k -2 \nu \mu \Delta \tilde{y}^k -  \left(\frac{1}{\alpha} +\mu^2 \right)\tilde{y}^k \right|^2 dxdt.
\end{equation}
\end{thrm}
\begin{proof}
The proof follows along the lines of the proof of Theorem~\ref{thm:errest}. Note that it holds
\begin{align}\label{eqthm}
 A_M^\mu((\tilde{y},\tilde{w}),(\tilde{y},\tilde{w})) &= \int_0^T \int_\Omega \tilde{y}_t^2 -2 \nu \mu |\nabla \tilde{y}|^2 + \left(\frac{1}{\alpha} + \mu^2 \right) \tilde{y}^2 + \tilde{w}^2 dxdt +\nonumber\\
&+  \int_\Omega \nu |\nabla \tilde{y}(T)|^2 - \mu |\tilde{y}(T)|^2 dx.
\end{align}
Using Green's formula, the definition of $\tilde{w}$ and Young's inequality, we can estimate the second summand in \eqref{eqthm} by
\begin{equation*}
 \begin{array}{r c l}
  \displaystyle\int_0^T \int_\Omega - 2 \nu \mu |\nabla \tilde{y}|^2 dxdt  =  \displaystyle\int_0^T \int_\Omega 2 \nu \mu \Delta \tilde{y} \; \tilde{y} \;dxdt = \int_0^T \int_\Omega -2  \mu \tilde{w} \tilde{y} \; dxdt \\
  \geq \displaystyle\int_0^T \int_\Omega -2  \mu |\tilde{w}| \; |\tilde{y}| \; dxdt
   \geq  \displaystyle\int_0^T \int_\Omega - 4 \delta \mu^2 \tilde{y}^2 - \frac{1}{4\delta} \tilde{w}^2 \; dxdt.
 \end{array}
\end{equation*}
With the choice $\delta:= \displaystyle\frac{1+2\alpha\mu^2}{8 \alpha \mu^2}$, it holds that $ -4\delta \mu^2 + \frac{1}{\alpha} + \mu^2 \geq 0$ and $ -\frac{1}{4\delta} + 1 \geq 0$.\\

Using the Poincar\'{e} inequality, we can estimate the last term in \eqref{eqthm} by
$$\int_\Omega \nu |\nabla y(T)|^2 - \mu |y(T)|^2 dx \geq \left( \frac{\nu}{c_p^2} - \mu \right) \|  y(T) \|_{L^2(\Omega)}^2$$
with Poincar\'{e} constant $c_p$. If $\mu \leq \nu / c_p^2$, then $\displaystyle\int_\Omega \nu |\nabla y(T)|^2 - \mu |y(T)|^2 dx \geq 0$. Thus, for $\mu \leq \nu / c_p^2$ it holds that
$$ A_M^\mu((\tilde{y},\tilde{w}),(\tilde{y},\tilde{w})) \geq ||| (\tilde{y},\tilde{w})|||^2.$$
With this, the a-posteriori error estimate follows in analogy to Theorem~\ref{thm:errest}.
\end{proof}

\section{Model predictive control}\label{secmpc}

In this section, we recall the MPC method to solve the optimal control problem \eqref{ocp}. For a comprehensive study we refer to e.g.\ \cite{GruP17,RawMD09}. The basic idea of the MPC approach is to split the (large) time horizon $[0 ,T]$ of the optimal control problem \eqref{ocp} 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 < \dots < t_m:= T$ with fixed time step size $\Delta t = T/m$. In Section \ref{secadaptivity}, we will propose a time adaptive method. Let us define a time domain with small time horizon as $[t_i, t_i^N]$ for $i=0,\ldots,m-1$, where $ 1 < N  \leq m+1 $ is a chosen positive natural number such that $t_i^N:=t_i+(N-1)\Delta t$. We denote the length of the {\em prediction horizon} by $\bar{T}$, i.e.\ $t_i^N - t_i = \bar{T} = (N-1)\Delta t$. In this section and in the sequel, $N$ denotes the number of time instances in each subinterval and $\bar{T}$ the size of the subinterval. Those values might not be constant as we will see in the next section.

The reduced cost functional over the domain $[t_i, t_i^N] \times \Omega$ which is considered in the $i$-th iteration of the MPC algorithm for $i=0,\dots,m-1$ reads as follows:
\begin{equation}\label{mpccost}
\hat{J}^N(u,t_i,y_i):=\int_{t_i}^{t_i^N} \ell (y_{[u,t_i,y_i]}(t),u(t))dt, \quad i=0,\ldots,m-1.
\end{equation}
Then, the open loop control problem in each level of the MPC method is
\begin{equation}\label{mpcocpreduced}
 \min_{u \in L^2(t_i,t_i^N;\Omega)} \hat{J}^N(u,t_i,y_i),  \quad i=0,\ldots,m-1.
\end{equation}

The MPC method works as follows. We start by solving the optimal control problem \eqref{mpcocpreduced} in $[t_0, t_0^N] \times \Omega$ and we store the optimal control $u^N$ on the first subinterval $[t_0,t_0 + \Delta 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_1^N]$ with $t_1=t_0+\Delta t$ where the initial condition $y_1$ is given by the optimal trajectory  $y_{[u^N,t_0,y_0]}(t)$ at $t=t_0+\Delta t$ using the optimal control $u^N(t)$ for $t\in(t_0,t_0+\Delta t]$. We store the feedback map $\phi^N$ which is defined as $\phi^N(y_{[u^N, t_0, y_0]}(t)) := u^N(t)$. We iterate this process. Note that \eqref{mpcocpreduced} is an open loop problem on a finite time horizon $[t_i,t_i^N]$ for $i=0, \dots, m-1$ which can be treated numerically by classical techniques, as in e.g.\ \cite{HinPUU09}. The Algorithm~\ref{Alg:NMPC} summarizes the MPC technique for an equidistant time stepping $\Delta t$ and a fixed time horizon $N$.

\begin{algorithm}[h]
\caption{Model predictive control (MPC)}
\label{Alg:NMPC}
\begin{algorithmic}
\REQUIRE Number of degrees of freedom in time $m$, finite horizon $N\in\mathbb N$, initial condition $y_0$, desired state $y_d$, source $f$, space domain $\Omega$, time domain $[0,T]$, constants $\nu,\mu,\alpha$.
\STATE Set $t_0 = 0$.
\FOR {$i=0,1,2,\ldots, m-1$}
\STATE Compute the optimal control over $[t_i,t_i^N]$ with $t_i^N = t_i+(N-1)\Delta t$, solving
\begin{equation}\label{MPC:FinHor}
u^N:= \arg\min_{ u\in L^2(t_i,t_i^N;\Omega)} \hat{J}^N(u,t_i,y_i).
\end{equation}
\STATE Define the MPC feedback value $\phi^N(y_{[u^N,t_i,y_i]}(t)):=u^N(t)$, $t\in(t_i,t_i+\Delta t]$.
\STATE Compute the associated state $y^N=y_{[\phi^N,t_i,y_i]}(t)$ by solving \eqref{heatdepletion} on $(t_i,t_i+\Delta t]$.
\STATE Set $y_{i+1}=y^N(t_i+\Delta t)$ and $t_{i+1}=t_i+\Delta t$.
\ENDFOR
\end{algorithmic}
\end{algorithm}

In general, it is well-known that the larger the prediction horizon, the better the feedback law one can obtain, see e.g.\ \cite{GruP17}. 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 \cite[Chapter 6]{GruP17}).

\section{Time adaptivity in MPC}\label{secadaptivity}
In this section, we propose time adaptive techniques within MPC. In Algorithm \ref{Alg:NMPC} 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:

\begin{center}
 {\em  How to choose a suitable time discretization for the prediction horizon in each level of the MPC?}
\end{center}

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 \textit{optimal} state.

\noindent Therefore, we propose different approaches to combine the computation of an adaptive time grid using the a-posteriori estimate provided in Section \ref{sec2:timedisc} within a MPC framework. The idea of adaptivity leads to different combinations using the error estimate \eqref{est-thm31mu}. 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 \cite{GSSnew}.\\

\paragraph{\bf (1) Offline approach:  Compute an adaptive time grid a-priori before starting the MPC} We compute an adaptive time grid $\{t_i\}_{i=0}^{m+N-1}$ according to \eqref{est-thm31mu}, 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~\ref{rem:heuristic} we will use a very coarse spatial mesh to solve \eqref{semitimediscrete_depletion} 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 \textit{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~\ref{fig:MPCoffline}, where we use a chosen number of time instances $N$ in the subintervals $[t_i,t_{i+N-1}]$ for $i=0,\dots,m-1$. Note that this procedure allows prediction horizons with different lengths $t_{i+N-1} - t_i = \bar{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~\ref{Alg:NMPC_off}.\\

 \begin{figure}[h]
   \centering
    \begin{tikzpicture}[scale=1.5]
       \draw[gray, line width=0.5mm] (-2.5,0) -- (5,0);

    \fill[gray] (-2.5,-0.1) rectangle (-2.55,0.1);
    \fill[gray] (5,-0.1) rectangle (5.05,0.1);

         \fill[gray] (-2,-0.1) rectangle (-2.05,0.1);
         \fill[blue] (-1,-0.1) rectangle (-1.05,0.1);
         \fill[blue] (-0.25,-0.1) rectangle (-0.3,0);
         \fill[purple] (-0.25,0) rectangle (-0.3,0.1);
         \fill[blue] (-0.1,-0.1) rectangle (-0.15,0);
        \fill[purple] (-0.1,0) rectangle (-0.15,0.1);
         \fill[blue] (0,-0.1) rectangle (0.05,0);
         \fill[purple] (0,0) rectangle (0.05,0.1);
         \fill[blue] (0.1,-0.1) rectangle (0.15,0);
         \fill[purple] (0.1,0) rectangle (0.15,0.1);
         \fill[blue] (0.3,-0.1) rectangle (0.35,0);
         \fill[purple] (0.3,0) rectangle (0.35,0.1);
         \fill[purple] (2.1,-0.1) rectangle (2.15,0.1);
         \fill[gray] (2.3,-0.1) rectangle (2.35,0.1);
         \fill[gray] (2.5,-0.1) rectangle (2.55,0.1);
         \fill[gray] (4,-0.1) rectangle (4.05,0.1);
         \fill[gray] (4.5,-0.1) rectangle (4.55,0.1);

          \draw[line width=0.5mm] (-2.5,-0.2) node {\color{gray}$t_0=0$};
    \draw[line width=0.5mm] (5,-0.2) node {\color{gray}$t_m=T$};
    \draw[blue,line width=0.5mm] (-1,-0.6) -- (0.3,-0.6);
    \draw[blue,line width=0.5mm] (-1,-0.6) -- (-1,-0.45);
    \draw[blue,line width=0.5mm] (0.3,-0.6) -- (0.3,-0.45);
    \draw[line width=0.5mm] (-1,-0.3) node {\color{blue}$t_i$};
    \draw[line width=0.5mm] (0.5,-0.3) node {\color{blue}$t_{i+N-1}$};
      \draw[line width=0.5mm] (-0.2,0.3) node {\color{purple}$t_{i+1}$};
        \draw[line width=0.5mm] (2.2,0.3) node {\color{purple}$t_{i+N}$};
    \draw[line width=0.5mm] (-0.3,-0.85) node {\footnotesize \color{blue} prediction horizon};
     \draw[purple,line width=0.5mm] (-0.25,0.5) -- (2.1,0.5);
     \draw[purple,line width=0.5mm] (-0.25,0.35) -- (-0.25,0.5);
    \draw[purple,line width=0.5mm] (2.1,0.35) -- (2.1,0.5);
     \draw[line width=0.5mm] (1,0.7) node {\footnotesize \color{purple} shifted prediction horizon};
     \draw[line width=0.5mm] (2.8,0.2) node{$\cdots$};
       \end{tikzpicture}
       \caption{Scheme of the offline approach: For a given number of time instances $N$ we compute on a prescribed adaptive grid the MPC  subproblem in the interval $[t_i,t_{i+N-1}]$.}
       \label{fig:MPCoffline}
       \end{figure}
\noindent

 \begin{algorithm}[h]
\caption{Model predictive control (MPC) with offline adaptive grid }
\label{Alg:NMPC_off}
\begin{algorithmic}
\REQUIRE Number of degrees of freedom in time $m$, finite horizon $N\in\mathbb N$, initial condition $y_0$, desired state $y_d$, source $f$, space domain $\Omega$, time domain $[0,T]$, constants $\nu,\mu,\alpha$, $t_{m+N-1} > T$.
\STATE Compute an adaptive time grid $\{t_i\}_{i=0}^{m+N-1}$ using \eqref{est-thm31mu} with $t_0=0,t_m=T$.
\FOR {$i=0,1,2,\ldots,m-1$}
\STATE Compute the optimal control over $[t_i,t_{i+N-1}]$ solving
\begin{equation}\label{MPC:FinHor_off}
u^N:= \arg\min_{ u\in L^2(t_i,t_{i+N-1};\Omega)} \hat{J}^N(u,t_i,y_i).
\end{equation}
\STATE Define the MPC feedback value $\phi^N(y_{[u^N,t_i,y_i]}(t))=u^N(t)$,  $t\in(t_i,t_{i+1}]$.
\STATE Compute the associated state $y^N=y_{[\phi^N,t_i,y_i]}(t)$ by solving \eqref{heatdepletion} on $(t_i,t_{i+1}]$.
\STATE Set $y_{i+1}=y^N(t_{i+1})$.
\ENDFOR
\end{algorithmic}
\end{algorithm}

\paragraph{\bf (2) Online approach: Compute an adaptive time grid for each time horizon in the MPC iteration} For a given prediction interval $[t^i_0,t^i_{N-1}]$ at each MPC iteration $i$, we make use of the a-posteriori error estimation {\color{black} \eqref{est-thm31mu}} for the state to compute an adaptive time grid within the current time horizon. Note that $t_0^0 := 0$ is the initial time. The scheme is visualized in Figure~\ref{fig:MPConline}.

  \begin{figure}[htbp]
   \centering
    \begin{tikzpicture}[scale=1.8]
       \draw[gray,line width=0.5mm] (-1.3,0) -- (5,0);

         \fill[blue] (-0.8,-0.1) rectangle (-0.85,0.1);
          \fill[blue] (0.1,-0.1) rectangle (0.15,0);
          \fill[purple] (0.1,0) rectangle (0.15,0.1);
          \fill[purple] (0.46,-0.1) rectangle (0.51,0.1);
          \fill[blue] (0.52,-0.1) rectangle (0.57,0.1);
          \fill[purple] (0.8,-0.1) rectangle (0.85,0.1);
          \fill[blue] (0.87,-0.1) rectangle (0.92,0.1);
          \fill[blue] (1.6,-0.1) rectangle (1.65,0.1);
         \fill[purple] (1.8,-0.1) rectangle (1.85,0.1);
          \fill[purple] (2.8,-0.1) rectangle (2.85,0.1);
          \fill[blue] (3,-0.1) rectangle (3.05,0.1);
          \fill[purple] (3.92, -0.1) rectangle (3.97,0.1);
          \draw[line width=0.5mm] (-0.9,-0.25) node {\color{blue} $t^i_0$};
          \draw[line width=0.5mm] (3,-0.25) node {\color{blue} $t^i_{N-1}$};
          \draw[blue,line width=0.5mm] (-0.8,-0.5) -- (3.02,-0.5);
          \draw[blue,line width=0.5mm] (-0.8,-0.5) -- (-0.8,-0.35);
          \draw[blue,line width=0.5mm] (3.02,-0.5) -- (3.02,-0.35);
          \draw (1,-0.7) node {\footnotesize \color{blue} prediction horizon};
          \draw[purple,line width=0.5mm] (0.1,0.55) -- (3.92,0.55);
          \draw[purple,line width=0.5mm] (0.1,0.4) -- (0.1,0.55);
          \draw[purple,line width=0.5mm] (3.92,0.4) -- (3.92,0.55);
          \draw (2,0.7) node {\footnotesize \color{purple} shifted prediction horizon };
          \draw (4, 0.25) node{\color{purple} $t^{i+1}_{N-1}$};
          \draw (0.15, 0.25) node{\color{purple} $t_1^i = t^{i+1}_{0}$};
       \end{tikzpicture}
       \caption{Scheme of the online approach: The blue color refers to the grid at iteration $i$ starting at time $t^i_0$ till $t^i_{N-1}$ whereas the red color refers to the next MPC level $i+1$. We note that the only guaranteed overlap of the time grids is for second time instance at iteration $i$ which corresponds to the first time instance at iteration $i+1$.}\label{fig:MPConline}
       \end{figure}

 For a given number of degrees of freedom $N$ the algorithm distributes the time instances according to the error estimation {\color{black}\eqref{est-thm31mu}} in subintervals of the same length $t^i_{N-1}-t^i_0 = \bar{T}$. The resulting adaptive time grid in each iteration level $i$ is related to the \textit{optimal} state of a corresponding open loop problem. Again, we assume that the heuristic assumptions of Remark~\ref{rem:heuristic} hold true which enables an efficient computation. The approach is summarized below in Algorithm~\ref{Alg:NMPC_on}.

 \begin{algorithm}[H]
\caption{Model predictive control (MPC) with online adaptive grid }
\label{Alg:NMPC_on}
\begin{algorithmic}
\REQUIRE Number of time instances in each subinterval $N$, prediction horizon $\bar{T}$, initial condition $y_0$, desired state $y_d$, source $f$, space domain $\Omega$, time domain $[0,T]$, constants $\nu,\mu,\alpha$.
\STATE Set $t_0^0=0, y_0^0=y_0$ and $i=0$.
\WHILE{$t_0^i < T$}
\STATE  Compute an adaptive grid $\{t^i_{j}\}_{j=0}^{N-1}$ in $[t_0^i, t_{N-1}^i:=t_0^i+\bar{T}]$ using \eqref{est-thm31mu}.
\STATE Compute the optimal control over $[t_0^i,t_{N-1}^i]$ solving
\begin{equation}\label{MPC:FinHor_on}
u^N:= \arg\min_{ u\in L^2(t_0^i,t_{N-1}^i;\Omega)} \hat{J}^N(u,t_{0}^i,y_0^i).
\end{equation}
\STATE Define the MPC feedback value $\phi^N(y_{[u^N,t_0^i,y_0^i]}(t))=u^N(t)$, $t\in(t_{0}^i,t_1^i]$.
\STATE Compute the associated state $y^N=y_{[\phi^N,t^i_0,y_0^i]}(t)$ by solving \eqref{heatdepletion}  on $(t_0^i,t_1^i]$.
\STATE Set $t_0^{i+1}:=t_1^i$, $y_0^{i+1}:=y^N(t_0^{i+1}), i \leftarrow i+1$.
\ENDWHILE
\end{algorithmic}
\end{algorithm}

 \begin{rmk}[Warm start]
  In order to make computations even more efficient, the information of the previous MPC iteration can be used as a \textit{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 {\color{black} one can use} as initial control the one computed at the previous step.
 \end{rmk}

\begin{rmk}[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 {\color{black} \eqref{est-thm31mu}} is used in an adaptive cycle for the {\color{black}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.
\end{rmk}

\section{Numerical example}\label{secnumerics}
In the following numerical tests, we compare the different adaptive time schemes within MPC proposed in Section~\ref{secadaptivity}. We will investigate numerically the influence of the time grid on the approximation quality of the MPC methods. In particular, we compare Algorithms~\ref{Alg:NMPC}, \ref{Alg:NMPC_off} and \ref{Alg:NMPC_on}.

In all numerical examples, the considered spatial domain is the open interval $\Omega = (0,1)$ and the temporal domain is $[0 ,T] = [0,1]$. In order to solve the mixed form \eqref{2ord_decoupled_depletion}, we introduce a partitioning of the space-time domain into regular orthotopes and use $\mathbb{Q}_1$ space-time finite elements for discretization, where $\mathbb{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 $\Delta x = 1/5$. \\
For the solution of the MPC open loop subproblems, we use an implicit Euler scheme for the temporal discretization and use $\mathbb{P}_1$ finite elements for the spatial discretization. The optimal control problem is solved {\color{black} with a direct solver}. We take as fine spatial resolution an equidistant discretization with $\Delta x = 1/100$. All coding is done in \textsc{Matlab R2019}a

\subsection{Test 1: Solution with a layer at $t=0.5$}
In this numerical test, we consider the optimal control of \eqref{heat}. The setting for this test example is taken from \cite[Example 5.2]{GonHZ12}, with the following choices: $\nu = 1$ in \eqref{heat} and $\alpha = 1$ in \eqref{cost}. The example is built such that the exact optimal solution $(y,u)$ {\color{black}to \eqref{ocpreduced}} is given by
$$ y(t,x) = \sin (\pi x) \text{atan} ((t-1/2)/\varepsilon), \quad u(t,x) = -\sin (\pi x) \sin (\pi t).$$
The initial condition is $y_0 (x) = \sin(\pi x) \text{atan}(-1/(2\varepsilon))$. The functions $f$ and $y_d$ are chosen accordingly as
\begin{align*}
\centering
f(t,x) & = \sin(\pi x) \left( \varepsilon/(t^2 - t + \varepsilon^2 + 1/4) + \pi^2 \text{atan}((t-1/2)/(\varepsilon)) + \sin(\pi t) \right),\\
y_d(t,x)& = \sin(\pi x) \left( \text{atan}((t-1/2)/(\varepsilon)) + \pi \cos(\pi t) - \pi^2 \sin (\pi t) \right).
\end{align*}
For small values of $\varepsilon$ (we use $\varepsilon = 10^{-3}$), the state $y$ develops a very steep gradient at $t = 0.5$, which can be seen in Figure~\ref{fig:ytrue}.

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.3]{test1ytrue.pdf}
\caption{Test 1: True optimal state solution}
\label{fig:ytrue}
\end{figure}

We now compare Algorithms~\ref{Alg:NMPC}, \ref{Alg:NMPC_off} and \ref{Alg:NMPC_on}. For an exemplary visualization, let us consider the following choices in Algorithm~\ref{Alg:NMPC_off}: $m=45, N=9$. For the online approach in Algorithm~\ref{Alg:NMPC_on} we pick $\bar{T}=0.2, N=9$ as comparison. The adaptive approaches are compared with the equidistant MPC Algorithm~\ref{Alg:NMPC} with $N=9, m=45$, i.e.\ $\bar{T}=0.2$. The numerical state solutions of the controlled problem with the different MPC approaches are shown in Figure~\ref{fig:off_y_MPC}. We can see that the standard MPC Algorithm~\ref{Alg:NMPC} with equidistant time grid fails
whereas using either of the Algorithms~\ref{Alg:NMPC_off} or \ref{Alg:NMPC_on}, it is possible to capture the layer at $t=0.5$ and the solutions comply much better with the true state solution

\begin{figure}[htbp]
\centering
  \includegraphics[scale=0.3]{test1y_equi_offline.pdf} \hspace{-0.4cm}
 \includegraphics[scale=0.3]{test1y_ad_offline.pdf} \hspace{-0.4cm}
    \includegraphics[scale=0.3]{test1y_ad_online.pdf}
 \caption{Test 1: MPC state solution $y$ using a uniform time discretization (left), offline adaptive approach (middle) and online adaptive scheme (right)}
 \label{fig:off_y_MPC}
\end{figure}

Let us now provide more details about the grids we obtained with the methods proposed. In the offline approach (Algorithm \ref{Alg:NMPC_off}), the temporal grid is computed before starting the MPC algorithm. The adaptive grids using a coarse and a fine spatial resolution for solving \eqref{2ord_decoupled} are shown, respectively, in the middle and right panel of Figure~\ref{fig:offlinegrid}. As already discussed in \cite{AllGH18,GonHZ12} and recalled in Remark~\ref{rem:heuristic}, we observe numerically that the temporal and spatial resolution decouple and we are able to compute cheaply the temporal grid by means of a coarse spatial resolution for the solution of \eqref{2ord_decoupled}.

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.3]{test1unigrid.pdf} \hspace{-0.4cm}
\includegraphics[scale=0.3]{test1coarsegrid.pdf} \hspace{-0.4cm}
\includegraphics[scale=0.3]{test1finegrid.pdf}
 \caption{Test 1: Uniform space-time grid with fine spatial resolution (left), offline adaptive grid with coarse (middle) and fine (right) spatial resolution}\label{fig:offlinegrid}
\end{figure}

Examples of the corresponding equidistant and adaptive temporal intervals are compared in the top and bottom panels of Figure~\ref{fig:offlineadgrids}. 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 \ref{fig:offlineadgrids} it is clear how the length of the prediction horizon $\bar{T}_j$ is different for each iteration of the MPC method.

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.3]{test1timegrid5_ad_offline.pdf} \hspace{-0.4cm}
\includegraphics[scale=0.3]{test1timegrid30_ad_offline.pdf} \hspace{-0.4cm}
\includegraphics[scale=0.3]{test1timegrid35_ad_offline.pdf}\hspace{-0.4cm}
\includegraphics[scale=0.3]{test1timegrid5_equi_offline.pdf} \hspace{-0.4cm}
\includegraphics[scale=0.3]{test1timegrid30_equi_offline.pdf} \hspace{-0.4cm}
\includegraphics[scale=0.3]{test1timegrid35_equi_offline.pdf}
\caption{Test 1: Adaptive time horizons according to the offline approach (top), uniform time horizons according to the standard MPC approach (bottom), MPC iteration levels $i=5$ (left), $i=30$ (middle), $i=35$ (right)}\label{fig:offlineadgrids}
\end{figure}

The online adaptive grid with a coarse and a fine spatial resolution is shown in the middle and right panel of Figure \ref{fig:onlinegrid}. Again, we observe that the time adaptivity is very insensitive with respect to the spatial resolution.

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.3]{test1spacetimeuniformonline.pdf}\hspace{-0.4cm}
\includegraphics[scale=0.3]{test1coarsegridonline.pdf}\hspace{-0.4cm}
\includegraphics[scale=0.3]{test1finegridonline.pdf}
 \caption{Test 1: Uniform space-time grid with fine spatial resolution (left), online adaptive grid with coarse (middle) and fine (right) spatial resolution}\label{fig:onlinegrid}
\end{figure}

The corresponding adaptive time horizons are shown in the top panels of Figure~\ref{fig:onlineadgrids}. As a comparison, the uniform time horizons of the same lengths are shown in the bottom panels of Figure~\ref{fig:onlineadgrids} using the same number of degrees of freedom in each interval. It is clear that the a-posteriori error estimate \eqref{est-thm31} leads to a time grid associated with the open loop optimal state which benefits the accuracy of the control problem.

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.3]{test1timegrid13_ad_online.pdf}\hspace{-0.4cm}
\includegraphics[scale=0.3]{test1timegrid15_ad_online.pdf}\hspace{-0.4cm}
\includegraphics[scale=0.3]{test1timegrid24_ad_online.pdf}\\
\includegraphics[scale=0.3]{test1timegrid13_equi_online.pdf}\hspace{-0.4cm}
\includegraphics[scale=0.3]{test1timegrid15_equi_online.pdf}\hspace{-0.4cm}
\includegraphics[scale=0.3]{test1timegrid24_qui_online.pdf}
\caption{Test 1: Adaptive time horizons according to the online approach (top), uniform time horizons according to the standard MPC approach (bottom), MPC iteration levels $i=13$ (left), $i=15$ (middle), $i=24$ (right)}\label{fig:onlineadgrids}
\end{figure}

Finally, we provide an error analysis for the computation of the approximate solutions using Algorithms~\ref{Alg:NMPC}, \ref{Alg:NMPC_off} and \ref{Alg:NMPC_on} {\color{black}for different choices of degrees of freedom in time and prediction horizons}. {\color{black}  For this,} we compute the error between the analytical optimal state solution to \eqref{ocpreduced} and its numerical approximation using the different MPC approaches measured in the $L^2(0,T;\Omega)-$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 \ref{fig:costoff} we show the error for Algorithms~\ref{Alg:NMPC} and \ref{Alg:NMPC_off}. 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 {\color{black}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~\ref{Alg:NMPC}) 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 \ref{fig:costoff} we see how the solution converges increasing the number of degrees of freedom.  As expected the larger the prediction horizon the better the approximation.

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.3]{test1abserr_equi_offline.pdf}
   \includegraphics[scale=0.3]{test1abserr_ad_offline.pdf}
  \caption{Test 1: $L^2-$error for the offline approach with equidistant (left) and adaptive (right) time grid} \label{fig:costoff}
  \end{figure}

In Figure \ref{fig:timeoff} we compare the computational time (measured in seconds) for Algorithms~\ref{Alg:NMPC} and  \ref{Alg:NMPC_off}. 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.

  \begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.3]{test1time_equi_offline.pdf}
   \includegraphics[scale=0.3]{test1time_ad_offline.pdf}
  \caption{Test 1: Computational time in seconds for an offline grid with equidistant (left) and adaptive (right) time grid}
  \label{fig:timeoff}
  \end{figure}

 {\color{black} To conclude our comments on the offline approach we show the results of the terms $ \|y - y_d\|_{L^2(0,T;\Omega)}^2$ in the top panels and $\|u\|_{L^2(0,T;\Omega)}^2$ in the bottom panels of Figure~\ref{fignew:cost} 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\|_{L^2(0,T;\Omega)}^2 = 26.8197$ and $\|u_{true}\|_{L^2(0,T;\Omega)}^2 = 0.25$ (``Analytical''). Moreover, Figure~\ref{fignew:costexemplary} compares the evolution of $\|y(t)-y_d\|_{L^2(\Omega)}^2$ and $\|u(t)\|_{L^2(\Omega)}^2$ over time comparing the MPC state and control solution with the true solution evaluated at the given time grid.}

  \begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.3]{test1yyd_equi_offline.pdf}
   \includegraphics[scale=0.3]{test1yyd_ad_offline.pdf}\\
    \includegraphics[scale=0.3]{test1u_equi_offline.pdf}
   \includegraphics[scale=0.3]{test1u_ad_offline.pdf}
  \caption{{\color{black}Test 1: Value of the tracking term $\|y-y_d\|_{L^2(0,T;\Omega)}^2$ (top) and control cost $\|u\|_{L^2(0,T;\Omega)}^2$ (bottom) comparing the MPC solution for different horizon lengths using the equidistant (left) and offline adaptive approach (right) with the optimal solution}}
  \label{fignew:cost}
  \end{figure}

     \begin{figure}[htbp]
   \centering
   \includegraphics[scale=0.3]{test1yydovertime_equi_offline.pdf}
    \includegraphics[scale=0.3]{test1yydovertime_ad_offline.pdf}\\
       \includegraphics[scale=0.3]{test1uovertime_equi_offline.pdf}
   \includegraphics[scale=0.3]{test1uovertime_ad_offline.pdf}
   \caption{{\color{black}Test 1: Value of the tracking term $\|y(t)-y_d(t)\|_{L^2(\Omega)}^2$ (top) and control cost $\|u(t)\|_{L^2(\Omega)}^2$ (bottom) comparing the MPC solution for $m=45,N=9$ using the equidistant (left) and offline adaptive approach (right); the vertical lines indicate the utilized time grid}}
   \label{fignew:costexemplary}
   \end{figure}

In Figure \ref{fig:coston} we compare the $L^2$-error for Algorithms~\ref{Alg:NMPC} and \ref{Alg:NMPC_on}. We fixed the prediction horizon $\bar{T}$ 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 $\bar{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.

  \begin{figure}[htbp]
  \centering
     \includegraphics[scale=0.3]{test1abserr_equi_online.pdf}
  \includegraphics[scale=0.3]{test1abserr_ad_online.pdf}
     \caption{Test 1: $L^2-$error for an online approach with equidistant (left) and adaptive (right) time grid}
  \label{fig:coston}
\end{figure}

In Figure \ref{fig:timeon} we compare the {\color{black} computational} time in seconds of Algorithm~\ref{Alg:NMPC_on} 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 {\color{black} $0.0872$ computed in $25.87$s} whereas, with the adaptive approach, to get an error of {\color{black}$0.0216$ we needed $16.06s$}. This shows that our method is more accurate and also more efficient computationally without any a-priori knowledge of the control problem.

  \begin{figure}[htbp]
\centering
     \includegraphics[scale=0.3]{test1time_equi_online.pdf}
  \includegraphics[scale=0.3]{test1time_ad_online.pdf}
     \caption{Test 1: {\color{black}Computational} time in seconds for the online approach with equidistant (left) and adaptive (right) time grid. }
  \label{fig:timeon}
\end{figure}

\subsection{Test 2: State equation with depletion term}

In this numerical test, we consider the optimal control of \eqref{heatdepletion} with $\mu > 0$. The regularization parameter in the cost is chosen as $\alpha = 1$ and the desired state is given by
$$y_d(t,x) = 10\sin(\pi x)\text{exp}(-((t-0.5)/\varepsilon)^2)  - x(x-1)-2\nu(t-1)-\mu x(x-1)(t-1).$$
The initial condition for the state is
$$y_0(x) \equiv 10\sin(\pi x)\text{exp}(-1/(4\varepsilon^2))$$
and the source term in the state equation is chosen in order to get as exact optimal {\color{black} solution $(y,u)$} {\color{black}for the original control problem}:
$$y(x,t) = 10\sin(\pi x)\text{exp}(-((t-0.5)/\varepsilon)^2),\quad {\color{black}u(x,t) = x(x-1)(t-1).}$$
Let us note that the Poincar\'{e} constant $c_p$ and the first eigenvalue $\lambda_1$ of the Laplace-Dirichlet operator are related by $\lambda_1 = 1/c_p^2$ (see e.g \cite[Proposition 8.4.3]{AttBM14}). For $\Omega = (0,1)$, the first eigenvalue $\lambda_1$ of the Laplace-Dirichlet operator is given by $\lambda_1 = \pi^2$ (see e.g. \cite[Proposition 8.5.2]{AttBM14}). Then, since Theorem~\ref{thm:errest_depl} is applicable if $\mu \leq \nu / c_p^2$, for this setting it requires $\mu \leq 0.1 \cdot \pi^2$. In our numerical simulation we set $\nu =0.1, \mu = 3, \varepsilon = 10^{-2}$. Thus, we show the unstable case for $\mu=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 $\mu\leq \nu\pi^2$ as required in Theorem ~\ref{thm:errest_depl}.

The optimal solution is shown on the left panel of Figure \ref{fig2:opt}. As one can see we have a layer at $t=0.5$. Then we show the results of Algorithm \ref{Alg:NMPC_off} choosing $N=5,m=45$ with an equidistant grid (middle panel of Figure \ref{fig2:opt}) and an adaptive grid (right panel of Figure \ref{fig2:opt}). It is easy to see that the adaptive grid catches better the behavior of the optimal solution.

\begin{figure}[htbp]
\centering
  \includegraphics[scale=0.3]{test2ytrue.pdf}\hspace{-0.4cm}
   \includegraphics[scale=0.3]{test2y_equi_offline.pdf}\hspace{-0.4cm}
      \includegraphics[scale=0.3]{test2y_ad_offline.pdf}
      \label{fig2:opt}
\caption{Test 2: Optimal {\color{black}state} solution (left), MPC solution using an equidistant {\color{black}time} grid (middle), MPC solution with an offline {\color{black}adaptive} approach (right).}
\end{figure}

The grid for this {\color{black} setting} is shown in Figure \ref{fig2:grid} (right) and, as expected, our adaptive approach refines consistently in a neighborhood of the layer.

\begin{figure}[htbp]
 \centering
  \includegraphics[scale=0.3]{test2grid_equi_offline.pdf}
 \includegraphics[scale=0.3]{test2grid_ad_offline.pdf}
 \label{fig2:grid}
  \caption{Test 2, $\mu = 3$: Equidistant (left) and offline adaptive grid (right) with $m = 45, N=5$}
\end{figure}

The quality of our results are also confirmed by the $L^2-$error plots in Figure \ref{fig2:err}. Under the settings considered we can see that the error with an equidistant grid is of order, {\color{black}around, $O(10^0)$ whereas with an adaptive grid {\color{black}we can reach} $O(10^{-1})$. Clearly decreasing $\Delta x$ we can improve the accuracy of our approximation.

 \begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.3]{test2abserr_equi_offline.pdf}
   \includegraphics[scale=0.3]{test2abserr_ad_offline.pdf}
 \caption{Test 2, $\mu = 3$: Value of the $L^2-$error for the MPC solution trajectory with different number of degrees of freedom in the whole time domain $[0,T]$ and increasing time horizon lengths; using the offline approach}
 \label{fig2:err}
 \end{figure}

Finally, we show the error behavior of the online approach in Figure \ref{fig2:err2}. For page limitation we do not show a complete analysis of the results with Algorithm \ref{Alg:NMPC_on}. However, considerations are very similar to the previous test.

 \begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.3]{test2abserr_equi_online.pdf}
   \includegraphics[scale=0.3]{test2abserr_ad_online.pdf}
 \caption{Test 2, $\mu = 3$: Value of the $L^2-$error for the MPC solution trajectory with different number of degrees of freedom within the prediction horizon and increasing time horizon lengths; using the online approach}
  \label{fig2:err2}
 \end{figure}

\section{Conclusions and Outlook}\label{secconcl}

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.

\bibliographystyle{plain}