A Certified Model Reduction Approach for Robust Parameter Optimization with PDE Constraints

TL;DR

This paper presents a novel certified model reduction method for robust parameter optimization constrained by PDEs, effectively handling uncertainty and reducing computational costs through adaptive surrogate modeling.

Contribution

It introduces an adaptive model order reduction technique with certification for robust PDE-constrained optimization under uncertainty, improving efficiency over traditional methods.

Findings

The proposed method effectively approximates worst-case scenarios in uncertain PDE models.

Numerical results demonstrate significant computational savings and accuracy.

The approach provides certified reduced order models for complex PDE-based optimization.

Abstract

We investigate an optimization problem governed by an elliptic partial differential equation with uncertain parameters. We introduce a robust optimization framework that accounts for uncertain model parameters. The resulting non-linear optimization problem has a bi-level structure due to the min-max formulation. To approximate the worst-case in the optimization problem we propose linear and quadratic approximations. However, this approach still turns out to be very expensive, therefore we propose an adaptive model order reduction technique which avoids long offline stages and provides a certified reduced order surrogate model for the parametrized PDE which is then utilized in the numerical optimization. Numerical results are presented to validate the presented approach.

Figures

Tables

Table 1. Table 1: Comparison of the results obtained by the optimization and robust optimization using different approximation orders.

	$V_{\mathrm{pm}}$	%	$p$	$E_0$ ($90$)	$E_0^{\text{worst}}$
Init	$133.00$	$100$	$(19.00,7.00,7.00)$	$30.3483$	$29.8873$
Nom. Opt	$62.36$	1$47$	$(21.08,2.96,6.62)$	$30.3483$	$29.8873$
Rob. Lin	$66.46$	1$50$	$(21.10,3.15,6.64)$	$30.5817$	$30.2322$
Rob. Quad	$67.95$	1$51$	$(21.10,3.22,6.65)$	$30.6995$	$30.3487$
Rob. Shift	$67.94$	1$51$	$(21.10,3.22,6.65)$	$30.6992$	$30.3483$

Table 2. Table 2: Performance of the SQP method and computational cost.

	CPU (s)	Iter	PDE solves
Nom. Opt	$38.76$	$18$	$72$
Rob. Lin	$42.71$	$8$	$64$
Rob. Quad	$46.47$	$6$	$72$
Rob. Shift	$88.51$	$7$	$84$

Table 3. Table 3: Performance of POD with nominal optimization

iter	$E_0(p,{𝐮}_h\left(p\right),90)$	$E_0(p,{𝐮}_h^{\ell }\left(p\right),90)$	$V$	$\#it$	$\underset{\left|\mu \right|\le 1}{\max }{\Delta }_{u_{\mu }^n}^{\ell }\left(p^k\right)$	$\ell$
0	30.3483	30.3483	133.00
1	30.3764	30.3483	x62.82	8	3.26e-02	4
2	30.3483	30.3483	x62.36	4	2.97e-06	8

Table 4. Table 4: Performance of POD with linear robust optimization

iter	$E_0(p,{𝐮}_h\left(p\right),90)$	$E_0(p,{𝐮}_h^{\ell }\left(p\right),90)$	$V$	$\#it$	$\underset{\left|\mu \right|\le 2}{\max }{\Delta }_{u_{\mu }^n}^{\ell }\left(p^k\right)$	$\ell$
0	30.3483	30.3483	133.00
1	30.7319	30.5729	x69.41	8	7.12e-01	8
2	30.5817	30.5816	x66.46	4	9.83e-05	16

Table 5. Table 5: Performance of POD with quadratic robust optimization

iter	$E_0(p,{𝐮}_h\left(p\right),90)$	$E_0(p,{𝐮}_h^{\ell }\left(p\right),90)$	$V$	$\#it$	$\underset{\left|\mu \right|\le 3}{\max }{\Delta }_{u_{\mu }^n}^{\ell }\left(p^k\right)$	$\ell$
0	30.3483	30.3483	133.00
1	30.8112	30.7332	x68.12	6	3.32e-02	8
2	30.6997	30.6995	x67.96	5	8.07e-05	16

Table 6. Table 6: Performance of POD with quadratic robust optimization with shift

iter	$E_0(p,{𝐮}_h\left(p\right),90)$	$E_0(p,{𝐮}_h^{\ell }\left(p\right),90)$	$V$	$\#it$	$\underset{\left|\mu \right|\le 3}{\max }{\Delta }_{u_{\mu }^n}^{\ell }\left(p^k\right)$	$\ell$
0	30.3483	30.3483	133.00
1	30.7882	30.7384	x67.59	9	4.05e-02	8
2	30.7908	30.7913	x67.74	9	7.11e-04	15
3	30.7934	30.7933	x67.87	9	7.29e-05	16

Table 7. Table 7: Comparison of the results obtained by POD coupled with the optimization and robust optimization using different approximation orders.

	$V_{\mathrm{pm}}$	%	$p$	$E_0$ ($90$)	$E_0^{\text{worst}}$
Init	133.00	100	(19.00 7.00 7.00)	30.34	29.88
Nom. Opt	61.24	46%	(21.25 2.88 6.87)	30.23	29.77
Rob. Lin	66.71	50%	(20.76 3.21 6.14)	30.56	30.21
Rob. Quad	67.48	51%	(21.21 3.17 6.92)	30.76	30.29
Rob. Shift	67.86	51%	(21.31 3.18 6.96)	30.79	30.32

Table 8. Table 8: Performance of the POD method and computational cost.

	CPU (s)	PDE solves	Reduced PDE solves
Nom. Opt	$7.7$	$8$	$48$
Rob. Lin	$11.3$	$24$	$110$
Rob. Quad	$14.2$	$36$	$132$
Rob. Shift	$20.9$	$48$	324

Equations

$$\left\{\begin{array}[]{l}\hbox{For }p\in\mathcal{D}\hbox{ find }\tilde{u}\in X(\Omega(p))\hbox{ s.t}\\ \tilde{a}(\tilde{u},v;p)=\tilde{f}(v;p,\phi),\quad\forall v\in X(\Omega(p)),\end{array}\right.$$

$$\tilde{a}(w,v;p)=\int_{\Omega(p)}\nu\nabla w\cdot\nabla v\,\mathrm{d}x\quad\mbox{and}\quad\tilde{f}(v;p,\phi)=\int_{\Omega(p)}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tilde{f}(\phi)}v\,\mathrm{d}x,$$

$$H^1(\Omega (p)):=\left\{f\in L^2(\Omega (p)):D^{\alpha }f\in (L^2(\Omega (p)){)}^2,|\alpha |\le 1\right\},$$

$$L^2(\Omega (p)):=\left\{f:\Omega (p)\to \mathbb{R},{\int }_{\Omega (p)}f(x{)}^2\mathrm{d}x<\infty \right\},$$

$$H_0^1(\Omega (p)):=\left\{f\in H^1(\Omega (p)):f\equiv 0\; a.e.\; on\partial \Omega (p)in\; the\; sense\; of\; traces\right\}.$$

$${\mathcal{T}}_q(x;p):{\Omega }_q(\bar{p})$$

$$x$$

$$\left\{\begin{array}[]{l}\hbox{For }p\in\mathcal{D}\hbox{ find }u\in X(\bar{\Omega})\hbox{ s.t.}\\ a(u,v;p)=f(v;p,\phi),\quad\forall v\in X(\bar{\Omega}),\end{array}\right.$$

$$a(u,v;p)=\sum _{q=1}^{Q_a}\sum _{i,j=1}^2{\hat{\Phi }}_q^{i,j}(p)a_q^{i,j}(u,v),$$

$$a_{q}^{i,j}(w,v):=\int_{\Omega_{q}(\bar{p})}\dfrac{\partial w}{\partial x_{i}}\dfrac{\partial v}{\partial x_{j}}\,\mathrm{d}x,\quad\mbox{for}\quad q=1,\ldots,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}Q_{a}}\mbox{ and }i,j=1,2.$$

$$f(v;p,\varphi )=\sum _{q=1}^{Q_f}{h}_q(\varphi ){\hat{\Phi }}_q^f(p)f_q(v)$$

$$f_{q}(v):=\int_{\Omega_{q}(\bar{p})}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(}f{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\circ\mathcal{T}_{q}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0})}v{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}|\det J_{\mathcal{T}_{q}}|}\,\mathrm{d}x,\quad\mbox{for}\quad q=1,\ldots,Q_{f}.$$

$$\alpha (p):=\underset{w\in X(\bar{\Omega }),w\ne 0}{\inf }\frac{a(w,w;p)}{\Vert w{\Vert }_{X(\bar{\Omega })}^2},$$

$$\gamma (p):=\underset{v\in X(\bar{\Omega }),v\ne 0}{\sup }\underset{w\in X(\bar{\Omega }),w\ne 0}{\sup }\frac{a(w,v;p)}{\Vert w{\Vert }_{X(\bar{\Omega })}\Vert v{\Vert }_{X(\bar{\Omega })}}.$$

$$\left\{\begin{array}[]{l}\hbox{For }p\in\mathcal{D}\hbox{ find }u^{1}_{i}\in X(\bar{\Omega)}\hbox{ s.t.}\\ a(u^{1}_{i},v;p)=\frac{\partial f}{\partial p_{i}}(v;p)-\frac{\partial a}{\partial p_{i}}(u,v;p),\quad\forall v\in X(\bar{\Omega}),\,i=1,\ldots,N_{p}.\end{array}\right.$$

$$a(u_{\mu}^{n},v;p)=\frac{\partial^{\mu}f}{\partial p^{\mu}}(v;p)-\sum_{1\leq|\kappa|\leq n}\binom{\mu}{\kappa}\frac{\partial^{\kappa}a}{\partial p^{\kappa}}(u^{n-|\kappa|}_{\mu-\kappa},v;p)\quad\forall v\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X(\bar{\Omega})},$$

$$\left\{\begin{array}[]{l}\hbox{For }p\in\mathcal{D}\hbox{ find }u_{h}\in X_{h}(\bar{\Omega})\hbox{ s.t.}\\ a(u_{h},v_{h};p)=f(v_{h};p),\quad\forall v_{h}\in X_{h}(\bar{\Omega}).\end{array}\right.$$

$$\left\{\begin{array}[]{l}\hbox{For }p\in\mathcal{D}\hbox{ find }{\bf u}_{h}(p)\in\mathbb{R}^{N}\hbox{ s.t.}\\ {\bf K}(p){\bf u_{h}}(p)={\bf f}(p),\end{array}\right.$$

$$\mathbf{K}(p)=\sum _{q=1}^Q\sum _{i,j=1}^2{\hat{\Phi }}_q^{i,j}(p){\mathbf{K}}_q^{i,j},$$

$$\mathbf{K}(p){\mathbf{u}}_{h,i}^1(p)={\tilde{\mathbf{f}}}_i^1,\text{mbox}fori=1,\dots ,N_p,$$

$${\tilde{\mathbf{f}}}_i^1=\frac{\partial \mathbf{f}}{\partial p_i}(p)-\frac{\partial \mathbf{K}}{\partial p_i}(p){\mathbf{u}}_h\text{mbox}fori=1,\dots ,N_p,$$

$$\mathbf{K}(p){\mathbf{u}}_{h,\mu }^n=\frac{{\partial }^{\mu }\mathbf{f}}{\partial p^{\mu }}(p)-\sum _{1\le |\kappa |\le n}\left(\genfrac{}{}{0ex}{}{\mu }{\kappa }\right)\frac{{\partial }^{\kappa }\mathbf{K}}{\partial p^{\kappa }}(p){\mathbf{u}}_{h,\mu -\kappa }^{n-|\kappa |}.$$

$${\mathbf{u}}_{h,\mu }^n(p):=\frac{{\partial }^{\mu }{\mathbf{u}}_h}{\partial p^{\mu }}(p)$$

$$\begin{array}[]{rl}\displaystyle\min_{p\in\mathbb{R}^{N_{p}},{\bf u}_{h}\in\mathbb{R}^{N}}&\tilde{g}_{0}({\bf u}_{h},p,\hat{\phi}),\\ \mbox{subject to (s.t.)}&\tilde{g}_{i}({\bf u}_{h},p,\hat{\phi})\leq 0,\quad i=1,\ldots,N_{g},\\ &e({\bf u}_{h},p,\hat{\phi})=0,\end{array}$$

$$g_i(p,\hat{\varphi }):={\tilde{g}}_i({\mathbf{u}}_h,p,\hat{\varphi })\text{mbox}fori=0,\dots ,N_g.$$

$$\begin{array}[]{rl}\displaystyle\min_{p\in\mathbb{R}^{N_{p}}}&g_{0}(p,\hat{\phi}),\\ \mbox{s.t.}&g_{i}(p,\hat{\phi})\leq 0,\quad i=1,\ldots,N_{g}.\end{array}$$

$$\mathcal{U}_{k}=\left\{\phi\in\mathbb{R}^{N_{\phi}}\,\big{|}\,\|D^{-1}(\phi-\hat{\phi})\|_{k}\leq 1\right\},$$

$$\mathcal{U}_{\infty}=\{\phi\in\mathbb{R}^{N_{\phi}}\,\big{|}\,\underline{\phi}\leq\phi\leq\overline{\phi}\}.$$

$${\phi }_i(p):=\underset{\varphi \in {\mathcal{U}}_k}{\max }{g}_i(p,\varphi ),i=0,\dots ,N_g.$$

$$\begin{array}[]{rl}\varphi_{i}(p):=\displaystyle\max_{\phi\in\mathbb{R}^{N_{\phi}}}&g_{i}(p,\phi)\\ \mbox{s.t.}&\|D^{-1}(\phi-\hat{\phi})\|_{k}\leq 1,\quad i=0,\ldots,N_{g}.\end{array}$$

$$\begin{array}[]{rl}\displaystyle\min_{p\in\mathbb{R}^{N_{p}}}&\varphi_{0}(p)\\ \mbox{s.t.}&\varphi_{i}(p)\leq 0,\quad i=1,\ldots,N_{g}.\end{array}$$

Full text

††journal: Advances in Computational Mathematics

1 Introduction

Parameter optimization governed by partial differential equations (PDEs) is a well-studied topic due to its relevance in industrial applications. If the model is considered to be perfect one can perform standard techniques to solve these type of problems (see e.g. [20, 39] and the references therein). However, its numerical approximation can be very expensive due to the dimension of the discretized PDE. For this reason, in the last decade, model order reduction techniques were introduced and successfully applied in the context of PDE constrained optimization. Model order reduction works in a Galerkin projection framework, where the basis functions are non-local and built upon information of the underlying system. Although a detailed description of these contributions goes beyond the scope of this paper, we want to mention Proper Orthogonal Decomposition (POD) and the reduced basis (RB) method. The former works in general situations such as time-dependent problems, and parametric steady and unsteady equations (see e.g. [42] for a presentation of the method within different applications). The latter is mainly used in the context of e.g. many-query scenarios (see e.g. [37]) for parametric steady problems where the basis functions are selected by means of a greedy algorithm. It is also possible to combine RB and POD in the so-called POD-Greedy algorithm [18] for parametric unsteady models. The strength of these methods is the presence of an a-posteriori error estimator which certifies the quality of the surrogate models. Model order reduction has been applied successfully to parameter optimization with PDE constrains in e.g. [13, 31]. However, we note that all these works strongly rely on the offline-online decomposition. In the current work, we propose to build the surrogate model towards the optimization. Similar approaches were studied in e.g. [35, 43], but our approach does not require offline-online decomposition and the reduced order model gets updated to desired accuracy during the optimization process.

For the purpose of this work we assume that the PDE which governs our system is given but material imperfections are present. These are due to e.g. manufacturing and introduce uncertainty to the model. This problem arises in many real life applications. One way to include this uncertainty into the optimization process is through robust optimization. In this case no probabilistic model of the uncertainty is required. Instead, a deterministic approach is applied by assuming that the uncertainty is restricted to a bounded uncertainty set. Using the notion of a robust counterpart the associated original uncertain optimization problem is reformulated. The solution obtained in this way stays feasible for all realizations from the uncertainty set and, at the same time, yields the best worst-case value of the objective function. Techniques which resolve the uncertainty by means of stochastic optimization e.g. [9, 25, 30, 40] are an alternative to the presented approach. These methods depend on sampling the uncertain parameters and hence can become prohibitive expensive in the context of PDE constrained optimization. For a general discussion of robust optimization we refer to e.g. [6, 7, 8].

The idea is to utilize a suitable approximation of the robust counterpart, e.g. [11, 22, 26, 38, 44] or exploit specific properties e.g. [22]. We investigate the approximation of the robust counterpart using a quadratic model. This allows us to reformulate the robust optimization problem as a mathematical program with equilibrium constraints (MPEC). This approach has been investigated in [26, 38, 3] in the context of PDE constrained optimization problems, where [26] forms the starting point for the present investigations. Our model will be a linear elliptic parametric equation where an affine decomposition will be applied to work with a reference parameter. The worst-case problem leads a non-linear optimization problem with a $\min-\max$ formulation. After approximating the inner maximization derivative based optimization is applied. In the setting of our work we will utilize a sensitivity based approach since we assume to deal with only a few parameters. This approach is computationally expensive because it requires the solution of multiple PDEs.

In this work the focus is on solving the robust counterpart in an efficient way utilizing a POD-based reduced order method. More specifically, we propose a certified robust optimization procedure which, following an idea proposed in [1], successively adapts the POD model during the iterative adjustment of the optimal robust design through appropriate enrichment of the snapshot set and avoids expensive offline stages. The method is certified by an a-posteriori error estimator for the state variable and the sensitivities. Therefore, error estimators for the state and the sensitivities are required. A generalized error estimator is derived to cover these needs. Here we build upon the robust optimization framework developed in [26], which for the convenience of the reader is summarized in Section 3, and which we have now complemented with a moving expansion point approach to further improve the accuracy of the utilized approximations, see Section 3.2.3 for details. Our approach to model reduction complements the method proposed in [26]. In fact, in [26] the model reduction in the application is performed by a POD-greedy procedure to reduce the degrees of freedom related to the rotation angle of a quasistatic model of an electrical machine. This POD surrogate model is then used to speed up the robust optimization. In this paper, we generate the reduced order model by applying POD to a snapshot set resulting from state and sensitivities corresponding to previous optimization iterates. The quality of the obtained optimal solution is evaluated by a posteriori error estimators. As long as the error is above a specified level, the snapshot set is updated by the state and sensitivities corresponding to the current iterate of the optimization variables. We demonstrate the performance and efficiency of our approach for robust design of electrical machines, where we use the setting of [26] in the static case, i.e., we do not consider rotation.

The paper is organized as follows: in Section 2 we present the mathematical model and in Section 3 we present the nominal and the robust optimization problem. Section 4 focuses on model order reduction for the optimization problem and in Section 5 we illustrate the effectiveness of the discussed methods by numerical examples.

2 Elliptic Parametric PDEs

We deal with an abstract linear stationary equation governed by a uniformly continuous and uniformly coercive bilinear form $\tilde{a}$ stemming e.g. from a parametric elliptic PDE. The parameter sets are denoted by $\mathcal{D}\subset\mathbb{R}^{N_{p}}$ and $\mathcal{U}_{k}\subset\mathbb{R}^{N_{\Phi}}$, where parameters $p\in\mathcal{D}$ enter the problem via a $p$-dependent regular bounded domain $\Omega(p)\subset\mathbb{R}^{2}$ which give rise to a $p$-dependent real and separable Hilbert space $X(\Omega(p))$ and dual space $X^{\prime}(\Omega(p))$. The parameters $\Phi\in\mathcal{U}_{k}$ will model the uncertainty in the problem. We note that for our abstract setting we impose the by now standard assumptions on parameter separability which are formulated for example in [17]. Those assumptions are in particular met for the application considered in Section 5. Our abstract problem then reads

[TABLE]

where $\tilde{u}$ is the unknown variable. In our application uncertainty only enters through the right hand side, where $\tilde{a}$ and $\tilde{f}$ are given by

[TABLE]

where $\nu\geq\bar{\nu}>0$ is an isotropic material coefficient (see Section 3 for more details.). Note that under these assumptions problem (1) admits a unique solution.

The function space $X(\Omega(p))$ is such that H^{1}_{0}(\Omega(p))\subset{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X(\Omega(p))}\subset H^{1}(\Omega(p)), with

[TABLE]

where $f$ is a measurable function, $\alpha\in\mathbb{N}_{0}^{2}$, $D^{\alpha}f$ denotes the weak $\alpha^{th}-$ partial derivative of $f$ and the functional spaces $L^{2}(\Omega(p))$, $H^{1}_{0}(\Omega(p))$ are defined as follows:

[TABLE]

[TABLE]

In order to obtain a computationally fast model and to avoid re-meshing when the parameter changes we require parameter separability, see, e.g., [17], and assume that the domains $\Omega(p)$ in (1) can be transformed to a fixed domain $\Omega(\bar{p}):=\cup_{q=1}^{Q}\Omega_{q}(\bar{p})$, where $\bar{p}$ denotes the reference parameter. Without loss of generality, we assume that the domain of interest $\Omega(p)$ can be decomposed in $Q$ non-overlapping triangles and the transformation on each triangle is linear, whereas piecewise-linear and continuous over the whole domain according to:

[TABLE]

for $q=1,\ldots,Q,$ where $C_{q}(p)\in\mathbb{R}^{2\times 2}$ and $d_{q}(p)\in\mathbb{R}^{2}$. From here onwards we write $\bar{\Omega}$ instead $\Omega(\bar{p})$ whenever this is appropriate.

As it is shown, e.g., in [37], the linear parametric map $\mathcal{T}_{q}(p)$ and its Jacobian $J_{\mathcal{T}_{q}}$ allow the definition of the bilinear and linear forms on the reference domain $\bar{\Omega}$. Then, problem (1) can be rewritten as

[TABLE]

where the assumptions in problem (1) hold true. The bilinear form $a(\cdot,\cdot;p)$ can be expressed with an affine linear decomposition:

[TABLE]

such that $\hat{\Phi}_{q}^{i,j}:\mathcal{D}\rightarrow\mathbb{R}\hbox{ for }q=1,\ldots,Q_{a},i,j=1,2,$ is a function depending on $p$ and $a_{q}^{i,j}$ are the parameter independent continuous bilinear forms. In our example we have

[TABLE]

This allows us to improve the computational efficiency in the evaluation of $a(u,v;p)$; the discrete approximation of the components $a_{q}^{i,j}(u,v)$ can be then computed once. The same affine decomposition is applied to $f$ with

[TABLE]

with $h_{q}(\phi)$ being nonlinear differentiable functions depending on the uncertain parameter and

[TABLE]

The assumption that the problem is affine dependent on the parameter $p$ is the key for the purpose of model reduction as we will see in Section 4. We note that without loss of generality, in the present work, we consider the same number of subdomains $Q=Q_{a}=Q_{f}$ for the (bi)linear forms $a$ and $f$.

Furthermore, we define the coercivity constant that will be a key ingredient in the certification of the model reduction algorithm, by

[TABLE]

and the continuity constant as

[TABLE]

To ease the notations, from here onwards, we drop the dependence on $\phi$ (e.g. $h_{q}(\phi):=1$ in (5)).

For the purpose of the optimization problem we will compute the sensitivities ${u}^{1}_{i}:=\frac{\partial u(p)}{\partial p_{i}}$ which are obtained by the derivative with respect to the parameters from equation (3) leading the following linear problem:

[TABLE]

We note that due to the affine decomposition the $i$-th partial derivatives of $a(\cdot,\cdot;p)$ and $f(\cdot;p)$ are given by the derivatives of coefficient functions $\hat{\Phi}_{q}^{i,j}$ and $\hat{\Phi}_{q}^{f}$, $q=1,\ldots,Q$, $i,j=1,2$ and can be computed analytically. More generally, for a multiindex $\mu\in\mathbb{N}_{0}^{N_{p}}$, $|\mu|=n$, the $\mu$-th sensitivity can be computed by the following proposition.

Proposition 2.1**.**

Let the coefficient functions $\hat{\Phi}_{q}^{i,j}(p)$ and $\hat{\Phi}_{q}^{f}(p,\phi)$, $q=1\ldots,Q$, $i,j=1,2$ be $n$-times differentiable with respect to $p$. Then the solution $u$ is differentiable with respect to $p$ and the sensitivities u_{\mu}^{n}=\frac{\partial^{\mu}u}{\partial p^{\mu}}\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X(\bar{\Omega})}, $|\mu|=n$, satisfy the sensitivity equation

[TABLE]

where $\binom{\nu}{\kappa}=\prod_{i=1}^{N_{p}}\binom{\mu_{i}}{\kappa_{i}}=\frac{\mu!}{\kappa!(\mu-\kappa)!}$ denotes the binomial coefficient for multiindices.

The proof of the proposition follows from the general Leibniz rule for the $\mu$-th derivative applied to (4) and is omitted here. We can easily see that for $n=0$ we find the state equation (3) and for $n=1$ the first sensitivity equation (8).

Next, we introduce an high dimensional finite element (FE) discretization of our model problem in the space X_{h}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\bar{\Omega})}\subset X{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\bar{\Omega})}. The discrete problem then reads

[TABLE]

For the discrete problem we use the ansatz $u_{h}=\sum_{i=1}^{N}({\bf u}_{h})_{i}\varphi_{i}$, where $\varphi_{i}$ are suitable FE ansatz functions. Problem (10) then is equivalent to the linear system

[TABLE]

where ${\bf K}(p)\in\mathbb{R}^{N\times N}$ is the stiffness matrix $\left({\bf K}(p)\right)_{ij}=a(\varphi_{j},\varphi_{i};p)$, $1\leq i,j\leq N$, and the right hand side ${\bf f}(p)\in\mathbb{R}^{N}$ is obtained by $\left({\bf f}(p)\right)_{i}=f(\varphi_{i};p)$, $1\leq i\leq N$. The FE system matrix keeps the dependency on the parameter $p$ and we apply the affine decomposition to it in the following way

[TABLE]

where ${\bf K}_{q}^{i,j}$, $q=1,\ldots,Q$ and $i,j=1,2$, are the system matrices on the $Q$ sub-domains. We note again that the $p$ dependency of ${\bf K}$ appears only in the weight functions $\hat{\Phi}_{q}^{i,j}$ which are easy and inexpensive to evaluate. The same applies for the right hand side. Similarly, the discretized version of equation (8) reads

[TABLE]

with

[TABLE]

where the subindex $i$ indicates as above the derivative with respect to the $i$-th parameter $p_{i}$. Note that these derivatives are easy to compute due to the previously introduced affine decomposition (12). The derivative of the matrix ${\bf K}(p)$ and the vector ${\bf f}(p)$ are given by the derivatives of the functions $\hat{\Phi}_{q}^{i,j}(p)$ and $\hat{\Phi}_{q}^{f}(p,\phi)$, $q=1,\ldots,Q$, $i,j=1,2$, respectively. For the general sensitivity equations (9) we get with a multiindex $\mu\in\mathbb{N}_{0}^{N_{p}}$, $|\mu|=n$,

[TABLE]

In the rest of the work we will continue to use the following compact notation when appropriate:

[TABLE]

in order to denote the $\mu$-th derivative of ${\bf u}_{h}$ with respect to $p$ for a multindex $\mu\in\mathbb{N}_{0}^{N_{p}}$, $|\mu|=n$. Moreover, we sometimes use the fact that the state is the zeroth order sensitivity, i.e., ${\bf u}_{h,\mu}^{0}(p)={\bf u}_{h}(p)$ with $|\mu|=0$. Similarly, we will adopt ${\bf f}^{n}_{\mu}$ and ${\bf K}^{n}_{\mu}$.

From now on we will focus on the discrete version (11) of the problem (1). All further steps are analogous in the continuous setting.

3 Optimization problem

This section is devoted to the optimization and the robust optimization problem. We will formulate an optimization problem governed by a parametrized PDE, then investigate the nominal optimization and its robust counterpart. By the robust formulation, uncertainties in model parameters are taken into account. This is done by utilizing a worst-case formulation. For an efficient realization, approximation techniques of different orders are investigated. Here we closely follow [26] and repeat the ingredients for the convenience of the reader.

3.1 The nominal optimization

We first introduce the nominal optimization problem. In this setting all model parameters $\phi$ are fixed to a nominal value $\hat{\phi}$. The optimization problem then reads as:

[TABLE]

where $p\in\mathbb{R}^{N_{p}}$ is the optimization variable and $\hat{\phi}\in\mathcal{U}_{k}\subset\mathbb{R}^{N_{\phi}}$ is the fixed model parameter. The function $\tilde{g}_{0}$ is the objective function and $\tilde{g}_{i}$ are the $N_{g}\in\mathbb{R}$ inequality constraints with $\tilde{g}_{i}:\mathbb{R}^{N}\times\mathbb{R}^{N_{p}}\times\mathbb{R}^{N_{\phi}}\mapsto\mathbb{R}$ for $i=0,\ldots,N_{g}$. Furthermore, $e$ is an equality constraint governed by the discretized parametric PDE in (11), i.e., e({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bf u}_{h},p,\hat{\Phi})=0\Longleftrightarrow K(p,\hat{\Phi})u_{h}=f(p,\hat{\Phi}). We assume $\tilde{g}_{i}$, $i=0,\ldots,N_{g}$, and $e$ to be continuously differentiable. It directly follows from the coercivity of the bilinear form $a$ that the parametric PDE for every admissible parameter $p$ admits a unique solution $u_{h}$, and that the Jacobian \frac{\partial e}{\partial{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bf u}} is boundedly invertible, compare for example [20]. Further, let $e({\bf u}_{h},p,\phi)=0$ have a unique solution ${\bf u}_{h}={\bf u}_{h}(p)$ for every admissible $p$. Then we can introduce the reduced objective and constraint functions

[TABLE]

The reduced optimization problem associated with (15) then reads

[TABLE]

Note that the reduced objective and the reduced constraints $g_{i}(p)$, $i=0,\ldots,N_{g}$, are continuously differentiable due to the previous assumptions. Since $e$ is continuously differentiable with invertible Jacobian {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\partial}{\partial{\bf u}}e({\bf u},p,\phi)}, the implicit function theorem guarantees that also ${\bf u}_{h}$ is continuously differentiable with respect to the parameter. In the further investigation, we focus on the reduced formulation.

3.2 The robust optimization

After having introduced the nominal optimization problem we will now introduce the robust optimization framework, where the model parameter $\phi$ is considered to be uncertain. In the presented setting, the uncertain parameter is a model parameter but the framework is not limited to this setting. We will next formulate a robust version of (16). For this, we first have to make some assumptions on the uncertain parameter and introduce the uncertainty set

[TABLE]

where $\hat{\phi}$ is a nominal value and $D\in\mathbb{R}^{N_{\phi}\times N_{\phi}}$ an invertible scaling matrix. We assume that the uncertain parameter remains within the given bounded set $\mathcal{U}_{k}$. For different choices of $k$ we get different uncertainty sets. The two most commonly used ones are $k=\{2,\infty\}$. The case $k=\infty$ is a special case and allows a representation using upper and lower bounds $\overline{\phi},\underline{\phi}$, respectively. For this we set $D=\mathrm{diag}((\overline{\phi}-\underline{\phi})/2)$ and $\hat{\phi}=(\underline{\phi}+\overline{\phi})/2$. The uncertainty set can then be written as

[TABLE]

Next we want to focus on the robust optimization problem utilizing the worst-case formulation [7, 11, 44]. Using the uncertainty set we can define the worst-case function as

[TABLE]

For an interpretation, for every fixed parameter $p$ the function $\varphi_{i}:\mathcal{U}_{k}\rightarrow\mathbb{R}$ is the worst-case value of the function $g_{i}$ with $\phi\in\mathcal{U}_{k}$. Using the definition of the uncertainty set, the worst-case function can be rewritten as

[TABLE]

Using the worst case function we will formulate the robust optimization problem. The robust counterpart of (16) is given by

[TABLE]

A solution $p$ of (19) is robust against uncertainties in the parameter $\phi$ and is hence referred to as robust optimal solution. Note that the solution is feasible for (16) for all $\phi\in\mathcal{U}_{k}$ and optimal with respect to the objective function.

The introduced robust optimization problem (19) is of bi-level structure and hence difficult to solve. Thus it is required to develop tailored methods to solve the problem. For general nonlinear problems in [7], it is proposed to employ approximations. The approximated robust counterpart of (19) is then introduced as

[TABLE]

where $\hat{\varphi}_{i}$ can be computed more efficiently compared to (18) and is referred to as the approximated worst-case function. First and second order approximations $\hat{\varphi}_{i}$ of the worst-case function $\varphi_{i}$ will be considered. The first order approximation for the general nonlinear case was investigated in [11, 44] while the second order approximation is a recent approach [26] and can be seen as a modification of [38].

3.2.1 Linear approximation of the robust counterpart

In the first order approach suggested in [11, 44], a linearization of the worst-case function is carried out. For this a nominal value $\hat{\phi}$ for the uncertain parameter is chosen and the approximated worst-case function is then given by the first order Taylor expansion

[TABLE]

where $\delta_{i}=\phi-\hat{\phi}$. The solution of the resulting convex optimization problem can be given analytically in the form

[TABLE]

with $\|D\cdot\|_{k^{*}}$ the dual norm of $\|D^{-1}\cdot\|_{k}$ for $k^{*}=k/(k-1)$ with $k^{*}=1$ for $k=\infty$. Note that this is a standard result for scaled Hölder norms. Utilizing these results, the linear approximated robust counterpart reads as

[TABLE]

We refer to (21) as the linear approximated robust counterpart of (16). Due to the norm the objective function and the inequality constraints are non-differentiable if the term inside the norm becomes zero. For the case $k=\infty$, the remedy is to introduce slack variables, i.e.,

[TABLE]

For the case $k=2$, the problem (21) can be reformulated by using second order cone constraints that can be handled effciently for example by interior point techniques. But usually the solution of (21) is sufficiently far away from points of nonsmoothness and in practice standard algorithms for smooth nonlinear problems work well without any modification. Nevertheless, one has to verify that the robust optimal solution does not lie in a non-differentiable point [11, 38].

The required derivatives for the approximations in this approach and the optimization can be computed using either the adjoint or the sensitivity approach [20]. Note that the adjoint approach is desirable when the number of uncertain parameters is large. For a detailed discussion about the different approaches we refer the reader to [11].

While computationally attractive, the first order approximation can suffer from an inaccurate approximation. Hence the influence of the uncertain parameter $\phi$ might be described insufficiently, which was already observed in [12]. This leads to the development of second order methods [26, 38].

3.2.2 Quadratic approximation of the robust counterpart

In the second order approach, a quadratic approximation of the worst case function is utilized. This idea has been investigated in [26, 38]. We will apply the method directly to the reduced problem (16). Since this approach is recent, we will give a short overview of the strategy.

The second order approximation is generated using the second order Taylor expansion of $g_{i}(p,\cdot)$ around $\hat{\phi}$, i.e.,

[TABLE]

where $\alpha_{i}(p,\hat{\phi})=g_{i}(p,\hat{\phi})$, $c_{i}(p,\hat{\phi})=\nabla_{\phi}g_{i}(p,\hat{\phi})$, and $H_{i}(p,\hat{\phi})=\nabla_{\phi\phi}g_{i}(p,\hat{\phi})$. As before, we have set $\delta_{i}=\phi-\hat{\phi}$. The corresponding approximated worst-case function is

[TABLE]

For the quadratic approximation we do not have a closed solution as in the linear case. Fortunately, the problem exhibits a well known structure. In the case $k=2$, the problem corresponds to a trust region problem and is well-studied. The exact solutions of (23) are characterized as follows.

Theorem 3.1**.**

For given $p$, the vector $\delta_{i}\in\mathbb{R}^{N_{\phi}}$ is a global solution of the trust region problem (23) if and only if there exists a Lagrange multiplier $\lambda_{i}\in\mathbb{R}$ satisfying

[TABLE]

with -H({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}p},\hat{\phi})+\lambda_{i}\mathbb{D} positive semidefinite and $\mathbb{D}=D^{-\top}D^{-1}$.

The proof of this theorem can be found in [10, 20]. By adding a square to the norms in the constraints, we obtain the equivalent but differentiable formulation of the quadratic approximated robust counterpart of (16) by

[TABLE]

where $A\preceq 0$ denotes that $A$ is a negative semidefinite matrix. The optimization variables are $p\in\mathbb{R}^{N_{p}}$, $\delta_{0},\ldots,\delta_{N_{g}}\in\mathbb{R}^{N_{\phi}}$, and $\lambda_{0},\ldots,\lambda_{N_{g}}\in\mathbb{R}$. The semidefinite constraint can be reformulated using the smallest eigenvalues as was outlined in [38] or be treated, e.g., by interior point techniques. Problem (24) is a mathematical program with equilibrium constraints (MPEC). To solve these problems numerically one has to pay attention to the complementarity constraint. It turns out that SQP method under relatively mild assumptions and few modifications [15, 27] can solve these type of problems efficiently. In [26] the required changes are outlined in details for the given problem.

The computation of the derivatives required for the quadratic approximation can be done by different approaches. Depending on the number of parameters, the sensitivity or the adjoint approach should be chosen. In this work, it is assumed that the number of uncertain parameters is small and the sensitivity based approach is utilized.

3.2.3 Improving the approximation by moving the expansion point

The quadratic approximation is often a notable improvement over the linear approximation. Nonetheless, Taylor approximations generally are accurate only locally, that is, in a neighborhood of the expansion point $\hat{\phi}$, and so the second-order model can sometimes be a poor predictor for the effects of the uncertain parameters, especially when the uncertainty set is relatively large and the expansion point is far away from the worst case. In these cases, the accuracy can be improved by moving the expansion point of the Taylor approximation towards a presumed maximum of the modeled function. This can be implemented with the aid of an additional set of variables $\hat{\phi}_{i}\in\mathbb{R}^{N_{\phi}}$, $i=0,\ldots,N_{g}$, for the individual expansion points of the models. The second-order Taylor expansion of g_{i}(p,\,\mathchoice{{\hbox to4.99995pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{1.5}{\raisebox{0.0pt}{\displaystyle\cdot}}}}\hss}}}{{\hbox to4.99995pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{1.5}{\raisebox{0.0pt}{\textstyle\cdot}}}}\hss}}}{{\hbox to2.91666pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{1.5}{\raisebox{0.0pt}{\scriptstyle\cdot}}}}\hss}}}{{\hbox to1.80553pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{1.5}{\raisebox{0.0pt}{\scriptscriptstyle\cdot}}}}\hss}}}\,) around $\hat{\phi}_{i}$ evaluated near the center $\hat{\phi}$ of the uncertainty set is given by

[TABLE]

where $\hat{\delta}_{i}=\hat{\phi}_{i}-\hat{\phi}$ as well as $H_{i}(p,\hat{\phi}_{i})=\nabla_{\phi\phi}g_{i}(p,\hat{\phi}_{i})$, and additionally

[TABLE]

Since a good choice for $\hat{\phi}_{i}$ is often not known a priori, it can be adaptively chosen in the course of the optimization. A schematic algorithm that does so is stated in Algorithm 1, which is based on a basic version of the SQP method by Powell, see [34].

The algorithm generates a sequence of optimization variables $x^{k}\coloneqq(p^{k},(\delta_{i}^{k})_{i},(\lambda_{i}^{k})_{i})$ and a sequence of expansion points $\hat{\Phi}^{k}\coloneqq((\hat{\phi}_{i}^{k})_{i})$, where $k\in\mathbb{N}$ is the iteration counter. The feasibility is eventually enforced with an exact $\ell_{1}$-penalty function, denoted here by $P_{\rho}^{1}(x;\hat{\Phi})$, with a penalty parameter $\rho>0$. Recall that, if we write (24) in the form $\min_{x}\,\{\,f(x;\hat{\Phi}):h_{\leq}(x;\hat{\Phi})\leq 0,\,h_{=}(x;\hat{\Phi})=0\,\}$ for brevity, the penalty function is defined by

[TABLE]

where the maximum is to be understood componentwise, and \lVert\,\mathchoice{{\hbox to4.99995pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{1.5}{\raisebox{0.0pt}{\displaystyle\cdot}}}}\hss}}}{{\hbox to4.99995pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{1.5}{\raisebox{0.0pt}{\textstyle\cdot}}}}\hss}}}{{\hbox to2.91666pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{1.5}{\raisebox{0.0pt}{\scriptstyle\cdot}}}}\hss}}}{{\hbox to1.80553pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{1.5}{\raisebox{0.0pt}{\scriptscriptstyle\cdot}}}}\hss}}}\,\rVert_{1} denotes the $1$-norm. We assume that a sufficiently large $\rho$ is known in advance, but in practise it is usually chosen dynamically.

The quadratic program in line 2 is obtained from (24) at the current iterate in the standard way, i.e., by linearization of the constraints and by a quadratic model of the objective function with a positive definite approximation of the Hessian of the associated Lagrangian. The novelty is that the problem functions of (24) are changed between iterations by replacing $\alpha_{i}(p,\hat{\phi})$, $c_{i}(p,\hat{\phi})$ and $H_{i}(p,\hat{\phi})$ by $\alpha_{i}^{s}(p,\hat{\phi}_{i}^{k})$, $c_{i}^{s}(p,\hat{\phi}_{i}^{k})$ and $H_{i}(p,\hat{\phi}_{i}^{k})$ before line 2.

There is some freedom in how the updates $\Delta\hat{\phi}_{i}^{k}$ are determined. One option is to perform a single step of a projected gradient method to maximize g_{i}(p^{k},\hat{\phi}_{i}^{k}+\,\mathchoice{{\hbox to4.99995pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{1.5}{\raisebox{0.0pt}{\displaystyle\cdot}}}}\hss}}}{{\hbox to4.99995pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{1.5}{\raisebox{0.0pt}{\textstyle\cdot}}}}\hss}}}{{\hbox to2.91666pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{1.5}{\raisebox{0.0pt}{\scriptstyle\cdot}}}}\hss}}}{{\hbox to1.80553pt{\hss\leavevmode\hbox{\set@color\raisebox{0.0pt}{\scalebox{1.5}{\raisebox{0.0pt}{\scriptscriptstyle\cdot}}}}\hss}}}\,). In this case, we choose $\hat{\phi}_{i}^{k+1}$ to be the projection of $\hat{\phi}_{i}^{k}+\sigma_{k}\nabla_{\phi}g_{i}(p^{k},\hat{\phi}_{i}^{k})$ onto the uncertainty set, where $\sigma_{k}>0$ is an adequate step length. For details, we refer to Section 2.2.2 of [20]. This strategy tries to increase the accuracy of the quadratic models by gradually moving the expansion points towards a local maximum of the modeled functions.

Since we change the MPEC formulation (24) in between two SQP iterations, standard convergence results do not apply. The major issue is that we cannot use the common monotonicity argument about the sequence of penalty function values, which is used to show that the penalty function values converge. A simple way to fix this is given in the next lemma.

Lemma 3.2**.**

Assume that the sequence of expansion points $(\hat{\Phi}^{k})_{k}$ converges as $k\to\infty$. Furthermore, assume that

[TABLE]

If the sequence $(x^{k})_{k}$ has an accumulation point, then

[TABLE]

Proof.

Let $x^{*}$ be an accumulation point of $(x^{k})_{k}$ and let $K\subseteq\mathbb{N}$ be an infinite subset such that $x^{k}\xrightarrow{K}x^{*}$. The continuity of $P_{\rho}^{1}$ implies that $(P_{\rho}^{1}(x^{k};\hat{\Phi}^{k}))_{k}$ is a Cauchy sequence. Hence, for any $\varepsilon>0$ there is some $N(\varepsilon)\in\mathbb{N}$ such that

[TABLE]

Noting $P_{\rho}^{1}(x^{k_{2}};\hat{\Phi}^{k_{2}})-P_{\rho}^{1}(x^{k_{1}};\hat{\Phi}^{k_{1}})=\sum_{k=k_{1}}^{k_{2}-1}P_{\rho}^{1}(x^{k+1};\hat{\Phi}^{k+1})-P_{\rho}^{1}(x^{k};\hat{\Phi}^{k})$, we obtain

[TABLE]

Since $P_{\rho}^{1}(x^{k+1};\hat{\Phi}^{k})-P_{\rho}^{1}(x^{k};\hat{\Phi}^{k})\leq 0$ by line 3 of Algorithm 1, the summability assumption shows that $0>\sum_{k=k_{1}}^{k_{2}-1}P_{\rho}^{1}(x^{k+1};\hat{\Phi}^{k})-P_{\rho}^{1}(x^{k};\hat{\Phi}^{k})>-\infty$ as $k_{2}\xrightarrow{K}\infty$. This implies $P_{\rho}^{1}(x^{k+1};\hat{\Phi}^{k})-P_{\rho}^{1}(x^{k};\hat{\Phi}^{k})\to 0$ for $k\to\infty$, from which the claim follows. ∎

The lemma can be used to apply standard convergence results to Algorithm 1. For example, we can argue analogously to the proof of Theorem 1 in [34] in order to show that every accumulation point of $(x^{k})_{k}$ is a KKT point of (24), under the additional, usual assumptions of SQP methods like uniform boundedness of the Hessian approximations of the Lagrangian. In this proof, the existence of a limit of the sequence $(P_{\rho}^{1}(x^{k};\hat{\Phi}^{k}))_{k}$, which is established by Lemma 3.2, implies that every accumulation point of $(x^{k})_{k}$ is feasible. The rest of the proof is simple to modify, so that we skip the details.

4 Proper Orthogonal Decomposition for Parametrized problems

Numerical approximation for robust optimization problems can be expensive since it involves the solution of several PDEs. Furthermore, the sensitivity approach enlarges the number of PDEs and it increases the computational costs of its approximation. For this reason, in this section, we introduce a model order reduction technique to reduce the complexity of the problem.

Here, we focus on the POD method for the approximate solution of the parametrized equation (11). As before, ${\bf u}_{h}(p)\in\mathbb{R}^{N}$ is the model vector associated to the FE solution of (11) for a given parameter $p\in\mathcal{D}\subset\mathbb{R}^{N_{p}}$ and ${\bf u}^{n}_{h,\mu}(p)\in\mathbb{R}^{N}$ are the sensitivities of order $n$ according to (14). For this purpose let $\{p^{j}\}_{j=1}^{m}$ be some points in $\mathcal{D}$ and let ${\bf u}_{h}(p^{j}),{\bf u}^{n}_{h,\mu}(p^{j})$ denote the corresponding solutions to (11), (14) for $p^{j}$. We define the snapshot set $Y\in\mathbb{R}^{n\times(mn_{max}+1})$

[TABLE]

including the states and the sensitivities up to order $n_{max}$ and determine a POD reduced space $\mathcal{V}^{\ell}:=\mbox{span}\{\psi_{1},\ldots,\psi_{\ell}\}$ of rank $\ell$ by solving the following minimization problem

[TABLE]

where we set ${\bf u}^{0}_{h,\mu}(p^{j})={\bf u}_{h}(p)$ with $|\mu|=0$ for brevity, $n_{max}$ is the maximum order of considered sensitivies, $\beta_{j}$ are nonnegative weights, $\delta_{ij}$ denotes the Kronecker symbol, ${\bf W}$ is a symmetric positive definite $N\times N$ matrix and $\psi_{i}\in\mathbb{R}^{N}$. The weighted inner product used is defined as follows: {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\langle{\bf u},{\bf v}\rangle_{\bf W}}:={\bf u}^{\top}{\bf W}{\bf v}.

It is well-known (see [16]) that problem (4) admits a unique solution $\{\psi_{1},\ldots,\psi_{\ell}\}$, where $\psi_{i}$ denotes the $i-$th eigenvector of the self-adjoint linear operator $\mathcal{R}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n},$ i.e., $\mathcal{R}\psi_{i}=\lambda_{i}\psi_{i}$ with $\lambda_{i}\in\mathbb{R}$ non-negative, where $\mathcal{R}$ is defined as

[TABLE]

Furthermore, the error in (4) can be expressed as

[TABLE]

where $d$ is the rank of the snapshots matrix $Y$.

4.1 POD approximation for state and sensitivities

We briefly recall how to generate the reduced order modeling by means of POD. Suppose we have computed the POD basis $\{\psi_{i},\ldots,\psi_{\ell}\}$ of rank $\ell$ according to the minimization problem (4). For the weight matrix we choose

[TABLE]

where $\bar{p}$ is a fixed reference parameter and $\bf M$ denotes the mass matrix. Then, ${\bf W}$ is the matrix associated to the discrete $H^{1}-$norm. We define the POD ansatz for the state as ${\bf u}_{h}^{\ell}(p):=\sum_{i=1}^{\ell}({\bf\bar{u}}^{\ell})_{i}\psi_{i}.$ This ansatz in (11) leads to an $\ell-$dimensional linear system for the unknowns $\{{\bf\bar{u}}_{i}\}_{i=1}^{\ell}$, namely

[TABLE]

Here the entries of the stiffness matrix are given by ${(\bf K^{\ell}})_{ij}=\psi_{j}^{\top}{\bf K}(p)\psi_{i}$ for $1\leq i,j\leq\ell$. The right hand side is given by ${(\bf f^{\ell}})_{i}=\psi_{i}^{\top}{\bf f}(p)$, $1\leq i\leq\ell$. We recall that due to the previously introduced affine decomposition this projection has to be computed only once and the system matrix can be written as

[TABLE]

where $({\bf K}^{q,\ell})_{ij}=\psi_{i}^{\top}{\bf K}^{q}\psi_{j}$ for $1\leq i,j\leq\ell$ and $q=1,\ldots,Q$. The same structure can be used for the right hand side. Note that this is very important in order to obtain an efficient reduced order model, since the system can be set up for different values of $p$ without the need of the original high dimensional matrices and right hand sides. The reduced order state equation reads:

[TABLE]

In an analogous way we obtain the general reduced sensitivity equation from (14). We need to make an ansatz for the sensitivities ${\bf u}^{n,\ell}_{h,\mu}$ and project the system onto the subspace spanned by the POD basis. In the present work we use the same basis functions for the state and the sensitivity variables and to achieve a better approximation property of the reduced order sensitivities we add the solution of the sensitivity equation to the snapshots set. Note that then the stiffness matrix in the reduced sensitivity equation is the same as in the reduced state equation, so that only the right-hand side needs to be projected.

4.2 A-posteriori error estimations

A-posteriori error estimators have a crucial role in model order reduction. They provide a certification of the surrogate model without the need of the computation of the truth solution. In the present work, we consider as truth solution the finite element approximation. If the transformation on a subdomain is given by (2), the bilinear form reads

[TABLE]

and a lower bound for the coercivity constant $\alpha(p)$ in (6) can be computed from the following estimate

[TABLE]

where {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\lambda_{min}} denotes the smallest eigenvalue of the matrix $\left(C_{q}(p)^{-1}\nu_{q}C_{q}(p)^{-T}\right)$ and $\bar{p}$ is the fixed reference parameter. An upper bound for the continuity constant $\gamma(p)$ in (7) can be computed with the same approach.

Then, we can derive an error bound for the reduced state and sensitivity equations (see [37] for more details) in terms of the reduced residual of the aforementioned equations.

For this purpose we define the residual for equation (14) as follows:

Definition 4.1**.**

Let u^{n,\ell}_{h,\mu}\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X_{h}(\bar{\Omega})}, $n=0,\ldots,n_{max}$, be the reduced order solution of (14). We define the residual

[TABLE]

Then, we have

Theorem 4.2**.**

Let u_{h,\mu}^{n}\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X_{h}((\bar{\Omega})} be the solution to (14) for $|\mu|=n$ and u^{n,\ell}_{h,\mu}\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X_{h}(\bar{\Omega})} be the corresponding reduced solution of (14). Then, the following inequality holds

[TABLE]

where $\gamma_{\kappa}$ is the continuity constant of the $\kappa$-th derivative of the coercive bilinear form.

Proof. We denote the error by $e_{u^{n}_{\mu}}:=u^{n}_{h,\mu}-u^{n,\ell}_{h,\mu}$. With (9) we find that

[TABLE]

Now we set $v=e_{u^{n}}$ and obtain

[TABLE]

By applying Cauchy-Schwarz and using the coercivity of $a$ as well as the continuity of $\frac{\partial^{\kappa}a}{\partial p^{\kappa}}$ we find

[TABLE]

Dividing by \alpha(p)\|e_{u^{n}_{h,\mu}}\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X_{h}(\bar{\Omega})}} leads to (31) first for $n=0$ and by inductively using the bound \|u_{h,\mu-\kappa}^{n-|\kappa|,\ell}-u_{h,\mu-\kappa}^{n-|\kappa|}\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X_{h}(\bar{\Omega})}}\leq\Delta^{\ell}_{u_{\mu-\kappa}^{n-|\kappa|}} for $n\leq n_{max}$. $\Box$

Remark 4.1**.**

We note that (31) is a generalization of well-known error bounds for state and first order sensitivity equations, compare e.g. [31], and [13] for time-dependent problems.

Remark 4.2**.**

If only the linear form $f$ depends on the parameter $p$, i.e., $a(u,v)=f(v;p)$, the error estimator (31) reads

[TABLE]

We note that the solution $u(\mu)$ in this case lives in a $Q_{f}$- dimensional linear subspace, where $Q_{f}$ is defined in (5), see e.g. [17]. Then, for $\ell\geq Q_{f}$ the error with the reduced model is zero and no error bound is needed. We note that sensitivity equations and error estimators for the uncertainties $\phi$ are analogous to (30)-(31) and are left to the reader.

4.3 The POD method for optimization problem

In this section we explain how to solve the parametrized optimization problem. In our construction the POD spaces depend on the parameter $p$. Since the solution to the parameter optimization problem is not known in advance the POD space has to be adapted/enriched during the parameter optimization procedure. To achieve this goal we here propose a certified extension of the approach suggested in [1] to our robust setting, where the error bound, $\Delta_{u_{\mu}^{n}}^{\ell}(p)$ introduced in Section 4.2, helps in the selection of the snapshot sets. The algorithm works as follows: we start with a very coarse parameter sample choosing only one parameter $p^{0}$ and solve the full problem together with the sensitivity equations associated to this parameter. Then, we compute the POD basis functions and perform the reduced optimization procedure. At the end of the process we find a new parameter $p^{1}$ which is an approximation of the optimal desired design, we update the parameter set $\mathcal{D}=\{p^{0},p^{1}\}$, solve the full problem and the sensitivity equations related to the new parameter $p^{1}$. Then, we enlarge the snapshots set and compute new POD basis functions. We iterate this process until the stopping criterion is reached. The procedure is summarized in Algorithm 2.

In our simulations this approach turned out to be very efficient since it avoids long pre-computations. In this way we are able to update the snapshot set and the POD basis functions. Our update contains information on the optimization problem and it improves the quality of our surrogate model. Note that every reduced optimization problem contains the a-posteriori error for the state and sensitivity equations introduced in Section 4.2.

5 Numerical tests

In our numerical example we consider an optimal design problem for a permanent magnet synchronous machine. We start by introducing the model and geometry under consideration. We consider a three-phase six-pole permanent magnet synchronous machine (PMSM) with one buried permanent magnet per pole. The geometry is shown in Figure 1. The goal of the design optimization is to change the size and location of the permanent magnet such that the material of the magnet is minimized while maintaining a lower bound on the electromotive force. We consider a description using three parameters: $p_{1}$ the width, $p_{2}$ the height and $p_{3}$ the central perpendicular distance between the rotor and the surface of the magnet. The region around the permanent magnet (Figure 1 red box) is decomposed into twelve triangles (Figure 1 blue lines). The introduced triangulation of the parametrized domain allows to perform the affine linear decomposition as introduced in (12) with $Q_{a}=12=Q_{f}$.

PMSMs can be described sufficiently accurate by the magneto-static approximation of Maxwell’s equation. The governing parametrized equation is given by

[TABLE]

with boundary conditions

[TABLE]

where $\nu$ is the reluctivity, $J_{src}$ is the source current density and $H_{pm}$ the field of the permanent magnets (PM). We note that equation (32) fits into the abstract formulation (1). In the 2D planar setting together with a finite element method for the magnetic vector potential, we lead to the discrete form given by the linear model presented in (11). To extract performance values the loading method is used to exploit the frequency domain [36]. To obtain quantities like the electromotive force (EMF) a Fourier analysis of the magnetic vector potential around the inner surface of the stator is carried out. This can be written as a linear function $E_{0}=\mathbb{E}^{\top}{\bf u}_{h}$. More details on the configuration we adopt in the present work can be found in [2, 21, 32].

Let us next formulate the optimization problem. We start by introducing the nominal optimization problem. The goal of the optimization is to minimize the required material for the permanent magnet while maintaining a lower bound on the EMF $E_{0}$. In the mathematical model this leads to a cost function of the form

[TABLE]

where $E_{0}^{d}$ is the desired lower bound on the EMF and $\rho\in\mathbb{R}^{+}$ a weight parameter. Additionally, we have the constraints

[TABLE]

The upper and lower bounds for the parameters and the first inequality are due to the parametrization of the geometry and the restriction that the permanent magnet has to stay within the red box in Figure 1 (left). The last inequality is a design restriction that avoids that the corner of the permanent magnet comes too close to the rotor surface. Note that we right away use the reduced formulation as introduced in (16). To obtain a smooth formulation we introduce a slack variable. We reformulate the nominal optimization problem by using the variable $x=(p,\xi)$ as

[TABLE]

In this form the optimization problem fits exactly into the framework of (16).

Next let us introduce the uncertainty. For our numerical example we assume that the magnetic field angle in the permanent magnet is uncertain [29]. This can be due to manufacturing imprecision. In the nominal optimization the magnetic field is aligned perfectly, i.e. the field angle is $90$, see Figure 1 (right plot). In practice this can not be met and a deviation is to be expected. The field angle enters the model problem (32) nonlinearly through the right hand side, in particular in the term $H_{pm}$. In our numerical example we allow a field angle $\phi$ in the range $[85,95]$. Following the definition of the uncertainty set in Section 3.2 we get

[TABLE]

with $k\in\{2,\infty\}$, i.e. \hat{\phi}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}90} and $D=5$ in (17). Using these settings we can now solve the linear and quadratic approximation of the robust counterpart (21) and (24). In the case of the linear approximation we choose $k=\infty$ and for the quadratic approximation we set $k=2$. Note that in the presented setting these two uncertainties sets are identical. The different choices are for technical reason as outlined in the derivation.

Before presenting the numerical results let us give a short overview of the numerical strategy utilized to solve the optimization problems. The computations are carried out in MATLAB. To solve the nominal and robust optimization problems an SQP method with Armijo-backtracking strategy using a $\ell_{1}$-penalty function is used [28]. The Hessian is computed via BFGS updates. Alternatively also routines like fmincon in MATLAB can be used obtaining similar results. The derivative of $g_{i},\,i=0,\ldots,8$ are computed using the sensitivity approach [20]. Also the derivatives with respect to the uncertain parameter $\phi$ are computed using this approach. The structure of the sensitivity equations are as outlined in (13).

5.1 Results obtained by the finite element approximation

We start by presenting the numerical results utilizing the finite element approximation. Piecewise linear and continuous finite elements are used to discretize equation (32) leading to a system with $61013$ degrees of freedom.

The initial geometry configuration corresponding to $p=(19,7,7)$ is shown in Figure 1 (left) together with the corresponding magnetic vector potential (right). From the magnetic vector potential we extract the EMF which we will use as the desired value $E_{0}^{d}=30.34$ in our optimization problem.

In Table 1 we show the results obtained in the optimization. In column $V_{pm}$ the volume of the permanent magnets is given. Note that through the optimization a significant reduction in size is achieved. The ratio is given in percent in the second column. For the nominal optimization a reduction by $53\%$ is obtained and in the robust case reductions of $50\%$ (linear case) and $49\%$ (quadratic case). The respective parameters are given in the third column. Lastly, in column four and five the EMF $E_{0}$ for the field angle $\phi=90$ and the worst-case EMF $E_{0}^{worst}$ are given. It can be seen that the uncertainty in the magnetic field of the permanent magnet has an impact on the performance. In the robust optimization this influence is incorporated. Hence it can be seen that in the case of the quadratic approximation (Rob. Quad and Rob. Shift) very good results can be achieved: the worst-case EMF stays above the target value of $30.34$. Shifting the expansion point (Rob. Shift) helps to avoid being overly conservative, so that the EMF target is precisely met, which allows to reduce the volume slightly. In the case of the nominal optimization a significant decay can be observed. In the case of the linear approximation the target can not be met since the approximation is not accurate enough. It can also be observed that by performing only a nominal optimization the worst-case can decrease compared to the initial configuration. In the presented case the difference is small but can become more significant in different settings.

The approximation quality and the behavior of the EMF is shown in Figure 2 for the different optimal designs. In the left plot the behavior of the EMF for different permanent magnet field angles $\phi$ is shown for the initial and the nominal optimal configuration. It can be seen that the target $E_{0}^{d}$ is only reached in $\phi=90$. In the middle plot the linear approximation and the actual EMF are compared. It can be seen that the approximation is not accurate enough to determine the worst-case. The right plot corresponds to the quadratic approximation. It can be seen that the EMF is approximated very well by the quadratic model. Hence also good results in the robust optimization using this approximation can be expected. We note that we do not include a plot for the quadratic approximation with a shifted expansion point because the difference is not visible from the standard quadratic approximation.

In Figure 3 the magnetic vector potentials are shown for the three optimal designs. To conclude we have a look at the computational expenses for the optimization. These are summarized in Table 2, where the computational time in seconds, the number of iterations and the number of required PDE solves are compared. While the computational time in this example is still very low it can be seen that the robust optimization is more expensive. Especially the number of PDE solves increases significantly compared to the nominal optimization. Hence we will use model order reduction to reduce the computational costs. This will be outlined in the next section.

5.2 Results obtained by the reduced order model

Let us start our analysis of the reduced order model by dealing with the approximation of the state and the sensitivity variables. To this aim we select $125$ parameters chosen as follows:

[TABLE]

and we compute an enlarged snapshot set with state and sensitivities for each parameter in $\mathcal{D}_{train}$. Altogether, we collect $500$ snapshots. We note that we do not need to compute a reduced order model with respect to the slack variable $\xi$ since it is not a variable that is involved in the PDE.

In Figure 4, we present an error analysis to check the quality of our a-posteriori error discussed in Section 4. We compare the error between the POD solution and the high dimensional approximation for state and sensitivities with the error bound presented in Section 4.2. With the W-norm introduced in Section 4, we set

[TABLE]

where $\mathcal{D}_{test}$ is chosen in the centers of the boxes of $\mathcal{D}_{train}$ as follows:

[TABLE]

As one can see model order reduction is able to reach only an error of size $10^{-4}$ (with a large amount of basis functions) and we would need to increase the number of basis functions to achieve higher accuracy. We also want to emphasize that a standard POD approach, with only state snapshots in $\mathcal{D}_{train}$, is able to reach an accuracy of order $10^{-2}$ with $50$ basis functions.

In Figure 5, we show that it is possible to reach an accuracy of order $10^{-6}$ or higher with only 4 POD basis functions if we compute the snapshots with respect to one parameter in $\mathcal{D}_{train}$ and then compute the error in a ball of radius 0.1, centered in the chosen parameter. This is our motivation to introduce Algorithm 1 for the successive enrichment of the POD model. In fact, we start our algorithm with only one parameter and we require the combined snapshots since in any other case we will not have enough data to generate a surrogate model from the simulation and the sensitivities. We note that if we add the sensitivities, the POD basis functions are improved and the parameter domain can be better explored. We refer the interested reader to [19] for more details.

Let us now draw our attention to the optimization problem and its performances when combined with model order reduction. Table 3 shows the convergence of Algorithm 2 for the nominal optimization. The number of updates of the snapshot set is given in the first column. Since our goal is to reach a desired electromotive force, we will also compare $E_{0}(p,{\bf u}_{h,\mu}^{0,\ell}(p),90)$ obtained from the POD, and $E_{0}(p,{\bf u}_{h,\mu}^{0}(p),90)$ obtained from the full simulation. As one can see with the surrogate-based optimization obtain the same results as with the full model. The fifth column presents the number of iterations needed in each sub-optimization problem and the sixth column presents the worst case of the error estimate $\Delta_{u_{\mu}^{n}}^{\ell}(p)$ for $0\leq|\mu|=n\leq 1$. In the stopping criterium of Algorithm 2 we choose $tol=1e-4$. In the last column we show the number of basis functions which also represents the dimension of the reduced problem. The dimension of the reduced space corresponds to the rank of current snapshot matrix. We note that the results of the algorithm are close to what we have shown in Table 2. The CPU time will be discussed at the end of the section.

Next we analyze the performance of POD for the robust optimization problem. In Table 4 we present the results of POD with a linear approximation of the robust counterpart with $n_{max}=2$. As one can see after two iterations we are able to recover the same results obtained without model order reduction, obtaining also a good approximation for the electromotive force (compare column 2 and column 3 in Table 4).

Then, we present the results of the surrogate models for the solution of the robust optimization problem with quadratic approximation with $n_{max}=3$ in Table 5 and with quadratic approximation with moving expansion point in Table 6. The same considerations discussed in the linear example hold true. We note that the number of POD basis functions used in each optimization in this case is larger due to a richer snapshot set.

Also in the quadratic case we obtain the same approximation quality as in the full model at lower computational costs (see Table 5 and Table 6).

To summarize, all the performances of the surrogate models are presented in Table 7. Again, we emphasize the results agree with Table 1. Strictly speaking, none of the solutions obtained with reduced models are robust for the high-fidelity problem. This is demonstrated by the last column of Table 7. All of the numbers should be, but none of them is greater than or equal to $30.34$. The relative worst-case infeasibility violation ranges from $1.9\,\%$ for the nominal solution to less than $0.07\,\%$ for the improved quadratic model. Seeing that the reduced models are able to reduce the volume more than the high-fidelity models, this is to be expected. In many applications, an infeasibility as low as $0.07\,\%$ in the worst case presumably is acceptable; if not, a post-optimization analysis is inevitable.

The computational costs of the POD approach are summarized in Table 8. As one can see we obtain a reduction of the amount of PDEs solved (compare Table 2), keeping the same accuracy attained in the high dimensional problem. For instance, in the quadratic robust optimization with model order reduction we only have to perform $36$ PDE solves in contrast to $72$ PDE solves in the optimization with the high-fidelity finite element model. The POD approach for quadratic approximation with shift presents a speed up of factor $5$. Although we solve more (reduced) PDEs the CPU time is reduced due to the number of full dimensional PDEs to be solved in the full setting, compare Tables 2 and 8. Finally, we note that the CPU time involves both offline and online stage, so that the speed up is considerable.

6 Conclusion

In this paper we present a new approach which combines model order reduction to robust optimal control. We investigate a parameter optimization problem governed by a parametric elliptic partial differential equation with uncertain parameters. We introduce a robust optimization framework that accounts for uncertain model and optimization parameters. By shifting the expansion point of the utilized Taylor approximation, we are able to model the nonlinear effects of the uncertain parameter more accurately and improve the robust solutions slightly. The resulting non-linear optimization problem has a bi-level structure due to the min-max formulation. We propose an adaptive model order reduction approach which avoids long offline stages and provides a certified reduced order surrogate model for the parametrized PDE which then is utilized in the numerical optimization. The presented numerical results clearly illustrate the validity and performance of the presented approach.

7 Acknowledgements

The authors wish to acknowledge the support of the German BMBF in the context of the SIMUROM project (grant number 05M2013), and the support of the German Research Foundation in the context of SFB 805.