4902
Ind. Eng. Chem. Res. 2001, 40, 4902-4911
A Novel Nonuniform Discrete Time Formulation for Short-Term Scheduling of Batch and Continuous Processes Kyu-Hwang Lee, Heung Il Park, and In-Beum Lee* Department of Chemical Engineering, Pohang University of Science and Technology, San 31 Hyoja-Dong, Pohang, 790-784, Korea
We propose a novel nonuniform time discretization method (NUDM) formulation for a shortterm scheduling problem of batch and continuous operations. It aims to reduce the size of resultant MILP models by avoiding the use of a high dimensioned task-time matching binary variable defined in previous NUDM approaches. We devised several efficient constraints, with subdivided binary variables assigned to the start, process, and end event of each task. Several types of examples show that the proposed formulation is both general and efficient. For instance, we discuss fixed and variable processing time problems, and a sequence-dependent setup-time problem. A simple optimization algorithm tailored for a periodic scheduling problem is also proposed. 1. Introduction In chemical processing industries, it is essential to enhance process flexibility and production efficiency due to changing market demands and a variety of customer requirements. One common approach concerns shortterm scheduling. Recently, to increase profitability of overall plant operation, scheduling has been enhanced from simply shortening completion times to minimizing the loss of potential sales and possible waste of resources, such as manpower, materials, and so on. Scheduling is now being reinvestigated as a process of organizing, choosing, and timing resource usage, to carry out the activities necessary to produce desired outputs at desired times, while satisfying a large number of time and relationship constraints among activities and resources (Morton and Pentico1). It gives rise to a number of research works on the resourceconstrained short-term scheduling model for a batch and/or continuous process and the optimization algorithm for its solution. Kondili et al.2 proposed a state-task network (STN) representation, by which the resource-constrained shortterm scheduling problem is formulated as a mixed integer linear programming (MILP) model. It is based on a uniform time discretization method (UDM), which discretizes the time horizon into a number of intervals of uniform duration and assumes that the processing time is fixed. Although this formulation can clearly represent many complexities that appear in practical situations, it has some difficulties due to time discretization. That is, the equal time slot may require a significant large number of equations when the greatest common divisor of processing times is considerably smaller than the longest processing time. An approximation of processing times in order to avoid an increase in model size may cause an infeasible or suboptimal solution. To circumvent these drawbacks, an alternative approach leads to a nonuniform time discretization method (NUDM), which divides the time horizon into several nonuniform time slots only when used for establishing material balances. As a result of * To whom all correspondence should be addressed. Tel: +82-54-279-2274. Fax: +82-54-279-2699. E-mail: iblee@ postech.ac.kr.
avoiding the use of redundant binary variables, the NUDM formulation may generate smaller sizes and may need shorter computational times than the conventional UDM formulation. Zhang and Sargent3 proposed a general formulation based on NUDM and considered the processing time of a task as not only fixed but also variable time that is a linear function of batch size. Their formulation includes a binary variable that assigns a starting event for each task, and another binary variable to determine an ending event of each task. Because of the nonconvexity in a mixed integer nonlinear programming (MINLP) formulation, a simplified algorithm is proposed. Schilling and Pantelides4 proposed a general mathematical formulation with tighter timing constraints than previous formulations, such as that of Zhang and Sargent.3 As continuous time formulations cause large MILP problems with large integrality gaps, they proposed a special branch-and-bound algorithm, where branching takes place on both continuous and binary variables. They extended their formulation to optimal periodic scheduling problems.5 Castro et al.6 presented an improved formulation, a relaxation of the problem presented by Schilling.7 By allowing finite storage within the processing tasks resource equipments of the involved raw materials and/or products, their approach leads to less degenerate mathematical models. Mockus and Reklaitis8 used NUDM incorporating a continuous time variable, a distinct binary variable that represents a birth or death of unit at each time interval, and a three-dimensional binary variable that connect a start event to an end event. They formulated a shortterm scheduling problem as MINLP and proposed a global optimization algorithm (GOA).10 However, they reported that its performance is problem-dependent according to their computational results. By introducing the assumption that a task operates only once during a predetermined time interval, they modified their previous work for reducing of model size (Mockus and Reklaitis9). The authors compared their NUDM with UDM from computational results with several examples. As the comparison results between the UDM and the NUDM approaches are problem-dependent,
10.1021/ie000513e CCC: $20.00 © 2001 American Chemical Society Published on Web 09/25/2001
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001 4903
they did not tell exactly which was better. However, they showed that the NUDM approaches could be improved as discretization lengths shorten. These kinds of NUDM formulations used high dimensional binary variables that, for example, were defined with a dimension of time multiplied by time in order determine the ending of a task for feasible material balances. These cause a large size of a model and may be difficult to apply to practical problems. To solve practical problems, Ierapetritou and Floudas11-13 proposed a simplified model. The authors proposed a binary variable dimension reduction method, which was made by separating the unit and task assign variables. In addition, their formulation can treat variable processing time as a linear function of batch size. With some sequencing constraints generated from recipes, authors showed an efficient solution in many cases. The large numbers of different classes of constraints to describe resource utilization were called sequence constraints by Ierapetritou and Floudas.11 Castro et al.6 reported that their formulation could be less accurate, and their solutions were infeasible in some cases because of these constraints. Although NUDMs provide a chance to reduce model size, the model efficiency may be inferior to that of UDM in some cases. The reason is that supplementary binary variables with a high dimension may be required to assign the task event and the time axis in the material balances. The objective of this paper is to develop a novel NUDM that addresses the short-term scheduling problem of batch and continuous processes. A MINLP formulation is presented by virtue of the introduction of new binary variables and the effective employment in task-unit-time relationships. They correspond to subdivided events, such as the start, process, and end of a task. In section 2, several important features of a problem considered in this study are described. Basic assumptions and definitions of new binary variables are presented. A novel mathematical formulation is constructed in section 3. Several examples from the literature include a sequence-dependent setup time problem and a periodic scheduling problem. They are considered and their computational results are illustrated in section 4. By changing the number of event points, the computational results are shown in section 5. Finally, in section 6, the performance of the proposed formulation is discussed and concluded. 2. Problem Statement and Characteristics The problem under consideration is a short-term scheduling of batch and/or continuous processes represented by STN. Given (i) process recipes, (ii) available units and their capacity limits, (iii) available storage limits for each material, and (iv) the time horizon under consideration, the objective is to determine the optimal sequence of task and operation size of each task taking place in each unit within restricted resources. Several important features of the problem are as follows. (i) Nonuniform Time Discretization. The time horizon is divided into Nmax time grids. The interval duration will be treated as a variable and determined optimally in optimization calculation. The total number of time grids is determined by an iterative procedure, which was used in the literature.4,6,11 (ii) Variable Processing Time. To handle various operation cases, the processing time of each task is
described as a linear function of operating size. For example, in the case of a batch operation with fixed duration, a variable term of the processing time equals zero. In continuous operation, a variable term is set as an inverse of a processing rate and a constant term is a minimum run length.5 (iii) Sequence-Dependent Setup Time and Cost. Equipment may require cleaning or some extra preparation between operations. If all cleaning times and costs of different units are assumed to be the same or small enough to be neglected, they can be included within a processing time and cost term, respectively. However, they will greatly change by the production sequence in reality. For instance, in a dye manufacturing plant, a unit may require extensive cleaning time or cost if a dark dye is to be processed before a lighter one. By contrast, little or no cleaning may be required in the reverse situation.9 In this case, a sequencedependent setup time must be considered, and can be classified into a fixed setup time and a variable setup time that depends on previous operation size. (iv) Mixed (Batch and Continuous) Operation. In addition to batch operations, continuous or mixed operation can also be modeled with the proposed formulation. The key difference between batch and continuous operation can be found in representation of a processing time and a change of the inventory level. To take it into account, a separate quantity of each state is devised. In this study, the following assumptions are made, which are almost same as those given by Ierapetritou and Floudas.11 (i) The material is transferred instantaneously from states to tasks and from tasks to states. (ii) No preemptive operation is allowed in a batch operation. (iii) All data are deterministic and fixed over the time horizon of interest. For UDM, the duration of a time interval is predetermined and the processing time is fixed, so the end point of a task can be easily obtained. Therefore, a task can be assigned to a time interval by using a binary variable. However, in NUDM, a variable time interval should be determined. Furthermore, when a processing time is employed as a function of the operation size, it is difficult to match several events of a task to a time axis. Moreover, to identify an end time, it has been indispensable for previous NUDMs to introduce another high-dimensioned binary variable, which leads to an increase in the total model size. In this study, as shown below, a task is partitioned into three separate events, such as start, process and end. For each event, a binary variable is defined as follows.
WSi,j,n ) 1 if task i starts in unit j at the time grid n 0 otherwise WPi,j,n ) 1 if task i is in processing in unit j at the time grid n 0 otherwise WEi,j,n ) 1 if task i ends in unit j between the time grid n - 1 and n 0 otherwise
{
{ {
The detailed definition of these binary variables can be
4904
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001
Figure 1. Definition of binary variables.
explained in Figure 1, where a Gantt chart of task i in unit j starting at time grid n + 1 and ending up at n + 3 is illustrated. The variable WSi,j,n+1, which corresponds to the starting event at time grid n + 1, is equal to 1, while the others are equal to 0. In the same fashion, only one variable WEi,j,n+3 has the value of one at time grid n + 3. It should be noted that the variable WPi,j,n+2 equals 1 because a task is being processed, but the time interval does not contain the starting or ending event. These subdivided binary variables can explicitly describe all events in a task without the additional use of a task-time matching binary variable; therefore, a more effective treatment of the task-unit-time relationship can be accomplished with a reduced problem size.
subject to
material balances of each state binary variable correlations operation size correlations nonuniform time discretization correlations sequence-dependent setup times.
Material Balances. A quantity of state s at time grid n is determined by the amounts of batch/continuous operation, a demand and a delivery of state s at the same time grid like eq 1. b c STs,n ) STs,n + STs,n - Ds,n + Rs,n, ∀ s ∈ S,n ∈ N (1)
Equation 2 represents the initial inventory correlations separated from batch and continuous operations. b c + STs,0 ∀s∈S STs,0 ) STs,0
b b STs,n ) STs,n-1 -
∑ ∑ Fcs,iBSi,j,n + ∑ ∑ Fps,iBEi,j,n,
b i∈Iin s j∈Ji
b i∈Iout s j∈Ji
(2)
As shown in the above equations, the proposed formulation considers the quantity of two parts in each state; one is a quantity from a batch operation and the other is from a continuous operation. The reason for variable separation can be found in the example of Figure 2. If the quantities at time grid n and n + 1 are expressed with the conventional state definition, they might be simply STs,n ) 3 and STs,n+1 ) 3. However, by classifying the quantities according to the sort of operations, they can be represented as STs,n ) 3, STs,n+1 ) 3, b c ) - 2, STcs,n ) 0, and STs,n+1 ) 2. STbs,n ) 3, STs,n+1 Although there is no change in the net inventory level at each grid, it can be found that not only a quantity of each state but also a quantity variation of each state between adjacent time grids can be described, i.e., the
∀ s ∈ S,n ∈ N (3)
∑ ∑ Fcs,iBSi,j,n-1 + ∑ ∑ Fps,iBEi,j,n,
c i∈Iin s j∈Ji
The mathematical model for a short-term scheduling of batch and/or continuous processes can be formulated, as follows:
profit
inventory level increase in continuous operation with the slope of 2 amounts per time interval. With a moderate increase in the variable number, a more precise description of the inventory level can be acquired.
c STcs,n ) STs,n-1 -
3. Mathematical Formulation
maximize
Figure 2. Inventory level example for material balance.
c i∈Iout s j∈Ji
∀ s ∈ S,n ∈ N (4)
Equations 3 and 4 represent material balances for batch and continuous operation, respectively. In each equation, the quantity of state s at time grid n is equal to that of time grid n - 1 updated by any quantity produced or consumed between the time grid n - 1 and n. In a batch operation, the production quantity is equal to an end quantity of input task into a state and a consumption quantity is equal to a start quantity of the output task from that state. However, a “batch” size of a continuous operation means an accumulated quantity between time grid n and n + 1. That is, if a continuous task starts at time grid n with a batch size BSi,j,n, the batch size does not mean a quantity used at time grid n but a preassigned quantity that will be used between time grid n and n + 1. Therefore, the consumption term is determined by the quantity of the batch size that was previously started, as denoted BSi,j,n-1. Binary Variable Correlations. The constraints among events in a task, such as start, process, and end, are described with the binary variables defined in the previous section. From the definition, three events of a task cannot happen at a time grid, simultaneously. Instead, only one task can start or remain in processing at a time grid (eq 5). Similarly, only one task can stay in processing or end (eq 6). When there is no task being performed in a time grid, no binary variable will be assigned, which is expressed by the inequalities of eq 5 and 6.
(WSi,j,n + WPi,j,n) e 1, ∑ i∈I
∀ j,n
(5)
(WPi,j,n + WEi,j,n) e 1, ∑ i∈I
∀ j,n
(6)
j
j
Equation 7 expresses the correlation of binary variables between two adjacent grids. That is, if a task starts or is in processing in a time grid, the possible event in the next grid is a processing or ending. It should be noted that if a task starts or is in processing event
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001 4905
at a time grid, the same task can start at the next grid, whenever it ends at the next grid.
WSi,j,n-1 + WPi,j,n-1 ) WPi,j,n + WEi,j,n, i ∈ Ij, ∀ j,n > 1 (7) Equations 8 and 9 express the fact that a task cannot be processed or end up at the initial time grid, respectively. With the same reason, eqs 10 and 11 are introduced for the final grid.
WPi,j,1 ) 0, i ∈ Ij, ∀ j
(8)
WEi,j,1 ) 0, i ∈ Ij, ∀ j
(9)
max WSi,j,n ) 0, i ∈ Ij, ∀ j
(10)
Figure 3. Time correlation example.
cleaning time, the setup time term will be neglected.
τj,n ) SEj,n +
(Ri,jWSi,j,n + βi,jBSi,j,n), ∑ i∈I
∀ j,n (17)
j
max ) 0, i ∈ Ij, ∀ j< WPi,j,n
(11)
In the case of a continuous task, eq 12 represents the prevention of the continuous task in processing across a time interval. That is, at each grid, continuous tasks should start or end in order to calculate the material balance defined in eq 1.
WPi,j,n ) 0, i ∈ Icj ,∀ j,n
(12)
Operation Size Correlation. When the operation size at the end time of a task is to be calculated in previous NUDM approaches, the end point should be deduced through high-dimensioned binary variables that represent interrelationships between the task and the order of the time grid. To avoid the use of the largesized binary variables, we directly define an operation size variable for each event of a task in the same way as was done for the task-unit-time binary variables, i.e., BSi,j,n, BPi,j,n, and BEi,j,n, which corresponds to WSi,j,n, WPi,j,n, and WEi,j,n, respectively. The following equation constrains the start operation size and the end size in material balance, which is the same correlation of eq 7. That is, if a task was previously started or processed with some quantity, the task should be processed or completed with the same quantity at the next grid.
BSi,j,n-1 + BPi,j,n-1 ) BPi,j,n + BEi,j,n, i ∈ Ij, ∀ j,n > 1 (13) The following three equations represent the upper and lower limits of each operation size, where the quantity of material undergoing a task i in unit j should be bounded by its maximum and minimum capacity, respectively. If a task is not performed (that is, the binary variable is zero), its operation size is set to zero. min Bi,j WSi,j,n
e BSi,j,n e
max Bi,j WSi,j,n,
i ∈ Ij, ∀ j, ∀ n (14)
min max Bi,j WPi,j,n e BPi,j,n e Bi,j WPi,j,n, i ∈ Ij, ∀ j, ∀ n (15) min max WEi,j,n e BEi,j,n e Bi,j WEi,j,n, i ∈ Ij, ∀ j, ∀ n Bi,j (16)
The processing time in each unit corresponds to the sum of a setup time and an operating time; that is, a linear function of operating size. If there is no need for a
Demand Correlation for Intermediate Due. Because of contractual obligations, it is often necessary to deliver a certain amount of Ds,n to customers at multiple Tn’s. After subdividing a problem into several subproblems on the bases of due dates, the time when demands arise should be within the end time at the previous time grid and a task’s start time. Although the number of event points between intermediate due dates should be determined by optimization procedures, there is no systematic procedure to determine it. Therefore, to determine the number of event points, iterative procedures are used.
TSj,n g Tn, i ∈ Ij, ∀ j,n
(18)
TEj,n-1 e Tn, i ∈ Ij, ∀ j,n > 1
(19)
Nonuniform Time Discretization Correlations. With previous constraints, the start, process, or end event of each task and its operation size can be determined. To combine them into a time axis via variable time intervals, the following constraints on the duration of each time interval are presented. Equation 20 means that the total sum of interval durations should be smaller than time horizon concerned. Equation 21 shows the precedence of end time of each unit. In eq 22, a start time of each unit is defined as the sum of previous interval durations. Equation 23 means that the end time of each unit should be larger than the sum of its starting time and processing time.
∑
∆Tn e H
(20)
n 1 TSj,n )
∑ ∆Tn′,
∀ j,n < nmax
(21) (22)
n′ 1 (24)
j
To reduce a searching space in the optimization step, several logical relations between time relevant variables are included as constraints. According to our computational studies, the CPU times and the number of nodes calculated without the following constraints are multiplied by more than two times. Equation 25 expresses that the remaining time duration at time grid n should be equal to or greater than the sum of processing times of tasks that will be assigned after a time grid n. When a task which has started at time grid n and ends at time grid n + 1 (that is, no processing event), the time duration between a time grid n and n + 1 should be equal to or greater than the processing time of the task. Otherwise, it will be relaxed.
H-
∑ ∆Tn′ g ∑ n′en n < n′
1 (29) Each unit at each time grid can be utilized by only one task, as shown in eq 30.
WTi,j,n e 1, ∑ i∈I
∀ j,n
(30)
j
∀ j,n < nmax - 1 (26)
j
Sequence-Dependent Setup Time Correlations. We assumed that a setup time includes not only preparation time for a starting task but also the cleaning time of the used unit in the previous operation. As mentioned previously, there can be two kinds of sequence-dependent setup times; a fixed setup time and a variable setup time. While the former is based on the relative orders of processing tasks, the latter is dependent not only on the relative orders but also on the previous operation size. In previous researches, the relative sequence is determined with the help of highly dimensioned binary variables, which are defined for all relative orders in each unit at each time grid. To describe the relative order, we must identify the task, which is recently carried out in a given unit from the reference time grid of interest, which may be called a “recently performed task (RPT)”. When there is no idle time between successive tasks, the task finished is defined as an RPT. However, if there is an idle time, we must deduce the RPT from the time precedence relation among previous performed tasks. To accomplish this clearly and efficiently, a binary variable is introduced as follows, which will require a lower dimension compared to one in previous works.
WTi,j,n )
{
1 if task i is a RPT in unit j at time grid n 0 otherwise
With this variable, the problem of an idle time can be simply treated by making the value of WTi,j,n be the same as that of WTi,j,n-1 as the time grid proceeds. When an end event occurs at a time grid, the corresponding
In a fixed setup time case, the setup time can be obtained from eq 31. That is, if task i is an RPT and task i′ starts in a unit, then the setup time is the same as the given setup time between production i and i′ in unit j, Si,i′,j. It can be rewritten in a linearized form, as in eqs 32 and 33.
SEj,n ) Si,i′,j where WTi,j,n ) 1 and WSi′,j,n ) 1, i ∈ Ij, i′ ∈ Ij, ∀ j,n (31) SEj,n e Si,i′,j + U(2 - WTi,j,n - WSi′,j,n), i ∈ Ij,i′ ∈ Ij, ∀ j,n (32) SEj,n g Si,i′,j - U(2 - WTi,j,n - WSi′,j,n), i ∈ Ij,i′ ∈ Ij, ∀ j,n (33) Otherwise, in a variable setup time case, we need to know the previously operated size. With the same concept of WTi,j,n, a new continuous variable BTi,j,n is defined for it. BTi,j,n is a batch size of an RPT. If a task ends at time grid n, BTi,j,n is set to the end size of the task like eq 34. Equation 35 means that when there is a task operated at time grid n-1 and its event does not end, BTi,j,n is the same value of BTi,j,n-1.
BTi,j,n ) BEi,j,n where WTi,j,n ) 1 and WEi,j,n ) 1, i ∈ Ij, ∀ j,n (34) BTi,j,n ) BTi,j,n-1 where WTi,j,n ) 1 and WEi,j,n ) 0, i ∈ Ij, ∀ j, n > 1 (35) Equations 34 and 35 can be linearized like eqs 36 and 37. Because of eq 27, it is impossible to obtain a case of WTi,j,n ) 0 and WEi,j,n ) 1. Therefore, the term of 1 -
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001 4907
WTi,j,n - WEi,j,n is sufficient in eq 37.
-U(2 - WTi,j,n - WEi,j,n) e BTi,j,n - BEi,j,n e U(2 WTi,j,n - WEi,j,n), i ∈ Ij, ∀ j,n (36) - U(1 - WTi,j,n - WEi,j,n) e BTi,j,n - BTi,j,n-1 e U(1 - WTi,j,n - WEi,j,n), i ∈ Ij, ∀ j,n > 1 (37) Representing the cleaning cost as a linear function of operation size, the following equation can be written in the case of WTi,j,n ) 1 and WSi′,j,n ) 1.
-U(2 - WTi,j,n - WSi′,j,n) e CCLi,j,n - CCAi,i′,j CCBi,i′,jBTi′,j,n e U(2 - WTi,j,n WSi′,j,n), i ∈ Ij,i′ ∈ Ij, ∀ j,n (38) To consider the cleaning time, Mockus and Reklaitis9 have introduced a four-dimensioned binary variable whose number increases with a square function of the number of tasks. However, with the above variables and constraints, the number of binary variables is a linear function of the number of tasks. Objective Function. The objective is to maximize a profit that can be expressed as follows:
profit ) value of products - cost of feedstocks - cost of utilities and manpowers - cost of cleaning value of products ≡
max ) ∑s cs(∑n Ds,n + STs,n
(39)
This expression includes both the value of material left in storage at the end of the time horizon and the value of material delivered to customers during the horizon. If material of state s is an intermediate and does not need to be stored, it will be subjected to the negative setting of the value of that storage (cs) or the restriction constraint at the quantity of state s.
cost of feedstocks ≡
∑s cs(STs,0 + ∑n Rs,n)
(40)
The cost of feedstocks is obtained from both the initial inventory quantities and the input quantities that are given as a parameter during the horizon.
cost of utilities and manpowers ≡ u (CAi,j
∑n
∑j ∑ i∈I j
u WSi,j,n + CBi,j
∑n BSi,j,n)
(41)
The cost of utilities and manpowers is determined by the operation size. As in previous studies, it is assumed as a linear function of operation size. The cost of cleaning is defined by the sum of cleaning cost, which is calculated in eq 38
cost of cleaning ≡
CCLi,j,n ∑n ∑j ∑ i∈I
(42)
j
4. Computational Results To illustrate the relative merits of the proposed formulation, we considered several examples. Using
Figure 4. State-task network for examples 1 and 2. Table 1. Data for Example 1 and 2 (a) Tasks And Units processing time for example 1 (mean processing time for example 2)
unit
capacity
suitability
heater reactor 1 reactor 2 still
100 50 80 200
heating reaction 1,2,3 reaction 1,2,3 separation
1.0 2.0, 2.0, 1.0 2.0, 2.0, 1.0 P1:1 for prod 2, 2 for IntAB P2:1.95 for prod 2, 2 for IntAB
(b) States state
storage capacity
initial amount
price
feed A feed B feed C hot A intAB intBC impureE prod1 prod2
unlimited unlimited unlimited 100 200 150 200 unlimited unlimited
unlimited unlimited unlimited 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 10.0 10.0
Table 2. Comparison of Computational Results of Example 1 method P1 UDM NUDM proposed P2 UDM NUDM proposed
no. ofvariables no. of no. of CPU binary continuous constraints nodes time (s) 86 465 165 1550 642 165
218 1351 401 3962 1846 401
226 2588 1011 4358 3571 1064
96 146 177 2101
0.27 12.5 5.4 >1 h
6945
96
AMPL/CPLEX 7.0.0 as a solver on a IBM-compatible, PENTIUM III 933 MHz machine running Windows ME, the solutions of the proposed mathematical formulation are obtained, except in example 2. Moreover, in example 2, a PENTIUM III 450 MHz machine is used in order to compare performances with the reported results of Castro et al. Example 1. This example is adopted from the paper of Kondili et al. In this example, two different products are produced through five processing stages: heating, reaction 1, 2, and 3, and separation. The basic recipe is presented in Table 1 and Figure 4. P1 indicates the problem where the processing time of the separation task for an intermediate AB is one, hence its greatest common divisor (GCD) of the processing times is 1.0. P2 is the modified problem by Mockus and Reklaitis, where the processing time of the separation task for intermediate AB is 1.95, so its GCD is 0.05. The time horizon under consideration is 10 h. The comparison of computational results with previous methods is summarized in Table 2. The proposed MILP formulation involves 165 binary variables, 401
4908
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001 Table 4. Data for Example 3 (a) States
Figure 5. Gantt chart of example 2.
states
capacity limits
feed int prod1 prod2
unlimited 5000 unlimited unlimited
Ri,j ) 200
(b) Costs βi,j ) 0.6
no. of event no. of binary variables no. of variables no. of constraints obj MILP obj relaxed LP nodes CPU time (s)
Castro et al.
Ierapetritou and Floudas
5 80
6 89
5 40
209(177)
226
375
376
297
507
1704.23 1480.50 22(31) 0.66(1.0)
1804.36 1480.06 60 0.32
2191.11 1480.05 747 117
CM s ) 0.18
time
Schilling and Pantelides
5 56(88)
10 8
(c) Demands
Table 3. Comparison of Computational Results of Example 2 proposed approach
prices 5
states
4
prod1 prod2
200 50
150
260
units
size
units suitability
processing times
374
unit1 unit2 unit3
1500 1000 1000
task1 task2 task3
1 for P1, 0.9 for P2 1 for P1, 1.1 for P2 1 for P1, 0.8 for P2
1732.36 1503.15 51 0.28
continuous variables and 1011 constraints for P1 and 165 binary variables, and 401 continuous variables and 1064 constraints for P2. The number of searching nodes for each problem in AMPL/CPLEX required 177 and 6945, respectively. It can be noted that UDM shows a better efficiency than NUDM in P1, where the duration of time interval is relatively large. However, in P2, the efficiency of UDM is not so good because the number of grid points increases from 11 to 201, as the GCD decreases from 1 to 0.05. While NUDM of Mockus and Reklaitis8 failed in solving P2 problem, our formulation provides the optimal solution. As the solution of P2 problem by UDM cannot be obtained within reasonable times in our system, we used the number of nodes reported in the paper of Mockus and Recklitis.8 Example 2. This example is presented to consider the variable processing time case, where all data are the same as in example 1 except that the processing times vary within 33% around the mean values and the time horizon of interest is 8 h. It has been dealt with in works of Schilling and Pantelides,4 Ierapetritou and Floudas,11 and Castro et al.6 For comparison purposes, the results of each formulation are quoted from Castro et al.6 The number of variables and examined nodes are shown in Table 3. The computational result of the proposed formulation is illustrated with a Gantt chart in Figure 5. The proposed formulation requires 88 binary variables, 177 continuous variables, and 376 constraints. However, since the three kinds of binary variables (WSi,j,n,WPi,j,n, and WEi,j,n) are dependent on each other, WEi,j,n is set to a continuous variable that has a range of [0,1]. In doing so, the formulation involves 56 binary variables, 209 continuous variables and 376 constraints. Moreover, the number of nodes for optimal solution is reduced from 31 to 22. According to Castro et al.,6 although a better objective function value was reported by Ierapetriotu and Floudas,11 their solution is infeasible since it violates the time horizon under consideration. When the number of events is 5, the proposed formulation requires more time than that of Castro et al.,6 but slightly improved objective function value can be obtained.
6
7
10
11
300
400 200
100
12 100
(d) Units, Tasks
Example 3. The intermediate due date case is investigated with the example of Mockus and Reklaitis.8,9 All data are shown in Table 4. While the GCD of processing times in P1 is one, it equals 0.1 in P2. As shown in Table 5, UDM fails to obtain a solution due to the increase of the binary variable number from 36 to 335. By contrast, both NUDMs can obtain the optimal solution without increase of model sizes. Note, however, that the proposed formulation is more efficient than that of Mockus and Reklaitis8,9 in terms of the number of searching nodes and the size of the model. Example 4. A sequence-dependent setup-time problem is considered, using the example of Mockus and Reklaitis.9 The single-unit sequencing problem is that given amounts (75 kg each) of three products (white, gray and black dye) is produced within a fixed time period (6 h). Each product requires the execution of a single task with a processing time of 1 h. After utilization, the unit requires a different cleaning time (60, 90, 75, and 66 min, respectively) in preparation for the next task. As already mentioned in the paper of Mockus and Recklitis,9 Table 6 shows that the number of variables and constraints of the UDM formulation will dramatically grow up as the discretization length is reduced. However, the NUDM formulations including the proposed formulation remain unchanged in the size of models. In this example, there is no need to define a processing and ending of each task because production is performed in a single unit, and the proposed formulation produces a smaller sized mode than the previous ones. Example 5. An example of Schilling and Pantelides5 is examined in order to show that the proposed formulation can be applied to a periodic scheduling problem whose objective is to maximize a profit per cycle length. They proposed an MINLP formulation and tackled the complexity caused by the nonlinear objective function with multiple MILP reduction. To alleviate computational difficulties in solving the MINLP model, we devised a simple algorithm based on the proposed formulation. If the objective function can be assumed a convex function, we can formulate the problem in an MILP form by fixing the length of the cycle because its feasible value can be roughly predicted in practical situations.
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001 4909 Table 5. Comparison of Computational Results of Example 3 no. of variables P1 P2
method
binary
continuous
no. of constraints
UDM NUDM proposed UDM NUDM proposed
36(33) 218(168) 53 335(305) 218(168) 53
101(93) 364(99) 142 940(860) 364(99) 142
88 838 348 1173 838 348
A value in parentheses comes form Mockus and
no. of nodes
CPU time (s)
5270 731 4
0.98 1.3 0.071
4849/10836 6
0.15
Reklaitis5.
Table 6. Comparison of Computational Results of Example 4 no. of variables cleaning discretization no. of method time (min) length (min) binary continuous constraints UDM NUDMa proposed UDM NUDMa proposed UDM NUDMa proposed UDM NUDMa proposed a
60
60
90
30
75
15
66
6
60 36 21 109 36 21 209 36 21 509 36 21
95 45 16 174 45 16 334 45 16 814 45 16
289 57 49 777 57 49 2432 57 49 12 797 57 49
Figure 6. Gantt chart of example 5. Table 8. Comparison of Computational Results of Example 5 no. of variables Schilling and Pantelides5 proposed
binary
continuous
no. of constraints
objective function value
80 (avg)
429 (avg)
432 (avg)
28.72
56
128
313
28.77
NUDM means the formulation by Mockus and Recklitis.9
Table 7. Durations and Limits of Size of Example 5 reactor 1 reactor 2 still
reaction 1 reaction 2 reaction 3 reaction 1 reaction 2 reaction 3 separation
duration 13.333 + 0.27(size) 13.333 + 0.27(size) 6.667 + 0.133(size) 13.333 + 0.17(size) 13.333 + 0.17(size) 6.667 + 0.0833(size) 6.667 + 0.0333(size)
size 20-80 20-80 20-80 30-80 30-80 30-80 40-80
An initial inventory STs,0, which was considered as a parameter in a previously suggested formulation, is regarded as a variable. In addition, from the definition of cyclic operation, the amount of the final inventory is set equal to that of the initial inventory, as follows. max STs,0 ) STs,n , ∀s
(43)
If we redefine a demand quantity, which has been defined as a parameter in the previous formulation, as a variable that denotes a produced quantity of a final product during a cycle, a final product state also can satisfy the eq 43. On the basis of these facts, an algorithm to apply the proposed formulation to a periodic scheduling problem is as follows: (i) as the first step, assume the minimum and maximum values of cycle length; (ii) solve the MILP maximizing a profit within a fixed cycle length; (iii) find the cycle length with a golden section search method; (iv) repeat procedures ii and iii until the cycle length converges within a tolerance range. The STN of Example 5 is the same as Example 1 without a heating task. The duration of each task and its size of limit are shown in Table 7. With a reduced sized model, we can obtain a good solution using this simple approach. Table 7 shows that the formulation of Schilling and Pantelides gives the objective function value of 28.72 with 80 binary variables, but the better
value of 28.77 can be obtained with only 56 binary variables in the proposed formulation. The Gantt chart of the solution is shown in Figure 6. As Schilling and Pantelides proposed the hybrid branch-and-bound algorithm that can be parallelized, it is difficult to directly compare our formulation with theirs. However, from the viewpoint of the size of models, the proposed model and procedure can obtain the similar solution with smaller size of a model. 5. Number of Event Points According to the results of previous literature or several examples shown above, the size of model formulated based on UDM approaches is dramatically increased when the number of event points is increased. One reason is that unused time slots may be defined and require more time to obtain the optimal solution. The NUDM approaches were developed in order to overcome these disadvantages. However, as there is no systematic approach to determine the number of time slots, we follow many researchers in using an iterative procedure. To obtain a solution, the number of event points is increased until the objective function value does not change. When the number of event points is less than the optimal one, the obtained solution cannot guarantee its global optimality. As the number of events points increases, to obtain a solution becomes harder. This tradeoff is one of the major disadvantages of NUDM approaches. After applying iterative procedures to Example 2 that Castro et al.6 used for examining the performances by changing the number of events points, computational results of the proposed formulation are shown in Table 8. More than five time slots are required to obtain the optimal solution, since the objective function value does not change. Although the CPU time of Castro et al. is less than ours at five event points, total CPU times of their formulation to guarantee its optimality is at least 5.23 s and ours is 4.17 s. Moreover, as the number of event points grows up, the proposed formulation shows
4910
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001
Table 9. Computational Statistics for Several Event Points event points 3
4
5
6
7
methods
Pa
Cb
P
C
P
C
P
C
P
C
binary continuous constraints relaxed obj obj nodes CPU (s)
24 94 159 520 520 0 0.04
25 87 90 520 520 5 0.25
40 151 267 966 866 1 0.16
47 146 175 975 866 15 0.25
56 209 376 1704 1480 22 0.66
80 226 297 1804 1480 60 0.32
72 267 485 1829 1480 357 3.31
122 324 446 2296 1480 1666 4.41
88 325 504 1829 1480 2203 15.5
173 440 622 2775 1480 41 024 143
a P means the proposed formulation. b C means the formulation of Castro et al.
better performances than one of Castro et al.6 Moreover, the number of binary variables, the relaxed objective function value and computational time of their formulation are more increased than ones of the proposed formulation. These results can be expected, because their formualtion uses high dimensional binary variable defined by time multiplied by time, to represent the processing of a task. In the proposed formulation, as the dimensionality of binary variables is reduced, much better performances are expected when the number of event points increases. 6. Conclusions We proposed that a novel nonuniform time discretization formulation for short-term scheduling of batch and continuous processes. In previous NUDM formulations, one or two kinds of high dimensional binary variables are defined to assign events to the relevant time interval. In this paper, we introduce a new set of binary variables that interrelates a task, a unit, and a time grid. On the basis of this, a general MINLP formulation is proposed and its linearized formulation requires smaller number of variables comparing previous ones. Several practical aspects in the scheduling problem can be taken into account such as a batch and continuous operation, a sequence-dependent setup time and an intermediate due date. Additionally, it is shown that the proposed formulation can be applicable to a periodic scheduling problem with a sequential search based approach. Several examples are considered to illustrate the efficiency of the proposed formulation. According to our computational studies, when the number of event points increase, the proposed formulation shows better performance than the previous formulations. Nomenclature Indices i ) task j ) unit n ) time grid s ) state Sets Ij ) a set of tasks which can be performed in unit j Ijb ) a set of batch tasks which can be performed in unit j Ijc ) a set of continuous tasks which can be performed in unit j Iin s ) a set of tasks producing material of state s Iout ) a set of tasks receiving material from state s s Ji ) a set of units in which a task i can be performed
Jib ) a set of units in which a batch task i can be performed Jci ) a set of units in which a continuous task i can be performed N ) a set of time grids Parameters Ri,j ) a constant term of processing time of task i at unit j βi,j ) a variable term of processing time of task i at unit j Fcs,i ) an input proportion of task i from state s Fps,i ) an output proportion of task i to state s max Bi,j ) a maximum capacity of unit j when used for performing task i min Bi,j ) a minimum capacity of unit j when used for performing task i Cs ) a price factor of state s u CAi,j ) a constant term of operating cost of task i at unit j u CBi,j ) a variable term of operating cost of task i at unit j CCAi,i′,j ) a constant term of cleaning cost of task i′ after task i in unit j CCBi,i′,j ) a variable term of cleaning cost of task i′ after task i in unit j Ds,n ) a demand quantity of state s at the beginning of time grid n H ) a time horizon of interest Rs,n ) an input quantity of resource of state s at time grid n Si,i′,j ) a setup time when task i′ is operated after task i in unit j Smax ) a maximum quantity of state s s Smin ) a minimum quantity of state s s STs,0 ) an initial quantity of state s Tn ) a time grid n corresponding to an intermediate due date Variables ∆Tn ) Duration of time interval between grid n and n + 1 τj,n ) a processing time of unit j at time grid n BEi,j,n ) an ending batch size of task i in unit j at time grid n BPi,j,n ) a batch size in processing of task i in unit j at time grid n BSi,j,n ) a starting batch size of task i in unit j at time grid n CCLi,j,n ) a cleaning cost of task i in unit j at time grid n SEj,n ) a setup time of unit j at time grid n b STs,0 ) an initial quantity of state s related with batch operations c STs,0 ) an initial quantity of state s related with continuous operations b STs,n ) a quantity of state s related with batch operations at time grid n c STs,n ) a quantity of state s related with continuous operations at time grid n STs,n ) a quantity of state s at time grid n, Ssmin e STs,n e Smax s TEj,n ) an end time of unit j at time grid n TSj,n ) a start time of unit j at time grid n WEi,j,n ) a binary variable that assigns the ending of task i in unit j between the time grid n - 1 and n WPi,j,n ) a binary variable that assigns the processing of task i in unit j at time grid n WSi,j,n ) a binary variable that assigns the starting of task i in unit j at time grid n
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001 4911 WTi,j,n ) a binary variable that assigns a RPT by task i in unit j at time grid n
Acknowledgment This work was supported by the international research program of the KOSEF (Grant No. 995-1100-0022) and the Brain Korea 21 project. Literature Cited (1) Morton, T. E.; Pentico, D. W. Heuristic Scheduling Systems: with Applications to Production Systems and Project Management; Wiley: New York, 1993. (2) Kondili, E.; Pantelides, C. C.; Sargent, R. W. H. A General Algorithm for Short-term Scheduling of Batch Operations-I. MILP Formulation. Comput. Chem. Eng. 1993, 17, 211. (3) Zhang, X.; Sargent, R. W. H. The Optimal Operation of Mixed Production Facilities-A General Formulation and Some Approaches for the Solution. Comput. Chem. Eng. 1995, 20, 897. (4) Schilling, G.; Pantelides, C. C. A Simple Continuous-time Process Scheduling Formulation and a Novel Solution Algorithm. Comput. Chem. Eng. 1996, 20, S1221. (5) Schilling, G.; Pantelides, C. C. Optimal Periodic Scheduling of Multipurpose Plants. Comput. Chem. Eng. 1999, 23, 635. (6) Castro, P.; Barbosa-Po´voa, A. P. F. D.; Matos, H. An Improved RTN Continuous-Time Formulation for the Short-Term
Scheduling of Multipurpose Batch Plants, Ind. Eng. Chem. Res. 2001, 40, 2059. (7) Schilling, G. Optimal Scheduling of Multipurpose Plants. Ph.D. Thesis, University of London, London, U.K., 1997. (8) Mockus, L.; Reklaitis, G. V. Mathematical Programming Formulation for Scheduling of Batch Operations Based on Nonuniform Time Discretization. Comput. Chem. Eng. 1997, 21, 1147. (9) Mockus, L.; Reklaitis, G. V. Continuous Time Representation Approach to Batch and Continuous Process Scheduling. 1. MINLP Formulation. Ind. Eng. Chem. Res. 1999, 38, 197. (10) Mockus, L.; Reklaitis, G. V. Continuous Time Representation Approach to Batch and Continuous Process Scheduling. 2. Computational Issues. Ind. Eng. Chem. Res. 1999, 38, 204. (11) Ierapetritou, M. G.; Floudas, C. A. Effective Continuoustime Formulation for Short-Rerm Scheduling. 1. Multipurpose Batch Processes. Ind. Eng. Chem. Res. 1998, 37, 4341. (12) Ierapetritou, M. G.; Floudas, C. A. Effective Continuoustime Formulation for Short-Term Scheduling. 2. Continuous and Semi-Continuous Processes. Ind. Eng. Chem. Res. 1998, 37, 4360. (13) Ierapetritou, M. G.; Hene´, T. S.; Floudas, C. A. Effective Continuous-Time Formulation for Short-Term Scheduling. 3. Multiple Intermediate Due Dates. Ind. Eng. Chem. Res. 1999, 38, 3446.
Received for review May 24, 2000 Accepted August 1, 2001 IE000513E
