Quantum-Based Analytical Techniques on the Tackling of Well Placement Optimization

: The high dimensional, multimodal, and discontinuous well placement optimization is one of the main di ﬃ cult factors in the development process of conventional as well as shale gas reservoir, and to optimize this problem, metaheuristic techniques still su ﬀ er from premature convergence. Hence, to tackle this problem, this study aims at introducing a dimension-wise diversity analysis for well placement optimization. Moreover, in this article, quantum computational techniques are proposed to tackle the well placement optimization problem. Diversity analysis reveals that dynamic exploration and exploitation strategy is required for each reservoir. In case studies, the results of the proposed approach outperformed all the state-of-the-art algorithms and provided a better solution than other algorithms with higher convergence rate, e ﬃ ciency, and e ﬀ ectiveness. Furthermore, statistical analysis shows that there is no statistical di ﬀ erence between the performance of Quantum bat algorithm and Quantum Particle swarm optimization algorithm. Hence, this quantum adaptation is the main factor that enhances the results of the optimization algorithm and the approach can be applied to locate wells in conventional and shale gas reservoir.


Introduction
In oil and gas industry, to maximize productivity, optimization of well placement is one of the most important issues in the field development process. The optimization techniques, based on contemporary research works on this field, can be indexed in three major categories. They are (i) classical methods, (ii) non-classical methods, and (iii) hybrid methods [1]. In the early novel research endeavors, focuses were mainly shown on classical methods to tackle well placement optimization problem and among them mixed-integer programming (MIP) [2], gradient-based finite difference method [3], multivariate interpolation algorithms [4], simultaneous perturbation stochastic approximation method [5,6], and steepest ascent method [7] are significant. However, the gradient-based techniques might get trapped in local optima, which is one of the biggest bottlenecks of the gradient-based techniques. Moreover, the calculation of the gradients is complex. Hence,

Governing Equations
An oil reservoir's dynamic behavior can be described by coupled spatiotemporal differential equations. In the case of a two-phase oil reservoir, these are [38]: where, µ f is the viscosity of phase f (f = 0 for oil and f = w for water), ε is the porosity, ρ f is the density, B f is the formation volume factor, K is the absolute permeability tensor, S f is the saturation, kr f is the relative permeability, P f is the pressure, and q f is the flow from (+ve for out, −ve for in) the reservoir. In mathematical analysis and modeling, the relative permeability of each phase is generally described using: where a and b present the exponents in Corey's correlation, S wr and S or are the residual water and oil saturations, respectively, kr 0 w and kr 0 o are the end-point relative permeabilities for water and oil, respectively.
The pressure and saturation of the oil and water phases are interrelated by the equation: S 0 + S w = 1; where P c is capillary pressure, P o and P w are oil and water pressure, and S 0 and S w presents oil and water saturation.
Then, the flow rate at the wellbore can be defined using the following well production equation (IPR or Inflow Performance Relation): q = Ψ (kr/µB(P − BHP). (6) where Ψ is the connection transmissibility factor for wells and the Bottom Hole Pressure (BHP) is determined by the Vertical Flow Performance (VFP) curve: Among them, ∆P presents the frictional pressure drop through well tubing, ∆p is the pressure drop due to the acceleration, L represents the depth of the well, and ρ is the density of the well production, and THP is tubing head pressure. For a given well chain and termination type, the total flow, the water/oil ratio (gas/oil ratio), and the inlet/outlet pressure determine the pressure drop in Equation (7).
The prime motivation behind well placement optimization is to make sure that the expenditure remains minimum while maximizing the net present reaches the maximum value. Well placement optimization, in general, can be formulated as: Max R(u n ) (8) R(u n ) = NPV(u n ) (9) Subjected to: LB ≤ u n ≤ UB∀ n (0, 1, 2, 3 . . . . . . N − 1), where u n represents well coordinates, NPV presents net present value, LB and UB are lower bound and upper bound of the reservoir, respectively. Net present value (NPV) changes randomly with the change of coordinates of the well location. Eclipse simulation was used to calculate the cumulative oil production, cumulative gas production, cumulative water production's value based on the coordinates of the well location. The variables used in Equation (11) are depicted from Reference [1]. Hence, NPV for a reservoir model can be formulated as: where P g denotes gas price, Q w presents cumulative water production, D is the discount rate, Q g is cumulative gas production, T is the number of years passed since the production has started, CAPEX is the capital expenditure, P O presents oil price, C w denotes cost per unit volume of produced water, OPEX is the operational expenditure, and Q O is cumulative oil production. Equation (1) constitutes an exact conventional as well as shale gas reservoir model, which cannot be solved analytically. Hence, by using complex commercial simulators like ECLIPSE, CMG to solve Equation (1) to Equation (7), production data Q O and Q g in Equation (11) can be found. Figure 1 depicts the general flow chart to acquire NPV value.
In general, Figure 1 shows the flow of the program to search the maximum net present value. At the start of the program, the parameters of the algorithm are initialized. Then, the program inserts a generated location in the data file. After that, the program calls eclipse and eclipse provide run summary file for a given location. The program will extract production data Q O and Q g from run summary file and calculate NPV. Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 25 In general, Figure 1 shows the flow of the program to search the maximum net present value. At the start of the program, the parameters of the algorithm are initialized. Then, the program inserts a generated location in the data file. After that, the program calls eclipse and eclipse provide run summary file for a given location. The program will extract production data and from run summary file and calculate NPV.

Methodology
Metaheuristic algorithms are stochastic and non-deterministic. There are various types of metaheuristic methods like local search, simulated annealing, Tabu search, variable neighborhood search, population-based or trajectory-based search, etc. Besides, among all the Metaheuristic search algorithms, one of the most popular algorithms are-gravitational search algorithm (GSA), particle swarm optimization algorithm (PSO), crow search algorithm (CSA), genetic algorithm (GA), differential evaluation (DE), and bat algorithm (BA), etc. Quantum-based optimization techniques were applied in several complex engineering applications using quantum parallelism mechanisms. For multimodal optimization applications, the study of Ross indicated that quantum behaved algorithm is superior to existing metaheuristic algorithms [39,40]. Another study found that quantum computation can manage highly non-linear multimodal optimization problems naturally [39]. Again, PSO searches linearly. On the other hand, QPSO's next position depends entirely on the probabilistic approach [41]. In QBA, the position of each bat depends on the mean best position. The quantumbehaved bat and mean best position lead local search to jump out of local positions [42]. Hence, QBA can easily avoid local optima. Since these properties are important to tackle multimodal well placement optimization problem, QPSO and QBA are discussed in the following section.

Methodology
Metaheuristic algorithms are stochastic and non-deterministic. There are various types of metaheuristic methods like local search, simulated annealing, Tabu search, variable neighborhood search, population-based or trajectory-based search, etc. Besides, among all the Metaheuristic search algorithms, one of the most popular algorithms are-gravitational search algorithm (GSA), particle swarm optimization algorithm (PSO), crow search algorithm (CSA), genetic algorithm (GA), differential evaluation (DE), and bat algorithm (BA), etc. Quantum-based optimization techniques were applied in several complex engineering applications using quantum parallelism mechanisms. For multimodal optimization applications, the study of Ross indicated that quantum behaved algorithm is superior to existing metaheuristic algorithms [39,40]. Another study found that quantum computation can manage highly non-linear multimodal optimization problems naturally [39]. Again, PSO searches linearly. On the other hand, QPSO's next position depends entirely on the probabilistic approach [41]. In QBA, the position of each bat depends on the mean best position. The quantum-behaved bat and mean best position lead local search to jump out of local positions [42]. Hence, QBA can easily avoid local optima. Since these properties are important to tackle multimodal well placement optimization problem, QPSO and QBA are discussed in the following section.

Quantum Particle Swarm Optimization Algorithm
In 1995, a novel search approach was developed and proposed by Kennedy and Eberhart [43], which is known as PSO and was stochastic. In a search space, the particles of PSO can fly and change their position with respect to time. At first, the initialization of individuals is done through a random process.
To update the velocity of each particle, the following equation is utilized where the individual i particles velocity is denoted by V i for iteration k. For defining a particle's acceleration, c 1 and c 2 are introduced as acceleration constant, w is used for weight vector, rand 1 and rand 2 are the random number ranging from 0 to 1, for individual i position is represented by x k i in iteration k, pbest k i and gbest k denotes the best position and global best of individual i in k iteration. In the search space, by using the following equation, each particle's next position will be updated: Sun et al. [44] implemented the principle of quantum mechanics with basic PSO. The QPSO algorithm not only tackles the drawbacks but also preserves the good features of the PSO algorithm [45], and thus the incorporation of improved search capability in addition to fast convergence is possible [46]. Another feature of QPSO is that it has only one parameter that needs to be tuned. The flowchart of QPSO is depicted in Figure 2. Figure 2 shows that at first the positions are randomly generated. Then, the cost function is calculated for each particle. Based on the results, mbest and gbest are updated. The search process will continue unless the termination criteria are met.
In 1995, a novel search approach was developed and proposed by Kennedy and Eberhart [43], which is known as PSO and was stochastic. In a search space, the particles of PSO can fly and change their position with respect to time. At first, the initialization of individuals is done through a random process.
To update the velocity of each particle, the following equation is utilized where the individual i particles velocity is denoted by for iteration k. For defining a particle's acceleration, are introduced as acceleration constant, w is used for weight vector, are the random number ranging from 0 to 1, for individual i position is represented by in iteration k, and denotes the best position and global best of individual i in k iteration. In the search space, by using the following equation, each particle's next position will be updated: Sun et al. [44] implemented the principle of quantum mechanics with basic PSO. The QPSO algorithm not only tackles the drawbacks but also preserves the good features of the PSO algorithm [45], and thus the incorporation of improved search capability in addition to fast convergence is possible [46]. Another feature of QPSO is that it has only one parameter that needs to be tuned. The flowchart of QPSO is depicted in Figure 2. Figure 2 shows that at first the positions are randomly generated. Then, the cost function is calculated for each particle. Based on the results, mbest and gbest are updated. The search process will continue unless the termination criteria are met. Unlike the behavior of the Newtonian particle that is described by both the position ( ⃗) and velocity ( ⃗) in traditional PSO, in case of QPSO to define the quantum state of a particle, the aid of wave function Ψ( ⃗, s) will be taken. By expanding the wave function in the three-dimensional space, we can get, Unlike the behavior of the Newtonian particle that is described by both the position ( → x ) and velocity ( → v ) in traditional PSO, in case of QPSO to define the quantum state of a particle, the aid of wave function Ψ ( → x , s) will be taken. By expanding the wave function in the three-dimensional space, we can get, where Rdxdydz is a probabilistic measure indicating the appearance of that particle at specific time s at (x, y, z) point. It should be noted that the square of the wave function (|Ψ | 2 ) is the probability density of the particle at that specific point satisfying the following expression: Appl. Sci. 2020, 10, 7000 7 of 25 Time variant 3D wave function Ψ ( → x ,s), which can be interpreted by (14) or (15), satisfies the following time-dependent Schrödinger equation [47]: In this case,h andĥ represents constant of Planck and Hamiltonian operator, respectively. Considering the particle in a potential field V with mass V(s) we can write: Assuming each particle as QPSO, Sun and associates [47] treat the particle with no spin in D dimensional Hilbert space with a given energy. Hence, their states can be characterized by the wave function. Now for D dimensional quantum space having a population that consists of k particles, ith particle's location can be defined by X i = (x i1 , x i2 , . . . , x iD ). For the ith particle, the best solution's position in previous step, i.e., pbest could be denoted by Qi = (Q i1 , Q i2 , . . . , Q iD ). Similarly, in search space of all the particles, the best particle's position, i.e., gbest can be written as Q g = (Q g1 , Q g2 , . . . , Q gD ). Using Monte-Carlo method, the quantum state of the particle's position could be expressed as [45]: where i = 1, 2, . . . , n, dimension d = 1, 2, . . . , D; u is defined as a random number ranging [0, 1]; q id denotes local attractor of ith particle on d dimension that could be expressed as [48]: where ϕ is defined as a random number, which is distributed uniformly ranging [0, 1]. A numerical value, L, originated from the current position of particle and best position of an individual that could be written as L = 2·β q id − x id . Now the quantum state of the particle's position in (18) will be expressed as: where contraction expansion (CE) is referred to as β, the only QPSO parameter earlier mentioned. The adaptive CE coefficient could be expressed by [48]: where t max and t are expressed for the maximum number and the current number of iterations, respectively. In the traditional PSO algorithm, premature convergence is common as a position, as well as the velocity of the particle, is directly used. However, to encounter this issue of PSO, in the QPSO algorithm Sun et al. [44] proposed mbest. Position of the particle is denoted by mbest and the particle's best position can be expressed as: Appl. Sci. 2020, 10, 7000 8 of 25 where particle i's best possible position is denoted by Q i . If mbest is used, then Equation (20) is transformed as follows:

Quantum-Behaved Bat Algorithm
The basic Bat Algorithm is formulated by using idealized rules. The three rules are, (i) The use of echolocation capability to measure distance and to sense the difference between their prey (food) and background barriers are common to every bat. (ii) If the bats are in x i position having velocity v i , with f min fixed frequency as well as λ 0 varying wavelength, then loudness A 0 will be used for the search of food while flying randomly. Automated adjustment of the wavelength of their emitted pulses as well as adjustment of the rate of pulses r [0, 1] can be done by bats which is dependent on target proximity, (iii) It is presumed that loudness changes from a positive large value A 0 to a minimum value A min that is constant.
Assuming that the solution is not known, to initialize bats randomly, in the search space of each dimension both lower and upper bound is used. The common solution will be generated by using: For ith bat in the jth dimension X 0 and X m denotes the lower and upper bound, respectively. Considering this, frequency, velocity, and position of bats can be formulated in the following manner: where α is random vector ranging [0, 1]; f i , f min as well as f max are pulse frequency, minimum frequency, and maximum frequency, respectively. Moreover, v t i , v t−1 i , x t i , x t−1 i , and g t denote ith bat's velocity for t iteration, ith bat's velocity at (t − 1) iteration, ith bat's position at t iteration, ith bat's position at (t − 1) iteration, and global best location currently found by bats, respectively.
For the generation of new solutions for every bat, a random natured local walk is used once a suitable solution is selected from the recent best solutions. Hence, a new position may be expressed as follows: where ε is denoted for expressing a random number ranging [−1, 1] and A t expresses average loudness of all the bats for total t iterations.
Loudness A i and pulse rate denoted by r i are the parameters that control the flow of the bat algorithm. For every iteration, loudness A i and pulse rate r i can be found using the following expression: where A t i , A t+1 i , r 0 i and r t+1 i denoted for ith bat's loudness for tth iteration, ith bat's loudness for (t + 1)th iteration, ith bat's initial pulse rate, and ith bat's pulse rate for (t + 1)th iteration, respectively. Both ∆ and γ are values that are constant ranging [0, 1] and positive, respectively.
The Quantum-behaved bat algorithm (QBA) is the improved variant of the original bat algorithm. Pioneered in 2010 by Xin-she Yang [49], the original bat algorithm was based on the echolocation or bio-sonar capability of bats. The flowchart of QBA is depicted in Figure 3. The flowchart in Figure 3 shows that after random initialization loudness, frequency and plus rate are updated. After that, the search equation is employed to update search location. A local search is incorporated into the search technique. This process will continue unless the termination criteria are satisfied.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 25 Unlike the traditional bat algorithm, for QBA, a new position is assessed differently, and it can be expressed by the equations below: where (0, ) denotes a Gaussian distribution that has zero mean and standard deviation of . and are global best position in d dimension for ith bat's at t + 1 iteration. is ith bat's loudness for tth iteration. To make sure that the standard deviation always remains positive, ε is introduced.
The global best in the swarm of bats can be treated as attractor using the Equation (23) and replacing by . The position of the bat with quantum behavior can be expressed by the following equation: where is the ith bat's position when it is in d dimension for tth iteration. By taking into consideration every bat's compensation, capability for Doppler Effect that is self-adaptive in nature reforms the updating formulas as per the following expression (23)(24).
It should be mentioned that for updating the value of velocity as well as for controlling the inheriting rate of the previous velocity, the inertia weight parameter (w) should be added. The Unlike the traditional bat algorithm, for QBA, a new position is assessed differently, and it can be expressed by the equations below: where j 0, σ 2 denotes a Gaussian distribution that has zero mean and standard deviation of σ 2 . x t+1 id and g t d are global best position in d dimension for ith bat's at t + 1 iteration. A t i is ith bat's loudness for tth iteration. To make sure that the standard deviation σ 2 always remains positive, ε is introduced.
The global best g t d in the swarm of bats can be treated as attractor using the Equation (23) and replacing q id by g t d . The position of the bat with quantum behavior can be expressed by the following equation: where x t id is the ith bat's position when it is in d dimension for tth iteration. By taking into consideration every bat's compensation, capability for Doppler Effect that is self-adaptive in nature reforms the updating formulas as per the following expression (23) and (24).
It should be mentioned that for updating the value of velocity as well as for controlling the inheriting rate of the previous velocity, the inertia weight parameter (w) should be added. The compensation rate C should be varied for each bat. Since the velocity of the sound through the normal air is assumed as 340 m/s the expression could be reformed as follows: where f id is used to denote ith bat's frequency at dimension d; v t−1 g is denoted for (t − 1)th iteration of the globally best position's the velocity and C i is ith bat's respective positive number ranging [0, 1]. To keep the computation simple, we have assumed that C = 0 leads to the fact that bats will not be able to compensate for Doppler effect in echoes but in the case of C = 1, the bat will be fully capable for compensating Doppler effect in the echoes.

Advantage and Disadvantage
From No Free Lunch theorem (NFL) it is known that no single algorithm can be best for all problems [50]. After analyzing the major attributes of QBA and QPSO algorithm, three key points for its success can be underscored in the following: • Better performance can be obtained from quantum-behaved algorithms compared to PSO, GA, CSA, GSA, and DE to carry out highly nonlinear, multi-modal optimization problem as quantum-behaved algorithms has the inherent capability to increase diversity in its population.

•
As PSO and GA update their location depending not only on personal best information but also on explicit global best, premature convergence is common in theses algorithms and to avoid this, mean best is used in this algorithm.
According to the No Free Lunch theorem (NFL), a single algorithm cannot be the best solution for every problem [50]. Table 1 delineates the overall advantage and disadvantages of these algorithms in brief.

Experimental Setting
To validate the proposed methods that are discussed previously, 2 case studies are conducted. Each algorithm ran 16 times. In case study 1, we have used 100 iterations and case study 2 we have used 30 iterations. In case study 1, the search space is 19 × 28 × 5 blocks are present. Again, in case study 2, 10 × 10 × 3 grid blocks are present. Hence, in case study 1, the search space is larger than case study 2. So, algorithms will require more effort to find optimum value in case study one. For PUNQ S3 reservoir 100 iterations is used in the following paper [51]. The detailed description of test case studies are shown in Figures 4 and 5. Again, Figure 6b shows that in case study 2, after a certain iteration the algorithms have reached their optimal result and did not improve further. Furthermore, after a certain point, the curve shows a flat line. The flat line in Figure 6b proves that the algorithms have converged and will not improve further. Hence, 30 iteration is sufficient for case study 2. In each trial, for case study 1, 100 iterations are conducted with 20 particles. A total of 2000 function evaluations are considered as a stopping criterion. In each trial, for case study 2, 30 iterations are conducted with 5 particles. A total of 150 function evaluations are considered as a stopping criterion. The parameters that were used in these algorithms and the economic parameters that were used in Equation (11) to conduct the experiment are listed in Tables 2 and 3.

Description of Case Studies
To select the test cases, the reservoirs that are not identical are considered. The surface of the reservoir can be highly non-smooth or it can multimodal for well placement optimization problem. The PUNQ S3 reservoir can be highly multimodal [26]. So, it will test capacity to tackle the multimodal optimization problem. Again, SPE-1 is fairly simple. So, the case studies will test the capacity to tackle highly multimodal optimization problem and simple multimodal optimization problem.
PUNQ-S3 is a real field-based model and had been used for Elf Exploration production. The detailed characteristics of this reservoir model are analyzed in the literature [58] and shown in Figure  4. The PUNQ-S3 consists of 19 × 28 × 5 blocks of the grid. To evaluate the proposed technique, the authors have chosen 4 vertical wells that need to be optimized, where each well has ( , ) coordinates. Thus, the overall number of the variables that need to be optimized in this experiment is 2 × 4.   To conduct the second case study, the first SPE model is adopted. The first SPE model is formulated based on a 3D black oil reservoir simulation. Hence, this type of synthetic reservoir model is exploited in this study. The detailed characteristics and specification of this reservoir model are discussed in the literature [59] and shown in Figure 5. There are 10 × 10 × 3 grid blocks in the First SPE model. Two vertical wells locations are optimized in this case study. Since, each well is defined by (x, y) coordinates, 2 × 2 variables are optimized in this case study.
(a) Initial pressure (b) Initial Gas saturation (c) Initial Oil saturation (d) Initial Water saturation

Convergence Analysis
To analyze the performance of different algorithms, the convergence curve and convergence speed is an important tool. So, plots of the mean net present value of algorithms versus iterations number are illustrated in Figure 6. As per the illustration of Figure 6a, QPSO algorithm surpassed all the other algorithms. The QBA provided 2nd best results in case study 1. The 3rd and 4th best net present values are achieved by the DE and PSO algorithm, respectively. However, it should be   To conduct the second case study, the first SPE model is adopted. The first SPE model is formulated based on a 3D black oil reservoir simulation. Hence, this type of synthetic reservoir model is exploited in this study. The detailed characteristics and specification of this reservoir model are discussed in the literature [59] and shown in Figure 5. There are 10 × 10 × 3 grid blocks in the First SPE model. Two vertical wells locations are optimized in this case study. Since, each well is defined by (x, y) coordinates, 2 × 2 variables are optimized in this case study.

Convergence Analysis
To analyze the performance of different algorithms, the convergence curve and convergence speed is an important tool. So, plots of the mean net present value of algorithms versus iterations number are illustrated in Figure 6. As per the illustration of Figure 6a, QPSO algorithm surpassed all the other algorithms. The QBA provided 2nd best results in case study 1. The 3rd and 4th best net present values are achieved by the DE and PSO algorithm, respectively. However, it should be

Proposed -QBA
The maximum and minimum inertia weight (w max and w min ) 0.9 and 0.5 The maximal and minimal frequency (f max and f min ) 1.5 and 0 Gamma, γ 0.9 Delta, δ 0.99 The frequency of updating the loudness and emission pulse rate, G 10 The maximum and minimum contraction expansion coefficient (β max and β min ) 1 and 0.5 The maximal and minimal loudness 2 and 1 The maximal and minimal pulse rate 1 and 0 The maximum and minimum compensation rate for Doppler effect (C max and C min ) 1 and 0.9 The maximum and minimum probability of habitat selection 0.9 and 0.6

4.
Proposed -QPSO  Figure 6b, QBA has exhibited superior performance to provide a better net present value compared to other algorithms in case study 2. After QBA, QPSO and PSO algorithms have gained better results than other methods. However, unlike the previous case study, GA performed the worst with respect average net present value. Furthermore, the CSA, GA, BA, and GSA algorithms again failed to obtain satisfactory NPV in this study. Figure 6 suggests that, unlike PSO and BA algorithm, modification in the search expression in PSO and BA with quantum computation helped avoid local optimum in convergence. Furthermore, the quantum computation in the PSO and BA made it more vibrant across both case studies. Hence, in a multimodal problem, QPSO and QBA is more dynamic than PSO and BA and is further validated by the exploration and exploitation graphs shown in Figure 6.

Performance Evaluation and Statistical Analysis
To assess the discussed algorithms performance, several benchmark criteria are considered for this problem [6,60]. Thus, the author's preferred benchmark for the evaluation are discussed in the following subsection: Effectiveness is the average value that is in between consecutive tests where the best solution is expressed with respect to the global optima as a percentage is effective and it is simply a measure of performance. Mathematically, effectiveness can be defined as follows:

Description of Case Studies
To select the test cases, the reservoirs that are not identical are considered. The surface of the reservoir can be highly non-smooth or it can multimodal for well placement optimization problem. The PUNQ S3 reservoir can be highly multimodal [26]. So, it will test capacity to tackle the multimodal optimization problem. Again, SPE-1 is fairly simple. So, the case studies will test the capacity to tackle highly multimodal optimization problem and simple multimodal optimization problem.
PUNQ-S3 is a real field-based model and had been used for Elf Exploration production. The detailed characteristics of this reservoir model are analyzed in the literature [58] and shown in Figure 4. The PUNQ-S3 consists of 19 × 28 × 5 blocks of the grid. To evaluate the proposed technique, the authors have chosen 4 vertical wells that need to be optimized, where each well has (x, y) coordinates. Thus, the overall number of the variables that need to be optimized in this experiment is 2 × 4.
To conduct the second case study, the first SPE model is adopted. The first SPE model is formulated based on a 3D black oil reservoir simulation. Hence, this type of synthetic reservoir model is exploited in this study. The detailed characteristics and specification of this reservoir model are discussed in the literature [59] and shown in Figure 5. There are 10 × 10 × 3 grid blocks in the First SPE model. Two vertical wells locations are optimized in this case study. Since, each well is defined by (x, y) coordinates, 2 × 2 variables are optimized in this case study.

Convergence Analysis
To analyze the performance of different algorithms, the convergence curve and convergence speed is an important tool. So, plots of the mean net present value of algorithms versus iterations number are illustrated in Figure 6. As per the illustration of Figure 6a, QPSO algorithm surpassed all the other algorithms. The QBA provided 2nd best results in case study 1. The 3rd and 4th best net present values are achieved by the DE and PSO algorithm, respectively. However, it should be mentioned that both GA, GSA, BA, and CSA have suffered the problem to converge in local optimum point. Overall, in the first case study, QPSO and QBA have achieved better solutions. As seen from Figure 6b, QBA has exhibited superior performance to provide a better net present value compared to other algorithms in case study 2. After QBA, QPSO and PSO algorithms have gained better results than other methods. However, unlike the previous case study, GA performed the worst with respect average net present value. Furthermore, the CSA, GA, BA, and GSA algorithms again failed to obtain satisfactory NPV in this study. Figure 6 suggests that, unlike PSO and BA algorithm, modification in the search expression in PSO and BA with quantum computation helped avoid local optimum in convergence. Furthermore, the quantum computation in the PSO and BA made it more vibrant across both case studies. Hence, in a multimodal problem, QPSO and QBA is more dynamic than PSO and BA and is further validated by the exploration and exploitation graphs shown in Figure 6.

Performance Evaluation and Statistical Analysis
To assess the discussed algorithms performance, several benchmark criteria are considered for this problem [6,60]. Thus, the author's preferred benchmark for the evaluation are discussed in the following subsection: Effectiveness is the average value that is in between consecutive tests where the best solution is expressed with respect to the global optima as a percentage is effective and it is simply a measure of performance. Mathematically, effectiveness can be defined as follows: where f (p) is defined as a solution for p, p * is denoted for the solution for global optima, pˆi is referred to as the best solution for ith trial, and N is trial number for each algorithm. Efficiency is the parameter that is an indicator of the algorithm's speed at which it reached a specific performance level using a number of distinctive evaluations to obtain a solution that is at least 98 percent of the best value in the experiment.
where additionally, L 98 i is the number of unique function evaluations required to find solution q such that f (q) _ 0:98f (ˆpi) for trial i (for minimization) and M is the total number of function evaluations per trial.
Besides these two criteria, statistical data like standard deviation, average, and min-max are collected in experimental trials and the results are shown in Tables 4 and 5. 5.14 × 10 9 3.72 × 10 9 5.09 × 10 9 5.13 × 10 9 5.30 × 10 9 5.03 × 10 9 5.33 × 10 9 Min 2.83 × 10 9 3.43 × 10 9 2.43 × 10 9 3.01 × 10 9 3.38 × 10 9 2.10 × 10 9 4.38 × 10 9 4.29 × 10 9 Average 3.33 × 10 9 4.07 × 10 9 3.24 × 10 9 3.67 × 10 9 4.26 × 10 9 3.28 × 10 9 4.77 × 10 9 4.67 × 10 9 Standard deviation 2.62 ×   The mean value and standard deviation indicate the robustness of the algorithm. The results of Case Study 1 show that the QPSO algorithm is optimal on four criteria. However, QBA was better on 2 other criteria. The box plot results are shown in Figure 7a indicates that QPSO and QBA have reached almost similar results compared to other algorithms and QPSO has the lowest standard deviation.  Again, the result of case study 2 in Table 5 shows that the QPSO and QBA algorithm is superior on three and five criteria, respectively. However, PSO and DE also achieved the same maximum value, their average value is lower than QBA and QPSO. The box plot results in Figure 7b illustrates that QBA has the lowest standard deviation compared to other algorithms, and QPSO has the secondlowest standard deviation. The main results of the study of case 2 are somewhat tantamount to the previous study as it shows quantum-based algorithms are a better choice considering all 6 criteria.

Wilcoxon's Rank-Sum Test
To validate whether the results of QPSO and QBA are statistically different from the other six algorithms, a non-parametric statistical test, the Wilcoxon rank-sum test [61,62], with a significance level of 0.05, is used. Table 6 shows the p-value one tail, p-value two tails, and Z values for the Wilcoxon rank sums test. A p-value less than 0.05 and Z value higher than 1.96 is required to reject the null hypothesis. From Table 6 it can be inferred that there is a statistical difference between the performance of QPSO and QBA with other algorithms except in one case. In case study 1 and 2, the performance of the QPSO and QBA are not statistically different. Therefore, it can be inferred that quantum computation can bring improved results from the other six algorithms. However, QPSO and QBA's proved to be statistically different from its primary algorithm.  Again, the result of case study 2 in Table 5 shows that the QPSO and QBA algorithm is superior on three and five criteria, respectively. However, PSO and DE also achieved the same maximum value, their average value is lower than QBA and QPSO. The box plot results in Figure 7b illustrates that QBA has the lowest standard deviation compared to other algorithms, and QPSO has the second-lowest standard deviation. The main results of the study of case 2 are somewhat tantamount to the previous study as it shows quantum-based algorithms are a better choice considering all 6 criteria.

Wilcoxon's Rank-Sum Test
To validate whether the results of QPSO and QBA are statistically different from the other six algorithms, a non-parametric statistical test, the Wilcoxon rank-sum test [61,62], with a significance level of 0.05, is used. Table 6 shows the p-value one tail, p-value two tails, and Z values for the Wilcoxon rank sums test. A p-value less than 0.05 and Z value higher than 1.96 is required to reject the null hypothesis. From Table 6 it can be inferred that there is a statistical difference between the performance of QPSO and QBA with other algorithms except in one case. In case study 1 and 2, the performance of the QPSO and QBA are not statistically different. Therefore, it can be inferred that quantum computation can bring improved results from the other six algorithms. However, QPSO and QBA's proved to be statistically different from its primary algorithm.

Exploration and Exploitation Analysis
Dimension wise analysis sheds light on the internal behavior of the algorithms [63]. In this context, it is important to have a strong diversification to better exploit the search space. A small degree of diversification implies local convergence. To measure the diversity in each dimension, the authors have considered the equation proposed in Reference [64]: where x ij is the dimension j of the ith swarm, median(x j ) refers to the median of dimension j in the whole population, and n is the total number of the population. So, the average diversity can be calculated using Equation (40): where Div is the diversity measurement of the whole population in an iteration. Hence, Exploration and Exploitation can be measured using the following equation: where, Div max presents the maximum diversity of whole populations in one run. Without revealing the behavior of swarms in the iterative process, it is difficult to understand the end result and the value of the objective function. Therefore, to understand the end results of these techniques or convergence curve in Figure 6, this study initiates a graphical illustration to show the algorithms' exploration and exploitation behavior over the whole iteration. From Figure 8, it can be inferred that a higher NPV can be achieved with an adequate exploration and development ratio. On the other hand, these figures clearly indicate a gradual change in intensity during exploration and exploitation during the iterative process.

Exploration and Exploitation Analysis
Dimension wise analysis sheds light on the internal behavior of the algorithms [63]. In this context, it is important to have a strong diversification to better exploit the search space. A small degree of diversification implies local convergence. To measure the diversity in each dimension, the authors have considered the equation proposed in Reference [64]: where is the dimension of the ith swarm, ( ) refers to the median of dimension in the whole population, and n is the total number of the population. So, the average diversity can be calculated using Equation (40): where is the diversity measurement of the whole population in an iteration. Hence, Exploration and Exploitation can be measured using the following equation: where, presents the maximum diversity of whole populations in one run. Without revealing the behavior of swarms in the iterative process, it is difficult to understand the end result and the value of the objective function. Therefore, to understand the end results of these techniques or convergence curve in Figure 6, this study initiates a graphical illustration to show the algorithms' exploration and exploitation behavior over the whole iteration. From Figure 8, it can be inferred that a higher NPV can be achieved with an adequate exploration and development ratio. On the other hand, these figures clearly indicate a gradual change in intensity during exploration and exploitation during the iterative process.
(a) CSA for case study 1 (b) GA for case study 1   The graphs shown in Figure 6 shows that while QPSO and QBA has a high average exploration rate. On the other hand, CSA, GSA, and GA had a relatively low exploration rate. In case Study 1, the overall exploration rate of QBA and QPSO is 74% and 60%, respectively, and a quantum behaved algorithm offered a better optimal solution. Again, among all algorithms, the least average exploration group contains BA, GA, CSA, and GSA and their exploration rate was 6.6%, 8.2%, 17.4%, and 6.4%, respectively. This shows that the average diversity obtained by the other algorithms is much less than the quantum-behaved algorithms. Furthermore, it can be noted that BA, GA, CSA, and GSA obtained the least average NPV value and quantum-behaved algorithms offer a better optimal solution. Again, PSO and DE obtained 27.4% and 54.6% exploration rate, respectively.
On the other hand, in case study 2, the average diversity of BA, CSA, and DE was 85.1%, 83.3%, and 85.2%, respectively. In this case, BA provides the highest exploration rate. However, BA and CSA were unable to achieve better optimization results. Hence, it showcases that the higher exploration The graphs shown in Figure 6 shows that while QPSO and QBA has a high average exploration rate. On the other hand, CSA, GSA, and GA had a relatively low exploration rate. In case Study 1, the overall exploration rate of QBA and QPSO is 74% and 60%, respectively, and a quantum behaved algorithm offered a better optimal solution. Again, among all algorithms, the least average exploration group contains BA, GA, CSA, and GSA and their exploration rate was 6.6%, 8.2%, 17.4%, and 6.4%, respectively. This shows that the average diversity obtained by the other algorithms is much less than the quantum-behaved algorithms. Furthermore, it can be noted that BA, GA, CSA, and GSA obtained the least average NPV value and quantum-behaved algorithms offer a better optimal solution. Again, PSO and DE obtained 27.4% and 54.6% exploration rate, respectively.
On the other hand, in case study 2, the average diversity of BA, CSA, and DE was 85.1%, 83.3%, and 85.2%, respectively. In this case, BA provides the highest exploration rate. However, BA and CSA were unable to achieve better optimization results. Hence, it showcases that the higher exploration rates are not necessary to optimize the well location. Again, PSO, GSA, and GA achieved 39%, 21.93%, 14.2% exploration rates. Moreover, PSO was able to achieve 3rd best average NPV. Similarly, 1st and 2nd best average NPV provider QBA and QPSO achieved 66.4%, 71% exploration, respectively. Therefore, the correct exploration-exploitation ratio should be carefully calculated according to the search field. In general, a roughly balanced response indicates that the algorithm is more efficient. This seems to be one of the explanations for the successful results of The QPSO and QBA algorithm.
Another important finding in this study is the successful search is produced through exploration mechanisms with the effective exploration of mechanisms and appropriate balanced responses. In some cases, the two metaheuristic algorithms offer very different performances to the quality of the solution, but the balanced responses are similar (i.e., in case study 2, DE and BA's Exploration percentage = 85, exploitation percentage = 15). It simply means that there are not enough balanced solutions to produce successful performance. Therefore, operators must consider the correct diversity conditions in a balanced response to produce viable options.
This study has limitations, as it is primarily focused on optimization methods. Uncertainty quantification and history matching have not been considered. The primary emphasis is on areas such as utilize quantum computation for reservoirs, enhancing efficiency, and increasing the NPV performance of well positioning. The emphasis in this analysis is only on the algorithm for optimization.

Conclusions
In this study, quantum computation is implemented for well placement optimization and their application to well placement optimization problem is investigated. Experimental results suggested that the QBA and QPSO algorithm can find a better solution than popular algorithms. It is true that due to strong attraction to local attractor QBA has a strong local search capacity that helped it to tackle this multimodal problem more efficiently. It can be inferred that quantum computation can bring improved results compared to other six algorithms. This study showcases that higher exploration rates are not necessary to optimize the well location. However, the correct exploration-exploitation ratio is should be carefully calculated in the task area. In general, a roughly balanced response indicates that the algorithm is more efficient. This seems to be one of the explanations for the successful results of The QPSO and QBA algorithm. The main conclusions of this paper are as follows: No global search algorithm can perform well on all kinds of reservoirs. Different reservoir requires different strategy.
Quantum techniques increase exploration rate in the search technique, which helps to find better results.
There is no statistical difference in the results of QBA and QPSO. The particles are volatile and diverse owning to the local attractor points. Quantum techniques are less susceptible to premature convergence and less likely to be stuck in local optima.
High exploitation rate can be linked with lower NPV value. Diversity analysis reveals that dynamic exploration and exploitation strategy is required for each reservoir.
Future optimization studies should focus on using the global search algorithm with a local search approach because it may have the advantage of solving well placement optimization problem successfully. Lots of research papers have used proxy models in recent years to replace actual reservoir simulators, and these models have been found to minimize runtime. The accuracy of this alternative model, therefore, depends on its range of sampling. Future research in this area may focus on improving the reliability of this technique.

Rand random ∆P
The frictional pressure drop through well tubing ∆p The pressure drop due to the acceleration L The depth of the well r Pulse rate λ varying wavelength Greek Symbols γ Gamma ρ The density of the well production δ Delta β The contraction expansion coefficient µ f The viscosity of phase ε The porosity ρ f The density S f The saturation kr f The relative permeability q f The flow from the reservoir Subscripts min Minimum max Maximum o oil w water