Bi-objective optimization of identical parallel machine scheduling with flexible maintenance and job release times

This paper investigates an identical parallel machine scheduling problem with flexible maintenance and job release times and attempts to optimize two objectives: the minimization of the makespan and total tardiness simultaneously. A mixed-integer programming (MIP) model for solving small- scale instances is presented first, and then a modified NSGA- Ⅱ (M-NSGA- Ⅱ ) algorithm is constructed for solving medium- and large-scale instances by incorporating several strategies. These strategies include: ( ⅰ ) the proposal of a decoding method based on dynamic programming, ( ⅱ ) the design of dynamic probability crossover and mutation operators, and ( ⅲ ) the presentation of neighborhood search method. The parameters of the proposed algorithm are optimized by the Taguchi method. Three scales of problems, including 52 instances, are generated to compare the performance of different optimization methods. The computational results demonstrate that the M-NSGA- Ⅱ algorithm obviously outperforms the original NSGA- Ⅱ algorithm when solving medium- and large-scale instances, although the time taken to solve the instances is slightly longer. that without neighbourhood search. These results fully show that the use of neighbourhood search can improve the search quality of the solution, improve the distribution and quality of the solution, and requires a longer time to obtain a better solution.


Introduction
Parallel machine scheduling (PMS) is an important branch of production scheduling (Cheng & Sin, 1990). In the literature on PMS, most studies assume that machines are available at all times (Liu et al., 2020;Kim et al., 2019). However, in real-world manufacturing, machines need to be maintained and hence may become unavailable during a certain period, which makes it highly important to jointly schedule production and maintenance tasks to simultaneously improve system availability and system throughput. The problem of scheduling jobs with maintenance reduces to the problem often referred to in the literature as scheduling with machine availability constraints. Lee and Chen (2000) first researched a parallel machine scheduling problem (PMSP) with maintenance in which each machine must be maintained once in the planning horizon, and they proposed branch and bound (B&B) algorithms based on the column generation approach. To date, many studies on various PMS with different preventive maintenance (PM) strategies have been carried out. Unfortunately, most research focuses on single-objective optimization, especially the makespan, and few studies focus on multiple objective optimization. Berrichi et al. (2009) researched the PMSP with machine reliability and simultaneously minimized the makespan and machine system unavailability. Based on this research, Wang and Liu (2015) proposed a multi-objective PMSP with two kinds of resources (machines and moulds) and flexible PM activities on resources, to simultaneously minimize the makespan, the machine system unavailability, and the mould system unavailability. With the development of intelligent manufacturing, customer satisfaction plays an increasingly vital role in both the manufacturing and service industries. On-time delivery is becoming increasingly important. Therefore, the bi-objective to simultaneously minimize the makespan and total tardiness is studied in this paper.
In addition, most studies on the joint optimization of production scheduling and PM usually assume that jobs reach zero time.
However, in real-world manufacturing, it is very common for jobs to arrive at different times. Many PMSP studies consider job release times without considering integration optimization with PM (Kramer et al., 2020;Abedi et al., 2015;Nessah et al., 2008;Yalaoui & Chu, 2006). Cui and Lu (2017) studied the single machine scheduling problem with machine flexible maintenance and job release times. Chen et al. (2021) studied two identical parallel machines with flexible maintenance to minimize the makespan without considering the job's release time. To the best of our knowledge, the PMSP with PM and job release times has not been documented in the literature. In this research, we extended the number of machines from references (Cui & Lu, 2017;Chen et al., 2021) to multiple machines and considered job release times to simultaneously minimize the makespan and total tardiness. Our objective is to determine the assignment of jobs to the parallel machines, the decision of the maintenance activities, and the optimal sequence of jobs and PM on each machine. The MIP model and an M-NSGA-II are successfully customized to solve the considered problem.
The rest of the paper is organized as follows: the next section surveys the related literature. A problem description is presented in Section 3. The MIP model is proposed in Section 4. Section 5 provides details about the modified NSGA-II. The experimental results are demonstrated in Section 6, and finally, Section 7 includes conclusions and discusses future research.

PMSP with preventive maintenance
PMSPs incorporated with PM can be divided into different models according to different criteria. The basic criterion is the number of maintenance tasks, including single maintenance (Xu et al., 2009;Wang & Wei, 2011) or multiple maintenance (Hashemian et al., 2014;Li et al., 2017;Xu et al., 2008). In addition, maintenance may be performed on one of the parallel machines or on all parallel machines. Some earlier studies assumed that maintenance is only performed on some parallel machines (Tan et al., 2011) or only once within the planning horizon (Lee & Chen, 2008;Yoo & Lee, 2016). Recent research has generally assumed that all parallel machines have multiple maintenance, and the related research is divided into three kinds: fixed interval maintenance, flexible maintenance, and variable maintenance in the form of maintenance activities. Fixed interval maintenance involves periodically performing fixed time maintenance according to a fixed time interval. Yoo and Lee (2016) studied the PMSP in which each machine requires maintenance activity once within a given time window. Li et al. (2017) studied the PMSP in which each machine is subject to periodic maintenance. Xu et al. (2009) considered two PMSPs in which one machine is periodically unavailable. Flexible maintenance is fixed time maintenance within a flexible time window or upper threshold. Xu et al. (2008) studied a PMSP with ɛ-almost periodic maintenance activities to minimize the makespan. Sun and Li (2010) studied a similar problem and assumed that the largest consecutive processing time for each machine cannot exceed an upper limit T. Chen et al. (2021) studied two identical PMSPs with flexible maintenance to minimize the makespan. Variable maintenance is a kind of maintenance in which the maintenance interval and time depend on the condition of the machines and/or jobs. Wang and Wei (2011) studied PMSPs with deteriorating maintenance activity, that is, when delaying maintenance increases the time required to perform it. Wu et al. (2020) considered the dispatchingdependent deterioration of machines and machine-health-dependent production rates and proposed a dynamic dispatching and PM model to minimize the weighted long-run average waiting costs of MTO systems. Moradi and Zandieh (2010) studied the joint production and maintenance problem of parallel machines, and they concluded that the availability of machines is related to the failure rate and repair rate.

Multi-objective optimization of PMSP
Instead of a single objective, multiple but conflicting objectives are often considered when a manager executes production scheduling. Thus, a Pareto front consisting of a set of non-dominated solutions is required so that the manager can choose one of the alternatives from the set while preparing the production plan. However, compared with a single objective, such as minimizing the makespan, multi-objective optimization of PMSPs is rarely studied. According to the number of objectives, this problem can be divided into two types: bi-objective and multi-objective. For bi-objective optimization, in the context of green manufacturing, minimizing the energy or power consumption and cost is an important objective; therefore, the objective is to minimize the power and the makespan Anghinolfi et al., 2021) as well as the tardiness penalty and power cost (Fang & Lin, 2013). In addition, the bi-objective to minimize the makespan and the total weighted earliness and tardiness of jobs (Abedi et al., 2015) as well as the makespan and the system unavailability (Berrichi et al., 2009), Moradi and Zandieh (2010)) can also be found in the literature. For multiple objective optimizations, based on the literature (Berrichi et al., 2009), the bi-objective is extended to three; Wang and Liu (2015) proposed a multi-objective PMSP with two kinds of resources (machines and moulds) and with flexible preventive maintenance activities on resources, and presented an integrated optimization method with NSGA-II adaptation. Bandyopadhyay and Bhattacharya (2013) proposed a modified version of NSGA-II with a new mutation algorithm for a PMSP with three objectives. Cochran et al. (2003) proposed a two-stage multipopulation genetic algorithm to solve PMSPs with multiple objectives.

Solution methodology of the PMSP
The solutions to PMSP generally include two decisions: one assigns jobs to machines, and the second sorts of the jobs assigned to the machine. These two decisions can be made in parallel, i.e., allocation first and then sequencing, or sequencing can be performed during assignment. The solution method for PMSP in the literature mainly includes the exact method, rule-based heuristics, the meta-heuristic algorithm and the intelligent evolutionary algorithm.
The exact method can obtain the optimal solution to the problem, but it is only suitable for small-scale instances. The exact methods proposed for solving the PMSP are mainly the B&B algorithm (Nessah et al., 2008;Yalaoui & Chu, 2006) and mathematical models (Wu & Wang, 2018;Naderi & Roshanaei, 2020). To improve the efficiency of the exact models, the lower bounds and the dominance properties have been incorporated into the B&B algorithm. In addition, Alhadi et al. (2020) presented a polynomial-time approximation scheme to generate an approximate Pareto Frontier that minimizes the maximum lateness and makespan. Rule-based heuristics are widely used in the real-world manufacturing industry (Lee, 2018). The rules proposed for the problem Pm// include the LPT, MULTIFIT, COMBINE, and LISTFIT (Gupta & Ruiz-Torres, 2001); those proposed for the problem Pm// ∑ include the EDD, SPT, TPI, MDD, KPM, Minimum Slack and BHG (Biskup et al., 2008); and those proposed for the problem Pm// ∑ are the QB and QBP procedures (Schaller & Valente, 2018). Considering the job's release time, Akturk and Ozdemir (2001) proposed some rules, namely, the ATC (apparent tardiness cost), X-RM (X-dispatch ATC), KZRM (Kanet and Zhou approach to ATC), and COVERT (weighted cost over time), for the problem 1/r / ∑ . Meta-heuristics and intelligent evolutionary algorithms are the focus of the existing research. Kayvanfar et al. (2017) proposed an intelligent water drop algorithm for PMSP with controllable processing times. Dipak and Gupta (2018) proposed an improved cuckoo search algorithm to minimize the makespan for identical PMSP. Abdeljaoued et al. (2020) proposed a simulated annealing meta-heuristic for PMS under resource constraints.

Problem statement
This paper investigates the problem that a set of jobs = { , ⋯ , ⋯ , } is to be scheduled on a set of identical machines = { , ⋯ , ⋯ , }, and the objectives are to simultaneously minimize the makespan and total tardiness. The problem can be denoted as Pm/ , nr, FPM/( , ∑ ), where "Pm" represents the identical parallel machines, " " represents jobs released at different times, "nr" represents jobs that are non-resumable, and "FPM" represents flexible preventive maintenance, which means that the continuous processing time of the machines cannot exceed the maintenance threshold UT. Additionally, we assume that ≤ for all jobs. The maintenance time is denoted by . The Gantt chart of a schedule for this problem is shown in Fig. 1. As shown in Fig. 1 . is the PM of machine , and [ ] is the makespan of machine . Proof. The sequencing problem of each machine in Pm/ , nr, FPM/( , ∑ ) is equivalent to the bin-packing problem. Because both the bin-packing problem and the P2// problem have been proven to be NP-Hard, the problem P2/nr, FPM/ is NP-Hard. In addition, the problem 1/ / ∑ is NP-Hard (Nessah et al. (2008)); thus, Pm/ / ∑ is NP-Hard. The problem Pm/ , nr, FPM/( , ∑ ) extends the number of identical parallel machines from 2 to m and considers the constraint of the job's release time and the objective of minimizing total tardiness; thus, problem Pm/ , nr, FPM/( , ∑ ) is NP-Hard.
Definition: Full loading (FL) means that the cumulative processing time of the jobs during two consecutive PMs, or the sum processing time of the jobs in a batch, is equal to the maintenance threshold UT. When = , batch is fully loaded; otherwise, batch is not fully loaded, and the waste time is equal to − .

MIP model
Based on the following assumptions, an MIP model for the problem studied in this paper is established: (1) each job has a release time; (2) the processing time, release time and due date of the job are known in advance; (3) the maintenance threshold and maintenance time are known in advance; (4) the time allocated to any machine for processing the same job is the same; (5) each machine can only process one job at any given time and shall maintain continuous processing, which shall not be pre-empted and interrupted by other jobs; and (6) there is no job importance, that is, the weight is not considered. The input parameters and decision variables are introduced before developing the proposed model.

Input Parameters
set of jobs set of machines , job index from set machine index from set the processing time of job the release time of job the due date of job the completion time of job the start time at the position of machine the cumulative processing time at the position of machine the completion time of machine in position k threshold of the cumulative processing time that the machine can process continuously time required for maintenance a sufficiently large positive integer, it is safe to set = ∑ + * Decision variables equal to 1 if job is processed at the position of machine , 0 otherwise equal to 1 if maintenance is conducted after the position of machine , 0 otherwise

Mathematical model
. (1) ≥ + * , ∀ = 1,2,3, … , = 2,3, … , In this model, the objectives are to minimize and ∑ as shown in Eq. (1) and Eq. (2). Constraints (3) and (4) ensure that each job can only be processed at one position of one machine and that each position of one machine can only be occupied by one job. Constraint (5) specifies that no PM is required after the last position of each machine. Constraints (6) and (7) limit the start time of each machine at each position. Constraint (8) defines the completion time of each machine at each position. Constraint (9) describes the cumulative processing time of each machine at the first position. Constraint (10) defines the cumulative processing time of each machine at each position. Constraint (11) ensures that the cumulative processing time of each machine at each position is not larger than UT. Constraint (12) ensures that the job is arranged in a more forward position. Constraint (13) defines the maximum completion time. Constraint (14) defines the completion time of each job. Constraint (15) defines the tardiness of each job.

M-NSGA-Ⅱ algorithm
For the proposed problem, an M-NSGA-II algorithm is presented by incorporating three strategies. (ⅰ) A decoding method based on dynamic programming (DP) is constructed. (ⅱ) Dynamic probability crossover and mutation are designed. (ⅲ) Neighbourhood search method is proposed. The general operation of the M-NSGA-II algorithm is shown in Algorithm 1.
Algorithm 1:M-NSGA-Ⅱ algorithm 1 Create an initial population of size using a random rule 2 Decode the individual using the method based on DP, sort the population based on the fast-non-dominated sorting method 3 Create offspring population by applying the dynamic probability crossover and mutation operators as well as elitist operators 4 While stopping criteria are not verified, do Create population R = ∪ of size 2 and construct the different of R using the fastnon-dominated sorting procedure Put of according to the crowding procedure Renew population with neighbourhood search method Create offspring population by selection, crossover, and mutation with dynamic probability End While

Chromosome representation
According to the decision of the proposed problem, the presentation of the chromosome is a vector with length n + m − 1, where n is the number of jobs, and m is the number of machines. An example is [ ] where the jobs allocated to different machines are distinguished by 0, and the string of numbers divided by 0 represents the sequence of the jobs processed on one machine. [ ] , although the job allocation and sequence on each machine have been determined according to the position of 0 in the chromosome, the maintenance or the job-grouping batch decision still needs to be made. For the identical PMSP without considering the job's release time to minimize , a near-optimal solution can be obtained by a maintenance decision made based on the FL principle as shown in Appendix 1. However, it is not true for the same problem considering the job release times. An example with five jobs on a single machine is used to explain, and the jobs data is shown in Table 1. Assume the solution π = {J , J , J , J , J }, UT = 10, and mt = 2. The Gantt chart of the solution obtained based on the FL and non-FL principle is shown in Fig. 2. The solution obtained without considering the FL is better than that obtained with the FL. The reason is that machine have to be idle for a period to be fully loaded. For the identical PMSP considering the job release times, the maintenance decision is a dynamic decision, and we propose a decoding method based on dynamic programming (DP) process (Zhou et al., 2018)). In each stage, only the different job-grouping batch schemes of a newly added job and the currently sorted job are compared, and the best optimal solution is selected as the final decoding result.

Neighbourhood Search operators
For the identical PMSP, it is easy to generate redundant solutions in the population because of the symmetry of the machines. For example, the solution = { , , , 0, , , } is equal to the solution = { , , , 0, , , }. Too many redundant solutions will lead to premature convergence of the algorithm. To improve the diversity of the population, we firstly need to identify the redundant solution, then to renew it with a neighbourhood search method. If chromosome and in the population is different, then the difference between the two chromosomes is ( , ) ≻ 0, otherwise ( , ) = 0. If the difference between chromosome and other chromosomes in the population is equal to 0, the solution corresponding to chromosome is redundant solution. Because each chromosome is composed of m segments, firstly the lengths of the m segments are compared. If they are different, the chromosomes are different; otherwise, the difference is determined by comparing the genes. The Eq. (19) and Eq. (20) respectively give the calculation method of the difference of each segment and the difference of the whole chromosome between two chromosomes and with the same length.
where and is gene segments on chromosomes and , respectively; is the length of the gene segment, and and is the locus in the gene segment. Let represent the number of redundant solutions in the population, and represents the density of the redundant solutions; then, = / . If the is larger than the threshold, a neighbourhood search method is used to renew the individuals in the redundant solutions population. If the new solution is better than the original solution, it replaces the original solution; otherwise, the original solution is kept. Because the two objectives are optimized, two kinds of four neighbourhood search operators are designed. The average makespan ̅ of the maximum makespan MaxC and the minimum makespan , the average ∑ of the maximum total tardiness M ∑ , and the minimum total tardiness ∑ in the redundant solution's population are calculated firstly, then the neighbourhood search operators are used.

Neighbourhood search guided by makespan improvement
If the objective of the individual satisfies > ̅ , one of the two M-neighbourhood search methods is randomly selected. The purpose of the M-neighbourhood search method is to improve the of the solution, including M-Insert and M-Swap. M-Insert: Select a batch with waste time p from machine with the minimum completion time, search for a job with a processing time equal to or less than in machine with the maximum completion time, and insert it into the selected batch on machine .
M-Swap: The job with a longer processing time selected from machine is swapped with the job with a shorter processing time selected from machine . The difference in processing time between the swapped jobs should be less than the difference in the completion time of the two machines.

Neighbourhood search guided by total tardiness improvement
If the objective of the individual satisfies ∑ > ∑ , one of the two T-neighbourhood search methods is randomly selected. The purpose of the T-neighbourhood search method is to improve the ∑ of the solution, including T-Insert and T-Swap.
T-Insert: Select the job with the shortest processing time from among the delayed jobs and insert them one by one into the position before the job with a larger due date time. T-Swap. Select the job with the longest processing time from among the delayed jobs and swap its position with the job whose subsequent due date is smaller than the job one by one.

Crossover operator
A fixed crossover probability = 0.5 is set for the NSGA-II algorithm, and a dynamic crossover probability p = 1.5 − × p is designed for the M-NSGA-Ⅱ algorithm. g represents the number of iterations at present, and G represents the total number of iterations to be performed.
A two-point crossover method is adopted. A random chromosome in the first Pareto front F of the population after fast non-dominated sorting and a random chromosome in the other front are selected as parents. A random gene in [1, + − 1] is selected from each parent for crossover, and the other genes are obtained by mapping. If the offspring is better than the parent , it will be retained. Otherwise, the original parent is retained. An example of crossover is shown in Fig. 3.  Fig. 3. Example of the crossover operation

Mutation operator
A fixed mutation probability = 0.05 is set for the NSGA-Ⅱ algorithm, and a dynamic mutation probability p = × p is designed for the M-NSGA-Ⅱ algorithm. Two methods, MO1 and MO2, are used for the mutation operation of the first 50% generation population and the last 50% generation population, respectively. MO1 randomly selects an individual, selects a gene fragment of length ⌊ / ⌋, and reverses the sequence of jobs in this gene fragment. MO2 selects an individual and randomly selects two genes to exchange. If the randomly selected gene is 0, reselection is required. An example of the mutation operation is shown in Fig. 4.

Parameter tuning and computational experiments
The NSGA-Ⅱ and the M-NSGA-Ⅱ algorithms were written using the Dev C++ language, and the MIP model was constructed on the IBM ILOG CPLEX Optimization Studio Ver. 12.7.1 platform. The three algorithms' simulation experiments were all run on an Intel(R) Xeon(R) E5-2650 v4 @2.20 GHz CPU and 64 GB RAM workstation.

Performance measures
This paper compares different algorithms in terms of operation efficiency and effectiveness. The operation efficiency is directly expressed by the operating time. Four measures are used to compare the operation effectiveness or different aspects of the non-dominated fronts obtained by MIP, the NSGA-Ⅱ algorithm and the M-NSGA-Ⅱ algorithm, including the number of Pareto solutions , the C metric, the inverted generational distance (IGD) and the distance ∆ between the different solutions. The C metric is used to show differences between two Pareto fronts (Wang & Liu, 2015). The ( , ) states the percentage of a solution in B dominated by at least one solution of A. It is calculated by Eq. (21): where < means that a dominates b or b is dominated by a. It can be seen from Equation (21) where is the Euclidean distance between solutions and , and Ω is the approximate Pareto front, and is acquired by merging the Pareto solutions obtained by all algorithms. The smaller the ( , Ω) is, the better of the algorithm.
∆ is used in the literature , and it is defined by Eq. (23): where is the Euclidean distance between two consecutive solutions in the obtained non-dominated set A, and ̅ is the average of all distances . The smaller the ∆ is, the better of the algorithm.

Parameter tuning
Parameters largely affect the performance of the M-NSGA-Ⅱ algorithm. Therefore, it is necessary to tune the parameters using appropriate methods. The Taguchi method is a strong stochastic technique from among the design of experiments (DOE) methods. It aims to improve quality and reduce cost when designing a process or production, and its main idea is to combine different factors and their degrees using vertical arrays and factorial designs. Then, the analysis is conducted with fewer numbers of experiments . The parameters that affect the performance of the M-NSGA-II algorithm include the population number Popsize, the number of iterations Gensize, the crossover rate p , and the mutation rate p . The levels of these four parameters were determined according to the preliminary experiments as shown in Table 2. A set of 13 typical problem instances with three scales (small, medium, and large) were generated to perform the parameter tuning experiments, and a decoding method based on FL was adopted. The signal-to-noise (S/N) ratio of the parameters was analysed to determine the parameters of the M-NSGA-II algorithm. Using Popsize, MaxGen, p and p as control factors where each factor has 3 levels, the number of orthogonal experiments is L (3 ), the objective function is used as the response variable, and the number of experiments is 13×9×10=1170. The results for the small-scale problem instances are shown in Table 3. Because there are two objectives, the same data normalization method used in reference (Yue et al., 2019) is applied. Through normalization of the S/N at different levels of different response variables, the comprehensive S/N ratio at the final parameter level was obtained, and Fig. 5 shows the mean S/N for the small-scale problems. Due to the two objectives of the minimization functions, the level of the parameter at which the minimum value of the evaluation criteria is obtained is preferred. The parameters of NSGA-II for different sets of problem instances were also selected using the Taguchi method. The optimum level of parameters for each algorithm for each instance set is illustrated in Table 4.   2019)). The experiments were conducted using the machine parameters, job parameters and PM activity parameters. The machine parameters include the number of machines ; the job parameters include the number of jobs , processing time , release time and due date ; and the PM activity parameters include the threshold and maintenance time . The instances and the range of experimental parameter are shown in Table 5.

Comparison results for different decoding methods
To validate the performance of the decoding method based on DP, we solved the problem instances using the DP-based and FL-based decoding methods. For each instance, the method was run 10 times by considering randomness, and the average values of the performance metrics were recorded. The fronts obtained by the two methods mentioned above were combined, and all the non-dominated solutions were used to form set Ω to compute the IGD. The results are shown in Table 6.  The results show that the DP-based decoding method performs better than the FL-based decoding method in terms of the C metric and IGD, especially for large-scale problem instances composed of many jobs. There is no significant difference between the two decoding methods in terms of the ∆ metric. The Nd metric obtained by the FL-based decoding method is slightly larger than that obtained by the DP-based decoding method. This result shows that the solution set obtained by the DP-based decoding method is more diversity and better approximates the Pareto front. In this paper, the Nd metric obtained by the NSGA-II and M-NSGA-II algorithms is limited since minimizing the makespan indirectly leads to minimizing the total tardiness, which is why the Nd metric of the FL-based decoding method slightly outperforms the DP-based decoding method. The DP-based decoding method searches for the job-grouping solution from the first position to the current position when adding a job at each stage to make maintenance decisions to obtain the best solution, which also limits the number of Pareto solutions. Thus, the DP-based decoding method requires more time than the FL-based decoding method.

Comparison results for different neighbourhood search methods
To validate the performance of the proposed neighbourhood search method, the M-NSGA-II algorithm with and without the neighbourhood search method was used to solve the instances. The results are shown in Fig. 6. The results show that the M-NSGA-II method with neighbourhood search can obtain better non-dominated solutions compared to the method without neighbourhood search in terms of the C metric and IGD. This advantage becomes more evident as the size of the problem increases. However, there is no significant difference between the two methods in terms of the ∆ metric and N . Since neighbourhood search adds further exploration and searching of the solution space, the time required is longer than that without neighbourhood search. These results fully show that the use of neighbourhood search can improve the search quality of the solution, improve the distribution and quality of the solution, and requires a longer time to obtain a better solution.

Comparison results for different algorithms
To compare the performance of MIP, M-NSGA-II, and NSGA-II, all three methods were used to solve the small-scale instances, and the M-NSGA-II and NSGA-II algorithms were used to solve the medium-and large-scale instances. The results are shown in Table 7 and Table 8.  As seen in Table 7, the MIP method can obtain the optimal Pareto solution for all 12 instances, although the time needed for solving the problem greatly increased with the increase in the size of the problem. The NSGA-Ⅱ algorithm slightly outperforms the M-NSGA-II algorithm in terms of the IGD, the C metric and N . There is no significant difference between the two methods in terms of ∆ . The NSGA-Ⅱ algorithm is more efficient than the MIP and M-NSGA-II algorithms, and this advantage becomes more evident as the size of the problem increases. It can be seen from Table 8 that the M-NSGA-Ⅱ algorithm outperforms the NSGA-Ⅱ algorithm in terms of the IGD and C metrics for medium-and large-scale problem instances, and this advantage also becomes more evident as the size of the problem increases. However, there is no significant difference between the two methods in terms of N and ∆ . The efficiency of the NSGA-II algorithm is better than that of M-NSGA-II because M-NSGA-II uses the DP-based decoding method and the neighbourhood search method.

Influence of PM parameters
To judge the influence of the PM parameters on the NSGA-Ⅱ and M-NSGA-Ⅱ algorithms, their four metrics were analysed using variance analysis (the confidence interval was 95%), and the results are shown in Table 9.
As shown in Table 9, for the different UT and , all the indicators for the M-NSGA-II algorithm are not significantly different. In contrast, for the NSGA-II algorithm, the N indicator for small instances and the IGD indicator for medium instances are significantly different for parameter . This result shows that the M-NSGA-Ⅱ algorithm is stable and effective for all instances.

Conclusions
In this paper, an identical PMSP with machine flexible maintenance and job release times has been investigated to minimize the makespan and total tardiness simultaneously. A MIP model was first established to obtain the exact Pareto fronts for smallscale instances. To tackle the medium-and large-scale instances, an M-NSGA-Ⅱ algorithm was constructed with three strategies. After tuning the M-NSGA-Ⅱ parameters using the DOE and the Taguchi method, computational experiments were conducted. The results show that the M-NSGA-Ⅱ algorithm obviously outperforms the original NSGA-Ⅱ algorithm when solving medium-and large-scale instances, although the time required to solve these instances is slightly longer.
In future research, more realistic parallel machine manufacturing environments such as unrelated parallel machines should be investigated. The maintenance period depends on the machine's reliability, and the objective is to minimize the energy and power costs. Furthermore, considering the development of Big Data and Intelligent Manufacturing, parallel machine dynamic scheduling and intelligent scheduling methods based on reinforcement learning should be developed.