Generation of a Pareto front for a bi-objective water flooding optimization problem using approximate ensemble gradients

Conflicting objectives are frequently encountered in most real-world problems. When dealing with conflicting objectives, decision makers prefer to obtain a range of possible optimal solutions from which to choose. In theory, methods exists that can produce a range of possible solutions, some of which are “Pareto Optimal”. The application of these methods to solve bi-objective production optimization problems is increasing. A recent paper introduced a method to find points on the boundary of the objective function space by solving a constrained optimization problem using adjoint gradients. In this work, we investigate the applicability of using ensemble optimization (EnOpt) (which relies on approximate ensemble gradients instead of exact adjoint-based gradients) to generate points along a “Pareto” front with acceptable computational effort.. Moreover, we investigate the applicability of this approximate gradient technique to solve constrained optimization problems using the augmented Lagrangian method. Finally, we compare the performance of this bi-objective optimization method to a traditional weighted sum method for bi-objective water flooding optimization of two different synthetic reservoir models. The two objectives used in this work are, undiscounted (0%) net present value (NPV), representing long-term targets and highly discounted (25%) NPV, representing short-term operational targets. The controls are inflow control valve (ICV) settings over time for one model and water injection rate controls for the other. The effect of different starting points and the computational efficiency of the constrained optimization method are also investigated. & 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
A majority of studies and applications of life-cycle water flooding optimization using a model-based approach have focused on a single objective optimization with emphasis being placed on the theoretical understanding and practical application of the optimization methodology. Life-cycle optimization essentially aims to find a strategy which optimizes long-term reservoir management targets, but life-cycle optimization is often at the expense of operationally significant short-term targets. Thus, there is a need to solve a bi-objective problem to obtain a strategy that accounts for the two objectives because the long-term perspective is usually in conflict with the short-term targets which are decided by operational constraints, contractual obligations etc. Van Essen et al. (2011) introduced a hierarchical optimization framework to solve such a multi-objective optimization problem. This was motivated by the observation in, e.g., Jansen et al. (2009) that the objective function space consists of many redundant degrees of freedom which can be exploited to optimize a secondary objective. This hierarchical structure provides a single optimal strategy which incorporates multiple objectives. However, decision makers usually prefer to have multiple strategies to choose from, especially when dealing with conflicting objectives. Isebor and Durlofsky (2014) applied an evolutionary algorithm to generate points along a "Pareto" front for a bi-objective water flooding problem. The main pitfall of this approach was the computational effort required to obtain the points on a Pareto front. Also they did not compare the front generated with any other method used to generate Pareto fronts to check if the front obtained was Pareto optimal. Liu and Reynolds (2014) applied the normal boundary intersection method (NBI) first introduced in Das and Dennis (1998) to a bi-objective water flooding problem with and without geological uncertainty. Liu and Reynolds (2014) showed that the NBI method is computationally more efficient than the method of Isebor and Durlofsky (2014) and produces better solutions than the traditional weighted sum method. The NBI method involves solving a series of constrained optimization sub-problems. In Liu and Reynolds (2014), these constrained optimization problems were solved using an augmented Lagrangian method using an adjoint formulation to compute the gradients. The adjoint formulation, an overview of which can be found in Jansen (2011) and references therein, is a computationally efficient method which requires access to the simulator source code to implement. Most commercial simulators either do not have a fully developed adjoint code or access to the source code is not permissible. This has led to an increase in the application of various approximate gradient based techniques which are computationally less efficient but use the simulator as a black-box, and are more flexible. Do and Reynolds (2013) provided theoretical connections between various existing approximate gradient techniques which use an ensemble of perturbed controls to estimate a gradient. One such approximate gradient technique introduced in Lorentzen et al. (2006) and thereafter in its current form by Chen et al. (2009) is the ensemble optimization (EnOpt) method. Recently many studies have used EnOpt for life-cycle production optimization problems. Fonseca et al. (2014) applied EnOpt to solve a bi-objective optimization problem using the hierarchical structure proposed by Van Essen et al. (2011). Additionally there has been an increase in the number of applications of different evolutionary algorithms to solve either a bi-objective optimization problem, Isebor and Durlofsky (2014) etc., or for history matching applications, as detailed in Liu and Reynolds (2014). In this work we investigate the applicability of the EnOpt technique to generate points along a "Pareto" front with acceptable computational effort. A secondary aim is the application of EnOpt to solve constrained optimization problems using the augmented Lagrangian method. Note that Fonseca et al. (2014) consider hierarchical optimization (using EnOpt), in which case an a-priory choice is made which of the two objectives is most important. Here we consider bi-objective optimization (using EnOpt) based on the Pareto front approach which provides freedom to the decision maker to choose the relative importance of each of the two objectives, as will be explained in more detail below.

Theory
This section investigates the applicability of the use of approximate ensemble gradients to calculate points on a Pareto front for bi-objective production optimization problems.

Objective functions
We first define the objective functions followed by an overview of EnOpt. We apply the usual expression for Net Present Value (NPV) as objective function J: , , / k t where q o k , is the average oil production rate in m 3 /day for time step k, q wp k , is the water production rate in m 3 /day for time step k, q wi k , is the water injection rate in m 3 /day for time step k, r o is the sale price of oil in $/ m 3 , r wp is the cost of water produced in $/ m 3 , r wi is the cost of water injected in $/m 3 , Δt k is the length of the k th time step in days, b is the discount factor, t k is the cumulative time in days corresponding to time step k, and τ t is the reference time period for discounting, typically one year (i.e. 365 days). In this work the two objective functions are: (1), representing the long-term objective ("recovery optimization").

Ensemble optimization (EnOpt)
In this section, we outline the standard formulation of the EnOpt algorithm as proposed by Chen et al. (2009). We take u to be a single control vector containing all the control variables to be optimized. This vector has N components where N is equal to the product of the controllable well parameters (number of well settings like bottom hole pressures, rates or valve settings) and the number of control time steps. Chen et al. (2009) sample the initial mean control vector from a Gaussian distribution while, at later iteration steps, the final control vector of the previous iteration is taken as the mean control. However, the initial controls can also be chosen by the user, as will be done in our experiments. The vector of controls is given by, where the counter i preempts the use of multiple control vectors, and where u i is assumed to be a random vector which has a mean u and covariance matrixC, i.e. u i $ N(u,C). Then an ensemble of M independent samples of N(u,C) are generated as, i.e., each z i is a vector of independent standard random normal deviates, andC 1/2 is any square root ofC. In our examples˜= C L 1/2 , where L is the lower triangular matrix in the Cholesky decomposition ofC. We truncate any element of the ensemble of controls outside of the set of bounds to the bound value. Then, the sample mean is computed as To estimate the gradient, a mean-shifted ensemble matrix is defined as In the present paper, we use as the search direction in a steepest ascent algorithm an approximation to the gradient, rather than the approximation of the smoothed gradient that is used in standard EnOpt. The approximate gradient is where the superscript † indicates the Moore-Penrose pseudo inverse, which is conveniently computed using a singular value decomposition (SVD); see, e.g., Strang (2006). Do and Reynolds (2013) demonstrated that it is akin to what is known as a 'simplex gradient', Conn et al. (2009). They also provided theoretical connections between various ensemble methods such as simultaneous perturbation stochastic approximation (SPSA), simplex gradient and EnOpt. Moreover, they proposed a modification to the gradient formulation which uses the current control vector ℓ u and the corresponding objective function value ℓ J to calculate the control and objective function anomalies ΔU and Δj: where the superscript ℓ is the optimization iteration index. In this work, we have used (Eqs. (9) and 10) in combination with Eq. (8).
We note that many authors also use single and double smoothed versions of Eq. (8) which can be obtained by pre-multiplication (either once or twice) of the gradient estimate by the covariance matrix used to generate the ensemble of controls. For the original derivation of EnOpt, including an extension to use the method under geological uncertainty (as represented by an ensemble of geologically different reservoir models), we refer to Chen et al. (2009). Here we restrict our application to "deterministic EnOpt", i.e. to optimization under the assumption that the geology is known.

Update rules
The approximate gradient g from Eq. (8) can be used in any gradient-based optimization algorithm. In our study we use the simple steepest ascent scheme given by where the superscript ℓ is the iteration index, and α ℓ is a step length in the direction of the approximate gradient. Note that we scaled the gradient by its infinity norm and choose a step length to be 10% of the difference between the maximum and minimum values of the controls. We allowed for a maximum of five backtracking steps, each time reducing the step size with a factor of one half if the objective function J decreases from one iteration to the next. If after the five back-tracking steps we still do not find an increase in J we accept the current control strategy and continue with the optimization until a convergence criteria is satisfied.
Optimal update schemes and their corresponding parameters are typically case-dependent and more sophisticated line search algorithms are sometimes beneficial. However, we chose to use a relatively simple update strategy to facilitate the comparison of the various multi-objective optimization schemes which form the key subject of our paper.

Multi-objective optimization
Most real world problems have multiple objectives that need to be satisfied. Usually these objectives are in conflict with each other, i.e. one must accept decreases in one objective to achieve increases in another objective. The process of optimizing systematically and simultaneously a collection of objective functions is called multi-objective optimization. In theory, there exist many methods to solve a multi-objective problem and recently there has been an increased focus on finding methods to solve multi-objective problems in the reservoir simulation community. These objectives are usually defined as long-term (life-cycle) objectives from a reservoir engineering viewpoint and short-term objectives from a production engineering/operational constraints viewpoint. Van Essen et al. (2011) showed that these two objectives may be in conflict with each other and suggested the use of a hierarchical framework for multi-objective optimization. An alternative to hierarchical bi-objective optimization (in which the primary objective is considered more important than the secondary objective), is regular bi-objective optimization in which there is no predefined preference for one of the objectives. Isebor and Durlofsky (2014), and Liu and Reynolds (2014) have introduced methodologies to generate the 'Pareto front' i.e. a range of possible solutions for a decision maker for a regular bi-objective reservoir optimization problem. Isebor and Durlofsky (2014) presented their methodology using a hybrid evolutionary algorithm, PSO-MADS, and reported results which were obtained with a significant computational effort. Liu and Reynolds (2014) presented a method using adjoint gradients which was shown to be computationally much more efficient. We use, in this work, the methods introduced in Liu and Reynolds (2014) and investigate their applicability in combination with the EnOpt method.
A point is defined as "Pareto optimal" if at that point the value of one objective function cannot be increased unless the value of a second objective function is decreased or, in other words, a control set is Pareto optimal if there does not exists any other control set which achieves better objective function solutions. Liu and Reynolds (2014) provide details of the commonly used theoretical definitions to determine whether points are non-dominated, i.e. Pareto optimal, and lie on a Pareto front.

Weighted sum method
The life-cycle waterflooding problem is inherently a long-term optimization problem as shown in Van Essen et al. (2011) and short-term goals are sacrificed to achieve the optimal long-term targets. A traditional technique to balance two conflicting objectives is the weighted sum method, see Marler and Arora (2004), which aims to optimize a weighted objective function that combines both objectives in a single function according to where J ws is the weighted sum objective function constructed from the long-term and short-term objective functions J 1 and J 2 with w 1 and w 2 as weighting factors. Liu and Reynolds (2014), among others, showed that the biggest drawback of this method in finding solutions on a Pareto curve is that the solutions tend to be concentrated on one part of the curve, i.e., the solutions generated are not evenly distributed along the Pareto front. Another disadvantage is that the weighted sum method cannot obtain points on the concave part of the Pareto front, see, for example, Fig. 1 of Liu and Reynolds (2014).

Adjusted weighted sum method
To overcome the difficulties of the weighted sum method, Liu and Reynolds (2014) proposed an adjusted weighted sum formulation where the weights w 1 and w 2 are now replaced bỹ and in this case maximizingJ ws corresponds to maximizing J 2 . Liu and Reynolds (2014) found that choosing decreasing w 1 from 1 to 0.1 in increments of 0.1, computing the corresponding values ofw 1 andw 2 , and maximizingJ ws for each of thesew 1 ,w 2 values tended to result in points that were well distributed along the Pareto front when maximizingJ ws , whereas Eq. (12) with the same set of w 1 values did not generate a well-distributed Pareto front.

Normal Boundary Intersection (NBI) method
In order to overcome the disadvantages of the weighted sum method, Das and Dennis (1998) proposed a technique, the Normal Boundary Intersection (NBI) method, to find points on the boundary of a feasible set starting from points along the "utopia line" which is defined as the line in the objective function space that connects the optimum solutions for the individual objective functions. The boundary points are then found by optimizing the magnitude of a unit normal to the utopia line in the objective function space. A detailed description of the NBI method can be found in Das and Dennis (1998) and, for petroleum engineering applications, in Liu and Reynolds (2014). The NBI method is motivated by the fact that the Pareto front must coincide with a part of the boundary of the feasible region. The disadvantage of NBI is that boundary points may or may not be Pareto optimal, i.e., may or may not lie on the Pareto front. However once optimal design vectors * u 1 , …., * u n are generated it is easy to check if each point is non-dominated by any other, which must be the case if (J 1 ( * u n ), J 2 ( * u n )) is a point on the Pareto front; see definitions in Reynolds (2014, 2016). The following is a brief description of the method as described in Liu and Reynolds (2014) for bi-objective water flooding optimization problems. For the two objective functions denoted by J 1 and J 2 the NBI procedure is repeated for different points along the utopia line. The general formulation for NBI is given by and where * u 1 and * u 2 are the optimal control strategies obtained for the individual optimizations of J 1 and J 2. The line segment that connects ( * ) j u 1 and ( * ) j u 2 in the objective space as is defined as the utopia line. To solve the equality constrained optimization problem as described in Eq. (16), Liu and Reynolds (2014) used the augmented Lagrangian method. While there exist several techniques to solve constrained optimization problems, we too have applied the augmented Lagrangian method because the main purpose of this work is to investigate the ability of approximate gradient techniques like EnOpt to generate solutions along a Pareto front. A by-product of this work is the demonstration of the applicability of an approximate gradient technique to solve constrained optimization problems. The augmented Lagrangian method (Nocedal and Wright, 2006), used to solve the different NBI sub-problems, is based on the augmented Lagrangian function which is defined by  (2014), because they were using adjoint gradients, calculated the gradient of the Lagrangian function with respect to u in terms of the gradients of objective functions J 1 and J 2 with respect to u. In this work, because we use approximate ensemble gradients, we calculate the gradient of the Lagrangian function directly using Eq. (8). The following is a brief algorithmic description of the NBI method as implemented in our case.

Tracking the Pareto front using NBI
The NBI method as implemented by Liu and Reynolds (2014) choses as a starting point a combination of the optimal control sets * u 1 and * u 2 depending on the weight factors chosen. Due to the nonlinearity of the problem these initial points usually have objective function values that do not lie exactly on the utopia line. The NBI problem does not necessarily require the starting points to be on or close to the utopia line, so we propose to generate points on the Pareto front by starting from a point on the front which has already been obtained with different values of β 1 and β 2 . This is akin to "tracking" a front. In the results section, we discuss the advantages/disadvantages of using this method of generating solutions on a front.

Hierarchical switching method
Van Essen et al. (2011) introduced a hierarchical optimization scheme to achieve multi-objective production optimization, which prioritizes the objective functions. The optimization of the secondary objective function J 2 is constrained by a maximum allowable change in the primary objective function. Thus the primary objective function J 1 will remain close to its optimal value. The ordering of the different objective functions is the prerogative of the user, thus secondary objectives can be implemented as primary objectives and vice versa. This hierarchical scheme is especially attractive in the presence of redundant degrees of freedom in the primary objective function. Van Essen (2011) proposed two different varieties of the hierarchical scheme: one requiring the computation of the Hessian matrix of the objective function with respect to the controls, and one which, in a more pragmatic fashion, alternatingly optimizes the short and long-term objectives while maintaining the first objective function value close to its initial maximum. Based on the results in Fonseca et al. (2014), we use the hierarchical switching method with EnOpt for the optimization where details of the implementation are provided. This method optimizes the objectives alternatingly with the use of a switching function according to where Ω 1 and Ω 2 are switching functions for J 1 and J 2 that take on values of 1 and 0 or vice versa: . 20 Here ε is the threshold value and J 1 * is the value of the primary objective at the optimal solution achieved during life cycle optimization. We will compare the results obtained from hierarchical switching optimization to the other methods presented above. The advantage of using a hierarchical switching method is that a user can decide the maximum allowable decrease in the primary objective value which is practically impossible to know when using the weighted sum method. However, with this hierarchical method, only a single control set is generated which may or may not be Pareto optimal since no other information is available for comparison. Fig. 1 depicts an overview of various optimization methods available to solve multi-objective reservoir optimization problems.
In the current paper we compare the methods indicated with numbers 1-5 with the aid of two numerical examples.

Reservoir model
Advances in technology have led to an increase in the application of inflow control valves (ICVs) to regulate flow rates and maintain pressure in the reservoir. We consider a control problem where ICV settings of injection and production wells in a 3D synthetic reservoir model, from JOA (2007), are manipulated to optimize waterflooding over the producing life of the reservoir. The model, illustrated in Fig. 2(a), consists of 25 Â 32 Â 5 ¼4000 grid blocks. The approximate size of each grid block is 110 Â 90 Â 20 m, so that the reservoir volume is 2.5 Â 3.5 Â 0.1 km 3 . The geological structure consists of uplifted blocks, separated by faults. The reservoir is produced using an inverted five-spot well pattern, i.e. four producers at the corners of the grid with an injector in the center. The reservoir is divided into five layers with different horizontal permeabilities. Fig. 2(b) is the top view of the transmissibility multipliers used for this model and the white cells are grid blocks that are inactive. There is a sealing fault on the North-Western side of the block, close to producer 1. Table 1 lists the reservoir and fluid properties of the model. A Corey model with exponents equal to 2 for both oil and water is used for the relative permeabilities where the connate water saturation is 0.2, the residual oil saturation is 0.3 and the end point relative permeabilities to oil and water are 0.8 and 0.4 respectively. Capillary pressure effects are not included. The wells penetrate all five layers with one ICV in each layer. The producing life of the reservoir is divided into 15 optimization control steps, each of which is one year (365 days) in duration, and there are 25 controls per control step which results in a total of 15 Â 25 ¼375 controls to be optimized. Water is injected at a constant pressure of 300 bars and the production wells are operated at a minimum pressure of 15 bars. We used an oil price r o ¼ 130 $/m 3 , water production costs r wp ¼25 $/m 3 , and water injection costs r wi ¼6 $/m 3 . Well index multipliers were used to model the ICVs in the simulator with bounds of 1 Â 10 À 4 and 1, where the finite lower bound was chosen to avoid numerical problems in the simulator associated with a zero lower bound. For the simulation of the model we used a commercial fully implicit finite difference black oil simulator (Eclipse, 2011).

Results for NBI
In this section we compare the use of the NBI method without tracking (method 3 in Fig. 1) and with tracking (method 4 in Fig. 1). Following Liu and Reynolds (2014) In Eq. (21), ( * ) j u 1 ¼ [9.1060 Â 10 9 , 3.3522 Â 10 9 ] T and ( * ) j u 2 ¼[8.7086 Â 10 9 , 4.4759 Â 10 9 ] T . The optimization is not dependent on the choice of n. The solution of this equation gives n¼ [2.822, 1] T , which is the same for all the different starting points used in this work. Solving multiple NBI sub-problems for different choices of weight combinations, we obtain the solutions shown in Fig. 3. The black circles are obtained for starting points based on the first step of the NBI algorithm presented previously. The objective of this initialization is to obtain a starting point on or close to the utopia line. Due to the non-linearity of the problem the objective function values achieved for this initial guess are never on the utopia line, but always slightly above the line. Using the solutions already obtained we also test the applicability of finding solutions which satisfy the constraints starting from points (control sets) that have previously satisfied the constraints. This is akin to "tracking" points along a front. The red circles in Fig. 3 are the points achieved when the tracking process begins from β 1 ¼0.1 (point A in Fig. 3). We observe that for most of the weight combinations, the tracking procedure achieve solutions that dominate the solutions represented by the black circles. Since there is no preference to choose from which end the tracking begins, we also began the tracking from β 1 ¼0.9 (point B in Fig. 3), to obtain the solutions shown by blue circles in Fig. 3. We observe that in this case the tracking procedure achieves solutions that dominate the solutions from the other two initialization procedures for all the points. Additionally this tracking procedure is computationally more efficient as is discussed later. Thus, different initial guesses for a given value of β 1 can have a significant impact on the solutions achieved with the bi-objective optimization algorithm. Besides the different starting points, all other algorithmic details are exactly the same for the three different sets of points generated. The gradients are estimated with an ensemble size equal to 30 with a perturbation size equal to 0.001. Table 2 provides the objective function values for 11 different optimum points (black circles) along a boundary front. We observe that for the β 1 ¼0.4 case, we obtain a 0.8% decrease in the primary, long-term objective function from its optimal value (9.1060 Â 10 9 $) and an approximately 22% increase in the secondary, short-term objective function. For the solution obtained with initial guesses based on the front tracking procedure for β 1 ¼0.1 we observe a 0.8% decrease in the primary objective to achieve a 33% increase in the secondary objective. There is only a 5% difference in the primary objective function values between the optimal strategies for the two objective functions J 1 and J 2 , i.e. the first and last points in Table 2. Thus, for the objective functions chosen in this study, we do not expect to observe major increases in primary (long-term) objective for minor decreases in the secondary (short-term) objective, indicating that there may exist fewer redundant degrees of freedom in the short-term objective function.
Independent of the method used to generate the initial guess of a given β 1 , the approximate Pareto front generated with NBI (Fig. 3) shows that one can obtain a sharp increase in the secondary objective function for a very minimal decrease in the primary objective. Moreover, it seems that for this case that the Pareto front consists of two branches; a near-horizontal one near the optimal secondary objective and a near-vertical one near the optimal primary objective. Fig. 4(a) is an illustration of the evolution of the Lagrangian function through the iteration process. The sharp drop in the value of the Lagrangian function corresponds to an update (decrease) of the penalty parameter μ in the augmented Lagrangian method. In most of the cases we observe that we generally perform 5 outer loop iterations in which we update the penalty parameter for 3 iterations and the vector of Lagrange multipliers λ for the remaining two iterations. Fig. 4(b), right-side plot, shows the constraint violation throughout the optimization process. Note that the constraint violation must be less than the given tolerance specified in the optimization algorithm to obtain convergence for the outer loop of the augmented Lagrangian algorithm. When the inner loop  Fig. 3. Boundary points achieved using the NBI method starting from different initial guesses. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) Table 2 Objective function values of the black circles in Fig. 3.
Long-term objective ( Â 10 9 ) $ Short-term objective ( Â 10 9 ) $ converges and the constraint violation is sufficiently small the algorithm converges. Thus it is possible that there are multiple points at which the constraint violation is satisfied, however there is only one point at which both the inner and outer loop's stopping criteria are satisfied. Fig. 5 depicts a comparison of the optimal control settings for two different strategies for the highest permeability layer in producer 2. The blue line is the life-cycle strategy, i.e. one end point of the utopia line, while the red line is the NBI strategy (blue circles in Fig. 3) obtained for weight combinations of β 1 ¼0.1 and β 2 ¼ 0.9, i.e. the strategy that achieved a 33% increase in the short-term objective for a 0.8% decrease in the long-term objective. From Fig. 5 we observe that significantly different strategies can be achieved by performing bi-objective optimization. For the red curve, the ICV setting is almost fully open for the first 10 years with lower setting values towards the end of the producing time period, which is in line with the emphasis on increasing the shortterm increase in NPV. For the optimal life-cycle strategy, the same ICV is almost closed for four or the first five years and then is fully open through most of the remaining producing life, in order to virtually maintain the goal of life-cycle optimization. The optimized control settings for other valves are similar to the trend shown in Fig. 5. Fig. 6 shows the saturation distribution in layer 4 after 4 years of production for the different optimal strategies whose controls are compared in Fig. 5. We see that the optimal life-cycle strategy, being less aggressive, sweeps a much smaller area with less water being injected, while the optimal NBI strategy, i.e. the one for β 1 ¼0.1 using the front tracking procedure, is more aggressive, i.e more is water injected and more oil is displaced and produced.

Comparison of weighted sum techniques
In this section we compare the adjusted sum method (method 1 in Fig. 1) and the adjusted weighted sum method (method 2 in Fig. 1). Liu and Reynolds (2014) showed cases where the adjusted weighted sum method produces a significantly better spread of solutions compared to the traditional weighted sum technique. Tables 3 and 4 provide the solutions for the various weight combinations used where we observe, as reported in Liu and Reynolds (2014), that the adjusted weighted sum technique provides a better spread of solutions, and in particular gives a better representation of the front near the optimal long-term NPV which is the most important part of the front. Fig. 7 provides a visual comparison of the solutions obtained with the two different methods. Note: The stopping criterion used to achieve this set of points is exactly the same as the stopping criterion used for the inner loop in the augmented Lagrangianbased NBI method albeit the objective functions are different.  the NBI method. The results here are very interesting: for w 1 ¼0.9 the solutions obtained with either method do not dominate each other while for w 1 ¼ 0.8 and w 1 ¼0.7 we observe that the adjusted weighted sum method achieves solutions which slightly dominate the solutions obtained with NBI. However, for the other weight combinations the solutions obtained with NBI dominate. It is difficult to know why this behavior is observed and it could be either case dependent or gradient quality dependent. However, Liu and Reynolds (2016) also find that the NBI method generally gives a better representation of the front than is obtained with the weighted sum method. Fig. 9 is a comparison of the optimization path for the different methods with the weight combination w 1 ¼0.7 and w 2 ¼ 0.3. The original NBI and the weighted sum have the same starting point, however they have very different paths. The adjusting of the weights in the adjusted weighted sum method leads to a significantly different starting point for the optimization. All the optimization results shown here are influenced not only by the gradient quality, but for the NBI method, also by the choices of the initial penalty parameter μ and Lagrange multipliers λ. Using a larger ensemble size for the gradient estimate (Fonseca et al., 2015) could lead to smoother optimization paths to solve the individual sub-problems and possibly better solution points, however for computational reasons this has not been investigated. Fig. 10 is an illustration of the total number of simulations performed for each of the methods including the two initial optimization runs to obtain the utopia line. The original NBI method was computationally the most expensive one with approximately 28,000 total simulations, while for the NBI tracking method about half the number of simulations required for the original NBI method were needed to achieve better solutions. Both the weighted sum variants were computationally much more efficient, similar to the results obtained in Liu and Reynolds (2014). We note that ensemble methods to generate approximate gradients, such Table 3 Solutions for different weight combination using the weighted sum method.

Comparison of weighted sum and NBI
w 1 w 2 Long-term objective ( Â 10 9 ) $ Short-term objective ( Â 10 9 ) $  Table 4 Solutions for different weight combination using the adjusted weighted sum method.
w 1 w 2 Long-term objective ( Â 10 9 ) $ Short-term objective ( Â 10 9 ) $  as EnOpt, are well suited for embarrassingly parallel processing. Moreover, when introducing geological uncertainty, in the form of different geological realizations, the computational load of ensemble methods becomes relatively less disadvantageous as compared to adjoint-methods; see Fonseca et al. (2015). We also note that the use of ensemble methods to compute approximate gradients involves various user-defined choices such as the ensemble size and the definition of the co-variance matrix (which may be used to enforce temporal correlations between the controls). Moreover a wide variety of ensemble gradient formulations is available; see, e.g., Fonseca et al. (2015) for an overview. In the present study we did not address these aspects, and we refer to Fonseca et al. (2015) for further information on the effects of the various formulations and user-defined parameters.

Hierarchical switching optimization method
In this section we compare the results of the NBI and weighted sum methods with those obtained with hierarchical switching (method 5 in Fig. 1). The switching method optimizes the objectives alternatingly, while staying within a maximum allowable decrease ε in the primary objective. The choice of ε is user dependent. Thus, the user has to a-priori decide the maximum allowable acceptable decrease in the optimal primary objective function value. Fig. 11 plots the optimization path where a maximum decrease of 0.3% in the primary objective is allowed (red curve). We see that we achieve approximately a 10% increase in the secondary objective. The values obtained are similar to using a weight combination of w 1 ¼ 0.7 and w 2 ¼0.3 for either the NBI method or the adjusted weighted sum method. However the solution is a non-dominated point when compared to all solutions obtained with the NBI and adjusted weighted sum methods. If the optimization is repeated for a 1% allowable decrease in the primary objective we observe that we achieve a 20% increase in the secondary objective function (black dotted curve). This solution however is dominated by the solutions from the other two methods. The hierarchical method only provides a single strategy which may or may not be Pareto optimal. A tracking procedure like the one implemented for the NBI method could be used with this method to generate a front. Alternatively we can use the primary objective function as a constraint while optimizing a secondary objective, i.e., a lexicographic approach; for details see Liu and Reynolds (2016).

Reservoir model
The model, illustrated in Fig. 12, which was introduced in Van Essen et al. (2009), represents a channelized depositional system in the form of a discrete permeability field modeled with 60 Â 60 Â 7¼ 25,200 grid cells of which 18,553 cells are active. A detailed description along with fluid and geological properties of the standardized version of this "egg model" is given in Jansen et al. (2014). The model is produced using eight peripheral water injection wells (blue) and four production wells (red) which are completed in all seven layers. No capillary pressures are included, the reservoir rock is assumed to be incompressible, and the liquids slightly compressible. The controls to be optimized are the water injection rates with a maximum allowable injection rate per well fixed at 79.5 m 3 /day and a minimum rate of 0 m 3 /day. The producers are operated at a minimum bottom hole pressure of 395 bars without rate constraints. The producing life of the reservoir is divided into 40 control steps of 90 days each, and therefore the control vector u has N ¼8 Â 40¼ 320 elements. The NPV parameters used are r o ¼ 126 $/m 3 , r wp ¼ 19 $/m 3 , and r wi ¼5 $/m 3 . The two objectives used are the same as in the previous example. We use a commercial fully implicit black oil simulator (Eclipse, 2011) for the reservoir simulations.  Primary Objective Function (10 9 $) Secondary Objective Function (10 9 $) epsilon = 0.3% Utopia Line epsilon = 1% NBI from 1 = 0.9 Adjusted Weighted Sum

Results and comparison
For this somewhat larger and more complex example we only compare the use of the NBI methods without and with tracking (methods 3 and 4 in Fig. 1). The normal vector n is obtained, similar to the example reported above, by setting the second component of n to 1 and solving Eq. (21) Fig. 14, which satisfy the stopping criteria of the augmented Lagrangian function and the constraint violation. The black circles are obtained for starting points which aim to start on the utopia line. Again, due to the non-linearity of the problem, the objective function values of the starting points are never on the utopia line. The spread in the points found using NBI is more continuous compared to the solutions achieved in the previous example. Again instead of solving the sub-problems from a starting point close to or on the utopia line we aim to "track" points along a front. The blue circles in Fig. 13 are the points achieved when the tracking process begins from β 1 ¼0. The results illustrated in Fig. 13 seem to suggest that there exists different fronts in the objective function space. We observe that till β 1 ¼0.7 the points seem to lie on a line with a certain slope for both the original NBI method as well as the NBI tracking method. From β 1 ¼0.6 onwards the points seem to align themselves along a line with a completely different slope for both the methods. The points obtained with NBI tracking seem to always dominate the solutions obtained with the original NBI method. Thus, similar to the previous example, different starting points of the optimization have significant impact on the solutions achieved. Besides the different starting points, all other algorithmic details are exactly the same for the three different sets of points generated. The gradients are estimated with an ensemble size equal to 30 with a perturbation size equal to 0.01.
With only a 0.7% decrease in the primary (long-term) objective function, we have found a solution that achieves an approximately 19% increase in the secondary (short-term) objective function for the black open circles, i.e. the original NBI method for β 1 ¼0.7. For the same weight combination, the NBI tracking method finds a solution for which a 0.3% decrease in the long-term objective leads to a 20% increase in short-term gains. Additionally, for the blue circles, i.e., NBI with tracking, we observe that for a 1.3% decrease in the primary objective, we can achieve an even more significant increase of 38% in the secondary objective. This last result corresponds to the β 1 ¼0.4 solution.
The results from the previous example illustrated that the adjusted weighted sum method produced a much better spread in solution points compared to the weighted sum method. Thus for this example we compare the NBI solutions with solutions from the adjusted weighted sum technique. We observe, as shown in Fig. 14, that the solutions achieved by the NBI tracking method dominate the solutions from the adjusted weighted sum method. The solutions from the original NBI method also dominate the solutions from the adjusted weighted sum method. A comparison of the plot of remaining oil saturation for the top layer after 3 years of production illustrates the difference in the strategies; see Fig. 15. The optimal long-term strategy is less aggressive, as significantly less area is swept by injected water, while the NBI tracking solution for β 1 ¼0.4, i.e. a 38% increase in short-term gains for a 1.3% decrease in long-term gains, is a more aggressive strategy as larger areas of the reservoir have been swept.
The control strategies that resulted in the saturation plots shown in Fig. 15 are illustrated in Fig. 16. The control settings have only been displayed for the first 3 years to highlight the differences between the strategies. We observe that the optimal NBI strategy injects much more water compared to the optimal longterm strategy. A similar trend in the control strategies is seen in the other wells. Fig. 17 is a comparison of the computational efficiency of the different methods. The original NBI method requires the highest computational effort similar to the results reported for the previous example. We needed about 21,000 simulations to achieve the 11 points for the original NBI method while we needed approximately 14,000 simulations when using the NBI to track the boundary front. The adjusted weighted sum method required less than 8000 simulations to find the 11 points, and is thus the computationally most efficient one, though the solutions achieved are far from Pareto optimal compared to the solutions achieved by the NBI method.

Discussion
The differences between the results of the adjusted weighted sum and the traditional weighted sum methods for this example  are far from significant because the adjustment of the weights does not lead to very different weight combinations, i.e., scaling of the problem is not as important for this problem as it was for the 5 spot ICV problem and the problems investigated in Liu and Reynolds (2014). For this example, the difference between the optimal primary long-term objective function values is 14%, which is much higher than in the previous example, while the difference between the short-term objective function values is 48%.

Conclusions
-Approximate gradient techniques like EnOpt can be used to generate solutions for a bi-objective optimization problem which may lie on a Pareto front although the computational costs are significant. -Tracking the Pareto front using NBI is a computationally more efficient method and produces better solutions for the decision maker to choose from compared to the original NBI method.
Different starting point have a significant impact on the optimal solutions achieved. This is observed for the two different example problems investigated. -The adjusted weighted sum produces a more even distribution of solutions and is marginally computationally more efficient compared to the traditional weighted sum technique for the ICV control problem. -For some weight combinations, the NBI method produces solutions which dominate solutions obtained by the weighted sum variants and vice versa for the ICV control problem. -A hierarchical switching method provides a single solution which satisfies the maximum allowable decrease in the objective function value. For the ICV problem, the single solution point was non-dominated when compared to all the points obtained from the NBI and weighted sum method.

Acknowledgements
We thank the member companies of TUPREP for their support.