Ind. Eng. Chem. Res. 2005, 44, 9129-9145
9129
Mixed-Time Representation for State-Task Network Models Christos T. Maravelias* Department of Chemical and Biological Engineering, University of WisconsinsMadison, Madison, Wisconsin 53706
A new mixed-time representation for STN-based scheduling models is proposed: the time grid is fixed, but processing times are allowed to be variable and span an unknown number of time periods. The proposed representation handles batch and continuous processes, and it accounts linearly for holding, backlog, and utility costs. It handles release and due dates at no additional computational cost, and it accounts for variable processing times. Computational improvements for the solution of the proposed mixed-integer programming formulation are also developed. 1. Introduction The purpose of the paper is to present a model that can be readily used for the simultaneous optimization of scheduling and supply chain management (SCM) problems. A number of issues that are usually neglected in standalone scheduling models become important when scheduling is integrated with SCM: • Smooth demand patterns at the customer-facing nodes are transformed into high peaks of demand at the production facilities, due to the inventory policies applied at the intermediate nodes (warehouses, distribution centers, etc.) to account for demand uncertainty. • Demand cannot always be satisfied on time. Backlog costs are introduced to account for unmet demand. The modeling of holding and backlog costs, therefore, is important. • Shipments are usually made at fixed time intervals (e.g., at the end of each week). This implies that scheduling models have to handle multiple fixed due dates at low computational cost. • Features such as seasonality of demand, max-ship life of final products, and scheduled maintenance and shut-downs require long planning horizons. In problems with long planning horizons, the tradeoff between changeover costs and holding costs (i.e., the lot-sizing problem) becomes very important. To account for scheduling decisions when we optimize the supply chain of process industries, therefore, current scheduling models should be amended to simultaneously handle lot-sizing, holding, backlog, changeover, and shipment costs and multiple intermediate due dates, while solution methods should be enhanced to effectively handle long time horizons. Several mixed-integer programming (MIP) scheduling models have been proposed in the literature (see refs 1-5 for reviews). The state task network (STN) formulation6 and the equivalent resource task network (RTN) formulation7 are two powerful modeling tools that have been refined by several researchers8-19 and extensively used for planning and scheduling problems. In discretetime STN/RTN models, the time horizon is divided into N periods (intervals) of equal duration, common for all units, and tasks begin and finish exactly at a time point. The disadvantage of discrete-time models is that variable processing times can be handled only as discrete * Corresponding author. Tel.: (608) 265-9026. Fax: (608) 262-5434. E-mail:
[email protected].
approximations and at high computational cost. In continuous-time models, the time horizon is divided into time periods of unequal and unknown duration. Continuous-time models account for variable processing times, but since the time points are not fixed, holding and backlog costs cannot be calculated linearly, due dates are very expensive to model, and the number of intervals needed to accurately represent the optimal solution is unknown. To solve problems where the objective is the minimization of total cost (i.e., holding + backlog + changeover costs) with intermediate due dates, therefore, a time representation that accounts for variable processing times, it linearly handles holding and backlog costs and accounts for due dates at no additional cost, is needed. Maravelias and Grossmann20 showed that the assumption of a fixed-time grid is not necessarily coupled with the requirement for constant processing times. In this paper, we exploit this observation and propose a mixed-time representation where the time grid is fixed but the processing time of the tasks can be variable. The proposed representation handles batch, continuous, and semicontinuous processes. It accounts linearly for holding, backlog, and utility costs, it handles due dates at no additional computational cost, and it accounts for variable processing times. In other words, it combines the modeling advantages of both discrete- and continuous-time models. The continuous-time STN model of Sundaramoorthy and Karimi16 for the scheduling of batch processes that is used for the derivation of the new time representation is presented in Section 2. The new MIP formulation for the scheduling of batch processes is presented in Section 3, while the modifications for the scheduling of continuous and semicontinuous processes are presented in Section 4. The modeling of special features, such as changeover times, shared storage tanks, orders with due dates, and holding and backlog costs, is presented in Section 5. A set of strong valid inequalities that lead to tighter formulations for the scheduling of batch processes are presented in Section 6. A reformulation for the scheduling of continuous processes is presented in Section 7. The application of the proposed approach is illustrated via three example problems in Section 8. 2. Continuous-Time STN Model Sundaramoorthy and Karimi16 proposed a continuoustime STN MIP model for the scheduling of batch
10.1021/ie0500117 CCC: $30.25 © 2005 American Chemical Society Published on Web 11/01/2005
9130
Ind. Eng. Chem. Res., Vol. 44, No. 24, 2005
processes that relies on the novel idea of time balances. Their model is the first not to include a big-M type of constraints, leading to a tighter formulation. A modified version of this model is used for the derivation of the mixed-time representation. Note, however, that any continuous-time STN model could have been used instead. A common, continuous partition of the time horizon is used to account for all possible plant configurations and resource constraints: N + 1 time points {0, 1, 2, ..., N} are defined, and the time horizon is partitioned into N periods {1, 2, ..., N}, where period n starts at time point n - 1 and finishes at time point n. A task idle_j is defined to represent idle time on equipment unit j, and task decoupling is used. Assignment constraints are expressed through binary Zjn for unit j and binaries Wsin, Wpin and Wfin for task i. Binary Zjn ) 1 if a task starts on unit j at time point n. Binary Wsin ) 1 if task i starts at time point n, binary Wpin ) 1 if task i is being processed at time point n, and binary Wfin ) 1 if task i finishes at or before time point n. The duration of task i that starts at time point n is denoted by Din, while the time remaining for the completion of the task that is being processed at unit j at time point n is denoted by δjn. The batch size of task i that starts at, is being processed at, and finishes at or before time point n is denoted by Bsin, Bpin, and Bfin, respectively. The inventory level of state s at time point n is denoted by Ssn, and the amount of resource r consumed by various tasks at time point n by Rrn. Time point n corresponds to time Tn, while the duration of period n (i.e., the period that starts at Tn-1 and finishes at Tn) is denoted by ∆tn. 2.1. Assignment Constraints. Since idle time is modeled as an additional task, a unit j is always allocated to a task. This implies that, whenever a task finishes, another task starts (eqs 1 and 2). Equation 3 ensures that no task can start/finish when another task is being processed, and eq 4 is a balance on the status of task i. The constraints in eqs 5 and 6 hold at the beginning and end, respectively, of the scheduling horizon.
Zjn )
∑ Wsin
∀j, ∀n < N
(1)
∑ Wfin
∀j, ∀n > 0
(2)
∀j, ∀n
(3)
i∈I(j)
Zjn )
i∈I(j)
∑
Wpin + Zjn ) 1
storage capacity for state s, and I(s)/O(s) is the set of tasks that consume/produce state s
BMIN Wsin e Bsin e BMAX Wsin i i
∀i, ∀n
(7)
BMIN Wpin e Bpin e BMAX Wpin i i
∀i, ∀n
(8)
(4)
Wfi0 ) Wpi0 ) 0
∀i
(5)
WfiN ) WpiN ) 0
∀i
(6)
2.2. Batch-Size Constraints and Material Balances. Equations 7, 8, and 9 are used to impose lower and upper bounds on the batch sizes Bsin, Bpin, and Bfin, respectively, while eq 10 is the batch-size balance that enforces variables Bsin and Bfin to be equal for the same batch of task i. Equation 11 is the mass balance equation and the capacity constraint for state s at time point n, where Fis is the stoichiometric coefficient of state is the s in task i (negative for consumed states), SMAX s
∑ FisBfin + ∑ FisBsin e
Ssn ) Ss,n-1 +
i∈O(s)
i∈I(s)
∀s, ∀n (11)
SMAX s
2.3. Resource Constraints. The amount RrnI of renewable resource r required by task i that starts at time point n is calculated by eq 12, where γir and µir are the fixed and variable requirements, respectively, of task i for resource r. The same amount RrnO is “released” when task i finishes and is calculated in eq 13. The total amount RBrn of resource r required at time n by batch processes is calculated and bounded not to by eq 14: exceed the maximum availability RMAX r
RIirn ) γirWsin + µirBsin
∀i, ∀r, ∀n
(12)
RO irn ) γirWfin + µirBfin
∀i, ∀r, ∀n
(13)
B RBrn ) Rrn-1
∑i ROirn + ∑i RIirn e RMAX r
∀r, ∀n (14)
2.4. Timing Constraints. The duration, Din, of task i that starts at time point n (i.e., at time Tn) is calculated via eq 15, where Ri and βi are the fixed and variable terms, respectively. The time remaining for completion, δjn, of the task being processed on unit j during period n is bounded by the constraints in eqs 16 and 17.
Din ) RiWsin + βiBsin δjn gδjn-1 +
∀i, ∀n > 0
(9)
∀i, ∀n > 0 (10)
Bsin-1 + Bpin-1 ) Bpin + Bfin
∑ Din-1 - ∆tn
∀i, ∀n ∀j, ∀n > 0
(15) (16)
i ∈I(j)
i∈I(j)
Wsin-1 + Wpin-1 ) Wpin + Wfin
∀i, ∀n
Wfin e Bfin e BMAX Wfin BMIN i i
δjn e
∑ (aiWpin + βiBpin)
∀j, ∀n > 0 (17)
i ∈I(j)
Equations 18-20 define the start and the end of the time horizon and enforce an ordering among time points, where H is the fixed scheduling horizon:
T0 ) 0
(18)
TN ) H
(19)
∆tn ) Tn - Tn-1
∀n>0
(20)
2.5. Objective Function. Various objective functions can be accommodated; here, we consider the maximiza-
Ind. Eng. Chem. Res., Vol. 44, No. 24, 2005 9131
Figure 1. Alternative time representations for STN/RTN-based models.
Figure 2. Zero-wait storage policy in the proposed mixed-time representation.
tion of revenues, where ζs is the price of product s, and sales occur only at the end of the horizon:
max Z )
∑s ζsSsN
and 24, respectively.
Zjn, Wsin∈{0,1}, Wpin, Wfin∈[0,1],
(21)
B Bsin, Bpin, Bfin, Ssn, RIirn, RO irn , Rrn , Din, δjn g 0 (25)
Finally, variables Wpin and Wfin need not be defined as binaries. The continuous-time MIP model (M1) consists of the constraints in eqs 1-22.
Model (M2) may yield slightly suboptimal solutions because tasks are not forced to finish exactly at a time point and equipment units may remain idle at the end of some time intervals. Nevertheless, model (M2) exhibits a number of modeling advantages: • Compared to continuous-time models, it accounts linearly for holding, backlog, and utility costs and handles release and due dates at no additional computational cost. • Compared to discrete-time models, it does not require the introduction of additional binary variables for the modeling of tasks with variable processing times. The proposed model can be readily extended to continuous and semicontinuous processes, with very good computational results. The only modeling limitation of model (M2) is that it cannot always handle zero-wait storage policy between two batch processes, because the end time of a task is not required to coincide with a point of the time grid (Figure 2). The gap between batches A and B in Figure 2a implies that the zero-wait requirement is not satisfied. Zero-wait policy, however, is usually applied between two continuous processes, where the output of the first process is continuously consumed by the second process. This case is effectively handled: tasks A and B are processed simultaneously and the zero-wait requirement is satisfied (Figure 2b). In terms of computational effectiveness, it seems that model (M2) “carries the computational curses” of both discrete and continuous representations, namely, the large number of time intervals and the poor linear programming (LP) relaxation gaps, respectively. Nevertheless, as will be shown in Section 6, the fact that the time grid is fixed allows us to apply a number of computational improvements for the solution of (M2).
Zjn, Wsin∈{0,1}, Wpin, Wfin∈[0,1], B Bsin, Bpin, Bfin, Ssn, RIirn, RO irn, Rrn, Tn, ∆tn, Din, δjn g 0 (22) 3. Mixed-Time Representation Maravelias and Grossmann20 showed that discretetime models can be derived from continuous-time models if the following two restrictions are applied: (A) Fixed-time grid: the scheduling horizon is divided into time intervals of equal and fixed (known) duration; (B) Fixed processing times: the processing times of all tasks are constant multiples of the duration of the time interval ∆t of the fixed-time grid. Their analysis showed that the assumption of fixedtime grid is not necessarily coupled with the assumption of fixed processing times. In this paper, we employ a fixed-time grid but we allow processing tasks to have variable duration, spanning an unknown number of time periods. Note that the opposite (variable-time grid, constant processing times) is a special case of continuous-time models. A task is forced to start exactly at a time point but is allowed to finish anywhere within a time interval. A schematic of the different time representations is shown in Figure 1. The restriction of the fixed-time grid means that the time Tn at which time point n occurs, and the length ∆tn of period n, are fixed parameters, defined via eqs 23 and 24, and not variables subject to the constraints in eqs 18-20:
∆tn ) ∆t ) H/N Tn ) n(H/N) ) n∆t
∀n ∀n
(23) (24)
The proposed mixed-time MIP model (M2), therefore, results directly from model (M1) and consists of eqs 1-17, 21, and 25, where ∆tn and Tn are given by eqs 23
4. Modeling of Continuous Processes The proposed model can be readily extended to handle continuous processes. In continuous processes, input/ output states are continuously consumed/produced, which means that the mass balance equation should
9132
Ind. Eng. Chem. Res., Vol. 44, No. 24, 2005
account for material produced/consumed throughout the execution of a task and not just at the beginning and the end of a task. To model this feature, we define a new variable Bin which is the amount of task i being processed during period n. Note that Bin is equal not to the batch size of task i but rather to the “per period production” of task i; e.g., if continuous task i starts at time point n, its production rate is 40 kg/h, its processing time is 2 h, and ∆t ) 0.5 h, we will have
Bin′ ) 40 (kg/h) × 0.5 (h) ) 20 (kg), n < n′ e n + 4 If task i was a batch process, we would have
Bsin ) Bpin+1 ) Bpin+2 ) Bpin+3 ) Bfin+4 ) 2 (h) × 40 (kg/h) ) 80 (kg) While the assignment constraints in eqs 1-6 remain exactly the same, the batch-size constraint for continuous processes is written as follows,
BMIN Wpin + Wfin e Bin e BMAX (Wpin + Wfin) i i ∀i, ∀n (26)
RCirn ) γir(Wpin + Wfin) + µirBin
where the resource requirement may be overestimated during the last period, enforcing feasibility. The resource constraint for plants with both batch and continuous processes is given by eq 33, where RBrn is given by eq 14:
Rrn )
∑ RCirn + RBrn e RMAX r
(27)
) rMAX ∆t BMAX i i
(28)
∀r, ∀n
(33)
i∈IC
In continuous processes, the duration Din of task i does not depend on the per period production Bin. This means that we do not need to define variables Din and δin because the batch size and the duration of a task are decoupled. Restrictions on the minimum and maximum production length (run time) can be, instead, and enforced through variables Wsin and Wfin. If DMIN i are the minimum and maximum production DMAX i length, respectively, of task i, we can directly calculate ) and maximum (NMAX ) number of the minimum (NMIN i i intervals in which task i can be carried out as follows:
and BMAX are the where is a small number and BMIN i i minimum and maximum amounts, respectively, of task i that can be processed in ∆t units of time:
) rMIN ∆t BMIN i i
∀ i, ∀ r, ∀ n > 0 (32)
) DMIN /∆t NMIN i i
(34)
) DMAX /∆t NMAX i i
(35)
The restrictions on the production length of continuous task i can be enforced via eqs 36 and 37. A task can start in time point n only if it finishes between points n and n + NMAX : + NMIN i i θeNMAX i
/rMAX is the minimum/maximum processing where rMIN i i rate of continuous task i. Note that we use a small lower bound on the production during the last time period, because task i may finish before time point n (i.e., it may finish anytime between Tn-1 and Tn). Furthermore, by enforcing ∆t] if Wfin ) 1, we exclude solutions where Bin∈[, rMAX i Wfin ) 1 but Bin ) 0. If it is sufficient to consider runs whose length is an exact multiple of ∆t, eq 26 is replaced by eq 29:
Similarly, a task can finish in period n only if it starts between points n - NMAX and n - NMIN : i i
(Wpin + Wfin) e Bin e BMAX (Wpin + Wfin) BMIN i i ∀i, ∀n (29)
The MIP model (M3) for the scheduling of continuous processes consists of eqs 1-6, 21, 26 or 29, 30, 32, 33, 36, 37, and 38.
The mass balance equation for continuous processes is written as follows:
Zjn, Wsin∈{0,1}, Wpin, Wfin∈[0,1],
Ssn ) Ss,n-1 +
∑
FisBin e
SMAX s
∀s, ∀n (30)
i∈O(s)∪I(s)
Note that, if a plant consists of batch and continuous processes, the mass balance equation reads
Ssn ) Ss,n-1 +
∑
FisBfin +
i∈IB∩O(s)
∑
∑
i∈IB∩I(s) FisBin e SMAX s
FisBsin + ∀s, ∀n (31)
i∈IC∩(O(s)∪I(s))
where IB and IC are the sets of batch and continuous processes, respectively. The amount RCirn of resource r consumed by continuous task i during period n depends on the amount Bin processed during period n and not the total batch size,
Wsin e
∑
Wfin+θ
∀i, ∀n
(36)
θgNiMIN
θeNMAX i
Wfin e
∑
Wsin-θ
∀i, ∀n
(37)
θgNMIN i
Bin, Ssn, Rrn, RCirn g 0 (38) Batch processes where feeds/products are added/removed at different (but fixed) time points can be modeled as described in Kondili et al.6 Semicontinuous processes with no restrictions on their processing time can be modeled either as continuous process (if the width ∆t is approximately equal to the processing time) or as batch processes with fixed duration. 5. Extensions 5.1. Changeover Times. In discrete-time STN models,6 changeovers can be modeled as separate tasks whose duration is a multiple of the duration ∆t of the period of the time grid. The modeling of changeovers as separate tasks means that an entire period is assigned to only one changeover task even if the actual
Ind. Eng. Chem. Res., Vol. 44, No. 24, 2005 9133
changeover time is smaller than ∆t. This may lead to suboptimal solutions if the duration of changeovers is substantially smaller than ∆t. To overcome this limitation, we consider the changeover from task i′ to task i as part of task i, and we take it into account when calculating the processing time Din of batch task i∈IB and the per period production Bin of continuous task i∈IC. The variable Xi′in is used to model the changeover from task i′ to task i* i′ at time point n, where tasks i and i′ are nonidle tasks and Xiin ) 0, ∀i, ∀n:
Xi′in g Wfi′n + Wsin - 1 ∀j, ∀i ∈I(j), ∀i′∈I(j), i * i′, ∀n > 0 (39)
∑
Xi′in eWfi′n
∀j, ∀i′∈I(j), ∀n > 0
(40)
∑
Xi′in eWsin
∀j, ∀i∈I(j), ∀n > 0
(41)
i∈I(j)\{i′}
new binary variable Yjsn that is equal to 1 if state s is stored in tank j during period n. If JT is the set of shared is the capacity of storage tank storage tanks and SMAX j j∈JT, the constraints in eqs 44 and 45 are added:
∑ Yjsn e 1
Din )
∑
σi′iXi′in + (RiWsin + βiBsin)
i′∈I(j)\{i}
∀i, ∀n (42)
where σi′i is the changeover time from task i′ to task i. If changeover times are smaller than ∆t, no adjustments are needed for the scheduling of batch processes. For continuous processes, we need only adjust the production during the first period of operation:
∑i′ σi′iXi′in-1) e Bin e (∆t Wpin - ∑σi′iXi′in-1) rMAX i i′
rMIN (∆tWpin i
∀i, ∀n (43)
Changeovers larger than ∆t can also be handled, but additional modifications are needed. In the proposed approach, a time period is not allocated to a changeover task whose duration is smaller than ∆t, which in turn means that equipment units are better utilized. The disadvantage, however, is that the states used as inputs in a batch task i should be available before the changeover to task i begins, and not before task i begins. In general, it is expected that the proposed approach will be better when changeover times are small compared to ∆t. Note that the two representations can be used simultaneously. Another interesting feature is that we can model changeovers from and to idle, because idle time is modeled as an additional task. This issue is often neglected, although it is very common in the process industry (e.g., tasks where heating to operating conditions are necessary). Note that, if changeovers from and to idle are important, the duration of continuous tasks should be exactly equal to a multiple of ∆t, i.e., the constraint in eq 29 should be used to bound the per period production of continuous task i. 5.2. Shared Storage Tanks. Shared storage tanks can be handled as equipment units. For each state s that can be stored to shared tank j (i.e., s∈S(j)), we define a
(44)
∀j∈JT, ∀s∈S(j), ∀n
Yjsn Ssn e SMAX j
(45)
Equation 44 does not allow more than one state to be stored in tank j∈JT at any time, and the constraint in eq 45 bounds the inventory of state s at time point n. If state s can be stored in more than one tank, eq 45 is expressed for variable Ssjn, which denotes the amount of state s stored in tank j during period n, and constraint 46 is added,
Ssn )
∑
∀s∈S*, ∀n
Ssjn
(46)
j∈JT(s)
i′∈I(j)\{i}
Variable Xi′in always gets integer values, due to eqs 3941, and it is defined as continuous variable Xii′n∈[0,1]. The duration of batch task i∈IB is now given by eq 42 instead of eq 15,
∀j∈JT, ∀n
s∈S(j)
where S* is the set of states for which there is no dedicated storage and JT(s) is the set of storage tanks in which state s can be stored. 5.3. Orders/Deliveries with Intermediate Due/ Release Times. Let K be the set of all orders and deliveries, and K(s) the set of orders and deliveries of state s. The due/release time of order/delivery k∈K is TDk, and the amount due/delivered is ADk (ADk is negative if the amount is due and positive if it is delivered). The net amount of state s due/delivered at time point n is denoted by Psn. Since the time grid is fixed, we can preassign orders/deliveries to time points. If the due/release time cannot be assigned exactly to a time point (e.g., the due time is 3.75 h and we do not want to use a time grid with ∆t ) 0.25 h), we assign it to the previous/next time point to ensure feasibility; i.e., if Tn e TDk < Tn+1, then order k∈K(s) is assigned to time point n, and Psn is equal to ADk. If there are multiple orders and/or deliveries for state s at time point n, Psn is equal to the sum of the corresponding AD k. The only modification needed to account for orders/deliveries with fixed due/release dates, therefore, is in the mass balance equation:
Ssn ) Ss,n-1 +
∑
∑
FisBfin +
i∈(O(s)∩IB)
FisBin + Psn e
∑
i∈(I(s)∩IB) SMAX s
FisBsin + ∀s, ∀n (47)
i∈(I(s)∪O(s)∩IC)
Note that, in continuous-time models, it is not possible to preassign due/release times to time points because we do not know the number of periods that will be needed to represent the optimal solution. Even if the total number of periods was known, it would not be possible to a priori determine how many periods are located between two due dates and thus preassign a due date to a time point. To model due dates, therefore, we need to use disjunctive programming with additional binary variables, and the resulting continuous-time models are computationally inefficient. 5.4. Holding and Backlog Costs. In practical applications, orders cannot always be met on time, and methods for the minimization of lateness are required. This feature is indirectly modeled through the introduction of a penalty term for backlogged demand. Delays can also be observed in the delivery of raw materials
9134
Ind. Eng. Chem. Res., Vol. 44, No. 24, 2005
and intermediate products, especially when interconnected plants are simultaneously optimized and an order in one plant is a delivery in another. In such a case, however, we can penalize the delayed shipment only once as an unmet order. The net amount Psn of state s that is shipped to/from another node of the supply chain at time point n is the sum of deliveries and orders, Psn ) P+ sn - Psn
∀s, ∀n
(48)
where P+ sn and Psn are the amounts delivered and due, respectively, at time point n. In this paper, we focus on single plants, and deliveries are, therefore, handled as fixed inputs, i.e., the amount of raw material or intermediate s delivered at time point n is always equal to P+ sn. To model unmet demand, we introduce the variable SD sn, which is equal to the amount of state s actually shipped to another node of the supply chain. Hence, if we replace eq 48 into the mass balance eq 47 and use variable SD sn instead of Psn, we obtain the following general mass balance equation:
Ssn ) Ss,n-1 +
∑
∑
i∈(O(s)∩IB)
FisBfin +
∑
FisBsin +
i∈(I(s)∩IB)
FisBin + Psn+ - SD sn e Cs
∀s, ∀n (49)
i∈(I(s)∪O(s)∩IC) Note that, if backlogs are not allowed, i.e., SD sn ) Psn, eq 49 reduces to eq 47. The backlogged demand SBsn of state s at time n is calculated in eq 50,
B D + PSBsn ) Ss,n-1 sn - Ssn
∀s, ∀n
∀s (51)
is the upper bound on the shipments of where P-MAX s state s, and ND(s) is the subset of time points that shipments are allowed. B If CH s and Cs are the holding and backlog costs of state s, and Cii′C is the changeover cost from task i to task i′, the objective function for the minimization of cost is written as follows:
min Z )
6. Valid Inequalities for the Scheduling of Batch Processes The fact that the scheduling horizon is divided into intervals of fixed duration enables us to calculate bounds on the number of intervals a task can span, something that is not possible, in principle, with conand longest tinuous-time models. The shortest DMIN i DMAX processing time of task i can be calculated from i eq 15:
DMIN ) Ri + βiBMIN i i
∀i
(53)
) Ri + βiBMAX DMAX i i
∀i
(54)
Since the duration of a time period is known, we can and maximum NMAX calculate the minimum NMIN i i number of time intervals a task i can span via eqs 34 and NMAX are used next to and 35. Parameters NMIN i i tighten model (M2). 6.1. Addition of Tightening Assignment Constraints. A task that starts at time point n should finish and n + NMAX , and similarly a task between n + NMIN i i that finishes at time point n should start between n and n - NMIN : NMAX i i
(50)
Note that shipments are usually made only at specified times, which means that if demand is not met on time it is backlogged until the next shipment date. There is an additional constraint, therefore, on variables SD sn, D -MAX SD n∈ND(s) ∧ SD sn e Ps sn ) 0n∉N (s)
tions with backlogged demand that could be avoided, while large values yield large integrality gaps. Nevertheless, backlogs are very common in practice, and they should be modeled. In Section 8, we elaborate how the proposed models can be used for the solution of such problems.
∑s ∑n (CHs Ssn + CBs SBsn) + ∑j ∑ ∑ ∑n Cii′CXii′n
(52)
i∈I(j)i′∈I(j)
The minimization of cost for given orders with fixed due dates is a very common practical problem, especially when we simultaneously optimize production, transportation, and inventory levels. In general, the introduction of backlogs makes scheduling problems more difficult, because infeasible points become (expensive) feasible points. An originally tight scheduling problem with very few feasible solutions, for example, becomes a loosely constrained problem with a practically infinite number of feasible solutions. It is, therefore, very important to correctly choose backlog cost CBs . Small values for CBs can potentially result to optimal solu-
θeNMAX i
∑
Wsin e
Wfin+θ
∀i, ∀n
(36)
Wsin-θ
∀i, ∀n
(37)
θgNMIN i θeNMAX i
Wfin e
∑
θgNMIN i
For fixed processing times (i.e., NMAX ) NMIN ) τi ) Ri), i i eqs 32 and 33 reduce to
Wfin ) Wsin-τi
∀i, ∀n gτi
which holds true in discrete-time STN models with constant processing times. Equations 1-6 enforce integer feasibility, but yield LP-relaxed solutions where the sum (Wsin + Wsi′n+ν) for two tasks i and i′, that can be assigned to unit j, is . In any greater than 1 for a time point ν|ν < NMIN i integer feasible solution, however, a task i′ cannot start on unit j if another task i has started less than NMIN i periods earlier, i.e., we always have (Wsin + Wsi′n+ν) e . More generally, a task 1, for i∈I(j), i′∈I(j) and νn. This means that the constraint in eq 55 is NMIN i satisfied by any integer feasible solution (i.e., it is a valid inequality), and its addition gives a tighter LP-relaxation: n′en
∑ ∑
Wsin′ e 1
∀j, ∀n
(55)
i∈I(j)n′>n-NMIN i
Note that, if processing times are fixed, we have NMIN i
Ind. Eng. Chem. Res., Vol. 44, No. 24, 2005 9135
) τi and the constraint in eq 55 reduces to the assignment constraint of the discrete-time STN model of Shah et al.:8 n′en
∑ n′>n-τ ∑ Wsin′ e 1
i∈I(j)
∀j, ∀n
(56)
i
Another feature of the LP-relaxed solutions of model (M2) is that the sum of tasks that finish on unit j before time point n is larger than the sum of tasks that have . In an started on unit j before time point n - NMIN i integer solution, however, a task i∈I(j) can finish at n or earlier. This means only if it has started at n - NMIN i that, in any integer solution, eq 57 must be satisfied:
∑ Wfin′ e ∑
n′en
Wsin′′
∀i, ∀n
(57)
n′′en-NMIN i
6.2. Addition of Tightening Constraints for Din. Let us assume that a task i is performed exactly once in an integer feasible solution of model (M2). This means that variable Din will be greater than or equal to DMIN i for one n∈{0, 1, ..., N - NMIN }, or that task i will span i periods. Hence, there must be at least at least NMIN i 1) binary variables Wpin that are equal to 1. If (NMIN i task i is performed more than once, then there must be at least
(NMIN i
- 1)
∑n Wsin
binaries Wpin that are equal to 1. In LP-relaxed solutions of model (M2), however, variables Wsin often obtain small fractional values that lead to values of Din variables that are smaller than DMIN , and the sum of i Wpin variables is smaller than
∀j
(60)
The assumption that idle tasks last exactly one period does not exclude feasible solutions, because idle tasks can appear consecutively. Hence, it leads to tighter formulations without introducing new variables and constraints, and it allows us to fix variables Wpidle_j,n and Wfidle_j,n:
Wpidle_j,n ) 0,
Wfidle_j,n ) Wsidle_j,n-1 ∀j, ∀n > 0 (61)
In the formulation of Sundaramoorthy and Karimi,16 the constraint in eq 17 enforces the processing time remaining on unit j to be smaller than the duration of this task (since Bsin ) Bpin′ for time points n′>n, for which Wsin ) Wpin′ ) 1) if a task i∈I(j) is being processed on unit j. A task i∈I(j), however, can only be processed in period n if it has started before time point n, which means that the remaining processing time at time point n can be, at most, equal to the duration of the task minus ∆t:
δjn e
∑ (aiWpin + βiBpin) - ∆t ∑ Wpin
i∈I(j)
∀j, ∀n
i∈I(j)
(62)
When no tasks are being processed on unit j at time point n (i.e., Wpin ) 0 and Bpin ) 0, ∀i∈I(j)), eq 62 is redundant because both the LHS (left-hand side) and the RHS (right-hand side) are equal to 0. When a task i∈I(j) is being processed, however, the RHS of eq 62 is smaller than the RHS of eq 17. Hence, eq 62 is tighter than eq 17. Alternatively, the RHS of eq 17 can be g ∆t, ∀i∈I(j): tightened via eq 63, if DMIN i
δjn e
∑n Wsin
(NMIN - 1) i
∑ (aiWpin + βiBpin) - ∆t ∑ Wsin-1
i∈I(j)
i∈I(j)
∀j, ∀n > 0 (63)
The addition of the valid inequality in eq 58, therefore, excludes these LP-relaxed solutions:
- 1)∑Wsin ∑n Wpin g(NMIN i n
Didle_j,n ) ∆tWsidle_j,n
∀i
(58)
The constraint in eq 59 is another valid inequality that interconnects variable Din with variables Wsin and Wpin, excluding a number of LP-relaxed solutions and gives an upper bound on the sum of processing times Din:
∑ Din′e∆t(n′gn ∑ (Wsin′ + Wpin′+1))
n′gn
(59) ∀i, ∀n < |N| - NMIN i 6.3. Tightening of Time Balance Constraints. Sundaramoorthy and Karimi16 proposed the novel idea of balance constraints for the processing time δjn remaining on unit j at the end of time period n. Their idea is novel in that eqs 16 and 17 replace big-M timematching constraints that were used in previous continuous-time STN models,9-15 and leads to tighter MIP formulations. To have a tight upper bound when unit j is idle, we can further assume that an idle task always has a duration equal to ∆t (i.e., Ridle_j ) ∆t, ∀j):
The constraint in eq 63 expresses the condition that the remaining time on unit j should be less than the duration of the task i∈I(j) being processed on unit j during period n minus ∆t, if task i started at n - 1. If g ∆t started at n - 1, eq 63 a task i∈I(j) with DMIN i holds as equality. If an i∈I(j) started before n - 1 and it is processed at time point n, eq 63 holds as inequality. If no task is being processed on unit j at time n, then no task started on unit j at n - 1, and both variable δjn and the RHS of eq 57 are equal to 0. The constraint in eq 63, therefore, is a valid inequality that is tighter than the constraint in eq 17 and yields a tighter formulation. 6.4. Addition of New Time Balance Variables and Constraints. We can further tighten our model by using the time remaining for a specific task i∈I(j) carried out on unit j. Variable δ h in denotes the remaining processing time of task i at time point n, and it is bounded via eqs 64 and 65, which are similar to eqs 16 and 62, respectively:
h in-1 - ∆t δ h in gDin-1 + δ δ h in e (RiWpin + βiBpin) - ∆tWpin
∀i, ∀n > 0
(64)
∀i, ∀n > 0 (65)
Finally, we can add eq 66 or use it instead of eq 65:
9136
Ind. Eng. Chem. Res., Vol. 44, No. 24, 2005
δ h in e (RiWpin + βiBpin) - ∆tWsin-1
∀i, ∀n > 0 (66)
6.5. Computational Effect of Valid Inequalities. In Sections 6.1-6.4, we proposed a number of valid inequalities. The addition of any single inequality to model (M2) leads to a tighter formulation. The number of these constraints, however, is large, and if all of them are added to (M2), it becomes more difficult to solve the resulting LP. Moreover, some of these originally valid inequalities become redundant in the presence of the other valid inequalities. In this section we study which inequalities are more effective. Example Ex1a, which is based on the network of Figure A1 in Appendix A and is a modified example from Kondili et al.,6 is used for the analysis (data in Appendix A). The objective is the maximization of revenues from sales for a time horizon of 8 h. The duration of the time interval is ∆t ) 0.5 h. First, we study eight different models: • (M2.1): (M2) + eq 36 • (M2.2): (M2) + eq 37 • (M2.3): (M2) + eq 55 • (M2.4): (M2) + eq 57 • (M2.5): (M2) with eq 62 instead of eq 17 • (M2.6): (M2) with eq 63 instead of eq 17 • (M2.7): (M2) + eq 64 and eq 65 • (M2.8): (M2) + eq 64 and eq 66 where (M2) consists of eqs 1-17, 21, 25, 60, and 61. All models were implemented in GAMS 21.3, and solved with CPLEX 9.0 on a Pentium M at 1.7 GHz running Windows XP. Models are solved to 0.5% relative optimality, using the default GAMS/CPLEX algorithmic options. A resource limit of 3600 CPU seconds is also used. All models consist of 260 binary variables. Models (M2) and (M2.1)-(M2.6) have 1 110 continuous variables; models (M2.7) and (M2.8) have 1 246 continuous variables. The number of constraints and the solution statistics are given in Table 1. As shown, the addition of any of the valid inequalities yields tighter formulations, leading to fewer nodes and shorter solution times. Models (M2.1) and (M2.4) appear to be the tightest, while models (M2.1), (M2.2), and (M2.4)-(M2.6) are computationally more effective. To further study the effect of the constraints in eqs 36, 37, 57, 62, and 63, we solve the same example using six new models (M2.9)-(M2.14), all of which have 260 binary and 1 110 continuous variables. The number of constraints and the solution statistics are given in Table 2. • (M2.9): (M2) + eqs 36 and 37 • (M2.10): (M2) + eqs 36, 37, and 57 • (M2.11): (M2) + eqs 36, 37, and 62 instead of eq 17 • (M2.12): (M2) + eqs 36, 37, and 63 instead of eq 17 • (M2.13): (M2) + eqs 36, 37, and 62 and 63 instead of eq 17 • (M2.14): (M2) + eqs 57, and 62 and 63 instead of eq 17 • (M2.15): (M2) + eqs 62 and 63 instead of eq 17 • (M2.16): (M2) + eqs 36, 37, 57, and 62 and 63 instead of eq 17 Note that the addition of the tightening constraints closes the integrality gap by more than 50%, from 23.6% for model (M2) to 11.2% for model (M2.16). Model (M2.13) appears to be the most effective. Models (M2.11), (M2.13), and (M2.16) yield tight formulations, while models (M2.10), (M2.12), (M2.13), and (M2.16) are effective in terms of computational requirements. There-
fore, we choose to use models (M2.10)-(M2.13) and (M2.16) to solve a series of examples: • Ex1b: Same process network/data as Ex1a; H ) 8 h; ∆t ) 0.33 h. • Ex1c: Same process network/data as Ex1a; H ) 12 h; ∆t ) 0.5 h. • Ex2a: Modified example from Maravelias and Grossmann21 (Example 2 in Appendix A); H ) 18 h; ∆t ) 0.5 h. • Ex2b: Same process network/data as Ex2a; H ) 18 h; ∆t ) 0.33 h. • Ex2c: Same process network as Ex2a; H ) 18; ∆t ) 0.25 h. The integrality gap (%), the number of nodes required, and the computational requirements for models (M2), (M2.10)-(M2.13), and (M2.16) for all instances are given in Tables 3, 4, and 5, respectively. Models (M2.11), (M2.13), and (M2.16) yield the smallest integrality gap, with the latter being the tightest. Models (M2.13) and (M2.16) usually require the fewest nodes. To fully study the effectiveness of the proposed inequalities, more examples need to be solved. For the purpose of this paper, however, we will use model (M2.16), which consists of eqs 1-16, 21, 25, 32, 33, 52, and 55-57. We will refer to this MIP model as model (M2*). 6.6. More Valid Inequalities. Despite the improvements described above, the scheduling of batch processes remains a very hard problem, and further computational improvements are needed in order to be able to solve problems of practical interest. Due to the fixed-time grid, the proposed time representation enables us to develop a large number of strong valid inequalities that can potentially enhance the solution of model (M2). Some of these inequalities are the following, n′en+NMAX i
Din e
∑
Tn′Wfin′ - TnWsin
n′gn+NMIN i
∀i, ∀n < N - NMIN (67) i θeNMAX i
Wsin e
∑
Wfi,n+θ -
(
θgNMIN i
)
ν7 200.0a >7 200.0a
(M2.10) (M2.11) (M2.12) (M2.13) (M2.16) 3.8 41.2 3 446.6 18.5 596.5 5 544.8
4.9 42.9 2 130.5 16.4 645.8 2 506.0
2.9 57.1 1 771.9 27.7 660.3 1 078.2
2.6 53.5 1 501.2 17.7 1641.9 6 478.6
3.8 43.9 1 871.4 20.7 573.9 2 377.0
Optimality gap after 7 200 CPU s: Ex1c ) 5.4%; Ex2b ) 1.8%; Ex2c ) 2.1%. a
and the time needed to solve each LP increases. These constraints, therefore, can be better used within a branch-and-cut scheme, where only violated inequalities are added. We are currently developing such a branchand-cut scheme that will be presented in another paper. 7. Reformulation for the Scheduling of Continuous Processes If we assume that the duration of a continuous process is a multiple of ∆t, and that there are no restrictions on the processing time of tasks, we can reformulate model (M3) using only one type of binary variable. More specifically, we can express all constraints using only binary variable Win, which is 1 if continuous task i is processed during period n and 0 otherwise. Note that, in continuous processes, variable Win refers to a time period and not a point. This implies that, if Wsin ) 1 at time point n, variable Win ) 0 for period n (period between points n - 1 and n) but Win )
∑ Win ) 1
∀j, ∀n
(70)
i∈I(j)
and the batch-size constraint is written as follows,
Win e Bin eBMAX Win BMIN i i
∀i, ∀n
(71)
Note that a lower bound on the production over the last because production lengths period is now equal to BMIN i are assumed to be multiples of ∆t. The mass balance equation remains the same (eq 30), and the resource constraint is written as follows:
Rrn )
∑i (γirWin + µirBin) e RMAX r
∀r, ∀n (72)
The reformulated MIP model (M3*) for the scheduling of continuous processes consists of the constraints in eqs 21, 30, and 70-73.
Win∈{0,1}, Bin, Ssn, Rrn g 0
(73)
Note that changeover constraints can be included in (M3*):
Xii′n g Win-1 + Wi′n - 1 ∀j, ∀i∈I(j), ∀i′∈I(j), i′ * i, ∀n > 0 (74)
∑ Xii′n eWin-1
∀j, ∀i∈I(j), ∀n > 0
(75)
i′∈I(j)
∑ Xii′n eWi′n
∀j, ∀i∈I(j), ∀n > 0
(76)
i∈I(j)
8. Examples and Computational Results Example 1 is used in Section 8.1 for qualitative comparisons with continuous-time representations. The application of the proposed model for the scheduling of continuous processes with holding, backlog, and changeover costs is illustrated in Sections 8.2 and 8.3. 8.1. Batch Processes. The maximization of production for the STN of Example 1 is considered. Tasks with both constant and variable processing times are included: tasks H and R1 have fixed processing times,
9138
Ind. Eng. Chem. Res., Vol. 44, No. 24, 2005
Table 6. Model and Solution Statistics for Maravelias and Grossmann (M&G)15 and (M2*) for Ex1a M&G15 # bin. var’s # cont. var’s # const. LP-relaxation objective CPU s nodes
(M2*)
N)6
∆t ) 0.5
∆t ) 0.33
∆t ) 0.25
112 766 1 667 2 370.6 1 890.4 4.5 1 406
260 1 110 2 383 1 964.0 1 766.7 3.6 150
388 1 630 3 519 1 948.0 1 834.2 4.8 81
516 2 150 4 655 2 017.1 1 837.0 43.9 1 099
while tasks R2, R3, and S have variable processing times, varying from 75% to 125% of the nominal processing time (used in Kondili et al.6). The proposed model is compared with the continuous-time STN model of Maravelias and Grossmann15 (M&G). The model and solution statistics for the model of M&G15 and model (M2*) are given in Table 6. The globally optimal solution for this example is obtained by the continuous-time model with six time periods (seven time points), and it has an optimal objective value of $1 890.40. Note that, in order to verify the optimality of this solution, we had to also solve the problem with seven and eight periods, and the computational requirements for the solution of the two additional MIPs were significantly higher than the CPU time reported. In particular, a limitation of all continuous-time MIP models is that the minimum number of time points/periods/events required to represent the optimal solution is unknown. In practice, the number is determined through an iterative approach during which the number of periods is increased until there is no further improvement. This, however, has two drawbacks. First, there is no guarantee that the necessary number of periods is N if a solution with the same objective function is found when N + 1 periods are used. There are several examples where there is no improvement from N to N + 1, but there is an improvement from N to N + 2. Second, the time to prove optimality (i.e., to solve the MIPs with N + 1 and N + 2 periods) is significantly larger than the time to find the optimal solution. In literature, however, it is common practice to report the time required to find the optimal solution and not to prove optimality. In this paper, we follow this convention. As expected, model (M2*) yields suboptimal solutions, and its computational performance worsens as the time grid becomes finer. Nevertheless, a relatively good solution with an objective value of $1 766.70 is obtained using ∆t ) 0.5 h, and a very good solution with an objective value of $1 837.00 is obtained with ∆t ) 0.25 h. The Gantt chart of the optimal solution obtained using the continuous-time representation is shown in Figure 3, while the Gantt chart of the best solution
Figure 3. Optimal solution of Example 1a using a continuoustime STN model.
obtained using model (M2*) with ∆t ) 0.5 h is shown in Figure 4. Note that, in the solution obtained by the continuous-time model, reactors RI and RII are fully
Figure 4. Optimal solution of Example 1a using (M2*) with ∆t ) 0.5. Table 7. Solution Statistics of Models Maravelias and Grossmann (M&G)22 and (M2*) for Example 1d M&G19 LP-relaxation objective CPU s nodes
(M2*)
N)9
∆t ) 0.5
∆t ) 0 0.33
∆t ) 0 0.25
2 119.1 1 490.4 292.5 51 239
1 564.0 1 366.7 1.8 79
1 548.0 1 434.2 3.7 63
1 617.1 1437.0 43.8 1 190
utilized, except for a small period of time at the end of the horizon due to luck of intermediate BC. In the solution obtained using model (M2*), however, units are not always fully utilized: there is idle time in RI between the two consecutive batches of R2 and in RII between the two consecutive batches of R3. To illustrate how due dates are handled by different time representations, we solved two modified instances of Example 1. In Example 1d, we add four orders of 10 kg of product P1 with intermediate due dates, and in Example 1e, we add three orders of 20 kg of P1 with intermediate due dates (see Table A1.3 in Appendix A). The scheduling horizon in Examples 1d and 1e is 8 and 10 h, respectively. In both examples, we try to maximize the revenues from additional sales at the end of the horizon. Using the continuous-time STN model of Maravelias and Grossmann as modified to account for intermediate due dates,22 we obtain a solution with an objective value of $1 490.40 with nine time periods. The proposed model yields an optimal solution with an objective ranging between $1 366.70 with ∆t ) 0.5 h and $1 437.00 with ∆t ) 0.25 h. The solution statistics for both models are given in Table 7. Example 1d has been “fabricated” to yield an optimal solution that is practically the same as the solution of Example 1a. The Gantt chart of the optimal solution of the continuous-time model of M&G22 is shown in Figure 5 and is almost identical to the optimal solution of
Figure 5. Optimal solution of Example 1d using a continuoustime STN model.
Example 1a (Figure 3). Nevertheless, the computational requirements for the continuous-time model increase significantly (from 4.5 to 292.5 CPU seconds) because new time points are now required to accurately describe the decrease in the inventory level of product P1 at t ) 4, 5, 6, and 7. The computational requirements of model (M2*), on the contrary, remain practically the same because orders are preassigned to time points. The only difference is that we now use eq 47 instead of eq 11 to express the mass balance of product P1. Example 1e is different in that the amounts due are larger and the time horizon is longer. The best solution
Ind. Eng. Chem. Res., Vol. 44, No. 24, 2005 9139 Table 8. Solution Statistics of Model Maravelias and Grossmann (M&G)22 for Example 1e LP-relaxation objective CPU s nodes optimality gap (%) a
N)7
N)8
N)9
N ) 10
N ) 11
2 421.8 1 783.6 75.6 18 230
