Multi-objective optimization for oil-gas production process based on compensation model of comprehensive energy consumption using improved evolutionary algorithm

This paper establishes an error compensation multi-objective optimization model of oil-gas production process for optimizing these production indices, including overall oil production, overall water production and comprehensive energy consumption per ton of oil. In order to reduce the error between the model output and the actual value of comprehensive energy consumption per ton of oil, combining the mechanism model with least squares support vector machine (LS-SVM) error model optimized by Bayesian optimization algorithm (BOA), a hybrid model is established to predict the comprehensive energy consumption, in which the mechanism model is used to describe the overall characteristics of oil-gas production process, and LS-SVM error model is established to compensate the mechanism model error. Then, in order to improve the performance of Pareto non-dominated solutions, an improved non-dominated sorting genetic algorithm-II with multi-strategy improvement (IMS-NSGA-II) is proposed to solve the error compensation multi-objective optimization model. Finally, the effectiveness and superiority of the the proposed optimization method are verified by the experiment results on some stand test problems and the optimization problem for the oil-gas production process in a block of an oil production operation area.


Introduction
Although technological progress has increased the diversity of energy alternatives, fossil energy is still the main energy in the modern society. Oil-gas is a kind of finite and scarce resource, and modern society has a strong dependence on it. In fact, the world's demand for petroleum continues to grow, but with the continuous exploitation of petroleum, the world's petroleum reserves continue to decline. In addition, with the continuous exploitation of oil, the industry of oil-gas production generally faces with some challenges such as the relative decline of formation energy, the increase of comprehensive water content, the gradual decline of surface gathering system load, the low efficiency and the gradual increase of the cost per ton of oil. These problems have caused the waste of resources and environmental pollution in the oil industry, and seriously affected the sustainable development of the oil industry. Therefore, in order to improve production, save energy, and make the production process safer and more efficient, advanced modeling, control and optimization technology have become extremely important and indispensable in the oil-gas production process (Al-Fehdly et al., 2019). Due to the complexity and uncertainty of oil-gas production system (Foss, 2012;Kosmidis et al., 2005), the modeling and optimization of oil-gas production process is a formidable challenge for comprehensive optimization of the process (Codas et al., 2012;Lei et al., 2019). Since the energy consumption is usually an important objective of optimization of oil-gas production process, and the energy consumption of pumping system and oil-gas gathering system accounts for a very large proportion of energy consumption in oil-gas production process, many scholars have proposed the optimization models and energy-saving methods of relevant systems.
In the optimization of sucker-rod pumping system, Zheng and Deng (2007) build an optimization model with the minimum energy consumption as objective function. Comprehensively considering the constraints of oil well productivity coordination, oil well production allocation, pumping unit bearing capacity, sucker rod string strength and motor power utilization, and taking the maximum system efficiency as the objective function, Dong et al. (2008a) propose the corresponding optimization simulation method of suction parameters for sucker-rod pumping. In order to overcome some shortcomings of single well research, Dong et al. (2010) also put forward the overall optimization method of an oil production block with the minimum overall input power as the objective function and the overall oil production as the constraint. The actual application shows the superiority of this optimization method. Considering the influence of plunger motion law, valve ball motion law, oil well inflow performance, etc, Wang et al. (2018) establish the dynamic parameter simulation model of pumping wells and study the methods to improve the system efficiency. For gas injection wells, Jun et al. (2017) use Monte Carlo simulation method to select the optimal gas recovery scheme, then genetic algorithm is used to find the optimal development scheme. In order to overcome the difficulties brought by the low oil prices and provide meaningful and valuable information for operators. Rui et al. (2018) develop an integrated evaluation model by analyzing a large number of actual historical data (from the subsurface reservoir to infield flow line).
For the optimization of oil-gas gathering and transferring system, some optimization methods have been reported (Chebouba et al., 2009;Cui et al., 2008;Dong et al., 2017;He et al., 2007;Liu and Chen, 1999;Tan et al., 2016;Wu, 2011;Wu et al., 2014), these methods are mainly to pursue the minimum production operation cost and energy consumption (Cheng et al., 2020). For example, in order to minimize the capital expenditure such as pipeline cost, power consumption cost and heat energy cost, Liu and Chen (1999) establish a mathematical optimization model of oilfield surface pipeline system by using fuzzy optimization method. He et al. (2007) establish the optimization model of gathering and transportation for electrical heating about ramiform net, in which the lowest investment and operating cost are taken as the objective functions. The experiment results show that the method can reduce the total investment and heating cost of gathering and transferring system. According to the characteristics of energy consumption increasing year by year in the gathering and transmission system, taking the tree-like water-mixed gathering and transportation network as the research object, and taking the water-mixed quantity, water-mixed temperature and water-mixed pressure of the network as the parameter variables, the genetic algorithm is used to optimize the parameters (Cui et al., 2008) to make the energy consumption of the network operation significantly lower than before optimization. For heavy oil pipeline heating transportation, Wu (2011) establishes the corresponding energy consumption model, which significantly reduces the energy consumption. Balancing the maximum operation benefits and the maximum transmission capacity, Wu et al. (2014) built an optimization model of natural gas pipelines by using the weight sum method, and according to the nonlinear characteristics of the optimization model, they propose an improved particle swarm optimization (PSO) algorithm to solve the optimization model. In view of the complex composition and behavior of gathering and transportation system, based on a large number of data accumulated in the process of oil and gas gathering and transportation, Tan et al. (2016) establish time series prediction model of production energy consumption in gathering and transportation system, and study the correlation between energy consumption factors. Wang et al. (2019) analyze the gelation behavior of waxy crude oil emulsion, and establish its dynamic model to provide a certain reference for waxy crude oil production and gathering and transportation optimization. For large-scale oil-gas gathering system,  establish a high-dimensional mixed integer nonlinear optimization mathematical model including pipeline design parameters, and propose a modified MPSO algorithm with global search ability for optimization solution. Aiming at the increase of water content in the old areas, the existing gathering and transportation modes cannot meet the requirements of energy saving in oilfield production, Zhao et al. (2019) derive the temperature drop and pressure drop models of oil, gas and water three-phase flow, and determine the operating parameters of single well gathering, branch connection gathering and annular gathering for four wells.
The above optimization methods for oil-gas production process have achieved remarkable results in practice, but they often only consider a single objective and lack of the comprehensive consideration of oil-gas production process. In the actual production process, oil recovery enterprises may have different production requirements during different periods, i.e., at the initial stage of production, the main goal is to maximize the oil production within a specified range of comprehensive energy consumption. At the middle and later stage, the main goal is to minimize the comprehensive energy consumption within a specified range of oil production. Additionally, the requirements of the overall oil production during different periods are not the same. So, it is necessary to consider some production indices simultaneously to satisfy the different requirements of the overall oil production, the energy consumption and so on. Ideally, these indices are expected to be the best at the same time. However, due to the conflict among these indicators, it is impossible for all the indices to reach the optimal values simultaneously. For example, when an oil well has a certain amount of oil production, it needs more input power of pumping system for continuing to increase oil production, which will consume more energy consumption of oil-gas gathering and transportation system. Ultimately, the more oil output tends to consume more energy. Additionally, when the output changes in a certain range, the increase of output will also lead to the increase of water production. Therefore, the improper allocation of these indices may lead to poor product quality, difficult to ensure the continuity and stability of production, and ultimately affect the economic benefits of oil recovery enterprises. In addition, due to the complexity of oil-gas production process, there may be a large error between the prediction value and the actual value of objective function. Therefore, it is necessary to establish a more accurate objective function model in the comprehensive optimization of oilgas production process.
To address the above problems, this paper establishes an error compensation multiobjective optimization model for the integrated oil-gas production process. In order to improve the prediction accuracy of objective function model, a LS-SVM error compensation method optimized by Bayesian optimization algorithm (BOA) is put forward to compensate the mechanism model error. Then, an improved non-dominated sorting genetic algorithm-II with multi-strategy improvement (IMS-NSGA-II) is proposed to solve the multi-objective optimization model. Finally, the oil-gas production process of one block in an oil recovery operating area is taken as the research background. The simulation results show the effectiveness and superiority of the established optimization model and algorithm.

Oil-gas production process
In this paper, rod pumping wells are taken as production wells, and the flowchart of the oilgas production process is shown in Figure 1. From Figure 1, it can be seen that a typical integrated oil-gas production system consists of reservoirs, production wells, manifolds, surface flow lines and surface facilities (Jahanshahi et al., 2020). First, the mixture of oil, gas and water is transported from the reservoir to each oil well, then to the oil well platform (manifold), and then to the surface facilities through the surface flow lines for oil, gas and water separation, processing and recycling. In the entire production system, the oil wells are the important component of production, and also the main production equipment to transport oil, gas and water multiphase flow from reservoir to the surface. In actual production, the oil production, water production and comprehensive energy consumption are usually important production indices, in which the comprehensive energy consumption mainly includes the power consumption in sucker-rod pumping system, the power consumption and the natural gas consumption in gathering system.
Multi-objective optimization model of oil-gas production process based on error compensation of comprehensive energy consumption This section is organized as follows: The first part is the objectives of the optimization model for the integrated oil-gas production process: including overall oil production, overall water production, comprehensive energy consumption per ton of oil. Then, the second part describes the multi-objective optimization model. Finally, the error compensation of comprehensive energy consumption per ton of oil is given in the third part: including hybrid model structure based on error compensation, LS-SVM error model based on Bayesian optimization algorithm.
Objectives of the optimization model for the integrated oil-gas production process Overall oil production. Overall oil production is the sum of the oil production of each well in the block platform, which is expressed as follows: where Q 0 is the overall oil production (m 3 /d), Q j is the actual flow rate (m 3 /d) of the j -th oil well, w 0j is the water content (%), and N is the number of rod-pumped wells. For each oil well, the wellbore model consists of two parts, including the multiphase flow from the reservoir into the wellbore and the flow rate in the pipe, they are represented by the inflow performance relation (IPR) curve of oil wells and outflow characteristic curve of pump respectively. When the static pressure is lower than the saturation pressure, the inflow performance relationship of the reservoir is expressed by Vogel equation (2) (Mccoy and Podio, 2001).
(2) where q max is the maximum liquid production ( p wf ¼ 0 ) of oil well (m 3 /d), p wf is the bottom hole flow pressure, p r is the reservoir pressure (Pa), p sat is defined as the saturation pressure (Pa). When p r > p sat , equation (2) can be expressed as follows (Patton and Goland, 1980): The inflow performance of oil well reflects the oil supply capacity of the reservoir to the oil well, which depends on reservoir pressure, geological properties of well surroundings, and the wellbore itself, etc.
Due to the oil flow resistance between reservoirs and well bores, the pressure of oil production well bore changes slowly with time with time, and the following equation can be obtained (Horne, 1998).
where B, l, k c , h, c t , r B and / are experimentally parameters related to the reservoir and the oil well. On the other hand, for rod pumping wells, the actual flow rate of the j -th well is calculated as follows: In equation (5), D j ; S j ; n j ; a j represent the pump diameter (m), the stroke (m), the stroke times (min À1 ) and the pump efficiency (%) of the j-th rod-pumped well, respectively. g S ; g F ; g L ; g V represent the effective stroke coefficient of oil well, fullness coefficient, leakage coefficient, and volume coefficient of gas dissolution crude oil under submerging pressure, which can be calculated as follows (Dong et al., 2008a;Dong et al., 2008b): In the above formula, S pmup is the effective stroke length of plunger (m), K is the clearance coefficient, R denotes the gas liquid ratio of pump suction (m 3 /m 3 ), p s and p d represent submerging pressure (Pa) and discharge pressure (Pa) shown in equation (7). A p is the cross section area of plunger pump (m 2 ), k is the index of natural gas polytropic process, w 0 is the water content (%), DQ is the liquid clearance leakage between the plunger and pump barrel (m 3 ). B ops , B wps are the oil volume coefficient and the water volume coefficient under the pump suction, which have a certain relationship between p d and reservoir temperature t L .
where q w is the water density (t/m 3 ), q 0 is the crude oil density (t/m 3 ), p c is the casing pressure (MPa), p o is the oil pressure (MPa) and H d is the working fluid level.
When the oil well is producing under stable working conditions, the flow rate from the reservoir into the wellbore is equal to the actual production flow rate of the oil well, which is to meet the oil well production coordination conditions as follows.
In addition, the relation between the bottom hole flow pressure and the submerging pressure is given as follows: where L Z is the central reservoir depth (m), L i is the pump depth (m) of the j -th rod-pumped well, q j is the density of oil water mixture (t/m 3 ). So far, equation (9) can be substituted into equation (8) to determine the nonlinear relationship between flow pressure, submerging pressure, reservoir characteristics, crude oil physical properties and suction parameters. When other relevant production parameters are determined and measured, the flowing pressure and the submerging pressure can be determined according to different strokes and stroke times. Therefore, the actual oil flow rate Q j can be obtained, and finally the overall oil production can be calculated as shown in equation (1).
Overall water production. Overall water production of the block is expressed as follows: where Q w is the overall water production (m 3 /d).
Comprehensive energy consumption per ton of oil. The comprehensive energy consumption per ton of oil is the energy consumption for producing one ton of crude oil, described as follows: where y is the comprehensive energy consumption per ton of oil (kg/t), Q DW , R and R 1 are the low calorific value of fuel (kJ/m 3 ), the electrical energy conversion coefficient (kJ/kW. h) and the equal value to the heat of per kilogram standard coal (kJ/kg), respectively. W 1 represents the power consumption in the sucker-rod pumping system (kW. h/t), W 2 and Q g represent the power consumption (kW.h/t) and the natural gas consumption (m 3 /t) in the gathering system, respectively. W 1 ; W 2 and Q g can be calculated by the following equations.
where W 1j ¼ P j t is the electricity consumption of the j-th oil well (kW. h), P j is the input power of suction system (kW), t is the operation time of the corresponding oil well in one day (h), P j is a function of suction parameters when the management level and the technical equipment are invariant (Dong et al., 2010). Q p and Q T are the total pressure loss (kJ/h) and heat loss (kJ/h) of gathering pipelines, respectively. And they can be calculated by the related equations (Dong et al., 2017;. When reservoir characteristics, crude oil physical property, technical equipment and management level are certain, Q T and Q p can be considered as the function of Q j , wellhead heating temperature T wj ( w) and pump head H j (m), g T and g p are the equipment efficiency.

Multi-objective optimization model
According to references (Dong et al., 2008a;Dong et al., 2008b;Dong et al., 2017;He et al., 2007;Zheng and Deng, 2007), the decision variables of the optimization model are selected as follows: where S ¼ S 1 ; S 2 ; Á Á Á; S N ½ 5 n 1 ; n 2 ; Á Á Á; n N ½ , T w 5 T w1 ; T w2 ; Á Á Á; T wN ½ and 5 H 1 ; H 2 ; Á Á Á; H N ½ represent the stroke, the stroke times, the wellhead heating temperature and the pump head vector of the rod-pumped wells, respectively. N is the number of rod-pumped wells. So, a multi-objective optimization model of oil-gas production process is established, which is composed of the decision variables in equation (13) and the objective functions described by equations (1), (10), (11), under the corresponding constraints, it is formulated as follows: where S jmin and S jmax are the minimum and maximum values of the j -th oil well stroke. n jmin and n jmax are the minimum and maximum values of the stroke times. T wmin and T wmax are the minimum and maximum values of heating temperature, respectively. H min and H max are the minimum and maximum limits of pump head vector, respectively. T j is the temperature of the oil entering the oil transfer station (production platform) from the oil gathering pipeline and T oil is the maximum allowable value of T j . T o ( C) is the temperature of oil flowing from station and T omax ( C) is the maximum allowable value of T o . p qi (Pa) is the wellhead back pressure and ½p o is the upper limit of p qi (Pa), p z (Pa) is the pressure of oil getting into station and ½p 1 (Pa) is the minimum allowable value of p z .
Error compensation of comprehensive energy consumption per ton of oil Hybrid model structure based on error compensation. According to the optimization model in equation (14), the expression of objective function y is complex. Because the mechanism model of the whole oil-gas production process is established under some assumptions and simplifications, there may be a large error between the output of the mechanism model and the actual value. In order to reduce the error, we propose a hybrid model based on the mechanism model and the LS-SVM error model, its structure is shown in Figure 2. In this hybrid model, the mechanism model is used to describe the overall characteristics of oil-gas production process, and then a LS-SVM error model based on Bayesian optimization algorithm is established to compensate the model error which can't be described by the mechanism model, the output of hybrid model is as follows: where y c is the prediction output of hybrid model, y is the output of the mechanism model for the comprehensive energy consumption per ton of oil, andê is the prediction error between the output of the mechanism model and the actual value. Through the analysis of oil-gas production process, the input vector X is composed of the stroke and the stroke times of each oil well, oil heating temperature and pump head constitute, and the comprehensive energy consumption per ton of oil is taken as the output of the model.
LS-SVM error model based on Bayesian optimization algorithm. For the convenience of description, we select fx i ; e i g L i¼1 as the sample set, where x i ¼ X, e i ¼ y Actual À y, L is the number of samples. The optimization problem of LS-SVM error model is described as follows : where w is the weight vector, c is the regularization parameter (also known as the penalty factor), 1 i is the error variables, b is the bias term. To solve this optimization problem, Lagrange function is constructed as follows: where a i are Lagrange multipliers. The solution of equation (17) can be determined by partial differentiation to w, b, 1 i and a i : After removing the variables w and 1 i , the equation is obtained as follows: where e ¼ ½e 1 ; e 2 ; Á Á Á; e L T , Z ¼ 1; 1; Á Á Á; 1 ½ , a ¼ a 1 ; a 2 ; Á Á Á; a L ½ T and X ij ¼ Kðx i ; x j Þ . Because there are various forms of kernel function Kðx; x i Þ, the most widely used Gaussian function is selected in this paper.
where r is kernel parameter, it denotes the width of the RBF which controls the regression ability. So far, the final expression of error estimation model of comprehensive energy consumption error is obtained as follows:ê During the training process of LS-SVM error model, the selection of penalty factor c and kernel function parameter r are important to improve the accuracy and generalization ability of the model. In order to reduce the influence of uncertainty in the actual production process on the model, and further improve the accuracy and stability of the model, a Bayesian optimization algorithm (BOA) (Shahriari et al., 2016) is used to obtain the optimal values of c and r by maximizing the posterior probability of the hyperparameters. Bayesian inference is constructed with three levels of inference (Xie et al., 2019): 1. First level of inference: Let H as the model space, D as the sample training set, for a given hyperparameter k ¼ 1=c, if the training samples are independent and identically distributed, and the prior distribution pðwjk; HÞ of the parameter w obeys the normal distribution, then the posterior probability of w can be obtained as follows: where e i is the regression error of each training sample, the optimal value w best is obtained by maximizing equation (22).

Second level of inference: Let
where A ¼ kr 2 E 1 þ r 2 E 2 , E best 1 and E best 2 are the values of E 1 and E 2 at w best , respectively. r is the vector differential operator, and then the optimal value k best is obtained by maximizing lnpðkjD; HÞ .
3. Third level of inference: Assuming that the prior probability of all sample data pðHÞ is a flat distribution, and then the posterior probability pðHjDÞ of parameter r can be obtained as follows pðHjDÞ / pðDjk best ; HÞ r Then, the optimal kernel parameter r best can be calculated by maximizing lnpðHjDÞ . A multi-objective optimization model based on error compensation of comprehensive energy consumption is constructed by replacing y in equation (14) with the output of hybrid model y c in equation (15). Other objective function models can also be compensated by the same method, and then a multi-objective optimization model based on error compensation is built, that is error compensation multi-objective optimization model.
Optimization solution of error compensation multi-objective optimization model using IMS-NSGA-II algorithm NSGA-II algorithm with multi-strategy improvement (IMS-NSGA-II) Through the above analysis, the multi-objective optimization model of oil-gas production process is a typical complex multi-objective optimization problem with multiple variables and multiple constraints. At present, the main methods to solve this kind of multi-objective optimization problems are multi-objective evolutionary algorithms (MOEAs), including dominant relation-based MOEAs, index-based MOEAs and decomposition-based MOEAs. The basic idea of dominant relation-based MOEAs is to find all non-dominant individuals from the current population by using Pareto-based fitness allocation strategy, such as non-dominated sorting genetic algorithm-II (NSGA-II) (Deb et al., 2002), Strength Pareto Evolutionary Algorithm2(SPEA2) (Zitzler et al., 2001) and Pareto Envelop-based Selection Algorithm (PESA-II) (Corne et al., 2001). NSGA-II is one of the most popular evolutionary algorithms, it is an improved version of "NSGA", its main feature is the fastnon-dominated sorting for ranking solutions at its selection stage, and an effective elite strategy is adopted. Compared with previous evolutionary algorithms, NSGA-II has lower computational complexity and better convergence (Noori-Darvish and Tavakkoli-Moghaddam, 2012). NSGA-II has been successfully applied to many practical complex projects (Maity et al., 2018;Wang et al., 2017;Yang et al., 2018). In recent years, many scholars have introduced the differential evolution algorithm (Rainer and Price, 1997) into the field of multi-objective problems, forming a multi-objective differential evolution (MODE). Some scholars have also integrated DE algorithm into NSGA-II, using crossover and mutation operators instead of simulated binary crossover operators to improve the convergence performance. However, many traditional methods do not consider how to maintain the diversity of Pareto optimal solutions, and are prone to converge prematurely, when the control parameters are not set properly, the algorithm is easy to fall into local optimum in the late evolution period (Guo et al., 2012). In addition, the traditional DE algorithm cannot use the gradient information to further explore the local optimal solution after obtaining the optimal solution of a certain generation of population. Therefore, through the above analysis, based on the NSGA algorithm, this paper introduces the chaos idea to avoid falling into local optimum, and uses the gradient information near the optimal solution to further improve the convergence ability of the algorithm, and proposes an improved NSGA-II algorithm with multi-strategy improvement (IMS-NSGA-II). The specific ideas are as follows: In IMS-NSGA-II, a chaotic sequence is first introduced to initialize the population to maintain the diversity of the initial population. Then, the gradient operator is introduced into the crossover operator and mutation operator to improve the local search speed of NSGA-II algorithm. Finally, during the evolutionary process, when an improved population is selected, the cyclic crowding distance is adopted to control the distribution of Pareto non-dominated solutions, so as to further improve the diversity and uniformity of the Pareto non-dominated solutions. Through these improved operations, IMS-NSGA-II algorithm not only guarantees the uniformity of initial population and the diversity of Pareto non-dominated solutions, but also improves the convergence of the algorithm to a great extent.
Population initialization based on chaotic sequence. Population initialization plays an important role in the accuracy and diversity of NSGA-II algorithm. In the standard NSGA-II algorithm, the uniformly distributed random number generator is usually used for crossover and mutation operations, but there are still some shortcomings in maintaining horizontal diversity and obtaining the Pareto-front with high uniformity (Guo et al., 2010). Chaotic motion has the characteristics of regularity, randomness and ergodicity. This process is not convergent, but bounded, and it is extremely sensitive to external parameters and initial values. Therefore, the population initialization based on chaotic sequence can traverse the search space within a certain range according to its own law, which significantly improves the accuracy and diversity of the solution obtained by the algorithm (Chen et al., 2014;Lu et al., 2013). Therefore, chaotic initialization is first introduced to enrich the diversity of the initial population, thus laying the foundation for global optimization.
In this paper, the chaotic sequence generated by Skew Tent mapping is used for population initialization (Hasler and Maistrenko, 1997;Long et al., 2019), and its mathematical model is expressed as follows: When u 2 (0,1) and x k 2 [0,1], equation (25) is in chaotic state. In equation (25), if x k is the j -dimensional variable, the generated chaotic variables are mapped to the decision variables, and the j -th decision variable x ji of the i -th initialized individual X 0 i ¼ ðx 0 1i ; Á Á Á ; x 0 ji ; Á Á Á ; x 0 ni Þ is obtained as follows: where x jmax and x jmin are the upper and lower limits of the j -th decision variable, respectively. m is the dimension of decision variable vector, N P is the population size.
The gradient-based operator. As a random search algorithm, the convergence speed of evolutionary algorithm cannot be ignored. Since the negative gradient direction is the fastest direction to reduce the value of the objective function, in order to improve the optimization speed, we introduce the gradient operator into NSGA-II algorithm, and combine it with the crossover and mutation operator to form the comprehensive operator (Yu et al., 2011), and then use the comprehensive operator to generate a new generation of population, expressed as follows: where x Ã i ! 0; i ¼ 1; 2; 3, equation (27) means that the new solutions are generated by the crossover operator, the mutation operator and the gradient-based operator with the probability x Ã 1 , x Ã 2 and x Ã 3 respectively. For gradient operator, in practical engineering problems, it is often impossible to get the analytical expressions directly, so it is necessary to estimate the gradient of function. The standard approach of gradient estimation is the finite difference (FD) method, i.e., fðxÞ ¼ ½f 1 ðxÞ; f 2 ðxÞ; Á Á Á; f m ðxÞ T ; x 2 R n represents a multi-objective function, x ¼ ðx 1 ; x 2 ; . . . ; x n Þ is a n -dimensional decision vector, the one-sided FD method used for the j -th objective gradient estimation requires ( n þ 1) function evaluations, as follows: where e i represents the unit vector in the i -th direction, and c is the step size. FD method is usually expensive in function evaluation. In contrast, the calculation time of Simultaneous Perturbation (SP) method is less (Zhou, 2017), so the SP method is used in the paper to estimate the gradient of the j -th objective function, and the i -th component of gradient is expressed as follows: where D is a n -dimensional vector of random perturbations, and D i is the i -th component of D . For each component of D, a simple choice is used to have a Bernoulli distribution AE1 and the probability of each AE1 outcome is 0.5. The step size c at the k th iteration (denoted by c k ) is given as follows: c k ¼ c 0 =ðk þ 1Þ c , the values of c 0 and c are set to 0.001 and 1/6, respectively. Then the search direction can be obtained by using the gradient operator in equation (30).
where the search direction e is the normalized vector of d, and the convex cone combination of negative gradient direction of each objective function, which makes each selected point moves to the Pareto front along a common descent direction of some objectives, k j ðj ¼ 1; 2; Á Á Á; mÞ 2 (0,1). Obviously, the direction d in equation (30) has a higher chance to minimize all the objectives (or at least one objective) than the variation direction in the whole space by the pure random genetic operator, which can also reduce the invalid trial times of crossover and mutation, and improve the local search ability.
After the e is determined, a new solution x kþ1 can be created by the gradient-based operator.
where t k is the step size of the k -th generation.
Circular crowded sorting strategy. In order to further improve the diversity of Pareto nondominated solutions and overcome the shortcomings of traditional crowded sorting, and to make the algorithm can jump out of local optimum and prevent premature convergence when solving complex multi-objective optimization problems. In this paper, a circulating crowded ranking algorithm is adopted (Luo et al., 2010). The specific steps are as follows: • Step 1: Calculate the crowding distance of each individual (solution), and delete the individual with the smallest crowding distance in the population • Step 2: Among the remaining individuals, recalculate the crowding distance, and continue to delete the individuals with the smallest crowding distance in the population.

Optimization solution process of error compensation multi-objective optimization model using IMS-NSGA-II algorithm
The design variables in the multi-objective optimization model of oil-gas production process include discrete variables and continuous variables. Here, we use binary encoding for discrete variables S and real encoding for continuous variables n, T w and H . The flow chart of proposed IMS-NSGA-II algorithm for solving the multi-objective optimization problem of oil-gas production process is shown in Figure 3.
The specific steps of optimization solution of error compensation multi-objective optimization model using IMS-NSGA-II algorithm are as follows: • Step 1: The measured data in actual production and the data generated by mechanism model of the comprehensive energy consumption per ton of oil are collected to establish the hybrid model of comprehensive energy consumption per ton of oil. • Step 2: Hybrid model of comprehensive energy consumption per ton of oil is introduced into the multi-objective optimization model of oil-gas production process, and then the multi-objective optimization model based on error compensation is built, that is the error compensation multi-objective optimization model.

Simulation results and analysis
Firstly, the production process with six oil wells in an oil recovery operation area is taken as background to verify the effectiveness of the established optimization model and the proposed IMS-NSGA-II algorithm. The basic data of oil well pumping equipment is shown in Table 1, where pumping unit model CYJ10-3-53HB and motor model Y250M-6 are used for oil wells numbers 11-22, 11-17, 11-25 and 11-20, and pumping unit model CYJ14-5-73HB and motor model YQ280M2-8 are used for wells numbers 11-23 and 11-27. Besides the basic data of the oil well, some production data, such as reservoir characteristic parameters, wellhead parameters and physical characteristic parameters of crude oil and natural gas, are also needed for model verification (Dong et al., 2010). According to the actual situation, some production parameters are fitted or selected as follows: B ¼1.21 m 3 /m 3 , l ¼10 mPa.s, k c ¼20 md, h ¼23 m, c t ¼8.32Âl0 À6 psi, r B ¼0.5m, / ¼0.21. oil reservoir pressure p r ¼15.3 Mpa, reservoir temperature t L ¼50.23 C, saturation pressure of crude oil p sat ¼10.5 Mpa, water density q w ¼1 t/m 3 , crude oil density q 0 ¼0.863 t/m 3 , crude oil viscosity l 0 ¼2.1 mPaÁs, water viscosity l w ¼1.005 mPaÁs, natural gas relative density d gs ¼0.75, natural gas polytropic process index k ¼1.1. For the maximum liquid production of oil well q max , under the condition of known reservoir pressure and saturation pressure, if a certain liquid production of oil well q 0 and its corresponding flow pressure p wf are measured, they are substituted into equation (2) or equation (3) to obtain the corresponding maximum liquid production of the oil well.  In oil well production, some production data such as oil pressure p o , casing pressure p c and the water content, are not constant and can be measured every day. Table 2 shows the corresponding production data of each oil well in the block on a certain day. As the water content will not change greatly in a short term and only fluctuate in a small range, the water content can be regarded as a constant in the simulation experiment.
According to the basic parameters of pumping units, oil pipelines and the actual production requirements, the parameters in constraints are set as follows: The stroke constraints of rod pumped well are s jmax ¼ 6 m and s jmin ¼ 1 m, respectively. The stroke times constraints of rod pumped well are n jmax ¼ 12 min À1 and n jmin ¼ 2 min À1 , respectively. The allowable temperature of crude oil getting into station T oil ¼ 35 C and the allowable temperature of crude oil out of station T omax ¼ 80 c, the pressure of oil gathering pipeline getting into the oil transfer station and the processing station is 0.2 MPa.
Parameters of IMS-NSGA-II are set as follows: the population size N p ¼ 100, the maximum evolutionary generation G max ¼ 150, distribution indexes for real encoding crossover operator and mutation operator of mu ¼ 20 and mum ¼ 20, crossover probability x Ã 1 ¼ 0:8, mutation probability x Ã 2 ¼ 1=n ( n is the dimension of decision vector), the probability of gradient-based operator x Ã 3 ¼ 1, u ¼ 0.5, the step size t k ¼ 0.1(gen-k þ 1)/gen, gen ¼ 100.
Firstly, 200 groups of data are collected form the actual production process, and the data are generated by the mechanism model of comprehensive energy consumption. Among them, 150 groups are taken as training samples and 50 groups are taken as testing samples. During training LS-SVM error model, the Bayesian optimization is applied to obtain the optimal values of ( c, r ), the values of c and r are determined to be 31.26 and 25.37, respectively. Then, the trained LS-SVM error model is used to compensate the error of mechanism model, and the hybrid model is compared with the mechanism model, the comparison results are shown in Figure 4. The performance comparison of the two models is shown in Table 3.
From Table 3, it can be seen that the average absolute error (MAE), root mean square error (RMSE) and mean absolute percentage error (MAPE) of the hybrid model are smaller than those of the mechanism model. IMS-NSGA-II algorithm is used to solve the multi-objective optimization models before and after error compensation. In order to further validate the reliability of optimization results, the manipulated variables are adjusted according to the optimization results, and the corresponding actual values of objectives are obtained. Part of Pareto non-dominated solutions are listed in Table 4, where y represents the comprehensive energy consumption per ton of oil. From the comparison results in Table 4, it can be seen that the corresponding actual values obtained by solving error compensation multi-objective optimization model are better than that those obtained by solving multi-objective optimization model before error compensation, such as solution 1, the actual value (40.7745, 135.1560, 74.3573) is better than the value (40.1320, 135.9883, 77.6981). It shows that the method not only improves the overall oil production rate, but also reduces the comprehensive energy consumption and the overall water production to a certain extent. Moreover, due to the higher prediction accuracy of the error compensation multi-objective optimization model, the optimization results obtained by the multi-objective optimization model based on the error compensation are closer to the actual values than that obtained by the uncompensated model, which can reduce the impact of uncertainty and greatly improves the reliability of optimization results.
Then, in order to verify the superiority of the proposed algorithm, IMS-NSGA-II algorithm is tested on the standard test problems including ZDT1 and ZDT2 (Zitzler et al., 2000), and compared with other algorithms, such as the original NSGA-II algorithm, SPEA2 and hybrid multi-objective intelligent algorithm (C-NSGA-II) (Yan et al., 2018) for multi-objective optimization.
The parameters of these algorithms are set as follows: For these algorithms, the population size is 100, the dimension of decision variables of test problems ZDT1 and ZDT2 is n ¼ 30. For IMS-NSGA-II algorithm, the parameters are set as above, and the other operations are consistent with those used in the relevant literatures (Deb et al., 2002;Yan et al., 2018;Zitzler et al., 2001).
In this paper, the hypervolume indicator I H (Zitzler et al., 2000) and spatial evaluation method  are applied to evaluate the performance of the algorithm. Hypervolume indicator I H is defined as follows: I H ðAÞ of a solution set A X can be where kðSÞ is the Lebesgue measure of a set S and ½f 1 ðaÞ; r 1 Â . . . Â ½f m ðaÞ; r m is the m -dimensional hypercuboid composed of all points weakly dominated by the point a rather than by the reference point. When the value of I H ðAÞ is higher, the convergence is better. These algorithms are run on each test problem for 10 times, respectively, and each run selects a different initial population. The maximum number of function evaluations of all algorithms is 10000. The optimization results of a random run on ZDT1 and ZDT2 are shown in Figures 5 and 6, respectively. The corresponding hypervolume indicator I H and D of these algorithms are listed in Tables 5 and 6, respectively. In order to show the significant difference between these algorithms, we also give the optimization results in Figures 5 and 6, they are obtained by these algorithms that the number of function evaluations is 5000.
From Figures 5 and 6 and Tables 5 and 6, it can be seen that the performance of IMS-NSGA-II on the problems: ZDT1 and ZDT2 is better than NSGA-II, SPEA2 and slightly better than C-NSGA-II in terms of the hypervolume indicator I H and the diversity metric D . The results also indicate that the Pareto non-dominated solutions obtained by IMS-NSGA-II have better convergence and diversity than other algorithms when these solutions are distributed along the real Pareto-optimal set. In addition, when solving these problems, it not only doesn't significantly increase the running time of these gradient-based operators, but also reduces the average running time used by IMS-NSGA-II. Finally, both NSGA-II algorithm and IMS-NSGA-II algorithm are used to solve the error compensation multi-objective optimization model. The two algorithms have been run  10 times on the multi-objective optimization problem that the number of evolutionary generations is 100. The optimization results of a random run are shown in Figure 7. The corresponding hypervolume indicator I H and D of the two algorithms are listed in Tables 7  and 8, respectively. It can be seen from Figure 7 that the Pareto non-dominated solution set obtained by IMS-NSGA-II algorithm are superior to that obtained by the traditional NSGA-II algorithm in terms of convergence and diversity. So, it is considered that IMS-NSGA-II algorithm can provide more superior solutions for decision-making departments in oil production enterprise, and the decision makers can choose more reasonable oil production programs to satisfy the different requirements of actual production.
From Table 7, it can be seen that each I H obtained by IMS-NSGA-II algorithm is larger than that obtained by NSGA-II algorithm. So, it further proves that the Pareto nondominated solutions obtained by the IMS-NSGA-II algorithm is closer to the real Pareto front and has better convergence. In Table 8, it can be seen that D of the Pareto  non-dominated solution set obtained by IMS-NSGA-II algorithm is smaller, which demonstrates that IMS-NSGA-II algorithm has a better diversity and uniformity. In addition, it can be seen from Tables 7 and 8 that the maximum value, the minimum value and the average value of hypervolume indicator I H obtained by IMS-NSGA-II algorithm for 10 runs are 0.5137, 0.4896 and 0.4990, respectively. The maximum value, the minimum value and the average value of D obtained by IMS-NSGA-II algorithm are 19.9750, 17.2326 and 18.30825, respectively. And for each run, the values of I H and D obtained by IMS-NSGA-II algorithm are better than that obtained by NSGA-II algorithm. Thus, all the results indicate that the improved algorithm is better than NSGA-II algorithm and has good robustness.

Conclusions
In this paper, we propose an error compensation multi-objective optimization model for the integrated oil-gas production process. Then an improved NSGA-II (IMS-NSGA-II) algorithm is presented to solve the optimization model. When IMS-NSGA-II algorithm is test on the problems: ZD1 and ZD2, the experimental results indicate the effectiveness and superiority of the proposed IMS-NSGA-II algorithm. Besides, the application results on oilgas multi-objective optimization demonstrate that the proposed optimization method can achieve better production indices. It not only satisfies the requirements of the actual oil-gas production process, but also reduce the impact of uncertainty and improves the reliability of the optimization results. Thus, the proposed optimization method is a more effective tool for optimization and control of oil-gas production process. These simulation results also imply that the proposed optimization method has a good application protest in dealing with the optimization problems of other similar complex production processes.