Regular Article
Modelbased production optimization under geological and economic uncertainties using multiobjective particle swarm method
Faculty of Petroleum Engineering, Amirkabir University of Technology, 424 Hafez Ave, Tehran 158754413, Iran
^{*} Corresponding author: drbdabir@aut.ac.ir
Received:
11
January
2021
Accepted:
28
June
2021
Optimization of the waterflooding process in the oilfields is inherently subject to several uncertainties arising from the imperfect reservoir subsurface model and inadequate data. On the other hand, the uncertainty of economic conditions due to oil price fluctuations puts the decisionmaking process at risk. It is essential to handle optimization problems under both geological and economic uncertainties. In this study, a Paretobased MultiObjective Particle Swarm Optimization (MOPSO) method has been utilized to maximize the shortterm and longterm production goals, robust to uncertainties. Some modifications, including applying a variable in the procedure of leader determination, namely crowding distance, a corrected archive controller, and a changing boundary exploration, are performed on the MOPSO algorithm. These corrections led to a complete Pareto front with enough diversity on the investigated model, covering the entire solution space. Net Present Value (NPV) is considered the first goal that represents the longterm gains, while a highly discounted NPV (with a discount rate of 25%) has been considered shortterm gains since economic uncertainty risk grows with time. The proposed optimization method has been used to optimize water flooding on the Egg benchmark model. Geological uncertainty is represented with ensembles, including 100 model realizations. The kmeans clustering method is utilized to reduce the realizations to 10 to reduce the computing cost. The Pareto front is obtained from Robust Optimization (RO) by maximizing average NPV over the ensembles, as the conservative production plan. Results show that optimization over the ensemble of a reduced number of realizations by the kmeans technique is consistent with all realizations’ ensembles results, comparing their cumulative density functions. Furthermore, 10 oil price functions have been considered to form the economic uncertainty space. When SNPV and LNPV are optimized, considering uncertainty in oil price scenarios, the Pareto front’s production scenarios are robust to oil price fluctuations. Using the robust Pareto front of LNPV versus SNPV in both cases, one can optimize production strategy conservatively and update it according to the current reservoir and economic conditions. This approach can help a decisionmaker to handle unexpected situations in reservoir management.
© M.M. Moshir Farahi et al., published by IFP Energies nouvelles, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
B : Boundary in MOPSO algorithm
BF: Boundary factor in MOPSO algorithm
C_{r} : Rock’s compressibility
C_{w} : Water’s compressibility
J_{RO} : Robust objective function
k : Number of means in kmeans method
K_{ro0} : Relativepermeability’s endpoint to oil
K_{rw0} : Relativepermeability’s endpoint to water
M : Number of problem constraints
NDSW: Worst NonDominated Solution
NDSE: Efficient NonDominated Solution
N_{e } : Number of realizations
NPV_{RO}: Robust Net Present Value
N_{geo} : Number of geological realizations
N_{eco} : Number of economic (oil price) realizations
P_{res} : Reservoir’s pressure
pr_{o} : Oil price – production
pr_{wp} : Water price – production
pr_{wt} : Water price – injection
Q_{o,k} : Oil production rate at k’s timestep
Q_{wp,k} : Rate of water rate at timestep kproduction
Q_{wi,k}: Rate of water rate at timestep kinjection
SaF: Saturation factor in MOPSO algorithm
SF: Shrink Factor in MOPSO algorithm
s ^{(t)} : Clusters in kmeans method
S_{wir} : Irreducible watersaturation
S_{wi} : Initial watersaturation
S_{or} : Residual oilsaturation
t_{k } : total time in k’s timestep
x_{p } : Data points in kmeans method
α_{i } : Randomly selected coefficients in ARMA model
β_{i } : Randomly selected coefficients in ARMA model
Δt_{k}: Difference of successive timesteps
τ_{ref} : Discounting referenced time
1 Introduction
One of the main concerns of the decisionmaking process in reservoir management is highlevel uncertainties and debatable parameters. This effect is a crucial challenge, especially in modelbased optimization. Using various realizations of the geological model and/or oil price function identifies such uncertainties for the investigated fields. Meanwhile, sometimes, it is needed to optimize two or more possibly conflicting objective functions in reservoir management. Therefore, MultiObjective Optimization (MOO) is a practical tool for decisionmaking that should be considered while identifying uncertainties.
Waterflooding optimization has been investigated in many kinds of research in petroleum literature, which focused mainly on singleobjective problems [1–3]. Accordingly, the field lifecycle’s economic performance has been reflected in the form of Longterm Net Present Value (LNPV) [4–6]. Although, neglecting the shortterm goals in an intended optimal design generally leads to a low profit in the shortterm compared to that of a reactive strategy. A limited constraint, e.g., maximum allowed watercut, is the only controller for designing the water injection scenario in a reactive strategy, selecting which conflicts with LNPV maximization in a reservoir. However, shortterm production gains are included in very recent studies [7–11]. One can generally use an annual discount rate in NPV calculation to consider the shortterm financial targets. Accordingly, the simultaneous optimization of LNPV and Shortterm NPV (SNPV) via an appropriate MOOhandling algorithm to form a Pareto front provides a valuable mechanism that helps the decisionmaker tradeoff between different goals [10, 12].
Modelbased optimization principally suffers from a high level of uncertainties due to either questionable geological models or fluctuating oil prices [13–15]. Instead of a unique model, various models may match the real reservoir’s available data, e.g., logs, well tests, or production data. Such uncertainty should be taken into account in the proposed optimization workflow [16]. The significant difference between geological and economic uncertainties is that economic uncertainty has a dynamic nature and develops with time. In this situation, maximizing the shortterm goals becomes more important. Thus, in this condition, developing uncertain optimal plans are inevitable [17]. Robust Optimization (RO) has been proposed to consider the effects of such uncertainties [18–20]. The average value of the NPVs belong to all possible scenarios is the objective function of a RO problem. Using RO, the obtained conservative enough optimal plan reduces the risk of uncertainties.
Many multiobjective production optimization algorithms were used in previous studies. Some of them utilized gradientbased methods [21]. Some others used the derivativefree techniques for production optimization [22]. These approaches are distinct by three selection methods including; criterionbased, aggregationbased, and Paretobased. In a standard design to handle multiobjective problems, a singleobjective optimization method like aggregation selection [23] may be used to allocate weight factors to each objective function. However, selecting the weight factors is challenging since a wrong selection yields an improper solution. Furthermore, there is no guarantee to find truly nonconvex results relying on a Pareto front [24]. Regarding these challenges, the selection algorithms accompanying Pareto have been introduced as more efficient approaches. Populationbased evolution algorithms can employ the Paretoranking approach, between which the Particle Swarm Optimization (PSO) [25] is nominated because of a faster rate of convergence.
To the author’s knowledge, MOO under uncertainty, particularly in the form of Pareto front formation, is not adequately addressed in the literature. In the limited reported case, the main focus is on the field model updating, i.e., historymatching problems. Reference [26] has utilized MOPSO to the reservoir historymatching while quantifying uncertainty. In another study, using a historical dataset, it is applied Spherical Search Optimizer to predict the oil production from Tahe oil field in China [27]. Also, in production optimization researches, most of the studies have been focused on hierarchical optimization [7–9, 17, 19]. In a work by [10], an oil field modelbased multiobjective development problem was solved by performing three multiobjective methods of adjoint (weightedsum method and gradientbased method), Genetic Algorithm, and PSO. They found that MOPSO is more appropriate for the applied mediumscale problems by comparing these approaches since it converges higher related to other approaches. Furthermore, it provides efficient recovery of the Pareto front.
Recently, some researchers have been studied MOO under uncertainty in a Pareto front form [28–32]. The main focus of the published RO methods has been on incorporating geological uncertainty [18, 33]. However, the inclusion of economic uncertainties, e.g., oil price fluctuation, has not been well addressed yet in MOO problems. Considering such uncertainties is a vital key for decisionmakers to manage the risks and make reasonable decisions. As a representative work in case of oil pricerelated uncertainty, for example, reference [34] correlated data values with a changeable oil price. Although they implicitly discussed the impact of uncertain economic parameters, an unknown oil price explicitly affects the NPV in shortterm and longterm investigations. In recent years, some studies focused on Utilizing various oil price scenarios to perform RO in Enhanced Oil Recovery (EOR) problems as a technique in reservoir risk management [9, 35], by considering the longterm and shortterm objectives separately. In addition, another study proposed a robust design for oil production optimization under demand and market uncertainty [20].
This study has explained a process for intelligent production control of a reservoir employing a modified version of MOPSO to maximize the longterm and shortterm profit simultaneously while considering various reservoirrelated uncertainties. The archive controller has been changed, and a dynamic boundary search has been applied as the main modifications to the MOPSO algorithm. Moreover, the selecting approach of leaders has been modified by taking into account the crowding distance concept. These changes have yielded a complete Pareto front while satisfying diversity and uniformity. The modified MOPSO has been employed to find a nondominated answers collection, i.e., the points lying on a Pareto front. In addition, being needless of access to simulator source code and having a high convergence rate make the MOPSO a good optimizer for multiobjective optimizations.
LNPV as the first objective function represents the longterm gains. Since economic uncertainty risk grows with time, a highly discounted NPV (25%) is considered as shortterm gains, SNPV. The proposed approach has been practiced on an eggshaped standard model. To characterize the geological uncertainty space, ensemble members of 100 model realizations and 10 selected ones through the kmeans approach have been used. The obtained results show that the produced RO design can help achieve a fast solution to financially manage the reservoir, enabling the reservoir engineer to solve the problem with a higher conservative viewpoint concerning geological models’ uncertainties.
Furthermore, in this study, we have exercised a framework to optimize the economic goals under oil price uncertainty applying MOPSO. The improved Paretobased MOPSO approach has been applied to maximize the shortterm and longterm economic goals, robust to oil price uncertainties. This approach’s application has also been verified by performing it on the Egg model as an evaluation. Several oilprice functions have been considered to form the economic uncertainty space. By optimizing the average SNPV and LNPV of all oil price scenarios, the advised production scenario is obtained to be robust to oilprice fluctuations. The results in both optimization cases under uncertainty illustrate that MOPSO appropriately provides logical Pareto fronts resulting in optimal LNPV and SNPV pairs. It provides the option of choosing the optimal production strategy between various available scenarios under either geological or oil pricerelated uncertainties.
2 Optimization problem for reservoir management
Field existing criteria and production constraints identify the field’s shortterm objectives that improve decisionmaking in the reservoir lifecycle. Optimizing such criteria may conflict with the ultimate preference operating scenario [29]. It is generally intended to find a solution set in which a significant promotion in the shortterm objectives yields a small reduction to the field’s lifecycle target [7, 19]. To start the analysis, a common statement of the NPV as a field intended objective function, which is under waterflooding process, provided by [29], is called:
in which, the prices and rates are respectively represented by pr [in USD/m^{3}] and Q [in m^{3}/day]. d [1/year], is a parameter representing a discounting factor. Δt_{k} and t_{k} [day] are difference of timesteps and the cumulative time after timestep k, respectively. A fixed value of 365 days is considered for discounting reference time, τ_{ref}. The values of the parameters for each problem are also presented in Table 1. The values of oil/water/gas rate (Q) for calculating the objective functions are obtained from reservoir simulation in each run.
In equation (1), the longterm and shortterm objective functions are achieved by setting d equal to respectively zero (d = 0) and considerable discount rate, e.g., 25% annually (d = 0.25). The first is an undiscounted version of NPV, i.e., LNPV that is commonly utilized as the first target to optimize production operation in field management problems [7]. Ideally, it is favorable to regain the primary capital in any field development project as quickly as feasible. That is why it is required to optimize the shortterm objective function, i.e., SNPV, concurrently with maximizing LNPV. The assumed discount rate corresponding with SNPV significantly considers the production in the field’s shorttarget points of view compared with the project’s future, having no economic meaning. However, this stimulates the shortterm achievement of implementing the production scenario [5]. Considering the abovementioned cases of short and longterm targets, the following two objectives are defined:
The longterm achievement: LNPV (d = 0 in Eq. (1)).
The shortterm achievement: SNPV (d = 0.25 in Eq. (1)).
To handle the optimality concept simultaneously for LNPV and SNPV in the multicriterion optimization problem of waterflooding in a reservoir, these two objective functions are considered incompatible purposes; the achievement of one may conflict with the other. The planned optimization problem’s main intention is to perform an optimal approach that contemplates the tradeoff among these objectives since growing one purpose makes another deteriorate. The original form of MOO is employed to be utilized in the waterflooding problem, as:
J and X are respectively the objective function, including f contradicting criteria, and the vector of controlling parameters, comprising n_{p} parameters. The control variables are intended to be defined so they satisfy the h vector, consisting of y timedependent problem limitations. These constraints should be defined by case. To obtain a set of results that tradeoff among the contradicting purposes, the Paretodominance principles [5, 36] can be performed. The related concepts to apply this theory are described in the literature in detail [26, 37–39]. For instance, the dominance concept addresses a situation between X_{1} and X_{2}, X_{2} < X_{1}, where if and only if: 1) X_{2} is not better than X_{1} in all J_{M}(X), for each M∈[1, …, f], and 2) X_{1} is notably higher than X_{2} in J_{M}(X), at the minimum one M∈[1, …, f]. The set of solutions is obtained as a Paretooptimized when it is nondominated, respecting the achievable range. Corresponding objective function sets obtain the Paretooptimized front. Evolutionary algorithms, like MOPSO because of its convergence rate to reliable solutions, are suitable candidates to handle the described MOO problem. Regarding various available uncertainties, e.g., the geological or marketrelated uncertainties, RO techniques can be modified and utilized in conjunction with the optimization algorithm. This can be performed in a framework specified by “MOO under uncertainty”.
3 Optimization procedure based on the MOPSO
Moore and Chapman [40] defined the utilization of particle swarm to multicriterion optimization problems for the first time. This technique was mentioned then as MOPSO by [41] and improved in various fields of studies [25, 26, 42, 43]. In a standard form of this algorithm developed by [25], the optimization problem results are obtained by applying the Paretodominance concept. A mutation operator is utilized that benefits the search scheme. Besides that, to guide future generations’ search, the nondominated outputs will be saved in a secondary population reservation. Based on the Paretoranking approach, MOPSO is exercised according to the illustrated steps in Figure 1 to handle MOO problems. This figure represents a summary of the detailed algorithm provided in [25, 44].
Fig. 1 Original implementing approach for MOPSO. 
In this algorithm, the best positions that are met by particles are stored; thus, all the previously obtained nondominated results are saved. Two steps in this method, namely “Secondary Population generation” and “Hypercube generation”, are utilized to choose/discard repository segments. They are supposed significant steps in the external repository generated by performing the algorithm. Secondary population generation is an archive controller, and hypercube generation is a scheme of adjusting the grid to include the points lying out of the current grid’s boundary in each iteration [25, 44]. Indeed, it yields adjusting the solutions’ boundary. The mutation step shows the mutation operator, using which, by adding some mutation (unconsciousness), leads to a change in the particles’ flight and makes the high amplitude of search coverage. When the approach is trapped in a local optimum due to the particles’ zero velocities, the mutation effect is of great importance to revise the search coverage [26, 44]. Repository updating is implemented in the next step. This approach includes the leader selection operator that is addressed as a scheme to define the positions that, saving them for checking the Paretodominance, is essential in the original version of MOPSO. In the condition that two stations have no superiority over each other, one may prefer to select randomly between them to obtain the leader [25] or select one of them based on the density determination. The latter one is referenced as a scheme to select the leader based on enhancing the crowd’s particle variety and adjacency. For example, in this class, two methods based on the neighborhoods of particles include sorting the nondominated results based on the crowding distance descending sequence by [26] and sorting them respecting niche count ascending sequence by [38]. Both of these approaches can be utilized as alternatives for the random selection of the leaders.
We utilized an adapted version of the MOPSO algorithm for the waterflooded optimization problem, dealing with MOO more efficiently, as given in Algorithm A (Appendix A). In this approach, an improvement proposed by [44] is employed that includes performing a mechanism based on the changing boundary to amplitude exploration space to find new solutions. It is used to guarantee to keep the population diversity during the algorithm’s entire running without losing the exploitation characteristic. By implementing this approach, the particle boundaries exploration is appropriately increasing/decreasing, yielding a search/utilization balance in the exploration scheme. Moreover, in the employed improved method, while the archiving scheme is modified to decrease the archive controller computational cost, the Pareto front’s diversity is satisfying yet. The implemented adjustments have been explained in [44] in detail. More details about the utilized approach are presented in Appendix A.
4 Clustering and realization selection
Production optimization under uncertainty requires a long run time, especially in the case of using the revolutionary approaches, e.g., the MOPSO, together with commercial simulators. To reduce the consumption time in optimization under uncertainty, employing a clustering approach to choose representative realizations subsets is unavoidable. Various methods have been investigated and utilized in this context in the literature. For instance, using flow characteristics in the method of clustering by kmedoids [45], implementing a starting control plan to rank the obtained NPVs [32], employing some static parameters based on simulation to perform kmeans clustering [46], and a combination strategy of permeabilitybased and flowbased quantities [47] are introduced previously in the literature. We select to use the kmeans algorithm to select a realizations subset representative of all geological models. Having an uncomplicated concept, simple implementation, and reasonable computations are the reasons for this preference [48, 49]. However, to ensure this approach’s results are qualified enough, we have checked the results of the selected representative realizations subset regarding features like those of the obtained NPV’s and the results of simulations [50, 51].
The algorithm of kmeans clustering is based on a prototype distribution clustering idea. In this technique, one can find the userspecified number of clusters, which are represented by their centroids (means). If a point is closer to a cluster’s centroid, it is considered to be in that particular cluster. With a k number initial set of means, kmeans method obtains the most reliable centroids by shifting among (1) specifying data points (x_{p}) (here, grids’ permeability) to clusters (s^{(t)}) on the basis of the present centroids whose mean is calculated by the least squared of Euclidean distance, which is intuitively the nearest mean. It can be represented by:
and, (2) choosing centroids based on the current assignment of data points to clusters, thus, calculating the new means of the data points in the new clusters, i.e.,:
The algorithm has converged when no change in the assignments has appeared. Therefore, the objective function to meet the final result is the minimization of the squared distance sum between each point and the group’s centroid. It is similar to the average dot product maximization, V[i] · V[j], if V[i] and V[j] belong to the same cluster, for unit vectors, V. Figure 2 schematically illustrates the process. There are several advantages for the kmeans method, from which fast computation, excellent performance on large datasets, and needing only one input parameter (k) are of great importance [52].
Fig. 2 Schematic of kmeans algorithm, dots are training examples, and triangles are cluster centroids. (a) Original dataset, (b) random initial cluster centroids, (c–f) illustration of running two iterations of kmeans. In each iteration, each training example is assigned to the closest cluster centroid; then each cluster centroid is moved to the mean of the points assigned to it. 
5 Robust optimization
In an ideal case, nominal optimization is utilized for production optimization. However, in the presence of uncertainties, not only there is not a guarantee that such a solution is optimum, it may be infeasible, too. In this respect, Van Essen et al. [18] implemented RO by considering model uncertainty in the optimization framework. RO is a conservative approach to deal with uncertainty, especially in green fields at the start of reservoir life. Using an ensemble of geological or economic model realizations, one can explicitly quantify the uncertainty. In this study, in the RO framework, objective function, also known as robust objective or mean optimization, is defined as the average of NPV over the ensemble of geological/economic models [18]:
where, J_{RO} is the robust objective function, N_{e} is the number of realizations, and J_{i} is the calculated objective function of ith realization.
5.1 Case study
The adjusted MOPSO in the form shown in Algorithm A (Appendix A) is performed on a standard oil reservoir model experiencing waterflooding process, namely Egg model. This benchmark has been previously applied for MOO problems in some works reported in the literature. As the first effort, without considering any uncertainty, we practice the adjusted MOPSO for simultaneous maximization of SNPV and LNPV to show the effectiveness of the approach compared to previous works in the literature. Then, RO is applied for handling geological uncertainty by implementing MOPSO for the mentioned MOO problem. SNPV and LNPV are desired to be maximum, while 100 realizations of the geological model are respected. The number of realizations is shown to be reduced in an effective way to decrease the required computational time, without a significant missing the accuracy of the solutions. Also in another investigation, oil price uncertainty is considered and analyzed based on the proposed approach results. In all cases, wells’ rates of injection are considered control parameters in the investigated optimization problems. The proposed MOPSO illustrated in Algorithm A (Appendix A) is completely compatible with a conventional reservoir simulator that is necessary for NPV computation at the optimization period. A zero rate of discount is assigned in equation (1) for LNPV calculation. The LNPV is an indicator of the cash flow in the lifecycle of the reservoir in cumulative form. For calculation of SNPV, however, a value of 0.25 is assigned for the rate of discount, to make it more significant than the LNPV which reflects the future exploitation. Since regaining the initial funding in each development plan as early as possible, is desirable, it is ideal to satisfy shortterm purposes simultaneously with meeting the longterm objective function. In this regard, we may select a high annual discount rate of 25% (b = 0.25) in equation (1) to reach SNPV that accentuates the significance of maximizing shortterm economic goals. Although NPV with a 25% discount rate may have no meaning in the economic sense, it weights the shortterm production remarkably more than the production in the following years.
5.2 Egg model
The egg model introduced by [18] is an oil reservoir model that is eggshaped and threedimensional. It composes 25 200 blocks with 18 533 ones being active that indeed build that eggshaped structure. Some channels belong to twelve wells, consisting of eight injectors and four producers, which are set in the model. The locations of the wells are depicted in Figure 3.
Fig. 3 Eggshaped model depicting the position of injectors and producers. 
Injection rates of injectors are desired to be adjusted in a way that maximizes both LNPV and SNPV in a related MOO problem. The range for these optimization parameters is set to be 0–79.5 m^{3}/day. Every 900 days are considered one optimization interval; four of them are required to cover the whole simulation run of 3600 days. On the sum, eight injectors multiplied by four timesteps means that 32 variables are defined to assign optimum value through performing the algorithm of optimization. Table 1 summarized the fluid properties, operational conditions, and geological details of the eggshaped model, as mentioned in [53]. Also, the values of oil and operational prices used in equation (1), are listed in Table 1.
5.3 Modelbased optimization: single realization
To verify the adapted form of MOPSO method as described in Algorithm A (Appendix A), first, the optimization algorithm has been applied on the base case Egg model in the case of considering no uncertainty. The applied parameters to run the algorithm is defined in Table 2.
Optimizer’ parameters of MOPSO – Base case.
It is evident from Figure 4 that MOPSO well makes the Pareto fronts in various iterations and in relatively higher ones, it spans the same range that is specified by the Utopia line. The Utopia line is obtained based on running a routine singleobjective optimization approach to maximize LNPV and SNPV, separately [21]. The search space, together with the ultimate optimum Pareto front, obtained by applying the modified MOPSO, is illustrated in Figure 5. From the diagram, the Pareto front’s final result, with a desired concave shape, is nearly formed the whole area of expected results, previously demonstrated by the Utopia line. In the case of using the approach of random deletion of leaders in the repository in the MOPSO, the diversity of solutions is not met, especially when the method has a high convergence rate. By considering the crowding distance in the deletion process, the Pareto front’s diversity is kept, and it is formed concave. Therefore, here a concavedown shape is desired, representing a complete Pareto front.
Fig. 4 Paretofronts of various iterations achieved through adjusted MOPSO – Basecase. 
Fig. 5 Pareto front, and the search space of particles, after 100 iterations of adjusted MOPSO. 
The obtained results are consistent with previous works by [2, 7], showing maximization of the shortterm objectives is permitted by an additional degree of freedom through optimization of reservoir lifetime. For this problem, to increase SNPV by about 33%, it is required to only decrease LNPV by about 10% that is an appropriate output based on their analyses’ results. The number of required objective function calls for this algorithm to converge to a reasonable solution set is compared with three previous methods practiced by [54], as shown in Table 3. These methods are consisting of: Original NBI method, NBI tracking method, and Adjusted Weighted Sum method, the details of which are defined in [21]. The number of necessary simulation runs for MOPSO is relatively low in comparison with most of them. For a decisionmaker, the developed Pareto front gives an effective scheme that makes it possible to tradeoff between practical water injection plans related to the different situations, but it also allows to deal with various constraints for goals of a short or long period of times.
Comparison of the number of simulation calls required for various methods – base case Egg model.
5.4 Production optimization under uncertainty
The geological model’s uncertainty is always present because the data from core samples or well logs and seismic outputs may contain an error. On the other hand, economic parameters such as oil prices, interest rates, political risks, and unpredictable operational conditions make the uncertainties inevitable. One way to deal with these issues is to use an ensemble of geological/economic realizations, and utilizing the mean of objective functions over them leads to a robust solution [18]. To solve the optimization problem in a multiobjective scheme, we consider the average of SNPVs and LNPVs, belong to all realizations in the ensemble, as objective functions that are robustly related to uncertainties. We considered two challenging types of uncertainty. The first one is the geological model uncertainty that focuses on uncertain permeability maps in this problem. It is quantified by considering 100 (and a reduced number of 10) model realizations and directly affects production/injection plan and indirectly on the objective function (NPV). The second one is the oil price uncertainty which is quantified by regarding 10 different oil price functions. This parameter directly and explicitly affects the NPV calculation. In the last section of the paper, also both uncertainties are considered. Therefore, the effect of uncertainty on this problem is that a range of solutions (the Pareto front) is obtained with a different set of controlling parameters (injection rate).
5.4.1 Optimization under geological uncertainty – RO
To respect the geological uncertainty in Egg model, 100 different model realizations, obtained by 100 different distributions of absolute permeability, are applied. The economic parameters including oil price r_{o}, and water production cost r_{w}, are chosen as the same as the data in Table 1. An average NPV can be determined by respecting all the geological realizations. It is the basis for the RO approach, as the objective function that is intended to be maximized:
where N_{geo} is the number of geological realizations and NPV^{i} is the corresponding objective function to the ith realization. The applied parameters to run the algorithm for this problem are defined in Table 4. It should be noted that, using the optimizer code and the simulator in a 64 bit PC with an Intel Core i7 processor and 16 GB RAM running Windows 7, each iteration lasts about 9000 s.
The MOPSO optimizer variables for production optimization under geological uncertainty.
Figure 6 shows the Pareto fronts obtained for each model realization. Also, the Pareto front resulted from RO is represented by the dotted line in this figure. The Pareto front obtained from averaging on all realizations is robust versus geological uncertainty. Therefore, the decisionmaker can utilize this Pareto front resulted from robust optimization as a conservative scenario to select the optimum injection plan. Such a strategy plays a significant role in the decisionmaking process by reducing the risks from reservoir model uncertainty. As it is clear from this figure, the Pareto front incline in the diagram is sharp near the maximum LNPV, which states that a small change in LNPV leads to a considerable change in SNPV. It has the advantage of making it possible to increase SNPV to a satisfying value while missing only a negligible performance related to the field’s lifecycle target.
Fig. 6 The obtained Pareto fronts using the ensemble of 100 realizations; the dotted line is the robust Pareto front. 
5.4.2 RO with reduced number of realizations
In this work, uncertainty has been considered for the permeability map while other system parameters, like porosity and NettoGross ratio, are considered to be known and constant. In order to use the kmeans for clustering the 100 realizations, we have employed permeability information as feature vectors, captured by permeability distribution in 25 200 grid blocks of the model. However, in the case of having a large number of realizations, Principal Component Analysis (PCA) representation of fields of permeability can be used as proposed by [47] to display these fields (geology) with respect to a few features.
Our goal for permeability clustering is choosing several subsets of realizations that represent the original 100 sets of model ones. The computation of the distance matrix is performed regarding the Euclidean distances between the generated vectors of feature through the kmeans algorithm. Then, dividing the realizations into 10 groups (clusters) by the algorithm is performed by minimizing the sum of square distances belonging to withincluster.
In Figure 7, 10 field permeability realizations, selected through performing the kmeans method, are presented.
Fig. 7 Ten permeability maps selected by the kmeans approach. 
The Pareto fronts resulted from implementing the optimization approach on the 10 selected models by the kmeans method are represented in Figure 8. The same parameters of the algorithm as mentioned in Table 4 are applied for practicing kmeans approach. The robust Pareto front of these models is also depicted, showing an acceptable consistency with the selected realizations.
Fig. 8 the obtained Pareto fronts using the ensemble of 10 selected realizations by kmeans method; the dottedline is RO Pareto front. 
5.4.3 Comparison of using the reduced realizations with all realizations
To check the precision of using the static parameter for clustering (reservoir permeability) with flow responses, a detailed comparison of the deviation of the two Pareto fronts, the diagram of Cumulative Density Function (CDF) for both long and shortterm income of robust solutions are presented in Figure 9, separately. It is clear that after 50 iterations, the obtained CDFs of both SNPV and LNPV of the 10 selected realizations are good estimations for the ones of the whole realizations. The outputs indicate that by utilizing a reasonable clustering method, one can optimize the problem to achieve a conservative estimation in minimum time from the proposed injection strategy.
Fig. 9 Comparison of the CDFs for SNPV (a) and LNPV (b) after 50 iterations for 10 selected realizations with that of 100 realizations for production optimization of the Egg model. 
As it is stated before, reservoir heterogeneity plays an important role in waterflooding performance. In the previous part, we implement the reservoir permeability as a static clustering parameter. For another comparison, in this part, the Recovery Factor (RF) of the 10 selected realizations is calculated and compared with those of all 100 realizations. It should be noted that the recovery factor is a flow response affected by permeability. To calculate the RF, we used three injection profiles: i) using the maximum allowable rate of injection (79.5 m^{3}/day), ii) using the injection profile which belongs to the maximum SNPV in Pareto front of robust optimization of 100 realizations (Fig. 6), iii) using the injection profile which belongs to the maximum LNPV in Pareto front of robust optimization of 100 realizations (Fig. 6). In Figure 10, the CDF diagram of RF belong to the three mentioned cases is depicted with their corresponding 10selected realizations based on RF. As it is clear, a relative consistency existed between the samples and ensemble from an RF point of view. It should be noted that, because the high injection rate in the first case of Figure 10 decreases the effect of permeability and heterogeneity, the range of RF is narrow and the CDF diagram partially deviates from the CDF which belongs to all realizations.
Fig. 10 Comparison of the CDFs for RF of all 100 realizations and their corresponding 10selected realizations resulted from (a) using the maximum allowable injection rate, (b) using injection profile of best LNPV in Figure 7, (c) using injection profile of best SNPV in Figure 6. 
A graphical comparison of the Pareto fronts obtained from the robust problem optimization based on 10 and 100 realizations is depicted in Figure 11. For a comparison to show the applicability of the selected models based on permeability, we have used the RF of all 100 realizations as a dynamic clustering parameter to include the effect of sweep efficiency. To select the new clusters, the same injection profiles as the last part is used for flow simulation. Ten realizations from each injection profile are selected by performing the kmeans method using the RF as a clustering parameter. Figure 11 illustrates the comparison of robust optimization of the newly selected realizations with the previous Pareto fronts. As shown in this figure, the Pareto front belongs to 100 realizations, and that of 10 realizations obtained from static clustering is relatively consistent; the Pareto front from the maximum injection rate scenario has highly deviated, as expected. It means that using the static clustering parameter (permeability) acceptably reflects the waterflooding performance in this case study. As shown in this figure, regarding the computational time that may affect the optimality of problemsolving, one can replace the whole realizationsrelated robust optimal plan with one derived only from the clustered models. It leads to acceptable results that are comparable with the outputs obtained from the whole realizations.
Fig. 11 Comparison of the obtained Pareto front after 50 iterations of robust optimization on different set of realizations of the Egg model (the definitions of RF cases: (a) using the maximum allowable injection rate, (b) using injection profile of best LNPV in Fig. 7, (c) using injection profile of best SNPV in Fig. 6). 
5.4.4 Production optimization under economic uncertainty
Economic parameters used in NPV formulation in equation (1), e.g., oil price and interest rate, are not predictable and always change with time. Since conventional petroleum reservoirs’ lifecycle spans 10 to 100 years, these parameters’ changes are unknown and generate a new source of uncertainty. For instance, in the current year, the Coronavirus pandemic drastically reduces the oil price by 50–80%.
As mentioned before, the inclusion of economic uncertainties, e.g., oil pricerelated issues, has not been well addressed in the petroleum engineering literature of MOO problems. In this section, the oil price uncertainty is explicitly included in the robust MOO framework to form the Pareto front [55]. A finite number of scenarios can represent varying oil prices. Other economic parameters, including waterinjection/production costs, are considered to be fixed.
Oil price scenarios
Prediction of the expected energy prices in the future is performed through different approaches. For example, the European Union and the French government use Prospective Outlook on Longterm Energy Systems (POLES) [56–57]. Siraj et al. utilized a simplified AutoRegressiveMovingAverage model (ARMA) [58] to make oil price timeseries and represent the impact of fluctuating oil prices [14]. We also use this model to consider oil price uncertainty. The ARMA model is defined as:
where ω_{k} is a whitenoise sequence, α_{i} and β_{i} are the randomly selected coefficients [8]. In an approach by [17], 10 scenarios were generated using 471 USD/m^{3} as base oil price of, as shown in Figure 12.
Optimization under oil price uncertainty – robust scenario
In this section, similar to optimization under geological uncertainty, RO under oil price uncertainty has been performed. In this case, it is assumed that the true model is the base case realization (Fig. 3), and different oil price scenarios were considered to obtain the conservative Pareto front. To do this, SNPV and LNPV for each price scenario are calculated using the injection profile. Then the average NPV is considered over the ensemble set of oil price scenarios and is maximized as an objective function. The average NPV defined over the ensemble of oil prices is obtained by replacing N_{geo} in equation (6) with N_{eco} as the economic realizations:
The applied parameters of the algorithm are mentioned in Table 5.
The MOPSO optimizer variables for production optimization under economic uncertainty.
In Figure 13, the results of optimization corresponding to each oil price scenario are illustrated in the form of the Pareto front to compare them with the result of RO. Also, the Pareto front belongs to the meanvalue of oil price and that of RO are obtained. As it is clear, the Pareto front resulted from RO is acceptably similar to the optimization based on the oil price meanvalue. Any point on the robust/mean Pareto front gives the injection procedure corresponding to the selected LNPVSNPV pairs, which can be considered a lowrisk scenario related to economic uncertainties.
Fig. 13 The obtained Pareto fronts using the ensemble of 10 oil price scenarios, mean oil price, and RO. 
Discussion: comparison between optimum injection scenarios
Figure 14a shows the graph of NPV versus optimization time, depicted for both the best LNPV and SNPV points belong to the Pareto front of the meanvalue oil price. The best LNPV means the maximum obtained LNPV in the Pareto front (the point on the right side of the graph), and the best SNPV means the maximum obtained SNPV in the Pareto front (the point on the left side of the graph), and their corresponding injection profile. As represented in this figure, in the best SNPV scenario, production from the reservoir after four years is no longer economical. In that time, due to the increasing costs of water injection and production, it is needed to reoptimize the injection scenario. In that case, continuing the reservoir production with the current condition and performed optimization is costeffective only by increasing the oil price prospective. At the same time, the best longterm scenario is a conservative one that makes the production plan costeffective for a minimum of seven years. The increasing oil price perspective enhances the analysis control. These analyses show the main contributions in the decisionmaking process, which include using the oil price meanvalue scenario, assigning the optimum scenario based on economic, technical, and operational assumptions, and updating it for any required time.
Fig. 14 NPV changes with simulation time for a hypothesis case of using the injection profiles belong to (a) meanvalue scenario and (b) robust scenario while both are utilizing the meanvalue oil price to generate objective functions through time. 
In Figure 14b, the previous diagram is reproduced for the injection profile resulted from RO. To generate this figure, the robust scenario’s injection profile and the meanvalue oil price are utilized. As it is evident, using this injection scheme, the charts for the best SNPV show a more sensible form for stable reservoir production. Using this scenario, in the case of improvement of future oil prices, one can hope to increase profitability. However, if the oil price decreases, the scenario is conservative enough to cover longterm profitability.
5.4.5 Optimization under both geological and economic uncertainty–robust scenario
This section considers both geological and oil price uncertainty, and RO under these uncertainties have been performed. We applied 10 oil price scenarios on the 10 model realization selected by kmeans with permeability as a clustering parameter. The average NPV is considered over the ensemble set of oil price and geological scenarios and maximized as an objective function. The average NPV defined over these ensembles is obtained by considering N_{geo} and N_{eco} as the geological and economic realizations as:
The applied parameters of the algorithm are mentioned in Table 5. For performing optimization under uncertainty, considering 10 oil price functions, 10 geological realizations (those selected by kmeans approach), a population of 30, and 30 maximum iterations, 90 000 simulations are needed.
In Figure 15, the optimization results corresponding to applying the oil price scenarios on each geological realization and the robust scenario regarding all uncertainties are illustrated in the form of the Pareto front. Each Pareto front in this figure is robust to economic uncertainty. The red one (the robust Pareto front) is robust to coupled economic and geological uncertainties and therefore is the most conservative optimum scenario. Any point on the robust Pareto front gives the injection procedure corresponding to the selected LNPVSNPV pairs, which can be considered a lowrisk scenario related to both geological and economic uncertainties.
Fig. 15 The obtained Pareto fronts from RO using the ensemble of 10 oil price scenarios on 10 selected realizations. 
Figure 16 depicts the proposed injection scenarios for the two endpoints of the obtained Pareto fronts in various analyses performed on the case study. This may provide useful information for the decisionmaking process respecting the comparison of different cases. For instance, the comparison of the best SNPV point respecting the 100 geological realizations with that point obtained by 10 geological ones shows that injection rates are higher for most of the wells, in the former case. Comparison of these two cases’ Pareto fronts on the other hand (Fig. 9), illustrates that their two values on the best SNPV points are not so far from each other, but their corresponding LNPV points are. It means that by increasing water injection rates from a scenario similar to case 3 to a one similar to case 2 (Fig. 16) not only yields ignorable improvements in SNPV, it also leads to missing considerable LNPV because of a longtime increasing in water injection rates. As another example, assume one decides to gain a maximum LNPV according to the Pareto front of either meanvalue optimization or RO that are approximately the same as shown in Figure 13. Thus, the provided analyses propose two different scenarios for the decisionmaker as depicted by case 5 and case 6 in Figure 16. It is defined that case 5 may be a better choice since it provides more slight and stable injection rates in comparison with the aggressive injection profiles in case 6. It is worth mentioning that the corresponding SNPV in the maximum LNPV is approximately the same for both cases as shown in Figure 13.
Fig. 16 Optimized injection scenarios [m^{3}/day] in Egg model obtained by MOPSO bluecolored lines: the best LNPVs, and redcolored lines: the best SNPVs. 
6 Conclusion
Geological and economic uncertainties in modelbased optimization are the main difficulties in maximizing shortterm financial goals. In this study, using a modified multiobjective approach in the optimization process, namely MOPSO, with explicit consideration of such uncertainties (i.e., RO) resulted in the robust prediction of longterm economic goals while allowing the right balance with shortterm objectives. In the process of selection/deletion of leaders in the repository of MOPSO, the crowding distance calculation was included. Also, the changing boundary exploration was employed in the MOPSO approach. Such modifications facilitate the optimizer to produce a complete Pareto front having an acceptable convergence speed and being needless of access to simulator source code. A benchmark eggshaped model was utilized to show the effectiveness of such an adjusted form of MOPSO. From various investigations on that model, the following conclusions can be deduced:
A practice of MOPSO on the model without considering any uncertainty, as a basecase, and to simultaneously update LNPV and SNPV, shows its effectiveness in comparison with some previous similar efforts via different algorithms in the literature.
Using the RO while performing MOPSO under the geological model uncertainty with 100 model realizations yields a Pareto front showing that the SNPV considerably increases without discrediting the LNPV. Also, utilizing the kmeans method, 10 geological model realizations were selected that are representatives of those 100 original ones. Results from optimization under uncertainty using these 10 model realizations show that the formed Pareto front is acceptably applicable as the problem solution. This approach yields decreasing the required computational time since it considerably reduces the number of simulator calls through running the optimization algorithm in the performed RO.
Respective of technical, political, and operational changes that affect future oil prices, attention to economic uncertainty is necessary. Considering the economic uncertainty in robust MOO explicitly, using 10 oil price profiles, the robust Pareto front with taking into account oil price fluctuations is derived. Then, employing the corresponding injection scenario introduces a conservative plan, which can cover the economic risks. Therefore, like the optimization under geological uncertainty, the proposed optimization method provides a framework that helps a decisionmaker select the best production scenario considering the oil price uncertainty.
Using other controls, such as geological or inversion parameters as an optimization variable is one aspect that one can consider in future work. Also, considering alternative geological realizations as an inversion parameter to use seismic data in history matching problems maybe an interesting idea to search. Moreover, to work with real reservoir data, utilizing the proxy models in production optimization is proposed to reduce the computational cost. Furthermore, using other optimization methods such as StoSAG or gradientbased optimization techniques is one aspect to consider in future studies to compare the efficiency of different methods from computational cost point of view.
Appendix A
We utilized an adapted version of the MOPSO algorithm for the waterflooded optimization problem, dealing with MOO more efficiently, as given in Algorithm A (Appendix A). In this approach, an improvement proposed by [44] is employed that includes performing a mechanism based on the changing boundary to amplitude exploration space to find new solutions. It is used to guarantee to keep the population diversity during the algorithm’s entire running without losing the exploitation characteristic. By implementing this approach, the particle boundaries exploration is appropriately increasing/decreasing, yielding a search/utilization balance in the exploration scheme. Moreover, in the employed improved method, while the archiving scheme is modified to decrease the archive controller computational cost, the Pareto front’s diversity is satisfying yet. The implemented adjustments have been explained in [44] in detail.
Algorithm A The adopted MOPSO algorithm improvement utilized for operation optimization in waterflooding problems. 
To start evaluating process in the initialization step, a random population (X_{i}), which includes controlling parameters within its boundary (B), is chosen:
Now, the objectives (SNPV/LNPV) are calculated using equation (1) and the best answers are stored. In the archive controller step (the 4^{th} step), generating hypercubes belong to the search space in a way that each particle is located by means of a coordinate system namely the hypercubes, according to its goal function value, and utilized to select external archive (repository) members through segmentation. Regarding all the particles are not saturated, the initial saturation count is set at zero. Zain et al. [44] fully described the utilized modifications in detail. After finding the best initial positions, in each succeeding iteration firstly new boundary is generated:
Then offspring population is constructed inside the new boundary by:
The next step is to determine the probability of each variable to go to its boundary by:
Now, respecting the new population evaluation of the objective functions (SNPV and LNPV), storage of the nondominated ones, updating repository by deleting dominated solutions, and updating hypercubes are performed, respectively. Updating memory of particles occurred by replacing the previous best location with the current best location of each particle. Two main components of the repository employed by MOPSO are the archive controller and the grid. The first is a decisionmaker tool for the addition/deletion of solutions. The second presents the objective space partitions to save the answers. This method is useful to decrease the computational cost; meanwhile, the archive controller requires adding or reducing the repository member. After performing mutation and finally selecting the best positions, the current boundary should be shrunk by shrink factor, as:
Then the BF is initialized by a random value, C5. If the objective function evaluation be more than half of its predetermined limit, the boundary factor is calculated as:
where C6 is a random value. It should be noted that we consider the value of half of repository size, REP.SIZE/2; if for REP.SIZE/2 iteration(s) the members in repository reach to maximum repository size, we consider the BF as C5, again, and recalculate the boundary factor. Moreover, if REP.SIZE/2 iteration(s) equals to saturation factor, we calculate the BF with equation (A6). Finally, we update the boundary by if:
In these formulas, B is the boundary in MOPSO algorithm, B_{i} is the lower and higher boundaries absolute difference divided by two, BF is boundary factor in MOPSO algorithm, X_{i} is initial population, X″ is population after performing mutation, X_{max} and X_{min} are higher and lower boundary of control variable respectively. Random values of C1–C6 in this algorithm are sampled from [0,1], and the constant value of A is set equal to 0.1. Two factors, namely Saturation Factor (SaF) and Shrink Factor (SF), with default values of 10 and 0.97 [44] are responsible for controlling the algorithm convergence trend. SaF that is a measure to define if a particle of a crowd considered saturated, is employed in Step #18, the procedure of updating population; and in the 23rd step, the measurement of boundary factor. The rate of shrinking for each particle’s exploration boundary is represented by SF in the MOPSO algorithm respecting each dimension. Mutation replaces the old population in this algorithm. Times of implementing the mutation procedure, however, is changed in the modified algorithm related to the original MOPSO; i.e., it is implemented only on the old population, as a substitute for employing it following the new population generation [44]. Additionally, in leader selection at Step #11 (Algorithm A), the crowding distance approach in a manner practiced by [26] has been employed as another technique of modification. Considering a specific point, i, the two points on either side of point i are selected. The crowding distance is the average of the distances between each of these points and i, and is computed for each objective. Areas with a larger crowding distance are particularly favored for the selection of a local leader in particular.
Each time new members are admitted, the archive controller in the original MOPSO entails determining each repository member’s domination. This process yields high computational costs, mainly during the repository size becomes more extensive in any iteration. Additionally, similar answers or particles in the group may be admitted into the finitesized repository, leading to quick overfilling. When new nondominated solutions are obtained in the modified approach, the archive controller will add them to the repository. Nonetheless, if the repository size surpasses the maximum permitted, some repository members will be dismissed by the archive controller [44].
In the original approach, the grid density decides which member to be excluded – the higher the density, the more chance of members to be excluded. In the modified approach, the parameter which controls the replacement depends on the Euclidean distance in the objective space between any member of the repository and the newest entitled member. In the applied approach, the old population is replaced utilizing mutation. The standard mutation in the original MOPSO is performed with a change in the time of its execution [44]. Therefore, it happens just on the old population alternately after making the new population, as described in Algorithm A. More details of the mutation process are defined in Zain et al. [44].
References
 Sarma P., Aziz K., Drlofsky L.J., (2005) Implementation of adjoint solution for optimal control of smart wells, in: SPE Reservoir Simulation Symposium Houston, Texas, USA. Paper SPE 92864. [Google Scholar]
 Jansen J.D., Douma S.D., Brouwer D.R., Van den Hof P.M.J., Bosgra O.H., Heemink A.W. (2009) Closed Loop Reservoir Management, in: SPE Reservoir Simulation Symposium, 2–4 February, The Woodlands, Texas, USA. https://doi.org/10.2118/119098MS. [Google Scholar]
 Chen Y., Oliver D.S. (2010) EnsembleBased ClosedLoop Optimization Applied to Brugge Field, SPE Reserv. Evaluation Eng. 13, 1, 56–71. https://doi.org/10.2118/118926PA. [Google Scholar]
 Jansen J.D. (2011) Adjointbased optimization of multiphase flow through porous mediaa review, Comput. Fluids. 46, 1, 40–51. [Google Scholar]
 Isebor O.J., Durlofsky L.J. (2014) Biobjective optimization for general oil field development, J. Pet. Sci. Eng. 119, 123–138. https://doi.org/10.1016/j.petrol.2014.04.021. [CrossRef] [Google Scholar]
 da Cruz Schaefer B., Sampaio M.A. (2020) Efficient workflow for optimizing intelligent well completion using production parameters in realtime, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 75, 69. [CrossRef] [Google Scholar]
 Van Essen G.M., Van den Hof P.M.J., Jansen J.D. (2011) Hierarchical longterm and shortterm production optimization, SPE J. 1, 191–199. https://doi.org/10.2118/124332PA. [CrossRef] [Google Scholar]
 Siraj M.M., Van den Hof P.M.J., Jansen J.D. (2015) Handling risk of uncertainty in modelbased production optimization: A robust hierarchical approach, 2nd IFAC Workshop Autom. Control Offshore Oil Gas Prod. Florianpolis, Brazil 48, 6, 248–253. [Google Scholar]
 Siraj M.M., Van den Hof P.M., Jansen J.D. (2016) Robust optimization of waterflooding in oil reservoirs using risk management tools, in: Proceedings of the 11th IFAC Symposium on Dynamics and Control of Process Systems, pp. 133–138. [Google Scholar]
 Fu J., Wen X.H. (2017) Modelbased MOO methods for efficient management of subsurface flow, SPE J. 22, 6, 1984–1998. https://doi.org/10.2118/182598PA. [CrossRef] [Google Scholar]
 von Hohendorff Filho J.C., Schiozer D.J. (2020) Influence of well management in the development of multiple reservoir sharing production facilities, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 75, 70. [CrossRef] [Google Scholar]
 Bagherinezhad A., Boozarjomehry Bozorgmehry R., Pishvaie M.R. (2016) Multicriterion based well placement and control in the waterflooding of naturally fractured reservoir, J. Pet. Sci. Eng. 149, 675–685. https://doi.org/10.1016/j.petrol.2016.11.013. [Google Scholar]
 Schiozer D.J., dos Santos A.A.S., Santos S.M.G., Hohendorff Filho J.C. (2019) Modelbased decision analysis applied to petroleum field development and management, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 74, 46. https://doi.org/10.2516/ogst/2019019. [CrossRef] [Google Scholar]
 Gallardo E., Deutsch C.V. (2019) Decision making in the presence of geological uncertainty with the meanvariance criterion and stochastic dominance rules, SPE Reserv. Evaluation Eng. 23, 1, 1094–6470. [Google Scholar]
 Vincent P., Schaaf T. (2019) Reservoir and economicuncertainties assessment for recoverystrategy selection using stochastic decision trees, SPE Reserv. Evaluation Eng. 22, 4, 1094–6470. [Google Scholar]
 DubosSallée N., Fourno A., ZarateRada J., Gervais V., Rasolofosaon P.N.J., Lerat O. (2020) A complete workflow applied on an oil reservoir analogue to evaluate the ability of 4D seismics to anticipate the success of a chemical enhanced oil recovery process, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 75, 18. [CrossRef] [Google Scholar]
 Siraj M.M., Van den Hof P.M.J., Jansen J.D. (2017) Handling geological and economic uncertainties in balancing shortterm and longterm objectives in waterflooding optimization, SPE J. 22, 4, 1313–1325. https://doi.org/10.2118/185954PA. [CrossRef] [Google Scholar]
 Van Essen G., Zandvliet M., Vanden Hof P.M.J., Bosgra O., Jansen J.D. (2009) Robust waterflooding optimization of multiple geological scenarios, SPE J. 14, 01, 202–210. [Google Scholar]
 Chen C., Li G., Reynolds A.C. (2012) Robust constrained optimization of short and longterm net present value for closedloop reservoir management, SPE J. 17, 3, 849–864. [Google Scholar]
 Beiranvand H., Ghazanfari M., Sahebi H., Pishvaee M.S. (2018) A robust crude oil supply chain design under uncertain demand and market price: A case study, Oil Gas Sci. Technol.– Revue d’IFP Energies nouvelles 73, 66. [Google Scholar]
 Fonseca R. (2015) A modified gradient formulation for ensemble optimization under geological uncertainty. PhD thesis, Delft University of Technology. [Google Scholar]
 Isebor O.J. (2013) Derivativefree optimization for generalized oil field development, Unpublished PhD thesis, Stanford University. [Google Scholar]
 Zitzler E., Deb K., Thiele L. (2000) Comparison of multiobjective evolutionary algorithms: Empirical results, Evol. Comput. 8, 2, 173–195. [Google Scholar]
 Das I., Dennis J.E. (1997) A closer look at drawbacks of minimizing weighted sums of objectives for Pareto set generation in multicriteria optimization problems, Struct. Optim. 14, 1, 63–69. [CrossRef] [Google Scholar]
 Coello Coello C.A., ToscanoPulido G., SalazarLechuga M. (2004) Handling multiple objectives with particle swarm optimization, IEEE Trans. Evol. Comput. 8, 3, 256–279. [CrossRef] [Google Scholar]
 Mohamed L., Christie M., Demyanov V. (2011) History matching and uncertainty quantification: Multiobjective particle swarm optimisation approach, Society of Petroleum Engineers, Vienna, Austria. https://doi.org/10.2118/143067MS. [Google Scholar]
 Alalimi A., Pan L., Alqaness M.A.A., Ewees A.A., Wang X., Abd Elaziz M. (2021) Optimized random vector functional link network to predict oil production from Tahe oil field in China, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 76, 3. [CrossRef] [Google Scholar]
 Fonseca R., Leeuwenburgh O., den Hof P.V., Jansen J. (2014) Ensemblebased hierarchical multiobjective production optimization of smart wells, Comput. Geosci.: Model. Simul. Data Anal. 18, 3–4, 449–461. [Google Scholar]
 Bouzarkouna Z., Ding D.Y., Auger A. (2012) Well placement optimization with the covariance matrix adaptation evolution strategy and metamodels, Comput Geosci. 16, 75–92. https://doi.org/10.1007/s1059601192542. [CrossRef] [Google Scholar]
 Fonseca R.M., Leeuwenburgh O., Della Rossa E., Van den Hof P.M.J., Jansen J.D. (2015) Ensemblebased multiobjective optimization of onoff control devices under geological uncertainty, SPE Reserv. Evaluation Eng. 18, 4, 1094–6470. [Google Scholar]
 Yasari E., Reza M., Khorasheh F., Salahshoor K. (2013) Application of multicriterion robust optimization in waterflooding of oil reservoir, J. Pet. Sci. Eng. 109, 1–11. https://doi.org/10.1016/j.petrol.2013.07.008. [CrossRef] [Google Scholar]
 Yasari E., Pishvaie M.R. (2015) Paretobased robust optimization of waterflooding using multiple realizations, J. Pet. Sci. Eng. 132, 18–27. https://doi.org/10.1016j.petrol.2015.04.038. [Google Scholar]
 Abellan A., Noetinger B. (2010) Optimizing Subsurface Field Data Acquisition Using Information Theory, Math Geosci. 42, 603–630. https://doi.org/10.1007/s1100401092856. [CrossRef] [Google Scholar]
 Wen T., Ciaurri D.E., Thiele M., Ye Y., Aziz K. (2014) How much is an oil price forecast worth in reservoir management?, in ECMOR XIV14th European Conference on the Mathematics of Oil Recovery, Catania, Sicily, Italy. https://doi.org/10.3997/22144609.20141900. [Google Scholar]
 Li H., Dang C., Mirbozorg A., Yang C., Nghiem L. (2019) Robust Optimization of ASP Flooding Under Oil Price Uncertainty, in: SPE Reservoir Simulation Conference, 10–11 April, Galveston, Texas, USA. [Google Scholar]
 Yu P.L. (1974) Cone convexity, cone extreme points, and nondominated solutions in decision problems with multiobjectives, J. Optim. Theory Appl. 14, 319–377. https://doi.org/10.1007/BF00932614. [CrossRef] [Google Scholar]
 Engelbrecht A. (2005) Fundamentals of Computational Swarm Intelligence, John Wiley & Sons, Chichester, England, UK. [Google Scholar]
 Reyes M., Coello C. (2006) Multiobjective particle swarm optimisers: A survey of the stateoftheart, Int. J. Comput. Intell. Res. 2, 3, 287–308. [Google Scholar]
 Deb K. (2009) Multiobjective optimisation using evolutionary algorithms (reprinted version), John Wiley & Sons, Chichester, England, UK. [Google Scholar]
 Moore J., Chapman R. (1999) Application of particle swarm to multiobjective optimization, Department of Computer Science and Software Engineering, Auburn University. [Google Scholar]
 Coello Coello C.A., Lechuga M.S. (2002) MOPSO: A proposal for multiple objective particle swarm optimization, in: Proc. Congr. Evolutionary Computation (CEC’2002), May, Honolulu, HI vol. 1, pp. 1051–1056. [Google Scholar]
 Xue B., Zhang M., Browne W.N. (2013) Particle swarm optimization for feature selection in classification: A multiobjective approach, IEEE Trans. Cybern. 43, 6, 1656–1671. [PubMed] [Google Scholar]
 Zheng Y.J., Ling H.F., Xue J.Y. (2014) Population classification in fire evacuation: A multiobjective particle swarm optimization approach, IEEE Trans. Evol. Comput. 18, 1, 70–81. [CrossRef] [Google Scholar]
 Zain M.Z.B.M., Kanesan J., Chuah J.H., Dhanapal S., Kendall G. (2018) A multiobjective particle swarm optimization algorithm based on dynamic boundary search for constrained optimization, Appl. Soft Comput. 70, 680–700. [CrossRef] [Google Scholar]
 Shirangi M.G., Mukerji T. (2012) Retrospective optimization of well controls under uncertainty using kernel clustering, Monterey, California, USA. [Google Scholar]
 Wang H., Echeverria Ciaurri D., Durlofsky L.J., Cominelli A. (2012) Optimal well placement under uncertainty using a retrospective optimization framework, SPE J. 17, 1, 112–121. [CrossRef] [Google Scholar]
 Shirangi M.G., Durlofsky L.J. (2016) A general method to select representative models for decision making and optimization under uncertainty, Comput. Geosci. 96, 109–123. [Google Scholar]
 Celebi M.E., Kingravi H.A., Vela P.A. (2013) A comparative study of efficient initialization methods for the kmeans clustering algorithm, Expert Syst. Appl. 40, 200–210. [CrossRef] [Google Scholar]
 Adeniran A., Adebayo A., Salami H., Yahaya M., Abdulraheem A. (2019) A competitive ensemble model for permeability prediction in heterogeneous oil and gas reservoirs, Appl. Comput. Geosci. 1, 100004. https://doi.org/10.1016/j.acags.2019.100004. [CrossRef] [Google Scholar]
 Meira L.A., Coelho G.P., Santos A.A.S., Schiozer D.J. (2015) Selection of representative models for decision analysis under uncertainty, Comput. Geosci. 88, 67–82. [CrossRef] [Google Scholar]
 Rahim S., Li Z., Trivedi J. (2015) Reservoir geological uncertainty reduction: An optimizationbased method using multiple static measures, Math. Geosci. 47, 4, 373–396. [Google Scholar]
 Anitha P., Patil M.M. (2019) RFM model for customer purchase behavior using Kmeans algorithm, J. King Saud Univ. – Comput. Inf. Sci. https://doi.org/10.1016/j.jksuci.2019.12.011 [Google Scholar]
 Jansen J.D., Fonseca R.M., Kahrobaei S., Siraj M.M., Van Essen G.M., Van den Hof P.M.J. (2014) The Egg model – a geological ensemble for reservoir simulation, Geosci. Data J. 1, 192–195. [CrossRef] [Google Scholar]
 Fonseca R., Reynolds A.C., Jansen J.D. (2016) Generation of a Pareto front for a biobjective water flooding optimization problem using approximate ensemble gradients, J. Petrol. Sci. Eng. 147, 249–260. [CrossRef] [Google Scholar]
 Soares J., Vale Z., Borges N., Lezama F., Kagan N. (2017) Multiobjective robust optimization to solve energy scheduling in buildings under uncertainty, in: International Conference on Intelligent System Application to Power Systems, September 17–21, San Antonio, Texas, USA, IEEE [Google Scholar]
 Criqui P. (2001) POLES: Prospective outlook on longterm energy systems. Information document, LEPIIEPE, Grenoble, France. http://web.upmfgrenoble.fr/lepiiepe/textes/POLES8p01.pdf. [Google Scholar]
 Lapillonne B., Chateau B., Criqui P., Kitous A., Menanteau P., Mima S., Gusbin D., Gilis S., Soria A., Russ P., Szabo L., Suwa W. (2007) World energy technology outlook – 2050 – WETOH2, PostPrint halshs00121063, HAL. [Google Scholar]
 Ljung L. (1999) System identification – theory for the user, PrenticeHall. [Google Scholar]
All Tables
Comparison of the number of simulation calls required for various methods – base case Egg model.
The MOPSO optimizer variables for production optimization under geological uncertainty.
The MOPSO optimizer variables for production optimization under economic uncertainty.
All Figures
Fig. 1 Original implementing approach for MOPSO. 

In the text 
Fig. 2 Schematic of kmeans algorithm, dots are training examples, and triangles are cluster centroids. (a) Original dataset, (b) random initial cluster centroids, (c–f) illustration of running two iterations of kmeans. In each iteration, each training example is assigned to the closest cluster centroid; then each cluster centroid is moved to the mean of the points assigned to it. 

In the text 
Fig. 3 Eggshaped model depicting the position of injectors and producers. 

In the text 
Fig. 4 Paretofronts of various iterations achieved through adjusted MOPSO – Basecase. 

In the text 
Fig. 5 Pareto front, and the search space of particles, after 100 iterations of adjusted MOPSO. 

In the text 
Fig. 6 The obtained Pareto fronts using the ensemble of 100 realizations; the dotted line is the robust Pareto front. 

In the text 
Fig. 7 Ten permeability maps selected by the kmeans approach. 

In the text 
Fig. 8 the obtained Pareto fronts using the ensemble of 10 selected realizations by kmeans method; the dottedline is RO Pareto front. 

In the text 
Fig. 9 Comparison of the CDFs for SNPV (a) and LNPV (b) after 50 iterations for 10 selected realizations with that of 100 realizations for production optimization of the Egg model. 

In the text 
Fig. 10 Comparison of the CDFs for RF of all 100 realizations and their corresponding 10selected realizations resulted from (a) using the maximum allowable injection rate, (b) using injection profile of best LNPV in Figure 7, (c) using injection profile of best SNPV in Figure 6. 

In the text 
Fig. 11 Comparison of the obtained Pareto front after 50 iterations of robust optimization on different set of realizations of the Egg model (the definitions of RF cases: (a) using the maximum allowable injection rate, (b) using injection profile of best LNPV in Fig. 7, (c) using injection profile of best SNPV in Fig. 6). 

In the text 
Fig. 12 Oil price functions for economic uncertainty characterizationdata from Siraj et al. [17]. 

In the text 
Fig. 13 The obtained Pareto fronts using the ensemble of 10 oil price scenarios, mean oil price, and RO. 

In the text 
Fig. 14 NPV changes with simulation time for a hypothesis case of using the injection profiles belong to (a) meanvalue scenario and (b) robust scenario while both are utilizing the meanvalue oil price to generate objective functions through time. 

In the text 
Fig. 15 The obtained Pareto fronts from RO using the ensemble of 10 oil price scenarios on 10 selected realizations. 

In the text 
Fig. 16 Optimized injection scenarios [m^{3}/day] in Egg model obtained by MOPSO bluecolored lines: the best LNPVs, and redcolored lines: the best SNPVs. 

In the text 
Algorithm A The adopted MOPSO algorithm improvement utilized for operation optimization in waterflooding problems. 

In the text 