Stochastic Programming Approach for the Planning of Offshore Oil or

3078

Ind. Eng. Chem. Res. 2009, 48, 3078–3097

Stochastic Programming Approach for the Planning of Offshore Oil or Gas Field Infrastructure under Decision-Dependent Uncertainty Bora Tarhan,† Ignacio E. Grossmann,*,† and Vikas Goel‡ Department of Chemical Engineering, Carnegie Mellon UniVersity, Pittsburgh, PennsylVania 15213, and ExxonMobil Upstream Research Company, Houston, Texas 77098

The planning of offshore oil or gas field infrastructure under uncertainty is addressed in this article. The main uncertainties considered are in the initial maximum oil or gas flowrate, the recoverable oil or gas volume, and the water breakthrough time of the reservoir, which are represented by discrete distributions. Furthermore, it is assumed that these uncertainties are not immediately realized, but are gradually revealed as a function of design and operation decisions. To account for these decision-dependent uncertainties, we propose a multistage stochastic programming model that captures the complex economic objectives and nonlinear reservoir behavior and simultaneously optimizes the investment and operating decisions over the entire planning horizon. The proposed solution algorithm relies on a duality-based branch-and-bound method involving subproblems as nonconvex mixed-integer nonlinear programs. Several examples involving nonlinear reservoir models are presented to illustrate the application of the proposed method. 1. Introduction Oil and gas field exploration and production operations consist of four major steps: exploration, appraisal, development, and production. In each step, many decisions have to be made that affect the overall performance of the operation. In the beginning of the exploration phase, many uncertainties exist, and depending on the decisions made at each step, the uncertainty reduces gradually. Unfortunately, many crucial decisions related to the planning of offshore oil or gas field infrastructure have to be made in the early steps when the uncertainty level is high. The quality of these decisions affects the overall profitability of the operation. The objective of this work is to develop and solve a model to optimize the decisions about planning of offshore oil or gas field infrastructure under gradual uncertainty resolution. For the sake of simplicity, we refer throughout the article to “oil and/or gas” as “oil” for short, but the presented model and solution algorithm applies to both of them equally. An oil field consists of several reservoirs, each of which contains a number of potential wells (Figure 1). These potential wells can be drilled and exploited for oil using different facilities. Previous work on the planning of oil field infrastructure can be classified into strategic, tactical, and operational approaches according to their time frames; whether uncertainty is considered; and if so, the way uncertainty is handled, i.e., simulationor optimization-based solution methods. Haugland et al.1 proposed simultaneous mixed-integer linear programming (MILP) models for oil field design and production planning. They demonstrated that their model can be solved by commercial solvers for only small instances because of the computational burden, thus necessitating the use of customized algorithms that take advantage of the special problem structure. Iyer et al.2 addressed the optimal planning and scheduling of offshore oil field infrastructure investment and operations. They proposed a multiperiod mixed-integer linear programming model that optimizes the planning and scheduling of investment and operating decisions. These decisions are the selection of * To whom correspondence should be addressed. E-mail: grossmann@ cmu.edu. † Carnegie Mellon University. ‡ ExxonMobil Upstream Research Company.

reservoirs to develop, selection of well sites and the production rates from wells of each time period. The model incorporates the nonlinear reservoir performance through piecewise linear approximations. van den Heever and Grossmann3 proposed a multiperiod mixed-integer nonlinear programming (MINLP) model for oil field infrastructure planning, for which they developed a bilevel decomposition method. In contrast to Iyer and Grossmann,4 they explicitly incorporated a nonlinear reservoir model into their formulation. In both articles, one of the major assumptions was that there was no uncertainty in the parameters. van den Heever et al.5 and van den Heever and Grossmann6 extended the work to handle complex economic objectives such as royalties, tariffs, and taxes. They incorporated these complexities into their model through disjunctions and converted the resulting disjunctive problem into a convex MINLP by obtaining the convex hull of each disjunction and using convex underestimators for nonconvexities in the constraints. Lin and Floudas7 investigated the long-term planning problem for gas field development. They optimized the investment and operating decisions using an MINLP model that incorporated nonlinear reservoir models and complex economic calculations using a continuous time formulation. They also proposed a twostage algorithm for solving the MINLP model in which they solved a simplified model in the first stage and the detailed model in the second stage. Carvalho et al.8 presented reformulations of the MILP model by Tsarbopoulou9 and a bilevel decomposition algorithm for solving the larger instances of the problem. They used design cuts to avoid subsets and supersets as proposed by Iyer and Grossmann4 and reported a significant improvement in solution time through the use of decomposition and design cuts in the algorithm. Furthermore, Carvalho et al.10 extended their work to handle multiple reservoirs. Behrenbruch11 presented a qualitative discussion on the development planning process, feasibility studies, cost estimating, economic evaluation, and risk analysis for offshore oil field. Guidelines for evaluating the feasibility of offshore petroleum projects in terms of appraisal programs, field development, and facility options were provided. This work emphasized the multidisciplinary nature of the problem, including the need to

10.1021/ie8013549 CCC: $40.75  2009 American Chemical Society Published on Web 02/23/2009

Ind. Eng. Chem. Res., Vol. 48, No. 6, 2009 3079

consider the correct geological model and to incorporate flexibility into the decision process. Jonsbraten12 presented an MILP model for optimal development of an oil field under oil price uncertainty. The model predicts decisions related to both design and operation while maximizing the expected net present value of the project. The author used a progressive hedging algorithm to solve the problem. Lund13 considered the value of flexibility in offshore petroleum projects through a stochastic dynamic programming model. The author assumed uncertainty in the oil price and recoverable oil volume and illustrated the change in the objective value with respect to the changes in the flexibility of possible decisions. Similarly, Begg et al.14 extended the value of the information concept to assess the value of flexibility in managing uncertainty in oil and gas investments. They used decision trees or asset models to find the optimum set of decisions and presented three examples to show the relationships among uncertainty, flexibility, value of information, and value of flexibility. Cullick et al.15 proposed the integration of a global search optimization algorithm, a finite-difference reservoir simulation, and economics. In the optimization algorithm, new decision variables were generated using metaheuristics, and uncertainties were handled through simulations for fixed design variables. They presented examples including multiple oil fields with uncertainties in the reservoir volume, fluid quality, deliverability, and costs. Aseeri et al.16 addressed financial risk management in the planning and scheduling of offshore oil infrastructures. They introduced uncertainty, risk management, and budgeting constraints into the model proposed by Iyer and Grossmann.4 A sampling average algorithm was used to overcome the numerical difficulties, and the results were compared with optimum results found using upper-bound risk curves. Ligero et al.17 proposed a new methodology to assess the value of information based on the simulation of several reservoir models. They showed the complexity of evaluating the value of information in the appraisal and development phases. Goel and Grossmann18 considered the gas field development problem under uncertainty in the size and quality of reserves where decisions on the timing of field drilling were assumed to yield an immediate resolution of the uncertainty. Linear reservoir models, which provide a reasonable approximation for gas fields, were used. In their solution strategy, the authors used a relaxed problem to predict upper bounds and solve multistage stochastic programs for a fixed scenario tree for finding lower bounds. Goel et al.19 later proposed a branch-and-bound algorithm for solving the corresponding disjunctive/mixed-integer programming model where the lower bounds are generated by Lagrangean duality. Ulstein et al.20 presented a model for the tactical planning of Norwegian petroleum production. The model maximized the net income before taxes from the production and sale of petroleum products. Different cases with demand variations, quality constraints, and system breakdowns were considered. The model was solved for different scenarios, and solutions were compared with the base-case scenario. The benefit of the model is its ability to identify feasible ways to satisfy the demand for varying network configurations. Ozdogan and Horne21 emphasized the value of time-dependent information for achieving better decisions in terms of reduced uncertainty and expected net present value. They used history-matched realizations and pseudohistory concepts to simulate and optimize the well placement decisions. The

application of this approach is rather limited because of the computational load of history-matching steps. This limits the application to fields with small numbers of wells. In addition to extensive mathematical programming and deterministic global optimization algorithms, many works have used a combination of reservoir modeling, economics, and decision making under uncertainty through simulation/optimization frameworks as a solution method.14,15,22,23 Another important aspect of this work is decision-dependent uncertainty. Most previous studies including uncertainty have considered exogenous uncertainty where stochastic processes are independent of decisions (e.g., demands), whereas problems in which stochastic processes are affected by decisions are said to have endogenous uncertainty. Decisions can affect stochastic processes in two different ways. One is that decisions can change the probability functions. Studies in which such an effect was considered include those by Ahmed,24 Vishwanath et al.,25 and Held et al.26 Ahmed24 presented examples on network design, server selection, and facility location problems with decision-dependent uncertainties, showing that these programs can be reformulated as MILP problems and solved by LP-based branch-and-bound algorithms. Vishwanath et al.25 addressed a network problem having endogenous uncertainty in survival distributions. The problem was a two-stage stochastic program in which first-period investment decisions were made for changing the survival probability distribution of arcs after a disaster. The aim was to find the investments that would minimize the expected shortest path from source to destination after a disaster. Held et al.26 worked on another instance in which the problem included endogenous uncertainty in the structure of a network. In each stage of this problem, an operator tried to find the shortest path from a source to a destination after some of the nodes in the network had become blocked. The aim was to maximize the probability of stopping the flow of goods or information in the network. Another way the decisions can impact the stochastic process is that they can affect the resolution of uncertainty or the time uncertainty resolves. Such decisions have been considered in Pflug,27 Jonsbraten et al.,28 and Goel and Grossmann.18 Pflug27 addressed the problems in the context of discrete-event dynamic systems in which the underlying stochastic process depended on the optimization decisions. Jonsbraten et al.28 addressed endogenous uncertainty where project decisions reduced the uncertainty and proposed an implicit enumeration algorithm in which decisions that affect the uncertain parameter values are made at the first stage. Tarhan and Grossmann29 considered the synthesis of process networks. They assumed that there were uncertainties in the yields of the processes and that these uncertainties resolved gradually over time depending on the investment and operating decisions. The problem was modeled as a disjunctive program and then converted to an MILP using big-M transformations. To solve the model, the authors used a duality-based branchand-bound procedure in which each scenario subproblem was solved independently during the subgradient optimization for calculating upper bounds and heuristics were used for generating feasible solutions. In this article, we address the planning of offshore oil field infrastructure involving endogenous uncertainty in the initial maximum oil flowrate, recoverable oil volume, and water breakthrough time of the reservoir, where decisions affect the resolution of these uncertainties. We extend the work of Goel and Grossmann18,30 and Goel et al.19 but with three major differences. The first is that the model focuses on a single field consisting of several reservoirs rather than multiple fields.

3080 Ind. Eng. Chem. Res., Vol. 48, No. 6, 2009

Figure 1. Configuration of fields, reservoirs, and wells.

However, more detailed decisions such as the number, type, and construction decisions for infrastructure are considered. The second difference is that nonlinear, rather than linear, reservoir models are considered. Nonlinear reservoir models are important because Goel and Grossmann18 focused on gas fields, for which linear models are an adequate approximation, but for oil fields, nonlinear reservoir models are required. Finally, the third difference is that the resolution of uncertainty is gradual over time instead of being resolved immediately. This article is organized as follows: In section 2, we present the statement of the problem under consideration, and then in section 3, we discuss the representation of the planning horizon and uncertainty. In sections 4 and 5, respectively, the optimization model and the proposed solution approach are explained in detail. Section 6 discusses the implementation details and the results found by the proposed solution approach.

and the installation schedule of these facilities, as well as selection of the types of wells and the drilling schedule for the wells. Operating decisions are amount of oil production for each time period given the limitations of the reservoirs. The goal is to capture the complex economic tradeoffs that arise from the investment and operating decisions in order to maximize the expected net present value of the project. Four major assumptions about the state of the reservoir are included in the problem. First, there is no free gas in the reservoir. Second, there is a strong aqueous support that creates sufficient pressure in the reservoir. As an extension of the first two assumptions, we assume there is no need for enhanced recovery, i.e., no need for injection of gas or water into the reservoir. Finally, all wells in one reservoir are identical, which leads to the same maximum oil flowrate for all wells in the reservoir given the cumulative oil extracted from that reservoir. These assumptions are made for the sake of simplicity, and both the model and the solution algorithm are flexible enough to incorporate more complex reservoir models. Figure 5 presents the nonlinear reservoir model. It shows the flowrate of oil and water from a single well versus the cumulative recovered oil. The maximum oil flowrate can be represented as a nonlinear function of the cumulative amount of oil recovered from the reservoir (eq 1) as shown in Figure 5. During oil recovery from a reservoir, the liquid coming from the well contains not only oil but also water, and the relative rates of these liquids can be characterized by the water-to-oil ratio. This ratio is approximated using a nonlinear function of the cumulative oil recovered and the recoverable oil volume of the reservoir (eq 2). The water flowrate can be calculated by multiplying the oil flowrate by the water-to-oil ratio as in eq 3.

2. Problem Statement In this article, the problem that we consider is the design and planning of an offshore oil field infrastructure (Figure 2) over a specified planning horizon. Specifically, we consider a field consisting of several reservoirs where a number of wells can be drilled and exploited for oil in every reservoir during the planning horizon. The field infrastructure can be composed of floating production storage and offloading (FPSO) (Figure 3) and/or tension leg platform (TLP) (Figure 4) facilities. An FPSO facility can be either a small FPSO, converted from a retiring oil tanker, or a large FPSO, a newly constructed grassroots facility. An FPSO facility can process, store, and offload the processed oil to other tankers. Processing means separating the oil and water that comes out of the well. Unlike an FPSO facility, a TLP facility cannot process oil; it has only drilling and oil recovering capabilities. TLP and FPSO facilities can be connected to each other through risers and subsea pipelines. Oil recovered from TLP facilities is pumped to FPSO facilities through these pipes. Each facility has a construction cost and a lead time between the construction decision and the actual startup. There are two options for drilling wells: Each well can be drilled as either a subsea or a TLP well. Drilling ships are used to drill subsea wells, so there is no need to have a facility present to drill a subsea well. Unlike subsea wells, each TLP well must be drilled by a TLP facility. For economic reasons, each type of well is drilled in groups consisting of a fixed number of wells. To recover oil from a well, the well must be connected to a facility. A subsea well must be connected to an FPSO facility, whereas a TLP well must be connected to a TLP facility. The problem involves making investment and operating decisions over the planning horizon. Investment decisions involve selection of the number, type, and capacity of facilities

∀t, ∀s, ∀r

cum,r,s oilr,s t e f(oilt-1 )

(

r,s worr,s t ) R

cum,r,s oilt-1

RECr,s

)

(1)

βr

(),r,s water(),r,s ) worr,s t t oilt

∀t, ∀s, ∀r

(2)

∀t, ∀s, ∀r

(3)

Uncertain parameters in the system are the initial maximum oil flowrate, the recoverable oil volume, and the water breakthrough time of the reservoir. Because we assumed that all wells within one reservoir are identical, we consider the uncertainty in the initial maximum oil flowrate for only a single well in the reservoir. The three uncertain parameters affect the nonlinear reservoir behavior in Figure 5 in different ways. The initial maximum oil flowrate affects the point where the maximum oil flowrate starts on the vertical axis in Figure 5. The amount of recoverable oil in the reservoir affects where the oil flowrate becomes zero, i.e., the point where the curve that represents the oil flowrate intersects the horizontal axis. The third uncertain parameter, the water breakthrough time, is the time elapsed until the water flowrate exceeds a prespecified value. We assume that the water breakthrough time is correlated with the scalar R in the nonlinear function representing the water-to-oil ratio (eq 2). Therefore, the uncertainty in the water breakthrough time can be represented by the uncertainty in the parameter R. If R has a high value, the water-to-oil ratio increases rapidly, as does the water flowrate (eq 3), making the water breakthrough time short. If R has a low value, the waterto-oil ratio and water flowrate increase in a moderate way, making the water breakthrough time longer. To consider the gradual resolution of uncertainties, we make the following assumptions:


Figure 2. Typical oil field infrastructure.

(a) We assume that an appraisal program is available to resolve the uncertainty in initial maximum oil flowrate. The appraisal program consists of drilling a number, N1, of wells. This appraisal program not only gives the actual value for the initial maximum oil flowrate, but also provides the posterior probabilities for recoverable oil from the reservoir depending on the outcome of the initial maximum oil flowrate. Therefore, at the end of the appraisal program, we resolve the uncertainty only in the initial maximum oil flowrate, but not in the volume of recoverable oil in the reservoir or the water breakthrough time. (b) We assume that the uncertainty in the recoverable oil volume can be resolved in two ways. One is to drill a total of more than N2 wells (N2 > N1), and the other is to produce from that reservoir for a duration of N3 years. Either of these actions has the same effect of resolving the uncertainty in the volume of recoverable oil. (c) The uncertainty in the water breakthrough time is resolved independently of the drilling decisions; it is affected only by the production. Specifically, production from the reservoir for N4 years resolves the uncertainty in the water breakthrough time. Therefore, it is possible to resolve the uncertainties in the initial maximum oil flowrate and recoverable oil by drilling wells before production from the reservoir, but there must be production to resolve the uncertainty in the water breakthrough time. Note that these assumptions on uncertainty resolution are completely flexible and can be modified depending on one’s assessments of how the uncertainties can be resolved.

Figure 3. Typical FPSO facility. (Reproduced from http://offalnews. blogspot.com/2007/03/producing-projects-terra-nova.html.) (accessed June 9, 2008).

Figure 4. Typical TLP facility. (Reproduced from http://en.wikipedia.org/ wiki/Image:Mars_Tension-leg_Platform.jpg.) (accessed June 9, 2008).

Figure 5. Nonlinear reservoir model.

3. Representation of Planning Horizon and Uncertainty In this work, the planning horizon is discretized into time periods and the probability distributions of uncertain parameters are discrete because these specifications allow us to represent the stochastic process by scenario trees. Figure 6 is a standard representation of a scenario tree having one uncertain parameter with two discrete values in two time periods, which leads to four scenarios. Uncertain parameters ξ1 and ξ2 are revealed at the end of the first and second time periods, respectively, generating four equally probable scenarios, given that the parameters have a uniform distribution. In a standard scenario tree, each node for time period t represents a possible state for that time period. Each arc represents the possible transition from one state in time period t to another state in time period t + 1. A path from the root node to a leaf node represents a scenario. Thus, a scenario is a combination of possible uncertain parameters in each of the time periods. The set of time periods that have the same amount of information defines a stage. Problems that have more than two stages are called multistage stochastic programs.

Figure 6. Scenario tree with uncertain parameters ξ1 and ξ2.

Figure 7 is an alternative representation of the scenario tree in Figure 6, proposed by Ruszczynski.31 In this representation, each scenario is represented by a set of unique nodes. The horizontal lines connecting nodes in time period t mean that these nodes are indistinguishable and have the same amount of information in that time period. The horizontal lines reduce the tree in Figure 7 to that in Figure 6. For modeling of the problem in this work, the scenario trees will be considered according to the representation given by Ruszczynski31 (see Figure 7) in order to explicitly handle the nonanticipativity constraints that give rise to different tree structures (see Goel and Grossmann30).


revst ) C1,tprodst

∀t, ∀s

(6)

The total cost at every time period t is the sum of capital and operating expenditures ∀t, ∀s

costst ) capexst + opexst Figure 7. Alternative scenario tree with uncertain parameters ξ1 and ξ2.

(7)

Capital expenditures at every time period t come from drilling a group r,s of TLP wells (dtwr,s t ) and/or subsea wells (dswt ) and building TLP (btst ), small FPSO (bsfst ), and large FPSO (blfst ) facilities during that time period capexst )

∑C

∑C

+

r,s 2,tdtwt

r,s 3,tdswt

r∈R

+ C4,tbtst +

r∈R

C5,tbsfts + C6,tblfts

∀t, ∀s (8)

The operating expenditure at every time period t is a linear function of the amount of oil produced ∀t, ∀s

opexst ) C7,tprodtotal,s t

Because of limitations on the maximum number of drilling rigs on a TLP facility, the number of TLP wells to be drilled (dtwr,s t ) during time period t by each TLP facility is limited by some upper bound (Utw)

Figure 8. Graphical representation of the aggregate infrastructure.

4. Mathematical Model In this section, we present the optimization model for the planning of offshore oil field infrastructure under uncertainty. This model optimizes investment decisions such as the numbers of wells to drill and facilities to build and operational decisions such as the oil production rate from the reservoirs to maximize the expected net present value. As explained in the previous sections, the interaction between decisions creates many tradeoffs and makes the decision-making process highly complex. The aggregate infrastructure is shown graphically in Figure 8. The variables in the model can be classified as decision, state, and recourse variables. Decision variables are related to decisions that are made at the beginning of each time period t and scenario s (e.g., numbers of wells to drill and facilities to build). Recourse variables are related to decisions made after partial or full resolution of uncertainty (e.g., oil production from reservoirs). State variables are the variables that are calculated automatically when decision and recourse variables at previous time periods are selected (e.g., water-to-oil ratio, number of subsea/TLP wells available for production). The sequence of events is as follows: Decision variables are implemented at the beginning of time period t. This is followed by the resolution of uncertainty. Recourse decisions are made at the end of period t. The sets and indices, continuous variables, integer variables, binary/logic variables, and parameters used in this work are defined in the Nomenclature section. Given these sets, variables, and parameters, the model (P) is as follows. The objective is to maximize the expected net present value of the project given by max enpv )

∑ P npv s

s

(4)

s∈S

The net present value of each scenario s can be calculated using revenues and costs for that scenario as npvs )

∑

Dt(revst - costst )

∀s

(9)

(5)

t∈T

The revenue at every time period t can be calculated using the total amount of oil produced and the sale price of oil at each time period

∑ dtw

r,s t

∀t, ∀s

e Utwntst

(10)

r∈R

Similarly, the number of subsea wells to be drilled (dswr,s t ) during time period t is limited by some upper bound (Usw)

∑ dsw

r,s t

e Usw

∀t, ∀s

(11)

r∈R

A group of TLP wells (Gtw) in reservoir r can be drilled and connected to TLP facilities during time period t (ntwr,s t > 0) only if at least one TLP facility has been built before that time period (ntst > 0). Also, there is a limit on the maximum number of wells to be drilled and connected per TLP facility (Nwmax)

∑ ntw

r,s t

e Nwmaxntst

∀t, ∀s

(12)

r∈R

Oil and water can be extracted from a well in reservoir r and pumped to TLP facilities during time period t only if there is a TLP well available during time period t ∀t, ∀s, ∀r

oiltlp,r,s + watertlp,r,s e M1ntwr,s t t t

(13)

There is an upper bound on the total number of TLP and subsea wells that can be drilled in each reservoir r,s w,r ntwr,s t + nswt e U

∀t, ∀s, ∀r

(14)

Oil can be recovered from reservoir r during time period t if there is a subsea well available for production during time period t oilfpso,r,s + waterfpso,r,s e M2nswr,s t t t

∀t, ∀s, ∀r

(15)

The big-M values M1 and M2 in eqs 13 and 15 can be taken as the maximum value of the sum of the maximum oil flowrate and the maximum water flowrate. These flowrates are well-defined curves in Figure 5. The rate of oil recovered from a well in reservoir r during time period t is bounded above by a linear or nonlinear function


of the cumulative amount of oil produced from the same reservoir up to time period t. Also, the oil flowrate must be 0 if there is no well ready for production (bwr,s t ) 0). oilr,s t

e

cum,r,s 2 γr,s 1 (oilt-1 )

+

+

cum,r,s γr,s 2 oilt-1

∀t, ∀s, ∀r (16)

r,s γr,s 3 bwt

Oil can be recovered from reservoir r during time period t if there is a well ready for production during time period t ∀t, ∀s, ∀r

r,s r,s bwr,s t e ntwt + nswt

(17)

The amount of oil pumped to TLP and FPSO facilities during time period t can be calculated by the numbers of available r,s TLP (ntwr,s t ) and subsea (nswt ) wells for production and the oil flowrate per well from the reservoir r,s oiltlp,r,s ) oilr,s t t ntwt r,s oilfpso,r,s ) oilr,s t t nswt

)

) prodtotal,s t

The number of wells and facilities available to operate during time period t (e.g., ntst ) can be calculated by summing the available wells and facilities during the previous time period t - 1 (e.g., ntst-1) and those that started to be built earlier (e.g., btst-τ1) and are expected to be ready at time period t

(

REC

)

s s + bsft-τ nsfst ) nsft-1 2

∀t, ∀s

(29)

∀t, ∀s

(30)

(19)

r,s r,s ntwr,s t ) ntwt-1 + dtwt-τ4

∀t, ∀s, ∀r

(31)

r,s r,s nswr,s t ) nswt-1 + dswt-τ5

∀t, ∀s, ∀r

(32)

∀t, ∀s, ∀r ∀t, ∀s, ∀r

∀t, ∀s, ∀r

(20)

There is an upper bound on the number of each type of facility (TLP, small FPSO, large FPSO) that can be built during the entire planning horizon

(21)

(22)

∀t, ∀s, ∀r (23)

cum,r,s oilcum,r,s ) oilt-1 + δt(oiltlp,r,s + oilfpso,r,s ) t t t

) and TLP The total amount of oil recovered by FPSO (oilfpso,r,s t ) facilities must be smaller than the total oil processing (oiltlp,r,s t capacity of the small FPSO (CAPsfoil) and large FPSO (CAPlfoil) facilities CAPsfoilnsfts + CAPlfoilnlfts g

fpso,r,s t

+ oiltlp,r,s ) t

∀t, ∀s (24)

r∈R

ntst e Ut

+ waterfpso,r,s + oiltlp,r,s + watertlp,r,s ) t t t

∀t, ∀s

r∈R

(25) The amount of production from reservoir r during time period t is calculated by summing the oil extracted by FPSO (oilfpso,r,s ) and t ) TLP facilities (oiltlp,r,s t fpso,r,s + oiltlp,r,s ) prodr,s t ) δt(oilt t

(33)

∀t, ∀s

(34)

nlfts e Ulf

∀t, ∀s

(35)

Equations 36-40 relate the state variables and the binary/logic variables (wr,s (.),t) that will be used to determine the indistinguishable scenario pairs at every time period. The total number of wells is r,s is 1/true. less than N1 if and only if binary/logic variable w1,t r,s is 1/true. Similarly, the total is less than N2 if and only if w2,t r,s r,s wr,s 1,t S (ntwt + nswt eN1 - 1)

∀t, ∀s, ∀r

(36)

r,s r,s wr,s 2,t S (ntwt + nswt eN2 - 1)

∀t, ∀s, ∀r

(37)

Production from reservoir r has occurred for less than N3 years if r,s is 1/true, and it has occurred for less than N4 years and only if w3,t r,s is 1/true if and only if w4,t

( (∑ t-1

wr,s 3,t S

∑ bp

∀t, ∀s, ∀r

(38)

bpr,s τ e N4 - 1

∀t, ∀s, ∀r

(39)

t-1

wr,s 4,t S

) )

e N3 - 1

r,s τ

τ)1

where bpr,s t is given by

CAPsfliqnsfts + CAPlfliqnlfts g fpso,r,s t

∀t, ∀s

nsfts e Usf

τ)1

The total liquid (oil and water) flowrate coming into the FPSO facilities, including the liquid from TLP facilities, must be smaller than the total liquid capacity of small FPSO (CAPsfliq) and large FPSO (CAPlfliq) facilities

∑ (oil

(28)

∀t, ∀s, ∀r

The cumulative volume of oil extracted from reservoir r up to the end of time period t is calculated as

∑ (oil

∀t, ∀s

s s ntst ) ntt-1 + btt-τ 1

The water-to-oil ratio in reservoir r is a nonlinear function of the cumulative oil production from the reservoir and the recoverable oil )

(27)

r

s s + blft-τ nlfst ) nlft-1 3

fpso,r,s waterfpso,r,s ) worr,s t t oilt

worr,s t

∀t, ∀s

r,s t

(18)

tlp,r,s worr,s t oilt

cum,r,s βr oilt-1 r,s R r,s

∑ prod

∀t, ∀s, ∀r

The water flowrate from reservoir r to TLP or FPSO facilities in time period t can be calculated using the water-to-oil ratio and the oil flowrate for wells in that reservoir watertlp,r,s t

The total amount of production during time period t is calculated by summing the production over all reservoirs

∀t, ∀s

(26)

r,s bpr,s t S (prodt gε)

∀t, ∀s, ∀r

(40)

Equations 41-47 are used to find the scenarios (s,s′) that are indistinguishable in time period t depending on the decisions made up to that time period. Actually, eqs 41-47 are a shortened version of the original equations, which are defined for both (s,s′) and (s′,s). The reduction was carried out using the proposition.30 Two scenarios (s,s′) that differ only in initial productivities (i.e., element of the set M1) will be indistinguishable at time period t if and only

3084 Ind. Eng. Chem. Res., Vol. 48, No. 6, 2009 Table 1. Sets That Differ Only in the Specified Parameters M1 initial maximum oil flowrate recoverable oil volume water breakthrough time

M2

×

×

M3

×

M4

M5

× ×

× ×

M6

M7

× ×

× × ×

able until there is production for N4 periods. Given restrictions on production such as the availability of wells and facilities and delays for building infrastructure, it can be shown that these scenarios must be indistinguishable for a certain number of years. Such nonanticipativity constraints are also included in the initial nonanticipativity constraints.

if for each reservoir that distinguishes those scenarios [i.e., element r,s is true (eq 41). Equations 42-47 can be interpreted of D(s,s′)] w1,t similarly. Note that each possible pair of scenarios (s,s′) will be an element of exactly one of the sets Mk, k ) 1,..., 7, that refer to the sets of scenarios that differ in the initial maximum oil flowrate, recoverable oil, and water breakthrough time of the reservoir (see Table 1). zs,s′ S t S zs,s′ t

∧ wr,s 1,t

r∈D(s,s′)

∀t, t g 2, ∀(s, s′) ∈ M1, s < s′ (41) ∀t, t g 2, ∀(s, s′) ∈ M2, s < s′

r,s ∧ wr,s 2,t ∧ w3,t

r∈D(s,s′)

(42) S zs,s′ t S zs,s′ t

zs,s′ S t

zs,s′ S t

zs,s′ S t

∧ wr,s 4,t

r∈D(s,s′)

∀t, t g 2, ∀(s, s′) ∈ M3, s < s′ (43)

r,s r,s ∧ wr,s 1,t ∧ w2,t ∧ w3,t

r∈D(s,s′)

∀t, t g 2, ∀(s, s′) ∈ M4, s < s′ (44)

r,s ∧ wr,s 1,t ∧ w4,t

r∈D(s,s′)

∀t, t g 2, ∀(s, s′) ∈ M5, s < s′(45)

r,s r,s ∧ wr,s 2,t ∧ w3,t ∧ w4,t

r∈D(s,s′)

∀t, t g 2, ∀(s, s′) ∈ M6, s < s′ (46)

r,s r,s r,s ∧ wr,s 1,t ∧ w2,t ∧ w3,t ∧ w4,t

r∈D(s,s′)

∀t, t g 2, ∀(s, s′) ∈ M7, s < s′ (47)

Note that it is sufficient to consider only eqs 41-43 that run over M1, M2, and M3 instead of constraints over all subsets of M. Assume scenario pairs (s1,s2) in M1 and (s1,s3) in M2 are indistinguishable. Using the relation in Table 1, it can be inferred that scenario pair (s2,s3), which is in M4, will be indistinguishable as well. Similarly, all scenario pairs in subsets M4-M7 can be inferred by the pairs in sets M1-M3, making it possible to reduce the constraints in eqs 41-47 to only those in eqs 41-43. The disjunctive constraint shows the decisions that are the same for time period t if s and s′ are indistinguishable at the beginning of time period t. By the same argument, eq 48 is also run over subsets M1, M2, and M3.

[ ] zs,s′ t

btst ) bts′t

bsfts ) bsfs′t blfts ) blfs′t

r,s′ dtwr,s t ) dtwt

∨ [¬zs,s′ t ]

∀r

r,s′ ∀r dswr,s t ) dswt ∀t g 2, ∀(s, s′) ∈ M1 ∪ M2 ∪ M3, ∀(s, s′, t) ∈ NC (48)

Equations 49-53 represent the initial nonanticipativity constraints. Initial nonanticipativity constraints include the firstperiod nonanticipativity constraints and some subsequent-time nonanticipativity constraints that appear as a result of the gradual resolution of uncertainty. For instance, the scenario pairs that differ only in the water breakthrough time will be indistinguish-

∀(s, s′, t) ∈ NI, ∀(s, s′, t) ∈ NI

btst ) bts′t

(49)

bsfts ) bsft s′

∀(s, s′, t) ∈ NI

(50)

blfts ) blfts′

∀(s, s′, t) ∈ NI

(51)

r,s′ dtwr,s t ) dtwt

∀(s, s′, t) ∈ NI, ∀r

(52)

r,s′ dswr,s t ) dswt

∀(s, s′, t) ∈ NI, ∀r

(53)

The model given by eqs 4-53 corresponds to a mixed-integer nonlinear/disjunctive programming model. The most direct way of solving this problem is to convert it into a mixed-integer nonlinear model (see Appendix B) by transforming the logic propositions into linear inequalities and disjunctions into mixedinteger constraints (e.g., see Raman and Grossmann32). Also, note that the model has been generated using the scenario tree approach proposed by Ruszczynski31 so that each scenario is represented by a set of unique variables where nonanticipativity constraints relate these variables for different scenarios. Also, most of the constraints are linear, except for a few relating the water-to-oil ratio to the cumulative oil production and the water flowrate, oil flowrate, and water-to-oil ratio. These constraints are nonconvex, which necessitates the use of a global optimization algorithm as explained in the next section. 5. Solution Strategy The size of model P increases exponentially with the number of scenarios and time periods, making this nonconvex MINLP model extremely difficult to solve in full space for real-size problems using commercial solvers. The model P is composed of subproblems connected through initial (eqs 49-53) and conditional (eqs 41-48) nonanticipativity constraints. Each subproblem is a nonconvex mixed-integer nonlinear program and can be solved independently when these nonanticipativity constraints are relaxed. Therefore, we propose a duality-based branch-and-bound algorithm that takes advantage of the special problem structure. In the following subsections, the upper and lower bounding procedures used at each node of the branch-and-bound tree, the branching scheme, and the proposed solution algorithm are presented. 5.1. Upper Bounding Procedure. In the proposed algorithm, the upper bound at the root node of the branch-and-bound tree is found by optimizing model PLR 0 in which the logic constraints (eqs 41-47) and disjunction (eq 48) have been removed and the initial nonanticipativity constraints (eqs 49-53) are dualized as follows (PLR 0 )

s,s′ s,s′ s,s′ s,s′ s,s′ φLR 0 (λbt,t, λbsf,t, λblf,t, λdtw,r,t, λdsw,r,t) )

max

∑ P npv s

s∈S

∑λ

s

+

∑

s,s′ [λbt,t (btst - bts′t ) +

(s,s′,t)∈NI

s,s′ s,s′ λbsf,t (bsfst - bsfs′t ) + λblf,t (blfst - blfs′t ) +

s,s′ r,s dtw,r,t(dtwt

- dtwr,s′ t ) +

r∈R

subject to eqs 5-40.

∑λ

s,s′ r,s dsw,r,t(dswt

r∈R

- dswr,s′ t )] (54)

Ind. Eng. Chem. Res., Vol. 48, No. 6, 2009 3085 s,s′ λbt,t ,

s,s′ λbsf,t ,

s,s′ λblf,t ,

s,s′ λdtw,r,t ,

s,s′ λdsw,r,t

The parameters and represent the Lagrange multipliers corresponding to constraints 49-53, respectively. To find the tightest upper bound generated by model P0LR at the root node, we consider the Lagrangean dual problem P0LD that minimizes the model P0LR in the space of multipliers (PLD 0 )

φn ) max

∑ P npv s

s

(56)

s∈S

subject to eqs 4-40 and ∀(s, s′, t) ∈ NI ∪ NnC

btst ) bts′t

(57)

bsfst ) bsfs′t

∀(s, s′, t) ∈ NI ∪ NnC

(58)

blfst ) blfs′t

∀(s, s′, t) ∈ NI ∪ NnC

(59)

r,s′ dgtwr,s t ) dgtwt

∀(s, s′, t) ∈ NI ∪ NnC, ∀r

(60)

r,s′ dgswr,s t ) dgswt

∀(s, s′, t) ∈ NI ∪ NnC, ∀r

(61)

The upper bound at node n is generated by optimizing model PnLR in which the nonanticipativity constraints (eqs 57-61) are dualized as follows φLR n ) max

∑ P npv + ∑ s

s∈S

s

s s′ [λs,s′ bt,t(btt - btt ) +

(s,s′,t)∈NI∪NCn

s,s′ s,s′ (bsfst - bsfs′t ) + λblf,t (blfst - blfs′t )+ λbsf,t

s,s′ r,s dtw,r,t(dtwt

- dtwr,s′ t ) +

r∈R

∑λ

s,s′ r,s dsw,r,t(dswt

- dswr,s′ t )]

r∈R

(62) subject to eqs 4-40. As for the root node, to find the tightest upper bound, we again consider the Lagrangean dual problem PnLD that is a minimization of the model PnLR in the space of multipliers (PLD n )

well type TLP subsea facility type TLP small FPSO large FPSO

λ

(Pn)

∑λ

construction lead time (years)

LR s,s′ s,s′ s,s′ s,s′ s,s′ φLD 0 ) min φ0 (λbt,t, λbsf,t, λblf,t, λdtw,r,t, λdsw,r,t) (55)

To find the upper bounds at any node n other than the root node in the branch-and-bound tree, model Pn is considered. In model Pn, in addition to the initial nonanticipativity constraints that are added to the model regardless of any decision, some conditional nonanticipativity constraints coming from the relaxed disjunction are also included. The conditional nonanticipativity constraints that apply in node n are included in the dynamic set NnC. The selection of which conditional nonanticipativity constraints to include in set NnC and some necessary cuts to be added to Pn are discussed in section 5.2. The model Pn, not including any of such cuts, is given by

(PLR n )

Table 2. Construction Lead Time (Years)

LR s,s′ s,s′ s,s′ s,s′ s,s′ φLD n ) min φn (λbt,t, λbsf,t, λblf,t, λdtw,r,t, λdsw,r,t) (63) λ

s,s′ s,s′ s,s′ s,s′ s,s′ where λbt,t , λbsf,t , λblf,t , λdtw,r,t , and λdsw,r,t represent the Lagrange multipliers corresponding to constraints 57-61. For fixed values of the multipliers, both models P0LR and PnLR can be decomposed into independent scenarios, each of which is represented by a nonconvex MINLP subproblem. Because both models are relaxations of the model P for any fixed values of the Lagrange multipliers, they yield a valid upper bound if all of the nonconvex subproblems are globally optimized. The validity of these upper bounds is proved in Appendix A.

1 1 1 2 4

Table 3. Representation of Scenarios Using Uncertain Parameters

scenario

initial productivity per well (kbd)

reservoir size (Mbbl)

water breakthrough time

1 2 3 4 5 6 7 8

10 10 20 20 10 10 20 20

300 300 300 300 1500 1500 1500 1500

5 2 5 2 5 2 5 2

Minimization of the Lagrangean dual models P0LD and PnLD in the space of multipliers is performed by the subgradient method proposed by Fisher.33 5.2. Branching. In general, at any node in the branch-andbound tree, the solution of the Lagrangean dual might not satisfy the dualized nonanticipativity constraints (eqs 57-61) or the conditional nonanticipativity constraints in the relaxed disjunction (eq 48) inferred by decisions. In this case, new branches are generated from the current node by considering the violations in the dualized nonanticipativity constraints or the relaxed disjunction. This branching procedure is explained in detail.30 5.3. Lower Bounding Procedure. The lower bound at each node is found by a heuristic, based on the solution found by the upper bound. Usually, the solution found by the upper-bound generation does not satisfy the nonanticipativity constraints. In this case, a feasible solution is generated using a rolling horizon approach.34 Starting from the first period, the probability-weighted averages of variables in indistinguishable scenarios are found (similar to that performed in branching; see step 7 in the algorithm). Assuming that these variables are fixed and considering the resolution of uncertainty depending on the fixed decisions, the decisions for the next period are found iteratively by calculating the probability-weighted averages of the variables in indistinguishable scenarios. This approach continues until all scenarios are distinguishable or the end of the planning horizon is reached, whichever comes first. Then, the model is solved in full space with these fixed decisions using an outer approximation procedure. This solution yields a feasible set of decisions on the numbers of wells to drill and facilities to build, the times to drill and build them, and the oil production for every time period. 5.4. Solution Algorithm (SP-GO). For convenience in the presentation, generic variable names are used in the explanation of the algorithm. Two criteria are taken into consideration for selecting the generic variable names. The first is that a variable can be discrete or continuous, and the second is that a variable can be a decision, state, or recourse variable. The discrete decision variables are denoted by d, whereas discrete state variables denoted by y. Continuous variables are denoted by w, z, or x: w for a decision variable, z for a state variable, and x for a recourse variable. Note that there is no discrete recourse variable. Based on the above definitions, the proposed algorithm is presented below. P denotes the list of open nodes, each having an upper bound φnUB found by the Lagrangean dual problem, and φLB represents


the objective value of the best feasible solution obtained so far. Given these variables, the major steps of the algorithm are as follows: Step 1. Initialization: φLB ) -∞. Step 2. Bound Generation: For the root node, generate upper (φ0UB) and lower (φLB) bounds with solutions (dˆ,yˆ,wˆ,zˆ,xˆ) and (dj,yj,w j ,zj,xj), respectively. [Details of this step (steps a-g) are explained below.] Step 3. Branching: Branch on the dualized nonanticipativity constraints or disjunctions that are violated by the solution (dˆ,yˆ,wˆ,zˆ,xˆ) of the relaxed problem P0 and generate two children nodes. Step 4. Bound Generation: For each child node n, generate upper (φnUB) and lower (φLB) bounds with solutions (dˆ,yˆ,wˆ,zˆ,xˆ) and (dj,yj,w j ,zj,xj), respectively. (Details of this step for each child node are similar to step 2). Step 5. Termination: If P ) L, stop. The current best solution is optimal. Otherwise, repeat steps 6-8. Step 6. Node Selection: Select and delete node n from P based on the best bound. Step 7. Branching: Branch on the dualized nonanticipativity constraints or disjunctions that are violated by the solution (dˆ,yˆ,wˆ,zˆ,xˆ) of the relaxed problem Pn. Generate two children nodes; add them to P. Step 8. Bound Generation: For each child node n, generate upper (φnUB) and lower (φLB) bounds with solutions (dˆ,yˆ,wˆ,zˆ,xˆ) and (dj,yj,w j ,zj,xj), respectively. (Details of this step for each child node are similar to step 2.) Note that steps 1-4 are a part of the initialization where the bounds for the root node and two children nodes are found. The algorithm iterates steps 5-8 during the rest of the algorithm. Steps 2, 4, and 8 share the same procedure for finding bounds. The steps of this common procedure (steps a-g) are explained below. Step a. Set iteration i ) 0. Step b. While i is less than or equal to max_iteration, i ) i + 1, repeat steps b-g. Step c. Generate Upper Bound: For fixed multipliers, use the global optimizer for each subproblem to solve PnLR to obtain solution (dˆ,yˆ,wˆ,zˆ,xˆ) with objective function value φˆ . Step d. Update Upper Bound: Update the upper bound by φUB ) min{φUB,φˆ }. Step e. Generate Lower Bound: Find a heuristic solution based on (dˆ,yˆ,wˆ,zˆ,xˆ) to generate a feasible solution (dj,yj,w j ,zj,xj) with j. objective value φ Step f. Update Lower Bound: Update the lower bound by j }. Delete from P all problems P′ with φˆ (P′) φLB ) max{φLB,φ LB eφ . Step g. Update Multipliers: Update multipliers using subgradient method. 6. Results In this section, we present results for specific example problems to show the effectiveness of the proposed solution algorithm, SP-GO. (For the generic problem statement, see section 2.) In the first example (section 6.1), planning decisions are optimized for a single reservoir for 10 years where uncertainty is represented using eight different scenarios. The oil flowrate from a single well in the reservoir is assumed to be given by a linear function of cumulative oil production from the same reservoir. In the second example problem (section 6.2), unlike in example 1, the oil flowrate from a single well in the reservoir is approximated using a nonlinear function of the cumulative production from the reservoir. Note that, in both

Figure 9. Maximum oil flowrate behavior under uncertainties in initial productivities and reservoir sizes.

Figure 10. Water flowrates for eight different scenarios. Table 4. Model Size for Example 1

integer variables binary variables continuous variables constraints

individual scenario

full-space model

100 10 121 271

800 600 869 8088

examples, the reservoir model is nonlinear because of the nonlinearity caused by the water-to-oil ratio and the water flowrate. All reported results were obtained on a Pentium IV, 3.20 GHz Windows machine. Also, we employed AIMMS 3.8 for implementing the solution algorithm using solvers CPLEX 10.1 for MILP, CONOPT 3.14 and SNOPT 6.1 for NLP, BARON 7.5.3 for the global optimization of MINLP, and AOA (AIMMS Outer Approximation Module) for local optimization of the MINLP. 6.1. Example 1. In this example, our aim is to optimize the infrastructure planning decisions for an offshore oil field having a single reservoir for 10 years. The planning decisions are the number, capacity, and installation schedule of FPSO/TLP facilities; the number and drilling schedule of subsea/TLP wells; and the oil production profile over time. There are certain specifications about the actions that can be taken during the planning process. Subsea and TLP wells are drilled in groups of three. Because of the limitations on the number of drilling rigs, the maximum number of subsea wells that can be drilled each year is 12. Also, each TLP facility can drill at most 6 TLP wells each year. A maximum of 30 TLP wells can be connected to each TLP facility. The construction time delays for different types of wells and facilities are listed in Table 2. The uncertainties in the initial maximum oil flowrate, the size of the reservoirs, and the water breakthrough time (see section 2 for the correlation between the water breakthrough time and

Ind. Eng. Chem. Res., Vol. 48, No. 6, 2009 3087 Table 5. Expected Value Solution Found for Example 1

parameter R) are represented by discrete distributions consisting of two values (high and low). These uncertainties are incorporated into the model using eight different scenarios. Each scenario is a unique combination of the possible values of uncertain parameters (see Table 3). In example 1, each scenario is assumed to be equally likely. We assume a linear relation between the maximum oil flowrate and the cumulative oil production from the reservoir as shown in Figure 9. Therefore, combination of the two different values (low and high) of initial productivities and the two different values (low and high) of reservoir sizes gives four different lines representing the maximum oil production flowrate. For each maximum oil flowrate line, there can be two possible curves for the water flowrate depending on the value of the water breakthrough parameter. This leads to eight different water flowrate curves, each of which corresponds to a scenario, as shown in Figure 10.

The specific values used for describing the uncertainty resolution are given as follows: The appraisal program is completed when a total of three wells are drilled (N1 ) 3) in one reservoir. As explained before, this appraisal program not only gives the actual value for the initial maximum oil flowrate, but also provides the posterior probabilities of reservoir sizes depending on the outcome. The uncertainty in reservoir size can be resolved if either a total of nine (N2 ) 9) or more wells are drilled or production is made from that reservoir for a duration of 1 year (N3 ) 1). Uncertainty in the water breakthrough time is resolved after 1 year of production from the reservoir (N4 ) 1). In Table 4, a comparison of model statistics is given for the full-space model and the individual scenarios for 10 time periods. We should note that the size of the full space model increases exponentially as a result of the increase in the number of binary variables for representing uncertainty resolution and

3088 Ind. Eng. Chem. Res., Vol. 48, No. 6, 2009 Table 6. Stochastic Programming Solution Found for Example 1

Table 7. Comparison of Capital Investments number of units in solutions expected value approach cost per unit ($ × 106) TLP facilities small FPSO facilities large FPSO facilities TLP wells subsea wells total capital expenditure ($ × 109)

250 700 1200 20 30

min 2 5 0 12 12 4.6

the nonanticipativity constraints that relate the decisions in indistinguishable scenarios.

average 2.8 6.5 0 57 20.3 7

max 5 8 0 150 33 10.8

stochastic programming approach (SP-GO) min 1 2 0 12 9 2.16

average 2.6 4.6 0.25 57.7 20.6 5.9

max 6 6 1 159 33 11.07

A simplified approach for finding a feasible solution to this planning problem is to use expected values, in which all discrete


Figure 11. Resulting branch-and-bound tree for example 1.

decisions for the entire planning horizon are first optimized using one-point estimates (mean, most likely, etc.) of the uncertain parameters. Then, starting from the first time period, discrete decisions are fixed successively, resolution of uncertainties are observed, and the rest of the planning horizon is reoptimized throughout the planning horizon. The expected solution was found by using mean values of the initial productivity, reservoir size, and breakthrough time parameter (15, 900, 3.5). The expected value solution (see Table 5) proposes building five small FPSO and two TLP facilities and drilling nine subsea wells in the first year. These decisions resolve the uncertainty in the initial productivity and reservoir size. Depending on the values of reservoir size and initial productivity, different decisions are implemented. This expected value approach gives an objective function value of $5.81 × 109. As seen in Table 6, the optimal stochastic programming solution yields an expected net present value of $6.37 × 109, which is higher than the expected value solution ($5.81 × 109). The solution proposes building two small FPSO and one TLP facilities and drilling nine subsea wells in the first year. Uncertainty in the initial oil flowrate and reservoir size is resolved after nine wells have been drilled. For scenarios 5 and 6, the solution proposes building four more small FPSO facilities, one large FPSO facility, and five TLP facilities and drilling 12 subsea and three TLP wells. For scenarios 7 and 8, the solution proposes building six more small FPSO facilities and one TLP facility and drilling six subsea wells. For solving each subproblem globally in a reasonable time frame, we customized the options for solving different scenarios depending on the time elapsed in different parts of the solver BARON (i.e., preprocessing, local search, bounds tightening, marginals testing, probing). This leads to a savings of up to 85% in the solution time compared to the default settings. Figure 11 shows the branch-and-bound tree generated by the proposed algorithm, which found the optimal solution with an expected net present value of $6.37 × 109. The root node is denoted by 0, and the rest of the nodes are numbered according to the order of exploration. At the root node, the Lagrangean problem yields an upper bound of $6.49 × 109 after five iterations of subgradient optimization, and the heuristic for finding a feasible solution yields a lower bound of $6.33 × 109. In the optimal solution of the Lagrangean problem at the root node, the highest violations in the initial nonanticipativity constraints are drilling groups of subsea wells. Therefore, the feasible region is separated into two, which are represented by nodes 1 and 2. The best feasible solution proposed by the algorithm is guaranteed to be within a 9% optimality gap, and the proposed branch-and-bound algorithm required 23 h because, at each node, 40 MINLP problems were solved to global optimality. A total of seven nodes were traversed, and the best feasible solution

Figure 12. Comparison of net present values.

was found at node 5. In this instance, the modified options of the solver BARON35 were used for solving each nonconvex MINLP with a 10% gap. Figure 12 compares the net present values of the expected value solution (patterned columns) and the stochastic programming solution (black columns) over the eight scenarios in the 10-year period. The added value of stochastic programming is due to the conservative initial investment strategy compared to the expected value solution strategy. The stochastic programming approach considers all eight scenarios before making the initial investment; therefore, it proposes building two small FPSO facilities instead of five and one TLP facility instead of two. Also, it builds more facilities and drills wells only after it determines that the reservoir size is 1500 Mbbl. A comparison of net present values of the expected value solution and stochastic programming shows that the added value of stochastic programming comes from handling the downside risk much better than the expected value solution. Table 7 compares the minimum, maximum, and average capital investment decisions made by the expected value and stochastic programming approaches over all scenarios. The variation in the amount of capital investment is higher in the stochastic programming approach because a smaller investment is made in unfavorable scenarios, whereas larger investments are made in favorable scenarios. Furthermore, the average investment proposed by stochastic programming is less than that from the expected value solution. Therefore, although the stochastic programming solution proposes smaller investments, it creates a higher expected net present value. 6.2. Example 2. In this example problem, the maximum oil flowrate from a single well in a reservoir is approximated using a quadratic function of the cumulative production from the

3090 Ind. Eng. Chem. Res., Vol. 48, No. 6, 2009 Table 8. Expected Value Solution Found for Example 2

reservoir. The basic idea is to investigate the sensitivity of the solution to this oil flowrate curve. γr,s 1

γr,s 3

1-λ-µ ) r,s 2 λ(1 - λ) (REC )

γr,s 2 )

γr,s 3 λ2 + µ - 1 RECr,s λ(1 - λ)

∀s, ∀r

∀s, ∀r

(64)

(65)

γ1r,s, γ2r,s, and γ3r,s are as given in eq 16. λ and µ are parameters specifying the third point (λRECr,s, µγ3r,s) that the oil flowrate passes in addition to the two prespecified points (0, γ3r,s) and (RECr,s,0). It is also more realistic to assume a convex curve such that the oil flowrate initially decreases more rapidly. This necessitates a third point satisfying λ + µ < 1. For λ ) 0.5 and

µ ) 0.25, the maximum oil flowrates from a single well in a reservoir for different scenarios are shown in Figure 13. Given the eight scenarios over the 10-year horizon, the problem data used in example 1, and the specifications for the nonlinear oil profile, the expected value solution (see Table 8) proposes building four small FPSO and three TLP facilities in the first year and starting to drill 12 TLP wells in year 2. These decisions resolve the uncertainty in the initial productivity and size in year 3. The production starts in year 3, and after all the scenarios become distinguishable in year 4, different decisions are implemented for the rest of the time horizon to maximize the net present value. This expected value approach gives an objective function value of $3.76 × 109. As seen in Table 9, the optimal stochastic programming solution yields an expected net present value of $4.59 × 109


Figure 13. Maximum oil flowrate behavior under uncertainties in initial productivities and reservoir sizes.

versus $3.76 × 109 from the expected value solution. The solution proposes building two small FPSO and one TLP facilities and drilling nine subsea wells initially. Drilling these nine wells will resolve the uncertainty in the initial oil flowrate

Figure 14. Comparison of net present values.

Table 9. Stochastic Programming Solution Found for Example 2 scenario year

1

2

3

4

1

5

6

7

8

build: 2 small FPSO 1 TLP drill: 9 subsea wells

2

build: 4 small FPSO 5 TLP drill: 12 subsea wells

drill: 3 TLP wells 3

drill: 6 TLP wells

build: 4 small FPSO 3 TLP drill: 6 TLP wells 6 subsea wells drill: 6 subsea wells 18 TLP wells

drill: 12 subsea wells 27 TLP wells

4 drill: 12 subsea wells 6 TLP wells


drill: 6 TLP wells


drill: 24 TLP wells

drill: 30 TLP wells

build: 1 small FPSO drill: 3 TLP wells

5

drill: 6 TLP wells

drill: 6 TLP wells

drill: 6 TLP wells

drill: 6 TLP wells

drill: 30 TLP wells

drill: 36 TLP wells

drill: 24 TLP wells

drill: 24 TLP wells

6

drill: 6 TLP wells

drill: 6 TLP wells

drill: 24 TLP wells

drill: 36 TLP wells

drill: 15 TLP wells


drill: 21 TLP wells

drill: 36 TLP wells

drill: 18 TLP wells

drill: 24 TLP wells

drill: 3 TLP wells

drill: 15 TLP wells

8.56

9.17

7 8

drill: 9 TLP wells

9-10 NPV ($×109) ENPV($×109)

0.83 4.59

0.83

1.86

1.87

6.69

6.89

Table 10. Comparison of Capital Investments number of units in solutions expected value approach

TLP facilities small FPSO facilities large FPSO facilities TLP wells subsea wells total capital expenditure ($ × 109)

cost per unit ($ × 106)

min

average

250 700 1200 20 30

3 4 0 39 0 4.33

3.50 5.00 0.00 77.25 5.63 6.08

max 4 5 0 111 18 7.26

stochastic programming approach (SP-GO) min

average

max

1 2 0 21 9 2.34

3 4.1 0 72 21.4 5.72

6 6 0 165 33 9.99


and reservoir size. For scenarios 5 and 6, the solution proposes building four small FPSO facilities and five TLP facilities and drilling 12 subsea wells. For scenarios 7 and 8, the solution proposes building four small FPSO facilities and three TLP facilities and drilling six TLP wells and six subsea wells. The best feasible solution proposed by the algorithm after 120 h is guaranteed to be within a 12% optimality gap. The reason for such a gap is that each subproblem is solved with a 10% gap, so in the worst case, the upper bound generated will be 10% better. A total of 14 nodes were traversed, and the best feasible solution was found after 90 h at node 6. Notice that, in this instance, some of the default options of the solver BARON (i.e., preprocessing, local search, bounds tightening, marginals testing, probing) were modified for solving each scenario subproblem. This leads to a 1-orderof-magnitude reduction in solution time compared to the default options. Figure 14 compares the net present values of the expected value solution (patterned columns) and the stochastic programming solution (black columns) over the eight scenarios in the 10-year horizon. As for example 1, the added value of stochastic programming is due to the conservative initial investment strategy compared to the expected value solution strategy. Comparing the result of example 2 with that of example 1 shows that the improvement by stochastic programming on the expected value solution is higher (22% vs 10%) when the conditions for extracting oil from the reservoir become harder. Note that, in the nonlinear case, it becomes harder to extract oil because the decrease in the oil flowrate is much faster than in the linear case. Table 10 compares the capital investment decisions made by the expected value and stochastic programming approaches for example 2. Although the average capital investments are similar in the two approaches, the stochastic programming approach gives lower minimum and higher maximum investments. This smaller investment is due to the conservative initial investment strategy proposed by the stochastic programming approach. Also, stochastic programming proposes higher investments only when the uncertainty parameters are found to be favorable. Comparing the solution quality of the two examples, the stochastic programming solution is clearly more robust with respect to the uncertainty in the shape of the curve representing the maximum oil flowrate versus the cumulative production from the reservoir. The first-year decisions proposed by stochastic programming are the same in the two examples, whereas the expected value solution decisions change. Stochastic programming is generating solutions that not only provide higher expected net present values, but also offer more robust solutions with respect to the uncertainties that were not included explicitly in the uncertainty space (e.g., uncertainty in shape of maximum oil flowrate). Although this result was not included in the objective explicitly, this finding supports the advantage of using stochastic programming. 7. Conclusions We have presented in this article a multistage stochastic programming model for the planning of offshore oil field infrastructure under uncertainty where the uncertainties, namely, the initial maximum oil flowrate, recoverable oil volume, and water breakthrough time of the reservoir, reduce gradually as a function of design and operating decisions. The proposed model considers possible investment strategies for reducing uncertainties before making the major facility investments. The proposed model is a disjunctive/mixed-integer nonlinear programming

model that is converted into an MINLP using big-M transformations. We have proposed a duality-based branch-and-bound algorithm that enables solving the scenario subproblems independently. Feasible solutions are generated using the infeasible solutions coming from individual scenarios. Implementation issues and results for two examples were presented to illustrate the application of the proposed method (SP-GO). In examples 1 and 2, the proposed algorithm found solutions that are nearly 10% and 22% better, respectively, than the solutions found by the expected value approach. The reason for these improvements is that stochastic programming incorporates the uncertainty and value of information analysis directly into the mathematical model and proposes decisions that hedge against possible outcomes of uncertain parameters. The solution times in both examples are rather long. In a forthcoming article, we will address how to improve the computational efficiency.36 Acknowledgment The authors acknowledge financial support from ExxonMobil Upstream Research Company and partial support from the National Science Foundation under Grant CTS-0521769. Nomenclature Sets and Indices D(s,s′) ) set of reservoirs that differentiate scenarios s and s′ M1 ) set of scenario pairs (s,s′) that differ only in initial productivities M2 ) set of scenario pairs (s,s′) that differ only in reservoir sizes M3 ) set of scenario pairs (s,s′) that differ only in water breakthrough times M4 ) set of scenario pairs (s,s′) that differ only in initial productivities and reservoir sizes M5 ) set of scenario pairs (s,s′) that differ only in initial productivities and water breakthrough times M6 ) set of scenario pairs (s,s′) that differ only in reservoir sizes and water breakthrough times M7 ) set of scenario pairs (s,s′) that differ in initial productivities, reservoir sizes and water breakthrough times NI ) set of initial nonanticipativity constraints NC ) set of conditional non-anticipativity constraints NCn ) set of conditional nonanticipativity constraints at node n R ) set of reservoirs r ) reservoir r ∈ R S ) set of scenarios s ) scenario s ∈ S T ) planning horizon t, τ ) time periods t ∈ T Continuous Variables capexts ) capital expenditures in time period t in scenario s costts ) total of capital and operating expenditures in time period t in scenario s enpv ) expected net present value of the project npvs ) net present value of the project in scenario s oiltr,s ) oil flowrate from one well in reservoir r in time period t in scenario s oilcum,r,s ) cumulative amount of oil extracted from reservoir r until t the end of time period t in scenario s oiltfpso,r,s ) amount of oil extracted from reservoir r and pumped to FPSO facilities in time period t in scenario s oilttlp,r,s ) amount of oil extracted from reservoir r and pumped to TLP facilities in time period t in scenario s opexts ) operating expenditures in time period t in scenario s

Ind. Eng. Chem. Res., Vol. 48, No. 6, 2009 3093 ) amount of oil produced from reservoir r in time period t in scenario s prodttotal,s ) total amount of oil produced in time period t in scenario s revst ) revenue coming from the sale of produced oil in time period t in scenario s watertfpso,r,s ) amount of water extracted from reservoir r and pumped to FPSO facilities in time period t in scenario s watertlp,r,s ) amount of water extracted from reservoir r and pumped t to TLP facilities in time period t in scenario s Integer Variables blfts ) number of large FPSO facilities built in time period t in scenario s bsfts ) number of small FPSO facilities built in time period t in scenario s btts ) number of TLP facilities built in time period t in scenario s dswr,s t ) number of subsea wells drilled in reservoir r in time period t in scenario s; it takes a value multiple of Gsw dtwtr,s ) number of TLP wells drilled in reservoir r in time period t in scenario s; it takes a value multiple of Gtw nlfts ) number of large FPSO facilities available in time period t in scenario s nsfts ) number of small FPSO facilities available in time period t in scenario s nswtr,s ) number of subsea wells available (ready to operate) in reservoir r in time period t in scenario s; it takes a value multiple of Gsw ntts ) number of TLP facilities available in time period t in scenario s ntwtr,s ) number of TLP wells available (ready to operate) in reservoir r in time period t in scenario s; it takes a value multiple of Gtw Binary/Logic Variables bptr,s ) true if there is production from reservoir r in time period t in scenario s bwr,s t ) true if there is a well ready for extracting oil from reservoir r in time period t in scenario s r,s ) true if and only if the number of wells drilled in reservoir r w1,t up to time period t in scenario s is less than N1 r,s w2,t ) true if and only if the number of wells drilled in reservoir r up to time period t in scenario s is less than N2 r,s w3,t ) true if the number of years of production is less than the amount of time needed to resolve uncertainty in reservoir size (i.e., N3) at the beginning of time period t in scenario s wr,s 4,t ) true if the number of years of production is less than amount of time needed to resolve uncertainty in water breakthrough time (i.e., N4) at beginning of time period t in scenario s zs,s′ t ) true if scenario s and s′ are indistinguishable at the beginning of time period t Parameters C1,t ) price of oil at time period t C2,t ) capital cost of drilling a TLP well connecting it to a TLP facility at time period t C3,t ) capital cost of drilling a subsea well at time period t C4,t ) capital cost of building a TLP facility at time period t C5,t ) capital cost of building a small FPSO facility at time period t C6,t ) capital cost of building a large FPSO facility at time period t C7,t ) cost of production per unit of product at time period t CAPlfliq ) maximum capacity of a large FPSO facility for liquid (oil and water) CAPlfoil ) maximum capacity of a large FPSO facility for oil prodtr,s

) maximum capacity of a small FPSO facility for liquid (oil and water) CAPsfoil ) maximum capacity of a small FPSO facility for oil Dt ) discounting factor at time period t Gsw ) number of subsea wells to be drilled together as one group Gtw ) number of TLP wells to be drilled together as one group N1 ) number of wells needed to complete appraisal program N2 ) total number of wells needed to be drilled to resolve uncertainty in reservoir size N3 ) number of years of production needed to resolve uncertainty in reservoir size N4 ) number of years of production needed to resolve uncertainty in water breakthrough time Nwmax ) maximum number of wells that can be connected to a TLP facility Ps ) probability of scenario s RECr,s ) amount of recoverable oil in reservoir r in scenario s Usf ) upper bound on the number of small FPSO facilities to build Usw ) upper bound on the number of subsea wells to drill in one time period Ut ) upper bound on the number of TLP facilities to build Ulf ) upper bound on the number of large FPSO facilities to build Utw ) upper bound on the number of TLP wells to drill in one time period per TLP facility Uw,r ) upper bound on the number of wells to drill in reservoir r Rr,s ) linear coefficient in the term for calculating water-to-oil ratio β′ ) exponent in the term for calculating water-to-oil ratio γ1r,s ) quadratic-term coefficient in the equation representing the oil flowrate from a well in reservoir r in scenario s γr,s 2 ) linear-term coefficient in the equation representing oil flowrate from a well in reservoir r in scenario s γ3r,s ) constant coefficient of the equation representing oil flowrate from a well in reservoir r in scenario s; also represents initial oil flowrate from a well in reservoir r in scenario s δt ) operating days per period τ() ) time delays CAPsfliq

Appendix A: On the Validity of the Upper Bound from Lagrangean Duality A stochastic program can be represented by eqs A1-A4, where the objective function is the sum of convex or nonconvex functions over all scenarios S. Constraint set A2 represents a nonlinear, nonconvex feasible space for scenario s. Constraint set A3 represents nonanticipativity constraints linking decisions in different scenarios. max xs

∑ f (x ) s

s

(A1)

s∈S

gs(xs) e 0 xs - xs′ ) 0 xs ∈ R

n

∀s ∈ S

(A2)

∀(s, s′) ∈ N

(A3)

∀s ∈ S

(A4)

Variables for different scenarios are linked through the nonanticipativity constraints A3, and if they are dualized, the model can be decomposed into independent subproblems for each scenario s. Proposition 1 If the nonanticipativity constraints A3 are dualized, the Lagrangean dual function is convex in the space of the multipliers. Proof Dualizing nonanticipativity constraints A3 yields


max xs

∑ f (x ) + ∑ s

λs,s′(xs - xs′)

s

s∈S

(A5)

or equal to the sum of individual objective function values of global optimal solutions (i.e., wait-and-see solution).

(A6)

Lemma 2

(s,s′)∈N

∀s ∈ S

gs(xs) e 0

Then, the model can be written as the sum of independent subproblems max xs

∑ [f (x ) + x ( ∑ λ s

s

s

s,s′

s∈S

∑λ

-

s

Stochastic Programming Approach for the Planning of Offshore Oil or

Recommend Documents