A Mixed-Integer convex formulation for production optimization of gas-lifted oil fields with routing and pressure constraints

PROCESS SYSTEMS ENGINEERING

A Mixed-Integer convex formulation for production optimization of gas-lifted oil fields with routing and pressure constraints

M. A. S. Aguiar^*,**; E. Camponogara; T. L. Silva

Department of Automation and Systems Engineering, Federal University of Santa Catarina C. P. 476, 88040-900 Florianópolis - SC, Brazil

ABSTRACT

Production optimization of gas-lifted oil fields under facility, routing, and pressure constraints has attracted the attention of researchers and practitioners for its scientific challenges and economic impact. The available methods fall into one of two categories: nonlinear or piecewise-linear approaches. The nonlinear methods optimize simulation models directly or use surrogates obtained by curve fitting. The piecewise-linear methods represent the nonlinear functions using a convex combination of sample points, thereby generating a Mixed-Integer Linear Programming (MILP) problem. The nonlinear methods rely on compact models, but can get stuck in local minima, whereas the piecewise-linear methods can reach globally optimal solutions, but their models tend to get very large. This work combines these methods, whereby piecewise-linear models are used to approximate production functions, which are then composed with convex-quadratic models that approximate pressure drops. The end result is a Mixed-Integer Convex Programming (MICP) problem which is more compact than the MILP model and for which globally optimal solutions can be reached.

Keywords: MINLP; MILP; Mixed-Integer Convex Programming; Oil Production Optimization.

INTRODUCTION

With the increasing demand for fossil energy, the oil industry has looked for new technologies in hardware and software to enable production optimization of oil fields. These are evolving technologies often referred to as smart fields (Williams and Webb, 2007; Moisés et al., 2008). However, before this concept is transferred to the oil fields, significant challenges in science and technology should be overcome. To this end, this paper addresses the problem of optimizing production of oil fields operated with artificial lifting and subject to facility, routing, and pressure constraints.

Large oil fields have several production wells spread over a wide area. The production of clusters of wells is concentrated in a manifold and then transferred by a flow line to a separator, where the multiphase flow is split into gas, oil, and water. In such fields, the reservoir internal pressure may not be sufficiently high to raise oil naturally to the surface, requiring the use of artificial lifting techniques. A common method is gas-lift, which consists in injecting high-pressure gas at the bottom of the well-bore to reduce the counter pressure and thereby induce a multiphase flow to the surface. In such systems, oil production is a function of the lift-gas injected and the pressure at the well head, which is also related to the separator nominal pressure through pressure drops in the production lines. Here arises the challenge of modeling the production function of wells and pressure drop in pipelines, which can be rather complex.

The production optimization problem consists of a Mixed-Integer Nonlinear Programming (MINLP) problem for which the direct application of standard algorithms may not be possible or effective. This problem is only known conceptually because the well-production functions and pressure-drop relations are not explicitly available. The proposed solution approach uses two-dimensional piecewise-linear models for the well-production functions that depend on the lift-gas rate and manifold pressure, while continuous convex functions are used to approximate the three-dimensional pressure-drop functions. The continuous convex approximation has the advantage of being quite compact when compared to three-dimensional piecewise-linear models. In the end, the production optimization problem is approximated as a Mixed-Integer Convex Programming (MICP) problem which can be tackled with off-the-shelf solvers. Computational experiments are performed to compare the MICP formulation and a Mixed-Integer Linear Programming (MILP) formulation obtained by piecewise-linearizing the pressure drops, which is more precise but demands a large number of sample points. The MICP formulation is further compared with an MINLP formulation obtained from MICP by imposing a precise pressure-balance in the flow lines. The computational experiments assess the trade-off between solution speed and the degree of approximation across formulations. A simulation analysis is carried out to compare the mean errors of the field variables predicted with the formulations, against the variables obtained with a multiphase-flow simulator. Because the true MINLP problem is only known conceptually, this analysis serves the purpose of comparing the effectiveness of the approximations and identifying how complex the underlying models need to be to represent the production process satisfactorily.

In what follows, the production optimization problem is first formalized in mathematical notation. Background is presented on multidimensional piecewise-linear (PWL) and quadratic programming (QP) models, which are then used in the problem formulation. Computational and simulation results are reported and discussed. Finally, the paper ends by offering some concluding remarks and suggesting directions for future research.

PROBLEM DEFINITION

In several offshore fields, oil production relies on artificial-lifting methods to compensate for low reservoir pressure and high depth reservoirs, notably the ones located off the coast under the sea bed. Among the artificial-lifting methods, gas-lift is a widely applied technique for its desirable features that include relatively low installation and maintenance costs, robustness for using few mechanical components, and efficiency. It works by injecting high pressure gas at the bottom of the production tubing to induce a pressure gradient from the reservoir all the way up to the surface facilities.

When the operating conditions on the surface are kept constant or change slowly, the modeling of well production can be carried out using Well Performance Curves (WPCs), which relate the production of oil, gas, and water to the lift-gas injection rate. Several works in the technical literature (Buitrago et al., 1996; Alarcón et al., 2002; Camponogara and Nakashima, 2006; Camponogara and de Conto, 2009; Misener et al., 2009; Codas and Camponogara, 2012) have addressed the problem of optimizing production for pre-determined surface conditions.

However, when surface conditions change frequently due to routing operations, failure of equipment, and shutting operations of wells, the standard WPC modeling needs to be extended to consider the pressure at the manifolds which concentrate production from the wells (Kosmidis et al., 2004). Some works from the literature take into account pressure balance constraints (Litvak et al., 1997; Kosmidis et al., 2004; Bieker, 2007; Gunnerud and Foss, 2010; Silva et al., 2012; Codas et al., 2012; Silva and Camponogara, 2014), with some approaches using nonlinear functions and others relying on piecewiselinear models to approximate pressure drops through pipelines.

Unlike the preceding works, we suggest the representation of the pressure drops using convex quadratic functions that are adjusted to measured or simulated data, while approximating well production functions with multidimensional piecewise-linear models. This hybrid approach has the advantage of ensuring convexity of the continuous relaxation of the pressuredrop approximation when combined with flow equations given in piecewise-linear form. (The composition of a convex function with an affine function is convex). Thus, any enumeration scheme such as branch-and-bound applied to the resulting approximation formulation will solve convex relaxations and therefore the global optimum can be reached. On the other hand, the composition of a convex pressuredrop function with a general and even concave production function will not necessarily be a convex function, thereby making global optimization computationally hard (Boyd and Vandenberghe, 2004). Notice that the well production function is typically concave for a varying lift-gas injection rate, provided that the well-head pressure is kept at the nominal value and the well does not have a kick-off rate.

The definition of the Production Optimization Problem (POP) is based on the following parameters:

N = {1,..., N} is the set of oil wells, N being their number, and N_m ⊆N is the subset of wells that have a flow line connected to manifold m;

M ={1,... , M} is the set of manifolds, M being the number of manifolds, and M_n ⊆ M is the subset of manifolds that can receive production from well n , which is then transferred to a single separator;

H = {o,g,w} has the multiphase flows: oil (o), gas (g) and water (w);

q_i^max is the lift-gas limit supplied by the compressors;

q_ⁿi,^min and qⁿ_i^max are operational limits for lift-gas injection into well n;

p^m,S is the nominal pressure of the separator that receives production from manifold m;

q^n,L and q^n,U are the lower and upper bound on the flow rate of well n;

variables:

q_iⁿ is the lift-gas rate allocated to well n;

y_n is 1 when well n is producing, and 0 otherwise;

z_n,m takes value 1 if the production from well n is directed to manifold m , and 0 otherwise;

q^n,m_h^, is the flow of phase h ∈H sent from well n to manifold m and q^n,m = (q_h: h ?H) is a vector with all phase flows. The gas flow rate received by the production manifold is the sum of the lift-gas injected into well n (Inj) and the gas from the reservoir (R): = + ;

q^mΣ _n∈Nm is the total flow received from the wells connected to manifold m for all phases;

p^m is the pressure of manifold m; and functions:

g is the profit function of the oil and gas production that accumulates in a manifold;

c is the cost function for lift-gas injection into a well;

q^n,m (p^m, q_iⁿ) is the production function of well n if connected to manifold m as a function of the manifold pressure p^m and the lift-gas injection rate q_iⁿ;

Δp^m (q^m) is the pressure drop in the flow line connecting manifold m to its dedicated separator as a function of the multiphase flow q^m.

Then the production optimization problem considering a group of gas-lifted wells, routing decisions about well-manifold connections, and lift-gas, separator flow handling, and pressure constraints can be conceptually formalized as an MINLP problem:

The problem aims to maximize the objective f which is composed of a gain function g (financial gains from oil and gas flows represented by q^m discounted the cost of water treatment before discharge) and the cost c incurred by compression and injection of lift-gas in each well.

The restriction given by Eq. (1b) says that there is a limited amount of lift-gas. Eq. (1c) says that the gas injection in an individual well is limited by an upper and lower bound if the well is producing when y_n=1, otherwise this injection is zero when y_n=0. Eq. (1d) ensures that each well must be connected to a single manifold when the well is producing. Eq. (1e) defines the multiphase flow q^m,n from well n to manifold m as a function of the manifold pressure p^m and the lift-gas injection q_iⁿ. Eq. (1f) restricts the multiphase flow of a well n to be within lower and upper bounds when the well is operating (y_n=1). Eq. (1g) ascertains that the inflow into a manifold m is precisely the sum of the multiphase flows from the connected wells. Eq. (1h) establishes processing capacity for the separators. Eq. (1i) ensures the pressure balance in the flow lines, namely that the difference between the manifold pressure and its separator is precisely the pressure drop through the flow line.

The MINLP problem given in Eq. (1) is a conceptual definition of the production optimization problem, since the well-production function q^n,m(p^m, q_iⁿ) and the pressure drop Δp^m(q^m) are not known explicitly. Correlations can be found in the literature and are often available in simulation software for approximating pressure relations in oil production systems (Beggs and Brill, 1973; Litvak and Darlow, 1995). Although one could conceivably model q^n,m and Δp^mwith these correlations, the resulting MINLP formulation would be highly nonlinear and complex, rendering a global optimization problem, which is a challenge for existing algorithms and software. This work is half-way between the approaches that use piecewise-linear models and those that rely on nonlinear correlations. Instead, the well-production functions will be approximated with piecewise-linear functions and the pressure drops with convex functions.

PIECEWISE-LINEAR AND QUADRATIC CONVEX APPROXIMATIONS

Two favorite strategies for solving the production optimization problem are nonlinear programming methods and MILP strategies which rely on models obtained by piecewise linearizing the nonlinear functions. While the former method is prone to get stuck in local minima, the latter method can lead to very large MILP problems.

This work suggests a hybrid approach that approximates the pressure drop functions with convex functions and the well-production functions with piecewise-linear models. Such an approach renders the approximation problem an MICP program which is far more compact than the MILP approximation and which can be solved efficiently up to optimality.

Convex Combination Model

Among the MILP models for piecewise linearization available in the literature (Vielma et al., 2010), this work uses the Aggregated Convex Combination (CC) model to approximate the well-production function q^n,m. Let f: → be a continuous function defined over a compact domain ⊆^d. According to Vielma et al. (2010), f is piecewise-linear if and only if there exists a family of polytopes , such that P∈ P =, , {C_P} P∈ and further:

Let V(P) be the set of vertices of polytope P and () be the set of all vertices. The CC model assigns weighting variables to each vertex v∈(). Thus a graph point is represented by (x,f(x)) = (v,f(v)) where and = 1. The CC model is given by:

where (v):={P∈:v V∈(P)} is the set of polytopes that contain vertex v.

The constraints (3a) represent a graph point (x, f (x)) as the convex combination of the vertices and their associated function values. Equations (3b) ensure that the λ weights define a convex combination of the vertices and function values. Equations (3c) limit convex combinations to a single polytope, and further guarantee that only the weighting variables associated with vertices of the active polytope can be non-zero.

Convex Quadratic Model

A convex quadratic approximation of a nonlinear function f:ⁿ → consists of

where Q is a positive definite matrix denoted by Q 0, b is a vector, and c is a constant. Such quadratic functions will be used later to approximate pressure-drop relations, since their behavior is dominated by the effects of friction, which depends on the square of the flow speed.

PROBLEM APPROXIMATION

This section begins by presenting the piecewiselinear modeling of well production, followed by the convex quadratic approximation of pressure drop and the synthesis of such models. These models are then combined to obtain the MICP formulation for oil production optimization.

Piecewise-Linear Approximation of Well Production Functions

The multiphase flow function q^n,m (p^m, q_iⁿ) depends on the pressure of the manifold m to which this well is connected, p^m, and the lift-gas injection rate q_iⁿ. Using the CC model, the piecewise-linear approximation of well production is formulated as follows:

For all n ∈N:

having the following extra parameters:

^n,m and ^n,m and are the sets of breakpoints for the lift-gas rate and manifold pressure when well n is connected to manifold m, respectively;

^n,m is the set of polytopes with vertices in ^n,m × ^n,m;

^n,m (q_i,p_r) = {P∈^n,m :(q_ip_r) ∈V(P)} has the polytopes that contain vertex (q_i,p_r);

p^m,max is the maximum manifold pressure; and extra variables:

is the weighting variable of a breakpoint pair in ^n,m × ^n,m;

is a binary variable for each polytope P∈^n,m which assumes value 1 if the convex combination is limited to polytope P . According to constraints (5g)-(5i), only the vertices of P can be part of the convex combination defining the gas injection into well n and manifold pressure p^m;

^n,m is the piecewise-linear approximation of q^n,m.

Notice that the injection bound (1c) and the wellproduction restriction (1f) are implicitly imposed by the piecewise-linear approximation, i.e., the infeasible points do not belong to the domain of the PWL approximation function.

Convex Quadratic Approximation of Pressure Drop Functions

The approximation of the pressure drop in the flowline from a manifold m to its separator is approximated with a convex quadratic function. Thus, for all m ∈:

where Q_m 0. Being an equality, Eq. (6c) induces a nonconvex and discontinuous set of candidate solutions, which results in a nonconvex approximation of the production optimization problem. However, this equality may be approximated by two convex inequalities:

with being a convex underestimation and being a concave overestimation of the pressure drop Δp^m. Replacing (6b)-(6c) with (7) would lead to a relaxation of the true MINLP. Notice that a problem R:Z_R =max{f_R (x):x∈X(R)} is a relaxation for a problem P:z_P =max{f_R (x): x∈X(R))} if X (R)⊇X(P) and f_R(x) > f_P(x) for all x ∈X(P). Provided that the piecewise-linear models of the production functions are precise, the optimal solution to the MICP arising from this replacement would induce an upper bound on the objective. The ability to produce relaxations based on under and overestimation within a reduced domain of the decision space, possibly using convex and concave functions as given in (7), would allow the application of a spatial branch and bound strategy.

The physical behavior in the oil production system is such that, for a given pressure difference (p^m-p^m,S) the flow q^mwill be as large as possible and, thereby, so will be the pressure drop Δp^m(q^m). This means that constraint (7b) may not be bounding and the approximation of (6b)-(6c) can be carried out only with (7a). Since the effectiveness of this approximation should be assessed by simulation of the oil field, we will consider a general single-sided approximation of the form:

where

is a convex function that meets one of the following cases:

1) Underestimation: if (q^m) < Δp^m (q^m) for all q^m, then (6b)-(6c) can be approximated by:

meaning that the resulting MICP formulation will be a relaxation for the true problem, provided that the piecewise-linear functions for well production can be regarded as precise models and the objective function is not modified.

2) Overestimation: if (q^m) > Δp^m (q^m) for all q^m, then (6b)-(6c) is approximated by:

3) Estimation: when (q^m) (q) neither underestimates, nor overestimates the true pressure drop, an approximation for (6b)-(6c) is given by:

The excess ε^m may be necessary to ensure feasibility of the MICP formulation when is an overestimation. Notice that the excess ε^m can be nullified when < p^m - p^m,S otherwise ε^m should be driven to zero. One possibility is to introduce in the objective a penalty factor on the excess, such as with µ > 0. As µ tends to zero the penalty increases, drawing ε^m towards the origin unless a feasible solution with ε^m = 0 does not exist. The strategy of introducing an excess variable is also viable for the case of overestimation.

To obtain the MICP formulation, the choice among the cases above should be based on an analysis of the approximation quality and validated through simulation.

Curve Fitting for Quadratic Approximation of Pressure Drop

Two relevant issues are how convex-quadratic approximations for the pressure drop are computed and whether or not such approximations are satisfactory. To resolve these issues, we suggest solving a Semi-Definite Program (SDP) minimizing the error with respect to a set of pressure-drop points obtained from field data or multiphase-flow simulators. The synthesis of the convex-quadratic approximation for Δp^m consists in solving the following problem:

where ={1,...,K} is the set of indexes of the pressure-drop sample points (q^m,1,Δp^m,1),...,(q^m,K, Δp^m,K). Notice that CF can be easily recast as an SDP by linearizing the objective function with the introduction of slack variables and linear inequalities, namely by expressing the l₁-norm of the vector of approximation errors as a system of linear inequalities. Reformulations of CF in SDP are briefly discussed in the Appendix.

The choice of convex approximation, the l₁-norm, and the relative error in the objective function was not arbitrary but rather the result of experimental analyses. For instance, the pressure-drop approximations obtained by minimizing concave functions were consistently worse than the minimization of convex functions. The experiments also revealed that the error should be normalized because pressure drops can vary drastically from low to very high values.

The analyses leading to these findings will be presented later in the paper.

Piecewise-Linear Convex-Quadratic Approximation

By piecewise-linearizing the well-production q^n,m function using the CC model and representing the pressure drop Δp^m with a convex-quadratic function, the oil production optimization problem is approximated by an MICP problem:

The approximation error for the well-production functions is controlled by the number and location of breakpoints. On the other hand, the pressure-drop error depends on the quadratic models obtained by SDP optimization, but also on the penalty factor µ. Let w(µ) be a vector with the solution to POP(µ). Our solution strategy consists in solving a problem sequence {POP(m_k)} with decreasing µ_k until approximation errors ε^mare sufficiently small. At iteration k, problem POP(µ_k) is solved starting from w(µ_k-1) which is feasible for POP(µ_k).

Notice that the (µ_k) formulation encompasses the three cases considered for . If is an underestimation, then ε^m will be driven to zero by the suggested iterative strategy and POP(µ_k) will induce an upper bound when ε^m =0 for all m. In this case, the excess variables and the penalty factors can be removed altogether, allowing POP(µ_k)µ_k) to be solved only once since it will be unaffected by the choice of µ. On the other hand, when overestimates the true pressure drop, POP(µ_k)µ_k) will be a handy approach to find a nearly optimal solution when the problem becomes infeasible without the excess variables.

Fully Piecewise-Linear Approximation

The MILP approximation for the oil production optimization problem is obtained from POP(µ_k) eliminating the excess variables ε^m, along with the penalty factor in the objective function, and replacing (13n) with a three-dimensional piecewiselinear formulation using the CC model:

with the parameters:

is the fluid flow of oil (kº), gas (k^g) , and water (k^W) of the vertex (kº, k^g, k^W) in the domain of Δ p^m;

m is the set of polytopes whose union defines the domain of the pressure-drop function;

^m (kº,k^g,k^W) ⊆ ^m is the set of polytopes that contain vertex (kº,k^g,k^W);

() is the se of all vertices appearing in the set of polytopes ^m and the additional variables:

is a binary variable associated to polytope P∈^m which assumes the value 1 when the convex combination is confined to P;

Ω^m_kº,k^g,k^W is a variable with the weight associated to vertex (kº,k^g,k^W) ∈();

y^m is an auxiliary binary variable which takes on the value 1 when manifold m receives production.

As before, the bound constraint (13m) is implicitly imposed by the domain of the PWL approximation of the pressure-drop relation.

COMPUTATIONAL ANALYSIS

This section evaluates the MICP formulation computationally. First, a synthetic oil field is instantiated in a commercial multiphase-flow simulator honoring the characteristics of real-world oil fields. Breakpoints for well-production functions and pressure drops in flow lines are obtained by sampling the simulator, which are used to yield piecewise-linear models for q^n,m and quadratic approximations for Δp^m Finally, the performance and solution quality of the MICP formulation is compared with the MILP formulation and an MINLP formulation obtained from MICP.

Oil Field Scenario

Inspired in a scenario from (Kosmidis et al., 2004), a synthetic oil field was put together for the purpose of computational and simulation analysis (Silva et al., 2012). This field consists of two separators with a limited separation capacity and an operational pressure of 300 psia (p^m,S).

Figure 1 illustrates the structure of the gathering system of the synthetic oil field. Separator 1 is connected to an adjoint manifold by a flow line of 100 m, whereas separator 2 is connected to another manifold by a flow line of 50 m. All of the 16 gas-lifted oil wells can be routed to one of the manifolds; however, wells 1 to 8 are closer to manifold 1, while wells 9 to 16 are closer to manifold 2. The wells are limited in the amount of fluids they can handle and lift-gas injection rates should be within bounds when the wells are producing.

The oil field was modeled with the Pipesim simulator, which allowed us to obtain breakpoints for piecewise-linearization of the well-production functions and pressure-drop points for the synthesis of quadratic approximations.

Approximation Analysis of Pressure Drop

Several strategies were tested in order to find the model which fits the sample data best. The minimization of the l₁-, l₂ -, and l_∞ -norm of the relative and absolute error vectors were considered. Because the visual analysis was not possible, the concaveand convex-quadratic fittings of the sample data (pressure drop) were analyzed. The different strategies considered for pressure-drop approximation with quadratic functions are shown in Table 1.

To find the approximations according to the strategies of Table 1, the corresponding variations of the curve fitting problem CF , given in Eq. (12), were implemented in YALMIP (Löfberg, 2004) and solved with the SDP solver SeDuMi version 1.3 (Peaucelle et al., 2002). The curve fitting problems were solved in a workstation running Ubuntu Linux with 64 bits. Eq. (12c) was replaced with Q_m = Q^T_m= 0 for the concave approximation, while Eq. (12a) was modified depending on the error type and norm. In all, 729 breakpoints consisting of the combination of 9 gas-flow rates (in the range from 137,777 to 1,240,000 m³/d), 9 oil-flow rates (from 400 to 3,600 m³/d), and 9 water-flow rates (from 92.2 to 830 m³/d) were sampled from Pipesim.

Intending to compare the curve fitting strategies, we established four quantitative indexes consisting of the maximum and mean values for the absolute (ª) and relative error (^%) induced by the resulting quadratic formulations. Table 2 shows these quantitative indexes with absolute errors given in psia and relative errors given in percentage.

The experiments are characterized by a tuple 〈m₁,m₂, m₃〉 where m₁ ∈ {convex, concave} indicates the curvature, m₂∈ {l₁, l₂, l_∞}indicates the error vector norm, and m₃ {minª, min^%} indicates minimization of absolute or relative error.

Analyzing the error indexes, it can be noticed that the strategies with concave curvature yield a worse fit than their convex counterpart: 〈concave, m₂, m₃〉 is worse than 〈convex m₂, m₃〉 for all m₂ and m₃.

Further, the convex strategies that minimize relative error tend to induce a better fit than the convex strategies that minimize absolute error: 〈convex,m₂ ,minª〉 is worse than 〈convex ,m₂, min^%〉 for all m₂. The remaining issue is the choice of the error norm m2 for% the strategy 〈convex,m₂,min^%〉 . Because the mean errors obtained with the strategy 〈convex,l_∞, min^%〉 are relatively high, the remaining choices are l1 and l₂. The strategy was chosen for the computational experiments because l1 produced lower mean error in absolute and relative terms.

Computational Analysis

The MILP and MICP formulations were programmed in AMPL and solved with CPLEX 11 in an Intel Core 2 Quad at 2.93 GHz, running on a 64-bit Linux workstation, with 4 GB of RAM. An MINLP formulation was obtained from MICP by imposing the pressure balance at equality on the flow lines (i.e., replacing inequality (13n) with an equality) was also programmed in AMPL, but solved with the global solver SCIP 3.0 (Achterberg et al., 2008; Achterberg, 2009) on the same workstation. All experiments ran with a time limit of 10 000 seconds (≈2.8 hours).

The lift-gas availability was varied in three different situations: the compressor has full capacity in the first case (High), half in the second (Medium), and only the capacity for maximizing the production of a single well (Low). This variation in compressing capacity aims to assess the relative performance of the formulations under disparate operational conditions.

The quality of approximation varies according to two different resolutions: a moderate number of polytopes in the PWL functions (Moderate), and a considerably high number of polytopes (Fine). The Moderate resolution has 18 polytopes (squares) for the WPC curve (6 breakpoints for injection rate, and 3 for manifold pressure), and 125 polytopes (cubes) for the pressure-drop function (5 breakpoints for each phase flow). The Fine resolution has 66 polytopes (squares) for WPC curves (11 breakpoints for injection rate, and 6 for manifold pressure), and 10³ polytopes (cubes) for pressure-drop functions (10 breakpoints for each phase flow). The goal for varying the resolution is to evaluate the trade-off between the quality of approximation and the relative performance of the formulations.

For the sake of simplicity, we assume that =_n, the number of breakpoints for qⁿ_i is = K and the number of breakpoints for p^m is =R in the well-production functions, which imply =(K-1)(R-1) for all m and n. Further, the number of breakpoints in each dimension (oil, gas, and water) for Δp^m is W which implies =(W-1)³ for all m. Table 3 gives the size of the conceptual MINLP, MILP and MICP formulation as a function of these parameters, along with the actual size of the Moderate and Fine instance. The size of the MICP and the MINLP formulation obtained by imposing an equality on the pressure balance are the same.

For solving POP(µ) , an iterative process was followed in which the penalty factor was initialized as µ⁰ =1, and updated as µ^k+1 =µ^k/10, until ^m became lower than 10^-6. Table 4 shows the execution time in seconds and the final dual gap of the solutions obtained with these formulations. For MICP, the table also provides the total number of iterations performed and the maximum slack on the pressure constraint (13n), i.e.,

The computational experiments revealed that the production optimization problem may not be efficiently solved with the MILP and MINLP formulation:

The MILP formulation reached the global optimum with the Moderate resolution for all compression capacities. With a Fine resolution, it was solved to optimality for high capacity, but failed to close the primal-dual gap with medium and low compression capacities.

The MINLP formulation was solved more efficiently than the MILP formulation for the Moderate resolution, but the primal-dual gap could not be closed with the Fine resolution.

On the other hand, globally optimal solutions were found relatively quickly using the MICP formulation for all instances, arguably due to the gradient information provided by the convexquadratic approximations and the reduced number of binary variables.

It can also be noticed that the approximation of the pressure-drop equation, = p^m - p^m,S with inequality (13n) was very accurate. This behavior corroborates the hypothesis that this inequality is binding, with the slack in the inequality being essentially zero.

Simulation Analysis

This section presents an approach for analyzing and reducing mean errors between the values of process variables calculated with the simulator and their predictions obtained with the MILP and MICP formulations. The solutions obtained with the MICP and MINLP formulation were essentially identical. This approach called off-line simulator-optimizerapplication loop is shown in Figure 2.

Aiming to reduce the discrepancy between optimization models and the process simulation model, the off-line loop applies four steps:

Step 1: The application gives the simulator an initial resolution quality.

Step 2: The nonlinear functions are sampled in the multi-phase flow simulator.

Step 3: The optimizer receives the sample breakpoints as inputs and finds a solution providing liftgas rates, well-manifold routes, flow rates, and pressure predictions.

Step 4: The lift-gas rates and well-manifold routes obtained by the optimizer are injected into the simulator. Then, the values calculated by the simulator and the optimizer predictions are given as inputs to the application. A mean error is calculated for the current resolution and the application decides whether a new iteration is necessary for the optimizer predictions to match simulator values. If the resolution quality needs to be improved, a new iteration starts from step 1.

The simulation analysis evaluated the mean errors of the scenario with low compression capacity for both optimization models (MILP and MICP) considering two resolution iterations: Moderate and Fine resolutions. Tables 5 and 6 present the simulation values and the relative errors of the optimization model predictions, for Moderate and Fine resolutions respectively. The MILP and MICP formulations produced the same well-manifold routes. The differences arise in the multiphase flows handled by the manifolds and their pressures.

The mean of the mean errors of the variables presented in Table 5 is 8.96% for MILP-Moderate and 9.53% for MICP-Moderate. MILP-Moderate is slightly more precise than MICP-Moderate. The predictions for water flows in manifold 2 were considerably different than the values calculated by the simulator (28.51% and 31.23%). This might be at-tributed to an insufficient number of breakpoints to represent the nonlinear characteristics of well-performance and pressure-drop curves.

According to Table 6, the simulator-relative errors of the Fine resolution are relatively low when compared to the Moderate resolution. The mean of the mean errors is 1.02% for MILP-Fine and 3.5% for MICP-Fine. Once the resolutions of the MILP and MICP are improved, the simulator-relative errors decrease and good solutions are reached with both formulations.

The MILP-Fine formulation is more precise and achieves a representation quite similar to the real process according to the simulator.

The computational and simulation analyses elicit the following remarks:

Performance: The computational analysis shows that an improvement in the resolution quality has a significant impact on the performance of the MILP formulation. When the lift-gas delivered by the compression unit is constrained to medium or low, the MILP formulation could not reach the global optimum with the Fine resolution. On the other hand, an improvement in the resolution quality does not slow down the MICP formulation significantly. This result was expected since the MICP formulation has the advantage of being more compact-the three-dimensional functions (pressure drops) are approximated with a quadratic function.

Process Representation: The analysis indicates that both MILP and MICP formulations reach better solutions when the approximation resolution is improved (the number of sampled breakpoints is increased). In order to reduce the optimizer prediction errors, some iterations of the off-line simulator-optimizer-application loop will be necessary to determine a sufficient number of breakpoints for a satisfactory approximation of the nonlinear functions.

Extrapolation of the Sampled Region: The MICP formulation was bound to the same sampled region of the MILP formulation. However, the MICP approach is able to extrapolate out of the sampled region within which the convex-quadratic approximation was fitted. It might not be a good approximation, but at least it gives us a lead about the existence of better operating points.

SUMMARY

This work proposed a mixed-integer convex programming formulation for oil production optimization of gas-lifted oil fields subject to operational, routing, and pressure constraints. The MICP formulation arises from the piecewise-linearization of the well-production functions and their composition with convex-quadratic approximation functions of the pressure drops. Besides being relatively compact, the convexity of the MICP formulation allows for optimization solvers to reach globally optimal solutions. The piecewise-linear models arose from the convex combination of well-production points sampled in the domain of lift-gas injection and manifold pressure, whereas the convex-quadratic models were obtained by solving a semi-definite program for fitting the model to pressure-drop points sampled from a multiphase-flow simulator.

A computational analysis performed in a synthetic oil field modeled in a multiphase-flow simulator showed that the MICP formulation is more efficient than the MILP and MINLP formulation.

A simulation analysis performed in the same oil field showed that the MILP and MICP formulations reach better solutions when the approximation resolution is improved. Although the MILP formulation is more precise, it might not be efficient when the approximation resolution is fine. The MICP formulation approximates the three-dimensional functions (pressure drops) with a single quadratic function resulting in a more compact formulation. This MICP approach has the advantage of not being significantly slowed down when the resolution quality is improved.

Future research directions include:

the design of more detailed models for pressure drop, possibly including temperature and outlet pressure in the pipelines;

the development of piecewise-convex models, which could more precisely approximate complex functions while being more compact than piecewiselinear models; and

the integration with compressor scheduling models (Camponogara et al., 2011; Camponogara et al., 2012).

ACKNOWLEDGMENTS

This work was funded in part by a research contract with Petróleo Brasileiro S.A. (Petrobras).

NOMENCLATURE

Parameters

N number of wells M number of manifolds q _i ^max

maximum lift-gas rate which represents the compressor capacity

q _i ^n,min lift-gas injection lower boundi q _i ^n,max lift-gas injection upper bound p^m,S nominal pressure for separator q ^n,L

lower bound on well flow rate

q ^n,U

upper bound on well flow rate

q ^m,min

minimum manifold flow rate

q ^m,max maximum manifold flow rate p^m,max maximum manifold pressure Sets

set of wells

set of manifolds

_m subset of wells that can be connected to manifold m _n manifolds that well n can send its production

set of phase flows

set of polytopes

V(P) set of vertices of a polytope ( ) set of all vertices P(v) set of polytopes that contain vertex v ^n,m breakpoints for the lift-gas rate when well n is connect to manifold m ^n,m manifold pressure breakpoints when well n is connected to manifold m ^n,m

set of polytopes with vertices in ^n,m × ^n,m

^n,m(q_i, p_r) subset of polytopes that contain vertex (q_i, p_r) ^m set of polytopes in the domain of the pressure drop function ^m(kº, k^g, k^w) subset of polytopes containing vertex (kº, k^g, k^w) (^m)

set of all vertices in the domain of the pressure drop function

Variables q _i ⁿ gas injection rate for well n y_n well n activation Z_n,m connection of well n to manifold m

flow of phase h from well n to manifold m q ^m total multiphase flow handled by manifold m p^m

pressure in manifold m

λ_v weight for vertex v in a piecewise-linear approximation y_P polyhedron selection in piecewise-linear approximation

weighting variable for vertex (q_i ,p_r)

binary variable indicating that a polytope P is active

approximation of well production

approximation of pressure drop in the flowline of manifold m ^m

excess variable for approximation error of pressure drop

weighting variable for vertex (kº,k^g, k^w) y^m auxiliary variable indicating that a manifold m is receiving production Functions g objective profit function C objective cost function q ^n,m production function of well n when connected to manifold m Δp^m

pressure drop in the flow line of manifold m

Boyd, S. and Vandenberghe, L., Convex Optimization. Cambridge University Press (2004).

Submitted: October 22, 2012 ; Revised: April 23, 2013 ; Accepted: April 28, 2013

This is an extended version of the manuscript presented at EngOpt 2012 - The 3rd International Conference on Engineering Optimization -Rio de Janeiro, Brazil, 1-5 July, 2012.

APPENDIX

SEMI-DEFINITE REFORMULATIONS

The strategies for finding quadratic approximations are a variation of problem CF given in Equation (12), which can be easily recast as a standard SDP by replacing the objective with:

The l_∞ - norm for approximation error can be handled in a similar manner by simple manipulation of problem CF . It suffices to introduce inequality (A1b), introduce the inequalities , and replace the objective with min .

The representation of the square of the approximation error (l₂ -norm) requires the use of Schur complement since quadratic functions are not explicitly allowed in SDP. The synthesis of convexquadratic models using the l₂ -norm of relative errors can be cast as an SPD program:

According to Schur complement, Eq. (A2c) is equivalent to:

Notice that the fitting of concave-quadratic functions is achieved by using Q_m 0 as the condition on the Hessian. Approximations that minimize absolute error are obtained by eliminating the denominator Δp^m,k from the equations that calculate error.

The gas-flow rates are of the order of 10⁶ and when squared reach values of the order of a trillion of m³/d, which are far larger than the water- and oil-flow rates. Thus, normalization of the flow rates was required before using the SDP solver to find a quadratic approximation.

The normalization is obtained with a diagonal matrix Tm as follows:

where is the normalized multiphase-flow rate and Tm,o(Tm,g, Tm,w) is the maximum oil-flow rate (gas-flow and water-flow rate). Solving problem CF using the normalized flow-rates yields a normalized quadratic approximation function given by . Thus,

by defining, Q_m = , and c_m = . Further, it can be shown that Q_m 0 if and only if _m0, provided that T_m is non-singular.

• Achterberg, T., SCIP: Solving constraint integer programs. Mathematical Programming Computation, 1(1), pp. 1-41 (2009).
• Achterberg, T., Berthold, T., Koch, T. and Wolter, K., Constraint integer programming: A new approach to integrate CP and MIP. In Proceedings of the 5th International Conference on Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems, pp. 6-20 (2008).
• Alarcón, G. A., Torres, C. F. and Gómez, L. E., Global optimization of gas allocation to a group of wells in artificial lift using nonlinear constrained programming. ASME Journal of Energy Resources Technology, 124 (4), pp. 262-268 (2002).
• Beggs, D. and Brill, J., A study of two-phase flow in inclined pipes. Journal of Petroleum Technology, 25(5), pp. 607-617 (1973).
• Bieker, H. P., Topics in Offshore Oil Production Optimization Using Real-time Data. Ph.D. Thesis, Norwegian University of Science and Technology (2007).
• Buitrago, S., Rodruiguez, E. and Espin, D., Global optimization techniques in gas allocation for continuous flow gas lift systems. In Proceedings of the SPE Gas Technology Symposium (1996).
• Camponogara, E. and Nakashima, P. H. R., Solving a gas-lift optimization problem by dynamic programming. European Journal of Operational Research, 174(2), pp. 1220-1246 (2006).
• Camponogara, E. and Conto, A. M., Lift-gas allocation under precedence constraints: MILP formulation and computational analysis. IEEE Transactions on Automation Science and Engineering, 6(3), pp. 544-551 (2009).
• Camponogara, E., Castro, M. P., Plucenio, A. and Pagano, D. J., Compressor scheduling in oil fields: Piecewise-linear formulation, valid inequalities, and computational analysis. Optimization and Engineering, 12, pp. 153-174 (2009).
• Camponogara, E., Nazari, L. F. and Meneses, C., A revised model for compressor design and scheduling in gas-lifted oil fields. IIE Transactions, 44, pp. 342-351 (2012).
• Codas, A. and Camponogara, E., Mixed-integer linear optimization for optimal lift-gas allocation with well-separator routing. European Journal of Operational Research, 217, pp. 222-231 (2012).
• Codas, A., Campos, S. R. V., Camponogara, E., Gunnerud, V. and Sunjerga, S., Integrated production optimization of oil fields with pressure and routing constraints: The Urucu field. Computers & Chemical Engineering, 46, pp. 178-189 (2012).
• Gunnerud, V. and Foss, B., Oil production optimization-a piecewise linear model, solved with two decomposition strategies. Computers & Chemical Engineering, 34(11), pp. 1803-1812 (2010).
• Kosmidis, V. D., Perkins, J. D. and Pistikopoulos, E. N., Optimization of well oil rate allocations in petroleum fields. Industrial & Engineering Chemistry Research, 43(14), pp. 3513-3527 (2004).
• Litvak, M. and Darlow, B., Surface network and well tubing-head pressure constraints in compositional simulation. In Proceedings of the SPE Reservoir Simulation Symposium, San Antonio, USA (1995).
• Litvak, M. L., Clark, A. J., Fairchild, J. W., Fossum, M. P., MacDonald, C. J. and Wood, A. R. O., Integration of Prudhoe Bay surface pipeline network and full field reservoir models. In Proceedings of the SPE Annual Technical Conference and Exhibition, San Antonio, USA (1997).
• Löfberg, J., YALMIP: A toolbox for modeling and optimization in MATLAB. In Proceedings of the IEEE International Symposium on Computer-Aided Control System Design, Taipei, Taiwan (2004).
• Misener, R., Gounaris, C. E. and Floudas, C. A., Global optimization of gas lifting operations: A comparative study of piecewise linear formulations. Industrial & Engineering Chemistry Research, 48(13), pp. 6098-6104 (2009).
• Moisés, G. V., Rolim, T. A. and Formigli, J. M., Gedig: Petrobras corporate program for digital integrated field management. In Proceedings of the SPE Intelligent Energy Conference and Exhibition, Amsterdam, The Netherlands (2008).
• Peaucelle, D., Henrion, D., Labit, Y. and Taitz, K., User's guide for SeDuMi interface 1.04. Technical Report, LAAS, CNRS (2002).
• Silva, T. L., Codas, A. and Camponogara, E., A computational analysis of convex combination models for multidimensional piecewise-linear approximation in oil production optimization. In Proceedings of the IFAC Workshop on Automatic Control in Offshore Oil and Gas Production, Trondheim, Norway (2012).
• Silva, T. L., Camponogara, E., A computational analysis of multidimensional piecewise-linear models with applications to oil production optimization. European Journal of Operational Research, 232(3), pp. 630-642 (2014).
• Vielma, J. P., Ahmed, S. and Nemhauser, G., Mixedinteger models for nonseparable piecewise linear optimization: unifying framework and extensions. Operations Research, 58(1), pp. 303-515 (2010).
• Williams, C. and Webb, T., Shell strives to make smart fields smarter. Society of Petroleum Engineers (2007)

*

To whom correspondence should be addressed

Publication Dates

• Publication in this collection
  07 July 2014
• Date of issue
  June 2014

History

• Received
  22 Oct 2012
• Accepted
  28 Apr 2013
• Reviewed
  23 Apr 2013