# Numerical Optimal Control for Problems with Random Forced SPDE Constraints

Authors: R. Naseri; A. Malek
Journal: ISRN Applied Mathematics
 (2014)

---

## Abstract

A numerical algorithm for solving optimization problems with stochastic diffusion equation as a constraint is proposed. First, separation of random and deterministic variables is done via Karhunen-Loeve expansion. Then, the problem is discretized, in spatial part, using the finite element method and the polynomial chaos expansion in the stochastic part of the problem. This process leads to the optimal control problem with a large scale system in its constraint. To overcome these difficulties the adjoint technique for derivative computation to implementation of the optimal control issue in preconditioned Newton’s conjugate gradient method is used. By some numerical simulation, it is shown that this hybrid approach is efficient and simple to implement.


---

## Body

## 1. Introduction

Physical problems, in many cases, can be formulated as optimization problems. These problems are utilized to gain a more widely understanding of physical systems. Typically, these problems depend on some models which, in many cases, are deterministic. Uncertainty might plague everything from modeling assumptions to experimental data. As such, in order to accommodate for these uncertainties, many practitioners have developed stochastic models. In order to make sense of and solve these models, in addition to randomness of the models, some additional theory is required to the resulting optimization problems. This paper proposes an adjoint based approach for solving optimization problems governed by stochastic diffusion equations. In order to deal with the stochastic partial differential equations (SPDE) as a constraint in an optimization problem one may first solve the SPDE. Since the provided systems by this approach are random and nonlinear, such kind of problems are difficult to handle and very challenging. One of the most challenging examples in this area is the control of stochastic diffusion equations with random forcing [1]. Polynomial chaos expansion (PCE) [2] provides a good direction for solution of nonlinear SPDEs numerically. In the presence of the random forcing, PCE seems to be more accurate and efficient numerical method than Monte Carlo simulation. In fact, the PCE can be interpreted as the Fourier expansion in the probability space. Particularly, the aim of this work is the numerical solution for the distributed control problems involving the stochastic diffusion equations in the form:
(1)
-
∇
·
(
a
(
x
,
ω
)
∇
u
(
x
,
ω
)
)
=
z
(
x
)
in
D
×
Ω
(2)
u
(
x
,
ω
)
=
0
on
⁡
∂
D
×
Ω
,
where D is the spatial domain, ∂
D is the boundary of the spatial domain, Ω is probability space, x
∈
D, ω
∈
Ω, u is the solution of the SPDE, z is the source function, and a
(
x
,
ω
) is the permeability field of the problem. As an important assumption for the stochastic diffusion equation (1), it is assumed that the random coefficient a
(
x
,
ω
) satisfies the elliptic condition [3]; that is, there exists a constant a
min
⁡ such that
(3)
0
<
a
min
⁡
≤
a
(
x
,
ω
)
∀
(
x
,
ω
)
∈
D
×
Ω
.Karhunen-Loeve expansion (KLE) of correlated random functions is used to separate the random and deterministic parts in random coefficient [4, 5]. Therefore, a finite dimensional approximation a
M
(
x
,
ω
) is performed by truncating the Karhunen-Loeve expansion of the permeability field a
(
x
,
ω
). Then, (1) is approximated by
(4)
∇
(
a
M
(
x
,
ω
)
∇
u
)
=
z
.For approximation of the functionu with Nrandom variables and the Gaussian random variables of the highest degree K, the number of PCE coefficients are(
N
+
K
)
!
/
N
!
K
!, [6]. Now, weak formulation of (4) and then Galerkin finite element method are used to solve the SPDE problem approximately. Generally in nature, optimization with SPDE-constraint is infinite dimensional, large, and complex. There are two numerical approaches for solving this problem [2, 7]:(i)
discretize then optimize (DO): for problems that can be trivially discretized first,
(ii)
optimize then discretize (OD): for problems that are differentiable.Our intention is to work with discretizethen optimization algorithms for smooth functions. Thus, we follow the DO approach. Galerkin finite element method incorporated with polynomial chaos expansion is used for discretizing the SPDE. By computing gradient and Hessian of the Lagrangian augmented function, preconditioned Newton’s conjugate gradient method is used to compute the best right hand side vector (control values at grid points), where the corresponding solution of SPDE (state values) according to this vector, stays as near as to the desired value approximation, and has lowest cost in the objective function. The organization of this paper is as follows: in Section2, problem formulation in addition to KLE, PCE, and stochastic Galerkin method is presented. In Section 3, distributed control of random forced diffusion equation with some numerical examples is proposed.


## 2. Problem Formulation

We consider optimal control problems of the form(5)
min
⁡
u
∈
U
,
z
∈
ℤ
J
(
u
,
z
)
s
.
t
.
e
(
u
,
z
)
=
0
,
where J
:
H
0
1
(
D
)
⨂
L
2
(
Ω
)
→
ℝ is the objective function, e
:
U
=
H
0
1
(
D
)
⨂
L
2
(
Ω
)
→
ℤ
=
L
2
(
D
) is an operator which stands as the constraint, D is the spatial domain, Ω is the probability space, and ⨂ is the tensor product.To provide a concrete setting for our discussion, the following problem is considered as the constraint operator(6)
-
∇
·
(
a
(
x
,
ω
)
∇
u
(
x
,
ω
)
)
=
z
(
x
)
in
D
×
Ω
,
u
(
x
,
ω
)
=
0
on
⁡
∂
D
×
Ω
,
where a
(
x
,
ω
) is the correlated random fields. Hence, the solutions u
(
x
,
ω
) of the SPDE (6) are also random fields. It is assumed that the control space z
(
x
)
∈
L
2
(
D
). The weak form of the constraint function (6) can be defined for all v
∈
U as
(7)
〈
e
(
u
,
z
)
,
v
〉
U
*
,
U
=
∫
Ω
ρ
(
ω
)
∫
D
(
a
(
x
,
ω
)
∇
u
(
x
,
ω
)
·
∇
v
(
x
,
ω
)
d
x
d
ω
-
z
(
x
)
v
(
x
,
ω
)
)
d
x
d
ω
.It is assumed that the objective or cost functional is as(8)
J
(
u
,
z
)
=
E
[
1
2
∥
u
(
x
,
ω
)
-
u
^
(
x
)
∥
L
2
(
D
)
2
]
+
α
2
∥
z
(
x
)
∥
L
2
(
D
)
2
,
where E
(
·
) denotes the expected value, u
^ is a given desired function, and α is the regularization parameter. Obviously, controlling the solution of this problem leads to controlling the statistics of the solution, for example, the expected value of the solution, given by
(9)
∥
E
[
u
(
x
,
ω
)
-
u
^
(
x
)
]
∥
L
2
(
D
)
2
≤
E
[
∥
u
(
x
,
ω
)
-
u
^
(
x
)
∥
L
2
(
D
)
2
]
≤
∥
u
(
x
,
ω
)
-
u
^
(
x
)
∥
L
2
(
D
)
⊗
L
2
(
Ω
)
2
.So, the optimal control problem can be interpreted as follows: given the random fielda
(
x
,
ω
), minimize the objective function (8) over all z
∈
ℤ, subject to u
∈
U, such that for the given random field a
(
x
,
ω
) satisfies in the weak formulation (7) [1].


### 2.1. Basic Definitions

In (5), J and e are assumed to be continuously F-differentiable and that for each z
∈
ℤ the state equation e
(
u
,
z
)
=
0 possesses a unique corresponding solution u
(
z
)
∈
U. Similar to the deterministic PDEs, existence and uniqueness of the solution for these kind of SPDEs can be established using the Lax-Milgram theorem. Thus, there is a solution operator z
∈
ℤ
→
u
(
z
)
∈
U. In addition, it is assumed that e
u
(
u
(
z
)
,
z
)
∈
ℒ
(
U
,
ℤ
) is continuously invertible. Then, existence and uniqueness of the solution for the above optimal control problems as well as the continuously differentiability of u
(
z
) are ensured by the implicit function theorem [8].


### 2.2. Karhunen-Loeve Expansion (KLE)

Consider a random fielda
(
x
,
ω
), x
∈
D, with finite second order moment
(10)
∫
Ω
E
[
a
2
(
x
,
ω
)
]
d
x
<
∞
.Assume thatE
[
a
]
=
a
¯
(
x
). It is possible to expand a
(
x
,
ω
), for a given orthonormal basis {
ψ
k
} in L
2
(
D
), as a generalized Fourier series
(11)
a
(
x
,
ω
)
=
a
¯
(
x
)
+
∑
k
=
1
∞
a
k
(
ω
)
ψ
k
(
x
)
,
where
(12)
a
k
(
ω
)
=
∫
Ω
a
(
x
,
ω
)
ψ
k
(
x
)
d
x
,
k
=
1,2
,
…
are random variables with zero means. It is important, now, to find a special basis {
ϕ
k
} that makes corresponding a
k uncorrelated: E
[
a
i
a
j
]
=
0 for all i
≠
j. Denoting the covariance function of a
(
x
,
ω
) by R
(
x
,
y
)
=
E
[
a
(
x
,
ω
)
a
(
y
,
ω
)
], the basis functions {
ϕ
k
} should satisfy
(13)
E
[
a
i
a
j
]
=
∫
D
ϕ
i
(
x
)
d
x
∫
D
R
(
x
,
y
)
ϕ
j
(
y
)
d
y
=
0
,
i
≠
j
.Completion and orthonormality of{
ϕ
k
} in L
2
(
D
) follow that ϕ
k
(
x
) are eigenfunctions of R
(
x
,
y
):
(14)
∫
D
R
(
x
,
y
)
ϕ
j
(
y
)
d
y
=
λ
j
ϕ
j
(
x
)
,
j
=
1,2
,
…
,
where λ
j
=
E
[
a
j
2
]
>
0. Indeed, by choosing basis functions ϕ
k
(
x
) as the solutions of the eigenproblem (14), the random variables a
k
(
ω
) will be uncorrelated. In (14), denoting θ
k
=
a
k
/
λ
k, yields the following expansion:
(15)
a
(
x
,
ω
)
=
a
¯
(
x
)
+
∑
k
=
1
∞
λ
k
θ
k
(
ω
)
ϕ
k
(
x
)
,
where θ
k satisfy E
[
θ
k
]
=
0 and E
[
θ
i
θ
j
]
=
δ
i
j. In the case that a
(
x
,
ω
) is considered as a Gaussian process, a
k
(
ω
), k
=
1,2
,
… will be independent Gaussian random variables. The expansion (15) is known as the Karhunen-Loeve expansion (KLE) of the stochastic process a
(
x
,
ω
).Using the KLE (15), the stochastic process can be represented as a series of uncorrelated random variables. Since the basis functions ϕ
k
(
x
) are deterministic, the spatial dependence of the random process can be resolved by them. The convergence property of the KLE to the random process a
(
x
,
ω
) can be represented in the mean square sense
(16)
lim
⁡
N
→
∞
⁡
∫
D
E
|
a
(
x
,
ω
)
-
a
N
(
x
,
ω
)
|
2
d
x
=
0
,
where
(17)
a
N
=
a
¯
(
x
)
+
∑
k
=
1
N
λ
k
θ
k
ϕ
k
is a finite term KLE [4, 5].


### 2.3. Polynomial Chaos Expansion (PCE)

There are problems, as the solution of a PDE with random inputs, that the covariance function of a random processu
(
x
,
ω
), x
∈
D is not known. The solution of such problems can be represented using a polynomial chaos expansion (PCE) given by
(18)
u
(
x
,
ω
)
=
∑
k
=
1
∞
u
k
(
x
)
Ψ
k
(
ξ
)
,
where the functions u
k
(
x
) are deterministic coefficients, ξ is a vector of orthonormal random variables, and Ψ
k
(
ξ
) are multidimensional orthogonal polynomials that satisfy in the following properties:
(19)
〈
Ψ
1
〉
≡
E
[
Ψ
1
]
=
1
,
〈
Ψ
k
〉
=
0
,
k
>
1
,
〈
Ψ
i
Ψ
j
〉
=
h
i
δ
i
j
.The convergence property of PCE for a random quantity inL
2 is ensured by the Cameron-Martin theorem [6, 9]; that is,
(20)
〈
u
(
x
,
ω
)
-
∑
k
=
1
∞
u
k
(
x
)
Ψ
k
(
ξ
)
〉
⟶
L
2
0
.Hence, this convergence justifies a truncation of PCE to a finite number of terms,(21)
u
(
x
,
ω
)
=
∑
k
=
1
P
u
k
(
x
)
Ψ
k
(
ξ
)
,
where the value of P is determined by the highest degree of polynomial K, used to represent u, and the number N of random variables—the length of ξ—with the formula P
+
1
=
(
N
+
K
)
!
/
N
!
K
!. Generally, the value of N is the same as the number of uncorrelated random variables in the system or equivalently the truncation length of the truncated KLE. Typically, the value of K is chosen by some heuristic method. Indeed, in the case of K
=
1 and N random variables, the KLE is a special case of the PCE [9, 10].


### 2.4. Stochastic Galerkin

SupposeW
i
⊂
L
ρ
i
2
(
Ω
i
) with dimension ρ
i for i
=
1
,
…
,
M and V
⊂
H
0
1
(
D
) with dimension N. In addition, let {
ψ
n
i
}
n
=
1
p
i for i
=
1
,
…
,
M be basis of W
i and {
ϕ
i
}
i
=
1
N a basis of V. The finite dimensional tensor product space W
1
⨂
⋯
⨂
W
M
⨂
V can be defined as the space spanned by the functions {
ψ
n
1
1
,
…
,
ψ
n
M
M
ϕ
i
} for all n
1
∈
{
1
,
…
,
p
1
}
,
…
,
n
M
∈
{
1
,
…
,
p
M
} and i
∈
{
1
,
…
,
N
}. For simplification, let n denote a multi-index whose kth component n
k
∈
{
1
,
…
,
p
k
} and I denote the set of such multi-indices; then, the basis functions for the tensor product space have the form
(22)
v
n
i
(
x
,
ω
)
=
ϕ
i
(
x
)
Ψ
n
(
ω
)
,
where
(23)
Ψ
n
=
∏
k
=
1
M
ψ
n
k
k
(
ω
k
)
.Then, the Galerkin method is looking for a solutionu
^
∈
W
1
⨂
⋯
⨂
W
M
⨂
V such that for x
=
(
x
,
y
), all n
∈
I and i
=
1
,
…
,
N,
(24)
∫
Ω
ρ
(
ω
)
∫
D
a
(
x
,
ω
)
∇
u
h
P
(
x
,
ω
)
·
∇
v
n
i
(
x
,
ω
)
d
x
d
ω
=
∫
Ω
ρ
(
ω
)
∫
D
z
(
x
,
ω
)
v
n
i
(
x
,
ω
)
d
x
d
ω
,
whereρ
(
ω
) is the integration weight. Equation (24) is equivalent to the following equation, which results by plugging in the explicit formula for the basis functions
(25)
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
∫
D
a
(
x
,
ω
)
∇
u
h
P
(
x
,
ω
)
·
∇
ϕ
i
(
x
)
d
x
d
ω
=
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
∫
D
z
(
x
,
ω
)
ϕ
i
(
x
)
d
x
d
ω
for all n
∈
I and i
=
1
,
…
,
N.Now, considering the approximate solution as(26)
u
h
P
(
x
,
ω
)
=
∑
j
=
1
N
∑
m
∈
I
u
j
m
ϕ
j
(
x
)
Ψ
m
(
ω
)
and substituting (26) into (25) yields
(27)
∑
j
=
1
N
∑
m
∈
I
u
j
m
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
×
∫
D
a
(
x
,
ω
)
∇
ϕ
j
(
x
)
·
∇
ϕ
i
(
x
)
d
x
d
ω
=
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
∫
D
z
(
x
,
ω
)
ϕ
i
(
x
)
d
x
d
ω
for all m
∈
I and i
=
1
,
…
,
N. Let P
=
|
I
|
=
∏
k
=
1
M
p
k; then, a natural bijection between {
1
,
…
,
P
} and Ican be defined. Thus, (27) can be written as
(28)
∑
m
=
1
P
∑
j
=
1
N
u
m
j
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
×
∫
D
a
(
x
,
ω
)
∇
ϕ
j
(
x
)
·
∇
ϕ
i
(
x
)
d
x
d
ω
=
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
∫
D
z
(
x
,
ω
)
ϕ
i
(
x
)
d
x
d
ω
.Now, ifa and z have the following KLE
(29)
a
(
x
,
ω
)
=
a
0
(
x
)
+
∑
i
=
1
M
ω
i
a
i
(
x
)
,
z
(
x
,
ω
)
=
z
0
(
x
)
+
∑
i
=
1
M
ω
i
z
i
(
x
)
,then, we can rewrite the Galerkin projection (28) as
(30)
∑
m
=
1
P
∑
j
=
1
N
u
m
j
(
∑
k
=
1
M
(
K
k
)
i
,
j
(
K
0
)
i
,
j
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
d
ω
+
∑
k
=
1
M
(
K
k
)
i
,
j
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
ω
k
d
ω
)
=
(
z
0
)
i
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
d
ω
+
∑
k
=
1
M
(
z
k
)
i
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
ω
k
d
ω
,
where
(31)
(
K
k
)
i
,
j
=
∫
D
a
k
(
x
)
∇
ϕ
i
(
x
)
·
∇
ϕ
j
(
x
)
d
x
≅
∑
r
w
r
a
k
(
x
r
)
∇
ϕ
i
(
x
r
)
·
∇
ϕ
j
(
x
r
)
,
(
z
k
)
i
=
∫
D
z
i
(
x
)
ϕ
i
(
x
)
d
x
≅
∑
r
w
r
z
i
(
x
r
)
ϕ
i
(
x
r
)
for k
=
1
,
…
,
M  [3]. Equation (31) is approximated by quadrature rule in spatial domain. Now, assume that there exist functions Ψ
n, n
=
1
,
…
P, such that
(32)
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
d
ω
=
δ
m
n
,
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
ω
k
d
ω
=
C
k
n
δ
n
m
,
where, for m
=
n, the above integral is approximated using quadrature rule in random domain as follows:
(33)
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
n
(
ω
)
ω
k
d
ω
≅
∑
s
w
s
′
ρ
(
ω
s
)
Ψ
n
2
(
ω
s
)
ω
k
.
Then, (30) becomes
(34)
∑
m
=
1
P
∑
j
=
1
N
u
m
j
(
(
K
0
)
i
,
j
+
∑
k
=
1
M
C
k
n
(
K
k
)
i
,
j
)
δ
n
m
=
(
z
0
)
i
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
d
ω
+
∑
k
=
1
M
(
z
k
)
i
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
ω
k
d
ω
or, equivalently,
(35)
∑
j
=
1
N
u
m
j
(
(
K
0
)
i
,
j
+
∑
k
=
1
M
C
k
n
(
K
k
)
i
,
j
)
=
(
z
0
)
i
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
d
ω
+
∑
k
=
1
M
(
z
k
)
i
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
ω
k
d
ω
.In (34) and (35) we have
(36)
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
d
ω
=
0
;
for
n
≠
0
,
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
ω
k
d
ω
=
∫
Ω
ρ
(
ω
)
Ψ
0
(
ω
)
Ψ
n
(
ω
)
ω
k
d
ω
=
C
k
n
δ
0
n
=
C
k
0
.Finally, (35) can be considered as the following block diagonal system:
(37)
(
K
0
+
∑
k
=
1
M
C
k
1
K
k
⋯
0
⋮
⋱
⋮
0
⋯
K
0
+
∑
k
=
1
M
C
k
p
K
k
)
(
u
→
1
⋮
u
→
P
)
=
(
z
→
1
⋮
z
→
p
)
or briefly K
u
→
=
z
→.


## 2.1. Basic Definitions

In (5), J and e are assumed to be continuously F-differentiable and that for each z
∈
ℤ the state equation e
(
u
,
z
)
=
0 possesses a unique corresponding solution u
(
z
)
∈
U. Similar to the deterministic PDEs, existence and uniqueness of the solution for these kind of SPDEs can be established using the Lax-Milgram theorem. Thus, there is a solution operator z
∈
ℤ
→
u
(
z
)
∈
U. In addition, it is assumed that e
u
(
u
(
z
)
,
z
)
∈
ℒ
(
U
,
ℤ
) is continuously invertible. Then, existence and uniqueness of the solution for the above optimal control problems as well as the continuously differentiability of u
(
z
) are ensured by the implicit function theorem [8].


## 2.2. Karhunen-Loeve Expansion (KLE)

Consider a random fielda
(
x
,
ω
), x
∈
D, with finite second order moment
(10)
∫
Ω
E
[
a
2
(
x
,
ω
)
]
d
x
<
∞
.Assume thatE
[
a
]
=
a
¯
(
x
). It is possible to expand a
(
x
,
ω
), for a given orthonormal basis {
ψ
k
} in L
2
(
D
), as a generalized Fourier series
(11)
a
(
x
,
ω
)
=
a
¯
(
x
)
+
∑
k
=
1
∞
a
k
(
ω
)
ψ
k
(
x
)
,
where
(12)
a
k
(
ω
)
=
∫
Ω
a
(
x
,
ω
)
ψ
k
(
x
)
d
x
,
k
=
1,2
,
…
are random variables with zero means. It is important, now, to find a special basis {
ϕ
k
} that makes corresponding a
k uncorrelated: E
[
a
i
a
j
]
=
0 for all i
≠
j. Denoting the covariance function of a
(
x
,
ω
) by R
(
x
,
y
)
=
E
[
a
(
x
,
ω
)
a
(
y
,
ω
)
], the basis functions {
ϕ
k
} should satisfy
(13)
E
[
a
i
a
j
]
=
∫
D
ϕ
i
(
x
)
d
x
∫
D
R
(
x
,
y
)
ϕ
j
(
y
)
d
y
=
0
,
i
≠
j
.Completion and orthonormality of{
ϕ
k
} in L
2
(
D
) follow that ϕ
k
(
x
) are eigenfunctions of R
(
x
,
y
):
(14)
∫
D
R
(
x
,
y
)
ϕ
j
(
y
)
d
y
=
λ
j
ϕ
j
(
x
)
,
j
=
1,2
,
…
,
where λ
j
=
E
[
a
j
2
]
>
0. Indeed, by choosing basis functions ϕ
k
(
x
) as the solutions of the eigenproblem (14), the random variables a
k
(
ω
) will be uncorrelated. In (14), denoting θ
k
=
a
k
/
λ
k, yields the following expansion:
(15)
a
(
x
,
ω
)
=
a
¯
(
x
)
+
∑
k
=
1
∞
λ
k
θ
k
(
ω
)
ϕ
k
(
x
)
,
where θ
k satisfy E
[
θ
k
]
=
0 and E
[
θ
i
θ
j
]
=
δ
i
j. In the case that a
(
x
,
ω
) is considered as a Gaussian process, a
k
(
ω
), k
=
1,2
,
… will be independent Gaussian random variables. The expansion (15) is known as the Karhunen-Loeve expansion (KLE) of the stochastic process a
(
x
,
ω
).Using the KLE (15), the stochastic process can be represented as a series of uncorrelated random variables. Since the basis functions ϕ
k
(
x
) are deterministic, the spatial dependence of the random process can be resolved by them. The convergence property of the KLE to the random process a
(
x
,
ω
) can be represented in the mean square sense
(16)
lim
⁡
N
→
∞
⁡
∫
D
E
|
a
(
x
,
ω
)
-
a
N
(
x
,
ω
)
|
2
d
x
=
0
,
where
(17)
a
N
=
a
¯
(
x
)
+
∑
k
=
1
N
λ
k
θ
k
ϕ
k
is a finite term KLE [4, 5].


## 2.3. Polynomial Chaos Expansion (PCE)

There are problems, as the solution of a PDE with random inputs, that the covariance function of a random processu
(
x
,
ω
), x
∈
D is not known. The solution of such problems can be represented using a polynomial chaos expansion (PCE) given by
(18)
u
(
x
,
ω
)
=
∑
k
=
1
∞
u
k
(
x
)
Ψ
k
(
ξ
)
,
where the functions u
k
(
x
) are deterministic coefficients, ξ is a vector of orthonormal random variables, and Ψ
k
(
ξ
) are multidimensional orthogonal polynomials that satisfy in the following properties:
(19)
〈
Ψ
1
〉
≡
E
[
Ψ
1
]
=
1
,
〈
Ψ
k
〉
=
0
,
k
>
1
,
〈
Ψ
i
Ψ
j
〉
=
h
i
δ
i
j
.The convergence property of PCE for a random quantity inL
2 is ensured by the Cameron-Martin theorem [6, 9]; that is,
(20)
〈
u
(
x
,
ω
)
-
∑
k
=
1
∞
u
k
(
x
)
Ψ
k
(
ξ
)
〉
⟶
L
2
0
.Hence, this convergence justifies a truncation of PCE to a finite number of terms,(21)
u
(
x
,
ω
)
=
∑
k
=
1
P
u
k
(
x
)
Ψ
k
(
ξ
)
,
where the value of P is determined by the highest degree of polynomial K, used to represent u, and the number N of random variables—the length of ξ—with the formula P
+
1
=
(
N
+
K
)
!
/
N
!
K
!. Generally, the value of N is the same as the number of uncorrelated random variables in the system or equivalently the truncation length of the truncated KLE. Typically, the value of K is chosen by some heuristic method. Indeed, in the case of K
=
1 and N random variables, the KLE is a special case of the PCE [9, 10].


## 2.4. Stochastic Galerkin

SupposeW
i
⊂
L
ρ
i
2
(
Ω
i
) with dimension ρ
i for i
=
1
,
…
,
M and V
⊂
H
0
1
(
D
) with dimension N. In addition, let {
ψ
n
i
}
n
=
1
p
i for i
=
1
,
…
,
M be basis of W
i and {
ϕ
i
}
i
=
1
N a basis of V. The finite dimensional tensor product space W
1
⨂
⋯
⨂
W
M
⨂
V can be defined as the space spanned by the functions {
ψ
n
1
1
,
…
,
ψ
n
M
M
ϕ
i
} for all n
1
∈
{
1
,
…
,
p
1
}
,
…
,
n
M
∈
{
1
,
…
,
p
M
} and i
∈
{
1
,
…
,
N
}. For simplification, let n denote a multi-index whose kth component n
k
∈
{
1
,
…
,
p
k
} and I denote the set of such multi-indices; then, the basis functions for the tensor product space have the form
(22)
v
n
i
(
x
,
ω
)
=
ϕ
i
(
x
)
Ψ
n
(
ω
)
,
where
(23)
Ψ
n
=
∏
k
=
1
M
ψ
n
k
k
(
ω
k
)
.Then, the Galerkin method is looking for a solutionu
^
∈
W
1
⨂
⋯
⨂
W
M
⨂
V such that for x
=
(
x
,
y
), all n
∈
I and i
=
1
,
…
,
N,
(24)
∫
Ω
ρ
(
ω
)
∫
D
a
(
x
,
ω
)
∇
u
h
P
(
x
,
ω
)
·
∇
v
n
i
(
x
,
ω
)
d
x
d
ω
=
∫
Ω
ρ
(
ω
)
∫
D
z
(
x
,
ω
)
v
n
i
(
x
,
ω
)
d
x
d
ω
,
whereρ
(
ω
) is the integration weight. Equation (24) is equivalent to the following equation, which results by plugging in the explicit formula for the basis functions
(25)
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
∫
D
a
(
x
,
ω
)
∇
u
h
P
(
x
,
ω
)
·
∇
ϕ
i
(
x
)
d
x
d
ω
=
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
∫
D
z
(
x
,
ω
)
ϕ
i
(
x
)
d
x
d
ω
for all n
∈
I and i
=
1
,
…
,
N.Now, considering the approximate solution as(26)
u
h
P
(
x
,
ω
)
=
∑
j
=
1
N
∑
m
∈
I
u
j
m
ϕ
j
(
x
)
Ψ
m
(
ω
)
and substituting (26) into (25) yields
(27)
∑
j
=
1
N
∑
m
∈
I
u
j
m
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
×
∫
D
a
(
x
,
ω
)
∇
ϕ
j
(
x
)
·
∇
ϕ
i
(
x
)
d
x
d
ω
=
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
∫
D
z
(
x
,
ω
)
ϕ
i
(
x
)
d
x
d
ω
for all m
∈
I and i
=
1
,
…
,
N. Let P
=
|
I
|
=
∏
k
=
1
M
p
k; then, a natural bijection between {
1
,
…
,
P
} and Ican be defined. Thus, (27) can be written as
(28)
∑
m
=
1
P
∑
j
=
1
N
u
m
j
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
×
∫
D
a
(
x
,
ω
)
∇
ϕ
j
(
x
)
·
∇
ϕ
i
(
x
)
d
x
d
ω
=
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
∫
D
z
(
x
,
ω
)
ϕ
i
(
x
)
d
x
d
ω
.Now, ifa and z have the following KLE
(29)
a
(
x
,
ω
)
=
a
0
(
x
)
+
∑
i
=
1
M
ω
i
a
i
(
x
)
,
z
(
x
,
ω
)
=
z
0
(
x
)
+
∑
i
=
1
M
ω
i
z
i
(
x
)
,then, we can rewrite the Galerkin projection (28) as
(30)
∑
m
=
1
P
∑
j
=
1
N
u
m
j
(
∑
k
=
1
M
(
K
k
)
i
,
j
(
K
0
)
i
,
j
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
d
ω
+
∑
k
=
1
M
(
K
k
)
i
,
j
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
ω
k
d
ω
)
=
(
z
0
)
i
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
d
ω
+
∑
k
=
1
M
(
z
k
)
i
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
ω
k
d
ω
,
where
(31)
(
K
k
)
i
,
j
=
∫
D
a
k
(
x
)
∇
ϕ
i
(
x
)
·
∇
ϕ
j
(
x
)
d
x
≅
∑
r
w
r
a
k
(
x
r
)
∇
ϕ
i
(
x
r
)
·
∇
ϕ
j
(
x
r
)
,
(
z
k
)
i
=
∫
D
z
i
(
x
)
ϕ
i
(
x
)
d
x
≅
∑
r
w
r
z
i
(
x
r
)
ϕ
i
(
x
r
)
for k
=
1
,
…
,
M  [3]. Equation (31) is approximated by quadrature rule in spatial domain. Now, assume that there exist functions Ψ
n, n
=
1
,
…
P, such that
(32)
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
d
ω
=
δ
m
n
,
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
ω
k
d
ω
=
C
k
n
δ
n
m
,
where, for m
=
n, the above integral is approximated using quadrature rule in random domain as follows:
(33)
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
n
(
ω
)
ω
k
d
ω
≅
∑
s
w
s
′
ρ
(
ω
s
)
Ψ
n
2
(
ω
s
)
ω
k
.
Then, (30) becomes
(34)
∑
m
=
1
P
∑
j
=
1
N
u
m
j
(
(
K
0
)
i
,
j
+
∑
k
=
1
M
C
k
n
(
K
k
)
i
,
j
)
δ
n
m
=
(
z
0
)
i
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
d
ω
+
∑
k
=
1
M
(
z
k
)
i
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
ω
k
d
ω
or, equivalently,
(35)
∑
j
=
1
N
u
m
j
(
(
K
0
)
i
,
j
+
∑
k
=
1
M
C
k
n
(
K
k
)
i
,
j
)
=
(
z
0
)
i
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
d
ω
+
∑
k
=
1
M
(
z
k
)
i
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
ω
k
d
ω
.In (34) and (35) we have
(36)
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
d
ω
=
0
;
for
n
≠
0
,
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
ω
k
d
ω
=
∫
Ω
ρ
(
ω
)
Ψ
0
(
ω
)
Ψ
n
(
ω
)
ω
k
d
ω
=
C
k
n
δ
0
n
=
C
k
0
.Finally, (35) can be considered as the following block diagonal system:
(37)
(
K
0
+
∑
k
=
1
M
C
k
1
K
k
⋯
0
⋮
⋱
⋮
0
⋯
K
0
+
∑
k
=
1
M
C
k
p
K
k
)
(
u
→
1
⋮
u
→
P
)
=
(
z
→
1
⋮
z
→
p
)
or briefly K
u
→
=
z
→.


## 3. Distributed Optimal Control Problem

Using the stochastic Galerkin method for general objective functions, computing the derivatives becomes very tedious and rather messy. Consider the following optimal control problem:(38)
min
⁡
u
∈
U
,
z
∈
𝕫
J
(
u
,
z
)
=
1
2
E
[
∥
u
(
x
,
ω
)
-
u
^
(
x
)
∥
L
2
(
D
)
2
]
min
⁡
u
∈
U
,
z
∈
𝕫
J
(
u
,
z
)
=
+
α
2
∥
z
(
x
)
∥
L
2
(
D
)
2
s
.
t
.
-
∇
·
(
a
(
x
,
ω
)
∇
u
(
x
,
ω
)
)
=
z
(
x
)
in
D
×
Ω
u
(
x
,
ω
)
=
0
on
⁡
∂
D
×
Ω
.The solutions to linear elliptic SPDEs live in the spaceU
=
H
0
1
(
D
)
⨂
L
2
(
Ω
) and the control space ℤ
=
L
2
(
D
). Denoting the constraint equation by e
(
u
,
z
)
=
0, usually, it is possible to invoke the implicit function theorem to find a solution u
(
z
) to the constraint equation. This leads to an implicitly defined objective function J
^
(
z
)
:
=
J
(
u
(
z
)
,
z
). Now we focus on computing derivatives of J
^ [11]. For any direction s
∈
ℤ,
(39)
〈
J
^
′
(
z
)
,
s
〉
ℤ
*
,
ℤ
=
〈
J
u
(
u
(
z
)
,
z
)
,
u
z
(
z
)
s
〉
U
*
,
U
+
〈
J
z
(
u
(
z
)
,
z
)
,
s
〉
ℤ
*
,
ℤ
,
where U
* and ℤ
* are the dual space of U and ℤ, respectively. Now, computing the derivative of the constraint e with respect to z, in the direction s, and applying it to the implicit solution to e
(
u
(
z
)
,
z
)
=
0 yield
(40)
e
u
(
u
(
z
)
,
z
)
u
z
(
z
)
s
+
e
z
(
u
(
z
)
,
z
)
s
=
0
.
So u
z
(
z
) is
(41)
u
z
(
z
)
s
=
-
e
u
(
u
(
z
)
,
z
)
-
1
e
z
(
u
(
z
)
,
z
)
s
.
Hence
(42)
〈
J
^
′
(
z
)
,
s
〉
ℤ
*
,
ℤ
=
-
〈
J
u
(
u
(
z
)
,
z
)
,
e
u
(
u
(
z
)
,
z
)
-
1
×
e
z
(
u
(
z
)
,
z
)
s
(
u
(
z
)
,
z
)
-
1
〉
U
*
,
U
+
〈
J
z
(
u
(
z
)
,
z
)
,
s
〉
ℤ
*
,
ℤ
.Applying the adjoint to theu
z
(
z
) term,
(43)
〈
J
^
′
(
z
)
,
s
〉
ℤ
*
,
ℤ
=
-
〈
e
z
(
u
(
z
)
,
z
)
*
e
u
(
u
(
z
)
,
z
)
-
*
×
J
u
(
u
(
z
)
,
z
)
,
s
(
u
(
z
)
,
z
)
-
*
〉
ℤ
*
,
ℤ
+
〈
J
z
(
u
(
z
)
,
z
)
,
s
〉
ℤ
*
,
ℤ
=
-
〈
e
z
(
u
(
z
)
,
z
)
*
e
u
(
u
(
z
)
,
z
)
-
*
J
u
(
u
(
z
)
,
z
)
+
J
z
(
u
(
z
)
,
z
)
,
s
(
u
(
z
)
,
z
)
-
*
〉
ℤ
*
,
ℤ
.
Considering the adjoint state, μ
=
μ
(
z
)
∈
ℤ
*, as the solution to
(44)
e
u
(
u
(
z
)
,
z
)
*
μ
=
-
J
u
(
u
(
z
)
,
z
)
the derivative of J
^
′ becomes
(45)
J
^
′
(
z
)
=
e
z
(
u
(
z
)
,
z
)
*
μ
(
z
)
+
J
z
(
u
(
z
)
,
z
)
.
The weak form of the constraint function is defined as follows: for all v
∈
U,
(46)
〈
e
(
u
,
z
)
,
v
〉
U
*
,
U
=
∫
Ω
ρ
(
ω
)
∫
D
(
a
(
x
,
ω
)
∇
u
(
x
,
ω
)
·
∇
v
(
x
,
ω
)
-
z
(
x
,
ω
)
v
(
x
,
ω
)
)
d
x
d
ω
.Thus, the constraint functione
:
U
×
ℤ
→
U
* and Y
=
U
*. Since the constraint function e acts linearly with respect to u and z, thus the corresponding derivative of e with respect to z in the direction δ
z
∈
ℤ, for all v
∈
U, is given by
(47)
〈
e
z
(
u
,
z
)
δ
z
,
v
〉
U
*
,
U
=
-
∫
Ω
ρ
(
ω
)
∫
D
δ
z
(
x
,
ω
)
v
(
x
,
ω
)
d
x
d
ω
.Similarly, the derivative ofe with respect to u, for any direction δ
u
∈
U and for all v
∈
U, is given by
(48)
〈
e
u
(
u
,
z
)
δ
u
,
v
〉
U
*
,
U
=
∫
Ω
ρ
(
ω
)
∫
D
a
(
x
,
ω
)
∇
δ
u
(
x
,
ω
)
·
∇
v
(
x
,
ω
)
d
x
d
ω
.Both of these derivatives are symmetric, soe
u
(
u
,
z
)
*
=
e
u
(
u
,
z
) and e
z
(
u
,
z
)
*
=
e
z
(
u
,
z
). From here, the adjoint can be computed as
(49)
e
u
(
u
(
z
)
,
z
)
μ
=
-
J
u
(
u
(
z
)
,
z
)
.Using this to any directionv
∈
U,
(50)
〈
e
u
(
u
(
z
)
,
z
)
μ
,
v
〉
U
*
,
U
=
∫
Ω
ρ
(
ω
)
∫
D
a
(
x
,
ω
)
∇
μ
(
x
,
ω
)
·
∇
v
(
x
,
ω
)
d
x
d
ω
=
-
J
u
(
u
(
z
)
,
z
)
v
.By the chain rule, the derivative ofJ with respect to u in the direction δ
u is
(51)
〈
J
u
(
u
,
z
)
,
δ
u
〉
U
*
,
U
=
E
[
∥
u
(
x
,
ω
)
-
u
^
(
x
)
∥
L
2
(
D
)
2
×
〈
u
(
x
,
ω
)

-

u
^
(
x
)
,
δ
u
(
x
,
ω
)
〉
L
2
(
D
)

∥
u
(
x
,
ω
)
-
u
^
(
x
)
∥
L
2
(
D
)
2
]
.Similarly, the derivative ofJ with respect to z in the δ
z
∈
ℤ direction is
(52)
〈
J
z
(
u
,
z
)
,
δ
z
〉
ℤ
*
,
ℤ
=
α
〈
z
,
δ
z
〉
ℤ
.With these computations, one can compute the first derivative ofJ
^
(
z
) via the adjoint approach. The aim is to derive the discrete versions of the objective function, gradient, and Hessian times a vector calculation corresponding to the stochastic Galerkin solution technique for SPDE; that is, e
(
u
,
z
)
=
0. Suppose X
h
⊂
H
0
1
(
D
) is a finite dimensional subspace of dimension N and Y
h
⊂
L
2
(
Ω
) is a finite dimensional subspace of dimension P; then X
h
⨂
Y
h is a finite dimensional subspace of the state space U. Similarly, let ℤ
h
⊂
ℤ be a finite dimensional subspace of the control space (with dimension). For the stochastic Galerkin method, it is assumed that Y
h
=
ℙ
P
-
1, the space of polynomials with highest degree P
-
1, and X
h is any finite element space; here, X
h is the space of linear functions built on a given mesh 𝒯
h. We choose the system of polynomials that are ρ-orthonormal to be a basis for Y
h; that is, Y
h
=
span
⁡
{
Ψ
1
(
ω
)
,
…
,
Ψ
P
(
ω
)
}, where
(53)
∫
Ω
ρ
(
ω
)
Ψ
n
(
ω
)
Ψ
m
(
ω
)
d
ω
=
δ
m
n
.The discretized optimization problem in the stochastic Galerkin framework is(54)
min
⁡
z
h
∈
ℤ
h
⁡
J
^
h
P
(
z
h
)
=
1
2
E
[
∥
∑
k
=
1
P
∑
j
=
1
N
(
u
→
k
(
z
h
)
)
j
Ψ
k
(
ω
)
ϕ
j
(
x
)
-
u
^
(
x
)
∑
k
=
1
P
∑
j
=
1
N
(
u
→
k
(
z
h
)
)
j
∥
L
2
(
D
)
2
]
+
α
2
∥
z
h
∥
L
2
(
D
)
2
,
where
(55)
u
→
k
(
z
h
)
=
(
u
→
1
(
z
h
)
T
,
…
,
u
→
P
(
z
h
)
T
)
T
=
u
→
is the stochastic Galerkin solution to the state equation
(56)
(
K
11
⋯
K
1
P
⋮
⋱
⋮
K
P
1
⋯
K
P
P
)
︸
=
K
(
u
→
1
⋮
u
→
P
)
︸
=
u
→
=
(
E
[
Ψ
1
]
M
z
→
h
⋮
E
[
Ψ
P
]
M
z
→
h
)
=
(
δ
11
M
z
→
h
⋮
δ
1
P
M
z
→
h
)
since E
[
Ψ
k
]
=
E
[
1
Ψ
k
]
=
E
[
Ψ
1
Ψ
k
]
=
δ
1
k. The blocks of the above K matrix have the form
(57)
(
K
k
l
)
i
j
=
∫
Ω
ρ
(
ω
)
Ψ
k
(
ω
)
Ψ
l
(
ω
)
×
∫
D
a
(
x
,
ω
)
∇
ϕ
i
(
x
)
·
∇
ϕ
j
(
x
)
d
x
d
ω
.First, in order to compute the derivatives of the objective function, we attempt to compute the adjoint state [8]. Indeed, the adjoint state, in the infinite dimensional formulation, solves the following equations:
(58)
-
∇
·
(
a
(
x
,
ω
)
∇
μ
(
x
,
ω
)
)
=
-
(
u
(
x
,
ω
)
-
u
^
(
x
)
)
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
in
D
×
Ω
,
μ
(
x
,
ω
)
=
0
on
⁡
∂
D
×
Ω
.Thus, as in (37), the block system for (58) can be written in the form K
μ
→
=
F, where F
=
(
F
1
T
,
…
,
F
p
T
)
T and F
k is defined as
(59)
(
F
→
k
)
i
=
-
∫
Ω
ρ
(
ω
)
Ψ
k
(
ω
)
∫
D
(
u
(
x
,
ω
)
-
u
^
(
x
)
)
ϕ
i
(
x
)
d
x
d
ω
=
-
∫
Ω
ρ
(
ω
)
Ψ
k
(
ω
)
∫
D
(
∑
l
=
1
P
∑
j
=
1
N
(
u
→
l
)
j
ϕ
j
(
x
)
Ψ
l
(
ω
)
ϕ
i
(
x
)
-
u
^
(
x
)
ϕ
i
(
x
)
∑
j
=
1
N
(
u
→
l
)
j
)
d
x
d
ω
=
-
∑
l
=
1
P
∑
j
=
1
N
(
u
→
l
)
j
∫
Ω
ρ
(
ω
)
Ψ
k
(
ω
)
Ψ
l
(
ω
)
d
ω
×
∫
D
ϕ
j
(
x
)
ϕ
i
(
x
)
d
x
+
∫
Ω
ρ
(
ω
)
Ψ
k
(
ω
)
d
ω
∫
D
u
^
(
x
)
ϕ
i
(
x
)
d
x
=
b
→
i
δ
1
k
-
∑
l
=
1
P
∑
j
=
1
N
M
i
j
(
u
→
l
)
j
δ
k
l
=
b
→
i
δ
1
k
-
∑
j
=
1
N
M
i
j
(
u
→
l
)
j
,
in which b
→
i
=
{
∫
u
^
ϕ
i
} and M
i
j
=
{
∫
ϕ
i
ϕ
j
}.Thus,(60)
(
F
→
k
)
i
=
{
b
→
-
M
u
1
if
k
=
1
-
M
u
if
k
=
2
,
…
,
P
,
where b
=
[
b
→
1
,
…
,
b
→
N
] and the matrix M
=
(
M
i
j
) is of order N
×
N.Now, the derivative ofJ
^
(
z
) in the direction ϕ
j
(
x
) for all j
=
1
,
…
,
N, with the computed adjoint state, is
(61)
〈
J
^
′
(
z
)
,
ϕ
j
(
x
)
〉
ℤ
*
,
ℤ
=
〈
(
u
(
z
,
·
,
·
)
,
z
)
ϕ
j
e
z
(
u
(
z
,
·
,
·
)
,
z
)
*
μ
+
J
z
(
u
(
z
,
·
,
·
)
,
z
)
,
ϕ
j
(
x
)
(
u
(
z
,
·
,
·
)
,
z
)
*
〉
ℤ
*
,
ℤ
=
-
∫
Ω
ρ
(
ω
)
∫
D
μ
(
z
,
·
,
·
)
ϕ
j
(
x
)
d
x
d
ω
+
α
∫
D
z
(
x
)
ϕ
j
(
x
)
d
x
=
-
∑
k
=
1
P
∑
i
=
1
N
(
μ
→
k
)
i
∫
Ω
ρ
(
ω
)
Ψ
k
(
ω
)
d
ω
×
∫
D
ϕ
j
(
x
)
ϕ
i
(
x
)
d
x
+
α
∑
i
=
1
N
z
→
i
∫
D
ϕ
j
(
x
)
ϕ
i
(
x
)
d
x
=
-
∑
k
=
1
P
∑
i
=
1
N
(
μ
→
k
)
i
M
i
j
δ
1
k
+
α
∑
i
=
1
N
M
i
j
z
→
i
=
-
∑
i
=
1
N
(
μ
→
1
)
i
M
i
j
+
α
∑
i
=
1
N
M
i
j
z
→
i
.Therefore, the gradient can be computed as follows:(62)
J
^
′
(
z
)
=
M
(
α
z
→
h
-
μ
→
1
)
.Now, for any vectorv
∈
ℤ, it is possible to write v as
(63)
v
=
∑
k
=
1
N
v
k
ϕ
k
(
x
)
.In order to compute the Hessian times a vector, that is, multiply the Hessian of theJ
^ with a some vector v
∈
ℤ, the equation e
u
(
u
(
z
)
,
z
)
w
=
e
z
(
u
(
z
)
,
z
)
v for w must be solved. Similar to the adjoint computation, w solves the linear system K
w
→
=
G where G
→
=
(
G
→
1
T
,
…
,
G
→
P
T
) and G
→
k for k
=
1
,
…
,
P is
(64)
(
G
→
k
)
i
=
∫
Ω
ρ
(
ω
)
Ψ
k
(
ω
)
d
ω
∫
D
∑
j
=
1
N
v
j
ϕ
j
(
x
)
ϕ
i
(
x
)
d
x
=
∑
j
=
1
N
M
i
j
v
j
δ
1
k
for all i
=
1
,
…
,
N. Thus,
(65)
G
→
k
=
{
M
v
→
if
k
=
1
0
if
k
=
2
,
…
,
P
.Now, solvinge
u
(
u
(
z
)
,
z
)
q
=
J
u
u
(
u
(
z
)
,
z
)
w for q requires the solution to the linear system K
q
→
=
H, where
(66)
(
H
→
k
)
i
=
∑
l
=
1
P
∑
j
=
1
N
(
w
→
l
)
j
∫
Ω
ρ
(
ω
)
Ψ
k
(
ω
)
Ψ
l
(
ω
)
d
ω
×
∫
D
ϕ
j
(
x
)
ϕ
i
(
x
)
d
x
=
∑
j
=
1
N
M
i
j
(
w
→
k
)
j
or equivalently H
→
k
=
M
w
→
k for all k
=
1
,
…
,
P. Hence, in the direction ϕ
i for i
=
1
,
…
,
N, J
^
′′
(
z
)
v can be approximated by
(67)
〈
J
^
′′
(
z
)
v
,
ϕ
i
〉
ℤ
*
,
ℤ
≈
〈
J
^
h
p
′′
(
z
h
)
v
,
ϕ
i
〉
ℤ
*
,
ℤ
=
-
∫
Ω
ρ
(
ω
)
∫
D
∑
k
=
1
P
∑
j
=
1
N
(
q
→
k
)
j
ϕ
j
(
x
)
Ψ
k
(
ω
)
×
ϕ
i
(
x
)
d
x
d
ω
+
α
∫
D
v
(
x
)
ϕ
i
(
x
)
d
x
=
-
∑
k
=
1
P
∑
j
=
1
N
M
i
j
(
q
→
k
)
j
δ
1
k
+
α
∑
j
=
1
N
M
i
j
v
→
j
=
∑
j
=
1
N
M
i
j
(
α
v
→
j
-
(
q
→
1
)
j
)
.Therefore,(68)
J
^
′′
(
z
)
v
≈
M
(
α
v
→
-
q
→
1
)
.Having (62) and (68), we can use the iterative solvers to the Newton equation
(69)
J
^
′′
(
z
k
)
s
k
=
-
J
^
′
(
z
k
)
.Now, by using the preconditioned conjugate gradient (PCG) method, the Newton equation (69) is solved approximately. When we access to the sufficiently small residual of Newton system, the PCG method is truncated. In real world computation, it is possible to employ some globalization technique for Newton’s method. Considering the case of implementation and relatively low computational cost, line search techniques are popular choices in this way. Finding an optimal step size α
k and using this step size to generate the iterate z
k
+
1
=
z
k
+
α
k
S
k, are the mission of a line search algorithm. The Armijo condition (or sufficient decrease condition) that the step size is required to satisfy is
(70)
J
(
z
k
+
α
k
S
k
)
≤
J
(
z
k
)
+
c
α
k
J
′
(
z
k
)
S
k
,
where c
∈
(
0,1
) and typically is quite small, for example, c
=
10
-
4 [12, 13].


### 3.1. Kronecker Product Preconditioners

Let us to rewrite (37) in the form A
^
u
→
=
z
→, where
(71)
A
^
=
(
∑
k
=
0
M
G
k
⨂
K
k
)
.The strategy here is to introduceP
^
=
G
⨂
K
0 as a preconditioner, such that
(72)
G
=
argmin
{
H
∈
ℝ
P
×
P
:
∥
A
^
-
H
⨂
K
0
∥
F
}
,
where as mentioned in PCE, P is equal with the degree of polynomial chaos truncation, and ∥
·
∥
F denotes the Frobenius norm. The closed form of the solution can be written as follow [14]:
(73)
G
=
I
+
∑
P
tr
⁡
(
K
k
T
K
0
)
tr
⁡
(
K
0
T
K
0
)
G
k
.Here,tr
⁡
(
K
k
T
K
0
)
=
∑
i
=
1
N
q
[
K
k
]
i
,
i
[
K
0
]
i
,
i and hence, the coefficients in the above equality can be computed straightforward. In addition, since A
^ and K
0 are symmetric and positive definite, thus G and P
^
=
G
⨂
K
0 have also these properties [15, 16].Example 1.
Consider the following distributed optimal control problem(74)
min
⁡
u
,
z
⁡
J
(
u
,
z
)
=
E
[
1
2
∥
u
(
x
,
ω
)
-
u
^
(
x
)
∥
L
2
(
D
)
2
]
+
α
2
∥
z
(
x
)
∥
L
2
(
D
)
2
,
-
∇
·
(
a
(
x
,
ω
)
∇
u
(
x
,
ω
)
)
=
z
(
x
)
,
x
∈
D
,
ω
∈
Ω
,
where D
=
[
-
1,1
]
2, Ω
~
U
[
-
1,1
], the boundary conditions are given as
(75)
u
(
-
1
,
y
,
ω
)
=
u
(
1
,
y
,
ω
)
=
u
(
x
,
1
,
ω
)
=
u
(
y
,
-
1
,
ω
)
=
0
,
and the desired function is given by
(76)
u
^
(
x
,
y
)
=
exp
⁡
[
-
64
(
(
x
-
1
2
)
2
+
(
y
-
1
2
)
2
)
]
.The random fielda
(
x
,
ω
) is characterized by its mean and covariance function
(77)
E
[
a
]
=
a
-
=
10
,
R
(
x
1
,
x
2
)
=
e
-
|
x
1
-
x
2
|
,
x
1
,
x
2
∈
[
-
1,1
]
.Following Section2.2, the truncated KLE can be represented as:
(78)
a
(
x
,
y
,
ω
)
=
10
+
∑
j
=
0
N
λ
j
θ
j
ϕ
j
(
x
)
.The aim is to calculateu
opt and z
opt such that for all u
∈
H
0
1
(
D
)
⨂
L
2
(
Ω
) and all z
∈
L
2
(
D
)
(79)
J
(
u
opt
,
z
opt
)
≤
J
(
u
,
z
)
.The eigenpairs{
λ
j
,
θ
j
} in truncated KLE solve the integral equation
(80)
∫
-
1
1
e
-
|
x
1
-
x
2
|
ϕ
j
(
x
2
)
d
x
2
=
λ
j
ϕ
j
(
x
1
)
.For this special case of the covariance function, we have explicit expressions forϕ
j and λ
j, [6]. Let ω
j
even and ω
j
odd solve the equations
(81)
1
-
ω
j
even
tan
⁡
(
ω
j
even
)
=
0
,
ω
j
odd
+
tan
⁡
(
ω
j
odd
)
=
0
.Then the even and odd indexed eigenfunctions are given by(82)
ϕ
j
even
(
x
)
=
cos
⁡
⁡
(
ω
j
even
x
)
1
+
(
sin
⁡
(
2
ω
j
even
)
/
2
ω
j
even
)
ϕ
j
odd
(
x
)
=
sin
⁡
(
ω
j
odd
x
)
1
-
(
sin
⁡
(
2
ω
j
odd
)
/
2
ω
j
odd
)
and corresponding eigenvalues
(83)
λ
j
even
(
x
)
=
2
ω
j
even
2
+
1
λ
j
odd
(
x
)
=
2
ω
j
odd
2
+
1
.We chooseθ
=
(
θ
1
,
…
,
θ
N
)
T to be independent random variables uniformly distributed over the interval [
-
1,1
].First of all, we find numerical approximations ofω
j
even and ω
j
oddwith bisection method, and then, with the eigenpairs evaluated by this ω
j
even  and ω
j
odd, by choosing dimension N
=
3 and order K
=
5, which emphasize that P
=
55, we construct the KLE. Taking α
=
0.00005, Algorithm 1 is used to compute u
opt and z
opt. Illustration of computed optimal control and corresponding state in addition to the global matrix sparsity, standard deviation and expected values of the solution, adjoint values, and initial state approximation of the desired function is shown in Figure 1.Algorithm 1
(1) Givenz
0 and tol
>
0, set k
=
0
(2) ComputeK, M, A, u
→, and b
→
(3) Compute(
F
→
k
)
i
=
{
b
-
M
u
if
k
=
1
-
M
u
if
k
=
2
,
…
,
P
.
(4) SolveK
μ
→
=
F
(5) ComputeJ
^
′
(
z
k
)
=
M
(
α
z
k
-
μ
→
1
).
(6) ComputeG
→
1
=
M
v
→ and G
→
t
=
0 for t
=
2
,
…
,
P
(7) SolveK
w
→
=
G
=
(
G
→
1
T
,
…
,
G
→
P
T
)
T
(8) ComputeH
→
k
=
M
w
→
k for k
=
1
,
…
,
P
(9) SolveK
q
→
=
H
=
(
H
→
1
T
,
…
,
H
→
P
T
)
T
(10) ComputeJ
^
′
′
(
z
)
v
≈
M
(
α
v
→
-
q
→
1
)
(11) If∥
∇
J
^
(
z
k
)
∥
<
tol, stop.
(12) SolveJ
^
′
′
(
z
k
)
s
k
=
-
J
^
′
(
z
k
) using PCG method with preconditioner P
^
=
G
⨂
K
0:
(i) Setx
0
=
0.
(ii) Setr
0
←
J
^
′
′
(
z
k
)
x
0
+
J
^
′
(
z
k
);
(iii) SolveP
^
y
0
=
r
0 for y
0;
(iv) Setp
0
=
-
r
0, t
←
0;
(v) While∥
r
t
∥
>
tol
α
k
←
r
t
T
y
t
p
t
T
A
p
t
x
t
+
1
←
x
t
+
α
t
p
t
r
t
+
1
←
r
t
+
α
t
A
p
t
P
y
t
+
1
←
r
t
+
1
β
t
+
1
←
r
t
+
1
T
y
t
+
1
r
t
T
y
t
p
t
+
1
←
-
y
t
+
1
+
β
t
+
1
p
t
t
←
t
+
1
End (while)
(vi)s
k
=
x
t
+
1
(13) Perform Armijo line-search.
(i) Setα
k
=
1 and evaluate J
(
z
k
+
α
k
S
k
)
(ii) WhileJ
(
z
k
+
α
k
S
k
)
>
J
(
z
k
)
+
10
-
4
α
k
S
k
∇
J
^
(
z
k
) do
(i) Setα
k
=
α
k
/
2 and evaluate J
(
z
k
+
α
k
S
k
).
(14) Setz
k
+
1
=
z
k
+
α
k
S
k, k
←
k
+
1. Go to (2).Illustration of the optimal control and state achieved by solving Example1, taking u
^
(
x
,
y
)
=
exp
⁡
⁡
[
-
64
(
(
x
-
(
1
/
2
)
)
2
+
(
y
-
(
1
/
2
)
)
2
)
], N
=
3, K
=
5, and α
=
0.00005.
(a)
Desired values
(b)
Initial state
(c)
Approximated solution (state)
(d)
Control
(e)
Adjoint
(f)
Expected value
(g)
Standard dev.
(h)
Sparsity of global matrixn
z
=
1815980In Figure1(a), the desired function u
^
(
x
,
y
) is plotted for x
i
=
y
i
=
-
1
+
0.04
*
i,  i
=
0,1
,
2
,
…
,
50. Initial state is computed and plotted in Figure 1(b) by solving (37) where z
→
=
0. Approximated u
opt corresponding to the optimal control z
opt is calculated by Algorithm 1 and its graph is depicted in Figure 1(c). As it was expected, by using the proposed Algorithm 1, even by very far initial state, the approximated state solution reaches to the desired function. In Figure 1(d), the z
opt is depicted. Serious changes during the computation for initial z
→
=
0 happened reach to the optimal control function z
opt. The adjoint state is computed and plotted in Figure 1(e). As it is expected from (58), the adjoint state corresponding to the z
o
p
t must approximate the initial state (see Figures 1(b) and 1(c)). Figures 1(f) and 1(g) represent counter plots of the optimal state expectation and standard deviation. Finally, in Figure 1(h), the representation of coefficient matrix sparsity with the number of nonzero component in the computation is plotted in x
-
y plane.In Example1 we considered an exponential desired function, with v-sharp points. In the next example the smooth desired function and boundary conditions are considered. This is the case that Algorithm 1 is more efficient to use.Example 2.
Consider the optimal control of Example1, in which the right hand side function as well as boundary conditions is replaced by u
^
(
x
,
y
)
=
2
π
2
sin
⁡
(
π
x
)
sin
⁡
(
π
y
), N
=
3, and K
=
4, which emphasize that P
=
34. Taking α
=
0.00001, the illustration of computed optimal control and corresponding state in addition to the global matrix sparsity, standard deviation and expected values of the solution, adjoint values and initial state approximation of the desired function is shown in Figure 2.Illustration of optimal control and state achieved by solving Example2, taking u
^
(
x
,
y
)
=
2
π
2
sin
⁡
(
π
x
)
sin
⁡
(
π
y
), N
=
3, K
=
4, and α
=
0.00001.
(a)
Desired values
(b)
Mean values
(c)
Approximated solution (State)
(d)
Control
(e)
Adjoint
(f)
Expected value-SGS
(g)
Standard dev.-SGS
(h)
Sparsity of global matrixn
z
=
3046160


## 3.1. Kronecker Product Preconditioners

Let us to rewrite (37) in the form A
^
u
→
=
z
→, where
(71)
A
^
=
(
∑
k
=
0
M
G
k
⨂
K
k
)
.The strategy here is to introduceP
^
=
G
⨂
K
0 as a preconditioner, such that
(72)
G
=
argmin
{
H
∈
ℝ
P
×
P
:
∥
A
^
-
H
⨂
K
0
∥
F
}
,
where as mentioned in PCE, P is equal with the degree of polynomial chaos truncation, and ∥
·
∥
F denotes the Frobenius norm. The closed form of the solution can be written as follow [14]:
(73)
G
=
I
+
∑
P
tr
⁡
(
K
k
T
K
0
)
tr
⁡
(
K
0
T
K
0
)
G
k
.Here,tr
⁡
(
K
k
T
K
0
)
=
∑
i
=
1
N
q
[
K
k
]
i
,
i
[
K
0
]
i
,
i and hence, the coefficients in the above equality can be computed straightforward. In addition, since A
^ and K
0 are symmetric and positive definite, thus G and P
^
=
G
⨂
K
0 have also these properties [15, 16].Example 1.
Consider the following distributed optimal control problem(74)
min
⁡
u
,
z
⁡
J
(
u
,
z
)
=
E
[
1
2
∥
u
(
x
,
ω
)
-
u
^
(
x
)
∥
L
2
(
D
)
2
]
+
α
2
∥
z
(
x
)
∥
L
2
(
D
)
2
,
-
∇
·
(
a
(
x
,
ω
)
∇
u
(
x
,
ω
)
)
=
z
(
x
)
,
x
∈
D
,
ω
∈
Ω
,
where D
=
[
-
1,1
]
2, Ω
~
U
[
-
1,1
], the boundary conditions are given as
(75)
u
(
-
1
,
y
,
ω
)
=
u
(
1
,
y
,
ω
)
=
u
(
x
,
1
,
ω
)
=
u
(
y
,
-
1
,
ω
)
=
0
,
and the desired function is given by
(76)
u
^
(
x
,
y
)
=
exp
⁡
[
-
64
(
(
x
-
1
2
)
2
+
(
y
-
1
2
)
2
)
]
.The random fielda
(
x
,
ω
) is characterized by its mean and covariance function
(77)
E
[
a
]
=
a
-
=
10
,
R
(
x
1
,
x
2
)
=
e
-
|
x
1
-
x
2
|
,
x
1
,
x
2
∈
[
-
1,1
]
.Following Section2.2, the truncated KLE can be represented as:
(78)
a
(
x
,
y
,
ω
)
=
10
+
∑
j
=
0
N
λ
j
θ
j
ϕ
j
(
x
)
.The aim is to calculateu
opt and z
opt such that for all u
∈
H
0
1
(
D
)
⨂
L
2
(
Ω
) and all z
∈
L
2
(
D
)
(79)
J
(
u
opt
,
z
opt
)
≤
J
(
u
,
z
)
.The eigenpairs{
λ
j
,
θ
j
} in truncated KLE solve the integral equation
(80)
∫
-
1
1
e
-
|
x
1
-
x
2
|
ϕ
j
(
x
2
)
d
x
2
=
λ
j
ϕ
j
(
x
1
)
.For this special case of the covariance function, we have explicit expressions forϕ
j and λ
j, [6]. Let ω
j
even and ω
j
odd solve the equations
(81)
1
-
ω
j
even
tan
⁡
(
ω
j
even
)
=
0
,
ω
j
odd
+
tan
⁡
(
ω
j
odd
)
=
0
.Then the even and odd indexed eigenfunctions are given by(82)
ϕ
j
even
(
x
)
=
cos
⁡
⁡
(
ω
j
even
x
)
1
+
(
sin
⁡
(
2
ω
j
even
)
/
2
ω
j
even
)
ϕ
j
odd
(
x
)
=
sin
⁡
(
ω
j
odd
x
)
1
-
(
sin
⁡
(
2
ω
j
odd
)
/
2
ω
j
odd
)
and corresponding eigenvalues
(83)
λ
j
even
(
x
)
=
2
ω
j
even
2
+
1
λ
j
odd
(
x
)
=
2
ω
j
odd
2
+
1
.We chooseθ
=
(
θ
1
,
…
,
θ
N
)
T to be independent random variables uniformly distributed over the interval [
-
1,1
].First of all, we find numerical approximations ofω
j
even and ω
j
oddwith bisection method, and then, with the eigenpairs evaluated by this ω
j
even  and ω
j
odd, by choosing dimension N
=
3 and order K
=
5, which emphasize that P
=
55, we construct the KLE. Taking α
=
0.00005, Algorithm 1 is used to compute u
opt and z
opt. Illustration of computed optimal control and corresponding state in addition to the global matrix sparsity, standard deviation and expected values of the solution, adjoint values, and initial state approximation of the desired function is shown in Figure 1.Algorithm 1
(1) Givenz
0 and tol
>
0, set k
=
0
(2) ComputeK, M, A, u
→, and b
→
(3) Compute(
F
→
k
)
i
=
{
b
-
M
u
if
k
=
1
-
M
u
if
k
=
2
,
…
,
P
.
(4) SolveK
μ
→
=
F
(5) ComputeJ
^
′
(
z
k
)
=
M
(
α
z
k
-
μ
→
1
).
(6) ComputeG
→
1
=
M
v
→ and G
→
t
=
0 for t
=
2
,
…
,
P
(7) SolveK
w
→
=
G
=
(
G
→
1
T
,
…
,
G
→
P
T
)
T
(8) ComputeH
→
k
=
M
w
→
k for k
=
1
,
…
,
P
(9) SolveK
q
→
=
H
=
(
H
→
1
T
,
…
,
H
→
P
T
)
T
(10) ComputeJ
^
′
′
(
z
)
v
≈
M
(
α
v
→
-
q
→
1
)
(11) If∥
∇
J
^
(
z
k
)
∥
<
tol, stop.
(12) SolveJ
^
′
′
(
z
k
)
s
k
=
-
J
^
′
(
z
k
) using PCG method with preconditioner P
^
=
G
⨂
K
0:
(i) Setx
0
=
0.
(ii) Setr
0
←
J
^
′
′
(
z
k
)
x
0
+
J
^
′
(
z
k
);
(iii) SolveP
^
y
0
=
r
0 for y
0;
(iv) Setp
0
=
-
r
0, t
←
0;
(v) While∥
r
t
∥
>
tol
α
k
←
r
t
T
y
t
p
t
T
A
p
t
x
t
+
1
←
x
t
+
α
t
p
t
r
t
+
1
←
r
t
+
α
t
A
p
t
P
y
t
+
1
←
r
t
+
1
β
t
+
1
←
r
t
+
1
T
y
t
+
1
r
t
T
y
t
p
t
+
1
←
-
y
t
+
1
+
β
t
+
1
p
t
t
←
t
+
1
End (while)
(vi)s
k
=
x
t
+
1
(13) Perform Armijo line-search.
(i) Setα
k
=
1 and evaluate J
(
z
k
+
α
k
S
k
)
(ii) WhileJ
(
z
k
+
α
k
S
k
)
>
J
(
z
k
)
+
10
-
4
α
k
S
k
∇
J
^
(
z
k
) do
(i) Setα
k
=
α
k
/
2 and evaluate J
(
z
k
+
α
k
S
k
).
(14) Setz
k
+
1
=
z
k
+
α
k
S
k, k
←
k
+
1. Go to (2).Illustration of the optimal control and state achieved by solving Example1, taking u
^
(
x
,
y
)
=
exp
⁡
⁡
[
-
64
(
(
x
-
(
1
/
2
)
)
2
+
(
y
-
(
1
/
2
)
)
2
)
], N
=
3, K
=
5, and α
=
0.00005.
(a)
Desired values
(b)
Initial state
(c)
Approximated solution (state)
(d)
Control
(e)
Adjoint
(f)
Expected value
(g)
Standard dev.
(h)
Sparsity of global matrixn
z
=
1815980In Figure1(a), the desired function u
^
(
x
,
y
) is plotted for x
i
=
y
i
=
-
1
+
0.04
*
i,  i
=
0,1
,
2
,
…
,
50. Initial state is computed and plotted in Figure 1(b) by solving (37) where z
→
=
0. Approximated u
opt corresponding to the optimal control z
opt is calculated by Algorithm 1 and its graph is depicted in Figure 1(c). As it was expected, by using the proposed Algorithm 1, even by very far initial state, the approximated state solution reaches to the desired function. In Figure 1(d), the z
opt is depicted. Serious changes during the computation for initial z
→
=
0 happened reach to the optimal control function z
opt. The adjoint state is computed and plotted in Figure 1(e). As it is expected from (58), the adjoint state corresponding to the z
o
p
t must approximate the initial state (see Figures 1(b) and 1(c)). Figures 1(f) and 1(g) represent counter plots of the optimal state expectation and standard deviation. Finally, in Figure 1(h), the representation of coefficient matrix sparsity with the number of nonzero component in the computation is plotted in x
-
y plane.In Example1 we considered an exponential desired function, with v-sharp points. In the next example the smooth desired function and boundary conditions are considered. This is the case that Algorithm 1 is more efficient to use.Example 2.
Consider the optimal control of Example1, in which the right hand side function as well as boundary conditions is replaced by u
^
(
x
,
y
)
=
2
π
2
sin
⁡
(
π
x
)
sin
⁡
(
π
y
), N
=
3, and K
=
4, which emphasize that P
=
34. Taking α
=
0.00001, the illustration of computed optimal control and corresponding state in addition to the global matrix sparsity, standard deviation and expected values of the solution, adjoint values and initial state approximation of the desired function is shown in Figure 2.Illustration of optimal control and state achieved by solving Example2, taking u
^
(
x
,
y
)
=
2
π
2
sin
⁡
(
π
x
)
sin
⁡
(
π
y
), N
=
3, K
=
4, and α
=
0.00001.
(a)
Desired values
(b)
Mean values
(c)
Approximated solution (State)
(d)
Control
(e)
Adjoint
(f)
Expected value-SGS
(g)
Standard dev.-SGS
(h)
Sparsity of global matrixn
z
=
3046160


## 4. Conclusion

The solution of SPDE-constrained optimization problems is a recently challenging computational task. Here, we consider distributed control problems in which stochastic diffusion equation is the SPDE. Since the saddle point system extracted from using KKT optimality condition of the problem is a very large system and more expensive to solve, hence, we use the strategy of adjoint technique and preconditioned Newton’s conjugate gradient method, which iteratively solve the problem and has low computational cost. By separating the stochastic and deterministic parts of the SPDE using KLE and discretizing each part by WCE and Galerkin finite element method, respectively, we adjoint technique to compute the gradient and Hessian of the discretized optimization problem. By using preconditioned Newton’s conjugate gradient method the optimal control and state of the problem are calculated numerically. Two numerical examples are given to illustrate the correspondence between theoretical and numerical approaches.



---

*Source: 974305-2014-02-20.xml*