A Lagrangean Decomposition Heuristic for the Design and Planning of

Ind. Eng. Chem. Res. 2001, 40, 2857-2875

2857

A Lagrangean Decomposition Heuristic for the Design and Planning of Offshore Hydrocarbon Field Infrastructures with Complex Economic Objectives Susara A. van den Heever and Ignacio E. Grossmann* Department of Chemical Engineering, Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, Pennsylvania 15213

Sriram Vasantharajan† and Krisanne Edwards‡ Mobil Technology Company

In this work, we present a multiperiod mixed-integer nonlinear programming model for the long-term design and planning of offshore hydrocarbon field infrastructures with complex economic objectives. We show that when complex fiscal rules, such as tariff, tax, and royalty calculations are added to the model, substantial increases in the net present value of the project are obtained. However, including these calculations leads to more than an order of magnitude increase in the computational effort required. To address this problem, we propose a specialized heuristic algorithm that relies on the concept of Lagrangean decomposition. We present an iterative scheme where feasible solutions are postulated from the Lagrangean decomposition, and multipliers are updated through a subgradient method. Results show that for moderately sized models up to 2 orders of magnitude of reduction in solution time is obtained compared to the full-space solution. Larger problems that could not be solved in 5 days of computation in full space are solved in a few hours with the proposed method. Introduction Oil and gas exploration and production is a huge industry, with typical oil production being on the order of 300 000 barrels/day and gas production on the order of 1700 million ft3/day worldwide for one of the leading companies.1 Around 1000 wells are drilled annually to sustain this production, leading to typical annual property acquisition, exploration, and development costs of $8 billion and net annual exploration and production earnings on the order of $5 billion worldwide. An offshore hydrocarbon field infrastructure typically consists of production platforms (PPs), well platforms (WPs), and interconnecting pipelines (see Figure 1) from which oil and/or gas is produced. The construction of offshore hydrocarbon field infrastructures gives rise to a very complex problem that involves design and planning decisions, as well as fiscal considerations. Given the large economic impact of these decisions, there is a clear motivation to develop a systematic optimization method that can effectively identify optimal designs and plans that capture all of the complex economic tradeoffs. At the design level discrete decisions include which PPs, WPs, and pipeline connections to install over the entire project lifetime, which can be from several years to several decades. Continuous design decisions include the capacities of these items. To facilitate decision making at the planning level, the operating horizon is divided into a number of time periods, for example, months or years, and discrete and continuous decisions * To whom correspondence should be addressed. E-mail: [email protected]. Tel: (412) 268-2230. Fax: (412) 2687139. † Now with Optimal Decision LLC. ‡ Now with ExxonMobil Upstream Research Co.

Figure 1. An offshore hydrocarbon field infrastructure with a PP, WPs, and connecting pipelines.

are made in every period. Discrete planning decisions include when to install the PPs, WPs, and pipelines, while continuous decisions include the production profile of each WP in each time period. Nonlinear reservoir behavior is represented by simplified nonlinear relationships of, for example, reservoir pressure vs cumulative hydrocarbon produced and surface pressure vs reservoir pressure and deliverability. While in the past simple revenue and cost functions have been used,2-5 complexities that arise from royalty, tax, and tariff rules specified by the corresponding government of the country where the hydrocarbon is extracted from have not been addressed. Including complex economic specifications, for example, taxes, tariffs, and royalties, into an optimization model leads to more than an order of magnitude increase in solution time, as will be demonstrated in this paper. This is due to the fact that often the complex

10.1021/ie000755e CCC: $20.00 © 2001 American Chemical Society Published on Web 06/02/2001

2858

Ind. Eng. Chem. Res., Vol. 40, No. 13, 2001

fiscal rules as specified by governments cannot be modeled in a simple form and require additional discrete variables and complex mixed-integer constraints. Neglecting these rules leads to economically suboptimal solutions. This effect is magnified in the hydrocarbon industry, where even a 1% difference in the objective can lead to millions of dollars in savings or losses. In addition, taxes, royalties, and tariffs are present in other areas, such as process network planning and supply chain planning, making the problem challenging and relevant to other areas. In this work, we present a multiperiod mixed-integer nonlinear programming (MINLP) model for the design and planning of offshore hydrocarbon field infrastructures with complex objectives. This model is capable of dealing with complex rules that specify different tax, royalty, and tariff factors depending on some economic indicator, such as net revenue or gross revenue. The economic part of the model is general and can be applied to industries other than the hydrocarbon industry. Because of the large computational expense involved in solving the proposed model, a Lagrangean decomposition heuristic is applied to this problem. Background The discrete and continuous nature of the variables, together with the nonlinear reservoir behavior, requires an MINLP model, where the objective is to optimize a specific economic function. Given the computational expense of solving multiperiod MINLP models, these have traditionally been avoided in favor of mixed-integer linear programs (MILPs) or linear programs (LPs). In the past, decisions regarding the design and planning of offshore hydrocarbon fields have often been made separately to ease the computational burden.6-9 The consequence of not including all of the design and planning decisions in one model is that interactions between these variables are not taken into account, leading to infeasible or suboptimal solutions. Simultaneous models have only appeared recently.2-4,10,11 Initially, attempts were made to solve these models with commercial MILP solvers, but only small instances could be solved this way.2,3,11 In the last 2 or 3 years, specialized algorithms emerged but only for linear models.4,10 In all of the methods mentioned thus far, nonlinearities have been avoided by either assuming simplified linear behavior,3,9 fixing an estimated production profile,10 or using piecewise linear interpolation.4,11 The first assumption is not realistic, while the second may neglect important interactions and lead to suboptimal solutions. Linear interpolation leads to improved accuracy, but this is offset by the computational effort resulting from the large increase in the number of integer variables. It is, therefore, important to deal with nonlinearities directly. In previous work5 the authors addressed the issues of problem size and nonlinear reservoir behavior through a multiperiod disjunctive optimization model. A bilevel decomposition technique,12 logic-based outer approximation (OA) algorithm,13 and the aggregation of time periods were combined into an iterative aggregation/ disaggregation algorithm to deal with the problem size. Convex underestimators14 were employed to deal with nonconvexities arising from the nonlinear reservoir behavior. It was shown that, by application of the proposed method, more than an order of magnitude reduction in solution time could be achieved compared

to a commercial solver such as DICOPT++.15 However, complex economic rules were not included in the objective function, because they lead to an intractable model. This has motivated the current study of alternative models for economic calculations. Fiscal considerations such as taxes, tariffs, and royalties have been largely simplified in the past. One such simplification is to work with an average profit per unit production in order to maximize the income as a function of profit per barrel multiplied by the number of barrels of oil produced.3,7,9 Other typical objectives include only sales revenue, operating costs, and investment costs.2,4,6 Sullivan11 includes royalties and taxes in the objective function but as a fixed percentage of some economic measure, such as the revenue. In reality, however, these percentages might not be fixed but change according to the profitability of the project. Frair8 addresses the issue of special tax considerations where cost depletion is calculated with two different factors, depending on whether the net revenue is positive or negative. A big-M formulation was proposed. This means the constraint takes the form g(x) e M(1 y) where g(x) e 0 is the original constraint, M is a big-M parameter (upper bound), and y is a binary variable. If y equals 1, the right-hand side terms goes to zero and the original constraint is applied. If y equals 0, the constraint is relaxed as g(x) e M so as to not influence the optimization. The author stated that it is uncertain whether the inclusion of detailed taxes in the model is justified, because of the large increase in computational effort. Mixed-integer optimization models often consist of some “easy” constraints and some “hard” constraints that are largely responsible for the intractibility. A wellstudied method for dealing with this is Lagrangean relaxation. The basic idea behind Lagrangean relaxation is that the “hard” constraints are dualized in the sense that the constraints are transferred into the objective function in the form of a penalty function, thus leaving the remaining model either relatively easy to solve or decomposable into easy to solve subproblems. In addition, the optimal solution from the Lagrangean problem provides an upper bound to the original problem if the objective is maximized. The term “Lagrangean relaxation” comes from the fact that the penalty term consists of the sum of products of the deviation in the removed “hard” constraints multiplied by their associated Lagrange multipliers. Fisher gives an excellent discussion on Lagrangean relaxation and a review on its practical application.16,17 Lagrangean relaxation is often used as part of a branch and bound algorithm in which it replaces the LP relaxation to provide bounds or as part of a heuristic in order to obtain good feasible solutions. Examples of the former approach in the literature include a capacitated network design problem18 and a pooling problem,19 while examples of the latter approach include a midterm planning problem20 and the capacity planning of phased implementation of flexible manufacturing systems.21 Numerous other authors report on the use of Lagrangean relaxation, for example, as applied to stochastic problems,22,23 unit commitment in power plants,24 and scheduling problems,25,26 to name a few. Lagrangean decomposition27 is a special case of Lagrangean relaxation, particularly applicable where the original model consists of two or more structured constraint sets. Identical copies are made of variables

Ind. Eng. Chem. Res., Vol. 40, No. 13, 2001 2859

Figure 2. Superstructure for 16 WPs.

that occur in both constraint sets, and an equality is added to set the original variables equal to the copies. When this additional constraint is dualized, the model decomposes into a separate model for each constraint set, and these can be solved individually in order to improve on the computational effort required. The method of “variable splitting”28 is similar. Lagrangean decomposition has successfully been applied to, for example, a new plant location problem,29 a combined transportation and scheduling problem,30 and the routing of container ships.31 In the hydrocarbon field industry, each WP often has its own tax, tariff, and royalty structure. This makes the overall model ideally suited for Lagrangean decomposition, because the most time consuming parts of the model, namely, the complex fiscal rules for each platform, are disconnected and the only link between platforms comes from the mass and pressure balances. In this work, we propose an iterative algorithm based on Lagrangean decomposition and a heuristic procedure for finding feasible solutions. A subgradient updating scheme is used to determine the Lagrange multipliers at each iteration. In the next section we present the problem statement, followed by the model. Subsequently, we present an example to illustrate the impact of the complex objectives on both the net present value (NPV) and computational effort. This is followed by a brief background on Lagrangean decomposition, the concept of reformulation of the overall model, and the solution algorithm. Example 1 is revisited to show the improvements achieved with the proposed method. A much larger example follows to illustrate the computational performance on a previously unsolvable gasfield planning model, as well as the improvements in the NPV. Concluding remarks follow. Problem Statement In this paper, we consider the multiperiod design and planning of offshore hydrocarbon field infrastructures (see Figure 2) with the inclusion of complex royalty, tax, and tariff rules. The hydrocarbon fields that we apply our model to are gas-producing fields. The objective is to maximize the NPV, which consists of sales revenues, capital costs, operating costs, taxes, tariffs, and royal-

ties. Gas from each field is produced at the associated WP and sent through pipelines to the PP or to another WP. Where more than one outgoing connection is possible, at most one can be selected. A pipeline connects the PP to the shore where the gas is treated in an existing gas plant (not shown in the superstructure). Decisions associated with the design are valid for the whole project lifetime and are (a) whether to include each PP, WP, and pipeline connection, (b) the type of each WP (if another WP can be connected to it or not), and (c) the capacities of the PPs, WPs, and pipelines. The planning decisions are (a) whether to invest in each PP, WP, and pipeline connection in each time period, (b) the production profile of each WP in each time period, and (c) whether compression is required at each PP in each time period. On the economic side, the government specifies the fiscal rules regarding taxes, tariffs, and royalties. Typically, the amount of taxes, tariffs, and royalties to be paid to the government needs to be calculated as a percentage of one or more objective elements (OE), for example, gross revenue or net revenue. The specific percentage depends on the value of one or more cumulative economic indicators (CEI), for example, the cumulative production, cumulative investment cost, or a combination of these and other indicators. These transition values are expressed in terms of tiers. For example, if the net revenue falls in tier 1, x % taxes need to be paid, while if the revenue falls in tier 2, y % taxes need to be paid, etc. Inclusion of these rules leads to complex economic tradeoffs whose considerations are essential in obtaining an optimal solution. However, such an inclusion requires additional discrete variables and complex mixed-integer constraints, leading to a potentially large increase in computational effort. Our goal is to investigate an effective optimization method to address this issue. Note that the model in this paper is different from the one that was presented by the authors in previous work5 apart from the fact that complex objectives were not included in that model. Not only is this model for a different hydrocarbon field with different physical characteristics, but the design and planning decisions involved, as well as the problem size, are also different.

2860


Model Please refer to the Nomenclature section listed at the end of the paper. Reservoir Model. The gas production from each field, f, in each time period, t, has to be less than the deliverability of the field [see constraint (1)]. Deliverability is determined through a nonlinear relationship with the reservoir pressure, surface pressure, and bottom hole pressure as denoted by the function h1 in constraint (2). The reservoir pressure, in turn, is represented by a nonlinear relationship, h2, in constraint (3) with the cumulative production determined in constraint (4) and the initial pressure in the reservoir. The explicit forms of the nonlinear constraints h1 and h2 cannot be included because of confidentiality reasons, but it suffices to say, for the purpose of this work, that these are either quadratic or cubic functions that are derived from reservoir simulations.

∀ f, t

d

Qf,t e Q

f,t

Qdf,t ) h1(Prf,t;Psf,t;Pbhf,t) r

c

init

P f,t ) h2(Q f,t;P

f)

Qcf,t ) Qcf,t-1 + Qf,t∆t

(1) ∀ f, t

∀ f, t

Qwp′,wp,t ) ∑Qwp,wp′,t + ∑ wp′ wp′ Qwp,pp,t ∑ pp Qpp,t )

∑ wp

Qwp,pp,t

(4)

∀ wp, wp′, t (11)

Pchkwp′,wp,t + Pslkwp′,wp,t

Pinpp,t ) Poutwp,t - Pdropwp,ppQwp,pp,tDistwp,pp ∀ pp, wp, t (12)

Pchkwp,pp,t + Pslkwp,pp,t Poutwp,t ) Pinwp,t - Pdropwp(Qwp,t +

Poutpp,t ) Pinpp,t +

Qwp,wp′,t) ∑ wp′

∑ Pbstpp,cmp,t

Pbstpp,1,t e Pinpp,t(CR1 - 1)

(15)

∀ pp, t

(16)

Poutpp,1,t ) Pinpp,t + Pbstpp,1,t

∀ pp, t

(17)

(6) (7)

Qpp,1,t ) Qpp,2,t + Qbypp,2,t

∀ pp, t

(8)

Pshrt ) Poutpp,t + h3(Qshrt) + Pslkpp,t

∀t

Fuelpp,cmp,t )

Poutpp,2,t ) Poutpp,1,t + Pbstpp,2,t

(9)

(ii) Pressure Balances. The inlet pressure at each WP or PP must be equal for all incoming pipelines. From the reservoir, the inlet pressure at the WP equals the surface pressure minus the pipeline choke pressure as shown in constraint (10). From another WP, the inlet pressure at each WP or PP equals the outlet pressure from the attached WP minus the pressure drop in the pipeline minus the choke pressure as shown in constraints (11) and (12). The slack variable ensures that

∀ pp, t

Pbstpp,2,t e Poutpp,1,t(CR2 - 1)

∀ pp, t

(Qpp,t - ∑ Fuelpp,cmp,t)] ∑ pp cmp

(14)

For each compressor stage, the pressure boost is determined by the compressor ratio and inlet pressure at that stage [constraints (15) and (16)], while the outlet pressure at each stage is the sum of the inlet pressure and pressure boost [constraints (17) and (18)]. The pressure at the shore equals the outlet pressure at a PP plus a nonlinear function of the gas flow rate to the shore as shown in constraint (19). The slack variable ensures that the constraint is only enforced for PPs that are in existence in period t. In constraint (20) the fuel consumed at each compressor stage is calculated by multiplying a factor that includes molecular weight, fuel required per each HP power, and conversion factors, with the power at each compressor stage. The latter is calculated through a standard expression that includes the average compressibility, inlet temperature, compressor ratio, specific heat capacity, and adiabatic compressor efficiency.

Qpp,t ) Qpp,1,t + Qbypp,1,t

Qshrt ) (1 - shrink)[

∀ pp, t

cmp

∀ wp, t (5)

∀ pp, t

(10)

Pinwp,t ) Poutwp′,t - Pdropwp′,wpQwp′,wp,tDistwp′,wp -

(3)

∀ f, t

∀ wp, t

Pinwp,t ) Pswp,t - Pchkwp,t

∀ wp, t (13)

(2)

Surface Model. (i) Mass Balances. The production at each WP plus all flows from other WPs, wp′, into wp must equal the total flow out of that WP to other WPs and to the PP [see constraint (5)]. Similarly, constraint (6) states that the total flow at each PP equals the sum of all flows from WP’s connected to it. Mass balances at each compressor inlet [constraints (7) and (8)] dictate how much gas bypasses the compressor and how much passes through it. The total flow of gas to the shore equals the flow from each PP minus the fuel required to operate the compressor at each PP multiplied by a shrinkage factor to account for losses due to fuel, flare, and so forth [constraint (9)].

Qwp,t +

the pressure balance need not apply if a connection is not in existence, and the corresponding upper bounding constraint is given later by (23). The outlet pressure from each WP is a function of its inlet pressure as well as the pressure drop, while the outlet pressure of the PP is a function of its inlet pressure and pressure boost from the compressor as shown in constraints (13) and (14), respectively.

f

conv

[

[

∀ pp, t

(18) (19)

]

(k -k 1) - 1] (k -k 1)

Qpp,cmp,tZavgTin CRcmp ηadia

∀ pp, t

∀ pp, cmp, t (20)

(iii) Variable Upper Bound Constraints. These constraints set variables to zero if a structure is not in operation in a given time period. In constraint (21) flows across WPs, PPs, and compressors are set to zero if the structure is not in operation in period t, and flows in


pipelines are set to zero if the pipeline connection is not in operation in period t. The bypass stream at each compressor stage is set to zero if that compressor stage is on in period t in constraint (22). Constraint (23) sets the slack pressure to zero if the structure or pipeline is in operation, thus enforcing the pressure constraints (11), (12), and (19). If the structure or pipeline is not in operation (bon(.),t ) 0), the slack pressure can take on large values, thus relaxing the pressure constraint. In constraint (24) the pressure boost at any compressor stage is set to zero if that stage is not in operation in period t.

∀ (.), t

Q(.),t e M(.)bon(.),t

∀ pp, cmp, t

Qbypp,cmp,t e Mpp(1 - bonpp,cmp,t) P

slk (.),t

e

Mslk(.)(1

∀ (.), t

on

-b

(21)

(.),t)

Pbstpp,cmp,t e Mbstcmpbonpp,cmp,t

∀ pp, cmp, t

(22) (23) (24)

(iv) Logical Constraints. A WP, PP, compressor stage, and pipeline can switch on or be installed at most once [constraint (25)]. There can only be one outgoing connection from a WP, and this connection should come on in the same period as the WP is installed as shown in constraint (26). Constraints (27)-(29) specify that a WP or PP must be on for a connection to it to switch on, and any connection can only switch on with flow in one direction. The compressor can only switch on once the PP is on, and compressor stage 2 can only be switched on once stage 1 is on as shown in constraints (30) and (31).

∑t bsw(.),t e 1

∀ (.)

bswwp,wp′,t + ∑bswwp,pp,t ) bswwp,t ∑ wp′ pp

(25) ∀ wp, t

bswwp,pp,t e btypewp,I

∀ wp, pp, t

(35)

∀ pp, t, wp

(28)

b

pp,cmp,t

b

eb

sw

pp,t on

pp,2,t

eb

pp,1,t

(29)

∀ pp, cmp, t

(30)

∀ pp, t

(31)

Constraint (32) determines that a structure or connection is “on” or in operation if it has been switched on or installed in a previous period. These constraints force the “on” variables to take 0-1 values even though they can be declared as continuous.

b

(.),t )

∑

bsw(.),τ

∀ (.), t

bswwp,wp′,t e 1 - btypewp,III ∀ wp, wp′, t if Distwp,wp′ g 10 (36) ∀ wp, tp, t

bswwp,tp,t e bswwp,t

∀ wp, tp, t

bswwp,tp,t e btypewp,tp bswwp,tp,t g bswwp,t + btypewp,tp - 1 bonwp,tp,t e bonwp,t

(32)

τ)1

There are three types of platforms (see Figure 3), and each existing platform has to be assigned only one type as specified in constraint (33). The first type of platform, type I, allows connections from other platforms, i.e., gas inflow from other platforms, and is therefore more expensive. Any WP that

∀ wp, tp, t

∀ wp, tp, t

bonwp,tp,t e btypewp,tp bonwp,tp,t g bonwp,t + btypewp,tp - 1

t

on

(33) (34)

bswwp,pp,t e bonpp,t

on

∀ wp

∀ wp, wp′, t

(27)

sw

btypewp,tp ) ∑bswwp,t ∑ tp t bswwp′,wp,t e btypewp,I

∀ wp, t, wp′

∀ wp, wp′, t

connects to the PP must be a type I platform. These two specifications are represented by constraints (34) and (35). The other two types of WPs do not allow any connections from other platforms and can only be connected to another WP and not to a PP. The two types are distinguished by their respective maximum capacities. The platform with a larger capacity (type II) is more expensive that the one with a smaller capacity (type III), and the choice depends on the largest expected gas flow rate. Constraint (36) is applied only if two WPs are further than 10 miles apart and implies that the connecting WP is in this case not allowed to be of type III (the smaller capacity). Constraints (37)(42) link the type of WP with the 0-1 installation decisions and force the variables bswwp,tp,t and bonwp,tp,t to take 0-1 variables even though they are defined as continuous.

(26)

bswwp′,wp,t e bonwp,t

bswwp,wp′,t + bswwp′,wp,t e 1

Figure 3. WP types I-III.

∀ wp, tp, t ∀ wp, tp, t

(37) (38) (39) (40) (41) (42)

(v) Flow Constraints. If a new maximum flow rate to the shore is selected, it must be maintained for at least 10 years before the flow rate can drop again. If the flow rate in the current period is less than the flow rate in the previous period, the binary variable bt10 must equal zero, implying that constraint (43) will be strict and constraint (44) is relaxed. If a new flow rate results, constraint (43) forces the binary variable to take a value

2862


of 1 and this forces constraint (44) to be strict. Thus, the flow rate in the 10 time periods following an increase in the flow rate must be at least as large as the last maximum flow rate.

∀t

Qshrt e Qshrt-1 + Qshrmaxb10t

(43)

Qshrτ g Qshrt-1 - Qshrmax(1 - b10t) ∀ t, t < τ < t + 10 (44) Economic Model. (i) Revenue. The total revenue in each time period is the sum of revenues from all WPs, which is the sales price of gas in period t multiplied by the flow rate from that wp and the number of days in that period.

Rtott )

Rwp,t ) ∑∆tPtQwp,t ∑ wp wp

∀t

(45)

(CCwp,t + CCpipewp,t) + ∑(CCpp,t + ∑ wp pp cmp ∀t CC pp,t) + ∑CCdevwp,t wp

∑

CCwp,t )

(47)

∑

fwpcc[

τ s.t. t∈Tinv(τ)

(Cconwp,wp′ + ∑ wp′

Cmilewp,wp′Distwp,wp′)bswwp,wp′,t +

(Cconwp,pp + ∑ pp

Cmilewp,ppDistwp,pp)bswwp,pp,t] CCpp,t )

(46)

∀ wp, t

fccwpCcstwp,tpbswwp,tp,τ


CCpipewp,t )

∀ wp, t (48)

∑

f

∑

fccppCcstcmpppbsw1,pp,τ

cc

sw ppCcstppb pp,τ

∀ pp, t (49)


CCcmppp,t )


∀ cmp, pp, t (50) CCdevwp,t )

∑


OCtott ) OCgpt +

OCpp,t + ∑OCwp,t ∑ pp wp

∀t ∀t

OCgpt ) FOcstgpt + VOcstgpQshrt

∀ pp, t

OCpp,t ) FOcstpp,tbonpp,t OCwp,t )

(52) (53) (54)

FOcstwp,tpbonwp,tp,t + VOcstwpQwp,t ∑ tp

∀ wp, t (55)

(ii) Capital Costs. The total capital cost in each time period equals the sum of individual capital costs for all WPs, PPs, pipelines, compressors, and well development costs [constraint (46)]. Constraint (47) calculates the capital cost for installing WP type tp. All capital costs occur over an investment horizon, usually between 1 and 3 years, preceding the actual installation. The constraint implies that a fixed fraction of the total cost is applied in each of these preceding years. Constraint (48) calculates the capital cost for the pipeline out of each WP, while the capital cost for each pp and its associated compressor is calculated in constraints (49) and (50). The cost for well development equals the cost for drilling and completing a well times the number of wells to be drilled at each WP as shown in constraint (51). In this case, the estimated number of wells per platform, Nwellwp, is a known parameter.

CCtott )

fixed charge and a variable charge that is determined by the flow of gas to the shore as shown in constraint (53). A fixed operating cost is applied for each PP in constraint (54), while the operating cost of each WP in constraint (55) is made up out of a fixed component depending on the type of platform and a variable component based on the flow rate.

CdcwwpNwellwpbswwp,τ

∀ wp, t (51)

(iii) Operating Costs. The total operating cost in each time period is the sum of the operating costs of the gas plant, the PPs, and the WPs [constraint (52)]. For the gas plant, the operating cost is the sum of a

(iv) Taxes, Tariffs, and Royalties. Because of the complexity of the tax, tariff, and royalty specifications, we present a detailed description of the model together with the mathematical representation. The general model is shown for the royalties, but it is similar for the taxes and tariffs. The explicit forms of some functions and variables are not shown for confidentiality reasons. The first specification is that the royalties paid for each WP should be the greater amount as calculated from an objective element 1 (OE1) and the amount as calculated from an objective element 2 (OE2). For ease of presentation, we will specify OE1 to be the net revenue and OE2 the gross revenue, although these are in reality specified by the government involved. OE1 and OE2 can then be calculated as follows:

OE2,wp,t ) Rwp,t

∀ wp, t

OE1,wp,t ) OE2,wp,t g1(CCwp,t;OCwp,t;Taxwp,t;Tariffwp,t)

(56) ∀ wp, t (57)

Constraint (56) sets OE2 equal to the gross revenue for each WP in each time period, while OE1 (in this case the net revenue) is calculated in constraint (57) as the gross revenue minus some amount involving capital costs, operating costs, taxes, and tariffs as represented by g1. If we define Roy1 and Roy2 to be the royalties as calculated from OE1 and OE2 respectively, then

Roywp,t ) max (Roy1wp,t, Roy2wp,t)

∀ wp, t (58)

This constraint can be reformulated as

Roywp,t g Roy1wp,t

∀ wp, t

(58a)

Roywp,t g Roy2wp,t

∀ wp, t

(58b)

Constraints (58a) and (58b) will be enforced as equalities, because royalties are effectively minimized in the objective function (73), which maximizes the negative of the royalties. Roy1 and Roy2 are calculated according to a tier system (tiers are denoted by the index i) as the product of a tier factor and the respective OE. Each OE has its own associated set of tier factors to choose from, namely, f1i and f2i, and the choice is determined by the Tier CEI (Cumulative Economic Indicator) TCEIi,wp,t. The latter variable is calculated in constraint (59) as


Figure 4. Tier structure.

the sum of the CEI, CEIwp,t, and the cumulative economic element, CEEwp,t, multiplied by a tier fraction, fTCEIi. Constraint (60) determines the value of CEIwp,t as the difference of a cumulative function of the capital costs, operating costs, and royalties, g2t, and the cumulative OE2 up to the current period. CEEwp,t equals the cumulative positive CEIwp,t up to the current period as shown in constraint (61).

TCEIi,wp,t ) CEIwp,t + fTCEIiCEEwp,t

∀ i, wp, t (59)

t

CEIwp,t )

t

∑g2τ(CCwp,τ;OCwp,τ;Roywp,t-1) - τ)1 ∑OE2,wp,t τ)1

∀ wp, t (60)

t

CEEwp,t )

CEIwp,t+ ∑ τ)1

∀ wp, t

(61)

For each period, the tier factor from tier i is used to calculate the respective royalties if TCEIi,wp,t is negative, while TCEIi+1,wp,t is still positive. The tier fractions, fTCEIi, are defined in such a way that tier 1 will turn on first, then tier 2, etc., and a tier “turns on” when TCEIi,wp,t becomes negative. This is better explained in Figure 4, which shows typical outcomes of TCEIi as a function of time. Tier 1 is defined where all of the TCEIi are still positive and is the initial case. As time progresses and the cumulative revenues increase, CEI will become negative for certain values of g2, as seen from constraint (60). When CEI becomes negative enough, TCEI will become negative, and the next tier will turn on, namely, tier 2. At this point more royalties will be paid than in tier 1. Similarly, tiers 3 and 4 will turn on as soon as their respective TCEI values turn negative. It is not obvious which of OE1 or OE2 will yield the smallest amount of royalties to be paid, neither is it obvious whether a transition to the next tier will have a positive or negative effect on the overall objective function. There are several tradeoffs due to these rules, and the optimization model presented here can explicitly capture these tradeoffs. To illustrate the complexity of these rules, consider a tier transition over 3 time periods as shown in Table 1 for some fictitious numbers. The first column shows possible values for the variable CEI, while the next three columns show the corresponding values for CEE, TCEI2, and TCEI3 calculated through constraints (59) and (61). If we assume values for OE1, OE2, f1, and f2 as shown, then the

royalties paid will be the maximum of Roy1 and Roy2, as shown in the last column. In the first period, tier 1 is on, while a transition to tier 2 occurs in period 2 when TCEI2 becomes negative, and another transition to tier 3 occurs in period 3 when TCEI3 becomes negative. What is noticeable is the large increase in royalties paid with the transition between tiers 2 and 3. It might therefore be more beneficial for the overall objective to keep OE1, or the net revenue, lower in order to avoid the transition to tier 3. This is counterintuitive, because one would expect to want the highest possible net revenue, but higher net revenues might lead to much larger royalties to be paid, and these might offset the gain in revenues to result in a lower total profit when royalties are subtracted from the net revenues. A simple calculation shows this for the instance in Table 1. Even though period 2 has a lower net revenue of 90, the royalties are also lower, leading to a higher total profit of 83.7 as opposed to the 77.7 obtained in period 3. If the royalty calculations are not included in the optimization model, such tradeoffs would be ignored and the solution might turn out to be severely nonoptimal after the royalties are calculated. However, with the inclusion of royalties in the optimization model, these tradeoffs are considered explicitly as will be shown in example 1, and potentially suboptimal investment and production decisions are avoided. The sequence of calculations is better explained in Figure 5: Values of revenues, capital costs, operating costs, taxes, and tariffs up to the current period and royalties up to the previous period are obtained from the overall model for each WP. The OEs, CEE, and CEI for the current period are then determined as functions of these variables. From these, TCEI is computed in order to define the different tiers. Then, depending on which tier becomes first negative, the royalty contribution is calculated for each time period. The royalties are then the maximum of the contributions as calculated for each economic indicator. We need to emphasize that these rules can be adapted to each specific problem instance. The arrow on the left-hand side of the diagram indicates that the royalty decisions interact with the overall model; e.g., the capital investment might change depending on the royalties that need to be paid, leading to a very complex tradeoff. The conditional part of the model specifies which tier turns on depending on the values of TCEI and requires the use of discrete variables. Initially, this constraint can be written as follows:

If -∞ < TCEIi,wp,t e 0 and 0 e TCEIi+1,wp,t < ∞, then Roy1wp,t ) f1iOE1wp,t and Roy2wp,t ) ∀ i, wp, t (62) f2iOE2wp,t Constraint (62) specifies that if TCEI is negative for tier i but positive for tier i + 1, then the factors associated with tier i are to be used to calculate the royalties. The discrete mathematical formulation of these conditional constraints is not obvious. In a previous paper32 the authors discussed both a disjunctive and big-M formulation for this constraint and concluded that it is not clear which is better. For simplicity, we use a big-M formulation to obtain the results presented in this paper because it is much easier to implement than the disjunctive formulation and performs similar for the model under consideration. The big-M model to calculate Roy1 is as follows (similar for Roy2):

2864


Table 1. Tier Transition over 3 Time Periods period

CEI

CEE

TCEI2 (f ) 0.2)

TCEI3 (f ) 0.4)

OE1

OE2

f1

f2

Roy1

Roy2

royalties paid

1 2 3

140 -50 -150

140 140 140

168 -22 -120

196 6 -94

-50 90 105

100 120 140

0.05 0.07 0.26

0.02 0.03 0.05

-2.5 6.3 27.3

2.0 3.6 7.0

2.0 6.3 27.3

Roy1wp,t )

∑i RoyC1i,wp,t

RoyC1i,wp,t e M1ibToi,wp,t

∀ wp, t

(63)

∀ i, wp, t

(64)

∀ i, wp, t

RoyC1i,wp,t g -M1ibToi,wp,t

(65)

RoyC1i,wp,t e (f1i - f1,i-1)OE1wp,t + M1i(1 - bToi,wp,t) ∀ i, wp, t (66) RoyC1i,wp,t g (f1i - f1,i-1)OE1wp,t - M1i(1 - bToi,wp,t) ∀ i, wp, t (67) Constraint (63) defines the total royalty contributions from OE1 to be the sum of the individual contribution from each tier i. The tier factors are now written as the difference between the current tier factor and the previous tier factor [see constraints (64) and (65)], so as to have the effect of a cumulative contribution with each additional tier that turns on. Constraints (64) and (65) specify that if bToi,wp,t equals 0 for tier i, then the royalty contribution for that tier is set to zero for WP wp in time period t. Similarly, constraints (66) and (67) specify that if bToi,wp,t equals 1 for tier i, then the royalty contribution for that tier is set to the appropriate factor times the OE for WP wp in time period t. Switching between tiers is represented by constraint (68) which forces bToi,wp,t to take a value of 1 if TCEIi,wp,t is negative.

TCEIi,wp,t g -MLTCEIibToi,wp,t

∀ i, wp, t (68)

Finally, some logical constraints are added. Constraint (69) specifies that each tier can be switched on only once, while (70) states that TCEIi,wp,t must be on before TCEIi+1,wp,t can switch on. Constraint (71) specifies that a tier can only be on if it had already switched on, and this constraint forces the continuous variables bToi,wp,t to take only discrete values. Constraint (72) specifies that a tier can only be switched on if the corresponding WP has been installed.

∑t bTsi,wp,t e 1 bToi,wp,t g bTsi+1,wp,t

Figure 5. Royalty calculations.

∀ i, wp ∀ i, wp, t

(69) (70)

Figure 6. 9 WP superstructure for example 1.

same tier structure significantly reduces the size of the model and, hence, the computational effort. (v) Objective. The objective is to maximize the NPV of the project over the complete planning horizon. The NPV consists of the discounted revenues minus the sum of capital costs, operating costs, taxes, tariffs, and royalties, where the terms of Taxwp,t and Tariffwp,t can, in general, also be determined by constraints similar to those for the term Roywp,t.

max NPV ) t

bToi,wp,t )

∑bTsi,wp,τ

∀ i, wp, t

(71)

τ)1 t

bTsi,wp,t e

∑bwp,τsw τ)1

∀ i, wp, t

(72)

It should be noted that often each WP has a different tier structure, meaning that there is a separate set of economic variables and constraints for each WP. However, the proposed model is directly applicable to the case where several WPs have the same tier structure. In fact, grouping more than one WP together with the

∑t dt[Rtott - CCtott - COtott (Taxwp,t + Tariffwp,t + Roywp,t)] ∑ wp

(73)

Example 1 To illustrate the effect that the complex economic rules have on both the economic outcome and computational effort, we compare the solution to the proposed model including royalties, tariffs, and taxes to the solution without considering these in the objective function. The problem under consideration is the design and planning of an infrastructure (see Figure 6) which could potentially consist of up to 9 WPs, denoted by the


letters A-I, and 1 PP, and decisions need to be made for a horizon of 7 years (7 time periods). As shown in Table 2, the inclusion of complex fiscal rules in the optimization model leads to a completely different solution. When royalties, taxes, and tariffs are omitted, the solution is to invest in 5 WPs, namely, A-D and I, and time periods t1, t1, t1, t4, and t4, respectively. In contrast, when these rules are included, the solution is to invest in 6 WPs, namely, A-D, H, and I, in time periods t1, t1, t1, t4, t6, and t5, respectively. Also, the solution without the complex economics requires compression from period 5 onward, while the solution when these factors are considered does not require compression. The solution with complex economics included results in a NPV of $520 million, while the NPV when these rules are not included is $457 million. In the latter case, the NPV is calculated by optimizing the design and planning model without considering the fiscal rules, after which the taxes, royalties, and tariffs are calculated for the optimal design and plan and subtracted. This represents an increase in NPV of $63 million, or 13.8%, when complex economic rules are included. Qualitatively, the large difference in NPV can be explained by the fact that, without complex economic objectives, all decisions are based solely on the discounted costs and revenues. In this specific instance, it seems at first improbable that the additional investment in platform H when the complex objective is used would lead to a higher NPV, because intuitively additional royalties could be expected to only decrease the NPV further from when revenues and operating and investment costs were optimized. However, it turns out that inclusion of complex economics captures the nonobvious tradeoff that when platform H is added, this allows platform I to produce less gas, leading to a large reduction in the royalties that are to be paid for platform I. In other words, the increase in royalties and investment costs due to platform H is outweighed by the large decrease in royalties due to production from platform I. It should be noted that global optimality is not guaranteed for either model, because of the presence of nonconvexities. It is therefore theoretically possible to obtain a worse NPV with complexities than without if a poor suboptimal solution is obtained in the model with complexities and a near-optimal solution is obtained for the model without complexities. However, in the experience of the authors, this never occurred in any of the instances that we solved. Because the nonconvexities in the model are minimal and fairly well-behaved, the comparison of NPVs in Table 2 is justified. A comparison of the problem size is also shown in Table 2. A large increase in the number of constraints, continuous variables, and binary variables results when the fiscal rules are included in the model. When the royalties, taxes, and tariffs are not included, the model solves in 96 CPU s on an HP9000/C110 with the outerapproximation method15 as implemented in DICOPT++ (GAMS33). However, when these rules are included, the solution time increases to 10 530 CPU s on the same machine with the same solver. This is 2 orders of magnitude increase in solution time, which is due to nearly two-thirds of the binary variables belonging to the tax, tariff, and royalty formulation. From these results it is clear that the complex rules are an important issue from both an economical and a computational perspective.

Table 2. Comparison of Solution with and without Complex Economic Considerations without complex objective

with complex objective

WP

invest?

when?

invest?

when?

A B C D E F G H I NPV ($ mil.) constraints continuous vars. binary vars. solution time (CPU s)

yes yes yes yes no no no no yes

t1 t1 t1 t4 n/a n/a n/a n/a t4

yes yes yes yes no no no yes yes

t1 t1 t1 t4 n/a n/a n/a t6 t5

457 2217 1499 76 96

520 3059 1990 142 10530

Solution Method It is clear from example 1 that the solution method calls for some sort of aggregation and/or decomposition technique. The structure of the model suggests a decomposition method where each WP can be solved independently, because the complex economics are allocated to each individual WP and the platforms are only linked by the physical constraints such as mass balances and pressure balances. If these physical constraints can be removed, the model will decompose such that each WP can be solved separately. These thoughts point to the concept of Lagrangean decomposition. Background: Lagrangean Decomposition. Guignard and Kim27 present Lagrangean decomposition as a generalization of Lagrangean relaxation by suitably relaxing side constraints such that the model becomes decomposable over two or more subsets of constraints. Consider a general optimization model (P) of the form

ZP ) max cTx + dT(y1 + y2) s.t.

A1x + B1y1 e b1 A2x + B2y2 e b2 h1(x) e 0 h2(x) e 0 x g 0, y1, y2 ∈ {0, 1}

Model P is from the start decomposable into two subproblems over the binary variables, y1 and y2, but is linked through the continuous variables, x, which occur in all of the constraints. To apply Lagrangean decomposition, model P is reformulated by creating identical copies of the x variables to yield model RP:

ZRP ) max cTx + dT(y1 + y2) s.t.

A1x + B1y1 e b1 A2z + B2y2 e b2 h1(x) e 0 h2(z) e 0 x)z x, z g 0, y1, y2 ∈ {0, 1}

2866


It is clear that model RP will decompose into two separate subproblems if the constraint x ) z is removed. Thus, by dualizing these constraints according to Lagrangean relaxation, we obtain the following decomposable model (LRλ):

ZLRl ) max cTx + dT(y1 + y2) + λT(z - x) s.t.

A1x + B1y1 e b1 A2z + B2y2 e b2 h1(x) e 0 h2(z) e 0

Figure 7. Decoupling of the infrastructure.

x, z g 0, y1, y2 ∈ {0, 1} where λ is the vector of Lagrange multipliers corresponding to the dualized constraints. Model LRλ can now be decomposed into models P1 and P2, and these can be solved independently.

(P1): s.t.

ZP1 ) max cTx + dTy1 - λTx A1x + B1y1 e b1 h1(x) e 0 x g 0, y1 ∈ {0, 1}

and

(P2): s.t.

λk+1 ) λk + tk(zk - xk)

ZP2 ) max dTy2 + λTz 2

2 2

where tk is a scalar step size and zk and xk are optimal solutions to the Lagrangean problem (LRλ) with the multipliers set to λk. The choice of step size, tk, has been shown16 to have a large effect on the convergence to the optimal value of λ. Ideally, tk should converge to zero at a moderate rate and a formula for tk that works well in practice is the following:16

2

A z+B y eb h2(z) e 0

z g 0, y2 ∈ {0, 1} It is well-known17 that ZLRl ) ZP1 + ZP2 is an upper bound to ZP for any λ if the constraints are convex. This follows by assuming an optimal solution x*, y1*, and y2* to P and observing that

Z* ) cTx* + dT(y1* + y2*) e cTx* + dT(y1* + y2*) + λ(z* - x*) ) ZLRl In fact, in the case of equality constraints being dualized, the Lagrangean relaxation should theoretically yield the optimal solution to P for an optimal value of x, y, and λ, because λ(z* - x*) ) 0. Obtaining the tightest upper bound to ZP requires the solution of the Lagrangean dual (D):

ZD ) min ZLRl λ

respectively. Solving D can be difficult to implement and time-consuming, although Fisher17 reports on some algorithms for this purpose. A code for solving the dual was developed by Kiwiel,36 but this code is not widely available. Solving the dual to optimality is therefore often circumvented by using an iterative heuristic approach where LRλ is solved to generate upper bounds to P and a heuristic method is used to generate feasible solutions to P which are also lower bounds. We develop and apply such a heuristic in this work. As for the choice of multipliers, λ, to use in the subproblems, we use an updating scheme based on a subgradient method, because it works well in practice and is relatively easy to implement.17 In this method a sequence of multiplier values are generated by the rule

(D)

If all of the constraints are convex and all of the variables are continuous, the optimum of D will equal the optimum of P. However, a duality gap may exist in the presence of integer variables or other nonconvexities, which means that the optimal solution to the dual problem will be strictly greater than the true optimum of P. Guignard34 and Bazaraa et al.35 give comprehensive graphical interpretations of the duality gap in the case of integer variables and nonconvex constraints,

tk )

Rk(ZLRl(λk) - Z*) |zk - xk|2

In this formula Rk is a scalar, ZLRl(λk) is the optimal solution to the Lagrangean relaxation (LRλ) with multipliers set to λk, and Z* is the best known feasible solution to model P. The initial value of the scalar Rk is chosen between 0 and 2 and in subsequent iterations is halved each time that ZLRl(λk) fails to improve in a fixed number of iterations. Model Reformulation Concept. To apply the theory presented in the previous section, the model needs to be reformulated such that it is easily decomposable on removal of some constraints. The only constraints that link the WPs and PPs together are the mass balances and pressure balances. We cannot simply remove these equations, because this would eliminate too much essential detail at each platform. Instead, we propose duplicating the pipeline variables of flow and pressure, assigning one variable to each WP or PP, and setting them equal. Physically, this can be interpreted as cutting the pipelines in two and thus decoupling the structure (see Figure 7). Consider the pipeline connecting WP1 and WP3 in Figure 7. The pipeline is “cut” at the point where it just


exits WP1, and the pressure at this point, therefore, equals the outlet pressure of WP1, Pout1,t, while the flow rate equals Q1,3,t. Duplicates of these two variables are made, distinguished by the superscript “P”, and equality constraints are added to specify that these are equivalent. By “cutting” the pipeline at this point, the pressure drop equation for the pipeline now belongs to the constraint set for WP3 and will contain the variables Pout,P1,t and QP1,3,t, while the constraint set for WP1 will not contain the pipeline pressure drop equation. The mass and pressure balances at each WP remain intact, while the mass balance in the pipeline itself is taken care of through the equality of the duplicate variables. To further facilitate decomposability, all of the terms in the objective function need to be written separately for each platform and summed. Applying these ideas to the model presented previously yields first a reformulated model (RM), to which Lagrangean relaxation is applied to yield the decomposable model (LRM). Subsequently, LRM is decomposed into a number of Lagrangean subproblems (LRMs), which will be used in the heuristic decomposition algorithm presented next. The detailed reformulation and decomposition to derive RM, LRM, and LRMs is given in Appendix A. Algorithm. The optimal solution to the set of LRMs yields an upper bound to the NPV of the overall model for any choice of Lagrangean multipliers if the constraints are convex. We, therefore, develop the algorithm for the case where all constraints are convex. We develop an iterative heuristic solution procedure where the upper bound is obtained from the Lagrangean decomposition and the lower bound by postulating a feasible solution from the solutions of LRMs. Because the approach is a heuristic regardless of convexity, an optimal solution cannot be guaranteed. In the presence of nonconvex constraints, the proposed method is still applicable, but then upper bounding cannot be guaranteed. The algorithm as applied in this work is shown in Figure 8. Let model RM with the integer requirement on the binary variables relaxed be known as RDM. We obtain the initial multipliers, λ0, by solving RDM. After these initial values are obtained, the sequence of LRMs are solved. The solutions to these subproblems are subsequently used to postulate a feasible set of binary variables for the overall model. After a feasible set of binaries has been found, they are substituted into the original model RM, which is solved as an NLP to obtain a feasible solution. This feasible solution yields a lower bound and is used together with the solutions to LRMs to update the multipliers for the next iteration. Because there are no guarantees for convergence or optimality associated with this method, it is terminated after a prespecified iteration limit. In this work we set the iteration limit to 20 iterations, because this allows one to obtain a good solution in a reasonable amount of time. In practice, one can choose this number according to how much time can be spent on finding the solution depending on the average time per iteration. An additional termination criterion is to terminate after no improvement in the lower bound is obtained in a prespecified number of iterations. Note that the decomposability allows parallelization of the algorithm, although this was not applied in this work. Postulating a feasible set of binary variables from the solutions of the Lagrangean decomposition subproblems

Figure 8. Heuristic Lagrangean decomposition algorithm.

is not an obvious task. The combination of binary variables from the optimal solutions of the original subproblems will typically not yield a feasible solution, and these values need to be manipulated to a point at which a feasible solution results. To find a feasible solution, we first fix all of the binary variables to the values obtained from each subproblem and afterward remove infeasible connections and/or replace them with feasible ones. For example, if from subproblem 1 a WP connects to another WP in subproblem 2 that does not exist, then the connection is removed. Also, if a WP from subproblem 1 connects in time t to a WP from subproblem 2 that only starts operating in a later time t + Tshift, then the connection is replaced with one in time period t + Tshift and the switching of fiscal tiers is shifted with Tshift periods for the first WP. The detailed procedure we propose for obtaining a feasible set of binary variables is shown in Appendix B. After this procedure is implemented for generating solutions, no significant infeasibilities were encountered. The binary variables b10t do not cause infeasibilities when fixed, because the gas flow rate can always be manipulated to reach a maximum in any specified period. Also, the types of WPs do not cause infeasibilities because no additional connections are added in the above procedure. We need to mention that our procedure is problem specific and can easily be changed depending on the specific infrastructure and available knowledge. The relaxed problem RDM and model RM with fixed binaries are NLPs and solved with CONOPT2.37 The Lagrangean subproblems, LRMs, are MINLPs and solved with DICOPT++.15 MILP master problems are solved with CPLEX 6.638 and NLP subproblems with CONOPT2. All models are written in GAMS and run on either a HP9000/C110 workstation (examples 1 and 2) or a Pentium III 667 MHz machine (example 3). It should be noted that some convex approximations of the nonlinear constraints (2), (3), and (19) for the MILP master problems were attempted initially in order

2868


Figure 9. 9 WP superstructure before and after decomposition. Table 3. Example 1: Comparison between Models with and without Complex Objectives and the Proposed Method with complex objective without complex objective

full space

proposed algorithm

96 457.5 0

10530 520 13.8

676 514.6 12.5

solution time (CPU s) NPV ($ mil.) improvement (%)

to guarantee a global optimum for the MINLP subproblems, for example, convex envelopes14 and McCormick underestimators,39 but these estimations led to infeasible solutions in the corresponding NLP. It turns out that, by using the straightforward linearization technique in DICOPT++, infeasiblilities in the NLP part of the outer approximation algorithm did not occur frequently, and good solutions were obtained. Example 1 Revisited To demonstrate the effectiveness of the proposed method, we revisit example 1 and compare the NPV with the NPV obtained without complex objectives included. We also compare the computational time used by the proposed algorithm to that required for solving the model in full space. For the purpose of decomposition, it is useful to note that 6 of the WPs, A-F, have the same tier structure and can therefore be grouped together with the PP into one subproblem. The remaining WPs, G-I, are each represented in separate subproblems, giving a total of 4 subproblems (see Figure 9). Each subproblem is a multiperiod MINLP to be solved over a horizon of 7 time periods. From Table 3 it is clear that the proposed method shows significant improvements in solution time compared to the full-space solution. While the full-space model takes 10 530 CPU s to solve, the solution with the proposed method requires only 676 CPU s. This is more than an order of magnitude reduction in solution time. Even though the proposed algorithm does not find the optimal solution of $520 million for the NPV, it finds an NPV of $514.6 million after 20 iterations, much better than if complex economics are not included. This is an improvement in NPV of 12.5%, or $57.1 million. Table 4 shows the progression of bounds reported for every second iteration in the proposed Lagrangean decomposition algorithm. The algorithm is initialized

Table 4. Example 1: Sequence of Bounds from the Lagrangean Decomposition Algorithm iteration

Lagrangean soln ($ mil.)

NLP relaxation 2 4 6 8 10 12 14 16 18 20

872.3 651.8 592.7 692.7 615.5 640.8 615.1 612.6 615.6 619.1 614.5

postulated LB ($ mil.) * 475.6 497.7 475.6 510.2 476.2 514.6 482.0 514.6 475.5

with the NLP relaxation which yields an initial upper bound of $872.3 million in 28 CPU s. The next step in the algorithm is to solve the sequence of Lagrangean subproblems with multipliers set to those from the NLP relaxation, yielding a solution of $651.8 million. Note the much better quality of the solution obtained from the Lagrangean decomposition compared to the NLP relaxation. Next, a solution is postulated, which turns out to be infeasible for this first iteration (indicated by an asterisk). The algorithm continues to yield a lower bound of $475.6 million in the next iteration, and so forth. The best NPV of $514.6 million is obtained after 13 iterations, and the final gap between bounds for the proposed method is 13.2%. Table 5 shows the typical problem size for the NLP relaxation, the Lagrangean decomposition subproblems, and the lower bounding NLP in terms of the number of variables and constraints, as well as the average solution times per iteration. The number of binary variables reported for the NLP relaxation represents the actual number of binaries for the full-space model, even though they are relaxed here. For the lower bounding NLP, the size is reported for the best solution. The average time per iteration for solving one sequence of Lagrangean MINLP subproblems is 32 CPU s, and the average time to solve the postulated lower bounding NLP is 2.5 CPU s. The actual solution obtained by the proposed method is to invest in WPs A-D, H, and I in time periods t1, t1, t1, t4, t6, and t5, respectively. This is the same as that for the full-space solution of $520 million. The main cause for the discrepancy with the full-space solution is the timing of the switching of economic tiers for WP I. In the optimal solution, tiers 2 and 3 are reached in

Ind. Eng. Chem. Res., Vol. 40, No. 13, 2001 2869 Table 5. Example 1: Problem Size and Solution Times problem size model

WPs

constraints

continuous vars.

binary vars.

solution time (CPU s)

NLP relaxation Lagrangean sub1 Lagrangean sub2 Lagrangean sub3 Lagrangean sub4 Lower bounding NLP

9 6 1 1 1 b

2668 1456 583 661 571 2514

1875 1055 420 440 412 1733

142a 70 32 36 28 n/a

28 16 4 6 6 2.5

a

Number of binary variables for the full-space model. b The number of WPs depends on postulation.

solution to even the first MILP in the outer approximation algorithm used in DICOPT++. With the proposed algorithm, a good solution is obtained in just more than 1 h. The NPV from the proposed method, namely, $933.8 million, is $55.4 million higher than the NPV obtained without the complex objectives, and this represents an increase of 6.3%. As in the previous example, the NPV for the problem without the complex objective is obtained by optimizing the design and planning problem without including the fiscal rules, calculating the taxes, royalties, and tariffs for this optimal design and plan, and subtracting it to obtain the actual NPV. Example 3

Figure 10. Comparison of gas production for WP I from the proposed method and full-space solution, respectively. Table 6. Comparison of the Results for 11 Years with complex objective

solution time (CPU s) NPV ($ mil.)

without complex objective

full space

proposed algorithm

748 878.4

>100 000 n/a

4587 933.8

t6 and tier 4 in t7, while for the solution of the proposed method, tiers 2 and 3 are reached in t7 and tier 4 is never reached. This can be explained by looking at the gas production for WP I as shown in Figure 10. The lower production from the proposed method has the result that the switching of tiers is postponed. This has the effect that less taxes, royalties, and tariffs need to be paid. However, in this instance the resulting savings are offset by the loss of revenue from the lower production, leading to the lower NPV of $514.6 million. Example 2 To obtain the best NPV of $520 million in example 1, the model has to be solved in the full space and the associated solution time of 10 530 CPU s might be considered reasonable for a $6 million improvement over the NPV from the proposed method. However, when longer time horizons than 7 years are considered or when more WPs are added, the problem becomes intractable and it is especially here where the Lagrangean decomposition pays off. Consider example 1, but with the planning horizon extended to 11 years instead of 7. Table 6 shows the comparison of solution times and NPV. The model without complex objectives included solves in 748 CPU s, but when the taxes, royalties, and tariffs are added, more than 3 days pass without getting a

In this final example we consider the complete oilfield infrastructure (refer to Figure 2) consisting of 16 WPs (A-P) and 1 PP for a planning horizon of 15 years (15 time periods). We solve this example with the proposed method and compare the solution to that obtained without the inclusion of complex objectives. No comparison with the full-space model is shown, because it is intractable and no solution to even one MILP could be obtained after 5 days. However, the problem size for the full-space model is shown in Table 7, together with the problem sizes of the various subproblems. The times reported for the proposed method correspond to the iteration that yielded the best solution. For subproblems 7-11, the solution times are so short (∼5 CPU s) because the associated WPs can only be installed after year 9, and these subproblems are thus effectively solved for 5 periods instead of 15. Table 8 shows a comparison of the NPV and total solution time. The model without complex objectives yields a NPV of $1123.4 million, and the solution time is 26.8 CPU h. The proposed method yields a higher NPV of $1218.7 million, which represents an increase in NPV of 8.5% or $95.3 million, in 12.3 CPU h or 20 iterations. Table 9 shows the progression of bounds reported for every second iteration in the proposed Lagrangean decomposition algorithm. Note the convergence of the Lagrangean solutions to the final value of $1242.8 million corresponding to the final lower bound of $1218.7 million. The average time per iteration for solving one sequence of Lagrangean MINLP subproblems is 2204 CPU s, and the average time to solve the lower bounding NLP is 69 CPU s. The final selected infrastructure is shown in Figure 11. All of the platforms, except platforms J-L and O, are selected. For the platforms where a choice of connection was available, P connects to the PP and not to E, N and M connect to the PP and not to H, H connects to the PP and not to B, and G connects to H and not to B. The comparison of timing of investments, WP types, connections, and compression requirements is shown in

2870


Table 7. Example 3: Comparison of Problem Sizes problem size model

WPs

constraints

continuous vars.

binary vars.

solution time (CPU min)

without complex objective without complex objective (proposed) NLP relaxation Lagrangean sub1 Lagrangean sub2 Lagrangean sub3 Lagrangean sub4 Lagrangean sub5 Lagrangean sub6 Lagrangean sub7 Lagrangean sub8 Lagrangean sub9 Lagrangean sub10 Lagrangean sub11 lower bounding NLP

16

8385

5676

341

1607.8

16 6 1 1 1 1 1 1 1 1 1 1 b

12696 4291 1503 1424 1647 1391 1467 955 955 940 955 955 11951

8552 2999 1044 999 1154 977 1020 729 729 719 729 729 7793

759a 269 88 82 110 71 76 46 46 39 46 46 n/a

a

15.3 7.0 3.0 5.2 35.6 1.0 11.7 0.08 0.03 0.03 0.05 0.02 1.7

Number of binary variables for the full-space model. b The number of WPs depends on the postulation. Table 10. Comparison of Investment Timing for 16 WPs over 15 Years without complex objectives equipment

Figure 11. Infrastructure obtained with the proposed method. Table 8. Comparison of the Results for 16 WPs over 15 Years with complex objectives

solution time (CPU h) NPV ($ mil.)

without complex objectives

full space

proposed algorithm

26.8 1123.4

n/a n/a

12.3 1218.7

Table 9. Sequence of Bounds from the Lagrangean Decomposition Algorithm iteration

Lagrangean soln ($ mil.)

postulated LB ($ mil.)

NLP relaxation 2 4 6 8 10 12 14 16 18 20

2271.0 1611.8 1363.7 1295.7 1265.2 1254.8 1248.6 1246.2 1243.5 1243.3 1242.8

* 1164.4 * 1174.4 1185.0 1205.7 1217.2 1107.5 1218.7 1117.7

Table 10. In the case where complex objectives are not included, only 10 WPs are selected, namely, A-I and P, and both stages of the compressor are required from year 6 onward. In the case where complex objectives are included, 12 WPs are selected and compressor stage 1 switches on in year 9, while compressor stage 2 switches on in year 11, meaning that the compressor installation is not required until year 9.

PP compressor stage 1 compressor stage 2 WP A B C D E F G H I J K L M N O P

proposed method

time connected time connected period type to period type to 1 6

shore PP

6

PP

11 1 1 1 8 8 4 12 6 6

I I I I I III II I III

PP PP PP PP PP B H PP C

12

I

PP

1 1 1 7 11 4 14 9 4

I I I I I III III II III

PP PP PP PP PP B B B C

11 13

I I

PP PP

12

I

PP

1 9

shore PP PP

Figures 12 and 13 show the production schedule with and without the complex fiscal rules. The bar diagram shows how much each WP contributes to the total production in each time period. Finally, we need to emphasize that the proposed model and algorithm are not intended for one-time use in order to plan for the next 15 years. Instead, the model will be updated and resolved periodically as more information becomes available or as the effect of implemented decisions is seen. Remark. The model discussed in this work is deterministic in that all of the data are considered to be perfectly known. This means that uncertainties in, for example, future gas prices, demand, and reservoir characteristics are not considered explicitly. However, it is important to consider these uncertainties when making long-term decisions. To illustrate the point, we solve the model in example 3 with a 10% increase and decrease respectively in the gas selling price from year 7 onward to see how sensitive the solution is to small changes. Indeed, one way of experimenting with uncertainty in practice is to run multiple scenarios of the deterministic model with different data sets.


Figure 12. Production schedule from the proposed method.

Figure 13. Production schedule without complex objectives.

The results of our experiment are shown in Table 11. The WPs that are invested in are the same as those in the base case. However, the timing of the investment differs for platforms G-I, M, and N for a 10% increase in price, while the timing of the investment differs for platforms D, G, H, N, and P for a 10% decrease in price. Also, the type of WP H is different for the 10% decrease in price. These differences are for the most part not large, with the greatest difference in timing being for platform H in the case of a 10% decrease in price. However, gas prices can fluctuate with 50% or more, and significantly different solutions can be expected in

such cases. There is, therefore, a strong need for the development of methods for explicitly incorporating issues of uncertainty into the long-term planning of offshore hydrocarbon field infrastructures. The biggest challenge in achieving such a goal lies in dealing with the problem size, because one can expect another order of magnitude increase in solution time. Conclusions We have presented a model for the long-term design and planning of offshore hydrocarbon field infrastruc-

2872


Table 11. Effect of Gas Price Uncertainty on Solution 10% increase in gas price equipment

time connected time connected period type to period type to

PP compressor stage 1 compressor stage 2 WP A B C D E F G H I J K L M N O P

1 9

shore PP

11

1 9

PP

10

shore PP

1 1 1 8 8 4 13 5 5

I I I I I III II I III

PP PP PP PP PP B H PP C

1 1 1 7 8 4 11 11 6

I I I I I III II II III

PP PP PP PP PP B H B C

12 14

I I

PP PP

11 14

I I

PP PP

12

I

PP

11

I

PP

Acknowledgment The authors acknowledge financial support from Mobil Technology Co. and the Department of Energy under Grant DE-FG02-98ER14875. The authors also acknowledge valuable comments from Prof. Arthur Westerberg. Appendix A: Model Reformulation and Decomposition The reformulated model RM is as follows:

(Rwp,t - CCwp,t - CCpipewp,t ∑t dt[∑ wp

CCdevwp,t - OCwp,t - Taxwp,t - Tariffwp,t - Roywp,t) -

(CCpp,t + CCcmppp,t + OCpp,t)] ∑ pp

subject to (1)-(4), (7)-(10), (14)-(51), and (56)-(77) Qwp,t +

QPwp′,wp,t ) ∑Qwp,wp′,t + ∑Qwp,pp,t ∑ wp′ wp′ pp

∀ wp, t (5a)

QPwp,pp,t ∑ wp

∀ pp, t

(6a)

Pinwp,t ) Pout,Pwp′,t - Pdropwp′,wp,tQPwp′,wp,tDistwp′,wp,t Pchkwp′,wp,t + Pslkwp′,wp,t

∀ wp, wp′, t (11a)

Pinpp,t ) Pout,Pwp,t - Pdropwp,pp,tQPwp,pp,tDistwp,pp,t -

PP

tures with complex fiscal rules regarding royalties, taxes, and tariffs. We have shown that the inclusion of these complex objectives into the optimization model can lead to a substantial increase in the NPV of the project. It was also demonstrated that this inclusion leads to more than an order of magnitude increase in the computational effort required for solving the problem. To address the computational difficulties, we have proposed a heuristic Lagrangean decomposition algorithm. The results show that application of the proposed method leads to more than an order of magnitude decrease in solution time and, in the case of intractable instances, yields good solutions in relatively short times. It has also been shown that the NPV obtained by the proposed model and algorithm is significantly higher than the NPV if complex objectives are not included.

max NPV )

Qpp,t )

10% decrease in gas price

∀ pp, wp, t (12a)

Pchkwp,pp,t + Pslkwp,pp,t Poutwp,t ) Pinwp,t - Pdropwp,t(Qwp,t +

QPwp,wp′,t) ∑ wp′

∀ wp, t (13a)

∀ (.), t

(21a)

1 FOcstgpt + FOcstpp,tbonpp,t |PP|

∀ pp, t (54a)

QP(.),t e M(.)bon(.),t OCpp,t ) OCwp,t )

FOcstwp,tpbonwp,tp,t + VOcstwpQwp,t + ∑ tp VOcstgpQwp,t

Q(.),t ) QP(.),t

∀ wp, t (55a)

∀ (.) ∈ {(wp, pp); (wp, wp′)}, t (79)

Poutwp,t ) Pout,Pwp,t

∀ wp, t

(80)

The objective function is now completely separable for each WP and PP. The reservoir constraints (1)-(4) are already separated and do not need to change, while the compressor remains connected to its associated PP and these constraints remain the same. Pressure balances between the reservoir and WP, at the PPs and compressors, and at the shore remain the same because they can be allocated to either a WP or PP respectively [constraints (10) and (14)-(20)]. The big-M constraints (21)-(24), as well as the logical conditions (25)-(42), remain unchanged. Note that, in the case where a linking binary variable, such as bswwp,wp′,t, occurs in a logical constraint, the variable as well as the constraint is duplicated for both wp and wp′. The binary variable is not required to be equal because equality will theoretically be enforced by the big-M constraints (21) and (21a) together with the equivalence of the continuous flow variables. Flow constraints (43) and (44), as well as revenue calculation (45) and capital cost calculations (46)-(51), remain in tact. The capital cost for only the incoming pipelines is allocated to platforms. The model for taxes, tariffs, and royalties, represented by constraints (56)-(77), is completely decomposable and therefore remains as is. Constraints (5a) and (6a) are the modified versions of constraints (5) and (6). Note that the duplicate variables, QP(.),t, are used for the incoming flows instead of the original variables. Constraints (11a), (12a), and (13a) correspond to the original constraints (11)-(13) where the duplicate variables, QP(.),t and Pout,Pwp,t, are used for all incoming flows and pipeline sections. Constraint (21a) complements (21) by setting the duplicate variables to zero if a connection is off. Constraints (52)-(55), the operating cost calculations, are removed and replaced with constraints (54a) and (55a), where the fixed operating cost of the gas plant is distributed among the PPs and the variable operating cost is distributed among the WPs. This


distribution captures the tradeoff between revenues and costs for each individual platform more accurately than had the gas-plant costs been considered as a whole, while the total still equals the total operating cost of the gas plant as specified originally in (53). The additional constraints (79) and (80) are the socalled “linking” constraints that link the duplicate variables to their original counterparts. These constraints link the platforms together, and removal of them yields a completely decoupled structure. Model RM has the same form as model RP and Lagrangean relaxation can now be applied to it to yield model LRM:

max NPV )

(Rwp,t - CCwp,t - CCpipewp,t ∑t dt[∑ wp

CCdevwp,t - OCwp,t - Taxwp,t - Tariffwp,t - Roywp,t) -

(CCpp,t + CCcmppp,t + OCpp,t) + ∑ pp ∑ ∑λwp,wp′,t(QPwp,wp′,t - Qwp,wp′,t) + wp wp′ λwp,pp,t(QPwp,pp,t - Qwp,pp,t) + ∑ ∑ wp pp λwp,t(Pout,Pwp,t - Poutwp,t)] ∑ wp

subject to (1)-(4), (7)-(10), (14)-(51), and (56)-(77) (5a), (6a), (11a)-(13a), (21a), (54a), and (55a) where the λ’s are the Lagrangean multipliers corresponding to constraints (79) and (80). Model LRM is decomposable into up to N subproblems, where N equals the number of WPs plus the number of PPs. In general, each PP can be grouped in the same subproblem with one or more WPs, because there is no revenue allocated specifically to the PPs. This will prevent solutions where the PPs do not get chosen. In addition, more than one WP can be grouped into the same subproblem, depending on the time it takes to solve and the quality of the upper bound obtained. To decompose (LRM), we define the set SUB(s,wp,pp) to be the set of WPs and PPs belonging to subproblem s. The set of subproblems after decomposition is as follows (LRMs):

max NPVs )

∑t {dt[

CCpipewp,t -

∑

(Rwp,t - CCwp,t wp∈SUB(s,wp) CCdevwp,t - OCwp,t - Taxwp,t -

∑

Tariffwp,t - Roywp,t) -

∑

OCpp,t)} +

(CCpp,t + CCpp,tcmp +

pp∈SUB(s,pp)

∑

λwp,wp′,tQPwp,wp′,t -

wp∉SUB(s,wp)wp′∈SUB(s,wp′)

∑

∑

λwp,wp′,tQwp,wp′,t +

∑

∑

λwp,pp,tQPwp,pp,t -

∑

∑

wp′∉SUB(s,wp′)wp∈SUB(s,wp)

wp∉SUB(s,wp)pp∈SUB(s,pp)

λwp,pp,tQwp,pp,t + pp∉SUB(s,pp)wp∈SUB(s,wp) λwp,tPout,Pwp,t λwp,tPoutwp,t] wp∉SUB(s,wp) wp∈SUB(s,wp)

∑

∑

subject to (1)-(4), (7)-(10), (14)-(51), and (56)-(77) (5a), (6a), (11a)-(13a), (21a), (54a), and (55a) wp, pp ∈ SUB(s,wp,pp)

}

In the objective function only costs and revenues of WP’s and PP’s belonging to that subproblem are taken into account. In the additional Lagrangean terms, the duplicate variables (with superscript “P”) are taken into account only if they belong to a pipeline coming into a WP or PP belonging to subproblem s from a WP not in subproblem s. The duplicate pressure variable is only added if the associated WP does not belong to s but is attached by a pipeline to a WP that does belong to s. The original variables are added if the WP belongs to s and is connected by a pipeline to a WP or PP not belonging to s. All of the constraints are applied only if the WP or PP belongs to subproblem s as determined by the set SUB(s,wp,pp). Appendix B: Procedure for Postulating Feasible Solutions The procedure for postulating feasible solutions is as follows: Step 1. From the solutions (LRMs), s ) 1, ..., S, fix the binaries bswwp,t, bswwp,wp′,t, bswwp,pp,t, btypewp,tp, bswpp,t, bswpp,cmp,t, b10t, and bTsi,wp,t for all wp, pp ∈ SUB(s,wp,pp). Step 2. Eliminate infeasible existence of WPs. If ∑tbswwp,t ) 1, bswwp,wp′,t ) 1, and ∑tbswwp′,t ) 0, i.e., if wp exists and connects to wp′ but wp′ does not exist, then remove wp from the set of feasible WPs by setting bswwp,t, bswwp,wp′,t, btypewp,tp, bTsi,wp,t ) 0 for all time periods t. It is similar for connections to PPs. Step 3. Correct infeasible timing of connections. If bswwp,t ) 1, bswwp,wp′,t ) 1, and bswwp′,τ ) 1 for some τ > t, i.e., if wp switches on and connects to wp′ in period t while wp′ switches on only in a later period τ, then perform the following: (a) Set bswwp,t and bswwp,wp′,t ) 0 and bswwp,τ and sw b wp,wp′,τ ) 1, to allow wp to be installed in the same period as wp′ and not before. (b) Let Tshift ) |τ| - |t|. If an economic tier for wp switches on in period t′ (bTsi,wp,t′ ) 1), then it will have to switch on in a later period t′ + Tshift instead. Set bTsi,wp,t′ ) 0 and bTsi,wp,t′+Tshift ) 1. If t′ + Tshift > |T|, then set bTsi,wp,t ) 0 for all periods t. It is similar for connections to PPs. Nomenclature Sets and Indices T ) set of time periods t ) time period t ∈ T PP ) set of PPs pp ) production platform pp ∈ PP F ) set of fields f ) field f ∈ F WP ) set of WPs wp ) well platform wp ∈ WP (alias wp′) C ) set of compressor stages cmp ) compressor stage cmp ∈ C TP ) set of WP types tp ) WP type tp ∈ TP Tinv(t) ) set of time periods in investment horizon preceding period t I ) set of economic tiers

2874


i ) economic tier i ∈ I S ) Lagrangean relaxation subproblems s ) subproblem s ∈ S K ) set of iterations for heuristic procedure k ) iteration k∈ K SUB(s,wp,pp) ) set of wp’s and pp’s belonging to subproblem s Note: The indices f and wp can be used interchangeable since there is one WP per field.

Nomenclature Subscript Combinations Denoted by (.) (f) or (wp) ) variable associated with field f or wp (pp) ) variable associated with production platform pp (wp, wp′) ) variable associated with the pipe section between wp and wp′ (wp, pp) ) variable associated with the pipe section between wp and pp (pp, cmp) ) variable associated with compressor stage cmp at pp (wp, tp) ) variable associated with wp type tp Continuous Variables Q(.),t ) gas/oil flow rate in period t QP(.),t ) duplicate of Q(.),t for pipeline decoupling Qdf,t ) gas/oil deliverability of field f in period t Qcf,t ) cumulative gas/oil production of field f in period t Qbypp,cmp,t ) gas/oil flow rate bypassing compressor stage cmp at pp in period t Qshrt ) gas/oil flow rate at the shore in period t Psf,t ) surface pressure at field f in period t Prf,t ) reservoir pressure at field f in period t Pbhf.t ) bottom hole pressure at field f in period t Pin(.),t ) inlet pressure in period t Pout(.),t ) outlet pressure in period t Pout,P(.),t ) duplicate of Pout(.),t for pipeline decoupling Pchk(.),t ) choke pressure in pipeline in period t Pslk(.),t ) slack to enforce pressure balance or not in period t Pbstpp,cmp,t ) pressure boost for compressor stage cmp at pp in period t Pshrt ) pressure at shore in period t Fuelpp,cmp,t ) fuel required by compressor stage cmp at pp in period t CCtott ) total capital cost in period t CCwp,t ) capital cost of wp in period t CCpp,t ) capital cost of pp in period t CCpipewp,t ) capital cost of pipeline connection from wp in period t CCcmppp,t ) capital cost of compressor at pp in period t CCdevwp,t ) capital cost for well development at WP wp in period t OCtott ) total operating cost in period t OCwp,t ) operating cost for each WP and its associated pipeline in period t OCpp,t ) operating cost for each PP in period t OCgpt ) operating cost for the gas plant in period t Rtott ) total revenue in period t Rwp,t ) revenue from wp in period t Roywp,t ) royalties for wp in period t OE1,wp,t ) objective element 1 for wp in period t OE2,wp,t ) objective element 2 for wp in period t Roy1wp,t ) royalties as calculated from OE1 for wp in period t Roy2wp,t ) royalties as calculated from OE2 for wp in period t Saleswp,t ) sales revenue for wp in period t Taxwp,t ) taxes paid for wp in period t Tariffwp,t ) tariffs paid for wp in period t CEIwp,t ) cumulative economic indicator for wp in period t

TCEIi,wp,t ) tier i CEI for wp in period t CEEwp,t ) cumulative economic element for wp in period t RoyC1i,wp,t ) royalty contribution from OE1 when TCEIi,wp,t e 0 and TCEIi+1,wp,t g 0 Binary Variables bsw(.),t ) 1 if pp/wp/pipeline/compressor stage is installed and switched on in period t; 0 otherwise btypewp,tp ) 1 if wp is of type tp; 0 otherwise b10t ) 1 if a new maximum flow rate to the shore starts in period t; 0 otherwise bTsi,wp,t ) 1 if TCEIi,wp,t switches e 0 in t; 0 otherwise Binary Variables That Can Be Treated as Continuous bswwp,tp,t ) 1 if wp type tp is installed in period t; 0 otherwise bon(.),t ) 1 if pp/wp/pipeline/compressor stage is in operation in period t; 0 otherwise bToi,wp,t ) 1 if TCEIi,wp,t e 0 in t; 0 otherwise Parameters Pinitf ) initial pressure for field f Pdrop(.) ) expected pressure drop in a pipeline or over the structure Dist(.) ) distance between structures shrink ) shrinkage factor for flow in a pipeline toward the shore CRcmp ) compressor ratio for compressor stage cmp f conv ) conversion factor to ensure proper units for compressor calculation Zavg ) compressibility at average temperature and pressure ηadia ) adiabatic compressor efficiency k ) specific heat capacity Tin ) inlet temperature at compressor M ) upper bound on oil/gas flow rates Mslk ) big-M parameter for enforcing pressure constraints Mbstcmp ) upper bound on pressure boost in compressor stage cmp Ccst(.) ) capital cost Cdcwwp ) capital cost for drilling and completing a well at wp Nwellwp ) estimated number of wells to be drilled at WP wp fcc(.) ) fraction of capital cost occurring in each year of the investment horizon preceding installation Ccon(.) ) capital cost for a pipeline connection Cmile(.) ) per mile capital cost for a pipeline FOcst(.) ) fixed operating cost VOcst(.) ) variable operating cost Pt ) predicted price of oil/gas in period t ∆t ) number of days in time period t Qshrmax ) maximum possible flow rate to the shore f1,i ) royalty factor for OE1 in tier i f2,i ) royalty factor for OE2 in tier i fTCEIi ) factor determining TCEIi,wp,t dt ) discounting factor in period t M1i ) upper bound on the royalty contribution from OE1 in tier i MLTCEIi ) lower bound on TCEIi,wp,t MUTCEIi ) upper bound on TCEIi,wp,t Tshift ) shifting time for postulation of a feasible solution in Lagrangean heuristic

Literature Cited (1) ExxonMobil Financial and Operating Review, 1999. (2) Bohannon, J. A. Linear Programming Model for Optimum Development of Multi-Reservoir Pipeline Systems. J. Pet. Technol. 1970, 22, 1429. (3) Haugland, D.; Hallefjord, Å.; Asheim, H. Models for Petroleum Field Exploitation. Eur. J. Oper. Res. 1988, 37, 58. (4) Iyer, R. R.; Grossmann, I. E.; Vasantharajan, S.; Cullick, A. S. Optimal Planning and Scheduling of Offshore Oil Field

Ind. Eng. Chem. Res., Vol. 40, No. 13, 2001 2875 Infrastructure Investment and Operations. Ind. Eng. Chem. Res. 1998, 37, 1380. (5) van den Heever, S. A.; Grossmann, I. E. An Iterative Aggregation/Disaggregation Approach for the Solution of a Mixed Integer Nonlinear Oilfield Infrastructure Planning Model. Ind. Eng. Chem. Res. 2000, 39 (6), 1955. (6) Aranofsky, J. S.; Williams, A. C. The Use of Linear Programming and Mathematical Models in Underground Oil Production. Manage. Sci. 1962, 8, 394. (7) Attra, H. D.; Wise, W. B.; Black, W. M. Application of Optimizing Techniques for Studying Field Production Operations. 34th Annu. Fall Meeting SPE 1961, 82. (8) Frair, L. C. Economic Optimization of Offshore Oilfield Development. Ph.D. Dissertation, University of Oklahoma, Tulsa, OK, 1973. (9) Lee, A. S.; Aranofsky, J. S. A Linear Programming Model for Scheduling Crude Oil Production. Pet. Trans., AIME 1958, 213, 389. (10) Eeg, O. S.; Herring, T. Combining Linear Programming and Reservoir Simulation to Optimize Asset Value. SPE 37446 Presented at the SPE Production Operations Symposium, Oklahoma City, OK, Mar 9-11, 1997. (11) Sullivan, J. A Computer Model for Planning the Development of an Offshore Oil Field. J. Pet. Technol. 1982, 34, 1555. (12) Iyer, R. R.; Grossmann, I. E. A Bilevel Decomposition Algorithm for Long-Range Planning of Process Networks. Ind. Eng. Chem. Res. 1998, 37, 474. (13) Turkay, M.; Grossmann, I. E. Logic-Based MINLP Algorithms for the Optimal Synthesis of Process Networks. Comput. Chem. Eng. 1996, 20, 959. (14) Horst, R., Pardalos, P. M., Eds. Handbook of Global Optimization; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1995. (15) Viswanathan, J.; Grossmann, I. E. A combined penalty function and outer-approximation method for MINLP optimization. Comput. Chem. Eng. 1990, 14, 769. (16) Fisher, M. L. An Applications Oriented Guide to Lagrangean Relaxation. Interfaces 1985, 15 (2), 10. (17) Fisher, M. L. The Lagrangean Relaxation Method for solving integer programming problems. Manage. Sci. 1981, 27 (1), 1. (18) Holmberg, K.; Yuan, D. A Lagrangean Heuristic Based Branch-and-Bound Approach for the Capacitated Network Design Problem, Oper. Res. 2000, 48 (3). (19) Adhya, N.; Tawarmalani, M.; Sahinidis, N. V. A Lagrangean Approach to the Pooling Problem. Ind. Eng. Chem. Res. 1999, 38, 1956. (20) Gupta, A.; Maranas, C. D. A Hierarchical Lagrangean Relaxation Procedure for Solving Midterm Planning Problems. Ind. Eng. Chem. Res. 1999, 38, 1937. (21) Lim, S.; Kim, Y. Capacity planning for phased implementation of flexible manufacturing systems under budget restrictions. Eur. J. Oper. Res. 1998, 104, 175. (22) Carøe, C. C.; Schultz, R. Dual Decomposition in Stochastic Integer Programming. Oper. Res. Lett. 1999, 24 (1-2), 37. (23) Nowak, M. P.; Ro¨misch, W. Stochastic Lagrangean Relaxation applied to Power Scheduling in a Hydro-Thermal System

under Uncertainty, Preprint 98-24 of Institut fu¨r Mathematik an der Mathematisch-Naturwissenschaftlichen Fakulta¨t II der Humboldt-Universita¨t zu Berlin, 1998. (24) Feltenmark, S. On Optimization of Power Production. Doctoral Thesis, Department of Mathematics/Optimization and Systems Theory, Royal Institute of Technology, Sweden, 1997. (25) Hoitomt, D. J.; Luh, P. B.; Pattipati, K. R. A Practical Approach to Job-Shop Scheduling Problems. IEEE Trans. Rob. Autom. 1993, 9 (1). (26) Chang, S.; Liao, D. Scheduling Flixible Flow shops with no setup effects. IEEE Trans. Rob. Autom. 1994, 10 (2). (27) Guignard, M.; Kim, S. Lagrangean decomposition: A model yielding stronger Lagrangean bounds. Math. Program. 1987, 39, 215. (28) Jo¨rsten, K. O.; Na¨sberg, M.; Smeds, P. A. Variable SplittingsA New Lagrangean Relaxation Approach to some Mathematical Programming Models. Department of Mathematics Report LiTH-MAT-R-85-04; Linko¨ping Institute of Technology: 1985. (29) Marı´n, A.; Pelegrı´n, B. The return plant location problem: Modeling and resolution. Eur. J. Oper. Res. 1998, 104, 375. (30) Equi, L.; Giorgio, G.; Marziale, S.; Weintraub, A. A combined transportation and scheduling problem. Eur. J. Oper. Res. 1997, 97, 94. (31) Rana, K.; Vickson, R. G. Routing Container Ships Using Lagrangean Relaxation and Decomposition. Transp. Sci. 1991, 25 (3). (32) van den Heever, S. A.; Grossmann, I. E.; Vasantharajan, S.; Edwards, K. Integrating complex economic objective with the design and planning of offshore oilfield infrastructures. Comput. Chem. Eng. 2000, 24, 1049. (33) Brooke, A.; Kendrick, D.; Meeraus, A. GAMS: A User’s Guide, Release 2.25; The Scientific Press: South San Francisco, 1992. (34) Guignard, M. Lagrangean relaxation: A short course. Belg. J. OR: Special Issue Francoro 1995, 35, 3. (35) Bazaraa, M. S.; Sherali, H. D.; Shetty, C. M. Nonlinear Programming; John Wiley: New York, 1994. (36) Kiwiel, K. C. User’s Guide for NOA 2.0/3.0: A Fortran package for convex nondifferentiable optimization; Polish Academy of Sciences, Systems Research Institute: Warsaw, Poland, 1993/ 1994. (37) Drud, A. CONOPTsA Large Scale GRG Code. ORSA J. Comput. 1994, 6, 207. (38) ILOG CPLEX 6.6. User’s Manual; ILOG Inc.: 2000. (39) McCormick, G. P. Computability of global solutions to factorable nonconvex programs: Part IsConvex underestimating problems. Math. Program. 1976, 10, 147.

Received for review August 15, 2000 Revised manuscript received February 7, 2001 Accepted April 25, 2001 IE000755E

A Lagrangean Decomposition Heuristic for the Design and Planning of

Recommend Documents