A hybrid fluid master–apprentice evolutionary algorithm for large-scale multiplicity flexible job-shop scheduling with sequence-dependent set-up time

In this article, a large-scale multiplicity flexible job-shop scheduling problem (FJSP) with sequence-dependent set-up time is studied. In this problem, the large production demand for each type of job yields the large-scale multiplicity manufacturing feature. To address the problem, a hybrid fluid master–apprentice evolutionary algorithm (HFMAE) is presented to minimize the makespan. In the first step, a fluid relaxation initialization method (FRI) and an initialize procedure are proposed to obtain high-quality initial solutions. In the FRI, an online fluid tracking policy is presented to improve the assignment decision and the sequencing decision of operations. In the second step, an improved master–apprentice evolutionary method (IMAE) is presented based on the generated initial solutions. In the IMAE, two neighbourhood structures and three makespan estimation approaches are presented to accelerate the solution space search efficiency. Numerical results show that the proposed HFMAE outperforms the comparison algorithms in solving large-scale multiplicity FJSPs.


Introduction
As a classical scheduling problem, the job-shop scheduling problem is an NP-hard problem that determines the operation sequence on each machine (Goncalves, Mendes, and Resende 2005).The flexible job-shop scheduling problem (FJSP) considers both the operation sequence and machine assignment problems (Chaudhry and Khan 2016).Therefore, the FJSP is a more complex NP-hard problem than the job-shop scheduling problem (Pezzella, Morganti, and Ciaschetti 2008).
In this article, the large-scale multiplicity FJSP is a variant of the regular FJSP that allows multiple jobs for each job type.It contains I machines, indexed as i = 1, 2, 3 . . .I, R job types, indexed as r = 1, 2, 3 . . .R, and N r jobs (multiplicity) of each job type r, indexed as n = 1, 2, . . ., N r .In the problem, jobs are divided into R different types and the jobs belonging to the same type are identical.The number of jobs N r of a specific type r is called multiplicity (N r ≥ 1).Differing from the research problem, regular FJSPs consist of a single job per job type (N r ≡ 1).The total number of jobs N = R r=1 N r is a large value (N × I ≥ 1000), which leads to a large-scale scheduling problem (Sun et al. 2019).
© 2022 Informa UK Limited, trading as Taylor & Francis Group The number of jobs is large in manufacturing problems, which leads to a large exponential solution space and unacceptable computing time (Mandavi, Shiri, and Rahnamayan 2015).There are many large-scale multiplicity FJSPs in real production systems, such as motorcycle parts manufacturing lines and semiconductor manufacturing lines (Monch et al. 2011).The set-up time is the interval needed to adjust the settings on a machine to process the next products.The set-up time is mainly related to the operations of changing tools, fixing or releasing parts to machines, and adjusting the tools to the machines (Pan et al. 2017).It is significant to consider sequence-dependent set-up time as it exists in some real manufacturing lines for certain products, such as in automobile and furniture manufacturing, and the paint industry (Ruiz and Stutzle 2008;Cakmakci 2009;Gawronski 2012).However, no relevant literature was found that considered both sequence-dependent set-up time and large-scale multiplicity characteristics.Therefore, to fill the gap in the literature, this study considers the inclusion of sequence-dependent set-up time for FJSP with the large-scale multiplicity characteristic.To address this problem efficiently, a hybrid fluid master-apprentice evolutionary algorithm (HFMAE) is presented.In the HFMAE, an online fluid tracking policy is designed based on the proposed fluid model.Then, a fluid relaxation initialization method (FRI) and an initialize procedure are proposed to obtain high-quality initial solutions.Finally, two neighbourhood structures and three makespan estimation approaches are presented to accelerate the solution space search efficiency.
The main contributions of this article can be summarized as follows: (1) A large-scale multiplicity FJSP with sequence-dependent set-up time is presented and the mathematical model is established.(2) A fluid model and an online fluid tracking policy are proposed to obtain a high-quality initial solution.
(3) Two neighbourhood structures and three makespan estimation approaches are presented to accelerate the solution space search efficiency.(4) The HFMAE is proposed for large-scale multiplicity FJSP with sequence-dependent set-up time, and numerical examples are applied to verify the superiority of the proposed algorithm.
The remainder of this article is organized as follows.In Section 2, the related literature is reviewed.In Section 3, the large-scale multiplicity FJSP with sequence-dependent set-up time is described in detail.In Section 4, the HFMAE is proposed to solve the proposed problem.In Section 5, the performance of the proposed HFMAE on cases of different scale and multiplicity is verified.Conclusions are presented in Section 6.

Literature review
In this section, the literature related to the fluid approach for large-scale multiplicity job-shop scheduling problems, and the non-decomposition and decomposition approaches for large-scale FJSPs, is reviewed.For more details on the recent research on the FJSP, readers can refer to Cao, Shi, and Chang (2022).

The fluid approach
In recent years, multiple algorithms based on the fluid approach have been proposed to solve the large-scale multiplicity job-shop scheduling problem.Bertsimas and Gamarnik (1999) proposed a fluid relaxation algorithm to solve packet routing and large-scale multiplicity job-shop scheduling problems.Boudoukh, Penn, and Weiss (2001) proposed an approximation algorithm that approximates the large-scale multiplicity job-shop scheduling problem with some identical or similar jobs as a continuous and deterministic scheduling problem.Dai and Weiss (2002) presented an online scheduling algorithm based on the fluid optimal solution that kept the bottleneck machine busy by controlling its safety inventory.Nazarathy and Weiss (2010) presented a fluid approximation method that was asymptotically optimal with the increase in the number of different job types.
Many scholars have approximated the fluid optimal solution by constructing formulae and proposed some fluid approaches to solve large-scale multiplicity job-shop scheduling problems.Bertsimas, Gamarnik, and Sethuraman (2003) combined the fluid relaxation approach and fair queuing in communication networks to construct the fluid tracking formula.Gu, Lu, and Gu (2017) presented a tracking virtual algorithm to solve this problem by constructing a virtual fluid formula.Gu et al. (2018) proposed a policy to address the large-scale multiplicity job-shop scheduling problem with fuzzy processing time.Ding, Guan, and Zhang (2021) presented a fluid relaxation algorithm for largescale re-entrant flexible job-shop scheduling.However, the algorithm is rule based, and cannot obtain high-quality solutions for large-scale multiplicity FJSPs, especially when the sequence-dependent set-up time is considered.
In conclusion, many scholars have conducted a large amount of research on large-scale multiplicity job-shop scheduling based on fluid approaches.However, little research was found on large-scale multiplicity FJSPs with fluid approaches.

Non-decomposition and decomposition approaches
Scholars have proposed many approaches to solve large-scale FJSPs, mainly including nondecomposition and decomposition approaches.For non-decomposition approaches, researchers have usually developed improved heuristic algorithms to obtain a high-quality solution within a limited computing time.Sobeyko and Monch (2016) proposed an efficient iterative local search approach and a simulated annealing criterion for large-scale FJSPs with identical and unrelated parallel machines.Sobeyko and Monch (2017) designed an iterative scheme based on the variable neighbourhood search and proposed appropriate neighbourhood structures for the variable neighbourhood search to integrate process planning and scheduling.Yang et al. (2022) proposed an algorithm with a dynamic opposite learning strategy to solve the large-scale FJSP.
Moreover, some non-decomposition approaches have been proposed for solving scheduling problems with sequence-dependent set-up time.Rifai, Windras Mara, and Sudiarso (2021) proposed an improved version of the multi-objective adaptive large neighbourhood search for the sequencedependent distributed re-entrant permutation flow shop.Alimian et al. (2022) developed a rolling horizon heuristic algorithm and a novel mixed-integer nonlinear programming model for a parallelline capacitated lot-sizing problem with a sequence-dependent set-up time.Rifai, Windras Mara, and Norcahyo (2022) presented a two-stage heuristic procedure to solve the sequence-dependent job sequencing and tool switching problem.Rani, Jain, and Angra (2022) developed a simulation model to assess the effect of routing flexibility on order release policies in a flexible job shop considering sequence-dependent set-up time, within a stochastic and dynamic manufacturing environment.
However, the non-decomposition approach may not solve large-scale FJSPs with sequencedependent set-up time effectively as the scale of the problem increases (Mandavi, Shiri, and Rahnamayan 2015).Thus, many decomposition approaches have been presented to improve the solving efficiency.Decomposition approaches mainly include machine-based approaches, task-based approaches and grouping-based approaches (Sun et al. 2019).Zhai et al. (2014) constructed many subproblems by decomposing the large-scale scheduling problem.Thus, a large-scale problem can be addressed by solving all of the decomposed subproblems.Borisovsky, Eremeev, and Kallrath (2020) proposed a hybrid approach that combined a constructive mixed-integer linear programming-based heuristic, a genetic algorithm and a decomposition approach.In addition, some decomposition algorithms based on cooperative evolutionary methods and the distributed model were proposed to solve large-scale FJSPs (Sun et al. 2019;Meselhi et al. 2020).
In conclusion, most of these non-decomposition and decomposition approaches were developed based on the population.Population-based approaches need to manage a large population to maintain the diversity of solutions, which is time consuming.Ding et al. (2019) proposed a two-individualbased evolutionary algorithm named the master-apprentice evolutionary algorithm (MAE) for solving FJSPs.In this article, an improved master-apprentice evolutionary method (IMAE) is presented based on the MAE framework.In the IMAE, two neighbourhood structures and three makespan estimation methods are presented to accelerate the local search phase.The proposed IMAE needs two initial solutions.Therefore, an initializing procedure based on the encoding-decoding mechanism (Defersha and Rooyani 2020) is proposed to obtain another high-quality initial solution.Ultimately, the HFMAE for the large-scale multiplicity FJSP with sequence-dependent set-up time is proposed by combining the FRI and IMAE.

Problem description
Up to now, few studies on large-scale FJSP with sequence-dependent set-up time or other characteristics have been found.Defersha and Rooyani (2020) first proposed an efficient two-stage genetic algorithm to solve the large-scale FJSP with sequence-dependent set-up time.In this section, the large-scale multiplicity FJSP with sequence-dependent detached set-up time is described.
Owing to the complexity of the large-scale multiplicity FJSP with sequence-dependent set-up time, the following assumptions are presented: (1) All machines are available at time zero.
(2) Each operation can only be processed on one machine.
(3) Each machine can only process one operation at a time.
(4) Once one operation has been processed on a machine, it must not be interrupted.
(5) An operation cannot be processed until its previous operation is completed.
The notation used in this method is described as follows.Binary variable, 1 if o r n j is the tight post-operation of o rnj processed on the same machine and 0 otherwise

Notation
The mathematical model is described as follows: subject to: st r n j ≥ ct rnj + S rjir j , z rnjr n j = 1, x rnji = 1 (7) In the mathematical model, Equation (1) represents the objective function.Constraint (2) means that the makespan is equal to the completion time of the final operation.Constraint (3) indicates that each operation can only be processed by one machine from the set of available machines.Constraint (4) is applied to determine the start time of the first operation processed on each machine.The completion time of an operation is calculated as constraint (5).Constraint (6) ensures the precedence constraints of operations belonging to the same job.Constraint (7) means that each machine can process only one operation at each time and cannot process any operation during the set-up time.

Hybrid fluid master-apprentice evolutionary algorithm
In the proposed HFMAE, the FRI is presented to obtain a high-quality initial solution.Then, another high-quality initial solution is obtained with a proposed initialize procedure.Finally, the global best solution to this problem is obtained by improving the generated two initial solutions with the IMAE.The flowchart of the HFMAE is shown in Figure 1.

Proposed fluid model
In the fluid model, the different operations of multiple job types are divided into different classes, indexed as k = 1, 2, 3 . . .K. Each operation j of job type r can only be represented by a unique class k.In this model, Z krj indicates whether class k represents the operation o rj .Moreover, jobs are composed of fluid and discrete jobs can be represented as continuous fluid flows.Therefore, the number of jobs can be represented as the number of fluids, which is not required to be an integer.The processing time of each machine can be allocated to its available classes and each machine has a fixed proportion of times (real number between 0 and 1) to process different available classes within a unit of time.Based on the fluid model, the proportion of processing times of all machines allocated to their available classes within a unit of time can be obtained.

Model detail
The notation used in the fluid model is introduced as follows.The fluid model is described as follows.In the fluid model, the fluid number and the actual number of each class satisfy Equations ( 8)-( 14).Equations ( 8) and ( 9) mean the fluid number of each class at the initial time and time t.Equations ( 10) and ( 11) provide the actual total number of each class that has not been processed at the initial time and time t.Equations ( 12) and ( 13) mean the total fluid number of each class that has not completed processing by machine i at the initial time and time t.Equation ( 14) means the completion time of each class in the fluid model.

Notation
The objective is to minimize the maximum completion time of all classes, which is calculated with Equation ( 15): Equations ( 16)-( 19) provide the constraints of the fluid model.Constraint ( 16) means that each operation of each type of job can only be represented by a unique class.Constraint (17) shows that machine utilization is less than 100%.Constraint (18) represents the range of decision variables.Constraint (19) indicates that the processing rate is no more than the arrival rate if the initial fluid number of class k is equal to zero, which guarantees that the solution is feasible.

Translating the fluid model solution into decision policy
Based on the proposed fluid model, an online fluid tracking policy (Nazarathy and Weiss 2010) is presented to guide the job selection and machine assignment in detail.In the online fluid tracking policy, p k (t) represents the priority of class k at time t, which is calculated by Equation ( 20).The priority of each available machine of the selected class k is calculated with Equation ( 21).At each decision phase, a class k with the maximum p k (t) value is selected (if the maximum value is not unique, an arbitrary tie-breaking rule is used to select a class).From the set of available machines, the machine that has the maximum p ik (t) value is selected to process the selected class k.
Finally, the class selection and machine assignment at each decision phase can be finished with the online fluid tracking policy.

Fluid relaxation initialization method
In this subsection, the FRI is presented to obtain a high-quality initial solution based on the online fluid tracking policy.As shown in Algorithm 1, FRI first selects the class k with maximum p k (t) value, and then selects the machine i with maximum p ik (t) value to process the selected class k .Finally, all classes are processed and the initial solution S 1 of the HFMAE is obtained with FRI.

Initialize procedure
In this subsection, an initialize procedure is presented to obtain another initial solution for the HFMAE.In the initialize procedure, D operation sequences are represented as {w 1 , w 2 , . . .w D } are randomly generated.In the {w 1 , w 2 , . . .w D }, each operation of each operation sequence is represented as (r, n, j) which has not been assigned a machine.Then, an available machine is assigned to each operation of each operation sequence in the {w 1 , w 2 , . . .w D } with the encoding-decoding mechanism (Defersha and Rooyani 2020) and D operation sequences (feasible solutions) are obtained {w 1 , w 2 , . . .w D }. Finally, the best solution from {w 1 , w 2 , . . .w D } is selected as the initial solution S 2 .The initialize procedure is shown in Algorithm 2.

Improved master-apprentice evolutionary method
In the improved master-apprentice evolutionary method (IMAE), S 1 and S 2 are two initial solutions which are obtained from the FRI and the proposed initialize procedure.Then, S 1 and S 2 are improved with the proposed local search procedure.Furthermore, the path relinking operator (PR) (Ding et al. 2019) is adopted to improve the global search ability of IMAE.The detail of the IMAE is shown in Algorithm 3.
Algorithm 1: The procedure of the fluid relaxation initialization method

Input:
The problem data parameters.Obtain all u ik values by solving the fluid model with CPLEX.

Initialization:
The current time time now ← 0. The next idle time of each machine time next (i) ← +∞, ∀i ∈ M I .The state of each machine state(i) = 0, ∀i ∈ M I (all machines are idling at the initial time).while Update the machine set M now (machines are idling and have processable classes at time now ): ∀i ∈ M I , if state(i) = 0 and Update time now , time next (i), and state(i): The solution of FRI.
In Algorithm 3, g cycle denotes the cycle length of the iterations.It is assumed that (r, n, j) is the n 1 th operation processed on machine i 1 in S 1 and the n 2 th operation processed on machine i 2 in S 2 .
If n 1 = n 2 and i 1 = i 2 , the operation (r, n, j) is in the same position in S 1 and S 2 .If the total number of operations in the same position in S 1 and S 2 is more than 0.9 × J o (J o is the total number of operations), S 1 and S 2 are similar (S 1 ≈ S 2 ).Then, S 2 is obtained with the proposed initialize procedure.

Local search by tabu search
The tabu search (TS) algorithm was first proposed by Glover (1990), and has been successfully applied to various scheduling problems.A TS algorithm is proposed as the local search strategy in the IMAE.In the proposed TS, the length of the tabu list changes dynamically.The TS terminates when the permitted maximum step size with no improvement (iter stop ).A tuple ((r, n, j), i, MP rnj , MS rnj , i , MP rnj , MS rnj ) is the tabu element in the tabu list, where PM rnj (PM rnj ) and SM rnj (SM rnj ) denote the predecessor and successor of operation (r, n, j) on the machine i before (after) the move.

Neighbourhood structures and makespan estimate approaches
In this subsection, two neighbourhood structures are proposed for the operation sequence subproblem, and the machine assignment change is performed in the neighbourhood called N α−ρ 1 , which was proposed by Kemmoé-Tchomté, Lamy, and Tchernev (2017) for the machine assignment subproblem.For accelerating the TS phase, three makespan estimation approaches are proposed to estimate the makespan of the neighbourhood solutions.
The FJSP can be represented by the disjunctive graph model (DauzerePeres and Paulli 1997).To clearly describe the neighbourhood structures, a disjunctive graph G = (V, A, E) is used to represent the FJSP with sequence-dependent set-up time.V denotes the set of nodes, in which each node Algorithm 3: The procedure of the improved master-apprentice evolutionary method

Initialization:
The number of iterations g: g ← 0.
Initialize solutions S 1 and S 2 , the best solution in the current cycle S * c , the best solution in the previous cycle S * p , and the best solution S * : S 1 ← FRI, S 2 , S * c , S * p , S * ← initialize_procedure(S 1 ).Initialize solutions S 1 and S 2 : end while Output: The best solution S * .
represents an operation or one of the start node and end node.A denotes the set of conjunctive arcs (oriented).The conjunctive arcs represent processing sequence constraints of the operations which belong to the same job.Furthermore, the connection arcs between the first operation and the start node (the last operation and the end node) are also represented by the conjunctive arcs.The weight of a conjunctive arc (v1, v2) is the processing time of the operation v1 on the assigned machine.E denotes the set of disjunctive arcs (non-oriented).The disjunctive arcs connect two operations that are processed continuously on the same machine and connect the start node with the first operations Algorithm 4: The procedure of the tabu search

Input:
The minimum iteration size with no improvement of TS iter min .The minimum length of the tabu list length max .The current number of iterations g of IMAE.The increment coefficient H of iter stop which is related to g.

Initialization:
The permitted maximum iteration size with no improvement of TS iter stop .
The  of each machine.The weight of a disjunctive arc (v1, v2) comprises the processing time of operation v1 and the set-up time of operation v1 and v2.
A solution corresponds to an acyclic subgraph in which each disjunctive arc is oriented.In the disjunctive graph model, a solution can be noted as S(α, π).In S(α, π), α denotes a vector where each operation is assigned to only one available machine (a feasible assignment) and π denotes the operation sequence on each machine.A feasible solution graph can be represented as G(α, π) = (V, A, E(α, π)).E(α, π) is the set of disjunctive arcs which represent the processing order of the operations on each machine.The makespan of G(α, π) is the cost of a critical path.The critical path has the maximum cost in the disjunctive graph model.A solution for a problem with three jobs and five machines is shown in Figure 2. In Figure 2, a critical path is shown in bold arcs.
Notation B i represents the sequence of consecutive operations assigned on machine i. R(v) denotes the longest path from the start node to operation v, and Q(v) denotes the longest path from operation v to the end node.For a critical operation v, C max = R(v) + t vi + Q(v).t vi denotes the processing time of operation v on the assigned machine i.In Figure 2 The proposed two neighbourhood structures are denoted as , which extends the neighbourhood structure N π (González, Vela, and Varela 2015).
Proposition 1: Let S(α, π) be a feasible schedule and let {B i , u, B i , v, B i } be a sequence of consecutive operations in the solution graph G(α, π).All operations in {B i , u, B i , v, B i } are assigned to the same machine i, ∀i ∈ M I .B i , B i and B i are consecutive operation sequences and any of them may be null.
(a) Let S(α, π ) be a schedule created from S(α, π) by moving u just after v.Then, G(α, π ) has no cycles if the following condition holds: if u is the last operation of the job (JS u = Ø), or JS u = ∅ and Q(v) ≥ Q(JS u ), v = JS u .(b) Let S(α, π ) be a schedule created from S(α, π) by moving v just before u.Then, G(α, π )

has no cycles if the following condition holds: if v is the first operation of the job (JP
A similar result can be obtained while moving v just before u.The proof of Proposition 1 is the same as Proposition 2 proposed by González, Vela, and Varela (2015).Then, the following neighbourhood structures are defined.
Definition 1: N π −B 1 structure.Let operation u be a critical operation which included in the operation sequence B i , B i = {B i , u, B i }.In a neighbouring solution, u is moved to another position in B i .The sufficient condition of feasibility that is given in Proposition 1(a) holds.
Figure 3(a) illustrates a new feasible solution after moving the operation u just after v with the neighbourhood structure N π−B 1 .Equation ( 22) estimates the makespan of the feasible solutions obtained by the neighbourhood structure N π−B 1 .
Definition 2: N π −B 2 structure.Let operation v be a critical operation that is included in the operation sequence B i , B i = {B i , v, B i }.In a neighbouring solution, v is moved to another position in B i .The sufficient condition of feasibility that is given in Proposition 1(b) must hold.
Figure 3 shows a neighbourhood solution with the neighbourhood structure N π −B

2
. The makespan of the neighbourhood solution is estimated by Equation ( 23).
Definition 3: N α−ρ 1 structure.Let S(α, π) be a schedule, and let j be a critical operation in any critical path ρ of G(α, π).A neighbouring solution is obtained by moving j into a position in B i (i is a new machine, i ∈ M j ).
Figure 3(c) illustrates the neighbourhood solution by moving a critical operation j into the operation sequence of a machine with the neighbourhood structure N α−ρ 1 .Equation ( 24) estimates the makespan of the neighbourhood.

Computational experiments and results
To verify the superiority of the proposed HFMAE, 10 generated benchmark data sets and different scale cases are adopted and tested.For all algorithms, the maximum central processing unit (CPU) time of n × m × t (ms) is set as the stopping criterion, where n is the total number of jobs that need to be scheduled, m is the number of machines, and t is a fixed value.The termination allows more computational time for cases on a larger scale.The parameter t is set as 200 in the experiments.
All algorithms are performed 30 times independently on each case and implemented on a personal computer with Intel ® Core TM I5-9300 CPU @ 2.40 GHz and 8 GB RAM and coded with Python 3.6.

Comparison with optimal value
In this subsection, the HFMAE was tested on generated small-scale instances (data01-data08).The optimal values are obtained with KNITRO 9.0.1 solver.Table 1 shows the results of the computation on the generated small instances.Within the limit of 2 h of computation time, the KNITRO 9.0.1 solver can only optimally solve the instance with less than or equal to {N = 9, I = 10}.The best solution obtained by the HFMAE in data06 and data07 outperforms KNITRO, which has a high gap.Therefore, the exact method is deemed to be ineffective for solving problems with a larger number of jobs and machines.

Experiments based on benchmarks
Brandimarte (1993) develops 10 benchmarks (BRdata) for the FJSP.However, the standard benchmarks consist of a single job per job type.Hence, to create large-scale benchmark data sets with multiplicity, each job is duplicated several times (N r ) based on the 10 benchmarks.In this subsection, the 10 generated large-scale benchmark data sets (LSM-Mk01-LSM-Mk10) are solved with the proposed HFMAE.Two classical evolutionary algorithms [tabu search algorithm (TS) (Mastrolilli and Gambardella 2000) and genetic algorithm (GA) (Pezzella, Morganti, and Ciaschetti 2008)] and two state-of-the-art algorithms [two-stage genetic algorithm (2SGA) (Defersha and Rooyani 2020) and master-apprentice evolutionary algorithm (MAE) (Ding et al. 2019)] are used to compare with the proposed HFMAE.Defersha and Rooyani (2020) proved that the 2SGA is efficient in solving large-scale FJSPs.It can obtain a solution in less than 30 min, which was not achieved after 72 h with the regular genetic algorithm.Therefore, the 2SGA is suitable as the comparison algorithm while solving the large-scale multiplicity FJSPs in this article.Table 2 shows the computational results of all the algorithms on the generated benchmark data sets.
In Table 2, the Best and avg columns represent the best solutions and the average solutions obtained by each algorithm.The t (s) column is the mean CPU time, which is obtained by recording the CPU time of all algorithms in each run.The best solutions obtained by each algorithm are shown in bold.The #best, #even and #worse rows represent that the number of instances obtained by the HFMAE is better than, equal to and worse than the reference algorithms.
From Table 2, the HFMAE obtains better results for the best and average makespan than the reference algorithms in most of the instances.The 2SGA is the closest to the HFMAE, but HFMAE is superior to 2SGA in eight out of the 10 instances, the exceptions being instances LSM-Mk03 and LSM-Mk04.Moreover, the HFMAE outperforms the four reference algorithms in seven out of the 10 instances in terms of CPU time.The reason for this is that the HFMAE can construct a high-quality initial solution with the proposed online fluid tracking policy.Furthermore, the proposed neighbourhood structures N π −B 1 and N π−B 2 extend the solution space and improve the possibility to find better solutions.The presented makespan estimate equations can quickly estimate the neighbourhood solutions and improve the solution space search efficiency.Therefore, it can be concluded that the HFMAE outperforms the algorithms TS, GA, MAE and 2SGA in solving the generated benchmark data sets.
Moreover, to check whether the proposed algorithm outperforms all reference algorithms statistically, the Wilcoxon test is used to compare the HFMAE with the reference algorithms (García et al. 2008).The Wilcoxon is a non-parametric test method that is suitable for continuous optimization problems in multiple-problem analysis.Table 3 summarizes the results of applying the Wilcoxon test.The table displays the sum of rankings obtained in each comparison and the associated p-value.The HFMAE outperforms these four algorithms considering independent pairwise comparisons since the p-values are below α = 0.05.However, Wilcoxon's test performs individual comparisons between two algorithms.If they are included within the multiple comparisons, the p-value obtained is p = 0.007797 in benchmark data sets according to Table 3.Then, it can be deduced that the HFMAE outperforms the four reference algorithms with a level of significance α = 0.05.

Experiments on generated data sets
In this subsection, some cases are constructed to verify the performance of the HFMAE in cases with different scales and multiplicity.Ten cases of different scales, i.e. small-and large-scale instances, are generated.These cases contain the same number of job types.Different numbers of jobs and machines lead to different problem scales and multiplicity.The processing time of all operations obeys the uniform random distribution U (10, 100).The set-up time obeys the uniform random distribution U (10, 50).The problem scale varies from 10 jobs with 30 operations processed on five machines to 200 jobs with 1000 operations processed on 30 machines.Table 4 shows the computational results.Note that the column Best is the best solution obtained by the HFMAE and reference algorithms.The column RPD is the average relative percentage deviation of the objective.Equation ( 25) defines the calculation of RPD.Here, C avg is the average makespan and C best is the best solution.The algorithm has a better performance when a smaller RPD is obtained.
The results in Table 4 show that the HFMAE obtains the best solution in nine out of the 10 cases.TS, GA and MAE only obtain the best solution to Small-1 problems.2SGA obtains the best solution for Small-1 and Small-2.Therefore, the HFMAE has a much better performance than GA, 2SGA and MAE on the best solution.It is concluded that the HFMAE obtains better solutions than the reference algorithms in all problems except for Small-1, and the improvement increases as the problem scale increases.Furthermore, 2SGA has superior performance among the reference algorithms.
Figure 4(a) shows the interval plot for the RPD of the HFMAE and reference algorithms, with 95% confidence intervals.From Figure 4(a), the HFMAE obtains better results, especially for large-scale problems.Moreover, the HFMAE has smaller fluctuations in results, which means that it has better stability and robustness.The mean CPU time of all algorithms in different scale cases is shown in Figure 4(b).The HFMAE takes less CPU time than the four reference algorithms to converge to the best solution as the problem size increases.This is due to the addition of FRI to HFMAE, which leads to a better initial solution.When S 1 and S 2 are similar, the initialize procedure generates a new highquality solution to ensure the diversity of the search.Furthermore, the HFMAE maintains a good balance between exploitation and exploration capabilities by choosing an appropriate parameter g cycle for different scale problems.
The Wilcoxon test is again used for checking whether the HFMAE outperforms the reference algorithms statistically.The HFMAE outperforms these four algorithms considering independent pairwise comparisons, as the p-values are below α = 0.05.Moreover, the p-value obtained is p = 0.011066 in the generated data sets when considering them within the multiple comparisons.Therefore, it can be concluded that the HFMAE outperforms the four reference algorithms with a level of significance α = 0.05.

Sensitivity analysis
To further assess the extent of the proposed algorithm, sensitivity analysis is performed on a largescale instance with {N = 50, I = 20} to analyse the impact of set-up time and processing time on the objective.Figure 5 shows the results of the sensitivity analysis.
From Figure 5(a) and (c), it can be observed that the changes in mean set-up time generally have a similar impact on the objective to the changes in mean processing time.Moreover, Figure 5(b) shows that as the variance of set-up time grows larger, the makespan tends to be reduced.Figure 5(d) shows that the changes in the variance of processing time have a similar impact on the objective.These   results demonstrate that the HFMAE can capture the influence of set-up time and processing time on the whole performance of the system.Furthermore, the proposed algorithm has good stability and robustness under various set-up time and processing time levels.

Industrial application
In this subsection, five manufacturing data sets of real-world motorcycle parts manufacturing plants are selected to verify the effectiveness of the proposed algorithm.The proposed algorithm, HFMAE, is compared with 2SGA, MAE and the following scheduling rules: shortest processing time (SPT), longest processing time (LPT) and first come, first served (FCFS).These rules are applied to the actual production scenario.The results are shown in Table 5.The results in Table 5 show that the HFMAE obtains the best solution in all five examples, which indicates that the HFMAE outperforms the three scheduling rules and the reference algorithms.As depicted in Figure 6, the HFMAE has the lowest RPD value, the best average makespan and the smallest interval, which means that the HFMAE obtains better results for all problems and has better stability.Therefore, it can be concluded that the proposed HFMAE can improve production efficiency in realistic manufacturing scenarios.

Conclusion
In this article, the HFMAE is proposed to solve the large-scale multiplicity FJSP with sequencedependent set-up time.In the HFMAE, an online fluid tracking policy is designed based on the proposed fluid model, and then the FRI is presented to obtain a high-quality initial solution.To further improve the initial solution in the HFMAE, the IMAE is proposed, based on the master-apprentice evolutionary algorithm framework.In the IMAE, two neighbourhood structures and the makespan estimate methods are presented to accelerate the local search process.
The experimental results on 10 large-scale multiplicity benchmark data sets and 10 randomly generated cases indicate that the proposed HFMAE outperforms the classical and state-of-the-art algorithms in solution quality, and has better stability and robustness.The HFMAE requires a much shorter computational time, especially for large-scale problems.Furthermore, the Wilcoxon test indicates that the HFMAE outperforms the four reference algorithms statistically.The real-world implementation results on actual plant data indicate that the HFMAE can improves production efficiency in realistic manufacturing scenarios.In the future, more constraints will be considered in the fluid model, such as machine failure, lot streaming and new job arrivals.

Disclosure statement
No potential conflict of interest was reported by the authors.

Figure 2 .
Figure 2. A feasible schedule for a problem with three jobs and five machines.

Figure 3 .
Figure 3. Illustration of three neighbourhood structures.

Figure 4 .
Figure 4. (a) Interval plot (with 95% confidence intervals) of the relative percentage deviation (RPD) of the hybrid fluid master-apprentice evolutionary algorithm (HFMAE) and reference algorithms; (b) mean central processing unit (CPU) time (s) of all algorithms.2SGA = two-stage genetic algorithm; GA = genetic algorithm; MAE = master-apprentice evolutionary algorithm; TS = tabu search.

Figure 5 .
Figure 5. Impact of set-up time and processing time on the objective (makespan): (a) impact of set-up time means; (b) impact of set-up time variances; (c) impact of processing time means; (d) impact of processing time variances.
p k , p ik .else if M now = Ø and ∀i ∈ M I , state(i) = 0 then End the FRI.else Update the class set K now (classes have at least one available machine is idling): ∀i ∈ M now , if k ∈ K i and Q k (time now ) > 0, then k ∈ K now .Select the class k with maximum p k value from the class set K now : k ← {K now } max_p k .Select the machine i with maximum p ik value from the machine set M now ∩ M k to process the class k :

Algorithm 2: The initialize procedure Input:
The initial solution S 1 .The total number J o of operations in S 1 .Initialization: Randomly generate D operation sequences based on the operation sequence of initial solution S 1 : {w 1 , w 2 , . . .w D }, Set d ← 1. while d ≤ D do Allocate the available machine for each operation of the operation sequence w d : Set p ← 1. while p ≤ J o do Select the machine i which completes the pth operation (r, n, j) of the operation sequence w d soonest From the available machine set M rnj : for ∀i ∈ M rnj do Calculate the completing time c rnji of the operation (r, n, j) in machine i. end for Select the machine with the i ← {M rnj } min_c rnji , p ← p + 1. Update w d :(r, n, j, i ) ← (r, n, j).{w 1 , w 2 , . . .w D }. S best ← {w 1 , w 2 , . . .w D } min_makespan , S 2 ← S best .
Output:The initial solution S 2 .
length of the tabu list length list .The iteration size with no improving iter number .iter number ← 0, length list ← length min .iter stop ← iter min + H × g. while iter number < iter stop do Generate the neighbourhood solutions and estimate these solutions with the proposed makespan estimate methods.if the aspiration criteria are satisfied then iter number ← 0. Update the current solution and best solution with the solution which satisfies the aspiration criteria.else iter number ← iter number + 1. Select the best no tabu solution as the current solution.end if Update the tabu list length.length list ← length min + (length max − length min ) iter number

Table 1 .
Computational results in small-scale instances.

Table 2 .
Comparison of the hybrid fluid master-apprentice evolutionary algorithm (HFMAE) and reference algorithms on generated benchmark data sets.

Table 3 .
Wilcoxon test considering the benchmark data sets.

Table 4 .
Comparison of the hybrid fluid master-apprentice evolutionary algorithm (HFMAE) and reference algorithms on the generated data sets.

Table 5 .
Comparison of the hybrid fluid master-apprentice evolutionary algorithm (HFMAE) and reference algorithms on the manufacturing data sets.FCFS = first come, first served; SPT = shortest processing time; LPT = longest processing time; MAE = master-apprentice evolutionary algorithm; 2SGA = two-stage genetic algorithm; RPD = relative percentage deviation.