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

Abstract

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.

Keywords

Oil-gas production process LS-SVM error model Bayesian optimization algorithm error compensation NSGA-II with multi-strategy improvement (IMS-NSGA-II)

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, Liu et al. (1999) 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 oil-gas production process.

To address the above problems, this paper establishes an error compensation multi-objective 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 oil-gas production process is shown in Figure 1.

Figure 1.

Flowchart of integrated oil-gas production system.

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:

Q_{0} = \sum_{j = 1}^{N} Q_{j} (1 - w_{0 j})

(1)

where

Q_{0}

is the overall oil production (m³/d),

Q_{j}

is the actual flow rate (m³/d) of the

j

-th oil well,

w_{0 j}

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).

q^{I N} = q_{\max} {1 - 0.2 (\frac{p^{w f}}{p_{r}^{}}) - 0.8 {(\frac{p^{w f}}{p_{r}^{}})}^{2}}

(2)

where

q_{\max}

is the maximum liquid production (

p^{w f} = 0

) of oil well (m³/d),

p^{w f}

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):

q^{I N} = {\begin{cases} \frac{p_{sat}}{p_{r}^{} - p_{sat}} (p_{r}^{} - p^{w f}) p^{w f} \geq p_{sat} \\ p_{sat} + (q_{\max} - p_{sat}) {1 - 0.2 (\frac{p^{w f}}{p_{sat}}) - 0.8 {(\frac{p^{w f}}{p_{sat}})}^{2}} p^{w f} < p_{sat} \end{cases}

(3)

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).

p^{w f} - p_{r}^{} = \frac{141.2 q^{I N} B μ}{k_{c} h} (\frac{1}{2} [\ln \frac{0.00264 k t}{ϕ μ c_{t} r_{B}^{2}} + 0.80907])

(4)

where

B

μ

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:

\begin{array}{l} Q_{j} = 1440 π D_{j}^{2} S_{j} n_{j} α_{j} / 4 \\ α_{j} = η_{i S} η_{i F} η_{i L} η_{i V} \end{array}

(5)

In equation (5), $D_{j}, S_{j}, n_{j}, α_{j}$ represent the pump diameter (m), the stroke (m), the stroke times (min⁻¹) and the pump efficiency (%) of the $j$ -th rod-pumped well, respectively. $η_{S}, η_{F}, η_{L}, η_{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):

{\begin{cases} η_{S} = S_{pmup} / S \\ η_{F} = \frac{1}{R} {1 - \frac{K R}{1 + R {(p_{s} / p_{d})}^{1 / k}} [1 - {(p_{s} / p_{d})}^{1 / k}]} \\ η_{L} = \frac{A_{p} S η_{S} η_{F} η_{V} - Δ Q}{A_{p} S η_{S} η_{F} η_{V}} \\ η_{V} = 1 / [(1 - w_{0}) B_{ops} + w_{0} B_{wps}] \end{cases}

(6)

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³/m³), $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²), $k$ is the index of natural gas polytropic process, $w_{0}$ is the water content (%), $Δ Q$ is the liquid clearance leakage between the plunger and pump barrel (m³). $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}^{}$ .

{\begin{cases} p_{d} = p_{0} + [w_{0} ρ_{w} + (1 - w_{0}) ρ_{0}] g L \\ p_{s} = p_{c} + [w_{0} ρ_{w} + (1 - w_{0}) ρ_{0}] g (L - H_{d}) \end{cases}

(7)

where

ρ_{w}

is the water density (t/m³),

ρ_{0}

is the crude oil density (t/m³), 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.

q_{j}^{I N} = Q_{j}

(8)

In addition, the relation between the bottom hole flow pressure and the submerging pressure is given as follows:

p_{j}^{w f} = p_{s j}^{} + (L_{Z} - L_{j}) ρ_{j} g

(9)

where

L_{Z}

is the central reservoir depth (m),

L_{i}

is the pump depth (m) of the

j

-th rod-pumped well,

ρ_{j}

is the density of oil water mixture (t/m³). 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:

Q_{w} = \sum_{j = 1}^{N} Q_{j} w_{0 j}

(10)

where

Q_{w}

is the overall water production (m³/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:

y = [Q_{g} Q_{D W} + (W_{1} + W_{2}) R] / R_{1}

(11)

where

y

is the comprehensive energy consumption per ton of oil (kg/t),

Q_{D W}

R

and

R_{1}

are the low calorific value of fuel (kJ/m³), 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³/t) in the gathering system, respectively.

W_{1}, W_{2}

and

Q_{g}

can be calculated by the following equations.

{\begin{cases} W_{1} = \sum_{i = 1}^{N} W_{1 i} / Q_{0} \\ W_{2} = Q_{p} t / (η_{p} R Q_{0}) \\ Q_{g} = Q_{T} t / (η_{T} Q_{D W} Q_{0}) \end{cases}

(12)

where

W_{1 j} = 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; Liu et al., 1999). 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_{w j}

(°w) and pump head

H_{j}

(m),

η_{T}

and

η_{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; Liu et al., 1999; Zheng and Deng, 2007), the decision variables of the optimization model are selected as follows:

X = {S, n, T_{w}, H}

(13)

where

S = [S_{1}, S_{2}, \cdot \cdot \cdot, S_{N}]

n= [n_{1}, n_{2}, \cdot \cdot \cdot, n_{N}]

T_{w} = [T_{w 1}, T_{w}_{2}, \cdot \cdot \cdot, T_{w N}]

and

H= [H_{1}, H_{2}, \cdot \cdot \cdot, 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:

\begin{array}{l} {\max Q_{0}, \min Q_{w}, \min y} \\ s . t . {\begin{cases} S_{j}_{\min} \leq S_{j} \leq S_{j}_{\max} \\ n_{j}_{\min} \leq n_{j} \leq n_{j}_{\max} \\ T_{w \min} \leq T_{w j} \leq T_{w \max} \\ H_{\min} \leq H_{j} \leq H_{\max} \\ T_{j} \geq T_{oil}, T_{o} \leq T_{o \max} \\ p_{q i} \leq [p_{o}], p_{z} \geq [p_{1}] \end{cases} \end{array}

(14)

where

S_{j \min}

and

S_{j \max}

are the minimum and maximum values of the

j

-th oil well stroke.

n_{j \min}

and

n_{j \max}

are the minimum and maximum values of the stroke times.

T_{w \min}

and

T_{w \max}

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_{o \max}

(°C) is the maximum allowable value of

T_{o}

p_{q i}

(Pa) is the wellhead back pressure and

[p_{o}]

is the upper limit of

p_{q i}

(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.

Figure 2.

Structure of hybrid model.

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:

y_{c} = y + \hat{e}

(15)

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

\hat{e}

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 ${x_{i}, e_{i}}_{i = 1}^{L}$ 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 (Chen et al., 1999):

\min \underset{w, b, ς}{J (w, ς)} = \frac{1}{2} w^{T} w + \frac{1}{2} γ \sum_{i = 1}^{L} ς_{i}^{2} s . t . e_{i} = 〈 w, φ (x_{i}) 〉 + b + ς_{i} i = 1, 2, \cdot \cdot \cdot, L

(16)

where

w

is the weight vector,

γ

is the regularization parameter (also known as the penalty factor),

ς_{i}

is the error variables,

b

is the bias term. To solve this optimization problem, Lagrange function is constructed as follows:

L (w, b, ς; α) = J (w, b, ς) - \sum_{i = 1}^{L} α_{i} {〈 w, φ (x_{i}) 〉 + b + ς_{i} - e_{i}}

(17)

where

α_{i}

are Lagrange multipliers. The solution of equation (17) can be determined by partial differentiation to

w

b

ς_{i}

and

α_{i}

{\begin{cases} \frac{\partial L}{\partial w} = 0 \to w = \sum_{i = 1}^{L} α_{i} φ (x_{i}) \\ \frac{\partial L}{\partial b} = 0 \to \sum_{i = 1}^{L} α_{i} = 0 \\ \frac{\partial L}{\partial ς_{i}} = 0 \to α_{i} = γ ς_{i}, i = 1, 2, \cdot \cdot \cdot, L \\ \frac{\partial L}{\partial α_{i}} = 0 \to 〈 w, φ (x_{i}) 〉 + b + ς_{i} - e_{i} = 0, i = 1, 2, \cdot \cdot \cdot, L \end{cases}

(18)

After removing the variables $w$ and $ς_{i}$ , the equation is obtained as follows:

[\begin{array}{l} 0 Z^{T} \\ Z Ω + 1 / γ \end{array}] [\begin{array}{l} b \\ α \end{array}] = [\begin{array}{l} 0 \\ e \end{array}]

(19)

where

e = {[e_{1}, e_{2}, \cdot \cdot \cdot, e_{L}]}^{T}

Z = [1, 1, \cdot \cdot \cdot, 1]

α = {[α_{1}, α_{2}, \cdot \cdot \cdot, α_{L}]}^{T}

and

Ω_{i j} = 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.

K (x, x_{i}) = \exp (- {‖ x - x_{i} ‖}^{2} / 2 σ^{2})

(20)

where

σ

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:

\hat{e} = \sum_{i = 1}^{L} α_{i} K (x, x_{i}) + b

(21)

During the training process of LS-SVM error model, the selection of penalty factor $γ$ and kernel function parameter $σ$ 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 $γ$ and $σ$ by maximizing the posterior probability of the hyperparameters. Bayesian inference is constructed with three levels of inference (Xie et al., 2019):

First level of inference: Let H as the model space, D as the sample training set, for a given hyperparameter $λ = 1 / γ$ , if the training samples are independent and identically distributed, and the prior distribution $p (w | λ, H)$ of the parameter $w$ obeys the normal distribution, then the posterior probability of $w$ can be obtained as follows:

p (w| D, λ, H) \propto \exp (- \frac{1}{2} w^{T} w - \sum_{i = 1}^{L} e_{i}^{2})

(22)

where $e_{i}^{}$ is the regression error of each training sample, the optimal value $w_{best}$ is obtained by maximizing equation (22).

2. Second level of inference: Let $E_{1} = \frac{1}{2} w^{T} w$ , $E_{2} = \sum_{i}^{L} e_{i}^{2}$ , the posterior probability of hyperparameter λ is shown as follows:

p (λ | D • H) \propto \frac{{(2 π λ)}^{k / 2} \exp (- λ E_{1}^{best} - E_{2}^{best})}{{(\det A)}^{1 / 2}}

(23)

where $A = λ \nabla^{2} E_{1}^{} + \nabla^{2} E_{2}^{}$ , $E_{1}^{best}$ and $E_{2}^{best}$ are the values of $E_{1}^{}$ and $E_{2}^{}$ at $w_{best}$ , respectively. $\nabla$ is the vector differential operator, and then the optimal value $λ_{best}$ is obtained by maximizing $\ln p (λ | D, 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 (H | D)$ of parameter $σ$ can be obtained as follows

p (H | D) \propto \frac{p (D | λ_{best}, H)}{σ}

(24)

Then, the optimal kernel parameter $σ_{best}$ can be calculated by maximizing $\ln p (H | D)$ .

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 fast-non-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:

{\begin{cases} x_{k + 1} = x_{k} / φ, 0 \leq x_{k} \leq φ \\ x_{k + 1} = (1 - x_{k}) / (1 - φ), φ < r_{k} \leq 1 \end{cases}

(25)

When $φ$ ∈ (0,1) and $x_{k}$ ∈ [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_{j i}^{}$ of the $i$ -th initialized individual $X_{i}^{0} = (x_{1 i}^{0}, \dots, x_{j i}^{0}, \dots, x_{n i}^{0})$ is obtained as follows:

\begin{array}{l} x_{j i}^{} = x_{j \min} + x_{j k} (x_{j \max} - x_{j \min}) \\ j = 1, 2, \dots, m, i = 1, 2, \dots, N_{P} \end{array}

(26)

where

x_{j \max}

and

x_{j \min}

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:

ω_{1}^{*} o perato r_{C}^{*} \oplus ω_{2}^{*} o perato r_{M}^{*} \oplus ω_{3}^{*} o perato r_{G}^{*}

(27)

where

ω_{i}^{*} \geq 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

ω_{1}^{*}

ω_{2}^{*}

and

ω_{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), \cdot \cdot \cdot, f_{m} (x)]}^{T}, x \in R^{n}

represents a multi-objective function,

x = (x_{1}, x_{2}, \dots, 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:

g_{j i} (x) = \frac{f_{j} (x + c e_{i}) - f_{j} (x)}{c}

(28)

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:

g_{j i} (x) = \frac{f_{j} (x + c Δ) - f_{j} (x)}{c Δ_{i}}

(29)

where

Δ

is a

n

-dimensional vector of random perturbations, and

Δ_{i}

is the

i

-th component of

Δ

. For each component of Δ, a simple choice is used to have a Bernoulli distribution ±1 and the probability of each ±1 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)}^{γ}

, the values of

c_{0}

and

γ

are set to 0.001 and 1/6, respectively. Then the search direction can be obtained by using the gradient operator in equation (30).

{\begin{cases} \bar{e} = d / ‖ d ‖ \\ d = - \sum_{j = 1}^{m} λ_{j} \frac{g_{j} (x)}{‖ g_{j} (x) ‖} \end{cases}

(30)

where the search direction

\bar{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,

λ_{j} (j = 1, 2, \cdot \cdot \cdot, m)

∈ (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 $\bar{e}$ is determined, a new solution $x_{k + 1}$ can be created by the gradient-based operator.

x_{k + 1} = x_{k} + t_{k} \bar{e}

(31)

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 non-dominated 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.

Step 3: Determine the number of remaining solutions $N'$ . If $N' < N_{p}$ , return to Step 2, otherwise, skip to Step 4.

Step 4: End the selection.

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.

Figure 3.

Flow chart of optimization solution using IMS-NSGA-II algorithm.

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.

Step 3: IMS-NSGA-II algorithm is put forward for solving the error compensation multi-objective optimization model established in Step 2.

Step 4: The values of manipulated variables can be adjusted according to the optimization results obtained in Step 3.

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.

Table 1.

Basic data of rod pumping equipment.

Oil well number	Central reservoir depth (m)	Casing diameter (mm)	Oil tube diameter (mm)	Pump diameter (mm)	Sucker rod string Lr (m)×dr(mm)	Pump depth (m)
11-22	2167.3	127	89	44	22 × 436.4 + 19 × 1363.0	1899.4
11-17	2198.3	139.7	73	56	22 × 703.6 + 19 × 1247.16	1950.76
11-25	2236.8	124	62	38	22 × 795.9 + 19 × 1354.19	2150.09
11-23	2245.2	127	89	44	22 × 710.1 + 19 × 1390.62	2100.72
11-27	2183.39	127	89	44	22 × 725.6 + 19 × 1257.63	1983.23
11-20	2175.5	124	62	38	22 × 746.8 + 19 × 1253.32	2000.12

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³/m³, $μ$ =10 mPa.s, $k_{c}$ =20 md, $h$ =23 m, $c_{t}$ =8.32×l0⁻⁶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 $ρ_{w}$ =1 t/m³, crude oil density $ρ_{0}$ =0.863 t/m³, crude oil viscosity $μ_{0}$ =2.1 mPa·s, water viscosity $μ_{w}$ =1.005 mPa·s, natural gas relative density $δ_{g s}$ =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'$ and its corresponding flow pressure $p^{w f}$ 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.

Table 2.

Part of production data for oil wells.

Oil well number	casing pressure p_c (MPa)	oil pressure p_o (MPa)	Water content (%)
11-22	2.3	0.3	69
11-17	2.5	0.25	76.4
11-25	5.9	0.4	76
11-23	8.3	0.3	82
11-27	2.9	0.3	73
11-20	3.8	0.2	74.5

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_{j \max} = 6$ m and $s_{j \min} = 1$ m, respectively. The stroke times constraints of rod pumped well are $n_{j \max} = 12$ min⁻¹ and $n_{j \min} = 2$ min⁻¹, 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_{o \max} = 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 $m u = 20$ and $mum = 20$ , crossover probability $ω_{1}^{*} = 0. 8$ , mutation probability $ω_{2}^{*} = 1 / n$ ( $n$ is the dimension of decision vector), the probability of gradient-based operator $ω_{3}^{*} = 1$ , φ = 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 ( $γ$ , $σ$ ), the values of $γ$ and $σ$ 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.

Figure 4.

Comparison of prediction results of comprehensive energy consumption per ton of oil.

Table 3.

Performance comparison of the two models.

Model	MAE	RMSE	MAPE
mechanism model	1.2561	2.0615	1.9202
hybrid model	0.5108	0.8147	0.7926

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.

Table 4.

Comparison of optimization results obtained by the optimization models before and after error compensation.

Pareto non-dominated solutions	Multi-objective optimization	$Q_{0}$ (m³/d)	$Q_{w}$ (m³/d)	$y$ (kg/t)
1	Before the error compensation	38.5341	134.2653	74.6722
	Actual value	40.1320	135.9883	77.6981
	After the error compensation	40.1914	134.3697	73.5061
	Actual value	40.7745	135.1560	74.3573
2	Before the error compensation	35.3816	115.6281	67.2615
	Actual value	33.9282	113.6112	70.2185
	After the error compensation	33.7540	111.5943	67. 6088
	Actual value	34.3298	112.0596	68.3626
3	Before the error compensation	28.6029	90.7451	66.7500
	Actual value	26.2670	88.7538	64.0897
	After the error compensation	27.2575	89.5439	62.9282
	Actual value	26.7489	88.8115	63.3762
4	Before the error compensation	19.4133	66.0509	58.9031
	Actual value	17.5297	68.1072	56.7124
	After the error compensation	18.1671	65.8618	57.3297
	Actual value	18.9598	66.4186	56.6079
5	Before the error compensation	11.7958	43.3517	53.4215
	Actual value	10.1160	46.0211	56.4906
	After the error compensation	11.2061	45.7571	53.7122
	Actual value	10.4965	45.1350	54.3072

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 (Deb et al., 2000) 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 \subseteq X$ can be defined as the hypervolume of the space, which is dominated by the set $A$ and bounded by a reference point $r = (r_{1}, \dots, r_{m}) \in R^{m}$ , expressed as follows:

I_{H} (A) = λ (\underset{a \in A}{\cup} [f_{1} (a), r_{1}] \times \dots \times [f_{m} (a), r_{m}])

(32)

where

λ (S)

is the Lebesgue measure of a set

S

and

[f_{1} (a), r_{1}] \times \dots \times [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 $Δ$ 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.

Figure 5.

Pareto non-dominated solutions obtained by NSGA-II, SPEA2, C-NSGA-II and IMS-NSGA-II on ZDT1.

Figure 6.

Pareto non-dominated solutions obtained by NSGA-II, SPEA2, C-NSGA-II and IMS-NSGA-II on ZDT2.

Table 5.

Convergence metric, diversity metric and average running time on ZDT1.

Algorithm	Average $I_{H}$	Average $Δ$	Average running time (min)
Algorithm	Evaluations =10000	Evaluations =10000	Evaluations = 10000
NSGA-II	0.80215	0.42064	0.04819
SPEA2	0.80329	0.79345	0.34566
C-NSGA-II	0.81650	0.23498	0.03843
IMS-NSGA-II	0.82703	0.21523	0.02975

Table 6.

The convergence metric, the diversity metric and average running time on ZDT2.

Algorithm	Average $I_{H}$	Average $Δ$	Average running time (min)
Algorithm	Evaluations =10000	Evaluations =10000	Evaluations = 10000
NSGA-II	0.50256	0.49253	0.04896
SPEA2	0.50327	0.78326	0.36425
C-NSGA-II	0.52212	0.30126	0.04237
IMS-NSGA-II	0.53018	0.27250	0.03049

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 $Δ$ . 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 $Δ$ of the two algorithms are listed in Tables 7 and 8, respectively.

Figure 7.

Distribution chart of Pareto set obtained by NSGA-II and IMS-NSGA-II algorithm.

Table 7.

Hypervolume indicator $I_{H}$ obtained by the two optimization algorithms.

Running time	Hypervolume indicator $I_{H}$
Running time	NSGA-II	IMS-NSGA-II
1	0.4825	0.5016
2	0.4510	0.4903
3	0.4839	0.4981
4	0.4576	0.5137
5	0.4178	0.5092
6	0.3995	0.4896
7	0.4655	0.4958
8	0.4960	0.5062
9	0.4582	0.4934
10	0.4729	0.4919

Table 8.

Deviation metrics $Δ$ obtained by the two optimization algorithms.

Running time	Diversity measure $Δ$
Running time	NSGA-II	IMS-NSGA-II
1	27.7852	18.0334
2	25.0656	17.9582
3	22.8245	18.0861
4	24.2369	18.7307
5	25.6253	19.9750
6	28.7549	18.1925
7	25.9347	17.6280
8	24.1972	17.2326
9	25.3425	19.0859
10	29.2590	18.1601

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 non-dominated 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 $Δ$ 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 $Δ$ 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 $Δ$ 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 oil-gas 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.

Footnotes

Acknowledgments

The author(s) would like to thank the editor and the anonymous reviewers for their helpful comments and high-quality suggestions, we also thank for the support from the Liaohe Oilfield of China, providing us research and experimental conditions.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by National Natural Science Foundation of China (32001415, 31901400, 61673281, 61903264), China Postdoctoral Science Foundation (2018M631820), Liaoning Province Doctoral Research Initiation Fund Guidance Project (2019-BS-207), Liaoning province natural fund guidance plan project (2019-ZD-0720).

ORCID iD

Tan Liu

References

Al-Fehdly

ElMaraghy

Wilkinson

(2019) Carbon footprint estimation for oil production: Iraq case study for the utilization of waste gas in generating electricity. Procedia CIRP 80: 389–392.

Chebouba

Yalaoui

Smati

, et al. (2009) Optimization of natural gas pipeline transportation using ant colony optimization. Computers & Operations Research 36(6): 1916–1923.

Chen

Wen

Gong

, et al. (1999) Realization of robot motion based on a learning method. IEEE Transactions on Automation Control 44: 271–376.

Chen

Yuan

, et al. (2014) Design of a fractional order PID controller for hydraulic turbine regulating system using chaotic non-dominated sorting genetic algorithm II. Energy Conversion and Management 84: 390–404.

Cheng

Song

, et al. (2020) Study on weak link of energy utilization in oil transfer station system: Insights from energy level analysis method. Energy Reports 6: 1097–1105.

Codas

Campos

Camponogara

, et al. (2012) Integrated production optimization of oil fields with pressure and routing constraints: The Urucu field. Computers & Chemical Engineering 46: 178–189.

Corne

Jerram

Knowles

, et al. (2001) PESA-II: region-based selection inevolutionary multiobjective optimization. In: Spector

Goodman

, et al. (eds) Proceedings of the genetic and evolutionary computation conference (GECCO’2001). San Francisco, CA: Morgan Kaufmann Publishers, pp.283–290.

Cui

Wang

Zhang

, et al. (2008) Optimization of the operation efficiency of oil-gas gathering pipeline network based on genetic algorithm. Journal of Xi’an Shiyou University (Natural Science Edition) 23(6): 59–62.

Deb

Agrawal

Pratap

, et al. (2000) A fast-elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II. In: 6th International conference on parallel problem solving from nature, Paris, France. 18-20 September 2000, pp.849–858. Germany: Springer Verlag.

10.

Deb

Pratap

Agarwal

, et al. (2002) A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation 6(2): 182–197.

11.

Dong

Yan

, et al. (2008a) Simulating maximum of system efficiency of rod pumping wells. Journal of System Simulation 20: 3533–3537.

12.

Dong

Wang

Sun

, et al. (2017) Oil gathering parameters optimization method based on the lowest energy consumption of machine mining and oil gathering system. Journal of China University of Petroleum 41(5): 123–129.

13.

Dong

Wang

, et al. (2008b) Simulation models for optimization design of suction parameters for rod pumping system in directional wells. Acta Petrolei Sinica 29: 120–123.

14.

Dong

Zhang

, et al. (2010) Simulation and optimization method for the integral energy-saving suction parameters of rod pumping wells in an oil field. Acta Petrolei Sinica 31: 475–479.

15.

Foss

(2012) Process control in conventional oil and gas fields – Challenges and opportunities. Control Engineering Practice 20(10): 1058–1064.

16.

Guo

Wang

P Huang

, et al. (2010) Chaotic-NSGA-II: An effective algorithm to solve multi-objective optimization problems. In: IEEE international conference on intelligent computing and integrated systems (ICISS), Guilin, China, 22-24 October 2010, pp.20–23. USA: IEEE.

17.

Guo

Gui

, and Yang

(2012) An improved hybrid differential evolution algorithm used for multi-objective optimization of aluminum electrolysis. Journal of Central South University (Science and Technology) 43(1): 185–188.

18.

Hasler

Maistrenko

(1997) An introduction to the synchronization of chaotic systems: Coupled skew tent maps. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications 44(10): 856–866.

19.

Zhu

(2007) Gathering and transferring parameter optimization of 401 transfer station of Xingliu joint station. Journal of Daqing Petroleum Institute 31: 34–37.

20.

Horne

(1998) Modern Well Test Analysis. 2nd ed. Palo Also: Petroway Inc.

21.

Jahanshahi

Krishnamoorthy

És Codas

, et al. (2020) Plantwide control of an oil production network. Computers and Chemical Engineering 136: 1–12.

22.

Jun

Kang

Jeong

Lee

(2017) An efficient approach for optimizing full field development plan using Monte-Carlo simulation coupled with genetic algorithm and new variable setting method for well placement applied to gas condensate field in Vietnam. Energy Exploration & Exploitation 35(1): 75–102.

23.

Kosmidis

Perkins

Pistikopoulos

(2005) A strategy for simulation and optimization of gas and oil production. Computers & Chemical Engineering 29(7): 1523–1541.

24.

Lei

Weng

Luo

, et al. (2019) Achievements and future work of oil and gas production engineering of CNPC. Petroleum Exploration and Development 46(1): 145–152.

25.

Liu

Chen

(1999) Optimal parameters design of oilfield surface pipeline systems using fuzzy models. Information Sciences 120(1–4): 13–21.

26.

Liu

Zhao

Wang

(1999) An optimizing design and analysis of oil gas gathering of ring type mixing hot water. Acta Petrolei Sinica 20: 77–81.

27.

Long

Cai

Jiao

, et al. (2019) An improved grey wolf optimization algorithm. Acta Electronica Sinica 47(1): 169–175.

28.

Niu

Liu

, et al. (2013) A chaotic non-dominated sorting genetic algorithm for the multi-objective automatic test task scheduling problem. Applied Soft Computing 13(5): 2790–2802.

29.

Luo

Chen

Zhang

(2010) Improved NSGA-II algorithm with circular crowded sorting. Control and Decision 25(2): 227–231.

30.

Mccoy

Podio

(2001) Rotaflex efficiency and balancing. In: Proceedings of SPE production and operations symposium, Oklahoma City, Oklahoma, USA, 24-27 March 2001, pp.453–461. USA: Society of Petroleum Engineers.

31.

Maity

Roy

Maiti

(2018) A rough multi-objective genetic algorithm for uncertain constrained multi-objective solid travelling salesman problem. Granular Computing 46: 1–18.

32.

Noori-Darvish

Tavakkoli-Moghaddam

(2012) Minimizing the total tardiness and make span in an open shop scheduling problem with sequence-dependent setup times. Journal of Industrial Engineering International 8(1): 1–13.

33.

Patton

Goland

(1980) Generalized IPR curves for predicting well behavior. Petroleum Engineering Institute 52(7): 74–82.

34.

Rainer

Price

(1997) Differential evolution – A simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization 11(4): 341–359.

35.

Rui

Wang

Zhang

, et al. (2018) A realistic and integrated model for evaluating oil sands development with steam assisted gravity drainage technology in Canada. Applied Energy 213: 76–91.

36.

Shahriari B, Swersky K, Wang ZY, et al. (2016) Taking the Human Out of the Loop: A Review of Bayesian Optimization. Proceedings of the IEEE 104(1): 149–175.

37.

Tan

Xiang

Zhao

, et al. (2016) Energy consumption prediction and application in oil and gas gathering and transferring system production based on large data. Acta Petrolei Sinica 37(S2): 158–164.

38.

Wang

Huang

, et al. (2017) A NSGA-II based memetic algorithm for multiobjective parallel flowshop scheduling problem. Computers & Industrial Engineering 113: 185–194.

39.

Wang

Dong

Gan

, et al. (2018) Method for improving system efficiency by dynamic parameters simulation model of low-production pumping wells and the ways to improve system efficiency. Acta Petrolei Sinica 39(11): 1299–1307.

40.

Wang

Bai

Zhang

, et al. (2019) Investigation on gelation nucleation kinetics of waxy crude oil emulsions by their thermal behavior. Journal of Petroleum Science and Engineering 181: 106230.

41.

Jia

, et al. (2014) Optimal operation of trunk natural gas pipelines via an inertia-adaptive particle swarm optimization algorithm. Journal of Natural Gas Science and Engineering 21: 10–18.

42.

(2011) A study on optimal operation and energy saving of super heavy oil pipeline. In: International conference on computer distributed control and intelligent environmental monitoring (CDCIEM), Changsha, China, 19-20 February 2011, pp.400–403. USA: IEEE.

43.

Xie

Wang

Liu

, et al. (2019) Nonlinear joint PP-PS AVO inversion based on improved Bayesian inference and LSSVM. Applied Geophysics 16(1): 64–76.

44.

Yan

Zhang

, et al. (2018) The optimal operation strategy for hybrid combined cooling, heating and power microgrid with compressed air energy storage. Proceedings of the CSEE 40: 6925–6935.

45.

Yang

Cao

and (2018) Optimizing the impact of economic restructuring on Shanghai’s carbon emissions using STIRPAT and NSGA-II. Sustainable Cities and Society 40: 44–53.

46.

Chai

Luo

(2011) Multi-objective production planning optimization using hybrid evolutionary algorithms for mineral processing. IEEE Transactions on Evolutionary Computation 15(4): 487–513.

47.

Zhao

Yao

Liu

, et al. (2019) Optimization of energy saving gathering and transportation mode in a block of Tahe oilfield. Case Studies in Thermal Engineering 13: 100378–100376.

48.

Zheng

Deng

(2007) Research and application on designing method of sucker-rod pumping system with the least energy consumption. Acta Petrolei Sinica 28: 129–132.

49.

Zhou

Hou

, et al. (2017) Estimation of relative permeability curves using an improved Levenberg-Marquardt method with simultaneous perturbation Jacobian approximation. Journal of Hydrology 544: 604–612.

50.

Zitzler

Laumanns

Thiele

, et al. (2001) SPEA2: Improving the strength Pareto evolutionary algorithm. TIK-report 103: 1–22.

51.

Zitzler

Thiele

Deb

(2000) Comparison of multi-objective evolutionary algorithms: Empirical results. Evolutionary Computation 8(2): 173–195.