Applying hybrid support vector regression and genetic algorithm to water alternating CO2 gas EOR

Water alternating CO2 gas injection (WAG CO2) is one of the most promising enhanced oil recovery techniques. The optimization of this process requires performing many time-consuming simulations. In this paper, an intelligent hybridization based on support vector regression (SVR) and genetic algorithm (GA) is introduced for the WAG process optimization in the presence of time-dependent constraints. Multiple SVRs are used as dynamic proxy to mimic numerical simulator behavior in real time. Latin hypercube design (LHD) is applied to generate the proper runs to train the proxy and ten supplementary runs are randomly chosen to validate it. The goal of GA in this study is twofold. First, it is employed during the training of multiple SVRs to find their appropriate hyper-parameters. Second, once the training and validation of the dynamic proxy are done, the GA is coupled with it to find the optimum WAG parameters which maximize field oil production total (FOPT) subject to time-dependent water-cut constraint and some domain constraints. The task is formulated as a non-linear constrained optimization problem. A semi-synthetic WAG CO2 case is used to examine the reliability of the approach. The results show that the established dynamic proxy is fast and accurate in reproducing the simulator outputs. The hybridization proxy-GA is demonstrated to be reliable for the real-time optimization of the formulated WAG process. C © 2020 The Authors. Greenhouse Gases: Science and Technology published by Society of Chemical Industry and John Wiley & Sons, Ltd.


Introduction
F ossil fuels will remain the world's primary energy source over the next few decades. 1,2 Many oil reservoirs are beyond or at the end of their primary stage of recovery, leaving considerable oil Correspondence to: Ashkan Jahanbani Ghahfarokhi, Department of Geoscience and Petroleum, Norwegian University of Science and Technology, S. P. Andersens veg 15b, 7031 Trondheim, Norway. E-mail: ashkan.jahanbani@ntnu.no MN Amar et al. Original Research Article: Applying hybrid support vector regression and genetic algorithm to water alternating CO 2 EOR Water alternating gas (WAG) is increasingly being applied as an EOR technique in numerous oil fields in the world. 5 This technique is recognized as the process of injecting water and gas in many cycles. The main advantages of this process are the microscopic and macroscopic sweep efficiencies enhancement. 5 The microscopic sweep is enhanced by gas injection, which is effective in displacing the oil not reached by water, while the macroscopic sweep is enhanced by subsequent water slugs, which control the mobility. Around 80% of the implemented WAG projects in USA were fruitful. 6 In addition, Skauge et al. 7 by reviewing 59 WAG field applications reported that the average oil recovery increased by up to 10% of the original oil in place (OOIP). Among different injection gases in WAG processes, impure/pure CO 2 is the most effective choice, because of the oil recovery factor improvement and the possibility of mitigating CO 2 emissions by storing CO 2 in the reservoir. 8,9 Several papers on WAG and WAG CO 2 processes revealed noticeable increase in the recovery factor using this EOR technique. 10,11 Various key parameters have considerable impacts on the successful design of any WAG process. 5 These parameters include the initialization time, water and gas injection rates, the half-cycle injection time, the WAG ratio and slug size. Proper determination of these parameters is a complex task that needs many computationally expensive numerical simulations. Therefore, developing a proxy model seems to be an alternative to deal with the computational issues, as demonstrated in our previous works 12,13 . A proxy model consists of using a simple mathematical paradigm that replaces the simulator and generates outputs quickly without sacrificing the accuracy. 14,15 This approach is recognized by different terminologies in the literature such as surrogate and meta modeling. 16 Proxy models are proved to be effective in various reservoir simulation tasks related to forecasting, management and optimization purposes. Jalali et al. 17 proposed a surrogate model to perform uncertainty analysis on a coal bed methane reservoir. In their study, the authors applied artificial neural networks (ANN) to build the surrogate model after identifying the most impacting input parameters using key performance indicators (KPIs). Then, Monte Carlo simulation was performed on the surrogate model in order to assess the uncertainties associated with reservoir parameters. The results of their study showed that the implemented approach led to a fast and accurate estimation of the needed outputs. Yao et al. 18 established a proxy for modeling the production profiles in steam assisted gravity drainage (SAGD) processes. The proxy models in their study were based on techniques such as Box-Jenkins models, linear non-recursive models and recursive ARX (autoregressive with exogenous input) models. Their findings demonstrated the reliability of the proposed proxy models for alleviating the formulated problem. Amini et al. 19 applied a surrogate reservoir model (SRM) to a real case CO 2 sequestration project to quantify the operational and geological uncertainties involved within reasonable time. The proposed workflow allowed the prediction of the distribution of CO 2 and pressure throughout the reservoir with satisfactory accuracy and a short runtime. Mohammadi and Ameli 20 optimized the control parameters of a SAGD process, namely injection rate and pressure, production, injection and offset well heights, offset well pressure and its injection and production periods, using response surface model (RSM). Their RSM-based substitutional scheme exhibited reliable performance in terms of time and accuracy. You et al. 21 proposed a robust computational proxy for optimizing carbon dioxide-enhanced oil recovery (CO 2 -EOR) project under multi-criteria goals. The task consisted of finding the scenarios that maximize the CO 2 storage, oil recovery and the project's economic outcomes. To this end, the authors proposed combination of a surrogate model, which was built using ANN, with multi-objective optimizers. Their results showed that the provided computational framework was very robust and can be considered as a user-friendly decision-making tool for similar optimization problems. Nait Amar et al. 12 and Nait Amar and Zeraibi 13 combined back-propagation techniques for learning ANN models and the best-found model was employed as proxy for optimizing mono-and multi-objective WAG problems, respectively. The results showed noticeable accuracy and reasonable runtime. Redouane et al. 16,22 developed efficient adaptive workflow for optimizing well placement in real fractured and synthetic reservoirs. They claimed that the proposed steps led to an optimized number of runs to build a representative proxy for well placement problems.
Over the last two decades, there has been a surge of interest in the introduction of soft computing techniques to resolve practical problems in several fields of science and technology, including petroleum engineering. 23 29 It is characterized by the high generalization capability which results from the well-formulated mathematical learning concept. 30 It has been applied in various fields such as signal processing, 31 finance, 32 biology, 33 biomedicine, 34 and engineering. 35 Many articles shed light on the successful application of this tool in petroleum and reservoir engineering. Bian et al. 36 applied SVR coupled with grey wolf optimization to model wax disappearance temperature (WDT). The established paradigm can estimate WDT with an average absolute relative deviation (AARD) of 0.7128%. Esfahani et al. 37 established an accurate model for estimating natural gas density using least square support vector machine (LSSVM). The performance analyses revealed that their model outperformed the prior approaches. Ziaee et al. 38 implemented an SVR model for predicting the solubility of CO 2 in different polymers. Their findings demonstrated the reliability of the proposed model. Nait Amar and Zeraibi 39 modeled minimum miscibility pressure (MMP) of CO 2 -oil systems using SVR optimized by artificial bee colony (ABC). It was found that the SVR-ABC model can predict the MMP of the systems with high accuracy. One of the studies that addressed the application of SVR as proxy of numerical simulators comes from the work of Ahmadi et al., 1 where an SVR model was developed as a static proxy to investigate the efficiency of CO 2 injection. No previous study has employed SVR technique as time-dependent proxy in the investigation of dynamic processes.
Genetic algorithm (GA), which is a consolidated approach in evolutionary computation, has received growing attention in resolving optimization problems. 40 Well test characterization 41 , history matching 42,43 , and well placement optimization 16,44 are among the well-known reservoir engineering tasks where GA was successfully applied.
This article focuses on the development, evaluation and validation of a fast and robust approach to optimize WAG CO 2 process subjected to time-dependent constraints. This approach consists of hybridization of multiple SVRs with GA. Multiple SVRs are constructed to generate the outputs of the numerical simulator in real-time with high accuracy. Latin hypercube design (LHD) is the space filling method applied to specify the runs utilized in the building phase of the proxy. GA is used initially during the training of the multiple SVRs to find the optimum SVR hyper-parameters. After checking the reliability of the proxy, it is coupled with GA to investigate the optimum WAG parameters, namely the initialization time, the injection rates of water and gas, the half-cycle injection time of water and gas, WAG ratio and slug size, which maximize field oil production total (FOPT) with respect to water-cut time-dependent constraint and domain constraints. This paper differs from published papers and literature as: (1) it sheds light on the development of time-dependent SVR as a dynamic proxy for optimizing WAG CO 2 process; (2) the developed SVR proxy generates the parameters involved in the optimization in real-time; (3) the parameters included in the optimization are not limited to the ones frequently considered in WAG process (the injection rates of water and gas, half-cycle injection time, the WAG ratio and the slug size), but another vital factor is considered, namely the initialization time of the process.
In this paper the second section introduces the theoretical aspects of SVR and GA. The third section describes the reservoir model utilized in this study. Fourth and fifth sections illustrate the mathematical formulation of the studied WAG CO 2 process and the procedure of building the SVR-based proxy, respectively. Results are presented and discussed in the sixth section. The last section sheds light on the main findings of the work.

Theory Support vector regression (SVR)
SVR is an advanced machine learning methodology which was introduced by Vapnik 45 . SVR aims to approximate a function that computes the functional dependency between targets T = {t 1 , t 2 , . . . , t m } defined on R and inputs X = {x 1 , x 2 , . . . , x m }, where x i ∈ R n and m is the database size. This function can be expressed as: where w and b are the weight vector and bias term, respectively and ϕ(x) is a high dimensional mapping feature.
To obtain w and b, the following regularized risk function should be minimized: Original Research Article: Applying hybrid support vector regression and genetic algorithm to water alternating CO 2 EOR 1 2 w 2 represent the empirical error and the flatness degree of the function, respectively. The regularization parameter C > 0 is recognized as a penalty parameter, which shows the trade-off between the complexity of the model and the empirical error.
ε-insensitive loss function L( f (x i ) − t i ) is applied for computing the empirical error. This function is formulated as follows 46 : where ε is the error tolerance.
The optimum parameters are then obtained by formulating the constrained optimization problem as: i and ξ + i are non-negative slack terms. By introducing Lagrange multipliers, the constrained optimization function of (Eqn. 4) can be transformed into dual space, and the solution is 47 : where α i and α * i denote Lagrange multipliers and must satisfy the constraints 0 ࣘ α i , α * i ࣘC, and the term K(x i , x j ) is called the Kernel function. The latter aims to map the input space onto some higher dimensional space (feature space), allowing SVR more flexibility in modeling complex cases. Gaussian type, polynomial function and radial basis function (RBF) are among the well-known Kernel functions. 48 In this study, RBF is considered as the Kernel function and is expressed as follows: where γ is the Kernel parameter. The reliability of SVR paradigm depends on the proper determination of C, ε and the parameter γ of the Kernel function. Hence, it is necessary to optimize these parameters using robust algorithms to automatically achieve the proper combinations.

Genetic algorithm and constraints handling strategy
Genetic algorithm GA is recognized as the first nature-inspired algorithm which applies the so-called genetic operators for exploring and exploiting the search space. 40 As a population-based optimization technique, GA randomly generates an initial population of possible solutions and represents them in the form of chromosomes. A fitness function is considered as an assessment criterion to differentiate the quality of the chromosomes. Across the iterations, the chromosomes are subjected to the genetic operators including selection, elitism, crossover and mutation. Selection is performed to pick the individuals to be parents and generate new offspring. Elitism allows the survival of the fittest elements and their passage to the next generation. Crossover consists of exchanging parts among chromosomes to give new ones. Mutation allows the local random research aspect for GA by varying certain genes of chromosome. These operators are reiterated until a stopping condition is satisfied.
In the present work, GA is used to optimize the SVR hyper parameters and the constrained WAG CO 2 process (after building the proxy). These two applications are explained in the following sections.

Constraint handling
This study utilizes 'three feasibility rules' method 49   To implement these rules in GA, a linear ranking fitness function is considered in this work. It should be noted that the selection operator of GA is related to this fitness function.
Original Research Article: Applying hybrid support vector regression and genetic algorithm to water alternating CO 2 EOR MN Amar et al.

Latin hypercube design (LHD)
LHD is one of the most widely used sampling designs in the field of design of experiment, mainly intended for numerical and computational purposes 48,50 . It was introduced by McKay et al. 51 in 1979. LHD is characterized by simplicity of construction, good coverage of the design space and uniform distribution of the points generated on the factorial axes. These advantages are due to the principle of LHD, which is basically the distribution made by dividing the intervals of the design variables into a number of bins with the same selection probability, and for each variable projection in a bin, there will only be one sample to select.
For m samples to be generated and n variables, the number of LHD combinations that can be obtained is given by the following expression: Based on this formula, it can be deduced that the number of potential LHDs is very high. However, to select the LHD that can cover the design space with the best repartition, several criteria and methodologies were proposed in the literature. More details can be found in published papers 48,52 . In the present work, a combination of Morris and Mitchell criterion 53 and GA were implemented for generating the best LHD.

Characterization of the reservoir model
In the present work, real dynamic and Pressure-Volume-Temperature (PVT) data from one of the Algerian fields are implemented with a synthetic static model to build the numerical model. Compositions of the injected gas and the reservoir fluid are reported in Table 1. The PVT properties and fluid viscosity were modeled by applying the three-parameter Peng Robinson (PR) equation of state (EOS) and the Lorenz-Bray-Clark (LBC) correlation, respectively. A tuning procedure was applied to match the experimental PVT data. Comparison of the experimental data and the predicted values using the PR-EOS after the tuning are illustrated in Table 2 (bubble point pressure) and Fig. 1 (the volumetric  parameters). According to Table 2 and Fig. 1, very satisfactory agreements were achieved between the measurements and the predictions of the tuned PR-EOS. It is necessary to mention that during the tuning steps, a detailed composition (up to C36+) was     Figure 3 illustrates the positions of these wells. Eclipse 300 is the numerical simulator used in the development. It is worth highlighting that the studied reservoir contains 16.5 × 10 6 sm 3 of oil in place. The start time is the beginning of 2020, and the end of the simulations is fixed to the time needed to inject 1.2 PV of gas.

WAG CO 2 optimization problem formulation
As achieving a high FOPT value is a paramount goal in any WAG injection process, this parameter is specified as the objective function to maximize, as a function of the design parameters of WAG CO 2 process. The design parameters and the possible ranges are summarized as follows:    In practice, field water-cut (FWCT) is often restricted by an upper limit. Therefore, a maximum water-cut value of 92% is considered as the problem constraint.
Based on real cases that are implemented in many fields 10,11 , we have considered the same half cycle time and initialization time for the injection wells, while for the injection rates, we split the FWIR and FGIR equally among the four injection wells.

Proxy development
To build a powerful dynamic proxy (multiple SVRs) that can emulate the parameters involved as functions of time, and capture the complexity of the WAG CO 2 process, a proper handling of the outputs and inputs of insignificant number of runs is done. According to the WAG CO 2 optimization problem formulated in the fourth section, two main outputs are basically required (FWPR and FOPR) to obtain the others. It is clear that FOPT is deduced from FOPR and time, while FWCT depends on FWPR and FOPR (Eqn 14). Therefore, multiple SVRs based proxy aims to estimate FWPR and FOPR as functions of time, and the other parameters (FOPT and FWCT) are deduced accordingly.
Before proceeding to the development step, a database is needed. For this purpose, appropriate number of numerical simulator runs are selected. From the IT and HCT constraints (Eqns 9 and 12), 3 and 4 levels are attributed to these parameters, respectively. Concerning FWIR and FGIR factors, 12 levels are considered for each. It should be noted that parameter levels definition consists of dividing the interval of a design variable into equal size sub-intervals and inserting points in the middle. Twenty-four runs are attributed to each IT level, and this results in a total number of 72 runs (24 × 3). It is worth noting that the runs associated with each IT level cover the whole HCT levels with an equal number of 6 runs (6 × 4). Figure 4 summarizes the splitting procedure. LHD is applied to design these runs. Three supplementary runs which design the interactions between the limits of the variables are added. Accordingly, the total number of runs is 75. Based on the variables of the entire runs, the WAG ratio (defined as the ratio of the injected volume of water to gas) ranges from 0.47 to 3.9, while the slug size (the injected gas volume during HCT gas ) is between 0.028 and 0.29 pore volume injected.
For the generated runs, a broad database is created by extracting the corresponding inputs and outputs for each time-step. The gained datapoints are divided into two sets: one contains the points conforming to the time-steps with gas as the injected fluid; and the other set involves the points where water was the injected fluid. Then, for each implied outcome, two SVRs are built using the proper part: one englobes the periods with gas injection and the other covers the periods with water injection. The global SVR of a given parameter comprises the two SVRs. The parameters required in the formulated optimization task and the considered inputs to develop the multiple SVRs are illustrated in Table 5. In order to include the effect of HCT on the outputs, its value at the prior timestep (t − 1), is employed as an input.
To validate the multiple SVRs proxy and test its robustness with blind data (data not used for the training), ten additional scenarios of WAG CO 2 are chosen randomly.   As mentioned earlier, the accuracy of SVR is related to the proper combination of C, ε and the Kernel function parameter γ . Therefore, GA is implemented to optimize these parameters for each SVR (each SVR aiming at emulating the outputs needed in the WAG CO 2 problem optimization, i.e., SVR-FOPR and SVR-FWPR during gas and water injections). Figure 5 briefly describes the SVR-GA model. The employed GA parameters in the optimization of SVR hyper-parameters are shown in Table 6.

Evaluation of the SVR based proxy model
As mentioned in the fifth section, GA is coupled with SVR during the training phase. This first hybridization aims at finding the proper SVR hyper-parameters. To this end, the steps stated in the workflow of Fig. 5 were followed. Mean square error (MSE) was considered as the fitness function during the training phase of the proxy. MSE is expressed as follows: where O refers to the needed parameter, N represents the number of points and subscripts e and p refer to eclipse and proxy, respectively. The values of the multiple SVRs optimum hyper-parameters are shown in Table 7. These values exhibit the final values of the hyper-parameters gained by GA. In all the cases developed, regularization factor (C) tends to give the best performance towards its largest values, while Epsilon and Kernel parameters give better models with small and medium values, respectively. Furthermore, statistical and graphical error analyses were applied to assess the reliability of the established proxy and evaluate its accuracy in emulating the needed parameters as functions of time. Absolute relative prediction error (ARPE) of each SVR is calculated by applying the following equation: where O e is the eclipse output and O p is the predicted response from the proxy model. In addition, for an in-depth visualization of the results and the robustness of the proxy, different graphical analyses are used: (1) cross plots that compare the alignment of the proxy results versus ideal paradigm, shown by a unit-slope line, (2) a comparative data index plot which consists of representing the simulator outputs against those of the proxy as a function of data index and (3) the plots showing the variations of different needed parameters as functions of time.
The average, minimum and maximum absolute relative errors of the multiple SVRs during the training phase are shown in Table 8. It is worth noting that the values reported in this table are assessed in terms of runs (not single points). Figure 6 displays cross plots for rates (oil and water) and water-cut, while the reported cross plots of Fig. 7 are for the cumulative parameters (FWPT and FOPT). It can obviously be seen from these figures that very satisfactory distributions of the predicted results by the established proxy are achieved (nearby the unit-slope line) for all the parameters. In addition, for an in-depth assessment of proxy approximations and the deviation from the numerical simulator, the detailed statistical results through the ARPE (Table 8) reveal that very small APRE are obtained for all the outcomes. The ARPE values are 1.13 and 0.07% for FOPR and FOPT t_final (at the final time-step), respectively; 1.46 and 0.14% for FWPR and FWPT t_final (at the final time-step), respectively; and 0.27% for FWCT.
As mentioned previously, in order to demonstrate the robustness of the multiple SVRs proxy and to check its accuracy with blind data, ten realizations (not used in the training process) are designed. Table 9 illustrates the performance results of the test step (in terms of runs). Furthermore, Figs 8 and 9 depict cross plots for the needed parameters. According to these cross plots, the accuracy of the developed proxy is justified using this visual survey on blind data. In the same context, and as listed in Table 7 and the visual comparison reported in the cross plots of Figs 8 and 9, it can be stated that the multiple SVRs proxy model implemented in this study provides promising results with the data not used for training the proxy, where very low APRE are reported. This shows the         Figure 13. Optimization results using hybridization SVR-GA. generalization capability of the developed multiple SVRs proxy model. Figure 10 shows the comparison of FOPT in sample series plots for the training and test runs at the final time-step (equivalent to the injection of 1.2 PV of gas). As seen in Fig. 10, good agreement between the multiple SVRs proxy predicted values and the numerical simulator results is noticed.
To evaluate the robustness of the implemented proxy in generating the parameters in real-time, Fig. 11 compares the multiple SVRs proxy with the reservoir simulator (Eclipse 300), in terms of results of randomly selected simulation runs that are included in the training phase. In Fig. 11(a) and (b), the comparisons are made as functions of time for FOPR and FOPT, and for FWPR and FWPT, respectively. Figure 11(c) illustrates the comparison of WCT as a function of time. It can be seen in these figures that the multiple SVRs proxy can reproduce the results of the numerical simulator with high accuracy. Figure 12 reports the comparison of the multiple SVRs proxy results with blind runs of the numerical simulator. The comparisons of FOPR and FOPT as functions of time are depicted in Fig. 12(a), while the comparison of FWPR and FWPT as functions of time are shown in Fig. 12(b). Figure 12(c) illustrates the comparison of WCT as a function of time. According to these figures, satisfactory matches are obtained for all the parameters even with the blind data. The results confirm the efficiency of the developed multiple SVRs proxy model in forecasting the numerical simulator performance.

WAG CO 2 process optimization using the hybridization proxy GA
After the validation of the multiple SVRs proxy through the blind runs and checking its performance, it is coupled with GA to find the optimum parameters of the formulated problem shown in the fourth section. The first step in the optimization with SVR-GA is to codify the problem variables (their combinations represent WAG CO 2 ) in the form of chromosomes. The representation (IT, FWIR, FGIR, HCT) is followed. Then, an initial population of the chromosomes is generated randomly. The evaluation step is performed based on FOPT and the state of the constraint. Accordingly, a linear ranking based fitness function is implemented. After evaluation, the GA operators, i.e., selection, elitism, crossover and mutation are applied to bring out offspring. Finally, by checking the maximum number of iterations as a stopping criterion, the best fit individual is chosen as the optimal WAG CO 2 scenario. Table 10 presents the GA control parameters considered in the formulated WAG CO 2 optimization problem. Figure 13 shows the optimization results of GA in its best run. The optimum FOPT is 11.79129 × 10 6 sm 3 , and the corresponding design parameters are listed in Table 11. Table 12 demonstrates the simulator values and the errors of the parameters needed in the optimization, gained from the multiple SVRs proxy model results shown in Table 11. The outputs of the multiple SVRs proxy model corresponding to the optimum WAG CO 2 design parameters are in good agreement with the results of numerical simulations.
According to Table 11, optimum conditions occur when the WAG is initiated after six years, the gas is injected at low rate and water is injected at high rate. Additionally, the half-cycle time should be set to six months. These results show that a high WAG ratio and small slug sizes are favourable in this process.
Finally, the main advantage of establishing a proxy is the noticeable reduction in the simulation time. One reservoir simulation run in a machine including an Intel R Core TM i5-5200 2.20 GHz and 6 GB of RAM, takes 15 to 23 minutes (depending on convergence time), while the developed multiple SVRs proxy run takes only 1.086 seconds.

Conclusions
In this work, multiple SVRs were designed and applied as time-dependent proxy for optimizing a constrained WAG CO 2 process in real time. The proxy was trained utilizing a database, which was exploited from each time-step input/output system of insignificant number of runs. Blind runs were then considered to check the reliability with unseen datapoints. The established proxy provided the desired outputs by conserving the required accuracy within reasonable time.
The hybridization of GA with the multiple SVRs proxy allowed the investigation of optimum WAG design parameters with high accuracy and reasonable CPU time. A middle-time IT, gas and water injection HCT of 6 months, high water rate and low gas rate with high WAG ratio and small slug size are the proper operational conditions for the WAG CO 2 process formulated in this study.