Multi-Objective Optimization of Microemulsion Flooding for Chemical Enhanced Oil Recovery

— Microemulsion ﬂ ooding is one of the most effective methods of Chemical Enhanced Oil Recovery (CEOR), particularly for the production of residual oil trapped in unconventional reservoirs. A criticalstepforsuccessfulapplicationofthistechniqueistoachieveasuitableformulation.Previousstudieshavealmostfocusedonthetechnicalaspectswhileconsideringbothpracticalandeconomicmattersascon ﬂ icting objectives have been neglected. In the present paper, the formulation of microemulsion is optimizedbasedonthetrade-offbetweenscienti ﬁ cand ﬁ nancialresponsesusingahybridwork ﬂ owinwhich experimental design and arti ﬁ cial intelligence methodologies are composed. To appraise the ef ﬁ ciency of developed algorithm, a challenge case study is ﬁ rst evaluated and compared to previous approaches. Thereafter, the second case is examined in which a newly developed formulation of microemulsion for high temperaturecarbonatereservoirsisoptimized.Theoutcomesofthismulti-attributework ﬂ owarecompared to a single-objective algorithm.Theresults indicatetheoutstanding performance of the proposed approach formulti-objectiveoptimizationofmicroemulsionformulation.Eventually,thepossibleconcernsregardingtheapplicationofmicroemulsion ﬂ ooding in unconventional reservoirs are discussed.


INTRODUCTION
An immense amount of oil is remained unrecovered in the reservoirs after primary and secondary mechanisms of oil recovery. Generally, it exceeds more than half of the Original Oil in Place (OOIP) in almost cases [1,2]. This considerable potential is of great interest for Enhanced Oil Recovery (EOR) methods which fall into three main categories: thermal, gas and chemical methods. Generally, they involve the injection of fluids into the reservoir to generate favorable fluid properties or interfacial conditions, overcoming forces trapping the oil and boosting oil production substantially [3][4][5].
Microemulsion flooding is one of the high-performance techniques of Chemical Enhanced Oil Recovery (CEOR) methods. Not only does it have different advantages such as ultra low Interfacial Tension (IFT) and/or better sweep efficiency compared to water or surfactant flooding but it also can overcome some important drawbacks of conventional CEOR processes; for instance, high adsorption of surfactants and difficulties with control over in situ gelation of polymers in the reservoirs. The results of experimental as well as modeling studies demonstrate microemulsion flooding is one of the most promising solutions to achieve incremental oil recovery [6][7][8].
The optimization of influential parameters is of vital importance in the successful application of any CEOR technique particularly microemulsion flooding because not only does it improve its efficiency, but it also reduces the chemical costs and their possible harmful environmental impacts [9][10][11] owing to the fact microemulsion flooding is a complex and expensive operation [12]. For this purpose, different methodologies have been proposed in the literature.
A reservoir simulator with an economic model was used by Wu et al. [13] to optimize a chemical flood of a typical sandstone onshore oilfield. Zerpa et al. [10] developed a surrogate-based optimization methodology to estimate the optimal values for a set of design variables (concentration and slug size of chemicals). Then, Carrero et al. [11] employed this approach to study a global sensitivity analysis. Anderson et al. [14] optimized chemical flooding for a mixed-wet dolomite reservoir using sensitivity analysis of key parameters. Prasanphanich et al. [15] optimized Net Present Value (NPV) of chemical flooding using in-house simulators coupled with experimental design methodology. Douarche et al. [16] utilized response surface methodology to develop a surrogate model by which oil recovery from a surfactant-polymer flood was estimated and optimized. AlSofi and Blunt [17] simulated rheological behavior of polymer as a non-Newtonian fluid using streamline approach to evaluate polymer-flooding design.
Previous studies thoroughly discussed and optimized the technical aspects of CEOR approaches. They have mostly considered oil recovery or NPV as a singleobjective optimization process. Seldom have they paid attention to both practical and economic features. When a CEOR operation is going to be efficiently implemented, different aspects of technical and economic factors should be simultaneously analyzed. In this regard, the prediction of optimum values for control variables considering the multiple conflicting criteria is known as multi-objective optimization. In such problems, optimal decisions need to be taken in the presence of trade-off between two or more competing objectives [18,19]. Methodologies for multi-attribute optimization of CEOR or even other EOR approaches have rarely been reported in the literature.
In this study, a robust algorithm for multi-criteria optimization is developed based on the design of experiments and artificial intelligence techniques. The proposed workflow is comprised of two main stages: modeling and then optimization. The modeling operation establishes a mathematical relationship between design variables and corresponding responses of the process. For this purpose, Response Surface Methodology (RSM) is employed to attain precise objective functions. In the next stage, Particle Swarm Optimization (PSO) is coupled with Fuzzy Logic (FL), generating a multi-criteria optimization methodology to find the best solution of the problem. To confirm the efficacy of this algorithm, it is applied for modeling and multi-objective optimization of microemulsion flooding in different two case studies.

The Advantages of Microemulsion Flooding
Demonstrated the outstanding performance of CEOR processes have the results of numerous experimental and numerical studies. Consequently, many thriving CEOR projects have been reported worldwide. Although there is no doubt about the inherent efficiency of CEOR techniques [20][21][22][23], several issues have threaten its successful application in unconventional resources, such as high price of chemicals, harsh condition of oilfields, and possible formation damage.
In recent years, the shadow of bad economic conditions and low price of oil have fallen across the development of EOR projects. It matters a great deal to CEOR techniques because they are known as expensive approaches among other EOR operations. Therefore, decision makers should have multiple perspectives rather than single-objective insight in which one aspect (e.g. technical) is considered regardless of other features (e.g. economic). In this atmosphere, the application of multi-purpose algorithms paves the way for the optimization of control variables in terms of conflicting objectives, improving the efficacy of the operation and simultaneously managing the limited budget. In this paper, a workflow for multi-attribute optimization was successfully proposed in two case studies of microemulsion flooding. The application of this efficient methodology would be extended beyond different cases of CEOR or even other EOR techniques if sufficient data can be gathered.
The next concern is about the development of unconventional reservoirs. High potential of CEOR processes are always attractive to such reservoirs where thermal stability and ionic compatibility of chemicals are permanent restrictions. This article deals with microemulsions which are defined as thermodynamic-stable dispersion of an immiscible liquid into another one. The feasibility of the successful application of microemulsions in high temperature reservoirs can be vindicated by considering two concepts: spontaneous formation and thermodynamic stability of microemulsions. According to the second law of thermodynamic, the formation of microemulsions is modeled using Gibbs free energy: where DG F is the change in Gibbs free energy, g is IFT, DA is the change in surface area (A), T is temperature and DS is the variation of configurational entropy [24]. A colloidal dispersion is formed spontaneously when the variations of Gibbs free energy become negative. To fulfill the conditions, ultra-low interfacial tension due to the presence of an appropriate surfactants required [25]. Consequently, gDA < À TDS and then DG F < 0. Thus, the surfactant plays the crucial role for the microemulsion formation and the highenergy external sources for mixing two separate phases does not require. It saves the energy and therefore reduces the expenditure of formulation compared to macroemulsions and nanoemulsions which usually require external sources of energy such as ultrasonic waves or stirring [26].
Accordingly, thermodynamic stability means the free energy of microemulsions is lower than the free energy of the reference state (oil and water phases). It indicates the system tends to continue its stable position whenever the thermodynamic conditions (e.g. temperature) remain unchanged. When a microemulsion is formed in a particular temperature, the thermal stability of the mixture is guaranteed until the temperature is kept constant. Therefore, this concept can be used as an efficient mechanism to retain the unlimited stability of microemulsion slug at high temperature. In contrast, macroemulsions and nanoemulsions are thermodynamically unstable although they can be kinetically stable which defines as their tendency to persist in a meta-stable position for long time rather than a limitless period [27].
For CEOR purposes, injection fluids should be stable for at least several months. The injection of a thermodynamically stable microemulsion ensures that the chemicals can keep their original properties and work with highest efficiency during the time interval between injection and production wells. However, the instability always threatens other types of dispersions. This idea has been successfully developed in our previous study in which the formulations of the microemulsions were stable for several months while the phase behavior, characterization of optimum mixture, and oil recovery experiments in carbonate rocks were being done [28].
Another restriction of CEOR applications in unconventional reservoirs is the ionic compatibility between the chemicals and sea water or formation water. The mixture of surfactant should be compatible with the composition of these samples of water.
The last point is the evaluation of formation damage due to the injection of microemulsions as they are colloidal dispersions. Insofar as the microemulsions contain the droplets of the interior phase dispersed in the continuous phase, the size of particles should be set so that they can prevent the plugging of pores and throats which causes the formation damage.

Response Surface Methodology
RSM is comprised of a collection of statistical and mathematical techniques. It can be used when several associated factors influence a corresponding response or a set of responses in a process. Compared to trial-and-error approaches of experimental design in which a factor is changed while others are constant, RSM has different advantages. The effects of factors on the process response(s) can be studied through varying the levels of factors simultaneously. Hence, not only does RSM provide the optimum experimental design to cover all the sampling space with minimum required experiments (saving time and the budget) but it also can fit an appropriate model to the data [29]. In the other words, it can develop a functional relationship between a number of related parameters (as the inputs of the process) and the response(s) through statistically analyzing their interactions. Furthermore, it can be applied for the optimization of the process to determine its best performance. Therefore, it is an efficient protocol for Design of Experiments (DOE) [29][30][31][32][33].

Artificial Neural Network
Artificial Neural Network (ANN) is an information processing paradigm which is inspired by biological neural networks. It consists of highly inter-connected elements that are known as neurons [34]. As a parallel, distributive and adaptive system, it can be effectively applied in complex and non-linear problems in different engineering disciplines [35][36][37]. Nowadays, ANN is a computational tool in petroleum engineering studies because not only it is a model-free function predictor but also it does not require any detailed knowledge regarding the process [38][39][40].
Multi-Layer Perceptron (MLP) network in accompaniment with back propagation training algorithm is known as the most popular paradigm for the construction of the network. MLP structure comprises different layers (input, hidden and output). Different number of neurons is placed on each layer. The connectivity between each neuron of the current layer and those on the proceeding layer is established through direct links that have their own weights. In the training phase, these weights are changed to satisfy the training goal. The output layer collects the received signals from previous layers and creates the response of the network; in other words, ANN models the process. Theoretical outlines of ANN have been discussed elsewhere [41].

Particle Swarm Optimization
Particle Swarm Optimization (PSO) was introduced by Kennedy and Eberhart in 1995 [42]. The idea was originally derived from mimicking social behaviors exhibited by group animals such as bird flocking and fish schooling. It is a population-based stochastic search like other biologically inspired optimization paradigms [43]. PSO is an attractive optimization methodology because not only can it be simply implemented with low computational efforts owing to the fact that there are few parameters to be adjusted, but it also adopts real numbers as particles; in other words, it does not need binary encoding or evolutionary operations which are required for Genetic Algorithm (GA). Furthermore, PSO utilizes collaborative approach (particles have memory to maintain knowledge of appropriate solutions) compared to GA in which a competitive one is used. Thus, PSO is very successful to solve a wide variety of problems of diverse engineering as well as computational science in continuous search space [44][45][46].
PSO uses a set of search points as particles in a population which is known as swarm. Two variables characterize the particles: the particle position (x) which is representative of a potential solution to the problem and particle velocity (v). It is initiated with a swarm of particles as random solutions distributed in the search space. Each particle moves through search space and compares its fitness value at the current position to the best fitness value so far attained. The best solution so far achieved by a particle is known as the personal best, or pbest. Accordingly, when PSO proceeds, global best or gbest is defined as the best global experience among all of the pbest positions achieved so far. The gbest is the final solution of the optimization problem. The stopping criterion can be the number of maximum iterations which is defined by user.
The definitions of pbest and gbest are employed to update the ith particle velocity in the dth dimension: where w is inertia weight factor (controlling the effect of previous velocities on the current velocity), C 1 and C 2 coefficients are acceleration constants which help particles to accelerate towards favorable areas of the search space. The function of rand(0,1) generates a random number between 0 and 1.
According to the updated velocity (Eq. 2), each particle repositions as follows:

Fuzzy Logic
The theory of fuzzy logic was introduced by Zadeh in 1965 [47]. It has been developed to facilitate approximate reasoning by which human expert's reasoning process is mimicked; in other words, a methodology to compute with words is provided using fuzzy logic. It is made up of the fuzzy sets [48]. In a crisp set (traditional Aristotle's logic), an object either belongs to this set or it does not. In a fuzzy set (Eq. 4), however, the transition from membership to nonmembership is gradual [49][50][51][52], that is to say any object belongs to a fuzzy set to a certain degree: where F is fuzzy set, x is a generic element of X as a space of points, and m f ðxÞ is membership degree which quantifies the belonging grade of x to the fuzzy set [53][54][55].
The concept of partial truth in fuzzy sets is represented by membership function which maps each object of X to a membership value between 0 and 1. There are different types of membership functions, such as Triangular, Trapezoidal and Gaussian.

Data
Two case studies of microemulsion flooding were considered (Appendix A). The conditions and procedures to produce experimental data have been described elsewhere [28,56]. Table 1 contains the technical and economic information relevant to the chemicals used in their formulation. Moreover, the oil price was set to $ 50.

Case 1
The first mixture was a triglyceride microemulsion [57] demonstrating the outstanding performance for CEOR applications in sandstones [58]. It was utilized in the first case study to confirm the efficiency of proposed algorithm via the comparison of results obtained by current methodology with those extracted from previous study [57]. The original dataset contained 160 data each of which comprised four inputs: the concentration of surfactant, salinity, co-surfactant concentration, and water content in the microemulsion phase. Interfacial Tension (IFT) of microemulsion and Recovery Factor (RF) were the outputs. The range of input and output parameters is presented in Table 2.

Case 2
The second case study was a newly developed microemulsion formulation to be applied for CEOR in high temperature carbonate reservoirs [28]. Four inputs were included in each data: surfactant concentration, salinity, cosolvent concentration and volume percentage of initial oleic phase. The Solubilization Capacity (SC) of oil phase into aqueous phase and total Cost of Microemulsion Production (CoMP) were the corresponding responses. The domain of variables for this case study is tabulated in Table 3.

Algorithm for Modeling and Optimization
An algorithm for multi-objective modeling and optimizing has been developed ( Fig. 1). Two main stages are included in the procedure. The former provides objective functions through modeling of the microemulsion flooding using RSM and the latter optimizes the process as a multi-objective problem using PSO and FL. The detailed description of the algorithm is presented as follows:

Process Modeling
The current optimization algorithm requires precise objective functions. In this regard, the process was modeled to establish mathematical relationships between inputs and outputs of the process. For this purpose, RSM was employed by which statistically significant quadratic equations can be developed. In the first place, a series of n experiments was designed. RSM consists of different approaches. The selected designs were dissimilar for two case studies. The former used D-optimal design while the latter utilized Central Composite Design (CCD). Both plans are the most popular RSM methods [59]. Considering four effective factors (Tabs. 2 and 3) with six center-point runs (to calculate experimental error and reproducibility of data), Doptimal design proposes 50 runs (n = 50) while CCD suggests a 30-run plan (n = 30). In each run, an arrangement of four independent variables in different levels is set.
In the next step, the corresponding response (output) for each run should be provided. The routine approach to provide the response(s) is their measurement in the laboratory. For newly developed formulation of microemulsion in high temperature (case study 2), solubilization capacity of each run was already on offer from our previous study [28]. Accordingly, the cost of microemulsion formulation was calculated based on the information of Table 1. However, the pattern of experimental design (Doptimal) differs from that of original dataset (Factorial) [56] in case study 1. The determination of 2 Â 50 required responses experimentally necessitates an immense amount of time and considerable experimental expenditures. To try to deal with address this problem, a reliable solution was substituted whereby the responses were determined numerically via the development of a proxy model. For this purpose, artificial neural network as a data-dependent To implement ANN modeling, the original dataset of experimental data was randomly divided into two separate datasets: training (60% of the data) and evaluation (the remaining 40%) datasets. 50 runs were put in a new dataset called testing. Inputs (four independent variables) and outputs (two dependent variables) of each individual data in different datasets were normalized using mean and standard deviation of variables. The process of ANN modeling was divided into three phases: training, generalization and operation.
The training phase was implemented to transfer expertise and derive functionality between inputs and outputs of the process which was generally unknown. For this end, a multilayer back propagation network was selected. This network consisted of 3 layers: input, one hidden and output. Tangent sigmoid and purelin types were assigned as transfer functions of hidden and output layers, respectively. Furthermore, Levenberg-Marquardt was the default training algorithm and the default number of neurons in hidden layer was two times the number of control variables, namely eight neurons.The performance goal or the stopping criterion of this stage was Mean Square Error (MSE) of 1 Â 10 À3 . The training dataset was fed to this phase.
In the generalization phase, the optimum structure of the ANN was achieved based on the evaluation dataset. To that end, different number of neurons in the hidden layer from 1 to 15 as well as 10 training algorithms were evaluated which were as follows: Levenberg-Marquardt (trainnlm), BFGS quasi-Newton (trainbfg), Bayesian regularization (trainbr), Powell-Beale conjugate gradient (traincgb), Fletcher-Powell conjugate gradient (traincgf), Polak-Ribiere conjugate gradient (traincgp), Gradient descent with adaptive learning rate (traingda), Gradient descent with momentum (traingdm), Gradient descent with momentum and adaptive learning rate (traingdx).
Three statistical parameters of Mean Absolute Percentage Error (MAPE), Symmetric MAPE (SMAPE) and Root Mean Square Error (RMSE) quantified the accuracy of training and generalization phases: j y p À y e y e j; ð5Þ The workflow for modeling and multi-objective optimization.
jy p À y e j jy p jþjy e j 2 ; ð6Þ where n is the number of data in the dataset, y p is predicted response by ANN and y e is experimental data.
In the operation phase as the last step of surrogate modeling, ANN with optimized structure (the best performance) was utilized to extract the outputs of testing dataset. It was the end part of ANN modeling by which required responses of 50 runs in D-optimal design for case study 1 were provided.
There was another problem with the responses of case study 1. They were not independent as whatever IFT is reduced RF increases. Therefore, they could not be considered as conflicting objectives in multi-attribute optimization methodology. To address this issue, RF was replaced with an economic output which was Chemical Efficiency (CE). It was defined as the total oil revenue calculated by RF at oil price of $ 50 divided by the total cost of chemicals quantified based on the information of Table 1. IFT and CE were separately influenced by technical and economic features of the process, respectively. That is why they were taken into account as conflicting objectives for multi-attribute optimization of this case study.
For both case studies, responses were finally fed to RSM designs which analyzed the obtained data statistically through Analysis of Variance (ANOVA) to establish relationships between effective factors (inputs) and corresponding responses (outputs). In other words, RSM developed second-order equations as objective functions required for multipurpose optimization algorithm.

Multi-Criteria Optimization Methodology
Having objective functions obtained from previous stage, multi-objective optimization of microemulsion formulations was performed using PSO-FL methodology. For both case studies, low and high factorial points were selected as the extremes of each control variable in optimization process (Tab. 1) except to salinity that was set to a constant value because the salinity of the injection water in which chemicals are dissolved is constant. Insofar as the basic fluid for almost water-based EOR techniques is sea water, the salinity of both formulations in two case studies was fixed at 4.2 wt% which is the typical salinity of sea water [60]. Accordingly, the optimization goal was set to maximize the technical efficiency and minimize the cost of formulations, simultaneously.
In the first iteration, a swarm of 25 particles was randomly positioned in the search space. The constant inertia (w) of 0.7298 and acceleration coefficients (C 1 and C 2 ) of 1.49618 were assigned because empirical results have indicated that The effects of different transformations on the ANN prediction in training phase of case study 1: (a) interfacial tension and (b) recovery factor. they help better convergent behavior [61]. Then, the components of velocity factor were initiated. Following this, the corresponding responses of each individual particle were calculated using objective functions.
To solve a multi-attribute problem, multiple objectives can be combined into one single-objective function. Moreover, finding pbest of each iteration necessitates the definition of a unique objective function. For this purpose, PSO was coupled with fuzzy logic. First, fuzzy membership functions of two objectives were generated in which the mathematical objective functions are treated as constrains (Eq. 8). In other words, each fuzzy membership function represents a fuzzy region of acceptability.
where m f ðF k Þ is fuzzy membership function of k, the objective function (k = 1, for the first and k = 2, for the second response of each case study), F k is the kth response, F min k and F max k are minimum and maximum values for each objective function, respectively. Following this, the new objective function can be defined as a function called zeta: where z j is satisfaction factor and the subscribe j = 1 :25 denotes the particle number. For two case studies, m f ðF 1 Þ was practical (IFT or SC) while m f ðF 2 Þ was economic (CE or CoMP) fuzzy membership functions, respectively. The best solution of current iteration (pbest i ) has the maximum value of particles' zeta: In the next iteration, a new swarm containing random particles is generated whose positions and velocities are updated according to equations (1) and (2), respectively. Finally, the multi-criteria solution is gbest: where gbest is the best global solution and pbest i is the best local solution of ith iteration (i = 1:60).

Confirmation of Proposed Methodology in Case Study 1
Considering multiple criteria, different methodologies were employed to develop an efficient algorithm for modeling and optimization of the process. To evaluate this workflow, two case studies of microemulsion flooding were considered. As shown in Figure 1, the algorithm was commenced with RSM (D-optimal method) by which 50 experiments were designed. Then, the responses have been provided. As mentioned before, the responses were numerically calculated using ANN paradigm whose results are presented and discussed as follows:

Modeling of Microemulsion Flooding Using ANN
ANN was used to model the production of triglyceride microemulsion, performing as a numerical simulator to provide the responses of D-optimal design. The neurosimulation process was started with training of a multi-layer back propagation network using training dataset. Error analysis of this phase is presented in Figure 2. Error analysis of MAPE and SMAPE for IFT response were 2833.1% and 107.8%, respectively ( Fig. 2(a), blue column) while they were calculated as less than 2% for RF response (Fig. 2(b), blue column). ANN was not trained for IFT very well although the results were acceptable for RF.
Comparing the predicted values of IFT with real ones in training dataset, it appears the huge error of MAPE has stemmed from two points. First, we found ANN predicted several negative IFT which are conceptually impossible for such oleic/ aqueous systems. Second, the training precision was the worst for low and ultralow IFT. To improve the accuracy of neurosimulation, the original dataset should be analyzed owing to the fact that ANN is an information-processing methodology and it is strongly dependent on the original data. For this end, the spatial distribution of original responses is plotted in Figure 3 (a) using a 3D plot [62]. It reveals that IFT did not lie in a standard distribution pattern because there is a marked difference between the minimum and maximum values of IFT (0.0001 IFT 11.2530). By contrast, RF was expanded over a limited domain. Since IFT was varied over several orders of magnitude it could not be trained very well.
To modify the wide distribution of IFT, some transformations were examined. First, it was replaced with ln (IFT) to restrict its range. In addition, solubilization ratio (S) using Chun Huh equation and ln(S) were assigned as new transformations: where c is constant (usually 0.3). To compare the spatial distribution of IFT original dataset with its transformed forms, the Standard Deviation (SD) as a useful indicator was employed. It evaluates the dispersion of the dataset: where x i is the ith sample of the dataset, m denotes the mean of value, and n is the number of samples.
The values of standard deviation for original and transformed datasetsof IFTwereasfollows:IFT = 1.97,ln(IFT) = 2.39,S = 6.77 and ln(S) = 1.19. They implies ln(S) provided the best distribution among others as the variations of IFT was modified very well (À1.8 < ln(S) < 3.7); compared to original IFT (0.0001 IFT 11.2530). The spatial distribution of ln(S) is visualized in Figure 3(b) which demonstrates its uniform distribution and narrower range of variations over the search space.
The effects of different transformations on the performance of the trained network are presented in Figure 2. Replacing IFT with ln(S) improved the quality of training significantly. MAPE and SMAPE were considerably reduced from 2833.1% to 65.7% and from 107.8% to 49.8%, respectively.
It is noteworthy that error analysis based on just one parameter (either MAPE or RMSE) can be misleading [63][64][65]. For example, RMSE of IFT (without any transformation) was 0.1215, compared to ln(S) which was 0.4602. A brief glance confirms that the transformations did not affect serious improvement if only RMSE was considered. However, MAPE and SMAPE of IFT were extremely high compared to ln(S). That is why the efficiency of training was examined based on three error estimators.
Accordingly, the use of SMAPE beside MAPE is essential because MAPE can be varied between zero and infinity while the upper endpoint of SMAPE (Eq. 6) is restricted to 200% [66,67]. In the current case, MAPE of 2833.1% does not make sense for engineers but they can deal better with SMAPE of 104.8%. Although the application of ln(S) boosted the performance of ANN training, current values of MAPE (66.7%) and SMAPE (49.8%) indicate that accuracy of ANN modeling should be further improved. In the generalization phase, attempts were made to increase ANN accuracy.
The efficiency of the trained network to predict a newly presented unseen data (evaluation dataset) was gradually improved during generalization operation in terms of the training algorithm, number of neurons, and performance goal.
In this step, 10 training algorithms were tested. Figure 4 demonstrates that the most efficient algorithm is Bayesian regularization back propagation by which overfitting can be prevented [41]. MAPE of ln(S) decreased substantially from 43.6% to 1.86% and SMAPE declined from 39.1% to 1.87% when Levenberg-Marquardt algorithm (default algorithm in training phase) was replaced with Bayesian regularization. Similar to ln(S), SMAPE of RF was reduced from 0.74% to 0.05%.
Following this, the optimal number of neurons was selected. For this end, it was changed from 1 to 15 and the variations of RMSE were followed. Figure 5 shows that the optimum neuron number was 10 since RMSE of ln(S) and RF were the least. Finally, the performance goal of training was reduced to 1 Â 10 À5 to minimize the prediction errors.
Thanks to the successful transform function of original response (IFT) to ln(S) and efficient structure of the network optimized during generalization operation, the accuracy of the ANN predictions improved significantly as indicated in Table 4. The overall errors of responses have been almost reduced less than 1%, demonstrating the efficiency of ANN to be successfully used for the accurate prediction of unknown data (has not been introduced to the network). Graphically, the crossplots of ANN estimation and original data in training and generalization phases are shown in Figure 6.
Eventually, the trained network with optimum structure was ready to be used as the surrogate model to produce 50 pairs of responses required for RSM. It is called operation phase in which new inputs are introduced to the network and The spatial distribution of outputs in case study 1: (a) IFT and RF (b) ln(S) and RF.
ANN can estimate output(s). Here, it was utilized to calculate the required responses of RSM. 50 runs of experimental design as testing dataset were fed to ANN and the corresponding responses were provided.

Generation of Objective Functions Using RSM (D-Optimal Pattern)
The provided responses from the previous stage were imported into software to fit quadratic models using RSM, developing objective functions. Statistical analysis of proposed mathematical models was tabulated in Table 5. Based on the results of ANOVA, the F Values for ln(S) and CE models were 139.00 and 1949.15, respectively. Moreover, the values of Prob > F were less than 0.0001 which indicates that both quadratic models were significant in 99.9% confidence interval. Other parameters of this table confirmed the performance of fitted equations; for instance, R 2 and R 2 adj of ln(S) were 0.9823 and 0.9753 while these control factors were calculated as 0.9987 and 0.9982 for CE. They demonstrated that equations fitted numerical data very well. The quadratic models (Appendix B) were extracted as objective functions of optimization process.

Multi-Objective Optimization
Having fitted quadratic equations, multi-attribute optimization was initiated. In each iteration of PSO, there was a swarm of 25 particles and optimization algorithm was continued to 60 iterations. The effect of training algorithms on the ANN prediction in generalization phase, (a) ln(S) and (b) recovery factor.

Figure 5
The effect of neuron number in the hidden layer on the accuracy of ANN predictions. In each iteration, position and velocity of every particle were randomly determined. Then, corresponding responses were calculated using objective functions (Appendix B). Following this, zeta function (Eq. 9) was used to find satisfaction factor of each particle, facilitating multi-objective optimization. The maximum value of z was assigned as pbest of curriculum iteration. Proceeding with the optimization process, z was increased. The variations of two original responses (outputs) during the progress of PSO-FL methodology are demonstrated in Figure 7 which The crossplots of ANN predictions versus actual data for training phase of (a) ln(S) and (b) recovery factor, generalization phase of (c) ln(S), and (d) recovery factor. indicates that the IFT decreased while the RF increased simultaneously until they reached their optimal values. Optimum arrangement of control variables by which gbest was achieved, are tabulated in Table 6. The algorithm proposed the optimum formulation of microemulsion as follows: surfactant concentration of 1.04 wt%, salinity of 4.2 wt%, co-surfactant concentration of 2.93 wt% and water content of 94 wt%. Performing this formulation, the algorithm was predicted ultra-low IFT of 4.01 Â 10 À4 mN/ m and chemical efficiency up to 13.50%. The trend of improvement for original responses is ploted in Figure 7.
To confirm the results achieved by PSO-FL methodology, RSM as a multi-purpose optimization technique was employed. The results are presented in Table 6. It shows that the optimum formulation of microemulsion had ultralow IFT of 2.81 Â 10 À4 mN/m and chemical efficiency of 13.18%. The results of PSO-FL were in consistent with those of RSM. Compared two methodologies, one can easily prove the efficacy of PSO-FL. Moreover, the results of PSO-FL were compared to single-objective optimization of each individual response. It indicates the trade-off between conflicting objectives could not be extracted when singleobjective algorithms were used which may mislead the engineers in such multi-attribute problems.

Case Study 2 3.2.1 Development of Objective Functions Using RSM (CCD Plan)
According to CCD approach, 30 experiments were designed. Insofar as the required responses were exactly available in the literature, they were gathered and directly fed into the RSM. In other words, the implementation of neuro-simulation step in the workflow (Fig. 1) did not necessitate for this case study. Anyway, RSM employed different control statistical parameters to develop precise quadratic models fitted to experimental data. Based on the results of ANOVA, the F Value of responses for solubilization capacity and cost of chemicals were calculated as 46.41 and 39.08 while Prob > F were less than 0.0001 for both responses. These values indicate the developed models were statistically significant in 99.9% confidence level (p-value < 0.001). R 2 , R 2 adj were calculated as 0.9659 and 0.9451 for solubilization capacity response. On the other hand, they were assessed at 0.9598 and 0.9353 for the second output (the cost of microemulsion production). Other statistical parameters are presented in Table 5. The results of ANOVA table pointed to favorable fits of quadratic equations (Appendix C) to experimental data.

Multi-Objective Optimization
The fitted quadratic equations were used as objective functions in the PSO-FL algorithm whereby the trade-off between technical and economic objectives was considered. Table 7 shows the optimal formulation of microemulsion. When aqueous phase (60 v%) consisting of 3.48 wt% surfactant and 4.2 wt% NaCl was mixed with oleic phase (40 v%) comprising 86.73 v% biodiesel and 13.27 v% cosolvent, the solubilization capacity reached up to 8.69% and the cost of microemulsion formulation was simultaneously reduced to 9.93 $/bbl.
The prominent role of PSO-FL methodology to consider trade-off between technical and economic matters of a process (here microemulsion formulation) was revealed when the best estimation of PSO-FL workflow was compared to PSO algorithm as a singleobjective paradigm.
When the first response (SC) was solely considered, the optimal concentrations of surfactant and co-solvent were 5.62 wt% and 15.00 v%, respectively. Although the solubilization capacity increased from 8.69% (PSO-FL: SC&CoMP) to 11.30% (PSO: SC), 22.7% rise in the cost of formulation was observed (from 9.93 $/bbl to 12.18 $/ bbl). On the other hand, the cheapest formulation was achieved when the single-objective PSO algorithm was run just based on the cost of chemicals. In this case, the optimum concentrations of surfactant and co-solvent were predicted as 1.78 wt% and 5 v%, respectively. Although the formulation cost was substantially reduced by 33.6% (from 9.93 $/ bbl to 6.59 $/bbl), poor quality of formulation was observed because solubilization capacity decreased drastically from 8.69% to 3.16%. A remarkable difference between singleand multi-objective optimization algorithms was explicitly indicated. Therefore, competing objectives necessitate the application of multi-attribute methodologies unless trade-off decisions cannot be made.
In the background section, several issues were discussed regarding the application of microemulsion for CEOR applications. The last two ones were salinity and formation damage. There was a wide range of salinity in the experimental evaluations of microemulsions in both case studies (Tabs. 2 and 3) which include many compositions of The improvement of two original responses toward their optimum values during PSO-FL methodology in case study 1. sea water and even formation water. The perturbation plot of salinity in these domains is presented in Figure 8 in which the effect of salinity factor on technical responses is examined while other factors were set to their zero (middle) levels in this assessment. It indicates the salinity was not a deterministic factor because the variations of salinity did not significantly influence on the responses in both cases owing to the fact that nonionic surfactants (APG and Polysorbate 80) were used in the microemulsion mixture which is more compatible than ionic surfactants [28,56]. As a result, such formulations demonstrate high potential for CEOR operations in unconventional reservoirs; not only for unlimited stability at high temperature but also for high tolerance against high salinity.  Theoretically, the diameter of microemulsions is less than 100 nm. Generally, the exact size of particles in any dispersion is determined in terms of various parameters, such as the type and concentration of chemicals, thermodynamic conditions of the system, and the source of energy for mixing operation [68]. The Z-average diameter of the newly developed formulation was approximately 56 nm [28]. Table 8 compares it with throat size of different reservoirs. The average diameter of particles was almost smaller than throats in different porous media even some types of tight oil formations. Hence, such colloidal dispersions can be potential candidates for diverse range of reservoirs without any concern regarding the formation damage.

CONCLUSION
Better design of microemulsion formulation for CEOR purposes is dependent on the optimization of its formulation in the presence of multiple criteria. In this article, a hybrid algorithm was proposed for modeling and multi-criteria optimization of microemulsion formulation in two different case studies. It effectively combined the abilities of response surface methodology as well as artificial neural networks to model the process and then particle swarm optimization and fuzzy logic were coupled to optimize it. To confirm the effectiveness of this algorithm, it was applied on the triglyceride microemulsion in the first case study. Thereafter, it was used for a newly developed microemulsion for hightemperature carbonate reservoirs. Based on the results of these case studies, the following conclusions can be drawn: algorithms that can consider simultaneously two or more competing objectives are of great interest in the optimization of microemulsion formulation and CEOR processes in general because they regard the problem from different points of view. This article demonstrates that trade-off decisions can be made only when methodologies supporting multiple perspectives are used rather than uniobjective algorithms; different techniques were employed in the structure of the hybrid algorithm for the multi-objective optimization of microemulsion formulation in terms of technical and economic aspects. The workflow consisted of two main phases: modeling and optimization. To address the challenges of finding accurate objective functions as the outcome of modeling stage, different approaches were used: appropriate selection of RSM design, modification of the spatial distribution of responses, and optimization of ANN structure. Thus, the overall symmetric prediction error was reduced less than 2%. Furthermore, the results of analysis of variance in RSM indicated that developed quadratic equations as objective functions were fitted very well; the coupling of particle swarm optimization and fuzzy logic were found to be a powerful multi-attribute optimization method by which optimal microemulsion formulation was successfully achieved based on conflicting practical and financial objectives; the optimization results of microemulsion flooding were confined to lab scale in both case studies and cannot be extended to field scale as the available data had been gathered from experimental works due to the lack of field data. It does not reflected the limitation of proposed workflow for optimization of case studies in field-scale.
The introduced methodology can be successfully applied for modeling and multi-purpose optimization of any CEOR process in both lab and field scales if relevant database is collected; the microemulsion flooding demonstrates high potential to be used as an effective CEOR approach in different reservoirs.

ACKNOWLEDGMENTS
The authors would like to thank Dr. Z. Jeirani for provision of experimental data in case study 1. In addition, the authors are grateful to IOR Research Institute for the financial support. The authors are also thankful to Stat-Ease, Minneapolis for the provision of the Design Expert package.

APPENDIX A
The data of the first database is available through the reference [57]. The original data of the second case study is presented in Table A1.