Optimization of Oil Reservoir Models Using Tuned
Evolutionary Algorithms and Adaptive Differential
Evolution
Claus Aranha∗ , Ryoji Tanabe‡ , Romain Chassagne† and Alex Fukunaga‡∗ Graduate School of Systems and Information Engineering, University of Tsukuba, Ibaraki, Japan† Institute of Petroleum Engineering, Heriot-Watt University, Edinburgh, Scotland‡ Graduate School of Arts and Sciences, The University of Tokyo, Tokyo, Japan

Emails:{caranha@cs.tsukuba.ac.jp,Romain.Chassagne@pet.hw.ac.uk,fukunaga@idea.c.u-tokyo.ac.jp,rt.ryoji.tanabe@gmail.com}
Abstract—In the petroleum industry, accurate oil reservoir
models are crucial in the decision making process. One criticalstep in reservoir modeling is History Matching (HM), where theparameters of a reservoir model are adjusted in order to improve
its accuracy and enhance future prediction. Recent works applied
evolutionary algorithms (EAs) such as GA, DE and PSO for the
HM problem, but they have been limited to classical versionsof these algorithms. A significant obstacle to applying EAs to
HM is that each call to the fitness function requires an expensive
simulation, making it difficult to tune the control parameters
for EAs in order to obtain the best performance. We apply
and evaluate state-of-the-art, adaptive differential algorithms(SHADE and jDE), as well as non-adaptive evolutionary algorithms (standard DE, PSO) that have been tuned using standard
black-box benchmark functions as training instances. Both ofthese approaches result in significant improvements compared tostandard methods in the HM literature. We also apply fitness
distance correlation analysis to the search space explored by ouralgorithms in order to better understand the landscape of the
HM problem.

I.

I NTRODUCTION

The History Matching (HM) problem is of great importancein the Oil Industry. Given a model for an oil reservoir, history
matching of this model is the tuning of the parameters ofthe model such that the output values for the model (for oil,gas, water, etc) match those that have been measured in thepast history of the real reservoir. A good History Matching
result means that the model for the reservoir will be accurate,
and thus useful for informing and guiding future decision
making, e.g., when/how to develop an existing oil field. History
Matching is essentially a parameter optimization problem. The
goal is to find a set of input parameters which results in an
accurate forecast model. However, because it is very difficult
to predict the influence of any one parameter in the final
output, a black box optimization approach is used, where eachparameter set is evaluated in a model simulation. Dependingon the complexity of the model, this simulation can take hours
or even days for complex cases.
For these reasons, there is recent interest in the use of meta
heuristic methods for History Matching, such as Differential
Evolution (DE) [1], Particle Swarm Optimization (PSO) [2],Ant Colony Optimization (ACO) [3], etc (see section II).
However, most of these approaches have two limitations:(1) They use simple variants of the “classical” evolutionary

algorithms, and (2) they rarely report enough information toallow for the reproduction and improvement on their results.
To address these issues, in this paper we evaluate andcompare DE and PSO, two methods that have shown goodresults in History Matching in the Oil Engineering literature in
recent years, as well as jDE [4] and SHADE [5], two adaptive
differential evolution methods. In particular, SHADE is a stateof-the-art algorithm which has been shown to perform very
well in the recent CEC2013 black-box optimization competition [6]. We evaluate these algorithms on three field models:
two from real world oil fields, and one from a synthetic field.
It is well known that the performance of Evolutionary
Algorithms (EAs) is greatly influenced by control parametersettings, and obtaining good performance on a particular
problem can require a significant amount of control parametertuning. The standard approach to control parameter tuning is to
run the algorithm on the problem many times using different
sets of control parameters, and picking the best-performing
set. While this can be automated and made relatively efficient
using tools such as SMAC (Sequential Model-based Algorithm
Configurator) [7], the tuning process fundamentally remains
the same, except that auto-tuners seek to intelligently sample
promising combinations of control parameters.
However, control parameter tuning of EAs for the History Matching problem is difficult because each call to thefitness function (i.e., evaluating a parametrized version of thereservoir simulator) requires an expensive simulation, where
each simulation can take hours or even days. To overcome
this problem, we tune each candidate EA using the automatedtuning tool SMAC [7], where instead of using the actualHistory Matching problem as the training instance, we use
the set of standard CEC2014 black-box benchmark problem
instances as the training set and tune the algorithms given
the same number of fitness function evaluations that will be
used for the History Matching problem. Thus, we can quickly
generate a set of tuned control parameter heuristics that can beapplied to the History Matching problem. Our results (SectionIV) show that this tuning process obtains control parameter setsthat significantly outperform the standard control parametersfor History Matching suggested in the petroleum engineering
literature.
Our comparative evaluation of both the tuned and untuned
versions of all algorithms shows that adaptive differential

evolution, particularly SHADE, performs well on the History
Matching problem. We also use the individuals generated in
the comparison experiment to perform an analysis of the search
space structure of the HM problem using Fitness DistanceCorrelation (FDC) [8], [9]. This allows us to better understandthe characteristics of HM as a parameter optimization problem,comparing it to optimization benchmarks, and guiding thesearch for more efficient methods and strategies.
The rest of the paper is structured as follows. We firstdescribe the History Matching problem (Section II), including
a review of previous work and a description of our benchmark
oil field data. Section III briefly describes the adaptive DE
algorithms that are evaluated in this paper. Section IV presents
the experimental comparison of tuned and untuned variantsof standard algorithms (DE, PSO) and adaptive DE (jDE,SHADE). Section V analyzes the search space of the History
Matching problem using fitness distance correlation analysis.
Finally, Section VI discusses our results and concludes withdirections for future work.
II.

T HE H ISTORY M ATCHING P ROBLEM

The fate of huge investments in the oil and gas industry
rests on decision making. Reservoir decision making, suchas field development planning, can be improved by having amore accurate forecast of the future behavior of the reservoir.
However, because of the difficulty of directly observing a largepart of a field’s quantitative data, it is hard to calibrate the fieldmodel with realistic values regarding the spatial distributionsof its properties. Some of these parameters may be “known”
from sparse and noisy data obtained from the reservoir, but it
is extremely difficult to measure this data in a precise way forevery part of the field.
The usual approach is to make a preliminary estimation offield parameters, such as porosity, permeability and so on, andthen run the model in a reservoir simulator. The results of thispreliminary run seldom match with the data observed fromthe field, because of the uncertainty in the initial estimates.
Consequently, the parameters are iteratively adjusted until the
updated model matches the measured data to a sufficientdegree of precision. This iterative process is known as “History
Matching” (HM), and used to adjust reservoir parameters. The
value for these parameters is adjusted based on the output of
the simulator. This process of adjusting the model’s parameterbased on the dynamic observation of change in the output value
is an inverse problem, in comparison to the reservoir model
simulator, which solves the forward simulation problem.
A recent review on the History Matching has been madeby Oliver and Chen [10]. Usually the problem is approachedby some sort of optimization method, be it a gradient based
method [11], or a non-gradient, such as simulated annealing [12] or evolutionary algorithms [13]. The main challengesof History Matching are its non-linearity, the potentially verylarge number of variables, and the excessive computation time
(which could take days in the largest cases).
A. Problem Description
An instance of the History Matching problem consists of: amodel for the studied field, a set of input parameters, and a setof output observations specifying the behavior of the reservoir

in time. The exact set of output observations may vary for eachreservoir model, but they generally include oil, gas and waterproduction for each well in the reservoir.
The goal of History Matching is to find an input parameter
set that produces outputs from the simulator that agree with themeasured data for that particular reservoir. The relevance of an
input parameter set is given by the misfit value. The misfit isdefined as the difference between the measured data for thefield, and the observed data from the simulator. For a given
parameter set P , a set of outputs O from the simulator overtime, and a set of measurements M , the misfit is calculated as
Misfit(P ) =T |O |
1 XXq i(Oj − Mji )2
T i=1 j=1

(1)where T is the number of time steps in the simulation, and|O| is the number of outputs (oil, gas, water) for each time
step.
The problem includes two kinds of constraints: box constraints and physical constraints. In this work, we deal withthe box constraints by limiting all candidate solutions to a(normalized) feasible region [0,1], where 0 is the minimal accepted value for a parameter, and a value of 1 is the maximum
accepted value for that parameter. Physical constraints referto occasions when the simulation does not converge because
a particular combination of parameter values is physicallyinfeasible. To deal with these constraints, we abort a simulation
run that takes more than a maximum expected amount of time
to complete, assuming that the simulation will never converge
in this case. The parameter set associated with this run will
be given the maximum misfit value (worst utility). Note thatthe parameters and box constraints used for optimization are
a subset of the complete set of model parameters. This subset
selection is currently made by a petroleum engineer using
domain knowledge (increasing the number of parameters being
optimized is a direction for future work).
B. Metaheuristics for the History Matching Problem
Metaheuristics (including Evolutionary Algorithms) havebeen a popular approach for the History Matching problem
for the last decade (see [14] for a recent survey). Among
recent approaches, Monfared et al. applied a simple GA toa History Matching problem [15]. Xavier et al. [16] also used
GA, using real-valued crossover operators. Hajizadeh tested
a variety of metaheuristic algorithms in a series of papers,
including ACO [17] and DE [18], on the TS2N oil field. He
used standard versions of these algorithms, and found that both
performed well. Arash [19] also used both standard DE andPSO in a different History Matching problem, and found that
DE outperformed PSO. Reynolds [20] studied a BOA-PSO
hybrid for the history matching problem. His result showed that
PSO by itself outperformed BOA, and slightly underperformed
the hybrid on a small synthetic case.
Based on these three previous lines of work, we use PSO
and standard DE as baselines for comparison vs. adaptive
DE variants (jDE and SHADE) in our experiments below.
A different, multi-objective approach has been studied by
Verga [21] using SPEA2. Multi-objective approaches are ofinterest for History Matching because for certain decision

TABLE I.

D IFFERENCES IN EXPERIMENTAL PARAMETERS FOR EACH
MODEL

Parameter
Input Parameters
Output Parameters
Time per Evaluation

TABLE II.

Fig. 1.

Fig. 2.

problems a lower misfit on one set of outputs rather thananother may be favored, so a set of diverse solutions may
be desirable.
C. Oil Reservoir Models
In this paper we study 3 different reservoir models: theTeal South Reservoir (TS2N), the Northern Sea Field (NSF)
and the Watt Field (WATT). The first two models refer to realworld fields, while the third one is a synthetic case.
Teal South Reservoir (TS2N): The Teal South reservoir is areal reservoir located in the Gulf of Mexico. Teal South has asingle production well that penetrates 4500 ft into the sand. It
was initially over-pressured at 3096 psi. Monthly production
data of Oil, Gas and Water is listed from 1996 to 1999. There
are five geological layers in the model with uniform properties.
The TS2N is a real life example, albeit a simple one. A
visualization of the TS2N model is shown in figure 1.
Northern Sea Field (NSF): NSF is a real reservoir in the
North Sea, and is the most complicated among these threecases as it is a complex, real field. Each full field evaluationof the NSF takes 8-10 hours of computing time, so we use an
up-scaled (lower resolution) version of the field, following themethod described by Christie [22]. Further details about thisfield are currently confidential.
WATT Reservoir (WATT): The Watt Field is a semi-synthetic
study, which means based in part on real field data and in parton synthetic in order to describe a realistic field example. It
has been developed in Arnold [23]. Top structure and wireline
data is based on real field data however the fluid properties,
relative permeability, and capillary data are synthetic. The fielddevelopment plan is also synthetic resulting in an artificial
production response. The field has been developed througha set of 16 horizontal production wells located across the

WATT
12
8
1-4h

NSF
8
3
10-50min

PARAMETERS FOR THE TS2N MODEL

Parameter
Porosity
Water Oil Contract
Rock Compressibility
Permeability X (Layer 1)
Permeability X (Layer 2)
Permeability X (Layer 3)
Permeability X (Layer 4)
Permeability X (Layer 5)
Permeability Z (Layer 1)
Permeability Z (Layer 2)
Permeability Z (Layer 3)
Permeability Z (Layer 4)
Permeability Z (Layer 5)

TS2N field, the colors define its 5 layers

WATT Synthetic Field

TS2N
13
3
2-10s

Minimum
0.18
4880
7.7e-5
60
60
60
60
60
60
60
60
60
60

Maximum
0.377
4935
9.7e-5
460
460
460
460
460
460
460
460
460
460

central part of the reservoir and 5 horizontal and 2 verticalinjection wells around the edges. As compared to Teal South
this case is much more complex in terms of physics involved.
A visualization of the WATT model is shown in figure 2.
Table I summarizes the differences between each model,
from a simulation point of view: In this table, Input parametersis the dimensionality of the optimization problem that we areusing in this experiment. The input parameters used for the
TS2N and WATT field are described in tables II and III.
Output parameters is the number of observations that themodel outputs at each time step of the simulation, and is used
as |O| in equation 1.
Time per evaluation is how long one simulation (oneevaluation) of the model would take. This running time ismeasured on a single core of a computer with an AMD Opteron
2384 CPU and 32GB of RAM. During the experiments, weparallelize the execution of the evolutionary algorithms byhaving each individual run on a different core of the machine.
III.

A DAPTIVE E VOLUTIONARY A LGORITHMS

In this paper we evaluate jDE [4] and SHADE [5], twostate-of-art variants of Differential Evolution (DE) [1]. First,
let us briefly review standard DE. A DE population is represented as a set of real parameter vectors xi = (x1 , ..., xD ),i = 1, ..., N , where D is the dimensionality of the target
problem, and N is the population size. In each generation t,
a mutant vector v i,t is generated from an existing population
member xi,t by applying some mutation strategy. One popular
TABLE III.

PARAMETERS FOR THE WATT MODEL

Parameter
Porosity Layers 1-4
Porosity Layers 5-8
Porosity Layers 9-12
Porosity Layers 13-16
Porosity Layers 17-20
Porosity Layers 21-24
Porosity Layers 25-32
Porosity Layers 33-40
Permeability Multiplier (4 parameters)

Minimum
0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.01

Maximum
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
1

mutation strategy is rand/1, defined as follow:
v i,t = xr1 ,t + F · (xr2 ,t − xr3 ,t )(2)

The indices r1 , r2 , r3 are randomly selected from [1, N ] suchthat they differ from each other as well as i. The parameterF ∈ [0, 1] controls the magnitude of the differential mutation
operator. After generating the mutant vector v i,t , it is crossedwith the parent xi,t in order to generate trial vector ui,t . Binomial Crossover, the most commonly used crossover operator
in DE, is implemented as follows:
v j,i,G if rand[0, 1) ≤ CR or j = jranduj,i,G =(3)xj,i,G otherwise
Where rand[0, 1) denotes a uniformly selected random number
and jrand is a decision variable index uniformly selectedfrom [1, D]. CR ∈ [0, 1] is the crossover rate. After allof the trial vectors ui,t , 0 ≤ i ≤ N have been generated,each individual xi,t is compared with its corresponding trial
vector ui,t . Then the better vector is kept in the population.
Below, we say that a generation of trial vector is successful if
this replacement occurs. This process is repeated until sometermination criterion is met.
The performance of DE depends significantly on thecontrol parameters such as F and CR (defined above). In
recent years, there have been much research on adaptive DE
variants that automatically adapt the control parameter values
as the search progresses. In this paper, we focus on jDE andSHADE because jDE is among the most cited adaptive DE
algorithms, and SHADE is the highest ranking DE methodin the CEC2013 competition on on Real-Parameter Single
Objective Optimization [6].
jDE [4] assigns a different set of parameter values Fi
and CRi to each xi , which is used for generating the trial
vectors. Initially, the parameters for all individuals xi are set
to Fi = 0.5, CRi = 0.9. The control parameter values of thetrial vectors are inherited from their parents. However, eachparameter is randomly modified (within a pre-specified range)
with some probability, and modified parameters are kept for
the next generation only when a trial is successful.
SHADE [5] is an adaptive DE which based on JADE [24]algorithm, but uses a historical memory of successful parameter settings based parameter adaptation scheme. Successhistory based adaptation uses a historical memory MCR , MF
which stores a set of CR, F values that have performed well in
the past, and generate new CR, F pairs by directly sampling
the parameter space close to one of these stored pairs.
IV.

C OMPARISON OF A DAPTIVE DE WITH P REVIOUS
M ETHODS

In this section, we evaluate the performance of adaptiveDEs (SHADE and jDE) on three instances of the History
Matching problem: TS2N, WATT and NSF, described in Section II. We compare these adaptive DEs with standard DE
and PSO algorithms that have been previously been evaluatedin the petroleum engineering literature (see Section II-B 1 ).
Our choice of jDE and SHADE as our adaptive differential
evolution algorithms is explained in Section III.
1 We used the standard gbest-model PSO using inertia weight strategy [2].

These four algorithms are used to generate input parameter
configurations for the three fields. Because of the computational cost of one evaluation (up to 50 minutes in the case of
NSF, and up to 5 hours in the case of WATT), the evaluationbudget for all algorithm runs is set to 500 evaluations. As thisis an unusually low number, we perform parameter tuning ineach of the four algorithms to compensate for it. Additionally,
we evaluate a uniform random sampling of 1, 000 parameter
sets from each of the three fields as a control sample.
For DE, PSO, jDE, and SHADE, we evaluate two configurations:Standard configuration: A parameter configuration taken fromrecent literature. In the case of DE and PSO, the implementation follow suggestions for configuring these algorithms fromrecent work in the History Matching literature [17], [20]. jDE
and SHADE follow the implementation described in Brest [4]and Tanabe [5], respectively.
Tuned configuration: A parameter configuration tuned for thesmaller budget of 500 evaluations used in this experiment. The
tuned parameter values for each algorithm were calculatedusing the Sequential Model-based Algorithm Configuration
method (SMAC) [7].
For the exact parameter values in the standard and tuned
configuration, see Table IV. In the tables and figures, the tuned
variants of the algorithms are denoted by the “tuned” prefix(tunedDE, tunedPSO, tunedjDE and tunedSHADE).
A. Parameter Tuning
As explained in Section I, the standard approach to parameter tuning by repeatedly running each algorithm with various
control parameter settings on the History Matching problem isimpractical because each call to the History Matching fitness
function involves a very expensive call to the reservoir simulator. We therefore perform parameter tuning using a standard
set of black-box benchmark functions as the tuning/training
instances, based on the belief that tuning the algorithms toperform well across a broad set of black-box benchmark
problems should result in good performance on the actual
problem of interest (the HM problem).
All 4 algorithms were tuned using the Sequential Modelbased Algorithm Configuration, SMAC [7]. An algorithm
configurator takes as input an algorithm executable, a formal
description of the parameters for the algorithm, and a set oftraining problem instances. It searches the space of possible
parameter values by repeatedly generating a candidate configuration (e.g., by local search) and evaluating the configurationon the set of the training instances (or some intelligently
selected subset of training instances). The configuration withhighest expected utility on the training set is returned. SMAC
can be used to tune real-valued, integer-valued, categorical,
and conditional parameters [7].2
The evaluation function used by SMAC to assess thequality of a candidate EA configuration was the mean of thedifference between the solution found by the EA configurationand the optimal value for each benchmark function in the training set, consisting of 30 functions F1 ∼ F30 in 10 dimensions
2 We used the most recent version of SMAC downloaded from the authors’
website, http://www.cs.ubc.ca/labs/beta/Projects/SMAC/.

TABLE IV.

F OR PSO, DE, J DE, SHADE, THE DEFAULT CONTROLPARAMETER VALUES , AS WELL AS THE BEST PARAMETERS FOUND BY
TUNING THE ALGORITHM WITH SMAC USING 10- DIMENSIONAL
CEC2014 BENCHMARK PROBLEMS F1 ∼ F30 AS TRAINING PROBLEMS .
N / A MEANS THAT ITS CONTROL PARAMETER IS NOT INCLUDED IN THE
TUNED PARAMETER SETTINGS .

PSO

DE

jDE

SHADE

Parameters
N

Range[2, 50]

Default
25

Tuned
10

w

[0, 1]

0.729

0.462

c1

[0, 4]

1.494

0.627

c2

[0, 4]

1.494

2.516

N

[6, 50]25

6F[0.1, 1]

0.5

0.730

CR

[0, 1]

0.5

0.401

strategy

see the text

rand/1current-to-pbest/1

p[0, 0.2]

0.05

0.175

archive rate

[0, 2]

1.0

1.738

N

[6, 50]25

7F[0.1, 1]

0.5

0.859

CR

[0, 1]

0.5

0.130

strategy

see the text

rand/1

best/1p[0, 0.2]

0.05

n/a

archive rate

[0, 2]

1.0

n/a

τF

[0.05, 0.3]

0.1

0.063

τCR

[0.05, 0.3]

0.1

0.098

N

[6, 50]25

8F[0.1, 1]

0.5

0.856

CR

[0, 1]

0.5

0.224

strategy

see the text

current-to-pbest/1“

p[0, 0.2]

0.05

0.017

archive rate

[0, 2]

1.0

0.996

memory size

[1, 20]

10

9from the CEC2014 benchmarks [25]. The maximum numberof fitness evaluations for EA was set to 500, following theevaluation budget of the HM optimization problem. Each runof SMAC was limited to 3, 000 EA configurations. For each
EA, for each training scenario, SMAC was run 5 times, andwe selected the best result out of these 5 runs.
For each EA, the default control parameter values (from[2], [24], [5]), the ranges for the parameters, as well as thevalues found by SMAC, are shown in Table IV. Binomial
crossover was used for all DE variants. For standard DE,
7 possible mutation strategies (rand/1, rand/2, best/1, best/2,
current-to-best/1, current-to-best/2, and current-to-pbest/1 witharchive) could be selected.3 The “current-to-pbest/1 witharchive” strategy has control parameter p and archive rate;
these are modeled as conditional parameters [7] in SMAC.
The archive size |A| is set to round(archive rate ×N ).
B. Results
The results of executing the four algorithms (and theirtuned variants) in the three History Matching fields are summarized in Table V and Figure 3. Table V shows the misfit
value obtained by each method for each field. In the caseof the TS2N and NSF fields, these values are the mean of
20 and 10 runs, respectively. The standard deviation of this
3 For details on these mutation strategies, see the survey by [26].

mean is given between parenthesis. The last column lists thebest misfit obtained from 1000 parameter sets drawn from a
uniform distribution, as a control. Figures 3a and 3b show box
plots of the misfit obtained for each repetition in the TS2N
and NSF, respectively.
In the TS2N reservoir, the tuned versions of the algorithms
easily outperformed the standard versions. TunedSHADE performed best, slightly ahead of tunedDE (Wilcoxon test p =0.10). Note that tunedDE included not only different values
for CR and F , but also uses the current-to-pbest/1 mutation,
which uses a “memory” of successful individuals.
For the NSF field, we observe somewhat similar results,
with the tuned versions of DE, jDE and SHADE outperforming
standard DE and PSO. In this field, the difference between
these three tuned algorithms is unclear (Wilcoxon test p = 0.48
for the null hypothesis that tunedDE equals tunedSHADE). For
this field, tuned PSO exhibits extremely anomalous behavior,
which reason we have not been able to identify.
For the WATT field, because of its extremely long evaluation time, we could only perform a single repetition of the
optimization. We observe that this result generally agrees tothe other two fields: The optimization methods outperform thereference random search, and tunedDE,tunedjDE,SHADE andtunedSHADE present the lowest misfit.
Overall, on these History Matching instances we observe
the following:•

All algorithms performed better than the trivial control
algorithm (random sampling).

•

The tuned versions of the algorithms generally performed better than the untuned versions, indicatingthat standard benchmark problems are useful as training problems for tuning evolutionary algorithms for
the History Matching problem.

•

Among untuned algorithms, SHADE, an adaptive DE
consistently performs the best.

•

All of the tuned algorithms achieve roughly similar average performance on these 3 instances. However, thestability of the tuned algorithms differ significantly.
Among tuned algorithms, tuned SHADE appears tohave the most stable performance.

V. A NALYSIS OF THE H ISTORY M ATCHING P ROBLEM
S EARCH S PACE U SING F ITNESS D ISTANCE C ORRELATIONIn this section, we experimentally analyze the structure of
the search space of the History Matching problem using Fitness
Distance Correlation (FDC) [8], [9]. FDC is a landscapeanalysis method which analyzes the correlation between f (x),the fitness of a candidate vector x and d(x, x∗ ) the distance
between x and the optimal solution vector x∗ . Sampling andplotting {f (x), d(x, x∗ )} for many points in the search space(e.g., Figure 4) can result in an informative, visualization of
the structure of the search space.
EA are often applied to black-box optimization problems
such as the History Matching problem where the structure of
the search space is unknown a priori. However, if it is possible
to know/predict the structure of the search space a priori, it isTABLE V.
Field
TS2N
NSF
WATT

DE
104.84(4.66)
1417.84 (318)
16938.96

M EAN M ISFIT VALUES FOR EACH METHOD AND EACH FIELD ( STANDARD DEVIATION IN PARENTHESIS )
Tuned DE
95.09 (4.01)
1301.33 (509)
16471.63

PSO
98.49 (2.82)
1660.24 (812)16510.02

Tuned PSO
96.70 (4.69)
2445.19 (1087)
16472.89

jDE
97.42 (2.59)
1113.84 (70)16564.55

Tuned jDE
96.17 (5.18)
1049.35 (85)16469.51

Shade
95.66 (2.2)
1067.24 (78)16469.51

Tuned Shade
93.83 (3.46)
1059.97 (86)16469.06

Random best
108.622916.0217490.13(a) TS2N Model(b) NSF Model
Fig. 3. Box plot of the best misfit values for all repetitions of the TS2N and NSF models. This means 20 repetitions for the TS2N case, and 10 for the NSF
case. The vertical axis in each figure is the misfit.

possible to adapt the EA to that structure in order to searchefficiently. For example, the population could be set to be smallor large depending on whether the search space is unimodal
or multimodal, respectively. Instances of a class of problems
tend to have similar search space structures. Thus, a FDC
analysis can guide the development of search algorithms for
the History Matching problem. In addition, while there havebeen a number of landscape analysis on synthetic benchmark
functions (c.f. [27]), landscape analysis of real world problems
are rare. Thus, we believe that a landscape analysis of theHistory Matching problem is a valuable contribution.
We measured the FDC of the 3 reservoir model instances(TS2N, WATT and NSF) described in Section II. Figures 4(a)–(c) plot {f (x), d(x, x∗ )} for all candidate vectors X that were
evaluated (in all trials, by all algorithms) for each instance inthe experiments in Section IV. Since the true optimal solutionx∗ is unknown, xb , the candidate vector with the lowest (best)
fitness value that was found among all of the sampled points, isused as an approximation/proxy for x∗ , i.e., the fitness distanceanalysis is performed relative to a point that is a possible, local
optimum, instead of the true global optimum.
For comparison and reference, Figures 4(d)–(g) show the
FDC for the 10-dimensional Sphere, Ellipsoidal, Rastrigin, andRosenbrock functions, where the points in the FDC plots wereobtained by uniformly sampling the search space 105 times.
The Sphere and Ellipsoidal functions are unimodal, while
the Rastrigin and Rosenbrock functions are multimodal. Also,
the Ellipsoidal function is ill-scaled, while the Rosenbrockfunction is non separable. Detailed definitions and analysisof the properties of the standard functions can be found in[28]. Table VI shows the rF DC metric [9], computed using
the following formula:rF DC =

cF D
sF sD

(4)

(a) TS2N(b) WATT(d) Sphere (D = 10)(e) Ellipsoidal (D = 10)(c) NSF(f) Rastrigin (D = 10)(g) Rosenbrock (D = 10)

Fig. 4. Figure (a)–(c) plots the Fitness Distance Correlation (FDC) for the three History Matching instances (TS2N, WATT and NSF). All candidate vectors
x evaluated by all algorithms, all trials for the experiments in Section IV are represented by a point in the plots. The horizontal axis represents the Euclidean
distance between x and xb , the best vector found for TS2N in Section IV (among all trials of all algorithms), and the vertical axis represents the fitness valueof x. Figure (d)–(g) show the FDC for the 10-dimensional Sphere, Ellipsoidal, Rastrigin, and Rosenbrock functions, where the points in the FDC plots wereobtained by uniformly sampling the search space 105 times.
ncF D =1X¯(fi − f¯)(di − d)n i=1

TABLE VI.

where n = |X|, fi , di are the fitness and distance (from xb¯ sF , sD are theor x∗ ) of the i-th vector (1 ≤ i ≤ n), and f¯, d,means and standard deviations of fitness values and distances,respectively. As rF DC approaches 1, the search space has amore pronounced big valley structure [8], which is well suitedfor Evolutionary Algorithms. On the other hand, rF DC close
to 0 indicates that the global structure is a weak “needle in ahaystack”, and a negative value of rF DC indicates a deceptive
structure (e.g., multi-funnel landscape).
Figure 4 shows that the FDC plot of the History Matching
instances (TS2N, WATT, NSF) resembles that of the standard
Ellipsoidal benchmark function. In addition, Table VI showsthat for both the History Matching instances and the Ellipsoidal
function, rF DC is low. Thus, it can be inferred that TS2N has
a search space structure that is somewhat similar to that of the
Ellipsoidal function. This result may have a future application:
In general, it is desirable for the training used for parametertuning to be similar to the test problems for which the tuning
is ultimately intended. In Section IV-A, our parameter tuning
used all 30 instances of the CEC2014 benchmarks [25] as thetraining set. By using a different set of training instances which
focus on problems similar to the Ellipsoid function, we may
be able to find control parameter values for EAs that are better
than the parameters we found in Section IV-A.
VI.

T HE rF DC METRICS [9] FOR EACH FUNCTION .

(5)

D ISCUSSION AND C ONCLUSION

We applied two adaptive differential evolution algorithms,
jDE and SHADE, to three different Oil Reservoir History

TS2N

WATT

NSF

SphereEllipsoidal

Rastrigin

Rosenbrock

0.33

0.49

0.22

0.99

0.40

0.61

0.82

matching problem instances. We experimentally compared
tuned and standard configurations of these adaptive algorithms,
as well as non-adaptive standard metaheuristics that have previously used for History Matching (PSO and DE). Our results
show that the tuned algorithms significantly outperformed thestandard versions from the petroleum engineering literature.
Among the untuned algorithms, SHADE consistently performed best. Although all of the tuned algorithms performed
comparably, the tuned version of SHADE had the most stableperformance.
Standard control parameter tuning of Evolutionary Algorithms is impractical on the History Matching problem due tothe extremely expensive fitness function. However, we showedthat tuning the EAs using standard black-box benchmark
problems as a surrogate (which is orders of magnitude fasterthan the oil reservoir simulation) is successful. All of thealgorithms tuned using this method performed significantly
better than their default configurations. In particular, it is worth
noting that standard DE and PSO using parameters tuned
with this methodology significantly outperformed DE and PSO
using control parameter values suggested by previous works on
the History Matching problem.
We also analyzed the search space of the History Matching
instances using the candidate solution vectors generated during

our optimization runs. The TS2N, WATT and NSF models
were all shown to have similar search spaces that are visibly
similar to the standard Ellipsoidal benchmark function, andthey also have rF DC values which are similar to the Ellipsoidal
function, suggesting that future work should focus on training
optimization algorithms for the HM problem using functions
similar to Ellipsoidal function as a surrogate model.
This paper focused on the application of relatively simple, evolutionary algorithms to the HM problem. Applicationof more complex methods that use surrogate models (e.g.,
HCMA [29]) is a direction for future work. Another interesting
challenge is the addition of seismic constraints to both themodel and input parameters. This means using a combinationof measured data from the wells (as in this study) with a 3D
image from the ground. This might generate an even more
precise model, in exchange of a more difficult evaluation
function. In the current work, sometimes the simulation did
not converge due to an unphysical parameter combination.
These failed simulations cannot be evaluated, but take precious
computational time. It would be interesting to develop a way
to learn these unphysical combinations as problem constraints.

[8]

[9]
[10]

[11]

[12]
[13]

[14]

[15]

[16]

ACKNOWLEDGMENT
R. Chassagne is supported by the sponsors of the Edinburgh
Time-Lapse Project, Phase V (BG, BP, CGG, Chevron, ConocoPhillips, ENI, ExxonMobil, Hess, Ikon Science, Landmark,Maersk, Nexen, Norsar, Petoro, Petrobras, RSI, Shell, Statoil,
Suncor and Total). A. Fukunaga is supported by a KAKENHI
grant. Thanks to Schlumberger for use of the ECLIPSE™ industry standard reservoir simulator software.

[17]

[18]
[19]

[20]

A PPENDIX A: C ODE AND DATA AVAILABILITY
All the source code, parameter files and the data for the
TS2N and WATT fields, used in this paper, are available at
http://github.com/caranha/HistoryMatchingProject. The
data for the NSF reservoir is not publicly available. A reservoir
simulator is required to perform the fitness evaluation. The
experiments in this paper used the ECLIPSE™ reservoir
simulator. The use of other simulation software is possible,
but may require some adjustment of the available code.

[22]

R EFERENCES

[24]

[1]

[2]
[3]

[4]

[5]
[6]
[7]

R. Storn and K. Price, “Differential Evolution - A Simple and Efficient
Heuristic for Global Optimization over Continuous Spaces,” J. Global
Optimiz., vol. 11, no. 4, pp. 341–359, 1997.
Y. Shi and R. E. and, “A modified particle swarm optimizer,” in IEEE
CEC, 1998, pp. 69–73.
M. Dorigo, V. Maniezzo, and A. Colorni, “Ant system: optimization bya colony of cooperating agents,” IEEE Trans. on Systems, Man, andCybernetics, Part B, vol. 26, no. 1, pp. 29–41, 1996.
J. Brest, S. Greiner, B. Bošković, M. Mernik, and V. Žumer, “SelfAdapting Control Parameters in Differential Evolution: A Comparative
Study on Numerical Benchmark Problems,” IEEE Tran. Evol. Comput.,vol. 10, no. 6, pp. 646–657, 2006.
R. Tanabe and A. Fukunaga, “Success-History Based Parameter Adaptation for Differential Evolution,” in IEEE CEC, 2013, pp. 71–78.
CEC2013, “http://www.ntu.edu.sg/home/EPNSugan/index files/CEC2013/CEC2013.htm.”
F. Hutter, H. H. Hoos, and K. Leyton-Brown, “Sequential model-based
optimization for general algorithm configuration,” in Proc. of LION-5,
2011, p. 507523.

[21]

[23]

[25]

[26]

[27]

[28]

[29]

K. D. Boese, A. B. Kahng, and S. Muddu, “A new adaptive multistart technique for combinatorial global optimizations,” Oper. Res. Lett.,
vol. 16, no. 2, pp. 101–113, 1994.
T. Jones and S. Forrest, “Fitness distance correlation as a measure of
problem difficulty for genetic algorithms,” in ICGA, 1995, pp. 184–192.
D. S. Oliver and Y. Chen, “Recent progress on reservoir history
matching: A review,” Computational Geosciences, vol. 15, pp. 185–
221, 2010.
F. Zhang and A. Reynolds, “Optimization algorithms for automatic
history matching of production data,” in Proceedings of 8th European
Conference on the Mathematics of Oil Recovery, 2002.
M. Sen and P. Stoffa, “Global optimization methods in geophysicalinversion,” Elsevier Science, 1995.
R. W. Schulze-Riegert, J. K. Axmann, O. Haase, D. T. Rian, and Y. L.
You, “Evolutionary algorithms applied to history matching of complex
reservoirs,” SPE Reservoir Evaluation and Engineering, vol. 5, no. 2,
pp. 163–172, 2002.
B. Davies, “Limits of automatic history matching,” ser. Workshop on
conditioning Reservoir Models to Dynamic Data. Barcelona, Spain:
EAGE/EUROPEC Meeting, June 2010.
A. D. Monfared and A. Helalizadeh, “Automatic history matching using
the integration of response surface modeling with a genetic algorithm,”
Petroleum Science and Technology, vol. 30, no. 4, pp. 360–374, 2012.
C. R. Xavier, E. P. dos Santos, V. da Fonseca Vieira, and R. W. dos
Santos, “Genetic algorithm for the history matching problem,” Procedia
Computer Science, vol. 18, pp. 946–955, 2013.
Y. Hajizadeh, M. Christie, and V. Demyanov, “History matching withdifferential evolution approach; a look at new search strategies,” in SPE
Annual Conference and Exhibition, 2010.
——, “Ant colony optimization for history matching,” in SPE EUROPEC Annual Conference and Exhibition, 2009.
A. Mirzabozog, L. Nghiem, Z. Chen, and C. Yang, “Differential
evolution for assisted history matching process: Sagd case study,” in
SPE Heavy Oil Conference, 2013.
A. P. Reynolds, A. Abdollahzadeh, D. W. Corne, M. Christie, B. Davies,
and G. Willians, “A parallel boa-pso hybrid algorithm for history
matching,” in IEEE Congress on Evolutionary Computation, 2011, pp.
894–901.
F. Verga, M. Cancelliere, and D. Viberti, “Improved application ofassisted history matching techniques,” Petroleum Science and Engineering, vol. 109, pp. 327–347, 2013.
M. A. Christie, “Upscaling for reservoir simulation,” Society ofPetroleum Engineers, 1996.
D. Arnold, V. Demyanov, D. Tatum, M. Christie, T. Rojas, S. Geiger,
and P. Corbett, “Hierarchical benchmark case study for history matching, uncertainty quantification and reservoir characterisation,” Computers and Geosciences, vol. 50, pp. 4–15, 2013.
J. Zhang and A. C. Sanderson, “JADE: Adaptive Differential Evolution
With Optional External Archive,” IEEE Tran. Evol. Comput., vol. 13,
no. 5, pp. 945–958, 2009.
J. J. Liang, B. Y. Qu, and P. N. Suganthan, “Problem Definitions and
Evaluation Criteria for the CEC 2014 Special Session and Competition on Single Objective Real-Parameter Numerical Optimization,”
Zhengzhou Univ. and Nanyang Technological Univ., Tech. Rep., 2013.
S. Das and P. N. Suganthan, “Differential Evolution: A Survey of theState-of-the-Art,” IEEE Tran. Evol. Comput., vol. 15, no. 1, pp. 4–31,
2011.
C. L. Müller and I. F. Sbalzarini, “Global characterization of the CEC
2005 fitness landscapes using fitness-distance analysis,” in EvoApplications, 2011, pp. 294–303.
N. Hansen, S. Finck, R. Ros, and A. Auger, “Real-parameter black-box
optimization benchmarking 2009: Noiseless functions definitions,”
INRIA, Tech. Rep. RR-6829, 2009, updated February 2010. [Online].
Available: http://coco.gforge.inria.fr/bbob2012-downloads
I. Loshchilov, M. Schoenauer, and M. Sebag, “Bi-population CMAES agorithms with surrogate models and line searches,” in GECCO(Companion), 2013, pp. 1177–1184.