INTRODUCTION
Oil production rate from a production well decreases gradually with respect to time. After some period of time, the production rate may fall to a level that economically no longer profitable. This condition encourages engineers to apply some artificial lift methods to increase oil production and one of them that commonly used is gas lift.
During the lift process using gas lift, gas is injected at selected point into the tubing to lighten the fluid column along the tubing, in order to enable the fluid to be delivered to the surface. Oil production normally increases as gas injection increases. However, too much gas injection will cause slippage where gas phase moves faster than liquid. This condition leads to reduction in oil production. Identifying optimal gas injection that maximizing oil production is the main interest in gas lift optimization problem.
In an oil field, oil is produced from a cluster of gas lift wells. The schematic
flow of a cluster of gas lift wells is shown in Fig. 1. Oil
and gas that are produced from each well are delivered to a production manifold
and then collected in a separator.

Fig. 1: 
Schematic flow of a cluster of gas lift wells 
After separation process, oil is distributed to the oil pipe line network.
Some gas is reused for injection and the remaining gas is for sale. Gas available
for injection is usually very limited and should be allocated in optimal form
for each well. Here the optimization problem to solve is: How to share the injected
gas to each well in order to yields maximum total oil production rate?
The optimization problem mentioned above is a complicated problem to which
researcher have pay attention for long time. For each well, liquid production
rate is a nonlinear function of gas injection rate, which is not known explicitly.
In existing approaches, the optimization problem has been solved in three steps
of procedure. The first step, a set of data that relate gas injection and oil
production rate from each well are collected. The data are obtained from empirical
or numerical simulation data. The second step, a regression or interpolation
method is applied to estimate the nonlinear function which relates gas injection
rate to liquid production rate. Some functions commonly used for regression
are quadratic (Nishikori et al., 1995), combination
between quadratic and logarithmic (Alarcon et al.,
2002) and exponential function (Sukarno et al.,
2006). Some researchers apply interpolation method to estimate the nonlinear
function using piecewise linear function (Ray and dan Sarker,
2007a, b; Camponogara and Nakashima,
2006). The third step, the constrained optimization problem is solved numerically
using non linear programming methods or other methods such as genetic algorithm.
For each well, liquid production as a non linear function of gas injection is obtained implicitly from gas lift model. The mathematical model for gas lift problem could be written as a two parameter family of a non linear differential equation (ODE):
that represents the steady state flow (gas and liquid) equation along the tubing, with wellhead pressure:
and bottom hole pressure:
as the boundary conditions. The real valued function is
non negative:
q_{g} and q_{l }where 0≤q_{g}≤1 and 0≤q_{l}≤1
are gas injection rate and liquid production rate respectively. Here F is pressure
gradient function that depends on well geometry (such as depth and tubing diameter)
and fluid properties (such as density, viscosity) in each well. The existence
and uniqueness of liquid production rate as an implicit function of gas injection
rate has been shown by Saepudin et al. (2007).
Also, a computation scheme using genetic algorithm to find optimum gas injection
rate has been proposed by Saepudin et al. (2007)
and also by Sukarno et al. (2009) for dual gas
lift well system. The last computation approach mentioned has reduced the collecting
data and regression or interpolation procedure as required in the previous approaches.
In this study, development of mathematical model and computation scheme by
Saepudin et al. (2007) and Sukarno
et al. (2009) are proposed to solve the optimization problem in a
cluster of gas lift wells system.
MATHEMATICAL MODEL
In an oil production well, reservoir fluid, which will be assumed consists of oil and water, flow from reservoir through tubing to the surface. When gas lift is applied, the reservoir fluid will mix with the gas injection and a two phase fluid flow (liquid and gas) takes place along the tubing. In the simplest form, a mathematical model for gas lift problem could be expressed as a combination of fluid flow in the reservoir and along the tubing.
Single gas lift well model: Assuming injection point near the bottom
hole and reservoir fluid consists of oil and water, the mathematical model for
gas lift problem in normalized form could be written as the boundary value problem
(1)(3). The Eq. 1 is derived from the mechanical energy balance
equation (Economides et al., 1994):
where, the terms:
correspond to the pressure drop due to gravity, friction and acceleration,
respectively. Mixture density and
velocity
are function of pressure P, gas injection rate q_{g} and liquid production
rate q_{l}. Since the liquid and gas superficial velocities are given
by:
and the cross section area of the tubing is:
then Eq. 5 could be written as:
In this study, the two phase fluid flow model along the tubing will be given
by a simple homogeneous model with slip, known as driftflux model. Although
this model is not as accurate as mechanistic models, it has advantages of being
relatively simple, continuous and differentiable (Shi et
al., 2005).
The in situ average density
in Eq. 7 can be expressed as:
where, gas density p_{g} that depends on the pressure P is given by:
The average void fraction y_{g} is given by driftflux model of Zuber
and Findlay (Shi et al., 2005; Guet
and Ooms, 2006):
The parameter distribution C_{0} accounts for the effect of the non
uniform distribution of both velocity and concentration profiles. Driftflux
velocity U_{d} accounts for mean relative velocity between two phases.
There are several driftflux correlations offering procedures to compute C_{0}
and U_{d} some of them are written by Shi et
al. (2005). The friction factor f is a function of Reynold Number and
the pipe roughness. In this study, f is assumed constant.
In the reservoir, liquid production rate in steady state can be expressed by the following equation:
where, P_{wf} and P_{r} represent bottom hole pressure and reservoir pressure, respectively and J is called productivity index.
Let, L be the injection depth (since the injection point is assumed near the bottom hole), then the injection depth is considered equals to the well depth. The gas lift performance model for single well is given by the differential Eq. 7 with boundary conditions are given by:
and
where, the right hand side of Eq. 13 is given by the right hand side of Eq. 11. By scaling:
the Eq. 7 can be expressed as:
Therefore, the gas lift model can be written as the boundary value problem
(13), with the right hand sides are given by the right hand sides of Eq.
15, 12 and 13, respectively. Liquid
production as a function of gas injection q_{g}:
could be obtained implicitly from the gas lift model (13) and the graph of Eq. 16 is known as Gas Lift Performance Curve (GLPC).
Multi gas lift wells model: In many cases of gas lift system, oil is
produced from an oil field which consists of a group of gas lift wells as illustrated
in Fig. 1. Gas for injection is shared for each well through
an injection manifold in appropriate proportion. Considering complexity of the
problem, we assume here the system consists of N gas lift wells, the amount
of injected gas through the injection choke can be controlled as required by
each well, the injection point is near the bottom of the well, the separator
position is close enough to the production manifold (the pressure difference
between separator and manifold can be neglected), flow lines are horizontal
and the separator capacity is large enough (unlimited). The multi gas lift wells
model then could be written as:
The terms P1_{k} and P2_{k} represent pressure in the pipe along the flow line and along the tubing of kth well, respectively. The model is a development of the single gas lift well model (13). The pressure gradient F1_{k} of the initial value problem (1718) which represents gas and liquid flow model along the horizontal flow line of kth well, could be obtained from Eq .15 by dropping the gravity terms and is given by:
Equation 19 represents gas and liquid flow model along the
tubing for kth well, where the pressure gradient is given by:
The term Qg_{av} represents the maximum gas available for injecting the N gas lift wells.
Let J_{m} and Pr_{m} be productivity index and reservoir pressure of mth well satisfying:
Taking:
and
Equation 1721 could be written as:
In the sequels we drop tildes from Eq. 2731.
For each k = 1, 2,..., N, the gas lift performance function of the kth well
is given by:
where, Eq. 32 satisfies Eq. 2731.
OPTIMIZATION PROBLEM
In gas lift well system, optimization problem which commonly studied is: How to maximize total oil production with respect to various constraints such as maximum gas available Qg_{av} to supply injected gas for N gas lift wells and surface facility constraints such as separator capacity. Reservoir liquid is assumed consists of water and oil with proportion given by Water Cut (WC). Therefore oil production rate is expressed as:
Oil production maximization without considering separator capacity: In this case, oil maximization problem could be expressed as a constrained optimization problem:
Determining q_{g}∈D_{1} that maximizes:
where for each k = 1, 2,..., N, function φ_{k} (qg_{k})
satisfies Eq. 2731, with search space:
In case, where total gas available for injection Qg_{av} is large enough, then for each k = 1, 2,..., N, the value of qg_{k} could be chosen such that maximizing φ_{k} (qg_{k}). However, in most cases, Qg_{k} is limited. This is a complicated optimization problem, since for each k = 1, 2,..., N, the function φ_{k} (qg_{k}) is never known explicitly.
For each k = 1, 2..., N, let P2_{k} (z, P; gg_{k}, ql_{k})
be solution of Eq. 2731 and for any given,
the liquid production rate (32) could be obtained implicitly using Eq.
31. Therefore the constrained maximization problem (3435) could be expressed
as:
Determining q_{g}∈D_{1} that maximizes:
with constraints:
Further, solution of the constrained maximization problem (3637) is equivalent with solution of the following minimization problem:
Minimize
in search space:
Oil production maximization with considering separator capacity: In
reality, separator capacity is limited and total liquid production rate possibly
exceeds separator capacity. If this condition happened, gas injection rates
should be recalculated in order the total liquid production rate equals to separator
capacity. Since solutions those satisfy this requirement are not unique, the
solution is chosen that maximizes total oil production and minimizes total gas
injection required. The optimization problem then could be expressed as follow
Determining q_{g}∈D_{2} that maximizes Eq. 34 and minimizes:
with search space:
Using penalty function approach, solution of the optimization problem could be obtained by solving the following minimization problem:
in the domain:
for large enough value of λ_{1}. In the following section, a numerical
scheme is built to solve the optimization problems Eq. 3839
and 4243.
NUMERICAL SCHEME
For each k = 1, 2,..., N and for given qg_{k} and gl_{k} let
(1;
qg_{k}, ql_{k}) be pressure function P2_{k} (z; qg_{k},
ql_{k}) satisfying initial value problem Eq. 2930,
which will be computed numerically using RungeKutta method. Using penalty function
approach (Bazara et al., 1993; Coello,
2002), solution of the following minimization problem:
in domain:
converges to solution of Eq. 3839. And
solution of minimization problem:
in domain:
converges to solution of Eq. 4243, for
large enough λ_{1} and λ_{2}. Solutions of the optimization
problems are computed using genetic algorithm (Goldberg, 1989).
In practice, one prefers a rectangle search space. Using the following transformation:
the domain Eq. 51 could be expressed as a rectangle domain:
And using transformation:
the domain Eq. 47 could be expressed as a rectangle domain:
A numerical scheme based on genetic algorithm is constructed. The domain D_{θ1}
and D_{θ2} are chosen as search spaces to keep population in the
domain as genetic operators work. Computation procedure could be written as
follow:
Step 1: 
Generate initial population v_{1}, v_{2},...,
v_{r} which correspond to {(θ^{(j)}, q_{l}^{(j)}),
j = 1, 2,..., r}⊂D_{θ1} for unbounded separator case and
{(θ^{(j)}, α^{(j)}), j = 1, 2,..., r}⊂D_{θ2}
for bounded separator case 
Step 2: 
Using transformation Eq. 4850 for
unbounded separator case and using transformation Eq. 5257
for bounded separator case, the corresponding {(q_{g}^{(j)},
q_{l}^{(j)}), j = 1, 2,..., r} could be obtained 
Step 3: 
Using 4th order Runge Kutta method, could
be computed 
Step 4: 
The value of or
are
obtained 
Step 5: 
Applying genetic operators (crossover, mutation, selection), new population
are obtained 
Step 6: 
Return to step 2 until convergence criteria are satisfied 
COMPUTATIONAL RESULTS
Here, some numerical results are shown for hypothetic data which is given in
Table 1. Gas lift performance curve for all wells are depicted
in Fig. 2ad.
As we could see from Fig. 2, the GLPC are unimodal curves
where their peak points correspond to the maximum liquid or oil production rates
obtained by gas lift. Relationship between oil and liquid production rate is
given by Eq. 33. Therefore, larger liquid production rate
does not correspond to larger oil production rate. One could check these facts
by comparing Fig. 2a with 2c or comparing Fig.
2b with 2d, where well No. 2 produces the largest liquid
rate, but the oil production rate is less than well No. 1.
Table 1: 
Hypothetic well data 


Fig. 2: 
Gas lift performance curves: (a) liquid prod. vs. gas inj.
rate (normalized) , (b) liquid prod. vs. gas inj. rate (field units), (c)
oil prod. vs. gas inj. rate (normalized) and (d) oil prod. vs. gas inj.
rate (field units) 
If large amount of gas available and separator capacity is large enough, the gas injection could be chosen those maximize liquid production rate. In this case, maximum total oil production rate is 2985.013 STBD (which corresponds to 9327.989 STBD of liquid) and total gas injection required is 7.405 MMSCFD.
For unlimited separator capacity case, numerical simulation is conducted for varying total gas available and the results are written in Table 2 and 3.
We could see from Table 2 and 3, for unlimited
gas injection case, the computational results could be compared with exact solutions
and they show quite agreement.
Table 2: 
Optimum solutions for unlimited separator capacity case (normalized) 

Table 3: 
Converted computational results into field units 


Fig. 3: 
Maximum total oil production vs. total gas available: (a)
normalized and (b) in field units 
Table 4: 
Optimal solutions for limited separator capacity case (normalized) 

Table 5: 
Optimal solutions for limited separator capacity case (field
units) 

Term exact solutions here mean optimum solutions which are obtained by computing
optimal solution for individual well separately.
Relationship between total oil production and total gas injection available is shown in Fig. 3a and b. As we can see from the figure, if the total gas injection is less than 6 MMSCFD, the total oil production rate increases quickly as total gas injection increases and just increases slightly if the total gas injection is larger than 6 MMSCFD. This means that allocating so much gas for injection is not always economically profitable.
For limited separator capacity case, computational results and corresponding conversion into field unit are given in Table 4 and 5. We could see from Table 5 that the largest contribution for oil production comes from well No. 1 which means that the oil production is dominated by well production that produces the largest oil production rate (Fig. 2).
In numerical simulation, population size is 100, crossover probability is 0.9, mutation probability is 0.1 and maximum generation is 500.
DISCUSSION
In future study, optimization problem in gas lift wells system should consider more additional constraints such as gas injection flow from casing annulus into the tubing that leads to the production stability conditions. The flexibility and robustness of the proposed computational method are potential to handle such constraints.
CONCLUSIONS
Some conclusion remarks are obtained through this study.
• 
Oil optimization problem in a cluster of gas lift wells system
could be written in mathematical model as a constrained optimization problem
in a class of boundary value problems 
• 
Genetic algorithm with penalty function constraint handling approach could
overcome complexity of computation in finding optimal solution of gas lift
problem 
• 
This approach gives better quality prediction for optimum solution than
previous approach since all solutions are obtained from model, not from
regression or interpolation 
ACKNOWLEDGMENT
The authors thank the Research Consortium on Gas and Pipeline Network ITB (RC OPPINET) for providing relevant data and field information.