Ind. Eng. Chem. Res. 1988,27, 1840-1848
1840
w = wall 1,2 = reactions in eq 1 and 2
Literature Cited Beek, J. Advances in Chemical Engineering; Academic: New York, 1962;Vol. 111, p 234. Craig, D. F.; Burklow, B. W. Ammonia Plant Saf. 1980,22, 131. Davies, J.; Lihou, D. A. Chem. Process Eng. 1971,52,71. De Deken, J. C.; Devos, E. F.; Froment, G. F. ACS Symp. Ser. 1982, 196,181-197. Ergun, S. Chem. Eng. Prog. 1952,48,89. Filla, M. Chem. Eng. Sci. 1984,39, 159. Grover, S. S.Hydrocarbon Process. 1970,49(4),109. Haldor Topsoe, Communication 698,The Institution of Gas Engineers, 31st Autumn Research Meeting, 1965. Hottel, H. C.; Sarofim, A. F. Radiatiue Transfer; McGraw-Hik New York, 1967. Hyman, M. H. Hydrocarbon Process. 1968,47(7), 131. Kern, D. Q. Process Heat Transfer; McGraw-Hill: Kogakusha, Tokyo, 1950.
McGreavy, C.; Newman, M. W. Paper presented at the Institution of Electrical Engineers Conference on the Industrial Applications of Dynamic Modeling, Durham, 1969. Murray, A. P.; Snyder, T. S. Ind. Eng. Chem. Process Des. Dew. 1985, 24,286. Murty, C.V. S. Ph.D. Thesis, Indian Institute of Technology, Madras, 1987. Roesler, F. C. Chem. Eng. Sci. 1967,22, 1325. Shen, Cai-Da; Yu, Zun-Hong Shang-hai Hua Kung Hsueh Yuan Hsueh Pao 1979, 1-2,51. Singh, C. P. P.; Saraf,D. N. Ind. Eng. Chem. Process Des. Dew. 1979, 18(1),1. Stehlik, P.; Sika, J.; Bebar, L. Chem. Prum. 1986a,36(61), 342. Stehlik, P.;Sika, J.; Bebar, L. Chem. Prum. 1986b, 36(61), 454. Stehlik, P.; Sika, J.; Bebar, L. Chem. Prum. 1986c,36(61),512. Subramaniam, T.K. Hydrocarbon Process. 1967,46(9),169. Yu, Zun-Hong; Shen, Cai-Da; Pan, Hui-Qin; Cai, Guo-Qiang; Sun, Xing-Yuan Huagong Xuebao 1980,2,143.
Received for review November 30,1987 Revised manuscript received May 18, 1988 Accepted June 6,1988
Scheduling in Serial Multiproduct Batch Processes with Finite Interstage Storage: A Mixed Integer Linear Program Formulation Hong-ming K u and I f t e k h a r A. Karimi* Department of Chemical Engineering, Northwestern University, Evanston, Illinois 60208
Production scheduling of N products in an M-unit serial multiproduct batch process under the minimum makespan criterion is studied. A finite number of storage units may be present between consecutive processing units. Recurrence relations requiring O(NM) simple steps are presented for computing operation times of all products in a given production sequence. An optimal MILP formulation is developed for scheduling a set of batches with stable intermediates and is then extended to allow for unstable intermediates. A heuristic strategy is also proposed to reduce the number of binary variables in the above formulation by assigning products to specific neighborhoods of a sequence. This approximate formulation yields solutions within 1%of the optima when applied to a set of test problems. For years, continuous operations have been the most prevalent and sought-after mode in chemical processing. In recent years, however, there has been a renewed interest in batch processes for processing situations involved in the production of complex, high-value chemicals or the production of multiple products by sharing the use of process equipment. An important problem that arises in such operations is the scheduling of products to make efficient use of a multiprocessor production facility with limited resources. A general short-term scheduling problem in multiproduct batch processes consists of scheduling a set of product batches on a network of processors so as to optimize a suitable cost or performance measure. Normally, this problem requires specifications of a matrix of batch processing times for various operations, matrices of resource utilization requirements for each operation, a matrix of changeover times and costs for pairs of products for each operation, and a set of rules governing the transfer of batches between processing stages. Then scheduling involves the determination of a sequence in which the batches should be processed and of a corresponding time table for all operations to optimize a cost or system performance criterion. Because of the staged nature of process equipment in batch chemical plants, interstage storage is a vital com-
* Author t o whom correspondence should be addressed.
ponent of batch processes. Due to the unsteady-state nature of batch operations, the equipment utilization and the productivity are usually low in these processes. Storage installed between processing stages can help reduce idle times in these stages by freeing them to process other batches and thus increase the equipment utilization and the productivity of a multiproduct batch process. Intermediate storage can also be used to decouple periodic operation of adjacent batch or semicontinuous units, to minimize the effects of variation in process parameters, to moderate the effeds of equipment failures, and to isolate intermediates associated with different products. In addition, it may be useful to have interstage storage for various other purposes. Hence, various forms of storage are used in batch plants, and its absencelpresence significantly affects the scheduling of operations. The rules governing the transfer of batches between processing stages can be classified into four storage policies, depending upon the number of batches that can be stored in the interstage storage. These policies are (1)unlimited intermediate storage (UIS), (2) finite intermediate storage (FIS), (3) no intermediate storage (NIS), and (4)zero wait or no wait (ZW or NW). The interstage storage capacity is measured in terms of the number of units and not the physical size of storage since it is usually assumed that each storage unit can temporarily hold any product batch. In the UIS and FIS
0888-5885/88/ 2627-1840$01.50/0 0 1988 American Chemical Society
Ind. Eng. Chem. Res., Vol. 27, No. 10, 1988 1841
0 UNIT1
0
0 UNIT M
UNIT?
STORAGE 1
Z1 U N I T S
STORAGE 2
Z2 U N I T S
STORAGE
ZM-,
M-l
UNITS
Figure 1. FIS system configuration.
policies, an infinite (LN-1) and finite number of storage units is available between stages, respectively. In both the NIS and ZW policies, there is no storage between stages. While a batch, after its completion in a processing unit, may be held temporarily in the unit in the UIS, FIS, or NIS mode, it must be transferred immediately to the downstream unit in the ZW mode. A process employing a combination of the four storage policies at different stages is referred to as a MIS (mixed intermediate storage) system. Despite the significance of interstage storage in batch plants, virtually all of the existing work in scheduling focuses on multiproduct processes with the UIS policy (see Graves (1981) and Reklaitis (1982)). To date, there have been only two reports on exact algorithms for solving the multiproduct FIS problem optimally (Dutta and Cunningham, 1975; Wiede and Reklaitis, 1984a,b;Wiede et al., 1984). They also proposed a number of heuristic algorithms. In both studies, the objective was to minimize the total production time of all products or makespan. In the first report, Dutta and Cunningham (1975) proposed an exact algorithm based on dynamic programming to solve the two-unit FIS problem; because of its computational demand, two approximate algorithms were also presented. However, no computational comparisons were reported for any of these methods. In the second report, Wiede et al. (1984) presented a heuristic method for the two-unit FIS system. In an extension of this paper, Wiede and Reklaitis (1984a) considered the serial multiunit FIS system in which the available storage is not localized between two adjacent consecutive units but is shared by all processing units. In their subsequent work (Wiede and Reklaitis, 1984b), they considered the same system under the MIS policy. In either case, a branch and bound solution procedure was presented to solve the problem optimally, subject to a suboptimal completion time algorithm. In this paper, we investigate the problem of scheduling N products across an M-stage serial processing system with a single unit per stage. Interstage storage may or may not be present between any two consecutive units. Thus, this system operates under a combination of the UIS, FIS, and NIS policies. The difference in the way storage is handled between this work and that of Wiede and Reklaitis (1984a,b) is that in our work, storage units, if present, can only receive finished products from the upstream batch unit, while in the latter case the storage can receive finished products from any batch unit.
Problem Description An M-unit serial multiproduct batch plant consists of M batch units in series. As shown in Figure 1, zj (zj I0, = 1,M-1) storage units are present between batch units J and (j+l).In a multiproduct plant, all products pass through the series of M units; however, some units may be skipped for some products. A set of N product batches is to be produced, and the sequence of operations to be performed is known for each batch. The processing con-
ditions are known a priori for all the products, so the processing time, tij,of product i on unit j is specified. If a product skips a unit, its processing time will be specified as zero on that unit. The objective is to obtain the order in which the batches should be produced so as to maximize the productivity of the plant, Le., to minimize the total time required to produce all the batches or the makespan. We make the following assumptions regarding the multiproduct process: 1. All products are produced in the same order on each processing unit. Such schedules are known as permutation schedules. 2. The operation is nonpreemptive. That is, once started, the processing of a product on a unit may not be interrupted and then resume later on. 3. A unit may not process more than one product at a time, nor may a product be processed by more than one unit simultaneously. 4. The last processing unit operates under the UIS policy. That is to say, there is always storage available for a finished product upon its completion on the last unit. 5. The times to transfer batches between units and storage are negligible compared to the batch processing times and are included in them. Similarly, the times required to set up and clean units and storage when changing from one product to another are also negligible. The first assumption may seem restrictive because permutation schedules are not necessarily optimal for serial flowshops or serial multiproduct plants. However, as mentioned by Reklaitis (1982),permutation schedules are more convenient to monitor in practice and hence have enjoyed wide use. Nonpreemptive operation is also desirable since it results in a minimum number of setups.
Determination of Completion Times We begin by presenting recurrence relations for computing completion times and the makespan of a given product sequence. A product sequence is characterized by a permutation of integers: k,-k,-k *..-kN, where product kiis in the ith position in the sequence. Let Cij denote the time at which the ith product in the sequence leaves unit j . Note that C j j is not necessarily equal to the time at which product ki finishes processing on unit j . For instance, when the downstream unit is busy and the storage is full, product ki must be held temporarily in unit j and Cij would then be the time at which product ki actually leaves unit j . To compute Cij, consider a scenario in which product ki has just finished processing on unit j . The transfer of product ki out of unit j is dictated by two events. One is the situation in which either the downstream unit is free of a storage unit is free between units j and (j+l). In either case, product ki can leave unit j,immediately upon completion. The other situation is the one in which the downstream unit is busy and no storage unit is available. In the first event, the completion time of product ki on unit j is simply the time at which unit j starts processing product k j plus its processing time, tkj. But unit j cannot start processing product ki until it has processed the previous product, namely, product k(i-,),or until product ki has been processed by the upstream unit, namely, unit ( j - I ) , as shown in Figure 2. So we have (1) Cij = max [C(i-l)j, Ci(j-,)]+ t k i j In the second event, product ki must be held temporarily in unit j until a storage unit becomes available. Since no storage unit is currently available and the downstream unit is busy, products k+,), k(i-2),...,kowzj)must be occupying the storage facility and product k(j-zj-l)must be in unit
1842 Ind. Eng. Chem. Res., Vol. 27, No. 10, 1988
Ci(j-1) = max
C I(J-I1
I 1-1
UNIT j-i
I
I
I
cij I
UNIT
t
j
1-1
I
I
I
I
1
I- t )
j
I
1
UNIT j
Figure 2. Determination of Cij defined in eq 1.
n
UNIT J+1
UNIT J
PRODUCT KI
Figure 3. Determination of Ci, when storage is full.
(j+l)as shown in Figure 3. The time at which a storage unit becomes available is the time at which product 12(i?.) is discharged to the downstream unit, namely, unit (j+lf. But this discharge is possible only when the downstream -1) leaves unit unit becomes available; i.e., product (j+l).Therefore, Cij = C(i-zj-l)G+l) (2)
From eq 1 and 2, we obtain the following proposition. Proposition I. For an M-unit serial FIS system, the completion times Cij are given by Cij
= max
[C(i-l)j, Ci(j-l), C(i-z,-l)(j+l)-tkjI
i = 1,2,...,N
C(i-l)j-tki( j-1)I
+
tk;(j-1)
+ tkij
j = 1,2,...,M (3)
where Cij = 0, if i I O or j I O or j > M. Note that the above recurrence relations are also applicable to the UIS and NIS policies, since they are in fact special cases of the FIS policy. In the case of UIS, z j N-1; thus, the term C(i-zj-l)Ci+l)-tk*i in the max function is negative and can be omitted altogether. In the case of NIS, although the recurrence relations can easily be obtained by setting zj equal to zero, a closer examination indicates that one of the three terms in the max function of eq 3 is redundant. Setting zj = 0 in eq 3 yields (4) c i j = max [C(i-l)j, CiG-l), C(i-l)(j+l)-tki jl + t k i j
5
Substituting j = j-1 in the above equation gives
= max
[ci(j-l), c(i-1)(j+l)-tki jl + tki j
In terms of complexity, our recurrence relations are easy to code and take at the most O(NM) computational steps. The fact that the completion times can be easily and cheaply computed is critical, because most good heuristics require a large number of makespan calculations. As we will soon see, these recurrence relations form the basis behind the optimal MILP formulation of the scheduling problem. They constitute a single-pass procedure for computing the completion times of the products and do not require any storage feasibility test. By a single-pass procedure, we mean that, once a product in the sequence is scheduled, its completion times are not readusted during the scheduling of the products following it. This is in contrast to the procedure presented by Wiede and Reklaitis (1984a) for the multiunit FIS system in which the available storage is shared among all processing units. Finally, the following lemma is true regarding the recurrence relations of proposition I. Lemma I. The recurrence relations of proposition I yield the minimum makespan for a given permutation sequence of products.
t C‘
UNIT j-i
Ci(j-2),
But this implies that Ci(j-l) 2 C(i-l)p Consequently, the term C(i-l)jin eq 4 is redundant and can be omitted. Thus, for the NIS policy,
i
I
[C(i-l)(j-l),
Optimal MILP Formulation The FIS system under makespan minimization has been shown (Papadimitriou and Kanellakis, 1980) to be an NP-complete problem. Thus, it is not likely that one will ever find an efficient, i.e., polynomially bounded, exact algorithm to solve it. However, in order to characterize the problem, an effort to find an exact but exponentialtime algorithm is still warranted. Furthermore, small-size problems can be solved optimally within reasonable computation time by using an exponential-time algorithm. Most important of all, an exact algorithm can usually provide insights into finding good, efficient heuristics for solving larger problems. We formulate the makespan FIS problem as a mixed integer linear program (MILP). A MILP is an optimization problem whose objective function and constraints are all linear. It differs from a linear program (LP) in that some of its variables are integer or binary (0-1). The notation used in the formulation is the same as before. The continuous variables are Cij as defined previously, and new binary variables are defined as follows:
x..={ 1 IJ
if product i is in position j in the sequence
0 otherwise
Since the criterion is to minimize the makespan, the objective function is minimize CNM The first set of constraints ensures that a product is assigned to one and only one position in the processing sequence: N Xij j= 1
=1
i = 1,2,...a
The second set of constraints ensures that a position in the sequence is assigned to one and only one product: N
CXij = 1
i=l
j = 1,2,...J
(6)
The remaining constraints come directly from proposition I, and they ensure that completion times of a product
Ind. Eng. Chem. Res., Vol. 27, No. 10, 1988 1843 UNIT
Table I. Processing Times
units products 1 2 3 4 5 6
1 10 15 20 14 6 13
2
3
20 8
5 12 9 15
7
6 11 7
5 17
4 30 10 5 10 15 10
Cil
CtklXki k=l
20
- C ( i - z , - l ) ~2 0
2
~
,
I
6
UNIT
,
31
16
,
2 , 6
37
45
78
I A
52
2-
64
d5
,
4
&I
3
A
UNIT a
2
I
I
,
6
38
,
, 3
4
rL.2 il
32
32
72
37
3
31
74
42
Figure 4. UIS schedule of the example of Table I. UNIT1 1 5 1 E
0
‘
I;
I
16
31
(8) I
UNIT2
for units 2 to M-1 ( j = 2, M-1):
3
,
58
~
22
(7)
4
44
i7
N
- c(i-1)1 -
I
6
UNIT2
sequence are correctly calculated. Since the max functions in the recurrence relations are discontinuous, we replace them by multiple, continuous relations as follows: for unit 1 ( j = 1) Cil
’
1 0
1
1
‘
I 17
I
i
’
~
50
65
I 85
4 ~ + %AA
37 43
73
57
92
N cij
- c(i-l)j - 2t k j x k i k=l
20
(9) UNIT
A22
3
N
Czj
- C(i-z, -I)( j + l ) 2 0
(11)
UNIT
(12)
STORAGE
A t ‘ ;
t
4
;
5
i
E
6
A IO1
,
82
92
82
A92
2
,
3
lG2 107
for unit M: N CiM
- C(i-1)M -
tkMXki k=l
2 0
- ci(M-1) -
tkMXki
k=l
1’6, 72
Figure 5. FIS schedule of the example of Table I.
N CiM
I 58
2 0
(13)
where i = 1,2,3,....,N in eq 7-13. Note that since only one of the X k i ’ S for each i will be one, the term E f = l t k , X k i in the above constraints will pick up the processing time of an appropriate product. Note also that for any value of i, if one or no term in a max function of a recurrence relation is nonzero (recall that C i j = 0 for i I 0 or j I 0), then we have an equality constraint instead of an inequality constraint for that recurrence relation and the zero terms do not constitute any constraints. For given M and N, it is clear that the number of binary variables in the formulation is N2,while the number of continuous variables is N X M. The number of constraints depends on the storage policies present between processing units. However, in any case, the number of constraints in the formulation is always less than N(3M-2). Clearly, the above optimal formulation is applicable for solving problems with any mix of the UIS, FIS, and NIS policies. As an example, consider the processing times in Table I. By use of the MILP formulation on LINDO (linear, interactive, discrete optimizer), an optimization package, we scheduled the products in UIS, FIS, and NIS systems. In the case of the FIS system, a storage configuration of z1 = 0, z2 = 0, and z3 = 1is assumed. Optimal sequences of 5-1-2-6-4-3 (makespan 107), 5-1-4-6-2-3 (makespan 107),and 5-6-1-4-2-3 (makespan 111)were obtained repectively for the UIS, FIS, and NIS systems, and their corresponding Gantt charts are shown in Figures 4 , 5 , and 6. Note that the sequences obtained for the three systems are different, which shows that the storage has a significant impact on the optimal solutions obtained. Thus, a heuristic approach available for solving the FIS problem, in which a product sequence is first obtained using existing
UNIT! 1 ’ ; 0
I
E
6
~
29
19
‘
I
64
UNIT
UNIT
3
4
2s
49
A I
I
22
43
1
84
&-
UNIT2 17
’
I
~
43
5
1 37
55
A
w
+;
72
4 70
I
91
2
3
16
, 34
4
, ,
94
135
2
, 3 IO6
if
Figure 6. NIS schedule of the example of Table I.
scheduling procedures for the UIS and NIS/ZW systems and completion times for this sequence are determined later by using proposition I, may not be very good. The reason is that such an approach does not account for the effect of storage in its determination of a product sequence. There are several advantages in formulating the FIS problem as an MILP. Optimization packages such as LINDO are becoming widely available for solving MILP’s. Our experience with LINDO on a VAX 11/780 computer and MPOS (multipurpose optimization system) on a CYBER 180/845 computer indicates that these packages are quite reliable for solving small-size problems. For LINDO on the VAX 11/780, the maximum problem size that one can solve within reasonable computation time is roughly N X M = 30. The computation time per problem ranges from a few seconds to approximately 20 min. The MPOS package appears to be more efficient than LINDO because it has routines that exploit the sparsity of our MILP formulation. (See the Appendix for an example of MILP formulation.) Unfortunately, MPOS on the CYBER is severely restricted by its memory limitation, and we could
1844 Ind. Eng. Chem. Res., Vol. 27, No. 10, 1988
hardly solve any problem beyond N X M = 20. For the FIS system, the MILP approach can be more convenient than the branch and bound (BAB) method because of the wide availability of MILP packages. In the MILP approach, only a simple program needs to be coded for generating the objective function and all the constraints, whereas an entire branch and bound algorithm must be coded in the BAl3 approach. In addition, efficient lower bound estimates, which dictate the efficiency of the BAB method, are not available for the FIS system. As discussed in the next section, the MILP formulation can also be extended to include the ZW policy. Despite its attractiveness, the MILP approach for the FIS system also has a drawback. It is observed that the major work involved in solving an MILP comes from the L P solution rather than the integer solution. This is particularly evident when N is small and M is large, because the number of binary variables is independent of M while the number of constraints is not.
Optimal Formulation with ZW Blocks While storage plays an important role in batch processing by reducing idle times of processing units, some product batches cannot be stored at all between some stages in the processing sequence. These batches are usually unstable intermediates which must be processed by the next processing unit immediately upon completion on the current processing unit. For such processing stages, the ZW policy is most appropriate. Hence, a typical serial multiproduct process consists of a mix of all types of storage policies operating at various processing stages. Although the ZW policy was not incorporated into the recurrence relations presented earlier, the optimal MILP formulation can be easily extended to include ZW blocks embedded in the multiproduct process. First, we redefine completion time Cijto be the time at which the ith product finishes processing on unit j , rather than the time at which it leaves the unit. We will denote this new completion time by C’ij. Consider a scenario in which the (i-1)th product, rather than ith product as before, has just finished processing on unit j ; we wish to determine the earliest time at which the ith product can start processing on unit j . Once this start time is determined, Cij is simply the start time plus t k i j , by our new definition of the completion time. The start time of the ith product on unit j is determined by either the time at which it finishes processing on the previous unit, namely C:+,), or the time at which the (i-1)th product leaves unit j . The transfer of product out of unit j is in turn governed by two events. In the first event, a storage unit or the downstream unit is available to receive product k+.l), in which case the product can leave unit j immediately upon completion. This time is simply given by C’(i-l)j. In the second event, product must wait in unit j , since no free storage is available and the downstream unit is busy. Clearly, products k(i-z),k(i-3),..., k(i-l-zj)must be occupying the storage, and product (k(i-2*.)must be in unit (j+l)as shown in Figure 7. The time at which a storage unit becomes free is the time at which product (i-l-zj) is discharged from storage and starts processing on unit ( j + l ) . This start time of product k(i-l-z.)on unit (j+l)is From the preceding dissimply C’(l-l-~j)(j+l)-t~(i-l-zj~(j+l). cussion, we get the following proposition. Proposition 11. The completion times C’, for a serial FIS flowshop are C’.. ?I max [c’i( j-l), c)i-l)j, c’(i-l-z,)(j+l)-tk(i-l-z,)(j+l)l + t k i j (14)
0 Ki-2
UNIT
j
UNIT
.i+i
1
u PRODUCT Ki-l
0 Ki-l-z.
Figure 7. Determination of C,: when storage is full.
where C:, = 0 if i 5 0 or j 5 0 or j > M. Again, the third and the second term in the max function of eq 14 can be omitted for the UIS and the NIS policy, respectively. Note that C, is always equal to C:j unless product k , must be held temporarily in unit j upon its completion. Since the last product in the processing sequence never has to wait in the last unit due to the assumption of unlimited storage for finshed products, makespans computed from both sets of recurrence relations are always the same, i.e., CNM= ChM. Hence, the two sets of recurrence relations are equivalent as far as makespan computation is concerned. An optimal MILP formulation can again be developed using eq 14 based on the new definition of completion time. Additional constraints can now be imposed on this formulation to accommodate the ZW mode of operation between some stages. Let units P and Q, respectively, be the first unit and the last unit in a block of units operating under ZW. Since in the ZW mode no wait is allowed between any two consecutive processing units, we must impose the following constraints: (15) C’ll = C’l(l-l) + t k , ] i = 1,...,N , j = P+1,...,Q In terms of computational complexity, it is obvious that both sets (eq 3 and 14) of recurrence relations are equivalent, as far as systems without a ZW policy are concerned. One may then ask why we even bothered to present the first set of recurrence relations and the MILP formulation based on it, since the second set of recurrence relations is superior in that th MILP formulation based on it can easily accommodate the ZW policy. A distinct merit of the first MILP formulation is that its constraints are more sparse than those of the second formulation. This is due to the difference in the constraints associated with the third term in the max function of both sets. In the first formulation, the terms t k , ] cancel out, while in the second formulation they do not. So in contrast to the first formulation, this introduces the binary variables in the second formulation. We have mentioned that a large MILP can be solved more efficiently if it is sparse, provided than an MILP solver has routines to exploit the sparsity. Hence, from the standpoint of efficiency, a sparse MILP formulation is generally preferrable to one which is not.
Approximate Scheduling Algorithms Our discussion in the previous sections has indicated that, although the MILP formulation is exact, it is computationally feasible only for small-size problems. Our experience with LINDO shows that the maximum size which we could solve within reasonable computation time
Ind. Eng. Chem. Res., Vol. 27, No. 10, 1988 1845 on a VAX 11/780 is about N X M = 30. In chemical batch plants, it is not uncommon to encounter the need to schedule more than 10 products across a sequence of processing stages, which renders the optimal MILP approach or for that matter any optimal approach, computationally infeasible or unsuitable. Hence, there is a need to find efficient algorithms which can give good, suboptimal solutions to large problems in reasonable computation time. In general, two options are available for finding efficient, suboptimal algorithms. One option is to devise one’s own scheduling algorithm, taking into consideration the features unique to a system. The other is to modify and use the existing algorithms available for solving other systems. Since the second option is more economical, several workers have applied various existing UIS heuristics to solve different FIS configurations with some success. We investigate both of these options. First, we evaluate the performance of an effective, existing heuristic for the UIS system on the FIS system, and second, we present and evaluate an approximate algorithm based on the optimal MILP formulation.
RAES Heuristic RAES or rapid access with extensive search was proposed by Dannenbring (1977) for solving the multiunit UIS makespan problem. It is a two-phase heuristic with a recursive improvement strategy. In the first phase, known as RA, a pseudo-two-unit UIS problem is generated from the original M-unit problem as follows: M
t*il = z(M-j+l)tij j=1
M
t*iz = x j t i j i = 1,2,...A (16) j=1
where t*ij is the processing time of product i on unit j for the pseudo-two-unit problem. This two-unit problem is then solved optimally by using Johnson’s algorithm (1954) to obtain a good, initial sequence. Johnson’s algorithm (see Ku et al. (1987)) is a simple but exact procedure to solve the two-unit UIS system under the makespan criterion. In the second phase, known as ES, this initial sequence is improved by generating (N-1) sequences via pairwise interchanges of adjacent products in the sequence and selecting the best of those as the next solution. The procedure is repeated until no improvement is found. Even though RAES was used for the UIS system, the basic strategy can be applied to the FIS system as well. The only change required is that, instead of using the UIS recurrence relations to obtain makespans, the FIS relations are used. We evaluated the performance of RAES on small-size FIS problems for each possible combination of N = 3, 4, 5, 6 and M ranging from 3 to 10. For each combination, 50 test problems were solved. The processing times were generated randomly from a uniform distribution in the range 0-99, and the number of storage units between two processing units was randomly generated in the range 0-5. Average percentage deviations based on the optimal solutions obtained through complete enumeration are shown in Table I1 with an overall average deviation of 0.81%. In addition, the percentages of the proportion optimal, which are the proportion of solutions equaling the optima, were also computed for each combination, and their values are shown in Table 11. Table I1 clearly indicates that RAES is effective for solving small FIS problems. The reason is that, when N is small, the FIS system behaves like the UIS system (not much interaction between the products and the storage), and thus RA provides a near-optimum starting solution for the FIS
Table 11. Evaluation of products 3 4 5 6 av
RAES on Small-Size FIS Problems % dev.n
0.24 0.48 1.02 1.49 0.81
prop. optimal: % 91.8 82.3 63.3 50.0 71.8
a % deviation = [(heuristic makespan - optimal makespan)/optimal makespan1100. *Proportion optimal = (number of solutions equalling optima/number of test problems) 100.
system as well. However, this may not be true as N increases. Note that the deviations become larger and the proportion optima become smaller as N increases. This is an indication that RAES may not be a good heuristic for solving larger FIS problems. The reason is that the use of an existing algorithm such as RAES does not exploit the special features unique to the problem, which in our case are the constraints on storage. In fact, the pairwise interchange strategy in RAES is random from the standpoint of optimization. Hence, it is clear that one should be able to develop an approximate scheduling algorithm that is more effective than RAES for the FIS system.
Approximate MILP Heuristic The MILP formulation is a good candidate and could be very effective for solving the FIS problem heuristically, because it inherently incorporates all the storage usage within the formulation. In principle, it is possible to devise an approximate algorithm based on the exact one. This is often done to increase the efficiency of the algorithm, especially if its time requirement increases exponentially with the problem size. In the case of the FIS system, larger problems can be solved based on the optimal MILP formulation if the computation time can somehow be reduced at the expense of accuracy. As discussed below, we investigated a number of strategies to achieve this and used RAES as a basis for comparison. One strategy is to relax integer requirements and solve only for the linear programming (LP) solution. Since the solution obtained this way may be fractional, i.e., 0 I Xi, I 1, the X i j can be rounded off to get an integer solution or a permutation. A second strategy is to first obtain an LP solution, but then instead of rounding off the L P solution, we now assign each product to only those positions of the sequence where the X,’s of the LP solution are nonzeros. A third strategy to increase the efficiency of the MILP formulation is to solve a series of small MILP’s. This can be done be dividing all products into a few groups and then solving these groups as separate MILP problems. The final product sequence is simply a concatenation of the sequences obtained from solving these separate MILP’s. Unfortunately, preliminary testing of these methods gave comparable results to those from RAES. So their requirements of higher computation time are not warranted. An obvious extension of the last approach is to avoid confining products in groups but confine every product to only certain positions instead of all positions as required in the optimal formulation. This forces us to include all products in a single MILP formulation. For instance, if N = 5, product 1 may be assigned to positions 1, 2, and 3 in the sequence. Thus, the binary variables X I 4and X I 5 are omitted altogether in the formulation. Similarly, we could do the same for other products. The result is that the number of constraints remains unchanged in the formulation, while the number of binary variables decreases considerably and the formulation becomes ‘more sparse.
1846 Ind. Eng. Chem. Res., Vol. 27, No. 10, 1988 Table 111. Comparison of RAES and Approximate MILP
CPU time 70 dev." of RAES
NIM 613 613 613 813 10/3
operating modes NIS
mean
FIS UIS
1.5 1.2 2.1
range 0-10.1
mean 0.9 0-8.8 0.3 0-12 0.6 0-7.0 0.75 av improvement over RAES = 1.35%
1.8
FIS FIS
% dev. of MILP
% Deviation = [(heuristic makespan - optimal makespan)/optimal makespan1100.
Chart I
1
1 ' X 2 / x 3 ' X 3 4 ;
.
I
I
Nl
2
3
x x x x
x x x x x
possible positions 4 5 6
x x x x x
x x x x x
x x x x x
*
x x x x x
*
.
N
x x x x x x x x x x x x
This concept of allocating products to various neighborhoods of a sequence stems from our observation that RA provides a good starting solution for solving the UIS problem by assigning products to the neighborhoods of the their positions in the optimal solution. The major virtue of RA is that its structure is based on Johnson's algorithm, Le., products with small processing times on their first few units are scheduled early in the sequence, while products with small processing times on their last few units are scheduled later in the sequence. This means that products, depending on their processing times, belong to certain neighborhoods in the sequence. Since pairwise interchange is basically a neighborhood search strategy, it is not surprising that RAES has been the best procedure for solving the UIS problem. The approach is certainly approximate, and its performance will depend on our ability to obtain a good starting product sequence. For our purposes, a good, initial solution should be near-optimal and inexpensive to compute. RA seems to be the best choice, because it also provides a good, initial solution for the FIS problem when N is small. In addition, RA ranks second both in terms of effectiveness and efficiency among the five heuristics that generate a single solution in Dannenbring's (1977) evaluation. On the other hand, among the five, the most effective heuristic, called merging, ranks fourth in terms of computation time. Thus, RA offers the best compromise between effectiveness and efficiency. The approximate MILP algorithm can then be summarized as follows: (1)Obtain a starting product sequence from RA. (2) Assign each product from the RA sequence to the positions in Chart I. Each element in the matrix corresponds to an Xi, in the formulation. For example, product 1 from the RA sequence is assigned to positions 1 , 2 , and 3, while product 4 is assigned to positions 2, 3, 4, 5, and 6. (3) Formulate an MILP with all products assigned according to the above diagram and solve it. The number of binary variables can be reduced even further by not declaring the diagonal elements in the matrix as integer variables. The diagonal elements will eventually be forced to be integers because the X,'s must sum to 1. Thus, instead of having N 2 binary variables in
range 0-4.1 0-2.1 0-6.3 0-2.8
(9)
of
MILP mean 50 50 26 167 394
rangeb 21-82 16-103 11-61 35-357 76-866
*VAX 11/780- VMS.
the formulation, we now only have (4N-6) (N L 4). The method was tested on a set of 100 problems with 20 problems for each of the configurations shown in Table 111. The processing times were generated randomly from a uniform distribution in the range 0-99. For the FIS systems, the number of storage units between two processing units was generated randomly in the range 0-5. As a basis for comparison, we also solved the same problems by using RAES. The approximate MILP approach gave makespans which were on the average within 1%of the optimum makespans obtained via the exact formulation, while the RAES makespans had deviations (from the optimal) which were more than twice those of the approximate MILP. The maximum deviations of RAES were also more than twice larger than those obtained by the MILP approach. These observations are true even for the UIS systems for which RAES has been the best procedure. For N I M = 1013,the optimal solutions could not be obtained easily by the exact approach, and hence only an average percentage improvement of 1.4% over the RAES solution was calculated. As shown in Table 111, the average computation time of the approximate MILP ranges from half a minute to about 7 min using LINDO on a VAX 111780 computer. This represents an approximately 50% reduction in computation time as compared to the exact approach and thus extends the range of solvable problems using this approach to N X M = 60. The major work of the heuristic comes from obtaining the LP solutions, and considerable savings can be obtained if routines that exploit sparsity are employed. Finally, we observed that, of all the 100 test problems solved, not even a single solution obtained from RAES was better than that obtained from our heuristic. This clearly indicates the superiority of the approximate MILP heuristic over RAES. Conclusions In this work scheduling of N products through an M-unit serial batch process with no or finite intermediate storage was studied under the criterion of makespan minimization. Recurrence relations requiring O(NM)simple computational steps were presented for computing completion times of all products in a given product sequence. Based on these relations, a mixed integer linear programming formulation was developed for solving the problem optimally. The formulation allows any mix of the four storage operating policies, UIS, FIS, NIS, and ZW. An effective heuristic based on the MILP formulation was also developed. Because of the NP-complete nature of the FIS problem, the optimal MILP formulation is feasible for solving problems up to a size of N X M = 30 by using LINDO on a VAX 11/780 computer. However, more than 50% savings in computation time can be achieved by the approximate approach, thus increasing the range of solvable problems to N X M = 60. The approximate algorithm is quite effective as it gives solutions within 1% of the optima and performs much better than RAES on the average in
terms of the quality of solution. The MILP formulation presented in this paper is general and can be extended to accommodate semicontinuous units. The proposed approach, especially the approximate MILP, is effective and convenient for small-sizeproblems, as MILP packages such as LINDO are widely available. Work is currently under way to develop heuristic approaches that will be feasible for solving large problems with N X M > 60. The results presented in this work are useful in increasing equipment utilization and productivities of batch plants because most batch processes make use of various forms of intermediate storage. Acknowledgment We acknowledge the financial support from Northwestern University and IBM Corp. Nomenclature Cij = time at which the ith product leaves unit j C$j = time at which the ith product finishes processing on unit
i
M = number of processing units max [ ] = maximum of the quantities in the brackets N = number of products tij = processing time of product i on unit j t * i j = weighted processing time of product i on unit j in RA Xij = binary variable defined in the optimal formulation zi = number of storage units between units i and (i+l) Abbreviations BAB = branch and bound method FIS = finite intermediate storage LINDO = linear, interactive, discrete optimizer LP = linear programming MILP = mixed integer linear programming MIS = mixed intermediate storage MPOS = multipurpose optimization system NIS = no intermediate storage NP-complete = nondeterministic polynomial-time complete NW = no wait RA = rapid access heuristic RAES = rapid access heuristic with extensive search UIS = unlimited intermediate storage VMS = virtual memory system ZW = zero wait Appendix: Example of an MILP Formulation Consider a four-unit, three-product scheduling problem whose processing times are given below. Let z1 = 1,z2 = 23 = 0 with ZW mode of operation in effect between units two and three. An optimal MILP formulation to minimize the total time required for production is as follows: processing times for units 1 1 2 3 4 20 15 8 11 14 18 9 I 25 17 4 30 1 6
An optimal sequence of 3-2-1 was obtained from LINDO with a makespan of 85 and the following completion times:
3 2 1
I
products 1 2 3
completion times for units 2 3 4 1 6 23 27 57 I 31 45 63 72 51 66 74 85 I
products
,
I
1
I
I
;
Objective function: minimize C’34. Constraints. For binary variables, Xi1 + Xi2 + X13 = 1 Xzl+ X22 + X23 = 1 X31+ X32 + X33 = 1 XI1 + XZl + x31=1 Xi2 + X22 + X32 = 1 X13 + X23 + X33 = 1 For unit 1,
Literature Cited Dannenbring, D. G. “An Evaluation of Flowshop Sequencing Heuristics” Manag. Sci. 1977, 23, 1174. Dutta, S. K.; Cunningham, A. A. “Sequencing Two Machine Flowshops with Finite Intermediate Storage” Manag. Sci. 1975, 21, 989.
Graves, S. C. “A Review of Production Scheduling” Opns. Res. 1981, 29, 646. Johnson, S. M. “Optimal Two and Three Stage Production Schedules with Set-up Times Included” Naval Res. Logist. Q. 1954,1, 61. Ku, H. M.; Rajagopalan, D.; Karimi, I. A. “Scheduling in Batch Processes” Chem. Eng. B o g . 1987, Aug, 1. Papadimitriou, C. H.; Kanellakis, P. C. ‘Flowshop Scheduling with Limited Temporary Storage” J. Assoc. Comput. Mach. 1980,27, 533. ReMaitis, G. V. “Review of Scheduling of Process Operation” AIChE Symp. Ser. 1982, 78 (214), 119-133.
Ind. Eng. Chem. Res. 1988,27, 1848-1862
1848
Wiede, W., Jr.; Reklaitis, G. V. "Determination of Completion Times for Serial Multiproduct Processes, Part 2: Multiunit Finite Intermediate Storage System" Presented at AIChE Annual Meeting, San Francisco, 1984a; Paper 123c. Wiede, W., Jr.; Reklaitis, G. V. "Determination of Completion Times for Serial Multiproduct Processes, Part 3: Mixed Finite Intermediate Systems" Presented a t AIChE Annual Meeting, San Francisco, 1984b; Paper 123c.
Wiede, W., Jr.; Kuriyan, K.; Reklaitis, G . V. "Determination of Completion Timea for Serial Multiproduct Processes, Part 1: Two Unit Finite Intermediate Storage System" Presented at AIChE Annual Meeting, San Francisco, 1984; Paper 123c.
Received for review August 12, 1987 Revised manuscript received March 23, 1988 Accepted April 15, 1988
Understanding the Dynamic Behavior of Distillation Columns Sigurd Skogestad* and Manfred Morari California Institute of Technology, Chemical Engineering, 206-41, Pasadena, California 91125
T h e dynamic behavior of a distillation column is approximated with a two time constant model. The response to changes in the external flows is approximately first order with time constant 7 1 . This dominant time constant can be estimated by using a simple mixing tank model for the column. The response to changes in the internal flows is also first order, but ita time constant T~ is generally significantly smaller than 71. The condition number and the RGA are smaller a t high frequency than a t steady state. Most models presented in the literature do not take this into account. The two time constant model does predict this behavior, and 72 can be estimated by matching the RGA at high frequency. Finally, it is shown that the effect of nonlinearity is almost eliminated if logarithmic compositions In XB and In (1- y ~are ) used. In particular, this applies to the initial response which is of primary importance for feedback control. 1. Introduction This paper is mainly concerned with understanding the composition response of distillation columns (Figure 1). Modeling the dynamic response of distillation columns is conceptually straightforward. Consider the simplest case with binary separation, constant relative volatility ( a ) , 100% tray efficiency, constant molar flows, and constant vapor or liquid holdups. The dynamic behavior of a column with N theoretical trays and a total condenser is then described by N + 1 first-order differential equations (see Appendix). For tray i (which is not a feed tray, reboiler or condenser), Mi dxi/dt E Miii = Lxi+l + Vyi-1 - Lxi - V y i (1)
where yi =
CuXi/(l
+ ( a - l)xJ
(2)
Solving this set of differential and algebraic equations on a digital computer is also straightforward (although computer times may be excessive for columns with a large number of trays and components). An early discussion on this is found in Rosenbrock (1957). However, even though the governing equations themselves are well-known, the understanding of the overall behavior of this set of equations is not. I t is known that the response is dominated , may be estimated by by one large time constant, T ~which some "mixing tank" approach (Davidson, 1956; Moczek et al., 1963, 1965; Wahl and Harriot, 1970; Skogestad and Morari, 1987a). However, this dominant time constant generally does not apply when we make changes in the internal flows only (simultaneous increase in L and V). This will become evident from the discussion in section 3. Furthermore, because of the nonlinear VLE relationship (eq 2), the response turns out to be strongly nonlinear, in particular for high-purity separations, and the value of the dominant time constant T~ may change drastically with
* Address correspondence t o this author at Chemical Engineering, Norwegian Institute of Technology (NTH), N-7034 Trondheim, Norway. 0888-588518812627-1848$01.50/ 0
operating conditions (e.g., Skogestad and Morari (1987a)). The main objective of this paper is to gain insight into the dynamic behavior of distillation columns by developing simple analytical models. 'For feedback control, an accurate model of the plant is usually not needed. Rather, a model which includes the factors most important for feedback control (inverse responses, multivariable effects, and sensitivity to model uncertainty) is desired. Pinpointing these factors is most easily accomplished with simple analytical models. We stress that these models are by no means intended to replace nonlinear simulations. Also, we stress that it is desirable that these models be analytic. Consequently, the models are not intended to replace low-order models which may be obtained by employing standard model reduction techniques (including residualization, balanced realization, collocation methods, etc.). Internal and External Flows. Throughout the paper, we make use of the terms external and internal flows. A change in external flows is any change which changes the ratio D/B, for example, an increase in flux L with constant boilup V. An increase in internal flows is accomplished by a simultaneous increase in L and V ,while keeping the product rates, D and B , constant. Linear Input-Output Model. For feedback control, we need an overall model that describes the effect of the inputs (flows) on the outputs (product compositions, Y D and x g ) . All results in this paper (gains, RGA values, etc.) are with reflux ( L )and boilup ( V )as inputs (LV configuration). We assume perfect level control. Distillate (0) and bottom flows ( B ) are then manipulated to keep constant holdups in the condenser (accumulator) and the reboiler (column base). With the additional assumption of constant molar flows and with no feed disturbances (dF = 0 ) this yields dD=dV-dL dB=-dD (3) This does not imply that the LV configuration is the preferred choice for control purposes. The choice is made because L and V have a direct effect on compositions and their effect is therefore only weakly dependent on the 0 1988 American Chemical Society
