A rapid simulation method applicable to complex reservoirs and a real-time optimization method for layered injection and production

Abstract

The present invention discloses a rapid simulation and real-time optimization method for stratified injection and production applicable to complex reservoirs, belonging to the field of numerical simulation of unconventional reservoirs, and comprising the following steps: Step 1, constructing a GMM-XGB hybrid model based on the Gaussian mixture model and the XGBoost algorithm; Step 2, constructing a stratified connectivity network model based on a physical model; Step 3, performing automatic history matching by using the ES-MDA algorithm; Step 4, establishing a mathematical model of a real-time optimization objective function for stratified injection and production, and performing real-time optimization of stratified injection and production based on the SHADE algorithm. The present invention can improve the recovery rate of complex reservoirs, extend the stable production period of oilfields, and achieve efficient utilization of resources and maximization of economic benefits.

Description

Technical Field

[0001] The present invention belongs to the field of numerical simulation of unconventional reservoirs, and particularly relates to a rapid simulation method suitable for complex reservoirs and a real-time optimization method for layered injection and production. Background Art

[0002] Reservoirs containing fractures and faults have obvious geological complexities, which pose great challenges to their effective development and simulation calculations. Fractures will exacerbate the planar and vertical heterogeneity of the formation, develop water injection preferential channels on the plane, and increase the difference in layer utilization in the vertical direction. Therefore, it is a difficult problem to achieve balanced and efficient development of well groups in fractured reservoirs. At the same time, the presence of faults will lead to the partition of the reservoir or the loss of local connectivity, thereby affecting the oil recovery efficiency and water flooding effect. In three-dimensional space, due to the discontinuity and irregularity of fractures and faults, traditional reservoir simulation methods require expensive costs to depict fractures and faults, and then perform simulation calculations more accurately. Summary of the Invention

[0003] In order to solve the above problems, the present invention proposes a rapid simulation method suitable for complex reservoirs and a real-time optimization method for layered injection and production. A GMM-XGB hybrid model is constructed by integrating the Gaussian mixture model and the XGBoost algorithm to predict the injection and production liquid volumes of small layers. And the characterization of faults and fractures is embedded when constructing the flow network model, and a layered connectivity network model suitable for complex reservoirs is successfully constructed. Then, the ES-MDA algorithm considering formation constraints is used to automatically history match and invert the inter-well connectivity and the dynamic effects of fractures and faults. On this basis, an adaptive mechanism based on successful history and the SHADE algorithm with an archive mechanism are introduced, and the economic net present value is used as the optimization target, considering various engineering constraints, to achieve real-time optimization of layered injection and production.

[0004] The technical solution of the present invention is as follows:

[0005] A rapid simulation method suitable for complex reservoirs and a real-time optimization method for layered injection and production, comprising the following steps:

[0006] Step 1: Construct a GMM-XGB hybrid model based on the Gaussian mixture model and the XGBoost algorithm;

[0007] Step 2: Construct a layered connectivity network model based on a physical model;

[0008] Step 3: Use the ES-MDA algorithm for automatic history matching;

[0009] Step 4: Establish a mathematical model of the real-time optimization objective function for layered injection and production, and perform real-time optimization of layered injection and production based on the SHADE algorithm.

[0010] Further, in step 1, the GMM-XGB hybrid model uses the Gaussian mixture model to perform clustering analysis on the original sample set, obtains the probability distribution of the injection-production fluid volume of the small layer, integrates the probability distribution into the original sample set as the input of the XGBoost algorithm, and predicts and outputs the injection-production fluid volume of the small layer. The specific working process of the GMM-XGB hybrid model is as follows:

[0011] Step 1.1: Collect the vertical splitting data of the injection-production fluid volume and construct the original sample set;

[0012] The vertical splitting data of the injection-production fluid volume includes porosity, permeability, oil layer thickness, initial water saturation, formation coefficient, number of perforations, perforation thickness, permeability range, porosity range, bottom hole flowing pressure, single-well injection-production fluid volume, and water cut;

[0013] Step 1.2: Perform data preprocessing on the original sample set. The specific process is as follows:

[0014] First, perform data cleaning. For the missing values of the random missing type, according to the geodetic coordinates of the well points, use the natural neighbor interpolation method in the matlab software for interpolation. For the outliers, select deletion or transformation processing according to the actual situation. The transformation processing includes smoothing or scaling. For the categorical data, first, it is necessary to judge the number of base classes of the categorical data. If the number of base classes is small, use the one-hot encoding method for processing. If the number of base classes is large, use the label encoding or target encoding method for processing;

[0015] Then, perform normalization processing on the cleaned data using the Min-Max normalization method;

[0016] Finally, perform correlation analysis on the normalized sample set. For the linear relationship in the data, use the Pearson correlation coefficient for correlation analysis. For the sample set with strong non-linear relationship and non-normal distribution, use the Spearman rank correlation coefficient for correlation analysis;

[0017] Step 1.3: Use the Gaussian mixture model to perform clustering analysis on the preprocessed original sample set to obtain the probability distribution of the injection-production fluid volume of the small layer;

[0018] Step 1.4: Integrate the probability distribution of the injection-production fluid volume of the small layer with the preprocessed original sample set to obtain the integrated sample set, and divide the integrated sample set into a training set and a validation set;

[0019] Step 1.5: Input the integrated sample set into the XGBoost algorithm and use the Bayesian optimization method for training until the requirements are met, and finally output the injection-production fluid volume of all wells, all time steps, and all small layers.

[0020] Further, the specific process of step 1.3 is as follows:

[0021] Step 1.3.1: Based on the preprocessed original sample set, construct a Gaussian mixture model; the Gaussian mixture model assumes that the injection-production fluid volume data points come from multiple different Gaussian distributions, and the probability density function of the Gaussian distribution is:

[0022] ;

[0023] where is the probability density function of the Gaussian distribution; is the injection-production fluid volume data point; is the mean vector; is the covariance matrix; is the data variance; is the data dimension; is the exponential function with base e;

[0024] Step 1.3.2: Use the Expectation-Maximization algorithm to iteratively solve the established Gaussian mixture model until convergence. The specific process is as follows:

[0025] First, initialize the parameters of each Gaussian distribution in the Gaussian mixture model;

[0026] Then, calculate the responsiveness of all injection-production fluid volume data points. The responsiveness calculation formula is as follows:

[0027] ;

[0028] where is the responsiveness; is a binary indicator variable indicating whether the th injection-production fluid volume data point is generated by the th Gaussian component; is the prior probability of the th Gaussian component; is the mean vector of the th Gaussian component; is the covariance matrix of the th Gaussian component; is the index variable of the Gaussian component, used to traverse all Gaussian components in the mixture Gaussian model for summation operations; is the total number of Gaussian components in the mixture Gaussian model; is the prior probability of the th Gaussian component; is the mean vector of the th Gaussian component; is the covariance matrix of the th Gaussian component;

[0029] Update the parameters of each Gaussian distribution according to the calculated responsivity, and the update formula is as follows:

[0030] ;

[0031] ;

[0032] ;

[0033] where, , , are the mean vector, covariance matrix, and prior probability of the -th iteration of the -th Gaussian component respectively; is the responsivity at the -th iteration; is the total number of injection-production fluid volume data points; is the transpose symbol;

[0034] Calculate the likelihood function of the model, and the formula is as follows:

[0035] ;

[0036] where, is the likelihood function;

[0037] Finally, continuously iterate and update the model parameters until the parameters of the Gaussian mixture model converge. At this time, the training ends, and the final responsivity is output. The responsivity is the probability distribution of the injection-production fluid volume of the sub-layer.

[0038] Furthermore, in step 1.5, Bayesian optimization is a method for automatic hyperparameter tuning. By establishing a surrogate model, the hyperparameters are continuously adjusted according to the historical performance of the model, so as to find the optimal hyperparameter configuration; the specific process is as follows:

[0039] Step 1.5.1, Initialize the Bayesian optimization algorithm framework, and initialize the evaluation of the objective function and set the hyperparameter search space of the GMM-XGB hybrid model in the framework; the objective function uses the mean absolute error, and the formula is as follows:

[0040] ;

[0041] where, is the mean absolute error; is the true value of the -th sample; is the predicted value of the -th sample; is the total number of samples;

[0042] The hyperparameters include the learning rate, the maximum depth of the tree structure, the subsample ratio, and the number of tree structures; the search space for setting the learning rate is [0.01, 0.2], the search space for the maximum depth of the tree structure is [5, 10], the search space for the subsample ratio is [0.5, 1], and the search space for the number of tree structures is [100, 500];

[0043] Step 1.5.2, iterative training to find the optimal hyperparameters of the model within the search space; first, sample collection is carried out, that is, a set of initial points is selected as the initial hyperparameters in the hyperparameter combination space of the model to be optimized, and then the Bayesian optimization algorithm evaluates the next set of hyperparameters to be selected based on the existing hyperparameters and the objective function values; finally, the hyperparameter combination with the largest objective function is selected as the optimal hyperparameter combination;

[0044] Step 1.5.3, train the GMM-XGB hybrid model using the selected optimal hyperparameter combination on the training set, and calculate the mean absolute error between the predicted value and the true value of the GMM-XGB hybrid model on the validation set to obtain the objective function value;

[0045] Step 1.5.4, after obtaining the new hyperparameter-objective function pair, update the surrogate model, and the surrogate model will continuously adjust the estimation of the objective function distribution in the hyperparameter space according to the new data;

[0046] Step 1.5.5, determine whether the pre-set termination condition is satisfied; if the termination condition is not satisfied, return to Step 1.5.2 to perform sample collection again and continue the iteration; when the termination condition is satisfied, select a set of hyperparameters that optimize the objective function value from all the evaluated hyperparameter combinations, and use this set of optimal hyperparameters to retrain the GMM-XGB hybrid model on the training set until the pre-set number of training times is reached; after the training is completed, the final GMM-XGB hybrid model for prediction is obtained.

[0047] Furthermore, the specific process of Step 2 is as follows:

[0048] Step 2.1, divide the reservoir into several well points and connected pipelines to construct an inter-well connection network; in the inter-well connection network, each well point represents an injection well or a production well on a small layer, and the well points are connected to each other through connected pipelines; the connected pipelines are characterized by two parameters, the inter-well conductivity and the connected volume;

[0049] The conductivity calculation formula is as follows:

[0050] ;

[0051] Among them, is the conductivity between well and well ; is well Seepage area between and well and well ; Average value of permeability between and well and well ;

[0052] The calculation formula of the connected volume is as follows:

[0053] ;

[0054] Wherein, Connected volume between and well and well ; Average value of effective thickness between and well

[0055] Step 2.2: Discretize the fracture into a group of uniformly arranged high-permeability nodes; discretize the fault into two groups of uniformly arranged nodes, located on both sides of the fault respectively, and characterize the fluid flow capacity near the fault through the heterogeneous conductivity; embed the two types of nodes of the fracture and the fault into the inter-well connection network to construct a dynamic flow channel.

[0056] Step 2.3: Based on the hierarchical inter-well connection relationship, set multiple connection units in each layer and perform one-dimensional grid discretization, and then map them into a standardized model file according to the connection relationship of the grids to construct a physically driven hierarchical connection network model.

[0057] Furthermore, the specific process of the said Step 2.3 is as follows:

[0058] First, divide the reservoir into multiple small layers according to the geological characteristics of the reservoir, and each small layer corresponds to an independent physical model; in each layer, set multiple connection units based on geological data and well test data; for the flow characteristics in each layer, refine them by means of one-dimensional grid discretization, and each grid unit is equivalent to a part of the connection unit.

[0059] Convert the flow process of the reservoir into a discretized numerical calculation problem; the basic idea of grid discretization is to convert the continuous fluid flow process in each connection unit into a discrete system composed of multiple small units, and the permeability, volume and flow rate of each small unit are solved by numerical methods; the connection relationship between discrete units is mapped through a standardized model file to generate a hierarchical connection network model.

[0060] Furthermore, the specific process of the said Step 3 is as follows:

[0061] Step 3.1. Construct the history matching objective function as follows:

[0062] ;

[0063] where, is the history matching objective function value; is the model fitting parameter vector; is the model simulation calculation result; is the actual observed data; is the covariance matrix of the observation error;

[0064] Step 3.2. Select the conductivity between discrete grids, the volume of discrete grids, the heterogeneous conductivity and the number of characteristic points of fractures, the heterogeneous conductivity and the number of characteristic points of faults, and the relative permeability parameters as control variables to form the model fitting parameter vector;

[0065] ;

[0066] where, is the conductivity between grid and grid at the initial time; is the volume of grid at the initial time; is the heterogeneous conductivity of the characteristic points and characteristic point of the fracture on layer at the initial time; is the number of characteristic points of the fracture on layer at the initial time; is the heterogeneous conductivity of the characteristic points and characteristic point of the fault on layer at the initial time; is the number of characteristic points of the fault on layer at the initial time; , , , are different relative permeability parameters;

[0067] During the reservoir history matching process, while minimizing the objective function, physical constraint conditions are imposed, specifically including:

[0068] ;

[0069] where, is the conductivity of grid and grid ; is the volume of grid ; is the number of grids; is the total volume of the reservoir;

[0070] Step 3.3: Use the ES-MDA algorithm, combined with historical production data and physical constraints, to achieve automatic inversion and real-time update of fracture and fault parameters. The specific process is as follows:

[0071] Firstly, the initial heterogeneous conductivity distribution of fractures and faults is set in the reservoir numerical model, and the number of characteristic points is defined. Then, in the history fitting process, the actual production data of oil wells (daily oil production per well) is used as observation information to perform multiple data assimilation on the model. Through the iterative update mechanism of the ES-MDA algorithm, the matrix parameters as well as the parameters of fractures and faults are continuously adjusted.

[0072] Furthermore, the specific process of step 4 is as follows:

[0073] Step 4.1: Select the economic net present value as the objective function of real-time optimization of injection and production; the following is the mathematical model of the optimization objective function:

[0074] ;

[0075] in, is the economic net present value; is the controlled variable, representing the regulation scheme of stratified water injection volume and stratified liquid production volume during the optimization process; for oil prices; For the Time step well Oil production; Cost of sewage treatment; For the Time step well The water production; is the water injection cost; For the Time step well The amount of water injected; is the total number of time steps; is the number of producing wells; is the number of water injection wells;

[0076] Step 4.2: Set engineering constraints, including injection-production ratio constraint, bottom hole pressure constraint, maximum water cut constraint, and single well maximum liquid volume constraint. The specific constraints are as follows:

[0077] ;

[0078] In the formula, is the minimum injection-production ratio; is the maximum injection-production ratio; For the Time step well Liquid production or injection volume; For the th time step of the well Minimum bottom hole flowing pressure; For the th time step of the well Maximum bottom hole flowing pressure; For the th time step of the well Bottom hole flowing pressure; For the th time step of the well Water cut; For the th time step of the well Ultimate water cut; For the th time step of the well Liquid production or injection volume; For the th time step of the well Ultimate liquid production or ultimate injection volume;

[0079] Step 4.3. Introduce an adaptive mechanism based on successful history and the SHADE algorithm with an archive mechanism for reservoir injection-production optimization.

[0080] Furthermore, in the said Step 4.3, the specific process of the SHADE algorithm for reservoir injection-production optimization is as follows:

[0081] Step 4.3.1. Adjust the injection-production flow rate to generate a new injection-production plan;

[0082] Step 4.3.2. Combine different injection-production plans and calculate the economic net present value of each injection-production plan combination;

[0083] Step 4.3.3. Select the injection-production plan combination that meets the engineering constraint conditions and has the maximum economic net present value as the individual entering the next generation. The individual will perform Step 4.3.1 and Step 4.3.2 again according to the adaptive mechanism until the number of iterations is reached; the final injection-production plan combination after the iteration ends is the optimal injection-production plan.

[0084] Advantageous technical effects brought by the present invention: The present invention uses a Gaussian mixture model for clustering analysis. Based on the clustering analysis results, the XGBoost algorithm is used to train the data with fused clustering features to obtain a GMM-XGB hybrid model. Compared with the conventional XGBoost method, the GMM-XGB hybrid model of the present invention can better discover the split data distribution pattern, initialize the XGBoost model, and thus reduce the overfitting problem. When establishing a flow network model, faults and fractures are embedded, which can characterize the blocking effect of faults and the conduction advantage of fractures, thereby improving the simulation accuracy of the fluid flow path in complex reservoirs and optimizing the injection-production regulation effect. The SHADE algorithm has more advantages in finding the global optimal solution compared with the traditional DE algorithm. The model established by the present invention can improve the recovery rate of complex reservoirs, extend the stable production period of oilfields, and achieve the efficient utilization of resources and the maximization of economic benefits. BRIEF DESCRIPTION OF THE DRAWINGS

[0085] Figure 1 It is a flowchart of the longitudinal production split and real-time injection-production optimization method applicable to complex reservoirs of the present invention.

[0086] Figure 2 It is a flowchart of predicting the injection-production fluid volume of sub-layers of the present invention.

[0087] Figure 3 It is a schematic diagram of the layered connection network model in the embodiment of the present invention.

[0088] Figure 4 It is a schematic diagram of the historical fitting result when the first iteration of sub-layer 1 of production well 1 in the embodiment of the present invention.

[0089] Figure 5 It is a schematic diagram of the historical fitting result when the second iteration of sub-layer 1 of production well 1 in the embodiment of the present invention.

[0090] Figure 6 It is a schematic diagram of the historical fitting result when the third iteration of sub-layer 1 of production well 1 in the embodiment of the present invention.

[0091] Figure 7 It is a comparison chart of water cut curves after historical fitting in the block in the embodiment of the present invention.

[0092] Figure 8 It is a comparison chart of daily oil production curves after historical fitting of sub-layer 1 in the embodiment of the present invention.

[0093] Figure 9 It is an injection regime regulation diagram of sub-layer 1 of injection well 5 in the embodiment of the present invention.

[0094] Figure 10 It is an injection regime regulation diagram of sub-layer 2 of injection well 5 in the embodiment of the present invention.

[0095] Figure 11 This is the comparison chart of cumulative oil production before and after optimization in the embodiment of the present invention.

[0096] Figure 12 This is the comparison chart of water cut before and after optimization in the embodiment of the present invention. Detailed implementation manners

[0097] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific implementation manners:

[0098] As Figure 1 shown, a rapid simulation and real-time optimization method for layered injection-production applicable to complex reservoirs includes the following steps:

[0099] Step 1: Construct a GMM-XGB hybrid model based on the Gaussian mixture model and the XGBoost algorithm; the GMM-XGB hybrid model uses the Gaussian mixture model to perform clustering analysis on the original sample set to obtain the probability distribution of injection-production fluid volume, and integrates the probability distribution into the original sample set as the input of the XGBoost algorithm to predict and output the injection-production fluid volume of sub-layers. The specific working process of the GMM-XGB hybrid model is as follows:

[0100] Step 1.1: Collect longitudinal splitting data of injection-production fluid volume and construct an original sample set;

[0101] The longitudinal splitting data of injection-production fluid volume includes various factor data affecting the fluid volume of sub-layers, such as porosity, permeability, oil layer thickness, initial water saturation, formation coefficient, number of perforations, perforation thickness, permeability range, porosity range, bottom hole flowing pressure, single well injection-production fluid volume, and water cut.

[0102] Step 1.2: Perform data preprocessing on the original sample set, including data cleaning and correlation analysis;

[0103] When constructing the sample set, it is necessary to ensure the integrity and consistency of the data. If the data source of the sample set is relatively large, there may be some noise and inconsistencies, so data cleaning must be carried out to provide a reliable basis for subsequent modeling and analysis. For missing values, first judge the type of missing values, and select a suitable method for handling missing values according to the type of missing values. For outliers, first detect outliers and outliers, and select a processing method of deletion or transformation according to the actual situation. For the categorical data in the sample set, first judge the base class level of the categorical data to select the coding method. Perform normalization processing on the cleaned data to adjust the data scale and better meet the requirements of machine learning algorithms. In order to find variables strongly correlated with the sub-layer injection-production fluid volume data, it is necessary to perform correlation analysis on the normalized sample set. For the sample set with strong non-linear relationship and non-normal distribution, the Spearman rank correlation coefficient is used for correlation analysis. The specific process is as follows:

[0104] Step 1.2.1: Missing values are mainly divided into three types: Missing Completely at Random (MCAR), Missing at Random (MAR), and Missing Not at Random (MNAR). After analysis, the missing values in reservoir data are mainly due to missing during on-site testing or recording. Such missing often occurs for several parameters simultaneously, so it is mainly MAR. According to the geodetic coordinates of the well points, the natural neighbor interpolation ("natural") method in Matlab software is used for interpolation. This method is based on Delaunay triangulation and achieves an effective balance between linear and cubic interpolation.

[0105] Step 1.2.2: Outliers and extreme points are important factors affecting data quality, so effective detection and processing are required. First, use box plots to detect outliers in the sample set and identify values that deviate from the normal distribution. Specifically, calculate the first quartile, median, and third quartile. After sorting the data from smallest to largest, the data at the 25% position is the first quartile, the data at the 50% position is the median, and the data at the 75% position is the third quartile. The box range is the range of most of the data, and the formula for the box is:

[0106] (1);

[0107] where, is the box; is the third quartile; is the first quartile;

[0108] Outliers are usually data points outside the upper and lower edges. The specific formulas for the upper and lower edges are:

[0109] (2);

[0110] (3);

[0111] where, is the upper edge; is the lower edge;

[0112] For the detected outliers, deletion or transformation can be selected according to the actual situation. For example, for significant outliers, direct deletion may be required; for some marginal data, data smoothing or scaling can be selected to make it conform to the normal data distribution. The key to this step lies in balancing the rigor of outlier processing and data fidelity.

[0113] Step 1.2.3: Some features in the sample set may be categorical data, such as geological type, injection-production type, etc. First, it is necessary to determine the number of base classes of the categorical data. If the number of base classes is small, the One-Hot Encoding method can be used for processing; while if the number of base classes is large, it is more efficient to use the Label Encoding or Target Encoding method. Select an appropriate encoding method according to different data types and task requirements. The appropriate encoding method helps to improve the training effect and prediction ability of the model.

[0114] Step 1.2.4: Inconsistent scales of data may affect the performance of machine learning algorithms, especially when using distance-based algorithms. To solve this problem, data normalization is a very important step. Normalization reduces the dimensional difference between different features by mapping the data to a unified scale range (such as [0,1] or [-1,1]), enabling machine learning algorithms to converge better. During the normalization process, the commonly used method is Min-Max normalization. Min-Max normalization can accelerate model training, improve the stability and accuracy of the model. The specific formula is:

[0115] (4);

[0116] where, is the data in the th column of the normalized data set; all columns after normalization form the data set , which is a matrix containing all features and labels; is the data in the th column of the data set before normalization; is the minimum value in the th column of data; is the maximum value in the th column of data.

[0117] Step 1.2.5: In order to find variables strongly correlated with the small-layer injection-production fluid volume data, it is necessary to perform a correlation analysis on the normalized sample set. First, calculate the correlation between each feature to understand which variables have a significant relationship with the injection-production fluid volume. For the linear relationship in the data, the Pearson correlation coefficient can be used for analysis. However, for some data with strong non-linear relationships and non-normal distributions, the Pearson correlation coefficient may not effectively reflect the relationship between variables. Therefore, using the Spearman rank correlation coefficient for analysis is a better choice. The Spearman rank correlation coefficient can measure the monotonic relationship between variables, not just the linear relationship, thus capturing the dependence relationship between variables more comprehensively. The specific calculation formula is as follows:

[0118] (5);

[0119] Among them, is the Spearman rank correlation coefficient; is the difference in the rank values between the th sample and the injection-production fluid volume data pair of the small layer; is the total number of samples.

[0120] Step 1.3: Use the Gaussian mixture model to perform clustering analysis on the preprocessed original sample set to obtain the probability distribution of the injection-production fluid volume of the small layer. The specific process is as follows:

[0121] Step 1.3.1: Based on the preprocessed original sample set, construct a Gaussian mixture model (Gaussian Mixture Models, GMM). The Gaussian mixture model assumes that the injection-production fluid volume data points come from multiple different Gaussian distributions.

[0122] The probability density function of the Gaussian distribution is:

[0123] (6);

[0124] Among them, is the probability density function of the Gaussian distribution; is the injection-production fluid volume data point; is the mean vector; is the covariance matrix; is the data variance; is the data dimension; is the exponential function with base e;

[0125] Step 1.3.2: Use the Expectation-Maximization (EM) algorithm to iteratively solve the established Gaussian mixture model until convergence, aiming to maximize the log-likelihood function of the observed data. The specific process of the EM algorithm is as follows:

[0126] First, initialize the parameters of each Gaussian distribution in the Gaussian mixture model;

[0127] Then, calculate the responsiveness of all injection-production fluid volume data points. The responsiveness calculation formula is as follows:

[0128] (7);

[0129] Among them, is the responsiveness; is a binary indicator variable, indicating whether the th injection-production fluid volume data point is from the Generation of Gaussian components; is the prior probability of the -th Gaussian component, also known as the mixing weight; is the mean vector of the -th Gaussian component; is the covariance matrix of the -th Gaussian component; is the index variable of the Gaussian component, used to traverse all Gaussian components in the Gaussian mixture model for summation operations; is the total number of Gaussian components in the Gaussian mixture model; is the prior probability of the -th Gaussian component; is the mean vector of the -th Gaussian component; is the covariance matrix of the -th Gaussian component;

[0130] Update the parameters of each Gaussian distribution according to the calculated responsivity. The update formulas are as follows:

[0131] (8);

[0132] (9);

[0133] (10);

[0134] where , , are the mean vector, covariance matrix, and prior probability of the -th Gaussian component at the -th iteration respectively; is the responsivity at the -th iteration; is the total number of injection-production fluid volume data points; is the transpose symbol;

[0135] Calculate the likelihood function of the model according to the calculated mixing weights. The likelihood function represents the goodness of fit of the model to the data. Use the likelihood function to determine whether the parameters of the model have converged, that is, if the value of the likelihood function between two iterations is less than 10 -5 , it is determined that the parameters of the model have converged. The formula for the likelihood function is as follows:

[0136] (11);

[0137] where is the likelihood function;

[0138] Finally, continuously iterate and update the model parameters until the parameters of the Gaussian mixture model converge. At this point, the training ends, and the final responsivity is output. The responsivity is the probability distribution of the fluid injection and production volume of the small layer.

[0139] When using this method to process small-sample problems such as splitting, GMM clustering has significant advantages compared to directly using the XGBoost algorithm for training: From the perspective of data understanding and feature mining, the feature information of small-sample data is limited, and it may be difficult to fully capture the data patterns directly using the XGBoost algorithm. GMM can obtain the probability distribution and different clusters of the split data. Each cluster represents a subset of data with similar features, which helps to deeply understand the internal characteristics of the data. Through clustering, hidden features that are not easily detected can be mined, providing richer and more discriminative feature information for the XGBoost algorithm, enabling the model to more accurately learn the key patterns in the data; In terms of improving the generalization ability of the model, directly using the XGBoost algorithm for split data is prone to overfitting, that is, the model overfits the training data and performs poorly when facing new data. After GMM clustering obtains different clusters and the probability distribution of the split data, the data distribution within each cluster is relatively more concentrated and similar. Training the XGBoost algorithm separately within each cluster allows the model to learn and optimize more carefully according to the characteristics of different clusters. This way of cluster-based training enables the model to better adapt to the diversity of data, reduce the dependence on specific samples, prevent overfitting, and thus enhance the generalization ability of the model on new data, making it perform more stably when facing new data.

[0140] Step 1.4: The probability distribution of the fluid injection and production volume of the small layer output contains the potential distribution pattern of the fluid volume, while the constructed original sample set contains more features related to the fluid volume. Therefore, the probability distribution of the fluid injection and production volume of the small layer output by the Gaussian mixture model in the previous step is fused with the preprocessed original sample set to obtain a fused sample set, and the fused sample set is divided into a training set and a validation set. The purpose of fusing these two is to combine the feature information of the probability distribution with the feature information in the actual samples, enhancing the model's prediction ability for the fluid injection and production volume distribution. The fusion process mainly consists of using the probability distribution output by the Gaussian mixture model as an additional input feature and combining it with other feature data samples to form a new training data set.

[0141] To reduce overfitting, the present invention uses the method of ten-fold cross-validation for partitioning. Specifically, ten-fold cross-validation divides the fused sample set into 10 subsets. In each iteration, 9 subsets are used for training, and the remaining 1 subset is used as the validation set. This process is repeated 10 times, and each subset will be used as the validation set once. Finally, the GMM-XGB hybrid model will be evaluated on all validation sets to ensure the balanced performance of the model on different data subsets, and thus a more robust GMM-XGB hybrid model is obtained.

[0142] Step 1.5: Input the fused sample set into the extreme gradient boosting (XGBoost) algorithm for training. To improve the performance of the XGBoost algorithm, the Bayesian optimization method is specifically used for training until the requirements are met, and finally the injection and production fluid volumes of all wells at all time steps and all sub-layers are output. By training different decision trees multiple times and combining them with weights, the model performance is gradually improved. The specific process is as follows:

[0143] Step 1.5.1: Bayesian optimization is a method for automated hyperparameter tuning. By establishing a surrogate model and continuously adjusting the hyperparameters according to the historical performance of the model, the optimal hyperparameter configuration can be found. The key to Bayesian optimization is to construct a surrogate model (usually a Gaussian process, GP) to approximate the objective function. In the Bayesian optimization algorithm framework, the evaluation of the objective function and the setting of the hyperparameter search space of the GMM-XGB hybrid model are initialized. The commonly used objective function is the Mean Absolute Error (MAE), and the formula is as follows:

[0144] (12);

[0145] Where, is the mean absolute error, is the true value of the th sample, is the predicted value of the th sample, that is, the estimation of the true value by the GMM-XGB hybrid model.

[0146] Common hyperparameters in the GMM-XGB hybrid model include the learning rate, the maximum depth of the tree structure, the subsample ratio, the number of tree structures, etc. When initializing the model, it is necessary to set the search space of the hyperparameters: the search space of the learning rate is [0.01, 0.2], the search space of the maximum depth of the tree structure is [5, 10], the search space of the subsample ratio is [0.5, 1], and the search space of the number of tree structures is [100, 500].

[0147] Step 1.5.2: After the GMM-XGB hybrid model is initialized, iterative training is carried out to find the optimal hyperparameters of the model within the search space. First, sample collection is performed, that is, a set of initial points is selected as the initial hyperparameters in the space of model hyperparameter combinations to be optimized (search space). Then, the Bayesian optimization algorithm evaluates the next set of hyperparameters to be selected based on the existing hyperparameters and objective function values. The Bayesian optimization algorithm automatically balances exploitation and exploration through an acquisition function (calculating the expected improvement) and a Gaussian process (calculating the predictive variance). This balance makes the Bayesian optimization more efficient than grid search or random search, especially suitable for hyperparameter tuning scenarios with high computational costs.

[0148] Step 1.5.3: Then, the GMM-XGB hybrid model is trained and evaluated. The GMM-XGB hybrid model is trained using the selected optimal hyperparameter combination. The GMM-XGB hybrid model is trained on the training set, and its performance is evaluated on the validation set. The mean absolute error between the predicted values and the true values of the GMM-XGB hybrid model is calculated on the validation set to obtain the objective function value.

[0149] Step 1.5.4: After obtaining the new hyperparameter-objective function pairs, the surrogate model (usually a Gaussian process model) is updated. The surrogate model continuously adjusts its estimate of the distribution of the objective function in the hyperparameter space according to the new data, making the approximation of the objective function more accurate.

[0150] Step 1.5.5: Check whether the pre-set termination conditions are met. A common termination condition is set to reach the pre-set maximum number of iterations. If the termination conditions are not met, return to Step 1.5.2 to perform sample collection again and continue the iteration; when the termination conditions are met, select the set of hyperparameters that makes the objective function value optimal from all the evaluated hyperparameter combinations. Use this set of optimal hyperparameters to retrain the GMM-XGB hybrid model on the training set until the pre-set number of training times is reached. After the training is completed, the final model for prediction or analysis is obtained. Then, the performance of this model can be evaluated on an independent test set to verify the generalization ability of the model.

[0151] The GMM-XGB hybrid model trained through Bayesian optimization will achieve good results on the given training set and validation set, and can accurately predict the distribution of the injection and production fluid volumes of small layers. Finally, the model will output the injection and production fluid volumes of all wells at each time step, refined to the fluid volume distribution of each small layer. These output results will provide guidance for oilfield management and optimization, helping to reasonably allocate the injection and production fluid volumes and improve the development efficiency of the oilfield.

[0152] Step 2: Construct a hierarchical connected network model based on the physical model. Taking the injection-production fluid volume of small layers and other dynamic and static parameters of the reservoir as inputs, the "characteristic point-connected network" method is proposed. Fractures and faults are abstracted as a series of characteristic points and embedded into the connected network between injection and production wells, and a more practical geological structure-connected network is reconstructed. Through the hierarchical discretization technology, the coupling of the matrix-fracture-fault multi-scale flow network is realized. Then, all the connections in the network are discretized into grids and mapped into a two-dimensional rectangular coordinate grid and input into the simulator for solution. The specific process is as follows:

[0153] Step 2.1: First, the reservoir can be equivalent to an inter-well connected network formed by a series of well points and hierarchical inter-well connected pipes. The well points are connected to each other through the connected pipes to form an inter-well connected network. The complexity of this network can be flexibly adjusted by controlling the well spacing and well angle, so that the model can adapt to reservoirs of different scales and complexities. In this network, the inter-well pipe is regarded as a connected unit, which is characterized by two parameters: inter-well conductivity and connected volume. By calculating these parameters, the characteristics of fluid flow in the reservoir can be effectively captured.

[0154] In the process of reservoir development, studying the fluid flow characteristics of the reservoir and the connectivity between injection and production wells is a key link in optimizing oilfield production and improving recovery efficiency. In order to quickly and accurately simulate this process, the reservoir can be regarded as an inter-well connected network composed of a series of well points and hierarchical inter-well connected pipes through an equivalent method. This network model simplifies and abstracts the complex physical process of the reservoir, enabling the analysis of fluid flow and inter-well interaction in a simpler and more intuitive way.

[0155] Specifically, in this network model, the reservoir is divided into several well points and connected pipes. Each well point represents an injection well or a production well on a small layer, and these well points are connected through a series of connected pipes. Each inter-well pipe is equivalent to a flow channel, and fluid can be transferred between the well points connected by it. The characteristics of each pipe are characterized by two main parameters: inter-well conductivity and connected volume. Among them, the inter-well conductivity reflects the ability of fluid to flow in the pipe, and the connected volume characterizes the reserves between this well point and the surrounding well points. These two parameters are the core parameters for studying the fluid dynamics of the reservoir and inter-well flow.

[0156] The conductivity calculation formula is as follows:

[0157] (13);

[0158] Among them, is the conductivity between well and well ; is the well Seepage area between the well and the well , unit is ; is the average permeability between the well and the well , unit is ; is the distance between the well and the well , unit is

[0159] The calculation formula for the connected volume is as follows:

[0160] (14);

[0161] where is the connected volume between the well and the well ; is the average value of the effective thickness between the well and the well .

[0162] The connected pipe structure between well points is a model highly dependent on the specific reservoir conditions. The two geometric parameters of well spacing and well angle directly determine the topological structure of the network. By adjusting these parameters, the complexity of the model can be flexibly controlled, thereby adapting to reservoir conditions of different scales and geological conditions. For example, in a reservoir, if the distance between well points is relatively close and the angle is small, the connectivity between wells will be stronger, and the fluid exchange and conduction ability will also be enhanced; while in the case of a relatively large well spacing or a large well angle, the fluid flow and conduction effect will be relatively weak. Therefore, controlling the well spacing and well angle enables the model to better simulate the inter-well flow in reality, thereby providing a scientific basis for production optimization.

[0163] In addition, this inter-well connected network model can simplify complex reservoir conditions into a mathematical model, thereby providing support for subsequent real-time optimization of injection and production. By independently calculating each connected unit (inter-well pipe), the flow characteristics of fluids in the reservoir and the mutual influence between wells can be effectively captured.

[0164] Step 2.2. Traditional flow network models are usually roughly characterized only by equivalent permeability parameters. The "characteristic point-connected network" method is proposed, which abstracts fractures and faults into discrete characteristic points and embeds them into the inter-well connected network to construct a dynamic flow channel closer to the actual geological structure. By endowing the characteristic points with different conductivity attributes (such as the high permeability of fractures and the blocking effect of faults), the model can quantitatively describe the high-speed seepage effect of fractures on injection and production fluids and the blocking effect of faults on flow paths.

[0165] In a reservoir, fractures and faults are two important geological features that significantly affect the flow paths of oil and gas and the inter-well connectivity, thereby influencing the production performance of the reservoir. To address this issue, in reservoir simulation, the impacts of fractures and faults need to be particularly considered. The "characteristic point - connectivity network" method is proposed. By embedding characteristic points representing fractures and faults into the original inter-well connectivity network model, the fluid flow in the reservoir can be more accurately described, providing a basis for optimizing injection and production management.

[0166] Fractures are high-speed flow channels in a reservoir, and their role in the reservoir is very significant. Fractures can greatly enhance the regional permeability of the reservoir and become the main channels for crude oil flow. Due to the existence of fractures, the flow paths of fluids in the reservoir will change, and the flow velocity will increase significantly. In the original connectivity network, well points are connected by relatively low-permeability channels, while in a reservoir with fractures, as high-speed channels, fractures usually require special treatment.

[0167] Therefore, based on network dissection, fractures can be discretized into a set of uniformly arranged high-permeability nodes. By introducing these fracture nodes on the basis of the original well points and network connectivity pipes, the connectivity relationship between well points and fracture nodes is reconstructed. The high-permeability characteristics of fracture nodes mean that the flow velocity of fluids at these nodes is relatively fast. Therefore, the connectivity between these nodes is usually stronger, the heterogeneous conductivity is larger, and the fluid flow tends to be more concentrated. This method can simulate the impact of fractures on reservoir permeability relatively quickly and accurately, further improving the dynamic simulation and optimization capabilities of the reservoir.

[0168] Faults are another geological structure that affects the flow characteristics of a reservoir. They usually cause the separation or enclosure of the reservoir. The existence of faults often restricts the fluid flow between different regions in the reservoir, especially forming a certain barrier to the migration and accumulation of crude oil. In actual reservoirs, faults have different properties. Some faults are completely closed, resulting in no mutual fluid flow between the two sides; while some faults are partially conductive and have a certain permeability.

[0169] In a reservoir model considering the impact of faults, faults are usually regarded as a closed structure. In this invention, they are discretized into two sets of uniformly arranged nodes, located on both sides of the fault respectively, and the heterogeneous conductivity is used to characterize the magnitude of the fluid flow ability near the fault. Through this discretization method, the closure effect of faults on inter-well connectivity can be accurately reflected, and fault nodes are formed in the network. When reconstructing the connectivity relationship between well points and these fault nodes, the connectivity between well points and the nodes on both sides of the fault will be affected by the fault. For a closed fault, there is almost no fluid conduction between the well point and the fault, while for a conductive fault, there may be a certain flow connection between the well point and the fault, but usually the permeability difference of the fault needs to be considered.

[0170] By discretizing faults into nodes and reconstructing the connectivity between well points and faults, the flow barriers or restrictions caused by faults in the reservoir can be more realistically simulated. This is of great significance for reservoir production optimization, injection-production scheduling, pressure management, etc.

[0171] Step 2.3: Based on the inter-layer well connectivity relationship, multiple connectivity units are set within each layer and one-dimensional grid discretization is performed. Then, according to the grid connection relationship, it is mapped into a standardized model file to construct a physically-driven inter-layer connectivity network model. In the model construction, the coupling of the matrix-fracture-fault multi-scale flow network is realized through the hierarchical discretization technique: for the flow in the matrix, one-dimensional grid discretization is used to describe the conventional seepage between wells; for the flow in fractures and faults, the non-uniform conductivity parameters between characteristic points are used to characterize the flow conductivity of fractures and the occlusion characteristics of faults. This coupling mechanism not only retains the computational efficiency of the flow network model but also realizes the refined description of complex flow paths, solving the simulation deviation problem caused by the simplified flow paths in fractured reservoirs by traditional methods. Finally, the pressures, saturations, injection-production rates, water cut, etc. of each injection-production well are solved through a numerical simulator.

[0172] First, according to the geological characteristics of the reservoir, the reservoir is divided into multiple small layers, and each small layer corresponds to an independent physical model. These small layers are usually divided according to the layer series combination division results given in the field and the well point connectivity between each layer is considered. Within each layer, multiple connectivity units are set based on geological data and well test data. For the flow characteristics within each layer, refinement is carried out through one-dimensional grid discretization. Each grid unit is equivalent to a part of the connectivity unit and can more accurately simulate the fluid flow in this area.

[0173] The flow process of the reservoir is transformed into a discretized numerical calculation problem. The basic idea of grid discretization is to transform the continuous fluid flow process within each connectivity unit into a discrete system composed of multiple small units, and the properties such as permeability, volume, and flow rate of each small unit can be solved by numerical methods. The connection relationship between these discrete units can be mapped through a standardized model file to generate a grid model suitable for numerical solution, and this network model is the inter-layer connectivity network model.

[0174] For the flow in the reservoir matrix, a one-dimensional grid discretization method is adopted to discretize the well-to-well connectivity relationship in each layer into multiple one-dimensional grid cells; for the flow in the fractures, the fractures are abstracted as a series of characteristic points with high heterogeneity conductivity, and the conductivity of the fractures is characterized by the heterogeneity conductivity parameter. The connectivity network between the fracture characteristic points can simulate the high-speed flow of fluids in the fractures and reflect the rapid conduction effect of the fractures on the injected and produced fluids; for the flow near the faults, the faults are abstracted as characteristic points with a blocking effect, and the characteristics of fluid flow near the faults are characterized by the heterogeneity conductivity parameter. The introduction of the fault characteristic points can accurately depict the blocking effect of the faults on the flow path.

[0175] Through the layered discretization technology, the flow networks of the matrix, fractures, and faults are coupled. The characteristic points of the fractures and faults are embedded in the matrix network, and the fluid exchange between the matrix and the fractures and faults is realized through conductivity. The fracture characteristic points serve as high-speed channels for matrix flow, and the fault characteristic points play a blocking role in the fluid near the faults, avoiding the poor simulation effect caused by the simplified characterization of fractures and faults in previous methods.

[0176] This grid discretization not only helps to improve the simulation accuracy but also facilitates parameter adjustment during the calculation process. The connection relationship between each grid cell and its surrounding cells determines the flow path and velocity of the fluid, and the changes in these connection relationships will directly affect the production dynamics of the entire reservoir.

[0177] By mapping the discretized grid cells to a standardized model file, a physically driven layered connectivity network model can be finally constructed. This model not only performs layered processing on the reservoir but also embeds the fractures and faults existing in the actual complex reservoir, further reducing the resource cost brought by numerical simulation.

[0178] Step 3: Use the ES-MDA algorithm for automatic history matching. First, select the daily oil production of a single oil well as the fitting target. On the basis of selecting conductivity, grid volume, and relative permeability parameters as fitting parameters, add the heterogeneity conductivity and the number of characteristic points of the fractures and faults to enable adaptive correction of the dynamic effects of the fractures and faults. Construct a mathematical model of the automatic history matching objective function and add engineering constraints. Use the ES-MDA algorithm to iteratively solve the objective function to make the daily oil production simulated by the layered connectivity network model coincide with the actually observed daily oil production while conforming to the actual geological understanding. The specific process is as follows:

[0179] Step 3.1. In the automatic history matching of oil reservoirs, the construction of the objective function is one of the core links, which determines the effect and efficiency of the matching process. Traditional history matching methods often rely on manual parameter adjustment and gradually optimize model parameters through the trial-and-error method. To improve the automation level and accuracy of the matching, a history matching objective function is constructed based on the Bayesian Maximum A Posteriori Estimation (MAP) method, which can effectively integrate the physical model and prior knowledge, thereby achieving more accurate history matching.

[0180] In the process of oil reservoir history matching, the daily oil production of oil wells is usually used as the observed data to make the model output as close as possible to the observed data. The constructed history matching objective function is:

[0181] (15);

[0182] Where, is the value of the history matching objective function, that is, the difference between the model simulation calculation result and the actual observed data; is the model fitting parameter vector; is the model simulation calculation result; is the actual observed data; is the covariance matrix of the observation error.

[0183] By minimizing the value of the history matching objective function, the coincidence between the model calculation result and the actual observed data is realized.

[0184] Step 3.2. In the process of oil reservoir history matching, it is crucial to select appropriate control variables. In the oil reservoir model, the transmissibility between discrete grids, the volume of discrete grids, the heterogeneous transmissibility and characteristic points of fractures, the heterogeneous transmissibility and characteristic points of faults, and the relative permeability parameters are usually selected as control variables to form the model fitting parameter vector;

[0185] (16);

[0186] Where, is the model fitting parameter vector, which contains the transmissibility between all grids and the volume of the grids; is the transmissibility of grid and grid at the initial time, with the unit of ; is the volume of grid at the initial time, with the unit of ; are the characteristic points and characteristic point of the fracture at the initial time, at The heterogeneous conductivity on the layer, with the unit of ; is the characteristic point number of the fracture on the layer at the initial time; is the characteristic point of the fault at the initial time and the characteristic point on the layer, with the unit of ; is the characteristic point number of the fault on the layer at the initial time; , , , are different relative permeability parameters;

[0187] In the process of reservoir history matching, some reasonable physical constraints must be imposed while minimizing the objective function. This can ensure that the model parameters obtained in the matching process can not only match the observed data but also follow the physical characteristics of the actual reservoir. Considering the reservoir characteristic constraints on the variation range of model parameters, the specific constraints can include:

[0188] (17);

[0189] Among them, is the conductivity of grid and grid ; is the volume of grid ; is the number of grids; is the total volume of the reservoir, with the unit of .

[0190] Ensure that the conductivity between wells is always greater than 0, the grid volume is always greater than 0 and less than the sum of all grid volumes, and the sum of all grid volumes is equal to the total volume of the reservoir, , is always between 1 and 6, , is always between 0 and 1.

[0191] Step 3.3. During the development of complex reservoirs, fractures and faults have a crucial impact on fluid flow characteristics. However, due to the high uncertainty in the spatial distribution of fractures and faults and the difficulty in directly obtaining their parameters through conventional logging or seismic data, accurately characterizing and modeling them pose a great challenge. To overcome the challenge of difficult direct acquisition of fracture and fault parameters, the present invention adopts the ES-MDA (Ensemble Smoother with Multiple Data Assimilation) algorithm, combines production history data and physical constraint conditions, and realizes the automatic inversion and real-time update of fracture and fault parameters.

[0192] The ES-MDA algorithm is a multiple data assimilation technology based on the ensemble method, which can continuously optimize reservoir parameters during the history matching process, making the prediction results of the model gradually approach the actual production situation. First, set the initial heterogeneous conductivity distribution of fractures and faults in the reservoir numerical model, and at the same time define the characteristic points to describe the complex geometric characteristics of fractures and faults. Then, during the history matching process, use the actual daily oil production of individual wells in the production data of oil wells as the observation information to perform multiple data assimilations on the model. Through the iterative update mechanism of the ES-MDA algorithm, continuously adjust the matrix parameters and the parameters of fractures and faults, so that the fitting degree between the simulation results and the actual production data gradually improves.

[0193] Through the deep integration of data-driven and physical mechanisms, this method enables the reservoir model to adaptively adjust the dynamic influence of fractures and faults on fluid flow, thereby improving the prediction ability of the model for the injection-production response of complex reservoirs. Compared with the traditional manual adjustment method, the automatic inversion and real-time update mechanism of the ES-MDA algorithm not only improves the calculation efficiency, but also can more accurately capture the change characteristics of the internal structure of the reservoir, providing more reliable technical support for reservoir dynamic management and optimized development.

[0194] Step 4. Establish a mathematical model for the real-time optimization of layered injection-production objective function, and perform real-time optimization of layered injection-production based on the SHADE algorithm. First, select the economic net present value as the objective of real-time optimization of injection-production, and the injection-production fluid volumes of individual layers of individual wells as the optimization variables. Considering the actual engineering constraints, establish a mathematical model for the real-time optimization of layered injection-production objective function, introduce an adaptive mechanism based on successful history and the SHADE algorithm with an archive mechanism, significantly improve the solution efficiency and accuracy of the injection-production optimization model, and maximize the economic benefits by iteratively maximizing the optimization objective. The specific process is as follows:

[0195] Step 4.1: Real-time injection-production optimization aims to maximize the oil recovery rate and economic benefits of the oilfield by continuously monitoring the production status of the oilfield and adjusting the water injection and oil production strategies. It can also make the oilfield management more flexible and intelligent, providing strong support for the long-term stable development of the oilfield. This process involves dynamically adjusting the water injection volume, liquid production volume, and the working status of wells, so as to optimize the resource allocation and enhance the overall economic benefits while ensuring the long-term stable production of the oilfield. In the actual oilfield production process, with the progress of exploitation, the geological, physical properties of the reservoir and its fluid distribution will change, and traditional injection-production optimization methods are difficult to cope with these dynamic changes. Therefore, real-time injection-production optimization has become an effective means to address this challenge.

[0196] To achieve the above goals, it is reasonable to select the economic net present value (NPV) as the objective function for real-time injection-production optimization. The economic net present value is obtained by discounting future cash flows to the current moment, comprehensively considering the investment return and production cost, so as to obtain the economic benefits of the oilfield project. Using the economic net present value as the objective function can effectively balance the relationship between the increase in oil recovery rate and economic benefits, ensuring that while pursuing long-term stable production, the economic return is maximized. Through real-time optimization, the economic net present value can be used as a direct indicator to measure the optimization effect, providing a scientific decision-making basis for the sustainable development of the oilfield. The following is the established mathematical model of the optimization objective function:

[0197] (18);

[0198] Where, is the economic net present value, with the unit of yuan; is the control variable, representing the regulation scheme of the layered water injection volume and layered liquid production volume during the optimization process; is the oil price, with the unit of yuan; is the th oil production of well , with the unit of ; is the sewage treatment cost, with the unit of yuan; is the th water production of well , with the unit of ; is the th water injection volume of well , with the unit of ; is the total number of production time steps;

[0199] Step 4.2. During the reservoir development process, the complexity and variability of the reservoir make the formulation and optimization of the development plan challenging. To ensure the rationality and feasibility of the development plan, relevant parameters need to be reasonably constrained according to the specific situation of the reservoir in practical applications. Through these constraints, not only can the accuracy of reservoir numerical simulation and optimization process be improved, but also unrealistic results can be avoided in actual operations. The following are some common constraints on reservoir development parameters, which help to improve the production stability and economic benefits while ensuring the reservoir development effect.

[0200] A reasonable injection-production ratio can ensure the coordination between water injection and oil production activities in the reservoir, and avoid the phenomena of over-injection or over-production. An excessively high injection ratio may lead to water flooding problems in the reservoir, affecting the long-term production potential of the reservoir; an excessively low injection ratio may cause insufficient displacement in the reservoir, affecting the recovery efficiency.

[0201] The bottom-hole flowing pressure directly affects the production capacity of the reservoir and the working state of the well. Reasonable control of the bottom-hole flowing pressure can avoid problems such as wellbore collapse and downhole equipment damage, and at the same time ensure the stable production of the reservoir. An excessively high bottom-hole flowing pressure may cause excessive wellbore pressure, affecting the oil production effect; an excessively low bottom-hole flowing pressure may lead to a decrease in the pumping capacity of the well and unable to achieve the ideal production.

[0202] As the exploitation progresses, the proportion of water often gradually increases, which is also a common problem in the reservoir development process. An excessively high water cut means that more water is produced, resulting in difficult oil-water separation, deteriorated oil quality, and increased later treatment costs. Therefore, a maximum water cut constraint is usually set in reservoir development to prevent the water cut from exceeding a certain critical value, thereby ensuring the economy of the production process. This constraint is usually set based on the economic analysis of the oilfield, equipment capacity, and production requirements.

[0203] The single-well limit liquid volume constraint is set based on the pumping capacity of the well, the fluid transportation capacity of the wellbore, and the processing capacity of the wellhead facilities. An excessively high liquid volume may exceed the bearing capacity of the wellbore, resulting in equipment failures or decreased production efficiency. At the same time, an excessively low liquid volume may also lead to insufficient development of the reservoir and unable to achieve the optimal recovery effect.

[0204] The main constraints include injection-production ratio constraint, bottom-hole flowing pressure constraint, maximum water cut constraint, and single-well limit liquid volume constraint. The specific constraints are as follows:

[0205] (19);

[0206] In the formula, is the minimum injection-production ratio; is the maximum injection-production ratio; is the Total water injection volume for a time step, unit ; is the liquid production volume or water injection volume of well at the time step; is the total liquid production volume at the time step, unit ; is the minimum bottom-hole flowing pressure of well at the time step, unit ; is the maximum bottom-hole flowing pressure of well at the time step, unit ; is the bottom-hole flowing pressure of well at the time step; is the water cut of well at the time step; is the ultimate water cut of well at the time step; is the liquid production volume or water injection volume of well at the time step, unit ; is the ultimate liquid production volume or ultimate water injection volume of well at the time step, unit .

[0207] Step 4.3. In the optimization process of the traditional DE algorithm, the scaling factor (F) and crossover probability (CR) need to be manually set, and its fixed parameters are difficult to adapt to the dynamic characteristics of complex injection-production optimization problems. Aiming at the limitations of the traditional differential evolution algorithm (DE) in dealing with high-dimensional and non-linear optimization problems, such as being prone to falling into local optimum and slow convergence speed, an adaptive mechanism based on success history and the SHADE algorithm (Success-History based Adaptive Differential Evolution) with an archive mechanism are introduced. At the same time, considering engineering constraints to solve the optimization objective function to obtain the maximum economic net present value, the solving efficiency and accuracy of the injection-production optimization model are significantly improved, which helps to adjust the injection-production system in a timely manner according to different stages and real-time production data during the oilfield development process, solve problems such as interlayer interference and well-to-well connectivity, and improve production efficiency.

[0208] The SHADE algorithm belongs to one of the modern evolutionary algorithms. Its core idea is to adjust the search strategy by introducing an adaptive mechanism, thereby effectively improving the convergence speed and global search ability of the optimization algorithm. SHADE dynamically adjusts the mutation parameters by utilizing historical successful experiences (successful mutation operations) in each generation of the algorithm, which enables SHADE to exhibit high robustness in complex optimization problems. The SHADE algorithm dynamically adjusts the values of F and CR by introducing a historical success parameter adaptive mechanism, enabling it to automatically adapt to problem characteristics based on the feedback information during the optimization process. This adaptive mechanism not only reduces the dependence on manual parameter tuning but also significantly improves the robustness and convergence efficiency of the algorithm in injection-production optimization. The SHADE algorithm also introduces a historical memory archive (Memory Archive) to store the parameter information of historical successful individuals and generate new trial individuals based on this information. This mechanism enables the algorithm to make full use of historical optimization experiences, avoid ineffective searches, and thus accelerate the convergence speed. In injection-production optimization, this feature is particularly important because it can quickly locate high-potential solution regions and reduce the waste of computing resources.

[0209] The SHADE algorithm includes mutation operations, crossover operations, and selection operations; in reservoir injection-production optimization, the mutation operation can generate new solutions by adjusting injection-production flow rates, while the crossover operation can combine different injection-production schemes to obtain more possible solutions. In this way, the SHADE algorithm can generate new solutions that meet the optimization requirements while ensuring the constraint conditions.

[0210] Individuals in each generation will undergo fitness evaluation and be sorted according to the objective function value and the satisfaction of constraints. The fitness evaluation not only considers the objective function (maximizing the economic net present value) but also ensures that engineering constraint conditions (such as injection-production ratio, bottom-hole flowing pressure, maximum water cut, etc.) are met. Those solutions that meet the engineering constraint conditions and can maximize the objective function will be preferentially selected.

[0211] Based on the fitness evaluation results, the selection operation will select those individuals with higher fitness to enter the next generation. These individuals will undergo mutation and crossover operations according to the adaptive mechanism to generate new solutions. In reservoir injection-production optimization, the selection operation can select the best injection-production scheme by evaluating the objective function and constraint conditions, thereby guiding the actual operations in oilfield development.

[0212] To prove the feasibility and superiority of the present invention, the following specific embodiments are given.

[0213] In the embodiment of the present invention, a one-injection and four-production model with two layers containing fractures and faults is established by using a finite volume method reservoir numerical simulator. Its grid scale is 25*25*2, the model size is 100m*100m*20m, the reservoir porosity is 0.3, the initial water saturation of the reservoir is 0.25, the oil and water viscosities are 4mPa*s and 1mPa*s respectively. There is 1 injection well and 4 production wells in the reservoir. The injection well number is 5, and the production well numbers are 1, 2, 3, and 4. The production is simulated for 4770 days by using the reservoir numerical simulator, and the overall injection-production balance of the block is achieved.

[0214] Then, as Figure 2 shown, the method of the present invention is used to predict the injection-production fluid volume of sub-layers. The specific process is as follows: First, reservoir data is collected to construct a sample library for longitudinal splitting of injection-production fluid volume. The production and injection profiles in the data are used as label data, and porosity, permeability, oil layer thickness, initial water saturation, formation factor, number of perforations, perforation thickness, permeability range, porosity range, bottom hole flowing pressure, single-well injection-production fluid volume, and water cut are used as input features. Then, a Gaussian mixture model is constructed to output the probability distribution of the injection-production fluid volume of sub-layers. The reservoir data is fused with the probability distribution of the injection-production fluid volume of sub-layers to construct a GMM-XGB hybrid model, and the Bayesian optimization algorithm is used to train the model. Finally, the predicted injection-production fluid volume of sub-layers is output. In the embodiment of the present invention, the whole-well injection-production volume is split into two sub-layers based on the above process.

[0215] Then, a hierarchical connectivity network model is established using the same reservoir parameters. The reservoir is equivalent to an inter-well connectivity network formed by a series of well points and inter-layer well-to-well connecting pipes. The well points are connected to each other through connecting pipes to form an inter-well connectivity network. The complexity of this network can be flexibly adjusted by controlling the well spacing and well angle. In this network, the inter-well pipes are regarded as a connectivity unit, which is characterized by two parameters: inter-well conductivity and connectivity volume. Then, the fractures and faults are abstracted into discrete feature points and embedded into the inter-well connectivity network to reconstruct the connectivity relationship between the well points and the fracture and fault nodes. The specific processing method is as Figure 3 shown. First, five well nodes are used to represent five well points. The gray diamond structure is a fault, and the gray rectangular structure is a fracture. Then, the fracture and the fault are respectively abstracted into three fracture nodes and six fault nodes. Finally, the connectivity relationship between the five well nodes (with serial numbers from 1 to 5), the fracture nodes, and the fault nodes is reconstructed to obtain the required hierarchical connectivity network model. Finally, the coupling of the matrix-fracture-fault multi-scale flow network is realized through hierarchical discretization technology: for the flow of the matrix, one-dimensional grid discretization is used to describe the conventional inter-well seepage; for the flow of fractures and faults, the non-homogeneous conductivity parameters between feature points are used to characterize the diversion ability of fractures and the blocking characteristics of faults.

[0216] The single forward running time of the model is 32 s. Compared with the finite volume reservoir numerical simulator which takes 123.3 s, the calculation speed is increased by about 4 times. The ES-MDA algorithm is used for history matching. The population size is selected as 200, and the iteration is carried out 3 times. Considering the constraint conditions, during the history matching, the comparison images of the prior model calculation results and the actual observation data of all sub-layers are output. Here, only the comparison image of the prior model calculation result and the actual observation data of sub-layer 1 of production well 1 is shown. Figures 4 - 6 They are the history matching results at the 1st, 2nd, and 3rd iterations respectively. It can be seen that during the history matching, the daily oil production curves (corresponding to the prior model calculation results) obtained by each set of matching parameters envelope the daily oil production curves (corresponding to the actual observation data) obtained by the numerical simulator. Figures 4 - 6 In the figure, one population corresponds to a set of matching parameters. Therefore, the figure contains the daily oil production curves obtained by 200 sets of matching parameters corresponding to 200 populations, and some gray dotted lines overlap.

[0217] After history matching, the fitting curve of the block water cut is output as shown in Figure 7 It can be seen that the true value basically coincides with the fitting value output by the numerical simulator after history matching.

[0218] After history matching, the fitting curve of sub-layer 1 of production well 1 is output as shown in Figure 8 It can be seen that the daily oil production rate of production well 1 using the method of the present invention basically coincides with the daily oil production rate output by the numerical simulator.

[0219] Using the model after history matching of the present invention, the SHADE algorithm is combined with engineering constraints to optimize the layered injection and production. A total of 12 time steps are optimized, each time step is 30 days, and a total of 360 days are optimized. Here, as shown in Figure 9 and Figure 10 only the injection regime regulation diagrams of sub-layer 1 and sub-layer 2 of injection well 5 are shown. From Figure 9 and Figure 10 it can be seen that restricted by the 20% upper and lower limits of the single-well liquid volume of the injection well, the optimized injection rate is restricted between 0.8 and 1.2 times of the injection rate before optimization.

[0220] After the injection and production regimes of all wells and all layers are regulated, the overall optimization results are as shown in Figure 11 and Figure 12 It can be seen that Figure 11 and Figure 12 are the optimization results of the cumulative oil production and water cut respectively. It can be seen that the optimized cumulative oil production is significantly increased and the water cut is slightly decreased, thus verifying the reliability of the method of the present invention in complex reservoir simulation and optimization calculation.

[0221] Certainly, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples either. Changes, modifications, additions or substitutions made by those skilled in the art within the scope of the essence of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for rapid simulation and real-time optimization of stratified injection and production applicable to complex oil reservoirs, characterized in that: The steps include: Step 1: Construct a GMM-XGB mixture model based on the Gaussian mixture model and XGBoost algorithm; Step 2: Construct a hierarchical connected network model based on the physical model; the specific process is as follows: Step 2.1, the oil reservoir is divided into several well points and connecting pipes, and an inter-well connecting network is constructed; in the inter-well connecting network, each well point represents a water injection well or oil production well on a small layer, and the well points are connected to each other through connecting pipes; the connecting pipes are characterized by two parameters: inter-well conductivity and connecting volume; The conductivity calculation formula is as follows: ； in, Well and well The conductivity between Well and well The seepage area between Well and well The average value of the permeability between Well and well The distance between The formula for calculating the connected volume is as follows: ； in, Well and well The connected volume between them; Well and well The average value of the effective thickness between Step 2.2, discretize the fracture into a group of uniformly arranged high permeability nodes; discretize the fault into two groups of uniformly arranged nodes, located on both sides of the fault, and characterize the flow capacity of the fluid near the fault by heterogeneous conductivity; embed the fracture and fault nodes into the well-to-well connection network to construct a dynamic flow channel; Step 2.3, based on the hierarchical well connectivity relationship, multiple connectivity units are set in each layer and one-dimensional grid discretization is performed, and then the grid connectivity relationship is mapped into a standardized model file to construct a physically driven hierarchical connectivity network model; Step 3: Use ES-MDA algorithm for automatic history matching; Step 4: Establish a mathematical model of the objective function for real-time optimization of stratified injection and production, and perform real-time optimization of stratified injection and production based on the SHADE algorithm.

2. The rapid simulation and real-time optimization method for layered injection and production applicable to complex oil reservoirs according to claim 1, characterized in that: In step 1, the GMM-XGB hybrid model uses a Gaussian mixture model to perform cluster analysis on the original sample set to obtain the probability distribution of the small layer injection and production volume, and integrates the probability distribution into the original sample set as the input of the XGBoost algorithm to predict the output of the small layer injection and production volume; the specific working process of the GMM-XGB hybrid model is as follows: Step 1.1, collect the longitudinal splitting data of injection and production volume, and construct the original sample set; The vertical splitting data of injection and production volume include porosity, permeability, oil layer thickness, initial water saturation, formation coefficient, number of perforations, perforation thickness, permeability extreme difference, porosity extreme difference, bottom hole flowing pressure, single well injection and production volume and water content; Step 1.2: Preprocess the original sample set. The specific process is as follows: First, data cleaning is performed. For missing values ​​of random missing types, the natural neighbor interpolation method in MATLAB software is used for interpolation according to the geodetic coordinates of the well points. For outliers, deletion or transformation processing is selected according to the actual situation. Transformation processing includes smoothing or scaling. For categorical data, the number of base classes of categorical data needs to be determined first. If there are few base classes, one-hot encoding is used for processing; if there are many base classes, label encoding or target encoding is used for processing. Then, the cleaned data was normalized using the Min-Max normalization method; Finally, the normalized sample set was subjected to correlation analysis; for linear relationships in the data, the Pearson correlation coefficient was used for correlation analysis; for sample sets with strong nonlinear relationships and no normal distribution, the Spearman rank correlation coefficient was used for correlation analysis; Step 1.3, using Gaussian mixture model to perform cluster analysis on the preprocessed original sample set to obtain the probability distribution of small layer injection and production volume; Step 1.4, the probability distribution of the small layer injection and production volume is fused with the preprocessed original sample set to obtain a fused sample set, and the fused sample set is divided into a training set and a validation set; Step 1.5: Input the fused sample set into the XGBoost algorithm and use the Bayesian optimization method to train until the requirements are met, and finally output the injection and production volume of all small layers in all wells and all time steps.

3. The rapid simulation and real-time optimization method for layered injection and production applicable to complex oil reservoirs according to claim 2, characterized in that: The specific process of step 1.3 is as follows: Step 1.3.1: Based on the preprocessed original sample set, a Gaussian mixture model is constructed; the Gaussian mixture model assumes that the injection and production volume data points come from multiple different Gaussian distributions, and the probability density function of the Gaussian distribution is: ； in, is the probability density function of Gaussian distribution; is the data point of injection and production volume; is the mean vector; is the covariance matrix; is the data variance; is the data dimension; is an exponential function with base e; Step 1.3.2: Use the expectation maximization algorithm to iteratively solve the established Gaussian mixture model until convergence. The specific process is as follows: First, initialize the parameters of each Gaussian distribution in the Gaussian mixture model; Then, the responsiveness of all injection and production volume data points is calculated. The responsiveness calculation formula is as follows: ； in, is responsiveness; is a binary indicator variable, indicating Injection and production volume data points Whether Gaussian components are generated; For the The prior probability of the Gaussian components; For the The mean vector of the Gaussian components; For the The covariance matrix of the Gaussian components; is the index variable of the Gaussian component, which is used to traverse all Gaussian components in the mixed Gaussian model for summation; is the total number of Gaussian components in the mixed Gaussian model; For the The prior probability of the Gaussian components; For the The mean vector of the Gaussian components; For the The covariance matrix of the Gaussian components; The parameters of each Gaussian distribution are updated according to the calculated responsiveness. The update formula is as follows: ； ； ； in, , , Respectively The first iteration The mean vector, covariance matrix, and prior probability of the Gaussian components; For the Responsiveness at iterations; is the total number of data points of injection and production volume; is the transpose symbol; Calculate the likelihood function of the model, the formula is as follows: ； in, is the likelihood function; Finally, the model parameters are continuously updated iteratively until the parameters of the Gaussian mixture model converge. At this time, the training ends and the final responsiveness is output. The responsiveness is the probability distribution of the injection and production volume of the small layer.

4. The rapid simulation and real-time optimization method for layered injection and production applicable to complex oil reservoirs according to claim 3, characterized in that: In step 1.5, Bayesian optimization is an automated hyperparameter tuning method that establishes a proxy model and continuously adjusts hyperparameters based on the historical performance of the model to find the optimal hyperparameter configuration. The specific process is as follows: Step 1.5.1, initialize the Bayesian optimization algorithm framework, initialize the evaluation of the objective function in the framework and set the hyperparameter search space of the GMM-XGB hybrid model; the objective function uses the mean absolute error, and the formula is as follows: ； in, is the mean absolute error; For the The true value of samples; For the The predicted value of samples; is the total number of samples; Hyperparameters include learning rate, maximum tree depth, subsample ratio, and number of tree structures. The search space of learning rate is set to [0.01, 0.2], the search space of maximum tree depth is set to [5, 10], the search space of subsample ratio is set to [0.5, 1], and the search space of number of tree structures is set to [100, 500]. Step 1.5.2, iterative training, to find the optimal hyperparameters of the model in the search space; firstly, sample collection is performed, that is, a set of initial points are selected as the initial hyperparameters in the model hyperparameter combination space to be optimized, and then the Bayesian optimization algorithm evaluates the next set of hyperparameters to be selected based on the existing hyperparameters and objective function values; finally, the hyperparameter combination with the largest objective function is selected as the optimal hyperparameter combination; Step 1.5.3: Use the selected optimal hyperparameter combination to train the GMM-XGB hybrid model on the training set, calculate the mean absolute error between the predicted value and the true value of the GMM-XGB hybrid model on the validation set, and obtain the objective function value; Step 1.5.4: After obtaining the new hyperparameter-objective function pair, update the surrogate model. The surrogate model will continuously adjust its estimate of the distribution of the objective function in the hyperparameter space based on the new data. Step 1.5.5: Determine whether the pre-set termination condition is met; if the termination condition is not met, return to step 1.5.2 to collect samples again and continue iterating; when the termination condition is met, select a set of hyperparameters that optimize the objective function value from all evaluated hyperparameter combinations, and use this set of optimal hyperparameters to retrain the GMM-XGB hybrid model on the training set until the pre-set number of training times is reached; after the training is completed, the final GMM-XGB hybrid model for prediction is obtained.

5. The rapid simulation and real-time optimization method for layered injection and production applicable to complex oil reservoirs according to claim 4, characterized in that: The specific process of step 2.3 is as follows: First, the reservoir is divided into multiple small layers according to its geological characteristics, and each small layer corresponds to an independent physical model. In each layer, multiple connected units are set based on geological data and well test data. The flow characteristics in each layer are refined by one-dimensional grid discretization, and each grid unit is equivalent to a part of the connected unit. The flow process of the reservoir is converted into a discretized numerical calculation problem. The basic idea of ​​grid discretization is to convert the continuous fluid flow process in each connected unit into a discrete system composed of multiple small units. The permeability, volume and flow rate of each small unit are solved by numerical methods. The connection relationships between discrete units are mapped through standardized model files to generate a hierarchical connected network model.

6. The method for rapid simulation and real-time optimization of stratified injection and production applicable to complex oil reservoirs according to claim 5, characterized in that: The specific process of step 3 is as follows: Step 3.1: Construct the history matching objective function as: ； in, is the historical fitting objective function value; is the model fitting parameter vector; Calculate the results for the model simulation; is the actual observation data; is the covariance matrix of the observation error; Step 3.2, select the conductivity between discrete grids, the volume of discrete grids, the heterogeneous conductivity and the number of characteristic points of fractures, the heterogeneous conductivity and the number of characteristic points of faults, and the phase permeability parameter as control variables to form a model fitting parameter vector; ； in, The initial moment grid and Grid The conductivity between The initial moment grid Volume; is the characteristic point of the crack at the initial moment and feature points exist Inhomogeneous conductivity across the layer; For the initial moment The number of characteristic points of cracks on the layer; is the characteristic point of the fault at the initial moment and feature points exist Inhomogeneous conductivity across the layer; For the initial moment The number of characteristic points of the fault on the layer; , , , are different phase permeability parameters; In the process of reservoir history matching, physical constraints are imposed while minimizing the objective function, including: ； in, For Grid and Grid conductivity; For Grid Volume; is the number of grids; is the total volume of the reservoir; Step 3.3: Use the ES-MDA algorithm, combined with historical production data and physical constraints, to achieve automatic inversion and real-time update of fracture and fault parameters. The specific process is as follows: Firstly, the initial heterogeneous conductivity distribution of fractures and faults is set in the reservoir numerical model, and the number of characteristic points is defined. Then, in the history fitting process, the actual production data of oil wells (daily oil production per well) is used as observation information to perform multiple data assimilation on the model. Through the iterative update mechanism of the ES-MDA algorithm, the matrix parameters as well as the parameters of fractures and faults are continuously adjusted.

7. The method for rapid simulation and real-time optimization of stratified injection and production applicable to complex oil reservoirs according to claim 6, characterized in that: The specific process of step 4 is as follows: Step 4.1: Select the economic net present value as the objective function of real-time optimization of injection and production; the following is the mathematical model of the optimization objective function: ； in, is the economic net present value; is the controlled variable, representing the regulation scheme of stratified water injection volume and stratified liquid production volume during the optimization process; for oil prices; For the Time step well Oil production; Cost of sewage treatment; For the Time step well The water production; is the water injection cost; For the Time step well The amount of water injected; is the total number of time steps; is the number of producing wells; is the number of water injection wells; Step 4.2: Set engineering constraints, including injection-production ratio constraint, bottom hole pressure constraint, maximum water cut constraint, and single well maximum liquid volume constraint. The specific constraints are as follows: ； In the formula, is the minimum injection-production ratio; is the maximum injection-production ratio; For the Time step well The amount of liquid produced or water injected; For the Time step well Minimum bottom hole flowing pressure; For the Time step well Maximum bottom hole flowing pressure; For the Time step well Bottom hole pressure; For the Time step well The moisture content; For the Time step well The ultimate moisture content; For the Time step well The amount of liquid produced or water injected; For the Time step well The maximum liquid production or maximum water injection volume; Step 4.3: Introduce the adaptive mechanism based on success history and the SHADE algorithm with archive mechanism to optimize reservoir injection and production.

8. The method for rapid simulation and real-time optimization of stratified injection and production applicable to complex oil reservoirs according to claim 7, characterized in that: In step 4.3, the specific process of reservoir injection and production optimization by SHADE algorithm is as follows: Step 4.3.1, adjusting the injection and production flow rate to generate a new injection and production plan; Step 4.3.2, combine different injection and production schemes and calculate the economic net present value of each injection and production scheme combination; Step 4.3.3: Select the injection-production scheme combination that meets the engineering constraints and has the largest economic net present value as the individual to enter the next generation. The individual will perform steps 4.3.1 and 4.3.2 again according to the adaptive mechanism until the number of iterations is reached; after the iteration, the final injection-production scheme combination is the optimal injection-production scheme.

Images