A44
Ensemble Based Multi-Objective Production
Optimization of Smart Wells
R.M. Fonseca* (Delft University of Technology), O. Leeuwenburgh (TNO) & J.D. Jansen (Delft University of Technology)
SUMMARY
In a recent study two hierarchical multi-objective methods were suggested to include short-term targets in life-cycle production optimization. However this previous study has two limitations: 1) the adjoint formulation is used to obtain gradient information, requiring simulator source code access and an extensive implementation effort, and 2) one of the two proposed methods relies on the Hessian matrix which is obtained by a computationally expensive method. In order to overcome the first of these
limitations, we used ensemble-based optimization (EnOpt). EnOpt does not require source code access and is relatively easy to implement. To address the second limitation, we used the BFGS algorithm to obtain an approximation of the Hessian matrix. We performed experiments in which a water flood was optimized in a geologically realistic multi-layer sector model. The controls were inflow control valve settings at pre-defined time intervals. Undiscounted Net Present Value (NPV) and highly discounted NPV were the long-term and short-long-term objective functions used. We obtained an increase of approximately 14% in the secondary objective for a decrease of only 0.2-0.5% in the primary objective. The study demonstrates that ensemble-based multi-objective optimization can achieve results of practical value in a computationally efficient manner.
ECMOR XIII – 13th European Conference on the Mathematics of Oil Recovery
Introduction
Optimization of hydrocarbon production has traditionally been a mainly reactive process. Closed-loop reservoir management aims to convert this practice into a more proactive process by model based optimization. There are numerous methods that are applied to achieve model-based optimization of hydrocarbon recovery. Some of them use gradient-based techniques, where the gradient information is obtained with the aid of an adjoint formulation see e.g. Brouwer and Jansen (2004); Sarma et al. (2005, 2006); Jansen (2011). The adjoint approach is computationally very efficient but has the disadvantage that it requires access to the simulation code to be implemented. Chen et al. (2008) introduced an ensemble-based optimization method (EnOpt), which is computationally less attractive but does not require simulator access and has proven to achieve good results. We use this optimization technique for this work. Many recent model based optimization studies are focused on optimizing over the life of the reservoir, as is appropriate, for example for water flooding strategies that typically aim to prevent early water breakthrough at wells. However, operational decisions are generally based on short-term objectives of a project in terms of operational criteria, production contracts etc., and strategies to optimize such objectives are often be in conflict with optimal long-term strategies. Jansen et al. (2009) observed that significantly different optimized long-term water flooding strategies result in nearly equal values of the objective function, defined as net present value (NPV). They concluded that the life cycle optimization problem is ill-posed and contains redundant degrees of freedom (DOF’s). Thus, there exist multiple solutions to the optimization problem, and different initial starting points may lead to different solutions in an optimal subset of the decision variable space. This formed the basis of the multi-objective optimization approach of Van Essen et al. (2011). They suggested a hierarchical optimization structure to include secondary objectives into the life cycle optimization using the adjoint formulation. They observed a significant increase in short-term objectives with minimal change to the primary objective function. In this paper we investigate the applicability of the EnOpt instead of the adjoint method for multi-objective optimization. We also propose a modification for the hierarchical optimization which improves the computational efficiency of the algorithm. This papers aims to provide a practical and relatively easy to implement alternative to adjoint based multi-objective optimization. We first provide an overview of the theoretical aspects followed by the application to a 3D reservoir model.
Theory
Life Cycle Production Optimization. Life cycle optimization has been successfully applied to
various hydrocarbon recovery mechanisms. This requires at least one decision variable as well as a model that provides relatively accurate long-term predictions. The most widely implemented secondary recovery mechanism in the petroleum industry is water flooding. Water flooding is our choice of recovery mechanism for the following reasons.
• There are many decision variables involved in a water flooding strategy.
• The process is well understood and can be modeled accurately over long time intervals. Recent improvements in technology have led to an increase in the application of downhole chokes or inflow control valves to regulate flow rates and maintain pressure in the reservoir. Smart wells, i.e. wells with Inflow Control Valves (ICVs), are important tools to increase oil recovery and delay water production. ICVs are also used to achieve a uniform water front thus sweeping the reservoir more efficiently (Brouwer and Jansen, 2004). Changes in ICV settings can be used to adjust the inflow from individual zones, reducing the water or gas production, and increase the profitability of a reservoir. Thus it is important to find the optimum value for the settings of the ICVs so as to maximize their potential. In order to improve economic life cycle performance dynamic optimization has to be performed over the producing life of the reservoir due to the transient nature of the saturation
distribution. We consider a control problem where ICV settings can be manipulated to achieve the best possible objective function. The objective function J, the Net Present Value (NPV) is defined as
{
, , ,}
1 ( ) ( ) ( ) , (1 )k t K o k o wp k wp wi k wi k t k q r q r q r t bJ
τ = ⎛ ⎡ ⋅ − ⋅ ⎤ ⎡− ⋅ ⎤ ⋅ Δ ⎞ ⎣ ⎦ ⎣ ⎦ ⎜ ⎟ ⎜ ⎟ ⎜ + ⎟ ⎝ ⎠=
∑
(1)where q is the oil production rate in bbl/day, o k, qwp k, is the water production rate in bbl/day, qwi k, is the water injection rate in bbl/day, ro is the price of oil produced in $/bbl, r is the price of water wp produced in $/bbl, r is the price of water injected in $/bbl, wi Δ is the difference between consecutive tk time steps in days, b is the discount factor expressed as a fraction per year, tk is the cumulative time in days corresponding to time step k, and τt is the reference time for discounting (365 days).
Ensemble Optimization (EnOpt). This technique, first introduced by Chen et al. (2008), is a
stochastic gradient-based optimization method, which utilizes an ensemble of control vectors to estimate a gradient. This method approximates the gradient based on the sensitivity of the ensemble averaged over the objective function with respect to the controls. Distinct characteristics of the EnOpt method are (Chen et al., 2008; Chen, 2008):
• The search direction (gradient) is obtained from the ensemble.
• It can be applied to maximize the expected objective function based on multiple geological realizations.
• It is largely independent of simulator specifics, and requires minimal code development. • It has been shown to work with high-dimensional control vectors.
Approximating the gradient from the sensitivity of the ensemble enables the use of any type of control variables without modification to the existing algorithm. A detailed description of the mathematics for EnOpt is available in Chen et al. (2008). The following is a short overview of the EnOpt method. EnOpt uses an ensemble of controls, for a given reservoir model, to approximate the gradients of the objective function J with respect to the controls. The dynamic well control variables form a control vector 1 1 ( ,...,T T T) ( ,..., )T K u uN = = u u u , (2)
where N is the number of controls which can be rather large. Thus u is a ‘super vector’ with a number of elements N that may be as large as the number of time steps K times the number of control variables (ICVs). For a single geological realization, the resulting stochastic gradient with respect to the controls given by EnOpt is
1 1 1 ( )( ) , 1 M i i u i J J M g C− u u = = − − −
∑
(3) where 1 1 M i i J J M = =∑
, 1 ( ) 1 u M ≈ − T C U U ,and where M is the number of ensemble members of controls, Ji is the objective function for each ensemble member, U is a n by M matrix of the difference between the individual ensemble elements and their mean. Using the gradient g approximated from equation (3) in combination with a steepest ascent scheme where l is the iteration index, we update the set of controls according to
1
l+ =αl l + l
ECMOR XIII – 13th European Conference on the Mathematics of Oil Recovery
Thus an updated set of controls is obtained with αl as a step length along the gradient. The process is iterated until an optimum is achieved.
Multi-Objective Optimization The process of optimizing systematically and simultaneously a
collection of objective functions is called multi-objective optimization. There are various methods of multi-objective optimization such as weighted sum method, Pareto optimality, goal programming etc. Van Essen et al. (2011) introduced a hierarchical optimization procedure to solve the multi-objective production optimization problem. They proposed two methods with the adjoint formulation which are explained below. We have, in this work used this hierarchical approach to investigate its applicability with the EnOpt and exploit the redundant degrees of freedom in the problem.
Hierarchical Optimization. Van Essen et al. (2011) introduced a hierarchical optimization structure
to achieve multi-objective production optimization, which prioritizes the objective functions. A general formulation of this hierarchical optimization is described below. The optimization of the secondary objective function J2 is constrained by a maximum allowable change in the primary objective function. Thus the primary objective function J1 will remain close to its optimal value. A general formulation for hierarchical optimization is as follows
1: 2 1: 1 1 1 1 1 * 1 1 1: max ( ) , . . ( , , ) , 0,...., 1, ( , ) , 0,..., 1 , ( ) , K K k k k k k k k K J s t k K k K J J ε + + + + + = = − ≤ = − − ≤ u u f u x x 0 c u x 0 u (5)
where u is the control vector (input vector, ICV settings), x is the state vector (grid block pressures, saturations), f is a vector valued function that represents the system equations, x0 is the state vector representing the initial state of the reservoir, the subscript k indicates discrete times and K is the total number of time steps. The vector of inequality constraints c related to the system for e.g. capacity limitations etc. The parameter ε> has an appropriately small value compared to J0 1*. Solving the above equations requires the knowledge of J1* which is the optimized value of J1 obtained from the primary objective optimization. Thus the hierarchical optimization constrains the optimization of the secondary objective with respect to the primary objective function. The ordering of the different objective functions is not unique, thus secondary objectives can be implemented as primary objectives and vice versa. This hierarchical structure is attractive when there is a presence of redundant degrees of freedom in the primary objective function. To exploit these degrees of freedom we require the Hessian of the primary objective function. Some concepts detailing the need for this Hessian and methods to exploit the redundant degrees of freedom are explained in the appendix.
Approximate Hessian Van Essen et al. (2011) proposed the use of a finite difference scheme in
combination with the adjoint formulation to approximate the second order derivatives of the objective function. Without the adjoint formulation, to estimate a finite difference based Hessian we require
n×(n+1) function evaluations where n is the number of controls. Thus this method is computationally
infeasible for realistic reservoir models and large numbers of controls, due to the high number of function evaluations needed. To alleviate this short-coming we propose to use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm which approximates and updates the Hessian during the optimization of the primary objective function. This leads to a significant reduction in overall computational costs incurred during optimization of the secondary objective function.
BFGS Algorithm From an initial guess u0 and an approximate Hessian matrix, H0 the following steps are repeated until u converges.
1. Obtain a direction pk by solving: H pk k = −∇J( )u . k
2. Determine an acceptable step size αk in the direction found in the first step, and then update 1
k+ = k +αk k
3. Set sk =αkp . k
4. Compute yk = ∇J(uk+1)− ∇J( )u . k
5. Compute the updated Hessian as 1 ( )
T T T k k k k k k k k T T k k k k k + = + − y y H s s H H H y s s H s .
J(u) denotes the objective function to be minimized. Practically, H0 can be initialized with H0 = I, so that the first step will be equivalent to a gradient descent, but further steps are more and more refined by Hk, the approximation to the Hessian.
Hierarchical Optimization Algorithm The algorithm is a modification of the algorithm proposed by
(Van Essen et al., 2011). The modification is the implementation of the BFGS algorithm to approximate the Hessian matrix.
1. Find a single optimal strategy u* that maximizes the primary objective function J
1 and use u = u* with n = 0 as a starting point for the secondary optimization problem where n is the iteration index.
2. Use the approximated Hessian H at u* and perform a singular value decomposition to obtain the orthonormal basis B for the null-space of H. Project P onto the null-space according to
T
P = BB .
3. Find the gradient s for the secondary objective function J2.
4. Project this improving direction s onto the orthonormal basis B to obtain the projected direction d, such that d is an improving direction for J2 and does not affect J1. Thus d is
⋅
d = P s.
5. Update the control vector un using the projected direction d in the steepest-ascent method.
1
n+ = n+ ⋅τ
u u d ,
where τ is an appropriately small step size.
6. Update the Hessian H using the BFGS algorithm for the new set of controls. 7. Perform steps 2 to 6 until convergence of J2
The EnOpt algorithm has been used to approximate the gradient of the secondary objective function s. Further details on the algorithm sketched above are provided in the Appendix.
Switching Algorithm The hierarchical algorithm presented above is computationally cumbersome
and not feasible for realistic reservoir models having a large number of input parameters when implemented using a finite difference based Hessian. To overcome this short-coming, Van Essen et al. (2011) presented a practical alternative method to the hierarchical algorithm with the use of a switching function according to
1 1 2 2 bal
J = Ω ⋅ + Ω ⋅ , J J (6)
where Ω and 1 Ω are switching functions for J2 1 and J* that take on values of 1 and 0 or vice versa:
* 1 1 1 1 * 1 1 * 1 1 2 1 * 1 1 1 if ( ) 0 if 0 if ( ) 1 if J J J J J J J J J J ε ε ε ε ⎧ − > Ω = ⎨ − ≤ ⎩ ⎧ − > Ω = ⎨ − ≤ ⎩ (7)
Here ε is the threshold value as defined in the inequality constraint in eq. (5). J1* is the value of the primary objective at the optimal solution achieved during life cycle optimization. The gradient of Jbal with respect to the input parameters is then given by
ECMOR XIII – 13th European Conference on the Mathematics of Oil Recovery 1 2 1 1, 2 1, 1 1 1 ( ) ( ) . bal n n n n n dJ dJ dJ J J du + = Ω ⋅ du + + Ω ⋅ du + (8)
The use of a balanced objective function in the optimization will give improving directions for either
J1 or J2, thus switching between feasible and infeasible solutions. However, the convergence towards an optimal solution may be rather slow due to the switching between the different solutions. In order to improve convergence speed Van Essen et al. (2010) suggested the following adaptation, in which gradients of the secondary objective function are projected onto the null space of the optimal primary objective function. Thus applying the general formula of projection we have
1 1 1 1 1 T T dJ dJ dJ dJ d d d d − ⎛ ⎞ = ⎜⎜ ⎟⎠ ⎝ P u u u u . (9)
In the neighbourhood of the optimum, the complement of the gradient with respect to the primary objective function can be used as a first order approximation to the null space of the Hessian of this function. Thus the improved direction for the secondary objective is given by
2 [ ] T dJ d = ⋅ a I - P u . (10)
Hence the alternative switching search direction d for solving the hierarchical problem is
1 2 1 1 2 1 1 ( ) ( ) [ ] T T T n n n n n dJ dJ J J d d + + + = Ω ⋅ + Ω ⋅ ⋅ d I - P u u . (11)
The EnOpt is used to approximate the gradients for both the primary and the secondary objective functions. The two hierarchical methods presented above are tested on a geologically realistic sector model explained below.
Numerical example
3D Synthetic Reservoir Model A 3D synthetic reservoir model was used to test the multi-objective
optimization methods. The life cycle of the reservoir covers 15 years, or 5470 days. The model illustrated in figure 1 consists of 25 x 32 x 5 = 4000 grid blocks. The approximate size of the grid blocks is 110 x 90 x 20 meters. The model represents an area of 2.5 km by 3.5 km by 100 m. The geological structure consists of connected uplifted/offset blocks. There is a sealing fault on the North-Western side of the block, close to producer 1. The initial average reservoir pressure is 200 bars. Table 1 lists the geological and fluid properties used to describe the model. The reservoir is produced using an inverted 5 spot well pattern, i.e. 4 producers at the edges of the grid with an injector in the center of the grid. The reservoir is divided into 5 layers having different horizontal permeabilities, however the permeability in each individual layer is constant. The wells penetrate all the five layers with an ICV in every layer resulting in a total of 25 controls per time step. The producing life of the reservoir is divided into 15 time intervals of one year (365 days) each, which results in a total of 375 controls that are to be optimized. Water is injected at a constant pressure of 300 bars and the production wells are operated at a minimum of 15 bars. A Corey model with exponents equal to 2 for both oil and water is used for the relative permeability’s where the connate water saturation is 20%, the residual oil saturation is 30% and the end point relative permeabilities to oil and water are 0.8 and 0.4 respectively. No capillary pressures have been included in this model. The reservoir rock is incompressible. The simulator used in this study is ECLIPSE 100.
Figure 1 Illustration of the reservoir model used to test multi-objective optimization methods.
Table 1 Geological and flow properties of the synthetic reservoir model.
Property Values Units
Porosity 20% --
Permeability- (layer 1 – layer 5) 100-300-50-600-100 mD
Reservoir Pressure 200 bar at 1950 m
Density of oil 800 kg/m3
Density of water 1000 kg/m3
Temperature 77 0F
Oil compressibility 4e-5 1/bar
Water compressibility 4e-5 1/bar
Viscosity of oil 2 cP at 1 bar
Viscosity of water 0.5 cP at 1 bar
Life Cycle Optimization An optimal life-cycle strategy of ICV settings for the individual layers is
obtained by optimizing the NPV as described in equation 2, with ro = 130 USD/m3, rwp = 25 USD/m3,
rwi = 6 USD/m3. The discount rate b was set to 0 Figure 2 is an illustration of the optimization with undiscounted NPV as the objective function, which is equivalent to cumulative cash flow over the producing life of the reservoir. The optimal solution u* was obtained using a steepest ascent scheme with the line search method and an ensemble size of 100 samples. Well productivity index (PI) multipliers are used to model ICVs in ECLIPSE with bounds of 1e-4 and 1. The starting point for the optimization is an initial control vector having values equal to 1. Thus all the ICVs/chokes are open as
ECMOR XIII – 13th European Conference on the Mathematics of Oil Recovery
a starting strategy. The optimization was allowed to run for 80 iterations although there was no significant improvement in objective function value after 65 iterations as indicated in Figure 2. Additional iterations were performed to allow the BFGS algorithm to estimate a Hessian matrix that is as close to the true Hessian at the optimum as possible. The optimized value of the objective function is USD 8.902x109 $.
Figure 2 Life cycle optimization with undiscounted NPV as the objective function.
Objective Functions We use undiscounted NPV as the primary objective function. Every project
aspires to recover the initial investments as soon as possible. In an ideal case we would like to additionally meet such short-term economic objectives whilst still maintaining the life cycle objectives. Thus we choose a secondary objective function which highlights the importance of maximizing short-term production. The secondary objective function has the same cost structure as the primary objective function but with a very high discount rate b of 25%.
Unconstrained Optimization First, the secondary objective function is optimized without being
constrained by the primary objective function. This case serves as a comparison to the hierarchical structure explained in the theory. The optimization was performed with an ensemble size of 75 members and was allowed to run for 65 iterations. The results are illustrated in Figure 3. The starting point of the optimization is the optimal solution achieved during optimization of the primary objective. A decrease of 1.65 % is seen in the primary objective function to achieve an increase of 14.2% in the secondary objective function.
10 20 30 40 50 60 70 80 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 9x 10 9 Undiscounted NPV (USD) Number of Iteration
Primary Objective Function
Figure 3 Illustration of an unconstrained optimization of secondary objective (red) starting from the optimal strategy u* obtained from life cycle optimization.
Hierarchical Optimization Figure 4 illustrates the optimization of the secondary objective function
within the null space of the primary objective function. Since none of the singular values of the Hessian are typically exactly equal to zero, a cut-off criterion must be defined to estimate the null space of the Hessian matrix. We have used a cut-off criterion of σi /σ1< 2e-9, where σi are the singular values of the Hessian matrix, with the values arranged from largest to smallest. The resulting null space consists of 187 vectors and its dimension remains almost constant throughout the optimization. The value of the secondary objective function at the optimum of the primary objective is 3.35·109 $,which is the starting point for the hierarchical optimization. We achieve a value of 3.823·109 $ after completing the hierarchical optimization equivalent to a 14.2% increase in the secondary objective function at the price of a corresponding 0.52% decrease in the primary objective function. This decrease is much smaller than that obtained with the unconstrained optimization approach and clearly illustrates the advantage of using the multi-objective optimization approach. Another illustration of the impact of the hierarchical optimization structure is provided in Figure 5 which shows a comparison of the cash flow (i.e. undiscounted NPV) over time resulting from the optimal life-cycle strategy (green), the strategy resulting from hierarchical optimization (red) and from unconstrained optimization of the secondary objective (blue). An increase in the short to medium term cash flow is observed with the strategy obtained by hierarchical optimization when compared to the life-cycle optimization strategy.
20 40 60 8 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 9x 10 9 Undiscounted NPV (USD) Number of Iteration
Primary Objective Function
20 40 60 3.35 3.4 3.45 3.5 3.55 3.6 3.65 3.7 3.75 3.8 3.85 3.9x 10 9 25% Discounted NPV(USD) Number of Iteration
Secondary Objective Function
14.2% increase in secondary objective function
1.6% decrease in primary objective function
ECMOR XIII – 13th European Conference on the Mathematics of Oil Recovery
Figure 4 Comparison of the primary (blue) and secondary (red) objective functions obtained by the hierarchical optimization algorithm.
The multi-objective optimization is seen to be useful to increase cash flow in the initial stages of the project whilst maintaining the life cycle goals. The unconstrained optimization as expected has the best cash flow in the short term period of the economic life. However, as shown in the inset figure, compromises are made to the life-cycle targets. The inset plot in Figure 5 also shows that the hierarchical optimization achieves a solution which performs better in the long-term compared to the unconstrained optimization.
Switching Algorithm The mathematical formulation of the switching algorithm described in the
theory requires the definition of a criterion ε. The criterion used is ε = 0.003 J*. An ensemble size of 75 samples was used and the optimization was allowed to run for 50 iterations; however after 35 iterations the improvements in the secondary objective were minimal. Results are shown in figure 6. The algorithm achieves an increase of 14.17% in the secondary objective for a corresponding 0.21% decrease in the primary objective. The switching algorithm thus performs very favourably when applied to this model. Since the performance is dependent on parameters such as ensemble size used for gradient evaluation, step length, size of the trust region etc., fine-tuning of these parameters may lead to further improved performance of the switching algorithm.
10 20 30 40 50 8 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 9x 10 9 Undiscounted NPV (USD) Number of Iteration
Primary Objective Function
10 20 30 40 50 3.35 3.4 3.45 3.5 3.55 3.6 3.65 3.7 3.75 3.8 3.85 3.9x 10 9 25% Discounted NPV(USD) Number of Iteration
Secondary Objective Function
0.5% decrease in primary objective function 14.2 % increase in secondary
Figure 5 Plot showing the different cash flows over time for the different optimization strategies: life-cycle optimized strategy (green), hierarchical optimized strategy (red) and unconstrained secondary optimized strategy (blue).
Similar to Figure 5 for the hierarchical optimization, Figure 7 compares the cash flow over time for the strategy from switching algorithm with the unconstrained and life-cycle only optimization results. It is observed that after 500 days the cash flow with life cycle optimization is approximately 1 ·109 USD. However the control strategy obtained with switching optimization achieves a cash flow of 1.5 ·109 USD. This increase of 0.5·109 $ over 500 days which will enable the project to achieve the break-even point faster, thus making this strategy extremely attractive. Similar to the results obtained in Van Essen et al (2011) and the hierarchical optimization results shown above, the NPV at the end time of the unconstrained optimization (blue curve) is decreased by 1.6% compared to the NPV obtained by the switching algorithm (red curve) which decreases by only 0.2%. Finally, in Figure 8 the set of controls obtained by the optimization using the switching algorithm (red) is compared to the optimal set of control obtained after life cycle optimization (black). The control sets are fairly different, so it must be concluded that rather different control sets may achieve very similar results for the primary optimization while drastically improving the secondary optimization.
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 1 2 3 4 5 6 7 8 9x 10 9
Cash Flow over Time (USD)
Time(days)
Comparison of Primary Objective Cash Flow for different Control Strategies
Value of objective function J
1 for initial optimal input J * Value of objective function J
1 for hierarchical optimization of J2 Value of objective function J
1 for unconstrained optimization of J2
3500 4000 4500 5000 7.8 8 8.2 8.4 8.6 8.8 9x 10 9 Times (days)
Cash Flow over Time (USD)
ECMOR XIII – 13th European Conference on the Mathematics of Oil Recovery
Figure 6 Illustration of switching algorithm optimization with projection applied to the gradient showing a 14.3% increase in secondary objective function (red) for a corresponding 0.2% decrease in primary objective (blue).
Discussion
The results for the model used here show that the switching method achieves better results compared to the hierarchical method in that the same improvement in secondary objective was obtained with fewer iterations and with a smaller decrease in the primary objective. It is not possible at this point, however, to draw any definitive general conclusion regarding the comparative performance of these methods for more complex models. The choice for the different objective functions may have a large impact and significantly affect the scope for multi-objective optimization. The Hessian approximation with the BFGS algorithm has shown to achieve good results for this case but may perform differently in a different case. The cut-off criterion used to define the dimension of the null space was chosen somewhat arbitrarily in this work. However this criterion was found to be very important for the success of the hierarchical optimization (not shown).
5 10 15 20 25 30 35 8 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 9x 10 9 Undiscounted NPV (USD) Number of Iteration Primary Objective Function
5 10 15 20 25 30 35 3.35 3.4 3.45 3.5 3.55 3.6 3.65 3.7 3.75 3.8 3.85 3.9x 10 9 25% Discounted NPV(USD) Number of Iteration Secondary Objective Function
14.2% increase in secondary objective function
0.21% decrease in primary objective function
Figure 7 Comparison of the cash flow over time for the switching (red), unconstrained (blue) and optimal life-cycle (green) strategies.
Results Summary
The two multi-objective methods using the EnOpt method for gradient evaluation have shown an improvement in a secondary objective constrained to a maximum specified decrease in the primary objective. The results obtained with the hierarchical optimization and switching algorithms are very similar. The results also demonstrate the attractiveness of constrained optimization relative to unconstrained optimization. Table 2 highlights the comparison between the different methods used in this paper to achieve multi-objective optimization.
Table 2 Summary of objective function vales and changes for the two multi-objective optimization
algorithms. Optimization Method Primary Objective (109) Secondary Objective (106) % Increase in second obj. % Decrease in primary obj. Hierarchical 8.854 3.827 14.23 0.52 Switching 8.882 3.823 14.17 0.21 Unconstrained 8.727 3.825 14.2 1.65 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 1 2 3 4 5 6 7 8 9x 10 9
Cash Flow over Time (USD)
Time(days)
Comparison of Primary Objective Cash Flow for different Control Strategies
Value of objective function J
1 for initial optimal input J *
Value of objective function J
1 for switching optimization of J2
Value of objective function J
1 for unconstrained optimization of J2
3500 4000 4500 5000 7.8 8 8.2 8.4 8.6 8.8 9x 10 9 Time (days)
Cash Flow over Time (USD)
ECMOR XIII – 13th European Conference on the Mathematics of Oil Recovery
Figure 8 Comparison of the control strategy for the switching algorithm (red) with optimal control strategy obtained during life cycle optimization (black), for the individual ICVs.
Reactive Control
For the given oil price of 130 $/bbl and water production costs of 25 $/bbl, the economic feasibility threshold will be achieved at a water cut (WC) of 83%. Thus a reactive control strategy is defined with a well shut-off limit of 83% water cut (WC). The wells are operated with all the ICVs fully open until the WC limit has been reached, after which the wells are shut-in. As expected the reactive control strategy has a much lower NPV at the end of the life of the reservoir. Figure 9 is a comparison of the cash flow for the reactive control strategy to the life-cycle and the switching based optimized strategies. The strategy obtained with the switching method (red) achieves an improved short-term performance compared to optimized life-cycle strategy (green line), but it is not as good as the reactive control strategy (blue), which gives the best short-term performance. Thus the two multi-objective optimization methods presented improve the short-term/secondary multi-objectives but do not truly recover the best possible short-term/secondary objectives as obtained with the reactive strategy. However, Figure 9 does confirm the advantage of optimal life-cycle strategies in comparison to a reactive control based operational strategy when long-term objectives are important.
2737.5 5475 0.5 1 ICV1 Days INJ 2737.5 5475 0.5 1 ICV2 Days 2737.5 5475 0.5 1 ICV3 Days 2737.5 5475 0.5 1 ICV4 Days 2737.5 5475 0.5 1 ICV5 Days 2737.5 5475 0.5 1 ICV1 Days PROD1 2737.5 5475 0.5 1 ICV2 Days 2737.5 5475 0.5 1 ICV3 Days 2737.5 5475 0.5 1 ICV4 Days 2737.5 5475 0.5 1 ICV5 Days 2737.5 5475 0.5 1 ICV1 Days PROD2 2737.5 5475 0.5 1 ICV2 Days 2737.5 5475 0.5 1 ICV3 Days 2737.5 5475 0.5 1 ICV4 Days 2737.5 5475 0.5 1 ICV5 Days 2737.5 5475 0.5 1 ICV1 Days PROD3 2737.5 5475 0.5 1 ICV2 Days 2737.5 5475 0.5 1 ICV3 Days 2737.5 5475 0.5 1 ICV4 Days 2737.5 5475 0.5 1 ICV5 Days 2737.5 5475 0.5 1 ICV1 Days PROD4 2737.5 5475 0.5 1 ICV2 Days 2737.5 5475 0.5 1 ICV3 Days 2737.5 5475 0.5 1 ICV4 Days 2737.5 5475 0.5 1 ICV5 Days
Figure 9 Comparison of the cash flow over time for the switching (red), reactive control (blue) and optimal life-cycle (green) strategies.
Conclusions
• Compromises made to short-term targets during life cycle optimization can be partly corrected for with an ensemble based multi-objective optimization method.
• The EnOpt is a good alternative to achieve practical results when the adjoint formulation is not available for multi-objective optimization.
• In our numerical simulation examples, two constrained multi-objective methods showed a 14.2% improvement in the secondary objective function (NPV @ 25% discount rate) constrained to the primary objective function (NPV @ 0% discount rate). The results obtained with the hierarchical optimization algorithm are similar to those resulting from the switching algorithm, although for the cases investigated here, the switching algorithm was found to perform slightly better.
• The BFGS algorithm used to estimate the Hessian for the hierarchical method, is computationally attractive compared to a finite difference method especially when dealing with large control sets, and led to good results for the case reported here.
• A reactive control operational strategy achieves the best possible short-term performance, at the cost of a decreased long-term performance.
• Multi-objective optimization of ICV settings shows significant scope for improvement in short to medium term goals constrained to life cycle targets.
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 1 2 3 4 5 6 7 8 9x 10 9
Cash Flow over Time $
Time(days)
Comparison of Primary Objective Cash Flow for different Control Strategies
Value of objective function J
1 for initial optimal input J * Value of objective function J
1 for switching optimization of J2 Value of objective function J
1 for reactive control 83% WCT
3500 4000 4500 5000 7.4 7.6 7.8 8 8.2 8.4 8.6 8.8 9x 10 9 Time (days)
Cash Flow (USD)
ECMOR XIII – 13th European Conference on the Mathematics of Oil Recovery
Acknowledgements
This research was carried out within the context of the ISAPP Knowledge Centre. ISAPP (Integrated Systems Approach to Petroleum Production) is a joint project of TNO, Delft University of Technology, ENI, Statoil and Petrobras. The authors would also like to acknowledge Paul Egberts (TNO) for his advice and contributions during this work.
References
Brouwer, D.R. and Jansen, J.D. [2004] Dynamic Optimization of Water Flooding With Smart Wells Using Optimal Control Theory. SPE Journal, 9(4), 391-402. DOI: 10.2118/78278-PA.
Chen, Y., Oliver, D.S. and Zhang, D. [2008] Efficient ensemble-based closed-loop production optimization. Paper SPE 112873 presented at the SPE/DOE Improved Oil Recovery Symposium, Tulsa, USA, 19-23 April. DOI: 10.2118/112873-MS.
Chen, Y. [2008] Efficient ensemble based reservoir management. PhD thesis, University of Oklahoma, USA.
Jansen, J.D., Douma, S.G., Brouwer, D.R., Van den Hof, P.M.J., Bosgra, O.H. and Heemink, A.W. [2009] Closed-loop reservoir management. Paper SPE 119098 presented at 2009 SPE Reservoir Simulation Symposium, The Woodlands, Texas, U.S.A, 2-4 February. DOI: 10.2118/119098-MS.
Jansen, J.D. [2011] Adjoint-based optimization of multiphase flow through porous media – a review. Computers and Fluids, 46(1), 40-51. DOI: 10.1016/j.compfluid.2010.09.039.
Oliver, D.S., Reynolds, A.C. and Liu, N. [2008] Inverse theory for petroleum reservoir characterization and history matching. Cambridge University Press, Cambridge.
Sarma, P., Aziz, K. and Durlofsky, L.J. [2005] Implementation of Adjoint Solution for Optimal Control of Smart Wells. Paper SPE 92864 presented at the SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA, 31 January–2 February. DOI: 10.2118/92864-MS. Sarma, P., Chen, W.H., Durlofsky, L.J. and Aziz, K. [2006] Production Optimization With Adjoint
Models Under Nonlinear Control-State Path Inequality Constraints. SPE Reservoir Engineering and Evaluation, 11(2), 326-339. DOI: 10.2118/99959-PA.
Van Essen, G.M., Van den Hof, P.M.J and Jansen, J.D. [2011] Hierarchical Long and Short-Term Production Optimization. SPE Journal, 16(1), 191-199. DOI: 10.2118/124332-PA.
Van Essen, G.M., Van den Hof, P.M.J and Jansen, J.D. [2010] Determination of lower and upper bounds of predicted production from history-matched models. Proc. 12th European Conference
on the Mathematics of Oil Recovery ECMOR XII, Oxford UK, 6-9 September 2010. Appendix
Degrees of freedom. When dealing with an optimization problem, there may exist multiple sets of
control variables for which we achieve similar results. This set of different optimal control variables is an indication of the presence of redundant degrees of freedom (DOF) in the system. This was first noticed by Jansen et al. (2009) who showed that different solutions exist for the optimal control problem of maximizing the economic objective function over the producing life of the reservoir. A similar observation in a history matching setting was made by Oliver et al. (page 145). This existence of multiple solutions suggests, when the optimality of a primary objective function is reached not all DOFs of the control variable space are fixed. This implies that there may exist redundant degrees of freedom in the optimization problem. This conclusion formed the basis for development of the multi-objective optimization algorithm in Van Essen et al. (2011).
Null Space Van Essen et al. (2011) approximate a Hessian matrix to find these DOFs, which is an
Hessian matrix, a matrix of second order derivatives of the objective function. Consider an objective
function J and let u* be a control vector. IfΔu is a vector (with same length as u*) of small perturbations then a Taylor expansion around the vector u* is given by
3 2 1 ( ) ( ) ( ) ( ) ( ) ( ) . 2 T T J u*+ Δ ≈u J u* + ∇J u* ⋅ Δ +u Δu ⋅∇ J u* Δ +u O Δu (A-1) If u* is a (local) optimum of J and u* is in the interior (i.e. not on the boundary) of the feasible domain
for a constrained optimization then we can conclude that ( )
J
∇ u* =0 (A-2)
Substituting equation (A-2) into equation (A-1) we get
3 2 1 ( ) ( ) ( ) ( ) ( ) . 2 T J u*+ Δ ≈u J u* + Δu ⋅∇ J u* Δ +u O Δu (A-3) If we choose 2 [ ( )] Null J Δ ∈ ∇ * u u then 2 ( ) J ∇ * Δ =
u u 0 . Thus equation (A-4) reduces to
3
( ) ( ) ( )
J *+ Δ ≈J * +O Δu
u u u . (A-4)
Thus equation (A-4) implies that for a small perturbation 2 [ ( )]
Null J
Δ ∈ ∇ *
u u , the adjusted control vector (u*+ Δu will have an objective function value very close to the objective function value) J(u*) which is an optimal value. This means that we can make an update to the control vector that is in the null space of the primary objective function to improve the secondary objective function. This proves the need to find the Hessian matrix at the optimum of the primary objective function and the set of vectors that span its null space. A singular value decomposition is used to obtain the null space and orthonormal basis B used in the algorithm.
