Optimal Allocation Procedure for Gas-Lift Optimization - American

Jan 26, 2010 - based on an upper convex profile requirement reduces a high dimensional problem into ... rise to the nonlinear gas-lift optimization pr...
1 downloads 0 Views 1MB Size
2286

Ind. Eng. Chem. Res. 2010, 49, 2286–2294

Optimal Allocation Procedure for Gas-Lift Optimization Kashif Rashid* Schlumberger-Doll Research, 1 Hampshire Street, Cambridge, Massachusetts 02139

This paper presents an optimal allocation procedure (referred to as the Newton reduction method) which based on an upper convex profile requirement reduces a high dimensional problem into one of a single variable. This allocation procedure is used in a methodology designed for the treatment of the gas-lift optimization problem. This problem is encountered in the oil industry and specifically concerns the distribution of a fixed quantity of lift gas over a number of wells in order to increase the field production capability. The allocation procedure is applied by enforcing the notion of well separability and the availability of upper convex lift profiles for each well. However, due to the effects of interdependent wells, the procedure is performed iteratively until convergence on wellhead pressures. The approach reported is shown to be several times more efficient than traditional approaches in computational and time costs, while delivering comparable results. Introduction As a producing oil field matures, the declining reservoir pressures from continued hydrocarbon extraction make oil production from existing and new wells harder. To alleviate this problem in part, natural gas is often injected at high pressure from the casing into the well bore. This method of artificial lift is known as gas lift. As it is relatively inexpensive, easy to implement, and applicable over a broad range of conditions, it is a favored method of lift in many operating fields and not only because some or all of the gas produced by a field can be used as the source of lift gas.1 When natural gas is injected at high pressure into the well bore, it mixes with the produced fluids from the reservoir (see Figure 1), reducing the density of the fluid column and effectively lowering the bottom-hole pressure (Pbh). The increased pressure differential induced across the sandface (the connection point between the reservoir and the well, given by Pr - Pbh, where Pr is the in situ reservoir pressure) allows more fluid to flow to the surface. However, too much lift gas increases the frictional pressure drop and reduces the fluid production possible. Hence, although each well has a desirable lift-gas quantity, when the entire gathering network is considered, an optimal distribution must be made to account for the backpressure effects imposed by interconnected wells. This gives rise to the nonlinear gas-lift optimization problem. The broader production problem can include the activation state of wells and the control of subsurface chokes, among other network elements. In addition, the numerous parameters involved in the gas-lift design problem (i.e., choice of gas-lift valve size, number, and depths) may also be considered, but are often excluded from the optimization problem as the well configuration is largely fixed in a mature field. That is, lift-gas injection normally takes place at the deepest valve at the highest possible injection pressure given the well configuration and operating conditions. Lastly, well stability is assumed to be integrable; otherwise the optimization procedure would be limited and premature. In this work, only the continuous gaslift optimization problem is considered.

to separator constraints. Under these, and other operating constraints, it is necessary for engineers to optimally allocate the available lift gas to maximize the field oil production, revenue, or indeed profit. In order to do so, it is common practice to model the physical system using a multiphase flow simulator with data collected at the well site. The ensuing model is used for optimization purposes, and if the model is an accurate representation of the physical system, the optimal configuration can be applied directly to the real system, either manually or automatically in a closed loop. A gas-lift network model in a steady-state multiphase flow simulator comprises a description of the gathering network, the well configurations, the pressures or flow rates at boundary

Problem Definition and Prior Art A gas-lifted field is constrained by the amount of gas available for injection or additionally the produced gas permissible due * E-mail: [email protected].

Figure 1. Gas-lift well schematic.

10.1021/ie900867r  2010 American Chemical Society Published on Web 01/26/2010

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

conditions, the composition of the produced fluid in each well, the multiphase flow correlations employed, and the quantity of lift gas injected into each well.2 The last can be considered control Variables, while the preceding elements can be deemed as constant network parameters (F) with respect to a gas-lift optimization scenario. For a network with N wells, the intent is to optimally allocate a fixed amount of gas C, such that the oil production at the sink node (F) is maximized. Mathematically, this can be stated as max Fnw ) fonline(L;F)

(1)

N

s.t.

∑L

i

eC

i)1

L ∈ RN where L is the vector of gas-lift rates and Li is the ith component of L. The problem defined by eq 1 is a nonlinear constrained optimization problem in which each function evaluation requires a call to the network simulator. In the context of the present methodology, this is referred to as the online problem. The direct approach is to solve eq 1 using a suitable nonlinear constrained solver. This includes the use of standard deterministic gradient-based methods such as sequential quadratic programming (SQP)3 and the augmented Lagrangian method (ALM).4 These methods tend to be computationally expensive if numerical derivatives are required. However, they offer fast convergence from good starting conditions. Under linearization assumptions (of the objective function), the sequential linear programming (SLP) method can be employed.5 Derivative-free methods, including the polytope method, the genetic algorithm (GA), and other heuristic-based schemes, can also be used.6-8 Treatment of the well-rate and the gas-lift allocation problem together, in a broader definition of eq 1, was demonstrated using a mixed-integer nonlinear programming (MINLP)9 and a dynamic programming (DP) approach.10 Evidently, as each function evaluation with the direct approach is a call to the underlying network simulator, these approaches can be time-consuming and computationally costly, especially if the number of variables is great, numerical derivatives are required, and the simulation is expensive to run, as is often the case. However, as the network model performs a rigorous pressure and flow rate balance, the benefit is that a steady-state solution is returned, in contrast to methods in which the interaction of interconnected wells is neglected. If the wells are considered independently of one another, for given wellhead pressure conditions, the problem in eq 1 can be defined by the separable program: max Fsep ) foffline(L;P)

(2)

N

s.t.

∑L

i

max Fsep )

∑ Q (L ) i

i

(3)

i)1

N

s.t.

∑L

i

eC

i)1

L, Q ∈ RN where L and Q are the vectors of gas-lift rates and flow rates, respectively, and Li and Qi are the ith components of L and Q, respectively. The lift profiles for each well can be obtained from actual well step-rate tests conducted at the well site or from single well nodal analysis calculation. While the former is likely to be more accurate and representative of the actual behavior observed, the latter is more practical for fields with many wells and can also provide a family of curves with varying wellhead pressure. In the present work, eq 2 is referred to as the offline problem. Early methods developed for the gas-lift optimization problem arose from the solution of the separable problem eq 3. They employed gas-lift performance curves to describe the well behavior but neglected to treat the effect of gas injection on wellhead pressure, in addition to the back-pressure effects imposed by interconnected wells. Hence, such methods give only pseudo-steady-state solutions, though, in general, this limitation can be overcome by repeating the procedure with updated well properties. The equal-slope concept dictates that the sensitivity of well flow rate to a change in the gas injected should be the same for all wells. This infers that a reduction and redistribution of lift gas from any one well to another will not increase the total production rate. Initially, a solution yielding zero slopes was sought assuming an unlimited supply of lift gas. However, a better solution was defined as one in which the incremental gain is equal to the incremental cost of injection.11-15 If consistent revenue and cost factors are employed, a slope solution of 1,14 the economic slope,15 the point of zero incremental gain,11,12 and the zero-slope solution for maximum profit all give rise to the same solution (see Figure 3). In addition, as the reciprocal of the slope of the GLPC gives the incremental GOR (IGOR), methods based on establishing equal IGOR (the marginal GOR) are also based on the equal-slope concept.16,17 A number of methods utilizing the equal-slope concept have been developed.11,12,14,15,18-22 Largely heuristic in nature, these methods discretize the available lift gas and allocate it incrementally to the wells with favorable production gradients until an equal-slope solution is obtained. The optimal allocation method presented in this paper, the Newton reduction method (NRM), also returns an equal-slope solution. The use of standard NLP solvers has also been reported.23-26 With second-order polynomials used to fit the GLPC, a quasi-

eC

i)1

L, P ∈ RN where L and P are the vectors of gas-lift rates and wellhead pressures, respectively, and Li and Pi are the ith components of L and P, respectively. Here, the flow rate versus gas injection profile, known as the gas-lift performance curve (GLPC), is defined for each well (see Figure 2) and the objective function is given as the sum of all well flow rates. Hence eq 2 can be rewritten as

2287

N

Figure 2. Gas-lift performance curve at wellhead pressure (P).

2288

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

Figure 5. Representative family of lift performance curves.

Figure 3. Performance curve representation and solution schemes.

Figure 4. Gas-lift performance curves at wellhead pressure (P).

Newton method23 and the Lagrange multiplier method with a standard NLP solver were demonstrated.24,25 Although the derivative requirement increased the computational cost, the methods can exhibit superlinear convergence from good starting conditions. However, loss of accuracy is possible due to the quality of the fitted polynomial. To overcome this, a modified second-order polynomial was used to better fit the data obtained from nodal analysis with an SQP solver based on an approximate Hessian update procedure.26 SLP techniques were proposed27,28 assuming piecewise linear lift performance curves and linearization of constraints using first-order Taylor series expansion. The approach can be fast, as only first-order derivatives are required, but it has the limitation that the linear model may be a poor representation of the nonlinear system. Wells with nonsmooth gas-lift performance curves, with discontinuities, and indeed with noninstantaneous flow (NIF) can impede the progress of conventional gradient-based optimization methods (see Figure 4). Due to poor handling of such wells and simplifications made by curve fitting, suboptimal solutions are often obtained. For example, as the flat section of an NIF well returns a gradient that is of little use in the optimization procedure, a minimum lift-gas injection constraint is often implemented. However, this requirement leads to an inefficient allocation and a suboptimal solution. Conventional methods also tend to find only local solutions. To overcome these issues, the use of robust global search algorithms were proposed, including heuristic-based schemes and genetic algorithms (GAs).29-31 The heuristic scheme employed random search along with a clustering method29 and is reportedly able

to handle wells with irregular profiles. A simple GA was reported in ref 31, and a multiobjective scheme to additionally minimize lift gas use was reported in ref 30. The latter employed a variant of the “elitist” multiobjective NSGA-II algorithm32 in order to maintain population diversity, but assumed piecewise linear lift curves. The limitation of curve-based schemes for problem 2 is that the solution is only pseudo steady state, while the schemes for problem 1 can be computationally expensive. A means to overcome the latter is to employ a proxy of the real networksimulation-based objective function during optimization. The use of a trained neural network (NN) in conjunction with a GA has been demonstrated.33 This can be considered an indirect approach, and can help alleviate the computational cost to some extent. The work presented in this paper comprises an iterative offline-online procedure, which solves eq 2 and uses the optimal solution determined to call eq 1 for updated well properties. The offline problem merits a fast solution, while the online evaluation ensures fidelity to the original problem. Notably, results comparable to direct optimization are obtained in only a fraction of the number of function evaluations. Optimal Allocation Procedure Consider the availability of upper convex lift performance curves for each well in the network model, describing the relationship of oil flow rate with the amount of lift-gas injection for varying wellhead pressure (see Figure 5). These can be extracted from a network simulator (on a well-by-well basis) using nodal analysis sensitivity studies.2 The curves are nonlinear and show that too much gas injection has a detrimental effect on the oil production rate. The reduced density of the mixture in the well bore increases initial production, but as the gas injection rate increases, frictional pressure losses become more significant, causing the production to peak and then start declining. Note that the use of liquid-rate curves is appropriate only in the case where each well has the same water-cut. If the wells in the network model are assumed decoupled from one another, the original problem, eq 1, can be restated as problem 2, and specifically by problem 3, in which the behavior of each well is described by a single lift curve at a specific wellhead pressure, given by Qi ) f(Li ;Pi)

(4)

If it is assumed that all the available lift gas is used, the constrained problem eq 3 can be converted to an unconstrained problem by formation of a Lagrangian function including a penalty factor (λ ∈ R):

min M(L, λ) ) -Fsep + λ

(

N

∑L

i

-C

i)1

)

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

Li ) gi(λ;Pi)

(5)

Imposing Karush-Kuhn-Tucker (KKT) conditions for optimality, the following expressions are obtained, in which notably ∂Qi/∂Lj ) 0 for all i * j and i, j ∈ [1, ..., N]:

2289

(10)

It is insightful to note that the superimposition and summation of all inverse derivative profiles gives a quantity with the same units as the allocated quantity (i.e., the unit of lift gas). This quantity is defined by the following: N

dQi ∂M )+λ)0 ∂Li dLi ∂M ) ∂λ

(6)

N



Li - C ) 0

(7)

i)1

E(λ) is constrained by the total gas available for injection, C, so in practice E(λ) e C. However, if we treat C with a strict equality constraint, under the assumption that all the available lift gas is utilized, a residual function can be composed: R(λ) ) E(λ) - C

(8)

N

∑L

i

)C

(9)

Equation 9 simply treats the allocated lift gas as an equality with respect to the available lift gas, while eq 8 suggests that the slopes of the operating curves for each of the wells have the same value, λ, as indicated by the equal-slope concept. Observing that λ merely indicates the level of each derivative profile (see Figure 6), it is evidently bound between the highest and lowest possible derivative values (dQi/dLi) over all wells and must also satisfy eq 9. Typically, once a well is flowing, the lift performance curve is concave and the ensuing derivative profile is convex. Hence, it is possible to establish the inverse of the derivative profile of each well. This is shown in Figure 7 for one particular operating curve and is given mathematically in eq 10.

(12)

Hence, the optimal value of λ is obtained for the solution of eq 12 with R(λ) ) 0 and can be achieved using Newton’s method (eq 13). λnew ) λold -

i)1

Figure 6. Lift performance curve and derivative profile.

(11)

i

i)1

From eqs 6 and 7 we obtain eqs 8 and 9, respectively: dQi )λ dLi

∑L

E(λ) )

R(λ) R'(λ)

(13)

)

(14)

where R(λ) )

(

N

∑ g (λ;P ) i

i

i)1

-C

and R'(λ) )

dR ) dλ

∑ ( dλ (λ;P )) N

dgi

i

(15)

i)1

Hence, the original N-dimensional inequality constrained nonlinear problem is converted to a single variable equality constrained problem that can be solved using a classic rootfinding technique. This is referred to as the Newton reduction method (NRM).34,35 Note that the obtained solution, satisfying the KKT conditions, is the global solution due to the convexity assumptions imposed.4 Although the addition of constraints can complicate matters, primary, secondary, and, to a lesser extent, tertiary level constraints can be accommodated. Primary constraints concern elements of the operating curve (i.e., GLIR and oil rate for each well) and are managed using penalty imposition34 (see Appendix). Secondary constraints (e.g., liquid rate or GOR limits) concern those factors which are functions of the control variables allowing imposition of the equivalent primarily level constraints. Tertiary constraints, such as produced gas limits, can act at the manifold or sink level. Manifold constraints are not easily managed unless upstream wells can be identified and an offline representation is possible.36 Sink constraints, however, can be managed iteratively. That is, the optimal quantity of lift gas to allocate is established that ensures the sink level constraints are met. In general, the procedure works best for the allocation problem alone, or with primary and secondary constraints only.37 Finally, as the wells are effectively decoupled in the offline problem, well interaction is accounted for by a call to the online problem. This leads to the overall methodology for the gas-lift optimization problem. Methodology for Gas-Lift Optimization

Figure 7. Inverse derivative profile.

The methodology developed uses an iterative offline-online approach and embodies the optimal allocation procedure of the

2290

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

Table 1. Example of GLO Methodology Using a Two-Well Network Modela offline online err new iteration P1, psia P2, psia L1, MMscf/day L2, MMscf/day Q1, bbl/day Q2, bbl/day Foffline, bbl/day Pnew , psia P 1 2 , psia Fonline, bbl/day |P|∞, psia 1 2 3

250.00 156.30 156.42

250.00 156.47 156.47

1.1316 1.0923 1.0926

0.8687 0.9078 0.9076

1440.7 1666.2 1666.0

969.5 1168.6 1168.6

2410.2 2834.8 2834.6

156.30 156.42 156.42

156.47 156.47 156.47

2836.2 2837.2 2837.2

93.70 0.12 0.00

a Three iterations are required to solve the problem. C ) 2 MMscf/day and 250 psia is the default initial wellhead pressure for both wells. The same result is obtained with other starting pressures. Comparative solutions from other approaches are shown in the Appendix for this case.

last section. Online refers to a simulation run of the actual network model, and offline refers to the solution of the decoupled problem eq 3 employing the lift performance curves shown in Figure 5 in place of the network simulator. First, in a preprocessing step, the family of lift performance curves are extracted for each well from the network model2 and validated for use. In particular, it is necessary to ensure that the inverse derivative profiles are monotonically decreasing. This step must be performed only once and the time cost is proportional to the number of control wells in the network model. In general, the time taken to generate the family of curves for each well depends on the range and number of pressure and gas injection values selected for the sensitivity study. A considered choice must be made to cover the anticipated operating range of interest with sufficient detail, while minimizing the total expected preprocessing cost. Once completed, however, there can be a significant advantage in time-cost reduction compared to other approaches if the optimization problem is re-run with a change in available lift gas, operating conditions, or constraints. Thus, the real cost is considered to be the number of simulation calls necessary to obtain a solution. Once the preprocessing step is completed, the vector of wellhead pressures (P) is initialized. This dictates the current operating curve for each well. At subsequent iterations the updated wellhead pressure estimate is used to set the current operating curve. If the desired wellhead pressure does not match the family of curves stored, it is generated by interpolation. Additionally, each performance curve is described using a spline function, while ensuring the convexity of the inverse derivative profile. When the operating curves have been set, the lift gas is ˆ ) among the N wells according to the total optimally allocated (L lift gas constraint C, to maximize the objective function eq 3 using the optimal allocation procedure presented. Note that the offline constrained nonlinear problem may be solved using alternative approaches and indeed this is necessary in the presence of multiple midnetwork constraints.37 Once the allocation procedure is completed, the optimal set of gas-lift rates are applied to the real network model. The oil production rate from the online simulation run can be compared to the offline solution, though primarily the new wellhead pressures are sought (Pnew). Convergence of the offline-online procedure is assessed using the L-infinity (maximum absolute difference) norm between the “old” and “new” estimates of the wellhead pressure vector (eq 16). ε ) |Pnew - P| ∞

(16)

This metric proves more effective in practice than other norms as it considers only the largest pressure differential at each iteration. In comparison, the Euclidean norm gives rise to an average estimate. If the convergence test is within a specified ˆ , Pˆ, tolerance, typically 0.5 psia, the optimal results (λˆ , Lˆ, Q

and Fˆ) are returned. Otherwise, the procedure repeats with the updated pressure vector, Pnew, used to set the current operating curves. The methodology is summarized in algorithmic form in eq 17. step 1:

k)1

step 2: P k ⇐ Pini step 3: L ⇐ foffline(λ;P k) step 4: P k+1 ⇐ fonline(L;G) step 5: ε ) |P k+1 - P k | ∞ step 6: ε > εtols ;

P k ⇐ P k+1 ;

k ⇐ k + 1;

return step 3

step 7: ε e εtols ;

ˆ , Pˆ, Fˆ ; return λˆ , Lˆ, Q

stop

(17)

An example of the application of the methodology described above is presented in Table 1 for a network model comprising two gas-lifted wells with 2 MMscf/day of available lift gas. Three iterations of the offline-online procedure (steps 3-6) are required until convergence is achieved, given by a near zero pressure norm. It is worth noting that the norm on the objective function values (|F|∞, not given in Table 1) is nonzero and indicates the implicit model error present in the offline curvebased representation. Collectively, a solution is obtained with only three network simulation evaluations. Direct optimization with a derivative-free polytope method required 20 function evaluations to give the same solution (see the Appendix). Once a solution has been obtained, as a consequence of the original inequality constraint, the sensitivity of the objective function must be tested to assess true optimality. This is also true if an offline solution is obtained alone (see Figure 8). If the derivative of the total production versus available lift-gas profile (dF/dC) is negative when allocating all the available lift gas (Cmax), then the optimal lift gas to allocate (Copt) must be identified, or in the presence of permissible produced gas constraint, the highest amount of lift gas that can be allocated (Ccon) without violating the constraint must be identified. If the plot of the total production versus available lift gas is established a priori and stored along with the corresponding allocated lift rates, the optimal configuration can be obtained immediately with changing gas-lift availability, assuming no production capacity constraints are violated, and used for realtime operation. For example, if C g Copt, the solution for Copt would be applied and, if C < Copt, the solution for optimally

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

2291

Qo(bbl/day) Qg(Mscf/day)

(18)

η)

Figure 8. Ensuring solution optimality.

allocating C would be applied. With a produced gas constraint, the solution for Ccon would be used while C g Ccon (assuming Ccon < Copt) and the solution for allocating C optimally otherwise (see Figure 8). Suitable consideration can be made in the case with Ccon > Copt. For significant changes in operating conditions, including individual well operability, imposition of new constraints, or changes in well properties, the optimization problem would need to be re-run, with updated GLPC data, if necessary. 56-Well Curve-Based Model In Buitrago et al.29 a 56-well curve-based test case, comprising 10 noninstantaneously flowing (NIF) wells, was presented with 22 500 Mscf/day of lift gas available for allocation. This example is used to demonstrate the allocation procedure and also the treatment of NIF wells. NIF wells have a minimum gas injection requirement before the well is able to produce (see Figure 4). If the level of gas injection is less than this minimum amount, the well remains dead. As NIF wells can complicate many optimization procedures, they are often forced to be on, under the assumption that a producing well is desirable, or simply switched off by default. The former results in a large amount of lift gas being used simply to keep the well aliVe, while the latter ignores the production potential of the well completely. For example, in this case, a total production of 20 366.8 bbl/day is obtained with all NIF wells on and a production of 22 015.1 bbl/day if all NIF wells are off. Hence, while not optimal, switching the NIF wells off is preferable to forcing the wells to produce. As an optimal solution will depend on some combination of NIF well activation, the state of such wells must be considered, though it is often neglected. In this work, rather than resorting to the use of a discrete optimizer, a rank-and-shut scheme is adopted. First, the requirement that NIF wells produce (i.e., implementing minimum injection constraints) is relaxed. Consequently, this can give rise to a solution in which certain wells may receive gas but produce nothing, i.e., NIF well operation in the dead zone. This is made possible by extension of the inverse derivative profiles to ensure that a solution is always possible, much like the procedure for constraint handling shown in the Appendix. Second, a metric for well-lift efficiency, η (eq 18), is introduced and is defined as the ratio of oil production to the amount of gas injected:

This metric, the level of oil production per unit lift-gas injection, infers the value of allocating a quantity of gas to any particular well. By imposing a minimum threshold production efficiency, ηmin, wells with an inefficient production solution can be ranked for shutdown. Thus, starting with the default case, in which all the NIF wells are on, the optimal allocation procedure is performed and the lift efficiency of each NIF well is determined at the solution. The NIF well with the lowest lift efficiency is shut down. The procedure of allocation and well closure is repeated until shutdown of a particular NIF well commences to decrease the total production achieved. In Table 2 solutions are given progressing from the case with all NIF wells on to all switched off. In practice, the procedure would be halted at iteration 8 when a reduction in the overall objective value is observed (see Figure 9). Note also that well 50 receives no lift gas at iteration 1. However, at iteration 2, once well 55 has been shut down, it flows with a production efficiency of 386 stb/Mscf. The best solution obtained using the rank-and-shut procedure appears in Table 2 at iteration 7, where the blank cells indicate closed NIF wells and all non-NIF wells are open. This solution was validated using an exhaustive search of all possible combinations (210 ) 1024). As the solution can depend on the starting condition it is not, in general, provably optimal. However, as shown, starting from the default assumption that all NIF wells are on, a good solution is obtained in less than n + 1 iterations (where n is the number of NIF wells). In Table 3, the optimal NRM result is given alongside the cases with all NIF wells switched on and off. In addition, solutions from a global derivative-free search method (EXIN),29 an equal-slope-based scheme (EQSL),29 and a multiobjective evolutionary approach (RAY)30 are given. As the solutions reported vary due to the choice of interpolation method and the solution technique employed, for comparative purposes the last row of Table 3 provides the equivalent values from the present model. Notably, the optimal NRM solution has the highest objective value. It is also worth noting that the EXIN solution is with all NIF wells off, while the RAY solution is with all NIF wells on. Hence, although both methods reportedly handle NIF wells, neither has proven particularly effective or optimal. As curve-based schemes can give only pseudo-steady-state solutions, the allocation procedure presented here is combined with the update procedure for network simulation models, as per the gas-lift optimization methodology presented. 100-Well Network-Based Model In this section, the proposed offline-online method, referred to as GLO,38 is demonstrated on a large gathering network of 100 wells, comprising 10 submanifolds, each gathering from 10 production wells (see Figure 10) with 80 MMscf/day of lift gas available for allocation. As the emphasis is on the comparative cost and quality of solution obtained, for convenience, details of network topology, well configuration, and fluid compositions are omitted. However, in general, each well uses a black-oil definition, with varying water-cut and gas-oil ratio value. These properties are exemplified by the lift performance curves obtained during the preprocessing stage.

2292

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

Table 2. NIF Iterations for 56-Well Curve-Based Model with C ) 22 500 Mscf/daya iteration NIF wells ηwell55 ηwell50 ηwell51 ηwell52 ηwell54 ηwell47 ηwell49 ηwell53 ηwell48 ηwell56 shut well model (bbl/day)

1 10 324 n/a 510 590 542 767 859 832 631 1303 55 20 367

2 9

3 8

4 7

5 6

6 5

7 4

8 3

9 2

10 1

11 0

386 521 600 570 772 867 913 776 1359 50 20 918

528 600 596 774 872 947 977 1370 51 21 464

600 603 776 875 966 1070 1375 52 22 084

606 777 877 973 1093 1375 54 22 502

777 877 973 1090 1371 47 22 667

877 972 1086 1345 49 22 687

972 1085 1334 53 22 648

1083 1312 48 22 552

1271 56 22 318

n/a 22 015

a

The η metric is in stb/Mscf, and the blank cells indicate closed NIF wells. The optimal solution (22 686.6 bbl/day) appears at iteration 7, for an allocation of 22 499.9 Mscf/day. Non-NIF well closure is not considered, although the scheme could be extended to include all wells.

Figure 9. Production variation with NIF well closure. Table 3. 56-Well Curve-Based Modela method

total gas (Mscf/day) total oil (bbl/day) model (bbl/day)

NRM (on)

NRM (off)

NRM (opt)

EXIN29

EQSL29

RAY30

22 500

22 500

22 500

20 479

24 661

22 376

20 367

22 015

22 687

21 790

21 265

22 033

20 367

22 015

22 687

20 655

20 844

22 140

a In Table 3 of the Buitrago paper,29 the 56-well Qg sums for EXIN and EQSL are incorrectly stated as 20 453.9 and 22 508 Mscf/day, respectively. The correct values (presented above) are 20 479 and 24 661 Mscf/day, respectively. The RAY scheme uses linear interpolation.

In Table 4 the GLO solution is compared to the solution obtained from a direct optimization approach using the downhill simplex method, referred to as the polytope method.39 The latter is selected for comparative purposes as an example of an efficient derivative-free direct solution method. Notably, from a uniformly distributed starting point, comparable objective function values are obtained. However, while the polytope method required 369 function evaluations (and 153 min), GLO required only eight calls to the actual network model. Also, in the latter case, 25 of the total 30 min required for solution was spent during the preprocessing stage. This is important to note, as this cost is necessary only once and with changing conditions (e.g., gas availability) the problem can be resolved quickly without having to repeat the preprocessing stage. This will only be necessary if the

Figure 10. A 100-well gathering network schematic. Key: o, wells; •, manifold nodes; |, branches. v is the transport line to the delivery sink O. The boundary conditions are all pressure specified. Table 4. 100-Well Network Resultsa method

Fonline (bbl/day) preprocessing time (min) solution time (min) total time (min) simulator calls % abs difference a

polytope

GLO

27 438 n/a 153.6 153.6 369 0

27 365 25.0 5.0 30.0 8 0.27

N ) 100 wells and C ) 40 MMscf/day.

network model conditions change appreciably, requiring updated lift performance curves to be generated. Furthermore, if the optimal solutions are obtained for varying quantities of lift gas a priori, the method can be readily adopted in a real-time closed-loop system. Note that the solution of realworld large-scale integrated networks has been demonstrated using this approach.37,40 Hence, the proposed offline-online procedure is computationally efficient, while providing

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

2293

significantly more efficient compared to other approaches, achieving comparable results in only a fraction of the number of simulator evaluations. Although the method can handle operating constraints using penalty imposition, it is somewhat limited with respect to manifold level constraints and dealing with complicated looped networks. However, the proposed method can be extended to accommodate energy (horsepower) allocation for ESP-lifted wells, methanol allocation for chemical stimulation, and the control of downhole chokes for pressure and rate management, i.e., a number of production optimization scenarios based on constrained resource allocation. Finally, the approach is considered a step forward in solving large-scale production optimization scenarios efficiently and can aid the development of real-time integrated asset solutions. Acknowledgment Thanks to Luca Letizia, Daniel Lucas-Clements, Balraj Grewal, Andy Shand, Trevor Tonkin, Scott Raphael, and Andy Howell, without whom this work and the development of GLO38 would not have been possible. I am grateful to William Bailey, Benoıˆt Coue¨t, and the anonymous reviewers for their comments on earlier versions of the manuscript. Appendix Handling Primary Constraints.

Figure 11. GLPC modification for primary constraint handling. Table 5. Network Model Test Resultsa method

polytope

proxy-polytope

GLO

A two-well example of GLPC modification for primary constraint handling is shown in Figure 11.

2-Well Case Fopt (bbl/day) simulator calls

2837.2 20

2836.2 14

2837.2 3

5750.0 18

5762.5 4

4-Well Case Fopt (bbl/day) simulator calls

5764.2 59 26-Well Case

Fopt (bbl/day) simulator calls

45 913 179

45 922 88

45 905 4

a

GLO is the proposed method. Polytope is the downhill simplex method. Proxy-polytope uses a neural-network proxy model with the polytope solver. C ) 2, 4, and 45 MMscf/day for 2, 4, and 26-well cases.

comparably good solutions. Results from other cases are shown in the Appendix. Conclusion A methodology for gas-lift optimization has been presented based on an iterative offline-online procedure. The actual online network model is replaced by an offline curve-based approximation by enforcing the notion of well separability. An optimal lift-gas allocation is achieved using the Newton reduction method (NRM), which converts the original nonlinear constrained problem into one of a single variable with a strict equality. At the final solution each well has the same sensitivity to an incremental gain in lift gas. A curve-based test demonstrated the efficacy of the method for the pseudo-steady-state solution and the treatment of noninstantaneously flowing wells using an iterative rank-andshut procedure. The steady-state solution was demonstrated on a network model of 100 wells. The proposed method is

Network Model Test Results. Table 5 lists results from the network model tests for the 2-, 4, and 26-well cases. Nomenclature Acronyms ALM ) augmented Lagrangian method DP ) dynamic programming ESP ) electric submersible pump EQSL ) method presented in Buitrago et al.29 EXIN ) method by Buitrago et al.29 GA ) genetic algorithm GL ) gas lift GLR ) gas-liquid ratio GLV ) gas-lift valve GLIR ) gas-lift injection rate GLPC ) gas-lift performance curve GLO ) gas-lift optimization methodology presented GOR ) gas-oil ratio IGOR ) incremental GOR KKT ) Karush-Kuhn-Tucker conditions for optimality MGOR ) marginal GOR MINLP ) mixed-integer nonlinear programming NIF ) noninstantaneous flow NN ) neural network NRM ) Newton reduction method NSGA ) nondominated sorting GA32 RAY ) reference for method by Ray and Sarker30 SLP ) sequential linear programming SQP ) sequential quadratic programming

2294

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

Symbols C ) available lift gas foffline ) function representing separable curve-based model fonline ) function representing network-based model Fnw ) oil production from network-based model (fonline) Fonline ) oil production from network-based model (fonline) Fsep ) oil production from curve-based model (foffline) Foffline ) oil production from curve-based model (foffline) L ) single well GLIR (used in figures) Li ) ith component of L L ) vector of gas-lift rates (L ∈ RN) N ) number of gas-lifted wells n ) number of NIF wells from N Pbh ) bottom-hole pressure Pr ) reservoir pressure Pwh ) wellhead pressure P ) single wellhead pressure (used in figures) Pi ) ith component of P P ) vector of wellhead pressures (P ∈ RN) Q ) single well oil flow rate (used in figures) Qo ) single well oil flow rate Qg ) single well gas flow rate Qw ) single well water flow rate Qi ) ith component of Q Q ) vector of oil flow rates (Q ∈ RN) RN ) N-dimensional search space

Literature Cited (1) Brown, K. E. Overview of Artificial Lift Systems. JPT, J. Pet. Technol. 1982, 34 (10), 2384–2396; 9979-PA. (2) PipeSim Network Simulator; Schlumberger SIS: Abingdon, Oxon, U.K., 2007. (3) Dutta-Roy, K.; Kattapuram, J. A New Approach to Gas-Lift Allocation Optimization. Presented at the SPE Western Regional Meeting, California, June 1997; paper no. 38333-MS. (4) Fletcher, R. Practical Methods of Optimization; John Wiley and Sons: New York, 2000. (5) Nadar, M. S.; Schneider, T. S.; Jackson, K. L.; McKie, C. J. N.; Hamid, J. Implementation of a Total System Production Optimization Model in a Complex Gas Lifted Offshore Operation. SPE Prod. Oper. 2008, 5– 13; 103670-PA. (6) Carroll, J. A.; Horne, R. N. Multivariate Optimization of Production Systems. SPE 1992, 44, 7; 22847-PA. (7) Palke, M.R.; Horne, R. N. Nonlinear Optimization of Well Production Considering Gas Lift and Phase Behavior. Presented at the SPE Production Operations Symposium, Oklahoma, March 1997; paper no. 37428-MS. (8) Vazquez, M.; Suarez, A.; Aponte, H.; Ocanto, L.; Fernandez, J. Global Optimization of Oil Production Systems: A Unified Operational View. Presented at SPE ATCE, Louisana, October 2001; paper no. 71561-MS. (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–3527. (10) Camponogara, E.; Nakashima, P. H. R. Solving a gas-lift optimization problem by dynamic programming. Eur. J. Oper. Res. 2006, 174 (2), 1220–1246. (11) Simmons, W. E. Optimizing Continuous Flow Gas Lift WellssPart 1. Pet. Eng. 1972, 45 (8), 46–48. (12) Simmons, W. E. Optimizing Continuous Flow Gas Lift WellssPart 2. Pet. Eng. 1972, 44 (10), 68–72. (13) Mayhill, T. D. Simplified Method for Gas-Lift Well Problem Identification and Diagnosis. Presented at the SPE Fall AIME Meeting, Texas, October 1974; paper no. 5151-MS. (14) Redden, J. D.; Sherman, T. A. G.; Blann, J. R. Optimizing GasLift Systems. SPE Fall AIME Meeting, Texas, October 1974; pp 1-13, paper no. 5150-MS. (15) Kanu, E. P.; Mach, J.; Brown, K. E. Economic Approach to Oil Production and Gas Allocation in Continuous Gas-Lift. SPE JPT 1981, 33 (10), 1887–1892; 9084-PA.

(16) Stoisits, R. F.; Scherer, P. W.; Schmidt, S. E. Gas Optimization at the Kuparak River Field. SPE ATCE, Louisana, September 1994; pp 3542, paper no. 28467-MS. (17) Weiss, J. L.; Masino, W. H.; Starley, G. P.; Bolling, J. D. Largescale Facility Expansion Evaluations at the Kuparuk River Field. Presented at the SPE Regional Meeting, California, April 1990; pp 297-304, paper no. 20046-MS. (18) Everitt, T. A. Gas-Lift Optimization in a Large GOM Field. SPE ATCE, Louisana, September 1994; pp 25-33, paper no. 28466-MS. (19) Schmidt, Z.; Doty, D. R.; Agena, B.; Liao, T.; Brown, K. E. New Developments to Improve Continuous-Flow Gas-Lift Utilizing Personal Computers. SPE ATCE, Louisana, September 1990; pp 615-630, paper no. 20677-MS. (20) Eclipse Reference Manual; Schlumberger SIS: Abingdon, Oxon, U.K., 2006. (21) Chia, Y. C.; Hussain, S. Gas Lift Optimization Efforts and Challenges. SPE Asia Improved Oil Recovery Conference, Malaysia, October 1999; pp 2-9, paper no. 57313-MS. (22) GOAL Gas-Lift Optimizer; Schlumberger SIS: Abingdon, Oxon, U.K., 1998. (23) Nishikiori, N.; Redner, R. A.; Schmidt, Z. An Improved Method for Gas-Lift Allocation Optimization. SPE ATCE, Texas, October 1989; pp105-118, paper no. 19711-MS. (24) Mayo, O. Procedure Optimizes Lift-Gas Allocation. Oil Gas J. 1999, (March), 38–41. (25) Lo, K. K. Optimum Lift-Gas Allocations under Multiple Production Constraints. SPE, 1992; paper no. 26017-MS. (26) Alarcon, G. A.; Torres, C. F.; Gomez, L. E. Global Optimization of Gas Allocation to a Group of Wells in Artificial Lift using Nonlinear Constrained Programming. J. Energy Resour. Technol. 2002, 124 (4), 262–268. (27) Fang, W. Y.; Lo, K. K. A Generalized Well-Management Scheme for Reservoir Simulation. SPE ReserVoir Eng. 1996, 11 (2), 116–120; 29124-PA. (28) Handley-Schachler, S.; McKie, C.; Quintero, N. New Mathematical Techniques for the Optimization of Oil and Gas Production Systems. Presented at the SPE European Petroleum Conference, France, October 2000; paper no. 565161-MS. (29) Buitrago, S.; Rodriguez, E.; Espin, D. Global Optimization Techniques in Gas Allocation for Continuous Flow Gas Lift Systems. Presented at the SPE Gas Technology Symposium, Canada, April 1996; paper no. 35616-MS. (30) Ray, T.; Sarker, R. Multi-objective Evolutionary Approach to the Solution of Gas Lift Optimization Problems. IEEE Congress on Evolutionary Computation, British Columbia, Canada, July 2006; pp 3182-3188. (31) Martinez, E. R.; Moreno, W. J.; Moreno, J. A. Application of Genetic Algorithm on the Distribution of Gas Lift Injection. Presented at the SPE Latin American Petroleum Engineering Conference, Argentina, April 1994; paper no. 26993-MS. (32) Deb, K. Multi-objectiVe Optimization using EVolutionary Algorithms; John Wiley and Sons: New York, 2001. (33) Stoisits, R. F.; MacAllister, D. J.; McCormack, M. D.; Lawal, A. S.; Ogbe, D. O. Production Optimization at the Kuparak River Field Utilizing Neural Networks and Genetic Algorithms. Presented at the SPE Mid-continent Operations Symposium, Oklahoma, March 1999; paper no. 52177-MS. (34) Rashid, K. A Method for Optimal Lift-Gas Allocation and other Production Scenarios, Research Note; Schlumberger AbTC: Abindgon, U.K., December 2006. (35) Rashid, K. A Method for Gas-Lift Optimization. Presented at the SIAM Annual Meeting, Boston, May 2008. (36) Rashid, K. Complex Constraint Handling for Gas-Lift Optimization, Research Note; Schlumberger-Doll Research: Boston, MA, February 2008; OFSR-rn-2008-034-MM-C. (37) Gutierrez, F.; Hallquist, A.; Shippen, M.; Rashid, K. A New Approach to Gas-Lift Optimization using an Integrated Asset Model. Presented at the SPE International Petroleum Technology Conference, U.A.E., December 2007; paper no. 11594-MS. (38) AVocet Gas-Lift Optimization; Schlumberger SIS: Calgary, Canada, 2007. (39) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in C++, 3rd ed.; Cambridge University Press: Cambridge, U.K., 2002. (40) Moitra, S. K.; Chand, S.; Barua, S.; Adenusi, D.; Agarwal, V. A Field-wide Integrated Production Model and Asset Management System for the Mumbai High Field. Presented at the SPE Offshore Technology Conference, Texas, April 2007; paper no. 18678-MS.

ReceiVed for reView May 26, 2009 ReVised manuscript receiVed October 26, 2009 Accepted January 7, 2010 IE900867R