Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Parameter estimation to infer injector-producer relationships in oil fields: from hybrid constrained nonlinear optimization to inference in probabilistic graphical model
(USC Thesis Other)
Parameter estimation to infer injector-producer relationships in oil fields: from hybrid constrained nonlinear optimization to inference in probabilistic graphical model
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
PARAMETER ESTIMATION TO INFER INJECTOR-PRODUCER
RELATIONSHIPS IN OIL FIELDS: FROM HYBRID CONSTRAINED
NONLINEAR OPTIMIZATION TO INFERENCE IN PROBABILISTIC
GRAPHICAL MODEL
by
Hyokyeong Lee
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(COMPUTER SCIENCE)
August 2010
Copyright 2010 Hyokyeong Lee
Dedication
Dedicated to God
ii
Acknowledgements
I would like to express my deepest gratitude to my advisor, Dr. Ke-Thia Yao and Professor
Aiichiro Nakano for their invaluable guidance, suggestions, patience, and encouragement
throughout this research. Their never-ending support enabled me to keep standing up
through the process.
I would like to thank Professor Robert Neches, for his continuous support and encour-
agement. He was such a person who showed me what a truly great mentor is.
I would like to thank Professor Iraj Ershaghi and Olu Okpani for the great technical help
as domain experts in petroleum engineering.
I would like to thank Professor Stephan Haas for serving on my qualifying exam and PhD
defense.
I would like to thank Professor Ulrich Neumann for serving on my qualifying exam.
I would like to thank all the students from countries who showed me that each of us
should be independent for our lives.
I would like to thank my family for the quiet support.
I would like to thank God for giving me the chance to break through, sending me all the
above helpers, and teaching me that any difficulty can be overcome.
iii
Table of Contents
Dedication ii
Acknowledgements iii
List of Tables vi
List of Figures viii
Abstract xi
Chapter 1: Introduction 1
1.1 Motivation ... ... .... ... .... ... ... .... ... .... . 1
1.2 Contributions .. ... .... ... .... ... ... .... ... .... . 5
Chapter 2: Capacitance-Resistive Model and Problem Formulation 8
Chapter 3: Hybrid Constrained Nonlinear Optimization 11
3.1 RelatedWork .. ... .... ... .... ... ... .... ... .... . 11
3.2 BaseAlgorithm . ... .... ... .... ... ... .... ... .... . 14
3.3 HybridConstrainedNonlinearOptimization . . ... .... ... .... . 19
3.4 PerformanceEvaluation ... ... .... ... ... .... ... .... . 22
3.4.1 ExperimentalSetting. ... .... ... ... .... ... .... . 23
3.4.2 Results . ... .... ... .... ... ... .... ... .... . 29
Chapter 4: Inference in Probabilistic Graphical Model 53
4.1 RelatedWork .. ... .... ... .... ... ... .... ... .... . 53
4.1.1 Markov and Bayesian Networks Structure Learning and Parameter
Estimation .. .... ... .... ... ... .... ... .... . 53
4.1.2 Factor Graph Structure Learning and Parameter Estimation . . . 56
4.2 FactorGraphandtheSum-ProductAlgorithm ... .... ... .... . 58
4.3 GraphicalRepresentationoftheCRM . . ... ... .... ... .... . 61
4.4 StructureLearningofFactorGraphsforParameterEstimation . .... . 63
4.5 PerformanceEvaluation ... ... .... ... ... .... ... .... . 68
4.5.1 ExperimentalSetting. ... .... ... ... .... ... .... . 68
4.5.2 Results . ... .... ... .... ... ... .... ... .... . 71
iv
Chapter 5: Conclusions 93
Bibliography 96
v
List of Tables
3.1 Simulated synthetic oil fields: Field name, field type, field size, number
of time steps, number of parameters, number of bounds and number of
equalityconstraints.. .... ... .... ... ... .... ... .... . 23
3.2 MaximumnumberofiterationsforBasevsHCNO . .... ... .... . 25
3.3 PopulationsizeofcGA.... ... .... ... ... .... ... .... . 26
3.4 Number of convergences, means of number of major iterations and search
time ... .... ... .... ... .... ... ... .... ... .... . 30
3.5 Means of root mean squared error (RMSE) and root relative squared error
(RRSE) . .... ... .... ... .... ... ... .... ... .... . 31
3.6 Number of negative ranks (HCNO - Base) and significance by Wilcox
signed-rank test for the number of major iterations, search time, RMSE
andRRSE.... ... .... ... .... ... ... .... ... .... . 35
3.7 Means of search time, root mean squared error (RMSE), and root relative
squarederror(RRSE)ofcGA ... .... ... ... .... ... .... . 36
3.8 Number of negative ranks (HCNO - cGA) and significance by Wilcox
signed-ranktestforsearchtime,RMSEandRRSE . .... ... .... . 37
4.1 Field name, field type, size, number of time steps, and number of parameters 68
4.2 Running times (Sec) of HCNO and PGL . . . . . . . .... ... .... . 71
4.3 Testerror(RMSE) . . .... ... .... ... ... .... ... .... . 72
4.4 True and estimated π
ij
forField4 . .... ... ... .... ... .... . 75
vi
4.5 True and estimated τ
ij
forField4 . .... ... ... .... ... .... . 76
vii
List of Figures
3.1 Well locations and true higher permeability . . . . . .... ... .... . 27
3.2 ExampleoftheinjectionandproductionsignalsofField1 . ... .... . 28
3.3 Example of the injection and production signals of Field2, 3, and 4 . . . . 28
3.4 Observed and predicted production curves of producer1 and 2 with the
leasttesterrorsoftestcaseno. 1ofBasealgorithmforField1. . .... . 33
3.5 Observed and predicted production curves of producer3 and 4 with the
leasttesterrorsoftestcaseno. 1ofHCNOforField1 ... ... .... . 34
3.6 Observed and predicted production curves of producer1 and 2 with the
leasttesterrorsoftestcaseno. 1ofBasealgorithmforField2. . .... . 38
3.7 Observed and predicted production curves of producer3 and 4 with the
leasttesterrorsoftestcaseno. 1ofHCNOforField2 ... ... .... . 39
3.8 Observed and predicted production curves of producer1 and 2 with the
leasttesterrorsoftestcaseno. 1ofBasealgorithmforField3. . .... . 40
3.9 Observed and predicted production curves of producer3 and 4 with the
leasttesterrorsoftestcaseno. 1ofBasealgorithmforField3. . .... . 41
3.10 Observed and predicted production curves of producer5 and 6 with the
leasttesterrorsoftestcaseno. 1ofBasealgorithmforField3. . .... . 42
3.11 Observed and predicted production curves of producer7 and 8 with the
leasttesterrorsoftestcaseno. 1ofBasealgorithmforField3. . .... . 43
3.12 Observed and predicted production curves of producer1 and 2 with the
leasttesterrorsoftestcaseno. 1ofBasealgorithmforField4. . .... . 44
viii
3.13 Observed and predicted production curves of producer3 and 4 with the
leasttesterrorsoftestcaseno. 1ofBasealgorithmforField4. . .... . 45
3.14 Observed and predicted production curves of producer5 and 6 with the
leasttesterrorsoftestcaseno. 1ofBasealgorithmforField4. . .... . 46
3.15 Observed and predicted production curves of producer7 and 8 with the
leasttesterrorsoftestcaseno. 1ofBasealgorithmforField4. . .... . 47
3.16 True (grey) and estimated (black) connectivity map for Field1: Base vs
HCNO .. .... ... .... ... .... ... ... .... ... .... . 49
3.17 True and estimated time constant map for Field1: Base vs HCNO . . . . 49
3.18 TrueandestimatedconnectivitymapforField2: BasevsHCNO .... . 50
3.19 True and estimated time constant map for Field2: Base vs HCNO . . . . 50
3.20 TrueandestimatedconnectivitymapforField3: BasevsHCNO .... . 51
3.21 True and estimated time constant map for Field3: Base vs HCNO . . . . 51
3.22 TrueandestimatedconnectivitymapforField4: BasevsHCNO .... . 52
3.23 True and estimated time constant map for Field4: Base vs HCNO . . . . 52
4.1 FactorGraphforEq.4.6... ... .... ... ... .... ... .... . 60
4.2 2×2fieldanditsfactorgraph .. .... ... ... .... ... .... . 61
4.3 Well locations and true higher permeability (grey lines) of four fields . . . 70
4.4 Observed and estimated production curves of Producer1 and 2 in Field1 . 74
4.5 Observed and estimated production curves of Producer3 and 4 in Field1 . 77
4.6 ObservedandestimatedproductioncurvesofField2 .... ... .... . 78
4.7 ObservedandestimatedproductioncurvesofField2 .... ... .... . 79
4.8 Observed and estimated production curves of proudcer1 and 2 in Field3 . 80
4.9 Observed and estimated production curves of proudcer3 and 4 in Field3 . 81
4.10 Observed and estimated production curves of proudcer5 and 6 in Field3 . 82
ix
4.11 Observed and estimated production curves of proudcer7 and 8 in Field3 . 83
4.12 Observed and estimated production curves of proudcer1 and 2 in Field4 . 84
4.13 Observed and estimated production curves of proudcer3 and 4 in Field4 . 85
4.14 Observed and estimated production curves of proudcer5 and 6 in Field4 . 86
4.15 Observed and estimated production curves of proudcer7 and 8 in Field4 . 87
4.16 True higher permeability and estimated connectivity in Field1 . . .... . 88
4.17 True higher permeability and estimated connectivity in Field2 . . .... . 89
4.18 True higher permeability and estimated connectivity in Field3 . . .... . 90
4.19 True higher permeability and estimated connectivity in Field4 . . .... . 90
4.20 True higher permeability and estimated time constant in Field1 . .... . 91
4.21 True higher permeability and estimated time constant in Field2 . .... . 91
4.22 True higher permeability and estimated time constant in Field3 . .... . 92
4.23 True higher permeability and estimated time constant in Field4 . .... . 92
x
Abstract
In petroleum community, the oil field optimization, i.e., minimizing the operational cost
and maximizing the oil recovery, is a challenging problem. One of the popular oil re-
covery techniques is waterflooding, which injects water into the oil reservoir to extract
oil. Thus, the knowledge about injector-producer relationships (IPRs), i.e., which injec-
tors contribute to the production of which producers, is a key for field optimization. The
difficulty associated with field optimization is that the underlying structure of oil fields is
unknown and it continuously changes over time. Recently, a capacitance-resistive model
(CRM) has been proposed to investigate the IPRs. The CRM is a predictive model that
predicts production rates given water injection rates. It consists of two sets of parameters,
connectivity and time constant. The connectivity parameter quantifies the contribution
of an injector to the production of a producer and the time constant parameter does
the degree of fluid storage between the two wells. Three constraints are posed on the
two sets of parameters: (1) the sum of the connectivity parameters of an injector to all
other producers should be 1; (2) individual connectivity parameters should fall between 0
and 1 (this implies a nonnegativity constraint as well); and (3) nonnegativity constraint
that the time constant parameters should be positive real numbers. Estimating parame-
ters is mapped to a parameter estimation problem for continuous nonlinear systems with
xi
constraints. The challenge is that an oil reservoir consists of hundreds of injection and
production wells, with which the number of CRM parameters increasing quadratically.
We propose two different approaches for the parameter estimation problem. First, we
propose a new constrained nonlinear optimization method. The parameter estimation
problem for the CRM is mapped to a large-scale constrained nonlinear optimization
problem. Nonlinear optimization is challenging due to the nonconvex property of the
objective function and/or constraints. Identifying a global optimum is extremely hard
even without constraints. The curse of dimensionality makes the problem more difficult
to solve: as the number of parameters increases, the number of local optima is likely to
increase. Moreover, there is no analytical solution for constrained nonlinear parameter
estimation problems due to the constraints. For the large-scale optimization problem,
we developed a hybrid constrained nonlinear optimization (HCNO) method, which is
based on sequential quadratic programming (SQP) with a quasi-Newton method in line
search framework. The constrained nonlinear time constant parameters are estimated by
the SQP so that the constrained nonlinear system is converted to a constrained linear
system. And then, the constrained linear connectivity parameters are estimated by the
constrained linear least squares (CLLS) method using a projected Hessian. The coupled
optimization is repeated until convergence. The experimental results show that HCNO
outperforms the base SQP method in terms of search time and prediction accuracy, i.e.,
faster convergence to better solutions.
xii
Second, we developed a factor graph-based probabilistic inference approach. The pro-
posed method is based on greedy forward search and belief propagation in a factor graph
for the learning of the graph structure and parameter estimation. To initiate the node
and edge search, an initial factor graph is created by introducing an inductive bias, the
locality principle. The graph structure is learned on the fly guided by belief discrepan-
cies, estimation error, and correlation analysis iteratively. At each round, the partially
constructed graph is the input for the learning process at the next round. The observed
data is exploited as evidences by local propagation. To reduce the estimation error caused
by the discretization of the continuous search spaces for the CRM parameters, the ex-
pected values are computed for the parameters. This learning process is repeated until
a stopping criterion is satisfied. Experimental results show that the elapsed time for the
message passing and the node and edge search is less than the search time of HCNO
for a large number of parameters, with the prediction accuracy comparable to that of
the optimization method. It is shown that the messages converge to a stable equilibrium
as time goes on in the loopy belief probabilistic system as advocated in [46] and the
approximation results in the high prediction accuracy.
xiii
Chapter 1
Introduction
1.1 Motivation
A major task of every oil company is to optimize field performance, i.e., maximizing oil
production and reducing operational cost. A difficulty associated with the optimal oper-
ation of an oil field is that its underlying structure is unknown and continuously changes
over time. Due to this uncertainty, field engineers have been struggling to optimize oil
recovery. One of the popular oil recovery techniques is waterflooding, which injects wa-
ter into the oil reservoir to extract oil. Much research has been done in investigating
IPRs. An Extended Kalman Filter (EKF) was used to forecast IPRs between multiple
injectors and a single producer [38]. IPRs were assumed to be non-stationary, i.e., IPRs
may change constantly due to changing bottom hole pressures, work-overs and natural
and/or man-made induced geomechanical effects. Their forecasts of IPRs are nonlinear
and adaptive functions of the data, where the nature of the nonlinearities is established by
the mathematics of the EKF. Production was predicted by the use of neural network (NN)
approach as a time series tool to represent the field as a black box [52]. NN was applied to
1
the data of mature field, where numerical reservoir simulation setup can be very complex
requiring several months to build a flow model. A probabilistic approach was proposed
using uncertainty quantification through the history matching process [4]. Probabilis-
tic multi-realizations of the reservoir model took into account static uncertainty using
JACTA, a software used to deal with geological uncertainties and determine OOIP eval-
uation. The impact of static uncertainties on the history match was calculated by means
of reservoir simulation and the impact of dynamic uncertainty on the history match was
calculated using experimental designs. Evolutionary algorithms were applied to history
matching [54]. Evolutionary strategies used continuous and discrete parameters. Continu-
ous parameters are permeabilities, transmissibilities, and relative permeabilities. Discrete
parameters represent different realizations of permeability maps or the placement of bar-
riers.
Recently, the capacitance-resistive model (CRM) has been employed to predict produc-
tion rates given water injection rates [49][50][56][57]. The CRM is a nonlinear model that
predicts production rates given water injection rates. The CRM consists of two parameter
types that represent the relationship of each injector-producer pair: one quantifies the
connectivity made by fractures or high permeability layers, and the other represents a time
constant that quantifies the degree of fluid storage between two wells. Three constraint
types are imposed on the two parameters: (1) the sum of the connectivity parameters
of an injector to all other producers should be 1; (2) individual connectivity parameters
π should fall between 0 and 1 (this implies a nonnegativity constraint as well); and (3)
the nonnegativity constraint that the time constant τ parameters should be positive real
2
numbers. The parameter estimation of CRM is expected to identify and quantify the
injector-producer (interwell) connectivity of multiple injectors-multiple producers and to
provide production estimation. The challenge is that an oil reservoir consists of hun-
dreds of injection and production wells, thus the number of CRM parameters increases
quadratically. Thus, parameter estimation for the CRM is a large-scale parameter esti-
mation problem for constrained nonlinear systems. In this research, we have developed
two different methods for the parameter estimation problem by mapping it to constrained
nonlinear optimization problem and inference in probabilistic graphical models.
Unfortunately, that there are no effective methods for solving a general nonlinear pro-
gramming problem. Estimating as few as ten variables can be extremely challenging,
so estimating hundreds of variables can be intractable. Therefore, several different ap-
proaches have been taken for the general nonlinear programming problem [45][8]. Global
optimization methods are used for problems with a small number of variables in which
search time is not critical. The goal of global optimization is to find the true global solu-
tion, thus the compromise is efficiency. The worst-case complexity of global optimization
methods grow exponentially with the problem sizes. While the true global solution is
found, it is not typical. Global optimization is used for small-scale problems only if it is
highly possible to find the true global solution. Local optimization methods do not always
find the global solution while they seek only a local solution at which the objective func-
tion is smaller than at all other feasible points in its vicinity. Local optimization does not
seek the optimal solution, which optimizes the objective over all feasible points. Instead,
a locally optimal point is sought, which optimizes the objective function among feasible
3
points. A large portion of the research on general nonlinear programming has focused on
local optimization methods. In consequence, they are well developed. In addition to the
fact that the true globally optimal solution is not found, there are several disadvantages
of local optimization methods. The methods require an initial guess for the optimization
variables. Little information is known about how far from globally optimal solution the
local solution is. Local optimization methods are sensitive to algorithm parameter values,
which need to be adjusted for a particular problem. In spite of those disadvantages, the
methods are widely applicable since they only require differentiability of the objective
and constraint functions, the methods are fast, can handle large-scale problems. Hence,
we developed a hybrid constrained nonlinear optimization (HCNO) method, which is a
gradient-based method for large-scale constrained nonlinear optimization [32]. HCNO is
on the basis of sequential quadratic programming in line search framework. The exper-
imental results showed that the method outperformed the conventional SQP, i.e., faster
convergence to better solutions. However, there are two major shortcomings in such con-
tinuous nonlinear optimization methods. First, the performance of those methods highly
depends on the initial guess due to the presence of multiple local optima as already men-
tioned. Second, they assign values to all parameters even if not every parameter needs to
be estimated. Thus, its computational cost becomes unnecessarily high.
To address these problems, we developed a factor graph-based probabilistic inference
approach to estimate the CRM parameters called PETROGRAPH Learning (PGL). A
factor graph represents a factorization of a complicated global function so that the fac-
torization can be exploited by algorithms that deal with the global function. Due to the
4
unknown IPRs, estimating parameters of the CRM is mapped to learning both the graph
structure and parameters in the graph from observable data. The search space of true
factor graphs is exponential, thus an exhaustive enumeration of all the possible factor
graphs is infeasible. The inputs to the problem are a set of injection rates over time, a set
of production rates over time, joint probabilities (JPs), and discrete search space. Given
the inputs, an inductive bias is introduced to construct an initial factor graph. The initial
factor graph is created based on the locality principle such that a producer is influenced
by its closest injector. Message passing is started on the initial factor graph, and then
nodes and edges are added to the graph at every round by investigating the discrepancy
among messages, estimation error, and the correlation of injection rates and the residual
production rates between the observed and estimated production rates. The observed data
is exploited as evidences by local propagation. To reduce the estimation error caused by
the discretization of the continuous search spaces for the CRM parameters, the expected
values are computed for the estimated parameters. This procedure is repeated until a
stopping criterion is satisfied. The outputs are a factor graph and the estimated values
for the connectivity and the time constant parameters of the CRM.
1.2 Contributions
The contributions of this thesis are the development of two different methods for the
large-scale parameter estimation for the constrained nonlinear system, CRM, using two
different approaches, a hybrid constrained nonlinear optimization algorithm based on
SQP in line search framework and structure learning of factor graph and inference in the
5
probabilistic graphical model.
• A Hybrid constrained nonlinear optimization (HCNO)
The nonlinear problem is decomposed into constrained nonlinear and constrained
linear optimization problems. A constrained nonlinear optimization algorithm is ap-
plied only to the estimation of the nonlinear parameters, and a convex optimization
algorithm is applied to estimating the linear parameters. SQP with a line search
framework is applied to the constrained nonlinear parameter estimation problem
such that the constrained nonlinear parameters are estimated first so that the non-
linear system is converted to a linear system, then the constrained linear parameters
are solved. To the best of our knowledge, no work has been done on applying the cou-
pled parameter estimation using SQP with a line search framework to constrained
nonlinear parameter estimation problems. HCNO is tested on synthetic oil field
data to estimate the parameters of a capacitance-resistive model, where the result
shows improvements in both search time and prediction accuracy.
• Structure learning of factor graph and inference in the probabilistic graphical model
guided by belief discrepancies, estimation error, and residual correlation analysis
The inputs to the problem are a set of injection rates over time, a set of production
6
rates over time, joint probabilities (JPs), and discrete search space. Given the in-
puts, an inductive bias is introduced to construct an initial factor graph. The initial
factor graph is created based on the locality principle such that a producer is influ-
enced by its closest injector. Message passing is started on the initial factor graph,
and then nodes and edges are added to the graph at every round by investigating
the discrepancy among messages, estimation error, and the correlation of injection
rates and the residual production rates between the observed and estimated pro-
duction rates. The observed data is exploited as evidences by local propagation. To
reduce the estimation error caused by the discretization of the continuous search
spaces for the CRM parameters, the expected values are computed for the estimated
parameters. This procedure is repeated until a stopping criterion is satisfied. The
outputs are a factor graph and the estimated values for the connectivity and the
time constant parameters of the CRM. The experimental results show that running
time is reduced compared with that of HCNO the large number of parameters and
the prediction accuracy is comparable to that of HCNO. We design a probabilistic
graphical model and demonstrate that messages converge to a stable equilibrium in
the loopy probabilistic system.
7
Chapter 2
Capacitance-Resistive Model and Problem Formulation
The capacitance-resistive model (CRM) relies on signal processing techniques, in which
the oil reservoir is viewed as a system that converts input signals into output signals
[49][50][56][57]. Water injection rates are treated as input signals, and oil production
rates are the reservoir response or output signals. Interwell connectivity and response
delay constitute the unknown system parameters, and production rates can be predicted
by fitting the unknown system parameters to the past observations. CRM models the
reservoir flow behavior in accordance with interactions between well pairs. For a given
CRM and a time-series dataset of injection and production rates, the estimated total
production rate of producer j at time t
n
is
ˆ q
j
(t
n
)=
L
i=1
n
k=1
(1− e
−Δt
k
/τ
ij
)π
ij
i
i
(k)e
−(tn−t
k
)/τ
ij
, (2.1)
where L is the number of injectors, Δt
k
is the time interval between t
k
and t
k−1
, π
ij
is
the fraction of injection rate of injector i flowing toward producer j, i
i
(k) is the injection
rate of injector i at time step k,and τ
ij
is a time constant. Our goal is to estimate values
8
of the two sets of parameters, π
ij
and τ
ij
(i = 1,...,L and j = 1,...,J), which minimize the
prediction error between the actual and the estimated production rates (J is the number
of producers). Thus, the objective function is defined as
f(x)=
J
j=1
N
n=1
(q
j
(t
n
)− ˆ q
j
(t
n
))
2
, (2.2)
where x =(π
11
,π
12
,...,π
LJ
,τ
11
,τ
12
,...,τ
LJ
), q
j
(t
n
) is the actual production rate of pro-
ducer j at time t
n
and N is the total number of time steps. The sum of square errors is
used as the measure of training error. Now, a constrained nonlinear optimization problem
can be formulated as follows: Minimize f (x) subject to constraints,
J
j=1
π
ij
=1, (2.3)
π
lb
≤ π
ij
≤ π
ub
, (2.4)
τ
lb
≤ τ
ij
≤ τ
ub
, (2.5)
where i = 1,2,...,L, j = 1,2,...,J, π
lb
=0 and π
ub
= 1. The objective is to minimize the
prediction error, Eq. 2.2, subject to the equality constraint, Eq. 2.3, and the upper and
lower bounds, Eqs. 2.4 and 2.5. Although the exact lower and upper bounds of τ
ij
are not
readily available, they can be approximated [49][50]. In general, τ
lb
is supposed to take a
very small value and τ
ub
is a reasonably large value. In our experiment, τ
lb
and τ
ub
are
set to 0.0001 and 100, respectively. By doing so, the bounds in Eqs. 2.4 and 2.5 provide
a feasible search region. The CRM parameters are evaluated based on only injection and
production history. Once the model parameters are estimated, the IPRs can be estimated
9
using the estimated model parameters and the prediction of production rates can be
made.
10
Chapter 3
Hybrid Constrained Nonlinear Optimization
3.1 Related Work
Nonlinear optimization is challenging due to the nonconvex property of the objective
function and/or constraints. Identifying a global optimum is extremely hard even with-
out constraints. The curse of dimensionality makes the problem more difficult to solve: as
the number of parameters increases, the number of local optima is likely to increase. No
universal algorithm is available for this problem, but there are numerous algorithms
broadly categorized into gradient-based methods and direct search methods. Gradient-
based methods use first or second derivatives of the objective function with respect to
the parameters. Examples of the gradient-based algorithms are Newtons method, quasi-
Newton methods, and sequential quadratic programming (SQP) [8][45][18][5]. On the
other hand, direct search methods do not rely on derivatives. Examples include the
Nelder-Mead method, genetic algorithms, differential evolution, simulated annealing, and
11
the particle swarm algorithm. This paper focuses on SQP methods for nonlinear optimiza-
tion problems with linear constraints. SQP is one of the most popular and robust algo-
rithms for constrained nonlinear continuous optimization. At each iteration, SQP solves a
quadratic programming (QP) subproblem subject to a linearization of the constraints and
defines the search direction as the solution of the QP subproblem. The objective function
of the subproblem is a quadratic approximation to the Lagrangian function. When SQP
algorithms employ the exact Hessian, the convergence rate is quadratic. When a quasi-
Newton approximation to the Hessian is used, the convergence rate is superlinear. SQP
algorithms have two shortcomings. First, they require the QP subproblem to be solvable
per iteration. Second, the solutions of the sequential QP subproblems may be unbounded,
which causes the sequence to diverge. Despite those drawbacks, SQP algorithms have been
used to solve a remarkably large number of important practical problems and are very
effective for solving small or medium-size nonlinear problems. When the solution to the
QP subproblem is used to compute the search direction, most modern SQP algorithms
employs active-set methods. The number of iterations required to solve QP subproblems
using an active-set method grows quickly with the number of parameters. A variety of
work has been done to reduce the SQP or QP iterations in applying SQP methods to solv-
ing large-scale constrained nonlinear optimization problems. The QP subproblems were
solved by an interior point method that could be prematurely halted by a trust region
constraint so as to achieve the overall effectiveness of the procedure [6]. A semidefinite
QP solver was utilized, which was based on a limited-memory quasi-Newton approx-
imation to the Hessian of the Lagrangian. A reduced-Hessian algorithm was used for
12
solving the QP subproblems [42][21]. For the necessity to deviate from the standard ap-
proach when solving large problems, more flexible SQP methods were analyzed by using
a particular augmented Lagrangian merit function but without requiring a minimizer of
the QP subproblems [43]. A SQP algorithm generated descent directions for an l
1
plus
log-barrier merit function and used a line-search to obtain a sufficient decrease of this
function [48]. The QP subproblems were solved by an interior-point method using an
inexact Newton approach, iterating to an accuracy just sufficient to produce a descent di-
rection in the early stages and tightening the accuracy as a solution was approached. For
practical performance in large-scale problems, the progress of the iterations was assessed
by a constrained merit function [7]. Since the exact computation of steps in a SQP can be
prohibitively expensive for very large problems, an exact penalty function was employed
to measure sufficient progress made by inexact step toward a solution of the nonlinear
program [10]. The QP subproblems were solved by an interior-point method only ap-
proximately by terminating the inner iterative process during the first early stage and
tightening the accuracy of the solutions of the quadratic problems as the iteration pro-
gresses [35]. These SQP algorithms estimated all the parameters as a whole solely by
con-strained nonlinear optimization, which incurs enormous complexity for large-scale
nonlinear problems with linear constraints. Many applications, e.g., CRM, however, can
be decomposed into subproblems, some of which are linear. In such cases, it is possible
to reduce the complexity by utilizing linear optimization methods.
13
3.2 Base Algorithm
The SQP algorithm is implemented in a line search framework with quasi-Newton Hes-
sian approximations used for the step computation. For the implementation, we utilize
constrained nonlinear optimization (fmincon) and constrained linear least squares (lsqlin)
functions in Matlab. We first present a base constrained nonlinear optimization algorithm
followed by the HCNO algorithm. Our base constrained nonlinear algorithm uses SQP
with a quasi-Newton method in a line search framework. At each iteration, the SQP
solves a QP subproblem whose solutions are a vector of search directions d
k
and a vec-
tor of Lagrange multipliers λ
k
for solution x
k
. The first step is to represent the equality
constraints and bounds, Eqs. 2.3-2.5 as the Jacobian matrix A of constraints
a
T
m
x = b
m
(m∈ ξ), (3.1)
a
T
m
≥ b
m
(m∈ ζ), (3.2)
where ξ and ζ are finite sets of indices, m, to which equality and inequality constraints
on a
m
apply, respectively. The bounds are converted to a set of inequality constraints
[17]
Ix≥ B
l
, (3.3)
−Ix≥−B
u
, (3.4)
14
where I is the identity matrix. The Jacobian matrices for all the constraints are con-
structed as
A =
⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ A
eq
I
−I
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , (3.5)
B =
⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ B
eq
B
l
−B
u
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , (3.6)
where A
eq
is the Jacobian for the equality constraints, B
eq
=(1,...,1)
T
, B
l
=(π
lb
,...,τ
lb
,
...)
T
andB
b
=(−π
ub
,...,−τ
ub
,...)
T
. The QP subproblem has the form: minimize
1
2
p
T
Gp+
g
T
p subject to
c
T
m
x
k
p + c
T
m
x
k
=0 (m∈ ξ), (3.7)
c
T
m
x
k
p + c
T
m
x
k
≥0(m∈ ζ), (3.8)
where p is a search direction, G is the approximated Hessian matrix of Lagrangian func-
tion L(x,λ), g is the function gradient, and c
T
m
x
k
= a
T
m
− b
m
. Having formulated the
QP subproblem, we now need to find a solution that satisfies the first-order optimality
conditions, Karush-Kuhn-Tucker (KKT) conditions. The projection-based active-set al-
gorithm constructs an active set W from the constraint set at each iteration [45][18]. For
15
the Hessian update, the quasi-Newton approximation of the Hessian of the Lagrangian
function at the k
th
iteration is
H
k+1
= H
k
−
H
k
s
k
s
T
k
H
k
s
T
k
H
k
s
k
+
y
k
y
T
k
y
T
k
s
k
, (3.9)
where s
k
= x
k+1
- x
k
and y
k
=
x
L(x
k+1
,λ
k+1
)−
x
L(x
k
,λ
k
). Specifically, we use the
Broyden-Fletcher-Goldfarb-Shanno (BFGS) method for the update [45][9]. The BFGS
update retains positive-definiteness of the approximated Hessian to prevent QP subprob-
lems from being poorly posed. For the progress evaluation due to the constraints, a merit
function is employed to measure the right balance of the two conflicting goals, i.e., reduc-
ing the objective function and satisfying the constraints. The merit function [47] is used
as a penalty function in finding the appropriate step size α
k
at every iteration such that
the merit function is minimized:
ψ(x)= f(x)+
me
i=1
r
i
c
i
(x)+
m
i=me+1
r
i
max{0,c
i
(x)}, (3.10)
where m
e
is the last index of the equality constraints, m is the last index of the inequality
constraints, and r
i
is the penalty parameter. Finally, the solution to the QP subproblem
gives a vector p
k
at the k
th
iteration and a new iteration is made as
x
k+1
= x
k
+ α
k
p
k
, (3.11)
where α
k
is a step length. The pseudocode of the base constrained nonlinear optimization
is given in Algorithm 1.
16
Algorithm 1. Base−caseconstrainednonlinearoptimization
1. Input :
2. Data: injection and production rates
3. init π, init τ: initial guess
4. Aeq, Beq: equality constraints
5. B π, B τ: bounds
6. Output :
7. est π new:estimateof π
8. est τ new:estimateof τ
9. Begin
10. est π old←{}
11. est τ old←{}
12. est π new ← init π
13. est τ new ← init τ
14. [A τ,B τ]← construct constraints(B τ)
15. iter ← 0
16. while∼done
17. g new ← computeGradient(Data, est π old, est τ old, est π new, est τ new)
18. s ← est τ new−est τ old
19. y ← g new−g old
20. H ← updateHessian(H , s, y)
21. est π old ← est π new
17
22. est τ old ← est τ new
23. g old ← g new
24. [λ τ,d τ,W τ]← solve QPsub(itest π old, est τ old, A τ, B τ, W τ, H )
25. meritFcnDone ← 0
26. while∼ itmeritFcnDone
27. α← updateStepSizebyMeritFcn(est π old,est τ old,α)
28. est π new ← est π old + α∗d
29. est τ new ← est τ old + α∗d
30. meritFcnDone ← checkMeritFncStoppingCriteria(Data,est π new,est τ new)
31. end while
32. iter ← iter +1
33. done ← checkStoppingCriteria()
34. end while
In Algorithm 1, computeGradient() computes the function gradient, updateHessian()
updates the approximation to the Hessian of the Lagrangian function, and solve QPsub()
solves QP subproblems producing the Lagrangian multipliers (λ), search direction (d)
and updated working set of constraints (W ). The function updateStepSizebyMeritFcn()
updates the step size at each iteration of the merit function, checkMeritFncStoppingCri-
teria() checks the stopping criteria of the merit function and checkStoppingCriteria()
checks any of the following stopping criteria: (1) the function gradient is less than a
tolerance; (2) the magnitude of search direction is less than a tolerance; (3) the difference
of two solutions of two consecutive iterations is less than a tolerance; (4) the first-order
18
optimality measure is less than a tolerance; and (5) maximum number of iterations is
exceeded.
3.3 Hybrid Constrained Nonlinear Optimization
In contrast to the base algorithm in section 3.1, HCNO is a two-step algorithm. Rather
than finding a solution for all the parameters together, the parameters are separated into
two types, linear and nonlinear. As we can see in Eq. (1), π
ij
are linear parameters and
τ
ij
are nonlinear parameters. The main idea of HCNO is to apply the base-case algorithm
in section 3.1 only to estimating the nonlinear parameters, whereas the linear parameters
are estimated by the constrained linear least squares (CLLS) method. Because of the
constraints, the linear parameter estimation problem cannot be solved by simple analyt-
ical formula. Therefore, coupled two-step optimization is performed per iteration of the
SQP algorithm. The algorithm works as follows. The nonlinear optimization minimizes
the objective function, f ( τ)( τ is the vector of τ parameters), subject to the bounds, Eq.
(2.5). The bounds are converted to a set of inequality constraints and the corresponding
QP subproblem is to minimize
1
2
subject to Eq. (3.8). The approximate Hessian is only
for the nonlinear parameters τ. With the resulting τ from the solution of the QP sub-
problem in the merit function evaluation, a constrained linear least squares problem is
then formulated as follows: Minimize
1
2
Cx −d
2
2
subject to
J
j=1
π
ij
=1, (3.12)
0≤ π
ij
≤ 1, (3.13)
19
where
C =
⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ q
1,1
(t
1
) ··· q
L,1
(t
1
)
.
.
.
q
L,1
(t
N
) ··· q
L,1
(t
N
)
q
1,2
(t
1
) ··· q
L,2
(t
1
)
.
.
.
q
1,2
(t
N
) ··· q
L,2
(t
N
)
.
.
.
q
1,J
(t
1
) ··· q
L,J
(t
1
)
.
.
.
q
1,J
(t
N
) ··· q
L,J
(t
N
)
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,
q
i,j
(t
n
)=
n
k=1
(1-e
−t
k
/τ
ij
)i(k)e
−(tn−t
k
)/τ
ij
,
d =(q
1
(t
1
),q
1
(t
2
),...,q
1
(t
n
),q
2
(t
1
),...,q
2
(t
N
),...,q
J
(t
1
),...,q
J
(t
N
))
T
.
The nonlinear parameters τ
ij
are now replaced by the estimated values from the QP sub-
problem, so that estimating π
ij
is a linear least squares problem with the reduced number
of parameters by half subject to the constraints. The problem can be solved by quadratic
programming which handles the constraints. Specifically, we adopt quadratic program-
ming using a projected Hessian as in [18]. The CLLS problem is simply to minimize the
2-norm of the residual of the objective function w.r.t. subject to the linear constraints,
Eqs. (3.12) and (3.13). The bounds are converted to a set of inequality constraints by a
similar method used in the base algorithm, so that a constraint matrix is constructed. The
quadratic problem for the CLLS is to minimize 0.5x
T
Gx +c
T
x subject to Eqs. (3.1) and
(3.2), where G = C
T
C and c = C
T
d. The pseudocode of HCNO is given Algorithm2,
20
where solve conlinlsq() is an iterative optimization function.
Algorithm2. HCNO
1. Input :
2. Data: injection and production rates
3. init π, init τ: initial guess
4. Aeq, Beq: equality constraints
5. B π, B τ: bounds
6. Output :
7. est π new:estimateof π
8. est τ new:estimateof τ
9. Begin
10. est π old←{}
11. est τ old←{}
12. est π new ← init π
13. est τ new ← init τ
14. [A τ,B τ]← construct constraints(B τ)
15. iter ← 0
16. while∼done
17. g new ← computeGradient(Data, est π old, est τ old, est π new, est τ new)
18. s ← est τ new−est τ old
19. y ← g new−g old
20. H ← updateHessian(H , s, y)
21
21. est τ old ← est τ new
22. g old ← g new
23. [λ τ,d τ,W τ]← solve QPsub(est τ old, A τ, B τ, W τ, H )
24. meritFcnDone ← 0
25. while∼ itmeritFcnDone
26. α← updateStepSizebyMeritFcn(est τ old, α τ)
27. est τ new ← est τ old + α τ∗d τ
28. [A π, B π]← construct constraints(Aeq, Beq, B π)
29. [est π new, λ π]← solve conlinlsq(est π new, est τ new, A π, B π)
30. meritFcnDone ← checkMeritFncStoppingCriteria(Data, est π new,
est τ new, W τ)
31. end while
32. iter ← iter +1
33. done ← checkStoppingCriteria()
34. end while
3.4 Performance Evaluation
HCNO was compared with two other algorithms, the base algorithm and the constrained
genetic algorithm (cGA). We used the genetic algorithm toolbox in Matlab (ga)for the
constrained genetic algorithm. Performance of the two algorithms was measured in terms
of prediction accuracy, search time, and locating connectivity. Real oil-field data is not
appropriate for this purpose, due to the unavailability of the underlying structure, e.g.,
22
locations of high permeability. The prediction accuracy of the production rates does not
guarantee the accuracy of the estimated values of the parameters that are used to infer
the interwell connectivity. Therefore, numerically simulated synthetic field data was used
for the evaluation, which provided the full geophysical information.
3.4.1 Experimental Setting
Data
CMG simulator, a commercial oil field simulation software, was used to generate the data
for four fields. Table 3.1 shows the information about the four synthetic fields ordered by
the ascending number of wells. Figure 3.1 shows the maps of the fields with the locations
of fractures, where (red) dots stand for injectors, (black) squares for producers, (grey)
bars for locations of higher permeability, and the x-y axes are the relative x-y coordinates
of wells from the origin (0,0) on the grid. Figure 3.2 shows an example of injection and
Table 3.1: Simulated synthetic oil fields: Field name, field type, field size, number of time
steps, number of parameters, number of bounds and number of equality constraints
Field name Type Size Time steps Parameters Bounds Equality
Field1 Homogeneous 251
Field2 Heterogeneous
5× 4
3041
40 80 5
Field3 Homogeneous 160
Field4 Heterogeneous
8× 8
512
128 256 8
production signals for Field1. As we can see, the injection rates are pulse-like signals. Fig-
ure 3.3 shows an example of the injection and production signals for Field2, Field3 and
23
Field4 in which all the signals have high variation instead of relatively simple pulse-like
signals.
Search Time and Test Error
For Field1 and Field2 a set of 20 initial guesses were randomly generated for the base
algorithm and HCNO, and this same set of initial guesses was used for all Field1 and
Field2 test cases. Another set of 20 initial guesses was used for the Field3 and Field4 test
cases. For cGA, 20 initial populations were randomly generated for Field1 and Field2, and
another set of 20 initial populations was generated for Field3 and Field4. Of the injection
and production data, 80% was used for the training data, and the remaining 20% of the
data was used to test prediction accuracy. For each field, the same initial guesses were
used for the base algorithm and HCNO. The two sets of initial guesses were also used
for HCNO for the comparison with cGA. All runs were made on a PC with a 2.6GHz
Intel Core2 Duo CPU and 2GB of RAM. For the statistical significance test, a Wilcoxon
signed-rank test for paired samples was performed, since the distributions of the results
were not normal but multimodal.
Number of Iterations and Population Size
Table 3.2 shows the optimization settings for the maximum numbers of iterations per
test case for the comparison of the base algorithm versus HCNO. The number of major
iterations is the number of iterations of the outmost loop to achieve convergence. The
number of minor iterations is the number of iteration to solve a QP subproblem. We have
found that the maximum number of minor iterations has negligible effect on the number
24
of major iterations, search time and test error. Thus for each test case, it was set to the
number of constraints for the field. MinorLsq is the number of iterations to solve a QP
problem for the constrained linear least squares and was set to the number of constraints
for the linear parameters. The maximum number of iterations was determined by follow-
ing rule. As mentioned, the maximum numbers of minor iterations and MinorLsq were
set to the number of constraints. The number of major iterations affected the quality of
solution. Because HCNO converged after a certain number of major iterations, the max-
imum number of major iterations of HCNO could be found after multiple trials. Once
the number of major iterations for HCNO was found, it was used as the baseline for the
subsequent test cases. The observation from the experiment was that it was very time-
demanding to find the number of major iterations for the base algorithm to converge, so
we stopped at reasonable point per field. Table 3.3 shows the experimental settings for
Table 3.2: Maximum number of iterations for Base vs HCNO
Field name No
Base HCNO
Major Minor Major Minor MinorLsq
Field1
110
85 10 40 45
220
350
4 100
Field2
110
85 10 40 45 250
3 100
Field3
1 200
264 200 128 136
2 500
Field4
1 200
264 200 128 136 2 500
2 1000
25
the comparison with cGA. The way that the number of generations and population size
for cGA were determined was by exploring their search space as we did in determining
the number of major iterations. The observation from the exhaustive search was that the
number of generations did not significantly contribute to the quality of solution while
the population size did. Thus, we tried multiple test cases by gradually increasing the
population size with a minimum of at least the number of parameters to be estimated.
In Table 3.3, the number of generations of each test case is 100.
Table 3.3: Population size of cGA
Field name No Population size
Field1
150
2 100
3 150
Field2
1 100
2 200
Field3
1 200
2 300
Field4
1 200
2 300
26
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(a) 5× 4 homogeneous field
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(b) 5× 4 heterogeneous field
0 200 400 600 800 1000 1200 1400 1600
0
200
400
600
800
1000
1200
1400
1600 Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(c) 8× 8 homogeneous field
0 200 400 600 800 1000 1200 1400 1600
0
200
400
600
800
1000
1200
1400
1600 Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(d) 8× 8 heterogeneous field
Figure 3.1: Well locations and true higher permeability
27
0 50 100 150 200 250 300
0
500
1000
1500
2000
Time Step
Injection Rate
(a) Injection rate
0 50 100 150 200 250 300
800
900
1000
1100
1200
1300
1400
1500
Time Step
Production Rate
(b) Production rate
Figure 3.2: Example of the injection and production signals of Field1
0 500 1000 1500 2000 2500 3000 3500
0
500
1000
1500
2000
2500
3000
3500
Time Step
Injection Rate
(a) Injection rate
0 500 1000 1500 2000 2500 3000 3500
200
300
400
500
600
700
Time Step
Production Rate
(b) Production rate
Figure 3.3: Example of the injection and production signals of Field2, 3, and 4
28
3.4.2 Results
Table 3.4 and 3.5 show the number of converged runs and means of the number of major
iteration, search time, and residual error of the test data. In Table 3.4, convergence is
the number of converged runs that satisfied any convergence criterion before reaching
the maximum number of major iterations. In Table 3.5, the root mean squared error
(RMSE) was measured to show the error with the same magnitude as the fluid quantity
being predicted because the magnitudes of the production rates of all the fields are not
exactly the same. As additional information, the root relative squared error (RRSE) is
measured to show the error relative to the training data and of the same magnitude as
the quantity predicted. Table 3.6 shows the result of the statistical significance test for
number of major iterations, search time, RMSE and RRSE. Negative rank is the number
of cases in which the HCNO result is less than the base result, and P indicates the
significance value. The four key findings are as follows.
• HCNO took less time to converge than the base algorithm did. The means of the
number of converged runs and search time in Table 3.5 show that HCNO converged
for most of the test cases. The significance test of search time confirms that the
search time of HCNO is significantly less than that of the base algorithm (Table
3.6). The base algorithm spent more time but failed to converge most of the time. It
turned out that CLLS contributed to the faster convergence of HCNO.
• HCNO converged to better solutions than the base algorithm did. RMSE and RRSE
of HCNO are significantly less than that of the base algorithm (Table 3.6), which
was confirmed by the significance test in Table 3.6. It implies that the solutions
29
Table 3.4: Number of convergences, means of number of major iterations and search time
Field name No
Convergence Mean
Major iterations Search time (sec)
Base HCNO Base HCNO Base HCNO
Field1
1
020
10
1.7
2.08
0.77
220 3.94
3 50 10.31
4 100 21.50
Field2
1
018
10
2.5
5.30
19.68 2 50 25.51
3 100 49.86
Field3
1
120
191.2
102.85
213.05
262.49
2 476.2 588.45
Field4
10
20
200
139.25
267.92
1137.53 2 0 500 694.15
3 20 915.85 1200.62
found by HCNO were local optima which were closer to global optima than those
found by the base algorithm.
• As the number of parameters increased, both algorithms needed more major itera-
tions and search time to converge. The means of major iterations and search time of
Field1 and Field2 were less than those of Field3 and Field4 (Table 3.4). This obser-
vation demonstrates the curse of dimensionality in optimization. However, HCNO
converged after spending less time than the base algorithm did.
• The search time and test error of HCNO are significantly less than those of cGA. Be-
cause the fundamental algorithms of HCNO and cGA are different, the same analysis
30
Table 3.5: Means of root mean squared error (RMSE) and root relative squared error
(RRSE)
Field name No
RMSE RRSE
Base HCNO BASE HCNO
Field1
1 1696.39
139.66
86.83
1.11
2 981.86 49.41
3 240.08 9.27
4 141.81 4.43
Field2
1 52743.65
22.63
42.19
0.30 2 8191.93 6.55
3 1391.29 1.11
Field3
1 4310.34
78.09
3.95
0.16
2 3369.99 3.08
Field4
1 3130.89
27.58
1.31
0.04 2 237.51 0.10
3 107.59 0.04
as we did for Base vs HCNO is inappropriate. Thus, general criteria such as search
time and test error are used. As Table 3.7 shows, the search time and test error
of cGA are significantly higher than those of HCNO. An interesting observation is
that the search time of cGA for Field4 is significantly less than that of HCNO but
the test error indicates that the quality of solution of cGA is significantly lower
than that of HCNO and the increase of population size did not reduce the error. An
exceptional observation is that the increase of population size for Field2 resulted in
a slightly higher test error for cGA.
In addition to those key findings, we make the following observations. The maximum
number of minor iterations did not make a significant difference in reducing the number
31
of major iterations, search time and test error. Several tests were done by increasing the
maximum number of minor iterations to see if it reduced the number of major iterations,
search time and test error. The result showed that the number of minor iterations made a
negligible difference, i.e., 0.005 ∼ 0.01 seconds reduction in search time and no improve-
ment in the number of converged runs, number of major iterations and test error. In
addition to the cases presented in Table 3.4 ∼ 3.8, we tried an even larger field consist-
ing of 13 injectors and 12 producers, in which the number of parameters was 312, the
number of bounds was 624 and the number of equality constraints was 13. For this case,
the base algorithm failed by diverging for some initial guesses soon after a few major it-
erations. For some initial guesses, the algorithm quickly stopped by satisfying a stopping
criterion but the test error was too high. Though HCNO did produce a converged result
for this case, we have not included the result due to the unavailability of the base data for
comparison. Even though the number of parameters and constraints was not very large,
the significant degree of nonlinearity of the objective function and the complex nature of
the data made the optimization problem very hard. We have also noted that the injection
signals with high variation make the optimization harder than the pulse-like injection
signals. Figs. 3.4 ∼ 3.15 compare the observed production rates (solid curves) and the
predicted production rates (dotted curves) as a function of time for the run with the
least test errors among the 20 runs for each producer of Field1 ∼ Field4 respectively for
both base and HCNO algorithms. In Figs. 3.4 ∼ 3.15, vertical dotted lines separate the
training and test data sets. Overall, the production prediction curves by HCNO
fit the observed production rates better than the base algorithm.
32
0 50 100 150 200 250 300
0
500
1000
1500
2000
2500
3000
Test
Producer1
(a) Base
0 50 100 150 200 250 300
200
400
600
800
1000
1200
1400
1600
Test
Producer1
(b) HCNO
0 50 100 150 200 250 300
0
500
1000
1500
2000
2500
3000
Test
Producer2
(c) Base
0 50 100 150 200 250 300
200
400
600
800
1000
1200
1400
1600
Test
Producer2
(d) HCNO
Figure 3.4: Observed and predicted production curves of producer1 and 2 with the least
test errors of test case no. 1 of Base algorithm for Field1
33
0 50 100 150 200 250 300
−200
0
200
400
600
800
1000
1200
1400
1600
Test
Producer3
(a) Base
0 50 100 150 200 250 300
0
200
400
600
800
1000
1200
1400
1600
Test
Producer3
(b) HCNO
0 50 100 150 200 250 300
0
200
400
600
800
1000
1200
1400
1600
Test
Producer4
(c) Base
0 50 100 150 200 250 300
0
200
400
600
800
1000
1200
1400
1600
Test
Producer4
(d) HCNO
Figure 3.5: Observed and predicted production curves of producer3 and 4 with the least
test errors of test case no. 1 of HCNO for Field1
34
Table 3.6: Number of negative ranks (HCNO - Base) and significance by Wilcox signed-rank test for the number of major
iterations, search time, RMSE and RRSE
Field name No
Major iterations Search time (sec) RMSE RRSE
Negative ranks P Negative ranks P Negative ranks P Negative ranks P
Field1
120 <.05 20 <.05 20 <.05 20 <.05
220 <.05 20 <.05 20 <.05 20 <.05
320 <.05 20 <.05 16 <.05 20 <.05
420 <.05 20 <.05 8 >.05 20 <.05
Field2
118 <.05 0 <.05 20 <.05 20 <.05
220 <.05 14 <.05 20 <.05 20 <.05
320 <.05 20 <.05 20 <.05 19 <.05
Field3
120 <.05 4 >.05 20 <.05 20 <.05
220 <.05 20 <.05 20 <.05 20 <.05
Field4
120 <.05 1 <.05 20 <.05 20 <.05
220 <.05 1 <.05 20 <.05 19 <.05
320 <.05 11 >.05 19 <.05 18 <.05
35
Table 3.7: Means of search time, root mean squared error (RMSE), and root relative
squared error (RRSE) of cGA
Field name No
Search time (sec) RMSE RRSE
cGA HCNO cGA HCNO cGA HCNO
Field1
1 12.18
0.77
2923.90
139.66
21.06
1.11 2 34.15 1927.35 13.88
3 101.60 521.89 3.75
Field2
1 139.19
19.68
20522.11
22.63
16.41
0.30 2 420.11 23314.33 18.65
Field3
1 230.85
262.49
3770.22
178.09
9.76
0.16
2 678.90 3759.33 9.73
Field4
1 333.63
1137.53
10728.38
27.58
12.67
0.04
2 812.73 10611.31 12.53
36
Table 3.8: Number of negative ranks (HCNO - cGA) and significance by Wilcox signed-rank test for search time, RMSE and
RRSE
Field name No
Search time (sec) RMSE RRSE
Negative ranks P Negative ranks P Negative ranks P
Field1
120 <.05 20 <.05 20 <.05
220 <.05 18 <.05 18 <.05
320 <.05 13 <.05 13 <.05
Field2
120 <.05 20 <.05 20 <.05
220 <.05 20 <.05 20 <.05
Field3
16 >.05 20 <.05 20 <.05
220 <.05 20 <.05 20 <.05
Field4
11 <.05 20 <.05 20 <.05
24 <.05 20 <.05 20 <.05
37
0 500 1000 1500 2000 2500 3000 3500
0
1000
2000
3000
4000
5000
Test
Producer1
(a) Base
0 500 1000 1500 2000 2500 3000 3500
1000
1500
2000
2500
3000
3500
4000
4500
5000
Test
Producer1
(b) HCNO
0 500 1000 1500 2000 2500 3000 3500
0
200
400
600
800
1000
1200
1400
1600
1800
Test
Producer2
(c) Base
0 500 1000 1500 2000 2500 3000 3500
50
100
150
200
250
300
350
400
450
500
Test
Producer2
(d) HCNO
Figure 3.6: Observed and predicted production curves of producer1 and 2 with the least
test errors of test case no. 1 of Base algorithm for Field2
38
0 500 1000 1500 2000 2500 3000 3500
0
500
1000
1500
2000
2500
3000
Test
Producer3
(a) Base
0 500 1000 1500 2000 2500 3000 3500
200
300
400
500
600
700
Test
Producer3
(b) HCNO
0 500 1000 1500 2000 2500 3000 3500
0
500
1000
1500
2000
2500
3000
3500
4000
4500
Test
Producer4
(c) Base
0 500 1000 1500 2000 2500 3000 3500
1000
1500
2000
2500
3000
3500
4000
4500
Test
Producer4
(d) HCNO
Figure 3.7: Observed and predicted production curves of producer3 and 4 with the least
test errors of test case no. 1 of HCNO for Field2
39
0 50 100 150 200
400
600
800
1000
1200
1400
1600
1800
Test
Producer1
(a) Base
0 50 100 150 200
400
600
800
1000
1200
1400
1600
1800
Test
Producer1
(b) HCNO
0 50 100 150 200
100
200
300
400
500
600
700
800
900
Test
Producer2
(c) Base
0 50 100 150 200
200
300
400
500
600
700
800
900
Test
Producer2
(d) HCNO
Figure 3.8: Observed and predicted production curves of producer1 and 2 with the least
test errors of test case no. 1 of Base algorithm for Field3
40
0 50 100 150 200
500
1000
1500
2000
Test
Producer3
(a) Base
0 50 100 150 200
600
800
1000
1200
1400
1600
1800
2000
2200
Test
Producer3
(b) HCNO
0 50 100 150 200
0
200
400
600
800
1000
1200
Test
Producer4
(c) Base
0 50 100 150 200
300
400
500
600
700
800
900
1000
1100
Test
Producer4
(d) HCNO
Figure 3.9: Observed and predicted production curves of producer3 and 4 with the least
test errors of test case no. 1 of Base algorithm for Field3
41
0 50 100 150 200
400
600
800
1000
1200
1400
1600
1800
2000
Test
Producer5
(a) Base
0 50 100 150 200
600
800
1000
1200
1400
1600
1800
2000
Test
Producer5
(b) HCNO
0 50 100 150 200
200
300
400
500
600
700
800
900
1000
1100
Test
Producer6
(c) Base
0 50 100 150 200
300
400
500
600
700
800
900
1000
1100
Test
Producer6
(d) HCNO
Figure 3.10: Observed and predicted production curves of producer5 and 6 with the least
test errors of test case no. 1 of Base algorithm for Field3
42
0 50 100 150 200
200
400
600
800
1000
1200
1400
1600
1800
Test
Producer7
(a) Base
0 50 100 150 200
400
600
800
1000
1200
1400
1600
1800
Test
Producer7
(b) HCNO
0 50 100 150 200
100
200
300
400
500
600
700
800
900
Test
Producer8
(c) Base
0 50 100 150 200
200
300
400
500
600
700
800
900
Test
Producer8
(d) HCNO
Figure 3.11: Observed and predicted production curves of producer7 and 8 with the least
test errors of test case no. 1 of Base algorithm for Field3
43
0 100 200 300 400 500 600
0
200
400
600
800
1000
1200
Test
Producer1
(a) Base
0 100 200 300 400 500 600
300
400
500
600
700
800
900
1000
1100
Test
Producer1
(b) HCNO
0 100 200 300 400 500 600
150
200
250
300
350
400
450
500
550
Test
Producer2
(c) Base
0 100 200 300 400 500 600
150
200
250
300
350
400
450
500
Test
Producer2
(d) HCNO
Figure 3.12: Observed and predicted production curves of producer1 and 2 with the least
test errors of test case no. 1 of Base algorithm for Field4
44
0 100 200 300 400 500 600
1000
1500
2000
2500
3000
3500
4000
Test
Producer3
(a) Base
0 100 200 300 400 500 600
1500
2000
2500
3000
3500
4000
Test
Producer3
(b) HCNO
0 100 200 300 400 500 600
50
100
150
200
250
300
350
400
450
500
Test
Producer4
(c) Base
0 100 200 300 400 500 600
150
200
250
300
350
400
450
500
Test
Producer4
(d) HCNO
Figure 3.13: Observed and predicted production curves of producer3 and 4 with the least
test errors of test case no. 1 of Base algorithm for Field4
45
0 100 200 300 400 500 600
100
200
300
400
500
600
700
800
900
Test
Producer5
(a) Base
0 100 200 300 400 500 600
300
400
500
600
700
800
900
Test
Producer5
(b) HCNO
0 100 200 300 400 500 600
600
800
1000
1200
1400
1600
1800
2000
2200
Test
Producer6
(c) Base
0 100 200 300 400 500 600
600
800
1000
1200
1400
1600
1800
2000
Test
Producer6
(d) HCNO
Figure 3.14: Observed and predicted production curves of producer5 and 6 with the least
test errors of test case no. 1 of Base algorithm for Field4
46
0 100 200 300 400 500 600
500
1000
1500
2000
2500
3000
Test
Producer7
(a) Base
0 100 200 300 400 500 600
1200
1400
1600
1800
2000
2200
2400
2600
2800
3000
Test
Producer7
(b) HCNO
0 100 200 300 400 500 600
50
100
150
200
250
300
350
400
Test
Producer8
(c) Base
0 100 200 300 400 500 600
100
150
200
250
300
350
400
Test
Producer8
(d) HCNO
Figure 3.15: Observed and predicted production curves of producer7 and 8 with the least
test errors of test case no. 1 of Base algorithm for Field4
47
Estimation of Permeability
The connectivity parameter π and time constant parameter τ indicate the locations of
high permeability. Figures 4.16 and 4.23 show the estimated connectivities (left) and
time constants (right) between each pair of injector and producer in Field1∼ Field4. The
(green) solid lines indicate the estimated connectivities and time constants. The magni-
tude of the connectivity is proportional to the thickness and length of the (green) solid
lines. For homogeneous fields with constant permeability, the magnitude of the connec-
tivity is inversely proportional to the distance between the injector and producer. If
boundary effects are ignored for the homogeneous fields (Field1 and Field3), the lengths
of the lines from an injector to producers at the same distance are expected to be the
same. For the heterogeneous fields (Field2 and Field4), the connectivity map is expected
to locate particularly longer and thicker (orange) striped bars between injectors and pro-
ducers which have direct communication paths. Well pairs with low connectivity caused
by low permeability, should be connected by shorter and thinner solid lines. Longer and
thicker solid lines are supposed to be located between wells for which the permeability is
particularly high. Injectors having such direct communication path with particular pro-
ducers have very weak influence on other producers. Overall, both algorithms generated
more accurate connectivity maps compared to the constant maps. For Field1, injector 3
(Inj3) is supposed to have the same strength of connectivity with all four producers, but
the map of the base algorithm does not reproduce the relationship. For Field2, HCNO
located the two strong connectivities at (Inj1, Prod1) and (Inj3, Prod4) in both of the
maps, but the base algorithm did not represent the relationship correctly.
48
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(a) Base
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(b) HCNO
Figure 3.16: True (grey) and estimated (black) connectivity map for Field1: Base vs
HCNO
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(a) Base
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(b) HCNO
Figure 3.17: True and estimated time constant map for Field1: Base vs HCNO
49
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(a) Base
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(b) HCNO
Figure 3.18: True and estimated connectivity map for Field2: Base vs HCNO
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(a) Base
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(b) HCNO
Figure 3.19: True and estimated time constant map for Field2: Base vs HCNO
50
0 200 400 600 800 1000 1200 1400 1600
0
200
400
600
800
1000
1200
1400
1600 Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(a) Base
0 200 400 600 800 1000 1200 1400 1600
0
200
400
600
800
1000
1200
1400
1600 Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(b) HCNO
Figure 3.20: True and estimated connectivity map for Field3: Base vs HCNO
0 200 400 600 800 1000 1200 1400 1600
0
200
400
600
800
1000
1200
1400
1600 Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(a) Base
0 200 400 600 800 1000 1200 1400 1600
0
200
400
600
800
1000
1200
1400
1600 Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(b) HCNO
Figure 3.21: True and estimated time constant map for Field3: Base vs HCNO
51
0 200 400 600 800 1000 1200 1400 1600
0
200
400
600
800
1000
1200
1400
1600 Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(a) Base
0 200 400 600 800 1000 1200 1400 1600
0
200
400
600
800
1000
1200
1400
1600 Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(b) HCNO
Figure 3.22: True and estimated connectivity map for Field4: Base vs HCNO
0 200 400 600 800 1000 1200 1400 1600
0
200
400
600
800
1000
1200
1400
1600 Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(a) Base
0 200 400 600 800 1000 1200 1400 1600
0
200
400
600
800
1000
1200
1400
1600 Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(b) HCNO
Figure 3.23: True and estimated time constant map for Field4: Base vs HCNO
52
Chapter 4
Inference in Probabilistic Graphical Model
4.1 Related Work
4.1.1 MarkovandBayesianNetworksStructureLearningandParameter
Estimation
Significant amount of research has been done in structure learning of Markov network
(MN) or Markov random field (MRF). MN structure was learned from data based on
the use of L
1
regularization on the weights of the log-linear model, which had the ef-
fect of biasing the model towards solutions where many of the parameters were zero
[34]. The formulation converted the MN learning problem into a convex optimization prob-
lem in a continuous space, which could be solved using efficient gradient methods. The
Dynamic Grow-Shrink Inference-based Markov network learning algorithm (DGSIMN)
[16]. DGSIMN works by conducting a series of statistical conditional independence tests
toward the goal of restricting the number of possible structures to one. Thus, inferred
that structure as the only possibly correct one. DGSIMN selected the locally optimal
test that was expected to increase the state of the knowledge about the structure the
53
most. Ensemble-of-trees (ET) model was used for the MN structure learning [37]. The
key feature of the ET model is that although it needs to sum over exponentially many
spanning trees, its partition function as well as data likelihood has a closed form. The
ET model is versatile to different MNs such as discrete or continuous, Gaussian or non-
Gaussian, low or high tree width, and so on. The ET model tends to fit data using
as small number of spanning trees as possible, thus it provides a natural regularization
mechanism for finding a sparse representation of a MN. In constrast to those methods
that employ L
1
-norm sparsity regularization, the ET model does not need to be set any
regularization. It was shown that the problem of learning a maximum likelihood MN
given certain observed data could be reduced to the problem of identifying a maximum
weight low-treewidth graph under a given input weight function [26]. They developed the
constant factor approximation algorithm for this problem. An automatic construction of
a MN was considered as an alternative to learning classification rules [55]. They presented
a procedure for constructing a MN from a sample set of observations, which minimized
the Kullback-Leibler cross-entropy between the network under construction and the set of
samples in a stepwise fashion until a preset threshold was reached. The Kullback-Leibler
cross-entropy allowed to choose a full joint probability distribution in the absence of com-
plete information about the exact distribution.
It has been shown that learning BNs from data is NP-hard [12], [51]. But approximate
methods have been proposed to learn BNs in an acceptable amount of time. Chicker-
ing [11] showed that among those BNs where each node had at most K parents, the
search problem of identifying a BN that had a relative posterior probability grater a
54
given constant was NP-complete when the Bayesian metric, called the BDe metric, was
used. Instead of a single task, a BN structure learning algorithm for a number of related
tasks, e.g., gene expression data for a number of related species, was proposed [44]. Heuris-
tic search was used to find a high scoring set of structures, where the score for a set of
structures was computed in a principled way. Given a data set D and a NB structure G,
P(G—D) was treated as a score so that one could search for a high scoring network using
the heuristic search in [25]. Bayesian networks structures were learned based on memory
binary particle swarm optimization method and the Minimum Description Length (MDL)
principle [36]. The purpose of the added memory feature was to prevent and overcome
premature convergence by providing particle specific alternate target points to be used
at times instead of the best current position of the particle. Due to the extremely large
search space in methods that cast learning BN as an optimization problem, the search
space was restricted [17]. The iterative algorithm, SparseCandidate algorithm, restricted
the parents of each variable to belong to a small subset of candidates. A network that
satisfied these constraints was searched, and the learned network was then used for se-
lecting better candidates for the next iteration. The structure of the BN was learned from
the complete data using K2 algorithm combined with genetic algorithm (GA) [58]. For
a method of sensor planning for a mobile robot localization problem, they presented the
causal relation between local sensing results, actions, and belief of the global localization
using a BN. A BN structure was searched based on Bayesian score, named the K2 algo-
rithm, to learn the causal relation between local environment information, robot action,
and global localization. The K2 algorithm used a constraint of ordering of nodes to re-
duce the search space and a GA was employed to search the best ordering. The recursive
55
autonomy identification (RAI) algorithm was proposed for constraint-based (CB) BN
structure learning. Conditional independence (CI) test was sequentially applied to learn
the structure. RAI algorithm learned edge direction and structure decomposition into
autonomous sub-structures that complied with the Markov property. The GA philosophy
was used for searching among alternative structure [31]. It was assumed that an ordering
between the nodes of the network structures. Next, the ordering assumption was released
by using a repair operator which converted illegal structure into legal ones. The structure
of a BN was learned by conducting a search through the space of equivalence classes of
BN using Ant Colony Optimization (ACO) [13].
4.1.2 Factor Graph Structure Learning and Parameter Estimation
Factor graphs subsume both MNs and BNs in the sense that every MN or BN can be
written as a factor graph of the same size [1]. Empirical entropy estimates were used
to select an approximate Markov blanket (MB) for each variable, and then the param-
eter estimation algorithm was used to estimate parameters and identify which factors
were likely to be irrelevant. It was shown that the class of factor graphs with bounded
factor size and bounded connectivity could be learned in polynomial time and polyno-
mial number of samples with the assumption that the data was generated by a network
in the class. Mirowski and LeCun [40] proposed a method for training dynamic factor
graphs (DFG) with continuous latent state variables. A DFG included factors modeling
joint probabilities between hidden and observed variables, and factors modeling dynam-
ical constraints on hidden variables. Given a training observation sequence, the optimal
state sequence was found by minimizing the energy using a gradient-based minimization
56
method. Second, the parameters of the model were updated using a gradient-based proce-
dure so as to decrease the energy. The two steps were repeated over all training sequences.
Steepest descent was used in estimating parameters in factor graphs [14]. By showing
that steepest descent was combined with the sum-product algorithm, they exposed that
steepest descent would be elegantly combined with other summary propagation algo-
rithms. Korl and Loeliger [28] estimated the autoregressive parameters using the sum-
product algorithm and factor graphs. A state-space form was derived for the AR model,
and then the messages were propagated over the corresponding factor graph. A paral-
lel inference algorithm on large factor graphs was proposed in the distributed memory
setting of computer clusters [22][23]. Graph-partitioning was incorporated to distribute
the work. The vertex with highest residual was updated because updating a vertex with
no residual wastes process cycles. A parallel implementation was designed for the loopy
belief propagation algorithm for factor graphs [39]. They found three main parameters
susceptible to be defined by the user; scheduling policies, stopping criteria, and initial
message values. The scheduling policies are a rule-based scheduling such that a node
sends messages to the neighbor nodes after receiving a fixed number of messages or re-
ceiving messages from all the nodes identified in SndSet. The stopping criteria are that
the algorithm stops after certain number of iterations or a maximum given number of
messages. The user can manually define the function values for each factor node. Elec-
tromyographic (EMG) signals were decomposed into its components using factor graph
57
and the sum-product algorithm [27]. The discrete-time model of EMG signals were rep-
resented by a factor graph and the sum-product algorithm operated by passing messages
forward and backward in the loopy factor graph.
4.2 Factor Graph and the Sum-Product Algorithm
When algorithms must deal with complicated global functions of many variables, they
often exploit the manner in which the given functions factor as a product of simpler local
functions so that computational efficiency can be derived. The key concept that we are
exploiting is the generalized distributive law (GDL) [2]. Let f (x,y,w)and g(x,z)begiven
real-valued functions, where x,y,z,and w are variables taking values in a finite set A
with q elements. Suppose we compute tables of the values of α(x, w)and β(y)thatare
defined as follows:
α(x,w)=
y,z∈A
f(x,y,w)g(x,z) (4.1)
α(x,w) is obtained by marginalizing out the variables y and z from the function f (x,y,z)·
g(x,z). The total number of arithmetic operations required for the computation of α(x,w)
is 2q
4
, i.e., q
4
summations and q
4
multiplications. Using factorization can simplify the
computation. Because of the distributive law, the sum in Eq. 4.1 factors
α(x,w)=(
y∈A
f(x,y,w))(
z∈A
g(x,z)). (4.2)
58
We compute tables of the functions
α
1
(x,w)=
y∈A
f(x,y,w), (4.3)
α
2
(x)=
z∈A
g(x,z), (4.4)
which requires a total of q
3
+ q
2
operations. Then we compute the q
2
values of α(x,w)
using Eq. 4.2
α(x,w)= α
1
(x,w)α
2
(x). (4.5)
which requires q
2
multiplications. Thus, the total number of operations required to com-
pute α(x,w) can be reduced from 2q
4
to q
3
+2q
2
. It clearly shows that computational
efficiency can be achieved by using the factorization.
A factor graph is a graphical model that visualizes the structure of the factorization
of a complicated global function and dependency among variables [30][3]. A factor graph
is a standard graphical language that represents a mathematical relation such as the is
an argument of relation between variables and local functions. Computational efficiency
can be derived by exploiting the factorization of the global function into a product of
simpler local functions. Let g(x
1
,x
2
,x
3
,x
4
,x
5
) be a function of five variables, and suppose
that g can be expressed as a product of five factors.
g(x
1
,x
2
,x
3
,x
4
,x
5
)= f
A
(x
1
)f
B
(x
2
)f
C
(x
1
,x
2
,x
3
)f
D
(x
3
,x
4
)f
E
(x
3
,x
5
). (4.6)
Figure 4.1 is the corresponding factor graph to Eq. 4.6. The sum-product algorithm
59
Figure 4.1: Factor Graph for Eq. 4.6
is a generic message-passing algorithm operating in factor graphs that computes either
exactly or approximately various marginal functions derived from the global function. For
example, the summary for x
2
by borrowing the notation from [30] is
g
2
(x
2
)=
∼x
2
g(x
1
,...,x
5
), (4.7)
where g
2
is the marginal function associated with g(x
1
,...,x
5
) and the summation is over
all the variables but x
2
. The summary for x
2
can be rewritten as the product of messages
from neighboring nodes and each message consists of messages from neighboring nodes
and each message consists of sum-of-products
g
2
(x
2
)= μ
f
B
→x
2
(x
2
)μ
f
C
→x
2
(x
2
), (4.8)
where μ
f
B
→x
2
(x
2
)=
∼x
2
f
B
(x
2
)= f
B
(x
2
)and μ
f
C
→x
2
(x
2
)=
x
2
(f
A
(x
1
)f
C
(x
1
,x
2
,x
3
)·
(
∼x
3
f
D
(x
3
,x
4
))(
∼x
3
f
E
(x
3
,x
5
))).
Several algorithms developed artificial intelligence are instances of the sum-product algo-
rithm, i.e., forward/backward algorithm, the Viterbi algorithm, the iterative turbo coding
60
algorithm, Pearl’s belief propagation algorithm for Bayesian networks, the Kalman filter,
and certain fast Fourier transform (FFT) algorithms. A message-passing schedule in a
factor graph is a specification of messages to be passed during each clock tick [29][30]. A
variety of message-passing schedules are possible. Flooding schedule calls for a message
to pass in each direction over each edge at each clock tick. In a serial schedule, at most
one message is passed anywhere in the graph at each clock tick.
4.3 Graphical Representation of the CRM
This section describes how the objective function and the constraints, Eq. 2.2∼ 2.5, are
represented by a factor graph. As a simple illustrative example, suppose that an oil field
consists of 2 injectors and 2 producers (2× 2 field) with the two injectors influencing the
two producers as shown in Figure 4.2, where a circle with an arrow represents an injector,
a triangle represents a producer, and dotted lines stand for the true connectivity. The
(a) 2×2field
Tau11
F11
Tau21
F21
Pi11 Pi21
Y11
F1
Y21
Tau12
F12
Tau22
F22
Pi12 Pi22
Y12
F2
Y22
Fe1 Fe2
(b) Factor graph
Figure 4.2: 2× 2 field and its factor graph
corresponding factor graph consists of three types of variable nodes and three types of
61
factor nodes. The continuous search spaces of π
ij
and τ
ij
are lower and upper bounded,
and they are discretized by certain intervals. Pi
ij
is thevariablenodefor themarginal
probability distribution of π
ij
, Tau
ij
is the variable node for the marginal probability
distribution of τ
ij
. Y
ij
is the marginal probability distribution for sets of possible pro-
duction rates for producer j, that are generated by the contribution of injector i.Each
factor node is assigned a joint probability table consisting of all the neighboring variable
nodes. Fe
i
is assigned the joint probability table for the equality constraint, Eq. (3), and
computes marginal probabilities of π
ij
by forcing the equality constraint. F
ij
is assigned
the joint probability table consisting of the three neighboring variable nodes and com-
putes marginal probabilities of the variables. F
j
is assigned the joint probability table for
Y
ij
and computes marginal probability of the variable. The observed production rates are
exploited as evidences in the way that the joint probability distribution for F
j
is com-
puted by Boltzman distribution that takes root mean squared error (RMSE) between the
observed and the estimated production rates as input. Variable nodes for injection and
production rates are not included in the factor graph because they are observable. Instead,
injection rates are implicitly included in F
ij
and production rates are implicitly included
in F
j
. Joint probability distributions for Fe
i
and F
ij
are uniform distribution. The size
of the joint probability table for F
j
increases as the number of injectors connected to
producer j increases and discretization intervals become smaller. To maintain the table
of a reasonable size, we keep only top k entries with relatively higher probabilities based
on the observation that the probabilities of most of the entries are very close to 0.
62
4.4 Structure Learning of Factor Graphs for Parameter
Estimation
The dynamic construction of factor graph and parameter estimation consist of two phases,
belief propagation over the factor graph using the sum-product algorithm and insertion of
nodes and edges into the factor graph. The two phases are incorporated in Algorithms 1
and 2. To start the message passing, the initial factor graph is created by connecting each
producer to its closest injector based on the locality principle even if these connections
are not true ones. π
ij
and τ
ij
of the truly non-existing edges eventually converge to 0 and
the upper bound, respectively. Thus, our method identifies the truly non-existing edges as
well, as will be shown in the experimental results. Once the messages converge after the
first round, injector and producer pairs that have missing edges need to be identified. A
naive approach is to try the exhaustive combinatorial search and to perform the message
passing repeatedly until a stopping criterion is satisfied. This method will not work for
a large number of injectors and producers. To avoid such excessive computations, we
employ the greedy forward edge search that consists of three steps:
• Step1: For every injector i and producers connected to the injector, determine
whether the most probable value for π
ij
identified by the message from Fe
i
to π
ij
is not the same as that identified by the message from F
ij
to π
ij
.
• Step2: For every producer j, determine whether there are producers whose estima-
tion errors are not less than the error threshold.
63
• Step3: Among injectors selected in Step1 and producers selected in Step2, determine
whether there are injector-producer pairs whose injection rates are positively and
significantly correlated with the difference between the observed and the estimated
production rates.
When there is any missing edge between injector i and producers, a discrepancy is ob-
served such that the most probable value for π
ij
determined by the message from Fe
i
to
π
ij
is different from that determined by the message from F
ij
to π
ij
,where j represents
the indices of all producers that were connected to injector i. Algorithm 1 is the pseu-
docode for identifying those injectors.
Algorithm 1. findFreeInjectors(PiDiscrete, msgFe
i
ToPi
ij
, msgF
ij
ToPi
ij
, L, J)
1. out = {}
2. for i =1to L
3. for each p connected to injector i
4. [maxFe
i
,idx1] = max(msgFe
i
ToPi
ij
(i,p,:))
5. [maxF
ij
,idx2] = max(msgF
ij
ToPi
ij
(i,p,:))
6. if PiDiscrete(idx1) > PiDiscrete(idx2)
7. out ← out
{i}
8. end if
9. end for
10. end for
11. return out
64
When the high estimation error is observed in a producer, it implies that the producer
has missed contributions from injectors that are supposed to be connected to the pro-
ducer. Thus, producers that have missed connection with injectors can be detected by
looking into the estimation errors of those producers. Now, we have found injectors and
producers that could be pairs. Still, we do not know which injectors should be connected
to which producers. The last observation is that producers with high estimation errors
have residuals between the observed and the estimated production rates. The idea is that
the residual indicates the lack of the contributions from injectors for the producer, thus
a positive correlation (≥ 0.3) between the set of injection rates of injector i and the
residual is used to determine the injector-producer pairs. Algorithm 2 shows the pseu-
docode for the construction of factor graph and parameter estimation for estimating IPRs.
Algorithm 2. PGL
1. Input :
2. injData,q,JPs
3. P: all producer ids
4. PiDiscrete,TauDiscrete: discrete values
5. Output :
6. PiGuess,TauGuess,G
7. Begin
8. freeInjectors←{},freeProducers←{}
9. numRound ← 0, done ← FALSE
10. [connInjIds,G]← getClosestInjectorInitG(P)
65
11. while∼done
12. if numRound > 0
13. changed ← FALSE
14. for each pp in freeProducers
15. if∼isempty(freeInjectors)
16. [maxCorr,injId]= max(corr(freeInjectors,injData, q(pp)- ˆ q(pp)))
17. if maxCorr ≥ TOL CORR
18. connInjIds(pp)← connInjIds(pp)
{injId}
19. changed ← TRUE
20. end if
21. initializeMessages(connInjIds, pp)
22. end if
23. end for
24. if∼changed
25. break
26. end if
27. end if
28. [mpPi
ij
, mpTauij, msgFe
i
ToPi
ij
, msgF
ij
ToPi
ij
]← SumProductFlooding(JPs,
G,TOL DIFF)
29. [PiGuess,TauGuess]← computeExpectation(PiDiscrete, TauDiscrete, mpPi
ij
,
mpTau
ij
)
30. freeInjectors ← findFreeInjectors(PiDiscrete, msgFe
i
ToPi
ij
, msgF
ij
ToPi
ij
, L, P)
31. freeProducers ← findFreeProducers(PiGuess, TauGuess, injData, q,TOL ERR)
66
32. update G by PiGuess and TauGuess
33. ˆ q ← estimateProduction(PiGuess, TauGuess, injData)
34. numRound ← numRound +1
35. if max(calcTrainingError(ˆ q,q)) < TOL ERR
36. done ← TRUE
37. end if
38. end while
At each round, the factor graph is updated by adding new π
ij
, τ
ij
, Y
ij
and F
ij
,and
updating Fe
i
and F
j
. The updated factor graph is used for the message passing at the
next round. Accordingly, the messages for the updated factor graph are initialized for
the next round. findFreeProducers() finds producers whose estimation errors are not less
than TOL ERR. The convergence criterion of the sum-product algorithm is that the
maximum difference between two consecutive marginal probabilities of all π
ij
and those
of all τ
ij
is not larger than the difference threshold (TOL DIFF). The exit condition is
any of the two criteria: (1) the maximum estimation error among producers is less than
the error threshold (TOL ERR); (2) there is no change between two consecutive sets of
selected injectors. The goodness-of-fit of the estimated factor graph and the parameters is
determined by the estimation error. Since the factor graph is a connected graph, flooding
message passing scheme is used.
67
4.5 Performance Evaluation
The performance of the proposed method, PGL, was compared with that of HCNO in [32]
for synthetic oil field data. The performance was measured by RMSE for the test data
(injection and production rates), running time, and accuracy of locating injector-producer
connectivity. Because the underlying structure, i.e., information about permeability, of
real fields is unknown, synthetic data was used for validating the algorithm. For the
purpose of quantitatively comparing PGL and HCNO methods, we generated oil field
data based on the CRM. In this way, we can eliminate uncertainty in deriving π
ij
and τ
ij
parameters from either real oil field or reservoir simulation data. 80% data was used for
training and 20% was used as the test data. The joint probability tables for the factor
nodes were generated using the training data.
4.5.1 Experimental Setting
We used four oil fields of different types. Table 4.1 shows field name, type, size (number
of injectors × number of producers), number of time steps, and number of parameters.
Figure 4.3 shows the map of the four fields. The (red) circles stand for injectors and
Table 4.1: Field name, field type, size, number of time steps, and number of parameters
Field Type Size TS Param
Field1 Homogeneous 5× 4 160 40
Field2 Heterogeneous 5× 4 160 40
Field3 Homogeneous 8× 8 160 128
Field4 Heterogeneous 8× 8 160 128
68
the (black) squares do for producers. The grey lines in Field2 and Field4 indicate the
higher permeability while the permeability in Field1 and Field3 is even over the field. In
the homogeneous fields (Field1 and Field3), the contribution (i.e., connectivity) of an
injector is evenly distributed for all the neighboring producers and the injector has almost
no influence on distant producers. The connectivity between the injector-producer pairs
with the higher permeability is stronger than that of any other pairs. π
lb
is 0 and π
ub
is
1, and the discretization interval for π
ij
is 0.1 and that for τ
ij
is set to 4. τ
lb
and τ
ub
are
set to 0.0001 and 20, respectively. For HCNO, randomly generated 10 initial guesses were
used. A set of initial guesses was used for Field1 and Field2, and the other set of initial
guesses was used for Field3 and Field4.
69
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(a) Field1
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(b) Field2
0 500 1000 1500 2000 2500 3000
0
500
1000
1500
2000
2500
3000
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(c) Field3
0 500 1000 1500 2000 2500 3000
0
500
1000
1500
2000
2500
3000
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(d) Field4
Figure 4.3: Well locations and true higher permeability (grey lines) of four fields
70
4.5.2 Results
The performance of the two algorithms are compared by test error and running time. The
running time of HCNO is the search time until convergence, and that of PGL is the
message passing and edge search time until an exit condition is satisfied. Table 4.2 shows
the running times of HCNO and PGL. HCNO(mean) is the average running time of the
10 running times and HCNO(min) is the smallest one among them. As the number of
Table 4.2: Running times (Sec) of HCNO and PGL
Field HCNO(mean) HCNO(min) PGL
Field1 3.26 2.48 10.16
Field2 2.60 2.24 2.68
Field3 802.32 645.72 261.07
Field4 867.79 1007.64 69.54
wells increased, the running time of PGL became less than that of HCNO for larger
number of parameters. For even larger fields, the scalability of HCNO and PGL can be
analyzed as follows. The complexity of HCNO is determined by the convergence rate,
the number of quadratic programming (QP) subproblem iterations, and the number of
constrained linear least squares iterations. Let n be the number of injectors and m be the
number of producers. In practical situations, m is on the order of n,andweset m = n in
the asymptotic analysis. The number of QP subproblem iterations per SQP iteration is
mainly determined by the number of inequality constraints because an active-set method
is used to solve the QP subproblem. Thus, the number of QP subproblem iterations per
SQP iteration is 2n
2
. The number of constrained linear least squares iterations is 2n
2
+
71
n. The number of SQP iterations is mainly determined by the convergence rate. quasi-
Newton methods typically converge Q-superlinearly, thus the number of SQP iterations
is determined by the convergence rate. The worst case complexity of HCNO is thus ln
2
,
where l is the number of SQP iterations determined by the convergence rate. Therefore,
the complexity of HCNO is O(ln
2
). On the other hand, PGL begins with disjoint O(n)
subgraphs, and for practical sparse graphs in which the number of edges is O(n), the
iterative learning of the factor graph structure should terminate after O(k) iterations
assuming that a new edge per injector is added per iteration, where k is the maximum
degree of any node. For real physical large reservoirs, we expect k to be a small constant
k
n. Consequently, the complexity of PGL can be made O(kn). We thus expect that
PGL is scalable to the large problem size of real oilfields. It should also be noted that
PGL based on a locality principle and message passing is readily amenable to parallel
computing. Table 4.3 shows the test errors of HCNO and PGL that correspond to the
result in Table 4.2. The test errors of HCNO and PGL are comparable and are orders
Table 4.3: Test error (RMSE)
Field Base(mean) HCNO(mean) HCNO(min) PGL
Field1 8308.90 24.64 16.70 18.76
Field2 9208.90 44.70 19.80 82.93
Field3 1402.56 44.01 19.29 43.92
Field4 552.57 147.64 84.08 33.86
of magnitude smaller than the most commonly used base SQP algorithm. The test error
of HCNO is not necessarily less than that of PGL. For Field4, PGL shows lower error
72
than the lowest error by HCNO. It was found that the quality of the HCNO solution was
highly dependent on the initial guess, while that of the PGL solution was affected by
the discretization interval, i.e., smaller discretization interval improved the quality of the
solution. Poorly estimated large search space resulted in poor solutions and longer search
time for HCNO, while the smaller discretization interval resulted in longer preprocessing
time for generating the joint probability tables for PGL. The results of π
ij
of Field4
are shown in Table 4.4 and 4.5. The estimated π
ij
by HCNO is the one with the least
test error among the 10 results. In Table 4.4, 0.00 indicates that the value was not 0
but very close to 0. Overall, the two methods estimated the parameters accurately. The
observation is that (Inj3, Prod5) and (Inj5, Prod7) are not truly connected but PGL
selected the injectors for the initial connection because they were the closest injectors to
those producers. Finally, PGL assigned 0 to π
35
and π
57
as shown in the table. We did not
show the result for the time constant parameters (τ
ij
) but τ
35
and τ
57
were determined to
be τ
ub
by PGL. Overall, the time constant parameters estimated by the two methods were
not as accurate as the connectivity parameters (Table 4.5). In contrast to the connectivity
parameters, PGL assigned invalid values to the time constant parameters for the two well
pairs (Inj3, Prod5) and (Inj5, Prod7). HCNO assigned invalid values to the time constant
parameters of most of the truly non-existing connections. Figure 4.4 ∼ 4.15 show the
observed and estimated production curves of all the producers in the four fields. The
(blue) dashed curves are for observed production rates and the (red) dotted curves are
for the estimated production rates.
73
0 50 100 150 200
0
500
1000
1500
2000
2500
3000
Test
Producer1
(a) HCNO
0 50 100 150 200
0
500
1000
1500
2000
2500
Test
Producer1
(b) PGL
0 50 100 150 200
0
500
1000
1500
2000
2500
3000
Test
Producer2
(c) HCNO
0 50 100 150 200
0
500
1000
1500
2000
2500
Test
Producer2
(d) PGL
Figure 4.4: Observed and estimated production curves of Producer1 and 2 in Field1
74
Table 4.4: True and estimated π
ij
for Field4
prod1 prod2 prod3 prod4 prod5 prod6 prod7 prod8
True π
ij
Inj1 0.5 - 0.5 - - - - -
Inj2 0.33 0.33 - 0.34 - - - -
Inj30.1 - - 0.9 - ---
Inj4 - 0.33 - 0.33 - 0.34 - -
Inj5 -- -- 0.2 - - 0.8
Inj6 -- --- 0.9 - 0.1
Inj7 -- --- - 0.4 0.6
Inj8 -- --- 0.5 - 0.5
HCNO
Inj1 0.51 0.00 0.46 0.00 0.02 0.00 0.00 0.00
Inj2 0.34 0.31 0.00 0.33 0 0.00 0.00 0.01
Inj3 0.08 0.00 0.01 0.91 0.00 0.00 0.00 0.00
Inj4 0.00 0.29 0.02 0.32 0.01 0.35 0.01 0.01
Inj5 0.00 0 0.00 0 0.15 0.00 0.00 0.85
Inj6 0.00 0.03 0.01 0.00 0.00 0.88 0.00 0.08
Inj7 0.00 0.03 0.00 0.01 0.03 0.00 0.39 0.54
Inj8 0 0.00 0 0 0.00 0.52 0.00 0.48
PGL
Inj1 0.5 - 0.5 - - - - -
Inj2 0.3 0.3 - 0.4 - - - -
Inj30.1 - - 0.9 0 ---
Inj4 - 0.39 - 0.31 - 0.3 - -
Inj5 -- -- 0.2 - 0 0.8
Inj6 -- --- 0.9 - 0.1
Inj7 -- --- - 0.4 0.6
Inj8 -- --- 0.5 - 0.5
75
Table 4.5: True and estimated τ
ij
for Field4
prod1 prod2 prod3 prod4 prod5 prod6 prod7 prod8
True π
ij
Inj1 10 - 10 - - ---
Inj2 0.16 0.16 - 0.16 - - - -
Inj3 0.19 - - 1 - - - -
Inj4 -16 - 16 -16 - -
Inj5 -- -- 18 -- 2
Inj6 -- --- 1 - 19
Inj7 -- --- - 17 9
Inj8 -- --- 10 - 10
HCNO
Inj1 11.66 20 9.88 1.31 11.65 14.29 0.18 12.38
Inj2 4.12 3.28 4.33 6.56 0.45 20 3.60 14.00
Inj3 6.70 18.03 13.89 3.39 11.55 20 0.18 20
Inj4 4.80 16.88 10.12 19.51 11.37 18.32 10.53 13.68
Inj5 13.02 20 14.52 9.85 17.45 5.70 11.51 3.80
Inj6 15.22 20 20 19.93 20 3.55 20 20
Inj7 20 19.65 2.47 10.50 15.29 13.24 16.57 13.63
Inj8 11.06 15.02 10.26 16.11 13.38 13.99 6.60 11.97
PGL
Inj1 10.99 - 10.00 - - - - -
Inj2 3.09 2.24 - 8.93 - - - -
Inj3 8.71 - - 2.95 10.00 - - -
Inj4 - 18.65 - 17.01 - 15.65 - -
Inj5 - - - - 18.25 - 10.00 3.82
Inj6 - - - - - 3.39 - 14.38
Inj7 - - - - - - 16.82 14.08
Inj8 - - - - - 14.54 - 12.71
76
0 50 100 150 200
0
500
1000
1500
2000
2500
3000
Test
Producer3
(a) HCNO
0 50 100 150 200
0
500
1000
1500
2000
2500
Test
Producer3
(b) PGL
0 50 100 150 200
0
500
1000
1500
2000
2500
3000
Test
Producer4
(c) HCNO
0 50 100 150 200
0
500
1000
1500
2000
2500
Test
Producer4
(d) PGL
Figure 4.5: Observed and estimated production curves of Producer3 and 4 in Field1
77
0 50 100 150 200
0
500
1000
1500
2000
2500
Test
Producer1
(a) HCNO
0 50 100 150 200
0
500
1000
1500
2000
2500
Test
Producer1
(b) PGL
0 50 100 150 200
200
400
600
800
1000
1200
1400
Test
Producer2
(c) HCNO
0 50 100 150 200
200
400
600
800
1000
1200
Test
Producer2
(d) PGL
Figure 4.6: Observed and estimated production curves of Field2
78
0 50 100 150 200
200
400
600
800
1000
1200
1400
1600
1800
2000
Test
Producer3
(a) HCNO
0 50 100 150 200
200
400
600
800
1000
1200
1400
1600
1800
2000
Test
Producer3
(b) PGL
0 50 100 150 200
0
500
1000
1500
2000
2500
3000
3500
Test
Producer4
(c) HCNO
0 50 100 150 200
0
500
1000
1500
2000
2500
3000
3500
Test
Producer4
(d) PGL
Figure 4.7: Observed and estimated production curves of Field2
79
0 50 100 150 200
0
200
400
600
800
1000
1200
1400
1600
1800
Test
Producer1
(a) HCNO
0 50 100 150 200
0
200
400
600
800
1000
1200
1400
1600
1800
Test
Producer1
(b) PGL
0 50 100 150 200
100
200
300
400
500
600
700
800
900
1000
Test
Producer2
(c) HCNO
0 50 100 150 200
100
200
300
400
500
600
700
800
900
1000
Test
Producer2
(d) PGL
Figure 4.8: Observed and estimated production curves of proudcer1 and 2 in Field3
80
0 50 100 150 200
0
500
1000
1500
2000
Test
Producer3
(a) HCNO
0 50 100 150 200
0
500
1000
1500
2000
Test
Producer3
(b) PGL
0 50 100 150 200
0
200
400
600
800
1000
1200
1400
1600
1800
Test
Producer4
(c) HCNO
0 50 100 150 200
0
200
400
600
800
1000
1200
1400
1600
1800
Test
Producer4
(d) PGL
Figure 4.9: Observed and estimated production curves of proudcer3 and 4 in Field3
81
0 50 100 150 200
0
200
400
600
800
1000
1200
1400
1600
1800
Test
Producer5
(a) HCNO
0 50 100 150 200
0
200
400
600
800
1000
1200
1400
1600
1800
Test
Producer5
(b) PGL
0 50 100 150 200
200
400
600
800
1000
1200
1400
1600
Test
Producer6
(c) HCNO
0 50 100 150 200
200
400
600
800
1000
1200
1400
1600
Test
Producer6
(d) PGL
Figure 4.10: Observed and estimated production curves of proudcer5 and 6 in Field3
82
0 50 100 150 200
0
200
400
600
800
1000
1200
Test
Producer7
(a) HCNO
0 50 100 150 200
0
200
400
600
800
1000
1200
Test
Producer7
(b) PGL
0 50 100 150 200
200
400
600
800
1000
1200
1400
1600
Test
Producer8
(c) HCNO
0 50 100 150 200
200
400
600
800
1000
1200
1400
1600
Test
Producer8
(d) PGL
Figure 4.11: Observed and estimated production curves of proudcer7 and 8 in Field3
83
0 50 100 150 200
0
200
400
600
800
1000
1200
1400
1600
Test
Producer1
(a) HCNO
0 50 100 150 200
0
200
400
600
800
1000
1200
1400
1600
Test
Producer1
(b) PGL
0 50 100 150 200
100
200
300
400
500
600
700
800
900
1000
Test
Producer2
(c) HCNO
0 50 100 150 200
0
200
400
600
800
1000
1200
Test
Producer2
(d) PGL
Figure 4.12: Observed and estimated production curves of proudcer1 and 2 in Field4
84
0 50 100 150 200
100
200
300
400
500
600
700
800
900
Test
Producer3
(a) HCNO
0 50 100 150 200
100
200
300
400
500
600
700
800
900
Test
Producer3
(b) PGL
0 50 100 150 200
0
500
1000
1500
2000
2500
Test
Producer4
(c) HCNO
0 50 100 150 200
0
500
1000
1500
2000
2500
Test
Producer4
(d) PGL
Figure 4.13: Observed and estimated production curves of proudcer3 and 4 in Field4
85
0 50 100 150 200
50
100
150
200
250
300
350
400
Test
Producer5
(a) HCNO
0 50 100 150 200
50
100
150
200
250
300
350
400
Test
Producer5
(b) PGL
0 50 100 150 200
0
500
1000
1500
2000
2500
3000
Test
Producer6
(c) HCNO
0 50 100 150 200
0
500
1000
1500
2000
2500
3000
Test
Producer6
(d) PGL
Figure 4.14: Observed and estimated production curves of proudcer5 and 6 in Field4
86
0 50 100 150 200
100
200
300
400
500
600
700
800
Test
Producer7
(a) HCNO
0 50 100 150 200
100
200
300
400
500
600
700
800
Test
Producer7
(b) PGL
0 50 100 150 200
0
500
1000
1500
2000
2500
3000
Test
Producer8
(c) HCNO
0 50 100 150 200
0
500
1000
1500
2000
2500
3000
Test
Producer8
(d) PGL
Figure 4.15: Observed and estimated production curves of proudcer7 and 8 in Field4
87
The estimation curves visually show that both of the two methods estimated the
parameters accurately. Figs. 4.16 ∼ 4.19 show the true higher permeability and connec-
tivity estimated by the two methods for the four fields. The black lines represent the
estimated π
ij
. The length and width of the black lines indicate the magnitude of
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(a) HCNO
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(b) PGL
Figure 4.16: True higher permeability and estimated connectivity in Field1
the connectivity parameters. As shown in the figures, the true permeability was also well
estimated by the two methods. The contributions of the locally connected injectors were
evenly distributed among the neighboring producers and the higher permeability across
wells was also well estimated. Figs. 4.20 ∼ 4.23 show the true higher permeability and
time constant estimated by the two methods for the four fields. The black lines represent
the estimated τ
ij
. The length and width of the black lines indicate the magnitude of the
time constant parameters. As the figures show, HCNO assigned values to non-existing
connections because it assumes a complete connectivity between a set of injectors and a
88
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(a) HCNO
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(b) PGL
Figure 4.17: True higher permeability and estimated connectivity in Field2
set of producers. Therefore, the HCNO took longer search time than PGL did for larger
number of parameters as shown in Table 4.2.
89
0 500 1000 1500 2000 2500 3000
0
500
1000
1500
2000
2500
3000
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(a) HCNO
0 500 1000 1500 2000 2500 3000
0
500
1000
1500
2000
2500
3000
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(b) PGL
Figure 4.18: True higher permeability and estimated connectivity in Field3
0 500 1000 1500 2000 2500 3000
0
500
1000
1500
2000
2500
3000
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(a) HCNO
0 500 1000 1500 2000 2500 3000
0
500
1000
1500
2000
2500
3000
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(b) PGL
Figure 4.19: True higher permeability and estimated connectivity in Field4
90
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(a) HCNO
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(b) PGL
Figure 4.20: True higher permeability and estimated time constant in Field1
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(a) HCNO
0 500 1000 1500 2000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
Inj1 Inj2
Inj3
Inj4 Inj5
Prod1
Prod2 Prod3
Prod4
(b) PGL
Figure 4.21: True higher permeability and estimated time constant in Field2
91
0 500 1000 1500 2000 2500 3000
0
500
1000
1500
2000
2500
3000
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(a) HCNO
0 500 1000 1500 2000 2500 3000
0
500
1000
1500
2000
2500
3000
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(b) PGL
Figure 4.22: True higher permeability and estimated time constant in Field3
0 500 1000 1500 2000 2500 3000
0
500
1000
1500
2000
2500
3000
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(a) HCNO
0 500 1000 1500 2000 2500 3000
0
500
1000
1500
2000
2500
3000
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
Inj1 Inj2
Inj3 Inj4
Inj5 Inj6
Inj7 Inj8
Prod1 Prod2
Prod3 Prod4
Prod5 Prod6
Prod7 Prod8
(b) PGL
Figure 4.23: True higher permeability and estimated time constant in Field4
92
Chapter 5
Conclusions
The parameter estimation for the constrained nonlinear system has been mapped to the
two different research problems, the constrained nonlinear optimization and inference in
probabilistic graphical model. We have designed a hybrid algorithm for constrained non-
linear optimization using coupled linear and nonlinear optimization. Rapid reduction of
the function gradients has been achieved by separating the nonlinear and linear optimiza-
tion subproblems. The experimental results show that the hybrid method improves the
optimization process by improving the convergence percentage, by reducing the search
time and by lowering the residual error. Providing solutions estimated by CLLS to the
nonlinear optimization process in the SQP and a line search framework has enabled the
constrained nonlinear optimization to converge faster by helping reduce the function gra-
dient. Finally, we have demonstrated the effectiveness of the HCNO algorithm by applying
it to an important reservoir modeling problem in petroleum engineering.
We have developed the graph structure learning and parameter estimation method using
factor graph and the sum-product algorithm. Instead of starting with a fully connected
93
graph for the given constrained nonlinear system, the graph structure was built incre-
mentally by learning from belief discrepancies, estimation error, and correlation analy-
sis. As the number of wells increased, the running time of PGL became less than that of
HCNO. The prediction accuracy of PGL was comparable to that of HCNO despite of the
error caused by the finite discrete values for the two sets of parameters. The three regu-
larizations prevented the greedy forward search from including truly non-existing edges in
the factor graph. We have demonstrated that the approach of the probabilistic inference
in the loopy graph structure was successful in estimating parameters of the constrained
nonlinear model.
In this study, we dealt with the large-scale parameter estimation problem for the con-
strained nonlinear system. The hybrid optimization method significantly reduced the
search time by dimensionality reduction. Not only the nonlinear system but also the con-
straints were represented by the graph. To the best of our knowledge, this is the first time
that constraints were represented in a graphical model and those constraints were forced
through the local message passing. When loops are present in a network, messages usually
oscillate. Such oscillation happens in non-probabilistic systems but do not normally occur
in probabilistic networks due to the stochastic nature of the link matrices. We designed a
probabilistic graphical model and demonstrated that messages converged to a stable equi-
librium in the loopy probabilistic system. It was also demonstrated that approximation
by inference in a probabilistic graphical model resulted in comparable prediction accuracy
to the solution by the optimization method. Learning in the probabilistic graphical model
for the IPRs opened the problem that structure learning of factor graphs and inference
94
in a probabilistic graph model could be used for the parameter estimation for large-scale
constrained nonlinear system, for which no analytical solutions exist and nonlinear pro-
gramming was the only method. The learned factor graph visualized the estimated IPRs
in a sophisticated way.
Future research could be parallel implementations of PGL. Since the flooding message
passing is inherently parallel processing, parallel computing techniques can be applied
to the message passing of PGL. As data size increases, it is impractical to solve them
on a single computer, especially given limited computer memory. A parallel implementa-
tion of PGL is expected to reduce the running time significantly than that of the hybrid
optimization method. Another problem for which PGL could be estimating parameters
for nonlinear and multiple regression (NMR). NMR can be solved by inference in prob-
abilistic graphical model. NMR is a simpler instance of the CRM since constraints are
not considered. The parameters to be estimated would be the linear coefficients given
independent and observed data and the lower and upper bounds for the exponents.
95
Bibliography
[1] P. Abbeel, D. Koller, A. Y. Ng. (2005). Learning Factor Graph in Polynomial Time
& Sample Complexity. Uncertainty in Artificial Intelligence.
[2] S. M. Aji and R. J. McEliece. (2000). The Generalized Distributive Law.IEEE Trans.
Information Theory. 46(2):325-343.
[3] C. Bishop. (2006). Pattern Recognition and Machine Learning. Springer.
[4] L. Boer, J. Poncet, P. Biver, P. Henriquel, V. Marlot. (2006). A New Technique
To Achieve History Match Using a Probabilistic Approach. SPE Annual Technical
Conference and Exhibition.
[5] P. T. Boggs and J. W. Tolle. (1995). Sequential Quadratic Programming. Acta Nu-
merica. 4:1-51.
[6] P. T. Boggs, A. J. Kearsley, and J. W. Tolle. (1999a). A Practical Algorithm for Gen-
eral Large Scale Nonlinear Optimization Problems. SIAM J. Optimization. 9(3):755-
778.
[7] P. T. Boggs, A. J. Kearsley, and J. W. Tolle. (1999b). A Global Convergence Anal-
ysis of an Algorithm for Large-Scale Nonlinear Optimization Problems. SIAM J.
Optimization 9(4):833-862.
[8] S. Boyd and L. Vandenberghe. (2004). Convex Optimization. Cambridge University
Press.
[9] C. G. Broyden. (1970). The Convergence of a Class of Double-rank Minimization
Algorithms. J. Institute of Mathematics and Its Applications 6:76-90.
[10] R. H. Byrd, F. E. Curtis, and J. Nocedal. (2008). An Inexact SQP Method for
Equality Constrained Optimization. SIAM J. Optimization. 19(1):351-369.
[11] D. M. Chickering. (1996). Learning Bayesian Networks is NP-complete. In D. Fisher
and H. J. Lenz, eds.:Learning from Data: Artificial Intelligence and Statistics V.
Springer Verlag.
[12] P. Dagum and M. Luby. (1993). Approximating Probabilistic Inference in Bayesian
Belief NetworksisNP-hard. Artificial Intelligence. 60(1):141-154.
[13] R. Daly and Q. Shen. (2009). Learning Bayesian Network Equivalence Classes with
Any Colony Optimization. J. Artificial Intelligence Research. 35(1):391-447.
96
[14] J. Dauwels, S. Korl, and H.-A. Loeliger. (2005). Steepest Descent on Factor Graphs.
IEEE ITSOC Information Theory Workshop on Coding and Complexity.
[15] J. L. Devore. (1995). Probability and Statistics for Engineering and the Sciences. 4th
ed. Duxbury Press.
[16] P. Gandhi, F. Bromberg, D. Margaritis. (2008). Learning Markov Network Structure
using Few Independence Tests. SIAM Conf. Data Mining.
[17] N. Friedman, I. Nachman, D. Pe´ er. (1999). Learning Bayesian Network Structure
from Massive Datasets: The ”Sparse Candidate” Algorithm. Conf. Uncertainty in
Artificial Intelligence.
[18] P. E. Gill, W. Murray, and M. H. Wright, Practical Optimization. Academic
Press,1981.
[19] P. E. Gill, W. Murray, and M. A. Saunders. (1984). Procedures for Optimization with
a Mixture of Bounds and General Linear Constraints. ACM Trans. Mathematical
Software. 10(3):282-298.
[20] P. E. Gill, W. Murray, M. A. Saunders and M. H. Wright. (1986). Considerations
of Numerical Analysis in a Sequential Quadratic Programming Method. Numerical
Analysis. 1230:46-62.
[21] P. E. Gill, W. Murray, and M. A. Saunders (2002). SNOPT: An SQP Algorithm for
Large-Scale Constrained Optimization. SIAM J. Optimization. 12(4):979-1006.
[22] J. E. Gonzalez, Y. Low, C. Guestrin, D. O’Hallaron. (2009a). Distributed Parallel
Inference on Large Factor Graphs. Uncertainty in Artificial Intelligence.
[23] J. Gonzalez, Y. Low, C. Guestrin. (2009b). Residual Splash for Optimally Paralleliz-
ing Belief Propagation. J. Machine Learning Research 5:177-184.
[24] J. Han and M. Kamber. (2006). Data Mining: Concepts and Techniques. 2nd ed.
Morgan Kaufmann.
[25] D Heckerman. (1999). A Bayesian Method for the Induction of Probabilistic Networks
from Data. Machine Learning.
[26] D. Karger and N. Srebro. (2001). Learning Markov Networks: Maximum Bounded
Tree-Width Graphs. ACM SIAM Symposium on Discrete Algorithms.
[27] V. M. Koch and H.-A. Loeliger. (2004). Decomposition of Electromyographic Signals
by Iterative Message Passing in a Graphical Model. Conf. IEEE Engineering in
Medicine and Biology Society.
[28] S. Korl, H.-A. Loeliger. (2004). AR Model Parameter Estimation: From Factor
Graphs to Algorithms. Int’l Conf. Acoustics, Speech, and Signal Processing.
97
[29] F. R. Kschischang and B. J. Frey. (1998). Iterative Decoding of Compound Codes
by Probability Propagation in Graphical Models. IEEE J. Selected Areas in Com-
munications. 16(1):219-230.
[30] F. R. Kschischang, B. J. Frey, and H.-A. Loeliger. (2001). Factor Graphs and the
Sum-Product Algorithm. IEEE Trans. Information Theory 47(2):498-519.
[31] P. Larra˜ naga, M. Poza, Y. Yurramendi, R. H. Murga, C. M. H. Kuijpers (1996).
Structure Learning of Bayesian Networks by Genetic Algorithms: A Performance
Analysis of Control Parameters. IEEE Trans. Pattern Analysis and Machine Intel-
ligence. 18(9):912-926.
[32] H. Lee, K. Yao, O. O. Okpani, A. Nakano, I. Ershaghi. (2010). Hybrid Constrained
Nonlinear Optimization to Infer Injector-Producer Relationships in Oil Fields. Int’l
J. Computational Science 4(1):1-22.
[33] H. Lee, K. Yao, A. Nakano. (2010b). Probabilistic Petrograph Parameter Estimation
via Targeted Restructuring of Graphs. Conf. Uncertainty in Artificial Intelligence.
Under review.
[34] S. Lee, V. Ganapathi, D Koller. (2006). Efficient Structure Learning of Markov
Networks using L
1
-Regularization. Conf. Neural Information Processing Systems.
[35] F. Leibfritz and E. W. Sachs. (1999). Inexact SQP Interior Point Methods and Large
Scale Optimal Control Problems. SIAM J. Optimization. 38(1):272-293.
[36] X. Li, S. Wang, X. He. (2006). Learning Bayesian Networks Structures Based on
Memory Binary Particle Swarm Optimization. T. D. Wang et al. (Eds.): SEAL,
LNCS 4247.
[37] Y. Lin and S. Zhu. (2009). Learning Sparse Markov Network Structure via Ensemble-
of-Trees Models. Conf. Artificial Intelligence and Statistics.
[38] F. Liu and J. M. Mendel. (2007). Forecasting Injector-Producer Relationships from
Production and Injection Rates Using an Extended Kalman Filter. SPE Annual
Technical Conference and Exhibition.
[39] A. Mendiburu, R. Santana, J. A. Lozano, E. Bengoetxea. (2007). A Parallel Frame-
work for Loopy Belief Propagation. Workshop on Parallel Bioinspired Algorithms in
conjunction with Genetic and Evolutionary Computation Conference.
[40] P. Mirowski, Y LeCun. (2009). Dynamic Factor Graphs for Time Series Modeling.
European Conf. Machine Learning and Principles and Practice of Knowledge Dis-
covery in Databases.
[41] K. P. Murphy, Y. Weiss, and M. I. Jordan. (1999). Loopy Belief Propagation for
Approximate Inference: An Empirical Study. Uncertainty in Artificial Intelligence.
[42] W. Murray. (1997).Sequential Quadratic Programming Methods for Large-Scale
Problems. Computational Optimization and Applications. 7(1):127-142.
98
[43] W. Murray and F. J. Prieto. (1995). A Sequential Quadratic Programming Algorithm
using an Incomplete Solution of the Subproblem. SIAM J. Optimization. 5(3):590-
640.
[44] A. Niculescu-Mizil and R. Caruana. (2007). Inductive Transfer for Bayesian Network
Structure Learning. Conf. Artificial Intelligence and Statistics .
[45] J. Nocedal and S. J. Wright. (1999). Numerical Optimization. Springer.
[46] J. Pearl. (1988). Probabilistic Reasoning in Intelligent Systems: Networks of Plausi-
ble Inference. Morgan Kaufmann.
[47] M. J. D. Powell. (1978). A Fast Algorithm for Nonlinearly Constrained Optimization
Calculations. Numerical Analysis. G.A.Watson ed. Lecture Notes in Mathematics.
Springer Verlag. 630:144-157.
[48] R. W. H. Sargent and M. Ding. (2000). A New SQP Algorithm for Large-Scale
Nonlinear Programming. SIAM J. Optimization, 11(3):716-747.
[49] M. Sayarpour, E. Zuluaga, C. S. Kabir, and L. W. Lake. (2007). The Use of
Capacitance-Resistive Models for Rapid Estimation of Waterflood Performance and
Optimization. SPE Annual Technical Conference and Exhibition.
[50] M. Sayarpour. (2008). Development and Application of Capacitance-Resistive Mod-
els too Water/CO2 Floods. PhD Thesis.
[51] S. E. Shimony. (1994). Finding Maps for Belief Networks is NP-hard. Artificial In-
telligence. 68(2):399-410.
[52] L. C. F. Silva, R. C. M. Portella, A. A. Emerick, N. F. F. Ebecken. (2007). Predic-
tive Data-Mining Technologies for Oil-Production Prediction in Petroleum Reservoir.
SPE Latin American and Caribbean Petroleum Engineering Conference.
[53] G. Strang. (1988). Linear Algebra and Its Applications. 3rd ed. Brooks Cole.
[54] R. W. Schulze-Riegert, J. K. Axmann, O. Haase, D. T. Rian, Y.-L. You. (2001).
Evolutionary Algorithms Applied to History Matching of Complex Reservoirs. SPE
Reservoir Simulation Symposium.
[55] S. K. M. Wong and Y. Xiang. (1994). Construction of a Markov Network from Data
for Probabilistic Inference. Workshop on Rough Sets and Soft Computing.
[56] A. A. Yousef, P. Gentil, J. L. Jensen, and L. W. Lake (2005). A Capacitance Model
to Infor Interwell Connectivity From Production and Injection Rate Fluctuations.
SPE Annual Technical Conference and Exhibition.
[57] A. A. Yousef. (2006). Investigating Statistical Techniques to Infer Interwell Connec-
tivity from Production and Injection Rate Fluctuations. PhD Thesis.
[58] H. Zhou and S. Sakane. (2003). Learning Bayesian Network Structure from Environ-
ment and Sensor Planning for Mobile Robot Localization. IEEE Conf. Multisensor
Fusion and Integration for Intelligent Systems.
99
Abstract (if available)
Abstract
In petroleum community, the oil field optimization, i.e., minimizing the operational cost and maximizing the oil recovery, is a challenging problem. One of the popular oil recovery techniques is waterflooding, which injects water into the oil reservoir to extract oil. Thus, the knowledge about injector-producer relationships (IPRs), i.e., which injectors contribute to the production of which producers, is a key for field optimization. The difficulty associated with field optimization is that the underlying structure of oil fields is unknown and it continuously changes over time.Recently, a capacitance-resistive model (CRM) has been proposed to investigate the IPRs. The CRM is a predictive model that predicts production rates given water injection rates. It consists of two sets of parameters, connectivity and time constant. The connectivity parameter quantifies the contribution of an injector to the production of a producer and the time constant parameter does the degree of fluid storage between the two wells. Three constraints are posed on the two sets of parameters: (1) the sum of the connectivity parameters of an injector to all other producers should be 1
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Inverse modeling and uncertainty quantification of nonlinear flow in porous media models
Asset Metadata
Creator
Lee, Hyokyeong
(author)
Core Title
Parameter estimation to infer injector-producer relationships in oil fields: from hybrid constrained nonlinear optimization to inference in probabilistic graphical model
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Computer Science
Publication Date
08/04/2010
Defense Date
05/06/2010
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
belief discrepancies,capacitance-resistive model,constrained linear least squares,constrained nonlinear optimization,dimensionality reduction,factor graph and the sum-product algorithm,hybrid optimization,injector-producer relationship,large-scale constrained nonlinear optimization,line search,locality principle,OAI-PMH Harvest,quasi-Newton method,sequential quadratic programming,structure learning of factor graphs
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Nakano, Aiichiro (
committee chair
), Haas, Stephan (
committee member
), Neches, Robert (
committee member
), Yao, Ke-Thia (
committee member
)
Creator Email
songof@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m3272
Unique identifier
UC1191449
Identifier
etd-Lee-3808 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-367302 (legacy record id),usctheses-m3272 (legacy record id)
Legacy Identifier
etd-Lee-3808.pdf
Dmrecord
367302
Document Type
Dissertation
Rights
Lee, Hyokyeong
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
belief discrepancies
capacitance-resistive model
constrained linear least squares
constrained nonlinear optimization
dimensionality reduction
factor graph and the sum-product algorithm
hybrid optimization
injector-producer relationship
large-scale constrained nonlinear optimization
line search
locality principle
quasi-Newton method
sequential quadratic programming
structure learning of factor graphs