1546
Ind. Eng. Chem. Res. 2008, 47, 1546-1555
Simultaneous Batching and Scheduling in Multistage Multiproduct Processes Arul Sundaramoorthy and Christos T. Maravelias* Department of Chemical & Biological Engineering, UniVersity of Wisconsin-Madison, Madison, Wisconsin 53705
In this article, we present a mixed-integer programming formulation for the simultaneous batching and scheduling in multiproduct multistage processes. The proposed sequence-based formulation addresses limitations of existing approaches where batching and scheduling decisions are carried out sequentially. To account for batching decisions, we use additional batch-selection and batch-size variables and introduce demand-satisfaction and unit-capacity constraints. Assignment constraints are active only for the subset of batches that are selected, and sequencing is carried out between batches that are assigned on the same processing unit. We also propose an alternate formulation to handle sequence-dependent changeover costs. Finally, to enhance the computational performance of the model, we present methods that allow us to fix a subset of sequencing variables and we develop a class of tightening inequalities based on time windows. 1. Introduction Scheduling is a decision-making process that concerns the allocation of limited resources to competing tasks over time with the goal of optimizing one or more objectives.1 The typical scheduling decisions are the assignment of production tasks to processing units, and the sequencing and timing of tasks on each unit. In the chemical industry, most batch facilities are multistage processes, where batches are processed in a sequence of stages without splitting, merging, or recycling. Multiproduct batch scheduling applications are prevalent in the production of pharmaceuticals, specialty chemicals, food processing, and microelectronics industries where product diversification and demand fluctuations are common. Although several methods have been proposed over the last two decades for the scheduling of such processes,2-5 there are classes of problems that have not been addressed effectively. Specifically, in existing approaches batching decisions are typically decoupled from scheduling decisions, thus often leading to suboptimal solutions. Further, existing methods consider fixed processing times, an assumption that is not always valid if batching and scheduling are considered simultaneously. The goal of this article is the development of a new mixedinteger programming (MIP) formulation for the simultaneous batching and scheduling of multiproduct multistage batch processes under variable processing times. The article is organized as follows. In section 2, we review multistage scheduling methods in the process engineering literature, provide motivation for our work, and state the problem at hand. In section 3, we present the mathematical formulation, followed by discussions on an alternate formulation for changeover costs and the development of tightening constraints in section 4. In section 5, we discuss how the proposed formulation can be simplified to solve standalone scheduling problems. Finally, five illustrative examples are presented in section 6, and computational results are presented in section 7. 2. Background and Problem Description 2.1. Multistage Batch Scheduling. In multistage batch processes, customer orders are usually divided into several * To whom correspondence should be addressed. Tel.: +1-608- 2659026. Fax: +1-608-262-5434. E-mail:
[email protected].
batches, that are processed through all stages sequentially, where each stage has one or more parallel units. The sequential processing of batches implies that batches can be treated as discrete entities moving through the process. Many researchers6-15 have exploited this structure and proposed mathematical formulations that are less general but computationally more effective than network-based formulations. Existing MIP formulations can be tagged on the basis of two approaches of modeling batch sequencing: (i) slot-based and (ii) sequence-based. In slot-based formulations,6,7,15 the available production time of each unit is divided into a number of ordered slots of unknown length. Batches are allotted to these slots, thus enforcing sequencing and timing restrictions. The solution time apparently depends on the number of slots postulated for each unit. On the other hand, sequence-based formulations8-12 use sequencing binaries to achieve the sequencing and timing of batches. Sequence-based formulations can further be classified into formulations with global or local (immediate) sequencing variables. Finally, Castro and Grossmann13 and Castro et al.14 presented multiple-time-grid continuous-time formulations based on resource-task-network representation.16 To solve large instances, Harjunkoski and Grossmann17 and Maravelias18 developed hybrid solution techniques. While in both approaches the original problem is decomposed into assignment and sequencing subproblems, the former uses constraint programming to solve the sequencing subproblem, whereas the latter uses special-purpose algorithms for the same. Reviews on multistage batch scheduling are given by Pekny and Reklaitis,3 Pinto and Grossmann,2 Shah,4 Kallrath,19 and Me´ndez et al.5 2.2. Simultaneous Batching and Scheduling. In the context of multistage batch processes, the batching and scheduling are viewed as two separate problems. Orders are first divided into batches based on the unit capacities, and these batches are used as inputs in the scheduling model. This two-step approach is effective when the demand is fixed and the parallel units have similar capacity at each stage. However, situations may often arise where the objective is the maximization of production or sales in addition to meeting customer demand for a given scheduling horizon. Further, capacity expansions and/or retrofit projects could potentially lead to processes with parallel units of dissimilar capacities. Hence, there exist multiple ways an order can be divided into batches. Therefore, batching and
10.1021/ie070944y CCC: $40.75 © 2008 American Chemical Society Published on Web 02/09/2008
Ind. Eng. Chem. Res., Vol. 47, No. 5, 2008 1547
Figure 1. Multistage batch plant with parallel nonidentical units per stage processing multiple orders of different release and due times.
scheduling decisions should be considered simultaneously.20,21 Furthermore, one common assumption in all existing multistage optimization methods is that the processing times are constant irrespective of the corresponding batch sizes. If batch sizes are fixed, the assumption is valid. However, if batching is a part of the optimization, then the processing times can also vary with the batch sizes. To address the above-mentioned limitations, we present a general global-sequence MIP formulation for the simultaneous batching and scheduling in multistage multiproduct processes with variable processing times. 2.3. Problem Statement. Given a multistage multiproduct batch plant (Figure 1) with: (i) a set of orders i∈I with demand qi, release/due time ri/di, and price πi; (ii) a set of units j∈J with minimum/maximum capacity bjmin/ bjmax, fixed/proportional processing time constants τijF/τijP, and fixed/proportional processing cost constants cijF/cijP; (iii) a set of stages k∈K with parallel (identical or nonidentical) units j∈Jk at stage k (J ) J1 ∪ J2 ...∪ J|K|); (iv) sets of forbidden units JFi for order i and forbidden paths (j,j′)∈FP for all orders (j∈JAik ) Jk\JFi is the set of allowable units of order i in stage k); ch (v) changeover time τii′j from a batch of order i to a batch of order i′ on unit j; the goal is to determine: (i) the number and size of batches for each customer order; (ii) the assignment of batches to units; (iii) the sequencing of assigned batches in each unit. We further assume that (a) preemption is not allowed, (b) all batches visit all stages, and (c) there is unlimited storage for intermediates between stages. Note that an order can have multiple products and a product can be demanded by multiple orders. However, we use index i to denote an order of a particular product. If an order has more than one products, then it is divided into multiple single-product orders with the same release and due times. 3. Mathematical Formulation The proposed formulation has three levels in the hierarchy of decisions: (i) the selection and sizing of batches, (ii) the assignment of batches to processing units, and (iii) the sequencing and timing of batches in each unit. The first level corresponds to the batching problem, whereas the second and the third levels comprise the scheduling problem.
Figure 2. Calculation of parameters bˆ max ˜ max ˆ min for order i. i , b i , and b i
3.1. Number of Batches. To enable the simultaneous optimization of batching and scheduling decisions, we introduce index l to denote the lth batch for order i∈I. The number of batches required to meet an order i depends on the demand qi and the capacities of units that can process the order. To this max end, the minimum (lmin i ) and maximum (li ) numbers of batches for order i are calculated using eqs 1 and 2:
limin ) qi/bˆ imax
∀i
(1)
limax ) qi/bˆ imin
∀i
(2)
) min[max(bmax ˆ min ) max[min (bmin where bˆ max i j )] and b i j )] are k∈K j∈JAik
k∈K j∈JAik
the maximum and minimum feasible batch sizes for order i respectively (Figure 2). Further, we define a set Li ) {1, 2,... lmax i } of potential batches for order i∈I. Note that, when the processing times are constant, limax can be calculated by,
limax ) qi/b˜ imax
∀i
(2a)
where b˜ max ) min[min (bmax i j )] is the largest batch size for k∈K j∈JAik
order i that can be processed in every allowable unit (Prasad and Maravelias21). This is because there is no incentive to divide a batch that can be processed on one unit into two or more batches when processing times are fixed. If processing times are variable, however, the division of an order into batches of size close to can potentially lead to a better solution. Therefore, eq 2 bˆ min i provides the absolute upper bound on the number of batches for order i. Nevertheless, eq 2a can be used to provide a more realistic estimate of lmax i , which is typically sufficient to obtain the optimal solution as will be discussed in section 7.
1548
Ind. Eng. Chem. Res., Vol. 47, No. 5, 2008
3.2. Batching and Assignment. If Bil denotes the size of the lth batch of order i∈I (hereafter referred to as batch (i,l)), then the sum of all batch sizes must be greater than or equal to the demand qi.
∑ Bil g qi
∀i
If Zil denotes the selection of batch (i,l), then batch size Bil is constrained by,
∀i,l∈Li
(4)
If batch (i,l) is selected for an order i, then it must be assigned to only one unit in each stage,
Zil )
∑
j ∈JAik
Xilj
∀i,l ∈Li,k
∑
j ∈JAik
B h ilj
ik
F P ch (τ i′j Xi′l′j + τ i′j B h i′l′j + τ ii′j Xi′l′j) i′k
H(1 - Yili′l′k)
∀ (i,l,i′l′) ∈ IL,k (10a)
∀i,l ∈Li,k
Note that eq 10a can be used under the assumption that the changeover time between any two orders i and i′ is such that ch ch min ch min F P min e τii′′j + τi′′j + τi′′i′j , ∀ i′′,j, where τi′′j ) τi′′j + τi′′j bˆ i′′ is τii′j the minimum processing time of order i′′ in unit j. The difference in the finish times of batch (i,l) at any two successive stages must be greater than or equal to its processing time at the latter stage.
Til(k + 1) g Tilk +
(6)
∑
(τ Fij Xilj + τ Pij B h ilj)
j∈JAi(k + 1)
(5)
where binary variable Xilj denotes the assignment of batch (i,l) to unit j. To model variable processing times, we disaggregate batch h ilj for the size of batch (i,l) processed size Bil into variables B on unit j.
Bil )
∑ j∈JA ∩ JA
(3)
l∈Li
Zil e Bil e bˆ max Zil bˆ min i i
Ti′l′k g Tilk +
∀ i,l ∈ Li,k < |K| (11) 3.4. Additional Constraints. The release time ri and due time di of each order should be respected,
Tilk g riZil + Tilk e diZil -
∑ ∑
(τ Fij Xilj + τ Pij B h ilj)
∑ ∑ k′>k j∈JA
(τ Fij Xilj + τ Pij B h ilj) + LTil
k′ek j ∈ JAik′
ik′
∀ i,l ∈ Li,k (13)
In addition, the size of a batch should satisfy the capacity limits on any assigned unit.
Xilj e B h ilj e bmax Xilj bmin j j
∀i,l ∈Li,j
(7)
Because batch (i,l) is assigned to only one unit in stage k (from eq 5), only one disaggregated variable B h ilj is allowed to be nonzero. Also, note that eq 4 becomes redundant in the presence of eqs 5-7. Further, numerical ordering in the selection of batches is enforced for symmetry breaking purposes.
Zi(l + 1) e Zil
∀i,l ∈ Li
(8)
3.3. Sequencing and Timing. Variable Yili′l′k is used to denote the pairwise sequencing between batches (i,l) and (i′,l′); that is, it is equal to 1 if batch (i,l) is processed before (but not necessarily immediately before) batch (i′,l′) on a unit at stage k. Equation 9 activates a sequencing binary when two batches are processed on the same unit in stage k,
Xilj + Xi′l′j - 1 e Yili′l′k + Yi′l′ilk ∀ (i,l,i′l′) ∈IL: i e i′,k,j ∈JAik ∩ JAi′k (9) where IL ) {i,i′,l ∈ Li,l′ ∈ Li′ : (i * i′) ∨ ((i ) i′) ∧ (l * l′))} is the set of all pairs of batches (i,l) and (i′,l′) that can be sequenced on a unit. Note that eq 9 is expressed only for i e i′ because both Yili′l′k and Yi′l′ilk are included. If batch (i,l) precedes batch (i′,l′), then the finish time Ti′l′k of batch (i′,l′) in stage k is constrained by,
Ti′l′k g Tilk +
∑
F (τ i′j Xi′l′j
+
P τ i′j B h i′l′j)
- H(1 - Yili′l′k)
j∈JAi′k
∀ (i,l,i′l′) ∈ IL,k (10) where the processing time of batch (i′,l′) has fixed and proportional terms. When sequence-dependent changeover times are present, eq 10 becomes
∀ i,l ∈ Li,k (12)
where LTil is a positive penalty variable if an order is satisfied after its due date. Note that LTil can be fixed to zero if the violation of due dates is not acceptable. Equation 13 is not needed if no due dates for orders exist. Equations 12 and 13 also serve as lower and upper bounds, respectively, on the finish times of batches at each stage. The constraints in eqs 14 and 15 ensure that forbidden paths and assignments are not violated.
∀ i,l ∈ Li,(j,j′) ∈ FP
Xilj + Xilj′ e Zil Xilj ) 0
∀ i,l ∈ Li,j ∈ JFi
(14) (15)
Further, eq 16 is added for symmetry breaking purposes.
Bi(l+1) e Bil
∀ i,l ∈ Li
(16)
For each order i, the selection binaries Zil are fixed to 1 for lelmin i .
Zil ) 1
∀ i,l e lmin i
(17)
Finally, redundant variables (i.e., variables for l ∉Li and (i,l,i′,l′) ∉IL) are fixed to zero via eq 18, and integrality and non-negativity constraints are expressed in eq 19.
Zil, Xilj, Bil, B h ilj, Tilk, LTil ) 0 ∀i,l ∉Li,j,k ∀ (i,l,i′l′) ∉ IL,k (18) Yili′l′k ) 0 Zil, Xilj, Yili′l′k ∈ {0,1}
Bil, B h ilj, Tilk, LTil g 0
(19)
Equations 3 and 5-19 form the general global-sequence MIP model (M1) for the simultaneous batching and scheduling of multistage multiproduct batch plants with variable processing times.
Ind. Eng. Chem. Res., Vol. 47, No. 5, 2008 1549
3.5. Objective Function. The model is flexible to accommodate a range of objective functions. The minimization of lateness, earliness, and processing costs are given in eqs 20, 21, and 22, respectively.
∑ LTil
(20)
∑ (di - Til|K|)
(21)
Figure 3. Differences between global- and local-sequence variables.
∑ (c Fij Xilj + c Pij Bh ilj) i,l∈L ,j∈JA
(22)
XFilj and XLilj that take the value of 1 if batch (i,l) is assigned first and last to unit j, respectively (Figure 3).
min
i,l∈Li
min
i,l∈Li
min
i
i
MS g Til|K| ∀i,l ∈ Li
(23)
min MS
(24)
∑
Wi′l′ilj + XFilj ) Xilj
∀i,l ∈ Li,j
(9a)
∑
Wili′l′j + XLilj ) Xilj
∀i,l ∈ Li,j
(9b)
i′,l′∈Li′
where JAi is the set of allowable units for order i. For the minimization of makespan MS in eq 24, we should enforce that all tasks are completed before MS.
i′,l′ ∈Li′
At most one order can be the first (last) in the sequence on each unit.
Finally, the objective of profit maximization is given in eq 25.
max
∑ πiBil - i,l∈L∑,j∈JA
i,l∈Li
i
(c Fij Xilj
+
c Pij B h ilj)
Ti′l′k g Tilk +
∀i
(3a)
i
The upper bound on production is also used to calculate the maximum number of batches. Equation 2, therefore, is replaced by eq 2b.
limax ) qimax/bˆ imin
(9c)
∑ XLilj e 1
∀j
(9d)
The sequencing and timing of batch (i′,l′) that immediately follows batch (i,l) on unit j is achieved through eq 10b.
4.1. Maximum Demand. The solutions to problems where the objective is the maximization of profit tend to favor the production of few products. Although mathematically correct, these solutions do not accurately represent reality, where more balanced production schedules are desirable. To address this shortcoming, we can introduce an upper bound qmax on the i production of product i. In this case, the total production of product i should satisfy eq 3a instead of eq 3.
Bil e qimax ∑ l∈L
∀j
i,l∈Li
(25)
i
4. Extensions and Computational Enhancements
qi e
∑ XFilj e 1
i,l∈Li
∀i
ik
F P ch (τ i′j Xi′l′j + τ i′j B h i′l′j + τ ii′j Wili′l′j) - H i′k
(1 -
∑
j∈JAik∩JAi′k
Wili′l′j)
∀(i,l,i′l′) ∈ IL,k (10b)
Integrality and non-negativity constraints are expressed in eq 19a.
Zil, Xilj, Wili′l′j ∈{0,1} XFilj, XLilj ∈[0,1] Bil, B h ilj, Tilk, LTil g 0 (19a) Equations 3, 5-8, 9a-9d, 10b, 11-18, and 19a constitute the alternate local-sequence MIP formulation (M2) for the simultaneous batching and scheduling of multistage multiproduct batch plants with variable processing times and sequencedependent changeover times and costs. The objectives of cost minimization and profit maximization can now be expressed via eqs 22a and 25a respectively,
(2b)
4.2. Sequence-Dependent Changeovers. Global-sequence models are on average better than local-sequence or slot-based models in terms of computational performance. However, these models cannot be used to address problems with sequencedependent changeover costs because the global-sequence variables are active for all of the pairs of batches assigned on the same unit. To address this limitation, we develop an alternate formulation for the sequencing and timing of batches in the same unit (eqs 9 and 10). We first define a new local-sequence variable Wili′l′j that is equal to 1 if batch (i,l) immediately precedes batch (i′,l′) on unit j. Unlike global-sequence variables Yili′l′k, the local-sequence variables Wili′l′j are active only for adjacent batches that are assigned to a specific unit j∈Jk. Further, we define variables
∑ j∈JA ∩JA
min
(c Fij Xilj + c Pij B h ilj) + ∑ ∑ i,l∈L ,j
max
γii′jWili′l′j (22a)
(i,l,i′,l′)∈IL,j
i
πiBil - ∑ (c Fij Xilj + c Pij B h ilj) - ∑ ∑ i,l∈L i,l∈L ,j i
i
(i,l,i′,l′)∈IL,j
γii′jWili′l′j (25a)
where γii′j is the changeover cost between orders i and i′ on unit j. 4.3. Fixing of Sequencing Binaries. To enhance the computational performance of the proposed models, we present some tightening procedures based on how early (late) a batch can start (end) at each stage. Let estik, eftik, lstik, and lftik denote the earliest start time, earliest finish time, latest start time, and latest finish time of any batch of order i in stage k, respectively (Appendix A). Then, some of the sequencing (local or global) variables can be fixed to zero based on the values of eftik and
1550
Ind. Eng. Chem. Res., Vol. 47, No. 5, 2008
lstik.
ened local-sequence MIP formulation (M2*) consists of eqs 3, 5-8, 9a-9d, 10b, 11-18, 19a, 26a, 27a, and 28b.
∀(i,l,i′,l′) ∈ IL,k : eftik > lsti′k
Yili′l′k ) 0 Wili′l′j ) 0
∀(i,l,i′,l′) ∈ IL,j ∈ Jk, k : eftik > lsti′k
(26) (26a)
Further, given a set of feasible predecessors14 Pi′k (Appendix A1), the variables Yili′l′k and Wili′l′j of orders that are not in the set for order i′ can be fixed to zero.
Yili′l′k ) 0 Wili′l′j ) 0
∀(i,l,i′l′) ∈ IL,k,i ∈/ Pi′k
(27)
∀(i,l,i′l′) ∈ IL,j ∈ Jk, k,i ∈/ Pi′k
(27a)
4.4. Tightening Inequalities. Given the release/due times of a subset of orders and the corresponding processing times, one can determine time windows (defined by the latest finish time and earliest start time among all the orders of the subset) within which all batches of this subset should be processed on each unit. Clearly, in any feasible solution the sum of the processing times of the batches that are assigned to a unit cannot exceed the length of the time window (assuming LTil ) 0).
(τ Fij Xilj + τ Pij B h ilj) e ubwjn - lbwjn ∑ i∈IC ,l∈L jn
∀j,n ∈ Nj (28)
i
where ICjn is the subset of orders considered for the nth time window in unit j, ubwjn (lbwjn) is the upper (lower) bound on the time window for the corresponding set of orders, and Nj is the set of inequalities generated for unit j. The details of the calculation of the time windows are given in Appendix B. Note here that we calculate unit-specific time windows taking into account forbidden unit-order assignments and not stage-specific time windows as in Maravelias18 and Prasad and Maravelias.21 If the objective is the minimization of makespan, then eq 28 should be modified as follows.
(τ Fij Xilj + τ Pij B h ilj) e MS - tailjn - lbwjn ∑ i∈IC ,l∈L jn
i
∀k,j ∈ Jk,n ∈ Nj (28a)
where tailjn is the minimum time required for an order i∈ICjn to be completed after being processed at stage k (Appendix A). In the presence of sequence-dependent changeover times, the local-sequence formulation allows us to add changeover terms in the LHS of eq 28.
(τ Fij Xilj + τ Pij B h ilj) + ∑ ∑ i∈IC ,l∈L jn
i
(i,l,i′,l′)∈IL, i,i′∈ICjn
ubwjn - lbwjn
ch τ ii′j Wili′l′j e
∀k,j ∈ Jk,n ∈ Nj (28b)
It is important to note here that although the global-sequence formulation handles changeover times, it cannot be used for the development of the valid inequalities because additional changeover terms will be incorrectly included in the second summation in the LHS of eq 28b. Thus, in addition to the modeling of changeover costs, the local-sequence model enables us to derive these tightening inequalities that, as we will show later, are very effective. The strengthened global-sequence MIP formulation (M1*) consists of eqs 3, 5-19, 26, 27, and 28, whereas the strength-
5. Scheduling Formulation Clearly, optimization decisions in scheduling problems (batchunit assignment, sequencing, and timing) are a subset of the optimization decisions in the proposed batching-scheduling model. Hence, the proposed model can be reduced to a scheduling model for a fixed set of batching decisions. If M is the fixed set of batches to be processed, and Bm the size of batch m∈M, then this is achieved as follows: (i) Set Definition. For every batch m∈M, we define an order i∈I and set Li ) {1} ∀i. Because all variables and constraints are defined only for l ) 1, this is equivalent to dropping index l from all variables and constraints. Sets J and K remain the same. (ii) Variable Fixing. Because the size of batch m is fixed to Bm, we can fix variables Bil ) Bm ∀i. In addition, the variable h ilj) can be replaced by a fixed processing time (τFijXilj + τPijB order-unit specific processing time τijXilj, where τij ) τFij + τPij Bm. The new term can be used to replace the variable processing ) lmax ) time in timing constraints (eqs 10-13). Because lmin i i 1, eq 17 implies that all variables Zil are fixed, that is, variables Zil can be replaced by 1 in all constraints. (iii) Constraint Elimination. Obviously, the selected batches and corresponding sizes should meet the demand. Thus, eq 3 can be removed. Because batch sizes are known, eq 4 can be used prior to optimization to find what units a batch can be assigned to. If the capacity bounds of a unit j are not appropriate for batch i, then unit j can be removed from set Ji of units suitable for order i, and the corresponding eq 7 can be removed. Equation 8 is also removed because variables Zil are fixed to 1. Finally, eq 16 is removed and eq 17 is used to fix Zil to 1. Thus, if we drop index l for simplicity (because l ) 1 in all variables and constraints), the reduced scheduling formulation for a given set of batching decisions consists of eqs 29-35.
∑
∀i,k
Xij e 1
(29)
j∈JAik
∀i e i′,k,j ∈ JAik ∩ JAi′k (30)
Xij + Xi′j - 1 e Yii′k + Yi′ik Ti′k g Tik +
∑
∀i * i′,k
τi′jXi′j - H(1 - Yii′k)
(31)
j∈JAi′k
Ti′k g Tik +
∑ j∈JA ∩JA ik
ch (τi′jXi′j + τ ii′j Xi′j) - H(1 - Yii′k) i′k
∀i * i′,k (31a)
Ti(k + 1) g Tik +
∑
τijXij
∀i,k < |K|
(32)
j∈JAi(k + 1)
Tik g ri +
∑ ∑
τijXij
∀i,k
(33)
k′ek j∈JAik′
Tik e di -
∑ ∑
τijXij + LTi
∀i,k
(34)
k′>k j∈JAik′
Xij + Xij′ e 1
∀i,(j,j′) ∈ FP
(35)
where set JAik is updated to contain only units with appropriate e Bm e bmax capacity bounds, that is, bmin j j . If only the number of batches but not the batch sizes is fixed, then the scheduling problem can be solved by setting lmin ) i ) 1, fixing Z ) 1 ∀i, l∈L and Z ) 0 ∀i, l ∉ L . lmax il i il i i
Ind. Eng. Chem. Res., Vol. 47, No. 5, 2008 1551 Table 1. Process Data for Examples 1-2a example 1 bmin j
example 2
bmax j
bmin j
bmax (kg) j
stages
units
1
J1 J2
10 20
30 40
10 20
30 40
2
J3 J4 J5
20 25
35 50
20 15 15
35 30 25
a
(kg)
(kg)
(kg)
Forbidden Assignments: JFC ) {J3} for example 2.
6. Examples In this section, we consider five examples. The first example illustrates the advantage of the simultaneous optimization of batching and scheduling decisions over the sequential approach. The second one illustrates the need to model variable processing times when batching and scheduling are considered simultaneously. Examples 3-5 highlight the application and performance of our models. All examples were modeled using GAMS 22.4 and solved using CPLEX 10.1 on a desktop computer with a 2.8 GHz Pentium 4 processor and 512 MB RAM running Windows XP. The data for the first two examples are given in Tables 1 and 2, whereas those for the last three examples are given in Supporting Information. We solved examples 1-4 using model M1* and example 5 using M2*. In all of the examples, we used lmax calculated by eq 2a. i 6.1. Example 1. Three orders have to be processed through two stages with two units per stage. Order A has a demand of 30 kg, and orders B and C have demand of 40 kg each. The objective is to minimize the makespan. The processing time depends on the unit and varies with batch size. Note that the bigger orders with 40 kg of demand can be processed as single batches on units J2 and J4. Because all orders are within capacity limits, one would be tempted to assign a priori each order to one batch and then solve for an optimal schedule. Figure 4a shows the Gantt chart of the solution obtained from such a sequential approach (example 1a) with a makespan of 17.2 h. We then solved the same instance without fixing any Zil (example 1b) in model M1*. The Gantt chart of the obtained solution, with a makespan equal to 14.5 h, is shown in Figure 4b. Because unit J2 of stage 1 is the bottleneck, its load is shared by unit J1 by dividing order B into two batches. Although one could argue that such division of orders is intuitive in sequential batching methods, which order to split, and in what sizes are optimization decisions that can easily become non-intuitive in larger instances. Note that the makespan in part b of Figure 4 is smaller, even though the total processing time of order B increases. A similar observation can also be made when we use lateness minimization for this example. 6.2. Example 2. This example considers three orders that are processed through two stages with two units in the first stage and three units in the second. The demands for orders A, B,
Figure 4. Gantt charts for example 1.
Figure 5. Gantt charts for example 2.
and C are 30, 70, and 55 kg, due at 10, 25, and 20 h, respectively. The objective is the minimization of earliness. First,
Table 2. Demands, Release/Due Times, and Constants for Processing Times for Examples 1 and 2 fixed term τFij (h) orders
qi (kg)
ri (h)
A B C
30 40 40
0 0 0
A B C
30 70 55
0 0 0
di (h)
10 25 20
J3
proportional term τPij (h/kg)
J1
J2
J4
2.500 2.500 2.500
2.000 2.000 2.000
Example 1 0.889 2.000 0.889 2.000 0.889 2.000
2.500 2.500 2.500
2.000 2.000 2.000
Example 2 1.111 2.667 1.111 2.667 2.667
J5
1.333 1.333 1.333
J1
J2
J3
J4
J5
0.083 0.083 0.083
0.100 0.100 0.100
0.089 0.089 0.089
0.080 0.080 0.080
-
0.083 0.083 0.083
0.100 0.100 0.100
0.111 0.111
0.178 0.178 0.178
0.267 0.267 0.267
1552
Ind. Eng. Chem. Res., Vol. 47, No. 5, 2008
Figure 6. Gantt chart of optimal solution for example 3.
Figure 7. Gantt chart of optimal solution for example 4.
we solved the example assuming that the processing time is irrespective of the batch size (example constant as τFij + τPijbmax j 2a). The Gantt chart of the solution is shown in Figure 5a with an earliness of 5 h. Note that both batches of order B take place on unit J3 even though units J4 and J5 are idle from 20 to 25 h. This is because the batch sizes of order B are greater than the capacities of units J4 or J5. Even if one of the batches is divided into two smaller batches and processed on units J4 and J5, the earliness cannot be reduced from 5 h because the constant processing time of order B in both of the units is 8 h. When we solved the same instance using variable processing times (example 2b), we obtained a solution with an earliness of 1.56 h (Figure 5b). In this case, the second batch of order B is divided into two batches (with shorter processing times) and then assigned to units J4 and J5, leading to a better solution. 6.3. Example 3. In this example, we consider 12 orders with release and due times and a batch plant with two stages and three units per stage. The objective is to meet customer demands at the minimum production cost. The Gantt chart of an optimal solution is shown in Figure 6. In the optimal solution, the 12 orders are divided into 24 batches. The model was solved to optimality in 320.6 CPU s. 6.4. Example 4. Because batching decisions were parted from the optimization model, previously proposed methods did not consider profit maximization until very recently.21 In this example, we consider the same problem instance for profit maximization as in Prasad and Maravelias.21 In addition, we introduce maximum demands for all orders while we treat the
Figure 8. Gantt chart of optimal solution for example 5. Gray boxes represent changeover times between batches of different orders. Table 3. Model Statistics for Examples 3-5 model M1/M2
model M1*/M2*
variables
variables
example binary continuous constraints binary continuous constraints 3 4 5a 5b
1680 1253 571 571
249 214 169 169
4529 3437 572 572
1364 1253 571 571
249 214 169 169
4712 3443 596 596
original order quantities as minimum demands. The calculation of a minimum and maximum number of batches for each order is now straightforward using eqs 1 and 2b. The objective is to maximize the net profit (reVenue from sales - total production cost). The Gantt chart of an optimal solution for a 60 h horizon is given in Figure 7. A total of 20 batches are processed over the given horizon. The optimal solution of $2,056 was obtained in 228.5 CPU s.
Ind. Eng. Chem. Res., Vol. 47, No. 5, 2008 1553 Table 4. Solution Statistics for Examples 3-5 model M1/M2
model M1*/M2*
example
optimal objective
objective value
LP-relaxation
CPU (s)
nodes
objective value
3 4 5a 5b
3037.00 2056.00 2960.00 2960.00
3037.00 2025.00 2960.00 2960.00
3028.00 2234.67 1277.21 1277.21
10 000.0a 10 000.0b 30.6 78.2
385 479 846 091 9757 22 914
3037.00 2056.00 2960.00 2960.00
a
LP-relaxation 3028.00 2168.51 1277.21 1277.21
CPU (s)
nodes
320.6 228.5 5.5 5.7
23 486 20 918 1366 1321
0.22% gap after 10 000 CPU s. b 3.16% gap after 10 000 CPU s.
Table 5. Comparison of Solution Statistics for the Examples min{lmax +1, lmax* } i i
lmax i example
optimal objective
LP-relaxation
CPU (s)
nodes
LP-relaxation
CPU (s)
nodes
1a 1b 2a 2b 3 4 5a 5b
17.20 14.50 5.00 1.56 3037.00 2056.00 2960.00 2960.00
14.67 11.35 0.00 0.00 3028.00 2168.51 1277.21 1277.21
0.1 0.4 0.1 0.8 320.6 228.5 5.5 5.7
0 614 79 595 23 486 20 918 1366 1321
14.67 11.35 0.00 0.00 3028.00 2168.51 1271.04 1271.04
0.1 0.8 0.3 1.5 3407.0 2315.4 9928.2 9550.8
0 1,328 121 658 114 767 65 831 1 408 252 1 496 356
6.5. Example 5. Finally, we illustrate the application of model M2* to cases where sequence-dependent changeover costs are significant. We consider a batch plant that processes six orders through two stages with two units per stage. In addition to the demand, due and release times, prices, and costs, we are also given sequence-dependent changeover times and costs. The objective is to minimize the sum of production and changeover costs while meeting customer demands. The Gantt chart of the minimum cost solution is shown in Figure 8, where nine batches are processed to meet the given orders. The model was solved in 5.5 CPU s to give an optimal cost of $2,960. 7. Computational Results In this section, we discuss the effectiveness of the tightening methods presented in sections 4.3 and 4.4 and the effect of the on solution times. To study the effect of choice of lmax i tightening constraints, we solved examples 3-4 using models M1 and M1* and example 5 using models M2 and M2*. To analyze the impact of changeover times, we consider two ch instances of example 5: one with the changeover times τii′j given in the Supporting Information (example 5a), and the other ch (exwith changeover times increased by 50%, i.e. 1.5* τii′j ample 5b). All other data in examples 5a and 5b are identical. The model and solution statistics are given in Tables 3 and 4, respectively. In terms of model statistics, we first note that model M1* has fewer variables than model M1 in example 3. This is because the time windows in the example are tighter, thus leading to variable fixing via the methods described in subsection 4.3. Second, the number of constraints in the strengthened formulations M1*/M2* is larger due to the addition of the valid inequalities (section 4.4). For the maximization of profit, we note that the addition of one valid inequality for set I for each unit reduces the computational requirements, whereas the addition of valid inequalities for all possible subsets of orders does not yield significantly better formulations. Thus, the number of additional constraints in example 4 is only six. In terms of computational statistics, we first note that the proposed methods strengthen the MIP formulation: in example 4 (profit maximization), the integrality gap is reduced from 8.69 to 5.47%. Further, the proposed methods improve the computational performance of the formulation significantly. Examples 3 and 4 were not solved to optimality even after 10 000 CPU s
using models M1 and M2, whereas they were solved to optimality in 320.6 and 228.5 CPU s, respectively, using the strengthened formulations. We also observe that the reduction in solution time due to the addition of the tigthening constraints is larger in example 5b than in example 5a. This is because changeover times are 50% larger in example 5b, which means that the inequalities in eq 28b are tighter. In other words, the proposed computational enhancements increase as the problem becomes more tightly constrained. The solutions presented in section 6 were obtained using calculated by eq 2a. Although these solutions are optimal, lmax i in this section we study how the solution time increases as we choose larger values for lmax i . As explained in section 3.1, the largest possible number of batches can be calculated by eq 2; let this number be denoted by lmax* . If the number of batches i in the optimal solution of an for order i is equal to lmax i new equal to example, then we re-solve this example with (lmax i ) max max* min{li +1, li }. Table 5 provides a comparison of solution and min{lmax + 1, statistics between the examples with lmax i i max* li } using models M1*/M2*. First, we observe that the calculated by eq 2a) solutions obtained in section 6 (using lmax i are optimal. Second, we observe that the computational requirements are significantly larger if the number of potential batches is increased. Therefore, it is important to estimate a tight but sufficient limax. 8. Conclusions The simultaneous optimization of batching and scheduling decisions in multistage processes has received little attention in the literature, although it is clear that it can potentially lead to solutions that are better than the ones obtained by sequential approaches. The challenge in addressing this problem is the modeling of batch-selection and sizing decisions, which in turn leads to variable processing times, a feature that current formulations do not account for. In this article, we presented an MIP formulation for the simultaneous batching and scheduling in multistage, multiproduct batch processes with variable processing times. Batching decisions are modeled via the introduction of additional selection binary variables and continuous sizing variables. Scheduling decisions are modeled via assignment and sequencing variables and constraints that are activated only if the corresponding batches are selected.
1554
Ind. Eng. Chem. Res., Vol. 47, No. 5, 2008
Processing times are modeled as functions of selection and batch-size variables. For the sequencing of batches on the same unit, we used a global-sequence formulation, which was found to perform better, on average, than slot-based or local-sequence formulation. Nevertheless, we developed an alternate localsequence formulation that accounts for sequence-dependent changeover costs. Furthermore, we showed how the proposed model can be used to address standalone scheduling problems. To address larger instances we also presented methods that allow us to fix certain sequencing variables and developed strengthening inequalities based on time-window information. The proposed methods were shown to yield significant improvements, enabling the solution of instances of medium size. Another advantage of the proposed approach is its inherent structure: selection binaries activate assignment binaries, which in turn activate sequencing binaries. In other words, information flows from the higher (selection) level to the lower (sequencing) level. This structure can be exploited in decomposition methods, where the original problem is decomposed into several subproblems that can be solved in parallel. Nomenclature Indices/Sets i,i′,i′′∈I ) orders j,j′∈J ) units k,k′∈K ) stages l,l′∈Li ) batches n∈Nj ) tightening inequalities for unit j m∈M ) batches in standalone scheduling problems Subsets JFi ) forbidden units for order i FP ) forbidden path for orders between units j and j′ IAj ) orders that are allowed to process in unit j IL ) pairs of batches (i,l) and (i′,l′) that can be sequenced ICjn, IC′ ) set of orders used in the preprocessing algorithm JAik ) units in stage k that are allowed to process order i JAi ) units that are allowed to process order i Jk ) units in stage k Kj ) stage of unit j Pi′k ) feasible predecessors for order i′ in stage k
ubii′j ) upper bound on start time of order i assuming that order i′ follows i in unit j πi ) price of product associated with order i τFij/τPij ) fixed/proportional constant of the processing time of order i in unit j τmin ij ) minimum processing time of order i in unit j ch ) changeover time for transition from a batch for i to a τii′j batch for i′ in unit j γii′j ) changeover cost between orders i and i′ in unit j Binary Variables Wili′l′j ) sequencing of batch (i,l) immediately before (i′,l′) on unit j Xilj ) assignment of batch (i,l) to unit j Yili′l′k ) sequencing of batch (i,l) before (i′,l′) on the same unit in stage k Zil ) selection of batch l of order i PositiVe Variables Bil ) size of batch (i,l) B h ilj ) size of disaggregated batch (i,l) in unit j LTil ) slack variable for late completion of batch (i,l) MS ) makespan Tilk ) finish time of batch (i,l) in stage k XFilj∈[0,1] ) assignment of batch (i,l) first in the sequence in unit j XLilj∈[0,1] ) assignment of batch (i,l) last in the sequence in unit j Appendix A. Parameters for Fixing Binaries and Feasible Predecessors. The earliest start time (estik), earliest finish time (eftik), latest start time (lstik), and latest finish time (lftik) of any batch of order i in stage k are calculated as shown below:
estik ) ri + eftik ) ri +
min (τ min ij )
(A1)
min min (τ min ij ) ) estik + min (τ ij ) (A2)
k′ e k j∈JAik′
j∈JAik
lstik ) di -
min (τ min ∑ ij ) j∈JA k′ g k
(A3)
ik′
lftik ) di -
Parameters max bmin ) minimum/maximum allowable batch size in unit j j /bj max bˆ i ) maximum feasible batch size of order i bˆ min ) minimum feasible batch size of order i i b˜ max ) largest batch size for order i that can be processed in i every allowable unit cFij/cPij ) fixed/proportional constant of the processing cost of order i in unit j estik/eftik ) earliest start/finish time for order i in stage k H ) a big number lbwjn/ubwjn ) lower/upper bound on time window for nth valid inequality of unit j limin/limax ) minimum/maximum number of batches for order i limax* ) largest possible number of batches for order i calculated by eq 2 lstik/lftik ) latest start/finish time for order i in stage k qi ) demand amount for order i qmax ) maximum demand for order i i ri/di ) release/due time for order i tailjn ) the sum of minimum processing times at stages later than stage Kj.
∑
∑
k′ < k j∈JAik′
min min (τ min ∑ ij ) ) lstik + min (τ ij ) j∈JA j∈JA k′ > k ik′
(A4)
ik
In addition, for each order in each stage a set of feasible predecessors Pi′k can be defined through eq A5 (Castro et al.14).
Pi′k ) {i∈I : estik e max (ubii′j)} ∀i′,k j∈JAik
(A5)
where ubii′j is the upper bound on the start time of order i in unit j assuming that it is followed by order i′,
ubii′j ) min (di τ min ij , di′ -
∑ min k′ > K j′∈JA
∑
min (τ ij′ )-
ik′
j
min min min (τ i′j′ ) - τ i′j - τ min ij ) ∀i,i′,j (A6)
k′ > Kj j′∈JAi′k′
and Kj is the stage associated with unit j. For makespan minimization, we use parameter tailjn on the RHS of tightening constraint (eq 28a).
tailjn ) min
∑
min min (τ ij′ )
i∈ICjn k′ > K j′∈JAik′ j
(A7)
Ind. Eng. Chem. Res., Vol. 47, No. 5, 2008 1555
Appendix B. Preprocessing Algorithms for Valid Inequalities. The preprocessing algorithm in Table B1 calculates the upper (ubwjn) and lower (lbwjn) bounds of time windows and tailjn by considering various subsets of orders ICjn for each unit j. The subsets are created by first removing the orders with the earliest starting time and then by removing those with the latest finish times. In this way, all possible combinations of time windows are covered. Table B1. Preprocessing Algorithm for Valid Inequalities for k∈{1,2,...|K|}, for j∈Jk, Initialize auxiliary set: IC′ ) I ∈IAj Reset cut number: n ) 0 Do while IC′ * L, Initialize active set: ICjn ) IC′ Do while ICjn * L and max (lftik) ) max (lftik), i∈ICjn
Update cut number: n ) n + 1 Calculate parameters: lbwjn ) min (estik)
i∈IC′jn
i∈ICjn
ubwjn ) max (lftik) i∈ICjn min tailjn ) min ∑ min (τij′ ) i∈ICjn k′>k j′∈JAik′
Update ICjn by removing the order(s) with the earliest start time: ICjn ) ICjn{i∈I: esti′k ) lbwjn} Update IC′ by removing the order(s) with the latest finish time: IC′ ) IC′\{i∈I: lftik ) ubwjn} Calculate the number of inequalities for unit j: Nj ) n
Acknowledgment The authors would like to acknowledge financial support from the American Chemical Society - The Petroleum Research Fund under Grant PRF #44050-G9, and Dr. Pradeep Prasad for his contributions to this work. Supporting Information Available: Tables of data for examples 3-5: process data; fixed and proportional constants of processing time; processing costs for orders; and demands, prices, release, and due times. This material is available free of charge via the Internet at http://pubs.acs.org. Literature Cited (1) Pinedo, M. Scheduling: Theory, Algorithms, and Systems; New York: Prentice Hall, 2002. (2) Pinto, J. M.; Grossmann, I. E. Assignment and Sequencing Models for the Scheduling of Process Systems. Ind. Eng. Chem. Res. 1998, 81, 433-466. (3) Pekny, J. F.; Reklaitis, G. V. Towards the Convergence of Theory and Practice: A Technology Guide for Scheduling/Planning Methodology. In Proceedings of the Third International Conference on Foundations of Computer-Aided Process Operations; Pekny, J. F.; Blau, G. E., Eds.: AIChE Symp. Ser. 1998, 94, 91-111.
(4) Shah, N. Single- and Multi-Site Planning and Scheduling: Current Status and Future Challenges. AIChE Symp. Ser. 1998, 94, 75. (5) 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-946. (6) Pinto, J. M.; Grossmann, I. E. A Continuous Time Mixed Integer Linear Programming Model for Short Term Scheduling of Multi-Stage Batch Plants. Ind. Eng. Chem. Res. 1995, 34, 3037-3051. (7) Pinto, J. M.; Grossmann, I. E. An Alternate MILP Model for ShortTerm Scheduling of Batch Plants with Preordering Constraints. Ind. Eng. Chem. Res. 1996, 35, 338-342. (8) Hui, C. H.; Gupta, A. A Novel MILP Formulation for Short-Term Scheduling of Multistage Multiproduct Batch Plants. Comput. Chem. Eng. 2000, 24, 1611-1617. (9) Hui, C. H.; Gupta, A.; Meulen, H. A. J. A Novel MILP Formulation for Short-Term Scheduling of Multistage Multiproduct Batch Plants with Sequence-Dependent Constraints. Comput. Chem. Eng. 2000, 24, 27052717. (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-711. (11) Gupta, S.; Karimi, I. A. Scheduling a Two-Stage Multiproduct Process with Limited Product Shelf Life in Intermediate Storage. Ind. Eng. Chem. Res. 2003, 42, 490-508. (12) Gupta, S.; Karimi, I. A. An Improved MILP Formulation for Scheduling Multi-Product, Multi-Stage Batch Plants. Ind. Eng. Chem. Res. 2003, 42, 2365-2380. (13) Castro, P. M.; Grossmann, I. E. New Continuous-Time MILP Model for the Short-Term Scheduling of Multi-Stage Batch Plants. Ind. Eng. Chem. Res. 2005, 44, 9175-9190. (14) Castro, P. M.; Grossmann, I. E.; Novais, A. Q. Two New Continuous-Time Models for the Scheduling of Multi-Stage Batch Plants with Sequence Dependent Changeovers. Ind. Eng. Chem. Res. 2006, 45, 6210-6226. (15) Liu, Y.; Karimi, I. A. Scheduling Multistage, Multiproduct Batch Plants with Nonidentical Parallel Units and Unlimited Intermediate Storage. Chem. Eng. Sci. 2007, 62, 1549-1566. (16) Pantelides, C. C. Unified Frameworks for the Optimal Process Planning and Scheduling. In Proceedings on the Second Conference on Foundations of Computer Aided Operations, Cache Publications: New York, 1994; pp 253-274. (17) Harjunkoski, I.; Grossmann, I. E. Decomposition Techniques for Multi-Stage Scheduling Problems Using Mixed-Integer and Constraint Programming Methods. Comput. Chem. Eng. 2002, 25, 1533-1552. (18) Maravelias, C. T. A Decomposition Framework for the Scheduling of Single- and Multi-Stage Processes. Comput. Chem. Eng. 2006, 30, 407420. (19) Kallrath, J. Planning and Scheduling in the Process Industry. OR Spectrum 2002, 24, 219-250. (20) Prasad, P.; Maravelias, C. T.; Kelly, J. D. Optimization of Aluminum Smelter Casthouse Operations. Ind. Eng. Chem. Res. 2006, 45, 7603-7617. (21) Prasad, P.; Maravelias, C. T. Batch Selection, Assignment and Sequencing in Multi-Stage Multi-Product Processes. Comput. Chem. Eng. 2008, in press (DOI: 10.1016/j.compchemeng.2007.06.012).
ReceiVed for reView July 11, 2007 ReVised manuscript receiVed November 23, 2007 Accepted December 4, 2007 IE070944Y
