Single machine batch processing problem with release dates to minimize total completion time

Article history: Received July 18 2017 Received in Revised Format July 29 2017 Accepted August 19 2017 Available online August 2


Introduction
Batch processing (BP) problem has been attracting many academics and practitioners because of its intensive applications in the real-world industry.In BP, a resource is able to process more than one job as a batch at the same time.BP problems can be observed in different industries such as cutting machines in textile industry, testers in semiconductors manufacturing, and ovens in metalworking.A similar problem wherein jobs must be grouped into batches and contemporarily sequences of jobs and batches must be managed to satisfy a certain objective is also ascribable to the fields denoted as scheduling with tool changes (Costa et al., 2016) and group scheduling (Costa et al., 2017).
Based on the required processing time for each batch, the traditional BP problems are classified into three categories as: (1) the processing time of a batch is equal to the sum of the processing times of the jobs assigned to the batch; (2) the processing time of a batch is equal to the maximum processing time of the jobs assigned to the batch; (3) the processing time of a batch is equal to a fixed processing time to process the jobs assigned to the batch.There is also another method of classifying the BP problems based on the capacity of batches: a) the number of jobs assigned to a batch is limited by the maximum number of jobs that can be assigned to a batch; b) the number of jobs assigned to a batch depends on a size capacity for each batch based on one of the jobs attributes such as weight or volume; c) the jobs assigned to a batch must respect both conditions explained in (a) and (b).
Most of the research performed in BP problems have been focused on category (2) as for processing time of each batch and categories (a) or (b) as for the capacities of the batches.Chandru et al. (1993) propose an exact procedure and several heuristics for solving the 1|batch| j C  problem for categories (2) and (a).Uzsoy (1994) tackles both 1|batch|Cmax and 1|batch| j C  problems and proposes an exact procedure based on Branch and Bound (B&B) algorithm for the 1|batch| j C  problem.Moreover, he devises several heuristics to solve those problems for categories (2) and (b) and proves that both problems are NP-hard.Another research of Uzsoy and Yang (1997) deals with the 1|batch| j j w C  problem for categories (2) and (a).They propose several heuristics and an exact approach based on B&B algorithm.Jolai and Dupont (1997) approach the 1|rj,batch|Cmax problem for categories (2) and (b) and propose a number of heuristics for the problem.Lee (1999) studies the 1|rj,batch|Cmax problem for categories (2) and (a).He proposes several polynomial and pseudo-polynomial time algorithms for a few particular instances and develop efficient heuristics for the general problem.Liu and Yu (2000) approach the 1|rj,batch|Cmax problem for categories (2) and (a) and propose a pseudo-polynomial algorithm for the case with a constant number of release dates and a greedy heuristic for the general case.Dupont and Dhaenens-Flipo (2002) address the 1|rj,batch|Cmax problem for categories (2) and (b) and develop an exact procedure based on B&B method for the research problem.Chang and Wang (2004) develop a heuristic algorithm for the 1|rj, batch| j C  problem for categories (2) and (b).Melouk et al. (2004) and Damodaran et al. (2006) develop several meta-heuristic algorithms based on simulated annealing (SA) and genetic algorithm (GA) for the 1|batch|Cmax problem for categories (2) and (b).Damodaran et al. (2007) cope with the 1|batch|Cmax problem for categories (2) and (c) and develop a meta-heuristic algorithm based on SA to solve the problem.Parsa et al. (2010) approach the 1|batch|Cmax problem for categories (2) and (b) and propose an exact approach based on Branch and Price (B&P) algorithm.They compare their proposed algorithm with the B&B method proposed by Dupont and Dhaenens-Flipo (2002) and show that the B&P algorithm is better than the B&B algorithm.Xu et al. (2012) solve the 1|rj,batch|Cmax problem for categories (2) and (b) through a mixed integer linear programming (MILP) model.They develop an effective lower bound method, a heuristic algorithm, and an ant colony optimization (ACO) algorithm for solving the mentioned BP problem. Lee and Lee (2013) develop several heuristics for the 1|batch|Cmax problem for categories (2) and (b).Jia and Leung (2014) provide an improved max-min ant system algorithm for the 1|batch|Cmax problem for categories (2) and (b).Zhou et al. (2014) approach the 1|rj,batch|Cmax problem for categories (2) and (b) and propose various heuristics to solve the problem.Al-Salamah (2015) devises an artificial bee colony method to minimize the 1|batch|Cmax problem for categories (2) and (b).Li et al. (2015) consider the 1|batch, dj=d| ( Accordingly to the literary contributions mentioned before, there exists a large amount of research dealing with the BP problem for categories (a) and (b).Conversely, BP problems for category (c) have been weakly investigated by the body of literature so far.In this research, a single machine BP problem with release dates to minimize the total completion time (1|rj,batch| j C  ) is approached for category (c).
The processing time of each batch is as the one defined by category (2) and the capacity of batches is defined by category (c).The application of the proposed research problem can be observed in the burnin operations in the final testing step of integrated circuits in semiconductor manufacturing (Uzsoy, 1994).The burn-in operations ensure that no faulty product is accepted.In burn-in operations, the integrated circuits are placed in an oven (a batch processing machine) at a fixed temperature for a long period of time.Each circuit (job) may have a different burn-in time (processing time).The jobs are loaded onto boards and then, the boards are put into an oven.The number of boards that can process at the same time is described as oven capacity.Thus, the boards must be split into batches.The processing time of a batch is defined by the longest processing time among all jobs in that batch.The processing time in burnin operations is too lengthy compared to other testing operations.Thus, they form a bottleneck in the final stage.The minimization of the total completion time would ensure the increase of the throughput.
In this research, it is assumed that n jobs are assigned to a machine to be processed.The machine can process the jobs arranged into batches.All jobs are not available at the beginning of the planning horizon.
The machine has a size capacity and also can process at most a limited amount of jobs at the same time.
Each job has a different size.The processing time of each batch is equal to the maximum processing time of jobs assigned to the batch.The goal is to minimize the total completion time of jobs.
The rest of this paper is organized as follows.In Section 2 a mathematical model is developed for the proposed research problem.In Section 3 several meta-heuristic algorithms are proposed to heuristically solve the problem.The specifications of the test problems used to compare the performance of the proposed algorithms, also including parameter settings and solution time are explained in Section 4. In Section 5 numerical results and comparison analysis involving the proposed metaheuristic algorithms are presented.In Section 6 conclusions and future research areas are discussed.

The Mathematical Model
In this section, an original MILP model is developed for the proposed research problem.In this model, a batch is considered to be active whether it has at least one job and the maximum number of batches is equal to the number of jobs.Indexes, parameters, decision variables, and the whole mathematical model are as follows: The model The objective, as presented by Eq. ( 1), is to minimize the total completion time.Constraint (2) is incorporated into the model to ensure that each job is assigned to one batch.Constraint (3) guarantees that the total size of the jobs assigned to a batch does not exceed the batch capacity.Constraint (4) assures that the number of jobs assigned to a batch cannot be greater than the number of jobs to be assigned to each batch.Constraint (5) allows the processing time of each batch is equal to or greater than the processing time of the jobs assigned to that batch.Constraint (6) ensures that the starting processing time of a batch is greater than or equal to the arrival time of all jobs assigned to that batch.According to Constraint (7) the single-machine cannot process more than one batch at a time.Constraint (8) states the completion time of a job is equal to the completion time of the batch the job is assigned to.Constraint (9) ensures at least one job is allocated to an active batch.Constraint (10) guarantees that all active batches are ordered consecutively.Since Uzsoy (1994) proved that the 1|batch| j C  problem for categories (2) and (b) is an NP-hard problem, the single machine batch processing problem under investigation, with release dates to minimize of the total completion time (1|rj,batch| j C  ) for categories (2) and (c), can be classified as an NP-hard problem too.As a result, meta-heuristic algorithms are needed to solve large-sized issues.

Meta-heuristic algorithms
Basically, two different classes can be considered to categorize meta-heuristic algorithms: single-solution based algorithms and population based algorithms.In the first class, algorithms try to make an improvement on a single candidate solution at each iteration.Tabu search (TS), simulating annealing, variable neighborhood search are examples of algorithms belonging to this class.In the second class, algorithms try to make an improvement on multiple candidate solutions at each iteration.Genetic Algorithms (GAs), Ant Colony Optimization (ACO), particle swarm intelligence (PSO) fall under this class.In this research, due to the significant performance of TS and PSO in solving relevant BP scheduling problems (Liao and Huang (2011), Damodaran et al. (2012) and Damodaran et al. (2013)), both of them have been taken into account for solving the research problem under investigation.The following sections deal with the adopted algorithms in depth.

Tabu Search
Tabu Search (TS) was originally developed by Glover and Laguna (1999) and basically consists of a neighborhood search method that keeps track of the up-to-now search path to avoid to be trapped into local optima or to try an explorative search of the solution domain (Costa et al., 2015).The proposed TS was inspired to that of Liao and Huang (2011) but, differently from the original one, in this paper an enhanced two-level TS algorithm was devised.In the first level (inside level), the best assignment of jobs to batches is examined.In the second level (outside level), the best sequence of batches is investigated.
The relationship between the aforementioned phases is that once the inside level search is carried out to assign jobs to batches, the search process is switched to the outside level.In fact, whenever the inside level search stopping criterion is met, the best obtained job assignment is considered and the search strategy changes to the outside level.The outside search stops when the related stopping criterion is met.The best obtained solution, which combines a sequence of batches and the allocation of jobs to each batch, is the final solution.The peculiarities of the proposed TS are as follows.

Initial Feasible Solution
In this research, an effective procedure inspired to the heuristic algorithm called DYNA (Jolai & Dupont, 1997) is proposed to generate the initial feasible solution to be handled by the inside level.Then, the best solution obtained by the inside level becomes the initial feasible solution for the subsequent phase, i.e., the outside level.The algorithm generating the initial solution works as follows: Step 1: First, the jobs are ordered based on the Earliest Completion Time (ECT) sorting rule.In words, jobs are ordered in nondecreasing order of their expected completion time (r j + p j ), with ties broken by r j .
Step 2: The first job in the ECT list is assigned to the first batch.
Step 3: The subsequent job is scheduled according to the following criteria: Step 3.1: The job on the top of the list of remaining jobs is assigned separately to the existing batches having enough space.Step 3.2: If the job on the top of the list of remaining jobs cannot be allocated to any batch, a new batch is created and the job is assigned to it.
Step 4: Then the total completion time is calculated for all possible combinations.
Step 5: The state with the minimum total completion time among all possible states is selected.
Step 6: All jobs have been scheduled.Go to Step7.Otherwise go to Step3.

Neighborhood generation mechanism
After an initial feasible solution is generated, the neighborhood search at the inside level is performed by two moves, the former being a swap move and the latter being an insert move: Swap move.The job in ath position of bth batch is swapped with the job in kth position of gth batch.if bth and gth batches have enough space the move is feasible.Otherwise the move is unfeasible and it is rejected from the neighborhood.
Insertion move.The job in ath position of bth batch is removed and inserted into kth position (empty position) of gth batch.if the gth batch has enough space the move is feasible.Otherwise the move is unfeasible and it is rejected from the neighborhood.
Let us consider a BP problem with n=5 jobs and two batches.The maximum number of items per batch is N=3 and B=10 is the machine capacity.The seed solution is [1-2-3-4-5] while the size of each job s j as well as the batch i each job is allocated to, i.e., X ji , are reported in Table 1.Six distinct moves are allowed by applying the swap operator as follows: (2,1), (2,4), (2,5), (3,1), (3,2), (3,5).Move (3,1) is not feasible as the obtained batch size (6+5) exceeds the provided maximum capacity B=10.As far as the insertion operator is concerned, since N=3, the following three different moves may be generated: (2,3,1), (2,3,4), (2,3,5).Only the last move is feasible under the batch size viewpoint.As a result, just six out of nine moves can be considered as candidates to be the next seed.

Table 1
An example of seed solution As for as the outside level, the neighborhood is generated according to a regular adjacent pairwise interchange method; thus, β-1 alternative moves may be produced, where β is the current amount of batches.Hence, β-1 additional moves have to be considered for selecting the next seed solution.

The Tabu List (TL)
TL is a list of the characteristics of forbidden moves.It prevents the cycling back to formerly visited solutions by storing the characteristics of these moves for a certain period.At each iteration, the best neighbor solution, i.e., that one with the best objective function value, is selected and the move that generated that candidate solution is compared with the moves stored in the TL.Whether the move does not exist in the TL, the corresponding solution is selected as a new seed solution that, subsequently, will be subject again to the neighborhood generation mechanism.Otherwise, the move associated to the next best solution is taken into account.If a tabu move allows a better objective function value than the best global one found so far, the tabu restriction is ignored and the solution related to that move is selected as the new seed solution.The best objective function value found so far is denoted as the aspiration criterion.

Stopping Criterion
The aforementioned process is repeated until the stopping criterion based on the maximum time designated for each problem is met.

The PSO
The PSO algorithm, introduced by Kennedy and Eberhart (1995), is a population-based algorithm to solve continuous optimization problems.It starts with a set of initial solutions (particles).Each particle has its own position and velocity vector.The position vector 1 2 , ,...,  represents the position of the ith particle in the tth iteration where t ij x illustrates the jth dimension of the n-dimensional position vector.The velocity vector 1 2 , ,...,  represents the velocity of the ith particle in the tth iteration where t ij v illustrates the jth dimension of the n-dimensional velocity vector.The dimension of search space is equal to the number of jobs n in this research.The successful behavior of each particle affects the behavior of other particles, so particles move toward areas with better objective function values according to the two following targets: P best : The best position visited by the ith particle at the tth iteration is shown as P best , Since particles move towards better positions during the search process, the velocity of each particle changes on the basis of P best and G best vectors at each iteration.The velocity of each particle at each iteration is updated by: where w is the inertia weight which controls the impact of the previous velocity vector at the tth iteration while c 1 and c 2 are two constants called acceleration coefficients.Also, r 1 and r 2 are random numbers drawn from a uniform distribution U[0, 1].The position of each particle varies as velocity changes.As a result, the position of the ith particle at the (t + 1)th iteration is updated as follows:

PSO encoding
The original PSO algorithm is used to solve the continuous optimization problems.Arabameri and Salmasi (2013), and Tadayon and Salmasi (2013) apply the smallest position value (SPV) rule that is an increasing order mechanism to transform a continuous PSO to a discrete one.In other words, the SPV method allows transforming a real-encoded solution into a permutation one, that is a sequence of jobs.
For instance, assume that t i X =[2.3,-1.5,-2.03,1.16] is the position vector of the ith particle at the tth iteration whose dimension is equal to four (i.e., four jobs).Since -2.03 is the smallest value in t i X , the third member of the position vector, that is J 3, is located in the first position of the sequence.The second smallest value is -1.5, so J 2 will be the second digit of the sequence, and so on.Following the same fashion, the final permutation sequence will be: J 3 -J 2 -J 4 -J 1 .

Initial solutions
In this research, three members (particles) of the initial population are created by the following dispatching rules: shortest processing time (SPT) (Baker & Trietsch, 2009), earliest release date (ERD) (Baker and Trietsch, 2009), and ECT.The rest of the members are randomly generated.

Updating Particles
Based on Poli et al. (2007), if a multiplier χ is inserted into Eq.( 11), the convergence speed as well as the effectiveness of PSO algorithm may be improved.The appropriate value of χ is calculated through the following Equations: To satisfy Eq. ( 13), the value of 2.05 is considered for both c 1 and c 2 .Therefore, according to Eq. ( 14),  is equal to 7298.As a result, the new expression for the velocity vector of each particle is: Eq. ( 12) is regularly used to update the position of each particle.As mentioned in section 3.2.1, the SPV rule is employed to assess the quality of the solution related to each particle and to update the values of P best i and G best , if applicable.

Improving the performance of PSO
Pros and cons characterize the particle swarm optimization algorithm (Coello Coello et al., 2007).The main advantage of PSO consists of its ability in exploitation while it may give weak result in exploring the solution space.In order to boost the performance of PSO, a series of hybrid PSO algorithms, which combine PSO with other performing search techniques, have been developed by literature so far (Sha & Hsu, 2006;Xia & Wu, 2006;Chen et al., 2013;Gao et al., 2014;Javidrad & Nazari, 2017).Notably, Arabameri and Salmasi (2013) demonstrated as the performance of PSO may be significantly improved if it is combined with a proper neighborhood search.In this research, two different hybrid methods have been considered to boost the performance of the regular PSO, the former being based on GA (PSO-GA), the latter being focused on TS (PSO-TS).

PSO-GA.
The GA, just like the PSO, is a population-based stochastic algorithm that evolves through a number of initial solutions called chromosomes.At each iteration, a series of basic operators called selection, crossover and mutation are performed on the population with the aim of performing both exploration and exploitation (Michalewicz, 1996;Lin & Kang, 1999).In this research, the selection of offspring is randomly executed on the population in order to strengthen the diversification phase.An arithmetic crossover operator based on Simon (2013) was employed to produce two chromosomes called offspring by combining two parent chromosomes.In this method, the genes related to two randomly selected individuals (ith and jth particles) are exchanged and modified through Eq. ( 16) and Eq. ( 17), where α l is an l-dimensional vector of random weights between 0 and 1. 1 kl X and 2 kl X represent the offspring particles at the kth iteration.The part of population on which the crossover operator is executed at each iteration depends on the crossover rate p c .

 
The goal of mutation is to maintain the population diversity, thus avoiding any premature convergence.It makes a new offspring up from a single parent and assists the algorithm to avoid to be trapped in local optima.A single parent (ith particle) is selected randomly to perform a Gaussian mutation procedure based on Simon (2013).Then, the jth element of the selected particle (X i ) is varied according to Kuo and Han (2011), as follows: where N(0,1) is a number from a standard normal distribution and rand[0,1] is a random number form a uniform distribution in the range [0,1].The portion of population in which mutation procedure is applied to depends on the mutation rate (p c ).A pseudo-code of the proposed PSO-GA is reported in the following paragraph.

Pseudo-code of PSO-GA
Step 1. Initialize control parameters and create a swarm with P particles.Do While (stopping criterion is not met) Step 2. Update position and velocity vectors for i= 1 to P for j=1 to n Step 3. Arithmetic crossover operator for k=1 to nc Select two particles randomly (ith and jth particles)

end for end for
Step 4. Gaussian mutation operator for q=1 to nm Select one single parent randomly (ith) and change the jth element of it Step 5. New particles Merge all the newly generated particles yielded by the crossover, the mutation, and PSO operators.Then, select the best P particles in terms of objective function value.
Step 6. Update the Pbest and Gbest vectors for Gbest = Pbests end if end for update inertia weight loop PSO-TS.Actually, three different versions of PSO-TS, hereinafter denoted with the subscripts a, b and c, have been developed as explained below: a) In the first version, whenever the G best is updated at a certain iteration, it is given to TS as a seed solution and TS is performed for a finite amount of time.If a better solution is found, it is considered as the new G best .

b)
In this version, whenever the P best is updated at a certain iteration, it is given to TS as a seed solution and TS is performed for a finite amount of time.If a better solution is found, it is considered as the new P best .c) In the last scenario, PSO and TS are hierarchically applied to the problem at hand.The regular PSO is executed in the first phase until a switching criterion is satisfied.Subsequently, the best solution found by the PSO is handled by the TS algorithm until a time-based stopping criterion is met.
It is worthy to point out that the aforementioned versions of hybrid PSO are powered by the same twolevels tabu search described in Section 3.1.

Calculating the objective function
Based on empirical studies, an appropriate assignment of jobs to the batches has a considerable effect on the value of the objective function.Therefore, the proposed heuristic algorithm as in Section 3.1.1was employed to calculate the total completion time of each solution.

Computational experiments
In order to compare the performances of the proposed meta-heuristic algorithms, a set of test problems have been generated, randomly.The whole set of test problems are solved by the five proposed algorithms, namely TS, PSO-GA, PSO a , PSO b , and PSO c , and subsequently the performances of these algorithms have been compared.The test problems specifications, parameter setting, the adopted solution time as well as the obtained numerical results are dealt with the subsequent sub-sections.

Test problems
To evaluate the performance of the proposed algorithms, a wide range of test problems has randomly been generated on the basis of the following five different factors: the number of jobs (n), the maximum number of jobs in a batch (N), the maximum capacity of the machine (B), the size of jobs (s j ) and the processing time of jobs (p j ).Three different classes of problems (small, medium, large) have been provided for each factor, as depicted in Table 2. Hence, the total amount of scenario problems to be investigated is 3 5 = 243.
Two replicates have been considered for each class; thus, a total amount of 486 test problems have been solved by means of each algorithm.As far as the release dates are concerned, they are drawn from the uniform discrete distribution U[ 0,( / ) max( ) ] and then rounded to the closer integer value; max(p j ) is the maximum value of processing time among the n jobs.

Setting of control parameters
Both effectiveness and efficiency of metaheuristic algorithms can be enhanced by giving appropriate values to their control parameters.As for TS algorithm the most influencing parameter to be chosen is the TL size.In this paper an empirical formula for choosing the TL size has been conceived making use of a trial-and-error approach involving an extended number of instances.Table 3 shows how the TL size should be set for both the first and the second level, on the basis of the expected neighborhood size.In fact, especially for the first level of TS, the neighborhood size (NS1) is not a priori known as it depends on how the jobs are allocated to the batches, conforming to the other parameters, i.e., N, B, s j .On the other hand, as concerns the second level of TS, the neighborhood size (NS2) is equal to β-1, where β is the current number of batches.To calculate both NS1 and NS2 the heuristic described in Section 3.3.1 has been employed.Whether the TL size in Table 3 does not yield any integer number, the obtained result must be properly rounded down.The control parameters pertaining to the other proposed algorithms (i.e., PSO a , PSO b , PSO c , PSO-GA) have been tuned through the response surface methodology (RSM) method.RSM is a collection of statistical and mathematical techniques used for developing, improving, and optimizing processes.In RSM, there are various input variables (control parameters) that can potentially influence some performance measures or quality characteristics that are called response (the objective function of the algorithm).Usually, in the RSM approach, it is convenient to transform the original variables into coded variables x 1 ,x 2 ,...,x l , which are usually defined to be dimensionless with mean zero and the same spread or standard deviation (Raymond et al., 2009).Eq. (190 shows how to transform an original variable to a coded one, where X i and x i represent the actual variable and the coded variable, respectively.Thus, the response function can be written as in Eq. ( 20).
where l represents the number of input variables.The goal is to gain the most suitable level of the algorithm parameters to optimize the response value.The main hypothesis is that the independent variables are continuous and controllable by experiments with negligible errors.Finding a proper approximation for the true functional relationship between independent variables and the response surface is required (Kwak, 2005).Raymond et al. (2009) proposed a second-order model, as follows: where y is the predicted response, β 0 is the model constant, β j is the linear coefficient, β jj is the quadratic coefficient, and β ij is the interaction coefficient.Shi and Eberhart (1999) stated that a large value of inertia weight w simplifies probing new positions (global search), while a small value of inertia weight simplifies a local search.As a result, a suitable balancing between local and global search can be achieved by adaptively reducing the inertia weight linearly according to Eq. ( 22), where t indicates the time of the current iteration and T refers to the total time considered for each test problem, respectively.Also, w start and w end are the starting and the finishing values for the inertia weight, respectively.
  Basically, NP, w start , and w end , p c , and p m are the control parameters considered as the input variables for PSO-GA.The input variables of the other PSO-based algorithms (namely PSO a , PSO b , and PSO c) are NP, w start , w end and TL size.According to Eq. ( 19), the PSO-GA parameters are denoted by X 1 (NP), X 2 (w start ), X 3 (w end ), X 4 (p c ), X 5 (p m ), while the other control parameters are denoted by X 1 (NP), X 2 (w start ), X 3 (w end ), X 4 (TL size).Each control parameter has been varied at three levels, conforming to low (-1), center (0) and high level (+1), as depicted in Table 4 and Table 5.As concerns the NP parameter, three scenarios depending on the problem size have been taken into account.The Box-Behnken experimental design was employed to handle a family of efficient three-level designs fitting the second-order response surfaces (Raymond et al., 2009).The number of experiments is equal to 2l(l-1)+n c , where n c is the number of the central points.As a result, 48 and 32 experiments have been performed for tuning PSO-GA and the different PSO-based algorithms, respectively.

Computational time
The exit criterion of the proposed metaheuristic algorithms is based on the computational time (CT) and, conforming to Gohari and Salmasi (2015), it has been parametrized as a function of the number of jobs (n), as follows: , where CT max and n up respectively are the maximum allowed computational time and the upper value of the range related to the number of jobs; as for example, in the range [5,20], n up is equal to 20.CT max has been set to 30, 90, and 180 seconds for small-, medium-, and large-sized problems, respectively.Hence, the values of δ for small-, medium-, and large-sized problems are equal to 1.5, 1.8 and 1.8, respectively.All the proposed algorithms were coded in C++ programming language.The whole set of experiments was performed on a PC with 2.8 GHz CPU and 2 GB RAM.

Numerical results and comparison analysis
In order to compare the different metaheuristics, a Relative Percentage Deviation (RPD) performance indicator as in Eq. ( 24) has been taken into account.It is worth pointing out that, for most small-size problems, the RPDs have been computed on the basis of the global optima obtained by solving the mathematical problem.As the reader can notice, in the last two rows both grand averages (ave) and standard deviations (st.dev) highlight the effectiveness of PSOc and PSO-GA over the other metaheuristics.In addition, the small values of RPDs confirm that each algorithm has been properly calibrated and designed for the problem at hand, though TS is strongly less performing than the other competitors.Table 11 reports the RPD medians for each algorithm, with respect to the three different classes of problems.PSO-GA and PSOc achieve the best results, regardless of the specific class of problem.Notably, for small and medium sized issues, medians related to PSO-GA and PSOc are comparably close to zero.Conversely, a slight advantage for PSO-GA emerges for solving the case of larger sized issues.The weakness of TS comes to light again, especially when the size of the problems increases.Similarly, the performances of PSOa and PSOb drastically deteriorate with the problem size though they are able to assure a median equal to zero for the smaller sized issues.

Conclusions and Future Research
In this research, we have approached a single machine batch processing problem with release dates to minimize the total completion time (1|r j , batch| j C  ).A mathematical model has been proposed to solve the research problem optimally.The proposed research problem is known to be NP-hard; hence, several meta-heuristic algorithms based on TS and PSO with two different approaches have been developed to heuristically solve the problem.Since the normal probability plot of residuals were not normally distributed, a non-parametric test (Kruskal-Wallis) was employed to compare the performance of the proposed algorithms.The results show that there was a significant statistical difference among the metaheuristic algorithms.Medians and Z rank demonstrate that PSO a and PSO-GA assure the most promising performances compared to the other metaheuristics.Therefore, in order to find the algorithm with the best efficiency between PSO c and PSO-GA, a post-hoc pairwise comparison based on a non-parametric Mann-Whitney U test was employed.The result of Mann-Whitney test indicates that PSO-GA had a better performance than PSO c .Since the proposed research problem was investigated in a single machine environment, the obtained results in this research can be applied in other environments such as parallel machines and flow shop.A lower bounding method can be provided to evaluate the performance of the proposed meta-heuristic algorithms for the future research.


(2) and (b) and develop a hybrid GA by combining GA with a heuristic algorithm managing the batches.Parsa et al. (2016) introduce a hybrid meta-heuristic algorithm based on the maxmin ant system for 1|batch| j C for categories for (2) and (b).

Fig. 1 .Fig. 2 .Fig. 3 .
Fig. 1. 3D surface plot involving NP and pc for PSO-GA Fig. 2. 3D surface plot involving NP and pm for PSO-GA The best position visited by all particles at the tth iteration is called G best ,

Table 2
Settings about the test problems

Table 3
The different values of TL size for different number of jobs

Table 4
The levels of the values of each parameter for PSO-GA

Table 5
The levels of the values of each parameter for PSO , PSO , and PSO

Table 6
Tuned value of PSO-GA's parameters

Table 8
Tuned value of PSO b 's parameters

Table 9
Tuned value of PSO c 's parameters Whereas, for both medium and larger-sized issues, each RPD was computed by exploiting the relative local optimum value among the different competing algorithms.As for the smaller-sized class of problems, the global optima have been achieved by ILOG CPLEX (version 12.2).Actually, the MILP model was able to optimally solve only 74 smaller-sized test problems out of 162; thus, the RPD values of the remaining 88 test problems have been computed making use of the local optima instead of the global ones.Table10shows the average values of RPDs over the two replicates per each scenario problem, related to the smaller-sized issues.

Table 10
Small-sized problems: average RPD results over the two replicates