# Model Reduction for Large Scale Systems

Tim Keil  and  Mario Ohlberger  Mathematics Münster, Westfälische Wilhelms-
Universität Münster, Einsteinstr. 62, D-48149 Münster, Germany.
https://www.wwu.de/AMM/ohlberger<EMAIL_ADDRESS>
###### Abstract.

Projection based model order reduction has become a mature technique for
simulation of large classes of parameterized systems. However, several
challenges remain for problems where the solution manifold of the
parameterized system cannot be well approximated by linear subspaces. While
the online efficiency of these model reduction methods is very convincing for
problems with a rapid decay of the Kolmogorov n-width, there are still major
drawbacks and limitations. Most importantly, the construction of the reduced
system in the offline phase is extremely CPU-time and memory consuming for
large scale and multi scale systems. For practical applications, it is thus
necessary to derive model reduction techniques that do not rely on a classical
offline/online splitting but allow for more flexibility in the usage of
computational resources. A promising approach with this respect is model
reduction with adaptive enrichment. In this contribution we investigate
Petrov-Galerkin based model reduction with adaptive basis enrichment within a
Trust Region approach for the solution of multi scale and large scale PDE
constrained parameter optimization.

###### Key words and phrases:

PDE constraint optimization and reduced basis method and trust region method.

The authors acknowledge funding by the Deutsche Forschungsgemeinschaft for the
project Localized Reduced Basis Methods for PDE-constrained Parameter
Optimization under contract OH 98/11-1 and by the Deutsche
Forschungsgemeinschaft under Germany’s Excellence Strategy EXC 2044 390685587,
Mathematics Münster: Dynamics – Geometry – Structure.

##  1\. Introduction

Model order reduction (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 on recent
development we refer e.g. to [4]. 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 spaces that are spanned from
suitably selected particular solutions, called snapshots. For time-dependent
problems, the POD-Greedy method [9] defines the Gold-Standard. As RB methods
rely on so called efficient offline/online splitting, they need to be combined
with supplementary interpolation methods in case of non-affine parameter
dependence or non-linear differential equations. The empirical interpolation
method (EIM) [3] and its various generalizations, e.g. [6], are key
technologies with this respect. While RB methods are meanwhile very well
established and analyzed for scalar coercive problems, there are still major
challenges for problems with a slow convergence of the Kolmogorov N-width
[13]. Such problems in particular include pPDEs with high dimensional or even
infinite dimensional parameter dependence, multiscale problems as well as
hyperbolic or advection dominated transport problems. Particular promising
approaches for high dimensional parameter dependence and large or multiscale
problems are localized model reduction approaches. We refer to [5] for a
recent review of such approaches, including the localized reduced basis
multiscale method (LRBMS) [15]. Several of these approaches have already been
applied in multiscale applications, in particular for battery simulation with
resolved electrode geometry and Buttler-Volmer kinetics [7]. Based on
efficient localized a posteriori error control and online enrichment, these
methods overcome traditional offline/online splitting and are thus
particularly well suited for applications in optimization or inverse problems
as recently demonstrated in [16, 14]. In the context of PDE constrained
optimization, a promising Trust Region (TR) – RB approach that updates the
reduced model during the trust region iteration has recently been studied in
[17, 11, 2]. In the latter two contributions a new non-conforming dual (NCD)
approach has been introduced that improves the convergence of the adaptive TR-
RB algorithm in the case of different reduced spaces for the corresponding
primal and dual equations of the first order optimality system.

While these contributions were all based on Galerkin projection of the
respective equations, we will introduce a new approach based on Petrov-
Galerkin (PG) projection in the following. As we will demonstrate in Section 2
below, the PG-reduced optimality system is a conforming approximation which
allows for more straight forward computation of derivative information and
respective a posteriori error estimates. In Section 3 we evaluate and compare
the Galerkin and Petrov-Galerkin approaches with respect to the error behavior
and the resulting TR-RB approaches with adaptive enrichment. Although the
convergence of the TR-RB method can be observed for both approaches, the
results demonstrate that the PG approach may require further investigation
with respect to stabilization.

##  2\. Petrov-Galerkin based model reduction for PDE constrained optimization

In this contribution we consider the following class of PDE constrained
minimization problems:

(P) (P.a) $\displaystyle\min_{\mu\in\mathcal{P}}\mathcal{J}(u_{\mu},\mu),$
$\displaystyle\text{with
}\mathcal{J}(u,\mu)=\Theta(\mu)+j_{\mu}(u)+k_{\mu}(u,u),$ subject to
$u_{\mu}\in V$ being the solution of the _state – or primal – equation_ (P.b)
$\displaystyle a_{\mu}(u_{\mu},v)=l_{\mu}(v)$ $\displaystyle\text{for all
}v\in V,$

where $\Theta\in\mathcal{P}\to\mathbb{R}$ denotes a parameter functional.

Here, $V$ denotes a real-valued Hilbert space with inner product
$(\cdot\,,\cdot)$ and its induced norm $\|\cdot\|$,
$\mathcal{P}\subset\mathbb{R}^{P}$, with $P\in\mathbb{N}$ denotes a compact
and convex admissible parameter set and
$\mathcal{J}:V\times\mathcal{P}\to\mathbb{R}$ a quadratic continuous
functional. In particular, we consider box-constraints of the form

$\mathcal{P}:=\left\\{\mu\in\mathbb{R}^{P}\,|\,\mu_{\mathsf{a}}\leq\mu\leq\mu_{\mathsf{b}}\right\\}\subset\mathbb{R}^{P},$

for given parameter bounds
$\mu_{\mathsf{a}},\mu_{\mathsf{b}}\in\mathbb{R}^{P}$, where “$\leq$” has to be
understood component-wise.

For each admissible parameter $\mu\in\mathcal{P}$, $a_{\mu}:V\times
V\to\mathbb{R}$ denotes a continuous and coercive bilinear form,
$l_{\mu},j_{\mu}:V\to\mathbb{R}$ are continuous linear functionals and
$k_{\mu}:V\times V\to\mathbb{R}$ denotes a continuous symmetric bilinear form.
The primal residual of (P.b) is key for the optimization as well as for a
posteriori error estimation. We define for given $u\in V$,
$\mu\in\mathcal{P}$, the primal residual $r_{\mu}^{\textnormal{pr}}(u)\in
V^{\prime}$ associated with (P.b) by

(1) $\displaystyle r_{\mu}^{\textnormal{pr}}(u)[v]:=l_{\mu}(v)-a_{\mu}(u,v)$
$\displaystyle\text{for all }v\in V.$

Following the approach of _first-optimize-then-discretize_ , we base our
discretization and model order reduction approach on the first order necessary
optimality system, i.e. (cf. [2] for details and further references)

(2a) $\displaystyle r_{\bar{\mu}}^{\textnormal{pr}}(\bar{u})[v]$
$\displaystyle=0$ $\displaystyle\text{for all }v\in V,$ (2b)
$\displaystyle\partial_{u}\mathcal{J}(\bar{u},\bar{\mu})[v]-a_{\mu}(v,\bar{p})$
$\displaystyle=0$ $\displaystyle\text{for all }v\in V,$ (2c)
$\displaystyle(\partial_{\mu}\mathcal{J}(\bar{u},\bar{\mu})+\nabla_{\mu}r^{\textnormal{pr}}_{\bar{\mu}}(\bar{u})[\bar{p}])\cdot(\nu-\bar{\mu})$
$\displaystyle\geq 0$ $\displaystyle\text{for all }\nu\in\mathcal{P}.$

From (2b) we deduce the so-called _adjoint – or dual – equation_

(3) $\displaystyle
a_{\mu}(q,p_{\mu})=\partial_{u}\mathcal{J}(u_{\mu},\mu)[q]=j_{\mu}(q)+2k_{\mu}(q,u_{\mu})$
$\displaystyle\text{for all }q\in V,$

with solution $p_{\mu}\in V$ for a fixed $\mu\in\mathcal{P}$ and given the
solution $u_{\mu}\in V$ to the state equation (P.b). For given $u,p\in V$, we
introduce the dual residual $r_{\mu}^{\textnormal{du}}(u,p)\in V^{\prime}$
associated with (3) as

(4) $\displaystyle
r_{\mu}^{\textnormal{du}}(u,p)[q]:=j_{\mu}(q)+2k_{\mu}(q,u)-a_{\mu}(q,p)$
$\displaystyle\text{for all }q\in V.$

###  2.1. Petrov-Galerkin based discretization and model reduction

Assuming $V_{h}^{\textnormal{pr}},V_{h}^{\textnormal{du}}\subset V$ to be a
finite-dimensional subspaces, we define a Petrov-Galerkin projection of (P)
onto $V_{h}^{\textnormal{pr}},V_{h}^{\textnormal{du}}$ by considering, for
each $\mu\in\mathcal{P}$, the solution $u_{h,\mu}\in V_{h}^{\textnormal{pr}}$
of the _discrete primal equation_

(5) $\displaystyle a_{\mu}(u_{h,\mu},v_{h})=l_{\mu}(v_{h})$
$\displaystyle\text{for all }v_{h}\in V_{h}^{\textnormal{du}},$

and for given $\mu,u_{h,\mu}$, the solution $p_{h,\mu}\in
V_{h}^{\textnormal{du}}$ of the _discrete dual equation_ as

(6) $\displaystyle
a_{\mu}(q_{h},p_{h,\mu})=\partial_{u}\mathcal{J}(u_{h,\mu},\mu)[q_{h}]=j_{\mu}(q_{h})+2k_{\mu}(q_{h},u_{h,\mu})$
$\displaystyle\forall q_{h}\in V_{h}^{\textnormal{pr}}.$

Note that the test space of one equation corresponds with the ansatz space of
the other equation. In order to obtain a quadratic system, we require $\dim
V_{h}^{\textnormal{pr}}=\dim V_{h}^{\textnormal{du}}$. Note that a Ritz-
Galerkin projection is obtained, if the primal and dual discrete spaces
coincide, i.e. $V_{h}^{\textnormal{pr}}=V_{h}^{\textnormal{du}}$.

Given problem adapted reduced basis (RB) spaces
$V_{r}^{\textnormal{pr}}\subset
V_{h}^{\textnormal{pr}},V_{r}^{\textnormal{du}}\subset
V_{h}^{\textnormal{du}}$ of the same low dimension $n:=\dim
V_{r}^{\textnormal{pr}}=\dim V_{r}^{\textnormal{du}}$ we obtain the reduced
versions for the optimality system as follows:

  * •

PG-RB approximation for (2a): 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

(7a) $\displaystyle a_{\mu}(u_{{r},\mu},v_{r})=l_{\mu}(v_{r})$
$\displaystyle\qquad\text{for all }v_{r}\in V_{r}^{\textnormal{du}}.$

  * •

PG-RB approximation for (2b): 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_

(7b) $\displaystyle
a_{\mu}(q_{r},p_{{r},\mu})=\partial_{u}\mathcal{J}(u_{{r},\mu},\mu)[q_{r}]=j_{\mu}(q_{r})+2k_{\mu}(q_{r},u_{{r},\mu})$
$\displaystyle\forall q_{r}\in V_{r}^{\textnormal{pr}}.$

We define the PG-RB reduced optimization functional by

(8)
$\displaystyle{{\hat{\mathcal{J}}_{r}}}(\mu):=\mathcal{J}(u_{{r},\mu},\mu)$

with $u_{{r},\mu}\in V_{r}^{\textnormal{pr}}$ being the solution of (7a). We
then consider the _RB reduced optimization problem_ by finding a locally
optimal solution $\bar{\mu}_{r}$ of

($\hat{\textnormal{P}}_{r}$)
$\displaystyle\min_{\mu\in\mathcal{P}}{{\hat{\mathcal{J}}_{r}}}(\mu).$

Note that in contrast to the NCD-approach that has been introduced in [11, 2],
we do not need to correct the reduced functional in our PG-RB approach, as the
primal and dual solutions automatically satisfy
$r_{\mu}^{\textnormal{pr}}(u_{{r},\mu})[p_{{r},\mu}]=0$. Actually, this is the
main motivation for the usage of the Petrov-Galerkin approach in this
contribution. In particular this results in the possibility to compute the
gradient of the reduced functional with respect to the parameters solely based
on the primal and dual solution of the PG-RB approximation, i.e.

(9)
$\displaystyle\big{(}\nabla_{\mu}{{\hat{\mathcal{J}}_{r}}}(\mu)\big{)}_{i}$
$\displaystyle=\partial_{\mu_{i}}\mathcal{J}(u_{{r},\mu},\mu)+\partial_{\mu_{i}}r_{\mu}^{\textnormal{pr}}(u_{{r},\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 PG-RB primal and dual reduced solutions of (7a) and (7b), respectively.
Note that in [17], the formula in (9) was motivated by replacing the full
order functions by their respective reduced counterpart which resulted in an
inexact gradient for the non-conforming approach. In the PG setting, (9)
instead defines the true gradient of ${{\hat{\mathcal{J}}_{r}}}$ without
having to add a correction term as proposed in [11, 2]. This also holds for
the true Hessian which can be computed by also replacing all full order
counterparts in the full order Hessian. In this contribution we however only
focus on quasi-Newton methods that do not require a reduced Hessian.

###  2.2. A posteriori error estimation for an error aware algorithm

In order to construct an adaptive trust region algorithm we require a
posteriori error estimation for (7a), (7b) and (8). For that, we can utilize
standard residual based estimation. For an a posteriori result of the gradient
$\nabla_{\mu}{{\hat{\mathcal{J}}_{r}}}(\mu)$, we also refer to [11].

######  Proposition 2.1 (Upper error bound for the reduced quantities).

For $\mu\in\mathcal{P}$, let $u_{h,\mu}\in V_{h}^{\textnormal{pr}}$ and
$p_{h,\mu}\in V_{h}^{\textnormal{du}}$ be solutions of (5) and (6) and let
$u_{{r},\mu}\in V_{r}^{\textnormal{pr}}$, $p_{{r},\mu}\in
V_{r}^{\textnormal{du}}$ be a solution of (7a), (7b). Then it holds

  1. (i)

$\|u_{h,\mu}-u_{{r},\mu}\|\leq\Delta_{\textnormal{pr}}(\mu):=\alpha_{\mu}^{-1}\,\|r_{\mu}^{\textnormal{pr}}(u_{{r},\mu})\|$,

  2. (ii)

$\|p_{h,\mu}-p_{{r},\mu}\|\leq\Delta_{\textnormal{du}}(\mu):=\alpha_{\mu}^{-1}\big{(}2\gamma_{k_{\mu}}\;\Delta_{\textnormal{pr}}(\mu)+\|r_{\mu}^{\textnormal{du}}(u_{{r},\mu},p_{{r},\mu})\|\Big{)}$,

  3. (iii)

$|\hat{\mathcal{J}}_{h}(\mu)-{{\hat{\mathcal{J}}_{r}}}(\mu)|\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}}$

where $\alpha_{\mu}$ and $\gamma_{k_{\mu}}$ define the inf-sup stability
constant of $a_{\mu}$ and the continuity constant of $k_{\mu}$, respectively.

###### Proof.

For $(i)$ and $(ii)$ we rely on the inf-sup stability of $a_{\mu}$ on
$V_{h}^{\textnormal{pr}}$ and $V_{h}^{\textnormal{du}}$ and proceed
analogously to [17]. For $(iii)$, we refer to [11] since
$r_{\mu}^{\textnormal{pr}}(u_{{r},\mu})[p_{{r},\mu}]=0$. ∎

It is important to mention that the computation of these error estimators
include the computation of the (parameter dependent) inf-sup constant of
$a_{\mu}$ which involves an eigenvalue problem on the FOM level. In practice,
cheaper techniques such as the successive constraint method [10] can be used.
For the conforming approach $V_{h}^{\textnormal{pr}}=V_{h}^{\textnormal{du}}$,
the inf-sup constant is equivalent to the coercivity constant
$\underline{a_{\mu}}$ of $a_{\mu}$ which can be cheaply bounded from below
with the help of the min-theta approach, c.f. [8].

###  2.3. Trust-Region optimization approach and adaptive enrichment

Error aware Trust-Region - Reduced Basis methods (TR-RB) with several
different advances and features have been extensively studied e.g. in [17, 11,
2]. They iteratively compute a first-order critical point of problem (P). For
each outer iteration $k\geq 0$ of the TR method, we consider a model function
$m^{(k)}$ as a cheap local approximation of the quadratic cost functional
$\mathcal{J}$ in the so-called trust-region, which has radius $\varrho^{(k)}$
and can be characterized by the a posteriori error estimator. In our approach
we choose

$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. Thus, we can use
$\Delta_{{{\hat{\mathcal{J}}_{r}}}}(\mu)$ for characterizing the trust-region.
We are therefore interested in solving the following error aware constrained
optimization sub-problem

(10)
$\displaystyle\min_{\widetilde{\mu}\in\mathcal{P}}{{\hat{\mathcal{J}}_{r}}}^{(k)}(\widetilde{\mu})\quad\text{
s.t. }$
$\displaystyle\frac{\Delta_{\hat{\mathcal{J}}}(\widetilde{\mu})}{{{\hat{\mathcal{J}}_{r}}}^{(k)}(\widetilde{\mu})}\leq\varrho^{(k)},\quad\widetilde{\mu}:=\mu^{(k)}+s\in\mathcal{P}$
$\displaystyle\text{ and
}r_{\tilde{\mu}}^{\textnormal{pr}}(u_{\tilde{\mu}})[v]=0\,\text{ for all }v\in
V.$

In our TR-RB algorithm we build on the algorithm in [11]. In the sequel, we
only summarize the main features of the algorithm and refer to the source for
more details. We initialize the RB spaces with the starting parameter
$u_{\mu^{(0)}}$, i.e.
$V^{\textnormal{pr},(0)}_{r}=\big{\\{}u_{h,\mu^{(0)}}\big{\\}}$ and
$V^{\textnormal{du},(0)}_{r}=\big{\\{}p_{h,\mu^{(0)}}\big{\\}}$. For every
iteration point $\mu_{k}$, we locally solve (10) with the quasi-Newton
projected BFGS algorithm combined with an Armijo-type condition and terminate
with a standard reduced FOC termination criteria, modified with a projection
on the parameter space $\mathrm{P}_{\mathcal{P}}$ to account for constraints
on the parameter space. Additionally, we use a second boundary termination
criteria for preventing the subproblem from spending too much computational
time on the boundary of the trust region. After the next iterate $\mu_{k+1}$
has been computed, the sufficient decrease conditions helps to decide whether
to accept the iterate:

(11)
$\displaystyle{{\hat{\mathcal{J}}_{r}}}^{(k+1)}(\mu^{(k+1)})\leq{{\hat{\mathcal{J}}_{r}}}^{(k)}(\mu_{\text{\rm{AGC}}}^{(k)})$
$\displaystyle\text{ for all }k\in\mathbb{N}$

This condition can be cheaply checked with the help of an sufficient and
necessary condition. If $\mu_{k+1}$ is accepted, we use the parameter to
enrich the RB spaces in a Lagrangian manner, i.e.
$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}\\}.$
For more basis constructions, we refer to [2] where also a strategy is
proposed to skip an enrichment. After the enrichment, overall convergence of
the algorithm can be checked with a FOM-type FOC condition

$\|\mu^{(k+1)}-\mathrm{P}_{\mathcal{P}}(\mu^{(k+1)}-\nabla_{\mu}\hat{\mathcal{J}}_{h}(\mu^{(k+1)}))\|_{2}\leq\tau_{\text{\rm{FOC}}},$

where the FOM quantities are available from the enrichment. Moreover, these
FOM quantities allow for computing a condition for possibly enlargement of the
TR radius if the reduced model is better than expected. For this, we use

(12)
$\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 the described algorithm the following convergence result holds.

######  Theorem 2.2 (Convergence of the TR-RB algorithm, c.f. [2]).

For sufficient assumptions on the Armijo search to solve (10), every
accumulation point $\bar{\mu}$ of the sequence
$\\{\mu^{(k)}\\}_{k\in\mathbb{N}}\subset\mathcal{P}$ generated by the above
described TR-RB algorithm is an approximate first-order critical point for
$\hat{\mathcal{J}}_{h}$, i.e., it holds

(13)
$\|\bar{\mu}-\mathrm{P}_{\mathcal{P}}(\bar{\mu}-\nabla_{\mu}\hat{\mathcal{J}}_{h}(\bar{\mu}))\|_{2}=0.$

Note that the above Theorem is not affected from the choice of the reduced
primal and dual equations as long as Proposition 2.1 holds. In fact, the
result can also be used for the TR-RB with the proposed PG approximation of
these equations.

Let us also mention that it has been shown in [2] that a projected Newton
algorithm for the subproblems can enhance the convergence speed and accuracy.
Moreover, a reduced Hessian can be used to introduce an a posteriori result
for the optimal parameter which can be used as post processing. In the work at
hand, we neglect to transfer the ideas from [2] for the PG variant.

##  3\. Numerical experiments

In this section, we analyze the behavior of the proposed PG variant of the TR-
RB (BFGS PG TR-RB) algorithm. We aim to compare the computational time and the
accuracy of the PG variant to existing approaches from the literature. To this
end, we mention the different algorithms that we compare in this contribution:
Projected BFGS FOM: As FOM reference optimization method, we consider a
standard projected BFGS method, which uses FOM evaluations for all required
quantities.
Non-conforming TR-RB algorithm from [11] (BFGS NCD TR-RB): For comparison, we
choose the above described TR-RB algorithm but with a non conforming choice of
the RB spaces, which means that (7a) and (7b) are Galerkin projected equations
where the test space coincides with the respective ansatz space. This requires
to use the NCD-corrected reduced objective functional
$\hat{\mathcal{J}}_{r}(\mu)+r_{\mu}^{\textnormal{pr}}(u_{{r},\mu})[p_{{r},\mu}]$
and requires additional quantities for computing the Gradient.

For computational details including the choice of all relevant tolerances for
both TR-RB algorithms we again refer to [11], where all details can be found.
We only differ in the choice of the stopping tolerance
$\tau_{\text{\rm{FOC}}}=10^{-6}$. Also note that, as pointed out in Section
2.2, we use the same error estimators as in [11] which include the coercivity
constant instead of the inf-sub constant. The source code for the presented
experiments can be found in [12], which is a revised version of the source
code for [11] and also contains detailed interactive jupyter-
notebooks111Available at https://github.com/TiKeil/Petrov-Galerkin-TR-RB-for-
pde-opt..

###  3.1. Model problem: Quadratic objective functional with elliptic PDE
constraints

For our numerical evaluation we reconsider Experiment 1 from [2]. We set the
objective functional to be a weighted $L^{2}$-misfit on a domain of interest
$D\subseteq\Omega$ with a weighted Tikhonov term, i.e.

$\displaystyle\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,$

with a desired state $u^{\text{d}}$ and desired parameter $\mu^{\text{d}}$. We
added the constant term $1$ to verify that $\mathcal{J}>0$. With respect to
the formulation in (P.a), we have
$\Theta(\mu)=\frac{\sigma_{d}}{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}}$,
$j_{\mu}(u)=-\sigma_{d}\int_{D}u^{\text{d}}u$, and
$k_{\mu}(u,u)=\frac{\sigma_{d}}{2}\int_{D}u^{2}$. From this general choice we
can construct several applications. One instance is to consider the stationary
heat equation as equally constraints, i.e. the weak formulation of the
parameterized equation

(14) $\begin{split}-\nabla\cdot\left(\kappa_{\mu}\nabla
u_{\mu}\right)&=f_{\mu}\hskip 58.0pt\text{in }\Omega,\\\
c_{\mu}(\kappa_{\mu}\nabla u_{\mu}\cdot n)&=(u_{\text{out}}-u_{\mu})\hskip
20.0pt\text{on }\partial\Omega.\end{split}$

with parametric diffusion coefficient $\kappa_{\mu}$ and source $f_{\mu}$,
outside temperature $u_{\text{out}}$ and Robin function $c_{\mu}$. From this
the bilinear and linear forms $a_{\mu}$ and $l_{\mu}$ can easily be deduced
and we set the parameter box constraints
$\mu_{i}\in[\mu_{i}^{\text{min}},\mu_{i}^{\text{max}}].$

As an application of (14), we use the blueprint of a building with windows,
heaters, doors and walls, which can be parameterized according to Figure 1. We
picked a certain domain of interest $D$ and we enumerated all windows, walls,
doors and heaters separately.

Figure 1. Parameterization based on a blueprint of a building floor, see [11]
for details. Numbers indicate potential parameters, where $i.$ is a window or
a heater, $\underline{j}$ are doors, and $k|$ are walls. The green region
illustrates the domain of interest $D$.

We set the computational domain to
$\Omega:=[0,2]\times[0,1]\subset\mathbb{R}^{2}$ and we model boundary
conditions by incorporating all walls and windows that touch the boundary of
the blueprint to the Robin function $c_{\mu}$. All other diffusion components
enter the diffusion coefficient $\kappa_{\mu}$, whereas the heaters are
incorporated as a source term on the right hand side $f_{\mu}$. Moreover, we
assume an outside temperature of $u_{\text{out}}=5$. For our discretization we
choose the conforming case for the FOM, i.e.
$V_{h}^{\textnormal{pr}}=V_{h}^{\textnormal{du}}$, and use a mesh size
$h=\sqrt{2}/200$ which resolves all features from the given picture and
results in $\dim V_{h}^{\textnormal{pr}}=80601$ degrees of freedom.

For the parameterization, we choose 2 doors, 7 heaters and 3 walls which
results in a parameter space of 12 dimensions. More details can be viewed in
the accompanying source code.

###  3.2. Analysis of the error behavior

This section aims to show and discuss the model reduction error for the
proposed Petrov–Galerkin approach. Furthermore, we compare it to the Galerkin
strategy from the NCD-corrected approach, where we denote the solutions by
$u_{{r}}^{G}$, $p_{{r}}^{G}$ and the corrected functional by
${{\hat{\mathcal{J}}_{r}}}^{\text{NCD}}$, as well as to the non corrected
approach from [17] whose functional we denote by $\hat{J}_{r}$. For this, we
employ a standard goal oriented Greedy search algorithm as also done in [11,
Section 4.3.1] with the relative a posteriori error of the objective
functional
$\Delta_{{{\hat{\mathcal{J}}_{r}}}}(\mu)/{{\hat{\mathcal{J}}_{r}}}(\mu)$. As
pointed out in Section 2.2, due to
$V_{h}^{\textnormal{pr}}=V_{h}^{\textnormal{du}}$, we can replace the inf-sup
constant by a lower bound for the coercivity constant of the conforming
approach. It is also important to mention that for this experiment, we have
simplified our objective functional $\mathcal{J}$ by setting the domain of
interest to the whole domain $D\equiv\Omega$. As a result the dual problem is
simpler which enhances the stability of the PG approach. The reason for that
is further discussed below. Figure 2 shows the difference in the decay and
accuracy of the different approaches. It can clearly be seen that the NCD-
corrected approach remains the most accurate approach while the objective
functional and the gradient of the PG approach shows a better approximation
compared to the non corrected version. Clearly the PG approximation of the
primal and dual solutions are less accurate and at least the primal error
decays sufficiently.

$20$$40$$60$$10^{-5}$$10^{-1}$$10^{3}$$20$$40$$60$greedy extension stepgreedy
extension step(A) functional $\hat{\mathcal{J}}_{h}$ and gradient
$\nabla\hat{\mathcal{J}}_{h}$(B) primal and dual solution

$|\hat{\mathcal{J}}_{h}-\hat{J}_{r}|$$|\nabla\hat{\mathcal{J}}_{h}-\tilde{\nabla}\hat{J}_{r}|$$|u_{h,\mu}-u^{G}_{{r},\mu}|$$|p_{h,\mu}-p^{G}_{{r},\mu}|$$|\hat{\mathcal{J}}_{h}-{{\hat{\mathcal{J}}_{r}}}^{\text{NCD}}|$$|\nabla\hat{\mathcal{J}}_{h}-\nabla{{\hat{\mathcal{J}}_{r}}}^{\text{NCD}}|$$|u_{h,\mu}-u_{{r},\mu}|$$|p_{h,\mu}-p_{{r},\mu}|$$|\hat{\mathcal{J}}_{h}-{{\hat{\mathcal{J}}_{r}}}|$$|\nabla\hat{\mathcal{J}}_{h}-\nabla{{\hat{\mathcal{J}}_{r}}}|$

Figure 2. Evolution of the true reduction error in the reduced functional and
gradient and its approximations (A) and the primal and dual solutions and its
approximations (B), during adaptive greedy basis generation. Depicted is the
$L^{\infty}(\mathcal{P}_{\textnormal{val}})$-error for a validation set
$\mathcal{P}_{\textnormal{val}}\subset\mathcal{P}$ of $100$ randomly selected
parameters, i.e. $|\hat{\mathcal{J}}_{h}-\hat{J}_{r}|$ corresponds to
$\max_{\mu\in\mathcal{P}_{\textnormal{val}}}|\hat{\mathcal{J}}_{h}(\mu)-\hat{J}_{r}(\mu)|$,
and so forth.

While producing the shown error study, we experienced instabilities of the PG
reduced primal and dual systems. This is due to the fact that for very complex
dual problems, the test space of each problem very poorly fits to the
respective ansatz space. In fact, the decay of the primal error that can be
seen in Figure 2(B) can not be expected in general and is indeed a consequence
of the simplification of the objective functional, where we chose
$D\equiv\Omega$. The stability problems can already be seen in the dual error
and partly on the primal error for basis size 20. In general it can happen
that the reduced systems are highly unstable for specific parameter values.
Instead, a sophisticated Greedy Algorithm for deducing appropriate reduced
spaces may require to build a larger dual or primal space by adding
stabilizing snapshots, i.e. search for supremizers [1]. Since the aim of the
work at hand is not to provide an appropriate Greedy based algorithm, we
instead decided to reduce the complexity of the functional where stability
issues are less present.

###  3.3. TR-RB algorithm

We now compare the PG variant of the TR-RB with the above mentioned NCD-
corrected TR-RB approach. Importantly, we use the original version of the
model problem, i.e. we pick the domain of interest to be defined as suggested
in Figure 1. In fact, we accept the possibility of high instabilities in the
reduced model. We pick ten random starting parameter, perform both algorithms
and compare the averaged result. In Figure 3 one particular starting parameter
is depicted and in Table 1 the averaged results for all ten optimization runs
are shown.

$0$$200$$400$$10^{-6}$$10^{-3}$$10^{0}$time in seconds [s]$\displaystyle\|\overline{\mu}-\mu^{(k)}\|^{\text{rel}}_{2}$BFGS NCD TR-RB [11]BFGS PG TR-RB Figure 3. Relative error decay w.r.t. the optimal parameter $\bar{\mu}$ and performance of selected algorithms for a single optimization run with random initial guess $\mu^{(0)}$ for $\tau_{\text{FOC}}=10^{-6}$. | runtime[s] |  |  iterations $k$ |  |
---|---|---|---|---|---
| avg. (min/max) | speed-up | avg. (min/max) | rel. error | FOC cond.
FOM BFGS | 6955 (4375/15556) | – | 471.44 (349/799) | 3.98e-5 | 3.36e-6
TR NCD-BFGS | 171 (135/215) | 40.72 |  12.56(11/15) | 4.56e-6 | 6.05e-7
TR PG-BFGS | 424 (183/609) | 16.39 |  17.56(11/22) | 4.62e-6 | 7.57e-7
Table 1. Performance and accuracy of the algorithms for ten optimization runs
with random initial guess $\mu^{(0)}$ and $\tau_{\text{\rm{FOC}}}=10^{-6}$.

It can be seen that the PG variant is a valid approach and converges
sufficiently fast with respect to the FOM BFGS method. Surely, it can not be
said that this stability issues do not enter the performance of the proposed
TR-RB methodology but regardless of the stability of the reduced system, we
note that the convergence result in Theorem 1 still holds true. However, the
instability of the reduced systems clearly harms the algorithm from iterating
as fast as the NCD-corrected approach. One reason for that is that the trust
region is much larger for the NCD-corrected approach, allowing the method to
step faster. We would like to emphasize that the depicted result in Figure 3
is neither an instance of the worst nor the best performance of the PG
approach but rather an intermediate performance. The comparison highly depends
on the starting parameter, and the structure of the optimization problem. As
discussed above, the suggested PG approach can potentially benefit from more
involved enrichment strategies that account for the mentioned stability
issues.

Last but not least, it is important to mention that the above experiment
showed weaknesses of the chosen projected BFGS approach as FOM method as well
as for the TR-RB subproblems which has been extensively studied in [2].
Instead, it is beneficial to choose higher order optimization methods, such as
projected Newton type methods, where also the PG variant can be applied to
neglect the more involved contributions from the NCD-corrected approach that
are entering the Hessian approximation.

##  4\. Concluding remarks

In this contribution we demonstrated, how adaptive enrichment based on
rigorous a posteriori error control can be used within a Trust-Region-Reduced-
Basis approach to speed up the solution of large scale PDE constrained
optimization problems. Within this approach the reduced approximation spaces
are tailored towards the solution of the overall optimization problem and thus
circumvent the offline construction of a reduced model with good approximation
properties with respect to the whole parameter regime. In particular, we
compared a new Petrov-Galerkin approach with the NCD-Galerkin approach that
has recently been introduced in [2]. Our results demonstrate the benefits of
the Petrov-Galerkin approach with respect to the approximation of the
objective functional and its derivatives as well as with respect to
corresponding a posteriori error estimation. However, the results also show
deficiencies with respect to possible stability issues. Thus, further
improvements of the enrichment strategy are needed in order to guarantee
uniform boundedness of the reduced inf-sup constants.

## References

  * [1] F. Ballarin, A. Manzoni, A. Quarteroni, and G. Rozza.  Supremizer stabilization of pod–galerkin approximation of parametrized steady incompressible navier–stokes equations.  International Journal for Numerical Methods in Engineering, 102(5):1136–1161, 2015.
  * [2] S. Banholzer, T. Keil, L. Mechelli, M. Ohlberger, F. Schindler, and S. Volkwein.  An adaptive projected newton non-conforming dual approach for trust-region reduced basis approximation of pde-constrained parameter optimization, 2020.  arXiv 2012.11653.
  * [3] M. Barrault, Y. Maday, N. C. Nguyen, and A. T. Patera.  An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations.  C. R. Math., 339(9):667–672, 2004.
  * [4] P. Benner, M. Ohlberger, A. Patera, G. Rozza, and K. Urban, editors.  Model reduction of parametrized systems, volume 17 of MS&A. Modeling, Simulation and Applications.  Springer, Cham, 2017.  Selected papers from the 3rd MoRePaS Conference held at the International School for Advanced Studies (SISSA), Trieste, October 13–16, 2015\.
  * [5] A. Buhr, L. Iapichino, M. Ohlberger, S. Rave, F. Schindler, and K. Smetana.  Localized model reduction for parameterized problems, pages 245–306.  De Gruyter, 2021.
  * [6] M. Drohmann, B. Haasdonk, and M. Ohlberger.  Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation.  SIAM Journal on Scientific Computing, 34(2):A937–A969, 2012.
  * [7] J. Feinauer, S. Hein, S. Rave, S. Schmidt, D. Westhoff, J. Zausch, O. Iliev, A. Latz, M. Ohlberger, and V. Schmidt.  MULTIBAT: Unified workflow for fast electrochemical 3d simulations of lithium-ion cells combining virtual stochastic microstructures, electrochemical degradation models and model order reduction.  Journal of Computational Science, 31:172–184, 2019.
  * [8] B. Haasdonk.  Reduced basis methods for parametrized PDEs – A tutorial introduction for stationary and instationary problems.  Technical report, 2014.  Chapter to appear in P. Benner, A. Cohen, M. Ohlberger and K. Willcox: ”Model Reduction and Approximation: Theory and Algorithms”, SIAM.
  * [9] B. Haasdonk and M. Ohlberger.  Reduced basis method for finite volume approximations of parametrized linear evolution equations.  Mathematical Modelling and Numerical Analysis, 42(2):277–302, 2008\.
  * [10] D. Huynh, D. Knezevic, A. Patera, and H. Li.  Methods and apparatus for constructing and analyzing component-based models of engineering systems, 2015.  US Patent 9,213,788.
  * [11] T. Keil, L. Mechelli, M. Ohlberger, F. Schindler, and S. Volkwein.  A non-conforming dual approach for adaptive trust-region reduced basis approximation of pde-constrained optimization, 2020.  arXiv 2006.09297. To appear in: ESAIM: M2AN (2021).
  * [12] T. Keil and M. Ohlberger.  Software for Model Reduction for Large Scale Systems, url: https://doi.org/10.5281/zenodo.4627971, Mar. 2021.
  * [13] M. Ohlberger and S. Rave.  Reduced basis methods: Success, limitations and future challenges.  Proceedings of the Conference Algoritmy, pages 1–12, 2016.
  * [14] M. Ohlberger, M. Schaefer, and F. Schindler.  Localized model reduction in PDE constrained optimization.  International Series of Numerical Mathematics, 169:143–163, 2018\.
  * [15] M. Ohlberger and F. Schindler.  Error control for the localized reduced basis multiscale method with adaptive on-line enrichment.  SIAM Journal on Scientific Computing, 37(6):A2865–A2895, 2015.
  * [16] M. Ohlberger and F. Schindler.  Non-conforming localized model reduction with online enrichment: towards optimal complexity in PDE constrained optimization.  In Finite volumes for complex applications VIII—hyperbolic, elliptic and parabolic problems, volume 200 of Springer Proc. Math. Stat., pages 357–365. Springer, Cham, 2017.
  * [17] E. Qian, M. Grepl, K. Veroy, and K. Willcox.  A certified trust region reduced basis approach to PDE-constrained optimization.  SIAM J. Sci. Comput., 39(5):S434–S460, 2017.