A Novel Surrogate-Assisted Multi-Objective Well Control Parameter Optimization Method Based on Selective Ensembles

Abstract

Multi-objective optimization algorithms are crucial for addressing real-world problems, particularly with regard to optimizing well control parameters, which are often computationally expensive due to their reliance on numerical simulations. Surrogate-assisted models help to reduce this computational burden, but their effectiveness depends on the quality of the surrogates, which can be affected by candidate dimension and noise. This study proposes a novel surrogate-assisted multi-objective optimization framework (MOO-SESA) that combines selective ensemble support-vector regression with NSGA-II. The framework’s uniqueness lies in its adaptive selection of a diverse subset of surrogates, established prior to iteration, to enhance accuracy, robustness, and computational efficiency. To our knowledge, this is the first instance in which selective ensemble techniques with multi-objective optimization have been applied to reservoir well control problems. Through employing an ensemble strategy for improving the quality of the surrogate model, MOO-SESA demonstrated superior well control scenarios and faster convergence compared to traditional surrogate-assisted models when applied to the SPE10 and Egg reservoir models.

1. Introduction

One of the most essential steps prior to production and after drilling and completion is numerical simulation-based well control optimization. This refers to finding an optimum scheme that could achieve higher productivity and profitability [1,2]. The net present value (NPV), which displays the comprehensive performance of developing the reservoir, is usually considered as the single objective for optimizing well control schemes [3]. However, oilfield development is complicated, as well control optimization problems are generally related to multiple objectives, including optimizing profitability and oil recovery while reducing operation fees and investments [4]. Moreover, these well control-related objectives may conflict with each other [5]. Two main methods, namely decomposition and the direct method, were utilized to search for the optimal solution of multi-objective problems. The drawback of the decomposition method, such as weighted aggregation, is the selection of the correct weight setting. In comparison, when the direct method is used, the optimum is reached with no human interference. By considering all objectives simultaneously, multi-objective optimization provides a compromise solution set for decision-makers. Individuals in the compromise solution set are non-dominated to each other, and together, they comprise the Pareto front (PF) [6]. To that end, multi-objective optimization, which simultaneously takes into account more than one objective, is more suitable for well control optimization [7].

In previous studies on addressing well control multi-objective optimization problems, many researchers have paid attention to two main techniques [8,9], namely gradient-based and gradient-free techniques, to obtain high-quality solutions that are represented by diversity and convergence. The gradient-based strategy is one of the traditional techniques, and it is also perfect in the field of computational efficiency during the process of numerical simulation [10]. Within the realm of gradient-based methods, the adjoint strategy, which relies on the simulation model and the reverse adjoint equation, results in neglecting the dimension of optimal variables. In this way, the adjoint technique is an efficient way to address well control optimization problems with practical reservoirs. Nevertheless, obtaining the adjoint gradient requires revising the code of simulation software to acquire the Jacobian matrix, which limits the application of the adjoint method [11,12]. Furthermore, the adjoint method has the limitation of falling into the local optimal while missing the global optimal [13]. Meanwhile, the gradient-free algorithm has also been extensively utilized in multi-objective reservoir optimization problems in recent years because of its broad suitability [14,15]. Specifically, Yasari et al. (2013) presented a robust multi-objective well control parameters optimization framework that considered three specified NPVs as objectives [16]. Luo et al. (2015) established a multi-objective well control parameters optimization framework with the particle swarm optimization (PSO) algorithm and estimation of distribution algorithm [17]. Fu and Wen (2017) compared the multi-objective particle swarm optimization and PSO in production optimization problems [18]. However, evaluating a lot of objectives to search for the optimum solution increased the computational burden and reduced the calculation efficiency. These studies have achieved impressive performance in the optimization of well control parameters, but the simulations that run during the optimization process are time-consuming. According to the grid size of the numerical model, it takes hours to run the simulation process. A natural concern that emerges is how to simultaneously reduce the computational load and maintain accuracy.

To lighten the computational cost and increase the optimization efficiencies during the well control optimization process, Proper Orthogonal Decomposition (POD), data-driven model and machine learning-based surrogate methods are the most employed methods. POD is brief in its computation, simple to execute, and offers a perfect approximation for complex numerical simulation systems [19]. Using the POD model, the optimization process could be more efficient. Similar to the traditional gradient technique, POD also requires the simulator code to be modified, which is difficult for a commercial simulator. Without any prior knowledge of geological and numerical models [20], a data-driven model could serve as a surrogate for well control optimization based on history matching of production data [21]. There are three main types of data-driven models, and these include the capacitance-resistance model, interwell numerical simulation model, and its front tracking [22,23,24]. Meanwhile, the main challenge associated with the data-driven model is that a lot of production data are required. As one of the most widely used approaches, machine learning theory-based surrogate methods can flexibly model any inputs and outputs relationship [2] with enough training samples, such that the expensive computation of the simulator could be lessened during the optimization process. Many machine learning approximation methods, such as the Gaussian process [25], artificial neural network [26], convolutional neural networks [27], support vector regression [23,28], and radial basis function [29], have been successfully employed in establishing surrogate models for computationally expensive problems. However, these are single or double models applied as surrogates, such that the quality of the model is difficult to guarantee. Meanwhile, the accuracy of machine learning surrogates is highly affected by the number of optimal variables and the scale of the training samples. With that in mind, the first task is to enhance the quality and robustness of the models with limited training samples. It is worth mentioning that the surrogate ensemble strategy, which constructs several single surrogate models, is suitable for improving the quality of the optimal solutions and has not been utilized in the field of solving well control optimization problems.

The goal of this study is to propose a novel surrogate model-based well control optimization framework with machine learning and intelligent algorithm, which could make full use of training samples to guide the multi-objective optimization search. To that end, several support vector regression (SVR) models are established before iteration and a simple, yet diverse, subset is adaptively selected, based on ensemble learning during the optimization process, to develop the surrogate ensemble while reducing the computational cost. More importantly, this is the first application of selective ensemble techniques with multi-objective optimization to reservoir well control problems.

2. The Multi-Objective Well Control Optimization Model

In well control optimization, the model typically involves identifying optimal variables and defining objective functions to improve production efficiency and profitability. In this study, the optimal variables include key well control parameters such as injection rates and bottom hole pressures, which directly influence reservoir performance. To assess the effectiveness of these strategies, two main objective functions are established. The first is NPV, a classic economic indicator used to evaluate the profitability of production by accounting for revenues, operating costs, and discounting future income. NPV serves as the primary objective in this study, aiming to maximize the financial returns from the reservoir. The second objective is cumulative oil production (COP), which measures the total volume of oil extracted and focuses on maximizing recovery efficiency. Because NPV and COP can sometimes conflict—where maximizing one may impact the other—this is a multi-objective optimization problem. The goal is to explore trade-offs between profitability and recovery, ultimately finding a balance that ensures both economic success and efficient reservoir depletion. The relationships between decision variables and objectives are formed as follows [2]:

$$J_1\left(u\right)={\displaystyle \sum _{n=1}^{N_t}\frac{\mathsf{\Delta }{t}_n}{{\left(1+d\right)}^{\frac{t_n}{365}}}\left[{\displaystyle \sum _{j=1}^P\left(r_o\cdot q_{o,j}^n-r_w\cdot q_{w,j}^n\right)-{\displaystyle \sum _{j=1}^I\left(c_{wi}\cdot q_{wi,j}^n\right)}}\right]}$$

(1)

$$J_2\left(u\right)={\displaystyle \sum _{n=1}^{N_t}{\displaystyle \sum _{j=1}^P{q}_{o,j}^n}}$$

(2)

The specific meanings of each symbol of Equations (1) and (2) are listed in the Nomenclature Section.

3. Related Method and Algorithm

3.1. Selective Surrogate Ensembles

Ensemble learning, as a surrogate management technique for machine learning, constructs a lot of base models and combines them as a strong model with more accuracy and robustness than a single model [30]. The two mainstream ensemble methods, namely boosting and bagging, have different strategies. Bagging involves minimizing variance using a parallel ensemble while boosting is minimizing bias by a sequential ensemble [31]. Moreover, during the iteration process, the bias of surrogates is less important. For this reason, bagging was adopted in this study to establish the ensemble. Meanwhile, subset learners were selected to improve the quality of the surrogate ensemble. Figure 1 consists of four steps that show the workflow of selective bagging.

As can be seen in Figure 1, each datapoint in N subsets (S₁, S₂, …, S_N) is randomly selected by a bootstrap that resamples the initial dataset. The traditional way of deciding the size of each subset is half-sampling. To make subsets more representative, a probability-based sampling strategy was employed in this study. Each point in the initial dataset was set to have a probability of 50% to be selected for each subset. This means the size of subsets is not fixed, which promotes the diversity of the ensembles. Then, the N SVR models were trained based on subsets S₁, S₂, …, and S_N separately. The specific process of establishing the surrogate ensemble is elaborated in Algorithm 1.

Algorithm 1 Pseudocode of establishing the surrogate ensemble

Input: D-the dataset, N-the number of models

1: for i from 1 to N

2:  Empty Si

3:  for individual in D

4:    if U(0,1) < 50%

5:      Add the individual to Si

6:    end if

7:  end for

8: end for

9: for i from 1 to N

10:  Establish an SVR model Mi using Si

11: end for

Output: subset dataset (S1, S2, …, SN) and model set (M1, M2,…, MN)

Another challenge related to selective ensembles relates to the type and method employed to select the SVR models for combinations. Within limits, the accuracy of surrogate ensemble improves with the increase in the selective ensemble size [32]. This indicates that the model selection strategy is essential for reducing the size of the ensemble without decreasing its quality. To enhance the accuracy and robustness of the surrogate ensemble, a random selection strategy was proposed for selecting SVR models in each generation. More specifically, fixed Q models were selected from the N models in each generation. Choosing the best individual (g_b) in the current generation, fitness values (F₁, F₂, …, F_N) of g_b with N models (M₁, M₂, …, M_N) were evaluated. Subsequently, the N SVR models were sorted based on the fitness values and separated into equal Q groups. To this end, one SVR for each data group was randomly selected and the surrogate ensemble for the next iteration was established. Algorithm 2 details the model selection strategy of the selective ensemble.

Algorithm 2 Pseudocode of model selection strategy

Input: Q-the scale of ensemble after selection, gb-the best individual of current generation, (M1, M2,…, MN)-the model set

1: if initial generation

2:  Select Q models randomly for model set

3: else

4:   Calculate fitness values (F1, F2, …, FN) of gb with (M1, M2, …, MN)

5:   Sort N SVR models by (F1, F2, …, FN)

6:   Divide N sorted SVR into equal Q groups

7:   for i = 1: Q

8:     Randomly select a model from group i and establish surrogate model

9:     end for

10: end if

Output: SVR models

A model set with eight models (M₁, M₂, …, M₈) is taken as an example to display the principle of selection strategy. Let the fitness value of g_b be F₁ = 1.1, F₂ = 1.4, F₃ = 1.6, F₄ = 1.5, F₅ = 1.7, F₆ = 1.8, F₇ = 1.2, F₈ = 1.3, respectively. To construct a surrogate ensemble with four models for the next iteration, these models are ranked according to fitness values and divided into four groups (F₁, F₇), (F₈, F₂), (F₄, F₃), and (F₅, F₆). Then, one model is chosen at random from every group. This way, F₁, F₈, F_4, and F₆ could be an output case.

3.2. Optimization Algorithm (NSGA-II)

NSGA-II, as the most traditional and perfect multi-objective optimization algorithm, was introduced by Deb [6]. Its workflow includes a specific sequence of cross, mutative, and non-dominated sorting. The key distinction of NSGA-II from other algorithms lies in its sorting strategy, where individuals are ranked based on crowding distance. Additionally, the elitism strategy applied in NSGA-II ensures that the best individuals are preserved for the next iteration, promoting faster convergence. Given its strong performance, NSGA-II was selected as the algorithm for this study.

3.3. Specific Workflow of MOO-SESA

The main framework of MOO-SESA, including two essential strategies, namely model selection and non-dominated sorting, is displayed in Figure 2. Firstly, the decision variables and objectives were determined as the initial conditions. Next, the initial individuals were produced by LHS and related objective values were calculated by the simulator to construct the dataset. Then, a selective ensemble SVR model was constructed to reduce the computational burden based on the dataset and ensemble strategy. This could make full use of limited samples and lead to more efficient and accurate surrogates. During each iteration, the surrogate ensemble was reconstructed. Finally, the NSGA-II was performed as the optimizer based on non-dominated sorting. The selective surrogate ensemble enhances the quality of prediction accuracy and releases the computational burden.

4. Case Study

To evaluate the performance of the MOO-SESA framework, two synthetic reservoir models are employed and the results are compared with those from the classical surrogate-based optimization framework that integrates a single static surrogate model with an optimization algorithm (SVR-NSGA-II). The first reservoir model was a two-dimensional extract from the SPE10 model, while the second was a three-dimensional Egg model based on Jansen et al. (2014) [33]. Simulations were conducted using the MATLAB Reservoir Simulation Toolbox (MRST) as described by Lie (2019) [34].

4.1. Case 1

To assess the performance of the MOO-SESA framework in this scenario, a two-dimensional heterogeneous reservoir model derived from a portion of layer four of the SPE10 model was utilized. This model is configured in a five-spot pattern, comprising one central injector and four producers located at the corners of the grid. The grid consists of 50 × 50 × 1 blocks, with each cell measuring 6 × 6 × 10 m. The arrangement places the injector at the center of the model, optimizing the injection and production process. The logarithmic permeability field for this setup is depicted in Figure 3.

In this case study, the decision variables are the injection rate of the injector and the bottom hole pressures (BHP) of the four producers. The injection rate ranges from 0 to 150 m³/day, and the BHP for the producers varies between 5 MPa and 25 MPa. The total production period is 600 days, divided into six intervals of 100 days each. This setup results in 30 decision variables.

Table 1. The values of parameters for calculating NPV.

Parameters	2	Value
Water disposal	rw	USD 3/STB
Water injection	cwi	USD 2/STB
Oil revenue	ro	USD 60/STB
Annual discount rate	d	0.05

provides the economic parameters used to calculate the net present value (NPV).

The accuracy and computational burden of the surrogate ensemble depends significantly on the size of the selected models [35]. To determine the optimal ensemble size, we examined the RMSE [36] and R² [37] indexes for ensemble sizes up to 100. Both NPV and COP were tested to identify the appropriate ensemble size and validate the performance of the proposed selective surrogate ensemble strategy. For this analysis, 200-well control samples were generated using Latin hypercube sampling, and the objectives were evaluated by a simulator to construct the dataset. Subsequently, 100 subsets of 150 samples each were generated via bootstrap sampling, and 100 SVR models were built using Algorithm 1. Models were randomly selected from the set to form surrogate ensembles, and their accuracy was blind tested using 50 samples.

As shown in Figure 4, RMSE decreased as the ensemble size increased, stabilizing when the size exceeded 40, while R² continued to improve, achieving significant precision around an ensemble size of 4. However, the computational burden increased linearly with ensemble size. Based on these observations, an ensemble size of 40 was chosen, providing a balance between sufficient accuracy and manageable computational load.

Based on the dataset, bootstrap sampling was performed for 50 subsets with 150 samples, and 50 SVR models were built to establish a surrogate ensemble for optimization. The specific parameter settings of the optimization process are shown in

Table 2. Values of parameters used in NSGA-II.

Parameters	Value
Number of individuals	50
Iteration steps	100
Probability of mutative	0.005
Probability of cross	0.65

.

Figure 5 shows the Pareto fronts of the SVR-NSGA-II and the MOO-SESA. The gray hollow circles display the trajectories of intermediate solutions of the MOO-SESA. The red solid circles represent the PF of the MOO-SESA. The blue solid circles are the PF of the SVR-NSGA-II. For the same number of iterations, the MOO-SESA framework produced a PF with superior diversity and convergence compared to the SVR-NSGA-II method. Furthermore, MOO-SESA resulted in an increased Pareto front compared to that achieved using SVR-NSGA-II. To that end, the optimization results of MOO-SESA were better than SVR-NSGA-II. It is worth mentioning that the major calculation burden of MOO-SESA in this example is only 200 times the numerical simulation run plus the very small amount of time required to train 50 SVR models. However, without the surrogate model, the computational burden runs 5000 simulations in this case, which is nearly 25 times slower than the MOO-SESA method. Meanwhile, the cost of the surrogate modeling and selection process is nearly 5% of the runtime of the numerical simulation model. For a decision maker, any solutions of well control schemes in the PF could be a perfect result.

For the purpose of further assessing the effectiveness of the MOO-SESA method, two optimization schemes—maximizing NPV (case 1) and maximizing COP (case 2)—were selected for detailed analysis. The detailed schemes for both cases are illustrated in Figure 6, while their corresponding NPV and COP values are shown in Figure 7. Case 1 achieved an NPV of USD 9.69 × 10⁶ with a COP of 7.36 × 10⁴ m³, whereas case 2 resulted in an NPV of USD 7.52 × 10⁶ and a COP of 9.02 × 10⁴ m³. This inverse relationship between NPV and COP clearly demonstrates the key characteristic of multi-objective optimization, where improving one objective often leads to a reduction in another.

4.2. Case 2

In this case, Egg model, which is a small synthetic reservoir model proposed by Jansen et al. (2014) [33], is selected as the 3D model to validate the quality and efficiency of the proposed MOO-SESA method. The grid block of this model is 60 × 60 × 7, totaling 25,200 grids, of which the active cells are 18,553. Because there are no active grids at the edge of the model, the model appears as an egg. For detailed and more specific parameters of the Egg model, one can refer to the work published by Jansen et al. (2014) [33]. The detailed permeability field is shown in Figure 8.

In this scenario, the production wells were configured to operate at a constant pressure of 38 MPa. Consequently, the decision variables are the injection rates of the eight injectors, each bounded between 0 m³/day and 70 m³/day. The total production period spanned 3600 days, divided into twenty intervals of 180 days each. This setup resulted in 160 decision variables. The economic parameters for calculating the net present value (NPV) were the same as those used in the first case.

To determine the optimal ensemble scale, RMSE and R² metrics were used to evaluate the accuracy of ensembles of up to 100 members. The study focused on NVP and COP, with decision variables totaling 160, significantly more than in the 2D case. Using Latin hypercube sampling (LHS), 1000-well control samples were generated and evaluated by the simulator to create the dataset. Subsequently, 100 SVR models were built from bootstrap samples of 950, and another 50 samples were used for blind testing. The RMSE and R² trends with varying ensemble sizes are depicted in Figure 9. As shown in Figure 9, RMSE decreased as the ensemble size increased, but stabilized within a narrow range once the ensemble size exceeded 48. Similarly, R² improved with larger ensemble sizes, achieving considerable precision around an ensemble size of 35. However, the linear increase in computational complexity with larger ensemble sizes remained a challenge. Based on these findings, an ensemble size of 48 was selected to balance sufficient accuracy with a relatively modest computational burden.

Using bootstrap sampling, 60 subsets with 950 samples were generated to build 60 SVR models, forming the surrogate ensemble for optimization. The optimizer’s parameters were consistent with those in the first case. As illustrated in Figure 10, the Pareto fronts obtained by the MOO-SESA method outperformed those of the SVR-NSGA-II, providing better diversity and convergence within the same number of iterations. The MOO-SESA method also significantly reduced the computation time, with ensemble surrogate models running in under 10 s, compared to over 300 s for the numerical simulation on an Intel Core i7-7800 CPU. This approach makes high-dimensional, computationally expensive optimization feasible, offering decision-makers a broader range of solutions tailored to specific objectives, such as maximizing economic benefits or recovery ratios.

To further assess the performance of the MOO-SESA method, we selected two extreme scenarios—one maximizing NPV (case 1) and the other maximizing COP (case 2)—to illustrate the details of the optimized solutions. The well control schemes for cases 1 and 2 during each production interval are depicted in Figure 11.

The specific NPV and COP curves are displayed in Figure 12. Cases 1 and 2 yielded NPVs of USD 9.14 × 10⁷ and USD 7.68 × 10⁷, respectively. Meanwhile, for COP, cases 1 and 2 produced 3.45 × 10⁵ m³ and 4.71 × 10⁵ m³, respectively. In essence, the maximum NPV corresponded to the minimum COP, and vice versa, illustrating a fundamental characteristic of multi-objective optimization (MOO). This means that improving one objective inevitably leads to a decline in one or more other objectives.

5. Conclusions

This study introduced the MOO-SESA framework, an innovative optimization method designed to address the computational challenges of well control optimization problems. Our comparative analysis, which was applied to both the SPE10 2D and Egg 3D reservoir models, demonstrated that MOO-SESA significantly outperformed SVR-NSGA-II across multiple metrics. Specifically, the MOO-SESA framework delivered more efficient and superior well control scenarios with enhanced diversity and convergence of the Pareto front. This led to faster convergence rates, greater accuracy in capturing the trade-offs between competing objectives, and ultimately better decision-making for well control optimization. The results confirmed that the MOO-SESA framework’s adaptability and superior performance are not limited to specific reservoir types or dimensions; it proved effective in both 2D and 3D models. This makes the MOO-SESA framework highly versatile and applicable across a range of reservoir simulation challenges, from simple, low-dimensional problems to more complex, high-dimensional ones. By significantly reducing computational costs while improving optimization quality, MOO-SESA holds tremendous promise for more efficient and robust decision-making in petroleum reservoir management. Moreover, the MOO-SESA framework holds potential for broader applications within petroleum engineering, including well placement and fracturing optimization. Furthermore, integrating MOO-SESA with real-time optimization frameworks is also worth investigating. Future research will focus on developing more accurate and efficient surrogate models and optimizing algorithms that are better suited to approaching high-dimensional, computationally complex multi-objective reservoir optimization challenges.