Transcript GLOBAL AND LOCAL OPTIMISATION FOR PARAMETER ESTIMATION OF OIL RESERVOIRS S. Gómez and L.

GLOBAL AND LOCAL OPTIMISATION
FOR PARAMETER ESTIMATION
OF OIL RESERVOIRS
S. Gómez and L. Castellanos
UNAM, México
O. Gosselin and T. Fincham
TOTAL-FINA-ELF, GRC-UK
A. Gavana
AGIP, Italy
IMA, Workshop January 2003
Objectives
Develop Global Optimisation Software to have:
• Accurate Production Forecast of oil reservoirs
• Speed and Robustness to get a useful tool
Within the Tunneling method Frame:
Local Optimisation methods
Scaling and regularising
Direct problem
To predict the production of a reservoir
•It is necessary to solve a system of PDE equations
in 3-D, with initial and boundary conditions
to simulate the flux, to obtain, for each production well:
the state variables (presure and saturation)
Bottom-Hole pressure (bhp)
Water cut (wct)
Gas-oil ratio (gor)
•In this work this system is solved by the oil simulator
ECLIPSE (like a black box)
Direct problem
The flux equations depend non-linearly on certain coefficients
p(x,y,z) that characterize the porous media
p(x,y,z) : permeability and porosity of the porous media
We call the solution of the flux equations Mod (p)
•In general the direct problems:
Cause
Effect
p(x,y,z)
pres & sat
Simulation
Forecast
Inverse problem
Characterize the reservoir: Find the properties of the
porous media p(x,y,z), which are not known,
using available information on the solution of the model
Mod(p) (pres. and sat.) measured on the wells
data (x, y, z, t )
Inverse: known
effect
ijk
not-known
cause
In principle p has to satisfy
Mod ( p ) = data
Then, to be able to predict production, we need
to find first the characteristics of the porous media,
permeability and porosity of the whole area, solving
an inverse problem.
Once the inverse problem is solved, we can solve
the flux equations, for future periods of time (prediction) .
These parameters can be obtained solving an
optimisation problem: History Matching
minimising a least squares error criteria or missfit
function, so that the simulation model of the reservoir
approximates the measured data closely
Find the 'optimal' parameters
p : permeability and porosity

1
data Mod p 
min f ( p)  
p
2 j i
2
2
subject to
known bounds on the parameters
pmin  p  pmax

1
data Mod p 
m in f ( p)  
p
2 j i
2
2
the data are the measured values of observable quantities
bottom-hole pressure (bhp)
water-cut
gas-oil ratio
Mod ( p ) is their calculated values from a reservoir
simulation with a certain set of values of parameters p
The index
j
The index
i
runs over the number of wells
of the measurements for each well
 denotes the normalisation factor for a given observable data.
This inverse problem is ill-posed :
I. The problem is non-convex
many local solutions
II. Non-continuous dependence on the data:
Regularisation is needed to get the best
possible approximation, preventing unbounded
magnification of the data-error
The parameterization of p(x,y,z) adds
uncertainty in the problem formulation
Regularisation
•Tikhonov, accurate solutions only with optimal parameter 
and to find is expensive. Theory well developed for linear
problems, but for non-linear it is still an open question
•Multiscale Optimisation (G. Chavent), very good
numerical results in seismic inversion (Chavent, Clement,
Gomez), and in water reservoirs (Gomez, Perez, Alvarez)
but difficult to implement with Tunneling + black box simulator
•Here, Conjugate Gradient + L-curve to get stable directions
within the T-Gauss-Newton Local method.
In order
• to get many local optimal solutions
•to deal with non-continuity
• to deal with uncertainty
THE GOAL is to get the set of stable optimal
solutions that produce a good match to the data :
alternative SCENARIOS of production
Global Optimisation
Local optimisation methods only find one
local minimum that may not fit the data or
that may not give the best fit.
We use our global optimisation method:
The Tunneling Method
f(x)
tunnel
tunnel
*
1
x
*
2
x
*
3
x
x
The Tunneling Method
• obtains a sequence of local minima with
decreasing values of the missfit function
   
 
f x  f x  f x
*
1
*
2
*
G
• It has two phases:
Minimisation Phase .- Finds a local minimum x*
Tunnelisation Phase .- Finds a point x in another valley
Minimisation phase
Given x0 find a local minimum x* with f*
We can use any local optimisation method
( gradient - based ) with bound constraints
Here we test
 Limited Memory Quasi-Newton LBFGSB (Nocedal et al)
 Truncated Gauss-Newton TRON (More et al)
Tunneling Phase
Once a local minimum x* has been obtained with f *,
to be able to tunnel from one valley to another,
that is
Find x0  B such that
T( x0 ) = f ( x0 ) - f *  0
using the same local gradient-type method, it
1) creates a pole to destroy the local minimum
2) solves an inequality problem T(x)  0
The Tunneling Phase
f(x)
f*
x*
f ( x)
x
f*
f ( x)
x
f *
0
Creates a pole to destroy the minimum
f ( x)
f*
f *
f ( x)
0
x
T (x)
T (x)

0
= x° solution
X*
x
To place a pole at x* we use one of the following
1. Classical Tunneling function,
Tc (x) 
f (x)  f (x*)
x  x*
*
0
2. Exponential Tunneling function,
 * 
Te ( x )   f ( x )  f ( x*)  exp
 0
 x x* 
Where * is the strength of the pole
 . is the squared euclidean norm
Minima at the same level
T ( x) 
f ( x)  f *
l
 || x  x* ||
i 1

To find minima within an interval of F*
F*
x*
Exponential Tunneling
effect of lambda : * = (0.1, 0.4, 0.9)
T(x)
50
40
30
20
10
x
0
-10
-3
-2
-1
0
1
2
3
Effect of 
15
10
5
0
-5
-10
-3
-2
-1
0
1
2
Classical Tunneling
Exact minimum at x= -0.127153
* = 1(blue), 2, 3
3
Effect of  : When the minimum is not exact x* = -.2
15
10
5
0
-5
-10
-3
-2
-1
0
1
2
3
For 1, in the direction where f(x) < f* the pole does not destroy
the minimum, a descent direction is not possible for T.
* has to be increased.
It has already shown to be very efficient to solve highly
difficult problems, both academic and real, such as
• optimal control of chemical process
• optimal structure of molecules (GRBS 94, 97)
• chemical multiphase equilibrium (NGL 2001 a & b)
• polimer reactors (D 2003)
It does not find local minima with larger objective function
than the best founded so far f(x*). This characteristic
makes it faster than other global methods like
Simulated annealing (NGL 2001)
Genetic algorithms (GCSQ 2001)
Tunneling in parallel
• After a local minimum x* has been found,
start the tunneling search from a point in
the neighborhood of x*
x = x* + r
Sequential :
r = random direction
if the search is unsuccesful, start
again from another initial point x
Parallel: start the search in parallel
in random directions (GCCS 2002)
Tunneling in parallel
Why LM-Quasi Newton
vs T-Gauss-Newton
•Our least squares problem could be highly non-linear
( there is an inverse operator in f )
•LM Quasi-Newton approximates the inverse of the
complete Hessian
JT J 

k 1,m
rk  2 rk
•T Gauss-Newton uses the approximation
JT J
with J given by the simulator.
We test if this approximation works for our problem
Limited Memory Quasi-Newton: with Linear search
Stopping Conditions:
   
 FRTOL
m ax f x  , f x  ,1
f x k  f x k 1
k 1
k
 
proj g x k

 PGTOL
Truncated Gauss-Newton: with Trust region (TRON)
Stopping Conditions:
 

 

f x k  f x k 1
 FRTOL
k
f x
 
g xk
2
 
 PGTOL * g x 0
2
Another reason to use Truncated Gauss-Newton
To find the descent direction s , solve the linear system:
(J
k 1 T
) J
k 1 k
s  J
k 1
r
k 1
using the Conjugate Gradient method
Stopping the CG iteration before the error in the matrix and
in the right hand side propagates, has a regularisation effect.
To know when to stop, we use our algorithm to find the
corner of the L-curve automatically (CGG 2002). The Trust
region size is modified accordingly.
TRON has the option to precondition with Incomplete Cholesky
Scaling
Characteristic of the History Matching problem:
Very flat valleys
•When an ill-posed problem is discretised, the
resulting Hessians are ill-conditioned
Scaling to make valleys sharper
GSCQ 2001
Scaling
If the original bounds are in
xi  li , ui 
then a new variable y will be in

yi  lin , uin
x  D yc
where D is a diagonal matrix with elements
and c is a vector with
u i  li
u in  l in
li uin  ui lin
uin  lin
•The function in y is not altered h(y) = f(x), but
the gradient g y  D g
and the Hessian H y  D H D

Scaling
This implies that the Truncated Gauss-Newton
linear system is preconditioned
D J J s  D J r
T
T
•Other traditional preconditioners like ICC may not work:
the singular values of the Hessians decay to zero.
The PUNQ Test Case
•The PUNQ problem is a benchmark for History Matching
and risk-analysis methodologies.
•It is a dynamical reservoir model based on a real
West Africa field, which has been discretised using a
19x28x5 corner-point grid, with 1,721 active cells.
•A history period, simulating 8 years of production from
six wells located close to the gas-oil contact (GOC),
was generated by using geostatistical distributions of
porosities and permeabilities.
•Gaussian noise has been added to the collected well data
to reproduce a real measurement process.
•Then, 8 years of forecast with five additional infilling wells
have been simulated.
The PUNQ Test Case
•The choice of history matching parameters are
pore volume, and horizontal transmissibility multipliers.
• The complexity of the H-M identification problem
was avoided by adopting a set of parameters, based on
the Gradzone (sensitivity) Analysis (GRC-UK), in which he
available a-priori geological information was included.
•This analysis leads to 30 history matching parameters,
which contains a signature of the geological model. We also
use a set of 10 history matching parameters, 5 and 5 :
one multiplier for each property was assigned to every layer.
The PUNQ Test Case
A multiplier parameter means that although the gridcells in
one zone can have a different pore volume, their values
relative to each other remain constant. This way to
parameterize the problem is a source of uncertainty.
To restrict the evolution of the system into a physically
reasonable region, simple bounds, acting as perfectly
absorbing surfaces, are imposed to the parameters :
0.1  xi  3.0, i 1,...,n
Original results on a simpler synthetic case without data-error
Phase
f*
fn-ev
Min1
Tun1
Min2
Tun2
Min3
Tun3
Min4
Tun4
Min5
0.0936
0.0887
0.0867
0.0742
0.0046
0.0045
0.0035
0.0010
0.0009
107
36
44
76
65
178
34
128
32
•Total cost
x*
fmax
1.011 0.986 1.0134 0.115 1.366
6 E+3
0.998 1.002 1.0820 0.097 0.763
1 E+4
0.991 1.001 1.0800 1.024 0.752
5 E+2
1.004 0.995 0.9280 0.797 1.076
7 E+2
1.007 0.991 0.9560 0.931 1.201
600 / 600
fn / grad
PUNQ-10 without noise
initial point = ip
Method
# Min
x* = (1, 1, ..., 1)
fun / grad
f*
frtol
Tun-TGN no-prec
2
30 / 21
0.00002
0.5x10-3
Tun-TGN
1
14 / 10
0.0000005

Tun-LBFGSB
5
166 / 166
0.0020
2.2x10-3
Tun-LBFGSB
2
229 / 229
0.0001
2.2x10-6
prec
No-error, the global min (zero residual) gives
the best approximation
PUNQ-10 without noise
initial point = ub
Method
# Min
x* = (1, 1, ..., 1)
fun / grad
f*
frtol
Tun-TGN no-prec
1
23 / 23
0.00005
0.5x10-3
Tun-TGN
1
27 / 19
0.00005

Tun-LBFGSB
10
408 / 408
0.00217
2.2x10-3
Tun-LBFGSB
4
483 / 483
0.0001
2.2x10-6
prec
PUNQ-10 with noise
initial point = (1, 1, ..., 1)
x* = ?
Method
# Min
fun/grad
Tun-TGN+noprec+ no-sca
4
211 / 105
0.309
0.5x10-3
Tun-TGN+
6
493 / 302
0.329
“
9
440 / 281
0.295
“
2
248 / 160
0.315
“
Tun-TGN+noprec+ sca [0.7,1.3]
5
278 / 170
0.288
“
prec+ sca [0.7,1.3]
3
387 / 236
0.317
“
Tun-LBFGSB
7
873 / 873
0.343
2.2x10-3
Tun-LBFGSB
7
599 / 599
0.317
2.2x10-6
prec + no-sca
Tun-TGN+noprec+ sca[1, 2]
Tun-TGN+
Tun-TGN+
prec+ sca[1,2]
f*
frtol
GOR PRO1
200
150
Not-scaled
100
50
0
0
5
10
15
20
Years
minimum 1
minimum 2
minimum 3
minimum 4
actual
GOR PRO-1
250
GOR (sm3/sm3)
GOR (sm 3/sm 3)
250
200
150
Scaled [0.7, 1.3]
100
50
0
0
5
10
15
20
Years
minimum 1
minimum 2
minimum 3
minimum 4
minimum 5
actual
The best minimum ?
GOR PRO-4
Well WCT PRO-11
1
160
0.9
140
GOR (sm 3/sm 3)
0.8
0.7
0.6
0.5
0.4
120
100
80
60
0.3
40
0.2
20
0.1
0
0
0
0
5
10
15
5
20
10
15
20
Years
Years
Minimum 1
Minimum 2
Minimum 3
Minimum 4
Minimum 5
actual
minimum 1
minimum 2
minimum 3
minimum 4
minimum 5
actual
•WCT is the most sensitive
(to error propagation?, to the formulation?):
an intermediate minimum gives the best approx.
•GOR the least sensitive
Extreme Cases
GOR PRO-15
200
80
70
60
50
40
30
20
10
0
180
GOR (sm 3/sm 3)
GOR (sm 3/sm 3)
GOR PRO-5
160
140
120
100
80
60
40
20
0
5
10
15
20
0
0
Years
minimum 1
minimum2
minimum3
5
10
15
20
Years
minimum4
minimum5
actual
Perfect cases :
GOR Pro1, Pro11, Pro12
minimum 1
minimum 2
minimum 3
minimum 4
minimum 5
actual
The only bad case
Scenarios: BHP
BHP - PRO-4
BHP - PRO-12
260
260
240
220
200
180
160
140
120
100
240
220
200
180
160
140
120
0
5
10
15
20
100
0
minimum 1
minimum2
minimum3
minimum4
minimum5
actual
2
4
minimum 1
6
minimum2
8
10
minimum3
12
minimum4
14
minimum5
16
18
actual
BHP - PRO-1
BHP - PRO-15
260
240
240
220
220
200
200
180
180
160
160
140
140
120
120
100
100
0
2
4
minimum 1
6
minimum2
8
minimum3
10
12
minimum4
14
minimum5
16
actual
18
0
5
minimum 1
minimum2
10
minimum3
15
minimum4
minimum5
20
actual
Forecast for the Field Properties
FWPT: Field Total Water Production
Punq10, FWPT
900000
800000
REAL
TUN_1
TUN_2
TUN_3
TUN_4
TUN_5
700000
FWPT (SM3)
600000
500000
400000
300000
200000
100000
0
8
9
10
11
12
13
TIM E (Years)
14
15
16
17
Forecast Field Gas Total Production
Punq10, FGPT
350000
REAL
TUN_1
TUN_2
TUN_3
TUN_4
TUN_5
300000
FGPT (SM3)
250000
200000
150000
100000
8
9
10
11
12
13
TIM E (Years)
14
15
16
17
Forecast Field Oil Total production
Punq10, FOPT
4000000
3500000
REAL
TUN_1
TUN_2
TUN_3
TUN_4
TUN_5
FOPT (SM 3)
3000000
2500000
2000000
1500000
8
9
10
11
12
13
TIM E (Years)
14
15
16
17
The error between the real case and the
calculated forecast for the 5 minima found.
The error has been calculated on 6 field properties
FOPT, FGPT, FWPT, FGOR, FGPR, FWCT
Min1
Min2
Min3
Min4
Min5
FOPT
FGPT
FWPT
FGOR
FGPR
FWCT
ERROR
6.241
8.620
8.226
8.483
5.202
6.655 7.449
7.716 11.150
7.790 10.356
9.140 9.652
5.016 6.157
6.679
8.476
8.957
9.071
5.573
6.275
7.085
7.589
9.152
4.830
6.265
7.645
8.294
8.577
5.315
6.594
8.449
8.535
9.013
5.349
Global Optimisation as a regularisation strategy ?
•The degree of approximation ( a local minimum with good
match to the data and forecast ), would be the regularisation
parameter.
•Then, if the data-error   0,
p  p
*
k
*
G (zero-residual)
•If there is no error in the data, the global minimum is the
zero-residual exact fit and gives the best approximation
•If there is error, an intermediate minimum gives
the best possible approximation (H-M & Forecast)
Conclusions
•The Tunneling Method + TGN has proven to be capable of
finding improved solutions
•Tunnel-TGN is much faster than Tunnel-Lbfgsb:
The number of function evaluations can be reduced up to 4
times, and gradient evaluations up to 6 times.
•Tunnel-TGN method, is attracted to the global minimum
when there is no-noise in the data, from some starting points,
whereas Tunnel-Lbfgsb is attracted to several local minima,
from every starting point.
Conclusions
•The Tunnel-TGN method, is able to find a set of local minima
with good HM, on a reasonable amount of time.
•These good local solutions, produce a set of scenarios of
production, where each production variable is well predicted
by one of these minima.
•Scaling (Diagonal Preconditioner) clearly improves the
forecast for noisy data
Conclusions
•The sensitivity of each state variable (pres. & sat.)
to the degree of attainable parameter identification
(the best approximation) may not be the same :
different local minimum may forecast better different
state variable
Only the scenarios
(given by local minima with good H-M),
can be useful as a Forecast decision tool
The code TUNEL we are using in this work
has recently been developed:
•Castellanos-Gómez 2000 - 2002
It has the option to use any local method and to scale
•It can be used to solve ANY global optimisation
problem, if the objective function is continuously
differentiable and the gradient is available
Future
•The regularisation properties of the Conjugate
Gradient iteration within the non-linear Newton
deserves further research
•The Parallel Tunneling method should make the
forecast much faster. The code is ready.
susanag @ servidor.unam.mx
Truncated Gauss-Newton
1
2
m in f x   r x  2
2
s.t. xmin  x  xmax
x x
k
k 1
r x  : R n  R m
s
k
1 k 1
k 1
s  arg min m s   J s  r
2
k
Trust Region:
k
s.t .
s k  k 1
2
2