An improved meta-heuristic approach for solving identical parallel processor scheduling problem

Abstract

This article considers the problem of scheduling n jobs on m identical parallel processors where an optimal schedule is defined as one that produces minimum makespan (the completion time of the last job) and total tardiness among the set of schedules. Such a problem is known as identical parallel processor makespan and total tardiness problem. In order to minimize makespan and total tardiness of identical parallel processors, improved versions of particle swarm optimization and harmony search algorithm are proposed to enhance scheduling performance with less computational burden. The major drawback of particle swarm optimization in terms of premature convergence at initial stage of iterations is avoided through the use of mutation, a commonly used operator in genetic algorithm, by introducing diversity in the solution. The proposed algorithm is termed as particle swarm optimization with mutation. The convergence rate of harmony search algorithm is improved by fine tuning of parameters such as pitch adjusting rate and bandwidth for improving the solution. The performance of the schedules is evaluated in terms of makespan and total tardiness. The results are analyzed in terms of percentage deviation of the solution from the lower bound on makespan. The results indicate that particle swarm optimization with mutation produces better solutions when compared with genetic algorithm and particle swarm optimization in terms of average percentage deviation. However, harmony search algorithm outperforms genetic algorithm, particle swarm optimization, and particle swarm optimization with mutation in terms of average percentage deviation. In certain instances, the solution obtained by harmony search algorithm outperforms existing clonal selection particle swarm optimization.

Introduction

Scheduling is a decision-making process primarily concerned with allocation of operations on machines in such a manner that some process performance measures such as makespan, tardiness, earliness, and flow time can be minimized.^1,2 The scheduling procedures available today are not only quite realistic but also possess fast solution-generating capability. It finds widespread applications in both manufacturing and service industries. ¹ The parallel processor scheduling is one of the scheduling scenarios that arise when a set of n jobs available at time 0 needs to be processed on m identical parallel processors (IPPs; machines). Each job has to be processed without interruption on one of the available m machines with processing time p_ij where jth job is processed on machine i. Each machine can process only one job at a time. Furthermore, the jobs cannot be processed more than one machine at the same time. The objective is to find a job schedule that can minimize the makespan and total tardiness. The problem of scheduling n jobs on m IPPs where an optimal schedule is defined as one that gives the smallest makespan (the completion time of the last job) and minimum total tardiness among the set of schedules is known as IPP makespan and total tardiness problems (IPPMST). Using the standard three field notations, makespan is normally denoted as $P\parallel C_{max}$ and total tardiness as $P\parallel \sum T_j$ , where P designates IPPs and the jobs that are not constrained, C_max represents makespan, and ΣT_j denotes total tardiness. Such a scheduling scenario is treated as NP-hard problem, and the exact solution for small size problems has been reported using traditional operation research methods such as dynamic programming and mixed-integer linear programming.^3–8 However, the computational time and cost increase in a non-polynomial manner for a reasonable size problem. In practical situation, large-scale problems are dealt with a view to obtain quick solutions at less computational effort. As exact models cannot address this issue, a large number of heuristics have been proposed for solving IPPMST.^9–15 Although heuristics generate good solutions, the solution quality is far away from optimal solution. Therefore, nature-inspired algorithms like tabu search (TS), simulated annealing (SA), genetic algorithm (GA), and particle swarm optimization (PSO) are proposed recently to generate approximate solutions close to the optimum with considerably less computational time.^13,15–17 Among these algorithms, PSO is considered as a robust algorithm for its computational efficiency than other algorithms, and also it needs very few parameters for tuning the algorithm. In PSO, the whole algorithm rests on two equations, namely, velocity and position updation rules. Normally, swarm-based optimization algorithms trap at local optima where the algorithm loses its diversity after few iterations. To avoid the premature convergence, mutation operator has been borrowed from GA to be embedded with PSO to avoid being trapped in local optima. Furthermore, this algorithm has been compared with the latest search-based algorithm known as harmony search algorithm (HSA).

In this article, smallest position value (SPV) rule proposed by Tasgetiren et al., ¹⁸ originally borrowed from the random key representation of Bean, ¹⁹ is used for problem representation in PSO for sequencing the jobs in IPPMST. Moreover, a non-uniform mutation operator is incorporated when it is observed that search procedure starts losing diversity and tends toward premature convergence. ²⁰ The results obtained from particle swarm optimization with mutation (MPSO) have further been validated using HSA with SPV rule problem representation integrating Mahdavi et al.’s ²¹ fine tuning schemes.

Literature review

In order to solve identical parallel machine scheduling problem, Lawler ² has proposed a transportation problem model when the job processing times are identical. Dogramaci ³ has proposed a dynamic programming algorithm to minimize total weighted tardiness. Elmaghraby and Park ⁴ have used a branch and bound (B&B) technique to minimize some penalty functions of tardiness considering identical due dates and processing times. Later, this method is modified by Barnes and Brennan, ⁵ Azigzoglu and Kirca, ⁶ and Yalaoui and Chu ⁷ to reduce the computational burden. Tanaka and Araki ⁸ have proposed a solution strategy based on B&B algorithm with Lagrangian relaxation technique to obtain a tight lower bound (LB). Most of the exact methods are limited to solve generalized IPPMST because they consider either equal processing times or common due dates or both. Therefore, a number of heuristic methods have been proposed to solve the problem. To minimize mean tardiness on a single machine, Wilkerson and Irwin ²² have presented a heuristic method based on neighborhood search. This method has been extended to solve the parallel machine case by Baker and Scudder. ²³ Kim et al. ²⁴ have considered the problem of determining the allocation and sequence of jobs on parallel machines for the objective of minimizing total tardiness. Kumar et al. ²⁵ have proposed a SA-fuzzy logic approach to select the optimal weighted earliness–tardiness combinations in a non-identical parallel machine environment. Dogramaci and Surkis ²⁶ have implemented a heuristic method based on list algorithms designed for parallel machine problem. Ho and Chang ²⁷ have proposed a new measurement technique called traffic congestion ratio combining processing time and due dates to obtain traffic priority index for job assignment. Simon and Andrew ²⁸ and Patrizia et al. ²⁹ have compared various heuristics on similar instances to study the effectiveness of heuristics. With the development of evolutionary computing techniques, meta-heuristics procedures are used to solve IPPMST problems to generate approximate solutions close to the optimum with considerably less computational effort. GA has been applied by Bean ¹⁹ to solve IPPMST where a representation technique called random keys is proposed to maintain feasibility from parent to offspring in GA. With the objective of workflow balancing, Rajakumar et al. ³⁰ used GA to solve parallel machine scheduling problem. Chaudhry and Drake ³¹ have presented a GA approach to minimize total tardiness of a set of tasks for identical parallel machines and worker assignment to machines. Koulamas ³² has proposed a polynomial decomposition heuristic and a hybrid SA heuristic for the problem. Kim and Shin ³³ have implemented TS algorithm that schedules jobs on parallel machines when release times and due dates of jobs are known a priori and at sequence-dependent setup times. Bilge et al. ³⁴ have developed a TS approach to solve the IPPMST with sequence-dependent setups on uniform parallel machines. Anghinolfi and Paolucci ³⁵ have proposed a hybrid meta-heuristic approach that integrates several features of TS, SA, and variable neighborhood search (VNS) for solving IPPMST.

For minimizing the makespan, Graham ⁹ has used well-known longest processing time (LPT) rule for assigning the list of jobs to the least loaded machines. Coffman et al. ¹⁰ have proposed the MULTIFIT heuristic and proved the close relation between the bin-packing problem and maximum completion time problem. Lee and Massey ¹¹ have proposed COMBINE heuristic that utilizes the LPT rule for an initial solution and improves the schedule by MULTIFIT algorithm. Gupta and Ruiz Torres ¹² have proposed LISTFIT heuristic for minimizing makespan in identical parallel machines with new LB). Kurz and Askin ³⁶ have investigated a hybrid flexible flow line environment with identical parallel machines and non-anticipatory sequence-dependent setup times with an objective to minimize the makespan.

Min and Cheng ¹³ have proposed GA for solving large-scale identical parallel machine scheduling problem with too many jobs and machines to minimize makespan. Kennedy and Eberhart ³⁷ have proposed PSO algorithm having the capability of updating of position and velocity of particles and quick convergence to obtain near-optimal solutions. Tasgetiren et al. ¹⁸ have proposed a PSO for permutation flow shop scheduling problem. Liao et al. ³⁸ and Liu et al. ³⁹ have also used PSO in flow shop scheduling problem. Niu et al. ⁴⁰ have proposed a clonal selection algorithm to improve the swarm diversity and avoid premature convergence. Geem et al. ⁴¹ have proposed HSA, a nature-inspired algorithm, mimicking the improvisation of music players. The harmony in music is analogous to the optimization solution vector, and the musician’s improvisations are analogous to the local and global search schemes in optimization techniques. It has been successfully applied to various optimization problems such as sequential quadratic programming, structural optimization, and so on.^21,42–44 Zhang et al. ⁴⁵ have addressed the dynamic job shop scheduling problem with random job arrivals and machine breakdowns. Renna ⁴⁶ has adopted a pheromone-based approach in creating schedules for the job shop environment in cellular manufacturing systems. Rui et al. ⁴⁷ have developed bidirectional convergence ant colony algorithm to solve the integrated job shop scheduling problem with tool flow in flexible manufacturing system (FMS). Kim and Lee ⁴⁸ have considered the scheduling problem in hybrid flow shops with parallel machines at each serial production stage where each job may visit each stage several times, called the re-entrant flows. Marimuthu et al. ⁴⁹ have addressed the problem of making sequencing and scheduling decisions for n jobs in m machine flow shops with a lot-sizing constraint. Pan et al. ⁵⁰ have presented a novel multiobjective particle swarm optimization (MOPSO) algorithm for solving no-wait flow shop scheduling problems with makespan and maximum tardiness criteria. Rossi and Dini ⁵¹ have presented a GA for generalized job shop problem-solving, and the generalization includes feeding times, sequences of setup-dependent operations, and jobs with different routings among work centers including “multi-identical” machines. Wang and Brunn ⁵² have proposed an effective GA for job shop sequencing and scheduling where a simple heuristic rule is adapted and embedded into the GA to avoid the production of unfeasible solutions.

Scheduling on IPPs

The scheduling of IPPs is regarded as one of the most-studied classical scheduling problems that involve the assignment of multiple jobs onto the machines. Among the performance measures considered, makespan, total tardiness, and earliness are commonly used to test the schedules. The problem can be stated as follows: there are n jobs waiting for being scheduled on m IPPs so as to minimize makespan and total tardiness. The basic assumptions need to be considered for solving IPPMST problems are as follows:

1. All jobs can be processed on any of the parallel machines.

2. Each of the parallel machines can process at most one job at a time.

3. Each job has to be processed without interruption on one of the machines.

4. All jobs are available to be processed at time 0.

5. No job can be processed by more than one machine.

6. The processing time and the due date of the jobs, the number of jobs, and the number of machines are given and fixed.

7. No downtime is considered, and setup time is included in processing time.

A set N = {J₁, J₂, …, J_n} of n jobs are to be scheduled on a set M = {M₁, M₂, …, M_m} of m identical parallel machines. The first objective is to find the optimal schedule S = {S₁, S₂, …, S_m} where S_j is a subset of jobs assigned to machine M_j such that max{C₁(S), C₂(S), …, C_m(S)} = C_max(S) is minimum where

$$C_j\left(S\right)=\sum _{p_i\in S_j}{p}_i$$

(1)

The following notations have been used for finding total tardiness: i is the index of machine, i = 1, 2, …, m, where m is the number of machines; $\Pi =\{S_1,S_2,\dots ,S_j,\dots ,S_{k_i}\}$ is the sequence of jobs which are assigned to machine i; S_j is the index of job, where j = 1, 2, …, k_i and k_i represents the number of total jobs to be processed by machine i; $p_{S_j}$ is the processing time of job S_j; $S{T}_{S_j}$ is the starting time of job S_j; $c_{S_j}$ is the completion time of job S_j; and $d_{S_j}$ is the due date of job S_j.

The IPP total tardiness can be represented mathematically in the following manner

$$S{T}_{S_1}=0$$

(2)

Constraint (2) describes the starting time of the machine is equal to 0

$$c_{S_1}=S{T}_{S_1}+p_{S_1}$$

(3)

Constraint (3) denotes the completion time of the job in the sequence which is the summation of starting time and the processing time of the particular job in the sequence

$$S{T}_{S_j}=c_{S_{j-1}},j>1$$

(4)

$$c_{S_j}=S{T}_{S_j}+p_{S_j}$$

(5)

The tardiness time of job S_j is calculated using the following equation

$$T_{S_j}=max(0,c_{S_j}-d_{S_j})$$

(6)

The second objective of this article is to minimize total tardiness which is defined using equation (7)

$$T=\sum _{i=1}^m\sum _{j=1}^{k_i}{T}_{S_j}$$

(7)

The mixed-integer programming approach for IPP machine scheduling problem formulated above can be solved by B&B technique, but computational efforts are prohibitive for large-scale problems.

PSO

PSO is inspired by the social behavior of bird flocking and fish schooling. It is a population-based, bio-inspired optimization algorithm originally introduced by Kennedy and Eberhart. ³⁷ In PSO, particles are called as potential solutions; each particle adjusts its position based on its own experience as well as the experience of a neighboring particle by utilizing the best known position in the search space as well as the entire swarm’s best known position. By adapting the global and local exploration capabilities, PSO has maintained its flexibility toward reaching the global optimal solutions. The particles in PSO have their own memory that enables to retain good solutions among all particles. But in case of GA, the previous knowledge of the problem has been destroyed once the population changes and chromosomes share their information with each other. In PSO, previous best value is called as P_best of the particle and the best value of all the particles among the P_best in the swarm is called as G_best. Compared with GA, all the particles in PSO tend to converge to the best solution quickly even in the local version in most cases, and it does not require any kind of classical optimization methods such as gradient descent and quasi-Newton to solve optimization problem. Due to its simplicity, easy implementation, and quick convergence, PSO has gained much attention among the researchers and has been successfully applied to a wide range of applications in almost all fields of engineering.^18,20

The basic elements of PSO are summarized as follows: The particle $X_i^t$ denotes the ith particle in the swarm at iteration t and is represented by D number of dimensions as $X_i^t=[x_{i1}^t,x_{i2}^t,\dots ,x_{\mathrm{iD}}^t]$ where $x_{\mathrm{ij}}^t$ is the position value of the ith particle with respect to the jth dimension (j = 1, 2, …, D). The population pop^t is the set of ρ particles in the swarm at iteration t, $\mathrm{po}{p}^t=[X_1^t,X_2^t,\dots ,X_{\rho }^t]$ . The permutation of the jobs is represented by ${\pi }_i^t=[{\pi }_{i1}^t,{\pi }_{i2}^t,\dots ,{\pi }_{\mathrm{in}}^t]$ , where ${\pi }_{\mathrm{ij}}^t$ is the assignment of job j of the particle i in the permutation at iteration t. The velocity of the particle is represented as $V_i^t=[v_{i1}^t,v_{i2}^t,\dots ,v_{\mathrm{iD}}^t]$ , where $V_{\mathrm{ij}}^t$ is the velocity of the particle i at iteration t with respect to the jth dimension. The inertia weight w^t is a parameter to control the impact of the previous velocities on current velocity. For each particle, the personal best is defined as $P_i^t=[p_{i1}^t,p_{i2}^t,\dots ,p_{\mathrm{iD}}^t]$ , where $p_{\mathrm{ij}}^t$ is the position value of the ith personal best with respect to the jth dimension. The global best is defined as $G^t=[g_1^t,g_2^t,\dots ,g_D^t]$ , where $g_{\mathrm{ij}}^t$ is the position value of the ith global best with respect to the jth dimension. The velocity and position of the particles are updated in equations (8) and (9)

$$v_{\mathrm{ij}}^t=w^{t-1}{v}_{\mathrm{ij}}^{t-1}+c_1{r}_1(p_{\mathrm{ij}}^{t-1}-x_{\mathrm{ij}}^{t-1})+c_2{r}_2(g_j^{t-1}-x_{\mathrm{ij}}^{t-1})$$

(8)

$$x_{\mathrm{ij}}^t=x_{\mathrm{ij}}^{t-1}+v_{\mathrm{ij}}^t$$

(9)

where c₁ and c₂ are social and cognitive parameters and r₁ and r₂ are uniform random numbers between (0, 1)

$$w={w_t}_{-1}\times \beta$$

(10)

where β is the decrement factor. The updation of equations (8) and (9) may lead the particles moving toward compound vector of global and local optimal solutions.

Proposed PSO algorithm

Solution representation

The solution representation is one of the important issues when designing PSO. To maintain the bonding between the problem domain and PSO particles for IPPMST, D number of dimensions is represented for n number of jobs. The particle $X_i^t=[x_{i1}^t,x_{i2}^t,\dots ,x_{\mathrm{iD}}^t]$ corresponds to the continuous position values for n number of jobs in IPPMST. The particle does not have the capability of doing permutation by itself. Hence, Tasgetiren et al. ¹⁸ introduced SPV rule to determine the permutation implied by the position values $x_{\mathrm{ij}}^t$ of particle $X_i^t$ as shown in Table 1.

Table 1.

SPV rule.

Dimension, j	1	2	3	4	5	6
$x_{\mathrm{ij}}^t$	2.13	1.21	−1.20	−2.92	0.97	−1.91
$v_{\mathrm{ij}}^t$	1.78	−0.34	−1.8	1.54	−0.81	0.6
Jobs, ${\Pi }_{\mathrm{ij}}^t$	4	6	3	5	2	1

SPV: smallest position value. The bold values shows the smallest position value.

The smallest position in the table is $x_{i4}^t=-2.92$ . So, the dimension j = 4 is assigned to the first job ${\Pi }_{i1}^t=4$ in the permutation ${\Pi }_i^t$ . Similarly, next smallest values are arranged in sequence. The different permutations at each iteration t can be obtained by updating the position of each particle. The initial sequence of jobs is ${\Pi }_{\mathrm{ij}}^t=[4,6,3,5,2,1]$ . After calculating the sequence of jobs by SPV rule, the jobs have been assigned to the machines with the same sequence. Let p_j be the processing time of n number of jobs as shown in Table 2.

Table 2.

Processing times for the jobs.

Jobs j	1	2	3	4	5	6
pj	6	8	10	9	12	7

For example, there are two machines (m = 2) and six jobs (n = 6). From the SPV rule, the fourth job occupies the first position in the initial sequence. The corresponding processing time of the particular job has been taken from Table 2 and placed in the first machine as shown in Table 3. Similarly, the sixth job in the initial sequence occupies the second position. So, the corresponding processing time for the sixth job from Table 3 has been placed in the second machine. Both the machines have been assigned with two different jobs. Among these machines, the machine that completes the assigned job first is assigned with a new job from the sequence in Table 1. In this case, second machine assigned with job number 6 completes the job first as processing time for the job is 7 units of time. The third position in the sequence is job number 3 having the processing time of 10 units of time. Hence, the second machine is occupied by the job number 3. This process is repeated until all the jobs are assigned to machines. From Table 3, it is evident that job completion time of first machine is 27 units of time and that of second machine is 25 units of time. Hence, the makespan has been considered as the maximum of the completion time in all the machines. In this case, makespan is considered as 27 units of time. From Table 3, p_ij denotes the processing time of the jth jobs in ith machine where 9(4) represents fourth job is processed having 9 units of time.

Table 3.

Assignment of jobs using SPV rule.

Machines mi	1	2
pij	9(4)	7(6)
	12(5)	10(3)
	6(1)	8(2)

SPV: smallest position value.

Mutation operator

Mutation operator adopted from GA tends to alter one or more gene values in a chromosome from the initial state. This leads to the generation of new gene value that is being added to the gene pool. Moreover, GA can be able to attain better solution than the previous one through these newly generated genes. Hence, mutation is an important operator of the genetic search as it helps to avoid local optimal by preventing the population of chromosomes. In case of PSO, the lack of population diversity among the swarms is considered as a factor for the convergence on local minima. Hence, the PSO has been incorporated with mutation operator to enhance its global search capacity and to improve performance. In this article, a non-uniform mutation operator by Michalewicz ²⁰ has been incorporated in PSO which helps to mutate some particles selected randomly from the swarm. If the number of iterations with no diversity in the solution exceeds, mutation is carried out. This operator works by changing a particle position dimension ( $x_{\mathrm{iD}}$ )

$$\mathrm{Mutate}\left(x_{\mathrm{iD}}\right)=\{\begin{array}{c}\hfill x_{\mathrm{iD}}+\mathrm{delta}(t,U-x_{\mathrm{iD}}):\mathrm{rb}=1\hfill \\ \multicolumn{1}{c}{x_{\mathrm{iD}}+\mathrm{delta}(t,x_{\mathrm{iD}}-L):\mathrm{rb}=0}\end{array}$$

(11)

where t is the current iteration number, U is the upper bound value of the particle dimension, L is the LB value of the particle dimension, rb is the randomly generated bit, and delta (t, y) returns the value in the range [0:y]

$$\mathrm{delta}(t,y)=y(1-r^{{(1-t/T)}^b})$$

(12)

where r is the random number generated from the uniform distribution range [0:1], T is the maximum number of iteration, and b is the tunable parameter which is set as 5 by Michalewicz. ¹⁹ The pseudo-code for the proposed MPSO is given as follows:

• Initialize the parameters, including swarm size, maximum number of iteration, w, c₁, c₂

• While (termination condition i.e. maximum Iteration)

• Do

• t = 0

• Initialize particle’s position and velocity stochastically;

• Apply SPV rule;

• Evaluate each particle’s fitness using the objective function;

• Initialize p_best position;

• Initialize g_best position with the particle with lowest fitness in the swarm;

• t = t + 1;

• Update the velocity of the particles by equation (8);

• Update the position of the particles by equation (9);

• Apply SPV rule;

• Evaluate each particle’s fitness using the objective function;

• Perform mutation if t < (T_max* PMUT) where T_max is maximum number of iteration and PMUT is probability of mutation;

• Find the new g_best and p_best values by comparison;

• Update g_best of the swarm and p_best of each particle;

• end do

• end

The pseudo-code has been further explained with the introduction of flow chart as shown in Figure 1.

Figure 1.

Flow chart for MPSO algorithm.

HSA

Geem et al. ⁴¹ introduced an interesting meta-heuristic algorithm called HSA. It is population based and mimics the musical process of searching for a better state of harmony such as jazz improvisation. It has its own advantages over traditional optimization techniques by having great power in global search and its simplicity in both concept and implementation. Due to its advantages and novelty, HSA has received increasing attention among the researchers and has been vigorously applied to various engineering optimization problems. One of the finest qualities of HSA is the identification of high-performance regions of the solution space at a reasonable time. However, while performing the local search for numerical application, it is not able to withstand at its best. Hence, fine tuning the HSA characteristics becomes inevitable. In order to improve such fine tuning characteristics and convergence rate, Mahdavi et al. ²¹ have proposed an improved harmony search which has the fine tuning feature of mathematical techniques having the potential of outperforming traditional HSA. This algorithm uses harmony memory considering rate (HMCR) and pitch adjusting rate (PAR) for finding the solution vector in the search space. The optimization procedure of HSA is given in the following steps.

1. Initialization of the harmony memory (HM);

2. Improvisation of a new HM;

3. Updation of the HM;

4. Checking for the stopping criteria. Otherwise repeat steps 2 and 3.

Initialization of HM

A set of initial solutions of harmony memory size (HMS) is generated to build the HM, which is represented by a matrix of two dimensions where rows consists of a set of solutions x_i (population size), while columns consists of the variables of each solution (jobs). Each solution xⁱ is considered as a one-dimensional array. The size of the array can be designed by the maximum number of jobs considered in the problem instance. HM is represented as given in equation (13)

$$\mathrm{HM}=\left[\begin{array}{ccccc}\hfill x_1^1\hfill & \hfill x_2^1\hfill & \hfill \dots \hfill & \hfill x_{K-1}^1\hfill & \hfill x_K^1\hfill \\ \multicolumn{1}{c}{x_1^2}& \hfill x_2^2\hfill & \hfill \dots \hfill & \hfill x_{K-1}^2\hfill & \hfill x_K^1\hfill \\ \multicolumn{1}{c}{.}& \hfill .\hfill & \hfill .\hfill & \hfill .\hfill & \hfill .\hfill \\ \multicolumn{1}{c}{.}& \hfill .\hfill & \hfill .\hfill & \hfill .\hfill & \hfill .\hfill \\ \multicolumn{1}{c}{.}& \hfill .\hfill & \hfill .\hfill & \hfill .\hfill & \hfill .\hfill \\ \multicolumn{1}{c}{x_1^{\mathrm{HMS}-1}}& \hfill x_2^{\mathrm{HMS}-1}\hfill & \hfill \dots \hfill & \hfill x_{K-1}^{\mathrm{HMS}-1}\hfill & \hfill x_K^{\mathrm{HMS}-1}\hfill \\ \multicolumn{1}{c}{x_1^{\mathrm{HMS}}}& \hfill x_2^{\mathrm{HMS}}\hfill & \hfill \dots \hfill & \hfill x_{K-1}^{\mathrm{HMS}}\hfill & \hfill x_K^{\mathrm{HMS}}\hfill \end{array}\right]$$

(13)

Improvisation of new HM

The divergence and convergence of the search in HSA has been maintained during this process. HMCR and PAR are considered as the main parameters for convergence or divergence of the search capability by Mahdavi et al. ²¹ Here, the new solutions are constructed stochastically using one of the following three operators: (1) memory consideration (based on the HMCR), (2) random consideration (based on 1 − HMCR), and (3) pitch adjustment (based on the PAR). For the memory consideration, the value of the first decision variable ( $x{\prime }_1$ ) has been taken from any of the values as specified in HM range ( $x{\prime }_1-x{\prime }_1^{\mathrm{HMS}}$ ). The other decision variable values ( $x{\prime }_2,x{\prime }_3,\dots ,x{\prime }_K$ ) have also been taken in the same manner

$$x{\prime }_i\leftarrow \{\begin{array}{c}\hfill \begin{array}{cc}\hfill {x\prime }_i\begin{array}{cc}\hfill \in \hfill & \hfill \{x_i^1,x_i^2,\dots ,x_i^{\mathrm{HMS}}\}\hfill \end{array}\hfill & \hfill \mathrm{with}\begin{array}{cc}\hfill \mathrm{probability}\hfill & \hfill \mathrm{HMCR}\hfill \end{array}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}\hfill {x\prime }_i\in {X\prime }_i\hfill & \hfill \mathrm{with}\begin{array}{cc}\hfill \mathrm{probability}\hfill & \hfill (1-\mathrm{HMCR})\hfill \end{array}\hfill \end{array}}\end{array}$$

(14)

Each and every component obtained by the memory consideration is analyzed to determine whether the pitch should be adjusted

$$\mathrm{Pitch\; adjusting\; decision\; for}{x^{\prime }}_i\leftarrow \{\begin{array}{cccc}Yes& with& probabilty& \mathrm{PAR}\\ No& with& probabilty& (1-\mathrm{PAR})\end{array}$$

(15)

If the pitch adjustment decision for $x{\prime }_i$ is Yes, then $x{\prime }_i$ is replaced as follows

$$x{\prime }_i\leftarrow x{\prime }_i\pm \mathrm{rand}\left(\right)*\mathrm{bw}$$

(16)

where bw represents an arbitrary distance bandwidth and rand() represents random number between 0 and 1.

By considering HM, pitch adjustment has been applied to each variable of the new HM vector. These improvisation steps are almost similar to reproduction in GA which uses crossover and mutation operators. Thus, the improvisation utilizes the full HM to construct the solutions, but in case of GA, new chromosomes are generated by crossover of two parents or mutation. While considering the operators in GA, it needs to be carefully designed for highly constrained problems for reaching out its optimal solutions.

Updation of HM

If the new harmony vector, ( $x\prime =x{\prime }_1,x{\prime }_2,\dots ,x{\prime }_K$ ), is better than the worst in the HM, judged in terms of the objective function value, the new harmony is included in the HM, and the existing worst harmony is removed from the HM.

Check for stopping criteria

If the maximum number of improvisations (stopping criterion) is satisfied, computation is terminated. Otherwise steps 2 and 3 will be repeated.

Proposed HSA

In HM vector, each member in the population pool is converted to job schedule using SPV rule–based problem representation discussed in section “Solution representation” for getting the initial sequence of jobs and assignment of jobs to the machines to calculate the objective function. The traditional HSA only uses fixed values for both PAR and bw. These values have been adjusted in the initialization step and not been changed during new generations. This has been considered as the main drawback of traditional HSA. Hence, the need has been raised to sort out these issues. However, by fine tuning the optimized solution vectors using the parameters such as PAR and bw, convergence rate has been increased, and this leads to optimal solution. Still there has been a probability of having a dispute by assigning small and larger values for PAR and bw or vice versa. Hence, the parameter values need to be selected judiciously. Fesanghary et al. ⁴² have developed one of such fine tuning parameters for PAR and bw for reducing the dispute among the values given to their parameters. It has been reported that using equation stated below for PAR and bw, these issues have been sorted out

$$\mathrm{PAR}\left(\mathrm{gn}\right)=\mathrm{PA}{\mathrm{R}}_{min}+\frac{(\mathrm{PA}{\mathrm{R}}_{max}-\mathrm{PA}{\mathrm{R}}_{min})}{\mathrm{NI}}\times \mathrm{gn}$$

(17)

where PAR represents the pitch adjusting rate for each generation, $\mathrm{PA}{R}_{min}$ represents the minimum pitch adjusting rate, $\mathrm{PA}{R}_{max}$ denotes the maximum pitch adjusting rate, NI denotes the number of solution vector generations, and gn represents the generation number.

For bandwidth, tuning has been done using the following equation

$$\mathrm{bw}\left(\mathrm{gn}\right)=b{w}_{max}\exp (c.\mathrm{gn})$$

(18)

where $c=\mathrm{Ln}(b{w}_{min}/b{w}_{max})/\mathrm{NI}$ , bw(gn) denotes the bandwidth for each generation, $b{w}_{max}$ denotes the maximum bandwidth, and $b{w}_{min}$ represents the minimum bandwidth.

Results and discussions

The computation study aims to analyze the performance of MPSO and HSA to minimize makespan and total tardiness for the IPP scheduling. The algorithms are coded in MATLAB R2010a and executed in Intel^® Core^™ i5 CPU M430 at 2.27 GHz with 4 GB RAM. The benchmark problems have been taken from Tanaka and Araki ⁸ and are also available at https://sites.google.com/site/shunjitanaka/pmtt. Fisher ⁵³ has proposed the standard method for generating these problems. The integer processing times p_j (1 ⩽ j ⩽ n) of these problems have been generated from the uniform distribution between [1, 100]. The total processing times have been computed by $P=\sum _{j=1}^n{p}_j$ , the due dates d_j (1 ⩽ j ⩽ n) have been generated using the uniform distribution $[P(1-\tau -R/2)/m,P(1-\tau +R/2)/m]$ , where n is the number of jobs, m is the number of machines, and $\tau$ is the tardiness factor, and due date range is changed by n = 20, m = {2, 3, 4, 5, 6, 7, 8, 9, 10}, τ = {0.2, 0.4, 0.6, 0.8, 1.0}, and R = {0.2, 0.4, 0.6, 0.8, 1.0}. For every combination of m, n, $\tau$ , and R, five problems have been generated. There are five characteristics which are used to represent the problem: number of jobs, number of machines, tardiness, due date range, and position of the instance. For example, the notation of problem 20_03_04_08_001 represents the first problem of three machines, 20-job, τ = 0.4, and R = 8. Among these benchmark problems, 30 problems have been randomly selected for testing and analysis of proposed algorithms. MPSO parameter settings are as follows: the inertial weight w is set as 0.6, both social and cognitive parameters are set as 1.04,and the decrement factor β is set as 0.7. The maximum number of iterations taken for MPSO is 15,000. HSA parameter settings are as follows: HMCR value as 0.9, PAR_min as 0.4, PAR_max as 0.9, bw_min as 0.0001, and bw_max as 1.0. The maximum number of iteration taken for HSA is 5000. For calculating GA, single-point crossover operation and swap mutation operation have been adopted. The crossover rate and the mutation rate are set as 0.85 and 0.1, respectively.

The MPSO algorithm has been tested with GA and PSO. In all the cases, the percentage deviation (PD) of MPSO is better than GA and PSO. HSA has been tested with the same benchmark problem where the PD value of HSA is better than GA, PSO, and MPSO. The results are represented in terms of PD of the solution from the LB reported by Gupta and Ruiz Torres. ¹² The LB has been calculated using the equation given below

$$\mathrm{LB}=\max \{\underset{1⩽i⩽n}{\max }{p}_i;P/m\}$$

(19)

$$\mathrm{PD}=\frac{\mathrm{Best}\mathrm{Makespan}-\mathrm{LB}}{\mathrm{LB}}\times 100$$

(20)

The average percentage deviation (APD) of the proposed MPSO has been tested with GA and PSO in IPPMST. The APD can be calculated as follows

$$\mathrm{APD}=\frac{\sum _L^I\mathrm{PD}\left(L\right)}{I}$$

(21)

where I is the total number of problems and L stands for index of the problem.

From Table 4, it is observed that the calculated value of APD for GA is 0.2686 and PSO is 0.1385. The APD value obtained by using MPSO is 0.0851. Based on APD comparisons, MPSO turns out to be superior to GA and PSO, whereas the performance of GA is worst among all algorithms. The APD value for HSA is 0.0523. It clearly shows HSA is superior to GA, PSO, and MPSO.

Table 4.

APD values for makespan.

m	n	Benchmark problem	LB	GA	PSO	MPSO	HSA	PD of GA = (LB − GA)/LB	PD of PSO = (LB − PSO)/LB	PD of MPSO = (LB − MPSO)/LB	PD of HSA = (LB − HSA)/LB
3	20	20_03_04_08_001	397	433	415	411	402	0.0907	0.0453	0.0353	0.0126
4	20	20_04_02_02_001	298	327	307	304	299	0.0973	0.0302	0.0201	0.0034
4	20	20_04_02_02_002	235	272	251	242	236	0.1574	0.0681	0.0298	0.0043
4	20	20_04_02_02_003	321	335	324	323	322	0.0436	0.0093	0.0062	0.0031
5	20	20_05_02_02_001	238	272	251	246	239	0.1429	0.0546	0.0336	0.0042
6	20	20_06_02_02_001	199	224	205	204	200	0.1256	0.0302	0.0251	0.0050
7	20	20_07_02_02_001	170	191	178	174	172	0.1235	0.0471	0.0235	0.0118
8	20	20_08_02_02_001	149	171	159	156	152	0.1477	0.0671	0.0470	0.0201
9	20	20_09_02_02_001	132	157	141	137	134	0.1894	0.0682	0.0379	0.0152
10	20	20_10_02_02_001	119	145	132	131	125	0.2185	0.1092	0.1008	0.0504
10	20	20_10_02_02_002	95	121	110	106	102	0.2737	0.1579	0.1158	0.0737
10	20	20_10_02_02_003	128	170	151	148	141	0.3281	0.1797	0.1563	0.1016
10	20	20_10_02_02_004	99	134	119	112	111	0.3535	0.2020	0.1313	0.1212
10	20	20_10_02_02_005	100	142	116	110	106	0.4200	0.1600	0.1000	0.0600
10	20	20_10_02_04_001	119	160	142	131	125	0.3445	0.1933	0.1008	0.0504
10	20	20_10_02_04_002	95	132	117	102	99	0.3895	0.2316	0.0737	0.0421
10	20	20_10_02_04_003	128	169	151	148	141	0.3203	0.1797	0.1563	0.1016
10	20	20_10_02_04_004	99	132	119	112	110	0.3333	0.2020	0.1313	0.1111
10	20	20_10_02_04_005	100	130	116	109	106	0.3000	0.1600	0.0900	0.0600
10	20	20_10_02_06_001	119	154	142	131	125	0.2941	0.1933	0.1008	0.0504
10	20	20_10_02_06_002	95	140	118	110	101	0.4737	0.2421	0.1579	0.0632
10	20	20_10_02_06_003	128	159	140	133	130	0.2422	0.0938	0.0391	0.0156
10	20	20_10_02_06_004	99	130	119	112	110	0.3131	0.2020	0.1313	0.1111
10	20	20_10_02_06_005	100	142	116	109	106	0.4200	0.1600	0.0900	0.0600
10	20	20_10_02_08_001	119	142	131	127	125	0.1933	0.1008	0.0672	0.0504
10	20	20_10_02_08_002	95	137	117	102	99	0.4421	0.2316	0.0737	0.0421
10	20	20_10_02_08_003	128	173	151	148	141	0.3516	0.1797	0.1563	0.1016
10	20	20_10_02_08_004	99	132	119	112	110	0.3333	0.2020	0.1313	0.1111
10	20	20_10_02_08_005	100	126	116	109	106	0.2600	0.1600	0.0900	0.0600
10	20	20_10_02_10_001	119	159	142	131	125	0.3361	0.1933	0.1008	0.0504
Average percentage deviation								0.2686	0.1385	0.0851	0.0523

APD: average percentage deviation; LB: lower bound; GA: genetic algorithm; PSO: particle swarm optimization; MPSO: particle swarm optimization with mutation; HSA: harmony search algorithm; PD: percentage deviation.

The improvement rate in APD using MPSO can be defined as follows

$$\mathrm{Improvement\; rate}(\%)=\frac{\left({\mathrm{APD}}_{\mathrm{GA},\mathrm{PSO}}-{\mathrm{APD}}_{\mathrm{MPSO}}\right)}{{\mathrm{APD}}_{\mathrm{GA},\mathrm{PSO}}}$$

(22)

The improvement rate of MPSO with respect to GA is 68.3% and with respect to PSO is 38.5%. The convergence curve for makespan with MPSO is presented in Figure 2.

Figure 2.

Convergence curve for makespan with MPSO for the problem 20_10_02_08_001.

Similarly, the improvement rate for HSA can be defined as follows

$$\mathrm{Improvement\; rate}\left(\%\right)=\frac{\left({\mathrm{APD}}_{\mathrm{GA},\mathrm{PSO},\mathrm{MPSO}}-{\mathrm{APD}}_{\mathrm{HSA}}\right)}{{\mathrm{APD}}_{\mathrm{GA},\mathrm{PSO},\mathrm{MPSO}}}$$

(23)

The improvement performance of HSA with respect to GA is 80.5%, with respect to PSO is 62.2%, and with respect to MPSO is 39.9%. The convergence curve for makespan with HSA is presented in Figure 3. The results of GA, PSO, MPSO, and HSA for makespan and their corresponding APD values are shown in Table 4.

Figure 3.

Convergence curve for makespan with HSA for the problem 20_10_02_08_001.

The same benchmark problem has been tested with total tardiness, and same instances of the problems have been compared with clonal selection PSO (CSPSO) reported by Niu et al. ⁴⁰ and also the optimal solution (OPT) reported by Tanaka and Araki. ⁸ The results for total tardiness are shown in Table 5. The comparison with instances by Niu et al. ⁴⁰ indicates that the total tardiness value by HSA is either exactly matching or better for eight instances out of nine instances. For instances 20_04_02_02_002 and 20_10_02_02_004, total tardiness value by CSPSO are 105 and 180, respectively, and same values are obtained by HSA. For the instance 20_10_02_02_002, the average total tardiness value by CSPSO is reported as 136.6, whereas the corresponding value by HSA is 136. For the instance 20_10_02_04_002, total tardiness values by CSPSO and OPT are 86.5 and 85, respectively, whereas HSA results in total tardiness of 85. For instances 20_10_02_04_004, 20_10_02_06_002, and 20_10_02_06_004, total tardiness values by HSA are better than CSPSO. However, in case of the instance 20_10_02_08_004, CSPSO results in better total tardiness value as compared to HSA. This is to be noted that total tardiness value produced by HSA is similar or superior to GA, PSO, and MPSO. Out of 30 instances, the total tardiness value for 9 instances is comparable with results produced by HSA. The convergence curves for the total tardiness are shown in Figures 4 and 5. PSO, originally proposed by Kennedy and Eberhart, ³⁷ is mainly used in continuous optimization problems due to quick solutions. However, it has an inherent drawback of premature convergence at initial stage of iterations. The limitations of the PSO can be overcome by the use of mutation, a commonly used operator in GA, through introduction of diversity in the solution. The proposed algorithm is termed as MPSO. MPSO is found to perform better than GA and PSO, in both continuous and combinatorial optimization problems. However, a newly developed meta-heuristics known as HSA is simple in concept and easy for implementation. HSA needs only few parameters to be tuned like PSO. Since HSA is a population-based meta-heuristic, the multiple harmonic groups can be used in parallel. This parallelism technique usually leads to better establishment with higher efficiency. Hence, fine balance of intensification and diversification as well as the good combination of parallelism with elitism is the key to the success of HSA and, in fact, to the success of any meta-heuristic algorithm. Thus, HSA performs better than GA, PSO, and MPSO.

Table 5.

Total tardiness value comparison.

m	n	Benchmark problem	OPT	CSPSO	GA	PSO	MPSO	HSA
3	20	20_03_04_08_001	240	–	460	448	432	250
4	20	20_04_02_02_001	144	–	175	158	148	144
4	20	20_04_02_02_002	105	105	122	109	108	105
4	20	20_04_02_02_003	151	–	210	192	161	157
5	20	20_05_02_02_001	149	–	181	163	156	149
6	20	20_06_02_02_001	161	–	193	178	175	168
7	20	20_07_02_02_001	164	–	201	188	178	167
8	20	20_08_02_02_001	177	–	210	192	181	177
9	20	20_09_02_02_001	181	–	205	188	182	181
10	20	20_10_02_02_001	195	–	224	206	201	196
10	20	20_10_02_02_002	135	136.6	152	144	139	136
10	20	20_10_02_02_003	202	–	233	220	206	203
10	20	20_10_02_02_004	180	180	202	195	182	180
10	20	20_10_02_02_005	157	–	200	184	174	160
10	20	20_10_02_04_001	165	–	191	176	167	167
10	20	20_10_02_04_002	85	86.5	129	115	87	85
10	20	20_10_02_04_003	151	–	193	178	162	152
10	20	20_10_02_04_004	162	165.8	186	171	165	162
10	20	20_10_02_04_005	108	–	158	138	135	118
10	20	20_10_02_06_001	161	–	189	174	164	163
10	20	20_10_02_06_002	35	36.7	69	58	39	36
10	20	20_10_02_06_003	122	–	165	147	127	125
10	20	20_10_02_06_004	157	159.4	188	172	159	157
10	20	20_10_02_06_005	83	–	133	125	118	98
10	20	20_10_02_08_001	170	–	223	204	177	173
10	20	20_10_02_08_002	24	24.3	71	63	36	24
10	20	20_10_02_08_003	120	–	156	138	128	122
10	20	20_10_02_08_004	164	169.6	220	199	178	172
10	20	20_10_02_08_005	78	–	128	102	88	81
10	20	20_10_02_10_001	177	–	241	216	204	184

OPT: optimal solution; PSO: particle swarm optimization; CSPSO: clonal selection particle swarm optimization; GA: genetic algorithm; HSA: harmony search algorithm; MPSO: particle swarm optimization with mutation.

Figure 4.

Convergence curve for total tardiness using MPSO for the problem 20_10_02_08_001.

Figure 5.

Convergence curve of total tardiness using HSA for the problem 20_10_02_08_001.

Conclusion

In this article, two new meta-heuristic approaches have been presented for minimizing makespan and total tardiness in IPP scheduling problems. In the first approach, PSO has been modified using a mutation operator and job sequence using SPV rule. In the second approach, HSA has been improved using the Mahdavi et al. ²¹ parameter in PAR and bw, and the sequence of job has been obtained through SPV rule. The scheduling performance by both the algorithms has been compared with GA and PSO. It is found that the performance of MPSO is superior to GA and PSO. However, it is to be noted that HSA outperforms GA, PSO, and MPSO. Both the algorithms happen to be novel in the field of scheduling. The study provides a framework for practicing engineers and managers to schedule large-scale IPP with less computational burden. However, tuning of algorithmic parameters is vital to obtain best schedules. Both the algorithms are extremely simple to implement. Although the proposed algorithms have been tested in IPPMST, the approach can be extended to flow shop, job shop, flexible low shop, flexible job shop, FMS, and resource-constrained project scheduling problem. The mapping of continuous position values of PSO or HM vector of HSA to schedule can be made with application of SPV rule with little modification. In case of continuous optimization problems, for example, process optimization, both the proposed algorithms can be applied without the need of application of SPV rule. However, combinatorial problems require the generation of combination by the use of SPV rule. Application of the proposed algorithms to combinatorial case makes problems more complex than continuous counterpart, and preservation of diversity of solution becomes an important issue. In this study, the processing time is assumed to be deterministic and machines are perfectly reliable. In future, the algorithms may be modified to incorporate stochastic processing times and machine unreliability for solving IPPMST. Although the study applies SPV rule for generation of job schedule, other mapping methods may be explored in conjunction with the proposed algorithm to study the performance of the algorithms in IPPMST.