Artificial intelligence for well placement optimization

Abstract

A method of optimizing a set of oil producing wells to be drilled in a reservoir includes receiving a base reservoir simulation representative of a geological oil reservoir, clustering the reservoir, and generating, in a first process, a possible pool of well positions in the reservoir, where each well position satisfies plausibility and respects operational constraints. The method further includes, in a second process, optimizing the pool of well positions, the second process maximizing a dynamic variable of the reservoir, and outputting the optimized pool of wells to a user. In embodiments, the second process further includes training a DNN surrogate model to forecast the production of the reservoir; and using the trained DNN surrogate model, optimizing the pool of wells.

Description

[0001] Patent Application Under the Patent Cooperation Treaty (“PCT”) ARTIFICIAL INTELLIGENCE FOR WELL PLACEMENT OPTIMIZATION CROSS-REFERENCE TO OTHER APPLICATIONS This application clams the benefit of United States Provisional Patent Application No. 63 / 543,614, filed on October 11, 2023, for “ARTIFICIAL INTELLIGENCE FOR WELL PLACEMENT OPTIMIZATION,” the entire disclosure of which is hereby incorporated by reference. TECHNICAL FIELD The present invention is directed to the application of artificial intelligence techniques to the optimal placement of oil wells in a reservoir. In particular, the present invention utilizes two loop optimizations. A first one to optimize the definition of possible well locations from an initial screening, and a second loop optimization to optimize final locations from a given set of possible well locations. BACKGROUND OF THE INVENTION Most of the worlds easily developable crude oil and natural gas reservoirs have likely been discovered and the annual volume of crude oil and natural gas discoveries is declining. The IEA forecast global demand to be 3.2 million barrels per day higher in 2030 than in 2023. Although 50,000 crude oil and natural gas reservoirs have been discovered only ten percent of these oil reservoirs have a material impact on global supply. Therefore, it is important for operators to optimize their field development plans, specifically where to locate wells in the reservoir to maximize production, improve capital expenditure (CAPEX) and reduce operating expenses (OPEX). This responsibility involves multidisiplinary teams of reservoir engineers (RE) and geoscientists attempting to develop optimal reservoir models that place wells based on complicated, heterogeneous geological maps and systems that depends on physical equations to characterize the interaction of the rock and the fluids. Traditionally, this problem has been approached by manual update of the model over multiple interations, where REs locate the best wells based on their assembled knowledge and experience combined with forecasting tools to evaluate the solution. However, due to the vast space to explore and limited timeline to make investment decisions, these techniques present many limitations. For this reason, leveraging state-of- the art computation hardware, algorithmic optimization and mathematical techniques can be used to solve this problem. With respect to optimization of a field development plan (“FDP”) for a hydrocarbon reservoir, it should be noted that, in general, all optimization algorithms require a first stage where the search space needs to be explored globally. After this first exploration stage, the solution found is often incrementally improved until some optimization stopping criterion is satisfied. The selection of an optimization algorithm has a strong dependency on the type of problem to solve, the size of the search space (i.e., the number of variables to optimize), its robustness to introduce constraints of different types, and the search / exploration capabilities. However, conventional optimization techniques for well placement present several common drawbacks: First, optimizers depend on computationally expensive finite volume reservoir simulators, where a simulation for large models can take hours or days. Additionally, optimizers potentially need to call the reservoir simulator millions of times to reach a local / global optimum. The combination of these issues results in most optimizers being impracticle for well placement, or, even if they can be used, the solution is often sub-optimal. Second, to be useful for well placement optimization (“WPO”), an algorithm must be created that generates well locations that are geologically plausible, in other words, respecting all geological constrains such as, for example, drilling in a high productivity rock type. In other words, the optimization algorithms must provide solutions that drilling engineers and geoscientists consider realistic enough to be implemented directly in the reservoir. Many optimizations – although theoretical - are not realistic, because they do not respect geological constraints. Finally, another challenge facing current optimization techniques is the mangement of realistic operational constraints such as minimum and maximum distances between wells, maximum rate to injection, etc. From a mathematical point of view, these non-linear constraints are complicated to manage in the optimizer and affect the convergence of the problem. This is due mainly to multiple local optima, and computational intensity, described below. One of the main challenges of nonlinear constrains is that they may have multiple local optima, which are not necessarily the global optimum. This makes it difficult to guarantee that the solution found is the best possible one. For example, when trying to find the lowest point on a complex surface, an algorithm may get stuck in a local valley that is not the deepest one globally. Moreover, nonlinear programming (NLP) can be computationally intensive, especially for large-scale problems. The algorithms required to solve NLP problems often involve complex calculations and iterations compared to linear programming. As detailed below, what is needed in the art is a reliable methodology that overcomes these problems with conventional optimizers, and that can determine an optimum field development plan, specifically optimizing the well locations while honoring the constraints to mazimize investment return. SUMMARY OF THE INVENTION In embodiments, systems and methods for the optimization of field development plans, specifically for well locations in real reservoirs, are presented. In embodiments, a method includes two loop optimization techniques. In a first loop, an optimization is made from an initial screening to define possible well locations. In a second loop, final well locations are optimized from a given set of possible well locations, where a dynamic variable, such as, for example, Net Present Value (NPV) is maximized, and all operational constraints respected. BRIEF DESCRIPTION OF THE DRAWINGS Fig.1 illustrates well trajectory parameterization, in accordance with various embodiments. Fig.2 illustrates projection of a solution into a feasible region, in accordance with various embodiments. Fig.3 is an example process flow chart for an example well optimization method, in accordance with various embodiments. Fig.4 is an example process flow chart for an example smart sampling sub-process, in accordance with various embodiments. Fig.5 is an example system diagram for surrogate modelling, in accordance with various embodiments. Fig.5A is an example process flow chart for training a deep neural network, in accordance with various embodiments. Fig.6 illustrates a block diagram of a computer device suitable for practicing the present disclosure, in accordance with various embodiments. Fig.7 illustrates an example computer-readable storage medium having instructions configured to implement software implementations of and / or practice various aspects of the present disclosure. DETAILED DESCRIPTION OF THE INVENTION For reservoir engineers (RE), optimizing well locations is not an easy task. Due to the large number of variables to optimize (e.g. dozens of wells, duration of water and gas cycle), the number of possible scenarios could be greater than 1024, which is even more than the number of stars in our galaxy. In addition, the RE has to wait hours -- and sometimes days -- to run a complex reservoir simulation for a single FDP scenario. After that, the RE has to analyze and manually determine the next scenario (which is not easy due to the complexity of the objective function to optimize) and follow this workflow until they reach a target deadline (e.g., six months). Reservoir simulation in porous media is a sophisticated computational technique used to model and predict fluid flow behavior in subsurface reservoirs. It combines principles of physics, mathematics, and computer science to create digital representations of oil and gas reservoirs, allowing engineers to analyze and forecast fluid movement through porous rock formations over time. The process involves integrating geological, petrophysical, and fluid models into a dynamic simulation that solves complex flow equations based on Darcy's law and conservation principles. These simulations typically use finite difference methods to discretize the reservoir into a grid of cells, enabling detailed 3D modeling of fluid interactions. Reservoir simulators can be black-oil, compositional, or thermal, depending on the level of detail required in fluid behavior representation. This powerful tool is essential for optimizing FDPs, making informed operational decisions, predicting reservoir performance, and evaluating enhanced oil recovery techniques, ultimately improving resource extraction efficiency and reducing costs in the oil and gas industry. The main idea is as follows. Imagine there are 100 possible well locations, and it is necessary to select 20 wells. The total number of possible combinations of the desired 20 wells is 5.359833704 x 10^20. Additionally, due to reservoirs usually being highly heterogeneous, with significant variations in porosity and permeability, it is difficult to determine the optimal combination (considering interference between wells) without a reservoir simulator tool. Consequently, humans have to manually run scenarios, with each scenario taking minutes or hours. It is very challenging for humans to perform this optimization manually. Usually, this convneitonal methodology is sub-optimal because the number of scenarios that the RE is able to evaluate is small, and the solutions are biased, which does not guarantee the improvement of the FDP iteration by iteration. For that reason, well location optimization is highly challenging to solve manually, and computational / mathematical techniques are necessary to address this problem. Decision-making for production strategy has been addressed formally from a mathematical simulation-based optimization perspective,

[0018] . For example, one paper applied optimization using Genetic Algorithms (GA) to China’s Pubei oil field

[0019] . Another example optimized the Eileen West End Areain Greater Prudhoe Bay, operated by British Petroleum, using for EOR such as WAG

[0014] , where the key parameters evaluated were injection volume, injection rate, WAG ratios and WAG sequencing or WAG cycle number. Although the formulation for well location is relatively new, several optimization techniques can be adapted to solve the target problem. Nonetheless, it should be understood that, in general, all optimization algorithms require a first stage where the search space needs to be explored globally. After this first exploration stage, the solution found is often incrementally improved upon until some optimization stopping criterion is satisfied such as maximum number of simulations, maximum computational time, no improvement of the objective function, etc. The selection of the algorithm has a strong dependency on the type of problem being solved, the size of the search space (i.e., the number of variables to optimize), its robustness to introduce constraints of different types, and search / exploration capabilities. In what follows, it is proposed to solve these drawbacks with a new optimization workflow based on artificial intelligence. In embodiments, the workflow manages the drawbacks separately. First, a two-loop AI optimization workflow may be used, where the first workflow provides a list of realistic candidate well locations (which may include, for example, hundreds of wells), and the second, an optimization workflow optimizes the best possible locations given a well candidate and a range of target wells (e.g., dozens). Secondly, in embodiments, a new surrogate modeling technology may be used that enables running simulations orders of magnitude faster than traditional reservoir simulation techniques -- while conserving the same level of accuracy. Finally, a new formulation may be used to manage complex and non-linear constraints in the AI optimizer algorithms. Methodology Optimizing well locations from a simulation-based optimization perspective is very challenging due to the large number of possible well locations to study, potentially greater than 1050possible combinations, as well as the dependency on an expensive reservoir simulator that limits the number of possible scenarios that can be evaluated. In embodiments, these drawbacks may be resolved using two different workflows that run sequentially, as follows: 1. Optimization to obtain possible well locations (high-loop optimization). Usually, well optimization algorithms have to evaluate possible well locations that are known, based on the properties of the reservoir, to not provide good results. Thus, in embodiments, the original action spaces may be pruned to obtain a set of solutions that could be optimal, which may then be further explored in a fine optimization workflow. In embodiments, to carry out this initial screening, an optimization problem may be used that provides a set of well trajectories to optimize in a reservoir-simulation (low-loop) first optimization workflow. 2. Actual well location optimization (low loop optimization). Once an initial well-trajectory candidate is generated, the goal is to then obtain the final optimal well-trajectory that maximizes a dynamic property, such as, for example, Net Present Value (NPV), or, for example, recovery factor or oil production. In embodiments, this problem may be formally addressed from a mathematical simulation-based optimization perspective. However, such a perspective depends on a highly time-consuming commercial reservoir simulator (such as, for example, INTERSECT, ECHELON, T-NAVIGATOR, or the like) that limits the exploration to a small number of scenarios. For this reason, in embodiments, surrogate modeling (described in detail below) may be implemented. Surrogate modeling based on deep learning has been shown to yield incredible results, allowing reservoir simulations to be accelerated several orders of magnitude

[0013] , [9]. However, known surrogate models present several limitations due to the fact that they are often not very accurate for the prediction of oil production, and they do not introduce physical elements. In this disclosure, this problem is solved, by using a new deep learning architecture that is inspired by the physcis of reservoir simulation, and which can also be trained based on the residual of the equation. It is noted that residuals in differential equations measure how well an approximate solution satisfies the original equation by quantifying the error when the approximation is substituted back into the equation. In addition, in embodiments, the disclosed methodology may implement a mechanism to manage non-linear operational constraints, which is a key element in optimizing realistic field development plan scenarios. Both workflows may become more challenging when introducing geological uncertainty. Reservoir simulations require input from geological models (e.g., maps of porosity, permeability, etc.). These maps are typically generated using geostatistics and limited input data, meaning that information from a small number of wells is used to create the geological model. Consequently, there are a set of equiprobable geological models that can capture the reservoir's behavior. From a mathematical perspective for well location optimization, uncertainty is introduced by calculating the expected value of an objective function, such as, for example, Net Present Value (NPV). This expected value is computed by simulating every production strategy that our optimization algorithm evaluates for all possible realizations, and subsequently calculating the average. This technique is computationally intensive, inasmuch as running simulations for multiple realizations can be expensive, especially for complex reservoir models. To address this drawback, in embodiments a new algorithm is implemented that the inventors term "smart sampling." The idea behind smart sampling is to determine the smallest set of realizations that allows for a similar expected value of the objective function to be obtained compared to the full set of realizations representing the uncertainty of the geological model. Each of the two workflows are described next. 1. Optimization to Obtain Possible Well Locations The main objective of a first high-loop optimization is to carry out an initial screening to obtain a set of candidate wells that are both plausible and that respect operational constraints. It is noted that the high-loop optimization goals are to provide an initial list of optimal well positions as input data for the lower loop optimization that is responsible for finding the final optimal well locations based on a reservoir-simulation optimization perspective. In embodiments, a first task is to automatically obtain a set of possible well-location trajectories. In embodiments, a well-location trajectory in 3D, y, may be parameterized using six different variables, y ∈ R6. These include three spatial coordinates for the initial position of the well (x, y, z), and then a well length, l, an elevation from the horizontal α, and an azimuth β (angle from true North), for the well, all as shown in Fig.1. However, such a formulation presents several limitations: First, a well trajectory in a reservoir simulator and / or surrogate must be provided in (i, j, k) cell coordinates which are integer variables. Second, we need to transform this straight-line formulation into trajectories that follow the irregular shape of the layers that usually characterize a pilar gridding mesh. These challenges are approached using a new function φ : R6→ I3, which may be used in embodiments to obtain ˆy from y; ˆy = φ(y). (For clarity, is intended to be the standard “y with a hat over it”; however, due to typographical limitations, in this disclosure the hat generally precedes the variable, and the same usage occurs for a “y bar” variable, where the bar preceds the variable letter, as “￣y”). This function φ projects x in cell coordinates space into the cell coordinates grid and provides a realistic trajectory that is aligned with the pillar gridding mesh that was given. Once a mathematical parametrization for a realistic well trajectory is obtained, in embodiments, the next task is to obtain an optimal set of candidate wells. In embodiments, this may be addressed from a mathematical optimization perspective. Here a geological / dynamic property may be maximized that provides a realistic and optimal screening for potential well candidates. In embodiments, for example, the normalized opportunity index (“OI”) may be used as this property. An expression of this is found, for example, in

[0017] . The OI provides an initial value for the amount of oil, and the mobility ratio associated with a well trajectory. Therefore, in embodiments, the optimization problem may be formulated as: Max OI (φ(x, y, z, l, α, β)) Eq. (1) x,y,z,l,α,β such that: x ≤ x ≤￣x and lower limit for each variable is given based on reservoir knowledge, x, y, z, l, α, β ∈ RNw,1and Nw is the number of wells which is given. However, it is noted that this formulation could collapse in a single trajectory and provide a set of wells that has low variability. For this reason, in embodiments, the loss function may be modified, introducing a new term. This new term may, in embodiments, be composed of two different components: first, the distance between the different wells in ˆy is calculated, dist(ˆy) ∈RNw,Nw, and subsequently the standard deviation of upper triangular distance matrix is calculated. Finally, Equation (1) may be modified to read as follows: Max λ1 ∗ OI(φ(x, y, z, l, α, β)) + λ2 ∗ std(distUT(φ(x, y, z, l, α, β))) x,y,z,l,α,β Eq. (2) such that: x ≤ x ≤￣x y ≤ y ≤￣y z ≤ z ≤￣z where λ1 and λ2 are the weights of the objective function, OI is normalized by dividing by the total distance of the well, and std is normalized using the following equation (std − stdmin) / (stdmax− stdmin). It is important to note that, in alternate embodiments, seismic attributes, for example, may be used instead of the OI. To do that one replaces the OI in Equation (2) with the seismic attribute map, such as impedance. However, in real fields the number of wells, lengths, azimuth, etc. may be different depending on the area of the reservoir, which may have different geological and / or petrophysical properties. To account for this, in embodiments, the optimization problem expressed in Equation (2) may be split into sub-optimization problems for each of these regions, obtaining the following Equation (3): max(xi, yi, zi, li, αi, βi)λ1 ∗ OI (φ(xi, yi, zi, li, αi, βi)) + λ2 ∗ std(distUT(φ(xI, yi, zi, li, αi, βi)))∀ i, i = 1, ...,Nc Eq. (3) such that: xi ≤ xi ≤￣xi ∀ i, i = 1, ...,Nc xi, yi, zi, li, αi, βi∈ RNwi,1and Nwiis the number of wells for the cluster i. In addition, uncertainty could be introduced in this problem, replacing the value of the objective function in Equations (1), (2) and (3) with the expected value of the objective function, and calculating the expected value such as was described above. In this case, smart sampling is not necessary; this is because the evaluation of the objective function may be performed very rapidly. In embodiments, different algorithms may be used to compute the clusters of the different sub-regions of the reservoir for Equation (3). For example, in some embodiments k-means clustering may be used. K-means clustering is a method of vector quantization, coming originally from signal processing, that aims to partition n observations into k clusters, in which each observation belongs to the cluster with the nearest mean (cluster centers or cluster centroid), thereby serving as a prototype of the cluster

[0010] . In embodiments, x and y coordinates as well as depth may be used in order to calculate clustering with k-means. However, k-means presents clustering with salt-and-pepper noise [1], which could be a drawback for the optimization problem of the Equation (3). For that reason, in embodiments, post-processing such as median filters to remove this noise may be used. It is important to note that alternatively, other post-processing techniques could be used, such as, for example, a pre-trained DNN that allows for the removal of noise

[0020] . Alternatively, other algorithms for clustering may also be used such as, for example, deep learning for image segmentation. In embodiments, to determine if a cluster may be used as a possible area for well location, a tolerance tol may be defined. This tolerance may be compared with the average opportunity index for each cluster to determine if wells should be drilled in this area. Finally, Algorithm 1, provided below, provides a complete method for obtaining the initial set of optimal well locations. Moreover, it is noted that a free derivative algorithm may preferably be used to solve the optimization problem of Equation (3). For example, the covariance matrix adaptation evolution strategy (CMA-ES) [6] algorithm may be used, which has shown high-quality results obtained in similar problems

[0012] . Additionally, the number of iterations necessary to solve CMA-ES is higher than 1MM, and thus, in embodiments, high- performance computation (HPC) may preferably be used, as well as a thoroughly optimized implementation of the function φ. Algorithm 1 - Algorithm to optimize possible well locations Given a reservoir realization m; Given number of clusters Nc; Calculate 3D cluster C using k-means, C ← kmeans(m); Given C, user provides constraints of Equation (3) s well as number of wells per cluster, Nwi; Given tolerance C, tol Given λ1 and λ2 x, y, z, α, β, l ← ∅ for i = 1 to Ncdo if average(Ci) > tol then Calculate xi, yi, zi, αi, βi, lisolving Equation (3) by CMA-ES x, y, z, α, β, l ← xi, yi, zi, αi, βi, li end if end for 2. Second Workflow: Optimization for Well Locations In embodiments, the objective of the second workflow is, given a set of well trajectories defined by W ∈ RNw= ¯NP+ ¯NIand characterized by x, y, z, l, α, β, a subset of producers and injectors p, q ⊂ W is obtained such that a dynamic variable (e.g., NPV) is maximized. It is noted that in other embodiments, another objective function may be used, instead of NPV, as noted above. In embodiments, this problem may also be addressed from a mathematical simulation-based optimization perspective. In particular, the problem may be formulated using the following Equation (4): max NPV (p, q) Eq. (4) p,q such that: dpp≤ d(i, j) ≤ d￣pp, ∀ i ∈ P, ∀ j ∈ P − {i} dip ≤ d(i, j) ≤ d￣ip, ∀ i ∈ I, ∀ j ∈ P dii≤ d(i, j) ≤ d￣ii, ∀ i ∈ I, ∀ j ∈ I − {i} sum(x) ≤￣Np-i ∀ i = 1, ... , Nc sum(y) ≤￣N I-i ∀ i = 1, ... , Nc where p ∈ {0, 1}¯NPis the binary vector for producer wells, and q ∈ {0, 1} ¯NI is the binary vector for injector wells; d (., .) is the Euclidean distance between two wells;￣NPis the maximum number of producer wells to drill, and￣NIis the maximum number of injectors to drill. In addition, I is the set of injector wells for a given y, and P is the set of producer wells for a given x. Finally, in embodiments, NPV may, for example, be calculated by: NPV = SUM from t=1 to Nt[[Q0 to tn(p,q)∗ Po− Nw∗ Pw] / [(1+λ)tn / 365]*[Δtn]] Eq. (5) (other calculations may also be utilized). Here, Qoindicates the field-wide oil production rate in stock tank barrels per day (STB / d), and Nw indicates the total number of wells drilled. Oil price, the cost to drill a well, and the annual discount rate are denoted by Po, Pw, and λ respectively. The quantities tnand Δtnindicate time and the size of the time step in days, respectively. One of the main challenges of this problem is to calculate Qo(p, q). Usually, Qois calculated using a computationally expensive reservoir simulator. However, this poses an important drawback, as optimization techniques to efficiently solve Equation (4) require hundreds of thousands of runs of the simulator, which is computationally prohibitive. In order to address this drawback, in embodiments, the use of a new surrogate modeling approach based on deep learning may be used. The example surrogate modeling approach can run simulations at a rate that is orders of magnitudes faster than traditional simulators, but still obtain a similar level of accuracy for well rates as that of traditional simulators. In embodiments, surrogate modeling offers a significant advantage with respect to traditional simulators, in that we can calculate the gradients or Jacobian matrix of the outputs (e.g., Qo) with respect to the input data (p, q) using automatic differentiation (AD) tools that are included in many machine learning libraries, such as, for example, Pytorch or TensorFlow. This allows one to use different gradient optimization algorithms, such as, for example, the Broyden-Fletcher-Goldfarb-Shanno algorithm (BFGS) to solve Equation (4). However, gradient optimization presents limitations in managing binary variables and non- linear constraints, as the derivative for binary variables can trap the optimization algorithm in local minima. Thus, in embodiments, for this specific formulation, free derivative optimization may be used to solve Equation (4). To do this, CMA-ES may be used, in embodiments. However, one of the main challenges in CMA-ES optimization is the introduction of non-linear operational constraints that are found in Equation (4). Although it is possible to introduce penalty constraints or an augmented Lagrangian method ([3],

[0011] ) here, this is challenging because the majority of the time the constrained term could direct the search in the optimization to suboptimal solutions. For that reason, in embodiments, a new operator for CMA-ES may be implemented that is applied after mutation and crossover. The main idea of such an operator is to project a solution that is unfeasible (a solution that does not respect the constraints) to the closest feasible solution in a feasible region. This is illustrated, for example, in Fig.2, where an unfeasible solution 201 is projected into the feasible region 230, at a minimum distance. This results in projected solution 210, which is in the feasible region 230, and which is the closest solution to the unfeasible solution 201 that satisfies or respects the constraints. Algorithm 2, provided below, illustrates how these projections may be carried out to respect the various distance constraints between two wells, in accordance with various embodiments. In general, these constraints are created or imposed by reservoir engineers or geo-scientists. Algorithm 2 Projection operator algorithm Given a solution x Given a feasible region characterized by the constraints of Equation (4) for i ∈ D do for j ∈ D − {i} do if dpp≤ d(i, j) ≤ dppthen Calculate the closest well to i that is feasible for ∀ j ∈ D − {i} → r Replace i by r end if end for end for where D is the set of drilled wells D ⊆ W. The inventors have observed that the projection operator provides better solutions the traditional penalty functions. Thus, penalty functions are generally not recommended as the primary optimization technique due to several drawbacks. They can cause numerical instability by creating severe slope changes at constraint boundaries, making it difficult for algorithms to converge. Further, their effectiveness is highly sensitive to the choice of penalty parameter, often requiring an iterative approach that increases computational cost. Penalty methods also risk getting trapped in local extrema, necessitating multiple optimization runs from different starting points. Additionally, they may not be equally effective for all types of constraints. Due to these limitations, practitioners often prefer alternative approaches such as barrier methods, exact penalty methods, interior point methods, or Lagrangian methods, which can provide more stable and efficient optimization processes. However, in embodiments, these drawbacks may be solved with the disclosed projection operator, as described above. Fig.3 is an example process flow chart for a method of optimizing an FDP, in accordance with various embodiments. With reference to Fig.3, process flow begins at block 301, where a base reservoir simulation case is uploaded. This may be, for example, a single case, or, for example, a set of cases with different maps of porosity and permeability that characterize the geological uncertainty. From block 301, process flow may optionally continue to block 302 or may skip block 302 and proceed directly to block 303. If the optional path to block 302 is chosen, then a smart sampling sub-process is performed. The objective of the smart sampling here is to obtain a subset of K realizations that provide identical or similar probability distributions to those of the N realizations, thus simplifying the computations. Details of this sub-process are described below with reference to Fig.4 From block 302, if that path is taken, process flow moves to block 303, where the reservoir is clustered. Alternatively, as noted, process flow may move directly from block 301 to block 303, and not perform smart sampling. Block 303 represents the first workflow described above, with reference to Equations (1), (2) and (3) and Algorithm 1. The output of this clustering operation is a pool of wells, generated at block 304. From block 304, there is another optional choice. This is to either proceed to blocks 305 followed by block 306, to train a surrogate model, or to proceed directly to block 307, without surrogate model training, and proceed directly to block 307, where the FDP is optimized. Block 307 invloves the second workflow described above, with reference to Equations (4) and (5) and Algorithm 2. If blocks 305 and 306 are not implemented, then instead of training a surrogate model, a reservoir simulator or other forecasting took may be used, to predict oil / gas / water production / injection. However, if surrogate training is used, then with reference to block 305, in embodiments, any sampling technique may be used to generate a sample of cases to train a surrogate model. For example, in some embodiments, Latin Hypercube Sampling (“LHS”) may be used to obtain the sample of cases. It is noted that the sampling must respect operational constraints such as, for example, minimum / maximum distance between wells, and / or maximum / minimum production or injection. In embodiments, what is sampled is, for example, well location. In either option, process flow ends at block 307. For the case where an examplary embodiment performs smart sampling at block 302 of Fig.3, then Fig.4 is an example process flow chart for that smart sampling processing, in accordance with various embodiments. This sub-process is next described. With reference to Fig.4, process flow begins at block 401, where a set of N reliazations of the reservoir is uploaded. The N relizations may be equiprobable, for example, or, alternatively, they may have different probabilities. Also at block 401, in addition, a number K, which is the desired number of relizations in a subset of N, is uploaded. As noted above, the objective here is to obtain a subset of K realizations that provide identical or similar probability distributions to those of the N realizations, thus simplifying the computations. From block 401 process flow continues to block 402, where, for each of the N realizations, a DNN pretraining for image classification or other deep learning algorithm that allow feature extraction for image or 3D volumes, is perfomed. The DNN may be, for example, C3D. From block 402, process flow moves to block 403, where all features of the intermediate layers are extracted, and a Gram Matrix is computed. The input data may be, for example, sedimentological or petrophysical properties, or, for example, a combination of them, or, even any other property tha a user may decide to use. From block 403, process flow moves to block 404, where a distance matrix is computed. The distance matrix represents a comparison between all of the N realizations. Thus, the principal diagonal of the distance matrix is equal to zero. From block 404, process flow moves to block 405, where a multi-dimensional scaling process is performed, using a multi- dimensional scaling algorithm to represent all of the N realizations in a 2D map. Alternatively, in other embodiments, another projection algorithm may be used here, such as, For example, T-sne. From block 405, process flow moves to block 406, where a clustering algorithm is used, such as, for example, K-means, to determine the subset of K realizations. K is a subset of the original N realizations. Finally, from block 406, process flow moves to block 407, where probabilities are computed, summing the probabilities of all realizations of the cluster. In embodiments, the probabilities are calculated for the number of realizations in the cluster. For example, assuming N is equal to 10 and it is desired to ultimately have 3 clusters (K = 3). The first cluster could have, for example, two realizations, the second cluster seven realizations and the last cluster one realization. Accordingly, in this example, the probability will be 20%, 70% and 10%, respectively for the first, second and third clusters. This probability may then be used later to calulate the expected value, using: Expected value= prob*NPV_realization Process flow then ends at block 407. As noted above with reference to Fig.3, blocks 305 and 306, in embodiments, a surrogate model may be trained to perform the optimization at block 307 of Fig.3, to determine the best FDP. Details of surrogate modelling, in accordance with various embodiments, are next described. It is noted that Fig.5 is an example architecture or system diagram for surrogate modelling, in accordance with such embodiments. Thus, in what follows, repeated reference will be made to the example architecture of Fig.5 when describing the surrogate modelling methodology. Surrogate Modeling for Well-Based Reservoir Simulation 1 Introduction In embodiments, the following methods may be used for well placement in reservoir simulations. The methods include using a surrogate model that incorporates graph convolutional networks (GCNs), a Gated Recurrent Unit (“GRU”), and an attention mechanism. All three of these processing units are shown in Fig.5, which is an example having four layers of GCN, a GRU layer, and an Attention Layer, as shown. GCNs are neural networks designed to capture the spatial relationships in data structured as a graph. In the context of well placement, each well may be viewed as a node in the graph, and the connections between wells can represent the relationships or dependencies between them, such as flow paths. GCNs enable a model to learn the influence of neighboring wells and their properties, which is crucial for predicting well performance, as wells often affect each other in reservoir simulations. GRU is a type of recurrent neural network (RNN) used to model sequential data. In well production prediction, where time-series data is involved, the GRU helps capture the temporal dependencies between past production values and future behavior. GRUs are particularly useful because they can handle long sequences of data efficiently without the risk of vanishing gradients, making them well-suited for predicting well production over time. The attention mechanism (“Attention Layer” in Fig.5) allows the model to focus on important parts of the input data when making predictions. In the context of this task, the attention mechanism can weigh the importance of certain wells or time steps when predicting future production. This is particularly useful when some wells have a more significant impact on the overall production than others, or when certain time intervals in the production history are more informative. Thus, in embodiments, an example architecture may include multiple layers of graph-based processing (e.g., GCNs) configured to learn the influence of neighboring wells and their completions. The number of these GCN layers is not limited to any specific value, and can be adjusted based on the complexity of the reservoir model or other simulation requirements. Additionally, as regards the GRU layer, in alternate embodiments, any layer that can capture sequential and temporal dependencies can be used in place of the GRU layer. The design is not restricted to GRUs and may include other types of recurrent layers or neural units, such as Long Short-Term Memory (LSTM) units or temporal convolution layers, depending on the specific needs of the system. In embodiments, the attention mechanism may be used to highlight the most relevant and important information at each stage of processing. However, this is not limited to traditional attention layers; any component capable of filtering out key features or important time steps may be applied. In general, the architecture is flexible and not limited to the aforementioned components. Alternative layers or mechanisms that perform equivalent functions—such as capturing spatial dependencies, processing sequential data, or emphasizing relevant information—can be used as substitutes, ensuring that the system can adapt to a wide range of scenarios and configurations. In embodiments, the input may include well locations, petrophysical properties, well completions (in the form of ijk coordinates) and well controls. The model is trained to predict time-series well production by incorporating spatial dependencies between wells and temporal dynamics in production data. For a simplistic example, the inputs may be well locations: `[(x1, y1), (x2, y2), (x3, y3)]`, where (x, y) are the spatial coordinates of the wells, petrophysical properties: [(0.25, 100), (0.20, 80), (0.30, 110)], where each pair represents the porosity and permeability at each well location and control at each time: [(1, 100), (2, 213), (3, 173)] where each pair represent bhp control for instance at time steps 1, 2 and 3. In embodiments, the output may be time-series field production: [(1, 180,), (2, 190,), (3, 211)], where each pair represents the field production at time steps 1, 2 and 3, respectively. This example shows how the model takes spatial and temporal inputs (well locations, completions, and controls) and predicts time-varying well production, integrating both spatial dependencies between wells and the dynamics of production over time. Thus, in embodiments, the architecture of the surrogate model is designed to enhance prediction accuracy by efficiently processing the connectivity of wells and integrating key features. 2 Use of Latin Hypercube Sampling (LHS): In embodiments, LHS may be used to ensure uniform sampling of each dimension bypartitioning the design space into ^^^^ equally probable intervals. For a set of parameters ^^^^ ={ ^^^^1, ^^^^2, … , ^^^^ ^^^^} (such as distances ^^^^ ^^^^ ^^^^, ^^^^ ^^^^ ^^^^, ^^^^ ^^^^ ^^^^ ), LHS selects one value from each intervalwithout repetition. This may be mathematically expressed as: ^^^^( ^^^^)^^^^^^^^ − 1 + ^^^^(0,1)^^^^ − ^^^^ where ^^^^^^^^and ^^^^^^^^are the ^^^^^^^^, and ^^^^(0,1) is a uniform random variable. 3. Exemplary methodology to account for distance constraints In embodiments, methods are aimed at generating uniform sampling based on (subject to) predefined real world constraints for the distances between wells. In embodiments, these constraints may, for example, be set for three main variables: • Distance between producers� ^^^^^^^^ ^^^^� • Distance between injectors ^^^^^^^^ ^^^^� •Distance between injectors(^^^^ ^^^^ ^^^^)In embodiments, the well locations may be represented by the set ^^^^ ={^^^^1, ^^^^2, … , ^^^^^^^^}, where each well ^^^^^^^^has a spatial position ^^^^^^^^= ( ^^^^^^^^, ^^^^^^^^, ^^^^^^^^). Moreover, in embodiments, ^^^^ may be defined the set of producer wells, and ^^^^ as the set of injector wells, where ^^^^ ⊆ ^^^^ and ^^^^ ⊆ ^^^^. 3.1 Satisfaction of Distance Constraints: For any two producers ^^^^^^^^, ^^^^^^^^∈ ^^^^ : ^^^^ � ^^^^ , ^^^^� =� ^^^^^^^^ ^^^^^^^^ ^^^^ ^^^^ ^^^^ ^^^^− ^^^^^^^^� ≥ ^^^^min where ^^^^^^^^ ^^^^^^^^ ^^^^is the Euclidean ^^^^minis the minimum allowable distance between two producers. For any injector ^^^^^^^^∈ ^^^^ and producer ^^^^^^^^∈ ^^^^ : ^^^^� ^^^^^, ^^^^^^^^� =� ^^^^^^ ^^^^^^^^ ^^^^ ^^^^^^^^^− ^^^^^^^^� ≥ ^^^^ where ^^^^ is the distance an a^^^^ ^^^^^^^^ ^^^^and ^^^^minis the minimum allowable distance between the injector and producer. For any two injectors ^^^^^^^^, ^^^^^^^^∈ ^^^^ : ^^^^^^^^ ^^^^( ^^^^^^^^, ^^^^^^^^) = ‖ ^^^^^^^^− ^^^^^^^^‖ ≥ ^^^^m^^^^ ^^^^in where ^^^^^^^^ ^^^^is the distanceminthe minimum allowable distance between injectors. For any injector ^^^^^^^^∈ ^^^^ and producer ^^^^^^^^∈ ^^^^, the distance between injectors and producers must also respect an upper limit: ^^^^^^^^ ^^^^≤ ^^^^� ^^^^ , ^^^^� =� ^^^^ − ^^^^^^ ^^^^^^^^ ^^^^ ^^^^ ^^^^ ^^^^^^^^^^� ≤ ^^^^ where: • ^^^^^^^^ ^^^^is the distance between an injector and a producer, • ^^^^^^^^ ^^^^minis the minimum allowable distance between injectors and producers, • ^^^^^^^^ ^^^^maxis the maximum allowable distance between injectors and producers. In embodiments, if the distance between an injector and a producer falls outside this range (i.e., the distances are smaller than the minimum allowable or, for injectors and producers, greater than a maximum allowable), then in such embodiments, the well configuration is replaced through an iterative replacement process. Regarding the constraints, these are typically set by reservoir engineers or geoscientists based on both physical and operational considerations of the reservoir. The minimum distance between producers, ^^^^^^^^ ^^^^min, is typically set to avoid interference between production wells. If two production wells are placed too close, they may start competing for the same resources, reducing their efficiency and leading to suboptimal production. In a typical oil reservoir, a minimum distance of 300 to 500 meters may be set between producers to ensure each well drains its section of the reservoir efficiently without depleting nearby wells. ^^^^^^^^ ^^^^minconstraint is often set to ensure the injected fluids (water or gas) have sufficient time and distance to sweep the hydrocarbons effectively toward the producers. Placing injectors too close to producers may lead to premature breakthrough, where the injected fluid reaches the producer too quickly without sweeping the reservoir efficiently. For a waterflood operation, a minimum distance of around 500 to 1000 meters between injectors and producers is common to ensure that the injected water has enough room to push oil toward the production wells. 3.2 Iterative Process: In embodiments, for each sampled configuration, the distance constraints may then be checked, as follows: For producers: ^^^^^^^^ ^^^^min≤ ^^^^^^^^ ^^^^� ^^^^^^^^, ^^^^^^^^�. For injectors and producers: ^^^^^^^^ ^^^^≤ ^^^^^^^^ ^^^^^^^^ ^^^^� ^^^^^^^^, ^^^^^^^^� ≤ For injectors: ^^^^^^^^ ^^^^ min≤ ^^^^^^^^ ^^^^( ^^^^^^^^, ^^^^^^^^. If any of the violated for the current sample, the following steps, for example, may be taken: Producers: First, producer wells are identified that are outside of the chosen sample that meet all the distance criteria with the other producers. These wells are selected to replace the ones that violate the constraints. Injectors: Second, after all producers meet the distance criteria, the injectors are checked. If any injector violates the distance constraints with producers or other injectors, new injectors from outside the sample are selected to satisfy the constraints. In embodiments, this process of checking and replacing wells continues iteratively until a valid well configuration is found where all wells meet the distance criteria. The final configuration of compliant wells is denoted as ^^^^∗. 4 Surrogate Model Architecture In embodiments, a surrogate model leverages a Graph Convolutional Network (“GCN”) to handle the spatial dependencies between wells, a GRU for temporal processing, and an attention mechanism (“AM”) for focusing on key time steps. As noted, each of these modules is shown in the architecture of Fig.5. 4.1 Input and Output Details In embodiments, the input to the model is a graph where the nodes represent wells and the edges represent the spatial relationships between the wells. This input is shown on the leftmost side of Fig.5, where a set of four wells is represented. The features of the nodes include petrophysical properties and ijk completions of the wells. The ijk features are first through the GCN and then concatenated with the other features. As noted, the ultimate output of the model is the time-series production prediction for each well over a specified number of days. Although not shown in Fig.5, this ultimate output is what is output by the FFNN module at the far right of the figure. 4.2 Graph Convolutional Network (GCN) As shown in FIG.5, in embodiments, an input to the GCN may consist of, for example: The graph connectivity, which connects all wells and also connects the top and bottom completions of each well. This connectivity is based on the well completions (ijk coordinates). As used herein “completion” means a specific location of reservoir where one drills with an ijk index of cells. Typically a reservoir is discretized into cells. The features of the nodes (wells), including petrophysical properties and well completions (ijk). In embodiments, the ijk well completions 503 are first passed to the GCN layer 505, and then the output from the GCN is concatenated at 510 with the rest of the features (petrophysical properties) 501. This concatenation significantly improves the accuracy of the model. In embodiments, each layer of the GCN may transform the input features as follows: ^^^^( ^^^^+1)= ^^^^� ^ˆ^^^ ^^^^( ^^^^)^^^^( ^^^^)� where: ^^^^( ^^^^)is the feature matrix at layer ^^^^, ^ˆ^^^ is the normalized adjacency matrix representing the graph connectivity of the wells, ^^^^( ^^^^)is the learnable weight matrix of the ^^^^-th layer, the ReLU activation function. In embodiments, the final output from the GCN after four layers (523, 525, 527 and 529, in Fig.5) captures the spatial relationships between the wells. This is the “Global add Pool” 530. 4.3 GRU and Attention Mechanism In embodiments, the GRU 540 handles the time-series data related to well controls. The input to the GRU 540 is the well control data over time 570, while the hidden state ℎ0is initialized at Global add Pool 530 using the output of the GCN layer 529 projected to a higher dimension, as follows: ℎ0= ReLU� ^^^^proj^^^^gcn� where ^^^^projis a learnable projection matrix, and ^^^^gcnis the final GCN output. In embodiments, the GRU 540 may then process the temporal sequence using the following dynamics: ^^^^^^^^= ^^^^( ^^^^^^^^^^^^^^^^+ ^^^^^^^^ℎ^^^^−1+ ^^^^^^^^) ^^^^ ^^^^ = ^^^^(^^^ ^^^^^ ^^^^ ^^^^ + ^^^^ ^^^^ℎ ^^^^−1 + ^^^^ ^^^^)ℎ˜ ^^^^ = tanh(^^^^ℎ ^^^^ ^^^^ + ^^^^ℎ(^^^ ^^^^^ ⊙ ℎ ^^^^−1)+ ^^^^ℎ)ℎ^^^^ = (1 − ^^^^ ^^^^) ⊙ ℎ ^^^^−1 + ^^^^ ^^^^ ⊙ ℎ˜ ^^^^where ^^^^^^^^is the input at time step ^^^^, and ℎ^^^^is the hidden state at time step ^^^^. In embodiments, an attention mechanism may then be applied, shown as “Attention Layer” 550 in Fig.5, to the GRU Layer 540 outputs to focus on key time steps, as follows: ^^^^ = ^^^^^^^^ℎ^^^^, ^^^^ = ^^^^^^^^ℎ^^^^, ^^^^ = ^^^^^^^^ℎ^^^^The attention may be computed, for example, as:tention( ^^^^, ^^^^, ^^^^) = softmax�^^^ ^^^^At^ ^^^^� ^^^^ where ^^^^^^^^is the dimension of the key. In the attention mechanism, the "key" refers to the representation used to determine how relevant a particular time step or feature is in relation to others. The dimension of the key, ^^^^^^^^, is essentially the number of features in the key vectors. It defines the size of each key vector that will be used in the dot product computation to determine the attention scores. 4.4 Fully Connected Layers and Output In embodiments, the output of the attention mechanism 550 is then passed through fully connected layers to produce the final production predictions: ^ˆ^^^ = where ^^^^1, ^^^^2, and (weights) in that are learned during the training process. In the context of neural networks, "learnable" means that these matrices start with random values and are adjusted (or trained) using backpropagation and gradient descent. As shown, the output of the attention layer 550 is fed to the final architectural stage, which is, for example, a feed-forward neural network (FFNN) 560, which consists of multiple fully connected layers. FFNN 560 receives the transformed data from the attention mechanism and applies the weight matrices to map this data to the desired output ( ^ˆ^^^). In embodiments, the number of layers in the FFNN is not limited to any specific value and may be adjusted depending on the complexity of the task or data. Furthermore, the final stage is not restricted to a traditional FFNN; any neural network structure or layer type capable of processing the output of the attention mechanism to generate the final prediction may be employed. This may include, for example, convolutional layers, recurrent layers, or hybrid architectures combining different types of neural layers. As noted, the final stage of the architecture is flexible and can be tailored to the specific requirements of the task, allowing for a wide range of configurations beyond just fully connected layers. 5 Training Process In embodiments, the model may be trained to minimize the following loss function, which is a weighted mean absolute error (MAE):1^^^^ where: ^^^^ is the number of samples, ^ˆ^^^^^^^is the predicted production for sample ^^^^, ^^^^^^^^is the actual production for sample ^^^^, ^^^^ is the number of days (time steps). It is understood, however, that the loss function is not limited to MAE and may alternatively include other types of loss functions commonly used in machine learning and deep learning. These may include, for example, mean squared error (MSE), cross-entropy loss, Huber loss, or any custom loss function that balances prediction accuracy with other considerations, such as, for example, regularization or constraint enforcement. Additionally, the loss function may be adapted to include multiple terms or components, allowing for the incorporation of additional objectives, such as penalties for constraint violations or rewards for specific behaviors (for example, higher accuracy for critical outputs). The weighting of these terms may be adjusted based on the desired behavior of the model. In embodiments, the training may be performed, for example, using an Adam (Adaptive Moment Estimation) optimizer with a learning rate scheduler. The optimizer updates the model parameters using the gradient of the loss function. The Adam optimizer is an algorithm used to update the parameters (weights) of a neural network during training. It combines the advantages of two other popular optimizers: gradient descent with momentum and RMSProp. Adam adapts the learning rate for each parameter individually, allowing faster convergence and better results in many cases. In embodiments, a learning rate scheduler may adjust the learning rate during training, often reducing it as training progresses. This helps the model to converge smoothly and avoid overshooting the optimal point in the loss landscape. At the beginning of training, a larger learning rate helps the model explore the parameter space quickly, but as training progresses, a smaller learning rate allows for fine-tuning around the optimal values, improving accuracy. Fig. 5A illustrates an example process flow 580 for a method of optimizing an FDP, in accordance with various embodiments. With reference to Fig. 5A, process flow begins at block 581, where a Well Pool is generated according to Algorithm 1 and block 304 in Fig.3. Continuing with reference to Fig.5A, from block 581, process flow moves to block 582, where a Sampling Training Dataset sub-process is performed. In this sub-process, a design of experiment (DOE) approach is used to generate a training dataset that incorporates distance constraints between injectors, producers and between themselves, as also shown in block 305 in Fig. 3. This ensures that the generated well configurations are consistent with operational rules and physical limitations. The sampled dataset represents a variety of well placement scenarios and serves as input to the deep neural network (DNN). From block 582, process flow moves to block 583, where the training DNN operation is performed. In this step, also shown as block 306 in Fig. 3, the DNN is trained using the sampled training dataset. The training process involves iteratively updating the model’s weights to minimize a loss function, allowing the DNN to learn the relationship between input features, such as, for example, well locations and petrophysical properties, and the predicted well production outcomes. Algorithm 3, described below, demonstrates this process in more detail. After the training process is complete, process flow continues to block 584, where the DNN is considered Ready for Prediction. At this stage, the trained DNN may be used to make predictions for well configurations that were not included in the training dataset. In embodiments, the DNN provides a computational advantage by generating production predictions much faster than traditional reservoir simulations, allowing for rapid evaluations of potential well placements. From block 584, process flow moves to block 585, where the Optimization of FDP is carried out. In this process, also shown as block 307 in Fig.3, the trained DNN is used in an inference mode to evaluate different well configurations from the well pool generated in block 581. The DNN computes an objective function, such as maximizing oil production or maximizing net percent value and determines the optimal set of well locations that best satisfy the objectives of the FDP. Algorithm 3 - Algorithm to Train Deep Neural Network (DNN) Given a training dataset Dtrain; Given user-defined parameters: number of epochs E and batch size B; Initialize deep neural network model DNN with random weights; Initialize optimizer (e.g., Adam) and loss function; for epoch e=1 to E do Shuffle the training dataset Dtrain; Divide Dtrain into batches of size B; for each batch b in Dtrain do Reset the gradients of the DNN; Perform forward pass to compute predicted output ^�^^^ for the input batch b; Calculate loss L( ^^�^^, y) where y is the true output for batch b; Perform backpropagation to compute gradients with respect to the DNN weights; Update weights of the DNN using the optimizer; end for end for Exaplanation of Algorithm 3: In embodiments, the DNN may be trained using an iterative process in which the model's weights are updated to minimize the difference between predicted values and actual target values in the training dataset. The number of epochs (i.e., the number of times the model processes the entire training dataset) and the batch size (i.e., the number of samples processed at once) are user-defined. 1. Epochs and Batches: In some embodiments, the training dataset may be divided into batches of a size specified by a user. For each epoch, the model processes the batches sequentially. This process may be repeated for a user-defined number of epochs. - For each epoch: The DNN processes all batches within the training dataset. For each batch within the dataset, the following operations are performed: 2. Forward Pass (DNN Prediction): For each batch: - Gradients are reset before processing, ensuring that no gradient accumulation from previous batches occurs. - The batch is passed through the layers of the DNN, wherein the DNN computes a predicted output. This computation, referred to as a forward pass, involves passing input data through various neural network layers (such as fully connected layers, convolutional layers, or other architectures), with the final layer producing the predicted output. 3. Loss Function Computation: The predicted output of the DNN is compared to the actual target values corresponding to the batch. A loss function (such as mean absolute error (MAE), mean squared error (MSE), or other suitable loss functions) is computed. This loss function quantifies the discrepancy between the predicted values and the target values. The computed loss serves as the objective function to be minimized during training. 4. Backpropagation and Gradient Descent: Following the computation of the loss function, backpropagation is applied. Backpropagation computes the gradients of the loss function with respect to each weight in the DNN. The optimizer, such as the Adam optimizer, then uses the gradients to adjust the model’s parameters (i.e., weights and biases). Gradient descent or another optimization technique is employed to update the parameters in a direction that reduces the loss, thereby improving the model's accuracy over time. 5. Batch Completion: After processing a batch, the DNN updates the weights based on the calculated gradients and moves to the next batch. This procedure is repeated until all batches in the dataset have been processed, thereby completing one epoch. 6. Epoch Completion: At the end of each epoch, the DNN's weights have been updated based on the data processed in that epoch. The training process continues for the number of epochs specified by the user, with the model weights being incrementally improved after each epoch. 7. Training Completion: Once the user-defined number of epochs is completed, the training process is finalized. The DNN is then ready for deployment and use in making predictions on unseen data. The trained model, along with its learned weights, may be saved in a format suitable for production use. Fig.6 is a block diagram of a computer device suitable for practicing the present disclosure, in accordance with various embodiments. As shown, computer device 600 may include one or more processors 602, and system memory 604. Each processor 602 may include one or more processor cores, and may include a hardware accelerator 605. An example of the hardware accelerator 605 may include, but is not limited to, programmed field programmable gate arrays (FPGA). Each processor 602 may include a DNN Training Module 621, and a DNN Processing Module 625. Further, each processor 602 may also include an FDP Optimizer 631, a Smart Sampling module 633, a Well Pool Generator 635, a Reservoir Clustering module 637, and a Surrogate Model Generation module 639. These latter modules and their functions were described above with reference to Figs.3, 4, 5, and 5A. Computer device 600 may also include system memory 604. In embodiments, system memory 604 may include any known volatile or non-volatile memory. Additionally, computer device 600 may include mass storage device(s) 606, input / output device interfaces 608 (to interface with various input / output devices, such as, mouse, cursor control, display device (including touch sensitive screen, and so forth) and communication interfaces 610 (such as network interface cards, modems and so forth). In embodiments, communication interfaces 610 may support wired or wireless communication, including near field communication. The elements may be coupled to each other via system bus 612, which may represent one or more buses. In the case of multiple buses, they may be bridged by one or more bus bridges (not shown). In embodiments, system memory 604, including mass storage device(s) 606 may be employed to store a working copy and a permanent copy of the executable code of the programming instructions of an operating system, one or more applications, and / or various software implemented components of the various modules and processing functions of Figs.3, 4, 5, and 5A collectively referred to as computational logic 622. The programming instructions implementing computational logic 622 may comprise assembler instructions supported by processor(s) 602 or high-level languages, such as, for example, C, that can be compiled into such instructions. In embodiments, some of computing logic may be implemented in hardware accelerator 605. In embodiments, part of computational logic 622, e.g., a portion of the computational logic 622 associated with the runtime environment of the compiler may be implemented in hardware accelerator 605. The permanent copy of the executable code of the programming instructions or the bit streams for configuring hardware accelerator 605 may be placed into permanent mass storage device(s) 606 and / or hardware accelerator 605 in the factory, or in the field, through, for example, a distribution medium (not shown), such as a compact disc (CD), or through communication interfaces 610 (from a distribution server (not shown)). The number, capability and / or capacity of these elements 602-639 may vary, depending on the intended use of example computer device 600, e.g., whether example computer device 600 is a smartphone, tablet, ultrabook, a laptop, a server, a set-top box, a game console, a camera, and so forth. The constitutions of these elements 602-639 are otherwise known, and accordingly will not be further described. Figure 7 illustrates an example computer-readable storage medium having instructions configured to implement all (or portion of) software implementations of concatenator 510, one or more GCN layers 505 and 523-529, Global add Pool 530, GRU Layer 540, Attention Layer 550, and FFNN 560 of Fig.5, as well as Well Pool 581, Sampling Training Dataset 582, Training DNN 583, and Optimization of FDP 585, all of Fig.5A, and / or practice (aspects of) processes shown in Figs.3, Fig.4 and Fig.5A, earlier described, in accordance with various embodiments. As illustrated, computer-readable storage medium 702 may include the executable code of a number of programming instructions or bit streams 704. Executable code of programming instructions (or bit streams) 704 may be configured to enable a device, e.g., computer device 600, in response to execution of the executable code / programming instructions (or operation of an encoded hardware accelerator 605), to perform (aspects of) processes performed by concatenator 510, one or more GCN layers 505 and 523-529, Global add Pool 530, GRU Layer 540, Attention Layer 550, and FFNN 560 of Fig.5, and / or practice (aspects of) processes illustrated in Figs.3, 4 and 5A. In alternate embodiments, executable code / programming instructions / bit streams 704 may be disposed on multiple non-transitory computer-readable storage medium 702 instead. In embodiments, computer-readable storage medium 702 may be non-transitory. In still other embodiments, executable code / programming instructions 704 may be encoded in transitory computer readable medium, such as signals. Referring back to Fig.6, for one embodiment, at least one of processors 602 may be packaged together with a computer-readable storage medium having some or all of computing logic 622 (in lieu of storing in system memory 604 and / or mass storage device 606) configured to practice all or selected ones of the operations earlier described with reference to Figures 3, 4 and 5A. For one embodiment, at least one of processors 602 may be packaged together with a computer-readable storage medium having some or all of computing logic 622 to form a System in Package (SiP). For one embodiment, at least one of processors 602 may be integrated on the same die with a computer-readable storage medium having some or all of computing logic 622. For one embodiment, at least one of processors 602 may be packaged together with a computer-readable storage medium having some or all of computing logic 622 to form a System on Chip (SoC). For at least one embodiment, the SoC may be utilized in, e.g., but not limited to, a hybrid computing tablet / laptop. Illustrative examples of the technologies disclosed herein are provided below. An embodiment of the technologies may include any one or more, and any combination of, the examples described below. REFERENCES CITED IN THE DESCRIPTION: [1] R. H. Chan, C.-W. Ho, and M. Nikolova. Salt-and-pepper noise removal by median-type noise detectors and detail-preserving regularization. IEEE Transactions on image processing, 14(10):1479–1485, 2005. [2] J. Christensen, E. Stenby, and A. Skauge. Review of wag field experience. society of petroleum engineers, 2001. [3] P. Dufoss´e and N. Hansen. Augmented lagrangian, penalty techniques and surrogate modeling for constrained optimization with cma-es. In Proceedings of the Genetic and Evolutionary Computation Conference, pages 519–527, 2021. [4] D. Erbas, M. Dunning, T. M. Nash, D. Cox, J. A. Stripe, and E. Duncan. Magnus wag pattern optimization through data integration. In SPE Improved Oil Recovery Conference?, pages SPE–169167. SPE, 2014. [5] F. Forouzanfar, G. Li, and A. C. Reynolds. A two-stage well placement optimization method based on adjoint gradient. In SPE Annual Technical Conference and Exhibition?, pages SPE–135304. SPE, 2010. [6] N. Hansen, S. D. M¨uller, and P. Koumoutsakos. Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (cma-es). Evolutionary computation, 11(1):1–18, 2003. [7] J. D. Head and M. C. Zerner. A broyden—fletcher—goldfarb—shanno optimization procedure for molecular geometries. Chemical physics letters, 122(3):264–270, 1985. [8] O. J. Isebor, D. E. Ciaurri, and L. J. Durlofsky. Generalized field-development optimization with derivative-free procedures. Spe Journal, 19(05):891–908, 2014. [9] Y. D. Kim and L. J. Durlofsky. Convolutional–recurrent neural network proxy for robust optimization and closed-loop reservoir management. Computational Geosciences, 27(2):179–202, 2023.

[0010] K. Krishna and M. N. Murty. Genetic k-means algorithm. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), 29(3):433–439, 1999.

[0011] R. Lu, F. Forouzanfar, and A. C. Reynolds. Bi-objective optimization of well placement and controls using stosag. In SPE Reservoir Simulation Conference?, page D012S014R006. SPE, 2017.

[0012] A. Miyagi, Y. Akimoto, and H. Yamamoto. Well placement optimization for carbon dioxide capture and storage via cma-es with mixed integer support. In Proceedings of the Genetic and Evolutionary Computation Conference Companion, pages 1696–1703, 2018.

[0013] J. Navr´atil, A. King, J. Rios, G. Kollias, R. Torrado, and A. Codas. Accelerating physics-based simulations using end-to-end neural network proxies: An application in oil reservoir modeling. Frontiers in big Data, 2:33, 2019.

[0014] M. Panda, J. G. Ambrose, G. Beuhler, and P. McGuire. Optimized eor design for the eileen west end area, greater prudhoe bay. SPE Reservoir Evaluation & Engineering, 12(01):25–32, 2009.

[0015] S. D. Rahmawati, C. H. Whitson, and B. Foss. A mixed-integer non-linear problem formulation for miscible wag injection. Journal of Petroleum Science and Engineering, 109:164–176, 2013.

[0016] M. Sarabian, P. Ruiz Mataran, and R. Rodriguez Torrado. Finite pinn net: Physics- informed deep convolutional neural networks for learning 3d transient darcy flows in heterogeneous porous media. Bulletin of the American Physical Society, 2022.

[0017] R. R. Torrado, D. E. Ciaurri, U. T. Mello, and S. E. Droz. Fast reservoir performance evaluation under uncertainty: opening new opportunities. In SPE Annual Technical Conference and Exhibition, page D011S007R005. SPE, 2013.

[0018] R. R. Torrado, G. D. Paola, A. F. Perez, A. Fuenmayor, M. S. De Azevedo, and S. Embid. Optimize a wag field development plan, use case of carbonate ultra-deep water reservoir. In SPE Europec featured at EAGE Conference and Exhibition, pages SPE– 174344. SPE, 2015.

[0019] D. Yang, Q. Zhang, H. Cui, H. Feng, and L. Li. Optimization of multivariate production- injection system for water-alternating-gas miscible flooding in pubei oil field. In SPE Western Regional Meeting, pages SPE–62856. SPE, 2000.

[0020] S. Yu, J. Ma, and W. Wang. Deep learning for denoising. Geophysics, 84(6): V333– V350, 2019.

Claims

WHAT IS CLAIMED:

1. An apparatus for computing, comprising: an input interface, configured to receive: a base reservoir simulation representative of a geological oil reservoir; a processor, the processor configured to: cluster the reservoir; generate, in a first process, a possible pool of well positions in the reservoir, where each well position satisfies plausibility and respects operational constraints; and optimize, in a second process, the pool of well positions, the second process maximizing a dynamic variable of the reservoir; and an output interface, configured to output the optimized pool of wells to a user. 2 The apparatus for computing of claim 1, wherein the output optimized pool of well positions specifies a set of producer wells and injector wells that is a subset of the possible pool of well positions generated in the first process. 3 The apparatus for computing of claim 1, wherein the base reservoir simulation is a set of cases with different maps of porosity and permeability that characterize geological uncertainty of the reservoir. 4 The apparatus for computing of claim 1, wherein the processor is further configured, prior to clustering the reservoir, to perform a smart sampling sub-process. 5 The apparatus for computing of claim 5, wherein the smart sampling sub-process outputs a subset of K realizations that provide similar probability distributions to those of an original set of N realizations.

6. The apparatus for computing of claim 5, wherein the smart sampling sub-process determines the smallest set of realizations that allows for a similar expected value of the reservoir objective function to be obtained compared to the full set of realizations.

7. The apparatus for computing of claim 1, wherein the clustering includes computing different clusters of the reservoir, each cluster representing a sub-region of the reservoir, using a clustering algorithm. 8 The apparatus for computing of claim 7, wherein the clustering algorithm is K-means. 9 The apparatus for computing of claim 8, wherein the processor is further configured to post-process the output of K-means clustering to remove salt-and-pepper noise. 10 The apparatus for computing of claim 9, wherein the post processing is performed using at least one of: median filters, or a pre-trained DNN to remove the noise. 11 The apparatus for computing of claim 1, wherein the dynamic variable of the reservoir to be optimized in the second process is net present value (“NPV”). 12 The apparatus for computing of claim 11, wherein NPV is calculated as: SUM from t=1 to Nt[[Q0 to tn(p,q)∗ Po− Nw∗ Pw] / [(1+λ)tn / 365]*[Δtn]], where: Qo indicates the field-wide oil production rate in stock tank barrels per day (STB / d); Nwindicates the total number of wells drilled; PoPw, and λ respectively indicate: oil price, the cost to drill a well, and the annual discount rate and, tnand Δtnindicate time and the size of the time step in days. 13 The apparatus of computing of claim 1, wherein the dynamic variable is recovery factor or oil production.

14. The apparatus for computing of claim 1, wherein the operational constraints include at least one of: distance between two wells, maximum number of injectors to drill, maximum number of producers to drill, distance between two producer wells, distance between two injector wells and distance between a producer well and an injector well.

15. The apparatus for computing of claim 1, wherein, in the second process, the processor is further configured to test a potential solution well x for feasibility, and if it is not feasible, to replace x with the closest well to x that is feasible, given pre-determined feasibility criteria for the reservoir.

16. The apparatus for computing of claim 1, wherein the processor is further configured, as part of the second process, to: train a surrogate model DNN to forecast the production of the reservoir; and use the trained surrogate model to optimize the pool of wells.

17. A method of optimizing a set of oil producing wells to be drilled in a reservoir, comprising: receiving a base reservoir simulation representative of a geological oil reservoir; clustering the reservoir; generating, in a first process, a possible pool of well positions in the reservoir, where each well position satisfies plausibility and respects operational constraints; and optimizing in a second process, the pool of well positions, the second process maximizing a dynamic variable of the reservoir; and outputting the optimized pool of wells to a user.

18. The method of claim 17, wherein the output optimized pool of well positions specifies a set of producer wells and injector wells that is a subset of the possible pool of well positions generated in the first process.

19. The method of claim 17, wherein the base reservoir simulation is a set of cases with different maps of porosity and permeability that characterize geological uncertainty of the reservoir.

20. The method of claim 17, further comprising, prior to clustering the reservoir, performing a smart sampling sub-process.

21. The method of claim 20, wherein the smart sampling sub-process outputs a subset of K realizations that provide similar probability distributions to those of an original set of N realizations.

22. The method of claim 21, wherein the smart sampling sub-process determines the smallest set of realizations that allows for a similar expected value of the reservoir objective function to be obtained compared to the full set of realizations.

23. The method of claim 17, wherein the clustering includes computing different clusters of the reservoir, each cluster representing a sub-region of the reservoir, using a clustering algorithm.

24. The method of claim 23, wherein the clustering algorithm is K-means.

25. The method of claim 24, wherein further comprising post-processing the output of the K- means clustering to remove salt-and-pepper noise.

26. The method of claim 25, wherein the post processing includes using at least one of: median filters, or a pre-trained DNN to remove the noise.

27. The method of claim 17, wherein the dynamic variable of the reservoir to be optimized in the second process is net present value (“NPV”).

28. The method of claim 27, wherein NPV is calculated as: SUM from t=1 to Nt [[Q0 to tn(p,q)∗ Po− Nw∗ Pw] / [(1+λ)tn / 365]*[Δtn]], where: Qoindicates the field-wide oil production rate in stock tank barrels per day (STB / d); Nwindicates the total number of wells drilled; Po, Pw, and λ respectively indicate: oil price, the cost to drill a well, and the annual discount rate; and, tnand Δtnindicate time and the size of the time step in days.

13. The apparatus of computing of claim 1, wherein the dynamic variable is recovery factor or oil production.

29. The method of claim 17, wherein the operational constraints include at least one of: distance between two wells, maximum number of injectors to drill, maximum number of producers to drill, distance between two producer wells, distance between two injector wells and distance between a producer well and an injector well.

30. The method of claim 17, further comprising, in the second process, testing a potential solution well x for feasibility, and if it is not feasible, replacing x with the closest well to x that is feasible, given pre-determined feasibility criteria for the reservoir.

31. The method of claim 17, further comprising, as part of the second process: training a DNN surrogate model to forecast the production of the reservoir; and using the trained DNN surrogate model, optimizing the pool of wells.

Images