Raw material release time control for complex make-to-order products with stochastic processing times

ARTICLE IN PRESS

Int. J. Production Economics 103 (2006) 371–385 www.elsevier.com/locate/ijpe

Raw material release time control for complex make-to-order products with stochastic processing times Dong-Ping Song International Shipping & Logistics Group, The Business School, University of Plymouth, Plymouth PL4 8AA, UK Received 29 September 2003; accepted 14 July 2005 Available online 9 February 2006

Abstract This paper considers a make-to-order manufacturing system producing complex products with multiple resource constraints, multiple levels of product structures, and stochastic processing times. The objective is to find the optimal raw material release times by minimising the expected sum of work-in-progress holding costs, product earliness costs and product tardiness costs. An evolution strategy-based planning tool is developed, in which the performance measure is evaluated by stochastic discrete event simulation. The tool is applied to an illustrative example that includes wide, deep and complex product structures and to a case study using data obtained from an industrial company that manufactures capital goods. It is demonstrated that the solutions produced perform substantially better than those produced by the backward scheduling based on infinite capacity, and also significantly better than those produced by the same search method based on deterministic processing times. More benefit can be achieved if the system has a higher degree of uncertainty. r 2006 Elsevier B.V. All rights reserved. Keywords: Manufacturing; Complex product structure; Raw material release time; Evolution strategy; Genetic algorithms

1. Introduction For products that have complex product structures and operate on a make-to-order basis, materials are required to perform different assembly activities during different time periods. An assembly operation cannot be started until all (or most) of its subassemblies and components are ready. For example, a turbine-generator has components with long sequences of manufacturing processes, such as blades, followed by major assemblies such as rotors, casings and bearings. The final product is assembled and tested before delivery to site. It is therefore Tel.: +44 1752232442; fax: +44 1752232249.

E-mail address: [email protected].

desirable to design material release (or ready) times as late as possible, but not too late, in order to reduce work-in-progress (WIP) holding costs and product earliness and tardiness costs. This is a difficult problem in stochastic situations, where uncertainties in production processes are accumulating and interacting due to the precedence constraints, resource constraints and assembly coordination requirements. The importance of material release time planning in stochastic manufacturing systems has attracted much attention in the last two decades. Similar terminologies used in the literature include planned dispatch time (Yano, 1987; Gong et al., 1994), job release time (Liao, 1992; Song et al., 2001; Elmaghraby, 2001), material and job ready time (Hodgson et al., 1997),

0925-5273/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ijpe.2005.07.006

ARTICLE IN PRESS 372

D.-P. Song / Int. J. Production Economics 103 (2006) 371–385

activity delay (Buss and Rosenblatt, 1997), material release time (HomemdeMello et al., 1999; Hasan and Spearman, 1999), and activity start time (Elmaghraby et al., 2000). A lot of previous research work is focused on single-machine systems. For example, Liao (1992) considered the optimal release time for the arriving jobs in a single-machine system with stochastic processing times. Hodgson et al. (1997) investigated the problem of determining simultaneously material and job ready times and production sequence in a single facility and a flow-shop environment subject to the given service level constraints. Elmaghraby et al. (2000) extended the above work by including the objective of maximising an economic gain. A dynamic programming model is applied to a single facility to determine the optimal activity/job start times. Elmaghraby (2001) further applied dynamic programming to optimise a common release time for all jobs and the job sequence in a system with a single or multiple parallel machines. Several authors considered serial production lines. Yano (1987) and Gong et al. (1994) investigated the problem of determining optimal planned lead-times in a stochastic serial production system. The planned lead-times are used to determine the batch dispatch time at each stage. They pointed out that the primary reason for holding back a batch is that the cost of holding it at subsequent stages makes it uneconomical to dispatch it immediately to the next stage despite potential penalty costs for tardiness. HomemdeMello et al. (1999) considered setting job release times in a stochastic production flow line in which the sequence of jobs is fixed. A Monte Carlo simulation-based optimisation method was proposed to minimise an expected value objective function that contains terms for tardiness and flow time costs. Hasan and Spearman (1999) studied a serial production line operating in maketo-order mode. Perturbation analysis and stochastic approximation was applied to find the optimal material release times for a fixed sequence of jobs. Some work focused on other planning issues in multistage assembly systems. For example, Pongcharoen et al. (2002) applied genetic algorithms (GA) to schedule products with many stages of operations by focusing on creating operation sequences with deterministic processing times. Hicks (2004) developed a genetic algorithm tool to design the layout of manufacturing facilities for capital goods companies that produce complex products. Lagodimos et al. (2004) presented heuristic algo-

rithms to schedule operation sequences and lot quantities in a multistage fabrication shop. Little work has been reported on material release time control for systems with general product structures. Buss and Rosenblatt (1997) showed that delaying the start of an activity from its earliest start time could indeed increase the expected net present value of a project. However, their analytical treatment was limited to the Markov activity networks and implied no resource constraints. Song et al. (2001) applied perturbation analysis technique to optimise planned job release times in a stochastic assembly system with general product structures where operation sequences on resources are fixed. The majority of the above work was limited to simple systems such as single machine and serial production line. Some of them assumed the fixed operation sequences on resources. Little work has been done on material release time planning in general make-to-order systems with complex product structures, stochastic processing times and multiple resource constraints. The main reason is the extreme complexity of the problem. This paper describes the development of an evolution strategy (ES)-based planning tool to seek the ‘‘optimal’’ raw material release times. An eventdriven simulation model is embedded to evaluate the objective function, which is defined as the expected total cost that consists of WIP holding costs, product earliness costs and product tardiness costs. Random search methods are often the alternatives to solve difficult optimisation problems, e.g. Aytug et al. (2003) reviewed the applications of GA (ES) to solve production and operations management problems. The rest of the paper is organised as follows. The mathematical formulation of the problem is given in Section 2. An event-driven simulation model, which is used to evaluate the performance measure, is described in Section 3. The application of the ES method to the raw material release time optimisation problem is presented in Section 4. A numerical example that covers a variety of scenarios and a case study using data obtained from an industrial company are given in Sections 5 and 6, respectively, to demonstrate the effectiveness of the method. Use of the model is discussed in Section 7. Finally, conclusions are made in Section 8. 2. Problem formulation In make-to-order manufacturing systems, products are highly customised and may have quite

ARTICLE IN PRESS D.-P. Song / Int. J. Production Economics 103 (2006) 371–385

different product structures. The manufacturing of a product involves a variety of items such as raw materials (components), work-in-progress, subassemblies, assemblies and final products. For simplicity, each item is termed as a job and it has only one operation, i.e. item, job and operation are used interchangeably in this paper. Therefore, the jobs without a preceding job represent raw materials, the jobs without a successive job represent final products and others represent work-in-progress. The relationships of different jobs are subject to precedence constraints, resource constraints and assembly co-ordination requirements. The job precedence constraints are represented by a tree-type product structure, where one job may have more than one immediately preceding job but has at most one immediately succeeding job. Clearly, the raw materials represent the leaf nodes in the product structure and all other jobs lie in the path from a raw material to a final product. The resource constraints are caused by the fact that some jobs may be processed on the same resource but each resource can only perform one job at any time. Assembly co-ordination requires that all jobs that go to the same assembly should be synchronised. Otherwise, early arrived jobs have to wait for other jobs to perform the same assembly activity and thus incur extra holding costs. Uncertainties in process times interact through the assembly activities and their effect on later stage is cumulative. It makes the job arrive at each resource randomly. If two or more jobs compete for service, a priority rule is then used to select the next one. Since each job lies in the path from a raw material to a final product, the production process can be controlled and co-ordinated by carefully assigning release times to the raw materials. The raw materials are not allowed to be processed prior to the assigned release times. The objective is to find the optimal raw material release times that minimise the expected sum of holding costs and tardiness cost. Notation is defined as follows: J F

L

Bi

the entire job set the set of the first job in each branch in each product structure, i.e. the set of jobs that have no preceding job. It is also termed the raw material set the set of the last job in each product structure, i.e. the set of jobs that have no successive job. It is also termed the product set the set of jobs which immediately precede (before) the job i in the product structure.

f(i)

r(i) j(i)

si xi ai ci di hi pi

373

This represents the assembly structural precedence among jobs the job that immediately follows the job i in the product structure. This also represents the structural precedence among jobs the resource on which the job i is processed the job which immediately precedes the job i on the resource r(i). This reflects the operational precedence among jobs that is the result of processing at each resource the planned release (or start) time of raw material iAF, which is a decision variable the processing time of job iAJ, which is a known random variable the actual processing start time of job iAJ the processing completion time of job iAJ the due date of product iAL the unit time holding cost for job iAJ the unit time tardiness cost (or penalty) for product iAL

For a given set of raw material release times, the evolution of the production process can be described by {ai} and {ci} subject to the following constraints: ai ¼ maxðfsi g [ fcjðiÞ gÞ;

i 2 F;

ai ¼ maxðfcjðiÞ g [ fcj : j 2 Bi gÞ; c i ¼ a i þ xi ;

i 2 J,

(1) i 2 J\F ;

(2) (3)

where J\F is the difference set between J and F. On the right-hand side (RHS) of Eq. (1), the first term represents that the raw material i cannot start processing prior to the planned release time si; and the second term describes that the raw material i cannot start processing on the resource r(i) until the job j(i) is completed. The second term on the RHS of (2) represents that a WIP job or a final product job cannot start processing until all its subassembly jobs are finished. Let s be a vector composed of all raw material release times. The performance measure is defined by hX V ðsÞ ¼ E h ða
ci Þ i2J\L i f ðiÞ X þ h maxðd i
ci ; 0Þ i2L i i X þ p maxðc
d ; 0Þ . ð4Þ i i i2L i Three terms on the RHS of (4) represent the WIP holding cost, product earliness cost and product

ARTICLE IN PRESS 374

D.-P. Song / Int. J. Production Economics 103 (2006) 371–385

tardiness cost, respectively. The problem is to find the optimal s by minimising the cost function V(s). Two points need to be emphasised. Firstly, the actual start time of a raw material may be later than its planned release time, which is reflected in Eq. (1). Secondly, the job sequence on each resource is undetermined because of the randomness of processing times. Due to the extreme complexity of the system, it is difficult to obtain the explicit form of the performance measure with respect to the decision variables {si, iAF}. This paper presents an ES method to seek the ‘‘optimal’’ solution. An event-driven simulation model is used to evaluate the performance measure, which is described in the next section first.

There are two types of events that enable system state transitions, namely the release of a raw material and the completion of a job. The material release event triggers a transition of the raw material state from ‘‘not ready’’ to ‘‘waiting’’ if the corresponding machine is busy, or to ‘‘processing’’ if the machine is idle. In the later case, the corresponding machine’s state changes to ‘‘busy’’. The job completion event may trigger the following state transitions:

 

 3. An event-driven simulation model

 For stochastic complex systems, discrete eventdriven simulation is often used to estimate a performance measure by performing averaging over multiple samples or replications. This section outlines an event-driven simulation model by identifying the system states, events and state transitions and illustrating how they evolve over time. The system state can be represented by the machine state and job state. Each machine (i.e. resource) has two states, busy or idle. Each raw material job has four states:

   

   

Not ready—the current time is less than its release time; Waiting—the current time is greater than its release time but the corresponding machine is busy; Processing—the raw material is being processed on a machine; Completion—the raw material has been processed on a machine. For any other job apart from raw materials, it also has four states with slightly different definitions: Not ready—at least one of its preceding jobs has not been finished by the current time; Waiting—all its preceding jobs have been finished by the current time but the corresponding machine is busy; Processing—the job is being processed on a machine; Completion—the job has been processed on a machine.

From ‘‘processing’’ to ‘‘completion’’ for the current job; From ‘‘not ready’’ to ‘‘waiting/processing’’ for the immediately successive job if the immediately successive job has no unfinished preceding jobs and the corresponding machine is busy/idle; From ‘‘waiting’’ to ‘‘processing’’ for the next job that is waiting for the current machine; From ‘‘busy’’ to ‘‘idle’’ for the current machine if no more job is in its queue.

The evolution of the system state and the active event list in the simulation is briefly described as follows. Suppose that the raw material release times {si, iAF} are given. Initially, all the raw material release events are put into the active event list. The first occurring event in the active event list is the one with the earliest release time (say sj). This event triggers a transition of the system state at time sj as follows: from ‘‘not ready’’ to ‘‘processing’’ for the raw material j and from ‘‘idle’’ to ‘‘busy’’ for the machine r(j). In other words, the raw material with the earliest release time will be processed exactly at its planned release time. The actual processing time of raw material j is generated following a known probability distribution. The active event list is updated by removing this raw material release event and adding its completion event. Then, search for the next occurring event in the active event list. It could be another raw material release event or the completion event. The system state is updated according to the type of the next event. For example, if it triggers a transition to the state ‘‘processing’’ and there is more than one job waiting for the machine, a specified priority rule is used to select the next job. The processing time for the selected job is generated from a known probability distribution. The active event list is updated by removing the occurred event and adding the new completion event. The process keeps progressing until all jobs are finished, i.e. the active event list becomes empty.

ARTICLE IN PRESS D.-P. Song / Int. J. Production Economics 103 (2006) 371–385

c1

5 Machine: jobs

1

3

1000: 1, 2, 3 1001: 4, 5

2

(i) The first occurring event in the active event list is the release of raw material 1 at time t ¼ s1 ¼ 0:5, which triggers the transition of the system state to {processing, not ready, not ready, not ready, not ready, busy, idle}. The active event list is updated to {release of raw material 2, release of raw material 3, comple-

c3 job 3

a2

a3

c4

c5

job 4

1001

0

c2

job 2

a1

Fig. 1. An illustrative example.

To have a clearer view of the above procedure, a simple example is given to illustrate how the system state and the active event list evolve over time. In Fig. 1, a single product is produced using two machines (numbered as 1000 and 1001). There are total 5 jobs including 3 raw materials (numbered as 1, 2, and 3) and 2 assembly operations (numbered as 4 and 5). Three raw materials are required to be processed on machine 1000 and two assembly operations are to be processed on machine 1001. Now assuming that the decision variables (i.e. the planned release times for raw materials) are made as follows: s1 ¼ 0:5, s2 ¼ 1:0 and s3 ¼ 2:0. Suppose that a sample of random processing times for five jobs has been generated such as: x1 ¼ 0:6, x2 ¼ 0:4, x3 ¼ 0:5, x4 ¼ 0:8, and x5 ¼ 0:5. Then, the Gantt chart (describing the start, duration and completion of all operations) of the actual processing on two machines can be determined by the event-driven simulation model and is shown in Fig. 2, where the length of each bar represents the corresponding job’s processing time (i.e. xi). The following part describes the details about how the simulation model determines the actual processing on machines. The system state is given by the set {job 1’s state, job 2’s state, job 3’s state, job 4’s state, job 5’s state, machine 1000’s state, machine 1001’s state}. There are total eight events including three raw material release events and five job completion events. Initially at time t ¼ 0, the active event list consists of all the raw material release events, i.e. {release of raw material 1, release of raw material 2, release of raw material 3}. The initial system state is {not ready, not ready, not ready, not ready, not ready, idle, idle}. The system state and the active event list are evolving as follows:

job 1

1000 machine

4

375

a4 s1

s2

0.5

1

job 5 a5

s3 1.5 t

2

2.5

3

Fig. 2. The Gantt chart of job processing on machines over time t.

(ii)

(iii)

(iv)

(v)

(vi)

tion of job 1}, obtained by removing the ‘‘release of raw material 1’’ event and adding the ‘‘completion of job 1’’ event. The next occurring event in the active event list is the release of raw material 2 at time t ¼ s2 ¼ 1. This event changes job 2’s state from ‘‘not ready’’ to ‘‘waiting’’ and the system state becomes {processing, waiting, not ready, not ready, not ready, busy, idle}. The active event list is updated to {release of raw material 3, completion of job 1}. The next event is the completion of job 1 occurring at t ¼ c1 ¼ 1:1. Machine 1000 starts processing job 2 immediately after finishing job 1 because job 2 is waiting. The system state becomes {completion, processing, not ready, not ready, not ready, busy, idle} and the active event list is updated to {release of raw material 3, completion of job 2}. The following event is the completion of job 2 occurring at t ¼ c2 ¼ 1:5. At this time, both raw materials 1 and 2 have completed and are ready to assembly on machine 1001. On the other hand, machine 1000 becomes idle because raw material 3 has not been released. Thus, the system state is changed to {completion, completion, not ready, processing, not ready, idle, busy}. The active event list becomes {release of raw material 3, completion of job 4}. The next event is the release of raw material 3 at time t ¼ s3 ¼ 2, which triggers the transition of the system state to {completion, completion, processing, processing, not ready, busy, busy}. The active event list becomes {completion of job 4, completion of job 3}. The next event is the completion of job 4 occurring at t ¼ c4 ¼ 2:3. The system state becomes {completion, completion, processing, completion, not ready, busy, idle} and the active event list is updated to {completion of job 3}.

ARTICLE IN PRESS D.-P. Song / Int. J. Production Economics 103 (2006) 371–385

376

(vii) There is only one event in the active event list, which occurs at t ¼ c3 ¼ 2:5. The system state is changed to {completion, completion, completion, completion, processing, idle, busy} and the active event list is updated to {completion of job 5}. (viii) The last event, the completion of job 5, occurs at t ¼ c5 ¼ 3, the system state becomes {completion, completion, completion, completion, completion, idle, idle} and the active event list becomes empty. It should be pointed out that due to the randomness of processing times the occurring sequence of the events may be different for different samples. For example, if the processing time of job 1 takes more than 1.5 (which implies that c1 4s2 and c1 4s3 ), then both job 2 and job 3 have already been released and are waiting for machine 1000 at time t ¼ c1 , in which case a priority rule should be employed to select the next one. Clearly, the simulation model yields {ai, ci, iAJ} dynamically and they satisfy Eqs. (1)–(3). By repeating the procedure K times for the fixed decision vector s, the cost defined in (4) can be estimated by averaging over K replications. In the next section, an ES method is applied to optimise the raw material release times {si, iAF}. 4. ES method Similar to GA, ES is a random search method that takes over the principle of biological evolution for the optimisation of technical systems (Back, 1996; Aytug et al., 2003). ES uses real variable to represent a solution and is thus more suitable for parameter optimisation. In the system under consideration a solution has been represented by a nonnegative vector s. Fig. 3 illustrates the ES procedure employed in this paper. The ES algorithm involves Initialisation

Selection

Recombination

No Mutation Optimal solution

Yes

Termination Adjustment Reproduction

Evaluation

Fig. 3. The evolution strategy procedure.

quite a few control parameters, which are introduced in the following sub-sections. 4.1. Initialisation In each generation, ES manipulates a set of individuals (solutions) that is called the parent generation, denoted by PP. Its size is denoted by jPPj. The initial parent population can be randomly or specifically generated. In this paper, for each individual in the initial parent population, the release time of each raw material is generated uniformly between zero and the product’s due date. Preliminary experiments show that using heuristic solutions (e.g. incorporating knowledge of mean processing times and product structures) in the initial population can improve the search speed, but does not guarantee a better solution. Let n denote the generation. 4.2. Selection and recombination From the parent population, two individuals are selected randomly, which will be used to generate a descendant. The offspring population is denoted by PO subject to jPOj4jPPj. It is assumed that all individuals in the parent population have the same probability to be selected. Other mechanism can also be applied, e.g. the parents with a higher fitness value (i.e. lower cost) have a higher probability to be chosen. However, preliminary experiments show there is no significant improvement using this fitness value-based probability selection rule. This may be explained by the fact that in the reproduction procedure (see later sub-sections) the jPPj best individuals have already been selected from total jPOj offspring. In addition, multiple parents could be selected to produce a descendant. This reflects the multisexual mechanism in the nature. Our tests show that multisexual mechanism can improve the convergence rate but not much. A descendant is generated by randomly copying the elements from two selected parents. That is, each gene of the descendant is randomly obtained from the corresponding gene in one of its parents. The number of genes remains to be jFj, where jFj represents the number of elements in the set F. This operator is the same as crossover in GA. An example of global uniform recombination for two parents at generation n is shown in Fig. 4, where si(n,j) denotes the planned release time for raw material i of the jth individual in the parent population at the nth generation.

ARTICLE IN PRESS D.-P. Song / Int. J. Production Economics 103 (2006) 371–385

377

Parent 1

s1(n,1)

s2(n,1)

s3(n,1)

... ...

s|F|-2(n,1)

s|F|-1(n,1)

s|F|(n,1)

Parent 2

s1(n,2)

s2(n,2)

s3(n,2)

... ...

s|F|-2(n,2)

s|F|-1(n,2)

s|F|(n,2)

Descendant

s1(n,2)

s2(n,1)

s3(n,1)

... ...

s|F|-2(n,2)

s|F|-1(n,1)

s|F|(n,2)

Fig. 4. Global uniform recombination.

4.3. Mutation and adjustment For each descendant, its genotype differs slightly from that of its parents. The deviations refer to individual genes, and are random and independent of each other. For example, for a descendant s(n,k), a Normal random vector z is generated such that its element zi  N(0,sn2), where z has the same size as s and sn2 is the variance for mutation at the nth generation. Then, the mutation operator is done by setting s(n,k) ¼ s(n,k)+z. After mutation, each descendant s(n,k) is adjusted by s(n,k) ¼ max(s(n,k), 0) to avoid negative release times. 4.4. Evaluation For each fixed descendant s(n,k), it is required to run K replications using the event-driven simulation model presented in the Section 3. The performance measure (4) can be approximated by averaging the costs over K replications. By the Law of Large Numbers, the approximated cost converges to (4) almost surely as K tends to infinity. However, it should be pointed that the total number of runs at each generation is the product of the offspring population size and the replication number, i.e. jPOjK. Take into account the fact that multiple generations (say N) are required, this gives rise to total jPOjK N runs to be required in the ES procedure. Therefore, K usually cannot be very large due to the computational limit. 4.5. Reproduction Reproduction is a process in which the best jPPj individuals are selected from the offspring population and they will form the parent population at the next generation. The remaining offspring are

simply discarded. Another method often used in this procedure is to select the best jPPj individuals from the current parent parents and the offspring population. That means, a parent is allowed to be alive for multiple generations. 4.6. Termination criteria Termination rules or stop criteria should be defined. Let nc denote the number of consecutive generations in which no improvement is achieved for the cost function. Let Nc be the maximum number of nc to be allowed in the procedure. Let N be the maximum number of generations. In this paper, it is assumed that the optimisation procedure is terminated if nc 4N c or n4N. The best solution up to now, denoted by s*, is then returned as the ‘‘optimal’’ solution. The standard deviation for mutation sn is reduced by a constant factor a if nc is greater than a predetermined number (denoted by Ns). In other words, Ns is a threshold value to trigger the reduction of the standard deviation sn. It has been observed that this can improve the search quality significantly compared with the fixed standard deviation. At the end of the ES procedure, the ‘‘optimal’’ solution s* provides schedulers or managers with the information about when each raw material should be released in order to minimise the cost function defined in (4). However, it should be pointed out that multiple high-quality local optima may exist due to the combinatorial nature of the problem. 5. Computational experiments This section gives an illustrative example, which has nine different product structures including wide

ARTICLE IN PRESS D.-P. Song / Int. J. Production Economics 103 (2006) 371–385

378

(ac), deep (df) and complex (gi) structures (See Fig. 5). Each product structure corresponds to one product and therefore nine products are to be produced. It involves total 91 jobs including 64 raw materials and 27 assembly operations. The jobs are numbered consecutively from 1 to 91 as shown in Fig. 5. The assembly operations are processed on four assembly machines (numbered by 1000–1003) and the raw materials are processed on eight nonassembly machines (numbered by 2000–2007). The relationships between jobs and machines are given in Table 1. These relationships indicate which job should be processed on which machine, but do not specify the processing sequences. The mean processing time of each job is uniformly obtained from the range of 1–10 days. The number of replications for evaluating each solution s is 100, i.e. K ¼ 100, which is used in stochastic discrete event simulation. Assuming that the holding cost coefficient hi equals the total time that has been spent on this job and on all preceding jobs multiplied by a constant 0.01. This 1 2

reflects that holding costs increase with added value. It is assumed that the product tardiness penalty coefficient is twice the corresponding product earliness penalty coefficient, i.e. pi ¼ 2hi , i 2 L. All products have the same due date 180 days. The simulation program is written in C-language and run in Pentium IV personal computers. A range of scenarios with different uncertainty levels are tested. In the literature, normal distribution and exponential distribution are commonly used to model the uncertainty of processing times. Four levels of uncertainty are considered in this paper. In the first three levels, the processing times follow normal distributions with left truncated at zero; and the coefficient of variation v (defined as the standard deviation divided by the mean) takes three values 0.1, 0.2 and 0.5, respectively. In the forth level of uncertainty, the processing times follow Exponential distributions determined by their means. Clearly, the variance of the processing time, which indicates the degree of uncertainty, increases from levels 1 to 4.

6

3

4

5

7

8

15 13

...

16

14

17

18

8 components (a)

29

...

30

31

16 components

(b)

(c) 46 47 37

50

49

32 38 33

39

51

52

34 41

40

53

54

36

35

42 44 (d)

43

55

45

57

(e)

56 58

(f) 80

59 61

60 62 (g)

63

67

64

65

68 66

72

73

69 74

75

84

71

70 76

82

81

77

78

85

79

(h)

(i) Fig. 5. Various product structures.

83 86

89 90

87 91

48

88

ARTICLE IN PRESS D.-P. Song / Int. J. Production Economics 103 (2006) 371–385 Table 1 Relationships between jobs and machines

Assembly machine

Machine

Jobs

1000

1, 15, 33, 55, 68, 69, 81, 89 38, 40, 42, 49, 61 6, 51, 59, 60, 84 32, 37, 46, 47, 53, 67, 71, 80, 82

1001 1002 1003 Non-assembly machine

2000 2001 2002 2003 2004 2005 2006 2007

16, 24, 27, 48, 52, 76, 86 2, 9, 18, 22, 65, 73, 83 3, 10, 11, 13, 36, 44, 54, 78, 90 30, 39, 43, 56, 62, 66, 77, 87 5, 23, 25, 29, 34, 45, 64, 70 4, 8, 14, 19, 20, 35, 50, 57, 63, 74, 79, 85, 88 7, 17, 21, 31, 41, 58, 91 12, 26, 28, 72, 75

In the simulation model, a priority rule should be specified for each machine. Before introducing the priority rules, two definitions are given, which are slightly different from those used in job shops.





The remaining expected processing time of operation i is defined as the cumulative mean processing times along the path from the operation i to its final P product. Mathematically, it is equal to x¯ i þ nk¼1 x¯ ik , where x¯ i represents the mean processing time of operation i and (i, i1, i2, y, in) is the operation list along the path from i to its final product satisfying i1 ¼ f ðiÞ and ikþ1 ¼ f ðik Þ for k ¼ 1, 2, y, n
1, and in 2 L. The P operation i’s due date is defined as d in
nk¼1 x¯ ik , where (i, i1, i2, y, in) is the operation list along the path from i to its final product.

In this paper, four priority rules are examined separately. First-come first-served (FCFS) selects the operation that has the earliest arriving time. Earliest operation due date (EOD) selects the operation that has the smallest operation due date. This rule is tailored from the well-known EDD rule. Shortest processing time (SPT) selects the operation that has the shortest mean processing time at that machine. Least remaining work (LRW) selects the operation that has the least remaining expected processing time.

379

To investigate the effectiveness of the proposed method, two other methods are implemented to make comparison. The first is the traditional backward scheduling based on mean data and infinity capacity (often used in MRP), denoted by BSIC. The second is the same ES method by treating the operation processing times as constants (based on their means), where the performance measure is evaluated by a single sample since the system is regarded as deterministic. This search method is denoted by ES-det. The solution obtained by BSIC or ES-det is then re-evaluated in stochastic situation by averaging over K replications. Intuitively, the BSIC designs raw material release times without considering resource capacity constraints and uncertainties in processing times, while ES-det takes into account the resource constraints but ignores the uncertainties. The ES includes a number of control parameters such as jPPj (parent population size), jPOj (offspring population size), s0 (initial standard deviation for mutation), a (reducing factor of the standard deviation), and Ns (the threshold value to trigger standard deviation reduction). The performance of ES may be affected by the control parameters. They are often chosen based on pilot runs (Aytug et al., 2003). This paper uses ES-det to perform pilot runs and the results show that the choice of parameters is related to the product due dates, the processing time, the number of raw materials and the total number of operations. For example, if jPOj is too small compared with the number of raw materials, then the search procedure may terminate early. This is possibly due to the lack of diversity of solutions at each generation. The parent population is the elite part of the offspring population, therefore jPPj should not be too small (lacking the diversity) or too large (not representing the elite part) compared to jPOj. If s0 is too small compared with the average processing time, then the search procedure appears very slow. This is possibly because the changes from parent generation to offspring generation are too small to make significant improvement. Based on the pilot runs, the following recommendations have been made: the offspring population size should be close to or larger than the number of raw materials; the parent population size is about 1/101/2 of the offspring population size; the initial standard deviation (s0) can be selected from the interval (d/jJj/2, 5d/jJj), where d is the largest product due date; the reducing factor of the standard deviation (a) should not be too small (e.g. X0.5); the threshold value to

ARTICLE IN PRESS D.-P. Song / Int. J. Production Economics 103 (2006) 371–385

380

Therefore, to complete each block (i.e. one priority rule) in Table 3, it requires total CPU time 4  32(10+680)/3600E25 h. It can be seen from Table 3 that ES-sto is substantially better than BSIC in all cases with the percentage of cost reduction between 63% and 80%. This indicates that designing raw material release times using BSIC is inappropriate for complex products. Interestingly, this result is almost independent of the degree of uncertainty. In other words, even for the systems with very low degree of uncertainty, the BSIC still produces very bad plans. On the other hand, ES-sto is also significantly better than ES-det with the percentage of cost reduction between 7% and 38%, where larger reduction can be achieved in situations with higher degrees of uncertainty (e.g. Normal distribution with v ¼ 0:5 and exponential distribution). This implies that using ES-det designing raw material release times appears to be satisfactory if the systems have low degrees of uncertainty. To an extreme situation, if the system is deterministic, then ES-sto and ES-det become the same. Table 3 also shows that on the whole the cost increases as the uncertainty degree increases for all four-priority rules. This is in agreement with the intuitions that higher uncertainty degree leads to more WIP waiting time and less reliable product delivery. In addition, Table 3 also shows that the costs corresponding to different priority rules are quite close, which means that the effects of different

trigger standard deviation reduction (Ns) can be 1/101/5 of Nc. In this experiment, two levels for each control parameter are experimented (see Table 2). The parameters of the termination criteria in ES are set by N c ¼ 50 and N ¼ 20000= j PO j. Since the number of runs of the ES equals the product of the offspring population size, the number of generations and the number of replications, this product is held constant (i.e. j PO j N K ¼ 20 000). The total number of combinations of parameters and levels is 25 ¼ 32. For each given scenario (with fixed uncertainty level) using a specified priority rule, the performance of ES is estimated by averaging the ‘‘optimal’’ costs over 32 combinations. The average costs under three methods for four levels of uncertainty and four priority rules are given in Table 3, where ES-sto denotes the ES method described in Section 4 that uses K replications to estimate the expected cost. The last two columns in Table 3 express the percentage of the cost reduction achieved by ES-sto compared to BSIC and ES-det. The CPU time per run is less than 0.06 s for BSIC, about 10 s for ES-det and 680 s for ES-sto. Table 2 Control parameters and levels for ES method Parameters

jPPj

jPOj

s0

a

Ns

Level 1 Level 2

10 20

60 100

1 10

0.5 0.9

5 10

Table 3 Costs in different scenarios with different priority rules BSIC

ES-det

ES-sto

% reduction (BSIC)

% reduction (ES-det)

FCFC

v ¼ 0:1 v ¼ 0:2 v ¼ 0:5 Exp

438 443 484 548

94 108 170 279

87 95 122 186

80 79 75 66

8 12 28 33

EOD

v ¼ 0:1 v ¼ 0:2 v ¼ 0:5 Exp

390 397 442 525

95 114 178 301

87 94 123 190

78 76 72 64

8 17 31 37

SPT

v ¼ 0:1 v ¼ 0:2 v ¼ 0:5 Exp

417 419 462 507

94 110 167 274

88 95 122 187

79 77 74 63

7 14 27 32

LRW

v ¼ 0:1 v ¼ 0:2 v ¼ 0:5 Exp

457 452 506 557

111 128 188 299

99 97 126 191

78 79 75 66

11 24 33 36

ARTICLE IN PRESS D.-P. Song / Int. J. Production Economics 103 (2006) 371–385 Table 4 Optimal raw material release times for different product structures Product 6 (flat structure)

Product 46 (deep structure)

Product 67 (complex structure)

r.m. i

si

r.m. i

si

r.m. i

si

7 8 9 10 11 12 13 14 —

163.55 122.65 163.64 166.95 163.15 160.18 145.38 154.90 —

48 50 52 54 56 57 58 — —

163.51 143.76 151.60 133.78 119.59 93.40 126.05 — —

70 72 73 74 75 76 77 78 79

157.07 156.65 144.67 127.15 156.92 122.67 118.48 121.52 104.15

priority rules are not very significant. This may be explained by the fact that complex products require jobs to be done during different time periods and therefore the opportunity of competition for the same resource is relatively lower than that in a job shop. It is interesting to investigate the relationship between the ‘‘optimal’’ release times achieved by the ES procedure and the system parameters such as due dates, product structures and resource capacity. As an example, Table 4 gives the ‘‘optimal’’ raw material (r.m.) release times for products 6, 46 and 67, which correspond to flat, deep and complex product structures in Fig. 5(b, f, h). The results in Table 4 are obtained by the ES-sto method (with jPPj ¼ 10, jPOj ¼ 100, s0 ¼ 10, a ¼ 0:9, and N s ¼ 10) for the case v ¼ 0:2 using the FCFS priority rule. Table 4 indicates that the earliest release time and the latest release time of raw materials are 122.65 and 166.95 for product 6; are 93.40 and 163.51 for product 46; are 104.15 and 157.07 for product 67. This reflects the effects of different product structures, e.g. flat structure tends to require the release of raw materials in a relatively narrow range of time so that they can reach the assembly operation at similar time; deep structure tends to release raw materials in a relatively wide range of time and roughly demonstrates a hierarchical structure (i.e. raw materials for later assembly operations are released later); while complex structure releases raw materials in a middle range of time compared with flat structure and deep structure. It should be pointed out that the release times of raw materials also depend on their processing times, product due

381

dates and the workloads on their corresponding machines. Generally speaking, raw materials with larger processing times, tighter product due dates or using machines with higher workloads should be released earlier. For example, in Table 4 the raw material 57 is released at time 93.40, which is significantly earlier than other raw materials. This can be partly explained by the facts that raw material 57 locates in the bottom of a deep product structure and is performed on machine 2005 that has the largest workload (see Table 1). Similarly, raw materials 8 and 14 are released much earlier than most of other raw materials of product 6, because they are also processed on the heavily loaded machine 2005. However, the explicit relationship between release times and the system parameters is almost surely problem-dependent. Nevertheless, the tool developed in this paper provides the flexibility for planners or managers to investigate various ‘‘what-if’’ scenarios.

6. A case study In this section, the methods are applied to a case study using data obtained from a capital goods company manufacturing steam turbines (Song, 2001). The case study has three products. The product structure information is given in Table 5, where the width of a product structure is defined as the number of raw materials and the depth of a product structure is defined as the maximum number of operations from a raw material to the final product. As an example, the detailed product structure of product No 3 (coded by 228) is shown in Fig. 6. The system involves total 239 operations (including 47 raw materials and 13 assembly operations) performed on 17 non-assembly machines and 2 assembly machines. The delivery date is fixed at 900 days from zero start time for all products. The mean operation times and the resource capacities are obtained from the company. The cost Table 5 Product structure information in the case study Product

Width

Depth

Assembly operations

Total operations

No 1 No 2 No 3

4 34 9

9 8 20

3 3 7

17 109 113

ARTICLE IN PRESS D.-P. Song / Int. J. Production Economics 103 (2006) 371–385

382

228

237

229

230

234

231

226 :15

246

239

243

247

240 :11

7 opers

240:1

241:1

244:9

9 opers

244:1

248:8

8 opers

. ..

11 opers

241:7

. ..

236:1

242

. ..

235:1

16 opers

238

. ..

233:1

10 opers

. ..

232:1

236 :16

235 :10

. ..

12 opers

. ..

226:1

12 opers

. ..

. ..

5 opers

233 :12

232 :12

245

248:1

Fig. 6. Product structure for the product 228 (Song, 2001).

coefficients are set in the same way as that in Section 5. The cost unit is GB(£)  1000. To investigate the effects of uncertainty, three methods (BSIC, ES-det and ES-sto) are applied to the system with varying degrees of uncertainty as such in Section 5. Again, four priority rules are included. Preliminary experiments show that if the initial standard deviation for mutation (i.e. s0) is too small compared to the due date, it may fail to find a good solution. Based on the results in Song (2001) and the experiments in Section 5, the control parameters for ES-det and ES-sto are selected as follows: j PP j¼ 10, j PO j¼ 100, s0 ¼ 10, a ¼ 0:9, and N s ¼ 10. The termination criteria are the same as those in Section 5. The costs under three methods and the cost reduction percentages achieved by ESsto are given in Table 6. The CPU time per run is about 34 s for ES-det and 2980 s for ES-sto in Pentium IV personal computers. Table 6 verifies the effectiveness of the ES-sto method. Compared with BSIC, the ES-sto can reduce the cost by around 80%. Compared with ES-det, significant improvement can be achieved in the situations with higher degrees of uncertainty, e.g. the ES-sto reduces the cost by 15–25% in the third and forth levels of uncertainty. For a given scenario with the fixed degree of uncertainty, the reduction percentages are quire similar for different priority rules.

7. Discussion The main purpose of this paper is to describe an ES-based planning tool and show that it works. The tool can help a manager to determine when the raw materials should be ready in manufacturing systems producing complex products with deterministic/ stochastic processing times. Suppose a manager has a number of products with known product structures, due dates, resource (machine) capacity, and specified priority rules for each machine, then a single run using a set of pre-selected control parameters can produce an ‘‘optimal’’ solution. For the example shown in the experimentation section, it takes 10 s to obtain an ‘‘optimal’’ solution if the processing time is deterministic and takes 680 s if the processing times follow known probability distributions. 7.1. Objective functions This paper defines a cost function that takes into account the costs of holding stock (workin-progress and products) and the costs of tardiness. Sometimes, the planner seeks to minimise the makespan or some function of tardiness (e.g. service level). By replacing the cost in (4) with the specified performance measure, the tool can be applied to

ARTICLE IN PRESS D.-P. Song / Int. J. Production Economics 103 (2006) 371–385

383

Table 6 Costs and reduction percentages in the case study BSIC

ES-det

ES-sto

% reduction (BSIC)

% reduction (ES-det)

FCFC

v ¼ 0:1 v ¼ 0:2 v ¼ 0:5 Exp

2328 2330 2350 2367

461 485 577 756

424 449 489 583

82 81 79 75

8 7 15 23

EOD

v ¼ 0:1 v ¼ 0:2 v ¼ 0:5 Exp

2330 2332 2350 2357

491 513 601 767

452 448 490 574

81 81 79 76

8 13 19 25

SPT

v ¼ 0:1 v ¼ 0:2 v ¼ 0:5 Exp

2335 2337 2355 2364

476 498 585 756

421 451 482 583

82 81 80 75

11 10 18 23

LRW

v ¼ 0:1 v ¼ 0:2 v ¼ 0:5 Exp

2335 2353 2397 2454

467 482 590 786

447 460 504 590

81 80 79 76

4 5 15 25

solve such problems without changing the algorithmic approach. This is a common attribute of random search methods, in which the simulation model and the search algorithm do not depend on the form of performance measure. 7.2. Priority rules Four commonly used priority rules are tested in this paper. In practice, there are many other priority rules and different machines may adopt different rules. The tool is flexible in terms of adopting priority rules. The experimentation and the case study show that although different priority rules lead to different solutions, their effects are not very significant. However, this result may not be generalised to general cases. 7.3. Speed and quality The ES-sto method requires a population of solutions, multiple generations, and multiple replications for evaluating each solution; it is therefore time-consuming for very large systems, especially when multiple scenarios (e.g. different priority rules, different resource capacities) are involved. How to improve the search speed is interesting. The most intuitive way is to start with a better set of initial solutions. For example, a two-stage ES method can be easily constructed as follows: at the first stage, the ES-det is used to produce a set of solutions for

the deterministic system; at the second stage, the ES-sto is used to produce an ‘‘optimal’’ solution for the stochastic system, in which the initial solutions are obtained from the results of the first stage. Our preliminary experiments show that such combined method could save up to half of the computational time compared with the ES-sto using random initial solutions. However, it is noteworthy that the search speed of ES is usually quick at the early stage (when the solutions are not good) but will slow down greatly at later stage. When multiple scenarios are required, parallel computing is one way to run multiple experiments simultaneously. In addition, reducing the number of samples or introducing tighter termination criteria can also reduce the computational time, but this may sacrifice the quality of solutions. In the absence of the capacity to perform the required ES-sto application, the ES-det may be applied by incorporating uncertainty into deterministic processing times (e.g. the expected processing time plus the standard deviation). Particularly, if the system’s uncertainty level is low, ES-det is recommended. 7.4. Optimality and lower bounds It has been well recognised that the random search methods including ES often find a sub-optimal solution instead of the global optimum. Due to the extreme complexity of the problem, there is no analytical result available. A lower bound of the

ARTICLE IN PRESS 384

D.-P. Song / Int. J. Production Economics 103 (2006) 371–385

cost function is useful to judge the quality of a solution. Let o denote a sample process, W(s, o) denote the sample cost function and s*(o) is the optimal raw material release times to minimise W(s, o). Then, it follows V(s) ¼ minsEW(s, o)X E minsW(s, o) ¼ EW(s*(o), o). Therefore, if the optimal solution s*o) can be found for any given sample process o, then this inequality provides a lower bound for the cost function V(s). Clearly, finding s*(o) is equivalent to solving an optimisation problem with deterministic processing times. If the operation sequences on each resource are fixed, Song (2001) presented an exact algorithm to find s*(o) under appropriate assumptions on costs. However, for the problem under consideration, it is difficult and no tight lower bounds have been reported. Therefore, it remains an open question to identify the closeness of a solution to the optimal one. 7.5. Validity to general cases The experimentation uses nine different product structures with mean processing times generated randomly from [1,10] and actual processing times follow normal or exponential distributions. In the case study, both the product structures and the mean processing times were obtained from an industrial company. Clearly the tool is applicable to other product mix and other processing time distribution and can identify the scopes of improvement. Note that ES-sto considers stochastic nature and other constraints as well while BSIC and ES-det only take into account part of constraints, it is believed that the effectiveness of ES-sto and its ability to achieve more benefit for systems with a higher degree of uncertainty are valid for general cases, although the exact quantity may be problem dependent. 7.6. Probability distributions of processing times A possible shortcoming is the assumption of known distributions of processing times. In maketo-order or engineer-to-order systems this can be difficult because of the inherent uniqueness of product, assembly and components. However, in many such industries, there are significant similarities between the products and operations. Historical data can therefore be used to provide process time distributions. On the other hand, this assumption can be relaxed to some extent. For example, instead of using a particular distribution, more

general distribution family can be used to represent stochastic operation times in the stochastic discrete event simulation (Song et al., 2002). 8. Conclusions This paper presents an ES method to optimise the raw material release times in make-to-order systems producing complex products with stochastic processing times and multiple resource constraints. The performance measure, which is the expected total cost consisting of work in progress holding costs and product earliness and tardiness costs, is evaluated by a simulation model through averaging over multiple replications. Compared with the backward scheduling based on infinity capacity, or the same search method but ignoring the randomness and using the expected values, the proposed method achieves significantly lower costs in a range of scenarios (with different uncertainty degrees) at a cost of much more computational time. It is found that more benefit can be achieved if the system has a higher degree of uncertainty. An illustrative example covering a range of product structures and a case study using industrial data demonstrate the results. The results confirm the intuition that it is beneficial to take into account the uncertainty in advance when designing the raw material release times. The experimentation and the case study also reveal that the effects of four priority rules on the cost are not very significant for the systems under consideration. Further research includes speeding up the search process and identifying the effects of uncertainties at different positions in the product structures on performance measures. An important but difficult issue is to establish tight lower bounds of performance measures, which can be used to judge the quality of a solution. In addition, extending the ES method to solve the re-planning and rescheduling problems in make-to-order systems is interesting. Acknowledgments The author would like to thank referees’ for careful reading and helpful comments, which significantly improved the presentation of the paper. References Aytug, H., Khouja, M., Vergara, F.E., 2003. Use of genetic algorithms to solve production and operations management

ARTICLE IN PRESS D.-P. Song / Int. J. Production Economics 103 (2006) 371–385 problems: A review. International Journal of Production Research 41 (17), 3955–4009. Back, T., 1996. Evolutionary Algorithms in Theory and Practice: Evolution Strategies, Evolutionary Programming, Genetic Algorithms. Oxford University Press, Oxford. Buss, A.H., Rosenblatt, M.J., 1997. Activity delay in stochastic project networks. Operations Research 45 (1), 126–139. Elmaghraby, S.E., 2001. On the optimal release time of jobs with random processing times with extensions to other criteria. International Journal of Production Economics 74, 103–113. Elmaghraby, S.E., Ferreira, A.A., Tavares, L.V., 2000. Optimal start times under stochastic activity durations. International Journal of Production Economics 64, 153–164. Gong, L.G., Kok, T., Ding, J., 1994. Optimal lead-times planning in a serial production system. Management Science 40 (5), 629–632. Hasan, C.N., Spearman, M.L., 1999. Optimal material release times in stochastic production environments. International Journal of Production Research 37 (6), 1201–1216. Hicks, C., 2004. A genetic algorithm tool for designing manufacturing facilities in the capital goods industry. International Journal of Production Economics 90, 199–211. Hodgson, T.J., Russell, E.K., Stanfield, P.M., 1997. Ready-time scheduling with stochastic service times. Operations Research 45 (5), 779–783.

385

HomemdeMello, T., Shapiro, A., Spearman, M.L., 1999. Finding optimal material release times using simulation-based optimisation. Management Science 45 (1), 86–102. Lagodimos, A.G., Mihiotis, A.N., Kosmidis, V.C., 2004. Scheduling a multi-stage fabrication shop for efficient subsequent assembly operations. International Journal of Production Economics 90, 345–359. Liao, C.J., 1992. Optimal release time on a stochastic single-machine. International Journal of Systems Science 23 (10), 1693–1701. Pongcharoen, P., Hicks, C., Braiden, P.M., Stewardson, D.J., 2002. Determining optimum genetic algorithm parameters for scheduling the manufacturing and assembly of complex products. International Journal of Production Economics 78, 311–322. Song, D.P., 2001. Stochastic models in planning complex engineer-to-order products. Ph.D. Thesis, University of Newcastle upon Tyne, UK. Song, D.P., Hicks, C., Earl, C.F., 2001. Setting planned job release times in stochastic assembly systems with resource constraints. International Journal of Production Research 39 (6), 1289–1301. Song, D.P., Hicks, C., Earl, C.F., 2002. Product due date assignment for complex assemblies. International Journal of Production Economics 76 (3), 243–256. Yano, C.A., 1987. Setting planned lead-times in serial production systems with tardiness costs. Management Science 33 (1), 95–106.