An efficient evolutionary grey wolf optimizer for multi-objective flexible job shop scheduling problem with hierarchical job precedence constraints

Journal Pre-proofs An efficient evolutionary grey wolf optimizer for multi-objective flexible job shop scheduling problem with hierarchical job precedence constraints Zhenwei Zhu, Xionghui Zhou PII: DOI: Reference:

S0360-8352(20)30014-0 https://doi.org/10.1016/j.cie.2020.106280 CAIE 106280

To appear in:

Computers & Industrial Engineering

Received Date: Revised Date: Accepted Date:

16 May 2019 11 November 2019 8 January 2020

Please cite this article as: Zhu, Z., Zhou, X., An efficient evolutionary grey wolf optimizer for multi-objective flexible job shop scheduling problem with hierarchical job precedence constraints, Computers & Industrial Engineering (2020), doi: https://doi.org/10.1016/j.cie.2020.106280

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier Ltd.

An efficient evolutionary grey wolf optimizer for multi-objective flexible job shop scheduling problem with hierarchical job precedence constraints An efficient evolutionary grey wolf optimizer for multi-objective flexible job shop scheduling problem with hierarchical job precedence constraints Authors: Zhenwei Zhu ([email protected]), Xionghui Zhou (corresponding author, [email protected]) Affiliation: National Engineering Research Center of Die and Mold CAD, Institute of Forming Technology and Equipment, Shanghai Jiao Tong University (all) Address: National Engineering Research Center of Die and Mold CAD, Shanghai Jiao Tong University, 1954 Huashan RD. Xuhui District, Shanghai 200030, China Declarations of interest: none

Abstract Concentrated on the production scheduling of complex products that are assembled by multiple and multilevel manufactured parts, this paper studies the flexible job shop scheduling problem with job precedence constraints (FJSSP-JPC). Distinguished from traditional scheduling model that only considers sequential operation precedence constraints defined by process routings, FJSSP-JPC takes additional hierarchical job precedence constraints defined by Bills-of-Materials (BOMs) of final products into account. A mixed integer programming mathematical model is formulated to describe FJSSP-JPC. To represent feasible solutions that satisfy the hybrid precedence constraints, a novel three-vector encoding scheme and a job precedence repair mechanism based on binary tree are elaborated. Subsequently, this paper proposes an efficient evolutionary multi-objective grey wolf optimizer (EMOGWO) to tackle FJSSP-JPC with minimizing the objectives of makespan, maximum machine workload and total machine workload simultaneously. The algorithm involves an improved social hierarchy and a diverse leader strategy to enhance the convergence speed and population diversity separately. Statistical experiments demonstrate that EMOGWO overwhelms other competing algorithms in aspects of cardinality, convergence and diversity metrics on the majority of test instances.

Keyword Flexible Job Shop Scheduling, Job Precedence Constraint, Assembly Job Shop Scheduling, Multiple Tree-structure Constraint Optimization, Multi-objective Optimization

1 Introduction For modern manufacturing enterprises, efficient production scheduling plays a vital role in enhancing the productivity and profitability. Scheduling is a decision-making process that allocates limited resources to a collection of tasks over time with the basic aim to ensure an effective use of resources (Branke et al., 2016). Owing to the growing consumer demand for product variation, extensive researches have been devoted to the flexible job shop scheduling problem (FJSSP) encountered in the flexible manufacturing. As an extension of the classical job shop scheduling problem (JSSP), FJSSP assumes that each job consists of a sequence of consecutive operations and each operation can be processed on one machine selected from a set of alternative machines rather than a prescribed machine. In the real-world production scheduling, many industrial products are complex assemblies made up of multiple and multilevel manufactured parts that are organized as Bills-Of-Materials (BOMs). As a result, multiple scheduling jobs become connected in tree structures by precedence constraints that correspond to part dependencies (Xu & Nagi, 2013). However, traditional FJSSP only considers sequential precedence constraints between operations within each job, while jobs are independent. Aimed at tackling the complex assembly production scheduling, this paper studies a variant of FJSSP with considering the job precedence constraints (called FJSSP-JPC). As illustrated in Fig. 1, there are tree-structure job precedence constraints and line-structure operation precedence constraints in FJSSP-JPC. To reflect the multiple assessment requirements in the practical production, this paper takes into account three objectives, namely, makespan, maximum machine workload and total machine workload. The studied multi-objective flexible job shop scheduling problem with job precedence constraints is abbreviated as MOFJSSP-JPC. The additional job precedence constraints increase the complexity of addressed problem dramatically. One crucial challenge lies in satisfying the hybrid sequential and hierarchical precedence constraints simultaneously to guarantee the feasibility of solutions. Furthermore, FJSSP-JPC is a strongly non-deterministic polynomial-time hard (NP-hard) problem since FJSSP has been proved to be a NP-hard problem (Brandimarte, 1993). Manifold metaheuristics have been extensively developed to cope with NP-hard combinatorial (discrete) optimization problems with low computational effort and time. Therefore, another critical challenge is developing an efficient metaheuristic to solve the addressed problem. Inspired from the social hierarchy and hunting mechanism of grey wolves in nature, Grey Wolf Optimizer (GWO) (Mirjalili et al., 2014) is one of the recently

proposed popular swarm intelligence metaheuristics. Plenty of researches and applications demonstrate that GWO is superior at solution quality and convergence speed (Faris et al., 2018). GWO has been successfully applied to diversified engineering fields such as heat and power dispatch problem (Jayakumar et al., 2016), unmanned combat aerial vehicle path planning problem (Zhang, Zhou et al., 2016) and welding scheduling problem (Lu et al., 2016). The impressive characteristics of GWO including easy implementation, favorable convergence and few parameters motivate us to tailor a new multi-objective discrete optimization algorithm based on GWO called evolutionary multi-objective grey wolf optimizer (EMOGWO) for solving MOFJSSP-JPC. The remainder of this paper is structured in the following manner. Section 2 discusses the related works about constrained and multi-objective flexible job shop scheduling. Section 3 formally describes the studied scheduling problem and gives a mathematical model. Section 4 illustrates the three-vector encoding scheme and the active decoding scheme. Section 5 introduces the proposed EMOGWO algorithm in detail. In section 6, abundant experimental studies are provided and analyzed. Section 7 provides some concluding remarks and directions for future work.

Fig. 1 The precedence constraints for a FJSSP-JPC example (operation 1 precedes 2; job B, C and D precede A)

2 Literature review 2.1 Constrained flexible job shop scheduling In view of various manufacturing settings for different kinds of production scheduling applications, the variants of FJSSP with additional constraints have come

to the foreground. For example, Rossi (2014) took sequence-dependent setup times and transportation times into consideration. Gao and Pan (2016) investigated the resource constrained FJSSP in which extra resources like labor and fixtures are required by machines to perform operations. Salido et al. (2017) and Wu & Sun (2018) addressed an extension of FJSSP where machines can work at different speed levels. Lu et al. (2017) worked on FJSSP with controllable processing times in which the processing time can be controlled by allocating additional resources. Devassia et al. (2018) introduced a FJSSP with resource recovery constraints where resources are available in batches. Khoukhi et al. (2017) studied the FJSSP with machine unavailability constraints in consideration of preventive maintenance activities. Bożek and Werner (2018) researched the FJSSP with lot streaming and lot sizing of variable sublots. To deal with the time uncertainty caused by fluctuating manufacturing environments, Sun et al. (2019) and Lin et al. (2019) modeled the processing time in FJSSP as a triangular fuzzy number, while Xie & Chen (2018) represented the uncertain processing time as an interval grey number. Considering the frequent real-time events in the practical production scheduling, some researchers were dedicated to flexible job shop rescheduling problem or called dynamic FJSSP. Gao et al. (2018) studied the dynamic FJSSP when newly arrived priority jobs must be inserted into the existing schedule. Sreekara Reddy et al. (2018) addressed the dynamic FJSSP under the condition of random machine breakdown. Zadeh et al. (2018) focused on the dynamic FJSSP with variable processing times which may be altered due to the nature of production processes or the varying machine settings. Concerning the job precedence constraints, most studies on FJSSP assume that jobs are independent, while some researchers took assembly operations into account and focused on dependent jobs. This variant of FJSSP with precedence constraints among jobs usually appears in the assembly job shop scheduling problem (AJSSP). Generally, AJSSP can be classified into two categories, two-stage AJSSP and hybrid AJSSP. In the two-stage AJSSP, the first stage is a fabrication stage that performs ordinary operations, while the second stage is an assembly stage that specializes in assembly operations. Researchers including Tian et al. (2013), Liao et al. (2015), Du et al. (2016), Zhang and Wang (2018) devoted major efforts to solving the two-stage AJSSP. In essence, the first stage is the traditional FJSSP with only sequential operation precedence constraints while the second stage contains only hierarchical precedence constraints. As a result, the two types of line-structure and tree-structure constraints in two individual stages can be handled separately. By contrast, the hybrid AJSSP blends ordinary operations and assembly operations, causing the mixture of hierarchical and sequential precedence constraints. The increased complexity leads to few literatures on the hybrid AJSSP in some sense. Dileeplal (2012) developed a multi-objective genetic algorithm based on a permutation encoding scheme and devised precedence preserving evolutionary operators to resolve multi-objective AJSSP. However, this study allots each operation to a fixed machine and omits the machine flexibility. To settle single-objective AJSSP, Zou et al. (2018) proposed a level-based method that generates the initial solutions level by level in the precedence tree, and implemented level barrier-based crossover and mutation operators for

evolutionary optimization of single-objective AJSSP. But this level-based method is quite complicated and tedious considering that the encoded solution must be decomposed into multiple levels. Zhu et al. (2019) provided a semantic graph method based on Neo4j to model the scheduling problem and developed an ant colony optimization algorithm that directly manipulates the job shop data in graphs. However, this study focuses on the single-objective FJSSP-JPC and the non-coding graph method can only be applicable to the ant colony optimization algorithm. In summary, the above literature review demonstrates that the absence of a concise and efficient solution representation methodology hinders the development of metaheuristics for FJSSP-JPC. The solution representation methodology which embodies the encoding scheme, decoding scheme and relevant operators of generating new solutions ought to maintain the feasibility of solutions by obeying the hybrid precedence constraints. In addition, an effective metaheuristic for tackling the unsettled MOFJSSP-JPC remains to be researched.

2.2 Multi-objective optimization A multi-objective optimization problem (MOP) can be defined as follows: min f ( X )  min[ f1 ( X ), f 2 ( X ),..., f m ( X )] X  ( x1 , x2 ,..., xn )  R n

(2)

where fi(X) is the i-th sub-objective function, X is a vector of solution, Rn is the decision variable space. Given that A, B ∈ Rn, a solution A is said to dominate another solution B (denoted by A≺B) if and only if fi(X) ≤ fi(Y) for each i ∈ {1, ..., m} and fj(X) < fj(Y) for at least one index j ∈ {1, ..., m}. A solution X* is Pareto optimal solution if there is none solution X ∈ Rn that dominates X*. The corresponding objective function f(X*) is called Pareto optimal front vector. A set of Pareto optimal solutions is called optimal Pareto set (PS*), while the corresponding set of Pareto optimal front vectors is called optimal Pareto front (PF*). In general, the ideal optimal Pareto front contains a large amount of points. Consequently, the prime goal of MOP is to acquire a superior Pareto front (PF) with a limited number of points that should be as close as possible to PF* and uniformly spread as well. In real-world production scheduling, the schedule evaluation usually incorporates manifold performance indicators such as the maximum completion time called makespan, the maximum machine workload, the total machine workload, the total energy consumption (Zhang et al., 2017), the total additional resource consumption (Lu et al., 2017) and the total weighted tardiness (Huang & Yu, 2017). Hence, multi-objective flexible job shop scheduling problem (MOFJSSP) has been extensively researched in recent years. Kurdi (2017) categorized the approaches for solving MOFJSSP into three types: weighting approach, lexicographical order-based approach and Pareto-based approach. The weighting approach converts MOFJSSP to a single-objective FJSSP by assigning fixed weights to each objective function such as (Karthikeyan et al., 2015; Meng, Zhang, & Fan, 2016). The lexicographical order

based approach optimizes separately each objective in turn according to the importance rank such as (González-Rodríguez, Vela, & Puente, 2010). The first two approaches rely heavily on the users’ priori knowledge about the objective preference or weight assignment, which is always absent in practice. By contrast, the Pareto-based approach is more desirable since it can yield multiple non-dominated solutions in a single run for users to choose. Combining with Pareto-based approach, the multi-objective metaheuristics can search multiple regions of the solution space and classify the solutions simultaneously based on dominance relationship in each iteration to gain multiple non-dominated solutions finally (Pholdee et al., 2017). Varieties of metaheuristics have been developed to address all kinds of MOFJSSP. Due to the discrete characteristic, the multi-objective evolutionary algorithm (MOEA) has gained much attention. Shen & Yao (2015) and Zhang & Chiong (2016) resolved the MOFJSSP by using the multi-objective genetic algorithm integrated with problem-specific local search strategies. Zhao et al. (2017) provided an improved multi-objective evolutionary algorithm based on decomposition for solving MOFJSSP. Gong & Deng et al. (2018) presented a memetic algorithm based on NSGAII to solve the MOFJSSP with worker flexibility. Recently, there is a tendency to transform the swarm intelligent algorithm which mainly aims at continuous optimization problems into a discrete version. Lu et al. (2017) developed a multi-objective discrete virus optimization algorithm to tackle the MOFJSSP with controllable processing times. Gong, Han, and Sun (2018) presented a hybrid multi-objective discrete artificial bee colony algorithm that incorporates crossover and local search operators. Sreekara Reddy et al. (2018) proposed an evolutionary based multi-objective teacher learning-based optimization algorithm to address MOFJSSP with considering machine breakdown. Qin et al. (2019) designed a hybrid multi-objective discrete grey wolf optimization algorithm embedded with tabu search for MOFJSSP in the casting production. Luo, Zhang, and Fan (2019) proposed a multi-objective grey wolf optimizer that integrates continuous real-number encoding and discrete integer encoding to solve MOFJSSP with variable processing speeds. Dai et al. (2019) developed an enhanced genetic algorithm combining with particle swarm optimization and simulated annealing algorithm to address an energy-efficient MOFJSSP with transportation constraints. Besides, Huang and Yu (2017) proposed a hybrid ant colony optimization algorithm combined with tabu search for addressing the MOFJSSP with equal-size lot-splitting. García-León (2019) presented a general local search approach to acquire Pareto fronts for MOFJSSP based on disjunctive graph. As the famous No Free Lunch (NFL) theorem (Wolpert & Macready, 1997) has logically proved, there is no metaheuristic best suited for all single-objective or multi-objective optimization problems. In other words, a particular metaheuristic may generate very promising results on a set of problems, but the same algorithm may show poor performance on a different set of problems. This motivates researchers to modify current approaches or invent brand new metaheuristics for tackling various optimization problems persistently. Thus, it is worthwhile to devise a fresh metaheuristic for MOFJSSP-JPC.

3 Problem formulation 3.1 Problem description In brief, FJSSP-JPC can be described as follows. There are a set of gn job groups G = {G1, G2, …, Ggn} and a set of mn machines M = {M1, M2, …, Mmn}. A job group Gi takes charge of manufacturing one type of final product and is composed of jni multilevel jobs {Ji,1, Ji,2, …, J i , jni } that correspond to the tree-structure BOM. A job Ji,j takes charge of producing one type of intermediate part and is composed of a sequence of oni,j consecutive operations {Oi,j,1, Oi,j,2, …, Oi , j ,oni , j } in accordance with the process routing. An operation Oi,j,k must be processed on a machine selected from a set of alternative machines AMi,j,k. Distinct from traditional FJSSP, all jobs are connected in a hierarchical job precedence tree for each job group. Meanwhile, operations are connected in a linear operation sequence for each job. An example of FJSSP-JPC is given in Table 1 and Fig. 2. The scheduling problem can be decomposed into two sub-problems: (1) assigning each operation to an appropriate machine and (2) sequencing the assigned operations on each machine. The aim is to minimize three objectives concurrently: (1) the maximum completion time called makespan, (2) the total workload of all machines and (3) the maximum machine workload. Moreover, the following assumptions are considered for this addressed problem: (1) Each operation must be processed without interruption. (2) All machines are available at time zero and all jobs are released at time zero. (3) Transportation times between operations and setting up times of machines are negligible.

Fig. 2 An example of the job precedence constraints (J2,1 and J2,2 precede J2,3, that is, the operations of J2,3 must be processed after finishing all operations of J2,1 and J2,2,) Table 1 An example of FJSSP-JPC Index Group Job Operation Machine Processing time 1 O1,1,1 [M1,M2] [2,3] G1 J1,1 2 O1,1,2 [M3,M4] [2,3]

3 4 5 6 7 8 9 10 11 12 13 14

J1,2 J1,3 J1,4 J1,5 J2,1 G2

J2,2 J2,3

O1,2,1 O1,2,2 O1,3,1 O1,4,1 O1,4,2 O1,4,3 O1,5,1 O2,1,1 O2,1,2 O2,2,1 O2,2,2 O2,3,1

[M1,M2] [M3,M4] [M5,M6] [M1,M2] [M3,M4] [M5,M6] [M5,M6] [M1,M2] [M5,M6] [M3,M4] [M5,M6] [M5,M6]

[3,4] [3,4] [4,5] [4,5] [4,5] [2,3] [2,3] [3,4] [2,3] [3,4] [2,3] [4,5]

3.2 Mathematical modelling The studied FJSSP-JPC is formulated into a mixed integer programming (MIP) model. The notations and constant parameters used throughout this paper are firstly given below. Then the decision variables, subjected constraints and objective functions are mathematically depicted. 1) Notations Gi the ith job group, i = 1,2,…, gn Ji,j the jth job of job group Gi, j = 1,2,…, jni Oi,j,k the kth operation of job Ji,j, k = 1,2,…,oni,j Ml the lth machine, l = 1,2,…,mn 2) Constant parameters gn the total number of job groups jni the total number of jobs in job group Gi oni,j the total number of operations in job Ji,j mn the total number of machines pti,j,k,l the processing time of operation Oi,j,k on machine Ml, 1, if operation Oi , j ,k can be processed on machine M l mai , j ,k ,l   0, otherwise AM i , j ,k  {M l | mai , j ,k ,l  1, l  1, 2,.., mn} , a set of alternative machines for Oi,j,k

1, if job J k ,m precedes J k ,n jpk ,m ,n   0, otherwise PJ i , j  {J i ,q | jpi ,q , j  1, q  1, 2,..., jni } , a set of all precedent jobs over job Ji,j j 1

i 1 jnq

p 1

q 1 p 1

Index(Oi , j ,k )  k   oni , p   onq , p , the fixed ID of operation Oi,j,k

jni

oni   oni , p , the total number of operations in job group Gi p 1

gn

jnq

on   onq , p , the total number of operations q 1 p 1

3) Decision variables 1, if operation Oi , j ,k is processed on machine M l X i , j , k ,l   0, otherwise Yi , j ,k ,r , s ,t

1, if O i , j ,k is the adjacent predecessor of Or , s ,t on the same machine   1, if O i , j ,k is the adjacent successor of O r , s ,t on the same machine  0, otherwise

sti,j,k the start time of operation Oi,j,k eti,j,k the end time of operation Oi,j,k eti,j the end time of job Ji,j 4) Constraints mn

eti , j ,k  sti , j ,k   X i , j ,k ,l pti , j ,k ,l , Oi , j ,k

(3)

eti , j  eti , j ,oni , j , J i , j

(4)

sti , j ,k  eti , j ,k 1 , (Oi , j ,k 1 , Oi , j ,k )

(5)

sti , j ,1  eti ,m jpi ,m , j , ( J i , j , J i ,m )

(6)

X i , j ,k ,l  mai , j ,k ,l , (Oi , j ,k , M l )

(7)

l 1

mn

X l 1

i , j , k ,l

 1, Oi , j ,k

( sti , j ,k  etr , s ,t ) X i , j ,k ,l X r , s ,t ,lYi , j ,k ,r , s ,t (Yi , j ,k ,r , s ,t  1) ( str , s ,t  eti , j ,k ) X i , j ,k ,l X r , s ,t ,lYi , j ,k ,r , s ,t (Yi , j ,k ,r , s ,t  1)  0, (Oi , j ,k , Or , s ,t , M l ) Yi , j ,k ,r , s ,t  Yr , s ,t ,i , j ,k  0, (Oi , j ,k , Or , s ,t )

(8)

(9) (10)

5) Objectives Minimize f1  max{eti , j ,k | i  1, 2,..., gn; j  1, 2,..., jni ; k  1, 2,..., oni , j } gn

jni oni , j

Minimize f 2    (eti , j ,k  sti , j ,k ) i 1 j 1 k 1

(11) (12)

gn

jni oni , j

Minimize f3  max{  X i , j ,k ,l pti , j ,k ,l | l  0,1,..., mn}

(13)

i 1 j 1 k 1

Constraints (2)(3) calculate the end time of operations and jobs. Constraint (4) defines the operation precedence constraint and ensures that the posterior operation cannot be started until its prior operation is completed. Constraint (5) defines the job precedence constraint and ensures that the subsequent job cannot be started before its precedent jobs are totally finished. Constraint (6) describes that each operation must be processed on the given alternative machines. Constraint (7) depicts that each operation must be processed on one and only one machine. Constraint (8) guarantees that each machine processes only one operation at each time, that is, the processing time periods of two adjacent operations cannot overlap. Objective (10) is the maximum completion time called makespan. Objective (11) is the total workload of all machines. Objective (12) is the maximum machine workload.

4 Solution representation As a foremost precondition of applying the metaheuristic, an effective and efficient solution representation method must be devised to encode solutions to digits for the convenience of searching the solution space and decode digits to practical solutions for calculating objectives.

4.1 Encoding scheme The pivotal and tough point of designing an appropriate encoding scheme for FJSSP-JPC lies in satisfying the hybrid constraints simultaneously especially tree-structure job precedence constraints. This paper proposes a novel encoding scheme which is composed of three parts namely job sequencing vector (JSV), grouped operation sequencing vector (GOSV) and machine assignment vector (MAV). All these three vectors are integer arrays with a size equal to the total number of operations under consideration. In addition, a GOSV repair mechanism based on binary tree is implemented to eliminate the violations of job precedence constraints. The job sequencing vector represents the processing permutation of jobs. In JSV, the job-based representation designates the jobs of each job group as the same job group index. The kth occurrence of a job group index corresponds to the kth job of this group in GOSV. The grouped operation sequencing vector represents the processing sequence of operations and is divided into a two-dimensional array based on job groups. In GOSV, the operation-based representation designates the operations of each job as the same job index. For each group in GOSV, the kth occurrence of a job index refers to the kth operation of this job. The machine assignment vector represents the machine allocation for each operation. In MAV, the machine index elements are arranged according to the fixed IDs of operations. Supposing that JSV[i] is the kth occurrence of element j in JSV. Meanwhile,

GOSV[j][k] is the nth occurrence of element m in GOSV[j]. As a result, the ith processed operation in the whole sequence would be Oj,m,n. If the fixed ID of Oj,m,n is p and MAV[p] is equal to q, then Oj,m,n would be operated on the machine Mq. To depict this encoding scheme concretely, a solution representation example is presented in Fig. 3. In order to guarantee the satisfaction of job precedence constraints, a job precedence repair mechanism is devised. The repair mechanism adjusts the GOSV by using the binary tree sort and in-order traversal approaches. As for the original GOSV in Fig. 3, some violations of job precedence constraints occur (e.g. O1,3,1 of J1,3 is processed before O1,2,1, O1,2,2 of J1,2 and O1,1,1, O1,1,2 of J1,1, while none operation of J1,3 could be processed until all operations of J1,2 and J1,3 are finished as restricted in Fig. 2). Therefore, the GOSV repair mechanism is triggered. For each job group, a matching binary tree is constructed. The root node is the first element in the corresponding GOSV partition. Then the binary tree is constructed and branched in accordance with the operation sequence and job precedence constraints of this job group. Given that element q is posterior to element p in GOSV. If Jq precedes Jp, node q goes to the left branch of node p. Otherwise, node q goes to the right branch of node p. At last, the binary sort tree is traversed in-order (left, root, right) to acquire a repaired GOSV permutation that totally conforms to the hybrid precedence constraints.

Fig. 3 An example of encoded solution representation for the illustrative FJSSP-JPC

instance in Fig.2 and Table 1

4.2 Decoding scheme During the decoding process, a realistic schedule is restored by determining the specific start time and end time for each operation on its assigned machine in line with the encoded vectors. Once the machine assignment is determined, the total machine workload and the maximum machine workload remain constant regardless of the operation permutation. Instead of arranging each operation on its assigned machine in strict conformity with the encoded sequencing vectors, the widely used active schedule strategy (Gao et al., 2014) is implemented to compress the makespan without affecting workload objectives. Each operation can only be started after all its precedent operations are completed and ought to be scheduled as early as possible. To schedule one operation, all idle time intervals of its assigned machine are searched and the earliest time interval that can accommodate the processing time is selected to insert this operation. Detailed decoding procedures are listed in Algorithm 1. After applying this to the encoded solution in Fig. 3, the decoded schedule is displayed as Gantt chart in Fig. 4. Algorithm 1 Pseudocode of active decoding algorithm 1: Let Il be the idle time intervals set of machine Ml, Il = [[0, inf],] 2: Let jck be the number of scheduled job indexes in job group Gk, jck = 0 3: Let ocm,n be the number of scheduled operations in job Jm,n, ocm,n = 0 4: for element x in JSV do 5: y = GOSV[x][jcx+1] 6: z = ocx,y+1 7: u = Index(Ox,y,z) 8: v = MAV[u] 9: // Assign operation Ox,y,z to machine Mv minst = max(etx,y,z-1, max(et f , g ) ) // earliest start time 10: J PJ f ,g

x,y

11: for interval [ist, est] in Iv do 12: if max(minst, ist) + ptx,y,z,v <= est then 13: stx,y,z = max(minst, ist) // start time 14: Update Iv, break 15: end if 16: end for 17: etx,y,z = stx,y,z + ptx,y,z,v // end time 18: jcx += 1, ocx,y += 1 19: if (z == onx,y) then etx,y = etx,y,z 20: end for Notes: inf denotes a positive infinite number. ist and est denote the lower bound and upper bound of the idle time interval separately. stx,y,z and etx,y,z denote the start time

and end time of operation Ox,y,z. ptx,y,z,v denotes the processing time of operation Ox,y,z on machine Mv. etx,y denotes the end time of job Jx,y. onx,y denotes the number of operations in job Jx,y. PJx,y denotes the precedent jobs set of job Jx,y.

(a)

(b) Fig. 4 Gantt charts for the encoded solution in Fig. 3， (a) original decoding, (b) active decoding

4.3 Conversion of continuous encoding Different from the proposed EMOGWO and other evolutionary algorithms which always encode solutions to integer vectors, swarm intelligence algorithms usually encode solutions to real-number vectors. For the comparative study of different algorithms in the subsequent section, the continuous encoding scheme is briefly introduced. The continuous encoding representation consists of three real number vectors equivalent to JSV, GOSV and MAV. These real-number vectors of which the value range is [-1, 1] need to be converted to integer vectors for decoding. In other words, the continuous encoding scheme needs to be transformed to the discrete encoding before decoding. For the conversion of JSV, the largest position value (LPV) rule (Yuan, Xu, & Yang, 2013) is employed to construct a new job permutation by sorting the old job permutation in the descending order of their position values as

illustrated in Fig. 5. The conversion of GOSV is the same as JSV except that GOSV is partitioned into groups (see Fig. 6). For the conversion of MAV, each real number is mapped to an integer machine index by reference to the formulation MAV [i ]  1 (14) mi  AM i [round (  (| AM i | 1)  1)] 2 where i is the position index, AMi denotes the alternative machines set for the operation whose fixed ID is i, and round is a function that returns the nearest integer. An example of MAV conversion is presented in Fig. 7.

Fig. 5 Conversion from continuous JSV to discrete job permutation

Fig. 6 Conversion from continuous GOSV to discrete operation permutation

Fig. 7 Conversion from continuous MAV to discrete machine assignment

5 Proposed algorithm 5.1 The original GWO As a novel metaheuristic based on swarm intelligence, Grey Wolf Optimizer (GWO) (Mirjalili et al., 2014) mimics the social leadership and hunting mechanism of grey wolves in nature. In order to simulate the social hierarchy, the first three best solutions are considered as the alpha (α) wolf, beta (β) wolf and delta wolf (δ) respectively, while the rest are omega (ω) wolves. During the hunting (optimization) process, ω wolves update their positions under the guidance of three leading wolves, i.e., α, β and δ. Additionally, the hunting prey (optimal solution) behavior involves three main steps, searching for prey, encircling prey and attacking prey. In GWO, each wolf represents a solution that is encoded as a continuous position vector. To mathematically model the encircling behavior, the following equations are defined: D | C  X p (t )  X (t ) |

(15)

X (t  1)  X (t )  A  D

(16)

A  2a  r1  a

(17)

C  2  r2

(18)

where t denotes the current iteration, Xp is the positon vector of the prey and X is the position vector of a wolf. Vector D represents the distance between the wolf and prey. A and C are coefficient vectors. a is linearly decreased from 2 to 0 during the iteration. r1 and r2 are random vectors in range [0,1]. It is assumed that alpha, beta and delta wolves hold better knowledge about the potential position of prey. As a result, the wolves update their positions by following the leading wolves. The relevant equations are formulated as: D | C1  X   X |, D | C2  X   X |, D | C3  X   X |

(19)

X 1  X   A1  D , X 2  X   A2  D , X 3  X   A3  D

(20)

X1  X 2  X 3 (21) 3 In summary, GWO begins with a random population of grey wolves. The three leading wolves α, β and δ represent the first three best solutions gathered so far. α, β and δ estimate the position of the prey, and other wolves updates their positions randomly around the prey. The exploration and exploitation are guaranteed by the adaptive values of a and A. |A| > 1 forces wolves to diverge from the prey (exploration), while |A|<1 forces wolves to converge towards the prey (exploitation). X (t  1) 

Furthermore, parameter C also contributes to exploration by providing a more random behavior when defining the distance. The multi-objective version of GWO, MOGWO (Mirjalili, et al., 2016), integrates two new components to adapt to multi-objective optimization. The first component is a fixed-sized external archive that stores and retrieves the best non-dominated solutions so far during optimization. The archive employs an adaptive grid mechanism to improve the diversity of non-dominated solutions. The second component is a leader selection strategy that assists to choose the three leading wolves from the archive.

5.2 The proposed EMOGWO Using real-number vectors to imitate the positions of wolves, the original GWO and MOGWO mainly aim to solve the continuous optimization problems. However, the addressed FJSSP-JPC is a combinatorial optimization problem. In this paper, an evolutionary multi-objective grey wolf optimizer (EMOGWO) is adapted from GWO to tackle the multi-objective discrete optimization problems like MOFJSSP-JPC. A flowchart in Fig. 8 outlines the proposed EMOGWO. On account of the Pareto dominance nature of MOP, the optimal result is usually not a single but a set of non-dominated solutions. Thus, for each wolf, a non-domination rank is calculated by the fast non-dominated sorting approach in NSGAII (Deb et al., 2002). The whole population can be categorized into wolf packs (Pareto fronts) with different non-domination ranks in accordance with the Pareto dominance relationship. A crowding distance indicator is also involved to measure the crowdedness neighboring each solution in each ranked Pareto front. The crowding distance estimates the perimeter of the cuboid formed by using the nearest neighbors as the vertices with normalized objectives. A solution with the lower non-domination rank and the larger crowding distance is preferred. The fast non-dominated sorting approach can be concisely sketched as follows. First, for each solution p, compute (1) domination count np, the number of solutions that dominates the solution p, and (2) domination set Sp, a set of solutions that the solution p dominates. All non-dominated solutions constitute the first non-domination rank set with a zero domination count. Then, for each solution p with np = 0, visit each solution q from the domination set Sp and reduce its dominated count nq by one. Any solution q whose domination count becomes zero will be put into the second non-domination rank set. This procedure continues until all non-domination ranks are identified. Let M be the number of objectives and N be the size of population. The overall computational complexity of EMOGWO is O(MN2), which is the same as MOGWO and NSGAII. But compared with these two metaheuristics, EMOGWO is more concise with only three adjustment parameters, namely, the population size, the neighbor probe probability and the termination criterion.

Fig. 8 Flowchart of the proposed EMOGWO

5.3 Social hierarchy As core ideas of GWO, the social hierarchy and the prey search under the guidance of leading wolves are retained in EMOGWO. As mentioned earlier, the wolf population is partitioned into different ranked wolf packs based on the Pareto dominance relationship. The three supreme non-dominated wolf packs with the first three non-domination ranks are equivalent to the three leading wolf packs from which alpha, beta and delta wolves are selected separately. To acquire a uniformly spread-out Pareto optimal front, the leading wolf is selected from the corresponding leading wolf pack by using a binary tournament selection operator based on the crowding distance comparison. Instead of sharing the same alpha, beta and delta wolves at each iteration, each wolf randomly chooses its own leading wolf to follow. This diverse leader strategy makes important sense for enriching the population diversity. Furthermore, the leading wolf ought to rank no worse than the following wolf. The leader pilots the follower for exploring the solution space by applying the crossover operator. In order to avoid stagnation in local optima, the unimproved wolves can randomly probe the neighborhood with a certain probability λ by applying the mutation operator. The prey search can be formulated as below, If the wolf is at the first non-domination rank, X new  Crossover ( X a , X )

(22)

If the wolf is at the second non-domination rank, Crossover ( X a , X ), if rand  0.5 X new   Crossover ( X  , X ), otherwise

(23)

If the wolf is at the third or worse non-domination rank, X new

Crossover ( X a , X ), if rand  0.33   Crossover ( X  , X ), else if rand < 0.66  Crossover ( X  , X ), otherwise

(24)

Afterwards, if the guidance search fails to improve the solution, there is a probability λ for each wolf to probe the neighborhood.  Mutate( X ),if X  X new and rand   '  X new  X new , otherwise

(25)

where rand denotes a uniform random value in [0,1], Crossover and Mutate are the search operators introduced in the succeeding section.

5.4 Operators of searching for prey The crossover operator takes advantage of the prey knowledge owned by leader to guide the prey search of follower. With regard to the job sequencing vector, the precedence crossover similar to the improved precedence operation crossover (Zheng, Wang, & Wang, 2014) is utilized to yield new JSV. Firstly, the job group index set is randomly divided into two distinct non-empty sets S1 and S2. Thereafter, the elements of leader JSVl that belong to S1 are copied to new JSVnew with maintaining the original positions. Finally, the remaining elements are copied from the follower wolf JSVw to fill JSVnew without changing the permutation. A JSV crossover example is illustrated in Fig. 9, where S1 = {2, 3} and S2 = {1, 4}. For the grouped operation vector, new GOSV is filled group by group. The crossover operator for each OSV group is the same as JSV. A GOSV crossover example is described in Fig. 10, where S1 = {{2, 3, 5}, {1, 3}} and S2 = {{1, 4}, {2}}. For the machine assignment vector, new MAV is filled by the elements in MAVl and MAVw randomly. A MAV crossover example is presented in Fig. 11, where the random numbers are [0.4, 0.7, 0.1, 0.8, 0.6, 0.3, 0.2, 0.4, 0.8, 0.9, 0.2, 0.1, 0.7, 0.3]. The overall crossover algorithm is detailed in Algorithm 2.

Fig. 9 Crossover for job sequencing (JSV)

Fig. 10 Crossover for operation sequencing (GOSV)

Fig. 11 Crossover for machine assignment (MAV)

Algorithm 2 Pseudocode of Crossover (X1, X2) 1: X1 = [JSV1, GOSV1, MAV1], X2 = [JSV2, GOSV2, MAV2] 2: Xnew = [JSVnew, GOSVnew, MAVnew] 3: // Crossover for JSV 4: Divide the job group index set randomly into two non-empty sets S1, S2 5: Copy all elements that belong to S2 in JSV2 to an empty vector C2 steadily 6: for i = 1 to on do 7: if JSV1[i] belongs to S1 then 8: JSVnew[i] = JSV1[i] 9: else 10: JSVnew[i] = C2.pop(0) 11: end for 12: // Crossover for GOSV 13: for i = 1 to gn do 14: GOSVnew[i] = crossover_like_jsv(GOSV1[i], GOSV2[i]) 15: end for 16: // Crossover for MAV 17: for i = 1 to on do 18: if random(0,1) < 0.5 then 19: MAVnew[i] = MAV1[i] 20: else 21: MAVnew[i] = MAV2[i] 22: end for 23: return Xnew Notes: on denotes the total number of operations. gn denotes the total number of job groups. random (0,1) generates a uniform random value in [0,1]. The mutation operator contributes to averting the premature convergence by probing the neighborhood. With respect to the job sequencing vector, the swap approach is utilized to exchange a pair of different elements at different positions. For the grouped operation sequencing vector, the swap approach is also exploited to exchange a pair of different elements in the same random group. For the machine assignment vector, the neighbor is generated by just randomly changing an assigned machine to another alternative machine. The details are depicted in Algorithm 3. Algorithm 3 Pseudocode of Mutate (X) 1: X = [JSV, GOSV, MAV], Xnew = X = [JSVnew, GOSVnew, MAVnew] 2: // Swap for JSV 3: Randomly choose i1, i2 in range [1, on] with i1 ≠ i2, JSV[i1] ≠ JSV[i2] 4: JSVnew[i1] = JSV[i2], JSVnew[i2] = JSV[i1] 5: // Swap for GOSV 6: Randomly choose a number j in range [1, gn] 7: Randomly choose k1, k2 in range [1, onj] with k1 ≠ k2, GOSV[j][k1]

≠GOSV[j][k2] 8: GOSVnew[j][k1] = GOSV [j][k2], GOSVnew[j][k2] = GOSV [j][k1] 9: Repair GOSVnew[j] by using the binary tree approach 10: // Mutation for MAV 11: Randomly choose a number p in range [1, on] 12: Get the operation Ox,y,z with Index(Ox,y,z) = p 13: Randomly choose a machine Mq from AMx,y,z with q ≠ MAV[p] 14: MAVnew[p] = q 15: return Xnew Notes: gn denotes the total number of job groups. onj denotes the total number of operations of job group Gj. on denotes the total number of operations under consideration. AMx,y,z denotes the alternative machines set of operation Ox,y,z.

5.5 Elitist preservation For the sake of accelerating convergence, the elitism is introduced into EMOGWO by preserving the elite solutions. After searching for prey, the newly generated solutions are merged with the original solutions to constitute a combined population from which a new elitist population for the next iteration is carefully chosen. Since the previous and current population are both included, the elitism can be ensured. The procedures of distinguishing the elitist are detailed as follows. Firstly, the combined population is classified and sorted according to the recalculated non-domination rank in ascending order. Then each ranked non-domination set is sorted by the recalculated crowding distance in descending order. Finally, the best ranked non-domination set in the remaining combined population is picked out to fill the elitist population iteratively until no more sets can be accommodated. If the picked set exceeds the residual capacity of elitist population, the required amount of solutions are picked out from this sorted set to finish filling. In this way, the population is continually improved by the solutions with as low non-domination rank and large crowding distance as possible. The emphasis on larger crowding distance can preserve diversity among solutions of the same non-dominated rank.

6 Experimental studies This section aims to validate the effectiveness of the proposed EMOGWO by conducting experimental studies on FJSSP-JPC instances. In the following subsections, experimental instances, parameter settings and performance metrics are elaborated at first. Then the effects of improvement strategies including neighbor probe probability, improved social hierarchy and diverse leader strategy are analyzed. Subsequently, EMOGWO is compared with not only well-known multi-objective algorithms: NSAGII and MOGWO, but also the newest algorithms: HMOGWO and EGA. Finally, the reasons for the superiority of EMOGWO are discussed.

6.1 Experiment setup The competing algorithms are implemented in Python with the assistance of jmetalPy (Benitez-Hidalgo et al., 2019) package, which is an object-oriented Python-based framework for multi-objective optimization. The numerical experiments are conducted on a personal computer equipped with 2.80GHz Intel Core i5-8400 processor and 16G RAM. Fifteen FJSSP-JPC instances are constructed and categorized into three problem scales as shown in Table 2. Each experiment is conducted 20 independent times on each test instance for each algorithm. These instances are adapted from the scheduling problem encountered in the hybrid multi-category and multi-model bicycle production, which is exactly the studied FJSSP-JPC. The bicycle is typical assembly product that is assembled by bike frame, pedal and wheels. The pedal and wheel are also assembled by manufactured parts and outsourced standard parts. The workshop consists of machines and workers like metal workers, welders, paint sprayers and assemblers that conduct specific production processes respectively like cutting tubes, welding, painting and assembling. The workshop can produce different categories of bikes such as touring bike, mountain bike, and road bike as well as their relevant distinct components. Each category of product or part has several models that share the same production processes. Table 2 Detailed configuration of test instances Category Problem Machines Job groups Jobs Operations S1 5 4 24 32 S2 5 6 36 48 small S3 6 6 40 56 scale S4 6 8 48 64 S5 7 8 50 70 M6 8 10 63 81 M7 8 10 68 92 medium M8 9 12 80 104 scale M9 9 12 88 120 M10 10 12 90 125 L11 11 15 100 136 L12 11 16 108 155 large L13 12 18 125 172 scale L14 12 20 136 185 L15 12 20 150 200 It is well known that parameter settings have a significant impact on the algorithm performance. Hence, pilot experiments have been conducted to determine the ultimate parameter configuration incorporating literature reference. For the proposed EMOGWO algorithm, there are two control parameters including population size and neighbor probe probability. The pilot experiment on the neighbor probe probability is provided in Section 6.3. As a common parameter for

metaheuristics, the appropriate population size is directly given due to its non-distinctive property and the limited space. To make a fair comparison, all algorithms use the same maximum number of function evaluations (NFEs) as the termination criterion. The population size is set to 100, 150 and 200 for small, medium and large scale problems respectively, while the maximum NFEs is set to 50000, 120000 and 200000 separately. The neighbor probe probability of EMOGWO is set to 0.1. The remaining parameter configuration is listed below. For NSGAII, the following parameters are chosen:  Offspring size = 100\150\200  Crossover probability = 0.9  Mutation probability = 0.1 For MOGWO, the following parameters are chosen:  Archive size = 100\150\200  Grid inflation parameter = 0.1  Number of grids per each dimension = 10 For HMOGWO, the following parameters are chosen:  Crossover probability = 0.9  Mutation probability = 0.1  Exploration control parameter = 0.03 For EGA, the following parameters are chosen:  Crossover probability = 0.9  Self-adaptive leaning factors = (0.3, 0.3)  Iteration number for the current temperature = 10  Temperature threshold = 50 In order to mitigate the impact of randomness on metaheuristic performance, a Wilcoxon signed rank test with the significance level of 0.05 is used for assessing the significant difference between the results acquired by different metaheuristics (Derrac et al., 2011). The sign “+” or “-” indicates that the proposed EMOGWO performs significantly better or worse than its compared algorithm. The sign “=” denotes that there is no significant difference between EMOGWO and its opponent. The optimal results are highlighted with bold in the following tables.

6.2 Performance metrics To quantify the quality of solutions acquired by different multi-objective optimization algorithms, researchers have been primarily defined three aspects of performance metrics (Riquelme et al., 2015). Cardinality metrics measure the number of obtained non-dominated solutions. Accuracy metrics estimate the convergence of obtained approximation Pareto front, in other words, the closeness to the theoretical Pareto optimal front. Diversity metrics assess the distribution and spread of a Pareto front. Accordingly, the below three commonly used metrics are computed to quantitatively compare EMOGWO and its rivals. (1) Two Set Coverage (C). This metric is a binary cardinality indicator that can be

defined as:

C ( A, B) 

| x  B | y  A : y x | |B|

(26)

where A and B are two approximation Pareto fronts. C(A,B) calculates the fraction of solutions in B that are dominated by at least one solution in A. If C(A,B) is large and C(B,A) is small, then A is better than B in a sense. (2) Generational Distance (GD). This metric is an accuracy indicator that measures how far the obtained approximation Pareto front (PF) is from the Pareto optimal front (PF*). This unary measurement calculates the average Euclidean distance between each point in PF and the nearest neighbor in PF*. It is formulated as:

GD 

 min d ( x, y)

xPF

2

yPF *

| PF |

(27)

where |PF| is the size of PF. d(x,y) denotes the Euclidean distance between two points. A lower GD reflects a better convergence performance. (3) Inverted Generational Distance (IGD). This metric is an inverse version of GD but measures both the diversity and convergence of a PF. The PF with a lower IGD value is desirable and has an excellent comprehensive performance. The metric is defined as: IGD 



xPF *

min d ( x, y ) yPF

| PF * |

(28)

Note that the theoretical Pareto optimal front (PF*) is unknown. Hence, all non-dominated solutions acquired by all competing algorithms in all independent runs constitute the reference Pareto optimal front (PF*) for each instance. Besides, data normalization is carried out for eliminating the negative impact of different objective scales as follows,

fi ' 

fi  fi min , i  1, 2,3 fi max  fi min

(29)

where fimin and fimax are the minimum and maximum of the ith objective values among all gained solutions respectively.

6.3 Effectiveness analysis of improved strategies In this section, the sensitivity study on the neighbor probe probability is carried out firstly. Then, EMOGWO is compared with its two variants, EMOGWO1 in which only alpha wolves are regarded as leaders, and EMOGWO2 in which the same leaders are shared by all wolves in each iteration. If the neighbor probe probability λ is too large, the algorithm executes the mutation operator more frequently and fails to explore the solution space sufficiently

by the crossover operator. If λ is too small, the algorithm is prone to fall into premature convergence without enough mutation to escape from the local optima. For determining an appropriate value for the neighbor probe probability, a series of experiments are conducted in which the parameter λ increases from 0 to 1 with a step size of 0.1. Due to the limited space, three problems from different scales are selected to display the statistical comparison results of GD and IGD metrics as box plots in Fig. 12. In each box, the five horizontal lines represent the minimum, first quartile, median, third quartile and maximum respectively from bottom to up, and the cross points denote the outliers. These box plots reveal that GD values show a slowly rising trend while IGD values keep a relatively stable tendency with the increase of the neighbor probe probability λ. That is to say, a smaller parameter λ holds a positive impact on the convergence performance. Thus, the parameter λ = 0.1 is recommended in this paper. In order to investigate the effectiveness of the improved social hierarchy strategy, EMOGWO is compared to its variant EMOGWO1 in which only the alpha wolves from the first non-dominated rank are selected to guide the wolf population. The statistical results of all metrics are recorded in Table 3. It can be observed that EMOGWO outperforms EMOGWO1 in terms of the set coverage metric on 10 out of 15 instances. Meanwhile, EMOGWO is also superior on 11 out of 15 instances with regard to the GD metric. As for the IGD metric, EMOGWO wins on all 15 instances. In general, EMOGWO can still obtain better solutions with higher quality especially greater diversity, although there is no dramatic difference between EMOGWO and EMOGWO1. Consequently, it can be concluded that the improved social hierarchy promotes a good balance between exploration and exploitation in some extent since the leaders are chosen from different non-domination ranks. The efficiency of the diverse leader strategy is also verified by comparing EMOGWO with is variant EMOGWO2 in which all wolves share the same alpha, beta and delta leaders in each iteration. As displayed in Table 4, the comparison results demonstrate obviously that EMOGWO overwhelms EMOGWO2 on all the test instances with statistical confidence in terms of set coverage, GD and IGD metrics. It signifies that the diverse leader strategy can boost a significant enlargement on both the quality and diversity of solutions. Moreover, GD and IGD values of EMOGWO are more stable, which reflects that the diverse leader strategy can also reinforce the algorithm stability. More specifically, the individual leader selection can enlarge the search space of each iteration to achieve a faster convergence speed and more uniform distribution.

(a)

(b)

(c)

(d)

(e)

(f) Fig. 12 Sensitivity study of the neighbor probe probability parameter, (a) GD metric for problem S1, (b) IGD metric for problem S1, (c) GD metric for problem M6, (d) IGD metric for problem M6, (e) GD metric for problem L11, (f) IGD metric for problem L11

Table 3 Comparison results of all metrics between EMOGWO (A) and EMOGWO1 (B) Problem

Set Coverage C(A,B) C(B,A)

GD (mean/std) EMOGWO

EMOGWO1

IGD (mean/std) win

EMOGWO

EMOGWO1

win

S1

0.168

0.280

3.71E-03/ 9.96E-04 4.55E-03/ 2.05E-03

=

2.44E-02/ 3.49E-03 2.59E-02/ 3.35E-03

=

S2

0.303

0.269

3.60E-03/ 1.18E-03 3.64E-03/ 1.11E-03

=

2.35E-02/ 6.20E-03 2.55E-02/ 8.69E-03

=

S3

0.295

0.252

3.90E-03/ 1.13E-03 4.66E-03/ 1.76E-03

=

2.63E-02/ 6.12E-03 3.14E-02/ 8.39E-03

=

S4

0.335

0.455

4.13E-03/ 1.25E-03 4.05E-03/ 1.09E-03

=

3.39E-02/ 4.80E-03 3.39E-02/ 5.28E-03

=

S5

0.195

0.568

5.46E-03/ 1.09E-03 5.37E-03/ 2.44E-03

=

3.45E-02/ 6.17E-03 3.67E-02/ 1.65E-02

=

M6

0.507

0.369

2.67E-03/ 6.05E-04 2.95E-03/ 6.08E-04

=

2.86E-02/ 4.34E-03 3.29E-02/ 5.80E-03

+

M7

0.508

0.328

2.24E-03/ 5.09E-04 2.33E-03/ 3.33E-04

=

2.73E-02/ 4.29E-03 2.92E-02/ 3.09E-03

=

M8

0.667

0.244

4.12E-03/ 1.45E-03 4.39E-03/ 1.18E-03

=

4.80E-02/ 1.84E-02 5.29E-02/ 1.57E-02

=

M9

0.432

0.367

2.16E-03/ 4.89E-04 2.18E-03/ 4.53E-04

=

3.29E-02/ 5.12E-03 3.36E-02/ 5.28E-03

=

M10

0.392

0.356

2.16E-03/ 4.53E-04 2.21E-03/ 5.61E-04

=

3.12E-02/ 4.84E-03 3.28E-02/ 6.11E-03

=

L11

0.508

0.301

1.58E-03/ 4.26E-04 1.72E-03/ 4.29E-04

=

2.40E-02/ 3.04E-03 2.67E-02/ 4.46E-03

+

L12

0.292

0.448

1.52E-03/ 3.61E-04 1.49E-03/ 3.44E-04

=

2.84E-02/ 5.75E-03 3.16E-02/ 6.10E-03

=

L13

0.401

0.415

1.61E-03/ 3.97E-04 1.69E-03/ 2.94E-04

=

2.60E-02/ 3.65E-03 3.02E-02/ 5.68E-03

+

L14

0.563

0.202

1.69E-03/ 4.70E-04 1.91E-03/ 4.04E-04

=

2.71E-02/ 4.16E-03 3.33E-02/ 6.67E-03

+

L15

0.422

0.340

1.65E-03/ 3.33E-04 1.56E-03/ 2.67E-04

=

2.50E-02/ 4.97E-03 2.72E-02/ 5.61E-03

=

Table 4 Comparison results of all metrics between EMOGWO (A) and EMOGWO2 (B) Proble m

Set Coverage C(A,B) C(B,A)

GD (mean/std) EMOGWO

EMOGWO2

IGD (mean/std) win

EMOGWO

EMOGWO2

win

S1

0.543

0.097

4.17E-03/ 1.17E-03 6.26E-03/ 1.62E-03

+

2.57E-02/ 3.35E-03 4.51E-02/ 1.46E-02

+

S2

0.492

0.126

3.02E-03/ 1.00E-03 4.84E-03/ 1.58E-03

+

2.30E-02/ 5.50E-03 3.75E-02/ 1.18E-02

+

S3

0.770

0.023

3.17E-03/ 9.12E-04 5.96E-03/ 1.57E-03

+

2.42E-02/ 5.34E-03 4.65E-02/ 1.27E-02

+

S4

0.628

0.166

3.88E-03/ 1.22E-03 5.71E-03/ 1.76E-03

+

3.20E-02/ 3.86E-03 4.35E-02/ 1.01E-02

+

S5

0.796

0.036

4.19E-03/ 8.40E-04 9.61E-03/ 3.06E-03

+

3.02E-02/ 5.01E-03 7.02E-02/ 2.47E-02

+

M6

0.915

0.032

2.10E-03/ 4.35E-04 3.50E-03/ 8.47E-04

+

2.73E-02/ 3.55E-03 4.49E-02/ 7.90E-03

+

M7

0.772

0.077

1.89E-03/ 4.27E-04 3.13E-03/ 9.36E-04

+

2.54E-02/ 4.00E-03 4.03E-02/ 1.11E-02

+

M8

0.985

0.004

3.67E-03/ 1.22E-03 7.35E-03/ 1.66E-03

+

4.33E-02/ 1.59E-02 1.17E-01/ 3.07E-02

+

M9

0.962

0.009

2.01E-03/ 4.78E-04 4.61E-03/ 1.15E-03

+

2.91E-02/ 4.07E-03 7.11E-02/ 2.09E-02

+

M10

0.957

0.008

2.06E-03/ 4.53E-04 3.93E-03/ 8.20E-04

+

2.97E-02/ 4.51E-03 6.01E-02/ 1.22E-02

+

L11

0.941

0.017

1.42E-03/ 3.70E-04 3.05E-03/ 7.20E-04

+

2.35E-02/ 2.67E-03 5.01E-02/ 1.16E-02

+

L12

0.962

0.006

1.41E-03/ 3.57E-04 2.99E-03/ 6.68E-04

+

2.70E-02/ 5.96E-03 6.51E-02/ 1.84E-02

+

L13

1.000

0.000

1.47E-03/ 4.10E-04 3.95E-03/ 9.28E-04

+

2.45E-02/ 3.54E-03 8.17E-02/ 1.70E-02

+

L14

0.996

0.000

1.58E-03/ 3.86E-04 4.55E-03/ 6.89E-04

+

2.68E-02/ 4.54E-03 9.65E-02/ 1.62E-02

+

L15

0.995

0.000

1.60E-03/ 3.27E-04 4.30E-03/ 7.81E-04

+

2.45E-02/ 5.17E-03 8.06E-02/ 1.64E-02

+

6.4 Performance comparison with other MOP algorithms To further validate its effectiveness, EMOGWO is compared with four multi-objective optimization metaheuristics including NSGAII (Deb et al., 2002), MOGWO (Mirjalili et al., 2016), HMOGWO (Luo et al., 2019) and EGA (Dai et al., 2019). NSGAII is one of the most popular and classical multi-objective evolutionary algorithms. MOGWO is a multi-objective version of GWO introduced in previous section 5.1. HMOGWO is a newly reported metaheuristic adapted from GWO that integrates continuous real-number encoding and discrete integer encoding for solving MOFJSSP with variable processing speeds. EGA is also a recently reported metaheuristic that incorporates genetic algorithm, particle swarm optimization and simulated annealing algorithm for solving MOFJSSP with transportation constraints. These competing algorithms are re-implemented and customized to tackle the studied FJSSP-JPC. NSGAII uses the same encoding scheme, decoding scheme, crossover operator and mutation operator as EMOGWO. MOGWO adopts the continuous encoding scheme in section 4.3. HMOGWO uses a hybrid solution representation method that encodes JSV and OSV to real-number vectors and MAV to integer vectors. EGA utilizes the same discrete encoding scheme as EMOGWO and applies the local search to the non-dominated solution with the largest crowding distance in each iteration. Specified parameter settings for every algorithm refer to section 6.1. All the statistical results on set coverage, GD and IGD metrics are separately presented in Table 5-7. With regard to the set coverage metric in Table 5, it is obvious that EMOGWO excels at acquiring a larger quantity of non-dominated solutions. Especially compared with MOGWO, HMOGWO and EGA, EMOGWO retains an absolute predominance of the cardinality measurement for the vast majority of test instances where C(A,B) is approximate to 1 and C(B,A) is approximate to 0. C(A,B) = 1 implies that all solutions of the opponent (B) are dominated by some solutions in EMOGWO (A), while C(B,A) = 0 means that none solution in EMOGWO can be dominated by someone in the opponent. Concerning the GD metric in Table 6, EMOGWO overwhelms other rivals with significant statistical confidence on almost all test instances. It indicates that EMOGWO has an excellent convergence performance which benefits from the valid social hierarchy and diverse guidance search. As for the IGD metric in Table 7, EMOGWO outperforms other opponents significantly on the small scale instances S1-S5, the medium scale instances M6-M10, and a large scale instance L11. Additionally, EMOGWO is still competitive and only outperformed by NSGAII on the large scale instances L12-L15. Incorporating the set coverage and GD metrics on the instances L12-L15, it can be summarized that EMOGWO exhibits overwhelming advantages in the cardinality and convergence along with the competitive ability in maintaining the diversity for the large scale problems. The reason is that the social hierarchy facilitates a faster convergence by focusing on exploring the more promising region, but sacrifices the diversity slightly in the

meantime. Compared with NSGAII, the superiority of EMOGWO lies in the leader guided search mode based on social hierarchy which proves to be efficient in enhancing the solution quality and convergence performance. The discrete encoding scheme used by EMOGWO may be more appropriate for the studied combinatorial optimization problem than the continuous encoding scheme used by MOGWO, since many neighboring continuous real-number vectors are mapped to the same discrete integer number vector which cause the invalid search. HMOGWO breaks up the linkage between operation sequencing and machine assignment by optimizing them separately which causes the waste of computational effort, while EMOGWO optimizes them concurrently. Compared with EGA, EMOGWO takes advantage of diverse leader solutions from better non-domination ranks instead of identical best solution in case of stagnation into local optima. On the whole, EMOGWO maintains a better balance between exploitation and exploration. To visualize the performance of competing algorithms, the Pareto fronts obtained by each algorithm from one independent run for three different scale problems are plotted in the three-dimension coordinate space as shown in Fig. 13. The 3D scatter points are further projected onto the three two-dimension coordinate planes. It reveals that the three objectives under consideration are in conflict. For example, the improvement of makespan may cause the deterioration of maximum and total machine workload. Besides, these scatter plots are consistent with some conclusions derived from the above computational statistical results. The Pareto front attained by EMOGWO occupies the forefront and is distributed widely as well as uniformly. Ultimately, it can be concluded that EMOGWO apparently outperforms the competing algorithms.

Table 5 Comparison results of set coverage (C) metric between EMOGWO, NSGAII, MOGWO, HMOGWO, and EGA. EMOGWO (A) vs. NSGAII EMOGWO (A) vs. MOGWO EMOGWO (A) vs. HMOGWO EMOGWO (A) vs. EGA Proble (B) (B) (B) (B) m C(A,B) C(B,A) C(A,B) C(B,A) C(A,B) C(B,A) C(A,B) C(B,A) S1 0.368 0.124 0.711 0.032 0.879 0.005 0.500 0.118 S2 0.476 0.172 0.925 0.019 0.979 0.010 0.845 0.016 S3 0.638 0.069 0.928 0.000 0.976 0.008 0.793 0.023 S4 0.409 0.264 0.922 0.006 0.952 0.000 0.974 0.000 S5 0.488 0.209 0.868 0.007 0.994 0.000 0.864 0.000 M6 0.657 0.198 0.988 0.000 0.991 0.003 0.995 0.000 M7 0.662 0.158 0.990 0.000 0.959 0.002 0.954 0.002 M8 0.819 0.113 1.000 0.000 0.986 0.000 0.970 0.000 M9 0.708 0.187 1.000 0.000 0.986 0.000 1.000 0.000 M10 0.611 0.187 1.000 0.000 0.980 0.002 1.000 0.000 L11 0.597 0.200 1.000 0.000 0.991 0.002 1.000 0.000 L12 0.556 0.241 1.000 0.000 0.964 0.010 0.952 0.000 L13 0.408 0.254 1.000 0.000 0.996 0.000 0.994 0.000 L14 0.507 0.217 1.000 0.000 0.980 0.002 0.995 0.000 L15 0.588 0.190 1.000 0.000 0.998 0.000 1.000 0.000

Table 6 Mean, standard deviation (std) value and Wilcoxon test (win) of GD metric obtained by EMOGWO, NSGAII, MOGWO, HMOGWO, and EGA EMOGWO NSGAII MOGWO HMOGWO EGA Problem mean/std mean/std win mean/std win mean/std win mean/std win S1 3.00E-03/ 8.43E-04 4.05E-03/ 1.64E-03 + 6.69E-03/ 1.28E-03 + 1.04E-02/ 2.09E-03 + 7.41E-03/ 2.69E-03 + S2 2.63E-03/ 8.65E-04 3.35E-03/ 7.73E-04 + 7.96E-03/ 1.59E-03 + 8.17E-03/ 1.05E-03 + 7.07E-03/ 2.25E-03 + S3 2.88E-03/ 7.91E-04 4.50E-03/ 8.18E-04 + 7.05E-03/ 1.74E-03 + 7.98E-03/ 1.92E-03 + 6.08E-03/ 2.12E-03 + S4 3.33E-03/ 1.05E-03 4.29E-03/ 1.36E-03 = 7.93E-03/ 1.93E-03 + 8.90E-03/ 1.49E-03 + 1.09E-02/ 1.87E-03 + S5 3.56E-03/ 7.29E-04 5.27E-03/ 1.07E-03 + 9.31E-03/ 1.67E-03 + 1.46E-02/ 1.62E-03 + 1.29E-02/ 1.86E-03 + M6 1.96E-03/ 4.62E-04 2.53E-03/ 4.58E-04 + 5.50E-03/ 9.61E-04 + 4.94E-03/ 1.01E-03 + 7.60E-03/ 1.24E-03 + M7 1.72E-03/ 4.12E-04 2.39E-03/ 3.73E-04 + 4.67E-03/ 8.54E-04 + 3.91E-03/ 5.32E-04 + 6.69E-03/ 1.31E-03 + M8 3.19E-03/ 1.15E-03 4.71E-03/ 1.22E-03 + 9.24E-03/ 1.74E-03 + 1.42E-02/ 1.66E-03 + 1.13E-02/ 2.38E-03 + M9 1.82E-03/ 4.21E-04 2.79E-03/ 5.89E-04 + 6.17E-03/ 1.12E-03 + 5.92E-03/ 1.53E-03 + 8.31E-03/ 1.12E-03 + M10 1.75E-03/ 3.58E-04 2.24E-03/ 3.25E-04 + 6.95E-03/ 1.22E-03 + 4.54E-03/ 1.03E-03 + 8.17E-03/ 1.52E-03 + L11 1.18E-03/ 2.98E-04 1.71E-03/ 2.89E-04 + 4.64E-03/ 8.61E-04 + 3.62E-03/ 9.14E-04 + 6.92E-03/ 9.68E-04 + L12 1.21E-03/ 3.31E-04 1.66E-03/ 3.86E-04 + 4.61E-03/ 6.18E-04 + 2.90E-03/ 6.09E-04 + 6.32E-03/ 1.20E-03 + L13 1.19E-03/ 3.13E-04 1.43E-03/ 2.98E-04 = 6.48E-03/ 9.19E-04 + 3.22E-03/ 7.54E-04 + 7.48E-03/ 7.12E-04 + L14 1.45E-03/ 4.15E-04 1.74E-03/ 2.85E-04 + 8.23E-03/ 1.04E-03 + 4.06E-03/ 3.99E-04 + 9.51E-03/ 1.01E-03 + L15 1.38E-03/ 2.87E-04 1.81E-03/ 4.18E-04 + 7.91E-03/ 1.30E-03 + 3.99E-03/ 7.94E-04 + 1.05E-02/ 1.03E-03 +

Table 7 Mean, standard deviation (std) value and Wilcoxon test (win) of IGD metric obtained by EMOGWO, NSGAII, MOGWO, HMOGWO, and EGA EMOGWO NSGAII MOGWO HMOGWO EGA Problem mean/std mean/std win mean/std win mean/std win mean/std win S1 2.21E-02/ 2.86E-03 2.47E-02/ 4.44E-03 + 4.74E-02/ 1.02E-02 + 8.07E-02/ 1.67E-02 + 3.55E-02/ 9.86E-03 + S2 2.30E-02/ 4.69E-03 2.65E-02/ 5.14E-03 + 5.23E-02/ 1.01E-02 + 6.99E-02/ 7.63E-03 + 4.67E-02/ 1.70E-02 + S3 2.46E-02/ 4.42E-03 2.72E-02/ 8.72E-03 = 5.83E-02/ 1.67E-02 + 5.60E-02/ 1.27E-02 + 4.27E-02/ 1.58E-02 + S4 3.07E-02/ 3.64E-03 3.70E-02/ 8.43E-03 + 5.89E-02/ 1.02E-02 + 6.56E-02/ 9.19E-03 + 8.56E-02/ 1.38E-02 + S5 3.15E-02/ 4.04E-03 3.77E-02/ 6.99E-03 + 7.43E-02/ 2.06E-02 + 1.24E-01/ 2.37E-02 + 9.77E-02/ 2.29E-02 + M6 2.42E-02/ 3.24E-03 2.82E-02/ 3.19E-03 + 6.69E-02/ 1.14E-02 + 5.52E-02/ 9.03E-03 + 8.97E-02/ 1.26E-02 + M7 2.55E-02/ 3.24E-03 2.85E-02/ 2.78E-03 + 5.79E-02/ 1.09E-02 + 4.86E-02/ 8.78E-03 + 8.10E-02/ 1.95E-02 + M8 4.04E-02/ 1.18E-02 4.88E-02/ 1.01E-02 + 1.24E-01/ 1.86E-02 + 1.90E-01/ 2.49E-02 + 1.40E-01/ 3.89E-02 + M9 2.73E-02/ 3.93E-03 3.10E-02/ 5.08E-03 + 9.95E-02/ 1.16E-02 + 9.97E-02/ 1.79E-02 + 1.14E-01/ 1.80E-02 + M10 2.73E-02/ 3.95E-03 3.04E-02/ 3.17E-03 + 9.69E-02/ 1.26E-02 + 6.62E-02/ 1.13E-02 + 1.07E-01/ 2.10E-02 + L11 2.06E-02/ 2.29E-03 2.34E-02/ 2.86E-03 + 8.05E-02/ 1.21E-02 + 5.87E-02/ 8.85E-03 + 9.77E-02/ 1.40E-02 + L12 2.70E-02/ 6.63E-03 2.37E-02/ 3.26E-03 = 1.04E-01/ 2.06E-02 + 6.15E-02/ 1.30E-02 + 1.15E-01/ 2.41E-02 + L13 5.81E-02/ 1.37E-02 2.63E-02/ 9.73E-03 - 2.04E-01/ 3.02E-02 + 8.63E-02/ 2.14E-02 + 1.68E-01/ 3.70E-02 + L14 3.20E-02/ 6.45E-03 2.70E-02/ 3.48E-03 - 1.58E-01/ 1.66E-02 + 6.50E-02/ 1.23E-02 + 1.69E-01/ 3.11E-02 + L15 3.36E-02/ 8.59E-03 2.70E-02/ 5.65E-03 - 1.66E-01/ 2.42E-02 + 6.89E-02/ 8.59E-03 + 1.77E-01/ 2.20E-02 +

(a)

(b)

(c) Fig. 13 Pareto fronts obtained by different multi-objective algorithms on three different scale problems, (a) S1, (b) M6, (c) L11

6.5 Discussions The primary reasons for the superiority of EMOGWO can be summarized as follows. Firstly, the effective social hierarchy contributes to the overwhelming advantage in the cardinality and accuracy of non-dominated solutions. The heuristic search process which is guided by leading solutions utilizes the embedded knowledge of leaders and focuses on exploring the more promising region in the solution space to boost the search efficiency. Secondly, the diverse leader strategy that selects a random individual leader from better non-domination ranks to follow for each solution significantly enriches the search diversity. As a result, the algorithm achieves an efficient tradeoff between convergence speed and search diversity. Thirdly, the elitist strategy that preserves the elite solutions with better non-domination rank and larger crowding distance further strengthens the convergence and still maintains the diversity simultaneously. This paper presents two different encoding schemes namely discrete encoding based on integer vectors and continuous encoding based on real-number vectors. As a combinatorial optimization problem, FJSSP consists of operation sequencing sub-problem and machine assignment sub-problem, which can be represented by discrete operation identifier permutation vector and discrete machine identifier selection vector respectively. Therefore, the continuous encoding scheme needs to be

transformed to discrete encoding scheme for decoding by mapping real-number vectors to integer number vectors. The many-to-one mapping relationship between continuous encoding and discrete encoding leads to a decline in search efficiency. For example, a continuous encoded solution with JSV = [-0.66, 0.59, 0.81, 0.73, 0.16, 0.11, 0.43, -0.02, -0.33, 0.43, -0.11, -0.28, 0.15, 0.53], GOSV = [[-0.39, -0.22, 0.59, 0.12, -0.88, 0.25, -0.66, -0.31, 0.48], [0.24, 0.48, -0.24, 0.34, -0.56]], and MAV = [0.12, -0.06, -0.48, 0.36, 0.08, 0.87, -0.55, -0.44, 0.92, -0.21, 0.22, -0.16, -0.67, 0.39] is equivalent to the illustrative continuous encoded solution in Fig. 5-7 on account of the mapping to the same encoded integer vectors. The proposed EMOGWO algorithm employs the discrete encoding scheme, while the multi-objective version of canonical GWO called MOGWO employs the continuous encoding scheme. Numerical experiments reflect that EMOGWO outperforms MOGWO comprehensively for the studied problem. The comparison result demonstrates that the discrete encoding is superior to the continuous encoding in some sense.

7 Conclusions and future work To accommodate the need for the production scheduling of complex assembly products, this paper studies the flexible job shop scheduling problem with job precedence constraints (FJSSP-JPC). In FJSSP-JPC, there are hierarchical precedence constraints between jobs as well as sequential precedence constraints between operations. The job precedence constraint trees correspond to the BOMs of final products which are made up by multilevel intermediate produced parts. To reflect the multiple assessment requirements in the practical production, three objectives are optimized concurrently including makespan, maximum machine workload and total machine workload. This study can also provide constructive enlightenment on other similar tree-structure constraint optimization problems. In conclusion, the main contributions of this paper can be outlined as follows: (1) To formally describe FJSSP-JPC, a mixed integer programming mathematical model is formulated. (2) To satisfy the hybrid hierarchical job precedence constraints, sequential operation precedence constraints and machine match constraints simultaneously, this paper designs a novel three-vector encoding scheme and a job precedence repair mechanism based on binary tree. The integer vector and real-number vector encoding schemes are both investigated. (3) Adapted from GWO, a new evolutionary multi-objective grey wolf optimizer (EMOGWO) is proposed to tackle the combinatorial optimization problems like MOFJSSP-JPC. An improved social hierarchy is employed to boost the convergence speed. A diverse leader strategy is utilized to enhance the population diversity. (4) The superiority of EMOGWO is validated by the comparison with other well-known and recent multi-objective optimization algorithms involving NSGAII, MOGWO, HMOGWO and EGA. The experimental results demonstrate that EMOGWO not only overwhelms other rivals significantly in aspects of cardinality and convergence metrics but also exhibits the excellence in diversity metric.

With the advent of green manufacturing, the future research may concentrate on energy-efficient FJSSP-JPC by considering and modeling energy efficiency-related constraints and objectives as well as developing energy consumption reduction strategies (Gao & Huang et al., 2019). Another future research direction is to design efficient ensembles of various strategies like local search operators that can be embedded into EMOGWO to enhance the algorithm performance (Gao & Cao et al., 2019).

References Benitez-Hidalgo, A., Nebro, A. J., Garcia-Nieto, J., Oregi, I., & Del Ser, J. (2019). jMetalPy: a Python Framework for Multi-Objective Optimization with Metaheuristics. arXiv:1903.02915. Bożek, A., & Werner, F. (2018). Flexible job shop scheduling with lot streaming and sublot size optimisation. International Journal of Production Research, 56(19), 6391-6411. Brandimarte, P. (1993). Routing and scheduling in a flexible job shop by tabu search. Annals of Operations Research, 41(3), 157-183. Branke, J., Su, N., Pickardt, C. W., & Zhang, M. (2016). Automated design of production scheduling heuristics: A review. IEEE Transactions on Evolutionary Computation, 20(1), 110-124. Dai, M., Tang, D., Giret, A., & Salido, M. A. (2019). Multi-objective optimization for energy-efficient flexible job shop scheduling problem with transportation constraints. Robotics and Computer-Integrated Manufacturing, 59, 143-57. Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. (2002). A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation, 6(2), 182-197. https://doi.org/10.1109/4235.996017 Derrac, J., García, S., Molina, D., & Herrera, F. (2011). A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms. Swarm and Evolutionary Computation, 1(1), 3-18. Devassia, J. V., Salazar-Aguilar, M. A., & Boyer, V. (2018). Flexible job-shop scheduling problem with resource recovery constraints. International Journal of Production Research, 56(9), 3326-3343. Dileeplal, J. (2012). Multi-objective assembly job shop scheduling using genetic algorithm and tabu search. Dissertion. Cochin University of Science And Technology. Du, H., Liu, D., & Zhang, M. H. (2016). A hybrid algorithm based on particle swarm optimization and artificial immune for an assembly job shop scheduling problem. Mathematical Problems in Engineering, 2016(2), 1-10. Faris, H., Aljarah, I., Al-Betar, M. A., & Mirjalili, S. (2018). Grey wolf optimizer: a review of recent variants and applications. Neural Computing and Applications, 30(2), 413-435. https://doi.org/10.1007/s00521-017-3272-5 Gao, K., Cao, Z., Zhang, L., Chen, Z., Han, Y., & Pan, Q. (2019). A review on swarm

intelligence and evolutionary algorithms for solving flexible job shop scheduling problems. IEEE/CAA Journal of Automatica Sinica, 6(4), 904-916. https://doi.org/10.1109/JAS.2019.1911540 Gao, K. Z., Suganthan, P. N., Pan, Q. K., Chua, T. J., Cai, T. X., & Chong, C. S. (2014). Pareto-based grouping discrete harmony search algorithm for multi-objective flexible job shop scheduling. Information Sciences, 289, 76-90. Gao, L., & Pan, Q. K. (2016). A shuffled multi-swarm micro-migrating birds optimizer for a multi-resource-constrained flexible job shop scheduling problem. Information Sciences, 372, 655-676. Gao, K., Yang, F., Zhou, M., Pan, Q., & Suganthan, P. N. (2018). Flexible Job-Shop Rescheduling for New Job Insertion by Using Discrete Jaya Algorithm. IEEE Transactions on Cybernetics, 1-12. https://doi.org/10.1109/TCYB.2018.2817240 Gao, K., Huang, Y., Sadollah, A., & Wang, L. (2019). A review of energy-efficient scheduling in intelligent production systems. Complex & Intelligent Systems. https://doi.org/10.1007/s40747-019-00122-6 García-León, A. A., Dauzère-Pérès, S., & Mati, Y. (2019). An efficient Pareto approach for solving the multi-objective flexible job-shop scheduling problem with regular criteria. Computers & Operations Research, 108, 187-200. Gong, D., Han, Y., & Sun, J. (2018). A novel hybrid multi-objective artificial bee colony algorithm for blocking lot-streaming flow shop scheduling problems. Knowledge-Based Systems, 148, 115-130. Gong, X., Deng, Q., Gong, G., Liu, W., & Ren, Q. (2018). A memetic algorithm for multi-objective flexible job-shop problem with worker flexibility. International Journal of Production Research, 56(7), 2506-2522. González-Rodríguez, I., Vela, C. R., & Puente, J. (2010). A genetic solution based on lexicographical goal programming for a multiobjective job shop with uncertainty. Journal of Intelligent Manufacturing, 21(1), 65-73. Huang, R. H., & Yu, T. H. (2017). An effective ant colony optimization algorithm for multi-objective job-shop scheduling with equal-size lot-splitting. Applied Soft Computing, 57, 642-656. Pholdee, N., Bureerat, S., & Yıldız, A. R. (2017). Hybrid real-code population-based incremental learning and differential evolution for many-objective optimisation of an automotive floor-frame. International Journal of Vehicle Design, 73(1-3), 20-53. Jayakumar, N., Subramanian, S., Ganesan, S., & Elanchezhian, E. B. (2016). Grey wolf optimization for combined heat and power dispatch with cogeneration systems. International Journal of Electrical Power & Energy Systems, 74, 252-264. Karthikeyan, S., Asokan, P., Nickolas, S., & Page, T. (2015). A hybrid discrete firefly algorithm for solving multi-objective flexible job shop scheduling problems. International Journal of Bio-Inspired Computation, 7(6), 386-401. El Khoukhi, F., Boukachour, J., & El Hilali Alaoui, A. (2017). The “Dual-Ants Colony”: A novel hybrid approach for the flexible job shop scheduling problem with preventive maintenance. Computers & Industrial Engineering, 106,

236-255. Kurdi, M. (2017). An improved island model memetic algorithm with a new cooperation phase for multi-objective job shop scheduling problem. Computers & Industrial Engineering, 111, 183-201. Liao, C. J., Lee, C. H., & Lee, H. C. (2015). An efficient heuristic for a two-stage assembly scheduling problem with batch setup times to minimize makespan. Computers & Industrial Engineering, 88, 317-325. Lin, J., Zhu, L., & Wang, Z. J. (2019). A hybrid multi-verse optimization for the fuzzy flexible job-shop scheduling problem. Computers & Industrial Engineering, 127, 1089-1100. Lu, C., Xiao, S., Li, X., & Gao, L. (2016). An effective multi-objective discrete grey wolf optimizer for a real-world scheduling problem in welding production. Advances in Engineering Software, 99, 161-176. Lu, C., Li, X., Gao, L., Liao, W., & Yi, J. (2017). An effective multi-objective discrete virus optimization algorithm for flexible job-shop scheduling problem with controllable processing times. Computers & Industrial Engineering, 104, 156-174. Luo, S., Zhang, L., & Fan, Y. (2019). Energy-efficient scheduling for multi-objective flexible job shops with variable processing speeds by grey wolf optimization. Journal of Cleaner Production, 234, 1365-1384. Meng, Q., Zhang, L., & Fan, Y. (2016). Approach of hybrid GA for multi-objective job-shop scheduling. International Journal of Modeling Simulation & Scientific Computing, 7(4), 1643006. Mirjalili, S., Mirjalili, S. M., & Lewis, A. (2014). Grey Wolf Optimizer. Advances in Engineering Software, 69, 46-61. Mirjalili, S., Saremi, S., Mirjalili, S. M., & Coelho, L. d. S. (2016). Multi-objective grey wolf optimizer: A novel algorithm for multi-criterion optimization. Expert Systems with Applications, 47, 106-119. Nebro, A. J., Durillo, J. J., Garcia-Nieto, J., Coello, C. A. C., Luna, F., & Alba, E. (2009). SMPSO: A new PSO-based metaheuristic for multi-objective optimization. In 2009 IEEE Symposium on Computational Intelligence in Multi-Criteria Decision-Making(MCDM), pp. 66-73. Riquelme, N., Von Lücken, C., & Baran, B. (2015). Performance metrics in multi-objective optimization. In 2015 Latin American Computing Conference (CLEI), pp. 1-11). Qin, H., Fan, P., Tang, H., Huang, P., Fang, B., & Pan, S. (2019). An effective hybrid discrete grey wolf optimizer for the casting production scheduling problem with multi-objective and multi-constraint. Computers & Industrial Engineering, 128, 458-476. Rossi, A. (2014). Flexible job shop scheduling with sequence-dependent setup and transportation times by ant colony with reinforced pheromone relationships. International Journal of Production Economics, 153(4), 253-267. Salido, M. A., Escamilla, J., Barber, F., & Giret, A. (2017). Rescheduling in job-shop problems for sustainable manufacturing systems. Journal of Cleaner Production,

162, S121-S132. Shen, X., & Yao, X. (2015). Mathematical modeling and multi-objective evolutionary algorithms applied to dynamic flexible job shop scheduling problems. Information Sciences, 298, 198-224. Sreekara Reddy, M.B.S. ,Ratnam, C., Rajyalakshmi, G., & Manupati, V. K. (2018). An effective hybrid multi objective evolutionary algorithm for solving real time event in flexible job shop scheduling problem. Measurement, 114, 78-90. Sun, L., Lin, L., Gen, M., & Li, H. (2019). A hybrid cooperative coevolution algorithm for fuzzy flexible job shop scheduling. IEEE Transactions on Fuzzy Systems, 27(5), 1008-1022. Tian, Y., Liu, D., Yuan, D., & Wang, K. (2013). A discrete PSO for two-stage assembly scheduling problem. The International Journal of Advanced Manufacturing Technology, 66(1-4), 481-499. Wolpert, D. H., & Macready, W. G. (1997). No free lunch theorems for optimization. IEEE Transactions on Evolutionary Computation, 1(1), 67-82. Wu, X., & Sun, Y. (2018). A green scheduling algorithm for flexible job shop with energy-saving measures. Journal of Cleaner Production, 172, 3249-3264. doi:https://doi.org/10.1016/j.jclepro.2017.10.342. Xie, N., & Chen, N. (2018). Flexible job shop scheduling problem with interval grey processing time. Applied Soft Computing, 70, 513-524. Xu, J., & Nagi, R. (2013). Solving assembly scheduling problems with tree-structure precedence constraints: A Lagrangian relaxation approach. IEEE Transactions on Automation Science and Engineering, 10(3), 757-771. Yuan, Y., Xu, H., & Yang, J. (2013). A hybrid harmony search algorithm for the flexible job shop scheduling problem. Applied Soft Computing, 13(7), 3259-3272. Zadeh, M. S., Katebi, Y., & Doniavi, A. (2019). A heuristic model for dynamic flexible job shop scheduling problem considering variable processing times. International Journal of Production Research, 57(10), 3020-3035. Zhang, Q., & Li, H. (2007). MOEA/D: A multiobjective evolutionary algorithm based on decomposition. IEEE Transactions on Evolutionary Computation, 11(6), 712-731. Zhang, S., Zhou, Y., Li, Z., & Pan, W. (2016). Grey wolf optimizer for unmanned combat aerial vehicle path planning. Advances in Engineering Software, 99, 121-136. Zhang, R., & Chiong, R. (2016). Solving the energy-efficient job shop scheduling problem: a multi-objective genetic algorithm with enhanced local search for minimizing the total weighted tardiness and total energy consumption. Journal of Cleaner Production, 112, 3361-3375. Zhang, Y., Wang, J., & Liu, Y. (2017). Game theory based real-time multi-objective flexible job shop scheduling considering environmental impact. Journal of Cleaner Production, 167, 665-679. Zhang, S., & Wang, S. (2018). Flexible assembly job-shop scheduling with sequence-dependent setup times and part sharing in a dynamic environment: constraint programming model, mixed-integer programming model, and

dispatching rules. IEEE Transactions on Engineering Management, 65(3), 487-504. Zhao, F., Chen, Z., Wang, J., & Zhang, C. (2017). An improved MOEA/D for multi-objective job shop scheduling problem. International Journal of Computer Integrated Manufacturing, 30(6), 616-640. Zheng, X., Wang, L., & Wang, S. (2014). A novel fruit fly optimization algorithm for the semiconductor final testing scheduling problem. Knowledge-Based Systems, 57, 95-103. Zhu, Z., Zhou, X., & Shao, K. (2019). A novel approach based on Neo4j for multi-constrained flexible job shop scheduling problem. Computers & Industrial Engineering, 130, 671-686. Zou, P., Rajora, M., & Liang, S. Y. (2018). A new algorithm based on evolutionary computation for hierarchically coupled constraint optimization: methodology and application to assembly job-shop scheduling. Journal of Scheduling, 21(5), 545-563.

Acknowledgement The authors want to thank the editor and anonymous referees.

Funding This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Author contributions Zhenwei Zhu: Conceptualization, Methodology, Software, Formal analysis, Investigation, Resources, Data Curation, Writing - Original Draft, Visualization. Xionghui Zhou: Writing - Review & Editing, Supervision, Project administration.