Solving A Gas-Lift Optimization Problem by Dynamic Programming

European Journal of Operational Research 174 (2006) 12201246
www.elsevier.com/locate/ejor
O.R. Applications
Solving a gas-lift optimization problem

by dynamic programming
Eduardo Camponogara *, Paulo H.R. Nakashima
Departamento de Automacao e Sistemas, Universidade Federal de Santa Catarina, Caixa Postal 476, Florianopolis, SC 88040-900, Brazil
Received 17 September 2003; accepted 7 February 2005
Available online 6 June 2005
Abstract
Gas lift is a costly, however indispensable means to recover oil from high-depth reservoirs that entails solving the
gas-lift optimization problem, GOP, often in response to variations in the dynamics of the reservoir and economic oscillations. GOP can be cast as a mixed integer nonlinear programming problem whose integer variables decide which oil
wells should produce, while the continuous variables allocate the gas-compressing capacity to the active ones. This
paper extends the GOP formulation to encompass uncertainties in the oil outow and precedence constraints imposed
on the activation of wells. Recursive solutions are suggested for instances devoid of precedence constraints, as well as
instances arising from precedence constraints induced by forests and general acyclic graphs. For the rst two classes,
pseudo-polynomial algorithms are developed to solve a discretized version of GOP, while the most general version
is shown to be NP-Hard in the strong sense. Numerical experiments provide evidence that the approximate algorithms
obtained by solving the recurrences produce near-optimal solutions. Further, the family of cover inequalities of the
knapsack problem is extended to the gas-lift optimization problem.
2005 Elsevier B.V. All rights reserved.
Keywords: Dynamic programming; Production; Oil industry; Gas-lift optimization; Polyhedra
1. Introduction
To counter the mounting pressure from economic markets and the ever-increasing demand for nonrenewable resources, the oil industry is investing in the development of improved, more ecient technologies
in all of its sectors. In particular, the gas-lift operation of oil elds is one of many production processes
whose performance can be improved. Because the internal pressure in high-depth or depleted reservoirs
*
Corresponding author. Tel.: +55 48 331 7688; fax: +55 48 331 9934.
E-mail addresses: camponog@das.ufsc.br (E. Camponogara), phrn@das.ufsc.br (P.H.R. Nakashima).
0377-2217/$ - see front matter 2005 Elsevier B.V. All rights reserved.
doi:10.1016/j.ejor.2005.03.004
E. Camponogara, P.H.R. Nakashima / European Journal of Operational Research 174 (2006) 12201246
1221
can force the ow of only a fraction of its oil to the surface, the use of articial means becomes imperative
to lift the oil, especially so for deep reservoirs found o-shore. Two examples of articial lifting are submerged pumps and continuous injection of gas. Although the former can in principle recover most of
the oil, its operating costs are excessively high for todays oil prices, not mentioning the potential of an
unfavorable energy trade and other technical hindrances. The gas-lift technique, on the other hand, harnesses the reservoirs gas by injecting natural gas into the production tubing so as to reduce the weight
of the oil column, thereby elevating the mix of oil, gas, and water to the surface. Fig. 1 depicts a small cluster of wells operated by a continuous gas-injection system, consisting of an oil reservoir, a chamber for
compressed gas, and a pool of electric compressors. The issues and problems related to the gas-lift production of an oil eld abound, ranging from the design of the gas-compressing unit and the gas pipeline,
through the modeling of the production ow as a function of the gas injection, to the allocation of compressed gas among wells. Despite recent technological advances, the operating processes can be further improved in a number of ways, especially in the development of more accurate models and the design of
algorithmic solutions to the problems thereof, all with the aim of increasing performance and the degree
of automation.
Herein, the focus is on the allocation of compressed gas to the wells of a cluster, taking into account
issues that play a part in the performance of the process, such as the limits on gas injection, the nonlinear
and uncertain nature of the outow [18], and logic constraints on the activation of the wells. Roughly stated, the problem amounts to rst deciding which wells will produce and second, determining the rate of gasinjection for the active ones. This problem, hereafter referred to as the gas-lift optimization problem, GOP,
is instantiated and solved continuously over time, often after the models of well outow are recalibrated
and as economic factors oscillate. It consists of a mixed-integer nonlinear optimization problem, whose discrete variables indicate which wells are active and which are not, whose continuous variables gauge the rate
of gas injected into the active wells, and whose constraints spell out the physical constraints.
The problem of concern has been the object of investigation in the literature. Grothey and McKinnon [6]
formulated the problem as a mixed-integer nonlinear program (MINLP) and analyzed the performance of a
specially-tailored branch-and-bound algorithm, a Benders decomposition approach whose sub-problems
are MINLPs, and also a Lagrangian relaxation procedure. The results reported from computational experiments revealed that the approximation approaches could nd near-optimal solutions relatively quickly,
while the branch-and-bound algorithm failed at the large instances.
Wang et al. [19] tackled a similar problem, in which they converted the MINLP into a mixed-integer linear programming problem (MILP) by piecewise-linearizing the nonlinear constraints and then applying
meta-heuristics. The same authors [20] later augmented the problem formulation to account for the ow
interactions among wells, but this time the discrete decisions regarding the activation of wells were treated
Fig. 1. A cluster of oil wells.
1222
in an ad hoc manner to allow the application of standard, nonlinear programming algorithms. In these papers, the authors described promising results for instances with a small number of wells, showing that the
latter model induces higher performance in their experimental scenarios.
Martinez et al. [10] applied genetic algorithms to the problem of distributing the ow of compressed gas
to various wells.
An extension of the equal-slope method [9] was proposed by Nishikiori et al. [15], which is based on a
quasi-Newton nonlinear optimization method to nd the optimal injection rates for a group of wells. A
drawback of this method is that it needs a good initial estimation of the injection rates to guarantee its convergence. Besides, convergence cannot be guaranteed if there exists wells that do not respond immediately
to gas-lift. More recently, Alarcon et al. [1] developed a sequential quadratic programming (SQP) algorithm
to replace the quasi-Newton method proposed by Nishikiori et al. [15]. Besides being less prone to failure
due to starting points that are far from the optimal solution, the SQP method has quadratic convergence
near locally optimal solutions and can handle additional constraints. Nevertheless, the SQP method resorts
to an ad hoc rule to handle wells that do not respond instantaneously to the injection of lift-gas or that have
nonlinear jumps, potentially converging to sub-optimal solutions.
This paper extends our preceding work on the application of dynamic programming [2] to approximately
solve the gas-lift optimization problem. Specically, it incorporates families of logic constraints that can be
utilized to express operational constraints on the gas pipelines, which impose further restrictions on the well
activations, and it also admits multiple outow models to hedge against uncertainties in the identication
processes. The algorithms proposed in the literature, in particular those cited above, could in principle be
applied to the problem at hand. A potential impediment to these approaches is the nonlinearities introduced by the multiple outow models, in the case of considering the worst case scenario, which can be naturally handled by dynamic programming. Although the proposed dynamic-programming algorithms only
solve GOP approximately, it compares favorably to the other methods with respect to some criteria. First,
the dynamic-programming approach solves a family of problems, namely one for each amount of gas supply, thereby accounting for the potential failure of one or more compressors. Second, it can be implemented
without the need of sophisticated commercial software to solve MILP and MINLP problems, making it
economically viable to automate the production processes of small- to medium-sized clusters.
The paper is organized as follows. Section 2 formalizes the gas-lift optimization problem, makes a few
assumptions, and describes some of its properties. Section 3 proposes a dynamic-programming algorithm
for the scenario devoid of precedence constraints on the activation of wells; further, it delivers a family of
valid inequalities, which can be extended to facets of the underlying polyhedron by means of lifting, and
reports numerical results from computational experiments. Section 4 focuses on the scenario where the precedence-constraint graph forms a forest, whose nodes are wells and whose arcs model the precedence constraints. Section 5 follows the same pattern of the preceding sections to present a DP algorithm for the case
of general, acyclic precedence-constraint graphs. Section 6 elaborates on extensions to the models and algorithms developed herein. Finally, Section 7 nalizes the paper by drawing some general conclusions and
suggesting directions for future research.
2. Problem formulation
The general gas-lift optimization problem of concern in this paper is of the form:
P G:
J Maximize
Subject to:
N
X
po cno pg cng pw cnw qno qni
n1
qno qni 2
N
X
pi qni
1:0
n1
Qn ;
n 1; . . . ; N ;
1:1
N
X
qni 6 qmax ;
1223
1:2
n1
ln y n 6 qni 6 un y n ; n 1; . . . ; N ;
y i 6 y j ; 8i; j 2 EG;
1:3
1:4
y n 2 f0; 1g;
1:5
n 1; . . . ; N ;
where
N is the number of oil wells;

qni is the rate at which gas is injected into well n;
yn takes on value 1 if the n-th well is activated, and 0 otherwise;
Qn is a set of polynomials that approximate the output ow of well n as a function of the gas injection,
namely Qn fqno;k qni an0;k y n an1;k qni an2;k qni 2 an3;k qni 3 : k 1; . . . ; kng, where k(n) denotes the
number of approximators and the set fani;k g contains the coecients of the polynomials;
qmax is an upper bound on the capacity of the gas-compressing station;
ln and un are lower and upper bounds, respectively, on the amount of gas that can be injected into well n;
po, pg, and pw are respectively the prevailing market prices of oil and gas, and the cost to treat water
before discharge;
pi is the cost incurred to compress gas before injection, typically encompassing the costs incurred by the
maintenance of the compressors and the electric energy to operate them;
cno ; cng , and cnw are respectively the fractions of oil, gas, and water that make up the output ow from well
n, thus cno cng cnw 1; and
G = (V, E) is a directed, acyclic graph whose vertices correspond to the wells, V = {1, . . . , N}, and whose
arcs express constraints on the operation of the wells, specically (i, j) 2 E means that well j must be
active if well i is to be activated.
The elements of Qn are obtained by applying system-identication [8,16,17] processes at well n. In short,
the output ow of a well is measured for a series of gas-injection rates to produce a set of operating points
that can be tted into a polynomial curve, possibly by minimizing the mean-squared error, which then
serves as a model to predict the behavior of the well.
If information about the reliability of these measurements is produced by the underlying identication
processes, in the form of probabilities, then the optimization of P(G) can take place over the expected values of the output ow. This problem can be cast as
P E G:
J E Maximize
fE
N
X
cn qno
n1
Subject to:
N
X
pi qni
2:0
n1
2
qno an0 y n an1 qni an2 qni an3 qni ;

N
X
qni 6 qmax ;
n 1; . . . ; N ;
2:1
2:2
n1
ln y n 6 qni 6 un y n ;
yi 6 yj;
n 1; . . . ; N ;
8i; j 2 EG;
y n 2 f0; 1g;
n 1; . . . ; N ;
2:3
2:4
2:5
where ani Eani;k is the expected value of the i-th coecient of the polynomial qno , and cn
po cno pg cng pw cnw .
1224
In the absence of data on the reliability of the models reported by the identication processes, the optimization can account for the worst case scenario, rendering the following version of the gas-lift optimization problem:
P W G:
J W Maximize
fW
N
X
cn qno
n1
Subject to:
N
X
pi qni
3:0
n1
qno minfqno qni 2 Qn g;

N
X
qni 6 qmax ;
n 1; . . . ; N ;
3:1
3:2
n1
ln y n 6 qni 6 un y n ;
yi 6 yj;
n 1; . . . ; N ;
8i; j 2 EG;
y n 2 f0; 1g;
n 1; . . . ; N .
3:3
3:4
3:5
Proposition 1. PE(G) is NP-Hard.

Proof. The knapsack problem (KP) is reducible to PE(G). Let IKP be an instance of KP consisting of N
items, where cn is the value of the n-th item and wn is its weight, and whose knapsack capacity is b. Generate
an instance IE of PE(G) as follows: set G = (V, ;); set qmax = b; set ln = un = wn for all n 2 V; set an0 cn
and an1 an2 an3 0 for each n; set cno 1 and cng cnw 0; and dene po = 1 while pg = pw = pi = 0. An
optimal solution z fy n ; qni ; qno g to IE yields an optimal solution to IKP, namely the n-th item goes into the
sack if and only if yn = 1. h
Corollary 1. PW(G) is NP-Hard.
Assumption 1. Each qno 2 Qn is a concave function such that qno qni P 0 for qni 2 ln ; un .
Assumption 2. po cno pg cng > pw cnw for each well n, which implies that cn > 0.
The assumptions reect what is typically observed, or can be satisfactorily approximated, in practice: (i)
the outow increases as the gas-injection rises until a saturation point is reached, dwindling as the injection
surpasses saturation owing to the dominating eects of the friction on the production tubing, which is captured by the concavity assumption on qno ; and (ii) the percentage of oil and gas production, co and cg, combined with their market prices should surpass the water treatment cost for the well to be economically
viable, which is the condition imposed by the second assumption.
Corollary 2. fE is a concave function of qi q1i ; . . . ; qNi within the feasible region of PE(G).
Proof. Because cn is a positive constant for each n, fn cn qno pi qni is concave and therefore
fE = f1 + + fn is a concave function of qi. h
Other than polynomials, the outow behavior of an oil well n can be approximated with functions that
n n
n n
have exponential terms, such as the functions ^
qno qni an0 2 ea1 qi an2 ea3 qi [13] and ~qno qni an0 an1 qni
2
an2 qni an3 lnqni 1 [1]. Simple conditions can be imposed on the parameters ani to render ^qno and ~qno con2
2
cave functions: for ^
qno it suces to enforce an0 ; an2 P 0, while 2an2 6 minfan3 =1 ln ; an3 =1 un g guarann
tees concavity of ~
qo . It is worth mentioning that when these functions are a better t to experimental data,
the DP algorithms to be developed hereafter can be easily modied to use ^qno or ~qno .
1225
This paper delivers DP algorithms to solve PE(G) for three classes of precedence-constraint graphs,
namely the class devoid of precedence constraints, PE(;), the class whose precedence constraints form a forest, PE(F), and the broader class comprising general acyclic graphs, PE(G). Owing to the couplings between
oil wells that arise from the precedence constraints, the algorithms become more intricate as E(G) varies
from the empty set, through forests, to acyclic graphs. We favored the study of PE(G) not only because
the algorithmic solution is simpler than its PW(G) counterpart, but also because the DP algorithms for
the former can be readily extended to tackle instances of the latter.
3. Solving P(;)
The only constraint binding the wells in PE(;) is the capacity of the gas-compressing station. If one knew
an optimal allocation of the compressed gas to the oil wells, PE(;) would be decomposed into N concave
optimization problems whose optimal solutions could be easily obtained. The solution procedure proposed
herein is not dissimilar to this two stage procedure, the dierence being that the gas allocation phase is
approximated with standard units of gas and eciently performed by dynamic programming. Let
d = qmax/M be the standard unit of compressed gas, where M is the total number of standard units that
make up the capacity of the compressing station. Assuming that ln = 0 and un = qmax for n = 1, . . . , N,
the number of feasible allocations is precisely the number of ways of placing M indistinguishable balls into
N boxes.
After denoting
the number of feasible allocations by T(M, N), one can deduce that
M N 1
1 N 1
T M; N
. Noticing that T M; N P MN
we conclude that the direct and exhausN 1

N 1
tive enumeration of all possible allocations is intractable for all, but the coarse discretizations of the
gas-compression capacity. To this end, dynamic programming can be applied to nd a solution to PE(;)
if this problem is further subjected to discrete allocations of standard units of gas. For the sake of simplicity, problem PE(;) will be referred to as P within this section. Having introduced the above terminology,
the version of P restricted to the allocation of discrete units of gas-compression capacity can be expressed
as
P M : J M Maximize
fE
N
X
cn qno
n1
Subject to:
N
X
pi qni
4:0
n1
2
qno an0 y n an1 qni an2 qni an3 qni ;

ln y n 6
M
X
qni
6 minfun ; wn dgy n ;
n 1; . . . ; N ;
n 1; . . . ; N ;
wn M;
4:1
4:2
4:3
n1
y n 2 f0; 1g;
n 1; . . . ; N ;
wn 2 f0; . . . ; Mg;
n 1; . . . ; N ;
4:4
4:5
where d = qmax/M is a standard unit of gas, M is the total number of standard units available, and wn is the
number of standard units secured for well n.
The similarities between P(G) and knapsack problems allow us to extend combinatorial properties of the
latter to the domain of the former, in particular families of facet-dening inequalities. In what follows we
extend the knapsack cover inequalities and their lifting procedure to P(;) [14,21]. Let P be the polyhedron
of the feasible solutions to P(;) in which (2.1) has been substituted for the variables qno , i.e.,
P fqi ; y 2 RN BN : subject to constraints families 2.2 and 2.3g, where B f0; 1g. Throughout
the text N f1; . . . ; N g will denote the set of well indices.
1226
Proposition 2. If ln < un 6 qmax for each n, then dimP 2N .

Proof. Let sn ^
qni ; ^y n be the element of P obtained by setting y n ^y n , yk = 0 for all k 5 n, qni ^qni , and
k
qi 0 for all k 5 n. Let S = {sn(ln, 1), sn(un, 1): n = 1, . . . , N}. Clearly, S is a set of 2N linearly independent
elements of P. Since 0 2 P, it follows that P has 2N + 1 afnely independent solutions. h
In what follows, we assume that the conditions of Proposition 2 are satised.
P
P
Denition 1. A set C {1, . . . , N} is a cover if n2C ln > qmax and for each n 2 C, k2Cfng lk < qmax .
Proposition 3. If C is a cover, then
X
y n 6 jCj 1
n2C
is a valid inequality for P. Further, the face F C fqi ; y 2 P :

inequality has dimension dim(FC) = 2jCj 1.
n2C y n
jCj 1g induced by the cover
P
P
P
Proof (Validity). If n2C y n jCj for some solution qi ; y 2 P, then n2C qni P n2C ln > qmax , but this
means that the solution P
is infeasible.
(Dimensionality). Let Pn2C pn y n ln qni 6 po be a general valid inequality for P inducing the maximal
n
face F p fqi ; y 2 P :
n2C pn y n ln qi po and yn = 0 for all n 2 N Cg such that FC Fp. The
developments below prove that FC is maximal by showing that the inequality that induces FC differs from
the one inducing Fp by a positive multiplicative constant.
First, we show that ln = 0 for all n 2 C. For each n 2 C perform the steps: set ^qi ; ^y 0; 0, dene
^qki lk for all k 2 Cn. Clearly ^qi ; ^y 2 F C . Let
Cn = C {j} for some j 2 C, j 5 n, and set ^y k 1 and P
n
~
qi ; ~y ^
qi ; ^y but this time set ~
qi minfun ; qmax k2Cn fng lk g. Because ^qi ; ^y ; ~qi ; ~y 2 F C , the
following equations must hold for ^
qi ; ^y ; ~
qi ; ~y 2 F p :
X
X
pk ^y k
lk ^
qki po ;
k2C n
k2C n
k2C n
pk ~y k
lk ~
qki po .
k2C n
P
Subtracting the rst equation from the second, we deduce that ln ~qni q^ni 0. Since k2Cn lk < qmax , it follows that ~
qni ^
qni > 0 and therefore ln = 0. After applying the above reasoning to each element of C, one
concludes that ln = 0 for each n 2 C.
Second, we show that pn = pk for each pair n, k 2 C. Let Cn = C {k} and dene ^qi ; ^y as follows:
^
qji lj and ^y j 1 for each j 2 Cn. Let Ck = C {n} and dene ~qi ; ~y as follows: ~qji lj and ~y j 1
for each j 2 Ck. Clearly, ^
qi ; ^y ; ~
qi ; ~y 2 F C and therefore the equalities below must hold if
^
qi ; ^y ; ~
qi ; ~y 2 F p :
X
pj ^y j po ;
j2C n
pj ~y j po .
j2C k
Subtracting the latter from

we obtain pn ^y n pk ~y k 0 ) pn pk . It follows from
P the former equality,
P
the solution P
q^i ; ^y that j2Cn pj ^y j po ) j2Cn pj po ) jCj 1pn po . By setting po = jCj 1, we
N
N
n
deduce that
n2C y n 6 jCj 1 induces a maximal face of P \ fqi ; y 2 R B : y n qi 0 for n 2
N Cg. Therefore dim(FC) = 2jCj 1. h
1227
According to Proposition 3, the valid inequality induced by a cover C would be facet-dening if the
cover comprised all of the wells, i.e., C N. One way to strengthen the cover inequality is by augmenting
the cover with the wells whose lower bounds on gas injection are equal to, or higher than, the lower bound
of each well of C. More precisely, the extended cover EC C [ fn 2 N C : ln P lk for all k 2 Cg
induces the extended cover inequality:
X
y n 6 jCj 1
6
n2EC
which is valid for P. The extended

P cover inequalityP(6) is valid because any element i of E(C) C can replace any element j of C, that is, n2Cfjg[fig ln P n2C ln > qmax . The extended cover inequality is said to
dominate the cover inequality because any point violating (5) also violates (6), but a point violating (6) may
not violate (5).
Another way is the principle called lifting [14,21], whereby a valid inequality of
PC P \ fqi ; y 2 RN BN : y n 0; n 2 N Cg can be extended to a valid inequality of P. Starting
with an ordered sequence s = (j1, . . . , jL) of the items from N C, the lifting procedure adds, one at a
time, the variables corresponding to the elements of s toPthe inequality. To add the variables of the rst
well, for instance, the procedure
P produces the inequality n2C y n aj1 y j1 6 jCj 1 that is valid for PC[fj1 g ,
where aj1 jCj 1 maxf n2C y n : qi ; y 2 PC[fj1 g ; y j1 1g. The resulting inequality induces a face
FC[{j1} of PC[fj1 g such that dim(FC[{j1}) P dim(FC) + 1. If the problem solved to compute aj1 attains
j
its maximum at a solution such that qi 1 > lj1 , then it can be shown that dimF C[fj1 g dimF C 2.
The steps outlined above can be systematically applied to the remaining elements of s until the lifted cover
inequality:
X
n2C
yn
L
X
ajk y jk 6 jCj 1
k1
is obtained, which is valid for P and may induce one of its facets under the conditions outlined above.
There is a simple way of computing lower bounds bn for the lifting factors an of the elements n 2 N C.
For k 2 {0, . . . , jCj}, let C(k) C be such that jC(k)j = k and min{ln: n 2 C(k)} P max{ln: n 2 C C(k)}
if k > 0, in other words, C(k) is the P
subset of C with the k elements that require the highest amount of liftgas for activation. Let also qC 0 n2C0 ln for any C 0 C. For n 2 N C, we can dene the lower bound
bn as:
bn Maximize
Subject to:
k 2 f0; . . . ; jCjg
ln P qCk.
Thus, by replacing the lifting factors an in (7) with their lower bounds bn we obtain the strengthened, extended cover inequality:
X
X
yn
bn y n 6 jCj 1
8
n2C
n2ECC
which is also valid for P. Notice that (8) may be stronger than (6) but weaker than (7).
3.1. Recurrences
By limiting the gas-compression capacity to m 2 {0, . . . , M} standard units, we can recast PM in the following recursive form:
1228
Pe m;n:
F m;n Maximize
Subject to:
fn cn qno pi qni
9:0
qno an0 y n an1 qni an2 qni 2 an3 qni 3 ;
9:1
ln y n 6 qni 6 minfun ; mdgy n ;
9:2
y n 2 f0; 1g;
9:3
P m;N : J m;N F m;N ;

P m;n: J m;n Maximize
Subject to:
10:0
F wn ;n J mwn ;n1
wn 2 f0; . . . ; mg.
11:0
11:1
In the above formulation, Jm,n denotes the value of the optimal solution to problem PM restricted to wells
n, n + 1, . . . , N and when m standard units of lift-gas are available. From the recurrences (10.0), (11.0),
(11.1), it follows that PM PM,1 and consequently JM = JM,1. A dynamic programming algorithm can
be readily obtained from the above recurrences. To nd the optimal solution to PM, the recursive procedure
begins by lling in the entries of Jm,N, which entails solving Pe m;N for each m 2 {0, . . . , M} as explained below. Then, the procedure computes the entries of Jm,N1 for each m 2 {0, . . . , M} by calculating the optimal
combined value of F wN 1 ;N 1 J mwN 1 ;N ; wN 1 2 f0; . . . ; mg, whereby problem Pe wN 1 ;N 1 is solved. The procedure continues backwards until it reaches well 1, at which point JM,1 is computed and the optimal solution to PM is thereby obtained.
Problem Pe m;n corresponds to computing the optimal gas injection rate into well n, if any, when the
amount of lift-gas available is md. If md < ln, then the objective value is Fm,n = 0 and well n is not activated.
Otherwise, an optimal injection rate can be eciently computed by considering the points at which fn can
take its maximum. Let:
zm,n be the unconstrained maximum of fn; and
Qm,n = {0, ln, min(un, md)} [ {zm,n: ln 6 zm,n 6 min(un, md)}, which will be shown to be the set of values
qni at which fn may attain its maximum if m standard units of gas are available to well n.
Proposition 4. fn attains its maximum only if qni 2 Qm;n .
Proof. There are two cases: (case 1) if yn = 0, then qni 0 and therefore qni 2 Qm;n ; (case 2) if yn = 1, then
due to the concavity of fn and the convexity of the feasible set, ln 6 qni 6 minfun ; mdg, it follows from the
rst-order optimality conditions that fn attains its maximum either at the extremes of the feasible set or at
the unconstrained maximum zm,n, if the latter is feasible. Hence the claim follows. h
When the available lift-gas equals or exceeds the lower bound, md P ln, an optimal injection rate is obtained by calculating the elements of Qm,n {0}, computing the value of fn qni for each qni 2 Qm;n , and taking its maximum as Fm,n. Whichever amount of lift-gas is available, the optimal rate qni is an element of Qm,n
as stated by Proposition 4.
In what follows, we make some comments about the behavior of Jm,n as m increases and give the steps of
the dynamic programming algorithm.
Remark 1. Jm,n is a nondecreasing function of m.
Remark 2. Let J1 = JM be the value of an optimal solution to PM for a discretization level M = M1. Let
J2 = JM for a discretization level M = M2. If M1 = k M2 for any integer k P 1, then J1 P J2.
1229
Proof. Let P1 PM be the gas-lift optimization problem for the discretization level M = M1. Let P2 PM
be the problem for the discretization level M = M2. The claim follows from the fact that the feasible region
of P1 contains the feasible region of P2. h
DP Algorithm for PM
(1)
(2)
(3)
(4)
(5)
Initialization
For m = 0, . . . , M
Compute Jm,N
Recursion
For n = N1, . . . , 1
For m = 0, . . . , M
Compute J m;n maxfF wn ;n J mwn ;n1 : wn 0; . . . ; mg
The DP algorithm solves PM in H(NM2) steps and demands H(MN) units of memory to store the table
of values Jm,n.
3.2. Computational results
In an eort to assess the quality of the approximation PM to P and the performance of the DP algorithm,
we conducted a number of numerical experiments whose problem instances possess the properties of realworld scenarios [4]. The DP algorithm was implemented in C/C++ programming languages in a Intel 800MHz Pentium-III workstation, equipped with 196 Mb of memory and the GNU Linux operating system.
Table 1 depicts the results for problem instances ranging from 6 to 48 wells and three levels of discretization. The quality of the approximation is measured by the relative distance of the objective induced by the
approximation, JM, to an upper bound computed with continuous relaxation and concave programming,
b
J . By ensuring the concavity of the outow functions qno within the interval [0, un] of gas injection, it was
possible to apply MINOS solver version 5.51 [12] to obtain the upper bounds.
To evaluate the merit of the valid inequalities derived from knapsack covers, we augmented the continuous
relaxations with cuts and computed upper bounds e
J with concave programming. A heuristic was implemented
to iteratively search for strengthened extended covers, as given by (8), that cut fractional solutions. The numerical analysis, as shown in Table 1, indicates that the inequalities failed to reduce the upper bound b
J in some
instances and, when they succeeded, the reductions were small and often occurred in instances of low availability of lift-gas which invariably forced the shuto of many wells. On average, the cover-derived inequalities
induced an upper-bound decrease of about 1% in absolute terms but roughly 25% in relative terms.
All in all, the computational results show that the DP algorithm is relatively fast and, perhaps more importantly, the experiments reveal that the approximated solutions are near-optimal for M suciently large.
4. Solving P(F)
Besides the capacity of the gas-compressing station, the wells in P(F) are coupled by constraints on their
activations which are imposed by the arcs of F. To optimally break P(F) into a set of decoupled concave
problems, one for each well, one would have to know an optimal gas distribution that meets the precedence
constraints. Our solution approach deals with the issue of gas allocation by rst, discretizing the capacity of
the gas-compressing station into M units of equal amount, d = qmax/M, and second, handling the precedence constraints through the topological order induced by F. As in the scenario devoid of precedence constraints, the explicit enumeration of all the feasible allocations can be intractable even in the presence of
1230
Table 1
Quality of approximated solutions to P
Prob. id
# of wells
JM
CPU time (seconds)
Bound b
J
b
J J M = b
J (%)
Bound w/cuts e
J
e
J J M = e
J (%)
10
50
100
920.23
976.66
977.87
0.00
0.02
0.07
978.01
5.91
0.14
0.02
978.01
5.91
0.14
0.02
10
100
300
625.23
625.23
628.36
0.00
0.06
0.44
635.65
1.64
1.64
1.15
634.64
1.48
1.48
0.99
20
989.17
0.10
989.17
0.00
989.17
0.00
100
200
324.21
324.21
0.07
0.26
349.10
7.13
7.13
328.53
1.31
1.31
12
40
70
1000
1926.82
1939.08
1947.08
0.03
0.08
13.66
1949.15
1.15
0.52
0.11
1949.15
0.15
0.52
0.11
12
40
100
1000
1308.63
1311.50
1314.23
0.03
0.16
13.43
1317.47
0.67
0.45
0.25
1317.47
0.67
0.45
0.25
12
50
2079.46
0.03
2079.47
0.00
2079.46
0.00
12
10
500
675.39
691.87
0.01
3.35
710.47
4.94
2.62
697.32
3.14
0.78
24
50
70
280
2150.61
2177.01
2181.89
0.10
0.18
2.33
2224.81
3.34
2.15
1.93
2221.27
3.18
2.03
1.77
10
24
50
100
800
2904.06
2944.07
2946.69
0.10
0.33
18.98
2972.11
2.29
0.95
0.86
2972.11
2.29
0.95
0.86
11
24
70
560
1173.83
1177.51
0.15
8.76
1199.84
2.17
1.86
1195.04
1.77
1.47
12
48
96
400
800
4400.78
4482.05
4495.51
0.65
10.05
39.49
4585.47
4.03
2.26
1.96
4578.36
3.88
2.10
1.82
13
48
100
200
800
5982.86
6072.18
6083.52
0.71
2.63
40.11
6116.58
2.19
0.73
0.54
6116.58
2.19
0.73
0.54
14
48
100
300
700
2380.55
2391.13
2398.32
0.70
5.41
28.73
2465.79
3.46
3.02
2.74
2458.82
3.18
2.75
2.46
Mean
9.56
3.60
2.63
precedence constraintsfor instance, the number of feasible activations of the wells, disregarding the availability of compressed gas, remains exponential in the number of wells for a forest F that forms a star. In
particular, a call to algorithm C#(F) given below computes the number of well activations that meet the
constraints induced by a forest F in linear time.1 When referring to a graph G, hereafter
d+(u) = {v: (u, v) 2 E(G)} denotes the set of vertices directly accessible from u, d(u) = {v: (v, u) 2 E(G)}
1
Algorithm C#(F) can be augmented in a straightforward manner to enumerate all of the feasible activations.
1231
denotes the set of vertices that can reach u directly, d+(u) = jd+(u)j is the out-degree of vertex u,
d(u) = jd(u)j is its in-degree, S(u) = {v: there is a path in G from u to v} is the set of descendants of
u, S+(u) = S(u) {u}, and root(G) = {u: d(u) = 0} is the set of roots of G.
Algorithm C#(F)
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
If jV(F)j = 1 then
return 2
If jroot(F)j = 1 then
return 1 + C#(F[V(F) root(F)])
w
1
For each u 2 root(F) do
w
w*C#(F[S(u)])
return w
Proposition 5. C#(F) correctly computes the number of well activations that meet the precedence constraints
imposed by a directed forest F.
Proof (By induction in the cardinality n of V(F)). For n = 1 the number reported by C# is clearly correct.
For n P 2, two cases have to be considered: (case 1) if jroot(F)j = 1, then C#(F) encompasses the activation of the root of F plus all of the feasible activations of F[V(F)root(F)], when the root is not active,
which is correctly computed by induction hypothesis; (case 2) if jroot(F)j = m P 2, then let
F = T1 [ [ Tm, where Tj, j = 1, . . . , m, are subtrees rooted at the vertices that appear in root(F);
because the wells appearing in a subtree are independent of those
Q appearing in other subtrees, it follows
from induction hypothesis that the number of activations is mj1 C # T j , which is precisely what algorithm C# computes. h
Because C#(F) visits exactly once each vertex of F, we conclude that the algorithms running time belongs
to H(n) where n = jV(F)j.
4.1. Recurrences
Along the lines of the preceding DP algorithm, PE(F) can be approximated by means of a discretization
of the gas-compression capacity. In the developments that follow, P is a shorthand for PE(F) and PM denotes the discretized version of P, where d = qmax/M is the standard unit of gas allocation and M is the
total number of units available. Before stating the recurrences, we provide short descriptions of the subproblems that are recursively solved by the DP algorithm:
P ym;R is the discretized version of P restricted to the wells that are descendants of the elements of R V(F)
within forest F, when the total amount of lift-gas available is m standard units, and such that the trees
rooted at the vertices of R are disjoint; further, each well r 2 R must be activated if y = 1 or else it may
or may not be activated (i.e., yr P y for all r 2 R), and the available compressing-capacity of m standard
units has to be shared among the elements of R and their descendants; J ym;R is the value of the objective
function of an optimal solution to P ym;R .
n
P ym;n
is the discretized version of P restricted to well n and its descendants in F, S+(n), where the activan
tion of well n is determined by yn, and there is m standard units of gas for all of the wells in S(n); J ym;n
is
yn
the value of the objective function of an optimal solution to P m;n .
1232
According to the above terminology, one can conclude that PM is equivalent to P 0M;rootF J M J 0M;rootF
where the latter can be solved by computing the following recursions:
n
P ym;n
:
n
J ym;n
Maximize
Subject to:
P ym;R:
n
fn cn qno pi qni J mw
n ;d n
qno
an0 y n
ln y n 6 qni 6
an1 qni
an2 qni 2
minfun ; wn dgy n ;
wn 2 f0; . . . ; mg;
12:0
an3 qni 3 ;
J ym;R 0 If R ;;
J ^ywn ;n J ymwn ;Rfng : If R 6 ;; pick any n 2 R
J ym;R Maximize
Subject to: wn 2 f0; . . . ; mg;
^y 2 fk 2 B : k P yg.
12:1
12:2
12:3
13:0
13:1
13:2
13:3
It is worth mentioning that the objective-function value of an infeasible problem is negative innity. For
instance, if m units is not sucient to activate all the wells in the set S(n), then P 1m;n is infeasible and
J 1m;n 1.
4.2. DP algorithm
To solve the recursions (12.0)(13.3), the topological order imposed by the edges of F can be used to
n g of sub-problems can be solved rst for the wells at the leaves of F,
advantage. Namely, the set fP ym;n
and then for those at higher levels while respecting the topological order. The algorithm outlined below
behaves much like in the way just suggested. To facilitate the specication of the algorithm, we denote
y
n : m 0; . . . ; Mg by Py n and, similarly, the set fP y
the set of sub-problems fP ym;n
m;R : m 0; . . . ; Mg by PR .
n
y
y
yn
yn
And extending this notation, we have Jn fJ m;n : m 0; . . . ; Mg and JR fJ m;R : m 0; . . . ; Mg.
DP Algorithm for PM(F)
Initialization
Let K = (k1, . . . , kN) be a permutation of {1, . . . , N}
obeying a topological order over F
Recursion
(2)
For k = N, . . . , 1
(3)
n
kk
(4)
If d+(n) = 0 then
(5)
Solve the set of sub-problems Pynn for yn 2 {0, 1} to obtain Jynn
Else-If
(6)
Let d+(n) = {r1, . . . , rd}, where d = d+(n)
(7)
Let J1frd g
J1rd
(8)
For t = d 1 down to 1
(9)
Compute J1frt ;...;rd g J1rt J1frt1 ;...;rd g by solving P1frt ;...;rd g
(10)
Let J0frd g
maxfJ0rd ; J1rd g
(11)
For t = d 1 down to 1
(12)
Compute J0frt ;...;rd g maxfJ0rt ; J1rt g J0frt1 ;...;rd g by solving P0frt ;...;rd g
(13)
Solve Pynn to obtain Jynn for yn 2 {0, 1}
End-If
Termination
(14) By following the steps outlined in the body of the Else-If above,
solve the set of sub-problems PyrootF to obtain JyrootF for y 2 {0, 1}
(1)
1233
Fig. 2. Illustration of the DP algorithm for P(F). The graph on the left denes the precedence constraints on the activation of wells.
For a topological order K = (1, 2, 3, 4, 5, 6), the sets of sub-problems solved and tables generated by the algorithm are depicted in the
tree on the right, which indicates the order in which the problems are solved. For instance, 5 : Pf5;6g ! Jf5;6g means that sub-problem
set Pyf5;6g was the 5th to be solved in order to obtain Jyf5;6g , for y 2 {0, 1}.
At line 1, the topological order K implies that for every i < j, ki 62 S(kj).
n
n
At line 5, the solution of sub-problem set Pynn entails solving P ym;n
and obtaining J ym;n
for m = 0, . . . , M.
1
At line 9, the algorithm solves P m;frt ;...;rd g for each m 2 {0, . . . , M} by: (i) deciding an optimal split of m
units into wrt units to the subtree of wells rooted at rt and the remaining m wrt units to the subtrees whose
roots are the elements of {rt+1, . . . , rd}; and (ii) fetching the solutions to the sub-problems P 1wrt ;rt and
P 1mwrt ;frt 1;...;rd g , namely J 1wrt ;rt and J 1mwrt ;frt 1;...;rd g , which must be available because the sub-problems are
solved in a reverse topological order. The values of J 1m;frt ;...;rd g are therefore computed as
J 1wrt ;rt J 1mwrt ;frt 1;...;rd g .
The way the algorithm solves the recurrences is illustrated in Fig. 2 for a cluster of 6 wells and a simple
precedence-constraint graph.
Two relevant issues are: how many tables does the algorithm compute? how many computational steps
does it take to ll them? If we think of the algorithm as operating on tables that are stored in a persistent
memory, and whose values are lled as the algorithm progresses, then upper bounds on the algorithms running time and its memory usage can be established. Considering
PN the number of tables, the algorithm will
compute 2N tables Jynn (one for each n and yn) and 2 n1 d n tables JyR , two for each n and
Rnt frt ; . . . ; rd g fr1 ; . . . ; rd g d n, t = 1, . . . , d+(n), in addition to 2jroot(F)j tables JyR where
R root(F).P Because each table has H(M) entries, altogether the tables will occupy
H2MN Nn1 d n jrootF j 2 H4MN oat-point units of memory.
Considering the number of computational steps, we can discern the solution of sub-problems between
leaf and nonleaf nodes of F. For a leaf node n, the computational cost to obtain Jynn is of the order
H(M). The computational cost associated with a nonleaf node n, on the other hand, is divided between
the calculation of Jyd n and Jyn . The former is equivalent to solving a problem PM(;) with d+(n) superwells,2 where the objective functions of the superwells are induced by the values of tables already computed, namely Jyk for k 2 d+(n); therefore, the computation of Jyd n terminates in H(d+(n)M2) steps. The
computation of the latter, Jyn , can be performed in H(M2) steps. Adding up all of the computational steps,
one can conclude that the algorithms running time is of the order H(4NM2).
2
A superwell consists of a whole subtree of wells with root at a particular child of vertex n.
1234

We implemented the DP algorithm outlined above on the same computer platform, using the same Unixbased operating system and C/C++ programming language. The algorithm was put to the test in a number
of experiments, ranging from clusters of 636 oil wells and with varying structures of precedence constraints. Table 2 reports gures of these experiments, comprising lower bounds obtained with the DP algorithm, its running time, upper bounds for the problem instances computed from continuous relaxation of
Table 2
Quality of approximated solutions to P(F)
Prob. id
# of wells
JM
CPU time (seconds)
Bound b
J
b
J J M = b
J (%)
Bound w/cuts e
J
e
J J M = e
J (%)
10
80
320
709.30
756.61
757.00
0.00
0.06
0.67
757.43
6.35
0.11
0.06
757.43
6.35
0.11
0.06
10
80
320
650.26
680.17
680.17
0.01
0.05
0.43
686.52
5.28
0.93
0.93
686.52
5.28
0.93
0.93
10
50
300
334.38
350.27
350.33
0.01
0.04
0.38
361.51
7.51
3.11
3.09
359.62
7.02
2.60
2.58
12
20
80
160
1695.18
1786.31
1787.24
0.03
0.21
0.67
1787.71
5.16
0.08
0.03
1787.71
5.16
0.08
0.03
12
20
80
320
1266.86
1287.20
1290.06
0.03
0.15
2.15
1315.96
3.73
2.19
1.97
1300.29
2.57
1.01
0.79
12
50
80
320
686.25
700.35
700.63
0.11
0.34
1.98
743.20
7.66
5.76
5.73
719.36
4.60
2.64
2.60
24
50
200
800
2050.28
2083.66
2101.30
1.85
27.09
425.88
2177.94
5.86
4.33
3.52
2177.94
5.86
4.33
3.52
24
50
200
800
2899.06
2944.07
2946.69
2.00
28.42
444.04
2967.41
2.31
0.79
0.70
2967.41
2.31
0.79
0.70
24
80
200
800
2187.44
2237.62
2244.39
4.37
26.02
609.40
2358.57
7.27
5.15
4.84
2339.96
6.52
4.37
4.08
10
36
100
400
800
3454.53
3543.07
3550.51
45.07
690.10
2746.23
3668.09
5.82
3.41
3.20
3668.09
5.82
3.41
3.20
11
36
100
400
800
3799.32
3884.45
3909.86
45.32
692.44
2758.98
4000.47
5.03
2.90
2.26
4000.47
5.03
2.90
2.26
12
36
200
400
600
3372.76
3386.50
3392.76
169.92
605.94
1354.80
3536.14
4.62
4.23
4.05
3518.45
4.14
3.75
3.57
Mean
296.82
3.61
3.11
1235
the discrete variables with a convex-programming solver Jb and from the continuous relaxation augmented with cuts e
J , and the relative objective distance between lower and upper bounds for a few levels
of discretization.
The cuts are from the families of valid inequalities to be derived in the next section for the polyhedron of
P(G) from knapsack-cover inequalities, which are also valid for the polyhedron of P(F). In the same vein of
the experiments for P, we implemented a heuristic that iteratively searches for strengthened extended covers
that cut the current fractional solution, quitting only when a cut cannot be found after several attempts. If a
cut is found, the corresponding inequality is introduced in the constraint set and the process is repeated.
For some instances cuts could not be found or the upper bound could not be reduced. For the remaining
instances the results depicted in the table show that the reduction on the upper bound was small on average,
about 0.5% in absolute terms but roughly 14% in relative terms.
Considering everything, the numerical experiments indicate that the DP algorithm is eective at approximating the solution to instances of P(F) of the size considered in this research.
5. Solving P(G)
This section addresses the full version of the gas-lift optimization problem, in which the precedence constraints are expressed by an acyclic graph G. The solution approach is not unlike the ones proposed for the
discretized versions of P(;) and P(F): rather than solving P(G), we solve its discretized form PM(G) in which
the gas-compression capacity is divided into M standard units, each with an amount of d = qmax/M. This
allows us to apply dynamic programming to the problem of distributing the gas-compression capacity
among the wells, but without violating the constraints on the activation of the wells. Roughly speaking,
the algorithm works by breaking PM(G) into sub-problems that are solved recursively and whose solutions
are later integrated to obtain a solution to PM(G), whereby each sub-problem is restricted to a reduced set
of wells and whose solution is stored in a data structure for future retrieval. The remainder of the section
presents some results on the structure of the constraint polyhedron obtained by removing the well outputow equations from the constraint set, describes a recursive algorithm for PM(G), and then reports results
from numerical experiments.
5.1. Insights into the structure of the problem
Let PG be the polyhedron of the feasible solutions to P(G) in which (2.1) has been substituted for the
variables qno , i.e., PG fqi ; y 2 RN BN : subject to constraints 2.22.4g. For the developments
P
herein, the graph terminology introduced in the preceding section has to be extended. Let kn k2Sn lk
be the minimum amount of gas necessary to activate well n and, in turn, all of its descendants in G. Also
let W(C) = [n2CS(n) be the set of all wells that have to be activated for the activation of those in
C N f1; . . . ; N g.
Proposition 6. If G is an acyclic graph, ln < un < qmax for all n, and kn < qmax for all n 2 root(G), then PG
is full dimensional, i.e., dimPG 2N .
Proof. Let L = (v1, . . . , vN) be a topological order for G. In what follows we show that unit vectors for qni
and yn can be obtained systematically, from elements of PG, in the order reverse to the one imposed by L.
More formally, the procedure nds unit vectors for the variables associated with well vN, then it nds unit
vectors for well vN1, and so forth until well v1 is reached. The steps of this procedure are given below.
Step 1. Let U = ; be the set of the indices of the wells for whose variables unit vectors have been
obtained. Set k = N.
1236
Step 2. Set n = vk and let ~z ~

qi ; ~y 2 PG be the solution obtained by setting ~y j 1 and q~ji lj for
each j 2 S(n), whereas ~y j 0 and ~
qji 0 for each j 62 S(n). Let ^z ^qi ; ^y 2 PG be the solution
obtained by setting ^y j 1 and ^
qji lj for each j 2 S+(n), while setting ^y n 1 and ^qni
minfun ; qmax kn ln g, and ~y j 0 and ~
qji 0 for each j 62 S(n). Notice that the solutions ~z and ^z agree
in all of their entries, with the exception of the gas-injection at well n, more precisely ~qni < ^qni because
kn < qmax. By induction hypothesis, we can assume that unit vectors have been obtained for all of the
nonzero entries of ~z other than those corresponding to well n, since those wells are preceded by well n
in the topological order. By linear combination of these unit vectors and the solutions ~z and ^z, we
can obtain unit vectors for the variables qni and yn. Therefore, U can be extended to encompass well
n, U
U [ {n}. If k = 1 then stop, or else decrement k by 1 and repeat from the beginning of Step 2.
At this point U = {1, . . . , N}, meaning that unit vectors have been obtained from elements of PG for all
of the variables. Since the origin belongs to PG, we conclude that PG contains 2N + 1 afnely independent elements, implying that dimPG 2N . h
From this point on, we assume that the conditions set down by Proposition 6 are met.
Denition 2. A set C {1, . . . , N} is a cover for PG if:
P
S(m) \ S(n) = ; for all m, n 2 C, m 5 n;
Pn2C kn > qmax ; and
n2C kn lk < qmax for each k 2 C.
Proposition 7. If C is a cover, then the cover inequality
X
y n 6 jCj 1
14
n2C
P
n
is valid for PG and, further, the face F C fqi ; y 2 PG :
n2C y n jCj 1, y n qi 0 for all
n 62 W(C)} induced by the cover inequality has dimension dim(FC) = 2jW(C)j 1.
P
P
P
Proof (Validity). If n2C y n jCj for some solution qi ; y 2 PG, then n2W C qni P n2C kn > qmax , but
this means that the solution is infeasible.
(Dimensionality) Let:
X
pn y n ln qni 6 po :
15
n2W C
general valid inequality for PG inducing the maximal face F p fqi ; y 2 PG :

ln qni po and y n qni 0 for all n 2 N W Cg such that FC Fp. Below we outline
a proof that FC is maximal by showing that inequality (14) (which induces FC) differs from inequality
(15) (which induces Fp) by a positive multiplicative constant. The demonstration consist of showing that
pn = 1 for all n 2 C, pn = 0 for all n 2 W(C) C, ln = 0 for all n 2 W(C), and p0 = jCj 1.
First, we show that ln = 0 for n 2 C. Let ~z ~y ; ~qi be obtained by setting ~y t 1 and ~qti lt for all
t
tP2 W(C {k}),
P where k 2 C {n}, and ~y t 0 and ~qi 0 for all t 2 S(k). Clearly, ~z 2 F C and
t
~
k
<
q
.
Similarly,
let
^
z
be
identical
to ~z with the exception that
max
t2W C qi
t2Cfkg t
P
n
t
^
qi g. For ~z and ^z to belong to Fp, (15) must hold at equality. Subtracting
qi minfun ; qmax t2W Cfng ~
(15) evaluated at ~z from (15) evaluated at ^z, and taking into account the fact that ^qni > ~qni , one concludes
that ln = 0. Repeating the above steps for the remaining elements of the cover, leads us to verify that
ln = 0 for all n 2 C.
be
P a
n2W C pn y n
1237
Second, we show that ln = 0 for n 2 W(C) C, where n 2 S+(l) for some l 2 C. Let ~z ~y ; q~i be
t
obtained by setting ~y t 1 and ~
qti lt for all t 2PW(C {k}),
Pwhere k 2 C {l}, and ~y t 0 and ~qi 0 for
t
all t 2 S(k). It follows that ~z belongs to FC and t2W
P C ~qi t2Cfkg kt < qmax . Let ^z 2 F C be identical to ~z,
except for ^
qni which takes on value minfun ; qmax t2W Cfng ~qti g. Notice that (15) must hold at equality for
~z and ^z because ~z; ^z 2 F p . By subtracting (15) evaluated at ~z from (15) evaluated at ^z, and using the fact that
^
qni > ~
qni , one can deduce that ln = 0. Repeating the above steps for the remaining wells, leads us to conclude
that lP
n = 0 for all n 2 W(C) C. At this point, we have discovered that the inequality inducing Fp is of the
form n2W C pn y n 6 p0 .
Third, we show that pn = 0 for all n 2 W(C) C. This, however, is shown by induction and exploiting
the structure imposed by the precedence constraints. Let L = (v1, . . . , vjW(C)Cj) be a permutation of the
elements of W(C) C such that vj precedes vk in a topological order of G if and only if j > k, in other
words, the elements of W(C) C are arranged in L in the order reverse to a topological order. It can be
shown by induction in l 2 {1, . . . , jW(C) Cj} that pv1 0 (induction basis), pv2 0 (induction step), and so
forth until the claim is demonstrated. We omit the induction basis and show only the induction step. For
some l P 2, let n = vl where n 2 S(k) for some k 2 C. Let ~z ~y ; ~qi be obtained by setting ~y t 1 and
t
~
qti lt for all t 2 W(C {k}) [ S(n), while ~y t 0 and
Because ~y k 0 and
P ~qi 0 for all
Pt 2 S(k) S(n). P
~y t 1 for all t 2 S(n) it follows that ~z 2 F C . Further, t2W C pt ~y t t2W CSk pt ~y t t2Sn pt ~y t p0 for
~z to belong to Fp. By induction, it follows that pt = 0 for all t 2 S(n) {n}. Thus
X
pt ~y t pn ~y n p0 .
16
t2W CSk
Now, let ^z 2 F C be equal to ~z, but with the exception that ^y t q^ti 0 for all t 2 S(n). Eq. (16) evaluated at ^z
must hold for ^z to belong to Fp. By subtracting (16) evaluated at ^z from (16) evaluated at ~z, it follows that
pn = 0. Repeating the above steps for l + 1 until l = jW(C) Cj, we can deduce that pn = 0 for all
n 2 W(C) C, leading us to conclude that the inequality (15) inducing Fp is of the form:
X
pn y n 6 p0 .
17
n2C
Fourth, it is easy to verify that pn = pk = p for all n, k 2 C, n 5 k. Just produce two elements of FC, one
activating well n but not k, say ~z, and another activating k but not n, say ^z. To have ~z; ^z 2 F p , inequality (17)
must hold at equality if evaluated for ~z and ^z. By subtracting (17) evaluated at ~z from (17) evaluated at ^z, we
conclude that pn = pk. Consequently, p0 = (jCj 1)p.
By setting p0 = jCj 1 and pn = 1 for all n 2 C, we deduce that inequality (14) induces a maximal
face of PG \ fqi ; y 2 RN BN : y n qni 0 for all n 2 N W Cg. Therefore dim(FC) = 2jW
(C)j 1. h
According to Proposition 7, the inequality induced by a cover C would be a facet of PG if W C N.
By applying the principle of lifting, one can extend a cover inequality for PU G
PG \ fqi ; y 2 RN BN : y n 0 for all n 2 N U g, where U = W(C), to a valid inequality of PG.
As in the lifting procedure outlined for cover inequalities of P (see Section 3), the variables associated with
the elements of N U would be lifted one at a time. Let s = (j1, . . . , jL) be an ordered sequence of
Pthe elements of N U . To lift the rst element of this sequence, j1, and produce the inequality
n2C y n
aj 1 y j 1 6 P
jCj 1 that is valid for PU [fj1 g G, the procedure would compute aj1 jCj 1
maxf n2C y n : qi ; y 2 PU [fj1 g G; y j1 1g by solving a problem related to P(G). The resulting inequality induces a face F U [fj1 g of PU [fj1 g G such that dimF U [fj1 g P dimF U 1. If there exists an optimal
j
solution to the problem of calculating aj1 such that qi 1 > lj1 , then it should not be dicult to show that
the induced face has dimension dimF U [fj1 g dimF U 2. Repeating this procedure for the remaining
elements of s may produce a stronger inequality, possibly a facet-inducing inequality for PG referred
to as lifted cover inequality:
1238
yn
n2C
L
X
ajk y jk 6 jCj 1.
18
k1
The inequality (18) depends as much on the cover C as on the order in which the elements of N W C are
lifted. Because the computation of the lifting factors is hard, one can attempt to obtain lower bounds bn for
the factors an. In what follows we describe a simple way of computing lower bounds. A subset EC N is
an extended cover of C if:
1. C E(C);
2. S(m) \ S(n) = ; for all m, n 2 E(C), m 5 n; and
3. kn P max{kj : j 2 C} for each n 2 E(C) C.
Then, it can be shown that the extended cover inequality:
X
y n 6 jCj 1
19
n2EC
is valid for PG. Notice that in (19) bn = 1 for all n 2 E(C) C while bn = 0 for all n 2 N EC. As in
the preceding developments, we can obtain higher lower bounds for the lifting factors. For k 2 {0, . . . , jCj},
let C(k) C be such that jC(k)j = k and min{kn: n 2 C(k)} P max{kn: n 2 C C(k)} if k > 0, that is, C(k)
is the subset of C P
encompassing the wells that need the highest minimum amount of lift-gas to be activated.
Let also qC 0 n2C0 kn for any C 0 C. For n 2 E(C) C, we can dene a lower bound bn as
bn Maximum
Subject to:
k 2 f1; . . . ; jCjg
kn P qCk.
By replacing the lifting factors an in (18) with bn for n 2 E(C) C we obtain the strengthened, extended
cover inequality:
X
X
yn
bn y n 6 jCj 1
20
n2C
n2ECC
which is valid for PG.

5.2. Recurrences
Before stating the recurrences, we provide a short description of the problems to be solved recursively by
the algorithm, and the tables produced thereof. The recurrences are centered around solving two subproblems:
Pm(H) is the version of the original problem restricted to the wells appearing in V(H), the precedence
constraints imposed by E(H), and a total amount of m standard units of lift-gas, where H is a subgraph
of G; Jm(H) denotes the objective value of the optimal solution to Pm(H), whereas JH
fJ m H : m 0; . . . ; Mg is the table with the objective value for all levels of gas-compression capacity;
PH denotes the set of sub-problems {Pm(H): m = 0, . . . , M}.
P rm H is a more stringent version of Pm(H) in which well r 2 root(H) should be active; J rm H denotes
the objective value of the optimal solution, which can take on value 1 if the amount md is not sucient to activate r and all of its descendants; Jr H fJ rm H g denotes the table with the objective values for all levels of gas-compression capacity, while Pr H fP rm H g denotes the set of corresponding
sub-problems.
Having introduced the above notation, we are in a position to state the recurrences.
X
P m H : J m H Maximize
fn qni if EH ;
1239
21:0
n2V H
Subject to: wn 2 f0; . . . ; mg;

X
wn m;
n 2 V H ;
21:1
21:2
n2V H
ln y n 6 qni 6 minfun ; dwn gy n ;

J m H
y n 2 f0; 1g;
r
maxfJ m H ; J m H V H
P rm H : J rm H Maximize
n 2 V H ;
21:3
n 2 V H ;
frgg for any r 2 rootH ;
if EH 6 ;;
fn qni if S H r V H
21:4
21:5
22:0
n2V H
Subject to:
wn 2 f0; . . . ; mg;
X
wn m;
n 2 V H ;
22:1
22:2
n2V H
ln 6 qni 6 minfun ; dwn g;

J rm H
Maximize
Subject to:
J rk H S H r
n 2 V H ;
J mk H V H S H r
22:3
if S H r 6 V H
k 2 f0; . . . ; mg.
22:4
22:5
In the above recurrences, SH(r) denotes the set of descendants of node r within graph H. The algorithms
that solve the recurrences use two global data structures, namely a dictionary D [3] to store tables JH
and Jr H , and a topological order T = (t1, . . . , tN) of G which guides the recursive solution of P(G).
The pseudo-codes of the algorithms SolveJ(H) and SolveJ(H, r), which recursively compute tables JH
and Jr H respectively, are described below. To solve P(G), one has to set up the global data structures
and invoke a call to SolveJ(G). At termination, a solution to P(G) can be obtained from table JG and
the remaining ones stored in D.
Algorithm SolveJ(H)
If JH 2 D then
Return JH
End-If
(3) If E(H) = ; then
(4)
Solve Pm(H) for m = 0, . . . , M to obtain JH
(5)
Store JH in D
(6)
Return JH
Else
(7)
Let r
tk, where k is the smallest index such that tk 2 V(H)
(Notice that r 2 root(H))
SolveJ H ; r
(8)
Jr H
(9)
JH r
SolveJ H r where Hr = H[V(H) {r}]
(10)
Compute JH by setting J m H maxfJ rm H ; J m H r g for m = 0, . . . , M
(11)
Store JH in D
(12)
Return JH
End-If
(1)
(2)
1240
Algorithm SolveJ(H, r)
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
(10)
(11)
(12)
(13)
If Jr H 2 D then
Return Jr H
End-If
If V(H) = SH(r) then
Solve P rm H for m = 0, . . . , M to obtain Jr H
Store Jr H in D
Return Jr H
Else
Let H1 = H[SH(r)]
Jr H 1
SolveJ H 1 ; r
Let H2 = H[V(H) SH(r)]
JH 2
SolveJ H 2
Compute Jr H by setting J rm H maxfJ rk H 1 J mk H 2 : k 0; . . . ; mg
for m = 0, . . . , M
Store Jr H in D
Return Jr H
End-If
The way the above algorithms operate is illustrated in Figure 3 for a simple instance consisting of 6 wells.
The graph with precedence constraints on the activations of wells is shown at the root of the sub-problem tree.
The algorithm follows the topological order T = (3, 1, 2, 4, 5) to decompose problems into sub-problems.
From the augmentations introduced by bookkeeping tables already computed3 and solving a problem
equivalent to P(;) whenever E(H) = ;, one can conclude that the algorithms running time is bounded
by O(jDjNM2) steps, where jDj is the number of tables stored in D. Despite these augmentations, a call
to SolveJ(G) can take exponential time in certain problem instances. Let G be a precedence constraint graph
where V = {1, . . . , 2k + 1} and E = {(1, k + 2), (2, k + 3), . . . , (k, 2k + 1)} [ {(k + 1, k + 2), . . . , (k + 1,
2k + 1)}. For the topological order T = (1, 2, . . . , 2k + 1), SolveJ(G) will compute X(2k) distinct tables
and, consequently, jDj 2 X(2k). On the other hand, the running time would be linear in the number of wells
if the topological order had well k + 1 at the head. This illustrates that the algorithms performance will
depend as much on G as on its topological order.
In light of the above discussion, one would wonder whether there exists a set of recurrences for PM(G)
that could be solved in pseudo-polynomial time, or else if this problem is NP-Hard in the strong sense [5]. A
problem is NP-Complete (NP-Hard) in the strong sense if the existence of a pseudo-polynomial algorithm
would deliver a polynomial time algorithm for the class of NP-Complete problems. The question raised
above is answered by the following proposition, which shows that PM(G) is NP-Hard in the strong sense
by means of a reduction from the partially ordered knapsack problem [5,7].
Proposition 8. PM(G) is NP-Hard in the strong sense.
Proof. An instance IK of the partially ordered knapsack problem is given by a nite set U, a partial order
0
on U, for each u 2 U a size s(u) 2 Z+ and a value v(u) 2 ZP
+, and positive integers
P b and k. Is there U U
0
0
0
0
such that if u 2 U and u u, then u 2 U , and such that u2U 0 su 6 b and u2U 0 vu P k? The algorithmic answer to this question constitutes the partially ordered knapsack problem. An instance IK can be
3
Without the data structure D, SolveJ(G) would compute H(2N) tables for problem instances that could force multiple solutions of
the same recurrences, such as the instances in which E(G) = ;.
1241
Fig. 3. An illustration of the recursive algorithms that solve P(G). The precedence-constraint graph G has 6 nodes as illustrated at the
root of the sub-problem tree. The topological order of G used to guide the solution is T = (3, 1, 2, 4, 5). Each node indicates a set of
sub-problems solved by the algorithm: P1; 2; 4 means the sub-problem set PGf1; 2; 4g.
reduced to an instance IP of PM(G) as follows. Assume w.l.o.g. that U = {1, . . . , N}. Let M = b, qmax = b,
and for each n 2 U set ln = un = s(u) and an0 vu while an1 an2 an3 0. Dene cno 1 while cng cnw 0
for each n, and set po = 1 while pg = pw = pi = 0. Further, let G = (V, E) where V = U and (u,u 0 ) 2 E if
and only if u 0 u. Then a solution SP to IP has objective value greater than or equal to k if and only if
there exists U 0 U with the desired properties, in other words, u 2 U 0 if and only if yu = 1. h
In spite of the unfavorable worst-case performance of the recursive algorithm, it produced solutions to
representative instances within time limits comparable to the running time of the algorithm for P(F). The
experimental evidence supporting this claim will be reported in a moment.
As a nal note, we would like to mention that the number of well activations that meet the precedence constraints imposed by the directed, acyclic graph G can be computed by an algorithm C#(G) that has the same
structure of the one proposed to solve P(G). Because this algorithm C#(G) would in essence have the same
recurrences, but without having to actually ll the tables, its running time would be bounded by O(jDj).
Noticing that the number of feasible well activations of a forest F is not smaller than the number of well activations of an acyclic graph G if E(F) E(G), one can alternatively use the well-activation counting algorithm
of the preceding section to compute an upper bound for the number of feasible activations for G.
The algorithm proposed to solve the recurrences (21.0)(21.5), (22.0)(22.5) was implemented in
our computer platform, using C/C++ programming language and the GNU Linux operating system. To
1242
simplify the programming eort, we used the data structures provided by LEDA library [11], mostly lists
and dictionaries. The quality of the approximation PM(G) to P(G) as well as the performance of the proposed algorithm were assessed in a number of experiments, with varying number of wells and structure of
the precedence constraints imposed by the edges of G. The results of these experiments appear in Table 3,
including parameters of the problem instances, lower bounds produced by the approximated solutions,
upper bounds, and the computational time.
Table 3
Quality of approximated solutions to P(G)
Prob. id
# of wells
JM
CPU time (seconds)
Bound b
J
b
J J M = b
J (%)
bound w/cuts e
J
e
J J M = e
J (%)
10
80
320
709.30
756.68
757.00
0.01
0.05
0.65
757.43
6.36
0.10
0.06
757.43
6.36
0.10
0.06
20
80
320
628.29
633.35
633.41
0.03
0.05
0.53
633.65
0.85
0.05
0.04
633.65
0.85
0.05
0.04
50
200
500
324.93
324.95
324.96
0.03
0.24
1.48
339.56
4.31
4.30
4.30
336.86
3.54
3.54
3.54
12
40
80
320
1774.67
1786.31
1787.58
0.08
0.29
3.94
1787.71
0.73
0.08
0.07
1787.71
0.73
0.08
0.07
12
40
80
320
1097.57
1112.83
1123.74
0.60
0.21
3.01
1144.33
4.09
2.75
1.80
1144.33
4.09
2.75
1.80
12
40
80
320
677.19
677.24
677.25
0.03
0.17
1.72
703.35
3.72
3.71
3.71
701.71
3.49
3.48
3.48
24
80
320
640
3646.45
3699.74
3703.53
1.47
21.25
84.17
3707.21
1.64
0.20
0.10
3707.21
1.64
0.20
0.10
24
80
320
640
1810.07
1820.37
1821.95
1.37
18.22
73.29
1942.58
6.82
6.29
6.21
1931.45
6.28
5.75
5.67
24
80
320
640
1226.30
1226.84
1226.90
1.20
15.98
63.08
1340.75
8.54
8.50
8.49
1336.69
8.26
8.22
8.21
10
36
80
200
800
1110.66
1117.85
1119.08
3.74
17.31
280.03
1182.99
6.11
5.51
5.40
1174.92
5.47
4.86
4.75
11
36
80
200
800
674.27
674.33
674.32
3.05
14.17
206.21
700.23
3.70
3.70
3.70
693.24
2.74
2.73
2.73
12
36
80
200
800
4078.28
4204.99
4263.75
4.42
22.17
556.81
4312.17
5.43
2.49
1.12
4312.17
5.43
2.49
1.12
Mean
36.14
3.47
3.19
1243
Table 4
Running time and memory usage
Edges (m)
Running time (seconds)

Nodes (n)
Number of sub-problems
Sub-problems (jDj)
24
36
48
56
Mean
24
36
48
56
Mean
dn/8e
dn/4e
dn/3e
dn/2e
d2n/3e
n
2n
3n
4n
5n
dn2/3e
dn2/2e
12.3
20.9
16.7
30.2
23.0
27.9
47.5
45.4
34.9
32.4
24.7
22.8
17.2
97.3
153.8
164.7
116.9
130.4
171.8
188.3
133.9
107.4
61.2
57.2
23.8
475.9
473.3
946.7
1297.2
1220.9
504.0
557.4
451.1
264.5
114.5
101.7
135.2
1783.9
2836.6
4948.6
2295.1
1890.1
1379.0
1105.3
729.6
451.1
141.6
134.5
47.1
594.5
870.1
1522.6
933.1
817.3
525.6
474.1
337.4
213.9
85.5
79.1
36
110
141
213
323
375
232
145
92
77
46
35
144
645
1275
1657
2348
3455
1211
671
358
258
77
58
240
3360
4702
7119
13498
15055
16291
2654
1189
555
119
77
696
6786
19455
31765
43549
57633
13874
4998
1928
1142
124
88
279
2725
6393
10189
14930
19130
7902
2117
892
508
92
65
Mean
28.3
116.7
535.1
1485.9
541.7
152
1013
5405
15170
5435
Upper bounds b
J were calculated from the continuous relaxation of P(G), which falls in the class of concave programming, and using MINOS solver [12].4 Aiming to reduce the upper bounds b
J , we implemented
a heuristic that randomly and systematically searches for a cover C and extended cover E(C) that induce a
strengthened, extended cover inequality (20) cutting o the current fractional solution. The heuristic iterates
until it nds a cutting plane or it fails to nd one for several iterations. The upper bound obtained from the
concave program augmented with cutting planes is denoted by e
J . The results reported in Table 3 show that
the upper bound reduction yielded by the strengthened, extended cover inequalities was small: on average,
the reduction was about 0.30% in absolute terms and 8% in relative terms.
In an eort to assess the experimental running time of the recursive algorithms, we solved 10 randomly
generated instances of problems of varying number of wells. The average running time and number of subproblems solved in these experiments appear in Table 4. The discretization level adopted for these experiments was M = 400. The experiments show that the size of the sub-problem tree is small for precedence
constraint graphs with small and high numbers of edges, but it tends to increase if the edge set has intermediate size. In the worst scenarios (56 wells and jEj = dn/2e), the algorithm took 1.4 hours on average to
generate the sub-problem tree, which is still acceptable for day to day operations. The use of virtual memory was necessary to store the data structure D, forcing the frequent retrieval of data from the hard disk
and thereby slowing down the computations. One way of expediting the solution is to increase the random
access memory. Another way is to use a dual processor workstation and solve sub-problems in parallel.
All in all, the experimental evidence indicate that the proposed algorithm can produce near-optimal solutions relatively quickly for small- to medium-sized instances.
6. Extensions
The developments reported heretofore can be extended in a number of ways and without substantial effort to tackle more general problems. Below, we sketch two extensions that include the multimodel outow
version of P(G) and problem instances whose precedence constraint graph G contains cycles.
4
The polynomial describing the output ow of each well n as a function of the gas injection, qno qni , induces a concave function in the
interval [0, un], this way ensuring that the continuous relaxation is a concave programming problem.
1244
6.1. Solving PW
Thus far, this paper has focused on versions of the gas-lift optimization problem in which the outow of
each well n is modeled by a single function, qno qni , and the maximization of the expected performance of a
cluster of wells, when estimates on the likelihood of each element of Qn are known. To hedge against the
worst-case scenario, the engineer may prefer to maximize the overall performance of an oil eld subject to
the minimum outow possible for each well, that is to say, qno minfqno qni 2 Qn g. Despite the additional
complications, it is relatively easy to adapt the algorithms designed to solve PM(;), PM(F), and PM(G) to
tackle their worst-case versions. For problem PM(;), for instance, one only need to change the last line of
the algorithm, replacing F wn ;n by the minimum over the outow functions qno 2 Qn . This change incurs a
penalty in the algorithms running time, which will be of the order H(LNM2) where L =
max{jQnj: n = 1, . . . , N}. Similarly, the algorithm proposed to solve PM(F) can be extended to solve its
n : just let f take on the value corresponding to the minimum among
worst-case version by augmenting P ym;n
n
the outow functions appearing in Qn. The running time of the resulting algorithm will belong to
H(LNM2). In summary, the worst-case algorithms will be L times slower than their single outow and expected performance counterparts.
6.2. Handling cyclic graphs of precedence constraints
Consider the discretized version PM(G) of a gas-lift optimization problem whose well-activation precedence graph G contains cycles. In what follows, we outline a procedure to solve such an instance of P(G)
which invokes the algorithm designed to solve instances whose precedence constraint graph is acyclic. First,
compute the strongly connected components C1, . . . , CK of G and model those containing more than one oil
well as superwells. Because each strongly connected component Ck will have all of its wells active or inactive, its corresponding superwell will behave like a single well. Second, solve a family of problems similar to
P(;) for each superwell, where all of its wells are active and for each m 2 {0, . . . , M}, thereby producing a
table with the solutions which is subsequently stored in the dictionary D. Third, solve a problem PM(G 0 )
where the vertices of G 0 correspond to the elements of {C1, . . . , CK} and whose edge set induces an acyclic
graph. Besides the computation of the connected components and the pre-calculation of tables, one has to
treat superwells dierently while applying the algorithm presented in Section 5 to solve PM(G 0 ): namely
look up a value in a table rather than solving a problem associated with a superwell. From the above outline of the procedure, it follows that the procedures running time is not worse but possibly faster than the
algorithm that tackles acyclic precedence constraint graphs.
7. Closing remarks
Nearly all sectors of the world economy rely on the oil industry. As the demand for fossil fuels continues
to rise, the industry searches for means to better run its production, processing, and transportation systems.
This paper contributed to the state-of-the-art of the oil industrys production systems by more formally
addressing the gas-lift optimization problem, GOP, the problem of deciding upon the activation of oil wells
and the amount of gas injection that renders a near-optimal overall performance. Besides introducing more
rigor to this problem, whose literature is somewhat sketchy, the paper shed light into the nature of the
problem and delivered simple, yet eective algorithms for small- to medium-sized instances. The algorithms
tackled discretized forms of GOP, a nonlinear mixed integer optimization problem, for varying complexity
of precedence constraints on the activation of wells. Despite the approximate nature of the solution approach, the computational evidence indicates that it can yield near-optimal solutions. While pseudo-poly-
1245
nomial algorithms were designed for the simplest versions of the problem, the most complex one was shown
to be NP-Hard in the strong sense.
Our future research agenda includes the integration of the optimization software with process control
systems and the application of this technology to real-world instances. We also intend to look into the
use of other algorithms, such as branch-and-cut with piecewise-linearization, and extend the models and
algorithms to handle more general versions of GOP. This research eort will be conducted in collaboration
with the Brazilian Oil Company.
Acknowledgements
This research was funded in part by Agencia Nacional de Petroleo, Brazil, under grant PRH 34/aciPG/
ANP, and by a grant from PRPG/FunPesquisa, Universidade Federal de Santa Catarina. We extend our
thanks to Dash Optimization for providing a free license of the Xpress-MP optimization software.
References
[1] G.A. Alarcon, C.F. Torres, L.E. Gomez, Global optimization of gas allocation to a group of wells in articial lift using nonlinear
constrained programming, ASME Journal of Energy Resources Technology 124 (4) (2002) 262268.
[2] E. Camponogara, P.H.R. Nakashima, Applying dynamic programming to a gas-lift optimization problem, in: Proceedings of the
2nd Brazilian Conference on Research and Development in Petroleum and Gas, Rio de Janeiro, Brazil, 2003.
[3] T.H. Cormen, C.E. Leiserson, R.L. Rivest, Introduction to Algorithms, MIT Press, Cambridge, MA, 1990.
[4] K. Dutta-Roy, J. Kattapuram, A new approach to gas-lift allocation optimization, in: Proceedings of SPE Western Regional
Meeting, Long Beach, California, June 1997, SPE Paper 38333, pp. 685691.
[5] M.R. Garey, D.S. Johnson, Computers and Intractability: A Guide to the Theory of NP-Completeness, W.H. Freeman and
Company, New York, NY, 1979.
[6] A. Grothey, K. McKinnon, Decomposing the optimization of a gas lifted oil well network, Technical Report MS-00-005,
Department of Mathematics and Statistics, University of Edinburgh, Edinburgh, Scotland, March 2000.
[7] O.H. Ibarra, C.E. Kim, Scheduling for maximum prot, Technical Report 75-2, University of Minnesota, Minneapolis, MN,
1975.
[8] J.-N. Juang, Applied System Identication, Prentice Hall, 1993.
[9] E.P. Kanu, J. Mach, K.E. Brown, Economic approach to oil production and gas allocation in continuous gas lift, Journal of
Petroleum Technology 33 (1981) 11871892.
[10] E.R. Martnez, W.J. Moreno, J.A. Moreno, R. Maggiolo, Application of genetic algorithm on the distribution of gas lift injection,
in: Proceedings of the 3rd SPE Latin American/Caribbean Petroleum Engineering Conference, Buenos Aires, Argentina, April
1994.
[11] K. Mehlhorn, S. Naher, LEDA: A Platform for Combinatorial and Geometric Computing, Cambridge University Press, 1999.
[12] B.H. Murtagh, M.A. Saunders, MINOS 5.5 Users Guide, Systems and Optimization Laboratory, Department of Operations
Research, Stanford University, 1998.
[13] P.H.R. Nakashima, E. Camponogara, Otimizacao da alocacao de gas de injecao para um conjunto de pocos de petroleo operando
via gas-lift, in: Anais do Congresso Brasileiro de Automatica, 2004 (in Portuguese).
[14] G.L. Nemhauser, L.A. Wolsey, Integer and Combinatorial Optimization, John Wiley & Sons, 1988.
[15] N. Nishikiori, R.A. Redner, D.R. Doty, Z. Schmidt, An improved method for gas lift allocation optimization, ASME Journal of
Energy Resources Technology 117 (1995) 8792.
[16] A. Plucenio, Stabilization and optimization of an oil well network operating with continuous gas lift, in: Proceedings of the SPE
International Student Paper Contest at the SPE Annual Technical Conference and Exhibition, San Antonio, TX, USA, 2002.
[17] A. Plucenio, Automation of oil wells production operating by continuous gas lift, Masters thesis, Federal University of Santa
Catarina, Florianopolis, Brazil, 2003 (in Portuguese).
[18] R.F. Stoisits, E.C. Batesole, J.H. Champion, D.H. Park, Application of nonlinear modelling of rigorous representation of
production facilities in reservoir simulation, in: Proceedings of the 67th Annual Technical Conference and Exhibition of the SPE,
Washington, DC, USA, 1992.
1246
[19] P. Wang, M. Litvak, K. Aziz. Optimization of production from mature elds, in: Proceedings of the 17th World Petroleum
Congress, Rio de Janeiro, Brazil, 2002.
[20] P. Wang, M. Litvak, K. Aziz, Optimization of production operations in petroleum elds, in: Proceedings of the SPE Annual
Technical Conference, Richardson, USA, 2002.
[21] L.A. Wolsey, Integer Programming, John Wiley & Sons, 1998.

Solving A Gas-Lift Optimization Problem by Dynamic Programming

Uploaded by

Copyright:

Available Formats

Solving A Gas-Lift Optimization Problem by Dynamic Programming

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Solving A Gas-Lift Optimization Problem by Dynamic Programming

Uploaded by

Copyright:

Available Formats

European Journal of Operational Research 174 (2006) 12201246

Solving a gas-lift optimization problem

Fig. 1. A cluster of oil wells.

po cno pg cng pw cnw qno qni

N is the number of oil wells;

qno an0 y n an1 qni an2 qni an3 qni ;

qno minfqno qni 2 Qn g;

Proposition 1. PE(G) is NP-Hard.

we conclude that the direct and exhausN 1

qno an0 y n an1 qni an2 qni an3 qni ;

Proposition 2. If ln < un 6 qmax for each n, then dimP 2N .

is a valid inequality for P. Further, the face F C fqi ; y 2 P :

jCj 1g induced by the cover

Subtracting the latter from

which is valid for P. The extended

qno an0 y n an1 qni an2 qni 2 an3 qni 3 ;

ln y n 6 qni 6 minfun ; mdgy n ;

P m;N : J m;N F m;N ;

CPU time (seconds)

4.3. Computational results

CPU time (seconds)

Step 2. Set n = vk and let ~z ~

general valid inequality for PG inducing the maximal face F p fqi ; y 2 PG :

which is valid for PG.

Subject to: wn 2 f0; . . . ; mg;

ln y n 6 qni 6 minfun ; dwn gy n ;

ln 6 qni 6 minfun ; dwn g;

CPU time (seconds)

Running time (seconds)

You might also like