A Surrogate Assisted Quantum-Behaved Algorithm for Well Placement Optimization

The oil and gas industry faces difficulties in optimizing well placement problems. These problems are multimodal, non-convex, and discontinuous in nature. Various traditional and non-traditional optimization algorithms have been developed to resolve these difficulties. Nevertheless, these techniques remain trapped in local optima and provide inconsistent performance for different reservoirs. This study thereby presents a Surrogate Assisted Quantum-behaved Algorithm to obtain a better solution for the well placement optimization problem. The proposed approach utilizes different metaheuristic optimization techniques such as the Quantum-inspired Particle Swarm Optimization and the Quantum-behaved Bat Algorithm in different implementation phases. Two complex reservoirs are used to investigate the performance of the proposed approach. A comparative study is carried out to verify the performance of the proposed approach. The result indicates that the proposed approach provides a better net present value for both complex reservoirs. Furthermore, it solves the problem of inconsistency exhibited in other methods for well placement optimization.


I. INTRODUCTION
Well placement is a boring process used to bring oil to the surface and placing wells in an appropriate location involves optimization techniques. Well placement optimization is a difficult task in the oil and gas industry as it creates inconsistency in the cost functions [1], [2]. It is also challenging due to the heterogeneities of reservoirs [1]. Reservoir heterogeneity is the variation of reservoir properties in space and time [3]. The surface of the search field in well placement optimization changes with the changes of reservoir heterogeneity. Furthermore, reservoirs such as PUNQ-S3 [4], [1], [5] and SPE-1 [6], have different properties and produce dynamic search spaces [7]. Hence, it is desirable to develop an effective algorithm to deal with complicated optimization problems.
The associate editor coordinating the review of this manuscript and approving it for publication was Chun-Hao Chen .

A. RELATED WORKS
The optimization algorithms used in the well placement optimization problems can be categorized into three main sections: (i) Traditional, (ii) Non-traditional, and (iii) Hybrid Techniques. Researchers initially applied traditional techniques such as the simultaneous perturbation stochastic approximation method (SPSA) [5], [6], mixed integer programming (MIP) [8], steepest ascent method [9], multivariate interpolation algorithms [10], and the finite difference method [11] to optimize the well placement problem. These traditional techniques have the vulnerability to entrap in local optima, as they use gradient information. Thus, it can be inferred that gradient-based techniques are inappropriate for well location optimization [12], [5]. In contrast, non-traditional or gradient-free techniques perform better than traditional techniques [13]- [15]. Various non-traditional techniques have been implemented such as Bat Algorithm (BA) [16], covariance matrix adaptation evolution strategy (CMA-ES) [17], firefly (FF) [18], differential evolution (DE) [19], particle swarm optimization (PSO) [20], [12], and genetic algorithm (GA) [21] to solve the well placement optimization problem. These techniques are derivative-free and provide a preferable solution for the optimization problem compared to the traditional techniques [22]. To obtain a better solution, niching techniques with Crow Search Algorithms (CSA) are used in well placement optimization problems [23]. However, the convergence capability of these techniques is poor. Furthermore, the problem of local optima for well placement optimization still abates the performance of these optimization algorithms [24]. Again, the particle swarm optimization algorithm is also incorporated with a novel weighting scheme [25]. The evaluation process used seven references data sets with different characteristics and complexity. The findings confirm that the proposed method produced the best results. Additionally, a novel populationbased optimization approach, the Aquila Optimizer (AO) is proposed in [26]. Experimental results demonstrate the superiority of the AO algorithm compared to well-known metaheuristic methods. Moreover, a study validated the Sine Cosine Algorithm's (SCA) success against related algorithms with a series of statistical tests [27]. However, the SCA does not have the ability to address the complexity of multimodal search space.
To improve the whale optimization algorithm (WOA), researchers combined multi-swarm and chaotic strategies to obtain optimized parameters and selected feature simultaneously for support vector machine (SVM) [28]. The results show that the CMWOAFS-SVM outperforms all other competitors. Another study proposes a variant of WOA, which incorporates two techniques at once [29]. The proposed EWOA (Evolutionary geography-based Whale Optimization Algorithm) has not been investigated in dynamic landscapes.Furthermore, Wang et al. [30] seek the optimal kernel extreme learning machine (KELM) using the chaotic moth-flame optimization (CMFO) approach. This technique performed better than the kernel extreme learning machine (KELM) models based on the GA, PSO, and MFO. Again, a fruit fly optimization (FOA) algorithm is used to optimize a KELM [31]. To avoid the limitation, an improved FOA is introduced by incorporating the Slime mould, Elite opposition-based learning, and levy flight algorithms. The proposed algorithm has a reliable trade-off between exploitation and exploration strategy. Double adaptive weight mechanism is introduced in the Moth flame optimization (MFO) to train kernel extreme learning machine (KELM). The proposed algorithm shows superior performance than other compared algorithms [32]. Table 1 illustrates the recent and compared algorithms that have been used for optimization in well placement optimization problem. It can be seen that a few well-known metaheuristics algorithms are used in experiments. Also, in many cases, experimenters only use primary algorithms for comparison purposes [24], [33]. Again, for better exploitation strategy, researchers incorporated local search techniques in the global search algorithms [13]. However, the local search algorithm's performance depends on the initialization [34]. Hence, investigators have implemented a hybrid algorithm incorporating non-traditional techniques for a better solution [35], [34]. The hybrid strategy, based on the best features of different algorithms, seeks a suitable solution to well placement optimization [36], [37], [17]. For instance, Dong et al. [38] proposed a hybrid of PSO to avoid the local optima, as the primary PSO algorithm can find a solution to a limited extent. Nwankwor et al. [24] used a combined HPSDE algorithm to determine optimal locations. They concluded that the hybridization of stand-alone DE and PSO algorithms performed better than stand-alone algorithms. Isebor et al. [22] combined two well-known search methods: the Mesh Adaptive Direct Search (MADS) and the PSO approach. Analysis demonstrates that the performance of the hybrid algorithm is superior compared to PSO and MADS. Humphries et al. [35] used a combination of PSO and generalized pattern search (GPS) strategy. Siddiqui et al. [39] conducted a comparison of CMA-ES, DE, and PSO in which DE performed better than PSO and CMA-ES.

B. RESEARCH GAPS AND MOTIVATIONS
Though researchers have mostly applied non-traditional [12], [16] and hybrid techniques [12]- [16] to resolve the well placement optimization problem, room for improvement remains. These techniques often fail to provide a better solution and faster convergence in different reservoirs [46]. Nevertheless, a better solution and faster convergence for a multimodal well placement optimization problem are still the dominant issues [5]. In the oil and gas industry, each reservoir has different sizes and search spaces for different well placement problems. Additionally, the surface may be non-smooth, or it may contain local optima. Metaheuristic techniques also require parameters tuning for different well placement problems. Hence, to provide better results in different search spaces, parameters tuning are required. However, the well placement optimization problems are computationally expensive. A single function evaluation requires one reservoir simulation, which is demanding in CPU time [47]. Thus, due to additional computational challenges, parameters tuning are difficult and researchers compare their work with few metaheuristic techniques [42].
In many studies only one reservoir is used. Hence, the performance of an algorithm is determined based on the results of one search space. In the oil and gas industry, each reservoir has different sizes and search spaces for different well placement problems. Additionally, the surface may be non-smooth, or it may contain local optima. However, in well placement optimization problem different reservoir will have different search space. The problem of this approach is that the algorithms parameter can be tuned to perform in one search space. Also, this process requires rigorous tuning of parameters. Again, the well placement optimization problems are computationally expensive. Therefore, in many sudies different reservoirs are not considered in the experiments [48]. For example, the PUNQ-S3 [4], [1], [5] and SPE-1 reservoirs [6] can be highly multimodal and both reservoirs are not considered in the same study. An ambiguity persisted when researchers used different reservoirs for evaluation in different studies. For example, DE performed better than CMA-ES in a specific study [39]. Conversely, CMA-ES performed better than DE in another study [49]. Hence, it can be observed that the algorithm's parameters setting in one study cannot be used in another study as it may not provide a better solution for a different reservoir [50]. Every reservoir requires a different exploration-exploitation strategy to explain this phenomenon [6]. Thus, the challenge is to obtain a better result in the different reservoirs using the same search algorithm.
To solve the limitation, this study considers an ensemble approach. Previous work demonstrates that quantumbased techniques, such as the quantum bat algorithm (QBA) and quantum particle swarm optimization algorithm (QPSO) performed better for well placement optimization [7], [51]. Moreover, quantum computation can manage highly nonlinear multimodal optimization problems [51], and PSO also works linearly. In contrast, the probabilistic approach can determine the QPSO's next position [52]. In QBA, researchers use the mean best approach to avoid local optima [53]. However, the QBA and QPSO techniques are better for specific reservoirs [6]. A single algorithm-based approach uses the same search approach. It may cause an algorithm to follow a similar trajectory. In turn, this may cause the algorithm to get stuck in local optima. However, the incorporation of different methods can improve the ability to find the best solution in complex conditions with complex areas [54]. Integrating several strategies using the appropriate adaptation mechanism allows an algorithm to select the appropriate strategy during optimization [55]. This integration can support search strategies with a variety of skills, improving the algorithm's performance. For example, a search strategy can find promising undiscovered areas. Using other search strategies can further improve the algorithm's performance [56].

C. RESEARCH CONTRIBUTIONS
A summary of the contributions of this study is as follows: • We combined different approaches with surrogate assistance which provides better solutions and faster convergence of each primary technique.
• A large set of algorithms are adopted for performance evaluation and two different reservoirs are considered for the evaluation.
• An ensemble approach of QPSO and QBA with surrogate assistance is proposed and implemented along with an approximation technique. It provides a better solution and faster convergence for the multimodal well placement optimization problems.
• Experimental evaluations were carried out to verify the proposed approach. The results demonstrate that the proposed method is more efficient and effective compared to existing optimization methods for well placement. To the best of our knowledge, this study is the first attempt that successfully applies ensemble-based optimization techniques to different complex reservoirs. The proposed search strategy provides the best solution in a dynamic search environment. The ensemble approach combines different methods and adjusts its strategy based on the success of its components. Furthermore, the ensemble approach does not require parameters tuning. Instead it utilizes the availability of diverse approaches at different stages and alleviates computationally intensive parameters tuning [57]. Finally, the ensemble strategy provides an effective tool to implement multiple search techniques suited to different reservoirs [57], [58].

II. PROBLEM FORMULATION
The Net Present Value (NPV) estimates the economic effect of a certain well location to extract oil/gas for a period. Certain well locations have effects on well NPV. We, therefore, propose using optimization techniques to find optimal well locations which provide maximum NPV. Figure 1 shows a common search technique to find the highest NPV for the well placement optimization problem. In this study, NPV is the objective function for well placement optimization. Equation (1) expresses NPV and considers oil, gas and water production, injection costs, oil sale prices, drilling cost, water production cost, and gas sale prices: (1) CAPEX designates the capital expenditure, Q w represents cumulative water production, P O signifies oil price, C W indicates the cost of produced water, Q O symbolizes cumulative oil production, OPEX stands for the operational expenditure, T denotes the numerical value of years that have passed, and D is the discount rate. The goal of well placement optimization is to maximize NPV and minimize expenditure. This research aims to optimize the location by maximizing production. In each iteration, the vectors containing all well positions in the PUNQ-S3 reservoir and SPE-1 reservoir are changed. For example, in the case study, investigators can place a well anywhere. After locating the well, the NPV or total production of the corresponding location is calculated. Therefore, an algorithm will try to change the position using the search technique. In each iteration, a new position is calculated and stored at its corresponding NPV. When the maximum number of iterations has been achieved, the algorithm displays the maximum NPV. The formulation of well placement optimization is the maximization of NPV based on well locations: Subjected to: where UB and LB represent the upper bound and the lower bound of the reservoir, respectively, NPV depicts net present value, and u n presents well coordinates.

III. PROPOSED SURROGATE ASSISTED QUANTUM-BEHAVED ALGORITHM
In the proposed approach an ensemble approach and proxy model are employed. The ensemble approach consists of QBA and QPSO. Additionally, a radial basis approximation technique is incorporated later. Figure 2 shows the concept of the proposed Surrogate Assisted Quantum-behaved Algorithm. It demonstrates the application of the QPSO and QBA techniques to the sample in the multimodal search space. Then the samples are evaluated and stored their corresponding fitness values. The sampling points and their corresponding fitness values allowed us to create a radial basis function (RBF) model. After solving the RBF model, its vertex has been identified. The approximate models are used to find better solutions rather than evaluate time-consuming cost functions. The following section gives an overview of QBA, QPSO, radial basis approximation, and the proposed approach. In the proposed methodology (refer to Figure 2), the QBA is used as a component of an ensemble approach with QPSO. Due to QBA's high exploration rate, it is used to evaluate the search space [6]. Yang et al. [59] proposed the original Bat Algorithm (BA), and they constructed the BA using three rules. The first rule states that usage of echolocation capability in every bat is similar, and echolocation capability can realize the distances between various background barriers and prey (food). In the second rule, bats in the x i position having velocity v i with varying wavelength λ 0 and fixed frequency f min use loudness A 0 to search for food. Depending on targeted proximity, adjustment of the rate of pulses and adjustment of the wavelengths in their emitted pulses is performed automatically. In rule three, they assumed that it could change loudness A 0 from a large positive value to a minimum value A min . The primary bat algorithm (BA) offers a fast convergence and straightforward implementation. However, the BA tends to get trapped into local optima points while optimizing the multimodal function. A study of bat trajectories reveals that, as the variety declines, many bats are restricted to the best local solutions. Also, the bats are guided by the best solution now available. However, if the best solution is categorized as a local point, the bats are then misguided. Furthermore, the BA has no mechanism for jumping out of local optima. Hence, to tackle the difficulties in boosting population variety and preventing premature convergence, quantum behavior is incorporated in the bat algorithm. The bats are guided by the present global solution in the early search stage, and the mean best position is employed during the later search. The formula, which is used to update location, is based on the Monte Carlo method [53]. Figure 3 depicts the flowchart of QBA. It indicates that frequency and pulse rates are updated after the random initialization of bats.
The upper and lower bounds are used to initialize the bat's position. The following equation determines the common solution: where X ij denotes the position of the j th dimension of the i th bat, X 0 and X m denote the upper and lower bounds, respectively, and rand is a random number between 0 and 1. This scenario leads to the following formula where we considered the bats' frequency, velocity, and position, respectively: where α represents a random vector ranging [0,1], f i represents the pulse frequency, f min is the minimum frequency, and f max is the maximum frequency. Furthermore, g t refers to the global best position of bats. In the following equations (9)(10)(11), the Doppler effect is considered. Moreover, for each bat, the compensating rate C is considered. As in normal air, the velocity of the air is 340 m/s, and the reformed equations (6)(7)(8) stand as: where f id represents the i th bats frequency at the dimension d, v t−1 g and v t g represent the velocity for the global best position at the (t − 1) th and the t th iteration, and C i refers to the number ranging [0,1] for the bat's position. ε is introduced so σ 2 , the standard deviation, remains positive. Furthermore, w stands for weight, x t id denotes the position in the d dimension for the i th bat at the t iteration, x t−1 id denotes the position in the d dimension for the i th bat at the t − 1 iteration, v t id denotes the velocity in the d dimension for the i th bat at the t iteration, v t−1 id denotes the velocity in the d dimension for the i th bat at the t-1 iteration, g t d denotes the position in the d dimension for the global best of the t iteration.
In QBA, the following equation can express a new position: where j(0,σ 2 ) denotes a gaussian distribution with zero mean and a standard deviation of σ 2 . x t+1 id is the i th bat's position in the d dimension at the t + 1 iteration, and g t d is the global best position in the d dimension at the t + 1 iteration. A t i is the i th bat's loudness.
Equation (12) shows the global best g t d is an attractant. Hence, the following equations express the position of the Quantum-behaved bat: where u is a random number. β is the contraction coefficient, mbest d is the average of personal best in the d dimension, and x t id is the i th bat's position in the d dimension for the t iteration.
After formalization of a new solution for every bat, we selected multiple solutions and used a random local nature walk. The new position for local search, therefore, was: where A t denotes the average loudness of bats, ε is used to denote a random number, x o is the present location, and x n is the new position after the local search.
In each iteration, the following equations can update the loudness A i and pulse rate r i : where A t+1 i denotes the i th bat's loudness in the (t + 1) th iteration and A t i denotes the i th bat's loudness in t th iteration. γ and are constant values. r 0 i denotes the i th bat's initial pulse rate, and r t+1 i denotes the i th bat's pulse rate at the (t + 1) th iteration.

B. QUANTUM PARTICLE SWARM OPTIMIZATION (QPSO)
In well placement optimization, different reservoirs have different properties. Thus, the search space will be different for each case [6]. To address this problem the QPSO is used in parallel with QBA as a component of an ensemble approach. Furthermore, the implementation of multiple search techniques is suited to different reservoirs [57].
Sun et al. [60] proposed an algorithm with the adaptation of the quantum mechanics principle for the basic PSO algorithm. There are certain dissimilarities between QPSO and PSO. PSO's current position is guided based on the personal and global best. On the other hand, QPSO follows a purely probabilistic scheme in which the next position is drawn from a probability distribution. In QPSO, current position is guided by mean best. Figure 4 illustrates the flowchart of QPSO.
The distinction between QPSO and traditional PSO is that QPSO follows quantum behavior in all particles, and all other versions of PSO follow classical Newtonian dynamics.
Instead of the position and the velocity, a wave function ( x, s) describes the particle's state in the Quantum-behaved Algorithm. The QPSO algorithm effectively removes the drawbacks and preserves the benefits PSO provides.
In PSO, the following equation updates the velocity of each particle: where x k i and V k i represent the i th individual's position and the velocity for iteration k, respectively. w is the weight vector, rand 1 and rand 2 are the random numbers, c 1 and c 2 are acceleration constants, pbest k i is the personal best of the individual i, and gbest k denotes the best position in the k iteration.
Each particle's new position is calculated using the following equation: In QPSO, for a population of k particles with d dimensions, . . , x id ) denotes the ith particle's location. Qi = (Q i1 , Q i2 , . . . , Q id ) denotes the i th particle's personal best, i.e., pbest. Similarly, Q g = (Q g1 , Q g2 , . . . , Q gd ) describes the global best position, i.e., gbest. q id , expressed as [64], denotes the local attractor of the i th particle on the d dimension: where ϕ is a random number. Sun et al. [60] proposed the mean best position (mbest) to avoid local optima. The mbest is calculated with the following equation (21): where mbest is the average position of all particles, and n is the number of particles. The following equation updates the i th particle's next position on the d dimension: where u is a random number and β represents the contraction coefficient and is expressed by [61]: where t max is the maximum number of iterations and t is the current number of iteration.

C. RADIAL BASIS FUNCTION APPROXIMATION
In optimization, an approximation technique is used to accelerate the search process for computationally expensive problems. The radial basis function (RBF) approximation technique is used to determine the optimal point based on all the particles' locations. After combining different approaches in the last stage, the approximation technique is employed to seek better solutions and faster convergence. In optimization, if N is the number of populations with d dimensions, then the input layer consists of the N × d matrix, and the output layer consists of the N ×1 matrix. Using RBF, we approximated the function u(x) as a linear combination of N radial functions.
The following equation expressed this [62]: where N denotes the data points number, λ j are the coefficients that need to be determined, and ∅ indicates the RBF. Thin Plate Spline (TPS) and Multi Quadrics (MQ) are considered advantageous for scattered data estimations [63]. For this reason, TPS is used in this study. The following equation defines a mth order TPS: ∅ x, x j = ∅(r j ) = r 2m j log r j , m = 1, 2, 3 . . . . . . , (25) where r j = x − x j denotes the Euclidean norm. Since ∅ is continuous, higher-order partial differential operators require a higher-order TPS. In the second-order equation, we utilized m = 2 as an assurance of the least C 2 continuity for u.

D. PROPOSED SURROGATE ASSISTED QUANTUM-BEHAVED ALGORITHMS
The key feature of the proposed approach is that it concurrently searches the solution space through two strategies, solutions, or individuals. Figure 5 and Algorithm 1 illustrate the flowchart and pseudocode of the proposed method. The proposed approach provides a framework for exchanging knowledge and immersive learning between algorithms with different search behaviors. Initially the population is subdivided into two groups. These two groups are used in two different search techniques such as QPSO and QBA to find a new position.
The new positions are evaluated and stored with the corresponding fitness values. Considering the position vector of the entire population as the input layer and the corresponding fitness values as the output layer, a proxy model with the TPS-RBF approximation technique is created. Then the optima of the approximate model are sought by utilizing the QBA technique. Finally, the optimal solutions of the approximate model are evaluated in the primary reservoir and the global best location is updated after comparing it with the current global best location. The approach of the proposed technique is below: Step 1: Subdivide the populations into two groups m and N-m.
Step 2: For the first m population, update search location X i,itr+1 using the QPSO algorithm.
Step 3: For the first m + 1 to nth population update, search location X i,itr+1 using the QBA algorithm.
Step 4: Evaluate the fitness function value.
Step 5: Apply TPS-RBF for the current population to generate an approximate model.
Step 6: Locate optima for the approximate model.
Step 7: Update personal best and global best location.

E. ADVANTAGES AND DISADVANTAGES OF PROPOSED TECHNIQUE
The characteristics of the proposed technique are: (1) the proposed approach spontaneously subdivides its population into two groups, enabling it to perform better than existing algorithms to address a nonlinear, multimodal optimization problem, (2) the PSO, BA, and GA have the disadvantage of premature convergence, (3) the proposed approach overcomes this limitation because it does not update its location based on the personal best information, and there is no explicit global best either, and (4) it functions as a QBA and QPSO, giving this technique the advantages of these two algorithms [64].

IV. RESULTS AND DISCUSSION
In this study, Eclipse, a numerical simulator, and MATLAB were used. The simulator provided production data for specific well placement. All simulations ran on a PC with an i7-7500U CPU @2.70 GHz (4 CPUs), 3.2 GHz, and 8 GB RAM. Two different case studies were considered. Case study 1 used the SPE-1 reservoir model, and case study 2 used the PUNQ-S3 reservoir model. In case study 1, the number of iterations and particles were 100 and 20, respectively, for all algorithms. In case study 2, the number of iterations and particles were 30 and 5, respectively, for all algorithms [6]. Table 2 lists the parameters of algorithms. Table 3 lists the economic parameters used to evaluate the objective function, as depicted in Equation 1.

A. CASE STUDY 1
In this case study, a 3D simulation of a black oil reservoir is used to develop the SPE-1 model. [65,66] describes detailed properties and specifications of the reservoir model, as shown in Figure 6(a). The SPE-1 model has a 10 × 10×3 grid block. As (x, y) is the coordinate for the wells, this case study optimizes the variable 2 × 2 for two wells. The dataset of this reservoir can be found at: https://www.spe.org/web/csp/datasets/set01.htm.

B. CASE STUDY 2
In this case study, the PUNQ-S3 model is developed by utilizing a real field used by Elf Exploration Production to test methods for quantifying uncertainty assessments. The PUNQ-S3 has 19×28×5 grid blocks. The details of the reservoir model can be found in [67]. Four vertical wells for optimization are considered. Hence, this experiment optimizes the 2 × 4 variable. Figure 6(b) shows a detailed description of Case Study 2. The data set of this reservoir can be found at: https://www.imperial.ac.uk/earth-science/research/researchgroups/perm/standard-models/eclipse-dataset/.

C. PERFORMANCE CRITERIA
Clerc [68] revealed that in trial a mean value runs alone, and it might be inadequate to measure the performance of an algorithm. The researchers utilized the graphical representation of the convergence curve with average value versus function evaluations. The standard deviation provides the consistency of the algorithm. As the evaluation process of the algorithms is a prime concern for this task, the researchers considered several criteria [68], [69]. These criteria are below: Effectiveness is a simple, important measure of performance. This is a measure of the average value between tests of the best solution found as a percentage of the global optimum or,f where f (p) refers to the solution of p, p * denotes the global optimum solution, p ∧ i represents the best solution found in trial i in N number of trials for each algorithm.   Efficiency, another crucial criterion, indicates the speed at which the algorithm reaches a performance level utilizing a unique evaluations number required to find a proper solution, at least 98% of the best solution found, on average FIGURE 6. Initial pressure, oil, and gas saturation properties under different case studies [6].  between tests or,L where L 98 i refers to the unique function evaluations number that is essential to calculate q as f (q) > 0.98f (p * ) for trial i (for maximization) and M denotes the function evaluations gross number per trial.

D. CONVERGENCE ANALYSIS
Each algorithm is run 30 times and their average convergence curve is shown in Figure 7. Figure 7(a) shows that, in case study 1, the proposed technique provided superior results compared to other compared algorithms for finding better NPV. Additionally, this study established that the second-best algorithm is QPSO, and QBA achieved the third best NPV. GA and PSO were trapped in local optima. Case study 1 showed the proposed algorithm has faster convergence and achieved the highest NPV. Figure 7(b) shows that the proposed technique acquired a better NPV in the PUNQ-S3 reservoir model. QBA was the second-best algorithm. However, The GA, CSA, and GSA algorithms could not provide a satisfactory NPV in either case study. It has been observed that the performance of the stand-alone algorithm was inconsistent [6]. Additionally, Figure 7 shows that the proposed algorithm reached convergence after 10 and 60 iterations. However, other algorithms required more iterations to achieve convergence. Therefore, it is worth mentioning that the proposed approach provided faster convergence than other algorithms in both case studies. Furthermore,  Table 4 lists the minimum, maximum, average, standard deviation, efficiency, and effectiveness of NPV from trials. It also shows that the proposed algorithm is better in four criteria. The BA, however, had superior efficiency. The reason for this phenomenon is that efficiency is calculated concerning its own best solution. Despite having a lower NPV than other algorithms, the BA achieved higher efficiency.
The QBA provided the maximum value, but the standard deviation was higher. Furthermore, the proposed technique had the highest average with the lowest standard deviation. Table 5 demonstrates that the proposed algorithm is better in five criteria. Figure 8 shows the convergence curve of the best performed result. Figure 9 shows the proposed algorithm provided the 2 nd lowest standard deviation compared to other algorithms. It can be inferred that the proposed algorithm's performance is better than in both of the case studies.

V. LIMITATIONS
Although researchers contributed to many areas in well placement selection, improving the reservoir proxy model, and NPV of well placement are the main interests of this study. To enhance results and extract maximum NPV from input data, a new type of algorithm is employed in this study. This study attempts to optimize well positions. In the oil and gas sectors, however, investigators must optimize history matching and well management parameters. In this analysis, well controls remain fixed. Nevertheless, the need exists for optimal well controls [72]. Furthermore, deciding the location of oil wells and operating settings (for example, infusion/recuperation rates for heterogeneous supplies) poses difficult challenges and has an impact on underground oil recovery and monetary value. The optimization of well position is an integer-based problem.
Moreover, researchers usually optimize the well position first, and they optimize well control settings with the fixed VOLUME 10, 2022   optimal well location [73]. Optimization requires a thorough sensitivity analysis. Additionally, only two case studies are used in this study. Furthermore, uncertainty analysis is not considered in this study. The Monte Carlo Simulation (MCS) is the most often used method for dealing with uncertainty problems. The key drawback of this method is that it converges slowly, which means it is expensive to compute. To maintain a variety of uncertainties, a standard MCS needs a few hundred runs, which is impractical for very large and complicated models. Furthermore, the findings of a MCS are highly vulnerable to distribution assumptions. Even if the mean and variance are the same, the outcomes can vary significantly due to different distributions [74]. To address these shortcomings, future research can integrate rigorous architecture optimization based on polynomial chaotic expansion [75] and Sparse Grid-based Polynomial Chaos (SGPC) [76]. Investigators can also consider info-gap decision theory as an alternative to MCS.

VI. CONCLUSION AND FUTURE DIRECTIONS
In this study, the QBA algorithm and QPSO were implemented in parallel for the investigation of well placement optimization. The performance is investigated on two separate reservoirs. As a standalone technique, QPSO's performance was better in the PUNQ-S3 reservoir than other stand-alone techniques. The QBA's performance was also better than other stand-alone techniques for the SPE-1 reservoir. Hence, this study implemented a Surrogate Assisted Quantum-behaved Algorithm and exploited the different search techniques.
The experimental results show the proposed technique can enhance the search technique and provide a better solution than other algorithms. Concluding remarks are as follows: • Due to the same search pattern, a stand-alone search algorithm cannot perform well.
• Quantum-based metaheuristic techniques are less likely to be stuck in local optima and less susceptible to premature convergence for well placement optimization.
• In both case studies, the proposed approach's standard deviation is lower than other existing state-of-the-art algorithms. Hence, the proposed approach can provide better solutions for well placement optimization problem.
• QBA performed well in case study 2 and QPSO performed well in case study 1.
• The ensemble strategy effectively solved the well placement optimization problem by providing better results for both case studies. The conclusion of the study is that proxy modelbased optimization techniques provided better results.
Furthermore, ensemble approaches of algorithms effectively address dynamic search space. Future work on well placement optimization can utilize other methods, such as Elephant Herding, Monarch Butterfly, the Kidneyinspired Algorithm, the Pity Beetle Algorithm, the Spotted Hyena Optimizer, Thermal Exchange Optimization, the Grasshopper Optimization Algorithm, and the Grey Wolf Optimizer.
Moreover, the main focus was the optimization algorithm. The history matching is not considered in this study since it is in another stage of oil production. As part of ongoing research, classical benchmark functions and other reservoirs will be included to evaluate the performance of the proposed algorithm. A sophisticated deep learning and data mining method to address reservoir uncertainty modeling is the clear way of the future [77].
The authors identified four challenges for well placement optimization: i) The general closed-loop workflow of proxy simulation for uncertainty optimization. ii) An approximation approach to generate a surrogate model that is computationally feasible and quantifies the uncertainty set of all chosen models. iii) Carrying out optimizations while considering complexity. iv) Perform risk identification and decision-making under decision-makers' attitudes and expectations. For metaheuristic algorithms, investigation of the effects of decision variables, limits, and internal parameters is critical. Problem-related customization, such as algorithm parameter tuning, is also important, and generating a diverse population to prevent local optima is a challenge. Strong diversity can help to prevent problems caused by local optima. To escape local optima, future research should consider incorporating levy flight, chaotic maps, and other techniques into current metaheuristic algorithms. Finally, researchers should consider a large search space to find the optimal solution. VOLUME