© M.M. Moshir Farahi et al., published by IFP Energies nouvelles, 2021

Nomenclature

B : Boundary in MOPSO algorithm

BF: Boundary factor in MOPSO algorithm

C : Random variable

C_o : Oil’s compressibility

C_r : Rock’s compressibility

C_w : Water’s compressibility

d : Discount factor

H : Problem constraint

J_RO : Robust objective function

k : Number of means in k-means method

K_ro0 : Relative-permeability’s endpoint to oil

K_rw0 : Relative-permeability’s endpoint to water

M : Number of problem constraints

$m_i^{\left(t\right)}$ : Means in cluster t

NDS-W: Worst Non-Dominated Solution

NDS-E: Efficient Non-Dominated Solution

N_p : Population’s number

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_c : Capillary pressure

Phi: Porosity

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 time-step

Q_wp,k : Rate of water rate at time-step k-production

Q_wi,k: Rate of water rate at time-step k-injection

r_o : Oil price function

SaF: Saturation factor in MOPSO algorithm

SF: Shrink Factor in MOPSO algorithm

s ^(t) : Clusters in k-means method

S_wir : Irreducible water-saturation

S_wi : Initial water-saturation

S_or : Residual oil-saturation

t_k : total time in k’s time-step

V: Unit vector

X : Control variable

X_i : Initial population

x_p : Data points in k-means method

α_i : Randomly selected coefficients in ARMA model

α₀ : Base oil price

β_i : Randomly selected coefficients in ARMA model

Δt_k: Difference of successive time-steps

τ_ref : Discounting referenced time

ω_k : White-noise sequence

1 Introduction

One of the main concerns of the decision-making process in reservoir management is high-level uncertainties and debatable parameters. This effect is a crucial challenge, especially in model-based 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, Multi-Objective Optimization (MOO) is a practical tool for decision-making that should be considered while identifying uncertainties.

Water-flooding optimization has been investigated in many kinds of research in petroleum literature, which focused mainly on single-objective problems [1–3]. Accordingly, the field life-cycle’s economic performance has been reflected in the form of Long-term Net Present Value (LNPV) [4–6]. Although, neglecting the short-term goals in an intended optimal design generally leads to a low profit in the short-term compared to that of a reactive strategy. A limited constraint, e.g., maximum allowed water-cut, is the only controller for designing the water injection scenario in a reactive strategy, selecting which conflicts with LNPV maximization in a reservoir. However, short-term production gains are included in very recent studies [7–11]. One can generally use an annual discount rate in NPV calculation to consider the short-term financial targets. Accordingly, the simultaneous optimization of LNPV and Short-term NPV (SNPV) via an appropriate MOO-handling algorithm to form a Pareto front provides a valuable mechanism that helps the decision-maker trade-off between different goals [10, 12].

Model-based 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 short-term 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 multi-objective production optimization algorithms were used in previous studies. Some of them utilized gradient-based methods [21]. Some others used the derivative-free techniques for production optimization [22]. These approaches are distinct by three selection methods including; criterion-based, aggregation-based, and Pareto-based. In a standard design to handle multi-objective problems, a single-objective 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 non-convex results relying on a Pareto front [24]. Regarding these challenges, the selection algorithms accompanying Pareto have been introduced as more efficient approaches. Population-based evolution algorithms can employ the Pareto-ranking 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., history-matching problems. Reference [26] has utilized MOPSO to the reservoir history-matching 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 model-based multi-objective development problem was solved by performing three multi-objective methods of adjoint (weighted-sum method and gradient-based method), Genetic Algorithm, and PSO. They found that MOPSO is more appropriate for the applied medium-scale 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 decision-makers to manage the risks and make reasonable decisions. As a representative work in case of oil price-related 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 short-term and long-term 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 long-term and short-term 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 long-term and short-term profit simultaneously while considering various reservoir-related 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 non-dominated 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 multi-objective optimizations.

LNPV as the first objective function represents the long-term gains. Since economic uncertainty risk grows with time, a highly discounted NPV (25%) is considered as short-term gains, SNPV. The proposed approach has been practiced on an egg-shaped standard model. To characterize the geological uncertainty space, ensemble members of 100 model realizations and 10 selected ones through the k-means 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 Pareto-based MOPSO approach has been applied to maximize the short-term and long-term 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 oil-price 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 oil-price 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 price-related uncertainties.

2 Optimization problem for reservoir management

Field existing criteria and production constraints identify the field’s short-term objectives that improve decision-making in the reservoir life-cycle. 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 short-term objectives yields a small reduction to the field’s life-cycle target [7, 19]. To start the analysis, a common statement of the NPV as a field intended objective function, which is under water-flooding process, provided by [29], is called:

$$\mathrm{NPV}=\sum _{k=1}^T\mathrm{}\left(\frac{\left\{\left[\left(Q_{\mathrm{o},k}\right)·p{r}_{\mathrm{o}}-\left(Q_{\mathrm{wp},k}\right)·p{r}_{\mathrm{wp}}\right]-\left[\left(Q_{\mathrm{wi},k}\right)·p{r}_{\mathrm{wi}}\right]\right\}·\Delta t_k}{{\left(1+d\right)}^{\frac{t_k}{{\tau }_{\mathrm{ref}}}}}\right)$$(1)

in which, the prices and rates are respectively represented by pr [in USD/m³] and Q [in m³/day]. d [1/year], is a parameter representing a discounting factor. Δt_k and t_k [day] are difference of time-steps and the cumulative time after time-step 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 long-term and short-term 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 short-term objective function, i.e., SNPV, concurrently with maximizing LNPV. The assumed discount rate corresponding with SNPV significantly considers the production in the field’s short-target points of view compared with the project’s future, having no economic meaning. However, this stimulates the short-term achievement of implementing the production scenario [5]. Considering the above-mentioned cases of short- and long-term targets, the following two objectives are defined:

• The long-term achievement: LNPV (d = 0 in Eq. (1)).

• The short-term achievement: SNPV (d = 0.25 in Eq. (1)).

To handle the optimality concept simultaneously for LNPV and SNPV in the multi-criterion optimization problem of water-flooding 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 trade-off among these objectives since growing one purpose makes another deteriorate. The original form of MOO is employed to be utilized in the water-flooding problem, as:

$$\begin{array}{cc}\mathrm{Optimize}:& J\left(X\right)=\left[J_1\right(X),J_2(X),\dots ,J_f(X\left)\right]\\ \mathrm{subject} \mathrm{to}:& \\ & \begin{array}{c}{X}_{\min }\le X_i\le X_{\max };i=\mathrm{1,2},\dots ,n_{\mathrm{p}}\\ h_{\mathrm{ii}}\left(X\right)\le 0 \mathrm{or} h_{\mathrm{ii}}\left(X\right)=0;\mathrm{ii}=\mathrm{1,2},...,y.\end{array}\end{array}$$(2)

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 time-dependent problem limitations. These constraints should be defined by case. To obtain a set of results that trade-off among the contradicting purposes, the Pareto-dominance 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₁ and X₂, X₂ < X₁, where if and only if: 1) X₂ is not better than X₁ in all J_M(X), for each M∈[1, …, f], and 2) X₁ is notably higher than X₂ in J_M(X), at the minimum one M∈[1, …, f]. The set of solutions is obtained as a Pareto-optimized when it is non-dominated, respecting the achievable range. Corresponding objective function sets obtain the Pareto-optimized 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 market-related 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 multi-criterion 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 Pareto-dominance concept. A mutation operator is utilized that benefits the search scheme. Besides that, to guide future generations’ search, the non-dominated outputs will be saved in a secondary population reservation. Based on the Pareto-ranking 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].

In this algorithm, the best positions that are met by particles are stored; thus, all the previously obtained non-dominated 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 Pareto-dominance, 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 non-dominated 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 water-flooded 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 k-medoids [45], implementing a starting control plan to rank the obtained NPVs [32], employing some static parameters based on simulation to perform k-means clustering [46], and a combination strategy of permeability-based and flow-based quantities [47] are introduced previously in the literature. We select to use the k-means 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 k-means clustering is based on a prototype distribution clustering idea. In this technique, one can find the user-specified 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 $m_1^{\left(1\right)}\dots m_k^{\left(1\right)}$ means, k-means 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:

$$s_i^{\left(t\right)}=\left\{x_p:{\Vert x_p-m_i^{\left(t\right)}\Vert }^2\le {\Vert x_p-m_j^{\left(t\right)}\Vert }^2\forall j,1\le j\le k\right\},$$(3)

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.,:

$$m_i^{(t+1)}=\frac{1}{\left|S_i^{\left(t\right)}\right|}\sum _{x_j\in S_i^{\left(t\right)}}\mathrm{}{x}_j.$$(4)

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 k-means method, from which fast computation, excellent performance on large datasets, and needing only one input parameter (k) are of great importance [52].

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]:

$$J_{\mathrm{RO}}=\frac{1}{N_e}\sum _{i=1}^{N_e}\mathrm{}{J}_i,$$(5)

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 water-flooding 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 life-cycle 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 short-term purposes simultaneously with meeting the long-term 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 short-term economic goals. Although NPV with a 25% discount rate may have no meaning in the economic sense, it weights the short-term 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 egg-shaped and three-dimensional. It composes 25 200 blocks with 18 533 ones being active that indeed build that egg-shaped 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.

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³/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 time-steps 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 egg-shaped model, as mentioned in [53]. Also, the values of oil and operational prices used in equation (1), are listed in Table 1.

5.3 Model-based 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.

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 single-objective 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 concave-down shape is desired, representing a complete Pareto front.

The obtained results are consistent with previous works by [2, 7], showing maximization of the short-term 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 decision-maker, the developed Pareto front gives an effective scheme that makes it possible to trade-off 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.

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 multi-objective 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:

$$\mathrm{NP}{\mathrm{V}}_{\mathrm{RO}}=\frac{1}{N_{\mathrm{geo}}}\sum _{i=1}^{N_{\mathrm{geo}}}\mathrm{}\mathrm{NP}{\mathrm{V}}^i,$$(6)

where N_geo is the number of geological realizations and NPVⁱ 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.

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 decision-maker 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 decision-making 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 life-cycle target.

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 Net-to-Gross ratio, are considered to be known and constant. In order to use the k-means 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 k-means algorithm. Then, dividing the realizations into 10 groups (clusters) by the algorithm is performed by minimizing the sum of square distances belonging to within-cluster.

In Figure 7, 10 field permeability realizations, selected through performing the k-means method, are presented.

The Pareto fronts resulted from implementing the optimization approach on the 10 selected models by the k-means method are represented in Figure 8. The same parameters of the algorithm as mentioned in Table 4 are applied for practicing k-means approach. The robust Pareto front of these models is also depicted, showing an acceptable consistency with the selected realizations.

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 short-term 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.

As it is stated before, reservoir heterogeneity plays an important role in water-flooding 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³/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 10-selected 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.

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 k-means 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 water-flooding performance in this case study. As shown in this figure, regarding the computational time that may affect the optimality of problem-solving, one can replace the whole realizations-related 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.

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’ life-cycle 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 price-related 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 water-injection/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 Long-term Energy Systems (POLES) [56–57]. Siraj et al. utilized a simplified Auto-Regressive-Moving-Average model (ARMA) [58] to make oil price time-series 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:

$$r_{o_k}={\alpha }_0+\sum _{i=1}^6\mathrm{}{\alpha }_i{r}_{o_{k-1}}+\sum _{i=1}^6\mathrm{}{\beta }_i{\omega }_{k-1},$$(7)

where ω_k is a white-noise sequence, α_i and β_i are the randomly selected coefficients [8]. In an approach by [17], 10 scenarios were generated using 471 USD/m³ 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:

$$\mathrm{NP}{\mathrm{V}}_{\mathrm{RO}}=\frac{1}{N_{\mathrm{eco}}}\sum _{i=1}^{N_{\mathrm{eco}}}\mathrm{}\mathrm{NP}{\mathrm{V}}^i.$$(8)

The applied parameters of the algorithm are mentioned in Table 5.

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 mean-value 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 mean-value. Any point on the robust/mean Pareto front gives the injection procedure corresponding to the selected LNPV-SNPV pairs, which can be considered a low-risk scenario related to economic uncertainties.

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 mean-value 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 re-optimize the injection scenario. In that case, continuing the reservoir production with the current condition and performed optimization is cost-effective only by increasing the oil price prospective. At the same time, the best long-term scenario is a conservative one that makes the production plan cost-effective for a minimum of seven years. The increasing oil price perspective enhances the analysis control. These analyses show the main contributions in the decision-making process, which include using the oil price mean-value scenario, assigning the optimum scenario based on economic, technical, and operational assumptions, and updating it for any required time.

In Figure 14b, the previous diagram is re-produced for the injection profile resulted from RO. To generate this figure, the robust scenario’s injection profile and the mean-value 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 long-term 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 k-means 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:

$$\mathrm{NP}{\mathrm{V}}_{\mathrm{RO}}=\frac{1}{N_{\mathrm{eco}}}\sum _{i=1}^{N_{\mathrm{eco}}}\mathrm{}\left(\frac{1}{N_{\mathrm{geo}}}\sum _{i=1}^{N_{\mathrm{geo}}}\mathrm{}\mathrm{NP}{\mathrm{V}}^i\right).$$(9)

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 k-means 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 LNPV-SNPV pairs, which can be considered a low-risk scenario related to both geological and economic uncertainties.

Figure 16 depicts the proposed injection scenarios for the two end-points of the obtained Pareto fronts in various analyses performed on the case study. This may provide useful information for the decision-making 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 long-time increasing in water injection rates. As another example, assume one decides to gain a maximum LNPV according to the Pareto front of either mean-value optimization or RO that are approximately the same as shown in Figure 13. Thus, the provided analyses propose two different scenarios for the decision-maker 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.

6 Conclusion

Geological and economic uncertainties in model-based optimization are the main difficulties in maximizing short-term financial goals. In this study, using a modified multi-objective approach in the optimization process, namely MOPSO, with explicit consideration of such uncertainties (i.e., RO) resulted in the robust prediction of long-term economic goals while allowing the right balance with short-term 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 egg-shaped 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 base-case, 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 k-means 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 decision-maker 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 gradient-based 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 water-flooded 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.

To start evaluating process in the initialization step, a random population (X_i), which includes controlling parameters within its boundary (B), is chosen:

$$X_i=X_{\min }+C1.B_i,i=1,...,n,B_i=\mathrm{abs}(X_{\max }-X_{\min })/2,B_i^0=B_i.$$(A1)

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:

$$X_{\min }^{\mathrm{\prime }}=X_i-B_i, X_{\max }^{\mathrm{\prime }}=X_i+B_i.$$(A2)

Then offspring population is constructed inside the new boundary by:

$$X_i^{\mathrm{\prime }}=X_{\min }^{\mathrm{\prime }}+C2·\mathrm{abs}\left(X_{\max }^{\mathrm{\prime }}-X_{\min }^{\mathrm{\prime }}\right).$$(A3)

The next step is to determine the probability of each variable to go to its boundary by:

$$\begin{array}{l}\mathrm{if} X_i^{\mathrm{\prime }}<A.C3.\mathrm{abs}\left(X_{\max }-X_{\min }\right)+X_{\min } \mathrm{then},X_i^{\mathrm{\prime }}=X_{\min },\\ {\mathrm{else} X}_i^{\mathrm{\prime }}>X_{\max }-A·C4·\mathrm{abs}\left(X_{\max }-X_{\min }\right), X_i^{\mathrm{\prime }}=X_{\max }.\end{array}$$(A4)

Now, respecting the new population evaluation of the objective functions (SNPV and LNPV), storage of the non-dominated 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 decision-maker 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:

$$B_i=B_i·\left(\mathrm{SF}\right).$$(A5)

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:

$$\mathrm{BF}=A·C6,$$(A6)

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 $B_i=B_i^0\left(\mathrm{BF}\right)$ if:

$$B_i<\left(\mathrm{BF}\right)·\mathrm{abs}\left(X_{\max }-X_{\min }\right)+X_{\min }.$$(A7)

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 finite-sized repository, leading to quick overfilling. When new non-dominated 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

All Tables

All Figures