Scheduling Gasoline Blending Operations from ... - ACS Publications

[AIChE J.2010, 56, 441–465] to avoid solving mixed-integer nonlinear ... While they considered transitions in blending pipelines, they did not ensur...
0 downloads 0 Views 2MB Size
ARTICLE pubs.acs.org/IECR

Scheduling Gasoline Blending Operations from Recipe Determination to Shipping Using Unit Slots Jie Li and I. A. Karimi* Department of Chemical and Biomolecular Engineering, National University of Singapore, 4 Engineering Drive 4, Singapore 117576 ABSTRACT: Gasoline can yield 60%70% of the total profit of a refinery. Integrated scheduling of gasoline blending and order delivery operations can increase profit by avoiding ship demurrage, improving customer satisfaction, minimizing quality giveaways, avoiding costly transitions and slop generation, and reducing inventory costs. In this paper, we develop a novel continuous-time formulation using unit slots for an integrated scheduling of gasoline blending operations incorporating features such as multipurpose tanks, parallel nonidentical blenders, minimum run lengths, blender setups, tank changeovers, piecewise-constant profiles for blend component qualities and feed rates, simultaneous receipt/delivery by tanks, etc. We employ and revise the schedule adjustment procedure of Li et al. [AIChE J. 2010, 56, 441465] to avoid solving mixed-integer nonlinear programming (MINLP) problems and illustrate the effectiveness of our new formulation on 14 industry-scale examples.

1. INTRODUCTION The petroleum refining industry constantly faces challenges such as fluctuating product demands, ever-changing crude prices, strict environmental regulations, etc. Optimization can play (and has played) a key role in managing oil refineries, as the refiners seek to exploit all potential cost-saving opportunities to survive and thrive in a low-margin business. Oil refineries have used optimization (especially successive linear programming or SLP) for planning refinery operations for a long time. Refinery operations1 involve three main segments, namely, crude oil storage and processing, intermediate processing, and product blending and distribution. While the first involves upstream operations from ship arrivals to crude distillation,2 the last involves downstream operations from the blending of hydrocarbon fractions to the storage and delivery of final products. Gasoline is one of the most profitable products of a refinery and can account for as much as 60%70% of total profit. A refinery typically blends several gasoline cuts or fractions from various processes to meet its customer orders of varying specifications. However, the large numbers of orders, delivery dates, blenders, blend components, tanks, quality specifications, etc. and nonlinear blending make this problem highly complex and nonlinear. If done manually, it can result in suboptimal schedules and costly quality giveaways. Therefore, the scheduling of gasoline blending and delivery is crucial. Scheduling such operations using advanced optimization techniques such as mixed-integer programming (MIP) can help avoid ship demurrage, improve order delivery and customer satisfaction, minimize quality giveaways, reduce transitions and slop generation, exploit low-quality cuts, and reduce inventory costs. The early work related to this problem focused on the optimal blending of various intermediate fractions and some additives to meet product quality specifications.3,4 Glismann and Gruhn5 developed a two-level decomposition approach to integrate short-term scheduling with blend-recipe optimization. They first solve a nonlinear programming (NLP) problem to obtain blend recipes and product quantities, and then they use mixed-integer r 2011 American Chemical Society

linear programming (MILP) to obtain a schedule for the blending operations. Pinto et al.1 and Joly and Pinto6 developed MIP models for scheduling blending operations involving fuel oil and asphalt. They used linear blending correlations for key elements, but they did not address nonlinear product properties such as sulfur. While they considered transitions in blending pipelines, they did not ensure constant blending rates, nor did they enforce minimum run lengths for blend runs. Jia and Ierapetritou7 proposed a continuous-time event-based MILP formulation for scheduling gasoline blending and order delivery operations simultaneously. While their model allows features such as multipurpose product tanks, one product tank delivering multiple orders, and multiple product tanks delivering one order, it lacks some key features, such as multiple parallel nonidentical blenders, variable recipes, and product specifications. However, as shown by Li et al.,8 it gives infeasible solutions and allows a product tank to hold several products at a time. Mendez et al.9 presented both discrete-time and slot-based continuous-time models for the simultaneous optimization of blending and short-term scheduling. Their model allows parallel identical blenders and determines optimal blend recipes. However, (i) it does not ensure constant blend rates or minimum run lengths; (ii) it does not allow multipurpose product tanks or setup times for blender switching; or (iii) it does not integrate order delivery time windows with blending operations. Furthermore, it uses nonlinear blending correlations for some specifications. Recently, Li et al.8 reported a slot-based continuous-time formulation using process slots10 for the simultaneous treatment of recipe, blending, storage, and order scheduling. They incorporated several operational features such as multiple-purpose product tanks, parallel nonidentical blenders, one blender Received: November 16, 2010 Accepted: May 20, 2011 Revised: April 24, 2011 Published: May 20, 2011 9156

dx.doi.org/10.1021/ie102321b | Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

(3) I component tanks, their initial inventories, limits on their holdups, flow profiles of feeds into the tanks, and limits on the flows out of the tanks; (4) P products and specification limits on their property indices; (5) B blenders, the products that each blender can process, lower limits on the blend times of these products, and limits on their blending rates; (6) J Product tanks, the products that each tank can store, limits on their holdups, the products and holdups at time zero, and delivery (lifting) rates for various products; (7) O orders, their constituent products, amounts, and delivery time windows; and (8) Revenues from product sales, component costs, inventory costs (for components and products), and demurrage costs for orders. Determine: (1) The blenders that each component tank should feed over time, and their feed rates; (2) The products that each blender should produce over time, and their production rates; (3) The products that each product tank should receive over time, from which blender, and at what flow rates; (4) The orders that each product tank should deliver over time and their amounts; and (5) The inventory profiles of various tanks (component and product). Assuming: (1) Flow rate profile of each component from the upstream process is piecewise constant; (2) Component quality profile is also piecewise constant; (3) Mixing in each blender is perfect; (4) Changeover times between products are negligible for product tanks; (5) Changeover times between product runs on blenders are product- and sequence-independent; (6) Each order involves only one product; as discussed earlier, each multiproduct order can be decomposed into several single-product orders; (7) Each order is completed during the scheduling horizon; and (8) Each order is single-product, i.e., multiproduct orders have been broken up into individual one-product orders. Subject to the operating rules: (1) A blender can process, at most, one product at any time. Once it begins processing a product, it must operate for some minimum time before it can switch to another product. (2) A blender can feed, at most, one product tank at any time. In addition to being the industry practice, this helps to decrease the number of tanks in use and increase their utilization. Allowing: (1) A component tank may receive and feed components at the same time; (2) A component tank may feed some or all blenders simultaneously; (3) Multiple component tanks may feed a blender at the same time; (4) A product tank may receive and feed products at the same time;

charging, at most, one product tank at a time, minimum run lengths, etc. They ensured the blending rate to be constant during a run, and they developed a schedule adjustment procedure to avoid solving nonconvex MINLPs. However, the use of process slots makes their model slow on larger examples, because of the low RMIPs and slow improvement of RMIPs. In addition, they assumed unlimited component inventories over the scheduling horizon to avoid more slots. This may not always hold in practice. In this paper, we develop a multigrid continuous-time formulation to improve the efficiency of treating recipe, blending, storage, and order scheduling simultaneously. We achieve this by employing unit slots11 in place of process slots. In addition to including all the operational features of the model of Li et al.,8 we allow limited component inventories, simultaneous receipt/ delivery by product tanks, and blender setup times. We also revise the schedule adjustment procedure of Li et al.8 for unit slots. We use 14 examples from Li et al.8 to evaluate our model and algorithm.

2. PROBLEM STATEMENT Consider Figure 1, which shows a gasoline blending unit (GBU) in a typical refinery. It involves I blend components (i = 1, 2, ..., I) in exactly I dedicated component tanks (i = 1, 2, ..., I), B blenders (b = 1, 2, ..., B), J multipurpose product tanks (j = 1, 2, ..., J), and some lifting ports. It makes P grades of gasoline, which we call products (p = 1, 2, ..., P). At time zero, GBDU has O orders (o = 1, 2, ..., O) with time windows [DDLo , DDU o , where o = 1, 2, ..., O] to fulfill during the scheduling horizon [0, H]. An order may involve multiple gasoline products, but we can break each multiproduct order into several single-product orders. Thus, each order involves a single product. A delivery of order o after DDU o incurs a demurrage (DMo, given in dollars per unit time of delay). The blend components are either gasoline fractions from various processes in a refinery (such as atmospheric distillation, fluid catalytic cracking unit (FCCU), and isoforming unit (IFU)) or additives (such as methyl tert-butyl ether (MTBE) and butane), which enhance the octane rating or act as corrosion inhibitors or lubricators. Each component i has a prefixed, distinct, and known quality or specification, The quality of a blend component is specified in terms of various property indices such as octane number (ON), Reid vapor pressure index (RVPI), sulfur index (SI), viscosity index (VI), etc. The flow profiles, over time, of various blend components into respective component tanks are known a priori. The components feed blend components into the blenders in various proportions to make various products over time. The blenders blend products one at a time in a semicontinuous mode. When a blender switches from one product to another, a minimum setup time (τb) is required. A product tank is assigned to receive product from each blend run. A product tank may hold different products over time. Subsequently, the product tanks load the product orders into vehicles or ships at appropriate times for delivery. We refer to this part as order scheduling. The GBDU operation involves decisions such as recipe determination, allocation of component tanks to blenders, assignment of product tanks to products, and scheduling of blending, transfer, and order loading. Thus, the gasoline blending and scheduling problem addressed in this paper can be stated as follows: Given: (1) A scheduling horizon [0, H]; (2) I components and profiles of their property indices over time; 9157

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Figure 1. Schematic of gasoline blending and distribution.

(5) A product tank may deliver multiple orders at the same time; and (6) Multiple product tanks may deliver an order at the same time. Minimizing: The total operating cost includes material (component), demurrage, transition, and back-order costs. The above problem description differs in three aspects from the one addressed by Li et al.8 First, we allow a product tank to receive and deliver products at the same time; Li et al.8 assumed that a tank cannot feed and deliver at the same time. Second, we do not assume unlimited supply of components from the component tanks over the scheduling horizon. In practice, the scheduling of a particular blend may be delayed due to the unavailability of one or more necessary components or an inability to match the quality from the existing supply of components. We now develop a continuous-time, multigrid11 MILP formulation for the above problem. In addition to the distinguishing problem features mentioned above, the main difference between our new formulation and that of Li et al.8 is the use of unit slots instead of process slots. The advantages of using unit slots are well-established11 in the literature. Our contribution is to extend the multigrid methodology of Susarla et al.,11 which was for multipurpose batch plants, to semicontinuous plants. Susarla et al.11 have detailed the differences and similarities between slotbased and event-based formulations. For the ease of exposition, we make two simplifying assumptions, which we relax later. First, we consider the single-period (SP) scenario, in which the flows of all components into the respective component tanks are constant over the entire scheduling horizon. Second, a product tank cannot receive and deliver products simultaneously. Throughout this paper, each variable is defined with specific ranges of its indices, and each constraint, unless otherwise stated, is written for all valid values of the indices of its constituent variables.

3. MILP FORMULATION We define component tanks (i), blenders (b), and product tanks (j) as units (n). For each unit n, we divide the scheduling horizon [0, H] into K (k = 1, 2, ..., K) ordered contiguous unit slots (see Figure 2). We define the period before time zero as slot zero (k = 0). Let Tnk (k = 0, 1, ..., K; Tn0 g 0, TnK e H) denote the

end time of slot k on unit n, where n becomes i for a component tank, b for a blender, and j for a product tank. Slot k on unit n starts at Tn(k1) and ends at Tnk. In contrast to Li et al.,8 the slots are asynchronous, so Tik, Tbk, and Tjk may vary with the unit and are partially independent.11 We define the following sets in this paper: BP = {(b, p) | blender b can process product p} BJ = {(b, j) | blender b can feed product tank j} JP = {(j, p) | product tank j can hold product p} JO = {(j, o) | product tank j may deliver order o} PI = {(p, i) | product p may use component i} OP = {(o, p) | order o is for product p} Since the slots on each unit n are chronologically ordered, we have Tnk g Tnðk1Þ

ðfor 1 e k e KÞ

ð1Þ

We assume that blending begins always from the start of a slot, but may end at any time within the slot. In other words, the idle time, if any, is always toward the end of the slot. Unless otherwise indicated, an index takes all its legitimate values in all the expressions or constraints in our formulation. 3.1. Blending and Storage. At any time, a blender must be either running or idle. However, we assume that it is always connected to a product tank. To facilitate this, we define a dummy product tank (j = 0). Thus, we have J real product tanks (j = 1, 2, ..., J) and one dummy (j = 0) product tank. If a blender is not running, then we connect it to the dummy tank (j = 0), otherwise it must be connected to a real product tank. Furthermore, we assume that a product tank always has a product. If it is empty, then we assume that it is holding the last product that it had. Chart 1 defines the binary variables used to model the blending and storage operations in the GBDU. Li et al.8 also used the aforementioned variables. They treated vbjk (j = 1, 2, ..., J) as binary and vb0k as 01 continuous, and they proved that xbpk would be binary automatically and, hence, can be treated as 01 continuous. Note that we must assign proper values for xbp0, based on the product that blender b was processing before time zero, and ujp0, based on what was inside tank j before time zero. For a formulation based on unit slots, Susarla et al.11 argued that, if a processing unit and a storage unit are exchanging material, then the unit slots in which this exchange occurs must have the same index on both the unit 9158

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

and the storage. On that basis, we use the following constraints from Li et al.8 in our formulation: vb0k þ



∑ vbjk ¼ 1 j:ðb, jÞ ∈ BJ

ujpk ¼ 1

ðfor 1 e k e KÞ

ðfor 1 e j e J, 1 e k e KÞ

xebk þ xbpk þ xbpðkþ1Þ e 2

ðfor ðb, pÞ ∈ BP, 0 e k < KÞ ð7Þ

ð2Þ

uejk g ujpk  ujpðkþ1Þ ðfor ðb, pÞ ∈ JP, 0 e k < K, 1 e j e JÞ

ð3Þ

ð8aÞ

p:ðj, pÞ ∈ JP

vb0k þ

∑ xbpk ¼ 1 p:ðb, pÞ ∈ BP

xbpk g ujpk þ vbjk  1

ðfor 1 e k e KÞ

ðfor ðb, pÞ ∈ JP, 0 e k < K, 1 e j e JÞ

ðfor ðb, pÞ ∈ BP, ðb, jÞ ∈ BJ,

ðj, pÞ ∈ JP, 1 e j e J, 1 e k e KÞ ujpk g xbpk þ vbjk  1

uejk g ujpðkþ1Þ  ujpk

ð4Þ

B K1

ð5aÞ

ðfor ðb, pÞ ∈ BP, ðb, jÞ ∈ BJ,

ðj, pÞ ∈ JP, 1 e j e J, 1 e k e KÞ

∑ ∑ xebk g NP  B b¼1 k¼1

ð5bÞ

xebk g xbpk  xbpðkþ1Þ

ðfor ðb, pÞ ∈ BP, 0 e k < KÞ ð6aÞ

xebk g xbpðkþ1Þ  xbpk

ðfor ðb, pÞ ∈ BP, 0 e k < KÞ ð6bÞ

ð8bÞ ð9Þ

where NP is the number of distinct products that must be processed during the scheduling horizon. Since we will impose a penalty for product transitions in a product tank, we do not force uejk = 0. If a blender b is idle during slot k, then it cannot receive from any component tank during k. Otherwise, it must receive at least one component. Therefore, we write vb0k þ yibk e 1 vb0k þ

I

∑ yibk g 1

i¼1

ðfor 1 e k e KÞ ðfor 1 e k e KÞ

ð10aÞ ð10bÞ

3.2. Blend Scheduling and Product Quality. Let BLbk be the time for which blender b processes a product p (p > 0) during its slot k. Li et al.8 did not need this variable. Clearly, BLbk must equal slot length, unless the current blend run is ending during slot k, or the blender is idle. In other words,

BLbk e Tbk  Tbðk1Þ

ðfor 1 e k e KÞ

BLbk þ Hðvb0k þ uebk Þ g Tbk  Tbðk1Þ

ð11aÞ

ðfor 1 e k e KÞ ð11bÞ

If the blender is idle during slot k, then BLbk must be zero. BLbk e Hð1  vb0k Þ

Figure 2. Schematic of unit-slot design.

ðfor 1 e k e KÞ

ð11cÞ

Chart 1. Binary Variables Used to Model the Blending and Storage Operations in the GBDU

9159

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Next, following Li et al.,8 we define RLbk (0 e k e K) as the length of the current run on blender b at the end of slot k, if the run does not end during slot k; otherwise, RLbk = 0. 8 > > < current run length ðif the current run of blender b does not end during its slot kÞ RLbk ¼ > > :0 ðotherwiseÞ Thus, RLb0 = 0 if a run has ended at time zero; otherwise, it is the current run length at time zero. To compute RLbk, we write RLbk e RLbðk1Þ þ BLbk

ðfor 1 e k e KÞ

RLbk e Hð1  xebk Þ

ðfor 0 e k e KÞ

ð12aÞ

Lastly, Qbk must respect the limits (FLb and FU b ) on the processing rates of blender b. FbL 3 BLbk e Qbk e FbU BLbk



ðfor ðb, pÞ ∈ BP, 1 e k e KÞ I

∑ RLLbpxbpk p:ðb, pÞ ∈ BP

RLbðk1Þ þ BLbk þ RLLb ð1  xebk Þ g

ðfor 1 e k < KÞ RLbðK1Þ þ BLbK g

P



p:ðb, pÞ ∈ BP

where RLLb ¼

max

p:ðb, pÞ ∈ BP

RLLbp xbpK

ð13bÞ

I

J



i¼1

qibk Fi θis g

I



i¼1

I

ðfor 1 e k e KÞ

ð14aÞ

ðfor 1 e k e KÞ

ðfor 1 e k < KÞ

ð14cÞ

Blending requires components and must ensure on-spec products. Let qibk be the volume that blender b receives from component tank i during its slot k; let Gbjk be the volume that it feeds into product tank j; and let Qbk be the volume that it processes. Then, a mass balance around blender b gives us Gbjk ∑ qibk ¼ j:ðb, ∑ i¼1 jÞ ∈ BJ

I

ð19bÞ

!

  qibk Fi θLps  Mb θLps  minðθLps Þ Fmax ð1  xbpk Þ p

ð20aÞ

ðfor ðb, pÞ ∈ BP, 1 e k e KÞ

ð14bÞ

Qbk ¼

p

ðfor ðb, pÞ ∈ BP, 1 e k e KÞ

  RLLbp

STbk e STbðk1Þ þ Tbk  Tbðk1Þ  BLbk

STbk g τb ½xebk  vb0ðkþ1Þ 



i¼1

ð19aÞ

  qibk θis e Qbk θUps þ Mb max ðθUps Þ  θUps ð1  xbpk Þ

ð13aÞ

When a blender switches from one product to another, it requires some minimum setup time (τb). Since we allow a blender to be idle across several consecutive slots after it ends a blend run, we must compute the cumulative idle time between successive blend runs. To this end, we define STbk as the accumulated idle time of a blender b by the end of slot k. We compute this using the following constraints: STbk e Hðxebk þ vb0k Þ

p

i¼1

ð12bÞ

P

ð18a; bÞ

The above allows the blending rate to vary from slot to slot during a blend run. This is not realistic, and we follow the procedure of Li et al.8 to address it later. For gasoline quality, Li et al.8 identified nine most commonly used indices and their additive bases. They proposed the following constraints to ensure on-spec products:   I L L L qibk θis g Qbk θps  Mb θps  minðθps Þ ð1  xbpk Þ

Then, to ensure a minimum length (RLLbp) for each blend run, we demand

ðfor 1 e k e KÞ

ðfor 1 e k e KÞ ð15a; bÞ

If blender b receives nothing from component tank i or is idle, or it does not feed product tank j in slot k, then qibk and Gbjk must be zero. qibk e Mb yibk

ðfor 1 e k e KÞ

ð16Þ

Gbjk e Mb vbjk

ðfor 1 e k e KÞ

ð17Þ

I

∑ qibk Fi θis e i∑¼ 1 qibk Fi i¼1

!

  θUps þ Mb max ðθUps Þ  θUps Fmax ð1  xbpk Þ p

ð20bÞ

ðfor ðb, pÞ ∈ BP, 1 e k e KÞ

  Qbk rpiL  Mb rpiL  minðrpiL Þ ð1  xbpk Þ e qibk e Qbk rpiU p   þ Mb max ðrpiU Þ  rpiU ð1  xbpk Þ p

ðfor ðb, pÞ ∈ BP, 1 e k e KÞ

ð21a,bÞ

where θis is the known blending index for a property s of component i, Fi is the density of component i, [θLps,θU ps] are the desired limits on property s of product p, and Fmax is the maximum possible density among all products. Equations 19 are for volume-based indices, eqs 20 are for weight-based indices, and the volume fraction of component i in product p satisfies limits [rLpi,rU pi]. 3.3. Order Scheduling. Li et al.8 defined one binary variable (zjok) for order scheduling operations: (   1 if product tank j delivers order o during its slot k zjok ¼ 0 ðotherwiseÞ ðfor ðj, oÞ ∈ JO, 1 e j e J, 1 e k e KÞ

They presented the following constraints that are applicable in our new model as well: J

K

∑ ∑ zjok g 1 j:ðj, oÞ ∈ JO k ¼ 1

U U Here, Mb = min[HFU b , max j VPj ] and VPj is the capacity of tank j.

9160

ð22Þ

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Figure 3. An example schedule to illustrate the delivery starts of O1, O2, and O3 by PT-101 at different times in the same slot.

vbjk þ zjok e 1 ðfor ðb, jÞ ∈ BJ, ðj, oÞ ∈ JO, 1 e j e J, 1 e k e KÞ

slot k on tank j. ð23Þ

tsjok g Tjðk1Þ ðfor ðj, oÞ ∈ JO, 1 e j e J, 1 e k e KÞ

zjok e ujpðk1Þ ðfor ðp, jÞ ∈ PJ, ðj, oÞ ∈ JO, ðo, pÞ ∈ OP, 1 e j e J, 1 e k e KÞ

tsjok þ

ð24Þ

ðfor ðp, jÞ ∈ PJ, ðj, oÞ ∈ JO, ðo, pÞ ∈ OP, 1 e j e J, 1 e k e KÞ

ð25Þ zjoðkþ1Þ þ uejk e 1 ðfor ðj, oÞ ∈ JO, 1 e j e J, 0 e k < KÞ

tsjok g DDLo zjok

ð26Þ

ðfor ðj, oÞ ∈ JO, 1 e j e J, 1 e k e KÞ ð27Þ

J



K

∑ DQjok ¼ TQo

do g tsjok þ

ð29Þ

where DQjok denotes the volume of order o delivered by tank j during slot k, TQo is the required amount for order o, and DRU j is the maximum cumulative rate for all orders. While a product tank may deliver multiple orders during a slot, Li et al.8 required them to begin simultaneously. We allow a delivery to start and end at any time during the slot. For instance, consider Figure 3. PT-101 begins delivering O1 at t = 2 and O2 at t = 4 in the same slot that starts at t = 1 and ends at t = 6. To know these times precisely, we define tsjok as the time at which product tank j begins delivering order o in slot k. Since the delivery rate (DRjo) is constant, the end time for this delivery is tsjok þ DQjok/DRjo. We force this delivery to be entirely within

DQjok  DDUo  Hð1  zjok Þ DRjo

ðfor ðj, oÞ ∈ JO, 1 e j e J, 1 e k e KÞ

ð28Þ

j:ðj, oÞ ∈ JO k ¼ 1

ð31Þ

If it ends after DDU o , then we have a delivery delay:

DQjok e DRjU ½Tjk  Tjðk1Þ  ∑ o:ðj, oÞ ∈ JO ðfor 1 e j e J, 1 e k e KÞ

ð30bÞ

Equations 30 allow order delivery to be intermittent from a tank, which is not realistic. We discuss this issue later. In addition to eq 30, we must ensure that a delivery does not begin before its own time window [DDLo , DDU o ].

DQjok e TQo zjok ðfor ðj, oÞ ∈ JO, 1 e j e J, 1 e k e KÞ

DQjok e Tjk DRjo

ðfor ðj, oÞ ∈ JO, 1 e j e J, 1 e k e KÞ

zjok e ujpk

ð30aÞ

ð32Þ

Finally, the delivery delay (do) must be bounded by (H DDU o ). 3.4. Slot Timings on Units. When using unit slots in dealing with shared resources such as inventories, the main challenge is to monitor the resource usage/production by the various sharing units accurately. The time instances corresponding to the flows in/out of a resource must be ordered chronologically, so that we can get a correct resource profile. To this end, we assume that a blending activity is always left justified (in time) within a slot. In other words, a blend activity must begin at the slot start, but may end before the slot end with an idle time. We also recall that each component tank has a constant inflow of material from upstream units in a singleperiod scenario that we are considering for now. Now, as explained in the Appendix, whenever a blender begins a new blend run in a slot, and a component tank is feeding that blender, then we require that the tank have a matching slot with 9161

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Tbðk1Þ þ BLbk e tik þ H½1  yibk þ yibðkþ1Þ 

ðfor 1 e k < KÞ

ð34cÞ Tbðk1Þ þ BLbk g tik  H½1  yibk þ yibðkþ1Þ 

ðfor 1 e k < KÞ

ð34dÞ While eq 33 applies only when a component tank begins to feed a blender during a blend run or a new run starts, eq 34 is effective only when a component tank ends its feed during a blend run or a blend run ends. If a blend run continues across successive slots with the same component tanks, then no new slots begin. Lastly, the synchronizations in eqs 33 and 34 are necessary, because of the continuous nature of blending operations. This is in contrast to the instantaneous flows that could be reasonably assumed in a batch process.11 For product tanks, the situation is simpler, because they cannot receive and deliver at the same time. So flow into and flow out from a product tank can never occur in the same slot. Therefore, if a blender feeds tank j during slot k, then the start (end) of a slot k on a product tank j must precede (succeed) the start (end) of slot k on the blender. Tjk e Tbk þ H½1  vbjðkþ1Þ  ðfor ðb, jÞ ∈ BJ, 1 e j e J, 0 e k < KÞ

ð35aÞ

Tjðkþ1Þ g Tbk þ BLbðkþ1Þ  H½1  vbjðkþ1Þ  ðfor ðb, jÞ ∈ BJ, 1 e j e J, 0 e k < KÞ

3.5. Inventory Balance. The timing constraints enable us to write the following inventory balances for component and product tanks:

Figure 4. Flowchart for the schedule adjustment procedure.

VCik ¼ VCiðk1Þ þ Fi ½Tik  Tiðk1Þ  

the same index and start time. Tbk g Tik  H½2  xebk  yibðkþ1Þ 

ðfor 0 e k < KÞ ð33aÞ

Tbk e Tik þ H½2  xebk  yibðkþ1Þ 

ðfor 0 e k < KÞ ð33bÞ

Similarly, we demand that, if a tank begins feeding a blender during a blend run, then also it should have a matching slot with the same index and start time. In other words, Tbk g Tik  H½1  yibðkþ1Þ þ yibk 

ðfor 0 e k < KÞ ð33cÞ

Tbk e Tik þ H½1  yibðkþ1Þ þ yibk 

ðfor 0 e k < KÞ ð33dÞ

To check the inventory in a tank during or at the end of a blend run, we define an intermediate point tik in each slot, i.e., Ti(k1) e tik e Tik. We then demand that the ends of each blend run and tank feed during a blend run should match the corresponding points on the tank. Tbðk1Þ þ BLbk e tik þ Hð2  xebk  yibk Þ

ðfor 1 e k e KÞ

ð34aÞ Tbðk1Þ þ BLbk g tik  Hð2  xebk  yibk Þ

ð35bÞ

ðfor 1 e k e KÞ

B

∑ qibk b¼1

ðfor 1 e k e KÞ VPjk ¼ VPjðk1Þ þ

B

ð36Þ

O

∑ Gbjk  o∑¼ 1 DQjok b¼1

ðfor ðj, oÞ ∈ JO, 1 e k e KÞ

ð37Þ

where Fi is the constant feed rate to tank i from upstream units, VCik (VCLi e VCik e VCU i ) is the inventory of component i at Tik, ) is the holdup in product tank j at Tjk. VPjk (0 e VPjk e VPU j Equation 37 computes inventory levels only at slot points. Since a blend run may end before the slot, we must also compute the inventory level at tik. Therefore, we use VCLi e VCiðk1Þ þ Fi ½tik  Tiðk1Þ  

B

∑ qibk e VCUi b¼1

ðfor 1 e k e KÞ

ð38a; bÞ

Note that eq 36 neglects the inflows into the component tanks after the end of schedule activities. In other words, eq 36 does not compute the inventory at the end of the horizon if the schedule ends before H. Lastly, a product transition cannot occur on a product tank j, unless the tank holdup at Tjk is zero. Thus, VPjk e VPUj ð1  uejk Þ

ð34bÞ 9162

ðfor 1 e k < KÞ

ð39Þ

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Table 1. Order Data for (a) Examples 19 and (b) Examples 1014 (a) Order Data for Examples 19

o

p

1

2

Amount

Delivery Rate

( 103 bbl)

( 103 bbl/h)

Delivery Window

3 4, 5 6 7, 8 9 1 2 3 4, 5 6 7, 8 9

1

2

3

4, 5

6

7, 8

9

O1 P1 11 10 10 10 10 10 10 3 5 5

5

5

5

5 [0,12]

[0,24]

[0,24]

[0,24]

[0,24]

[0,24]

[0,24]

O2 P2 3

3

3

3

3

3

3 3 3 3

3

3

3

3 [0,12]

[0,24]

[0,24]

[0,24]

[0,24]

[0,24]

[0,24]

O3 P2 60

3

3

3

3

3

3 5 3 3

3

3

3

3 [24,48]

[24,50]

[0,24]

[0,24]

[0,24]

[0,24]

[0,24]

O4 P1 125 10 10 10 10 10 10 5 5 5

5

5

5

5 [24,72]

[24,50]

[0,24]

[0,24]

[0,24]

[0,24]

[0,24]

O5 P2 3 O6 P1 -

3 3 3 3 3 3 3 3 3 10 10 10 10 10 10 - 5 5

3 5

3 5

3 5

3 [24,48] 5 -

[48,72] [48,72]

[24,48] [24,48]

[24,48] [24,48]

[24,48] [24,48]

[24,48] [24,48]

[24,48] [24,48]

O7 P2 -

3

- 3 3

3

3

3

3

-

[48,120]

[24,48]

[24,48]

[24,48]

[24,48]

[24,48]

O8 P1 - 120 100 100 90 100 100 - 5 5

5

5

5

5

-

[118,190]

[118,190]

[118,190]

[118,190]

[118,190]

[118,190]

O9 P2 -

- 3 3

3

3

3

3

-

[144,168]

[144,168]

[144,168]

[144,168]

[144,168]

[144,168]

O10 P4 - 150 150 150 150 150 150 - 5 5

5

5

5

5

-

O11 P3 -

-

20 20 45 60 60 - - 5

5

5

5

5

-

-

[144,168]

[144,168]

[144,168]

[144,168]

[144,168]

O12 P2 -

-

30 30 30 20 20 - - 5

5

5

5

5

-

-

[24,48]

[24,48]

[24,48]

[24,48]

[24,48]

O13 P4 O14 P3 -

-

-

60 45 60 60 - - 10 15 15 20 - - -

5 5

5 5

5 5

5 5

-

-

-

[0,56] [48,72]

[0,56] [48,72]

[0,56] [48,72]

[0,56] [48,72]

O15 P2 -

-

-

20 15 20 20 - - -

4

4

4

4

-

-

-

[0,72]

[0,75]

[0,72]

[0,72]

O16 P2 -

-

-

-

10 20 20 - - -

-

4

5

5

-

-

-

-

[48,72]

[48,72]

[48,72]

O17 P1 -

-

-

-

10 10 10 - - -

-

5

5

5

-

-

-

-

[48,96]

[48,72]

[48,72]

O18 P1 -

-

-

-

10 10 10 - - -

-

5

5

5

-

-

-

-

[48,96]

[48,72]

[48,72]

O19 P2 -

-

-

-

-

60 60 - - -

-

-

5

5

-

-

-

-

-

[0,50]

[0,50]

O20 P2 -

-

-

-

-

40 40 - - -

-

-

5

5

-

-

-

-

-

[144, 168]

[144,168]

O21 P1 O22 P5 -

-

-

-

-

-

30 - - 40 - - -

-

-

-

5 5

-

-

-

-

-

-

[96,120] [144,168]

O23 P3 -

-

-

-

-

-

20 - - -

-

-

-

5

-

-

-

-

-

-

[144,168]

3

3 3

3 3

3 3

3

3

3

3

[150.5,185.5] [150.5,185.5] [150.5,185.5] [150.5,185.5] [150.5,185.5] [150.5,185.5]

(b) Order Data for Examples 1014 Amount ( 103 bbl)

Product o

12 the rest 10

Delivery Rate ( 103 bbl/h)

Delivery Window

11

12

13

14

10

11

12

13

14

10

11

12

13

14

O1 P1

P1

10

10

10

10

10

5

5

5

5

5

[0,24]

[0,24]

[0,24]

[0,24]

[0,24]

O1 P1 O2 P2

P1 P2

10 3

10 3

10 3

10 3

10 3

5 3

5 3

5 3

5 3

5 3

[0,24] [0,24]

[0,24] [0,24]

[0,24] [0,24]

[0,24] [0,24]

[0,24] [0,24]

O3 P2

P2

3

3

3

3

3

3

3

3

3

3

[0,24]

[0,24]

[0,24]

[0,24]

[0,24]

O4 P1

P1

10

10

10

10

10

5

5

5

5

5

[0,24]

[0,24]

[0,24]

[0,24]

[0,24]

O5 P2

P2

3

3

3

3

3

3

3

3

3

3

[24,48]

[24,48]

[24,48]

[24,48]

[24,48]

O6 P1

P1

10

10

10

10

10

5

5

5

5

5

[24,48]

[24,48]

[24,48]

[24,48]

[24,48]

O7 P2

P2

3

3

3

3

3

3

3

3

3

3

[24,48]

[24,48]

[24,48]

[24,48]

[24,48]

O8 P1

P1

100 100 100 100 100

5

5

5

5

5

[118,190]

[118,190]

[118,190]

[118,190]

[118,190]

O9 P2 O10 P4

P2 P4

3 3 3 3 3 150 150 150 150 150

3 5

3 5

3 5

3 5

3 5

O11 P3

P3

60

60

60

60

60

5

5

5

5

5

[144,168]

[144,168]

[144,168]

[144,168]

[144,168]

O12 P2

P2

20

20

20

20

20

5

5

5

5

5

[24,48]

[24,48]

[24,48]

[24,48]

[24,48]

O13 P4

P4

60

60

60

60

60

5

5

5

5

5

[0,56]

[0,56]

[0,56]

[0,56]

[0,56]

O14 P3

P3

20

20

15

20

20

5

5

5

5

5

[48,72]

[48,72]

[48,72]

[48,72]

[48,72]

O15 P2

P2

20

20

20

20

20

4

4

4

4

4

[0,72]

[0,72]

[0,72]

[0,72]

[0,72]

O16 P2

P2

20

20

20

20

20

5

5

5

5

5

[48,72]

[48,72]

[48,72]

[48,72]

[48,72]

O17 P1 O18 P1

P1 P1

10 10

10 10

10 10

10 10

10 10

5 5

5 5

5 5

5 5

5 5

[48,72] [48,72]

[48,72] [48,72]

[48,72] [48,72]

[48,72] [48,72]

[48,72] [48,72]

O19 P2

P2

60

60

60

60

60

5

5

5

5

5

[0,50]

[0,50]

[0,50]

[0,50]

[0,50]

O20 P2

P2

40

40

40

40

40

5

5

5

5

5

[144,168]

[144,168]

[144,168]

[144,168]

[144,168]

[144,168] [144,168] [144,168] [144,168] [144,168] [150.5,185.5] [150.5,185. 5] [150.5,185.5] [150.5,185.5] [150.5,185.5]

9163

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Table 1. Continued (b) Order Data for Examples 1014 3

o

Delivery Rate ( 103 bbl/h)

Amount ( 10 bbl)

Product 12 the rest 10

11

12

13

14

10

11

12

13

14

Delivery Window 10

11

12

13

14

O21 P5

P1

30

30

30

30

30

5

5

5

5

5

[96,120]

[96,120]

[96,120]

[96,120]

[96,120]

O22 P5

P5

40

40

40

40

40

5

5

5

5

5

[144,168]

[144,168]

[144,168]

[144,168]

[144,168]

O23 P3

P3

20

20

20

20

20

5

5

5

5

5

[144,168]

[144,168]

[144,168]

[144,168]

[144,168]

O24 P5 O25 P5

P5 P5

6 20

6 20

6 20

6 20

6 20

3 5

3 5

3 5

3 5

3 5

[96,120] [144,168]

[96,120] [144,168]

[96,120] [144,168]

[96,120] [144,168]

[96,120] [144,168]

O26 P3

P1

-

10

30

10

10

-

4

4

4

4

-

[0,76]

[144,168]

[0,76]

[0,76]

O27 P3

P4

-

20

20

20

20

-

5

4

5

5

-

[120,144]

[72,96]

[120,144]

[120,144]

O28 P4

P1

-

25

3

25

25

-

5

3

5

5

-

[120,144]

[72,96]

[120,144]

[120,144]

O29 P4

P5

-

10

15

10

10

-

5

3

5

5

-

[120,144]

[96,120]

[120,144]

[120,144]

O30 P1

P4

-

15

15

15

15

-

5

3

5

5

-

[120,144]

[96,120]

[120,144]

[120,144]

O31 P2

P1

-

-

15

15

15

-

-

5

5

5

-

-

[96,120]

[120,144]

[120,144]

O32 P5 O33 P1

P1 P4

-

-

20 20

20 20

20 20

-

-

2 5

5 5

5 5

-

-

[96,120] [0,76]

[144,168] [144,168]

[144,168] [144,168]

O34 P3

P4

-

-

20

20

20

-

-

5

5

5

-

-

[120,144]

[168,192]

[168,192]

O35 P3

P5

-

-

30

30

30

-

-

5

5

5

-

-

[120,144]

[168,192]

[168,192]

O36 -

P2

-

-

-

3

3

-

-

-

3

3

-

-

-

[168,192]

[168,192]

O37 -

P1

-

-

-

10

10

-

-

-

5

5

-

-

-

[168,192]

[168,192]

O38 -

P1

-

-

-

40

40

-

-

-

5

5

-

-

-

[168,192]

[168,192]

O39 -

P4

-

-

-

10

10

-

-

-

5

5

-

-

-

[168,192]

[168,192]

O40 O41 -

P5 P1

-

-

-

10 -

10 15

-

-

-

5 -

5 5

-

-

-

[168,192] -

[168,192] [168,192]

O42 -

P2

-

-

-

-

20

-

-

-

-

3

-

-

-

-

[168,192]

O43 -

P3

-

-

-

-

15

-

-

-

-

5

-

-

-

-

[144,168]

O44 -

P5

-

-

-

-

20

-

-

-

-

4

-

-

-

-

[168,192]

O45 -

P4

-

-

-

-

10

-

-

-

-

5

-

-

-

-

[96,120]

3.6. Scheduling Objective. The total operating cost of GBDU from Li et al.8 is given by

TC ¼

I

B

K

B

K1

∑ ∑ ∑ ci qibk þ b∑¼ 1 k∑¼ 1 CBb xebk i¼1 b¼1 k¼1

þ

J

K1

O

∑ ∑ CTj uejk þ o∑¼ 1 DMo do j¼1 k¼1

ð40Þ

where ci is the price of component i ($ per unit volume), CBb is the cost of transition on blender b ($ per occurrence), CTj is the cost of transition in product tank j ($ per occurrence), and DMo is the demurrage cost of order o ($ per unit time). This completes our multigrid formulation for GBDU operations. It comprises eqs 140. It has three main limitations. First, it allows the blending rate to vary from slot to slot and order delivery to be discontinuous; these are undesirable in practice. Second, it assumes constant feed rates for components, or it is a single-period formulation. Third, it does not allow a product tank to receive and deliver product at the same time. For the first limitation, we revise the adjustment procedure of Li et al.8 slightly to accommodate unit slots. The second limitation or the extension to multiple periods is readily addressed as done by Li et al.8 for their model, so we simply state the modified equations without elaboration. However, the third limitation requires a more-detailed consideration. We begin with the first limitation.

4. SCHEDULE ADJUSTMENT A schedule adjustment procedure is necessary to eliminate the slot-to-slot variability of blend rates. The main modification required in the procedure of Li et al.8 concerns the blend length variable (BLbk), which their formulation did not require. Given a solution from our basic model, we compute the correct values for run lengths (CRLbk), volumes (CCQbk), and total volume (TCQbk) processed by a blender in a run as follows: CRLbk ¼ CCQbk ¼ 0

ðif xebk ¼ 1Þ

ð41a,bÞ

CRLbk ¼ CRLbðk1Þ þ ½BLbk 

ðif xebk ¼ 0Þ

ð42aÞ

CCQbk ¼ CCQbðk1Þ þ ½Qbk 

ðif xebk ¼ 0Þ

ð42bÞ

TCQbk ¼ 0

ðif xebk ¼ 0Þ

TCQbk ¼ CCQbðk1Þ þ ½Qbk 

ðif xebk ¼ 1Þ

ð43aÞ ð43bÞ

where [BLbk] and [Q bk] are the values of BL bk and Q bk , repectively. Using the above parameters, the blending rate (R bk) for each blending run at the slot where it ends is given by Rbk ¼ 9164

CCQbðk1Þ þ Qbk CRLbðk1Þ þ BLbk

for k with xebk ¼ 1 and vb0k ¼ 0

ð44Þ

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

1

P3

P3

P2

P2

P1

-

-

-

-

-

C1

C2

C3

C4 C5

C6

C7

C8

C9

tank

PT-101

PT-102

PT-103

PT-104

PT-105

PT-106

PT-107 PT-108

PT-109

PT-110

9165

PT-111

CT-101

CT-102

CT-103

CT-104 CT-105

CT-106

CT-107

CT-108

CT-109

C9

C8

C7

C6

C4 C5

C3

C2

C1

P1

P1

P4

P1 P1

P2

P2

P4

P2

P3

P3

214

Initial Product

13.84

49.34

27.30

19.53

28.29 10.59

59.44

34.58

27.38

-

-

-

-

-

20.20

28.49

14.08

0.00

90.20

1

44.58

49.47

46.91

19.53

44.44 10.59

59.44

67.90

26.46

12.36

60.00

23.96

13.79 12.36

57.59

28.49

25.00

14.08

0.00

30.00

214

Initial Stock ( 103 bbl)

200

100

100

200

250 200

250

200

200

-

-

-

-

-

100

100

150

100

100

1

250

250

250

250

300 200

300

300

250

150

150

200

200 150

150

200

200

150

150

150

214

Capacity ( 103 bbl)

Table 2. Product and Component Tank Data for Examples 114

P2, P3

P2, P3

P2, P3

28

-

-

-

-

-

-

-

-

-

-

-

-

-

P1,P2

-

-

-

-

-

-

-

-

P1,P4

P1,P4

P1,P4

P1,P4 P1,P4

P2,P3

P2, P3

P1,P2 P2P4

P2,P3

P2,P3

P2,P3

1

13, 14

12

-

-

-

-

-

-

-

-

P1,P4

P1,P4

P1,P4

P1,P4 P1,P4

P2,P5

P2,P5

P2P4

-

-

-

-

-

-

-

-

P1,P4

P1,P4

P1,P4

P1,P4 P1,P4

P2,P5

P2,P5

P2P4

-

-

-

-

-

-

-

-

P1,P4

P1,P4

P1,P4

P1,P4 P1,P4

P2,P3,P5

P2,P3,P5

P2P5

P2, P3, P5 P2, P3, P5 P2, P3, P5

P2, P3, P5 P2, P3, P5 P2, P3, P5

P2, P3, P5 P2, P3, P5 P2, P3, P5

911

Storable Products

1.00

0.00

0.00

1.00

1.20 1.00

1.20

1.00

1.00

-

-

-

-

-

15

15

15

15

15

1

1.00

0.00

0.00

0.50

1.00 0.80

1.00

0.50

1.00

20

20

20

20 20

20

20

20

20

20

20

25

1.00

0.00

0.00

0.50

1.00 0.80

1.00

0.50

1.00

25

25

25

25 25

25

25

25

25

25

25

6

1.00

0.00

0.00

0.50

1.00 0.80

1.00

0.50

1.00

20

20

20

20 20

20

20

20

20

20

20

7, 8

1.00

0.00

0.00

0.50

1.00 0.80

1.00

0.50

1.00

30

30

30

30 30

30

30

30

30

30

30

9, 10

1.00

0.00

0.00

0.50

1.00 0.70

1.00

0.50

1.00

30

30

30

30 30

30

30

30

30

30

30

11

1.00

0.00

0.00

0.50

1.00 0.80

1.00

0.50

1.00

30

30

30

30 30

30

30

30

30

30

30

12

Maximum Delivery Rate/Feed Rate ( 103 bbl/h)

1.00

0.50

0.50

0.50

1.00 0.70

1.00

0.50

1.00

30

30

30

30 30

30

30

30

30

30

30

13, 14

Industrial & Engineering Chemistry Research ARTICLE

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Table 3. Component and Product Property Indices for Examples 114 RBN comp./ product

1

28

RVI

SI

BI

AROI

914

28

914

28

914

28

914

28

914

C1

86.50

86.50

86.50

140.47

140.47

80.00

80.00

0.78

0.78

25.00

25.00

C2

103.66

103.66

103.66

68.92

68.92

40.00

40.00

0.98

0.98

31.70

31.70

C3

111.35

111.35

111.35

87.68

87.68

0.00

0.00

1.20

1.20

48.00

48.00

C4

113.93

113.93

113.93

51.47

51.47

5.00

5.00

0.00

0.00

0.00

0.00

C5

94.50

94.50

94.50

175.59

175.59

0.00

0.00

0.10

0.10

0.00

0.00

C6

118.16

118.16

118.16

19.91

19.91

0.08

0.08

0.01

0.01

0.00

0.00

C7

144.68

144.68

144.68

12.55

12.55

7.50

7.50

0.01

0.01

0.05

0.05

C8 C9

150.66 92.50

150.66 92.50

150.66 92.50

110.59 436.34

110.59 436.34

2.00 30.00

2.00 30.00

0.25 0.09

0.25 0.09

19.20 24.00

19.20 24.00

P1

[110.45,þ¥]

[110.45,þ¥]

[110.45,þ¥]

[15, 170]

[15, 170]

[0,45]

[0,45]

[0,0.86]

[0,0.86]

[0, 35.00]

[0, 35.00]

P2

[111.95,þ¥]

[111.95,þ¥]

[111.95,þ¥]

[15, 170]

[15, 170]

[0, 50]

[0, 50]

[0,0.92]

[0,0.92]

[0, 36.00]

[0, 36.00]

P3

[108.97,þ¥]

[108.97,þ¥]

[108.97,þ¥]

[15, 170]

[15, 170]

[0,44]

[0,42]

[0,0.94]

[0,0.94]

[0, 42.00]

[0, 42.00]

P4

-

[103.24,þ¥]

[103.24,þ¥]

[15, 170]

[15, 170]

[0, 50]

[0, 50]

[0,0.90]

[0,0.90]

[0, 40.00]

[0, 40.00]

P5

-

-

[115.01,þ¥]

-

[15, 170]

-

[0, 48]

-

[0,0.93]

-

[0, 40.00]

OI

DNI

FLI

OXI

4

58

914

48

914

48

914

46

7, 8

911

1214

C1 C2

1.00 23.80

1.00 23.80

1.00 23.80

1.49 1.33

1.49 1.33

3.45 6.25

3.45 6.25

0.25 0.75

0.25 0.75

0.25 0.75

0.25 0.75

C3

0.85

0.85

0.85

1.22

1.22

2.36

2.36

2.00

2.00

2.00

2.00

C4

0.00

0.00

0.00

1.58

1.58

3.56

3.56

1.25

1.25

1.25

1.25

C5

0.40

0.40

0.40

1.50

1.50

1.96

1.96

0.08

0.08

0.08

0.08

C6

0.72

0.72

0.72

1.44

1.44

3.65

3.65

0.00

0.00

0.00

0.00

C7

0.00

0.00

0.00

1.15

1.15

2.96

2.96

0.00

0.00

0.00

0.00

C8

0.15

0.15

0.15

1.35

1.35

5.46

5.46

18.20

18.20

18.20

18.20

C9 P1

0.06 [0, 20.00]

0.06 [0, 20.00]

0.06 [0, 20.00]

1.61 [1.19, 1.67]

1.61 [1.19, 1.67]

7.95 [1.4, 7.60]

7.95 [1.4,7.60]

0.85 [0,1.85]

0.85 [0,2.80]

0.85 [0,2.80]

0.85 [0, 2.80]

P2

[0, 18.00]

[0, 18.00]

[0, 18.00]

[1.20, 1.67]

[1.20, 1.67]

[1.4, 7.25]

[1.4, 7.25]

[0,1.90]

[0,2.75]

[4, 7.25]

[0, 2.75]

P3

[0, 20.00]

[0, 20.00]

[0, 20.00]

[1.18, 1.67]

[1.18, 1.67]

[1.4, 7.20]

[1.4, 7.20]

[0,2.10]

[0,2.90]

[0, 2.90]

[0, 2.90]

P4

[0, 18.00]

[0, 18.00]

[0, 18.00]

[1.19, 1.67]

[1.19, 1.67]

[1.4, 7.50]

[1.4, 7.50]

[0,2.00]

[0,2.70]

[0, 2.70]

[0, 2.70]

P5

-

-

[0, 20.00]

-

[1.20, 1.67]

-

[1.4,7.40]

-

-

[0, 3.00]

[0, 3.00]

Then, we set R bk for all slots within each run to be the above value. Thus, if a run spans slots 36 inclusive, then we set R bk = R b6 for k = 35. Now, to obtain a realistic schedule with the constant blend rates computed above, we fix xbpk, xebk, and vb0k and follow the procedure of Li et al.8 This gives us the following reduced formulation, where [xbpk], [xebk], and [vb0k] respectively denote the optimal values of xbpk, xebk, and vb0k obtained from the full model. Equations 11a11c become BLbk ¼ 0 BLbk ¼ Tbk  Tbðk1Þ

for ðb, kÞ with ½vb0  ¼ 1

ð45aÞ

Equations 12 and 13 become RLbk e RLbðk1Þ þ BLbk for ðb, kÞ with ½xebk  ¼ ½vb0k  ¼ 0 ð46Þ RLbk ¼ 0 for ðb, kÞ with ½vb0k  ¼ 1 or ½xebk  ¼ 1 and ½vb0k  ¼ 0

RLbðk1Þ þ BLbk g

P



p¼1

ð47Þ RLLbp ½xbpk  for ðb, kÞ with ½xebk 

¼ 1 and ½vb0k  ¼ 0, ðb, pÞ ∈ BP, 1 e k < KÞ ð48aÞ RLbðk1Þ þ BLbk g

P

RLLbp ½xbpk  ∑ p¼1

for ðb, kÞ with ½vb0k 

¼ 0, ðb, pÞ ∈ BP, k ¼ KÞ

for ðb, kÞ with ½xebk  ¼ ½vb0k  ¼ 0 ð45bÞ

ð48bÞ

Equations 18 becomes BLbk e Tbk  Tbðk1Þ

Qbk ¼ 0

for ðb, kÞ with ½xebk  ¼ 1 and ½vb0k  ¼ 0

Qbk ¼ Rbk 3 BLbk

ð45cÞ 9166

for ðb, kÞ with ½vb0k  ¼ 1 for ðb, kÞ with ½vb0k  ¼ 0

ð49aÞ ð49bÞ

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Table 4. Allowable Composition Ranges for Components in Products of Examples 114 C1

C2

product

1

28

914

1

P1

[0,0.22]

[0,0.22]

[0,0.22]

[0.1,1]

P2

[0,0.24]

[0,0.24]

[0,0.24]

[0.1, 1]

P3

[0,0.25]

[0,0.25]

[0,0.25]

[0.1, 1]

P4

-

[0,0.24]

[0,0.24]

P5

-

-

[0, 0.30]

28

1

28

1

28

914

1

24

57

8

912

13, 14

[0.10, 1]

[0.10, 1]

[0,1]

[0,1]

[0,1]

[0,0.40]

[0,0.40]

[0,0.40]

[0,0.40]

[0,0.40]

[0,0.40]

[0.10, 1]

[0.10, 1]

[0,1]

[0,1]

[0,1]

[0,0.45]

[0,0.45]

[0,0.45]

[0,0.45]

[0,0.45]

[0,0.45]

[0.10, 1]

[0.10, 1]

[0,1]

[0,1]

[0,1]

[0,0.43]

[0,0.43]

[0,0.43]

[0,0.43]

[0,0.43]

[0,0.43]

-

[0.10, 1]

[0.10, 1]

-

[0,1]

[0,1]

[0,0.44]

[0,0.44]

[0,0.44]

[0,0.44]

[0,0.44]

-

-

[0.15, 1]

-

-

[0,1]

-

-

-

[0,0.40]

[0,0.40]

C6

914

C4

914

C5 product

C3

1

28

914

-

C7 1

C8

28

914

1

C9

28

914

1

28

914

P1

[0, 0.25] [0, 0.25] [0, 0.25] [0, 0.20] [0, 0.20] [0, 0.20] [0, 0.25] [0, 0.25] [0, 0.25] [0, 0.30] [0, 0.30] [0, 0.30] [0, 0.15] [0, 0.15] [0, 0.15]

P2 P3

[0, 0.25] [0, 0.25] [0, 0.25] [0, 0.22] [0, 0.22] [0, 0.22] [0, 0.25] [0, 0.25] [0, 0.25] [0, 0.30] [0, 0.30] [0, 0.30] [0, 0.18] [0, 0.18] [0, 0.18] [0, 0.25] [0, 0.25] [0, 0.25] [0, 0.18] [0, 0.18] [0, 0.18] [0, 0.25] [0, 0.25] [0, 0.25] [0, 0.30] [0, 0.30] [0, 0.30] [0, 0.20] [0, 0.20] [0, 0.20]

P4

-

P5

-

[0, 0.25] [0, 0.25] -

-

[0, 0.25]

[0, 0.20] [0, 0.20]

-

-

[0, 0.20]

-

[0, 0.25] [0, 0.25] -

∑ vbjk ¼ 1  ½vb0k 

j¼1

ðfor ðb, jÞ ∈ BP, 0 < j e J, 0 < k e KÞ ð50Þ

I

∑ qibk θis g

ð51aÞ ðfor ðb, pÞ ∈ BP, ðb, jÞ ∈ BJ, ðp, jÞ ∈ PJ, 0 < j e J, 0 < k e KÞ

ð51bÞ



I



i¼1

[0, 0.17]

 Qbk θLps

 Mb

ð52aÞ

ð55bÞ

 θLps



minðθLps Þ p

ð1  ½xbpk Þ ð56aÞ

  qibk θis e Qbk θUps þ Mb max ðθUps Þ  θUps ð1  ½xbpk Þ p

I

ð56bÞ

I

ðθLps ÞgFmax ð1  ½xbpk Þ ∑ qibkFi θis g ði∑¼ 1 qibk Fi ÞθLps  Mb fθLps  min p i¼1

I



i¼1

yibk

ðfor ðb, kÞ with ½vb0k  ¼ 0, 1 e k e KÞ

-

ðfor ðb, kÞ with ½vb0k  ¼ 0, ðb, pÞ ∈ BP, 0 < k e KÞ

I

i¼1

[0, 0.16] [0, 0.16]

-

ðfor ðb, kÞ with ½vb0k  ¼ 0, ðb, pÞ ∈ BP, 0 < k e KÞ

ujpk g ½xbpk  þ vbjk  1

1  ½vb0k  e

-

[0, 0.30]

ðfor ðb, kÞ with ½vb0k  ¼ 0, ðb, pÞ ∈ BP, 0 < k e KÞ

ðfor ðb, pÞ ∈ BP, ðb, jÞ ∈ BJ, ðp, jÞ ∈ JP, 0 < j e J, 0 < k e KÞ

yibk e 1  ½vb0k  ðfor ðb, kÞ with ½vb0k  ¼ 1, 1 e k e KÞ

-

ðfor ðb, kÞ with ½xebk  ¼ 1, 1 e k e KÞ

i¼1

½xbpk  g ujpk þ vbjk  1

[0, 0.30] [0, 0.30]

-

Tbðk1Þ þ BLbk g tik  Hð1  yibk Þ

Fixing the values of xbpk, xebk, and vb0k, eqs 2, 5, 10, 14a, 14c, 1921, 33a, 33b, 34a, 34b, and 40 become J

-

[0, 0.25]

I



qibk Fi θis e

i¼1

!   qibk Fi θUps þ Mb maxðθUps Þ  θUps Fmax ð1  ½xbpk Þ p

ðfor ðb, kÞ with ½vb0k  ¼ 0, ðb, pÞ ∈ BP, 0 < k e KÞ

ð52bÞ

ð57aÞ

ð57bÞ

Qbk rpiL  Mb frpiL  minðrpiL Þgð1  ½xbpk Þ e qibk e Qbk rpiU p

STbk ¼ 0 ðfor ðb, kÞ with ½xebk  ¼ 0 and ½vb0k  ¼ 0, 1 e k e KÞ ð53aÞ

þ Mb fmax ðrpiU Þ  rpiU gð1  ½xbpk Þ p ðfor ðb, kÞ with ½vb0k  ¼ 0, ðb, pÞ ∈ BP, 0 < k e KÞ ð58a; bÞ

STbk g τb ðfor ðb, kÞ with ½xebk  ¼ 1 and ½vb0ðkþ1Þ  ¼ 0, 1 e k < KÞ ð53cÞ

TC ¼

Tbk g Tik  H½1  yibðkþ1Þ  ðfor ðb, kÞ with ½xebk  ¼ 1, 0 e k < KÞ

þ ð54aÞ

Tbk e Tik þ H½1  yibðkþ1Þ  ðfor ðb, kÞ with ½xebk  ¼ 1, 0 e k < KÞ

ð54bÞ

Tbðk1Þ þ BLbk e tik þ Hð1  yibk Þ ðfor ðb, kÞ with ½xebk  ¼ 1, 1 e k e KÞ

I

ð55aÞ

B

K

B

K1

∑ ∑ ∑ ci qibk þ b∑¼ 1 k∑¼ 1 CBb ½xebk  i¼1 b¼1 k¼1 J

K1

O

∑ ∑ CTjuejk þ o∑¼ 1 DMo do j¼1 k¼1

ð59Þ

The revised model (RSPM) comprises eqs 1, 3, 8, 14b, 1517, 2232, 33c, 33d, 34c, 34d, 3539, and 4159, whose solution will guarantee constant blend rates with applicable limits and minimum run lengths. As noted by Li et al.,8 the schedule from RSPM may still show intermittent delivery of orders. If a delivery spans contiguous slots, then we can revise the schedule to simply deliver the 9167

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Table 5. Blender and Economic Data for Examples 114 RLb,0, h and CQb,0 ( 103 bbl) blender

14

5

6, 7

812

Minimum and Maximum Blending Rate ( 103 bbl/h)

13, 14

1

2

36

7

8, 9

1012

Allowable Product

13, 14

1

27

8

912 13, 14

B-101 0 and 0 10 and 150 0 and 0 0 and 0 0 and 0 2.015 1.515 1.520 1.525 1.525 1.530 1.530 P1P3 P1P4 P1P4 P1P5 P1P5 B-102

-

-

-

B-103

-

-

-

0 and 0 0 and 0 -

0 and 0

-

-

-

-

-

-

-

-

1.525 1.530 1.530 -

-

1.525

-

-

-

-

P1P4 P1P5 P1P5 -

-

P1P5

Minimum Blend Run Length (h) P1 17

P2

812

P3

13, 14

17

812

13

14

17

P4

812

13, 14

27

P5

811

12

13, 14

911

12

13, 14

B-101

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

5

5

B-102

-

6

6

-

6

6

6

-

6

6

-

6

6

6

6

5

5

B-103

-

6

-

-

6

6

-

-

6

-

-

-

6

-

-

5

Cost ($/bbl) Example

C1

C2

C3

C4

C5

C6

C7

C8

C9

Example

demurrage cost ( 103 $/h)

scheduling horizon (h)

1 214

20 20

40 24

55 30

45 25

25 22

50 27

58 50

60 50

30 22.5

1 214

2.1 2.5

72 192

Transition Cost ( 103 $/instance) Example PT-101 PT-102 PT-103 PT-104 PT-105 PT-106 PT-107 PT-108 PT-109 PT-110 PT-111 transition cost in blender ( 103 $/instance) 1

9.8

9.8

14.5

9.8

98

-

-

-

-

-

214

14.5

14.5

19

19

19

14.5

19

14.5

19

14.5

complete order at a constant rate.8 Figure 4 shows the complete schedule adjustment procedure.

5. MULTIPERIOD EXTENSION Following the approach of Li et al.,8 we modify eqs 19 and 20 to ensure on-spec products for the piecewise-constant profiles of the multiperiod scenario. I

∑ qibk θist g

i¼1

Qbk θLps

 Mb fθLps



minðθLps Þgð1  xbpk Þ p

ðfor ðb, pÞ ∈ BP, ðt, kÞ ∈ TK, 0 < k e KÞ

ð60aÞ

I

ðθUps Þ  θUps gð1  xbpk Þ ∑ qibk θist e Qbk θUps þ Mb fmax p i¼1 ðfor ðb, pÞ ∈ BP, ðt, kÞ ∈ TK, 0 < k e KÞ I

ð60bÞ

ðfor ðb, pÞ ∈ BP, ðt, kÞ ∈ TK, 0 < k e KÞ



i¼1

20

Furthermore, eqs 36 and 38 change as follows: VCik ¼ VCiðk1Þ þ

B

∑ Fit ½Tik  Tiðk1Þ   b∑¼ 1 qibk t ∈ TK ðfor 1 e k e KÞ

ViL

e VCiðk1Þ þ



t ∈ TK

Fit ½tik  Tiðk1Þ  

B

∑ qibk e

b¼1

ðfor 1 e k e KÞ

qibk Fit θist e ð

I



i¼1

ð61aÞ

qibk Fit ÞθUps þ Mb fmax ðθUps Þ  θUps gFmax t ð1  xbpk Þ p

ðfor ðb, pÞ ∈ BP, ðt, kÞ ∈ TK, 0 < k e KÞ

ð61bÞ

where TK = {(t, k) | slot k is in period t}, θist denotes the known blending index for a property s of component i during period t, Fit is the density of component i during period t, and Fmax is the t maximum possible density among all products during period t.

ð62Þ

ViU ð63a; bÞ

where Fit is the flow rate of component i in period t. With this, our multiperiod model (MPM) comprises eqs 118, 2135, 37, 3940, and 6063. For revised MPM (RMPM), we only need to change eqs 60 and 61, as follows:   I ∑ qibkθist g Qbk θLps  Mb θLps  minðθLps Þ ð1  ½xbpk Þ p

i¼1

ðfor ðb, kÞ with ½vb0k  ¼ 0, ðb, pÞ ∈ BP, ðt, kÞ ∈ TK, 0 < k e KÞ

ð64aÞ

I

ðθLps ÞgFmax ∑ qibk Fit θist g ði∑¼ 1 qibk Fit ÞθLps  Mb fθLps  min t ð1  xbpk Þ p i¼1

I

14.5

I



i¼1

qibk θist e Qbk θUps þ Mb fmax ðθUps Þ  θUps gð1  ½xbpk Þ p

ðfor ðb, kÞ with ½vb0k  ¼ 0, ðb, pÞ ∈ BP, ðt, kÞ ∈ TK, 0 < k e KÞ I

I

∑ qibk Fit θist g i∑¼ 1 qibk Fit

i¼1

!

ð64bÞ

  θLps  Mb θLps  minðθLps Þ Fmax t ð1  ½xbpk Þ p

ðfor ðb, kÞ with ½vb0k  ¼ 0, ðb, pÞ ∈ BP, ðt, kÞ ∈ TK, 0 < k e KÞ ð65aÞ !   I I qibk Fit θist e qibk Fit θUps þ Mb max ðθUps Þ  θUps Fmax t ð1  ½xbpk Þ



i¼1



i¼1

p

ðfor ðb, kÞ with ½vb0k  ¼ 0, ðb, pÞ ∈ BP, ðt, kÞ ∈ TK, 0 < k e KÞ ð65bÞ 9168

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Table 6. Periods, Durations, Slots, and Feed Rates to Component Tanks for Examples 114 Feed Flow Rate to Component Tank ( 103 bbl/h) example

period

duration

slot

CT-101

CT-102

CT-103

CT-104

CT-105

CT-106

CT-107

CT-108

CT-109

1

1

40

1, 2

1.2

0.8

1.2

1.2

0.5

0.8

0.0

0.0

1.0

2

32

3, 4

0.8

0.6

0.6

0.8

0.5

0.6

0.5

0.5

0.0

1

60

1, 2

1.2

0.8

1.2

1.2

0.5

0.8

0.0

0.0

1.0

2

132

3

0.8

0.6

0.6

0.8

0.5

0.6

0.5

0.5

0.0

1

100

1, 2

1.2

0.8

1.2

1.2

0.5

0.8

0.0

0.0

1.0

2

92

3

0.8

0.6

0.6

0.8

0.5

0.6

0.5

0.5

0.0

5

1

100

13

1.2

0.8

1.2

1.2

0.5

0.8

0.0

0.0

1.0

6

2 1

92 120

4, 5 1, 2

0.8 1.2

0.6 0.8

0.6 1.2

0.8 1.2

0.5 0.7

0.6 0.8

0.5 0.0

0.5 0.0

0.0 1.0

2

72

3, 4

0.8

0.6

0.6

0.8

0.5

0.6

0.5

0.5

0.0

1

80

13

1.2

0.8

1.2

1.2

0.7

0.8

0.0

0.0

1.0

2

70

46

0.8

0.6

0.6

0.8

0.5

0.6

0.5

0.5

0.0

3

42

7

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

1

80

1, 2

1.2

0.8

1.2

1.2

0.5

0.8

0.0

0.0

1.0

2

70

3, 4

0.8

0.6

0.6

0.8

0.5

0.6

0.5

0.5

0.0

3 1

42 80

5, 6 1, 2

0.0 1.0

0.0 0.5

0.0 1.0

0.0 1.0

0.0 0.5

0.0 0.5

0.0 0.0

0.0 0.0

0.0 1.0

2

70

3, 4

0.8

0.6

0.6

0.8

0.5

0.6

0.5

0.5

0.0

3

42

5, 6

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

1

80

1, 2

1.0

0.5

1.0

1.0

0.8

0.5

0.0

0.0

1.0

2

60

35

0.8

0.6

0.6

0.8

0.5

0.6

0.5

0.5

0.0

3

52

6, 7

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

1

80

13

1.0

0.5

1.0

1.0

0.7

0.5

0.0

0.0

1.0

2 3

60 52

46 7

0.8 0.5

0.6 0.5

0.6 0.5

0.8 0.5

0.5 0.5

0.6 0.5

0.5 0.0

0.5 0.0

0.0 0.5

1

50

13

1.0

0.5

1.0

1.0

0.8

0.5

0.0

0.0

1.0

2

50

46

0.8

0.6

0.6

0.8

0.5

0.6

0.5

0.5

0.0

3

50

79

0.5

0.5

0.5

0.5

0.5

0.5

0.0

0.0

0.5

4

42

1012

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

1

50

13

1.0

0.5

1.0

1.0

0.7

0.5

0.5

0.5

1.0

2

50

47

0.8

0.6

0.6

0.8

0.5

0.6

0.5

0.5

0.0

3 4

50 42

810 1113

0.5 0.0

0.5 0.0

0.5 0.0

0.5 0.0

0.5 0.0

0.5 0.0

0.0 0.0

0.0 0.0

0.5 0.0

1

50

13

1.0

0.5

1.0

1.0

0.7

0.5

0.5

0.5

1.0

2

50

47

0.8

0.6

0.6

0.8

0.5

0.6

0.5

0.5

0.0

3

50

810

0.5

0.5

0.5

0.5

0.5

0.5

0.0

0.0

0.5

4

42

1113

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

2 3, 4

7

8

9

10

11

12

13

14

Thus, RMPM comprises eqs 1, 3, 8, 14b, 1517, 2232, 33c, 33d, 34c, 34d, 35, 37, 39, 4155, 58, 59, and 6265. So far, we do not allow a product tank to receive and deliver simultaneously. Next, we address the aforementioned third limitation.

As we did for xebk and uejk, we write the following to compute and force zejok to be binary:

6. SIMULTANEOUS RECEIPT/DELIVERY BY PRODUCT TANKS If a product tank can receive and deliver simultaneously, then we need one more 01 continuous variable (zejok) to denote the end of the current delivery run as follows:

zejok g zjoðkþ1Þ  zjok

(

zejok ¼



1 if tank j ends current delivery run of order o in slot k 0 ðotherwiseÞ

ðfor ðj, oÞ ∈ JO, 1 e j e J, 1 e k < KÞ



zejok g zjok  zjoðkþ1Þ

ðfor ðj, oÞ ∈ JO, 1 e j e J, 1 e k < KÞ

ð66aÞ ðfor ðj, oÞ ∈ JO, 1 e j e J, 1 e k < KÞ

ð66bÞ zejok þ zjok þ zjoðkþ1Þ e 2

ðfor ðj, oÞ ∈ JO, 1 e j e J, 1 e k < KÞ

ð67Þ If tank j does not end delivering order o in slot k, then it must continue in (k þ 1). In other words, the delivery in slot k must end exactly at Tjk. Similarly, if the delivery during slot k does not end, then the delivery in slot (kþ1) must begin 9169

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Table 7. Solution Statistics of Various Algorithms/Codes for SPM for Examples 114a order Example 1 Example 2 Example 3

5 10 12

number of

number of

number of

cost

total CPU

number of

number of slots

discrete variables

continuous variables

constraints

( 103 $)

time (s)

switching operations

unit slot

4

136

429

1577

5149.73

0.89

2

process slotc

4

100

269

1098

5149.73

0.72

2

unit slot

3

202

733

2743

3658.11

2.12

1

process slotc

4

252

598

2847

3658.11

6.7

1

model

unit slot

3

246

799

3047

3159.12

process slotc

4

315

633

3199

3159.12

18.8

2.62

1

1

3 4

285 372

925 723

3551 3823

4556.67 4556.67

4.20 39.4

1 1

Example 4

15

unit slot process slotc

Example 5

15

unit slot

3

285

925

3563

4556.67

2.41

1

process slotc

4

372

723

3867

4556.67

6.94

1

Example 6

18

unit slot

4

476

1320

5372

5213.88

3024

process slotc

5

576

950

5438

5213.88

10801

2 2

unit slot

5

672

1745

7323

8100.35

10802

3

Example 7

20

process slotc

6

775

1197

7175

8100.35

10801

3

Example 8

20

unit slot process slotc

4 5

595 682

1548 1164

7327 7664

8080.35 8080.35

10819 10801

2 2

Example 9

23

unit slot

4

635

1651

8083

10573.65

10814

4

process slotc

8

1251

1890

13922

10640.05

14400b

6 4

Example 10 Example 11

25 30

unit slot

7

1265

2913

15001

11286.10

14400b

process slotc

9

1503

2209

16492

11326.10

10817

6

unit slot

8

1644

3729

19679

13248.58

46800b

4

process slotc

12

2328

3216

24755

13437.25

46800b

10 5 13

Example 12

35

unit slot process slotc

9 16

2237 3846

4994 4963

26135 39520

14764.86 15475.34

46800b 46800b

Example 13

40

unit slot

12

3378

7222

43139

15646.15

118800b

6

process slot

17

4414

5966

51775

15879.39

108282

13

unit slot

12

3664

7877

46529

17737.39

118800b

9

17

4830

6439

57657

18249.07

118800b

13

c

Example 14

45

c

process slot a

Note: CPU time limit for SPM is set at 10 800 s for Examples 110, 36000 s for Examples 11 and 12, 108000 s for Examples 13 and 14. CPU time limit for RSPM is set at 3600 s for Examples 110, and 10800 s for Examples 1114. b Reached total CPU time limit. c Process-slot from Li et al.8

tsjok e Tbðk1Þ þ Hð2  vbjk  zjok Þ ðfor ðj, oÞ ∈ JO, 1 e k e KÞ

exactly at Tjk. Equation 31, along with the following expressions, ensures this. tsjok þ

DQjok g Tjk  H½1  zjok þ zejok  DRjo

ðfor ðj, oÞ ∈ JO, 1 e j e J, 1 e k e KÞ

tsjok þ ð68aÞ

tsjoðkþ1Þ e Tjk þ H½1  zjoðkþ1Þ þ zejok  ðfor ðj, oÞ ∈ JO, 1 e j e J, 1 e k e KÞ

ð68bÞ

When a product tank can receive and deliver products at the same time, the flows into/out of the tank can occur in the same slot. Then, eq 36 cannot ensure a correct inventory profile. We need additional constraints. Since we define tsjok as the start of delivery of o by tank j in slot k and the end of delivery is given by tsjok þ DQjok/DRjo. Now, if a blender b is feeding a tank j during its slot k, and tank j is delivering order o at the same time, then the start/end times of delivery must match the start/end times of blender during slot k. This gives us the following additional timing constraints for product tank j: tsjok g Tbðk1Þ  Hð2  vbjk  zjok Þ ðfor ðj, oÞ ∈ JO, 1 e k e KÞ

tsjok þ

ð69aÞ

DQjok g Tbðk1Þ þ BLbk  Hð2  vbjk  zjok Þ DRjo ðfor ðj, oÞ ∈ JO, 1 e k e KÞ

ð69bÞ

ð70aÞ

DQjok e Tbðk1Þ þ BLbk þ Hð2  vbjk  zjok Þ DRjo ðfor j, oÞ ∈ JO, 1 e k e KÞ ð70bÞ

Thus, the model equations for the simultaneous receipt and delivery by product tanks are SPM: eqs 122, 2440, and 6670. RSPM: eqs 1, 3, 8, 14b, 1517, 22, 2432, 33c,33d, 34c34d, 3539, 4159, and 6670. MPM: eqs 118, 21, 22, 2435, 37, 39, 40, 6063, and 6670. RMPM: eqs 1, 3, 8, 14b, 1517, 22, 2432, 33c, 33d, 34c, 34d, 35, 37, 39, 4155, 58, 59, 6265, and 6670. Note that allowing simultaneous receipt and delivery increases the number of slots considerably, and it makes the problem much more difficult. 9170

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Table 8. Solution Statistics of Various Algorithms/Codes for MPM for Examples 114a order Example 1 Example 2 Example 3

5 10 12

Example 4

15

Example 5

15

Example 6 Example 7

18 20

Example 8

20

Example 9

23

Example 10

25

model

30

number of

number of

number of

total CPU

cost

number of

slots

discrete variables

continuous variables

constraints

time (s)

( 103 $)

switching operations

unit slot

4

136

429

1703

1.25

5149.73

2

process slotc

5

130

328

1384

6.89

5149.73

2

unit slot

3

202

733

2683

process slotc

5

329

729

3600

0.69 50.7 1.16

3658.11

1

3658.11

1

unit slot

3

246

799

2987

3159.12

1

process slotc

5

411

769

4048

80.0

3179.12

2

unit slot process slotc

3 5

285 486

925 877

3479 4880

1.97 170

4556.67 4556.67

1 1

unit slot

5

531

1461

6121

46.1

4556.67

1

process slotc

5

486

877

4892

56.9

4556.67

1

unit slot

4

476

1320

5404

1037

5213.88

2

process slotc

6

712

1114

6631

10801

5213.88

2

unit slot

7

986

2381

10337

10814

8100.35

3

process slotc

9

1219

1725

10922

14400b

8682.34

3

unit slot process slotc

6 8

949 1159

2244 1779

11491 12464

14170 10831

8082.85 8123.90

2 4

unit slot

6

1015

2391

12635

10817

10573.65

4

process slotc

9

1423

2107

15708

10833

10636.40

5

unit slot process slotc

Example 11

number of

unit slot

7

1265

2913

15469

11735

11306.10

5

10

1685

2436

18368

10802

11349.54

6

7

1419

3289

17179

36075

13248.58

4

process slotc

14

2742

3720

28957

42342

13348.88

5

Example 12

35

unit slot process slotc

12 17

3050 4099

6572 5258

34151 42026

46800 36467

14809.36 15316.24

7 9

Example 13

40

unit slot

13

3677

7801

46771

118800b

15626.15

5

process slot

18

4686

6302

54857

108005

16187.58

5

unit slot

13

3989

8508

49097

118800b

17737.39

9

process slotc

17

4830

6439

57657

118800b

18678.48

5

c

Example 14

45

a

Note: CPU time limit for SPM is set at 10800 s for Examples 110, 36000 s for Examples 11 and 12, 108000 s for Examples 13 and 14. CPU time limit for RSPM is set at 3600 s for Examples 110, and 10800 s for Examples 1114. b Reached total CPU time limit. c Process-slot-based methodology from Li et al.8

7. NUMERICAL EVALUATION Tables 16 show the data for 14 examples (Examples 114) from Li et al.8 These examples vary in structure, size, scale, and complexity to mimic real-life industry scenarios. Example 1 involves nine components, five product tanks, and a scheduling horizon of 72 h (3 days). Examples 214 have 9 components, 11 product tanks, and a horizon of 192 h (8 days). Examples 17 have one blender, Examples 812 have two, and Examples 13 and 14 have three blenders. Example 1 has three products, Examples 28 have four products, Examples 913 have five products, and Example 14 has six products. The number of orders varies from 10 to 45 in various examples. Li et al.8 used nine product properties: octane number (ON), Reid vapor pressure index (RVPI), sulfur index (SI), benzene index (BI), aromatics index (AI), olefin index (OI), specific gravity index (SGI), flammability index (FI), and oxygenates index (OXI). We considered one property index in Example 1, five in Examples 2 and 3, and nine in Examples 414. The blenders are idle at time zero except in Example 5, where the blender is processing product 1. We used CPLEX 12.1.0 on Dell Precision PWS690 computer (Intel Xeon CPU 3.00 GHz, 16.0 GB memory) running Windows XP to solve the 14 examples.

Tables 7 and 8 present the model statistics for unit-slot and process-slot-based SPM and MPM. Unit-slot-based SPM (USSPM) requires fewer slots in Examples 214, compared to process-slot-based SPM (PSSPM). For instance, Example 12 requires 9 unit slots, compared to 16 process slots. While both USSPM and PSSPM attain the optimal solution in Examples 18, the former solves faster than the latter. For instance, USSPM requires 4.20 CPU s in Example 4 versus 39.4 CPU s for PSSPM. In Example 6, PSSPM fails to prove the optimal solution, even after 3 days of CPU time. While USSPM takes 3024 s to give a solution of $5213.88K, PSSPM requires 10 801 CPU s for the same. Examples 714 are unsolvable to optimality by both models within the allocated CPU times in Table 7. However, USSPM obtains better solutions than PSSPM in all cases. For instance, USSPM gives a cost of $10 573.65K, compared to $10 640.05K for PSSPM in Example 9. In Example 14, USSPM gives $17 737.39K, compared to $18 249.07K for PSSPM. In addition, we also observe that the solutions from USSPM involve, at most, nine product transitions (i.e., on blenders and in product tanks). However, USSPM always gives the same or fewer switches than PSSPM for all examples. For instance, USSPM gives four transitions for Example 10, whereas PSSPM gives six. Similar observations hold true for MPM models 9171

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

Figure 5. (a) Blending schedule for Example 12 (35 orders) from RSPM. (b) Order delivery schedule for Example 12 (35 orders) from RSPM.

as well, except for Example 14, where unit-slot-based (USSPM) gives nine transitions, compared to five for process-slot-based MPM (PSMPM). For illustration, Figure 5 gives the final schedule from SPM for Example 12 (35 orders).

8. CONCLUSIONS We developed a novel multigrid continuous-time formulation for addressing recipe, blending, storage, and order

scheduling simultaneously. Our new formulation improves on our previous work 8 in several respects. By using unit slots instead of process slots, it is able to solve larger problems and is computationally more efficient. In addition, it avoids the assumption of unlimited component inventories, allows blender setup times and simultaneous receipt/delivery by tanks, and revises the schedule adjustment procedure of Li et al.8 for unit slots. 9172

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

to acknowledge Prof. Vladimir Mahalec (McMaster University) for his valuable suggestions on some practical features of real blending systems, which we have tried to include in this paper.

’ NOTATION Sets

BP = set of pairs (blender b, product p) such that b can process p BJ = set of pairs (blender b, product tank j) such that b can feed j PJ = set of pairs (product p, product tank j) such that j can hold p JO = set of pair (product tank j, order o) such that j can deliver o Subscripts

i = component and its tank b = blender p = product o = order j = product tank s = gasoline property specification k = slot t = period Superscripts Figure A1. An example to indicate why, if a blender begins a blend run in a slot, each tank feeding that blender must have a matching slot with the same index and start time. It also illustrates why a check on the tank inventory is needed at the end of the blend run.

’ APPENDIX Consider the example in Figure A1, where CT-101 has a capacity of 1.4 kbbl and feeds B-101 only. The feed rate to CT-101 from upstream units is 1.5 kbbl/h, whereas the consumption rate of B-101 is 3.5 kbbl/h. If BT-101 withdraws from CT-101 during [4, 5] h, then the inventory levels are 0.0 kbbl at 3 h, 1.0 kbbl at 6 h, and 1.5 kbbl at 4 h, and 0.5 kbbl at 5 h. Clearly, we must compute and check the inventory levels at 3, 4, 5, and 6 h. However, this cannot be done by using [3,6] h as one slot on CT-101 and [4,5] h as another on B-101 (Figure A1a). The duration [3, 6] h must be broken up into several slots. [3, 4] h can be merged into the preceding slot. A slot must begin on both B-101 and CT-101 at 4 h. Finally, [4, 6] h can be considered a single slot on both B-101 and CT-101, and we must check the inventory level at 5 h using a continuous variable, as illustrated in Figure A1b. This will guarantee no storage violation for CT-101. From this example, we conclude that whenever a blender begins a blend run in a slot, any component tank feeding that blender must have a matching slot with the same index and the same start time. Moreover, we must check the inventory level at the end of each blend run. ’ AUTHOR INFORMATION Corresponding Author

*Tel.: þ65 6516-6359. Fax: þ65 6779-1936. E-mail: cheiak@ nus.edu.sg.

’ ACKNOWLEDGMENT We would like to acknowledge financial support for this work from The Agency for Science, Technology, and Research (A*Star), under Grant No. 052 116 0074. We would also like

U = upper limit L = lower limit Parameters

Ci = price of component i ($ per unit volume) CBb = cost a transition on blender b ($ per occurrence) CTj = cost of transition in product tank j ($ per occurrence) DDLo = start delivery time of order o DDU o = due date of order o DMo = demurrage cost ($ per unit time) of order o DRjo = delivery rate of product tank j to order o fttc = constant feed rate of component i from upstream units in period t = limits on blending rate of blender b FL/U b Fbk = blending rate of blender b in slot k Fi = constant feed rate of component i from upstream units = limits on the fraction of component i in product p rL/U pi RLLbp = minimum run length of blender b for product p TCQbk = volume of blender b processed in slot k TQo = total amount of order o = capacity limits of component tank i VCL/U i VPU j = maximum capacity of product tank j θL/U ps = limits on property indices of product p θis = blending index for property s of component i Fi = density of component i Fmax = maximum density among all products τb = setup time for blender b Binary Variables

ujpk = 1 if product tank j holds product p during slot k vbjk = 1 if blender b feeds product tank j during its slot k yibk = 1 if blender b receives from component tank i during its slot k zjok = 1 if product tank j is delivering order o during slot k 01 Continuous Variables

uejk = 1 if product tank j switches products at the end of its slot k xbpk = 1 if blender b processes product p during its slot k xebk = 1 if blender b ends its current blend run during its slot k 9173

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174

Industrial & Engineering Chemistry Research

ARTICLE

zejok = 1 if product tank j ends current delivery run of order o during its slot k Continuous Variables

CQbk = volume processed of blender b if it does not end its run during slot k do = demurrage for order o DQjok = volume of order o delivered by product tank j during slot k Gbjk = volume that blender b feeds product tank j during slot k qibk = volume of component i used by blender b during slot k Qbk = volume processed in blender b during slot k RLbk = length of the current run of blender b at the end of slot k STbk = accumulated idle time of a blender b by the end of slot k Tnk = end time of slot k on unit n, where n becomes i for a component tank, b for a blender, and j for a product tank TC = total operating cost ($) VCik = inventory of component tank i at the end of slot k VPjk = inventory of product tank j at the end of slot k

’ REFERENCES (1) Pinto, J. M.; Joly, M.; Moro, L. F. L. Planning and Scheduling Models for Refinery Operations. Comput. Chem. Eng. 2000, 24, 2259–2276. (2) Li, J.; Li, W. K.; Karimi, I. A.; Srinivasan, R. Improving the Robustness and Efficiency of Crude Scheduling Algorithms. AIChE J. 2007, 53, 2659–2680. (3) Dewitt, C. W; Lasdon, L. S.; Waren, A. D.; Brenner, D. A.; Melhem, S. A. OMEGA: An Improved Gasoline Blending System for Texaco. Interfaces 1989, 19, 85–101. (4) Rigby, B.; Lasdon, L. S.; Waren, A. D. The Evolution of Texaco’s Blending Systems: From OMEGA to Starblend. Interfaces 1995, 25, 64–83. (5) Glismann, K.; Gruhn, G. Short-Term Scheduling and Recipe Optimization of Blending Processes. Comput. Chem. Eng. 2001, 25, 627–634. (6) Joly, M; Pinto, J. M. Mixed-Integer Programming Techniques for the Scheduling of Fuel Oil and Asphalt Production. Trans. Inst. Chem. Eng., Part A 2003, 81, 427–447. (7) Jia, Z. Y.; Ierapetritou, M. Mixed-Integer Linear Programming Model for Gasoline Blending and Distribution Scheduling. Ind. Eng. Chem. Res. 2003, 42, 825–835. (8) Li, J.; Karimi, I. A.; Srinivasan, R. Recipe Determination and Scheduling of Gasoline Blending Operations. AIChE J. 2010, 56, 441–465. (9) Mendez, C. A.; Grossmann, I. E.; Harjunkoski, I.; Kabore, P. A. Simultaneous Optimization Approach for Off-Line Blending and Scheduling of Oil-Refinery Operations. Comput. Chem. Eng. 2006, 30, 614–634. (10) Liu, Y.; Karimi, I. A. Scheduling Multistage, Multiproduct Batch Plants with Nonidentical Parallel Units and Unlimited Intermediate Storage. Chem. Eng. Sci. 2007, 62, 1549–1566. (11) Susarla, N.; Li, J.; Karimi, I. A. A New Approach to Scheduling Multipurpose Batch Plants using Unit Slots. AIChE J. 2010, 56, 1859–1879.

9174

dx.doi.org/10.1021/ie102321b |Ind. Eng. Chem. Res. 2011, 50, 9156–9174