Multigoal-oriented optimal control problems with nonlinear PDE constraints

TL;DR

This paper develops an adaptive solution strategy for multigoal-oriented optimal control problems constrained by nonlinear PDEs, using a dual-weighted residual method to balance errors and improve accuracy.

Contribution

It introduces a combined a posteriori error estimator for multiple quantities of interest in nonlinear PDE-constrained control problems, enabling adaptive mesh refinement.

Findings

Effective error balancing between discretization and nonlinear iteration.

Enhanced accuracy in control solutions through adaptive mesh refinement.

Numerical examples demonstrate the method's efficiency and robustness.

Abstract

In this work, we consider an optimal control problem subject to a nonlinear PDE constraint and apply it to the regularized $p$-Laplace equation. To this end, a reduced unconstrained optimization problem in terms of the control variable is formulated. Based on the reduced approach, we then derive an a posteriori error representation and mesh adaptivity for multiple quantities of interest. All quantities are combined to one, and then the dual-weighted residual (DWR) method is applied to this combined functional. Furthermore, the estimator allows for balancing the discretization error and the nonlinear iteration error. These developments allow us to formulate an adaptive solution strategy, which is finally substantiated via several numerical examples.

Figures

Tables

Table 1. Table 1: Example 1: I effs subscript 𝐼 effs I_{\text{effs}} for I ​ ( u , q ) := ∫ Ω u ​ ( x ) 2 ​ q ​ ( x ) 2 ​ 𝑑 x assign 𝐼 𝑢 𝑞 subscript Ω 𝑢 superscript 𝑥 2 𝑞 superscript 𝑥 2 differential-d 𝑥 I(u,q):=\int_{\Omega}u(x)^{2}q(x)^{2}dx .

$\alpha$	0.01		0.1		1		10
$l$	$I_{\text{eff}}$	DOFs	$I_{\text{eff}}$	DOFs	$I_{\text{eff}}$	DOFs	$I_{\text{eff}}$	DOFs
0	0.88	275	0.69	275	0.89	275	0.90	275
1	0.94	326	0.79	506	0.91	565	0.93	568
2	0.94	381	0.68	759	0.91	832	0.92	845
3	0.99	561	0.66	1 266	0.92	1 367	0.93	1 451
4	1.05	719	0.63	2 084	0.91	2 246	0.93	2 385
5	1.06	1 151	0.50	3 013	0.89	3 115	0.93	3 263
6	1.13	1 856	0.59	5 031	0.92	5 072	0.95	5 444
7	1.05	2 419	0.55	8 137	0.94	8 367	0.97	8 865
8	1.11	3 363	0.36	12 498	0.94	11 880	0.97	12 479
9	1.12	5 691	0.56	20 690	0.95	17 591	0.98	19 357
10	1.15	7 852	0.47	33 247	0.95	31 035	0.99	32 970
11	1.13	10 752	0.38	50 864	0.95	45 721	0.99	47 850
12	1.14	17 094	0.56	84 368	0.96	72 636	0.99	78 502
13	1.19	25 916	0.44	135 166	0.96	126 711	0.99	133 541
14	1.14	35 482	0.39	207 466	0.96	184 754	1.00	192 946
	$I(u,q)$	DOFs	$I(u,q)$	DOFs	$I(u,q)$	DOFs	$I(u,q)$	DOFs
$\infty$	0.2316036	1 326 503	0.07069658	2 127 499	0.1502366	1 996 755	0.1635741	2 107 007

Table 2. Table 2: Example 3: Comparison of Newton’s Method with adaptive stopping rule (AN) and the classical, non-adaptive, Newton method (FN); I ​ t F 𝐼 subscript 𝑡 𝐹 It_{F} : number of iterations for FN, I ​ t A 𝐼 subscript 𝑡 𝐴 It_{A} : number of iterations for AN, | 𝒯 h , F | subscript 𝒯 ℎ 𝐹 |\mathcal{T}_{h,F}| : number of elements in the adaptive mesh resulting from using FN, | 𝒯 h , A | subscript 𝒯 ℎ 𝐴 |\mathcal{T}_{h,A}| : number of elements in the adaptive mesh resulting from using AN, I eff,F subscript 𝐼 eff,F I_{\text{eff,F}} : I eff subscript 𝐼 eff I_{\text{eff}} for FN, I eff,A subscript 𝐼 eff,A I_{\text{eff,A}} : I eff subscript 𝐼 eff I_{\text{eff}} for AN, I eff,c,F subscript 𝐼 eff,c,F I_{\text{eff,c,F}} : I eff,c subscript 𝐼 eff,c I_{\text{eff,c}} for FN, I eff,c,A subscript 𝐼 eff,c,A I_{\text{eff,c,A}} : I eff,c subscript 𝐼 eff,c I_{\text{eff,c}} for AN

$l$	$I{t}_F$	$I{t}_A$	$\left|{𝒯}_{h,F}\right|$	$\left|{𝒯}_{h,A}\right|$	$I_{\text{eff,F}}$	$I_{\text{eff,A}}$	$I_{\text{eff,c,F}}$	$I_{\text{eff,c,A}}$
$0$	$4$	$3$	$116$	$116$	$0.882$	$0.882$	$0.882$	$0.882$
$1$	$5$	$3$	$137$	$137$	$0.829$	$0.830$	$0.829$	$0.830$
$2$	$5$	$1$	$158$	$158$	$0.942$	$0.933$	$0.942$	$0.941$
$3$	$6$	$2$	$215$	$215$	$0.918$	$0.912$	$0.918$	$0.918$
$4$	$6$	$2$	$347$	$347$	$1.023$	$1.019$	$1.023$	$1.022$
$5$	$6$	$2$	$494$	$494$	$1.078$	$1.068$	$1.078$	$1.076$
$6$	$7$	$2$	$800$	$800$	$1.069$	$1.067$	$1.069$	$1.069$
$7$	$12$	$2$	$\mathrm{1\hspace{0.17em}283}$	$\mathrm{1\hspace{0.17em}283}$	$1.091$	$1.087$	$1.091$	$1.090$
$8$	$17$	$2$	$\mathrm{1\hspace{0.17em}898}$	$\mathrm{1\hspace{0.17em}895}$	$1.089$	$1.093$	$1.089$	$1.090$
$9$	$6$	$2$	$\mathrm{2\hspace{0.17em}966}$	$\mathrm{2\hspace{0.17em}957}$	$1.089$	$1.090$	$1.090$	$1.090$
$10$	$9$	$2$	$\mathrm{4\hspace{0.17em}802}$	$\mathrm{4\hspace{0.17em}790}$	$1.085$	$1.088$	$1.085$	$1.083$
$11$	$5$	$2$	$\mathrm{7\hspace{0.17em}097}$	$\mathrm{7\hspace{0.17em}091}$	$1.089$	$1.077$	$1.089$	$1.093$
$12$	$4$	$2$	$\mathrm{11\hspace{0.17em}153}$	$\mathrm{11\hspace{0.17em}099}$	$1.096$	$1.088$	$1.097$	$1.095$
$13$	$2$	$2$	$\mathrm{18\hspace{0.17em}206}$	$\mathrm{18\hspace{0.17em}140}$	$1.128$	$1.129$	$1.122$	$1.123$
$14$	$2$	$2$	$\mathrm{27\hspace{0.17em}341}$	$\mathrm{27\hspace{0.17em}269}$	$1.151$	$1.152$	$1.136$	$1.137$

Table 3. Table 3: Example 3: Comparison of Newton’s Method with adaptive stopping rule (AN) and the classical, non-adaptive, Newton method (FN); η k , F subscript 𝜂 𝑘 𝐹 \eta_{k,F} : iteration error estimate for FN, η k , F subscript 𝜂 𝑘 𝐹 \eta_{k,F} : iteration error estimate for AN.

$l$	Error in $I_{𝔈,F}$	Error in $I_{𝔈,A}$	${\eta }_{k,F}$	${\eta }_{k,A}$
$0$	$5.98\cdot {10}^0$	$5.98\cdot {10}^0$	$-1.99\cdot {10}^{-5}$	$-3.2\cdot {10}^{-6}$
$1$	$2.93\cdot {10}^2$	$2.75\cdot {10}^2$	$1.22\cdot {10}^{-2}$	$-4.19\cdot {10}^{-2}$
$2$	$2.92\cdot {10}^0$	$2.98\cdot {10}^0$	$-5.18\cdot {10}^{-5}$	$2.46\cdot {10}^{-2}$
$3$	$1.23\cdot {10}^0$	$1.24\cdot {10}^0$	$6.72\cdot {10}^{-5}$	$7.71\cdot {10}^{-3}$
$4$	$4.70\cdot {10}^{-1}$	$4.73\cdot {10}^{-1}$	$9.75\cdot {10}^{-6}$	$1.67\cdot {10}^{-3}$
$5$	$2.93\cdot {10}^{-1}$	$2.96\cdot {10}^{-1}$	$4.44\cdot {10}^{-6}$	$2.59\cdot {10}^{-3}$
$6$	$2.40\cdot {10}^{-1}$	$2.40\cdot {10}^{-1}$	$1.99\cdot {10}^{-6}$	$3.48\cdot {10}^{-4}$
$7$	$1.35\cdot {10}^{-1}$	$1.36\cdot {10}^{-1}$	$-2.41\cdot {10}^{-6}$	$5.21\cdot {10}^{-4}$
$8$	$8.13\cdot {10}^{-2}$	$8.09\cdot {10}^{-2}$	$-3.38\cdot {10}^{-6}$	$-2.73\cdot {10}^{-4}$
$9$	$6.00\cdot {10}^{-2}$	$6.00\cdot {10}^{-2}$	$1.94\cdot {10}^{-5}$	$-2.64\cdot {10}^{-5}$
$10$	$3.77\cdot {10}^{-2}$	$3.77\cdot {10}^{-2}$	$3.65\cdot {10}^{-6}$	$-2.03\cdot {10}^{-4}$
$11$	$2.07\cdot {10}^{-2}$	$2.11\cdot {10}^{-2}$	$-8.56\cdot {10}^{-7}$	$3.45\cdot {10}^{-4}$
$12$	$1.57\cdot {10}^{-2}$	$1.58\cdot {10}^{-2}$	$2.61\cdot {10}^{-5}$	$1.17\cdot {10}^{-4}$
$13$	$9.41\cdot {10}^{-3}$	$9.45\cdot {10}^{-3}$	$-5.67\cdot {10}^{-5}$	$-5.96\cdot {10}^{-5}$
$14$	$5.09\cdot {10}^{-3}$	$5.11\cdot {10}^{-3}$	$-7.89\cdot {10}^{-5}$	$-7.85\cdot {10}^{-5}$

Equations

$$\begin{array}{rl}\underset{(u,q)}{\min }J(u,q)& u\in U,q\in Q,\\ \text{s.t. }A(u,q)=0& \text{ in }{V}^{\ast },\end{array}$$

$$A(S(q),q)=0,\forall q\in Q.$$

$$\underset{q}{\min }j(q),q\in Q,$$

$$j^{\prime }(\overline{q})(\delta q)=0\forall \delta q\in Q.$$

$$\mathcal{L}(u,q,z):=J(u,q)-A(u,q)(z),\forall u\in U,q\in Q,z\in V.$$

$$J_u^{\prime }(\bar{u},\bar{q})(\delta u)-A_u^{\prime }(\bar{u},\bar{q})(z)(\delta u)={\mathcal{L}}_u^{\prime }(\bar{u},\bar{q},\bar{z})(\delta u)$$

$$J_q^{\prime }(\bar{u},\bar{q})(\delta q)-A_q^{\prime }(\bar{u},\bar{q})(\bar{z})(\delta q)={\mathcal{L}}_q^{\prime }(\bar{u},\bar{q},\bar{z})(\delta q)$$

$$-A(\bar{u},\bar{q})(\delta z)={\mathcal{L}}_z^{\prime }(\bar{u},\bar{q},\bar{z})(\delta z)$$

$${\mathcal{A}}_p:W_0^{1,p}(\Omega )\times (W_0^{1,p}(\Omega ){)}^{\ast }\mapsto (W_0^{1,p}(\Omega ){)}^{\ast },$$

$${\mathcal{A}}_p(u,q)(v):=$$

$$\underset{(u,q)}{\min }J(q,u)u\in U,q\in Q$$

$$\text{s.t. }{\mathcal{A}}_p(u,q)(v)=0,$$

$$J(q,u)=\frac{1}{2}\Vert u-{\bar{u}}^d{\Vert }_{L^2(\Omega )}^2+\frac{\alpha }{2}\Vert q-{\bar{q}}^d{\Vert }_{L^2(\Omega )}^2,$$

$$Q_{DG}^r:=\{v_h\in L^{\infty }(\Omega ):v_{h|K}\in Q_r(K),\forall K\in {\mathcal{T}}_h\},$$

$$\begin{array}{rl}\underset{(u_h,q_h)}{\min }J(u_h,q_h)& u_h\in U_h,q_h\in Q_h,\\ \text{s.t }A(u_h,q_h)=0& \text{ in }{V}_h^{\ast }.\end{array}$$

$$A(S_h(q_h),q_h)=0\forall q_h\in Q_h.$$

$$\underset{q_h}{\min }{j}_h(q_h)q_h\in Q_h.$$

$$j_h^{\prime }({\bar{q}}_h)(\delta q_h)=0\forall \delta q_h\in Q_h.$$

$$\mathcal{L}(u_h,q_h,z_h):=J(u_h,q_h)-A(u_h,q_h)(z_h),\forall u_h\in U_h,q_h\in Q_h,z_h\in V_h.$$

$$\begin{array}{rl}{J}_u^{\prime }({\bar{u}}_h,{\bar{q}}_h)(\delta u_h)-A_u^{\prime }({\bar{u}}_h,{\bar{q}}_h)({\bar{z}}_h)(\delta u_h)={\mathcal{L}}_u^{\prime }({\bar{u}}_h,{\bar{q}}_h,{\bar{z}}_h)(\delta u_h)& =0\forall \delta u_h\in U_h,\\ J_q^{\prime }({\bar{u}}_h,{\bar{q}}_h)(\delta q_h)-A_q^{\prime }({\bar{u}}_h,{\bar{q}}_h)({\bar{z}}_h)(\delta q_h)={\mathcal{L}}_q^{\prime }({\bar{u}}_h,{\bar{q}}_h,{\bar{z}}_h)(\delta q_h)& =0\forall \delta q_h\in Q_h,\\ -A({\bar{u}}_h,{\bar{q}}_h)(\delta z_h)={\mathcal{L}}_z^{\prime }({\bar{u}}_h,{\bar{q}}_h,{\bar{z}}_h)(\delta z_h)& =0\forall \delta z_h\in V_h,\end{array}$$

$$I(S(\bar{q}),\bar{q})-I(S_h({\tilde{q}}_h),{\tilde{q}}_h)=i(\bar{q})-i({\tilde{q}}_h)+i({\tilde{q}}_h)-i_h({\tilde{q}}_h).$$

$$j^{\prime \prime }(\overline{q})(\delta q,\overline{p})=-i^{\prime }(\overline{q})(\delta q)\forall \delta q\in Q.$$

$$\begin{array}{rl}i(\overline{q})-i({\tilde{q}}_h)& =\frac{1}{2}\rho ({\tilde{q}}_h)(\overline{p}-{\tilde{p}}_h)+\frac{1}{2}{\rho }^{\ast }({\tilde{q}}_h,{\tilde{p}}_h)(\overline{q}-{\tilde{q}}_h)+\rho ({\tilde{q}}_h)({\tilde{p}}_h)+{\mathcal{R}}^{(3)},\end{array}$$

$$\rho ({\tilde{q}}_h)(\cdot )$$

$${\rho }^{\ast }({\tilde{q}}_h,{\tilde{p}}_h)(\cdot )$$

$${\mathcal{R}}^{(3)}:=\frac{1}{2}{\int }_0^1[i^{\prime \prime \prime }({\tilde{q}}_h+se)(e,e,e)+j^{\prime \prime \prime \prime }({\tilde{q}}_h+se)(e,e,e,{\tilde{p}}_h+s{e}^{\ast })+3j^{\prime \prime \prime }({\tilde{q}}_h+se)(e,e,e^{\ast })]s(s-1)ds,$$

$$m(\overline{x})-m({\tilde{x}}_h)=$$

$$=$$

$$m^{\prime }(\overline{x})(e_x)=i^{\prime }(\overline{q})+j^{\prime \prime }(\overline{q})(e,\overline{p})+j^{\prime }(\overline{q})(e^{\ast })=0$$

$$m(\overline{x})-m({\tilde{x}}_h)=\frac{1}{2}{m}^{\prime }({\tilde{x}}_h)(e_x)+{\mathcal{R}}^{(3)}.$$

$$i(\overline{q})-i({\tilde{q}}_h)=$$

$$=$$

$$\rho ({\tilde{q}}_h)(\overline{p}-{\tilde{p}}_h)=-j^{\prime }({\tilde{q}}_h)(\overline{p}-{\tilde{p}}_h)=-J_q^{\prime }(S({\tilde{q}}_h),{\tilde{q}}_h)(\overline{p}-{\tilde{p}}_h)-J_u^{\prime }(S({\tilde{q}}_h),{\tilde{q}}_h)(S^{\prime }({\tilde{q}}_h)(\overline{p}-{\tilde{p}}_h)).$$

Full text

1 Introduction

Optimal control problems with nonlinear PDE constraints have been studied for a long time in many works. In particular, employing the (regularized) $p$-Laplacian (see e.g., [29, 22, 34, 46]) as a nonlinear constraint of an optimal control problem was considered for instance in [17].

In many applications, however, not the entire solution is of interest, but only parts or certain quantities of interest, so-called goal functionals. In the past, often a single goal functional was analyzed. However, it may be of interest to control multiple goal functionals simultaneously [33, 32, 48, 28, 35, 42]. In this paper, these three topics are combined: optimal control, the regularized $p$-Laplacian as a numerical example of a quasi-linear PDE constraint, and multiple goal-oriented a posteriori error estimation.

In the following, we briefly refer to studies that treat parts of the three topics. Optimal control problems (specifically, a priori estimates and optimality conditions) with quasi-linear (as the $p$-Laplacian can be classified) elliptic PDE constraints were considered in [16, 18, 15]. More recently, the extension to optimal control with parabolic PDEs was discussed in [9] and [14].

Optimal control problems with (single) goal functionals were investigated in [6, 40, 5, 50, 52, 43]. The $p$-Laplacian and a posteriori error estimates were considered in [36, 12, 20, 13], and, more specifically, for goal functional evaluations, we refer to [34, 44, 25]. To estimate goal functionals, we adopt the dual-weighted residual (DWR) method [7, 8] in which an adjoint problem is solved to obtain (local) sensitivity measures that are used for mesh refinement. As is well-known, using a gradient-based approach for the numerical solution of optimal control problems, the same adjoint problem as for the DWR error estimator can be employed. For this reason, it is natural to combine gradient-based optimization with adjoint-based error estimation.

We are specifically interested in an extended DWR version in which the discretization and (linear/nonlinear) iteration error are balanced [39, 44, 37]. As localization technique we employ integration by parts as done in [8] or, for residual based error estimates, in [49]. The extension of [44] to multiple goal functionals was recently undertaken in [25].

Three major aims constitute the main contents of this paper: first, the design of a framework for goal-oriented error estimation for optimal control subject to a nonlinear PDE and balancing the discretization and nonlinear iteration error (Section 3). From the optimization point of view, we carefully revisit the important elements for the DWR estimator for optimization problems. The main result in this respect is the a posteriori error representation for the reduced optimal control system for an abstract problem formulation. The second aim is the extension to the simultaneous control of multiple goal functionals (Section 4). As a third goal, based on our theoretical developments, we carefully design an adaptive solution algorithm (Section 5). The performance of our algorithms are investigated in terms of the usual quality measures of convergence behavior and effectivity indices in Section 6. The latter one measures the quality of our proposed error estimator in comparison to (known) true errors, which are computed on sufficiently refined meshes.

We summarize the outline of this work as follows: In Section 2, the problem setting is introduced. Next, in Section 3, the dual-weighted residual method for the reduced optimization problem is formulated. The multi-goal approach is then introduced in Section 4. Our algorithmic developments to solve the multiple goal-functional optimal control problem are derived in Section 5. In Section 6, we present several numerical examples that demonstrate the performance of our approach. Therein, we study different Tikhonov regularization parameters, we perform mesh refinement studies, and consider different goal functionals. In Section 7, we summarize the key outcomes of this work.

2 The Optimal Control Problem

In this section, we define an abstract problem formulation and collect some properties that we will rely on when deriving the a posteriori error estimates.

2.1 The Abstract Problem Formulation

Let $U$ and $Q$ be Banach spaces. We would like to find a control $\overline{q}\in Q$ and an associated state $\overline{u}\in U$ such that the pair $(\overline{u},\overline{q})$ is a local minimizer of some given cost functional $J(u,q)\colon U\times Q\to\mathbb{R}$, where $u$ and $q$ have to fulfill the so called state equation $A(u,q)=0$ with nonlinear differential operator $A$ acting between Sobolev spaces. More precisely, the arising PDE-constrained optimization problem reads as follows:

[TABLE]

for some operator $A\colon U\times Q\mapsto V^{*}$, where $V^{*}$ denotes the dual space of some Banach space $V$. For the theoretical findings in this paper, we assume that, for each $q\in Q$, the PDE is uniquely solvable. More precisely, we assume the following:

Assumption 1**.**

Let there exist a unique mapping $S\colon Q\mapsto U$ which is implicitly defined by

[TABLE]

Moreover, we assume that $S$ is twice continuously Fréchet differentiable.

Without further mention, we also assume the existence of a at least one global minimizer for Problem (1). For instance, we refer to [47] for general theorems on existence of solutions for problems with linear and semilinear state equations. Moreover, let $A$ and $J$ be smooth enough for all operations occurring in the next Section.

With the help of the so called control-to-state mapping $S$, we reformulate (1) as an unconstrained optimization problem

[TABLE]

where $j(q):=J(S(q),q)$. Here, we will also assume sufficient smoothness in order to derive all further estimates.

2.2 First Order Necessary Optimality Conditions

It is clear, that under our implicit smoothness assumptions, the first order necessary optimality conditions for a locally optimal control $\bar{q}\in Q$ for Problem (2.1) are given by

[TABLE]

For completeness and further use, we rewrite these conditions for the non-reduced formulation with the help of the well-known Lagrange approach. We define the Lagrangian $\mathcal{L}\colon U\times Q\times V\mapsto\mathbb{R}$ for this problem as follows

[TABLE]

To shorten notation, we consider the abbreviation $B^{\prime}_{\zeta}:=\frac{\partial}{\partial\zeta}B$ for the partial derivatives of some operator $B$. The first order necessary optimality conditions for (1) are then given by

[TABLE]

Moreover, $\bar{u}=S\bar{q}$ denotes the optimal state associated with $\bar{q}$, and $\bar{z}=(S^{\prime}(\bar{q}))^{*}J^{\prime}_{u}(\bar{u},\bar{q})$ the associated adjoint state. In order for the Newton algorithm to work, and for the error estimator we need the following assumption.

Assumption 2**.**

We assume that $A^{\prime}_{u}=\mathcal{L}^{\prime\prime}_{uz}$ is invertible.

2.3 An Example: the Regularized $p$-Laplacian and Tracking-type Cost Functional

Let us finish this section by defining $A$ for a concrete example (i.e., a PDE) that motivates our numerical studies. To this end, a (regularized) $p$-Laplace equation for $p\neq 2$ is considered, even though, it does not necessarily fit into the theory setting. For details, we refer to [22, 34, 46] and the references therein regarding the (regularization of) the $p$-Laplace equation. We consider the following setting: Let $\Omega\subset\mathbb{R}^{d}$ be open and bounded with $C^{1}$ boundary, and let $p\in(\frac{2d}{2+d},\infty)$. Then we define

[TABLE]

by the identity

[TABLE]

for $u,v\in U:=W^{1,p}_{0}(\Omega),$ $f\in V^{*}$, where $\langle\cdot,\cdot\rangle$ is the usual notation for duality pairings. Note that in this example, we have $U=V$.

Let $u\in U$ be the state, and $q\in Q$, e.g., $Q=L^{2}(\Omega)$, be the control variable. Then our optimal control problem is given by

[TABLE]

with the tracking-type cost functional

[TABLE]

with $\alpha>0$ and given $f\in U^{*}$, $\bar{u}^{d}\in L^{2}(\Omega)$ and $\bar{q}^{d}\in L^{2}(\Omega)$.

3 The Dual Weighted Residual Method for the Reduced System

We now formulate the DWR method for the reduced optimal control system and develop a posteriori error estimators. The presentation is kept as general as possible so that the extension to multiple goal functionals outlined in Section 4 can easily be incorporated. Firstly, we briefly outline the important elements of the discretization.

3.1 Discretization

The method of choice, which will be used in the numerical examples, is the finite element method [19, 11, 31]. However, the algorithms presented in this work can also be adapted to other discretization techniques where adaptivity can be accomplished, like isogeometric analysis, the virtual element method, or finite cell methods. For the spaces $U_{h}=V_{h}$, we use continuous tensor product finite elements $Q_{c}^{r}$ ;see, for instance, [19]. For $Q_{h}$ we use discontinuous tensor product finite elements $Q_{DG}^{r}$. Let $\mathcal{T}_{h}$ be a subdivision (triangulation) of the domain $\Omega$ into quadrilateral elements such that $\bigcup_{K\in\mathcal{T}_{h}}\overline{K}=\overline{\Omega}$ and $K\cap K^{\prime}=\emptyset$ for all $K,K^{\prime}\in\mathcal{T}_{h}$ where $K\neq K^{\prime}$. Furthermore, let $\psi_{K}$ be a multilinear mapping from the reference element $\hat{K}=(0,1)^{d}$ to the element $K\in\mathcal{T}_{h}$. We define the space $Q_{DG}^{r}$ as

[TABLE]

with $Q_{r}(K):=\{v_{|\hat{K}}\circ\psi_{K}^{-1}:\,v(\hat{x})=\prod_{i=1}^{d}(\sum_{\beta=0}^{r}c_{\beta,i}\hat{x}_{i}^{\beta}),\,c_{\beta,i}\in\mathbb{R}\}$. The use of these finite dimensional spaces leads to a conforming discretization for Example 2.3. We point out that the conforming discretization is needed in order to keep Theorem 3.5 valid. The discretized abstract model problem reads as follows: Find $u_{h}\in U_{h}$ and $q_{h}\in Q_{h}$ such that they are a local solution pair of

[TABLE]

Assumption 3**.**

There exists a unique discrete mapping $S_{h}\colon Q_{h}\mapsto U_{h}$, which is implicitly defined by

[TABLE]

As for its continuous counterpart, we assume that it is twice continuously Fréchet differentiable.

Using the discrete mapping $S_{h}$, we can reformulate Problem (7) as the unconstrained optimization problem: Find $q_{h}\in Q_{h}$ such that it solves

[TABLE]

Similar to Section 2.2, we also provide the discrete version of the first order necessary optimality conditions. If $\bar{q}_{h}\in Q_{h}$ is a local solution, then these conditions are given by

[TABLE]

We will also use the non-reduced formulation with the help of the Lagrange-approach, with

[TABLE]

The discrete first order necessary optimality conditions for (7) are then given by

[TABLE]

where $\bar{u}_{h}=S_{h}(\bar{q}_{h})$ and $\bar{z}_{h}=(S_{h}^{\prime}(\bar{q}_{h}))^{*}J^{\prime}_{u}(\bar{u}_{h},\bar{q}_{h})$.

3.2 Error Representation for the Reduced System

We are now interested in an error estimator for a quantity of interest $I\colon U\times Q\mapsto\mathbb{R}$. Let $\overline{q}$ be an optimal control of Problem (2.1) with associated optimal state $\bar{u}=S\bar{(}q)$. While we are interested in $I(\overline{u},\overline{q})$, we can only compute an approximation $I(\tilde{u}_{h},\tilde{q}_{h})$ of this value. Note that we assume, for most of what follows, that $\tilde{u}_{h}:=S_{h}(\tilde{q}_{h})$ is exactly solved by means of the solution operator $S_{h}$ for the discrete state equation, cf. Section 3.1. To estimate this error, we apply the previously mentioned DWR method (e.g., [8]) to the first order optimality conditions of our reduced system.

Defining $i(q):=I(S(q),q)$ as well as $i_{h}(q):=I(S_{h}(q),q)$, the error between $I(S(\bar{q}),\bar{q})$ and $I(S_{h}(\tilde{q}_{h}),\tilde{q}_{h})$ can be split into

[TABLE]

Therefore, $i_{h}$ still corresponds to our "true" quantity of interest, but computed with the discrete solutions $\tilde{q}_{h}$ and $S_{h}(\tilde{q}_{h})$. We start by estimating the first part of the error, which actually has a practical relevance: if some approximate control $\tilde{q}_{h}$ is computed and applied in a practical situation, then the corresponding physical system will produce a "true" state $\tilde{u}:=S(\tilde{q}_{h})$ instead of an approximation $\tilde{u}_{h}=S_{h}(\tilde{q}_{h})$.

As a first result, we formulate a theoretical error estimator, where we need the adjoint problem to the first order optimality conditions, which is given by: Find $\overline{p}\in Q$ such that

[TABLE]

Assumption 4**.**

We assume that (11) has a unique solution.

Theorem 3.1** (Error Representation for Reduced System).**

Let us assume that $j\in\mathcal{C}^{4}(Q,\mathbb{R})$ and $i\in\mathcal{C}^{3}(Q,\mathbb{R})$. If $\overline{q}$ solves (2.1) and $\overline{p}$ solves (11) for $\overline{q}\in Q$, then, for arbitrary fixed $\tilde{q}_{h}\in Q$ and $\tilde{p}_{h}\in Q$, we find:

[TABLE]

where

[TABLE]

and the remainder term satisfies

[TABLE]

with $e:=\overline{q}-\tilde{q}_{h}$ and $e^{*}:=\overline{p}-\tilde{p}_{h}$.

Proof.

The proof follows the same idea as in [44, 25] but is stated for completeness of presentation. Define $e$ and $e^{*}$ as above and let $\overline{x}$, $x$, $\tilde{x}_{h}$ be defined as $\overline{x}:=(\overline{q},\overline{p})$, ${x}:=({q},{p})$, $\tilde{x}_{h}:=(\tilde{q}_{h},\tilde{p}_{h}),$ as well as $m(x):=i(q)+j^{\prime}(q)(p)$. Furthermore, let $e_{x}$ be defined $e_{x}:=\overline{x}-\tilde{x}_{h}$. By the fundamental theorem of calculus as well as the trapezoidal rule, we observe that

[TABLE]

By carefully inspecting $\frac{1}{2}\int\limits_{0}^{1}m^{\prime\prime\prime}(\tilde{x}_{h}+se_{x})(e_{x},e_{x},e_{x})s(s-1)\text{ d}s$, it follows that it coincides with $\mathcal{R}^{(3)}$. Additionally, we can deduce that

[TABLE]

due to (3) and (11). Combining (13) and (14) results in the following identity

[TABLE]

Therefore, using again (3) as well as (12), we get

[TABLE]

where we have applied (15). This proves the theorem after verifying that $m^{\prime}(\tilde{x}_{h})(e_{x})=\rho(\tilde{q}_{h})(\overline{p}-\tilde{p}_{h})+\rho^{*}(\tilde{q}_{h},\tilde{p}_{h})(\overline{q}-\tilde{q}_{h})$. ∎

Remark 3.2**.**

One objective of this representation, in addition to the fact that for instance $\tilde{u}=S(\tilde{q}_{h})$ is not readily available exactly, is to obtain indicators for local adaptivity. By inspecting the primal part of the error estimator $\rho(\tilde{q}_{h})(\overline{p}-\tilde{p}_{h})$, we observe that

[TABLE]

Since it is not clear how to localize $S^{\prime}(\tilde{q}_{h})(\overline{p}-\tilde{p}_{h})$, we do not follow this path to compute the error indicators, but prove a localizable error estimator in a similar fashion in Theorem 3.5, which makes use of (5) as well.

For another idea, we consider the adjoint problem to the first order optimality conditions for the Lagrangian defined in (4): Find $(\overline{v},\overline{p_{2}},\overline{y})\in U\times Q\times V$ such that

[TABLE]

where the argument in the partial derivatives is always given by $(\overline{u},\overline{q},\overline{z})$.

Assumption 5**.**

We assume that (16) has a unique solution.

In order to obtain the variables $\overline{v}$ and $\overline{y}$ with the help of the solution of the reduced adjoint problem (11), the following lemma is useful.

Lemma 3.3**.**

If $\bar{q}\in Q$ with associated state $\bar{u}=S(\bar{q})$ is a local solution of (2.1), and $\bar{p}$ solves (11), then $\bar{v}=S^{\prime}(\bar{q})\bar{p}$, $\bar{p}_{2}=\bar{p}$, and $\bar{y}$ given by (21) solve (16).

Proof.

Let $\mathfrak{p}\in Q$ be arbitrary. Using the definition of the reduced functionals, we obtain

[TABLE]

and

[TABLE]

Furthermore, with the definition of the solution operator, we obtain from (2) that

[TABLE]

and

[TABLE]

By subtracting (19) from (17), it follows that

[TABLE]

Further, from (18) we get

[TABLE]

Thus $\overline{p}_{2}=\overline{p}$ and $\overline{v}=S^{\prime}(\overline{q})\overline{p}$ satisfy the third line in (16).

To proceed, we note that $\overline{q}$, $\overline{u}=S(\overline{q})$ and $\overline{z}$ solves (5), thus we have that $\mathcal{L}^{\prime}_{u}(S(\overline{q}),\overline{q},\overline{z})=0$. This leads to

[TABLE]

Now, we define $\overline{y}$ by the first line of (16), we get

[TABLE]

With this, we can rewrite (20) as

[TABLE]

Now, we can use the definition of $\overline{p}$, $\mathcal{L}^{\prime\prime}_{uz}=(\mathcal{L}^{\prime\prime}_{zu})*,\mathcal{L}^{\prime\prime}_{zq}=(\mathcal{L}^{\prime\prime}_{qz})*$, the formula for $S^{\prime}(\overline{q})$ and the representation of $i^{\prime}(\overline{q})$ to get

[TABLE]

and the second line in (16) follows. ∎

Lemma 3.3 allows to obtain $\overline{p}=\bar{p}_{2}$ by solving the reduced adjoint equation (11). Then, $\overline{v}$ can be computed by solving the tangent equation

[TABLE]

which is the last row of (16). Using this solution, we can deduce $\overline{y}$ from the first row of (16).

An analogue to (16) on the discrete level is given by: Find $(\tilde{v}_{h},\tilde{p}_{h},\tilde{y}_{h})\in U_{h}\times Q_{h}\times V_{h}$ such that

[TABLE]

where the arguments in the partial derivatives are given by $(\tilde{u}_{h},\tilde{q}_{h},\tilde{z}_{h})$.

Remark 3.4**.**

If (22) is considered at the linearization point $\bar{q}_{h},\bar{u}_{h},\bar{z}_{h}$, then Lemma 3.3 holds also true for the discrete problem, i.e. if $\overline{p}_{h}\in Q_{h}$ solves

[TABLE]

then $\tilde{p}_{h}=\bar{p}_{h}$. This can be shown by the same proof replacing $S$ by $S_{h}$.

Similar as explained above, the variables $\tilde{v}_{h}$ and $\tilde{y}_{h}$ can be deduced from the knowledge of $\bar{p}_{h}$ and the discrete version of Lemma 3.3.

Theorem 3.5** (Localizable Error Representation for Reduced System).**

Let us assume that $j\in\mathcal{C}^{4}(Q,\mathbb{R})$ and $i\in\mathcal{C}^{3}(Q,\mathbb{R})$. Let $\overline{q}$ be a local solution of (2.1), with $\overline{\xi}=(\overline{u},\overline{q},\overline{z})$ the corresponding KKT-triplet given by (5), and let the triple $\overline{\xi}^{*}=(\bar{v},\overline{p},\bar{y})\in U\times Q\times V$ solve (16). Moreover, let $\tilde{q}_{h}\in Q_{h}$ be an arbitrary fixed discrete control, and let $\tilde{\xi}_{h}^{*}=(\tilde{p}_{h},\tilde{v}_{h},\tilde{y}_{h})$ be the solution to (10) and the first and last row of (22) at the linearization point $\tilde{\xi}_{h}=(\tilde{u}_{h},\tilde{q}_{h},\tilde{z}_{h})$ with $\tilde{u}_{h}=S_{h}(\tilde{q}_{h})$ and $\tilde{z}_{h}=S^{\prime}_{h}(\tilde{q}_{h})^{*}J^{\prime}_{u}(\tilde{u}_{h},\tilde{q}_{h})$. Then we have the error representation

[TABLE]

where

[TABLE]

and the remainder term

[TABLE]

with $\tilde{e}_{\xi}=\overline{\xi}-\tilde{\xi}_{h}$, $\tilde{e}_{\xi}^{*}=\overline{\xi}^{*}-\tilde{\xi}_{h}^{*}$.

Proof.

The proof follows a similar structure as the proof of Theorem 3.1. Let $\overline{x}:=(\overline{\xi},\overline{\xi}^{*})$, $\tilde{x}_{h}:=(\tilde{\xi}_{h},\tilde{\xi}^{*}_{h})$. For $x=(\xi,\xi^{*})=(u,q,z,\xi^{*})$ we define $\mathcal{M}(x):=I(\xi)+\mathcal{L}^{\prime}(\xi)(\xi^{*})=I(u,q)+\mathcal{L}^{\prime}(\xi)(\xi^{*})$. It holds that

[TABLE]

where $e_{x}=\overline{x}-\tilde{x}_{h}$. By carefully inspecting $\mathcal{M}^{\prime\prime\prime}(\tilde{x}_{h}+se_{x})(e_{x},e_{x},e_{x})$ it follows that

[TABLE]

since $\mathcal{M}^{\prime\prime}_{\xi^{*}\xi^{*}}=0$ and $\mathcal{M}^{\prime}_{\xi^{*}}(\tilde{x}_{h}+se_{x})(e_{x})=\mathcal{L}^{\prime}(\tilde{\xi}_{h}+s(\tilde{e}_{\xi}))(\tilde{e}_{\xi}^{*})$. Thus, (25) gives

[TABLE]

For the part $\mathcal{M}^{\prime}(\overline{x})(e_{x})$ of (25), we can deduce that

[TABLE]

since $\overline{\xi}^{*}$ solves (16) and $\overline{\xi}$ solves (5). Finally, relation (25) reduces to the following identity

[TABLE]

Therefore, we get

[TABLE]

Furthermore, we can deduce that $\mathcal{M}^{\prime}(\tilde{x}_{h})(e_{x})=\mathcal{L}^{\prime}(\tilde{\xi}_{h})(\tilde{e}_{\xi}^{*})+I^{\prime}(\tilde{\xi}_{h})(\tilde{e}_{\xi})+\mathcal{L}^{\prime\prime}(\tilde{\xi}_{h})(\tilde{e}_{\xi},\tilde{\xi}^{*}_{h})$. Gathering the results from above, we obtain, noting that $\tilde{\xi}_{h}=(\tilde{u}_{h},\tilde{q}_{h},\tilde{z}_{h})=(S_{h}(\tilde{q}_{h}),\tilde{q}_{h},\tilde{z}_{h})$

[TABLE]

Straightforward calculations show

[TABLE]

and

[TABLE]

∎

Let us end this section with some further observations.

Remark 3.6**.**

Note that if $\tilde{q}_{h}=\bar{q}_{h}$, then $(\tilde{v}_{h},\tilde{p}_{h}=\bar{p}_{h},\tilde{y}_{h})$ in fact solve (10), cf. Remark 3.4, and consequently $j^{\prime}_{h}(\tilde{q}_{h})(\tilde{p}_{h})=0$.

Remark 3.7**.**

From numerical experiments for the regularized $p$-Laplacian computed in [27], we can deduce that $\mathcal{R}^{(3)}$ can be neglected on sufficiently refined meshes.

An identity also observed in [51], is the following:

Proposition 3.1**.**

If $I=J$ and $j^{\prime\prime}(q)$ is injective, then we have $(\overline{v},\overline{p},\overline{y})=(0,0,\overline{z})$.

Proof.

Since $J$ is the cost functional and $(\overline{u},\overline{q})$ is a local minimizer of our optimization problem the first order necessary condition is given by $j^{\prime}(\overline{q})=0$. Therefore the adjoint equation reads as

[TABLE]

If $j^{\prime\prime}(\overline{q})$ is injective, then $\overline{p}=0$. From the tangent equation

[TABLE]

we can deduce that $\overline{v}=0$. Finally the optimality system reduces to

[TABLE]

From this follows that $\overline{y}=\overline{z}$, which completes the proof. ∎

3.3 The Parts of the Error Estimator

We now briefly discuss the two main parts of the error estimator:

[TABLE]

where the first part refers to the iteration error, and the second term denotes the discretization error to be defined in the following. We recall that $\eta_{h,k}^{(2)}$ is designed to estimate $i(\overline{q})-i_{h}(\tilde{q}_{h})$ given in (24).

The iteration error estimator

The iteration error estimator

[TABLE]

can be used as stopping rule for the nonlinear solver like for Newton’s method as in [44, 25, 25] and Algorithm 1 presented in Section 5.

The discretization error estimator

Of course the exact solution of the optimal control problem in formula (24) are not known. They can either be replaced by a (patch-wise) higher order polynomial interpolation or by approximations on enriched spaces [8, 4].

The discretization error estimator using the solutions $(u_{h}^{(2)},q_{h}^{(2)},z_{h}^{(2)})$ and $(v_{h}^{(2)},p_{h}^{(2)},y_{h}^{(2)})$ on enriched spaces reads as

[TABLE]

The replacement is justified if a strengthened saturation assumption is fulfilled as shown in [27] for both the nonlinear state equation and the goal functionals.

We briefly recall that the localization can be performed in three ways: classical integration by parts yielding the strong problem formulation [8], a filtering approach employing the weak problem formulation [10], or a partition-of-unity using again the weak form of the problem [45]. All three techniques are analyzed (theoretically and computationally) with respect to their effectivity in [45]. In the theoretical analysis, a discrete version of Lemma 3.3 is necessary to justify that $(v_{h}^{(2)},p_{h}^{(2)},y_{h}^{(2)})$ is indeed a solution in the enriched spaces.

4 Extension to Multiple Goal Functionals

In Section 3, we discussed how the DWR method works for one functional. However, for some problems, several functional evaluations would be of interest. Let us consider $N$ goal functionals $I_{1},I_{2},\ldots,I_{N}$ for some $N\in\mathbb{N}$. One possibility would be to compute the error estimators separately as described in Section 3. However, we would have to solve the adjoint problem $N$ times, leading to high computational cost. There are several ways to tackle this problem as for example discussed in [33, 32, 42, 1] and more recently in [35, 28, 25, 26, 27].

Adopting the techniques presented in [25], we try to combine the functionals to one, and apply the DWR method for one functional to it. In the following section, we consider $\overline{u}$, $\overline{q}$ as the solution of (1), and $\tilde{u}_{h}$, $\tilde{q}_{h}$ as some approximations. To construct the combination, we introduce a so called error weighting function:

Definition 4.1** (Error weighting function [25]).**

Let $M\subseteq\mathbb{R}^{N}$. We say that $\mathfrak{E}:(\mathbb{R}^{+}_{0})^{N}\times M\mapsto\mathbb{R}^{+}_{0}$ is an error-weighting function if $\mathfrak{E}(\cdot,m)\in\mathcal{C}^{1}((\mathbb{R}^{+}_{0})^{N},\mathbb{R}^{+}_{0})$ is strictly monotonically increasing in each component and $\mathfrak{E}(0,m)=0$ for all $m\in M$.

As in [25], let $\vec{I}(\cdot):=(I_{1}(\cdot),I_{2}(\cdot),\ldots,I_{N}(\cdot))$ mapping from $\bigcap_{i=1}^{N}\mathcal{D}(I_{i})\subset U\times Q\mapsto\mathbb{R}^{N}$. Furthermore, we define $|\cdot|_{N}:\mathbb{R}^{N}\mapsto(\mathbb{R}_{0}^{+})^{N}$ as the component-wise absolute value. This allows us to construct the error function $I_{\mathfrak{E}}$ as follows:

[TABLE]

Remark 4.2**.**

The error functional $\tilde{I}_{\mathfrak{E}}$ is constructed in a way, that avoids error cancellation between two or more functionals. For a more detailed discussion, we refer the reader to [25, 27].

Remark 4.3**.**

The quantity (30) is not computable, since it depends on $\vec{I}(\overline{u},\overline{q})$, which is not known. However, we can use a higher order polynomial approximation to approximate this quantity, as done in [33, 25, 27], where consequences of the replacement are discussed in [27].

The resulting error weighting functional is given by

[TABLE]

where $u_{h}^{(2)}$, $q_{h}^{(2)}$ denote the solutions on enriched finite element spaces.

Remark 4.4**.**

We notice that, for the choice $\mathfrak{E}(x,m):=\sum_{\ell=1}^{N}\frac{x_{\ell}}{|m_{\ell}|}$, we obtain the same combined functional as in [28] up to sign. The same holds for [33, 32] in the case of linear problems. This choice is used in our numerical examples.

Remark 4.5**.**

Finally, the method explained in Section 3 is applied to $I_{\mathfrak{E}}$ instead of $I$ to achieve a control of the errors in all functionals at once, as algorithmically illustrated in Section 5.

5 Algorithmic Details

In this section, we briefly recapitulate the algorithmic techniques to solve the optimal control problem with multiple goal functionals that we have outlined in the previous sections. The algorithms for the forward problem including multiple goal functionals evaluations were derived in [25]. Therein, the goal functionals were estimated using the DWR method (thus an adjoint approach). Hence, the extension to optimal control using a gradient-based approach is straightforward. The implementation of the following algorithms is done in the open-source library DOpElib [23, 30]. For a general overview of optimization algorithms, we refer to [41, 38]. First, we present the reduced Newton method described in Algorithm 1.

Remark 5.1**.**

The parameter $\gamma$ is chosen as $10^{-2}$ in the numerical experiments.

Remark 5.2**.**

In [30], we specifically used DOpE::ReducedNewtonAlgorithm::ReducedNewtonLineSearch to obtain the line search parameter $\alpha^{k}$.

Remark 5.3**.**

*The arising linear problems $j_{h}^{\prime\prime}(q^{l,k}_{h})(v_{h},p^{l,k}_{h})=j_{h}^{\prime}(q^{l,k}_{h})(v_{h})$ and

$j_{h}^{\prime\prime}(q^{l,k}_{h})(v_{h},p^{l,k}_{h})=(i_{\mathfrak{E},h}^{(k)})^{\prime}(q^{l,k}_{h})(v_{h})$ were solved by using the algorithm

DOpE::ReducedNewtonAlgorithm::SolveReducedLinearSystem implemented in [30].*

With the help of Algorithm 1, we can now state the final Algorithm 2 used in this paper.

Remark 5.4**.**

In Algorithm 2 in Step 8, we use Dörfler marking with $\theta=0.5$ as marking strategy [24].

Remark 5.5**.**

The reduced discrete cost functional $j_{h}$ on the space $Q_{h}^{l,(2)}$ is constructed by means of the corresponding discrete solution operator on the enriched space.

Remark 5.6**.**

To solve the linear systems arising form the forward state equation, we use the sparse direct solver UMFPACK [21].

6 Numerical examples

In the current section, we provide some numerical examples demonstrating the performance of the theoretical arguments and algorithms developed previously. The implementation is done in DOpElib [23, 30] using the finite elements from deal.II [3, 2]. However, large parts of the programming are new . For this reason, we first present a linear example with a single goal functional, which has been already studied in the literature. In the second example, we then consider the $p$-Laplacian and again the case of a single goal functional. In Example 3, we study several nonlinear goal functionals that are simultaneously controlled. The quality of our results will be measured by effectivity index which is given by

[TABLE]

whereas the primal and adjoint effectivity indices are defined by

[TABLE]

and

[TABLE]

Notice that we do not apply the absolute value to the contributions. Hence, we also estimate the sign of the error.

6.1 Example 1: linear Laplacian, single goal functional

In this first numerical test, we consider a standard linear example, which is implemented, for instance, in DOpElib[23, 30][OPT/StatPDE/Example1, Section 6.1.1]. The main purpose is to validate our novel programming code against known findings. The domain is $\Omega:=(0,1)^{2}$. The right-hand side forces of the PDE are f(x,y):=\big{(}20\pi^{2}\text{sin}(4\pi x)-\alpha^{-1}\text{sin}(\pi x)\big{)}\text{sin}(2\pi y). The given control is $q^{d}:=0$, and the desired state is u^{d}:=\big{(}5\pi^{2}\text{sin}(\pi x)+\text{sin}(4\pi x)\big{)}\text{sin}(2\pi y). The regularization is chosen as $\alpha=10^{-2}$.

The problem statement is as follows: Find $(\overline{u},\overline{q})\in H^{1}_{0}(\Omega)\times L^{2}(\Omega)$ such that it is a minimizer of

[TABLE]

with the constraints

[TABLE]

The exact minimizer of the problem is known, and given by $\overline{u}(x,y)=\text{sin}(4\pi x)\text{sin}(2\pi y)$ and $\overline{q}(x,y)=\alpha^{-1}\text{sin}(\pi x)\text{sin}(2\pi y)$. First of all, we use $I=J$, so the cost functional as quantity of interest. Here, the exact value is given by J(\overline{u},\overline{q})=\frac{1}{8}\big{(}25\pi^{4}+\alpha^{-1}\big{)}.

In the Figures 1 and 2, the effectivity index $I_{\text{eff}}$ and the error are both shown against the number of degrees of freedom (DOFs). For the single error parts, primal and adjoint estimators, the effectivity indices show significant differences from the asymptotically expected value. Combining both parts, then yields an optimal $I_{\text{eff}}=1$. Convergence of adaptive and uniform mesh refinement are shown in Figure 2.

1e-050.00010.0010.010.1110100101001000100001000001e+06DOFsError in $I_{\mathfrak{E}}$ (adp.)Estimated ErrorError in $I_{\mathfrak{E}}$ (uni.)$\mathcal{O}(\text{DOFs}^{-1})$

1e-050.00010.0010.010.1110100101001000100001000001e+06DOFsError in $I_{\mathfrak{E}}$ (adp.)Estimated ErrorError in $I_{\mathfrak{E}}$ (uni.)$\mathcal{O}(\text{DOFs}^{-1})$

In this second part of the example, we apply the method to a quantity that is different to the cost functional. We are interested in $I(u,q):=\|u\|_{L^{1}(\Omega)}$. The exact value is given by $I(\overline{u},\overline{q})=4\pi^{-2}$. The corresponding numerical findings are displayed in the Figures 3 and 4. We observe excellent effectivity indices in Figure 3. Optimal convergence rates also in comparison with uniform mesh refinement are observed in Figure  4.

-2.5-2-1.5-1-0.500.511.522.5101001000100001000001e+06DOFsI_{eff}$$I_{effp}$$I_{effa}1

1e-050.00010.0010.010.1110101001000100001000001e+06DOFsError in $I_{\mathfrak{E}}$ (adp.)Estimated ErrorError in $I_{\mathfrak{E}}$ (uni.)$\mathcal{O}(\text{DOFs}^{-1})$

6.2 Example 2: $p$-Laplacian, single goal functional

We now proceed to nonlinear state equations and consider the example PDE provided in Section 2.3. Here, $\Omega$ (and the initial mesh) and $u^{d}$ are given in Figure 5. Furthermore, $q^{d}=1$, $p=4$, $\varepsilon=1$ and $f=0$. In particular, we investigate various regularization parameters $\alpha$. The goal functional $I(u,q)$ is given by $I(u,q):=\int_{\Omega}u(x)^{2}q(x)^{2}dx$.

In Table 1, we obtain, for $\alpha=0.01,\ldots,10$, effectivity indices in the range of $0.88$ to $1.30$, which are excellent findings in view of the nonlinear behavior of the state equation and the geometric singularities introduced by the domain. In the case of $\alpha=0.1$, we obtain a $I_{\text{eff}}$ in the range of $0.36$ to $0.79$, which might be affected by cancellation effects from adding the different contributions to the error estimator. The exact value of the functionals was approximated by one additional $p$ and $h$ refinement, and is given in the last line of Table 1 corresponding to $l=\infty$, with additional information on the number of DOFs used to compute this values.

In the Figures 6 and 7, the final meshes for different $\alpha$ are shown. For $\alpha=10^{-2}$, we observe very localized mesh refinement, while, for larger $\alpha$, the mesh is still locally refined, but in a somewhat uniform behavior. The states and controls on these final meshes are displayed in the Figures 8 and 9.

6.3 Example 3: $p$-Laplacian, multiple goal functionals

In this third example, we proceed to multiple goal functionals. The setup is the same as in Example 2, but with a single $\alpha=0.01$ and multiple goal functionals:

• •

$I_{1}(u,q)=\frac{1}{2}\int_{\Omega}(u-\overline{u}^{d})^{2}dx\approx 1.15760$,

• •

$I_{2}(u,q)=\frac{1}{2}\int_{\Omega}(q-\overline{q}^{d})^{2}dx\approx 21.3305$ ,

• •

$I_{3}(u,q)=\int_{([4,5]\times\mathbb{R})\cap\Omega}udx\approx-0.236288$ ,

• •

$I_{4}(u,q)=\int_{[1,\frac{25}{4}]\times[2,\frac{5}{2}]}qdx\approx 0.328042$ ,

• •

$I_{5}(u,q)=\frac{1}{2}\int_{\Omega}u^{2}q^{2}dx\approx 0.231615$.

The geometry alongside with the goal functionals $I_{3}$ and $I_{4}$ is illustrated in Figure 10.

0.60.70.80.911.11.21.31.4100100010000100000DOFsI_{eff}$$I_{effp}$$I_{effa}1

0.00010.0010.010.111010010001001000100001000001e+06DOFsError in $I_{\mathfrak{E}}$ (adp.)Estimated ErrorError in $I_{\mathfrak{E}}$ (uni.)$\mathcal{O}(\text{DOFs}^{-1})$

1e-050.00010.0010.010.11101001000100100010000100000DOFsI_{1}$$I_{2}$$I_{3}$$I_{4}$$I_{5}$$I_{\mathfrak{E}}

1e-050.00010.0010.010.11101001000100100010000100000DOFs\eta_{h}^{(2)}$$\rho_{u}$$\rho_{p}$$\rho_{z}$$\rho_{v}$$\rho_{q}$$\rho_{y}

The reference values are computed on a fine grid ($716\,792$ DOFs for $q$ + $730\,199$ DOFs for $u$) , which is obtained by 12 adaptive refinements for $J_{\mathfrak{E}}$ followed by two uniform h-refinements and one uniform p-refinement. Our findings are displayed in the Figures 11, 12, 13 and 14. In Figure 11 the calculated effectivity indices are excellent in view of the nonlinearities of the domain, state equation and multiple goal functionals. Curves of the errors and estimators are shown in the Figures 12, 13 and 14. Here, the combined functional (as expected) bounds all single functionals. In Figure 12, we observe that adaptive refinement pays off in delivering the same error as uniform mesh refinement, but with a lower computational cost. The convergence rates are the same, which lies in the fact that the control is chosen in such a way that a sufficiently smooth final solution is obtained. Finally, we compare the adaptive stopping rule used in Algorithm 1 with the standard stopping rule, which is used in the DOpElib [23, 30] algorithm DOpE::ReducedNewtonAlgorithm::Solve with the absolute residual nonlinear_global_tol = 1.e-7 and relative residual nonlinear_tol= 8.e-5. Since the discretization error estimate is not given for $l=0$ in Algorithm 1, we use $\eta_{h}^{l-1}=10^{-5}$.

We abbreviate the first algorithm with AN (Adaptive Newton) and the second algorithm with FN (Full Newton). In Table 2, we monitor that the $I_{\text{eff}}$ show a pretty similar behavior even for the adaptive stopping rule. Even though we need $1-3$ iterations in case of the adaptive stopping rule compared to $2-17$ iterations for the standard stopping rule, which is illustrated in Table 2 as well. Furthermore, we want to notice that the refined meshes for both algorithms coincide exactly up to $l=7$. For $l=8$, it is exactly one element, which is refined additionally in the case of FN. If we compare the corrected effectivity indices

[TABLE]

for the two stopping rules, we observe that they coincide even more after the correction.

In Table 3, the comparison between the estimated iteration error and the real error in the combined functional is shown. The ratio between $\eta_{k}$ and the error mimics the choice of $\gamma$ in Algorithm 1 for our adaptive stopping rule, whereas there is almost no correlation for the standard stopping rule.

7 Conclusions

In this work, we developed a novel a posteriori multiple goal-oriented error estimation for optimal control problems subject to a nonlinear state equation. The error estimator also serves for balancing the discretization and nonlinear iteration error. The overall optimization problem is solved via a reduced approach in which the state equation is eliminated by a control-to-state solution operator. In Section 3.2, the theoretical results yield an a posteriori estimate for a single goal functional. The extension to multiple goal functionals was made in Section 4. Based on these theoretical aspects, the algorithmic details were worked out in the following section. Three numerical examples were investigated. In the first example, our approach was tested against configurations known in the literature. The Examples 2 and 3 are more advanced by considering the regularized $p$-Laplacian as nonlinear state equation. The main criterion whether the proposed error estimator works sufficiently well is given by the effectivity index. In the numerical examples, values around one were obtained. These are excellent findings in view of the challenging nature of the underlying problem configuration; namely domain (corner) singularities, quasi-linear state equations within an optimal control setting, and finally multiple nonlinear goal functionals. Ongoing work considers the extension to elasticity and more practical applications.

8 Acknowledgments

This work has been supported by the Austrian Science Fund (FWF) under the grant P 29181 ‘Goal-Oriented Error Control for Phase-Field Fracture Coupled to Multiphysics Problems’ and the DFG-SPP 1962 ‘Non-smooth and Complementarity-based Distributed Parameter Systems: Simulation and Hierarchical Optimization’ within the project ‘Optimizing Fracture Propagation Using a Phase-Field Approach’ under grant numbers NE1941/1-1 and WO1936/4-1.