Optimal control for multiscale equations with rough coefficients

TL;DR

This paper presents a novel numerical homogenization approach using Rough Polyharmonic Splines for optimal control problems governed by multiscale elliptic equations with rough coefficients, achieving efficient solutions without regularity or scale-separation assumptions.

Contribution

It introduces the GRPS method for multiscale elliptic equations in optimal control, enabling efficient repeated solutions with minimal fine-scale computations.

Findings

Optimal convergence rates independent of coefficient regularity

One-time fine-scale pre-computation reduces iterative costs

Validated through numerical experiments

Abstract

This paper concerns the convex optimal control problem governed by multiscale elliptic equations with arbitrarily rough $L^\infty$ coefficients, which has important applications in composite materials and geophysics. We use one of the recently developed numerical homogenization techniques, the so-called Rough Polyharmonic Splines (RPS) and its generalization (GRPS) for the efficient resolution of the elliptic operator on the coarse scale. Those methods have optimal convergence rate which do not rely on the regularity of the coefficients nor the concepts of scale-separation or ergodicity. As the iterative solution of the OCP-OPT formulation of the optimal control problem requires solving the corresponding (state and co-state) multiscale elliptic equations many times with different right hand sides, numerical homgogenization approach only requires one-time pre-computation on the fine scale and the following iterations can be done with computational cost proportional to coarse degrees of freedom. Numerical experiments are presented to validate the theoretical analysis.

Equations

$$\begin{array}[]{cl}&\min\limits_{u\in K\subset L^{2}(\Omega_{U})}g(y)+h(u)\\ \text{subject to}&\left\{\begin{array}[]{cc}-\operatorname{div}(a(x)\nabla y)=f+Bu,&\text{ in }\Omega,\\ y=0,&\text{ on }\partial\Omega.\end{array}\right.\end{array}$$

$$(u,v{)}_Y$$

$$(u,v{)}_U$$

$$a_{\min }|\xi {|}^2\le {\xi }^{T}a(x)\xi \le a_{\max }|\xi {|}^2,\forall \xi \in R^d/\{0\},x\in \bar{\Omega }.$$

$$\Vert g^{\prime }(u)-g^{\prime }(v){\Vert }_Y\le L_g\Vert u-v{\Vert }_Y,\forall u,v\in Y,$$

$$\Vert h^{\prime }(u)-h^{\prime }(v){\Vert }_U\le L_h\Vert u-v{\Vert }_U,\forall u,v\in U,$$

$$(g^{\prime }(u)-g^{\prime }(\tilde{u}),u-\tilde{u}{)}_Y\ge M_g\Vert u-\tilde{u}{\Vert }_Y^2,\forall u,\tilde{u}\in Y,$$

$$(h^{\prime }(u)-h^{\prime }(\tilde{u}),u-\tilde{u}{)}_U\ge M_h\Vert u-\tilde{u}{\Vert }_U^2,\forall u,\tilde{u}\in U,$$

$$a(y,v):={\int }_{\Omega }a(x)\nabla y\cdot \nabla v,\text{ for }y,v\in V.$$

$$\begin{array}[]{cc}&\min\limits_{u\in K\subset U}g(y)+h(u)\\ &\text{subject to}\quad a(y,v)=(f+Bu,v),\forall v\in V.\end{array}$$

$$\begin{array}[]{cc}a(y,v)=(f+Bu,v)_{Y},&\forall v\in V,\\ a(p,q)=(g^{\prime}(y),q)_{Y},&\forall q\in V,\\ (h^{\prime}(u)+B^{*}p,\tilde{u}-u)_{U}\geq 0,&\forall\tilde{u}\in K\subset U.\end{array}$$

$$a(z,w)=(\rho ,w),\text{for any }w\in V,$$

$$\begin{array}[]{cc}&\min\limits_{u_{h}\in K_{h}\subset U_{h}}g(y_{h})+h(u_{h})\\ \text{subject to}&a(y_{h},v_{h})=(f+Bu_{h},v_{h}),\forall v_{h}\in V_{h}.\end{array}$$

$$\begin{array}[]{cc}a(y_{h},v_{h})=(f+Bu_{h},v_{h}),&\forall v_{h}\in V_{h},\\ a(p_{h},q_{h})=(g^{\prime}(y_{h}),q_{h}),&\forall q_{h}\in V_{h},\\ (h^{\prime}(u_{h})+B^{*}p_{h},\tilde{u}_{h}-u_{h})_{U}\geq 0,&\forall\tilde{u}_{h}\in K_{h}\subset U_{h}.\end{array}$$

$$\Vert y-y_h{\Vert }_{1,\Omega }+\Vert p-p_h{\Vert }_{1,\Omega }+\Vert u-u_h{\Vert }_{0,{\Omega }_U}\le C(h_U+h)(\Vert y{\Vert }_2+\Vert p{\Vert }_2+\Vert u{\Vert }_1).$$

$$\{\begin{array}{rl}-\nabla \cdot (a(x,x/\epsilon )\nabla z^{\epsilon })=\rho ,& x\in \Omega ,\\ z^{\epsilon }=0,& x\in \partial \Omega .\end{array}$$

$$\Vert y-y_h{\Vert }_{1,\Omega }+\Vert p-p_h{\Vert }_{1,\Omega }+\Vert u-u_h{\Vert }_{0,{\Omega }_U}\le C(h+\sqrt{\epsilon /h}).$$

$$\{\begin{array}{rl}-\nabla \cdot (a(x)\nabla z)=\rho (x),& x\in \Omega ,\\ z=0,& x\in \partial \Omega .\end{array}$$

$$V_i:=\{\varphi \in V|{\int }_{\Omega }\varphi (x){\psi }_j(x)\text{dx}={\delta }_{i,j},j=1,\dots ,N\}.$$

$$\{\begin{array}{l}\text{Minimize }\Vert \varphi {\Vert }^2\\ \text{Subject to }\varphi \in V_i\end{array}$$

$$\left\{\begin{array}[]{lr}-\text{div}(a\nabla z)=\rho(x)&x\in\Omega\\ z=0&x\in\partial\Omega\\ \end{array}\right.$$

$$\Gamma (x,y)={\int }_{{\Omega }^2}G(x,z)\Lambda (z,z^{\prime })G(y,z^{\prime })\text{dz}{\text{dz}}^{\prime }.$$

$$\mathbb{E}[z(x)|\Psi ]=\sum _{i=1}^N{\Psi }_i{\varphi }_i(x),\text{ where }{\varphi }_i(x):=\sum _{j=1}^N{\Theta }_{ij}^{-1}{\int }_{\Omega }\Gamma (x,y){\psi }_j(y)\text{dy}$$

$$\Vert z{\Vert }^2=\Vert z_{\mathrm{I}}{\Vert }^2+\Vert z-z_{\mathrm{I}}{\Vert }^2.$$

$$\Vert \nabla z-\nabla z_{\mathrm{I}}{\Vert }_{L^2(\Omega )}\le C_{\mathrm{I}}H\Vert \rho {\Vert }_{L^2(\Omega )}.$$

$$\Vert z-z_H{\Vert }_{H_0^1(\Omega )}\le C_{\mathrm{NH}}^{g}H\Vert \rho {\Vert }_{L^2(\Omega )}.$$

$${\int }_{\Omega \cap (B(x_i,r){)}^c}\nabla {\varphi }_i^{t}a\nabla {\varphi }_i\le \exp (1-\frac{r}{CH}){\int }_{\Omega }\nabla {\varphi }_i^{t}a\nabla {\varphi }_i.$$

$$V_i^{\text{loc}}:=\{\varphi \in V|{\int }_{{\Omega }_i}\varphi (x){\psi }_j(x)\text{dx}={\delta }_{i,j},j=1,\dots ,N\}.$$

$$\{\begin{array}{l}\text{Minimize }\Vert \varphi {\Vert }_{{\Omega }_i}^2\\ \text{Subject to }\varphi \in V_i^{\text{loc}}\end{array}$$

$$\Vert {\varphi }_i-{\varphi }_i^{\mathrm{loc}}{\Vert }_a\le C\exp (-\frac{r}{2lH}).$$

$$\Vert z-z_H^{\text{loc}}{\Vert }_{H_0^1(\Omega )}\le C_{\mathrm{NH}}H\Vert \rho {\Vert }_{L^2(\Omega )}.$$

$$\begin{array}[]{cc}a(y_{H},v_{H})=(f+Bu_{H},v_{H}),&\forall v_{H}\in V_{H},\\ a(p_{H},q_{H})=(g^{\prime}(y_{H}),q_{H}),&\forall q_{H}\in V_{H},\\ (h^{\prime}(u_{H})+B^{*}p_{H},\tilde{u}_{H}-u_{H})_{U}\geq 0,&\forall\tilde{u}_{H}\in K_{H}\subset U_{H}.\end{array}$$

Full text

{AMS}

49J20, 65N15, 65N30, 74Q05

1 Introduction

We consider the following convex optimal control problem (OCP) governed by elliptic partial differential equations with rough coefficients $a(x)\in L^{\infty}(\Omega)$

[TABLE]

where $\Omega$ ($\Omega_{U}$) is a bounded convex polygon in $\mathbb{R}^{d}$ ($d=2,3$) with Lipschitz boundary $\partial\Omega$ ($\partial\Omega_{U}$). We adopt the standard notation $W^{m,p}$ for Sobolev spaces on $\Omega$ with norm $\|\cdot\|_{m,p}$ and seminorm $|\cdot|_{m,p}$ [1]. We define $W^{m,p}_{0}:=\{\psi\in W^{m,p}:|\psi|_{\partial\Omega}=0\}$. For $p=2$, we denote $H^{m}(\Omega)=W^{m,2}(\Omega)$, and $\|\cdot\|_{m}=\|\cdot\|_{m,2}$. Take $\|\cdot\|=\|\cdot\|_{0,2}$. In the rest of the paper, we take $V=H^{1}_{0}(\Omega)$ as the state space, $U=L^{2}(\Omega_{U})$ as the control space, and $Y=L^{2}(\Omega)$ as the observation space.

The admissible set $K\subset U$ for the control variable $u$ is a closed convex set, for example, we can take $K=\{u\in L^{2}(\Omega_{U})|u\geq 0\}$, or $K=\{u\in L^{2}(\Omega_{U})|\int_{\Omega_{U}}u\geq 0\}$, or $K=\{u\in L^{2}(\Omega_{U})|a\leq u\leq b\}$. $y\in V$ is the state variable. $g(y)+h(u)$ is the objective functional. $g(\cdot)$, $h(\cdot)$ are two convex differentiable functions on $Y$ and $U$, respectively. $B:U\to L^{2}(\Omega)$ is a continuous linear operator such that $\|Bw\|_{Y}\leq C_{B}\|w\|_{U}$ for any $w\in U$, for example, $(Bu)(x)=c(x)u(x)$, where $c(x)$ is a density factor. Assume that the forcing term $f(x)\in L^{2}(\Omega)$.

Define the following inner products,

[TABLE]

Denote $\|\cdot\|_{Y}$ and $\|\cdot\|_{U}$ as the corresponding norms induced by $(\cdot,\cdot)_{Y}$ and $(\cdot,\cdot)_{U}$, respectively. We may drop the subscript $Y$ and $U$ if no confusion arises.

We impose the following assumptions on the coefficients and data.

(A1) The coefficients matrix $a(x)=(a_{ij}(x))$ is a symmetric $d\times d$ matrix which satisfies the uniformly elliptic condition, i.e.,

[TABLE]

where $a_{\min}$ and $a_{\max}$ are positive constants. $\kappa=a_{\max}/a_{\min}$ is the contrast.

(A2) Let $g^{\prime}(\cdot)$ be the Gâteaux derivative of the functional $g(\cdot)$ on $Y$, $g^{\prime}\in Y^{\prime}=Y$. Similarly, let $h^{\prime}(\cdot)$ be the Gâteaux derivative of the functional $h(\cdot)$ on $U$, $h^{\prime}\in U^{\prime}=U$. Assume that $g^{\prime}$ and $h^{\prime}$ are Lipschitz continuous, i.e.,

[TABLE]

where $L_{g}$ and $L_{h}$ are positive constants.

(A3) The functionals $g(\cdot)$ and $h(\cdot)$ are uniformly convex,

[TABLE]

where $M_{g}$ and $M_{h}$ are positive constants.

Optimal control plays an increasingly important role in many engineering branches, and efficient numerical methods are essential to its successful application [23, 40]. Over the past 30 years, finite element method (FEM) has become one of the most widely used numerical methods for optimal control problems. For the optimal control of elliptic or parabolic equations: a priori error estimates [23, 22, 45], a posteriori error estimates [31, 4], and some superconvergence results [8, 9] for the FEM methods have been developed in the literature. However, those error estimates require the $H^{2}$ regularity of the solutions.

The optimal control problems governed by partial differential equations with rough coefficients (such as permeabilities in reservoir modelling) have become a great challenge, owing to the lack of regularity of the coefficients $a(x)$ and therefore the solutions $y(x)$ and $p(x)$ ($p(x)$ is the co-state variable in (8)). Even if $a(x)$ is a smooth but highly oscillatory function in the form of $\displaystyle a(x,\frac{x}{\varepsilon})$ with a small parameter $\varepsilon\ll 1$ (such as material properties of composite materials), conventional FEMs based on piecewise polynomial basis will not be effective [24]. To be more precise, the convergence of conventional FEMs relies on the $H^{2}$-regularity of the solution, but the prefactor of the error is of the order $O(1/\varepsilon)$. Thus, conventional FEMs require prohibitively small mesh size $h<\varepsilon$ to yield good numerical approximations for the optimal control problem.

Numerical homogenization for problems with multiple scales have attracted increasing attention in recent years. If the coefficient $a(x)$ has structural properties such as scale separation and periodicity, together with some regularity assumptions (e.g., $a(x)\in W^{1,\infty}$), classical homogenization [29, 26] can be used to construct efficient multiscale computational methods and have been applied to optimal control problems, such as multiscale asymptotic expansions method[6, 7, 30], multiscale finite element method (MsFEM) [24, 12, 10, 11], and heterogeneous multiscale method (HMM) [42, 21].

For multi-scale PDEs with non-separable scales and high-contrast coefficients which appear in many applications such as reservoir modelling and damage in composite materials, the coefficients do not have structural properties such as periodicity/scale separation. Numerical homogenization with non-separable scales concerns approximation of the solution space of such problems by a (coarse) finite dimensional space, instead of focusing on the classical issue of the homogenization limit. Fundamental questions for numerical homogenization are: How to approximate the high dimensional solution space by a low dimensional approximation space with optimal error control, and furthermore, how to construct the approximation space efficiently, for example, whether its basis can be localized on a coarse patch of size $O(H)$. Several novel approaches for numerical homogenization and their rigorous error analysis have been developed recently, such as: rough polyharmonic splines (RPS) method [39, 35, 36, 33]; Localized Orthogonal Decomposition (LOD) method [34, 19]; and Generalized multiscale finite element method (GMsFEM) [44, 20].

In the context of optimal control, homogenization based methods have been applied to problems governed by multiscale PDEs with separable scales [10, 30, 11]. To the best of our knowledge, few literature concerns the optimal control with non-separable scales, which is of great importance for applications.

The purpose of this work is to obtain the convergence result for the solution of optimal control problems governed by multiscale elliptic equations with rough coefficients using the so-called generalized rough polyharmonic splines (GRPS) method [33]. Those optimal control problems often arise in the optimal design of composite materials and the control of water injection in reservoir simulation. The GRPS approximation space is generated by an interpolation basis minimizing an appropriate energy norm subject to certain constrains. The resulting approximation space leads to a quasi-optimal $H^{1}$-accurate approximation of the solution space together with quasi-optimal localization properties. The GRPS approximation can be cast as a Bayesian inference problem under partial information [35].

The paper is organized as follows: We introduce the optimal control problem in Section § 2. We formulate the GRPS method and show its numerical properties in Section § 3. In Section § 4, we present the convergence analysis for the solution of optimal control problem using GRPS method. Numerical algorithms and results are given in Section § 5 to complement and justify the theoretical analysis. We conclude the paper in § 6.

Notation

For a finite set $A$, we will use $\#A$ to denote the cardinality of $A$. For a measurable set $\tau$, we use $|\tau|$ to denote its measure.

The symbol $C$ denotes generic positive constant that may change from one line of an estimate to the next. The dependencies of $C$ will normally be clear from the context or stated explicitly.

2 Formulation of the Optimal Control Problem

In this section, we present the equivalent formulations of the optimal control problem (1), as well as the corresponding finite element formulations and the error estimates.

Recall that $V=H^{1}_{0}(\Omega)$ is the state space, and $U=L^{2}(\Omega_{U})$ is the control space, we define the following energy product,

[TABLE]

The weak formulation of the optimal control problem (OCP) (1) is: find $(y,u)\in V\times U$ such that

[TABLE]

It is well know that (see, e.g, [23, Theorem 1.46] and [28]) the convex optimal control problem (1) has a unique solution $(y,u)$, and that $(y,u)$ is the solution of (1) if and only if there exists a co-state $p\in V$ such that the triple $(y,u,p)$ satisfies the following optimality conditions (OCP-OPT)

[TABLE]

where $B^{*}$ is the adjoint operator of $B$. In (8), the first equation is satisfied by the state function $y$ in (1), the second equation is satisfied by the co-state function $p$, and the third one is a variational inequality satisfied by the control function $u$.

Remark 2.1**.**

For some $K$ the variational inequality in (8) admits an analytical solution. One such as example is the admissible set $K=\{u\in L^{2}(\Omega_{U})|\int_{\Omega_{U}}u\geq 0\}$, the variational inequality in (8) is equivalent to $h^{\prime}(u)=-B^{*}p+\max\{0,\overline{B^{*}p}\}$, where $\bar{p}=\int_{\Omega_{U}}p/\int_{\Omega_{U}}1$ is the integral average of $p$ on $\Omega_{U}$. Another example is $K=\{u\in L^{2}(\Omega_{U})|a\leq u\leq b\}$, and the variational inequality is equivalent to $h^{\prime}(u)=\max\{a,\min(b,-B^{*}p)\}$.

Both $y$ and $p$ are solutions of the elliptic problem of the following form,

[TABLE]

where $\rho\in Y=L^{2}(\Omega)$.

The finite element approximation ${\rm(OCP)}_{h}$ for the optimal control problem (1) can be obtained by the restriction of $U$ and $V$ to their finite dimensional subspaces $U_{h}$ and $V_{h}$, respectively: find $(y_{h},u_{h})\in V_{h}\times U_{h}$ such that

[TABLE]

Again, the discretized control problem in (10) has a unique solution $(y_{h},u_{h})$ if and only if there is a co-state $p_{h}\in V_{h}$ such that the triplet $(y_{h},u_{h},p_{h})$ satisfies the following optimality conditions $\textrm{ (OCP-OPT)}_{h}$ [32]:

[TABLE]

Let $\mathcal{T}^{h}$ ($\mathcal{T}^{h}_{U}$) be a quasi-uniform triangulation of $\Omega$ ($\Omega_{U}$). Note that the regularities of the control variable $u$ is lower than the regularity of the state variable $y$ and the co-state variable $p$. We can choose the piecewise constant finite element space over $\mathcal{T}^{h}_{U}$ as $U_{h}$, and $K_{h}:=K\cap U_{h}$. If weak solutions of elliptic problem (9) admit $H^{2}$ regularity, namely, $\|z\|_{H^{2}}\leq C\|\rho\|_{Y}$, $V_{h}$ can be taken as the conforming piecewise linear finite element space over $\mathcal{T}^{h}$. We have the following theorem for the error estimates of $(y_{h},u_{h},p_{h})$ [32, Theorem 4.1.1].

Theorem 1**.**

Let $(y,u,p)$ be the solution of (8), and $(y_{h},u_{h},p_{h})$ be the finite element solution of (11). Assume that $u\in H^{1}(\Omega_{U})$, $y,p\in H^{2}(\Omega)$, it holds true that

[TABLE]

For optimal control governed by multiscale equations, we assume $a(x)=a(x,x/\varepsilon)\in W^{1,\infty}$ in the scale separation case. The Dirichlet problem

[TABLE]

has the following regularity estimate $\displaystyle\|z^{\varepsilon}\|_{2}\leq\frac{C}{\varepsilon}\|\rho\|$ for $\rho\in L^{2}(\Omega)$. Therefore, the error bound in Theorem 12 becomes $\displaystyle O(\frac{h_{U}+h}{\varepsilon})$, and one may need to take extremely small $h\ll\varepsilon$ and $h_{U}\ll\varepsilon$ in order to obtain accurate solution with conventional piecewise polynomial FEM spaces, and the computational costs will become prohibitive.

Multiscale numerical methods such as multiscale finite element method (MsFEM) [24, 25] can be used to efficiently capture the large scale components of the solution on the coarse grid. In [11, Theorem 4.11] and [27], MsFEM was applied to solve the control problem governed by multiscale PDEs with separable scales, a priori error estimates can be obtained as follows,

[TABLE]

However, numerical homogenization methods based on concepts such as scale separation and periodicity/ergodicity cannot be applied directly to problems with nonseparable scales. In the next section, we will introduce the concept of numerical homogenization that does not rely on the classical assumptions of homogenization theory such as scale separation and ergodicity, but only on the compactness of the solution space [38].

3 Numerical Homogenization and Generalized Rough Polyharmonic Splines

It is clear that a good approximation space for the elliptic equation (14) is important for the accurate solution of the optimal control problems (1) with arbitrarily rough coefficients $a(x)\in L^{\infty}$,

[TABLE]

This is closely linked to the problem of numerical homogenization with non-separable scales. Up to now, numerical homogenization has become a large field. It is motivated by the fact that standard methods, such as finite-element method with piecewise linear elements [3] can perform arbitrarily badly for PDEs with rough coefficients. Although some numerical homogenization methods such as multiscale finite element methods [24, 43, 25, 17], heterogeneous multiscale methods [41, 18] are directly inspired from classical homogenization concepts such as periodic homogenization and scale separation [29], one of the main objectives of numerical homogenization is to achieve a numerical approximation of the solution space of (14) with arbitrary rough coefficients (i.e., in particular, without the assumptions found in classical homogenization, such as scale separation, ergodicity at fine scales and $\varepsilon$-sequences of operators). Methods such as harmonic coordinates [37], generalized multiscale finite element [44, 14], flux-norm based approaches [5] have been proposed for numerical homogenization with arbitrary rough coefficients.

We first formulate the problem of numerical homogenization. Suppose $\varepsilon$ is the smallest scale of the elliptic problem, let $H$ be an artificial scale determined by available computational power and/or desired precision, $N$ is the corresponding dof such that $H\sim N^{-1/d}$, and the scales $\varepsilon\ll H\ll 1$. The goal of numerical homogenization is to construct a finite dimensional space $V_{H}$ ($V_{N}$), and to find an approximate solution $z_{H}\in V_{H}$, such that,

• •

$z_{H}$ has guaranteed error estimate in certain norm $\|\cdot\|$, e.g., $\|z-z_{H}\|\leq CH^{\alpha}$ or equivalently $\|z-z_{N}\|\leq CN^{-\alpha/d}$, with $C$ independent of $\varepsilon$ (and contrast) and $\alpha$ is the optimal convergence rate.

• •

$V_{H}$ is constructed via precomputed subproblems which are optimally localized and can be solved in parallel, also those subproblems do not depend on the forcing term and boundary condition (analog of cell problems in classical homogenization).

The basis of numerical homogenization needs to be pre-computed. The computation of each basis will be independent and the support of each basis needs to be localized on a small patch. The possibility to compute such bases on localized sub-domains of the global domain without loss of accuracy is therefore a problem of practical importance. We refer to [13, 16, 2, 38, 34] for recent localization results for divergence-form elliptic PDEs.

We will introduce the generalized rough polyharmonic splines (GRPS) in [39, 33]. Given $N$ measurement functions $\psi_{i}$, $i=1,\dots,N$, define

[TABLE]

The GRPS basis is given by the solution of the following constrained minimization problem which is strictly convex and has a unique minimizer $\phi_{i}\in V_{i}$,

[TABLE]

for an appropriate norm $\|\cdot\|$.

Remark 3.1**.**

The GRPS basis can be given by the Bayesian Formulation in [35]. We can ask the following question: Given $N$ observables $\Psi_{i}=\int_{\Omega}z(x)\psi_{i}(x)\textrm{dx}$ of the solution $z(x)$ of $-\operatorname{div}(a\nabla z)=\rho$, $i=1,\dots,N$, what is the best guess of $z(x)$? The answer can be given by the following procedure of randomization and conditioning.

Randomization: Put a prior on $z(x)$, e.g. $z(x)$ is given by

[TABLE]

where $\xi(x)$ is a Gaussian field with covariance function $\Lambda(x,y)$, therefore $z(x)$ is a Gaussian field with covariance function

[TABLE] 2. 2.

Conditioning: take the conditional exception $\mathbb{E}[z(x)|\Psi]$, we have

[TABLE]

and $\Theta_{i,j}:=\int_{\Omega^{2}}\psi_{i}(x)\Gamma(x,y)\psi_{j}(y)\textrm{dx}\textrm{dy}$ is the covariance matrix of $\Psi$. The $\phi_{i}$ given by the conditional expectation in (16) is exactly the one given in the variational formulation in (15) (by choosing an appropriate $\|\cdot\|$).

*this Bayesian framework can be further generalized to the game/decision theoretical framework as in [36]. *

In fact, we have some flexibility to choose $\|\cdot\|$ in (15) (corresponding to $\Lambda(x,y)$ in the Bayesian framework). For example, it can be taken as $\|\operatorname{div}a(x)\nabla u\|_{L^{2}(\Omega)}$ which is used in the rough polyharmonic splines (RPS) paper [39], or the energy norm $\|u\|_{a}$ as in [36, 33].

We also have different choices for the measurement function $\psi_{i}$, for example,

• •

$\psi_{i}=\delta(x-x_{i})$ (point value observables), together with the norm $\|\operatorname{div}a(x)\nabla z\|$, we recover the RPS basis in [39];

• •

$\psi_{i}$ as characteristic functions of patches in a coarse triangulation $\mathcal{T}_{H}$, observables are patch averages of $z$. We refer to this basis as the GRPS basis in the current paper;

• •

$\psi_{i}$ as characteristic functions of edges in a coarse triangulation $\mathcal{T}_{H}$, observables are edge averages of $z$ ;

• •

For a quasi-interpolation operator $\mathfrak{I}_{H}$ (Clement, Oswald, etc.), there exists $\{\psi_{i}\}$ such that $\mathfrak{I}_{H}(z)(x_{i})=\int_{\Omega}\psi_{i}z\textrm{dx}$, which recovers the LOD approach by Peterseim et.al [34].

The GRPS approach naturally induces a two-level decomposition: Define the coarse space as $V_{H}:={\rm span}\{\phi_{i}\}$, the fine space can be naturally defined as (given the full space $V=H^{1}_{0}(\Omega)$ or $V_{h}$),

$V_{f}:=\{v\in X|\int\psi_{i}v\textrm{dx}=0,\forall i\}.$

The $\langle\cdot,\cdot\rangle$ inner product is induced by the norm $\|\cdot\|$. Then $V_{H}\perp V_{f}$ with respect to the $\langle\cdot,\cdot\rangle$ product. Furthermore, we have the following optimal recovery property of $z(x)$ in $V_{H}$. Let $z_{\rm I}:=\sum_{i}(\int_{\Omega}\psi_{i}z\textrm{dx})\phi_{i}(x)$ be the interpolant of u in $V_{H}$, we have

[TABLE]

The following properties of GRPS is important for the proof of the convergence for the optimal control problem, the proof of those results can be found in [33, 39, 36].

Theorem 1**.**

[Optimal Approximation Property of $V_{H}$]:

[TABLE]

This is true for the following constructions

• •

RPS basis, using higher order Poincare inequality $\|\nabla v\|_{L^{2}(\Omega)}\leq CH\|{\rm div}a\nabla v\|_{L^{2}(\Omega)}$ for $v\in V_{f}$.

• •

GRPS basis, using Poincare inequality for $v\in V_{f}$.

• •

LOD construction, approximation property for quasi-interpolation operator $\mathfrak{I}_{H}$: $H^{-1}\|v-\mathfrak{I}_{H}v\|_{L^{2}(T)}+\|\nabla(v-\mathfrak{I}_{H}v)\|_{L^{2}(T)}\leq C\|\nabla v\|_{L^{2}(\omega_{T})}.$

Note that in the variational formulation (15), the minimization is done for functions defined on whole $\Omega$, and we call the corresponding basis global basis.

Theorem 2**.**

For the finite element solution $z_{H}\in V_{H}$, where $V_{H}$ is the space of global basis. We have

[TABLE]

Of course, it is not preferable to use global basis in practical computation. The good thing is, the global basis has the following exponential decay property which can be proved by using Cacciopoli like argument for harmonic functions [34, 36].

Theorem 3**.**

We have the following exponential decay property of global basis.

[TABLE]

The exponential decay property opens an avenue for the local approximation of global basis. Let $\Omega_{i}\subset\Omega$, introduce

[TABLE]

The local basis is given by the solution of the following constrained minimization problem which is strictly convex and has a unique minimizer $\phi^{\textrm{loc}}_{i}\in V^{\textrm{loc}}_{i}$,

[TABLE]

We have the following properties of the localized basis.

Theorem 4**.**

[Truncation error for the localized basis]: If ${\rm diam}(\Omega_{i}:={\rm supp}(\phi^{\rm loc}_{i}))=r$, then

[TABLE]

and the convergence of the FEM with localized basis

Theorem 5**.**

[Accuracy of FEM with Localized Basis]: If $\displaystyle(B(x_{i},CH\ln\frac{1}{H})\cap\Omega)\subset\Omega_{i}$, $z_{H}^{\textrm{loc}}$ is the FEM solution in ${\rm span}\{\phi_{i}^{\textrm{loc}}\}$, the space of localized basis, then

[TABLE]

4 Convergence Analysis for Optimal Control Problem

In the rest of the paper, we use the localized numerical homogenization basis $V_{H}:=\rm span\{\phi_{i}^{loc}\}_{i\in\mathcal{N}}$ which satisfies the Theorem 5.

Consider the following optimal control problem

[TABLE]

where $U_{H}$ is the piecewise constant finite element space over $\mathcal{T}^{H}_{U}$, and $K_{H}:=K\cap U_{H}$.

4.1 A priori error estimates

Fixed the control approximation $u_{H}\in V_{H}$, define the auxiliary solutions $(y(u_{H}),p(u_{H}))\in V\times V$ which are the solutions of the following equations:

[TABLE]

We have the following lemma for the accuracy of auxiliary solutions.

Lemma 1**.**

Let $y(u_{H})$ and $p(u_{H})$ be the solutions of (23), $y$ and $p$ be the finite element solutions of (8) in $V$. It holds true that

[TABLE]

*where $C$ depends on $a_{\min}$, $a_{\max}$, $d$, $C_{B}$, $L_{g}$.

Proof 4.1**.**

By (8) and (23), we have

[TABLE]

which implies that

[TABLE]

The following lemma bounds the accuracy of the approximate solution of (22).

Lemma 2**.**

Let $y(u_{H})$ and $p(u_{H})$ be the solutions of (23), $y_{H}$ and $p_{H}$ be the finite element solutions of (22) in $V_{H}$. It holds true that

[TABLE]

where $\displaystyle C=[\frac{a_{\max}}{a_{\min}}(1+d/\pi)]^{1/2}C_{\mathrm{I}}$. And,

[TABLE]

*where $\displaystyle C=\max\{(\frac{a_{\max}}{a_{\min}}(1+d/\pi)\sqrt{2}+\frac{1}{\sqrt{2}})C_{NH},\frac{L_{g}}{a_{\min}}(1+d/\pi)\sqrt{2}\}$, $d=\mathrm{diam}(\Omega)$.

Proof 4.2**.**

The first inequality (24) is due to usual finite element estimate, Poincaré inequality, and (17),

[TABLE]

where $y(u_{H})_{\mathrm{I}}$ is the interpolation of $y(u_{H})$ in $V_{H}$.

By Poincaré inequality, Lipschitz property (3) of $g^{\prime}(\cdot)$, (22), (23), we obtain

[TABLE]

where $(p(u_{H}))_{\mathrm{I}}$ is the interpolation of $p(u_{H})$ in $V_{H}$.

Take $\displaystyle\varepsilon=\frac{1}{2}\frac{a_{\min}}{1+d/\pi}$, we have

[TABLE]

where the constants $C$ depends on $a_{\min}$, $a_{\max}$, $d$, $L_{g}$, $C_{I}$ and $C_{\mathrm{NH}}$, but not on $H$.

Define the averaging projection([32]) $\Pi_{H}:K\rightarrow K_{H}$ as,

[TABLE]

Note that $K_{H}\subset K$. We have that if $v\in H^{1}(\Omega_{U})$ (see e.g.[32, 15]),

[TABLE]

Moreover, if $u\in H^{1}(\Omega_{U})$, $p\in H^{1}(\Omega)$, assume that $h^{\prime}(u)$ is Lipschitz continuous, we have

[TABLE]

Lemma 3**.**

Let $(y,p,u)$ be the solution of equation (8), and $(y_{H},p_{H},u_{H})$ be the finite element solution of (22). Assume that $u\in H^{1}(\Omega_{U})$, it holds true that

[TABLE]

*where $C$ is a constant depending on $a_{\min}$, $a_{\max}$, $L_{h}$, $L_{g}$, $C_{B}$, $C_{I}$ and $C_{\mathrm{NH}}$, but not on $H$ and $H_{U}$. *

Proof 4.3**.**

It follows from (8) that

[TABLE]

Similarly,

[TABLE]

The convexities of $g$ and $h$ imply that,

[TABLE]

Combining (29), the optimality conditions (8) and (22), the estimate (28), we have

[TABLE]

where $\varepsilon$ is a arbitrarily positive constant.

By Lemma 1, we have

[TABLE]

Thus, if we choose $\displaystyle\varepsilon=\frac{M_{h}}{2+C^{2}}$, it follows from (27), (31) and (30) that,

[TABLE]

Using (24), (25), Lemma 1, and choosing $H$ sufficient small such that $CL_{g}^{2}H^{2}\leq 1$, we have

[TABLE]

Thus, it follows from (32) and (33) that,

[TABLE]

*By choosing $H$ sufficient small such that $\displaystyle 1-C(C_{B}^{2}+L_{g}^{2})H^{2}\geq\frac{1}{2}$, we conclude the proof of the theorem. *

Combining Lemmas 1, 2, 3, we have the following a priori error estimates.

Theorem 4**.**

Let $(y,p,u)$ be the solution of equation (8), and $(y_{H},p_{H},u_{H})$ be the finite element solution of (22). Assume that $u\in H^{1}(\Omega_{U})$, then it holds true that

[TABLE]

Proof 4.4**.**

Note that

[TABLE]

Lemma 1 leads to

[TABLE]

Lemma 2 leads to

[TABLE]

Lemma 3 leads to

[TABLE]

*Combining all those estimates, we conclude the proof by taking $\displaystyle CL_{g}H\leq\frac{1}{2}$. *

5 Algorithm and Numerical Experiments

5.1 Numerical algorithm

To solve the optimal control problem (22), we will introduce the following projection algorithm.

Define the projection operator $P^{b}_{K}:U\rightarrow K$: for $w\in U$, find $P^{b}_{K}w\in K$ such that

[TABLE]

which is equivalent to the inequality

[TABLE]

It is clear that $P^{b}_{K}$ is well-defined for any closed convex subset $K\subset U$. For example, when $K=\{v\in L^{2}(\Omega_{U})|\int_{\Omega_{U}}v\geq 0\}$, for any $w\in U$, we have

[TABLE]

where $\bar{w}=\int_{\Omega_{U}}w/\int_{\Omega_{U}}1$ denotes the average of $w$ on $\Omega_{U}$. The formulas for other important cases can be found in [32].

We have the following lemma.

Lemma 1**.**

For the solution $u$ of

[TABLE]

or the equivalent optimality condition

[TABLE]

where $J(u)$ is a convex functional of $u$. It holds that

[TABLE]

Furthermore, for any $u$, $v\in U$

[TABLE]

Proof 5.1**.**

[TABLE]

therefore $u=P^{b}_{K}(u-\rho J^{\prime}(u))$.

Furthermore,

[TABLE]

Hence,

[TABLE]

The convergence results for this algorithm are given in Theorem 8.2.1 and Remark 8.2.1 of [32].

Let $V_{h}\subset H^{1}_{0}(\Omega)$ be the fine scale finite element space associated with fine mesh triangulation $\mathcal{T}^{h}$, and $U_{h}\subset L^{2}(\Omega)$ is the piecewise constant finite element space associated with $\mathcal{T}^{h}_{U}$. $V_{H}$ and $U_{H}$ are the corresponding coarse mesh finite element spaces introduced in § 3. We have the following iterative algorithm to solve (22).

Now, we are in the position to prove the convergence of the above algorithm.

Theorem 2**.**

The triplet $(y^{(n)},u^{(n)},p^{(n)})$ in Algorithm 1 converges to the triplet $(y_{H},u_{H},p_{H})$ in (22). To be more precise, if we take $\rho^{(n)}=\rho$ with $0<\rho<1$ such that $1-2\rho M_{h}+2\rho^{2}L_{h}^{2}<\delta^{2}/2$ and $2\rho C^{2}_{B}L_{g}+2\rho^{2}C_{B}^{4}L_{g}^{2}<\delta^{2}/2$, with $0<\delta<1$, we have

[TABLE]

Proof 5.2**.**

By (37) and (22), we have the following equation for $y^{(n)}-y_{H}$ and $p^{(n)}-p_{H}$,

[TABLE]

Therefore, we have

[TABLE]

where the constant $C$ only depends on $a_{\min}$, $a_{\max}$, and $\mathrm{diam}(\Omega)$.

For simplicity of notation, we refer to $P^{b}_{K}$ as the projection from $U_{H}$ to $K_{H}$, and let $\rho^{(n)}=\rho$. Therefore, $u^{(n+1)}=P^{b}_{K}(u^{(n)}-\rho(h^{\prime}(u^{(n)})+B^{*}p^{(n)}))$, and $u_{H}=P^{b}_{K}(u_{H}-\rho(h^{\prime}(u_{H})+B^{*}p_{H}))$. Hence,

[TABLE]

Take $0<\rho<1$ such that $1-2\rho M_{h}+2\rho^{2}L_{h}^{2}<\delta^{2}/2$ and $2\rho C^{2}_{B}L_{g}+2\rho^{2}C_{B}^{4}L_{g}^{2}<\delta^{2}/2$, with $\delta<1$, we have

[TABLE]

*Therefore, $\|y^{(n)}-y_{H}\|_{1}+\|p^{(n)}-p_{H}\|_{1}+\|u^{(n)}-u_{H}\|\leq C\delta^{n}\|u_{H}\|$. *

We can also write down the matrix form of the Algorithm 1 in the following. Denote

[TABLE]

where $i,j=1,2,\ldots,N$; $k,l=1,2,\ldots,m$. The iterative scheme for gradient descent algorithm is as follows:

Therefore, the main computational cost of the iterative algorithm is to solve the state and control equations.

5.2 Numerical Results

In this section, we present two numerical examples to verify the error estimates presented in the previous sections. We consider the elliptic optimal control problem with rough coefficients,

[TABLE]

and the constraint set $K=\{v\in L^{2}(\Omega)|\int_{\Omega}v\geq 0\}$. The dual equation of the state equation is given by,

[TABLE]

For simplicity, we use Dirichlet boundary condition and set $f=1$. The continuous linear operator is chosen $B=I$, where $I$ is the identity operator.

We test the numerical methods for two different types of diffusion coefficients. The first is a multiscale trigonometric function, and the second is the SPE10 benchmark for reservoir simulation (http://www.spe.org/web/csp/).

5.2.1 Multiscale Trigonometric Example

For the first example, $a(x)$ is a scalar function given by the following expression,

[TABLE]

where $\displaystyle\varepsilon_{1}=\frac{1}{5},\varepsilon_{2}=\frac{1}{13},\varepsilon_{3}=\frac{1}{17},\varepsilon_{4}=\frac{1}{31},\varepsilon_{5}=\frac{1}{65}.$ The coefficient $a(x)$ is highly oscillatory with non-separable scales, as shown in Figure 1.

We take the domain as the unit square $\Omega_{U}=\Omega=[0,1]\times[0,1]$. The regular coarse mesh $\mathcal{T}_{H}$ is obtained by first subdividing $\Omega$ uniformly into $N_{c}\times N_{c}$ squares, then each square can be partitioned into two triangles along the $(1,1)$ direction. We can further refine the coarse mesh uniformly by dividing each triangle into four similar subtriangles. We refine $\mathcal{T}_{H}$ $J$ times to obtain the fine mesh $\mathcal{T}_{h}$, therefore, $H=2^{J}h$. Let us refer to Figure 2 for an illustration of the partition. The degrees of freedom of global RPS basis and the global GRPS basis are $(N_{c}-1)^{2}$ and $2N_{c}^{2}$, respectively.

We compute localized RPS basis $\phi_{i}^{l}$ on localized sub-domains $\Omega_{i}^{l}$ defined by adding $l$ layers of coarse triangles around coarse node $x_{i}$. More precisely $\Omega_{i}^{1}$ is the union of triangles $T\in\mathcal{T}_{H}$ with $x_{i}$ as a common node, and $\displaystyle\Omega_{i}^{l+1}:=\cup_{T\cap\Omega_{i}^{l}\neq\emptyset}T$.

Similarly, we compute localized GRPS basis $\phi_{i}^{l}$ on localized sub-domains $\bar{\Omega}_{i}^{l}$ defined by adding $l$ layers of coarse triangles around coarse triangle $T_{i}$. More precisely $\bar{\Omega}_{i}^{1}$ is the union of triangles $T\in\mathcal{T}_{H}$ with $T_{i}$ as common node, common edge, and $\displaystyle\bar{\Omega}_{i}^{l+1}:=\cup_{T\cap\bar{\Omega}_{i}^{l}\neq\emptyset}T$. We refer to Figure 3 for an illustration of the patches $\Omega_{i}^{l}$ and $\bar{\Omega}_{i}^{l}$ of the RPS basis and GRPS basis with $l=1,2,3$.

To better understand the exponential decay properties of the RPS basis and GRPS basis, we plot the global RPS basis and the localized RPS basis in Figure 4. The coarse mesh has $H=1/32$, and the fine mesh has $h=1/256$ (namely, $N_{c}=32$, $J=3$). Figure 4(a) shows the shape of global basis $\phi_{i}$ centered at $(0.5,0.5)$ in the $\log_{10}$ scale. Figure 4(b) shows a slice of $\phi_{481}$ along the x-axis. Figure 5 shows the localized basis $\phi_{i}^{loc}$ for the same node for various levels of localization(i.e., for $l=1,4,7$). It is consistent with the exponential decay and localization results.

We compute the optimal control problem using Algorithm 1 with RPS space $V_{H}$ or GRPS space $\bar{V}_{H}$.

Figure 6 shows the relative error of $y_{H}$ and $p_{H}$ using GRPS $\bar{V}_{H}$, in $L^{2}$, $H^{1}$ and $L^{\infty}$ norm in the $\log_{10}$ scale as a function of the number of layers $l$.

Figure 7(a), (b) and (c) show the relative errors of $\|y-y_{H}\|_{H^{1}(\Omega)}$, $\|p-p_{H}\|_{H^{1}(\Omega)}$ and $\|u-u_{H}\|_{L^{2}(\Omega)}$ in $\log_{10}$ scale with respect to the number of layers $l=1,\cdots,9$, and (d), (e), (f) show the relative error with respect to coarse dof. Note that the fine mesh is fixed with $h=1/256$.

Figure 8(a)(b) shows that: with fixed $H$, when the number of layers increase , $\|y-y_{H}\|_{1}+\|p-p_{H}\|_{1}+\|u-u_{H}\|$ decreases first and then saturates; as $H$ decrease, and $l\sim 1/H\log(1/H)$, we have the optimal convergence rate $O(H)$. This is consistent with the Theorem 4.

For comparison with GRPS solution, we show in Figure 9 the relative errors of RPS solutions, for each variable $y_{H}$, $p_{H}$ and $u_{H}$, and show in Figure 10 the combined error of the RPS solutions of the optimal control problem. It seems the numerical performance of RPS and GRPS are similar, and GRPS is a little more stable.

5.2.2 SPE 10

The purpose of this section is to show the performance of GRPS method for a more practical example, namely, the SPE10 benchmark problem (http://www.spe.org/web/csp/). SPE10 is the latest industry benchmark problems (SPE10) from the Society of Petroleum Engineers (SPE). The physical domain is a cube with dimension $[0,220]\times[0,60]\times[0,85]$. The coefficients $a(x,y,z)$ are given as piecewise constant numerical values rather than a continuous function. We select the layer 39 with respect to the height dimension and treat them as coefficients of two dimensional problems. The contrast of the coefficients is $1.74956\times 10^{7}$. The domain is chosen as $\Omega=[0,2.2]\times[0,0.6]$. We first uniformly divide $\Omega$ with the coarse mesh size $H$, then we further refine the mesh with the fine mesh size $h$. In the numerical experiment, we choose $H=1/25,1/50,1/100,1/200$, respectively. The fine mesh size is fixed as $h=1/800$. We illustrate the SPE10 coefficients in Figure 11.

The errors of GRPS solutions are shown in Figures 12 and 13.

6 Conclusions

In this paper, we have introduced the generalized rough polyharmonic splines (GRPS, including RPS) method for the efficient solution of optimal control problem govern by multiscale elliptic equation with rough coefficients. We have derived rigorous error estimates and the numerical experiments complement well with the theoretical analysis.

We plan to apply this strategy to optimal control problem with more general control conditions, such as boundary control, as well as optimal control problem for elasticity and Stokes equations with rough coefficients.

The numerical homogenization based method in this paper can provide coarse scale accuracy of the optimal control solution. If fine scale accuracy is desired, numerical homogenization may help design efficient preconditioners, and furthermore, efficient mutlgrid/multilevel method such as Gamblet based multigrid method [36, 46], can be used for the efficient resolution of the optimal control governed by multiscale problems.