A Comparison between Evolutionary Metaheuristics and Mathematical Optimization to Solve the Wells Placement Problem


The Wells Placement Problem (WPP) consists in choosing well locations within an oil reservoir grid to maximize the reservoir total oil production, subject to distance threshold between wells and number of wells cap constraints. A popular approach to WPP is Genetic Algorithms (GA). Alternatively, WPP has been approached in the literature through Mathematical Optimization. Here, we conduct a computational study of both methods and compare their solutions and performance. Our results indicate that, while GA can provide near-optimal solutions to instances of WPP, typically Mathematical Optimization provides better solutions within less computational time.


Share and Cite:

G. AlQahtani, A. Alzahabi, E. Kozyreff, I. Farias and M. Soliman, "A Comparison between Evolutionary Metaheuristics and Mathematical Optimization to Solve the Wells Placement Problem," Advances in Chemical Engineering and Science, Vol. 3 No. 4A, 2013, pp. 30-36. doi: 10.4236/aces.2013.34A1005.

1. Introduction

Consider an oil reservoir for which oil production totals over the entire reservoir area are given (or estimated). We are interested in the Wells Placement Problem (WPP), i.e., the problem of maximizing the total oil production of the reservoir by selecting locations to build at most a certain number of vertical wells at a given minimum distance apart from one another.

WPP (and variations of it) has been the subject of a number of studies in the last decades [1-9]. The main approaches used to find solutions have been Genetic Algorithms (GA), which are a special type of Evolutionary Metaheuristics, and Integer Programming (IP), a subclass of Mathematical Optimization.

To the best of our knowledge, the earliest study of wells placement optimization using IP was the work conducted in 1974 by Rosenwald and Green [8], who developed a numerical optimization framework to select optimal positions of wells. Their objective was to minimize the difference between actual and desired flow curves by choosing optimum well locations from predefined location sets. In 1997, Vasantharajan and Cullick [9] introduced the term reservoir quality to represent reservoir connectivity adjusted by tortuosity, which was the surface area surrounding the flow by the volume contained. The well site selection process was performed through IP, with the goal of maximizing reservoir quality. Later, in 1999, Cullick et al. [4] improved the formulation of the problem by reducing its number of constraints.

On the other side, GA is the most widely used method to solve WPP optimization [1]. Guyaguler et al. [6] developed a hybrid method of GA linked with polytope search and coupled with a proxy generated by kriging. The objective function in this work was to maximize net present value by finding the optimum injection well locations. The output of this work indicates that using kriging as a proxy with the hybrid optimization framework reduced the computational expense considerably. This work was extended by Badru and Kabir [2] by including horizontal wells and linking the hybrid GA with experimental design to find the impact of uncertainty on recovery and the number of wells required. In 2006, Bangerth et al. [3] compared GA with four algorithms designed to solve WPP. GA was also compared with Covariance Matrix Adaptation-Evolution Strategy by Ding [5] to optimize well locations for different types of wells with calculated productivity index as the objective function. Onwunalu et al. [7] used GA combined with statistical proxy to establish a relationship between the variables and the objective function using unsupervised learning techniques.

Some of the articles cited above report computational results using either GA or IP, but we are not aware of a study that compares both. In this paper, we use GA and IP to solve instances of WPP, and report the computational results of those tests.

Specifically, our test reservoir sits on an grid, where each of the sites in the grid represents a possible location for a well. For each, where and, a production total is associated with the site, and each pair of sites chosen needs to be at least grid units apart from each other. Finally, the maximum number of well locations selected, , is predetermined.

2. Evolutionary Metaheuristics

Evolutionary Metaheuristics use ideas of biology and Darwin’s theory of evolution, such as reproduction, mutation, and “survival of the fittest”, to provide solutions to problems that arise in other fields of study. Because one of the most used types of evolutionary metaheuristics in petroleum engineering is genetic algorithms, we developed a GA to tackle our problem of interest (for details on the structure and use of GA in general, see [10]).

2.1. Definition of Terms and Variables

Our GA starts with an initial population that is improved through crossover, local search, mutation, and the addition of new members to the population, which we call intruders. After a generation, the best fit individuals of the population are selected, the others are discarded, and a new generation begins. The algorithm runs for a predetermined number of generations, , and returns the fittest individuals of the final population. The fitness of an individual, here, is the total oil production associated with it.

An individual is represented by an -array of ordered pairs (which we call coordinates), and denoted as

where, , represent the locations on the grid of the wells defined by. Associated with, the value is the total oil production that results from choosing this individual, given by Equation (1).


We denote the production totals matrix by, i.e., an element of row and column of is, and we define We let be the sorted vector of the elements of and its corresponding coordinates:





2.2. The Genetic Algorithm

The high level structure of our GA is as follows:

At each generation, is the number of individuals in the initial population, is the number of crossover operations performed, and is the number of intruders created. The value is used in the sub-routine LOCAL-SEARCH, and will be explained later.

Lines 1 - 5 initialize the variables of the algorithm, which are the individuals of the population, at zero. Line 6 creates the initial population, and lines 7 - 12 perform evolutionary sub-routines for generations. Finally, line 13 returns the best fit individual found.

Throughout the algorithm, and for all individuals considered, it is verified whether the distance of all locations in the grid associated with the individual satisfy the minimum distance constraint. When it does, we say that the individual is feasible. This verification is done using the following sub-routine.

This routine takes as input an individual, , and returns 1 if the individual is feasible, or 0 otherwise. The procedure simply computes the Euclidean distance between every pair of coordinates that define the individual, and checks whether the distance is at least.

The routine INITIAL-POPULATION can be divided into two main blocks: lines 1 - 12 and 13 - 29. The first block creates a “greedy” solution to the problem by choosing the coordinates of, the first individual of the population, based on the vector of sorted values of. The second block creates the rest of the initial population by selecting random locations on the grid (using the function, which returns a random element of a finite set), with the condition that each individual created is feasible.

We note that, for large values of, it may be hard, and perhaps even impossible, to find locations on the grid such that the minimum distance between each pair of locations is. Our algorithm makes attempts to find such locations, and, in the case that it is not successful, some of the coordinates of the individual being created remain as.

The sub-routine CROSSOVER performs the reproduction of individuals. In this routine, two individuals of the population are selected at random (lines 2 - 3), and the crossover point, represented by, is also selected at random (line 4). In lines 5 - 12, the crossover occurs, and two new individuals are created. These new individuals are then checked for feasibility, and, in the case that they are not feasible, they are “erased” from the population by setting all of their coordinates to (lines 13 - 21). The procedure is repeated times, and the routine creates new individuals. 

Our GA then performs a local search in the grid, in the geographical sense, to try to improve every individual of the population.

In line 2, a coordinate of individual is chosen at random, and, in case it is not, a local search in the reservoir grid, within a square of area centered at, is performed to try to find a location in the search area such that, if substituted in place of, will maintain feasibility of the individual and increase its associated value the most. The routine performs this local search for, and returns their updated values.

The next genetic operator is mutation, which, again, is performed for every individual of the population.

This routine picks a coordinate from individual at random and replaces it with another coordinate picked at random in the entire grid. If the replacement keeps feasible and improves the value of, then is updated; otherwise, remains unchanged.

Our last genetic operator, INTRUDERS, creates individuals that enter the population “from outside”. The intention of this procedure is to add variability to the population thus far created.

This routine is similar to INITIAL-POPULATION. The main differences are in lines 2 and 6, where the choice of the coordinates of the new individuals is more selective: they are chosen among the top values of the matrix. Thus, this routine brings to the population not only variability, but also members with high fitness.

Our final sub-routine selects and preserves the best fit individuals from the entire population, which will become the initial population of the next generation.

The sorting procedure in line 1 is, in fact, a renumbering of the individuals based on the values,. As a clarifying example, suppose we give as input to this line, , and , with, , and. Then the output will be, , and. Lines 2 - 7 simply erase variables and keep which are the surviving individuals that will start a new generation. (In case, then there will not be another generation, and WPP-GA will return, the individual with the best fitness found after generations.)

A careful analysis of WPP-GA and its sub-routines shows that the algorithm terminates, but not necessarily with an optimal solution to the problem, i.e., an individual whose fitness is the highest possible.

3. Mathematical Optimization

Mathematical Optimization approaches quantitative problems using tools such as linear algebra, calculus, and graph theory. In a sense, it is a more sophisticated method than Evolutionary Metaheuristics, and, in practice, usually requires bigger and more structured algorithms to solve a given problem.

The main advantage of using mathematical optimization, and, in our particular case, IP, is that a solution may be proven to be optimal, as opposed to GA, which does not guarantee optimality of the best solution found.

It is common practice, in IP, to formulate problems by defining an objective function to be maximized (or minimized), subject to constraints that define the problem in question. Here, the formulation of WPP is as follow:

where a variable represents the location in the oil reservoir grid, and is equal to 1 if the location is chosen, or 0 otherwise. Hence, a feasible solution to this problem is a -vector of 0’s and 1’s that satisfies all constraints defined in the formulation shown. Note that this is a completely different way of representing a solution, in comparison with our GA. Nonetheless, the values of the objective function, here, and the fitness of the individuals, in the GA, are comparable, since they have the same meaning.

To solve our instances with the IP approach, we use a commercial solver based on branch-and-cut that can handle large-scale optimization problems. For details on branch-and-cut, the reader is referred to [11].

4. Computational Tests

We ran our computational tests in the Texas Tech High Performance Computing Center’s 3.0 GHz CPUs with 16 GB RAM nodes. For the tests using IP, we used CPLEX 12 with default settings, except for the search strategy, which was set to Traditional Branch-and-Cut, and the relative and absolute gap tolerances, which were set to zero. Every test had a time limit of 3600 seconds.

Our test grid has, and we set the values of to, and. For each of these values, we tried to solve WPP with set to, up to the point at which either the time limit was reached, or the value of was not reached, i.e., neither method was able to find a solution that represented well locations. For the GA, we set, , , , and for all tests.

The values of the matrix were obtained with the commercial reservoir simulator Eclipse, using the Quality Maps approach (see [12] for details). Some main characteristics of our heterogeneous and anisotropic reservoir are: initial pressure of, porosity of 22%, and horizontal permeability average of with standard deviation of. The distance between two closest grid points is, and the thickness of the reservoir is.

Table 1 shows the results obtained using both GA and IP. The columns “Best Sol’ n” displays the value of the best solutions found by both methods. (Note: in GA, this value is called the fitness of the best individual, whereas in IP it is called the objective function value of the best solution. We will use the latter expression in our analysis.) “Time” shows the computational time, in seconds, required for the test to finish. We note that, because of the settings used in CPLEX, all solutions obtained using IP are provably optimal precisely when the time is less than 3600 seconds. The column “Card” shows the cardinality of the solution found, i.e., the number of well locations associated with that solution. Finally, “Gap” refers to the relative percentage difference between the GA and IP best solutions, and is computed using Equation (2).


From Table 1, it is clear that the problem gets harder as the values of and increase. On the IP section of the table, this is reflected in the time required to solve the instances. On the GA side, both the computational time and the cardinality of the best solution show the increase in the difficulty of the problem. For instance, for and, the best solution that GA can find has cardinality 117, although solutions with greater cardinality (and better objective function value) do exist, as the IP results show.

Overall, IP was capable of finding better solutions than GA, and, most of the time, was faster. In every test, the best solution found by IP had an objective function value at least as good as the one found by the GA. In 8 cases, both methods found an optimal solution to the problem, and in 11 other cases, the gap between the two methods was less than 1%. This shows that GA can be effective for instances that are less challenging, but, for the harder instances, IP was notably more successful than GA.

The merit of being faster and finding better solutions is due not only to the differences in approaching the problem (IP vs. GA), but also a result of the algorithms used to solve the instances. While our GA has less than 500 lines of code, CPLEX is the product of professional software developers who worked on it for decades. The seemingly small values in the last columns reflect this contrast of algorithms: a gap as small as 0.1% can be very difficult to close, and thus finding a sub-optimal solution can be a much easier task than finding an optimal one. Moreover, the fact that CPLEX can guarantee the optimality of a solution is a feature that a GA does not have. It is only in comparison studies such as the present study, that one can visualize the effectiveness of a GA, in finding optimal solutions. Finally, we note that CPLEX is a multi-purpose solver that can handle a vast number of different problems, while our GA was developed specifically to tackle WPP.

5. Conclusions and Further Research

In this study, we compared the performance of GA and IP to solve instances of the problem of placing vertical wells in an grid that sits on a reservoir, with the

Table 1. Summary of results using GA and IP.

objective of extracting the maximum amount of oil from the reservoir. Our results indicate that GA can be effective for easier instances, but it lacks performance for harder ones. In comparison, CPLEX (which takes the IP approach) always found a solution at least as good as our GA, and most of the time faster. Moreover, IP has the advantage of finding provably optimal solutions, while GA is not able to guarantee that a solution is optimal.

The GA presented here can, of course, be modified and have its performance enhanced, either by actually changing the algorithm, or by simply adjusting its parameters. Similarly, the solution via IP can be made faster by considering different formulations of the problem and tuning some parameters on CPLEX. To maintain both approaches simplely, we chose not to vary those parameters, to provide a basic GA, and to use a straightforward formulation for our IP approach.

As further research, we suggest the study of non-conventional wells placement, a more challenging problem, since it involves wells that are not necessarily vertical, and thus may require a 3D-grid as the set of possible locations for the wells.

6. Acknowledgements

The authors acknowledge the High Performance Computing Center at Texas Tech University at Lubbock for providing resources that have contributed to the research results reported within this paper. URL: http://www.hpcc.ttu.edu We thank Dr. Jerome Onwunalu’s for guidance and support, and for providing a GA code used in preliminary tests.

This research was partially supported by the Office of Naval Research (ONR) through grant N000141310041 to de Farias.

Conflicts of Interest

The authors declare no conflicts of interest.


[1] G. AlQahtani, R. Vadapalli, S. Siddiqui and S. Bhattacharya, “Well Optimization Strategies in Conventional Reservoirs,” Proceedings of SPE Saudi Arabia Section Technical Symposium and Exhibition, Al-Khobar, 8-11 April 2012, 13 Pages.
[2] O. Badru and C. S. Kabir, “Well Placement Optimization in Field Development,” SPE Annual Technical Conference and Exhibition, Denver, 5-8 October 2003, 9 Pages.
[3] W. Bangerth, H. Klie, M. F. Wheeler, P. L. Stoffa and M. K. Sen, “On Optimization Algorithms for the Reservoir Oil Well Placement Problem,” Computational Geosciences, Vol. 10, No. 3, 2006, pp. 303-319.
[4] A. S. Cullick, S. Vasantharajan and M. W. Dobin, “Determining Optimal Well Locations from a 3D Reservoir Model,” US Patent No. US6549879 B1, 2003.
[5] D. Y. Ding, “Optimization of Well Placement Using Evolutionary Algorithms,” SPE EAGE Annual Technical Conference and Exhibition, Rome, 9-12 June 2008, p. 912.
[6] B. Guyaguler, R. N. Horne, L. L. Rogers and J. J. Rosenzweig, “Optimization of Well Placement in a Gulf of Mexico Waterflooding Project,” SPE Annual Technical Conference and Exhibition, Vol. 5, No. 3, 2000, pp. 110-118.
[7] J. Onwunalu, M. Litvak, L. J. Durlofsky and K. Aziz, “Application of Statistical Proxies to Speed up Field Development Optimization Procedures,” International Conference and Exhibition, Abu Dhabi, 3-6 November 2008, 14 Pages.
[8] G. W. Rosenwald and D. W. Green, “A Method for Determining the Optimum Location of Wells in a Reservoir Using Mixed-Integer Programming,” SPE Journal, Vol. 14, No. 1, 1974, 12 Pages.
[9] S. Vasantharajan and A. S. Cullick, “Well Site Selection Using Integer Programming Optimization,” Proceedings of IAMG’97, 22-26 September 1997, pp. 421-426.
[10] R. Haupt and S. Haupt, “Practical Genetic Algorithms,” John Wiley and Sons, Hoboken, 2004.
[11] G. L. Nemhauser and L. A. Wolsey, “Integer and Combinatorial Optimization,” John Wiley and Sons, Hoboken, 1988.
[12] P. S. Cruz, R. N. Horne and C. V. Deutsch, “The Quality Map: A Tool for Reservoir Uncertainty Quantification and Decision Making,” SPE Annual Technical Conference and Exhibition, Vol. 7, No. 1, 2004, 9 Pages.

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.