Global Optimization of Gas Lifting Operations: A Comparative Study of

6098

Ind. Eng. Chem. Res. 2009, 48, 6098–6104

Global Optimization of Gas Lifting Operations: A Comparative Study of Piecewise Linear Formulations Ruth Misener, Chrysanthos E. Gounaris, and Christodoulos A. Floudas* Department of Chemical Engineering, Princeton UniVersity, Princeton, New Jersey 08544-5263

Continuous gas lifting is the process of increasing oil well production by injecting compressed natural gas, called “lift gas”, into the production tubing of an oil well [Buitrago et al. Presented at the SPE Gas Technology Symposium, Calgary, Alberta, Canada, 1996]. This paper considers the problem of optimizing the distribution of a limited supply of lift gas to wells in an oil field using piecewise linearization techniques. Four modeling approaches, proposed by Nemhauser and Woolsey [Integer and Combinatorial Optimization; : New York, 1988], Foudas [Nonlinear and Mixed-Integer Optimization: Fundamentals and Applications; : New York, 1995], Sherali [Oper. Res. Lett. 2001, 28, 155.], and Keha et al. [Oper. Res. Lett. 2004, 32, 44.], are presented and the gas lifting problem is solved using each method. Each of the four frameworks is sufficient to solve the problem to global optimality and the method presented by Keha et al. has the best computational performance. The gas lifting problem is used within the context of larger problems such as well scheduling in oil fields [Kosmidis et al. Comput. Chem. Eng. 2005, 29, 1523.], and this proposed work can be used to make one of the key parts of the well scheduling problem more efficient. 1. Introduction and Literature Review 1.1. Gas Lifting Introduction. Continuous gas lifting is the process of injecting compressed natural gas, called “lift gas”, into the production tubing of an oil well.1 The lift gas reduces the weight of the column of fluid, increasing the pressure gradient between the reservoir and well fluid that pushes the fluid to the surface.7 The fluid from the well, which includes oil, gas, and water, is treated to separate water. The oil is sold, and the remaining gas is either sold or recycled for use in another gas lifting operation.8 The relevant optimization problem at hand is gas allocation to a field of oil wells. Assuming that reservoir dynamics are constant over the period of a few days,9 engineers construct gas requirement curves for each well that relate gas injection (qGAS,i) to oil production (qOIL,i).10 The index i ∈{1,..., NWELLS} represents the varying response of each individual well to gas lift. The gas requirement curves account for individual well factors such as reservoir pressure, composition, and depth. The correlation qGAS,i vs qOIL,i can be calculated by interpolating between the entries of a hydraulic table.9 An interesting feature of these gas requirement curves is that they are piecewise linear functions, that is, they are defined by vertices (qGAS,i,k, qOIL,i,k) where i ∈{1,..., NWELLS} and k ∈{1,..., NVERTICES}. The scope of this paper is the solution of the gas lifting problem using four methods for optimizing over piecewise linear functions. The four approaches applied are a classic formulation which is presented by both Nemhauser and Woolsey2 and Floudas3 and three recently introduced methods by Floudas,3 Sherali,4 and Keha et al.5 All four methods can solve the gas lifting problem to global optimality and in this paper we investigate their comparative computational performance. 1.2. Literature Review. Multiple groups have studied well scheduling and planning on oil and gas fields6,9,11-16 Given the price of natural gas and the compression cost of injecting lift gas into an oil well, Saputelli et al.13 identify the well-scheduling subproblem of lift gas allocation as an important supervisory control operation. Because of its immediate relevance to oil * To whom all correspondence should be addressed. E-mail: [email protected]. Tel.: (609) 258-4595. Fax: (609) 258-0211.

production, interest in economically allocating lift gas to oil field wells has been an area of research since 1981 when Kanu et al.17 computed the gas requirement curve (qOIL,i vs qGAS,i) for six wells, calculated the gas injection necessary for optimal profit in each well, and found the optimal gas supply to the test case under unlimited lift gas conditions. Kanu et al.17 also considered the case of limited lift gas for the six-well test case. Later, Buitrago et al.1 considered insufficient lift gas supply to reach economic optimality, effectively equating maximum oil production and maximum profit. This condition of limited lift gas, which was adopted by Mason et al.,10 will also be used in this paper. Past attempts to solve this problem have encountered varying success. Kanu et al.17 optimized a six-well test case using an algorithm that requires continuously differentiable, concave gas requirement curves. The assumption of continuously differentiable gas requirement curves was also made by Nishikiori et al.,18 who locally optimized the gas lift problem using a quasiNewton method. Fang and Lo19 transformed the gas requirement curves into a piecewise-linear representation and optimized gas lift problems for concave curves. Observing that many gas

Figure 1. Illustration of the Buitrago et al. literature test case.1

10.1021/ie8012117 CCC: $40.75  2009 American Chemical Society Published on Web 12/15/2008

Ind. Eng. Chem. Res., Vol. 48, No. 13, 2009

6099

Table 1. Gas Lifting Nomenclature indices data methods

variables

parameters

i k (qGAS,i,k, qOIL,i,k) CM LSM CHM SSM QOIL,TOTAL λ ∈ [0, 1]NWELLNVERTICES λL ∈ [0, 1]NWELLNVERTICES-1 λR ∈[0, 1]NWELLNVERTICES-1 y ∈ {0, 1}NWELLNVERTICES-1 QGAS ∈ RNWELLNVERTICES-1 z ∈ {0, 1}NWELL Ri,k βi,k QGAS,LIMIT QGAS,MIN

requirement curves are not necessarily concave (some wells require a minimum gas supply threshold to produce oil), Buitrago et al.1 used a multistart, stochastic algorithm to distribute lift gas to 6- and 56-well test cases. But, as the results of Mason et al.10 show, the algorithm of Buitragoet al.1 rarely attained global optimality. To improve the results of the Buitrago et al.1 stochastic algorithm, Mason et al.10 introduced a “discrete gradient method” (DGM) that adapts standard descent methods for the case of piecewise linear functions. The DGM has been used to solve a variety of optimization problems.20 Mason et al.10 employ four realizations of DGM to solve the problems posed by Kanu et al.17 and Buitrago et al.1 and their own 58-well test case. Although this method is introduced as a “global optimization algorithm”, the four variants of DGM produce different local solutions. Additionally, the success of all four algorithms depends on the initialization point. Camponogara and Nakashima8,21-23 considered the gas lifting problem. Their first solution technique, outlined in Camponogara and Nakashima21 and Nakashima and Camponogara,23 uses a dynamic programming algorithm. But, as Camponogara and Nakashima explain,22 their dynamic programming method produces solutions that are only “close to optimal”. Their second algorithm, described in Camponogara and Nakashima,8,22 begins with the piecewise linear method introduced by Sherali4 and then constructs lifting inequalities to make the problem run more quickly. Camponogara and Nakashima8,22 claim that the lifting inequalities are effective in attaining global optimality, but these papers neither use standard literature test cases, nor do they report the test cases they constructed. In an insightful generalization of the gas lifting problem, Kosmidis et al.6,9 point out that it is often best to integrate factors like pressure drop across well tubing and merging of oil well lines into oil well scheduling optimization. Kosmidis et al.6,9 account for these factors by adding a choke setting variable that controls well tubing pressure, thereby controlling the oil flow rate. Kosmidis et al.6,9 also add additional constraints to better represent multiphase flow in well tubing. To address this class of well scheduling problems, Kosmidis et al.6 outline a solution procedure for simultaneously optimizing over well, manifold (well stream mixers), and separator (water removal facilities) nodes. Kosmidis et al.6 solve the gas lifting problem described in this report for many pressure drops. The method they employed to solve the piecewise linear problem corresponds to the classic algorithm, presented in section 2.2.1. In this paper, we study four methods to solve the gas lifting problem, suggesting improvements to the piecewise linear algorithm of Kosmidis et al.6,9

well index i ranges from 1 to NWELLS index of vertex points between 1 and NVERTICES vertex pairs determined using a reference table classic method (section 2.2.1) linear segmentation method (section 2.2.2) convex hull method (section 2.2.3) special structure method (section 2.2.4) objective to maximize convex combination weight used in CM and SSM Convex combination weight used in CHM convex combination weight used in CHM binary decision variable for CM, LSM, and CHM continuous variable used in LSM binary variable determines if well is activated intercept of a line segment used as in LSM slope of a line segment used as in LSM constraint on the availability of lift gas constraint on the lower bound of the gas rate

2. Problem Formulation 2.1. Problem Statement. The data sets for this problem appear as coordinates (qGAS,i,k, qOIL,i,k) where i ∈{1,..., NWELLS} denotes the well index and k ∈{1,..., NVERTICES} represents the vertex point. Figure 1 illustrates the data given by the six-well1 literature test case. Notice that Well 6 in Figure 1 does not instantaneously begin producing oil when lift gas is injected. This lagging phenomena is a source of problem nonconvexity that forces us to use mixed-integer linear programming rather than linear programming.19 If we were to assume that pathological wells such as Well 61 do not exist in the field, then we could use the linear programming techniques of Fang and Lo19 to address the gas lifting problem. Our goal is to choose the gas injection level for each well so as to maximize oil production. Although this study assumes limited gas supply, note from Figure 1 the importance of accounting for the diminishing returns of injecting large quantities of lift gas under unlimited supply conditions. Other work in the literature, like the studies by Camponogara and Nakashima8,21,22 and Kosmidis et al.9 assume that there is enough lift gas available but that the cost of compressing the lift gas will linearly shift the objective function by the volume of injected gas. Under this assumption, the objective function is the following: max(QOIL,TOTAL - CCOMPRESSIONQGAS,TOTAL) QGAS

(1)

rather than: max QOIL,TOTAL QGAS

(2)

This difference in the objective function will change the global optimum, but the addition of linear terms will not significantly change the performance of the four mathematical frameworks this study compares. Furthermore, the objective function in eq 2 is the one used in the literature test cases of Buitrago et al.1 and Mason et al.,10 allowing us to confirm the solutions. Table 1 defines the nomenclature used in this paper. 2.2. Problem Formulation via Four Approaches. The four approaches presented in sections 2.2.1-2.2.4 are logically equivalent, but they use different constraint representations with the aim of producing shorter solution times. 2.2.1. Classic Method. The classic method, which serves as a baseline for designing the other three approaches, is described by Nemhauser and Woolsey2 and Floudas.3 Its formulation, depicted in Figure 2, introduces λ ∈ [0, 1]NWELLNVERTICES, a continuous variable that serves as a convex combination weight, and y ∈ {0, 1}NWELLNVERTICES-1, a binary variable that restricts the position of nonzero convex combination weights.

6100


Figure 2. Classic method using binary decision variable y to activate the appropriate segment and convex combination weight λ to find the optimal point along that segment.

Figure 4. Convex hull method using binary decision variable y to activate the appropriate segment and left (λL) and right (λR) convex combination weights to interpolate between segment end points.

used as a continuous variable that ranges along a well’s active linear segment. Rather than directly describing the objective function as a convex combination of the vertices, the objective function of the linear segmentation method is composed of a sum of linear functions: NWELL NVERTICES

max QOIL,TOTAL )

y,QGAS

∑ ∑ i

(Ri,k)(yi,k) + (βi,k)(QGAS,i,k)

k

(10) where Ri,k and βi,k are defined as Ri,k ) qOIL,i,k - βi,kqGAS,i,k

(11)

Figure 3. Linear segmentation method uses binary decision variable y to activate the appropriate segment and parameters R and β to represent the intercept and slope, respectively, of the segments.

qOIL,i,k+1 - qOIL,i,k , qGAS,i,k+1 - qGAS,i,k ∀ i ) 1, ..., NWELL, k ) 1, ..., NVERTICES-1

(12)

The classic method features the objective function as a convex combination of the oil exiting at each vertex point:

The variable QGAS,i,k will be nonzero on at most one line segment, so QGAS,i,k is constrained by binary variable yi,k:

NWELLS NVERTICES

max QOIL,TOTAL ) y,λ

∑ i

∑

(λi,k)(qOIL,i,k)

(3)

k

(qGAS,i,k)(yi,k) e QGAS,i,k e (qGAS,i,k+1)(yi,k), ∀ i ) 1, ..., NWELL, k ) 1, ..., NVERTICES-1(13) NVERTICES-1

where λi,k is constrained by λi,1 e yi,1 ,

βi,k )

∀ i ) 1, ..., NWELL

(4)

∀ i ) 1, ..., NWELL, k ) 2, ..., NVERTICES-1 (5) ∀ i ) 1, ..., NWELL

λg0

yi,k ) 1,

∀ i ) 1, ..., NWELL

(6)

∀ i ) 1, ..., NWELL, k ) 1, ..., NVERTICES-1 (15)

Additionally, the sum of the gas injected into each well is constrained by the field limit: NWELL NVERTICES

(7)

∑ ∑ i

NVERTICES-1

∑

yi,k ) 1,

∀ i ) 1, ..., NWELL ; yi,k ∈ {0, 1} (8)

k)1

Since one of the challenges of gas lifting operations is the limited supply of lift gas, we add an extra constraint that limits the supply of natural gas to an upper bound: NWELL NVERTICES

∑ ∑ i

(λi,k)(qGAS,i,k) e QGAS,LIMIT

(14)

k)1

yi,k ∈ {0, 1},

λi,k e yi,k-1 + yi,k, λi,NVERTICES e yi,NVERTICES-1 ,

∑

(9)

k

2.2.2. Linear Segmentation Method. The linear segmentation method, depicted in Figure 3, is presented by Floudas.3 It uses binary variable y ∈ {0, 1}NWELLNVERTICES-1 to restrict each well to a single linear segment, but now QGAS ∈ RNWELLNVERTICES-1 is

QGAS,i,k e QGAS,LIMIT

(16)

k

2.2.3. Convex Hull Method. The convex hull representation, a modified version of the classic method presented in section 2.2.1, was introduced by Sherali.4 In this formulation, which is illustrated in Figure 4, the left and right convex combination weights, λL ∈ [0, 1]NWELLNVERTICES-1 and λR ∈ [0, 1]NWELLNVERTICES-1, are continuous variables and y ∈ {0, 1}NWELLNVERTICES-1 is a binary variable. The interesting feature of the convex hull method is that, for piecewise-linear functions, the convex hull of the constraint polytope (eqs 19-22) is equal to the linear programming relaxation of the constraint polytope. Keha et al.5 explains that this equality makes the formulation locally ideal. Local ideality is stated as follows:


PLP ) {(λ, y)|(λ, y) satisfies (19)-(22)} ) conv{PLP ∩ {(λ, y)|y ∈ {0, 1}}} ) PC (17)

∑

λi,k ) 1,

NWELL NVERTICES y,λL,λR

∑ ∑ i

(λRi,k)(qOIL,i,k) + (λLi,k)(qOIL,i,k)

k

(18) R L and λi,k can be described by where λi,k

λRi,k + λLi,k ) yi,k,

∀ i ) 1, ..., NWELL, k ) 1, ..., NVERTICES-1 (19)

NVERTICES-1

∑

yi,k ) 1,

∀ i ) 1, ..., NWELL

(20)

k)1

(λR, λL) g 0

(21)

yi,k ∈ {0, 1}

(22)

As in the other methods, the total lift gas injected into the field cannot exceed the field limit introduced by the following constraint: NWELL NVERTICES

∑ ∑ i

(λRi,k)(qGAS,i,k) + (λLi,k)(qGAS,i,k) e QGAS,LIMIT

k

(23) 2.2.4. Special Structure Method. The final representation of the problem, depicted in Figure 5, was created using the method of Keha et al.5 This formulation, like the convex hull representation (section 2.2.3), is derived from the classic method (section 2.2.1). The special structure representation recognizes the fact that any two λi,k values should be consecutive but does this without introducing binary variable yi,k. So, although the only variable used in the special structure representation is the convex combination weight λ ∈ [0, 1]NWELLNVERTICES, the variable λ is constrained using a special structure. The special structure method has a similar objective function to the classic formulation: NWELL NVERTICES

max QOIL,TOTAL ) λ

∑ ∑ i

(λi,k)(qOIL,i,k)

(24)

k

but the constraints become

Figure 5. Special structure method uses only convex combination weights λ as variables and takes advantage of the property that at most two consecutive variables are activated.

∀ i ) 1, ..., NWELL

(25)

k)1

The objective function is similar to the one defined in eq 3, but now, the convex weights are disaggregated, that is, they are envisioned as left and right end points: max QOIL,TOTAL )

6101

NVERTICES

λg0

(26)

λ ∈ SOS2

(27)

where SOS2, or special ordered set 2, stipulates that, ∀ i ) 1,..., NWELL, at most two λi,k can be positive and, if two λi,k are positive, they must be consecutive. Keha et al.24 starts from the constraint set composed of eqs 25-27 and designs a branch-and-cut algorithm to enforce the SOS2 property that was first proposed by Beale and Tomlin.25,26 However, implementation of the branch-and-cut algorithm is unnecessary if an elaborate mixed-integer linear programming (MILP) solver that directly admits SOS2 variable sets (e.g., GAMS/CPLEX27,28) is used. The last constraint in the special structure method is the limited availability of lift gas: NWELL NVERTICES

∑ ∑ i

(λi,k)(qGAS,i,k) e QGAS,LIMIT

(28)

k

3. Computational Results 3.1. Literature Test Cases. The four formulations described in sections 2.2.1-2.2.4 CM classic method (section 2.2.1) LSM linear segmentation method (section 2.2.2) CHM convex hull method (section 2.2.3) SSM special structure method (section 2.2.4) were used to solve each of the following test cases: • 6-well test case with 4600 Mscf/day supply limit from the work of Buitrago et al.1 • 56-well test case with 22 500 Mscf/day supply limit from the work of Buitrago et al.1 • 58-well test case with 50 000 Mscf/day supply limit from the work of Mason et al.10 Table 2 compares the results of solving the four mathematical formulations in GAMS27 using the solver CPLEX28 (version 9.0.2) on a Pentium 4 running Linux. Because each of the four methods is formulated as a mixed-integer linear program (MILP), GAMS/CPLEX certifies global optimality for every test case regardless of initialization point. Table 2 also includes the best local solution attained by the DGM optimization algorithms of Mason et al.10 Note that Mason et al.10 actually reported at least four trial runs for each test case, no one of which is consistently the top performer. The “root node” rows of Table 2 refer to the optimal solution of the continuous relaxation of the problem. The root nodes overestimate the true objective value because the problem has been posed as a maximization of oil production. 3.2. Combination of Cases. Camponogara and Nakashima8,22 report their use of minimal lifting cover inequalities to solve gas lifting problems that “obey the patterns typically encountered in practice.” They approximate their well performance curves with 19 linear segments and then perform computational experiments to compare the Sherali4 formulation (CHM in section 2.2.3) to their lifting technique. The reason that Camponogara and Nakashima8,22 need these lifting inequalities is that, even with a test case containing only 32 wells, their implementation of the convex hull method requires as many as 8050 nodes to attain optimality. Their implementation of a 128 well test case visits as many as 2 × 107 nodes. Camponogara and Nakashima8,22 reported that they needed complex machinery to solve the gas lifting problem using ILOG CPLEX (version 9.0) on a Pentium 4 running Linux (our

6102


Table 2. Comparison of the Four Formulations to the Best “Discrete Gradient Method” Trial Runs Using Literature Test Cases Mason et al.

CM

LSM

CHM

SSM

3665.6 3669.1 46 5 0.01

3665.6 3669.1 12 1 0.00

22632.0 22637.4 163 2 0.02

22632.0 22637.5 104 2 0.01

119592.3 119592.3 165 2 0.04

119592.3 119592.3 105 0 0.00

1

6-well test case with 4600 Mscf/day supply limit objective (B/day) root node (B/day) no. of iterations nodes visited CPU (s)

3665.6

3665.6 3665.6 43 0 0.00

183 0.09

3665.6 3669.1 71 1 0.00

56-well test case with 22500 Mscf/day supply limit1 objective (B/day) root node (B/day) no. of iterations nodes visited CPU (s)

21213.7

22632.0 22637.5 175 3 0.04

10608 709.2

22632.0 22892.4 514 16 0.04

58-well test case with 50000 Mscf/day supply limit10 objective (B/day) root node (B/day) no. of iterations nodes visited CPU (s)

115951.0

119592.3 119592.3 167 0 0.01

1345 11.52

Table 3. 200-Well Case (Complication of a 126-Well Concatenation of Literature Test Cases) obj val (B/day)

root node val (B/day)

no. iter

no. nodes

119592.3 119592.3 516 6 0.02

Table 4. 200-Well Case (Complication of a 126-Well Concatenation of Literature Test Cases) with an Enforced Lower Bound on the Gas Flow Rate, QGAS,MIN(i) ) 100

CPU (s) obj val (B/day)

supply limit: 3000 Mscf/day CM LSM CHM SSM

66188.5 66188.5 66188.5 66188.5

CM LSM CHM SSM

82436.4 82436.4 82436.4 82436.4

CM LSM CHM SSM

92124.1 92124.1 92124.1 92124.1

CM LSM CHM SSM

165851.6 165851.6 165851.6 165851.6

66188.5 66188.5 66188.5 66188.5

731 1236 191 481

0 0 0 0

0.06 0.06 0.06 0.03

1 5 5 2

0.15 0.15 0.10 0.03

0 0 8 0

0.06 0.10 0.11 0.03

0 0 0 0

0.06 0.12 0.06 0.02

17 18 22 13

0.28 0.26 0.15 0.05

supply limit: 7000 Mscf/day 82436.4 82436.4 82436.4 82436.4

708 1368 233 468


707 1429 250 461


698 1480 356 429


186537.3 186537.3 186537.3 186537.3

186543.4 186597.2 186543.4 186543.4

965 1780 450 454

computational results are reported using CPLEX version 9.0.2 on a Pentium 4 running Linux). A new case study consisting of 200 wells verifies that the four piecewise linear representations presented in sections 2.2.1-2.2.4 are sufficient to solve complicated versions of the gas lifting problems. The 200-well case has similarities to a 126-well concatenation of the four literature test cases: the 6- and 56-well cases from Buitrago et al.,1 the 58-well case from Mason et al.10 and the 6-well case from Kanu et al.17 The well performance curves are approximated with as many as 11 linear segments. To further complicate the case, almost half (94) of the 200 wells in the new test case exhibit nonconvexity in the form of the noninstantaneous oil production described in section 2.1. This new test case study is reported in the Supporting Information which is available online.


no. iter

no. nodes

CPU (s)


66188.5 66188.5 66188.5 66188.5

CM LSM CHM SSM

82285.7 82285.7 82285.7 82285.7

66188.5 66188.5 66188.5 66188.5

126 1093 87 79

1 0 0 0

0.21 0.12 0.07 0.03

12 8 15 0

0.36 0.34 0.30 0.03

2 3 4 0

0.16 0.26 0.20 0.02

22 8 16 0

0.19 0.32 0.22 0.02

212 272 206 18

0.69 1.35 0.91 0.08


737 1326 200 102


91905.4 91905.4 91905.4 91905.4

CM LSM CHM SSM

165728.9 165728.9 165728.9 165728.9

CM LSM CHM SSM

186467.7 186467.7 186467.7 186467.7

91905.4 91905.4 91905.4 91905.4

150 1337 154 101


344 1516 307 176


705 2365 871 280

Table 3 analyzes the performance of the four formulations with respect to different lift gas supply limits. All three methods continue to quickly compute the optimal distribution of the lift gas. Notice the close correspondence between root node and objective values in Table 3. For example, the CHM returns 82 436.43 B/day at the root node and 82 436.41 B/day as a final objective value for the case of a 7000 Mscf/day supply limit. These four formulations do not always solve the problem at the root node, but the root nodes are within 1% of the true objective value. The root node values of the CM, CHM, and SSM are identical because the variation between the formulations is their method of constraining convex combination weights. Since constraining the convex combination weights is

Ind. Eng. Chem. Res., Vol. 48, No. 13, 2009 Table 5. 200-Well Case (Complication of a 126-Well Concatenation of Literature Test Cases) with an Enforced Lower Bound on the Gas Flow Rate, QGAS,MIN(i) ) 500 obj val (B/day)


no. iter

no. nodes

CPU (s)


65272.0 65272.0 65272.0 65272.0

CM LSM CHM SSM

79201.9 79201.9 79201.9 79201.9

CM LSM CHM SSM

88117.2 88117.2 88117.2 88117.2

65272.0 65272.0 65272.0 65272.0

147 885 111 97

0 0 0 0

0.04 0.1 0.05 0.01

6 20 17 0

0.17 0.34 0.27 0.01

57 131 217 288

0.59 1.00 1.00 0.81

1837 210 2029 72

4.63 0.98 5.05 0.28

3730 224 836 109

9.12 1.29 2.79 0.33


184 1084 211 117


1325 2467 1036 950


161756.3 161756.3 161756.3 161756.3

CM LSM CHM SSM

182965.8 182965.8 182965.8 182965.8

162034.3 162090.1 162034.3 162034.3

7846 1996 5300 906


9522 2232 2720 772

Table 6. 200-Well Case (Complication of a 126-Well Concatenation of Literature Test Cases) with an Enforced Lower Bound on the Gas Flow Rate, QGAS,MIN(i) ) 1000 obj val (B/day)


no. iter

no. nodes

CPU (s)


62426.1 62426.1 62426.1 62426.1

CM LSM CHM SSM

75268.0 75268.0 75268.0 75268.0

CM LSM CHM SSM

82407.7 82407.7 82407.7 82407.7

CM LSM CHM SSM

151224.2 151224.2 151224.2 151224.2

62426.1 62426.1 62426.1 62426.1

147 752 113 92

0 0 0 0

0.04 0.10 0.04 0.05

4 4 2 2

0.16 0.22 0.22 0.10

9006 21 59 97

25.00 0.32 0.39 0.39

59 154 1513 29

0.28 0.69 3.91 0.13

2119 1403 4717 237

4.99 4.21 10.45 0.35


227 907 211 142


41418 988 399 832

6103

than 20 nodes are needed to reach global optimality. In every case, the special structure method visits the fewest number of nodes and takes least CPU time. 3.3. Gas Rate Lower Bound. To explore other possible constraints on the problem, this section explores the addition of lower bounds on the gas flow rate. A new binary variable, zi ∈ {0, 1}, ∀ i ) 1,..., NWELL, activates when QGAS,i * 0. The additional constraints are QGAS,i g QGAS,MINzi, ∀ i ) 1,..., NWELL. Tables 4-6 give solution data for the 200-well case with lower bounds QGAS,MIN ) 100, QGAS,MIN ) 500, and QGAS,MIN ) 1000. As expected from the additional constraints, the objective values in Tables 4-6 are somewhat lower than those reported in Table 3. Also, the computational requirements grow because of the additional binary variables. As the supply limit grows, the four formulations often require more computational time. All four frameworks are somewhat similar in their performance, but notice that both the classic method and the convex hull method occasionally function less efficiently than the other two. In Table 6, the classic formulation needs 25 s to complete (at supply limit 10 000 Mscf/day and QGAS,MIN(i) ) 1000). As this is the algorithm used in Kosmidis et al.,6 we suggest that their frequent calls to a piecewise linear algorithm would be more robust if they used another method. The special structure method from Keha et al.5 consistently out-performs the three other formulations. 4. Conclusion This paper focused on a comparative study of four modeling approaches for addressing the gas lifting problem. Each of the four formulations is sufficient to solve the problem to global optimality. However, the tests reported in this paper reveal that the classic formulation may be the least efficient of the four. We recommend using one of the other methods for challenging gas lifting problems that employ piecewise linear representations. It is also noteworthy to point out that given the power of the four formulations, piecewise linearization techniques may be used to optimize other functions, such as the piecewise underestimators recently proposed by Wicaksono and Karimi29 for bilinear problems. Piecewise linearization techniques may also be used for applications that involve bilinear monomials, trilinear monomials, and higher order nonlinearities that appear in separation synthesis, heat exchanger networks, reactor networks, and design and scheduling.30-34 Acknowledgment The authors gratefully acknowledge support from the National Science Foundation. R.M. is thankful for her National Science Foundation Graduate Research Fellowship.


478 1568 3265 304

Supporting Information Available: 200-Well test case described in section 3.2 as a supplement to our study. This material is available free of charge via the Internet at http:// pubs.acs.org.


173389.3 173389.3 173389.3 173389.3

173822.2 174060.5 173822.2 173822.2

4717 4745 8706 661

done by either using or avoiding binary variables yi,k, the CM, CHM, and SSM linear programming relaxations should be the same. The algorithm performance presented in Table 3 is representative of a whole range of supply limits. Notice that no more

Literature Cited (1) Buitrago, S.; Rodrı´guez, E.; Espin, D. Global Optimization Techniques in Gas Allocation for Continuous Flow Gas Lift Systems. Presented at the SPE Gas Technology Symposium, Calgary, Alberta, Canada, 1996; SPE Number 35616. (2) Nemhauser, G. L.; Woolsey, L. A. Integer and Combinatorial Optimization; J. Wiley: New York, 1988. (3) Floudas, C. A. Nonlinear and Mixed-Integer Optimization: Fundamentals and Applications; Oxford University Press: New York, 1995.

6104


(4) Sherali, H. D. On mixed-integer zero-one representations for separable lower-semicontinuous piecewise-linear functions. Oper. Res. Lett. 2001, 28, 155. (5) Keha, A. B.; de Farias, I. R., Jr.; Nemhauser, G. L. Models for representing piecewise linear cost functions. Oper. Res. Lett. 2004, 32, 44. (6) Kosmidis, V. D.; Perkins, J. D.; Pistikopoulos, E. N. A mixed integer optimization formulation for the well scheduling problem on petroleum fields. Comput. Chem. Eng. 2005, 29, 1523. (7) Bahadori, A.; Ayatollahi, Sh.; Moshfeghian, M. Simulation and Optimization of Continuous Gas Lift in Aghajari Oil Field. Presented at the SPE Gas Technology Symposium, Kuala Lumpur, Malaysia, 2001; SPE Number 72169. (8) Camponogara, E.; Nakashima, P. Optimizing gas-lift production of oil wells: piecewise linear formulation and computational analysis. IIE Trans. 2006, 38, 173. (9) Kosmidis, V. D.; Perkins, J. D.; Pistikopoulos, E. N. Optimization of Well Oil Rate Allocations in Petroleum Fields. Ind. Eng. Chem. Res. 2004, 43, 3513. (10) Mason, T. L.; Bagirov, A.; Klunder, J.; van Berkel, J. Development and Application of the Discrete Gradient Method for Non-Smooth Nonconvex Production Optimization. Presented at the International Oil Conference and Exhibition, Veracruz, Mexico, 2007; SPE Number 108496. (11) Lin, X.; Floudas, C. A. A Novel Continuous-Time Modeling and Optimization Framework for Well Platform Planning Problems. Optim. Eng. 2003, 4, 1–2. (12) Saputelli, L.; Nikolaou, M.; Economides, M. J. Self-Learning Reservoir Management Presented at the SPE Annual Technical Conference and Exhibition, Denver, CO, 2003; SPE Number 84064. (13) Saputelli, L.; Nikolaou, M.; Economides, M. Real-time reservoir management: A multiscale adaptive optimization and control approach. Comput. Geosci. 2006, 10 (1), 61–96. (14) Goel, V.; Grossmann, I. E. A stochastic programming approach to planning of offshore gas field developments under uncertainty in reserves. Comput. Chem. Eng. 2004, 28 (8), 1409–1429. (15) Goel, V.; Grossmann, I. E.; El-Bakryb, A. S.; Mulkay, E. L. A novel branch and bound algorithm for optimal development of gas fields under uncertainty in reserves. Comput. Chem. Eng. 2006, 30 (6-7), 1076– 1092. (16) Barragań-Hernańdez, V.; Va´zquez-Romań, R.; Rosales-Marines, L.; Garcıá-Sańchez, F. A strategy for simulation and optimization of gas and oil production. Comput. Chem. Eng. 2005, 30 (2), 215–227. (17) Kanu, E. P.; Mach, J.; Brown, K. E. Economic Approach to Oil Production and Gas Allocation in Continuous Gas Lift. J. Pet. Technol. 1981, 33, 1887. (18) Nishikiori, N.; Redner, R. A.; Doty, D. R.; Schmidt, Z. An Improved Method for Gas Lift Allocation Optimization. J. Energy Resourc. Technol. 1995, 117 (2), 87–92. (19) Fang, W. Y.; Lo, K. K. A Generalized Well-Management Scheme for Reservoir Simulation. SPE ReserV. Eng. 1996, 11 (2), 116–120; SPE Number 29124.

(20) Mason, T. L.; Emelle, C.; van Berkel, J.; Bagirov, A. M.; Kampas, F.; Pinte´r, J. D. Integrated Production System Optimization Using Global Optimization Techniques. J. Ind. Manage. Optim. 2007, 3, 257. (21) Camponogara, E.; Nakashima, P. Solving a gas-lift optimization problem by dynamic programming. Eur. J. Oper. Res. 2006, 174, 1220. (22) Camponogara, E.; Nakashima, P. Optimal Allocation of Lift-Gas Rates Under Multiple Facility Constraints: A Mixed Integer Linear Programming Approach. J. Energy Resourc. Technol. 2006, 128, 280. (23) Nakashima, P.; Camponogara, E. Optimization of Lift-Gas Allocation Using Dynamic Programming. IEEE Trans. Syst., Man., Cybernetics, Part A. 2006, 36, 407. (24) Keha, A. B.; de Farias, I. R., Jr.; Nemhauser, G. L. A Branch-andCut Algorithm Without Binary Variables for Nonconvex Piecewise Linear Optimization. Oper. Res. 2006, 54, 847. (25) Beale, E. M. L.; Tomlin, J. A. Special Facilities in a General Mathematical Programming System for Non-convex Problems Using Ordered Sets of Variables. In Proceedings of the Fifth International Conference on Operational Research; Lawrence, J., Ed: Tavistock Publishing: London, 1970; pp 447-454. (26) Tomlin, J. A. Special ordered sets and an application to gas supply operations planning. Math. Prog. 1988, 42:1-3, 69–84. (27) Brooke, A.; Kendrick, D.; Meeraus, A. GAMS: A User’s Guide; GAMS Development Corporation: Washington, D.C., 2005. (28) ILOG CPLEX 9.0.2 User’s Manual; ILOG, INC.; Mountain View, CA, Tavistock Publishing: London, 2005. (29) Wicaksono, D. S.; Karimi, I. A. Piecewise MILP Under- and Overestimators for Global Optimization of Bilinear Programs. AIChE J. 2008, 54, 991. (30) Floudas, C. A.; Grossmann, I. E. Synthesis of Flexible Heat Exchanger Networks With Uncertain Flowrates and Temperatures. Comput. Chem. Eng. 1987, 11:4, 319–336. (31) Floudas, C. A.; Anastasiadis, S. H. Synthesis of General Distillation Sequences with Several Multicomponent Feeds and Products. Chem. Eng. Sci. 1988, 43 (9), 2407–2419. (32) Ciric, A. R.; Floudas, C. A. A Retrofit Approach of Heat Exchanger Networks. Comput. Chem. Eng. 1989, 13 (6), 703–715. (33) Kokossis, A. C.; Floudas, C. A. Optimization of Complex Reactor Networks-II: Nonisothermal Operation. Chem. Eng. Sci. 1994, 49 (7), 1037– 1051. (34) Lin, X.; Floudas, C. A. Design, Synthesis and Scheduling of Multipurpose Batch Plants via an Effective Continuous-Time Formulation. Comput. Chem. Eng. 2001, 25, 665–674.

ReceiVed for reView August 6, 2008 ReVised manuscript receiVed November 7, 2008 Accepted November 10, 2008 IE8012117

Global Optimization of Gas Lifting Operations: A Comparative Study of

Recommend Documents