Using Detailed Scheduling To Obtain Realistic Operating Policies for a

Ind. Eng. Chem. Res. 1997, 36, 1717-1726

1717

Using Detailed Scheduling To Obtain Realistic Operating Policies for a Batch Processing Facility Matthew H. Bassett* DowElanco, 9330 Zionsville Road, Indianapolis, Indiana 46268

Joseph F. Pekny and Gintaras V. Reklaitis School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907

With many industries changing to a just-in-time manufacturing environment, the ability to reliably meet due dates and/or determine delivery/promise dates is gaining importance. In order to achieve this, the ability to determine operating policies based on realistic production plans is paramount. For these policies to be reasonable, the uncertain nature of the process must be addressed. These uncertainties include, but are not limited to, processing time fluctuations, equipment reliability/availability, process yields, demands, and manpower changes. In this paper we present a framework for addressing uncertainty by means of Monte Carlo sampling. This framework is not limited to a specific model, but results obtained show that the need of including sufficient detail is best satisfied by a scheduling model. A number of stopping criteria are derived, and the framework is utilized to obtain operating policies for two industrially-based examples. 1. Introduction To minimize inventory costs, many chemical processing industries are turning to just-in-time manufacturing. This change in philosophy requires planners and schedulers to delay production as long as possible such that intermediates and final products are stored for the minimum amount of time. To reliably meet due dates or predict delivery/promise dates in this environment, realistic operating policies are a necessity. In order to determine these policies, realistic production planning is necessary. For production planning to be realistic, it needs to take into account a number of uncertainties which can effect throughput. These can include, but are not limited to, processing time uncertainties, equipment reliability/availability, process yields, demands, and manpower fluctuations. Much of the production planning literature has focused on the use of aggregate production planning (APP) models to determine production quotas (see the review by McKay et al., 1995). In general, APP models determine lot sizes for all products within each aggregate interval and are therefore formulated as linear programs (LPs). Some researchers have looked at including a number of uncertain parameters within APP; however, little work has been done to include these parameters for detailed resource constrained scheduling problems (RCSP) associated with the chemical engineering field. RCSPs, unlike APP models, make discrete decisions concerning task-unit assignments to meet production targets, thus necessitating the use of a mixed-integer linear programming (MILP) formulation. This discrete nature leads to difficulty in solving RCSPs deterministically for industrial-scale problems. However, if a RCSP solution can be obtained, it is apt to be much more realistic than the APP model solution for a number of reasons. First, no aggregating assumptions are made, which can lead to overestimating the production capacity. For example, a change in processing times can have far-reaching effects due to precedence relationships within recipes. The RCSP takes these relationships into account, while the APP model misses these interactions due to the aggregation of time. Second, a unit being unavailable 20% of the time simply S0888-5885(96)00470-8 CCC: $14.00

reduces the total time available on the unit in the APP model, whereas within the RCSP it is important to know exactly when the downtimes occur. This is due to the fact that an equipment failure during a critical utilization period has a much greater effect on the production plan than if the unit is underutilized during this time. In this paper we present a framework for including uncertain parameters within a general APP or RCSP model. This is achieved by using Monte Carlo sampling to generate random instances, determining a schedule/ production plan for each instance, and generating distributions of aggregate properties to infer a number of operating policies. These aggregate properties may be the time necessary to meet demands, the success/ failure to meet a due date, maximum intermediate inventory levels over the horizon, etc. It should be stressed that this framework does not determine a specific robust schedule but determines operating policies which are robust. The policies determined using these distributions could include realistic production lead times, possible maintenance protocols, and strategic inventory levels. We can determine these operating policies by addressing the following issues: (a) determining the probability of meeting given due dates for specified demands with a given confidence or (b) determining a promise/delivery date for a new demand such that it is missed a given percentage of the time with a given confidence. We will compare the results obtained using a RCSP model to those achieved when a APP model is used within the same framework. The progression of the paper is as follows. Section 2 will present a number of model formulations which are used throughout to remainder of the paper. Section 3 will detail the Monte Carlo framework, with specific emphasis on possible stopping criteria and data analysis methods. Section 4 will discuss the implications of assumptions made in the framework as it is applied in practice versus theory. Section 5 will present a number of industrially-based examples followed by some concluding remarks in section 6. 2. Formulations 2.1. Uniform Discretization Model Formulation. Uniform discretization models (UDMs) for scheduling © 1997 American Chemical Society

1718 Ind. Eng. Chem. Res., Vol. 36, No. 5, 1997

problems were first proposed as early as 1959 (see Bowman, 1959; Manne, 1960). In terms of industrialscale problems, very little was done with UDMs at the time, but with the increased capabilities of modern computers, interest in these models has reemerged. Within the chemical engineering community, UDMs and possible solution methods were investigated by Kondili et al. (1988). Since then, significant work utilizing UDMs has been published in the chemical engineering literature due to the general applicability of the model (Sahinidis and Grossmann, 1991; Kondili et al., 1993; Pekny and Zentner, 1993; Shah et al., 1993; Pantelides, 1994; Zentner et al., 1994). The UDM formulation discretizes the planning horizon into uniform time intervals based on the greatest common denominator of the processing times. In this way, all the dynamics of the scheduling problem can be captured in the fewest intervals. Variables are generated within each interval to model the inventory levels (asjt′), task-unit assignments (wijt′), production levels (bijt′), and resource purchases (qbst′). These variables are used to form constraint families which define feasible schedules. One important constraint family limits a unit to processing only one task at a time:

∑ i∈I

∑

pijtWijt e Ht ∑ i∈I

s

i

∑ ∑ i∈I j∈J

s

∀s ∈ S

∀t ∈ T (5)

min Vmax ij Wijt g Bijt g Vij Wijt ∀i ∈ I

∀j ∈ Ji

∀t ∈ T (6)

i

where

∀s ∈ S

Asjt ) asjt′ Bijt )

bij′ ∑ t′∈t

A second constraint in the UDM formulation maintains a material balance from one time interval to the next for every resource. This constraint ensures that the amount of material in inventory at the end of the interval is equal to the amount in inventory at the beginning of the interval plus material additions due to purchases and completed tasks minus material consumed due to demands and tasks started during the interval:

QBst )

∑

j∈Js

∑ ∑ i∈I j∈J

bij′ fconsume + dst′ ) ijs

j∈Js

i

∑ ∑bij′-p i∈I j∈J

∑

asjt′-1 + qbst′ +

ijt′

fcreate ijs

∀s ∈ S

∀t′ ∈ T′ (2)

i

The final important UDM constraint, considered for purposes of this paper, limits the production level to be between a given upper and lower limit when a task is assigned to a unit: min Vmax ij wijt′ g bijt′ g Vij wijt′ ∀i ∈ I

∀j ∈ Ji

∀ t′ ∈ T′ (3)

The major problem with this model is that the greater the number of intervals, tasks, units, and/or resources, the larger the dimensionality of the formulation generated. Most importantly, the number of binary assignment variables increases, leading to greater combinatorial complexity. It is the inability to effectively address the combinatorial complexity of the UDM formulation which is its greatest drawback, limiting its usefulness to smaller problems unless highly customized algorithms are employed (Bassett et al., 1996a). However, the combinatorial complexity is what gives reliable production estimates; thus, any solution methodology must address it in some manner. 2.2. Aggregate Planning Problem Model Formulation. From the UDM formulation, a more compact but less exact formulation can be easily constructed

(4)

Bijt-pijtfcreate ijs

x|(pijt′+x>t′ and xet′)

asjt′ +

∀t ∈ T

Asjt + ∑ ∑ Bijtfconsume + Dst ) ∑ Asjt-1 + QBst + ∑ ijs j∈J i∈I j∈J j∈J

Dst )

wijx e 1

∀j ∈ J

j

∀t′ ∈ T′ (1)

j

∀j ∈ J

(Bassett et al., 1996a). This formulation is referred to as an aggregate formulation because it is generated by summing over all the UDM variables and constraints for a predetermined set of discrete time intervals (in other words, aggregating the UDM formulation). This yields

Wijt )

∀j ∈ Js ∀i ∈ I

dst′ ∑ t′∈t qbst′ ∑ t′∈t

wijt′ ∑ t′∈t Ht )

1 ∑ t′∈t

∀ t′ ) t

∀j ∈ Ji

∀s ∈ S ∀s ∈ S

∀i ∈ I

∀j ∈ Ji

∀t ∈ T

(7) (8) (9) (10) (11) (12)

Severe assumptions are made in constructing the APP model formulation. The first is that any demands which are due during the time covered by an aggregate interval are assumed to be due at the end of the interval. Since this may artificially extend due dates which do not coincide with aggregate interval end points, it is important to ensure that aggregation is done systematically to reduce this effect. In a similar manner, arrival and/or purchase of resources occurs only at the beginning of aggregate intervals. Finally, due to the loss of sequencing characteristics, changeover costs and/or times cannot be modeled within the aggregate formulation. This can lead to vastly overestimating production capacity (underestimating production costs) when changeover times (costs) play a significant role in processing. This overestimation is acceptable when planning for highly dedicated facilities or simple process recipes since these effects are typically minimal. However, when equipment and/or task sharing play a large role in the facility, a more detailed approach is needed. This consideration led to the approach presented in the next section, in which many of the details of the UDM formulation are obtained at the reduced expense of the APP model formulation. 2.3. Reverse Rolling Window Approach. In Bassett et al. (1996b), a rolling window approach to solving large-scale scheduling problems was presented. This formulation has the desirable property of merging both the UDM and APP model formulations into a single hybrid formulation. As shown in Figure 1, the reverse rolling window approach (reverse RWA) starts with a

Ind. Eng. Chem. Res., Vol. 36, No. 5, 1997 1719

Figure 2. Monte Carlo framework.

Figure 1. Reverse rolling window with crossover variables and LP region.

detailed UDM formulation window at the end of the horizon and the APP model formulation everywhere else. This is the opposite of other rolling horizon/ window approaches but was shown to reduce both solution time and infeasibilities. As the solution progresses, the window shifts, leaving behind a fixed locally feasible solution. In this way, a solution to the global problem is obtained from a sequence of local solutions. In order to reduce the possibility of obtaining a locally infeasible solution to a globally feasible problem, an LP crossover region is introduced to mesh the aggregate and detailed regions. Within this region, the UDM formulation is generated but integrality is not enforced. Due to the lack of integrality, a feasible schedule is not obtained for the crossover region, but this region is a more faithful representation than the aggregate region since process causality is taken into account. This allows a smooth transition between detailed and aggregate regions and was shown to improve solution algorithm effectiveness. 3. Monte Carlo Framework The inclusion of uncertain processing times and/or equipment availability within the formulations presented can be difficult. Most recently, Ierapetritou and Pistikopoulos (1996) theoretically addressed process time and demand uncertainty through the use of scenarios for both detailed scheduling and lot-sizing constraints within a plant design formulation. A scenario contains a discrete value for all uncertain variables within a given time interval and its associated probability. Thus, if there are 10 uncertain variables at two levels each, a total of 1024 (210) scenarios will be needed. The scheduling constraints and variables corresponding to each of these scenarios are combined into a single composit problem formulation whose objective function includes the sum of the performance terms corresponding to each scenario weighted by its probability of occurrence. Obviously, this quickly leads to a very large formulation for even small problems, as is bore out by the size of the example problems presented in which only demand uncertainty for at most two products was considered. Additionally, determining what scenarios to generate can be a complex problem in and of itself, especially if the data given is continuous. To overcome the problem of both scenario generation and formulation intractability, a Monte Carlo framework can be utilized. Within this framework, an itera-

tion consists of randomly sampling each uncertain parameter from its given distribution and then solving the associated problem. Then from the data collected for a number of iterations, trends and distributions can be determined and statistically analyzed. Lee and Malone (1995) used Monte Carlo sampling and simulated annealing for production sequencing under uncertainties in demand, due date, processing time, and production cost. Thompson et al. (1993) performed a Monte Carlo simulation for a lot-sizing problem to test a number of operating policies with uncertainty in demands, equipment availability, and production costs. Within the Monte Carlo framework there are four major considerations, as shown in Figure 2: parameter sampling scheme, solution methodology, stopping criteria, and data analysis. In order to fully understand the framework, each of these will be discussed in detail. 3.1. Sampling. It is important that the uncertain parameters are randomly sampled from their probability distributions. However, whole-sale sampling of all parameters can lead to a computationally intractable problem size. Thus, some suitable assumptions are needed for more efficient sampling of processing times and equipment downtimes, which are the uncertainties being dealt with in this paper. 3.1.1. Processing Times. In the interest of reasonable formulation sizes, only integer processing times are considered. This restriction is necessary due to the use of the greatest common factor (GCF) for the discretization of time in the UDM formulation. Allowing arbitrary real values would require discretizations on the order of the number of significant digits, thus enormously increasing problem size. For problems with short processing times, this limitation could force drastic deviations from the average to occur. For example, if the average processing time is 2, the minimum symmetric deviation is (1, which is a (50% change. The deviation in this case could be reduced to (25% by using a GCF of 1/2 at the expense of doubling the formulation size. For the examples presented later, we found that the use of a GCF of 1 was not detrimental to the results obtained; however, it is ultimately up to the modeler to determine the discretization which is acceptable. In section 4, we will discuss further the implications of this assumption. 3.1.2. Equipment Downtime. In the absence of detailed failure data, specific assumptions are needed in order to determine the downtime profiles for each equipment item. For the present purposes, failures for each piece of equipment are assumed to follow an exponential distribution and once a failure occurs the equipment downtimes are assumed to follow a lognormal distribution. The exponential distribution is


used because it has been shown to mimic failur patterns for numerous systems (Lipson and Sheth, 1973). The log-normal distribution allows downtimes to only be positive, thus removing the possibility of negative downtimes. In order to use these distributions, an equipment “half-life” and mean repair rate are needed. These values are called the mean time to failure (mttf) and the mean time to repair (mttr), respectively. Depending on the values chosen, different downtime profiles emerge for the same overall downtime. Ideally, actual equipment reliability data can be obtained to help determine the values for these parameters before runs are made. Within the RCSP formulation, the method for including equipment downtimes is to set the assignment variables (wijt′) equal to zero for all tasks i when unit j is not available. For cases when the APP model formulation is used, either the downtimes calculated can be aggregated over the specific aggregate regions or the overall percent downtime can be evenly spread over all aggregate intervals. The former leads to slightly more accurate solutions depending on the level of aggregation used. The latter, on the other hand, has the advantage of not requiring any assumptions on the half-life and repair rate for the equipment. In either case, the downtime is modeled by reducing the right-hand side of eq 4 by the amount of time which the equipment is unavailable. For this study a priori knowledge of when equipment will be down is assumed. This means that the schedule generated is able to smartly reroute production to avoid unavailable equipment. In section 4 we will discuss the implications of this assumption on both theory and practice. 3.2. Solution Method. Depending on the size of the facility as well as the scope of the solution required, one of the three formulations presented could be used. For problems with minimal task precedence relationships, the APP model formulation may be sufficient. If the horizon or recipe network are small enough, using just the detailed UDM formulation may be tractable. However, for the vast majority of industrially relevant planning problems, the use of the reverse RWA is necessary. With the reverse RWA, task precedence relationships as well as practically sized problems can be included. In order to obtain realistic operating policies, it is desirable to minimize lead times required to meet demands, thereby maximizing the production rate. To minimize lead time, an objective function is used which, for the case of the reverse RWA, strictly monotonically increasingly penalizes production the further from the demand that it occurs:

Min

Cijt′bijt′) ∑ ∑ (∑ CijtBijt + ∑ i∈I j∈J t∈T t′∈t

(13)

i

with the requirement that the modeler predefines the cost parameters such that they satisfy the following conditions:

Cijt1 > Cijt2

∀t1 < t2

(14)

Cijt′1 > Cijt′2

∀t′1 < t′2

(15)

Cijt′1 > Cijt2

∀t′1 < (t′2 ∈ t2)

(16)

Cijt1 > Cijt′2

∀(t′1 ∈ t1) < t′2

(17)

This objective forces as much production close to the demand as possible. Then the “makespan” is given by the production finish date minus the start time of the first production batch. This objective does not necessarily assure the absolute minimum makespan is obtained, but in practice it has given reasonable results. Of course, the modeler is free to modify the objective function based on the problem data. For example, it may be useful to utilize a declining future value approach applied to all time-dependent parameters as the forcing function. In order to allow the possibility that production may take longer than the original horizon, negative time can be utilized so that infeasible solutions are not obtained. 3.3. Stopping Criteria and Data Analysis. In any Monte Carlo calculation procedure, the issues of when to stop the sampling/solution process and how to analyze the solution information which is obtained are intimately connected. In general, the nature, quality, and, more specifically, confidence level of the information which will be extracted from the solution data will dictate the number of required instances of the sampling/ solution process. Since the output for each iteration is the makespan for that instance, the stopping criteria(s) must obviously utilize this information. There are a number of different stopping criteria which ca be used depending of the answer being sought. These criteria can utilize the mean standard deviation or confidence interval for the makespan, the probability of violating (meeting) a due date, or even just the number of iterations. For the case where the due date is given, the probability of due date violation is the obvious metric to utilize since this is what we wish to determine in the end. The upper limit for the probability of due date violation (which we will refer to as the failure probability, Pu) can be calculated by utilizing a simple equation derived from the binomial distribution (Pieruschka, 1963):

Pu )

1 N-f 1+ (f + 1)F(R,2(f + 1),2(N-f))

(18)

The use of this equation does not require any knowledge about the distribution which due date violations follow. All that is needed is the number of iterations (N) performed so far, the number of due date violations (f) (makespan > due date), and the confidence level (R) with which the probability is required. From these values, the F-distribution is calculated and finally the failure probability is determined. The failure probability is a constant for a given system and scheduling/ planning approach. The use of eq 18 allows this value to be statistically estimated with a selected confidence for a given sample. The larger the sample, the more likely that the value calculated is the value for the population. A possible stopping criteria could be when the calculated failure probability is less than some preselected tolerance level. However, depending on the due date, the actual failure probability may be a value higher than the tolerance level. Figure 3 shows both possibilities. If a due date violation rate of no more than 10% is acceptable, the algorithm will stop for both the 95 and 99% confidence levels. However, if a due date violation rate of no more than 5% is required, the algorithm will not stop even after 1000 iterations since the failure


the makespan. These values are easily calculated by eqs 19-21. Unless we have prior knowledge of the Ni

∑ Valuek

Meani )

Std Devi )

Figure 3. Calculated failure probability.

Figure 4. Change in calculated failure probability with calculated bounds for 99 and 95% confidence levels.

probability for the given system is ∼6% for both confidence levels. Instead of continuing indefinitely, it is desirable to stop when the calculated failure probability has plateaued. The possible change in the calculated failure probability can be bounded by using the data collected so far and postulating future occurrences. Consider performing one additional iteration. The solution of this iteration will either meet its due date or not. A due date violation will increase the calculated failure probability by the greatest amount possible, while meeting the due date will reduce the calculated failure probability by the greatest amount. Taking the difference between these future probabilities and the current probability gives a measure of the maximum and minimum slope at the current iteration. Based on this information, a possible stopping criteria is when these slopes are within a chosen tolerance level. From Figure 4 it can be seen that the bounds calculated are very tight since the result of the next iteration is known a priori to be either a met or missed due date. For the data presented, if a tolerance of (0.5% is acceptable, the algorithm will stop after 248 iterations, giving a failure rate of 7.96% for a 99% confidence level, or after 236 iterations, giving a failure rate of 6.56% for a 95% confidence level. The failure probability cannot be used as a stopping criteria if the due date is not fixed, since a due date is needed to determine failures. However, stopping criteria can be constructed based on the mean, the standard deviation, and/or the confidence interval for

x

k)1

Ni

(19)

Ni

(Valuek - Meani)2 ∑ k)1

Conf Inti ) (

Ni - 1 Std DevitR/2;Ni-1

xNi

(20)

(21)

facility, the mean and standard deviation are unknown quantities, but what we do know is that these values will converge to constants if enough data are collected. The confidence interval, on the other hand, always approaches zero with enough data points since it is a function of 1/xNi. So a possible stopping criteria could be when the confidence interval is within a preselected tolerance. However, like the case for the failure probability, there comes a point of diminishing returns when there is no need to continue (maximum and minimum slopes are within tolerance). In order to determine the extreme (bounding) slopes, assumptions must be made concerning the maximum and minimum makespans which may be encountered during the solution of the problem. For the case where only processing times are uncertain this can be accomplished by solving two problems, one where all processing times are maximum and one where all processing times are minimum. This yields exact bounds for the maximum and minimum possible makespans, respectively. However, when equipment reliability is uncertain, this approach underestimates the extreme makespans, leading to incorrect bounds for the slopes. One way to improve the makespan estimates is to utilize the mean and standard deviation for the data already collected to postulate future extreme values. By assuming that makespans can be fit by a normal distribution (the validity of which will be tested in the examples presented later), we can calculate the number of standard deviations needed (zR/2) to cover a preselected percentage of the distribution. For example, by using (1.96 standard deviations, we cover 95% of a normal distribution described by the data collected so far. Using these maximum and minimum makespans, bounds can be calculated for the maximum and minimum change in mean, standard deviation, and confidence interval. The relative change in the mean for an additional iteration is given by:

% ∆Mean ≡

Meani+1 Valuei+1 - Meani (22) -1) Meani Meani(Ni + 1)

So the maximum % ∆Mean is achieved when Valuei+1 is maximized and the minimum % ∆Mean when Valuei+1 is minimized. If Valuei+1 is assumed to equal Meani ( zR/2 × Std Devi, then eq 22 can be further simplified to:


% ∆Mean ) (

zR/2 × Std Devi Meani(Ni + 1)

(23)

Continuing with the standard deviation, the relative change for an additional iteration is given by:

% ∆Std Dev ≡

x

Std i+1 -1) Std Devi

1-

(Valuei+1 - Meani)2 1 + - 1 (24) Ni N × Std Dev 2 i

i

The minimum % ∆Std Dev is achieved when Valuei+1 is equal to Mean, making eq 24

min % ∆Std Dev )

x

1 -1 Ni

1-

(25)

Figure 5. Percent change in mean with calculated bounds for 99 and 95% confidence levels.

The maximum % ∆Std Dev occurs when Valuei+1 is as far from Meani as possible. If Valuei+1 is assumed to equal Meani ( zR/2 × Std Devi, then eq 24 can be simplified to:

max % ∆Std Dev )

x

zR/22 - 1 1+ -1 Ni

(26)

The confidence interval is related to the standard deviation by eq 26. In general, the % ∆Conf Int is calculated by:

% ∆Conf Int ≡

Conf Inti+1 -1) Conf Inti

tR/2;Ni Std Devi+1 tR/2;Ni Std Devi Std Devi tR/2;Ni-1

x

Ni - 1 (27) Ni + 1

Figure 6. Percent change in standard deviation with calculated bounds for 99 and 95% confidence levels.

Utilizing eqs 25 and 26, this can be simplified to:

min % ∆Conf Int )

tR/2;Ni tR/2;Ni-1

x

Ni - 1 - 1 (28) Ni + 1

and

max % ∆Conf Int )

tR/2;Ni tR/2;Ni-1

x

Ni + zR/22 - 1 -1 Ni + 1 (29)

Figures 5-7 show how the bounds progress for each of the previous criteria for an arbitrary run of 1000 iterations. The bounds for the % ∆Mean are symmetric and vary with zR/2 as expected from eq 23. The standard deviation and confidence interval lower bounds are very tight since they are only a function of the number of iterations and the confidence level. The upper bounds can be tightened by reducing the number of standard deviations covered by the extreme makespans, but this can lead to even more values falling outside these bounds due to incomplete coverage. The greater the coverage, the fewer the number of points outside the predicted range, but the longer the time needed to achieve a given tolerance. An acceptable coverage is determined by the needs of the application. So, three additional stopping criteria available require the relative change in the mean, standard deviation,

Figure 7. Percent change in confidence interval with calculated bounds for 99 and 95% confidence levels.

and/or the confidence interval for the makespan to be within acceptable tolerances. Once the algorithm stops, a due date still needs to be determined which meets the reliability criteria within the confidence required. This is easily done by solving eq 18 iteratively for f with N, Pu, and R fixed. Once f is found, we determine the due date which leads to f due date violations based on the data which have been collected. 3.4. Parallelizability. We utilize a master slave paradigm to achieve parallelization of the Monte Carlo framework. The master spawns the slaves, collects


results (calculated makespans), and checks stopping criteria. The slaves randomly select processing times and determine equipment availability patterns, solve the planning/scheduling problem thus specified, and then return the results to the master process. Due to the independence of each iteration, it is easy to parallelize the framework to achieve almost linear speedup. The only bottleneck to perfectly linear speedup is the time needed to start up the slaves and to pass results from the slaves to the master. For this work a messagepassing library called parallel virtual machine (PVM) was used (Geist et al., 1994). This package allows the user to have processes running on different architecture machines, thereby utilizing all of the computational resources available within a heterogeneous environment. The modification of the code to utilize PVM is very straightforward and took less than one person-day once the serial code was working correctly. 4. Implications of Framework in Practice and Theory There are two important issues associated with the framework just presented which need to be clarified. The first matter concerns discretization width and its implication on uncertain parameters. As alluded to in section 3.1.1, in this framework only fluctuations in processing times which are on the order of the discretization width are considered. For most large problems, this means that only significant upsets to the plant not minor variations are accommodated. In theory, if minor variations are important, this framework can deal with them but at the expense of longer solution times. In practice, studies by Honkomp (1995) showed that minor variations in processing times had minimal effects on the overall solution of a scheduling problem. The second issue involves the use of a priori knowledge concerning equipment availability. If the problem were solved only once, the use of this perfect knowledge would be unacceptable, but not if the problem is solved many times, with each instance a randomly generated set of perfect information. Taking this to the extreme, in theory, a solution (schedule) is obtained for the set of all possible outcomes. Based on the information from these schedules, operating policies (e.g., due dates) are determined and operation of the facility is initialized. When an upset occurs, the facility is rescheduled. Since all the possible outcomes for the problem are known, this new schedule must be contained in the set of known outcomes. With each new upset, a new schedule is obtained which is always contained in the set of possible outcomes. So in the end, a schedule which was obtained without perfect knowledge is contained in the set of outcomes based on perfect knowledge. The framework presented generates this set of known outcomes. Of course, in practice, we do not determine the set of all known outcomes but collect data for a random subset and use statistics to bridge the gap. In addition, rescheduling after every upset is impractical, so use of this framework in practice will give results which differ slightly from reality. With an understanding of the implications presented, the framework presented can be a strong tool for determining a number of operating policies. 5. Examples All the runs were performed on a parallel virtual machine made up of 21 HP workstations having a total

Figure 8. Recipe for industrial case study. Table 1. Task-Unit Assignments and Processing Times for Industrial Case Study unit task U1, U2, U3 U4, U5, U6 U7, U8, U9 T1 T2 T3 T4 T5 T6 T7 T8 T9 T10 T11 T12 T13 T14 T15

U10

U11, U12

U(7,17) U(1,5) U(2,12) U(1,11) U(7,17) U(1,5) U(1,11) U(6,16) U(2,12) U(1,5) U(7,17) U(1,3) U(3,13) U(1,7) U(1,5)

of 22 processors. The processors used were (with approximate relative speeds) two 9000/770’s (1.00), one 9000/755 (0.66), one 9000/735 (0.80), and 18 9000/715’s (0.33). This gives a cumulative computing power of ∼9.5 9000/770’s, so results were obtained 9.5 times faster by running in parallel than running on a single 9000/770 processor. 5.1. Example 1. This example is based on an industrial case study. The recipe, possible task-unit assignments, and processing times are given in Figure 8 and Table 1. The final product is generated by a three-stage process consisting of 15 tasks on 12 units utilizing 22 resources. Each stage produces a storeable intermediate: I1, I2, and P1. In addition, many of the raw materials can be recovered, stored, and recycled to reduce material costs. As can be seen in Table 1 there is a high degree of sharing involved in this process such that many tasks can be performed on the same unit as well as many units being able to perform the same task. The first part presents the case when the due date is given, while the second looks at determining the due date. 5.1.1. Single Product Due Date Reliability. For this example we will determine the reliability of meeting a demand of 175 000 units of P1 within 15 weeks (2520 h) of starting production with a 99% confidence. In addition to the processing time uncertainties, the equipment is assumed to be down approximately 10% of the time. Since no specific information concerning the typical failure and repair rates for the equipment are given, two Monte Carlo runs are performed. The first run assumes short, frequent breakdowns for all equipment (mttf ) 72 h, mttr ) 8 h), while the second run assumes long, infrequent breakdowns (mttf ) 216 h, mttr ) 24 h). For each randomly generated instance, the problem was solved using both the APP model and the reverse


Figure 9. Progression of Pu.

RWA. The UDM formulation by itself was not used due to the high dimensionality of the problem being solved. By using the other two formulations, we can directly compare the results obtained without bias. For the APP model aggregate intervals of 24 h were used with equipment downtime aggregated for each specific interval, not averaged over the horizon. The reverse RWA utilized a window with a MILP region of 24 h and no LP region. Each subproblem within the reverse RWA was solved to within at least 1% of optimality. The stopping criteria used was that the slope of the failure probability for both solution methods had to be within (0.25%. For this example, a total of 442 iterations were needed for the first run and 441 iterations for the second, with each run taking ∼7 h to complete. The progression of the calculated failure probability is shown in Figure 9. The most obvious point of these plots is that the APP model vastly underestimates the likelihood of missing the due date (only ∼10% of the time), while the reverse RWA, by determining a realistic production profile, shows that the due date is very likely to be missed much of the time (between ∼40% and ∼60%). Another point to notice from the plots is that the downtime profile had very little effect on the failure probability for the APP model, while it had a significant effect for the reverse RWA (∼20% difference between profiles). This can be explained by the discrete nature of batches within the UDM formulation of the reverse RWA. Within the UDM formulation, when a unit fails, any task running on the unit is lost and must be run again once the unit is repaired. The more often a unit goes down, the more often batches have to be restarted. This is reflected in a higher failure probability for the case with frequent, short breakdowns. Since the APP model does not treat batches as discrete entities, only the capacity of the batch is reduced when a unit fails, leading to possibly large overestimations in the production capacity. As long as the total downtime is the same, the failure probability is unaffected by the difference in profiles. From this example we were able to see that, for this facility, it would be unwise to determine operating policies based on the APP model. 5.1.2. Single Product Promise Date Determination. Since we have shown that the original due date in the previous example is highly unreliable, the obvious next step is to determine a promise date which is

Figure 10. Histograms for makespan for the reverse RWA.

Figure 11. Progression of mean and standard deviation for the reverse RWA.

reliable. In this case, we would like a promise date which can be met 95% of the time with 99% confidence. The same two runs were performed as in the previous example, but this time only the reverse RWA was used since the APP model is unsatisfactory for this problem. The stopping criteria used was when the maximum changes in both the mean and standard deviation of the makespan were no more than (0.25%. For this stopping criteria, a total of 1105 iterations was needed for both runs, taking ∼18 h of computing time each. As shown in Figure 10, the makespan appears to be normally distributed, so using this assumption within the stopping criteria is valid. From Figure 11 it is obvious that both the mean and standard deviation had settled down well before 1105 iterations. So a less stringent stopping criteria would have allowed the algorithm to terminate sooner without loss of detail, thus reducing the overall solution time. Also, both the mean and standard deviation are larger for the case where there are short, frequent breakdowns. The reason for this follows from the discrete nature of batches discussed in the previous example: more frequent breakdowns lead to more interrupted batches which lead to longer makespans. Given N ) 1105, R ) 0.99, and Pu ) ∼5%, eq 18 can be solved iteratively, giving a value of f ) 39. From the data collected, this corresponds to a promise date of 3334 h for the case of short, frequent breakdowns and a promise date of 3021 h for long, infrequent breakdowns.


Figure 12. Recipe information for example 2.

Figure 13. Transportation information for example 2. Table 2. Demand Distribution for a Multisite Problem warehouse product

W1

W2

P1 P2 P3

120 000 50 000 20 000

60 000 40 000 70 000

The differences in due dates and reliabilities obtained for different equipment downtime profiles emphasize the importance of having good information about equipment breakdowns. It also helps to underscore the benefits of preventive maintenance. The longer between breakdowns, the greater the reliability and, in turn, the shorter the necessary makespan. So an operating policy may be that when a unit goes down, it may be better to leave it down longer to fix more than just what caused the failure. This extended downtime may help to reduce the frequency of future downtimes. 5.2. Example 2. The previous example dealt with only a single demand. Multiple demands, while possibly taking more time, do not significantly increase the complexity of the problem. One method is to prioritize the demands and solve the problem for each demand sequentially. As each demand is added, the makespan is determined and further production is prevented during the time each unit has already been utilized. Another method is to solve the entire problem all at once and “integrate” over the production profile(s) to determine makespans necessary to meet each demand. The latter method is the approach utilized for this example. This example involves a multisite problem based on industrial data first presented in Bassett et al. (1996a). The problem consists of two plants producing three products which are shipped to two warehouses (see Figures 12 and 13). Each product takes roughly five tasks to be produced, with two products being produced in one plant and all three in the other. For this example, since all the demands given in Table 2 have to be met within 1440 h with equal priority, missing any demand is unacceptable. 5.2.1. Multiple-Product Due Date Reliability. First, we would like to determine the probability of meeting all the demands, as well as only five of the six

demands, within 2 mo (1440 h) with a 99% confidence. As in the previous example, the equipment is assumed to be down ∼10% of the time with a mttf of 216 h and mttr of 24 h. Only the reverse RWA was used for this problem with an MILP region of 27 h and no LP region. Each subproblem was solved to within 1% of optimality. The stopping criterion used was that the slope of the failure probability for both cases had to be within (0.25%. After 261 iterations and ∼2.5 h, the algorithm stopped with a failure probability of ∼75% for meeting all of the demands and ∼74% for meeting five of the six demands. The reason these probabilities are so close is that since there are two plants shipping to both warehouses, it is not uncommon to have multiple demands satisfied at approximately the same time. Therefore if one demand is missed, probably two are missed. In other words, in order to meet five of six demands, all six demands are likely to be met. 5.2.2. Multiple-Product Promise Date Determination. Finally, we would like to determine a promise date for this problem at which time all of the demands can be met 95% of the time with 99% confidence. The stopping criterion used was when the maximum change in both the mean and standard deviation of the makespan was no more than (0.5%. For this stopping criterion, a total of 569 iterations was needed, requiring ∼5 h of computing time. Given N ) 569, R ) 0.99, and Pu ) ∼5%, eq 18 can be solved iteratively, giving a value of f ) 16. From the data collected, this corresponds to a promise date of 2075 h, which is ∼44% later than the original due date of 1440 h. 6. Conclusions In this paper we have shown that for highly shared batch facilities it may be necessary to use detailed scheduling to obtain realistic production plans rather than an aggregate production planning model. We presented a Monte Carlo framework to determine the reliability of given due dates or determine due dates which met a certain reliability when processing times and equipment reliability were uncertain. Using this framework, we also showed the importance of having good information concerning equipment reliability (section 5.1). Additionally, we derived a number of stopping criteria which may be used within the Monte Carlo framework. If information about intermediate storage is needed, these data are easily obtained by this framework as well. From each iteration we obtain a detailed schedule; therefore, we know storage levels for all resources at all times. Using all the iterations, we can obtain distributions of maximum and/or minimum storage levels for every time period and/or over the whole horizon. Distributions can also be obtained for other quantities such as equipment utilization or manpower. Nomenclature Sets I ) set of all tasks Ij ) set of all tasks which can run on unit j J ) set of all units Ji ) set of all units which can run task i Js ) set of all units which can store resource s S ) set of all resources T ) set of all aggregate time periods

1726 Ind. Eng. Chem. Res., Vol. 36, No. 5, 1997 T′ ) set of all discrete time periods Parameters Cijt ) production penalty to run task i on unit j at time t Dst ) demand for resource s during aggregate interval t dst′ ) demand for resource s at discrete time t′ consume ) fraction of task i on unit j which consumes fijs resource s create ) fraction of task i on unit j which creates resource fijs s Ht ) length of aggregate interval t pijt′ ) processing time for task i on unit j at discrete time t′ Vijmax ) maximum production quantity for task i on unit j Vijmin ) minimum production quantity for task i on unit j Variables Asjt ) quantity of resource s stored in unit j at aggregate time t asjt′ ) quantity of resource s stored in unit j at discrete time t′ Bijt ) production quantity for task i on unit j at aggregate time t bijt′ ) production quantity for task i on unit j at discrete time t′ QBst ) quantity of resource s purchased at aggregate time t qbst′ ) quantity of resource s purchased at discrete time t′ Wijt ) assignment variable for task i on unit j at aggregate time t wijt′ ) assignment variable for task i on unit j at discrete time t′

Literature Cited Bassett, M. H.; Doyle, F. J., III; Kudva, G. K.; Pekny, J. F.; Reklaitis, G. V.; Subrahmanyam, S.; Zentner, M. G.; Miller, D. L. Perspectives on model-based integration of process operations. Comput. Chem. Eng. 1996a, 20, 821. Bassett, M. H.; Pekny, J. F.; Reklaitis, G. V. Decomposition techniques for the solution of large scale scheduling problems associated with batch chemical processes. AIChE J. 1996b, 42, 3373. Bowman, E. H. The scheduling sequencing problem. Oper. Res. 1959, 7, 621. Geist, A.; Beguelian, A.; Dongarra, J.; Jiang, W.; Manchek, R.; Sunderam, V. PVM3 User’s Guide and Reference Manual; Oak Ridge National Laboratory: Oak Ridge, TN, 1994. Honkomp, S. Development of a model based testing framework for reactive scheduling under uncertainty. M.S. Thesis. School of Chemical Engineering, Purdue University, West Lafayette, IN, 1995.

Ierapetritou, M. G.; Pistikopoulos, E. N. Batch plant design and operation under uncertainty. Ind. Eng. Chem. Res. 1996, 35, 772. Kondili, E.; Pantelides, C. C.; Sargent, R. W. H. A geneal algorithm for scheduling batch operations. Third International Symposium on Process Systems Engineering, Sydney, Australia, 1988; p 62. Kondili, E.; Pantelides, C. C.; Sargent, R. W. H. A general algorithm for short-term scheduling of batch operations. Part IsMathematical Formulation. Comput. Chem. Eng. 1993, 17, 211. Lee, Y.-G.; Malone, M. F. A general treatment of uncertainties in batch process scheduling. AIChE 1995 Annual Meeting, Miami Beach, FL, Nov 1995; Paper 174b. Lipson, C.; Sheth, N. J. Statistical design and analysis of engineering experiments; McGraw-Hill Publishing Co.: New York, 1973. Manne, A. S. On the job-shop scheduling problem. Oper. Res. 1960, 8, 219. McKay, K. N.; Safayeni, F. R.; Buzacott, J. A. A review of hierarchical production planning and its applicability for modern manufacturing. Prod. planning control 1995, 6, 384. Pantelides, C. C. Unified framework for optimal process planning and scheduling. Proceedings of the Second Conference of the Foundations of Computer Aided Operations, Crested Butte, CO, 1993; p 253. Pekny, J. F.; Zentner, M. G. Learning to solve process scheduling problems: the role of rigorous knowledge acquisition frameworks. Proceedings of the Second Conference of the Foundations of Computer Aided Operations, Crested Butte, CO, 1993; p 275. Pieruschka, E. Principles of reliability; Prentice Hall: Englewood Cliff, NJ, 1963. Sahinidis, N. V.; Grossmann, I. Reformulation of multiperiod MILP models for planning and scheduling of chemical processes. Comput. Chem. Eng. 1991, 15, 255. Shah, N.; Pantelides, C. C.; Sargent, R. W. H. A general algorithm for short-term scheduling of batch operations. Part IIs computational issues. Comput. Chem. Eng. 1993, 17, 229. Thompson, S. D.; Wantanabe, D. T.; Davis, W. J. A comparative study of aggregate production planning strategies under conditions of uncertainty and cyclic product demands. Int. J. Prod. Res. 1993, 31, 1957. Zentner, M. G.; Pekny, J. F.; Reklaitis, G. V.; Gupta, J. N. D. Practical considerations in using model-based optimization for the scheduling and planning of batch/semicontinuous processes. J. Process Control 1994, 4, 259.

Received for review July 31, 1996 Revised manuscript received February 18, 1997 Accepted February 19, 1997X IE960470V

X Abstract published in Advance ACS Abstracts, April 1, 1997.

Using Detailed Scheduling To Obtain Realistic Operating Policies for a

Recommend Documents