Research on Low-Carbon-Emission Scheduling of Workshop under Uncertainty

Abstract

Focusing on the problems of uncertainty and carbon emissions in the manufacturing process, this paper studies the low-carbon-emission scheduling optimization problem. Firstly, the variations in workpiece processing time and delivery date are selected as the uncertainty factors. A low-carbon-emission scheduling model for uncertain job shops is constructed with the optimization objectives of the time index, carbon emission index, and robustness index. Secondly, an improved third-generation non-dominated sorting genetic algorithm (NSGA-III) is proposed. Based on the original NSGA-III algorithm, this algorithm introduces the state transition algorithm to perform state transformation, neighborhood sampling, selection update, and alternate rotation operations on the parent population, generating new candidate solutions. Finally, the scheduling model and the improved algorithm are applied to a workshop example. Through case study computation and result analysis, the feasibility and effectiveness of the model and algorithm in addressing the low-carbon-emission job shop scheduling problem under uncertainty are further verified.

1. Introduction

Increasing carbon emissions play a crucial role in global warming, especially in industrial fields. Industries produce 35% of carbon emissions in China, which hurts the environment [1,2]. Carbon emissions are not eco-friendly, and the associated repercussions are global warming, severe meteorological deviations, and serious environmental pollution [3]. With the development of people’s ecological awareness, low carbon emissions and energy efficiency in the manufacturing industry have been receiving much attention [4]. The workshop scheduling optimization method effectively reduces carbon emissions and relieves environmental pressure [5]. Various scheduling strategies result in varying levels of carbon emissions, and reasonable scheduling methods are essential to complex flexible shop-floors. To increase enterprise profitability and control carbon emissions in the manufacturing sector, it is critical to implement a low-carbon-emission scheduling optimization method for workshops.

Manufacturing enterprises must carry out high-efficiency, low-carbon-emission, and low-cost manufacturing to enhance their competitiveness. The scheduling optimization model is constructed to optimize the manufacturing workshop. Chen and Hao [6] applied the non-dominated sorting genetic algorithm (NSGA) to design a non-compact flow shop scheduling plan. They successfully solved the multi-objective optimization problem considering process connection, which carries profound theoretical and practical significance for enterprises, e.g., improving the scheduling of non-compact flow shops, the production efficiency, and the response to market situations. Han [7] proposed the improved NSGA-II to minimize the makespan, total energy consumption, and total costs of the flexible job shop cell scheduling problem to achieve flexible manufacturing. Sang et al. [8] established the many-objective job shop intelligent scheduling model with complex constraints. They proposed an improved intelligent decision optimization algorithm named NSGA-III-APEV based on NSGA-III, which can ensure that the decision-maker obtains a more suitable production scheduling scheme.

Various uncertainties in the production process can disturb the production system and affect the scheduling plan, so many studies are focused on workshop robustness to improve resistance to disturbances. Chaari et al. [9] proposed a genetic algorithm for robust hybrid flow shop scheduling to evaluate the quality of the robustness when faced with uncertainty. Allaoui et al. [10] introduced a robustness framework for the stochastic hybrid flow shop problem under uncertainty to find a robust solution that can resist changes in input data. Shen et al. [11] constructed the multi-objective stochastic flexible job shop scheduling problem model, considering three objectives: makespan, maximal machine workload, and robustness to uncertainties. Considering various uncertainties in the actual manufacturing process, Zhang et al. [12] proposed a flexible job shop scheduling problem with machine breakdown, with the objectives of maximum completion time and robustness. Meanwhile, there is more and more research on low-carbon scheduling problems. Ding et al. [13] designed a multi-objective NEH algorithm and a modified multi-objective iterated greedy algorithm to solve a permutation flow shop model to minimize the carbon emissions and makespan. Pan and Lei [14] studied the problem of distributed low-carbon parallel machines. Dong and Ye [15] constructed a two-stage re-entrant hybrid-flow-shop bi-level scheduling model and used the NSGA-III algorithm to lower carbon emissions and energy costs. Tong and Zhu [16] proposed that green manufacturing has become a hot topic globally, and they designed a customer-oriented method to support multi-task green scheduling in cloud manufacturing, intending to minimize the total energy consumption during scheduling.

According to the reviewed literature, most low-carbon optimization scheduling studies focus on single machines or flow shops. However, the workshop is a complex, flexible job shop that is full of complexity, uncertainty, multiple constraints, and a multi-objective nature. Further research needs to be carried out on complex flexible workshops. Gen et al. [17] used the EDA method (understanding the estimation of distribution algorithm) to solve the stochastic flexible job shop scheduling problem with a min–max regret version and processing time uncertainty. To solve the flexible job shop scheduling problem, Kacem et al. [18] designed a hybrid Pareto algorithm based on fuzzy logic to minimize the maximum completion time and the total machine load. Zheng et al. [19] developed a neighborhood search algorithm based on a multi-objective group to solve the fuzzy flexible job shop scheduling problem by minimizing the maximum fuzzy completion time and machine load. Chang et al. [20] proposed an overall DT-enabled real-time scheduling (DTE-RS) framework for complex product shop-floors to effectively reduce adverse impacts of the dynamic disturbances and minimize the makespan based on a global twin. Piroozfard et al. [21] established a job shop scheduling model with total delay time and carbon emissions as objective functions and proposed a new multi-objective imperialist competition algorithm to solve the problem. Lin et al. [22] extended the flow shop scheduling problem to optimize low carbon emissions and variable processing parameters. They proposed two scheduling methods for specific problems to describe the scheduling problem. Lu and Jiang [23] proposed a bi-population-based discrete bat algorithm to research the low-carbon job shop scheduling problem.

To minimize the negative impact of uncertain factors on production, ensure the stability of the production system, and realize low-carbon and efficient operation of the manufacturing system, based on an in-depth discussion of the uncertain factors in the job shop, combined with a relevant literature search, this paper selects two typical uncertain factors (uncertain processing time and uncertain delivery date) as the uncertain factors in the production process. The corresponding mathematical tools are used to describe the uncertain factors, and the uncertainty-oriented low-carbon scheduling model of the job shop is constructed. Considering the lack of convergence in the genetic algorithm [24], this paper proposes an improved NSGA-III based on the state transition algorithm (NSGA-III-ST), which is proposed as an extension of the NSGA-III by introducing the state transition algorithm and enhancing the initial clustering operator, while maintaining the properties of the NSGA-III. The optimization of the manufacturing system’s low-carbon, efficient, and orderly operation is realized by applying the model and algorithm to the actual production scenario.

The main contributions of the paper are as follows:

• A multi-objective low-carbon-emission scheduling model for workshops under uncertainty is established, which optimizes green indexes such as carbon emissions and economic indexes such as the robustness index and time-based index, using the processing time and order delivery time as uncertainty. By optimizing the workshop model, the manufacturing process can achieve low-carbon, economic, and efficiency goals.
• An improved NSGA-III-ST algorithm is proposed. The algorithm is a hybrid of the NSGA-III and the state transition algorithm. The test also validates the feasibility and validity.

Our work has critical academic significance and engineering value for multi-objective optimization problems in the context of low carbon goals. The remainder of the paper is organized as follows: the second section constructs the low-carbon-emission scheduling model, the third section proposes the NSGA-III-ST algorithm, the fourth section solves the optimization model, the fifth section discusses the algorithm’s advantages and disadvantages, and the sixth section presents the conclusions. The method flow chart is shown in Figure 1.

2. Low-Carbon-Emission Scheduling Optimization Modeling

2.1. Description and Handling of Uncertainty

There are many uncertainties in the manufacturing process, and the emergence of uncertainties will strongly influence the accuracy of workshop scheduling [25]. A workshop’s uncertainties can be divided into four categories based on different sources: the system’s inherent uncertainty, uncertainty of equipment, uncertainty of operation personnel, and external factors caused by sale behavior or purchasing behavior. Because the system’s intrinsic uncertainty is generally not the focus of research, and the uncertainties caused by personnel or management methods are challenging to measure and calculate, uncertainty research focuses on the equipment, which consists mainly of machine breakdown, order insertion, order cancellation, workpiece processing time, unknown workpiece arrival time, order delivery time and raw material delay, and so on. After the weights of uncertain factors and the objective function are determined by the expert scoring method, and sensitivity analysis is performed by an analytic hierarchy process, the importance of the workshop uncertainty is determined as the workpiece processing time, order delivery time, machine breakdown, raw material delay, unknown workpiece arrival time, order insertion, and order cancellation. Meanwhile, the uncertainty of the processing time leads to uncertainty of the order delivery time, which affects the satisfaction degree (SD), and the uncertainty of processing time affects the makespan and carbon emissions. Therefore, the processing and order delivery times are selected as the uncertainties of the low-carbon-emission workshop scheduling problem.

The processing time cannot be precisely measured in manufacturing workshops; fuzzy sets are frequently used to describe it [26,27]. The triangular fuzzy number expresses the fuzzy processing time. Furthermore, because trapezoidal fuzzy numbers are commonly used when considering advance and delay, a trapezoidal fuzzy number expresses the fuzzy delivery time.

(a)

Fuzzy processing time

The triangular fuzzy number $\tilde{A}=(a_1,a_2,a_3)$ is used to denote the fuzzy processing time of the workpiece, where $a_1$ denotes the shortest processing time of the process, $a_2$ is the most probable processing time, and $a_3$ denotes the longest processing time. The probability distribution function of the affiliation function ${\mu }_{\tilde{A}}$ is illustrated in Equation (1). When $a_2=(a_1+a_3)/2$, the affiliation degree value of the triangular fuzzy number is 1.

$${\mu }_{\tilde{A}}=\left\{\begin{array}{c}\frac{x-a_1}{a_2-a_1},a_1\le x\le a_2\\ \frac{x-a_3}{a_2-a_3},a_2\le x\le a_3\\ 0,others\end{array}\right.$$

(1)

Set the triangular fuzzy number $\tilde{B}=(b_1,b_2,b_3)$. The addition of triangular fuzzy numbers can be expressed as $\tilde{A}+\tilde{B}=\left(a_1,a_2,a_3\right)+\left(b_1,b_2,b_3\right)=(a_1+b_1,a_2+b_2,a_3+b_3)$, where $b_1$ denotes the shortest processing time of the process, $b_2$ is the most probable processing time, and $b_3$ denotes the longest processing time. The symbol ∨ denotes taking the more significant number of fuzzy numbers, and the operation of taking the larger one is $\tilde{A}\vee \tilde{B}=\left(a_1,a_2,a_3\right)\vee \left(b_1,b_2,b_3\right)\simeq (a_1\vee b_1,a_2\vee b_2,a_3\vee b_3)$. The chosen large operation is chosen to approximate the fuzzy number, because the fuzzy number obtained by the more extensive operation is not always the triangular fuzzy number.

(b)

Fuzzy delivery time

The trapezoidal fuzzy number $\tilde{C}=(c_1,c_2,c_3,c_4)$ is used to express the fuzzy delivery period. $c_1$ denotes the shortest delivery period, $c_2$ and $c_3$ are probable delivery periods, and $c_4$ denotes the most extended delivery period. And ${\mu }_{\tilde{c}}$ is displayed in Equation (2). From the affiliation function, it is known that when $c_2=c_3$, the trapezoidal fuzzy number is transformed into a triangular fuzzy number; when $[c_1,c_2]=[c_3,c_4]$, the trapezoidal fuzzy number is symmetric about $l=(c_2+c_3)/2$, and the trapezoidal fuzzy number is called the symmetric trapezoidal fuzzy number.

$${\mu }_{\tilde{A}}=\left\{\begin{array}{c}\frac{x-c_1}{c_2-c_1},c_1\le x\le c_2\\ 1,c_2\le x\le c_3\\ \frac{c_4-x}{c_4-c_3},c_3\le x\le c_4\\ 0,others\end{array}\right.$$

(2)

2.2. Problem Description

The problem of low-carbon-emission workshop scheduling under uncertainty is generally described as follows: there are n mutually independent workpieces $J=\left\{J_i|i\in \mathrm{1,2},\dots ,n\right\}$ to be manufactured on $m$ machines $M=\left\{M_k|k\in \mathrm{1,2},\dots ,m\right\}$. Each workpiece $J_i$ is manufactured by one process or more processes. Each process can be constructed from a machine selected from the machine set $M$. There is a sequence constraint among each process, and only after the previous process is completed can the following process be processed. There are uncertainties in the process, such as the processing time and delivery time, so the processing time and fuzzy delivery time are given to characterize the range of these uncertainties.

The model of this study is based on the following assumptions:

• The production task starts at 0 time, and the processing route of the workpiece has been determined.
• The machine can only process one workpiece at a time, and there are no interruptions.
• The machine is operational at the outset. It begins with the first workpiece and does not stop until all workpieces have been completed.
• The time it takes to replace a machine tooling fixture and the time it takes to handle a workpiece between different machines is neglected.
• The processing time of each process is not determined.
• Each workpiece can only be manufactured on one machine at a time.
• There are no constraint level constraints between different workpieces and processes.

2.3. Parameter Description

The description of each parameter of the model is shown in

Table 1. Description of each parameter of the model.

Parameters	Parameter Meaning
n	Total number of processed workpieces
m	Total number of processed machines
l	Total number of processing processes
i	Workpiece subscript
k	Machine subscript
j	Operation subscript
$J_i$	$\mathrm{Workpiece} i,i\in \left\{\mathrm{1,2},\cdots ,n\right\}$
$M_k$	$\mathrm{Machine} \mathrm{k}, k\in \left\{\mathrm{1,2},\cdots ,m\right\}$
$O_{ij}$	jth process of the workpiece i
$\widetilde{D_i}$	$\mathrm{Vague} \mathrm{delivery} \mathrm{time} \mathrm{of} \mathrm{workpiece} J_i$
$\widetilde{C_i}$	$\mathrm{Fuzzy} \mathrm{completion} \mathrm{time} \mathrm{of} \mathrm{workpiece} J_i$
$\widetilde{C_{max}}$	Maximum fuzzy completion time
$\widetilde{T_{ijk}}$	$\mathrm{Fuzzy} \mathrm{processing} \mathrm{time} \mathrm{of} \mathrm{operation} O_{ij}$  $\mathrm{on} \mathrm{machine} M_k$
$\widetilde{t_{ijk}}$	$\mathrm{Fuzzy} \mathrm{idle} \mathrm{time} \mathrm{of} \mathrm{operation} O_{ij}$  $\mathrm{on} \mathrm{machine} M_k$
$\widetilde{C_{ijk}}$	$\mathrm{Fuzzy} \mathrm{completion} \mathrm{time} \mathrm{of} \mathrm{operation} O_{ij}$  $\mathrm{on} \mathrm{machine} M_k$  $, \mathrm{and} \widetilde{C_{ijk}}=\widetilde{T_{ijk}}+\widetilde{t_{ijk}}$
$\widetilde{{ST}_{ijk}}$	$\mathrm{Fuzzy} \mathrm{start} \mathrm{time} \mathrm{of} \mathrm{operation} O_{ij}$  $\mathrm{on} \mathrm{machine} M_k$
$\widetilde{{ET}_{ijk}}$	$\mathrm{Fuzzy} \mathrm{end} \mathrm{time} \mathrm{of} \mathrm{operation} O_{ij}$  $\mathrm{on} \mathrm{machine} M_k$
$F_e$	Carbon emission factor of electricity
$F_c$	Carbon emission factor of coolant
$P_{ik}$	$\mathrm{Machining} \mathrm{power} \mathrm{of} \mathrm{machining} \mathrm{on} \mathrm{the} \mathrm{machine} M_k$  $\mathrm{for} \mathrm{workpiece} J_i$
$T_f$	$\mathrm{Effective} \mathrm{cycle} \mathrm{period} \mathrm{of} \mathrm{coolant} \mathrm{on} \mathrm{machine} M_k$
$L_{fk}$	$\mathrm{Amount} \mathrm{of} \mathrm{coolant} \mathrm{circulating} \mathrm{on} \mathrm{the} \mathrm{machine} M_k$
$x_{ijk}$	$\mathrm{Boolean} \mathrm{value}: \mathrm{when} \mathrm{the} \mathrm{operation} O_{ij}$$\mathrm{processes} \mathrm{on} \mathrm{machine} M_k$, the value is 1; otherwise, the value is 0
$x_{ijkt}$	$\mathrm{Boolean} \mathrm{value}: \mathrm{when} \mathrm{the} \mathrm{operation} O_{ij}$$\mathrm{processes} \mathrm{on} \mathrm{machine} M_k$ at the time t, the value is 1; otherwise, the value is 0

.

2.4. Objective Function

Based on the problem description and problem assumptions, the scheduling goals are a time-based index ($TBI$), sum carbon emissions ($SCE$), and robustness index ($\Delta \delta$) depreciation, which are shown as $F_1$, $F_2$, and $F_3$.

The relationship between the fuzzy completion time and fuzzy delivery time is not only ahead of the time or delay, so the completion time satisfaction ($A{I}_i$) is introduced, which is shown as $A{I}_i=\frac{area(\widetilde{C_i}\cap \widetilde{D_i})}{area\widetilde{C_i}}$. The index value represents the ratio of the cross area of $C_i$, and $D_i$ membership functions to the area of the whole membership function of $C_i$. The Average Weighted Satisfaction ($AWS$) and the Production Minimum Satisfaction ($PMS$) are introduced to comprehensively evaluate $A{I}_i$, where $AWS=\frac{1}{n}{\sum }_{i=1}^n{w}_{i}A{I}_i$, $PMS=\mathrm{m}\mathrm{i}\mathrm{n}(A{I}_1,A{I}_2,\dots ,A{I}_n$), and $w_i$ represents the weighting coefficient of SD, which reflects the importance of a customer’s requirements for delivery time, and ${\sum }_{i=1}^n{w}_i=1$. $TBI={\theta }_{1}AWS+{\theta }_{2}PMS$, where ${\theta }_1$ and ${\theta }_2$ represent the weighting coefficients of $AWS$ and $PMS$, and ${\theta }_1+{\theta }_2=1$. Therefore, the objective function $F_1$ is shown in Equation (3).

$$F_1=\mathit{min}(TBI)=\mathit{min}({\theta }_{1}AWS+{\theta }_{2}PMS)$$

(3)

The model only considers the fuzzy carbon emissions generated by power loss during the machining process $\widetilde{{CE}_p}$, machine tool no-load process $\widetilde{{CE}_r}$, coolant consumption production $\widetilde{{CE}_f}$, and the public electricity consumption $\widetilde{PCE}$ generated by the operation of public auxiliary facilities. The flow of coolant is the same when different parts are machined on the same machine, and the power of public auxiliary facilities $P_{sumRP}$ includes the lighting facilities’ power $P_{sumLED}$, management kanban power $P_{sumBoard}$, and ventilator power $P_{sumFan}$.

$$\widetilde{C{E}_p}={\sum }_{i=1}^n{\sum }_{k=1}^m{\sum }_{j=1}^l{F}_e{x}_{ijk}{P}_{ik}\widetilde{T_{ijk}}$$

(4)

$$\widetilde{C{E}_r}={\sum }_{i=1}^n{\sum }_{k=1}^m{\sum }_{j=1}^l{F}_e{x}_{ijk}{P}_{ko}+\widetilde{t_{ijk}}$$

(5)

$$\widetilde{C{E}_f}={\sum }_{i=1}^n{\sum }_{k=1}^m{\sum }_{j=1}^l\frac{\widetilde{T_{ijk}}}{T_f}{F}_c{L}_{fk}{x}_{ijk}$$

(6)

$$P_{sumRP}=P_{sumLED}+P_{sumBoard}+P_{sumFan}$$

(7)

$$\widetilde{PCE}=F_e\widetilde{C_{max}}{P}_{sumRP}$$

(8)

Therefore, the objective function $F_2$ is shown in Equation (9).

$$F_2=\min \left(SCE\right)=\mathrm{m}\mathrm{i}\mathrm{n}(\widetilde{C{E}_p}+\widetilde{C{E}_r}+\widetilde{C{E}_f}+\widetilde{PCE})$$

(9)

The uncertainty of the fuzzy completion time can be determined by the difference in the span between the shortest and longest times, which reflects robustness. So the $△{\delta }_1$ and $△{\delta }_2$ are introduced, where $△{\delta }_1={\mu }_1\left(a_1+a_3-2a_2\right)+{\mu }_2{a}_2$, $△{\delta }_2=\sigma (a_3-a_1)$, ${\mu }_1$, ${\mu }_2$, and $\sigma$ are the non-negative weights, and ${\sum }_{i=1}^2{\mu }_i=1$. Therefore, $F_3$ is expressed in Equation (10).

$$F_3=\mathit{min}\left(△\delta \right)=\mathit{min}(△{\delta }_1+△{\delta }_2)$$

(10)

The objective function of the model is shown in Equation (11).

$$F=\mathrm{o}\mathrm{p}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{l}\left(F_1,F_2,F_3\right)=(\min \left(TBI\right),\min \left(SCE\right),\min \left(\mathsf{\Delta }\delta \right))$$

(11)

The scheduling aims to select the best machine and sequence for machining the workpiece to achieve the optimization goals. According to references [28,29], the constraints of the model are as follows:

$$\sum _{k=1}^m{x}_{ijk}=1,\forall i,j$$

(12)

$$\sum _{i=1}^n\sum _{j=1}^l{x}_{ijkt}\le 1,\forall k,t$$

(13)

$$\widetilde{S{T}_{ijk}}\ge \widetilde{E{T}_{i(j-1)k}},\forall i,j,k$$

(14)

$$\widetilde{E{T}_{ijk}}-\widetilde{S{T}_{ijk}}=\widetilde{T_{ijk}},\forall i,j,k$$

(15)

where Equation (12) indicates that a process can only be operated on one machine; Equation (13) represents that a machine can only process one process at most at any time; and Equation (14) indicates the process sequence constraints of the same workpiece, which is that only after the completion of the previous process can the following process begin to be processed. Equation (15) means that the machine only processes a workpiece at a time without interruptions.

3. Methods

NSGA-III is designed to improve the multi-objective optimization problem by changing the selection [30]. Despite its high diversity, NSGA-III has low precision and poor searchability when solving specific problems. Therefore, the improved NSGA-III-ST algorithm is proposed.

The NSGA-III-ST algorithm introduces the state transition algorithm and improves the original clustering operator. The state transition algorithm is used in the NSGA-III-ST to generate candidate solutions for the state transition of random initialization populations. The candidate solutions are further compared, filtered, and iteratively optimized by the Pareto sorting, population association, and microhabitat operations of NSGA-III to obtain the Pareto optimal solution set. The state transition algorithm is set up with various state operators, each performing a different function to guarantee the characteristics of the NSGA-III, and it can avoid the issue of traditional optimization algorithms’ prematureness or stagnation. The algorithm flow is shown in Figure 2.

3.1. Crossover and Mutation

Since the POX (Precedence Operation Crossover) operator can ensure the solution’s legality and improve the algorithm’s search ability, the POX is selected as the crossover operator [30]. Taking the coding crossover of 3 workpieces and 3 machines, the crossover is shown in Figure 3, with Parent1 and Parent2 as parents and Children1 and Children2 as children.

The workpieces are randomly assigned to the non-empty subsets $J_1$ and $J_2$, with $J_1=\left\{2\right\}$ and $J_2=\left\{\mathrm{1,3}\right\}$. We then copy the dark codes belonging to Parent1 to Children1 and the dark codes belonging to Parent2 to Children2 in $J_1$ and keep the position. We copy the light codes belonging to Parent2 to Children1 and light codes belonging to Parent1 to Children2 in $J_2$ and keep the position.

Insertional mutations are selected among the mutation ways, randomly selecting a gene to insert into a random position.

3.2. Encoding and Decoding

The coding method is composed of n workpieces, $J_i$, and the number of occurrences of each workpiece is equal to the process number $O_i$. When processing starts, the workpieces are arranged from left to right for processing.

Table 2. Process and fuzzy processing time in 3 × 3 processing scene.

Workpiece	Process	Machine
		${\mathit{M}}_{\mathbf{1}}$	${\mathit{M}}_{\mathbf{2}}$	${\mathit{M}}_{\mathbf{3}}$
$J_1$	1	(2, 4, 6)	(4, 5, 8)	(8, 9, 10)
	2	(4, 6, 8)	(3, 7, 10)	(1, 3, 5)
	3	(4, 5, 8)	(1, 4, 7)	(5, 6, 7)
$J_2$	1	(2, 3, 5)	(6, 8, 10)	(7, 8, 10)
	2	(3, 6, 10)	(4, 5, 7)	(2, 5, 8)
$J_3$	1	(3, 5, 8)	(2, 4, 5)	(7, 9, 11)
	2	(6, 8, 9)	(4, 8, 10)	(1, 2, 4)
	3	(2, 3, 4)	(5, 7, 8)	(8, 10, 13)

provides a more explicit description.

We set the process code as [1,1,3,2,3,1,2,3], as shown in Figure 4. The first “1” represents the first process of workpiece $J_1$, and the second “1” represents the second process of workpiece $J_1$. The corresponding processing machine should be selected according to the processing time. From

Table 3. Pseudocode for sampling strategy.

$\mathrm{Input}: SE, Bes{t}_k, \delta , R_a$

$\mathrm{Output}: State$

1. $\mathrm{for} i\leftarrow 1, SE \mathrm{do}$ 2. $State(:,i)\leftarrow Bes{t}_k+\delta R_{a}Bes{t}_k$ 3. end for

, the machining operation $O_{11}$ of $J_1$ is processed on machine $M_1$, and the machining operation $O_{12}$ of $J_1$ is processed on machine $M_3$. In turn, the processing process can be known as $O_{11}\to O_{12}\to O_{31}\to O_{21}\to O_{32}\to O_{13}\to O_{22}\to O_{33}$, and the corresponding machine is $M_1\to M_3\to M_2\to M_1\to M_3\to M_2\to M_3\to M_1$.

After the process-based encoding, the chromosomes containing process sequencing and machine selection information are decoded. The process should be allocated according to the earliest possible completion time. We suppose that the workpiece processing time is less than the idle time of the machine. In that case, the current process is inserted into the idle time of the machine, so that each process can be allocated to a better machine to complete the processing.

3.3. State Transition Operator

The state transition operator includes rotation transformation, translation transformation, expansion transformation, and axesion transformation [31,32]. Each transformation process takes place as follows: Rotation transformation (RT) is calculated as $x_{k+1}=x_k+a\frac{1}{n\left|\right|x_k|{|}_2}{R}_r{x}_k$, where $a$ is a positive constant, called the rotation factor, and $x_{k+1}$ denotes a candidate solution, while $x_k$ is the current state, $R_r\in R^{n\times n}$ is the random matrix with its elements belonging to the range of $[-1, 1]$, and $\left|\right|x|{|}_2$ represents the Euclidean norm. Translation transformation (TT) realizes a function of linear search for the line connected from point $x_{k-1}$ to point $x_k$, with the origin $x_k$ and maximum search length β, where β is a positive constant; $R_t\in R$ is a random variable with the elements belonging to the range of $\left[\mathrm{0,1}\right]$, whose formula is $x_{k+1}=x_k+\beta R_t\frac{x_k-x_{k-1}}{\left|\right|x_k-x_{k-1}|{|}_2}$. Expansion transformation (ET) is calculated as $x_{k+1}=x_k+\gamma R_e{x}_k$, where $\gamma$ is a positive constant called the expansion factor. $R_e\in R^{n\times n}$ is a random diagonal matrix, and the elements in the matrix conform to Gaussian distribution. Axesion transformation (AT) is calculated as $x_{k+1}=x_k+\delta R_a{x}_k$, where $\delta$ is a positive constant, called the axesion factor, and $R_a\in R^{n\times n}$ is a random diagonal matrix, which obeys Gaussian distribution, and only one random index has value.

3.4. Neighborhoods and Sampling

The candidate solution $x_{k+1}$ is generated by the current state $x_k$ through the state transformation operator. At the same time, since each state transformation operator contains a random matrix transformation, an infinite number of $x_{k+1}$ can be generated and form a neighborhood, which can be denoted as $N_{x_k}^{operator}$, where operator represents one of the four operators. Taking AT as an example, the pseudocode of the sampling strategy is illustrated in Table 3.

Here, $Bes{t}_k$ denotes the current optimal solution, $SE$ represents the sample size, and State means the candidate solution set with $SE$, generated by the current optimal solution $Bes{t}_k$ through the axial operator under a given sampling method.

3.5. Select and Update

The selection of individuals mainly consists of two aspects: one is to select the optimal solution from the candidate state’s set state containing $SE$ candidate solutions, and the other is to compare the chosen optimal solution with the original $Bes{t}_k$ to update the population, as given by Equation (16). The greedy criterion operates on the solution in the process, where $newBest$ denotes the selected optimal solution.

$$Bes{t}_k=\left\{\begin{array}{c}newBest, if \left(newBest\right)<f\left(Bes{t}_k\right)\\ Bes{t}_{k+1}, others\end{array}\right\}$$

(16)

3.6. Alternate Rotation

The pseudocode for the step is illustrated in

Table 4. Pseudocode for alternating rotations.

$\mathrm{Input}: Bes{t}_0, SE, func, a_{min}, a_{max}, \alpha , \beta , \gamma , \delta , \theta$

$\mathrm{Output}: Best$

1.   $Bes{t}_0\leftarrow Random Initialization$ 2.   repeat 3.   $\mathrm{if} a<{a \mathrm{then}}_{min}$ 4.   $\alpha \to {\alpha }_{max}$ 5.   end if 6.   $Bes{t}_k\leftarrow \mathrm{expansion} (func, Bes{t}_k, SE, \beta , \gamma )$ 7.   $Bes{t}_k\leftarrow \mathrm{rotation} (func, Bes{t}_k, SE, \beta , \alpha )$ 8.   $Bes{t}_k\leftarrow \mathrm{axesion} (func, Bes{t}_k, SE, \beta , \delta )$ 9.   $\alpha \leftarrow \theta$ 10. until terminal condition 11. $Best\leftarrow Bes{t}_k$

, where $Bes{t}_0$ is the initial solution generated randomly. $a_{min}$ and ${\alpha }_{max}$ denote the value ranges of the rotation factor, respectively, and $\theta$ is a constant.

3.7. Improvement of Clustering Operator

To decrease computational complexity, the improved clustering operator replaces the vertical distance between the reference line and the individuals with the angle between the reference line and the individuals. We set $f^n\left(x\right)=({f_1}^n\left(x\right),{f_2}^n\left(x\right),\dots ,{f_i}^n\left(x\right))$ as the normalized target and let ${\lambda }^j$ be the reference point and $\theta$ the acute angle between $f^n\left(x\right)$ and ${\lambda }^j$. The niche operation is performed after taking $\theta$ as the minimum individual’s set and the corresponding reference point’s set.

3.8. Performance Test of the Algorithm

To better reflect the comprehensive performance of the algorithm, the three indexes Hypervolume (HV), Inverted Generational Distance (IGD), and Averaged Hausdorff Distance (Δp) are used to evaluate the performance of the algorithm [33,34].

(a)

HV

The HV calculates the volume covered by the Pareto front in the objective space, a comprehensive index. The larger the index is, the better the diversity and convergence of the obtained Pareto front are [35]. The formula is given by Equation (17).

$$\mathrm{HV}=\delta \left(U_{i=1}^{\left|Q\right|}{v}_i\right)$$

(17)

where $\delta$ denotes the Leberger measure. $\left|Q\right|$ denotes the numbers of the Pareto front obtained by solving the algorithm. $v_i$ represents the Hypervolume formed by the reference point and the ith solution in the resulting solution set.

(b)

IGD

The IGD is an inverse mapping of the convergence metric generation distance, representing the Euclidean distance between the solutions on the true Pareto front and the calculated non-dominated solution set. In addition to reflecting the algorithm’s convergence [36], it can also reflect the uniformity and extensiveness of the solution set’s distribution. A smaller value indicates that the solved Pareto front is closer to the true one. The formula is given by Equation (18).

$$\mathrm{IGD}(p^{\ast },D)=\frac{{\sum }_{v\in p^{\ast }}dis(v,D)}{\left|p^{\ast }\right|}$$

(18)

where $p^{\ast }$ denotes the true Pareto optimal solution set; D denotes the non-dominated solution set that is actually solved by the algorithm; and $dis(v,D)$ means the minimum Euclidean distance from the individual in $p^{\ast }$ to the corresponding point in D. $\left|p^{\ast }\right|$ represents the size of the true Pareto optimal solution.

(c)

$∆p$

$∆p$ is used to evaluate the average Hausdorff distance between the solution set obtained by the algorithm solution and the true Pareto front of a given multi-objective optimization problem. It is a comprehensive performance evaluation indicator combining improved GD and IGD indexes. The smaller the indicator is, the better the comprehensive performance of the algorithm is. The improved GD and IGD are denoted as GD_p and IGD_p, which are calculated in Equations (19) and (20), and Δp is calculated in Equation (21).

$$\mathrm{G}{\mathrm{D}}_p(S,P,q)=\left(\frac{1}{\left|S\right|}\sum _{x\in S}\underset{p\in P}{min}{d}^q\right(x,p){)}^{\frac{1}{q}}$$

(19)

$$\mathrm{IG}{\mathrm{D}}_p(S,P,q)=\left(\frac{1}{\left|P\right|}\sum _{p\in P}\underset{x\in S}{min}{d}^q\right(p,x){)}^{\frac{1}{q}}$$

(20)

$$∆p\left(S,P,q\right)=\mathrm{m}\mathrm{a}\mathrm{x}(G{D}_p\left(S,P,q\right),I{GD}_p\left(S,P,q\right))$$

(21)

where $x$ is the solution; $S$ is the solution set; $P$ is the Pareto approximation front; the size of $q$ affects the size of $∆p$, and the penalty of $∆p$ on discrete points is greater when $q$ is larger; and $d$ is the Euclidean distance.

3.9. Analysis of Test Results

To verify the superiority of the improved algorithm, MOEA/D, NSGA-II, NSGA-III, and the NSGA-III-ST algorithm are selected to perform simulation experiments on ZDT (ZDT1~ZDT4) and DTLZ (DTLZ1, DTLZ2, and DTLZ7) [37] test functions, which are common function sets. The four algorithms use simulated binary crossover and polynomial variation to achieve iterative population evolution, and the variation probability is $1/k$, where $k$ denotes the dimensionality of the corresponding decision variables. According to multiple tests, the rotation factor $\alpha$, maximum search length $\beta ,$ and scaling factor $\gamma$ are all set as 1, and the sample size is set as $SE=5$. Each algorithm runs 30 times independently, and the average value (Average) and standard deviation (St.dev) of the three indexes IGD, HV, and $∆p$ are taken to reflect the comprehensive performance of the algorithm. The optimal values of all indexes are marked in dark gray, and the suboptimal values are marked in light gray.

(a)

HV

It can be seen from

Table 5. The test results based on the HV.

Function	Algorithm	Average	St.dev
ZDT1	MOEA/D	7.1146 × 10−1	6.57 × 10−3
	NSGA-II	7.1913 × 10−1	2.34 × 10−4
	NSGA-III	7.2000 × 10−1	9.33 × 10−5
	NSGA-III-ST	7.2028 × 10−1	2.94 × 10−5
ZDT2	MOEA/D	4.1891 × 10−1	3.06 × 10−2
	NSGA-II	4.4379 × 10−1	1.89 × 10−4
	NSGA-III	4.4460 × 10−1	1.40 × 10−4
	NSGA-III-ST	4.4501 × 10−1	2.77 × 10−5
ZDT3	MOEA/D	5.9748 × 10−1	3.96 × 10−2
	NSGA-II	5.8232 × 10−1	2.44 × 10−2
	NSGA-III	5.8212 × 10−1	2.02 × 10−4
	NSGA-III-ST	5.8924 × 10−1	2.45 × 10−4
ZDT4	MOEA/D	6.9605 × 10−1	1.21 × 10−2
	NSGA-II	7.1726 × 10−1	1.70 × 10−3
	NSGA-III	7.1307 × 10−1	8.60 × 10−3
	NSGA-III-ST	7.1975 × 10−1	4.28 × 10−4
DTLZ1	MOEA/D	8.3994 × 10−1	1.53 × 10−3
	NSGA-II	8.2031 × 10−1	5.55 × 10−3
	NSGA-III	8.3962 × 10−1	1.65 × 10−3
	NSGA-III-ST	8.3884 × 10−1	9.24 × 10−3
DTLZ2	MOEA/D	5.5947 × 10−1	2.65 × 10−5
	NSGA-II	5.3180 × 10−1	4.47 × 10−3
	NSGA-III	5.5946 × 10−1	4.55 × 10−5
	NSGA-III-ST	5.5824 × 10−1	1.23 × 10−3
DTLZ4	MOEA/D	4.0454 × 10−1	1.35 × 10−1
	NSGA-II	5.0529 × 10−1	1.13 × 10−1
	NSGA-III	4.6308 × 10−1	1.28 × 10−1
	NSGA-III-ST	5.5790 × 10−1	1.48 × 10−3

that among the 7 test functions based on the HV, NSGA-III-ST obtains 4 optimal values and 1 suboptimal value, and the st.devs are smaller than other algorithms on some functions that achieve optimal values, such as ZDT1 and ZDT2; MOEA/D obtains 3 optimal values; NSGA-II obtains 2 suboptimal values; and NSGA-III obtains 5 suboptimal values. That means that NSGA-III-ST has a better comprehensive performance.

(b)

IGD

Similarly, the IGD is also tested, and the following results are obtained, as shown in

Table 6. The test results based on the IGD.

Function	Algorithm	Average	St.dev
ZDT1	MOEA/D	1.0569 × 10−2	8.96 × 10−3
	NSGA-II	4.7751 × 10−3	1.83 × 10−4
	NSGA-III	3.9120 × 10−3	1.20 × 10−5
	NSGA-III-ST	3.8896 × 10−3	3.89 × 10−6
ZDT2	MOEA/D	1.9229 × 10−2	3.01 × 10−2
	NSGA-II	4.9055 × 10−3	1.96 × 10−4
	NSGA-III	3.8716 × 10−3	4.50 × 10−5
	NSGA-III-ST	3.8098 × 10−3	2.69 × 10−6
ZDT3	MOEA/D	2.6567 × 10−2	1.30 × 10−2
	NSGA-II	7.3206 × 10−3	7.39 × 10−3
	NSGA-III	6.1249 × 10−3	2.05 × 10−4
	NSGA-III-ST	6.0997 × 10−3	1.82 × 10−4
ZDT4	MOEA/D	2.0121 × 10−2	1.08 × 10−2
	NSGA-II	5.3536 × 10−3	8.16 × 10−4
	NSGA-III	9.5565 × 10−3	1.29 × 10−2
	NSGA-III-ST	4.0198 × 10−3	1.67 × 10−4
DTLZ1	MOEA/D	2.0741 × 10−2	2.17 × 10−4
	NSGA-II	2.7259 × 10−2	1.51 × 10−3
	NSGA-III	2.0780 × 10−2	2.20 × 10−4
	NSGA-III-ST	2.1333 × 10−2	3.36 × 10−3
DTLZ2	MOEA/D	5.4467 × 10−2	1.33 × 10−6
	NSGA-II	6.9125 × 10−2	2.38 × 10−3
	NSGA-III	5.4480 × 10−2	7.48 × 10−6
	NSGA-III-ST	5.4668 × 10−2	6.80 × 10−4
DTLZ4	MOEA/D	3.8997 × 10−1	2.78 × 10−1
	NSGA-II	1.2666 × 10−1	2.23 × 10−1
	NSGA-III	2.6281 × 10−1	2.69 × 10−1
	NSGA-III-ST	5.5056 × 10−2	1.14 × 10−3

: NSGA-III-ST obtains 5 optimal values; MOEA/D obtains 2 optimal values; NSGA-II obtains 2 suboptimal values; and NSGA-III obtains 5 suboptimal values. It can be seen that NSGA-III-ST has a higher convergence than the other algorithms.

(c)

$∆p$

Meanwhile, $∆p$ is tested, and the following results are obtained, as shown in

Table 7. The test results based on the Δp.

Function	Algorithm	Average	St.dev
ZDT1	MOEA/D	1.0569 × 10−2	8.96 × 10−3
	NSGA-II	4.7751 × 10−3	1.83 × 10−4
	NSGA-III	3.9120 × 10−3	1.20 × 10−5
	NSGA-III-ST	3.8896 × 10−3	3.89 × 10−6
ZDT2	MOEA/D	1.9229 × 10−2	3.01 × 10−2
	NSGA-II	4.9055 × 10−3	1.96 × 10−4
	NSGA-III	3.8716 × 10−3	4.50 × 10−5
	NSGA-III-ST	3.8098 × 10−3	2.69 × 10−6
ZDT3	MOEA/D	2.6567 × 10−2	1.30 × 10−2
	NSGA-II	7.3206 × 10−3	7.39 × 10−3
	NSGA-III	6.1249 × 10−3	2.05 × 10−4
	NSGA-III-ST	6.0997 × 10−3	1.82 × 10−4
ZDT4	MOEA/D	2.0127 × 10−2	1.08 × 10−2
	NSGA-II	5.3536 × 10−3	8.16 × 10−4
	NSGA-III	9.5565 × 10−3	1.29 × 10−2
	NSGA-III-ST	4.0198 × 10−3	1.67 × 10−4
DTLZ1	MOEA/D	2.0725 × 10−2	1.56 × 10−4
	NSGA-II	2.7259 × 10−2	1.51 × 10−3
	NSGA-III	2.0780 × 10−2	2.20 × 10−4
	NSGA-III-ST	2.1333 × 10−2	3.36 × 10−3
DTLZ2	MOEA/D	5.4467 × 10−2	1.34 × 10−6
	NSGA-II	6.9125 × 10−2	2.38 × 10−3
	NSGA-III	5.4480 × 10−2	7.48 × 10−6
	NSGA-III-ST	5.4668 × 10−2	6.80 × 10−4
DTLZ4	MOEA/D	3.2778 × 10−1	2.70 × 10−1
	NSGA-II	1.2666 × 10−1	2.23 × 10−1
	NSGA-III	2.6281 × 10−1	2.69 × 10−1
	NSGA-III-ST	5.5056 × 10−2	1.14 × 10−3

: NSGA-III-ST achieves 5 optimal values; MOEA/D makes 3 optimal values; NSGA-II achieves 2 suboptimal values; NSGA-III achieves 5 suboptimal values. That means that NSGA-III-ST has better comprehensive performances.

The convergence of the four algorithms on the HV on the ZDT1 problem is drawn, as shown in Figure 5. It can be seen that the convergence speed of NSGA-III-ST has been greatly improved.

Although NSGA-III-ST has no optimal solution among the four algorithms in the DTLZ2 problem, where the Pareto optimal solution set is evenly distributed, the improved algorithm retains the evaluation method based on the reference vector in NSGA-III. It ensures the uniformity of the intersection point with the Pareto front. The HV and IGD indexes calculated by NSGA-III-ST are only slightly different from the optimal values, as shown in Figure 6.

In summary, NSGA-III-ST is superior to the other three algorithms in terms of comprehensive performances for HV, IGD, and $∆p$. Considering the convergence level and diversity, the NSGA-III-ST algorithm performs better for the Pareto front.

4. Results

4.1. Related Data on Workshop Scheduling Model

There are eight machines in a production workshop, and the set $\left\{M_i|i\in \mathrm{1,2},\dots ,8\right\}$ is used to denote the available processing machines set. The data on machine power and consumption are shown in

Table 8. Machine power and coolant consumption.

Machine	Processing Power (kW)	Idle Power (kW)	Coolant Usage (L)	Coolant Circulation Cycle (104 s)
$M_1$	11.5	2.45	600	120
$M_2$	12.5	1.82	600	120
$M_3$	11.5	1.50	400	90
$M_4$	12	1.58	400	90
$M_5$	10	1.41	400	90
$M_6$	6.5	0.45	350	86
$M_7$	7.5	0.71	350	86
$M_8$	10	1.80	350	86

.

The fuzzy processing time is obtained using previous data and staff experience, as shown in

Table 9. Fuzzy processing time on the machine for each product process (seconds).

Workpiece	Process	Machine
		${\mathit{M}}_{\mathbf{1}}$	${\mathit{M}}_{\mathbf{2}}$	${\mathit{M}}_{\mathbf{3}}$	${\mathit{M}}_{\mathbf{4}}$
$J_1$	1	(85, 90, 92)	(88, 95, 98)	(0, 0, 0)	(79, 82, 85)
	2	(0, 0, 0)	(35, 40, 45)	(32, 35, 38)	(0, 0, 0)
	3	(0, 0, 0)	(12, 15, 20)	(0, 0, 0)	(14, 16, 17)
$J_2$	1	(0, 0, 0)	(0, 0, 0)	(83, 85, 86)	(80, 81, 83)
	2	(50, 51, 53)	(45, 49, 51)	(0, 0, 0)	(0, 0, 0)
	3	(29, 31, 33)	(0, 0, 0)	(0, 0, 0)	(32, 35, 36)
	4	(0, 0, 0)	(0, 0, 0)	(19, 20, 22)	(0, 0, 0)
	5	(18, 19, 21)	(17, 19, 22)	(0, 0, 0)	(18, 22, 24)
	6	(0, 0, 0)	(10, 11, 12)	(0, 0, 0)	(8, 9, 11)
$J_3$	1	(81, 83, 85)	(0, 0, 0)	(79, 81, 83)	(0, 0, 0)
	2	(49, 51, 53)	(47, 49, 52)	(51, 53, 55)	(46, 48, 50)
	3	(0, 0, 0)	(32, 34, 36)	(0, 0, 0)	(33, 35, 36)
	4	(0, 0, 0)	(0, 0, 0)	(10, 11, 12)	(0, 0, 0)
$J_4$	1	(75, 76, 78)	(72, 73, 74)	(0, 0, 0)	(0, 0, 0)
	2	(0, 0, 0)	(38, 40, 42)	(38, 39, 42)	(36, 38, 39)
	3	(10, 12, 13)	(0, 0, 0)	(9, 10, 12)	(9, 10, 11)
$J_5$	1	(0, 0, 0)	(52, 55, 56)	(0, 0, 0)	(51, 53, 55)
	2	(0, 0, 0)	(0, 0, 0)	(0, 0, 0)	(0, 0, 0)
$J_6$	1	(81, 82, 83)	(0, 0, 0)	(79, 81, 82)	(0, 0, 0)
	2	(45, 46, 48)	(0, 0, 0)	(0, 0, 0)	(0, 0, 0)
	3	(10, 11, 12)	(0, 0, 0)	(9, 12, 13)	(8, 10, 12)
$J_7$	1	(0, 0, 0)	(71, 72, 73)	(70, 73, 76)	(0, 0, 0)
	2	(42, 44, 45)	(0, 0, 0)	(0, 0, 0)	(40, 42, 43)
	3	(0, 0, 0)	(23, 25, 26)	(22, 23, 24)	(0, 0, 0)
	4	(0, 0, 0)-	(0, 0, 0)	(8, 9, 11)	(0, 0, 0)
$J_8$	1	(81, 83, 85)	(79, 81, 83)	(0, 0, 0)	(0, 0, 0)
	2	(45, 46, 48)	(43, 45, 46)	(0, 0, 0)	(44, 45, 47)
	3	(21, 22, 24)	(0, 0, 0)	(19, 21, 22)	(0, 0, 0)
Workpiece	Process	Machine
		${\mathit{M}}_{\mathbf{5}}$	${\mathit{M}}_{\mathbf{6}}$	${\mathit{M}}_{\mathbf{7}}$	${\mathit{M}}_{\mathbf{8}}$
$J_1$	1	(0, 0, 0)	(92, 95, 98)	(84, 91, 96)	(0, 0, 0)
	2	(30, 33, 35)	(28, 30, 32)	(0, 0, 0)	(30, 31, 33)
	3	(0, 0, 0)	(0, 0, 0)	(12, 13, 15)	(0, 0, 0)
$J_2$	1	(85, 88, 90)	(89, 92, 95)	(0, 0, 0)	(79, 82, 85)
	2	(46, 48, 50)	(42, 45, 48)	(0, 0, 0)	(0, 0, 0)
	3	(0, 0, 0)	(30, 33, 35)	(31, 33, 35)	(0, 0, 0)
	4	(21, 22, 25)	(0, 0, 0)	(20, 25, 27)	(0, 0, 0)
	5	(0, 0, 0)	(19, 20, 22)	(20, 21, 22)	(0, 0, 0)
	6	(8, 10, 12)	(7, 8, 9)	(0, 0, 0)	(0, 0, 0)
$J_3$	1	(0, 0, 0)	(0, 0, 0)	(82, 83, 85)	(0, 0, 0)
	2	(0, 0, 0)	(0, 0, 0)	(0, 0, 0)	(0, 0, 0)
	3	(0, 0, 0)	(30, 35, 36)	(0, 0, 0)	(35, 36, 38)
	4	(9, 10, 12)	(8, 10, 12)	(0, 0, 0)	(9, 12, 13)
$J_4$	1	(0, 0, 0)	(70, 74, 78)	(71, 73, 75)	(0, 0, 0)
	2	(0, 0, 0)	(0, 0, 0)	(40, 41, 43)	(0, 0, 0)
	3	(0, 0, 0)	(8, 10, 12)	(0, 0, 0)	(11, 12, 13)
$J_5$	1	(49, 52, 54)	(0, 0, 0)	(0, 0, 0)	(0, 0, 0)
	2	(0, 0, 0)	(0, 0, 0)	(21, 22, 24)	(19, 21, 23)
$J_6$	1	(83, 85, 86)	(91, 93, 94)	(0, 0, 0)	(0, 0, 0)
	2	(44, 45, 46)	(48, 50, 51)	(52, 53, 55)	(0, 0, 0)
	3	(9, 11, 12)	(0, 0, 0)	(0, 0, 0)	(10, 12, 13)
$J_7$	1	(69, 71, 72)	(73, 75, 78)	(0, 0, 0)	(0, 0, 0)
	2	(0, 0, 0)	(39, 41, 42)	(43, 44, 46)	(38, 40, 42)
	3	(0, 0, 0)	(19, 20, 22)	(0, 0, 0)	(25, 26, 27)
	4	(8, 10, 12)	(0, 0, 0)	(10, 11, 12)	(0, 0, 0)
$J_8$	1	(0, 0, 0)	(0, 0, 0)	(0, 0, 0)	(82, 83, 85)
	2	(0, 0, 0)	(0, 0, 0)	(0, 0, 0)	(0, 0, 0)
	3	(20, 21, 23)	(0, 0, 0)	(18, 19, 21)	(0, 0, 0)

, where the value (0,0,0) indicates that the workpiece’s current process is not processed on that machine. Meanwhile, each workpiece has an internal delivery period, as shown in

Table 10. Fuzzy internal lead time for each product (seconds).

Workpiece	Delivery Time
$J_1$	(170, 180, 200, 210)
$J_2$	(210, 215, 231, 245)
$J_3$	(180, 198, 210, 220)
$J_4$	(150, 160, 175, 180)
$J_5$	(100, 125, 130, 145)
$J_6$	(160, 170, 180, 200)
$J_7$	(180, 190, 205, 210)
$J_8$	(150, 164, 182, 206)

. The delivery period is rounded to account for subsequent calculations.

Each machine in the workshop is equipped with a Management Kanban; lighting is provided by energy-saving fluorescent lamps of the T8 type, and axial-flow roof fans provide exhaust. The number and rated power of each utility involved in carrying out the batch production tasks are shown in

Table 11. Number and rated power of each public facility.

Name	Number	Rated Power (kW)	Total Power (kW)
Management Kanban	8	0.1	0.8
Energy-saving fluorescent lamp of T8 type	20	0.03	0.6
Axial flow fans	5	1.1	5.5

.

To compare with the optimized scheme, the scheduling scheme is intercepted for the same quantity at that time. In the original scheduling scheme, the carbon emission is 150.64 kg·CO₂, and the makespan is 260 s.

4.2. Solving of Model Based on NSGA-III-ST

The model is solved using the NSGA-III and NSGA-III-ST algorithms. We set the initial population size as 100, the variation probability as 0.8, and the number of iterations as 100. The simulation software used in this paper is Pycharm 2022.1.2, and the simulation experiment environment is AMD Ryzen 5 5600H with Radeon Graphics, 3.30 GHz. The RAM is 16G, with the 64-bit Windows 11-based operating system for x64 processors. The two algorithms run 10 times independently, and the results are illustrated in

Table 12. Example results of solving for the two algorithms.

Type	Parameter	NSGA-III			NSGA-III-ST
		TBI	$\mathit{S}\mathit{C}\mathit{E}$	$\mathit{\Delta }\mathit{\delta }$	TBI	$\mathit{S}\mathit{C}\mathit{E}$	$\mathit{\Delta }\mathit{\delta }$
Value	Best	0.08	136.13	197.32	0.03	133.51	195.71
	Average	0.90	142.71	216.48	1.78	135.81	206.32
Optimal	TBI minimization	0.08	145.43	236.35	0.03	143.39	231.31
	SCE minimization	2.52	136.13	219.83	0.30	133.51	228.76
	∆δ minimization	2.36	143.08	197.32	1.99	136.69	195.71

.

Table 12 shows that the best value (Best) obtained by the NSGA-III-ST algorithm is superior to those obtained by the NSGA-III algorithm, and the average of the $SCE$ index obtained by NSGA-III-ST is better than those of NSGA-III. Meanwhile, we run the NSGA-III-ST against the NSGA-III using the test functions DTLZ1-4, finding that the difference in running time between the two algorithms for the test functions is insignificant.

Figure 7 plots the Pareto solution sets obtained by the two algorithms. The pentagram represents the solution set obtained by the NSGA-III-ST, and the diamond denotes that of NSGA-III. Most of the diamond points are dominated by the pentagram points, which means that the ability of NSGA-III-ST to obtain Pareto solutions is outstanding. Considering that the NSGA-III-ST algorithm outperforms the NSGA-III in terms of solution performance, the NSGA-III-ST algorithm has a higher efficiency.

Three optimization results are taken out owing to the $TBI$ minimization, $SCE$ minimization, and $∆\delta$ minimization, as shown in Table 12. The fuzzy Gantt chart of the scheduling scheme is drawn based on the minimization of the three indexes, as shown in Figure 8, Figure 9 and Figure 10.

From Table 12, it can be seen that NSGA-III-ST makes better solutions in the Best parameter but achieves a larger average in the time index $TBI$, which means that the optimization scheme still has some room for improvement. In the optimization scheme based on the minimization of different indexes, the job is mainly concentrated in $M_4$, $M_5$, $M_6$, and $M_7$. The robust scheduling scheme generated by the NSGA-III-ST algorithm is feasible under the influence of uncertainties in terms of processing time and delivery period. The minimum scheduling scheme based on the NSGA-III-ST algorithm calculates the carbon emissions to be 133.51 kg CO₂ and the makespan to be 226 s, which means tha the carbon emissions are reduced, and the efficiency of the scheduling scheme is improved compared with the original scheme.

Therefore, the following improvements can be made: for equipment with relatively concentrated operations, replacing it with lower-energy-consuming equipment or considering adjusting processes; while meeting internal delivery requirements, the external delivery period should also be ensured, that is, meeting the customer’s needs in terms of both time and quality to improve customer satisfaction.

5. Discussion

Focusing on the problem of insufficient convergence of the NSGA-III algorithm, this paper introduces the state transition algorithm into the NSGA-III algorithm and proposes the NSGA-III-ST algorithm. The state transformation operator includes rotation, translation, expansion, and axesion transformations to improve the algorithm’s searchability. Neighborhood and sampling operations avoid traversing the elements in the neighborhood to improve the search speed. The selection and update operation selects the optimal solution from the state set of candidate solutions. The alternate rotation operation is used to realize the iterative optimization of the global optimal solution. The solution set generated by the above state changes is combined with the parent population. Then, the non-dominated sorting, population adaptive normalization, and niche preservation strategies are carried out. When the algorithm meets the iteration conditions, the optimal solution is output. The NSGA-III_ST algorithm is compared with the classical algorithm, and ZDT, UF, and DTLZ series benchmark functions and scheduling cases are selected for testing. The test results show the effectiveness of the proposed algorithm. However, the paper still has some limitations. Firstly, there are more uncertainties in the actual workshops, and the anti-interference capability of the model needs to be enhanced. The algorithm mainly focuses on updating and selecting candidate solutions, with insufficient research being carried out on reference point associations and microhabitat operations. For future research, equipment failure and the quantification of personnel operations, etc., should be considered, and in-depth research on reference point associations and microhabitat operations is suggested.

6. Conclusions

In this paper, to solve workshop scheduling problems under uncertainty, a low-carbon-emission workshop scheduling model under uncertainty and an improved NSGA-III-ST algorithm are proposed. The proposed model considers the scheduling objectives of a minimum time index, minimum total carbon emissions, and optimal robust indexes. The NSGA-III-ST algorithm is proposed by introducing the state transition algorithm and improving the original clustering operator based on NSGA-III, increasing the ability of the local search while maintaining the advantages of the NSGA-III algorithm. The model is applied to a real case. In the minimum scheduling scheme of the carbon emission index generated based on the NSGA-III-ST algorithm, the calculated carbon emission is 133.51 kg-CO₂, and the maximum completion time is 226 s. Compared with other scheduling schemes, the method has some advantages in reducing carbon emission and reducing the carbon emissions of automobiles by 11.3%, which improves, to a certain extent, the efficiency of the original scheduling scheme. The maximum completion time is reduced by 13.1%, which verifies the effectiveness of the model and algorithm.