Order release planning with clearing functions: A queueing-theoretical analysis of the clearing function concept

Int. J. Production Economics 131 (2011) 399–406

Contents lists available at ScienceDirect

Int. J. Production Economics journal homepage: www.elsevier.com/locate/ijpe

Order release planning with clearing functions: A queueing-theoretical analysis of the clearing function concept Hubert Missbauer  Department of Information Systems, Production and Logistics Management, University of Innsbruck, Universit¨ atsstrasse 15, A-6020 Innsbruck, Austria

a r t i c l e in fo

abstract

Article history: Received 1 June 2008 Accepted 7 September 2009 Available online 21 October 2009

We consider the problem of planning future order releases in hierarchical production planning and control systems. An established research direction is the clearing function concept: the planned material flow through a production unit is modelled by inventory balance equations for WIP and final products, and the consequences of the stochastic properties of the material flow are modelled by clearing functions, which is the functional relationship between the level of WIP and the maximum output of a work centre in a period. Using a transient M/M/1 model, our paper shows that the usual definition of a nonlinear clearing function suffers from substantial shortcomings concerning both the definition of the function and empirical estimation of its parameters. We propose an alternative transient clearing function and derive a procedure for its parameterization. & 2010 Published by Elsevier B.V.

Keywords: Production planning Order release Workload control Queueing theory

1. Problem description We consider the problem of determining future order releases to a production unit where detailed scheduling decisions within the production unit are performed locally at the shop floor level. This is an important decision problem in production planning and control, because order releases determine the level of work-inprocess (WIP) in the production unit and thus the average flow times and the possible output. These flow times, in turn, must be known in order to release the work at the appropriate times and to maintain high due-date performance. In the following we assume that this load-dependence of the flow times should not be considered by setting fixed lead times and determining the order releases based on these target lead times (as is usually done in traditional workload control), but by determining order releases and (variable) lead times simultaneously. That is, lead times should be adjusted dynamically according to the required output and the load situation in the shop. Models have been developed for this planning task that simultaneously determine future order releases, lead times and WIP levels. This paper aims at contributing to improvements of these models. The paper is organised as follows. Section 2 reviews relevant literature. Section 3 describes the basic structure of the aggregate order release planning models under consideration. Section 4 describes the clearing function as a model of the stochasticity of the material flow at the work centres and discusses its  Tel.: + 43 512 507 7503; fax: +43 512 507 2847.

E-mail address: [email protected] 0925-5273/$ - see front matter & 2010 Published by Elsevier B.V. doi:10.1016/j.ijpe.2009.09.003

shortcomings. Section 5 presents a queueing-theoretical analysis of the clearing function concept and derives a new transient clearing function together with rules for its parameterization based on empirical data. Section 6 provides some conclusions and research directions.

2. Related literature The roles of lead time management and of the order release function in hierarchical production planning and control systems where detailed scheduling is performed locally within the production units has been analyzed in several publications, mainly in the 1980s (Tatsiopoulos and Kingsman, 1983; Bertrand and Wortmann, 1981; Bertrand et al., 1990). The workload control concept that aims at controlling lead times by controlling order releases and thus the level of WIP in the production unit emerges ¨ from this idea (description in Bertrand et al., 1990; Zapfel and Missbauer, 1993). Workload control can be performed by order release mechanisms that aim at balancing work input and capacity for a short planning horizon and (mostly) given capacities. Reviews of this technique can be found in Bergamaschi et al. (1997), Land (2004), Stevenson and Hendry (2006). Due to their short time horizon, the planning capabilities of these order release mechanisms are very limited. An alternative approach (which might be complementary for performing multi-period planning) which is more related to linear programming models for production planning (Johnson and Montgomery, 1974, p. 187ff.) and Input/Output Control

400

H. Missbauer / Int. J. Production Economics 131 (2011) 399–406

(Wight, 1974) is the following: the planned order release volumes over time (at the level of work orders, products or aggregate products) are evaluated using a multi-period model of the material flow through the network of work centres that constitute the production unit. Appropriate mathematical programming techniques are used to optimise the order release volumes for the periods of the specified planning horizon. Methods for flow time prediction and control must be integrated into this model. If fixed lead times are assumed for all work that is released, the resulting models are largely state of the art (Hackman and Leachman, 1989; Voss and Woodruff, 2006). However, queueing theory and practical experience indicate that flow times strongly depend on the workload or capacity utilization of the facilities, which is determined by the release decisions. Hence the planned lead times should be adjusted accordingly. The development of order release planning models with load-dependent lead times is still an active area of research today. Two main approaches can be distinguished: In flow–timeoriented models, estimated flow times for the released orders are generated that depend on the order release plan, and hence on the capacity load. These estimated flow times can be derived from operational characteristics (e.g., relationship between capacity load and flow time) that represent long-term averages (e.g., Voss and Woodruff, 2006, p. 164ff.), or from a simulation or queueing model of the production unit, which requires iterative adjustment of order release planning and flow time prediction (Hung and ˜ o, 2003). In contrast, WIP-oriented models Leachman, 1996; Rian represent the flow of WIP through the work centres by inventory balance equations. Flow times are not represented explicitly, although they can be inferred from the WIP level at a certain time and the output in the future periods. For deterministic, continuous material flow it is possible to formulate and solve models of this type in continuous time (fluid relaxation; see Bertsimas and Sethuraman, 2002), but for order release, discrete-time models are by far more common (literature given below). In the following we pursue this latter direction and concentrate on the decisive research topic, namely the integration of the stochastic properties of the material flow into WIP-oriented order release models in discrete time.

3. Models for aggregate order release planning We consider a production unit consisting of M work centres (m= 1, y, M) that produces J products or product groups (j= 1, y, J), which may represent final products or components. Demand forecasts Djt for the products or product groups over the planning horizon t =1, y, T are available. These forecasts can be the customer demand (if the j denote final products) or the demand for components that results from the MRP system. The order releases Rjt are to be determined. As mentioned above, there are different approaches to perform order release planning in this context (for an overview, see Missbauer and Uzsoy, 2008). In the following we pursue an important approach that is at least close to practical applicability. The flow of work through the network of work centres that forms the production unit is modelled by means of inventory balance equations, complemented by a clearing function for each work centre that models the effect of the stochasticity. The basic formulation of the model for aggregate products with similar routings within the product groups is as follows: Variables (all measured in hours of work) Wjmt

work of product group j waiting at work centre m at the end of period t.

Ijt Rjt Xjmt Ajmt

finished goods inventory of product group j at the end of period t. released work of product group j in period t. output of product group j from work centre m in period t. input of product group j to work centre m in period t.

Parameters Djt p~ jim

zjimt

demand for product group j in period t (measured in hours of work). average work arriving at work centre m when one unit of product group j is finished at work centre i ( = transition probabilities pjim times the ratio of the operation times at work centres i and m). proportion of the output of product group j from work centre i to work centre m that arrives at m t periods after completion at i.

Generic model structure: Wjmt ¼ Wj;m;t1 þ Ajmt Xjmt

Ajmt ¼

M X 1 X

for all j; m; t

Xj;i;tt p~ jim zjimt þ

i¼1t¼0

1 X

ð1Þ

Rj;tt p~ j0m zj0mt

for all j; m; t

t¼0

ð2Þ This is the definition of the WIP at the work centres Wjmt. The sources of the input Ajmt are the output of the other work centres and the release of orders. The work input can be delayed due to transportation, etc., which is represented by the t in the subscripts. i=0 denotes order release. Ijt ¼ Ij;t1 þ

M X 1 X

Xj;m;tt p~ jm0 zjm0t Djt

for all j; t

ð3Þ

m¼1t¼0

This is the definition of the final product inventory. Work centre index 0 in p and z denotes completion. If the work is not aggregated to product groups and separate products or work orders are considered, the model structure is basically the same. Explicit representation of setups and lot sizing represents a major extension to the model which is not treated here. The objective function is an appropriate cost function, consisting of the holding cost for WIP and finished goods inventory and possibly other cost components. Note that this model is deterministic; it is analogous to fluid approximation in queueing theory (Kleinrock, 1976; Chen and Mandelbaum, 1994). Hence the model assumes that work that is scheduled to be available at a work centre i in period t can be processed (cleared) in period t up to the maximum capacity. This assumption usually overestimates the output due to the variability of the material flow. Clearing functions have been proposed to capture the effects of this stochastic behaviour; these are described and criticised in the next section.

4. Clearing functions in order release models In order to incorporate the stochastic effects into model (1)–(3), we must estimate the fraction of the planned available work at work centre i in period t that is processed (cleared) in period t. A clearing function can be defined as the expected or maximum output of a work centre i in period t as a function of some measure of WIP (e.g., average WIP or planned available work) in period t and of the capacity Cit of the work centre. In this

H. Missbauer / Int. J. Production Economics 131 (2011) 399–406

401

Initial WIP W0 = Ls (0)

Poisson input, E [Ak] = λ τ t Time of planning t = -(k-1) τ

t=0

k-th period No. 1

t=τ

Fig. 1. Integration of the M/M/1 model into the order release planning model.

paper we use the amount of available work as the WIP measure. This is termed the load of work centre i in period t (denoted by Lit) and is defined as the WIP at the beginning of period t plus the planned input in period t (Ait): Lit ¼ Wi;t1 þ Ait

ð4Þ

Hence the clearing function for work centre i is a functional relationship of the form Xit ¼ fi ðWi;t1 þ Ait ; Cit Þ

ð5Þ

Starting with Karmarkar (1989), concave, saturating clearing functions have been widely used in the literature (Srinivasan et al., 1988; Missbauer, 1998, 2002; Asmundsson et al., 2002, 2006, 2009) and can be regarded as a standard modelling technique. However, clearing functions of this type assume that all factors that influence the expected output in a period for a given expected load (e.g., the ratio between Wi,t  1 and Ait) are either kept constant or are determined by the expected load. Thus a clearing function can be derived, e.g., for a steady-state period of a single-stage queueing system or for the kth period of a singlestage queueing system starting at the origin. If, as is common, the clearing function is derived from simulated data, it reflects the system states that occurred in the simulation. Hence the use of this type of clearing functions in order release planning models constitutes a logical contradiction. Order release planning models determine the time-dependent input to the production units and hence the extent to which the (e.g., steady-state) assumptions on which the clearing function is based can hold. Combining this with a conventional clearing function is problematic. There is strong evidence that this in fact leads to shortcomings in the performance of the models (Missbauer, 2009). This is less important for long planning periods, because in this case it is more likely that the system reaches the state (e.g., steady-state) the clearing function is based on. A look at the historical development of the clearing function idea is helpful here: Clearing functions were introduced by Graves (1986). He assumed a linear (proportional) clearing function that is consistent with constant flow times (Graves, 1986, p. 524f.) and thus models production control practices that stabilize the flow times, e.g., variation of the processing speed, short-term capacity expansion, changes in worker assignment (see Graves, 1986, p. 524). Note that there is no specific reference to stochastic arrival and departure processes at the work centres; a proportional clearing function only states that in each planning period a constant proportion of the available work is processed. In contrast, the concave, saturating shape of the clearing function introduced by Karmarkar (1989) explicitly addresses the congestion behaviour of capacity-constrained production resources (Karmarkar, 1989, p. 109) suggested by queueing models. The reason for congestion when utilizationo1 is the stochastic nature of arrival and departure process. This is a major shift – from modelling production control practices to modelling the behaviour of stochastic systems – whose implications are not yet fully understood. Surely the long-term relationship between average

WIP and average output that results from the stochastic properties of the manufacturing system is well-known in the literature and is an important basis for production control, especially for order release mechanisms based on the workload control concept. However, optimization models for order release planning require models of the material flow over relatively short periods. The resulting difficulties in interpreting and handling saturating clearing functions intended to model the stochastic behaviour of production resources lead to substantial shortcomings in the research on this topic. We now present a new approach to overcome these shortcomings.

5. Definition and parameter setting of a transient clearing function 5.1. Queueing-theoretical analysis of the clearing function in a transient state We assume a single-server M/M/1 queueing system with service rate m = 1 and an initial WIP (orders in the queue and at the server) of Wk  1 at the beginning of period k of length t time units. An input of E[Ak]= l  t (with l denoting the arrival rate in period k) is expected. The expected output E[Xk] is to be calculated. Wk  1, Ak and Xk are the decision variables of the order release model. In general the actual WIP, input and output represented by these variables are random variables at the time of planning. For exposition we express WIP and input as number of customers, while the output Xk is measured in units of time (e.g., standard hours). The model is depicted in Fig. 1. We now derive the expression for the expected output E[Xk] and examine whether it can be modelled to an acceptable degree of accuracy as a function of a single variable, namely the available work (load) in period k. We define: prn(t)

conditional probability of having n customers in the system at time t given r customers in the system at time 0.

An analytical expression for prn(t) as a function of the service rate m and the arrival rate l for the M/M/1 model can be derived (see Cohen 1969, p.82ff. and p. 178; the expression is not given here). The quantity most interesting to us is pr0(t)=P{Ls(t)=0| Ls(0) = r} the probability of idleness at time t given r customers in the system at time 0. Ls(t) denotes the number in system at time t.1 In the following we let k-1= 0 without loss of generality; that is, period k is scaled as period 1 such that t= 0 at the start of the 1 In this paper the index for the periods (discrete time) is denoted as a subscript; continuous time is denoted in parenthesis. Since in this section we consider the period with number 1 (even if it is not the first period), the start of the period is at time t= 0, which means that Ls0 =Ls(0). In the formulas that specifically consider the data of period 1, we use the symbol Ls0.

402

H. Missbauer / Int. J. Production Economics 131 (2011) 399–406

6

Expected output E [X1]

5

4

3

2 Ls(0)=0 Ls(0)=1 Ls(0)=2 Ls(0) =3 Ls(0)=4 Ideal curve

1

0 0

1

2

3

4

5

6

7

8

9

10

11

Expected load E [L1] Fig. 2. Clearing functions for different deterministic initial WIP.

period. Then the expected output E[X1] for a given deterministic initial WIP Ls(0) becomes: Z t E½X1  ¼ t pLsð0Þ;0 ðtÞdt ð6Þ t¼0

In general Ls(0) is a random variable. Assuming no systematic bias, its planned value (given by the value of the corresponding decision variable in the optimization model) can be interpreted as the expected value of the distribution, denoted by E[Ls(0)]. We define pn(t) as the probability of n customers in the system at time t. Then we have: p0 ðtÞ ¼

1 X

pr ð0Þ:pr0 ðtÞ

0ot rt

ð7Þ

r¼0

The expected output, measured in units of time, is Z t E½X1  ¼ t p0 ðtÞdt

ð8Þ

t¼0

By substituting (7) for p0(t), this yields: Z t 1 X E½X1  ¼ t pr ð0Þ:pr0 ðtÞdt

ð9Þ

t¼0r¼0

The probability distribution of the number in system at the end of period k= 1 (time t) can be obtained analogously as: pn ðtÞ ¼

1 X

pr ð0Þ:prn ðtÞ

ð10Þ

r¼0

Thus the pn(0) (note that t = 0 is just our time scale and not the beginning of the planning horizon) can be calculated recursively from the pr(  t) analogous to (10); the recursion can continue to the beginning of the planning horizon. However, this is not an applicable procedure in the context of order release planning models. Therefore we explore which independent variables determine the expected output in a period k to a sufficient degree of accuracy. First we examine whether the one-dimensional clearing function (5) can be used or whether the composition of the expected load (the relative proportion of E[Ls0] and E[A1]) has a significant impact on the expected output, which would mean that (at least) a two-dimensional clearing of the form E½X1  ¼ f2 ðE½Ls0 ; E½A1 Þ

one-dimensional clearing function (5) does not provide the required precision. Fig. 2 depicts the expected output E[X1] (9) as a function of the expected load (the total available work) of the period, defined as E½L1  ¼ E½Ls0  þ l:t

ð12Þ

for the length of the period t = 5 and the deterministic number in system at the start of the period Ls0 = Ls(0) varying from 0 to 4. These are the clearing functions for deterministic initial state Ls0 and different proportions of initial WIP and work input. l is the control variable, and the output values for identical values of the expected load on the x-axis correspond to different l values for different values of the initial WIP (see (12)).2 The clearing functions differ substantially, showing that for a given expected load the expected output is substantially higher as the initial number in system Ls0 increases (and the expected input decreases accordingly). In this case (deterministic initial state) the variance of the load, expressed in orders, depends on the ratio of initial WIP to expected input which might partly explain the difference. Hence a two-dimensional clearing function (11) with initial WIP and expected input as two independent variables in the spirit of Andersson et al. (1981) can be expected to improve the accuracy of the clearing function estimate compared to the usual one-dimensional clearing function. If the initial WIP is measured in time units (e.g., hours of work), the clearing function can be derived analogously (see Missbauer, 2007). In this case the dependence of the expected output on the initial WIP is stronger than in Fig. 2 because if W0 ZC1 the output X1 =C1 (ideal curve in Fig. 2). Next we examine the relationship between E[X1] and the distribution pr(0)8r given by (9) for a given expected value P E½Ls0  ¼ r r:pr ð0Þ. This analysis should clarify whether a two-dimensional clearing function (11) is sufficient for a given input process (which determines the variability of the input A1) or whether it is necessary to include measures of the variability of Ls0 (e.g., the standard deviation) as independent variable(s). Note that the standard deviation of Ls0 is zero for the first period and presumably will increase in the later periods.

ð11Þ

must be applied. Since the prn(t) are nearly intractable analytically even for the M/M/1 model, only numerical analyses are possible. Therefore we use a numerical example to show that the

2 Due to the time-consuming computations, the numerical precision of the calculation is limited, which leads to the slight oscillations in the curves of Figs. 2 and 3.

H. Missbauer / Int. J. Production Economics 131 (2011) 399–406

Since complete independence of the expected output according to (9) from the standard deviation of Ls0 cannot be expected, it is important to specify the possible range of this standard deviation or the family of distributions pn(0). The following extreme or idealized cases are relevant:

 For the first period the number in system at the beginning is deterministic.

 If the system is in steady-state, the number in system follows the distribution: pn ð0Þ ¼ ð1rÞ:rn



for n ¼ 0; 1; . . .

ð13Þ

with r = l/m. If the input process is a Poisson process with time-dependent input rate (M(t)/M/1-system) and E[Ls(0)] is a fixed value, this can result from a variety of different patterns of l(t) up to t = 0. If we assume a step function as

lðtÞ ¼ 0 for t o tstart ¼ lstart for tstart rt r 0

ð14Þ

then the limit for  tstart-0 is a pulse input at t =0, and the number in system follows a Poisson distribution with pn ð0Þ ¼

ln n!

el

for n ¼ 0; 1; . . .

ð15Þ

Fig. 3 depicts the expected output E[X1] (9) as a function of the expected load (the total available work) of the period (12) for the length of the period t = 5, expected number in system at time zero E[Ls0] =4 and pn(0) according to the three distributions of Ls0 mentioned above: deterministic, Poisson and steady-state. Again l is the control variable that determines the expected load (12) and the expected output (9). Although the expected initial WIP and the expected work input (and hence the expected load) for given l are the same for all three situations (from (12)), the clearing functions differ substantially, which can only be due to the Ls0 distribution. The expected output for a given expected load decreases as the variability of Ls0 increasespfrom coefficient of variation CV = 0 ffiffiffi (deterministic) to CV ¼ 1= l ¼ 0:5 (Poisson distribution with mean l = 4) to CV =1.12 (steady-state distribution of the M/M/1 system with mean r/(1  r)= 4). The output from the Poisson distribution (which is the consequence of a pulse input) and the

403

output from the steady-state distribution differ by 6–22% which is clearly a significant impact. Thus the impact of the variability of the number in system at the beginning of the period and of the process leading to this variability should not be ignored. This variability is a consequence of the history of the process (work input and output before the beginning of the period under consideration), and it seems reasonable to consider it the essential variable to represent the impact of the history on the output of the period. From this analysis we can conclude the following: for a given stationary input process during a planning period, the expected output in this period can be modelled as a function of three separate variables:

 the expected initial WIP level,  the expected input during the period,  the probability distribution of the WIP level at the beginning of the period. It seems reasonable to express the distribution of the initial WIP level as the mean and the standard deviation (or the coefficient of variation) of this distribution. This yields a three-dimensional clearing function with the expected initial WIP, the expected input and an appropriate measure of the variability of the initial WIP as independent variables of the regression. It remains to be analyzed whether including higher moments of the initial WIP distribution is worthwhile, but it seems difficult to manage a clearing function with four or more independent variables. The inclusion of the variability of the initial WIP in a future period, estimated at the beginning of the planning horizon, as independent variable in a regression leads to a conceptual problem: It is not observable in empirical data. In order to analyze this problem we must review some theoretical issues of clearing functions and their empirical estimation.

5.2. Estimation of clearing functions from empirical data Clearing functions usually are estimated from empirical or simulated data as follows. For each period of the time horizon under consideration the WIP measure (average WIP level; load as defined above) and the total output in the period are recorded. Fig. 4 shows the data for a printing machine of a manufacturer of optical storage media.

6

Expected output E [X1]

5 4 3 2 1 0 0

2

4

6 Expected load E [L1]

8

Steady-state initial WIP distribution

Deterministic initial WIP

Poisson initial WIP distribution

Ideal curve

10

Fig. 3. Clearing functions for period 1 for different distributions of the initial WIP.

12

404

H. Missbauer / Int. J. Production Economics 131 (2011) 399–406

We assume that the initial WIP Wk  1 and the total input during period k, denoted Ak, are expressed in number of orders and are random variables following the joint distribution

Actual output Ideal shape

Output

fWk1 ;Ak ðw; yÞ probability of initial WIP Wk  1 =w and input Ak =y. The marginal distributions are fAk ðyÞ ¼

1 X

fWk1 ;Ak ðw; yÞ

ð17Þ

w¼0

denoting the probability that the input during the period k is y, and fWk1 ðwÞ ¼

fWk1 ;Ak ðw; yÞ

ð18Þ

y¼0

Load (available work) Fig. 4. Empirical data for the clearing function. Each data point is a combination of load (= total available work) and output of one day for a single machine. The ideal shape is the minimum of load and capacity.

giving the probability that the initial WIP at the start of period k is w. The expected output during period k is E½Xk  ¼

In contrast, a clearing function used in an order release model must use the estimates of the initial WIP Wk  1 and work input Ak (made at the time of planning, before uncertainty is realized) as independent variables for the calculation of the expected output (E[Xk]) in future periods. Hence the usual way to obtain a clearing function from empirical data does not yield a correct clearing function for the planning model; it systematically underestimates the stochastic effects by treating a point estimate of a random variable as a deterministic quantity, which leads to a biased result. It is shown below that this bias results in an overestimation of the output that can actually be obtained. We now present a way to compensate for this. The starting point is the following crucial insight: The empirical clearing function (regression over empirical data) is the conditional expectation of the output Xk given the observed (deterministic) values of the initial WIP and of the work input in a period. We assume that the stochastic arrival process at the work centre is determined by the structure of the production unit and its planning and scheduling system. Order release can change only the arrival rate. Based on these assumptions, the correct clearing function for a planning model (with estimated initial WIP and work input as independent variables) can be derived from the empirical regression function (with actual initial WIP and work input as independent variables) as follows. For the following analysis we assume a two-dimensional clearing function: Xk ¼ f ðWk1 ; Ak Þ

1 X

E½Xk jWk1 ¼ w; Ak ¼ y conditional expectation of the output in period k given that the initial WIP Wk  1 =w and the input during the period k Ak = y.

E½Xk jWk1 ¼ w; Ak ¼ y:fWk1 ;Ak ðw; yÞ

ð19Þ

w¼0y¼0

The question is how to determine the probability distribution fWk1 ;Ak ðw; yÞ. The following approximation seems a reasonable starting point:

 Assuming stochastic independence of initial WIP and input





ð16Þ

The empirical data used in the regression (Fig. 4) are the observed (and hence deterministic) values for Xsk, Ws,k  1 and Ask for different periods k and possibly different simulation runs s. The variability that is included in the empirical data (but not observed explicitly) is the distribution of the work input Ak within the period and, if WIP and input are measured in number of orders, the variability of the service times (machine breakdowns are not considered here). Appropriate functional forms for the regression function in the one-dimensional case Xk ¼ f ðWk1 þAk Þ can be derived from queueing theory (see Missbauer, 1998, 2002); for the two-dimensional function (16) there is little guidance available in the current literature. The starting point is the empirical two-dimensional clearing function (16), which yields the expected output for the specified values of Wk  1 and Ak. We define:

1 X 1 X



during the period, the marginal distributions (17) and (18) can be estimated independently. Considering stochastic dependence (which is likely since a workload control approach is applied) is a topic for future research. The probability distribution of the number of arriving orders depends on the arrival process, e.g., a Poisson distribution for Poisson arrival process. The parameter can be determined from the estimated input to the work centre in the optimization model. The probability distribution of the initial WIP represents the history of the process in prior periods. We interpret the variable Wk  1 in the optimization model as the expected initial WIP E[Wk  1]. The second moment of the distribution of the initial WIP cannot be observed and must be calculated from a model. Abate and Whitt describe approximations for the transient behaviour of the first and second moment of the number in system and of the workload in M/M/1 and M/G/1 systems with utilizationo1 (Abate and Whitt, 1987a, 1987b, 1988, 1994). These results, which cannot be described here in detail (see Missbauer, 2007), can be used to obtain appropriate approximations for the variance or coefficient of variation of Wk  1. Working out this procedure in detail is a topic for future research. Determining the distribution of the random variable Wk  1 requires assumptions about the family of distributions, since this information is not available from the first two moments. The properties of analytical WIP distributions for single-stage queueing models can provide some guidance here. Extending the procedure to include an estimation of the third moment, which would specify the family of distributions more precisely, must be considered a substantial complication, especially since less insights into the time-dependent behaviour of these higher moments seem to be available. Thus we do not pursue this direction.

Based on these assumptions, which must be refined and justified in future research, the marginal distributions fWk1 and fAk can be calculated. The assumed stochastic independence leads to fWk1 ;Ak ðw; yÞ ¼ fWk1 ðwÞ:fAk ðyÞ

ð20Þ

H. Missbauer / Int. J. Production Economics 131 (2011) 399–406

The expected output for period k can be calculated by substituting (20) and the empirical clearing function E½Xk jWk 1 ¼ w; Ak ¼ y into (19). 5.3. Numerical example We use the analytical results depicted in Fig. 2 (deterministic initial WIP) as the data for the empirical clearing function. These are the values for E½Xk jWk1 ¼ w; Ak ¼ y. We assume a steadystate distribution of the M/M/1 system for Wk  1: fWk1 ðwÞ ¼ PfWk1 ¼ wg ¼ ð1aÞaw

w ¼ 0; 1; . . .

ð21Þ

with E½Wk1  ¼

a

ð22Þ

1a

as the specified value Wk  1 for the corrected clearing function (i.e., the planned value Wk  1 of the optimization model), which yields:

a¼

E½Wk1  1E½Wk1 

ð23Þ

We assume a Poisson input process, hence: fAk ðyÞ ¼ PfAk ¼ yg ¼

ly y!

el

y ¼ 0; 1; . . .

ð24Þ

with E½Ak  ¼ l ) l ¼ E½Ak 

ð25Þ

The mechanism for the correction is as follows. We want to calculate the expected output E[Xk] as a function of the planned ( =expected) values for Wk  1 and Ak (again assuming that the value of the decision variable Wk  1 is equal to the expected value of the actual WIP at the start of period k). The specified value for E[Wk  1] yields a (23) and hence the distribution fWk1 ðwÞ (21). Likewise, the specified value for E[Ak] yields l and hence the

405

distribution fAk ðyÞ ((25) and (24)). The resulting joint distribution fWk1 ;Ak ðw; yÞ (20) is the basis for the correction according to (19). The numerical example shows the order of magnitude of the correction. We use the same model as in Fig. 2 and the twodimensional clearing function (16) that results from deterministic Wk 1 and Poisson-distributed Ak for m =5 [orders per period]. The resulting values for the expected output E[Xk] are used as empirical clearing function values (note that these are arbitrary data in this context). Fig. 5 depicts the plot of the empirical clearing function. Table 1 shows the relative difference between the empirical and the corrected clearing function. Table 1 shows that the corrected values are lower than the empirical values (which underestimate the stochasticity), suggesting that the application of the empirical regression as clearing function for order release planning is too optimistic, (i.e., it overestimates the output that can be expected). This is likely to produce releases that are too late, resulting in poor due-date performance. This consequence is confirmed by simulations. Selc- uk (2007) finds that an optimistic clearing function estimation (high expected output for given WIP) leads to low final product inventory, but increases duedate deviations and vice-versa (see the simulation results in Selc- uk, 2007, pp. 102–107). Asmundsson et al. (2009) report that in a high breakdown scenario the clearing functions that are obtained from simulation data by minimizing the sum of squared errors and then correcting the final segment of the outer linearization to account for average availability, ‘‘significantly overestimate the capability of the machine’’ (p. 153). Since this does not occur in the low breakdown scenario (Asmundsson et al., 2009, p. 153), future research must clarify to what extent the bias described in this paper can explain the results. The analysis presented here provides a sound theoretical basis for avoiding this by defining a clearing function that correctly models the properties of the manufacturing system.

6. Conclusions and research topics

5

2

4

0 0 2

2

k−

1

3 W

E [X k]

4

Clearing function models currently appear to be the only approach to address load-dependent lead times in order release planning that is at least close to practical applicability. The structure of optimization models based on clearing functions has been worked out extensively, but the definition of clearing functions and the estimation of their parameters are not as yet based on a sound theory. This paper describes a way to define and estimate a transient clearing function based on theoretical insights on the transient behaviour of work centres. It does not solve all issues that arise in this context. The main limitations that require further research are the following:

4 6 Ak

1

 The proposed approach must be worked out in detail, and tests

8 10

with empirical or simulated data must be performed.

0

 The proposed approach treats each period independently. The

Fig. 5. Plot of the assumed empirical clearing function. Due to the demanding computations, E[Xk] was set to the maximum output (= 5) beyond a certain point.

actual work input in a future period (its deviation from the expected or planned value) has no impact on the work input of

Table 1 Relative difference between the empirical and the corrected clearing function, defined as follows: (empirical value  corrected value)/empirical value.

Wk  1 = 0 Wk  1 = 1 Wk  1 = 2 Wk  1 = 3 Wk  1 = 4 Wk  1 = 5

Ak = 0

Ak = 1

Ak = 2

Ak = 3

Ak = 4

Ak = 5

Ak = 6

Ak = 7

Ak =8

Ak =9

Ak =10

Indeterminate 0.0732 0.17 0.248 0.303 0.343

0.0559 0.0909 0.164 0.224 0.266 0.295

0.0595 0.0967 0.154 0.201 0.233 0.254

0.0595 0.0969 0.143 0.179 0.202 0.217

0.056 0.0927 0.131 0.157 0.173 0.184

0.0828 0.111 0.138 0.155 0.165 0.177

0.0498 0.0829 0.107 0.12 0.127 0.148

0.0471 0.0776 0.0967 0.105 0.0998 0.125

0.0414 0.0703 0.0853 0.0905 0.0886 0.106

0.0345 0.0624 0.0744 0.098 0.0892 0.0904

0.0269 0.0547 0.0992 0.0815 0.0753 0.0785

406

 

H. Missbauer / Int. J. Production Economics 131 (2011) 399–406

the subsequent period(s). Since the stochastic nature of the work input is partly due to the uncertainty in the timing of the input (the total input over a longer time horizon is determined by order release), it can be expected that lower input in a certain period indicates a higher input in the next period(s). This autocorrelation is not considered here. The approach presented here deals only with the improvements of the clearing function. The optimization model incorporating these improvements remains to be formulated. A three-dimensional clearing function, possibly with case distinctions (e.g., overloaded/underloaded periods) might make linearization difficult if not impossible. Appropriate algorithms to perform the optimization must be found.

The desired result of this research program is a theoretically consistent way to handle load-dependent lead times in order release planning. The insights presented here indicate that this goal is within reach. References Abate, J., Whitt, W., 1987a. Transient behavior of the M/M/1 queue: starting at the origin. Queueing Systems: Theory and Applications 2 (1), 41–65. Abate, J., Whitt, W., 1987b. Transient behavior of regulated Brownian motion, II: non-zero initial conditions. Advances in Applied Probability 19 (3), 599–631. Abate, J., Whitt, W., 1988. Transient behavior of the M/M/1 queue via Laplace transforms. Advances in Applied Probability 20 (1), 145–178. Abate, J., Whitt, W., 1994. Transient behavior of the M/G/1 workload process. Operations Research 42 (4), 750–764. ¨ ¨ Andersson, H., Axsater, S., Jonsson, H., 1981. Hierarchical material requirements planning. International Journal of Production Research 19 (1), 45–57. Asmundsson, J.M., Rardin, R.L., Uzsoy, R., 2006. Tractable nonlinear production planning models for semiconductor wafer fabrication facilities. IEEE Transactions on Semiconductor Manufacturing 19, 95–111. Asmundsson, J., Rardin, R. L., Uzsoy, R., 2002. Tractable nonlinear capacity models for aggregate production planning. Working Paper, Laboratory for Extended Enterprises at Purdue, School of Industrial Engineering, Purdue University, West Lafayette, IN. Asmundsson, J.M., Rardin, R.L., Turkseven, C.H., Uzsoy, R., 2009. Production planning with resources subject to congestion. Naval Research Logistics 56, 142–157. Bergamaschi, D., Cigolini, R., Perona, M., Portoli, A., 1997. Order review and release strategies in a job shop environment: a review and a classification. International Journal of Production Research 35, 399–420. Bertrand, J.W.M., Wortmann, J.C., 1981. Production Control and Information Systems for Component-Manufacturing Shops. Elsevier, Amsterdam. Bertrand, J.W.M., Wortmann, J.C., Wijngaard, J., 1990. Production Control: a Structural and Design Oriented Approach. Elsevier, Amsterdam. Bertsimas, D., Sethuraman, J., 2002. From fluid relaxations to practical algorithms for job shop scheduling: the makespan objective. Mathematical Programming, Series A 92, 61–102.

Chen, H., Mandelbaum, A., 1994. Hierarchical modeling of stochastic networks, part I: fluid models. In: Yao, D.D. (Ed.), Stochastic Modeling and Analysis of Manufacturing Systems. Springer-Verlag, New York. Cohen, J.W., 1969. The Single Server Queue. North Holland, Amsterdam-London. Graves, S.C., 1986. A tactical planning model for a job shop. Operations Research 34 (4), 522–533. Hackman, S.T., Leachman, R.C., 1989. A general framework for modeling production. Management Science 35, 478–495. Hung, Y.F., Leachman, R.C., 1996. A production planning methodology for semiconductor manufacturing based on iterative simulation and linear programming calculations. IEEE Transactions on Semiconductor Manufacturing 9 (2), 257–269. Johnson, L.A., Montgomery, D.C., 1974. Operations Research in Production Planning, Scheduling and Inventory Control. John Wiley, New York. Karmarkar, U.S., 1989. Capacity loading and release planning with work-inprogress (WIP) and leadtimes. Journal of Manufacturing and Operations Management 2, 105–123. Kleinrock, L., 1976. Queueing Aystems. vol. II: Computer System Applications. John Wiley and Sons, New York-Sydney-Toronto. Land, M., 2004. Workload Control in Job Shops, Grasping the Tap. Labyrinth Publications, Ridderkerk. ¨ eine Neugestaltung von PPSMissbauer, H., 1998. Bestandsregelung als Basis fur Systemen. Physica, Heidelberg. Missbauer, H., 2002. Aggregate order release planning for time-varying demand. International Journal of Production Research 40 (3), 699–718. Missbauer, H., 2007. Meta-models of the transient behaviour of a work centre to improve clearing function models. Working Paper, Department of Information Systems, Production and Logistics Management, University of Innsbruck, Innsbruck, Austria. Missbauer, H., 2009. Models of the transient behaviour of production units to optimize the aggregate material flow. International Journal of Production Economics 118 (2), 387–397. Missbauer, H., Uzsoy, R., 2008. Optimization models of production planning problems. In: Kempf, K., Keskinocak, P., Uzsoy, R. (Eds.). Planning Production and Inventories in the Extended Enterprise: a State of the Art Handbook. Kluwer International Series in Operation Research and Management Science (forthcoming). ˜ o, G., 2003. Transient behavior of stochastic networks: application to Rian production planning with load-dependent lead times. PhD Thesis, School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, USA. Selc- uk, B., 2007. Dynamic performance of hierarchical planning systems: modeling and evaluation with dynamic planned lead times. PhD Thesis, Technische Universiteit Eindhoven, The Netherlands. Srinivasan, A., Carey, M., Morton, T.E., 1988. Resource pricing and aggregate scheduling in manufacturing systems. Unpublished paper, GSIA, CarnegieMellon University. Stevenson, M., Hendry, L.C., 2006. Aggregate load-oriented workload control: a review and a re-classification of a key approach. International Journal of Production Economics 104 (2), 676–693. Tatsiopoulos, I.P., Kingsman, B.P., 1983. Lead time management. European Journal of Operational Research 14, 351–358. Voss, S., Woodruff, D.L., 2006. Introduction to Computational Optimization Models for Production Planning in a Supply Chain, second ed Springer-Verlag, Berlin. Wight, O.W., 1974. Production and Inventory Management in the Computer Age. CBI Publishing Co., Boston. ¨ Zapfel, G., Missbauer, H., 1993. New concepts for production planning and control. European Journal of Operational Research 67, 297–320.