**1. Introduction**

The flexible job scheduling problem (FJSP) is an extension of the classic job scheduling problem (JSP) and is closer to the actual production environment. In the scheduling process of the FJSP, processing operations can be processed on all optional machines. The assignable machine expands the search range of feasible solutions and also increases the complexity and the difficulty of solving the problem. The FJSP is a complex NP-hard problem, and its solution time increases exponentially as the problem size increases.

With the increasingly prominent energy crisis and environmental pollution, manufacturing has gradually become one of the hot spots in modern manufacturing. Manufacturing has adopted a new sustainable manufacturing model that has attracted widespread attention from industry and academia. In the manufacturing process of an enterprise, workshop scheduling is an important factor in the manufacturing process. It not only affects the production efficiency and the economic benefits of the enterprise but also is closely related to the social responsibility of the enterprise. Therefore, it has important theoretical and practical significance to conducting research on flexible job scheduling with the goal of protecting the environment and saving energy.

In the past few decades, the single-objective FJSP (SO-FJSP), which has been extensively studied in the literature, has usually sought to minimize the total completion time [1–6]. However, many realistic scheduling problems often need to optimize multiple

**Citation:** Sun, X.; Wang, Y.; Kang, H.; Shen, Y.; Chen, Q.; Wang, D. Modified Multi-Crossover Operator NSGA-III for Solving Low Carbon Flexible Job Shop Scheduling Problem. *Processes* **2021**, *9*, 62. https://doi.org/10.3390/ pr9010062


Received: 4 December 2020 Accepted: 26 December 2020 Published: 29 December 2020

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

objectives at the same time, and these objectives usually conflict with each other. The main method used to solve the MO-FJSP is the MOEA, which can be roughly divided into two categories: the weighting method and the Pareto method. The weighting method solves the MO-FJSP by assigning different weights to each objective and transforming the multi-objective problem into a single-objective problem. The Pareto method solves the MO-FJSP based on the Pareto dominance relationship and generates a set of Pareto optimal solutions. As a Pareto method, the non-dominated sorting genetic algorithm-II (NSGA-II) [7] is an effective method to solve various multi-objective optimization problems in recent years. Z.-Q. Jiang et al. [8] used the NSGA-II algorithm, which optimizes mutation strategies to solve the multi-objective FJSP of strategy. Yuan Y. and Xu H. [9] proposed a new memetic algorithm that combines the memetic algorithm with the NSGA-II to solve the FJSP with the goal of minimizing completion time, total workload, and critical workload. Bandyopadhyay and Bhattacharyaput [10] proposed a modified NSGA-II with a new mutation algorithm for a parallel machine scheduling problem and proved the effectiveness of the algorithm. The research goal of all the improved algorithms is to solve the MO-FJSP more effectively.

The FJSP is widely present in various manufacturing processes. It has received extensive attention from researchers. A large number of research results have appeared [8–19]. However, in these studies, the objective function of the problem is rarely to minimize carbon emissions or total energy consumption. The FJSP with these objectives has not attracted attention. The existing FJSP research focuses on the relationship between carbon emissions or energy consumption and time [16–19]. The machine load is rarely optimized as a target problem. The completion time and the machine load are also two conflicting issues. The price of minimizing the total completion time is the long-term overload of high-performance machines. Therefore, it is necessary to optimize machine load as an objective.

This paper establishes an MO-FJSP model targeting carbon emissions, the completion time, and the machine load and modifies NSGA-III [20]. By studying the differences in exploration and developmental capabilities of different NSGA-III variants in the MO-FJSP decision space, indicator-based thought is introduced into NSGA-III, a genetic model of multiple populations and multiple crossover operators is established, and a new evolutionary mechanism is proposed. We apply this evolutionary mechanism to NSGA-III and propose the co-evolutionary NSGA-III (NSGA-III-COE). Then, calculation experiments are carried out on 27 well-known FJSP benchmark instances [21,22]. Quantities of experiments in this paper prove that the NSGA-III-COE achieves good results in solving the low carbon MO-FJSP and verifies the advantages and the competitiveness of the NSGA-III-COE in solving the low carbon MO-FJSP.

#### **2. Mathematical Modeling of the MO-FJSP**

Before building the mathematical model and the assumptions are listed below.

#### *2.1. Assumptions*

The machining process satisfies the following assumptions and constraints. These assumptions and constraints are common in the FJSP literature [5–12].


#### *2.2. Mathematical Model*

The mathematical formulas and constraints are as follows:

$$T = \min T\_{\max} = \min \left\{ \sum\_{p=1}^{N\_n} \left\{ T\_p \right\} \right\},\tag{1}$$

$$\mathcal{W} = \min \mathcal{W}\_{\max} = \min \left\{ \sum\_{h=1}^{N\_m} \{ \mathcal{W}\_h \} \right\}, \tag{2}$$

$$\mathcal{C} = \min \mathcal{C}\_{\text{max}} = \min \left\{ \sum\_{p=1}^{N\_{\text{tr}}} \sum\_{h=1}^{N\_{\text{tr}}} \sum\_{q=1}^{N\_p} \left( T\_{\text{sh}} \mathcal{C}\_{\text{sh}} + T\_{pq\text{h}} \mathcal{C}\_{pq\text{h}} \right) \right\}. \tag{3}$$

$$\mathcal{S}\_{pq}, \mathcal{T}\_{pqh} \ge 0, \mathcal{J}\_p \in \mathcal{J}; q = 1, 2, \dots, \mathcal{N}\_p; h = 1, 2, \dots, \mathcal{N}\_m \tag{4}$$

$$\sum\_{h=1}^{M\_{p\eta}} \sigma\_{p\eta h} = 1, l\_p \in f; q = 1, 2, \dots \cdot N\_{p\prime} \tag{5}$$

$$\mathcal{S}\_{p(q+1)} \ge \mathcal{S}\_{pq} + T\_{pq\text{h}},\\I\_p \in \mathcal{J}; \mathcal{M}\_{\text{li}} \in \mathcal{M}\_{pq}; q = 1, 2, \dots, N\_p - 1,\tag{6}$$

$$\sum\_{p=1}^{N\_{\rm tr}} \sum\_{q=1}^{N\_{\rm p}} T\_{p\rm qph} \sigma\_{p\rm qlr} \le \mathcal{W}\_{\rm h}, f\_p \in f; M\_{\rm h} \in \mathcal{M}\_{p\rm q}.\tag{7}$$

Objective (1) represents the objective function that minimizes the maximum completion time, objective (2) represents the objective function that minimizes the total machine load, and objective (3) represents the objective function that minimizes total carbon emissions during processing. Constraint (4) indicates that the start time and the processing time of the process are greater than or equal to 0, constraint (5) indicates that each process can only select one machine from the set of candidate processing machines, constraint (6) indicates that each job must be processed in the given order, and constraint (7) represents the total machine load.

#### *2.3. Chromosome Encoding*

The FJSP needs to select processing machines for each process and sort the processes allocated on each machine. According to the characteristics of the FJSP, this paper adopts the two-dimensional chromosome encoding method based on the combination of process coding and machine coding [9]. The following uses an FJSP instance to illustrate chromosome encoding. The processing time of an FJSP instance with three jobs, three machines, and seven processes is shown in Table 1.


**Table 1.** Machine processing schedule of the flexible job shop scheduling (FJSP) instance. '-' means that the process cannot be processed by this machine.

Table 2 is a set of randomly generated chromosome codes corresponding to the instances in Table 1. Pro is a process-based code used to determine the processing order, and Mac is a machine-based code used to allocate the processing machines for each process. Each column of genes can be interpreted as *Opq*,*Mh* , that is, Pro is the processing sequence

of the process *O*<sup>21</sup> → *O*<sup>11</sup> → *O*<sup>12</sup> → *O*<sup>22</sup> → *O*<sup>31</sup> → *O*<sup>13</sup> → *O*<sup>32</sup> , and the corresponding Mac is the machine on which the process is processed (*M*2, *M*1, *M*2, *M*2, *M*3, *M*3, *M*3). The gene (2, 2) in the first column of Table 2 can be interpreted as (*O*21, *M*2), and the gene (2, 2) in the fourth column can be interpreted as (*O*22, *M*2). That is, the first process of the second job is processed on machine 2, and the required processing time is *T*<sup>212</sup> = 1; the second process of the first job is processed on machine 2, and the required processing time is *T*<sup>222</sup> = 1 (the processing time is found in Table 2).

**Table 2.** A set of chromosome code representation instances.


#### *2.4. Chromosome Decoding*

Chromosome decoding allocates a period of time for each operation on the designated machine according to the sequence of the process in the chromosome. Take the FJSP instance shown in Table 1 as an instance to decode the chromosomes in Table 2. There are two different decoding schemes; the first is to assign the machine processing according to the sequence of the process chromosome Pro to the Mac chromosome. The scheduling Gantt chart is shown in Figure 1a. The second is when processing a process in the Pro chromosome, first obtain the machine selected in the Mac chromosome for the process, and then scan the machine from left to right to determine the idle time interval between the processing processes and insert the current process until the time period that can be processed is found. The second scheduling Gantt chart is shown in Figure 1b. The second decoding scheme allows process scheduling to search for the earliest available idle time interval on a specified machine, which can effectively reduce the production cycle. Therefore, this paper uses the second decoding scheme.

**Figure 1.** Chromosome decoding corresponds to the Gantt chart.

#### **3. Modification of the NSGA-III**

#### *3.1. Introduction to the NSGA-III*

The NSGA-III replaces the crowding distance selection operation in the NSGA-II with a reference point-based selection operation and uses well-distributed reference points to maintain the diversity of the population. This is the reason why this paper selects it. In addition, the NSGA-III is the most widely used MOEA in the existing literature. Next, we briefly describe the main procedures of the NSGA-III.

The NSGA-III first defines a set of reference points, then randomly produces an initial population containing *N* individuals, and then iterates until the termination condition is met. *Pi* is the population in the *i*th generation, and *Qi* is the population generated by *Pi* after the reproduction phase. In order to select *N* individuals from population *Ri*(*Ri* = *Pi* ∪ *Qi*) into the next generation, non-dominated sorting is used to divide the individuals in *Ri* into several different non-dominated layers (*F*1, *F*2, ···) and to add the non-dominated

layers into *Si* in order. *Si* is determined to be selected. It is the population of the (*i* + 1)th generation *Pi*+1. Assuming that *Fl* is the last non-dominated layer where the population size of *Si* is larger than N for the first time, use the reference point to find the optimal number of remaining *Pi*+<sup>1</sup> individuals in *Fl* and join the next generation population *Pi*+1.

#### *3.2. Study of NSGA-III Variants*

The original NSGA-III used simulated binary crossover (SBX) [23] to generate individual offspring. This section calls this method the NSGA-III-SBX to distinguish it from other NSGA-III variants studied in this section. In this paper, we introduce cycle crossover (CX) [24], order-based crossover (OBX) [25], order-crossover (OX) [26], partially mapped crossover (PMX) [27], and position-based crossover (PBX) [25] into the NSGA-III, replacing the original SBX operator and forming 5 NSGA-III variants, namely, NSGA-III-CX, NSGA-III-OBX, NSGA-III-OX, NSGA-III-PBX, and NSGA-III-PMX.

In order to study the search performance of all NSGA-III variant algorithms in the decision space, we use part of three sets of well-known FJSP benchmark instances, including ka3, ka4, and ka5 in the Kacem instance [22] and mk4, mk5, and mk7 in the BRdata instance [21], to conduct exploration and testing. These six instances are representative from simple to complex. In the experiments in this section, we use these six benchmark instances to explore the NSGA-III variants mentioned in this section and propose a modified NSGA-III based on the research results.

Table 3 lists the parameters used in this section to study the different variants of NSGA-III, and we use uniform parameter values for all variants. In preliminary research, we found that the widely used MOEAs are prone to fall into local optimum on FJSP. Some studies in the literature have found that a larger mutation probability can effectively help the population jump out of the local optimum [10]. Thus, we use a high mutation probability. In order to explore whether the performance of different variants is related to population size, we use two population sizes in our research, 200 and 300. When the number of iterations reaches the set maximum number of iterations, the algorithm is terminated. To ensure a fair comparison, for each benchmark instance, all variants are run independently with the same initial population 30 times, and the average of 30 experiments is taken for comparison.

**Table 3.** Parameter setting of the non-dominated sorting genetic algorithm-III (NSGA-III) variant algorithm.


In our experiments, we use the generational distance (GD) [28] and the inverted generational distance (IGD) [29] as evaluation indicators to evaluate the convergence of the non-dominated solution set and the comprehensive performance of the algorithm. They can be expressed as follows.

*GD*: Assuming that *P* is the solution set obtained by the algorithm and *P*∗ is a set of uniformly distributed reference points sampled from the Pareto front (PF), the *GD* of solution set *P* is defined as follows:

$$GD\left(P, P^\*\right) = \frac{1}{|P|} \sqrt{\sum\_{y \in P} \min\_{\mathbf{x} \in P^\*} \left(\text{dis}\left(\mathbf{x}, y\right)^2\right)}\,\tag{8}$$

*dis*(*x*, *y*) represents the Euclidean distance between point *y* in solution set *P* and point *x* in reference set *P*∗. *GD* only evaluates the convergence of the solution set. The smaller the *GD* value is, the better the convergence of the algorithm is.

*IGD*: Assuming that *P* is the solution set obtained by the algorithm and *P*∗ is a set of uniformly distributed reference points sampled from the Pareto front (PF), then the IGD value of solution set *P* is defined as follows:

$$IGD(P, P^\*) = \frac{1}{|P^\*|} \sum\_{x \in P^\*} \min\_{y \in P} dis(x, y)\_{\prime} \tag{9}$$

*dis*(*x*, *y*) represents the Euclidean distance between point *x* in reference set *P*∗ and point *y* in solution set *P*. If |*P*∗| is large enough to fully represent the Pareto front, then the *IGD* can comprehensively measure the convergence and the diversity of the solution set. If we want to obtain a smaller *IGD*, the solution set must be close enough to the Pareto front in the target space.

When calculating the *GD* and the *IGD*, a reference set is needed. Since the actual Pareto front of the benchmark instance is unknown, the reference set used in the calculation of the *GD* and the *IGD* in this paper is formed by collecting all the non-dominated solutions found during the runtime of all implemented algorithms.

Next, we compare the search behavior of different NSGA-III variants on the MO-FJSP. Simply put, we want to understand the algorithm search process for solution sets in the decision space as well as which algorithm is better at exploration and which algorithm is better at development. Knowing these can help us design more effective evolutionary mechanisms. First, we randomly initialize a population for each instance and then perform 200 and 300 iterations.

Figures 2 and 3 depicts the trajectories of the GDs obtained from the six NSGA-III variants as the number of iterations increases when the population sizes are 200 and 300, respectively. Figures 2 and 3 show that, although the population sizes are different, the same NSGA-III variant shows very similar convergence trends on these benchmark instances, and the convergence of different NSGA-III variants is significantly different due to the different complexities of the benchmark instances. As the complexity of the benchmark instance increases, NSGA-III-CX, NSGA-III-OBX, and NSGA-III-PBX show better convergence, and the GDs of these three algorithms decrease in a similar way during the evolutionary process. This phenomenon indicates that the initial population is a randomly distributed solution in the decision space, and then the optimal solution is searched continuously. NSGA-III-CX, NSGA-III-OBX, and NSGA-III-PBX can search for better solutions faster than other algorithms. The above experimental results show that the three NSGA-III variants, NSGA-III-CX, NSGA-III-OBX, and NSGA-III-PBX, can explore more optimal solutions more effectively in the decision space.

In order to study the development capabilities of different NSGA-III variants, we conduct a similar experiment. The difference from the previous experiment is that we replace the initial population with a population that is already closer to the Pareto front, which can be obtained through iteration by any multi-objective evolutionary algorithm. In this experiment, only the three NSGA-III variants with better exploration capabilities are used. The purpose is to study the abilities of NSGA-III-CX, NSGA-III-OBX, and NSGA-III-PBX to develop better solutions. The three NSGA-III variant algorithms use the population close to the Pareto front as the initial population to perform 100 iterations. As before, use two population sizes, 200 and 300.

Figures 4 and 5 depicts the trajectories of IGDs obtained from the three NSGA-III variants as the number of iterations increases when the initial population is close to the Pareto front when the population size is 200 and 300 respectively. We compare Figures 4 and 5 first. Similar to the previous experiment, the same NSGA-III variant showed very similar ability to develop better solutions when the population size was different. However, the situation in Figures 4 and 5 is very different from that in Figures 2 and 3. In Figures 4 and 5, from the beginning, as the number of iterations increases, a certain algorithm will reduce IGD faster, while the IGD of other algorithms will decrease more slowly. Since the initial population is a population closer to the Pareto frontier, an algorithm with a faster IGD decline has a better ability to develop better solutions in the decision space. In Figures 4e and 5e, on the

instance mk5, the NSGA-III-CX has the best ability to develop better solutions. NSGA-III-OBX has the best development capability on other instances. It is speculated from this that when faced with different decision spaces, the crossover operator with better development capabilities may change. The research in this section can help us design more effective evolutionary mechanisms to solve low carbon FJSP.

**Figure 2.** At population size of 200, the evolutionary trajectory of the generational distances (GDs) of the NSGA-III variants on the six FJSP instances (the average of the results of 30 independent runs; the initial population is a random population).

**Figure 3.** At population size of 300, the evolutionary trajectory of the GDs of the NSGA-III variants on the six FJSP instances (the average of the results of 30 independent runs; the initial population is a random population).

**Figure 4.** At population size of 200, the evolutionary trajectory of the inverted generational distance (IGD) for the NSGA-III variants on the six FJSP instances (the average of the results of 30 independent runs; the initial population is the population close to the Pareto front in the target space).

**Figure 5.** At population size of 300, the evolutionary trajectory of the IGD for the NSGA-III variants on the six FJSP instances (the average of the results of 30 independent runs; the initial population is the population close to the Pareto front in the target space).

#### *3.3. The NSGA-III-COE Proposal*

The research in the previous two sections shows that the NSGA-III variant using the three crossover operators CX, OBX, and PBX has better exploration capabilities than others in the decision space. When the initial population is close to the Pareto frontier, in most instances, NSGA-III-OBX has the best ability to develop better solutions in the decision space. However, in a few instances, NSGA-III-OBX does not have the best development capabilities. Therefore, which cross-operator has the best development capability is still uncertain. The above research aims to help us design a more effective evolutionary mechanism when solving the MO-FJSP.

Many studies have shown that exploring and developing strategies at the same time can find more useful information from the decision space in the process of finding a better solution. If we make full use of the three crossover operators of CX, OBX, and PBX, we can expect the algorithm to achieve better performance. This is the motivation for proposing the NSGA-III-COE algorithm. The effects of three different crossover operators are naturally integrated to improve the search ability of the decision space and maintain the diversity of the population. This is the main purpose of the NSGA-III-COE.

The NSGA-III-COE is the result of the combination of Pareto dominance and indicatorbased thought. In order to achieve our purpose, we decided to coevolve three subpopulations using CX, OBX, and PBX crossover operators. In the process of evolution, natural selection is carried out by simulating the evolution of biological populations to achieve the purpose of survival of the fittest. To achieve natural selection, a certain parameter is necessary to guide the evolution of the population. Therefore, we combine the indicator-based idea with NSGA-III and add the concept of indicator into NSGA-III to guide the natural selection of the population.

In order to propose the NSGA-III-COE algorithm, we introduce the set coverage (SC) [30]. Assuming that both set *A* and set *B* are obtained approximate solution sets, the numerator of formula (10) represents the number of solutions in which the solution in *B* is dominated by at least one solution in *A*, and the denominator represents the total number of solutions contained in *B*. The SC is the probability that the solutions in *B* is dominated by at least one solution in *A*.

$$\mathcal{C}(A, B) = \frac{|\{x \in B | \exists y \in A : y \text{ dominates } x\}|}{|B|},\tag{10}$$

Each subpopulation in the initial population has the same number of individuals. In the evolution process, the evolution of biological populations is simulated, and a small number of individuals are randomly exchanged in each iteration to increase the diversity of chromosomes in the decision space and to increase the amount of useful information in the decision space.

When the evolution reaches half of the maximum number of iterations, the SC indicator intervenes. The subpopulation size is adjusted every 10 generations according to the SC indicator. The SC indicator makes natural selection of the subpopulation based on the exploration and the development ability of the subpopulation in the decision space. Natural selection in the evolutionary process means increasing the size of superior subpopulations and reducing the size of disadvantaged subpopulations in order to achieve the survival of the fittest. Algorithm 1 describes the evolutionary mechanism of the NSGA-III-COE.

#### **Algorithm 1:** Evolutionary mechanism of the NSGA-III-COE.


#### **4. Experimental Results and Discussion**

In order to verify the advantages of the NSGA-III-COE in the MO-FJSP decision space exploration and development capabilities, a large number of computational experiments were carried out. These experiments were implemented by MATLAB programming and were tested on three sets of well-known benchmark instances, including 5 Kacem instances (ka1, ka2, ka3, ka4, ka5) [22], 10 BRdata instances (mk1–mk10) [21], and 12 BRdata instances (01a–12a) [22]. These collections cover most of the problem instances used in the FJSP literature. In fact, most of the existing research only considers a small part of them. In our experiment, 27 instances are used to comprehensively evaluate the algorithm we propose.

Table 4 lists the parameter settings of the algorithm, and we use uniform parameter values for the algorithm. For all instances, the maximum number of iterations is set to 300, and it is the same for all implemented algorithms. When the number of iterations reaches the set maximum number of iterations, the algorithm terminates to ensure fair comparison. For each instance, all algorithms run independently 30 times starting with the same initial population.

**Table 4.** Experimental parameter settings for the co-evolutionary NSGA-III (NSGA-III-COE) performance evaluation.


We are not sure whether all the nondominant solutions of the 27 benchmark instances collected in this paper enable *IGD* to more accurately reflect the overall performance of the algorithm on all instances. Therefore, this paper uses the hypervolume (*HV*) [30] indicator to evaluate the overall performance of the algorithms for all 27 instances.

Suppose *P* is the solution set obtained by the algorithm, and *q* = (*q*1, *q*2, ··· , *qm*) *<sup>T</sup>* is a reference point in the target space, which is dominated by all the target vectors in the solution set *P*. Then, the *HV* for reference point *q* refers to the volume of the target space dominated by solution set *P* and bounded by reference point *q*.

$$HV(P,q) = \text{volume}\left(\underset{p \in P}{\cup} [p\_{1\prime}q\_1] \times [p\_{2\prime}q\_2] \cdot \dots [p\_{m\prime}q\_m]\right) \tag{11}$$

Figure 6 illustrates the meaning of *HV* in a two-dimensional target space. The *HV* indicator calculation does not require a reference point set and can comprehensively reflect the convergence and the diversity of the solution set. The larger the *HV* is, the better the solution set obtained by the algorithm is and the better the overall performance of the algorithm will be.

**Figure 6.** Diagram of the hypervolume (*HV*).

### *4.1. Comparison of NSGA-III-COE and NSGA-III Variants*

In this section, we compare the NSGA-III-COE with NSGA-III-CX, NSGA-III-OBX, and NSGA-III-PBX to verify whether the population evolution mechanism proposed in this paper can integrate the effects of the three different crossover operators, enhance exploration and development capabilities in the decision space, and improve the performance of the algorithm. Table 5 shows the normalized average HVs obtained by running four algorithms 30 times independently on 27 benchmark instances. In addition, the features of the instance are also listed in the table. The first column indicates the name of the instance, and the second column indicates the size of the instance, where *Nn* indicates the number of processes, and *Nm* indicates the number of machines.

First, we analyze the three algorithms, NSGA-III-CX, NSGA-III-OBX, and NSGA-III-PBX. In 14 instances, the NSGA-III-OBX obtains the largest HV. In 11 instances, the NSGA-III-CX obtains the largest *HV*, and the NSGA-III-PBX obtains the largest *HV* in two instances. This phenomenon validates our conjecture when studying the developmental capabilities of NSGA-III variants. The same crossover operator shows different capabilities for developing better solutions when facing different MO-FJSPs.

Then, we add the NSGA-III-COE for analysis. Except for instances ka1 and mk2, the HVs of the NSGA-III-COE in the remaining 25 instances are all larger than those of the other NSGA-III variants. In instance ka1, the four algorithms have the same HVs. In instance mk2, NSGA-III-COE and NSGA-III-PBX both have the largest HVs. With the increase of instance complexity, the *HV* value of the NSGA-III-COE increases more and more obviously than other algorithms. This shows that the NSGA-III-COE performs best in all 27 instances, and as the complexity of the instance increases, the NSGA-III-COE shows better performance.


**Table 5.** Performance evaluation of NSGA-III-COE, NSGA-III-CX (cycle crossover), NSGA-III-OBX (order-based crossover), and NSGA-III-PBX (partially mapped crossover); the average HVs of 27 problem cases running independently 30 times. For each instance, the results which are better than the others are marked in bold (these have the largest *HV* value).

> In order to further verify whether the evolutionary mechanism proposed in this paper reaches our original design intention, Figure 7 shows the performances of the four algorithms of NSGA-III-COE, NSGA-III-CX, NSGA-III-OBX, and NSGA-III-PBX in instance mk2–mk10 at obtaining the non-dominated solution set in the same coordinate system obtained above. It can be seen from the figure that the non-dominated solution set obtained by the NSGA-III-COE is closer to the Pareto frontier than the other three algorithms and has good diversity. This indicates that the evolutionary mechanism proposed in this paper enhances search and development capabilities in the MO-FJSP decision space and speeds up the convergence speed while maintaining the diversity of the population, thus the algorithm obtains better performance.

**Figure 7.** The non-dominated solutions of NSGA-III-COE, NSGA-III-CX, NSGA-III-OBX, and NSGA-III-PBX on instances mk2–mk10.

#### *4.2. Comparison of the NSGA-III-COE with Widely Used MOEAs*

This section compares the NSGA-III-COE with the existing widely used MOEAs. NSGA-III, NSGA-II, NSGA-II/strengthened dominance relation(NSGA-II/SDR) [31], improving the strength pareto evolutionary algorithm (SPEA2) [32] and hypervolume estimation algorithm for multi-objective optimization (HypE) [33] are chosen as the comparison algorithms because they are all currently widely used MOEAs that can be directly applied to solve the MO-FJSP after simple modification, and some of them have been used many times in the previous FJSP literature. Ahmadi et al. [34] used the NSGA-II to solve the FJSP with random failures. Bandyopadhyay et al. [10] compared the calculation results with NSGA-II and SPEA2 when solving a parallel machine scheduling problem. Yuan Y et al. [9] also used the NSGA-II for comparison when solving the FJSP. NSGA-II is used very frequently in solving scheduling problems. Therefore, this paper adds the newer modified NSGA-II algorithm NSGA-II/SDR to the comparison algorithms.

Table 6 shows the average HVs for 30 independent runs on 27 instances of NSGA-III-COE, NSGA-III, NSGA-II, NSGA-II/SDR, SPEA2, and HypE. The table shows that the HVs of the NSGA-III-COE on the three instances of ka1, ka2, and mk1 are slightly larger than those of the other algorithms. On the remaining 24 instances, the HVs of the NSGA-III-COE are much larger than those of the other algorithms. As the complexity of the

benchmark instances increases, the advantages of the NSGA-III-COE become increasingly more obvious.

**Table 6.** Performance evaluation of the NSGA-III-COE and other comparison algorithms; the average HVs of 27 problem instances running 30 times independently. For each instance, the results which are significantly better than the others are marked in bold (these have much larger *HV* value than the other algorithms).


Figures 8–10 show the non-dominated solution set in the same coordinate system obtained by the NSGA-III-COE and the comparison algorithm used in this section on almost all benchmark instances. It can be seen from the figure that the performance of the original NSGA-III is relatively close to that of other comparison algorithms. However, the non-dominated solution sets obtained by all the comparison algorithms are far away from the Pareto front, which obviously falls into the local optimum. The non-dominated solution obtained by the NSGA-III-COE is far superior to other comparison algorithms. This shows that the population evolution mechanism proposed in this paper is not only conducive to the improvement of convergence speed but also improves the ability to jump out of the local optimum on the MO-FJSP and greatly improves the performance of NSGA-III in solving the MO-FJSP.

**Figure 8.** The non-dominated solutions of NSGA-III-COE, NSGA-III, NSGA-II, NSGA-II/SDR, SPEA2, and HypE on instances ka3, ka4, and ka5.

**Figure 9.** The non-dominated solutions of NSGA-III-COE, NSGA-III, NSGA-II, NSGA-II/SDR, SPEA2, and HypE on instances mk2–mk10.

**Figure 10.** The non-dominated solutions of NSGA-III-COE, NSGA-III, NSGA-II, NSGA-II/SDR, SPEA2, and HypE on instances 01a–12a.

We count the running time of six algorithms on 27 instances, and the results are shown in Figure 11. The running time of the NSGA-III-COE is slightly longer than that of the NSGA-III. Compared with these widely used MOEAs, the computational cost of the NSGA-III-COE is at a moderate level.

**Figure 11.** Average CPU time (in seconds) consumed by NSGA-III-COE and other comparison algorithms on 27 problem instances.

#### *4.3. Sensitivity Analysis of Population Size*

In this section, we associate the performance of the NSGA-III-COE and the other five MOEAs with the population size, which ranges from 30 to 500. We study the sensitivity of the NSGA-III-COE to population size on 27 FJSP benchmark instances.

We use *SR* (sum of ranks) to represent the performance of the algorithm:

$$SR = rank\_{\rm ka1} + \dots + rank\_{\rm ka5} + rank\_{\rm mk1} + \dots + rank\_{\rm mk10} + rank\_{01a} + \dots + rank\_{12\nu} \tag{12}$$

where *rank* represents the ranking of an algorithm among all algorithms for a benchmark instance. The ranking of the algorithm is based on the *HV* value of the non-dominated solution set. The larger the *HV* value is, the smaller the *SR* is, which means a higher ranking. For example, *rank*ka1 represents the ranking of an algorithm on instance ka1. *SR* represents the cumulative sum rankings of an algorithm on all instances.

From Figure 12, the following experimental observation results can be obtained.

**Figure 12.** The impact of population size on the performance of the NSGA-III-COE and other multiobjective evolutionary algorithms (MOEAs) on 27 benchmark instances. The figure shows the sum of rankings of 27 instances obtained by each algorithm when the population size is 30, 60, 100, 200, 300, 400, and 500.


Considering experimental results, running time, and sensitivity to population size, the NSGA-III-COE is a very competitive algorithm for solving low carbon FJSP.

#### **5. Conclusions**

In this paper, a low carbon MO-FJSP mathematical model is established to minimize total completion time, total carbon emissions, and total machine load. This mathematical model is very close to the real production environment and conforms to the concepts of green manufacturing and sustainable development. By understanding the existing literature on the MO-FJSP research, this paper introduces five crossover operators into the NSGA-III to produce different NSGA-III variants. Then, the ability of different NSGA-III variants to explore and develop better solutions in the decision space on some FJSP benchmark instances is studied. Through research and analysis of the experimental results, the indicator-based thought is introduced into NSGA-III, and a new co-evolutionary mechanism incorporated with multi-crossover operator and natural selection is proposed, which combines the capabilities of different crossover operators to make the algorithm obtain better performance. Subsequently, we introduce the new evolutionary mechanism into the NSGA-III and propose the NSGA-III-COE. Using the NSGA-III-COE to solve low carbon MO-FJSP, multiple experiments are done. The NSGA-III-COE achieves good results in solving the MO-FJSP.

In the experiment, we compare the NSGA-III-COE with five existing widely used MOEAs on 27 benchmark instances in the three\ test sets of the FJSP. Experimental results show that the NSGA-III-COE has a strong ability to optimize the low carbon MO-FJSP, and the computational cost of solving the problem is similar to that of widely used MOEAs. Compared with other widely used algorithms, the NSGA-III-COE algorithm has obvious advantages in convergence speed and the ability to jump out of local optimum. Especially, it shows better performance when dealing with complex problem cases. Since the solving time of FJSP increases exponentially with the increase of the problem scale, our research is very meaningful for solving the low carbon MO-FJSP.

The work done in this paper only shows that the proposed the NSGA-III-COE is effective for solving the MO-FJSP. In the future, we will further study the production scheduling problem, continue to research and improve the algorithm, and apply the algorithm to solve other multi-objective production scheduling problems.

**Author Contributions:** Conceptualization, Y.W. and H.K.; Software, Y.W. and X.S.; Validation, Y.W. and X.S.; Formal Analysis, Y.W. and Y.S.; Investigation, H.K. and Q.C.; Data Curation, Y.W. and X.S.; Writing—Original Draft Preparation, Y.W.; Writing—Review & Editing, Y.W. and D.W.; Visualization, Y.W. and Y.S.; Funding Acquisition, X.S., H.K., Y.S. and Q.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China, grant number 61663046, 1876166. This research was funded by Open Foundation of Key Laboratory of Software Engineering of Yunnan Province, grant number 2020SE308, 2020SE309. This research was funded by Science Research Foundation of Yunnan Education Committee of China, grant number 2020Y0004.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing is not applicable to this article.

**Conflicts of Interest:** The authors declare no conflict of interest.
