A Discrete Multi-Objective Rider Optimization Algorithm for Hybrid Flowshop Scheduling Problem Considering Makespan, Noise and Dust Pollution

Many optimization algorithms have been proposed to solve hybrid flowshop scheduling problem (HFSP). However, with the development of industry and society, labor right and labor safety have become important problem to consider in production scheduling. So the green HFSP considering makespan, noise and dust pollution becomes an urgent problem to be solved. In this paper, the rider optimization algorithm (ROA) is modified into the multi-objective rider optimization algorithm (MOROA) using Pareto archive and neighborhood sorting techniques. The Pareto archive and neighborhood sorting technology make the Pareto optimal solution set of MOROA have higher coverage and more solutions. Then MOROA is discretized into discrete MOROA (DMOROA) to solve the HFSP considering makespan, noise and dust pollution. DMOROA is tested on 10, 30 and 50 jobs HFSP considering makespan, noise and dust pollution. The test results are compared with two multi-objective algorithms to verify the performance of DMOROA. And the test results verify that the DMOROA is superior to the comparison algorithms in search accuracy, number of non-dominated solutions, diversity of solution set and stability. Therefore, DMOROA is effective in solving multi-objective HFSP considering makespan, noise and dust pollution.


I. INTRODUCTION
The hybrid flowshop scheduling problem (HFSP) is a kind of combinatorial optimization problem which integrates the traditional flowshop scheduling and parallel machine scheduling [1]. HFSP is a complex problem in flowshop scheduling and is difficult to be solved. Complex HFSP is not suitable to be solved by precise algorithm, so heuristic algorithm is put forward to solve HFSP [2]- [5]. The goal of HFSP is to promote the production efficiency. Therefore, HFSP usually uses minimum manufacturing time, total completion time, total delay or other criterions about production efficiency to establish the objective functions [6].
The associate editor coordinating the review of this manuscript and approving it for publication was Diego Oliva .
The single object HFSP considering efficiency has been researched widely. Reference [7] applied immune algorithm to solve HFSP of flexible printed circuit board and obtained good effect. But the algorithm was proposed for specific applications and not easy to generalize to other applications. Literature [8] proposed an artificial immune algorithm to minimize the makespan of HFSP, however the stability was poor. Reference [9] presented a tabu search optimization algorithm to deal with the multiprocessor HFSP. The algorithm was easy to fall into local optimum. Literature [10] used ant colony algorithm to solve HFSP. In reference [11], the iterative greedy algorithm and a variable block heuristic algorithm were fused to solve the HFSP with minimum total flow time, yet only a few solutions were given. In reference [12], an improved particle swarm optimization VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ algorithm was developed based on constructing heuristic algorithm to solve HFSP with minimum maximum completion time and total flow time. However the Pareto front end was not uniform enough. The multi-objective HFSP usually employs the economic indexes such as efficiency or energy consumption to establish the objective functions. Reference [13] presented a discrete multi-objective firefly algorithm to solve the HFSP considering makespan, total flow time and machine idle time. But there was no other comparison metrics except for coverage metric in the experiment. In reference [14], a multi-objective discrete artificial bee colony algorithm was proposed to solve the green HFSP with the maximum completion time and the minimum total energy consumption. However the coding method was complex and the efficiency was low. Literature [15] applied a multi-objective gray wolf optimization algorithm based on cellular automata and variable neighborhood search to optimize HFSP considering production efficiency, energy consumption and noise pollution. The experiment was comprehensive, but lacked the complex test problem. Reference [16] enhanced the multi-objective evolutionary algorithm with problem-specific heuristics to deal with the lot-streaming HFSP. This method was a new research direction and close to the reality, yet the algorithm had high time complexity. In literature [17], a discrete multi-objective firework algorithm using cluster analysis and fuzzy correlation entropy analysis was introduced to handle the HFSP with total production cost, maximum completion time, average running time and average idle time. Unfortunately, this research focused on static situation and could not deal with dynamic situation. Reference [18] developed a multi-objective multi-verse optimizer algorithm based on Levy flight and chaos local search to solve reentrant HFSP with minimum completion time, delay and energy consumption. The method had strong ability of machine scheduling, but its performance were not outstanding in the case of continuous scheduling. In literature [19], a multi-objective evolutionary algorithm was introduced to search the optimal solutions of hot-rolling scheduling problem in the compact strip production. Reference [20] employed a dual-objective evolutionary algorithm to solve stochastic hybrid flowshop deteriorating scheduling problem considering makespan and total tardiness.
With the development of economic globalization and the spread of people-oriented concept, more and more countries pay attention to labor right and labor safety, but the loss caused by occupational disease and injurie is still huge every year [21]- [23]. Occupational diseases mainly refer to the diseases caused by workers' exposure to dust, noise or radioactive substance in their occupational activities. According to a report released by the international labor organization in 2019 [24], about 6500 people die of occupational disease every day in the world, and the economic loss caused by occupational safety and health problem accounts for 4% of the global GDP [25], [26]. The report suggested that countries could use new technologies to reduce occupational risks and create new jobs through sustainable development and green economy to continuously improve occupational safety.
There are two modes of labor right protection. One mode is to protect labor right legally through state legislative intervention or regulation [27]. The other mode relies on the joint participation of labor and capital in enterprise decision-making to implement enterprise autonomy and labor self-protection [28], [29]. Due to the different levels of economic and social development in different countries, the levels of labor protection in national law are uneven. Therefore, labor protection needs a long time and a lot of financial investment, and it is not easy to implement. In this study, new decision-making model and optimization technique are applied to production planning and scheduling to improve the protection of labor right without any change of production process and equipment. This paper is devoted to solving the green HFSP considering makespan, noise and dust pollution.
Rider optimization algorithm (ROA) is a new heuristic algorithm proposed by D. Binu and B. S kariyapa in 2019. ROA optimizes the problem by assuming that a group of riders drive towards the goal to compete for the championship [30]. ROA has unique structure and advantage in dealing with complex problem. The ROA is modified into the multi-objective rider optimization algorithm (MOROA) using Pareto archive and neighborhood sorting techniques. The Pareto archive and neighborhood sorting technology make the Pareto optimal solution set of MOROA have higher coverage and more solutions. Then MOROA is discretized into discrete MOROA (DMOROA) to solve the HFSP considering makespan, noise and dust pollution. DMOROA and HFSP considering labor safety are both novel issues, and this paper can provide an enlightening idea for these fields.
The rest of this paper is arranged as follows. The section II gives the definition of HFSP considering makespan, noise and dust pollution. The section III describes the ROA and gives the pseudo code. The section IV describes how to change ROA to MOROA and gives the pseudo code. In section V, MOROA will be discretized into discrete MOROA and DMOROA is used to solve the HFSP considering makespan, noise and dust pollution. In the section VI, the simulation experiment is carried out, and the experiment results are analyzed according to the evaluation metrics. The section VII summarizes the whole paper and discusses the future work.

II. DEFINITION OF HFSP CONSIDERING MAKESPAN, NOISE AND DUST POLLUTION
The definition of HFSP considering makespan, noise and dust pollution is given by model hypothesis and mathematical equations [31]- [34]. The schematic diagram of HFSP is shown in Figure 1.

A. MODEL ASSUMPTIONS
The model assumptions for HFSP considering makespan, noise and dust pollution are described below.
(1) The n jobs are processed in s working stages hypothetically. There are m parallel machines in each stage. The speed of processing job, unit pollution of noise and dust may be different in each stage.
(2) In each stage, each job can only be processed by one machine, and each job must complete all stages. Moreover the machine of a certain stage cannot handle operation of other stages.
(3) All machines and jobs are available at 0 time, and jobs follow the principle of first come first process.
(4) The buffer area between two consecutive stages is infinite. The preparation time and the transmission time are included in the processing time.
(5) Each machine can only process one job at one time. The work that is carried out must be completed without interruption.
(6) If at least one machine is available in a stage, the stage shall be started immediately after the completion of the previous stage.

B. DESCRIPTION OF VARIABLES AND SYMBOLS
Variables and symbols of HFSP considering makespan, noise and dust pollution are described as follows.

C. ESTABLISHMENT OF MULTI-OBJECTIVE OPTIMIZATION MODEL
The mathematical model of HFSP considering makespan, noise and dust pollution is expressed as follows.
Finally, the model of HFSP considering makespan, noise and dust pollution is established as expression (4).

III. RIDER OPTIMIZATION ALGORITHM
Optimization algorithm can be roughly divided into two categories: natural heuristic algorithm and artificial computation algorithm [35], [36]. Natural heuristic algorithm seeks solution to optimization problem based on inspiration from biological or physical phenomena [37]. Artificial computation algorithm is a guiding way to solve complex problem based on human thought. In 2019, D. Binu and B. S kariyapa proposed a new computing method that is called virtual computing [38], [39]. Virtual computing follows an idea based on imagination to solve optimization problem. At the same time, they put forward ROA based on the idea of virtual computing. The ROA assumes that there are four groups of riders with the same number of riders moving towards the same goal to compete for the championship. The four groups of riders are the bypass riders, followers, overtakers and VOLUME 8, 2020 attackers. The bypass riders take the lead by making a detour. The followers follow the leader most of the time. The overtakers overtake by judging the leaders and their positions. The attackers use the fastest speed to rush to the target point. In the ROA, each rider adjusts his own direction, gear, accelerator, brake and other parameters in real time to update his position. At the end of the competition (usually to reach the maximum number of iterations), the leader wins. The mathematical representation of the ROA is described below. The schematic diagram of ROA is shown in Figure 2.
is the position vector of the i-th rider at time t. B, F, O and A respectively represent the number of four groups of riders: bypass riders, followers, overtakers and attackers.

B. PARAMETER INITIALIZATION
The ROA includes the following parameters or parameter matrices: direction, gear, accelerator, brake, success rate. Equation (11) gives the direction angle of the rider at time t. Equations (12), (13) and (14) give the initialization mode at time 0.
At time t, the gear of the i-th rider is selected from {0,1,2,3,4}, and the initial gear value of each rider is set to 0. The formula is equation (15).
The value range of accelerator and decelerator of each rider belongs to [0,1]. Equation (16) represents accelerator, which is initialized to 0. Equation (17) represents decelerator, which is initialized to 1.
The maximum speed limit of each gear is shown as expression (18).
E represents the number of gears. V i max is the maximum speed that the i-th rider can reach. X i U and X i L are the upper and lower bounds of the rider position values, respectively. T OFF is the end time of the competition or the maximum number of iterations.
The success rate is defined as the reciprocal of the distance from the current position to the leader position. X i is the current position of the i-th rider. L T is the leader position. The rider individual with the best fitness value is the leader. The equation of success rate is equation (20).

C. LOCATION UPDATE EQUATIONS 1) BYPASS RIDER LOCATION UPDATE
Bypass rider location is updated by equation (21).
where δ and β belong to the random numbers of 0-1 uniform distribution, η and ξ belong to the random integers of 1-R uniform distribution.

2) FOLLOWER LOCATION UPDATE
Follower location is updated by expression (22).
where X L (L, K ) represents the k-th dimension of the leader position. T t i,k represents the k-th dimension turning angle of the i-th rider at t time. d t i represents the distance that the rider will travel next. d t i is described by equation (23).
v t i refers to the speed of the i-th rider at time t, which depends on the rider parameters. V i max is the maximum speed that the i-th rider can reach. E t i , e t i , K t i and V E i respectively represent the gear, accelerator, brake and gear speed limit of the rider.
The dimension k that is determined during location update depends on the real-time probability P t ON . Equation (25) is used to calculate the real-time probability P t ON . P t ON is multiplied by the corresponding dimension index j. Then j * P t ON is compared with Q to determine the coordinate that is considered for the change of position. And this comparison process is described in equation (26).

3) OVERTAKER LOCATION UPDATE
Equation (27) represents the update equation of overtaker position, and the updated dimensions are selected according to equations (30) and (31). The difference between the leader and the overtaker is calculated according to equation (30), and the dimensions less than the average value are selected. Then, the direction of each rider is calculated according to the relative success rate of equation (29). Finally, equation (27) is employed to update the overtaker position according to the direction of travel, the position of the leading driver, the position of the current driver and the selected dimension.
D I t (i) is the direction indicator of the i-th rider at time t. D I t (i) is calculated using the relative success rate S R t (i). The relative success rate is the ratio of the success rate of the corresponding rider to the maximum success rate of rider. The relative success rate is calculated by equation (29), where r t (i) is the success rate of the i-th rider at time t.
In equation (30),l(i, j) is the distance between the current rider position and the leader position in the j-th coordinate. The value of l(i, j) that is smaller than its average value is selected from j according to equation (31). In the equation (31), µ i is the average value of the standardized distance vector of the i-th rider.

4) ATTACKER LOCATION UPDATE
The attacker follows a similar update process as the follower. However, the attacker will update all the values in the coordinate instead of only the selected values. The equation (32) is the update equation of attacker position.

D. RIDER PARAMETER UPDATE
In addition to the rider parameters shown in the above, here add a new parameter: the activity counter. The direction and gear are updated according to the activity counter. The activity counter is updated based on the success rate. The mathematical representation of the activity counter is given in equation (33).
The activity counter value of the i-th rider at time t+ 1 is A t+1 c (i), which depends on r t (i) and r i+1 (i). If r i+1 (i) > r t (i) then the A t+1 c (i) value is 1, otherwise 0. Equation (34) updates the direction matrix with the value of the activity counter.
Equation (35) updates the gear according to the value of the activity counter. Equations (36) and (37) determine the values of the accelerator and the brake based on the newly selected gear.

E. THE PSEUDO CODE OF RIDER OPTIMIZATION ALGORITHM
The pseudo code of ROA is described in algorithm 1.

IV. MULTI-OBJECTIVE RIDER OPTIMIZATION ALGORITHM
ROA was originally developed to solve the continuous singleobjective optimization problem. It cannot be directly used to solve the multi-objective optimization problem. The Pareto archive and neighborhood sorting technique are applied to change ROA into MOROA.

A. DEFINITION OF MULTI-OBJECTIVE OPTIMIZATION PROBLEM
In single-objective optimization, because the objective function is unitary, the solutions are easy to compare. For the minimization problem, if x < y, then the solution x is better than y. However, the multi-objective optimization problem has multiple objective values. The solutions in the solution space cannot be compared by using common relational operators. At this time, the concept of Pareto optimality is given: when a solution shows better or equal objective values on all objectives, and provides a better objective value in at least VOLUME 8, 2020 Algorithm 1 Rider Optimization Algorithm (1). Input: Random positions of the riders X t (2). Output: Leader X L (3). Begin (4). Initialize the population (5).
Initialize the rider parameters: Steering angle T, Gear E, Accelerator e, and Break K (6).
Find the success rate r t (7).
Rank the riders based on r t (14).
Select the rider having the maximum r t as the leader (15).
(1) Pareto dominance. Suppose there are two vectors x 1 , x 2 ∈ X . If the vector x 1 dominates vector x 2 , that is (x 1 x 2 ), then it can be expressed by equation (38).
(2) When the solution vector x * is Pareto optimality, it can be expressed as equation (39).
(3) A group of Pareto optimal solutions consist Pareto optimal set, the Pareto optimal set is expressed by equation (40).
(4) The set of objective function values corresponding to Pareto optimal solution set is called Pareto optimal front, which is expressed as equation (41).
When solving multi-objective problem, the goal of algorithm is to find a set of Pareto optimal solutions. This set of Pareto optimal solutions needs to be closest to the real Pareto front, large number and even distribution.

B. PARETO ARCHIVE AND NEIGHBORHOOD SORTING TECHNIQUE 1) PARETO ARCHIVE
The multi-objective heuristic algorithm will give a set of Pareto optimal solutions when the maximum number of iterations is reached. If the Pareto optimal solution set has a higher coverage rate and more solution number as far as possible, the Pareto optimal solution set can try to approach the real Pareto front end. MOROA applies Pareto archive technique to maintain the Pareto optimal solution set in the optimization process.
Pareto archive is an independent multi-dimensional solution space with predefined size. The non-dominated solutions generated on each iteration will be added to the Pareto archive according to certain rules. The non-dominated solutions will be carefully maintained in Pareto archive. The non-dominated solutions will be added to Pareto archive or delete from Pareto archive followed a series of rules, and these rules are listed below.
(1) Addition of non-dominated solutions. After the new non-dominated solutions are produced in the current iteration, they should be added to the Pareto archive for maintenance.
(2) Deletion of non-dominated solutions. The size of the Pareto archive is limited, and the solutions in the Pareto archive should be dominated by each other. In the following cases, the deletion program of non-dominated solutions is executed.
(a) After the new non-dominated solutions are added, the solutions in the Pareto archive no longer dominate each other. When the non-dominated solutions in the updated population are added to the Pareto archive, the solutions in the Pareto archive may no longer dominate each other. At this time, it is necessary to remove the no-dominated solutions from the archive.
(b) After the above procedures are executed, the number of solutions in the Pareto archive exceeds the upper limit of the Pareto archive. The excess non-dominated solutions need to be deleted. It is necessary to select the deleted solutions reasonably to ensure that the obtained Pareto front end is uniform and highly approximate to the real Pareto front end. Therefore MOROA uses neighborhood sorting technique to select the deleted solutions and ensure to get optimal Pareto solution set.

2) NEIGHBORHOOD SORTING TECHNIQUE
In order to make the Pareto front as evenly distributed as possible, Neighborhood sorting technique is used to delete the non-dominated solutions that exceed the upper limit of Pareto archive. And neighborhood sorting technique is applied to select the leaders and guide the population position update. Neighborhood sorting technique calculates the density around a non-dominated solution. The density of the solution is employed as the probability to determine whether the solution is deleted or selected as the leader.
Neighborhood sorting technique defines a hypersphere for each solution in the archive. The radius of the hypersphere is determined by the prior knowledge of the problem. The prior knowledge is related to the background of the practical problem. Each hypersphere will contain other solutions (0, 1 or multiple solutions) around the solution. The probability of each hypersphere is equation (42).
The c is a constant greater than 1, N i is the number of non-dominated solutions in the i-th hypersphere. The roulette mechanism is used to choose the solutions to be deleted, and the interval of roulette will be determine by probability P i . The equation (42) means that the algorithm has a higher probability to select deleted solutions from areas with higher density. The purpose of this choice is to make the non-dominated solutions distribute evenly in the Pareto front end. Therefore, the solutions exceeding the Pareto archive capacity will be deleted according to roulette mechanism and equation (42).
The neighborhood sorting technique will be used to choose the leader to guide the Pareto front end more evenly. Similarly, the hypersphere for each solution will be defined and the probability of each hypersphere is expressed by equation (43).
where c is a constant greater than 1, and N i is the number of solutions in the i-th hypersphere. In the roulette mechanism, the hypersphere with the lowest density assigns the maximum probability to be selected as the leader.  Figure 3 shows an example of the above two selection processes. The red circle represents the hypersphere with the highest density, and the black circle represents the hypersphere with the lowest density. The corresponding solutions of these two hyperspheres have the maximum probability of being deleted or selected as the leader respectively.

C. PSEUDO CODE OF MOROA
The pseudo code of MOROA is explained in algorithm 2.
Calculate the r t of riders (17).

V. DMOROA FOR HFSP CONSIDERING MAKESPAN, NOISE AND DUST POLLUTION
In literature [30], ROA has been employed to diagnose analog circuit faults, and proven to be effective in handling complex engineering problem. ROA has the advantages of fast search speed, high search accuracy and good at dealing with complex problem. Therefore, MOROA will be discretized into DMOROA. DMOROA is used to solve the HFSP considering makespan, noise and dust pollution. The data of industrial production is difficult to obtain and lacks of universality, so the experimental data is randomly generated. The randomly generated data is usually used to test HFSP [15]- [17]. Therefore, the value of matrix S, C, N, D, T is determined by random initialization in this paper. The initialization of matrix S and C should satisfy the conditions of equations (6) and (7).
Each HFSP is identified as ''number of jobs -number of stages -number of parallel machines in each stage''. Suppose that n jobs and s stages are needed, and m parallel machines  can be selected in each stage, then the HFSP can be named as ''n-s-m'' problem. A specific ''3-3-2'' problem will be applied as an example to illustrate the process of population initialization and discretization in DMOROA. The time sequence Gantt chart of problem ''3-3-2'' is Figure 4.

A. POPULATION INITIALIZATION
The search space of ''n-s-m'' problem is composed of n + n * s * m dimensions. The vector X = (x 1 , x 2 , . . . ,x n+n * s * m ), x i ∈ [0, 1] (i = 1,. . . ,n + n * s * m) represents the value of the job processing sequence and machine selection. For example, the search space of the ''3-3-2'' problem has 21 dimensions, and the vector X of the ''3-3-2'' problem is initialized and shown in Table 1.
The first n dimensions of vector X represent job processing sequence and form vector Q = (x 1 , x 2 , . . . , x n ). The n + 1 to n + n * s * m dimensions of vector X form machine selection matrix M . The 1-3 dimensions of the ''3-3-2'' problem form vector Q and list in the first row of Table2. The 4-21 dimensions of the ''3-3-2'' problem form matrix M and represent in Table 3.

B. POPULATION DISCRETIZATION
The discretization method of reference [13] is used to discretize the population. The components of vector Q are numbered in order from small to large. After numbering, the vector Q is discretized into vector Q . For ''3-3-2'' problem, the first row of Table 2 is converted to the second row of Table 2. The processing sequence of the jobs 1, 2, 3 is 1, 2, 3.
The smallest position value rule [27] is used to determine machine selection.  Table 3 represent that job n1 selects machines m1 or m2 in stage s1. The minimum value is marked as 1 and the other values are marked as 0, so (0.26744, 0.33730) becomes (1, 0). According to the smallest position value rule, matrix M is discretized into matrix Z . For ''3-3-2'' problem, the minimum position rule discretizes Table 3 into Table 4.

C. SOLVING HFSP CONSIDERING MAKESPAN, NOISE AND DUST POLLUTION WITH DMOROA
When MOROA is used to solve the HFSP, the initial population will be discretized to transformed MOROA into DMOROA. The DMOROA will be employed to solve HFSP considering the makespan, noise and dust pollution. The three objective functions are determined by equations (1), (2) and (3)     The parameters of equation (3) are determined by matrix D and T . For problem ''3-3-2'', the matrix T , N and D are Table 5, 6 and 7. DMOROA solves the ''3-3-2'' problem and obtain the optimal solutions that are listed in Table 8.
The principle of DMOROA for HFSP considering makespan, noise and dust pollution is described in algorithm 3.

VI. SIMULATION EXPERIMENT
The experiment results of HFSP considering makespan, noise and dust pollution are compared between DMOROA, non-dominated sorting genetic algorithm 2 (NSGA2) and multi-objective particle swarm optimization algorithm (MOPSO) to evaluate the performance of DMOROA.

A. EXPERIMENTAL PLATFORM
The operating environment of the simulation experiment will be illustrated. The machine model is Dawning 5000A. The server is configured with Xeon x5620 CPU (4 cores) 2, 24GB memory and 1 * 300 GB SAS hard disk. The operating system is Linux rhel5.6. The programming tool is Matlab 2012a (for Linux).

B. EXPERIMENTAL METHOD
In order to evaluate the performance of DMOROA, the different configuration data from the size and structure of the HFSP considering makespan, noise and dust pollution is needed. In addition, relevant performance evaluation metrics are needed and the level and scope of some evaluation metrics are shown in Table 9.
In the experiment, each algorithm uses the same coding strategy that is explained in Part V, and the parameter settings of each algorithm are shown in Table 10. Each algorithm will test 10 times for each problem and run 500 iterations each time. The results are compared with the evaluation metrics proposed in the next section.

C. EVALUATION CRITERIONS
For multi-objective HFSP, it is difficult to choose a single performance metric to compare the effectiveness and efficiency of different algorithms. Therefore, four performance criterions are proposed for convergence, diversity, distribution uniformity and stability. The four performance criterions are coverage metric, number of non-dominated solutions, spread and mean relative deviation index. VOLUME 8, 2020  (4). Initialize the population X (5). X is divided into job processing sequence vector Q and machine selection matrix M (6). vector Q is discretized to vector Q (7). Matrix M is discretized to matrix Z according to the smallest position value rule (8). Initialize the rider parameters: Steering angle T, Gear E, Accelerator e, and Break K (9). Initialize the Archive, find the leader (10). Find the success rate r t (11). while t< T OFF (12). for i=1 to R (13).
Calculate the r t of riders (21).
Return Archive (23). t=t+1 (24). end for (25). end while (26). Terminate  C(B, A). This means that it is necessary to compare the values of C(A, B) and C(B, A) to determine which non-dominated solution set has higher precision.

2) NUMBER OF NON-DOMINATED SOLUTIONS
Number of non-dominated solutions represents the number of non-dominated solutions obtained by the algorithm that is applied to test a HFSP. In multi-objective HFSP, producers often need to choose specific scheduling strategies according to the results of scheduling evaluation. Getting more Pareto optimal solutions can make the Pareto front more uniform and smooth, which is beneficial for producer to choose the most suitable scheduling strategy for the current situation.

3) SPREAD
The spread metric measures the diversity and the distribution uniformity of solution set. The spread is a more comprehensive metric to show the quality of solution set. The calculation equation of spread is equation (45).
where d i is the Euclidean distance between the i-th solution in one solution set and its closest solution in its neighbor,d is the average value of all d i . d e j is the Euclidean distance between the extreme solutions in the j-th objective and the boundary solutions of the obtained solution set. The obj_num is the number of objective.

4) MEAN RELATIVE DEVIATION INDEX
Mean relative deviation index is applied to measure the stability of the algorithm. The mean relative deviation index of each objective is calculated by equation (46).
In the equation (45), MRDI j represents the mean relative deviation index of the j-th target. obj * j represents the optimal 88536 VOLUME 8, 2020  value of the j-th target. obj i,j represents the j-th objective function value of the i-th non-dominated solution.

D. RESULTS AND DISCUSSION
The results of experiments will be analyzed to demonstrate the performance of DMOROA. Figure 5 to Figure 7 show the Pareto front end obtained by each algorithm on each problem. Figure 8 to Figure 10 are the comparison of the target values of each problem on the two-dimensional coordinate figure. Figure 5 and Figure 8 show the 10 workpiece problems, Figure 6 and Figure 9 show the 30 workpiece problems, Figure 7 and Figure 10 show the 50 workpiece problems. In Figure 5, the solutions of DMOROA obviously dominate the solutions of the comparison algorithms in ''10-3-3'', ''10-3-6'', ''10-3-9'' and ''10-5-6''. In the ''10-5-3'' problem, the solutions of DMOROA are slightly better than NSGA2, but the coverage of Pareto front end is wider and the diversity of solutions is better. In the 10-5-9 problem, a part of DMOROA solutions are excellent, which dominate the solutions of all comparison algorithms, but the other solutions are mixed with the solutions of MOPSO and dominated by the solutions of NSGA2.
The test results have been compared and analyzed, indicating that DMOROA has the following characteristics.
(1) For the multi-objective HFSP, in general, the higher complexity of the problem has, the higher search precision of DMOROA has.
(2) Compared with classical multi-objective heuristic algorithm, DMOROA can show excellent performance in solving the HFSP that have the same complexity.
(3) If DMOROA can obtain more non-dominated solutions, Pareto front of DMOROA has more diversity. But with the increase of the complexity of the HFSP, the number of non-dominated solutions of DMOROA will decrease.
(4) When HFSP has low complexity, DMOROA can obtained the uniform and diverse solution set easily. With the increase of the HFSP complexity, the solution set of DMOROA decreases the uniformity slightly.

VII. CONCLUSION
This paper studies green HFSP from the perspective of labor right, and proposes a HFSP considering makespan, noise and dust pollution. The MOROA is presented based on ROA, and MOROA is discretized into DMOROA. MOROA not only keeps the advantages of ROA, but also improves the ability of handling complex problem. DMOROA is introduced to solve HFSP considering makespan, noise and dust pollution.
The simulation experiment is carried out with 10, 30 and 50 jobs HFSP considering the makespan, noise and dust pollution between DMOROA, NSGA2 and MOPSO. The experimental results show that DMOROA is superior to NSGA2 and MOPSO in search accuracy, number of non-dominated solution and stability. Therefore DMOROA is suitable to solve HFSP considering makespan, noise and dust pollution. When the HFSP considering makespan, noise and dust pollution increases complexity, the number of non-dominated solutions and diversity of solution set are slightly reduced. In the future research, this problem is worth studying and solving. The green scheduling factors considered in this paper are limited. The factors such as energy efficiency and radioactive substance can be introduced as the optimization objectives in the further work [47].