Group shops scheduling with makespan criterion subject to random release dates and processing times

Computers & Operations Research 37 (2010) 152 -- 162

Contents lists available at ScienceDirect

Computers & Operations Research journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / c o r

Group shops scheduling with makespan criterion subject to random release dates and processing times Fardin Ahmadizar a,∗ , Mehdi Ghazanfari b , Seyyed Mohammad Taghi Fatemi Ghomi c a

Department of Industrial Engineering, University of Kurdistan, Pasdaran Boulvard, Sanandaj, Iran Department of Industrial Engineering, Iran University of Science and Technology, Narmak, Tehran 16846 13114, Iran c Department of Industrial Engineering, Amirkabir University of Technology, 424 Hafez Avenue, Tehran, Iran b

A R T I C L E

I N F O

Available online 18 April 2009 Keywords: Group shop scheduling problem Random release dates Random processing times Makespan Ant colony algorithm Simulation

A B S T R A C T

This paper deals with a stochastic group shop scheduling problem. The group shop scheduling problem is a general formulation that includes the other shop scheduling problems such as the flow shop, the job shop and the open shop scheduling problems. Both the release date of each job and the processing time of each job on each machine are random variables with known distributions. The objective is to find a job schedule which minimizes the expected makespan. First, the problem is formulated in a form of stochastic programming and then a lower bound on the expected makespan is proposed which may be used as a measure for evaluating the performance of a solution without simulating. To solve the stochastic problem efficiently, a simulation optimization approach is developed that is a hybrid of an ant colony optimization algorithm and a heuristic algorithm to generate good solutions and a discrete event simulation model to evaluate the expected makespan. The proposed approach is tested on instances where the random variables are normally, exponentially or uniformly distributed and gives promising results. © 2009 Published by Elsevier Ltd.

1. Introduction Manufacturing systems are subject to uncertain events due to both human and machine resource factors. This uncertainty may be caused by machine failure, worker unavailability, failure to take delivery of raw material, change in availability dates, change in the customer orders, etc. Consequently, considering a scheduling problem in a stochastic situation is more realistic than in a deterministic situation. A great deal of research has been conducted on scheduling problems. Many of these assume that parameters such as processing times are fixed and deterministic. In fact, deterministic scheduling problems relying on precise data have received considerable attention and a lot of papers exist concerning these problems. However, the uncertainty in parameters has not received enough attention and work remains to be done concerning stochastic scheduling problems in particular job shops and open shops, in which probability distributions are used for characterizing uncertainty. Usually, in stochastic scheduling models the number of jobs is known in advance. However, the processing times are assumed to be random variables from known probability distributions. Different

∗ Corresponding author. Tel./fax: 98 871 6660073. E-mail addresses: [email protected], [email protected] (F. Ahmadizar), [email protected] (M. Ghazanfari), [email protected] (S.M.T. Fatemi Ghomi). 0305-0548/$ - see front matter © 2009 Published by Elsevier Ltd. doi:10.1016/j.cor.2009.04.002

operations may have dissimilar processing time distributions. The release dates and due dates of jobs may also be random variables from known distributions. The actual realization of a random variable only becomes known when it actually occurs. Since such a problem is a stochastic optimization problem, the decision maker has to decide a strategy that optimizes the objective in some stochastic sense, most often in expectation. In stochastic scheduling, such forms of randomness can be formulated in several ways; for example, random breakdowns may be modeled as either an integral part of the processing times or a separate stochastic process [1]. In the literature, there are a few references where stochastic scheduling problems are considered. The most commonly considered problems in this field are those that assume machines to be completely reliable and processing times to be stochastic. Some recent works about stochastic job shop and open shop scheduling problems are first given. Singer [2] has studied a job shop with random processing times and individual due dates, where the expected total weighted tardiness must be minimized. He has proposed a heuristic that amplifies the expected processing times by a factor and measures the quality of a given policy using a risk averse penalty function combining the expected difference and the worst case difference from the optimal schedule. Some researchers have considered job shop scheduling problems with random processing times under various costs for not completing a job just on time, for minimizing the average value of total costs. The general procedure has been that heuristic decision-making rules are used in situations when several jobs

F. Ahmadizar et al. / Computers & Operations Research 37 (2010) 152 -- 162

are ready to be served on one machine [3–5]. Yoshitomi [6] has proposed a simulation-based genetic algorithm for solving a job shop with random processing times, in order to minimize the expected makespan. Yoshitomi and Yamaguchi [7] have enhanced this genetic algorithm by applying a new crossover operator. In these two algorithms, solutions having very high frequency through all generations are selected as good solutions. Tavakkoli-Moghaddam et al. [8] have considered a job shop with random operations, where the difference between the delivery and the completion times of jobs as well as related operational or idle cost of machines must be minimized. They have presented a hybrid method consisting of a neural network to generate an initial solution and a simulated annealing algorithm to improve the quality of the initial solution. Some authors have developed exact and heuristic algorithms for job shop scheduling problems with the makespan or mean flow time criterion subject to random processing times, which can take any real value between known lower and upper bounds [9–11]. Luh et al. [12] have presented an approach based on a combined Lagrangian relaxation and stochastic dynamic programming for a job shop considering uncertain arrival times, processing times, due dates, and part priorities, for the objective of minimizing expected job tardiness and earliness cost. Neumann and Schneider [13] have proposed heuristics based on aggregate schedules for a job shop scheduling problem with stochastic precedence constraints, where the expected makespan is to be minimized. Finally, Alcaide et al. [14] have studied open shop scheduling problems with random processing times, breakdowns and repair times, where the expected makespan must be minimized. They have developed a general dynamic procedure that converts any problem into a set of scheduling problems without breakdowns. This paper deals with stochastic group shop scheduling (SGSS) problems. The group shop scheduling (GSS) problem was first introduced in the context of a mathematical contest in 1997 [15]. The GSS problem frequently occurs in industrial environments and includes the job shop and the open shop scheduling problems. Consequently, the GSS problem is NP-hard and difficult to solve optimally. However, very few papers have concentrated on this issue [16,17]. The present paper is focused on the SGSS problem subject to stochastically distributed release dates and processing times, governed by known distribution functions. The objective is to find a schedule in the class of static policies to minimize the expected completion times of all jobs, or makespan. The problem is first formulated in a form of stochastic programming and then an approach is developed which is a combination of an ant colony optimization (ACO) algorithm and a heuristic algorithm to generate good solutions and a discrete event simulation model to evaluate the objective. The proposed approach is tested on large instances. The rest of the paper is organized as follows. In the next section, the SGSS problem is introduced and formulated. Section 3 provides a lower bound on the expected makespan, which may be used as a measure for evaluating the performance of a solution without simulating. The proposed approach is described in Section 4. Computational results and some concluding remarks are provided in Section 5. Finally, Section 6 gives a summary as well as future work. 2. Stochastic group shop scheduling problem In the following, a description of a general deterministic shop scheduling problem, called the GSS problem, including the job shop and the open shop scheduling problems is first provided and then the stochastic case is presented. The GSS problem can be formulated, in general, as a problem of processing a set of N jobs on a set of M machines. Job j which is a set of Oj operations has at most one operation to be processed on each machine. No job may be processed on more than one machine at a time and no machine can process more than one job at a time.

153

Moreover, preemption is not allowed. The machines are available and ready at time zero. The release date of job j is Rj , that is, the processing of this job cannot be started before that. Job j is processed during a time Pij by machine i. Setup times are sequence independent and included in processing times, and transportation times are negligible. Between two machines, jobs can wait in an unlimited buffer. In addition, the operations of each job j are partitioned into Gj groups, where the operations in the same group are unordered while operations in two distinct groups satisfy the predefined precedence constraint between the groups. If each group consists of only one operation, that is, for all j Gj = Oj , this GSS problem is equivalent to a job shop scheduling problem. If for each job, the number of groups is equal to 1, that is, for all j Gj = 1, the GSS problem is equivalent to an open shop scheduling problem. In the GSS problem, processing orders for all machines and all groups have to be determined in order to minimize a given optimality criterion such as makespan. In the SGSS problem, some classical assumptions for the deterministic case may be replaced by the following: processing times as well as release dates are not known in advance. The SGSS problem is concerned with finding a job schedule which optimizes the criterion in some stochastic sense, most often in expectation. To solve scheduling problems, two approaches are proposed, namely the static and the dynamic scheduling. In the static scheduling approach, the schedule is determined before the start of processing, while in the dynamic context it is constructed at each time that an event such as end of processing by a machine and arrival of a new job occurs. In this paper, the group shop static scheduling problem is considered with stochastically and independently distributed release dates and processing times. The release date Rj follows a distribution function given by FRj (t) with mean E(Rj ) and variance V(Rj ); assume that each job arrives at the shop independently and this differs from the situation that all jobs arrive according to a joint distribution considered in the queuing theory. Moreover, the processing time Pij follows a distribution function given by FPij (t) with mean E(Pij ) and variance V(Pij ). The objective is to find a job schedule which minimizes the expected makespan, E(Cmax ), where Cmax is the makespan, i.e., the maximum completion time. It is clear that if all jobs release after time zero, the makespan contains the period of time at the beginning of the production schedule in which all machines are kept idle. 2.1. Stochastic disjunctive graph representation The SGSS problem can be represented by a stochastic disjunctive graph G(ND,A) adapted from the one proposed by Roy and Sussman [18], where ND is the set of nodes and A is the set of arcs. The set of nodes ND contains one element for each operation, a source node and a sink one. The source is equivalent to a dummy operation Se preceding all other nodes and the sink to a dummy operation En succeeding all other nodes. Both dummy operations have zero processing time. The set of arcs A includes both conjunctive (directed) and disjunctive (undirected) arcs. The set of conjunctive arcs contains three types of arcs: the set of arcs CAO reflecting the predefined precedence constraints between successive groups of a job, the set of arcs CAR—corresponding to the dummy operation Se—reflecting the job release dates constraints and the set of arcs CAE corresponding to the dummy operation En. Furthermore, the set of disjunctive arcs contains two types of arcs: the set of arcs DAG which connect mutually unordered operations belonging to the same group and the set of arcs DAM which connect mutually unordered operations requiring the same machine for their execution. Each of the arcs has a random length. Each arc of CAO and CAE has a length equal to the processing time of the operation where the arc begins and each arc of CAR has a length equal to the release date of its corresponding

154

F. Ahmadizar et al. / Computers & Operations Research 37 (2010) 152 -- 162

arcs of CAO

arcs of DAG

arcs of CAR

arcs of DAM

arcs of CAE

P31

P21 P21 J1 m2

P31

J1 m3

J1 m1

J1 m4

R1

P41 R1

P11 R2

Se

J2 11111111P12 m1

J2 m3

P32

R3

En P33

P43

J3 m4

P43

J3 m1

P23

J3 m2

J3 m3

P13 Fig. 1. The stochastic disjunctive graph for the given instance

P31

J1 m3 R1

Se

P31 R2

J2 m1

R3

P21

J1 m2

J1 m1

P41

P12

J2 m3

P12 P43

J1 m4

P21 P32

P32

P43 J3 m4

P11

J3 m1

P13 P13

J3 m2

P23

En P33

J3 m3

Fig. 2. A feasible solution of the problem instance shown in Fig. 1.

job, while the length of each disjunctive arc depends on the decision. After a disjunctive arc is directed, it acts as an arc of CAO. In the disjunctive graph model, the scheduling problem is equivalent to selecting one arc in each disjunction. By fixing directions of all disjunctive arcs, the order of all operations belonging to the same group as well as the order of all operations requiring the same machine is determined and a complete solution is obtained. In addition, the resulting directed graph has to be acyclic. The expected length of the longest path from the source to the sink determines the expected makespan. The corresponding disjunctive graph for a simple instance of the SGSS problem is illustrated in Fig. 1. This group shop consists of three jobs (J1 , J2 and J3 ) and four machines (m1 , m2 , m3 and m4 ), where job J1 has four operations partitioned into two groups, job J2 has two operations partitioned into two groups and job J3 has four operations partitioned into three groups. Job J1 may be processed first on machine m2 and then on machine m3 or vice versa, i.e., the first two operations of this job belong to the same group, and subsequently the job may be processed first on machine m1 and then on machine m4 or vice versa, and so on. Fig. 2 illustrates a directed graph that corresponds to a feasible solution to the SGSS problem instance shown in Fig. 1 (these figures only show necessary arcs). Given a feasible solution, a new feasible schedule can be obtained by reversing some disjunctive arcs in such a way that the obtained graph is also acyclic.

2.2. Stochastic programming The SGSS problem is formulated as a stochastic programming problem based on its stochastic disjunctive graph representation. To present this formulation, let the decision variable Tij be the earliest possible starting time of the processing of job j on machine i which is denoted by operation (i,j). In this formulation, the makespan is equal to the earliest possible starting time of the dummy operation En. Moreover, the starting time of the dummy operation Se is equal to zero. To ensure that each operation belonging to a group of a job cannot start before all operations belonging to preceding groups of the same job are completed, a set of constraints is necessary as follows: Tij − Ti j  Pi j

for all (i , j) → (i, j) ∈ CAO.

(1)

With regard to the set of arcs DAG, to ensure that some ordering exists among operations belonging to the same group of a job, a set of disjunctive constraints is required as follows: Tij − Ti j  Pi j

or Ti j − Tij  Pij

for all (i , j) ↔ (i, j) ∈ DAG.

(2)

Furthermore, with regard to the set of arcs DAM, to ensure that some ordering exists among operations of different jobs requiring the same machine for their execution, a set of disjunctive constraints

F. Ahmadizar et al. / Computers & Operations Research 37 (2010) 152 -- 162

is required as follows: Tij − Tij  Pij

or Tij − Tij  Pij

for all (i, j ) ↔ (i, j) ∈ DAM.

Since the max function is convex [20], applying Jensen's inequality to the right side of the above equation gives (3)

To ensure that the processing of each job cannot start before its release date, a set of constraints is necessary as follows: Tij  Rj

for all (Se) → (i, j) ∈ CAR.

(4)

Finally, the following set of constraints is required to ensure that the makespan is equal to the completion time of all jobs: Cmax − Tij  Pij

for all (i, j) → (En) ∈ CAE.

(5)

The proposed model for the SGSS problem may then be expressed as: Model 1: Minimize the objective function subject to (1)–(5). In this stochastic programming problem, because of the random changeability of the release dates and processing times, each schedule can have a distribution of objective function values. Since the makespan is a random variable and solution methods for ordinary mathematical programming problems cannot be directly applied, a simulation-based optimization approach is proposed to minimize the expected makespan in Model 1. It is noteworthy that a feasible solution to the problem given in Model 1 can only be found by turning each of the disjunctive constraints in (2) and (3) into a conjunctive one, i.e., by selecting one constraint from each pair and deleting the other one, such that the selected constraints, together with the conjunctive constraints in (1), (4) and (5), do not conflict.

In the stochastic shop scheduling problem under consideration, a lower bound can be obtained for the expected makespan objective by replacing all the random release dates and processing times by their expected values and turning the stochastic group shop into a deterministic one. Let S denote the set of all feasible static, i.e., predetermined, schedules for the SGSS problem and assume that s s be the stochastic directed is a member of S. Let Gs (ND, A) and Cmax graph and the makespan corresponding to solution s, respectively. Recall that the resulting graph Gs (ND, A) has no cycles and all arcs in A are conjunctive. Define the deterministic directed graph Gˆ s (ND, A) in the same way as Gs (ND, A), but replacing each random variable by s denote the length of the longest path its expected value. Let Cˆ max from node Se to node En in graph Gˆ s (ND, A), which can be computed using the algorithm proposed by Bellman [19]. Proposition 1. In the SGSS problem with independent random release dates and processing times, the expected makespan satisfies the inequality given in (6): for all s ∈ S.

(6)

Proof. Assume that Cjs denotes the completion time of job j in graph Gs (ND,A) corresponding to a schedule s. Then, the expected makespan is  s )=E E(Cmax

 Max Cjs . j=1,2, ...,N

Let Tijs denote the earliest possible starting time of operation (i,j) in graph Gs (ND, A). Suppose that under schedule s the last operation of job j is processed by machine l. Clearly,  s )=E E(Cmax

 Max (Tljs + Plj ) . j=1,2, ...,N

s ) E(Cmax

Max [E(Tljs + Plj )].

j=1,2, ...,N

Let Tˆ ijs denote the earliest possible starting time of operation (i,j) in graph Gˆ s (ND, A). Now, let us suppose that E(T s )  Tˆ s . Then, lj

s ) E(Cmax

Max

j=1,2, ...,N

[E(Tljs ) + E(Plj )] 

Max

j=1,2, ...,N

[Tˆ ljs

lj

+ E(Plj )],

and assuming that Cˆ js denotes the completion time of job j in graph Gˆ s (ND, A) results s ) E(Cmax

s Max Cˆ js = Cˆ max .

j=1,2, ...,N

At this point, we show that inequality E(Tijs )  Tˆ ijs holds for not only the last operation of each job j but also each operation (i,j), and this completes the proof of the proposition. In view of the fact that under each schedule s each operation (i,j) immediately follows at most two other operations, all possible cases are considered. Case 1: If operation (i,j) is immediately preceded by the dummy operation Se only, then E(Tijs ) = E(Rj ) = Tˆ ijs . Case 2: If operation (i,j) is immediately preceded by operation (i , j) only, which is an operation specified in case 1, then E(Tijs ) = E(Tis j + Pi j ) = Tˆ is j + E(Pi j ) = Tˆ ijs .

3. Lower bound on the expected makespan

s s )  Cˆ max E(Cmax

155

Case 3: If operation (i,j) is immediately preceded by the dummy operation Se and operation (i, j ), which is an operation specified in case 1, then E(Tijs ) = E[Max(Rj , Tijs  + Pij )]  Max[E(Rj ), Tˆ ijs  + E(Pij )] = Tˆ ijs . Case 4: If operation (i,j) is immediately preceded by the dummy operation Se and operation (i, j ), which is an operation specified in case 2 or case 3, the inequality is demonstrable in the same way as case 3. Case 5: If operation (i,j) is immediately preceded by operations (i , j) and (i, j ), which are operations specified in case 1, then E(Tijs ) = E[Max(Tis j + Pi j , Tijs  + Pij )]

 Max[Tˆ is j + E(Pi j ), Tˆ ijs  + E(Pij )] = Tˆ ijs . Case 6: If operation (i,j) is immediately preceded by operation (i , j), which is an operation specified in case 1 and operation (i, j ), which is an operation specified in case 2, the inequality is demonstrable in the same way as case 5. It is clear that for each of the other possible cases, the inequality is provable in the same way as the previous cases.  In such a manner described in the next section, the above proposition can be used for making a judgment on the performance of each feasible solution to the SGSS problem. 4. Description of the proposed approach In a stochastic optimization context, if the expected objective value can be obtained analytically for each feasible solution, then the problem is considered to be a deterministic optimization problem which can be solved either analytically or numerically by mathematical programming problems solution methods [21]. On the other hand, the performance analysis in complex problems, where the objective

156

F. Ahmadizar et al. / Computers & Operations Research 37 (2010) 152 -- 162

function cannot be analytically expressed in terms of the model variables, must be considered via Monte Carlo simulation [22]. Consequently, in order to solve the SGSS problem, a simulation optimization approach is developed that is an efficient combination between an optimization algorithm to generate high-quality solutions and a discrete event stochastic simulation model to evaluate the expected makespan. Furthermore, the proposed optimization procedure is a hybrid of an ant colony algorithm and a heuristic algorithm. Recently, ACO algorithms which belong to the class of constructive metaheuristics have been used for solving the combinatorial optimization problems. The first research effort was made by Dorigo in 1992 [23], and an introduction to the ACO algorithms was then provided by Dorigo et al. [24]. ACO algorithms are nature-inspired, stochastic search procedures derived from the behavior of real ants which use pheromones as a communication medium and are capable of finding the shortest path from their nest to a food source without using visual cues. This characteristic of real ants is exploited in a colony of simple agents, called artificial ants, in order to solve the combinatorial optimization problems. In an ACO algorithm, each artificial ant probabilistically constructs a solution by starting with a null one and iteratively adding solution components until a complete solution is built. The construction of solutions by artificial ants is guided by taking into account the pheromone trails and the problemspecific heuristic information. At run-time, to reflect the agents' acquired search experience, the pheromone trails change dynamically according to a function of the values of the objective function for the previously generated solutions [25,26]. Some methods such as simulated annealing and genetic algorithm requiring only the evaluation of the objective function for complete solutions are often used for solving stochastic optimization problems. Unlike these methods, ACO algorithms require the performance evaluation of not only complete solutions for modifying the pheromone trails and updating the best constructed solution, but also partial solutions for assigning the heuristic information to solution components. Consequently, ACO algorithms have to be adapted significantly for a stochastic optimization context. In the remainder of this section, the general structure of the proposed approach is first provided and then the details are presented. 4.1. General structure of the approach In the developed approach, an ACO algorithm is applied to generate high-quality solutions to the SGSS problem. A heuristic algorithm is used to within the ACO framework to construct feasible solutions. Considering the lower bound on the expected makespan, a constructed solution that may improve the best solution constructed so far is specified; we refer to such a solution by a satisfactory solution. The performance of a satisfactory solution is evaluated by means of a discrete event simulation model. In other words, the performance of a constructed solution that is not satisfactory is evaluated without being simulated. A satisfactory solution is then compared with the current best one by applying a statistical test. In addition, the pheromone trails are dynamically updated and limited. The general structure of the approach can be represented as follows: Algorithm 1. Simulation-based ACO algorithm. Step 1. Set the parameters and initialize the pheromone trails. Step 2. While the termination condition is not met, do: Step 2.1. For each ant in the colony do: Step 2.1.1. Construct a complete solution using Algorithm 2; while constructing the solution, also update each pheromone trail that is compatible to the sequence by applying the local updating rule.

Step 2.1.2. If the solution is satisfactory, evaluate its performance using Algorithm 3, i.e., the simulation model; then, compare the solution with the best one constructed so far by applying the statistical test; and in case of an improved solution, update the best solution found and limit the pheromone trails. Step 2.2. Update the pheromone trails by applying the global updating rule; while updating the pheromone trails, also limit them. Step 2.3. If the restart condition is met, reset the pheromone trails. Step 3. Return the best solution found. 4.2. Pheromone model To apply an ACO algorithm to a combinatorial optimization problem, a pheromone model which forms a kind of adaptive memory of previously found solutions should properly be defined. To solve the SGSS problem, a pheromone model is used that differs from the relation-learning model developed for solving the GSS problem [16,27,28]. In the relation-learning model, pheromone trails are associated with pairs of related operations, where two operations are called related if they either belong to the same group or have to be processed by the same machine, while in the proposed model a pheromone trail is associated with the assignment of an operation to a position in its related group or machine. (1)g Let ij be the trail intensity of setting operation (i,j) in the po(2)m

sition g of its related group and ij denote the trail intensity of setting this operation in the position m of machine i. Each of them represents the desirability of processing (scheduling) operation (i,j) in the given position. In other words, a pheromone matrix is associated with each machine as well as each group given that a decision is required, i.e., when the number of its corresponding operations is more than one. Like most application of ACO algorithms, a fixed value is first assigned to all pheromone trails. Then, at run-time the trail intensities are modified by applying both the local and global updating rules, and to keep the algorithm from converging towards a solution, like the max–min ant system [29], the trail intensities are always limited between a lower bound min and an upper bound max . In addition, the maximum and minimum trail bounds are dynamically modified whenever the best solution obtained so far is updated. 4.3. Construction of a schedule In the developed approach, each artificial ant probabilistically builds a feasible solution according to a list scheduler algorithm [30], a heuristic procedure to solve shop scheduling problems, by starting with the null sequence and iteratively adding operations until a complete solution which is a permutation of all operations is constructed. At each construction iteration, an unscheduled operation can only be chosen if all of its predecessors are already in the partial solution; in that case the operation is said to be schedulable. Therefore, an ant constructs a feasible permutation of all operations step by step until all the groups and machines are scheduled. It is noteworthy that each feasible permutation of all operations is equivalent to one feasible solution—under each realization of the release dates and processing times—in the disjunctive graph model as well as Model 1 (there is generally a many-to-one mapping). At each iteration of a list scheduler algorithm an operation is often chosen by means of dispatching rules, while in the proposed algorithm this is done by applying a transition rule, the same as the pseudo-random

F. Ahmadizar et al. / Computers & Operations Research 37 (2010) 152 -- 162

proportional rule employed in the ant colony system [31], taking into account the pheromone trails and the heuristic information. Let Ps be a partial schedule and define the deterministic directed graph Gˆ Ps (ND, A) corresponding to Ps, where all arcs in A are conjunctive and each random variable is replaced by its expected value. At the start, in which Ps is empty, A consists of only CAO, CAR and CAE, and after appending at each construction iteration an operation from a set of allowed operations to Ps, the corresponding disjunctive arcs are directed and added to A. Let Tˆ ijPs denote the earliest possible starting time of schedulable operation (i,j) in graph Gˆ Ps (ND, A), which is equal to the length of the longest path from node Se to node (i,j). At the start of the algorithm, the set of schedulable operations consists of all operations (i,j) belonging to the first group of each job and for each of them it holds that Tˆ ijPs = E(Rj ). The heuristic algorithm applied to construct solutions to the SGSS problem can then be detailed as follows: Algorithm 2. Solution construction mechanism. Step 1. Define SO as the set of schedulable operations. Start by setting Ps as the null partial schedule. Step 2. If there is an operation (i,j) in SO with no related and unscheduled operations left, i.e., an operation that is the only unscheduled one belonging to its related group and also the only unscheduled operation among all the unscheduled ones requiring machine i, choose it (when two or more operations are candidates to be chosen, select one of them randomly) and go to Step 6. Otherwise, go to the next step. Step 3. Determine the earliest completion time ECT among all operations in SO as ECT = Min(i,j)∈SO [Tˆ ijPs + E(Pij )]. Determine one of the operations (i∗ , j∗ ) with minimal completion time ECT. Define RSO as the set of all operations (i,j) in SO such that at least one of the conditions i = i* , j = j* (i.e., belong to the same group as operation (i∗ , j∗ )) or Tˆ ijPs < ECT is satisfied. Step 4. Generate a random number q uniformly distributed in [0, 1]. If q  q0 , where q0 is a parameter between 0 and 1 determining the relative importance of exploitation versus exploration, select the unscheduled operation (i,j) in RSO as (7), and go to Step 6. Otherwise, go to the next step. (i, j) = arg max[(ij ) (ij ) ],

(7)

where the trail intensity ij represents the `desire' of placing operation (i,j) in the current position of the partial schedule, ij denotes the heuristic information assigned to this operation, and  and  are two positive parameters determining the relative importance of the pheromone trail and the heuristic information. Step 5. Select one of the unscheduled operations in RSO according to the following transition probabilities (exploration): pij = 

(ij ) (ij )

(i ,j )∈RSO (i j )



(i j )

:

if (i, j) ∈ RSO.

(8)

Step 6. Add the selected operation to Ps as early as possible and update SO. Return to Step 2 until a complete solution is constructed. The trail intensity ij employed in Steps 4 and 5 of the above algorithm is calculated as follows. If operation (i,j) is the only unscheduled one belonging to its related group and not the only unscheduled operation among all the unscheduled ones requiring machine i, then ij = (2)m , where m−1 operations requiring the same machine are ij already in the partial solution. If operation (i,j) is the only unscheduled operation among all the unscheduled ones requiring machine

157

i and not the only unscheduled one belonging to its related group, (1)g then ij = ij , where g−1 operations belonging to the same group are already in the partial solution. If operation (i,j) is not the only unscheduled one belonging to its related group and also not the only unscheduled operation among all the unscheduled ones requiring (1)g (2)m machine i, then ij = Min(ij , ij ), where g−1 operations belonging to the same group and also m−1 operations requiring the same machine are already scheduled. Note that the other case is considered by Step 2. Moreover, the heuristic information which is often needed for achieving a high algorithm performance [26] is assigned to each operation in RSO according to the earliest starting time dispatching rule, depending on the partial schedule, as follows:

ij =

1 1 + Tˆ ijPs

if (i, j) ∈ RSO.

:

(9)

At each iteration of Algorithm 2, taking the expected values into account, the set of schedulable operations is restricted and then an operation from the restricted set is selected by applying a rule exploiting both the expected values (the heuristic information) and the simulation results (the pheromone trails, as shown in Section 4.5). From Proposition 1, it follows that for all operations Tˆ ijPs  E(TijPs ), where, with regard to the random variables, TijPs denotes the earliest possible starting time of operation (i,j) depending on partial schedule Ps. As a result, each of the expected values may be multiplied by a factor in order to increase each Tˆ ijPs towards E(TijPs ). Since it is not practicable to assign an accurate factor to each random variable, one may use the same expansion factor for each of them. Clearly, in such a situation, the resulted algorithm is similar to Algorithm 2, except for an insignificant difference in the heuristic information calculation. 4.4. Performance evaluation of a constructed solution As previously mentioned, we propose to use a discrete event simulation model to evaluate the performance measure of a constructed solution to the SGSS problem. Since the input variables to the stochastic simulation model, i.e., the release dates and processing times, are random, the performance measure, i.e., the makespan, is also a random variable, which may be evaluated through its mean. Each execution of the simulation model does not directly yield the expected makespan of the solution but it gives only an estimation of the expected makespan. Consequently, to evaluate the estimation accuracy, we need to calculate a confidence interval on the expected criterion. By executing the simulation model, some statistically independent replications are successively realized. At each replication, a sample of the release dates and processing times according to the distribution functions is first chosen and the makespan of the solution is calculated. The mean of the obtained results for all replications is then an estimation of the expected makespan. Suppose that s X Cmax denotes the output of the simulation model corresponding to solution s. Proposition 2. In the SGSS problem under consideration, the output of the simulation model should satisfy the inequality stated in (10): s s X Cmax  Cˆ max

for all s ∈ S.

(10)

s

s ) and from Proposition Proof. Since X Cmax is an estimation of E(Cmax s s , it is clear 1 we know that E(Cmax ) is larger than or equal to Cˆ max that the inequality should be satisfied. 

In other words, if for a given number of replications we s s , the simulation has to be continued until the get X Cmax < Cˆ max

158

F. Ahmadizar et al. / Computers & Operations Research 37 (2010) 152 -- 162

inequality stated in (10) is satisfied. This happens due to the fact that, from the strong law of large numbers, as the number of observations increases, the estimation accuracy increases too, that is, s s X Cmax converges towards E(Cmax ). In stochastic simulation, a statistical test has to be conducted in order to compare solutions properly [32]. Let bs denote the best solution constructed up to now. A two-sample comparison is then performed to determine if there is a statistically significant difference s

bs

between the two estimates X Cmax and X Cmax , i.e., there is an improvement. Assuming that the number of observations for each of them is large enough, the following statistic is used: s

bs

X C − X Cmax , Z =  max Sd2bs Sd2s + ns nbs

(11)

where Sds and ns (Sdbs and nbs ) are, respectively, the estimated standard deviation and the number of all replications for solution s (bs). s

bs

Statistically, the two estimates X Cmax and X Cmax are not significantly different when Z ∈ [−Z/2 , Z/2 ], where 1− is the level of confidence and Z/2 (−Z/2 ) denotes the point that under the standard normal curve, its right (left) area is equal to /2; in this case, we keep the current best solution bs, because it cannot be considered that there is an improvement. By contrast, the two estimates are significantly different and hence can be compared when Z ∈/ [−Z/2 , Z/2 ]; in the case of Z > Z/2 , bs is significantly better than s and consequently solution s has to be rejected, i.e., the current best solution bs is kept, and in the case of Z < − Z/2 , s is significantly better than bs and consequently the current best solution has to be replaced by s. As seen, if it holds that

s X Cmax

s

bs  X Cmax , solution s is rejected, but it is possible that bs than X Cmax , nevertheless solution s is again rejected.

X Cmax is smaller Considering the above explanations, the following proposition can efficiently be used for specifying, without simulating, a solution that cannot improve the current best one, i.e., a solution that is not satisfactory:

Algorithm 3. Simulation model. Step 1. Set t = 1 and Nr as the number of replications at each iteration, that is, the minimum number of replications. Step 2. For rp = (t − 1)Nr + 1, . . . , tNr do: Step 2.1. Choose a sample of the release dates and processing times randomly, according to the distribution functions. s . Step 2.2. Determine the makespan of solution s, Xrp Step 3. Calculate the mean of the obtained results for all replications,  s s s ˆs that is, X Cmax = tNr rp=1 Xrp /(tNr). If it holds that X Cmax  Cmax , according to Proposition 2, stop. Otherwise, set t = t+1 and return to Step 2. Determining the number of replications is an important issue because a very large number may waste time, while a very small number may lead to an inaccurate estimation. The required number of replications can be determined in order to achieve a desired accuracy which is often expressed as the margin of error, that is, the plus or minus part of the confidence interval. In this study, to assure a good sampling, it has been desired to have the margin of error of the (1−)×100% confidence interval to be mostly < 3% of the output of the simulation model. The experimental analysis has shown the minimum number of replications about 200 to perform well, thus we set Nr = 200. 4.5. Updating of the pheromone trails After each solution construction iteration, the local updating rule is applied to each pheromone trail corresponding to the selected operation. For example, if at an iteration operation (i,j) has been added to the partial solution and this selection is compatible to set this operation in the position g of its related group and also to set it in the position m of machine i, the following pheromone updates are applied: (1)g (1)g = (1 −  )ij +  min , ij

Proposition 3. In the SGSS problem under consideration, if a solution s satisfies the inequality stated in (12), it cannot be better than the solution bs and therefore, it is not necessary to evaluate the performance measure of solution s using the simulation model: bs s  X Cmax . Cˆ max

(12)

Proof. Assume that we evaluate the performance measure of solus tion s using the simulation model, and let X Cmax be the output. From s

bs

Proposition 2 and (12) it is clear that X Cmax  X Cmax and consequently, according to the above statistical test, solution s cannot be better than the best solution constructed so far. Hence it follows logically that this solution can be rejected without being simulated.  In fact according to the above proposition, it is determined as quickly as possible, without wasting time, that a constructed solution is satisfactory or not. If the solution is not satisfactory, its performance evaluation no longer continues and it is directly rejected. Otherwise, its performance has to be evaluated by the simulation model which is time-consuming to determine whether there is an improvement. It should be noted that while building the solution, whether or not it is satisfactory, each pheromone trail that is compatible to the sequence is locally updated. The algorithm of the stochastic simulation model proposed for evaluating the expected makespan of a satisfactory solution s, i.e., bs , can a constructed solution s which satisfies inequality Cˆ s < X max

then be detailed as follows:

Cmax

(2)m (2)m = (1 −  )ij +  min , ij

(13)

where  , a parameter between 0 and 1, is the local pheromone trail evaporation rate. Depending on the manner of selecting operation (i,j), it is possible that only one or none of above two updates is essential. For example, if operation (i,j) has been selected in Step 2 of Algorithm 2, since none of the pheromone trails has been used in selecting it, then both of them are unnecessary. The effect of this updating is to decrease the desirability of constructed solutions in order to explore different areas by the next ants in the colony, i.e., one ant's solution will be chosen with a lower probability by other ants. Once all ants in the colony have built their solutions, the global updating rule is applied to the pheromone trails. First, the amount of each pheromone trail is modified as follows: (1)g (1)g = (1 − )ij , ij

(2)m (2)m = (1 − )ij , ij

(14)

where , a parameter between 0 and 1, is the global pheromone trail evaporation rate. Then, the amount of each pheromone trail that is compatible to the best solution constructed so far is modified as follows: (1)g (1)g = ij + , ij

(2)m (2)m = ij + , ij

(15)

where

 =

H bs X Cmax

,

(16)

F. Ahmadizar et al. / Computers & Operations Research 37 (2010) 152 -- 162

where H is a positive parameter. Recall that for some operations (i,j), it is possible that only one or none of the two updates stated in (15) is essential. According to the global updating rule, all of the trail intensities first evaporate and then only those corresponding to the current best solution increase in order to make the search more directed. As previously mentioned, the trail intensities are always limited between min and max . These bounds are calculated as follows:

max =

 , 

min = 0.01max .

(17)

At the start of Algorithm 1 and also each pheromone trails resetting, the current max is assigned to all pheromone trails (at the start, we set max = 1). Since  is not constant, max as well as min is modified whenever the best solution obtained so far is updated. At that time, all pheromone trails are checked and if a trail intensity is smaller than the new min or greater than the new max , it is set to min or max , respectively. Regarding the old and new amounts of max , only the first constructed solution may cause a trail intensity to be greater than the new max . The local updating rule does not cause a trail intensity to be smaller than min , but since the first phase of the global updating rule may cause such a situation, after this phase and before the second one all pheromone trails are checked and if a trail intensity is smaller than min , it is set to min . Moreover, the second phase of the global updating rule does not cause a trail intensity to be greater than max . As seen, the pheromone trails are based on the simulation results. 5. Computational results and conclusions To evaluate the performance of the proposed simulation-based ACO approach, we use the only real GSS instance, called whizzkids97 [15], which is large enough to consider. It consists of 15 machines, 20 jobs that are all available at time zero, and the total number of 197 operations partitioned into 124 groups. To define a release date for each job, an integer has been generated from the uniform distribution (5,30) for each of them. The obtained release dates for the first job through the last one are, respectively, 11, 8, 9, 6, 11, 23, 18, 28, 22, 11, 22, 13, 17, 10, 8, 14, 14, 9, 12 and 17. For transforming this GSS instance to several SGSS ones, the release dates and processing times are assumed to be random variables which are normally, uniformly or exponentially distributed. The mean of each of the random variables is set to be equal to the corresponding value in the GSS

159

instance. Since the normal distribution as well as the uniform distribution has two parameters, the standard deviation is given as the ratio to the mean; the ratio of the standard deviation to the mean is referred to as the variability. For simplicity, the distribution as well as the variability is assumed to be the same for all the random variables. The ratios of the standard deviation to the mean are selected as 0.05, 0.1, 0.15, 0.2, 0.25 or 0.3. In addition, two different levels of confidence are considered: 0.99 and 0.95, i.e.,  = 0.01 and 0.05. The proposed algorithms have been coded in Visual C++ 6.0 under Microsoft Windows XP operating systems, running on a Pentium IV, 2 GHz PC with 256 MB memory. In the preliminary experiment, different values of the numeric parameters have been considered and to make the performance of the proposed approach robust, each combination of the parameter values has been tested on some of the problem instances and/or some ones generated by repartitioning the operations of each job into new groups. Then, the following values have been superior and used for all further experiments of this paper: 25 ants in the colony, q0 = 0.95,  = 0.1,  = 2,  = 0.01,  = 0.05 and H = 1. The pheromone trails are reset if the best constructed solution is not updated, i.e., does not improve, within 500 successive iterations in Step 2 of Algorithm 1. Furthermore, Algorithm 1 terminates when the total number of iterations in Step 2 reaches 3000. Due to the stochastic nature of the ACO search algorithm, each of the problem instances has been tested for ten trials. The computational results of the proposed simulation optimization approach applied to the problem instances are shown in Tables 1–5. For each problem instance and each trial, we give the expected makespan (point estimation), the corresponding confidence interval (interval estimation) and the lower bound on the expected makespan. The results concerning average expected makespan, CPU time (in seconds) and number of solutions evaluated without being simulated are also reported for each problem instance. Considering the results shown in Tables 1–4, it is seen that under the same conditions as the variability grows, the average expected makespan as well as the width of the corresponding confidence interval increases, which is reasonable. Moreover, as the variability grows, the average number of solutions evaluated via the simulation model becomes larger and since this model is time-consuming, the average CPU time increases too. With regard to Proposition 3, it follows that the number of solutions that have to be evaluated using the simulation model depends on the difference between the expected makespan and its lower bound. If this difference is larger, the probability that a solution is specified as a satisfactory one increases and therefore the number of solutions evaluated without

Table 1 Computational results for the normal distribution with  = 0.01. Trial

Variability 0.05

1 2 3 4 5 6 7 8 9 10 Average Time (s) NS

0.1 a

b

566.6 –563 (565.8,567.4) 553.9–551 (552.9,555.0) 560.0–555 (559.2,560.8) 568.9–564 (568.0,569.9) 563.4–561 (562.5,564.4) 545.6–538 (544.8,546.4) 558.9–554 (558.0,559.8) 554.1–549 (553.2,555.0) 566.8–560 (565.8,567.9) 559.5–554 (558.5,560.6) 559.8 245.6 72 537

c

558.0–542 574.4–564 558.2–541 568.3–557 564.5–552 568.1–558 580.0–563 577.9–567 579.2–568 574.1–568 570.3 351.3 70 217

0.15 (556.3,559.6) (572.6,576.2) (556.6,559.9) (566.6,570.0) (562.7,566.2) (566.5,569.7) (578.3,581.8) (575.8,580.0) (577.4,581.1) (572.2,575.9)

577.9–551 581.0–556 577.2–550 577.5–549 582.4–554 585.6–565 569.8–543 590.6–560 580.5–554 569.9–542 579.2 440.7 68 403

NS: number of solutions evaluated without being simulated (on the average). a Point estimation of the expected makespan. b Lower bound on the expected makespan. c Interval estimation of the expected makespan.

0.2 (575.8,580.0) (578.8,583.3) (574.6,579.7) (574.8,580.1) (580.2,584.6) (583.1,588.2) (567.6,572.1) (588.2,593.0) (578.3,582.8) (567.5,572.2)

598.9–556 587.9–547 606.1–563 608.4–574 585.9–552 616.0–576 585.9–544 607.3–560 614.6–567 591.0–555 600.2 1008.0 55 901

0.25 (595.5,602.3) (584.8,591.0) (603.1,609.0) (604.8,612.1) (582.1,589.6) (612.3,619.7) (582.8,589.1) (603.7,610.9) (611.4,617.9) (587.9,594.1)

604.2–550 614.1–563 623.3–575 607.8–547 609.4–555 624.9–573 620.0–562 629.4–574 617.7–567 604.4–557 615.5 1232.4 48 761

0.3 (600.3,608.0) (609.9,618.4) (619.6,627.0) (603.9,611.7) (605.7,613.0) (620.2,629.6) (616.2,623.8) (625.6,633.3) (613.5,621.8) (600.4,608.4)

626.2–560 643.6–568 640.1–571 632.3–562 634.5–567 631.2–577 643.5–575 625.3–567 634.5–572 624.4–560 633.6 1654.1 38 858

(621.6,630.7) (638.7,648.6) (635.4,644.9) (628.0,636.6) (629.9,639.1) (626.6,635.8) (638.8,648.1) (620.4,630.2) (630.1,638.9) (620.4,628.5)

160

F. Ahmadizar et al. / Computers & Operations Research 37 (2010) 152 -- 162

Table 2 Computational results for the normal distribution with  = 0.05. Trial

1 2 3 4 5 6 7 8 9 10 Average Time (s) NS

Variability 0.05

0.1

558.5a –555b (557.8,559.2)c 554.6–550 (553.8,555.3) 565.0–559 (564.3,565.8) 561.5–555 (560.9,562.2) 558.7–555 (557.9,559.5) 555.7–552 (555.0,556.4) 562.4–558 (561.7,563.0) 550.5–547 (549.8,551.2) 564.0–562 (563.2,564.7) 566.8–561 (566.2,567.4) 559.8 225.7 73 184

564.2–552 559.1–548 561.0–552 569.1–553 570.6–553 573.0–561 564.2–553 568.7–555 564.3–551 571.2–553 566.5 349.2 70 221

0.15 (562.8,565.5) (557.6,560.5) (559.7,562.3) (567.7,570.4) (569.1,572.1) (571.7,574.2) (562.8,565.5) (567.4,570.0) (563.1,565.4) (570.0,572.3)

574.5–551 576.0–553 564.7–540 591.4–568 572.9–549 578.8–553 578.9–552 590.1–570 588.6–571 579.3–554 579.5 525.7 65 912

0.2 (572.7,576.4) (574.2,577.8) (562.9,566.4) (589.5,593.4) (571.0,574.8) (576.9,580.6) (576.9,580.9) (587.9,592.3) (586.6,590.5) (577.5,581.2)

591.0–561 599.3–561 602.5–567 597.1–555 610.9–570 594.2–561 617.0–578 597.2–561 601.6–558 613.4–570 602.4 744.8 61 268

0.25 (588.6,593.5) (596.9,601.7) (599.9,605.0) (594.6,599.6) (608.5,613.3) (591.6,596.8) (614.5,619.5) (594.9,599.5) (599.1,604.1) (610.6,616.1)

613.4–559 632.9–579 611.5–565 628.4–578 628.1–572 619.6–564 611.7–556 619.7–577 614.6–565 610.8–557 619.1 1116.4 52 465

0.3 (610.4,616.3) (629.8,636.0) (608.6,614.5) (625.4,631.3) (625.2,631.1) (616.1,623.2) (609.0,614.5) (616.6,622.9) (611.6,617.5) (607.7,613.9)

640.3–572 623.9–559 646.2–576 623.9–549 627.7–557 635.4–565 633.9–563 643.8–575 635.7–570 641.8–571 635.3 1690.2 38 038

(636.8,643.7) (620.6,627.2) (642.7,649.7) (620.3,627.4) (624.1,631.3) (631.8,639.0) (630.3,637.5) (640.1,647.5) (631.6,639.8) (638.5,645.0)

NS: number of solutions evaluated without being simulated (on the average). a Point estimation of the expected makespan. b Lower bound on the expected makespan. c Interval estimation of the expected makespan. Table 3 Computational results for the uniform distribution with  = 0.01. Trial

Variability 0.05

1 2 3 4 5 6 7 8 9 10 Average Time (s) NS

0.1 a

b

562.9 –557 (562.1,563.8) 554.7–550 (553.7,555.7) 559.8–554 (559.0,560.6) 543.4–537 (542.6,544.3) 547.4–544 (546.6,548.3) 558.4–553 (557.6,559.3) 547.9–544 (547.0,548.7) 560.4–554 (559.6,561.3) 562.0–557 (561.1,562.9) 551.2–546 (550.3,552.1) 554.8 186.5 73 278

c

582.5–565 561.8–542 572.9–558 582.5–568 569.4–552 583.9–563 584.0–569 580.9–569 578.9–566 565.5–548 576.2 241.3 70 171

0.15 (580.7,584.3) (560.2,563.4) (571.2,574.6) (580.6,584.3) (567.7,571.2) (582.3,585.6) (582.4,585.5) (579.2,582.6) (577.0,580.8) (564.1,567.0)

577.2–558 585.3–555 584.3–554 573.6–552 583.7–564 585.2–563 589.6–563 594.9–570 582.2–552 585.9–560 584.2 316.6 65 543

0.2 (574.8,579.5) (582.9,587.8) (581.9,586.7) (571.3,575.9) (581.5,586.0) (582.8,587.6) (587.1,592.1) (592.4,597.5) (579.9,584.6) (583.6,588.3)

591.1–549 597.9–574 600.6–566 601.4–559 594.8–557 599.6–562 601.1–558 589.0–552 604.4–557 599.6–567 598.0 407.8 59 749

0.25 (588.4,593.8) (594.4,601.5) (597.5,603.8) (598.4,604.5) (591.8,597.9) (596.8,602.5) (597.9,604.2) (585.9,592.1) (601.2,607.6) (596.4,602.8)

616.2–572 613.3–569 604.2–547 615.8–566 604.4–568 608.0–555 618.8–560 622.8–570 623.3–571 609.3–563 613.6 667.9 46 583

0.3 (612.7,619.8) (609.4,617.1) (601.0,607.5) (612.0,619.6) (601.0,607.8) (604.5,611.5) (615.7,622.0) (618.3,627.3) (619.0,627.6) (605.6,613.0)

637.1–570 628.2–567 632.1–558 628.7–563 614.9–548 622.4–559 637.2–573 639.9–565 645.4–578 636.2–579 632.2 797.5 38 170

(632.5,641.7) (623.5,632.9) (628.2,635.9) (624.5,632.9) (609.9,619.9) (617.5,627.3) (632.8,641.7) (635.7,644.1) (640.6,650.1) (631.5,641.0)

NS: number of solutions evaluated without being simulated (on the average). a Point estimation of the expected makespan. b c

Lower bound on the expected makespan. Interval estimation of the expected makespan.

Table 4 Computational results for the uniform distribution with  = 0.05. Trial

Variability 0.05

1 2 3 4 5 6 7 8 9 10 Average Time (s) NS

0.1 a

b

547.6 –545 (546.8,548.4) 566.2–566 (565.4,567.0) 553.3–549 (552.6,553.9) 565.5–561 (564.7,566.3) 560.4–557 (559.7,561.1) 578.0–573 (577.4,578.7) 561.6–556 (560.9,562.3) 564.0–560 (563.3,564.7) 536.6–530 (535.9,537.3) 563.6–558 (562.8,564.4) 559.7 190.5 73 008

c

576.7–564 545.4–531 577.3–562 562.4–551 585.1–574 571.2–560 579.5–564 567.7–558 570.4–555 578.8–567 571.5 233.1 69 893

0.15 (575.2,578.3) (544.2,546.6) (576.0,578.5) (561.1,563.8) (583.7,586.4) (570.0,572.5) (578.3,580.7) (566.3,569.2) (569.3,571.4) (577.4,580.2)

582.8–563 584.9–553 572.5–541 573.7–551 573.8–546 582.9–566 589.7–571 577.3–553 596.4–575 566.5–540 580.1 322.8 65 349

0.2 (580.9,584.7) (582.9,586.9) (570.8,574.2) (571.8,575.6) (572.0,575.6) (581.1,584.7) (587.7,591.8) (575.7,578.9) (594.4,598.5) (564.7,568.2)

607.1–568 603.7–564 606.1–568 611.4–573 599.4–569 595.9–557 590.0–551 606.4–578 610.8–578 588.5–549 601.9 450.9 57 184

0.25 (605.1,609.2) (601.2,606.3) (603.8,608.4) (608.9,613.8) (597.0,601.7) (593.4,598.4) (587.7,592.3) (604.2,608.6) (608.2,613.4) (586.2,590.8)

616.0–559 592.3–543 612.5–555 606.3–552 608.4–560 611.5–572 611.2–558 613.4–564 608.4–555 613.6–565 609.4 470.0 56 116

0.3 (613.1,619.0) (589.8,594.8) (609.9,615.1) (603.5,609.2) (605.7,611.0) (608.1,614.9) (608.5,613.9) (610.2,616.5) (605.4,611.3) (610.7,616.6)

610.5–541 621.5–545 639.9–573 624.1–565 639.3–573 615.5–547 628.3–572 626.5–548 627.8–563 626.6–558 626.0 769.4 38 487

(607.4,613.5) (618.4,624.6) (636.3,643.5) (621.2,627.0) (635.8,642.7) (612.2,618.8) (625.2,631.4) (623.3,629.7) (624.8,630.7) (623.6,629.7)

NS: number of solutions evaluated without being simulated (on the average). a Point estimation of the expected makespan. b c

Lower bound on the expected makespan. Interval estimation of the expected makespan.

being simulated becomes smaller. On the other hand, with regard to Proposition 1, if the random variables have larger variances, the difference between the expected makespan and its lower bound is

also larger—as shown in Tables 1–4. This is due to the fact that as the variability grows, the Jensen's inequality holds in a stronger sense, that is, E[Max( . . . )] becomes much larger than Max[E( . . . )].

F. Ahmadizar et al. / Computers & Operations Research 37 (2010) 152 -- 162 Table 5 Computational results for the exponential distribution. Trial

1 2 3 4 5 6 7 8 9 10 Average Time (s) NS

 0.01

0.05

908.9a –564b (889.7,928.1)c 943.8–598 (922.3,965.2) 950.5–603 (928.5,972.5) 934.2–588 (913.3,955.2) 932.5–599 (911.8,953.1) 920.6–583 (897.6,943.5) 908.8–577 (889.2,928.4) 938.1–584 (917.0,959.1) 957.8–612 (936.7,978.9) 917.7–589 (895.8,939.7) 931.3 1752.1 0

938.9–582 916.9–587 934.7–604 944.1–613 913.1–615 936.4–580 945.7–600 939.3–611 941.1–607 928.9–595 933.9 1775.8 0

(921.3,956.4) (901.9,932.0) (917.5,951.8) (926.0,962.2) (895.3,930.9) (919.7,953.1) (926.8,964.7) (921.9,956.7) (922.2,959.9) (911.4,946.4)

NS: number of solutions evaluated without being simulated (on the average). a Point estimation of the expected makespan. b c

Lower bound on the expected makespan. Interval estimation of the expected makespan.

Considering the results shown in Table 5, it is seen that since the standard deviation of an exponential random variable is equal to its mean, that is, the variability is equal to 1, all constructed solutions have been evaluated via the simulation model. As a result if the variability is less important, a greater number of solutions can be constructed and evaluated in the same CPU time. The results show that under the same conditions as the level of confidence decreases, the average width of the confidence interval of the expected makespan decreases, which is not surprising; while decreasing the level of confidence does not have a significant effect on the average expected makespan, CPU time and number of solutions evaluated without being simulated. Furthermore, considering the results shown in Tables 1–4, it can be seen that under the same conditions, changing the distribution of the release dates and processing times makes no clear effect. This could be due to the fact that in a complete solution the number of independent random variables forming the longest path from node Se to node En is large enough and hence, from the central limit theorem, the sum of these random variables, i.e., the makespan, has a probability distribution that is approximately normal, regardless of their distribution. Therefore, one may conclude that if for the normal distribution as well as the uniform distribution the ratio of the standard deviation to the mean is selected as 1, the results must be compatible to those shown in Table 5. Of course, because it takes a great amount of time to generate a random number from a normal distribution, under the same conditions the average CPU times presented by Tables 1 and 2 are longer than those presented by Tables 3 and 4. 6. Summary and future work This paper deals with a stochastic group shop scheduling problem. The group shop scheduling problem is a general formulation that frequently occurs in industrial environments and includes the other shop scheduling problems. Manufacturing systems are subject to uncertain events which may be due to both human and machine resource factors and hence, considering a scheduling problem in a stochastic situation is more practical than in a deterministic situation. We study the stochastic group shop scheduling problem subject to stochastically distributed release dates and processing times, governed by known distribution functions, where the objective is to find a job schedule in the class of static policies to minimize the expected makespan. The problem is first formulated in a form of stochastic programming and then a lower bound on the expected makespan is

161

proposed. To solve the problem efficiently, a simulation optimization approach is developed that is a combination between an ant colony optimization algorithm and a heuristic algorithm to generate high-quality solutions and a discrete event simulation model to evaluate the expected makespan. Considering the lower bound on the expected makespan, a generated solution that may improve the current best one is specified; the time-consuming simulation model is used only for evaluating the performance of such a solution. In other words, the lower bound on the expected makespan may be used as a measure for evaluating the performance of a solution without simulating, particularly when the variability of the parameters is not very high. The proposed approach is tested on large instances where the random variables are normally, exponentially or uniformly distributed. The computational results are given to demonstrate the suitable performance of the approach. Several extensions to this study could be investigated in future research. It would be interesting to apply other metaheuristics such as genetic algorithm and simulated annealing to generate good solutions to the problem; however, in such situations the lower bound on the expected criterion may not be employed in the same manner. Another extension would be to consider other objective functions, in particular those are related to the due dates of jobs. Moreover, the study of bicriteria and multicriteria stochastic group shop scheduling problems would be useful. References [1] Pinedo M. Scheduling: theory, algorithms and systems. 2nd ed, Englewood Cliffs, NJ: Prentice-Hall; 2002. [2] Singer M. Forecasting policies for scheduling a stochastic due date job shop. Int J Prod Res 2000;38(15):3623–37. [3] Golenko-Ginzburg D, Kesler Sh, Landsman Z. Industrial job-shop scheduling with random operations and different priorities. Int J Prod Econ 1995;40: 185–95. [4] Golenko-Ginzburg D, Gonik A. Using ``look-ahead'' techniques in job-shop scheduling with random operations. Int J Prod Econ 1997;50:13–22. [5] Golenko-Ginzburg D, Gonik A. Optimal job-shop scheduling with random operations and cost objectives. Int J Prod Econ 2002;76:147–57. [6] Yoshitomi Y. A genetic algorithm approach to solving stochastic job-shop scheduling problems. Int Trans Oper Res 2002;9:479–95. [7] Yoshitomi Y, Yamaguchi R. A genetic algorithm and the Monte Carlo method for stochastic job-shop scheduling. Int Trans Oper Res 2003;10:577–96. [8] Tavakkoli-Moghaddam R, Jolai F, Vaziri F, Ahmed PK, Azaron A. A hybrid method for solving stochastic job shop scheduling problems. Appl Math Comput 2005;170:185–206. [9] Lai T-C, Sotskov YN, Sotskova NY, Werner F. Optimal makespan scheduling with given bounds of processing times. Math Comput Modell 1997;26(3):67–86. [10] Lai T-C, Sotskov YN. Sequencing with uncertain numerical data for makespan minimization. J Oper Res Soc 1999;50:230–43. [11] Lai T-C, Sotskov YN, Sotskova N, Werner F. Mean flow time minimization with given bounds of processing times. Eur J Oper Res 2004;159:558–73. [12] Luh PB, Chen D, Thakur LS. An effective approach for job-shop scheduling with uncertain processing requirements. IEEE Trans Robot Autom 1999;15(2): 328–39. [13] Neumann K, Schneider WG. Heuristic algorithms for job-shop scheduling problems with stochastic precedence constraints. Ann Oper Res 1999;92: 45–63. [14] Alcaide D, Rodriguez-Gonzalez A, Sicilia J. A heuristic approach to minimize expected makespan in open shops subject to stochastic processing times and failures. Int J Flex Manuf Syst 2006;17(3):201–26. [15] http://www.win.tue.nl/whizzkids/1997. [16] Blum C, Sampels M. An ant colony optimization algorithm for shop scheduling problems. J Math Modell Algorithms 2004;3:285–308. [17] Liu SQ, Ong HL, Ng KM. A fast tabu search algorithm for the group shop scheduling problem. Adv Eng Software 2005;36:533–9. [18] Roy B, Sussmann B. Les problemes d'ordonnancement avec constraintes disjonctives. Technical report note, DS 9, SEMA, Paris, 1964. [19] Bellman RE. On a routing problem. Quart Appl Math 1958;16:87–90. [20] Hardy GH, Littlewood JE, Polya G. Inequalities. 2nd ed, Cambridge: Cambridge University Press; 1952. [21] Ho YC. On the numerical solution of stochastic optimization problem. IEEE Trans Autom Control 1997;42(5):727–9. [22] Rubinstein RY. Monte Carlo optimization, simulation and sensitivity of queuing networks. New York: Wiley; 1986. [23] Dorigo M. Optimization, learning and natural algorithm. PhD thesis, DEI, Politecnico di Milano, Italy, 1992 [in Italian]. [24] Dorigo M, Maniezzo V, Colorni A. The ant system: optimization by a colony of cooperating agents. IEEE Trans Syst Man Cybernet B 1996;26(1):29–41.

162

F. Ahmadizar et al. / Computers & Operations Research 37 (2010) 152 -- 162

[25] Dorigo M, Stutzle T. Ant colony optimization. Cambridge: MIT Press; 2004. [26] Dorigo M, Blum C. Ant colony optimization theory: a survey. Theoret Comput Sci 2005;344:243–78. [27] Blum C, Sampels M. Ant colony optimization for FOP shop scheduling: a case study on different pheromone representations. In: Proceedings of the 2002 congress on evolutionary computation (CEC'02), vol. 2. Los Alamitos, CA: IEEE Computer Society Press; 2002. p. 1558–63. [28] Blum C, Dorigo M. The hyper-cube framework for ant colony optimization. IEEE Trans Syst Man Cybernet B 2004;34(2):1161–72.

[29] Stutzle T, Hoos HH. MAX–MIN ant system. Future Gen Comput Syst 2000;16(8):889–914. [30] Giffler B, Thompson GL. Algorithms for solving production scheduling problems. Oper Res 1960;8(4):487–503. [31] Dorigo M, Gambardella LM. Ant colony system: a cooperative learning approach to the traveling salesman problem. IEEE Trans Evol Comput 1997;1(1): 53–66. [32] Banks J, Carson JS, Nelson BL. Discrete-event system simulation. 2nd ed, Englewood Cliffs, NJ: Prentice-Hall; 1996.