A global-local neighborhood search algorithm and tabu search for flexible job shop scheduling problem

General description of the GLNSA

Many recent works have proposed hybrid optimization algorithms that combine a population method with another local search method. These methods have produced satisfactory results, and new proposals continue to be developed to have equally efficient algorithms, with a simple implementation and lower complexity of execution. In these methods, the population part applies a series of information sharing and mutation operators to improve the optimization process. Usually, these methods are used serially to each individual to obtain a new position.

In the algorithm proposed in this work, a different strategy is taken, inspired by the operation of a cellular automaton, where the solution of a complex problem is achieved through simple operators’ cooperation. This inspiration has been proved before in the design of infinite-impulse filters (Lagos-Eulogio et al., 2017) and a concurrent layout and scheduling problem (Hernández-Gress et al., 2020), obtaining good results.

A cellular automaton is a discrete dynamic system made up of indivisible elements called cells, where each cell changes its state over time. The state change can depend on both the current state of each cell and its neighboring cells. With such simple dynamics, cellular automata can create periodic, chaotic, or complex global behavior (Wolfram, 2002; McIntosh, 2009; Adamatzky, 2010; Bilan, Bilan & Motornyuk, 2020).

In this work, the idea of a cellular-automaton neighborhood is an inspiration to propose a new algorithm that optimizes instances of the FJSP. The algorithm has S_n leading solutions or smart-cells, where each of them first performs a global search mainly focused on making modifications to the sequencing of operations, applying several operators to form a neighborhood, and selecting the best modification. Then, each smart-cell executes a local search focused on machine assignment, applying a TS with increasing iterations at each step of the algorithm. Thus, in each iteration, a global exploration search and a local exploitation search are performed to optimize instances of the FJSP, and hence, this process is called the global-local neighborhood search algorithm (GLNSA).

The novel part of the GLNSA is that instead of serially applying the sharing-information and mutation operators to each smart-cell to obtain a new solution as most recently proposed hybrid methods do, in the GLNSA, a neighborhood of new solutions is generated. Each neighbor is the product of applying a random information exchange or mutation operator to the original smart-cell. Then the best neighbor is selected to replace it, following a neighborhood structure inspired by cellular automata.

In this way, the contribution of GLNSA is to show that notions of cellular automata are helpful to inspire new forms of hybridization in the conceptualization of new optimization methods for the FJSP, obtaining satisfactory results as shown in the subsequent sections. The corresponding flow chart is described in Fig. 1.

The detailed explanation of the encoding and decoding of solutions, the neighborhood used, its operators for the global search, and the TS operators are presented in the following sub-sections. The general procedure of the GLNSA is in Algorithm 1.

Encoding and decoding of solutions

To represent a solution for an instance of the FJSP, we take the encoding with two strings (OS and MS) described in Li & Gao (2016), OS for operations and MS for machines.

The string OS consists of a permutation with repetitions where each job J_i appears n_i times. The string OS is read from left to right, and the j−th appearance of J_i indicates that the O_i,j operation of job J_i should be processed. This coding of the sequence of operations OS has the advantage of any permutation with repetitions producing a valid sequence so that the operators used in this work will always yield a feasible sequence.

The string MS consists of a string with a length equal to the number of the total operations. The string is divided into n parts, where the i−th part contains the machines assigned to J_i with n_i elements. For each i−th part of MS, the j−th element indicates the machine assigned to O_i,j.

Initially, each solution OS is generated at random, making sure that each job J_i appears exactly n_i times. For every operation, one of the possible machines that can process it is selected at random as well.

For each pair of OS and MS strings that represents a solution for an instance of the FJSP, their decoding is done using an active scheduling. For each operation O_i,j in OS and its assigned machine k in MS, its initial time s(O_i,j) is taken as the greater time between the completion time of the previous operation O_{i,j − 1} and the lesser time available in the machine k (not necessarily after the last operation programmed in that machine) where there is an available time slot, such that the processing time p_i,j,k is less than or equal to the size of this slot. The time the O_i,j operation is completed is called C(O_ij), and for j = 1, s(O_i,j − 1) = 0 for all 1 ≤ i ≤ n. An example with 3 jobs and 2 machines is shown in Fig. 2.

In Fig. 2 (A), there are 3 jobs, each with 2 operations. Almost every operation can be executed on the 2 available machines. On average, each operation can be executed by 1.83 machines; this is called the system’s flexibility. Note that the time at which each operation is performed on each machine may be different. In part (B), you can see how a smart-cell is encoded. It consists of two strings; the first OS is a permutation with repetitions, where each job appears 2 times. The second string MS contains the machines programmed for each operation, where the first 2 elements correspond to the machines assigned to the operations O_1,1 and O_1,2, the second block of 2 elements specify the machines assigned to operations O_2,1 and O₂, and so on. Finally, part (C) indicates the decoding of the smart-cell reading of the string OS from left to right. In this case, the O₂₂ operation, which is the fifth task programmed in OS, is actively accommodated, since there is a gap in the M₁ machine between O₁₁ and O₃₂, which is of sufficient length to place O₂₂ just after the preceding operation on the machine M₂ and without moving the operations already programmed in M₁. This scheduling gives a final makespan of 7 time units using the active decoding.

Selection method

The GLNSA uses elitism and tournament to refine the population, by considering the value of the makespan of each smart-cell. For elitism, a proportion E_p of smart-cells is taken with the best values of population’s makespan. Those solutions will remain unchanged in the next generation of the algorithm. The rest of the population members are selected using a tournament scheme, where a group of b smart-cells are randomly selected from the current population and competed with each other, and the smart-cell with the best makespan is selected to become part of the population to be improved using the global and local search operators described later. In this work, we take b = 2. This mixture of elitism and tournament allows a balance between exploring and exploiting the information in the population, keeping the best smart-cells, and allowing the other smart-cells with a good makespan to continue in the optimization process.

Neighborhood structure

In this section, two neighborhood structures are presented for the FJSP. The first neighborhood focuses mainly on the sequencing of operations and random assignments of machines to each operation, and the second neighborhood focuses on a local search neighborhood to improve the assignment of the machines once the sequence of operations has been modified.

For the second neighborhood, a simplification of the Nopt1 neighborhood and the makespan estimation explained in Mastrolilli & Gambardella (2000) are used to improve the GLNSA execution time.

Exploration neighborhood

For the global search neighborhood, well-known operators used in various task sequencing problems are used. These operators are used to optimize the sequence of operations. For the machine assignment, the mutation operator used in Li & Gao (2016) is applied. Each smart-cell generates l neighbors using one of the three possible operators (insertion, swapping, or PR) with probability α_I, α_S and α_P, respectively, to generate a variant of OS. For the machine assignment, the mutation operator with probability α_M is employed to generate another variant of MS. From these l neighbors, the one with the smallest makespan is chosen to be the new smart-cell. This neighborhood is exemplified in Fig. 3.

Insertion

The insertion operator consists of selecting two different positions, k₁ and k₂, of the string OS in a smart-cell to obtain another string OS′. For example, if k₁ > k₂ with OS = (O₁, … Ok₂ … O k₁ …), then we get the string OS′ = (O₁, … Ok₁ Ok₂ … Ok₁ − 1 …). This is analogous if k₁ < k₂. This is exemplified in Fig. 4.

Swapping

The swapping operator consists of selecting 2 different positions, k₁ and k₂, from the string OS to exchange their positions. For example, if k₁ < k₂ with OS = (O₁, … Ok₁ … Ok₂ …), then the string OS′ = (O₁, … Ok₂ … Ok₁ …). This is analogous if k₁ > k₂. This is exemplified in Fig. 5.

Insertion and swapping are operators classically used in metaheuristics for task scheduling problems, since they provide good quality solutions (Błażewicz, Domschke & Pesch, 1996; Cheng, Gen & Tsujimura, 1999; Deroussi, Gourgand & Norre, 2006).

Path relinking

Path relinking (PR) involves establishing a route between two different smart-cells for their strings OS and OS′. The route defines intermediate strings ranging from OS to OS′. To form this route, the first position k of OS is taken such that OS_k ≠ OS′_k. Then, the first position p is located after k such that OS_p = OS′_k. Once both positions have been found, the positions k and p are exchanged in OS to obtain a new string closer to OS′. The PR is repeated until OS′ is obtained. In the end, one of the intermediate solutions is randomly chosen.

The idea of the PR is to generate solutions that combine the information of OS and OS′ and fulfill two objectives in the GLNSA: if both smart-cells have similar machine strings, PR acts as a local search method that refines the strings of operations. Conversely, if both smart-cells have very different machine strings, PR works as an exploration method that generates new variants of one of the smart-cells, taking the other as a guide, as shown in Fig. 6. PR has already been used successfully in the FJSP, as shown in its application with different neighborhood variants in González, Vela & Varela (2015).

Mutation of the string of machines

For the MS string of each smart-cell, a mutation operator is applied that selects half of the positions in MS at random. For each position linked to an operation in OS, the assigned machine is changed for another selected at random so that it can perform the respective operation in OS.

Figure 7 shows an example of the mutation operator applied to a sequence MS. Three machines are selected corresponding to the operations O₁₁, O₂₁, and O₂₂ in OS. The new machine assignment is obtained by changing the machines in these positions to others that can also perform these operations.

Local search neighborhood

Simple local search methods are based on creating slight modifications of a solution to generate an improved one until the process can no longer achieve an improvement, or a certain number of repetitions is met. These methods tend to get trapped in local minima or regions where many solutions have the same cost. In order to make an algorithm more adaptable and capable of escaping local optima, various strategies have been proposed. One of the most effective ones so far is the TS.

TS is a metaheuristic developed by Glover & Laguna (1998) to guide a local search more effectively by incorporating an adaptive memory and a more intelligent exploration. It is one of the most used methods and offers satisfactory results in many optimization problems (Glover & Kochenberger, 2006). TS allows, during the search process, solutions with a worse cost to be accepted to provide a more significant search capacity. A short-term memory prohibits returning to solutions that have already been recently explored. They are again taken if they are outside the tabu threshold or meet an aspiration criterion (Chaudhry & Khan, 2016). In this way, the search does not have a cyclical behavior and produces new trajectories in the solution space.

In the GLNSA, TS allows generating solutions that may not improve the makespan of the original smart-cell as long as the operation and the machine selected to obtain the new solution are not forbidden. TS keeps a record of the operation and the machine selected in each movement and the threshold at which this movement will remain tabu, and the aspiration criterion allows a solution to be accepted even if it is tabu.

The implementation of TS is resumed using a simplification of the neighborhood structure proposed in Mastrolilli & Gambardella (2000) and the makespan estimation explained in the same work to reduce the computational time of TS. The general TS procedure is described in Algorithm 2.

The flow chart of the TS is shown in Fig. 8.

Critical path

The TS for this job uses the following definition of the critical path cp. To form cp, one of the operations with completion time equal to the makespan is selected randomly. Once this operation has been chosen, the operation that precedes it is selected either on the same machine or by the previous operation of the same job, and the one whose completion time is equal to the initial time of the current operation is selected. If both previous operations have the same completion time, one of them is selected randomly. This process is repeated until an operation with a start time of 0 is selected. These operations form a critical path cp of length q.

Simplified Nopt1 neighborhood

The neighborhood used for the local search is based on the one defined as Nopt1 in Mastrolilli & Gambardella (2000). Given a critical path cp, in the original neighborhood Nopt1, for each operation in cp, a set of preceding and succeeding operations can be found in each feasible machine, such that a new placement of the operation between these operations on the new machine optimizes the makespan. The calculation of Nopt1 depends on the review of the start and tail times of the operations programmed in each machine, implying an almost constant computational time. When this operation is performed several times (such as in a meta-heuristic algorithm), the computational time can increase considerably, especially if the number of jobs and machines is large and the system has high flexibility.

In this work, it is proposed to use a simplification of the Nopt1 neighborhood, where a search for the best position is not carried out to accommodate each operation in a new feasible machine, but its machine assignment is just changed in the string MS, preserving the position in the string OS. The idea is that the global search operations applied in the neighborhood based on cellular automata (insertion, swapping, and PR) are capable of finding this optimal assignment as the optimization algorithm advances, especially for systems with greater flexibility, where more machines are capable of processing the same operation. Thus, the objective is to avoid carrying out the operations that explicitly seek to accommodate an operation on a machine and directly take the position that is being refined by the global neighborhood’s operations.

There are more optimal job placement positions for systems with high flexibility, given the greater availability of feasible machines, so this simplification focuses on showing that good results can be obtained with a reduced process for this type of FJSP instance.

Parameters of the GLNSA

The following are the parameters of the proposed algorithm:

• Number of iterations for the whole optimization process: G_n.

• Number of smart-cells: S_n.

• Global neighborhood size: l.

• Probabilities α_I, α_S, and α_P for insertion, swapping, and PR, respectively, and probability α_M for mutation.

• Maximum number of stagnation iterations: S_b.

• Proportion of elitist solutions: E_p.

• Number of tabu iterations for every optimization iteration: T_n.

• Tabu threshold: T_u.

Parameter tuning

A preliminary study was carried out considering different values of the GLNSA parameters. They were applied to the same problem using a similar number of iterations to select the best parameters applied to instances of FJSP with high flexibility.

For the population size S_n, the number of iterations G_n, and the mutation probability α_M, the results presented in González, Vela & Varela (2015) and Li & Gao (2016) are taken as a basis, since they are recent works that show great effectiveness both in the makespan calculation as well as in the runtime for FJSP instances.

For G_n, values of 200 and 250 were tested. Also, while S_n is taken between 20 and 100 in González, Vela & Varela (2015) and defined as 400 in Li & Gao (2016), we tested S_n between 40 and 80.

The number of neighborslthat each smart-cell has in our algorithm to generate the global search neighborhood was tested with values between 2 and 3, in order to preserve a population close to 250 solutions at most (number of smart-cells by number of neighbors) and keep a computational execution close to the cited references. To form the global neighborhood, the probability combinations (α_I, α_S, α_P) with respective values (0.5,0.25,0.25), (0.25,0.5,0.25) and (0.25,0.25, 0.5) were tested. The mutation probability α_M was taken with values 0.1 and 0.2.

To control the stagnation limit S_b, the value proposed by Li & Gao (2016) is taken to test S_b with values 20 and 40, and for the elitist proportion of solutions, a value of E_p of 0.025 and 0.05 is considered.

Without a doubt, the TS is the most computationally expensive process that the proposed algorithm has. In González, Vela & Varela (2015), the execution of the TS is tested up to 10,000 iterations per solution. In Li & Gao (2016), this number goes up to about 80,000 times per solution, of course, with different ways of creating solutions and estimating the makespan. Our algorithm uses the makespan estimation developed in Mastrolilli & Gambardella (2000) during the TS to reduce the computational time of the optimization process.

In our algorithm, we take a point of view similar to Li & Gao (2016), using an increasing number of iterations of TS, as the number of iterations of the optimization process grows. For each iteration i, T_n * i TS iterations will be taken for each smart-cell, where T_n values of 1 and 2 are tested to keep a maximum TS close to 60,000 iterations per smart-cell. Altogether, it took 384 different combinations of parameters to tune the GLNSA.

For this preliminary tuning study, we took the problem la31 with more flexibility (known as part of the vdata set), initially proposed by Hurink, Jurisch & Thole (1994), and downloaded from https:/people.idsia.ch/monaldo/fjsp.html. From the tuning study, the parameters G_n = 250, S_n = 40, l = 2, α_I = 0.5, α_S = 0.25, α_P = 0.25, α_M = 0.1, S_b = 40; E_p = 0.025 and T_n = 1 were selected as the most appropriate values to apply the GLNSA.

The threshold used in the tabu list for each entry (operation/machine) consists of the sum of the length of the random critical path plus the number of feasible machines that can perform the corresponding operation. This criterion was also proposed in Mastrolilli & Gambardella (2000) and has been widely used by similar algorithms.

Figure 9 shows the convergence of the makespan by applying the GLNSA to the parameters indicated above on the instance la31-vdata. The implementation was developed in Matlab (the implementation characteristics and computational experiments are specified in detail in the next section), and the execution time was 56 s.

Comparative analysis of computational complexity between algorithms

Since it is difficult to compare the execution time of each of these algorithms since they were tested in different architectures, languages, and programming skills, these algorithms are compared based on their computational complexity $\mathcal{O}(o)$, taking as a reference the total number of operations $o={\sum }_{i=1}^n{\sum }_{j=1}^{n_i}{O}_{i,j}$ to be processed for an instance of the FJSP. To generalize the way of measuring each algorithm’s complexity, we will use the notation defined for the GLNSA parameters; the number of iterations will be denoted by G_n and the number of solutions handled by each algorithm will be defined by X, which in the case of the GLNSA is X = S_nl. For the local search each algorithm applies, the number of iterations is expressed with T_n, and for the number of machines, m will be used.

The ePSOGSA first uses a modified PSO to update each individual by rate or mutation, then applies a modified GSA to obtain an extended population and sort a final population. The ePSOGSA does not use an iterative local search based on the calculation of critical paths. The GRASP algorithm constructs an optimized Gantt chart and then applies a greedy local search that is quadratic concerning the number of operations. The HA has an architecture similar to the GLNSA proposed in this work, and first applies four types of genetic operators to each solution. Then, HA applies a TS that is supported by calculating the critical path and selecting different machines for each operation. The IJA uses a modified Java algorithm that applies three change operators to each solution and then, in the same cycle, a local search to all individuals based on the critical path and random exchange of critical blocks in the machines, obtaining a complexity very similar to HA. The SLGA uses two genetic operators for each individual and then uses greedy reinforcement learning based on selecting suitable individuals and two types of actions that also do not need to calculate a critical path. The SSPR algorithm also has a similar structure to HA and GLNSA. It first uses a TS for each individual and then uses PR between two pairs of solutions iteratively to obtain a new individual, which is again improved by a TS based on the critical path calculation. The SSPR also employs a diversification phase when all pairs of individuals have been selected for the PR. The TLBO algorithm is based on a teaching-learning algorithm that uses real coding that utilizes three learning rules and must review the feasibility of the solutions. It also applies a local search based on the exchange of critical operations and their change of machines. The TlPSO algorithm applies two modifications to the PSO in order to optimize the routing problem in the first stage, and within this stage, optimizes the machine assignment problem with another PSO. The last algorithm taken for comparison was the VNSGA, which first applies three genetic operators to each individual and then performs a local search based on the critical path of each solution where processing-time conditions are detected for four neighborhood types available for the machines assigned to each solution.

Finally, the algorithm proposed in this work (GLNSA) performs an elitist selection of solutions, then a search for l based on simple operators such as insertion, swapping, PR, and machine mutation on the set of S_n smart-cells and then a TS only on S_n smart-cells. Since S_n < X, the GLNSA performs less computation for the local search.

Table 1 presents the algorithms used for comparison in this work ordered depending on their complexity, taking as a reference that S_n < X <T_n, since in the algorithms that apply a local search, a high value of T_n is usually chosen for best results. We can see that the GLNSA has a competitive computational complexity compared to state-of-art algorithms recently proposed in the literature for the FJSP.

It should be noted that this analysis only considers the computational complexity required to manipulate and modify the scheduling and routing of the sequences of operations of an FJSP instance. The computational complexity for the makespan calculation is $\mathcal{O}(o^2)$ and the makespan estimation taking into account only critical operations is bounded by $\mathcal{O}(o)$. These processes are not contemplated in the analysis presented since all the algorithms use these operations, and we only focus on the study of the computational process that distinguishes each algorithm.

First experiment, Kacem instances

Table 2 presents the GLNSA results and the comparison with algorithms that present results for the Kacem benchmark. n represents the number of jobs and m the number of machines, in addition to the flexibility rate β of each instance. This benchmark dataset is characterized by its high flexibility and starts from instances with low dimensionality to instances with a greater number of jobs and machines.

For this benchmark, only the IJA and TLBO algorithms report complete results (Caldeira & Gnanavelbabu, 2019; Buddala & Mahapatra, 2019). In this experiment, the GLNSA obtains the best-known results for the makespan value in all cases, like the IJA, and improves the TLBO for the K5 instance with greater dimensionality. This experiment corroborates the excellent performance of the GLNSA for problems with a high rate of flexibility β.

Second experiment, Brandimarte instances

To compare the results of the nine algorithms in this benchmark, the relative percentage deviation (RPD) is defined in Eq. 7.

(7) $$RPD={\displaystyle \frac{BOV-BKV}{BOV}\times 100}$$

where BOV is the best value obtained by the algorithm, and BKV is the best-known value for each instance, in the Brandimarte benchmark, and the best-reported values are taken from Chen et al. (2020). Table 3 presents the results of the GLNSA and its comparison with the other algorithms. In this table, the rate β of each instance is presented, which varies from 0.15 to 0.35. All cases have partial flexibility. The best makespan obtained by any of the algorithms is indicated with a *. In this experiment, it can be observed that the GLNSA obtains the best makespan 6 times, only behind the ePSOGA and the SSPR, and shares the number of best results obtained with the GRASP, HA, and TlPSO algorithms.

Table 4 presents the average RPD of each algorithm in the 10 instances and the ranking of each algorithm, taking the average RPD as a reference. The GLNSA is seen to obtain the fifth position among the algorithms taken for comparison, which shows its competitiveness with state-of-the-art algorithms for this benchmark.

A non-parametric Friedman test was performed with the RPD values in all instances to corroborate whether there is a statistically significant comparison between the results obtained by all the algorithms (Derrac et al., 2011). The result of the Friedman test is observed in the last column of Table 4 with a value of p = 0.0003, which is less than the significant level of 0.05 to reject the null hypothesis that algorithms behave statistically similarly. The p value demonstrates significant differences in the performance of the 10 algorithms, showing the GLNSA’s competitiveness to optimize this benchmark.

Third experiment, Hurink instances with low flexibility

This experiment takes 43 instances from the Hurink benchmark (rdata) with low flexibility, whose rate β ≤ 0.4. The results of the GLNSA are compared taking the BKVs of every instance as reported in Li & Gao (2016). Table 5 shows the comparison of the GLNSA with the previous algorithms that report results for this benchmark (ePSOGA, HA, IJA, and SSPR); the ePSOGA only reports 15 results from the 43 instances, so only these results were taken to complete the comparative analysis. The results marked with * are the best obtained among the five algorithms.

Table 6 presents the average RPD of the algorithms that report complete results for the 43 instances and the ranking of each algorithm, taking the average RPD as a reference. The GLNSA is seen to obtain the fourth position among the algorithms taken for comparison, and the Friedman test obtains a value of p = 0.000000002, which is less than the significant level of 0.05, proving that there are significant differences in the performance of the 4 algorithms. The low competitiveness of the GLNSA with respect to the other algorithms can be explained by the results in instances la21 to la40, which correspond to the examples with higher dimensionality (number of jobs and machines) and lower flexibility. This is expected given the functioning of the GLNSA that is focused on resolving FJSP instances with a higher rate β.

However, if we analyze only the problems with greater flexibility (la01 to la15) to consider the ePSOGSA algorithm and do the same average RPD examination and Friedman test, we obtain the result shown in Table 7. For these instances with more flexibility, the GLNSA ranks second among the 5 algorithms with a value of p = 0.0339, which is less than the significant level of 0.05, demonstrating significant differences in performance of the 5 algorithms for these instances. The preceding results corroborate the efficiency of the GLNSA to optimize FJSP instances with high flexibility.

Fourth experiment, Hurink instances with greater flexibility

The fourth experiment takes 43 instances of the Hurink-vdata benchmark with rate β = 0.5, indicating that an operation can be processed by around half the machines, already implying a high degree of flexibility. The results of the GLNSA are compared with the algorithms taken before that report results for this benchmark (HA, IJA, and SSPR) in Table 8, taking the BKVs for statistical analysis again from Li & Gao (2016).

Table 9 presents the average RPD of the algorithms and their ranking. According to this average, the GLNSA obtains the third position among the algorithms taken for comparison. Two non-parametric Friedman tests were made, one among the 4 algorithms and the other only comparing the GLNSA with the HA, which are the ones that had the closest values.

In Table 9, a value of p = 0.0000003 is obtained, which is less than the significant level of 0.05, proving that there are significant differences in the performance of the 4 algorithms. However, when comparing only the GLNSA to the HA, a value of p = 0.5637 is obtained, which rejects a significant difference between both algorithms to optimize this benchmark dataset. This test verifies that the efficiency of the GLNSA is similar to the HA and is only significantly below the SSPR. These results confirm the observation that the simplified Nopt1 neighborhood works adequately for problems with greater flexibility. From the results of the four sets of experiments, the GLNSA has obtained results comparable with state-of-art algorithms and with a competitive computational complexity, especially for problems with a high rate of flexibility.