Crank-Nicolson time stepping and variational discretization of control-constrained parabolic optimal control problems

The purpose of this paper is to prove optimal a priori error bounds for the variational time semi-discretization of a generic parabolic optimal control problem, where the state in time is approximated with a Petrov Galerkin scheme. The key idea consists in choosing piecewise linear, continuous test functions and a discontinuous, piecewise constant Ansatz for the approximation of the state equation. With this Petrov Galerkin Ansatz variational discretization of the optimal control problem delivers a cG(1) time approximation of the optimal time semi-discrete adjoint state. The resulting time integration schemes for the state and the adjoint state are variants of the Crank-Nicolson scheme. Combining this setting with the supercloseness result of Corollary 4.3 for interval means we are able to prove second order in time convergence of the time discrete optimal control in Theorem 5.2. Moreover, the piecewise linear and continuous parabolic projection of the discrete state based on the the values of the discrete state on the dual time grid provides a second order convergent approximation of the optimal state without further numerical effort, see Theorem 5.3.

Our work is motivated by the work [MV11] of Meidner and Vexler, whoms technical results we use whenever possible. Under mild assumptions on the active set they show the same convergence order in time for the post-processed piecewise linear, continuous parabolic projection of the piecewise constant in time optimal control. This control is obtained by a Petrov Galerkin scheme with variational discretization of the parabolic optimal control problem. In comparison to their work, we switch Ansatz and test space in our numerical schemes.

Our work is novel in several aspects:

• •

  In (4.1) to the best of the authors knowledge a new, fully variational time-discretization scheme for the parabolic state equation is presented. It results in a Crank-Nicolson scheme with initial damping step for the nodal values of the state.

• •

  In Theorem 5.2 we provide an optimal error estimate for the time-discrete control, namely

  $$\left\|{\bar{u}}-{\bar{u}_{k}}\right\|_{L^{2}(0,T,\mathbb{R}^{D})}\leq Ck^{2},$$

  with ${\bar{u}}$ the optimal control, which is the solution of the optimal control problem ($\mathbb{P}$), and ${\bar{u}_{k}}$ denoting the optimal control obtained from the related discretized problem ($\mathbb{P}_{k}$) (see below). Here, $k$ denotes the grid size of the time grid. This result could be compared to [MV11, Theorem 6.2]. There, under mild assumptions on the structure of the active set w.r.t. ${\bar{u}}$, a similar bound is obtained for the post-processed parabolic projection of a piecewise constant optimal control. Our approach avoids such an assumption in the numerical analysis, and presents an error estimate for the variational-discrete optimal control.

• •

  In Theorem 5.3 we prove

  $$\left\|{\bar{y}}-\pi_{P_{k}^{*}}{\bar{y}_{k}}\right\|_{L^{2}(0,T,L^{2}(\Omega))}\leq Ck^{2},$$

  where $\pi_{P_{k}^{*}}{\bar{y}_{k}}$ denotes the piecewise linear and continuous parabolic projection of the discrete state based on the the values of the discrete state on the dual time grid defined in Section 4. Since these values are already known, the projection is for free.

• •

  Our approach demonstrates that variational discretization of [Hin05] through the choice of Ansatz and test space in a Petrov Galerkin approximation of parabolic optimal control problems offers the possibility to specify the discrete structure of variational optimal controls. Of course this feature applies also to other classes of PDE constrained optimal control problems.

In our note we only consider semi-discretization in time for two reasons. On the one hand, error estimation for standard spatial finite element approximations of the time-semidiscrete optimal control problem ($\mathbb{P}_{k}$) is along the lines of [MV11, Section 6.2]. On the other hand, we are interested in the approximation of optimal controls, which in a realistic time-dependent scenario only depend on time, see the possible definitions of the control operator $B$ below.

With $I:=(0,T)\subset\mathbb{R}$, $T<\infty$, and a fixed function $y_{d}\in{L^{2}(I,{L^{2}(\Omega)})}$, we consider the linear-quadratic optimal control problem

		$\displaystyle\min_{y\in Y,u\in{U_{\textup{ad}}}}J(y,u)=\frac{1}{2}\|y-y_{d}\|^{2}_{{L^{2}(I,{L^{2}(\Omega)})}}+\frac{\alpha}{2}\|u\|^{2}_{U},$		($\mathbb{P}$)
		$\displaystyle\text{s.t. }y=S(Bu,y_{0}).$

The state space $Y$ is given by

$$Y:=W(I):=\{v\in L^{2}(I,{H^{1}_{0}(\Omega)}),\partial_{t}v\in L^{2}(I,{H^{-1}(\Omega)})\}\hookrightarrow C([0,T],L^{2}(\Omega)),$$

and the operator $S:L^{2}(I,{H^{-1}(\Omega)})\times L^{2}(\Omega)\rightarrow W(I)$, $(f,\kappa)\mapsto y:=S(f,\kappa)$, denotes the weak solution operator associated with the parabolic problem

	$\displaystyle\partial_{t}y-\Delta y$	$\displaystyle=f$		$\displaystyle\text{in }I\times\Omega\,,$		(1.1)
	$\displaystyle y$	$\displaystyle=0$		$\displaystyle\text{in }I\times\partial\Omega\,,$
	$\displaystyle y(0)$	$\displaystyle=\kappa$		$\displaystyle\text{in }\Omega,$

i.e. for $(f,\kappa)\in L^{2}(I,{H^{-1}(\Omega)})\times L^{2}(\Omega)$ the function $y\in W(I)$ satisfies $y(0)=\kappa$ and

$$\int\limits_{0}^{T}\langle\partial_{t}y(t),v(t)\rangle_{{H^{-1}(\Omega)}{H^{1}_{0}(\Omega)}}+a(y(t),v(t))\,dt=\int\limits_{0}^{T}\langle f(t),v(t)\rangle_{{H^{-1}(\Omega)}{H^{1}_{0}(\Omega)}}\,dt\\ \quad\forall\,v\in L^{2}(I,{H^{1}_{0}(\Omega)}).$$

(1.2)

Here $\Omega\subset\mathbb{R}^{n}$, $n=2,3$, is a convex polygonal domain with boundary $\partial\Omega$, and for $y,v\in{H^{1}_{0}(\Omega)}$ we define

$$a(y,v):=\int\limits_{\Omega}\nabla y(x)\nabla v(x)\ dx.$$

In what follows several choices of control spaces are feasible. In particular we may choose as control space $U=L^{2}(I,\mathbb{R}^{D})$, $D\in\mathbb{N}$, and as the admissible set

$${U_{\textup{ad}}}=\left\{u\in U\,\left|\;a_{i}\leq u_{i}(t)\leq b_{i}\text{ a.e. in $I$, }i=1,\dots,D\right.\right\},$$

where $a_{i}$, $b_{i}\in\mathbb{R}$, $a_{i}<b_{i}$ $(i=1,\dots,D)$. In this case the control operator is given by

$$B:U\rightarrow L^{2}(I,{H^{-1}(\Omega)})\,,\quad u\mapsto\left(t\mapsto\sum_{i=1}^{D}u_{i}(t)g_{i}\right),$$

(1.3)

where $g_{i}\in{H^{-1}(\Omega)}$ are given functionals, whose regularity is specified in Assumption 1.1 below. Clearly, $B$ is linear and bounded. A further possible choice for the control space is $U={L^{2}(I,{L^{2}(\Omega)})}$ with

$${U_{\textup{ad}}}=\left\{u\in U\,\left|\;a\leq u(t,x)\leq b\text{ a.e. in }I\times\Omega\right.\right\},$$

where $a<b$ denote real constants. In this case the control operator $B$ is the injection from $U$ into $L^{2}(I,H^{-1}(\Omega))$. In both cases the admissible set ${U_{\textup{ad}}}$ is closed and convex. We build our exposition upon the practical more relevant first choice of time-dependent amplitudes as controls.

It is well known that the operator $S$ is well defined, i.e. for every $(f,\kappa)\in L^{2}(I,{H^{-1}(\Omega)})\times L^{2}(\Omega)$ a unique state $y\in W(I)$ satisfying (1.2) exists. Furthermore, it fulfills

$$\|y\|_{W(I)}\leq C\left\{\|f\|_{L^{2}(I,{H^{-1}(\Omega)})}+\|\kappa\|_{L^{2}(\Omega)}\right\}.$$

Now let $y\in Y$ denote the unique solution of (1.2), and let $v\in W(I)$. Then it follows from integration by parts for functions in $W(I)$, that with the bilinear form $A:W(I)\times W(I)\rightarrow\mathbb{R}$ defined by

$$A(y,v):=\int\limits_{0}^{T}-\langle\partial_{t}v(t),y(t)\rangle_{{H^{-1}(\Omega)}{H^{1}_{0}(\Omega)}}+a(y(t),v(t))\,dt+(y(T),v(T))_{{L^{2}(\Omega)}}\,,$$

(1.4)

the state $y$ also satisfies

$$A(y,v)=\int\limits_{0}^{T}\langle f(t),v(t)\rangle_{{H^{-1}(\Omega)}{H^{1}_{0}(\Omega)}}\,dt+(\kappa,v(0))_{{L^{2}(\Omega)}}\quad\forall\ v\in W(I).$$

(1.5)

Furthermore, $y$ is the only function in $Y$ which satisfies (1.5). In the next section we use the bilinear form $A$ to define our numerical approximation scheme for the state equation.

A lot of literature is available on optimal control problems with parabolic state equations. We refer to [HPUU09] for a comprehensive discussion, and also to [AF12, MV08a, MV08b, MV11, SV13] for the most recent developments related to optimal control with Galerkin methods in time.

The paper is organized as follows. In section 2 we briefly summarize the solution theory of the optimal control problem. In section 3 we analyse the regularity of the state and the adjoint state, which plays an important role in the time discretization. In section 4 the time discretization of state and adjoint state is discussed in detail. In section 5 we introduce variational discretization of the optimal control problem $(P)$ and prove second order convergence of the variational discrete controls in time. In section 6 we present numerical results which confirm our analytical findings.

It is well known that problem ($\mathbb{P}$) admits a unique solution $({\bar{y}},{\bar{u}})\in Y\times U$, where ${\bar{y}}=S(B{\bar{u}},y_{0})$. Moreover, using the orthogonal projection $P_{U_{\textup{ad}}}:L^{2}(I,\mathbb{R}^{D})\rightarrow{U_{\textup{ad}}}$, the optimal control is characterized by the first-order necessary and sufficient condition

$${\bar{u}}=P_{U_{\textup{ad}}}\left(-\frac{1}{\alpha}B^{\prime}\bar{p}\right),$$

(2.1)

where $(\bar{p},\bar{q})\in L^{2}(I,{H^{1}_{0}(\Omega)})\times{L^{2}(\Omega)}$ (here we use reflexivity of the involved spaces) denotes the adjoint variable which is the unique solution to

$$\int\limits_{0}^{T}\langle\partial_{t}\tilde{y}(t),\bar{p}(t)\rangle_{{H^{-1}(\Omega)}{H^{1}_{0}(\Omega)}}+a(\tilde{y}(t),\bar{p}(t))\,dt+(\tilde{y}(0),\bar{q})_{{L^{2}(\Omega)}}\\ =\int\limits_{0}^{T}\int\limits_{\Omega}({\bar{y}}(t,x)-y_{d}(t,x))\tilde{y}(t,x)\,dxdt\quad\forall\ \tilde{y}\in W(I).$$

(2.2)

Here, $B^{\prime}:L^{2}(I,{H^{1}_{0}(\Omega)})\rightarrow L^{2}(I,\mathbb{R}^{D})$ denotes the adjoint operator of $B$, which is characterized by

$$B^{\prime}q(t)=\left(\langle g_{1},q(t)\rangle_{{H^{-1}(\Omega)}{H^{1}_{0}(\Omega)}}\ ,\ \dots\ ,\ \langle g_{D},q(t)\rangle_{{H^{-1}(\Omega)}{H^{1}_{0}(\Omega)}}\right)^{T}.$$

(2.3)

Furthermore we note that for $v\in L^{2}(I,\mathbb{R}^{D})$ there holds

$$P_{U_{\textup{ad}}}(v)(t)=\left(P_{[a_{i},b_{i}]}(v_{i}(t))\right)_{i=1}^{D},$$

where for $a,b,z\in\mathbb{R}$ with $a\leq b$ we set $P_{[a,b]}(z):=\max\{a,\min\{z,b\}\}$.

Since ${\bar{y}}-y_{d}\in L^{2}(I,L^{2}(\Omega))$ in (2.2), we have $\bar{p}\in W(I)$, so that by integration by parts for functions in $W(I)$ we conclude from (2.2) (compare (1.5))

$$\int\limits_{0}^{T}-\langle\partial_{t}\bar{p}(t),\tilde{y}(t)\rangle_{{H^{-1}(\Omega)}{H^{1}_{0}(\Omega)}}+a(\tilde{y}(t),\bar{p}(t))\,dt+\\ (\tilde{y}(0),\bar{q})_{{L^{2}(\Omega)}}+(\tilde{y}(T),\bar{p}(T))_{{L^{2}(\Omega)}}-(\tilde{y}(0),\bar{p}(0))_{{L^{2}(\Omega)}}\\ =\int\limits_{0}^{T}\int\limits_{\Omega}({\bar{y}}(t,x)-y_{d}(t,x))\tilde{y}(t,x)\,dxdt\quad\forall\ \tilde{y}\in W(I),$$

(2.4)

so that the function $\bar{p}$ can be identified with the unique weak solution to the adjoint equation

	$\displaystyle-\partial_{t}\bar{p}-\Delta\bar{p}$	$\displaystyle=h$		$\displaystyle\text{in }I\times\Omega\,,$		(2.5)
	$\displaystyle\bar{p}$	$\displaystyle=0$		$\displaystyle\text{on }I\times\partial\Omega\,,$
	$\displaystyle\bar{p}(T)$	$\displaystyle=0$		$\displaystyle\text{on }\Omega,$

with $h:=\bar{y}-y_{d}$. Moreover, $\bar{q}=\bar{p}(0).$

However, in order to achieve $\mathcal{O}(k^{2})$-convergence we need more regularity, i.e., at least second weak time derivatives. From [MV11, Proposition 2.1] we have

From Lemma 3.1 we conclude that the optimal state ${\bar{y}}$ lives in $H^{1}({I},{L^{2}(\Omega)})$, and ${\bar{y}}(T)\in{H^{1}_{0}(\Omega)}$. Thus, by Lemma 3.2 and Assumption 1.1 the optimal adjoint state ${\bar{p}}$ is an element of $H^{2}({I},{L^{2}(\Omega)})$. It then follows from (2.3) that $B^{\prime}{\bar{p}}\in H^{2}({I},\mathbb{R}^{D})$. Furthermore, for $v\in W^{1,r}({I},\mathbb{R}^{D})$, $1\leq r\leq\infty$, one has

$$\|\partial_{t}P_{{U_{\textup{ad}}}}(v)\|_{L^{r}({I},\mathbb{R}^{D})}\leq\|\partial_{t}v\|_{L^{r}({I},\mathbb{R}^{D})},$$

so that (2.1) and $B^{\prime}{\bar{p}}\in H^{2}({I},\mathbb{R}^{D})$ imply ${\bar{u}}\in W^{1,\infty}({I},\mathbb{R}^{D})$.

Hence, using our Assumption 1.1, Lemma 3.2 is applicable to the solution of the state equation and one obtains the following result, see e.g. [MV11, Proposition 2.3].

Let $[0,T)=\bigcup_{{m}=1}^{M}I_{m}$, where the intervals $I_{m}=[t_{{m}-1},t_{m})$ are defined through the partition $0=t_{0}<t_{1}<\dots<t_{M}=T$. Furthermore, let $t_{m}^{*}=\frac{t_{{m}-1}+t_{m}}{2}$ for $m=1,\dots,{M}$ denote the interval midpoints. By $0=:t_{0}^{*}<t_{1}^{*}<\dots<t_{M}^{*}<t_{{M}+1}^{*}:=T$ we get the so-called dual partition of $[0,T)$, namely $[0,T)=\bigcup_{{m}=1}^{{M}+1}I_{m}^{*}$, with $I_{m}^{*}=[t_{{m}-1}^{*},t_{m}^{*})$. The grid width of the first (primal) partition is defined by the mesh-parameters $k_{m}=t_{m}-t_{m-1}$ and

$$k=\max_{1\leq{m}\leq{M}}k_{m}.$$

On these partitions we define the Ansatz and test spaces of our Petrov Galerkin scheme for the numerical approximation of the optimal control problem ($\mathbb{P}$) w.r.t. time. We set

$$P_{k}:=\left\{v\in C([0,T],{H^{1}_{0}(\Omega)})\,\left|\;{\left.\kern-1.2ptv\vphantom{\big{|}}\right|_{I_{m}}}\in\mathcal{P}_{1}(I_{m},{H^{1}_{0}(\Omega)})\right.\right\}\hookrightarrow W(I),$$

$$P_{k}^{*}:=\left\{v\in C([0,T],{H^{1}_{0}(\Omega)})\,\left|\;{\left.\kern-1.2ptv\vphantom{\big{|}}\right|_{I_{m}^{*}}}\in\mathcal{P}_{1}(I_{m}^{*},{H^{1}_{0}(\Omega)})\right.\right\}\hookrightarrow W(I)$$

and

$$Y_{k}:=\left\{v:[0,T]\rightarrow{H^{1}_{0}(\Omega)}\,\left|\;{\left.\kern-1.2ptv\vphantom{\big{|}}\right|_{I_{m}}}\in\mathcal{P}_{0}(I_{m},{H^{1}_{0}(\Omega)})\right.\right\}\,.$$

Here, $\mathcal{P}_{i}(J,{H^{1}_{0}(\Omega)})$, $J\subset\bar{I}$, $i\in\{0,1\}$, denotes the set of polynomial functions in time of degree at most $i$ on the interval $J$ with values in ${H^{1}_{0}(\Omega)}$. We note that functions in $P_{k}\cup Y_{k}$ can be uniquely determined by ${M}+1$ elements from ${H^{1}_{0}(\Omega)}$. Furthermore each function in $Y_{k}$ is also an element of $L^{2}(I,{H^{1}_{0}(\Omega)})$.
In what follows we frequently use the interpolation operators

1. 1.

   $\mathcal{P}_{Y_{k}}:L^{2}(I,{H^{1}_{0}(\Omega)})\rightarrow Y_{k}$

   $${\left.\kern-1.2pt\mathcal{P}_{Y_{k}}v\vphantom{\big{|}}\right|_{I_{m}}}:=\frac{1}{k_{m}}\int_{t_{{m}-1}}^{t_{m}}vdt\text{ for }m=1,\dots,{M},\text{ and }\mathcal{P}_{Y_{k}}v(T):=0$$

2. 2.

   $\Pi_{Y_{k}}:C([0,T],{H^{1}_{0}(\Omega)})\rightarrow Y_{k}$

   $${\left.\kern-1.2pt\Pi_{Y_{k}}v\vphantom{\big{|}}\right|_{I_{m}}}:=v\left(t_{m}^{*}\right)\text{ for }m=1,\dots,{M},\text{ and }\Pi_{Y_{k}}v(T):=v(T).$$

3. 3.

   $\pi_{P_{k}^{*}}:C([0,T],{H^{1}_{0}(\Omega)})\cup Y_{k}\rightarrow P_{k}^{*}$

   	$\displaystyle{\left.\kern-1.2pt\pi_{P_{k}^{*}}v\vphantom{\big{|}}\right|_{I_{1}^{*}\cup I_{2}^{*}}}$	$\displaystyle=v(t_{1}^{*})+\frac{t-t_{1}^{*}}{t_{2}^{*}-t_{1}^{*}}(v(t_{2}^{*})-v(t_{1}^{*})),$
   	$\displaystyle{\left.\kern-1.2pt\pi_{P_{k}^{*}}v\vphantom{\big{|}}\right|_{I_{m}^{*}}}$	$\displaystyle=v\left(t_{{m}-1}^{*}\right)+\frac{t-t_{{m}-1}^{*}}{t_{m}^{*}-t_{{m}-1}^{*}}(v(t_{m}^{*})-v(t_{{m}-1}^{*})),\text{ for }m=3,\dots,{M}-1,$
   	$\displaystyle{\left.\kern-1.2pt\pi_{P_{k}^{*}}v\vphantom{\big{|}}\right|_{I_{M}^{*}\cup I_{{M}+1}^{*}}}$	$\displaystyle=v(t_{{M}-1}^{*})+\frac{t-t_{{M}-1}^{*}}{t_{M}^{*}-t_{{M}-1}^{*}}(v(t_{M}^{*})-v(t_{{M}-1}^{*})).$

To apply variational discretization to ($\mathbb{P}$) we next introduce the Petrov-Galerkin scheme for the approximation of the states. For this purpose we extend the bilinear form $A$ of (1.4) from $W(I)$ to $W(I)\cup Y_{k}$, i.e. we consider $A$ as a mapping $A:W(I)\cup Y_{k}\times W(I)\rightarrow\mathbb{R}$. Then, according to (1.5) we for $(f,\kappa)\in L^{2}(I,{H^{-1}(\Omega)})\times L^{2}(\Omega)$ consider the time-semidiscrete problem: Find $y_{k}\in Y_{k}$, such that

$$A(y_{k},v_{k})=\int\limits_{0}^{T}\langle f(t),v_{k}(t)\rangle_{{H^{-1}(\Omega)}{H^{1}_{0}(\Omega)}}\,dt+(\kappa,v_{k}(0))_{{L^{2}(\Omega)}}\quad\forall\ v_{k}\in P_{k}.$$

(4.1)

Then $y_{k}\in Y_{k}$ is uniquely determined. This follows from the fact that with

$$y_{k}=\alpha_{M+1}\chi_{\{T\}}+\sum\limits_{i=1}^{M}\alpha_{i}\chi_{I_{i}},\quad\alpha_{i}\in{H^{1}_{0}(\Omega)}\text{ for }i=1,\dots,M+1,$$

the coefficients $\alpha_{i}$ for $i=2,\dots,M$ are determined by a Crank-Nicolson scheme with a (Rannacher) smoothing step [Ran84] for $\alpha_{1}$, and $\alpha_{M+1}$ is uniquely determined by $\alpha_{M}$.

Note that in all of the following results $C$ denotes a generic, strict positive real constant that does not depend on quantities which appear to the right of it.

The following stability result will be useful in the later analysis.

Note that the proof of [MV11, Lemma 5.3] is applicable in our situation since the initial value $\kappa$ is the same for both, the continuous problem (1.2) and (4.1).

Since the state is discretized by piecewise constant functions, we can only expect first order convergence in time for its discretization error. The following Lemma shows that a projected version of the discretized state converges second order in time to the continuous state. The benefit of this result will be discussed in the numerics section.