Novel Self-Adaptive Shale Gas Production Proxy Model and Its Practical Application.

Recently, production optimization has gained increasing interest in the petroleum industry. The most computationally intensive and critical part of the production optimization process is the evaluation of the production function performed by the numerical reservoir simulator. Employing proxy models as a substitute for the reservoir simulator is proposed for alleviating this high computational cost. In this study, a new approach to construct adaptive proxy models for production optimization problems is proposed. An adaptive difference evolution algorithm (SaDE) optimized least-squares support vector machine (LSSVM) is used as an approximation function, while training is performed using a self-adaptive response surface experimental design (SaRSE). SaDE selects the optimal hyperparameters of LSSVM during the training process to improve the prediction accuracy of the proxy model. Cross-validation methods are used in the recursive training and network evaluation phases. The developed method is used to optimize the production of block gas reservoir models. Computational results confirm that the developed adaptive proxy model outperforms traditional regression methods. It is further verified that when the experimental data are updated, the alternative model still has high prediction accuracy when performing the objective function evaluation. The results show that the proposed proxy modeling approach enhances the entire optimization process by providing a fast approximation of the actual reservoir simulation model with better accuracy.

Introduction

The development of onshore oil and gas has made rapid progress over the past decade, primarily due to the development of unconventional resources, particularly shale gas.[1−7] Although shale gas extraction has broken through the industrial capacity barrier, it is at the critical point of marginal benefits, making it difficult to achieve beneficial development. Shale gas production optimization is getting increasing attention in the oil industry. Production evaluation is one of the most critical and difficult processes in production optimization under complex geological and engineering conditions. Reasonable and accurate production evaluation is essential for production optimization.

Numerical reservoir simulation technology, developed by reservoir geologists, reservoir engineers, and reservoir numerical simulators, is a physics-driven modeling approach recognized as the standard decision-making approach in the petroleum industry.[8] Numerical reservoir simulations can effectively simulate discrete fracture networks, fluid seepage, and key physical effects that could be important in shale reservoirs. Zhang et al.[9] applied a porous system with bedrock grid to describe the flow of unsteady gas from bedrock to fracture interval and used numerical simulation software to simulate the effects of parameters such as porosity, permeability, bedrock-fracture coupling parameters, and fracture half-length spacing on shale gas production capacity. Hu et al.[10] simulated the sensitivity of production capacity to parameters such as the number of fracture bars, fracture conductivity, and fracture spacing, which provided guidance for designing fracture construction plans for shale gas wells. Numerical simulations with higher accuracy have been frequently used as a powerful tool for engineering design and optimization. However, hundreds or thousands of reservoir simulations need to be run to perform optimization tasks. A single run of a simulated reservoir model consisting of thousands or even millions of grid blocks can take several hours. In addition, the large number of control parameters exacerbates this numerical simulation-based optimization design problem.

Essentially, the proxy model is a data mining tool. Advances in data mining, especially in machine learning, made it possible to create machine learning models that outperformed many existing computer models. Hassani et al.[11] used the traditional quadratic, multiplicative, and redial basis function (RBF) models to approximate the cumulative outflow from the reservoir when drilling a new horizontal well, which significantly reduced the computational time. Schuetter et al.[12] employed simple regression and other advanced methods such as random forests (RFs), support vector regression (SVR), gradient propulsion machines (GBM), and multidimensional kriging to build prediction models for production indicators. Zhang et al.[13] proposed a new multipoint adaptive Gaussian process agent model for design domain (MAGPSM-IDD) and applied it to an integrated field and an actual reservoir, showing superior computational performance. Wang et al.[14] employed a deep neural network (DNN) and Sobol global sensitivity analysis to analyze 2919 wells (2780 multistage hydraulically fractured horizontal wells and 139 straight wells) in the Bakken Formation. A thermal coding method was also employed to process the classified data. A reliable DNN model was built by evaluating Xavier initialization, exit techniques, and batch normalization. The results show that the proposed DNN model can integrate directly into existing hydraulic fracturing design programs. Brantson et al.[15] applied back-propagation artificial neural network (BPANN), radial basis function neural network (RBNN), and generalized regression neural network (GRNN) as proxy models to predict the historical production decreasing trend of extralow porosity tight gas reservoirs and further validated the effectiveness of the proxy models.

Although many studies have shown the broad application of data mining techniques in various engineering disciplines, its application in petroleum engineering is still in its infancy. Moreover, the performance of production forecasting in unconventional shale reservoirs is not impressive.[16−18] An important reason is that reservoir properties, including permeability, porosity, gas adsorption, gas saturation distribution, and reservoir pressure, as well as fracture stimulation design, will have a severe impact on shale gas production. More significantly, the prediction performance of machine learning algorithms is influenced a lot by their hyperparameters, such as learning rate, initial weights and thresholds, activation function, number of neurons, and the number of layers of a multilayer neural network, as well as the characteristics of the training samples. Hyperparameter tuning or optimization is the significant process of automatically testing different configurations for training a machine learning model. Prabusankarlal et al.[19] combined the adaptive difference evolution algorithm (SaDE) with extreme learning machine (ELM), where the SaDE algorithm optimizes the hidden node learning parameters of ELM to improve the model prediction accuracy and generalization performance significantly.

We used a self-adaptive response surface experimental design as the data set. Response surface design is a statistical method that uses reasonable experimental design methods and obtains specific data through experiments. Using multiple quadratic regression equations to fit the functional relationship between factors and response values, we seek optimal process parameters by analyzing the regression equations to solve multivariate problems.[20−22] If the training sample is limited, then the proxy model will be challenging to extract discriminative features automatically and guarantee regression accuracy. As the number of training samples increases, all of these estimates approach the true regression accuracy, which is the accuracy of a proxy designed with the full knowledge of the actual situation. This study will use a self-adaptive response surface design scheme, which generates training samples with specified accuracy. This adaptive accuracy model ensures that the training samples have actual reservoir characteristics, thus avoiding unsatisfactory learning results due to errors in the training samples. In addition, we introduce the least-squares support vector machine (LSSVM) as a proxy model. The LSSVM is functionalized by SaDE to optimally select the unknown parameters of the LSSVM and compare the results with other optimization algorithms (genetic algorithm (GA), particle swarm algorithm (PSO), gray wolf algorithm (GWO)) to generate a reasonable production proxy model. Our approach efficiently (in a short time) automatically searches for the best machine learning algorithm and hyperparameter values for a given machine learning problem.

Constructing the Proxy Model

Methodology of Least-Squares Support Vector Machines

The support vector machine (SVM) is a typical machine learning method based on statistical theory. It is effective in solving problems with a few samples, nonlinearity, and high dimensionality. With the advancement of the research, SVM is being widely used for complex nonlinear modeling problems in various fields. The main disadvantage of SVMs is the dimensional disaster problem. LSSVM is an SVM-based improvement type proposed by Suykens[23] in 1999. The higher computational load of SVM is overcome by LSSVM, which uses a set of linear equations to solve the problem.

The basic principle of LSSVM is to classify the hyperplane that can correctly partition the data set and has the largest geometric interval for a given data set. As shown in Figure , ω·x + b = 0 is the classification hyperplane. For a linearly divisible data set, there are infinitely many such hyperplanes, but the classification hyperplane with the largest geometric interval is the only one.

Figure 1

Schematic diagram of basic principles of LSSVM.

This study utilizes the regression form of LSSVM as the proxy model. LSSVM can transform linearly inseparable data into linearly separable data using kernel functions to map them to a high-dimensional space. The function of the kernel function in LSSVM is shown in Figure .

Figure 2

Function of the kernel function in LSSVM.

For the training sample set S = {(x, ···, y)} with size n, consisting of the input value x ∈ R and the corresponding output values y ∈ R, its LSSVM regression model is determined as followswhere ω is the weight vector, b is the offset, and the nonlinear transformation φ(x) is the mapping from the low-dimensional space to the high-dimensional space. According to the principle of structural risk minimization, in the LSSVM, the squared error loss function in the optimization targets is selected so that the regression problem is transformed into a quadratic optimization problem as followswhere e is the slack variable and γ is the penalty. The constraint condition is as followsTo solve the optimization problem, the Lagrange function was introducedwhere α represents the Lagrange multiplier. The following formula could be obtained based on the Karush–Kuhn–Tucker (KKT) optimization conditionswhere K(x,x) refers to the kernel function. In this study, the radial basis function was selected as the kernel function of the LSSVM regression model for its excellent generalization ability and wide-range convergence domain. Its expression is as followswhere σ refers to the kernel function width; γ and σ have a significant degree of impact on the learning and generalization ability of LSSVM models. Therefore, optimization algorithms (e.g., PSO, GA, SaDE, GWO) are generally applied to select the hyperparameters γ and σ of the LSSVM model comprehensively and optimally.

The quadratic optimization problem can be transformed into the problem of linear equationThe ultimate LSSVM regression function was built as follows

Methodology of Parameter Optimization of the LSSVM

Differential evolutionary (DE) algorithm is an excellent swarm intelligence optimization algorithm after genetic algorithm (GA) and particle swarm algorithm (PSO). The DE algorithm has been extensively applied because of its simple structure, few control parameters, easy implementation using accurate number coding, and fast convergence speed.

Differential Evolution Algorithm

The basic procedures of the DE algorithm include initialization, mutation, crossover, and selection. Individuals are randomly generated in the search space by initialization; then, mutation and crossover generate new individuals and eventually select which individual will enter the next generation. This process will continue to repeat until its termination condition is reached.

Initialization

For minimizing the objective function concerning the parameter vector x ∈ R, the population containing NP parameter vectors evolves toward the global minimum. The ith individual is as followswhere i = 1, 2, ···, NP. G refers to the number of evolutionary generations.

Initialize to generate a set of NP parameter vectors x to cover the parameter space as much as possiblewhere xmin = [xmin1, xmin2, ···, xmin] and xmax = [xmax1, xmax2, ···, xmax] refer to the maximum and minimum parameters boundary, respectively.

Mutation

The existing vectors are converted into new mutant vectors by a self-organization scheme that takes the difference vectors of randomly selected population vectors. For each parameter vector x of the current generation, a mutation vector V is generated by some mutation strategy. Several popular mutation strategies (MT) are listed as followsMT1: DE/rand/1MT2: DE/rand/2MT3: DE/rand-to-best/2MT4: DE/current-to-rand/1where F refers to the positive amplification factor, which controls the scaling of the difference vector and is usually chosen in the range 0 ≤ F ≤ 2. K refers to the control parameter, which is randomly generated in the range 0 ≤ K ≤ 1. Indices r are r1, r2, r3, r4, and r5 mutually exclusive integers and are randomly generated within the range: [1, 2, ···, NP].

Among the different vector generation strategies mentioned above, MT1 has a strong exploration capability and is suitable for solving multimodal problems. However, this strategy usually exhibits a slow convergence rate. The two-difference vector-based strategies “MT2” and “MT3” have better perturbation effects than the one-difference vector-based strategy, but they require a higher computational cost. The MT3 strategy relies on the best solution found so far, has a fast convergence rate, and performs well for single-peaked problems. However, this strategy will likely fall into a local optimum for multipeaked issues, leading to premature convergence. The “MT4” strategy is rotationally invariant. The results show that the MT4 strategy is effective for solving multiobjective optimization problems.

Crossover

To increase the diversity of the perturbation parameter vectors, a crossover approach was used. At generation G, concerning each mutation vector V = [V1, V2, ···, V], the trial vector U = [U1, U2, ···, U] will be generated according to the following crossover equationwhere CR refers to the crossover rate, which controls the fraction of parameter values obtained from the mutation vector with positive values between [0,1], rand refers to the jth sought value of the uniform random number generator, which takes values in [0,1], and jrand is a chosen integer in [1,D] randomly and it ensures that one parameter in U is different from the target vector x at least.

Selection

For the minimization problem, each target vector x, and its corresponding trial vector U, a selection procedure is performed using the fitness function. The population with the lower fitness function value will be selected as the next-generation population.

Self-Adaptive Differential Evolution Algorithm

The control parameters used by the DE algorithm, such as CR, F, and NP, with the learning strategy, determine the algorithm’s performance. For achieving superior performance on a given problem, fine-tuning individual key control parameters is essential. Additionally, but more importantly, all available learning strategies must be tested in the mutation phase. A better solution is to use the SaDE algorithm that can adjust the control parameters and learning strategies during the evolutionary process through the self-adaptive mechanism automatically.

Among the three key control parameters of the DE algorithm, CR, F, and NP, NP is adopted as a user-defined value to deal with different dimensions. F corresponds to the speed of convergence, which should be maintained throughout the evolutionary process for local searches with smaller values of F and global searches with larger values of F, to produce influential mutant individuals. For different individuals in the current population, F can select any value in the (0, 2) range of the normal distribution with mean (0.5) and standard deviation (0.3) randomly.

In general, CR affects the convergence speed and robustness of the search process, where choosing a good CR can achieve better learning results under different learning strategies. Therefore, the value of CR can be dynamically adjusted within a particular generation interval by accumulating previous learning experiences. After several generations, the CR values are modified several times, using all recorded CR values, and then the mean of the CR normal distribution is recalculated based on the vector of successful trials. When the process is repeated with this new mean and standard deviation of the normal distribution, the range of correct CR values for the current problem has been automatically learned.

This work introduces a self-adaptive mechanism based on statistical experience and the roulette wheel selection method to select the better mutation strategy from the pool of variant strategy candidates at different stages of evolution. MT1, MT2, MT3, and MT4 are four candidate learning strategies that can be selected during evolution. The probabilities of applying the learning strategies to each individual in the current population are P1, P2, P3, and P4, respectively. The initial possibilities of all four strategies are equal (0.25). For a population of size NP, a probability vector of length NP can be randomly generated and uniformly distributed in the range [0, 1].

If the value of the lth element of the vector is less than or equal to P1, the strategy MT1 will be applied to the lth individual in the current aggregate or moved to the next available strategy. After evaluating all trial vectors, the number of new trial vectors generated by different strategies into the next-generation G is NS1, NS2, NS3, and NS4. The number of discarded trial vectors is NF1, NF2, NF3, and NF4. These values are accumulated within a specific number of G that is called a learning period. Then, the probability p of strategy l is updated as followsAt the end of the learning period, it will automatically update the probabilities of applying the four learning strategies and reset the values obtained in the previous period for NS1–NS4 and NF1–NF4. This process will gradually evolve as the most suitable learning strategy for the given problem.

The performance of the LSSVM is affected to a significant degree by two key hyperparameters (eqs –8): γ and σ. Therefore, these hyperparameters are to be optimized using the SaDE algorithm, and the optimization process is shown in Figure .

Figure 3

Flow diagram for LSSVM’s parameter optimization with SaDE.

The variables are chosen as the hyperparameters γ and σ to be optimized, with the selected fitness value function as followswhere y(x) and y(x) are the true and proxy data at the test point x, respectively, and R(γ,σ) refers to the root mean square of the prediction error (RMSE), which varies with the LSSVM parameters γ and σ. For proxy models, any error or correlation metric is obtained by validation in a test data set. The complexity of a proxy model requires multiple methods to assess its accuracy. The indicator for evaluating the performance of a regression model is generally the regression accuracy, which is also known as the model accuracy. In addition to RMSE, relative error (RE) is also a fundamental factor to be considered.

Generally, the relative error better reflects the confidence level of the calculation. It is defined by the following equation

Constructing Integrated Optimization Model

The utilization of hydraulic fracturing and horizontal drilling has facilitated the recovery of oil reserves from shale and other tight formations, resulting in horizontal well design and fracture optimization being the most critical aspects of low-permeability reservoir development design. The characteristics of horizontal well placement and fracturing design are analyzed, combined with the proxy model in Section , which establishes an integrated mathematical model for optimal strategy in this section.

Most conventional horizontal well design and fracturing optimization methods focus on finding the optimal design parameters to specific metrics, such as estimated ultimate recovery (EUR), cumulative gas production over a while, or net present value (NPV).[24−30] The economic analysis recommends the decision with maximum NPV (present value of net benefits, or benefits minus costs, over time) or maximum benefit–cost ratio (present value of benefits to present value of fees). The NPV is calculated to obtain the present value of future income generated by eliminating the use of standard economic techniques. NPV is the difference between the present value of cash inflows and outflows associated with a project, adjusted for inflation to its present value, the objective function most widely used in reservoir optimization. For the optimization problem of horizontal well placement and fracturing design in shale gas reservoirs, the objective function (NPV) of the optimization is calculated as followswhere xf and Lw denote the optimization variables fracture half-length and horizontal well length, respectively; Nt and t are the total year and the time step calculated to year t, respectively; b is the discount rate; Pg denotes the natural gas price; Qgas denotes the total gas production in the t-th year; Coperate is the operating cost of natural gas in the year t; ctax is the tax constant; Cground is the cost of ground works; and Cdrill and Cfracture denote drilling costs and fracturing costs, respectively.

Drilling costs are usually defined as a linear function of horizontal well lengthwhere Cd fix denotes the fixed cost of drilling a new horizontal well and C is the drilling penetration cost per lateral length of the horizontal well. Table provides the values for the fixed cost and penetration cost of the drilling process.

Table 1

Parameter for the NPV Calculation

parameter	values	unit	symbols
discount rate	8	%/100	b
gas price	1.3	CNY/m3	Pg
operating cost	2.64	106 CNY	Coperate
tax constant	0.3277	CNY/m3	ctax
cost of ground works	5	106 CNY	Cground
fixed cost of drilling	19.923	106 CNY	Cd fix
drilling cost per lateral length	2269	CNY/m	Cd
fixed fracturing cost	8.57	106 CNY	Cfix

In horizontal well multistage fracturing, large volumes of fracturing fluid and proppant are pumped into the formation. The proppant is pumped into the fractured rock to help improve the permeability of the formation fluids once the injected fracturing fluid has released the fracturing pressure. Therefore, the fracturing cost is primarily determined by the fracturing fluid and proppant. The following fracturing costs are definedwhere cf and cp are the fracturing fluid and proppant costs per fracturing section, respectively; nf, is the fracture section number, which is determined by the horizontal section length and fracture spacing; Cfix is the fixed fracturing cost per well; rp is the proppant concentration in the fracturing fluid; and Fin, is the total volume of fracturing fluid. Typically, xf and Fin, are nonlinearly related. This is partly because as the fracture half-length increases, the filtration loss increases due to the increased contact area between the fracturing fluid and the formation during injection and also because the injection pressure decreases as the fracture half-length increases under the influence of formation frictional resistance. This ultimately leads to the linear increase of injection volume and does not achieve the fracturing result of the linear increase of fracture half-length. In this study, according to the literature,[31,32] we regressed the relationship between Fin, and xf as a function of the following, based on the statistics of actual fracturing data

It is essential to emphasize that production forecasting and optimal design should be an integrated system. With the development of computer technology, researchers utilize optimization algorithms combined with proxy models to deal with production problems such as well placement and fracturing optimization design.[33−36] The optimization results tend to be more influenced by the optimization algorithm when the accuracy of the proxy model is sufficient.

Guyaguler et al.[37] proposed a hybrid genetic algorithm (HGA), which reduces the computational burden of performing many simulations by combining GA and general kriging algorithms. They obtained results consistent with observations from the Pompano field in the Gulf of Mexico and used the HGA to study the PUNQ-S3 reservoir. Ma et al.[38] applied gradient-based finite difference (FD), discrete simultaneous regressive stochastic approximation (DSPSA), and genetic algorithm (GA) to resolve the hydraulic fracture arrangement problem. Plaksina et al.[39] offer the algorithm that gives quantitative and qualitative measures of “goodness” of the optimal production plans, which can handle objectives effectively and produce the Pareto optimal solutions without requiring the user to assign weights to each aim inside an aggregate function. Li et al.[30] proposed a dynamic simplex interpolation-based alternating subspace (DSIAS) search method for the mixed-integer optimization problem of shale gas development projects and also presented a case study on the development of the Barnett Shale gas field and validated the optimization effect.

Results and Discussion

In this section, we apply the SADE-LSSVM model to actual blocks (the lower Silurian Longmaxi Formation from the southern Sichuan Basin). First, a production proxy model is built based on the adaptive response surface experimental design with the historical fitting of regional production data. The applicability of the proxy model was also verified, with the proxy model remaining efficient and accurate when the data set was updated.

Model Building and Experiment Designing

The central target formation in the study area is the Longmaxi shale gas formation, which has a reservoir thickness of 30 m, a reservoir depth of about 2500–3000 m, average porosity of 4%, an average permeability of 3 × 10–5 mD, low-pressure coefficient, and fault and microfracture developed. The reservoir is assumed to be homogeneous, with uniform fracture spacing and porosity, and permeability independent of stress. Only gas flows in the reservoir, and the effects of gas desorption are considered. Detailed reservoir information for this section is shown in Table .

Table 2

Parameters Used in Reservoir Simulation

parameter	value	unit
porosity	4	%
matrix permeability	3.6 × 10–5	mD
gas content saturation	65	%
natural fracture spacing	1	m
hydraulic fracture half-length	105	m
hydraulic fracture height	30	m
fracture conductivity	1.2	mD·m

Self-adaptive response surface experimental design (SaRSE) was applied based on six uncertain parameters with reasonable ranges of parameter values:

porosity (⌀)

matrix permeability (p)

fracture half-length (xf)

horizontal section length (lw)

gas content saturation (Sg)

SaRSE experimental data are determined by the response surface regression accuracy, with the number of experiments increasing until the specified accuracy is reached. Experimental data thus obtained can better reflect the actual reservoir characteristics and seepage patterns. Numerical simulation software was used to calculate the 15-year gas production (Cg). Some of the experimental design results are shown in Supporting Information Table S1.

Proxy Model Accuracy

Based on SaRSE, the response surface regression formula will be obtained. The results of adaptive response surface regression are shown in Figure . The response surface formula predicts results with low accuracy. Therefore, a new method is warranted to build the agent model.

Figure 4

Regression results of adaptive response surface.

Cross-validation is introduced to improve the model’s reliability. The total data set was divided into three data sets: training set, validation set, and test set. More specifically, first, 25% of the entire data set was randomly selected as the test set and excluded from the training and validation process. The remaining data set was then randomly divided into 10 parts, and a fivefold cross-validation method was used to use 60% of the data as the training set and 15% of the data to validate the model for each cycle.

Figure depicts the Cg for the training, validation, test, and total data sets. The results show that the error difference between the training and test data sets is slight, indicating that the training process is reliable (no overfitting).

Figure 5

SaDE-LSSVM comparison of predicted and target Cg for (a) training set, (b) validation set, (c) testing set, and (d) total set.

The predictions are in good agreement with the target values, and the accuracy is within the acceptable range. We compared the SADE-LSSVM model with the PSO-LSSVM, GA-LSSVM, and GWO-LSSVM models, showing the high prediction accuracy of the SADE optimized LSSVM, as shown in Figures –8. Table shows a more comprehensive set of error metrics.

Figure 6

PSO-LSSVM comparison of predicted and target Cg for (a) training set, (b) validation set, (c) testing set, and (d) total set.

Figure 8

GWO-LSSVM comparison of predicted and target Cg for (a) training set, (b) validation set, (c) testing set, and (d) total set.

Table 3

Comparison of LSSVM Model Errors (RMSE and RE)

proxy model		SaDE-LSSVM	PSO-LSSVM	GA-LSSVM	GWO-LSSVM
training set	RMSE	0.002264	0.005254	0.007137	0.005262
	RE	0.001350	0.002044	0.002488	0.002047
validation set	RMSE	0.002610	0.005651	0.007552	0.005659
	RE	0.001368	0.002677	0.003418	0.002680
testing set	RMSE	0.002955	0.007006	0.009437	0.007016
	RE	0.001429	0.002632	0.003381	0.002636
total set	RMSE	0.002503	0.005780	0.007808	0.005788
	RE	0.001372	0.002290	0.002856	0.002293

GA-LSSVM comparison of predicted and target Cg for (a) training set, (b) validation set, (c) testing set, and (d) total set.

Sensitivity Analysis and Optimization Results

A first sensitivity analysis was carried out according to the experimental design in Supporting Information Table 1. The F-value of analysis of variance is a fine indicator of the degree of influence of different parameters. Figure shows that all investigated parameters have different degrees of influence on Cg and NPV. The results showed that Cg and NPV exhibited the same behavior in response to parameter variations. Cg and NPV are sensitive to geological properties (⌀ and Sg) to a high degree. The effect of k on Cg is relatively small. Shale gas is mainly recovered by hydraulic fracturing with high inflow capacity, where changing the matrix permeability does not bring significant changes in production, which is the reason why hydraulic fracturing is used for shale gas development.

Figure 9

Degree of influence of geological and engineering parameters on (a) Cg and (b) NPV.

A second sensitivity study was then performed for each influencing parameter to determine the effect of individual parameters on Cg and NPV based on the experimental data in Supporting Information Table 1. The effect patterns of parameters on Cg and NPV are shown in Figure . Obviously, the increase of ⌀ and Sg can directly increase Cg with no cost constraint, which eventually makes the NPV show the same trend as Cg. However, an xf growth leads to a higher Cfracture while increasing Cg, which results in a different trend of Cg and NPV. For lw, both Cg and Cdrill are close to linear to lw, which leads to Cdrill not being able to constrain NPV, which makes Cg and NPV trend the same while the optimal lw is taken to the extreme value on the other hand.

Figure 10

Effects of (a) ⌀, (c) Sg, (e) xf, and (g) lw on Cg, and effects of (b) ⌀, (d) Sg, (f) xf, and (h) lw on NPV.

In this study, the SaDE-LSSVM proxy model with relatively high accuracy is applied to the single-well NPV optimization with the optimization variables xf and lw. The final optimization results are xf = 155 m and lw = 2000 m.

Practical Application of Proxy Model

The numerical simulation itself is a simplified mathematical model, so the LSSVM fitting is excellent. The reliability of the LSSVM method can be tested with actual field data. Therefore, to further validate the applicability of the proxy models, we developed the EUR proxy model and the absolute open flow (AOF) proxy model for the study area based on actual data. We select 16 real production wells data as training samples.

The characteristic variables are formation pressure, clump number, fracturing fluid volume, proppant dose, ⌀, k, Sg, rejection rate, fracturing pressure, and pumping stop pressure. Five wells were also selected as the test set to verify the prediction accuracy of the proxy model. Figure shows the accuracy of the SADE-LSSVM model and PSO-LSSVM model, respectively. A comprehensive set of error evaluations is shown in Table . The proxy model errors are within acceptable limits, indicating the reliability of using actual production data to build the proxy models.

Figure 11

Predicted results of QAOF and EUR under SaDE-LSSVM and PSO-LSSVM, respectively.

Table 4

Comparison of Errors in Predicting QAOF and EUR by LSSVM Models

proxy model			SaDE-LSSVM	PSO-LSSVM
QAOF	training set	RMSE	0.039221	0.022444
		RE	0.828875	0.425939
	testing set	RMSE	0.030152	0.019894
		RE	0.633160	0.812785
EUR	training set	RMSE	0.013228	0.054516
		RE	0.032447	0.042602
	testing set	RMSE	0.098707	0.100341
		RE	0.079433	0.100720

Conclusions

In this study, a proxy model of shale gas production, SADE-LSSVM, is developed using a swarm optimization algorithm and machine learning. Its errors are compared and evaluated, showing a high accuracy of the proxy model. The proxy model built by the swarm optimization algorithm combined with the machine learning approach still has strong applicability when the data set is updated. We obtained the following recognition points:

(1) The developed adaptive proxy modeling method has high accuracy, due to the fact that a sufficient number of data are selected with the self-adaptive response surface design in the developed method, which improves the accuracy of the proxy model.

(2) The developed proxy model can simulate the reservoir simulator response with acceptable accuracy, thus replacing numerical simulation in the production optimization process to satisfy the optimization process under multiple iterations and multiple control parameter variations, which is particularly important for production optimization applications.

(3) The proxy models can be easily updated when any new data sets are available. These data sets include geological data, drilling data, fracturing data, or other custom design-of-experiment data. This data-driven approach has great potential for application by quickly implementing proxy models on other shale reservoir data sets.