*Article* **A Multi-Objective Optimization Method for Flexible Job Shop Scheduling Considering Cutting-Tool Degradation with Energy-Saving Measures**

**Ying Tian 1,2, Zhanxu Gao 1,2,\*, Lei Zhang 3,4,\*, Yujing Chen 1,2 and Taiyong Wang 1,2**


**Abstract:** Traditional energy-saving optimization of shop scheduling often separates the coupling relationship between a single machine and the shop system, which not only limits the potential of energy-saving but also leads to a large deviation between the optimized result and the actual application. In practice, cutting-tool degradation during operation is inevitable, which will not only lead to the increase in actual machining power but also the resulting tool change operation will disrupt the rhythm of production scheduling. Therefore, to make the energy consumption calculation in scheduling optimization more consistent with the actual machining conditions and reduce the impact of tool degradation on the manufacturing shop, this paper constructs an integrated optimization model including a flexible job shop scheduling problem (FJSP), machining power prediction, tool life prediction and energy-saving strategy. First, an exponential function is formulated using actual cutting experiment data under certain machining conditions to express cutting-tool degradation. Utilizing this function, a reasonable cutting-tool change schedule is obtained. A hybrid energysaving strategy that combines a cutting-tool change with machine tool turn-on/off schedules to reduce the difference between the simulated and actual machining power while optimizing the energy savings is then proposed. Second, a multi-objective optimization model was established to reduce the makespan, total machine tool load, number of times machine tools are turned on/off and cutting tools are changed, and the total energy consumption of the workshop and the fast and elitist multi-objective genetic algorithm (NSGA-II) is used to solve the model. Finally, combined with the workshop production cost evaluation indicator, a practical FJSP example is presented to demonstrate the proposed optimization model. The prediction accuracy of the machining power is more than 93%. The hybrid energy-saving strategy can further reduce the energy consumption of the workshop by 4.44% and the production cost by 2.44% on the basis of saving 93.5% of non-processing energy consumption by the machine on/off energy-saving strategy.

**Keywords:** cutting-tool degradation; machine tool turning-on/off schedule; hybrid energy-saving strategy; multi-objective optimization; flexible job shop scheduling

**MSC:** 90B30; 90B35

#### **1. Introduction**

In the current industrial environment, the manufacturing industry, as an important part, consumes a lot of energy and resources in the process of product manufacturing [1]. The report on power consumption released by the China Electricity Council in 2022 showed that China's industrial electricity consumption accounted for 64.5% of the total social electricity consumption in the first 10 months of 2022, while manufacturing electricity

**Citation:** Tian, Y.; Gao, Z.; Zhang, L.; Chen, Y.; Wang, T. A Multi-Objective Optimization Method for Flexible Job Shop Scheduling Considering Cutting-Tool Degradation with Energy-Saving Measures. *Mathematics* **2023**, *11*, 324. https:// doi.org/10.3390/math11020324

Academic Editor: Javier Alcaraz

Received: 12 December 2022 Revised: 4 January 2023 Accepted: 6 January 2023 Published: 8 January 2023

**Copyright:** © 2023 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/).

consumption accounted for 76% of industrial electricity consumption. In addition, research shows that 99% of environment-related problems in mechanical processes are due to electrical energy consumption [2]. Therefore, the establishment of an energy-saving machining system is an urgent requirement to reduce environmental impacts and every manufacturing enterprise needs to focus on it.

Energy-saving strategies using new materials and technologies may require enterprises to transform and invest a lot in existing manufacturing systems, therefore enterprises are usually inclined to carry out energy-saving scheduling and management [3]. Through scientific matching of production tasks and machine tools, more accurate calculation of tasks sequencing, reduce idle time of machine tools, and reasonable selection of machine tools on/off time can improve energy efficiency [4]. Additionally, Guzman et al. indicated that a gap still exists in developing mathematical models to deal with scheduling problems. Novel modeling approaches should be developed to address and associate the parameters related to production and sustainability [5], among which Feng et al. integrated multiple optimization algorithms and apply edge artificial intelligence (AI) to smart green scheduling of sustainable flexible shop floors [6]. Guzman et al. provided a mixed integer linear programming (MILP) model to address the multi-machine CLSD-BPIM (a capacitated lot-sizing problem with sequence-dependent setups and parallel machines in a bi-part injection molding) [7]. Mula et al. proposed a matheuristic algorithm to optimize the job-shop problem, which combines a genetic algorithm with a disjunctive mathematical model to cut computational times, and the Coin-OR Branch and Cut open-source solver is employed [8]. Rakovitis et al. developed a novel mathematical formulation for the energy-efficient flexible job-shop scheduling problem using the improved unit-specific event-based time representation and proposed a grouping-based decomposition approach to efficiently solve large-scale problems [9]. Knowing that approximately 80% of the energy consumption of machine tools is attributed to non-processing operations, whereas the actual energy consumed by processing operations accounts for less than 20% [10]. If only relying on advanced algorithms to achieve further energy saving in the workshop, the effect is limited. Wu and Sun realized energy saving by changing the turning on/off time of machine tools and choosing different machining speeds [11]. Gong et al. effectively reduced the number of machine restarts and total energy consumption by changing the start time of operations on different machines [12]. Cheng et al. proposed machine tool on/off criterion criteria, speed-scaling policy and transportation optimization strategy, and applied them to manufacturing unit scheduling problems to achieve overall energy saving [13]. An et al. proposed a worn cutting-tool maintenance strategy that reduced the impact of cutting-tool degradation and the total energy consumption of cutting-tool maintenance in manufacturing workshops [14]. Setiawan et al. studied a shop rescheduling problem caused by the failure or reduced service life of cutting tools [15].

As can be seen from the aforementioned literature, on the one hand, most energysaving scheduling problems start from the perspective of improving the performance of the algorithm, which makes the optimization calculation of shop energy consumption more accurate. However, on the other hand, from the perspective of workshop system management to achieve energy saving, in order to further realize the energy saving of the manufacturing system, it is important to consider the contribution of the coupling relationship between the energy consumption of individual equipment and the energy consumption of the system to the actual production and optimization objectives; however, this was almost ignored in previous studies. As the basic energy consumption equipment in the manufacturing process [16], the energy consumption caused by machine tools cannot be ignored. However, accurate estimation of energy consumption is the basis for improving energy efficiency. In recent years, different modeling methods for machine tool energy consumption have been proposed, such as those by He et al. and He et al., who combined the tool machining path with the energy consumption model to improve machining efficiency [17,18]. Shailendra et al. established an empirical model between cutting parameters and energy consumption of end turning by experiments [19]. Haruhiko et al. proposed

an empirical model for predicting machine tool power consumption based on the power function between specific energy consumption and material removal rate [20]. In addition, as the direct implementer of machine tool cutting, tool changing and maintenance will directly affect the production schedule, in order to avoid the tool suddenly reaching the end of life, resulting in the conflict of resources, energy consumption increase and the extension of the makespan and other problems. T. Mikołajczyk et al. and Sun et al. established the prediction model of tool residual life based on the historical data of tool wear [21,22]. M. Castejo'n et al. and P.J. Bagga et al. constructed a tool wear image dataset to predict tool life using cluster analysis and an exponential model [23,24]. Shi et al., Zhang et al. and Muhammad et al. introduced tool wear into the energy consumption model to achieve accurate energy consumption modeling, which laid a foundation for the integration of tool life prediction and energy consumption model [25–27]. Figure 1 summarizes the above literature.

**Figure 1.** Analysis of energy-saving scheduling research status.

In summary, few researchers combine shop scheduling under low-carbon production with single-machine tool energy consumption and tool life prediction. This paper analyzes the coupling relationship between shop scheduling, single-machine tool energy consumption and tool life prediction, and organically integrates the three to achieve deeper shop consumption reduction. Firstly, the machining power model and the tool life model of the machine tool were established through the tool wear-cutting experiment. Then, the two models were integrated into the shop scheduling system to obtain the machining power of each production procedure and the tool change time of each machine tool in the shop scheduling process, so as to realize the precise modeling of energy consumption at the system level. In addition, on the basis of the machine tool turn-on/off strategy of the workshop, considering the relationship between the tool change time and the turn-on/off time of the machine tool, the tool change time is adjusted to further reduce the machining power and the makespan of the workshop, so as to reduce the production energy consumption of the workshop, as shown in Figure 2.

The remainder of the paper is organized as follows. Section 2 describes the FJSP, cutting-tool degradation model, and hybrid energy-saving strategy of cutting-tool change and machine tool turn-on/off. In Section 3, a multi-objective optimization model of flexible job shop scheduling is established that considers tool degradation and energysaving measures. Section 4 introduces the proposed NSGA-II algorithm and its specific improvements. Section 5 sets the optimization model parameters through data collection and the analysis of actual cases. The rationality, effectiveness, and practical effects of the proposed model and algorithm are analyzed through verification experiments. Section 6 presents the conclusions and directions for future study.

**Figure 2.** Energy-saving scheduling research route.

#### **2. Problem Description**

The relevant symbols are provided in this section. Then, the FJSP that considers cuttingtool degradation with energy-saving measures (FJSP–CTD–ESM) is described. First, the FJSP is described. Then, the calculation method of the cutting-tool life, dynamic machining power, and the hybrid energy-saving strategy of cutting-tool change and machine tool turn-on/off is proposed, which combines the cutting-tool degradation and machine tool turn-on/off effects.

#### *2.1. FJSP Description*

In the FJSP, there are *n* kinds of jobs *J* = {*Ji*}*i*=1,2,...,*<sup>n</sup>* and *k* machine tools *M* = {*Mm*}*m*=1,2,...,*k*, and each job *Ji* has *Si* preset sequence of operations *<sup>O</sup>* <sup>=</sup> - *oi*,*j j*=1,2,...,*Si*

(Li et al., 2012). At least one operation *oi*,*<sup>j</sup>* in *O* can be processed by different machine tools, with a corresponding difference in the processing time and efficiency for the same operation.

The following conditions should be met in the FJSP: (1) A machine tool cannot be assigned to two or more operations simultaneously. (2) Each job has the same processing priority: initially, all jobs can be processed. (3) There is no constraint relationship between different jobs. (4) The optional machine tools for the job have no priority relationship. (5) All job processing tasks are non-preemptive. (6) The processing power of the machine tool and degree of cutting-tool wear obey the law of tool degradation. (7) Once a process begins, it cannot be interrupted before completion. Changing the cutting tool and turning the machine tool on/off cannot be inserted into the machining process. (8) The conversion time between different jobs with the same machine tool as well as the transportation time between different stages of the same job are ignored.

The symbols used in this paper are defined in Table 1.

**Table 1.** The symbols used in this paper.



**Table 1.** *Cont.*

#### *2.2. Cutting-Tool Degradation Model*

In the FJSP, the degradation of the cutting tool reduces its machining capacity, leading to an increase in the machining power and the interruption of the process caused by the blunt cutting tool. If the cutting-tool wear is considered in advance during the scheduling process, the change in machining power caused by cutting-tool wear can be accurately predicted. This not only improves processing efficiency and reduces energy consumption but also prevents the cutting tool from becoming blunt.

This section introduces the machining power model and cutting-tool life model derived from the tool degradation model.

(1) Machining power model

From the point of view of the working state of the machine tool, the machine tool power *Pm* in the workshop production process can be divided into two parts, as shown in Equation (1). The first part is the dynamic machining power of the machine tool *Pdm*, which includes the spindle power of the machine tool in the workpiece-cutting process. The second part is the static power *Psm* of the machine tool, including the no-load power of the motor and the power of the numerical control, lighting, and cooling systems.

$$P\_m = P\_{dm} + P\_{sm} \tag{1}$$

*Psm* exhibits little change in the machining process; hence, it is regarded as a constant value.

The dynamic power model [28,29] proposed by Tian et al. and Tian et al. is divided into two parts: the initial dynamic power without tool wear, and the additional dynamic power caused by tool wear, as shown in Equation (2):

$$P\_{dm} = K\_1 n\_v \, ^{a\_1} f^{a\_2} a\_p \, ^{a\_3} a\_c \, ^{a\_4} + K\_2 t\_m n \, ^{a\_5} f^{a\_6} a\_p \, ^{a\_7} a\_c \, ^{a\_8} \tag{2}$$

#### (2) Cutting-tool life model

To determine the relationship between the cutting-tool life and different cutting parameters, a type of cutting-tool failure should be selected as the criterion. According to ISO 8688-2:1989 *Tool life testing in milling-part 2: end milling* (1989), the wear of an end milling cutter can be divided into rake-face wear and flank-face wear [30]. Because the flank-face wear is easy to measure, the blunt standard of the tool wear is often set according to the maximum allowable value of the flank-face wear (usually expressed as VB). In this study, the end of the end milling cutting tool's life was defined as having a maximum VB of 0.3 mm in one of all teeth (VBmax = 0.3). The cutting-tool life model of Tian et al. and Sun et al. is shown in Equation (3) [22,29]:

$$T\_m = \mathbb{K}\_3 \cdot n\_v^{b\_1} f^{b\_2} a\_e^{b\_3} a\_e^{b\_4} \tag{3}$$

As we all know, tool wear is produced in complex mechanical and thermal environments, and there will be different dullness criteria for different processing objects or different quality requirements. In this paper, according to the dullness criterion mentioned in the ISO standard, in other application scenarios with higher cutting quality requirements, this part of modeling needs to establish a dullness criterion that meets the quality requirements and build models under these standards. This article only provides such a solution.

#### *2.3. The Hybrid Energy-Saving Strategy*

Section 2.2 shows that the cutting-tool life is not only related to the material and specifications of the cutting tool but also to the cutting parameters. Therefore, the cutting tool remaining useful life (RUL) cannot be calculated directly from the processing time of different operations. This results in a unique tool-changing schedule that affects the makespan and machine tool turn-on/off schedule.

Three measures are proposed to solve this problem.

(1) The cutting-tool change strategy

In this study, the cutting tool is changed before it is damaged to ensure that it meets processing quality requirements, reduces the risk of accidents, and improves the reliability of the processing system. Therefore, if the remaining service life of the cutting tool is insufficient to support the next processing task in the schedule, the cutting tool is considered unavailable and changed before the start of the next processing task, as expressed by Equation (4).

$$T\_m - t\_m < PT \tag{4}$$

Owing to the different cutting-tool lives under different cutting parameters, the cuttingtool service time cannot be added directly. A normalized approach is adopted to deal with this problem, that is, the increase in cutting-tool wear caused by the processing task is obtained using the processing task time/cutting-tool life under the cutting parameters of the task. Then, the total cutting-tool wear *Wm* is used to determine whether the cutting tool has reached the end point of its service life, as expressed by Equation (5).

$$\mathcal{W}\_m + PT \;/\; T\_m < 1\tag{5}$$

Equation (5) defines that cutting-tool wear must be less than 1.

(2) The machine tool turn-on/off strategy

During the production process, if a machine tool remains idle for some time, it is sensible to turn it off to avoid wasting energy and reduce carbon emissions. Turning machine tools on/off leads to additional energy consumption and could also damage their performance and service life. Therefore, the no-load balance time *TRm* should be set to control when to turn the machine tool on/off, as expressed by Equation (6). Meanwhile, to reduce the damage caused by turning the machine tool on/off, the interval between

on/off times is controlled by setting the on/off security threshold time *Hm*, as expressed by Equation (7) [11].

$$\mathcal{S}T\_{mr} - \mathcal{C}T\_{m(r-1)} \supseteq \; ^\*T\_{\mathcal{R}m} \eta\_{mr} \tag{6}$$

$$\left| \eta\_{mr} \mathcal{C} T\_{m(r-1)} - \eta\_{mr'} \mathcal{S} T\_{mr'} \right| \ge H\_m \delta\_{m,r} \quad r > r' \tag{7}$$

Equation (6) shows that if the time interval between two subsequent processing tasks is greater than the no-load balancing time, the machine tool should be turned off. Equation (7) shows that if the machine tool is turned on/off twice or more, the interval between the time the machine tool is turned on and the next turn-off must exceed the on/off security threshold time of the machine tool.

(3) Hybrid energy-saving strategy of cutting-tool change and machine tool turn-on/off

If the cutting-tool change operation is separated from the machine tool turn-on/off operation, the single cutting-tool change operation not only increases the makespan but also the standby energy consumption of the machine tool. Sacrificing a small amount of cutting-tool processing capacity by advancing the timing of the cutting-tool change will shorten the makespan and reduce energy consumption. Here, we set a reasonable coefficient *CFm* of the cutting-tool capacity to control when to change the cutting tool. Equation (8) ensures that the reduced processing capacity of the machine tool *Mm* is within an acceptable range. The cutting-tool change operation can be carried out before *x* tasks to realize energy savings.

$$CF\_m \ge 1 - \mathcal{W}\_m + \sum\_{r=0}^{\chi} \left(\lambda\_{mr} PT\_{mr} / \ T\_m\right) \tag{8}$$

#### **3. Formulation of FJSP–CTD–ESM**

In this section, the energy footprint model is defined. Then, the optimization model of the FJSP–CTD–ESM is established.

#### *3.1. Energy Footprint Model*

The total energy consumption *Etotal* of the workshop consists of five parts, as shown in Equation (9). Figure 3 shows the energy consumption of the machine tool at different stages [1].

$$E\_{\text{total}} = E\_{\text{C}} + E\_{\text{R}e} + E\_{\text{ct}} + E\_{\text{S}} + E\_{\text{Add}} \tag{9}$$

**Figure 3.** Schematic diagram of energy consumption distributions for different running states.

The energy consumption of each part is analyzed in detail below.

(1) Processing energy consumption *Ec* of machine tools

The processing energy consumption *Ec* of process *oi*,*<sup>j</sup>* is closely related to the cutting time *PTijm* and the machine tool power *Pijm*, which varies with the machine tool type, cutting parameters, and cutting-tool wear. The energy consumption of cutting is calculated as Equation (10).

$$E\_{\mathbb{C}} = \sum\_{m=1}^{k} \sum\_{i=1}^{n} \sum\_{j=1}^{S\_i} \sum\_{r=1}^{l\_m} \gamma\_{ijnm} P\_{ijm} P T\_{ijm} \tag{10}$$

(2) Energy consumption *ERe* of turning machine tools on/off

The energy consumption when turning the machine tool on/off *ERe* is influenced by the type and performance of the machine tool; it has no relation with the processing operation. Thus, the energy consumption of turning the machine tool on/off *ERem* once was set as a constant value. The energy consumption of machine tool turn-on/off is calculated as Equation (11).

$$E\_{R\varepsilon} = \sum\_{m=1}^{k} \sum\_{r=1}^{l\_m} \eta\_{mr} E\_{R\varepsilon m} \tag{11}$$

(3) Total energy consumption *Ect* of changing the cutting tool

As the processing time increases, it is necessary to change the cutting tool before it becomes blunt. The energy consumption of cutting-tool change is calculated as Equation (12).

$$E\_{ct} = \sum\_{m=1}^{k} \sum\_{r=1}^{l\_m} \lambda\_{mr} P\_{ctm} t\_{ctm} \tag{12}$$

(4) Standby energy consumption *Es*

When the machine tool is idle and kept on between two processes, it consumes energy while on standby, which is expressed as:

$$E\_s = \sum\_{m=1}^{k} \sum\_{r=1}^{l\_m} P\_{\rm sm} (1 - \eta\_{mr}) \left( S T\_{mr} - C T\_{m(r-1)} \right) \tag{13}$$

(5) Additional energy consumption *EAdd* of the workshop

The energy consumption of the workshop results not only from machine tool-related processes but also from lighting, computer utilization, and other sources. In this study, the additional energy consumption is not examined in detail; hence, it is set to a constant and calculated as Equation (14).

$$E\_{Add} = P\_{Add} \* \mathbb{C}\_{max} \tag{14}$$

where *PAdd* = 9.65 *Kw*.

#### *3.2. Formulation of the FJSP–CTD–ESM Optimization Model*

In actual production, the total energy consumption is not the only indicator; the makespan, total load of the machine tool, and the total number of times the machine tool is turned on/off and cutting tools are changed also need to be considered. Therefore, the multi-objective optimization model proposed in this paper has four objectives: the makespan *f*<sup>1</sup> (min), the total energy consumption of *f*<sup>2</sup> (Kw·min), the total load of machine tools *f*<sup>3</sup> (min), and the total number of times the machine tools are turned on/off and cutting tools are changed *f*<sup>4</sup> (time), expressed as Equations (15)–(35).

$$\min F = [f\_1, f\_2, f\_3, f\_4] \tag{15}$$

where

$$f\_1 = \mathbb{C}\_{\max} = \max\_{m,r} \mathbb{C}T\_{mr} \tag{16}$$

$$f\_2 = E\_{total} = E\_c + E\_{Rc} + E\_{ct} + E\_s + E\_{Add}$$

$$= \begin{cases} \sum\_{m=1}^{k} \sum\_{i=1}^{n} \sum\_{j=1}^{S\_i} \sum\_{r=1}^{l\_m} \eta\_{ijmr} P\_{ijm} P\_{ijm} + \\\ \sum\_{m=1}^{k} \sum\_{r=1}^{l\_m} \eta\_{mr} E\_{Rcm} +\\\ \sum\_{m=1}^{k} \sum\_{r=1}^{l\_m} \lambda\_{mr} P\_{ctm} t\_{ctm} +\\\ \sum\_{m=1}^{k} \sum\_{r=1}^{l\_m} P\_{sm} (1 - \eta\_{mr}) \left( S T\_{mr} - C T\_{m(r-1)} \right) + P\_{Add} \* C\_{max} \end{cases} (17)$$

$$f\_3 = WL = \sum\_{m=1}^{k} \sum\_{i=1}^{n} \sum\_{j=1}^{S\_i} \sum\_{r=1}^{l\_m} \gamma\_{ijmr} PT\_{ijm} \tag{18}$$

$$f\_4 = G = \sum\_{m=1}^{k} \sum\_{r=1}^{l\_m} (\eta\_{mr} + \lambda\_{mr}) \tag{19}$$

Subject to

$$CT\_{ij} - PT\_{ij} \ge CT\_{i(j-1)}\tag{20}$$

$$\begin{cases} \mathcal{S}T\_{m(r+1)} - \mathcal{S}T\_{mr} > 0\\ \mathcal{S}T\_{m(r+1)} - \mathcal{C}T\_{mr} \ge 0 \end{cases} \tag{21}$$

$$\sum\_{m=1}^{k} \gamma\_{ijmr} = 1\tag{22}$$

$$\text{ST}\_{mr} - \text{CT}\_{m(r-1)} \stackrel{>}{\geq} T\_{\text{Rm}} \eta\_{mr} \tag{2.3}$$

$$\left| \eta\_{mr} \mathcal{C} T\_{m(r-1)} - \eta\_{mr'} \mathcal{S} T\_{mr'} \right| \ge H\_m \delta\_{m'} \quad r > r' \tag{24}$$

$$N\_{mr} < 1\tag{25}$$

$$\text{ST}\_{mr} - \text{CT}\_{m(r-1)} \ge t\_{\text{ctm}} \lambda\_{mr} \tag{26}$$

$$CF\_{\rm m} \ge 1 - \mathcal{W}\_{\rm m} + \sum\_{r=0}^{\infty} (\lambda\_{\rm mr} P T\_{\rm mr} / \ | \ T\_{\rm m}) \tag{27}$$

*γijmr* = 1, 0, if the operation *oi*,*<sup>j</sup>* is the r th processing task of the machine tool *Mm* otherwise (28)

$$\eta\_{mr} = \begin{cases} 1, \text{ if the machine tool } M\_{\text{ill}} \text{ is turned on/off before its rth processing task} \\ 0, & \text{otherwise} \end{cases} \tag{29}$$

*λmr* = 1, 0, if the machine tool *Mm* changes the cutting tool before its r th processing task otherwise

$$
\delta\_{\mathfrak{m}} = \begin{cases} 1, & \text{if the machine tool } M\_{\mathfrak{m}} \text{ is turned on/off twice or more} \\ 0, & \text{otherwise} \end{cases} \tag{31}
$$

$$ST \ge 0\tag{32}$$

$$CT > 0\tag{33}$$

$$PT \ge 0\tag{34}$$

$$\forall i \in [1, n], j \in [1, s], m \in [1, k], r \in [1, l\_m] \tag{35}$$

Constraint [20] indicates that an operation cannot begin unless the preceding operation was completed. Constraint [21] confirms that a machine tool cannot be assigned to two or more processes simultaneously. Constraint [22] ensures that the same process cannot be conducted by two or more machine tools. Constraint [23] states that the machine tool turn-on/off does not overlap with the processing task, and the on/off time must exceed the no-load balance time. If it is the first round of processing, *CTijm*<sup>0</sup> = 0; otherwise, *CTijm*<sup>0</sup> = *CTijmlm* . Constraint [24] indicates that the interval between the time the machine tool is turned on and the next turn-off must exceed the machine tool on/off security threshold time *Mm*. Constraint [25] implies that the actual degree of tool wear *Wmr* must be less than 1. Constraint [26] requires that the cutting-tool change does not overlap with the processing task, and the time interval between the two processes is greater than the defined cutting-tool change time. Constraint [27] shows that the cutting-tool change of *Mm* is carried out in advance to achieve energy savings under the reduced processing capacity of the machine tool within an acceptable range. Constraints [28–35] are the constraints of decision variables.

#### **4. Proposed NSGA-II**

In this section, a general framework to solve the FJSP–CTD–ESM is proposed; the NSGA-II is briefly introduced, and the motivation behind the NSGA-II to optimize the FJSP–CTD–ESM is analyzed. Specific improvement measures of the NSGA-II algorithm are also described.

#### *4.1. Framework*

The optimization method of integrating the cutting-tool degradation, hybrid energysaving strategy, and production scheduling determines the cutting-tool capability and the order and priority of production tasks by combining machine tools and scheduling. The scheduling scheme is implemented in two steps: First, the machining power model and cutting-tool life model based on tool degradation were added to the fast and elitist multi-objective genetic algorithm (NSGA-II) scheduling algorithm [31] to generate an initial scheduling scheme, which includes the cutting-tool change. Second, through the scheduling mechanism, the cutting-tool change and machine tool turn-on/off were arranged to generate the final scheduling scheme to achieve the goal of energy conservation.

#### *4.2. Details and Improvements in NSGA-II*

#### 4.2.1. Scheduling Mechanisms

In the FJSP–CTD–ESM, processing energy consumption is mainly determined by machining power, which is also affected by cutting parameters and machine tool type. Non-machining energy consumption is mainly determined by the running state of machine tools. Because the cutting parameters were determined, the allocation of the appropriate machine tool to the operation and choosing the suitable machine tool turn-on/off and cutting-tool change times are the main considerations of the scheduling scheme. Two scheduling mechanisms are proposed:

(1) Cutting-tool degradation mechanism

During the processing operation, the cumulative processing time of the cutting tool is calculated in advance and the machining power *Pm* and degree of cutting-tool wear *Wm* of machine tool *Mm* are then calculated using the tool degradation model. When the degree of cutting-tool wear is expected to be greater than 1 (*Wm* > 1) after the machine tool completes the next operation, the cutting tool will be changed before processing (*λmr* = 1), and the machining power *Pm* is recalculated. The cutting-tool degradation mechanism is shown in Figure 4.

(2) Hybrid mechanism of cutting-tool change and machine tool turn-on/off

Based on the scheduling algorithm, the hybrid mechanism of cutting-tool change and machine tool turn-on/off is added to determine when to turn the machine tool on/off and change the cutting tool, as shown in Figure 5.

Step 1: According to the start and end time of the processing task, determine whether the non-processing time is greater than the no-load balancing time *TRm*. If so, go to Step 2; otherwise, end.

Step 2: Determine whether there is an on/off operation before this non-processing stage. If so, go to Step 3. Otherwise, turn off the machine tool in this non-processing stage and go to Step 5.

Step 3: Compare whether the time difference between the start of the non-processing stage and the last turn-on time of the machine tool is greater than the defined on/off security threshold time *Hm*. If the difference exceeds *Hm*, turn off the machine tool at the beginning of the non-processing stage and go to Step 5; otherwise, go to Step 4.

Step 4: Set the last turn-on time of the machine tool to *t*<sup>1</sup> and the start time of the non-processing stage to *t*2. If *t*<sup>2</sup> − *t*<sup>1</sup> − *Hm* > *TRm*, turn off the machine tool at (*t*<sup>1</sup> + *Hm*) and go to Step 5; otherwise, end.

Step 5: Check whether the machine tools have both cutting-tool change and turnon/off operations. If so, go to Step 6; otherwise, end.

Step 6: Determine whether the degree of tool wear *Wm* at the nearest machine tool on/off position before the cutting-tool change operation is greater than the difference between 1 and the cutting-tool capacity coefficient. If so, combine the cutting-tool change and machine tool turn-on/off operations, and recalculate the optimization target value; otherwise, end.

**Figure 4.** Cutting-tool degradation mechanism.

156

**,QSXW**,QLWLDO6FKHGXOLQJVFKHPH ݇ ܴܶ݉ ܨܥ݉ 6HHRULJLQDO16\*\$,,'HEHWDODQG6HFWLRQ **IRU** ݅ LQUDQJH݇ ܶݐݏ݅ 7XUQRQ7LPH6FKHGXOLQJVFKHPH7KHWXUQRQVWDUWWLPHVHWRI ܯ݅ ܵݐ݅ 6WDQGE\7LPH6FKHGXOLQJVFKHPH7KHVWDQGE\WLPHVHWRI ܯ݅ ܵݐݏ݅ 6WDQGE\VWDUW7LPH6FKHGXOLQJVFKHPH7KHVWDQGE\VWDUWWLPHVHWRI ܯ݅ **IRU** ݆ LQUDQJHOHQܵݐ݅ ܴܶ݉ !@݆<݅ݐܵ **LI** @[݆<݅ݐݏܵ ݔ LI݅ݐݏܶHQXPHUDWHLQ ݔ IRU < HVVBOLVW/ **LI**/HVHBOLVW >@ **LI** ܵݐݏ݅<݆@PD[/HVVBOLVW! ܪ݉ 6FKHGXOLQJVFKHPH &KDQJH6FKHPH6FKHGXOLQJVFKHPH ܶݐݏ݅ ܵݐ݅ ܵݐݏ݅0RGLI\WKHVFKHGXOLQJVFKHPHWRFKDQJHWKHVWDQGE\VWDWHWRWKHWXUQRII VWDWHDQG6HFWLRQ ܴܶ݉ !@݆<݅ݐݏܵ ݉ܪ HVVBOLVW/[PD@݆>݅ݐܵ **HOLI** 6FKHGXOLQJVFKHPH &KDQJH6FKHPH6FKHGXOLQJVFKHPH ܶݐݏ݅ ܵݐ݅ ܵݐݏ݅ ܪ݉ 0RGLI\WKHVFKHGXOLQJVFKHPHWRFKDQJHWKHVWDQGE\VWDWHWRWKHWXUQ RIIVWDWHIURPPD[/HVVBOLVW ܪ݉ DQG6HFWLRQ **HOVH** 6FKHGXOLQJVFKHPH &KDQJH6FKHPH6FKHGXOLQJVFKHPH ܶݐݏ݅ ܵݐ݅ ܵݐݏ݅ **IRU** ݅ LQUDQJH݇ ܥܶݐ݅ &KDQJHWRRO7LPH6FKHGXOLQJVFKHPH7KHFKDQJHWRROWLPHVHWRI ܯ݅ ܥܶݐݏ݅ &KDQJHWRRO7LPH6FKHGXOLQJVFKHPH7KHFKDQJHWRROVWDUWWLPHVHWRI ܯ݅ ܶݐ݅ 7XUQRQRII7LPH6FKHGXOLQJVFKHPH7KHWXUQRQRIIWLPHVHWRI ܯ݅ ܶݐݏ݅ 7XUQRQ7LPH6FKHGXOLQJVFKHPH7KHWXUQRQVWDUWWLPHVHWRI ܯ݅ @< ݅ݐܶ DQG>@ ݅ݐܶܥ **LI IRU** ݆ LQUDQJHOHQܥܶݐݏ݅ ܹ݉ 2EWDLQ:HDUܶݐݏ݅<݆@ ܥܶݐݏ݅<݆@2EWDLQWRROZHDUDWWKHQHDUHVWWXUQRQRIISRVLWLRQEHIRUHWKHFXWWLQJWRROFKDQJH ݅ݐܶܥ ! ݅ݐܶ DQG ݉ܨܥ  ܹ݉**LI** 6FKHGXOLQJVFKHPH &KDQJH6FKHPH6FKHGXOLQJVFKHPH ܥܶݐ݅ ܥܶݐݏ݅ ܶݐ݅ ܶݐݏ݅ 0RGLI\ WKH VFKHGXOLQJ VFKHPH WR FRPELQH WKH FXWWLQJ WRRO FKDQJLQJZLWKWKHWXUQRQRIIDQG6HFWLRQ ݅ݐܶܥ ݅ݐܶ DQG ݉ܨܥ  ܹ݉**HOLI** 6FKHGXOLQJVFKHPH &KDQJH6FKHPH6FKHGXOLQJVFKHPH ܥܶݐ݅ ܥܶݐݏ݅ ܶݐ݅ ܶݐݏ݅0RGLI\WKHVFKHGXOLQJVFKHPHWRGHOHWHWKHWXUQRQRIIDQGEULQJ WKHFXWWLQJWRROFKDQJLQJWRWKLVSRVLWLRQLQDGYDQFHDQG6HFWLRQ

**2XWSXW**6FKHGXOLQJVFKHPH

**Figure 5.** Hybrid mechanism algorithm form of cutting-tool change and machine tool turn-on/off.

#### 4.2.2. Encoding and Decoding with Changing Cutting and Turn-on/off

In the FJSP, encoding and decoding are expressed in the form of chromosomes, which are divided into the process chromosome and the machine tool chromosome. We utilize the encoding and decoding methods proposed by Zhang et al. [32], as shown in Figure 6.


**Figure 6.** Chromosome encoding process.

In addition, changing the cutting tool and turning the machine tool on/off affects the scheduling scheme, including the following six scenarios. (1) Scenario 1: If the machine tool *Mm* needs to be turned off before *oi*,*j*+1, the non-processing time and the last turn-on time of *Mm* should be considered. When the non-processing time is greater than the no-load balancing time *TRm*, and the difference between the start time *CTmr* of the non-processing time and the last turn-on time *STmr* of *Mm* is greater than the defined on/off security threshold time *Hm*, the machine tool is directly turned off between *CTmr* and *STm*(*r*+1), as shown in Figure 7. (2) Scenario 2: If *CTmr* − *STmr*- < *Hm* and *STm*(*r*+1) − *STmr*- − *Hm* > *TRm*, turn off *Mm* between (*STmr*-+ *Hk*) and *STm*(*r*+1), as shown in Figure 8. (3) Scenario 3. If machine tool *Mm* needs a cutting-tool change before *oi*,*j*<sup>+</sup>1, the non-processing time must be considered. If (*STm*(*r*+1)-*CTmr*) is greater than the cutting-tool change time *tctm*, change the cutting tool in *CTmr*, as shown in Figure 9. (4) Scenario 4: If the machine tool *Mm* needs a cutting-tool change before *oi*,*j*+<sup>1</sup> and *STm*(*r*+1)-*CTmr* < *tctm*, *STm*(*r*+1) needs to move to the right to *ST*- *<sup>m</sup>*(*r*+1) to make *ST*- *m*(*r*+1) -*CTmr* = *tctm*, as shown in Figure 10. (5) Scenario 5: If

the cutting-tool change and machine tool turn-on/off occurs in sequence, these should be combined, as shown in Figure 11. (6) Scenario 6: If the conditions of scheduling mechanism 2 are met, two cases will occur. Case 1: the cutting-tool change is incorporated into the machine tool turn-on/off. Case 2: Advance the cutting-tool change and eliminate the machine tool turn-on/off, as shown in Figure 12.

**Figure 7.** Scenario 1: (**a**) Initial Gantt chart; (**b**) Adjusted Gantt chart.

**Figure 8.** Scenario 2: (**a**) Initial Gantt chart; (**b**) Adjusted Gantt chart.

**Figure 9.** Scenario 3: (**a**) Initial Gantt chart; (**b**) Adjusted Gantt chart.

**Figure 10.** Scenario 4: (**a**) Initial Gantt chart; (**b**) Adjusted Gantt chart.

**Figure 11.** Scenario 5: (**a**) Initial Gantt chart; (**b**) Adjusted Gantt chart.

**Figure 12.** Scenario 6: (**1**) Sufficient machine-off time; (**2**) Insufficient machine-off time. (**a**) Initial Gantt chart; (**b**) Adjusted Gantt chart.

#### 4.2.3. Crossover

A crossover is used to maintain population diversity and explore a new solution space. In the parent population, two types of crossover were performed according to chromosome types [32], with the crossover of the process chromosome occurring in odd positions, while that of the machine tool chromosome occurring in even positions, as shown in Figure 13.

**Figure 13.** Crossover: (**a**) Crossover of the process chromosome; (**b**) Crossover of the machine tool chromosome.

#### 4.2.4. Mutation

A mutation is the use of the solution space to generate different neighborhood solutions through different mutations to prevent the population from falling into a local optimum during the process of population evolution convergence and increase the diversity of the population. The mutation can be divided into the process chromosome mutation and the machine tool chromosome mutation [32]. In the process chromosome mutation, gene *I* at a random mutation site on a random parental chromosome is mutated into the random gene *I* - . At the same time, one of the genes *I* on the other sites of this chromosome is randomly selected to become gene *I*, and the corresponding machine tool chromosome is changed accordingly. In the machine tool chromosome mutation, the genes at two random mutation sites on a random parental chromosome are mutated into any gene within the allowed mutation range, as shown in Figure 14.


**Figure 14.** Mutation: (**a**) Crossover of the process chromosome; (**b**) Crossover of the machine tool chromosome.

#### **5. Case Study**

In this section, the parameter determination method in FJSP–CTD–ESM is briefly described, and the comparison experiment based on the FJSP is designed. In addition, to provide enterprises with a better basis for selecting scheduling schemes, this paper integrates production functional resources with job-shop scheduling and provides production cost indicators to evaluate scheduling schemes. Thereafter, the rationality and effectiveness of the cutting-tool degradation model and hybrid energy-saving strategy of cutting-tool change and machine tool turn-on/off are illustrated by the example results.

#### *5.1. Design of Experiments*

#### 5.1.1. Environment Setting

All the algorithms were run in Python 2.7 on a personal computer with an Intel (R) Core (TM) i7-9750H, 2.6 GHz CPU and 8 GB RAM.

#### 5.1.2. Model Parameter Determination

To obtain stable machining power and tool wear data, the workpiece (i.e., 45# steel with dimensions 60 × 60 × 35 mm) was processed by milling. The cutting mode was straight-line face milling. The cutting path of the experimental workpiece is shown in Figure 15. The basic properties of the machine tool and tool used in the experiment are shown in Table 2. To establish the relationship between the machining power of machine tools and machining parameters, and that between the tool life and machining parameters, the machining power of machine tools and tool life under different cutting parameters were obtained by an orthogonal experiment. By referring to the Concise Manual of Cutting Parameters, it is found that in general, when the high-speed end milling cutter face-milling 45 steel, the recommended cutting speed is 21~40 m/min, and the recommended range of

feed per tooth is 0.12~0.2 mm. Therefore, within the recommended range, we consider the workpiece conditions, machine conditions, test costs and other factors, and through a large number of tests select the cutting parameters, as shown in Table 3 for the experiment. All the power data of a single milling cutter in a stable cutting period were collected for each group of experiments. The experimental devices are shown in Figure 16.

**Figure 15.** The tool cutting path.

**Table 2.** Basic parameters of experimental devices.


**Table 3.** Experimental cutting parameters table.



**Table 3.** *Cont.*

**Figure 16.** Experimental device: (**a**) Tool wear detection device; (**b**) Clamp on power logger (PW3360A982-04); (**c**) Tool wear detection process; (**d**) Power acquisition wiring diagram; (**e**) CNC machining center (M1/2,M3/4,M5/6); (**f**) Positioning device and fixture.

During the milling process, cutting-tool wear was detected once per ten cuts according to the cutting route. When the end point of the cutting-tool wear was reached, machining was stopped and the tool was changed. Photos of the detected partial cutting-tool wear are shown in Figure 17.

**Figure 17.** Partial tool wear detection photos.

The multivariate linear regression method [33] was used to determine the parameters in the cutting-tool degradation model, as shown in Figure 18. The machining power and cutting-tool life models are shown in Table 4.

**Figure 18.** Multivariate linear regression method.


**Table 4.** Cutting-tool degradation model.

The no-load balance time *TRm* of the machine tool was determined by measurement experiments and the on/off security threshold time *Hm* was obtained from the equipment manual.

The no-load balance time *TRm* must meet the following condition:

$$\begin{cases} \begin{array}{c} T\_{Rm} \ge \frac{\mathcal{R}T\_{\text{mean}}}{\overline{R}T\_{\text{sum}}} \ast SF\_m \end{array} \end{cases} \tag{36}$$

where *SFm* = 1.2.

5.1.3. Job and Workshop Configuration Information

The flexible job shop has six machine tools and five types of workpieces. The machining process information on the machine tools and workpieces are listed in Tables 5 and 6, respectively. The axonometric drawing of the five workpieces is shown in Figure 19.

**Table 5.** Machine tool information.


**Figure 19.** Three-dimensional model drawing of the five workpieces.


**Table 6.** Workpiece information.

#### 5.1.4. Establishment of Production Cost Indicator

Reducing cost plays an important role in the profitability, survival and development of an enterprise. How to get the maximum benefit with the minimum cost is an important topic that enterprises and even the whole society face and need to study and solve. To provide a better basis for selecting scheduling schemes, this paper establishes a mathematical model of production cost to evaluate scheduling schemes [34,35], as shown in Equation (36).

$$\text{COST} = \omega\_c E\_{\text{total}} + \omega\_m \text{WL} + \omega\_l G + \omega\_l \text{C}\_{\text{max}} \tag{37}$$

Table 7 shows the specific unit cost components in the production cost indicator, which includes unit energy cost, unit operating cost of machine tool, machine tool turn-on/off loss cost and cost per unit of labor time.

**Table 7.** Related unit costs.


#### *5.2. Evaluation*

#### 5.2.1. Experiment Results

(1) Hybrid mechanism of cutting-tool change and machine tool turn-on/off

It can be seen from Section 2.3 that machining power is not only related to tool wear but also related to cutting parameters. Therefore, to reflect and highlight the relationship between machining power variations and tool wear, the relationship between machining power and tool wear and the accuracy of the power model was analyzed by using actual and simulation data under the premise of certain cutting parameters.

It can be seen from the model that the machining power is linearly related to the tool utilization time. By comparing the actual collected power data with the tool utilization time and power data predicted by regression, it was found that the errors were all within a controllable range. The maximum error of the *M*1,*M*2; *M*3,*M*4; and *M*5,*M*<sup>6</sup> models are 6.44%, 3.36%, 4.67%, respectively, as shown in Figure 20.

**Figure 20.** Comparison of machining power and tool utilization time data: (**a**) *M*1,*M*2; (**b**) *M*3,*M*4; (**c**) *M*5,*M*6.

(2) Scheduling algorithm results

To reflect the performance of the improved NSGA-II (INSGA-II), scheduling solutions were generated under the premise that the degrees of cutting-tool wear of machine tools *M*1-*M*<sup>6</sup> are 60%, 70%, 50%, 50%, 70% and 60%, respectively. Table 8 summarizes the experimental results and shows the number of Pareto solutions, target values (*Cmax*, *Etotal*, *WL*, *and G*) and production cost indicator for each Pareto solution. The production cost is between CNY 36 and 40. The Pareto solution of the example is shown in Figure 21, where the X-, Y-, and *Z*-axis represents the *WL*, *Cmax*, and *Etotal*, respectively; the color represents *G*. As seen in Figure 21, the solution space is sufficient and the Pareto front is well distributed. The makespan is between 27 and 31 min, the total energy consumption of the workshop is between 319 and 350 *Kw*·*min*, the total machine tool load is 80–90 min, and the total number of times machine tools were turned on/off and cutting tools were changed is from five to nine times. The above data show that the INSGA-II algorithm can balance the four target values. The decision-maker can choose the optimal compromise using the multicriteria decision-making method.

**Table 8.** Cutting-tool degradation model.


**Figure 21.** Pareto solutions.

Figure 22 is a Gantt chart of Pareto solutions for this example. The *X*-axis represents the time, the *Y*-axis represents the machine number, and each block represents an operation. *O* denotes the machine tool processing operation, e.g., the first block *O*5,1 of *M*<sup>3</sup> indicates that the first process of *J*<sup>3</sup> was processed on *M*3. *CT* denotes the cutting-tool change. *Idle* indicates that the machine tool is not in use, e.g., the gray block on *M*6. *On*/*Off* indicates that the machine tool is turned off, e.g., the green block on *M*1. *R*\_*C*\_*T* indicates that the machine tool is off; however, the cutting tool will be changed at the beginning or end of this period, e.g., the gray-green block on *M*2. *S* indicates that the machine tool is on standby, e.g., the light blue block on *M*5. The scheduling scheme can provide the appropriate cutting-tool change time and turn off the machine tools on standby when necessary to save energy.

**Figure 22.** Gantt chart of an instance.

5.2.2. Assessing the Effects of the Cutting-Tool Degradation Model

In this section, based on the cutting-tool degradation model, multi-objective scheduling optimization is conducted starting with new cutting tools. Five schemes are selected from the Pareto optimal solutions for comparison, as follows: Scheme 1 includes the minimum sum of the four target values (*Cmax*, *Etotal*, *WL*, and *G*) with the weight of [0.3, 0.1, 0.5, and 0.1], and its four target values and cost are CNY (26.64, 309.7, 92.96, and 7) and 37.01. Scheme 2 includes the minimum makespan, and its four target values and cost are CNY (26.64, 309.7, 92.96, and 7) and 37.01. Scheme 3 includes the minimum total energy consumption of the workshop, and its four target values and cost are CNY (26.64, 309.7, 92.96, and 7) and 37.01. Scheme 4 includes the minimum total load of machine tools, and its four target values and cost are CNY (28.95, 328.13, 81.93, and 7) and 36.73. Scheme 5 includes the least number of times that the machine tools are turned on/off and cutting tools are changed, and its four target values and cost are CNY (36.82, 408.87, 90.61, and 5) and 40.94. Figure 23 shows the Gantt chart of production scheduling of Scheme 1. Figures 24 and 25 show the simulation power curve and degree of cutting tool wear curve of each machine tool, respectively, in the production process of Scheme 1. These reflect the change in the machining power of each machine tool with the processing operation, verify the influence of the cutting-tool degradation model on the total energy consumption of machine tools, and highlight the necessity of the tool degradation model.

**Figure 23.** Gantt chart of scheme 1.

**Figure 24.** Simulation machining power curve.

**Figure 25.** Scatter chart of tool wear.

From the cutting-tool degradation model, all machine tools can provide feedback on the real power condition of the machining process. Taking machine tool *M*<sup>1</sup> as an example, the curve of simulated and actual machining power can be obtained after actual machining, as shown in Figure 26. The simulated and actual average machining power of the three processes are 58.46 W, 69.49 W, and 65.24 W; and 56.7 W, 68.9 W, and 69.5 W, respectively. The power errors of the three processes are 3.15%, 0.85%, and 6.14%, respectively. The main reason for the considerable fluctuation in the actual machining power is that the cutting direction is not constant during machining. Changes in the cutting direction lead to an instantaneous power decrease because no cutting occurs at that moment, and the spindle generates a large amount of instantaneous power when it just touches the workpiece. The reason for the large error in the simulation results of the third process is that during hole enlargement, the actual cutting width is larger than the given cutting width (1.5 mm) because the actual processing path is circumferential; this makes the actual machining power larger than the simulation power. However, overall, the errors are all within the acceptable range and the cutting-tool degradation model can be applied in scheduling planning. This is conducive to making the simulation closer to the actual production, making the scheduling scheme and scheduling results more practical, and it plays a key role in predicting the machining power of machine tools.

**Figure 26.** Comparison curve between simulated power and actual machining power.

5.2.3. Assessing the Effects of the Energy-Saving Strategies

Comparing the results of the five scheduling schemes in Section 2.3 shows that adopting the machine tool turn-on/off measure reduced the average total machine standby time in the five schemes by over 99.2%: from 77.18 min to 35.40 s. The total energy consumption of turning machine tools on/off increased from 0 *Kw*·*min* to 1.744 *Kw*·*min*, whereas the standby energy consumption of the machine tools decreased from 31.891 *Kw*·*min* to 2.044 *Kw*·*min*. In addition, the cost of energy consumption decreased by about CNY 0.36. Although the machine tool turn-on/off strategy slightly increased the energy consumption when these are turned on/off, it significantly reduced the standby energy consumption, as shown in Figure 27. Therefore, it can be concluded that the machine tool turn-on/off energy-saving strategy is very effective.

**Figure 27.** Influence of the machine tool on/off strategy on non-processing energy consumption.

To assess the effects of the hybrid energy-saving strategy, we compared the energy consumption and makespan distinction before and after its adoption. Figure 28 shows the scheduling scheme changes before and after the hybrid energy-saving measure was applied. The four target values and cost before and after optimization are CNY 30.47, 347.18, 89.93, 10 and 39.92 and CNY 28.95, 331.76, 89.93, 10 and 38.97, respectively. It can be seen from the Gantt chart that after adopting the energy-saving strategy, the cutting-tool change of *M*<sup>3</sup> was performed before *O*5,1, resulting in a 3.95% reduction in the makespan from 30.14 min to 28.95 min, a 4.44% reduction in the total energy consumption of the workshop from 347.18 *kW* ·*min* to 331.76 *kW* ·*min* and a 2.44% reduction in the production cost from CNY 39.92 to 38.97, equivalent to CNY 47.3 saved every 24 h. However, the cutting-tool life of

*M*<sup>3</sup> was calculated as 143.002 min using the cutting-tool life model, while the processing time of *O*5,1 was 8.16 min, thus reducing the processing capacity of the cutting tool *DW* by 5.7%, as shown in Figure 29.

**Figure 28.** Comparison of scheduling schemes before and after the hybrid energy-saving strategy: (**a**) Before the hybrid energy-saving strategy; (**b**) After the hybrid energy-saving strategy.

**Figure 29.** Total impact of the hybrid energy-saving strategy.

The energy-saving strategy mainly reduced the energy consumption of the workshop by changing the timing of the cutting-tool change while reducing the makespan; the total machine tool load and the total number of times the machine tools were turned on/off and cutting tools were changed were not affected. The changes in the total energy consumption of the workshop were analyzed, as shown in Figure 30. In the solutions before and after optimization, the additional energy consumption of the workshop *EAdd* was the highest. This is because the energy consumption of a large number of additional equipment such as lighting, air pumps, and air conditioners in the actual workshop was far greater than that of the machine tools. Because this equipment is continuously operated, its energy consumption is positively correlated with time. As shown in Figure 31, the power consumption of *M*<sup>3</sup> changed when the cutting tool was changed in advance. The red and green parts represent the energy consumed before and after optimization, respectively. On the whole, the processing energy consumption *Ec* of the optimized solution decreased.

**Figure 30.** Effect of the hybrid energy-saving strategy on energy consumption.

**Figure 31.** Power contrast curve of *M*<sup>3</sup> after adopting the hybrid energy-saving strategy.

However, the reduction in *Ec* is much smaller than that in *EAdd* in terms of energy consumption and has little influence on the total energy consumption. The decrease in *EAdd* is mainly attributed to the integration of the cutting-tool change and machine tool turn-on/off processes, which leads to the reduction in the makespan, thus reducing the energy consumption of additional equipment. However, the cutting-tool change and machine tool turn-on/off processes were not eliminated; hence, their energy consumption did not change. In summary, the hybrid energy-saving strategy effectively reduced energy consumption and optimized the makespan, making it vital within the acceptable degree of reduction in the cutting-tool processing capacity.

#### **6. Conclusions**

Production planning and scheduling are usually the most critical activities in intelligent manufacturing enterprises. In the manufacturing process, manufacturers not only need to use the minimum resources to meet the production demand with as little energy consumption and in as short a time as possible but also face the challenge of the lack of mutual responsibility between the scheduling system and the single machine, which often leads to a large deviation between the scheduling optimization results and the actual application. Therefore, a new FJSP–CTD–ESM method is proposed in this paper to provide strong support for intelligent manufacturing enterprises to reduce the time and energy consumption in the production process. Through analyzing the coupling relationship between shop scheduling, single-machine tool energy consumption and tool life prediction, and organically integrating the three to achieve deeper shop consumption reduction. The resulting effect is as follows:

(1) Cutting-tool degradation during shop scheduling was analyzed. Based on the experimental data, exponential regression models of the dynamic power and cutting-tool life were established under certain machining conditions, with an error of approximately 6.5%. (2) A dynamic cutting-tool change strategy by monitoring the RUL was proposed to change the cutting tool before it becomes blunt. This makes the optimization model closer to the real machining situation. (3) Oriented towards low-carbon production objectives, the conventional machine tool turn-on/off schedule can reduce the non-processing energy consumption by 93.5%. Integrating the cutting-tool change strategy into the conventional machine tool turn-on/off schedule further reduces the total energy consumption by 4.44% and production cost by 2.44%. It was proved that this hybrid energy-saving strategy effectively reduces the energy consumption of workshops and has great application prospects.

In terms of the defects in this study, the proposed model does not consider the constraints of transport, clamping, and assembly on shop scheduling. In addition, the establishment of the machining power model and the tool life model of each machine tool in the workshop needs to spend a lot of time on the cutting wear experiment (about 45 h), which brings a lot of work to the preparation of the early production. When the shop changes the machine tool or changes a different type of tool, the models need to be rebuilt. Based on the above limitations, some suggestions are recommended as follows: (1) To explore a fast method to obtain the machining power model and the tool life model and make these models have a certain universal applicability. (2) To integrate more practical constraints such as transport, clamping, assembly, random breakdown or rush orders into the optimization model. (3) To design an efficient solution algorithm to solve multi-objective and many-objective optimization problems.

**Author Contributions:** Conceptualization, Y.T. and Z.G.; methodology, Y.T. and Z.G.; software, Z.G.; validation, Z.G.; formal analysis, Z.G.; investigation, Y.T., Z.G., L.Z., Y.C. and T.W.; resources, Y.T., L.Z. and T.W.; data curation, Z.G. and Y.C.; writing—original draft preparation, Y.T. and Z.G.; writing—review and editing, Y.T., Z.G. and T.W.; visualization, Z.G.; supervision, Y.T. and Z.G.; project administration, Y.T., L.Z. and Z.G.; funding acquisition, Y.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 51975407.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank the editors and anonymous reviewers whose comments helped greatly improve this paper. The work is supported by the National Natural Science Foundation of China (Grant No.51975407).

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

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Optimal Design of Reverse Logistics Recycling Network for Express Packaging Considering Carbon Emissions**

**Jia Mao 1, Jinyuan Cheng 1, Xiangyu Li 1, Honggang Zhao 2,\* and Ciyun Lin <sup>1</sup>**


**Abstract:** With the development of China's express delivery industry, the number of express packaging has proliferated, leading to many problems such as environmental pollution and resource waste. In this paper, the process of reverse logistics network design for express packaging recycling is given as an example in the M region, and a four-level network containing primary recycling nodes, recycling centers, processing centers, and terminals is established. A candidate node selection model based on the K-means algorithm is constructed to cluster by distance from 535 courier outlets to select 15 candidate nodes of recycling centers and processing centers. A node selection model based on the NSGA-II algorithm is constructed to identify recycling centers and processing centers from 15 candidate nodes with minimizing total cost and carbon emission as the objective function, and a set of Pareto solution sets containing 43 solutions is obtained. According to the distribution of the solution set, the 43 solutions are classified into I, II, and III categories. The results indicate that the solutions corresponding to Class I and Class II solutions can be selected when the recycling system gives priority to cost, Class II and Class III solutions can be selected when the recycling system gives priority to environmental benefits, and Class III solutions can be selected when the society-wide recycling system has developed to a certain extent. In addition, this paper also randomly selects a sample solution from each of the three types of solution sets, conducts coding interpretation for site selection, vehicle selection, and treatment technology selection, and gives an example design scheme.

**Keywords:** express packaging; green reverse logistics; reduce carbon emissions; K-means algorithm; bi-objective model; NSGA-II algorithm

**MSC:** 93-10

#### **1. Introduction**

In recent years, to promote the high-quality development of green couriers, green transportation, green consumption, and other areas, the concept of saving has been deeply rooted in people's hearts; green low-carbon mode of production and lifestyle is accelerating the formation. The report of the 20th National Congress of the Communist Party of China proposed to "accelerate the green transformation of the development mode", and such measures are especially evident in the express delivery industry. The 2021 government work report clearly pointed out that we should strengthen the construction of urban and rural circulation systems, especially to speed up e-commerce and express into the countryside, and expand consumption at the county and township level, which is the eighth time "express" was included in the government work report. According to the National Bureau of Statistics data, China's 2011–2021 express business volume increased significantly (see Figure 1). At the same time, along with the in-depth implementation of the new development pattern strategy, China's express industry will certainly enter a new stage of development shortly, bringing a new round of growth in the volume of express business and its revenue [1–4].

**Citation:** Mao, J.; Cheng, J.; Li, X.; Zhao, H.; Lin, C. Optimal Design of Reverse Logistics Recycling Network for Express Packaging Considering Carbon Emissions. *Mathematics* **2023**, *11*, 812. https://doi.org/10.3390/ math11040812

Academic Editor: Yaping Ren

Received: 11 January 2023 Revised: 29 January 2023 Accepted: 3 February 2023 Published: 5 February 2023

**Copyright:** © 2023 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/).

**Figure 1.** The 2011–2021 express business volume.

The monitoring data of the State Post Bureau showed that the country received 569 million express parcels on the Double Eleven, an increase of 28.54% year-on-year. However, the rapid development of the express industry has also brought about many social and environmental problems; especially the environmental pollution, management chaos, and waste of resources caused by express packaging waste are increasingly obvious [5–7]. Data show that in China's megacities, the increment of express packaging waste has accounted for 93% of the increment of domestic waste, and in some large cities 85% to 90%. How to effectively manage courier waste has become an urgent environmental problem to be solved. The main problems are manifested in the following areas.


(5) The reverse logistics system of express packaging recycling is not sound, and systematic scientific planning guidance is missing, manifested by the lagging work of express packaging classification, the confusion of social recycling channels and the low degree of specialization of treatment methods, the inadequacy of institutional mechanisms of relevant enterprises and government departments, the lack of laws and regulations and policy support, and the low enthusiasm of consumer participation in recycling [13–16].

Promoting the green transformation of express packaging, solving the bottlenecks faced by the industry, and achieving the sustainable development of the express industry is a complex and long-term systemic project that requires the participation and joint efforts of experts from different fields and different industry sectors [17]. In response to the above-mentioned problems in the express delivery industry, scholars have studied express packaging recycling from different perspectives.

In terms of site selection, Harsaj proposed a fuzzy multi-objective optimization model that quantifies three aspects simultaneously: economic, environmental, and social, and solved it using an improved particle swarm algorithm (PSO), and finally validated the model with a medical syringe recycling system [18]. Gao proposed a reverse logistics network design scheme based on a forward logistics network, and constructed an optimization model based on a multi-objective scenario with the objective functions of maximizing the expected total monetary profit, minimizing the expected total carbon emission cost, and maximizing the expected total job creation, and transformed it into a single-objective model to finally obtain the Pareto-optimal solution, and finally validated the effectiveness by using tires as an example [19]. Nie studied the supply chain configuration problem, constructed a mixed integer linear programming model with minimizing carbon emissions as the objective function, solved it using dynamic programming algorithms, and finally carried out an example verification to achieve a balance between economic and social benefits [20]. Guo studied fresh food distribution, built a two-stage model, considered the total system cost and vehicle path, and solved using a genetic algorithm and particle swarm algorithm, which effectively reduced carbon emissions and total cost [21]. Reddy constructed a multi-level multi-period mixed integer linear programming model with profit maximization as the objective function, considering the effects of facility location, vehicle type, and return rate, and finally gave an example analysis [22].

From the perspective of recycling model selection research, Liang used the Internet as a bridge to realize the design of a virtual APP and recycling device from the perspective of consumer psychology, real consumption situation, and the current situation of packaging recycling, to form a complete express packaging recycling system [2]. Yang constructed a multi-agent express packaging waste recycling system including the government, individuals, and enterprises. Based on differential game theory, the behavioral characteristics of individuals and the optimal strategies of government and enterprises under the market-driven recycling model, government-driven recycling model, and cooperativedriven recycling model were explored [23]. Based on previous studies by scholars, the main research contributions of this paper are as follows [24].


#### **2. Problem Description**

In this paper, according to the economic development and administrative division of Changchun, seven district administrative units under Changchun are selected as the area under study (hereinafter collectively referred to as M area), including Nanguan District, Kuancheng District, Chaoyang District, Erdao District, Lvyuan District, Shuangyang District, and Jiutai District (not considering Gongzhuling District). In this paper, we design the reverse logistics network for express packaging in region M. The regional map of region M is shown in Figure 2.

**Figure 2.** Regional map of M area.

According to the relevant data from Changchun Bureau of Statistics, this paper obtains the 2018–2020 population figures for the seven administrative regions mentioned above, as shown in Appendix A Table A1.

The proportion of the population of each district to the total population of Changchun was obtained (see Appendix A Table A2), in which the average proportion of the municipal districts to the total population of the city from 2018 to 2020 was the weighted average, and the weights of each year from 2018 to 2020 were 0.5, 0.3, and 0.2 according to the principle that the closer the year, the greater the weight.

In this paper, Baidu map API was used to obtain the original data of the latitude and longitude of courier points in the M area of Changchun and obtain 535 final valid data points after eliminating individual invalid data points, including 99 in Nanguan District, 81 in Kuancheng District, 107 in Chaoyang District, 103 in Erdao District, 81 in Lvyuan District, 17 in Shuangyang District, and 47 in Jiutai District. The latitude and longitude of express points in Nanguan District are shown in Appendix A Table A3, and the rest of the areas are omitted.

The relevant statistical information of Changchun Postal Administration was checked to obtain the express business volume in Changchun from 2013 to 2020, as shown in Table 1.


**Table 1.** The 2013–2020 express business statistics table in Changchun.

To meet the design requirements of the subsequent recycling network, the data in Table 1 need to be used to forecast the express business volume of Changchun in the next few years and select the appropriate value as the design value. In this paper, the least squares method is used to fit a linear function and an exponential function to the express business volume data of Changchun from 2013 to 2020, and the obtained functional relationships are shown in Equations (1) and (2).

$$y = 2242.9 \text{x} - 540.05 \tag{1}$$

$$y = 2183e^{0.2971x} \tag{2}$$

In the above Equations (1) and (2), *x* = *Year* − 2012. The fitted image is shown in Figure 3, and it can be visually seen that the fitting effect of Equation (2) is better than the fitting effect of Equation (1). Using Equations (1) and (2), the prediction of express business volume for 2013–2020 and 2021–2025 is shown in Table 2 and Figure 3.

**Figure 3.** Fitting graph of express business volume.


**Table 2.** Express business volume forecast.

As we can see in Table 3, the prediction effect of Equation (2) is better than that of Equation (1), so the express volume of 103,855,000 pieces in 2025 (x = 13) predicted by Equation (2) is selected as the design value.

**Table 3.** The parameters given in this paper and their values.


#### **3. Algorithm Introduction**


The K-means algorithm belongs to unsupervised machine learning and is a common classical clustering algorithm with the advantages of simplicity, efficiency, and ease of implementation, but it is sensitive to the selection of the initial clustering centers, which can affect the accuracy and speed of convergence if not selected properly. The K-means problem in a general sense can be described as follows: given a dataset containing N elements, each of which is an M-dimensional real vector, the objective is to select K points as clustering centers and divide the N elements into K sets, where each set corresponds to a clustering center so that the sum of the squared distances of all elements to the corresponding clustering center is minimized, at which point the clustering center of each set is the mean point of each set element [25].

3.1.2. K Value Determination Method


#### *3.2. Introduction of NSGA-II Algorithm*

3.2.1. Introduction to Multi-Objective Optimization Algorithms

When there are two or more objective functions, it is called multi-objective optimization.

The general form of multi-objective optimization is shown in Equation (3).

$$\min F = \left[ f\_1(\mathbf{x}), f\_2(\mathbf{x}), f\_3(\mathbf{x}), \dots, f\_m(\mathbf{x}) \right]^T \tag{3}$$

In the above equation, *F*(*x*) is the multi-objective optimization result, *f*1(*x*), *f*2(*x*), *f*3(*x*),... *fm*(*x*) is the objective component, and m is the objective dimension.

#### 3.2.2. Introduction to NSGA-II Algorithm

NSGA-II algorithm is a commonly used multi-objective genetic algorithm with the advantages of lower computational complexity and better population goodness and diversity. Its core is the introduction of fast non-dominated sorting, crowding degree, and elite strategy, as shown below [26].

#### (1) Introduction to the fast non-dominated sorting method.

Let *ni* denote the number of individuals dominating individual in the population, and *Si* is the set of individuals dominated by individual *i*. Find all individuals with *ni*= 0 in the population, i.e., the number of individuals dominating individual *i* is 0, i.e., individual *i* is not dominated, and deposit the eligible individual *i* into the non-dominated set R1, which means the subdominated rank is 1.

For all individuals j in the current non-dominated set R1, iterate through the set *Sj* of the individuals it dominates. Since the individuals *j* dominating individual t have been deposited in the current non-dominated set R1, the *nt* of each individual *t* in the set *Sj* is subtracted by 1. That is, the number of individuals governing the solution of individual t is reduced by 1. If *nt* − 1 = 0 is satisfied, then individual t is deposited in the set H.

R1 is used as the first level of the set of non-dominated individuals, and the individuals in this set are only dominated by other individuals and not by any other individuals, and all individuals in this set are assigned the same non-dominated ranking level, and then the above grading operation is continued for the set H, and the corresponding ranking level is also assigned, until all individuals are graded, i.e., all individuals are assigned the corresponding ranking level.

(2) Crowding degree profile.

The crowding degree *id* denotes the density of individuals around a given point in a population of a given generation, and in practice, it is measured by the length of the largest rectangle around individual *i* that contains individual *i* but not other individuals, where the crowding degree of individuals on each rank boundary is +∞. According to the definition of crowding degree, it can be seen that the larger the crowding degree is, the better the individuals are. The specific algorithm text is not repeated.

(3) Introduction to the elite strategy.

The elite strategy is to prevent the elimination of good individuals in the population in each generation, and to mix all individuals in the parent and child generations, and then select them according to the rank of non-dominance sorting and the size of crowding degree to get the new generation population that meets the population size requirement, effectively avoiding the loss of good individuals in the parent population, and the execution process is shown in Figure 4 [27].

**Figure 4.** Elite strategy diagram.

As shown above, firstly, the parents *Pt* and *Qt* are merged to obtain a new population with two times the original population size, and the individuals in the new population are sorted non-dominantly and selected according to the rank of non-dominant rank from smallest to largest until the selection reaches the rank of *Zi*, so that the number of selected individuals plus the number of individuals in this rank is greater than the original population size for the first time, and then the *Zi* is calculated. *Zi* rank in the crowding degree of individuals, and according to the size of the crowding degree from the largest to the smallest selection, until the population number requirement is met, so that the new generation of parents *P*- *t*.

#### **4. Model Building**

#### *4.1. Modeling of Reverse Logistics Network in M Region*

This paper determines the third-party logistics-centered express packaging model considering government participation, i.e., in the subsequent design of this paper, it is assumed that a third-party logistics enterprise specializing in express packaging recycling will carry out unified recycling and processing of express packaging of each courier company.

#### 4.1.1. Network Level and Node Analysis

According to the development status of the M region, this paper designs a fourlayer recycling network, and the schematic diagram of express packaging reverse logistics network layers and nodes in the M region is shown in Figure 5.

Among them, the first level is the primary recycling layer; the nodes of this level are the 535 courier points acquired, responsible for the recovery of express packaging directly from consumers, mainly by the consumers themselves to return express packaging, supplemented by door-to-door service for recycling. The second level is the recycling level; the node of this level is the recycling center, which is responsible for collecting and storing the express packaging recovered by the courier points within a certain area and connects to the primary recycling level with the relevant carriers leased or purchased, and the transportation process is short-distance transportation. The third level is the processing layer; the node of this level is the processing center, which is responsible for classifying the

express packaging recovered by each recycling center and carrying out different technical treatments according to different categories, mainly including reuse treatment and transport to each recycling center, transport to the paper mill, and transport to landfill, with the relevant carrier leased or purchased to connect the recycling layer and the terminal layer, the transport process is medium and long-distance transport compared with the transport process of the first and second levels. The transportation process is medium to long distance compared to the first and second levels. The fourth level is the terminal level, the nodes of this level are the paper mill and the landfill, which are responsible for accepting the express packaging after sorting in the processing center, the longitude and latitude of the paper mill and the landfill are (43.910551◦ N, 125.420776◦ E) and (43.964575◦ N, 125.358692◦ E), and they are connected to the processing level by the relevant carriers leased or purchased. The process is also medium to long distance.

**Figure 5.** Schematic diagram of express packaging reverse logistics network hierarchy and nodes in M region.

4.1.2. Determination of Candidate nodes Based on the K-Means Algorithm

• Basic assumptions

(1) Euclidean distance is used in this paper to calculate the distance between the data and the center of clusters (center of mass).

(2) It is assumed that the influence of the Earth's surface on the distance calculation is negligible in the range of M region.



• Model Building

$$SSE = \sum\_{i=1}^{K} \sum\_{\mathbf{x}\_{i} \in \mathbf{S}\_{i}} \left| \left| \mathbf{x}\_{i} - \mu\_{i} \right| \right|^{2} \tag{4}$$

$$\min V = \sum\_{\mathbf{i}=1}^{\mathbf{K}\_0} \sum\_{\mathbf{x}\_{\mathbf{j}} \in \mathbf{S}\_{\mathbf{i}}} \left| \left| \mathbf{x}\_{\mathbf{j}} - \boldsymbol{\mu}\_{\mathbf{i}} \right| \right|^2 \tag{5}$$

	- (1) Determine K0 (see Figure 6)

**Figure 6.** Determining K0 in the K-means algorithm.

(2) K-means clustering (see Figure 7)

**Figure 7.** K-means algorithm specific steps.

#### *4.2. NSGA-II Algorithm-Based Node Siting Modeling*



• Symbol Description


• Objective function

In this paper, the objective function is set from two perspectives: economic and environmental.


$$
min F = f\_1 + f\_2 + f\_3 \tag{6}
$$

$$
min G = \mathcal{g}\_1 + \mathcal{g}\_2 \tag{7}
$$

$$f\_1 = \sum\_{i=1}^{N} \sum\_{j=1}^{N} X\_i X\_{ji} d\_{ji} V\_{ji} m\_0 t\_{rk} + \sum\_{i=1}^{N} X\_i d\_{im} V\_{im} m\_0 t\_{rk} + \sum\_{i=1}^{N} X\_i d\_{in} V\_{in} m\_0 t\_{rk} + \sum\_{i=1}^{N} \sum\_{j=1}^{N} X\_i X\_{ij} d\_{ij} V\_{ij} m\_0 t\_{rk} \tag{8}$$

$$f\_2 = \begin{cases} a\_j V\_i m\_\prime \, V\_i m \le V\_{\max} \, m\_0 \\ +\infty, \quad V\_{\max} \, m\_0 < \, V\_i \, V\_i \end{cases} \tag{9}$$

$$f\_3 = \begin{cases} \sum\_{\substack{i=1\\N\\i=1\end{r}}^N X\_i V\_i m \mathbb{C}\_1, \ i \in \{1, 2, 3\} \\\ \sum\_{i=1}^N X\_i V\_i m \mathbb{C}\_2, \ i \in \{4, 5, 6\} \\\ \sum\_{i=1}^N X\_i V\_i m \mathbb{C}\_3, \ i \in \{7, 8, 9, 10, 11, 12, 13\} \\\ n \sum\_{i=1}^N X\_i V\_i m \mathbb{C}\_4, \ i \in \{14, 15\} \end{cases} \tag{10}$$

$$\mathcal{G}\_{1} = \sum\_{i=1}^{N} \sum\_{j=1}^{N} X\_i X\_{ji} d\_{ji} V\_{ji} m\_0 L\_k \varepsilon\_{r\bar{k}} + \sum\_{i=1}^{N} X\_i d\_{i\bar{n}} V\_{i\bar{m}} m\_0 L\_k \varepsilon\_{\bar{k}\bar{r}} + \sum\_{i=1}^{N} X\_i d\_{i\bar{n}} V\_{i\bar{n}} m\_0 L\_k \varepsilon\_{r\bar{k}} + \sum\_{i=1}^{N} \sum\_{j=1}^{N} X\_i X\_{\bar{i}j} d\_{\bar{i}\bar{j}} V\_{\bar{i}\bar{j}} m\_0 L\_k \varepsilon\_{\bar{k}\bar{r}} \tag{11}$$

$$L\_k = L\_k^0 + \frac{L\_k^\* - L\_k^0}{N} m\_0 \tag{12}$$

$$\lg\_2 = V\_1 m\_0 c\_{arj} \tag{13}$$

• Constraints

$$X\_{l} = \begin{cases} \text{1, } \text{The ith candidate node is selected as the processing center} \\ \text{0, The ith candidate node is selected as the reciprocing center} \end{cases} \tag{14a}$$

$$X\_{ji} = \begin{cases} \text{1, } \text{The jth candidate node is assigned to the ith candidate node} \\ \text{0, The jth candidate node is not assigned to the ith candidate node} \end{cases} \tag{14b}$$

$$0 \le V\_i \le V\_{\max} \tag{15a}$$

$$V\_{\bar{i}} = V\_{\bar{i}\bar{i}} + V\_{\bar{j}\bar{i}} \tag{15b}$$

$$V\_i = V\_{in} + V\_{in} \tag{15c}$$

Equation (6) represents the cost, including transportation cost, construction cost, and treatment cost, in which, transportation cost includes three parts from recycling center to treatment center, treatment center to landfill and paper mill, and treatment center back to the recycling center, construction cost considers four types of areas A, B, C, and D and the cost per unit construction area is different for candidate nodes in different areas, and treatment cost considers alternative of three technologies, and the treatment cost of different technologies is different. Equation (7) represents carbon emissions, including carbon emissions from the transportation process and carbon emissions from the treatment process, where carbon emissions from the transportation process include three parts from recycling center to treatment center, from treatment center to landfill and paper mill, and from treatment center back to the recycling center, and carbon emissions from treatment process consider the three alternative technologies, and the carbon emissions generated by different technologies are different. The above Equation (15a) indicates that the storage volume of node *i* takes a range of values, Equation (15b) indicates that node *i* is equal to its storage volume plus the volume transported from point *j* to point *i*, Equation (15c) indicates that the volume transported out of node *i* is not greater than the storage volume of node *i*. Equations (15a) and (15b) indicate that the above two equations indicate the conservation of flow, the way to achieve each of the above constraints, especially the capacity constraint through the constraint matrix.

According to the relevant information and combined with the actual situation, the following values of the relevant parameters are given in this paper. The values of each variable are as follows in Table 3.

Take fuel consumption per unit distance at no load L0 <sup>1</sup>= 0.11 L/km, L<sup>∗</sup> <sup>1</sup>= 0.15 L/km at full load, carbon emission factor per unit fuel consumption ctr1 = 2.5 kg/L, freight per unit tr1 = 0.001 CNY/km/kg in vehicle type 1, fuel consumption per unit distance at no load L0 <sup>2</sup>= 0.1 L/km, L<sup>∗</sup> <sup>2</sup>= 0.15 L/km at full load. The carbon emission factor per unit of fuel consumption ctr2 is 2.8 kg/L, the freight cost per unit tr2 is 0.0008, and the electricity consumption per unit distance *L*<sup>0</sup> <sup>3</sup> is 0.1 in vehicle type 2. When vehicle type 3 is empty, L∗ <sup>3</sup>= 0.3, when fully loaded, the carbon emission factor per unit of electricity consumption ctr3 is 0.8, the freight cost per unit tr3 is 0.0015, and the electricity consumption per unit distance L0 <sup>4</sup>= 0.1. When vehicle type 4 is empty, L<sup>∗</sup> <sup>4</sup>= 0.3 when fully loaded, the carbon emission factor per unit of electricity consumption ctr4 = 0.1, unit freight tr4 is 0.0012.

#### **5. Algorithm Design**

#### *5.1. NSGA-II Algorithm Description*

(1) Chromosome coding

In this paper, we set each generation of the population containing n individuals, and each individual has only one chromosome, using a repeatable integer coding method, and the total length of the chromosome is 75. From the perspective of corresponding

functions, each chromosome can be divided into five parts, and the length of each part is 15. Specifically, the first part indicates the site selection code, the second part indicates the vehicle type selection code from the recycling center to the treatment center, the third part indicates the vehicle selection code from the processing center to the recycling center, the fourth part indicates the vehicle selection type code from the processing center to the paper mill and the landfill, and the fifth part is the processing technology selection code.

• First part: site selection code

For a generational population, the first part of the chromosome of the ith (1 ≤ *I* ≤ n, and *i* is an integer, the same below) individual, the *j*th (1 ≤ *j* ≤ 15, and *j* is an integer, the same below) position indicates *j* candidate nodes and the value *k* corresponding to the *j*th position indicates that the *j*th candidate node is assigned to node k, and *k* becomes the processing center, where *k* ∈ {1 ≤ *k* ≤ 15, and *k* is an integer}.

For example, Figure 8I represents the first part of the chromosome of the first individual of a generation population, whose length is 15, and the value corresponding to the first position is 1, indicating that node 1 is assigned to node 1, i.e., node 1 is selected as a processing center by the candidate node, and so on, and the value corresponding to the 15th position is 15, indicating that node 15 is assigned to node 15, i.e., node 15 is selected as a processing center.


**Figure 8.** (**I**,**II**) indicate the site selection code, (**III**,**IV**) indicate the vehicle type selection code from the recycling center to the processing center, the vehicle selection code from the processing center to the recycling center and the vehicle type selection code from the processing center to the paper mill and landfill, (**V**,**VI**) are the processing technology selection code.

Figure 8II represents the first part of the chromosome of the 7th individual of this generation population, whose length is 15, and the 1st, 2nd, and 3rd positions correspond to values all of 5, indicating that nodes 1, 2, and 3 are assigned to node 5, i.e., node 5 is selected as a processing center by the candidate node, and accordingly, nodes 1, 2, and 3 are selected as a recycling center by the candidate node, and so on, and the 12th, 13th, 14th, and 15th positions correspond to the value 13, indicating that nodes 12, 13, 14, and 15 are assigned to node 13, i.e., node 13 is selected as the processing center by the candidate node, and nodes 12, 14, and 15 are selected as the recycling center by the candidate node.

• Second, third, and fourth parts.

For the second part of the *i*th (1 ≤ *i* ≤ *n* and *i* is an integer, same below) individual chromosome of a generational population, the value *l* corresponding to the *j*th (16 ≤ *j* ≤ 30 and *j* is an integer, same below) position indicates the type of transport vehicle between the (*j*-15)th candidate node and the kth node corresponding to the (*j*-15)th position in the first part, where *l* ∈ {1 ≤ *l* ≤ 4 and *l* is an integer } and *k*∈{1 ≤ *k* ≤ 15 and *k* is an integer}.

For example, Figure 8III represents the second part of the chromosome of the first individual of a generational population which is the same population as the population of Figure 8I with a length of 15. The 16th position corresponds to a value of 1, indicating that the type of the transport vehicle between the first candidate node and the first candidate node corresponding to the first position in the first part is type 1, and so on, and the 30th position corresponds to the value of 4 indicates that the type of the transport vehicle between the 15th candidate node and the 15th candidate node corresponding to the 15th position in the first part is type 4.

Figure 8IV represents the chromosome of the 7th individual of this generation population with a length of 15, and the values corresponding to the 16th, 20th, 24th, and 28th positions are all 1, indicating that the type of transport vehicle between the 1st, 5th, 9th, and 13th candidate nodes and the 5th, 6th, 11th, and 13th candidate nodes corresponding to the first part is all type 1, and the values corresponding to the 17th, 21st, 25th, and 29th positions are 2, indicating that the vehicle types of the transport vehicles between the 2nd, 6th, 10th, and 14th candidate nodes and the corresponding 2nd, 6th, 10th, and 14th candidate nodes in the first part are all of type 2, and the remaining cases and so on.

• Fifth part.

For the *i*th (1≤ *i* ≤ n, *i* is an integer, the same below) individual chromosome of the fifth part of a generational population, the value *p* corresponding to the *j*th (61 ≤ *j* ≤ 75, *j* is an integer, the same below) position denotes the *p*th technique chosen for the *k*th node corresponding to the *j*-15th position of the first part, where *p*∈{1≤ *l* ≤ 3 and *p* is an integer} and *k* ∈{1 ≤ *k* ≤ 15 and *k* is an integer}. For example, Figure 8V represents the fifth part of the chromosome of the first individual of a generational population (this population is the same population as the population in Figure 8), whose length is 15, and the 61st position corresponds to a value of 1, indicating the *j*60th candidate node, i.e., the first candidate node uses technology 1, and so on, and the 75th position corresponds to a value of 1, indicating the *j*60th candidate node i.e., the 15th candidate node also adopts technique 1.

Figure 8VI represents the chromosome of the seventh individual of this generation population with a length of 15, and the values corresponding to the 61st, 64th, 65th, 69th, 72nd, and 73rd positions are all 1, indicating that the 1st, 4th, 5th, 9th, 12th, and 13th candidate nodes corresponding to the first part adopt technology 1, and the values corresponding to the 62nd, 66th, 70th, and 74th positions are all 2, indicating that the 2nd, 6th, 10th, and 14 candidate nodes adopt technique 2, and the values corresponding to the 63rd, 67th, 68th, 71st, and 75th positions are all 3, indicating that the 3rd, 7th, 8th, 11th, and 15th candidate nodes adopt technique 3.

(2) Population initialization

According to the above coding method, in this paper, let there be *n* = 50 individuals per generation of population, each individual has only one chromosome, and the length of each chromosome is 75, i.e., m = 75, the first part corresponds to the candidate processing center coding, each position randomly generates a repeatable integer *k* (1≤ *k* ≤15), the second, third, and fourth parts correspond to the vehicle type coding, each position randomly generates a repeatable integer *l* (1 ≤ *l* ≤ 4), the fifth part corresponds to the type of technology used, *p* ∈ {1 ≤ *l* ≤ 3, and *p* is an integer}.

(3) Adaptation degree function

In this paper, the fitness of each individual is calculated according to the rank size and crowding degree of the non-dominated stratum, specifically, the parents and children are merged, the new population after the merger is non-dominated stratified, and the crowding degree is calculated for all individuals, and finally, the individuals are selected according to the principle that priority is given to individuals with small non-dominated stratification rank and priority is given to individuals with large crowding degree in the same stratum until the population number is satisfied requirements, the method has a slightly different logic from the classical NSGA-II algorithm, but is identical in purpose and fully equivalent in effect [32–34].

(4) Crossover operation [35]

In this paper, the simulated binary crossover method is used to perform crossover operations on population chromosomes. Assuming that the children generated by the crossover of parents *xa* and *xb* are *ya* and *yb*, then for the *kth* position of children *ya* and *yb* we have.

$$\mathbf{x}\_d(k) = \frac{1}{2} [(1+\beta)\mathbf{x}\_d(k) + (1-\beta)\mathbf{x}\_b(k)] \tag{16}$$

$$y\_b(k) = \frac{1}{2} [(1 - \beta)\mathbf{x}\_a(k) + (1 + \beta)\mathbf{x}\_b(k)] \tag{17}$$

Among them,

$$\beta = \begin{cases} \, \, 2r^{\frac{1}{1+\eta}}, r \le 0.5\\ \, \, (2-2r)r^{-\frac{1}{1+\eta}}, r > 0.5 \end{cases} \tag{18}$$

In the above equation, *r* ∼ *U* [0, 1], *η* is a custom parameter, and the larger the value, the closer the offspring is to the parent.

In this paper, we take *η* = 20, and for crossover, the first half of individuals and the second half of individuals in each generation of the population are combined two by two, and when the number of individuals is odd, the last individual does not participate in the crossover.

(5) Variation operation

In this paper, the polynomial variation method is used to perform various operations on population chromosomes; specifically, the variation form is,

$$
\boldsymbol{v}\_{k}^{l} = \boldsymbol{v}\_{k} + \delta(\boldsymbol{u}\_{k} - \boldsymbol{l}\_{k}) \tag{19}
$$

Among them,

$$\delta = \begin{cases} \begin{bmatrix} 2u + (1 - 2u)(1 - \delta\_1)^{\eta\_m + 1} \end{bmatrix}^{\frac{1}{\eta\_m + 1}} - 1, \ u \le 0.5\\ 1 - \left[ 2(1 - u) + 2(u - 0.5)(1 - \delta\_2)^{\eta\_m + 1} \right]^{\frac{1}{\eta\_m + 1}}, \ u > 0.5 \end{cases} \tag{20}$$

In the above equation, *δ*<sup>1</sup> = (*vk* − *Ik*)/(*uk* − *Ik*), *δ*<sup>2</sup> = (*vk* − *vk*)/(*uk* − *Ik*), *u* is a random number in the interval [0, 1], *η<sup>m</sup>* is a custom parameter, and this paper takes *ηm*= 20.

#### *5.2. NSGA-II Algorithm Steps*

Step 1: encoding by repeatable integer coding.

Step 2: initialize the population and generate a population containing m individuals, each containing one chromosome, at this point, set as the initial population.

Step 3: non-dominated stratification of the individuals of the initial population.

Step 4: calculate the fitness of the individuals of the initial population based on the results of the non-dominated stratification in step 3.

Step 5: select a certain number of individuals in the initial population as the evolutionary generation 0 according to the calculation result in step 4, and all individuals of the initial population are selected as generation 0 in this paper.

Step 6: start evolution, take generation *i* (*i* ∈ {0 ≤ *i* < 50 and *k* is an integer}) as the parent, perform crossover and mutation operations on the individuals of the parent to generate the children corresponding to the parent of generation i, and fuse the individuals of the parent and children of generation *i*.

Step 7: decode the fused individuals from step 6 and calculate the objective function values of the fused individuals.

Step 8: non-dominated stratification of the fused individuals from step 6 and calculation of their crowding degree.

Step 9: according to the calculation result of step 8, select m individuals as the (*i* + 1)*th* generation according to the non-dominated stratification level from the lowest to the highest and when the same level according to the crowding degree from the largest to the smallest, and return to step 6 if the required number of evolutionary iterations is not satisfied.

Step 10: the (*i* + 1)*th* generation of individuals has been noted as the Pareto solution and the corresponding solution is the Pareto frontier solution [36–38].

The algorithm terminates.

The basic flow chart of the algorithm is shown in the following Figure 9.

**Figure 9.** The basic process of mutation operation.

#### **6. Analysis of Results**

*6.1. Candidate Points*

(1) Execution of the algorithm yields the SSE versus K plot, as shown in Figure 10.

**Figure 10.** SSE versus K.

From the above Figure 8, it can be seen that the inflection point appears between K = 5 and K = 10. It is known from the rule of elbow law that K0 = 15 should satisfy 5 ≤ K0 ≤ 10, but considering the actual demand, K0 can be expanded appropriately by the rule of on-demand selection law, and K0 = 15 is taken in this paper.

(2) Take K0 = 15, execute K-means algorithm, get 15 clustering centers, and use them as candidate processing centers, whose latitude and longitude information is shown in Table 4, and the location schematic is shown in Figure 11.

**Figure 11.** Schematic diagram of the location of the candidate nodes, paper mill, and landfill.


**Table 4.** Latitude and longitude of candidate nodes.

6.1.1. Classification of Candidate Nodes in Region M

Considering the economic development status of seven districts, this paper delineates four categories of regions A, B, C, and D, as shown in Figure 12.

**Figure 12.** Schematic diagram of candidate node sub-region location.

According to the above classification results, the regions to which each candidate node belongs are shown in Table 5.

**Table 5.** Classification regions and administrative regions to which each candidate node belongs.


#### 6.1.2. Determination of The Distance between Candidate Nodes

The latitude and longitude of 15 nodes can be known from Table 3, and the distance between any two points (*L*1, *N*1) and (*L*2, *N*2) is calculated using Equation (20) to obtain the distance between each candidate node, as shown in Table 6.

$$d = \sqrt{\left[111(N\_1 - N\_2)\right]^2 + \left\{111[E\_1 \cos(N\_1) - E\_2 \cos(N\_2)]\right\}^2} \tag{21}$$

**Table 6.** Distance between candidate nodes.


6.1.3. Determination of the Candidate Node Express Volume

According to the previous data, the design values were weighted according to the average ratio of the population in the seven districts of the M region to the total population of Changchun, and the results are shown in Table 7.

The classification results of the candidate nodes, Nanguan District has one class A and C regional nodes, respectively, candidate nodes 2 and 12, known from Table 1 Nanguan District express the business volume of 103,352,300 pieces. This paper assumes that the ratio of A and C regional nodes express business volume in Nanguan District is 7: 3, then the express business volume of candidate nodes 2 and 12 are 7234.65967 and 3100.56843 million pieces. According to the above rules, the express business volume of each candidate node is shown in Table 8.


**Table 7.** Express business volume by district.

**Table 8.** Express the business volume of each candidate node.


#### *6.2. Number of Iterations, Crossover Variance Probability Selection*

In this paper, we use Python 3.8.5 for algorithm implementation, in which some subfunctions directly call the functions of the Geatpy library, such as selection sub-functions and crossover sub-functions.

In this paper, we set the number of population individuals m = 100 and M = 3000, in the actual problem, the number of iterations will have a large impact on the performance of the whole algorithm, so in this paper, the two single objectives of total cost and carbon dioxide total emission, as shown in Figure 13, where the blue line indicates the average value and the red line indicates the minimum value, it can be seen that both the average value and the minimum value decrease gradually with the increase of the number of generations and tend to be stable, and both objective functions show good convergence, so it is feasible to take M = 3000.

**Figure 13.** Plot of the total cost–number of iterations, total carbon dioxide emissions–number of iterations c (M = 3000).

In this paper, we compare and analyze the relevant indicators for four cases with cross-variance probabilities of 0.9, 0.8, 0.7, and 0.6. Total cost vs. total carbon dioxide

emissions relationship diagram (cross-variance probability is 0.9, 0.8, 0.7, 0.6) as shown in Figure 14A–D and Table 9.

**Figure 14.** (**A**) Total cost vs. total carbon dioxide emissions relationship diagram (cross-variance probability is 0.9). (**B**) Total cost and total carbon dioxide emissions relationship diagram (crossvariance probability is 0.8). (**C**) Total cost versus total carbon dioxide emissions relationship diagram (cross-variance probability is 0.7). (**D**) Total cost and total carbon dioxide emissions relationship diagram (cross-variance probability is 0.6).



Note: After testing, each index fluctuates somewhat during repeated calculations, and the best value in each calculation process is taken.

Integrating the three indicators of non-dominated solution percentage, HV, and Spacings, and the final three generated images, this paper selects the case of M = 3000 and the probability of cross-variance is 0.7 for analysis, and pools the analysis results to give management suggestions. According to the results of the selection of the number of iterations, this paper runs the procedure at M = 3000 and the cross-variance probability is 0.7, and a total of 43 Pareto solutions are obtained, as shown in Figure 15, which can be more clearly seen in the Pareto frontier solutions [39,40].

**Figure 15.** Total cost and carbon dioxide emissions Pareto solution set.

For further analysis, the obtained Pareto solution sets can be classified into I, II, and III, as shown in Figure 16.

**Figure 16.** Classification of total cost and total carbon dioxide emissions (M = 3000) Pareto solution set.

Among them, the total cost of class I is at a low level and carbon emission is at a high level, the total cost and carbon emission of class II are both at a medium level, and the total cost of class III is at a high level and carbon emission is at a low level. For different development stages of recycling system construction and development, different solution sets of different regions can be selected as design solutions under the consideration of total cost and carbon emission only. When the whole society's recycling system has developed to a certain extent, Region III can be chosen.

In this paper, one point from each of the three categories I, II, and III is randomly selected for analysis, and the selected points are shown in Figure 17.

**Figure 17.** Schematic diagram of sample solution selection for the Pareto solution of total cost and total carbon dioxide emissions (M = 3000).

The coordinates of the point selected for Class I are (3.333 × 107, 2.160 × 105), which is noted as the sample solution for Class I, i.e., the total cost is 3.333 × <sup>10</sup><sup>7</sup> CNY and carbon emission is 2.160 × 105 kg. The first, second, third, and fourth parts of the chromosome corresponding to this point are coded as


According to the algorithm design rules, the above code is decoded to obtain the recycling center responsible for each processing center, the model used, and the technology used, as shown in Table 10 (for convenience, Table 11 notes the recycling center to the processing center as Section 1, the processing center back to the recycling center as Section 2, and the processing center to the landfill and paper mill as Section 3, the same as the following table). The location of each node is shown in Figure 18A.


**Table 10.** Class I sample solution analysis.

**Table 11.** Class II sample solution analysis. **Processing Center Recycling Center in Charge Road Section Technical Processing <sup>123</sup>** 372233 432233 594422 614233 7 4, 6 2, 2 1, 1 2, 1 2, 3 10 2 2 2 3 3 11 5 3 2 4 3 12 14 1 1 3 3 13 13 4 2 3 3 14 8, 10, 11, 12, 15 4, 1, 3, 4, 4 2, 1, 2, 2, 2 4, 4, 3, 3, 4 3, 2, 3, 3, 3

**Figure 18.** (**A**) Schematic diagram of the location of each node of the class I sample solution. (**B**) Schematic diagram of the location of each node of the class II sample solution. (**C**) Location diagram of nodes of Class III sample solutions.

The coordinates of the point selected for class II are (3.340 × <sup>10</sup>7, 2.023 × 105), which is noted as the sample solution for class II, i.e., the total cost is 3.340 × 107 *CNY* and the carbon emission is 2.023 × 105 kg. The first, second, third, and forth parts of the chromosome corresponding to this point are coded as


The above codes were decoded to obtain the recycling centers responsible for each processing center, the models used, and the technologies used, as shown in Table 11. The location of each node is shown in Figure 18B.

The coordinates of the point selected for class III are 3.363 × 107, 1.961 × 105 , which is noted as the sample solution for class III, i.e., total cost is 3.363 × 107 *CNY* and carbon emission is 1.961 × 105 kg. The first, second, third, and fourth parts of the chromosome corresponding to this point are coded as,


After the above encoding and decoding, the recycling centers, models, and technologies adopted by each processing center are obtained, as shown in Table 12. The location of each node is shown in Figure 18(C).


**Table 12.** Class III sample solution analysis.

Compare and analyze the sample solutions of class I, II, and III selected above: From the perspective of site selection, it can be seen that among the nodes selected by class I sample solutions, 2 nodes are located in region B, 5 nodes are located in region C, and 1 node is located in region D. Among the nodes selected by region B sample solutions, 1 node is located in region A, 3 nodes are located in region B, 5 nodes are located in region C, and 1 node is located in region D. Among the nodes selected by class III sample solutions, 2 nodes are located in region A, 4 are located in region B, 5 are located in region C, and 1 is located in region D, and the node selected by the sample solution of class III contains class II, and class II contains class I. From the perspective of the selected vehicle types, excluding some invalid codes such as starting and ending points being the same node, it can be seen that the transportation modes between nodes are widely selected, and all four types of vehicles have been applied. From the perspective of the selected technology, it can be seen that the same processing center can choose 2 or more kinds of processing technology, and each processing technology has an application. To sum up, after selecting the regional category, the specific location, vehicle type, and technology should be taken into consideration to optimize the whole system.

#### **7. Conclusions**

Green and low-carbon products are becoming increasingly popular, and green carbon reduction has become the mainstream way of consumption upgrading. This paper analyzes and summarizes the existing courier packaging recycling model, and establishes a new courier packaging recycling model based on the concept of sharing from a low-carbon perspective. From the perspective of engineering research, this paper proposes a complete set of reverse logistics network design process for express packaging, especially from the existing express network, and establishes a network optimization model by combining qualitative and quantitative analysis, which provides a certain technical reference value for similar projects. From the application value point of view, this paper defines the scope of region M. According to the population of each administrative region in region M, the design value is used to weight according to the population number to estimate the courier volume of each administrative region in region M. The location information of 535 courier points in region M was obtained and filtered. The courier packaging recycling mode adopted in region M was determined. This paper also randomly selects one sample solution from each of the three types of solution sets, conducts the coding interpretation of site selection, vehicle selection, and processing technology selection and gives an example design scheme. The express packaging recycling network constructed in this paper can avoid the waste of express packaging, reduce environmental pollution, and promote the sustainable development of social environment and economy for the social development of region M. For the express enterprises in region M, it can improve the utilization rate of express packaging, reduce the cost, actively assume social responsibility, and establish a good corporate image. There are shortcomings in this paper. Affected by the epidemic, there are large errors in the estimation of express business volume in each administrative region of M. The performance of the program written by the relevant algorithm is unstable, and the time complexity and space complexity are not considered, and the algorithm design and program writing should be further optimized. In future research, this aspect should be considered more comprehensively and carefully.

**Author Contributions:** Conceptualization, J.M. and J.C.; methodology, J.C.; software, J.C.; validation, X.L. and H.Z.; formal analysis, J.C. and C.L.; investigation, J.M.; resources, X.L.; data curation, J.C.; writing—original draft preparation, J.C.; writing—review and editing, H.Z.; visualization, C.L. and X.L.; supervision, C.L. and H.Z.; project administration, H.Z. and C.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors thank the editor and the anonymous referees for their helpful comments and critics, and Professor Guangdong Tian for helpful discussions and guidance.

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

#### **Appendix A**

**Table A1.** Total population of each district in Region M.


**Table A2.** The proportion of the population of each district in the total population of Changchun in M region.




**Table A3.** *Cont.*

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
