1. INTRODUCTION
Design optimization deals with exploring multiple instances of a design with the aim of identifying the best among them. In the presence of uncertainty, each of these instances entails the exploration of several nominally equivalent realizations intended to cover the range of anticipated scatter. In cases of interest to the present study, this uncertainty and the associated scatter are described as probabilistic measures on parameters or environmental conditions associated with each design instance. In the context of model-based optimization, each instance of the design involves the evaluation of a model, typically within a computational setting. Under conditions of uncertainty, the model is probabilistic, and each evaluation involves the exploration of the set of possible events, thus considerably increasing the number of function evaluations required to identify the optimal design. Clearly, the characterization of the probability measure on parameter space is a critical step in proceeding with any optimization under uncertainty (OUU) problem. That step itself can be cast as an iteration over a model, requiring the repetitive evaluation of an objective function or a likelihood function. This calibration or update step notwithstanding, the added computational burden is a key challenge when introducing uncertainty to a computational problem. For large-scale computational models, this challenge is on the critical path of credible model-supported decision making and design.
Typically, the objective function in an OUU problem is the expectation of a quantity of interest (QoI), Q, with respect to the probability measure of the random variables that describe the uncertainty in the problem. This expectation could, for instance, represent average performance, or the probability of some particular event, often associated with failure. Monte Carlo sampling techniques are standard tools for evaluating these expectations, and several optimization procedures have been integrated as sampling strategies resulting in algorithms such as stochastic gradient descent (Bottou, Reference Bottou2010), sample average approximation (Kleywegt, Reference Kleywegt, Shapiro and Homem-de Mello2002), Robbins–Monro (Robbins & Monro, Reference Robbins and Monro1951), the simultaneous perturbation stochastic approximation (Spall, 1952), and random search (Heyman & Sobel, Reference Heyman and Sobel2003). Although adapted to specific problems, and generally more efficient than standard sampling strategies, these algorithms can still be prohibitive for problems with large uncertainties or expensive function evaluations. In addition to uncertainty in the objective functions, uncertainty can also be present in the constraints. In such chance-constrained problems, the task of verifying the constraints incurs a significant computational burden, and tolerances on satisfying probabilistic constraints must be introduced. Robust optimization problems have been introduced to alleviate the computational burden by replacing the stochastic problem with a deterministic equivalent (Ben-Tal et al., Reference Ben-Tal, El Ghaoui and Nemirovski2009). This equivalence often pertains to low-order statistics and may fail to capture extreme behavior, relevant in many design situations.
Reservoir modeling and subsurface hydrology are two of the areas where design optimization problems are of particular interest and where the computational model involved is usually very expensive to solve as a result of the complex multiphase/multicomponent physics of flow and transport in porous media that give rise to coupled systems of equations and the fine scale discretization of the domain that is imposed on them. As such, the need for efficient methods in related design optimization problems is crucial, and an extended literature has been devoted to it. The objectives are typically to infer subsurface properties such as permeability and porosity or to optimize decision-making strategies such as oil production and soil remediation. The design problem that often arises in this case is that of determining the optimal locations to place the wells in an oil reservoir or the boreholes in an aquifer such that some oil production-related function or some information-based criterion is maximized (Tsilifis et al., Reference Tsilifis, Ghanem and Hajali2017). Variations on the oil well placement problem using optimization algorithms (Rian & Hage, Reference Rian and Hage1994; Beckner & Song, Reference Beckner and Song1995), including uncertainty (Güyagüler & Horne, Reference Güyagüler and Horne2001; Güyagüler, Reference Güyagüler2002; Yeten et al., Reference Yeten, Durlofsky and Aziz2002; Bangerth et al., Reference Bangerth, Klie, Wheeler, Stoffa and Sen2006), have explored well placement scenarios, details of the flow model, or specifics of the cost function.
One way to cope with the computational challenges described above is to replace the forward model with a cheap-to-evaluate surrogate such as polynomial chaos (PCE; Eldred, Reference Eldred2009; Keshavarzzadeh et al., Reference Keshavarzzadeh, Meidani and Tortorelli2016) or Gaussian process (Bilionis & Zabaras, Reference Bilionis and Zabaras2013). Once such a surrogate is available, the computational cost required to perform an optimization algorithm is reduced dramatically; however. there is no free lunch. Constructing the surrogate will often be an expensive procedure and can become even prohibitive for high-dimensional inputs. In this work, we develop an alternative framework based on polynomial chaos expansions that reduces further the cost of constructing the surrogate and therefore can be applicable in higher dimensional settings where standard nonintrusive methods would fail. The idea is based on constructing adapted polynomial chaos expansions that was first introduced in Tipireddy and Ghanem (Reference Tipireddy and Ghanem2014), which relies on rotating the Gaussian input of the chaos expansion and constructing a reduced expansion dependent on only a few Gaussian inputs, thus reducing the input dimension and requiring a signicantly smaller number of forward evaluations for the estimation of its coefficients. The novel contribution of the present work is in considering the optimization objective function as the QoI used in adaptation. In this manner, an adaptation is constructed for each iteration in the design space. This procedure permits the efficient evaluation of probabilistic operations (expectations or probabilities of failure) at each design point. We note that the proposed formalism is very different from sparse PCE procedures where a subset of the initial random variables is used in representing the solution. Substantial additional reduction is achieved in the present work by seeking direction in the parameter space along which the solution is concentrated. Typically, all of the original components are active along these directions and will thus contribute to the final solution.
The paper is structured as follows: Subsection 2.1 describes the design optimization objectives. In Subsection 2.2 we present the polynomial chaos properties and the basis adaptation concept along with the computational algorithm for its implementation. In Subsection 2.3 we describe how the polynomial chaos basis adaptation framework can be applied in a design optimization problem and a computational algorithm is also provided. In Section 3 we apply our methodology in a standard oil well placement problem that consists of a two-dimensional oil reservoir with random permeability where one production well and one injection well need to be placed such that threshold above which the cummulative production rate will be with a high confidence level, is maximized. Conclusions are presented in Section 4.
2. METHODOLOGY
2.1. Design optimization objectives
We present the general objective of the design optimization problem under uncertainty that we are investigating. We consider a scalar QoI
$Q({\bf p}\comma \,{\bf \theta} ) \in {\ssf R}$
that can be thought of as the outcome of a physical model or a computer code that depends on design parameters
${\bf p} \in {\ssf R}^n $
and some random parameters
${\bf \theta} \in {\ssf R}^m $
that are defined on some probability space
$\lpar {\Omega\comma \, {\ssf F}\comma \, P} \rpar $
and follow some probability distribution that is assumed to be known. Let
$F_Q (q) = P(Q({\bf p} \comma \,{\bf \theta} ) \le q)$
be the cummulative probability distribution of the quantity of interest evaluated at
${\bf p}$
. It is a common task to seek the estimation of the value q
α that provides a one-sided confidence at level α ∈ (0, 1) (typically small). Thus, denoting the cumulative distribution function with F
Q
and the quantile with q
α,
$$F_Q (q) = P(Q({\bf p}\comma \,{\bf \theta} ) \le q)\comma \,\quad q_{\rm \alpha} = F_Q^{ - 1} ({\rm \alpha} )\comma$$
and the objective is to find q
α such that
$Q({\bf p}\comma \,{\bf \theta} )$
is above q
α with probability 1 − α. Clearly, q
α depends on
${\bf p}$
. Motivated by this, we are interested in solving the optimization problem of the form
$${\bf p}^{\ast} = {\arg} \mathop {\max} \limits_{{\bf p} \in {\ssf X}} J({\bf p})\comma \,\quad J({\bf p}) = q_{\rm \alpha}\comma $$
where
${\ssf X} \subset {\ssf R} ^n $
is the domain of the design parameters (typically bounded). The dependence of q
α on
${\bf p}$
comes through the distribution F
Q
. For evaluation of the objective function subject to minimization, knowledge of this distribution is required. In all but the simplest cases, it will not be possible to compute F
Q
analytically and only an estimate of it will be feasible. This could typically be estimated based on empirical procedures, for instance, a histogram obtained from Monte Carlo (MC) sampling or similar sampling-based techniques such as kernel density estimation (Scott, Reference Scott2015). The main drawback in these approaches is that the number of MC estimates required for such estimations can be prohibitive in terms of the computational resources necessary to carry out such tasks. Among other characteristics that affect the required number of samples are the complexity of the computational model involved and the dimensionality and the distribution of
${\bf \theta} $
. These both affect the structure of the probability distribution that we seek to estimate. In the case where the computational model involved is replaced by a surrogate, the computational burden of MC sampling does not necessarily disappear as the cost of actually constructing the surrogate can be equally high, if not higher. Our ultimate goal in this work is to provide an alternative methodology that attempts to reduce the cost of estimating the distribution of F
Q
that is applicable in the case where the model uncertainties can be written as functions of Gaussian noise and the QoI admits a homogeneous chaos expansion.
2.2. Polynomial chaos basis adaptation
2.2.1 Homogeneous chaos
In the previous paragraph, we introduced uncertainty in our QoI through the random vector
${\bf \theta} $
whose distribution is assumed known. In order to develop our methodology, it is further necessary to assume that
${\bf \theta} $
can be written as a function of d uncorrelated Gaussian variables
${\bf \xi} \in {\ssf R} ^d $
and write
${\bf \theta} = {\bf g}({\bf \xi} )$
. Note that d can be considered to be infinite, but here we will only consider the finite case. Clearly Q will then be also written as a function of the same Gaussian variables as
$Q({\bf p}\comma \,{\bf g}({\bf \xi} ))$
. For the sake of simplicity we will denote
$Q({\bf p}\comma \,{\bf \xi} )$
, ignoring the dependence through
${\bf g}$
when there is no confusion. One necessary condition for the above assumption to hold is for
${\bf \theta} $
to have finite second moment and be measurable with respect to the σ-algebra
${\ssf F} ({\bf \xi} )$
, generated by
${\bf \xi} $
. Then, the Cameron–Martin theorem (Cameron & Martin, Reference Cameron and Martin1947) ensures that
${\bf \theta} $
admits a representation of the form
$${\bf \theta} = \sum\limits_{ \vert {\rm \alpha} \vert = 0}^p {\bf \theta} _{\rm \alpha} \Psi _{\rm \alpha} ({\bf \xi} )\comma \, \quad j = 1\comma \, \ldots\comma \, m\comma \,\quad {\bf \theta}\comma \, {\bf \theta} _{\rm \alpha} \in {\ssf R}^m\comma $$
that is, mean-square convergent as p → ∞, α = (α1, . . . , α
d
),
$\vert {\rm \alpha} \vert = \sum\nolimits_{i = 1}^d \!{\rm \alpha} _i $
, and
$\Psi _{\rm \alpha} ({\bf \xi} )$
are multidimensional Hermite polynomials of order α defined as
$$\Psi ({\bf \xi} ) = \prod \limits_{i = 1}^d {\rm \psi}_{\alpha _i} (\xi _i )\;.$$
Here,
${\rm \psi} _{{\rm \alpha} _i} ({\rm \xi} _i )$
are the normalized one-dimensional Hermite polynomials of order α
i
. Assuming that q has also finite second moment, it is clear that a similar representation can be obtained:
$$Q({\bf p}\comma {\;\bf \xi} ) = \sum\limits_{ \vert {\rm \alpha} \vert = 0}^p Q_{\rm \alpha} ({\bf p})\Psi _{\rm \alpha} ({\bf \xi} ).$$
The total number P of terms in the expansion is given by
$$P + 1 = \displaystyle{{(d + p)!} \over {d!p!}}.$$
That the functions {Ψα}α form an orthonormal basis on the space of square integrable random variables
$L^2 (\Omega\comma \, {\ssf F} ({\bf \xi} )\comma \,P)$
gives a convenient expession for the coefficients
$$Q_{\rm \alpha} ({\bf p}) = \langle Q({\bf p}\comma .)\comma \,\Psi _{\rm \alpha} \rangle \comma $$
which can be estimated using several techniques, the most popular of which are Galerkin projections (Ghanem & Spanos, Reference Ghanem and Spanos1991; Le Maître, Reference Le Maître, Reagan, Najm, Ghanem and Knio2002) and pseudo-spectral projections (Reagan et al., Reference Reagan, Najm, Ghanem and Knio2003) that are based on numerical integration using quadrature rules, known also as intrusive and nonintrusive approaches, respectively. In the application included in this work, we will focus on the second approach, which allows treating the computational model as a black box. Once the coefficients
$Q_{\rm \alpha} ({\bf p})$
are estimated, one can reconstruct the density function of Q through samples of
${\bf \xi} $
, as evaluation of the polynomial series is easy and fast and one can proceed in solving optimization problems of the form in Equation (2). However, the number of evaluations required for estimation of the coefficients grows exponentially fast as a function of the dimensionality of
${\bf \xi} $
and quickly becomes prohibitive or at least not cheaper, for the purpose of solving Equation (2), than direct sampling using the forward model. One way to address this issue, as we demonstrate in the next section, is to adapt the Gaussian basis of the expansion in a fashion that reduces the dimensionality and consequently requires fewer evaluations.
2.2.2. Basis adaptation
Let
${\bf A}:{\ssf R} ^d \to {\ssf R} ^d $
be an isometry on
${\ssf R} ^d $
, or equivalently, an orthogonal matrix (
${\bf AA}^T = {\bf I}_d $
) and define
${\bf \eta} = {\bf A}{\bf \xi} $
. It follows (Janson, Reference Janson1999) that the components of
${\bf \eta} $
are also uncorrelated standard Gaussian variables that span the same Gaussian Hilbert space with
${\bf \xi} $
and generate the same σ-algebra
${\ssf F} ({\bf \xi} )$
; therefore, a polynomial chaos expansion of any element
$Q \in L^2 (\Omega\comma \, {\ssf F} ({\bf \xi} )\comma \,P)$
can be taken with respect to the new germ
${\bf \eta} $
as
$$Q^{\bf A} ({\bf p} \comma\, {\bf \eta} ) = \sum\limits_{{\rm \beta\comma\beta} = 0}^p Q_{\rm \beta}^{\bf A} ({\bf p}){\rm \psi}_{\rm \beta} ^{\bf A} ({\bf \xi} )\comma $$
where we use the superscript on
$Q^{\bf A} $
and
$Q_{\rm \beta}^{\bf A} $
to indicate that the new basis has been obtained after applying rotation
${\bf A}$
on
${\bf \xi} $
and
${\rm \psi} _{\rm \beta} ^{\bf A} ({\bf \xi} ) = {\rm \psi} _{\rm \beta} ({\bf A\xi} )$
. Then, we have that
$Q^{\bf A} ({\bf p}\comma \,{\bf \eta} ) = Q({\bf p}\comma \,{\bf \xi} )$
almost surely and the coefficients of the new expansion can be written in terms of
$\{{ Q_{\rm\alpha} ({\bf p})} _{\rm \alpha = 0}^{\vert p \vert}\} $
as
$$Q_{\rm \beta} ^{\bf A} ({\bf p}) = \sum\limits_{ \vert {\rm \alpha} \vert = 0}^p Q_{\rm \alpha} ({\bf p})\langle {\rm \psi} _{\rm \alpha} \comma\ {\rm \psi} _{\rm \beta} ^{\bf A} \rangle. $$
We have
$\langle {\rm \psi} _{\rm \alpha} \comma \, {\rm \psi} _{\rm \beta} ^{\bf A} \rangle = 0$
for |α| ≠ |β|, while for |α| = |β| the inner products can be computed explicitly (Tsilifis & Ghansem, Reference Tsilifis and Ghanem2016). The motivation behind applying a change of basis as described above is that, by choosing an appropriate matrix
${\bf A}$
, the new expansion can concentrate the dominant effects of the original expansion to fewer coefficients of
$Q^{\bf A} $
and thus provide sparse representations depending on only a few components of
${\bf \eta} = ({\rm \eta} _1\comma \, \ldots \comma \, {\rm \eta} _d )$
, say,
${\bf \eta} _{d_0} = ({\rm \eta} _1\comma \, \ldots \comma \, {\rm \eta} _{d_0} )$
with d
0 < d, and therefore, we can discard terms of the adapted expansion via projection (see Tipieddy & Ghanem, Reference Tipireddy and Ghanem2014; Tsifilis & Ghanem, Reference Tsilifis and Ghanem2016, for details as well as several popular choices of
${\bf A}$
). If we let
${\ssf J} _{\,p \comma d} = {\vert{\rm \alpha}\vert} = ({\rm \alpha} _1\comma \, \ldots\comma \, {\rm \alpha} _d ): {\vert{\rm \alpha}\vert} \le {\vert p\vert} $
be the set of d-dimensional multi-indices of order up to P and take any
${\ssf J} _{\,p_0\comma d_0} $
for some p
0 < p, d
0 < d, which implies that
${\ssf J} _{\,p_0\comma d_0} \subset {\ssf J} _{\,p\comma d} $
, then for given
${\bf A}$
, the reduced decomposition
$$Q^{\bf A} ({\bf p} \comma\, {\bf \eta} _{d_0} ) = \sum\limits_{{\rm \beta} \in {\ssf J} _{\,p_0 \comma d_0}} Q_{\rm \beta} ^{\bf A} ({\bf p}){\rm \psi} _{\rm \beta} ({\bf \eta} _{d_0} )$$
approaches
$$Q^{\bf A} ({\bf p}\comma \,{\bf \eta} )\mathop = \limits^{[0pt] \;\;almost\ surely} Q({\bf p}\comma \,{\bf \xi} )$$
as p
0 → p and d
0 → d. One then can carefully choose p
0, d
0 such the reduced decomposition in Eq. (10) can be sufficiently close to
$Q({\bf p}\comma \,{\bf \xi} )$
.
Obtaining this expression requires the estimation of
$Q_{\rm \beta}^{\bf A} ({\bf p})$
,
${\rm \beta} \in {\ssf J} _{\,p_0\comma d_0} $
. In practice, using Eq. (9) does not deliver us from the computational cost of obtaining a full polynomial chaos expansion because the coefficients
$Q_{\rm \alpha} ({\bf p})$
are required. Instead, once an isometry
${\bf A}$
is chosen, one can compute the coefficients
$Q_{\rm \beta} ^{\bf A} ({\bf p})$
of the adapted lower dimensional expansion by following a nonintrusive implementation as described in Tipireddy and Ghanem (Reference Tipireddy and Ghanem2014). The computational benefits of doing so is that for
$d_0 = d$
, the number of quadrature points required for a nonintrusive estimation of the coefficients
$Q_{\rm \beta} ^{\bf A} ({\bf p})$
is much smaller than that of the original expansion while the adapted expansion serves as an accurate estimate of our QoI. The details of the procedure are summarized in Algorithm 1.
Nonintrusive implementation
Throughout this paper we restrict ourselves with obtaining
${\bf A}$
by using the Gaussian adaptation. This consists of first estimating the coefficients
$\{{ Q_{{\rm \varepsilon} _i} }\} _{i = 1}^d $
of the first-order polynomials
$ \{{ {\rm \psi} _{{\rm \varepsilon} _i} ({\rm \xi} _i )}\} _{i = 1}^d $
, where ε
i
= (0, … , 1, … , 0), that is, the multi-index with 1 at location i and 0 everywhere else. This is achieved using a low-level quadrature rule with very few points. Next, we define the vector
${\bf Q}_1 = [Q_{{\rm\varepsilon}_1}\comma \, \ldots\comma \, Q_{{\rm\varepsilon}_d} ]$
, and we set the first row of
${\bf A}$
, denoted with
${\bf A}_{1 \cdot} $
, to be
$${\bf A}_{1 \cdot} = \displaystyle{{{\bf Q}_1} \over { \vert \vert {\bf Q}_1 \vert \vert }}. $$
The matrix
${\bf A}$
can then be completed by performing a Gram–Schmidt procedure.
2.3. Design optimization using adapted polynomial chaos
The adaptation framework described in the previous paragraph can be implemented in order to obtain cheap surrogates of the QoI at different locations
${\bf p}$
, which can then be used for estimation of F
Q
and its inverse. The optimization problem described in Section 2.1 can then be approached as follows: first, one needs to employ an optimization algorithm to be used for maximization of
$J({\bf p})$
. At each iteration k, the adapted chaos expansion of the QoI needs to be computed at design parameters
${\bf p}_k $
, for a dimension value d
0 that achieves a satisfactory level of accuracy. Then F
Q
can be estimated with MC sampling, and
$F_Q^{ - 1} ( \cdot )$
can be evaluated at α, thus giving the values of the objective function
$J({\bf p})$
. The choice of d
0 is not known a priori, and it needs to be identified. On way to do this is to simply start by constructing a one-dimensional adapted expansion and then increase the dimension until some relative error of the consecutive adaptations falls below some threshold. Because the number of evaluations required increases exponentially as a function of the dimension, the evaluations required to compute all adaptations up to dimension d
0 will be still much smaller than those required for the d-dimensional one as long as d
0 remains significanly smaller than d. Algorithm 2 displays the pseudocode with all the steps required for one evaluation of
$J({\bf p})$
.
Evaluation of J(p).
The outside loop, where d 0 is being incremented, may seem to present an additional workload beyond standard OUU procedures. It is noted, however, that for values of d 0 well below the value of d, the computational burden associated with evaluations for dimensions up to and including d 0 is still a fraction of the cost for an evaluation at dimension d. While the numerical value of the tolerance, ε, can be selected to become smaller as d 0 is approached, in the present paper we fix ε to a predetermined value of 0.05.
3. APPLICATION: OIL WELL PLACEMENT PROBLEM
3.1. The oil reservoir model
To show the computational efficiency of the proposed method, we consider a multiphase flow in a petroleum reservoir. The governing equations of the corresponding physics describing black oil flow through a heterogeneous porous medium are the following (Chen et al., 2008):
$$\eqalign{&{\displaystyle{{\partial} \over {{\partial} t}}\left[ {\displaystyle{{{\rm \phi} {\rm \rho}_{W_s}} \over {B_w}} S_w} \right] + \nabla \cdot \left( {\displaystyle{{{\rm \rho}_{W_s}} \over {B_w}} \vec u_w} \right) - q_W = 0}\comma \cr & {\displaystyle{{\partial} \over {{\partial} t}}\left[ {\displaystyle{{{\rm \phi} {\rm \rho}_{O_s}} \over {B_o}} S_o} \right] + \nabla \cdot \left( {\displaystyle{{{\rm \rho}_{O_s}} \over {B_o}} \vec u_o} \right) - q_O = 0 \comma } \cr & {\displaystyle{{\partial} \over {{\partial} t}}\left[ {{\rm \phi} \left( {\displaystyle{{R_{{\rm so}} {\rm \rho}_{G_s}} \over {B_o}} S_o + \displaystyle{{{\rm \rho}_{G_s}} \over {B_g}} S_g} \right)} \right] + \nabla \times \left( {\displaystyle{{R_{\rm so} {\rm \rho}_{G_s}} \over {B_o}} u_o + \displaystyle{{{\rm \rho}_{G_s}} \over {B_g}} u_g} \right) } \cr & \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad - q_G = 0.}$$
where ϕ denotes the porosity;
$S_{\rm \alpha} \comma \,\mu _{\rm \alpha} \comma \,p_{\rm \alpha} \comma \,B_{\rm \alpha} \comma \, {\rm \rho} _{{\rm \alpha} _s}\comma \hbox{ and } q_{\rm \alpha} $
are the saturation, viscosity, pressure, formation volume factor, density, and volumetric rate at each phase α (water, gas, and oil); R
so is the gas solubility; u
w
, u
g
, and u
o
are the respective Darcy velocities obtained by solving Darcy's law for each phase
$$\eqalign{u_{\rm \alpha} &= - \displaystyle{{k_{r {\rm \alpha}} k({\bf x} \comma\, {\rm \omega})} \over {{\rm \mu} _{\rm \alpha}}} \lpar {\nabla p_{\rm \alpha} - {\rm \rho}_{\rm \alpha} {\rm \gamma} \nabla z} \rpar\cr {\rm \alpha} &= o \comma\, w \comma\, g \comma \quad {\bf x} \in \Gamma \comma \quad {\rm \omega} \in \Omega \comma}$$
where ρα and k
rα are the respective mass density and relative permeabilities of each α phase, γ is the magnitude of the gravitational acceleration, z is the depth, and
$k({\bf x})$
is the intrinsic permeability. The relative permeabilities k
rα are functions of saturation, and their values generally lie between 0 and 1. They are obtained either experimentally or by using analytical expressions such as the Corey (Reference Corey1994) and Stone (Reference Stone1970, Reference Stone1973) models. In the above equations, fluid injection wells and production wells are represented by point sources and sinks whenever they are present (Peaceman, Reference Peaceman2000). Uncertainty is introduced in our model by considering the absolute permeabilities
$k({\bf x})$
to be log-normally distributed, that is,
$\log k({\bf x}): {\ssf N} ({\bf m}({\bf x})\comma \,{\bf C}({\bf x}\comma \,{\bf x}))$
for all
${\bf x}$
. Here we use the SPE-10 data set (Christie & Blunt, Reference Christie and Blunt2001), which is discretized in a 60 × 220 × 85 regular Cartesian grid with real range being 600 × 2200 × 170 (ft.)3. We treat the 85 layers as different realizations of our two-dimensional domain to construct our sample mean and sample covariance matrix to be used as our
${\bf m}({\bf x})$
and
${\bf C}({\bf x}\comma \,{\bf x'})$
, respectively. Next we parameterize our random field by considering a Karhunen–Loeve expansion of
$k({\bf x})$
as
$$k({\bf x}) = {\bf m}({\bf x}) + \sum\limits_{i = 1}^{20} \sqrt {{\rm \lambda} _i} {\rm \phi} _i ({\bf x}){\rm \xi} _i\comma $$
truncated up to 20 terms. The mean and variance of the generated log-permeability field are shown in Figure 1, and two sample realizations drawn from the distribution are shown in Figure 2. The permeability values used in our simulations are measured in milidarcy. We should note that the specific nature of the distribution is not relevant for the applicability of the proposed method. Under general conditions on their moments, most data-driven distrbutions can be represented as polynomials of Gaussian random variables.
Mean (left) and variance (right) of the log-pemeability field.
Sample realizations of the log-permeability field.
3.2. Oil production as a QoI
The quantity of interest
$Q({\bf p}\comma \,{\bf \theta} )$
considered in this example is taken to be the cumulative oil production over T = 2000 days, where the design parameters
${\bf p} = (x_1^i\comma \, x_2^i\comma \, x_1^p\comma \, x_2^p )$
are the locations
$(x_1^i\comma \, x_2^i )$
and
$(x_1^p\comma \, x_2^p )$
of one injection and one production well, respectively, and
${\bf \theta} $
is taken here to be the random permeability field. The cumulative production is defined as
$$Q({\bf p} \comma {\bf \theta} ) = \int_E \int_T u_o ({\bf p} \comma\, {\bf \theta} ) \times d{\bf E}\,dt \comma $$
where
$u_o ({\bf p}\comma \,{\bf \theta} )$
is the Darcy velocity of the oil phase, E is the cross-sectional area of the wellbore, and
${\bf E} = E{ \hat n}$
is the vector area, where
${ \hat n}$
is the unit vector normal to the area. The objective function that we are interested in maximizing is
$$J({\bf p}) = q_{\rm \alpha} = F_Q^{ - 1} (\rm \alpha ).$$
3.3. Optimization results
Because the objective function considered here is complex in nature and its theoretical properties are not known, it can exhibit nonconvex behavior. In order to avoid exploring only an area around the local minimum, we start our convex optimization algorithm with an initial guess obtained from the brute force grid parameter study. For this study, the parameter space
${\bf p}$
is divided into four evenly spaced zones
$(x_{1g}^i = \{0 \comma \,200 \comma \,400 \comma\,600\}, x_{2g}^i = \{ 0 \comma \,733.3 \comma \,1466.7 \comma \,2200\}, x_{1g}^p = \{ 0 \comma \,200 \comma \,400 \comma \,600\}, x_{2g}^p = \{ 0 \comma \,733.3 \comma \,1466.7 \comma \,2200\})$
, and the objective function is evaluated at the corresponding grid design points
${\bf p}_g = (x_{1g}^i \times x_{2g}^i \times x_{1g}^p \times x_{2g}^p )$
. The values of the objective function evaluated at the 256 grid design points
${\bf p}_g $
are shown in Figure 3, highlighting with red, the ones that give oil production above 100 with probability 0.93. The possible locations of the wells where
$J({\bf x})$
are evaluated and the location of the well combinations resulting in maximum q
α are shown in Figure 4 (left and middle). The largest value of the objective function corresponds to the design point (0, 146.67, 60, 146.67), which we use as the initial guess to run the limited memory Broyden–Fletcher–Goldfarb–Shanno algorithm (Byrd et al., Reference Byrd, Lu, Nocedal and Zhu1995) that can perform bound constrained optimization, in order to find the optimal design point. The well locations generated from iterations of the algorithm and the resulting optimal design location of the injection and production wells (I and P labels, respectively) are shown in the Figure 4 (right). Figure 5 shows the water and oil saturation on the 2000th day at the optimal design point corresponding to the mean permeability field, and Figure 6 shows a sample realization of the water and oil saturation on the same day and at the same design point.
Values of q
α at the 256 design points. The values with red “
${\bf W} $
” are above 100.
Left: Injection and production well locations generating 256 possible combinations. Middle: 8 pairs of well locations resulting in q α > 100. Right: Pairs of well locations produced by running limited memory Broyden–Fletcher–Goldfarb–Shanno optimization.
Realization of the oil saturation and water saturation on the 2000th day at the optimal design point using mean permeability.
Sample realization of the oil and water saturation on the 2000th day at the optimal design point.
Finally, in order to perform a comparison of the computational needs of our approach with possible alternatives, it is necessary to look at the dimensionality of the adapted expansions used in the evaluations of the objective function. Figure 7 shows the histogram with the dimension of the adapted expansions used in the estimation of F
Q
(·) at all design points. We see that the dimensionality of the
${\bf \eta} $
components necessary for accurate estimation of F
Q
(·) does not exceed 7, and no more than 30 design points require dimension larger than 3. Furthermore, Figure 8 depicts the contribution of each ξ
i
to η1, … η7 at one of the design configuration where the dimension of the adapted dimension is 7. These contributions are obtained as the respective rows of the adaptation matrix
${\bf A}$
. This figure reveals that all 20 ξs are contributing to the adapted dimensions, showing that the quantity of interest depends on all 20 dimensions in the ξ space. A polynomial order of 3 in the reduced dimension d
0 is pursued in all subsequent calculations. The number of quadrature points, and consequently the number of forward evaluations required for the estimation of the coefficients of the adapted expansions of order 3 for different dimensionality of
${\bf \eta} _{d_0} $
, is shown in Table 1. Smolyak sparse quadrature rule at Level 2 is used for estimation of multidimensional integrals.
Dimension of the adapted chaos expansions used in all 256 design points.
Contribution of each ξ1 to η1, … , η.
Number of forward evaluations required for estimation of d-dimensional adapted chaos expansions for d = 1, … , 7, with polynomial order = 3 and level 2 Smolyak quadrature
3.4. Variability of the isometry
${\bf A}$
In this subsection, we investigate the dependence of the isometry
${\bf A}$
on the design
${\bf p}$
, and write
${\bf A}({\bf p})$
. This dependence is clearly a result of the dependence of the coefficients
$q_{\rm \alpha} ({\bf p})$
of
$q({\bf p}\comma \,{\bf \xi} )$
on
${\bf p}$
and particularly the ones corresponding to the first-order Hermite polynomials that are used in the first row of
${\bf A}$
. We perform a simple test that aims to provide some further insight on how
${\bf A}$
varies as a function of
${\bf p}$
. The motivation for doing so is that if this relation is known a priori to the experimenter, then the computation cost involved in the algorithms described above, which is related to the evaluation of
${\bf A}$
, could be omitted, thus resulting in even more efficient implentation. The idea is simply that
${\bf A}({\bf p})$
, or at least its first row
${\bf A}$
1, which is defined in Eq. (11), is often a smooth function of
${\bf p}$
, and therefore, we should be able to estimate it; for instance, we could evaluate
${\bf A}$
1 at several points and then apply some interpolation scheme. A brief theoretical discussion on the smoothness of the entries of
${\bf A}$
can also be found in (Tsifiis & Ghanem, Reference Tsilifis and Ghanem2016), where it is shown that if
$Q_\alpha ({\bf p})$
are square integrable functions of
${\bf p}$
, then the covariance kernels of η
i
s define Hilbert–Schmidt operators. Here, we only demonstrate that in our example this claim is true and the entries of
${\bf A}$
1 or some properties of them could be estimated. Further development of interpolation or estimating procedures is beyond the scope of the present paper.
To simplify the dependence on
${\bf p}$
, we fix the location of the injection well in our example to be at the point (1, 2200) and thus
${\bf p} = (x_1^p\comma \, x_2^p )$
becomes two dimensional. We compute the Gaussian part
$q_{\rm \alpha} ({\bf p})$
, |α| = 1 at 256 locations as shown in Figure 8 (right), and from that we find the angles
$${\rm \phi} = {\rm arccos}({\bf A}_{1 \cdot} ({\bf p}_o ) \times {\bf A}_{1 \cdot} ({\bf p})) \comma $$
that is, the angles between the row
${\bf A}$
1·(
${\bf p}$
) and
${\bf A}_{1 \cdot}({\bf p}_o )$
for
${\bf p}_o = (600\comma \,1)$
and all other
${\bf p}$
. There are shown in Figure 9 (left). We observe that the angle between the first rows of the rotations
${\bf A}({\bf p})$
and
${\bf A}({\bf p}_o )$
is increasing as the production well moves to the upper area of the domain, that is, it approaches the injection well.
Spatial distribution of rotation angle for a fixed injection well, as a function of the production well.
3.5. Effect of the confidence level on the optimal design
In this section we are interested in the effect of the confidence level (1−α) on the design. For this study, similarly to Section 3.3 the parameter space
${\bf p}$
is divided into four evenly spaced zones
$(x_{1g}^i = \{ 0 \comma \, 200 \comma \, 400 \comma \,600\} \comma \,x_{2g}^i = \{ 0 \comma \, 733.3 \comma 1466.7 \comma \, 2200\}\comma \,x_{1g}^p = \{ 0 \comma \, 200 \comma \, 400 \comma \,600\}\comma \, x_{2g}^p = \{ 0\comma \,733.3 \comma 1466.7 \comma \, 2200\} )$
, and the objective function is evaluated at the corresponding grid design points
${\bf p}_g = (x_{1g}^i \times x_{2g}^i \times x_{1g}^p \times x_{2g}^p )$
. Each of the grid points designated as design number 1, 2, … , 256. The values of the objective function evaluated at the 256 grid points for different confidence levels 30, 35, 40, … , 95. For each confidence level, the optimal design is picked where the objective function is maximum. Figure 10 (left) shows how the optimal design (here its index number) changes for different value of the confidence level. We observe that there are essentially only three possible optimal designs along the whole confidence level spectrum. The production and injection wells ({P
i
, I
i
}
i = 1,2,3) for these three optimal designs are shown in Figure 10 (right). Here, the optimal design one, two, and three corresponds to confidence levels 0.2 to 0.4, 0.45 to 0.65 and 0.70 to 0.95, respectively. Figure 11 (left) shows the quantile values q
α obtained at the optimal design point as a function of the confidence level, and Figure 11 (right) shows the probability density functions of q
α for the three optimal designs. It is worth mentioning here that the support and the shape of the probability density functions indicates that although a low confidence level will result most likely in a low minimum production q
α, there is still a positive probability that a large value (up to around 4000) of oil production might be achieved. As the confidence level increases, the support of the probability density function shrinks, indicating a lower upper bound for q
α that, however, can be achieved with higher probability, and thus the minimum oil production is steadily increasing while the associated risk is decreasing.
Confidence level versus the optimal design number (left) and the locations of the optimal designs on the domain (right).
Dependence of the optimal production level on confidence level (1–α) (left) and probability density function Q α of the three optimal designs (right).
4. CONCLUSIONS
We presented a polynomial chaos based approach to design optimization that makes use of the basis adaptation capabilities of the PCE in order to obtain highly efficient reduced expansions that offer an alternative to design optimization problems for a special case of objective functions that consist of failure probabilities and thus require full knowledge of the density function of the QoI. We applied our methodology to the oil well placement problem with an expensive-to-solve reservoir model, and our numerical results indicate that our approach can provide solutions that dramatically decrease the computational cost that would be required by computing a polynomial chaos expansion on the initial basis. Much is left to be explored and understood in this novel direction in design optimization under uncertainty based on stochastic expansions, including the behavior of the rotation matrix as a function of the design parameters as well possible constructions of the rotation matrix other than the Gaussian adaptation that is used here.
ACKNOWLEDGMENTS
This research was supported by the ScramJet-UQ project funded under DARPA EQUIPS Program.
Charanraj Thimmisetty is a Postdoctoral Associate at the Center for Applied Scientific Computing at Lawrence Livermore National Laboratory. Dr. Thimmisetty has a PhD in civil and environmental engineering from the University of Southern California (USC).
Panagiotis Tsilifis is a Postdoctoral Associate in the Department of Civil Engineering at USC, where he attained his PhD in mathematics.
Roger Ghanem is the Gordon S. Marshall Professor of Engineering Technology at USC. He obtained his PhD from Rice University. Dr. Ghanem's current research deals with probabilistic uncertainty quantification and data-based and model-based predictability.
