Efficient numerical methods for elliptic optimal control problems with random coefficient

1.   Introduction

It is well known that many physical and engineering problems could be described by the optimal control problems (OCPs) with partial differential equation (PDE) or random/stochastic partial differential equation (RPDE/SPDE) constraints, such as the optimal heat source problem, the optimal design of the aircraft, and the weather forecasting problem etc. For the deterministic PDE optimal control problem, there exists plenty research results, and we refer to [15,21,26] for more details. However, only few works have been developed for the RPDE/SPDE optimal control problem.

Generally speaking, the numerical methods for RPDE/SPDE constraint optimal control problems are divided into two categories: sample methods and projection methods. For the former, Cao et al. proposed an improved Monte Carlo method using compression variance technique to deal with optimal control problems under random Burger constraints in 2003 [7]. Based on the gradient error, Kouri et al. combined the trust region method and the stochastic collocation method to propose an adaptive collocation method, which is superior to the traditional Newtonian conjugate gradient method [17]. For the latter, Xiu developed a polynomial chaos expansion for the target functional, which transforms the original SPDE constraint optimal control problem to a nonlinear problem with finite coefficients [25]. However, no theoretical guarantee was given there. In 2011, Gunzburger and his collaborators transformed the SPDE constraint optimal control problem into a KKT nonlinear system, and then proved the existence of the solution of the KKT system. By using Karhunen Loeve (KL) expansion or appropriate orthogonal decomposition, they obtained a series of deterministic optimization systems numerically [13]. In 2013, Kunoth and Schwab used generalized polynomial chaos expansion (gPC) to handle the RPDE/SPDE optimal control problems [18]. Recently, Taylor approximation and variance reduction method was proposed by Chen et al. for PDE-constrained optimal control under uncertainty [8].

In this paper, we mainly concentrate on the efficient numerical methods for the optimal control problem constrained by the Poisson equation with random coefficient. The major challenges for dealing with the randomness are: (Ⅰ) How to deal with the uncertainties appeared in the constraint and the objective functional, and formulate a discretized optimization problem. (Ⅱ) How to solve the resulted optimization problem, and propose an efficient algorithm. We will analyze the challenges above one by one, and present the corresponding solutions.

Monte Carlo (MC) method [5,7], the stochastic collocation (SC) method [1,2], the stochastic Galerkin (SG) method [22,23], the polynomial chaos (PC) expansion [16], and so on are some of the commonly used methods for dealing with the first issue. Based on the generality and parallelizability of MC method, we will apply it to deal with the random sample space. Further, using the finite element method (FEM) to approximate the physical space, we shall obtain a discretized optimization problem. This is the most natural idea, but it will face huge difficulties in computation and storage when the high accuracy solution is required. The KL expansion is a spectral representation of the covariance function, which consists of the eigenvalues and eigenfunctions. This method is widely used in RPDE/SPDE fields when the eigenvalues decay exponentially or non-Gaussian spatially-dependent random with large correlation length, otherwise it will run into the curse of dimensionality [3]. As an improvement, the multi-modes expansion (MME) [11] will be considered to approximate the state variable firstly, and then carry out MC and FEM subsequently. By exploiting the relationships among the terms of MME expansion, we can transform the discretized optimization problem into a quadratic programming problem with deterministic linear constraints. In this paper, we will follow this approach to resolve the first challenge.

The second challenge is closed related to the first one. If the first challenge can be handled very well, the latter one will be much less difficult. For solving the resulted optimization problem, there exist fruitful algorithms, such as the gradient descent method, the conjugate gradient method, the Newton's method, the regularization method, and the alternating direction method of multipliers (ADMM) etc., the interested readers are referred to [13,21] and references therein for the rich literature. Because of the separability of state and control of the resulted optimization problem, ADMM is applied in this paper, which has been used in the optimal control problem with deterministic PDE constraints, and displayed satisfactory results [14,26]. Moreover, to the best of our knowledge, there is no complete convergence analysis for elliptic RPDE-constrained optimization problems. The existing results merely present the difference between the exact solution and the discretized solution [13,19], while the error for solving the optimization problem is ignored. Here, based on the convergence analysis of ADMM in [6] and the techniques for dealing with the former challenge, the global error estimate in terms of the deviation between objective functionals will be established for our proposed algorithms later.

For the clarity of the structure, the efficient algorithm proposed in this paper is developed by three steps. Firstly, the MC, the FEM, and the ADMM are used sequentially to formulate a primitive algorithm, Algorithm 1, which needs to compute $M+1$ inverses of $N\times N$ matrices, and store $\mathcal{O}(M+1)$ matrices (see Section 3), where $M$ and $N$ denote the sample number in probability space and the degree of freedom in physical space, respectively. To reduce the computation cost, we propose algorithm 2, where the MME technique is introduced. Using this algorithm, we only need to compute inverse matrix twice, but still need to store $\mathcal{O}(M+1)$ matrices (see Section 4). As the last step, by applying the recurrence relation on MME, the random term in the constraint is shifted to the objective functional. Then, based on Algorithm 2, we propose Algorithm 3, in which we only need to compute inverse matrix twice, and store $\mathcal{O}(1)$ matrices (see Section 5). Naturally, Algorithm 3 reduces the computation cost and storage cost significantly.

The outline of this paper is as follows. In Section 2, the model problem, the well-posedness of the solution, and the first order optimal condition of the model problem are introduced. Algorithm 1, 2, 3, and their convergence analyses are presented in Section 3, 4, and 5, respectively. In Section 6, numerical simulations are carried out to test the performance of the proposed algorithms. Finally, some conclusion remarks are given in Section 7.

2.   Optimal control problem with random coefficient

The model problem, the well-posedness of the studied model, and its corresponding first order optimal condition will be presented in this section.

2.1.   The model in concern

For easy of description, most of the notations and definitions used in this paper are adopted from the reference [9]. We consider the model problem on a polygon $D\subset \mathbb{R}^{2}$ with the boundary $\partial D$. $L^{p}(D)$, $L^{p}(\partial D)$, $H^{p}(D)$, and $H^{p}(\partial D)$ denote the Lebesgue integrable function spaces and Hilbert spaces on $D$ and $\partial D$, respectively ($1\leq p\leq+\infty$). In addition, we introduce a standard probability space $(\Theta, \mathcal{F}, \mathcal{P})$, and a Bochner space $L^{p}(\Theta, X)$ with the norms

where $X$ stands for a Hilbert space.

Using the notions above, the Poisson equation constraint optimal control problem with random coefficient and Dirichlet boundary condition could be written as

with

where $y_{d}\in L^{2}(D)$ and $u_{d}\in H^{2}(D)$ stand for the desirable state variable and the control variable respectively, which are given in advance. The random diffusion coefficient $\alpha(\xi, x)$ belongs to $L^{2}(\Theta, D)$, where $\xi$ represents a random sample in $\Theta$. The expectation of $\rho(\xi)$ is defined by $\mathbb{E}[\rho(\xi)] = \int_{\Theta}\rho(\xi)\, \mathrm{d}\mathcal{P}(\xi)$, and the boundary condition $g\in L^{2}(\partial D)$ is a known function, which is set to be $0$ in this paper for simplicity.

The corresponding variational formulation of the Poisson equation with random coefficient is given by: Finding $y\in L^{2}(\Theta, H_{0}^{1}(D))$, such that

with $a(y, v) = (\alpha (\xi, x)\nabla y, \nabla v)_{L^{2}(D)}$. In order to ensure the well-posedness of the variational formulation (2), the following assumptions on $\alpha$ are needed.

Assumption 1. ([3]) The random diffusion coefficient $\alpha(\xi, x)$ in $(OCP)$ is assumed to be a strongly measurable mapping from $\Theta$ to $L^{\infty}(D)$. And there exist constants $0<\alpha_{1}<\alpha_{2}<\infty$ such that random coefficient $\alpha(x, \xi)$ is uniformly elliptic, i.e., for every $\xi\in \Theta$, there holds

Lemma 2.1. ([2,4,12]) If Assumption 1 holds, then for any given right-hand-side function $u\in L^{2}(D)$, the variational equation (2) has a unique solution $y(\xi, x)\in L^{2}(\Theta, H_{0}^{1}(D))$. Furthermore, we have the following pointwise $H^{1}$-norm estimate

and the global $L^{2}$-error estimate

where $C(\alpha)$ is a positive constant depending on $\alpha_{1}$ and $\alpha_{2}$.

2.2.   First order optimal condition

For simplicity, the diffusion coefficient $\alpha(\xi, x)$ in (OCP) is set to be

other cases could be considered similarly. Here, $\sigma(\xi, x)\in L^{2}(\Theta, L^{\infty}(D))$ denotes the random process satisfying

and $\varepsilon$ stands for the magnitude of the perturbation. It is not difficult to obtain that the random coefficient $\alpha$ defined in (5) satisfies Assumption 1. Therefore, according to Lemma 2.1, there exists a linear operator $T: L^2(D)\rightarrow L^2(\Theta, H^1(D))$, such that

where $u$ and $y$ satisfy the variational equation (2). Further, it follows from the inequality (3) that the following estimate holds,

where $j = 2$ or $4$, and $C_{1}$ is a constant depending on the bound of $\alpha$.

Therefore, if the state variable space is $Y = L^{2}(\Theta, H^{1}_{0}(D))$, and the control variable space is $U = L^{2}(D)$, the optimal control problem (OCP) could be rewritten as

which can also be reformulated as

Here, $J$ is a functional from $L^{2}(D)$ to $\mathbb{R}$, and the existence and uniqueness results of the solution for (8) could be found in [16,19].

Finally, let $u^{*}$ be the optimal solution of (8). By applying the definition of the Fr$\acute{e}$chet derivative [9], the linearity of the expectation operator $\mathbb{E}$, and the solution operator $T$ defined in (6), the first order optimal condition of (P) can be written as

where $T^{*}$ is the adjoint operator of $T$.

Now, the optimal control problem (P), and its first order optimal condition have all been presented. In the remainder of the paper, we will propose three algorithms for solving the optimal control problem (P) incrementally, and compare them with each other.

3.   FEM-MC-ADMM for the model problem

The first algorithm studied in this paper is presented in this section. The motivation comes from the finite element method (FEM), the Monte Carlo approximation (MC), and the alternating directions method of multipliers (ADMM). We will carry out these techniques sequentially to solve the optimal control problem (P), and present an global error estimate in the form of the deviation between the objective functionals.

3.1.   FEM discretization

Let $\mathcal{K}_{h}$ be a regular triangulation of the spatial domain $D$, then the linear finite element function space could be defined by

where $\mathbb{P}_{1}$ is the space of polynomials of degree less than or equal to $1$. Further, let ${\bf R} = {\bf R}(x) = (\phi_1(x), \cdots, \phi_N(x))$ denotes a vector valued function consisting of all the basis functions of $V_{h}$, the state variable $y(\xi, x)\in L^{2}(\Theta, H^{1}(D))$ can be approximated by

which belongs to $Y_{h} = L^{2}(\Theta, V_{h})$. Here, the vector ${\bf y}(\xi) = (y_1(\xi), \cdots, y_N(\xi))^{T} \in \mathbb{R}^{N}$ for any fixed $\xi \in \Theta$.

Based on above notations, for a fixed sample $\xi\in \Theta$ and a given control variable $u\in L^{2}(D)$, the FEM approximation for the constraint equation of (OCP) shall be described as : Finding $y_h(\xi, \cdot)\in V_h$ such that

This means there exists an operator $T_{h}:L^{2}(D)\rightarrow L^{2}(\Theta, V_{h})\subset L^{2}(\Theta, H^{1}(D))$ such that $y_h = T_{h}(u)$. It is easy to see that $T_{h}$ is a linear operator.

Lemma 3.1. The pathwise FEM discretization (11) satisfies the following stability estimate: for any fixed sample $\xi\in \Theta$,

where $c(\alpha)$ depends on $\alpha$.

The proof, which will be omitted here, follows from the coerciveness and continuity of the bilinear form $a(u, v)$. The interested readers are referred to [15].

Now, we generalize the stability estimate of the deterministic problem to the random coefficient case. It is well known that the FEM approximation of (2) is: Finding $y_h\in L^{2}(\Theta, V_h)$ such that

By applying the Assumption 1 and the estimate (12), we can derive:

Lemma 3.2. For any control variable $u\in L^{2}(D)$, if $y_{h}$ is the FEM solution of (13), then we have

where $j = 2, 4$, and $C_{2}$ is a positive constant depending on $\alpha$ and $j$.

By using the estimates (4), (7), (14), the Theorem 4.3 in [10], and the triangle inequality, we can obtain the following lemma.

Lemma 3.3. ([11,19]) If the control variable $u\in L^{2}(D)$, then it holds

where $j = 2, 4$, and $C_{3}$ is a positive constant independent of the mesh size $h$.

Similar to (10), the control variable $u$ could be approximated by

where ${\bf u} = (u_1, \cdots, u_N)^{T} \in \mathbb{R}^{N}$. Substituting the state variable $y$ and the control variable $u$ in (P) by $y_h$ and $u_h$, respectively, we can get a stochastic optimization problem

3.2.   MC approximation

In this subsection, we mainly apply the standard MC method to approximate the expectation $\mathbb{E}$ in (P$_{1h}$). Let ${\{\xi_{j}\}_{j = 1}^{M}}$ ($M \gg 1$) be the identically distributed independent samples chosen from the probability space $(\Theta, \mathcal{F}, P)$. For any state variable $y\in L^{2}(\Theta, H^{1}(D))$, by using the law of large numbers and the central limit theorem, we can approximate the expectation $\mathbb{E}[y_h(\xi, x)]$ by

which satisfies the error estimate (cf. [4])

Substituting the MC approximation (17) into (P$_{1h}$), we can derive the following discretized optimization problem:

where

Similar to (P), the optimization problem (P$^M_{1h}$) has a unique solution $u_{M, h}^{*}\in V_{h}$, which satisfies the first order optimal condition

where $T_h^{*}$ is the adjoint operator of $T_h$, and $\mathbb{Q}_{h}$ is the $L^{2}$ projection operator from $L^2(D)$ to $V_{h}$.

Using the first order optimal conditions (9) and (20), we can get

Lemma 3.4. ([19]) Let $u^*$ and $u_{M, h}^*$ be the optimal controls of (P) and (P$^M_{1h}$), respectively, we have the following estimate

with

where the constant $C$ is independent of $M$ and $h$.

Further, applying the definitions of the functionals $F(y, u)$ and $F_{M, h}(y_{h}, u_{h})$, the operators $T$ and $T_h$, and the error estimate (21) in Lemma 3.4, we can get

Theorem 3.5. Let $(y^{*}, u^{*})$ and $(y_{M, h}^{*}, u_{M, h}^{*})$ be the optimal solutions of the optimal control problem (P) and it approximation (P$^M_{1h}$), respectively, then we have

where

For the sake of structural consideration, we put the proof in Appendix.

3.3.   ADMM algorithm for (P$_{1h}^M$)

In this subsection, we will introduce the ADMM to solve the optimization problem (P$^M_{1h}$), and present the error estimate of this algorithm. For these purposes, we need to rewrite the optimization problem (P$^M_{1h}$) in the matrix-vector form firstly. Replacing $y_h$ and $u_h$ in (19) by ${\bf R}{\bf y}(\xi)$ and ${\bf R}{\bf u}$, which have been defined in (10) and (16), the optimization problem (P$^M_{1h}$) can be rewritten as

where $\mathcal{F}_{M, h}({\bf y}(\xi), {\bf u}) = F_{M, h}({\bf{R}\bf y(\xi)}, {\bf{R}\bf u}) = F_{M, h}(y_{h}, u_{h})$, and

The augmented Lagrangian functional of the optimization problem $(\widetilde{\mathrm{P}}^M_{1h})$ is

where ${\bf \lambda}(\xi)\in \mathbb{R}^{N}$ is the Lagrange multipliers and $\beta$ denotes the penalty parameter. Here and hereafter, the notation $(\cdot, \cdot)$ represents the inner product of vectors, and $\|\cdot\|$ stands for the $l^{2}$ norm of a vector. It is clear that the optimization problem ($\widetilde{\mathrm{P}}^M_{1h}$) equals the saddle-point problem (cf. [21])

For the saddle-point problem (23), the ADMM is an efficient algorithm. We will apply ADMM to solve this problem, and the whole process can be described by

Notice that every iteration of the ADMM has the closed-form solution

with $j = 1, \dots, M$.

Now, we consider the convergence of the ADMM. It is easy to find that the stiffness matrixs $A_{1}(\xi_j)$, $j = 1, \dots, M$, and the mass matrix $B_{1}$ are all symmetric positive definite matrices, and $\mathcal{F}_{M, h}({\bf y}(\xi), {\bf u})$ is a convex functional. For the deterministic case under conditions above, the convergence result has been established in [6]. Applying the Assumption 1, it is not hard to obtain the following convergence result for the random case (cf. [20,24]).

Lemma 3.6. For the problem (23), after $t$ iterations, the iterate solution of the objective function generated by ADMM converges at the rate $\mathcal{O}(\frac{1}{t})$, that is

where $({\bf y}^{*}_{M}, \; {\bf u}^{*}_{M})$ is the theoretical solution of (23), and

Further, applying the Theorem 3.5, the identity

the Lemma 3.6, and the triangle inequality, we can get the global error estimate.

Theorem 3.7. A global error estimate for Algorithm 1 is as follows:

From Algorithm 1, we can see that the stiffness matrix $A_{1}(\xi)$ is a random one, which needs to be computed and stored for every sample $\xi_j$, $j = 1, \dots, M$. The cost to compute $A_{1}(\xi)$ and the mass matrix $B_1$ is $\mathcal{O}((M+1)N^{2})$. In addition, from the closed form solutions (24), we also need to calculate and save the inverse matric $\left(\gamma B_{1}+\beta B_{1}^TB_{1}\right)^{-1}$, and the inverse matrices $\left(B_{1}+\beta A_{1}(\xi_j)^T A_{1}(\xi_j)\right)^{-1}$ for all the samples. Thus, we must compute $M+1$ inverses of $N\times N$ matrices, which will cost a large amount of memory and computation time.

4.   MME-FEM-MC-ADMM for the model problem

To reduce the computational cost in Algorithm 1, we shall introduce the multi-modes expansion (MME), and propose a new algorithm. The main idea is to approximate the state variable by MME technique, and then transform the original constraint in Algorithm 1 and thus shift the randomness to the right hand side mass matrix. The advantage is that we can calculate the expectation of the right hand side random matrix in advance, then compute and store the inverse of the matrix only once, which could reduce the computational cost significantly.

4.1.   Multi-modes expansion

The multi-modes expansion for the state variable $y$ can be expressed as (cf. [11])

where $\varepsilon$ is the magnitude of the perturbation defined in (5), and $y_{q}\in L^{2}(\Theta, H_0^{1}(D))$. Substituting (27) into the state equation of the optimal control problem (OCP), and comparing the coefficients of the powers of $\varepsilon$, we could get the following Poisson equations in recursive form.

It is apparent that the relationship between $y_{0}$ and $u$ is deterministic, so we can solve the equations (28) by recursion for $q\geq 1$. The well-posedness of the MME solutions has been established in Theorem 3.1 of the reference [11].

Lemma 4.1. The equations (28) has a unique solution $y_{q}\in L^{2}(\Theta, H_{0}^{1}(D))$. Besides, for any $q\geq 0$, if $y_{q}\in L^{2}(\Theta, H^{2}(D))$, then the following energy estimate holds

where $c(\alpha)$ is given by (12) and is independent of $q$ and $\varepsilon$.

For the purpose of numerical solution, we truncate the expression $(27)$ after the first $Q$ terms, and denote it by

Lemma 4.2. ([11]) If the perturbation magnitude $\varepsilon <\min{\{1, c({\alpha})^{-\frac{1}{2}}\}}$, then we have

Lemma 4.1 shows that the problem (28) is well-posed, and Lemma 4.2 shows that the expression (30) is an approximation of the state variable $y$. Based on these facts, there must exist a linear operator $T^{Q}:L^{2}(D)\rightarrow L^{2}(\Theta, H^{1}(D))$ such that

Feng et al. [11] showed that the truncation error of MME after $Q$ terms is $\mathcal{O}(\varepsilon ^{Q})$.

Lemma 4.3. Under the assumptions of Lemma 4.1 and Lemma 4.2, it holds that

where $j = 2$ or $4$, $C_{4}$ is a positive constant depending on $\alpha$, $Q$ and $j$.

4.2.   FEM and MC approximation

In this subsection, we will apply FEM to discretize $y_{q} (q\geq 0)$ on physical space, and MC method to approximate the expectation $\mathbb{E}$ on probability space. Then the discretized optimization problem, its first order optimal condition, and the convergence analysis for the proposed approach will be presented sequentially.

For any fixed $\xi\in\Theta$, the FEM approximation of $y_{q}(\xi, x)$ in (28) can be written as

where ${\bf y_q}(\xi)\in \mathbb{R}^N$. Naturally, the truncated MME state $y^{Q}(\xi, x)$ could be approximated by

And the variational formulation of the constraint equation in the truncated MME form is given by: For each $\xi\in \Theta$, find $y^Q_h(\xi, x)\in L^2(\Theta, V_h)$ such that

for any $v_h(x)\in V_h$, where $q = 1, \dots, Q-1$. Similar to (31), there exists a linear operator $T_{h}^{Q}:L^{2}(D)\rightarrow L^{2}(\Theta, V_{h})\subset L^{2}(\Theta, H^{1}(D))$ such that

Lemma 4.4. (cf. [11]) Under the assumption of Lemma 4.2, for $j = 2, 4$ we have the following stability estimate

and the error estimate

where $C_5$ and $C_{6}$ are constants depending on $\varepsilon$ and $Q$.

Further, applying the triangle inequality, the estimates (32) and (38), we have the following error estimate.

Lemma 4.5. Under the conditions of Lemma 4.1 and Lemma 4.4, we have the following estimate

where $j = 2$ or $4$.

By arguments similar to the ones used to get the optimization problem (19), applying the expressions (16), (34) and (36), the original problem (P) could be approximated by

where ${F}_{M, h}$ has been given in (19). Similar to (P$_{1h}^{M}$), the optimization problem (P$_{2h}^{Q}$) has a unique solution $u_{M, h}^{Q*}\in V_{h}$, which satisfies the first order optimal condition

Using the first order optimal conditions (9) and (40), we have

Lemma 4.6. (cf. [19]) Under the conditions of Lemma 4.1-4.5, the following error estimate

holds, where

Similar to Theorem 3.5, we can get the error estimate on the difference of the objective functionals, as follows

Theorem 4.7. Let $(y^{*}, u^{*})$ and $(y_{M, h}^{Q*}, u_{M, h}^{Q*})$ be the optimal solutions of the optimal control problem (P) and its approximation (P$^Q_{2h}$), respectively, then we have

where

The proof is similar to Theorem 3.5, we omit it here.

4.3.   ADMM algorithm for (P$_{2h}^{Q}$)

In this subsection, the ADMM algorithm for solving the optimization problem (P$_{2h}^{Q}$) will be presented. In order to describe the algorithm clearly, we simplify (P$_{2h}^{Q}$) firstly.

Substituting the control variable $u$ in the variational formulation (35) by $u_h$ defined in (16), and using the definition (33) on $y_{q, h}$, the constraint of the optimization problem (P$_{2h}^{Q}$) could be rewritten as

with

Applying the definition of ${\bf y}^Q_{h}$ in (34), we know that

Further, combining above equality with the equations (43), the relation between state variable ${\bf y}^Q(\xi)$ and control variable ${\bf u}$ can be given by

where

Therefore, the optimization problem (P$_{2h}^{Q}$) can be reformulated as:

where the functional $\mathcal{F}_{M, h}$ is defined in (23).

The optimization problem ($\widetilde{\mathrm{P}}^Q_{2h}$) is equivalent to the saddle point problem

with the augmented Lagrangian function

where ${\bf \lambda}(\xi)\in \mathbb{R}^{N}$ is the Lagrange multipliers and $\beta$ denotes the penalty parameter.

Now, the ADMM algorithm for solving the saddle point problem (44) could be described by

The ADMM process has the following closed form solutions

with $j = 1, 2\dots, M$. Similar to Lemma 3.6, let $({\bf y}^{Q, *}_M, {\bf u}^{Q, *}_M)$ be the exact solution of the problem $(\widetilde{P}^Q_{2h})$, which satisfies $(y_{M, h}^{Q*}, u_{M, h}^{Q*}) = ({\bf{R}}{\bf y}^{Q, *}_M, {\bf{R}}{\bf u}^{Q, *}_M)$, then there exists a constant $\widetilde{K}_3$ such that

Furthermore, applying the estimate of Theorem 4.7 and triangle inequality, we obtain the following estimate:

Theorem 4.8. Let $({\bf y}^{Q, t}, {\bf u}^{t})$ denote the solution obtained after using FEM, MME, and MC approximation and $(y^{*}, u^{*})$ be the solution of (1). Then we have

where $F(\cdot, \cdot)$ and $F_{M, h}(\cdot, \cdot)$ are functionals defined by (1) and (19), respectively.

From the iterations (45), we can see that all the matrices are deterministic except the matrix $B_{2}(\xi)$. Fortunately, the matrix $\mathbb{E}_M [B_2(\xi)^T B_2(\xi)]$ could be obtained in advance for the given samples. Therefore, we need to compute only two inverse matrices, which could reduces the computational cost on the inverse matrices in the Algorithm 1 greatly. However, in addition to the cost of computing and storing $A_2$, we still need to compute and store the mass matrices $B_{2}(\xi_j)$ for all samples. Thus, the storage requirement does not decrease, it is still $\mathcal{O}((M+1)N^{2})$.

5.   SMME-FEM-MC-ADMM for the model problem

In this section, a new algorithm is proposed to overcome the disadvantages of the Algorithm 2, which can reduce the storage requirement significantly.

From the variational formulation (35), we know that $y_{q}$ depends on the preceeding term $y_{q-1}$ for $q\geq 1$. Therefore, we can establish a relationship between $y_0$ and $y_q$. Furthermore, applying the definition (30) of $y^{Q}$, there must exist an relation between $y_0$ and $y^Q$. Based on above facts, using the equations (34) and (43), we can obtain

where $P(\xi)$ is a matrix with random variables, defined as follows

The relationship (47) implies that ${\bf y}^Q(\xi)$ is determined by ${\bf y}_0$, so we need to consider only the constraint on ${\bf y}_0$. The optimization problem $(\widetilde{\mathrm{P}}^Q_{2h})$ shall be simplified as

where the functional $\mathcal{F}_{M, h}$ is defined in (23), $A_{3} = A_{2}$ and $B_{3} = B_{1}$ are two deterministic matrices defined in (43). This problem is a classical optimization problem, which is the same as

with the augmented Lagrangian function

where ${\bf \lambda}\in \mathbb{R}^N$ is Lagrange multiplier independent of the samples, and $\beta$ is the penalty parameter.

Now, we present the ADMM algorithm to solve the saddle point problem (48).

By straightforward calculations, we can get

It is not hard to see that Algorithm 3 has the same convergence property as Algorithm 2 (see Theorem 4.8). From the iterations (49), we can see that we don't need to store the information for all the samples. The only matrices needing to be remembered are $\mathbb{E}_{M}\big[P(\xi)\big]$, $A_3$, and $B_3$, which means the memory cost is $\mathcal{O}(N^{2})$. Moreover, we need to compute inverse matrices only twice. The above advantages are in sharp contrast with the Algorithm 1 and Algorithm 2. Hence, Algorithm 3 is the most efficient one both in computational and storage aspects.

6.   Numerical experiments

In this section, results of numerical experiments, which are carried out to test the performance of the proposed algorithms, will be reported. The convergence order of the FEM for all algorithms is verified firstly. And then, the ADMM convergence order in the form of the objective function is shown. Finally, to illustrate the advantage of Algorithm 3, we present the storage costs of the three algorithms. All the simulations are conducted using Matlab on a high-performance computer with $16$ kernels and $256$ GB RAM.

We consider the optimal control problem on the spatial domain $D = [-1, 1]^{2}$, and assume the random variable $\sigma(x, \cdot)$ associated with the random coefficient $\alpha(\xi, x)$ obeys the uniform distribution on the interval of $[-1, 1]$ for any fixed $x\in D$. The weighted parameter $\gamma$ of the objective functional is set to be $1$, the desired state and control are set to be

It is not difficult to find that $y_d$ and $u_d$ are the optimal state and control of the original problem (P) in the expected sense. In addition, the parameters associated with the ADMM are set to be $\beta = 1$, $\varepsilon^0_{i} = 10^{-4}$, $i = 1, 2, 3$, respectively.

To verify the convergence order of the FEM, we fix other parameters firstly. For simplicity, the sample number of the MC, the iteration steps of the ADMM, the magnitude of the perturbation defined in (5), and the truncated length of the MME are set to be $M = 1000$, $t = 20$, $\varepsilon = 10^{-4}$, and $Q = 1$, respectively. Other cases could be discussed in similar manners. Now, we carry out our algorithms on the nested meshes with $h = (1/2)^{i}, \; i = 1, 2, 3, 4, 5$, and present the results in Figure 1. We emphasize that Algorithm 1 is out of memory even with spatial meshsize $h = 1/32$. For aesthetics, we extend the line to the fifth point. From Figure 1, we can see that the FEM convergence order with respect to the control variable $u$ and the state variable $y$ are all of order $2$ for the proposed three algorithms, which coincide with the theoretical analysis in Sections 3-5. To illustrate the correctness of the algorithms intuitively, the images of the optimal control in theory, and the numerical solutions on the proposed algorithms with $h = 1/32$ are shown in Figure 2, where the numerical solution for Algorithm 1 is obtained with $h = 1/16$ for the memory limit.

Now, we will test the performance of ADMM. For this purpose, other parameters in the proposed algorithms need to be fixed in advance. Without loss of generality, we set $M = 1000$, $h = 1/32$, and $Q = 1$ for the latter two algorithms and the Algorithm 1 with $h = 1/16$. The ADMM errors for the proposed algorithms in the form of the deviation between functionals are shown in Figure 3. From the results, we can see that the convergence rate is superior to $\mathcal{O}(\frac{1}{t})$, which means that the numerical performance of ADMM is better than the theoretical case.

Finally, we analyze the CPU memory cost of the proposed algorithms, which is one of the main concerns in this paper. For the fixed MC parameter $M = 1000$, the MME perturbation magnitude $\varepsilon = 10^{-4}$ and the iteration number $t = 20$, the variation trends of the CPU memory cost with respect to the FEM parameter $h$ and the MME parameter $Q$ are given in Table 1. From Table 1, we can find three obvious facts: (1) For any one of the three algorithms, the memory cost increase as the parameter $h$ decreases. (2) For any one of the Algorithm 2 and 3, the memory cost does not have obvious change with respect to the parameter $Q$. (3) For any fixed $h$ and $Q$, the memory cost of Algorithm 3 is significantly lower than other two algorithms. The above facts agree with the theoretical analysis in Section 3, 4, and 5 very well.

In order to observe the storage changes of the three proposed algorithms with respect to the sample number $M$, we fix the FEM parameter $h = 1/16$, the stochastic perturbation parameter $\varepsilon = 10^{-4}$ and the ADMM iteration number $t = 20$. The CPU memory costs with different cases are presented in Table 2. From Table 2, we can see that (1) The memory costs of Algorithm 1 and 2 increase linearly with the sample number $M$, which roughly equals $\mathcal{O}((M+1)N^{2})$. (2) The storage cost of Algorithm 3 basically remain unchanged with the increasing of the sample number, that is approximately $\mathcal{O}(N^{2})$. Based on all the results exhibited in Table 1 and Table 2, we can conclude that Algorithm 3 is an efficient method for the optimal control problem with random coefficient.

7.   Conclusions

In this paper, three algorithms are proposed for the elliptic equation constraint optimal control problem with random coefficient. The Algorithm 1 is based on the FEM discretization, the MC approximation, and the ADMM iteration, which has global convergence. However, the computational cost and the storage cost are both very high. As an improvement, the Algorithm 2 is proposed, which reduce the computational cost by using the MME technique, but memory cost remain high. To overcome the disadvantage of Algorithm 2, Algorithm 3 is proposed, which reduce the memory cost significantly by transforming the constraint with random term to a deterministic one. Numerical simulations are carried out to verify the theoretical analysis, which illustrate the efficiency of the Algorithm 3 further. The applications of the proposed algorithms to other random optimal control problems will be discussed in our future works.

Acknowledgments

The work of H. Song was supported by the NSF of China under the grant No.11701210, the science and technology department project of Jilin Province under the grant No. 20190103029JH and 20200201269JC. the education department project of Jilin Province under the grant No. JJKH20180113KJ, and the fundamental research funds for the Central Universities. The work of X. Pang was supported by the NSF of China under the grant No. 11871245, 11771179, and 11726102. The work of J. Zhang was supported by the NSF of China under the grant No. 11971221 and 11731006, the Shenzhen Sci-Tech Fund under the grant No. JCYJ20170818153840322 and JCYJ20190809150413261, and Guangdong Provincial Key Laboratory of Computational Science and Material Design under the grant No. 2019B030301001. The authors also wish to thank the High Performance Computing Center of Jilin University, the Computing Center of Jilin Province, and the Key Laboratory of Symbolic Computation and Knowledge Engineering of Ministry of Education for essential computing support.

Appendix. The proof of Theorem 7.

Proof. Using the definitions of the $F(y, u)$ and $F_{M, h}(y_{h}, u_{h})$, and subtracting (19) from (1), we can get

For the first term $i$, we have

By using the MC error estimate (18), the triangle inequality, the estimate (7) with $j = 4$, and the inequality (52), it yields

Applying the Cauchy-Schwarz inequality, the triangle inequality, the MC error estimate (18), and the inequality

we can get

and

Taking $j = 2, 4$ in (7), (14), and (15), it follows from (52) that

Similarly, we obtain

When $M\gg 1$ and $h\ll 1$, we have

Next, we deal with the second term, applying the Cauchy-Schwarz inequality, triangle inequality, and the estimate (21), we obtain

Substituting $I_{1}$, $I_{2}$, and $II$ into (50), the conclusion is confirmed naturally.