Adaptive finite element approximation of sparse optimal control with integral fractional Laplacian

In this paper, we present and analyze a weighted residual a posteriori error estimate for an optimal control problem involving a nondifferentiable cost functional, a state equation with an integral fractional Laplacian, and control constraints. For a bounded Lipschitz domain $\Omega\subset R^{d},\Omega^{c}:=R^{d}\backslash\overline{\Omega}$, we consider

subject to

and the control constraints

Here parameters $\gamma>0$ and $\beta>0$. In the sequel, $y$ is state and $u$ is the control variable. The function $y_{d}\in L^{2}(\Omega)$ is referred to as desired state. To focus on the scenario of nondifferentiability, it is assumed that $a,b\in R$ satisfy the condition that $a<0<b$. We notice that the set $U_{ad}$, is a nonempty, bounded, closed and convex subset of $L^{2}(\Omega)$.

The introduction of nonsmooth regularization term $\beta\|u\|_{L^{1}(\Omega)}$ in the PDE-constrained optimization problems promotes sparsity in the solutions. This allows the control variable in the optimization process to tend towards zero in regions where it has negligible impact on the cost function, therefore minimizing the cost function. Sparse optimization is widely used in many practical applications, especially in the processing and analysis of high dimensional data, such as noise processing, machine learning, face recognition, etc.

Previous research has addressed the analysis of optimal control problems with a cost term containing $L^{1}(\Omega)$ norm [1, 2, 3, 4, 5, 6, 7, 8, 20]. For instance, in [7], the authors investigated the $L^{1}(\Omega)$ control problem constrained by a linear elliptic PDE, where the objective functional incorporated a regularization technique based on the $L^{2}(\Omega)$ control cost term. The authors analyzed the optimality conditions and proposed a semismooth Newton method that achieves local convergence with superlinear speed. Building upon this work, [8] provided a priori and posteriori error estimates through finite element analysis. Furthermore, in [5], the authors considered a semilinear elliptic PDE as the state equation and analyzed second-order optimality conditions. Additionally, the authors in [6] studied the sparse control problem with a fractional diffusion equation as the state equation. They analyzed a priori error estimate for the fully discrete case using finite element methods. More recently, Ot$\rm{\acute{a}}$rola et al. studied a sparse optimal control problem with a non-differentiable cost functional, where the state equations are Poisson’s problem and fractional diffusion equation, respectively in [9, 10]. In [9], the authors studied three different strategies for approximating the control variable, they proposed and analyzed a reliable and efficient a posteriori error estimate, and designed an adaptive strategy to achieve optimal convergence rates. In [10], the authors investigated an adaptive finite element method for sparse optimal control of fractional diffusion, taking into account the spectral definition of the fractional Laplacian operator.

In comparison to the priori error analysis of finite element approximations for PDE-constrained optimization, the design and analysis of a posteriori error estimate is not much. The initial work on reliable a posteriori error estimation for optimal control problems was presented in [11], followed by a series of related studies [12, 13, 14, 15, 16]. Residual-based a posteriori error estimates incorporating data oscillation were introduced in [17]. Later, a unified framework for the a posterior error analysis of linear quadratic optimal control problems with control constraints was established in [18], and the pure convergence of an adaptive finite element method for optimal control problems with variable divergence control was proved, that is, convergence without convergence rate. In [19], the authors rigorously proved convergence and quasi-optimality of AFEM for optimal control problem involving state and adjoint state variable.

However, to the best of our knowledge, no previous research has combined adaptive finite element methods (AFEMs) with integral fractional Laplacian sparse optimal control to address such problems. Therefore, in this paper, we focus on the adaptive finite element approximation for sparse optimal control with the integral fractional Laplacian. We outline and analyze the solution methodology for problem (1.1)-(1.2) based on the following considerations:

$\bullet$ The optimal control problem involving the fractional Laplacian operator can effectively simulate groundwater pollution [21], turbulent flow [22], and chaotic dynamics [23]. Unlike integer-order diffusion equations, the fractional Laplacian operator exhibits power-law decay, which accurately captures heavy-tailed power-law decay phenomena observed in these applications. Hence, studying fractional optimal control problem is essential.

$\bullet$ Objective functional by introducing the $L^{1}(\Omega)$ norm to control some specific physical quantities or locations, and the $L^{2}(\Omega)$ norm to maintain smoothness and continuity, can better solve the practical problems that need to control the optimal cost. The existence of $L^{1}(\Omega)$ term in the objective function requires us to derive the first-order condition using a subdifferential approach [5, 7, 8], which is different from distributed optimal control problems.

$\bullet$ Due to the non-locality, non-differentiability, and intrinsic constraints of the fractional Laplacian operator, by adopting adaptive strategies and a posteriori error estimate, we can identify singularities and refine the mesh accordingly, which can more effectively allocate computational resources and achieve higher accuracy with lower computational costs. One of the challenges in designing the a posterior error estimator is the nature of the residual, that is, it is not necessarily in $L^{2}(\Omega)$. We refer to [26] and introduce the weighted residual estimator, where the weights are given by the power of the distance to the grid skeleton

(see section 4).

$\bullet$ The adaptive finite element method is widely used, but there are not many convergence analyses of the algorithm. The optimal control problem we studied is a coupled system with nonlinear properties, which leads to the lack of orthogonality presented in [30] and brings difficulties to our convergence analysis. In order to address this issue, we refer to reference [28] and prove its quasi-orthogonality.

Recently, the only work on a posteriori error analysis for sparse optimal control constrained by fractional order equations, as in (1.1)-(1.2), is found in [10]. Compared with [10], this paper studies the integral definition of the fractional Laplacian operator, which plays an important role in the modeling of complex non-local and nonlinear phenomena such as diffusion, heat transfer, resistance and elasticity. And the main difference is that the convergence of the adaptive algorithm is also analyzed in this paper. In this paper, we use piecewise linear polynomial dispersion for the state variable and variational discretization for the control variable. We design a posterior error estimator that requires only discretization of the state variable and adjoint variable. Notably, in the $\frac{3}{2}<\alpha<2$, the residual does not satisfy the $L^{2}(\Omega)$-regularity. To address this issue, an additional weight based on the power of the distance from the mesh skeleton is included in the weighted residual estimator. An h-adaptive algorithm driven by the D$\rm{\ddot{o}}$rfler marking criterion based on the a posteriori error estimator is proposed and its convergence is proved.

The organization of the paper is as follows: In section 2, we introduce the symbols used and provide a brief overview of elements in convex analysis, along with the regularity of solutions to optimal control problem. In section 3, we analyze the first-order optimality conditions for the problem. In section 4, we introduce the finite element discretization of the optimal control problem (1.1)-(1.2) and design a weighted residual estimator. The core of our work is presented in sections 4 and 5. For the discretization introduced at the beginning of section 4, we first derive upper and lower bounds for the a posteriori error estimate of the finite element approximation for the optimal control problem. An h-adaptive algorithm driven by the posterior view error estimator based on D$\rm{\ddot{o}}$rfler labeling criterion is proposed. In section 5, we show that the sequence of approximations produced by the adaptive algorithm converges at the optimal algebraic rate. In section 6, a series of numerical examples are provided to demonstrate the effectiveness of our theoretical findings.

In this section we introduce some preliminaries about fractional Sobolev spaces, subdifferential and fractional Laplacian. For a bounded domain $\Lambda\subset R^{d},L^{2}(\Lambda)$ denotes the Banach spaces of standard 2-th Lebesgue integrable functions on $\Lambda$. For $s\in(0,1),$ $H^{s}(\Lambda)$ denotes the fractional Sobolev space. $H_{0}^{s}(\Lambda)$ is the subspace of $H^{s}(\Lambda)$ consisting of functions whose trace is zero on $\partial\Lambda$. Let $(\cdot,\cdot)$ and $\|\cdot\|$ denote the inner product and norm in $L^{2}(\Lambda)$, respectively. The seminorm $|\cdot|_{H^{s}(\Lambda)}$ and the full norm $\|\cdot\|_{H^{s}(\Lambda)}$ are denoted as follows

Moreover, we introduce the following space, which will be used in the weak formulation of state equation

Next, we will review some concepts with respect to subdifferentials from convex analysis that will be useful in our upcoming analysis. For details, please refer to [24]. Consider a real and normed vector space $G$. Suppose $\phi:G\rightarrow R\cup\{\infty\}$ be a convex and proper functional. Let $v\in G$ be such that $\phi(v)<\infty$. A subgradient of $G$ at $v$ is an element $v^{*}\in G^{*}$ that satisfies

Here, $\langle\cdot,\cdot\rangle_{G^{*},G}$ represents the duality pairing between $G^{*}$ and $G$. The set of all subgradients of $\phi$ at $\bar{v}$, denoted by $\partial\phi(\bar{v})$, refers to the subdifferential of $\phi$ at $\bar{v}$.

As $\phi$ is a convex functional, the subdifferential at any point $v$ within the effective domain of $\phi$ is not empty. Additionally, it is important to note that the subdifferential is monotone, i.e.,

Finally, we introduce the definition of fractional Laplacian:

Here $0<\alpha<2$, and

and ”p.v.” denotes the principal value of the integral:

where $B_{\epsilon}(v)$ is a ball of radius $\epsilon$ centered at $x$. The difference $y(x)-y(w)$ in the numerator of (2.3), which vanishes at the singularity, provides a regularization, which together with averaging of positive and negative parts allows the principal value to exist. A consequence of this definition is the mapping property (see [25]).

The weak formulation of state equation (1.2) reads: Find $y\in\widetilde{H}^{\frac{\alpha}{2}}(\Omega)$ such that

Here

We define

As the $H^{\frac{\alpha}{2}}(R^{d})$ seminorm is equivalent to the $H^{\frac{\alpha}{2}}(R^{d})$ norm on $\widetilde{H}^{\frac{\alpha}{2}}(\Omega)$ (see [29]), by Lax-Milgram theorem, the solution $y\in\widetilde{H}^{\frac{\alpha}{2}}(\Omega)$ exists and is unique.

For the state equation with the right hand $f$ we can define a linear and bounded solution operator $\mathcal{S}:L^{2}(\Omega)\longrightarrow\widetilde{H}^{\frac{\alpha}{2}}(\Omega)$ such that $y=\mathcal{S}f$. Moreover, the following regularity result holds for the state equation.

The weak formulation of the optimal control problem (1.1)-(1.2) reads:

subject to

Since the $J$ is strictly convex and weakly lower semicontinuous, this problem admits a unique optimal solution $(y,u)\in\widetilde{H}^{\frac{\alpha}{2}}(\Omega)\times L^{2}(\Omega)$.

In order to obtain optimality conditions for (3.2)-(3.3), we introduce the following adjoint state $p$ as follows:

Set

and $j_{2}(u)=\|u\|_{L^{1}(\Omega)}$. Then we obtain the reduced problem of (1.1):

Although the reduced cost functional (3.5) is nonsmooth, it consists in the sum of a regular part and a convex nondifferentiable term. Thanks to the structure, optimality conditions can still be established according to the following result.

For $a,b\in R$, we introduce a projection operator $\Pi_{[a,b]}:L^{2}(\Omega)\rightarrow U_{ad}$ defined by

Then we have the following projection formulas.

At end we present the following first order optimality conditions for above optimal control problems.

We begin by partitioning the domain $\Omega$ into simplices $K$ with size $h_{K}:=|K|^{\frac{1}{d}}$, forming a conforming partition $\mathcal{T}_{h}=\{K\}$. We then define $h_{\mathcal{T}_{h}}=\max\limits_{K\in\mathcal{T}_{h}}h_{K}$ and denote by $\mathbb{T}$ the collection of conforming and shape regular meshes that are refinements of an initial mesh $\mathcal{T}_{h_{0}}$. For $\mathcal{T}_{h}\in\mathbb{T},$ let $\mathbb{V}_{\mathcal{T}_{h}}$ be the finite element space consisting of continuous piecewise linear functions over the triangulation $\mathcal{T}_{h}$

For all elements $K\in\mathcal{T}_{h}$ and $k\in\mathbb{N}_{0}$, we introduce the $k$-th order element patch inductively by

The finite element approximation of the optimal control problem (3.2)-(3.3) can be characterized as

subject to

Here the admissible set of control $U_{ad}$ is not discretized, i.e., the so-called variational discretization approach. Similar to the continuous case we have the discrete first order optimality condition

where $\lambda_{\mathcal{T}_{h}}\in\partial j_{2}(u_{\mathcal{T}_{h}}).$ Next, we give the following discrete projection formula.

Similar to the continuous case, we can define a discrete control-to-state mapping $\mathcal{S}_{h}:L^{2}(\Omega)\longrightarrow\mathbb{V}_{\mathcal{T}_{h}}$. Set $\Theta(h)=\sup\limits_{f\in L^{2}(\Omega),\|f\|=1}\inf\limits_{\chi_{\mathcal{T}_{h}}\in\mathbb{V}_{\mathcal{T}_{h}}}\|\mathcal{S}f-\chi_{\mathcal{T}_{h}}\|_{\widetilde{H}^{\frac{\alpha}{2}}(\Omega)}.$ Then we have $\Theta(h)\ll 1$ for $h\in(0,h_{0})$ with $h_{0}\ll 1.$

To derive a posteriori error analysis we need to introduce the following auxiliary problems

Note that the residuals do not satisfy the $L^{2}(\Omega)$-regularity for $\frac{3}{2}<\alpha<2$. To address this issue, we require a weight function to measure the distance from the mesh skeleton. For a mesh $\mathcal{T}_{h}$, we introduce the weight function defined in [26]

We further define the corresponding weighted residual errors as follows

where

Then on a subset $\omega\subset\Omega$, we define the error estimators of the state and adjoint state by

Thus, $E_{y}(y_{\mathcal{T}_{h}},\mathcal{T}_{h})$ and $E_{p}(p_{\mathcal{T}_{h}},\mathcal{T}_{h})$ constitute the error estimators for the state equation and the adjoint state equation on $\Omega$ with respect to $\mathcal{T}_{h}$ as follows

Moreover we also need the Scott-Zhang operator([26]) $\Pi_{\mathcal{T}_{h}}:L^{2}(\Omega)\rightarrow\mathbb{V}_{\mathcal{T}_{h}}$ that satisfy the following properties