Scheduling Crude Oil Unloading, Storage, and Processing - Industrial

Ind. Eng. Chem. Res. 2002, 41, 6723-6734

6723

Scheduling Crude Oil Unloading, Storage, and Processing Li Wenkai and Chi-Wai Hui* Chemical Engineering Department, Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong

Ben Hua Chemical Engineering Research Institute, South China University of Technology, Guangzhou, P.R. China

Zhongxuan Tong Juijiang Petrochemical Complex, SINOPEC, Jiujiang, Jiangxi, P.R. China

This paper presents a solution algorithm and effective mathematical formulations for shortterm scheduling of crude oil unloading, storage, and processing with multiple oil types, multiple berths, and multiple processing units. The problem can be formulated as a mixed-integer nonlinear programming (MINLP) model. However, solving a MINLP model directly results in inconsistency in solution quality and time. This paper proposes a solution algorithm that iteratively solves two mixed-integer linear programming (MIP) models and a nonlinear programming (NLP) model, resulting in better quality, stability, and efficiency than solving the MINLP model directly. The proposed MIP formulations disaggregate the three-dimensional binary decision variables in traditional formulations using two-dimensional binary variables, significantly reduce the total number of binary variables required, and hence improve the overall solution time. For handling large-scale problems, heuristics are included in the formulations to further reduce the solution time. Finally, an industrial case study is used to demonstrate the effectiveness of the proposed formulation and solution algorithm. 1. Introduction Mathematical programming has been extensively studied and implemented for long-term plant-wide refinery planning. Some commercial software applied linear programming (LP) models, such as RPMS (Refinery and Petrochemical Modeling System1) and PIMS (Process Industry Modeling System2), have been developed for refinery production planning. It has been pointed out3 that this planning technology can be considered well developed, and startling further progress should not be expected. Although commercial software tailored for refinery short-term scheduling such as OmniSuite4 and PetroPlan5 are available, compared with long-term planning, short-term operation scheduling technology is much less mature. Furthermore, most of the research work in the short-term scheduling area is about the scheduling of batch processes.6 Far fewer papers have reported research on the scheduling of continuous multiproduct plants. However, intense competition requires that plants promptly respond to market requirements so that scheduling of continuous processes has received more and more attention.7 Shah8 adopted a mathematical programming model for crude oil scheduling. In his paper, the problem was decomposed into two smaller ones: a downstream problem and an upstream problem. The downstream problem was solved first, and the upstream problem was solved subsequently. Lee et al.9 addressed the problem of inventory management for a refinery that included an unloading schedule, a transferring schedule, and a charging schedule for crude oils. * To whom correspondence should be addressed. E-mail: [email protected].

Tu¨rkayetal.10 developedacomputerinterface,CRUDEOIL, to solve the problem of inventory management of a refinery according to the model proposed in Lee et al.9 Pinto et al.7 discussed planning and scheduling refinery operations which included planning models, crude oil scheduling, production and distribution scheduling, fuel oil/asphalt production scheduling, and LPG scheduling. Zhang et al.11 presented a method for overall refinery optimization through integration of the hydrogen network and the utility system with the material processing system. They adopted a decomposition approach to make the problem solvable. The material processing was optimized first using a LP model, and then supporting systems are used to minimize operation costs for a process condition determined from the LP model. In modeling of a refinery process, the calculation of crude oil mixing generates nonconvex bilinear constraints. Some researchers have proposed techniques to linearize these bilinear constraints. Sherali and Alameddine12 proposed a new reformulation-linearization technique (RLT) and imbedded it within a provably convergent branch-and-bound algorithm. This RLT process yields a LP problem whose optimal value provides a tight lower bound on the optimal value of the bilinear programming problem. Quesada and Grossmann13 applied the technique proposed by Sherali and Alameddine12 to model process networks consisting of one or several splitters, mixers, and linear process units that involve multicomponent streams. Lee et al.9 also applied the RLT to reformulate bilinear terms in mass balance equations into linear equations. Although the RLT constraints can sometimes replace bilinear constraints, this often leads to inconsistent solutions, in particular when handling storages of multiple oil types. New

10.1021/ie020130b CCC: $22.00 © 2002 American Chemical Society Published on Web 11/16/2002

6724

Ind. Eng. Chem. Res., Vol. 41, No. 26, 2002

Figure 1. Overview of the crude oil unloading and charging system.

modeling techniques and a solution algorithm are proposed in this paper that overcome the bilinear problem, avoiding solving of the problem as a MINLP but still keeping the consistency of the solution. These will be further discussed in the coming sections. 2. Problem Description The system configuration of the scheduling problem in this paper is illustrated in Figure 1. During a given horizon (7 days, 8 days, etc.), several vessels arrive at the docking station according to the crude supply plan. The number of berths at the docking station is different in different refineries. In this paper, example 1 has one berth, and example 2 has four berths. If all of the berths are busy, the vessel has to wait. After being allocated a berth, a vessel can unload its crude oil to storage tanks. Several vessels can unload their crude oils simultaneously if the number of berths is more than one. The crude oil is then transferred from the storage tanks to charging tanks. The function of the storage tanks and charging tanks is different in different refineries. Some refineries do not allow mixing of crude oils in the storage tanks; the mixing of crude oils can happen only in the charging tanks. This is the case of example 1. Some refineries have no charging tanks, so the mixing of crude oils is allowed in the storage tanks. This is the case of example 2. The crude oil mix, either from the storage tanks (if there are no charging tanks) or from the charging tanks, is then fed into crude distillation units (CDUs). The number of CDUs is different in different refineries. There is one CDU in example 1 and two CDUs in example 2. The data needed to start the scheduling are as follows: (i) Data on the number of vessels to arrive, the arrival date of each vessel, the crude amount and crude type in each vessel, the crude composition data if the number of crude oils is more than one in one vessel, the sea waiting cost, and the unloading cost of each vessel. (ii) Data on the docking station: the number of berths and the maximum and minimum flow rates of crude oil from the docking station to the storage tanks. (iii) Data on the storage tanks and/or charging tanks: the inventory cost, the initial inventory, the minimum and maximum capacities of each tank, the range of concentration of crude oils allowed in the tank, and the maximum and minimum flow rates from the storage tanks to the charging tanks and from the charging tanks or storage tanks to the CDUs. (iv) Data on CDUs: the production targets or product demands of each CDU and the changeover cost of each CDU.

(v) If the model is to maximize the overall profit of a refinery, product prices and crude oil costs are required. The scheduling procedure must then determine the following: (i) The allocation of berths to vessels if the number of berths is more than one. (ii) The sea waiting time and unloading duration for each vessel. (iii) The flow rate of crude oils from the vessels to the storage tanks, from the storage tanks to the charging tanks, and from the storage tanks or charging tanks to the CDUs. (iv) Inventory levels in the storage tanks and charging tanks. (v) The distribution of crude oils from the vessels to the charging tanks or storage tanks to meet the concentration requirements of each tank. The scheduling model minimizes the operating cost or maximizes the profit involved in the operation. The operating cost involves vessel unloading costs, sea waiting costs, tank inventory costs, and CDU transition costs. Some operating rules must be obeyed in the scheduling horizon. These rules will be described together with the formulations in section 4. Other assumptions used in this paper are as follows: (i) Crude oil from each vessel can be pumped to any storage tanks through a main pipeline. (ii) The amount of crude oil remaining in the pipeline is neglected. (iii) The time needed for changeover is neglected. (iv) Perfect mixing in the charging tanks or in storage tanks is assumed. 3. Discussion The difficulties in scheduling of crude oil unloading and charging concern primarily (a) the crude oil composition calculations for the storage and charging tanks which introduce bilinear terms that require the model to be solved as a MINLP, (b) the large number of discrete variables, and (c) the large number of constraints and continuous variables. To overcome these difficulties, some formulation techniques are proposed as follows. 3.1. Calculating the Crude Oil Composition in the Storage/Charge Tanks. The calculation of the concentration of key components in the charge tanks will illustrate the problem. The mass balance equation for key component k in charging tank j at time t can be written as

VBKjkt ) VBKOjk +

∑(∑FSBij(tt) × ESik -

ttet i∈I

FBCKjlk(tt)), ∑ l∈L

∀ j ∈ J, ∀ k ∈ K, ∀ t ∈ Stt (a)

The calculation of the concentration of key component EBjkt will generate nonconvex bilinear constraints:

VBKjkt ) EBjkt × VBjt

(b)

FBCKjlkt ) EBjkt × FBCjlt

(c)

EBMINjk e EBjkt e EBMAXjk

(d)

These constraints cannot be solved using a MIP model. Some researchers (e.g., Lee et al.9) suggest reformulation of the above constraints to avoid the bilinear terms.

Ind. Eng. Chem. Res., Vol. 41, No. 26, 2002 6725

They multiply inequality (d) by FBCjlt to avoid the bilinear term in (c) to obtain (e):

FBCjlt × EBMINjk e FBCKjlkt e FBCjlt × EBMAXjk (e) They then multiply inequality (d) by VBjt to avoid the bilinear term in (b):

VBjt × EBMINjk e VBKjkt e VBjt × EBMAXjk (f) Constraints (a)-(d) are then replaced by constraints (a), (e), and (f) to avoid all bilinear terms so that the problem becomes a MIP problem. Constraints (a), (e), and (f) do not explicitly calculate the key component concentration. These constraints alone are not rigorous enough to ensure that the concentration of the key components inside the charging tank is the same as the streamflow out from the tank. This often causes inconsistencies in the resulting solution. This will be discussed further in example 1. Another approach applies RLT to reformulate the bilinear constraints, as presented by Quesada and Grossmann.13 Constraints (g1)-(g4) are used to replace constraint (e):

FBCKjlkt g FBCMINjl × EBjkt + EBMINjk × FBCjlt - FBCMINjl × EBMINjk (g1) FBCKjlkt g FBCMAXjl × EBjkt + EBMAXjk × FBCjlt - FBCMAXjl × EBMAXjk (g2) FBCKjlkt e FBCMAXjl × EBjkt + EBMINjk × FBCjlt - FBCMAXjl × EBMINjk (g3) FBCKjlkt e FBCMINjl × EBjkt + EBMAXjk × FBCjlt - FBCMINjl × EBMAXjk (g4) Constraints (g1)-(g4) are tighter than constraint (e) and also provide a valid lower bound on the global optimum. However, using constraints (a), (f), and (g1)(g4) to replace (a)-(d) still cannot guarantee the consistency of the key component concentrations. In general, if a unit, such as a mixer, has no mass accumulation inside, then its bilinear terms in mass balance constraints can be replaced by individual component flows and reformulated into exact linear constraints. Nonlinear constraints are thus avoided. If mass accumulates in a unit, such as a charging tank, in some cases using individual component flows and solving the model with MIP will give incorrect results. To avoid possible errors, this paper proposes a new algorithm that properly handles the concentrations of key components or volume fractions of crude oils in charging (or storage) tanks without solving the MINLP. 3.2. Decomposing Tri-index Discrete Variables to Bi-index Discrete Variables. Some researchers14 have studied the issue of reducing the overall number of discrete decision variables by decomposing tri-index variables to bi-index variables. In this paper we give a more general formulation that can be used to perform this decomposition. Many cases exist which call for the use of tri-index binary variables. To denote that the crude oil mix in storage tank i charges crude oil to CDU l at time t, a tri-index variable CDilt is needed. To denote a vessel unloading its crude oil to tank i at time t, calls for the

Table 1. Value of XWvit for Different VIvi and VTvt Values VIvi

VTvt

XWvit

VIvi

VTvt

XWvit

1 1

1 0

1 0

0 0

1 0

0 0

tri-index variable XWvit. If the number of storage tanks and vessels is very large, tri-index variables such as XWvit will introduce a lot of discrete variables into the model. It is proposed that, instead of using tri-index binary variable XWvit, two bi-index binary variables, VIvi and VTvt, can represent vessel v unloading its crude oil to tank i and vessel v unloading at time t. Variable XWvit is kept in the new formulation as a continuous positive variable. The assignment of vessel v to tank i at time t is therefore reformulated as

XWvit g VIvi + VTvt - 1

(h1)

XWvit e VIvi

(h2)

XWvit e VTvt

(h3)

The possible values of VIvi, VTvt, and XWvit are listed in Table 1. It is clear that XWvit can only take the values 0 or 1. Using two bi-index binary variables to replace one triindex binary variable can reduce the number of discrete decision variables significantly. For example, if there are 17 vessels (v ) 17), 7 storage tanks (i ) 7), and a scheduling horizon of 7 (t ) 7), then defining XWvit as a binary variable will introduce 833 discrete variables into the model. However, defining VIvi and VTvt as binary variables to replace XWvit generates only 238 discrete variables. Thus, the number of discrete decision variables has been decreased by 595 or 71%. Such formulations can shorten the solution time significantly. 3.3. Changeover Calculations. In some refineries, feedstocks are premixed in charging tanks before distillation. Changeover (changing charging tanks) is required when a charging tank is running out of stock or when economic or operating conditions change. The number of changeovers should be minimized because each changeover incurs additional operating costs and may cause production losses during the changeover period. Because changeover is a discrete operation, additional binary variables might be required when the costs of changeover are considered in a mathematical model. The technique for handling the changeover in a mathematical model often depends on the definition of the changeover. Definition D0: If a CDU is charged by charging tank j1 at time t and charged by charging tank j2 at time t + 1, this constitutes changeover. This definition was used by Lee et al.9 and was modeled as

Zjj′lt g CDj′lt + CDjl(t-1) - 1

(i)

In constraint (i), changeover happens when CDU l changes the charging tank from j to j′ at time t, forcing the integer variable Zjj′lt to 1. This constraint, by its nature, creates a tremendous number of integer variables and constraints for a problem with multiple CDUs, making the model very time-consuming to solve. Definition D1: Changeover happens as soon as a charging tank starts to charge a CDU.

6726

Ind. Eng. Chem. Res., Vol. 41, No. 26, 2002

Table 2. Comparison of Constraints (i) and (j) single equations j)2 j)4 j)7

nonzero elements

constraint (i)

constraint (j)

difference

constraint (i)

constraint (j)

difference

341 615 1146

339 547 859

2 68 287

1493 2747 5093

1489 2491 3994

4 256 1099

This can be modeled using the following constraint:

Zjlt g CDjlt - CDjl(t-1)

(j)

Constraint (j) is much more general and simpler than (i) but covers the condition that (i) describes. For a problem with multiple CDUs, some CDUs occasionally might be stopped for a short period of time because of changes in product demand, maintenance, etc. Constraint (i) considers this temporary stop not a changeover, even though the charging tanks before and after the stop might have changed, which is indeed untrue. In the case where the CDUs never stop, constraints (i) and (j) should produce the same result, but constraint (j) should require many fewer constraints and binary variables, allowing it to outperform constraint (i). Table 2 lists a comparison of the number of “single equations” and “nonzero elements” in the model MIP I of example 1 in this paper when constraint (i) or (j) is used. Definition D2: Changeover happens whenever a crude oil tank begins or stops charging a CDU. There are two methods to handle this definition: (i) The Absolute Value Method (a General Method). The absolute value of CDjlt - CDjl(t-1) can be used to indicate the changeover as tank j begins or stops charging CDU l at time t. When the value is zero, no changeover occurs. Otherwise, the charging tank either begins or stops charging the CDU. One major drawback of this approach is that calculating absolute values using a direct function is not allowed in most MIP or MINLP solvers. To overcome this, the following formulation is commonly used:

Xjl(t+1) ) CDjl(t+1) - CDjlt

(k1)

Xjl(t+1) ) CD1jl(t+1) - CD2jl(t+1)

(k2)

Zjl(t+1) ) CD1jl(t+1) + CD2jl(t+1)

(k3)

CD1jl(t+1) e Yjl(t+1)

(k4)

CD2jl(t+1) e 1 - Yjl(t+1)

(k5)

Constraints (k1)-(k5) can be used to calculate the absolute value of CDjlt - CDjl(t-1) with the additional binary variables Yjlt, positive variables CD1jlt and CD2jlt ranging between 0 and 1, and free variable Xjlt. Although (k1)-(k5) can be used for calculating the absolute value of CDjlt - CDjl(t-1), a quite significant number of additional binary and continuous variables and constraints are required. (ii) A Simple Approach. Because CDjlt can only be 0 or 1, there is a better method to calculate Zjlt apart from the absolute value method:

Zjlt g CDjlt - CDjl(t-1) Zjlt g CDjl(t-1) - CDjlt

(L1)

If CDjl(t-1) ) 0 and CDjlt ) 1, constraint (L1) becomes Zjlt g 1 and constraint (L2) becomes Zjlt g -1. Because Zjlt is a positive variable and the objective function is to minimize the value of Zjlt, the value of Zjlt can become either 0 or 1. 4. General Mathematical Formulation 4.1. Objective Function. Operating cost ) unloading cost of crude vessels + sea waiting cost of crude vessels + inventory cost of storage tanks and charging tanks + changeover cost.

∑ (TLv - TFv + 1) × C_uldv +

minimize cost )

v∈Svv

∑ (TFv - TARRv) × C_seav + i∈Si ∑ t∈St ∑ 0.5(VSit +

v∈Svv

i

VSi(t-1)) × C_ivsi +

t

∑ ∑ 0.5(VBjt + VBj(t-1)) × j∈J t∈St t

C_inbj +

∑ ∑ ∑ C_setjl × Zjlt j∈J l∈L t∈St

(obj)

t

If there is no charging tank and the crude is charged directly from the storage tanks to the CDUs, then the inventory cost of charging tanks can be neglected and the changeover variable Zjlt should be replaced by Zilt. 4.2. Constraints. (a) Vessel Unloading Constraints. (i) Throughout the scheduling horizon, each vessel can arrive at the docking station only once.

∑

XFvt ) 1, ∀ v ∈ Svv

(1a1)

t∈Varrvt

(ii) Throughout the scheduling horizon, each vessel can leave the docking station only once.

∑

XLvt ) 1, ∀ v ∈ Svv

(1a2)

t∈Varrvt

(iii) Each vessel leaves the dock after its arrival.

∑

XLv(tt) g XFvt, ∀ (v, t) ∈ Varrvt

(1b)

tt∈Vttt(tt)

(iv) The unloading initiation time.

TFv )

∑

ORDt × DT × XFvt, ∀ v ∈ Svv (1c)

t∈Varrvt

In above constraint, ORDt means the relative position of time interval t in set T. (v) The unloading completion time.

TLv )

∑

ORDt × DT × XLvt, ∀ v ∈ Svv

(1d)

t∈Varrvt

(L2)

The constraints (L1) and (L2) are much simpler and more effective in calculating the absolute value of Zjlt.

(vi) If there are several berths at the dock, a vessel can enter any free berth. The unloading of a vessel begins at the same time or after the unloading of the

Ind. Eng. Chem. Res., Vol. 41, No. 26, 2002 6727

preceding vessel.

∑

XWvt e

∑

XFv(tt), ∀ (v, t) ∈ Varrvt

∑

XLv(tt), ∀ (v, t) ∈ Varrvt

(3a1)

ttet∩tt∈Varrv(tt)

XF(v+1)(tt) g XFvt,

tt∈Varr(v+1)(tt)∩Vttt(tt)

∀ (v, t) ∈ Varrvt, v < Nv (1e)

XWvt e

(3a2)

ttgt∩tt∈Varrv(tt)

∑

XLv(tt) -

∑

XFv(tt),

Nv is the total number of vessels that will arrive during the scheduling horizon. If there is only one berth, a vessel can dock only after the preceding vessel leaves. Thus, the above constraint should be replaced by

XWvt g

∑ tt∈Varr

For the case when a vessel is allowed to unload its crude oil to only one storage tank, defining bi-index variable XWvt and using constraints (3a1)-(3a3) are enough for determining if the vessel is being unloaded. XWvt is forced to be either 0 or 1 by knowing the initial and final day of the unloading. If the minimum unloading rate is 0 (FVSMINvi ) 0), (3a3) can be neglected. Without (3a3), XWvt can be any number between 0 and 1 when the vessel is docked. This does not affect the correctness of the result. For cases where vessels are allowed to unload to more than one storage tank, a tri-index variable XWvit is needed. Constraints (3a1)-(3a3) are replaced by (3a1′) and (3a2′), and constraints (3b)-(3e) are added to the model.

(v+1)(tt)∩Vttt,tt

XF(v+1)(tt) g XLvt, ∀ (v, t) ∈ Varrvt, v < Nv (1e′)

(b) Distribution of Vessels to Berths. These include constraints (2a1)-(2b3). If there is only one berth at the dock, these constraints should not be included in the model. (i) In 1 day, one berth can accommodate only one vessel.

∑

XDvdt e 1, ∀ d ∈ D, ∀ t ∈ Stt

(2a1)

v∈Varrvt

(ii) In 1 day, one vessel can moor at only one berth.

∑ XDvdt e 1,

ttgt∩tt∈Varrv(tt)

ttgt+1∩tt∈Varrv(tt)

∀ (v, t) ∈ Varrvt (3a3)

XWvit e

∑

XFv(tt), ∀ (v, i, t) ∈ Vitvit (3a1′)

∑

XLv(tt), ∀ (v, i, t) ∈ Vitvit (3a2′)

ttet∩tt∈Varrv(tt)

∀ (v, t) ∈ Varrvt

(2a2)

d∈D

XWvit e

ttgt∩tt∈Varrv(tt)

(iii) After arrival, one vessel should moor at one berth.

∑ ∑

XDvdt g 1, ∀ v ∈ Svv

(2a3)

d∈Dt∈Varrvt

XWvit e 1 -

(iv) The mooring time of a vessel can vary only between its minimum and maximum unloading time.

MinVtv e

∑ ∑

XDvdt e MaxVtv, ∀ v ∈ Svv

d∈Dt∈Varrvt

(2a4,2a5)

(v) A vessel can only moor at the berth between the unloading initiation time and completion time.

∑ XDvdt ) ttet ∑XFv(tt) + ttgt∩tt∈Varr ∑ XLv(tt) - 1,

d∈D

vt

(vi) To ensure that a vessel should stay at the same berth throughout its unloading period, another bi-index binary variable XVDvd is introduced. Constraint (2b3) ensures that the value of XVDvd equals unity at one berth for one vessel, and then constraint (2b2) ensures that XDvdt has the value of unity at the same berth during the unloading period.

∑

XDvdt, ∀ d ∈ D, ∀ v ∈ Svv (2b2) ∀ v ∈ Svv

∀ (v, i, t) ∈ Vitvit

(3b)

When crude oils are charged from charging tanks to CDUs, this operating rule is applied to the flow rate from storage tanks to charging tanks (constraint (5b1,5b2)), thus, (3b) should not be included in the model. (iii) A vessel can feed at most one storage tank.

∑

XWvit e 1, ∀ (v, t) ∈ Varrvt

(3c)

(iv) Constraints (3d1)-(3d3) are used to ensure that positive variable XWvit can take only the values 0 or 1.

XWvit g VIvi + VTvt - 1, ∀ (v, i, t) ∈ Vitvit

(3d1)

XWvit e VIvi, ∀ (v, i, t) ∈ Vitvit

(3d2)

XWvit e VTvt, ∀ (v, i, t) ∈ Vitvit

(3d3)

(v) At time t, a vessel can unload its crude oil only when the vessel has already arrived and moored at a berth.

t∈Varrvt

∑ XVDvd ) 1,

CDilt, ∑ l∈L

i∈Vitvit

∀ (v, t) ∈ Varrvt (2b1)

XVDvd × MaxVtv g

(ii) When a storage tank is charging a CDU, it cannot be fed.

XWvit e

∑ XDvdt,

∀ (v, i, t) ∈ Vitvit

(3e)

d∈D

(2b3)

d∈D

(c) Denoting When a Vessel Can Unload Its Crude Oil. These include constraints (3a1)-(3a8). (i) Unloading is possible between time TFv and TLv.

The upper bound of XWvit is 1. (d) Material Balance for Vessels. (i) Total crude oil volume in vessel v at time t ) initial crude oil volume in vessel v - crude oil transferred from vessel v to storage tanks up to time t.

6728

Ind. Eng. Chem. Res., Vol. 41, No. 26, 2002

∑ ∑FVSvi(tt),

(4a)

FBCMINjl × CDjlt e FBCjlt e FBCMAXjl × CDjlt, ∀ j ∈ J, ∀ l ∈ L, ∀ t ∈ Stt, FBCMINjl > 0 (6b1,6b2)

(ii) Crude oil transfer rate from vessel v to storage tank i at time t.

(iii) Volume capacity limitations for charging tank j at time t:

FVSMINvi × XWvit e FVSvit e FVSMAXvi × XWvit, ∀ (v, i, t) ∈ Vitvit ∩ FVSMINvi > 0 (4b1,4b2)

VBMINj e VBjt e VBMAXj

VVvt ) VV0v -

∀ (v, t) ∈ Varrvt

i∈Vitvi(tt)ttet

(iii) The initial crude oil volume of vessel v equals the volume of crude oil transferred from vessel v to storage tanks i during the scheduling horizon.

∑ ∑

FVSvit ) VV0v, ∀ v ∈ Svv

(4c)

i∈Vitvitt∈Vitvit

(g) Total Amount of Crude Oil Charged Should Meet the Demand of the CDUs. (i) The total amount of crude oil transferred from storage tanks to each CDU is restricted by the maximum and minimum demands of the CDU.

DMMINl e

∑ FSCilt e DMMAXl,

∀ t ∈ Stt, ∀ l ∈ L

i∈Stt

Note that, for cases with storage tanks and charging tanks and where no mixing is allowed in storage tanks, such as example 1 in this paper, XWvit should be replaced by XWvt in constraint (4b1,4b2). (e) Material Balance for Storage Tanks. (i) Total crude oil volume in storage tank i at time t ) initial crude oil volume in storage tank i + crude oil transferred from vessels to storage tank i up to time t - crude oil transferred from storage tank i to CDUs up to time t:

VSit ) VS0i +

(6b3)

∑ ∑FVSvi(tt) - ∑ ∑FSCil(tt), l∈L ttet

v∈Vitvi(tt)ttet

∀ i ∈ Sii, ∀ t ∈ Stt (5a) (ii) Crude oil transfer rate from storage tank i to CDU l at time t constraints.

(7a,7b)

(ii) If the total amount of crude oil charged has to meet an exact demand of a CDU, (7a) and (7b) are replaced by (7a′).

∑ ∑ FBCjlt ) DMj, l∈L t∈St

∀j∈J

(7a′)

t

(h) Material Balance of Crude Oil c in a Storage Tank. (i) Volume of crude oil c in storage tank i at time t ) initial volume of crude oil c in storage tank i + crude oil c transferred from vessels to storage tank i up to time t - crude oil c charged into CDUs up to time t:

VSCict ) VSC0ic +

∑(

∑

FVSvitEVvc -

ttet v∈V_Cvc∩Vitvi(tt)

FSCMINil × CDilt e FSCilt e FSCMAXil × CDilt, ∀ i ∈ Sii, ∀ t ∈ Stt, ∀ l ∈ L, FSCMINil > 0 (5b1,5b2)

∑

FSCCilc(tt)),

l∈L,(i,c)∈I_Cic

∀ (i, c, t) ∈ Ι_Cic ∩ Stt ∩ Si (8a) (ii) Flow rate of crude oil c from storage tank i to CDU l:

Note that when there are storage tanks and charging tanks in the problem, crude oil is transferred from the storage tanks to the charging tanks, and variables FSCilt, FSCMINil, and FSCMAXil in constraints (5a) and (5b1,5b2) should be replaced by FSBijt, FSBminij, and FSBmaxij; CDilt in (5b1,5b2) should be replaced by 1 ∑l∈LCDjlt, and index l replaced by j. (iii) Volume capacity limitations for storage tank i at time t:

VSMINi e VSit e VSMAXi (f) Material Balance for Charging Tanks. These include constraints (6a) and (6b1,6b2). If there is no charging tank, these constraints should not be included in the model. (i) Crude oil in charging tank j at time t ) initial crude oil in charging tank j + crude oil transferred from storage tanks to charging tank j up to time t - crude oil charged into CDUs up to time t:

VBjt ) VB0j +

∑ ∑ i∈I ttet

FSBij(tt) -

∑ ∑ l∈L ttet

FBCjl(tt),

∀ j ∈ J, ∀ t ∈ Stt (6a)

(ii) Crude oil transfer rate from charging tank j to CDU l at time t.

FSCilt × ESMINic e FSCCilct e FSCilt × ESMAXic, ∀ l ∈ L, ∀ (i, c, t) ∈ Sii ∩ Stt ∩ I_Cic, ESMINic > 0 (8b1,8b2) (iii) Constraints on the volume of crude oil c in storage tank i at time t:

VSit × ESMINic e VSCict e VSit × ESMAXic, ∀ (i, c, t) ∈ Sii ∩ Stt ∩ I_Cic, ESMINic > 0 (8c1,8c2) (iv) Constraints on the concentration of crude oil c in storage tank i at time t. The concentration limits in storage tank i at time t:

ESMINic e EBCict e ESMAXic, ∀ (i, c, t) ∈ Sii ∩ Stt ∩ I_Cic The flow rate of crude oil c from storage tank i should equal the total flow rate of crude oil from storage tank i multiplied by the concentration of crude oil c in storage tank i:

FSCCilct ) EBCict × FSCilt, ∀ l ∈ L, ∀ (i, c, t) ∈ Sii ∩ Stt ∩ I_Cic (8d1)

Ind. Eng. Chem. Res., Vol. 41, No. 26, 2002 6729

The volume of crude oil c in storage tank i should equal the total volume of crude oil in storage tank i multiplied by the concentration of crude oil c in storage tank i:

VSCict ) EBCict × VSit, ∀ (i, c, t) ∈ Sii ∩ Stt ∩ I_Cic (8d2) The sum of the concentrations of all crude oil types should equal unity:

∑

EBCict ) 1, ∀ (i, t) ∈ Sii ∩ Stt

(8d3)

c∈I_Cic

Constraints (8a)-(8d3) account for problems with crude oils mixing in storage tanks when the CDUs are directly charged from the storage tanks. For problems which concern key component balances, such as the sulfur content, or with charging tanks, the formulation in this section will require some modifications. (i) Crude Oil Charging Operating Rules. (i) A storage tank i can charge at most one CDU at any time t:

CDilt e 1, ∑ l∈L

∀ i ∈ Sii, ∀ t ∈ Stt

(9a1)

(ii) A CDU l can be charged by at most two storage tanks simultaneously at any time t:

∑ CDilt e 2,

∀ l ∈ L, ∀ t ∈ Stt

(9a2)

Figure 2. Procedure of new algorithm to solve the model.

i∈Sii

If a CDU can and must be charged by one tank at any time, then constraint (9a2) becomes (9a2′):

∑ CDilt ) 1,

∀ l ∈ L, ∀ t ∈ Stt

(9a2′)

i∈Sii

(iii) Constraints (9b1)-(9b3) ensure that continuous positive variable CDilt can only take the values 0 or 1.

CDilt g ILil + ITit - 1, ∀ l ∈ L, ∀ i ∈ Sii, ∀ t ∈ Stt (9b1) CDilt e ILil, ∀ l ∈ L, ∀ i ∈ Sii, ∀ t ∈ Stt (9b2) CDilt e ITit, ∀ l ∈ L, ∀ i ∈ Sii, ∀ t ∈ Stt (9b3) Note that if crude oils are charged from charging tanks to CDUs, index i in constraints (9a1)-(9b3) should be replaced by j; CDilt, ILil, and ITit should be replaced by CDjlt, ILjl, and ITjt, respectively. (iv) Changeover calculation formulation for definition D1:

Zjlt ) CDjlt - CDjl(t-1), ∀ l ∈ L, ∀ j ∈ J, ∀ t ∈ (Stt, t > 1) (9c) If changeover definition D2 is applied, then (9c) will be replaced by (9c1) and (9c2):

Zilt g CDilt - CDil(t-1), ∀ l ∈ L, ∀ i ∈ Sii, ∀ t ∈ (Stt, t > 1) (9c1) Zilt g CDil(t-1) - CDilt, ∀ l ∈ L, ∀ i ∈ Sii, ∀ t ∈ (Stt, t > 1) (9c2) 5. A New Algorithm To Solve the Mixed-Integer Nonlinear Problem Model The problem presented in this paper is indeed a mixed-integer nonlinear problem (MINLP). However,

the only nonlinearity is caused by calculating the concentration inside the storage and feed tanks, constraints (8d1) and (8d2). Instead of solving the problem as a MINLP, a new solution algorithm is proposed that solves the problem iteratively as a NLP and a MIP. Figure 2 shows the flow diagram of the solution procedure. Before discussion of the procedure, some explanations of the terms that are shown in Figure 2 are given: MIP I. The first MIP model for the problem. In this model, the bilinear term is linearized using the method in work by Lee et al.9 It includes all constraints listed in section 4 except (8d1)-(8d3). Although the model does not guarantee a consistent concentration inside and in the discharge from the tanks, the model generates a good starting point for the coming calculations. Iteration Condition I. The procedure at this point checks the consistency of the concentration values from MIP I. The consistency criteria are (i) whether the concentration of crude oils in the mixing tank is the same as that of the stream that flows out from the mixing tank and (ii) whether the sum of the concentrations of all crude oil types equals unity. If both (i) and (ii) are positive, an optimum solution has been found and the procedure is terminated. Otherwise, the procedure goes on to the next step. Integer Map Transfer I. The values of discrete variables and 0-1 continuous variables which were used in material balance constraints and CDU charging constraints determined by MIP I are sent to the corresponding parameters in the following NLP model. These variables are XFvt, XLvt, CDilt (or CDjlt), and XWvit (or XWvt). The name “transfer variables” will be used to denote these four variables. For example, the value of discrete variable XFvt determined by MIP I is transferred to parameter PXFvt in the NLP model. NLP. The constraints in the NLP are the same as those in MIP I except for the following:

6730

Ind. Eng. Chem. Res., Vol. 41, No. 26, 2002

(i) Constraints in MIP I that were used to determine the values of “transfer variables” are not included in the NLP model. (ii) The “transfer variables” in MIP I are fixed as parameters at the values that were calculated by MIP I. For instance, constraints (3a1)-(3e) in MIP I are used to determine the value of XWvit, so they are not included in the NLP model. XWvit used in constraints (4b1) and (4b2) is fixed at the present value. (iii) Nonrigorous constraints (8b1), (8b2), (8c1), and (8c2) used to determine the concentration of crude oils are replaced by rigorous constraints (8d1) and (8d2). The NLP model determines accurate concentration values using rigorous bilinear constraints. These concentration data are further used in the next step. In the NLP model, all of the discrete variables are fixed as parameters. Iteration Condition II. The iteration is terminated if the difference between the objective values of the NLP and the previous MIP run is smaller than an appropriate tolerance (the absolute relative difference). The tolerance used in this paper is 1.0 × 10-5. Concentration Transfer. If the previous check is not satisfied, indicating that the concentrations in the tanks strongly affect the objective values, then they may even affect the integer variables that were previously determined in MIP I. To verify this, the concentrations in the tanks are fixed at the values generated from the NLP, but the integer valuables that were previously fixed will be freed in the coming step. MIP II. A new MIP II model is generated from the previous step by fixing the tank concentrations but allowing all of the integer variables to vary. The new model uses PEBCict (or PEBjkt) to replace EBCict (or EBjkt) in (8d1) and (8d2) to construct two new pseudoNLP constraints. These two new constraints replace nonrigorous constraints (8b1), (8b2), (8c1), and (8c2) in MIP I to form MIP II. MIP II uses the concentration data determined by the NLP to determine a new integer map. Iteration Condition III. If the objective value of MIP II is better than that of the NLP model and any value of a “transfer variable” calculated by MIP II is different from the corresponding parameter, then the procedure continues to iterate. Otherwise, the procedure is terminated. Integer Map Transfer II. The values of the “transfer variables” determined by MIP II are fixed at the present values in the next NLP model. The proposed algorithm shown in Figure 2 first solves the problem using MIP I. The concentration values of the tanks are determined; however, they may not be consistent. If the concentration results are correct (by iteration condition I), then the final solution has already been obtained (such as in the case reported by Quesada and Grossmann13). Otherwise, the procedure goes to the next step, using integer map transfer I to transfer values of the “transfer variables” to the NLP as parameters. The NLP model is solved to obtain correct concentrations corresponding to the integer map generated from MIP I. The procedure then uses iteration condition II to decide whether the objective value of MIP I in the previous iteration and the objective value of the NLP in the current iteration are close enough. If yes, no better solution to the problem can be expected, and this stops the whole algorithm. If no, then these concentration data are transferred to MIP II as param-

Table 3. Solution Results of Example 1 in Three Cases case A (MIP) case A (MINLP) case A (new algorithm) case B (MIP) case B (MINLP) case B (new algorithm) case C (MIP) case C (MINLP) case C (new algorithm)

CPU time (s)

objective value

1.6 24.7 1.6 0.3 1.3 0.7 1.32 3.5 1.9

211.2 211.2 211.2 207.9 208.0 208.0 211.0 211.0 211.0

eters. In MIP II, the concentration values remain unchanged, but all of the other variables, including the “transfer variables”, are calculated again. After MIP II is solved, the procedure uses iteration condition III to decide whether a better solution to the problem is possible. If yes, integer map transfer II transfers the newest “transfer variables” to the next NLP model and continues the second iteration until the terminating condition of iteration condition II becomes true. In most cases tested in this study, one to three iterations were enough to get the final solution. 6. Examples Two examples with different problem sizes and complexities demonstrate the capabilities and effectiveness of the new algorithm and formulations proposed in this paper. The models were formulated in GAMS15 on a 933 MHz Pentium III PC. The code DICOPT in GAMS 2.25 was used for MINLP, OSL for MIP, and CONOPT for NLP. 6.1. Example 1. Example 1 is taken from the motivating example of Lee et al.9 In this example, two crude vessels are to arrive on days 1 and 5, and unloading for both vessels should be completed by day 8. There is one CDU, which has to process two kinds of mixed crude. Crude oils are initially unloaded to two storage tanks and then mixed into two charging tanks to make two crude oil mixes X and Y. Only one key component, the sulfur concentration, is used to characterize each crude oil. The sulfur concentration of these two crude oil mixes should be in the ranges of 0.0150.025 for X and 0.045-0.055 for Y. The demand for crude mix X and crude mix Y by the CDU is 1 million bbl each during the scheduling horizon. Example 1 was solved using MIP, MINLP, and the new algorithm separately. The results of the following three cases are shown in Table 3. In case A, all data are the same as those of Lee et al.9 All three algorithms (MIP, MINLP, and the new algorithm in this paper) resulted in the same solution no matter whether the RLT technique was used or not in the MIP model. From case A of Table 3, the solution time with the new algorithm is the same as that with MIP because an optimum solution was found in the first step of the new algorithm. The solution time with MINLP is much longer than the other two. Note that in case A the solution is slightly different from that of Lee et al.9 because some new formulations, such as the new definition of a changeover, have been applied in this paper. In case B, assume that the CDU is under maintenance and the demand of the CDU for crude mixes X and Y during the scheduling horizon is reduced from 1 to 0.4 million bbl. More crude oil is stored as the inventory in the storage tanks and the charging tanks to be processed in the next scheduling horizon. All other data of case B remain the same as those of case A. In

Ind. Eng. Chem. Res., Vol. 41, No. 26, 2002 6731 Table 4. Concentration Data Calculated by MIP for Case B time 1

2

3

4

5

6

7

8

VX 0.02 0.02 0.02 0.018 0.017 0.015 0.015 0.015 VF UNDF UNDF UNDF UNDF 0.025 0.025 0.015 0.015 Table 5. Concentration Data Calculated by MIP for Case C time 1

2

3

4

5

6

7

8

VX UNDF UNDF 0.055 0.065 0.047 0.063 0.015 UNDF VF 0.02 UNDF UNDF 0.015 0.12 0.015 0.11 0.015 Table 6. Concentration Data Calculated by the New Algorithm for Case B time 1

2

3

4

5

6

7

8

VX 0.02 0.02 0.02 0.018 0.018 0.018 0.018 0.018 VF UNDF UNDF UNDF UNDF 0.018 0.018 0.018 0.018 Table 7. Concentration Data Calculated by the New Algorithm for Case C time 1

2

3

4

5

6

7

8

VX UNDF UNDF 0.055 0.055 0.055 0.055 0.055 UNDF VF 0.02 UNDF UNDF 0.055 0.055 0.055 0.055 0.055

case C, the range of sulfur concentration to be processed in the CDU is widened. Thus, the sulfur concentration of the crude oil mix X changes from the range of 0.0150.025 to 0.015-0.12. All other data of case C also remain the same as those of case A. In cases B and C, the MIP approach failed to generate consistent sulfur concentrations in the tanks. The concentrations in charging tank j1 during the scheduling horizon calculated by MIP are listed in Table 4 (case B) and Table 5 (case C). The “VX” row of Tables 4 and 5 is the concentration of component k in the charging tank (VXjt ) VBKjkt/VBjt), and the “VF” row is the concentration of component k in the stream that flows out from the charging tank (VFjt ) FBCKjlkt/FBCjlt). “UNDF” in the tables means that at that time either the volume in the charging tank is 0 (VBjt ) 0) or the flow rate from the charging tank to the CDU is 0 (FBCjlt ) 0). Because the crude oil flow out from charging tank at time t has the same composition as the crude oil in charging tank at the end of time t - 1, VFjt must equal VXj(t-1). However, from Tables 4 and 5, it is clear that this requirement is not satisfied. The results calculated using MIP in both cases B and C are not correct. In case B, the objective value calculated by the MIP model is slightly lower because its concentration constraints are looser than the others. In case C, the objective values from the three algorithms are the same, indicating that the objective value is not sensitive to the concentration in the charging tanks. Both MINLP and the new algorithm get the correct results, but the solution time with the new algorithm is faster than that with MINLP. It takes one iteration for the new algorithm to find the optimal solution in both cases B and C. The concentration data calculated by the new algorithm are listed in Table 6 (case B) and Table 7 (case C). It is found that the concentration requirement has been satisfied. 6.2. Example 2. Example 2 is a real industrial example provided by Juijiang refinery in China with

Figure 3. Flow diagram for example 2.

some modifications in data to maintain the secrecy of the company. The flow diagram for this example is shown in Figure 3. In the example, 17 vessels are planned to moor at the docking station during the scheduling horizon of 7 days. The amount and type of crude oil in each vessel is different. The vessel data are given in Table 8. In Table 8, vessel 1 has 3000 tons of crude oil; it arrives at the docking station on day 1 and should stay at the docking station at least 1 day but at most 3 days. Vessel 1 only has one type of crude oil: C1. Vessel 17 arrives at day 7. It carries equal quantities of two types of crude oil: C2 and C6. The refinery docking station has four berths. Every day, at most four vessels can arrive at the docking station. One of the four berths needs to be allocated to each vessel. A vessel can feed at most one storage tank. A storage tank can be charged simultaneously by a maximum of four vessels. The company has seven storage tanks. Unlike some other refineries, the refinery has no charging tank. The number of crude oil types allowed in the vessels and storage tanks during the scheduling horizon is six (C1-C6). The crude oil types that are allowed to mix in each storage tank are different. In this example, only tank 6 is allowed to mix crude oils C2, C4, and C6. All of the other tanks can store only one specific type of crude oil. The storage tank information is given in Table 9a,b. Instead of calculation of the sulfur concentration of the feed oil, the concentrations of crude oils are calculated directly by their mass balance so that refinery operators know the exact volume of crude oils in each tank. This example consists of two CDUs that can be charged simultaneously from a maximum of two storage tanks. The daily maximum and minimum demands for crude by the CDUs are 8000 and 2000 tons/day for CDU 1 and 3000 and 1000 tons/day for CDU 2. The inventory cost for each storage tank is 0.05 Yuan/day/ton. Unloading cost and sea waiting cost for each vessel are 7000 and 5000 Yuan/day, respectively. The changeover cost is 1000 Yuan for each change. The maximum and minimum flow rates of crude oil from the berths to the storage tanks are 5000 and 0.0 tons/day, respectively. The maximum and minimum flow rates of crude oil from the storage tanks to the CDUs are 5000 and 0.0 tons/ day. The cost of raw materials (C_rawc) and profits of processing the crude oils (C_profitc) are given in Table

6732

Ind. Eng. Chem. Res., Vol. 41, No. 26, 2002

Table 8. Vessel Data for Example 2 amount of crude (tons)

arrival time

least stay time (days)

most stay time (days)

name of oil type 1

ratio of oil type 1

3000 3000 3000 3000 5000 5000 3000 3000 3000 5000 5000 3500 3500 3000 3000 3000 3000

1 1 2 2 3 3 3 3 4 5 5 5 5 6 6 7 7

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

c1 c1 c1 c1 c1 c1 c1 c1 c1 c2 c2 c2 c2 c1 c1 c4 c2

1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.50

v1 v2 v3 v4 v5 v6 v7 v8 v9 v10 v11 v12 v13 v14 v15 v16 v17

name of oil type 2

ratio of oil type 2

unloading cost (Yuan/day)

sea waiting cost (Yuan/day)

0.50

7000 7000 7000 7000 7000 7000 7000 7000 7000 7000 7000 7000 7000 7000 7000 7000 7000

5000 5000 5000 5000 5000 5000 5000 5000 5000 5000 5000 5000 5000 5000 5000 5000 5000

c6

Table 9. Tank Data for Example 2 and Data for the Mixed Oil Permitted in Tank 6 in Example 2 a. Tank Data for Example 2 inventory cost (Yuan/ton/day)

min capacity (tons)

max capacity (tons)

name of oil type 1

0.05 0.05 0.05 0.05 0.05 0.05 0.05

0.0 0.0 0.0 0.0 0.0 0.0 0.0

25 000.0 25 000.0 25 000.0 25 000.0 40 000.0 40 000.0 40 000.0

c1 c2 c3 c1 c1 c2 c3

tank 1 tank 2 tank 3 tank 4 tank 5 tank 6 tank 7

oil type 1 composition scope low high 1.00 1.00 1.00 1.00 1.00 0.00 1.00

oil type 1 initial inventory (tons)

1.00 1.00 1.00 1.00 1.00 1.00 1.00

5000.0 5000.0 7000.0 8000.0 8000.0 10000.0 10000.0

b. Data for the Mixed Oil Permitted in Tank 6 in Example 2 oil type 2 compostion scope low high

name of oil type 2 tank 6

c4

0.00

oil type 2 initial inventory (tons)

name of oil type 3

5000.0

c6

1.00

Table 10. Cost and Price Data for Example 2 C_rawc (Yuan/ton) C_profitc (Yuan/ton) C_prodc (Yuan/ton)

C1

C2

C3

C4

C5

C6

1130 1000 2130

1574 1000 2574

1418 1000 2418

1150 1000 2150

1150 1000 2150

1671 1000 2671

10. The prices of products that can be produced from the crude oils, C_prodc, are also given in that table. The objective function used in example 2 is different from the one listed in section 4.1. The prices of raw materials and products are added to the objective function to account for the actual profitability of the whole production site. The objective of the problem is to maximize the profit rather than minimize the cost. Thus, the objective function of example 2 becomes

∑

maximize profit )

∑

C_rawc

c∈C,t∈Stt

∑

C_prodc

c∈C,t∈Stt

FSCCilct -

i∈I_Cic,l∈L

∑

FVSvit × EVvc -

(v,i)∈VITvit∩V_Cvc

∑ (TLv - TFv + 1) × C_uldv - v∈Sv ∑ (TFv -

v∈Svv

v

TARRv) × C_seav -

∑ ∑ 0.5(VSit + VSi(t-1)) × i∈Si t∈St i

C_ivsi -

t

∑ ∑ ∑ C_setjl × Zjlt j∈J l∈L t∈St

(obj)

t

Model MIP I of example 2 involved 175 discrete variables, 1489 single continuous variables, and 3247 single constraints. Note that XFvt, XLvt, and XVDvd were

oil type 3 compostion scope low high 0.00

1.00

oil type 3 initial inventory (tons) 5000.0

Table 11. Summary of Computational Results for Example 2 model initialization MIP I iteration 1 NLP MIP II iteration 2 NLP MIP II iteration 3 NLP total MINLP a

objective value 105 048 575 80 489 297 94 761 497 94 882 608 95 904 658 95 904 657

CPU nodes iterations time (s) 2202 N/A 156 N/A 118 N/A

intermediate N/A noninteger

50000 257 10231 276 5730 204 50000

247.4 0.55 38.0 0.55 22.8 0.50 309.8 464a

The CPU time for MINLP is the sum of MIP I and MINLP.

defined as SOS1 variables and thus not included in discrete variables (if XFvt, XLvt, and XVDvd were included in discrete variables, the numbers of discrete variables and single continuous variables were 379 and 1285, respectively). To reduce the computing time, the relative optimality gap of the OSL solver was set to 0.1 with an upper iteration limit of 50 000. Table 11 summarizes the solution. When the problem is solved as a MINLP, no feasible solution was found within the iteration limit even after imposing initial values of variables from the MIP I solution. Using the new algorithm, the concentration values calculated by MIP I violate the criteria described in iteration condition I of section 5; thus, the algorithm continues to solve the NLP and MIP II models. The procedure continued for three iterations for convergence, requiring 309 s for the solution. Results of example 2 are summarized in Table 11-14.

Ind. Eng. Chem. Res., Vol. 41, No. 26, 2002 6733 Table 12. Berth Distribution Schedule time berth

1

d1 d2 d3 d4

v1

2 v4

v2

v3

3

4

5

6

7

v8 v6 v7 v5

v9

v12 v10

v14 v10 v13 v11

v15 v16 v17

v11

Table 13. CDUs Charging Schedule time CDU

1

2

3

4

5

6

7

l1

i4 i5 i2 i6

i1 i7 i2 i6

i1 i3 i2 i6

i3 i5 i2 i6

i4 i7 i2

i1 i4 i6

i3 i4 i2

l2

Table 14. Storage Tank Charging Schedule time tank

1

i1 i2 i3 i4 i5 i6 i7

v1, v2

2

3

4

5

6

7 v15

v10, v11, v13 v5, v6, v7 v9 v3, v4 v8

v14 v12

v16, v17

The changeover calculation formulations proposed in this paper are very important in example 2. The changeover formulation for definition D0 (constraint (i) in section 3.3) cannot be applied in example 2 because it introduces so many variables and constraints that the solver needs to allocate a huge amount of memory for it and finally fails to solve the whole model. 7. Conclusions In this paper, the short-term scheduling of crude oil unloading and charging has been studied. This MINLP problem involves bilinear constraints for calculating the crude oil concentration. A new algorithm was proposed to decompose the MINLP into a MIP and a NLP. The total number of discrete decision variables was reduced by decomposing tri-index discrete variables to bi-index discrete variables, hence reducing the solution time. The changeover cost in crude oil charging was efficiently calculated using the systematic method proposed in this paper. Finally, two examples, including an industrial example, are described and solved using the new formulations proposed in this paper. Acknowledgment The authors acknowledge financial support from the Research Grants Council of Hong Kong (Grants 6014/ 99P and 6036/98P) and the National Science Foundation of China (Grant 79931000) and financial and technical support from Jiujiang Refinery (Sinopec), Jiangxi, China. Definition of Sets and Parameters The notations used for the sets and parameters in the mathematical formulation are as follows. Part of the nomenclature in this paper is taken from Lee et al.9 (a) Indices c ) different types of crude oils d ) berths in a docking station i ) crude oil storage tanks

j, j′ ) crude oil charging tanks k ) key components of crude oil l ) crude distillation units t, tt ) time intervals v ) crude vessels (b) Sets C ) types of crude oils D ) number of berths in a docking station I ) crude oil storage tanks J ) crude oil charging tanks K ) key components of crude oil L ) crude distillation units T ) time intervals V ) crude vessels Sii ) number of storage tanks being scheduled Stt ) days being scheduled Svv ) vessels being scheduled V_Cvc ) oil type(s) that vessel v carries I_Cic ) oil type(s) in storage tank i Varrvt ) feasible unloading time for vessel v Vttt,tt ) time tt equal to or greater than time t V_Ivi ) feasible storage tank i for vessel v Vitvit ) feasible time period t for vessel v to unload to storage tank i (c) Parameters C_inbj ) inventory cost of charging tank j per unit time per unit volume C_ivsi ) inventory cost of storage tank i per unit time per unit volume C_profitc ) profit of processing one unit of crude oil c C_prodc ) overall product value of crude oil c C_rawc ) purchasing price of crude oil c C_seav ) sea waiting cost of vessel v per unit time interval C_setjj′l ) changeover cost for transition from crude mix j to j′ in CDU l C_setjl ) changeover cost for crude mix transition from charging tank j to CDU l C_setil ) changeover cost for crude mix transition from storage tank i to CDU l C_uldv ) unloading cost for vessel v per unit time interval DMj ) demand of crude mix j by CDUs during the scheduling horizon DMMINl ) minimum demand of crude oil by CDU l each day DMMAXl ) maximum demand of crude oil by CDU l each day DT ) length of each time interval ESik ) concentration of component k in the crude oil of storage tank i EVvc ) concentration of crude oil c in vessel v EBMINjk ) minimum concentration of component k in a crude mix of charging tank j EBMAXjk ) maximum concentration of component k in a crude mix of charging tank j ESMINic ) minimum concentration of crude oil c in storage tank i ESMAXic ) maximum concentration of crude oil c in storage tank i EB0jk ) initial concentration of component k in a crude mix of charging tank j FVSMINvi ) minimum crude oil transfer rate from vessel v to storage tank i FVSMAXvi ) maximum crude oil transfer rate from vessel v to storage tank i FSBMINij ) minimum crude oil transfer rate from storage tank i to charging tank j FSBMAXij ) maximum crude oil transfer rate from storage tank i to charging tank j

6734

Ind. Eng. Chem. Res., Vol. 41, No. 26, 2002

FSCMINil ) minimum crude oil charging rate from storage tank i to CDU l FSCMAXil ) maximum crude oil charging rate from storage tank i to CDU l FBCMINjl ) minimum crude oil charging rate from charging tank j to CDU l FBCMAXjl ) maximum crude oil charging rate from charging tank j to CDU l MaxVtv ) maximum unloading time of vessel v MinVtv ) minimum unloading time of vessel v ORDt ) sequential order of time interval t TARRv ) arrival time of crude vessel v VV0v ) initial volume of crude oil in crude vessel v VSMINi ) minimum crude oil volume of storage tank i VSMAXi ) maximum crude oil volume of storage tank i VS0i ) initial crude oil volume in storage tank i VSC0ic ) initial volume of crude oil c in storage tank i VBMINj ) minimum mixed crude oil volume of charging tank j VBMAXj ) maximum mixed crude oil volume of charging tank j VB0j ) initial crude oil volume in charging tank j VBK0jk ) initial volume of component k in charging tank j VBC0jc ) initial volume of crude oil c in charging tank j PCDjlt ) parameter corresponding to CDjlt PCDilt ) parameter corresponding to CDilt PXFvt ) parameter corresponding to XFvt PXLvt ) parameter corresponding to XLvt PEBjkt ) parameter corresponding to EBjkt PEBCict ) parameter corresponding to EBCict PXWvt ) parameter corresponding to XWvt PXWvit ) parameter corresponding to XWvit (d) Variables XFvt ) binary variable to denote if vessel v starts unloading at time t XLvt ) binary variable to denote if vessel v completes unloading at time t ILjl ) binary variable to denote if tank j charges CDU l ITjt ) binary variable to denote if tank j charges (to any CDU) at time t ILil ) binary variable to denote if tank i charges CDU l ITit ) binary variable to denote if tank i charges (to any CDU) at time t VIvi ) binary variable to denote if vessel v unloads (at least as one time) to tank i VTvt ) binary variable to denote if vessel v unloads at time t XVDvd ) binary variable to denote if vessel v moors at berth d Cost ) total cost of the refinery during the scheduling horizon CDilt ) 0-1 continuous variable, crude oil mix in storage tank i charges CDU l at time t CDjlt ) 0-1 continuous variable, crude oil mix in charging tank j charges CDU l at time t EBjkt ) concentration of key component k in charging tank j EBCict ) volume fraction of crude oil c in storage tank i at time t FBCjlt ) volumetric flow rate of crude oil mix from charging tank j to CDU l at time t FBCKjlkt ) volumetric flow rate of component k from charging tank j to CDU l at time t FSBijt ) volumetric flow rate of crude oil from storage tank i to charging tank j at time t FSCilt ) volumetric flow rate of crude oil from storage tank i to CDU l at time t FSCCilct ) volumetric flow rate of crude oil c from storage tank i to CDU l at time t

FVSvit ) volumetric flow rate of crude oil from vessel v to storage tank i at time t Profit ) total profit of the refinery during the scheduling horizon TFv ) unloading initiation time of vessel v TLv ) unloading completion time of vessel v VBjt ) volume of mixed oil in charging tank j at time t VBKjkt ) volume of component k in charging tank j at time t VSit ) volume of crude oil in storage tank i at time t VSCict ) volume of crude oil c in storage tank i at time t VVvt ) volume of crude oil in crude vessel v at time t XDvdt ) 0-1, a continuous variable to denote if vessel v moors at berth d at time t XWvt ) 0-1, a continuous variable to denote if vessel v is unloading its crude oil at time t XWvit ) 0-1, a continuous variable to denote if vessel v is unloading its crude oil to storage tank i at time t Yjlt ) binary variable for absolute value calculation Zjj′lt ) 0-1 continuous variable, transition from crude mix j to j′ at time t in CDU l Zilt ) 0-1 continuous variable, transition of storage tanks that charge CDU l at time t

Literature Cited (1) RPMS (Refinery and Petrochemical Modeling System) is a software product developed by Bonner & Moore Management Science, 2727 Allen Parkway, Suite 1200, Houston, TX 77019. (2) PIMS System Reference, version 11.0; AspenTech, Inc.: Cambridge, MA, 1999. (3) Pelham, R.; Pharris, C. Refinery operation and control: a future vision. Hydrocarbon Process. 1996, 75 (7), 89-94. (4) OmniSuitesUnified planning and Scheduling of Crude supply, Refining Product Blending and Distribution, Haverly Systems, Inc., http://www.haverly.com/product.htm. (5) PetroPlan, Ami Consultants Inc. of Texas,. http:// www.amiconsultants.com/index.html. (6) Reklaitis, G. V. Overview of scheduling and planning of batch process operations. Presented at NATO Advanced Study Institutesbatch process system engineering, Antalya, Turkey, 1992. (7) Pinto, J.; Joly, M.; Moro, L. Planning and Scheduling models for refinery operations. Comput. Chem. Eng. 2000, 24 (9-10), 2259-2276. (8) Shah, N. Mathematical Programming Techniques for Crude Oil Scheduling. Comput. Chem. Eng. 1996, 20 (Suppl.), S1227S1232. (9) Lee, H.; Pinto, J.; Grossmann, I. E.; Park, S. Mixed-Integer Programming Model for Refinery Short-Term Scheduling of Crude Oil Unloading with Inventory Management. Ind. Eng. Chem. Res. 1996, 35, 1630-1641. (10) Tu¨rkay, A.; Lee, H.; Pinto, J. M.; Grossmann I. E.; Park, S. CRUDEOILsMILP Model for Refinery Short-Term Scheduling of Crude Oil Unloading with Inventory Management, http:// egon.cheme.cmu.edu/Group/MINLPSOFTWARE/crudeoil.html. (11) Zhang, J.; Zhu, X. X.; Towler, G. P. A Simultaneous Optimization Strategy for Overall Integration in Refinery Planning. Ind. Eng. Chem. Res. 2001, 40, 2640-2653. (12) Sherali, H.; Alameddine, A. A New ReformulationLinearization Technique for Bilinear Programming Problems. J. Global Optim. 1992, 2, 379-410. (13) Quesada, I.; Grossmann, I. E. Global Optimization of Bilinear Process Networks With Multicomponent Flows. Comput. Chem. Eng. 1995, 19, 1219-1242. (14) Hui, C. W.; Gupta, A. A Bi-Index Continuous Time Model for Single-Stage Batch Scheduling with Parallel Units. Ind. Eng. Chem. Res. 2001, 40, 5960-5967. (15) Brooke, A.; Kendrick, D.; Meeraus, A. GAMSsA User’s Guide (Release 2.25); The Scientific Press: San Francisco, CA, 1992.

Received for review February 11, 2002 Accepted October 8, 2002 IE020130B