INTEGRATED PRODUCTION SYSTEM OPTIMIZATION USING GLOBAL OPTIMIZATION TECHNIQUES

. Many optimization problems related to integrated oil and gas production systems are nonconvex and multimodal. Additionally, apart from the innate nonsmoothness of many optimization problems, nonsmooth functions such as minimum and maximum functions may be used to model ﬂow/pressure controllers and cascade mass in the gas gathering and blending networks. In this paper we study the application of diﬀerent versions of the derivative free Discrete Gradient Method (DGM) as well as the Lipschitz Global Optimizer (LGO) suite to production optimization in integrated oil and gas production systems and their comparison with various local and global solvers used with the General Algebraic Modeling System (GAMS). Four nonconvex and non- smooth test cases were constructed from a small but realistic integrated gas production system optimization problem. The derivation of the system of equations for the various test cases is also presented. Results demonstrate that DGM is especially eﬀective for solving nonsmooth optimization problems and its two versions are capable global optimization algorithms. We also demonstrate that LGO solves successfully the presented test (as well as other related real-world) problems. discontinuous optimization problem.


1.
Introduction. It is well recognized that many optimization problems in the oil and gas industry are global optimization problems. They may have a large number of decision variables and the objective and/or constraint functions in these problems may contain nonsmooth functions such as maximum, minimum functions and if statements. The presence of such functions leads to nonsmooth and even to discontinuous optimization problem.
So far, different meta-heuristics such as genetic algorithms (GA) and simulated annealing (SA) methods [26] have mainly been used to tackle the nonconvexity in optimization problems from the oil and gas industry. Parallelization of these algorithms have been used to solve large scale problems [18]. Smoothing of discontinuities, which are usually present in many mechanical models, can improve the performance of these global optimization techniques [10], possibly with some degradation in solution accuracy.
However meta-heuristics have many drawbacks. First of all, they require a large number of the objective and constraint function evaluations which is not acceptable when their evaluations are expensive. This makes both simulated annealing and genetic algorithms suitable for only small problems. Second, meta-heuristics have difficulties in dealing with the complicated continuous constraint functions which are common for many optimization problems from the oil and gas industry. Third, meta-heuristics sometimes cannot locate a global solution with high accuracy: as a result they may produce only suboptimal solutions. For large production system optimization problems, the preferred "standard" algorithms are still the Sequential Linear Programming (SLP) [13,24] and Sequential Quadratic Programming (SQP) [9] methods. However, these methods also have their own limitations. In many integrated production system optimization involving gas blending such as in the liquified Natural Gas (LNG) processing, the simultaneous consideration of blending and fluid transmission can result in blending problems which are similar to the Haverly Pooling problem [14]. Problems of this type are generally nonconvex and can have multiple solutions [11]. In this situation, the local solvers may get stuck in one of the many stationary points.
Derivative free methods do not rely explicitly on the local model of the objective and/or constraint functions and therefore they are more effective to solve global and nonsmooth optimization than gradient or subgradient-based methods. In this paper we study the application of different versions of the derivative free Discrete Gradient Method (DGM) as well as the Lipschitz Global Optimizer (LGO) solver suite to production optimization in integrated oil and gas production systems and their comparison with various local and global solvers from the General Algebraic Modeling System (GAMS). It is well known that comparative studies of optimization algorithms can be problematic as there may be objections to the choice of parameters and procedures used (Dolan and Moré as cited in [17]). In this paper we analyze the performance of DGM and LGO applied to typical gas production optimization problems. In integrated production system modeling and optimization, the coupling of subsurface dynamic reservoir models to surface network and processing plants models is non-trivial, with many tie-in possibilities [5]. In this paper, we follow earlier coupling ideas [25] but the minimum and maximum functions are used to combine the pressure/rate transmission equations and gas blending constraints into a single optimization problem. Four test problems, with increasing complexity in the objective and constraints functions are used to illustrate this approach. We present the results of numerical experiments on these test problems.
The paper is organized as follows. Section 2 describes optimization models in the integrated production systems. The description of optimization problems are given in Section 3. Section 4 gives the brief description of algorithms used in this paper and Section 5 presents the results of numerical experiments. We discuss the results of numerical experiments and give some preliminary conclusions in Section 6.
2. Mathematical model of the gas blending process. We consider the process of the blending of gas produced from five fields F 1 , . . . , F 5 to supply different quality gas at six processing plants P 1 , . . . , P 6 through a converging-diverging gas gathering/distribution network, as shown in Figure 1.   This model has five fields, where the average reservoir pressures are fixed and the relationship between the average reservoir pressure and flowing bottom hole pressure (FBHP) can be expressed as a quadratic Inflow Performance Relationship (IPR) (see, for example, [7]). From the bottom of the well to the tubing head, the Vertical Flow Performance (VFP), which associate the pressure drop in vertical pipes with production, is assumed to be quadratic.
In order to find flow rates Q 1 , ..., Q 13 in the above gas network, a set of simultaneous equations is developed that relates the nodes' pressures with flow rates. The governing principles used for fluid flow in a network under steady state conditions are based on Kirchoff's law, and the fact that gas flows down a pressure gradient. It is assumed that flow controllers are present at the gathering centers GC 1 , GC 2 and GC 3 to dissipate energy in order to equalize incoming pressures and spilt flow.

2.2.
Gas inflow performance relationship. The determination of the flowing bottom hole pressure, P F BHP , given the average reservoir pressure, P RES and gas flow rate, Q, is based on the following Inflow Performance Relationship (IPR) equation: where A and F are Darcy and Non-Darcy flow coefficients, respectively. Like most IPR for gas-condensate reservoirs (see, [16] for other and newer IPRs), the dependence (1) is also quadratic. When the pressure is expressed in psia and gas flow rate expressed in MMscf/d (million standard cubic feet per day), the numerical values for flow coefficients of the five fields are shown in Table 1. Table 1. Flow coefficients A and F .

2.3.
Gas vertical flow performance. The vertical flow performance (VFP) for the wells is an adaptation of the Cullender-Smith equation [8]. The equation relates the flowing bottom flow pressure, P F BHP , to the flowing tubing head pressure, P F T HP , and the flow rate, Q, as follows: (2) The VFP flow coefficients for the five wells, for the same pressure and rate units as above, are shown in Table 2.   (1) and (2) imply that the relationship of P T HP with Q for each field is given by 2.5. Pressure drop along horizontal pipes. The pressure drop (P IN − P OUT ) associated with a flow, Q, in a horizontal pipe is calculated from the equation: The coefficient H for each pipe in the network is shown in Figure 2. Like most pressure drop equations [12,23], the above equation is also quadratic. The minimum outlet pressures for P 1 , . . . , P 6 are 500 psia, 550 psia, 600 psia, 580 psia, 560 psia and 670 psia, respectively, as shown in Figure 3. The Tubing Head Pressure (THP) for all wells must be greater or equal to 1543 psia.   The THP provides the energy to drive the gas to deliver at or above the minimum delivery pressure at a plant. Some of the energy is lost due to friction. To illustrate the pressure-rate calculation, we will consider production from F 1 , . . . , F 4 to P 1 and P 3 . The steps from F 1 to P 1 , given Q 1 , a decision variable, are illustrated below: After substitution, P 2 OUT 5 can be expressed as: For feasible flow, P 2 OUT 5 ≥ (500psia) 2 , the squared of the minimum delivery pressure at P 1 , therefore one condition for flow is: The condition for flow to the processing plant, P 3 , is more involved as we have to consider production from F 1 to F 4 and flow and blending at three gathering centers GC 1 , GC 2 and GC 3 .
The outlet pressures from the four fields at GC 3 are calculated as shown above and they are as follows: Then the pressure P GC 3 at GC 3 is: Thus, the condition for flow from GC 3 to P 3 is as follows: or alternatively, In addition to the above equations, at the gathering centers the following mole (and volume) conservation applies: 2.7. Gas composition. It is desired that gas delivered to processing plant P 6 must meet CO 2 and N 2 concentration limits while maximizing LPG content. For the blending calculation, we need the gas composition of all the five fields as they all contribute towards production to P 6 . The gas composition, in volume fraction, for the five fields are shown in Table 3.  Table 3. Field gas composition.
The concentrations of CO 2 , N 2 and LPG at GC 3 are as follows: (22) Thus at P 6 we get and likewise for [ 3. The optimization problems. We will consider four different optimization problems. These problems may differ from each other by both objective and constraint functions. Their general form is as follows: Here the function f may have one of the following forms: 1. The objective function for gas production 2. The objective function for the amount of LPG produced at P 6 : i is the sum of the amount of the propane and butane components. 3. The objective function for maximization of LPG production at plant P 6 and total condensate production: In production operations the upper and lower bounds are generally known and here the decision variables Q i , i = 1, . . . , 13 are in the range between 0 MMscf/day and 200 MMscf/day.
The feasible set X is described by the following constraint functions. There are seven pressure/rate constraint equations to deliver as per minimum outlet pressure requirements shown in Figure 3: There are three volume (mole) fraction constraints for CO 2 and N 2 at GC 3 and LP G P5 : In addition there are three linear equalities from (17) - (19) describing mole (and volume) conservation at the gathering centers: (40) It should be noted that functions ϕ 1 , . . . , ϕ 7 are concave, ϕ 8 , ϕ 9 are nonconvex and ϕ 10 is nonconcave.
We will consider the following problems.
Since f 1 is linear function and the set X 1 is convex the problem (1) is convex.
Despite the fact that f 1 is linear function, the set X 2 is nonconvex and therefore the problem (2) is nonconvex and it may have many local minimizers.
In this problem both the objective function and the feasible set X 2 are nonconvex and therefore the problem (3) is nonconvex.
The problem (4)   Here we use four different versions of DGM. Their main difference is the way they compute descent directions. Therefore we will describe here only an algorithm for the computation of descent directions in DGM. The description of DGM itself can be found in the above mentioned papers.
Let f be a locally Lipschitz continuous function defined on IR n . Let S 1 = {g ∈ IR n : g = 1}, G = {e ∈ IR n : e = (e 1 , . . . , e n ), |e j | = 1, j = 1, . . . , n}, Here S 1 is the unit sphere, G is the set of vertices of the unit hypercube in IR n and P is the set of univariate positive infinitesimal functions.
Then for given x ∈ IR n and z ∈ P we define a sequence of n + 1 points as follows: Definition 1. (see [1] - [4]) The discrete gradient of the function f at the point x ∈ IR n is the vector Γ i (x, g, e, z, λ, α) = (Γ i 1 , . . . , Γ i n ) ∈ IR n , g ∈ S 1 with the following coordinates: Remark 1. One can see from the definition that n − 1 coordinates of the discrete gradient are computed differently from its i-th coordinate. They are approximations to n−1 coordinates of a certain subgradient at the point x+λg and these coordinates dependent in particular on z ∈ P . One can take z very small number and fix it for all λ > 0 or z ∈ P can be any univariate positive infinitesimal function. Discrete gradients approximate subgradients for a broad subset of nonsmooth functions [3,4].

Algorithm 1.
An algorithm for the computation of the descent direction.
Step 2. Calculate the vector w k 2 = min{ w 2 : w ∈ D k (x)}. If then stop. Otherwise go to Step 3.
Step 3. Calculate the search direction by g k+1 = − w k −1 w k .
Step 4. If then stop. Otherwise go to Step 5.
Remark 2. Algorithm 1 is a terminating [3,4]. This fact is true for any values of λ > 0. Small values of λ > 0 give approximations to subgradients of f and in this case Algorithm 1 calculates local descent directions. Large values of λ > 0 do not give approximations to subgradients anymore however they still allow one to find descent directions from x and such directions are global descent directions. Algorithm 1 is capable of finding such directions even from local minimizers.
We will consider four different versions of the discrete gradient method: DGM1L, DGM1G, DGM2L and DGM2G. In DGM1L and DGM1G discrete gradients are computed using z(λ) = λ l where l ≥ 2 whereas z(λ) = 10 −12 for any λ > 0 in DGM2L and DGM2G. Only local descent directions are computed in DGM1L and DGM2L. In DGM1G and DGM2G both local and global descent directions are computed.
The idea behind both DGM1G and DGM2G is the same. We take any starting point and applying DGM1L (or DGM2L in DGM2G) compute the first local solution. Then applying the algorithm for the computation of global descent directions we find descent direction from this local minimizer and compute a new starting point for DGM1L (or DGM2L) and so on. This continues until the algorithm for the computation of global descent directions cannot find a descent direction from a local minimizer.

4.2.
The Lipschitz-continuous global optimizer. The Lipschitz-Continuous Global Optimizer (LGO) solver suite has been successfully applied to complex, large-scale models both in research and commercial contexts for over a decade. Detailed technical descriptions and user documentation can be found in [19] and [20] and elsewhere, including the peer review [6].
LGO solver suites is based on a seamless combination of a suite of global and local scope nonlinear solvers. Currently, LGO includes the following solver options: • Branch-and-bound (adaptive partition and sampling) based global search (LGO-BB); • Adaptive global random search (with single-start) (LGO-GARS); • Adaptive multi-start global random search (LGO-MS).
• Constrained local search based on the generalized reduced gradient method (GRG) (LGO-LS). The LGO global search methodology has been described in [19], [20] and Section 2 of [21] reviews several implementations. In all three global search modes the model functions are aggregated by an exact penalty (aggregated merit) function. By contrast, in the local search phase all model functions are considered and treated individually. The global search phases are equipped also with stochastic sampling procedures that support the usage of statistical bound estimation methods. All LGO search algorithms are derivative-free: specifically, in the local search phase central differences are used to approximate gradients. The compiler-based LGO solver suite is used as an option linked to various modeling environments [21]. In its core text I/O based version, the application-specific LGO executable program (that includes a driver file, the LGO solver library and the model function file) reads an input text file that contains application-specific information (model name, variable and constraint names, variable bounds and nominal values, and constraint types) as well as a few key solver options (global solver type, precision settings, resource and time limits). Upon completing the LGO run, a summary and a detailed report file are available. As can be expected, this LGO version has the lowest demands for hardware, it also is fastest and can be directly embedded into vertical and proprietary user applications.

GAMS local and global solvers.
For an independent verification and comparison, the four test problems were run through various local and global solvers available in GAMS. Brief descriptions of the solvers are as follows (see http://www.gams.com/solvers/solvers.htm, for more details): • SNOPT is a large scale Sequential Quadratic Programming (SQP) solver developed by Philip Gill (University of California at San Diego) and Walter Murray and Michael Saunders (Stanford University). • MINOS from the Systems Optimization Laboratory at Stanford University iteratively solves subproblems with linearized constraints and an augmented Lagrangian objective function. MSMINOS is MINOS with multi-start.
• KNITRO from Ziena Optimization, Inc is a software package for finding local solutions of continuous, smooth nonlinear optimization problems, with or without constraints. • PATHNLP solves an NLP by internally constructing the Karush-Kuhn-Tucker (KKT) system of first-order optimality conditions associated with the NLP and solving this system using the PATH solver for complementarity problems. • CONOPT from ARKI Consulting and Development in Denmark is a generalized reduced gradient solver. MSCONOPT is multi-start CONOPT.

Results of numerical experiments.
It is obvious from the construction that the above formulated optimization problems can have multiple solutions, which is what we expect in real world production operations. In production operations the upper and lower bounds are generally known and for the test cases, the decision variables are in the range between 0 MMscf/day and 200 MMscf/day. The starting point for each variable is set to the mid-point value of 100 MMscf/day. However, this selection may be an infeasible starting point which it is in these test cases.
LGO is available for the Mathematica (version 5.1 was used) platform as the MathOptimizer Professional (MOP) software package. MOP auto-converts the Mathematica model code into Fortran (the compiler used was Intel Fortran version 7.1) and then the external LGO engine solves the model. Results are reported back to the calling Mathematica document (see [22]).
No tuning of search parameters was made for the algorithms and all were run in the default option mode. All versions of DGM have no user tuneable parameters which is an important feature for deployment in the operating companies. 5.1. Problem 1: maximizing gas production. The solutions obtained by LGO and DGM are shown in Table 4.
One can see from Table 4 Table 5. Solutions from various GAMS local solvers with both single start and multi-start options for Problem 1.

5.2.
Problem 2: maximizing gas production with gas quality constraints. The solutions obtained by LGO solvers and different versions of DGM are presented in Table 7.
The solutions by the various local and global solvers of GAMS are presented in Tables 8 and 9. The solvers SNOPT, MINOS, KNITRO and PATHNLP failed to locate any feasible solution. However, the multi-start MINOS managed to converge to a feasible solution but no such improvement was observed for CONOPT (see Table 8).
From the results of the various GAMS and the LGO solvers, the consensus best value for the objective function is 276.6. DGM1L and DGM1G search modes did come up with better values but these were achieved at the expense of a single constraint violation (Equation 35) in the order of 10 −2 . The maximum constraint   Table 7.
LGO and DGM solutions for Problem 2.
violation for the other three LGO and GAMS solvers is about 10 −6 . This may explain why the DGM solutions are better.
5.3. Problem 3: plant P 6 LPG maximization. Results for LGO solvers and different versions of DGM are presented in Table 10.
One can see from Table 10 that In this test case, LGO gave better solutions compared to the two DGM versions. For the first three LGO solutions the constraint violations were all less than 1.1 · 10 −7 . The constraint violation for the LGO-BB was 4.5 · 10 −6 . It is intuitive that for this test case the best strategy would be to divert as much gas to Plant P 6 as possible. This is indeed the case for the three LGO algorithms.
Results for GAMS solvers are shown in Tables 11 and 12 and these results show that Problem 3 has many solutions and different solvers given the same starting point, which is the mid-point of the lower and upper bounds. Overall, MSCONOPT (multi-start CONOPT) gave the best solution, followed by LGO-MS. MINOS, MSMINOS and SNOPT failed to locate feasible solutions.   Table 9. Solutions obtained by GAMS global solvers for Problem 2.
5.4. Problem 4: plant P 6 LPG and total condensate production maximization. The solutions are as shown in Table 13. For this problem case, the DGM gave better solutions than the LGO. The DGM solution were some 50% better.
The solutions obtained by GAMS solvers are shown in Tables 14 and 15.
One can see that Problem 4 has multiple solutions and different solvers produce different solutions starting from the the same point. Overall, DGM gave the best solutions. SNOPT, MINOS and MSMINOS failed to locate feasible solutions.
6. Discussions and conclusions. "The future is gas", said Jeroen van der Veer, the Chief Executive of the Royal Dutch Shell plc, in a speech delivered on the 23 rd World Gas Conference held in Amsterdam on June 6, 2006. It is anticipated that over the next two decades, gas either transported in pipelines or as Liquified Natural Gas (LNG) will be the major business of many major oil and gas companies. Optimization algorithms may play important role to meet the challenges of gas  Table 11. Solutions from various GAMS local solvers with both single start and multi-start options for Problem 3.
production and blending systems. Most of optimization problems in these systems are nonconvex and nonsmooth global optimization problems. It is well known that gas optimization problems, as typified by the Haverly gas pooling problem [14], are nonconvex and can have multiple solutions [11]. Production behaviors in nature may have many discontinuities and so are generally nonsmooth and nonconvex. Therefore, depending on starting points, one of the local scope solvers (such as SLP and SQP) may end up at a local minimizer. Finding the best among the many local optima may lead to substantial business advantages. On the other hand, data from gas production systems may contain noise which makes the application of gradient-based methods impossible to solve in such problems. Since the derivative-free methods do not rely explicitly on the local models of the objective and/or constraint functions such methods are suitable for solving many gas blending optimization problems. In this paper we pursue the use of derivative-free techniques such as the LGO and the DGM.
In gas production optimization with blending, it is desirable to be able to take into account the flow splitting and pressure dissipation in gathering and distribution 0.000 0.000 0.000 0.000   Table 13.
LGO and DGM Solutions for Problem 4.
networks, as shown in Figure 1. We took advantage of the capability of the DGM and LGO to solve nonsmooth problems and introduced the use of the minimum function to model chokes and flow controllers in production systems. In the forward flow from the well to the processing plant, termed the "forward-pass", the condition of flow at any gathering center, is that the operating pressure there is the minimum of the inlet pressures feeding into the gathering center. Excess pressures of the other inlets at the gathering center are dissipated through chokes and valves. In the "backward-pass", given the minimum delivery pressures at the plant, the maximum function is used to cascade the pressures and rates back to the production wells. This approach allows us to incorporate the network pressure/flow rate constraints and the blending requirements into a single optimization problem.
The above test problems model real-world situations. In this study, we are interested in the range of engineering solutions rather than just the performance of individual algorithms. Starting from a simple gas optimization under pressure constraints we gradually add blending requirements which involve fractional terms into the constraints and objective functions. The CO 2 volume fraction constraint (35)   In Problem 3, we introduced a fractional objective function, to test the capability of the DGM to address cases where blending is part of the objective function as well; in this case the maximization of LPG production at a particular plant. As shown in Table 10, only DGM2G came close to those obtained by LGO. For the GAM solvers, there is a wide spread in feasible solutions, ranging from 6.9 to 14.21, indicating that this problem has multiple solutions. As in Problem 2, when the constraint (35) is expressed as a polynomial function, DGM2G gave an objective value of 8.20 with almost no constraint violation (< 10 −8 ). In this case, DGM did not perform as good as LGO.
In Problem 4, the objective function is a sum of condensate and LPG blending at P lant 6 . In this case, intuitively a good production strategy would be to produce to the other plants as well. Here DGM did outstandingly well, with the values of the objective function in the 124.35 to 131.17, but again while violating the constraint (35) in the order of 0.007. DGM2G found some other 20 stationary points, ranging from 41 to 105. When the constraint (35) is expressed as a polynomial function, DGM2G found the solution with the objective function value 55.85 with almost no constraint violations (< 10 −8 ). The consensus value for the objective function based on the GAMS and LGO solvers is 85.35. However, depending on what is an acceptable degree of constraint violation, there can be many solutions.
Overall, compare with the LGO and the GAMS solvers, the above results are indicative that the DGM2G is a capable global optimization algorithm, less so for DGM1G. To facilitate deployment (not too many knobs for the users) both algorithms were deliberately designed with no user tuneable parameters. The above experience with fractional objective functions and constraints suggests that the DGM is not as robust as LGO when dealing fractional functions.
As the dimension in all problems is small we do not give any detail analysis of other performance indicators.
A more detailed comparative study of the DGM using a large number of real world problems is the subject of our future research. Results with the above test problems demonstrate that LGO and DGM have their own strength; together they can be very useful in production optimization. For engineering applications it is better to have a broad spectrum of global optimization solvers, some methods are more appropriate than others for some special optimization opportunities. The LGO has now been incorporated into the Shell Hydrocarbon Field Planning Tool (HFPT). Comparative studies of HFPT/LGO against off-the-shelf 3 rd party integrated production system software have shown that LGO can come up with better alternative solutions. In the case of a optimal gas lift allocation problem to maximize production [15], the HFPT/LGO combination came up with better gas lift allocation and usage (see Figure 4), which has resulted in better oil production (see Figure 5). At the current oil price of some US$70/barrel, the difference observed may have significant impact on the field development strategy and revenue.