Hierarchical Approach for Short-Term Scheduling in Refineries

Apr 7, 2007 - ... the operation modes of some multipurpose tanks within aggregate tanks. ... A decision tree based decomposition method for oil refine...
1 downloads 0 Views 140KB Size
Download PDF
3656

Ind. Eng. Chem. Res. 2007, 46, 3656-3668

Hierarchical Approach for Short-Term Scheduling in Refineries Chunpeng Luo and Gang Rong* National Laboratory of Industrial Control Technology, Institute of AdVanced Process Control, Zhejiang UniVersity, Hangzhou 310027, China

This paper focuses on a hierarchical approach with two decision levels for short-term scheduling problems in refineries. The optimization model at the upper level and the heuristics and rules adopted in simulation system at the lower level are presented. The iteration procedure between the upper and the lower levels is also introduced. The optimization model is used to decide sequencing and timing of the operation modes of processing units and pipelines and to determine the quantities of materials consumed/produced by each operation mode of a unit. Only aggregate tanks are used at the upper level, and the simulation system at the lower level uses a heuristic to adjust the operation modes of some multipurpose tanks within aggregate tanks. The iteration procedure between the upper level and the lower level can be used, if needed, to find optimal solution under updated aggregate storage capacities. The simulation system also uses some heuristics and rules to arrange actual tanks to receive/send materials with logic operation constraints (operation rules) on tanks respected. The main purpose of our approach is to reduce the binary variables and the size of the optimization model caused by the representation of multipurpose tanks and operation rules and, at the same time, to determine the detailed material movements with operation constraints respected. Case studies illustrated the efficiency of the proposed approach. 1. Introduction In order to compete successfully in international markets and with global competition, oil refineries are increasingly concerned with improving the planning and scheduling of their operations to achieve better economic performance. Decisions are normally made with the help of mathematic programming technologies such as linear programming (LP), mixed integer linear programming (MILP), and nonlinear programming (NLP). In general, mathematical programming is suited and currently used for longterm production planning, and the availability of LP-based commercial software for refinery production planning, such as RPMS (Refinery and Petrochemical Modeling System) and PIMS (Process Industry Modeling System), has allowed the development of general production plans of the whole refinery.1 However, the short-term scheduling problem is still one of the most challenging problems in operational research due to the complexity of the scheduling operations and the corresponding process models. Many scheduling models for refineries have been proposed from different viewpoints in the literature. Jia and Ierapetritou2 developed a comprehensive mathematical programming model based on a continuous time formulation for the scheduling of oil-refinery operations and decomposed the overall problems into three domains: the crude-oil unloading and blending, the production unit operations, and the product blending and delivery. Reddy and Karimi3 presented the first complete continuous-time (MILP) formulation for the shortterm scheduling of crude-oil unloading operations. Go¨theLundgren and Lundgren4 formulated a mixed integer linear programming model for the scheduling problem with inventory capacities and modes of operation considered and discuss a number of modifications and extensions of the model. Glismann and Gruhn5 also developed an integrated approach to coordinate short-term scheduling of multiproduct blending facilities with nonlinear recipe optimization. Up to now, little literature has dealt with the scheduling problems of the whole refinery. The * To whom correspondence should be addressed. Tel.: +86-57187952248. E-mail: [email protected].

model will be very complicated and difficult to build and solve if all decisions in different levels of the scheduling of an entire refinery are simultaneously considered. The large number of continuous and binary variables, as well as nonlinear constraints involved in the scheduling model, may make MILP/MINLP fail to come up with a feasible production schedule at the right time. Especially the operation rules, such as specifying a minimum amount of time to allow the separation of the brine of crude oils, that just one tank receives the same kind of material at the same time, and that a tank cannot receive and send material at the same time, will lead to a large number of 0-1 variables that would make the model solution unattainable within a reasonable time. Moreover when a model is described with discrete time, the time slot duration must be short enough to honor all the operation rules and the binary variables during the whole scheduling horizon will increase dramatically. Pinto and Joly1 had to define a time slot duration of 15 min for crudeoil scheduling problems, which generated a discrete time mixed integer optimization model with 21 504 binary variables in a 3-4 days horizon. They pointed out that the solution of such a problem is far beyond the capabilities of the current optimization technology. The continuous-time modeling is particularly suitable for scheduling problems with activities ranging from a few minutes to several hours. Because of the possibility of eliminating a major fraction of the inactive event-time interval, the resulting mathematical programming problems are usually of much smaller sizes and require less computational efforts for their solution.6-8 However, due to the variable nature of the timings of the events, it becomes more challenging to model the scheduling processes and the continuous-time approach may lead to mathematical models with more complicated structures compared with their discrete-time counterparts.9 Some other authors have resorted to intelligent simulation systems to deal with scheduling problems. Paolucci10 used a simulation-based DSS to improve the effectiveness of the decisions on allocating crude-oil supply to port and refinery tanks. Chryssolouris and Papakostas11 used an integrated simulation-based approach to evaluate the performance of short-

10.1021/ie061354n CCC: $37.00 © 2007 American Chemical Society Published on Web 04/07/2007

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007 3657

term scheduling with tank farm, inventory, and distillation management. The use of a simulation model to provide decision support can allow the application of the public domain of heuristic knowledge, support what-if analysis, and be able to evaluate the performance of alternative solutions. Simulation systems can also have the ability to quickly react to stochastic disturbances in production by discrete event simulation. But simulation systems rely on the independent variables specified by users, e.g., flows into and out of a unit over time. The simulator can then calculate the dependent variables such as a tank’s holdup and quality at each time interval. The values of those independent variables are mainly depended on scheduler’s experience and cannot ensure optimal schedules. Mathematical programming approaches and simulation approaches all have their own advantages to deal with scheduling problems. In this study, we propose a scheduling system that uses a hierarchical approach with two decision levels. The upper level is a mathematical programming model based on discrete time and the lower level is an intelligent simulation system. The upper level considers the operations related to processing units and pipelines and decides sequencing and timing of the operation modes of these units. The upper level also determines quantities of materials consumed and produced by each operation mode of a unit. The mathematical model does not take the detailed material movements from/to individual tanks into account and only considers aggregate storage capacity for each material. The time slot is 24 h, and many detailed operation rules are ignored in this model. These all will reduce binary variables and the size of the mathematical model. The simulation system automatically controls the detailed material movements according to heuristics and operation rules. Although the upper level is based on discrete-time formulation, we adopt some methods to extend the model and allow an activity of a unit to occur at any time point of a period. In refineries, there are multipurpose tanks that can store several materials but can only store one at a time. This will undoubtedly increase the number of binary variables to represent the different material storage modes of these tanks in the optimization model. We adopt some slack variables to soften the constraints on the limitation of the maximum storage capacities of aggregate tanks in the mathematical model, and aggregate storage capacities are adjusted by changing the operation modes of multipurpose tanks within aggregate tanks according to a heuristic in the simulation system. A multipurpose tank can be adjusted to store another type of material when it is changed to another operation mode, and the storage capacities of the corresponding aggregate tanks will be changed. When the storage capacities of aggregate tanks are changed, it is necessary to recalculate the optimization model to find a new optimal solution. So, there is an iterative link between the upper lever and the lower lever in our approach. Scheduling includes both creating a predictive schedule and adapting an existing schedule when unforeseen events occur.12 The optimization model at the upper level can generate an optimal predictive schedule. And, the “what-if?” analysis of the simulation system at the lower level can enable schedulers to manipulate the schedule and see quickly the impact of changes and to provide anticipation of events by modeling the entire production processes. Plant managers may use it to review daily plant schedules. In this paper, the architecture of the scheduling approach and the termination conditions of our overall algorithm are first described in section 2. In section 3, the optimization model at the upper level of the proposed algorithm is introduced. The heuristics and rules that are used in the simulation system at

the lower level are introduced in section 4. Finally, two case studies are presented in section 5. 2. Architecture of the Scheduling Approach The upper level is a mathematical programming model based on discrete time formulation. The purpose of this mathematical programming model is to determine when, how long, and how much an operation (operation mode) is active and, also, to determine how many materials (raw materials, intermediate materials, and final products) are produced and consumed by the operations. In this model, only the operation modes of processing units, blenders, and pipelines are considered. The details of material movements between all units (including tanks) will be determined at the lower level with operation constraints respected. The optimization model at the upper level can be regarded as a capacitated lot-sizing problem involving material balance constraints, material supply and demand constraints, sequence constraints, capacity constraints, and quality specifications. Only an aggregate tank for each material is considered in this model; but, the maximum storage capacities of aggregate tanks can be violated, and the extra capacity requirements are punished by penalties. In our approach, the lower level is an intelligent simulation system. The variables that are calculated at the upper level are taken as inputs in the simulation system. These variables include the start time and the end time of the operation modes of each unit, the quantity of each material consumed/ produced by every operation mode of a unit, and the extra storage capacity requirements. An important role that the lower level plays is trying to eliminate the extra capacity requirements by adjusting the operation modes of the multipurpose tanks within some aggregate tanks that have redundant capacities. A heuristic is used to adjust the operation modes of these tanks in order to avoid a large number of 0-1 variables in mathematic model. When the simulation system adjusts the operation modes of some multipurpose tanks, the storage capacities of the corresponding aggregate tanks will be changed. The optimization model will be recalculated to find the optimal solution according to the updated storage capacities. So, there may be, if necessary, an iterative link between the upper lever and the lower level. After the termination of the iteration procedure, if an optimal solution is found, the simulation system will use heuristics and knowledge to determine the detailed material movements between all actual units with operation constraints respected. The upper bounds of the storage capacities of aggregate tanks are enlarged in the optimization model, because the storage capacities can be violated in this model. If the enlarged storage capacities cannot ensure a feasible solution in the optimization model, the scheduling problem will not find a feasible solution and the iteration procedure will terminate. If an optimal solution, without extra capacity requirements, is found in the optimization model, the solution will be a feasible solution to the scheduling problem, but it may not be the optimal solution to the scheduling problem. When there is no storage capacity reaching its upper bound in this solution, this indicates that all the current storage capacities are enough to ensure optimal solution and the current solution will be considered as the optimal solution to the scheduling problem; the iteration procedure will then terminate. When there are some storage capacities reaching their upper bound, this indicates that these storage capacities may be bottlenecks to the production processes. Although current storage capacities are not violated in the optimization model, this may be achieved by violating other

3658

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007

materials consumed and produced. The purpose of this model is to determine the start time and end time of each operation and the quantity of materials consumed/produced by each operation of a unit. The objective is to minimize production costs and some kind of penalties with capacity constraints of units and material storages respected and, at the same time, to ensure that material supply, products demands, and product quality specifications are satisfied. The simplified flow sheet of a refinery is illustrated in Figure 2. The figure consists of a set of rectangles/triangles that represent the units/tanks in the plant and “connections” which represent the flow paths between the equipment. The small squares connected to the equipment by short lines are the “ports” where flow enters or leaves the equipment. The pentagons represent the “perimeter” units (equivalent to pipelines) through which materials enter/leave the plant. The ports of perimeters are not shown in the figure for simplification. Supply orders or demand orders can be defined on the perimeters. The optimization model is described as follows. Mass balance at each port:

∑

Qs,u,m,t )

Qs′,u′,m′,s,u,m,t

(s′,u′,m′)∈CONMs,u,m

∀ u ∈ U, ∀ t ∈ TP, ∀ s ∈ ISu, ∀ m ∈ Mu (1)

Figure 1. Iteration procedure of our scheduling approach.

constraints, which are punished with smaller penalties in the objective function. So, it is possible to find better solutions by adjusting storage capacities. If the simulation system at the lower level cannot find redundant storage capacities to eliminate these bottlenecks or the value of the objective function is not better when storage capacities are changed, the iteration procedure will terminate and the current solution is considered as the optimal solution to the scheduling problem. If a solution, with maximum storage capacities of some aggregate tanks violated, is found in the optimization model, the extra capacity requirements must be eliminated by adjusting storage capacities to ensure that the solution is feasible to the scheduling problem. However if the simulation system at the lower level cannot find redundant storage capacities, the maximum storage capacity constraints at the upper level must be changed to hard constraints that cannot be violated. The optimal solution, found by the optimization model with hard constraints, will also be optimal to the scheduling problem, and the iteration procedure terminates. If the optimization model with hard constraints cannot find a feasible solution, the scheduling problem will have no feasible solution and the iteration procedure terminates. When the simulation system finds redundant storage capacities that can be utilized, the upper level will be recalculated according to the updated storage capacities and the iteration procedure continues. It is believed that the number of iterations will be acceptable because of the limit of the number of multipurpose tanks. The iteration procedure of our scheduling approach is showed in Figure 1. 3. Optimization Model at the Upper Level The optimization model is a MILP model based on discrete time. The scheduling horizon is 5-10 days with each period being 24 h. The decision variable yu,m,t represents whether unit u operates in mode m during period t. The operation modes of blenders determine final products, and the operation modes of pipelines determine materials being transferred from/to a plant. The operations modes of other processing units determine the

∑

Qs,u,m,t )

Qs,u,m,s′,u′,m′,t

(s′,u′,m′)∈CONMs,u,m

∀ u ∈ U, ∀ t ∈ TP, ∀ s ∈ OSu, ∀ m ∈ Mu (2)

Equality 1 relates the volumetric flow of inflow port s with every flow that enters that port. Equality 2 relates the volumetric flow of each outflow port s with every flow that leaves that port. Qs′,u′,m′,s,u,m,t represents the flow leaving port s′ of unit u′ in operation mode m′ and entering port s of unit u in operation mode m during period t. Qs,u,m,t represents the flow of port s of unit u in operation mode m during period t. Capacity constraint for each unit (tanks excluded):

QFu,m,t )

∑ Qs,u,m,t

∀ u ∈ NTU, ∀ t ∈ TP, ∀ m ∈ Mu (3)

s∈ISu

up yu,m,tQFlow u e QFu,m,t e yu,m,tQFu ∀ u ∈ NTU, ∀ m ∈ Mu, ∀ t ∈ TP (4)

∑

yu,m,t e 1 ∀ u ∈ NTU, ∀ t ∈ TP

(5)

m∈Mu

The charge size of each unit in operation mode m during period t is the sum of the volumetric flow of each inflow port of that unit, as showed in equality 3. Constraint 4 specifies that the minimum and maximum volumetric capacity must be satisfied if unit u operates in mode m during period t. If unit u does not operate in mode m during period t (yu,m,t ) 0), the charge size of the unit in mode m during period t will be zero. Constraint 5 specifies that these units can at most be in one operation mode during period t and may be idle during period t. Variety of the charge size of each unit (tanks excluded):

TQFu,t )

∑

QFu,m,t ∀ u ∈ NTU, ∀ t ∈ TP

(6)

m∈Mu

TQFu,t - TQFu,t-1 e ∆TQFu,t ∀ u ∈ NTU, ∀ t ∈ TP (7) TQFu,t-1 - TQFu,t e ∆TQFu,t ∀ u ∈ NTU, ∀ t ∈ TP (8) In refineries, schedulers are concerned about the fluctuation of

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007 3659

Figure 2. Simplified flow sheet of a refinery.

the charge size of each unit. Equality 6 calculates the charge size of unit u during period t. ∆TQFu,t is the difference of the charge size between two adjacent periods. Constraints 7 and 8 are used to avoid the negative value of ∆TQFu,t because the objective function is to minimize the value of ∆TQFu,t. The value of ∆TQFu,t can become either (TQFu,t - TQFu,t-1) or (TQFu,t-1 - TQFu,t) depending on which is positive. Volume of each outflow port in one operation mode of a unit (tanks excluded):

Qs′,u,m,t ) QFu,m,tYIEDs′,u,m ∀ s′ ∈ OSu, ∀ u ∈ NTU, ∀ m ∈ Mu, ∀ t ∈ TP (9) Yield expressions in our approach are based on a standard value, YIEDs,u.m, which is determined over average values obtained from plant data. In different modes of a processing unit, the materials produced have different yields. Material balance equation for each tank:

INVu,m,t ) INVu,m,t-1 + Qs,u,m,t - Qs′,u,m,t ∀ u ∈ TK, ∀ s ∈ ISu, ∀ s′ ∈ OSu, ∀ t ∈ TP (10) The inventory in each tank at the end of period t is equal to the inventory at the end of period t - 1 plus the amount of the flows entering the unit during period t minus the amount of flows leaving the unit during period t. The tanks used in this model are aggregate tanks. Each aggregate tank has only one operation mode, one inflow port, and one outflow port. The storage capacity of each aggregate tank can change by adjusting the operation mode of the multipurpose tanks included in them. Storage capacity constraint:

Material movements control:

Qs′u′,m′,s,u,m,t W TFs′,u′,s,u,tFLOWUPs′,u′,s,u ∀ (s′,u′,s,u) ∈ CONN, ∀ t ∈ TP (13) Sometimes the materials from one port of a unit can move to some different destinations, but they can only go to one place at a time and the different movements may have different priorities. In order to satisfy this kind of requirement and ensure that the higher priority movement is always first selected, different penalties are imposed on different movements in the objective function (the sixth term in the objective function) and the higher priority movement has the smaller penalty. But, when some destinations (tanks in particular) are close to the minimum or maximum capacity, the changeover between different movements may be very frequent. So in the objective function, the term ∑t∑s′,u′ (∑s,u∈CONs′,u′ TFs′,u′,s,u,tPENLFs′,u′,s,u) is used to minimize all the movements during the whole scheduling horizon. Binary variable TFs′,u′,s,u,t represents whether or not a movement from port s′ of unit u′ to port s of unit u is selected during period t. Constraint 13 ensures that the volumetric flow will be zero when the corresponding movement is not selected. Product property constraints:

PROs′′,u,m,p,t )

∑

∑

QFRs′,u′,m′,s,u,m,tPROs′,u′,m′,p

s∈ISu (s′,u′,m′)∈CONMs,u,m

∀ s′′∈ OSu, ∀ u ∈ BLEN, ∀ m ∈ Mu, ∀ t ∈ TP, ∀ p ∈ P (14) Qs′,u′,m′,s,u,m,t ) QFu,m,tQFRs′,u′,m′,s,u,m,t ∀ s′′ ∈ OSu, ∀ u ∈ BLEN, ∀ m ∈ Mu, ∀ t ∈ TP, ∀ p ∈ P (15)

Softened product storage capacity constraint:

low up PROs′′,u,m,p,t e PROs′′,u,m,p,t e PROs′′,u,m,p,t ∀ s′′ ∈ OSu, ∀ u ∈ BLEN, ∀ m ∈ Mu, ∀ t ∈ TP, ∀ p ∈ P (16)

+ INVu,m,t e INVup u + INVu,m,t ∀ u ∈ VTA, ∀ t ∈ TP (12)

low PROs′′,u,m,p,t QFu,m,t e

up INVlow ∀ u ∈ TA, ∀ t ∈ TP (11) u e INVu,m,t e INVu

Constraint 11 specifies that the minimum and maximum capacity of each aggregate tank must be satisfied if there are no other storage capacities that can be adjusted to extend the current storage capacity. Otherwise, the maximum storage capacity can be violated as showed by constraint 12.

∑

∑

s∈ISu (s′,u′,m′)∈CONMs,u,m up Qs′,u′,m′,s,u,m,tPROs′,u′,m′,p e PROs′′,u,m,p,t QFu,m,t

∀ s′′ ∈ OSu, ∀ u ∈ BLEN, ∀ m ∈ Mu, ∀ t ∈ TP, ∀ p ∈ P (17) These constraints are used to blenders only. The properties of

3660

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007

the product leaving the outflow port (each blender only has one outflow port) can be directly predicted by using a volumetric average as shown in eq 14. QFRs′,u′,m′,s,u,m,t is the volume fraction of flow Qs′,u′,m′,s,u,m,t in the flows that enter blender u during period t. The volume fraction variable QFRs′,u′,m′,s,u,m,t is linked to the volumetric flow variable Qs′,u′,m′,s,u,m,t and the charge size variable QFu,m,t through equality 15. Constraint 16 specifies that property p of the flow leaving the outflow port must satisfy the minimum and maximum specifications. Constraint 16 makes the model bilinear, and in order to preserve the linearity of the model, constraint 16 can be expressed in an alternative way by multiplying it by QFu,m,t. Then, constraints 14-16 can be replaced by constraint 17. For simplification, in this paper, we assume that the properties of the flows entering the blenders are constant. Generally speaking, one operation mode defined on a blender represents one grade of product produced in this mode or one kind of recipe for the product. The blending index for each of the different nonlinear qualities is used to convert each of the nonlinear qualities to a linear index that may be blended linearly by volume. Product supply and demand constraint: low DDu,m,t1,t2 e

up Qs,u,m,t e DDu,m,t1,t2 ∑ t1etet2

∀ s ∈ ISu, ∀ u ∈ PERI, ∀ m ∈ Mu (18)

low e DDu,m.t1,t2

up ∑ Qs,u,m,t e DDu,m.t1,t2 t1etet2

∀ s ∈ OSu, ∀ u ∈ PERI, ∀ m ∈ Mu (19)

The orders defined on the perimeter units specify the minimum and maximum quantities of demands/supplies and the time horizon of demands/supplies. Constraints 18 and 19 ensure that the sum of the volumetric flows entering/leaving perimeter u during the time horizon satisfy the minimum and maximum demands/supplies. An inflow perimeter has only one outflow port that supplies materials to the plant, and an outflow perimeter has only one inflow port that receives products from the plant. The different operation modes defined on perimeters represent the different grades of products or different raw materials transferred in these modes. Constraint 18 is for outflow perimeter, and constraint 19 is for inflow perimeter. Changeover between operation modes:

yu,m,t-1 + yu,m′,t e SMMu,m,m′,t + 1 ∀ u ∈ NTU, ∀ (m, m′) ∈ Mu, ∀ t ∈ TP, m * m′ (20) yu,m,t - yu,m,t-1 e SMu,m,t ∀ u ∈ NTU, ∀ m ∈ Mu, ∀ t ∈ TP (21) Whenever there is a change of the operation mode from m to m′, constraint 20 forces the corresponding variable SMMu,m,m′,t to take the value of 1. And whenever there is a change of the operation mode, constraint 21 forces the corresponding variable SMu,m,t to be 1. The change of operation mode is punished with penalties in the objective function to avoid the frequent fluctuation of production. When the penalties depend on the sequence of two different modes changed between, constraint 20 is used. Otherwise, constraint 21 is used. The positive penalties for SMMu,m,m′,t and SMu,m,t will always force them to take the value 0 when there is no changeover of operation modes and to take the value 1 when these are changes. So, the integer requirement on these two variables can be relaxed, and they become 0-1 continuous variables.

Target for the charge size of the processing units:

∑t TQFu,t ) TGQFu + TGQF+u - TGQF-u

∀ u ∈ NTU, ∀ m ∈ Mu (22)

Targeted values of the inventory in the storages: INVu,m,t ) TGINVu + TGINV+ u - TGINVu ∀ u ∈ TK, ∀ t ∈ ETP (23)

In order to improve the performance of the plant or ensure the completion of the monthly production plan, schedulers in refineries may want to control the charge size of some important processing units during the scheduling horizon and control the inventory at the end of the scheduling horizon. Constraints 22 and 23 are used to ensure that these targets are achieved as possible, and the deviations from these targets are punished with penalties in the objective function. Objective function:

min

∑ ∑ ∑

∑

s′∈ISu′ u′∈U m′∈Mu′ (s, u, m)∈CONMs′,u′,m′

∑t

CQs,u,mQs,u,m,s′,u′,m′,t +

∑u m∈NTU ∑ m′∈M ∑ ∑t PLSMMu,m,m′SMMu,m,m′,t + u

∑u m∈NTU ∑ ∑t ( ∑ ∑t ∑ s′,u′ s,u∈CON

PLSMu,mSMu,m,t +

∑ ∑t PLTu∆TQFu,t +

u∈NTU

TFs′,u′,s,u,tPENLFs′,u′,s,u) + s′,u′

∑ ∑ ∑

∑

s′∈OSu′ u′∈U m′∈Mu′ (s,u,m)∈CONMs′,u′,m′

∑t

PLQs,u,m,s′,u′,m′Qs,u,m,s′,u′,m′,t +

+ + PLI+ ∑ ∑ u INVu,m,t + ∑ PLTGQu(TGQFu + u∈VTA t u∈NTU + ) + PLTGIN TGQF∑ u u(TGINVu + u∈ΤΚ TGINVu ) + ∑ ∑ CINVuINVu,m,t t u∈ΤΚ

The main objective of the scheduling problem is to minimize the cost of products and all kind of penalties. The first term in the objective function is used to reduce the product cost. In this paper, we only consider the cost of the products produced by the blenders that are the sum of the costs of all the flows entering into the blenders. The second and third terms in the objective function use penalties to avoid frequent changeover of operation modes. The fourth term uses penalties to minimize the fluctuation of the charge size of each unit. The fifth and sixth terms are used to ensure that higher priority movement is always first selected and the number of all movements is minimum. The seventh term is used to force the variable + INVu,m,t to be zero when the current storage capacities are sufficient. The eighth term uses penalties to minimize the deviations of the charge size of processing units from targets, and the deviations of the inventories at the end of the scheduling horizon from targets are minimized by the ninth term. The tenth term is used to minimize the total inventory levels. This model is based on discrete time with a time period of 24 h, and this means that each operation mode can only start at the beginning of a period and end at the end of a period. In order to extend the model and allow the start and end of each operation mode at any point in a period, we use the methods

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007 3661 + mentioned by Go¨the-Lundgren.4 The term yu,m,t-1 means that operation mode m of unit u is extended from period t - 1 to t, represents the delay of the start of an operation mode and yu,m,t m of unit u. In refineries, when a unit is changed to another operation mode, the unit usually needs time to stabilize under the new operation conditions. The quality of the products will therefore fluctuate for some time after a changeover. The policies of how to deal with these products are depended on the sequence of the operation modes changed between. For example, when a unit is changed from 0# diesel mode to -10# diesel mode, the diesel produced can only move to tanks that store 0# diesel before the -10# diesel mode is stabilized. When the unit is changed form -10# diesel mode to 0# diesel mode, the diesel produced can first move to tanks that store -10# diesel for some time and then move to tanks that store 0# diesel immediately after the quality of the products produced cannot satisfy the specification of -10# diesel. So, we adopt another variable tu,m,t to represent the time that operation mode m can continue to stay stable during period t when the operation mode changes from m to another. The value of this variable depends on the value of variable SMMu,m,m′,t and parameter tu,m,m′. Here, tu,m,m′ represents the time that operation mode m can continue to stay stable when the operation mode changes from m to m′. Then, we modify constraint 4 to the following constraints:

+ (yu,m,t + yu,m,t-1 + tu,m,t - yu,m,t )QFlow u e QFu,m,t e (yu,m,t + + + tu,m,t yu,m,t )QFup yu,m,t-1 u ∀ u ∈ NTU, ∀ m ∈ Mu, ∀ t ∈ TP (24)

tu,m,t )

∑

SMMu,m,m′,ttu,m,m′/24

m′∈Mu

+ e yu,m,t-1

∀ u ∈ NTU, ∀ m ∈ Mu, ∀ t ∈ TP (25)

∑

SMMu,m,m′,t ∀ u ∈ Uu, ∀ m ∈ Mu, ∀ t ∈ TP

m′∈Mu

yu,m,t e

∑

(26)

SMMu,m′,m,t ∀ u ∈ NTU, ∀ m ∈ Mu, ∀ t ∈ TP

m′∈Mu

∑

m∈Mu

+ (yu,m,t-1 + tu,m,t)

(27) )

∑

yu,m,t

∀ u ∈ NTU, ∀ t ∈ TP

m∈Mu

(28)

Equality 25 calculates the time of operation mode m that can continue to stay stable during period t when there is a change of operation mode from m to another in that period. The unit of tu,m,m′ is hour, and the time period in our model is 1 day (24 h). So, tu,m,m′ is divided by 24 in equality 25. Constraint 26 specifies that operation mode m can be extended to the next period only if there is a change of operation mode from m to another mode in the next period. Constraint 27 ensures that operation mode m can be delayed only if this mode is started during period t. Constraint 28 ensures that the extended time equals the delayed time. 4. Simulation System at the Lower Level At the lower level, the simulation system serves five objectives: (1) to try to adjust the operation modes of multipurpose tanks within some aggregate tanks to change storage capacities of other aggregate tanks, (2) to construct a schedule including detailed material movements between units (process units, pipelines, and actual tanks), (3) to evaluate the

performance of alternative schedules obtained by the experience of schedulers, (4) to generate detailed scheduling orders to workshops, and (5) to evaluate the effects on production when production fluctuation and unpredicted events occur. 4.1. Adjustment of the Operation Modes of Some Actual Tanks. If a solution, with the maximum storage capacities of some aggregate tanks violated, is found in the optimization + in the solution), the first thing model (there is nonzero INVu,m,t this simulation system will do is to try to eliminate these extra + capacity requirements. The value of arg maxt∈TP INVu,m,t represents the degree to which the current maximum storage capacity of aggregate tank u cannot satisfy the demands of production. The reason we use this variable at the upper level is that its value can be eliminated by some polices. The policy we take is a heuristic trying to adjust the operation modes of some multipurpose tanks within other aggregate tanks that the storage capacities are redundant. In the simulation system, the heuristic sequences aggregate tanks in descending order on the values of arg maxt∈TP + INVu,m,t . The heuristic first selects aggregate tank u that has + and tries to adjust the maximum value of arg maxt∈TP INVu,m,t the operation modes of the multipurpose tanks within other aggregate tanks to eliminate the extra storage capacity requirement of this aggregate tank. If t is the first period that the maximum storage capacity of aggregate tank u is violated, multipurpose tank u′′ that can be used to extend the maximum capacity of u must satisfy the following requirements:

(1) It must at least be idle during period t - 1. In order to ensure enough time for multipurpose tank u′′ to be cleaned for another use during period t, tank u′′ should at least be idle during period t - 1. Whether tank u′′ can be idle during period t - 1 is judged by the following formulations. up up INVu′′ e INVu′,t′ - INVu′,m′,t′ ∀ u′′ ∈ u, t′ ) t - 2, and t′ ) t - 1 (29)

INIINVu′′ -

∑

low Qs′,u′,m′,t e INVu′′

t t - 1 (31) INVu′′

Formulation 31 specifies that the remaining storage capacity of aggregate tank u′ from period t to the end of the scheduling

3662

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007

horizon must still be sufficient after u′′ is removed from u′ during period t. Among the tanks that can be selected for aggregate tank u, if there are some tanks whose maximum storage capacity + , the tank with the smallest exceeds arg maxt∈TP INVu,m,t maximum storage capacity is selected. If there are no tanks with + maximum capacity bigger than arg maxt∈TP INVu,m,t , the tank with the biggest maximum storage capacity is first selected, and then, the tank with the biggest maximum storage capacity among the remaining tanks is selected until the sum of the maximum + . If the capacities selected is more than arg maxt∈TP INVu,m,t sum of the capacities of all the candidate tanks is less than arg + + maxt∈TP INVu,m,t , this means that arg maxt∈TP INVu,m,t cannot be fully eliminated and all the candidate tanks must be selected to change their operation modes. Once tank u′′ within aggregate tank u′ is selected to extend the maximum storage capacity of aggregate tank u, the heuristic will modify the values of INVup u up and INVu′ according to the storage capacity of u′′. If a tank in u′′ within aggregate tank u′ is adjusted for another use during period t, the maximum storage capacity of u′′ will be removed from u′ from period t to the end of the scheduling horizon. If there are no tanks that can be used to extend the storage capacity aggregate tank u, the heuristic will try to deal with the next + aggregate tank with smaller arg maxt∈TP INVu,m,t . When the adjustment for aggregate tank u finishes, the heuristic will select the next aggregate tank with smaller arg + to deal with until all the aggregate tanks with maxt∈TP INVu,m,t + are dealt with. At last, all the nonzero values of INVu′,m,t adjusted maximum storage capacities will feed back to the upper level. The storage capacity constraints on the aggregate tanks, whose maximum storage capacities cannot be extended any more, will be changed to hard constraints as constraint 11. Then, the optimization model will be recalculated in order to find the optimal solution according to the new conditions. The pseudocode for this heuristic is showed in Appendix A. 4.2. Determination of the Material Movements. When the iteration procedure finishes and an optimal solution is found, the simulation system will determine the detailed material movements between units about when, where, and how much the materials move. General speaking, the material movements from/to tanks can be well controlled by some heuristics. These heuristics reflect the experience of schedulers and the actual operations adopted by refineries. We found four policies to select tanks in refineries. The objectives of these four policies are the following:

(1) Minimizing the number of tanks used to stock a material

(2) Maximizing the quantity of materials in each tank that is not completely filled (3) Minimizing the number of tanks used to feed processing units (4) Minimizing the quantity of materials in each tank that is not completely emptied Five simple heuristics are used in our system for these objectives. The second and third heuristics are based on the work of Paolucci:10 Multipurpose Tank to Be Adjusted for Another Use First (MUL_OBJ_FRS) aims at emptying the multipurpose tank before its operation mode is adjusted:

NEXSELt ) multipurpose tank u ∀ u∈ u′, ∀ u′ ∈ TK

Maximum Available Capacity First (MAX_CAP_FRS) aims at using the minimum number of tanks to stock a material:

NEXSELt ) arg max[CAPu - INVu,t]

∀ u ∈ u′, ∀ u′ ∈ TK

Minimum Available Capacity First (MIN_CAP_FRS) gives priority to maximize the quantity of materials in each tank that is not completely filled, i.e., it aims at reducing the ‘‘fragmentation” of the free capacity of the storage system:

NEXSELt ) arg min[CAPu - INVu,t] ∀ u ∈ u′, ∀ u′ ∈ TK Maximum Available Volume First (MAX_VOL_FRS) aims at using the minimum number of tanks to feed units:

NEXSELt ) arg max[INVu,t] ∀ u ∈ u′, ∀ u′ ∈ TK Minimum Volume First (MIN_VOL_FRS) gives priority to minimize the quantity of materials in each tank that is not completely drawn, i.e., it aims at minimizing the quantity of materials in each tank:

NEXSELt ) arg min[INVu,t] ∀ u ∈ u′, ∀ u′ ∈ TK These five heuristics determine the next tank to be selected when the present tank is full or empty. The multipurpose tanks that will be adjusted for another use are always first selected to supply materials, so as to ensure that the inventory in these tanks are completely drawn before their operation modes are adjusted. The heuristics mainly determine that the next tank may be used to receive or supply materials. When a tank is selected by the heuristics to execute the current task, the simulation system will check if the new operation on this tank can satisfy the logic constraints (operation rules) on it. For example, according to the heuristic MAX_CAP_FRS, a tank with minimal inventory is the next tank to receive materials, but this tank has drawn materials out of it at that time and this tank has a logic constraint that it cannot be filled and drawn at the same time. Then, the simulation system has to select another tank to receive materials according the same heuristic. The operation rules in refineries can be represented as if (condition) wait (time) then (action) rules. These rules reflect the relationship between two adjacent operations of a unit and decide the start time of the next states and operations according to the current states and operations. Based on these rules, the simulation system can control the detailed start time and end time of each material movement. For example, the fill-draw-delay constraint can be stated as the following: IF (Tank1.state ) FillEnd) WAIT (1 h) THEN (Tank.state ) DrawBegin). This rule shows that the unit Tank1 can only supply materials 1 h later after the end of receiving materials. In order to determine the detailed start and end time point of an operation of an actual tank, the simulation system also need to know the flow rate of the materials entering into or leaving the aggregate tank that the actual tank belongs to. Because we allow that the operation of a unit can start and end at any time point in each period, this will result in unequal flow rates of each material in a period. So in the simulation system, we divide a period into several small time slots and make sure that the flow rate of a material is constant during a time slot. These time slots depend on the time points of the change of operation modes of each unit. For each aggregate tank, the simulation system first sequences the time point on the ascending values + + tu,m,t - yu,m,t in each period, where u is the of yu,m,t + yu,m,t-1 unit that sends or receives material to/from this aggregate tank.

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007 3663 Table 1. Capacity Data and Initial Inventory

Figure 3. Gantt chart for some materials and tanks.

Then, several time slots are obtained according to each of two adjacent time points. The simulation system calculates the flow rate of each flow from and to every aggregate tank in each period according to the following two formulations:

Fs,u,m,s′,u′,m′,t ) Qs,u,m,s′,u′,m′,t/(yu,m,t + + + tu,m,t - yu,m,t )/24 ∀ s ∈ OSu, ∀ u ∈ NTU, yu,m,t-1 ∀ m ∈ Mu, ∀ s′ ∈ ISu′, ∀ u′ ∈ TK, ∀ m ∈ Mu′, ∀ t ∈ TP (32)

FOs′,u′,m′,s,u,m,t ) Qs′,u′,m′,s,u,m,t/(yu,m,t + + + tu,m,t - yu,m,t )/24 ∀ s ∈ ISu, ∀ u ∈ NTU, yu,m,t-1 ∀ m ∈ Mu, ∀ s′ ∈ OSu′, ∀ u′ ∈ TK, ∀ m ∈ Mu′, ∀ t ∈ TP (33)

where Fs,u,m,s′,u′,m′,t and FOs′,u′,m′,s,u,m,t are the flow rate of each flow entering and leaving aggregate tank u′, respectively. The total flow rates in each slot are the sum of all the flow rates of the flows entering/leaving the aggregate tank in this slot. Then, the simulation system will calculate the start time and end time of the operations on each actual tank according to the heuristics, as well as the total flow rates. The simulation calculated six time points: the start time of filling a tank u-fill-start, the end time of filling of a tank u-fill-end, the start time of drawing from a tank u-draw-start, the end time of drawing from a tank u-draw-end, the start time of the fullness of a tank u-full, and the start time of the emptying of a tank u-empty. The pseudocode for this procedure is showed in Appendix B. These time points are important for the simulation system to control the logic operations according to the heuristics and rules. And, these time points are also important to determine the detailed start time and end times of each material that moves from/to a tank. Let us take the data in Figure 3 for example. Material s from two different units moves to tanks that are storing material s. Unit unit1 produces material s from time t2 to t5, and unit unit2 produces material s from time t3 to t6. Tank1 and tank2 receive material s from t1 to t4 and t4 to t7, respectively. Then, we will be clear that material s from unit1 moves to tank1 from time t2 to t4 and moves to tank2 from time t4 to t5. Also, it can be concluded that material s from unit2 moves to tank1 from time t3 to t4 and to tank2 from t4 to t6. When the timing of material movements, the timing of operation modes of units, and the quantity of each flow are determined by this simulation system, then detailed scheduling orders can be generated to instruct the operations of operators at workshops. This system can also be used to evaluate the effects on production when production fluctuation and unpredicted events occur. The user can view the simulation results by modifying the unit states, material inventories, and supply and demand information according to actual production situations. When

unit

LB

UP

INIVF

unit

LB

UP

TKNA TKFCC1 TKFCC2 TKFCC3 TKVR TKLNA TKHOGO TKMTBE TKHT TKG90 TKG93 TKG95 LIH P1 P2

0 1600 3400 0 0 0 0 800 5000 4536 4174 3000 0 0 0

2400 5400 15300 40500 10000 1000 20000 4000 50000 32000 28024 24000 6400 600 600

2076 5229 9528.5 20655 9200 900 17500 3952 17300 21733.3 17887.75 10000

DAQ CDU1 CDU2 FCC1 FCC2 FCC3 HT RAU CKU BLEN1 BLEN2 KER DO MTBE

0 312 350 40 45 100 180 20 125 100 100 0 0 0

600 375 562 100 125 208 282 68.49 250 400 400 600 600 600

INIVF

Table 2. Property and Cost Data of Materials Entering the Blenders unit

port

proper

cost

TKLNA TKHOGO TKMTBE FCC1 FCC2 FCC3

OUT OUT OUT GASO GASO GASO

80 95 110 88 87 86

100 118.75 137.5 110 108.75 107.5

infeasible information is found by this system, then the user can adjust unit operation modes and material movement paths according to his experience or resort to the optimization model at the upper level to find the optimal solution under the updated production situations. 5. Case Studies 5.1. Case 1. Case 1 is studied based on the flow sheet shown in Figure 2. The flow sheet is composed of an atmospheric distillation column (CDU1), an atmospheric distillation column combined with a vacuum distillation column (CDU2), three fluidized catalytic cracking units (FCC1, FCC2, and FCC3), a delaying coke unit (CKU), a hydrotreatment unit (HT), a reforming unit (RAU), two blenders (BLEN1 and BLEN2), twelve aggregate tanks, three inflow perimeters, and four outflow perimeters. Atmospheric distillation fractionates crude oil into the following hydrocarbon streams naphtha (NA), kerosene (KER), light diesel (LD), and atmospheric residue (AR). The vacuum distillation column fractionates the AR stream into two streams: vacuum gas oil (VGO) and vacuum residue (VR). The FCC unit produces a gasoline component and diesel oil. CKU produces a gasoline component, diesel oil, and wax oil. HT improves diesel oil quality by reducing the sulfur content. The 90# gasoline, 93# gasoline, 95# gasoline can be blended by each blender. Two different crude oils DAQ and LIH are supplied to CDU1 and CDU2, respectively. The minimum capacity (LB) and maximum capacity (UP) for each unit and the initial inventory (INIVF) for tanks are shown in Table 1. The unit of each minimum capacity, maximum capacity, and initial inventory for the tanks is cubic meters. Otherwise, the unit is cubic meters per hour. The actual tanks in TKG90, TKG93, and TKG95 are all multipurpose tanks. The octane value and cost of each material entering into the two blenders are shown in Table 2. Port OUT represents the outflow port of a tank. The unit and the corresponding port represent what material is transferred. For example, unit TKLNA and port OUT represent that material LNA (light naphtha) leaves from port OUT of aggregate tank TKLNA. Also, the combination unit FCC1 and port GASO represents the case where material GASO (gasoline component) leaves port GASO of unit FCC1.

3664

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007

Table 3. Different Penalties for Different Movements srunit

srport

deunit

deport

penalty

CDU1 CDU1 CDU1 CDU2 CDU2 CDU2

AR AR AR VGO VGO VGO

TKFCC1 TKFCC2 TKFCC3 TKFCC1 TKFCC2 TKFCC3

IN IN IN IN IN IN

2 1 0 0 1 2

Table 4. Specifications on Octane Value of Gasoline Products mode

LB

UP

G95 G93 G90

94.85 92.85 89.85

95.25 93.25 90.25

Table 5. Orders Defined on Each Perimeter unit

mode

ST

ED

LB (m3)

UP (m3)

DAQ LIH P1 P2 P2 KER P3 MTBE

DAQ LIH G90 G93 G93 KER DO MTBE

1 1 1 1 4 1 1 1

5 5 5 3 5 5 5 5

0 0 28800 23760 15840 0 0 0

60000 60000 28800 23760 15840 60000 60000 6000

The materials that leave from port AR of unit CDU1 and port VGO of unit CDU2 can go to either aggregate tanks of TKFCC1, TKFCC2, and TKFCC3. But, each material can only move to one tank at a time and the different movements have different priorities. The penalties for different movements are shown in Table 3. Port IN represents the inflow port of a unit. A smaller penalty implies a higher priority. For example, the movement from port AR of unit CDU1 to port IN of unit TKFCC1 has the highest priority (smallest penalty) for the material leaving port AR of unit CDU1. The two blenders (BLEN1 and BLEN2) and perimeters P1 and P2 all have three operation modes: G90 (90# gasoline mode), G93 (93# gasoline mode), and G95 (95# gasoline mode). The value of tu,m,m′ for the blenders is 0.5 h (30 min), and that of tu,m,m′ for P1 and P2 is 0. They can only transfer material to/from one aggregate tank in one operation mode. For example, BLEN1 in operation mode G90 can only send material to TKG90, and P2 in operation mode G93 can only receive material from TKG93.¢¢ The specifications on octane value of the products leaving the blenders in different modes are shown in Table 4. A scheduling horizon of 5 days is used for this case study. The penalty for inventory CINVu and the penalty for the extra storage capacity requirement are 1. The penalties for different flows are shown in Table 3. The other penalties are all 10 000. The targets for the charge size of CDU1 and CDU2 during the entire scheduling horizon are 45 000 m3 and 67 440 m3, respectively, and the targets for the inventories of TKG90, TKG93, and TKG95 at the end of the scheduling horizon are 6936, 6974, and 5400 m3, respectively. The orders defined on each perimeter unit are shown in Table 5. Each order specifies the minimum and maximum quantities of demands/supplies and the time horizon of demands/supplies. The unit and mode shown in Table 5 represent the materials transferred by the unit in that mode. ST and ED represent the start time period and end time period in which the materials can be transferred. LB and UP represent the minimum and maximum quantities that can be transferred between period ST and DT. The MILP models were implemented in LINGO 9. The calculations were performed on a Celeron(R) 2.40 GHz/512Mb RAM platform. The model involves 175 discrete variables, 5924

Table 6. Inventory in TKVR and the Extra Storage Capacity Requirement in Each Period period

inventory (m3)

extra (m3)

1 2 3 4 5

9593.312 9986.624 10379.94 10773.25 11166.56

0 0 379.936 773.248 1166.56

Table 7. Inventory in TKG90, TKG93, and TKG95 at the End of Each Period period unit

1 (m3)

2 (m3)

3 (m3)

4 (m3)

5 (m3)

TKG90 TKG93 TKG95

18773.84 14705 10000

15814.38 11522.25 10000

12854.92 6574 11765.5

9895.46 6574 8582.75

6936 6574 5400

Table 8. Multipurpose Tanks in Aggregate Tanks aggregate tank

actual tank

LB (m3)

UP (m3)

INIINV (m3)

TKG90 TKG90 TKG90 TKG90 TKG93 TKG93 TKG93 TKG95 TKG95 TKG95

Tank1 Tank2 Tank3 Tank4 Tank5 Tank6 Tank7 Tank8 Tank9 Tank10

1156 1110 1100 1170 1384 1395 1395 1000 1000 1000

8000 7900 8000 8100 9316 9308 9400 8000 8000 8000

7700 7763 3000 3270.3 8100 7345 2442.75 5000 3000 2000

continuous variables, and 5867 single constraints. The CPU time required for the solution of the MILP model is 58 s, and the value of the objective function is 5 195 766. The simulation system finds that only the maximum storage capacity of aggregate TKVR is violated during period 3, and there are no other bottleneck tanks. The inventory in TKVR and the extra capacity requirement are shown in Table 6. The flows leaving TKG90, TKG93, and TKG95 during period 1 are 5760, 7920, and 0 m3, respectively, and the inventory in these tanks at the end of each period is shown in Table 7. The data of the actual multipurpose tanks in TKG90, TKG93, and TKG95 are shown in Table 8. According to formulations 29, 30, and 31, Tank3, Tank4, and all the tanks in TKG93 can be used to extend the maximum storage capacity of TKVR. The maximum storage capacity of each of these tanks is bigger than the maximum extra capacity requirement of TKVR. Because Tank3 has the smallest maximum storage capacity, it is selected to extend the maximum capacity of TKVR. By this adjustment, the maximum storage capacity of TKVR and TKG90 change to 18 000 and 24 000 m3, respectively. The CPU time required for the solution of the updated MILP model is 61 s, and the value of the objective function is 5 193 446. There are no tanks with the maximum capacities violated found in the new solution, and the iteration procedure terminates. It is found that the inventory in each aggregate tank at the end of each period does not take any change; it may be caused by the tenth term in the objective function used to minimize the total inventories. In the new solution, only the operation modes of Blen2 change from G93 + for Blen2 in to G95 during period 3. The value of yu,m,t-1 operation mode G93 during period 3 is 0.6458, and this means that operation mode G93 is extended for about 15.5 h (0.6458 × 24). The value of yu,m,t-1 for Blen2 in operation mode G95 during period 3 is 0.6666, and this means that operation mode G95 is delayed for about 16 h (0.6666 × 24). The delayed time is equal to the extended time plus 0.5 h, that is, the time for G93 to be still stable. The amount of products produced by mode G93 during period 3 is 2971.75 m3, and the

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007 3665 Table 9. Flow Paths and the Amount of the Flows Leaving Port VGO of CDU2 amount during period destination

1

2

3

4

5

TKFCC1 TKFCC2

With Penalty on the Number of Total Flows 0 3641.76 3641.76 1245.72 3641.76 0 0 2396.04

3641.76 0

TKFCC1 TKFCC2

With No Penalty on the Number of Total Flows 2571 1158.24 3641.76 1158.24 1070.76 2483.52 0 2483.52

3641.76 0

Table 10. Inventory in TKFCC1 period1

period 2

period 3

period 4

period 5

2829

With Penalty on the Number of Total Flows 4070 5312 4158

5400

5400

With No Penalty on the Number of Total Flows 4158 5400 4158

5400

amount of products produced by mode G95 during the same period is 1765.5 m3. The flows leaving port AR of unit CDU1 all move to TKFCC3. The flow paths and the amount of the flows leaving port VGO of CDU2 are shown in Table 9, and the inventory in TKFCC1 at the end of each period is shown in Table 10. It can be seen that the total number of flows is smaller if penalties are imposed on the number of flows, and thus, the frequent change of flow paths is avoided. Each actual tank in TKG90, TKG93, and TKG95 cannot receive and send materials at the same time, and they can only send materials 1 h later after they receive materials. Heuristic MAX_CAP_FRS is used for them to receive materials, and heuristic MAX_VOL_FRS is used for them to send materials. The simulation for the operations of the actual tanks in tank TKG93 is shown by the Gantt chart in Figure 4. The gray strip represents the duration of a tank receiving materials, and the white strip represents the duration of a tank sending materials. 5.2. Case 2. Case 2 is used to test the computational time for the MILP model of a larger scheduling problem. In case 2, another four operation modes are added to each of unit CDU1, CDU2, FCC1, FCC2, and FCC3. This model involves 275 discrete variables, 9549 continuous variables, and 9979 single constraints. The CPU time required for the solution of this MILP model is 153 s, and the value of the objective function is 5 014 680. The number of 0-1 discrete variables in our optimization model is determined by the number of operation modes each unit possesses, the number of flow paths needed to be controlled by constraint 13, and the length of the scheduling horizon. In a real world refinery, some modes for units can be excluded according to monthly plans and demand and supply orders during a scheduling horizon, as well as the experience of schedulers. Doing this can reduce the 0-1 variables and the size of the optimization model. In order to compare with our model, a monolithic optimization model with the same scheduling horizon and time periods is

Figure 4. Gantt chart for the operations of the actual tanks in TKG93.

built based on case 2. We refer to the study of Shah13 to model the multipurpose tanks. Each multipurpose tank within TKG90, TKG90, and TKG95 has four operation modes G90, G93, G95, and VR. Moreover, in each period of this monolithic model, the flows leaving each outflow port of unit Blen1 and Blen2 can at most move to two of these multipurpose tanks and each of P1 and P2 can receive materials from at most two of these multipurpose tanks. This model involves 630 discrete variables, 10 018 continuous variables, and 11 342 single constraints. The optimal solution of this model cannot be found within 1 h. 6. Conclusions In this paper, a hierarchy scheduling approach is introduced for the short-term scheduling of refineries. The main optimization model at the upper level and the heuristics and rules adopted in simulation system at the lower level are presented. As discussed, several reasons motivate the use of a hierarchy scheduling approach for the short-term scheduling; one of which is the difficulty in the industrial practice to represent all the constraints in a mathematical formulation when mathematical programming is adopted. The other reason for this is that mathematic representations involve a lot of continuous, binary variables and a large number of nonlinear constraints that may lead to the failure of getting feasible solutions at the right time. The hierarchical approach proposed in this paper is able to model and optimize the entire production processes in refineries. The optimization model decides the start time and end time of each operation and the quantity of materials consumed/produced by each operation mode of a unit. At the same time, some rules and heuristics are introduced to automatically make decisions on the arrangement of tanks to store materials and to adjust the material storage capacities. The iterative procedure between the upper level and the lower level is used to find the optimal solution under new material storage capacities. The simulation system at the lower level also uses heuristics and rules to determine the detailed material movements between all units with operation constraints respected. The approach proposed not only can reduce the binary variables and the size of the optimization model but can also consider the details of material movements and operation rules. Acknowledgment This research is supported by the National Natural Science Foundation of China (60421002). Appendix A: Pseudocode for the Heuristic of Adjustment of the Operation Modes of Tanks up Each INVu,t′ ) INVup u , each INIINVu′,m′,t′ ) 0, and each low INVu′,m′,t′ ) 0. Sequence (u, t) into u-queue with descending values of arg + maxt∈TP INVu,m,t . For each (u, t) in u-queue Select multipurpose tanks u′′ and corresponding aggregate tanks u′ with redundant capacities into u-set For each (u′′, u′) in u-set For each t′ > ) t - 2 up up e INVu′,t′ - INVu′,m′,t′ and INIINVu′,m′,t′ + INIIN If INVu′′ low low + INVu′′ Vu′′,m′ - ∑tet′-1 s′, u′, m′, t e INVu′,m′,t′ Continue Else Break End If

3666

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007

Add (u′′, u′) to u-candidate End For End For Select tanks with maximum capacities more than arg maxt∈TP + INVu,m,t into u-set-more from u-candidate If u-set-more is nonempty Select tank (u′′, u′) with the smallest maximum capacity in u-set-more For each t′ > t - 1 up up up INVu,t′ ) INVu,t′ + INVu′′ End For For each t′ > t - 2 up up up INVu′,t′ ) INVu′,t′ + INVu′′ INIINVu′,m′,t′ ) INIINVu′,m′,t′ + INIINVu′′ low low low INVu′,m′,t′ ) INVu′,m′,t′ + INVu′′ End For Else sequence tanks on ascending maximum capacities into u-set-less sum ) 0 For each (u′′, u′) in u-set-less + If sum > arg maxt∈TP INVu,m,t Break up Else sum ) sum + INVu′′ For each t′ > t - 2 up up up INVu′,t′ ) INVu′,t′ + INVu′′ INIINVu′,m′,t′ ) INIINVu′,m′,t′ + INIINVu′′ low low low INVu′,m′,t′ ) INVu′,m′,t′ + INVu′′ End For End If End For For each t′ > t - 1 up up ) INVu,t′ + sum INVu,t′ End For End If End For Appendix B: Pseudocode for the Calculation of Six Time Points of the Operations of Tanks For each aggregate tank u ttt ) 0 For each t Select all flows received and sent by u into flow-set For each flow in flow-set + Fs′,u′,m′,s,u,m,t ) Qs,u,m,s′,u′,m′,t/(yu,m,t + yu,m,t-1 + tu,m,t yu,m,t)/24 + FOs′,u′,m′,s,u,m,t ) Qs′,u′,m′,s,u,m,t/(yu,m,t + yu,m,t-1 + tu,m,t yu,m,t)/24 End For Sequence the time point on the ascending value of yu,m,t + + yu,m,t-1 + tu,m,t - yu,m,t Calculate the length Tsl of each slot sl according to adjacent time points For each time slot sl TFsl ) 0 and TFOsl ) 0 For each flow during this time slot TFsl ) TFsl + Fs′,u′,m′,s,u,m,t TFOsl ) TFOsl + FOs,u,m,s′,u′,m′,t End For If TFsl ) 0 and TFOsl ) 0 ttt ) ttt + Tsl End If If TFOsl > 0 and there is no active tank Select a tank u′ to receive and supply material s

according to the heuristics and rules u′-fill-start ) ttt End If If TFsl > 0 and there is no active tank Select a tank u′′ to receive and supply material s according to the heuristics and rules u′′-draw-start ) ttt End If tt ) 0 While tt < Tsl If u′ and u′′ are the same INVu′,sl ) INVu′,sl + (TFOsl - TFsl)(Tsl - tt) up If INVT > INVu′ up tt′ ) Tsl - (INVu′,sl - INVu′ )/(TFOsl - TFsl) T-full ) ttt + (tt′ - tt) ttt ) T-full T-fill-end ) ttt up INVu′,sl ) INVu′ tt ) tt′ Next u′ according to the heuristics and rules u′-fill-start ) ttt Else If INVu′,sl < 0 tt′ ) Tsl + INVu′,sl/(TFOsl - TFsl) u′-empty ) ttt + (tt′ - tt) ttt ) u′-empty u′-draw-end ) ttt INVu′,sl ) 0 tt ) tt′ Next u′ according to the heuristics and rules u′-draw-start ) ttt End If End If If there is a tank u′ or u′′ that can be selected according to the heuristics and rules and T and T′ are not the same INVu′,sl ) INVu′,sl + TFOsl(Tsl - tt) INVu′′,sl ) INVu′′,sl - TFsl(Tsl - tt) up If INVu′,sl > INVu′ and INVu′′,sl g 0 up tt′) Tsl - (INVu′,sl - INVu′ )/ TFOsl u′-full ) ttt + (tt′ - tt) ttt ) T-full u′-fill-end ) ttt up INVu′,sl ) INVu′ INVu′′,sl ) INVu′′,sl - TFsl(tt′ - tt) tt ) tt′ Next u′ according to the heuristics and rules u′-fill-start ) ttt End If up and INVu′′,sl < 0 If INVu′,sl e INVu′ tt′) Tsl + INVu′′,sl/TFsl u′′-empty ) ttt +(tt′ - tt) ttt ) u′′-empty u′′-draw-end ) ttt INVT′ ) 0 INVu′,sl ) INVu′,sl + TFOsl(tt′ - tt) tt ) tt′ Next u′′ according to the heuristics and rules u′′-draw-start ) ttt End If up and INVu′′,sl < 0 If INVu′,sl > INVu′ up tt′ ) Tsl - (INVu′,sl - INVu′ )/TFOsl tt′′ ) Tsl + INVu′′,sl/TFsl If tt′ < tt′′ u′-full ) ttt + (tt′ - tt)

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007 3667

ttt ) u′-full u′-fill-end ) ttt up INVu′,sl ) INVu′ INVu′′,sl ) INVu′′,sl - TFsl(tt′ - tt) Tt ) tt′ Next u′ according to the heuristics and rules u′-fill-start ) ttt Else If tt′′< tt′ u′′-empty ) ttt + (tt′′ - tt) ttt ) u′′-empty u′′-draw-end ) ttt INVu′′,sl ) 0 INVu′,sl ) INVu′,sl + TFOsl(tt′′ - tt) tt ) tt′′ Next u′′ according to the heuristics and rules u′′-draw-start ) ttt Then u′-full ) ttt + (tt′ - tt) u′′-empty ) ttt + (tt′ - tt) ttt ) u′-full u′-fill-end ) ttt up INVu′,sl ) INVu′ u′′-draw-end ) ttt INVu′′,sl ) 0 tt ) tt′ Next u′ and u′′ according to the heuristics and rules u′-fill-start ) ttt u′′-draw-start ) ttt End If End If End If tt ) Tsl ttt ) ttt + tt End While End For End For End For Nomenclature Indices P ) property s, s′ ) port t ) time period u, u′ ) unit m,m′ ) operation mode Sets U ) processing units, perimeter units, and aggregate tanks Mu ) operation modes of unit u TP ) time periods being scheduled ISu ) inflow ports of unit u OSu ) outflow ports of unit u NTU ) units with tanks excluded TK ) aggregate tanks TA ) aggregate tanks in which maximum storage capacities cannot be violated VTA ) aggregate tanks in which maximum storage capacities can be violated CONN ) flow paths between two units CONs,u ) flow paths that connect port s of unit u CONMs,u,m ) flow paths that connect port s of unit u in operation m PERI) perimeter units

Parameters tu,m,m′ ) time of operation m stays stable when operation mode of u is changed from m to m′ QFlow u ) lower bound for charge size of unit u QFup u ) upper bound for charge size of unit u YIEDs,u.m ) standard yield of material leaving port s of unit u in operation m INVlow u ) minimum storage capacity of aggregate tank u INVup u ) maximum storage capacity of aggregate tank u low ) minimum value of property p for flow of port s of PROs,u,m,p unit u in operation m up PROs,u,m,p ) minimum value of property p for flow of port s of unit u in operation m low DDu,m,t1,t2 ) minimum demand/supply of material transferred by unit u in operation m from period t1 to t2 up DDu,m,t1,t2 ) maximum demand/supply of material transferred by unit u in operation m from period t1 to t2 TGQFu ) target for charge size of unit u during scheduling horizon TGINVu ) target for inventory in tank u at the end of scheduling horizon CQs,u,m ) cost of flow leaving port s of unit u in operation m PLSMMu,m,m′ ) penalty for unit u changing operation from m to m′ PLSMu,m ) penalty for unit u changing operation to m PLTu ) penalty for variance of charge size of unit u PENLFs′,u′,s,u ) penalty for number of flow paths from port s′ of unit u′ to port s of unit u PLQs′,u′m′,s,u,m ) penalty for flow from port s′ of unit u′ in mode m′ to port s of unit u in mode m PLI+ u ) penalty for extra storage capacity requirement of aggregate tank u PLTGQu ) penalty for deviation from target charge size of unit u PLTGINu) penalty for deviation from target inventory of tank u CINVu ) penlty for inventory in tank u FLOWUPs′,u′,s,u ) upper bound for flow rate of flow from port s′ of unit u′ to port s of unit u Variables Qs,u,m,t ) volumetric flow of port s of unit u in operation m during period t Qs,u,m,s′,u′,m′,t ) flow from port s of unit u in operation m to port s′ of unit u′ in operation m′ during period t QFu,m,t ) charge size of unit u in operation m during period t TQFu,t ) charge size of unit u during period t yu,m,t ) binary variable denoting that unit u is in operation m during period t ∆TQFu,t ) variance of charge size of unit u between period t and t - 1 INVu,m,t ) inventory in tank u in operation m at the end of period t + INVu,m,t ) extra storage capacity requirement of aggregate tank u during period t TFs′,u′,s,u,t ) binary variable denoting that there is flow from port s′ of unit u′ to port s of unit u during period t QFRs′,u′m′,s,u,m,t ) volume fraction of flow Qs′,u′m′,s,u,m,t in the flows entering u during period t SMMu,m,m′,t ) binary variable denoting a changeover from operation mode m to m′ between time period t - 1 and t SMu,m,t ) binary variable denoting the start of operation mode m during period t

3668

Ind. Eng. Chem. Res., Vol. 46, No. 11, 2007

+ yu,m,t-1 ) time of operation mode m of unit u extended to period t from t - 1 yu,m,t ) delay of the start of operation mode m of unit u during period t tu,m,t ) time of operation mode m that can continue to stay stable during period t PROs,u,m,p,t ) value of property p for the flow of port s of unit u in operation m during period t

Literature Cited (1) Pinto, J. M.; Joly, M. Planning and scheduling models for refinery operations. Comput. Chem. Eng. 2000, 24, 2259-2276. (2) Jia, Z.; Ierapetritou, M. Efficient short-term scheduling of refinery operations based on a continuous time formulation. Comput. Chem. Eng. 2004, 28, 1001-1019. (3) Reddy, P. C. P.; Karimi, I. A. A new continuous-time formulation for scheduling crude oil operations. Chem. Eng. Sci. 2004, 59, 1325-1341. (4) Go¨the-Lundgren, M.; Lundgren, J. T. An optimization model for refinery production scheduling. Int. J. Prod. Econ. 2002, 78, 255-270. (5) Glismann, K.; Gruhn, G. Short-term scheduling and recipe optimization of blending processes. Comput. Chem. Eng. 2001, 25, 627-634. (6) Jia, Z.; Ierapetritou, M. Refinery Short-Term Scheduling Using Continuous Time Formulation: Crude-Oil Operations. Ind. Eng. Chem. Res. 2003, 42, 3085-3097.

(7) Moro, L. F. L.; Pinto, J. M. Mixed-Integer Programming Approach for Short-Term Crude Oil Scheduling. Ind. Eng. Chem. Res. 2004, 43, 8594. (8) Jia, Z.; Ierapetritou, M. Mixed-Integer Linear Programming Model for Gasoline Blending and Distribution Scheduling. Ind. Eng. Chem. Res. 2003, 42, 825-835. (9) Floudas, C. A.; Lin, X. Continuous-time versus discrete-time approaches for scheduling of chemical processes: a review. Comput. Chem. Eng. 2004, 28, 2109-2129. (10) Paolucci, M.; Sacile, R. Allocating crude oil supply to port and refinery tanks: a simulation-based decision support system. Decision Support Syst. 2002, 33, 39-54. (11) Chryssolouris, G.; Papakostas, N. Refinery short-term scheduling with tank farm, inventory and distillation management: An integrated simulation-based approach. Eur. J. Oper. Res. 2005, 166, 812-827. (12) Henning, G. P.; Jaime, C. Knowledge-based predictive and reactive scheduling in industrial Environments. Comput. Chem. Eng. 2000, 24, 2315-2338. (13) Shah, N. Mathematical programming techniques for crude oil scheduling. Comput. Chem. Eng. 1996, 20, 1227-1232.

ReceiVed for reView October 21, 2006 ReVised manuscript receiVed March 10, 2007 Accepted March 12, 2007 IE061354N