2.1. Flexible Job-shop Scheduling Problem

The FJSP can be described as follows. There are n jobs J = {J₁, J₂,...,J_n} to be processed on m machines M = {M₁, M,...M_m}. Each job J_i has n_i operations {O_i,1, O_i,2,...,O_i,ni} to be processed according to a given sequence. Each operation O_i,j can be processed on any machine among a subset M_i,j ⊆ M. The FJSP is to solve the assignment of machines and the sequence of operations to minimize a certain scheduling objective, e.g., the makespan of all the jobs (C_max).

Moreover, the following conditions should be satisfied while processing. Each machine processes one operation at a given time. Each operation is assigned to only one machine. Once the process starts, it cannot be interrupt. All jobs and machines are available at the beginning. The order of the operations for each job is predefined and cannot be modified.

For explaining explicitly, an example of FJSP is shown in Table 1. There are 2 jobs and 4 machines, where the rows correspond to the operations and the columns correspond to the machines. Each element denotes the processing time of this operation on the corresponding machine, and the “-” means that an operation cannot be processed on the corresponding machine.

2.2. Standard Artificial Fish Swarm Algorithm

The artificial fish swarm algorithm (AFSA) is a population-based optimization algorithm, which is inspired from fish swarm behaviors. In an AFSA system, each artificial fish (AF) adjusts its behavior according to its current state and its environmental state, making use of the best position encountered by itself and its neighbors. The optimization of artificial fish swarm algorithm is conducted by four behaviors, i.e., preying, swarming, following, and moving.

Suppose X_i = (x₁,x₂,...,x_n) is the current position of artificial fish AF_i; Y_i = f(X_i) is the fitness function at position X_i. Visual is the visible distance of AF; try_number is the try times of preying behavior; Step is the maximum moving step of AF; δ is the crowd factor; n_f is the number of AFs within its visual. For AF_i, one target position X_j in its visual can be described by Eq. (1), Rand() is a function that generate random numbers in the interval [0,1]. Then the AF_i updates its state by using Eq. (2) when the updating condition is satisfied.

The four behaviors of AF are described as follows:

1. (1)

   Preying: The AF_i chooses a position X_j randomly within its visible region using Eq. (1). If Y_j < Y_i, it moves one step to X_j according to Eq. (2). Otherwise, it chooses another position and determines whether it satisfies the requirement Y_j < Y_i. If the requirement is still not satisfied after try_number times, AF_i executes the moving behavior.

2. (2)

   Swarming: Suppose X_C is the center position in the visible region, if the center has more food and low crowd degree as indicated by Y_C · n_f < Y_i · δ, then AF_i moves a step towards X_C according to Eq. (2). Otherwise, AF_i executes default preying behavior.

3. (3)

   Following: Suppose X_b is the best found position in the visible region, if the position X_b has high food consistence and low crowd degree as indicated by Y_b · n_f < Y_i · δ, then AF_i moves a step towards X_b according to Eq. (2). Otherwise, AF_i executes default preying behavior.

4. (4)

   Moving: AF_i choose a random position in its visual region and moves a step towards this direction. It is a default behavior of preying.

In AFSA, swarming and following are simulated in each generation. The AFs will choose the behavior to find the position with better fitness value, and the default behavior is preying. The flow chart of AFSA is shown in Fig. 1.

Fig. 1.

The framework of basic AFSA

3.1. Pre-principle and post-principle Mechanisms

The FJSP optimizes the objective function by adjusting the machine assignment and the operation sequence. The order of solving the two sub-problems may affect the optimal results. Thus, we propose a pre-principle arranging mechanism and a post-principle arranging mechanism to adjust machine assignment and operation sequence. The objective is to enhance the diversity of population and make the algorithm searching the feasible solutions more comprehensively. The proposed mechanisms work as follows:

1. (1)

   Pre-principle arranging mechanism (PrA): The machine assignment part is firstly adjusted for balancing the workload of each machine, and then the operation sequence part is adjusted for minimizing the total makespan.

2. (2)

   Post-principle arranging mechanism (PoA): The operation sequence part is firstly adjusted for balancing the completion time of each job, and then the machine assignment part is adjusted for minimizing the total workload of machines.

While implementing the AFSA, the population is divided into two sub-populations. The two sub-populations respectively adopts the PrA and PoA mechanisms. After the independent evolution for each sub-populations, they are recombined to an entire population for further evolution.

3.2. The preying behavior based on estimation of distribution

In an AFSA system, preying behavior usually tends to be blindfold since the selection of destination locations is achieved by a random process. To overcome this shortcoming, we propose to estimate the distribution of the individuals, and then use the distribution model to guide preying behavior. The estimation of distribution algorithm (EDA) can reduce the randomness of behavior and make the search move toward and converge to the promising regions in the solution space.

The EDA works as follows: (1) Select a set of promising individuals from the population according to the fitness value; (2) Estimate the probability distribution of the selected individuals according to a probabilistic model. The probability distribution is constructed by two matrixes, i.e., the machine probabilistic matrix and the operation probabilistic matrix; (3) Generate new individuals according to the estimated probability.

Let I denote the number of promising individuals selected from the current population, and D1l be the machine probabilistic matrix at the lth generation. Each entry dijkl of D1l means the probability of operation O_i,j being processed on machine M_k and it is determined by the following formula:

Let D2l denote the operation probabilistic matrix at the lth generation. Each entry dijlof D2l means the probability of job J_j being arranged in position i in the operation sequence and it is calculated by the following formula:

While implementing the preying behavior, a new individual is generated by sampling the two probabilistic matrixes. The machine assignment vector is generated through sampling the probabilistic matrix D1l. For each operation O_i,j, machine M_k is selected with a probability of dijkl. Similarly, the operation sequence vector is generated by sampling the probabilistic matrix D2l. Job J_j is selected with a probability of dijl to replace the j-th position of operation sequence vector.

3.3. Attracting behavior

In AFSA, each AF determines its next position according to current state and its environmental state within its visible region, which may limit the exploration ability and interaction with the other AFs outside its visible region. We propose an attracting behavior to enhance exploration ability.

The bulletin board that records the state of the current optimal individual is setup. For each AF, it reads the position information of the optimal individual from the bulletin board, and then it moves one step toward this direction.

Suppose X_i is the current position of AF_i, Y_i, is the fitness value. X_best is the position of the global optimal individual and its fitness value is Y_best. If Y_i > Y_best, then AF_i, moves towards X_best for a step according to Eq. (2). Otherwise, if Y_i > Y_best, which means that it is the optimal individual, then it executes default preying behavior.

3.4. Critical Path Local Search based on Public Factor

In this subsection, a local search procedure is presented to enhance the exploitation around the best solution obtained by the AFSA. The local search is based on the critical path, so it is applied to the schedule represented by the disjunctive graph rather than by the AF. Hence, when a solution is to be improved by the local search, it should be firstly decoded to a schedule represented by the disjunctive graph.

3.4.1. The disjunctive graph

A feasible solution of FJSP can be represented by the disjunctive graph G = (N,B ∪ C), where N is a set of nodes which includes all the operations and two dummy nodes: starting and terminating; B is a set of conjunctive arcs which represents precedence constraints in the same job; C is a set of disjunctive arcs which correspond to the precedence of operations processed on the same machine. The weight of each node is the processing time of corresponding operation. In a disjunctive graph, the longest path form the starting node to terminating node is called critical path, whose length determines the makespan of the schedule. Any operation on the critical path is called critical operation.

Take the problem shown in Table 1 for instance, a possible schedule represented by the disjunctive graph is showed in Fig. 2. S_G and E_G are respectively dummy starting and terminating nodes. The operations O_1,1 and O_2,3 are performed on M₁ successively, O_2,1 and O_1,2 are performed on M₃ successively, and O_2,2 is processed in M₂.

Fig. 2.

Illustration of framework of disjunctive graph

3.4.2. Neighborhood structure based on public factor

The disjunctive graph usually has more than one critical paths. Only changing the length of all the critical paths, the makespan can be changed. For obtaining a better schedule from the current one, lots of operations may be tried to move. The process is time consuming. So the public factor based critical path search in the neighborhoods is proposed. The public factor is used to identify the influence degree of the critical for all the critical paths. The public factor of an operation O_i,j can be defined by the formula:

(5)Pfi,j=Ni,j/Tcp

Where T_cp is the total number of critical paths, and N_i,j is the number of the critical paths including the operation O_i,j. Pf_i,j is in the interval [0,1]. Pf_i,j = 0 means this operation does not include any critical path, in other words, it is not a critical operation; Pf_i,j = 1 means this operation is included in all critical paths. For each critical operation, the higher value of public factor one operation possess, the greater impact for disjunctive graph it has while moving it. In the neighborhood structure, the critical operation with the highest public factor will be moved preferentially.

While moving an operation, the precedence constrains should be satisfied. For an operation O_i,j processed on M_k, we define SE_G(O_i,j) as the earliest start time and CE_G(O_i,j) = SE_G(O_i,j) + t_{i, j, k} as the earliest completion time. Similarly, denote the latest start time without delaying the makespan as SL_G(O_i,j) and the latest completion time as CL_G(O_i,j) = SL_G(O_i,j) + t_{i, j, k}. Let PJ_i,j = O_{i, j−1} be the precedent of operation O_i,j and FJ_i,j = O_{i, j+1} be the successor of operation O_i,j. Denote PM_i,j as the operation performed on M_k right before O_i,j and SM_i,j as the operation performed on M_k right after O_i,j. In the disjunctive graph G, the process of moving an operation O_i,j is to delete it from its current machine sequence by moving all its disjunctive arcs and then insert it at another available machine by adding disjunctive arcs. Let PO_l (l = 1,2, ..., x) be the critical operation to be moved, where x is the total number of critical operations in G. Let G⁻ be the disjunctive graph obtained by deleting the critical operation PO_l from G. For no increasing the makespan after inserting operation PO_l, we take C_max(G) as the makespan of G⁻ when we calculate the latest start time SL_G⁻ (O_i,j) for each operation in G⁻. If PO_l is inserted before O_i,j on machine M_k in G⁻, it could be started as early as CE_G⁻ (PM_i,j) and should be completed as late as SL_G⁻ (O_i,j) without delaying the required makespan C_max(G). In addition, PO_l needs to comply with the operation precedence constraints. So, the available idle time for inserting PO_l to machine M_k need to satisfy the following condition:

(6)max⁡{CEG−(PMi,j),CEG−(PJPOl)}+tPOl,k                      <min⁡{SLG−(Oi,j),SLG−(FJPOl)}

This moving process is repeated until a better schedule strategy is found or all critical operations have been tried to move. The procedure of local search is shown in Algorithm 1.

1.	Convert the feasible schedule to the disjunctive graph G
2.	Get all the critical operations in the G
3.	Calculate the Pf of all the critical operations and sort them in descending order to form a set {PO1, PO2,...POx}
4.	for i = 1 to x do
	Delete POi from G to get G−
	Calculate all the idle time intervals in G−
	if existing an available time interval for POi
	then
	Insert the operation POi to get Gi
	Return the disjunctive graph Gi
	end if
	end for
5.	Return the final disjunctive graph.

Algorithm 1.

The procedure of critical path local search

Figure 3 gives an illustration of the local search for the example in Table 1. In Fig. 3 (a), the only critical path is S_G → O_1,1 → O_2,1 → O_2,2 → O_2,3→ E_G, and the makespan is 17. The critical operation set is PO_l = {O_1,1, O_2,1, O_2,2, O_2,3}, and all the public factors for the critical operations are 1. So we can try to move the critical operations O_1,1, O_2,1, O_2,2, and O_2,3 successively. In this process, the algorithm preferentially selects the machine with the least processing time is and then judges the feasibility according to Eq. (6). Fig. 3 (b) show the disjunctive graph G₁ obtained by moving the critical operation O_1,1 to machine M₁. In this case, the makespan is reduced to 12. Fig. 3 (c) show the disjunctive graph G₂ obtained by moving the critical operation O_2,1 to machine M₃. IN this case, the makespan is 11. Following this, the critical operation O_2,2 is not satisfied the moving condition according to Eq. (6). Finally, the algorithm obtains the disjunctive graph Fig. 3 (d) by moving the critical operation O_2,3 to machine M₄. In this case, the makespan is reduced to 9.

Fig. 3.

Illustration of the local search

4.1. Representation and movement

In AFSA-ED, each AF represents a feasible solution of the problem. Each AF is expressed by two vectors: machine assignment vector and operation sequence vector, which correspond to the two sub-problems of the FJSP. The machine assignment vector is represented by a vector of N integer values and N is the total number of operations. Each element of vector denotes the machine selection of each operation and the value is the index of the array of alternative machine set. The operation sequence vector is an un-partitioned permutation with repetitions of job J_i (i = 1,2,...,n). The length of operation sequence vector equals to N. The index i of job J_i occurs n_i times in the vector, and the k-th occurrence of a job number refers to the k-th operation in the technological sequence of this job.

For the problem in Table 1, a representation of a feasible solution is shown in Fig. 4. The machine assignment vector is (1,2,2,2,3). If the operation O_2,1 is processed on {M₂, M₃, M₄}, then the corresponding element ‘2’ means that operation O_2,1 will be assigned to the second machine M₃. If the operation sequence is given as (2,1,2,2,1), then scanning the vector from left to right, the processing order of operations can be obtained: O_2,1 − O_1,1 − O_2,2 − O_2,3 − O_1,2.

Fig. 4.

Illustration of the representation of a solution

To address the discrete FJSP, the movement of an AF in the solution space is completed by learning the partial structure from the target AF. Fig. 5 gives an illustrative example of the movement process. Assume the i-th AF is (1,2,2,2,3;2,1,2,2,1) in the current generation, and its target AF is (2,2,1,3,2;1,2,2,1,2), then an indicator vector with the same length is produced by randomly filled with the elements of the set {0, 1}. The number of the element “1” is int[N2] in the first part of the indicator vector, whereas the second part of the indicator vector has only one position filled by “1”. The elements in the first part of the indicator vector decide the sources of the corresponding elements in new produced AF, namely from the parent or the target. The elements “1” in the second part of the indicator vector represents the starting position of the operation sequence segment from the target AF, and its length is also taken as int[N2], and the remainder is from the parent AF. Finally, m machine selections and operation orders are adjusted randomly according to the corresponding update formula, and m = int[step · rand()] as in the formula (2).

Fig. 5.

Illustration of the movement of an AF

4.2. Decoding of AF

For calculating the value of the makespan, each individual in the artificial fish swarm is decoded to the corresponding schedule sequence. The decoding is achieved by the process that assigns operations to the machines at their earliest possible starting time according to technological order of the jobs. It is worth noticing that the scheduling acquired in this way is semi-active. Then the active decoding is applied, which checks the possible blank time interval before appending an operation at the last position, and fills the first blank interval before the last operation to convert the semi-active schedule to an active one so that the makespan can be shorten³⁴.

4.3. Initialization

In this subsection, an integrated initialization algorithm is proposed for machine assignment initialization and operation sequence initialization. To generate the initial machine assignments, the following rules are applied:

1. (1)

   Global approach of localization (GAL)¹⁷.

2. (2)

   Local approach of localization (LAL).

3. (3)

   Random rule.

The LAL makes a random permutation for the positions of machines. Following this, for each operation, it selects the machine with minimum processing time in the alternative machine set, and updates the machine workload by adding this processing to the processing time of the remaining unarranged operations within the same job. Take the problem showed in Table 1 as an example, Table 2 shows a possible machine assignment obtained by using the LAL. The last four columns indicate the final assignments obtained by LAL. In the table, the items in bold type are the updated workload of machines, and the “-” means that the operation cannot be processed on the corresponding machine.

	M3	M2	M4	M1	M3	M2	M4	M1	M3	M2	M4	M1	...	M3	M2	M4	M1
O1,1	-	5	6	3	-	5	6		-	5	6		…	-	5	6
O1,2	4	-	5	6	4	-	5	9		-	5	9	…		-	5	9
O2,1	2	5	3	-	2	5	3	-	2	5	3	-	…		5	3	-
O2,2	5	1	3	1	5	1	3	1	5	1	3	1		7		3	1
O2,3	-	3	2	2	-	3	2	2	-	3	2	2		-	4		2

Table 2.

Initial assignments by LAL

The random rule executes by randomly selecting a machine from the alternative machine set for each operation. The GAL emphasizes the global workload among all the machines. The advantage of the LAL is that it obtains different initial assignments in different runs of the algorithm and emphasizes the workload among the set of machines within the same job. In addition, the random rule can increase the diversity of initial population. In our algorithm, the above three rules are used in a hybrid way. More specially, the initial machine assignments of 30% solutions in the population are generated by the GAL, 50% solutions by the LAL, and 20% solutions by random rule.

In our algorithm, the initial operation sequences are generated by the following three dispatching rules:

1. (1)

   Most time remaining (MTR). The job with the most remaining processing time will be arranged first.

2. (2)

   Most number of operations remaining rule (MOR). The job with the most remaining unprocessed operations has a high priority to be arranged.

3. (3)

   Random rule. Randomly generate the sequence of the operations. In particular, the above three rules are used in a hybrid way, then 20% of initial operation sequences are generated by random rule, 40% by the MTR, and 40% by the MOR.

4.4. The framework of AFSA-ED

The framework of AFSA-ED for solving FJSP is shown in Fig. 6. At the beginning of each generation, the entire population is randomly divided into two equal size sub-populations with one sub-population using machine based mechanism, and the other sub-population using operation based mechanism. Each AF evaluates by executing swarming, following, attracting behaviors. If the algorithm obtains a better solution using the three behaviors, it selects a better behavior to execute, otherwise, it executes the default improved preying behavior. When all individuals complete the searching process, the two sub-populations are combined into an entire population. Then the optimal individual of the population executes the critical path based on local search for further exploitation. If a new better individual is obtained, then the bulletin board is updated accordingly. The algorithm stops when the maximum iteration time is reached.

Fig. 6.

The framework of the AFSA-ED for the FJSP

5.1. Instances and parameters

To evaluate the performance of the AFSA-ED, we consider three sets of well-known benchmarks with 160 instances:

1. (1)

   BRdata: The data set includes 10 instances from Brandimarte¹¹. The number of jobs ranges from 10 to 20, the number of machines ranges from 4 to 15, and the flexibility of each operation ranges from 1.43 to 4.10.

2. (2)

   BCdata: The data set includes 21 instances from Barnes and Chambers³⁵, which were acquired from the classical JSP mt10 and the la24, la40 instances. The number of jobs ranges from 10 to 15, the number of machines ranges from 11 to 18, and the flexibility of each operation ranges from 1.07 to 1.30.

3. (3)

   HUdata: The data set includes 129 instances from Hurink et al. ³⁶, which were obtained from 3 classical JSP instances (mt06, mt10, mt20) by Fisher and Thompson and 40 classical JSP instances (la01–la40) by Lawrence. HUdata is divided into three subsets: Edata, Rdata, and Vdata. The number of jobs ranges from 6 to 30, the number of machines ranges from 5 to 15, and the flexibility of each operation ranges from 1.15 to 7.5.

The AFSA-ED is coded and implemented in matlab language on an Intel Core i5 2.53 GHz personal computer with 1GB of RAM. The algorithm runs 30 independent times for each instance from BRdata and BCdata, and runs 10 independent times for each instance from HUdata on account of the large number of instances in this data set. The computational results are compared with several performing algorithms from the existing literatures.

Each instance can be characterized by the following parameters: number of jobs (n), number of machines (m), number of operations (T₀) and the flexibility of problem (Flex). The parameters in the AFSA-ED include population size (AFs), the maximum iteration times, the try number, the step of AF (Step), the visual of AF (Visual), and the crowd factor (delta). In our experiment, the iteration times and the try number are taken as 40. The settings of the other parameters for each data set are listed in Table 3.

5.2. Computational results

The computational results for each data set are shown in this subsection. In the following tables, (LB,UB) denotes the lower and upper bounds³⁷. The LB of BRdata and BCdata instances are taken from Mastrolilli¹², while the LB of HUdata instances are computed by Jurisch.³⁸ C_m denotes the best makespan. AV denotes the average makespan. SD denotes standard deviation of the makespan. T_a is the average running time of the algorithm in terms of seconds. To illustrate the quality of the results obtained by the AFSA and the compared algorithms, the mean relative error (MRE) is also introduced. The relative error (RE) is calculated as follows: RE = (MS − LB)/LB × 100%, where MS is the makespan obtained by the corresponding algorithm.

5.2.1. Influence of proposed EDA process

To investigate the influence of proposed EDA process for preying behavior, the experiments on the BRdata are conducted by implementing the AFSA-ED and the AFSA without EDA, respectively. The parameters in the AFSA are taken as the same in AFSA-ED, and the two algorithms runs 30 independent times for each instance. The results are given in Table 4. Compared with AFSA, AFSA-ED outperforms AFSA in all 10 instances for the average makespan and the best makespan, and the SD values obtained by AFSA-ED are relative smaller than AFSA for most of the instances. However, the overall computation time of the AFSA-ED is slightly longer than AFSA because of the EDA process. Fig. 7 and Fig. 8 show the convergence curves in solving MK7 and MK10 by AFSA-ED and AFSA respectively. The experiment results show the validity of the proposed EDA process for preying behavior.

Instance	n × m	T0	Flex.	(LB,UB)	ABC			HHS			AFSA				AFSA-ED
					Cm	AV	SD	Cm	AV	SD	Cm	AV	SD	Ta	Cm	AV	SD	Ta
Mk01	10×6	55	2.09	(36,42)	40	40.00	0.00	40	40.00	0.00	40	41.25	0.65	1.39	40	40.00	0.00	2.59
Mk02	10×6	58	4.01	(24,32)	26	26.50	0.50	26	26.63	0.49	28	29.04	0.58	1.61	26	26.10	0.30	2.82
Mk03	15×8	150	3.01	(204,211)	204	204.00	0.00	204	204.00	0.00	204	204.00	0.00	1.00	204	204.00	0.00	1.12
Mk04	15×8	90	1.91	(48,81)	60	61.22	1.36	60	60.03	1.18	64	65.51	1.03	5.68	60	60.27	0.72	12.94
Mk05	15×4	106	1.71	(168,186)	172	172.98	0.14	172	172.80	0.41	177	177.69	0.47	4.75	172	172.40	0.48	10.49
Mk06	10×15	150	3.27	(33,86)	60	64.48	1.75	58	59.13	0.63	62	63.86	0.72	20.58	57	58.60	0.75	39.27
Mk07	20×5	100	2.83	(133,157)	139	141.42	1.20	139	139.57	0.50	142	143.03	0.95	28.47	139	140.63	0.71	57.63
Mk08	20×10	225	1.43	523	523	523.00	0.00	523	523.00	0.00	523	523.00	0.00	1.04	523	523.00	0.00	2.14
Mk09	20×10	240	2.53	(299,369)	307	308.76	1.63	307	307.00	0.00	310	310.75	0.54	4.06	307	307.13	0.43	9.40
Mk10	20×15	240	2.98	(165,296)	208	212.84	2.43	205	211.13	2.37	213	214.91	1.37	50.44	201	201.93	1.06	104.10

Table 4.

Results of ten BRdata instances

Fig. 7.

Convergence curves of MK07

Fig. 8.

Convergence curves of MK10

5.2.2. Results of BRdata problems

The AFSA-ED is first tested on ten instances of BRdata. Meanwhile, we make a comparison with ABC algorithm by Wang et al.³⁹ and HHS algorithm by Yuan et al.⁴⁰. These results are also given in Table 4. It can be seen that AFSA-ED, ABC and HHS obtain the same best result for instances Mk01-Mk05, and Mk07-Mk09. For Mk06 and Mk10, AFSA-ED obtains the values of 57 and 201, respectively. On the other hand, the ABC obtains the best values of 60 and 208, respectively, and HHS obtains the best values of 58 and 205, respectively. Compared with ABC, AFSA-ED outperforms ABC in all 10 instances for the average makespan, and the SD values of AFSA-ED are relative smaller than ABC. In comparison with HHS, AFSA-ED outperforms HHS in 7 out 10 instances for the average makespan. It is worth noting that AFSA-ED can obtain average value of 201.93 and SD value of 1.06 for instance Mk10, respectively. On the other hand, the HHS can obtain average value of 211.13 and SD value of 2.37, respectively.

In addition, Table 5 gives a detailed comparison in terms of the MRE of the best value and the MRE of average value. We compare AFSA-ED with the PVNS of Yazdani et al.⁴¹, the CDDS of Ben Hmida et al.³⁷, the BEDA of Wang et al.⁴², the ABC and the HHS. BKS represents the best known solution ever reported in the literature for each instance. It can be seen that AFSA-ED finds 9 best known solutions for 10 instances. The AFSA-ED obtains the MRE of the best value which is equal to 14.85%, while PVNS, CDDS, ABC, BDEA, and HHS are 16.43%, 14.98%, 16.19%, 16.07%, and 15.40%, respectively. For the MRE of average value obtained, AFSA-ED generates 15.64%, faced to 16.48% for HHS, 17.01% for PVNS, 18.55% for ABC, and 19.24 % for BDEA. However, faced to 15.34 % for the CDDS algorithm, the AFSA-ED Obtains a higher MRE. The Gantt chart of solution of MK06 obtained by AFSA-ED is showed in Fig. 9.

Instance	BKS	AFSA-ED		PVNS		CDDS
		Cm(AV)	MRE	Cm(AV)	MRE	Cm(AV)	MRE
Mk01	40	40 (40)	11.11(11.11)	40 (40)	11.11(11.11)	40 (40)	11.11(11.11)
Mk02	26	26(26.10)	8.33(8.75)	26(26.04)	8.33(8.50)	26(26)	8.33(8.33)
Mk03	204	204(204)	0.00(0.00)	204(204)	0.00(0.00)	204(204)	0.00(0.00)
Mk04	60	60(60.27)	25.00(25.56)	60(60.60)	25.00(26.25)	60(60)	25.00(25.00)
Mk05	172	172(172.40)	2.38(2.62)	173(173)	2.98(2.98)	173(173.5)	2.98(3.27)
Mk06	57	57(58.60)	72.73(77.57)	60(61)	81.82(84.85)	58(59)	75.76(78.79)
Mk07	139	139(140.63)	4.51(5.73)	141(141.2)	6.02(6.17)	139(139)	4.51(4.51)
Mk08	523	523(523)	0.00(0.00)	523(523)	0.00(0.00)	523(523)	0.00(0.00)
Mk09	307	307(307.13)	2.68(2.72)	308(308.8)	3.01(3.28)	307(307)	2.68(2.68)
Mk10	197	201(201.93)	21.82(22.38)	208(209.4)	26.06(26.91)	197(197.75)	19.39(19.85)
MRE		14.85(15.64)		16.43(17.01)		14.98(15.34)

Instance	BKS	ABC		BEDA		HHS
		Cm(AV)	MRE	Cm(AV)	MRE	Cm(AV)	MRE
Mk01	40	40 (40)	11.11(11.11)	40 (41.02)	11.11(13.94)	40 (40)	11.11(11.11)
Mk02	26	26(26.50)	8.33(10.24)	26(27.25)	8.33(13.54)	26(26.63)	8.33(10.96)
Mk03	20	204(204)	0.00(0.00)	204(204)	0.00(0.00)	204(204)	0.00(0.00)
Mk04	60	60(61.22)	25.00(27.54)	60(63.69)	25.00(32.69)	60(60.03)	25.00(25.06)
Mk05	172	172(172.98)	2.38(2.96)	173(173.38)	2.98(3.20)	172(172.80)	2.38(2.86)
Mk06	57	60(64.48)	81.82(95.39)	60(62.83)	81.82(90.39)	58(59.13)	75.76(79.18)
Mk07	139	139(141.42)	4.51(6.33)	139(141.55)	4.51(6.43)	139(139.57)	4.51(4.94)
Mk08	523	523(523)	0.00(0.00)	523(523)	0.00(0.00)	523(523)	0.00(0.00)
Mk09	307	307(308.76)	2.68(2.93)	307(310.35)	2.68(3.80)	307(307)	2.68(2.68)
Mk10	197	208(212.84)	26.06(28.99)	206(211.92)	24.85(28.44)	205(211.13)	24.24(27.96)
MRE		16.19(18.55)		16.07(19.24)		15.40(16.48)

Table 5.

Comparison between AFSA-ED and several existing algorithms on BRdata

Fig. 9.

Gantt chart of solution of MK06 obtained by AFSA-ED (makespan = 57)

To determine the statistical differences between the AFSA-ED and the compared algorithms, the Friedman test and Holm procedure are conducted. The results are presented in Table 6. It can be seen from the Friedman test results that the differences among the six algorithms are statistically relevant with 98% certainty. The AFSA-ED obtains the best overall rank. The holm procedure shows that the AFSA-ED obtains better results than the compared five algorithms, and the differences are statistically relevant with 97%, 66%, 81%, 86%, and 67% certainty, respectively.

Friedman test					Holm procedure
Algorithm	Rank	χ2	1 − p value	Diff.?	Algorithm	z	1 − p value
AFSA-ED	2.85	14.42	0.98	Yes	-	-	-
PVNS	4.50				AFSA-ED v.s. PVNS	1.97	0.98
CDDS	3.20				AFSA-ED v.s. CDDS	0.41	0.66
ABC	3.60				AFSA-ED v.s. ABC	0.89	0.81
BEDA	3.75				AFSA-ED v.s. BEDA	1.07	0.86
HHS	3.10				AFSA-ED v.s. HHS	0.29	0.67

Table 6.

Friedman test and Holm procedure of different algorithms

5.3. Discussions

In Table 11, we summarize the MRE of the best makespan obtained by our algorithm and other known algorithm for BRdata set, BCdata set and HUdata set. The first column shows the name of data set, the second column shows the number of problems for each set, the following nine columns show the MRE of AFSA-ED, GA of Chen¹⁵, TS of Mastrolilli ¹², GA of Pezzella¹⁷, CDDS, PVNS, HHS, HDE-N₁ and HDE-N₂. As can be seen in Table 11, HDE-N₂ achieves the state-of-the-art performance on three data sets. Except for HDE-N₂, AFSA-ED outperforms the other seven algorithms on BRdata and BCdata; it outperforms GA_Chen, GA_Pezzella, PVNS and HHS on three sub-data sets of HUdata, and CDDS on Edata and Rdata. TS works better than AFSA-ED on the HUdata, and CDDS works better than AFSA-ED on Vdata. In particular, HDE-N₂ is a time consuming algorithm, in which the computational time is about 5 times longer than AFSA-ED as shown in Table 7. In conclusion, the results show that the AFSA-ED is an effective algorithm for FJSP.