ROM-based multiobjective optimization of elliptic PDEs via numerical continuation

TL;DR

This paper introduces a ROM-based continuation method for multiobjective PDE-constrained optimization, enabling efficient computation of Pareto sets with multiple objectives by reducing computational costs through surrogate modeling.

Contribution

It combines reduced basis model reduction with a continuation approach to handle multiple objectives efficiently in PDE-constrained optimization.

Findings

Significant reduction in computational time for multiobjective PDE optimization.

Ability to handle an arbitrary number of objectives.

Effective surrogate modeling for expensive PDE evaluations.

Abstract

Multiobjective optimization plays an increasingly important role in modern applications, where several objectives are often of equal importance. The task in multiobjective optimization and multiobjective optimal control is therefore to compute the set of optimal compromises (the Pareto set) between the conflicting objectives. Since the Pareto set generally consists of an infinite number of solutions, the computational effort can quickly become challenging which is particularly problematic when the objectives are costly to evaluate as is the case for models governed by partial differential equations (PDEs). To decrease the numerical effort to an affordable amount, surrogate models can be used to replace the expensive PDE evaluations. Existing multiobjective optimization methods using model reduction are limited either to low parameter dimensions or to few (ideally two) objectives. In this article, we present a combination of the reduced basis model reduction method with a continuation approach using inexact gradients. The resulting approach can handle an arbitrary number of objectives while yielding a significant reduction in computing time.

Tables

Table 1. Table 1 : Number of basis functions for different error bounds

Error bound	Strong Greedy	Error Estimate
$\epsilon =1e-6$	24	56
$\epsilon =1e-5$	20	50
$\epsilon =1e-4$	16	40
$\epsilon =1e-3$	12	32
$\epsilon =1e-2$	12	26
$\epsilon =1e-1$	10	20

Table 2. Table 2 : Comparison of the performance of the exact continuation method, Strategy 1 and Strategy 2 for Example ( 46 ). The number of subproblems is split up in subproblems for the continuation and initialization (cf. Section 2.4 )

Algorithm	# Boxes	# Subproblems	Runtime (in seconds)
Exact cont.	$15916$	$18721+25$	$17501s$
Strategy 1	$21750$	$24490+25$	$1426s$
Strategy 2	$899$	$1027+225$	$276s$

Equations

$$\underset{u\in {\mathbb{R}}^n}{\min }J(u).$$

$$\alpha \in {\Delta }_k:=\left\{\alpha \in ({\mathbb{R}}^{\ge 0}{)}^k:\sum _{i=1}^k{\alpha }_i=1\right\}$$

$$DJ(\bar{u}{)}^{\top }\alpha =\sum _{i=1}^k{\alpha }_i\nabla J_i(\bar{u})=0.$$

$$F:{\mathbb{R}}^n\times ({\mathbb{R}}^{>0}{)}^k\to {\mathbb{R}}^{n+1},(u,\alpha )\mapsto \left(\begin{array}{c}{\sum }_{i=1}^k{\alpha }_i\nabla J_i(u)\\ 1-{\sum }_{i=1}^k{\alpha }_i\end{array}\right),$$

$$rk(DF(u,\alpha ))=n+1\forall (u,\alpha )\in \mathcal{M},$$

$$T_{(u,\alpha )}\mathcal{M}=ker(DF(u,\alpha )).$$

$$\mathcal{B}(r)$$

$${\mathcal{B}}_c(r):=\{B\in \mathcal{B}(r):P_c\cap B\ne \varnothing \}.$$

$$N(u,r):=\{B\in \mathcal{B}(r):B(u,r)\cap B\ne \varnothing \}$$

$${\mathcal{B}}^{\prime }(\bar{u},r)=\{B^{\prime }\in \mathcal{B}(r):B^{\prime }\in N(\bar{u},r),B^{\prime }\cap \bar{u}+p{r}_u(T_{(\bar{u},\bar{\alpha })}\mathcal{M})\ne \varnothing \}.$$

$$\underset{u\in B,\alpha \in {\Delta }_k}{\min }\Vert DJ(u{)}^{\top }\alpha {\Vert }_2^2$$

$$\theta (B)=0\iff B\cap P_c\ne \varnothing .$$

$$\underset{u\in {\mathbb{R}}^n}{\sup }{\Vert \nabla J_i(u)-\nabla J_i^r(u)\Vert }_2\le {ϵ}_i,i\in \{1,...,k\},$$

$$\Vert D{J}^r(\bar{u}{)}^{\top }\bar{\alpha }{\Vert }_2\le \sum _{i=1}^k{\bar{\alpha }}_i{ϵ}_i\le \Vert ϵ{\Vert }_{\infty }.$$

$${\Vert D{J}^r(\bar{u}{)}^{\top }\bar{\alpha }\Vert }_2$$

$$\le \sum _{i=1}^k{\Vert \nabla J_i^r(\bar{u})-\nabla J_i(\bar{u})\Vert }_2{\bar{\alpha }}_i\le \sum _{i=1}^k{\bar{\alpha }}_i{ϵ}_i\le {\Vert ϵ\Vert }_{\infty }$$

$$P_1^r$$

$$P_2^r$$

$$P_c\subseteq P_2^r\subseteq P_1^r\text{ and }{P}_c^r\subseteq P_2^r\subseteq P_1^r,$$

$$\underset{u\in {\mathbb{R}}^n}{\sup }\Vert \nabla {\tilde{J}}_i(u)-\nabla J_i^r(u){\Vert }_2\le {ϵ}_i\forall i\in \{1,...,k\}$$

$$\tilde{\alpha }$$

$$\nu$$

$$g(u)$$

$$\tilde{J}(u)$$

$$\Vert \nabla {\tilde{J}}_i(u)-\nabla J_i^r(u){\Vert }_2=\Vert \nabla g_i(u){\Vert }_2=\frac{{ϵ}_i}{{\tilde{\alpha }}^{\top }ϵ}\Vert \nu {\Vert }_2\le {ϵ}_i\forall u\in {\mathbb{R}}^n\text{ and }\forall i\in \{1,...,k\},$$

$$D\tilde{J}(\tilde{u}{)}^{\top }\tilde{\alpha }=\nu +\sum _{i=1}^k{\tilde{\alpha }}_i\nabla g_i(\tilde{u})=\nu -\sum _{i=1}^k{\tilde{\alpha }}_i\frac{{ϵ}_i}{{\tilde{\alpha }}^{\top }ϵ}\nu =0,$$

$$J^r:{\mathbb{R}}^2\to {\mathbb{R}}^2,u\mapsto \left(\begin{array}{c}(u_1-1{)}^2+(u_2-1{)}^4\\ (u_1+1{)}^2+(u_2+1{)}^2\end{array}\right).$$

$$A:=\left\{u\in {\mathbb{R}}^n:\underset{\alpha \in {\Delta }_k}{\min }\left(\Vert D{J}^r(u{)}^{\top }\alpha {\Vert }_2^2-({\alpha }^{\top }ϵ{)}^2\right)<0\right\}.$$

$${\alpha }^i\in {\text{argmin}}_{\alpha \in {\Delta }_k}\left(\Vert D{J}^r(u^i{)}^{\top }\alpha {\Vert }_2^2-({\alpha }^{\top }ϵ{)}^2\right).$$

$$\Psi :{\mathbb{R}}^n\times {\Delta }_k\to \mathbb{R},(u,\alpha )\mapsto \Vert D{J}^r(u{)}^{\top }\alpha {\Vert }_2^2-({\alpha }^{\top }ϵ{)}^2.$$

$$\bar{\alpha }\in {\text{argmin}}_{\alpha \in {\Delta }_k}\left(\Vert D{J}^r(\bar{u}{)}^{\top }\alpha {\Vert }_2^2-({\alpha }^{\top }ϵ{)}^2\right).$$

$$\underset{\alpha \in {\Delta }_k}{\min }\left(\Vert D{J}^r(u{)}^{\top }\alpha {\Vert }_2^2-({\alpha }^{\top }ϵ{)}^2\right)\le \psi (u)<0\forall u\in U$$

$${\mathcal{B}}_c^r(r):=\{B\in \mathcal{B}(r):B\cap P_2^r\ne \varnothing \}.$$

Full text

1 Introduction

The dilemma of deciding between multiple, equally important goals is present in almost all areas of engineering and economy. A prominent example comes from production, where we want to produce a product at minimal cost while simultaneously preserving a high quality. In the same manner, multiple goals are present in most technical applications, maximizing the velocity while minimizing the energy consumption of electric vehicles [24] being only one of many examples. These conflicting goals result in multiobjective optimization problems (MOPs) [9], where we want to optimize all objectives simultaneously. Since the objectives are in general contradictory, there exists an infinite number of optimal compromises. The set of these compromise solutions is called the Pareto set, and the goal in multiobjective optimization is to approximate this set in an efficient manner, which is significantly more expensive than solving a single objective problem. Due to this, the development of efficient numerical approximation methods is an active area of research, and methods range from scalarization [9, 14] over set-oriented approaches [8] and continuation [14] to evolutionary algorithms [7]. Recent advances have paved the way to new challenging application areas for multiobjective optimization such as feedback control or problems constrained by partial differential equations (PDEs); cf. [22] for a survey.

In the presence of PDE constraints, the computational effort can quickly become infeasible such that special means have to be taken in order to accelerate the computation. To this end, surrogate models form a promising approach for significantly reducing the computational effort. A widely used approach is to directly construct a mapping from the parameter to the objective space using as few function evaluations of the expensive model as possible, cf. [30, 6] for extensive reviews. In the case of PDE constraints, an alternative approach is via dimension reduction techniques such as Proper Orthogonal Decomposition (POD) [29, 18] or the reduced basis (RB) method [11]. In these methods, a small number of high-fidelity solutions is used to construct a low-dimensional surrogate model for the PDE which can be evaluated significantly faster while guaranteeing convergence using error estimates. In recent years, several methods have been proposed where model reduction is used in multiobjective optimization and optimal control. In [17] and [16], scalarization using the so-called weighted sum method was combined with RB and POD, respectively. In [1, 2], convex problems were solved using reference point scalarization and POD, and set-oriented approaches were used in [4, 3]. A comparison of both was performed in [23] for the Navier–Stokes equations.

In this article we combine an extension of the continuation methods presented in [14, 28] to inexact gradients (Section 2) with a reduced basis approach for elliptic PDEs (Section 3). To deal with the error introduced by the RB approach, we combine the KKT conditions for MOPs with error estimates for the RB method to obtain a tight superset of the Pareto set. For the example considered here, the proposed method yields a speed-up factor of approximately 63 compared to the direct solution of the expensive problem (Section 4). Additionally, our approach allows us to control the quality of the result by controlling the errors for each objective function individually.

2 A continuation method for MOPs with inexact objective gradients

In this section, we will begin by briefly introducing the basic concepts of multiobjective optimization upon which we will build in this article (see [9, 14] for detailed introductions). Afterwards, we will discuss the continuation method for MOPs and present two modifications of it that can deal with inexact gradient information.

2.1 Multiobjective optimization

The goal of multiobjective optimization is to minimize several conflicting criteria at the same time. In other words, we want to minimize an objective $J=(J_{1},...,J_{k}):\mathbb{R}^{n}\rightarrow\mathbb{R}^{k}$ that is vector valued. It maps the variable space $\mathbb{R}^{n}$ to the image space $\mathbb{R}^{k}$. In contrast to single-objective optimization (i.e., $k=1$), there exists no natural total order of the image space $\mathbb{R}^{k}$ for $k>1$. As a result, the classical concept of optimality has to be generalized:

Definition 2.1**.**

• (a)

$\bar{u}\in\mathbb{R}^{n}$* is called (globally) Pareto optimal if there is no other point $u\in\mathbb{R}^{n}$ such that $J_{i}(u)\leq J_{i}(\bar{u})$ for all $i\in\{1,...,k\}$ and $J_{j}(u)<J_{j}(\bar{u})$ for some $j\in\{1,...,k\}$.*

• (b)

The set $P$ of all Pareto optimal points is called the Pareto set. Its image under $J$ is the Pareto front.

The Pareto set is the solution of the multiobjective optimization problem (MOP)

[TABLE]

Constrained MOPs can be formulated analogously by restricting $u$ in Definition 2.1 to a subset $U\subseteq\mathbb{R}^{n}$. Similar to the scalar-valued case, if $J$ is differentiable, we can use the derivative of $J$ to obtain necessary conditions for Pareto optimality, the Karush-Kuhn-Tucker (KKT) conditions [14]:

Theorem 2.2**.**

Let $\bar{u}$ be a Pareto optimal point of (MOP). Then there exist multipliers

[TABLE]

such that

[TABLE]

For $k=1$, this reduces to the well-known optimality condition $\nabla J(\bar{u})=0$. If $J$ is non-convex, then the points satisfying (KKT) form a proper superset of the Pareto set $P$:

Definition 2.3**.**

If $\bar{u}\in\mathbb{R}^{n}$ and $\bar{\alpha}\in\Delta_{k}$ satisfy (KKT), then $\bar{u}$ is called Pareto critical with corresponding KKT vector $\bar{\alpha}$, containing the KKT multipliers $\bar{\alpha}_{i}$, $i\in\{1,...,k\}$. The set $P_{c}$ of all Pareto critical points is called the Pareto critical set.

When solving an MOP, an initial step can be to compute the Pareto critical set. This set possesses additional structure which can be exploited in numerical schemes. Introducing the function

[TABLE]

we see that Pareto critical points and their corresponding KKT vectors can be described as the zero level set of $F$. As shown by Hillermeier [14], this has the following implication:

Theorem 2.4**.**

Let $J$ be twice continuously differentiable.

• (a)

Let $\mathcal{M}:=\{(u,\alpha)\in\mathbb{R}^{n}\times(\mathbb{R}^{>0})^{k}:F(u,\alpha)=0\}$. If the Jacobian of $F$ has full rank everywhere, i.e.,

[TABLE]

then $\mathcal{M}$ is a $(k-1)$-dimensional differentiable submanifold of $\mathbb{R}^{n+k}$. The tangent space of $\mathcal{M}$ at $(u,\alpha)$ is given by

[TABLE]

• (b)

Let $(u,\alpha)\in\mathcal{M}$ such that (1) holds in $(u,\alpha)$. Then there is an open set $U\subseteq\mathbb{R}^{n}\times\mathbb{R}^{k}$ with $(u,\alpha)\in U$ such that $\mathcal{M}\cap U$ is a manifold as in (a). In other words, $\mathcal{M}$ locally possesses a manifold structure in all points satisfying (1).

Theorem 2.4 forms the basis for the continuation method we use in this article.

2.2 Continuation method with exact gradients

We only give a brief description of the method here and refer to [28] and [14] for details. By Theorem 2.4, the Pareto critical set is – except for the boundary – the projection of the differentiable manifold $\mathcal{M}\subseteq\mathbb{R}^{n}\times\mathbb{R}^{k}$ onto its first $n$ components. In [10] it has been shown that generically, this also holds for the first-order approximations, i.e., the projection of the tangent space of $\mathcal{M}$ yields the tangent cone of $P_{c}$. Given a Pareto critical point $\bar{u}\in P_{c}$, this means that we can find first-order candidates for new Pareto critical points in the vicinity of $\bar{u}$ by moving in the projected tangent space of $\mathcal{M}$. The idea of the continuation method is to do this iteratively to explore the entire Pareto critical set.

Instead of approximating $P_{c}$ by a set of points, we use a set-oriented numerical approach; cf. [28] for details. This has the key advantage that it is easy to check whether a certain part of the set has already been computed, which is difficult when working with points. Additionally, a covering of $P_{c}$ by boxes makes it easy to obtain (and exploit) its topological properties. In the approach, we evenly divide the variable space $\mathbb{R}^{n}$ into hypercubes or boxes $B$ with radius $r>0$:

[TABLE]

Remark 2.5**.**

For ease of notation and readability, we will only consider the case where points $u\in\mathbb{R}^{n}$ are contained in single boxes. In other words, we only consider the case where $u$ is in the interior of a box and not in the intersection of multiple boxes. Since this is the generic case, this has no impact on the numerical methods we will propose later. $\blacksquare$

For $u\in\mathbb{R}^{n}$ let $B(u,r)$ be the box containing $u$. We want to compute the subset of $\mathcal{B}(r)$ covering the Pareto critical set for a given radius $r$, i.e.,

[TABLE]

Since we are interested in a covering via boxes instead of an approximation via points, when moving in a tangent direction of the critical set, we will search for tangent boxes instead of single points. For $u\in\mathbb{R}^{n}$ let

[TABLE]

be the set of neighboring boxes of $B(u,r)$. Starting from a box $B(\bar{u},r)$ containing a critical point $\bar{u}$ with KKT vector $\bar{\alpha}$, we want to explore the neighboring boxes covering the projected tangent space at $\bar{u}$, i.e.,

[TABLE]

Here, $pr_{u}:\mathbb{R}^{n+k}\rightarrow\mathbb{R}^{n}$ is the projection of the tangent space onto the first $n$ components, i.e., the variable space. The typical situation is visualized in Figure 1.

As the tangent space of the Pareto critical set is only a linear approximation, a corrector step is required to verify that a given tangent box actually contains part of the Pareto critical set. This means that there has to be at least one $u\in B$ satisfying (KKT). To this end, for a box $B$, we consider the problem

[TABLE]

Let $\theta(B)$ be the optimal value of this problem. Then obviously

[TABLE]

In particular, if $\theta(B)=0$ and $(\bar{u},\bar{\alpha})$ is the solution of (PC-Box), then $\bar{u}$ is Pareto critical with corresponding KKT vector $\bar{\alpha}$. After solving (PC-Box) in each tangent box, all boxes with $\theta(B)=0$ are added to a queue and a new iteration of the method is started with the first element in the queue. The method stops when the queue is empty, i.e., when there is no neighboring box of the current set of boxes that contains part of the Pareto critical set. For the remainder of this article, we will refer to this method as the exact continuation method.

2.3 Continuation method with inexact gradients

Using ROM to solve the state equation of an MOP of an elliptic PDE will introduce an error in the objective functions and the corresponding gradients, which has to be taken into account in order to ensure Pareto criticality of the solution. We here present a method that calculates a tight superset of the Pareto critical set via numerical continuation, using upper bounds for the errors in the approximated gradients. Formally, we now assume that for each gradient $\nabla J_{i}$, we only have an approximation $\nabla J_{i}^{r}$ such that

[TABLE]

with upper bounds $\epsilon=(\epsilon_{1},...,\epsilon_{k})^{\top}\in\mathbb{R}^{k}$. Let $P_{c}$ and $P^{r}_{c}$ be the Pareto critical sets corresponding to $(\nabla J_{i})_{i}$ and $(\nabla J_{i}^{r})_{i}$, respectively. The following lemma shows how these error bounds translate to error bounds for the KKT conditions:

Lemma 2.6**.**

Let $\bar{u}\in\mathbb{R}^{n}$ be Pareto critical for $J$ with KKT vector $\bar{\alpha}\in\Delta_{k}$. Then

[TABLE]

Proof.

From the estimate

[TABLE]

we derive the claim. ∎

Remark 2.7**.**

Lemma 2.6 can be generalized to equality and inequality constrained MOPs using the constrained version of the optimality conditions from [14]. In this case, in the norm on the left-hand side of the inequality in Lemma 2.6, one additionally has to add a linear combination of the gradients of the equality and inequality constraints. $\blacksquare$

Lemma 2.6 shows that we have to weaken the conditions for Pareto criticality of the reduced objective function to obtain a superset of the actual Pareto critical set $P_{c}$. Formally, let

[TABLE]

$P_{1}^{r}$ was also considered in [21] in the context of descent directions, where the solution of $\min_{\alpha\in\Delta_{k}}\|DJ^{r}(u)^{\top}\alpha\|_{2}^{2}$ is the squared length of the steepest descent direction in $u$. The condition for a point being in $P_{1}^{r}$ only depends on the maximal error $\|\epsilon\|_{\infty}$ and can be seen as a relaxed version of the KKT conditions for the inexact objective function. In contrast to this, the condition in $P_{2}^{r}$ actually considers the individual error bounds. By Lemma 2.6,

[TABLE]

i.e., both $P_{1}^{r}$ and $P_{2}^{r}$ are supersets of $P_{c}$ and $P_{c}^{r}$ (the points $\bar{u}$ for which the inexact gradients satisfy (KKT)). In fact, $P_{2}^{r}$ is a tight superset of $P_{c}$ in the following sense:

Lemma 2.8**.**

Let $\tilde{u}\in P_{2}^{r}$. Then there is some continuously differentiable $\tilde{J}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{k}$ with

[TABLE]

such that $\tilde{u}$ is Pareto critical for $\tilde{J}$.

Proof.

Let

[TABLE]

Since $\tilde{u}\in P_{2}^{r}$ by assumption, we have $\|\nu\|_{2}\leq\tilde{\alpha}^{\top}\epsilon$. Thus

[TABLE]

and

[TABLE]

which proves the lemma. ∎

Lemma 2.8 shows that for each point $\tilde{u}$ in $P_{2}^{r}$, there is an objective function satisfying the error bounds (4) for which $\tilde{u}$ is Pareto critical. As a result, $P_{2}^{r}$ is the tightest superset of $P_{c}$ we can hope for if we only have the estimates in (4). The following example shows both supersets for a simple MOP (cf. [21]).

Example 2.9**.**

Let

[TABLE]

We consider the two error bounds $\epsilon^{1}=(0.2,0.05)^{\top}$ and $\epsilon^{2}=(0,0.2)^{\top}$. The corresponding supersets $P_{1}^{r}$ and $P_{2}^{r}$ are shown in Figure 2.

As $\|\epsilon_{1}\|_{\infty}=\|\epsilon_{2}\|_{\infty}=0.2$, $P_{1}^{r}$ is identical for both error bounds. Considering each component of $J^{r}$ individually, the critical points of $J_{1}^{r}$ and $J_{2}^{r}$ are located at $u^{1}=(1,1)^{\top}$ and $u^{2}=(-1,-1)^{\top}$, respectively. For $P_{2}^{r}$, we see that the difference between $P_{c}^{r}$ and $P_{2}^{r}$ becomes smaller the closer we get to the critical point of the objective function with the smaller error bound. This can be expected, as the influence (or weight) of $\nabla J^{r}_{i}(u)$ in the KKT conditions (KKT) becomes larger the closer $u$ is to $u^{i}$. In particular, in Figure 2(b), the difference between $P_{2}^{r}$ and $P_{c}^{r}$ at $(1,1)^{\top}$ becomes zero, as $\epsilon^{2}_{1}=0$. $\Diamond$

If we set $\epsilon_{i}=\|\epsilon\|_{\infty}$ for all $i\in\{1,...,k\}$, then $P_{1}^{r}=P_{2}^{r}$. Thus, we will from now on only consider $P_{2}^{r}$. As shown in the previous example, the “dimension” of $P_{2}^{r}$ is higher than the “dimension” of $P_{c}^{r}$. More precisely, $P_{2}^{r}$ contains the closure of an open subset of $\mathbb{R}^{n}$, which is shown in the following lemma:

Lemma 2.10**.**

Let $\nabla J^{r}_{i}$ be continuous for all $i\in\{1,...,k\}$. Let

[TABLE]

Then

1. (a)

$P_{2}^{r}$* is closed. In particular, $\overline{A}\subseteq P_{2}^{r}$.* 2. (b)

$A$* is open.*

Proof.

(a) The case $P_{2}^{r}=\emptyset$ is trivial, so we assume that $P_{2}^{r}\neq\emptyset$. Let $\bar{u}\in\overline{P_{2}^{r}}$. Then there is a sequence $(u^{i})_{i}\in P_{2}^{r}$ with $\lim_{i\rightarrow\infty}u^{i}=\bar{u}$. Consider the sequence $(\alpha^{i})_{i}\in\Delta_{k}$ with

[TABLE]

By compactness of $\Delta_{k}$, we can assume w.l.o.g. that there is some $\bar{\alpha}\in\Delta_{k}$ with $\lim_{i\rightarrow\infty}\alpha^{i}=\bar{\alpha}$. Let

[TABLE]

By our assumption, $\Psi$ is continuous. From $\Psi(u^{i},\alpha^{i})<0$ for all $i\in\mathbb{N}$ it follows that $\Psi(\bar{u},\bar{\alpha})\leq 0$, which yields $\bar{u}\in P_{2}^{r}$.

(b) The case $A=\emptyset$ is again trivial such that we assume $A\neq\emptyset$. Let $\bar{u}\in A$ with

[TABLE]

Let $\psi:\mathbb{R}^{n}\rightarrow\mathbb{R}$, $u\mapsto\|DJ^{r}(u)^{\top}\bar{\alpha}\|_{2}^{2}-(\bar{\alpha}^{\top}\epsilon)^{2}$. Then $\psi(\bar{u})<0$ and by our assumption, $\psi$ is continuous. Therefore, there is some open set $U\subseteq\mathbb{R}^{n}$ with $\bar{u}\in U$ such that $\psi(u)<0$ for all $u\in U$. Since

[TABLE]

we have $U\subseteq A$ such that $A$ is open. ∎

We will now present two strategies for the numerical computation of $P^{r}_{2}$. Analogously to the case with exact gradients, we will approximate $P_{2}^{r}$ via the box covering

[TABLE]

2.3.1 Strategy 1

The idea of our first method is to mimic the exact continuation method to calculate $\mathcal{B}_{c}^{r}$. For this, there are mainly two modifications we have to make:

By Lemma 2.10, $P_{2}^{r}$ is not a lower-dimensional object in $\mathbb{R}^{n}$, so it makes no sense to use tangent information to find first-order candidates as in (3). Instead, we have to consider all neighboring boxes. 2. 2.

The problem (PC-Box) has to be replaced by a problem that checks the defining inequality of $P_{2}^{r}$.

As a replacement for (PC-Box), we consider the following problem:

[TABLE]

Let $\theta_{\epsilon}(B)$ be the optimal value of this problem. Note that $\theta_{\epsilon}(B)<0$ is sufficient to verify that a box $B$ contains part of $P_{2}^{r}$. As a result, we do not need to solve ($\epsilon$PC-Box) exactly. For example, when using an iterative method for the solution of ($\epsilon$PC-Box), we can stop when the function value is negative. The above mentioned changes yield Algorithm 1.

Due to the loss of low-dimensionality of $P_{2}^{r}$, the formulation of the continuation method becomes much simpler. As a consequence, it is straightforward to show that Algorithm 1 yields the desired covering $\mathcal{B}_{c}^{r}(r)$.

When executing the exact continuation method directly using inexact gradients (i.e., forgetting about the inexactness) and comparing it to Algorithm 1 (with the same box radius), the former will generally be much faster than the latter. A suitable way to evaluate the run time is to compare the number of times Problems (PC-Box) and ($\epsilon$PC-Box) need to be solved, respectively, as they require the majority of the computing time and are equally difficult to solve. (Here, we assume that both problems are solved with equal precision.) For each box added to the collection $\mathcal{B}$ in either algorithm, one of these problems has to be solved. Consequently, the longer run time of Algorithm 1 is partly due to the fact that $P_{2}^{r}$ is a superset of $P_{c}^{r}$, which means that more boxes are required to cover $P_{2}^{r}$ than $P_{c}^{r}$. However, even if the error bounds $\epsilon$ are small such that $P_{2}^{r}$ and $P_{c}^{r}$ are almost equal, Algorithm 1 will be slower. This is due to the fact that instead of only the tangent boxes, all neighboring boxes have to be tested with ($\epsilon$PC-Box) in each loop of Algorithm 1. While this does not matter in the interior of $P_{2}^{r}$ (as all neighboring boxes are in fact in $P_{2}^{r}$ in that case), it is very inefficient at the boundary of $P_{2}^{r}$. This is the motivation for the second strategy.

2.3.2 Strategy 2

By Lemma 2.10, $P_{2}^{r}$ has the same dimension as the space of variables $\mathbb{R}^{n}$. This means that it can be described much more efficiently by its topological boundary $\partial P_{2}^{r}$. To be more precise, $\mathbb{R}^{n}\setminus\partial P_{2}^{r}$ consists of different connected components that lie either completely inside or completely outside $P_{2}^{r}$. So if we know $\partial P_{2}^{r}$, we merely have to test one point of each connected component if it is contained in $P_{2}^{r}$ or not to completely determine $P_{2}^{r}$. Therefore, the idea of our second strategy is to only compute $\partial P_{2}^{r}$.

Let

[TABLE]

This map is well-defined since $\Delta_{k}$ is compact, i.e., the minimum always exists. By Lemma 2.10, we have $\partial P_{2}^{r}\subseteq\varphi^{-1}(0)$. Our goal is to compute $\varphi^{-1}(0)$ via a continuation approach. To this end, we first have to show that $\varphi$ is differentiable. We will do this by investigating the properties of the optimization problem in (5), i.e., of the problem

[TABLE]

for

[TABLE]

This leads to the following result.

Theorem 2.11**.**

Let $\bar{u}\in\varphi^{-1}(0)$ such that (2.3.2) has a unique solution $\bar{\alpha}\in\Delta_{k}$ with $\bar{\alpha}_{i}>0$ for all $i\in\{1,...,k\}$. Let (2.3.2) be uniquely solvable in a neighborhood of $\bar{u}$. Then there is an open set $U\subseteq\mathbb{R}^{n}$ with $\bar{u}\in U$ such that $\varphi|_{U}$ is continuously differentiable.

Proof.

See Appendix A. ∎

For a standard continuation approach, we also have to show that $\varphi^{-1}(0)$ is a manifold. By the Level Set Theorem (cf. [19], Corollary 5.14), to properly show that $\varphi^{-1}(0)$ is a manifold in a neighborhood of some $\bar{u}\in\varphi^{-1}(0)$, we would have to show that $D\varphi|_{U}(\bar{u})\neq 0$ (cf. (A.4)). From the theoretical point of view, this poses a problem as there is no obvious way to achieve this. In practice however, we can test this by checking if the norm of $D\varphi|_{U}(\bar{u})$ is below a certain threshold. If this is the case, and if $\varphi^{-1}(0)$ is indeed not a manifold, we again have to consider all neighboring boxes as tangent boxes as in strategy 1. Otherwise, if $D\varphi|_{U}(\bar{u})\neq 0$, we can compute the tangent space $T_{\bar{u}}$ of $\varphi^{-1}(0)$ at $\bar{u}$ via

[TABLE]

Finally, in analogy to (PC-Box) and ($\epsilon$PC-Box), we will use the following problem to test if a box $B$ contains part of $\partial P_{2}^{r}$:

[TABLE]

The resulting continuation method is presented in Algorithm 2.

Remark 2.12**.**

Since for every evaluation of $\varphi$ the solution of the quadratic problem (2.3.2) has to be computed, ($\partial\epsilon$PC-Box) is significantly more difficult to solve than ($\epsilon$PC-Box). Additionally, we are looking for the points $u$ where $\varphi(u)=0$, i.e., where the problem (2.3.2) is not positive definite. This increased difficulty of Strategy 2 is compensated by the fact that far fewer boxes have to be checked with ($\partial\epsilon$PC-Box) than with ($\epsilon$PC-Box) in Strategy 1. 2. 2.

When all $\epsilon_{i}=\bar{\epsilon}$ are equal, $\varphi(u)=-\bar{\epsilon}^{2}$ for all $u\in P_{c}^{r}$, i.e., $\varphi$ is constant on the Pareto critical set $P_{c}^{r}$. This means that local solvers may fail to find a minimum of $\varphi$ when the box $B$ in ($\partial\epsilon$PC-Box) has a nonempty intersection with $P^{r}$. An obvious but expensive way to circumvent this problem is to start the local solver multiple times with different initial points. Alternatively, one can use sufficient conditions for a box $B$ containing part of $\varphi^{-1}(0)$ before actually solving ($\partial\epsilon$PC-Box). For example, by the intermediate value theorem, if there are two points in $B$ where $\varphi$ has different sign, we immediately know that $\varphi(u)=0$ for some $u\in B$. (But note that for this method, we still need to find a point in $\varphi^{-1}(0)\cap B$ to be able to calculate the tangent space of $\varphi^{-1}(0)$). 3. 3.

In practice, error bounds which are zero can cause problems for the stability of Strategy 2. For example, in Figure 2(b), the width of $P_{2}^{r}$ becomes arbitrarily small near $(1,1)^{\top}$. As a result, Strategy 2 may jump between different parts of the boundary and thus miss certain parts. Additionally, since the boundary of $P_{2}^{r}$ typically intersects the Pareto critical set $P_{c}$ in this case, ($\partial\epsilon$PC-Box) may be difficult to solve (as in 2.). Thus, in practice, one should use error bounds that are slightly larger than zero, even if the corresponding gradients are exact.

$\blacksquare$**

2.4 Globalization approach

Note that all algorithms presented in this section so far approximate either $P_{c}$, $P_{2}^{r}$ or $\partial P_{2}^{r}$ by starting in an initial point $u_{0}$ and then locally exploring in all (tangent) directions. Thus, if the set we want to approximate is disconnected, we can only compute the connected component that contains $u_{0}$. In the following, we will describe how we can solve this problem, i.e., how our methods can be globalized.

As mentioned earlier, an advantage of using boxes in the continuation method instead of points is the fact that it is easy to detect whether a region has already been explored. In particular, this allows us to start the continuation in multiple initial points at the same time, by simply adding all of them to the queue $Q$ in step 1 of Algorithms 1 or 2 (and initializing the covering $\mathcal{B}$ with the corresponding boxes). As a result, to globalize our methods, we merely have to find an initial set $U_{0}$ of points such that the intersection of $U_{0}$ with each connected component is nonempty.

For obtaining an initial set, we make use of the optimization problems that verify if a box contains part of the set we want to approximate, i.e., the problems (PC-Box), ($\epsilon$PC-Box) and ($\partial\epsilon$PC-Box). The idea is to consider a box covering as in (2) with large radius $R$ and then simply test each box for relevant points using these problems. Let $B_{0}$ be a compact superset of the set that we want to approximate (i.e., of $P_{c}$, $P_{2}^{r}$ or $\partial P_{2}^{r}$), e.g., a large outer box. For ease of notation, we assume that $B_{0}$ is a union of boxes in $\mathcal{B}(R)$. For the case of the Pareto critical set $P_{c}$, i.e., the globalization of the exact continuation method, the resulting method is presented in Algorithm 3. The corresponding globalization methods for Algorithm 1 and 2 are obtained by replacing (PC-Box) in step 3 by ($\epsilon$PC-Box) and ($\partial\epsilon$PC-Box), respectively.

The radius $R$ has to be chosen such that for each connected component, there is at least one box in our covering that only has an intersection with the desired component. In theory, $R$ can obviously become very small if two different connected components are very close to each other. In this case, Algorithm 3 becomes infeasible to use, as the number of boxes that have to be tested becomes too large. In practice however, the components are often sufficiently far apart such that a large radius is sufficient and only few boxes have to be considered.

For the globalization of the exact continuation method and Algorithm 1, we only have to take the non-connectivity of $P_{c}$ and $P_{2}^{r}$ into account. For Algorithm 2, an additional problem may arise since the boundary $\partial P_{2}^{r}$ does not necessarily need to be smooth. Non-smoothness of $\partial P_{2}^{r}$ is caused by points in which $\varphi$ is not differentiable. (By Theorem 2.11, these are points where the solution of ($\partial\epsilon$PC-Box) is not unique.) In these points, $\partial P_{2}^{r}$ does not posses a tangent space, and our method will be unable to continue. As a result, we have to ensure in the initialization of Algorithm 2 that we choose an initial point in $U_{0}$ on each smooth component of $\partial P_{2}^{r}$. Visually, these can be thought of as the faces of $P_{2}^{r}$.

We conclude this section with some remarks on the practical use of Algorithm 3.

Remark 2.13**.**

For MOPs with a high-dimensional variable space, Algorithm 3 quickly becomes infeasible due to the exponential growth of the number of boxes in $\mathcal{B}(R)$. For these cases, an initialization based on points instead of boxes should be used, for example by applying methods from global optimization to modified versions of (PC-Box), ($\epsilon$PC-Box) and ($\partial\epsilon$PC-Box), where $u$ is not constrained to a box $B$. 2. 2.

Instead of directly looping over all boxes in step 2 of Algorithm 3, in some cases it might be more beneficial to first execute a few steps of the subdivision algorithm (cf. **[8]**) to quickly discard boxes that are far away from the Pareto critical set.

$\blacksquare$**

3 Multiobjective optimization of an elliptic PDE using the RB method

In this section we will present a multiobjective (parameter) optimization problem of an elliptic advection-diffusion-reaction equation and show how the reduced basis method can be applied in view of the continuation method for inexact gradients from Section 2.3 (see Algorithms 1 and 2).

3.1 Multiobjective optimization of an elliptic PDE

Given a domain $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$, we consider the problem

[TABLE]

s.t.

[TABLE]

and the bilateral box constraints

[TABLE]

where $u=(u_{1},\ldots,u_{m})=(\kappa_{1},\ldots,\kappa_{m^{\prime}},c,r)\in\mathbb{R}^{m}$ is the parameter of dimension $m:=m^{\prime}+2$, $U_{\textsl{ad}}:=\{u\in\mathbb{R}^{m}\mid u_{a}\leq u\leq u_{b}\}$ is the admissible parameter set, and $y\in L^{2}(\Omega)=:H$ is the state variable.

The domain $\Omega$ is divided into $m^{\prime}$ pairwise disjoint subdomains $\Omega=\Omega_{1}\dot{\cup}\ldots\dot{\cup}\,\Omega_{m^{\prime}}$, such that $\kappa_{i}$ is the diffusion coefficient on $\Omega_{i}$. The vector field $b\in L^{\infty}(\Omega,\mathbb{R}^{d})$ is the given advection, whose strength and orientation can be controlled by the parameter $c\in\mathbb{R}$. Moreover, the reaction coefficient is given by the parameter $r>0$, and $f\in H$ is the inhomogeneity on the right-hand side of the equation. On the boundary we impose homogeneous Neumann boundary conditions.

The cost functions $\mathcal{J}_{1},\ldots,\mathcal{J}_{k-1}:H\times U_{\textsl{ad}}\to\mathbb{R}^{k}$ are of tracking type with respect to the desired states $y_{d}^{1},\ldots,y_{d}^{k-1}\in H$, and the cost function $\mathcal{J}_{k}:H\times U_{\textsl{ad}}\to\mathbb{R}^{k}$ measures the parameter cost.

Setting $V:=H^{1}(\Omega)$ and using the parameter-dependent bilinear form $a(u;\cdot,\cdot):V\times V\to\mathbb{R}$ defined by

[TABLE]

for all $u\in U_{\textsl{ad}}$ and $\varphi,\psi\in V$, and the linear functional $F:V\to\mathbb{R}$ given by $F(\varphi):=\langle f,\varphi\rangle_{H}$ for all $\varphi\in V$, we can write (EPDE) in its weak formulation as: Find $y\in V$ such that

[TABLE]

is satisfied. It is possible to show the unique solvability of (13) under some conditions on the parameter $u$.

Theorem 3.1**.**

There are $\kappa_{\textsl{min}}\in(0,\infty)^{m^{\prime}}$, $c_{\textsl{min}},c_{\textsl{max}}\in\mathbb{R}$ with $c_{\textsl{min}}<c_{\textsl{max}}$ and $r_{\textsl{min}}\in(0,\infty)$ such that (13) has a unique solution $y(u)\in V$ for every parameter $u=\left(\kappa,c,r\right)\in\mathbb{R}^{m}$ with $\kappa>\kappa_{\textsl{min}}$, $c_{\textsl{min}}<c<c_{\textsl{max}}$ and $r>r_{\textsl{min}}$.

Proof.

It is straightforward to show that for all parameters $u\in\mathbb{R}^{m}$ the bilinear form $a(u;\cdot,\cdot)$ and the linear functional $F$ are continuous, and that there are $\kappa_{\textsl{min}}\in(0,\infty)^{m^{\prime}}$, $c_{\textsl{min}},c_{\textsl{max}}\in\mathbb{R}$ with $c_{\textsl{min}}<c_{\textsl{max}}$ and $r_{\textsl{min}}\in(0,\infty)$ such that $a(u;\cdot,\cdot)$ is coercive for all $u=\left(\kappa,c,r\right)\in\mathbb{R}^{m}$ with $\kappa>\kappa_{\textsl{min}}$, $c_{\textsl{min}}<c<c_{\textsl{max}}$ and $r>r_{\textsl{min}}$. Now the Lax-Milgram Theorem can be applied to show the unique solvability of (13). ∎

With Theorem 3.1 in mind we can introduce the solution operator of the elliptic PDE.

Definition 3.2**.**

Define the set $U_{\textsl{eq}}:=(\kappa_{\textsl{min}},\infty)\times(c_{\textsl{min}},c_{\textsl{max}})\times(r_{\textsl{min}},\infty)$ with the constants from Theorem 3.1. Let $\mathcal{S}:U_{\textsl{eq}}\to V\hookrightarrow H$ be defined as the solution operator of (13), i.e., the function $y:=\mathcal{S}(u)$ solves the weak formulation (13) for any parameter $u\in U_{\textsl{eq}}$.

Remark 3.3**.**

In the following we suppose that it holds $U_{\textsl{ad}}\subset U_{\textsl{eq}}$. $\blacksquare$

Using the explicit dependence of the state $y$ on the parameter $u$ for all $u\in U_{\textsl{ad}}$, the essential cost functions $J_{1},\ldots,J_{k}:U_{\textsl{ad}}\to\mathbb{R}$ can be defined.

Definition 3.4**.**

For any $i\in\{1,\ldots,k\}$ let the essential cost function $J_{i}:U_{\textsl{ad}}\to\mathbb{R}$ be given by $J_{i}(u):=\mathcal{J}_{i}(\mathcal{S}(u),u)$ for all $u\in U_{\textsl{ad}}$.

For applying the continuation method from Section 2, which is based on Theorem 2.4, to solve this multiobjective parameter optimization problem, the cost functions $J_{1},\ldots,J_{k}$ need to be twice continuously differentiable. This is the statement of the next lemma.

Lemma 3.5**.**

The cost functions $J_{1},\ldots,J_{k}$ are twice continuously differentiable.

Proof.

It is clear that the cost function $J_{k}$ is twice continuously differentiable. Furthermore, it is possible to show that the solution operator $\mathcal{S}$ of (13) is twice continuously differentiable (this can be shown by rewriting (13) in the form $e(y,u)=0$ and then using the implicit function theorem, cf. [15, Section 1.6]). From this it immediately follows that the cost functions $J_{1},\ldots,J_{k-1}$ are twice continuously differentiable as well. ∎

For later use, we need an explicit formula for the gradients $\nabla J_{1},\ldots,\nabla J_{k}$. Therefore, we introduce the so-called adjoint equation for all $i\in\{1,\ldots,k-1\}$: Find $p\in V$ such that it holds

[TABLE]

With the same arguments as in Theorem 3.1 it is possible to show that (14) has a unique solution for all $u\in U_{\textsl{eq}}$.

Definition 3.6**.**

Denote by $\mathcal{A}_{i}:U_{\textsl{eq}}\to V\hookrightarrow H$ the solution operator of the adjoint equation (14) for all $i\in\{1,\ldots,k-1\}$.

Now a small computation shows that

[TABLE]

which yields

[TABLE]

for all $i\in\{1,\ldots,k-1\}$. Lastly, it is obvious that $\nabla J_{k}(u)=u$.

3.2 The reduced basis method

For computing the Pareto critical set of the problem (MPOP) by the exact continuation method introduced in Section 2.2, the problem (PC-Box) has to be solved numerous times. However, already one gradient evaluation of all cost functions $\nabla J_{1}(u),\ldots,\nabla J_{k}(u)$ involves the solution of one state and $k-1$ adjoint equations. Thus, using a finite element discretization for the weak formulations (13) and (14), which leads to large linear equation systems, is numerically very costly and time consuming. Therefore, the use of reduced-order modelling (ROM) is a common tool to lower the computational costs.

The idea of ROM is to use a low-dimensional subspace $V^{r}\subset V$ as a surrogate for the infinite-dimensional space $V$ in the weak formulations (13) and (14). Given a finite-dimensional reduced-order space $V^{r}\subset V$, the reduced-order state equation reads: Find $y^{r}\in V^{r}$ such that

[TABLE]

is satisfied.

With the same arguments as in Theorem 3.1 it can be shown that (22) has a unique solution for all $u\in U_{\textsl{eq}}$. Therefore, we can follow the procedure of Section 3.1 and introduce the solution operator $\mathcal{S}^{r}:U_{\textsl{eq}}\to V^{r}\subset V\hookrightarrow H$ of the ROM state equation (22) and consequently the ROM essential cost functions $J_{1}^{r},\ldots,J_{k}^{r}$, which are defined by $J_{i}^{r}(u):=\mathcal{J}_{i}(\mathcal{S}^{r}(u),u)$ for all $u\in U_{\textsl{ad}}$ and all $i\in\{1,\ldots,k\}$. Again, it can be shown that the functions $J_{1}^{r},\ldots,J_{k}^{r}$ are twice continuously differentiable so that they fit into the framework of Theorem 2.4. The gradient of the cost functions can also be displayed by the reduced-order adjoint equations

[TABLE]

for all $i\in\{1,\ldots,k-1\}$, whose solution operator we denote by $\mathcal{A}_{i}^{r}:U_{\textsl{eq}}\to V^{r}\subset V\hookrightarrow H$. With this definition it holds

[TABLE]

for all $i\in\{1,\ldots,k-1\}$. Moreover, we have $\nabla J_{k}^{r}(u)=u=\nabla J_{k}(u)$.

In this paper we use a particular model-order reduction technique, namely the reduced basis (RB) method (see e.g. [26, 13, 25]). In the RB method the snapshot space $V^{r}$ is spanned by solutions of the state equation and the adjoint equations to different parameter values $u\in U_{\textsl{ad}}$. The reduced basis is then given by an orthonormal basis $(\Phi_{1},\ldots,\Phi_{N})$ of the space $V^{r}$.

By using the RB method we introduce an error in the state equation, which transfers to the cost functions, its gradients and eventually to the Pareto critical set, which we want to compute. In Section 2.3 two strategies were presented to deal with the inflicted inexactness in the gradients of the multiobjective optimization problem. Both are based on the estimates (4) for the errors in the gradients of the cost functions. Thus, when applying the RB method we need to ensure these estimates. This is done by using the well-known greedy algorithm (cf. [5]). Given a sufficiently fine finite parameter training set $\mathcal{P}\subset U_{\textsl{ad}}$ new solution snapshots are computed until the error in the gradients of all cost functions is smaller than the predefined error tolerance for all parameters in $\mathcal{P}$. The parameter for the new snapshots is thereby chosen as the one for which the error in the gradient is the largest. The procedure is summarized in Algorithm 4.

3.3 Error estimation for the gradients

In the greedy procedure in Algorithm 4, the error between the full-order and the reduced-order gradients has to be evaluated. There are two strategies to do so.

The full-order gradients are computed and stored at the beginning of the greedy procedure. Therefore, in each greedy iteration, only the reduced-order gradients have to be computed and the error can be easily evaluated. Of course, this implies large computational costs at the beginning of the greedy procedure. This method is called strong greedy algorithm (cf. [12, 5]). 2. 2.

An a-posteriori error estimator for the errors in the gradient is used, which can be efficiently evaluated. This results in computational costs for the greedy algorithm, which only depend on the reduced-order dimension $N$.

To be able to follow the second strategy we introduce a rigorous a-posteriori error estimator for the error in the gradient of the cost functions.

Using the gradient representations (21) and (30) we can write for $i\in\{1,\ldots,k-1\}$

[TABLE]

Due to the bilinearity and the continuity of $a_{1},\ldots,a_{m}$ and the triangle inequality, we can further write

[TABLE]

for all $j\in\{1,\ldots,m\}$.

Therefore, we need a-posteriori error estimators for the state and the adjoint equations in order to be able to estimate the approximation error induced in the gradients. To this end, we use the following well-known estimators (cf. [26]).

[TABLE]

where the residuals $r_{\mathcal{S}}(u)$ and $r_{\mathcal{A}_{i}}(u)$ are given by

[TABLE]

For methods on how to estimate $\alpha(u)$ and to evaluate the terms $\left\|r_{\mathcal{S}}(u)\right\|_{V^{\prime}}$ and $\left\|r_{\mathcal{A}_{i}}(u)\right\|_{V^{\prime}}$ efficiently, we refer for example to [26].

Remark 3.7**.**

Since $J_{k}=J_{k}^{r}$, the gradients of the two functions also coincide, so that the $\nabla J_{k}$ is approximated exactly by $\nabla J_{k}^{r}$. $\blacksquare$

4 Numerical results

In this section we will numerically investigate the application of the continuation method presented in Section 2 to the PDE-constrained multiobjective optimization problem using the reduced basis method in Section 3.

For the discretization of the state and adjoint equations we used linear finite elements with 714 degrees of freedom.

4.1 Generation of the reduced basis

For investigating the generation of the reduced basis by the greedy algorithm in Algorithm 4, we consider the MPOP

[TABLE]

with $u=(\kappa_{1},\kappa_{2},c,r)$, $\Omega_{1}=(0,1)\times(0,0.5)$, $\Omega_{2}=(0,1)\times(0.5,1)$, and the admissible parameter set

[TABLE]

The reason for setting $c=0$ in this example is that the coercivity constant $\alpha(u)$ of the bilinear form $a(u;\cdot,\cdot)$ is explicitly given by $\alpha(u)=\min\{\kappa_{1},\kappa_{2},r\}$ for all $u\in U_{ad}$, so that we expect a good efficiency of the error estimator of both the state and adjoint equations.

This is verified by the results shown in Figure 3 (a), where the efficiency of the error estimator for both equations is shown for a given reduced basis for 1000 randomly chosen parameter values. However, the resulting efficiency of the error estimator for the error in the gradient is between $10^{3}$ and $10^{6}$ (see Figure 3 (b)) and thus not well suited for a greedy procedure, which depends on a good error estimation. The huge overestimation of the error estimator is mainly due to the use of the triangle inequality (31) and the continuity estimates (32), as can be seen in Figure 3 (b).

Compared to the strong greedy algorithm, we can see in Table 1 that this overestimation results in far more basis elements than actually needed to reach the given error bound. Since we want to investigate the influence of the error bounds in the estimate (4) on the problem, we want that the estimate (4) is satisfied sharply by the RB. Therefore, we will not use the error estimator to generate the basis, but instead use the strong greedy algorithm.

4.2 Application of the continuation methods to an MPOP

For the numerical investigation of the continuation method applied to a PDE-constrained multiobjective parameter optimization problem together with the use of the reduced basis method, we consider the MPOP

[TABLE]

with $u=(\kappa,c,r)$ and

[TABLE]

i.e., the reaction parameter $r$ is a constant so that we only optimize the diffusivity in the whole domain $\Omega$ and the strength and orientation of the advection field $b$. Thus, this can be seen as a problem with two parameters.

As described before, the reduced basis is generated by the strong greedy Algorithm 4, where the error bounds $\epsilon_{1},\ldots,\epsilon_{4}$ are chosen in accordance with the estimate (4). As a reference, the exact solution of (46) (via exact continuation and FEM discretization of the weak formulations) is shown in Figure 4.

Remark 4.1**.**

Since (46) is constrained to a box, we have to use a constrained version of the exact continuation method (cf. [14]) to calculate Pareto critical points that lie on the boundaries of (46). But note that for this example, all Pareto critical points on the boundary are also Pareto critical if we ignore the constraints. In other words, for each Pareto critical point $\bar{u}$ on the boundary, there is a sequence of Pareto critical point in the interior that converges to $\bar{u}$. By continuity of $DJ$, the gradients of the (active) inequality constraints in the KKT conditions can be ignored. As a result, we can treat (46) as an unconstrained problem that we only solve in a certain area. $\blacksquare$

As a first test, we will compare the time needed to compute the exact solution of (46) with the time needed for Strategy 1 and 2. For the error bounds we choose $\epsilon=(0.03,0.03,0.01,0.01)$ and for the box radius we choose $r=\frac{3-0.5}{2^{9}}\approx 0.0049$. The results are shown in Figure 5.

All three methods were implemented in Matlab. For the solution of the subproblems (PC-Box), ($\epsilon$PC-Box) and ($\partial\epsilon$PC-Box), the SQP-Algorithm of fmincon was used. (For increased stability during the continuation, each subproblem where the SQP-Algorithm found an optimal value larger than zero was restarted using the Interior-Point-Method and the Active-Set-Method of fmincon). The runtime, number of boxes and number of subproblems needed are shown in Table 2.

When comparing Strategy 1 and Strategy 2, we see that Strategy 2 needs about $20$ times fewer boxes and solutions of subproblems than Strategy 1. This is to be expected, since Strategy 2 only computes a covering of the boundary of $P_{2}^{r}$, i.e., of a lower dimensional set. When comparing the actual runtime, Strategy 2 is about $5$ times faster than Strategy 1, since the subproblems in Strategy 2 are more expensive to solve than the ones in Strategy 1 (cf. Remark 2.12). Finally, Strategy 2 is about $63$ times faster than the exact continuation method with FEM discretization, illustrating the large increase in efficiency we gain from our approach.

Although it is a lot quicker to use inexact gradients from ROM instead of the exact gradients via FEM, it is important to keep in mind that our methods are computing a superset of the actual Pareto critical set. For example, in Figure 5, the right side of the lower connected component is only approximated poorly by $P_{2}^{r}$. Therefore, we will now investigate the influence of the error bounds $\epsilon=(\epsilon_{1},\epsilon_{2},\epsilon_{3},\epsilon_{4})$ on $P_{2}^{r}$, by applying Strategy 2 with reduced bases for different values of $\epsilon$. Note that in all our tests we set $\epsilon_{4}=0.01$, although the error in the gradient of the fourth cost function is zero for all parameters. This is done to make the solution of ($\partial\epsilon$PC-Box) in line 7 of Algorithm 2 numerically stable (cf. Remark 2.12).

The results of our experiment can be seen in Figure 6. Generally, as expected, the boundary $\partial P_{2}^{r}$ encloses the Pareto critical set $P_{c}$ sharper and sharper for decreasing $\epsilon$. Moreover, we observe that it is crucial to choose an $\epsilon$ which is not too large: For the value $\epsilon=(0.1,0.1,0.1,0.01)$ the shape of the boundary $\partial P_{2}^{r}$ implies that the set $P_{2}^{r}$ is connected, i.e., we lose the topological information that the Pareto critical set actually consists of two connected components. Decreasing $\epsilon$ to $\epsilon=(0.0885,0.0885,0.0885,0.01)$ we are in the limit case in which the boundary $\partial P_{2}^{r}$ touches the box constraints at around $(2.3,1)$, so that this is the approximate $\epsilon$ for which we regain the basic topological information of a disconnected Pareto critical set.

If we compare the results for $\epsilon=(0.03,0.03,0.03,0.01)$, $\epsilon=(0.03,0.03,0.01,0.01)$ and $\epsilon=(0.03,0.01,0.01,0.01)$, the influence of changing one component of $\epsilon$ becomes obvious. For $\epsilon=(0.03,0.03,0.03,0.01)$ the set $\partial P_{2}^{r}$ encloses the set $P_{c}$ quite sharply at the upper connected component and at the left part of the lower connected component, where the second and third component of the corresponding KKT-multipliers $\alpha$ are small. On the other hand, in the right part of the lower connected component of $P_{c}$, where the second and third component of the corresponding KKT-multipliers are relatively large, the deviation of $\partial P_{2}^{r}$ to $P_{c}$ is still large. Consequently, first reducing $\epsilon_{3}$ and then also $\epsilon_{2}$ from 0.03 to 0.01 leads to a clearly visible sharper enclosing of this part of $P_{c}$.

5 Conclusion and outlook

In this article, we present a way to efficiently solve multiobjective parameter optimization problems of elliptic PDEs by combining the reduced basis method from PDE-constrained optimization with the continuation approach from multiobjective optimization, which computes a box covering of the Pareto critical set. Using the RB method in this setting introduces an error in the objective functions and their gradients that has to be considered when solving the MPOP. To this end, we require that the reduced basis guarantees error bounds for the gradients of the objective functions. These error bounds are then incorporated into the KKT optimality conditions for MOPs to derive a tight superset $P_{2}^{r}$ of the actual Pareto (critical) set. This superset can be computed using a straightforward modification of the continuation method for MOPs (Strategy 1). Since $P_{2}^{r}$ has the same dimensions as the variable space of the MOP, we afterwards present a second method that only computes the boundary $\partial P_{2}^{r}$ of $P_{2}^{r}$ (Strategy 2). We do this by showing that $\partial P_{2}^{r}$ can be written as the level set of a differentiable mapping, which again enables the use of a continuation approach to compute it. For constructing the reduced basis, we use a greedy procedure which incorporates, and thus ensures, the error bounds for the gradients of the objective functions.

Our numerical tests show that the presented a-posteriori error estimator for the error in the gradients is not well-suited for the application in a greedy procedure due to its bad efficiency. Therefore, a strong greedy algorithm is used to build the reduced basis. Concerning the solution of the MPOP we investigate two aspects: First, the runtimes of our methods are compared. In our case, Strategy 1 is about 13 times and Strategy 2 about $63$ times faster than the exact solution of the MPOP (via the classical continuation method with FEM discretization). Second, the influence of the error bound for the gradients of the objective functions is investigated. As expected, a smaller error bound leads to a tighter covering of the Pareto critical set. Moreover, we observe that single components of the error bound strongly influence the tightness of the covering in areas, in which the corresponding components of the KKT-multipliers are large. Thus, by individually adapting the single components of the error bound, we can nicely control the tightness of the covering.

For future work, there are some theoretical and practical aspects that should be investigated further:

• •

As mentioned in Remark 2.12, in certain situations there can be difficulties when solving the problem ($\partial\epsilon$PC-Box). In these situations, specialized methods that take these difficulties into account should be developed and used instead of standard methods for constrained optimization.

• •

If the number of objectives of the MPOP is larger than the number of variables, it may be possible to combine our approaches in this article with the hierarchical decomposition of the Pareto critical set presented in [10].

• •

The development of a more efficient a-posteriori error estimator for the error in the gradients of the objective functions would allow to use it in the greedy procedure. In that way, the expensive strong greedy procedure would be avoided in the offline phase. One way to do so might be the application of localized RB methods, see e.g. [20].

• •

As explained in the globalization approach in Section 2.4, we have to use multiple initial points to ensure that we find all connected components of $P_{2}^{r}$ (and faces of $\partial P_{2}^{r}$). Due to the local nature of the continuation method, this approach can potentially be parallelized, increasing the efficiency of our methods even more.

• •

If a decision maker is present with a certain preference, it may be worth to steer our continuation method in a direction that results from that preference instead of approximating the complete Pareto set. For the case with exact gradients, this was done in [27].

Acknowledgment

This research was funded by the DFG Priority Programme 1962 “Non-smooth and Complementarity-based Distributed Parameter Systems”.

Appendix A Proof of Theorem 2.11

To prove Theorem 2.11, we first have to investigate some of the properties of the optimization problem (2.3.2). This problem is quadratic with linear equality and inequality constraints. We will first investigate the uniqueness of the solution in the following lemma.

Lemma A.1**.**

Let $u\in\varphi^{-1}(0)$ and let $\alpha^{1}$ and $\alpha^{2}$ be two solutions of (2.3.2) with $\alpha^{1}\neq\alpha^{2}$. Then $\omega(\alpha)=0$ for all $\alpha\in span(\{\alpha^{1},\alpha^{2}\})$ and

[TABLE]

Proof.

For $c_{1},c_{2}\in\mathbb{R}\setminus\{0\}$ we have

[TABLE]

From $\omega(\alpha^{1})=\omega(\alpha^{2})=0$ it follows that $\epsilon^{\top}\alpha^{1}=\|DJ^{r}(u)^{\top}\alpha^{1}\|$ and $\epsilon^{\top}\alpha^{2}=\|DJ^{r}(u)^{\top}\alpha^{2}\|$. Let $\sphericalangle$ be the angle between $DJ^{r}(u)^{\top}\alpha^{1}$ and $DJ^{r}(u)^{\top}\alpha^{2}$. Then

[TABLE]

Assume $cos(\sphericalangle)\neq 1$ (i.e., $cos(\sphericalangle)-1<0$), $\|DJ^{r}(u)^{\top}\alpha^{1}\|\neq 0$ and $\|DJ^{r}(u)^{\top}\alpha^{2}\|\neq 0$. If we choose $c_{1}=t$ and $c_{2}=1-t$ for $t\in(0,1)$, then $t\alpha^{1}+(1-t)\alpha^{2}\in\Delta_{k}$ and $\omega(t\alpha^{1}+(1-t)\alpha^{2})<0$, which contradicts $u\in\varphi^{-1}(0)$. If $\|DJ^{r}(u)^{\top}\alpha^{1}\|=0$ or $\|DJ^{r}(u)^{\top}\alpha^{2}\|=0$ then (47) holds for $\bar{\alpha}=\alpha^{1}$ or $\bar{\alpha}=\alpha^{2}$, respectively. If $cos(\sphericalangle)-1=0$ then $DJ^{r}(u)^{\top}\alpha^{1}$ and $DJ^{r}(u)^{\top}\alpha^{2}$ are linearly dependent, so there are $\bar{c}_{1}$, $\bar{c}_{2}\in\mathbb{R}\setminus\{0\}$ such that $DJ^{r}(u)^{\top}\bar{\alpha}=0$ for $\bar{\alpha}=\bar{c}_{1}\alpha^{1}+\bar{c}_{2}\alpha^{2}$. In particular, in any case we must have $\omega(\alpha)=0$ for all $\alpha\in span(\{\alpha^{1},\alpha^{2}\})$. ∎

The previous lemma implies that for $k=2$, the solution of (2.3.2) for $u\in\varphi^{-1}(0)$ is non-unique iff $DJ^{r}(u)DJ^{r}(u)^{\top}-\epsilon\epsilon^{\top}=0$. For $k>2$, we can only have non-uniqueness if (47) holds. If we consider the dimensions of the spaces in (47), we see that in the generic case, it can only hold if

[TABLE]

i.e., if all gradients of the objectives are linearly dependent in $u$. This motivates us to assume that in general, the solution of (2.3.2) is unique for almost all $u\in\varphi^{-1}(0)$.

We will now investigate the differentiability of $\varphi$. Our strategy is to apply the implicit function theorem to the KKT conditions of (2.3.2) to obtain a differentiable function $\phi$ that maps a point $u\in\mathbb{R}^{n}$ onto the solution of (2.3.2) in $u$. This would imply the differentiability of $\varphi$ via concatenation with $\omega$. An obvious problem here is the fact that (2.3.2) has inequality constraints which, when activated or deactivated under variation of $u$, lead to non-differentiabilities in $\phi$. Note that an inequality constraint being active means that one component of $\alpha$ is zero, i.e., one of the objective functions has no impact on the current problem. Thus, for our theoretical purposes, if there is an active inequality constraint in (2.3.2) we will just ignore the corresponding objective function. This approach is strongly related to the hierarchical decomposition of the Pareto critical set (cf. [10]).

For the reasons mentioned above, we will now consider the case where the solution of (2.3.2) is strictly positive in each component. The following lemma shows a technical result that will be used in a later proof.

Lemma A.2**.**

Let $u\in\varphi^{-1}(0)$ and let $\bar{\alpha}\in\Delta_{k}$ be a solution of (2.3.2) with $\alpha_{i}>0$ $\forall i\in\{1,...,k\}$. Then $\bar{\alpha}$ is unique if and only if there is no $\beta\in\mathbb{R}^{k}\setminus\{0\}$ with $\omega(\beta)=0$ and $\sum_{i=1}^{k}\beta_{i}=0$.

Proof.

We will show that $\alpha$ is non-unique if and only if there is some $\beta\in\mathbb{R}^{k}$ with $\omega(\beta)=0$ and $\sum_{i=1}^{k}\beta_{i}=0$.

$\Rightarrow$: Let $\tilde{\alpha}$ be another solution of (2.3.2). Then, as in the proof of Lemma A.1, we must have $\omega(c_{1}\bar{\alpha}+c_{2}\tilde{\alpha})=0$ for all $c_{1},c_{2}\in\mathbb{R}$. This means we can choose $\beta=\bar{\alpha}-\tilde{\alpha}$.

$\Leftarrow$: Let $\beta\in\mathbb{R}^{k}$ with $\omega(\beta)=0$ and $\sum_{i=1}^{k}\beta_{i}=0$. Let $s>0$ be small enough such that $\bar{\alpha}+s\beta\in\Delta_{k}$. Then, as in (A), we have

[TABLE]

Since by assumption $\varphi(u)=0$ we must have $\omega(\bar{\alpha}+s\beta)=0$, so $\bar{\alpha}+s\beta$ is another solution of (2.3.2). ∎

To be able to use the KKT conditions of (2.3.2) to obtain its solution, we have to make sure that these conditions are sufficient. Since (2.3.2) is a quadratic problem, this means we have to show that the matrix in the objective $\omega$ is positive semidefinite.

Lemma A.3**.**

Let $u\in\varphi^{-1}(0)$ and let $\bar{\alpha}\in\Delta_{k}$ be the unique solution of (2.3.2) with $\bar{\alpha}_{i}>0$ $\forall i\in\{1,...,k\}$. Then $\omega(\beta)\geq 0$ for all $\beta\in\mathbb{R}^{k}$. In particular, $DJ(u)DJ(u)^{\top}-\epsilon\epsilon^{\top}$ is positive semidefinite.

Proof.

Assume there is some $\beta\in\mathbb{R}^{k}$ with $\omega(\beta)<0$, i.e., $\epsilon^{\top}\beta>\|DJ^{r}(u)^{\top}\beta\|$. We distinguish between two cases:

Case 1: $\sum_{i=1}^{k}\beta_{i}=0$: Similar to the proof of Lemma A.1 we get

[TABLE]

for all $s>0$. In particular, since $\bar{\alpha}$ is positive, there is some $\bar{s}>0$ such that $\bar{\alpha}+\bar{s}\beta\in\Delta_{k}$ with $\omega(\bar{\alpha}+\bar{s}\beta)<0$, which is a contradiction.

Case 2: $\sum_{i=1}^{k}\beta_{i}\neq 0$. W.l.o.g. assume that $\sum_{i=1}^{k}\beta_{i}=1$. Consider

[TABLE]

Then $\bar{\omega}(0)=0$ and $\bar{\omega}(1)<0$. By assumption we must have $\bar{\omega}(s)>0$ for all $s$ such that $\bar{\alpha}+s(\beta-\bar{\alpha})\in\Delta_{k}$. By continuity of $\bar{\omega}$ there must be some $s^{*}$ with $\bar{\omega}(s^{*})=0$. Let $\bar{\beta}:=\bar{\alpha}+s^{*}(\beta-\bar{\alpha})$. Using (A) we get

[TABLE]

for all $t\in(0,1)$, which is a contradiction. ∎

The previous results now allow us to prove Theorem 2.11.

Theorem 2.11

*Let $\bar{u}\in\varphi^{-1}(0)$ such that (2.3.2) has a unique solution $\bar{\alpha}\in\Delta_{k}$ with $\bar{\alpha}_{i}>0$ for all $i\in\{1,...,k\}$. Let (2.3.2) be uniquely solvable in a neighborhood of $\bar{u}$. Then there is an open set $U\subseteq\mathbb{R}^{n}$ with $\bar{u}\in U$ such that $\varphi|_{U}$ is continuously differentiable. *

Proof.

The KKT conditions for (2.3.2) are

[TABLE]

for $\lambda\in\mathbb{R}$ and $\mu\in\mathbb{R}^{k}$. By Lemma A.3 these conditions are sufficient for optimality. By our assumption there is an open set $U^{\prime}$ with $\bar{u}\in U^{\prime}$ such that the solution of (2.3.2) is unique and positive. Thus, on $U^{\prime}$, (A) is equivalent to

[TABLE]

for some $\lambda\in\mathbb{R}$. This system can be rewritten as $G(u,(\alpha,\lambda))=0$ for

[TABLE]

Derivating $G$ with respect to $(\alpha,\lambda)$ yields

[TABLE]

Let $\bar{\lambda}\in\mathbb{R}$ such that $G(\bar{u},(\bar{\alpha},\bar{\lambda}))=0$. (Note that uniqueness of $\bar{\alpha}$ implies uniqueness of $\bar{\lambda}$ here.) For $D_{(\alpha,\lambda)}G(\bar{u},(\bar{\alpha},\bar{\lambda}))$ to be singular, there would have to be some $v=(v^{1},v^{2})\in\mathbb{R}^{k+1}$ with

[TABLE]

and thus

[TABLE]

By Lemma A.2, this is a contradiction to the assumption that $\bar{\alpha}$ is a unique solution of (2.3.2). So $D_{(\alpha,\lambda)}G(\bar{u},(\bar{\alpha},\bar{\lambda}))$ has to be regular. This means we can apply the implicit function theorem to obtain open sets $U\subseteq U^{\prime}\subseteq\mathbb{R}^{n}$, $V\subseteq\mathbb{R}^{k+1}$ with $\bar{u}\in U$, $(\bar{\alpha},\bar{\lambda})\in V$ and a continuously differentiable function $\phi=(\phi_{\alpha},\phi_{\lambda}):U\rightarrow V$ with

[TABLE]

In particular,

[TABLE]

so $\varphi|_{U}$ is continuously differentiable. ∎

Remark A.4**.**

From the proof of Theorem 2.11 we can even derive an explicit formula for the derivative of $\varphi|_{U}$ in $\bar{u}$: First of all, the derivative of the implicit function $\phi$ is given by

[TABLE]

By applying the chain rule to (50) we obtain

[TABLE]

$\blacksquare$**