A multi-mode expansion method for boundary optimal control problems constrained by random Poisson equations

1.   Introduction

In recent years, optimization problems governed by PDEs with uncertain coefficients have been the subject of growing interest in the scientific community [11,3,22]. This subject lies at the interface of PDEs with uncertain coefficients, optimization in Banach spaces, and stochastic programming.

In this paper, we study an optimal control problem constrained by a Poisson equation with uncertain diffusion coefficients. The control is a determinate function in the Dirichlet boundary condition. Our goal is to design efficient algorithms for this stochastic boundary control problem.

Deterministic optimal control problems constrained by partial differential equations have been well developed and investigated for several decades. While, there are not as much papers devoted to the random optimal control problems. In contrary to the deterministic problems, it is more expensive to get the numerical solution for the random optimization problem. Since the parameters are uncertain, the resulting states are random fields, i.e., random variables. And the inclusion of the stochastic dimension brings in additional freedom in the cost functionals. Therefore, we need to discretize the governing PDE in space and approximate the variables in random field. Moreover, it needs special probability theories to analyze and handle the stochastic domain in the optimal control problem. Existent algorithms for this problem always need to solve a large number of deterministic PDEs at each optimization iteration, and the memory consumption is huge. For the random discretization, commonly used methods are mainly divided into two categories: projection-based methods, such as stochastic Galerkin (SG) method and generalized polynomial chaos (gPCs) method, etc.; and sample-based methods, such as Monte Carlo (MC) method and stochastic collocation (SC) method, etc.. Naseri and Malek apply the SG method combined with preconditioned Newton's conjugate gradient method to solve an SPDE-constrained optimization problem with random force [28]. Cao designs an efficient MC method with variance reduction technique to solve random Burgers equation-constrained optimization problem [9]. In [30], the authors apply the SC method together with a gradient descent for the SPDE constrained control problem. Kouri and his collaborators improve the efficiency of SC method by adopting adaptive sparse-grid collocation with a trust-region framework for the random discretization [23]. Currently, a domain decomposition techniques have been applied for the random optimal control problems in [21]. In [5], the authors presented a low-rank tensor method to discrete PDE-constrained optimization problem. And we refer to [1,24,25] for some theoretical results.

The SG method provides a solid mathematical framework for the analysis and algorithms design, but it is not always the most computationally efficient means of solving large problems. While, the stochastic collocation method is often more popular than the projection-based method for it only needs to solve a collection of decoupled deterministic problems instead of a stochastic problem. However, the collocation method usually be exploited to a large class of random PDE-constrained optimization problems at the cost of losing the non-intrusivity property(see [4,32]). Over the past several year, a new stochastic discretization method has been proposed by Feng and his colleges for the Stochastic PDEs, which shows higher efficiency than traditional methods, see [13,15,16,14]. Inspired by this idea, we developed a multi-mode representation preprocessing technique for the optimal control problems constrained by random Helmholtz equations in [26]. This technique combines the merits of projection-based methods and sample-based methods, which reduces the number of equations evaluations sharply. Meanwhile, it is robust and easy to implement. The computational complexity is similar to solve a deterministic PDE-constrained control problem and a mathematical expectation of random matrices.

The main focus of this paper is on the efficient discretize then optimize numerical method for solving the optimal control problem constrained by random Poisson equation with deterministic control. While most of existing articles for the random PDEs constrained control problem are only concerned with distributed parameter controls, it is not a realistic situation since such controls are not easy to implement in practice. We think that boundary controls are much more natural and only consider the model problems with these controls. In order to overcome the difficulties of unbearable computation and memory consumption, we first do some pre-processing for the model problem with a modified multi-mode expansion (MME) form. The key idea is to transfer the randomness from the coefficients in the state equations to the coefficients in the cost functional, in which the random field is much easy to handle. With the MME technique, we can convert the random state equation to a deterministic one. Then discrete the physical space by the finite element method (FEM). After approximating the expectation in the cost functional with Monte Carlo method, we get a deterministic PDE-constrained control problem which can be solved by an alternating direction method of multipliers (ADMM) fastly. The memory cost of the proposed algorithm is only $O(N^2)$. The theoretical analysis in this paper are more difficult than those in the random Helmholtz equations constrained control problems for the randomness of diffusion coefficients in the state equations. We prove the existence of optimal control and deduce the global error estimate of the algorithm.

The rest of the paper is organized as follows. In Section 2 we present an optimality problem constrained by random Poisson equations under unconstrained control (UC) condition and box-constrained control (BCC) condition. The existence of optimal solutions and the first order necessary conditions are also discussed in this part. In Section 3, we first deduce the multi-modes representation approximation scheme for the model problem, then discrete the physical space by FEM. After that, we transform the new random discrete optimal control problem into a deterministic optimization problem. This deterministic problem will be solved by ADMM and two-block ADMM for the UC and BCC cases respectively in section 4. Section 5 presents the global error estimates for two algorithms. And Section 6 provides numerical experiments which demonstrate the efficiency of our algorithms. The last section is devoted to some concluding remarks.

2.   Model problem

2.1.   Definitions

Let $D\subset \mathbb{R}^d(d = 1,2,3)$ be a bounded Lipschitz domain. Suppose $H^{m}(D) = \{v\in L^2(D)\; \big|\; \partial^\alpha v\in L^2(D)\; ,\; 0\leq |\alpha| \leq m\},\quad m\in \mathbb{N}^+$ is the Sobolev space on $D$ with the norm

Let $(\Omega,\mathcal{F},P)$ be a complete probability space, where $\Omega$ is the set of outcomes, $\mathcal{F}$ is the $\sigma$-algebra of events and $P$ is the probability measure. Let the Bochner space $L^p(V,\Omega)$ be a space of strongly measurable functions $w$ ranging in a separable Hilbert space $V$ with the norm

Before listing the model problem, a detailed description of the state equations and some necessary properties will be given. Suppose that the deterministic control variable $u(x)$ and the random state variable $y(x,\xi)$ satisfies the following random elliptic equations,

In the elliptic equations, $a(x,\xi)$ is the diffusion coefficient with a small random perturbation defined by:

Here, the random process $\eta\in L^2( L^\infty(D);\Omega)$ satisfies

and the parameter $\varepsilon$ represents the magnitude of the random perturbation.

To ensure the well-posedness of the forthcoming model problem, it requires the following assumption on the random diffusion coefficient $a$:

Assumption 2.1. There exist constants $0<a_-<a_+<\infty$ such that the random coefficient $a(x,\xi)$ satisfies:

The corresponding weak formulation of the random elliptic equation (1)$\sim$ (2) is given by: Find $y\in L^2(H_u^1(D);\Omega)$, such that

with $H_u^1(D): = \{y\in H^1(D)| y = u \; \text{on}\; \partial D\}$, $H_0^1(D): = \{y\in H^1(D)| y = 0 \; \text{on}\; \partial D\}$, $(v_1,v_2)_D = \int_{D}v_1(x,\xi)v_2(x,\xi)\,\mathrm{d}x$ and $\mathbb{E}[w(\xi)] = \int_{\Omega}w(\xi)\,\mathrm{d}P(\xi)$.

Lemma 2.2. Under the assumption 2.1, for every $u\in L^2(\partial D)$, the elliptic equaiton (1) admits a unique solution $y\in L^2(H_u^1(D);\Omega)$. Further more, if $u\in H^{1/2+\sigma}(\partial D)$, then $y\in L^2(H_u^{1+\sigma}(D);\Omega), 0<\sigma \leq 1$, and from the elliptic regularity theory it has the following regularity estimation,

where $C(D,a)$ is some constant dependent on domain $D$ and coefficient $a$.

Proof. The existence of the solution to the equation (5) can be deduced from the Lax-Milgram Theorem. And the regularity estimate (6) follows by the elliptic regularity theory together with the transposition method (see [27]).

Note $Y: = L^2(H^{1+\sigma}(D),\Omega)$ and $U: = H^{1/2+\sigma}(\partial D),0<\sigma \leq 1$ as the state space and control space respectively. Lemma 2.2 indicates that for all $u \in U$, there is exactly one random state $y = y(u)\in Y$. So we use the continuous mapping $S$ in the following to denote the "random weak solution operator" of equation (5), i.e.

From (6), the estimate of the operator $S$ is

In next part, the target random optimal control problem will be represented in the form of the "random weak solution operator". By $S$, we are able to formally do the FEM error analysis for the model problem.

2.2.   Random elliptic equation constrained control problem

With the definition of $S$ in (7), the optimal control problem with random elliptic constraint is given by

where

in which $y_d\in Y$ and $u_0\in U$ are given functions.

In this paper, two kinds of control constraint $U_{ad}$ are considered: one is the unconstrained case (UC) with $U_{ad}: = U$, and the other is the box-constrained case (BCC) with $U_{ad}: = \{u(x)\in U \mid u_a\leq u(x)\leq u_b \}$.

Firstly, the existence and uniqueness results for the problem $(P)$ are shown in the following theorem.

Theorem 2.3. Under the assumption (2.1), and $\alpha>0$, the problem $(P)$ has a unique pair of optimal solutions $(\bar{y},\bar{u})$.

Proof. Denote the feasible set by

Since $F\geqslant 0$ and $W$ is nonempty, the infimum

exists and then we can find a minimizing sequence $(y_n,u_n)\subset W$ with

The sequence $\{u_n\}$ is bounded, since $F(y_n,u_n)\geq \frac{\alpha}{2}\|u_n\|_{L^(\partial D)}^2$. $\{u_n\}$ is included in a complet space $H^{1/2+\sigma}(\partial D)$, so there exists a weakly convergent subsequence $\{u_{n_j}\}\subset \{u_n\}$ with $u_{n_j}\rightharpoonup \bar{u}$ as $j\rightarrow \infty$.

It is easy to know that $S$ is a bounded and linear operator, then

It follows from the closedness and convexity of the product space $Y\times U_{ad}$ that $(\bar{y},\bar{u})\in Y\times U_{ad}$, and $(\bar{y},\bar{u})$ satisfies the state equation $\bar{y} = S(\bar{u})$, so $(\bar{y},\bar{u})\in W$. Due to the expression of the cost function (10), $(y,u)\in Y\times U \rightarrow F(y,u)$ is continuous and convex. Hence, $F$ is weakly lower semicontinuous and using that $(\bar{y},\bar{u})$ is feasible yields

Thus, $(\bar{y},\bar{u})$ is the optimal solution of $(P)$. And with $\alpha >0$ and the linearity of the operator $S$, it is easy to find that $u \mapsto F(S(u),u)$ is strongly convex, which implies the uniqueness.

Next, we consider the optimality conditions for the optimal control problem $(P)$ under unconstrained case and box-constrained case. The problem $(P)$ can be rewrited in the following reduced form

where $G(u): = F(S(u),u) = \mathbb{E}\left[\frac{1}{2}\int_D |S(u)-y_d|^2 dx\right]+\frac{\gamma}{2}\int_{\partial D} u^2 dx$.

For the unconstrained case (UC): Under this condition, the optimal control satisfies

where $\bar{u}\in U$ is the optimal control and $G'(\bar{u})$ is the Fr$\acute{e}$chet derivative of $G$ at $\bar{u}$. With the definition of Fr$\acute{e}$chet derivative, the linearity of expectation operator $\mathbb{E}$ and the operator $S$, one can deduces the first order necessary condition

where $S^*$ is the adjoint of operator $S$.

For the box-constrained case (BCC): The bounded optimal control satisfies the following lemma,

Lemma 2.4 ([20]). Suppose $G: U_{ad}\mapsto \mathbb{R}$ is Gateaux-differentiable and $U_{ad}$ is nonempty and convex. If $\bar{u}$ is an optimal control of the problem

then $\bar{u}\in U_{ad}$ and $\bar{u}$ satisfies the variational inequality

With lemma 2.4 and the expression of $G'(u)$ in (12), the first order necessary condition of $(P)$ under box-constrained case can be summarized as

3.   Numerical approximation

To compute numerical solutions of the random optimal control problem $(P)$ in $(9)$, both a random space approximation and a physical space finite element discretisation are needed. In order to formulate a computationally efficient scheme, we introduce a multi-modes expansion (MME) preprocessing technique which will be described in detail later to approximate random variables. Do not like the stochastic collocation method or standard Monte Carlo method, the MME technique simplifies the random space without solving tedious and messy optimized system. The essential feature of MME is to transfer the random field from the coefficients to the right hand side functions. And the advantages of this technique will be more clear after the MME approximate optimization problem is discretized by FEM. Moreover, With an equation about the stochastic state variable the discrete stochastic PDE-constrained optimization problem can be recast as an almost deterministic optimization problem with a mathematical expectation of random coefficients in the cost functional. Next, using Monte Carlo method we can compute the mathematical expectation without spending too many computing resources.

3.1.   Multi-modes expansion approximation scheme(MME)

Recall that $y$ denotes the random state variable, and also it is the real solution to equations (1)-(2). And $y$ has the following multi-modes representation as a power series of the perturbation parameter.

Lemma 3.1 ([15]). Under $4\varepsilon \sqrt{C(D,a)}(1+k_0)<1$, where $C(D,a)$ is given by (6), the random state variable $y$ has the following multi-modes expansion in terms of powers of the perturbation parameter $\varepsilon$,

where $y_q\in Y, q = 0,1,2,\cdots$.

The proof of lemma 3.1 can be found in Feng's paper [15].

Actually, we are more interested in the finite dimensional truncation of the multi-modes expansion forms, which is

The truncation error of $y$ and $y^Q$ will be stated later.

Substituting the expansion of $y^Q$ for $y$ in the equation (1) and boundary condition (2), and match the coefficients of $\varepsilon^q$ order terms for $q = 0,1,2,\cdots,Q-1$ it follows that

with the boundary conditions for each mode function $y_q$ as

It is easy to find that the first equation in (16) is irrelevant to randomness, and $y_0(x)$ is determinate. And we remind that in the following text $y_0(x,\xi)$ always refers to $y_0(x)$.

Let $V = H_0^2(D)$, and relative finite element weak formulations to equations (15)-(17) are defined as follows: Find $y^Q(\cdot,\xi)\in Y$ such that

with $y_0\in\{v(x)\in H_0^2(D)|v(x) = u(x),\forall x\in \partial D\}$, $y_q\in Y, q = 1,\cdots, Q-1$ and

According to the weak formulations, the "MME solution operator" $S^Q:L^2(\partial D) \rightarrow L^2(H^2(D);\Omega)$ is defined by

where $u$ is the Dirichlet boundary function, and $y^Q$ is the relative multi-modes expansion truncation solution solved by the above weak formulations. Obviously, $S^Q$ is a linear operator respect to $u$.

Using lemma 2.2 we can get the following theorem about the convergence analysis between the real solution and the truncated solution.

Theorem 3.2. Assume that $\varepsilon<1$, and $u\in H^{1/2+\sigma}(\partial D)$ for $0<\sigma \leq 1$. Then

which can also be writed as

Similar proof can be found in [15].

Then, with the "MME solution operator" $S^Q$, we can rewrite the optimal control problem $(P)$ as

where the definition of functional $F$ can be found in (10).

3.2.   Finite element discretization

In this part, we deduce the variational discretization for problem $(P^Q)$.

Suppose $\mathcal{T}_{h} = \bigcup_{T\in \mathcal{T}_h} T$ is a triangulation on $D$ which consists of shape regular polygons. Let $h = \max\limits_{T\in \mathcal{T}_h} h_T$, where $h_{T}$ is the maximal diameter of the polygon $T$. Let the set of edges under $\mathcal{T}_h$ be denoted by $\mathcal{E}_h = \bigcup_{e\in \mathcal{E}_h} e$.

Associated with the triangulation, we consider the finite element space

where $\mathcal{P}_1$ denotes the linear polynomial space. The discrete control constraint set is given by

Denote by $\{x_l\}_{l = 1}^{N}$ and $\{\phi_j\}_{j = 1}^{N}$ the nodal points and the corresponding basis functions of $V_h$ respectively, where

For any $\xi\in \Omega$, the FEM interpolations of $y^Q(x,\xi)$ in state equations (18) is:

with $y_j(\xi) = y^Q(x_j,\xi)$ for $j = 1,2,\cdots,N$.

with $y_{0,j} = y_{0}(x_j)$ for $j = 1,2,\cdots,N$.

And for any $\xi\in \Omega$, the FEM interpolations of $y_q(x,\xi)$ in equation (19) is:

with $y_{q,j}(\xi) = y_{q}(x_j,\xi)$ for $j = 1,2,\cdots,N$ and $q = 1,\cdots,Q-1,$.

We also use the same basis functions to express the FEM interpolation of boundary control function $u_h(x)$, that is

where

Then the finite element scheme for the weak formulations (18)-(19) is given by: For each realization $\xi\in \Omega$, seeking $y_h^Q(\cdot,\xi)\in V_{h}$ such that

where $y_{0,h}\in W_{hu}$ and $y_{q,h}(x,\xi)(q = 1,2, \cdots,Q-1)\in V_h$ for fixed $\xi\in \Omega$ satisfy the following equations,

Notice that the discrete control variable $u_h$ is implicit in the solution space $W_{hu}$ that $y_{0,h}$ belongs to, which makes $y_h^Q(x,\xi)$ in (26) related to $u_h$. And (26)-(27) can be summarized as

where $S_h^Q: L^2(\partial D)\rightarrow Y_h$ is a random operator.

Next we focus on the regularity estimate of $S_h^Q$ and error estimate between $S$ and $S_h^Q$.

Lemma 3.3. From [15], we can summarize that:

(a) For any given $u_h\in U_h$, suppose $\{y_{q,h}\}_{q\geq 0}$ are a series of solutions solved by equations (27) with respect to $u_h$. Then the regularity estimates for $y_{q,h}$ are

(b) For any given $u_h\in U_h$, suppose $y^Q$ and $y_h^Q$ are the solutions in equations (18)-(19) and equations (26)-(27) relate to $u_h$ respectively, which means $y^Q = S^Q(u_h)$ and $y_h^Q = S_h^Q(u_h)$. Then we have the error estimate as

From the inequalities (28), we can deduce that

where $C_1 = \sum_{q = 0}^{Q-1} C(D,a)^{q+1}\varepsilon^{jq}$. In another word, it means

Combining theorem 3.2 with inequality (29), we get the error estimate between $S$ and $S_h^Q$.

Now, the stochastic optimal control problem $(P^Q)$ in (21) has a discretization scheme as

Here, we give the one-order necessary optimality conditions of ($P^Q_h$) which will be used later,

in which $\bar{u}_h$ is the optimum solution of the problem $(P^Q_h)$.

For simplicity of calculation, we rewrite the FEM functions to relative vector forms. Suppose there is a vector function $\textbf{R}(x) = [\phi_1(x),\phi_2(x),\cdots,\phi_N(x)]$. Then from (22) the function $y_h^Q$ in (26) can be written in the following form

where

Similarly, there have

where

and

where

And from (25) we set

with

By substituting (35)-(38) into equations (26)-(27) and eliminating $\textbf{R}(x)$, we get

with

where

Finally, it is necessary to mention that equations (39)-(40) are equivalent to equations (26)-(27) in different formats.

3.3.   A determinant optimal control problem derived from $(P^Q_h)$

In this part, we will show that the problem ($P^Q_h$) in (32) can be approximated to a discrete optimal control problem with determined state and control variables. The key to this transformation lies in equality between the random state variable $y^Q_h(x,\xi)$ and the determined finite dimensional variable ${\bf y_{0,h}}$. First, let us study on the relationship between ${\bf y^Q_h}(\xi)$ and ${\bf y_{0,h}}$.

In equations (40), due to the singularity of stiffness matrix $\textbf{A}$ it holds

Then bringing (41) into (39), we can compute that

where ${\bf H_Q}(\xi)$ is a random matrix defined as

Notice the equality (42) that the random field and physical space in ${\bf y^Q_h}(\xi)$ are separated into a random matrix ${\bf H_Q}(\xi)$ and a a non-random variable ${\bf y}_{\bf 0,h}$.

With the help of equality (42), we are ready to get the new optimal control problem. By substituting $y_h^Q(x,\xi) = \textbf{R}(x){\bf y^Q_h}(\xi) = \textbf{R}(x){\bf H_Q}(\xi){\bf y_{0,h}}$ and $u_h(x) = {\bf R}(x){\bf u_h}$ in the cost functional $F(y^Q_h(x,\xi),u_h(x))$ in (32), it has

where

In the above equation, randomness only exists in $\mathbb{E}\big[{\bf H_Q}(\xi)^T{\bf R}(x)^T{\bf R}(x) {\bf H_Q}(\xi)\big]$ and $\mathbb{E}\big[{\bf H_Q}(\xi)^T\big]$ which can be estimated approximately by Monte Carlo simulation at the very beginning.

Next we list some theoretical conclusions of Monte Carlo simulation in order to get further global error estimates of our algorithms.

Suppose $\{\xi_i\}_{i = 1}^{M}$ are arbitrary independent, identically distributed samples in the probability space $(\Omega, \mathcal{F}, P)$. For given $v(x,\xi)\in Y$, define an approximation of the expectation $\mathbb{E}[v(x,\xi)]$ as

And we can find the following estimate in [2],

Applying the Monte Carlo simulation (44) to the expanded functional $F_h({\bf y_{0,h}}, {\bf u_h})$, we get an equivalent optimal control problem $(\tilde{P}_h^Q)$ to $(32)$

where

and

Compared to the stochastic optimal control problem $(P^Q_h)$ in (32), the new problem $(\tilde{P}_h^Q)$ possesses a simpler composition with deterministic state variable, control variable, and state equation. Although these two problems have different formation, they are almost equivalent. Major change is the substitution of mathematical expectation for its Monte Carlo approximation. The whole transformation from $(P^Q_h)$ to $(\tilde{P}_h^Q)$ is carried out under the equality (42). And most importantly the control variables are identical, which leads to the similar first order necessary conditions between the two problems.

The following proposition reveals the first order necessary condition of problem $(\tilde{P}_h^Q)$, which is derived by substituting the Monte Carlo simulation $E_M$ for the mathematical expectation $\mathbb{E}$ in (33) and (34).

Proposition 1. Problem $(\tilde{P}_h^Q)$ admits a unique FEM optimal solution $(\bar{y}_{0,h},\bar{u}_h)$ with $\bar{y}_{0,h} = {\bf R}(x){\bf \bar{y}_{0,h}},\bar{u}_h = {\bf R}(x){\bf \bar{u}_h}$, and the optimum control $\bar{u}_h$ satisfies the following one order necessary condition:

4.   The ADMM for ($\tilde{P}_h^Q$) under unconstrained and box-constrained constrains

In this section, we shall solve the optimal control problem $(\tilde{P}_h^Q)$ in (46) by an alternating direction method of multipliers (ADMM). The origin of ADMM can be traced back to the 1950s for nonlinear variational problems [17]. In recent years, it is mainly applied to sparse signal recovery, image restoration, compressed sensing, machine learning and so on [29,31]. Besides, ADMM is an efficient numerical algorithm for convex optimization problem, which combines the benefits of dual decomposition and augmented Lagrangian methods. Relative theoretical results of ADMM can be found in [6,7,12,18,19].

The advantages of ADMM for solving the PDE-constrained optimization $(\tilde{P}_h^Q)$ can be summarized as follows: (1) ADMM enjoys a global convergence with a convergence rate of $O(\frac{1}{k})$, where $k$ denotes the number of iterations, compared with the local convergence method such as Newton's method or SQP method. (2) ADMM makes use of the character that the control variable and the state variable are decoupled in the cost functional, which will reduce the computational cost. (3) The stiff matrix of left hand side is invariant in each iteration, so that we only need to calculate the inverse matrices for one time in the whole process which reduces the computing time of solving the algebraic system sharply [33].

Next we will present two ADMM algorithms for unconstrained case and box-constrained case separately.

4.1.   Unconstrained case

Because the cost functional in (46) has variables separable structure, it can be divided into two irrelevant items. Let

Then the augmented Lagrangian function for problem $(\widetilde{P}^Q_{h})$ is given by

where ${\bf \lambda_h}\in V_h$ is the Lagrange multiplier and $\rho$ is the penalty parameter. Similarly, the problem $(\widetilde{P}^{Q}_{h})$ is equivalent to the saddle-point problem

Let $\tau$ be the tolerance error and $\tau^{k} = \|{\bf u}^{k+1}_{\bf h}-{\bf u}^{k}_{\bf h}\|+\|{\bf \lambda}^{k+1}_{\bf h}-{\bf \lambda}^{k}_{\bf h}\|$ be the iterative error at $k$-th step. Then the ADMM for solving (54) under unconstrained case is given by the following algorithm 1.

One of the advantages for applying ADMM to the optimal control problem $(\tilde{P}_h^Q)$ is that ${\bf y}^{k+1}_{\bf 0, h}, {\bf u}^{k+1}_{\bf h}, {\bf \lambda}^{k+1}_{\bf h}$ are easy to figure out in each iteration step. By simple calculation, we obtain

where ${\bf W} = E_M\Big[{\bf H_Q}(\xi)^T{\bf D}{\bf H_Q}(\xi)\Big]$, ${\bf w} = E_M\Big[{\bf H_Q}(\xi)^T{\bf c}\Big]$ and the matrix ${\bf D}$ and vectors ${\bf c}$ are defined as

4.2.   Box-constrained case

For the box-constrained case, there should bring in additional constraints ${\bf x_h}$ and an indicative function $\Omega({\bf x_h})$ as

Secondly, to get the ADMM form under the box-constrained case, we regard ${\bf y}_{\bf 0, h}$ and ${\bf x_h}$ as a global variable ${\bf \omega} = ({\bf y}_{\bf 0, h},{\bf x_h})$. Thus the discrete optimal control problem (46) under box-constrained case can be rewritten as

Suppose

Then the equivalent saddle-point problem is

Let $\tau_2^{k} = \|{\bf u}^{k}_{\bf h}-{\bf u}^{k-1}_{\bf h}\|+\| {\bf \lambda}^{k}_{\bf 1h}-{\bf \lambda}^{k-1}_{\bf 1h}\|+\| {\bf \lambda}^{k}_{\bf 2h}-{\bf \lambda}^{k-1}_{\bf 2h}\|$, therefore the two block ADMM algorithm for the saddle-point problem (54) is

The sub-optimal problems in algorithm 2 can be computed directly without iteration. Let $P_{{\bf U_{ad}^h}}(v)$ be the projection of $v$ onto ${\bf U_{ad}^h}$, i.e.

Then, the fomulas in each iterative step of algorithm 2 turns to be

where $\textbf{I}$ is the identity matrix.

We finish this part by giving the convergence rate of the general ADMM algorithm. Because algorithm 1 and algorithm 2 have essentially the same structures, they share common convergence results.

Lemma 4.1 ([8]). The sequence $({\bf y}^{k}_{\bf 0, h}, {\bf u}^{k}_{\bf h},{\bf \lambda}^{k}_{\bf h})$ in both algorithm 1 and 2 converges to the optimal solutions $({\bf \bar{y}_h}, {\bf \bar{u}_h},{\bf \bar{\lambda}_h})$. And ${\bf \lambda}^{k}_{\bf h}$ denotes the vector $(\lambda^{k}_{\bf 1h} ,\lambda^{k}_{\bf 2h})$ in the second algorithm. Let ${\bf v}^{k}_{\bf h} = ({\bf u}^{k}_{\bf h},{\bf \lambda}^{k}_{\bf h})$ and ${\bf \bar{v}}_{\bf h} = ({\bf \bar{u}}_{\bf h},{\bf \bar{\lambda}}_{\bf h})$. Both algorithms are contractive and they satisfy the following estimate

where ${\bf v}^{0}_{\bf h} = ({\bf u}^{0}_{\bf h},{\bf \lambda}^{0}_{\bf h})$ is the initial value of the ADMM.

Remark 1 ([8]). Let $\bar{{\bf u}}_{\bf h}$ be the real solution of the problem ($\tilde{P}_h$) in (46), and ${\bf u}^{\bf k}_{\bf h}$ be the numerical solution after $k$ steps iteration in both algorithm 1 and algorithm 2. Then in general case, Lemma 4.1 implies that the theoretical convergence rate for the ADMM is

where $\mathcal{H}$ is some given matrix related to $B$ and $\rho$ (concrete definition can be found in [8]). The computational results always show faster convergence rate than $\mathcal{O}(1/k)$ (see numerical simulations in Section 6).

5.   Convergence analysis

This section provides some error estimates for the proposed algorithms in the foregoing discussion. Suppose $\bar{u},{\bf \bar{u}_h},{\bf u^k_h}$ denote the optimum control for the problem $(P)$ in (9), the real solution for the problem $(\tilde{P}_h^Q)$ in (46), and the iterative solution for algorithm 1 and algorithm 2. The ultimate goal is to derive the global error estimate $\|\bar{u}-{\bf R}(x){\bf u^k_h}\|_{L^2(\partial D)}$ for algorithm 1 and algorithm 2. For this purpose, firstly we will show the discretization error estimate for $\|\bar{u}-\bar{u}_h\|_{L^2(\partial D)}$, where $\bar{u}_h = {\bf R}(x){\bf \bar{u}_h}$ is the corresponding FEM function of ${\bf \bar{u}_h}$. Then with the existing ADMM convergence rate (56), the global error estimate is derived.

5.1.   Discretization error estimates

Let us start with the unconstrained case. Similar error estimates have been discussed in our earlier article [26], so we draw the following conclusion directly.

Theorem 5.1. In algorithm 1 for the optimal control problem $(\tilde{P}_h^Q)$, let $M$ be the number of samples for MC simulation, $h$ be the FEM mesh size, $\varepsilon$ be the magnitude of the random perturbation, and $Q$ be the numbers of MME expansion. Then, algorithm 1 enjoys the following conclusion,

where $K_1$ and $K_2$ are constants independent of $M$ and $h$, and $0<\sigma \leq 1$.

Proof. With the help of first order necessary conditions (12) for $(P)$ and (47) for $(\tilde{P}_h^Q)$ under unconstrained case, lemma 3.3, the Monte Carlo error estimate in (45), we can get the estimate. For a rigorous proof the reader is referred to [26].

Next, the box-constrained case will be considered. For ease of understanding, we transcribe the one order necessary conditions (13) for $(P)$ and (48) for $(\tilde{P}_h^Q)$ under box-constrained case below,

To estimate $\|\bar{u}-\bar{u}_h\|_{L^2(\partial D)}$ from (58) and (59), common skill is to insert $u = \bar{u}_h$ in (58), and if we can also take $v_h = \bar{u}$ in (59), the estimation will be trivial to carry out by adding these two inequalities together. Unfortunately, this can not be realized because $\bar{u}$ does not belong to $U_{ad}^h$. Therefore, a discrete control $\hat{u}_h\in U_{ad}^h$ will be constructed to approximate $\bar{u}$. Especially, the variable $\hat{u}_h$ should satisfies $G'(\bar{u})(\bar{u}-\hat{u}_h) = 0$ for subsequent estimates where the functional $G$ is defined in (11).

Following the idea in [10], we construct the discrete control $\hat{u}_h\in U_{ad}^h$ with $\hat{u}_h(x_j) = u_{h,j}$ for every node $x_j\in \partial D$ as

where

The following lemma reveals that $\hat{u}_h$ satisfies some requirements.

Lemma 5.2 ([10]). $\hat{u}_h\in U_h$ defined by 60 obeys the following properties:

1. $\hat{u}_h\in U_h^{ad}$.

2. $G'(\bar{u})(\bar{u}-\hat{u}_h) = 0$.

3. The approximation property

where $p\geq 2$.

With the help of the function $\hat{u}_h$ we can now prove the following theorem.

Theorem 5.3. Suppose $M$ is the number of samples for MC simulation, $h$ is the FEM mesh size, $\varepsilon$ is the magnitude of the random perturbation, and $Q$ is the numbers of MME expansion in the optimal control problem $(\tilde{P}_h^Q)$ in (46). Then, for the problem $(\tilde{P}_h^Q)$ under box-constrained case it has the following conclusion,

where $p\geq 2$.

Proof. Insert $\bar{u}_h\subset U_{ad}$ and $\hat{u}_h\in U^{ad}_h$ into (13) and (48) respectively, then we have

For simplicity, we omit the subscript of the inner product $\langle \cdot,\cdot \rangle_{U',U}$ in the following. Then, by adding these two inequalities and re-ordering the result it gets

We divide the right hand of the above inequality in four terms and estimate them step by step. For the first term $\textbf{I}$, it follows from the estimate (45), the Cauchy-Schwarz inequality, and the triangle inequality that

For further estimation of $\textbf{I}$, we should use the following inequality

Under the above inequality, the estimate (30) when $j = 4$, and the approximation property (61) it can be deduced that

For the second term Ⅱ, it turns to be

Next, using the Cauchy-Schwarz inequality, the error estimate (31), the estimate of $S$ in (8), and (30) when $j = 2$, it is easy to check that

The rest is to estimate Ⅳ. Before proceeding further, let us see an important equation which comes from the second property of $\hat{u}_h$ in Lemma 5.4.

With the above equation, we are now turning to the following estimation.

where $0<\sigma\leq 1$ and $p\geq 2$.

Substituting the estimates for $I, II, III, IV$ into the inequality (65), we arrive at

With Young's inequality, it turns to be

By rearranging the above inequality, we get

Thus

which completes the proof.

5.2.   Global error estimate

Finally, we summarise the above conclusions and arrive at the following global error estimate of the control variables.

Lemma 5.4. Let $\bar{u}$ be the optimum control for the problem $(P)$ in (9) under either the unconstrained case (UC) or the box-constrained case (BCC). And ${\bf u^k_h}$ is the optimal iterative solution for both algorithm 1 and algorithm 2. It has the following global error estimate

where $C$ is a generic constant, and

with $0<\sigma \leq 1$ and $p\geq 2$.

The above lemma can be easily proved by theorem 5.3, theorem 5.3, remark 1 and the triangle inequality

where $\bar{u}_h = {\bf R}(x){\bf \bar{u}}_{\bf h}$ and $\mathcal{H}$ is some bounded matrix defined as

6.   Numerical simulations

In this section, we present some numerical tests to show the performance of algorithm 1 and algorithm 2. Let $D = [-1, 1]\times[-1, 1]\subset \mathbb{R}^2$, and consider

We are going to solve the following optimal control problems constrained by random Poisson equations,

where

In the following, we will show the convergence rates of FEM, ADMM, and MME for algorithm 1 and algorithm 2 separately. As it is commonly used in other articles, we choose the random variable to be uniformly distributed in $[-0.5,0.5]$. The magnitude of the random fluctuation $\varepsilon$ will set to be different values in the following tests. First, we will get the mean value of coefficient by Monte Carlo method with 1000 samples (which can guarantee that the error of MC is around $3.2\times 10^{-3}$) and discretize the problem by continuous piecewise linear finite elements. The resulting deterministic finite dimensional optimization problem is solved with ADMM. We use classical ADMM for the unconstrained case, and two-block ADMM for the box-constrained case.

Test 1: The convergence rates of FEM

We start by examining the convergence rates of FEM numerically. The FEM mesh sizes are set to be $h = 1/2,1/4,1/8,1/16,1/32$, and $1/64$ respectively. Note that for this problem we do not have an exact solution, therefore the error is computed with respect to a reference numerical solution derived under $h = 1/128$. Let the iterative steps $k\approx 200$ in both unconstrained case and box-constrained case to guarantee that the error of ADMM is about $10^{-4}$. The perturbation magnitude $\varepsilon = 0.1$, and we use three modes (i.e., Q = 3) in MME approximation. Figure 1 gives the FEM convergence rates of the numerical solutions ${\bf u}^{\bf k}_{\bf h}$ and ${\bf y}^{\bf k}_{\bf h}$ to the algorithm 1 (left) and algorithm 2 (right) correspondingly. From figure 1, we can find that the FEM convergence rate of $u$ in unconstrained case is approximate to order 2, and the FEM convergence rate of $u$ is about order 1.5, which are both higher than the theoretical results.

Test 2: The convergence rates of ADMM

Consider the convergence rates of ADMM with H-norm numerically in both cases. The perturbation magnitude $\varepsilon$ is setted to be 0.1. In figure 2, three modes in MME are used. And the FEM mesh size is fixed to be $h = 1/64$ to guarantee the ADMM errors $\|{\bf \bar{u}}_{\bf h}-{\bf u}^{\bf k}_{\bf h}\|_\mathcal{H}$ being the major error source. Figure 2 shows the ADMM errors for two algorithms with respect to the $k$-iteration in log-scale. In these results we can observe that the ADMM convergence rate is much faster than $O(\frac{1}{\sqrt{k}})$ in algorithm 1. This is possible, since the control in the model problem under unconstrained case is smooth. Comparing to algorithm 1, the ADMM convergence rate is lower in algorithm 2 because of the non-smoothness of the control variable $u$ under the box-constrained condition.

Test 3: The convergence rates of MME

The goal of this test is to validate the convergence rates of MME numerically with respect to the perturbation magnitude $\varepsilon$ and the truncated number $Q$. The mesh size of the FEM is $h = 1/64$, and the tolerance error is set to $10^{-4}$ (in which the iterative number of ADMM is about $k = 200$). First, we fix $Q = 3$ to observe the convergence rate about $\varepsilon$ alone. Figure 3 shows the convergence rates is about order 3 which is consistent with the theory. Then, fix $\varepsilon = 0.1$ to see the error of control variable $u$ under different modes. Table 1 gives the results in which we can find that the bigger the number of modes is, the lower is the error of $u$. But, it also cost more CPU time. So we just using three modes for the two algorithms.

7.   Conclusion

In this paper, we presented two numerical algorithms for solving optimal control problems constrained by random Poisson equations under unconstrained condition and box-constrained condition. Both algorithms adopt MME technique for the random space and FEM for physical space discretization. Then under a special state variable substitution (42), the random constrained control problem can be approximated to a deterministic discrete optimization problem. Further in algorithm 1 and algorithm 2, we use ADMM and two-block ADMM to solve the fully discretized problem respectively. A detailed analysis is carried out for the error estimation of the discrete optimal solution under some mild assumptions on the random input data. Because it only needs to solve a deterministic constrained control problem, our algorithms can reduce the memory consumption sharply, especially for high-dimensional and large-scale problems which require a great number of samples.

Acknowledgments

This work of Jingshi Li was supported by the Natural Science Foundation of Jiangsu Province (Grant No. BK20190766), the Startup Foundation for Introducing Talent of NUIST (No. 2018r022), and the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (No. 19KJB110018). The work of Jiachuan Zhang, Guoliang Ju and Juntao You was partially supported by the NSF of China No. 11971221 and 11731006, the Shenzhen Sci-Tech Fund No. JCYJ20170818153840322 and JCYJ20190809150413261, and Guangdong Provincial Key Laboratory of Computational Science and Material Design No. 2019B030301001.