Modeling and Optimization for Multi-Objective Nonidentical Parallel Machining Line Scheduling with a Jumping Process Operation Constraint

: This paper investigates the nonidentical parallel production line scheduling problem derived from an axle housing machining workshop of an axle manufacturer. The characteristics of axle housing machining lines are analyzed, and a nonidentical parallel line scheduling model with a jumping process operation (NPPLS-JP), which considers mixed model production, machine eligibility constraints, and fuzzy due dates, is established so as to minimize the makespan and earliness/tardiness penalty cost. While the physical structures of the parallel lines in the NPPLS-JP model are symmetric, the production capacities and process capabilities are asymmetric for different models. Different from the general parallel line scheduling problem, NPPLS-JP allows for a job to transfer to another production line to complete the subsequent operations (i.e., jumping process operations), and the transfer is unidirectional. The signiﬁcance of the NPPLS-JP model is that it meets the demands of multivariety mixed model production and makes full use of the capacities of parallel production lines. Aiming to solve the NPPLS-JP problem, we propose a hybrid algorithm named the multi-objective grey wolf optimizer based on decomposition (MOGWO/D). This new algorithm combines the GWO with the multi-objective evolutionary algorithm based on decomposition (MOEA/D) to balance the exploration and exploitation abilities of the original MOEA/D. Furthermore, coding and decoding rules are developed according to the features of the NPPLS-JP problem. To evaluate the effectiveness of the proposed MOGWO/D algorithm, a set of instances with different job scales, job types, and production scenarios is designed, and the results are compared with those of three other famous multi-objective optimization algorithms. The experimental results show that the proposed MOGWO/D algorithm exhibits superiority in most instances. of great signiﬁcance for production practice. This research has a substantive impact on improving the production efﬁciency of the workshop, and can signiﬁcantly enhance the production management level of enterprises, so as to increase the market competitiveness.


Introduction
Flow shop scheduling problems are the most common and widely studied problems in the manufacturing industry, and the examined problems are usually simplified versions of the real flow shop scheduling problem so as to reduce the difficulty of modeling and solving these problems. These oversimplified schemes often cannot perfectly solve such scheduling problems in the actual environment. In real-world manufacturing situations, some special constraints or uncertainties must usually be considered and handled, such as sequence-dependent setup times [1][2][3], the kinds of parallel machines under study [4,5], machine eligibility constraints [6][7][8][9], resource constraints [10,11], and fuzzy stochastic demand [12,13]. Considerations of these additional constraints and uncertainties make the developed scheduling models closer to real production scenarios, but also increase their scheduling complexity. Because of product iteration requirements and the diversified needs of customers, production systems must often address multivariety production on multiple production lines. This is very common in manufacturing enterprises such as those in the automobile industry and household appliance industry, as well as for construction machinery manufacturers. When coping with these situations, the equipment configurations of multiple lines may be different for meeting multivariety production. In this paper, we study this nonidentical parallel production line scheduling problem. To improve the machine utilization and shorten the waiting times of the jobs to be processed, a jumping process operation is often used. This means that if a certain process operation of a job is finished, it can move to another production line to complete the subsequent process operations. This jumping process operation is unidirectional. To solve this kind of scheduling problem, nonidentical parallel line scheduling with a jumping process operation (NPPLS-JP) is proposed; this is also in essence a parallel production line scheduling problem. Notably, the proposed NPPLS-JP problem is an NP-hard problem because of its complexity.
This NPPLS-JP problem is derived from the axle housing machining workshop of an axle manufacturer. Axle housing is an important part of axle production; it usually adopts make-to-stock (MTS) production and make-to-order (MTO) assembly. The machining of an axle housing is shown in Figure 1. The axle manufacturer adopts nonidentical parallel production lines for axle housing machining. There are two parallel production lines (A and B) in the axle housing machining workshop, and each production line is installed linearly, as shown in Figure 1. The physical structures of the two production lines are symmetrical. Each production line contains five stages corresponding to five operations, and the parallel machines at any stage of each line are identical. The corresponding stages of production lines A and B have similar functions, but the configurations of the machines in the different lines are different in order to meet the needs of multivariety mixed production. In this production situation, the production load of each stage on any line is not easy to balance via a simple scheduling scheme.
parallel machines under study [4,5], machine eligibility constraints [6][7][8][9], resource constraints [10,11], and fuzzy stochastic demand [12,13]. Considerations of these additional constraints and uncertainties make the developed scheduling models closer to real production scenarios, but also increase their scheduling complexity. Because of product iteration requirements and the diversified needs of customers, production systems must often address multivariety production on multiple production lines. This is very common in manufacturing enterprises such as those in the automobile industry and household appliance industry, as well as for construction machinery manufacturers. When coping with these situations, the equipment configurations of multiple lines may be different for meeting multivariety production. In this paper, we study this nonidentical parallel production line scheduling problem. To improve the machine utilization and shorten the waiting times of the jobs to be processed, a jumping process operation is often used. This means that if a certain process operation of a job is finished, it can move to another production line to complete the subsequent process operations. This jumping process operation is unidirectional. To solve this kind of scheduling problem, nonidentical parallel line scheduling with a jumping process operation (NPPLS-JP) is proposed; this is also in essence a parallel production line scheduling problem. Notably, the proposed NPPLS-JP problem is an NP-hard problem because of its complexity.
This NPPLS-JP problem is derived from the axle housing machining workshop of an axle manufacturer. Axle housing is an important part of axle production; it usually adopts make-to-stock (MTS) production and make-to-order (MTO) assembly. The machining of an axle housing is shown in Figure 1. The axle manufacturer adopts nonidentical parallel production lines for axle housing machining. There are two parallel production lines (A and B) in the axle housing machining workshop, and each production line is installed linearly, as shown in Figure 1. The physical structures of the two production lines are symmetrical. Each production line contains five stages corresponding to five operations, and the parallel machines at any stage of each line are identical. The corresponding stages of production lines A and B have similar functions, but the configurations of the machines in the different lines are different in order to meet the needs of multivariety mixed production. In this production situation, the production load of each stage on any line is not easy to balance via a simple scheduling scheme. According to the different vehicle models, there are eight types of axle housing products that can be processed on line A and line B with multivariety mixed model production; all types of products are processed through five operations, as shown in Figure 1. Two types of axle housing products have machine eligibility constraints, and the operation "combined machining I" can only be performed at the corresponding machines in production line B; the others can finish processing on any line independently. For different types of axle housing, the different processing times required for the same operation on the same machine and the different mixing ratios of the various axle housing types increase the complexity of the scheduling problems. It is difficult to balance all the stages of each production line, so a jumping process operation is adopted to address this problem by allowing a job to be processed on two production lines. For example, when According to the different vehicle models, there are eight types of axle housing products that can be processed on line A and line B with multivariety mixed model production; all types of products are processed through five operations, as shown in Figure 1. Two types of axle housing products have machine eligibility constraints, and the operation "combined machining I" can only be performed at the corresponding machines in production line B; the others can finish processing on any line independently. For different types of axle housing, the different processing times required for the same operation on the same machine and the different mixing ratios of the various axle housing types increase the complexity of the scheduling problems. It is difficult to balance all the stages of each production line, so a jumping process operation is adopted to address this problem by allowing a job to be processed on two production lines. For example, when the operation "combined machining I" is finished on line A, the job is transferred to line B to complete Symmetry 2021, 13, 1521 3 of 24 the subsequent machining operations; this is called the jumping process operation, and the stage "combined machining I" is called the jumping process point. A jumping process operation is unidirectional, which means that the axle housings processed on line A are allowed to be transferred to production line B, but not the other way around. Appropriate jumping process operations can reduce waiting times, improve the utilization rate of equipment, balance the production capacity of each stage, and ensure the due date of each order.
From the above description of a production system, it can be seen that the axle housing machining line scheduling problem has the following characteristics: (1) The configurations of the multiple production lines are similar but not the same. (2) Mixed multivariety production is adopted to organize production. (3) Several jobs with special types have machine eligibility constraints. (4) The jumping process operation is unidirectional in the production process. This problem can be regarded as a variant and extension of the flow shop scheduling problem or general parallel production line scheduling problems, and it involves four key decision-making processes, namely: job sequencing decisions, parallel line decisions, parallel machine decisions, and job jumping process operation decisions. It is obvious that the NPPLS-JP problem proposed in this paper is a rather complex scheduling optimization problem.
Because of the fierce competition in the market and the diversified needs of customers, the NPPLS-JP problem is widespread in manufacturing environments and has an important impact on the manufacturing efficiency of production systems. However, there is no relevant research on this topic in the existing literature, so the study in this paper is of exploratory and practical significance. In this paper, we establish a scheduling model for the NPPLS-JP problem, which involves multivariety mixed model production, multiline scheduling, and machine eligibility constraints. The objectives are to minimize the makespan and earliness/tardiness penalty in a production cycle. In the NPPLS-JP model, to more closely approximate the actual production environment, the due date of a production order denoted by the fuzzy earliness/tardiness penalty model [14], and a hybrid algorithm combining the grey wolf optimization algorithm and the multi-objective evolutionary algorithm based on decomposition (MOEA/D) are proposed to solve the NPPLS-JP problem.
The remainder of this paper is organized as follows. Section 2 reviews the literature relevant to multivariety mixed model production, parallel line scheduling, and multiobjective optimization. Section 3 gives a general statement of the NPPLS-JP problem and establishes a production-based, order-oriented, multi-objective scheduling model. Section 4 proposes the multi-objective grey wolf optimizer based on decomposition (MOGWO/D) for NPPLS-JP and describes the procedures in detail. Section 5 tests the performance of the proposed MOGWO/D algorithm by comparing it with three other famous multiobjective algorithms based on a set of designed test instances, and the experimental results are analyzed. Section 6 summarizes the research content and discusses the direction of future research.

Literature Review
Mixed model production refers to the production of a variety of products on a single production line to increase the flexibility of the line and meet the multivariety and smallbatch production demands. For multiproduct demands, mixed model production is widely adopted by manufacturing enterprises, especially in assembly workshops [15,16]. Because of the widespread application of mixed model production, many scholars are devoted to research in this field and have made many achievements [17][18][19][20][21].
Mcmullen et al. [22] studied the mixed model scheduling problem with the consideration of a setup time. They presented a bean search heuristic method to obtain an efficient front. Leu and Hwang [23] proposed a resource-constrained mixed-production flow shop scheduling system for mixed precast production task problems and developed a search method based on a genetic algorithm (GA) to minimize the output makespan under re-source constraints and mixed production. Wang et al. [24] studied final assembly line scheduling, which considers order scheduling and mixed model sequencing simultaneously, and combined the original artificial bee colony (ABC) algorithm with some steps of the GA and Pareto optimality to solve this problem. Bahman et al. [25] constructed a mixed-integer linear programming model with a tighter linear relaxation for a realistic automotive industry assembly line, including a set of specific requirements involving moving workers and limited workspace. Alghazi and Kurz [26] proposed a mixed model line-balancing integer program for mixed model assembly lines with the aim of minimizing the number of chemical workers; a constraint programming model was established to address larger assembly line balancing problems.
Parallel line scheduling problems are very common in both mass production and multiproduct production. Parallel line production can enhance the stability and flexibility of the production system and improve production efficiency. When one production line breaks down, not all production activities are stopped. All parallel lines may have the same number of processing stages, the pieces of equipment in the same stage are similar and can complete the same production processes, and every line can substitute for all other lines to produce all or some types of the desired products [27]. However, in some situations, the configurations of the machines in different lines are nonidentical; this situation is more convenient for the production of multiple products. For example, two types of equipment in different production lines can independently complete the operation of milling surfaces, but the machining accuracies are different.
Haq et al. [28] studied the line scheduling problem with multiple parallel processing in job shops. Each job can only be processed on a particular line and is not allowed to move between parallel lines. Meyr and Mann et al. [29] introduced a new solution approach to determine the lot-sizing and scheduling problem for parallel production lines with the consideration of scarce capacity; sequence-dependent setup times; and the deterministic, dynamic demands of multiple products. Mumtaz et al. [30] investigated the multiple assembly line scheduling problem for a printed circuit board (PCB) assembly and developed a hybrid spider monkey optimization approach with an improved replacement strategy to solve it. Rajeswari et al. [27] presented parallel flow line scheduling with a dual objective to minimize the tardiness and earliness of jobs. All parallel flow lines had similar sets, and the authors developed a hybrid algorithm that used a GA and particle swarm optimization (PSO) to incorporate greedy randomized adaptive search to address the problem. Ebrahimipour et al. [31] proposed linear programming and a bagged binary knapsack to address the multiple production line scheduling problem. Mumtaz et al. [32] developed a mixed-integer programming model for the multilevel planning and scheduling problem of parallel PCB assembly lines, and a hybrid spider monkey optimization (HSMO) algorithm was proposed.
Multi-objective optimization refers to a situation where more than one conflicting objective is to be optimized simultaneously. It is often impossible to obtain an optimal solution as in a single-objective optimization problem, but a set of tradeoff solutions, namely, nondominated solutions, can be used to choose the most suitable solution according to the actual requirements. Therefore, it is more applicable to the actual situation, which requires the consideration of multiple indicators affecting decision making. As it is conducive to obtaining an ideal decision making effect, multi-objective optimization has a wider application field and more practical value. MOEAs are global optimization algorithms based on populations that simulate the evolution process of natural organisms. Since the whole solution set can be obtained in one run, MOEAs have become some of the mainstream algorithms in multi-objective optimization. According to their selection mechanisms, MOEAs can be classified into three classes [33]: Pareto-based algorithms, indicator-based algorithms, and decomposition-based algorithms. The Pareto-based method was first proposed by Goldberg et al. [34] and has been widely studied since then. Many classical multi-objective optimization algorithms are based on the Pareto relation, and many have been proposed based on the Pareto dominance relationship, such as the famous SPEA [35], SPEA2 [36], and NSGA-II [37]. The indicator-based method uses performance evaluation indicators to guide the search process and the choice of solutions [38,39]. The MOEA/D was first proposed by Zhang and Hui in 2007 [40], and it is the most representative multi-objective optimization algorithm based on decomposition. Different from the classical multi-objective optimization algorithms, it decomposes an input multi-objective optimization problem (MOP) into a series of single-objective optimization subproblems by using a set of uniformly distributed weight vectors and optimizing these subproblems simultaneously. Since the MOEA/D was proposed, it has attracted increasing attention from scholars, improvements and applications for the MOEA/D are constantly emerging, and it has become one of the best multi-objective optimization algorithms.
Li and Landa-Silva et al. [41] proposed evolutionary multi-objective simulated annealing (EMOSA), which incorporates a simulated annealing algorithm and introduces an adaptive search strategy. The experimental results showed that the algorithm obtained a good effect in terms of solving the multi-objective knapsack problem and multi-objective salesman problem. Tan et al. [42] developed a multi-objective meme algorithm based on decomposition, which integrates a simplified quadratic approximation (SQA) into the MOEA/D as a local search operator to balance its local and global search strategies. Wang et al. [43] designed a multi-objective particle swarm optimization algorithm based on decomposition (MPSO/D). This algorithm adopts relevant measures to ensure that only one solution is present in each subregion in oder to maintain solution diversity, and the fitness is calculated by the crowding distance. Ke and Zhang et al. [44] proposed the MOEA/D-ACO algorithm, which incorporates ant colony optimization (ACO) into the MOEA/D; they then tested the performance of the proposed algorithm in 12 instances and obtained good results. Alhindi et al. [45] developed a hybrid algorithm called MOEA/D-GLS, which integrated guided local search (GLS) with the MOEA/D to promote the exploitation ability of the original MOEA/D. The experimental results showed that the proposed MOEA/D-GLS was superior to the original MOEA/D. Zhang et al. [46] proposed MOEA/D-EGO to address expensive MOPs. In this method, the input problem is decomposed into several subproblems, and a prediction model is established for each subproblem based on the evaluated points in order to reduce modeling costs and improve the prediction quality. Wang et al. [47] proposed adaptive replacement strategies by adjusting the problem size dynamically for the MOEA/D. This approach can balance the diversity and convergence of the MOEA/D.
However, according to the no free lunch theorem [48], no algorithm can solve all of the optimization problems in all of the fields. Because of the continuous emergence of new optimization problems, the existing algorithms cannot solve these new optimization problems well, so new algorithms or improved algorithms are needed. The MOEA/D algorithm exhibits good diversity in solving MOPs; its characteristics include simplicity, few parameters, and better result distributions. In this paper, the MOGWO/D, which incorporates the GWO into the MOEA/D, is proposed to solve the NPPLS-JP problem.

Problem Definition and Assumption
Suppose that O orders are processed in L production lines, the job types for each order are the same and those for different orders may be different, and the operations of each type of job are predetermined and similar. All production lines have the same number of stages, and the machine configurations may be different. If parallel machines exist in some stages of the production line, they are identical parallel machines. Some types of jobs may have machine eligibility constraints, namely some job operations must be carried out on specific machines at certain stages of some production lines. For any production line l, if s is set as the jumping process operation point, it means that after processing is completed in stage s, the job can be transferred to line l to continue the processing of the subsequent operations (l = l and l, l ∈ {1, 2, · · · , L}); the jumping process operation is unidirectional. The scheduling objectives are to minimize the makespan and the earliness/tardiness penalty cost.
In addition, there are usually several complicated constraints and perturbations in the real-world production environment. To prevent the loss of generality and reduce the computational complexity of the scheduling model, some modeling assumptions are given, as follows.
(1) The type, quantity, and due date of each order are known.
(2) The jobs in each order are of the same type.
(3) All the machines in the production system are available at the beginning. Meanwhile, to solve the NPPLS-JP problem more conveniently, the concept of a virtual production line is introduced. This means that if there is a jumping process operation point in the manufacturing process of a job of a certain type, all the stages in two different production lines that can complete the processing of jobs of this type are regarded as a new production line, namely, a virtual production line.
In a real-world production environment, a breach of the order due dates is not always unacceptable. In general, a few occurrences of tardiness are allowed, to achieve the smallest due date penalty cost across the total orders. Here, the fuzzy due date is used to deal with this situation. Trapezoidal fuzzy due date and triangular fuzzy due date are two common fuzzy due dates that have been investigated in the literature regarding scheduling problems [14,[49][50][51]. Most researchers choose the type of fuzzy due date depending on the research background and problem characteristics. In our research, neither early nor late completion were the best solutions for the automobile industries. Completing the production order in advance will increase the inventory cost, while the order delay will lead to customer penalty loss. Finishing and delivering the orders in a given period is the most feasible result. Therefore, the trapezoidal fuzzy due date and earliness/tardiness penalty cost model was adopted in this scheduling problem according to the real-world production requirement.
As shown in Figure 2a-b, (a) is the trapezoidal fuzzy due date and earliness/tardiness penalty cost model, (b) is the corresponding satisfaction model of the fuzzy due date. Model (a) shows each order has a corresponding trapezoidal fuzzy due date that is denoted by a trapezoidal fuzzy number o , and the earliness/tardiness penalty cost coefficients are α o and β o , respectively. A completion time before d 1 o it means that the orders are produced prematurely, and additional inventory costs are generated; if the completion time comes after d 4 o , the production order seriously violates the due date requirement. These two cases are both unacceptable, so the satisfaction is 0, and the maximal earliness/tardiness penalty is imposed. If the completion time is in the time interval , the due date penalty costs decrease and increase linearly, respectively, and the corresponding due date satisfaction is just the opposite. Only in the time interval are the completion times reasonable, and the due date penalty cost in this case is 0. Therefore, the production orders should be arranged optimally in terms of time according to different due dates and earliness/tardiness penalty costs.

Mathematical Modeling
A mathematical model for the proposed NPPLS-JP problem is established using the above notations, and the two objective functions, which minimize the makespan and earliness/tardiness penalty cost, are formulated as Equations (1) and (2), respectively.
where the calculation formulas of the completion time and due earliness/tardiness penalty for order o are formulated as Equations (3) and (4), respectively: The constraints are as follows:

Mathematical Modeling
A mathematical model for the proposed NPPLS-JP problem is established using the above notations, and the two objective functions, which minimize the makespan and earliness/tardiness penalty cost, are formulated as Equations (1) and (2), respectively.
where the calculation formulas of the completion time and due earliness/tardiness penalty for order o are formulated as Equations (3) and (4), respectively: The constraints are as follows: Symmetry 2021, 13, 1521 Among the above constraints, constraints (5)- (7) together define the jumping process operation. Constraint (5) defines the operations before the jumping process operation point (including the operation on the jumping process operation stage), which can only be completed in one production line. Constraint (6) restricts all of the operations after the jumping process operation point for each job that can only be processed on the same production line. Constraint (7) states the jumping process operation for any job may or may not occur after completing another jumping process operation, and the three constraints together guarantee that the jumping process operation can only occur once at most. Constraint sets (8) and (9) define the relationships between several decision variables. Constraint (10) states that the processing of each job operation must satisfy the machine eligibility constraints. Constraint (11) ensures that each operation of a job can only be processed on one machine. Constraint sets (12) and (13) together restrict the processing sequence of two jobs on a machine to only one possible result. Constraint (14) guarantees that one machine can only process one job at a time, which means that the completion time of the current job is longer than the sum of the completion time of the immediate predecessor job and the processing time of the current job. Constraint (15) ensures that a job is processed by only one machine at a time, that is, the completion time of the job operation is greater than the sum of the completion time of the immediate predecessor operation and the processing time of the current operation.

Proposed MOGWO/D Algorithm
The MOEA/D provides a general framework that allows any single objective to be applied to the subproblems of a MOP [41]. Compared with the other multi-objective optimization algorithms, such as the Pareto-based optimization algorithms, MOEA/D has less computational complexity, and its results have better diversity. In this section, we present a hybrid algorithm that combines GWO with MOEA/D, and uses the mechanism of searching for prey in the GWO algorithm to enforce a balance between exploration and exploitation. According to the characteristics of the NPPLS-JP problem, problemspecific encoding and decoding rules are given, and some main procedures in the proposed MOGWO/D are also stated in detail.

Original GWO
The GWO was inspired by the leadership hierarchy and hunting behaviors of grey wolves [52]. In GWO, initial populations are used to simulate the grey wolf group, which is divided into four hierarchies; the solutions with the best, second, and third fitness values are α, β, and δ, respectively, are utilized to find the optimal solution by simulating the hunting process of grey wolves. In the process of hunting, the location of the prey is unknown. Therefore, to simulate the hunting behavior of grey wolves and the prey behavior from the perspective of mathematical modeling, suppose that α, β, and δ are closest to the potential position of the prey. Under the guidance of α, β, and δ, the position vector is updated to approximate the optimal solutions in the search space.
The main procedure of wolf hunting includes encircling prey and hunting, and the mathematical models for grey wolves approaching and encircling their prey are as follows: Symmetry 2021, 13, 1521 In Equations (16) and (17) (18) and (19). By changing the value of the vector → A, the search process can be guided. When | → A| > 1, α, β, and δ diverge from each other, which is good for global search; when | → A| < 1, α, β, and δ converge to the prey, which contributes to the local search. The parameter → C is generated randomly to help grey wolves jump out of the local optima.

MOGWO/D Algorithm Framework
The proposed MOGWO/D is a hybrid algorithm that integrates GWO into MOEA/D. Similar to the original MOEA/D, the MOGWO/D algorithm decomposes the input multiobjective problem into a series of single-objective scalar optimization subproblems by utilizing a set of uniformly distributed weight vectors and a scalar function. Here, we use the Tchebycheff method to construct each subproblem; then, subproblem i can be described as follows [40]: The purpose is to minimize each singleobjective function g te (x i |λ i , z * ), and each subproblem uses the approach of the GWO to update its position vectors. It is worth noting that finding accurate reference points is difficult and time-consuming work, so the best objective values z = (z 1 , z 2 , · · · , z m ) T method is used in the initial population as the initial reference point, and to update the reference point over the course of iterations by generations.
As the normalization method for objectives is conducive to increase the uniformness of the obtained solutions when the input objectives are disparately scaled [40], we used a simple normalization method to replace f i and obtained a normalization Tchebycheff approach, as follows.
where z * is the reference point and z nad = z nad 1 , z nad 2 , · · · , z nad m is the nadir point in the objective space. In our calculation, z * is replaced by z in Step 2.3 of Algorithm 1, and the maximum value of f i (x) in the current population is the substitute for z nad i . This calculation strategy can meet the needs of the algorithm.

Input:
A multiobjective problem; A stopping criterion; A set of uniformly spread weight vectors λ 1 , λ 2 , · · · , λ N ; N : population size (equal to the number of the weight vectors or subproblems); T : neighborhood size; T : the number of position vectors in the neighborhood to be updated of a subproblem (where T < T).

Output:
External population, EP for short.
Step 1.3) Randomly generate an initial population x 1 , x 2 , · · · , x N or use a problem-specific approach. The objective of each position vector is calculated and labeled as Step 1.4) Initialize z = (z 1 , z 2 , · · · , z m ) T , z i = min f i x 1 , f i x 2 , · · · , f i x N , i = 1, 2, · · · m.
Step 1.5) Calculate g te x j λ i , z * for each j ∈ B(i), and the three best position vectors are labeled as x α , x β and x δ , respectively corresponding to weight vector λ i .
Step 2) Update: for i = 1, 2, · · · ,N, do Step 2.1) Randomly select T indexes k 1 , k 2 , · · · , k T from B(i), then yield a set of new position vectors x k 1 , x k 2 , · · · , x k T according to the Equations (20)-(23) by the guidance of Comparing the value of g te x λ i , z with g te x i α λ i , z , g te x i β λ i , z and g te x i δ λ i , z , x ∈ PS(i), then update x i α , x i β and x i δ with the three best position vectors of all.
Step 2.3) Update of z. For each j = 1, 2, · · · , m, if f j x i α < z j , set z j = f j x i α .
Step 2.4) Update of neighborhood. For each if the number of vectors exceeds the EP capacity, the kth nearest neighbor method is used as a truncation strategy.
If the vectors in EP are dominated by f x i α , remove from EP.

Step 3) Stopping criterion:
If the stopping criteria is satisfied, stop running and output EP. Otherwise, return to Step 2. Similar to the original MOEA/D, the MOGWO/D algorithm (Algorithm 1) also optimizes a number of scalar optimization subproblems simultaneously in one iteration, thus improving the optimization efficiency of the proposed algorithm.

Generate a Set of Uniform Weight Vectors
In Step 1.2 of Algorithm 1, a simplex-lattice design [53] is adopted to generate a set of uniformly distributed weight vectors λ i = λ i 1 , λ i 2 , · · · , λ i m , i ∈ N, and m is the dimensionality of the objective space. For each λ i , H is a predetermined positive integer determined according to the sizes of the problems, so, a total of C m−1 H+m−1 weight vectors are obtained. For each λ i , the Euclidean distance to any weight vector is calculated, defining a set B (i) = {i 1 , i 2 , · · · , i T }, in which λ i 1 , λ i 2 , · · · , λ i T are the indexes of the T closest weight vectors to λ i ; then, B (i) is called neighborhood of λ i (including λ i itself, as λ i is the closest weight vector to itself, of which the Euclidean distance is 0). At the same time, the response of x i to λ i also generates a neighborhood, and each individual in the neighborhood corresponds to each weight vector determined by B(i).

Encoding and Decoding
In the proposed MOGWO/D algorithm, the encoding method is similar to the original GWO, and the initial population is randomly generated from a uniform distribution. All position vectors in the initial population are continuous, but the scheduling solutions to the proposed combinatorial optimization problem are not, so a decoding approach is needed to convert the continuous position vectors to the scheduling solutions.
In the proposed NPPLS-JP problem, each position vector needs to include two pieces of information, a job permutation and a production line sequence, and the production line For convenience of expression, the first N position values are marked as Part 1, which corresponds to the job permutation, and the last N position values are marked as Part 2, which corresponds to the production line sequence.
Part 1 and Part 2 are decoded independently. The decoding methods for Part 1 and Part 2 are different because the machining operations of the jobs of some types have machine eligibility constraints. The selection of the production line involves the decoded information of Part 1, that is, the decoding process of Part 2 depends on the obtained jobs' permutation. To discretize the continuous position vectors, the ranked-order value (ROV) rule [54] is used in the decoding processes of Part 1 and Part 2.
For each position value in Part 1 of the position vector, the ROV rule is used to generate ROVs according to the position values in ascending order. If identical position values exist, the ROVs increase from left to right, and then the ROV permutation is obtained. Then, the N 1 smallest values are picked and all are assigned a value of V 1 . V 1 is the order number of order 1, and N 1 is the size of order 1; similarly, ROVs (N 1 + 1) to N 2 are picked and assigned V 2 . V 2 is the order number of order 2, and N 2 is the size of order 2. In the same way, the job permutation is obtained.
For Part 2, first, as in Part 1, the ROV sequence is obtained according to the position values in Part 2 in ascending order. The next step is different from Part 1. For each ROV in Part 2, the same ROV in Part 1 is found, the corresponding order number is obtained, and then its job type is determined. Then, the job type in the first column of the given line selection information table is found, and the corresponding number of optional lines in the second column is obtained. Next, the remainder of the current ROV divided by the number of available lines is obtained. Finally, the corresponding production line number in the third column is found according to the calculated remainder, and the production line number is assigned to the position in Part 2 corresponding to current ROV. In the same way, the production line sequence is obtained. This decoding method for Part 2 can prevent increases in the calculation cost due to invalid solutions caused by the selection of unusable production lines.
A simple example, as shown in Table 1, is used to demonstrate the decoding rules. The sizes of the three orders are two, three, and one, and the total number of jobs is 6six. The detailed decoding process for the example is shown in Figure 3.  Table 1, is used to demonstrate the decodin The sizes of the three orders are two, three, and one, and the total number of jobs The detailed decoding process for the example is shown in Figure 3.   Table 2.

Updating the Position Vectors
The Tchebycheff approach is adopted to decompose the input MOP into objective optimization subproblems and to optimize them simultaneousl neighborhood of each subproblem is defined based on the distances between their vectors. Adjacent subproblems have similar approximate solutions, so each subp is optimized, and only the information of neighboring subproblems is used. Fo generation, position vectors corresponding to the subproblems are updated, w is a positive integer smaller than the neighborhood size, and these position vec randomly selected from the neighborhood corresponding to the subproblem.
Each selected position vector of each subproblem is updated, and this inform . For the decoding process, the production line selection information used in Part 2 is given in Table 2.

Updating the Position Vectors
The Tchebycheff approach is adopted to decompose the input MOP into N singleobjective optimization subproblems and to optimize them simultaneously. The neighborhood of each subproblem is defined based on the distances between their weight vectors. Adjacent subproblems have similar approximate solutions, so each subproblem is optimized, and only the information of neighboring subproblems is used. For every generation, T position vectors corresponding to the subproblems are updated, where T is a positive integer smaller than the neighborhood size, and these T position vectors are randomly selected from the neighborhood corresponding to the subproblem.
Each selected position vector of each subproblem is updated, and this information is used in its neighborhood. First, position vectors are updated according to Equations (20)-(23) through the guidance of the current three best position vectors x α , x β , and x δ , and the position vector set PS = {x 1, x 2 , · · · , x T } of the subproblem is obtained. Second, calculating the g te value of every position vector in the GS(i) = PS ∪ x α , x β , x δ corresponding to λ i . Then, x α , x β , and x δ are updated with the best, second best, and third best position vectors in the GS(i) corresponding to the weight vector λ i , where x α is the optimal solution to the current subproblem. After that, the reference point is updated by calculating the objectives with x α ; if f j (x α ) < z j , then z j = f j (x α ). Third, the neighborhood of the current subproblem is updated with respect to each position vector to the weight vectors of the neighborhood; for each j ∈ B(i), if g te x α λ i , z * ≤ g te x j λ j , z * , x j = x α and FV j = f (x α ) simultaneously. The pseudocode for implementing the update process in one iteration is shown in Algorithm 2.  (2); updateZ(x α ); //update the reference point; updateNeighborhood(x α ); updateEP(x α ); } }

Updating the External Population (EP)
After obtaining the best position vector x α of every subproblem in each generation, the EP needs to be updated. If the condition that F(x α ) is not dominated by the individuals in the EP, F(x α ) is added to the EP, and the individuals that are dominated by F(x α ) are removed.
In the search process of the algorithm, excess individuals are added to the EP. Too many nondominated individuals are not of great significance for solving practical problems, but they increase the difficulty of the data analysis. Therefore, a special truncation strategy is used to maximize the retention of nondominated solution characteristics while maintaining the appropriate EP size. In this strategy, the kth nearest neighbor method [36] is used here to evaluate the individuals in the EP, and the calculation approach for the kth nearest neighbor distance is as follows.
where D(i) is the inverse function of the Euclidean distance from individual i to its k-th nearest neighbor, which is used to reflect the density information of the objective space. σ k i is the kth nearest Euclidean distance of individual i, where the value of k is calculated using Equation (27), and |P| and |EP| are the population size and external population size, respectively. The smaller the D(i) value is, the more scattered the solutions are.

Computational Experiments and Results Analysis
NPPLS-JP is a novel multiline scheduling problem derived from a real-world manufacturing workshop; it is an extension of regular flow shop scheduling problems that has no related research. Therefore, there are no benchmarks available for the proposed MOGWO/D algorithm. In this section, test experiments are designed to assess the performance of the proposed MOGWO/D algorithm by comparing the results obtained on the proposed NPPLS-JP problem with those of three other famous multi-objective optimization algorithms, i.e., NSGA-II [37], MOGWO [55], and MOPSO [56], in terms of three metrics. The results illustrate the effectiveness of the proposed MOGWO/D algorithm.

Evaluation Metrics
Different from single objective optimization problems, MOPs involve the simultaneous optimization of multiple conflicting objectives. An improvement of one objective results in the deterioration of another objective, so MOP algorithms usually obtain a set of tradeoff solutions in terms of the desired objectives, namely, nondominated solutions. There is no absolute optimum among these solutions, and fitness functions cannot evaluate their effectiveness, so a set of metrics is needed to evaluate the performance of multi-objective algorithms for solving MOPs. If the obtained Pareto front is closer to the Pareto optimal front, covering the extreme regions as much as possible, and the nondominated solutions are uniformly distributed in the obtained Pareto front, this means that the obtained results have better convergence and distribution effects. In this paper, the following three metrics are used: (1) Generational distance (GD) [57]. The GD is the most common multi-objective indicator for convergence. It is used to calculate the mean Euclidean distance between the obtained Pareto front and the Pareto optimal front. The calculation formula for the GD is as follows.
where d i is the Euclidean distance from point i of the obtained Pareto front to the closest point in the Pareto optimal front, and |OF| is the number of nondominated solutions in the obtained Pareto front; therefore, GD denotes the mean value of the closest distance from each point in the obtained Pareto front to the Pareto optimal front. A smaller GD value indicates that the obtained Pareto front is closer to the Pareto optimal front; namely, the obtained Pareto front has good convergence. When GD equals zero, the obtained Pareto front is located at the Pareto optimal front. (2) Inversed generational distance (IGD) [58]. This metric is a variant of the GD and is a comprehensive performance indicator. This metric represents the mean Euclidean distance from the points in the Pareto optimal front to the obtained Pareto front. The formulation of the IGD is as follows.
where |PF| denotes the number of points in the Pareto optimal front and D i is the Euclidean distance from point i in the Pareto optimal front to the closest point in the obtained Pareto front. A smaller IGD value indicates better convergence and diversity for the obtained Pareto front. In our experiments, the nondominated solutions obtained from all independent runs of the four algorithms on each instance are regarded as the Pareto optimal front of that instance.
(3) Spread (∆) [37]. ∆ is the diversity metric of the multi-objective optimization that can measure the distribution and spread of solutions. ∆ is calculated as follows.
where m is the number of objectives and n is the number of solutions in the obtained Pareto front. d e j is the minimum Euclidean distance from the nondominated solutions in the obtained Pareto front to the extreme point j of the Pareto optimal front, and d i is the Euclidean distance of the closest pairwise points in the obtained Pareto front, and d is the average value of d i . A smaller value of ∆ represents a better distribution and increased diversity. The calculation of ∆ is simple and does not require knowledge of the whole Pareto optimal front, it uses only the extreme objectives of the Pareto optimal front.

Instance Generation
In the ideal situation, the proposed MOGWO/D algorithm is suitable for all kinds of problems that meet the model definition of NPPLS-JP with different numbers of parallel production lines and jobs, different mixing ratios of job types, and different configurations of production lines. Here, by combining the production system and customers' demand data to generate a set of instances, the performance of the proposed MOGWO/D algorithm is tested. Therefore, we evaluated the effectiveness of the proposed algorithm in different production scenarios by varying the order quantity and order capacity. In the experiment in this section, we only considered the case with two parallel production lines, and gave three production scenarios with different configurations, four different numbers of orders, and three different order capacities; thus, a tally of 3 × 4 × 3 = 36 instances were generated by the combination of the three factors.

Type Available Quantity
For each order in the instances, the trapezoidal fuzzy number for the trapezoidal fuzzy due date is generated as follows. First, a uniform distribution U 0.8 × N J , 3 × N J is given, pt s , N o is the number of jobs in order o, S is the number of operations in the job, and pt s is the maximum processing time for each job operation in order o at any available machine, as shown in Table 3. Then, four integers are randomly taken from the uniform distribution U 0.8 × N J , 3 × N J . The last four time points are sorted in ascending order, as required for the trapezoidal fuzzy due date. The earliness/tardiness penalty coefficients of all orders are set to 5 and 15, respectively.
There are 36 instances to be tested using the proposed MOGWO/D algorithm and the three compared algorithms. For experimentation, the four algorithms are all coded in Java, and all instances are run on an HP Pavilion m4 notebook PC with a Windows 10 Professional 64-bit operating System, 8 GB of RAM, and an Intel Core i5 CPU at 2.60 GHz.

Parameter Settings
From the 36 instances generated above, we can see that the jobs in each instance have eight different scales: 20, 40, 60, 80, 120, 160, 180, and 240. For instances with different job sizes, different population sizes and numbers of iterations should be set so as to obtain the best optimization performance for the algorithms; thus, four different population sizes are given in Table 5. For fairness, MOGWO and MOPSO used the same encoding and decoding method as the proposed MOGWO/D algorithm. Because of the different mechanisms, the encoding and decoding methods of NSGA-II are slightly different. In the NSGA-II algorithm, Part 1 and Part 2 use a permutation encoding scheme [59] to obtain independent encodings, both of which are natural number permutations. The decoding scheme is similar to that of the proposed MOGWO/D algorithm; the only difference is that the individuals in the population are discrete, and these discrete values are natural number permutations. They can be regarded as ROVs, and then decoded according to the decoding rules in Section 4. Unlike the other three algorithms, the encodings in NSGA-II do not need to be discretized first.  In addition, in NSGA-II, a binary tournament is used to as a selection operator for Part 1 and Part 2, the partial mapped crossover (PMX) and swap operation are used as the crossover operator of the two parts, and the insert operation is used as the mutation operator of both parts. The other parameter settings of the four algorithms are also shown in Table 5.
Each instance is run 30 times independently for the proposed MOGWO/D algorithm and the other three comparison algorithms. Table 6 presents the means and standard deviations of the three metrics for the MOGWO/D algorithm-NSGA-II, MOGWO, and MOPSO. By comparing the GD, IGD, and Spread values of the four multi-objective algorithms, it can be seen in Table 6 that the proposed MOGWO/D algorithm has better results in the vast majority of instances. Specifically, with regard to the convergence metric GD, the MOGWO/D algorithm achieves the best results for 33 instances; it is inferior to MOGWO on the "2_8_10" and "2_8_20" instances and inferior to MOPSO on "3_4_30", but the differences are very small. For the spread metric, MOGWO/D algorithm achieves the optimal scheme in 34 instances; it is only inferior to MOGWO on "2_8_10" and inferior to NSGA-II on "3_4_30", but still better than MOPSO. Regarding the comprehensive metric IGD, 34 instances obtained the best metrics values with the MOGWO/D algorithm, only instance "2_ 8_10" was slightly better when using MOGWO, and NSGA-II is superior to MOGWO/D algorithm on "3_4_30". Furthermore, the standard deviations achieved for the three metrics for these instances by the MOGWO/D algorithm were better than those of the other three comparison algorithms in the vast majority of instances. It is clear that the MOGWO/D algorithm has a better optimization performance for solving NPPLS-JP problems, and meets the needs to solve such problems. This may be because of the better balance of the MOGWO/D algorithm between exploitation and exploration, as well as the strategy in which the three best solutions x α , x β , and x δ can be obtained from different levels of the grey wolves' social hierarchy. The experimental results show that the proposed algorithm is superior to other multi-objective optimization algorithms in solving the NPPLS-JP problem.

Experimental Results Analysis
To observe the experimental results more intuitively, Figure 5a-j gives the convergence curve of each instance in the instance set, for which the order capacity is 20, and each curve is the best running result according to the comprehensive performance metric IGD over 30 independent runs for the proposed MOGWO/D algorithm and the three other comparison algorithms. It can be seen clearly that the convergence curves of the proposed algorithm are better than those of the comparison algorithms through these curve graphs. Among them, only on instances "3_2_20" and "3_8_20" is the proposed algorithm similar to the comparison algorithms, and the other instances all yield much better convergence curves than those of the comparison algorithms, which proves the effectiveness of the proposed MOGWO/D algorithm in solving the NPPLS-JP problem.
The NPPLS-JP problem proposed in this paper comes from the actual demand of an axle housing machining workshop. This scheduling demand exists widely in a multivariety mixed production environment, so the proposed model and method are of great significance for production practice. This research has a substantive impact on improving the production efficiency of the workshop, and can significantly enhance the production management level of enterprises, so as to increase the market competitiveness.  The NPPLS-JP problem proposed in this paper comes from the actual demand of a axle housing machining workshop. This scheduling demand exists widely in multivariety mixed production environment, so the proposed model and method are o great significance for production practice. This research has a substantive impact o improving the production efficiency of the workshop, and can significantly enhance th production management level of enterprises, so as to increase the market competitiveness (f) 2_4_20; (g) 2_6_20; (h) 2_8_20; (i) 3_2_20; (j) 3_4_20; (k) 3_6_20; (l) 3_8_20.

Conclusions
In this paper, a multi-objective NPPLS-JP derived from the real-life axle housing machining workshop of an axle manufacturer is studied. In the established NPPLS-JP model, the structures of all parallel lines are symmetrical. However, because of the demands of multivariety mixed production, the process capabilities and production capacities of these parallel production lines are asymmetric, and some types of job operations must be processed on the specific lines. This situation greatly affects the production efficiency of the production system and increases the difficulty of scheduling. To make multivariety mixed production more efficient and to maximize the utilization of the production capacity, a jumping process operation is introduced into the proposed model, which is the largest difference relative to the other general parallel production line scheduling problems. In the NPPLS-JP model, the multiline scheduling, multivariety mixed production, machine eligibility constraints, and MOPs are involved, so it is an NP-hard scheduling problem. In view of this model, we propose a hybrid multi-objective optimization algorithm that incorporates the single-objective GWO into the MOEA/D. The basic idea is to compensate for the shortcomings of the original algorithms by the reasonable mixing of several algorithms to balance their exploration and exploitation of abilities. To verify the effectiveness of the proposed algorithm, a set of instances is designed, and comparative experiments are conducted using the MOGWO/D algorithm as well as three other famous multi-objective optimization algorithms. The experimental results demonstrate that the proposed algorithm is superior to the compared algorithms for solving the NPPLS-JP problem. Furthermore, the experiment also proves that algorithm mixing can improve the performance and expand the application field of the constitutive algorithms.
In future research, we could solve the NPPLS-JP problem under the condition of considering the sequence-dependent setup times, or we could design new metaheuristics to solve the NPPLS-JP problem. Another interesting research direction is to explore new problem-specific rules to improve the performance of the MOGWO/D algorithm in terms of solving the NPPLS-JP. We can also focus on utilizing the MOGWO/D algorithm to solve other workshop scheduling problems, such as job shop scheduling problems and regular flow shop scheduling problems.

Conflicts of Interest:
The authors declare that they have no competing interest.

Abbreviations
The mathematical modeling notations are listed as follows. J a set of jobs, J = 1, 2, · · · , N J i, i the index of jobs, i, i ∈ J N L the total number of production lines L a set of production lines, L = {1, 2, · · · , N L } l, l the index of production lines, l, l ∈ L N S the number of operations, which is equals to the number of stages S a set of operations, S = {1, 2, · · · , N S } s the index of operations, which is also the index of stages, s ∈ S s l the jumping process point of production line l, s l = (1, 2, · · · , S − 1) t the index of job types N t the number of job types k the index of machines M l,s the number of machines at stage s of production line l pt k,t,s the processing time of operation s for a job of type t on machine k