An Artificial Bee Colony Algorithm for the Job Shop Scheduling Problem with Random Processing Times

Due to the influence of unpredictable random events, the processing time of each operation should be treated as random variables if we aim at a robust production schedule. However, compared with the extensive research on the deterministic model, the stochastic job shop scheduling problem (SJSSP) has not received sufficient attention. In this paper, we propose an artificial bee colony (ABC) algorithm for SJSSP with the objective of minimizing the maximum lateness (which is an index of service quality). First, we propose a performance estimate for preliminary screening of the candidate solutions. Then, the K-armed bandit model is utilized for reducing the computational burden in the exact evaluation (through Monte Carlo simulation) process. Finally, the computational results on different-scale test problems validate the effectiveness and efficiency of the proposed approach.


Introduction
The job shop scheduling problem (JSSP) has been known as an extremely stubborn combinatorial optimization problem since the 1950s.In terms of computational complexity, JSSP is N P-hard in the strong sense [1].It models an important decision problem in contemporary manufacturing systems, so extensive research has been conducted on JSSP.Due to its high complexity, the recent research has focused on the meta-heuristic approaches, such as genetic algorithm (GA) [2], estimation of distribution algorithm (EDA) [3], tabu search (TS) [4], particle swarm optimization (PSO) [5], ant colony optimization (ACO) [6], artificial bee colony (ABC) algorithm [7].
However, most existing research is based on the standard JSSP model, which means it is impossible to directly apply these algorithms in real-world scheduling scenarios.Especially, we emphasize the following two points.
(1) Most research on JSSP has focused on the makespan criterion (i.e., minimizing the maximum completion time).However, in the make-to-order (MTO) manufacturing environment, due date related performances are apparently more relevant for decision makers, because the in-time delivery of goods is vital for maintaining a high service reputation.Therefore, the research that aims at minimizing lateness/tardiness in JSSP deserves more attention.(2) Most existing algorithms are designed for the deterministic JSSP, in which all the data (e.g., the processing times) are assumed to be fixed and precisely known in advance.In real-world manufacturing, however, the processing of operations is constantly affected by uncertain factors.Machine breakdowns, worker absenteeism, order changes, etc. can all lead to variations in the operation times.In this case, solving the deterministic JSSP will not result in a robust production schedule.Therefore, it is rewarding to focus more research effort on the stochastic JSSP.
In an attempt to overcome the drawbacks, in this paper we study the stochastic job shop scheduling problem (SJSSP) with the objective of minimizing maximum lateness (L max ) in the sense of expectation.Compared with the standard JSSP, the studied SJSSP moves a step closer to the practical scheduling features.The main contribution of this paper is an efficient optimization approach based on artificial bee colony for solving SJSSP.The algorithm first uses a quick-and-dirty performance estimate to roughly evaluate the candidate solutions.Then, the selected promising solutions undergo a more accurate evaluation through simulation.In this process, a new computing budget allocation policy is adopted in order to promote the computational efficiency.
The rest of the paper is organized as follows.Section 2 briefly reviews the recent publications on stochastic job shop scheduling and artificial bee colony algorithms.Section 3 provides the formulation of SJSSP and an introduction to the standard ABC principles.Section 4 introduces an estimate for the objective function (maximum lateness).Section 5 describes the design of ABC for SJSSP in detail.Section 6 gives the computational results.Finally, some conclusions are made in Section 7.

The Stochastic Job Shop Scheduling Problem
Due to the inevitable uncertainties in manufacturing systems, SJSSP is more realistic than its deterministic counterpart.In SJSSP, the number of jobs is usually known in advance, while the processing time of each operation is a random variable with known probability distribution.The due dates can be regarded as either fixed or random variables, depending on the frequency of changes in customer orders.The aim is to find a feasible schedule (i.e., sequence of operations) that minimizes the objective function in the sense of expectation.
In [8,9], simulation-based GAs are proposed for solving SJSSP with the expected makespan criterion.The individuals that appear through the generations with very high frequency are selected as good solutions.In [10], the authors study the SJSSP with just-in-time objective (completion before or after the due date will both result in a cost).Several decision-making rules are proposed for selecting a job when several jobs are competing for a machine.In [11], the goal is to minimize the processing time variations, the operational costs and the idle costs in SJSSP.A hybrid method is proposed, in which the initial solutions are first generated by a neural network and then improved by SA.Luh et al. [12] presents an algorithm based on Lagrangian relaxation and stochastic dynamic programming for SJSSP with uncertain arrival times, due dates, and part priorities.In [13], exact and heuristic algorithms are proposed to solve SJSSP with mean flow time criterion.The aim is to find a minimum set of schedules which contains at least one optimal schedule for any realization of the random processing times.Singer [14] presents a heuristic which amplifies the expected processing times by a factor and then applies a deterministic scheduling algorithm.Recently, Lei proposes a genetic algorithm for minimizing makespan in a stochastic job shop with machine breakdowns and non-resumable jobs [15].The processing time, the downtime and the normal operation time between successive breakdowns are all assumed to follow exponential distributions.In addition, a multi-objective genetic algorithm is proposed in [16] for stochastic job shop scheduling problems in which the makespan and the total tardiness ratio should be minimized.In [17,18], quantum genetic algorithms are proposed to solve SJSSP with expected makespan criterion.In [19], a simulation-based decision support system is presented for the production control of a stochastic flexible job shop manufacturing system.In [20], an algorithm based on computer simulation and artificial neural networks is proposed to select the optimal dispatching rule for each machine from a set of rules in order to minimize the makespan in SJSSP.
Generally, the most critical factor that affects the efficiency of solving a stochastic combinatorial optimization problem (SCOP) is the evaluation of solutions [21].Since the objective function is in the expectation form, the evaluation may not be so straightforward if no closed-form formula is available for calculating the expectation.Due to the complexity of the job shop model, we have to rely on Monte Carlo simulations to obtain an approximation for the expected L max .However, the simulation process is usually very time-consuming.To achieve a balance, we propose an artificial bee colony algorithm for solving SJSSP in this paper.Time-saving strategies have been designed and utilized in the searching process in order to ensure a satisfactory solution within reasonable computational time.

The Artificial Bee Colony Algorithm
The artificial bee colony (ABC) algorithm is a relatively new swarm intelligence based optimizer.It mimics the cooperative foraging behavior of a swarm of honey bees [22].ABC was initially proposed by Karaboga in 2005 for optimizing multi-variable and multi-modal continuous functions [23].The latest research has revealed some good properties of ABC [24][25][26].Especially, the number of control parameters in ABC is fewer than that of other population-based algorithms, which makes it easier to be implemented.Meanwhile, the optimization performance of ABC is comparable and sometimes superior to the state-of-the-art meta-heuristics.Therefore, ABC has aroused much interest and has been successfully applied to different kinds of optimization problems [27][28][29].
The ABC algorithm systematically incorporates exploration and exploitation mechanisms, so it is suitable for complex scheduling problems.For example, Huang and Lin [30] proposed a new ABC for the open shop scheduling problem, which is more difficult to solve than JSSP due to the undecided precedence relations between job operations.In this paper, we will use ABC as the basic optimization framework in the process of solving SJSSP.To our knowledge, this is the first attempt that ABC is applied to a stochastic scheduling problem.Of course, some problem-specific information should be embedded into the searching process of ABC in order to promote the overall optimization efficiency for SJSSP.

Formulation of the SJSSP
In an SJSSP instance, a set of n jobs {J j } n j=1 are to be processed on a set of m machines {M k } m k=1 under the following basic assumptions: (1) there is no machine breakdown; (2) no preemption of operations is allowed; (3) all jobs are released at time 0; (4) the transportation times and the setup times are all neglected; (5) each machine can process only one job at a time; (6) each job may be processed by only one machine at a time.
Each job has a fixed processing route which traverses all the machines in a predetermined order.The manufacturing process of job j on machine k is noted as operation O jk .Besides, a preset due date d j (describing the level of urgency) is given for each job j.The duration of an operation is influenced by many real-time factors such as the condition of workers and machines.Therefore, in the scheduling stage, the processing time of each operation is usually not known exactly.We assume that the processing times (p jk , for each O jk ) are independent random variables with known expectation (E (p jk )) and variance (var (p jk )).The objective function is the expected maximum lateness (L max ).If we use C j to denote the completion time of job j, then L max is defined as max n j=1 {C j − d j }.We consider L max because due date related performances are becoming very significant in the make-to-order manufacturing environment nowadays.
Like its deterministic counterpart, SJSSP can also be described by a disjunctive graph G(O, A, E) [31].O = {O jk |j = 1, . . ., n, k = 1, . . ., m} is the set of nodes.A is the set of conjunctive arcs which connect successive operations of the same job, so A describes the technological constraints in the SJSSP instance.E = ∪ m k=1 E k is the set of disjunctive arcs, where E k denotes the disjunctive arcs corresponding to the operations on machine k.Each arc in E k connects a pair of operations to be processed by machine k and ensures that the two operations should not be processed simultaneously.Initially, the disjunctive arcs do not have fixed directions, unlike the conjunctive arcs in A.
Under the disjunctive graph representation, finding a feasible schedule for the SJSSP is equivalent to orienting all the disjunctive arcs so that no directed cycles exist in the resulting graph.In this paper, we use σ to denote the set of directed disjunctive arcs which is transformed from the original E. Thus, if A ∪ σ is acyclic, the schedule corresponding to σ is feasible [32].

  
where the k-th row specifies the resolved disjunctive arcs related with machine k (k = 1, 2, 3).
Based on the disjunctive graph model, the discussed SJSSP can be mathematically formulated as follows: In this formulation, (x) + = max{x, 0}.t jk represents the starting time of operation O jk .k j denotes the index of the machine that processes the last operation of job j, so the completion time of job j is t jk j + p jk j .Constraint (a) ensures that the processing order of the operations from each job is consistent with the technological routes.Constraint (b) ensures that the processing order of the operations on each machine complies with the sequence specified by the schedule, σ.
When a feasible schedule σ is given, a minimum L σ max can be achieved for each realization of {p jk }.Since {p jk } are random variables, {t jk } and L σ max are also random variables.Therefore, the objective function is expressed in the form of expectation, and constraints (a-c) should hold almost everywhere (a.e.).The aim of solving SJSSP is to find a schedule σ with the minimum E (L σ max ).

Principles of the Artificial Bee Colony (ABC) Algorithm
In the ABC algorithm, the artificial bees are classified into three groups: the employed bees, the onlookers and the scouts.A bee that is exploiting a food source is classified as employed.The employed bees share information with the onlooker bees, which are waiting in the hive and watching the dances of the employed bees.The onlooker bees will then choose a food source with probability proportional to the quality of that food source.Therefore, good food sources attract more bees than the bad ones.Scout bees search for new food sources randomly in the vicinity of the hive.When a scout or onlooker bee finds a food source, it becomes employed.When a food source has been fully exploited, all the employed bees associated with it will abandon the position, and may become scouts again.Therefore, scout bees perform the job of "exploration", whereas employed and onlooker bees perform the job of "exploitation".In the algorithm, a food source corresponds to a possible solution to the optimization problem, and the nectar amount of a food source corresponds to the fitness of the associated solution.
In ABC, the first half of the colony consist of employed bees and the other half are onlookers.The number of employed bees is equal to the number of food sources (SN ) because it is assumed there is only one employed bee for each food source.Thus, the number of onlooker bees is also equal to the number of solutions under consideration.The ABC algorithm starts with a group of randomly generated food sources.The main procedure of ABC can be described as follows.
Step 1: Initialize the food sources.
Step 2: Each employed bee starts to work on a food source.
Step 3: Each onlooker bee selects a food source according to the nectar information shared by the employed.
Step 4: Determine the scout bees, which will search for food sources in a random manner.
Step 5: Test whether the termination condition is met.If not, go back to Step 2.
The detailed description for each step is given below.
(1) The initialization phase.The SN initial solutions are randomly-generated D-dimensional real vectors.Let x i = {x i,1 , x i,2 , . . ., x i,D } represent the i-th food source, which is obtained by where r is a uniform random number in the range [0, 1], and x min (2) The employed bee phase.In this stage, each employed bee is associated with a solution.She exerts a random modification on the solution (original food source) for finding a new solution (new food source).This implements the function of neighborhood search.The new solution v i is generated from x i using a differential expression: where d is randomly selected from {1, . . ., D}, k is randomly selected from {1, . . ., SN } such that k ̸ = i, and r ′ is a uniform random number in the range Once v i is obtained, it will be evaluated and compared to x i .If the fitness of v i is better than that of x i (i.e., the nectar amount of the new food source is higher than the old one), the bee will forget the old solution and memorize the new one.Otherwise, she will keep working on x i .(3) The onlooker bee phase.When all employed bees have finished their local search, they share the nectar information of their food source with the onlookers, each of whom will then select a food source in a probabilistic manner.The probability p i by which an onlooker bee chooses food source x i is calculated as follows: Obviously, the onlooker bees tend to choose the food sources with higher nectar amount.
Once the onlooker has selected a food source x i , she will also conduct a local search on x i according to Equation (2).As in the previous case, if the modified solution has a better fitness, the new solution will replace x i .(4) The scout bee phase.In ABC, if the quality of a solution cannot be improved after a predetermined number (limit) of trials, the food source is assumed to be abandoned, and the corresponding employed bee becomes a scout.The scout will then produce a food source randomly by using Equation (1).
To facilitate the understanding of the algorithm, a flow chart [33] is provided as Figure 2.

An Estimate for the Expected Maximum Lateness
In the studied SJSSP, the objective function is expressed as an expectation because the processing times are random variables.Meanwhile, due to the N P-hardness of JSSP, there does not exist a closed-form expression to calculate the expected objective value.Therefore, the evaluation of a given schedule is not a simple task.Usually, we can utilize the idea of Monte Carlo simulation, i.e., taking the average objective value in a large number of realizations as an estimate for the expectation.However, this definitely increases the computational burden, especially when used in an optimization framework (frequent solution evaluations are needed).
If there is no strict requirement on the evaluation accuracy, it is natural to think of a time-saving strategy: cut down the number of realizations.Furthermore, if we allow only one realization of the random processing times, then a rational choice is to use the mean (expected) value of each processing time.Then, we can show the following property, that is, such an estimate is a lower bound of the true objective value.
Theorem 1.Let σ denote a feasible schedule of the considered SJSSP instance.The following inequality must hold: where L σ max (random variable) is the maximum lateness corresponding to the schedule, and L σ max (constant value) is the maximum lateness in the case where each random processing time takes the value of its expectation.Proof.For the sake of simplicity, we will omit the superscript "σ" in the following proof because we are already focusing on the given schedule.We will denote the completion time of operation O jk by c jk (a random variable).Furthermore, let t jk (resp.c jk ) denote the starting time (resp.completion time) of operation O jk when each processing time is replaced by its expected value.Since A ∪ σ is acyclic, all the operations can be sorted topologically, yielding an ordered sequence {O In this sequence, the machine predecessors and the job predecessors of any operation O jk must be positioned before O jk .First we will prove that, for any operation O jk , E (t jk ) ≥ t jk .
We start from the first operation (i = 1) in the sequence.Because O [1] jk has no predecessor operations, we have E (t jk ) = t jk (actually t jk = 0).Then, the proof procedure will continue for i = 2, . . ., n × m.Suppose we have already proved E (t jk ) ≥ t jk for each operation before jk in the sequence, and without loss of generality, we assume that O jk has an immediate machine predecessor O j ′ k and an immediate job predecessor O jk ′ , then we have If O jk has no job predecessor, then we set c jk ′ = 0 in the above proof.Likewise, if O jk has no machine predecessor, then we set c j ′ k = 0. Therefore, the reasoning applies to each operation in the sequence.
Having proved E (t jk ) ≥ t jk , we can now move to the objective function (L max ).Let C j (resp.L j ) denote the completion time (resp.lateness) of job j when each random processing time takes its expected value.Meanwhile, k j represents the machine which processes the last operation of job j.Then, This completes the proof of E (L σ max ) ≥ L σ max .

The Proposed ABC Algorithm for Solving SJSSP
In this section, the proposed ABC is introduced in detail.First, we describe how ABC can be adapted to a discrete optimization problem like SJSSP.Second, we describe the method for comparing the quality (subject to estimation errors) of different solutions in the proposed ABC.Since the problem is stochastic, in the solution comparison process, we must decide how to allocate the available simulation replications.This problem is resolved by the model proposed in the third subsection.Finally, to facilitate the utilization of the simulation allocation model, we slightly revise the implementation of the local search procedure in standard ABC.

Adaptation to the Discrete Problem
Since the standard ABC is intended for continuous function optimization, we need to introduce some modifications in order to adapt it to SJSSP.

Encoding and Decoding
First, the encoding scheme for SJSSP is based on operation sequences instead of real number vectors.In particular, each solution (food source) is represented by a sequence of n × m operations, which can be transformed into a feasible schedule by the "schedule builder" presented below.This procedure functions by iteratively scanning the sequence and then picking out the first schedulable operation [34] in the sequence.In order to build an active schedule, each operation is inserted into the earliest possible idle period (i.e., time interval not occupied by the previously scheduled operations).
Input: A sequence of operations, π, with a length of n × m.
Step 1: Let σ be an empty matrix of size m × n.
Step 2: If π = ∅, output σ (and the corresponding L max if necessary) and terminate the procedure.
Otherwise, continue the following steps.Step 6: Delete operation O * from π. Go back to Step 2.

Initialization
Second, the solution initialization policy (Equation ( 1) in standard ABC) should be modified to comply with the solution encoding scheme.Here we use the famous ATC (apparent tardiness cost) rule [35] to produce the initial population.To generate the u-th solution, we first make a sample of the random processing times according to their distributions, denoted as {p (u) jk }.Then, the ATC rule is applied in a simulation-based scheduling procedure, where the operation with the largest Z jk value will be selected for processing when a machine is freed at time t.
where JS(O jk ) is the set of job successors of operation O jk .p (u) denotes the average processing time of the currently waiting operations in the machine's buffer.K is a scaling parameter (or called "look-ahead" parameter).W jk ′ is the estimated lead time of operation O jk ′ .Here it is assumed K = 2 and W jk ′ = 0.4p (u) jk ′ .

Neighborhood Structure
Third, the neighborhood search mechanism (Equation ( 2) in standard ABC) also needs to be modified accordingly.Let's recall that, in the case of deterministic JSSP, we must focus on the operations that belong to the blocks [36] of tardy jobs in order to improve the current solution.Therefore, the neighborhood search for SJSSP can be performed by using the expected values of processing times to determine the critical paths and blocks.Particularly, we first select two adjacent operations randomly from a block related with a tardy job, and then the SWAP operator is applied to exchange the positions of these two operations in the encoded solution.

The Comparison of Solutions
In order to approximate E (L σ max ), we need to implement the schedule σ under different realizations of the random processing times.When σ has been evaluated for a sufficient number of times, say τ , then its objective value can be approximated by 1 τ ∑ τ i=1 L σ max (i) (where L σ max (i) corresponds to the i-th realization).This is consistent with the idea of Monte Carlo simulation.
In the employed bee phase and the onlooker bee phase, the newly-generated solution must be compared with the original solution in order to determine which one should be kept.For deterministic optimization problems, this can be done by simply comparing the exact objective values of the two solutions.But in the stochastic case, the comparisons may not be so straightforward because we can only obtain approximated (noisy) objective values as mentioned above.In this study, we will utilize the following two mechanisms for comparison purposes.
(I) Pre-screening.Because L σ max is a lower bound for E (L σ max ) (Theorem 1), we can arrive at the following conclusion which is useful for the pre-screening of candidate solutions.Corollary 1.For two candidate solutions x 1 (the equivalent schedule is denoted by σ 1 ) and x 2 (the equivalent schedule is denoted by σ 2 ), if L σ 2 max ≥ E (L σ 1 max ), then x 2 must be inferior to x 1 and thus need not be considered.
(II) Hypothesis test.If the candidate solution has passed the pre-screening, then hypothesis test is used for comparing the quality of the two solutions.Suppose we have implemented n i simulation replications for solution x i whose true objective value is f Then, the sample mean and sample variance can be calculated by is the objective value obtained in the j-th simulation replication for solution x i .
Let the null hypothesis H 0 be "f (x 1 ) = f (x 2 )", and thus the alternative hypothesis H 1 is "f (x 1 ) ̸ = f (x 2 )".According to the statistical theory, the critical region of H 0 is where z α/2 is the value such that the area to its right under the standard normal curve is exactly α/2.Therefore, if f 1 − f 2 < Z (i.e., the null hypothesis holds), it is said that there exists no statistical difference between x 1 and x 2 .Otherwise, if

The Allocation of Simulation Replications
In the previous subsection, we assume that n i simulation replications have been performed for solution x i before the hypothesis test.But how should the value of n i be determined?We know that full-length simulation is very time-consuming, so a smaller value is desirable for n i .On the other hand, however, if n i is too small, the approximation for the objective values will be over-imprecise, which can mislead the hypothesis test.To achieve a balance, here we borrow an idea from the solution to the famous K-armed bandit problem.
The K-armed bandit problem can be described as K random variables, X i (1 ≤ i ≤ K), where X i represents the stochastic reward given by the i-th gambling machine.The distributions of X i are independent but generally not identical.The laws of the distributions and the expectations µ i for the rewards X i are unknown.The goal is to find a strategy that determines the next machine to play (based on the past plays and the received rewards) that maximizes the expected total reward for the player.
Since a player does not know which of the machines is the best, he can only guess the reward distributions from successive plays.Clearly, he has to make a trade-off between the following two choices: • The player tends to concentrate his efforts on the machines that gave the highest rewards in the past plays.Intuitively, this choice will maximize the potential gains.• The player also wants to try the machines which he has played very few times in the past.The reward distributions of these machines are quite unclear, but it is likely that they provide even higher rewards.
Therefore, it is necessary to strike a wise balance between exploiting the currently best machine and exploring the other machines to make sure that none of those is even better.
When allocating the simulation replications, we are actually facing a similar situation.There are a number of solutions waiting to be evaluated.For each solution, the objective value obtained from one simulation is a random variable following an unknown distribution.We want to find out the quality of these solutions with the minimum total number of simulation replications.The question is to decide which solution should get the next chance of simulation.Obviously, the solutions are analogous to the gambling machines, and the simulations are analogous to the gambler's tries.
Here we use an algorithm called UCB1 [38] for making the allocation decisions.If the number of candidate solutions is K, while the available computing resource is only enough to support a total of T replications of simulation (T > δK, where δ is the number of replications assigned at a time), then the algorithm can be adapted to our problem as follows.
Input: A set of K candidate solutions.In the adapted UCB1 (noted as A-UCB1) algorithm, ν k records the number of times that solution x k has been evaluated through simulation, while ν is the total number of simulation replications that have been assigned so far.The "reward" from trying x k is defined as the relative standard deviation of the simulated objective values.Under the effect of ρ, the definition of r k ensures that the simulation replications tend to be allocated to the solutions whose quality is still unclear (represented by the relatively large s k ).Meanwhile, the algorithm will also allocate simulation replications to the solutions which have been evaluated very few times (represented by the small ν k ).In sum, the A-UCB1 algorithm aims at promoting the computational efficiency in the evaluation of solutions.Based on extensive computational tests, the value of δ is set as ⌊T /100⌋.

Revised Implementation of the Local Search in ABC
In the employed bee phase and the onlooker bee phase, the artificial bees are actually performing the local search task.Considering the special features of the stochastic optimization problem, we slightly modify the implementation of this local search mechanism.In particular, we evaluate and compare the solutions group by group rather than one by one.The aim is to utilize the benefits of the A-UCB1 algorithm (for controlling the computational burden) and the hypothesis test (for increasing the diversity of the whole population).Our implementation is detailed as follows.Suppose the bees are already associated with each solution.
Step 1: Use A-UCB1 to allocate a total of T simulation replications to the solutions in the current population P .Save the estimated mean objective value and variance for each solution.
Step 2: Generate a new solution based on each original solution using the SWAP operator.Apply the pre-screening method (Corollary 1) to ensure the quality of each new solution [39].Denote the set of new solutions by P new .
Step 3: Use A-UCB1 to allocate a total of T simulation replications to the solutions in P new .Save the estimated mean objective value and variance for each solution.
Step 4: Perform hypothesis tests to determine the new population: (4.1) Sort all the solutions in P ∪ P new in non-decreasing order of the estimated objective values, yielding a sequence {x [1] , . . ., x [2SN ] }.Put x [1] into the ultimate population P * .Let i = 1, j = 2.
(4.2) Perform hypothesis test for x [j] and the i-th solution (the latest) in P * .If the null hypothesis holds, x [j] is rejected.Otherwise, x [j] is appended to P * and let i ← i + 1.
(4.4) If i = SN , the ultimate population has been determined.Otherwise, generate (SN − i) new solutions randomly, evaluate them (each one is evaluated with ⌊T /SN ⌋ simulation replications) and append them to P * .
Step 5: Compare P * with P and check whether each new solution has been accepted (i.e., in P * ).If accepted, let the corresponding bee fly to the new solution.
In the proposed ABC algorithm, the above procedure is used in place of the one-to-one greedy selection policy in standard ABC.This is in consideration of the specific features of stochastic optimization.Evaluating a group of solutions together is beneficial, because A-UCB1 will allocate more simulation replications to good solutions (with smaller f k ) to make sure about their true performance and avoid wasting time on inferior solutions (with larger f k ).The hypothesis test is used for maintaining an adequate level of diversity of the solution population.In Step 4.2, if x [j] does not significantly differ from the most recent solution in P * , it will not have a chance to enter P * .Finally, if the number of qualified solutions is less than SN , the rest of solutions will be generated randomly.

Generation of Test Instances
To test the effectiveness of the proposed algorithm, computational experiments are conducted on a number of randomly-generated test instances.In each instance, the route of each job is a random permutation of m machines.Three types of distribution patterns are considered for the random processing times: normal distribution, uniform distribution and exponential distribution.In all cases, the mean values (µ jk ) are generated from the uniform distribution U (1, 99).In the case of normal distributions (i.e., p jk ∼ N (µ jk , σ 2 jk )), the standard deviation is controlled by σ jk = θ × µ jk (θ describes the level of variability).In the case of uniform distributions (i.e., p jk ∼ U (µ jk − ω jk , µ jk + ω jk )), the width parameter is given by ω jk = θ × µ jk .In the case of exponential distributions (i.e., p jk ∼ Exp(λ jk )), the only parameter is given by λ jk = 1/µ jk .The due dates are obtained by a series of simulation runs which apply different priority rules (such as SPT, EDD, etc.) on each machine, and the due date of each job is finally set as its average completion time.This method can generate reasonably tight due dates.Meanwhile, the weight of each job is generated from the uniform distribution U (1, 10).The following computational experiments are conducted in Visual C++ 2010 on an Intel Core i5-750/3GB RAM/Windows 7 PC.

The Computational Results and Comparisons
Based on extensive computational tests (not listed due to space limitations), the settings of key parameters in the final ABC algorithm are: SN = 30, limit = 40 and T = 1000 (for each evaluation).The confidence level for the hypothesis test is α = 0.05.
In the following, we will use the proposed ABC to solve different-sized SJSSP instances.The results are compared with the hybrid optimization method PSO-SA [40], which utilizes SA (simulated annealing) as a local optimizer for PSO.In order to make the comparisons meaningful, we set a computational time limit for both algorithms.Normally, the time limit should be set according to the size of the instance.In our experiment, the time limit for solving a n × m instance is determined as t = 0.2 × (n × m).For example, the allowed computational time for a 20 × 10 instance is 40 sec.Each algorithm is run for 20 independent times on each SJSSP instance.
We report the best, average and worst objective values (exact evaluation [41] for the output solutions) in the 20 runs.Tables 1-3 report the results under normal distributions, Tables 4-6 report the results under uniform distributions, and Table 7 reports the results under exponential distributions.In addition, we have also performed Mann-Whitney U tests to statistically compare the computational results.Because each algorithm is run for 20 times, we have n 1 = n 2 = 20 in the Mann-Whitney U test.The null hypothesis is that there is no difference between the results of the two compared algorithms.Therefore, if the obtained U (the lesser of U 1 and U 2 ) is below 127, the null hypothesis can be rejected at the 5% level.Furthermore, if U < 105, the null can be rejected at the 1% level.The values of U are listed in Table 8.Based on the statistical tests, we can conclude that ABC is significantly more effective than the comparative method.In addition, the following comments can be made.
(1) According to the tables, the advantage of ABC over PSO-SA is greater when the variability level of processing times is higher (represented by larger θ or the case of exponential distribution).This can be attributed to the function of A-UCB1, which is responsible for the allocation of the limited computational resources.If θ is small, the objective value of a solution can be obtained with only a few replications.In this case, the A-UCB1 strategy is not significantly better than an equal allocation of the available replications (the case of PSO-SA).However, when the variability increases, the computational time becomes a relatively scarce resource.In order to correctly identify high-quality solutions, the limited replications should be allocated in an efficient manner rather than evenly.So in this case, the advantage of using A-UCB1 becomes evident.(2) According to the tables, ABC outperforms PSO-SA to a greater extent when solving the larger-scale instances.If the solution space is huge, PSO must rely on the additional local search module (SA) to promote the searching efficiency.However, SA is not so efficient as the inherent local search mechanism (employed and onlooker bees) in ABC, especially under tight time budgets.From another perspective, ABC systematically combines the exploration and exploitation abilities, and thus it works in a coordinated fashion.By contrast, the hybrid algorithm uses PSO for exploration and SA for exploitation, but the two algorithms have different search patterns.This may weaken the cooperation between PSO and SA in solving large-scale SJSSP.Therefore, ABC alone is more efficient than the hybrid algorithm.

Conclusions
In this paper, an artificial bee colony algorithm is proposed for solving the job shop scheduling problem with random processing times.The objective function is to minimize the expected maximum lateness.In view of the stochastic nature of the problem, two mechanisms are devised for the evaluation of solutions.First, a quick-and-dirty performance estimate is used to pre-screen the candidate solutions, so that the obviously inferior solutions can be eliminated at an early stage.Then, Monte Carlo simulation is applied to obtain a more accurate evaluation for the surviving solutions.In this process, a simulation budget allocation method is designed based on the K-armed bandit metaphor.This helps to utilize the limited computational time in an efficient manner.The computational results on a wide range of test instances reveal the superiority of the proposed approach.
The future research can be conducted from the following aspects: (1) It is worthwhile to consider other types of randomness in job shops, for example, the uncertainty in the processing routes of certain jobs.(2) It is worthwhile to investigate the method for discovering and utilizing problem-specific characteristics of SJSSP.This will make the ABC algorithm more pertinent to this problem.

Figure 1 (
a) shows the disjunctive graph for a 3 × 3 instance.The solid lines represent the conjunctive arcs while the dashed lines with bidirectional arrows express the disjunctive arcs.

Figure 1 .
Figure 1.A concrete example for the disjunctive graph representation of SJSSP.

d and x max d are the lower
and upper bounds for dimension d, respectively.

Figure 2 .
Figure 2. Flow chart of the ABC algorithm.

Step 3 :( 5 . 2 )
Find the first schedulable operation in the sequence π, denoted by O * .Step 4: Identify the machine required to process O * and denote it by k * .Record the expected processing time of O * as p * .Step 5: Schedule the operation O * : (5.1) Scan the Gantt chart of machine k * (which records the processing information of the already scheduled operations) from time zero and test whether O * can be inserted into each idle period [a, b], i.e., whether the following condition is met: max{a, C JP * } + p * ≤ b (where C JP * denotes the completion time of the immediate job predecessor of operation O * ).If the above inequality is satisfied for the idle interval between operation o 1 and o 2 on machine k * , then insert O * between o 1 and o 2 in the k * -th row of σ.Otherwise (no idle intervals can hold O * ), insert it at the back of the k * -th row of σ.Update the Gantt chart records for the starting time and completion time of operation O * .

Step 1 :
Perform δ simulation replications for each solution x k .Denote the mean objective value by f k , and denote the standard deviation by s k .Calculate the estimated "reward" as rk = s k /f k .Set ν k = δ (k = 1, . . ., K) and ν = δK.Step 2: Calculate a priority index for each solution x k : ρ k = rk + √ 2 ln ν ν k .Step 3: Perform δ additional simulation replications for the solution x k * with the maximum ρ value.Update f k * , s k * and rk * for this solution.Let ν k * ← ν k * + δ and ν ← ν + δ.Step 4: If ν < T , go back to Step 2. Otherwise, terminate the procedure.

Table 1 .
The computational results under normal distributions with θ = 0.1.

Table 2 .
The computational results under normal distributions with θ = 0.2.

Table 3 .
The computational results under normal distributions with θ = 0.3.

Table 4 .
The computational results under uniform distributions with θ = 0.1.

Table 5 .
The computational results under uniform distributions with θ = 0.2.

Table 6 .
The computational results under uniform distributions with θ = 0.3.

Table 7 .
The computational results under exponential distributions

Table 8 .
Mann-Whitney U tests on the computational results.