Ind. Eng. Chem. Process Des. Dev. 1986, 25,979-908
979
Scheduling and Sequencing of Batch Operations in a Multipurpose Plant Steven H. Rich and George J. Prokopakls' Department of Chemical Englneerlng and Applied Chemistry, Columbia University, New York, New York 10027
A mixed integer programming approach to schedulingand sequencing multiple production runs of N products on M processors over a finite time horizon is presented in this paper. The scheduling aspect determines the production quantities, while the Sequencing aspect establishes the job processing order at each unit. The model's flexibility
allows for a choice of scheduling/sequencing objectives. LINM), a library quality mixed integer programming computer package, is used to solve the problem. Two sample problems are considered. The approach yields good solutions with reasonable computational effort.
In recent years, plants involving batch operations have received their share of attention in the literature. In particular, researchers have developed models to optimize the design of multiproduct and multipurpose batch facilities with respect to the number and size of equipment units (Sparrow et al., 1975; Grossmann and Sargent, 1979; Knopf et al., 1982; Suhami and Mah, 1982). However, only limited attention has been paid to the problem of devising production control systems for existing facilities with batch operations (Prabhakar, 1974; Overturf et al., 1978; Maudedi and Rippin, 1980; Suhami and Mah, 1984). A production control system attempts to effectively coordinate the availability of labor, equipment, and materials. In a multipurpose plant, the system must consider numerous (short or intermediate range) planning activities including (a) forecasting the demand, (b) controlling the inventory control, (c) determining the quantities of products to be produced in each production run (scheduling), and (d) establishing process ordering a t each processor (sequencing). The scheduling problem concerns itself with assigning a set of limited resources (idle processors) to a set of competing tasks (product production). This problem is further complicated when the sequencing (timing) of the operations on the processors is considered in conjunction with the scheduling activity. As indicated in a recent review article on production planning (Graves, 1980, few models have been presented in the Operations Research/Management Science literature which incorporate both scheduling aspects of a production control system, though there is an abundance of work which studies these problems separately. This paper addresses a mathematical programming model for a multipurpose batch facility in which both scheduling and sequencing activities are simultaneously considered. A qualitative description of the problem is presented in the next section, followed by a detailed analysis of the model formulation. Next, an overview of the proposed solution methodology is presented, and finally, the application of the model to two sample problems is discussed in the last section.
Problem Description The model presented in this paper pertains to a multipurpose plant in which N products are produced. Each product requires a subset of the M processors in the plant. The processors in this subset are deployed simultaneously and production of a product on a processor cannot be preempted once processing has commenced. The batch
* Author to whom correspondence should be addressed. 0196-4305/86/1125-0979$01.50/0
size associated with a product-processor pair is fixed. There are two product types, saleables and intermediates. Saleables are products for which a due date has been specified. A product may be classified as both saleable and intermediate if it is used in the production of another product. Unlike a multiproduct plant in which the production of all saleable end products follows a fixed routing through the plant, the processing sequences for the saleable products in the multipurpose plant differ. Moreover, different products may be produced simultaneously throughout the plant, and processing sequences for successive production runs of a saleable product may vary. Precursor relationships among the products are specified by a network structure. The model is appropriate for a planning horizon of either short or intermediate range. Within the planning horizon, due dates, Dik, are specified for completion of the kth order of a saleable product i. The j t h production run of a product i represents the processing of Nil consecutive batches of that product. The duration of this task is given by tiNijwhere t, is the processing time per batch of product i. The maximum allowable number of production runs for a product is set equal to the number of orders of that product. The kth customer order for a saleable product i ( X i kunits) is met from the quantity of that product produced in production run k and/or from inventory. Demand for each product is known, though it is not necessarily uniform over the planning horizon. The saleable products demands, in conjunction with stoichiometric considerations, determine the demands for intermediates. The purpose of the model is to select the number of batches of saleable and intermediate products to be produced in a production run so as to satisfy customer demand over the planning horizon. In addition, the model specifies operation starting times, so that processing order of operations on any batch unit may be inferred. Binary variables are incorporated in the model to prevent simultaneous production of two products in a batch unit. There are several schedulingfsequencing criteria from which the decision maker can choose to measure the effectiveness of a solution. These range from minimizing the mean tardiness of saleable products to minimizing the idle time on a processor.
Model Formulation Network Structure. The network structure which constitutes the precursor relationships among the products is a set of node-disjoint, acyclic digraphs (a cycle-free set of vertices and directed arcs, grouped into subsets of vertices which are disconnected). Such a network structure 0 1986 American Chemical Society
980
Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 4, 1986
@ %
0 Figure 1. Example of the precedence network structure. Products 4,5 , and 7 are saleable. All others are intermediates.
is shown in Figure 1. Each digraph specifies the sequence of operations for the production of a saleable product. A node of a digraph represents the task of producing a product. A task duration and a set of processors are associated with each node i. An arc (ij)in a digraph represents a precedence relationship between tasks i and j ; i.e., task j cannot begin until the completion of task i. The digraphs must be assembly structures; that is, each node may have a t most one immediate successor, though the node may have multiple immediate predecessors. Therefore, a production sequence network in which a product is an intermediate for more than one saleable product, as shown in Figure 2a, must be modified prior to using the model. The general structure shown in Figure 2a can be decomposed into two assembly structures as shown in Figure 2b. Similar logic holds if a product can be classified as both saleable and intermediate. For complex assembly structures, the decomposition algorithm of Afentakis and Gavish (1985) can be used. The total number of nodes in the precedence network is N’. Model Assumptions. The following limiting assumptions are made to reduce the complexity of the problem: (a) Storage tanks of sufficient capacity are available to meet inventory requirements. (b) Products are not perishable over the length of time they are stored. (c) Inventory costs are not considered. (d) Changeover costs and changeover times are independent of job sequences on the processors. Moreover, setup times for consecutive batches of a product are considered to be incorporated in task durations. Model Parameters. The model requires the specification of the following parameters, classified into three groups: Group A: User (Internally) Controlled. (a) Due dates for the kth order of saleable product i, D& (also under customer influence). (b) The maximum number of production runs per saleable product i, R(i). In this paper, R(i) is set equal to the number of customer orders for saleable product i. The number of runs per intermediate is equal to the number of runs of saleable product i for which the intermediate is produced. (c) The processor batch size for product i, B,. (d) The objective function weight which gives the (relative) penalty cost for each unit o i time that the kth order of saleable product i is tardy, wik. (e) The initial inventory for product i, Ii,. Group B: Customer (Externally) Controlled. (a) The size of the kth customer order for saleable product i, x i , (expressed in the same units as Bi). The sum of the demands x i k must be met by the end of production run n of product i.
*
multi-copied portion of netNork
2‘
Figure 2. Decomposition of a general structure network (a, top) a general structure network; (b, bottom) decomposed general structure network.
Group C: Process Parameters. (a) The cycle time ti (setup + filling + holding + drainage time) for a processor (or processor set) required by product i. (b) The stoichiometric factors f i j , i.e., the number of units of product i requires to produce one unit of the immediate successor product j . There is an f i j associated with each arc (ij)in the precedence network. (c) The postprocessing time required to ready a batch of product i for production of its immediate successor product j , uiP If the set of nodes (i = il, i.,, ...,i, = j ) represents a path in the precedence network from node i to (saleable product) node j , then (ngl’fi,i,,,)x$=:Xjkrepresents the production requirement for intermediate product i which must be met by the end of the product’s nth production run. To simplify the notation, this expression is replaced by C;=IXik*
Decision Variables. The model determines values for the following sets of variables. (a) Start time, S i k , of k = 1,...,R(i). From these start product i, run it; i = 1,...8’; times, the job sequence on a processor may be inferred. (b) Number of batches, Nik,of product i produced (consecutively) in production run k ; i = 1,...,N’; k = 1,..$( i). Nik may take on any nonnegative value. A task duration is given as tiNik and the amount of product i produced in run k is BiNik. ( c ) Binary variables, Yikjk’, which resolve conflicts on processors. (d) Tardiness Tik+ of the kth order of saleable product i; iEA; k = 1,...,R(i). Tardiness is defined as the time by which the completion time Fikof a saleable product exceeds its due date Dik. If the due date is met, the tardiness is equal to zero. Hence, Tik+= max [o, F i k - Dikl. Therefore, the number of variables is 2xE1R(i) + &AR(i) + Bin. Objective Function. Different criteria may be used to evaluate a production schedule. The quantitative expression of these objectives follows. (i) Minimize mean weighted tardiness. It is equivalent to R(i)
min
CwikTik+ iEAk=l
(ii) Minimize maximum tardiness of saleable products.
Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 4, 1986
The objective function will read min 2 = max [Tik+]
This relationship is established by appending the following set of Constraints to the formulation: 2 1 Tik+ iEA k = 1,...$(i) Although these constraints allow 2 to take on any value greater than or equal to max [Tik+], 2 w i l l equal max [Tik’l since 2 is being minimized. (iii) Minimize mean completion time of saleable products. The objective function will read R(i)
min C C (Sik + Nikti) i€Ak=l
No additional constraints are required. (iv) Minimize makespan, the time required to produce all saleable products. The objective function is min Y’ = max [Sip(i)+ Ni3&] iEA
For each saleable product (iEA) the following set of constraints must be added: Y’ 2 Sis(i)+ NiR(i)ti (v) Minimize idle time on processor k. Greenberg (1968) shows that achieving this goal is equivalent to minimizing the completion time of the last job on processor k: min Gk = max + Nip&] GLJk
For all products (iEUk),the following set of constraints should be added: Gk
si,R(i) +
(vi) Minimize idle time on all processors. It follows from the preceding discussion that achieving this goal is equivalent to M
min CGk k=l
with the constraints
Si,R(i)+ Ni,R(i)ti k = 1,...N i E Uk Gk
Multiple Solutions. It is conceivable that the chosen objective function may attain its best value at two or more solutions; i.e., multiple optimal solutions exist. Although these solutions are considered to be equivalent based on a single criterion, one solution may prove to be superior to the others, when the solutions are judged with secondary criteria as well. To distinguish among the alternate optima, one might employ a preemptive priority goal programming scheme, in which the problem is sequentially solved with a ranked set of objective functions. The decision maker chooses a subset of the schedulinglsequencing criteria and prioritizes them depending on their relative importance. The problem is fnst solved with the objective function associated with the criteria deemed to be of greatest importance. Once the optimal objective value is found, an upper bound constraint of the form (objective function of priority i) 5 (optimal objective function value under priority i) is appended to the formulation, and this revised problem is then re-solved with a new objective function, namely the objective function corresponding to the criteria of next greatest importance. Subsequently, an upper bound constraint on this second objective function is added to the formulation. This process continues until it is determined that the current
981
optimal solution is unique or the chosen set of objective functions of secondary importance is exhausted. The last optimal solution found by this approach (assuming that it is unique) is optimal under the first criterion and dominates all other optimal solutions (to the initial objective) based on the secondary criteria. Constraint Sets. The constraints of the problem have been grouped in seven sets according to their function. I. Precedence Constraints. There are two types of precedence constraints: (a) the ones dictated by the physical process itself and (b) the ones dictated by the scheduling needs. Type a. Precedence between two jobs i and j is established through an arc in the precedence network: s j k 2 Sik + Nikti + UijN*k (i,j)Enetwork
k = 1,...,R( i)
where N*k = Nik if d l Nik batches of product i just produced are postprocessed together and Nak= Njk if Njk batches of product i are postprocessed. If Nik> Njk,then the remaining Nik- Njkbatches are stored and are postprocessed immediately prior to their use in producing product j in subsequent production runs. If Nib < Njk, then production of product j cannot begin until additional Njk - Nik batches of product i are taken from inventory. Whether Nik or Njk is chosen as N,k is a t the decision maker’s discretion. I t is assumed that conflicts will not arise a t the postprocessing equipment. Alternatively, if the postprocessing time is independent of the size of the batch, the last term of the right-hand side of the inequality may be replaced by a constant. The significanceof these constraints is that a run cannot commence until the intermediates (if any) are available. There are &JVwc,$(i) of these constraints in the formulation, where Nwc,is the number of arcs in the digraph associated with saleable product i. Type b. Precedence is established by production runs
Si,k+l 2 S i k + Nikti i = 1,...,N’ k = 1,...,R(i) - 1 That is, the (k + 1)stproduction run of a product cannot precede the kth. Such a constraint is unnecessary if R(i) = 1;i.e., only one production run of i should be scheduled. There is a total of [R(i)- 11constraints of this type. 11. Demand Constraints. These constraints take the form of
IiO
k
k
p=l
p=l
+ Bi x N i p 2 c xi,
i = 1,...,N’
k = 1,...,R(i)
They require that demand through the kth order of product i must be filled from inventory or production run k. There are CKIR(i) constraints of type 11. If k = 1, the demand constraint takes the form Iio &Nil L Xi,
+
which gives the lower bound constraint Nil 2 [Xi1 - Iiol/Bi Expressing Nil in terms of a new variable Nil’ Nil = Nil’ [Xi, - Iio]/Bi
+
and substituting for Nil throughout the formulation transforms the lower bound constraint to a nonnegativity constraint Nil’ 1 0, thereby reducing the size of the formulation by N’ constraints.
982
Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 4, 1986
111. Production Limitations Established by Reactant Availability. The production quantity of product j will be limited by the amount(s) of immediate predecessor product(s) i available. Consequently, for each arc (ij7 in the precedence network, a contraint of the following form is required: BjNjk 5 [Xik + l i k l ( l / f i j ) Substituting for I i k gives k
k-1
i = 1,...,R(i) - 1
Note that a constraint is unnecessary for k = R(i) since BjNjp(i)will be equated to the remaining unsatisfied demand not yet produced. Moreover, these constraints are not needed for digraphs in the network corresponding to single production run saleable products. There are CiEANarc,i(R(i) - 1) such constraints. IV. Due Date/Finishing Time Relationships [(si,
+ Nikt;) - Dik] - T i k + 5 0 i€A
1 if an n-level path connects nodes i and j = 0 otherwise
{
Define P* as L
k = 1,...,Hi)
These constraints equate tardiness to max [0, Fik - &]. Fikis the finishing time of the kth run of product i and
can be equated to S i k + Nikti. Note that if the due date is met, i.e., the finishing time is less than Dik, then Tik+ may take on any nonnegative value. However, Tik+will be set equal to zero, since a penalty cost is incurred whenever Tik+ is greater than zero. By similar logic, it follows that Tik+will equal the difference between the finishing time and the due date whenever the finishing time exceeds the due date. This set of constraints is unnecessary if one of the objective functions (iii)-(vi) is chosen. There are CiEAR(i) constraints of type IV. V. Disjunctive (either/or) Constraints. No two operations can be processed simultaneously in a reactor unit. For each pair of products that are in potential conflict, a pair of disjunctive constraints is introduced to determine job ordering in a processor. Two products are in potential conflict if the two products require the same processor k (i.e., i j E U k ) and either of the following conditions is met: (i) A directed path does not connect nodes i and j in the network; Le., a production precedence relationship does not exist between products i and j . (ii) A directed path does connect nodes i and j , and the precedence constraints of set I do not give a process ordering between i and j . This situation arises if the respective starting times of production run k of the product i and production run k'of product j are compared, where i # j and k > k'. If k 5 k', then process ordering is established by the relationships S ' k ' 1 s j k (by production run precedence) and S j k IS i k (by precedence network). In order to construct the disjunctive constraints, it is necessary to identify all directed paths in the precedence network in a systematic fashion. Some background information from graph theory will be useful a t this stage. Let P' be the N' X "adjacency matrix which gives all one-level paths (an n-level path is a directed path consisting of n arcs) in a precedence network. Its elements are 1 if arc (ij)Enetwork
L
P* = CP" = C(P')"
i = 1,...,R(i) - 1
(ij)Enetwork
(ij)Enetwork
In general, P"is the precedence matrix for n-level paths with elements
n=O
n=O
where Po is the N ' x N ' identity matrix and L is the number of arcs in the longest path in the network (PL+l is the null matrix). An element of P* gives the number of paths in a digraph connecting nodes i and j . If pi,* = 0, then no directed path exists between i and j . It follows from the relationship P* = Cf;=,(P')"that P* = (I - P1)-l. This last expression provides a means of determining P*. Introducing the connection matrix Q, defined as
the conditions for potential conflict can be restated as follows. Given the quadruple (i,pi,p'), a conflict may result on a process if (i) there exists a k such that ij€Uk and (ii) one of the following three mutually exclusive conditions are satisfied: qij = 0
i#j
qi,=l
p>p'
i#j
q;i=1
Pep'
Evaluation of the Q matrix (Q = P* when all digraphs are assembly structure networks) and consideration of the above criteria in a preprocessing step lead to the identification of the disjunctive constraints that need to be included. They assume the form s j k , - S i k + w(1 - Yik;k') 2 Nikti
where 1 if product i, run k precedes product j , run k ' 0 otherwise
{ and W is a large number (relative to any task duration yikjk'=
Nikti), specified by the decision maker, with upper bound R(i) N'
[CXikti
k=l ill
+ ( i j ) €network uij(xik/Bi)l
The lower bound on W is the smallest value for which the optimal starting times satisfy all disjunctive constraints. The rationale for these constraints is straightforward. If i precedes j , then it is required that s j k ! 1 S i k + Nikti. If j precedes i on processor k, then it is required that S i k 1 S,, + Njkd;. If &jk' = 1,then the first constraint of the 1S i k + Nikti, thus preventing simultaneous pair gives production of i and j . The second constraint of the pair gives S i k + W 2 s j k , + Njkd, which is trivially satisfied if W is properly set. Similar reasoning holds for Yikjk' = 0. The number of disjunctive constraints is twice the number of binary variables. The number of binary variables is primarily a function of (a) the demand for each processor, (b) the number of products, and (c) the number of production runs for each product. It is desirable to choose Was small as possible so as to tighten the formulation. However, it is impossible to know
Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 4, 1986 Q83
S u b j e c t to: Precedence: 'jk 'i,k+l
'ik
1.
Nikti 'ik
+
Nikti
j)Enetwork k=l,...,R(i)
V(i,
,...,
N' i=l k:l,...,R(i)-l
Demand: Vi; k=l.....R(i)
Availability: BjNjk
k. The following relationships exist = 0 implies
Yikjk'
T o t a l n u m b e r of c o n s t r a i n t s :
Figure. 3. Summary of the mathematical programming formulation.
the lower bound on W a priori. Low values of W are desirable since a judicious choice of W can immediately eliminate certain poor binary variable combinations, in that their associated operation starting times will violate a subset of the disjunctive constraints. Unfortunately, it cannot be asserted that a candidate solution, whose starting times do not satisfy all the disjunctive Constraints with a given value of W, is inferior to a solution which is feasible with respect to all the disjunctive constraints. Hence, one must be careful not to set W too low, so as to eliminate potentially favorable candidate solutions. To solve this dilemma, one might take W a s 1.2 times a hypothetical makespan in which all intermediate and saleable demand is processed without considering product precedence requirements. Alternatively, if a good candidate solution is known, W may be selected to be 1.3 times the smallest value for which all disjunctive constraints are satisfied with these starting times. A summary of the model is presented in Figure 3. The number of variables and constraints are also presented in Figure 3. As noted in the demand constraint section, an equality constraint of the form Nil = X i l / B iis required for all i such that R(i) = 1. Since X i l / B iis a constant, Nil may be substituted throughout the formulation. Hence, the number of continuous variables and the number of constraints are both reduced by the cardinality of the set (i(R(i)= 1). The presence of the binary variables yikjk' requires that an integer programming algorithm be employed to search for the optimal solution. A branch-and-bound solution
=0
= 1 implies Yikjk' = 1
Yimjkr
T o t a l n u m b e r of V a r i a b l e s :
Yimlk!
The first condition states that if run k' of product j precedes run k of product i on a processor, then necessarily consequent runs of run k will follow run k'of product j , i.e., S j k ' I S i k 5 The second condition states that if run k'of product j follows run m of product i, then run k'of product j must also follow all runs of product i which precede run m. If Yikjk' = 1, then Yimjk' is free to take on the value of 0 or 1. Likewise, Y i m j k , = 0 does not dictate the value of Yikjk'. The addition of the constraint
si,.
Yimjk'
I Yikjk'
covers all cases. Case B: Yikjk' and Y m p l k ' for qim = 1 and p I k . The condition qim= 1 implies that there exists a directed path between nodes i and m. Then = 0 implies
Yikjk' Ympjk'
ympjk!
=0
Yikjk'
=1
= 1 implies
The addition of the constraint Ymplk'
I Yikjk'
covers all cases. Case C: Yikjk' and yipjplfor k this case yields
> p and k' < p'. Similarly,
I Yipjpl
Yikjk'
These constraints can be aggregated so that relationships among groups of variables can be expressed in a more concise form. For example, the two constraints Y1363
I Y1263
Y1363
Yll63
can be written as
+ Y1163 Though equivalent to the disaggregated constraints in the 2Y1363
Y1263
984
Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 4, 1986
Table I. Data for Test Problem la obj. function
demand prod production run (batches) 4 1 10 4 2 20 5 1 10 prod
due date wt 25 1 50 10 15 1 processing time/batch ( t ; ) 0.5 1.0 0.5 1.0 1.5
0
Figure 4. Precedence network of test problem 1. Products 4 and 5 are saleable.
“All f,, = 1 and u,, = 0. Two processors: VI = [1,2,4]; U2= [3,5].
C k’ are listed. Since the precedence network is acyclic, qrJ qJ,cannot equal 2, for this would imply ql, = q,l = 1.
integer sense, the aggregated constraint is not equivalent in the continuous sense; for instance Yl363
=
Yll63
+
For each (i,k,j,k‘)tuple, a pair of disjunctive constraints will be required in the formulation. In certain cases, a YakJk’ variable may be suppressed and its associated disjunctive constraints eliminated if a conflict between product i, run k and product j , run k’is unlikely, e.g., DJ, >> Djk’. In other words, the two operations are likely to occur a t opposite ends of the planning horizon. The optimal sequence/ schedule should be checked to ensure that the operation starting times are feasible with respect to all suppressed binary variables Yak#. The LINDO computer code (Schrage, 1981) and its associated branch-and-bound capabilities were chosen to solve the mixed integer programming problem. A thorough discussion of the branch-and-bound solution technique can be found in any standard mathematical programming text (e.g., Bradley et al., 1977).
y2
74
Y1263
=
satisfies the aggregated constraint but not the two disaggregated constraints. The aggregated constraints, however, may represent potential savings in computation time since the size of the model decreases as a result of aggregation. Further testing is required to determine whether the loss of information due to aggregation is significant and whether aggregation reduces CPU time. Before the section on the model formulation is concluded, it should be emphasized that unlike constraint sets I-V, the sets VI and VI1 are optional and have no bearing to the optimal solution. Therefore, it is not necessary to include in the formulation every constraint which falls in sets VI and VII. However, their inclusion may expedite the search for an optimal solution, improving the efficiency of the mixed integer approach.
Test Cases The attributes of the approach will be presented through two test cases. The first is a small problem which will provide the basis for discussion on different objective functions and multiple solutions. The second test case was presented by Suhami and Mah (1984), who studied a problem similar to the one addressed in this paper. The computer studies were performed on a DEC 2060 computer. Test Problem 1. The problem is presented in Figure 4 and Table I. It consists of two saleable products (4 and 5), three intermediates (1-3) and two processors. The problem was analyzed under four objective functions, namely, minimization of (a) mean weighted tardiness, (b) maximum tardiness, (c) mean completion time, and (d) makespan. The basic problem is characterized by 25 continuous and 8 binary variables and by a total of 41 constraints, 16 of which are disjunctive. The solution under objectives (b) and (d) above requires a slighly augmented set of constraints, as discussed in a previous section. Optimal solutions are shown in Figure 5. A summary of the solutions is presented in Table 11. These results demonstrate the existence of multiple solutions. For example, Table I1 shows that the optimal solution to the mean weighted tardiness objective is also optimal under the second and fourth objectives. If minimizing maximum tardiness were the objective of primary
Solution Methodology The binary variables used in constraint sets V and VI1 are identified in a preprocessing step. The inputs to this step are the number of products, N , the number of processors, M , the product sets for each processor k , U,, the number of production runs for each product, R(i),and the matrix of one level paths in the precedence network, P’. The matrix (I- P1)-lis used to determine the Q matrix. A conflict matrix C of size N’ x N’ is developed from the product sets uk. An element caJis equal to one if there exists a k such that i, jE uk. The conflict matrix is symmetric. To locate binary variables, the preprocessor scans through all (ij)pairs for which j > i. If crj = 0, then a conflict cannot arise and the next pair (ij)is checked. If cv = 1, then qv and qJaare summed. Should this sum be zero, then all runs of product i and all runs of product j are in potential conflict. Consequently, binary variables are listed for each (i,k,j,k’) tuple. If qbJ+ qJ1= 1, the computer code determines which q is nonzero. If qrj = 1 and product i has more that one production run (R(i)> l),then all tuples ( i , k i , k ’ ) for which k > k’are listed. If qJ,= 1 and RG) > 1, then all tuples ( i , k j , k ? for which k
Table 11. Summary of Optimal Solutions of Test Problem 1 (Cf. Figure 4) mean weighted T,, TA,+ Tu tardiness objective minimize mean weighted tardiness 30 25 0 28013 0 33013 minimize max tardiness 30 30 5 40 0 40513 minimize mean completion time minimize makespan 30 31013 30 25 +
1 0 1 1 0
+
max tardiness 30 30 40 30
mean completion time 14513 15013 13513 17513
makesuan 75 80 90
75
-
Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 4, 1986
Proconor -1
NI1.P
IN21.10I 15
NZ.20
25
IN4l.lOI I5
N42-20
55
7
NII~P.5
I
985
N11 Y(1 Y ( 1 = 7 . 5 * I 5 '15 41.25 I5 18.75 52.5
Nzz.,5
N21- Z.5
11.35
33.75
-2
0
15
33 75
0
15
fl5
Proconor
propnor
-2
-2 0
35
I5
45
55
60
80
-
Rocanm -1
46
I95
515
-2
FTw=$wr
0
-2
0
15
60
PO
propnor -2 0
25
30
45
55
I5
25
30
45
55
60
90
70
75
Bigure 5. Optimal solutions to test problem 1 under different criteria; t, = 1. (a, top) Minimize mean weighted tardiness. (b, top middle) Minimize maximum tardiness. (c, bottom middle) Minimize mean completion time. This is also the optimal solution when mean weighted tardiness is minimized and all weights are equal. (d, bottom) Minimize makeapan.
importance, then it would prove to be myopic to limit oneself to the optimal solution found by LINDO,since it is clear that multiple optimal solutions exist. The goal programming approach discussed earlier can be effective in this situation, because it allows these multiple optimal solutions to be reevaluated under secondary criteria. For all problems discussed in this section, a value of 90 was selected for W, the coefficient of the binary variables in the disjunctive constraints. The optimal solution to the mean weighted tardiness problem could have been obtained with a value of W set as low as 55. However, an attempt to minimize mean completion time with W set to 55 would cause LINDO to bypass the optimal solution since it is "infeasible"; that is, the disjunctive constraint s 5 1 - s32 + W(l - Y3251) N32t3 is not satisfied, 0 - 60 + 55(1- 0) $20(1/2) = 10 (cf. Figure 5c) Therefore, the value of W must be chosen carefully so as not to risk the loss of good candidate solutions for the sake of a reduction in CPU time. One of the merits of the model formulation is that production run quantities need not be specified by the decision maker. The benefits that may be derived from this innovation are best shown by example. The problem presented in Figure 4/Table I is modified by changing the processing time of product 4 ( t 4 )to 0.25 time unit. The schedule/sequence which minimizes mean weighted tardiness is shown in Figure 6a. The optimal objective function value is 48.75. This solution may be impractical because some of the optimal production quantities are noninteger. Rounding these quantities to the nearest integer leads to the solution shown in Figure 6b with weighted tardiness of 59.5. Had the decision maker restricted production quantities to the intuitively appealing values of 10 batches for all Nil associated with the first run
Figure 6. Effect of prespecification of the batch sizes on test problem 1; t4 = 0.25. (a, top) Optimal solution (batches allowed to be nonintegers); Tc = 48.75. (b, middle) Production quantities in (a) rounded to the closest integer; Tc = 59.5. (c, bottom) Optimal solution with Nil = 10 and Ni2 = 20, i = 1-4; Tc = 122.5.
m
W
n
w
0
Figure 7. Precedence network of test problem 2. All products except 5 and 6 are saleable.
of product 4 and 20 batches for all Ni, associated with product 4, run 2, then the formulation would yield the solution given in Figure 6c. Note that weighted tardiness shows a marked increase to 122.5 because these preset values of 10 and 20 batches do not allow for effective simultaneous production of produds on the two processors. Test Problem 2. Suhami and Mah (1984) have developed a heuristic procedure to tackle the N product, M processor sequencing/scheduling problem similar to the one being studied in this paper. The test problem shown in Figure 7 and Table I11 comprised of 12 products (10 of which are saleable) and 43 processors was presented in their paper. Since each product required only one production run and task durations were specified, the production quantity set of decision variables, Nik, is not necessary. The formulation consisted of 2 constraints of type I, 10 constraints of type IV, 32 constraints of type V (used to resolve the 16 potential conflicts), and 4 relational constraints of type VII. There were 22 continuous and 16 binary variables. Minimization of mean weighted tardiness with equal weights was chosen as the objective function. In solving the problem, an optimal solution was found after 20 partitions (solution of 40 subproblems) and 134 simplex iterations. No additional partitions were required to complete the branch-and-bound tree. The solution was
Be6
InU. Eng. Chem. Process Des. Dev., Vol. 25, No. 4, 1986
Table 111. Data for Test Problem 2 prod processor set processing time due date 180 720 1 37-40,44 160 720 2 26, 27, 41, 43 608 1440 3 28-35 4 1, 2, 4-6, 8, 10, 11 552 1440 360 5 3,1519 6 20, 24, 25, 31-34, 36, 40, 43 48 7 2, 4-6,9, 12-14, 19 576 1440 8 20, 22, 24, 41, 42 200 1440 360 1440 9 2,4-9,19 40 2160 10 37, 39,40, 42 48 2160 11 20-23 120 2160 12 15-19
48
208 384424
496
1104
7
4
i0 I 48
.
,
228
1
.,,
5
I 12 p$:yp
,',
,
1128
I
. . ,; ",, ,
'I ', ,,
I
,
i
,
, , , , _ ,,,, 1824
16081
,lo
,I
7
4
I
9
I12
I608
7
4
0
9
360
B]
KEI:
S l u t Tire
I12
I608
1
is also added to the formulation to preserve optimality under the tardiness objective. The problem is then resolved with the schedule in Figure 8a as a good initial solution. The resulting optimal solution is presented in Figure 8b. This solution differs from the single-objective solution in Figure 8a in the left shifting of product 11 so that its production begins a t time 1560. The minimum makespan is 1608. In this case the application of goal programming could continue, since there are multiple optimal solutions satisfying the two criteria. Suhami and Mah's solution is shown in Figure 8c. I t is optimal under both the tardiness and makespan objectives. The effort of the heuristic approach to obtain the schedule in Figure 8c was limited to the solution of three linear programs. Judging from the number of linear programs (LP), the Suhami and Mah approach appears far more efficient (3 vs. 40 LPs). However, this basis of comparison is misleading, since their linear programs correspond to reconstructed problems, requiring in general a full search. On the other hand, the proposed approach solves slightly modified linear programs from partition to partition, resulting in a small number of simplex iterations per linear program. The number of linear programs that need to be solved by the heuristic procedure is bounded by a quantity which increases linearly with the size of the problem. The corresponding number of linear programs required by the proposed procedure is bounded by &U, where n is the number of variables needed to resolve the potential conflicts. However, our computational experience suggests that only a small fraction of the theoretical upper bound of the linear programs is actually necessary. For example,
18241
9
768
"All products requires an additional 336 time units for postprocessing testing. Tardiness is defined as max {(F,+ 336) - Di, 01.
obtained after 4.87 CPU s and is presented in Figure 8a. All due dates are met, with the exception of product 4, which is tardy by 24 time units. The proposed mathematical programming approach was also implemented on an IBM personal computer using PC-LINDO. The full model required 1 min and 50 s to obtain the solution in Figure 8a. The solution in Figure 8a indicates that since the "limiting" sequence is 4, 7,9, 12, multiple solutions exist which simply result in rearrangement of the rest of the schedule. In such a case, a decision maker might reevaluate the solutions based on a secondary criterion. Preemptive priorty goal programming may be applied using as a secondary criterion the minimization of makespan. The program is modified by changing the objective function and adding the associated set of constraints discussed earlier in the objective function section. The constraint CTi+I24 i€A
624
--+IO
50
+Piaimh T i w
Id'* Procersor
Figure 8. Optimal solutions to test problem 2. (a, top) Proposed mathematical programming approach; 7'+ = T7+= 24. (b, middle) Heuristic approach; T+= T7+= 24. (c, bottom) Goal programming applied to the proposed approach; the solution satisfied both criteria: minimum tardiness and minimum makespan.
although in test problem 2 there are over 130000 combinations ( n = 16) of binary variables that could have been searched, the proposed procedure searches over 40 combinations only. The primary differences between our method and that of Suhami and Mah's are 3-fold. First, in our approach production quantities are variables, not parameters of the model. Second, due dates for intermediate products need not be specified. Finally the proposed approach is an exact method and, as such, guarantees an optimal solution if the computer run is allowed to run to completion. Computational experience has shown that LINDO will find a good solution early in its search (regardless of the size of the problem), though the entire computer run may take a lengthy period of time either in fine tuning this good candidate solution (if suboptimal) in order to drive it to optimality or in completing the branch-and-bound tree to verify that this good solution is optimal. When computer runs appear to be taking excessive time, we have found it beneficial to stop the run, review the branch-and-bound search, and restart the run with certain binary variables fixed so that different areas of the branch-and-bound tree are forced to be searched. Extensive computer experimentation is necessary to establish a relationship between the size of the problem
Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 4, 1986
and the effort required by the proposed approach. Potential measures of problem structure size (complexity) include (1) number of products, (2) percentage of saleable products, (3) average number of orders per saleable product, (4) number of processors, (5) average number of products in a processor set, (6) average number of processors required by a product, (7) processor congestion, (8) variance of activity durations, (9) number of potential conflicts, (10) network connectivity (density), (11) distribution of demand and placement of due dates over the planning horizon. The importance of these measures is currently under investigation, and the results will be presented in forthcoming publications. These measures can also be used to determine for what problem structures heuristic techniques and exact methods are best suited. For example, if in test problem 2 the due dates of products 4, 7, 9, and 12 are changed to 720,720,1080, and 1440, respectively, and no postprocessing of saleable products is required, the mathematical programming approach would alter the schedule in Figure 8a by shifting the start times of these products to S, = 552, S, = 912, S12= 1032, and S7 = 1608. The total tardiness is now 888 which can be attributed to lateness in the production of product 7. The heuristic, however, would not alter its schedule in Figure 8c. The reasoning is as follows: The initialization phase of the algorithm orders products based on the period in which the product is due and for those products due in the same period, priority is given to products in increasing order of their processing times. Hence, the process ordering is initially 4,7,9, and 12, assuming a period of 360 time units. This ordering can be changed in the heuristic's improvement phase only if an idle processor is identified; that is, a gap in the schedule is located into which a product can be introduced. Because processors 2, 4, 5, 8, and 19 are utilized consistently throughout the planning horizon, no gaps will result in the schedule which will allow for ordering transpositions among products 4, 7, 9, and 12. Consequently, this initial process order remains, leading to a suboptimal solution of total tardiness of 984 (T7+= 408, T9+= 408, T12+ = 168). Therefore, if processor congestion is likely to occur due to heavy demand for certain processors over the planning horizon, the initial process ordering established by the heuristic is crucial because the lack of gaps in the schedule will reduce the possibilities for changes in process orderings during the improvement phase. The mathematical programming approach has the advantage of being able to search through all promising process orderings on a congested processor. Conclusions and Significance A mathematical programming formulation is proposed which simultaneously considers the scheduling and sequencing short-range planning activities in a production control system for a multipurpose plant with batch operations. The mixed (binary) integer program attempts to determine the number of batches of products to produce in each production run and the prrocess orderings on all plant equipment units. The schedule/sequence must be consistent with physical precursor relationships among the products, and quantities produced must satisfy customer demand requirements. In addition, the model resolves potential conflicts when two or more products require the same processing unit(s). The formulation offers many advantages to the user. Production run quantities are variables, not parameters of the model. Due dates for intermediate products need not be specified. There is a variety of objective functions from which to choose. Commercial mixed integer pro-
087
gramming software can be used to solve for the optimal schedule/sequence. A good suboptimal or an optimal solution is generally attainable for practically sized problems with reasonable computational effort. Finally, sensitivity analysis is readily implemented by simply altering parameters in the model formulation or adding/deleting constraints. Alternatively, the recent methods proposed by Schrage and Wolsey (1985) or Ohtake and Nishida (1985) can be used for parametric analysis. Future directions for research include devising a better means to handle general structure precedence networks and extending the formulation to additionally consider material, labor, and/or storage requirements and restrictions. Work is currently under way to determine which problem attributes have the greatest effect on computational effort. Nomenclature A = set of indexes corresponding to saleable products Bi = batch size for product i Bin = number of binary variables C = conflict matrix cij = element of conflict matrix C Dik = due date for kth order of saleable product i Fik = time of completion of production of saleable product i, production run k fi, = number of units of product i required to produce one unit of immediate successor product j G k = idle time on processor k Iio = initial inventory for product i l i k = number of batches of product i that are produced in excess of demand by the end of production run k M = number of processors N = actual number of products N' = total number of nodes in the precedence network, N'>N, the strict inequality holds whenever portions of the precedence network must be copied to accommodate either (a) a product which is both intermediate and saleable or (b) a product with multiple successors (Figure 2) NBIC,( = number of arcs associated with saleable product i Nik = number of batches of product i produced (consecutively) in production run k (decision variable) P" = precedence matrix for n-level paths P* = structural matrix of the network pijn = element of matrix P" pij* = element of matrix P* Q = incidence matrix of the network qij = element of matrix Q R(i) = number of customer orders for saleable product i; maximum number of production runs for product i S i k = start time of production of the kth run of product i (decision variable) Tik+ = tardiness of kth order of saleable product i (decision variable) = max [0, (Sik+ Nikti) - Dik] ti = processing time per batch of product i uk = the set of products which require processor k u;j = time required to postprocess a batch of product i to ready it for product j production W = coefficient of binary variables in the disjunctive constraints Wik = penalty cost (relative) per unit that product i, run k is tardy Xik = size of the kth customer order for saleable product i (expressed in the same units as Bi) Y' = makespan yikj = binary variable associated with the pair: product ilrun R - product j/run p 2 = maximum tardiness Literature Cited Afentakls, P.; Gavish, 6. Working Paper No. QM-8318,1985; University of Rochester, Graduate School of Management. Bradley, S. P.; Hax, A. C.; Magnanti, T. L. Applied Mafhemaflcal Programming; Addlson Wesley: Reading, MA; 1977. Graves, S. Oper. Res. lD81, 29(4), 646.
988
Ind. Eng. Chem. Process Des. Dev. 1986, 25, 988-996
Greenberg, H. H. Oper. Res. 1988, 76(2), 353. Grossmann. I. E.; argent. R. W. H. I n d . Eng. Chem. process Des. Dev. 1979, 78(2), 343. Knopf, F. C.; Okos, M. R.; G. V. Reklattls, 0. V. Ind. Eng. Chem. Process D e s . D e v . 1982, 27(1), 79. M a M i , A.; Rlppin, D. W. T. Chem. Eng. Rug. 1980, 76(4). 37. Ohtake, Y.; NisMda, N. Oper. R e s . Lett. 1985, 4(1), 41. Overturf. E. W.: Reklaltk. G. V.: Woods. J. M Ind, E m " . Chem. Rmess Des. &.'1978,'77(2), 166. ' Prabhakar, T. Manage. Sci. 1974, 27(1), 34. Schrage, L. Users Manual for LINDO, ScientlHc: New York, 1981 Schrage, L.; Woisey, L. Oper. Res. 1985, 33(5), 1008.
Sparrow, R. E.; Forder, G. J.; Rippin, D. W. T. Ind. Eng. Chem. Process D e s . D e v . 1975, 14(3), 197. Suhemi, I.; Mah, R. S. H. Ind. Eng. Chem. procesS Des. D e v . 1982, 27(1). 94. Suhaml, I.; Mah, R. S. H. Proceedings of the Second FOCAFD Conference; Westerberg, A. W., Chien, H. H., Eds.; American Institute of Chemical Englneers: New York, 1984.
Received for review March 12, 1985 Revised manuscript received March 17, 1986 Accepted May 20, 1986
Experimental and Mudel9ng Studks in Fixed-Bed Char Gasification Amltava Bhattacharya, Lyle Salam, M. P. Dudukovlb, and Babu Joseph" Chmicai Engineering Department, Washington Universlfy, St. Louis, Mlssowi 63 730
A laboratory benchscale reactor for the gasification of char and coal using steam and air is described. This setup is designed to generate data that can be used to valldate proposed gasifier models. A previously proposed gasifier model was adapted to simulate thls sasifler. Experhnental results on the gasification of Wyodac char are compared with model predictions. The model shows good agreement with experimental data in the initial stages of gasification. The discrepancy between the model and experimental data at later stages of gasification suggest that better modeling of the kinetics of chargas reactions is needed.
Fixed-bed coal gasification is a simple and reliable way of converting coal to gaseous fuels. Among a number of gasifiers of this type, the Lurgi gasifier is most commonly used commercially. Increased use of fixed-bed gasification led to model development. Yoon et al. (1978) developed a steady-state, homogenous model for the Lurgi gasifier. Radial variations were accounted for approximately by consideration of an adiabatic core and a boundary layer next to the wall. Cho (1980) and Desai (1975) developed one-dimensionalheterogeneous models for the steady-state case. Kim (1981) extended Cho's formulation to include the dynamics of a fixed-bed gasifier when it went from one steady state to another because of a disturbance in one of the inputs. Yu (1981) included radial effects in Yoon's formulation and also incorporated the dynamics. There have been very little attempts a t validating the above models primarily because of a lack of appropriate experimental data. A notable difference is a paper by Fbm and Thorsness (1978), who have modeled one-dimensionalflow in a permeable packed bed. They also report experimental results on product gas compositions and reaction zone propagation rates, though not on temperature profiles. Most comparisons so far have been with exit gas compositions of plant-scale gasifiers operated a t steady state. However, this is insufficient to establish the validity of the model since exit gas compositions are well defined by the constraints of overall steady-state mass balances and equilibrium relations (Yoon et al., 1978). Thus, it is possible to match these compositions without detailed modeling of the phenomena taking place in the interior of the gasifier. Yet, it is this interior behavior that is of interest from a design and operation standpoint. Some attempts have been made to obtain detailed experimental information on gasification and combustion in a fixed bed of coal. These experiments, however, have been done to support mainly the modeling of the underground
* Author to whom correspondence should be addressed. 0198-4305/86l1125-0908$01.50/0
coal gasification processes. Westbrook (1976) reported the construction of a vessel to investigate in situ combustion of Texas lignite. Successful experimental results were difficult to obtain because of the high temperatures involved. Work by Lu (1980) involved the construction of a similar tube for combustion experiments. The results were compared with a one-dimensional mathematical model. Data on large-scale gasifiers are either not easily measured for lack of good instrumentation or not available for proprietary reasons. On the other hand it is difficult to construct and operate a small-scale fixed-bed gasifier with continuous feed. An alternative here is to construct a "true" fixed-bed reactor without a continuous feed of coal. The operation in this case will be batch and hence dynamic in nature. Such an approach is followed in this work. The primary objective of the work was to develop an experimental database for the validation of proposed gasifier models. Preliminary results on the gasification of Wyodac char are provided in this article. The experimental data are compared with results obtained from a proposed model. Due to the cumulative nature of modeling errors in this transient test, we believe that this setup provides an excellent vehicle for validating existing and new proposed models for fixed-bed gasification.
Experimental Setup and Experimental Results This section briefly describes the process of char gasification as carried out in a bench-scale fixed-bed reactor which was designed and operated by Salam (1983). The overall process is similar to commercial-scale systems requiring similar unit operations. A schematic summary of the process is given in Figure 1. The solid reactant in this process is devolatilized Wyodac coal (char). The devolatilization process is carried out separately, and the devolatilized char is loaded into the gasification unit. Standard analyses (proximate and ultimate) of the coal and the char were performed, and these are presented in Tables I and 11. Char was used instead 0 1986 American Chemical Society
