An efficient method for smart well production optimisation

Journal of Petroleum Science and Engineering 69 (2009) 25–39

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / p e t r o l

Research paper

An efficient method for smart well production optimisation Daniel Chr. Doublet a, Sigurd I. Aanonsen b,⁎, Xue-Cheng Tai c a

Centre for Integrated Petroleum Research and Department of Mathematics, University of Bergen, Allégaten 41, N-5007 Bergen, Norway Centre for Integrated Petroleum Research and Department of Mathematics, University of Bergen, Allégaten 41, N-5007 Bergen, Norway Currently with StatoilHydro, Sandsliveien 90, N-5254 Sandsli, Norway c Department of Mathematics and Centre for Integrated Petroleum Research, University of Bergen, Allégaten 41, N-5007 Bergen, Norway b

a r t i c l e

i n f o

Article history: Received 10 January 2008 Accepted 7 June 2009 Keywords: production optimization smart wells augmented Lagrangian method

a b s t r a c t A method for dynamic optimisation of water flooding with smart wells is presented. The method finds the optimal injection and production rates for every well segment of the smart wells. We formulate the problem as a constrained optimisation problem and state this problem as an augmented Lagrangian saddle point problem. Comparisons are made with a more traditional optimal control method based on solving the adjoint systems of equations. In the examples tested the method obtains same maximum profit while using less computational effort. © 2009 Elsevier B.V. All rights reserved.

1. Introduction As the oil resources of the world are becoming increasingly difficult to recover, it has become more important to produce existing fields as efficiently as possible and to decrease the development and operating costs. This problem is an optimisation problem where we want to maximise some profit function. In this paper we propose a method for maximising the net present value (NPV) of an oil reservoir, by reducing water production and increasing oil recovery at the same time as we are delaying water breakthrough. Optimal control theory methods have been used by many authors to solve this problem using the adjoint method, see e.g. (Asheim, 1988; Brouwer and Jansen, 2004; Lien et al., 2006; Pallav et al., 2005; Ramirez, 1987; Sudaryanto and Yortsos, 2001; Zakirov et al., 1996). However, in all these approaches it is required to solve the state equations exactly for each new estimate of the control variables, which is computationally expensive. In this paper we formulate the optimisation problem as an augmented Lagrangian saddle point problem, and present a method for solving it by solving the Karush-Kuhn-Tucker (KKT) conditions for the augmented Lagrangian functional. By solving the KKT conditions sequentially, we avoid solving the state equations exactly for each new estimate of the control variables, thereby reducing the computational cost for finding a new estimate of the controls. Although the state equations are not fulfilled during the optimisation procedure, they will be so at convergence. Preliminary results have been presented earlier in (Doublet et al., 2006) and (Doublet et al., 2007). The use of the augmented Lagrangian functional for constrained optimisation, with, was introduced by Hestenes in (Hestenes, 1969) and

⁎ Corresponding author. E-mail address: [email protected] (S.I. Aanonsen). 0920-4105/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.petrol.2009.06.008

has been applied to optimal control problems later in (Ito and Kunisch, 1990; Kunisch and Tai,1997). In this paper we make use of the augmented Lagrangian functional, but the solution approach is somewhat different than the one used in (Ito and Kunisch, 1990) and (Kunisch and Tai, 1997). A proof of convergence for a quadratic optimisation problem with linear constraints is presented in the appendix, showing that it is possible to find a penalisation parameter that guarantees convergence in the quadratic case. The efficiency of the method is demonstrated through several numerical examples. We also present comparisons with the more traditional adjoint method.

2. Problem formulation We consider a rectangular, heterogeneous, two-dimensional, twophase (oil and water) reservoir with no-flow boundaries inspired by the model of Brouwer and Jansen (Brouwer and Jansen, 2004). The model is horizontal such that gravitational effects can be ignored. The capillary pressure is zero, and we use the Corey models to define the relative permeabilities. There is an injector along the left edge of the reservoir, and a producer along the right edge of the reservoir, as shown in Fig. 1. The two smart wells have several segments, indicated by black dots, such that the injection and the production can be controlled individually for each segment. Initially the reservoir is completely oil saturated, and at the start of operation water is injected in the segments on the left hand side of the reservoir, whilst we are producing from the segments at the right hand side of the reservoir. Initially, only oil is produced, but after a certain time we will start to produce both oil and water, and finally only water is produced. There is a profit associated with production of oil. However when producing both oil and water simultaneously we must remove the water from the oil, resulting in a cost associated with water production. The total production rate is fixed,

26

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

to number Ninj are assumed to belong to the injection well, while the segments from number 1 + Ninj to number Nprod + Ninj are assumed to belong to the production well. Moreover, the control variables must also satisfy the following inequality constraints, for n = 1,…,N, n

0≤vi ≤1; for i = 1; …; Ninj + Nprod :

ð4Þ

Since only water is injected, the liquid rate equals the water rate for the segments in the injector. Thus we have that, for n = 1,…,N, n

n

Vwi = vi Q; for i = 1; …; Ninj ;

ð5Þ

n is the water rate at injection segment i at time step n. In the where Vwi segments of the producer, however, the liquid rate equals the sum of the water and oil rates so that, for n = 1,…,N,

n

and we assume that we inject at the same rate as we produce. The profit may be controlled by adjusting the percentage of total injection/ production in each of the segments, and our objective is to maximise the NPV from an oil reservoir by adjusting the individual segments injection and production rates at distinct times. The reservoir flow equations are discretised using a standard cell centred grid with upstream weighting (Aziz and Settari, 1979) using backward Euler to approximate the time derivative. From the discretisation we derive a discrete time model of the two-phase conservation equations, n

n

n

e ðv̂ ; û ; û

n−1

Þ = 0; for n = 1; N N;

ð1Þ

where en − {eno,enw}T is the residual column vector corresponding to the discretised version of the two-phase flow equations (Aziz and Settari, 1979), the superscript n denotes the discrete time step, N is the total T T number of time steps, (ûn)T = {pn , Sn } is a column vector consisting of the pressures and the water saturations in all the grid cells at time step n, and v̂n is the column vector consisting of the control variables at time step n. The elements of v̂n are related to the water injection and liquid production rates in the different segments of the wells. Instead of using a well model, we will directly control water injection and liquid production at each segment such that the well segments are treated as source and sink terms. We assume that the total water injection rate equals the total liquid production rate and denote the total injection/production rate by Q. Letting vni denote the percent of total injection or production in segment i at time step n, we have the relations

Ninj

n

∑ vi = 1; for n = 1; …; N;

i=1

ð2Þ

and Nprod + Ninj

∑

i = 1 + Ninj

n

vi = 1; for n = 1; …; N;

n

n

Vwi + Voi = −vi Q; for i = 1 + Ninj ; …; Nprod + Ninj ;

Fig. 1. Reservoir.

ð3Þ

where Ninj is the number of segments in the injector and Nprod is the number of segments in the producer. All the segments from number 1

ð6Þ

n n and Voi are the water and oil rates, respectively, at where Vwi production segment i at time step n. The different phase production rates can be expressed as functions of the liquid rate and the fractional flow at the well segment so that

n

λwi n v Q for i = 1 + Ninj ; …; Nprod + Ninj ; + λnoi i

ð7Þ

λnoi n v Q; for i = 1 + Ninj ; …; Nprod + Ninj ; λnoi + λnwi i

ð8Þ

n

Vwi = −

λnwi

and n

Voi = −

where λnwi and λnoi are the water and oil mobilities, respectively, in the grid cell containing well segment i and at time step n. 2.1. Profit function Our aim is to maximise net present value, by controlling the individual segment rates during the entire production period. The net present value (NPV), J, is given as N

n

n

n

Jðv̂ ; û Þ = ∑ J ðv̂ ; û Þ; n=1

ð9Þ

with

n

n

n

"N

J ðv̂ ; û Þ = ΔxΔyh

# n n −Iw ⋅Vwi −Io ⋅Voi n Δt ; tn i = 1 + Ninj ð1 + b=100Þ prod

+ Ninj

∑

ð10Þ

where Vo and Vw are yearly produced volumes of water and oil, Io and Iw are the revenue of oil and the cost of water produced per volume expressed in /m3, respectively, b is the annual interest rate expressed in %, Δx and Δy are the dimensions of the grid cells in, respectively, horizontal and vertical direction, Δt n is the size of the nth time step, and tn = Σni = 1Δti is the time expressed in years at time step n. Since n n and Voi , are less than or equal the production rates of water and oil, Vwi to zero, Io is a positive constant and Iw is a negative constant. The objective is to maximise the function (9) by adjusting the percentage of total injection and production in each of the individual segments of the injection well and the production well, respectively. If we at some point produce so much water that J n in Eq. (10) is negative, then we set J n to zero. This is done to maximise the profit over the field life

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

time, while still letting the profitable production period be unknown. It should be noted that the method will reduce production from waterproducing connections to zero, or almost zero, when this is beneficial, and increase production from other connections correspondingly.

27

for n = 1,…,N. Since the last of these conditions yields en = 0, for n = 1,…, N, it is easy to see that the KKT conditions also hold for the augmented Lagrangian functional ⋆ ⋆ ⋆ ∂ Lc ðv̂ ; û ; η Þ = 0; ∂û n

2.2. A saddle point formulation inj + Nprod , Let us define the control variables at time step n as v̂n ={vni }N i=1 and let v̂={v̂n}nN= 1 and û={ûn}nN= 1 . We can now formulate our problem as a constrained minimisation problem where we want to solve

min −Jðv̂; û Þ

ð11Þ

v̂;û

ð18Þ

∂ ⋆ ⋆ ⋆ L ðv̂ ; û ; η Þ = 0; ∂v̂ n c ∂ ⋆ ⋆ ⋆ L ðv̂ ; û ; η Þ = 0; ∂ηn c

subject to n

n

n

e ðv̂ ; û ; û

n−1

Þ = 0; for n = 1; …; N:

ð12Þ

The corresponding Lagrangian is N

n

n

n

Lðv̂ ; û ; ηÞ = ∑ L ðv̂ ; û ; û

n−1

; η Þ;

n

n

n=1

n

ð13Þ

for n = 1,…,N. In order to find the solution of Eq. (11), we can thus solve the KKT system Eq. (17) or equivalently Eq. (18). In the recent years this problem has been solved using the adjoint method, see e.g. Brouwer and Jansen (Brouwer and Jansen, 2004). However, when solving this problem using the adjoint method, the forward problem, Eq. (1), has to be solved exactly for each new gradient. Solving the forward problem is time consuming, and we present here a method which does not require that the forward problem is solved exactly at each iteration. Thus we propose an algorithm for solving Eq. (11), Eq. (12) by solving the system of Eq. (18).

where 3. Optimisation method n

n

n

L ðv̂ ; û ; û

n−1

n

n

nT n

n

n

; η Þ = −J ðv̂ ; û Þ + η e ðv̂ ; û ; û

n−1

Þ;

ð14Þ

and η = {ηn}nN= 1 are the Lagrange multipliers. Furthermore, we define the augmented Lagrangian by N

Lc ðv̂ ; û ;ηÞ = ∑

n=1

n n n n−1 n Lc ðv̂ ; û ; û ; η Þ;

ð15Þ

Letting the subscript k be the outer iteration counter, we propose the following algorithm to solve the KKT conditions Eq. (18). Algorithm 1. KKT Optimisation Algorithm 1. Choose v̂0 = {v̂n0}nN= 1, û0 = {ûn0}nN= 1, c N 0. For k = 1,2,… do: 2. Find ηN k such that ∂

N

∂û N

where

N

N

N−1

N

Lc ðv̂ k−1 ; û k−1 ; v̂ k−1 ; ηk Þ = 0:

ð19Þ

3. For n = N − 1, N − 2,…1, find ηnk , such that n n n n−1 n Lc ðv̂ ; û ; û ;η Þ

n

n

n

= L ðv̂ ; û ; û

n−1

n

;η Þ

ð16Þ

c n n n n−1 T n n n n−1 + e ðv̂ ; û ; û Þ e ðv̂ ; û ; û Þ; 2

∂ n L ∂û n c +

and c N 0 is a penalisation constant. It is known that the solution of Eq. (11), Eq. (12) is a saddle point of Eq. (15) for a sufficiently large penalisation parameter c, see for example Bonnans et al. (Bonnans et al., 2006) pp. 275–279 for proof. Given that the set of constraint gradients is linearly independent at the solution of Eq. (11), Eq. (12), the following conditions hold at the solution of Eq. (11), Eq. (12), see e.g. Nocedal and Wright (Nocedal and Wright, 1999). 2.2.1. Karush Kuhn Tucker (KKT) conditions Suppose that {v̂⋆,û⋆} is a solution of Eq. (11). Then there exists a set of vectors η⋆ such that the following conditions hold at the point {v̂⋆, û⋆},

⋆ ⋆ ⋆ ∂ Lðv̂ ; û ; η Þ = 0; ∂û n ∂ ⋆ ⋆ ⋆ Lðv̂ ; û ; η Þ = 0; ∂v̂

∂ ⋆ ⋆ ⋆ Lðv̂ ; û ; η Þ = 0; ∂ηn

+ 1

n + 1

n + 1

n

n + 1

ðv̂ k−1 ; û k−1 ; û k−1 ; ηk

ð20Þ

Þ

∂ n n n n−1 n L ðv̂ ; û ; û ; η Þ = 0: ∂û n c k−1 k−1 k−1 k

4. Then, for n = 1,…,N, find v̂nk such that n

n

n

n

n−1

n

v̂ k = arg min Lc ðv̂ ; û k−1 ; û k−1 ; ηk Þ: n

ð21Þ

v̂

5. Then update ûk − 1 = {ûnk − 1}nN= 1 by n

!−1 ∂2 ∂ L ðv̂ ; û ; η Þ L ðv̂ ; û ; η Þ: ∂ηn c k k−1 k ∂ηn ∂û n c k k−1 k

n

û k = û k−1 −

ð22Þ

After finding ηk from Eqs. (19) and (20) in the algorithm, the first KKT condition is fulfilled. By studying the augmented Lagrangian given ð17Þ

by Eqs. (15) and (16) we observe that ∂ ∂û

n

Lc = N

∂ ∂û

and η and

n

∂

∂

N

Lc =

∂ ∂û

N

LNc and that

Lc is dependent on v̂N, ûN, ûN − 1

Lnc

+

∂

Lc , for n = 1,…,N− 1, is dependent on v̂n + 1, v̂n, ûn + 1, ûn,

∂û

n

∂û

n

Lnc + 1 ,

∂ ∂û

so that

∂û

N

ûn − 1, ηn + 1 and ηn. Thus the first KKT condition can be solved exactly ∂ for η if we start by finding ηN such that N Lc = 0 and then, for n =N− 1, ∂û

28

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

N − 2,…,1, we find ηn such that

∂ ∂û

n

Lc = 0. Furthermore we observe that,

when calculating v̂k from Eq. (21) in the algorithm above, we find v̂k such that

n v̂ i

∂ L ðv̂ ; û ; η Þ = 0; ∂v̂ n c k k−1 k for n = 1,…,N. Thus the second KKT condition is fulfilled after performing the calculations in Eq. (21). Finally in Eq. (22) of the algorithm, we are solving the third KKT condition by a Newtoniteration. In this way, the third KKT condition is fulfilled at convergence, and since we are enforcing the first and the second of the KKT conditions in Eqs. (19)–(21) the solution is found at convergence. The algorithm then becomes: Algorithm 2. KKT optimisation algorithm, in detail v̂0 = {v̂n0}nN= 1,

1. Choose 2. Find ηkN such that −

∂J N ∂û N

û0 = {ûn0}nN= 1,

c N 0. For k = 1,2,… do:

N

N

ðv̂ k−1 ; û k−1 Þ

N

N

N

N

N−1

T

+ ðηk + c⋅e ðv̂ k−1 ; û k−1 ; û k−1 Þ

∂eN ∂û

N

N

N−1

ð23Þ

ðv̂ k−1 ; û k−1 ; û k−1 Þ = 0; N

n ∂J n n n n n n n n−1 T ∂e n n n−1 ðv̂ ; û Þ + ðηk + c⋅e ðv̂ k−1 ; û k−1 ; û k−1 ÞÞ ðv̂ ; û ; û Þ ∂û n k−1 k−1 ∂û n k−1 k−1 k−1 n + 1

+ ðηk

+ c⋅e

n + 1

n + 1

n + 1

n

n

T

ðv̂ k−1 ; û k−1 ; û k−1 û k−1 ÞÞ

4. Then, for n = 1,…,N, we will find v̂nk , while holding ûnk − 1 constant, such that   T    c n n n n n n n n n−1 n n n n−1 v̂ k = arg min −J v̂ ; û k−1 + ηk + e v̂ ; û k−1 ; û k−1 e v̂ ; û k−1 ; û k−1 n 2 v̂

ð25Þ

n

3.1. Comparison with the adjoint method To evaluate the KKT algorithm, we compare with an optimal control theory method, the adjoint method, used previously in (Brouwer and Jansen, 2004; Lien et al., 2006; Sarma et al., 2005) and (Zakirov et al., 1996). Letting, as for the KKT algorithm, the subscript k be the outer iteration counter, the adjoint algorithm is as follows. Algorithm 3. The adjoint algorithm

n  −1   ∂e  n n n−1 n n n n−1 e v̂ k ; û k−1 ; û k−1 : n v̂ k ; û k−1 ; û k−1 ∂û

−

N N    ∂J  N N N T ∂e N N N−1 = 0: v̂ k−1 ; û k + ηk v̂ k−1 ; û k ; û k N N ∂û ∂û

ð29Þ

4. For n = N − 1, N − 2,…,1, find ηnk , such that −

n n   ∂J  n n nT ∂e n n n−1 n v̂ k−1 ; û k + ηk n v̂ k−1 ; û k ; û k ∂û ∂û n + 1  n + 1 n + 1 n n + 1T ∂e + ηk v̂ k−1 ; û k ; û k = 0: n ∂û

ð30Þ

ð26Þ

In this step Eq. (26) is evaluated by inserting ûk − 1 and v̂k directly into the flow (simulator) equations without actually running the reservoir simulator. When finding ηn, for n − 1,…,N in Eqs. (23) and (24) of the algorithm we need to solve a linear system of equations, and we do this with the GMRES method, as described in for example Golub and van Loan (Golub and van Loan, 1996) or Saad (Saad, 2000).

Table 1 Estimation of Cost per Iteration, where Mc is the number of grid cells and Mw = Ninj + Nprod is the total number of segments of injection and production type. KKT algorithm N linear systems of equations of size 2·Mc N non-linear systems of equations of size Mw N linear systems of equations of size 2·Mc Adjoint algorithm Eqs. (29) and (30) Eq. (28) Eq. (31)

ð28Þ

3. Find ηN k such that



Eqs. (23) and (24) Eq. (25) Eq. (26)

ð27Þ

1 where Hi = B− i , and Bi is the limited memory BFGS matrix as described in (Byrd et al., 1994, 1995). Finally, for the update in Eq. (26), we use as in Eqs. (23) and (24), the GMRES method to solve the linear system of equations.

5. And finally for n = 1,…,N update the state variables ûn by n



n ∂L  n n n−1 n = v̂ i −Hi cn v̂ ; û k−1 ; û k−1 ; ηk ; ∂v̂

n + 1

∂e n + 1 n + 1 n ðv̂ k−1 ; û k−1 ; û k−1 Þ = 0; ∂û n

ð24Þ

û k = û k−1 −

+ 1

1. Choose v̂0 = {v̂n0}nN= 1, û00 and c N 0. For k = 1,2,… do: 2. For n = 1,…,N, find ûnk such that   n n n n−1 = 0: e v̂ k−1 ; û k ; û k

3. For n = N − 1, N − 2,…,1, find ηnk , such that −

For the minimisation in Eq. (25), we use the LBFGS method which is a limited memory quasi-Newton method, see e.g. Byrd et al. (Byrd et al., 1994, 1995). This is done by finding a new estimate of v̂kn as

N linear systems of equations of size 2·Mc N non-linear systems of equations of size 2·Mc Updating the inverse of the LBFGS second derivative

5. Finally we find the gradient of −J as dJ  n  ∂Lc  n n v̂ v̂ ; û = dv̂ n k−1 ∂v̂ n k−1 k n  T ∂J n  n n n ∂e n n n−1 = 0; = − n v̂ k−1 ; û k + ηk n v̂ k−1 ; û k ; û k ∂v̂ ∂v̂ ð31Þ n

−

and use this gradient in athe LBFGS algorithm to find v̂k = {v̂nk }nN= 1. Eq. (28) in the algorithm is in fact the same as solving the forward problem, and when solving the non-linear systems of equations in Eq. (28), we use Newton's method so that we find ûnk as • Choose ûnk,0. For l = 1,2,… do until convergence: For n = 1,2,…,N do: n

n

û k;l = û k;l−1 −



 −1   ∂en  n n n n n n−1 e v̂ k−1 ; û k;l−1 ; û k;l ; n v̂ k−1 ; û k;l−1 ∂û

ð32Þ

where we solve the linear system of equations with the GMRES method. In order to find ηn for the adjoint algorithm we need to solve a linear system of equations in Eq. (29) and in Eq. (30), for which we use the GMRES method.

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

3.2. Estimation of cost per iteration Let us denote the number of grid cells by Mc, and let Mw = Ninj + Nprod denote the total number of segments of injection and production type. Now we investigate the computational cost of our two algorithms, and a summary is given in Table 1. In practice, the computational time required to perform the calculations in Eq. (25) in the KKT algorithm is negligible compared to that of Eq. (23), Eq. (24) and Eq. (26), since Mw bb Mc. Thus the dominating factor is the time required to solve 2·N linear systems of equations of size 2·Mc. For the adjoint method, the calculation of the derivative and the approximation to the Hessian in Eq. (31) is negligible in comparison with the work done in Eqs. (28), (29) and (30). In Eqs. (29) and (30) we are solving N linear systems of equations of size 2·Mc, and in Eq. (28) the computational time depends upon how many iterations of Newton's method which is required in order to solve the forward problem. For each iteration of Newton's method, we must solve N linear systems of equations of size 2·Mc. Our experience shows that one needs to do between 3 and 5 iterations of Newton's algorithm in order to reach convergence, giving a total of 4 N to 6 N linear systems of equations of size 2·Mc to solve in each outer loop of the algorithm. From this, one sees that the KKT algorithm is from two to three times faster per iteration compared to the adjoint algorithm. Furthermore it must be noted here that by one iteration of the adjoint algorithm we mean one iteration of the algorithm as it is stated in this paper.

29

One can ensure that the inequality constraints (35) are fulfilled for the variables vn1, vn2,…,vnNinj − 1 and vnNinj + 1, vnNinj + 2,…,vnNinj + Nprod − 1 n= 1,…, N by using a projected gradient method in the optimisation, but for the variables vnNinj and vnNinj + Nprod for n = 1,…,N, there is however no easy way of ensuring that the inequality constraints (35) are fulfilled. To overcome this difficulty, we will instead use a different approach, where we define a set of (Ninj + Nprod − 2)·N functions n

n

0 ≤ ξ1 ; N ; ξNinj

+ Nprod −2

≤ 1; for n = 1…N;

ð38Þ

inj + Nprod − 2 and denote ξn = {ξni }N . Now the vn are given in terms of ξn by i=1 the relations,

n

n

v1 = ξ1  n n n v2 = 1−ξ1 ξ2  n n  n n v3 = 1−ξ1 1−ξ2 ξ3 ⋮ n vNinj −1

=

n ξNinj −1 Ninj −1

vnNinj = ∏



i=1

Ninj −2



∏

i=1

ð39Þ

n 1−ξi

 1−ξni ;

and

3.3. The additional constraints on the control variables vn1

= ξnNinj   vn2 + Ninj = 1−ξnNinj ξn1 + Ninj   vn3 + Ninj = 1−ξnNinj 1−ξn1 +

In addition to the constraints given in Eq. (1), we have other equality constraints given in Eqs. (2) and (3), restated here Ninj

n

∑ vi = 1; for n = 1; …; N;

ð33Þ

i=1

Nprod + Ninj

∑

i = 1 + Ninj

n

vi = 1; for n = 1; …; N:

ð34Þ

 Ninj

ξn2 + Ninj ð40Þ

⋮ vnNprod

and

+ Ninj

+ Ninj −1

= ξnNprod

Nprod + Ninj −3 + Ninj −2

Nprod + Ninj −2

n

vNprod

+ Ninj



∏

=

i=1

∏

i=1



1−ξni



n 1−ξi :

In addition, we also have inequality constraints given in Eq. (4), n

0≤vi ≤1; for i = 1; …; Ninj + Nprod ;

ð35Þ

for n = 1,…,N. In order to deal with the equality constraints (33) and (34) it is possible to add them to the augmented Lagrangian functional Eq. (15). However we want to compare the convergence of the KKT algorithm to that of the adjoint method with the reservoir flow equations as constraints. Adding additional constraints to the augmented Lagrangian functional will further complicate the problem and the convergence will not be clear, if the constraints are handled in this way. Another method of handling the linear equality constraints (33) and (34), could be to do the optimisation with respect to the variables vn1, v2n,…,vNinjn − 1 and vNninj + 1, vNninj + 2,…,vNninj + Nprod − 1 for n = 1,…,N. Then we can calculate the remaining variables using the linear equality constraints Eqs. (33) and (34) so that n

Ninj −1

n

vNinj = 1− ∑ vi for n = 1; …; N i=1

Now we observe that, when defining the rates as in Eqs. (39) and (40), Eqs. (33) and (34) hold for all n

n

0 ≤ ξ1 ; N ; ξNinj

+ Nprod −2

≤ 1:

ð41Þ

Using ξn instead of vn, we have a optimisation problem without the linear equality constraints. However the inequality constraints in Eq. (38) are transformed into the inequality constraints in Eq. (41), which we handle with a projected gradient method. This is automated in the LBFGS code written by Nocedal and details regarding this, can be found in Byrd et al. (1995) and Byrd et al. (1994). To find the derivative, we have by the chain rule that for n = 1,…,N, n

N

inj ∂Lnc ∂Lnc ∂vj n = ∑ n n ; for 0 ≤ i ≤ Ninj −1; ∂ξi j = 1 ∂vj ∂ξi

ð42Þ

ð36Þ and

and n vNinj + Nprod

Nprod + Ninj −1

= 1−

∑

i = 1 + Ninj

n

n vi

for n = 1; …; N:

ð37Þ

∂Lc = ∂ξni

n n ∂Lc ∂vj n n ; for Ninj ≤ i ≤ Ninj + Nprod −2: + Ninj ∂vj ∂ξi

Ninj + Nprod

∑

j=1

ð43Þ

30

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

The derivative given by

∂Lnc , ∂vnj

we already know, and the other derivative is

8 0 > > > > > j−1   > > > ∏ 1−ξnk > > > > n < k=1 ∂vj n = ξj j−1  n > ∂ξni − ∏ 1−ξk > > > 1−ξni k = 1 > > > > > > 1 j−1  > n > :− ∏ 1−ξk 1−ξni k = 1

if i N j if i = j; and j; i b Ninj if i b j; and j; i b Ninj if i b j = Ninj

3.4. The additional constraints on the state variables Our state variables consists of pressures and saturations in all the grid cells. On the saturation there are some additional inequality constraints, stating that the saturations cannot be less than the residual saturation, Swr ≤ Si ≤ 1−Sor :

ð44Þ

When Eq. (22) is performed, in the KKT algorithm, the constraints (44) may be violated. If so is the case, the saturations are set to either Swr or 1 − Sor, depending on which constraint that is violated. That is, we do as follows 1. Find ûnk {pnk , Snk } from (22) 2. For i = 1,2,…,Mc do: 3. If Snki b Swr n

Ski = Swr : 4. If Snki N 1 − Sor n

Ski = Sor : Note that this procedure does not guarantee mass balance. However, we have done tests where the forward model (for which material balance is required) is run with the final optimal controls obtained using the optimization are used, and the results of such a run has always been close to the results from the KKT algorithm. Thus we do not believe that this is a big problem. 4. Numerical examples In this section we present numerical examples of profit optimisation done with the KKT method. Herein we also present comparisons with the more traditional adjoint method. In all of our examples we use synthetic reservoir models inspired by the models of Brouwer and Jansen (Brouwer and Jansen, 2004) with dimensions 450 m × 450 m × 10 m. Along the left hand side of the reservoir there is one injection well, with one segment in each grid cell along the well. We have the possibility to control the injection at each of the individual segments at each time step. Similarly, there is one production well along the right hand side of the reservoir, with one segment in each grid cell along the producer. For the producer we also have the possibility to control the individual production rates at each of the segments at each time step.

Table 2 Reservoir simulator constants—equal for all experiments. κro

κrw

eo

ew

Sor

Swr

μo

μw

ϕ

1.0

0.1

2.0

1.5

0.2

0.2

0.5

0.5

0.2

Table 3 The numerical experiments. Exp.

Type

Grid

Sim. time

Time steps

# Control var.

c

1 2 3 4 5 6

One channel One channel Two channels Two channels One channel Two channels

5×5 15 × 15 15 × 15 15 × 15 45 × 45 45 × 45

2 years 2 years 2 years 2 years 5 years 5 years

64 64 64 64 171 149

640 1920 1920 1920 15390 13410

1·108 1.25·107 1·107 1·107 4·106 5·106

As an initial guess for our rates, we use uniform injection and production, such that the individual injection and production rates are equal at all the grid cells which are penetrated by a well. For the KKT method we also need an initial guess for the state variables u. We find this initial u by solving the forward problem for the initial rates, and using the resulting state variables, u, as our initial guess for u. Initially the water saturation is equal to the residual water saturation of 0.2 everywhere in the reservoir, and the initial pressure is 190 bar, i.e. 19·106 Pa, in the entire reservoir. Other reservoir simulator constants are given in Table 2. For the revenue of oil produced, we set Io = 80 $/m3 and set the cost associated with water production to Iw = − 20 $/m3. Although the revenue constant is low comparing with the prices of today, it is set deliberately to this value since we wish to make a comparison with the study done by (Brouwer and Jansen, 2004). The annual interest rate is set to 20% in Eq. (10). Six experiments were performed as listed in Table 3. The permeability fields applied are plotted in Fig. 2. The first experiment is a small, simple test case with 25 grid cells. The reservoir has one high permeable channel. In this case, one can see intuitively that uniform injection and production over the smart wells will result in an early water breakthrough since the fluids will flow faster through the channel. This serves well as an initial test, since the optimisation algorithms should be able to avoid such a situation. The next 3 experiments are coarse versions of example 1, 2, and 3 in Brouwer and Jansen (Brouwer and Jansen, 2004), while in the last two we apply the same model size of 45 × 45 cells. In all examples, a uniform injection and production is expected to give early water break through and reduced profit, but we expect that it will be more difficult for the optimizer to find the optimum control in the cases with two channels. The main difference between the permeability fields in experiments 3 and 6 compared to the one for experiment 4 lies in the permeability contrast between the channels and the surroundings, which should give different optimal valve settings. One pore volume (i.e., true pore volume, not moveable pore volume) of water is injected in all cases, but the total rate is less in the last two experiments giving a longer total simulation time. In all our experiments, we have used a constant penalty parameter, c, throughout the iterations. It is possible to change this parameter from iteration to iteration. This may give faster convergence, but also introduces other uncertainties and makes it necessary to find a robust method to change the parameter in way that ensures fast convergence to the optimum of the objective function and at the same time honouring the constraints. The value of the penalty constant c in the augmented Lagrangian functional Eq. (15), for the six experiments are given in Table 3. The constant is found experimentally for each specific experiment on a trial and error basis, where we start with an initial value of c, and then monitor the development of the constraints en. If T the value of c is too small, we typically observe that the value of enen increases rapidly (after about five iterations) towards a very high value. When this happens, we multiply c by ten and try again. This process is repeated until we see a decrease of the constraints. Then we divide the penalty parameter by 2 until we start to experience the problem of non-convergence again. The final c-value is chosen as the smallest one which give a decreasing residual after about five iterations. The initial value is chosen arbitrarily as c = 106. This way

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

of determining c increases the amount of work by a negligible amount compared to the work involved in solving the optimal control problem. In addition, we see from Table 3 that the c-value for experiments 2, 3 and 4 is approximately the same, indicating that the same value for the parameter can be used for problems with the same number of grid cells. It is our experience that if one uses a too high value of the constant c, the KKT algorithm converges slower. On the other hand, if one uses a too small value for c, the constraints may not converge to zero such that the algorithm does not converge. To compare the efficiency of the KKT method with the adjoint method, we also plot the objective function against the number of times that we solve N linear systems of equations of size 2·Mc.

31

4.1. Results The behaviour of the KKT algorithm for all the experiments is shown in Figs. 3–12. Comparison with the adjoint method is shown for the first four examples. We expect that a comparison with the adjoint method would be similar for the last two experiments. In Figs. 3, 4, 7, 8, 11 and 12 we show the development of the objective function and the residual as a function of iterations of the KKT algorithm. Since the equation constraint is only fulfilled at convergence for the KKT method, it is not straightforward to compare the convergence behaviour of the KKT method with the adjoint method. In order to facilitate this comparison, we solve the forward

Fig. 2. Permeability fields, given in m2.

32

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

Fig. 3. Experiment 1: Objective function and equation residual for the 20 first iterations of the KKT algorithm. We see that the objective function, J, stabilizes when the residual converge.

Fig. 4. Development of the NPV for experiment 1. The NPV for the adjoint method is plotted as a dashed line and the true NPV, J(vm,ũm) for the KKT method is shown as a solid line. We have $ on the vertical axis and the number of times N linear systems of equations of size 2·Mc are solved on the horizontal axis. We see that the KKT method finds maximum after 20 units on the horizontal axis, while it takes 60 for the adjoint method, showing that the KKT method is faster than the adjoint.

Fig. 5. Experiment 2: Objective function and equation residual for the 30 first iterations of the KKT algorithm. We see that the objective function, J, stabilizes when the residual converge.

Fig. 6. Development of the NPV for experiment 2. The NPV for the adjoint method is plotted as a dashed line and the true NPV, J(vm,um) for the KKT method is shown as a solid line. We have $ on the vertical axis and the number of times N linear systems of equations of size 2·Mc are solved on the horizontal axis. We see that the KKT method finds maximum after 500 units on the horizontal axis, while it takes 800 for the adjoint method, showing that the KKT method is faster than the adjoint.

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

33

Fig. 7. Experiment 3: Objective function and equation residual for the 30 first iterations of the KKT algorithm. We see that the objective function, J, stabilizes when the residual converge.

Fig. 8. Development of the NPV for experiment 3. The NPV for the adjoint method is plotted as a dashed line and the true NPV, J(vm,ũm) for the KKT method is shown as a solid line. We have $ on the vertical axis and the number of times N linear systems of equations of size 2·Mc are solved on the horizontal axis. We see that the KKT method finds maximum after 100 units on the horizontal axis, while it takes 200 for the adjoint method, showing that the KKT method is faster than the adjoint.

Fig. 9. Experiment 4: Objective function and equation residual for the 60 first iterations of the KKT algorithm. We see that the objective function, J, stabilizes when the residual converge.

Fig. 10. Development of the NPV for experiment 4. The NPV for the adjoint method is plotted as a dashed line and the true NPV, J(vm,ũm) for the KKT method is shown as a solid line. We have $ on the vertical axis and the number of times N linear systems of equations of size 2·Mc are solved on the horizontal axis. We see that the KKT method finds maximum after 70 units on the horizontal axis, while it takes 150 for the adjoint method, showing that the KKT method is faster than the adjoint.

34

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

Fig. 11. Experiment 5: Objective function and equation residual for the 200 first iterations of the KKT algorithm. We see that the objective function, J, stabilizes when the residual converge.

Fig. 12. Experiment 6: Objective function and equation residual for the 200 first iterations of the KKT algorithm. We see that the objective function, J, stabilizes when the residual converge.

Fig. 13. Final estimated relative rate allocation for experiment 1. Injection well (top) and production well (bottom). Horizontal axis: discrete time steps. Vertical axis: Well segments.

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

35

Fig. 14. Final estimated relative rate allocation for experiment 2. Injection well (top) and production well (bottom). Horizontal axis: discrete time steps. Vertical axis: Well segments.

Fig. 15. Final estimated relative rate allocation for experiment 3. Injection well (top) and production well (bottom). Horizontal axis: discrete time steps. Vertical axis: Well segments.

Fig. 16. Final oil saturation distributions for experiment 3, showing the KKT algorithm to the left and the adjoint algorithm to the right.

36

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

Fig. 17. Difference in final oil saturation distributions for experiment 3, where the square root of the squared difference is shown.

problem using the latest set of control variables after each iteration of the KKT method. The resulting state variables given the rates vk will be denoted ũk. That is, ũk = {u|en(vnk , un, un + 1) = 0, for n = 1…N}. Correspondingly, Jk̃ = (vk, ũk) as opposed to Jk = J(vk, uk), which is the objective function calculated from a current iteration of the KKT method. For the different experiments, we plot Jk = J(vk, uk), and J̃k versus the number of iterationsqofffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the KKT algorithm. We also provide plots of the development of

∑Nn = 1 ‖en ðvk ; uk Þ‖2 = ð2⋅Mc NÞ to show the

convergence of the residual for the various experiments. Note that because the initial guess is taken from a forward model run, with an equation residual given by the convergence criteria of the forward model, the equation residual will always start at a low value. During the iterations, the equation constraint is not fulfilled, so it increases initially, before it converges to a low value as the KKT iterations proceed. Note also that the equation residual converges faster than the NPV in all cases. When the equation residual approaches zero, and it seems that these become Jk = J(vk, uk) approaches Jk̃ as expected, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi approximately equal at a value of ∑Nn = 1 ‖en ðvk ; uk Þ‖2 = ð2⋅Mc NÞ approximately equal to 10− 3 in all cases. For comparison, the convergence criteria used in the forward model corresponds to a value of qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ∑Nn = 1 ‖en ðvk ; uk Þ‖2 = ð2⋅Mc NÞ equal to 7·10− 4.

Fig. 19. Final estimated relative rate allocation for experiment 5. Injection well (top) and production well (bottom). Horizontal axis: discrete time steps. Vertical axis: Well segments.

We also see from the Figs. 5, 6, 9 and 10 that the optimal profit found with the KKT method is approximately the same as found with the adjoint method. However, the computational effort is less for the KKT method since it takes less computational effort to find a similar NPV value, as we see in for example Fig. 5 where the KKT method finds maximum after 20 units on the horizontal axis, while it takes 60 for the adjoint method. Notice that the NPV does not change initially for the adjoint method. This is due to the initial line search of the LBFGS optimisation routine used by the adjoint method. Although not clearly seen because of the scale, there will be several linear solves without any increase in the adjoint NPV for every iteration, corresponding to the average number of iterations applied in the outer (Newton) iteration loop in each time step of the forward solver. In these examples there is a large period with little changes at the end of the simulations, where only 1 or 2 Newton iterations are required. Thus, the gain in computational speed with the KKT method may be even larger in other cases. The final allocated relative rates found with the KKT method are shown in Figs. 13–20. Comparison with rates obtained with the adjoint method are shown for the first 4 experiments. The injection rates are shown in the top sub-figures, and the production rates are shown in bottom sub-figures. Time steps are on the horizontal axis, and well

Fig. 18. Final estimated relative rate allocation for experiment 4. Injection well (top) and production well (bottom). Horizontal axis: discrete time steps.Vertical axis: Well segments.

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

37

Let us consider the simplified problem min Jðv̂ ; ûÞ = min v̂;û

v̂;û

1 T 1 T T T T v̂ Hv̂ v̂v̂ + v̂ Hû û û + v̂ Hû v̂ û + gv̂ v̂ + gû û ; 2 2 ð45Þ

subject to Kðv̂ ; ûÞ = Nv̂ + Bû −b = 0;

ð46Þ

where Hv̂v̂∈Rnc × nc , Hûû∈Rns × ns , B∈Rns × ns , Hûv̂∈Rnc × ns , N∈Rnc × ns , gv̂∈Rnc , gû∈Rns , b∈Rns and K∈Rns and all of these vectors and matrices are independent of û and v̂. Furthermore, B is assumed to be non-singular. Let us define the augmented Lagrangian as Fig. 20. Final estimated relative rate allocation for experiment 6. Injection well (top) and production well (bottom). Horizontal axis: discrete time steps. Vertical axis: Well segments.

c T Kðv̂ ; ûÞ Kðv̂ ; û Þ; 2

T

Lc ðv̂ ; ûÞ = Jðv̂ ; ûÞ + Kðv̂ ; ûÞ η +

ð47Þ

where η∈Rns is the Lagrange multiplier. Now, the KKT conditions of optimality for Lc gives that segments are on the vertical axis. Note that the plots end earlier than what is indicated in Table 3. This is because it is not profitable to continue production. From these figures we see that the optimal rate distribution for the two methods are similar, but not the same, suggesting that there are several rate distributions giving the same optimal profit, corresponding to several maxima for the objective function. In Fig. 16 we have plotted the final oil saturation distributions obtained with the KKT method and the adjoint method for the two-channel case in example 3. The square root of the squared difference between the two final oil saturation distributions is shown in Fig. 17. Note that this difference is of order 10− 3, which is quite small.

5. Conclusions We have presented a method for solving a NPV maximisation problem based on solving the Karush Kuhn Tucker equations for the augmented Lagrangian functional. In the examples tested, the KKT method finds approximately the same NPV value as the adjoint method with less computational effort. The performance of the KKT method is dependent on the value of the penalty parameter, c, but in the examples tested we have observed that the penalty parameter is approximately the same for reservoirs with the same number of grid cells. This suggest that, once it has been found, we do not need to change it further, when the permeability field is changed in a closed-loop profit optimisation process.

Acknowledgements This work is financed by the Norwegian Research Council, Petromaks programme, project No. 163383/S30, and we are grateful for the support. Furthermore, we want to thank Raymond Martinsen for providing the forward simulator code.

2

T

T

Hv̂ v̂ + cN N 6 6 H T + cBT N 4 û v̂

Hû v̂ + cN B T

Hû û + cB B

N

B

N

T

T

B 0

3 3 2 T T v̂ gv̂ −cN b 76 7 7 6 76 û 7 = −6 g T −cBT b 7: 54 5 5 4 û η −b 32

ð48Þ

The second order conditions of optimality give that " z

T

Hv̂ v̂

Hû v̂

# z N 0;

HûT v̂ Hû û

ð49Þ

Solving the problem (Eqs. (45)–(46)) with the KKT algorithm gives the following iterations −1

ûk + 1 = −B ηk +

1

v̂ k + 1

ð−b + Nv̂ k Þ       = −B−T gûT −cBT b + HûT v̂ + cBT N v̂ k + Hû û + cBT B û k + 1  −1 h   = − Hv̂ v̂ + cN T N gv̂T −cNT b + Hû v̂ + cNT B û k + 1 + NT η k +

1

:

ð50Þ We may simplify this iteration by substituting ûk + 1 in the second equation by the first, giving

ηk +

1

= −B

−T

T

T

T

T

ðgû −cB b + ðHû v̂ + cB NÞv̂ k T

−1

−ðHû û + cB BÞB

ð51Þ

ð−b + Nv̂ k ÞÞ:

Furthermore we substitute ûk + 1 and ηk + 1 in the third equation by the first and the second equation, giving

Appendix A. The quadratic problem with linear constraints In order to understand the convergence properties of the optimisation method proposed, it makes sense to investigate the performance of the algorithm on a simplified problem. And we will thus consider a quadratic optimal control problem with linear constraints.

v̂ k

+ 1

 −1 T   T T T −1 = − Hv̂ v̂ + cN N gv̂−cN b− Hû v̂ + cN B B ð−b + Nv̂ k Þ   T −T T T T T ð52Þ gû −cB b + Hû v̂ + cB N v̂ −N B k   T −1 − Hû û + cB B B ð−b + Nv̂ k Þ :

½

ð

Þ

38

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39

Proof. Denoting an eigenvalue of M3 by α and the corresponding eigenvector by w, we have that

By simple algebraic manipulation this algorithm becomes û k +

= −B

1

ηk + 1 = −B v̂ k

−1

−T

ð−b + Nv̂ k Þ

ðgû

= ðHv̂ v̂ +

+ 1

T

T −T

+N B

ð53Þ −1

+ Hû û B

T

b + ðHû v̂ −ðHû û ÞB

T T −1 cN NÞ ½−gv̂ −Hû v̂ −1

Hû û B

T −T

T

−1

B

−1

b + ðHû v̂ B −1

+ cN N−N B

−1

Hû û B

T −T

b+N B T −T

N+N B

gû

NÞv̂ k Þ

T

ðHû v̂ B−1 N + NT B−T Hû v̂ + cNT N−NT B−T Hû û B−1 NÞw

ð65Þ

T

= αðHv̂ v̂ + cN NÞw:

T

Multiplying Eq. (65) on the right by wT, and can reformulate as

T

Hû v̂

wT Hv̂ v̂ w−2wT Hû v̂ B−1 Nw + wT NT B−T Hû û B−1 Nw

NÞv̂ k :

= ð1−αÞðwT Hv̂ v̂ wT+ c‖Nw‖2 Þ:

Now we define 2

3 ûk ω̂ k = 4 ηk 5; v̂ k

ð66Þ

By choosing ð54Þ

2

3 0 0 M1 M = 4 0 0 M2 5 0 0 M3

 z=

w

−1

B

w

 ;

ð67Þ

in the second order optimality conditions (49), we find that ð55Þ

T

−1

T

w Hv̂ v̂ w−2w Hû v̂ B

T

T −T

Nw + w N B

−1

Hû û B

Nw N 0;

ð68Þ

so that Eq. (66) gives that and T

3 β1 β = 4 β2 5; β3

ð56Þ

where −1

M1 = −B

−T

M2 = −B

−1

T

ðHû v̂ −ðHû û ÞB T

−1

−1

ðHû v̂ B

ð58Þ

NÞ; T −T

N+N B

T

T

T −T

Hû v̂ + cN N−N B

−1

Hû û B

−T

β2 = −B

−1

T

ðgû + Hû û B T

β3 = ðHv̂ v̂ + cN NÞ

−1

ð61Þ

bÞ;

T

ð−gv̂ −Hû v̂ B

−1

T −T T gû +

b+N B

T −T

N B

−1

Hû û B

T

2

w Hv̂ v̂ w + c‖Nw‖ N 0;

ð70Þ

and it follows that α b 1:

ð60Þ

b;

ð71Þ

Now, we look at the second part. Let ε N 0, and assume by contradiction that α ≤ −ε. Suppose, that for a sequence of c → ∞, there exists a sequence of corresponding unit eigenvectors wc→w ≠ 0 so that Eq. (66) holds wTc Hv̂ v̂wc −2wTc Hû v̂ B−1 Nwc + wTc NT B−T Hû û B−1 Nwc = ð1−αÞðwTc Hv̂ v̂ wTc + c‖Nwc ‖2 Þ:

bÞ:

ð62Þ

T

T

2

≥ð1 + εÞðwc Hv̂ v̂ wc + c‖Nwc ‖ Þ

2

ð63Þ

This algorithm will be convergent if the spectral radius of the matrix M, ρ(M)|max(eig(M)| is less than 1. From the definition of M in Eq. (55) and the iteration (63), we observe that ûk + 1 and ηk + 1 depend on v̂k and that they are independent of ûk and ηk, yielding that ρ(M) = ρ(M3). Thus, the algorithm is convergent if the spectral radius of M3 ρðM3 Þ = j maxðeigðM3 Þj b 1:

Lemma 1. The spectral radius of M has the following properties 1. ρ(M3) b 1 2. For ∀ε N 0, ρ(M3) ∈ (−ε,1) for c large enough.

ð72Þ

ð73Þ

If we divide Eq. (72) by c, and let c tend to its limit, we get that

Then we can write the system Eq. (53) as ω̂ k + 1 = Mω̂ k + β:

ð69Þ

NÞ;

ð59Þ −1

2

According to Finsler's lemma, see for example Bonnans et al. (Bonnans et al., 2006), p. 285 , the second order optimality conditions (49) gives that T

ð57Þ

N;

M3 = ðHv̂ v̂ + cN NÞ

β1 = B

T

ð1−αÞðw Hv̂ v̂ w + c‖Nw‖ Þ N 0:

2

ð64Þ

ð1 + εÞ‖Nw‖ ≤ 0;

ð74Þ

Hence it follows that Nw = 0, and Eq. (72) give when c tends to its limit T

T

T

w Hv̂ v̂ w≥ð1 + εÞðw Hv̂ v̂ w Þ;

ð75Þ

which yields the desired contradiction. By Lemma 1, we can always choose a sufficiently large c such that ρ(M3)∈(− 1,1) and thereby verifying condition (64) such we are ensured convergence. References Asheim, H., 1988. Maximization of water sweep efficiency by controlling production and injection rates. Paper SPE 18365 presented at the European Petroleum Conference, 16–19 October 1988, London, United Kingdom. Aziz, K., Settari, A., 1979. Petroleum Reservoir Simulation. Kluwer Academic Publishers. Bonnans, Frédéric, Charles, Gilbert J., Claude, Lemaréchal J., Sagastizábal, ClaudiaA., 2006. Numerical Optimization, Theoretical and Practical Aspects. Springer Verlag.

D.C. Doublet et al. / Journal of Petroleum Science and Engineering 69 (2009) 25–39 Brouwer, D.R., Jansen, J.-D., 2004. Dynamic optimization of waterflooding with smart wells using optimal control theory. SPE J. 9 (4), 391–402. Byrd, Richard H., Nocedal, Jorge, Schnabel, Robert B., 1994. Representations of quasinewton matrices and their use in limited memory methods. Math. Program. 63 (4), 129–156. Byrd, Richard H., Lu, Peihuang, Nocedal, Jorge, Zhu, Ciyou, 1995. A limited memory algorithm for bound constrained optimization. SIAM J. Sci. Stat. Comput. 16 (5), 1190–1208. Doublet, D.C., Martinsen, R., Aanonsen, S.I., Tai, X.-C., 2006. Efficient optimisation of production from smart wells based on the augmented lagrangian method. Presented at European Conferance on the Mathematics of Oil Recovery X, Amsterdam, The Netherlands. September. Doublet, D.C., Aanonsen, S.I., Tai, X.-C., 2007. Efficient history matching and production optimisation with the augmented lagrangian method. Presented at SPE Reservoir Simulations Symposium, 26-28 Feburary 2007, Houston, Texas, USA, 2007. Golub, G.H., Van Loan, Charles F., 1996. Matrix Computations. The Johns Hopkins University Press. Hestenes, M.R., 1969. Multiplier and gradient methods. J. Optim. Theory Appl. 4 (1), 303–320. Ito, Kazufumi, Kunisch, Karl, 1990. The augmented lagrangian method for parameter estimation in elliptic systems. SIAM J. Control Optim. 28 (1), 113–136.

39

Kunisch, K., Tai, X.-C., 1997. Sequential and parallel spitting methods for bilinear control problems in hilbert spaces. SIAM J. Numer. Anal. 34 (1), 91–118. Lien, M., Brouwer, R., Jansen, J.D., Mannseth, T., 2006. Multiscale regularization of flooding optimization for smart field management. Paper 99728-MS, in proceedings of the Intelligent Energy Conference and Exhibition, 11–13 April, Amsterdam, The Netherlands. Nocedal, Jorge, Wright, Stephen J., 1999. Numerical Optimization. Springer Verlag. Ramirez, W.F., 1987. Optimal control theory to enhance oil recovery, developments in petroleum science. Elsevier, Amsterdam. Saad, Yousef, 2000. Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics, Philadelphia, PA. Sarma, Pallav, Aziz, Khalid, Durlofsky, Louis J., 2005. Implementation of adjoint solution for optimal control of smart wells. Paper 92864-MS, in proceedings of the SPE Reservoir Simulation Symposium, 31 January–2 Feburary, The Woodlands, Texas. Sudaryanto, B., Yortsos, Y.C., 2001. Optimization of displacements in porous media using rate control. Paper SPE 71509 presented at the SPE Annual Technical Conference and Exhibition, 30 September–3 October 2001, New Orleans, Louisiana. Zakirov, Iskander, S., Aanonsen, Sigurd Ivar, Zakirov, Ernest, S., Palatnik, Boris, M., 1996. Optimizing reservoir performance by allocation of well rates. The European conference on the Mathematics of Oil Recovery, Leoben, Austrian.