An integrated MILP method for gathering pipeline networks considering hydraulic characteristics

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

Contents lists available at ScienceDirect

Chemical Engineering Research and Design journal homepage: www.elsevier.com/locate/cherd

An integrated MILP method for gathering pipeline networks considering hydraulic characteristics Bingyuan Hong a , Xiaoping Li a , Guojia Di a , Yu Li a,b , Xingshuai Liu c , Shilin Chen d , Jing Gong a,∗ a

National Engineering Laboratory for Pipeline Safety/MOE Key Laboratory of Petroleum Engineering/Beijing Key Laboratory of Urban Oil and Gas Distribution Technology, China University of Petroleum-Beijing, Fuxue Road No. 18, Changping District, Beijing, 102249, PR China b China National Oil & Gas Exploration and Development Co. Ltd, Bldg D, Guoji Invest Plaza, Fuchengmen North Street, Xicheng District, Beijing, 100034, PR China c China Petroleum Engineering & Construction Corp Beijing design branch, Xinxi Road No. 8, Shangdi, Haidian District, Beijing, 100085, PR China d CNOOC China Limited, Unconventional Oil & Gas Branch/China United Coalbed Methane Corporation Ltd, No.21B Jiuxianqiao Road, Chaoyang District, Beijing, 100015, PR China

a r t i c l e

i n f o

a b s t r a c t

Article history:

High investment highlights the importance of optimizing gathering pipeline networks. In

Received 9 December 2018

view of the problem of ignoring hydraulic characteristics in previous studies, an integrated

Received in revised form 13 August

mixed-integer linear programming (MILP) model for gathering pipeline networks is devel-

2019

oped in this study. Taking the minimum total construction cost as the objective function,

Accepted 19 August 2019

the proposed model considers the economic and technical constraints such as obstacles,

Available online 11 October 2019

three-dimensional terrain, pipeline topological structures, pipe diameters, wellhead back pressure, and pressure equipment. The ant colony optimization algorithm is used for route

Keywords:

optimization to provide parameters for the proposed model. Moreover, a piecewise method

Natural gas

is employed to linearize the nonlinear hydraulic equations. Then, the model is solved by

Gathering pipeline networks

the CPLEX solver to obtain the optimal solution integrally, including the optimal connection

Optimization

topology, the location of the central processing facility, the position of pressure equipment,

Topological structure

the diameter and detailed route of each pipeline, the pipeline flowrate and the node pressure.

Hydraulic characteristics

Finally, two real-world gas fields and a virtual oil field are taken as examples to demonstrate the feasibility and practicality of the model. The optimal results illustrate that this model

MILP

can be implemented as a decision-support tool to optimize gathering pipeline networks in the actual design process. © 2019 Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

1.

Introduction

1.1.

Background

With the economic development and population growth, global energy demand is expected to increase by 30% in 2035 (Kan et al., 2019). Oil and natural gas, as important energy

∗

resources, play a vital role in the development of the modern global economy (Yu et al., 2018b). The global proved reserves of natural gas increased by 21.76% between 2000 and 2008, and further increased by 9.99% in 2008–2016. At the same time, oil, coal, and natural gas production increased by 9.85%, 6.27%, and 16.29%, respectively, between 2008 and 2016, and the consumption increments were 9.76%, 5.77%,

Corresponding author. E-mail address: [email protected] (J. Gong). https://doi.org/10.1016/j.cherd.2019.08.013 0263-8762/© 2019 Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

Nomenclature Abbreviations ACO Ant colony optimization Cascade dendritic pipeline network. CDPN Chinese Yuan. CNY Central processing facility. CPF IDPN Insertion dendritic pipeline network. Linear programming. LP MILP Mixed-integer linear programming model. MINLP Mixed-integer nonlinear programming. Indexes and sets Set of the numbering of flowrate segments in A the piecewise method, denoted by the index a. Set of pipeline diameters, denoted by the D indices d. I, J Set of coordinate of all nodes in the studied area, denoted by the indices i and j. Set of coordinate of well nodes in the studied IW, JW area, IW ε I, JW ε J. Set of potential supercharging devices, denoted V by the indices v. K Set of direction of connecting nodes, denoted by the indices k, rk is the opposite direction of what k stands for. Parameters bwi,j Determining parameter of well sites. If node (i, j) is well site bwi,j = 1, otherwise, bwi,j = 0. Construction cost of CPF. cc cpd Construction cost of pipeline per unit length with diameter d. Construction cost of supercharging device type csv v. disi,j,k Connecting distance from node (i, j) to the adjacent node at its direction k. Flow coefficient of pressure drop equation for fd pipeline with diameter d. Parameter of pressure term in pressure drop m1 equation. m2 Parameter of flowrate term in pressure drop equation. obi,j Determining parameter of obstacles. If there are obstacles in node (i, j), obi,j = 1 and 0 otherwise. qpmaxd Upper bound of economic flowrate of pipeline with diameter d. qpmind Lower bound of economic flowrate of pipeline with diameter d. qsmaxv Upper bound of flowrate of supercharging device type v. qsminv Lower bound of flowrate of supercharging device type v. qtotal Total production of all well sites. Output of well site (i, j). qwi,j = 0 if (i, j) is not a qwi,j well site. qxa Value of flowrate at break points by piecewise linear approximation. qya Value of the m2th power to flowrate at break points by piecewise linear approximation. ppmax Maximum pipeline network operating pressure.

ppmin pcpfmax pcpfmin pwi,j psv

321

Minimum pipeline network operating pressure. Maximum requirement pressure of CPF. Minimum requirement pressure of CPF. Maximum wellhead pressure of well site (i, j). pwi,j = 0 if (i, j) is not a well site. Pressure provided by supercharging device v.

Binary decision variables Binary variable of CPF. If the node (i, j) is the BCi,j central processing facility in the gathering network, BCi,j = 1 and 0 otherwise. Binary variable of ordinary node. If the node (i, BNi,j j) is neither a well site nor a CPF, BNi,j = 1 and 0 otherwise. BPi,j,k,d Binary variable of pipeline construction. If a new gathering pipeline is constructed from node (i, j) to the adjacent node at its direction k with diameter d, BPi,j,k,d = 1 and 0 otherwise. BQAi,j,k,a Binary variable of the flowrate segment. If the flowrate delivered in the pipeline from node (i, j) to the adjacent node at its direction k in segment a for pressure drop equation, BQAi,j,k,a = 1 and 0 otherwise. Binary variable of supercharging device. If there BSi,j,v will be built a supercharging device type v in node (i, j), BSi,j,v = 1 and 0 otherwise. Positive decision variables F Total cost. The construction cost of CPF. FC The construction cost of pipelines. FP FS The construction cost of supercharging devices. Value of the m1th power to pressure of node (i, Pi,j m1 j). QPi,j,k Flowrate from node (i, j) to the adjacent node at its direction k. QPAi,j,k Approximate flowrate value from node (i, j) to the adjacent node at its direction k for pressure drop equation. Receiving flowrate at node (i, j). QRi,j QWi,j,k,a weights of break points by piecewise linear approximation.

and 16.36%, respectively (Chen et al., 2018). Although China’s natural-gas consumption was only 29.184 × 109 m3 in 2002, the consumption has increased with an average annual growth rate of 15.35% from 2002 to 2014, reaching 186.894 × 109 m3 and accounting for 5.4% worldwide in 2014 (Ding, 2018). The natural gas consumption of China in 2016 was 210.3 × 109 m3 , with 7.7% growth rate per annum (Wang et al., 2018b). The significant increase in oil and natural gas utilization necessitates the planning, design and construction of production fields to serve the ever-increasing demand (Crow et al., 2018). The gathering pipeline network is utilized to transport oil and natural gas during the production process in oil-gas fields and it is composed of well sites, pipelines, pressure equipment and central processing facilities (CPF) (He et al., 2018; Wei et al., 2016). There are two common connection structures of gathering pipeline networks, namely cascade dendritic pipeline network (CDPN) and insertion dendritic pipeline network (IDPN) (Zhang et al., 2017b; Zhou et al., 2014). For CDPN, as shown in Fig. 1(a), single wells are usually connected with the nearest well site.

322

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

Fig. 1 – Schematic diagram of gathering pipeline network. The well sites are connected to other adjacent well site, and the trunk lines are gathered to CPF in radial structure. Such a connection structure is advantageous for reducing the length and construction cost of the pipeline and it is applicable for the rolling development in the fields with larger areas, more wells and divided ground due to geographical conditions. For IDPN, as shown in Fig. 1(b), the gathering trunk line runs through the main production area of the oil-gas field. The flow line of each well site is not connected to each other like CDPN, but to be connected with the nearby trunk line in a shortest vertical distance. The connection mode is applicable for fields with evenly distributed wells and no need for supercharging. Generally, the gathering pipeline network is characterized by a massive and complex connection structure, which makes its construction cost account for a large proportion of the total investment of the entire oil-gas field (Rios-Mercado and Borraz-Sanchez, 2015). The high proportion of cost highlights the importance of optimizing the topological structure and gathering process parameters of the gathering pipeline network, which substantially affects the efficiency and profit of oil-gas field production. Traditionally, the design of gathering pipeline network has been manually performed by experienced engineers, taking various factors such as terrain conditions, environmental conservation, economic and technical constraints into account. However, since innumerable existences of connection types between all well sites, this is an iterative and very complex process highly dependent on the expertise of engineers, and it takes a long time to obtain a satisfactory result (de Lucena et al., 2014; Marcoulaki et al., 2012; Zhou et al., 2018). Therefore, it has been recognized that there is an urgent need for optimization techniques that not only reduces investment but also accelerates engineering design procedures. The purpose of this work is to establish an optimal design method of gathering pipeline network, which can be used to design a gathering network with high transport efficiency and low construction cost. This optimization problem contains several sub-problems, including the division of well groups (Wang et al., 2014, 2012), the topological structure of gathering pipeline networks (Hong et al., 2018; Rosa et al., 2018; Zhang et al., 2017b; Zhou et al., 2014, 2018), the location of facilities (Marcoulaki et al., 2012), the detailed shortest route of pipelines (Baeza et al., 2017; Baioco et al., 2017; de Lucena et al., 2014; Wang et al., 2018a; Zhou et al., 2017) and others (Jonsson and Magnusdottir, 2017).

Some of those problems are typical non-deterministic polynomial complete problems and often nonlinear coupled, which means that it is difficult to reach the optimal results in polynomial time (Alves et al., 2016; Mezura-Montes and Coello Coello, 2011). Due to the various pipe network structures, the three-dimensional terrain, the existence of obstacles, and the hydraulic characteristics of the pipeline network, these optimization problems are very complicated, massively computational, multidisciplinary and multi-objective.

1.2.

Related work

Many scholars have acknowledged the need for systematic methodologies to design optimal gathering pipeline networks and have carried out researches on different aspects of this issue (Rios-Mercado and Borraz-Sanchez, 2015). Gupta and Grossmann (2012) established a mixed-integer nonlinear programming (MINLP) model for optimal planning of offshore oil and gas field infrastructure, which obtained the maximizing total net present value and determined FPSO installation and connection with target fields. Marcoulaki et al. (2012) presented a nonlinear programming model for the routing and equipment design of main pipelines, while considering the process constraint and applying the stochastic optimization method to solve it. This model was extended to a broader space of pipeline configurations by adding to the effects of corrosion in their later work (Marcoulaki et al., 2014). Wang et al. (2012) proposed a mathematical model with the minimum connection cost between subsea wells and cluster manifolds to determine the layout of subsea production system and provided an efficient algorithm. In a follow-up study (Wang et al., 2014), the optimization of layout scenarios of cluster manifolds with pipeline end manifolds was discussed. Rodrigues et al. (2016) developed a 0–1 Linear Programming (LP) integrated model which sought to minimize the development costs of a given oil field, taking the whole system from the reservoir to the topside into consideration. Zhang et al. (2017b) came up with a Mixed-Integer Linear Programming (MILP) model to optimize the production well gathering system taking the minimum total investment as the objective function while considering terrain, obstacles, and other operational constraints. The optimal topological structure, position of the CPF, diameter and route of each pipeline were obtained through GUROBI solver. Zhang et al. (2017a) introduced a model that considered the gathering radius, economic flow

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

rate, terrain obstacles and production technique to optimize the offshore oilfield gathering system. Ramos Rosa et al. (Rosa et al., 2018) proposed a MILP model to design subsea production networks that accounts for the placement and number of manifold, the pipeline layout, the reservoir dynamics and multiphase flow. Hong et al. (2018) established an integrated model to optimize the layout of the subsea production system, taking seabed topography and obstacles into consideration, coupling the simulated annealing and Dijkstra algorithms to solve the model. In recent years, intelligent algorithms have been widely applied to address the design optimization of complex gathering networks. Lucena et al. (de Lucena et al., 2014) proposed to divide the design criteria of submarine pipeline routes into ¨ ¨ hard ¨ and presented a mathematic programming s¨ oftor ones, model to assess an optimal route. The different constraint handling methods were considered and the model was solved by genetic algorithm. Zhang (2015) investigated multiple connection modes of gathering networks and applied a genetic algorithm with hybrid encoding via a Hash table to solve the programming model for coalbed methane fields, obtaining the optimal topological connections. Dossary and Nasrabadi (Al Dossary and Nasrabadi, 2016) dealt with well placement optimization for maximum well productivity by using a metaheuristics algorithm known as the imperialist competitive algorithm, and compared the performance with that of the particle swarm optimization as well as the genetic algorithm. Baioco et al. (2017) introduced the advances in computational tools based on evolutionary algorithms for the synthesis and optimization of submarine pipeline routes and studied the influence of different on-bottom stability criteria for particular scenarios in shallow waters. Zhou et al. (2017) considered the route optimization of two-phase pipelines and adopted the general genetic algorithm and steady-state genetic algorithm to get the optimal route, respectively, with a focus on the influences of different parameters setting to the result. In the work of Baeza et al. (2017), two optimization approaches, namely ant colony optimization (ACO) method and Dijkstra algorithm, were compared in terms of their computational cost and accuracy for optimal ore concentrate pipeline routing in a complex terrain. Wang et al. (2018a) applied the A* algorithm to determine the route of natural gas pipeline networks. De Padua Agripa Sales et al. (de Padua Agripa Sales et al., 2018) presented a Monte Carlo simulation integrated with a genetic algorithm to obtain adequate solutions of the field layout design considering the uncertainties in well deliverability. A greedy algorithm was employed to evaluate the performance of the genetic algorithm. However, although the intelligence algorithms mentioned above can solve the optimization problem quickly, they could easily converge to the locally optimal solution as a result of the strong randomness (Fister et al., 2015; Xiang et al., 2014; Yang et al., 2007). Other studies of spatio-temporal energy system optimization such as CO2 , heating, water or hydrogen pipeline can also provide good references for the optimization of the gathering pipeline networks. Tian et al. (2017) proposed a robust optimization method for CO2 pipeline transportation design considering the multiple uncertainties and using the linear matrix inequalities. A stepwise method is presented to further improve the optimization performance. Kim et al. (2018) studied practical deployment of pipelines for the Carbon Capture and Storage (CCS) network in critical conditions using MINLP modelling and optimization. A rigorous cost model

323

was developed, and subdivided penalty factors were applied to consider the drastic change in geographical conditions. The effect of reservoir uncertainty and CCS policy options are analyzed. For urban energy systems, Morvaj et al. (2016) investigated the optimal design and operation of distributed energy systems, as well as optimal heating network layouts for different economic and environmental objectives. A MILP model was used for multi-objective optimization to minimize total cost and carbon emissions. Diogo et al. (2018) presented a model for optimal rehabilitation of separate sanitary sewer systems, taking into account network performance and wastewater treatment costs, average values of several input variables, and rates that can reflect the adoption of different predictable or limiting scenarios. Johnson and Ogden (2012) investigated long-term hydrogen pipeline planning and described a network optimization tool for identifying the lowest cost centralized production and pipeline transmission infrastructure within real geographic regions. In reality, field engineers often give the maximum back wellhead pressure and the pressure range required by the CPF (Drouven and Grossmann, 2016, 2017). Due to the frictional resistance of pipeline, the fluid pressure gradually decreases during the flow from well sites to the CPF (Kostowski and Skorek, 2012; Yu et al., 2018a). In order to ensure that the fluid can meet the pressure requirements of the CPF, it is essential to consider hydraulic characteristics in the design of the gathering pipeline network. The hydraulic characteristics are affected by the pipe diameter with different frictional resistance and the pressure devices installed at well sites (Abeysekera et al., 2016). Therefore, hydraulic constraints must be strictly considered in the optimization of gathering pipeline network. However, many of the previous studies have neglected the flow characteristics of pipelines or replaced the detailed hydraulic calculation with the gathering radius. As shown in Fig. 1, in the method of gathering radius, it is considered that when the distance from well site to CPF is less than the prescribed gathering radius, the whole gathering pipeline network satisfies the hydraulic and thermal constraints. This simplified method may lead to an imbalance of pipeline network pressure, especially under terrain relief conditions, such that the fluid cannot be successfully transported from wellhead to CPF, which is unrealistic in engineering.

1.3.

Contributions of this work

To address those issues, this study proposes a MILP optimization model for gathering pipeline networks with the minimum total construction cost as the objective function, considering two common gathering pipeline network connection forms, obstacles, three-dimensional terrain and other factors simultaneously. In addition, hydraulic constraints and pressure equipment constraints which were neglected by previous researches are taken into consideration to ensure that the optimal solutions meet the production requirements and are more in line with the actual situation. A piecewise method is employed to linearize the nonlinear equations in the model, thus the model can be solved by the MILP solver based on the branch-and-bound algorithm to obtain the global optimal solutions, including the optimal topology, the location of CPF, the location of pressure facilities, the size and route of each pipeline. Furthermore, two real-world gas fields in China and a virtual oil field are regarded as exam-

324

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

ples to verify the feasibility and practicality of the proposed model.

1.4.

Paper organization

The remainder of this paper is organized as follows. In Section 2, the problem is described, then the required data, model assumptions, the objective function, decision variables are stated. In Section 3, the development of the proposed optimization methodology for gathering pipeline networks is discussed. The details of the mathematical model and model solving methods are provided. Two real-world cases in China to illustrate our methodology are presented in Section 4. Finally, it is followed by conclusions and future work in Section 5.

2.

Problem description

In the construction of oil and gas fields, two kinds of boundary conditions must be observed. The first one is the data given by reservoir engineer, such as the well site position, the expected output and the maximum wellhead back pressure. The second one is that the pressure required of the CPF is usually within a limited range. The task of optimization design is to give an economical and safe construction scheme under these two kinds of boundary conditions, which can meet the hydraulic requirements with the lowest construction cost. The hydraulic characteristics refer to the relationship between fluid pressure and flowrate along the pipeline, obeying the law of conservation, including the continuity equation, momentum equation and energy equation (Pambour et al., ˘ 2014). The difference of physical 2016; Üster and Dilaveroglu, properties between oil and gas makes their hydraulic characteristics different. The oil is incompressible, so the fluid velocity does not change with the pressure when it flows through the diameter-determined pipe. However, due to the compressibility of the gas, the pressure drops during the flow, the density decreases, and the gas velocity increases. Therefore, when calculating the hydraulic characteristics of gas pipeline, the compressibility of the gas must be taken into account (Abeysekera et al., 2016). In addition to the hydraulic characteristics mentioned above, terrain and obstacles are also important considerations. In practice, many well sites are distributed over complex three-dimensional terrain, resulting in that the optimal distance between two nodes is not necessarily a straight line. Meanwhile, many existing buildings should be avoided when constructing pipelines. Hence, the routing between two nodes has to be optimized. We now state the planning problem described in this paper as follows. Given: The required information can roughly be divided into three parts.

(1) Geographic information: the detailed terrain of the study area, obstacle information and the location of well sites. (2) Economic parameters: pipe unit price (of different diameters), pressure equipment types (capacity and corresponding price of different classes), fixed cost of the CPF. (3) Technical parameters: the production rate and back pressure of each well site, the pressure required of the CPF.

Determine: (1) (2) (3) (4) (5)

The location of CPF. Connections and detailed routes of each pipeline. The pipeline diameter based on hydraulic characteristics. The installation location of pressure equipment. The flowrate of each pipeline and the operating pressure of each node. (6) The total construction cost. To enhance the solving efficiency, the following assumptions are made:

(1) This article focuses on two most common connection structures of gathering pipeline networks, namely IDPN and CDPN. Every well node can only choose one direction to deliver gas or oil to other nodes, except for the nodes confirmed as CPF. (2) There is only one CPF in gathering pipeline network and pressure equipment can only be installed in well site. (3) Ant colony optimization algorithm is used for minimizing the distance and without considering the influence of elevation changes. The elevation changes are also not taken into account in hydraulic calculation. (4) The hydraulic calculation in this paper refers to the single-phase steady-state model, without considering the variation of fluid parameters with time and multiphase flow.

3.

Methodology

The purpose of the optimization method for gathering pipeline networks is to establish a mathematical programming model and work out an optimal design scheme under given economic and technological parameters. As mentioned above, the optimal design is economically optimal under the premise of meeting hydraulic constraints, which means the pressure of the pipeline network should also be kept within allowable range. A methodology is proposed in this paper, which integrates the three-dimensional terrain, pipeline connection forms, obstacles, pressure equipment and hydraulic calculation into a MILP model simultaneously. The first step to apply our method is parameter input, including geographic parameters, economic parameters and technical parameters. Then, terrain parameters are gridprocessed and the corresponding parameters to the grid points are determined after generating a two-layer nonuniform grid. Subsequently, the ant colony optimization is used to determine the values of the parameters in the model that describe the pipeline lengths required to connect the nodes at the first layer coarse grid level. Furthermore, the piecewise approximation approach is employed to determine the values of the parameters in the model to eliminate the nonlinear term of hydraulic constraints. After parameter preparation, the proposed MILP model can be carried out to generate an optimal construction scheme of gathering pipeline network, which consists of the total construction cost, the integral topological structure, the location of CPF and pressure equipment, the detailed diameter and route of each pipeline, and the hydraulic state of the whole gathering pipeline network. The flow chart of the methodology is shown in Fig. 2 and the details of each step are illustrated in the following sections.

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

325

Fig. 2 – Flow chart of the methodology.

3.1.

Discrete grid division

According to the given location of the well site, we can obtain the corresponding geographic coordinates and elevation from Google Earth, import these data into MATLAB, and use cubic spline interpolation to derive the topographic data for the study area. The model adopts a double-layer non-uniform discrete grid by dividing the study area into a number of sub-areas. The first layer of coarse grid is used to determine the topology of the pipeline network (Rodrigues et al., 2016; Zhang et al., 2017a), and the second layer of fine grid is used to determine the detailed route between the two selected nodes. As shown in Fig. 3(a), the first layer grid is created based on the location of the well sites, and each node has eight connection directions. This eight-node-type mode is employed to establish the corresponding constraints and objective function for each sub-areas and grid nodes. Supposing fluid flow from node (xi ,yj ) to another adjacent node along its direction k, the terminal node can be defined as (xi ,yj, k). rk is the reverse direction of k, so (xi ,yj, rk) stands for the coordinate of the starting node which connects node (xi ,yj ) at the direction rk. Obviously, (xi ,yj, k) and (xi ,yj, rk) are the same node. As shown in Fig. 3(b), in the second layer grid, the detailed routes of the pipeline can be determined more finely on the three-dimensional terrain by using the ant colony optimization. Based on the accuracy of the terrain data obtained from Google Earth, we set the grid size to 50 m. It should be noted that the connection

routes determined by ant colony optimization are calculated in advance and then used as parameters for the optimization model.

3.2.

Ant colony optimization for optimal route

The ant colony optimization has been widely applied to various engineering optimization problems, particularly in path planning. Inspired by the foraging behavior of ants in finding the shortest path from their nest to food sources, this algorithm was first proposed by M. Dorigo (Dorigo, 1992). During the foraging process, ant colonies randomly choose their path, but leave a kind of chemical compositions called pheromone on their trails when they seek for food. The more pheromones are left on the path, the more possibly other ants will choose this path. Consequently, pheromone concentration on such a path will accumulate faster and attract more ants to follow the same route. Finally, the optimal path is selected through exchanges of information between ant individuals and mutual cooperation (Mokhtari and Rekioua, 2018). In this paper, we use the classical ant colony optimization to search the detailed routes between each connected node in three-dimensional terrain. The objective function is minimizing the distance in three-dimensional terrain. The calculation process of ACO is listed in Table 1. For more details on the spatial path optimization of ant colony algorithms, it is recommended to review the literature (Baeza et al., 2017; Jiang et al., 2015).

326

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

Fig. 3 – Two-layer non-uniform discrete grid. built from node (i, j) to the adjacent node at its direction k with diameter d, BPi,j,k,d =1 and 0 otherwise. It should be noted that the distance disi,j,k between a node (i, j) and the adjacent node is calculated previously by the ant colony optimization.

Table 1 – Calculation process of ACO method. Input topographic parameters Put the number of anm artificial ants at the start location (xstart,ystart) For in =1,2. . .,inm For an = 1,2. . .,anm Calculate the state transition probability in the surrounding 8 directions and use the roulette algorithm to select the next move point until the ant individual reaches the end point or no neighbor nodes are available for selection. End all ants complete the path search Update the pheromone matrix based on ant-cycle model. End Output the optimal route for gathering pipeline. L=



2

2

(xi − xi ) + (yi − yi ) + (zi − zi )

FP =

cpd disi,j,k BPi,j,k,d

(2)

i∈I j∈J k∈K d∈D

The expressions of the construction cost of CPF can be calculated by the following equation. FC =



ccBCi,j

(3)

i∈I j∈J

2

i∈I

anm is the artificial ants number of each generation; inm is the irritation times of ACO, I is the set of points in optimal path.

3.3.



BSi,j,v is the binary variable of pressure devices. If a new supercharging device of type v is built in node (i, j), BSi,j,v = 1 and 0 otherwise. Eq. (4) means that the construction cost of pressure devices should be added to the total cost. FS =

MILP model for optimal pipeline network

 

csv BSi,j,v

(4)

i ∈ IW j ∈ JW v ∈ V

According to the problem statement and parameter preparation in the previous section, a MILP model is proposed to address the optimal design of gathering pipeline networks considering the hydraulic characteristics. All the parameters are denoted with lower-case symbols, and all the decision variables are denoted with upper-case symbols.

3.3.2.

Node and pipeline constraints

Each node cannot belong to well site node (bwi,j ), ordinary node (BNi,j ) or the CPF node (BCi,j ) simultaneously. This constraint is expressed as: bwi,j + BNi,j + BCi,j = 1 i ∈ I, j ∈ J

3.3.1.

(5)

Objective function

The objective function is to work out the construction scheme with the lowest total cost under each given constraint (Eq. 1). The first term refers to the construction cost of pipelines, while the second represents the construction cost of CPF, and the last term is the construction cost of supercharging devices (compressors for gas fields or pumps for oil fields). min F = FP + FC + FS

(1)

The construction cost of pipelines depended on diameter and length should be added to the total cost. BPi,j,k,d is the binary variable of pipeline construction. If a new pipeline is

For a gathering pipeline network, there is generally only one CPF, and once the location is determined, it does not change over time (Zhou et al., 2019).



BCi,j = 1

(6)

i∈I j∈J

BPi,j,k,d represents the construction of pipelines from the node (i,j) to other directions, and BPi,j,rk,d represents the construction of pipelines from other directions to the node (i,j). We define that the direction of pipeline construction is consistent with the direction of fluid flow. There are only one route and

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

one flow direction between two nodes, either from the node (i,j) to other directions or from other directions to the node (i,j). In addition, each node can have up to one route connected with the next node, in other words, each node has at most one flow direction to the next node. This constraint is expressed as:



BPi,j,k,d +

d∈D



BPi,j,rk,d ≤ 1 i ∈ I, j ∈ J, k, rk ∈ K

(7)

d∈D

However, CPF is the terminal of fluid, so it must have no further connection.



BPi,j,k,d ≤ 1 − BCi,j i ∈ I, j ∈ J

(8)

k∈K d∈D

The well site node is the starting point for fluid and must have a further connection through a pipeline with a specific diameter in a certain direction, i.e.,



BPi,j,k,d ≥ bwi,j i ∈ IW , j ∈ JW

327

nodes (QPi,j,k ) and receiving flowrate (QRi,j ). This specific logical relationship is expressed as follows: qwi,j +

 rk ∈ K

QPi,j,rk =



QPi,j,k + QRi,j i ∈ I, j ∈ J

(12)

k∈K

As for the receiving flowrate, all the products are delivered from well sites to the CPF for processing, so on the one hand if the node is a CPF, the receiving flowrate should equal the total producing rate of the well sites, on the other hand if the node is not a CPF, the receiving flowrate must equal to 0. QRi,j = BCi,j qtotal i ∈ I, j ∈ J

(13)

For oil pipelines, the flowrate in pipeline should be kept in the efficient range by selecting appropriate pipeline diameter to meet the constraint of economic flowrate according to the Chinese code GB50350-2015 (China, 2015) and Zhang (Zhang et al., 2017a). QPi,j,k ≥ qpmind BPi,j,k,d i ∈ I, j ∈ J, k ∈ K, d ∈ D

(14)

(9) QPi,j,k ≤ qpmaxd + qtotal(1 − BPi,j,k,d ) i ∈ I, j ∈ J, k ∈ K, d ∈ D

k∈K d∈D

(15)

3.3.3.

Connection structure constraints

If the connection structure of gathering pipeline network is CDPN, which means well sites are connected with the nearest well site through branch pipelines, no special constraints are required. However, for IDPN, the well site node (bwi,j = 1) can only be taken as a starting point to deliver fluid to other nodes, which means that  no fluid flows from other nodes to the well site node ( BPi,j,rk,d = 0). rk ∈ K d ∈ D



BPi,j,rk,d ≤ 1 − bwi,j i ∈ IW , j ∈ JW

(10)

rk ∈ K d ∈ D

3.3.4.

Pipeline flow constraints

3.3.5.

As shown in Fig. 4(a), we define four flow parameters for each node. qwi,j is output of each well site. qwi,j = 0 if node (i, j) is not a well site. QPi,j,rk stands for inflow rate from adjacent node along direction rk to the node (i, j) via pipeline. QPi,j,k stands for outflow rate from node (i, j) to the adjacent node along direction k via pipeline. QRi,j represents receiving flowrate of node (i, j).If node (i, j) is not a CPF then QRi,j = 0. All production are transported from well sites to CPF for processing, therefore, if node (i, j) is a CPF, then QRi,j means the flowrate of oil or gas consumed by the CPF, which equals the total production of all well sites (qtotal =



qwi,j ).

i ∈ IW j ∈ JW

If fluid can flow from the node (i, j) to the adjacent node at its direction k, there must be a pipeline between those two nodes and the flowrate should be less than the total production of all well sites.



QPi,j,k ≤ qtotal

The relationships of (14) and (15) are activated only when the binary variable BPi,j,k,d = 1. When BPi,j,k,d = 0, Eqs. (14) and (15) are not binding as their right hand sides respectively become 0 and (qpmax + qtotal). For gas pipelines, because of the compressibility of gas, the economic flowrate is not suitable, and it is necessary to select the pipe diameter by hydraulic calculation. In addition, the gas flowrate control is often used to prevent erosion due to the gas velocity being too fast (Botros et al., 2018; Parsi et al., 2014), so the pipeline flowrate should be less than the corresponding upper limit, which can be expressed as constraint (15).

Equipment constraints

If the maximum wellhead back pressure of a well site is too large, the throttle valve at the wellhead can be used for throttling without the need for additional throttling equipment. However, if the maximum wellhead back pressure is too small to transport the fluid, a supercharger device is required at the well site. In addition, only a certain kind of device can be installed in well site nodes.



BSi,j,v ≤ bwi,j i ∈ I, j ∈ J

(16)

v∈V

There are different kinds of pressure equipment, and in virtue of their capacity limits, the corresponding equipment should be determined based on the real flowrate.



QPi,j,k ≥ qsminv BSi,j,v i ∈ IW , j ∈ JW , v ∈ V

(17)

k∈K

BPi,j,k,d i ∈ I, j ∈ J, k ∈ K

(11)

d∈D

The flowrate at each node should meet the conservation of flowrate. The sum of production flowrate (qwi,j ) and inflow from other nodes (QPi,j,rk ) equals to the sum of outflow to other

 k∈K

QPi,j,k ≤ qsmaxv + qtotal(1 − BSi,j,v ) i ∈ IW , j ∈ JW , v ∈ V (18)

328

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

Fig. 4 – Flow and pressure parameters for each node. (a) Flow parameters (b) Pressure parameters. The previous relationships of (17) and (18) are activated only when the binary variable BSi,j,v = 1. When BSi,j,v = 0, Eqs. (17) and (18) are not binding as their right hand sides respectively become 0 and (qpmax + qtotal).

3.3.6.

Hydraulic constraints

The hydraulic pressure drop equation for the pipeline from node (i,j) to the adjacent node at its direction k can be simplified as follows (Bagajewicz and Valtinson, 2014): Pm1 − Pm1 − fd disi,j,k QPi,j,k m2 = 0 i,j,k i,j,rk

(19)

i ∈ I, j ∈ J, k, rk ∈ K, d ∈ D

Generally, the oil gathering pipeline is operated within hydraulic smooth area (Liao et al., 2018); in this flow state, 0.25 m1 = 1, m2 = 1.75, fd = 0.0246 d4.75 . As for gas, the flow is usually in the region of quadratic resistance law (Wang et al., 2018b), ∗T m1 = 2, m2 = 2, fd = Z 2 5 . As can be seen from constraint (19), C0 D

the pressure drop equation is a nonlinear equation that would add the complexity of the optimal model. Therefore, the piecewise linear approximation is used to eliminate the nonlinear term QP. By using piecewise approximation, the constraint (19) can be transformed into (20) and (21).





Pm1 − Pm1 − fd disi,j,k QPAi,j,k ≥ ppmaxm1 BPi,j,k,d − 1 i,j i ,j i, i ∈ I, j, j ∈ J, k ∈ K, d ∈ D



− Pm1 − fd disi,j,k QPAi,j,k ≤ ppmaxm1 1 − BPi,j,k,d Pm1 i,j i ,j i, i ∈ I, j, j ∈ J, k ∈ K, d ∈ D

(20)

 (21)

It is important to note that pressure terms in constraint (20) and (21) are handled in a special way. It should be noted that ppmax is a suitably large number that will make the constraints non-binding when BPi,j,k,d = 0. As for oil, m1 equals 1 so the pressure term is linear and does not need to be linearized; as for gas, m1 equals 2 so the pressure term is nonlinear. However, in order to reduce the problem size (number of variables and number of constraints), the pressure terms in flow-pressuredrop relationship are considered as variables themselves, so that there is no need to linearize each pressure term. That is, Pij m1 are the variables (not Pij ) and the constraints are all linear with respect to Pij m1 . Obviously, once the solution of Pij m1 is obtained, the specific value of the pressure can be found by the square root operation. In order to facilitate understanding of the proposed model, we still use Pij m1 to represent the non-

Fig. 5 – Schematic diagram for piecewise linear approximation. linear pressure term, but the whole of Pij m1 is a variable (not Pij ). As shown in Fig. 5, the piecewise linear function is made of several break points (A∼F) and straight-line segments. Constraints (22)–(26) represent the process of obtaining the approximation QPAi,j,k of nonlinear items QPm2 . The QWi,j,k,a are new variables which can be interpreted as ‘weights’ to be attached to the break points. (22)–(24) determine the value of QPi,j,k and QPAi,j,k , while (25) and (26) guarantee that corresponding values of QPi,j,k and QPAi,j,k lie on same straight line segment. BQAi,j,k,a is equal to 1 when the flowrate QPi,j,k is in segment a. If there is no pipeline built from node (i, j) to the adjacent node at its direction k, then no segment can be chosen. For more details on piecewise linear approximation, it is recommended to Model Building in Mathematical Programming section 7.3, 9.3 (Williams, 2013). It should be noted that there are only |A|-1 intervals, leaving BQAi,j,k,|A| as a variable that does nothing. QPi,j,k =



QWi,j,k,a qxa i ∈ I, j ∈ J, k ∈ K

(22)

a∈A

QPAi,j,k =



QWi,j,k,a qya i ∈ I, j ∈ J, k ∈ K

(23)

a∈A

 a∈A

QWi,j,k,a = 1 i ∈ I, j ∈ J, k ∈ K

(24)

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

If a node (i, j) is CPF, the pressure of this node should be within the specified range.

QWi,j,k,a ≤ BQAi,j,k,a i ∈ I, j ∈ J, k ∈ K, a = 1 QWi,j,k,a ≤ BQAi,j,k,a−1 + BQAi,j,k,a i ∈ I, j ∈ J, k ∈ K, a = 2...|A| − 1

(25)

QWi,j,k,a ≤ BQAi,j,k,a−1 i ∈ I, j ∈ J, k ∈ K, a = |A|



|A|−1

BQAi,j,k,a =



BPi,j,k,d i ∈ I, j ∈ J, k ∈ K

(26)

d∈D

a=1

Creating a piecewise approximation creates its own optimization problem that aims to generate the most accurate approximation with the fewest number of break points. Obviously, the more break points are made, the more accurate of the approximation is, but the longer it takes to solve the proposed optimization model. In this paper, the range of flowrate is from min (qw(i,j) ) to qtotal, and the number of break points is determined by average absolute error, which is 1% in this paper based on the measurement accuracy of the field instrument. It should be noted that because of the existence of the constraint 26, when a node has no pipe connection, the flowrate of this node is 0, so 0 also needs to be a break point. However, there will be no other value of flowrate between 0 and minimum output of well site min (qw(i,j) ), so the calculation error will not be increased.

m1 Pm1 BCi,j i ∈ I, j ∈ J i,j ≥ pcpfmin

(31)

m1 Pm1 + ppmaxm1 (1 − BCi,j ) i ∈ I, j ∈ J i,j ≤ pcpfmax

(32)

Recall that whole of the term Pij m1 in Eqs. (27)–(32) is the variable in the model, which means that these equations are all linear. Furthermore, ppmax in Eqs. (30) and (32) is the same parameter used in Eqs. (20) and (21).

3.3.8.

Pipeline pressure constraints

In view of the installation of the pressure equipment, the pressure change caused by the pressure equipment needs to be considered. As shown in Fig. 4(b), we define three pressure parameters for each node. pwi,j is maximum permissible value of wellhead back pressure. pwi,j = 0 if node (i, j) is not a well site. psv stands for pressure provided by supercharging device v. Pi,j represents pressure of node (i, j). If node (i, j) is well site node, then the pressure should be less than the sum of maximum wellhead back pressure and the pressure provided by supercharging device. In addition, only a certain kind of devices can be installed in well site nodes. Therefore, it is expressed as (27) for oil flow (m1 = 1) and (28) for gas flow (m1 = 2). m1 Pm1 i,j ≤ pwi,j +



BSi,j,v psv m1 i ∈ IW , j ∈ JW

(27)

v

m1 Pm1 i,j ≤ pwi,j + 2pwi,j

+





BSi,j,v psv

v

BSi,j,v psv m1 i ∈ IW , j ∈ JW

(28)

Obstacles constraints

During the construction of the gathering pipeline network, some places where the terrain is too undulating or where other facilities already exist are regarded as obstacles (obi ,j = 1). The CPF and pipelines cannot be constructed in the area with obstacles. Therefore, the corresponding binary variables need to take the corresponding values. BCi ,j = 0 i ∈ I, j ∈ J



BPi ,j ,k,d = 0 i ∈ I, j ∈ J

(33) (34)

k∈K d∈D

3.3.9. 3.3.7.

329

Overall construction constraints

If there are already CPF, pipelines and devices in a gathering pipeline network, the corresponding binary variables need to take the corresponding values. BCi ,j = 1 i ∈ I, j ∈ J

(35)

BPi ,j ,k ,d = 1 i ∈ I, j ∈ J, k ∈ K, d ∈ D

(36)

BSi ,j ,v = 1 i ∈ I, j ∈ J, v ∈ V

(37)

3.4.

MILP model solution

The proposed model contains many constraints, including node and pipeline, pipeline flow, equipment, pipeline pressure, obstacles and overall construction. Combined with specific actual data and piecewise linearization method, the model can be transformed into a MILP mathematical model of which all the constraints and the objective function are linear. Therefore, the model can be solved by CPLEX, a commercial MILP model solver, which is based on the branch and bound algorithm (Yokoyama et al., 2015). Hence, the optimal construction scheme of gathering pipeline network can be obtained.

v

If a node (i, j) is connected to other nodes with a pipeline or a node (i, j) is CPF, then the pressure of this node (i, j) should be limited to the safe range, neither exceeding the maximum operating pressure nor lower than the minimum operating pressure. On the contrary, if a node (i, j) has no pipeline or is not CPF, then the pressure of this node (i, j) should be 0. m1 Pm1 ( i,j ≥ ppmin



BPi,j,k,d + BCi,j ) i ∈ I, j ∈ J

(29)

k∈K d∈D m1 ( Pm1 i,j ≤ ppmax

 k∈K d∈D

BPi,j,k,d + BCi,j ) i ∈ I, j ∈ J

(30)

4.

Case studies

In this section, the performance of the proposed optimization methodology is illustrated by its application to two gas fields in China and a virtual oil field. The first two examples deal with the optimization design for gas fields, and the variations of gathering pipeline network when the production rates decline are analyzed. Next, the MILP model is applied to an oil field to show the advantages of the proposed method for providing an optimal solution. The mathematical models were implemented in GAMS 24.1.3 as a programming environment, solved by CPLEX 12.5.1 as the MILP solver and executed on a workstation with an eight-core Intel Xeon E5-2609 (1.70 GHz).

330

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

Fig. 6 – Optimal gathering pipeline network of case study 1. (a) IDPN (b) CDPN.

4.1.

Case study 1

The first case aims to test the accuracy of hydraulic calculation results and the effectiveness of the optimization results, taking a gas field in Linxing area of Shanxi, China as an example. The data are obtained from a natural gas production enterprise that is responsible for the operating management for this gas field. There are five cluster well sites as shown in Fig. S1, including 25 different types of production wells, such as horizontal wells, directional wells and straight wells. The production of each well site varies from 3 × 104 m3 /d to 26 × 104 m3 /d and the specific production information of those five well sites is shown in Table S1. Only one CPF in this area and the construction cost of CPF is 20 × 106 CNY. Moreover, the required pressure into the CPF is 5.4 MPa, the maximum wellhead back pressure is 7.4 MPa and the maximum pipeline network operating pressure is 9 MPa. In addition, the fluid velocity of gas pipelines is controlled to be 4 m/s. There are six kinds of pipes with different diameters to choose from, and the specific information is shown in Table S2. This case for IDPN consists 34,490 constraints, 31,684 variables (13,235 binary variables of reduced MIP), and we obtained the solution in 3027 s with the gap of 5%. As for CDPN, it involves 34,485 constraints, 31,684 variables (13,769 binary variables of reduced MIP) and the optimal was reached in 772 s with the gap of 5%. The piecewise linear approximation was employed and the number of break points is 13 for both IDPN and CDPN. The optimized gathering pipeline network connection and the construction location of CPF are shown in Fig. 6 for IDPN and CDPN, respectively. The length of each pipeline is shown in Table S3. Compared with the actual pipeline construction cost of 6.28 × 106 CNY, the optimal pipeline construction cost can be effectively reduced by 38.4% for IDPN and 44.4% for CDPN. The total construction cost is shown in Fig.7. Compared with the total actual construction cost of 26.28 × 106 CNY, the optimal total construction cost can be effectively reduced by 9.2% for IDPN and 10.6% for CDPN. Due to the pressure drop caused by friction, the hydraulic check for the optimization pipeline network should be necessary and was implemented by PIPESIM, which is a well-known commercial multiphase flow simulation software. The Modelling interface of PIPESIM for IDPN and CDPN are shown in

Fig.S2. The results are shown in Fig. 8. The average error of hydraulic calculation for IDPN is 1.67%, and the average error for CDPN is 1.07%. Therefore, the calculation results of the model can be considered accurate. Furthermore, it can be concluded that the pressure of all the nodes calculated by the proposed model can meet the pressure limit, being greater than the inbound pressure and less than the maximum operating pressure of the pipeline. Moreover, the results of gas velocity in pipeline, as shown in Fig.S3, are close and smaller than the controlled velocity of 4 m/s. In summary, the model is proved to be applied to the optimization of gathering pipeline network with different structures under the premise of meeting hydraulic requirements.

4.2.

Case study 2

The second case presents another gas field in Linxing area of Shanxi, China, and investigates the variations of gathering pipeline network when the bottom-hole pressure and production rate decline. The data come from a research institution responsible for the planning and design of this gas field. In the early stage of gas field development, one CPF and two well sites (LX8 and LX9) have already been built. Now, nine well sites including 32 different types of production wells have been added due to the expansion of production. The specific production information of those well sites is shown in Table S4. Moreover, the required pressure into the CPF is 1.0 MPa, the maximum wellhead back pressure is 3.5 MPa, the maximum and minimum pipeline network operating pressure is 4.0 MPa and 1 MPa, respectively. In addition, the fluid velocity of gas pipelines is controlled to be 8 m/s. There are three kinds of pipes with different diameters to choose from, and the specific information is shown in Table S5. This case for IDPN consists 34,611 constraints, 37,443 variables (12,531 binary variables of reduced MIP), and we obtained the solution in 23,588 s with the gap of 7%. As for CDPN, it involves 34,599 constraints, 37,443 variables (13,656 binary variables of reduced MIP) and the optimal was reached in 28,794 s with the gap of 7%. The piecewise linear approximation was employed and the number of break points is 13 for both IDPN and CDPN. The new construction scheme are shown in Fig. 9 for IDPN and CDPN, respectively. The length of each pipeline is shown

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

331

Fig. 7 – Comparison of construction costs among onsite, IDPN and CDPN.

Fig. 8 – Well site pressure for case study 1. (a) IDPN (b) CDPN.

in Table S6. Compared with the actual pipeline construction cost of 20.76 × 106 CNY, the optimal pipeline construction cost can be effectively reduced by 12.8% for IDPN and 15.6% for CDPN.

In the later stages of gas field development, the amount of production and wellhead back pressure will decay. In this case, the gas may not have sufficient pressure to flow from wellhead to CPF, so it is necessary to change the hydraulic characteristics of the pipeline network by installing pressurized equipment.

332

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

Fig. 9 – Optimal gathering pipeline network of case study 2. (a) IDPN (b) CDPN.

Fig. 10 – Optimal gathering pipeline network of case study 3. (a) IDPN (b) CDPN.

Based on the above optimal pipeline network, we tested the performance of the model under this scenario. It is assumed that the production rate and pressure of well sites will change with time, as shown in Table S7. The corresponding class and cost of equipment are shown in Table S8. This case for IDPN

consists 34,830 constraints, 38,020 variables (8306 binary variables of reduced MIP), and we obtained the solution in 0.69 s with the gap of 0%. As for CDPN, it involves 34,819 constraints, 38,020 variables (8603 binary variables of reduced MIP) and the optimal was reached in 0.78 s with the gap of 0%.

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

333

Fig. 11 – Hydraulic conditions of oil pipelines for case study 3. (a) IDPN (b) CDPN. The implementation results, as shown in Table S9, indicate that the wellhead back pressure in some well sites is not enough and pressurized equipment is needed. More specifically, the pressurized equipment of class V1 is needed in well site LX59 and LX31 for IDPN, the cost is 3.57 × 106 CNY. As for CDPN, the pressurized equipment of class V1, V1, V1 and V2 are needed in well site LX19, LX59, LX58 and LX31, respectively. The cost is 8.15 × 106 CNY. According to the results, it is concluded that the proposed model can identify the boost point and select the appropriate class of pressurized equipment.

4.3.

Case study 3

Case study 3 shows the application of the proposed model in a virtual oil field. There are 10 cluster well sites and the production of each well site varies from 350 m3 /d to 460 m3 /d as shown in Table S10. Only one CPF in this area and the construction cost of CPF is 20 × 106 CNY. Moreover, the required pressure into the CPF is 1.0 MPa, the maximum wellhead back pressure is 5.0 MPa and the maximum pipeline network operating pressure is 7 MPa. In addition, the fluid velocity of oil pipelines is controlled to be 0.8 ∼2.0 m/s according to Chinese design code GB50350-2015. There are three kinds of pipes with different diameters to choose from, and the specific information is shown in Table S11.

This case for IDPN consists 12,990 constraints, 12,723 variables (4092 binary variables of reduced MIP), and we obtained the solution in 40,959 s with the gap of 7%. As for CDPN, it involves 13,000 constraints, 12,723 variables (3348 binary variables of reduced MIP) and the optimal was reached in 45,694 s with the gap of 7%. The piecewise linear approximation was employed and the number of break points is 10 for both IDPN and CDPN. The results are shown in Fig. 10 for IDPN and CDPN, respectively. The total construction cost of IDPN is 76.18 × 106 CNY, and CDPN is 73.91 × 106 CNY. We used the same method as in Case study 1 to verify the hydraulic conditions of the pipeline network obtained by the proposed model. The Modelling interface of PIPESIM for IDPN and CDPN are shown in Fig. S4, and the results for IDPN and CDPN are shown in Fig. 11. Compared with PIPESIM, the average error of hydraulic calculation for IDPN is 8.47%, and the average error for CDPN is 6.87%. The requirements of economic flowrate are also met, as shown in Fig. S5.

5.

Conclusion and future work

This paper focuses on two common connection structures of gathering pipeline networks, and a detailed procedure for optimization of gathering pipeline networks is presented. Considering obstacles, three-dimensional topography,

334

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

pipeline network structure, pipeline diameter, wellhead back pressure and pressure equipment, an integrated MILP model with the minimum total cost as the objective function is established. The ant colony optimization is used to determine the values of the parameters in the model that describe the detailed route of pipeline in three-dimensional terrain. More important, hydraulic characteristics are considered and linearized by a piecewise method to eliminate the nonlinear term of hydraulic constraints. The optimal solution is obtained integrally by CPLEX solver, including the optimal connection topology, the location of the central processing facility and pressure facilities, the diameter and route of each pipeline, pipeline flowrate and node pressure. The feasibility and practicality of proposed MILP model are confirmed with two real-world gas fields and a virtual oil field. The first two examples deal with the optimization design for gas fields, and the variations of gathering pipeline network when the production rates decline are analyzed. Example 3 shows the application of the proposed model in oil field optimization design. Compared with on-site data, the results illustrate that the model lowers the total construction cost, therefore, it can be used to provide theoretical and technical support for gathering pipeline networks optimization. There are other elements of more realistic oil-gas fields that are not catered for in this paper, such as gas injection inline and at the wells, large elevation changes, multiphase flow, the relationship between flowing bottom-hole pressure and production rate in each well, constraints on minimum flowing bottom-hole pressure and on erosional velocity. Moreover, in actual production, because of the rolling development method, it takes a long time to establish the whole gathering pipeline network. The consideration of the time factor in planning may find a more economical result in some cases. Therefore, further efforts will try to consider the above issues for gathering pipeline networks.

Acknowledgments This work was supported by the National Science and Technology Major Project of China (2016ZX05066005-001 and 2016ZX05028004-001), National Natural Science Foundation of China (51874323, 51534007), National Key Research and Development Plan of China (2016YFC0303704), and Science Foundation of China University of Petroleum - Beijing (C201602), all of which are gratefully acknowledged. Thanks to Shangfei Song and Huirong Huang for English language editing.

Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/ j.cherd.2019.08.013.

References Abeysekera, M., Wu, J., Jenkins, N., Rees, M., 2016. Steady state analysis of gas networks with distributed injection of alternative gas. Appl. Energy 164, 991–1002. Al Dossary, M.A., Nasrabadi, H., 2016. Well placement optimization using imperialist competitive algorithm. J. Pet. Sci. Eng. 147, 237–248. Alves, F.D.S., Miranda de Souza, J.N., Hemerly Costa, A.L., 2016. Multi-objective design optimization of natural gas

transmission networks. Comput. Chem. Eng. 93, 212–220. Baeza, D., Ihle, C.F., Ortiz, J.M., 2017. A comparison between ACO and Dijkstra algorithms for optimal ore concentrate pipeline routing. J. Clean. Prod. 144, 149–160. Bagajewicz, M., Valtinson, G., 2014. Computation of natural gas pipeline hydraulics. Ind. Eng. Chem. Res. 53, 10707–10720. Baioco, J.S., Albrecht, C.H., de Lima, B., Jacob, B.P., Rocha, D.M., Asme, 2017. Multi-objective optimization of subsea pipeline routes in shallow waters. Amer Soc Mechanical Engineers, New York. Botros, K.K., Jensen, L., Foo, S., 2018. Determination of erosion-based maximum velocity limits in natural gas facilities. J. Nat. Gas Sci. Eng. 55, 395–405. Chen, J., Yu, J., Ai, B., Song, M., Hou, W., 2018. Determinants of global natural gas consumption and import–export flows. Energy Econ., http://dx.doi.org/10.1016/j.eneco.2018.06.025 (in press). Mo.Ha.U.-R.Dot.Ps.Ro. China, GB 50350-2015 2015. Code for design of oil-gas gathering and transportation systems of oilfield. China Planning Press, Beijing. Crow, D.J.G., Giarola, S., Hawkes, A.D., 2018. A dynamic model of global natural gas supply. Appl. Energy 218, 452–469. de Lucena, R.R., Baioco, J.S., Leite Pires de Lima, B.S., Albrecht, C.H., Jacob, B.P., 2014. Optimal design of submarine pipeline routes by genetic algorithm with different constraint handling techniques. Adv. Eng. Softw. 76, 110–124. de Padua Agripa Sales, L., Pitombeira-Neto, A.R., de Athayde Prata, B., 2018. A genetic algorithm integrated with Monte Carlo simulation for the field layout design problem. Oil Gas Sci. Technol. 73, 16. Ding, S., 2018. A novel self-adapting intelligent grey model for forecasting China’s natural-gas demand. Energy 162, 393–407. Diogo, A.F., Barros, L.T., Santos, J., Temido, J.S., 2018. An effective and comprehensive model for optimal rehabilitation of separate sanitary sewer systems. Sci. Total Environ. 612, 1042–1057. Dorigo, M.J.P.T., 1992. Politecnico di Milano. In: Optimization, Learning and Natural Algorithms. Drouven, M.G., Grossmann, I.E., 2016. Multi-period planning, design, and strategic models for long-term, quality-sensitive shale gas development. AIChE J. 62, 2296–2323. Drouven, M.G., Grossmann, I.E., 2017. Mixed-integer programming models for line pressure optimization in shale gas gathering systems. J. Pet. Sci. Eng. 157, 1021–1032. Fister, I., Perc, M., Kamal, S.M., Fister, I., 2015. A review of chaos-based firefly algorithms: perspectives and research challenges. Appl. Math. Comput. 252, 155–165. Gupta, V., Grossmann, I.E., 2012. An efficient multiperiod MINLP model for optimal planning of offshore oil and gas field infrastructure. Ind. Eng. Chem. Res. 51, 6823–6840. He, G.X., Tang, D.D., Yin, B.B., Sun, L.Y., Ding, D.Q., Liang, Y.T., Liao, K.X., 2018. Comparison and analysis of drainage measures for draining accumulated water condensed from wet CBM and transported in surface gathering pipeline network. J. Nat. Gas Sci. Eng. 56, 281–298. Hong, C., Estefen, S.F., Wang, Y., Lourenco, M.I., 2018. An integrated optimization model for the layout design of a subsea production system. Appl. Ocean Res. 77, 1–13. Jiang, W.-Y., Lin, Y., Chen, M., Yu, Y.-Y., 2015. A co-evolutionary improved multi-ant colony optimization for ship multiple and branch pipe route design. Ocean. Eng. 102, 63–70. Johnson, N., Ogden, J., 2012. A spatially-explicit optimization model for long-term hydrogen pipeline planning. Int. J. Hydrogen Energy 37, 5421–5433. Jonsson, M.T., Magnusdottir, L., 2017. Minimizing visual effects and optimizing routes and locations for geothermal steam gathering system. Kan, S.Y., Chen, B., Wu, X.F., Chen, Z.M., Chen, G.Q., 2019. Natural gas overview for world economy: from primary supply to final demand via global supply chains. Energy Policy 124, 215–225. Kim, C., Kim, K., Kim, J., Ahmed, U., Han, C., 2018. Practical deployment of pipelines for the CCS network in critical

Chemical Engineering Research and Design 1 5 2 ( 2 0 1 9 ) 320–335

conditions using MINLP modelling and optimization: a case study of South Korea. Int. J. Greenhouse Gas Control. 73, 79–94. Kostowski, W.J., Skorek, J., 2012. Real gas flow simulation in damaged distribution pipelines. Energy 45, 481–488. Liao, Q., Liang, Y., Xu, N., Zhang, H., Wang, J., Zhou, X., 2018. An MILP approach for detailed scheduling of multi-product pipeline in pressure control mode. Chem. Eng. Res. Des. 136, 620–637. Marcoulaki, E.C., Papazoglou, I.A., Pixopoulou, N., 2012. Integrated framework for the design of pipeline systems using stochastic optimisation and GIS tools. Chem. Eng. Res. Des. 90, 2209–2222. Marcoulaki, E.C., Tsoutsias, A.V., Papazoglou, I.A., 2014. Design of optimal pipeline systems using internal corrosion models and GIS tools. Ind. Eng. Chem. Res. 53, 11755–11765. Mezura-Montes, E., Coello Coello, C.A., 2011. Constraint-handling in nature-inspired numerical optimization: past, present and future. Swarm Evol. Comput. 1, 173–194. Mokhtari, Y., Rekioua, D., 2018. High performance of maximum power point tracking using ant colony algorithm in wind turbine. Renew. Energy 126, 1055–1063. Morvaj, B., Evins, R., Carmeliet, J., 2016. Optimising urban energy systems: simultaneous system sizing, operation and district heating network layout. Energy 116, 619–636. Pambour, K.A., Bolado-Lavin, R., Dijkema, G.P.J., 2016. An integrated transient model for simulating the operation of natural gas transport systems. J. Nat. Gas Sci. Eng. 28, 672–690. Parsi, M., Najmi, K., Najafifard, F., Hassani, S., McLaury, B.S., Shirazi, S.A., 2014. A comprehensive review of solid particle erosion modeling for oil and gas wells and pipelines applications. J. Nat. Gas Sci. Eng. 21, 850–873. Rios-Mercado, R.Z., Borraz-Sanchez, C., 2015. Optimization problems in natural gas transportation systems: a state-of-the-art review. Appl. Energy 147, 536–555. Rodrigues, H.W.L., Prata, B.A., Bonates, T.O., 2016. Integrated optimization model for location and sizing of offshore platforms and location of oil wells. J. Pet. Sci. Eng. 145, 734–741. Rosa, V.R., Camponogara, E., Ferreira, V.J.M., 2018. Design optimization of oilfield subsea infrastructures with manifold placement and pipeline layout. Comput. Chem. Eng. 108, 163–178. Tian, Q., Zhao, D., Li, Z., Zhu, Q., 2017. Robust and stepwise optimization design for CO2 pipeline transportation. Int. J. Greenhouse Gas Control. 58, 10–18. ˘ Üster, H., Dilaveroglu, S¸., 2014. Optimization for design and operation of natural gas transmission networks. Appl. Energy 133, 56–69. Wang, B., Yuan, M., Zhang, H., Zhao, W., Liang, Y., 2018a. An MILP model for optimal design of multi-period natural gas transmission network. Chem. Eng. Res. Des. 129, 122–131. Wang, B., Yuan, M., Zhang, H., Zhao, W., Liang, Y., 2018b. An MILP model for optimal design of multi-period natural gas transmission network. Chem. Eng. Res. Des. 129, 122–131.

335

Wang, Y.Y., Duan, M.L., Feng, J.K., Mao, D.F., Xu, M.H., Estefen, S.F., 2014. Modeling for the optimization of layout scenarios of cluster manifolds with pipeline end manifolds. Appl. Ocean Res. 46, 94–103. Wang, Y.Y., Duan, M.L., Xu, M.H., Wang, D.G., Feng, W., 2012. A mathematical model for subsea wells partition in the layout of cluster manifolds. Appl. Ocean Res. 36, 26–35. Wei, L., Dong, H., Zhao, J., Zhou, G., 2016. Optimization model establishment and optimization software development of gas field gathering and transmission pipeline network system. J. Intell. Fuzzy Syst. 31, 2375–2382. Williams, H.P., 2013. Model Building in Mathematical Programming. John Wiley & Sons. Xiang, W.L., An, M.Q., Li, Y.Z., He, R.C., Zhang, J.F., 2014. An improved global-best harmony search algorithm for faster optimization. Expert Syst. Appl. 41, 5788–5803. Yang, X.M., Yuan, J.S., Yuan, J.Y., Mao, H.N., 2007. A modified particle swarm optimizer with dynamic adaptation. Appl. Math. Comput. 189, 1205–1213. Yokoyama, R., Shinano, Y., Taniguchi, S., Ohkura, M., Wakui, T., 2015. Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation. Energy Convers. Manage. 92, 92–104. Yu, W., Song, S., Li, Y., Min, Y., Huang, W., Wen, K., Gong, J., 2018a. Gas supply reliability assessment of natural gas transmission pipeline systems. Energy 162, 853–870. Yu, W., Wen, K., Min, Y., He, L., Huang, W., Gong, J., 2018b. A methodology to quantify the gas supply capacity of natural gas transmission pipeline system using reliability theory. Reliab. Eng. Syst. Saf. 175, 128–141. Zhang, H., Liang, Y., Wu, M., Qian, C., Li, K., Yan, Y., 2015. Study on the optimal topological structure of the producing pipeline network system of CBM fields. International Petroleum Technology Conference. Zhang, H., Liang, Y., Ma, J., Qian, C., Yan, X., 2017a. An MILP method for optimal offshore oilfield gathering system. Ocean. Eng. 141, 25–34. Zhang, H., Liang, Y., Zhang, W., Wang, B., Yan, X., Liao, Q., 2017b. A unified MILP model for topological structure of production well gathering pipeline network. J. Pet. Sci. Eng. 152, 284–293. Zhou, J., Gong, J., Li, X., Tong, T., Cheng, M., Zhang, S., 2014. Optimization of coalbed methane gathering system in China. Adv. Mech. Eng. 6, 147381. Zhou, J., Liang, G., Deng, T., 2018. Optimal design of star-tree oil-gas pipeline network in discrete space. J. Pipeline Syst. Eng. Pract., 9. Zhou, J., Liang, G., Deng, T., Gong, J., 2017. Route optimization of pipeline in gas-liquid two-phase flow based on genetic algorithm., pp. 1–9. Zhou, J., Peng, J., Liang, G., Deng, T., 2019. Layout optimization of tree-tree gas pipeline network. J. Pet. Sci. Eng. 173, 666–680.