A Multilevel Heterogeneous ADMM Algorithm for Elliptic Optimal Control Problems with L¹-Control Cost

Abstract

In this paper, elliptic optimal control problems with $L^1$-control cost and box constraints on the control are considered. To numerically solve the optimal control problems, we use the First optimize, then discretize approach. We focus on the inexact alternating direction method of multipliers (iADMM) and employ the standard piecewise linear finite element approach to discretize the subproblems in each iteration. However, in general, solving the subproblems is expensive, especially when the discretization is at a fine level. Motivated by the efficiency of the multigrid method for solving large-scale problems, we combine the multigrid strategy with the iADMM algorithm. Instead of fixing the mesh size before the computation process, we propose the strategy of gradually refining the grid. Moreover, to overcome the difficulty whereby the $L^1$-norm does not have a decoupled form, we apply nodal quadrature formulas to approximately discretize the $L^1$-norm and $L^2$-norm. Based on these strategies, an efficient multilevel heterogeneous ADMM (mhADMM) algorithm is proposed. The total error of the mhADMM consists of two parts: the discretization error resulting from the finite-element discretization and the iteration error resulting from solving the discretized subproblems. Both errors can be regarded as the error of inexactly solving infinite-dimensional subproblems. Thus, the mhADMM can be regarded as the iADMM in function space. Furthermore, theoretical results on the global convergence, as well as the iteration complexity results $o(1/k)$ for the mhADMM, are given. Numerical results show the efficiency of the mhADMM algorithm.

1. Introduction

Sparse optimal control problems are widespread in many areas, such as the placement of actuators [1], quantum spin systems [2], etc. [3,4]. It is well-known that adding the $L^1$-norm penalty can lead to sparse optimal control problems, i.e., the infinite dimensional control variable is localized in its domain of action. In this paper, we consider the elliptic optimal control problem with $L^1$-control cost and box constraints on the control:

$$\left\{\begin{array}{cc}\hfill \underset{(y,u)\in Y\times U}{\overset{}{min}}& J(y,u)=\frac{1}{2}\parallel y-y_d{\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\frac{\alpha }{2}{\parallel u\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\beta {\parallel u\parallel }_{L^1\left(\mathsf{\Omega }\right)}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& Ly=u+y_r\mathrm{in}\mathsf{\Omega },\hfill \\ \hfill & y=0\mathrm{on}\partial \mathsf{\Omega },\hfill \\ \hfill & u\in U_{ad}=\left\{v\left(x\right)\right|a\le v\left(x\right)\le b,\mathrm{a}.\mathrm{e}\mathrm{on}\mathsf{\Omega }\}\subseteq U,\hfill \end{array}\right.$$

(1)

where $Y:=H_0^1\left(\mathsf{\Omega }\right),U:=L^2\left(\mathsf{\Omega }\right),\mathsf{\Omega }\subseteq {\mathbb{R}}^n(n=2,3)$ is a convex, open and bounded domain with $C^{1,1}$- or polygonal boundary; the desired state $y_d\in H^1\left(\mathsf{\Omega }\right)$ and the source term $y_r\in H^1\left(\mathsf{\Omega }\right)$ are given; parameters $\alpha ,\beta >0$, $-\infty <a<0<b<+\infty$; L is the uniformly elliptic differential operator given by

$$Ly:=-\sum _{i,j=1}^n{\partial }_{x_j}\left(a_{ij}{y}_{x_i}\right)+c_{0}y,$$

(2)

where $a_{ij},c_0\in L^{\infty }\left(\mathsf{\Omega }\right),c_0⩾0,a_{ij}=a_{ji}$ and there is a constant $\theta >0$, such that

$$\sum _{i,j=1}^n{a}_{ij}\left(x\right){\xi }_i{\xi }_j⩾\theta {\parallel \xi \parallel }^2,\mathrm{a}.\mathrm{a}.x\in \mathsf{\Omega },\forall \xi \in {\mathbb{R}}^n.$$

(3)

Due to the $L^1$-control cost, the objective function $J(y,u)$ is non-smooth. This means that the structure of the control significantly differs from that of the smooth one. To solve elliptic optimal control problems with $L^1$-control cost, Stadler et al. proposed a semi-smooth Newton method to solve elliptic optimal control problems with $L^1$-control cost [1]. Porcelli et al. proposed a semi-smooth Newton method with a robust preconditioner for different formulations of the Newton equation [5]. However, solving Newton equations is expensive, especially when the discretization is at a fine level. Thus, some efficient first-order algorithms have received much attention in recent years. In [6], Schindele and Borzì proposed an inexact accelerated proximal gradient (APG) method in function space to solve elliptic optimal control problems with $L^1$-control cost. However, the efficiency of the APG method relies closely on the step length. The backtracking approach is applied to obtain the appropriate step length, but this greatly increases the computational cost. In [7], Song et al. proposed an inexact heterogeneous ADMM (ihADMM) algorithm. Different from the classical ADMM, two different weighted inner products are utilized to define the augmented Lagrangian function for two subproblems. Moreover, the ihADMM was applied to solve PDE-constrained optimization problems with $L^2$-control cost [8] and elliptic optimal control problems with pointwise box constraints on the state [9]. For more applications of ADMM-type algorithms to solve PDE-constrained optimization problems, see [10,11,12,13]. Inspired by the efficiency of ADMM-type algorithms on PDE-constrained optimal control problems, we consider using ADMM-type algorithms to solve (1).

In classical finite-element-based algorithms, problems are always discretized and computed at a fixed grid level. As the mesh becomes finer and finer, the scale of the discretized problems will be larger and the computation cost will increase, causing a bottleneck. Thus, it is essential to develop new approaches to solve optimal control problems in an accurate and computationally efficient way. The multigrid method is a modern field of research, which started in the early 1970s. It is well-known that the multigrid method is the optimal solution to many discretized partial differential equations [14]. Different from the classical multigrid method, Deuflhard proposed a cascadic multigrid method for elliptic partial differential equations, which solves problems from the coarse grid to the fine grid [15]. In [16], the multigrid method was applied to tackle infinite dimensional non-linear partial differential equations using Newton methods. Due to the efficiency of the multigrid method, it is a natural idea to construct the multigrid type numerical method for PDE-constrained optimization problems. The classical multigrid method for the optimal control problem is designed to solve the linear algebraic systems formulated on each step of the optimization algorithm; refer to [17] for an overview. In [18], Gong et al. proposed an adaptive multilevel correction method, which solves the optimal control problems from the coarsest mesh to the finest mesh. Chen et al. proposed a strategy of gradually refining the grid, and proposed a multilevel ADMM (mADMM) algorithm for PDE-constrained optimization problems with a $L^2$-control cost in [19]. The authors proved the global convergence and the iteration complexity results $o(1/k)$ for the mADMM.

In this paper, we focus on the first optimize, then discretize [20] approach. This approach provides the freedom to discretize subproblems using different discretization schemes. Motivated by the success of the multigrid method, we combined the multigrid method and the ihADMM algorithm to numerically solve (1). At the early stage of the whole iteration process of the algorithm, using a coarse grid will not make the precision worse, but will reduce the computational cost. Thus, we propose a strategy of gradually refining the grid. Specifically, we first introduced the iterative scheme of the iADMM in function space. Then, nodal quadrature formulas were used to approximately discretize the $L^1$-norm and $L^2$-norm to ensure that discretized subproblems have decoupled forms. Finally, appropriate methods, such as Krylov-based methods, are used to solve the discretized subproblems.

The main contribution of this paper can be summarized as follows:

• We propose an efficient multilevel heterogeneous ADMM (mADMM) algorithm for solving sparse optimal control problems with $L^1$-control cost. Specifically, we first apply the iADMM in the function space as an outer optimization method. Then, the subproblems in each iteration are discretized by the finite element method. Instead of fixing the grid size before the computation process, we apply the strategy of gradually refining the grid. Moreover, we apply nodal quadrature formulas to approximately $L^1$-norm and $L^2$-norm to overcome the difficulty that $L^1$-norm does not have a decoupled form and ensure z-subproblems have closed form solutions. Finally, we use appropriate methods to inexactly solve the subproblems in each iteration. The global convergence and the iteration complexity results for the mhADMM algorithm are also established.
• We present numerical experiments to illustrate the effectiveness of the mhADMM algorithm. In addition, we compare the performance of the mhADMM with two benchmark methods: the ihADMM and the classical ADMM. Compared to the two methods, the mhADMM has an evident advantage in terms of computational time. Moreover, the numerical results regarding iterations illustrate the mesh-independent performance of the mhADMM.

The paper is organized as follows. In Section 2, we first briefly review the iteration format of the inexact ADMM algorithm in function space, then introduce the finite element approximation. A strategy of gradually refining the grid is introduced and the mhADMM algorithm is proposed. Numerical computation of the subproblems is also discussed. In Section 3, we present the convergence results of the mhADMM algorithm. Numerical results are given in Section 4, and concluding remarks are drawn in Section 5.

2. An Multilevel Heterogeneous ADMM Algorithm

In this section, we propose an efficient convergent multilevel alternating direction method of multipliers (mhADMM). Moreover, we introduce the numerical computation of the subproblems in the mhADMM algorithm.

To apply ADMM to solve (1), we introduce an artificial variable z; then, we can rewrite (1) in an equivalent form:

$$\left\{\begin{array}{cc}\hfill \underset{(y,u,z)\in Y\times U\times U}{\overset{}{min}}& J(y,u,z)=\frac{1}{2}\parallel y-y_d{\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\frac{\alpha }{2}{\parallel u\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\beta {\parallel z\parallel }_{L^1\left(\mathsf{\Omega }\right)}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& Ly=u+y_r\mathrm{in}\mathsf{\Omega },\hfill \\ \hfill & y=0\mathrm{on}\partial \mathsf{\Omega },\hfill \\ \hfill & u=z,\hfill \\ \hfill & z\in U_{ad}=\left\{v\left(x\right)\right|a\le v\left(x\right)\le b,\mathrm{a}.\mathrm{e}\mathrm{on}\mathsf{\Omega }\}\subseteq U.\hfill \end{array}\right.$$

(4)

For the existence and uniqueness of the solution of the elliptic PDE involved in (1), where L is defined by (2), the following proposition holds.

Proposition 1. 

For every $u\in L^2\left(\mathsf{\Omega }\right)$ and $y_r\in H^1\left(\mathsf{\Omega }\right)$, the elliptic PDE involved in (1):

$$\begin{array}{cc}\hfill Ly& =u+y_r\mathrm{in}\mathsf{\Omega },\hfill \\ \hfill y& =0\mathrm{on}\partial \mathsf{\Omega },\hfill \end{array}$$

(5)

has a unique weak solution $y=y\left(u\right):=S(u+y_r)$, where $S:L^2\left(\mathsf{\Omega }\right)\to H_0^1\left(\mathsf{\Omega }\right)$ denotes the solution operator. Moreover, S is a well-defined continuous linear injective operator. The adjoint operator $S^{*}:H^{-1}\left(\mathsf{\Omega }\right)\to H_0^1\left(\mathsf{\Omega }\right)$ is also a continuous linear operator.

Proof. 

The weak formulation of (5) is given by the following:

$$\mathrm{Find}y\in H_0^1\left(\mathsf{\Omega }\right):a\left(y,v\right)={\left(u+y_r,v\right)}_{L^2\left(\mathsf{\Omega }\right)},\forall v\in H_0^1\left(\mathsf{\Omega }\right),$$

where $a:V\times V\to \mathbb{R}$ is a bilinear form:

$$a\left(y,v\right):={\int }_{\mathsf{\Omega }}\left(\sum _{i,j=1}^n{a}_{ij}{y}_{x_i}{v}_{x_j}+c_{0}yv\right)dx.$$

We know from the assumption (3) that $a\left(·,·\right)$ is symmetric, i.e., $a\left(y,v\right)=a\left(v,y\right)$ for all $y,v\in H_0^1\left(\mathsf{\Omega }\right)$. Thus, $\left(y,v\right):=a\left(y,v\right)$ defines a new inner product on $H_0^1\left(\mathsf{\Omega }\right)$. Then, the existence of a unique solution of (5) directly follows from the Riesz representation theorem. Moreover,

$$\begin{array}{c}\hfill {\theta \parallel y\parallel }_{H_0^1\left(\mathsf{\Omega }\right)}^2\le a\left(y,y\right)={\left(u+y_r,y\right)}_{L^2\left(\mathsf{\Omega }\right)}\le \parallel u+y_r{\parallel }_{L^2\left(\mathsf{\Omega }\right)}{\parallel y\parallel }_{L^2\left(\mathsf{\Omega }\right)}\le \parallel u+y_r{\parallel }_{L^2\left(\mathsf{\Omega }\right)}{\parallel y\parallel }_{H_0^1\left(\mathsf{\Omega }\right)},\end{array}$$

where $\theta$ is a constant, depending only on $\mathsf{\Omega }$, and we use the equivalence of norms in $H_0^1\left(\mathsf{\Omega }\right)$ in the first inequality. Then, we have

$${\parallel y\parallel }_{H_0^1\left(\mathsf{\Omega }\right)}^2\le \frac{1}{\theta }{\parallel u+y_r\parallel }_{L^2\left(\mathsf{\Omega }\right)}.$$

Thus, the solution operator S is well-defined and called the control-to-state mapping, which is a continuous linear injective operator. Since $H_0^1$ is a Hilbert space, the adjoint operator $S^{*}:H^{-1}\left(\mathsf{\Omega }\right)\to H_0^1\left(\mathsf{\Omega }\right)$ is also a continuous linear operator.    □

Since (4) is continuous and strongly convex, the solution of (4) exists and is unique. The following Karush–Kuhn–Tucker (KKT) conditions hold at the optimal solution of (4).

Theorem 1. 

(First-order optimality condition) $(y^{*},u^{*},z^{*})$ is the optimal solution of (4), if, and only if, adjoint state $p^{*}\in H_0^1\left(\mathsf{\Omega }\right)$ and Lagrange multiplier ${\lambda }^{*}\in L^2\left(\mathsf{\Omega }\right)$ exists, such that the following conditions hold in the weak sense:

$$\begin{array}{cc}\hfill & y^{*}=\mathcal{S}\left(u^{*}+y_r\right),\hfill \\ \hfill & p^{*}={\mathcal{S}}^{*}\left(y^{*}-y_d\right),\hfill \\ \hfill & \alpha u^{*}-p^{*}+{\lambda }^{*}=0,\hfill \\ \hfill & u^{*}=z^{*},\hfill \\ \hfill & z^{*}\in U_{ad},\hfill \\ \hfill & {〈-{\lambda }^{*},\tilde{z}-z^{*}〉}_{L^2\left(\mathsf{\Omega }\right)}+\beta \left(\parallel \tilde{z}{\parallel }_{L^1\left(\mathsf{\Omega }\right)}-{∥z^{*}∥}_{L^1\left(\mathsf{\Omega }\right)}\right)\ge 0,\forall \tilde{z}\in U_{ad}.\hfill \end{array}$$

Moreover, we have

$$u^{*}:={\mathsf{\Pi }}_{U_{ad}}\left(\frac{1}{\alpha }soft\left(p^{*},\beta \right)\right),$$

where the projection operator ${\mathsf{\Pi }}_{U_{ad}}\left(·\right)$ and the soft thresholding operator $\mathrm{soft}\left(·\right)$, respectively, are defined as follows:

$$\begin{array}{cc}\hfill {\mathsf{\Pi }}_{U_{ad}}\left(v\left(x\right)\right)& :=max\{a,min\{v\left(x\right),b\left\}\right\},\hfill \\ \hfill soft\left(v\right(x),\beta )& :=sgn\left(v\right(x\left)\right)\circ max\left(\right|v\left(x\right)|-\beta ,0).\hfill \end{array}$$

Using the solution operator S, (4) can be equivalently rewritten as the following reduced form:

$$\left\{\begin{array}{cc}\hfill \underset{(u,z)\in U\times U}{\overset{}{min}}& f\left(u\right)+g\left(z\right)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& u=z,\hfill \end{array}\right.$$

(6)

where

$$\begin{array}{cc}\hfill & f\left(u\right):=\frac{1}{2}\parallel S(u+y_r)-y_d{\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\frac{\alpha }{2}{\parallel u\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2,\hfill \\ \hfill & g\left(z\right):={\beta \parallel z\parallel }_{L^1\left(\mathsf{\Omega }\right)}+{\delta }_{U_{ad}}\left(z\right),\hfill \end{array}$$

${\delta }_{U_{ad}}$ denotes the indicator function of $U_{ad}$, i.e.,

$${\delta }_{U_{ad}}\left(z\right)=\left\{\begin{array}{cc}0,\hfill & z\in U_{ad},\hfill \\ \infty ,\hfill & z\notin U_{ad}.\hfill \end{array}\right.$$

Notice that (6) is a two-block separable convex optimization problem with linear equality constraints; the ADMM algorithm can be used to solve (6). The augmented Lagrangian function of (6) is defined as follows:

$$L_{\sigma }(u,z,\lambda ;\sigma )=f\left(u\right)+g\left(z\right)+{\langle \lambda ,u-z\rangle }_{L^2\left(\mathsf{\Omega }\right)}+{\displaystyle \frac{\sigma }{2}}{\parallel u-z\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2,$$

where $\lambda \in L^2\left(\mathsf{\Omega }\right)$ denotes the Lagrangian multiplier; $\sigma >0$ is the penalty parameter. Given the initial point $(u^0,z^0,{\lambda }^0)\in L^2\left(\mathsf{\Omega }\right)\times \mathrm{dom}\left({\delta }_{U_{ad}}(·)\right)\times L^2\left(\mathsf{\Omega }\right)$, parameters $\sigma >0$, $\tau \in \left(0,{\displaystyle \frac{\sqrt{5}+1}{2}}\right)$, the iteration format of ADMM in function space is as follows:

$$\left\{\begin{array}{c}\mathrm{Step}1:{\overline{u}}^{k+1}=\underset{u}{argmin}f\left(u\right)+{\langle {\overline{\lambda }}^k,u-{\overline{z}}^k\rangle }_{L^2\left(\mathsf{\Omega }\right)}+{\displaystyle \frac{\sigma }{2}}{\parallel u-{\overline{z}}^k\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2,\hfill \\ \mathrm{Step}2:{\overline{z}}^{k+1}=\underset{z}{argmin}g\left(z\right)+{\langle {\overline{\lambda }}^k,{\overline{u}}^{k+1}-z\rangle }_{L^2\left(\mathsf{\Omega }\right)}+{\displaystyle \frac{\sigma }{2}}{\parallel {\overline{u}}^{k+1}-z\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2,\hfill \\ \mathrm{Step}3:{\overline{\lambda }}^{k+1}={\overline{\lambda }}^k+\tau \sigma ({\overline{u}}^{k+1}-{\overline{z}}^{k+1}).\hfill \end{array}\right.$$

(7)

However, computing the exact solution of each subproblem is usually expensive and unnecessary. In [7], Song et al. propose an inexact ADMM (iADMM) in function space for (6). Krylov-based methods are used to solve the subproblems, which are equivalent to large-scale linear systems. We show the iterative scheme of the iADMM algorithm in Algorithm 1.

It is easy to see that the z-subproblem has a closed-form solution

$$\begin{array}{cc}\hfill z^{k+1}& =\underset{z}{argmin}{\displaystyle \frac{\sigma }{2}}\parallel z-{\displaystyle \frac{1}{\sigma }}(\sigma u^{k+1}+{\lambda }^k){\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\beta {\parallel z\parallel }_{L^1\left(\mathsf{\Omega }\right)}+{\delta }_{U_{ad}}\left(z\right)\hfill \\ \hfill & ={\mathsf{\Pi }}_{U_{ad}}\left(\frac{1}{\sigma }soft\left(\sigma u^{k+1}+{\lambda }^k,\beta \right)\right).\hfill \end{array}$$

For the global convergence results and the iteration complexity for Algorithm 1, we have the following theorem.

Algorithm 1 Inexact ADMM (iADMM) algorithm for (6)

Input: Choose the initial point $(u^0,z^0,{\lambda }^0)\in L^2\left(\mathsf{\Omega }\right)\times \mathrm{dom}\left({\delta }_{U_{ad}}(·)\right)\times L^2\left(\mathsf{\Omega }\right)$, parameters $\sigma >0$, $\tau \in \left(0,{\displaystyle \frac{\sqrt{5}+1}{2}}\right)$. Let ${\left\{{ϵ}_{k+1}\right\}}_{k=0}^{\infty }$ be a sequence satisfying ${\left\{{ϵ}_{k+1}\right\}}_{k=0}^{\infty }\subseteq [0,+\infty )$ and ${\displaystyle \sum _{k=0}^{\infty }{ϵ}_k<\infty }$. Set $k=0$. Output:$u^k,z^k,{\lambda }^k$. Step 1 Compute $u^{k+1}$ as an approximation solution of $$\begin{array}{c}\hfill \underset{u}{min}f\left(u\right)+{\langle {\lambda }^k,u-z^k\rangle }_{L^2\left(\mathsf{\Omega }\right)}+{\displaystyle \frac{\sigma }{2}}{\parallel u-z^k\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2\end{array}$$        such that the residual ${\delta }_u^{k+1}:=\nabla f\left(u^{k+1}\right)+{\lambda }^k+\sigma (u^{k+1}-z^k)$ satisfies        $\parallel {\delta }_u^{k+1}{\parallel }_{L^2\left(\mathsf{\Omega }\right)}\le {ϵ}_{k+1}.$ Step 2 Compute $z^{k+1}$ as follows: $$\begin{array}{c}\hfill z^{k+1}=\underset{z}{argmin}g\left(z\right)+{\langle {\lambda }^k,u^{k+1}-z\rangle }_{L^2\left(\mathsf{\Omega }\right)}+{\displaystyle \frac{\sigma }{2}}{\parallel u^{k+1}-z\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2.\end{array}$$ Step 3 Compute $${\lambda }^{k+1}={\lambda }^k+\tau \sigma (u^{k+1}-z^{k+1}).$$ Step 4 If a termination criterion is met, stop; otherwise, set $k:=k+1$ and go to Step 1.

Theorem 2. 

([7] Theorem 3) Let $(y^{*},u^{*},z^{*},p^{*},{\lambda }^{*})$ be the KKT point of (4), the sequence ${\left\{(u^k,z^k,{\lambda }^k)\right\}}_{k=0}^{\infty }$ is generated by Algorithm 1 with the associated state ${\left\{y^k\right\}}_{k=0}^{\infty }$ and adjoint state ${\left\{p^k\right\}}_{k=0}^{\infty }$; then, we have

$$\begin{array}{cc}\hfill & \underset{k\to \infty }{lim}\{\parallel u^k-u^{*}{\parallel }_{L^2\left(\mathsf{\Omega }\right)}+\parallel z^k-z^{*}{\parallel }_{L^2\left(\mathsf{\Omega }\right)}+\parallel {\lambda }^k-{\lambda }^{*}{\parallel }_{L^2\left(\mathsf{\Omega }\right)}\}=0,\hfill \\ \hfill & \underset{k\to \infty }{lim}\{\parallel y^k-y^{*}{\parallel }_{H_0^1\left(\mathsf{\Omega }\right)}+\parallel p^k-p^{*}{\parallel }_{H_0^1\left(\mathsf{\Omega }\right)}\}=0.\hfill \end{array}$$

Moreover, a constant $C_0$ only depends on the initial point $(u^0,z^0,{\lambda }^0)$ and the optimal solution $(u^{*},z^{*},{\lambda }^{*})$, such that, for $k\ge 1$,

$$\underset{1\le i\le k}{min}R(u^i,z^i,{\lambda }^i)\le \frac{C_0}{k},\underset{k\to \infty }{lim}(k·\underset{1\le i\le k}{min}R(u^i,z^i,{\lambda }^i))=0,$$

where the function $R:(u,z,\lambda )\to [0,\infty )$ is defined as

$$R(u,z,\lambda ):={\parallel \nabla f\left(u\right)+\lambda \parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+{\mathrm{dist}}^2(0,-\lambda +\partial g\left(z\right))+{\parallel u-z\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2.$$

2.1. The mhADMM Algorithm

To numerically solve (6), we consider the full discretization, in which both the state y and the control u are discretized by piecewise, linear, globally continuous finite elements. We introduce a family of regular and quasi-uniform triangulations $\left\{{\mathcal{T}}_h\right\}$ of $\overline{\mathsf{\Omega }}$, i.e., $\overline{\mathsf{\Omega }}={\bigcup }_{T\in {\mathcal{T}}_h}\overline{T}$. With each element $T\in {\mathcal{T}}_h$, we define the diameter of the set T by ${\rho }_T:=\mathrm{diam}T$ and let ${\sigma }_T$ denote the diameter of the largest ball contained in T. The grid size is defined by $h:={\max }_{T\in {\mathcal{T}}_h}{\rho }_T$. We suppose the following standard assumption holds (see [20]).

Assumption 1.

(Regular and quasi-uniform triangulations) There are two positive constants κ and τ such that

$$\frac{{\rho }_T}{{\sigma }_T}\le \kappa ,\frac{h}{{\rho }_T}\le \tau$$

hold for all $T\in {\mathcal{T}}_h$ and all $h>0$. Moreover, let us define ${\overline{\mathsf{\Omega }}}_h={\bigcup }_{T\in {\mathcal{T}}_h}\overline{T}$ and let ${\mathsf{\Omega }}_h\subseteq \mathsf{\Omega }$ and ${\Gamma }_h$ denote its interior and boundary, respectively. In the case that Ω is a convex polyhedral domain, we have $\mathsf{\Omega }={\mathsf{\Omega }}_h$. In the case that Ω is a domain with a $C^{1,1}$- boundary Γ, we assume that ${\overline{\mathsf{\Omega }}}_h$ is convex and all boundary vertices of ${\overline{\mathsf{\Omega }}}_h$ are contained in Γ, such that

$$|\mathsf{\Omega }\setminus {\mathsf{\Omega }}_h|\le c^{*}{h}^2,$$

where $|·|$ denotes the measure of the set, and $c^{*}>0$ is a constant.

Due to the homogeneous boundary condition of the state equation, we use

$$\begin{array}{cc}\hfill & Y_h:=\{y_h\in C\left(\overline{\mathsf{\Omega }}\right)|y_{h|T}\in {\mathcal{P}}_1,\forall T\in {\mathcal{T}}_h,y_h=0\mathrm{in}\overline{\mathsf{\Omega }}\setminus {\mathsf{\Omega }}_h\},\hfill \\ \hfill & U_h:=\{u_h\in C\left(\overline{\mathsf{\Omega }}\right)|u_{h|T}\in {\mathcal{P}}_1,\forall T\in {\mathcal{T}}_h,u_h=0\mathrm{in}\overline{\mathsf{\Omega }}\setminus {\mathsf{\Omega }}_h\},\hfill \end{array}$$

as the discretized state space and the discretized control space, respectively. ${\mathcal{P}}_1$ denotes the space of polynomials of degree less than or equal to 1. For the given regular and quasi-uniform triangulation ${\mathcal{T}}_h$ with nodes ${\left\{x_i\right\}}_{i=1}^{N_h}$, let ${\left\{{\varphi }_i\left(x\right)\right\}}_{i=1}^{N_h}$ be a basis of $Y_h$, $U_h$, which satisfies the following properties:

$${\varphi }_i\left(x\right)⩾0,{\parallel {\varphi }_i\left(x\right)\parallel }_{\infty }=1\forall i=1,...,N_h,\sum _{i=1}^{N_h}{\varphi }_i\left(x\right)=1.$$

Then, $u_h\in U_h$, $y_h\in Y_h$ can be represented in the following forms, respectively,

$$y_h=\sum _{i=1}^{N_h}{y}_i{\varphi }_i,u_h=\sum _{i=1}^{N_h}{u}_i{\varphi }_i,$$

where $y_i:=y_h\left(x_i\right)$, $u_i:=u_h\left(x_i\right)$. Moreover, $y_h$ can be expressed by $y_h\left(u\right)=S_h(u+y_r)$, where $S_h$ denotes the discretized version of the solution operator S. Let $U_{ad,h}$ denotes the discretized feasible set, which is defined by

$$U_{ad,h}:=U_h\cap U_{ad}=\left\{z_h=\sum _{i=1}^{N_h}{z}_i{\varphi }_i|a\le z_i\le b,\forall i=1,\cdots ,N_h\right\}\subset U_{ad}.$$

To overcome the difficulty that the discretized $L^1$-norm does not have a decoupled form, we choose the nodal quadrature formulas introduced in [21] to approximately discretized the $L^1$-norm:

$${∥z_h∥}_{L_h^1\left({\mathsf{\Omega }}_h\right)}=\sum _{i=1}^n\left|z_i\right|{\int }_{{\mathsf{\Omega }}_h}{\varphi }_i\left(x\right)\mathrm{d}x.$$

Moreover, in order to obtain a closed form solution for the z-subproblem, a similar quadrature formulae introduced in [8] is also used to discretize the $L^2$-norm in the z-subproblem:

$${∥z_h∥}_{L_h^2\left({\mathsf{\Omega }}_h\right)}={\left(\sum _{i=1}^n{\left(z_i\right)}^2{\int }_{{\mathsf{\Omega }}_h}{\varphi }_i\left(x\right)\mathrm{d}x\right)}^{\frac{1}{2}}.$$

For the given $y_r\in H^1\left(\mathsf{\Omega }\right)$ and $u\in L^2\left(\mathsf{\Omega }\right)$, the unique discretized state $y_h$ associated with u is can be expressed by $y_h\left(u\right)=S_h(u+y_r)$, where $S_h$ is the discretized version of the solution operator S. Then, we have the well-known error estimates:

Lemma 1 

([22], Theorem 4.4.6). For a given $u\in L^2\left(\mathsf{\Omega }\right)$, let y be the unique weak solution of the state Equation (5) and $y_h$ be the unique discretized state. Then, there exists a constant $c>0$ independent of h, u and $y_r$, such that

$$\parallel y-y_h{\parallel }_{L^2\left(\mathsf{\Omega }\right)}+h\parallel \nabla y-\nabla y_h{\parallel }_{L^2\left(\mathsf{\Omega }\right)}\le c{h}^2{(\parallel u\parallel }_{L^2\left(\mathsf{\Omega }\right)}+\parallel y_r{\parallel }_{L^2\left(\mathsf{\Omega }\right)}).$$

In particular, this implies $\parallel S-S_h{\parallel }_{\mathcal{L}(L^2,L^2)}\le c{h}^2$ and $\parallel S-S_h{\parallel }_{\mathcal{L}(L^2,H^1)}\le ch$.

To project the solution obtained on the coarser grid to the finer grid, we introduce the definition of the node interpolation operator $I_h$.

Definition 1. 

For a given regular and quasi-uniform triangulation ${\mathcal{T}}_h$ of Ω with nodes ${\left\{x_i\right\}}_{i=1}^{N_h}$, let ${\left\{{\varphi }_i\left(x\right)\right\}}_{i=1}^{N_h}$ denotes the standard set of nodal basis functions. The interpolation operator $I_h$ is defined as

$$I_{h}w\left(x\right):=\sum _{i=1}^{N_h}w\left(x_i\right){\varphi }_i\left(x\right)\mathrm{for}\mathrm{any}w\in C^0\left(\mathsf{\Omega }\right)\cap H^1\left(\mathsf{\Omega }\right).$$

For the interpolation error estimate, we have the following Theorem 3.

Theorem 3 

([19] Theorem 2). For all $w\in C^0\left(\mathsf{\Omega }\right)\cap H^1\left(\mathsf{\Omega }\right)$, we have

$$\parallel w-I_h{w\parallel }_{L^2\left(\mathsf{\Omega }\right)}\le c_{I}h{\parallel w\parallel }_{H^1\left(\mathsf{\Omega }\right)},$$

where $c_I$ is a constant independent of h.

At the early stage of the whole process, computing on the coarser grid can reduce the computation cost without making the precision worse. While as the iteration process proceeds, the iteration precision is supposed to increase. In this case, it is necessary to use the finer grid at the late stage. Thus, we apply the strategy of gradually refining the grid. In the initial iteration, we obtained a solution on the coarse grid, then projected the obtained solution to the finer grid. For the convenience of representing subproblems on different grids, we define

$$\begin{array}{cc}\hfill & f_{h_{k+1}}\left(u\right):=\frac{1}{2}\parallel S_{h_{k+1}}(u+I_{h_{k+1}}{y}_r)-I_{h_{k+1}}{y}_d{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}^2+\frac{\alpha }{2}{\parallel u\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}^2,\hfill \\ \hfill & g_{h_{k+1}}\left(z\right):=\beta {\parallel z\parallel }_{L_{h_{k+1}}^1\left({\mathsf{\Omega }}_{h_{k+1}}\right)}+{\delta }_{U_{ad},h_{k+1}}\left(z\right).\hfill \end{array}$$

Moreover, let

$${\displaystyle I_h{y}_r:=\sum _{i=1}^{N_h}{y}_r^i{\varphi }_i,I_h{y}_d:=\sum _{i=1}^{N_h}{y}_d^i{\varphi }_i}$$

denotes the $L^2$-projection of $y_r$ and $y_d$ onto $Y_h$, respectively. Then, we show the iterative scheme of the multilevel heterogeneous ADMM alternating direction method of multipliers (mhADMM) in Algorithm 2.

Algorithm 2 Multilevel heterogeneous ADMM (mhADMM) algorithm for (6)

Input: Choose the initial point $(u_{h_1}^0,z_{h_1}^0,{\lambda }_{h_1}^0)\in H^1\left(\mathsf{\Omega }\right)\times H^1\left(\mathsf{\Omega }\right)\times H^1\left(\mathsf{\Omega }\right)$, parameters $\sigma >0$, $\tau \in \left(0,{\displaystyle \frac{\sqrt{5}+1}{2}}\right)$. Let ${\left\{{ϵ}_{k+1}\right\}}_{k=0}^{\infty }$ be a sequence satisfying ${\left\{{ϵ}_{k+1}\right\}}_{k=0}^{\infty }\subseteq [0,+\infty )$ and ${\displaystyle \sum _{k=0}^{\infty }{ϵ}_{k+1}<\infty }$, mesh sizes ${\left\{h_{k+1}\right\}}_{k=0}^{\infty }$ of each level satisfying ${\displaystyle \sum _{k=0}^{\infty }{h}_{k+1}<\infty .}$ Set $k=0$. Output:$u_{h_k}^k,z_{h_k}^k,{\lambda }_{h_k}^k$. Step 1 Compute $u_{h_{k+1}}^{k+1}$ as an approximation solution of $$\begin{array}{c}\hfill \underset{u}{min}{f}_{h_{k+1}}\left(u\right)+{\langle {\lambda }_{h_{k+1}}^k,u-z_{h_{k+1}}^k\rangle }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}+{\displaystyle \frac{\sigma }{2}}{\parallel u-z_{h_{k+1}}^k\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}^2\end{array}$$        such that the residual ${\delta }_{u,h_{k+1}}^{k+1}:=\nabla f_{h_{k+1}}\left(u_{h_{k+1}}^{k+1}\right)+{\lambda }_{h_{k+1}}^k+\sigma (u_{h_{k+1}}^{k+1}-z_{h_{k+1}}^k)$ satisfies $\parallel {\delta }_{u,h_{k+1}}^{k+1}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}\le {ϵ}_{k+1}.$ Step 2 Compute $z_{h_{k+1}}^{k+1}$ as follows: $$\begin{array}{cc}\hfill z_{h_{k+1}}^{k+1}& =\underset{z}{argmin}{g}_{h_{k+1}}\left(z\right)+{\langle {\lambda }_{h_{k+1}}^k,u_{h_{k+1}}^{k+1}-z\rangle }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}+{\displaystyle \frac{\sigma }{2}}{\parallel u_{h_{k+1}}^{k+1}-z\parallel }_{L_{h_{k+1}}^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}^2.\hfill \end{array}$$ Step 3 Compute $${\lambda }_{h_{k+1}}^{k+1}={\lambda }_{h_{k+1}}^k+\tau \sigma (u_{h_{k+1}}^{k+1}-z_{h_{k+1}}^{k+1}).$$ Step 4 If a termination criterion is met, stop; otherwise, set $k:=k+1$ and go to Step 1.

It is easy to see that the z-subproblem has a closed-form solution:

$$\begin{array}{cc}\hfill z_{h_{k+1}}^{k+1}& =\underset{z}{argmin}{\displaystyle \frac{\sigma }{2}}\parallel z-{\displaystyle \frac{1}{\sigma }}(\sigma u_{h_{k+1}}^{k+1}+{\lambda }_{h_{k+1}}^k){\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}^2+\beta {\parallel z\parallel }_{L_{h_{k+1}}^1\left({\mathsf{\Omega }}_{h_{k+1}}\right)}+{\delta }_{U_{ad},h_{k+1}}\left(z\right)\hfill \\ \hfill & ={\mathsf{\Pi }}_{U_{ad},h_{k+1}}\left(\frac{1}{\sigma }soft\left(\sigma u_{h_{k+1}}^{k+1}+{\lambda }_{h_{k+1}}^k,\beta \right)\right).\hfill \end{array}$$

2.2. Numerical Computation of the Subproblems in Algorithm 2

To rewrite the subproblems into matrix-vector forms, we define the following matrices

$$\begin{array}{cc}\hfill & K_h:={\left(a({\varphi }_i,{\varphi }_j)\right)}_{i,j=1}^{N_h},\hfill \\ \hfill & M_h:={\left({\int }_{{\mathsf{\Omega }}_h}{\varphi }_i{\varphi }_{j}dx\right)}_{i,j=1}^{N_h},\hfill \\ \hfill & W_h:=diag{\left({\int }_{{\mathsf{\Omega }}_h}{\varphi }_i\left(x\right)\mathrm{d}x\right)}_{i,j=1}^{N_h},\hfill \end{array}$$

where $K_h$, $M_h$ and $W_h$ denote the finite element stiffness matrix, mass matrix and lump mass matrix, respectively.

For ${\displaystyle u_h=\sum _{i=1}^{N_h}{u}_i{\varphi }_i\in U_h}$, ${\displaystyle y_h=\sum _{i=1}^{N_h}{y}_i{\varphi }_i\in Y_h}$, let

$${\mathbf{u}}_h=(u_1,...,u_{N_h}),{\mathbf{y}}_h=(y_1,...,y_{N_h})$$

be the relative coefficient vectors, respectively. For ${\displaystyle I_h{y}_r=\sum _{i=1}^{N_h}{y}_r^i{\varphi }_i}$, ${\displaystyle I_h{y}_d=\sum _{i=1}^{N_h}{y}_d^i{\varphi }_i}$, let

$${\mathbf{y}}_{r,h}=(y_r^1,y_r^2,...,y_r^{N_h}),{\mathbf{y}}_{d,h}=(y_d^1,y_d^2,...,y_d^{N_h})$$

be the coefficient vectors, respectively. Let ${\mathbf{I}}_h$ denotes the vector version of the interpolation operator. Moreover, we define

$$\begin{array}{cc}\hfill f\left(\mathbf{u}\right)& :={\displaystyle \frac{1}{2}}\parallel K_h^{-1}{M}_h(\mathbf{u}+{\mathbf{y}}_{r,h})-{\mathbf{y}}_{d,h}{\parallel }_{M_h}^2+{\displaystyle \frac{\alpha }{2}}{\parallel \mathbf{u}\parallel }_{M_h}^2,\hfill \\ \hfill g\left(\mathbf{z}\right)& :=\beta \parallel W_h{\mathbf{z}\parallel }_1^2+{\delta }_{{[a,b]}^{N_h}}\left(\mathbf{z}\right).\hfill \end{array}$$

Then, the matrix-vector form of Algorithm 2 is given in Algorithm 3.

Algorithm 3 Matrix-vector form of the mhADMM algorithm

Input:$({\mathbf{u}}^0,{\mathbf{z}}^0,{\mathsf{\lambda }}^0)\in {\mathbb{R}}^{N_h}\times {[a,b]}^{N_h}\times {\mathbb{R}}^{N_h}$, parameters $\sigma >0$, $\tau \in \left(0,{\displaystyle \frac{\sqrt{5}+1}{2}}\right)$. Let ${\left\{{ϵ}_{k+1}\right\}}_{k=0}^{\infty }$ be a sequence satisfying ${\left\{{ϵ}_{k+1}\right\}}_{k=0}^{\infty }\subseteq [0,+\infty )$ and ${\displaystyle \sum _{k=0}^{\infty }\frac{{ϵ}_{k+1}}{\sqrt{\parallel M_{h_{k+1}}{\parallel }_2}}<\infty }$, mesh sizes ${\left\{h_{k+1}\right\}}_{k=0}^{\infty }$ of each iteration satisfy ${\displaystyle \sum _{k=0}^{\infty }{h}_{k+1}<\infty .}$ Set $k=0$. Output:${\mathbf{u}}_{h_k}^k,{\mathbf{z}}_{h_k}^k,{\mathsf{\lambda }}_{h_k}^k$. Step 1: Compute ${\mathbf{u}}_{h_{k+1}}^{k+1}$ as an approximation solution of $$\underset{\mathbf{u}}{min}f\left(\mathbf{u}\right)+\langle M_{h_{k+1}}{\mathsf{\lambda }}_{h_{k+1}}^k,\mathbf{u}-{\mathbf{z}}_{h_{k+1}}^k\rangle +{\displaystyle \frac{\sigma }{2}}{\parallel \mathbf{u}-{\mathbf{z}}_{h_{k+1}}^k\parallel }_{M_{h_{k+1}}}^2$$        such that the residual ${\delta }_{\mathit{u},h_{k+1}}^{k+1}:=\nabla f\left({\mathbf{u}}_{h_{k+1}}^{k+1}\right)+M_{h_{k+1}}{\mathsf{\lambda }}_{h_{k+1}}^k+\sigma M_{h_{k+1}}({\mathbf{u}}_{h_{k+1}}^{k+1}-{\mathbf{z}}_{h_{k+1}}^k)$        satisfies $\parallel {\delta }_{\mathit{u},h_{k+1}}^{k+1}\parallel \le \frac{{ϵ}_{k+1}}{\sqrt{\parallel M_{h_{k+1}}{\parallel }_2}}$. Step 2: Compute ${\mathbf{z}}_{h_{k+1}}^{k+1}$ as follows: $${\mathbf{z}}_{h_{k+1}}^{k+1}=\underset{z}{argmin}g\left(\mathbf{z}\right)+\langle M_h{\mathsf{\lambda }}_{h_{k+1}}^k,{\mathbf{u}}_{h_{k+1}}^{k+1}-\mathbf{z}\rangle +{\displaystyle \frac{\sigma }{2}}{\parallel {\mathbf{u}}_{h_{k+1}}^{k+1}-\mathbf{z}\parallel }_{W_h}^2.$$ Step 3: Compute $${\mathsf{\lambda }}_{h_{k+1}}^{k+1}={\mathbf{I}}_{h_{k+1}}{\mathsf{\lambda }}_{h_k}^k+\tau \sigma ({\mathbf{u}}_{h_{k+1}}^{k+1}-{\mathbf{z}}_{h_{k+1}}^{k+1}).$$ Step 4: If a termination criterion is met, stop; otherwise, set $k:=k+1$ and go to Step 1.

The u-subproblem at the kth iteration is equivalent to the following linear system:

$$M_{h_k}{K}_{h_k}^{-1}{M}_{h_k}(K_{h_k}^{-1}{M}_{h_k}({\mathbf{u}}_{h_k}^k+{\mathbf{y}}_{r,h_k})-{\mathbf{y}}_{d,h_k})+\alpha M_{h_k}{\mathbf{u}}_{h_k}^k+{\mathsf{\lambda }}_{h_k}^k+\sigma ({\mathbf{u}}_{h_k}^k-{\mathbf{z}}_{h_k}^k)=0,$$

(8)

where ${\mathbf{y}}_{h_k}^k:=K_{h_k}^{-1}{M}_{h_k}({\mathbf{u}}_{h_k}^k+{\mathbf{y}}_{r,h_k}),{\mathbf{p}}_{h_k}^k:=K_{h_k}^{-1}{M}_{h_k}({\mathbf{y}}_{d,h_k}-{\mathbf{y}}_{h_k}^k)$ denote the discretized state and the discretized adjoint state, respectively. Then (8) can be rewritten as:

$$\left. [\begin{array}{ccc}{M}_{h_k}& 0& K_{h_k}\\ 0& (\alpha +\sigma )M_{h_k}& -M_{h_k}\\ K_{h_k}& -M_{h_k}& 0\end{array}\right]\begin{array}{c}\left. [\begin{array}{c}{\mathbf{y}}_{h_k}^k\\ {\mathbf{u}}_{h_k}^k\\ {\mathbf{p}}_{h_k}^k\end{array}\right]\end{array}=\left. [\begin{array}{c}{M}_{h_k}{\mathbf{y}}_{d,h_k}\\ M_{h_k}(\sigma {\mathbf{I}}_{h_k}{\mathbf{z}}_{h_{k-1}}^{k-1}-{\mathbf{I}}_{h_k}{\mathsf{\lambda }}_{h_{k-1}}^{k-1})\\ M_{h_k}{\mathbf{y}}_{r,h_k}\end{array}\right].$$

(9)

We know from (9) that ${\mathbf{p}}_{h_k}^k=(\alpha +\sigma ){\mathbf{u}}_{h_k}^k-\sigma {\mathbf{I}}_{h_k}{\mathbf{z}}_{h_{k-1}}^{k-1}+{\mathbf{I}}_{h_k}{\mathsf{\lambda }}_{h_{k-1}}^{k-1}$. By eliminating the variable ${\mathbf{p}}_{h_k}^k$, (9) can be rewritten in the following reduced form without any additional computational cost:

$$\left. [\begin{array}{cc}{M}_{h_k}& (\alpha +\sigma )K_{h_k}\\ -K_{h_k}& M_{h_k}\end{array}\right]\begin{array}{c}\left. [\begin{array}{c}{\mathbf{y}}_{h_k}^k\\ {\mathbf{u}}_{h_k}^k\end{array}\right]\end{array}=\left. [\begin{array}{c}{M}_{h_k}{\mathbf{y}}_{d,h_k}+K_{h_k}(\sigma {\mathbf{I}}_{h_k}{\mathbf{z}}_{h_{k-1}}^{k-1}-{\mathbf{I}}_{h_k}{\mathsf{\lambda }}_{h_{k-1}}^{k-1})\\ -M_{h_k}{\mathbf{y}}_{r,h_k}\end{array}\right].$$

(10)

The equation system (10) can be solved by the generalized minimal residual (GMRES) with the preconditioned variant of modified hermitian and skew-hermitian splitting (PMHSS) preconditioner [23,24].

For the z-subproblem, there is a closed-form solution

$${\mathbf{z}}_{h_{k+1}}^{k+1}={\mathsf{\Pi }}_{{[a,b]}^{N_{h_{k+1}}}}\left(\frac{1}{\sigma }soft\left(\sigma {\mathbf{u}}_{h_{k+1}}^{k+1}+W_h^{-1}{M}_h{\mathsf{\lambda }}_{h_{k+1}}^k,\beta \right)\right).$$

3. Convergence Analysis

In this section, we establish the global convergence and the iteration complexity results in non-ergodic sense for the sequence generated by Algorithm 2. Before giving the convergence analysis, we first introduce the exact multi-level heterogeneous ADMM (mhADMM) algorithm. Each subproblem of the exact mhADMM algorithm is solved exactly. Given the initial point $(u^0,z^0,{\lambda }^0)\in H^1\left(\mathsf{\Omega }\right)\times H^1\left(\mathsf{\Omega }\right)\times H^1\left(\mathsf{\Omega }\right)$, parameters $\sigma >0$, $\tau \in \left(0,{\displaystyle \frac{\sqrt{5}+1}{2}}\right)$. The mesh sizes ${\left\{h_{k+1}\right\}}_{k=0}^{\infty }$ of each iteration satisfy ${\sum }_{k=0}^{\infty }{h}_{k+1}<\infty .$ Then, each iteration of the exact mhADMM has three main steps:

$$\left\{\begin{array}{c}\mathrm{Step}1:{\overline{u}}_{h_{k+1}}^{k+1}=\underset{u}{argmin}{f}_{h_{k+1}}\left(u\right)+{\langle {\overline{\lambda }}_{h_{k+1}}^k,u-{\overline{z}}_{h_{k+1}}^k\rangle }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}+{\displaystyle \frac{\sigma }{2}}{\parallel u-{\overline{z}}_{h_{k+1}}^k\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)},\hfill \\ \mathrm{Step}2:{\overline{z}}_{h_{k+1}}^{k+1}=\underset{z}{argmin}{g}_{h_{k+1}}\left(z\right)+{\langle {\overline{\lambda }}_{h_{k+1}}^k,{\overline{u}}_{h_{k+1}}^{k+1}-z\rangle }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}+{\displaystyle \frac{\sigma }{2}}{\parallel {\overline{u}}_{h_{k+1}}^{k+1}-z\parallel }_{L_{h_{k+1}}^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}^2,\hfill \\ \mathrm{Step}3:{\overline{\lambda }}_{h_{k+1}}^{k+1}={\overline{\lambda }}_{h_{k+1}}^k+\tau \sigma ({\overline{u}}_{h_{k+1}}^{k+1}-{\overline{z}}_{h_{k+1}}^{k+1}),\hfill \end{array}\right.$$

(11)

where ${\overline{\lambda }}_{h_{k+1}}^k:=I_{h_{k+1}}{\overline{\lambda }}_{h_k}^k$, ${\overline{z}}_{h_{k+1}}^k:=I_{h_{k+1}}{\overline{z}}_{h_k}^k$. Then, we use the following lemma to measure the gap between the solution sequence obtained by the ADMM in function space and the exact mhADMM algorithm in finite dimensional space.

Lemma 2. 

Let the initial point be $(z^0,{\lambda }^0)\in H^1\left(\mathsf{\Omega }\right)\times H^1\left(\mathsf{\Omega }\right)$. Let ${\left\{({\overline{u}}^k,{\overline{z}}^k,{\overline{\lambda }}^k)\right\}}_{k=0}^{\infty }$ defined in (7) be the sequence generated by the ADMM in function space and ${\left\{({\overline{u}}_{h_k}^k,{\overline{z}}_{h_k}^k,{\overline{\lambda }}_{h_k}^k)\right\}}_{k=0}^{\infty }$ defined in (11) be the sequence generated by the exact mhADMM algorithm. Then, for all $k⩾1$, we have

$$\begin{array}{cc}\hfill & \parallel {\overline{u}}^k-{\overline{u}}_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\le C_{u,k}{h}_k,\hfill \\ \hfill & \parallel {\overline{z}}^k-{\overline{z}}_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\le C_{z,k}{h}_k,\hfill \\ \hfill & \parallel {\overline{\lambda }}^{k-1}-{\overline{\lambda }}_{h_k}^{k-1}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\le C_{\lambda ,k}{h}_k,\hfill \end{array}$$

where $C_{u,k},C_{z,k},C_{\lambda ,k}$ are constants independent of $h_k$ and there is a constant C such that $C_{u,k}\le C$ for any $k\ge 1$. Thus, we have

$$\sum _{k=1}^{\infty }{\parallel {\overline{u}}^k-{\overline{u}}_{h_k}^k\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\le C\sum _{k=1}^{\infty }{h}_k.$$

Proof. 

We employ the mathematical induction to prove the conclusion. For $k=1$, we know from Theorem 3 that

$$\begin{array}{cc}\hfill \parallel {\overline{\lambda }}^0-{\overline{\lambda }}_{h_1}^0{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_1}\right)}& =\parallel {\lambda }^0-I_{h_1}{\lambda }^0{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_1}\right)}\hfill \\ \hfill & \le c_I{h}_1{\parallel {\lambda }^0\parallel }_{H^1\left({\mathsf{\Omega }}_{h_1}\right)}\hfill \\ \hfill & \le C_{\lambda ,1}{h}_1,\hfill \end{array}$$

where $C_{\lambda ,1}:=c_I{\parallel {\overline{\lambda }}^0\parallel }_{H^1\left({\mathsf{\Omega }}_{h_1}\right)}$ is a constant independent of h.

For u-subproblems, ${\overline{u}}^1$ and ${\overline{u}}_{h_1}^1$ satisfy the following optimality conditions, respectively,

$$\begin{array}{cc}\hfill & S^{*}[S({\overline{u}}^1+y_r)-y_d]+\alpha {\overline{u}}^1+{\lambda }^0+\sigma ({\overline{u}}^1-z^0)=0,\hfill \\ \hfill & S_{h_1}^{*}[S_{h_1}({\overline{u}}_{h_1}^1+I_{h_1}{y}_r)-I_{h_1}{y}_d]+\alpha {\overline{u}}_{h_1}^1+I_{h_1}{\lambda }^0+\sigma ({\overline{u}}_{h_1}^1-I_{h_1}{z}^0)=0.\hfill \end{array}$$

By subtracting the two equalities above, we have

$$\begin{array}{cc}\hfill & \left. [-(\alpha +\sigma )I-S_{h_1}^{*}{S}_{h_1}\right]({\overline{u}}^1-{\overline{u}}_{h_1}^1)\hfill \\ \hfill =& S^{*}S{\overline{u}}^1-S_{h_1}^{*}{S}_{h_1}{\overline{u}}^1+S^{*}S{y}_r-S_{h_1}^{*}{S}_{h_1}{I}_{h_1}{y}_r-S^{*}{y}_d+S_{h_1}^{*}{I}_{h_1}{y}_d+{\lambda }^0-I_{h_1}{\lambda }^0+\sigma I_{h_1}{z}^0-\sigma z^0.\hfill \end{array}$$

Then, we know from Lemma 2 in [19] that $\parallel {\overline{u}}^1-{\overline{u}}_{h_1}^1{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_1}\right)}\le C_{u,1}{h}_1$, where $C_{u,1}$ is a constant independent of $h_1$.

For z-subproblems, ${\overline{z}}^1$ and ${\overline{z}}_{h_1}^1$ satisfy

$$\begin{array}{cc}\hfill {\overline{z}}^1& ={\mathsf{\Pi }}_{U_{ad}}\left(\frac{1}{\sigma }soft\left(\sigma {\overline{u}}^1+{\lambda }^0,\beta \right)\right),\hfill \\ \hfill {\overline{z}}_{h_1}^1& ={\mathsf{\Pi }}_{U_{ad},h_1}\left(\frac{1}{\sigma }soft\left(\sigma u_{h_1}^1+{\lambda }_{h_1}^0,\beta \right)\right),\hfill \end{array}$$

respectively. Then, we know from the projection operator $\mathsf{\Pi }$ and the soft thresholding operator $\mathrm{soft}\left(·\right)$ are nonexpansive, such that

$$\begin{array}{cc}\hfill \parallel {\overline{z}}^1-{\overline{z}}_{h_1}^1{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_1}\right)}& =\parallel {\mathsf{\Pi }}_{U_{ad}}\left(\frac{1}{\sigma }soft\left(\sigma {\overline{u}}^1+{\lambda }^0,\beta \right)\right)-{\mathsf{\Pi }}_{U_{ad},h_1}\left(\frac{1}{\sigma }soft\left(\sigma {\overline{u}}_{h_1}^1+{\lambda }_{h_1}^0,\beta \right)\right){\parallel }_{L^2\left({\mathsf{\Omega }}_{h_1}\right)}\hfill \\ \hfill & \le \frac{1}{\sigma }{\parallel \sigma {\overline{u}}^1-\sigma {\overline{u}}_{h_1}^1+{\lambda }^0-{\lambda }_{h_1}^0\parallel }_{L^2\left({\mathsf{\Omega }}_{h_1}\right)}\hfill \\ \hfill & \le \parallel {\overline{u}}^1-{\overline{u}}_{h_1}^1{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_1}\right)}+\frac{1}{\sigma }{\parallel {\lambda }^0-{\lambda }_{h_1}^0\parallel }_{L^2\left({\mathsf{\Omega }}_{h_1}\right)}\hfill \\ \hfill & \le C_{u,1}{h}_1+\frac{1}{\sigma }{C}_{\lambda ,1}{h}_1\hfill \\ \hfill & =C_{z,1}{h}_1,\hfill \end{array}$$

where $C_{z,1}:=C_{u,1}+\frac{1}{\sigma }{C}_{\lambda ,1}$. Hence, the statement is true for $k=1$.

For $k>1$, we assume the statement is true for $\forall j\le k$. Then, for $j=k+1$, we have

$$\begin{array}{cc}\hfill \parallel {\overline{\lambda }}^k-{\overline{\lambda }}_{h_{k+1}}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}=& \parallel {\overline{\lambda }}^k-{\overline{\lambda }}_{h_k}^k+{\overline{\lambda }}_{h_k}^k-I_{h_{k+1}}{\overline{\lambda }}_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}\hfill \\ \hfill \le & \parallel {\overline{\lambda }}^k-{\overline{\lambda }}_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}+{\parallel {\overline{\lambda }}_{h_k}^k-I_{h_{k+1}}{\overline{\lambda }}_{h_k}^k\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}\hfill \\ \hfill \le & \parallel {\overline{\lambda }}^{k-1}-{\overline{\lambda }}_{h_k}^{k-1}+\tau \sigma ({\overline{u}}^k-{\overline{u}}_{h_k}^k)-\tau \sigma ({\overline{z}}^k-{\overline{z}}_{h_k}^k){\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\hfill \\ \hfill & +\parallel {\overline{\lambda }}_{h_k}^k-I_{h_{k+1}}{\overline{\lambda }}_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}\hfill \\ \hfill \le & \parallel {\overline{\lambda }}^{k-1}-{\overline{\lambda }}_{h_k}^{k-1}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}+\tau \sigma \left(\parallel {\overline{u}}^k-{\overline{u}}_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}+{\parallel {\overline{z}}^k-{\overline{z}}_{h_k}^k\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\right)\hfill \\ \hfill & +\parallel {\overline{\lambda }}_{h_k}^k-I_{h_{k+1}}{\overline{\lambda }}_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}\hfill \\ \hfill \le & (C_{\lambda ,k}+\tau \sigma C_{u,k}+\tau \sigma C_{z,k})h_k+c_I{h}_{k+1}{\parallel {\overline{\lambda }}_{h_k}^k\parallel }_{H^1\left({\mathsf{\Omega }}_{h_{k+1}}\right)}\hfill \\ \hfill \le & C_{\lambda ,k+1}{h}_{k+1},\hfill \end{array}$$

where $C_{\lambda ,k+1}:=C_{k+1}(C_{\lambda ,k}{h}_k+\tau \sigma C_{u,k}+\tau \sigma C_{z,k})+c_I{\parallel {\overline{\lambda }}_{h_k}^k\parallel }_{H^1\left({\mathsf{\Omega }}_{h_{k+1}}\right)}$ is a constant independent of $h_{k+1}$. In the last equality, we use the property ${\sum }_{k=0}^{\infty }{h}_k<\infty$; thus, there exists a constant $C_{k+1}$ such that $h_k<C_{k+1}{h}_{k+1}$.

For u-subproblems, ${\overline{u}}^{k+1}$ and ${\overline{u}}_{h_{k+1}}^{k+1}$ satisfy the following optimality conditions respectively,

$$\begin{array}{cc}\hfill & S^{*}[S({\overline{u}}^k+y_r)-y_d]+\alpha {\overline{u}}^{k+1}+{\overline{\lambda }}^k+\sigma ({\overline{u}}^{k+1}-{\overline{z}}^k)=0,\hfill \\ \hfill & S_{h_{k+1}}^{*}[S_{h_{k+1}}({\overline{u}}_{h_{k+1}}^{k+1}+I_{h_{k+1}}{y}_r)-I_{h_{k+1}}{y}_d]+\alpha {\overline{u}}_{h_{k+1}}^{k+1}+{\overline{\lambda }}_{h_{k+1}}^k+\sigma ({\overline{u}}_{h_{k+1}}^{k+1}-{\overline{z}}_{h_{k+1}}^k)=0.\hfill \end{array}$$

By subtracting the two equalities above, we have

$$\begin{array}{cc}\hfill & -\left. [(\alpha +\sigma )I+S_{h_{k+1}}^{*}{S}_{h_{k+1}}\right]({\overline{u}}^{k+1}-{\overline{u}}_{h_{k+1}}^{k+1})\hfill \\ \hfill =& S^{*}S{\overline{u}}^{k+1}-S_{h_{k+1}}^{*}{S}_{h_{k+1}}{\overline{u}}^{k+1}+S^{*}S{y}_r-S_{h_{k+1}}^{*}{S}_{h_{k+1}}{I}_{h_{k+1}}{y}_r-S^{*}{y}_d+S_{h_{k+1}}^{*}{I}_{h_{k+1}}{y}_d\hfill \\ \hfill & +{\overline{\lambda }}^k-{\overline{\lambda }}_{h_{k+1}}^k-\sigma ({\overline{z}}^k-{\overline{z}}_{h_{k+1}}^k).\hfill \end{array}$$

Then, we know from Lemma 2 in [19] that

$$\parallel {\overline{u}}^{k+1}-{\overline{u}}_{h_{k+1}}^{k+1}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}\le C_{u,k+1}{h}_{k+1},$$

where $C_{u,k+1}$ is a constant independent of $h_{k+1}$.

For z-subproblems, ${\overline{z}}^{k+1}$ and ${\overline{z}}_{h_{k+1}}^{k+1}$ satisfy

$$\begin{array}{cc}\hfill {\overline{z}}^{k+1}& ={\mathsf{\Pi }}_{U_{ad}}\left(\frac{1}{\sigma }soft\left(\sigma {\overline{u}}^{k+1}+{\lambda }^k,\beta \right)\right),\hfill \\ \hfill {\overline{z}}_{h_{k+1}}^{k+1}& ={\mathsf{\Pi }}_{U_{ad},h_1}\left(\frac{1}{\sigma }soft\left(\sigma {\overline{u}}_{h_{k+1}}^{k+1}+{\lambda }_{h_{k+1}}^k,\beta \right)\right),\hfill \end{array}$$

respectively. We know the projection operator $\mathsf{\Pi }$ and the soft thresholding operator $\mathrm{soft}\left(·\right)$ are nonexpansive, such that

$$\begin{array}{cc}\hfill \parallel {\overline{z}}^{k+1}-{\overline{z}}_{h_{k+1}}^{k+1}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}& \le \frac{1}{\sigma }{\parallel \sigma {\overline{u}}^{k+1}-\sigma {\overline{u}}_{h_{k+1}}^{k+1}+{\lambda }^k-{\lambda }_{h_{k+1}}^k\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}\hfill \\ \hfill & \le \parallel {\overline{u}}^{k+1}-{\overline{u}}_{h_{k+1}}^{k+1}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}+\frac{1}{\sigma }{\parallel {\lambda }^k-{\lambda }_{h_{k+1}}^k\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}\hfill \\ \hfill & \le C_{u,k+1}{h}_{k+1}+\frac{1}{\sigma }{C}_{\lambda ,k+1}{h}_{k+1}\hfill \\ \hfill & =C_{z,k+1}{h}_{k+1},\hfill \end{array}$$

where $C_{z,k+1}:=C_{u,k+1}+\frac{1}{\sigma }{C}_{\lambda ,k+1}$. Hence, the statement is true for $j=k+1$. We complete the whole proof of Lemma 2.    □

Similarly, we have the following lemma. Lemma 3 shows the gap between the sequence $(u^k,z^k,{\lambda }^k)$ generated by Algorithm 1 and the sequence $(u_{h_k}^k,z_{h_k}^k,{\lambda }_{h_k}^k)$ generated by Algorithm 2.

Lemma 3. 

Let the initial point be $(u^0,z^0,{\lambda }^0)\in H^1\left(\mathsf{\Omega }\right)\times H^1\left(\mathsf{\Omega }\right)\times H^1\left(\mathsf{\Omega }\right)$. Let ${\left\{(u^k,z^k,{\lambda }^k)\right\}}_{k=0}^{\infty }$ be the sequence generated by Algorithm 1 and ${\left\{(u_{h_k}^k,z_{h_k}^k,{\lambda }_{h_k}^k)\right\}}_{k=0}^{\infty }$ be the sequence generated by Algorithm 2. Then, for all $k⩾1$, we have

$$\begin{array}{cc}\hfill & \parallel u^k-u_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\le {\widehat{C}}_{u,k}(h_k+\parallel {\delta }_{u,h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}),\hfill \\ \hfill & \parallel z^k-z_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\le {\widehat{C}}_{z,k}(h_k+\parallel {\delta }_{u,h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}),\hfill \\ \hfill & \parallel {\lambda }^{k-1}-{\lambda }_{h_k}^{k-1}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\le {\widehat{C}}_{\lambda ,k}(h_k+\parallel {\delta }_{u,h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}),\hfill \end{array}$$

where ${\widehat{C}}_{u,k},{\widehat{C}}_{z,k},{\widehat{C}}_{\lambda ,k}$ are constants independent of $h_k$, and there exists a constant $\widehat{C}$ such that ${\widehat{C}}_{u,k}\le \widehat{C}$ for any $k\ge 1$. Thus we have

$$\sum _{k=1}^{\infty }{\parallel u^k-u_{h_k}^k\parallel }_{L^2\left(\mathsf{\Omega }\right)}\le \widehat{C}\sum _{k=1}^{\infty }\left(h_k+{\parallel {\delta }_{u,h_k}^k\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\right).$$

Proof. 

We employ the mathematical induction to prove the conclusion. The proof is similar to Lemma 2. We do not discuss this in detail here.    □

The total error of utilizing numerical methods to solve PDE-constrained optimal control problem consists of two parts: the discretization error and the iteration error. These two kinds of errors can be regarded as the error of inexactly solving infinite-dimensional subproblems. Thus, the mhADMM algorithm can be regarded as the iADMM algorithm in function space. Inspired by the results of Theorem 2, we have the following convergence results.

Theorem 4. 

Let $(y^{*},u^{*},z^{*},p^{*},{\lambda }^{*})$ be the KKT point of (1), ${(u_{h_k}^k,z_{h_k}^k,{\lambda }_{h_k}^k)}_{k=0}^{\infty }$ be the sequence generated by Algorithm 2 with the associated state ${\left\{y_{h_k}^k\right\}}_{k=0}^{\infty }$ and the adjoint state ${\left\{p_{h_k}^k\right\}}_{k=0}^{\infty }$. Then we have

$$\begin{array}{cc}\hfill & \underset{k\to \infty }{lim}\{\parallel u_{h_k}^k-u^{*}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}+\parallel z_{h_k}^k-z^{*}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}+\parallel {\lambda }_{h_k}^k-{\lambda }^{*}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\}=0,\hfill \\ \hfill & \underset{k\to \infty }{lim}\{\parallel y_{h_k}^k-y^{*}{\parallel }_{H_0^1\left({\mathsf{\Omega }}_{h_k}\right)}+\parallel p_{h_k}^k-p^{*}{\parallel }_{H_0^1\left({\mathsf{\Omega }}_{h_k}\right)}\}=0.\hfill \end{array}$$

Moreover, there exists a constant $\overline{C}$ that only depends on the initial point $(u^0,z^0,{\lambda }^0)$ and the optimal solution $(u^{*},z^{*},{\lambda }^{*})$ such that, for $k\ge 1$,

$$\underset{1\le i\le k}{min}{R}_{h_i}(u_{h_i}^i,z_{h_i}^i,{\lambda }_{h_i}^i)\le \frac{\overline{C}}{k},\underset{k\to \infty }{lim}(k·\underset{1\le i\le k}{min}{R}_{h_i}(u_{h_i}^i,z_{h_i}^i,{\lambda }_{h_i}^i))=0,$$

where $R_{h_i}:(u_{h_i}^i,z_{h_i}^i,{\lambda }_{h_i}^i)\to [0,\infty )$ is defined as

$$R_h(u_{h_i}^i,z_{h_i}^i,{\lambda }_{h_i}^i):=\parallel \nabla f_h(u_{h_i}^i)+{\lambda }_{h_i}^{i-1}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_i}\right)}^2+{\mathrm{dist}}^2(0,-{\lambda }_{h_i}^{i-1}+\partial g_h\left(z_{h_i}^i\right))+{\parallel u_{h_i}^i-z_{h_i}^i\parallel }_{L^2\left({\mathsf{\Omega }}_{h_i}\right)}^2.$$

Proof. 

Note that $(u_{h_k}^k,z_{h_k}^k,{\lambda }_{h_k}^k)$ can be regarded as the inexact solution obtained by Algorithm 1. The total error ${\delta }_u^k$ consists of two parts, the discretization error from gradually refining the grid and the iteration error from the inexactly solving the subproblems. Then, we know from the optimality conditions of the u-subproblem in Algorithm 1 that

$$S^{*}[S(u_{h_k}^k+y_r)-y_d]+\alpha u_{h_k}^k+{\lambda }_{h_{k-1}}^{k-1}+\sigma (u_{h_k}^k-z_{h_{k-1}}^{k-1})={\delta }_u^k.$$

(12)

Moreover, we know from the optimality condition of the u-subproblem in ADMM in function space, the optimality conditions of the u-subproblem in Algorithm 2 and the exact multi-level ADMM that

$$S^{*}[S({\overline{u}}^k+y_r)-y_d]+\alpha {\overline{u}}^k+{\overline{\lambda }}^{k-1}+\sigma ({\overline{u}}^k-{\overline{z}}^{k-1})=0,$$

(13)

$$S_{h_k}^{*}[S_{h_k}(u_{h_k}^k+I_{h_k}{y}_r)-I_{h_k}{y}_d]+\alpha u_{h_k}^k+{\lambda }_{h_k}^{k-1}+\sigma (u_{h_k}^k-z_{h_k}^{k-1})={\delta }_{u,h_k}^k,$$

(14)

$$S_{h_k}^{*}[S_{h_k}({\overline{u}}_{h_k}^k+I_{h_k}{y}_r)-I_{h_k}{y}_d]+\alpha {\overline{u}}_{h_k}^k+{\overline{\lambda }}_{h_k}^{k-1}+\sigma ({\overline{u}}_{h_k}^k-{\overline{z}}_{h_k}^{k-1})=0.$$

(15)

Then, we know from (12)–(15) that

$$\begin{array}{cc}\hfill {\delta }_u^k=& {\delta }_u^k-{\delta }_{u,h_k}^k+{\delta }_{u,h_k}^k\hfill \\ \hfill =& {\delta }_{u,h_k}^k+S^{*}S(u_{h_k}^k-{\overline{u}}^k)+(\alpha +\sigma )(u_{h_k}^k-{\overline{u}}^k)+({\lambda }^{k-1}-{\overline{\lambda }}^{k-1})-\sigma (z^{k-1}-{\overline{z}}^{k-1})\hfill \\ \hfill & +S_{h_k}^{*}{S}_{h_k}({\overline{u}}_{h_k}^k-u_{h_k}^k)+(\alpha +\sigma )({\overline{u}}_{h_k}^k-u_{h_k}^k)+({\overline{\lambda }}_{h_k}^{k-1}-{\lambda }_{h_k}^{k-1})-\sigma ({\overline{z}}_{h_k}^{k-1}-z_{h_k}^{k-1})\hfill \\ \hfill =& {\delta }_{u,h_k}^k+\underset{I_1}{\underbrace{[(\alpha +\sigma )I+S^{*}S]({\overline{u}}_{h_k}^k-{\overline{u}}^k)}}+\underset{I_2}{\underbrace{({\lambda }^{k-1}-{\lambda }_{h_k}^{k-1})}}-\underset{I_3}{\underbrace{\sigma (z^{k-1}-z_{h_k}^{k-1})}}\hfill \\ \hfill & +\underset{I_4}{\underbrace{({\overline{\lambda }}_{h_k}^{k-1}-{\overline{\lambda }}^{k-1})}}-\underset{I_5}{\underbrace{\sigma ({\overline{z}}_{h_k}^{k-1}-{\overline{z}}^{k-1})}}+\underset{I_6}{\underbrace{(S^{*}S-S_{h_k}^{*}{S}_{h_k})(u_{h_k}^k-{\overline{u}}_{h_k}^k)}}.\hfill \end{array}$$

(16)

For the term $I_1$ and $I_4$, we know from Lemma 2 that

$$\begin{array}{cc}\hfill \parallel I_1{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}& \le (\alpha +\sigma +\parallel S^{*}\parallel \parallel S\parallel )C_{u,k}{h}_k,\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \parallel I_4{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}& \le C_{\lambda ,k}{h}_k.\hfill \end{array}$$

(18)

For the term $I_5$, we know from Theorem 3 and Lemma 2 that

$$\begin{array}{cc}\hfill \parallel I_5{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}& =\sigma \parallel {\overline{z}}_{h_k}^{k-1}-{\overline{z}}_{h_{k-1}}^{k-1}+{\overline{z}}_{h_{k-1}}^{k-1}-{\overline{z}}^{k-1}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\hfill \\ \hfill & \le \sigma \parallel {\overline{z}}_{h_k}^{k-1}-{\overline{z}}_{h_{k-1}}^{k-1}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}+\sigma {\parallel {\overline{z}}_{h_{k-1}}^{k-1}-{\overline{z}}^{k-1}\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\hfill \\ \hfill & \le \sigma c_I{h}_k{\parallel {\overline{z}}_{h_{k-1}}^{k-1}\parallel }_{H^1\left({\mathsf{\Omega }}_{h_k}\right)}+\sigma C_{z,k-1}{h}_{k-1}\hfill \\ \hfill & \le \sigma c_I{h}_k{\parallel {\overline{z}}_{h_{k-1}}^{k-1}\parallel }_{H^1\left({\mathsf{\Omega }}_{h_k}\right)}+\sigma C_{z,k-1}{C}_k{h}_k\hfill \\ \hfill & \le c_5{h}_k,\hfill \end{array}$$

(19)

where $c_5:=\sigma c_I{\parallel {\overline{z}}_{h_{k-1}}^{k-1}\parallel }_{H^1\left({\mathsf{\Omega }}_{h_k}\right)}+\sigma C_{z,k-1}{C}_k$ is a constant, $h_{k-1}\le C_k{h}_k$.

For the term $I_2$, we know from Lemma 3 that

$$\parallel I_2{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\le {\widehat{C}}_{\lambda ,k}(h_k+\parallel {\delta }_{u,h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}).$$

(20)

For the term $I_3$, we know from Theorem 3 and Lemma 2 that

$$\parallel I_3{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\le c_3(h_k+\parallel {\delta }_{u,h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}),$$

(21)

where $c_3$ is a constant.

Finally, for the term $I_6$, we make use of the decomposition

$$\begin{array}{cc}\hfill \parallel u_{h_k}^k-{\overline{u}}_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}=& \parallel u_{h_k}^k-u^k+u^k-u^{*}+u^{*}-{\overline{u}}^k+{\overline{u}}^k-{\overline{u}}_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\hfill \\ \hfill \le & \parallel u_{h_k}^k-u^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}+\parallel u^k-u^{*}{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}+\parallel u^{*}-{\overline{u}}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}+{\parallel {\overline{u}}^k-{\overline{u}}_{h_k}^k\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\hfill \\ \hfill \le & {\widehat{C}}_{u,k}(h_k+\parallel {\delta }_{u,h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)})+C^{*},\hfill \end{array}$$

(22)

where ${\widehat{C}}_{u,k},C^{*}$ are constants in dependent of $h_k$. In the last equality, we used Lemma 2, Lemma 3, the convergence property of ADMM in function space and the inexact ADMM in function space. Then, we know from Proposition 1 and Lemma 1 that

$$\begin{array}{cc}\hfill \parallel I_6{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}=& \parallel (S^{*}S-S^{*}{S}_{h_k}+S^{*}{S}_{h_k}-S_{h_k}^{*}{S}_{h_k})(u_{h_k}^k-{\overline{u}}_{h_k}^k){\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\hfill \\ \hfill \le & \parallel (S^{*}S-S^{*}{S}_{h_k})(u_{h_k}^k-{\overline{u}}_{h_k}^k){\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}+{\parallel (S^{*}{S}_{h_k}-S_{h_k}^{*}{S}_{h_k})(u_{h_k}^k-{\overline{u}}_{h_k}^k)\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\hfill \\ \hfill \le & \parallel S^{*}\parallel \parallel S-S_{h_k}\parallel \parallel u_{h_k}^k-{\overline{u}}_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}+\parallel S^{*}-S_{h_k}^{*}\parallel \parallel S_{h_k}\parallel \parallel u_{h_k}^k-{\overline{u}}_{h_k}^k{\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}\hfill \\ \hfill \le & c_6{h}_k,\hfill \end{array}$$

(23)

where $c_6$ is a constant.

Then, we know from (16)–(23) that there are constants $C_1^{*}$ and $C_2^{*}$ such that

$$\begin{array}{c}\hfill \sum _{k=1}^{\infty }\parallel {\delta }_u^k{\parallel }_{L^2\left(\mathsf{\Omega }\right)}\le C_1^{*}\sum _{k=1}^{\infty }{\parallel {\delta }_{u,h_k}^k\parallel }_{L^2\left({\mathsf{\Omega }}_{h_k}\right)}+C_2^{*}\sum _{k=1}^{\infty }{h}_k.\end{array}$$

Moreover, the mesh sizes ${\left\{h_{k+1}\right\}}_{k=0}^{\infty }$ of each mhADMM iteration satisfy ${\displaystyle \sum _{k=0}^{\infty }{h}_{k+1}<\infty }$, the residuals of each mhADMM iteration satisfy ${\displaystyle \sum _{k=0}^{\infty }{\parallel {\delta }_{u,h_{k+1}}^{k+1}\parallel }_{L^2\left({\mathsf{\Omega }}_{h_{k+1}}\right)}\le \sum _{k=0}^{\infty }{ϵ}_{k+1}<\infty }$, thus we have

$$\sum _{k=1}^{\infty }{\parallel {\delta }_u^k\parallel }_{L^2\left(\mathsf{\Omega }\right)}<\infty .$$

Then, we know from Theorem 2 that the global convergence and the iteration complexity results $o(1/k)$ for Algorithm 2 are guaranteed.    □

4. Numerical Experiments

In this section, we illustrate the numerical performance of the mhADMM algorithm for the elliptic PDE-constrained optimization problems with $L^1$-control cost. For our numerical experiment, we used MATLAB R2021b with the FEM package iFEM [25] on a Thinkpad laptop with 2.8 GHz Intel Core i7 processor with 16GB of RAM.

In the mhADMM algorithm, the accuracy of a numerical solution is measured by the following residual. Let $ϵ$ be a given accuracy tolerance, we terminate the algorithm when $\eta <ϵ$, where $\eta :=\max \{{\eta }_1,{\eta }_2,{\eta }_3,{\eta }_4,{\eta }_5\},$ where

$$\begin{array}{cc}\hfill {\eta }_1& :=\frac{\parallel K_h\mathrm{y}-M_h\mathrm{u}-M_h{\mathrm{y}}_{\mathrm{r}}\parallel }{1+\parallel M_h{\mathrm{y}}_r\parallel },{\eta }_2:=\frac{\parallel M_h(\mathrm{u}-\mathrm{z})\parallel }{1+\parallel \mathrm{u}\parallel },{\eta }_3:=\frac{\parallel M_h(\mathrm{y}-{\mathrm{y}}_d)+K_h\mathrm{p}\parallel }{1+\parallel M_h{\mathrm{y}}_d\parallel },\hfill \\ \hfill {\eta }_4& :=\frac{\parallel \alpha M_h\mathrm{u}-M_h\mathrm{p}+M_h\lambda \parallel }{1+\parallel \mathrm{u}\parallel },{\eta }_5:=\frac{\parallel \mathrm{z}-{\mathsf{\Pi }}_{{[a,b]}^{N_h}}\left(soft\left(W_h^{-1}{M}_h\lambda ,\beta \right)\right)\parallel }{1+\parallel \mathrm{z}\parallel }.\hfill \end{array}$$

To present the finite element error estimates’ results, we introduce the experimental order of convergence, a brief EOC defined by

$$\mathrm{EOC}:=\frac{logE\left(h_1\right)-logE\left(h_2\right)}{log{h}_1-log{h}_2},$$

where $h_1,h_2>0,h_1\ne h_2$ denotes different grid sizes, E denotes the positive error functional

$$E\left(h\right):={∥u-u_h∥}_{L^2\left(\mathsf{\Omega }\right)}.$$

We note that if $E\left(h\right)=\mathcal{O}\left(h^{\gamma }\right)$ holds, then $\mathrm{EOC}\approx \gamma$.

As shown in Section 2.1, instead of using the standard piecewise linear and continuous finite elements, nodal quadrature formulas are used to approximately discretize the $L^1$-norm and $L^2$-norm in Algorithm 2. In both examples, the mhADMM algorithm, the ihADMM algorithm and the classical ADMM algorithm are employed to obtain numerical solutions of different grid sizes. For both numerical examples and all algorithms, we chose $(u^0,z^0,{\lambda }^0)=(0,0,0)$ as the initial values. The penalty parameter $\sigma$ was chosen as $\sigma =\alpha$. For the step length $\tau$, we chose $\tau =1.618$. We terminate the algorithms when the residual $\eta <{10}^{-6}$ with the maximum number of iterations set to 500.

In numerical experiments, we show the numerical results for different final mesh sizes. In

Table 1. The convergence behavior of the mhADMM, the ihADMM and the classical ADMM for Example 1.

Case	h	#dofs	E	EOC	Index	mhADMM	ihADMM	ADMM
1	$\frac{\sqrt{2}}{2^4}$	225	9.66 × 10${}^{-2}$	0.99	residual $\eta$	8.93 × 10${}^{-7}$	6.69 × 10${}^{-7}$	6.46 × 10${}^{-7}$
					time (s)	0.12	0.13	0.13
					#iter	20	16	18
2	$\frac{\sqrt{2}}{2^5}$	961	4.46 × 10${}^{-2}$	1.05	residual $\eta$	9.44 × 10${}^{-7}$	6.60 × 10${}^{-7}$	8.63 × 10${}^{-7}$
					time (s)	0.40	0.42	0.53
					#iter	20	18	24
3	$\frac{\sqrt{2}}{2^6}$	3969	1.49 × 10${}^{-2}$	1.23	residual $\eta$	3.30 × 10${}^{-7}$	7.37 × 10${}^{-7}$	9.40 × 10${}^{-7}$
					time (s)	2.16	2.67	8.51
					#iter	22	21	49
4	$\frac{\sqrt{2}}{2^7}$	16,129	4.92 × 10${}^{-3}$	1.32	residual $\eta$	5.57 × 10${}^{-7}$	5.54 × 10${}^{-7}$	8.87 × 10${}^{-7}$
					time (s)	14.86	18.35	262.56
					#iter	21	23	120
5	$\frac{\sqrt{2}}{2^8}$	65,025	1.65 × 10${}^{-3}$	1.37	residual $\eta$	7.02 × 10${}^{-7}$	4.61 × 10${}^{-7}$	9.97 × 10${}^{-7}$
					time (s)	114.33	200.43	7576.14
					#iter	20	25	257
6	$\frac{\sqrt{2}}{2^9}$	261,121	5.83 × 10${}^{-4}$	1.39	residual $\eta$	6.91 × 10${}^{-7}$	2.82 × 10${}^{-7}$	1.17 × 10${}^{-5}$
					time (s)	2457.43	3850.03	279,881.28
					#iter	20	27	500

and

Table 2. The convergence behavior of the mhADMM, the ihADMM and the classical ADMM for Example 2.

Case	h	#dofs	E	EOC	Index	mhADMM	ihADMM	ADMM
1	$\frac{\sqrt{2}}{2^4}$	225	9.69 × 10${}^{-1}$	1.19	residual $\eta$	7.12 × 10${}^{-7}$	7.12 × 10${}^{-7}$	8.04 × 10${}^{-7}$
					time (s)	0.25	0.30	0.33
					#iter	22	22	35
2	$\frac{\sqrt{2}}{2^5}$	961	5.72 × 10${}^{-1}$	0.97	residual $\eta$	7.71 × 10${}^{-7}$	7.77 × 10${}^{-7}$	9.37 × 10${}^{-7}$
					time (s)	0.95	1.10	1.66
					#iter	19	19	88
3	$\frac{\sqrt{2}}{2^6}$	3969	1.64 × 10${}^{-1}$	1.25	residual $\eta$	6.03 × 10${}^{-7}$	6.02 × 10${}^{-7}$	7.63 × 10${}^{-7}$
					time (s)	5.03	5.57	35.73
					#iter	20	20	198
4	$\frac{\sqrt{2}}{2^7}$	16,129	4.44 × 10${}^{-2}$	1.41	residual $\eta$	8.48 × 10${}^{-7}$	8.25 × 10${}^{-7}$	8.10 × 10${}^{-7}$
					time (s)	28.13	37.70	851.48
					#iter	20	20	454
5	$\frac{\sqrt{2}}{2^8}$	65,025	1.72 × 10${}^{-2}$	1.40	residual $\eta$	7.15 × 10${}^{-7}$	7.74 × 10${}^{-7}$	1.66 × 10${}^{-5}$
					time (s)	134.46	196.07	8659.08
					#iter	21	21	500
6	$\frac{\sqrt{2}}{2^9}$	261,121	-	-	residual $\eta$	9.09 × 10${}^{-7}$	9.01 × 10${}^{-7}$	6.17 × 10${}^{-5}$
					time (s)	1052.59	1972.46	151,623.99
					#iter	21	21	500

, h denotes the final mesh size, ‘#dofs’ denotes the dimension of the control variable on each grid level, and ‘iter’ represents the times of iteration. To guarantee the sequence ${\left\{{ϵ}_{k+1}\right\}}_{k=0}^{\infty }\subseteq [0,+\infty )$ satisfies ${\displaystyle \sum _{k=0}^{\infty }\frac{{ϵ}_{k+1}}{\sqrt{\parallel M_{h_{k+1}}{\parallel }_2^2}}<\infty }$, and the mesh sizes ${\left\{h_k\right\}}_{k=0}^{\infty }\subseteq [0,+\infty )$ of each iteration satisfy ${\displaystyle \sum _{k=0}^{\infty }{h}_{k+1}\le \infty }$, we choose ${ϵ}_{k+1}=\frac{C}{{(k+1)}^2}$, where C is a constant and $h_k=\frac{\sqrt{2}}{2^{k+3}},k\in \mathbb{Z},k\ge 1$ in both examples. Moreover, we would like to point out that, in the iteration of mhADMM algorithm, once the mesh size $h_k$ reaches the final mesh size h, we continue the iteration in the final mesh until the stopping criterion above is satisfied.

Before providing examples, we first introduce the following algorithm, which can help us formulate sparse optimal control problems.

According to the first-order optimality condition given in Theorem 1, it is easy to see that Algorithm 4 provides a construction strategy for problems with known optimal solutions $(y^{*},u^{*})$.

Algorithm 4 Construct the optimal control problem

Step 1 Choose $y^{*}\in H_0^1\left(\mathsf{\Omega }\right)$ and $p^{*}\in H_0^1\left(\mathsf{\Omega }\right)$ arbitrarily. Step 2 $$\begin{array}{ccc}& \hfill u^{*}:=& {\mathsf{\Pi }}_{U_{ad}}\left(\frac{1}{\alpha }soft\left(p^{*},\beta \right)\right)\hfill \\ & & =\left\{\begin{array}{c}min\left\{\frac{p^{*}-\beta }{\alpha },b\right\}, \mathrm{on}x\in \mathsf{\Omega }:p^{*}\left(x\right)>\beta ,\hfill \\ \max \left\{\frac{p^{*}+\beta }{\alpha },a\right\}, \mathrm{on}x\in \mathsf{\Omega }:p^{*}\left(x\right)<-\beta \hfill \\ 0, \mathrm{elsewhere}.\hfill \end{array}\right.\hfill \end{array}$$ Step 3 Set $y_r=S^{-1}{y}^{*}-u^{*}$ and $y_d=y^{*}-{\left(S^{*}\right)}^{-1}p$.

Example 1. 

Consider

$$\left\{\begin{array}{cc}\hfill \underset{(y,u)\in H_0^1\left(\mathsf{\Omega }\right)\times L^2\left(\mathsf{\Omega }\right)}{\overset{}{min}}J(y,u)& =\frac{1}{2}\parallel y-y_d{\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\frac{\alpha }{2}{\parallel u\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\beta {\parallel u\parallel }_{L^1\left(\mathsf{\Omega }\right)}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.-\mathsf{\Delta }y& =u+y_r\mathrm{in}\mathsf{\Omega },\hfill \\ \hfill y& =0\mathrm{on}\partial \mathsf{\Omega },\hfill \\ \hfill u& \in U_{ad}=\left\{v\left(x\right)\right|a\le v\left(x\right)\le b,\mathrm{a}.\mathrm{e}\mathrm{on}\mathsf{\Omega }\},\hfill \end{array}\right.$$

where $\mathsf{\Omega }={(0,1)}^2$, the parameters $\alpha =0.5$, $\beta =0.5$, $a=-0.5$, $b=0.5$. As this is a constructed problem, we set $y^{*}=sin\left(\pi x_1\right)sin\left(\pi x_2\right)$ and $p^{*}=2\beta sin\left(2\pi x_1\right)exp\left(0.5x_1\right)sin\left(4\pi x_2\right)$. Then, through Algorithm 4, we can obtain the optimal control solution $u^{*}={\mathsf{\Pi }}_{U_{ad}}\left(\frac{1}{\alpha }soft\left(p^{*},\beta \right)\right)$, the source term $y_r=S^{-1}{y}^{*}-u^{*}$ and the desired state $y_d=y^{*}-{\left(S^{*}\right)}^{-1}{p}^{*}$. Thus, we construct the example for which we know the exact solution.

We then test the mhADMM, the ihADMM and the classical ADMM for Example 1. The exact optimal control u and an example for the numerical optimal control obtained by mhADMM on the grid with $h=\frac{\sqrt{2}}{2^7}$ are shown in Figure 1.

In Table 1, we report the dimension of the control variable on each grid level, the error E of the control u, the EOC, the residual $\eta$, the computational time and the number of iterations obtained by the mhADMM, the ihADMM and the classical ADMM. As can be seen from the fourth and fifth column of Table 1, a high-subfigures are correct.precision solution does not improve the accuracy of the discretization error and the EOC. The computational time on the seventh, eighth and ninth columns show that the mhADMM is much faster than the ihADMM and the classical ADMM, especially when the discretization is at a fine level. The mhADMM algorithm can significantly reduce the computational cost and make the algorithm faster. This is mainly because the mhADMM adopts the strategy of gradually refining the grid, while the ihADMM and the classical ADMM compute the problem on a fixed grid size, which is their computational bottleneck. Moreover, the seventh column illustrates the mesh-independent performance of mhADMM; that is, the number of iteration of the mhADMM is independent of the discretization level. Above all, we can see that the mhADMM is much more efficient than the ihADMM and the classical ADMM.

Example 2. 

Consider

$$\left\{\begin{array}{cc}\hfill \underset{(y,u)\in H_0^1\left(\mathsf{\Omega }\right)\times L^2\left(\mathsf{\Omega }\right)}{\overset{}{min}}J(y,u)& =\frac{1}{2}\parallel y-y_d{\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\frac{\alpha }{2}{\parallel u\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\beta {\parallel u\parallel }_{L^1\left(\mathsf{\Omega }\right)}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.-\mathsf{\Delta }y& =u\mathrm{in}\mathsf{\Omega },\hfill \\ \hfill y& =0\mathrm{on}\partial \mathsf{\Omega },\hfill \\ \hfill u& \in U_{ad}=\left\{v\left(x\right)\right|a\le v\left(x\right)\le b,\mathrm{a}.\mathrm{e}\mathrm{on}\mathsf{\Omega }\},\hfill \end{array}\right.$$

where $\mathsf{\Omega }={(0,1)}^2$, the parameters $\alpha ={10}^{-4}$, $\beta ={10}^{-3}$, $a=-10$, $b=10$. The exact sparse solution of this problem is not known in advance. Instead, we use the numerical solutions computed on the grid with the grid size $h=\frac{\sqrt{2}}{2^9}$ as reference solutions.

As an example, Figure 2 presents the numerical optimal control for Example 2 on the grid with $h=\frac{\sqrt{2}}{2^7}$.

Table 2 shows the dimension of the control variable at each grid level, the error E of the control u, the EOC, the residual $\eta$, the computational time and the number of iterations for three methods. Similar to Example 1, the mhADMM still outperforms the ihADMM and the classical ADMM in terms of the computational time. The fourth and fifth column of Table 2 clearly show that a high-precision solution does not improve the accuracy of the discretization error and the EOC. Furthermore, the numerical results in the seventh column also illustrate the mesh-independent performance of our mhADMM algorithm. These results demonstrate that the mhADMM is more efficient than the ihADMM and the classical ADMM.

5. Conclusions

In this paper, we propose a new, efficient, multilevel, heterogeneous ADMM (mhADMM) algorithm for solving sparse elliptic PDE-constrained optimal control problems with $L^1$-control cost and box constraints on the control. Specifically, the inexact ADMM is first applied in the function space. Then, we propose the strategy of gradually refining the grid and employ the standard piecewise linear finite element to discretize the related subproblems appearing in each iteration of the inexact ADMM algorithm. Moreover, nodal quadrature formulas are utilized to approximately discretize the $L^1$-norm and $L^2$-norm to overcome the difficulty that the $L^1$-norm does not have a decoupled form. Finally, subproblems are solved by appropriate numerical methods. Theoretical results regarding the global convergence and iteration complexity are presented. In our numerical experiments, we show that the proposed mhADMM is superior to the ihADMM and the classical ADMM in terms of the efficiency.