6126
Ind. Eng. Chem. Res. 2008, 47, 6126–6139
Short-Term Scheduling of Multistage Batch Plants with Unlimited Intermediate Storage Pedro M. Castro* and Augusto Q. Novais Departamento de Modelac¸a˜o e Simulac¸a˜o de Processos, Instituto Nacional de Engenharia, Tecnologia e InoVac¸a˜o, 1649-038 Lisbon, Portugal
This paper presents a new multiple time grid, continuous-time, mixed integer linear program model for the short-term scheduling of multistage, multiproduct batch plants. The formulation is based on the resource-task network process representation and explicitly considers a single intermediate storage unit per stage to keep track of material transfer between processing units belonging to consecutive stages of production. It is compared to other existing unit-specific and sequence-based approaches through the solution of a few example problems taken from the literature. Alternative objective functions such as total cost, makespan, and total earliness minimization, as well as revenue maximization, are considered. The results show that the proposed formulation performs similarly to the one of Castro and Grossmann (2005), which is limited to problems involving a single product batch per stage. In addition, it can be up to 3 or 4 orders of magnitude faster than the one of Shaik and Floudas (2008), due to a lower integrality gap, the use of a smaller number of event points, and fewer constraints. 1. Introduction In order to remain competitive in the global marketplace, companies need to rely on flexible manufacturing facilities that can handle a large product portfolio and are able to respond quickly to demand fluctuations while keeping inventories low. Multipurpose, multiproduct batch plants, where equipment resources are shared among competing tasks, have replaced dedicated, single product, continuous plants. While more flexible, multipurpose plants are also much more difficult to operate since the equipment network is more complex. In fact, a particular unit (e.g., reactor, filter) can be suitable for the execution of tasks belonging to different stages of production. Multistage plants are a special case, where equipments are allocated to a single stage. Multistage batch processes are commonly used for the production of high-value, low-volume products such as specialty chemicals and pharmaceuticals, as well as in consumer product and food processing, metal casting, and microelectronics.1 The challenging optimization problems arising from flexible plants have motivated the research community to develop efficient scheduling approaches that can have a significant economic impact. Nowadays, optimal scheduling is viewed as an essential part of enterprise-wide optimization (EWO)2,3 and despite the major developments of the last two decades there are still enormous challenges to face to render the effective solution of large-scale problems a trivial task. In terms of the modeling challenge, it can be said that state-task/resource-task network4,5 (STN/RTN) models can in theory be used for the solution of virtually all classes of scheduling problems. However, there are some critical modeling issues that are strongly dependent on the characteristics of the problem being solved. The most important concerns the treatment of time. Dividing the time horizon into a fixed number of equal length intervals is probably the best option whenever a reasonable number of time intervals can be used to account for fairly accurate processing times. On the other hand, in the presence of variable duration tasks, such as continuous tasks or batch * To whom correspondence should be addressed. Tel.: +351210924643. Fax: +351-217167016. E-mail:
[email protected].
tasks with a duration dependent on the batch size, it may be necessary to apply a continuous-time formulation or at least a mixed-time approach that allows tasks to end between fixed interval boundaries.6,7 Even with continuous-time models there are still a few possibilities for dealing with multistage problems. First, if there is no batch mixing or splitting, i.e., if each batch moves as a discrete entity throughout the process, no time grid is necessary. Such sequence-based formulations can use immediate8 or global precedence9,10 sequencing variables. Second, one can rely on a single nonuniform time grid, which has significant advantages when dealing with shared resources like storage units with limited capacity, utilities, or manpower. The drawback is that tasks may span more than a single time interval, which has motivated the development of multiple time grid (also known as unit-specific) formulations that need fewer event points to attain global optimal solutions. Sequence-based scheduling approaches are typically used following the selection, for each particular product, of the number and size of batches to produce.1 However, decoupling batching from scheduling decisions may not be an easy task in processes with parallel units of dissimilar capacities.26 Furthermore, it often leads to suboptimal solutions since one cannot take advantage of batch mixing/splitting to capitalize on the different capacities from one stage to the next. Maravelias and co-workers1,26 have recently shown that considering batching and scheduling simultaneously, yields better solutions than those obtained by existing two-stage methods. Their formulations can be viewed as sequence-based explicit batching11 approaches since there are binary variables that identify whether a particular batch of a certain order has been selected, from a prepostulated set of potential batches. The formulations were first proposed1 with global precedence 5-index sequencing variables in problems involving fixed processing times and were later26 extended to cases with batch size dependent processing times and involving sequence dependent changeovers costs, through the use of immediate precedence sequencing variables. In this paper, we present a new multiple time grid formulation that can simultaneously determine the optimal set of batches (number and size) together with the assignment and sequencing of batches in processing units of multistage plants. It can be
10.1021/ie800194b CCC: $40.75 2008 American Chemical Society Published on Web 07/24/2008
Ind. Eng. Chem. Res., Vol. 47, No. 16, 2008 6127
Figure 1. RTN process representation for multistage multiproduct plant.
classified as an implicit batching11 approach since there are no integer variables to account for the number of batches processed in each unit. It can be viewed as an extension of our previous work,12 which dealt with a subclass of the problems considered here, more specifically, those involving a single product batch per stage. The main novelty is the definition of a single intermediate storage unit per stage to which an explicit time grid is associated to. The event points belonging to the time grid of a particular storage unit are then related to the corresponding event points of time grids of processing units feeding and receiving material from the storage unit. Such timing constraints, together with the material balances, guarantee that a product can only start to be produced in stage k + 1 after enough intermediate has been produced in stage k. The new formulation is thoroughly compared to its predecessor, to the also unit-specific, RTN formulation of Shaik and Floudas,13 and to the sequence-based approach of Harjunkoski and Grossmann,9 through the solution of ten example problems taken from the literature, for alternative objective functions. The remainder of the paper is structured as follows. In section 2 the problem under consideration is defined and the general RTN process representation is given. In section 3, fundamental concepts related to unit-specific time grid formulations are given. The conceptual differences between the three scheduling models considered are identified. Furthermore, we study their impact on a vital performance indicator such as the minimum number of event points required to find the global optimal solution, through a simple example. The proposed formulation is then given in section 4, while relevant features of the others are left for section 5. Detailed computational studies are performed in section 6, where the problems are separated into single and multiple batch problems. A couple of limiting aspects of the new formulation are identified there. Finally, the conclusions are left for section 7. 2. Problem Definition In this paper, the short-term scheduling problem of multistage, multiproduct batch plants is considered. Given are a set I of products that must follow a sequence of processing stages k ∈
Figure 2. Number of event points required by Shaik and Floudas13 (nt), Castro and Grossmann12 (Tt,m), and Castro at al.15 (Tt) approaches (|I| ) 2, |MPR| ) 6, |K| ) 3).
K to reach the condition of final products. A set m ∈ M of equipment units is available at the plant, which is divided into processing units (MPR) and storage units (MST, with |MST| ) |K| - 1). Each unit is allocated to a single stage, with set Mk including the processing units belonging to stage k. A particular processing unit can handle all products belonging to set Im. Storage units are shared among the products and have unlimited capacity. Processing times are assumed to be given by a constant Ri,m plus a term βi,m that is proportional to the batch size Sizei,m, which cannot exceed the maximum batch size, Max_Sizei,m. In general, multiple batches of product i can be produced and a different number of batches may be involved in different stages of production; i.e., batches can be split and/or mixed. As a special case, we will consider problems with fixed batch sizes and a single batch per product in every stage. In such problems, the product’s release ri and due dates di are assumed to be known and are enforced as hard constraints. The process representation is given in Figure 1 in the form of a resource-task network. The plant resources (circles) are the equipment units and the material states. The processing tasks, named DO_i_m, are represented as rectangles. Products are numbered from 1 to |I|, with I(|I|) representing the last element of the set. For simplicity, it is assumed that storage units are numbered after processing units, so the latter go from M1 to M(|MPR|) while the former from M(|MPR| + 1) to M(|M|). The fundamental difference between the two types of units is that
6128 Ind. Eng. Chem. Res., Vol. 47, No. 16, 2008
Figure 3. Underlying continuous-time grid for processing units.
Figure 4. Number of event points required by the new approach (|I| ) 2, |MPR| ) 6, |K| ) 3).
unlike processing units, storage units are always available. Thus, there is no need to model the consumption and production of storage resources and this is why in Figure 1 there are no arrows starting or ending in these resources. Resources and tasks are represented within the boundaries of the corresponding stage. As an example, in stage 2 (K2), product I1 can be processed on parallel processing units ranging from M(|M1| + 1) to M(|M1| + |M2|). The corresponding tasks consume material (I1,K1) and produce (I1,K2). Note that intermediate material states are directly associated to the stage where they are produced and are located inside the storage vessel linked to that same stage, e.g. (I1,K2) only exists inside M(|MPR| + 2). Four objective functions will be considered. The first threestotal cost, revenue, total earlinessscan be thought of as being product related since the time horizon (H) is fixed. In total cost minimization, the objective will be to find the lowest allocation cost of products to units, given by ci,m, so that all products can be produced. In revenue maximization, and given the product selling prices Vi, the objective will be to produce more of the most valuable and less time-consuming products. When minimizing total earliness, the objective will be to end product production as close as possible to their due dates. Finally, makespan minimization is a time-related objective since we aim to find the minimum production time for a given amount of products (∆i). In dealing with this objective, the productivity of the plant as a whole is also maximized. 3. Fundamental Concepts The new formulation to be presented in the next section can be classified14 as a monolithic, continuous-time approach relying on multiple (also known as unit-specific) time grids for event representation. It uses the resource-task network4 (see section 2), as a unified framework for process representation, to deal with the optimal set of batches, simultaneously with the allocation, sequencing, and timing of processing tasks. Nevertheless, in comparison to the seminal paper, the general concept of resource has been disaggregated into two specific types, equipment resources and material states, in order to have a better
correspondence between the model entities and the real plant entities. As a consequence, fewer structural parameters are required for this multistage case when compared to the more general multipurpose case,14,15 making the model less abstract and easier to understand. While unified frameworks were a major breakthrough in the early 1990s,4,5 the most critical issue in scheduling models is the handling of time. Time grid based, continuous-time models are the state of the art. Among these, there are two alternatives for time representation: (i) single time grid with global event points; (ii) multiple time grids with unit-specific event points. Well-known examples of the former type can be seen in work by Castro and co-workers,15,16 Maravelias and Grossmann,17 and Sundaramoorthy and Karimi,18 while the most cited unitspecific formulations can be found in papers of Floudas and co-workers13,19–21 and Giannelos and Georgiadis.22 A detailed computational study involving all these formulations for multipurpose batch plants can be seen in Shaik et al.,23 who have employed the same hardware and software. Comparative studies for multistage, multiproduct batch plants with24 and without12 sequence-dependent changeovers and for continuous plants with intermediate storage units25 can also be found in the literature. Overall, it can be said that single time grid formulations are more general since they can always find the global optimum for a sufficiently high number of event points. On the other hand, unit-specific approaches are not as reliable and can sometimes cutoff the global optimal solution from the feasible region or fail to rigorously account for intermediate storage requirements. Historical evidence of this can be found in papers from Floudas and co-workers13,20,21 that have extended the original one.19 However, the way in which event points of different time grids are related has not changed. Nevertheless, in some circumstances, multiple time grid formulations are significantly more powerful than their single time grid counterparts due to the requirement of fewer event points per grid to find global optimal solutions. Smaller, tighter models typically result that are easier to solve. The new multiple time grid formulation presented in this paper is suitable for multistage multiproduct plants under a shortterm mode of operation. It is restricted to an unlimited intermediate storage policy since even though there are variables that hold the amount of materials in storage at the unit’s event points, they are generally not accurate. More specifically, under a finite intermediate storage policy, it cannot be guaranteed that predefined storage capacity is not exceeded and that the resulting solutions are not in fact infeasible. 3.1. Relation of Event Points Belonging to Different Time Grids. In this section, we discuss the conceptual differences between the unit-specific formulation of Floudas and coworkers, based on their newly improved version,13 the multiple time grid formulation of Castro and Grossmann,12 and the proposed new approach. They all rely on the resource-task network process representation. Besides the type of network structure that can be handled, there is a clear distinction between the two published models, which can be realized following an analysis of the timing variables. While the model by Shaik and Floudas13 is called a unit-specific approach, the model does not rely on explicit time grids for equipment units. In fact, the model timing variables Ts(i,n) link event points n with tasks i. The excess resource balances then ensure that two tasks cannot start at the same event point in the same unit, while the timing (sequencing) constraints guarantee that the event points are properly spaced. The unit’s time grids are naturally implicit in these constraints
Ind. Eng. Chem. Res., Vol. 47, No. 16, 2008 6129
Figure 5. RTN for example P12 (duration in hours, size in tonnes).
and can be generated following the solution of the model. In general, a unit-generated time grid will feature only some of the event points. In view of the above, Floudas and co-workers models are better classified as task-specific event point based rather than unit-specific. On the other hand, the model by Castro and Grossmann12 is a true unit-specific approach, since one time grid is explicitly defined for each processing unit, with the timing variables Tt,m holding the time of event point t of the time grid linked to unit m. The timing variables then relate the times of two consecutive event points belonging to the same grid and ensure that a product can only start to be processed in stage k + 1 after going through stage k. What is interesting is that these conceptual differences influence the number of event points that is required to find the global optimal solution as can be seen in the simple example given in Figure 2, involving 2 products, 6 processing units, and 3 stages. In the example, every unit is processing exactly one task. Using the model by Shaik and Floudas,13 we can solve the problem by using 3 event points (n1-n3). Note that due to their sequencing constraints relating different tasks in different units and the excess resource balances, production of product i in stage k starts necessarily at least one event point after production in stage k - 1. Thus, I1 and I2 start to be produced at events points n1 in stage 1 (units M1 and M2), n2 in stage 2 (M3, M4), and n3 in stage 3 (M5, M6). If one uses the model by Castro and Grossmann,12 only two event points (Tt,m) are required per unit. It is particularly important at this point to emphasize that the last event point is not important since the relevant event points are only those where tasks are started. This model can be easily reformulated to use one less event point, so from now on it is the reformulated version that will be considered for comparative purposes. Its underlying time grid, which is also used by the new proposed model, is given in Figure 3. It should be noted that the location of event points may vary from unit to unit.
As can be seen in Figure 2, the first event points of units M1 and M2, namely T1,1 and T1,2, are placed at the origin (T1). The first event points of units M3 and M4 are located after the end of the task whose output is the product that they will be processing, i.e. T1,3 g T1,2 + RI2,2 and T1,4 g T1,1 + RI1,1 (assuming that the duration of the task is independent of the batch size). Note that all tasks can potentially start at the first event point, regardless of the stage where the unit in which they are processed is located. In Shaik and Floudas13 approach, neither second stage tasks can start at the first event point, nor can first stage tasks start at the last event point (in multistage problems). The corresponding binary variables can therefore be eliminated. Finally, Figure 2 also shows the event point location for the single continuous-time formulation of Castro at al.15 and it can be concluded that in such a case, a significantly higher number of event points (7) is required. For the objectives of total cost and makespan minimization, the new formulation can lead to a degenerate solution to that shown in Figure 2 with just 2 event points per unit (see Figure 4), while no solution can be found using just 1 event point. Thus, this simple example also serves to illustrate one important property of the reformulated formulation of Castro and Grossmann.12 It is a minimum event points formulation and it is not conceptually possible for a time grid based continuous-time formulation to solve a multistage, multiproduct problem, to global optimality by using fewer event points. While in Castro and Grossmann12 explicit timing variables TDi,k are used to ensure that product i in stage k + 1 starts to be processed after stage k has ended, in the new formulation this is done implicitly by means of the time grid of the intermediate storage units, M7 and M8, for stages 1 and 2, respectively. More specifically, the model constraints for M7 imply T1,7 g T1,2 + RI2,2 and T1,7 g T1,1 + RI1,1. Since M7 feeds units M3 and M4, then T1,3 g T1,7 and T1,4 g T1,7. From these constraints it should have been possible for I1 and I2 to start at the first event point in M1 and M2. However, since I1 is the bottleneck, M6 cannot wait for
6130 Ind. Eng. Chem. Res., Vol. 47, No. 16, 2008
I2 to end its production in M3 so that the tasks could also start at the first event point in units M5 and M6. As a consequence, one more event point is required, with I1 starting at the first event point in units M1, M4, and M6 and I2 starting at the second event point in M2, M3, and M5. The relevant constraints for M8 are all valid since T1,8 g T1,3 (note that I2 does not start at event point 1 of M3 so there is no need to consider the RI2,3 term on the RHS), T1,8 g T1,4 + RI1,4, T1,6 g T1,8, T1,5 g T1,8, T2,8 g T2,4, T2,8 g T2,3 + RI2,3, T2,5 g T2,8, and T2,6 g T2,8. This schedule requires a total of 5 event points (T1 - T5) when using a single time grid.15 The advantage of the new model over that given in Castro and Grossmann12 is that when relating the time grid of the storage unit with those of the preceding and succeeding processing units, we do not need to know which orders are being executed and no big-M constraints are required, which generally have a detrimental effect on the integrality gap and computational effort. But more importantly, the new model is no longer limited to the execution of a single batch per product per stage, and so it is considerably more general. 4. New Multiple Time Grid Continuous-Time Formulation (CN) The new continuous-time formulation has many similarities to the one by Castro and Grossmann.12 First, some model variables are the same. Binary variables Ni,m,t identify the execution of product i in process unit m starting at event point t. Variables Rm,t identify the excess amount of process unit m at time t. Although they can be defined as continuous variables, for performance reasons it is better to define them as binary variables. Continuous variables Tt,m give the time of event point t belonging to the time grid of unit m. Now, however, these timing variables refer not only to process units but also to storage units. By considering the intermediate storage units explicitly, it is possible to keep track of material inventory over time. Excess resource variables Si,m,t (m ∈ MST) are used for that purpose and give the amount of product i in intermediate storage unit m at event point t. Note that when compared to both the discrete and single time grid continuous-time formulations, also given in Castro and Grossmann,12 there has been an index change from Si,k,t. The other additions are the continuous extent variables ξi,m,t that account for the amount of material i produced in unit m at event point t. With these variables it becomes possible to model features like variable batch sizes and variable duration tasks, which are particularly relevant in problems involving multiple product batches. The model constraints are given next. Starting with the excess resource balance constraints, eq 1 relates the excess amount of unit m at time t to the occurrence of processing tasks starting at that same event point. If a task is executed, the unit becomes unavailable (Rm,t ) 0) for other tasks, thus ensuring that a particular equipment unit can only handle one product at a time. A fundamental property of the model is that it can assume, without loss of generality, that all tasks can last a single time interval. Thus, tasks consuming the equipment unit at t give it back to the system at t + 1. This is the reason why eq 1 looks slightly different from a typical multiperiod balance constraint, which features a term with the value of the variable at an earlier time (see eq 2). In eq 2, the production of product i in processing unit m′ belonging to the same stage of storage unit m (m′ ∈ Mmto) increases the availability of the material by the amount processed in m′, while production in processing unit m′′ belonging to the next stage (m′′ ∈ Mmfr), decreases inventory levels by the amount consumed in m′′. It is assumed that the
amount processed by the task is equal to the amount of material consumed, while the yield leading to the amount produced can be accounted for through parameter νji,m. Rm,t ) 1 -
∑N
∀ m ∈ MPR, t ∈ T
i,m,t
(1)
i∈Im
Si,m,t ) Si,m,t-1 +
∑
νi,m′ξi,m′,t -
to m′∈Mm
i∈Im′
∑
∀ i ∈ I, m ∈ MST, t ∈ T (2)
ξi,m′′,t
fr m′′∈Mm
i∈Im′′
The timing constraints relate the absolute times of consecutive event points in the same time grid and those of a particular event point in different time grids. Equation 3 forces the duration of time interval t of processing unit m to be greater that the duration of the task taking place. Equations 4 and 5 are new and together with eq 2 ensure that a product can only start to be produced in stage k + 1 after production in stage k has finished. Equation 4 states that if processing unit m sends material to storage unit m′, then the absolute time of event point t belonging to grid m′ must be greater than that of t in grid m plus the duration of the task taking place. Equation 5 then ensures that if processing unit m receives material from storage unit m′, the absolute time of event point t belonging to grid m′ must be lower than that of t in grid m. The timing variables are naturally bounded by the time horizon, H (see eq 6). H|t)|T| + Tt+1,m|t*|T| - Tt,m g
∑ (R
i,mNi,m,t +
i∈Im
∀ m ∈ MPR, t ∈ T
(3)
to ∀ m ′ ∈ MST, m ∈ MPR, m ∈ Mm′ , t∈T
(4)
βi,mξi,m,t) Tt,m′ g Tt,m +
∑ (R
i,mNi,m,t +
i∈Im
βi,mξi,m,t) Tt,m′ e Tt,m
∀m′ ∈M , m∈M , ST
Tt,m e H
PR
fr m ∈ Mm′ ,
t∈T
∀ m ∈ M, t ∈ T
(5) (6)
Finally, eq 7 and a few more objective function related constraints that are described in the following subsections complete the formulation. It is a capacity constraint that guarantees that the amount of material involved in the task does not exceed the maximum amount that can be handled by the processing unit (Vmmax). ξi,m,t e Vmmax Ni,m,t
∀ m ∈ MPR, i ∈ Im, t ∈ T
(7)
4.1. Objective Function Related Constraints. 4.1.1. Revenue Maximization. In problems involving several production alternatives and possibly multiple product batches over a fixed time horizon, the most suitable objective might be revenue maximization, eq 8. Note that we are accounting for the total amount produced in processing units belonging to the final production stage. max
∑ ∑ ∑ V ·ν i
i,m·ξi,m,t
(8)
m∈M|K| i∈Im t∈T
4.1.2. Makespan Minimization. In makespan minimization, the objective is to minimize the production time of a certain amount of products. MS being the variable representing the makespan, eq 9 replaces eq 3 and ensures that all tasks end before the makespan. Although not required, the formulation becomes tighter12 if eq 10 is used. When compared to the
Ind. Eng. Chem. Res., Vol. 47, No. 16, 2008 6131 12
original constraint, there is a second term inside the square brackets to account for the fact that the duration of the tasks is dependent on the batch size. Note that we cannot simply add the minimum size-dependent parameters of processing units belonging to subsequent stages (βi,m′) because it may be possible to divide the total amount among such stages processing units. Furthermore, we need to account for the units yields. The lower bounds on product production are achieved through eq 11.
∑ (R
MS|t)|T| + Tt+1,m|t*|T| - Tt,m g
i,mNi,m,t +
i∈Im
Tt,m e MS -
∑
i∈Im
(
βi,mξi,m,t)
[ (
Ni,m,t Ri,m +
ξi,m,t βi,m + νi,m
[(
∑
k′>k
)
k′′∈K
] )]
k 0
(34)
s + Ri,mNi,m,t e diNi,m,t + H(1 Ti,m,t Ni,m,t) ∀ m ∈ M|K|, t ∈ T, i ∈ Im,t (35) s + Ri,mNi,m,t + H(1 DDi e Ti,m,t Ni,m,t) ∀ m ∈ M|K|, t ∈ T, i ∈ Im,t (36) s + Ri,mNi,m,t - H(1 DDi g Ti,m,t Ni,m,t) ∀ m ∈ M|K|, t ∈ T, i ∈ Im,t (37)
min
∑ (d - DD ) i
i
(38)
i∈I
5.2. Castro and Grossmann12 Reformulated (CG). The multiple time grid formulation of Castro and Grossmann12 is not as general as the other two since, instead of keeping track of the inventory levels of intermediate materials, it assumes that there is a single transfer of material per product from one stage to the next. For this reason, it can only address problems involving a single batch per product on each stage. Variables TDi,k hold the transfer time for product i in stage k, which must be greater than the product finishing time in that stage (eq 39) and lower than its starting time in stage k + 1 (eq 40). Equations 41 and 42 improve the performance both for total cost minimization and makespan minimization, while eqs 43–45 replace them for total earliness minimization. For makespan minimization, it is also good to relate such variables with the makespan variable; see eq 46. Finally, eq 47 reduces the number of degenerate solutions by enforcing the allocation of tasks to event points with as low an index as possible.24 Note that this will occur even without that constraint for the objective of total earliness minimization (see section 4.1.4). The reformulated version of the model, which requires one less event point than the original one12 as discussed in section 3.1, is completed with eqs 1, 3 or 9 and 10 (without the terms featuring parameter βi,m), 6, 15–20, and either eq 12, 13, or min MS as the objective function.
Ind. Eng. Chem. Res., Vol. 47, No. 16, 2008 6133
TDi,k g Tt,m + Ri,mNi,m,t -
(
H 1-
∑N
i,m,t′
t′∈T t′gt
TDi,k-1 e Tt,m +
(
H 1-
∑N
)
∀ k ∈ K, k * |K|, m ∈ Mk, i ∈ Im, t ∈ T
t′et
TDi,k g ri +
(39)
)
∀ k ∈ K, k * 1, m ∈ Mk, i ∈ Im, t ∈ T (40)
∑∑ ∑
Ri,mNi,m,t
∀ i ∈ I, k ∈ K, k * |K|
t∈T k′∈K m∈Mk′ k′ek i∈Im
(41) TDi,k e di -
∑∑ ∑
Ri,mNi,m,t
∀ i ∈ I, k ∈ K, k * |K|
t∈T k′∈K m∈Mk′ k′>k i∈Im
(42) TDi,k g ri +
∑∑
∀ i ∈ I, k ∈ K, k ) 1
Ri,mNi,m,t
t∈T m∈Mk i∈Im
(43) TDi,k e di -
problem/model
optimum CG
P3C (|I| ) 8, |M | ) 6,|K| ) 2) P4C (|I| ) 8, |MPR| ) 6,|K| ) 2) P5C (|I| ) 10, |MPR| ) 6,|K| ) 3) P6C (|I| ) 10, |MPR| ) 6,|K| ) 3) P7C (|I| ) 12, |MPR| ) 8,|K| ) 3) P8C (|I| ) 12, |MPR| ) 8,|K| ) 3) P9C (|I| ) 15, |MPR| ) 6,|K| ) 2) P10C (|I| ) 10, |MPR| ) 8,|K| ) 4) PR
i,m,t′
t′∈T
Table 1. Single Batch Problems: Overview of Computational Performance (CPU s) for Total Cost Minimizationa
∑ ∑
Ri,mNi,m,t
∀ i ∈ I, k ∈ K, k ) |K|-1
t∈T m∈Mk+1
16 1063 89 886 61 664 88 154
0.12 0.16 6.26 3.6 54.9 189 7.91 1994
CN
SF
HG
0.14 0.3 1.16 20.5 18 3.16 177 30730b
0.22 0.14 1203 839 60000c 60000d 60000e 25369
0.12 0.11 0.16 0.11 3.86 2.00 15.7 363
a BPS ) best possible solution at the time of termination. MRL)maximum resource limit exceeded. NS ) no solution found. OM ) solver ran out of memory. SO ) suboptimal solution returned. b SO ) 155. c MRL, BPS ) 60. d MRL, NS, BPS ) 664. e MRL, SO ) 89, BPS ) 87.
Table 2. Single Batch Problems: Overview of Computational Performance (CPU s) for Makespan Minimization problem/model P3M P4M P5M P6M P7M P8M P9M P10M
optimum
CG
CN
SF
HG
793 793 1435 1394 1456 1655 235 252
0.62 0.39 6.11 8.3 521 79.6 60000a 56000b
153 5.31 31.7 0.17 4.11 10.6 60000c 5.34
634 355 5759 1801 60000d 60000e 60000f 60000g
0.41 0.86 1615 7.01 7200h 3557 26084i 9000j
a MRL, SO ) 236, BPS ) 232. b MRL, NS, BPS ) 252. c MRL, BPS ) 233.46. d MRL, NS, BPS ) 1200. e MRL, NS, BPS ) 1320. f MRL, SO ) 246, BPS ) 188.42. g MRL, BPS ) 178.99. h MRL, BPS ) 1455. i OM, SO ) 243, BPS ) 225.21. j MRL, SO ) 254, BPS ) 243.
i∈Im
(44) TDi,k g TDi,k-1 -
∑∑
Ri,mNi,m,t
problem/model
∀ i ∈ I, k ∈ K, k * 1, k * |K| (45)
t∈T m∈Mk i∈Im
TDi,k e MS -
Table 3. Single Batch Problems: Overview of Computational Performance (CPU s) for Total Earliness Minimization
∑∑ ∑
Ri,mNi,m,t
∀ i ∈ I, k ∈ K, k * |K|
t∈T k′∈K m∈Mk′ k′>k i∈Im
(46) Rm,t g Rm,t-1
∀ m ∈ MPR, t ∈ T, t * 1
(47)
6. Computational Results The performance of the three mathematical formulations is illustrated through the solution of 10 example problems (P3 to P12). These have been divided in two sets. The first set (P3 to P10) deal with problems involving a single product batch per stage and includes 8 examples taken from Castro and Grossmann12 (P3-P8 are originally from Harjunkoski and Grossmann9). Since in ref 12 three models other than those which are now considered have solved the problems, the same nomenclature will be used here to facilitate comparison. A character after the number identifies the objective function being considered, e.g., C for total cost minimization, E for total earliness, and M for makespan. In the second set of problems, a few product batches are involved so formulation CG can no longer be used. Problem P11 corresponds to example 1 of Shaik and Floudas13 and deals with the production of a single product in a multistage plant (examples 2 and 3 involve multipurpose plants that cannot be solved by CG neither CN). In P12 we have kept the same product data for I1 (νI1 ) 5 $/t) but added
P3E P4E P5E P6E P7E P8E P9E P10E
optimum
CG
CN
SF
HG
135 135 35 69 700 1380 228 184
2.01 0.73 8.08 65.6 295 20.6 8.31 45520
3.55 2.78 14.8 3.39 314 8.58 20.6 5.84a
57 359 3574 610 60000b 60000c 60000d 60000e
9.98 17.2 0.38 0.94 0.39 1.97 10835f 60000g
a SO ) 192 that does not improve with an increase in |T| (special case). b MRL, BPS ) 0. c MRL, NS, BPS ) 0. d MRL, SO ) 417, BPS ) 0. e MRL, NS, BPS ) 0. f OM, SO ) 239, BPS ) 17.46. g MRL, BPS ) 173.22.
another product, I2 (νI2 ) 6 $/t). The remaining data can be found in the RTN given in Figure 5. The MILP models were implemented and solved in GAMS 22.5, using CPLEX 10.2 as the solver. All problems were solved to optimality (1E-6 relative tolerance) unless otherwise stated. The computer used was a Pentium-4 3.4 GHz processor with 2 GB of RAM running Windows XP Professional. 6.1. Single Batch Problems. The single batch problems under consideration range from 8 products in 6 processing units in 2 stages (P3-P4) to 15 products in 6 units in 2 stages (P9) and 10 products in 8 units in 4 stages (P10). In order to achieve more general results, the sequence-based model of Harjunkoski and Grossmann9 (HG), as described in Castro and Grossmann,12 is considered together with the three unit-specific formulations presented in sections 4 and 5. The value of the optimal solution and the total computational efforts are given in Tables 1–3, one for each objective function. 6.1.1. Overview of Computational Performance. The model with global precedence sequencing variables of Harjun-
6134 Ind. Eng. Chem. Res., Vol. 47, No. 16, 2008 Table 4. Number of Event Points (|T|) Employed By the Different Models objective cost minimization
makespan minimization
earliness minimization
problem/model
CG
CN
SF
CG
CN
SF
CG
CN
SF
P3 P4 P5 P6 P7 P8 P9 P10
3 3 6 6 12 12 6 5
3 3 6 7 12 12 7 8
4 4 8 9 14 14 7 8
3 3 5 5 12 12 6 5
4 4 6 5 12 12 8 6
4 4 7 7 14 14 7 8
3 3 5 5 12 12 6 5
3 3 5 5 12 12 6 6
4 4 7 7 14 14 7 8
koski and Grossmann9 (HG) emerges as the best overall performer since it is the best in half of the examples. Furthermore, it needs to be solved just once to find the global optimal solution, which is a major advantage over time grid based models. In other words, since it does not rely on explicit time grids, no iterative search procedure over |T| is required. In terms of objective functions, HG is particularly powerful for total cost minimization, being outdone by the model of Castro and Grossmann12 (CG) only for P9C. Model CG is actually the best performer among the multiple time grid approaches even though it failed to find the global optimal solutions for problems P9M and P10M. For the latter, no feasible solution could be found even after leaving the computer running overnight. Such difficulties in solving the fourstage problem had already been reported for the original model12 and this is probably due to the larger number of big-M constraints involved (eqs 39 and 40). The new formulation (CN) is very competitive and, like HG, could always find feasible solutions to the problems. It is actually the best performer for makespan minimization. When comparing these results to those of our previous work,12 namely those from a constraint programming model, which is generally regarded as the best approach for makespan minimization, one concludes that the performances are identical. In particular, P10 M was solved by CN in just 5.34 s, while OPL Studio 3.7.1 took 16 679 s on a Pentium-4, 2.8 GHz PC. Failures in finding the global optimal solution happened in P10C, due to the use of an insufficient number of event points, and P10E, which is a special case to be discussed in section 6.1.5. The model by Shaik and Floudas13 comes in last place. While it could solve examples P3-P6 in an acceptable amount of time, problems P7-P8 were found to be mostly intractable since in four out of six instances no feasible solution could be attained. Nevertheless, SF did found the optimal solutions without proving optimality for P7C and P7E. Such difficulties for P7-P8 arise mainly because those problems feature a few forbidden paths that lead to the need of higher number of event points that would normally be expected. In problems P9-P10, there was one infeasible, three suboptimal, and two optimal solutions being returned, one of which without proving optimality. Nevertheless, SF performed better than CN in P10C and better than CG in P10M. 6.1.2. Event Points Used. The number of event points required to find the global optimal solution is an important performance indicator in time grid based continuous-time formulations. It directly influences the size of the mathematical problem since the large majority of variables and constraints feature index t, as can be seen in sections 4 and 5. Hence, it is reasonable to aim at models that use fewer event points and this has actually been the main motivation to develop unitspecific models as opposed to single time grid models. But is
this a good heuristic? Table 4 will help us to answer this question. The most important result from Table 4 is that Castro and Grosmann’s12 formulation is the one that requires the smallest number of event points. With the exception of P10M, we know that the |T| values listed in the CG columns are enough to find the global optimal solutions and that in 15 out of the 24 examples, such values are also sufficient for the new formulation (CN). In the remaining cases, one extra event point is usually sufficient. The two exceptions occurred with P9M, two event points above the minimum, and P10C. For the latter, the computational statistics (see Table 1) are given for 8 event points despite leading to a suboptimal solution (155), simply because it already takes a considerable amount of time to solve. In fact, to generate the global optimal solution obtained with CG, a total of 10 event points are actually required by the new formulation. The same problem also serves to illustrate that the widely used procedure for time grid based formulations of stopping the search for the global optimal solution when the solution is the same for two consecutive values of |T| is sometimes inappropriate. Note that when solving P10C with CN, we found a solution worth 156 for both 5 and 6 event points, and 155 for |T| ) 7, 8 and probably 9. The comparison to Shaik and Floudas13 formulation (SF) is more difficult to perform since the global optimal solutions could only be found in two-thirds of the problems. For the others, the values listed in Table 4 are one event point above those for which the problem becomes infeasible, which can be proved rather rapidly. Interestingly, we do not get a constant difference between the number of event points for CG and SF, neither between those for CN and SF. From the discussion in sections 3.1 and 5.1, we know that SF implicitly assumes that a product being produced in stage k + 1 cannot start before event point t + 1, if started at event point t in stage k. This artificial increase in the number of event points with the number of stages is naturally reflected in the results but the difference in the number of event points is not always equal to |K| - 1. In fact, the new CN formulation needs more event points than SF for P10C. This is evidence that the three models are indeed conceptually different. Nevertheless, for all instances solved of P11 and P12, the following condition applies: |T|SF ) |T|CN + |K| - 1. To help us answer the question posed at the beginning of this section, we identify the formulation using the smallest number of event points in bold. In 9 problem instances, CG used fewer event points than the others and when comparing the results given in Tables 1–3, it was indeed the best performer in 7 of them (CN did find better solutions for P9M and P10M). Thus, it can be concluded that the heuristic that ranks models based on the number of event points needed to find the global optimal solution, is a good one.
Ind. Eng. Chem. Res., Vol. 47, No. 16, 2008 6135 Table 5. Integrality Gap (%) objective cost minimization
makespan minimization
earliness minimization
problem/model
CG
CN
SF
HG
CG
CN
SF
HG
CG
CN
SF
HG
P3 P4 P5 P6 P7 P8 P9 P10
0.0 0.0 0.0 1.4 1.6 0.0 1.0 4.8
0.0 0.0 0.0 1.6 1.6 0.0 1.5 7.1
0.0 0.0 0.0 1.6 1.6 0.0 1.1 5.8
0.0 0.0 0.0 0.6 1.6 0.0 1.4 6.6
13.3 13.2 0.6 2.8 0.0 2.4 5.4 1.5
19.8 19.6 10.9 2.6 0.0 2.4 16.4 11.4
24.6 24.6 19.9 23.0 17.6 20.2 20.6 31.5
7.6 7.6 7.7 11.3 0.1 6.6 6.9 4.8
100 100 100 100 100 81.2 71.8 16.2
100 100 100 100 100 81.2 71.8 35.3
100 100 100 100 100 100 100 100
100 100 100 100 100 100 100 100
Table 6. Number of Model Entities Used by the Different Formulations entity discrete variables problem/model P3C P4C P5C P6C P7C P8C P9C P10C
continuous variables
constraints
CG
CN
SF
HG
CG
CN
SF
HG
CG
CN
SF
HG
159 159 384 384 1224 1224 576 440
159 159 384 448 1224 1224 672 704
165 165 396 460 1240 1240 582 464
103 103 193 193 292 292 300 260
186 186 441 441 1345 1345 628 511
205 205 553 645 1633 1633 827 1033
339 339 905 1047 2705 2705 1228 1105
120 120 224 224 329 329 331 301
273 273 723 723 2019 2019 775 893
146 146 353 405 875 875 346 633
1135 1177 3654 4442 15104 15248 8329 4373
384 384 601 601 1163 1163 1342 839
6.1.3. Integrality Gap. An important performance indicator for MILPs is the integrality gap (IG), which measures the difference between the optimal solution of the relaxed problem (linear) and the optimal solution of the MILP. In principle, the lower the gap, the fewer the number of branches in the search tree that will be searched before finding the optimal solution and proving optimality, and the lower the computational time. Table 5 lists the relative values, with respect to the optimal solution of the MILP for the 96 instances solved. It should be noted that for the Shaik and Floudas formulation13 there is no guarantee, for examples P7M, P8-P9, and P10E, that the value of |T| used is sufficient to find the optimum. Since increasing the number of event points typically increases the integrality gap, such values can be viewed as lower bounds. The most important result is that only in 20% of the cases the best performer did not feature the lowest integrality gap. Model CG is normally the tightest formulation and with the exception of P6M, P10C, and eventually P9C, we get the following rank among the unit-specific approaches: IGCG e IGCN e IGSF. The integrality gap usually increases when going from total cost to makespan and then to total earliness minimization. However, this does not necessarily mean that the computational performance decreases in that order as illustrated for CN with P9C, P9M, and P9E, which take 177, 60 000, and 20.6 s to solve (see Tables 1–3) for integrality gaps of 1.5, 16.4, and 71.8%, respectively. And the reason is that the number of event points used (7, 8, 6) also plays an important part. Finally, it is important to highlight that the new formulation and CG can lead to integrality gaps lower than 100% for total earliness minimization contrary to SF, HG, and also the single time grid formulation of Castro et al.15 6.1.4. Size of the MILPs. The four formulations use different sets of variables and constraints, so it is appropriate to compare the size of the mathematical problems. Table 6 lists the corresponding values for the objective of total cost minimization, where it can be seen that HG uses the smallest number of binary and continuous variables while the new formulation, CN, uses the fewest constraints. Note that HG uses binary assignment, yi,m, and sequencing variables, xi,i′,k (i′ > i), while the others rely on binaries Ni,m,t and Rm,t. As for the continuous variables,
HG only needs to account for the completion time of product i in stage k (ci,k), while the unit-specific models use Tt,m (CN, s CG), Ti,m,t (SF), Si,m,t (CN, SF), and TDi,k (CG). The number of constraints in HG exceeds the ones for CN and sometimes for CG, since 4-index constraints (i,i′,m,k), which are of the big-M type, are used to ensure that different products are not being produced on a particular unit at the same time. Among the multiple time grid models, CG and CN use the exact same number of binary variables whenever the same number of event points is used (refer also to Table 4), while the number for SF is slightly higher. Nevertheless, it is relevant to emphasize that, if variables Rm,t are defined as continuous, the same number of binary variables also results for SF even when using a different value of |T| (recall that the binary extent variables Ni,m,t have a different domain in SF; see eq 22). The differences start to appear when comparing the number of single variables and are clearly obvious in terms of number of constraints. Model CG uses the fewest variables essentially because it relies on variables TDi,k instead of Si,m,t to ensure a valid transfer of material between consecutive stages. Models CN and SF use more or less the same number of variables Si,m,t so the discrepancy between them is due to the fact that the timing variables in the latter are product dependent. Taking P7C s as an example, SF uses 1128 Ti,m,t variables, while CN uses 120 Tt,m variables. The proposed model requires the smallest number of constraints since even though it requires the material balances (eq 2) and timing constraints for the intermediate storage units (eqs 4 and 5), these have fewer indices in their domain than the transfer time constraints in CG (eqs 39 and 40). For problem P7C, such constraints total 405 for CN and 1416 for CG, with this difference being the large majority of that listed in Table 6. The number of constraints in SF is 1 order of magnitude larger than for CN mostly due to the need to relate the starting times of different tasks taking place on the same unit (eq 26 is disaggregated into 11 132 constraints for P7C). This is probably why SF has a high failure rate in P7-P8. 6.1.5. Limitation of New Formulation for Total Earliness Minimization. We have seen that the new formulation returns a suboptimal solution (192) for problem P10E. Optimality of such solution is attained rather rapidly for 6 event
6136 Ind. Eng. Chem. Res., Vol. 47, No. 16, 2008
Figure 6. Best solution obtained by new formulation for problem P10E (not the global optimum).
points, so it is possible to continue the search procedure. Interestingly, the same value is returned up to |T| ) 10, which is below the global optimum of 184. This is a clear indication that the new formulation cannot find the global solution, which requires just 5 event points by CG. In other words, every processing unit handles exactly 5 products. If, in model CN, one forces the optimal assignments of products to units without being concerned with the starting event point of the tasks, the resulting mathematical problem is infeasible up to 8 event points for which the solver finds a solution of 854. This is a rather poor value and it is not the real total earliness, which equals 211. This limitation of the new formulation, stemming from the way in which the objective of total earliness minimization is defined, had already been explained in section 4.1.4. What happens is that for the optimal assignments, the new formulation cannot allocate tasks to the final stage units from the first to the last time interval. In fact, no tasks start at event points 2, 3, and 5 for M7 or at event points 1, 4, and 7 for M8. As a consequence, the starting times of those intervals cannot be made equal to the time horizon (e.g., T1,8 ) 124 < 288 ) H) and the value from eq 13 becomes much larger than the total earliness. The best solution returned from the new formulation for P10E and 6 event points is given in Figure 6. The location of all event points in the time grids of the processing units is given. Notice that there are two units (M2 and M7) handling 6 products and that there are consecutive event points with the same absolute time, T5,6 ) T6,6 ) 218 and T5,8 ) T6,8 ) 288. The latter correspond to starting points of unused time intervals and since the time value is the horizon, in this case, the value of the objective function is the real total earliness. 6.2. Multiple Batch Problems. In this paper and in ref 12, the single batch problems have been solved by 7 alternative formulations so it is fair to say that the global optimal solutions are already known. On the other hand, multiple batch problems are significantly more complex and were the main motivation for the development of the proposed model. Problems involving a single (P11) or two products (P12) are considered for two objective functions. A few instances involving different time horizons (H) for revenue maximization, or different production goals (∆i) for makespan minimization, are tested. The results are given in Tables 7 and 8 and are grouped for |T|SF ) |T|CN + 2 since we have found out that this is the relation between the number of event points required by the two formulations to find a particular solution, as already highlighted in section 6.1.2.
For revenue maximization, the results in Table 7 for P11 show virtually the same performance for SF and CN. This is no surprise considering that the problem sizes and the integrality gap are the same. However, the two models react differently to the addition of a second product. The difference between the number of continuous variables and constraints increases and the major change occurs in the integrality gap, which is significantly lower for CN. Because of this, a better performance is achieved for the new formulation, which becomes more evident as the time horizon increases. Notice that for P12 (H ) 16), CN is 1 order of magnitude faster than SF, for the minimum number of event points that leads to the global optimal solution (e.g., 6 and 8, respectively), and 50 times as fast when |T| is increased by 1. When the objective function is changed to makespan minimization, we get the confirmation that the new formulation is tighter, even for the single product problem (P11) where we get a very small difference in integrality gap; see Table 8. Again, the overall computational performance for P11 is similar between CN and SF for |T|CN + 2, even though there is a higher degree of variability, most noticeable when ∆i ) 2000/|T| ) 13 for CN and ∆i ) 2000/|T| ) 15 for SF, for which the former is 5 times faster. For P12, one confirms that the new formulation is a better one since CN can find the global optimal solution, 29.748 h, in just 16.2 s as opposed to 3.5 h by SF (almost 3 orders of magnitude slower). Furthermore, the solution is confirmed in less than 4 min for a single increase in |T|, while the same is not true for SF, even after 16 h. 6.2.1. Intermediate Storage Levels. Variables Si,m,t hold the amount of product i in intermediate storage unit m at event point t but can these values be used to generate the intermediate storage profiles? To help us answer this question, Figure 7 gives the optimal schedule for problem P12 when minimizing makespan with the new formulation. Starting with product I2 (in red), which in the first stage is processed mostly in M1, the first 100 tonnes are produced at the 3 h mark. Also starting at t ) 1 but in unit M2, the first batch of product I1 in stage 1 ends at 3.333 h, which is the time of the first event point of storage unit M6 (see Figure 8). At the same time, 150 tonnes of I1 start to be processed in stage 2 (unit M3), so no storage is required for I1 at state K1 (see Figure 5). The values of the excess resource variables Si,m,t for unit M6 (amount in storage) are shown in Figure 8 as small line segments. It can be seen that SI1,M6,1 ) 0 and SI2,M6,1 ) 100. The second event point of
Ind. Eng. Chem. Res., Vol. 47, No. 16, 2008 6137 Table 7. Computational Results for Problems P11-P12 for Revenue Maximization problem
model
|T|
discrete variables
single variables
constraints
RMIP ($)
MIP ($)
CPU s
nodes
P11 (H ) 8)
CN SF CN SF CN SF CN SF CN SF CN SF CN SF CN SF CN SF CN SF CN SF CN SF
2 4 3 5 4 6 5 7 7 9 8 10 3 5 4 6 4 6 5 7 6 8 7 9
20 30 30 40 40 50 50 60 70 80 80 90 45 55 60 70 60 70 75 85 90 100 105 115
49 59 73 81 97 103 121 125 169 169 193 191 109 136 145 175 145 175 181 214 217 253 253 292
47 59 70 80 93 101 116 122 162 164 185 185 91 149 121 196 121 196 151 243 181 290 211 337
2000 2000 2493.1 2493.1 3890 3890 4229.1 4229.1 5956.1 5956.1 6120.4 6120.4 2679.2 3308.82 2929.2 3378.88 4146.1 4679.27 4504.99 5059.66 6034.79 6639.58 6300.33 6777.01
1840.2 1840.2 1840.2 1840.2 3463.6 3463.6 3463.6 3463.6 5038.1 5038.1 5038.1 5038.1 1861.96 1861.96 1861.96 1861.96 3487.11 3487.11 3487.11 3487.11 5200.47 5200.47 5200.47 5200.47
0.08 0.06 0.11 0.09 0.09 0.08 0.16 0.19 0.97 0.91 9.72 8.12 0.33 0.7 2.06 6.11 0.28 1.52 3.23 23.3 2.41 35.4 13.9 718
0 0 102 97 21 15 350 358 3261 2983 38542 30034 472 1877 5485 17079 257 2531 4751 52538 2212 64110 20407 914207
P11 (H ) 12)
P11 (H ) 16)
P12 (H ) 8)
P12 (H ) 12)
P12 (H ) 16)
Table 8. Computational Results for Problems P11-P12 for Makespan Maximization problem
model
|T|
discrete variables
single variables
constraints
RMIP (h)
MIP (h)
CPU s
nodes
P11 (∆i ) 2000)
CN SF CN SF CN SF CN SF CN SF CN SF CN SF CN SF
12 14 13 15 20 22 21 23 22 24 10 12 11 13 12 14
120 130 130 140 200 210 210 220 220 230 150 160 165 175 180 190
289 280 313 302 481 456 505 478 529 500 361 410 397 449 433 488
337 270 365 291 561 438 589 459 617 480 352 480 387 527 422 574
25.401 25.358 25.093 25.064 51.585 51.362 50.175 50.06 49.638 49.594 29.04 26.302 27.32 26.24 26.833 26.223
27.881 27.881 27.881 27.881 52.433 52.433 52.072 52.072 52.072 52.072 30.116 30.116 29.748 29.748 29.748 29.748
5.94 11.6 75.2 365 9.22 9.14 6.62 14.7 853 397 1.2 80.6 16.2 12776 230 60000a
13860 30837 168617 590916 24395 26500 12375 36340 1511062 892915 713 58342 9598 2732012 190282 450000
P11 (∆i ) 4000)
P12 (∆i ) 1000)
a
MRL, BPS ) 27.345.
Figure 7. Optimal schedule for problem P12 (∆i ) 1000).
M6 occurs at T2,M6 ) 6.665 h, and the values of the excess resource variables are SI1,M6,2 ) 150 and SI2,M6,2 ) 0. One could be tempted to say that between those two event points the levels of I1 have increased while those of I2 have decreased, but that is not true for I2 as it can be seen in Figure 8.
Figure 8. Inventory profiles in intermediate storage unit M6 for problem P12 (∆i ) 1000).
It is clear that there is first an increase in the amount of I2 caused by the end of its second batch in unit M1, which occurs at the 6 h mark. At that point, the level raises to 200 tonnes and stays there until the start of the first batch of I2 in stage 2,
6138 Ind. Eng. Chem. Res., Vol. 47, No. 16, 2008
which happens at T2,M3 ) 6.665. That exact same amount is then consumed, which brings the intermediate level of I2 in stage 1 down to zero. In general, we get this behavior whenever the ending time of a production task in stage k does not coincide with the starting time of the production task linked to the same product in stage k + 1. From Figure 8 we can see that this can happen relatively frequently so we cannot rely solely on variables Si,m,t to generate the intermediate storage profiles. The most important consequence from the above observation is that maximum intermediate storage levels cannot be enforced. For this example, we could have used for stage 1, values of 150 and 100 tonnes for I1 and I2, respectively, and come up with the schedule given in Figure 7, which leads to maximum storage levels that equal twice those values, as can be seen in Figure 8. Thus, the new formulation can strongly underestimate storage requirements and this is the reason why the paper has focused on unlimited intermediate storage. Modeling other storage policies like finite or no intermediate storage will be the subject of future work. 7. Conclusions This paper has presented a new continuous-time formulation for the optimal short-term scheduling of multistage, multiproduct batch plants under an unlimited intermediate storage policy. It relies on multiple time grids, one per equipment unit, where units are divided into processing and storage. One storage unit is defined for each production stage but the last, in order to model material transfers between consecutive stages of production. The use of explicit time grids for intermediate storage units is the main novelty of this paper. They open the way to deal with problems featuring multiple product batches per stage and to take advantage of batch splitting and/or mixing. The performance of the MILP formulation has been illustrated with a set of test problems taken from the literature that were divided into single and multiple batch problems. The former set had already been thoroughly examined and therefore was the obvious choice to validate the new approach. The results have shown this new approach to be competitive when compared to the multiple time grid approach of Castro and Grossmann,12 which is limited to single batch problems. In particular, it was a better performer for the objective of makespan minimization. The drawback is that the new formulation may need more event points to find the global optimal solution, which normally leads to a higher computational effort. In addition, since the event points of interacting time grids are related differently, the elegant objective function definition for total earliness minimization is no longer entirely appropriate. The new formulation has also been compared to one of the most recent versions of the unit-specific event based approach13 from the prolific group of Floudas at Princeton. For the single batch problems, the proposed approach is undoubtedly much better having been outdone in only 2 of the 24 examples. In particular, for the larger instances, it was frequent to achieve 3 to 4 orders of magnitude reduction in computational time when using the new formulation. And the reason for this was found in the typical performance indicators. The new formulation needs a smaller number of event points to find global solutions, has a lower integrality gap, and generates mathematical problems with significantly fewer constraints, for approximately the same number of binary variables. The larger number of constraints in Shaik and Floudas model13 is primarily due to the need to relate the starting times of different tasks that take place on the same unit. Such results have been confirmed in the multiple batch problems. Nevertheless, it is fair to admit that the new
formulation is limited to multistage plants with unlimited intermediate storage, while the unit-specific event based models of Floudas and co-workers can handle the more general multipurpose network structure and other storage policies. As for other limitations of the proposed approach, the characteristics of the most difficult problem instances give an idea on how large a problem can be in order to be effectively solved. We recognize that real scheduling problems may go well beyond the model’s capabilities. Efficient decomposition approaches for large-size problems with embedded sequence-based models similar to the one evaluated in this paper can be found in the literature. The same is not true for time grid based approaches so this is an obvious research challenge. Future work will widen the scope of the proposed approach not only in terms of problem size but also in terms of plant structure. Concerning the incorporation of sequence-dependent changeovers within the proposed model, which is often a very important issue in multiproduct plants, the reader is directed to the extension of this paper.27 Literature Cited (1) Prasad, P.; Maravelias, C. T. Batch Selection, Assignment and Sequencing in Multi-Stage Multi-Product Processes. Comput. Chem. Eng. 2008, 32, 1114. (2) Grossmann, I. Enterprise-wide Optimization: A New Frontier in Process Systems Engineering. AIChE J. 2005, 51, 1846. (3) Varma, V. A.; Reklaitis, G. V.; Blau, G. E.; Pekny, J. F. Enterprisewide modeling & optimization- An overview of emerging research challenges and opportunities. Comput. Chem. Eng. 2007, 31, 692. (4) Pantelides, C. C. Unified Frameworks for the Optimal Process Planning and Scheduling. In Proceedings of the Second Conference on Foundations of Computer Aided Operations; Cache Publications: New York, 1994; pp 253. (5) Kondili, E.; Pantelides, C. C.; Sargent, R. A General Algorithm for Short-Term Scheduling of Batch Operations - I. MILP Formulation. Comput. Chem. Eng. 1993, 17, 211. (6) Maravelias, C. T. Mixed-time Representation for State-Task Network Models. Ind. Eng. Chem. Res. 2005, 44, 9129. (7) Westerlund, J.; Ha¨stbacka, M.; Forssell, S.; Westerlund, T. MixedTime Mixed-Integer Linear Programming Scheduling Model. Ind. Eng. Chem. Res. 2007, 46, 2781. (8) Gupta, S.; Karimi, I. A. An Improved MILP Formulation for Scheduling Multiproduct, Multistage Batch Plants. Ind. Eng. Chem. Res. 2003, 42, 2365. (9) Harjunkoski, I.; Grossmann, I. Decomposition Techniques for Multistage Scheduling Problems using Mixed-integer and Constraint Programming Methods. Comput. Chem. Eng. 2002, 26, 1533. (10) Me´ndez, C. A.; Henning, G. P.; Cerda´, J. An MILP Continuoustime Approach to Short-term Scheduling of Resource Constrained Multistage Flowshop Batch Facilities. Comput. Chem. Eng. 2001, 25, 701. (11) Castro, P.; Erdirik-Dogan, M.; Grossmann, I. E. Simultaneous Batching and Scheduling of Single Stage Batch Plants with Parallel Units. AIChE J. 2008, 54, 183. (12) Castro, P. M.; Grossmann, I. E. New Continuous-Time MILP Model for the Short-Term Scheduling of Multistage Batch Plants. Ind. Eng. Chem. Res. 2005, 44, 9175. (13) Shaik, M.; Floudas, C. A. Unit-specific event-based continuoustime approach for short-term scheduling of batch plants using RTN framework. Comput. Chem. Eng. 2008, 32, 260. (14) Me´ndez, C. A.; Cerda´, J.; Grossmann, I. E.; Harjunkoski, I.; Fahl, M. State-of-the-art Review of Optimization Methods for Short-Term Scheduling of Batch Processes. Comput. Chem. Eng. 2006, 30, 913. (15) Castro, P. M.; Barbosa-Po´voa, A. P.; Matos, H. A.; Novais, A. Q. Simple Continuous-time Formulation for Short-Term Scheduling of Batch and Continuous Processes. Ind. Eng. Chem. Res. 2004, 43, 105. (16) Castro, P.; Barbosa-Po´voa, A. P. F. D.; Matos, H. An Improved RTN Continuous-Time Formulation for the Short-term Scheduling of Multipurpose Batch Plants. Ind. Eng. Chem. Res. 2001, 40, 2059. (17) Maravelias, C. T.; Grossmann, I. E. New General Continuous-Time State-Task Network Formulation for Short-Term Scheduling of Multipurpose Batch Plants. Ind. Eng. Chem. Res. 2003, 42, 3056.
Ind. Eng. Chem. Res., Vol. 47, No. 16, 2008 6139 (18) Sundaramoorthy, A.; Karimi, I. A. A simpler better slot-based continuous-time formulation for short-term scheduling in multipurpose batch plants. Chem. Eng. Sci. 2005, 60, 2679. (19) Ierapetritou, M. G.; Floudas, C. A. Effective continuous-time formulation for short-term scheduling. 1. Multipurpose batch processes. Ind. Eng. Chem. Res. 1998, 37, 4341. (20) Janak, S. L.; Lin, X.; Floudas, C. A. Enhanced Continuous-Time Unit-Specific Event-Based Formulation for Short-Term Scheduling of Multipurpose Batch Processes: Resource Constraints and Mixed Storage Policies. Ind. Eng. Chem. Res. 2004, 43, 2516. (21) Shaik, M.; Floudas, C. A. Improved Unit-Specific Event-Based Continuous-Time Model for Short-Term Scheduling of Continuous Processes: Rigorous Treatment of Storage Requirements. Ind. Eng. Chem. Res. 2007, 46, 1764. (22) Giannelos, N. F.; Georgiadis, M. C. A simple new continuoustime formulation for short-term scheduling of multipurpose batch processes. Ind. Eng. Chem. Res. 2002, 41, 2178. (23) Shaik, M.; Janak, S.; Floudas, C. Continuous-Time Models for Short-Term Scheduling of Multipurpose Batch Plants: A Comparative Study. Ind. Eng. Chem. Res. 2006, 45, 6190.
(24) Castro, P. M.; Grossmann, I. E.; Novais, A. Q. Two New Continuous-Time Models for the Scheduling of Multistage Batch Plants with Sequence Dependent Changeovers. Ind. Eng. Chem. Res. 2006, 45, 6210. (25) Castro, P. M.; Novais, A. Q. Optimal Periodic Scheduling of Multistage Continuous Plants with Single and Multiple Time Grid Formulations. Ind. Eng. Chem. Res. 2007, 46, 3669. (26) Sundaramoorthy, A.; Maravelias, C. T. Simultaneous Batching and Scheduling in Multistage Multiproduct Processes. Ind. Eng. Chem. Res. 2008, 47, 1546. (27) Castro, P. M.; Novais, A. Q. Scheduling Multistage Batch Plants with Sequence Dependent Changeovers. Submitted for publication in AIChE J.
ReceiVed for reView February 1, 2008 ReVised manuscript receiVed April 30, 2008 Accepted May 8, 2008 IE800194B
