---
abstract: |
    In this contribution we propose and rigorously analyze new variants of adaptive Trust-Region methods for parameter optimization with PDE constraints and bilateral parameter constraints. The approach employs successively enriched Reduced Basis surrogate models that are constructed during the outer optimization loop and used as model function for the Trust-Region method. Each Trust-Region sub-problem is solved with the projected BFGS method. Moreover, we propose a non-conforming dual (NCD) approach to improve the standard RB approximation of the optimality system. Rigorous improved a posteriori error bounds are derived and used to prove convergence of the resulting NCD-corrected adaptive Trust-Region Reduced Basis algorithm. Numerical experiments demonstrate that this approach enables to reduce the computational demand for large scale or multi-scale PDE constrained optimization problems significantly.

    **Keywords**: PDE constrained optimization, Trust-Region method, error analysis, Reduced Basis method, model order reduction, parametrized systems, large scale problems

    **AMS Mathematics Subject Classification**: 49M20, 49K20, 35J20, 65N30, 90C06
author:
- |
    Tim Keil$^\dag$, Luca Mechelli$^\ddag$, Mario Ohlberger$^\dag$,\
    Felix Schindler[^1], Stefan Volkwein[^2]
bibliography:
- 'bibliography.bib'
title: 'A non-conforming dual approach for adaptive Trust-Region Reduced Basis approximation of PDE-constrained optimization[^3] '
---

Introduction {#sec:introduction .unnumbered}
============

We are concerned with the development and rigorous analysis of novel efficient model order reduction methods for parameter optimization constrained by coercive variational state equations using the [*first optimize, then discretize*]{} approach. The methods are based on successive enrichment of the underlying reduced order models within the framework of Trust-Region optimization. Optimization problems constrained by partial differential equations (PDEs) arise in many fields of application in engineering and across all sciences. Examples of such problems include optimal (material) design or optimal control of processes and inverse problems, where parameters of a PDE model are unknown and need to be estimated from measurements. The numerical solution of such problems is very challenging as the underlying PDEs have to be solved repeatedly within outer optimization algorithms and the dimension of the parameters that need to be optimized might be very high or even infinite dimensional. PDE constrained optimization problems have been of interest for many decades. Classically, the underlying PDE (forward problem) is approximated by a high dimensional full order model (FOM) that results from discretization, e.g., by the Finite Element or Finite Volume method. Hence, the complexity of the optimization problem directly depends on the numbers of degrees of freedom (DOF) of the FOM. Mesh adaptivity has been advised to minimize the number of DOFs; see, e.g., [@BKR2000; @MR2556843; @Clever2012; @MR2905006; @MR1887737; @MR2892950] and the references therein.

**Model order reduction for PDE constrained optimization and optimal control.** A more recent approach is the usage of model order reduction (MOR) methods in order to replace the FOM by a surrogate reduced order model (ROM) of possibly very low dimension. MOR is a very active research field that has seen tremendous development in recent years, both from a theoretical and application point of view. For an introduction and overview we refer to the monographs and collections [@BCOW2017; @BOP+2017; @HRS2016; @QMN2016]. A particular promising model reduction approach for parameterized partial differential equations (pPDEs) is the Reduced Basis (RB) method that relies on the approximation of the solution manifold of pPDEs by low dimensional linear approximation spaces that are spanned from suitably selected particular solutions, called snapshots. A posteriori error estimation for solutions of the ROM with respect to the FOM is the basis for efficient Greedy algorithms to select the snapshots in a quasi-optimal way [@BCD+2011; @Haasdonk2013]. Alternatively, construction of reduced bases using proper orthogonal decomposition (POD) may be used [@GV17]. The construction of a reduced basis and the respective projected ROM is generally called the offline phase, whereas evaluating the ROM is called online phase.

There exists a large amount of literature using such reduced order surrogate models for optimization methods. A posteriori error estimates for reduced approximation of linear-quadratic optimization problems and parametrized optimal control problems with control constraints were studied, e.g., in [@Dede2012; @GK2011; @KTV13; @NRMQ2013; @OP2007]. In [@dihl15] an RB approach is proposed which also enables an estimation on the actual error on the control variable and not only on the gradient of the output functional. Certified Reduced Basis methods for parametrized elliptic optimal control problems with distributed controls were studied in [@KTGV2018]. With the help of an a posteriori error estimator, ROMs can be constructed with respect to a desired accuracy but also with respect to a local area in the parameter set [@EPR2010; @HDO2011]. For very high dimensional parameter sets, simultaneous parameter and state reduction has been advised [@HO2014; @HO2015; @LWG2010]. However, constructing a reduced order surrogate for a prohibitively expensive forward problem can also take a significant amount of computational resources. To remedy this, it is beneficial to use optimization methods that optimize on a local level of the control variable, assuming the surrogate only to be accurate enough in the respective parameter region. Hence, we require an approach which goes beyond the classical offline/online decomposition. Recently, RB methods have been advised with a progressive construction of ROMs [@BMV2020; @GHH2016; @ZF2015]. Also localized RB methods that are based on efficient localized a posteriori error control and online enrichment [@BEOR2017; @OS2015] overcome traditional offline/online splitting and are thus particularly well suited for applications in optimization or inverse problems [@OSS2018; @OS2017].

**Trust-Region reduced order models for second-order methods.** Trust-Region (TR) approaches are a class of optimization methods that are advantageous for the usage of locally accurate surrogate models. The main idea is to solve optimization sub-problems only in a local area of the parameter set which resolves the burden of constructing a global RB space. The problem that might occur is the fact that during this minimization one usually moves away from the original parameters on which the reduced order model was built, and the quality of the reduced model cannot be guaranteed anymore. For that reason, a priori and a posteriori error analysis are required to ensure accurate reduced order approximations for the optimization problem; cf. [@GV17; @HV08; @KTV13]. In [@AFS00; @SS13] a TR approach was proposed to control the quality of the (POD) reduced order model, referred to as TR-POD, a meanwhile well-established method in applications; cf. [@BC08; @CAN12].

TR methods ensure global convergence for locally convergent methods. In each iteration of the TR algorithm the nonlinear objective is replaced by a model function which can be optimized with much less effort; cf. [@CGT00; @NW06]. One suitable choice for the model is a reduced order discretization of the objective (e.g., by utilizing a second-order Taylor approximation). To ensure convergence to stationary points the accuracy of the model function and of its gradient have to be monitored. In [@RoggTV17] a posteriori error bounds are utilized to monitor the approximation quality of the gradient. We also refer to [@GGMRV17], where the authors utilize basis update strategies to improve the reduced order approximation scheme with respect to the optimization goal. The TR strategy can be combined with second-order methods for nonlinear optimization: with the Newton method to solve the reduced problem and with the SQP method for the all-at-once approach; cf. [@HV06].

Constraints on the control and the metric for the Trust-Region radius can affect the convergence of the method. For an error-aware TR method, the TR radius is directly characterized by the a posteriori error estimator for the cost functional of the surrogate model. Thus, the offline phase of the RB method can completely be omitted since the RB model can be adaptively enriched during the outer optimization loop. With this procedure the surrogate model eventually will have a high accuracy around the optimum of the optimization problem, ignoring the accuracy of the part which the outer (and inner) optimization loop does not approach at all. Error aware TR-RB methods can be utilized in many different ways. One possible TR-RB approach has been extensively studied in [@QGVW2017] for linear parametric elliptic equations, which ensures convergence of the nonlocal TR-RB. Note that the experiments in [@QGVW2017] are for up to six dimensional parameter sets without inequality constraints. In [@YM2013], the TR framework is combined with an efficient RB error bound for defining the Trust-Region in the design optimization of vibrating structures using frequency domain formulations. **Main results.** In this contribution we present several significant advances for adaptive Trust-Region Reduced Basis optimization methods for parameterized partial differential equations:

-   For the model function in the TR-RB approach, we follow a non-conforming dual (NCD) approach by choosing as model function the Lagrangian associated to the optimization problem. This permits more accurate results in terms of approximation of the optimal solution;

-   we provide efficiently computable a posteriori error estimates for all reduced quantities for different choices of the cost functional and its (approximate) gradient;

-   we rigorously prove the convergence of the TR-RB method with bilateral inequality constraints on the parameters;

-   we devise several new adaptive enrichment strategies for the progressive construction of the Reduced Basis spaces;

-   we demonstrate in numerical experiments that our new TR-RB methods outperform existing model reduction approaches for large scale optimization problems in well defined benchmark problems.

**Organization of the article.** In Section \[sec:problem\] we introduce the PDE constrained optimization problem and state first- and second-order optimality conditions. These serve as a basis for the full order discretization derived in Section \[sec:mor\]. Moreover, in Section \[sec:mor\] we introduce different strategies of model order reduction for the full order model and derive rigorous a posteriori error estimates for all equations, functionals, and gradient information. Section \[sec:TRRB\_and\_adaptiveenrichment\] is devoted to the derivation of Trust-Region – Reduced Basis methods and the presentation of the convergence analysis of the adaptive TR-RB algorithm. In addition, we discuss in detail several variants of new TR-RB algorithms that differ in their respective reduced gradient information as well as in the enrichment strategies for the construction of the corresponding reduced models. All variants are thoroughly analyzed numerically in Section \[sec:num\_experiments\], where we consider three well defined benchmark problems. We also compare with selected state of the art optimization methods from the literature.

Problem formulation {#sec:problem}
===================

Given $\mu_\mathsf{a},\mu_\mathsf{b}\in\mathbb{R}^P$ with $P \in {\mathbb{N}}$ we consider the compact and convex admissible parameter set $${\mathcal{P}}:= \left\{\mu\in\mathbb{R}^P\,|\,\mu_\mathsf{a} \leq \mu \leq \mu_\mathsf{b} \right\} \subset {\mathbb{R}}^P,$$ where $\leq$ is understood component-wise. Let $V$ be a real-valued Hilbert space with inner product $(\cdot \,,\cdot)$ and induced norm $\|\cdot\|$. We are interested in efficiently approximating PDE-constrained parameter optimization problems with the quadratic continuous cost functional $${\mathcal{J}}: V \times {\mathcal{P}}\to {\mathbb{R}}, \quad  (u,\mu) \mapsto {\mathcal{J}}(u, \mu) = \Theta(\mu) + j_\mu(u) + k_\mu(u, u),$$ where $\Theta: {\mathcal{P}}\to \mathbb{R}$ denotes a parameter functional and, for each $\mu \in {\mathcal{P}}$, $j_\mu \in V'$ is a parameter-dependent continuous linear functional and $k_\mu: V \times V \to {\mathbb{R}}$ a continuous symmetric bilinear form. To be more precise, we consider the following constrained minimization problem: [ $$\tag{P}
    \label{P}$$ ]{}

$$\min_{(u,\mu)\in V\times{\mathcal{P}}} {\mathcal{J}}(u, \mu), 
    \tag{P.a}\label{P.argmin}$$

subject to $(u,\mu)$ satisfying the *state – or primal – equation* $$\begin{aligned}
    a_\mu(u, v) = l_\mu(v) && \text{for all } v \in V,
    \tag{P.b}\label{P.state}
    \end{aligned}$$

where, for each $\mu \in {\mathcal{P}}$, $a_\mu: V \times V \to {\mathbb{R}}$ denotes a continuous and coercive symmetric bilinear form and $l_\mu \in V'$ denotes a continuous linear functional. For given $u \in V$, $\mu \in {\mathcal{P}}$, we introduce the primal residual $r_\mu^{\textnormal{pr}}(u) \in V'$ associated with by $$\begin{aligned}
r_\mu^{\textnormal{pr}}(u)[v] := l_\mu(v) - a_\mu(u, v) &&\text{for all }v \in V.
\label{eq:primal_residual}\end{aligned}$$ The primal residual plays a crucial role for a posteriori error analysis and for sensitivities of solution maps.

The Lagrange functional for is given by $\mathcal{L}(u,\mu,p) = {\mathcal{J}}(u,\mu) + r_\mu^{\textnormal{pr}}(u)[p]$ for $(u,\mu)\in V\times{\mathcal{P}}$ and for $p\in V$.

To apply RB methods efficiently, we require the parametrization of the problem to be separable from $V$ throughout the work. This separability is a standard assumption for RB methods and can be circumvented by using empirical interpolation techniques [@BMNP2004; @CS2010; @DHO2012].

\[asmpt:parameter\_separable\] We assume $a_\mu$, $l_\mu$, $j_\mu$, $k_\mu$ to be parameter separable with $\Xi^a, \Xi^l, \Xi^j, \Xi^k \in {\mathbb{N}}$ non-parametric components $a_\xi: V \times V \to {\mathbb{R}}$ for $1 \leq \xi \leq \Xi^a$, $l_\xi \in V'$ for $1 \leq \xi \leq \Xi^l$, $j_\xi \in V'$ for $1 \leq \xi \leq \Xi^j$ and $k_\xi: V \times V \to {\mathbb{R}}$ for $1 \leq \xi \leq \Xi^k$, and respective parameter functionals $\theta_\xi^a, \theta_\xi^l, \theta_\xi^j, \theta_\xi^k: {\mathcal{P}}\to \mathbb{R}$, such that $$\begin{aligned}
    a_\mu(u, v) &= \sum_{\xi = 1}^{\Xi^a} \theta_\xi^a(\mu)\, a_\xi(u, v)
    , &&&&& l_\mu(v) &= \sum_{\xi = 1}^{\Xi^l} \theta_\xi^l(\mu)\, l_\xi(v),
    \end{aligned}$$ and analogously for $j_\mu$ and $k_\mu$.

Due to Assumption \[asmpt:parameter\_separable\], all quantities which linearly depend on $a_\mu$, $l_\mu$, $j_\mu$ and $k_\mu$ (such as ${\mathcal{J}}$ or the primal residual) are also separable w.r.t. the parameter. Since we will use a Lagrangian ansatz for an explicit computation of derivatives, we require some notation that we use throughout this paper.

A note on differentiability {#sec:differentiability}
---------------------------

If ${\mathcal{J}}: V \times {\mathcal{P}}\to {\mathbb{R}}$ is Fréchet differentiable w.r.t. each $\mu \in {\mathcal{P}}$, for each $u \in V$ and each $\mu \in {\mathcal{P}}$ there exists a bounded linear functional $\partial_\mu {\mathcal{J}}(u, \mu) \in \mathbb{R}^P$, such that the Fréchet derivative of ${\mathcal{J}}$ w.r.t. its second argument in the direction of $\nu \in \mathbb{R}^P$ is given by $\partial_\mu {\mathcal{J}}(u, \mu) \cdot \nu$ (noting that the dual space of $\mathbb{R}^P$ is itself). We refer to $\partial_\mu {\mathcal{J}}(u, \mu)$ as the derivative w.r.t. $\mu$. In addition, for $u \in V$, $\mu \in {\mathcal{P}}$ we denote the *partial derivative* of ${\mathcal{J}}(u, \mu)$ w.r.t. the $i$-th component of $\mu$ by $d_{\mu_i}{\mathcal{J}}(u, \mu)$ for $1 \leq i \leq P$. Note that $d_{\mu_i}{\mathcal{J}}(u, \mu) = \partial_\mu {\mathcal{J}}(u, \mu)\cdot e_i $, where $e_i \in {\mathbb{R}}^P$ denotes the $i$-th canonical unit vector. Furthermore, we denote the *gradient* w.r.t. its second argument – the vector of components $d_{\mu_i}{\mathcal{J}}(u, \mu)$ – by the operator $\nabla_\mu {\mathcal{J}}: V \times {\mathcal{P}}\to {\mathbb{R}}^P$. Similarly, if ${\mathcal{J}}$ is Fréchet differentiable w.r.t. each $u \in V$, for each $u \in V$ and each $\mu \in {\mathcal{P}}$ there exists a bounded linear functional $\partial_u {\mathcal{J}}(u, \mu)\in V'$, such that the Fréchet derivative of ${\mathcal{J}}$ w.r.t. its first argument in any direction $v \in V$ is given by $\partial_u {\mathcal{J}}(u, \mu)[v]$. We refer to $\partial_u {\mathcal{J}}(u, \mu)$ simply as the derivative w.r.t. $u$. If ${\mathcal{J}}$ is twice Fréchet differentiable w.r.t. each $\mu \in {\mathcal{P}}$, we denote its *hessian* w.r.t. its second argument by the operator ${\mathcal{H}}_\mu {\mathcal{J}}: V \times {\mathcal{P}}\to {\mathbb{R}}^{P \times P}$.

We treat $a$, $l$, $j$ and $k$ in a similar manner, although, for notational compactness, we indicate their parameter-dependency differently. For instance, interpreting the bilinear form $a$ as a map $a: V \times V \times {\mathcal{P}}\to {\mathbb{R}}$, $\left(u, v, \mu\right) \mapsto a_\mu(u, v)$, we denote the Fréchet derivatives of $a$ w.r.t. the first, second and third argument of said map in the direction of $w\in V$, $\nu \in \mathbb{R}^P$ by $\partial_u a_\mu(u, v)[w] \in {\mathbb{R}}$, $\partial_v a_\mu(u, v)[w] \in {\mathbb{R}}$ and $\partial_\mu a_\mu(u, v)\cdot \nu \in {\mathbb{R}}$, respectively. Similarly, interpreting the linear functional $l$ as a map $l: V \times {\mathcal{P}}\to {\mathbb{R}}$, $(v, \mu) \mapsto l_\mu(v)$, we denote the Fréchet derivatives of $l$ w.r.t. the first and second argument of said map in the direction of $w \in V$, $\nu \in \mathbb{R}^P$ by $\partial_v l_\mu(v)[w] \in {\mathbb{R}}$ and $\partial_\mu l_\mu(v) \cdot \nu \in {\mathbb{R}}$, respectively. We omit the word Fréchet when referring to the derivatives of ${\mathcal{J}}$, $a$, $l$, $j$ and $k$, in order to simplify the notation, unless it is strictly necessary to specify it.

We apply this notation for Fréchet and partial derivatives for functionals and bilinear forms throughout this manuscript. Note that we denote the derivatives w.r.t. the symbol of the argument in the original definition of the functional or bilinear form, not w.r.t. the symbol of the actual argument, i.e. we use $\partial_u {\mathcal{J}}(u_\mu, \mu)$ for the derivative w.r.t. the first argument, not $\partial_{u_\mu} {\mathcal{J}}(u_\mu, \mu)$ or $\partial_v a_\mu(u, p)$ for the derivative w.r.t. the second argument, not $\partial_p a_\mu(u, p)$. Note also that, due to Assumption \[asmpt:parameter\_separable\], we can exchange the order of differentiation w.r.t. $V$ and $\mathbb{R}^P$, i.e. $\partial_u\big(\partial_\mu {\mathcal{J}}(u, \mu)\cdot\nu\big)[w] = \partial_\mu\big(\partial_u {\mathcal{J}}(u, \mu)[w]\big)\cdot\nu$.

\[asmpt:differentiability\] We assume $a_{\mu}$, $l_{\mu}$ and ${\mathcal{J}}$ to be twice Fréchet differentiable w.r.t. $\mu$. This obviously requires that all parameter-dependent coefficient functions in Assumption [[\[asmpt:parameter\_separable\]]{}]{} are twice differentiable as well.

\[prop:gateaux\_wrt\_V\] Due to the (bi-)linearity of $a$, $l$, $j$ and $k$, we can immediately compute their derivatives w.r.t. arguments in $V$. For $u, v \in V$, $\mu \in {\mathcal{P}}$, the derivatives of $a$, $l$ and ${\mathcal{J}}$ w.r.t. arguments in $V$ in the direction of $w \in V$ are given, respectively, by $$\partial_u a_\mu(u, v)[w] = a_\mu(w, v),\hspace{1em}\partial_v a_\mu(u, v)[w] = a_\mu(u, w), \hspace{1em} \partial_v l_\mu(v)[w] = l_\mu(w), \hspace{1em}\partial_u {\mathcal{J}}(u, \mu)[w] = j_\mu(w) + 2 k_\mu(w, u).$$

We compute the partial derivatives of $a$ and $l$ w.r.t. the parameter by means of their separable decomposition.

\[rem:gateaux\_wrt\_P\] For $\mu \in {\mathcal{P}}$, $u, v, \in V$ the derivatives of $a$ and $l$ w.r.t. $\mu$ in the direction of $\nu \in \mathbb{R}^P$ are given by $$\begin{aligned}
    \partial_\mu a_\mu(u, v) \cdot \nu &= \sum_{\xi = 1}^{\Xi^a} \big(\partial_\mu \theta_\xi^a(\mu) \cdot \nu\big)\, a_\xi(u, v)&&\text{and}&&&
    \partial_\mu l_\mu(v) \cdot \nu &= \sum_{\xi = 1}^{\Xi^l} \big(\partial_\mu \theta_\xi^l(\mu) \cdot \nu\big)\, l_\xi(v),
    \notag
    \end{aligned}$$ respectively, if $u, v$ do not depend on $\mu$. We also introduce the following shorthand notation for the derivative of functionals and bilinear forms w.r.t. the parameter in the direction of $\nu \in \mathbb{R}^P$, e.g. for $\mu \in {\mathcal{P}}$ we introduce $$\begin{aligned}
    \partial_\mu l_\mu &\cdot \nu \in V'&&& v \mapsto \big(\partial_\mu l_\mu \cdot \nu\big)(v) &:= \partial_\mu l_\mu(v)\cdot \nu&&\text{and}
    \notag\\
    \partial_\mu a_\mu &\cdot \nu \in V \times V \to {\mathbb{R}}&&& u, v \mapsto \big(\partial_\mu a_\mu \cdot \nu\big)(u, v) &:= \partial_\mu a_\mu(u, v)\cdot \nu,
    \notag
    \end{aligned}$$ and note that $\partial_\mu l_\mu$ and $\partial_\mu a_\mu$ are continuous and separable w.r.t. the parameter, owing to Assumption [\[asmpt:parameter\_separable\]]{}.

The bilinear form $a_\mu(\cdot\,,\cdot)$ is continuous and coercive for all $\mu\in{\mathcal{P}}$. Thus we can define the bounded solution map $\mathcal{S}:{\mathcal{P}}\to V$, $\mu \mapsto u_\mu:= \mathcal S(\mu)$, where $u_\mu$ is the unique solution to for a given $\mu$. The Fréchet derivatives of $\mathcal{S}$ are a common tool for RB methods and optimization, e.g., for constructing Taylor RB spaces that consist of the primal solution as well as their sensitivities (see [@HAA2017]) or for deriving optimality conditions for (see [@HPUU2009]).

\[prop:solution\_dmu\_eta\] Considering the solution map $\mathcal S:{\mathcal{P}}\to V$, $\mu \mapsto u_\mu$, its Fréchet derivative $d_{\nu} u_\mu \in V$ w.r.t. a direction $\nu\in\mathbb{R}^P$ is the unique solution of $$\begin{aligned}
\label{eq:primal_sens}
    a_\mu(d_{\nu}u_\mu, v) = \partial_\mu r_\mu^{\textnormal{pr}}(u_\mu)[v] \cdot \nu &&\text{for all } v \in V.
    \end{aligned}$$

We refer to [@HPUU2009] for the proof of this result.

Optimal solution and optimality conditions {#sec:first_order_optimality_conditions}
------------------------------------------

In this section, we discuss the existence of an optimal solution for problem . Then, we characterize a locally optimal solution through first- and second-order optimality conditions. Throughout the paper, a bar indicates optimality.

\[Thm:existence\_optimalsolution\] Problem admits an optimal solution pair $(\bar u, \bar \mu)\in V\times{\mathcal{P}}$, where $\bar u:= u_{\bar \mu}$ is the solution of for the parameter $\bar\mu$.

Note that the quantities involved in problem satisfies Assumption 1.44 in [@HPUU2009]. Thus the existence follows from [@HPUU2009 Theorem 1.45].

Let us introduce the reduced cost functional ${\hat{\mathcal{J}}}: {\mathcal{P}}\mapsto \mathbb{R},\,\mu\mapsto {\hat{\mathcal{J}}}(\mu) := {\mathcal{J}}(u_\mu, \mu)= {\mathcal{J}}( \mathcal S(\mu),\mu)$. Then problem is equivalent to the so-called reduced problem $$\begin{aligned}
\min_{\mu \in {\mathcal{P}}} {\hat{\mathcal{J}}}(\mu).
\tag{$\hat{\textnormal{P}}$}\label{Phat}\end{aligned}$$

1.  Since $r^{\textnormal{pr}}_\mu(u_\mu)[p] = 0$ for any $p\in V$, it follows that ${\hat{\mathcal{J}}}(\mu) = \mathcal L(u_\mu,\mu,p)$ for any $p\in V$.

2.  The cost functional ${\hat{\mathcal{J}}}$ is in general non-convex, thus the existence of a unique minimum for ${\hat{\mathcal{J}}}$ (and thus of ${\mathcal{J}}$) can not be guaranteed.

3.  Let $(\bar u, \bar \mu) \in V \times {\mathcal{P}}$ be a local optimal solution to with $\bar u := u_{\bar \mu}$ the solution of the primal equation for the parameter $\bar\mu$. Then the following constraint qualification holds true: For any $f\in V'$ there exists a pair $(u,\mu) \in V\times {\mathbb{R}}^P$ solving $$\begin{aligned}
                    a_{\bar \mu}(u, v) - \partial_\mu r_{\bar \mu}^{\textnormal{pr}}(\bar u)[v] \cdot \mu = f(v) &&\text{for all } v \in V.
            \end{aligned}$$

4.  Theorem [[\[Thm:existence\_optimalsolution\]]{}]{} does not provide any solution method.

One can derive first-order necessary optimality conditions in order to compute candidates for a local optimal solution of . We refer to [@HPUU2009 Cor. 1.3] for a proof of the following result:

\[prop:first\_order\_opt\_cond\] Let $(\bar u, \bar \mu) \in V \times {\mathcal{P}}$ be a local optimal solution to . Moreover, let Assumption [[\[asmpt:differentiability\]]{}]{} hold true. Then there exists a unique Lagrange multiplier $\bar p\in V$ such that the following first-order necessary optimality conditions hold:

\[eq:optimality\_conditions\] $$\begin{aligned}
        r_{\bar \mu}^{\textnormal{pr}}(\bar u)[v] &= 0 &&\text{for all } v \in V,
        \label{eq:optimality_conditions:u}\\
        \partial_u {\mathcal{J}}(\bar u,\bar \mu)[v] - a_{\bar\mu}(v,\bar p) &= 0 &&\text{for all } v \in V,
        \label{eq:optimality_conditions:p}\\
        (\nabla_\mu {\mathcal{J}}(\bar u,\bar \mu)+\nabla_{\mu} r^{\textnormal{pr}}_{\bar\mu}(\bar u)[\bar p]) \cdot (\nu-\bar \mu) &\geq 0 &&\text{for all } \nu \in {\mathcal{P}}. 
        \label{eq:optimality_conditions:mu}
        \end{aligned}$$

Note that resembles the state equation . From we deduce the *adjoint – or dual – equation* with unique solution $p_{\mu} \in V$ for a fixed $\mu \in {\mathcal{P}}$, i.e. $$\begin{aligned}
a_\mu(v, p_\mu) = \partial_u {\mathcal{J}}(u_\mu, \mu)[v]
= j_\mu(v) + 2 k_\mu(v, u_\mu)&&\text{for all } v \in V,
\label{eq:dual_solution}\end{aligned}$$ given the solution $u_\mu \in V$ to the state equation . From we observe that the variable $\bar p$ of the optimal triple solves the dual equation for $\bar \mu$. Similarly to the primal solution, we can consider the dual solution map $\mathcal A:{\mathcal{P}}\to V$, $\mu \mapsto \mathcal A(\mu) := p_\mu$, where $p_\mu$ is the solution of for the parameter $\mu$. In particular, $\bar p = p_{\bar \mu}$. For given $u, p \in V$, we also introduce the dual residual $r_\mu^{\textnormal{du}}(u, p) \in V'$ associated with by $$\begin{aligned}
r_\mu^{\textnormal{du}}(u, p)[v] := j_\mu(v) + 2k_\mu(v, u) - a_\mu(v, p)&&\text{for all }v \in V.
\label{eq:dual_residual}\end{aligned}$$ In addition, from the dual equation , we obtain the following formulation for the dual sensitivities.

\[prop:dual\_solution\_dmu\_eta\] Considering the dual solution map $\mathcal A:{\mathcal{P}}\to V$, $\mu \mapsto p_\mu$, its directional derivative $d_{\eta} p_\mu\in V$ w.r.t. a direction $\eta \in {\mathcal{P}}$ is the solution of $$\label{eq:dual_sens}
    \begin{split}
    a_\mu(q, d_{\eta} p_\mu) &= -\partial_\mu a_\mu(q, p_\mu)\cdot \eta + d_\mu \partial_u {\mathcal{J}}(u_\mu, \mu)[q] \cdot \eta        = \partial_\mu r_\mu^{\textnormal{du}}(u_\mu, p_\mu)[q] \cdot \eta + 2 k_\mu(q, d_{\eta}u_\mu)
    \end{split}$$ for all $q \in V$, where the latter equality holds for quadratic ${\mathcal{J}}$ as in .

Note that $\mathcal A$ is well defined because the bilinear form $a_{\mu}(\cdot \, , \cdot)$ is continuous and coercive. For a proof of the other claims we refer to [@HPUU2009], for instance.

Furthermore, we can compute first-order derivatives of ${\hat{\mathcal{J}}}$.

\[prop:grad\_Jhat\] For given $\mu\in{\mathcal{P}}$, the gradient of ${\hat{\mathcal{J}}}$, $\nabla_\mu{\hat{\mathcal{J}}}: {\mathcal{P}}\to {\mathbb{R}}^P$, is given by $$\begin{aligned}
    \nabla_{\mu}{\hat{\mathcal{J}}}(\mu) &= \nabla_{\mu}{\mathcal{J}}(u_{\mu}, \mu)+\nabla_{\mu}r_\mu^{\textnormal{pr}}(u_{\mu})[p_{\mu}]= \nabla_\mu \Theta(\mu) + \nabla_\mu j_\mu(u_\mu)
    +  \nabla_\mu k_\mu(u_\mu, u_\mu) + \nabla_\mu l_\mu(p_\mu) - \nabla_\mu a_\mu(u_\mu, p_\mu).
    \end{aligned}$$

This follows from , , and , cf. [@HPUU2009].

The proof of Proposition [[\[prop:grad\_Jhat\]]{}]{} relies on the fact that both $u_\mu$ and $p_\mu$ belong to the same space $V$; cf. *[@HPUU2009]*. In particular, for any $\mu\in{\mathcal{P}}$, we have $\nabla_\mu {\hat{\mathcal{J}}}(\mu) = \nabla_\mu \mathcal L(u_\mu,\mu,p_\mu)$.

For $\bar\mu$ satisfying the first-order necessary optimality conditions, we have that $\bar \mu$ is a stationary point of the cost functional ${\hat{\mathcal{J}}}$. Thus, $\bar \mu$ can be either a local minimum, a saddle point or a local maximum of the cost functional ${\hat{\mathcal{J}}}$ (and obviously the same relationship occurs between $(\bar u,\bar \mu)$ and ${\mathcal{J}}$). We thus consider second-order sufficient optimality conditions in order to characterize local minima of the functional ${\hat{\mathcal{J}}}$, requiring its hessian.

\[prop:hessian\_Jhat\] The hessian of ${\hat{\mathcal{J}}}$, ${\hat{\mathcal{H}}}_\mu := {\mathcal{H}}_\mu {\hat{\mathcal{J}}}: {\mathcal{P}}\to {\mathbb{R}}^{P \times P}$, is determined by its application to a direction $\nu\in \mathbb{R}^P$, given by $$\begin{aligned}
    {\hat{\mathcal{H}}}_\mu(\mu) \cdot \nu
    = \nabla_\mu\Big(&\partial_u {\mathcal{J}}(u_\mu, \mu)[d_{\nu}u_\mu] + l_\mu(d_{\nu}p_\mu) - a_\mu(d_{\nu}u_\mu, p_\mu) - a_\mu(u_\mu, d_{\nu}p_\mu)
    \notag\\
    &+\big(\partial_\mu {\mathcal{J}}(u_\mu, \mu)+ \partial_\mu l_\mu(p_\mu)- \partial_\mu a_\mu(u_\mu, p_\mu)\big)\cdot \nu\Big),
    \notag
    \end{aligned}$$ where $u_\mu, p_\mu \in V$ denote the primal and dual solutions, respectively. For a quadratic ${\mathcal{J}}$ as in the above formula simplifies to $$\begin{aligned}
    {\hat{\mathcal{H}}}_\mu(\mu) \cdot \nu
    = \nabla_\mu\Big(&j_\mu(d_{\nu}u_\mu) + 2k_\mu(d_{\nu}u_\mu, u_\mu) + l_\mu(d_{\nu}p_\mu) - a_\mu(d_{\nu}u_\mu, p_\mu) - a_\mu(u_\mu, d_{\nu}p_\mu)
    \notag\\
    &+\big(\partial_\mu {\mathcal{J}}(u_\mu, \mu)+ \partial_\mu l_\mu(p_\mu)- \partial_\mu a_\mu(u_\mu, p_\mu)\big)\cdot \nu\Big).
    \notag
    \end{aligned}$$

See, e.g., [@HPUU2009] for the first part. The second one follows from Remark \[prop:gateaux\_wrt\_V\].

\[prop:second\_order\] Let Assumption [[\[asmpt:differentiability\]]{}]{} hold true. Suppose that $\bar \mu\in {\mathcal{P}}$ satisfies the first-order necessary optimality conditions . If ${\hat{\mathcal{H}}}_\mu(\bar \mu)$ is positive definite on the *critical cone* $\mathcal C(\bar\mu)$ at $\bar\mu\in{\mathcal{P}}$, i.e., if $\nu \cdot({\hat{\mathcal{H}}}_\mu(\bar \mu)\cdot \nu) > 0$ for all $\nu\in\mathcal C(\bar\mu)\setminus\{0\}$, with $$\begin{aligned}
    \mathcal C(\bar\mu):= \big\{\nu\in\mathbb{R}^P\,\big|\, \exists \mu\in{\mathcal{P}},\,c_1>0: \nu = c_1(\mu-\bar \mu),\, \nabla_\mu {\hat{\mathcal{J}}}(\bar\mu)\cdot \nu = 0  \big\},
    \end{aligned}$$ then $\bar \mu$ is a strict local minimum of .

For this result we refer to [@CasTr15; @Nocedal06], for instance.

High dimensional discretization and model order reduction {#sec:mor}
=========================================================

We first discretize the optimization problem as well as the corresponding optimality conditions using the classical Ritz-Galerkin projection into a possibly high dimensional approximation space $V_h \subset V$, such as conforming Finite Elements. Note that we restrict ourselves to a conforming approximation for simplicity and that we do not further specify the choice of $V_h$, as neither impacts the analysis below. Based on this idea, we then derive different ways for the ROM using the Reduced Basis method with possibly different reduced primal and dual state spaces. Thus, the resulting ROM optimality system will in general not be equivalent to a Ritz-Galerkin projection of the FOM one onto a reduce space $V_{{r}}\subset V_h$. For this reason, we will introduce a non-conforming dual-corrected (NCD-corrected) approach; cf. Section \[sec:ncd\_approach\].

FOM for the optimality system {#sec:problem_fom}
-----------------------------

For the discretization of the optimization problem we assume that a finite-dimensional subspace $V_h \subset V$ is given and obtain the FOM for the optimality system of by Ritz-Galerkin projection of equations onto $V_h$. In particular, we have for each $\mu \in {\mathcal{P}}$ the solution $u_{h, \mu} \in V_h$ of the *discrete primal equation* $$\begin{aligned}
{a_{\mu}}(u_{h, \mu}, v_h) = {l_{\mu}}(v_h) &&\text{for all } v_h \in V_h,
\label{eq:state_h}\end{aligned}$$ and hence ${r_{\mu}}^{\textnormal{pr}}(u_{h, \mu})[v_h] = 0$ for all $v_h \in V_h$, $\mu \in {\mathcal{P}}$. We also have for each $\mu \in {\mathcal{P}}$ the solution $p_{h, \mu} \in V_h$ of the *discrete dual equation* $$\begin{aligned}
{a_{\mu}}(v_h, p_{h, \mu}) = \partial_u {\mathcal{J}}(u_{h, \mu}, \mu)[v_h] = {j_{\mu}}(v_h) + 2 {k_{\mu}}(v_h, u_{h, \mu}) &&\text{for all } v_h \in V_h,
\label{eq:dual_solution_h}\end{aligned}$$ and hence ${r_{\mu}}^{\textnormal{du}}(u_{h, \mu}, p_{h, \mu})[v_h] = 0$ for all $v_h \in V_h$, $\mu \in {\mathcal{P}}$. Similarly, the *discrete primal sensitivity equations* for solving for $d_{\nu} u_{h, \mu} \in V_h$ as well as *discrete dual sensitivity equations* for solving for $d_{\nu} p_{h, \mu} \in V_h$ at any direction $\nu \in \mathbb{R}^P$ follow analogously to Propositions \[prop:solution\_dmu\_eta\] and \[prop:dual\_solution\_dmu\_eta\]. Furthermore, ${\hat{\mathcal{J}}}$ is approximated by the *discrete reduced functional* $$\begin{aligned}
{\hat{\mathcal{J}}}_h(\mu) := {\mathcal{J}}(u_{h, \mu}, \mu) = \mathcal L(u_{h,\mu},\mu,p_h) && \text{for all } p_h\in V_h,
\label{eq:Jhat_h}\end{aligned}$$ where $u_{h, \mu} \in V_h$ is the solution of and we formulate the discrete optimization problem $$\begin{aligned}
\min_{\mu \in {\mathcal{P}}} {\hat{\mathcal{J}}}_h(\mu).
\tag{$\hat{\textnormal{P}}_h$}\label{Phat_h}\end{aligned}$$ Further, $\bar\mu_h$ denotes a locally optimal solution to satisfying the first- and second-order optimality conditions.

Since $u_{h,\mu}$ and $p_{h,\mu}$ belong to the same space $V_h$, Propositions [[\[prop:first\_order\_opt\_cond\]-[\[prop:grad\_Jhat\]]{},[\[prop:hessian\_Jhat\]]{}-\[prop:second\_order\]]{}]{} from Section [[\[sec:problem\_fom\]]{}]{} hold for the FOM as well, with all quantities replaced by their discrete counterparts.

As usual in the context of RB methods, we eliminate the issue of “truth” by assuming that the high dimensional space $V_h$ is accurate enough to approximate the true solution.

\[asmpt:truth\] We assume that the primal discretization error $\|u_\mu - u_{h, \mu}\|$, the dual error $\|p_\mu - p_{h, \mu}\|$, the primal sensitivity errors $\|d_{\mu_i} u_\mu - d_{\mu_i} u_{h, \mu}\|$ and the dual sensitivity errors $\|d_{\mu_i} p_\mu - d_{\mu_i} p_{h, \mu}\|$ are negligible for all $\mu \in {\mathcal{P}}$, $1 \leq i \leq P$.

To define suitable ROMs, in what follows, we assume that we have computed problem adapted RB spaces $V_{{r}}^{\textnormal{pr}}, V_{{r}}^{\textnormal{du}}\subset V_h$, the construction of which is detailed in Section \[sec:construct\_RB\]. We stress here that $V_{{r}}^{\textnormal{pr}}$ and $V_{{r}}^{\textnormal{du}}$ might not coincide, this will imply further discussions of the RB approximation of the optimality system .

ROM for the optimality system – Standard approach {#sec:standard_approach}
-------------------------------------------------

Given a RB space $V_{{r}}^{\textnormal{pr}}\subset V_h$ of low dimension $n := \dim V_{{r}}^{\textnormal{pr}}$ and dual RB space $V_{{r}}^{\textnormal{du}}\subset V_h$ of low dimension $m := \dim V_{{r}}^{\textnormal{du}}$, we obtain the RB approximation of state and adjoint equations as follows:

\[eq:optimality\_conditionsRB\]

-   RB approximation for : For each $\mu \in {\mathcal{P}}$ the primal variable $u_{{{r}}, \mu} \in V_{{r}}^{\textnormal{pr}}$ of the *RB approximate primal equation* is defined through $$\begin{aligned}
            {a_{\mu}}(u_{{{r}}, \mu}, v_{{r}}) = {l_{\mu}}(v_{{r}}) &\qquad \text{for all } v_{{r}}\in V_{{r}}^{\textnormal{pr}}.
            \label{eq:state_red}
            \end{aligned}$$

-   RB approximation for : For each $\mu \in {\mathcal{P}}$, $u_{{{r}}, \mu} \in V_{{r}}^{\textnormal{pr}}$ the dual/adjoint variable $p_{{{r}}, \mu} \in V_{{r}}^{\textnormal{du}}$ satisfies the *RB approximate dual equation* through $$\begin{aligned}
            {a_{\mu}}(q_{{r}}, p_{{{r}}, \mu}) = \partial_u {\mathcal{J}}(u_{{{r}}, \mu}, \mu)[q_{{r}}] = {j_{\mu}}(q_{{r}}) + 2 {k_{\mu}}(q_{{r}}, u_{{{r}}, \mu}) &&\text{for all } q_{{r}}\in V_{{r}}^{\textnormal{du}}.
            \label{eq:dual_solution_red}
            \end{aligned}$$

Analogously to Propostion \[prop:solution\_dmu\_eta\], we define the *RB solution map* $\mathcal{S}_{{r}}:{\mathcal{P}}\to V_{{r}}^{\textnormal{pr}}$ by $\mu \mapsto u_{{{r}}, \mu}$ and analogously to Propostion \[prop:dual\_solution\_dmu\_eta\] the *RB dual solution map* $\mathcal{A}_{{r}}:{\mathcal{P}}\to V_{{r}}^{\textnormal{du}}$ by $\mu \mapsto p_{{{r}}, \mu}$, where $u_{{{r}}, \mu}$ and $p_{{{r}}, \mu}$ denote the primal and dual reduced solutions of and , respectively. To approximate , we introduce the *RB reduced functional* by $$\begin{aligned}
{\hat{J}}_{{r}}(\mu) := {\mathcal{J}}(u_{{{r}}, \mu}, \mu)={\mathcal{J}}(\mathcal S_r(\mu), \mu),&&
\text{where $u_{{{r}}, \mu} \in V_{{r}}^{\textnormal{pr}}$ is the solution of \eqref{eq:state_red}}
\label{eq:Jhat_red}\end{aligned}$$ instead of ${\hat{\mathcal{J}}}_h$ and the problem of finding a locally optimal solution $\bar \mu_{{r}}$ of $$\begin{aligned}
\min_{\mu \in {\mathcal{P}}} {\hat{J}}_{{r}}(\mu).
\label{eq:Phat_red_uncorrected}\end{aligned}$$ Now, a solution to the optimality system is approximated by the RB triple $(u_{{{r}},\bar\mu_{{r}}},\bar\mu_{{r}},p_{{{r}},\bar\mu_{{r}}})$.

As proposed in [@QGVW2017], for computing an approximation of the gradient of ${\hat{J}}_{{r}}$, the gradient from Propostion \[prop:grad\_Jhat\] can be utilized by replacing $u_{\mu}$ and $p_{\mu}$ with their RB counterparts. However, it can not be guaranteed in general that the computed gradient is the actual gradient of ${\hat{J}}_{{r}}$, if $V^{\textnormal{pr}}_{{r}}$ and $V^{\textnormal{du}}_{{r}}$ are chosen to be different. To see this, we consider first the Lagrangian and note that, for $1\leq i\leq P$ and all $p\in V$, it holds $${\hat{J}}_{{r}}(\mu)  = \mathcal{L}(u_{{{r}},\mu},\mu,p), \qquad
    \label{Gradient_lagrangian}
     \big(\nabla_\mu {\hat{J}}_{{r}}(\mu)\big)_i  =  \partial_u \mathcal{L}(u_{{{r}},\mu},\mu,p)[d_{\mu_i}u_{{{r}},\mu}] 
     + d_{\mu_i} \mathcal{L}(u_{{{r}},\mu},\mu,p).$$ Now, following [@QGVW2017], we define the *inexact gradient* ${\widetilde\nabla}_\mu {\hat{J}}_{{r}}:{\mathcal{P}}\to {\mathbb{R}}^P$ by $$\label{naive:red_grad}
\big({\widetilde\nabla}_\mu {\hat{J}}_{{r}}(\mu)\big)_i := d_{\mu_i}{\mathcal{J}}(u_{{{r}}, \mu}, \mu) + d_{\mu_i} r_\mu^{\textnormal{pr}}(u_{{{r}}, \mu})[p_{{{r}}, \mu}] = d_{\mu_i} \mathcal{L}(u_{{{r}},\mu},\mu,p_{{{r}},\mu})$$ for all $1 \leq i \leq P$ and $\mu \in {\mathcal{P}}$, where $u_{{{r}}, \mu} \in V_{{r}}^{\textnormal{pr}}$ and $p_{{{r}}, \mu} \in V_{{r}}^{\textnormal{du}}$ denote the primal and approximate dual reduced solutions of and , respectively. With the superscript $\sim$ we stress that ${\widetilde\nabla}_\mu {\hat{J}}_{{r}}(\mu)$ is not the actual gradient of ${\hat{J}}_{{r}}$, but its approximation. Choosing $p= p_{{{r}},\mu}\in V^{\textnormal{du}}_{{r}}$ in and considering lead to $$\begin{aligned}
\big(\nabla_\mu {\hat{J}}_{{r}}(\mu)\big)_i & = \partial_u \mathcal{L}(u_{{{r}},\mu},\mu,p_{{{r}},\mu})[d_{\mu_i}u_{{{r}},\mu}] + {\widetilde\nabla}_\mu {\hat{J}}_{{r}}(\mu).\end{aligned}$$ Note that, in general, it does not hold that $\partial_u\mathcal L(u_{{{r}},\mu},\mu, p_{{{r}},\mu})= 0$, since is not the dual equation with respect to the optimization problem , cf. [@HPUU2009 Section 1.6.4], which would only be true if $V^{\textnormal{du}}_{{r}}\subseteq V^{\textnormal{pr}}_{{r}}$. Thus, defines only an approximation of the true gradient of ${\hat{J}}_{{r}}$ with the choice made in [@QGVW2017]. This introduces an additional approximation error in reconstructing the solution of the optimality system , which is well visible in our numerical experiments (see Section \[sec:mmexc\_example\]): the standard RB approach leads to a significant lack in accuracy, requiring additional steps to enrich the RB space and cover this gap. We therefore propose to add a correction term to ${\hat{J}}_{{r}}$ based on the previous remarks.

ROM for the optimality system – NCD-corrected approach {#sec:ncd_approach}
------------------------------------------------------

Following the primal-dual RB approach for linear output functionals [@HAA2017 Section 2.4], it is more suitable to add a correction term to the output functional for which improved error estimates are available. We seek to minimize the Lagrangian corresponding to problem . A similar approach, in the context of adaptive finite elements, can be found in [@BKR2000; @Rannacher2006]. We utilize to extend the primal-dual RB approach of [@HAA2017 Section 2.4] to quadratic output functionals and define the *NCD-corrected RB reduced functional* $$\begin{aligned}
{{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu) := \mathcal{L}(u_{{{r}},\mu},\mu,p_{{{r}},\mu}) = {\hat{J}}_{{r}}(\mu) + r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[p_{{{r}},\mu}]
\label{eq:Jhat_red_corected}\end{aligned}$$ with $u_{{{r}}, \mu} \in V_{{r}}^{\textnormal{pr}}$ and $p_{{{r}}, \mu} \in V_{{r}}^{\textnormal{du}}$ the solutions of and , respectively. Note that ${{{{\hat{\mathcal{J}}}_{{r}}}}}$ coincides with the functional ${\hat{J}}_{{r}}$ in if $V_{{r}}^{\textnormal{du}}= V_{{r}}^{\textnormal{pr}}$. We then consider the *RB reduced optimization problem* of finding a locally optimal solution $\bar \mu_{{r}}$ of $$\begin{aligned}
\min_{\mu \in {\mathcal{P}}} {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu).
\tag{$\hat{\textnormal{P}}_{{r}}$}\label{Phat_{{r}}}\end{aligned}$$

Computing the actual gradient of ${\hat{\mathcal{J}}}_{{r}}$ results in the next proposition, proved following [@HPUU2009 Section 1.6.2].

\[prop:true\_corrected\_reduced\_gradient\_adj\] The $i$-th component of the true gradient of ${{{{\hat{\mathcal{J}}}_{{r}}}}}$ is given by $$\begin{aligned}
    \big(\nabla_\mu {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)\big)_i & = d_{\mu_i}{\mathcal{J}}(u_{{{r}},\mu},\mu) + d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[p_{{{r}},\mu}+w_{{{r}},\mu}] - d_{\mu_i}r^{\textnormal{du}}_\mu (u_{{{r}},\mu},p_{{{r}},\mu})[z_{{{r}},\mu}] 
    \end{aligned}$$ where $u_{{{r}}, \mu} \in V_{{r}}^{\textnormal{pr}}$ and $p_{{{r}}, \mu} \in V_{{r}}^{\textnormal{du}}$ denote the RB approximate primal and dual solutions of and , $z_{{{r}},\mu} \in V_{{r}}^{\textnormal{du}}$ solves $$\label{A2:z_eq}
    a_\mu(z_{{{r}},\mu},q) = -r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[q] \quad \forall q\in V^{\textnormal{du}}_r$$ and $w_{{{r}},\mu} \in V_{{r}}^{\textnormal{pr}}$ solves $$\label{A2:w_eq}
    a_\mu(v,w_{{{r}},\mu}) = r_\mu^{\textnormal{du}}(u_{{{r}},\mu},p_{{{r}},\mu})[v]-2k_\mu(z_{{{r}},\mu},v), \quad \forall v\in V^{\textnormal{pr}}_r.$$

A posteriori error analysis {#sec:a_post_error_estimates}
---------------------------

A posteriori error estimates are required for controlling the accuracy of the reduced order model. In addition, we also use them for the error aware TR method (which is explained in Section \[sec:TR\]). We derive a posteriori error estimates for all reduced terms that we need for the TR method. Moreover, we suggest further advances for the reduction of sensitivities and gradients. From a model reduction perspective, these error estimates need to be computed efficiently such that the time for the evaluation for many parameters can be neglected. Note that Assumption \[asmpt:parameter\_separable\] is crucial for this, since it allows to precompute most of the required terms. For any functional $l \in V_h'$ or bilinear form $a: V_h \times V_h \to {\mathbb{R}}$, we denote their dual or operator norms $\|l\|$ and $\|a\|$ by the continuity constants ${\gamma_{l}}$ and ${\gamma_{a}}$, respectively. The same consideration applies for the norm $\|\cdot\|$ in $V_h'$ of the residuals. For $\mu \in {\mathcal{P}}$, we denote the coercivity constant of ${a_{\mu}}$ w.r.t. the $V_h$-norm by $\underline{{a_{\mu}}} > 0$.

### Standard RB estimates for the optimality system {#sec:estimation_for_the_standard_and_corrected_reduced_order_model}

We start with the residual based a posteriori error estimation for the primal variable, which is a standard result from RB theory and has extensively been used in the literature. For a proof, we refer to [@rozza2007].

\[prop:primal\_rom\_error\] For $\mu \in {\mathcal{P}}$ let $u_{h, \mu} \in V_h$ be the solution of and $u_{{{r}}, \mu} \in V_{{r}}^{\textnormal{pr}}$ the solution of . Then it holds $$\begin{aligned}
    \|u_{h, \mu} - u_{{{r}}, \mu}\| \leq \Delta_{\textnormal{pr}}(\mu)  := \underline{{a_{\mu}}}^{-1}\, \|{r_{\mu}}^{\textnormal{pr}}(u_{{{r}}, \mu})\|.
    \notag
    \end{aligned}$$

For the reduced dual problem, a similar idea can be used to derive the following estimation, accounting for the fact that $p_{{{r}}, \mu}$ is not a Galerkin projection of $p_{h, \mu}$. For a proof, we refer to [@QGVW2017 Lemma 3].

\[prop:dual\_rom\_error\] For $\mu \in {\mathcal{P}}$, let $p_{h, \mu} \in V_h$ be the solution of and $p_{{{r}}, \mu} \in V_{{r}}^{\textnormal{du}}$ the solution of . Then it holds $$\begin{aligned}
    \|p_{h, \mu} - p_{{{r}}, \mu}\| &\leq \Delta_{\textnormal{du}}(\mu) := \underline{{a_{\mu}}}^{-1}\big(2 {\gamma_{{k_{\mu}}}}\;\Delta_{\textnormal{pr}}(\mu) + \|{r_{\mu}}^{\textnormal{du}}(u_{{{r}}, \mu}, p_{{{r}}, \mu})\|\Big).
    \notag
    \end{aligned}$$

In the next proposition we state the result of the standard approach from [@QGVW2017 Theorem 4]. Furthermore we show an improved version by using, in contrast to [@QGVW2017], the NCD-corrected reduced functional, which results in an optimal higher order a posteriori upper bound without lower order terms.

\[prop:Jhat\_error\]

(i) With the notation from above, we have for the standard RB reduced cost functional $$\begin{aligned}
            |{\hat{\mathcal{J}}}_h(\mu) - &{\hat{J}}_r(\mu)| \leq \Delta_{{\hat{J}}_{{r}}}(\mu) :=  \Delta_{\textnormal{pr}}(\mu) \|{r_{\mu}}^{\textnormal{du}}(u_{{{r}}, \mu}, p_{{{r}},\mu})\| + \Delta_{\textnormal{pr}}(\mu)^2 {\gamma_{{k_{\mu}}}} + \big|r_{\mu}^{\textnormal{pr}}(u_{{{r}}, \mu})[p_{{{r}},\mu}]\big|.
            \notag
            \end{aligned}$$

(ii) Furthermore, we have for the NCD-corrected RB reduced cost functional (or equivalently for the Lagrangian for any $p\in V_h$) $$\begin{aligned}
            |{\hat{\mathcal{J}}}_h(\mu) - &{{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)| = |\mathcal{L}(u_{h,\mu},\mu,p)-\mathcal{L}(u_{{{r}},\mu},\mu,p)| \leq \Delta_{{{{{\hat{\mathcal{J}}}_{{r}}}}}}(\mu) :=  \Delta_{\textnormal{pr}}(\mu) \|{r_{\mu}}^{\textnormal{du}}(u_{{{r}}, \mu}, p_{{{r}},\mu})\| + \Delta_{\textnormal{pr}}(\mu)^2 {\gamma_{{k_{\mu}}}}.
            \notag
            \end{aligned}$$

We refer to [@QGVW2017 Theorem 4] for a proof of $(i)$. Regarding $(ii)$, using the shorthand $e_{h, \mu}^{\textnormal{pr}}:= u_{h, \mu} - u_{{{r}}, \mu}$ and $a_\mu(e_{h, \mu}^{\textnormal{pr}},p_{{{r}},\mu}) = r_{\mu}^{\textnormal{pr}}(u_{{{r}}, \mu})[p_{{{r}},\mu}]$ lead us to $$\begin{aligned}
        |{\hat{\mathcal{J}}}_h(\mu) - {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)| &= |j_{h, \mu}(e_{h, \mu}^{\textnormal{pr}}) + {k_{\mu}}(u_{h, \mu}, u_{h, \mu}) - {k_{\mu}}(u_{{{r}}, \mu}, u_{{{r}}, \mu}) - a_\mu(e_{h, \mu}^{\textnormal{pr}},p_{{{r}},\mu})| \notag\\
        &= |r_{\mu}^{\textnormal{du}}(u_{{{r}}, \mu}, p_{{{r}},\mu})[e_{h, \mu}^{\textnormal{pr}}] - 2{k_{\mu}}(u_{{{r}}, \mu},e_{h, \mu}^{\textnormal{pr}})
        + {k_{\mu}}(u_{h, \mu}, u_{h, \mu}) - {k_{\mu}}(u_{{{r}}, \mu}, u_{{{r}}, \mu})|  \notag \\
        &\leq \|r_{\mu}^{\textnormal{du}}(u_{{{r}}, \mu}, p_{{{r}},\mu})\|\; \|e_{h, \mu}^{\textnormal{pr}}\| +{\gamma_{{k_{\mu}}}} \; \|e_{h, \mu}^{\textnormal{pr}}\|^2,
        \notag
        \end{aligned}$$ where we used the definition of the dual residual in the second equality and Cauchy-Schwarz for the inequality. The assertion follows by using Propostion \[prop:primal\_rom\_error\].

\[continuity\_of\_estimator\] The estimator $\Delta_{{{{{\hat{\mathcal{J}}}_{{r}}}}}}(\mu)$ is continuous w.r.t. $\mu$, since the Riesz-representative of the residual is continuous.

For the inexact and NCD-corrected gradient, we derive the following a posteriori estimators.

\[prop:grad\_Jhat\_error\]

(i) For the inexact gradient ${\widetilde\nabla}_\mu {\hat{J}}_{{r}}(\mu)$ from the standard-RB approach , we have $$\begin{aligned}
            \big\|\nabla_\mu {\hat{\mathcal{J}}}_h(\mu) - {\widetilde\nabla}_\mu {\hat{J}}_{{r}}(\mu)\big\|_2 &\leq \Delta_{{\widetilde\nabla}{\hat{J}}_{{r}}}(\mu) = \big\|\underline{\Delta_{{\widetilde\nabla}{\hat{J}}_{{r}}}(\mu)}\big\|_2 \quad\quad\text{with}
            \notag\\
            \big(\underline{\Delta_{{\widetilde\nabla}{\hat{J}}_{{r}}}(\mu)}\big)_i  
            := 2\Delta_{\textnormal{pr}}(\mu) \|u_{{{r}}, \mu}\|\; {\gamma_{d_{\mu_i} {k_{\mu}}}}
            &+ \Delta_{\textnormal{pr}}(\mu)\big({\gamma_{d_{\mu_i} {j_{\mu}}}} + {\gamma_{d_{\mu_i} {a_{\mu}}}}\; \|p_{{{r}}, \mu}\|\big)
            \notag\\
          + \Delta_{\textnormal{du}}(\mu)\big({\gamma_{d_{\mu_i} {l_{\mu}}}} &+ {\gamma_{d_{\mu_i} {a_{\mu}}}}\; \|u_{{{r}}, \mu}\|\big)
            + \Delta_{\textnormal{pr}}(\mu)\; \Delta_{\textnormal{du}}(\mu)\; {\gamma_{d_{\mu_i} {a_{\mu}}}}
            + (\Delta_{\textnormal{pr}})^2(\mu)\; {\gamma_{d_{\mu_i} {k_{\mu}}}}.
            \notag
            \end{aligned}$$

(ii) For the gradient $\nabla_{\mu} {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)$ of the NCD-corrected reduced functional, computed with the adjoint approach from Definition [[\[prop:true\_corrected\_reduced\_gradient\_adj\]]{}]{}, we have $$\begin{aligned}
            &\big\|\nabla_\mu {\hat{\mathcal{J}}}_h(\mu) - \nabla_{\mu} {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)\big\|_2 \leq \Delta^*_{\nabla {\hat{\mathcal{J}}}_{{r}}}(\mu) 
        = \big\|\underline{\Delta^*_{\nabla {\hat{\mathcal{J}}}_{{r}}}(\mu)}\big\|_2 \quad\quad\text{with} \\
          \big(\underline{\Delta^*_{\nabla {\hat{\mathcal{J}}}_{{r}}}(\mu)}\big)_i &:=
            2\Delta_{\textnormal{pr}}(\mu) \|u_{{{r}}, \mu}\|\; {\gamma_{d_{\mu_i} {k_{\mu}}}}
            + \Delta_{\textnormal{pr}}(\mu)\big({\gamma_{d_{\mu_i} {j_{\mu}}}} + {\gamma_{d_{\mu_i} {a_{\mu}}}}\; \|p_{{{r}}, \mu}\|\big)\\
            &+ \Delta_{\textnormal{du}}(\mu)\big({\gamma_{d_{\mu_i} {l_{\mu}}}} + {\gamma_{d_{\mu_i} {a_{\mu}}}}\; \|u_{{{r}}, \mu}\|\big)
            + \Delta_{\textnormal{pr}}(\mu)\; \Delta_{\textnormal{du}}(\mu)\; {\gamma_{d_{\mu_i} {a_{\mu}}}}
            + (\Delta_{\textnormal{pr}})^2(\mu)\; {\gamma_{d_{\mu_i} {k_{\mu}}}} \\
            &+ ({\gamma_{d_{\mu_i}l_{\mu}}} + {\gamma_{d_{\mu_i}a_{\mu}}} \| u_{{{r}},\mu} \|) \underline{{a_{\mu}}}^{-1} \big( \|r_\mu^{\textnormal{du}}(u_{{{r}},\mu},p_{{{r}},\mu})\| 
            + 2{\gamma_{k_\mu}} \underline{{a_{\mu}}}^{-1} \|r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})\| \big) \\
            &+ \underline{{a_{\mu}}}^{-1} \|r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})\| \big( {\gamma_{d_{\mu_i}j}}+ 2{\gamma_{d_{\mu_i}k}} \|u_{{{r}},\mu}\|
            + {\gamma_{d_{\mu_i}a}} \|p_{{{r}},\mu}\| \big).
            \end{aligned}$$

\(i) For $\Delta_{{\widetilde\nabla}{\hat{J}}_{{r}}}(\mu)$, we have $$\begin{aligned}
        \big(\nabla_\mu {\hat{\mathcal{J}}}_h(\mu) - {\widetilde\nabla}_{\mu} {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)\big)_i
        &=d_{\mu_i}{\mathcal{J}}(u_{h, \mu}, \mu) - d_{\mu_i}{\mathcal{J}}(u_{{{r}}, \mu}, \mu)  + d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{h,\mu})[p_{h,\mu}] - d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[p_{{{r}},\mu}].
        \end{aligned}$$ Regarding the first contribution, we obtain with $\|u_{h, \mu}\|_h \leq \|e_{h, \mu}^{\textnormal{pr}}\|_h + \|u_{{{r}}, \mu}\|_h$ $$\begin{aligned}
        \big|d_{\mu_i}{\mathcal{J}}(u_{h, \mu}, \mu) - d_{\mu_i}{\mathcal{J}}(u_{{{r}}, \mu}, \mu)\big| &=
        | d_{\mu_i}j_{h, \mu}(e_{h, \mu}^{\textnormal{pr}}) +  d_{\mu_i}k_{\mu}(e_{h, \mu}^{\textnormal{pr}}, u_{{{r}}, \mu}) +  d_{\mu_i}k_{\mu}(u_{h, \mu}, e_{h, \mu}^{\textnormal{pr}})|
        \\
        &\leq \Delta_{\textnormal{pr}}(\mu)\Big( {\gamma_{d_{\mu_i} {j_{\mu}}}} + {\gamma_{d_{\mu_i} {k_{\mu}}}}\big(2 \|u_{{{r}}, \mu}\| + \Delta_{\textnormal{pr}}(\mu)\big)\Big).
        \end{aligned}$$ For the other contributions we refer to [@QGVW2017 Theorem 5].

\(ii) For the adjoint estimator $\Delta^*_{\nabla_\mu {\hat{\mathcal{J}}}_{{r}}}$, we have $$\begin{aligned}
        \partial_\mu {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)\cdot\nu & = \partial_\mu {\mathcal{J}}(u_{{{r}},\mu},\mu)\cdot\nu + \partial_\mu r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[p_{{{r}},\mu}+w_{{{r}},\mu}]\cdot\nu
        - \partial_\mu r_\mu^{\textnormal{du}}(u_{{{r}},\mu},p_{{{r}},\mu})[z_{{{r}},\mu}]\cdot\nu
        \end{aligned}$$ and thus $$\begin{aligned}
        \big(\nabla_\mu {\hat{\mathcal{J}}}_h(\mu) - \nabla_{\mu} {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)\big)_i
        &=d_{\mu_i}{\mathcal{J}}(u_{h, \mu}, \mu) - d_{\mu_i}{\mathcal{J}}(u_{{{r}}, \mu}, \mu)  + d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{h,\mu})[p_{h,\mu}] - d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[p_{{{r}},\mu}] \\
        & \qquad \qquad - d_{\mu_i} r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[w_{{{r}},\mu}] - d_{\mu_i} r_\mu^{\textnormal{du}}(u_{{{r}},\mu},p_{{{r}},\mu})[z_{{{r}},\mu}]
        \end{aligned}$$ The first line is equal to the estimator $\Delta_{{\widetilde\nabla}{\hat{J}}_{{r}}}(\mu)$, the first term of the second line can be estimated by $$\begin{aligned}
        d_{\mu_i} r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[w_{{{r}},\mu}] &\leq {\gamma_{d_{\mu_i}l_{\mu}}}\|w_{\mu}\| + {\gamma_{d_{\mu_i}a_{\mu}}} \| u_{{{r}},\mu} \| \|w_{\mu}\| 
        \end{aligned}$$ The second term can analogously be estimated by $$\begin{aligned}
        d_{\mu_i} r_\mu^{\textnormal{du}}(u_{{{r}},\mu},p_{{{r}},\mu})[z_{\mu}] \leq {\gamma_{d_{\mu_i}j}} \|z_{\mu}\| + 2{\gamma_{d_{\mu_i}k}} \| z_{\mu} \| \|u_{{{r}},\mu}\|
        + {\gamma_{d_{\mu_i}a}} \|z_{\mu}\| \|p_{{{r}},\mu}\|.
        \end{aligned}$$ We also have $$\begin{aligned}
        \underline{{a_{\mu}}} \|w_{\mu}\|^2 &\leq {a_{\mu}}(w_{\mu},w_{\mu}) = r_\mu^{\textnormal{du}}(u_{{{r}},\mu},p_{{{r}},\mu})[w_{\mu}]-2k_\mu(z_{\mu},w_{\mu}) 
        \leq \|r_\mu^{\textnormal{du}}(u_{{{r}},\mu},p_{{{r}},\mu})\| \|w_{\mu}\| + 2{\gamma_{k_\mu}} \|z_{\mu}\| \|w_{\mu}\|
        \end{aligned}$$ which gives $$\|w_{\mu}\| \leq \underline{{a_{\mu}}}^{-1} \left( \|r_\mu^{\textnormal{du}}(u_{{{r}},\mu},p_{{{r}},\mu})\| + 2{\gamma_{k_\mu}} \|z_{\mu}\| \right).$$ For $z_\mu$ we estimate $$\begin{aligned}
        \underline{{a_{\mu}}} \|z_{\mu}\|^2 &\leq {a_{\mu}}(z_{\mu},z_{\mu}) = -r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[z_\mu]\leq \|r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})\| \|z_{\mu}\|.
        \end{aligned}$$ Summing all together gives the assertion.

In a view of Section \[sec:ncd\_approach\], we emphasize that the estimator for the NCD-corrected gradient does not show a better approximation of the FOM gradient since more terms are added to the standard estimate. Propositon \[prop:Jhat\_error\].(ii) suggests that there exist an estimator of higher order which we derive in the following section.

### Sensitivity based approximation and estimation {#sec:sens_estimation}

We elaborate a better estimator for the NCD-corrected gradient by using sensitivities of the reduced primal and dual solutions. In addition, approximated sensitivities that are computed from the FOM sensitivities suggest an even better approximation of the FOM gradient.

We define the derivatives of the primal and dual solution maps associated with in direction $\nu \in \mathbb{R}^P$ as the solutions $d_{\nu} u_{{{r}}, \mu} \in V_{{r}}^{\textnormal{pr}}$ and $d_{\nu} p_{{{r}}, \mu} \in V_{{r}}^{\textnormal{du}}$ of $$\begin{aligned}
a_\mu(d_{\nu} u_{{{r}}, \mu}, v_{{r}}) &= \partial_{\mu} r_\mu^{\textnormal{pr}}(u_{{{r}}, \mu})[v_{{r}}]\cdot\nu &&\text{for all } v_{{r}}\in V_{{r}}^{\textnormal{pr}}\text{ and}\label{eq:true_primal_reduced_sensitivity}\\
a_\mu(q_{{r}}, d_{\nu} p_{{{r}}, \mu}) &= -\partial_\mu a_\mu(q_{{r}}, p_{r, \mu})\cdot\nu + d_\mu\partial_u {\mathcal{J}}(u_{{{r}}, \mu}, \mu)[q_{{r}}]\cdot\nu\notag\\
&= \partial_\mu r_\mu^{\textnormal{du}}(u_{{{r}}, \mu}, p_{{{r}}, \mu})[q_{{r}}]\cdot\nu + 2k_\mu(q_{{r}}, d_{\nu} u_{{{r}}, \mu}) &&\text{for all } q_{{r}}\in V_{{r}}^{\textnormal{du}},
\label{eq:true_dual_reduced_sensitivity}\end{aligned}$$ respectively, analogously to Propositions \[prop:solution\_dmu\_eta\] and \[prop:dual\_solution\_dmu\_eta\], where the last equality holds for quadratic functionals as in . With these sensitivities we can compute the same gradient of the NCD-corrected RB reduced functional from Propostion \[prop:true\_corrected\_reduced\_gradient\_adj\] in a different manner.

\[prop:true\_corrected\_reduced\_gradient\] The $i$-th component of the true gradient of ${{{{\hat{\mathcal{J}}}_{{r}}}}}$, $\nabla_\mu{{{{\hat{\mathcal{J}}}_{{r}}}}}:{\mathcal{P}}\to {\mathbb{R}}^P$, is given by $$\begin{aligned}
    \big(\nabla_\mu {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)\big)_i & = d_{\mu_i}{\mathcal{J}}(u_{{{r}}, \mu}, \mu) + d_{\mu_i} r_\mu^{\textnormal{pr}}(u_{{{r}}, \mu})[p_{{{r}}, \mu}] 
    + r_\mu^{\textnormal{pr}}(u_{{{r}}, \mu})[d_{\mu_i}p_{{{r}}, \mu}]
    + r_\mu^{\textnormal{du}}(u_{{{r}}, \mu},p_{{{r}}, \mu})[d_{\mu_i} u_{{{r}}, \mu}] 
    \end{aligned}$$ for all $1 \leq i \leq P$ and $\mu \in {\mathcal{P}}$, where $u_{{{r}}, \mu} \in V_{{r}}^{\textnormal{pr}}$ and $p_{{{r}}, \mu} \in V_{{r}}^{\textnormal{du}}$ solve , $d_{\mu_i}u_{{{r}}, \mu} \in V_{{r}}^{\textnormal{pr}}$ and $d_{\mu_i}p_{{{r}}, \mu} \in V_{{r}}^{\textnormal{du}}$ denote the derivatives of RB primal and dual solution maps from and .

It follows from the chain rule and Remark \[prop:gateaux\_wrt\_V\]; see [@HPUU2009 Section 1.6.1].

Note that the sensitivity based gradient is mathematically equivalent to the one in Propostion \[prop:true\_corrected\_reduced\_gradient\_adj\], but the second only requires to solve and once, because they can be reused for every component $\mu_i$, whereas the computation of the gradient in Propostion \[prop:true\_corrected\_reduced\_gradient\] requires to solve and for each $1\leq i\leq P$; cf. [@HPUU2009].

In terms of numerical approximation w.r.t. the FOM functional, we note that, e.g., a solution $d_{\mu_i} u_{{{r}}, \mu} \in V_{{r}}^{\textnormal{pr}}$ of does not necessarily need to be a good approximation of the FOM version $d_{\mu_i} u_{h, \mu} \in V_h$ even though $u_{h, \mu}$ is contained in $V_{{r}}^{\textnormal{pr}}$ since the high dimensional sensitivities are not generally contained in the respective reduced space (c.f. Propostion \[prop:primal\_solution\_dmui\_error\]).

To remedy this we could compute the FOM sensitivities for all canonical directions and either include them in the respective primal and dual space (thus forming Taylor RB spaces) or distribute all directional sensitivities to problem adapted RB spaces for the primal and dual sensitivities w.r.t. all canonical directions: $V_{{r}}^{{\textnormal{pr}},d_{\mu_i}}, V_{{r}}^{{\textnormal{du}},d_{\mu_i}} \subset V_h$. Thus, we again commit a variational crime.

\[prop:red\_sens\] Considering the reduced primal and dual solution maps ${\mathcal{P}}\to V_{{r}}^{\textnormal{pr}}$, $\mu \mapsto u_{{{r}}, \mu}$ and ${\mathcal{P}}\to V_{{r}}^{\textnormal{du}}$, $\mu \mapsto p_{{{r}}, \mu}$, respectively, where $u_{{{r}}, \mu}$ and $p_{{{r}}, \mu}$ are the solutions of and , we define their approximate partial derivatives w.r.t. the $i$th component of $\mu$ by ${\tilde d_{\mu_i}} u_{{{r}}, \mu} \in V_{{r}}^{{\textnormal{pr}},d_{\mu_i}}$ and ${\tilde d_{\mu_i}} p_{{{r}}, \mu} \in V_{{r}}^{{\textnormal{du}},d_{\mu_i}}$, respectively, as solutions of the sensitivity equations $$\begin{aligned}
    {a_{\mu}}({\tilde d_{\mu_i}} u_{{{r}}, \mu}, v_{{r}}) &= \partial_\mu {r_{\mu}}^{\textnormal{pr}}(u_{{{r}}, \mu})[v_{{r}}] \cdot e_i&&\text{for all } v_{{r}}\in V_{{r}}^{{\textnormal{pr}},d_{\mu_i}},
    \label{eq:red_sens_pr} \\
    {a_{\mu}}(q_{{r}}, {\tilde d_{\mu_i}} p_{{{r}}, \mu}) &= \partial_\mu {r_{\mu}}^{\textnormal{du}}(u_{{{r}}, \mu}, p_{{{r}}, \mu})[q_{{r}}] \cdot e_i 
    + 2{k_{\mu}}(q_{{r}}, {\tilde d_{\mu_i}} u_{{{r}}, \mu})&&\text{for all }q_{{r}}\in V_{{r}}^{{\textnormal{du}},d_{\mu_i}}.
    \label{eq:red_sens_du}
    \end{aligned}$$ Similarly, we denote the approximate partial derivatives in direction $\nu \in \mathbb{R}^P$ by ${\tilde d_{\nu}} u_{{{r}}, \mu}$ and ${\tilde d_{\nu}} p_{{{r}}, \mu}$, respectively, defined by substituting $e_i$ with $\nu$ above.

Following Propositions \[prop:solution\_dmu\_eta\] and \[prop:dual\_solution\_dmu\_eta\] we would obtain ${\tilde d_{\mu_i}} u_{{{r}}, \mu} = d_{\mu_i} u_{{{r}}, \mu}$, if $V_{{r}}^{{\textnormal{pr}},d_{\mu_i}} =V_{{r}}^{{\textnormal{pr}}}$ and ${\tilde d_{\mu_i}} p_{{{r}}, \mu} = d_{\mu_i} p_{{{r}}, \mu}$, if $V_{{r}}^{{\textnormal{du}},d_{\mu_i}} =V_{{r}}^{{\textnormal{du}}}$. Moreover, the approximate partial derivatives depend on the choice of the corresponding reduced approximation spaces.

\[prop:red\_approx\_grad\] We define the approximate gradient ${\widetilde\nabla_{\mu}}{{{{\hat{\mathcal{J}}}_{{r}}}}}:{\mathcal{P}}\to {\mathbb{R}}^P$ of ${{{{\hat{\mathcal{J}}}_{{r}}}}}$ by $$\begin{aligned}
    \big({\widetilde\nabla_{\mu}} {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)\big)_i :=  d_{\mu_i}{\mathcal{J}}(u_{{{r}},\mu}, \mu) &+ d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[p_{{{r}},\mu}] 
    + r_\mu^{\textnormal{pr}}(u_{{{r}}, \mu})[{\tilde d_{\mu_i}}p_{{{r}}, \mu}] + r_\mu^{\textnormal{du}}(u_{{{r}}, \mu},p_{{{r}}, \mu})[{\tilde d_{\mu_i}} u_{{{r}}, \mu}], 
    \end{aligned}$$ for $1 \leq i \leq P$, where $u_{{{r}},\mu} \in V_{{r}}^{\textnormal{pr}}$, $p_{{{r}},\mu} \in V_{{r}}^{\textnormal{du}}$ denote the reduced primal and dual solutions and ${\tilde d_{\mu_i}} u_{{{r}}, \mu} \in V_{{r}}^{{\textnormal{pr}},d_{\mu_i}}$ and ${\tilde d_{\mu_i}} p_{{{r}}, \mu} \in V_{{r}}^{{\textnormal{du}},d_{\mu_i}}$ denote the solutions of and .

Both gradients from Definition \[prop:true\_corrected\_reduced\_gradient\] and Propostion \[prop:red\_approx\_grad\] yield higher order estimate. To show this, we first derive error estimates for the reduction error of the reduced sensitivities from and as well as for and . For $v_h \in V_h$, the residuals of the equation in Propostion \[prop:solution\_dmu\_eta\] and Propostion \[prop:dual\_solution\_dmu\_eta\] for the canonical directions are respectively given by $$\begin{aligned}
{r_{\mu}}^{{\textnormal{pr}},d_{\mu_i}}(u_{h, \mu}, d_{\mu_i} u_{h, \mu})[v_h] &:= d_{\mu_i} {r_{\mu}}^{\textnormal{pr}}(u_{h, \mu})[v_h] - {a_{\mu}}(d_{\mu_i} u_{h, \mu}, v_h),
\label{sens_res_pr}\\
{r_{\mu}}^{{\textnormal{du}},d_{\mu_i}}(u_{h, \mu}, p_{h, \mu}, d_{\mu_i} u_{h, \mu}, d_{\mu_i} p_{h, \mu})[v_h] 
&:=  d_{\mu_i} {r_{\mu}}^{\textnormal{du}}(u_{h, \mu}, p_{h, \mu})[q_{{r}}] + 2{k_{\mu}}(q_{{r}}, d_{\mu_i} u_{h, \mu}) - {a_{\mu}}(q_h, d_{\mu_i} p_{h, \mu}) .
\label{sens_res_du}\end{aligned}$$

\[prop:primal\_solution\_dmui\_error\] For $\mu \in {\mathcal{P}}$ and $1 \leq i \leq P$, let $d_{\mu_i}u_{h, \mu} \in V_h$ be the solution of the discrete version of and $d_{\mu_i}u_{{{r}}, \mu} \in V_{{r}}^{{\textnormal{pr}},d_{\mu_i}}$ be the solution of . We then have $$\begin{aligned}
    \|d_{\mu_i}u_{h, \mu} - d_{\mu_i}u_{{{r}}, \mu}\| &\leq \Delta_{d_{\mu_i}{\textnormal{pr}}}(\mu) := \underline{{a_{\mu}}}^{-1}\Big({\gamma_{d_{\mu_i} {a_{\mu}}}} \Delta_{\textnormal{pr}}(\mu) + \|{r_{\mu}}^{{\textnormal{pr}},d_{\mu_i}}(u_{{{r}}, \mu}, d_{\mu_i}u_{{{r}}, \mu})\| \Big).
    \notag
    \end{aligned}$$

Using the shorthand $d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}:= d_{\mu_i}u_{h, \mu} - d_{\mu_i}u_{{{r}}, \mu}$, we obtain $$\begin{aligned}
    \underline{{a_{\mu}}}\, \|d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}\|^2 &\leq {a_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}, d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}) 
    = {a_{\mu}}(d_{\mu_i}u_{h, \mu}, d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}) -\, {a_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}) \\
    &= d_{\mu_i} {r_{\mu}}^{\textnormal{pr}}(u_{h, \mu})[d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}] - {a_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}) \\
    &= d_{\mu_i} {r_{\mu}}^{\textnormal{pr}}(u_{h, \mu})[d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}] - d_{\mu_i} {r_{\mu}}^{\textnormal{pr}}(u_{{{r}}, \mu})[d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}] + d_{\mu_i} {r_{\mu}}^{\textnormal{pr}}(u_{{{r}}, \mu})[d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}] - {a_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}) \\   
    &= - d_{\mu_i} {a_{\mu}}(e_{h, \mu}^{\textnormal{pr}}, d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}) + {r_{\mu}}^{{\textnormal{pr}},d_{\mu_i}}(u_{{{r}}, \mu}, d_{\mu_i}u_{{{r}}, \mu})[d_{\mu_i} e_{{{r}}, \mu}^{\textnormal{pr}}] \\
    &\leq {\gamma_{d_{\mu_i} {a_{\mu}}}} \;\|e_{h, \mu}^{\textnormal{pr}}\| \;\|d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}\|  + \|{r_{\mu}}^{{\textnormal{pr}},d_{\mu_i}}(u_{{{r}}, \mu}, d_{\mu_i}u_{{{r}}, \mu})\|\; \|d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}\|
    \end{aligned}$$ using the coercivity of ${a_{\mu}}$ in the first inequality, the definition $d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}$ in the first equality, Propostion \[prop:solution\_dmu\_eta\] applied to $u_{h, \mu}$ in the second equality, the definition of the discrete sensitivity primal residual in the third equality and the continuity of $d_{\mu_i} {a_{\mu}}$ in the last inequality.

We emphasize that the same result can be shown for ${\tilde d_{\mu_i}} u_{{{r}}, \mu}$ by replacing $d_{\mu_i} u_{{{r}}, \mu}$ and using the equation instead of . We call the resulting error estimator $\Delta_{{\tilde d_{\mu_i}}{\textnormal{pr}}}(\mu)$.

\[prop:dual\_solution\_dmui\_error\] For $\mu \in {\mathcal{P}}$ and $1 \leq i \leq P$, let $d_{\mu_i}p_{h, \mu} \in V_h$ be the solution of the discrete version of and $d_{\mu_i}p_{{{r}}, \mu} \in V_{{r}}^{{\textnormal{pr}},d_{\mu_i}}$ be the solution of . We then obtain $$\begin{aligned}
    \|d_{\mu_i}p_{h, \mu} - d_{\mu_i}p_{{{r}}, \mu}\| \leq \Delta_{ d_{\mu_i}{\textnormal{du}}}(\mu)\hspace*{5cm}\text{with}\\
    \Delta_{ d_{\mu_i}{\textnormal{du}}}(\mu) := \underline{{a_{\mu}}}^{-1}\Big(
    2 {\gamma_{d_{\mu_i} {k_{\mu}}}} \; \Delta_{\textnormal{pr}}(\mu) +  {\gamma_{d_{\mu_i} {a_{\mu}}}} \; \Delta_{\textnormal{du}}(\mu) + 2 {\gamma_{{k_{\mu}}}}  \; \Delta_{d_{\mu_i}pr}(\mu) 
    + \| {r_{\mu}}^{{\textnormal{du}},d_{\mu_i}}(u_{{{r}}, \mu}, p_{{{r}}, \mu}, d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i}p_{{{r}}, \mu}) \| \Big).
    \end{aligned}$$

Using the shorthand $d_{\mu_i} e_{h, \mu}^{\textnormal{du}}:= d_{\mu_i}p_{h, \mu} - d_{\mu_i}p_{{{r}}, \mu}$ and $e_{h, \mu}^{\textnormal{du}}:= p_{h, \mu} - p_{{{r}}, \mu}$, we obtain $$\begin{aligned}
    \underline{{a_{\mu}}}\;&\|d_{\mu_i} e_{h, \mu}^{\textnormal{du}}\|^2 \leq {a_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, d_{\mu_i} e_{h, \mu}^{\textnormal{du}})
    = \underbrace{{a_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, d_{\mu_i}p_{h, \mu})}_{= d_{\mu_i} {r_{\mu}}^{\textnormal{du}}(u_{h, \mu}, p_{h, \mu})[d_{\mu_i} e_{h, \mu}^{\textnormal{du}}] 
        + 2{k_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, d_{\mu_i}u_{h, \mu})} - {a_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, d_{\mu_i}p_{{{r}}, \mu})
    \\
    &= d_{\mu_i} {r_{\mu}}^{\textnormal{du}}(u_{h, \mu}, p_{h, \mu})[d_{\mu_i} e_{h, \mu}^{\textnormal{du}}] + 2{k_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, d_{\mu_i}u_{h, \mu})
    - d_{\mu_i} {r_{\mu}}^{\textnormal{du}}(u_{{{r}}, \mu}, p_{{{r}}, \mu})[d_{\mu_i} e_{h, \mu}^{\textnormal{du}}] \\
    &\qquad - 2{k_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, d_{\mu_i}u_{{{r}}, \mu})
    + d_{\mu_i} {r_{\mu}}^{\textnormal{du}}(u_{{{r}}, \mu}, p_{{{r}}, \mu})[d_{\mu_i} e_{h, \mu}^{\textnormal{du}}] + 2{k_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, d_{\mu_i}u_{{{r}}, \mu}) 
     - {a_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, d_{\mu_i}p_{{{r}}, \mu}) \\
    &= d_{\mu_i} {j_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}) + 2 d_{\mu_i}{k_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, u_{h, \mu}) - d_{\mu_i}{a_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, p_h)
     - d_{\mu_i}{j_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}) + 2 d_{\mu_i}{k_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, u_{{{r}}, \mu})\\
    &\qquad 
    - d_{\mu_i}{a_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, p_{{r}}) +\! 2{k_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, d_{\mu_i}u_{h, \mu}) \!-\! 2{k_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, d_{\mu_i}u_{{{r}}, \mu}) \\
    &\qquad + {r_{\mu}}^{{\textnormal{du}},d_{\mu_i}}(u_{{{r}}, \mu}, p_{{{r}}, \mu}, d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i}p_{{{r}}, \mu})[d_{\mu_i} e_{h, \mu}^{\textnormal{du}}] \\
    &= 2 d_{\mu_i}{k_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, e_{h, \mu}^{\textnormal{pr}}) - d_{\mu_i}{a_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, e_{h, \mu}^{\textnormal{du}})
     + 2{k_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{du}}, d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}) \\
    &\qquad
    + {r_{\mu}}^{{\textnormal{du}},d_{\mu_i}}(u_{{{r}}, \mu}, p_{{{r}}, \mu}, d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i}p_{{{r}}, \mu})[d_{\mu_i} e_{h, \mu}^{\textnormal{du}}]\\
    &\leq \big(2 {\gamma_{d_{\mu_i} {k_{\mu}}}} \; \|e_{h, \mu}^{\textnormal{pr}}\| +  {\gamma_{d_{\mu_i} {a_{\mu}}}} \; \|e_{h, \mu}^{\textnormal{du}}\|  \big)\|d_{\mu_i} e_{h, \mu}^{\textnormal{du}}\| \\
    &\qquad + 2 {\gamma_{{k_{\mu}}}}  \; \|d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}\| \; \|d_{\mu_i} e_{h, \mu}^{\textnormal{du}}\|
    + \| {r_{\mu}}^{{\textnormal{du}},d_{\mu_i}}(u_{{{r}}, \mu}, p_{{{r}}, \mu}, d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i}p_{{{r}}, \mu}) \| \; \|d_{\mu_i} e_{h, \mu}^{\textnormal{du}}\|
    \end{aligned}$$ using the coercivity of ${a_{\mu}}$ in the first inequality, the definition of $d_{\mu_i} e_{h, \mu}^{\textnormal{du}}$ in the first equality, Propostion \[prop:dual\_solution\_dmu\_eta\] applied to $p_{h, \mu}$ in the second equality, the definition of the dual residual in in the third equality and continuity of all parts in the last inequality.

Again, the same result holds for ${\tilde d_{\mu_i}} p_{{{r}}, \mu}$ if we replace $d_{\mu_i} p_{{{r}}, \mu}$ and use instead of . The resulting error estimator is then called $\Delta_{{\tilde d_{\mu_i}}{\textnormal{du}}}(\mu)$.

Using the residual based a posteriori error estimates for the primal sensitivities, we are able to state two a posteriori error bounds on the model reduction error of the true gradient and the approximated gradient of the NCD-corrected functional.

\[prop:grad\_Jhat\_error\_sens\]

(i) For the gradient $\nabla_{\mu} {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)$ of the NCD-corrected RB reduced functional, computed with sensitivities according to Proposition [[\[prop:true\_corrected\_reduced\_gradient\]]{}]{}, we have $$\begin{aligned}
          \big\|\nabla_\mu {\hat{\mathcal{J}}}_h(\mu) - \nabla_{\mu} {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)\big\|_2 &\leq \Delta_{\nabla{\hat{\mathcal{J}}}_{{r}}}(\mu) 
        = \big\|\underline{\Delta_{\nabla{\hat{\mathcal{J}}}_{{r}}}(\mu)}\big\|_2 \quad\quad\text{with} \\
          \big(\underline{\Delta_{\nabla{\hat{\mathcal{J}}}_{{r}}}(\mu)}\big)_i := {\gamma_{d_{\mu_i} {k_{\mu}}}} \, \left(\Delta_{\textnormal{pr}}(\mu)\right)^2 &+ {\gamma_{{a_{\mu}}}} \, 
            \Delta_{d_{\mu_i}pr}(\mu) \,\Delta_{\textnormal{du}}(\mu)  + \| {r_{\mu}}^{{\textnormal{du}},d_{\mu_i}}(u_{{{r}}, \mu}, p_{{{r}}, \mu}, d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i}p_{{{r}}, \mu}) \| \,
            \Delta_{\textnormal{pr}}(\mu).
            \end{aligned}$$

(ii) Furthermore, we have for the approximate gradient from Definition [[\[prop:red\_approx\_grad\]]{}]{} $$\begin{aligned}
            \big\|\nabla_\mu {\hat{\mathcal{J}}}_h(\mu) - {\widetilde\nabla_{\mu}} {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)\big\|_2 &\leq \Delta_{\widetilde\nabla{\hat{\mathcal{J}}}_{{r}}}(\mu) 
        = \big\|\underline{\Delta_{\widetilde\nabla{\hat{\mathcal{J}}}_{{r}}}(\mu)}\big\|_2 \quad\quad\text{with} \\
        \big(\underline{\Delta_{\widetilde\nabla{\hat{\mathcal{J}}}_{{r}}}(\mu)}\big)_i := {\gamma_{d_{\mu_i} {k_{\mu}}}} \, \left(\Delta_{\textnormal{pr}}(\mu)\right)^2 &+ {\gamma_{{a_{\mu}}}} \, 
            \Delta_{{\tilde d_{\mu_i}}{\textnormal{pr}}}(\mu) \,\Delta_{\textnormal{du}}(\mu) + \| {r_{\mu}}^{{\textnormal{du}},d_{\mu_i}}(u_{{{r}}, \mu}, &p_{{{r}}, \mu}, {\tilde d_{\mu_i}}u_{{{r}}, \mu}, {\tilde d_{\mu_i}}p_{{{r}}, \mu}) \| \,
            \Delta_{\textnormal{pr}}(\mu).
            \end{aligned}$$

\(i) To prove the first assertion, we use $r_\mu^{\textnormal{pr}}(u_{h, \mu})[d_{\mu_i}p_{{{r}}, \mu}] = 0$ and $r_\mu^{\textnormal{du}}(u_{h, \mu},p_{h, \mu})[d_{\mu_i} u_{{{r}}, \mu}]=0$ to obtain $$\begin{aligned}
        \big(\nabla_\mu {\hat{\mathcal{J}}}_h(\mu) &- \nabla_\mu {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)\big)_i 
        =d_{\mu_i}{\mathcal{J}}(u_{h, \mu}, \mu) - d_{\mu_i}{\mathcal{J}}(u_{{{r}}, \mu}, \mu) + d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{h,\mu})[p_{h,\mu}] - d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[p_{{{r}},\mu}]  \\
        &= \partial_\mu {j_{\mu}}(e_{h, \mu}^{\textnormal{pr}}) + \partial_\mu {k_{\mu}}(u_{h, \mu}, u_{h, \mu}) - \partial_\mu {k_{\mu}}(u_{{{r}}, \mu},u_{{{r}}, \mu})  
         + d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{h,\mu})[p_{h,\mu}] - d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[p_{{{r}},\mu}] \\
        &\qquad + \underset{=(*)}{\underbrace{r_\mu^{\textnormal{pr}}(e_{h, \mu}^{\textnormal{pr}})[d_{\mu_i}p_{{{r}}, \mu}]}} 
        + \underset{=(**)}{\underbrace{r_\mu^{\textnormal{du}}(e_{h, \mu}^{\textnormal{pr}},e_{h, \mu}^{\textnormal{du}})[d_{\mu_i} u_{{{r}}, \mu}]}}.
        \end{aligned}$$ For the last two residual terms we have $$\begin{aligned}
        (*) &= {l_{\mu}}(d_{\mu_i}p_{{{r}}, \mu}) - {l_{\mu}}(d_{\mu_i}p_{{{r}}, \mu}) 
        - {a_{\mu}}(e_{h, \mu}^{\textnormal{pr}},d_{\mu_i}p_{{{r}}, \mu}) \\ 
        &= - {a_{\mu}}(e_{h, \mu}^{\textnormal{pr}},d_{\mu_i}p_{{{r}}, \mu}) 
        + d_{\mu_i} r_\mu^{\textnormal{du}}(u_{{{r}}, \mu},p_{{{r}}, \mu})[e_{h, \mu}^{\textnormal{pr}}] 
        + 2{k_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}, e_{h, \mu}^{\textnormal{pr}}) 
        - d_{\mu_i} r_\mu^{\textnormal{du}}(u_{{{r}}, \mu},p_{{{r}}, \mu})[e_{h, \mu}^{\textnormal{pr}}] 
        - 2{k_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}, e_{h, \mu}^{\textnormal{pr}}) \\
        &= {r_{\mu}}^{{\textnormal{du}},d_{\mu_i}}(u_{{{r}}, \mu}, p_{{{r}}, \mu}, d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i}p_{{{r}}, \mu})[e_{h, \mu}^{\textnormal{pr}}] 
        - d_{\mu_i} r_\mu^{\textnormal{du}}(u_{{{r}}, \mu},p_{{{r}}, \mu})[e_{h, \mu}^{\textnormal{pr}}] 
        - 2{k_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}, e_{h, \mu}^{\textnormal{pr}}) 
        \end{aligned}$$ and $$\begin{aligned}
        (**) &= {j_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}) - {j_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}) 
        + 2{k_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}, e_{h, \mu}^{\textnormal{pr}}) - {a_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}, e_{h, \mu}^{\textnormal{du}}). 
        \end{aligned}$$ Thus, by summing both terms we have $$\begin{aligned}
        (*) + (**)= 
        {r_{\mu}}^{{\textnormal{du}},d_{\mu_i}}(u_{{{r}}, \mu}, p_{{{r}}, \mu},d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i}p_{{{r}}, \mu})[e_{h, \mu}^{\textnormal{pr}}] 
        - \underset{=(***)}{\underbrace{d_{\mu_i} r_\mu^{\textnormal{du}}(u_{{{r}}, \mu},p_{{{r}}, \mu})[e_{h, \mu}^{\textnormal{pr}}]}} - {a_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}, e_{h, \mu}^{\textnormal{du}}) 
        \end{aligned}$$ and for $(***)$ it holds $$\begin{aligned}
        (***) = d_{\mu_i}{j_{\mu}}(e_{h, \mu}^{\textnormal{pr}}) + 2 d_{\mu_i}{k_{\mu}}(e_{h, \mu}^{\textnormal{pr}}, u_{{{r}},\mu}) - d_{\mu_i}{a_{\mu}}(e_{h, \mu}^{\textnormal{pr}},p_{{{r}},\mu}).
        \end{aligned}$$ Combining $(*)$, $(**)$ and $(***)$ with the previous result, we have $$\begin{aligned}
        \big(\nabla_\mu {\hat{\mathcal{J}}}_h(\mu) &- \nabla_\mu {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)\big)_i 
         = \partial_\mu {j_{\mu}}(e_{h, \mu}^{\textnormal{pr}}) + \partial_\mu {k_{\mu}}(u_{h, \mu}, u_{h, \mu}) - \partial_\mu {k_{\mu}}(u_{{{r}}, \mu},u_{{{r}}, \mu})  \\
        &\qquad - d_{\mu_i}{j_{\mu}}(e_{h, \mu}^{\textnormal{pr}}) - 2 d_{\mu_i}{k_{\mu}}(e_{h, \mu}^{\textnormal{pr}}, u_{{{r}},\mu}) + d_{\mu_i}{a_{\mu}}(e_{h, \mu}^{\textnormal{pr}},p_{{{r}},\mu}) + d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{h,\mu})[p_{h,\mu}]  \\
        &\qquad - d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[p_{{{r}},\mu}] 
         +  {r_{\mu}}^{{\textnormal{du}},d_{\mu_i}}(u_{{{r}}, \mu}, p_{{{r}}, \mu},d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i}p_{{{r}}, \mu})[e_{h, \mu}^{\textnormal{pr}}] 
        - {a_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}, e_{h, \mu}^{\textnormal{du}})  \\
        & = d_{\mu_i} {k_{\mu}}(e_{h, \mu}^{\textnormal{pr}}, e_{h, \mu}^{\textnormal{pr}}) 
        +  {r_{\mu}}^{{\textnormal{du}},d_{\mu_i}}(u_{{{r}}, \mu}, p_{{{r}}, \mu},d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i}p_{{{r}}, \mu})[e_{h, \mu}^{\textnormal{pr}}] \\
        & \qquad + \underset{=(****)}{\underbrace{d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{h,\mu})[p_{h,\mu}] - d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[p_{{{r}},\mu}] 
                +  d_{\mu_i}{a_{\mu}}(e_{h, \mu}^{\textnormal{pr}},p_{{{r}},\mu}) -  {a_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}, e_{h, \mu}^{\textnormal{du}})}}.
        \end{aligned}$$ Further, we have $$\begin{aligned}
        d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{h,\mu})[p_{h,\mu}] &- d_{\mu_i}r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[p_{{{r}},\mu}] = d_{\mu_i} {l_{\mu}}(e_{h, \mu}^{\textnormal{du}}) - 
        d_{\mu_i} {a_{\mu}}(u_{h, \mu}, p_{h, \mu}) + d_{\mu_i} {a_{\mu}}(u_{{{r}}, \mu}, p_{{{r}}, \mu}) \\
        & =  {a_{\mu}}(d_{\mu_i}u_h, e_{h, \mu}^{\textnormal{du}}) + d_{\mu_i} {a_{\mu}}(u_h, e_{h, \mu}^{\textnormal{du}}) 
         - d_{\mu_i} {a_{\mu}}(u_{h, \mu}, p_{h, \mu}) + d_{\mu_i} {a_{\mu}}(u_{{{r}}, \mu}, p_{{{r}}, \mu}),
        \end{aligned}$$ where we used the discretized version of in the second equality. Inserting this into $(*)$ gives $$\begin{aligned}
        (***\,*) =&  {a_{\mu}}(d_{\mu_i}u_h, e_{h, \mu}^{\textnormal{du}}) -  {a_{\mu}}(d_{\mu_i}u_{{{r}}, \mu}, e_{h, \mu}^{\textnormal{du}}) \\ & +  
        \underset{=0}{\underbrace{d_{\mu_i} {a_{\mu}}(u_h, e_{h, \mu}^{\textnormal{du}}) - d_{\mu_i} {a_{\mu}}(u_{h, \mu}, p_{h, \mu}) 
                + d_{\mu_i} {a_{\mu}}(u_{{{r}}, \mu}, p_{{{r}}, \mu}) 
                +d_{\mu_i}{a_{\mu}}(e_{h, \mu}^{\textnormal{pr}},p_{{{r}},\mu}) }} 
        = {a_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}, e_{h, \mu}^{\textnormal{du}}).
        \end{aligned}$$ In total, we have $$\begin{aligned}
        & &\big(\nabla_\mu {\hat{\mathcal{J}}}_h(\mu) - \nabla_\mu {{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)\big)_i 
         = d_{\mu_i} {k_{\mu}}(e_{h, \mu}^{\textnormal{pr}}, e_{h, \mu}^{\textnormal{pr}}) 
        + {a_{\mu}}(d_{\mu_i} e_{h, \mu}^{\textnormal{pr}}, e_{h, \mu}^{\textnormal{du}}) +  {r_{\mu}}^{{\textnormal{du}},d_{\mu_i}}(u_{{{r}}, \mu}, p_{{{r}}, \mu},d_{\mu_i}u_{{{r}}, \mu}, d_{\mu_i}p_{{{r}}, \mu})[e_{h, \mu}^{\textnormal{pr}}] 
        \end{aligned}$$ which proofs the assertion.\
(ii) The estimate follows analogously to (i), by replacing $d_{\mu_i} u_{{{r}}, \mu}$ and $d_{\mu_i} p_{{{r}}, \mu}$ with ${\tilde d_{\mu_i}} u_{{{r}}, \mu}$ and ${\tilde d_{\mu_i}} p_{{{r}}, \mu}$, respectively.

To conclude, $\Delta_{\widetilde\nabla{\hat{\mathcal{J}}}_{{r}}}(\mu)$ and $\Delta_{\nabla{\hat{\mathcal{J}}}_{{r}}}(\mu)$ both decay with second order (cf. Section \[sec:estimator\_study\]). We also point out, that $\Delta_{\widetilde\nabla{\hat{\mathcal{J}}}_{{r}}}(\mu)$ is an improved estimator which can be used to replace the poor estimator $\Delta^{{{r}},*}_{\nabla_\mu {\hat{\mathcal{J}}}_{{r}}}(\mu)$. However both higher order estimators $\Delta_{\widetilde\nabla{\hat{\mathcal{J}}}_{{r}}}(\mu)$ and $\Delta_{\nabla{\hat{\mathcal{J}}}_{{r}}}(\mu)$ come with the price of computing the dual norm of the sensitivity residuals in and for each direction which aggravates the computational complexity.

The Trust-Region Method and adaptive enrichment strategies {#sec:TRRB_and_adaptiveenrichment}
==========================================================

To solve problem we apply the TR method, which iteratively computes a first-order critical point of . At each iteration $k\geq 0$, we consider a so-called model function $m^{(k)}$, which is a cheaply computable approximation of the quadratic cost functional ${\mathcal{J}}$ in a neighbourhood of the parameter $\mu^{(k)}$, i.e., the Trust-Region. Therefore, for $k\geq 0$, given a TR radius $\delta^{(k)}$, we consider the TR minimization sub-problem $$\label{TRsubprob}
\min_{s\in \mathbb{R}^P} m^{(k)}(s) \, \text{ subject to } \|s\|_2 \leq \delta^{(k)},\, \widetilde{\mu}:= \mu^{(k)}+s \in{\mathcal{P}}\text{ and } r_{\tilde{\mu}}^{\textnormal{pr}}(u_{\tilde{\mu}})[v]= 0 \, \text{ for all }  v\in V.$$ Under suitable assumptions on $m^{(k)}$, problem admits a unique solution $\bar s^{(k)}$, which is used to compute the next iterate $\mu^{(k+1)} = \mu^{(k)} + \bar s^{(k)}$.

The Trust-Region Reduced Basis Method {#sec:TR}
-------------------------------------

Slightly different from [@AFS00; @QGVW2017], we choose as model function the NCD-corrected RB reduced functional ${{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}$ defined in (\[eq:Jhat\_red\_corected\]), i.e. $m^{(k)}(\cdot)= {{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}(\mu^{(k)}+\cdot)$ for $k\geq 0$, where the super-index $(k)$ indicates that we use different RB spaces $V_{{r}}^{*, (k)}$ in each iteration. As indicated in Proposition \[prop:Jhat\_error\] and shown in our numerical experiments below, ${{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}$ converges to ${\hat{\mathcal{J}}}$ with higher order in comparison to the standard RB reduced functional from , which has been considered in [@QGVW2017]. We initialize the RB spaces using the initial guess $\mu^{(0)}$, i.e. setting $V^{{\textnormal{pr}},0}_{{r}}= \left\{u_{h,\mu^{(0)}}\right\}$ and $V^{{\textnormal{du}},0}_{{r}}= \left\{p_{h,\mu^{(0)}}\right\}$. At every iteration $k$ we may – depending on the a posteriori estimates – enrich the obtained space using the computed parameter $\mu^{(k+1)}$; for further details see Section \[sec:construct\_RB\]. Possible sufficient and necessary conditions for convergence, dependent on the approximate generalized Cauchy point (AGC) $\mu^{(k)}_\text{\rm{AGC}}$ (see Definition \[Def:AGC\]), are given in [@QGVW2017]. In contrast to [@QGVW2017], we consider additional bilateral parameter constraints in . In particular, the presence of these inequality constraints requires a review of the proof of convergence for the TR-RB algorithm. In [@QGVW2017], the convergence is based on the results contained in [@YM2013], where the authors consider an equality-constrained optimization problem. We state first how our method differs from the one in [@QGVW2017], then we prove the convergence of this modified algorithm. According to [@QGVW2017], the inexact RB version of problem is $$\label{TRRBsubprob}
\min_{\widetilde{\mu}\in{\mathcal{P}}} {{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}(\widetilde{\mu}) \text{ s.t. }  \frac{\Delta_{{{{{\hat{\mathcal{J}}}_{{r}}}}}}(\widetilde{\mu})}{{{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}(\widetilde{\mu})}\leq \delta^{(k)},$$ where $\widetilde{\mu}:= \mu^{(k)}+s$, the equality constraint $r_{\tilde\mu}^{\textnormal{pr}}(u_{\tilde{\mu}})[v]= 0$ is hidden in the definition of ${{{{\hat{\mathcal{J}}}_{{r}}}}}$ and the inequality ones are concealed in the request $\widetilde{\mu}\in {\mathcal{P}}$. As also remarked in [@QGVW2017], the projected BFGS method [@Kel99], which we use in order to solve , computes the ACG point $\mu^{(k)}_\text{\rm{AGC}}$ in the first iterate and generates a sequence $\{\mu^{(k,\ell)}\}_{\ell=1}^L$ where $L$ is the last BFGS iteration. In what follows, $\mu^{(k,1)}:= \mu^{(k)}_\text{\rm{AGC}}$ and the TR iterate $\mu^{(k+1)}:=\mu^{(k,L)}$. Throughout the paper the index $k$ refers to the current outer TR iteration, $\ell$ refers instead to the inner BFGS iteration. Note that $L$ may be different for each iteration $k$, but we will indicate it only when strictly necessary in order to simplify the notation. To describe the projected BFGS method in details, we define $$\begin{aligned}
\label{eq:General_Opt_Step_point}
\mu^{(k,\ell)}(j):= {\text{P}}_\mathcal{P}(\mu^{(k,\ell)} + \kappa^j d^{(k,\ell)}) \in{\mathcal{P}}&& \text{for } j\geq 0,\end{aligned}$$ where $\kappa\in(0,1)$, $d^{(k,\ell)}\in\mathbb{R}^P$ is the chosen descent direction at the iteration $(k,\ell)$ and the projection operator ${\text{P}}_{\mathcal{P}}: \mathbb{R}^P\rightarrow {\mathcal{P}}$ is defined as $$\begin{aligned}
({\text{P}}_{\mathcal{P}}(\mu))_i:= \left\{ \begin{array}{ll}
(\mu_\mathsf{a})_i & \text{if } (\mu)_i\leq (\mu_\mathsf{a})_i, \\
(\mu_\mathsf{b})_i & \text{if } (\mu)_i\geq (\mu_\mathsf{b})_i, \\
(\mu)_i & \text{otherwise,}
\end{array} \right. 
&& \text{for } i=1,\ldots,P.\end{aligned}$$ Note that the operator ${\text{P}}_{\mathcal{P}}$ is Lipschitz continuous with constant one; cf. [@Kel99]. For computing the descent direction $d^{(k,\ell)}$ we follow the projected BFGS algorithm reported in [@Kel99 Section 5.5.3]. Furthermore, we enforce respectively an Armijo-type condition and the additional TR constraint on ${{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}$ $$\begin{aligned}
\label{Armijo}{{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}(\mu^{(k,\ell)}(j)) - {{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}(\mu^{(k,\ell)}) & \leq  -\frac{\kappa_{\mathsf{arm}}}{\kappa^j} \big\| \mu^{(k,\ell)}(j)-\mu^{(k,\ell)}\big\|^2_2, \\
\label{TR_radius_condition} \frac{\Delta_{{{{{\hat{\mathcal{J}}}_{{r}}}}}}(\mu^{(k,\ell)}(j))}{{{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}(\mu^{(k,\ell)}(j))} & \leq \delta^{(k)},\end{aligned}$$ by selecting $\mu^{(k,\ell+1)} = \mu^{(k,\ell)}(j^{(k,\ell)})$ for $\ell\geq 1$, where $j^{(k,\ell)}<\infty$ is the smallest index for which the conditions - hold for some $\kappa_{\mathsf{arm}}\in(0,\frac{1}{2})$, generally $\kappa_{\mathsf{arm}}=10^{-4}$; cf. [@QGVW2017]. Moreover, we use as termination criteria for the optimization sub-problem

\[Termination\_crit\_subproblem\]

$$\label{FOC_subproblem}
        \big\|\mu^{(k,\ell)}-{\text{P}}_{\mathcal{P}}(\mu^{(k,\ell)}-\nabla_\mu {\hat{\mathcal{J}}}_{{r}}^{(k)}(\mu^{(k,\ell)}))\big\|_2\leq \tau_\text{\rm{sub}}
        \hspace{2em} \text{or}$$

$$\label{Cut_of_TR}
        \beta_2\delta^{(k)} \leq \frac{\Delta_{{\hat{\mathcal{J}}}_{{r}}^{(k)}}(\mu)}{{\hat{\mathcal{J}}}_{{r}}^{(k)}(\mu)} \leq \delta^{(k)},$$

1em

where $\tau_\text{\rm{sub}}\in(0,1)$ is a predefined tolerance and $\beta_2\in(0,1)$, generally close to one. Condition is used to prevent that the optimizer spends much time close to the boundary of the Trust-Region, where the model is poor in approximation; cf. [@QGVW2017]. Note that, without the projection operator ${\text{P}}_{\mathcal{P}}$, conditions - coincide with the ones in [@QGVW2017], apart from using the NCD-corrected RB reduced functional. Furthermore, in addition to [@QGVW2017], we consider a condition which allows enlarging the TR radius. A drawback of the TR algorithm proposed in [@QGVW2017] is that the TR radius may be significantly shrunk at the beginning, i.e. when the TR model is poor in approximation. Afterwards, even if the RB space is enriched, i.e. the approximation of the TR model function is improved, the TR radius is kept small. Thus, one misses the local second-order rate of convergence of the BFGS method. More precisely, if $\mu^{(k,\ell)}$ is close to the locally optimal solution $\bar \mu^{(k)}$ of the TR sub-problem, we want to make full BFGS steps, which gives us faster convergence. The possibility to enlarge the TR radius at each iteration will also decrease the number of outer iterations needed to converge. As a condition for enlarging the radius we check whether the sufficient reduction predicted by the model function ${{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}$ is realized by the objective function, i.e. we check if $$\label{TR_act_decrease}
\varrho^{(k)}:= \frac{{\hat{\mathcal{J}}}_h(\mu^{(k)})-{\hat{\mathcal{J}}}_h(\mu^{(k+1)})}{{{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}(\mu^{(k)})-{{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}(\mu^{(k+1)})}  \geq  \eta_\varrho$$ for a tolerance $\eta_\varrho \in [3/4,1)$. Condition seems costly because of the evaluation of the FOM cost functional ${\hat{\mathcal{J}}}_h$, but, after the enrichment of the RB space, the quantities in the numerator of are cheaply accessible, since one has already solved the FOM to generate the new snapshots for the RB space enrichment. Note that this also implies that we can cheaply evaluate the FOM gradient $\nabla_\mu {\hat{\mathcal{J}}}_h(\mu^{k+1})$ in case of an enrichment. This knowledge will be used for the stopping criterium in the outer loop of the algorithm. Finally, let us define the AGC point for our constrained case.

\[Def:AGC\] At iteration $k$, we define the AGC point as $$\mu_\text{\rm{AGC}}^{(k)}:= \mu^{(k,0)}(j^{(k)}_c)=  {\text{P}}_\mathcal{P}(\mu^{(k,0)} + \kappa^{j^{(k)}_c} d^{(k,0)}),$$ where $\mu^{(k,0)}:= \mu^{(k)}$, $d^{(k,0)}:= -\nabla_\mu {{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}(\mu^{(k,0)})$ and $j^{(k)}_c$ is the smallest non-negative integer $j$ for which $\mu^{(k,0)}(j)$ satisfies - for $\ell=0$.

We refer to Algorithm \[Alg:TR-RBmethod\] for the proposed TR-RB algorithm.

Set $k=0$ and $=$

Convergence analysis {#sec:TR_convergence_analysis}
--------------------

In order to guarantee the well-posedness (because of ) and the convergence of the method, we make the following assumption

\[asmpt:bound\_J\] The cost functional ${\mathcal{J}}(u,\mu)$ is strictly positive for all $u\in V$ and all parameters $\mu\in{\mathcal{P}}$.

Note that this assumption is not too restrictive, since the boundedness from below is a standard assumption in optimization to guarantee the existence of a solution for the minimization problem. If a global lower bound for the cost functional is also known, one can add a sufficiently large constant, without changing the position of its local minima and maxima. Another important request, pointed out in [@QGVW2017; @YM2013], is that an error-aware sufficient decrease condition $$\begin{aligned}
\label{Suff_decrease_condition}
{{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k+1)}(\mu^{(k+1)})\leq {{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}(\mu_\text{\rm{AGC}}^{(k)}) && \text{ for all } k\in\mathbb{N}\end{aligned}$$ is fulfilled at each iteration $k$ of the TR-RB algorithm. As in [@QGVW2017; @YM2013], we consider cheaply computable sufficient and necessary conditions for in Algorithm \[Alg:TR-RBmethod\] (Step \[Suff\_condition\_TRRB\] and Step \[Nec\_condition\_TRRB\], respectively). The TR-RB algorithm rejects, then, any computed point which does not satisfy . One may be concerned of the fact that Algorithm \[Alg:TR-RBmethod\] may be trapped in an infinite loop where every computed point is rejected and the TR radius is shrunk all time. We point out that this never happened in our numerical tests. Anyway, we consider a safety termination criteria, which is triggered when the TR radius is smaller than the double machine precision. To prove convergence of Algorithm \[Alg:TR-RBmethod\], in what follows, we then assume that this situation can not occur.

\[asmpt:rejection\] For each $k\geq 0$, there exists a radius $\widetilde{\delta}^{(k)}>\tau_{\rm{\text{mac}}}$ for which a solution of is such that $\eqref{Suff_decrease_condition}$ is verified, where $\tau_{\rm{\text{mac}}}= 2.22\cdot10^{-16}$ is the double machine precision.

\[Lemma:AGC\_search\] Let Assumptions [[[\[asmpt:parameter\_separable\]]{}]{}–[[\[asmpt:rejection\]]{}]{}]{} hold true. The search of the AGC point defined in Definition [[\[Def:AGC\]]{}]{} takes finitely many iterations at each step $k$ of the TR-RB Algorithm.

We want to prove that there exists an index $j_c^{(k)}<\infty$ for each $k\geq 0$, for which $\mu_\text{\rm{AGC}}^{(k)}= \mu^{(k,0)}(j^{(k)}_c)$ satisfies - for $\ell=0$. From [@Kel99 Theorem 5.4.5] (and the subsequent discussion) we conclude that for all $k\in\mathbb{N}$ there exists a strictly positive index $j^{(k)}_1\in\mathbb{N}$ such that $\mu^{(k,0)}(j)$ satisfies for $j\geq j^{(k)}_1$ and $\ell=0$. If $k=0$, by construction we have that $\Delta_{{\hat{\mathcal{J}}}_{{r}}^{(0)}}(\mu^{(0)}) = 0$. Therefore, there exists a sufficiently large (but finite) index $j^{(0)}_2\in\mathbb{N}$ such that $\mu^{(0,0)}(j)$ satisfies for all $j\geq j^{(0)}_2$ and $\ell=0$. The reason relies on the continuity w.r.t. $\mu$ of the error estimator $\Delta_{{\hat{\mathcal{J}}}_{{r}}^{(k)}}(\mu)$ (cf. Remark \[continuity\_of\_estimator\]) and of the cost functional ${{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k)}(\mu)$ for all $k\in\mathbb{N}$. Hence there exists $j^{(0)}_c=\max(j^{(0)}_1,j^{(0)}_2)<\infty$, for which $\mu^{(0,0)}(j)$ satisfies - for $\ell=0$. If $k\geq 1$, since the model has been enriched, i.e. $\Delta_{{\hat{\mathcal{J}}}_{{r}}^{(k)}}(\mu^{(k)}) = 0$, we can show the claim arguing as we did for $k=0$. Note, in fact, that we increase the iteration counter only when $\mu^{(k)}$ is accepted at iteration $k-1$ and, thus, when the RB model is enriched at this parameter.

\[Thm:convergence\_of\_TR\] Let the hypotheses of Lemma [[\[Lemma:AGC\_search\]]{}]{} be verified. Then every accumulation point $\bar\mu$ of the sequence $\{\mu^{(k)}\}_{k\in\mathbb{N}}\subset {\mathcal{P}}$ generated by the TR-RB algorithm is an approximate first-order critical point for ${\hat{\mathcal{J}}}_h$ (up to the chosen tolerance $\tau_\text{\rm{sub}}$), i.e., it holds $$\label{First-order_critical_condition}
    \|\bar \mu-P_{\mathcal{P}}(\bar \mu-\nabla_\mu {\hat{\mathcal{J}}}_h(\bar \mu))\|_2 \leq \tau_\text{\rm{sub}}.$$

The set ${\mathcal{P}}\subset \mathbb{R}^P$ is compact. Therefore there exists a sequence of indices $\left\{k_i\right\}_{i\in\mathbb{N}}$ such that the sub-sequence $\{\mu^{(k_i)}\}_{i\in\mathbb{N}}$ converges to a point $\bar \mu\in {\mathcal{P}}$. It remains to show that $\bar \mu$ is an approximate first-order critical point. At first, note that once the RB space is enriched at a point $\mu^{(k)}$, we have $\Delta_{{\hat{\mathcal{J}}}_{{r}}^{(k)}}(\mu^{(k)})= 0$. Hence, also $q^{(k)}(\mu^{(k)}) = 0$ holds, where $$\begin{aligned}
    q^{(k)}(\mu) := \frac{\Delta_{{\hat{\mathcal{J}}}_{{r}}^{(k)}}(\mu)}{{\hat{\mathcal{J}}}_{{r}}^{(k)}(\mu)} && \text{for all } k\in\mathbb{N},\, \mu\in{\mathcal{P}}.
    \end{aligned}$$ Note that both the estimator $\Delta_{{\hat{\mathcal{J}}}_{{r}}^{(k)}}$ and $q^{(k)}$ are uniformly continuous on ${\mathcal{P}}$ for all $k\geq 0$. This follows directly from Remark \[continuity\_of\_estimator\] and the Heine-Cantor theorem. When the model is enriched at a parameter $\mu^{(k)}$, from the uniform continuity of $q^{(k)}$ it follows that for all $\varepsilon>0$ there exists an $\eta^{(k)}>0$ (depending on $\varepsilon$) such that $\|\mu^{(k)}-\mu\|_2< \eta^{(k)}$ implies $$\big|q^{(k)}(\mu)-\underset{= 0}{\underbrace{q^{(k)}(\mu^{(k)})}}\big| < \varepsilon.$$ Furthermore, due to the convergence of the sub-sequence $\left\{\mu^{(k_i)}\right\}_{i\in\mathbb{N}}\subset{\mathcal{P}}$, we have that there exists a sufficiently large constant $I>0$ and a constant $\gamma>0$ such that $\|\mu^{(k_i)}-\mu^{(k_{i+1})}\|_2< \gamma < \eta^{(k_i)}$ for all $i \geq I$. Then we have $$\begin{aligned}
    \label{unif_cont_estimator_step}
    q^{(k_i)}(\mu^{(k_{i+1})}) = \frac{\Delta_{{\hat{\mathcal{J}}}_{{r}}^{(k_i)}}(\mu^{(k_{i+1})})}{{\hat{\mathcal{J}}}_{{r}}^{(k_i)}(\mu^{(k_{i+1})})}<\varepsilon && \text{ for all } i\geq I.
    \end{aligned}$$ We want to prove that $q^{(k_{i+1}-1)}(\mu^{(k_{i+1})})< \beta_2 \delta^{(k_{i+1}-1)}$ for all $i\geq I$, such that the unique solution $\mu^{(k_{i+1})}$ to (for $k=k_{i+1}-1$) is not triggering the termination criteria . Note that $\varepsilon$ in can be chosen appropriately (which implies a certain $\eta^{(k_i)}$ for all $i\in\mathbb{N}$ and thus a sufficiently large index $I$, of course). Since the RB space is enriched at each iteration of Algorithm \[Alg:TR-RBmethod\], we especially have that $\Delta_{{\hat{\mathcal{J}}}_{{r}}^{(k_{i+1}-1)}}(\mu^{(k_{i+1})})\leq \Delta_{{\hat{\mathcal{J}}}_{{r}}^{(k_i)}}(\mu^{(k_{i+1})})$. Using , we find that $\Delta_{{\hat{\mathcal{J}}}_{{r}}^{(k_{i+1}-1)}}(\mu^{(k_{i+1})})\leq \Delta_{{\hat{\mathcal{J}}}_{{r}}^{(k_i)}}(\mu^{(k_{i+1})})< \varepsilon {\hat{\mathcal{J}}}_{{r}}^{(k_i)}(\mu^{(k_{i+1})})$. Hence, $q^{(k_{i+1}-1)}(\mu^{(k_{i+1})})< \beta_2 \delta^{(k_{i+1}-1)}$ holds for $$\varepsilon = \beta_2 \delta^{(k_{i+1}-1)}\frac{{\hat{\mathcal{J}}}_{{r}}^{(k_{i+1}-1)}(\mu^{(k_{i+1})})}{{\hat{\mathcal{J}}}_{{r}}^{(k_i)}(\mu^{(k_{i+1})})}$$ and for all $i\geq I$. This shows that from a certain iteration $I$, we are far enough from the boundary of the Trust-Region for all $i\geq I$, so that does not affect the projected BFGS algorithm. Thus, must hold for $\mu^{(k_{i+1})}= \mu^{(k_{i+1}-1,L^{(k_{i+1}-1)})}$ for $i\geq I$. Hence, we have proved that each $\mu^{(k_{i+1})}$ is an approximate first-order critical point for ${{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k_{i+1}-1)}$ (up to the chosen tolerance $\tau_\text{\rm{sub}}$) for all $i\geq I$, which yields to $$\begin{aligned}
        \|\mu^{(k_{i+1})}-P_{\mathcal{P}}(\mu^{(k_{i+1})}-\nabla_\mu{{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k_{i+1}-1)}(\mu^{(k_{i+1})})\|_2 \leq \tau_\text{\rm{sub}}, && \text{for all } i\geq I.
    \end{aligned}$$ Moreover, taking into account the RB method properties and the fact that $V_h$ is a finite dimensional space, there exists a constant $I_\nabla>0$ sufficiently large, such that $\nabla_\mu {{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k_i)}(\mu)= \nabla_\mu {\hat{\mathcal{J}}}_h(\mu)+\epsilon^{(k_i)}$ for all $\mu$ in a neighborhood of $\bar\mu$ and for $i\geq I_{\nabla}$, with $\epsilon^{(k_i)}\to 0$ as $i\to\infty$. Thus, exploiting the continuity of the projection operator and assuming $i\geq \max(I,I_\nabla)$, we have that $$\begin{aligned}
    \tau_\text{\rm{sub}} & \geq \|\mu^{(k_{i+1})}-P_{\mathcal{P}}(\mu^{(k_{i+1})}-\nabla_\mu{{{{\hat{\mathcal{J}}}_{{r}}}}}^{(k_{i+1}-1)}(\mu^{(k_{i+1})})\|_2 \\
    & = \|\mu^{(k_{i+1})}-P_{\mathcal{P}}(\mu^{(k_{i+1})}-\nabla_\mu{\hat{\mathcal{J}}}_h(\mu^{(k_{i+1})})+\epsilon^{(k_{i+1}-1)})\|_2 \to \|\bar\mu - P_{\mathcal{P}}(\bar\mu-\nabla_\mu{\hat{\mathcal{J}}}_h(\bar \mu))\|_2.
    \end{aligned}$$ Hence, the accumulation point $\bar \mu$ is an approximate first-order critical point (up to the tolerance $\tau_\text{\rm{sub}}$).

\[rmk:local\_minimum\_TR\] What remains to prove is that $\bar \mu$ is a local minimum of ${\hat{\mathcal{J}}}_h$ (or rather a sufficiently close approximation of a local minimum). Exploiting the sufficient decrease condition, one can easily show by contradiction that $\bar \mu$ is not a maximum of ${\hat{\mathcal{J}}}_h$. It can still be a saddle point as well as a local minimum. In the numerical experiments, to verify that the computed point $\bar \mu$ is actually a local minimum, we employ the second-order sufficient optimality conditions after the algorithm terminates.

Construction of RB spaces {#sec:construct_RB}
-------------------------

In an enrichment step of the outer loop of the TR-algorithm \[Alg:TR-RBmethod\] for $\mu \in {\mathcal{P}}$, we assume to have access to the primal and dual solutions $u_{h, \mu}, p_{h, \mu} \in V_h$ and consider two strategies to enrich the RB spaces.

(a) *Lagrangian RB spaces*: \[enrich:lag\] We add each FOM solution to the RB space that is directly related to its respective reduced formulation, i.e. for a given $\mu \in {\mathcal{P}}$, we enrich by $
        V^{{\textnormal{pr}},k}_{{{r}}} = V^{{\textnormal{pr}},k-1}_{{{r}}} \cup \{u_{h,\mu}\},  
        V^{{\textnormal{du}},k}_{{{r}}} = V^{{\textnormal{du}},k-1}_{{{r}}} \cup \{p_{h,\mu}\}.
        $

(b) *Single RB space*: We add all available information into a single RB space, i.e. $V^{{\textnormal{pr}}, k}_{{{r}}} = V^{{\textnormal{du}}, k}_{{{r}}} = V^{{\textnormal{pr}},k-1}_{{{r}}} \cup \{u_{h,\mu}, p_{h,\mu}\}$. According to Section \[sec:mor\] this results in ${{{{\hat{\mathcal{J}}}_{{r}}}}}(\mu)$ being equal to the standard RB reduced functional from . \[enrich:single\]

These strategies for constructing RB spaces have a significant impact on the performance and accuracy of the TR-RB method. Note that offline computations for the construction of RB models scale quadratically with the number of basis functions in the RB space. Thus, Lagrange RB spaces in are computationally beneficial compared to at a potential loss of accuracy of the corresponding RB models (since less information is added). Moreover, different spaces as in destroy the duality of state and adjoint equations, cf. Section \[sec:standard\_approach\].

Trust-Region variants based on adaptive enrichment strategies {#sec:high_order_TR}
-------------------------------------------------------------

A major contribution of this article is to introduce and analyze variants of adaptive TR-RB methods with projected BFGS as sub-problem solver for efficiently computing a solution of the optimization problem . In terms of performance we need to account for all computational costs, including traditional offline and online costs of the algorithms. The proposed methods mainly differ in terms of the model function and its gradient information. Following Section \[sec:TR\], we propose a TR method which adaptively builds an RB space along the path of optimization (see Algorithm \[Alg:TR-RBmethod\]). From a MOR perspective this diminishes the offline time of the ROM significantly since no global RB space (with respect to the parameter domain) has to be built in advance. We enrich the model after the sub-problem of the TR method has been solved. We distinguish three different approaches:

1.  *standard approach:* Following Section \[sec:standard\_approach\], the standard approach for the functional is to replace the FOM quantities by their respective ROM counterpart, i.e. we consider the map $\mu \mapsto {\hat{J}}_{{r}}(\mu)$ from . Gradient information can be computed by reducing the corresponding FOM gradient which results in ${\widetilde\nabla}_\mu {\mathcal{J}}(u_{{{r}}, \mu}, \mu)$ from . Consequently this approach does not allow for using a higher order estimate but $\Delta_{{\hat{J}}_{{r}}}(\mu)$.

2.  *semi NCD-corrected approach:* A first correction strategy is to replace the functional by the NCD-corrected RB reduced functional ${{{{\hat{\mathcal{J}}}_{{r}}}}}$ from but stick with the inexact gradient of the standard approach. This allows to use the higher order estimator for the functional, i.e. $\Delta_{{{{{\hat{\mathcal{J}}}_{{r}}}}}}(\mu)$.

3.  *NCD-corrected approach:* We propose to consider the NCD-corrected RB reduced functional ${{{{\hat{\mathcal{J}}}_{{r}}}}}$ from and its actual gradient according to Propostion \[prop:true\_corrected\_reduced\_gradient\_adj\]. Note that we only need to solve two additional equations, independently of the dimension $P$ of the parameter space.

For the basis construction, we may use variants or from Section \[sec:construct\_RB\]. Note however, that by using the basis enrichment , all approaches $1.$ - $3.$ are equivalent. Using variant with BFGS is inspired from [@QGVW2017]. However, our algorithms differ from the TR-RB approach in [@QGVW2017] since we are working with the NCD-corrected reduced cost functional (in 2) and its actual gradient (in 2 and 3). Note that the presence of inequality constraints, which are missing in [@QGVW2017], implies a projection-based optimization algorithm. In addition, we stress that, differently from [@QGVW2017], we take advantage of the proposed condition for enlarging the TR radius and of a stopping criterium independent from the RB a posteriori estimates, as presented in Section \[sec:TR\].

Note that we do not use the sensitivity based quantities from Section [[\[sec:sens\_estimation\]]{}]{} although they suggest the highest numerical accuracy w.r.t. the FOM optimality system. However, for the experiments in Section \[sec:num\_experiments\], additional computational cost for computing FOM sensitivities will not pay off in the TR-RB algorithm, especially for high-dimensional parameter spaces.

Numerical experiments {#sec:num_experiments}
=====================

We present numerical experiments to demonstrate the adaptive TR-RB variants from Section \[sec:high\_order\_TR\] with both RB constructions from Section \[sec:construct\_RB\] for quadratic objective functionals with elliptic PDE constraints as in , and compare them to state-of-the art algorithms from the literature. We also validate the higher-order a posteriori error estimates from Section \[sec:a\_post\_error\_estimates\] numerically. We consider two setups: first, the elliptic thermal fin problem from [@QGVW2017 Sec. 5.1.1] (where the correction term of the proposed NCD-corrected approach vanishes) in Section \[sec:thermalfin\]. Second, we consider a more challenging optimization problem in Section \[sec:mmexc\_example\], including a detailed analysis of the a posteriori error estimates from Section \[sec:a\_post\_error\_estimates\]. All simulations have been performed with a pure `Python` implementation based on the open source MOR library pyMOR [@milk2016pymor], making use of pyMORs builtin vectorized `numpy/scipy`-based discretizer for the FOM and generic MOR algorithms for projection and orthonormalization (such as a stabilized Gram-Schmidt algorithm) to effortlessly obtain efficient ROMs. The source code to reproduce all results (including detailed interactive `jupyter`-notebooks[^4]) is available at [@Code]. All experiments are based on the same implementation (including a reimplementation of [@QGVW2017]) and were performed on the same machine multiple times to avoid caching or multi-query effects. Timings may thus be used to compare and judge the computational efficiency of the different algorithms.

We consider stationary heat transfer in a bounded connected spatial domain $\Omega \subset {\mathbb{R}}^2$ with polygonal boundary $\partial\Omega$ partitioned into a non-empty Robin boundary $\Gamma_\text{\rm{R}} \subset \partial\Omega$ and possibly empty distinct Neumann boundary $\Gamma_\text{\rm{N}} = \partial\Omega\backslash \Gamma_\text{\rm{R}}$, and unit outer normal $n: \partial\Omega \to {\mathbb{R}}^2$. We consider the Hilbert space $V = H^1(\Omega) := \{v \in L^2(\Omega) \,|\, \nabla v \in L^2(\Omega) \}$ of weakly differentiable functions and, for an admissible parameter $\mu \in {\mathcal{P}}$, we seek the temperature $u_\mu \in V$ as the solution of $$\begin{aligned}
  -\nabla\cdot\big(\kappa_\mu \nabla u_\mu \big) &= f_\mu \ \text{in }\Omega, &
  \kappa_\mu \nabla u_\mu\cdot n &= c_\mu (u_\text{\rm{out}} - u_\mu) \ \text{on }\Gamma_\text{\rm{R}},
\label{eq:heat_equation}
  & \kappa_\mu\nabla u_\mu \cdot n &= g_\text{\rm{N}} \ \text{on }\Gamma_\text{\rm{N}}\end{aligned}$$ in the weak sense, with the admissible parameter set, the spatial domain and its boundaries and the data functions $\kappa_\mu \in L^\infty(\Omega)$, $f_\mu \in L^2(\Omega)$, $c_\mu \in L^\infty(\Gamma_\text{\rm{R}})$ and $u_\text{\rm{out}} \in L^2(\Gamma_\text{\rm{R}})$ defined in the respective experiment. The bilinear form $a$ and linear functional $l$ in are thus given for all $\mu \in {\mathcal{P}}$ and $v, w \in V$ by $$\begin{aligned}
  a_\mu(v, w) := \int_\Omega\kappa_\mu \nabla v \cdot \nabla w {{\operatorname{d}{\hspace{-1.25pt}} x}}+ \int_{\Gamma_\text{\rm{R}}}\hspace{-5pt}c_\mu\, vw{{\operatorname{d}{\hspace{-1.25pt}} s}}&&\text{and}&& l_\mu(v) := \int_\Omega f_\mu\,v{{\operatorname{d}{\hspace{-1.25pt}} x}}+ \int_{\Gamma_\text{\rm{R}}}\hspace{-5pt}c_\mu\, u_\text{\rm{out}}v{{\operatorname{d}{\hspace{-1.25pt}} s}}+ \int_{\Gamma_\text{\rm{N}}}\hspace{-5pt}g_\text{\rm{N}}v{{\operatorname{d}{\hspace{-1.25pt}} s}}.\end{aligned}$$ For the FOM we fix a fine enough reference simplicial or cubic mesh and define $V_h \subset V$ as the respective space of continuous piecewise (bi-)linear Finite Elements.

Since the inner product and norm have a big influence on the computational efficiency of the a posteriori error estimates as well as their sharpness, we use the mesh-independent energy-product $(u, v) := a_{\check{\mu}}(u, v)$ for a fixed parameter $\check{\mu} \in {\mathcal{P}}$, which is a product over $V$ due to the symmetry, continuity and coercivity of the bilinear form for each example below. Owing to this choice of the product, we may use the $\min$-theta approach from [@HAA2017 Prop. 2.35] to obtain lower bounds on coercivity constants and the $\max$-theta approach from [@HAA2017 Ex. 5.12] to obtain upper bounds on continuity constants, each required for the a posteriori error estimates. Compared to the more general Successive Constraint Method [@pat07], this approach yields quite sharp estimates and is computationally more efficient, both offline and online. Due to Assumption \[asmpt:parameter\_separable\] and the bi-linearity of the objective functional, we may carry out the preassembly of all high-dimensional quantities after each enrichment, which is well-known for RB methods [@HAA2017 Sec. 2.5]. We would like to point out that while the more accurate and stable preassembly of the estimates from [@buhr14] is readily available in pyMOR, the slightly cheaper standard preassembly of the estimates was sufficient for our experiments.

For all experiments, we use an initial TR radius of $\delta^0 = 0.1$, a TR shrinking factor $\beta_1=0.5$, an Armijo step-length $\kappa=0.5$, a truncation of the TR boundary of $\beta_2 = 0.95$, a tolerance for enlarging the TR radius of $\eta_\varrho = 0.75$, a stopping tolerance for the TR sub-problems of $\tau_\text{\rm{sub}} = 10^{-8}$, a maximum number of TR iteration $K = 40$, a maximum number of sub-problem iteration $K_{\text{\rm{sub}}}= 400$ and a maximum number of Armijo iteration of $50$. We also point out that the stopping tolerance for the FOC condition $\tau_\text{\rm{FOC}}$ is specified in each experiment.

State of the art optimization methods {#sec:state_of_the_art_methods}
-------------------------------------

We compare our proposed methods to the following ones from the literature:

**Adaptive TR-RB with BFGS sub-problem solver and Lagrangian RBs [@QGVW2017]**: We consider the same method as in [@QGVW2017], where the authors used the standard functional and gradient from Section \[sec:standard\_approach\]. Furthermore, no enlarging strategy has been used for the TR-radius and no projection for parameter constraints has been considered. Importantly, the authors did not take advantage of the fact, that the full order FOC condition in line 23 of Algorithm \[Alg:TR-RBmethod\] is cheaply available after an enrichment step. Instead they used the reduced FOC condition plus the estimator for the gradient of the cost functional $\|{\widetilde\nabla}_\mu {\hat{J}}_{{r}}(\mu^{(k+1)}))\|_2 + \Delta_{{\widetilde\nabla}_\mu {\hat{J}}_{{r}}}(\mu) \leq \tau_\text{\rm{FOC}}$ in line 23. Note that this approach has multiple drawbacks. First, the evaluation is more costly due to the estimator. Second, it is less accurate and third, it can prevent the TR-RB from converging in case the estimator is not able to be small enough (for instance governed by high constants or numerical issues in the estimator).

**FOM projected BFGS**: We consider a standard projected BFGS method, which uses FOM evaluations of the forward model to compute the reduced cost functional and its gradient. We restrict the maximum number of iterations by $400$.

Model problem 1: Elliptic thermal fin model problem {#sec:thermalfin}
---------------------------------------------------

We consider the six-dimensional elliptic thermal fin example from [@QGVW2017 Sec. 5.1.1] and refer to Figure \[fig:fin\] for the problem definition. The purpose of this experiment is to show the applicability of the proposed algorithms and to compare them to the one proposed in [@QGVW2017]. For all runs we prescribe the same desired parameter $\mu^\text{d} \in {\mathcal{P}}$ by randomly drawing $k_1, \dots, k_4$ strictly within ${\mathcal{P}}$ and by setting $k_0 = 0.1$ and $\text{Bi} = 0.01$, to artificially mimic the situation where parameter constraints have to be tackled. Defining $T^\text{d} := q(u_{\mu^\text{d}})$ where $u_{\mu^\text{d}} \in V$ is the solution of associated with the desired parameter and where $q(v) := \int_{\Gamma_\text{\rm{N}}} v {{\operatorname{d}{\hspace{-1.25pt}} s}}$ for $v\in V$ denotes the mean temperature at the root of the fin, we consider a cost functional ${\mathcal{J}}(u, \mu) = \Theta(\mu) + j_\mu(u) + k_\mu(u, u)$ as in with $\Theta(\mu) := (\|\mu^\text{d} - \mu\| / \|\mu^\text{d}\|)^2 + {T^\text{d}}^2 + 1$, $j_\mu(v) := - T^\text{d}\,q(v)$ and $k_\mu(v, w) := 1/2\,q(v)\,q(w)$. We would like to point out that the authors in [@QGVW2017] dropped the ${T^\text{d}}^2 + 1$ term from the definition of $\Theta$, which we re-add to ensure Assumption \[asmpt:bound\_J\]. This constant term does not change the position of local minima and the derivatives of the cost functional. However, this makes the Trust-Region radius shrink especially at the beginning, slowing down the TR-RB methods. This does not affect the comparison among the TR-RB methods, since all suffer from this issue. Note that for this particular example, the proposed NCD-correction term vanishes, see Remark \[rem:fin\]. For the FOM, we generate an unstructured simplicial mesh using pyMORs `gmsh` (see [@gmsh]) bindings, resulting in $\dim V_h = 77537$.

(-.5,0) – (-.5,.75); (-.5,.75) – (-3.,.75); at (-1.8,0.86) [[$k_1$]{}]{}; at (-1.8,1.86) [[$k_2$]{}]{}; at (-1.8,2.86) [[$k_3$]{}]{}; at (-1.8,3.86) [[$k_4$]{}]{}; at (1.8,0.86) [[$k_1$]{}]{}; at (1.8,1.86) [[$k_2$]{}]{}; at (1.8,2.86) [[$k_3$]{}]{}; at (1.8,3.86) [[$k_4$]{}]{}; at (0,2) [[$k_0$]{}]{}; (-3.,.75) – (-3.,1) – (-.5,1); (-.5,1) – (-.5,1.75) – (-3.,1.75) – (-3.,2) – (-.5,2); (-.5,2) – (-.5,2.75) – (-3.,2.75) – (-3.,3) – (-.5,3); (-.5,3) – (-.5,3.75) – (-3.,3.75) – (-3.,4) – (3,4); (.5,0) – (.5,.75) – (3.,.75) – (3.,1) – (.5,1); (.5,1) – (.5,1.75) – (3.,1.75) – (3.,2) – (.5,2); (.5,2) – (.5,2.75) – (3.,2.75) – (3.,3) – (.5,3); (.5,3) – (.5,3.75) – (3.,3.75) – (3.,4); (-.5,0) – (.5,0) node \[midway,above\] [$\Gamma_\text{\rm{N}}$]{}; (-.5,.75) – (-.5,1); (-.5,1.75) – (-.5,2); (-.5,2.75) – (-.5,3); (-.5,3.75) – (-.5,4); (.5,.75) – (.5,1); (.5,1.75) – (.5,2); (.5,2.75) – (.5,3); (.5,3.75) – (.5,4); (3.2,0.45) – (3.2,0.75); (3.2,1.3) – (3.2,1.); (0.5,0.45) – (3,0.45) node\[midway, fill=white\] [$L$]{}; (3,0.75) – (3,1) node\[midway,right\] [$t$]{};

Starting with ten different randomly drawn initial parameters $\mu^{(0)}$, we measure the total computational runtime, the number of TR iterations $k$ and the error in the optimal parameter for all combinations of adaptive TR algorithms from Section \[sec:TRRB\_and\_adaptiveenrichment\] and choice of RB spaces from Section \[sec:construct\_RB\], as well as for the state of the art methods from the literature from Section \[sec:state\_of\_the\_art\_methods\].

                                                                               av. (min/max) runtime\[s\]   speed-up   av. (min/max) iter.       rel. error             FOC cond.
  --------------------------------------------------------------------------- ---------------------------- ---------- --------------------- --------------------- ---------------------
  FOM proj. BFGS                                                                967.86 (176.69/3401.06)        –         111.20 (25/400)     $3.13\cdot 10^{-3}$   $1.19\cdot 10^{-2}$
  TR-RB from [@QGVW2017]                                                          68.06 (43.28/88.21)        10.40         7.20 (8/13)       $1.34\cdot 10^{-6}$   $4.31\cdot 10^{-5}$
  1(a) TR-RB with $V_{{r}}^{\textnormal{pr}}\neq V_{{r}}^{\textnormal{du}}$       44.56 (34.22/74.96)        21.72         8.80 (8/11)       $3.08\cdot 10^{-6}$   $4.64\cdot 10^{-5}$
  1(b) TR-RB with $V_{{r}}^{\textnormal{pr}}= V_{{r}}^{\textnormal{du}}$          43.86 (34.09/74.35)        22.07         8.70 (8/10)       $3.37\cdot 10^{-6}$   $6.40\cdot 10^{-5}$

  :  [Performance and accuracy of selected algorithms for the example from Section \[sec:thermalfin\] for ten optimization runs with randomly initial guesses $\mu^{(0)}$: averaged, minimum and maximum total computational time (column 2) and speed-up compared to the FOM variant (column 3); average, minimum and maximum number of iterations $k$ required until convergence (column 4), average relative error in the parameter (column 5) and average FOC condition (column 6).]{} []{data-label="table:fin"}

All considered optimization methods converged (up to a tolerance), but we restrict the presentation to the most informative ones (all results can be found in the accompanying code). As we observe from Table \[table:fin\], the ROM based adaptive TR-RB algorithms vastly outperform the FOM variant, noting that the computational time of the ROM variants includes all offline and online computations. Figure \[fig:fin\_timing\_mu\_d\] details the decay of the error decay in the optimal parameter during the optimization for a selected random initial guess. We observe that the choice of the RB enrichment does not impact the performance of the algorithm for this example too much, see Remark \[rem:fin\]. Also methods 2(\[enrich:lag\]) and 3(\[enrich:lag\]) show a comparable computational speed (not shown). We also observe that the method from [@QGVW2017] requires more time and more iterations on average, variants 1 are still faster due to the enlarging of the TR radius and of the use of a termination criterium which does not depend on a posteriori estimates, which may force additional TR iterations.

\[rem:fin\] It is important to notice that this model problem is not suitable to fully demonstrate the capabilities of the NCD-corrected approach. The reason is that the choice of the functional is a misfit on only the root edge of the thermal fin, plus a Tikhonov regularization term. Since the root of the thermal fin is also the source of the primal problem, the dual solutions $p_{{{r}}, \mu}$ of the reduced dual equation are thus linearly dependent on the respective primal solutions $u_{{{r}}, \mu}$ and the correction term $r_\mu^{\textnormal{pr}}(u_{{{r}},\mu})[p_{{{r}},\mu}]$ for the NCD-corrected RB reduced functional from vanishes. In general, for quadratic objective functionals, this is not the case and all variants with correction terms thus waste unnecessary computational time.

Model problem 2: stationary heat distribution in a building {#sec:mmexc_example}
-----------------------------------------------------------

For these experiments we consider as objective functional a weighted $L^2$-misfit on a domain of interest $D \subseteq \Omega$ and a weighted Tikhonov term comparable to design optimization, optimal control or inverse problems, i.e. $$\begin{aligned}
\mathcal{J}(v, \mu) = \frac{\sigma_D}{2} \int_{D}^{} (v - u^{\text{d}})^2 + \frac{1}{2} \sum^{M}_{i=1} \sigma_i (\mu_i-\mu^{\text{d}}_i)^2 + 1,
\label{eq:mmexc_example_J}\end{aligned}$$ with given desired state $u^{\text{d}} \in V$ and parameter $\mu^{\text{d}} \in {\mathcal{P}}$ and weights $\sigma_D, \sigma_i$ specified further below. With respect to , we thus have $\Theta(\mu) = \frac{1}{2} \sum^{M}_{i=1} \sigma_i (\mu_i-\mu^{\text{d}}_i)^2 + \frac{\sigma_D}{2} \int_{D}^{} u^{\text{d}} u^{\text{d}} + 1$, $j_{\mu}(u) = -\sigma_D \int_{D}^{} u^{\text{d}}u$ and $k_{\mu}(u,v) = \frac{\sigma_D}{2} \int_{D}^{} uv$.

[0.60]{} ![image](Pictures/full_diffusion_with_big_numbers_with_D.png){width="\textwidth"}

Motivated by ensuring a desired temperature in a single room of a building floor, we consider blueprints with windows, heaters, doors and walls, yielding parameterized diffusion, forces and boundary values as sketched in Figure \[ex1:blueprint\].[^5] For simplicity we omit a realistic modeling of temperature and restrict ourselves to academic numbers of the diffusion and heat source quantities. We seek to ensure a desired temperature $u^\text{d}=18$ and set $\mu^\text{d}_i = 0$. For the FOM discretization we choose a cubic mesh which resolves all features of the data functions derived from Figure \[ex1:blueprint\], resulting in $\dim V_h = 80601$ degrees of freedom. We consider a ten-dimensional parameter example with three wall sets $\{1|,2|,3|,8|\}$, $\{4|,5|,6|,7|\}$ and $\{9|\}$ and seven heater sets, $\{1,2\}$, $\{3,4\}$ and $\{5\}$, $\{6\}$, $\{7\}$, $\{8\}$ and $\{9,10,11,12\}$ (each set governed by a single parameter component). The set of admissible parameters is given by ${\mathcal{P}}= [0.025,0.1]^3\times[0,100]^7$ and we choose $\sigma_D= 100$ and $(\sigma_i)_{1\leq i\leq 10} = (10\sigma_w,5\sigma_w,\sigma_w,2\sigma_h,$ $2\sigma_h,\sigma_h,\sigma_h,\sigma_h,\sigma_h,4\sigma_h)$ in , with $\sigma_w = 0.05$ and $\sigma_h= 0.001$. The choice of $\sigma_i$ is related to the measure of the walls and how many heaters are considered in each group. The other components of the data functions are fixed and thus not directly involved in the optimization process. Briefly, the diffusion coefficient of air and inside doors is set to $0.5$, of the outside wall to $0.001$, of outside doors $\underline{8}$ and $\underline{9}$ to $0.01$ and of windows to $0.025$. For the energy product, we choose $\check{\mu}= (0.05,0.05,0.05,10,10,10,10,10,10,10)$.

We use this setup to inspect different TR-RB algorithms in Section \[sec:mmexc\_opt\_results\], but also to study the a posteriori error estimates from Section \[sec:a\_post\_error\_estimates\] in the following section.

### Numerical validation of the a posteriori error estimates {#sec:estimator_study}

To study the performance of the a posteriori error estimates proposed in Section \[sec:a\_post\_error\_estimates\], we neglect the outer-loop optimization and simply use a goal oriented adaptive greedy algorithm [@HDO2011] with basis extension (a) from Section \[sec:construct\_RB\] to generate a ROM which ensures that the worst relative estimated error for the reduced functional and its gradient over the adaptively generated training set and a randomly chosen validation set is below a prescribed tolerance of $\tau_{\text{\rm{FOC}}}= 5 \cdot 10^{-4}$. In particular we first ensure $\Delta_{\hat{J}_{{r}}}(\mu)/{\hat{J}}_{{r}}(\mu) < \tau_{\rm{FOC}}$ for $\Delta_{\hat{J}_{{r}}}$ from Proposition \[prop:Jhat\_error\].i and continue with $\Delta_{{\widetilde\nabla}{\hat{J}}_{{r}}}(\mu)/\|{\widetilde\nabla}{\hat{J}}_{{r}}(\mu)\|_2 < \tau_{{\widetilde\nabla}{\hat{J}}}$ for $\Delta_{{\widetilde\nabla}{\hat{J}}_{{r}}}$ from Proposition \[prop:grad\_Jhat\_error\].i, cf. [@QGVW2017 Algorithm 2]. Let us mention that the goal for $\Delta_{\hat{J}_{{r}}}$ is fulfilled after $24$ basis enrichments and we have $\Delta_{{\widetilde\nabla}{\hat{J}}_{{r}}}(\mu)/\|{\widetilde\nabla}{\hat{J}}_{{r}}(\mu)\| < 4.84$ after $56$ basis enrichments, where we artificially stop the algorithm since the associated computational effort is already roughly 17 hours, demonstrating the need for the proposed adaptive TR-RB algorithm studied in the next section.

As we observe from Figure \[fig:estimator\_study\], the error of the NCD-corrected terms is of several orders of magnitude smaller than the corresponding terms of the standard approach. It can also be seen that the (computationally more costly) sensitivity bases quantities, i.e. ${\widetilde\nabla}{{{{\hat{\mathcal{J}}}_{{r}}}}}$, show the best error. However, all estimators for the corrected and sensitivity based quantities show a worse effectivity, hinting that there is still room for improvement.

### Optimization results {#sec:mmexc_opt_results}

Similar to Section \[sec:thermalfin\], starting with ten different randomly drawn initial parameters $\mu^{(0)}$, we measure the total computational runtime, the number of TR iterations $k$ and the error in the optimal parameter for all combinations of adaptive TR algorithms from Section \[sec:TRRB\_and\_adaptiveenrichment\] and choice of RB spaces from Section \[sec:construct\_RB\], as well as for the state of the art methods from the literature from Section \[sec:state\_of\_the\_art\_methods\].

All algorithms converged (up to a tolerance) to the same point $\bar\mu$ and it was verified a posteriori that this point is a local minimum of ${\hat{\mathcal{J}}}$, i.e. it satisfies the second-order sufficient optimality conditions. The value of $\bar\mu$ in order to compute the relative error was calculated with the FOM projected Newton method for a FOC condition tolerance of $10^{-12}$ and, thanks to the choice of the cost functional weights, the target $u^\text{d}$ is approximate by $\bar u$ with a relative error of $1.7\cdot 10^{-6}$ in $D$. We consider the same setup for two different stopping tolerances $\tau_{\text{\rm{FOC}}}= 5\cdot 10^{-4}$ and $\tau_{\text{\rm{FOC}}}= 10^{-6}$ to demonstrate that the performance (both in terms of time and convergence) of the methods vastly depends on the choice of $\tau_{\text{\rm{FOC}}}$.

  <span style="font-variant:small-caps;">**(A) Result for $\tau_{\text{\rm{FOC}}}= 5\cdot 10^{-4}$**</span>    av. (min/max) runtime\[s\]   speed-up   av. (min/max) iter.       rel. error             FOC cond.
  ----------------------------------------------------------------------------------------------------------- ---------------------------- ---------- --------------------- --------------------- ---------------------
  FOM proj. BFGS                                                                                                 332.57 (196.51/591.85)        –          44.30 (30/60)      $1.40\cdot 10^{-3}$   $1.80\cdot 10^{-4}$
  TR-RB from [@QGVW2017]                                                                                         117.87 (70.29/166.31)        2.82        10.10 (6/14)       $5.46\cdot 10^{-4}$   $1.41\cdot 10^{-4}$
  1(a) TR-RB with $V_{{r}}^{\textnormal{pr}}\neq V_{{r}}^{\textnormal{du}}$                                       91.50 (47.07/230.29)        3.63         8.30 (5/10)       $2.01\cdot 10^{-3}$   $2.04\cdot 10^{-4}$
  1(b) TR-RB with $V_{{r}}^{\textnormal{pr}}= V_{{r}}^{\textnormal{du}}$                                          78.65 (54.69/114.36)        4.23         6.90 (5/9)        $2.53\cdot 10^{-4}$   $8.23\cdot 10^{-5}$
  2(a) TR-RB semi NCD-corrected                                                                                   79.47 (63.38/94.28)         4.18         8.50 (7/10)       $5.98\cdot 10^{-5}$   $1.02\cdot 10^{-4}$
  3(a) TR-RB NCD-corrected                                                                                        71.84 (50.38/87.16)         4.63         7.40 (5/9)        $1.09\cdot 10^{-3}$   $6.12\cdot 10^{-5}$

  :  [Performance and accuracy of selected algorithms for two choices of $\tau_\text{\rm{FOC}}$ for the example from Sec. \[sec:mmexc\_opt\_results\] for ten optimization runs with random initial guess, compare Table \[table:fin\].]{} []{data-label="table:MP2_com"}

  <span style="font-variant:small-caps;">**(B) Result for $\tau_{\text{\rm{FOC}}}= 10^{-6}$**</span>    av. (min/max) runtime\[s\]   speed-up   av. (min/max) iter.       rel. error             FOC cond.
  ---------------------------------------------------------------------------------------------------- ---------------------------- ---------- --------------------- --------------------- ---------------------
  FOM proj. BFGS                                                                                          409.28 (317.25/637.55)        –          57.00 (49/71)      $2.82\cdot 10^{-6}$   $3.35\cdot 10^{-7}$
  TR-RB from [@QGVW2017]                                                                                  614.81 (566.66/671.97)       0.66        40.00 (40/40)      $8.46\cdot 10^{-7}$   $8.44\cdot 10^{-8}$
  1(a) TR-RB with $V_{{r}}^{\textnormal{pr}}\neq V_{{r}}^{\textnormal{du}}$                               165.48 (92.26/417.24)        2.47        15.30 (10/40)      $3.29\cdot 10^{-6}$   $5.43\cdot 10^{-7}$
  1(b) TR-RB with $V_{{r}}^{\textnormal{pr}}= V_{{r}}^{\textnormal{du}}$                                   86.39 (62.68/124.43)        4.74         7.80 (6/10)       $3.52\cdot 10^{-6}$   $3.03\cdot 10^{-7}$
  2(a) TR-RB semi NCD-corrected                                                                            90.37 (80.97/102.60)        4.53         9.80 (9/11)       $8.12\cdot 10^{-7}$   $2.26\cdot 10^{-7}$
  3(a) TR-RB NCD-corrected                                                                                 88.24 (58.18/108.90)        4.64         8.90 (6/10)       $2.65\cdot 10^{-6}$   $2.73\cdot 10^{-7}$

  :  [Performance and accuracy of selected algorithms for two choices of $\tau_\text{\rm{FOC}}$ for the example from Sec. \[sec:mmexc\_opt\_results\] for ten optimization runs with random initial guess, compare Table \[table:fin\].]{} []{data-label="table:MP2_com"}

From Table \[table:MP2\_com\], we observe that all proposed TR-RB methods speed up the FOM projected BFGS method with the NCD-corrected approach outperforming the others, since the gradient used is the true one of the model function ${\hat{\mathcal{J}}}_{{r}}$. Moreover, independently of the model function, the algorithm from [@QGVW2017] is much slower, demonstrating the positive impact of the suggested improvements on enlarging the TR radius and on the termination criterium based on cheaply available FOM information (instead of relying on an a posteriori estimate), also visible in the number of outer TR iterations. Comparing our proposed TR variants in terms of iterations, it is more beneficial to consider a single RB space, i.e. $V^{\textnormal{pr}}_{{r}}= V^{\textnormal{du}}_{{r}}$. While enrichment (a) is more costly and the time-to-ROM-solution is slightly larger, the richer space seems to allow for better approximations of ${\hat{\mathcal{J}}}_h$.

All methods approximate the optimal parameter $\bar\mu$ with a small relative error and reach the desired tolerance for the FOC condition. However, in view of the resulting relative error in Table \[table:MP2\_com\] and Figure \[Fig:\_FOC\_cond\], we observe that the choice $\tau_\text{\rm{FOC}}= 5\cdot10^{-4}$ is not sufficiently small for this model problem. In fact, we observe for most of the variants, that this choice for the tolerance $\tau_\text{\rm{FOC}}$ does not guarantee an adequately low relative error in approximating $\bar\mu$ and affects the timings by stopping the method too early. We conclude that the choice $\tau_{\text{\rm{FOC}}}= 10^{-6}$ instead results in a valid optimum of all variants (up to a tolerance of $10^{-6}$). Importantly, for this choice of $\tau_\text{\rm{FOC}}$, we point out that the variant from [@QGVW2017] only stopped because we restricted the maximum number of iterations to $40$, although the FOC condition dropped under the depicted tolerance of $10^{-6}$.

[0.45]{}

[0.45]{}

This is caused by the fact that in [@QGVW2017] the a posteriori estimate, which is summed to the FOC condition, cannot get numerically small enough, showing the limit of the proposed stopping criterium in [@QGVW2017]. From Figure \[Fig:\_FOC\_cond\](B) we conclude that the NCD-corrected approaches 2(\[enrich:lag\]) and 3(\[enrich:lag\]) outperform the standard ROM variant 1(\[enrich:lag\]), which also reached the maximum number of iterations for one of the ten samples. Consequently, the NCD-correction entirely resolves the issue of the variational crime (introduced by splitting the reduced spaces), since it shows roughly the same performance as variant 1(\[enrich:single\]). However, looking at the minimum and maximum number of computational time in Table \[table:MP2\_com\], variant 3(\[enrich:lag\]) shows a less volatile and more robust behavior.

Conclusion
==========

In this work we proposed and analyzed several variants of new adaptive Trust-Region Reduced Basis methods for parameterized partial differential equations. First, we proved convergence of the modified algorithm in case of additional bilateral constraints on the parameter set, making this method more appealing for real-world applications. Second, the use of a NCD-corrected RB reduced functional improves the RB approximation compared to the standard approach, and enables the possibility of using an exact gradient in the case of separate RB spaces (each variant accompanied by rigorous a posteriori error estimates). This approach turns out to be the most reliable in terms of computational time and accuracy, outperforming the existing TR-RB method. Furthermore, the proposed cheaply-computable criteria for enlarging the TR radius and for terminating the iterations ensure a faster convergence. In future works we are interested in considering the projected Newton method to replace the projected BFGS method used in this contribution. This leads to additional effort on developing a posteriori estimates for the RB approximation of the hessian and of the optimal parameter. In addition, we are interested in combining the proposed TR-RB algorithm with localized RB methods for large-scale applications.

[^1]: Mathematics Münster, Westfälische Wilhelms-Universität Münster, Einsteinstr. 62, D-48149 Münster, [{tim.keil,mario.ohlberger,felix.schindler,}@uni-muenster.de]({tim.keil,mario.ohlberger,felix.schindler,}@uni-muenster.de)

[^2]: Department of Mathematics and Statistics, University of Konstanz, D-78457 Konstanz, [{luca.mechelli,stefan.volkwein}@uni-konstanz.de]({luca.mechelli,stefan.volkwein}@uni-konstanz.de)

[^3]: [**Funding:**]{} The authors acknowledge funding by the Deutsche Forschungsgemeinschaft for the project [*Localized Reduced Basis Methods for PDE-constrained Parameter Optimization*]{} under contracts OH 98/11-1; SCHI 1493/1-1; VO 1658/6-1. Tim Keil, Mario Ohlberger and Felix Schindler acknowledge funding by the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy EXC 2044 390685587, Mathematics Münster: Dynamics – Geometry – Structure.

[^4]: Available at <https://github.com/TiKeil/NCD-corrected-TR-RB-approach-for-pde-opt>.

[^5]: See <https://github.com/TiKeil/NCD-corrected-TR-RB-approach-for-pde-opt> for the definition of the data functions.
