Guaranteed lower bounds for cost functionals of time-periodic parabolic optimization problems

TL;DR

This paper introduces a new technique for deriving guaranteed lower bounds for cost functionals in time-periodic parabolic optimal control problems, enhancing the ability to estimate and control solutions accurately.

Contribution

The paper presents a novel method for obtaining computable lower bounds for two cost functionals in time-periodic parabolic control problems, complementing existing upper bounds and enabling two-sided estimates.

Findings

Derived guaranteed lower bounds for cost functionals.

Validated bounds through numerical experiments.

Established a basis for adaptive time-periodic optimization schemes.

Abstract

In this paper, a new technique is shown for deriving computable, guaranteed lower bounds of functional type (minorants) for two different cost functionals subject to a parabolic time-periodic boundary value problem. Together with previous results on upper bounds (majorants) for one of the cost functionals, both minorants and majorants lead to two-sided estimates of functional type for the optimal control problem. Both upper and lower bounds are derived for the second new cost functional subject to the same parabolic PDE-constraints, but where the target is a desired gradient. The time-periodic optimal control problems are discretized by the multiharmonic finite element method leading to large systems of linear equations having a saddle point structure. The derivation of preconditioners for the minimal residual method for the new optimization problem is discussed in more detail. Finally, several numerical experiments for both optimal control problems are presented confirming the theoretical results obtained. This work provides the basis for an adaptive scheme for time-periodic optimization problems.

Tables

Table 1. Table 1. Example 1 . Minorant 𝒥 0 ⊖ subscript superscript 𝒥 symmetric-difference 0 \mathcal{J}^{\ominus}_{0} , majorant 𝒥 0 ⊕ subscript superscript 𝒥 direct-sum 0 \mathcal{J}^{\oplus}_{0} and their efficiency indices computed on grids of different sizes.

grid	$t^{\text{sec}}$	${𝒥}_0^{\ominus }$	$I_{\text{eff}}^{{𝒥}_0^{\ominus }}$	${𝒥}_0^{\oplus }$	$I_{\text{eff}}^{{𝒥}_0^{\oplus }}$	$I_{\text{eff}}^{𝒥,0}$	$I_{\text{eff}}^{{ℳ}_1,0}$
$16\times 16$	0.02	1.13e+05	0.90	1.26e+05	1.01	1.12	1.53
$32\times 32$	0.07	1.14e+05	0.90	1.27e+05	1.00	1.11	1.47
$64\times 64$	0.24	1.14e+05	0.90	1.27e+05	1.00	1.11	1.44
$128\times 128$	1.16	1.14e+05	0.90	1.27e+05	1.00	1.11	1.43
$256\times 256$	4.51	1.14e+05	0.90	1.27e+05	1.00	1.11	1.42

Table 2. Table 2. Example 1 . Minorant 𝒥 1 ⊖ subscript superscript 𝒥 symmetric-difference 1 \mathcal{J}^{\ominus}_{1} , majorant 𝒥 1 ⊕ subscript superscript 𝒥 direct-sum 1 \mathcal{J}^{\oplus}_{1} and their efficiency indices computed on grids of different sizes.

grid	$t^{\text{sec}}$	${𝒥}_1^{\ominus }$	$I_{\text{eff}}^{{𝒥}_1^{\ominus }}$	${𝒥}_1^{\oplus }$	$I_{\text{eff}}^{{𝒥}_1^{\oplus }}$	$I_{\text{eff}}^{𝒥,1}$	$I_{\text{eff}}^{{ℳ}_1,1}$
$16\times 16$	0.02	4.27e+05	0.90	4.74e+05	1.00	1.11	6.18
$32\times 32$	0.07	4.31e+05	0.90	4.79e+05	1.00	1.11	6.01
$64\times 64$	0.25	4.32e+05	0.90	4.80e+05	1.00	1.11	5.93
$128\times 128$	1.13	4.32e+05	0.90	4.80e+05	1.00	1.11	5.90
$256\times 256$	4.66	4.32e+05	0.90	4.80e+05	1.00	1.11	5.89

Table 3. Table 3. Example 1 . Overall minorant 𝒥 ⊖ superscript 𝒥 symmetric-difference \mathcal{J}^{\ominus} and overall majorant 𝒥 ⊕ superscript 𝒥 direct-sum \mathcal{J}^{\oplus} , their parts, and their efficiency indices computed on a grid of size 256 × 256 256 256 256\times 256 .

mode	$t^{\text{sec}}$	${𝒥}^{\ominus }$	$I_{\text{eff}}^{{𝒥}^{\ominus }}$	${𝒥}^{\oplus }$	$I_{\text{eff}}^{{𝒥}^{\oplus }}$	$I_{\text{eff}}^{𝒥}$	$I_{\text{eff}}^{{ℳ}_1}$
$k=0$	4.51	1.14e+05	0.90	1.27e+05	1.00	1.11	1.42
$k=1$	4.66	4.32e+05	0.90	4.80e+05	1.00	1.11	5.89
$k=2$	4.75	1.79e+05	0.90	1.99e+05	1.00	1.11	1.64
$k=3$	4.81	6.10e+04	0.90	6.74e+04	1.00	1.11	1.84
$k=4$	4.74	7.68e+03	0.91	8.42e+03	1.00	1.10	11.06
$k=6$	4.72	2.05e+02	0.90	2.29e+02	1.00	1.11	2.08
$k=8$	4.81	1.97e+01	0.93	2.19e+01	1.04	1.12	1.22
overall ($N=3$)	–	2.86e+06	0.90	3.17e+06	1.00	1.11	2.08
overall ($N=8$)	–	2.86e+06	0.90	3.17e+06	1.00	1.11	2.09

Table 4. Table 4. Example 2 . Overall minorant 𝒥 ⊖ superscript 𝒥 symmetric-difference \mathcal{J}^{\ominus} and overall majorant 𝒥 ⊕ superscript 𝒥 direct-sum \mathcal{J}^{\oplus} , their parts, and their efficiency indices computed on a grid of size 256 × 256 256 256 256\times 256 .

mode	$t^{\text{sec}}$	${𝒥}^{\ominus }$	$I_{\text{eff}}^{{𝒥}^{\ominus }}$	${𝒥}^{\oplus }$	$I_{\text{eff}}^{{𝒥}^{\oplus }}$	$I_{\text{eff}}^{𝒥}$	$I_{\text{eff}}^{{ℳ}_1}$
$k=0$	4.55	3.20e+05	0.90	3.56e+05	1.00	1.11	1.43
$k=1$	4.53	1.02e+06	0.90	1.14e+06	1.00	1.11	3.19
$k=2$	4.62	2.57e+05	0.90	2.85e+05	1.00	1.11	3.17
$k=3$	4.80	6.07e+04	0.91	6.69e+04	1.00	1.10	1.82
$k=4$	4.75	1.99e+04	0.91	2.19e+04	1.00	1.10	1.55
$k=6$	4.83	4.05e+03	0.92	4.38e+03	1.00	1.08	1.32
$k=8$	4.90	1.30e+03	0.93	1.40e+03	1.00	1.07	1.20
$k=10$	4.72	5.40e+02	0.94	5.75e+02	1.00	1.06	1.14
overall ($N=6$)	–	6.34e+06	0.90	7.06e+06	1.00	1.11	2.09
overall ($N=10$)	–	6.34e+06	0.90	7.06e+06	1.00	1.11	2.09

Table 5. Table 5. Example 3 . Minorants, majorants, and their efficiency indices as well as the efficiency indices of the combined norm computed on a grid of size 256 × 256 256 256 256\times 256 .

mode	$t^{\text{sec}}$	${𝒥}_k^{\ominus }$	$I_{\text{eff}}^{{𝒥}_k^{\ominus }}$	${𝒥}_k^{\oplus }$	$I_{\text{eff}}^{{𝒥}_k^{\oplus }}$	$I_{\text{eff}}^{𝒥,k}$	$I_{\text{eff}}^{{ℳ}_1,k}$
$k=0$	6.44	5.71e+04	0.92	8.20e+04	1.32	1.44	6.48
$k=1$	6.45	9.30e+04	0.92	1.31e+05	1.30	1.41	3.59
$k=3$	6.46	1.06e+04	0.95	1.35e+04	1.20	1.27	2.86
$k=5$	6.39	3.90e+03	0.97	4.63e+03	1.15	1.19	3.16
$k=7$	6.46	2.01e+03	0.98	2.28e+03	1.11	1.13	3.28
$k=9$	6.55	1.23e+03	0.98	1.35e+03	1.08	1.10	3.30
$k=11$	6.48	8.23e+02	0.99	8.89e+02	1.07	1.08	3.23
$k=21$	6.44	2.28e+02	1.00	2.37e+02	1.03	1.04	3.88
$k=41$	6.46	5.99e+01	1.00	6.19e+01	1.03	1.03	4.81
$k=81$	6.42	1.53e+01	1.00	1.63e+01	1.06	1.06	3.19

Table 6. Table 6. Example 4 . Minorant 𝒥 ~ 0 ⊖ subscript superscript ~ 𝒥 symmetric-difference 0 \tilde{\mathcal{J}}^{\ominus}_{0} , majorant 𝒥 ~ 0 ⊕ subscript superscript ~ 𝒥 direct-sum 0 \tilde{\mathcal{J}}^{\oplus}_{0} and their efficiency indices computed on grids of different sizes.

grid	$t^{\text{sec}}$	${\widetilde{𝒥}}_0^{\ominus }$	$I_{\text{eff}}^{{\widetilde{𝒥}}_0^{\ominus }}$	${\widetilde{𝒥}}_0^{\oplus }$	$I_{\text{eff}}^{{\widetilde{𝒥}}_0^{\oplus }}$	$I_{\text{eff}}^{\widetilde{𝒥},0}$	$I_{\text{eff}}^{{\widetilde{ℳ}}_1,0}$
$16\times 16$	0.02	8.92e+03	0.99	9.85e+03	1.09	1.10	2.04
$32\times 32$	0.06	9.24e+03	0.99	9.95e+03	1.07	1.08	1.97
$64\times 64$	0.24	9.32e+03	0.99	9.97e+03	1.05	1.07	1.94
$128\times 128$	1.04	9.34e+03	0.98	9.98e+03	1.05	1.07	1.93
$256\times 256$	4.39	9.35e+03	0.98	9.98e+03	1.05	1.07	1.92

Table 7. Table 7. Example 4 . Minorant 𝒥 ~ 1 ⊖ subscript superscript ~ 𝒥 symmetric-difference 1 \tilde{\mathcal{J}}^{\ominus}_{1} , majorant 𝒥 ~ 1 ⊕ subscript superscript ~ 𝒥 direct-sum 1 \tilde{\mathcal{J}}^{\oplus}_{1} and their efficiency indices computed on grids of different sizes.

grid	$t^{\text{sec}}$	${\widetilde{𝒥}}_1^{\ominus }$	$I_{\text{eff}}^{{\widetilde{𝒥}}_1^{\ominus }}$	${\widetilde{𝒥}}_1^{\oplus }$	$I_{\text{eff}}^{{\widetilde{𝒥}}_1^{\oplus }}$	$I_{\text{eff}}^{\widetilde{𝒥},1}$	$I_{\text{eff}}^{{\widetilde{ℳ}}_1,1}$
$16\times 16$	0.02	3.22e+04	0.95	3.51e+04	1.03	1.09	1.52
$32\times 32$	0.06	3.39e+04	0.96	3.58e+04	1.02	1.05	1.33
$64\times 64$	0.26	3.45e+04	0.97	3.60e+04	1.01	1.04	1.25
$128\times 128$	1.08	3.48e+04	0.97	3.62e+04	1.01	1.04	1.21
$256\times 256$	4.31	3.51e+04	0.98	3.65e+04	1.01	1.04	1.19

Table 8. Table 8. Example 4 . Overall minorant 𝒥 ~ ⊖ superscript ~ 𝒥 symmetric-difference \tilde{\mathcal{J}}^{\ominus} and overall majorant 𝒥 ~ ⊕ superscript ~ 𝒥 direct-sum \tilde{\mathcal{J}}^{\oplus} , their parts, and their efficiency indices computed on a grid of size 256 × 256 256 256 256\times 256 .

mode	$t^{\text{sec}}$	${\widetilde{𝒥}}^{\ominus }$	$I_{\text{eff}}^{{\widetilde{𝒥}}^{\ominus }}$	${\widetilde{𝒥}}^{\oplus }$	$I_{\text{eff}}^{{\widetilde{𝒥}}^{\oplus }}$	$I_{\text{eff}}^{\widetilde{𝒥}}$	$I_{\text{eff}}^{{\widetilde{ℳ}}_1}$
$k=0$	4.39	9.35e+03	0.98	9.98e+03	1.05	1.07	1.92
$k=1$	4.31	3.51e+04	0.98	3.65e+04	1.01	1.04	1.19
$k=2$	4.42	9.23e+03	0.63	1.57e+04	1.06	1.70	1.65
$k=3$	4.43	2.85e+03	0.58	5.06e+03	1.03	1.78	1.08
$k=4$	4.44	5.89e+02	0.98	6.47e+02	1.07	1.10	5.91
$k=6$	4.36	1.74e+01	0.97	2.82e+01	1.58	1.62	2.32
$k=8$	4.33	2.99e-01	0.23	1.37e+00	1.05	4.59	1.97
overall ($N=6$)	–	2.09e+05	0.89	2.45e+05	1.04	1.17	2.46
overall ($N=8$)	–	2.09e+05	0.89	2.45e+05	1.04	1.17	2.46

Table 9. Table 9. Example 5 . Overall minorant 𝒥 ~ ⊖ superscript ~ 𝒥 symmetric-difference \tilde{\mathcal{J}}^{\ominus} and overall majorant 𝒥 ~ ⊕ superscript ~ 𝒥 direct-sum \tilde{\mathcal{J}}^{\oplus} , their parts, and efficiency indices of them and the combined norm computed on a grid of size 256 × 256 256 256 256\times 256 .

mode	$t^{\text{sec}}$	${\widetilde{𝒥}}^{\ominus }$	$I_{\text{eff}}^{{\widetilde{𝒥}}^{\ominus }}$	${\widetilde{𝒥}}^{\oplus }$	$I_{\text{eff}}^{{\widetilde{𝒥}}^{\oplus }}$	$I_{\text{eff}}^{\widetilde{𝒥}}$	$I_{\text{eff}}^{{\widetilde{ℳ}}_1}$
$k=0$	4.31	2.63e+04	1.00	2.79e+04	1.06	1.06	1.36
$k=1$	4.30	8.49e+04	1.00	8.60e+04	1.02	1.01	1.00
$k=2$	4.39	2.08e+04	0.98	2.21e+04	1.04	1.06	2.83
$k=4$	4.35	1.58e+03	0.96	1.75e+03	1.06	1.11	1.74
$k=6$	4.36	2.93e+02	0.87	3.53e+02	1.05	1.20	1.63
$k=8$	4.37	8.95e+01	0.82	1.20e+02	1.10	1.34	1.08
$k=10$	4.34	3.27e+01	0.71	5.23e+01	1.14	1.60	1.19
overall ($N=6$)	–	5.23e+05	1.00	5.43e+05	1.04	1.04	2.00
overall ($N=10$)	–	5.22e+05	1.00	5.42e+05	1.04	1.04	2.00

Table 10. Table 10. Example 6 . Minorants, majorants, and their efficiency indices as well as the efficiency indices of the combined norm computed on a grid of size 256 × 256 256 256 256\times 256 .

mode	$t^{\text{sec}}$	${𝒥}_k^{\ominus }$	$I_{\text{eff}}^{{𝒥}_k^{\ominus }}$	${𝒥}_k^{\oplus }$	$I_{\text{eff}}^{{𝒥}_k^{\oplus }}$	$I_{\text{eff}}^{𝒥,k}$	$I_{\text{eff}}^{{ℳ}_1,k}$
$k=0$	8.12	1.55e+01	0.79	2.05e+01	1.04	1.32	2.23
$k=1$	8.03	7.70e+00	0.86	1.13e+01	1.26	1.47	2.87
$k=3$	8.51	1.46e+01	0.93	1.73e+01	1.10	1.18	2.56
$k=5$	8.25	1.44e+01	0.99	1.65e+01	1.14	1.15	2.59
$k=7$	8.08	8.92e+00	0.98	9.80e+00	1.08	1.10	1.18
$k=9$	8.35	4.51e+00	0.90	5.14e+00	1.03	1.14	1.26
$k=11$	8.36	2.57e+00	0.96	3.13e+00	1.17	1.22	1.51
$k=21$	8.51	1.36e+00	0.99	2.27e+00	1.65	1.67	3.16
$k=41$	7.78	4.54e+00	0.86	6.26e+00	1.18	1.38	3.09
$k=81$	7.81	3.19e+00	0.79	7.09e+00	1.75	2.22	3.44

Equations

$$\sigma {\partial }_{t}y-\nabla \cdot (\nu \nabla y)$$

$$y$$

$$y(0)$$

$$H^{1,0}(Q)$$

$$H^{0,1}(Q)$$

$$H^{1,1}(Q)$$

$$\Vert u{\Vert }_{1,1}$$

$$u(\mathbf{x},t)=u_0^c(\mathbf{x})+\sum _{k=1}^{\infty }\left(u_k^c(\mathbf{x})\cos (k\omega t)+u_k^s(\mathbf{x})\sin (k\omega t)\right)$$

$$\begin{array}{r}{u}_0^c(\mathbf{x})=\frac{1}{T}{\int }_0^{T}u(\mathbf{x},t)dt,u_k^c(\mathbf{x})=\frac{2}{T}{\int }_0^{T}u(\mathbf{x},t)\cos (k\omega t)dt,u_k^s(\mathbf{x})=\frac{2}{T}{\int }_0^{T}u(\mathbf{x},t)\sin (k\omega t)dt.\end{array}$$

$$\displaystyle\big{\|}\partial^{1/2}_{t}u\big{\|}^{2}:=|u|_{0,\frac{1}{2}}^{2}:=\frac{T}{2}\sum_{k=1}^{\infty}k\omega\|\boldsymbol{u}_{k}\|_{\Omega}^{2}$$

$$|u{|}_{1,\frac{1}{2}}^2=T\Vert \nabla u_0^c{\Vert }_{\Omega }^2+\frac{T}{2}\sum _{k=1}^{\infty }\left(k\omega \Vert {\mathbf{u}}_k{\Vert }_{\Omega }^2+\Vert \nabla {\mathbf{u}}_k{\Vert }_{\Omega }^2\right)=\Vert \nabla u{\Vert }^2+\Vert {\partial }_t^{1/2}u{\Vert }^2$$

$$\begin{array}{rl}{u}^{\perp }(\mathbf{x},t)& :=\sum _{k=1}^{\infty }(-{\mathbf{u}}_k^{\perp }{)}^T\cdot (\cos (k\omega t),\sin (k\omega t){)}^T.\end{array}$$

$$\begin{array}{r}⟨\sigma {\partial }_t^{1/2}u,{\partial }_t^{1/2}v⟩=⟨\sigma {\partial }_{t}u,v^{\perp }⟩,⟨\sigma {\partial }_t^{1/2}u,{\partial }_t^{1/2}{v}^{\perp }⟩=⟨\sigma {\partial }_{t}u,v⟩,\forall u\in H_{per}^{0,1}(Q)\forall v\in H_{per}^{0,\frac{1}{2}}(Q),\end{array}$$

$$⟨\sigma {\partial }_{t}u,u⟩=0$$

$$⟨\sigma {\partial }_t^{1/2}u,{\partial }_t^{1/2}{u}^{\perp }⟩=0$$

$${\int }_Q\kappa {\partial }_t^{1/2}{u}^{\perp }d\mathbf{x}dt=-{\int }_Q{\partial }_t^{1/2}{\kappa }^{\perp }ud\mathbf{x}dt\forall \kappa ,u\in H_{per}^{0,\frac{1}{2}}(Q).$$

$$\mathcal{J}(y,u):=\frac{1}{2}\Vert y-y_d{\Vert }^2+\frac{\lambda }{2}\Vert u{\Vert }^2$$

$$\mathcal{J}(y,u)=T{\mathcal{J}}_0(y_0^c,u_0^c)+\frac{T}{2}\sum _{k=1}^{\infty }{\mathcal{J}}_k({\mathbf{y}}_k,{\mathbf{u}}_k),$$

$$\displaystyle\int_{Q}\Big{(}y\,z-\nu\nabla p\cdot\nabla z+\sigma\partial_{t}^{1/2}p\,\partial_{t}^{1/2}z^{\perp}\Big{)}\,d\boldsymbol{x}\,dt=\int_{Q}y_{d}\,z\,d\boldsymbol{x}\,dt,\qquad\forall z\in H^{1,\frac{1}{2}}_{0,per}(Q),$$

$$\displaystyle\int_{Q}\Big{(}\nu\nabla y\cdot\nabla q+\sigma\partial_{t}^{1/2}y\,\partial_{t}^{1/2}q^{\perp}+\lambda^{-1}p\,q\Big{)}\,d\boldsymbol{x}\,dt=0,\qquad\forall q\in H^{1,\frac{1}{2}}_{0,per}(Q).$$

$$\displaystyle\int_{\Omega}\big{(}\boldsymbol{y}_{k}\cdot\boldsymbol{z}_{k}-\nu\nabla\boldsymbol{p}_{k}\cdot\nabla\boldsymbol{z}_{k}+k\omega\,\sigma\boldsymbol{p}_{k}\cdot\boldsymbol{z}_{k}^{\perp}\big{)}\,d\boldsymbol{x}=\int_{\Omega}\boldsymbol{y_{d}}_{k}\cdot\boldsymbol{z}_{k}\,d\boldsymbol{x},\quad\forall\boldsymbol{z}_{k}\in\mathbb{V},$$

$$\displaystyle\int_{\Omega}\big{(}\nu\nabla\boldsymbol{y}_{k}\cdot\nabla\boldsymbol{q}_{k}+k\omega\,\sigma\boldsymbol{y}_{k}\cdot\boldsymbol{q}_{k}^{\perp}+\lambda^{-1}\boldsymbol{p}_{k}\cdot\boldsymbol{q}_{k}\big{)}\,d\boldsymbol{x}=0,\quad\forall\boldsymbol{q}_{k}\in\mathbb{V}.$$

$$\displaystyle\int_{\Omega}\big{(}y_{0}^{c}\cdot z_{0}^{c}-\nu\nabla p_{0}^{c}\cdot\nabla z_{0}^{c}\big{)}\,d\boldsymbol{x}=\int_{\Omega}{y_{d}^{c}}_{0}\cdot z_{0}^{c}\,d\boldsymbol{x},\quad\forall z_{0}^{c}\in V,$$

$$\displaystyle\int_{\Omega}\big{(}\nu\nabla y_{0}^{c}\cdot\nabla q_{0}^{c}+\lambda^{-1}p_{0}^{c}\cdot q_{0}^{c}\big{)}\,d\boldsymbol{x}=0,\quad\forall q_{0}^{c}\in V.$$

$$\mathcal{J}(y(v),v)\le {\mathcal{J}}^{\oplus }(\alpha ,\beta ;\eta ,\mathbf{\tau },v)\forall v\in L^2(Q),$$

$$\mathbf{\tau }\in H(\text{div},Q):=\{\mathbf{\tau }\in [L^2(Q){]}^d:\nabla \cdot \mathbf{\tau }(\cdot ,t)\in L^2(\Omega )\text{ for a.e. }t\in (0,T)\},$$

$${\int }_{\Omega }\nabla \cdot \mathbf{\tau }wd\mathbf{x}=-{\int }_{\Omega }\mathbf{\tau }\cdot \nabla wd\mathbf{x}\forall w\in V$$

$$\begin{array}{rl}{\mathcal{J}}^{\oplus }(\alpha ,\beta ;\eta ,\mathbf{\tau },v):=& \frac{1+\alpha }{2}\Vert \eta -y_d{\Vert }^2+\gamma (\Vert {\mathcal{R}}_2(\eta ,\mathbf{\tau }){\Vert }^2+\frac{C_F^2}{\beta }\Vert {\mathcal{R}}_1(\eta ,\mathbf{\tau },v){\Vert }^2)+\frac{\lambda }{2}\Vert v{\Vert }^2,\end{array}$$

$$|y(v)-\eta {|}_{1,\frac{1}{2}}\le \frac{1}{\underline{{\mu }_1}}\left(C_F\Vert {\mathcal{R}}_1(\eta ,\mathbf{\tau },v)\Vert +\Vert {\mathcal{R}}_2(\eta ,\mathbf{\tau })\Vert \right),$$

$${\mathcal{R}}_1(\eta ,\mathbf{\tau },v):=\sigma {\partial }_t\eta -\nabla \cdot \mathbf{\tau }-v\text{ and }{\mathcal{R}}_2(\eta ,\mathbf{\tau }):=\mathbf{\tau }-\nu \nabla \eta .$$

$$\underset{\begin{array}{c}\eta \in H_{0,per}^{1,1}(Q),\mathbf{\tau }\in H(\text{div},Q)\\ v\in L^2(Q),\alpha ,\beta >0\end{array}}{\inf }{\mathcal{J}}^{\oplus }(\alpha ,\beta ;\eta ,\mathbf{\tau },v)=\mathcal{J}(y(u),u),$$

$$\begin{array}{r}\mathcal{J}(y(u),u)\le {\mathcal{J}}^{\oplus }(\alpha ,\beta ;\eta ,\mathbf{\tau },v)\forall \eta \in H_{0,per}^{1,1}(Q),\mathbf{\tau }\in H(\text{div},Q),v\in L^2(Q),\alpha ,\beta >0.\end{array}$$

$$\tilde{\mathcal{J}}(y,u):=\frac{1}{2}\Vert \nabla y-{\mathbf{g}}_d{\Vert }^2+\frac{\lambda }{2}\Vert u{\Vert }^2$$

$$\tilde{\mathcal{L}}(y,u,p)=\tilde{\mathcal{J}}(y,u)-{\int }_Q\left(\sigma {\partial }_{t}y-\nabla \cdot (\nu \nabla y)-u\right)pd\mathbf{x}dt.$$

Full text

This research was supported by the Academy of Finland, grant 295897.

1. Introduction

In this work, we derive fully computable, guaranteed lower bounds (minorants) for cost functionals of parabolic optimization problems with given time-periodic conditions. The fully computable upper bounds (majorants) for one of the cost functionals were derived in [26]. The second cost functional is new and for this one both upper and lower bounds are presented. The motivation for the second problem lies in applications, where the target is the gradient or flux instead of the state function. Lower bounds for cost functionals of time-periodic parabolic optimal control problems have not been discussed yet. However, optimal control problems are highly important for different applications (see e.g. [36] as well as the original work [29]). These applications also include time-periodic problems, see for instance [1] and [15] considering problems in electromagnetics and biochemistry, respectively. For these types of problems the multiharmonic (or harmonic-balanced) finite element method (short MhFEM) is a natural choice. The functions are approximated by truncated Fourier series in time and by the finite element method (FEM) in space – more precisely, the Fourier coefficients. We refer to the application of this discretization technique already in [45] as well as later in [3, 4, 5, 9] for non-linear time-harmonic eddy current problems. Moreover, time-periodic optimal control problems and the MhFEM were discussed in e.g. [16, 17, 27] and [18, 19]. Recent works on robust preconditioners for time-periodic parabolic and eddy-current optimal control problems are discussed in [28] and [2], respectively. In this work, a standard finite element discretization is used with continuous, piecewise linear finite elements and a regular grid as discussed, e.g., in [6, 7, 43]. However, the method is wider applicable using also other finite elements or also finite differences (for instances, if the given domain is geometrically rather non-complex).

Functional a posteriori estimation provides a useful machinery to derive computable and guaranteed quantities for the desired unknown solution, see, e.g., [39, 12] on parabolic problems. Recent works on new estimates for parabolic problems and parabolic optimal control problems can be found in [32] and [42], respectively. A posteriori estimates of functional type for elliptic optimal control problems can be found in [10, 11, 40, 31]. First functional type estimates for inverse problems, which are related to optimal control problems, can be found in [41, 8]. Moreover, recent results on guaranteed computable estimates for convection-dominated diffusion problems are presented in [33]. In [26], majorants for one cost functional of a time-periodic parabolic optimal control and for the corresponding optimality system were presented. This work presents the corresponding minorants for this cost functional using the new technique presented in [44], which makes use of ideas derived by Mikhlin [34] but generalized for the class of optimal control problems. We mention here that [30] presents a different approach for the derivation of a lower bound for a class of elliptic optimal control problems.

We extend the analysis in this paper and consider a second cost functional with respect to the same parabolic time-periodic boundary value problem. In the second optimal control problem, the target is a given desired gradient. Problems of that type have been earlier discussed in [11]. The results on computable lower bounds together with the upper bounds lead to two-sided estimates which can be used to derive majorants for the discretization error in state and control. These majorants and minorants provide a new formulation of the optimization problems since they can, in principle, be used as objects of direct minimization on their difference. The majorants and minorants can be used in order to derive an adaptive scheme in time and in space. The linearity of the optimal control problems leads to decoupled problems in the Fourier modes including decoupling in the majorants and minorants, which is shown in this work. The overall estimators provide the modes/mode numbers which are necessary for computations. The problems for the different modes can then be computed on different grids, for which estimators in space can be used with respect to the finite element discretized Fourier coefficients. Altogether we derive a space-time adaptive method. Its idea has been for the first time introduced in [26] and has been called the adaptive multiharmonic finite element method (AMhFEM).

In this work, robust preconditioners for the preconditioned minimal residual (MinRes) method (see [37]) are discussed for the second optimization problem, which are new for this case. Also the practical performance of the AMLI preconditioner MinRes solver is presented in various numerical experiments for both optimization problems. For additional numerical tests regarding the AMLI precondititioned MinRes solver used in this work and its performance for different cases of given data in time-periodic parabolic optimal control problems, we refer the reader to [27].

The work is arranged in the following sections: In Section 2, the two types of cost functionals are presented including some preliminary results. We denote the problems by optimization problem I and II. Moreover, the former result on the majorant for problem I is summarized there. The new minorant for optimization problem I is derived in Section 3 followed by the discussion of the majorant and minorant for optimization problem II in Section 4. In Section 5, the MhFEM for both optimization problems is presented. Robust preconditioners for applying the preconditioned MinRes method on the problems discretized by the MhFEM are presented in Section 6. Section 7 discusses detailed a set of various numerical experiments for both optimization problems I and II. A few final remarks are drawn in Section 8.

2. Time-periodic parabolic PDE, the two cost functionals and preliminary results

We denote by $\Omega\subset\mathbb{R}^{d}$ with possible dimensions $d=\{1,2,3\}$ the spatial bounded Lipschitz domain with boundary $\Gamma:=\partial\Omega$. Also we denote by $Q:=\Omega\times(0,T)$ the space-time domain and $\Sigma:=\Gamma\times(0,T)$ its lateral surface, where $(0,T)$ is the given time interval. The optimization problems are both subject to the following parabolic PDE with given homogeneous Dirichlet boundary conditions and time-periodical condition:

[TABLE]

The function $y$ is the state and $u$ will be the control function. We assume that the coefficient functions $\sigma$ and $\nu$ are positive and bounded satisfying $0<\underline{\sigma}\leq\sigma(\boldsymbol{x})\leq\overline{\sigma}$ and $0<\underline{\nu}\leq\nu(\boldsymbol{x})\leq\overline{\nu}$ for $\boldsymbol{x}\in\Omega$ with constants $\underline{\sigma}$, $\overline{\sigma}$, $\underline{\nu}$ and $\overline{\nu}$. The time-periodic problems in this paper are motived by real-life applications such as computational electromagnetics, where these parameters correspond to the reluctivity and conductivity being usually piecewise constant because of the various materials of the electrical devices.

2.1. Preliminaries

In the following, we present a proper functional space setting for time-periodic problems which starts by defining the Hilbert spaces

[TABLE]

(see also [23, 24]). For instance, the norm in $H^{1,1}$ is given by

[TABLE]

In the following, we skip the subindex for the norms and inner products in $L^{2}(Q)$ and denote them by $\|\cdot\|$ and $\langle\cdot,\cdot\rangle$. For $L^{2}(\Omega)$ and $H^{1}(\Omega)$, we denote them by $\|\cdot\|_{\Omega}$ and $\langle\cdot,\cdot\rangle_{\Omega}$ and $\|\cdot\|_{1,\Omega}$ and $\langle\cdot,\cdot\rangle_{1,\Omega}$, respectively.

Time-periodic functions which are at least from the space $L^{2}$ can be naturally represented by Fourier series as

[TABLE]

for $\omega=2\pi/T$ being the frequency, $T$ the period and with the Fourier cofficients

[TABLE]

We define the norm in Fourier space as follows

[TABLE]

as well as the spaces H^{0,\frac{1}{2}}_{per}(Q)=\{u\in L^{2}(Q):\big{\|}\partial^{1/2}_{t}u\big{\|}<\infty\}, H^{1,\frac{1}{2}}_{per}(Q)=\{u\in H^{1,0}(Q):\big{\|}\partial^{1/2}_{t}u\big{\|}<\infty\}, $H^{1,\frac{1}{2}}_{0,per}(Q)=\{u\in H^{1,\frac{1}{2}}_{per}(Q):u=0\mbox{ on }\Sigma\}$, where $\boldsymbol{u}_{k}=(u_{k}^{c},u_{k}^{s})^{T}$, $k\in\mathbb{N}$. We also introduce the orthogonal vector $\boldsymbol{u}_{k}^{\perp}=(-u_{k}^{s},u_{k}^{c})^{T}$. The inner products (including also a $\sigma$-weighted version) in these spaces are defined by $\langle\partial^{1/2}_{t}u,\partial^{1/2}_{t}v\rangle:=\frac{T}{2}\sum_{k=1}^{\infty}k\omega\langle\boldsymbol{u}_{k},\boldsymbol{v}_{k}\rangle_{\Omega}$, and $\langle\sigma\partial^{1/2}_{t}u,\partial^{1/2}_{t}v\rangle:=\frac{T}{2}\sum_{k=1}^{\infty}k\omega\langle\sigma\boldsymbol{u}_{k},\boldsymbol{v}_{k}\rangle_{\Omega}$. The $H^{1,\frac{1}{2}}_{per}(Q)$-seminorm is defined by

[TABLE]

and the corresponding norm is $\|u\|_{1,\frac{1}{2}}^{2}:=\|u\|^{2}+|u|_{1,\frac{1}{2}}^{2}$. Finally, we also define the product $\langle\kappa,\partial_{t}^{1/2}u\rangle:=\frac{T}{2}\sum_{k=1}^{\infty}(k\omega)^{1/2}\langle\boldsymbol{\kappa}_{k},\boldsymbol{u}_{k}\rangle_{\Omega}$ as well as the orthogonal Fourier series

[TABLE]

Using this notation we can prove that $(u^{\perp})^{\perp}$ = -u, \big{\|}u^{\perp}\big{\|}=\big{\|}u\big{\|} and \big{\|}\partial^{1/2}_{t}u^{\perp}\big{\|}=\big{\|}\partial^{1/2}_{t}u\big{\|} for all $u\in H^{0,\frac{1}{2}}_{per}(Q)$ and also that $\|\boldsymbol{u}_{k}^{\perp}\|_{\Omega}=\|\boldsymbol{u}_{k}\|_{\Omega}$. Briefly, we recall from [25], the following identities

[TABLE]

as well as orthogonality relations

[TABLE]

Note that the following identity is valid (in Fourier series sense)

[TABLE]

Friedrichs’ inequality in $Q$ can be proved by using standard Friedrichs’ inequality on the Fourier coefficients with respect to the spatial domain $\Omega$. We have that $\|\nabla u\|^{2}\geq\frac{1}{C_{F}^{2}}\|u\|^{2}$, where $C_{F}>0$ is Friedrichs’ constant.

In the following, the parameter $\lambda>0$ denotes the regularization or cost parameter.

2.2. Optimization problem I

In the first case, we want to minimize the following cost functional with respect to the unknown state $y$ and control $u$:

[TABLE]

subject to the time-periodic boundary value problem (1)–(3). The given desired state $y_{d}\in L^{2}(Q)$ does not have to be time-periodic. It only has to be from the space $L^{2}(Q)$. The cost functional $\mathcal{J}$ defined in (9) can be written as

[TABLE]

where $\mathcal{J}_{0}(y_{0}^{c},u_{0}^{c})=\frac{1}{2}\|y_{0}^{c}-{y_{d}^{c}}_{0}\|_{\Omega}^{2}+\frac{\lambda}{2}\|u_{0}^{c}\|_{\Omega}^{2}$ and $\mathcal{J}_{k}(\boldsymbol{y}_{k},\boldsymbol{u}_{k})=\frac{1}{2}\|\boldsymbol{y}_{k}-{\boldsymbol{y}_{d}}_{k}\|_{\Omega}^{2}+\frac{\lambda}{2}\|\boldsymbol{u}_{k}\|_{\Omega}^{2}$. In [26], the corresponding optimality system is derived, which is given in weak formulation as follows: Given $y_{d}\in L^{2}(Q)$, find $y,p\in H^{1,\frac{1}{2}}_{0,per}(Q)$ such that

[TABLE]

The reduced optimality system (10)–(11) i.a. results from using the condition $u=-\lambda^{-1}p$, since no box constraints are imposed on the control function $u$ in this paper. This also leads to the space of admissible controls being given by $H^{1,\frac{1}{2}}_{0,per}(Q)$.

Remark 1**.**

Since the optimality system is derived first, and then the discretization is performed later, we can say that the first optimize, then discretize approach is applied here (as discussed, e.g., earlier in [14]).

Let us define $V:=H^{1}_{0}(\Omega)$ and $\mathbb{V}:=V\times V$. Expanding all functions into Fourier series in (10)–(11) together with using the orthogonality of the cosine and sine functions yields the following problems corresponding to the $k$th and 0th Fourier coefficients and which are decoupled due to the linearity of the optimal control problem: Find $\boldsymbol{y}_{k},\boldsymbol{p}_{k}\in\mathbb{V}$ so that

[TABLE]

and $y_{0}^{c},p_{0}^{c}\in V$ so that

[TABLE]

Both problems (12)–(13) and (14)–(15) are uniquely solvable (see [27]).

2.3. Majorant for

cost functional (9)

Here, the results of [26] on upper bounds for optimization problem I are summerized, since they are needed later to derive the two-sided estimate for the cost functional (9), which deepens and extends the a posteriori error analysis for optimization problem I. Let $y=y(v)$ be the corresponding state to an arbitrary control $v$. The following upper bound can be proved:

[TABLE]

for arbitrary $\alpha,\beta>0$, $\eta\in H^{1,1}_{0,per}(Q)$ and

[TABLE]

where, for any $\boldsymbol{\tau}\in H(\text{div},Q)$, the identity

[TABLE]

is valid. The guaranteed and fully computable majorant is given by

[TABLE]

where $\underline{\mu_{1}}=\frac{1}{\sqrt{2}}\min\{\underline{\nu},\underline{\sigma}\}$ and we have set $\gamma:=\frac{(1+\alpha)(1+\beta)C_{F}^{2}}{2\alpha\underline{\mu_{1}}^{2}}$. The parameters $\alpha,\beta>0$ have been introduced in order to obtain a quadratic functional by applying Young’s inequality.

Remark 2**.**

The arbitrary functions $\eta\in H^{1,1}_{0,per}(Q)$ and $v\in L^{2}(Q)$ can be taken later as the approximate solutions of the optimal control problem (9) subject to (1)–(3) and $\boldsymbol{\tau}\in H(\emph{div},Q)$ represents the image of the flux $\nu\nabla\eta$. Note again that first optimize, then discretize is applied in this paper, see also Remark 1.

For the derivation of (16), the following estimate for the approximation error has been used:

[TABLE]

where

[TABLE]

The derivation of estimate (17) can be found in [25]. The function $\mathcal{J}^{\oplus}(\alpha,\beta;\eta,\boldsymbol{\tau},v)$ is a sharp upper bound on $\mathcal{J}(y(v),v)$ for arbitrary but fixed $v$ as well as on the optimal value $\mathcal{J}(y(u),u)$

[TABLE]

since the infimum is attained for the optimal control $u$, its corresponding state $y(u)$ and its exact flux $\nu\nabla y(u)$, and for $\alpha$ going to zero. Therefore, we have the estimate

[TABLE]

2.4. Optimization problem II

In the second case, we want to minimize the following cost functional with respect to the unknown state $y$ and control $u$:

[TABLE]

subject to the time-periodic boundary value problem (1)–(3), where $\boldsymbol{g}_{d}\in[L^{2}(Q)]^{d}$ is the given desired gradient. The optimality system can analogously be derived as for optimization problem I using the Lagrange functional

[TABLE]

Optimality equations $\partial_{u}\tilde{\mathcal{L}}(y,u,p)=0$ and $\partial_{p}\tilde{\mathcal{L}}(y,u,p)=0$ are similar. Equation $\partial_{y}\tilde{\mathcal{L}}(y,u,p)=0$ is different. The optimality conditions are given in weak form as follows: Given $\boldsymbol{g}_{d}\in[L^{2}(Q)]^{d}$, find $y,p\in H^{1,\frac{1}{2}}_{0,per}(Q)$ such that

[TABLE]

The optimality systems corresponding to every mode $k$ are analogously derived as for optimization problem I (similar to (12)–(13) and (14)–(15)). In Section 4, we will derive new two-sided bounds for optimization problem II.

3. Guaranteed lower bounds leading to

two-sided bounds for optimization problem I

In this work, we complement the guaranteed upper bounds for the discretization error in state and control of minimizing cost functional $\mathcal{J}$ defined in (9) subject to (1)–(3). This is done by obtaining fully computable lower bounds for $\mathcal{J}$ following the technique from [44] (derived for elliptic problems) leading to two-sided bounds for the cost functional (9).

3.1. Minorant for cost functional (9)

Let $y=y(u)$ be the optimal state corresponding to the optimal control function $u\in L^{2}(Q)$, which is connected with the adjoint state $p=p(u)$ by the identity $u=-\lambda^{-1}p(u)$. Then $y=y(u)$ is the solution of the variational formulation

[TABLE]

of the boundary value problem (1)–(3) (see also (10)–(11)). For any arbitray function $\eta\in H^{1,1}_{0,per}(Q)$, one can obtain that

[TABLE]

Since $\frac{1}{2}\|y-\eta\|^{2}\geq 0$ and using the identity $u=-\lambda^{-1}p(u)$, $\mathcal{J}$ can be estimated from below by

[TABLE]

For $\eta\in H^{1,1}_{0,per}(Q)$, let $p_{\eta},\tilde{p}_{\eta}\in H^{1,\frac{1}{2}}_{0,per}(Q)$ be the solutions to the equations

[TABLE]

Remark 3**.**

Note that we assumed that $\eta\in H^{1,1}_{0,per}(Q)$ according to the derivation of the majorant, but so far the assumption $\eta\in H^{1,\frac{1}{2}}_{0,per}(Q)$ would be enough.

Adding and subtracting $p_{\eta}$ in (25) together with $\frac{1}{2\lambda}\|p-p_{\eta}\|^{2}\geq 0$ yields the estimate

[TABLE]

Using equation (26) and identity (8) leads to the estimate

[TABLE]

By applying equations (24) and (27), it follows that

[TABLE]

We introduce now the arbitrary function $\zeta\in H^{1,1}_{0,per}(Q)$. Note that at the moment $\zeta\in H^{1,\frac{1}{2}}_{0,per}(Q)$ would be enough, but the higher regularity in time will be needed in another step. This goes along with the higher regularity assumption on $\eta$ (see Remark 3). Since $\frac{1}{2\lambda}\|p_{\eta}-\zeta\|^{2}\geq 0$, we have that

[TABLE]

Now we add and subtract $\lambda^{-1}\tilde{p}_{\eta}\zeta$ in the last integral as well as use equation (27) again. Moreover, we exploit the fact that we have assumed that $\eta\in H^{1,1}_{0,per}(Q)$, hence, we can apply the identities (5). Altogether these steps yield the estimate

[TABLE]

In the following, we need to estimate the last integral of this expression in order to formulate a computable lower bound for the cost functional. For that let us first prove a computable upper bound for the error in the adjoint state, which is presented in the following theorem. Note that here we will need the higher regularity assumption (in time) on $\zeta$.

Theorem 1**.**

Let $y_{d}\in L^{2}(Q)$ be given. Let $p_{\eta}\in H^{1,\frac{1}{2}}_{0,per}(Q)$ solve (26) for an arbitrary $\eta\in H^{1,1}_{0,per}(Q)$. The following estimate holds:

[TABLE]

for any $\zeta\in H^{1,1}_{0,per}(Q)$, where $\underline{\mu_{1}}=\frac{1}{\sqrt{2}}\min\{\underline{\nu},\underline{\sigma}\}$, $\mathcal{R}_{3}(\zeta,\boldsymbol{\rho},\eta)=\eta-y_{d}+\nabla\cdot\boldsymbol{\rho}+\sigma\partial_{t}\zeta$ and $\mathcal{R}_{4}(\zeta,\boldsymbol{\rho})=\boldsymbol{\rho}-\nu\nabla\zeta$ with $\boldsymbol{\rho}\in H(\emph{div},Q)$ and $C_{F}>0$.

Proof.

Let us consider the adjoint equation (26). Adding and subtracting $\zeta\in H^{1,1}_{0,per}(Q)$ in the left-hand side of the equation leads to

[TABLE]

Next we introduce the auxiliary variable $\boldsymbol{\rho}\in H(\text{div},Q)$. Together with using that $\zeta\in H^{1,1}_{0,per}(Q)$ as well as applying Cauchy–Schwarz’ and Friedrichs’ inequalities, the following estimate for the right-hand side of (31) can be obtained:

[TABLE]

Using the boundedness of the coefficients $\sigma$ and $\nu$, the orthogonality relations (7) and applying that $\|z^{\perp}\|=\|z\|$, we can prove the estimate from below for the left-hand side of (31), which is

[TABLE]

First, we estimate the supremum from below with choosing $z=(p_{\eta}-\zeta)+(p_{\eta}-\zeta)^{\perp}$, for which

[TABLE]

using the orthogonality relations (7). Next, applying the second equation in (7) gives

[TABLE]

and applying the first equation in (7) as well as using the identity $((p_{\eta}-\zeta)^{\perp})^{\perp}=-(p_{\eta}-\zeta)$ gives

[TABLE]

leading to the estimate

[TABLE]

where $\underline{\mu_{1}}=\frac{1}{\sqrt{2}}\min\{\underline{\nu},\underline{\sigma}\}$. Combining now both estimates together with $\|\nabla(p_{\eta}-\zeta)\|\leq|p_{\eta}-\zeta|_{1,\frac{1}{2}}$ we finally derive (30). ∎

Now we have all the tools in order to estimate the last term of (29). We obtain as follows

[TABLE]

leading to

[TABLE]

where $\boldsymbol{\tau},\boldsymbol{\rho}\in H(\text{div},Q)$ and we have used equation (27), relations (6)–(7), Cauchy–Schwarz’ and Friedrichs’ inequalities, estimate (30) and that $\eta\in H^{1,1}_{0,per}(Q)$.

Finally, we obtain the following estimate from below:

[TABLE]

where the (fully computable) minorant is given by

[TABLE]

Theorem 2**.**

The supremum of the minorant $\mathcal{J}^{\ominus}$ as defined in (34) is attained for the minimum of the cost functional (9) subject to (1)–(3), which is equivalent to the optimal value of the optimality system (10)–(11) as follows

[TABLE]

Proof.

The estimate is sharp for the exact control $u$, $\eta=y(u)$, $\zeta=p(u)$, $\boldsymbol{\tau}=\nu\nabla y(u)$ and $\boldsymbol{\rho}=\nu\nabla p(u)$. Hence, $\mathcal{J}^{\ominus}(y(u),p(u),\nu\nabla y(u)),\nu\nabla p(u))=\frac{1}{2}\|y-y_{d}\|^{2}+\frac{1}{2\lambda}\|p\|^{2}=\frac{1}{2}\|y-y_{d}\|^{2}+\frac{\lambda}{2}\|u\|^{2}=\mathcal{J}(y(u),u)$. ∎

Remark 4**.**

It has been shown in [10] that if we choose finite dimensional subspaces that are limit dense in the spaces of the exact solution $y(u)$, of its adjoint $p(u)$ and of their fluxes, $\nu\nabla y(u)$ and $\nu\nabla p(u)$, which are $H^{1,1}_{0,per}(Q)$ and $H(\emph{div},Q)$, and we choose sequences of the functions $(\eta,{\zeta},\boldsymbol{\tau},\boldsymbol{\rho})$ in these finite dimensional subspaces, so for instance $(\eta_{i},{\zeta}_{i},\boldsymbol{\tau}_{i},\boldsymbol{\rho}_{i})$, then they converge (let $i\rightarrow\infty$) to the exact solution, its adjoint and their fluxes. The corresponding majorants $\mathcal{J}^{\oplus}_{i}$ and minorants $\mathcal{J}^{\ominus}_{i}$ converge to the exact value of the cost functional $\mathcal{J}$. As shown in [10], majorants can be used in order to produce sequences of state and controls with values of the cost functional being as close to the exact cost functional value as one needs it.

3.2. Guaranteed upper bounds for the discretization errors of the control

and the state

Here, we present a posteriori error estimates for control and state measured by the norm $|||u-v|||^{2}:=\frac{1}{2}\|y(u)-y(v)\|^{2}+\frac{\lambda}{2}\|u-v\|^{2}$. The next theorem was proved for the elliptic case (together with control constraints) in [44]. The norm $|||\cdot|||$ can be represented in terms of the state and the adjoint state (instead of the control), since there are no control constraints imposed. Hence, $u=-\lambda^{-1}p(u)$, $v=-\lambda^{-1}p(v)$, and $|||u-v|||^{2}=\frac{1}{2}\|y(u)-y(v)\|^{2}+\frac{1}{2\lambda}\|p(u)-p(v)\|^{2}$.

Theorem 3**.**

We obtain the following identity for an arbitrary control $v\in L^{2}(Q)$:

[TABLE]

Proof.

Together with $u=-\lambda^{-1}p(u)$ and $v=-\lambda^{-1}p(v)$ the difference can be computed as

[TABLE]

Since the adjoint states $p(u),p(v)\in H^{1,1}_{0,per}(Q)$ fulfill (10)–(11) for the corresponding states $y(u)$, $y(v)\in H^{1,1}_{0,per}(Q)$, we obtain

[TABLE]

This proves now the equality relation (36) by applying the equations (10)–(11) for $(u,y(u),p(u))$ as well as $(v,y(v),p(v))$. ∎

Using the result of Theorem 3, we can derive the majorant for the discretization errors of control and state in the norm $|||\cdot|||$.

Theorem 4**.**

The functional

[TABLE]

for an arbitrary $v\in L^{2}(Q)$, $\eta,\zeta\in H^{1,1}_{0,per}(Q)$, $\boldsymbol{\tau},\boldsymbol{\rho}\in H(\emph{div},Q)$ and $\alpha,\beta>0$, and where $\underline{\mu_{1}}=\frac{1}{\sqrt{2}}\min\{\underline{\nu},\underline{\sigma}\}$, is a majorant for the discretization error

[TABLE]

Proof.

Applying (36) together with (20) and (33) yields estimate (37). ∎

Proposition 1**.**

The infimum of the majorant $\mathcal{M}^{\oplus}(\alpha,\beta;\eta,\zeta,\boldsymbol{\tau},\boldsymbol{\rho},v)$ in (37) is attained for the minimum of the optimization problem I, which is equivalent to the solution of the optimality system (10)–(11) $(v=u$, $\eta=y(u)$, ${\zeta}=p(u)=-\lambda u$, $\boldsymbol{\tau}=\nu\nabla y(u)$, $\boldsymbol{\rho}=\nu\nabla p(u))$ as follows

[TABLE]

Proof.

The majorant $\mathcal{M}^{\oplus}(\alpha,\beta;y(u),p(u),\nu\nabla y(u),\nu\nabla p(u),u)=\,\frac{\alpha}{2}\|y(u)-y_{d}\|^{2}$ equals zero if $\alpha$ goes to zero. ∎

The majorant $\mathcal{M}^{\oplus}$ is a sharp, guaranteed and fully computable upper bound for the control-state error in $|||\cdot|||$. However, it overestimates the $L^{2}$-norm $|||\cdot|||$ which is of order $h^{2}$, since the majorant only decreases with order $h$. Following the idea from [44] a weighted $H^{1}$-norm is introduced decreasing with the same order as the majorant. For that, we define the norm $|||u-v|||_{1}^{2}:=\frac{1}{2}\|y(u)-y(v)\|^{2}+\frac{\lambda\underline{\mu_{1}}^{2}}{2C_{F}^{2}}|y(u)-y(v)|_{1,\frac{1}{2}}^{2}$.

Theorem 5**.**

The following estimate:

[TABLE]

is valid for an arbitrary control function $v\in L^{2}(Q)$ and with $\mathcal{R}_{1}(\eta,\boldsymbol{\tau},v)$ and $\mathcal{R}_{2}(\eta,\boldsymbol{\tau})$ defined as in (18).

Proof.

Let the parameter $\delta>0$ be arbitrary but fixed. We add and subtract $\eta$, apply triangle inequality on $\frac{\underline{\mu_{1}}^{2}}{C_{F}^{2}\delta}|y(u)-y(v)|_{1,\frac{1}{2}}^{2}$ and obtain

[TABLE]

where we further add and subtract $v$. Then we apply triangle inequality two times leading to

[TABLE]

Together with (36) this yields

[TABLE]

We see that the choice $\delta=2/\lambda$ finally provides estimate (38). ∎

This theorem directly leads to the following two results presented in Propositions 2 and 3.

Proposition 2**.**

The following error majorant for any control $v\in L^{2}(Q)$ is obtained:

[TABLE]

with

[TABLE]

where $\underline{\mu_{1}}=\frac{1}{\sqrt{2}}\min\{\underline{\nu},\underline{\sigma}\}$, $\alpha,\beta>0$ as well as arbitrary $\eta,\zeta\in H^{1,1}_{0,per}(Q)$ and $\boldsymbol{\tau},\boldsymbol{\rho}\in H(\emph{div},Q)$.

Proposition 3**.**

The infimum of the majorant $\mathcal{M}_{1}^{\oplus}(\alpha,\beta;\eta,\zeta,\boldsymbol{\tau},\boldsymbol{\rho},v)$ in (39) is attained for the minimum of the optimization problem I being equivalent to solving the optimality system (10)–(11) $(v=u$, $\eta=y(u)$, ${\zeta}=p(u)=-\lambda u$, $\boldsymbol{\tau}=\nu\nabla y(u)$, $\boldsymbol{\rho}=\nu\nabla p(u))$ as follows

[TABLE]

and $\alpha$ going to zero.

4. Two-sided bounds for optimization problem II

In this section, we analogously derive the majorants and minorants for the second cost functional, however, skipping details which are similar to the derivation in the case of optimization problem I.

4.1. Majorant for cost functional (21)

Adding and subtracting $\nabla\eta$ in the cost functional $\tilde{\mathcal{J}}(y(v),v)$, applying the triangle inequality as well as using the estimate

[TABLE]

we conclude that

[TABLE]

Together with (17) this leads to the estimate

[TABLE]

where again $\underline{\mu_{1}}=\frac{1}{\sqrt{2}}\min\{\underline{\nu},\underline{\sigma}\}$ as well as $\mathcal{R}_{1}(\eta,\boldsymbol{\tau},v)$ and $\mathcal{R}_{2}(\eta,\boldsymbol{\tau})$ are defined as in (18). Finally, introducing parameters $\alpha,\beta>0$ and applying Young’s inequality, we can reformulate the estimate such that the right-hand side is given by a quadratic functional as follows

[TABLE]

where

[TABLE]

The infimum of the majorant (40) is attained for the minimum of the cost functional (21) subject to (1)–(3), which is equivalent to the optimal value of the optimality system (22)–(23). Analogously to (19), we can show that

[TABLE]

and that

[TABLE]

4.2. Minorant for cost functional (21)

Let us derive now the minorant. For any $\eta\in H^{1,1}_{0,per}(Q)$, we have that

[TABLE]

for all $v\in L^{2}(Q)$. The first norm is again greater or equal to zero, together with the identity $v=-\lambda^{-1}p(v)$, we can estimate $\tilde{\mathcal{J}}$ from below by

[TABLE]

Note that Remark 3 applies here as well. For $\eta\in H^{1,1}_{0,per}(Q)$, let $p_{\eta},\tilde{p}_{\eta}\in H^{1,\frac{1}{2}}_{0,per}(Q)$ be the solutions to the equations

[TABLE]

Deriving the minorant for the second minimization functional uses similar ideas now as presented for problem I (see Subsection 3.1). However, in the following we present the main steps which are important for problem II. Adding and subtracting $p_{\eta}$ together with $\frac{1}{2\lambda}\|p-p_{\eta}\|^{2}\geq 0$ yields

[TABLE]

Appying equation (43), identity (8) (analogously to (28)) and then using equations (24) and (44) provides the estimate

[TABLE]

Together with introducing an arbitrary function $\zeta\in H^{1,1}_{0,per}(Q)$, following (29), applying Theorem 1 and using (32), we finally derive the estimate

[TABLE]

with the fully computable minorant

[TABLE]

The following theorem is analogous to Theorem 2 and so is its proof. Also Remark 4 can be applied for optimization problem II.

Theorem 6**.**

The supremum of the minorant $\tilde{\mathcal{J}}^{\ominus}$ as defined in (46) is attained for the minimum of the cost functional (21) subject to (1)–(3), which is equivalent to the optimal value of the optimality system (22)–(23) as follows

[TABLE]

4.3. Guaranteed upper bounds for the discretization errors of the control

and the state

We present the majorants for the control-state errors measured in the norm $||||u-v||||^{2}:=\frac{1}{2}\|\nabla y(u)-\nabla y(v)\|^{2}+\frac{\lambda}{2}\|u-v\|^{2}=\frac{1}{2}\|\nabla y(u)-\nabla y(v)\|^{2}+\frac{1}{2\lambda}\|p(u)-p(v)\|^{2}$. We obtain the identity $||||u-v||||^{2}=\tilde{\mathcal{J}}(y(v),v)-\tilde{\mathcal{J}}(y(u),u)$ for an arbitrary control $v\in L^{2}(Q)$ yielding the error majorant

[TABLE]

for $\alpha,\beta>0$, arbitrary $\eta,\zeta\in H^{1,1}_{0,per}(Q)$, $\boldsymbol{\tau},\boldsymbol{\rho}\in H(\text{div},Q)$ with

[TABLE]

where $\underline{\mu_{1}}=\frac{1}{\sqrt{2}}\min\{\underline{\nu},\underline{\sigma}\}$. The infimum of this majorant is attained for the minimum of the optimization problem II being equivalent to solving the optimality system (22)–(23) $(v=u$, $\eta=y(u)$, ${\zeta}=p(u)=-\lambda u$, $\boldsymbol{\tau}=\nu\nabla y(u)$ and $\boldsymbol{\rho}=\nu\nabla p(u))$ as follows

[TABLE]

Analogously, defining $||||u-v||||_{1}^{2}:=\frac{1}{2}\|\nabla y(u)-\nabla y(v)\|^{2}+\frac{\lambda\underline{\mu_{1}}^{2}}{2C_{F}^{2}}|y(u)-y(v)|_{1,\frac{1}{2}}^{2}$ we derive

[TABLE]

where now

[TABLE]

All other results similar to Propositions 2 and 3 follow completely.

5. Multiharmonic finite element (MhFE) discretization

The desired state $y_{d}$ and desired gradient $\boldsymbol{g}_{d}$ belong to $L^{2}(Q)$ and $[L^{2}(Q)]^{d}$, respectively. So they can be represented as Fourier series having Fourier coefficients from $L^{2}(\Omega)$. Moreover, we assume that $\eta$ and $\zeta$ approximating the exact state $y$ and adjoint state $p$, respectively, as well as the vector-valued functions $\boldsymbol{\tau}$ and $\boldsymbol{\rho}$ are given by truncated Fourier series. We also have the multiharmonic time derivative defined by $\partial_{t}\eta(\boldsymbol{x},t)=\sum_{k=1}^{N}\left(k\omega\,\eta_{k}^{s}(\boldsymbol{x})\cos(k\omega t)-k\omega\,\eta_{k}^{c}(\boldsymbol{x})\sin(k\omega t)\right)$ as well as the gradient and divergence by

[TABLE]

For the numerical treatment, we truncate the Fourier series expansions of all appearing functions at an index $N\in\mathbb{N}$ creating Fourier series approximations of the functions. Next, we approximate the Fourier coefficients $\boldsymbol{y}_{k}=(y_{k}^{c},y_{k}^{s})^{T}\in\mathbb{V}$ and $\boldsymbol{p}_{k}=(p_{k}^{c},p_{k}^{s})^{T}\in\mathbb{V}$ of the unknown state and adjoint state functions by finite element functions $\boldsymbol{y}_{kh}=(y_{kh}^{c},y_{kh}^{s})^{T}\in\mathbb{V}_{h}$ and $\boldsymbol{p}_{kh}=(p_{kh}^{c},p_{kh}^{s})^{T}\in\mathbb{V}_{h}$. The finite element space $\mathbb{V}_{h}=V_{h}\times V_{h}\subset\mathbb{V}$ with $V_{h}=\mbox{span}\{\phi_{1},\dots,\phi_{n}\}$ and $\{\phi_{i}(\boldsymbol{x}):i=1,2,\dots,n_{h}\}$ is conforming. We have defined $n=n_{h}=\mbox{dim}V_{h}=O(h^{-d})$ and $h$ is the discretization parameter. The basis of the finite element space $V_{h}$ on the triangulation $\mathcal{T}_{h}$, which is regular, consists of piecewise linear and continuous elements.

5.1. Optimization problem I

The MhFE discretization (see also [27] for more details on the multiharmonic finite element analysis for such an optimal control problem) yields the system of linear equations having a saddle point structure corresponding to the decoupled problems (12)–(13) as follows

[TABLE]

The functions $y_{kh}^{c}(\boldsymbol{x})=\sum_{j=1}^{n}y_{k,j}^{c}\phi_{j}(\boldsymbol{x})$, $y_{kh}^{s}(\boldsymbol{x})=\sum_{j=1}^{n}y_{k,j}^{s}\phi_{j}(\boldsymbol{x})$, $p_{kh}^{c}(\boldsymbol{x})=\sum_{j=1}^{n}p_{k,j}^{c}\phi_{j}(\boldsymbol{x})$ and $p_{kh}^{s}(\boldsymbol{x})=\sum_{j=1}^{n}p_{k,j}^{s}\phi_{j}(\boldsymbol{x})$ are defined by the corresponding nodal function values $\underline{y}_{k}^{c}=(y_{k,j}^{c})_{j=1,\dots,n},\,\underline{y}_{k}^{s}=(y_{k,j}^{s})_{j=1,\dots,n},\,\underline{p}_{k}^{c}=(p_{k,j}^{c})_{j=1,\dots,n},\,\underline{p}_{k}^{s}=(p_{k,j}^{s})_{j=1,\dots,n}\in\mathbb{R}^{n}$. We have defined the stiffness matrix $K_{h,\nu}$ as well as the mass matrices $M_{h}$ and $M_{h,\sigma}$ by their entries

[TABLE]

The right-hand side can be obtained by computing the vectors

[TABLE]

Note that for the case $k=0$, hence for (14)–(15), we obtain

[TABLE]

Solving all the systems of linear equations finally lead to the contributions for computing the MhFE approximations $y_{Nh}(\boldsymbol{x},t)$ and $p_{Nh}(\boldsymbol{x},t)$ given by $y_{Nh}(\boldsymbol{x},t)=y_{0h}^{c}(\boldsymbol{x})+\sum_{k=1}^{N}(y_{kh}^{c}(\boldsymbol{x})\cos(k\omega t)$ $+y_{kh}^{s}(\boldsymbol{x})\sin(k\omega t))$, $p_{Nh}(\boldsymbol{x},t)=p_{0h}^{c}(\boldsymbol{x})+\sum_{k=1}^{N}(p_{kh}^{c}(\boldsymbol{x})\cos(k\omega t)+p_{kh}^{s}(\boldsymbol{x})\sin(k\omega t))$. For some proper fast solvers for these systems we refer to [17, 21, 27]. Both, majorant (16) and minorant (34) of the cost functional $\mathcal{J}$ can be computed by choosing the MhFE approximations $y_{Nh}$, $p_{Nh}$ and $u_{Nh}=-\lambda^{-1}p_{Nh}$ as $\eta$, $\zeta$ and $v$, respectively. The arbitrary functions $\boldsymbol{\tau}$ and $\boldsymbol{\rho}$ can also be represented in form of multiharmonic functions $\boldsymbol{\tau}_{Nh}=\boldsymbol{\tau}_{0h}^{c}(\boldsymbol{x})+\sum_{k=1}^{N}(\boldsymbol{\tau}_{kh}^{c}(\boldsymbol{x})\cos(k\omega t)+\boldsymbol{\tau}_{kh}^{s}(\boldsymbol{x})\sin(k\omega t))$ and $\boldsymbol{\rho}_{Nh}=\boldsymbol{\rho}_{0h}^{c}(\boldsymbol{x})+\sum_{k=1}^{N}(\boldsymbol{\rho}_{kh}^{c}(\boldsymbol{x})\cos(k\omega t)+\boldsymbol{\rho}_{kh}^{s}(\boldsymbol{x})\sin(k\omega t))$. Hence, majorant (16) and minorant (34) have a multiharmonic structure too. The linearity of the problem again yields the decoupling of the problems introducing $\alpha_{k},\beta_{k}>0$ and resulting into majorants $\mathcal{J}_{k}^{\oplus}$ and minorants $\mathcal{J}_{k}^{\ominus}$ corresponding to every Fourier mode. We start with the majorants $\mathcal{J}_{0}^{\oplus}=\mathcal{J}_{0}^{\oplus}(\alpha_{0},\beta_{0};y_{0h}^{c},p_{0h}^{c},\boldsymbol{\tau}_{0h}^{c})$ and $\mathcal{J}_{k}^{\oplus}=\mathcal{J}_{k}^{\oplus}(\alpha_{k},\beta_{k};\boldsymbol{y}_{kh},\boldsymbol{p}_{kh},\boldsymbol{\tau}_{kh})$ together with defining the parameter $\gamma_{k}:=((1+\alpha_{k})(1+\beta_{k})C_{F}^{2})/(2\alpha_{k}\underline{\mu_{1}}^{2})$. We have that

[TABLE]

and

[TABLE]

Defining $\boldsymbol{\alpha}_{N+1}=(\alpha_{0},\dots,\alpha_{N+1})^{T}$ and $\boldsymbol{\beta}_{N}=(\beta_{0},\dots,\beta_{N})^{T}$, we can write the overall majorant as

[TABLE]

Here, the terms are ${\mathcal{R}_{1}}^{c}_{0}=\nabla\cdot\boldsymbol{\tau}_{0h}^{c}-\lambda^{-1}p_{0h}^{c}$, ${\mathcal{R}_{2}}^{c}_{0}=\boldsymbol{\tau}_{0h}^{c}-\nu\nabla y_{0h}^{c}$,

[TABLE]

and ${\mathcal{R}_{2}}_{k}=\boldsymbol{\tau}_{kh}-\nu\nabla\boldsymbol{y}_{kh}=({\mathcal{R}_{2}}^{c}_{k},{\mathcal{R}_{2}}^{s}_{k})^{T}=(\boldsymbol{\tau}_{kh}^{c}-\nu\nabla y_{kh}^{c},\boldsymbol{\tau}_{kh}^{s}-\nu\nabla y_{kh}^{s})^{T}$. The truncation’s remainder term $\mathcal{E}_{N}:=\|y_{d}-{y_{d}}_{N}\|^{2}=\frac{T}{2}\sum_{k=N+1}^{\infty}\|\boldsymbol{y_{d}}_{k}\|_{\Omega}^{2}=\frac{T}{2}\sum_{k=N+1}^{\infty}\left(\|{y_{d}^{c}}_{k}\|_{\Omega}^{2}+\|{y_{d}^{s}}_{k}\|_{\Omega}^{2}\right)$ can always be computed with any desired accuracy, since $y_{d}$ is known (see also [25]). The minorant (34) can be written as

[TABLE]

where $\mathcal{J}_{0}^{\ominus}=\mathcal{J}_{0}^{\ominus}(y_{0h}^{c},p_{0h}^{c},\boldsymbol{\tau}_{0h}^{c},\boldsymbol{\rho}_{0h}^{c})$ and $\mathcal{J}_{k}^{\ominus}=\mathcal{J}_{k}^{\ominus}(\boldsymbol{y}_{kh},\boldsymbol{p}_{kh},\boldsymbol{\tau}_{kh},\boldsymbol{\rho}_{kh})$ are given by

[TABLE]

with ${\mathcal{R}_{3}}^{c}_{0}=\nabla\cdot\boldsymbol{\rho}_{0h}^{c}+y_{0h}^{c}-{y_{d}}_{0}^{c}$, ${\mathcal{R}_{4}}^{c}_{0}=\boldsymbol{\rho}_{0h}^{c}-\nu\nabla p_{0h}^{c}$, and

[TABLE]

with ${\mathcal{R}_{3}}_{k}=k\omega\,\sigma\boldsymbol{p}_{kh}^{\perp}+\nabla\cdot\boldsymbol{\rho}_{kh}+\boldsymbol{y}_{kh}-{\boldsymbol{y}_{d}}_{k}=({\mathcal{R}_{3}}^{c}_{k},{\mathcal{R}_{3}}^{s}_{k})^{T}=(-k\omega\,\sigma p_{kh}^{s}+\nabla\cdot\boldsymbol{\rho}_{kh}^{c}+y_{kh}^{c}-{y_{d}}_{k}^{c},k\omega\,\sigma p_{kh}^{c}+\nabla\cdot\boldsymbol{\rho}_{kh}^{s}+y_{kh}^{s}-{y_{d}}_{k}^{s})^{T}$ and ${\mathcal{R}_{4}}_{k}=\boldsymbol{\rho}_{kh}-\nu\nabla\boldsymbol{p}_{kh}=({\mathcal{R}_{4}}^{c}_{k},{\mathcal{R}_{4}}^{s}_{k})^{T}=(\boldsymbol{\rho}_{kh}^{c}-\nu\nabla p_{kh}^{c},\boldsymbol{\rho}_{kh}^{s}-\nu\nabla p_{kh}^{s})^{T}$.

Remark 5**.**

For any index $\bar{N}\in\mathbb{N}$, $\bar{N}>N$, the truncated remainder term

[TABLE]

is a fully computable lower bound for the remainder term $\mathcal{E}_{N}$. For any given $y_{d}\in L^{2}(Q)$, this provides an arbitrarily tight lower bound for the minorant (72), which is in return optimal.

The fluxes $\boldsymbol{\tau}_{0h}^{c}$, $\boldsymbol{\rho}_{0h}^{c}$ and $\boldsymbol{\tau}_{kh}$, $\boldsymbol{\rho}_{kh}$ for all $k=1,\dots,N$, denoted by $\boldsymbol{\tau}_{kh}=R_{h}^{\text{\tiny{flux}}}(\nu\nabla\boldsymbol{y}_{kh})$ and $\boldsymbol{\rho}_{kh}=R_{h}^{\text{\tiny{flux}}}(\nu\nabla\boldsymbol{p}_{kh})$ are reconstructed by lowest-order Raviart-Thomas elements mapping $L^{2}$-functions to $H(\text{div},\Omega)$, see [38] as well as [25, 26], leading to $\boldsymbol{\tau}_{Nh}=R_{h}^{\text{\tiny{flux}}}(\nu\nabla y_{Nh})$ and $\boldsymbol{\rho}_{Nh}=R_{h}^{\text{\tiny{flux}}}(\nu\nabla p_{Nh})$. We minimize $\mathcal{J}^{\oplus}$ with respect to the positive parameters $\alpha_{k}$ and $\beta_{k}$ leading to the optimized $\boldsymbol{\alpha}_{N+1}$ and $\boldsymbol{\beta}_{N}$. Finally, the multiharmonic majorant (71) and minorant (72) lead to upper and lower bounds for $\mathcal{J}$ which are guaranteed and computable. The computation of the infimum of $\mathcal{J}^{\oplus}$ and the supremum of $\mathcal{J}^{\ominus}$ provide the minimum of $\mathcal{J}$, see also [26].

5.2. Optimization problem II

For the second problem, we only summarize the main results and changes. The MhFE discretization leads to the following discrete problem:

[TABLE]

and

[TABLE]

Now, the right-hand side vectors are computed by

[TABLE]

We summarize the discrete majorant (40) and minorant (46) of the cost functional $\tilde{\mathcal{J}}$ computed by choosing the MhFE approximations for all the (arbitrary) functions. Defining now $\tilde{\mathcal{E}}_{N}:=\|\boldsymbol{g}_{d}-{\boldsymbol{g}_{d}}_{N}\|^{2}=\frac{T}{2}\sum_{k=N+1}^{\infty}\|{\boldsymbol{g}_{d}}_{k}\|_{\Omega}^{2}=\frac{T}{2}\sum_{k=N+1}^{\infty}\left(\|{\boldsymbol{g}_{d}^{c}}_{k}\|_{\Omega}^{2}+\|{\boldsymbol{g}_{d}^{s}}_{k}\|_{\Omega}^{2}\right)$ as truncation’s remainder term, we write the majorant (40) as

[TABLE]

where $\tilde{\mathcal{J}}_{0}^{\oplus}=\tilde{\mathcal{J}}_{0}^{\oplus}(\alpha_{0},\beta_{0};y_{0h}^{c},p_{0h}^{c},\boldsymbol{\tau}_{0h}^{c})$ and $\tilde{\mathcal{J}}_{k}^{\oplus}=\tilde{\mathcal{J}}_{k}^{\oplus}(\alpha_{k},\beta_{k}s;\boldsymbol{y}_{kh},\boldsymbol{p}_{kh},\boldsymbol{\tau}_{kh})$ are given by

[TABLE]

and

[TABLE]

The minorant (46) can be written as

[TABLE]

where $\tilde{\mathcal{J}}_{0}^{\ominus}=\tilde{\mathcal{J}}_{0}^{\ominus}(y_{0h}^{c},p_{0h}^{c},\boldsymbol{\tau}_{0h}^{c},\boldsymbol{\rho}_{0h}^{c})$ and $\tilde{\mathcal{J}}_{k}^{\ominus}=\tilde{\mathcal{J}}_{k}^{\ominus}(\boldsymbol{y}_{kh},\boldsymbol{p}_{kh},\boldsymbol{\tau}_{kh},\boldsymbol{\rho}_{kh})$ are given by

[TABLE]

and

[TABLE]

6. Robust preconditioners for the minimal residual method

The saddle point systems (61), (68), (87) and (94) can be solved by using the preconditioned MinRes method, see [37]. A convergence result for the preconditioned MinRes method can be found in [13] stating that the convergence rate of the preconditioned MinRes method depends on the condition number of the preconditioned system. The derivation of preconditioners for problems (61) for $k=1,\dots,N$ and (68) for $k=0$ have already been presented and discussed in [17, 27] given by

[TABLE]

respectively, where $D_{k}=\sqrt{\lambda}K_{h,\nu}+k\omega\sqrt{\lambda}M_{h,\sigma}+M_{h}$ and $D_{0}=M_{h}+\sqrt{\lambda}K_{h,\nu}$. In [17], preconditioners are derived following the technique in [46] based on operator interpolation theory.

In this section, we present new robust preconditioners for the problem matrices in (87) for $k=1,\dots,N$ and in (94) for $k=0$ in order to solve optimization problem II. Here, we assume that $\sigma$ and $\nu$ are constant, which we also choose in the numerical results in Section 7. Hence, $M_{h,\sigma}=\sigma M_{h}$ and $K_{h,\nu}=\nu K_{h}$. The block-diagonal preconditioners are practically implemented by the version of the algebraic multilevel iteration (AMLI) method from [20]. The AMLI preconditioned MinRes solver is robust and of optimal complexity which is proved in [21]. This can be also observed in the numerical results’ section. Let us consider the saddle point system (87) for $k\geq N$. The derivation of preconditioners for system (94) in case of $k=0$ is completely analogous. We refer the reader to [17] on details for the derivation of the preconditioners based on Schur complements.

Defining the matrices and vectors

[TABLE]

$\underline{f}=(K_{h}{\underline{g}_{d}^{c}}_{k},K_{h}{\underline{g}_{d}^{s}}_{k})^{T}$, $\underline{y}=(\underline{y}_{k}^{c},\underline{y}_{k}^{s})^{T}$ and $\underline{p}=(\underline{p}_{k}^{c},\underline{p}_{k}^{s})^{T}$ leads to the following problem structure

[TABLE]

with the symmetric and positive definite matrices $A$ and $C$. We define the negative Schur complements $S=C+BA^{-1}B^{T}$ and $R=A+B^{T}C^{-1}B$ yielding the preconditioners for $k\geq N$ as follows

[TABLE]

The negative Schur complements are given by

[TABLE]

and

[TABLE]

Let us define

[TABLE]

and

[TABLE]

Then $S$ and $R$ can be written as $S=\text{diag}(\tilde{D}_{k}^{S},\tilde{D}_{k}^{S})$ and $R=\text{diag}(\tilde{D}_{k}^{R},\tilde{D}_{k}^{R})$, respectively.

Remark 6**.**

Both Schur complement preconditioners $\tilde{\mathcal{P}}$ and $\tilde{\mathcal{Q}}$ in (112) can be chosen for computations of optimization problem II leading to fast and robust convergence rates, see [22] and [35].

Inserting now $A$ and $C$ from (104) and $S$ and $R$ from (113) and (114), respectively, the Schur complement preconditioners in the form of $\tilde{\mathcal{P}}$ as presented in (112) in the case of $k\geq N$ as well as for $k=0$ are given by

[TABLE]

Analogously in the form of $\tilde{\mathcal{Q}}$, they are given by

[TABLE]

for $k\geq N$ and $k=0$ for the saddle point systems (87) and (94), respectively. We refer again to [17] for further details on the preconditioners’ derivation. For the numerical experiments of this work, we simplify the preconditioner such that $\tilde{\tilde{\mathcal{P}}}_{k}=\text{diag}(\tilde{\tilde{D}}_{k}^{S},\tilde{\tilde{D}}_{k}^{S},\lambda^{-1}\tilde{\tilde{D}}_{k}^{S},\lambda^{-1}\tilde{\tilde{D}}_{k}^{S})$ and $\tilde{\tilde{\mathcal{P}}}_{0}=\text{diag}(\tilde{\tilde{D}}_{0}^{S},\lambda^{-1}\tilde{\tilde{D}}_{0}^{S})$, where $\tilde{\tilde{D}}_{k}^{S}=\nu\sqrt{\lambda}K_{h}+(1+k\omega\sigma\sqrt{\lambda})M_{h}$ and $\tilde{\tilde{D}}_{0}^{S}=M_{h}+\nu\sqrt{\lambda}K_{h}$ in a similar form as (101) in order to apply the AMLI preconditioned MinRes method analogously as for optimization problem I.

7. Numerical results

In this section, we present numerical results performed in C++ for computing the minorants and majorants of the two optimal control problems with cost functionals (9) and (21) and for different cases of given data denoted by Examples 1–6. For the numerical results on the majorants of optimization problem I, we refer to [26] and for the numerical analysis including convergence results, to [17, 27]. However, all the numerical experiments on the majorants of optimization problem II as well as all on minorants for both problems I and II are new. We perform numerical experiments for the same three cases (applied on the problem data) as for problem I but in Examples 4–6 they are applied on the desired gradients.

The domain $\Omega$ is the two-dimensional unit square $\Omega=(0,1)^{2}$ with Friedrichs’ constant $C_{F}=1/(\sqrt{2}\pi)$. The triangulation of the domain is performed regular leading to a uniform grid. The finite elements are chosen as described in Section 5. The coefficients are chosen to be $\nu=\sigma=1$. In Examples 1, 2, 4 and 5, the cost parameter is chosen to be $\lambda=0.1$, $T=2\pi/\omega$ and $\omega=1$. In Examples 3 and 6, $\lambda=0.01$, $T=1$ and $\omega=2\pi$. The MhFE approximations for $\eta$, $\zeta$ and $\boldsymbol{\tau},\boldsymbol{\rho}$ are chosen as well as the fluxes are reconstructed by $RT^{0}$-extensions (lowest-order standard Raviart-Thomas) resulting in averaged fluxes which are now from $H(\text{div},\Omega)$, see also [25, 26] for further details. The grid sizes range between 16$\times$16 and 256$\times$256 as well as 512$\times$512 to obtain finer grid solutions for Examples 3 and 6 as a reference for the exact solution. The preconditioned MinRes iteration was stopped after 8 iteration steps in all computations using the AMLI preconditioner with 4 inner iterations. The numerical experiments for Examples 1, 2, 4 and 5 were performed on a laptop with Intel(R) Core(TM) i5-6267U CPU @ 2.90GHz processor and 16 GB 2133 MHz LPDDR3 memory. The numerical experiments for Examples 3 and 6 were performed on a CPU server with a Tumbleweed distribution having 64 cores and 1 terabyte memory in order to provide enough memory for computing the finer grid solutions in addition. All computational times $t^{\text{sec}}$ are presented in seconds and include the CPU times also needed to derive the majorants and minorants. However, we want to highlight that these times are much smaller compared to the rest. The computational times of Examples 3 and 6 exclude the computation of the solution on the finer grid (512$\times$512).

7.1. Numerical results for optimization problem I

7.1.1. Example 1

The desired state is a time-periodic and time-analytic function

[TABLE]

which is however not time-harmonic. Note that the exact state function for this example is

[TABLE]

The truncation index for the multiharmonic approximations is chosen as $N=8$ here. Table 1 presents for different grid sizes CPU times $t^{\text{sec}}$, values for the majorants $\mathcal{J}^{\oplus}_{0}$ and minorants $\mathcal{J}^{\ominus}_{0}$ as defined in (69) and (73), and efficiency indices $I_{\text{eff}}^{\mathcal{J}^{\oplus}_{0}}=\mathcal{J}^{\oplus}_{0}/\mathcal{J}_{0}$, $I_{\text{eff}}^{\mathcal{J}^{\ominus}_{0}}=\mathcal{J}^{\ominus}_{0}/\mathcal{J}_{0}$ and $I_{\text{eff}}^{\mathcal{J},0}=\mathcal{J}^{\oplus}_{0}/\mathcal{J}^{\ominus}_{0}$. Here, $\mathcal{J}_{0}=\mathcal{J}_{0}(y_{0}^{c},u_{0}^{c})=\frac{1}{2}\|y_{0}^{c}-{y_{d}^{c}}_{0}\|_{\Omega}^{2}+\frac{\lambda}{2}\|u_{0}^{c}\|_{\Omega}^{2}$ as introduced in Subsection 2.2. In Table 2, the numerical results for the Fourier mode $k=1$ are presented including $\mathcal{J}^{\oplus}_{k}$, $\mathcal{J}^{\ominus}_{k}$ as defined in (70) and (74) and the corresponding efficiency indices $I_{\text{eff}}^{\mathcal{J}^{\oplus}_{k}}=\mathcal{J}^{\oplus}_{k}/\mathcal{J}_{k}$, $I_{\text{eff}}^{\mathcal{J}^{\ominus}_{k}}=\mathcal{J}^{\ominus}_{k}/\mathcal{J}_{k}$ and $I_{\text{eff}}^{\mathcal{J},k}=\mathcal{J}^{\oplus}_{k}/\mathcal{J}^{\ominus}_{k}$. Moreover, we present the efficiency indices for $\mathcal{M}_{1}^{\oplus}$ given for the modes by

[TABLE]

The error norms for the modes are given by

[TABLE]

and

[TABLE]

leading the representation

[TABLE]

with the remainder term $\mathcal{F}_{N}:=\frac{T}{2}\sum_{k={N+1}}^{\infty}|||\boldsymbol{y}_{k}|||_{1,k}^{2}$. For the numerical experiments, we can estimate the efficiency index for $\mathcal{M}_{1}^{\oplus}$ from above by estimating (121) from below ignoring the remainder term $\mathcal{F}_{N}$ leading to the overall efficiency index for $\mathcal{M}_{1}^{\oplus}$

[TABLE]

The corresponding majorants $\mathcal{M}_{1,0}^{\oplus}=\mathcal{M}_{1,0}^{\oplus}(\alpha_{0},\beta_{0};y_{0h}^{c},p_{0h}^{c},\boldsymbol{\tau}_{0h}^{c},\boldsymbol{\rho}_{0h}^{c})$ and $\mathcal{M}_{1,k}^{\oplus}=\mathcal{M}_{1,k}^{\oplus}(\alpha_{k},\beta_{k};$ $\boldsymbol{y}_{kh},\boldsymbol{p}_{kh},\boldsymbol{\tau}_{kh},\boldsymbol{\rho}_{kh})$ are given by

[TABLE]

and also $\mathcal{M}_{1,k}^{\oplus}=\mathcal{J}_{k}^{\oplus}-\mathcal{J}_{k}^{\ominus}+\frac{3\lambda}{4C_{F}^{2}}\left(C_{F}\|{\mathcal{R}_{1}}_{k}\|_{\Omega}+\|{\mathcal{R}_{2}}_{k}\|_{\Omega}\right)^{2}$. Table 3 sums up the numerical results for Example 1 by presenting the minorants, majorants and efficiency indices on a grid of size $256\times 256$ for all $k$ up to $N=4$ and then for $k=6$ and $k=8$ (since their results were similar as for $k=5$ and $k=7$).

For $N=3$ and $N=8$, the truncation’s remainder terms can be precomputed and are given by $\mathcal{E}_{3}=63694.86$ and $\mathcal{E}_{8}=106.06$, respectively. Since the overall efficiency indices in Table 3 stay all in approximately the same range, we observe that the method is robust. However, the efficiency indices for the combined error norm $I_{\text{eff}}^{\mathcal{M}_{1}}$ indicate that the modes $k=1$ and $k=4$ are the most significant to represent the solution by its multiharmonic approximation. Comparing the last two lines of Table 3 shows that the value for representing the cost functional of the exact solution is already sufficiently accurate for a truncation index $N=3$. One of the reasons for this is that the remainder term $\mathcal{E}_{N}$ can be precomputed exactly.

7.1.2. Example 2

We choose the time-analytic, not time-periodic, desired state function

[TABLE]

having as exact solution the state function

[TABLE]

The approximations by the MhFEM are computed for the truncation index $N=10$. For Example 2, it suffices to present here only the overall results as in Table 3, now presented in Table 4. We compare the overall values of majorants and minorants for different truncation indices $N=6$ and $N=10$, for which the corresponding truncation’s remainder terms are given by $\mathcal{E}_{6}=44094.84$ and $\mathcal{E}_{10}=10597.20$. One can see from the last two lines that the truncation index $N=6$ suffices already to provide an accurate enough approximate solution. Also efficiency indices being around 1 show that the majorants and minorants perform well for that example.

7.1.3. Example 3

The desired state is chosen to be a space-time non-smooth function

[TABLE]

Here, we denote by $\chi$ the space-time characteristic function. The desired state’s Fourier cofficients are analytically computed and given by ${y_{d}^{c}}_{0}(\boldsymbol{x})=\chi_{[\frac{1}{2},1]^{2}}(\boldsymbol{x})/2$ and

[TABLE]

The desired state has jumps in space and time. In this example, the exact solution cannot be precomputed analytically. Hence, as approximation for it, we use a MhFE representation on the finer grid with size $512\times 512$. Note also that the Fourier coefficients are zero for the even Fourier modes besides for $k=0$. We can observe in Table 5 that the efficiency indices, especially regarding the combined norm, are similar for higher modes. We have computed the results up to $N=11$ and also added the results for the modes $k=21$, $k=41$ and $k=81$ to the table as examples. The majorants of the combined norm stay approximately in the same range for $k\geq 1$. The majorants and minorants for the cost functional are close to 1, which demonstrates their efficiency also in this numerical example, where the given data has jumps in space and time.

7.2. Numerical results for optimization problem II

We compute the numerical results for the three same cases as for problem I but now applied on the desired gradient.

7.2.1. Example 4

We set the desired gradient to be time-periodic and time-analytic

[TABLE]

The exact solution for the state function is given by (120). Moreover, we present the efficiency indices for $\tilde{\mathcal{M}}_{1}^{\oplus}$ given for the modes by

[TABLE]

The error norms for the modes are given by

[TABLE]

and

[TABLE]

leading the representation

[TABLE]

with the remainder term $\tilde{\mathcal{F}}_{N}:=\frac{T}{2}\sum_{k={N+1}}^{\infty}||||\boldsymbol{y}_{k}||||_{1,k}^{2}$. For the numerical experiments, the efficiency index for $\tilde{\mathcal{M}}_{1}^{\oplus}$ from above by estimating (127) from below ignoring the remainder term $\tilde{\mathcal{F}}_{N}$ leading to the overall efficiency index for $\tilde{\mathcal{M}}_{1}^{\oplus}$ given by

[TABLE]

The corresponding majorants are given by $\tilde{\mathcal{M}}_{1,0}^{\oplus}=\tilde{\mathcal{M}}_{1,0}^{\oplus}(\alpha_{0},$ $\beta_{0};y_{0h}^{c},p_{0h}^{c},\boldsymbol{\tau}_{0h}^{c},\boldsymbol{\rho}_{0h}^{c})=\tilde{\mathcal{J}}_{0}^{\oplus}-\tilde{\mathcal{J}}_{0}^{\ominus}+\frac{3\lambda}{4C_{F}^{2}}(C_{F}\|{\mathcal{R}_{1}}^{c}_{0}\|_{\Omega}+\|{\mathcal{R}_{2}}^{c}_{0}\|_{\Omega})^{2}$ and $\tilde{\mathcal{M}}_{1,k}^{\oplus}=\tilde{\mathcal{M}}_{1,k}^{\oplus}(\alpha_{k},$ $\beta_{k};\boldsymbol{y}_{kh},\boldsymbol{p}_{kh},\boldsymbol{\tau}_{kh},\boldsymbol{\rho}_{kh})$ $=\tilde{\mathcal{J}}_{k}^{\oplus}-\tilde{\mathcal{J}}_{k}^{\ominus}+\frac{3\lambda}{4C_{F}^{2}}\left(C_{F}\|{\mathcal{R}_{1}}_{k}\|_{\Omega}+\|{\mathcal{R}_{2}}_{k}\|_{\Omega}\right)^{2}$. We present the numerical results for the modes $k=0$ and $k=1$ for different grid sizes in Tables 6 and 7. The efficiency indices for the majorants and minorants are very close to 1.00. Also the efficiency indices for $\tilde{\mathcal{M}}_{1,0}$ show a good accuracy.

Table 8 compares the results for different Fourier modes up to $N=8$ computed on a grid of size $256\times 256$. Here, the overall minorants, majorants and efficiency indices are presented, where the remainder terms $\mathcal{E}_{N}$ for $N=8$ and also $N=6$ have been precomputed exactly. The values of the efficiency indices vary for different modes $k$. For example, the results for $\tilde{\mathcal{M}}_{1,4}$ indicate that the mode $k=4$ is essential to represent the solution accurately. The values for $I_{\text{eff}}^{\mathcal{J}^{\ominus}_{7}}$ and $I_{\text{eff}}^{\mathcal{J}^{\ominus}_{8}}$ indicate that the minorants require a different, higher refinement for a more accurate representation of the overall solution. An adaptive scheme is the natural choice. On the other hand, in these cases the majorants give a good representation for the cost functional. Finally, comparing the last two lines of Table 8 again shows that the overall value for representing the cost functional of the exact solution is already sufficiently accurate for a truncation index $N=6$.

7.2.2. Example 5

We choose the non time-periodic but time-analytic desired gradient

[TABLE]

leading to the time-analytic, but not time-periodic exact state (124). We compute the MhFE approximation of the desired gradient and solve the systems (61) and (68) for modes up to $N=10$ on a $256\times 256$-mesh and present the results in Table 9. The remainder terms for $N=6$ and $N=10$ are $\mathcal{E}_{6}=4796.54$ and $\mathcal{E}_{10}=1149.65$, respectively. The efficiency indices for the overall majorant and minorant show that a truncation index of $N=6$ already gives a sufficiently accurate approximation for the overall cost functional. Note that the efficiency index for $\tilde{\mathcal{M}}_{1,2}$ indicates that the mode $k=2$ is essential for the multiharmonic approximation giving an accurate representation of the solution.

7.2.3. Example 6

We set the space-time non-smooth desired gradient

[TABLE]

Also the coefficients of the Fourier expansion associated with $\boldsymbol{g}_{d}$ can be found analytically. They are as in Example 3 given by (126) for each direction of the gradient (128). Again the exact solution cannot be computed analytically and hence we use its MhFE approximations on the finer mesh of size $512\times 512$ as a reference. Table 10 presents the results for modes up to truncation index $N=11$ as well as for $k=21$, $k=41$ and $k=81$ analogously to Example 3. The results reflected by the efficiency indices show the good representation by using the minorants and majorants, especially, considering the efficiency indices in the last two columns of Table 10. This again demonstrates the efficiency of the minorants and majorants for data having jumps in space and time but now for optimization problem II.

8. Conclusions and outlook

In this work, the a posteriori error analysis started in [26] has been extended now by deriving new lower bounds, called minorants, for the cost functional leading to an upper estimate for the error norm of the state and control or equivalently in state and adjoint state. These lower bounds are guaranteed and computable. Together with using the results from [44] as well as [16] one can apply the method also to time-periodic optimal control problems, where box constraints are being imposed on the Fourier coefficients of the control. The estimates are derived for two different cost functionals, where the second one is now new in this context.

Since in the linear case the problems are decoupled, the solutions on the Fourier coefficients could easily be computed on grids of different sizes depending on the accuracy needed, which could be exactly determined by using the a posteriori estimates presented in this work leading to an adaptive method in time. Together with the adaptive finite element method we then obtain a space-time adaptive method, the adaptive multiharmonic finite element method, which we call AMhFEM, as mentioned for the first time in [26].

In this work, a first derivation of preconditioners for the MinRes method for the second optimization problem has been presented as well as a preconditoner for applying AMLI has been suggested. Several numerical tests for optimization problem I and II have been presented showing the efficiency of the upper and especially – with regard to the article – lower bounds for the cost functionals in practice.

Acknowledgments

The author gratefully acknowledges the financial support by the Academy of Finland under the grant 295897. The author would like to thank the anonymous referees for their valuable comments improving the article.