**1. Introduction**

Scheduling has been certified to be critical in manufacturing for improving the productivity of the manufacturing system and utilization of equipment, as well as shorting the manufacturing cycle [1]. Traditional machining scheduling typically envisions product processing phases where jobs are independent of each other and are not sequence constrained. In the conventional processing scheduling method, the two processes of processing and assembly are separated from each other, which will lead to the destruction of the original parallel relationship between processing and assembly [2]. Generating a processing sequence on the base of the assembly process would reduce the frequency of this phenomenon. Nevertheless, because the process is complicated (such as product hierarchy, number of workpieces, and processing steps), the resolution of these problems considered processing sequence in machining systems is more difficult than traditional machining scheduling. Moreover, it is an NP-hard problem [3].

Temporally, the completion time of workpieces is affected by the hierarchy of the product tree. For example, a product that includes 10 levels is processed from the bottom up. If each level is processed when the previous level is finished, without a doubt, it will prolong the entire production cycle. Therefore, batching the workpieces and making the scheduling scheme in the production of the multi-level structures have practical significance for reducing the completion time.

**Citation:** Li, N.; Feng, C. Research on Machining Workshop Batch Scheduling Incorporating the Completion Time and Non-Processing Energy Consumption Considering Product Structure. *Energies* **2021**, *14*, 6079. https:// doi.org/10.3390/en14196079

Received: 3 August 2021 Accepted: 18 September 2021 Published: 24 September 2021

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

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

1

The manufacturing industry has consumed large amounts of energy in the process of transforming resources into products or services, leading to many environmental problems [4]. The realization of low-carbon manufacturing is extremely important for improving the sustainability of the manufacturing industry [5]. Energy-efficient scheduling, which can be approved in the manufacturing industry, has achieved energy conservation and emissions reduction [6]. Therefore, the scheduling scheme should not only consider the rationality of the process sequence for the workpieces but also reduce the energy consumption of the production system.

The processing energy consumption (PEC) and non-processing energy consumption (NPEC) of the machine are the two components of the production energy consumption of the workshop. PEC stands for the energy consumption of the machine at the processing stage, which is related to the processing power and processing time of the machine. NPEC is the sum of standby energy consumption, conversion energy consumption, and handling energy consumption. Compared with the energy consumed by the machine in other operating phases, the equipment consumes less energy when processing workpieces, especially in mass production, which generally only accounts for approximately 10% of the total energy consumption [7].

Most of the energy consumption in production is generated by auxiliary operations. Here, auxiliary operations are defined as operations that are not directly involved in processing but indispensable in the production process, such as those for equipment standby, state conversion, and workpiece handling. Compared with the energy consumption generated by the processing phase, the energy consumption generated by auxiliary operations can be large. Therefore, if the focus of energy conservation is on developments in processing and energy-saving equipment [8], the considerable energy-saving potential of auxiliary operations will be ignored. In addition, relative to changing a processing technology or researching and developing more energy-saving processing, an optimized workshop scheduling scheme can provide good application value with a low investment [9]. Therefore, in a production system based on processing sequences of workpieces, comprehensively considering the energy consumption composition in the production process, optimizing the allocation of workshop resources, and formulating reasonable scheduling arrangements will be more conducive to reducing energy consumption and improving efficiency in manufacturing enterprises.

Gray wolf optimization [10] (GWO) is a new intelligent optimization algorithm proposed in recent years. Compared with the genetic algorithm (GA) and particle swarm optimization (PSO), GWO algorithm results are more competitive [11]. At present, the gray wolf algorithm has been widely applied in thermodynamics [12], power systems [13], energy and fuels [14], cloud technology [15], and workshop scheduling [16–18]. Lu [16] embedded genetic operators into the multi-objective GWO to enhance the searchability of the algorithm. Qin [17] used the improved multi-objective gray wolf algorithm to solve the casting shop scheduling to minimize the production cycle, total production cost, and total delivery delay. Lu [18] added a random search model based on traditional GWO search to enhance global search capability. Although GWO has been successfully used in many different types of production environments, there is limited literature on GWO to solve energy-saving scheduling problems in a machine-shop, especially to optimize auxiliary production energy consumption. Therefore, we extended the single-objective GWO to the multi-objective GWO to consider completion time and total energy consumption minimization.

Reducing energy consumption through NPCE optimization and minimum completion time are major design goals. The research motivation and research problem will be clearer in the discussion of background research, and then the related mathematical model will be introduced, followed by the MOGWO algorithm and how it is applied to optimize the bi-objective scheduling problem; finally, a presentation of a case study will demonstrate the model and algorithm in the case of two different energy consumption optimization objectives and different algorithms.

In the rest of this paper, the current research progress will be introduced in Section 2 and a bi-objective model based on the research problem is established in Section 3. In Section 4, the working procedure of MOGWO of this optimization problem is described. Section 5 conducts case analysis and comparison. Finally, Section 6 summarizes this article.

#### **2. Literature Review**

The scheduling research of multi-level structure products was first studied for an assembly workshop [19,20]. Li et al. [21] developed four batch strategies to solve three different multi-level product assembly problems. Lu et al. [22] studied tree-like products scheduling problems in an assembly workshop, aiming to minimize the assembly completion time. Wan et al. [23] proposed a visual modeling and scheduling model for assembly processing based on a workflow for the assembly process of complex products and designed heuristic scheduling rules. Suharyanti et al. [24] investigated the optimal lot size of complex products in the job shop. Batching can effectively reduce product production cycles. However, the scheduling research on multi-level structure products is not enough to solve the problem of collaboration production because sequence constraint scheduling should be compatible with the assembly sequence and has higher complexity.

A scheduling problem that comprehensively considers the optimization of the energy consumption with the traditional targets (completion time/production cost) is complex, and it is particularly important to carry out in-depth research on this [25]. At present, a large amount of energy waste occurring in processes that have nothing to do with equipment processing operations has been found. Dahmus and Gutowski [7] analyzed the energy consumption of machining and proved that the actual processing operations only accounted for a small part of the total energy consumption, whereas auxiliary operation energy consumption accounted for 30–50%. The idle time of the machine occupied 16% of the total completion time [26]. Thus, the NPEC is large. Based on the knowledge regarding total energy consumption, it has become a trend to decompose the total energy consumption and reduce energy waste through optimization of the scheduling scheme. Wang [27] simulated a processing process and classification of energy consumption by product quality. In general, the total energy consumption can divide into PEC and NPEC, where NPEC, as indicated above, refers to the energy consumption generated by auxiliary operations such as equipment start-up, shutdown, and idling.

In recent years, more and more research of NPEC has been conducted thoroughly. Luan et al. [28] studied the energy consumption of non-cutting status and established an accurate power model to accurately predict the power of the feed motion. Liu et al. [29] improved the machine utilization rate by 8.2% by optimizing the processing sequence for the workpieces. Peng et al. [30] considered standby energy consumption. Wu et al. [31] studied a renewable energy scheduling problem of a flow shop and established a multi-objective renewable energy power supply model, intending to reduce the processing and idle energy consumption during processing. Gilles et al. [32] investigated the impacts of batch production on energy consumption and order completion time. It was considered that batching could effectively reduce the number of conversions and equipment standby time, thereby reducing the conversion energy consumption and standby energy consumption. Che et al. [25] used a clustering algorithm to determine whether a shutdown operation was required between two tasks to reduce the standby energy consumption and/or optimally sort the processing tasks. Liu et al. [29] integrated scattered short standby periods into a long standby time and judged whether off or assigned other tasks. Wang et al. [33] considered the power changes in the standby state and processing states of the machine. However, the above research mostly focused on the conversion and standby energy consumption separately. In the research of NPEC, they should be considered more systematically.

In the research articles above, most of the influences of workpiece handling on the energy consumption of the workshop were ignored. However, it is more realistic to consider the scheduling and optimization of such handling. The impact of handling on energy consumption has generally not been considered [34–36]. The handling of workpieces

between machines is also part of auxiliary processing, and its energy consumption should belong to the NPEC in the system. Different from previous studies, the NPEC would be determined based on the energy consumption of the equipment standby, conversion, and workpiece handling. In addition, a machine may perform many processing tasks, and the frequent shutdown and start-up of the equipment will increase energy consumption. Therefore, this part of the energy consumption is generally not considered.

Based on the above discussion, the research on multi-level product production and energy-conscious scheduling in the machining workshop is limited. The existing workshop scheduling methods cannot meet the requirements of production and energy saving because each mode has its characteristics. The research on NPEC should be considered when proposing a scheduling method that is different from the traditional optimal energy consumption scheduling. Therefore, the scheduling model of the machining workshop is more complicated than that of the traditional workshop. To bridge this gap, a scheduling method was introduced specifically for machining workshops to minimize completion time and energy consumption.

#### **3. Problem Statement**

#### *3.1. Related Product Structure*

In real-world production, a product is a combination of a group of parts/components with order constraints and can be represented using a tree structure diagram. The processing and assembly of products are conducted according to the tree structure. The edges in the tree represent the assembly constraint relationships, the leaf nodes are the parts/components, and the root node is the final product. Each level's leaf node can be called a child node of its upper-level node, and the processing starts from the lowest-level leaf node. It is composed of *n* parts/components with order constraints. The production of components includes several processes, each process is handled by a machine (in M), and the alternative processing equipment for different processes can be identical.

According to the structural characteristics and commonalities between products, the literature [37] has generally summarized product structures into three types: flat, tall, and complex. The corresponding product structure trees are shown in Figure 1. The flat type is a single-layer product and is directly assembled from first-level parts into products. The tall type has multiple levels of parts/components, and each sub-workpiece contains at most two nodes. The complex type is a multi-layer composite of flat- and tall-type tree structures, in which at least one parent node contains more than two child nodes, as shown in Figure 1. The nodes of each tree are arranged hierarchically, where level 0 represents the complete product, and level 1 is the hierarchical arrangement of parts. For example, the structure of the B product is divided into 1–3 levels, and the branch nodes under it are called components, such as *P*1 and *P*2. The production sequence is as follows: first produce the workpieces *J*3 and *J*4, and then the superior *P*2 and *J*2 for production, and so on.

Based on the coupling relationship(s) between the parts/components in the product structure, batch production and handling are conducted, in which each type of workpiece is divided into equal batches, and the numbers and size of the sub-batch of each product are determined, as well as the sub-batch production and handling sequence. A batch scheduling problem based on the product structure will face difficulties caused by the coordination of the processing times between the workpieces. The arrangement of the production sequence of each sub-batch to meet the coupling sequence of the products is very important. For example, in the mass production of the workpieces in a machining workshop, each type of workpiece can be divided into an equal number of sub-batches. According to the coupling relationship between the workpieces and the production time of the workpiece, the lower-level workpieces are produced first, and then the start processing time of the upper-level workpieces will be later than the next level of workpieces, minimizing machine standby while optimizing handling equipment to reduce NPEC.

**Figure 1.** Three types of product structure.

*3.2. Problem Definition and Assumptions*

The problem can be expressed as follows. Related sets and decision variables are shown in the Appendix A.

The following assumptions are used in the scheduling.

• The time required for loading and unloading workpieces is ignored.

#### *3.3. Problem Formulation*

3.3.1. Objectives

During workpieces processing, there are three main states of equipment: the processing state, conversion state, and standby state. The equipment power and state vary with time [25], as shown in Figure 2. *Sim* and *Cim* represent the start and completion times of the ith operation on device *m*, respectively, and *Omj*(1)*i*(1) represents the *Oji* processing on equipment *m*.

**Figure 2.** Gantt chart of *m* and the distribution of corresponding power.

The total energy consumption of the production system is composed of the PEC (*Ep*) and NPEC (*En*) during the processing of components and parts. Of these, the NPEC consists of the state transition energy consumption (*Es*), the standby energy consumption (*Ew*), and the handing energy consumption (*Ed*) (including the energy consumption of the moving parts/components between machines and moving to the assembly workshop when finished).

•The energy consumption criterion *E*

It consists of the PEC and NPEC. Specifically, this energy consumption metric is given as follows:

$$E = E\_P + E\_n \tag{1}$$

The PEC consumed by all workpieces owing to the processing, as shown in Equation (2), as follows

$$E\_p = \sum\_{j=1}^{m} \sum\_{k=1}^{B\_j} \sum\_{i=1}^{O\_j} \sum\_{m=1}^{w} \alpha\_j T B\_{j\text{kim}} \tag{2}$$

The NPEC consists of the equipment standby energy consumption (*Ew*), conversion energy consumption (*Es*), and handling energy consumption (*Ed*). The equation is given by Equation (3), as follows:

$$E\_n = E\_w + E\_S + E\_d \tag{3}$$

Among them, *Ew* refers to the idling state where the equipment is on non-stop, and no parts/components are processing on the machine. The formula for the energy consumption during the standby period is shown as follows:

$$E\_{w} = \sum\_{k=1}^{B\_{j}} \sum\_{i=1}^{u} \sum\_{m=1}^{w} \left( S B\_{j(k+1)im} - C B\_{j\text{kim}} \right) P\_{m}^{w} \tag{4}$$

In the above, *SBj*(*k*+1)*m* represents the start time of the (*k* + 1)th sub-batch of *j* on the processing equipment *m*, represents the kth sub-batch on equipment *m* completion time, and *Pwm* is the standby power of the machine.

*Es* represents the energy consumption of the equipment state transition. In Equation (5), *αj* − *αj*- is the absolute value of the difference in energy consumption involved in switching power owing to processing different types of workpieces (*j*→*j*-); *Rjjm* is a 0–1 variable. If the processing of the workpiece on equipment *m* is different from the workpieces to be processed, *Rjjm* = 1; otherwise, *Rjjm* = 0.

$$E\_s = \sum\_{j=1}^{n} \sum\_{k=1}^{B\_j} \sum\_{i=1}^{O\_j} \sum\_{m=1}^{M} \mathcal{R}\_{j\text{kim}} \mathcal{R}\_{\vec{j}\vec{l}} \iota\_m \left| a\_j - a\_{\vec{j}} \iota \right| \tag{5}$$

The workshop handling energy consumption is related to the sub-batch quantity Bpj and sub-batch *Qpjk* of *j*, selected handling equipment *Hh*, distance *dmm*- between the equipment, and distance *dmP* between the equipment and assembly workshop P. Ed represents the energy consumption during handling. The transportation of workpieces includes two parts: one part comprises transporting the current sub-batch process to the next processing machine after the completion of the current sub-batch process, and the other comprises transporting it to assembly workshop P after the last process of the sub-batch process is completed. The two parts of energy consumption are described in detail as follows. This can be expressed using Equation (6).

$$E\_d = E\_{jkinum'}^{h} + E\_{jkmP}^{h} \tag{6}$$

After the process *Oji* is processed on equipment m, the workpiece will be transported to the next process *Oj*(*i*+1). The energy consumption *Ehjkimm*- of the transportation equipment *Hh* at the selected equipment m- is shown in Equation (7).

$$\begin{array}{l} E\_{jklmnw'}^{h} = \sum\_{j=1}^{n} \sum\_{k=1}^{B\_{j}} \sum\_{i=1}^{O\_{j}} \sum\_{m=1}^{W} \sum\_{k=1}^{H} S\_{jklh}^{1} H\_{jklmh}^{1} m\_{jk} P\_{j}^{k} t\_{jklmnw'}^{h} \\ = \sum\_{j=1}^{n} \sum\_{k=1}^{B\_{j}} \sum\_{i=1}^{O\_{j}} \sum\_{m=1}^{W} \sum\_{h=1}^{H} S\_{jklh}^{1} H\_{jikmh}^{1} \left[ \frac{Q\_{jkl}}{s\_{jh}} \right] P\_{j}^{h} \frac{d\_{mm'}}{V\_{h}} \end{array} \tag{7}$$

Here, *njk* represents the number of pieces of handling equipment required for the kth sub-batch of *j*, and *thjkimm*- indicates the time that *Hh* moves *j* from process *Oji* of equipment *m* to *m'. Sβjkih* is a 0–1 variable. If *Hh* was selected to handle the ith process of the kth sub-batch of *j*, then *Sβjkih* = 1; otherwise, *Sβjkih* = 0. *Hβjkimh* is also a 0–1 variable. *β* can take two values, 1 and 2; *β* = 1 indicates that the workpieces are transported between equipment; *β* = 2 indicates that the workpieces are transported from the last piece of equipment *m* to assembly shop P. If the ith process of the kth sub-batch of *j* is transported by *Hh*, then *Hβjkimh* = 1; otherwise, *Hβjkimh* = 0.

After the workpieces *j* are processed, they are transported from equipment *m* to assembly workshop P by *Hh*. The energy consumption *EhjkmP* of Hh is given by Equation (8), as follows:

$$E\_{jkmP}^{h} = \sum\_{j=1}^{n} \sum\_{k=1}^{B\_j} \sum\_{h=1}^{H} X\_{jkO\_j} S\_{jkh}^2 H\_{jkmh}^2 \mu\_{jh}^h P\_j^h t\_{jkmP}^h = \sum\_{j=1}^{n} \sum\_{k=1}^{B\_j} \sum\_{h=1}^{H} X\_{jkO\_j} S\_{jkh}^2 H\_{jkmh}^2 \left[ \frac{Q\_{jh}}{S\_{jh}} \right] P\_j^h \frac{d\_{hP}}{V\_h} \tag{8}$$

In the above, *t h jkmP* represents the time taken by *Hh* to transport the kth sub-batch of *j* from *m* where the last process is located to the assembly workshop P. *XjkOj m* is a 0–1 variable. If the last process of the kth sub-batch of *j* is completed on equipment *m*, then *XjkOj m* = 1; otherwise, *XjkOj m* = 0.

• The makespan criterion *C*

The makespan criterion is defined as the maximum time for the last sub-batch for the equipment in the workshop to be processed and transported to the assembly workshop; *CBjBj Ojm* represents the completion time of the last process *Oj* of the last sub-batch *Bj* of the workpiece *j* on the equipment *m*. In addition, *XjkOj mS*2*l jBj OjhH*2*jBj Ojh* represents the completion of the processing of the part *j* after choosing the handling equipment *h* and transporting it to the assembly workshop. The calculation is shown in Equation (9), as follows.

$$\mathcal{C} = \sum\_{j=1}^{u} \sum\_{m=1}^{M} \mathcal{C}B\_{j\overline{B}\_jO\_{j^\*}m} + \sum\_{j=1}^{u} \sum\_{m=1}^{M} \sum\_{h=1}^{H} X\_{j\overline{B}\_jO\_{j^\*}h} S\_{j\overline{B}\_jO\_{j^\*}h}^2 H\_{j\overline{B}\_jO\_{j^\*}h}^2 \frac{d\_{mP}}{V\_h} \tag{9}$$

The multi-objective model is shown in Equations (10) and (11).

$$f\_1 = \min \mathbb{C}\_{\text{max}}\tag{10}$$

$$f\_2 = \min \to \tag{11}$$

3.3.2. Constraints

> Two issues should be considered in the scheduling:

The batches of workpieces should satisfy the condition that the sum of the divided sub-batch batches is equal to the processing quantity of the workpieces; moreover, the number of divided sub-batches should not exceed the total quantity of workpieces.

$$\begin{cases} \ Q\_{\bar{j}} = \sum\_{k=1}^{B\_{\bar{j}}} Q\_{\bar{j}k} \\ \ 2 \le B\_{\bar{j}} \le Q\_{\bar{j}k} \end{cases} \tag{12}$$

The workpieces are split into equal batches. If the number of batches is not an integer, it is rounded down, and the remaining workpieces comprise a single batch.

$$Q\_{jk} = \left\{ \begin{array}{c} \lfloor Q\_j \div B\_j \rfloor \mid k \le B\_{j-1} \\\ \lfloor Q\_j - \lfloor Q\_j \div B\_j \rfloor \end{array} \right. \\ \left. \right. \tag{13} \\ \left. \begin{array}{c} \text{(13)} \\\ \text{(13)} \end{array} \right. \tag{14}$$

If the processing type of the current workpieces is the same as that of the workpieces being processed, there is no need to switch the state of the machine; otherwise, the state needs to be switched.

$$R\_{j\bar{j}'m} = \begin{cases} \ 0, P\_{jtk\dot{t}\dot{t}\dot{m}} = NP\_{j\dot{k}\dot{m}}\\ 1, P\_{jtk\dot{t}\dot{t}m} \neq NP\_{j\dot{k}\dot{m}} \end{cases} \tag{14}$$

The scheduling is performed according to the sub-batch priority relationships of the parts/components in the product structure tree. The process starts processing time at the completion time of the last sub-batch process in the selected equipment, and the sub-batch workpieces after the previous process are completed and transported according to the maximum value between the device moments.

$$\begin{cases} \begin{aligned} \label{eq:1} \mathcal{S}B\_{j k O\_{j}, \mathfrak{m}} = \mathcal{C}B\_{j k o\_{j}, \mathfrak{m}} = 0 \\ \begin{aligned} \label{eq:1} \mathcal{S}B\_{j k \text{im}} = \max\{\mathcal{C}B\_{\mathfrak{m}}^{n-1}, \mathcal{C}B\_{j k (i-1) \mathfrak{m}'} + t\_{j k \text{im} \mathfrak{m}'}^{h} \end{aligned} \end{cases} \end{cases} \tag{15}$$

The production process of the same sub-batch should not be interrupted.

$$CB\_{j\bar{k}im} = SB\_{j\bar{k}im} + TB\_{j\bar{k}im} \tag{16}$$

The part production sequence should meet the requirement that the production time of the lower-level parts/components is earlier than the start-up time of the upperlevel parts/components; that is, the start-up time of the first process in the n-level parts/ components sub-batch should be later than the (*n* + 1) level parts/components.

$$S\_{j\le 1} \ge S\_{j(n+1)1} \tag{17}$$

Among them, *Sjn*1 is the start time of the first process of the *n*-level parts/components, and *Sj*(*n*+1)1 is the start time of the first process of the first sub-batch of the (*n* + 1) parts/components *j*.

For two adjacent processes for the same workpiece, the processing sequence constraints between the processes need to be met, and the next process can only be conducted after the previous process is completed and the workpiece is transported to the selected equipment *m*- for the start of the next process.

$$SB\_{j\bar{k}(i+1)m'} \ge R\_{j\bar{j}'m} R\_{j\bar{k}im} + TB\_{p\bar{j}ko\_{\bar{j}i}} + t\_{j\bar{k}imm'}^{\bar{h}} \tag{18}$$

In each process, *Oji* can only select one piece of machine for processing.

$$\sum\_{m=1}^{w} MP\_{j\text{inr}} = 1\tag{19}$$

One type of handling equipment is selected for each handling instance.

$$\sum\_{h=1}^{H} H\_{jkinh}^{\beta} = 1\tag{20}$$

#### **4. Multi-Objective Gray Wolf Optimization Algorithm**

*4.1. Basic Gray Wolf Optimization Algorithm*

Mirjalili [10] proposed the gray wolf optimizer (GWO) in 2014. The core of the algorithm is to manage an optimization problem by imitating the hunting process of a gray wolf population. Owing to its balance of local and global search capabilities, convergence speed, and depth balancing, it has attracted widespread attention since its proposal.

The basic idea is that *α* wolf was chosen to be the most suitable plan, and *β* wolf and *δ*wolf were the second and third optimal plans. The rest are *ω*. *α* and *β* are the guiders of hunting, followed by *δ* and *ω* wolf. The equation for simulating the corresponding behaviors is defined as follows.

$$
\stackrel{\rightarrow}{D} = \left| \stackrel{\rightarrow}{\mathbf{C}} \cdot \stackrel{\rightarrow}{X}\_p(\mathbf{t}) - \stackrel{\rightarrow}{X}(\mathbf{t}) \right| \tag{21}
$$

Among them, → *D* represents the distance to the prey, *t* indicates the current iteration, → *Xp* represents the prey's position vector, and → *X* represents the wolf's position vector. The coefficient vectors are represented by → *A* and → *C*, and the formula is:

$$
\overrightarrow{A} = 2\overrightarrow{a} \cdot \overrightarrow{r\_1} - \overrightarrow{a} \tag{22}
$$

$$
\stackrel{\rightarrow}{C} = 2\stackrel{\rightarrow}{r\_2} \tag{23}
$$

In the search process, *a* linearly decreases from *2* to *0* and is used to emphasize detection and discovery of prey. → *R*1 and →*r*2 are selected in the range [0,1] randomly. *α*, *β* are the guider in hunting, and *δ* wolves also can join the hunting. The location of the prey (optimal) is unknown. Simulating the hunting behavior of gray wolves, *α, β*, and *δ* wolves are assumed to be more familiar with the potential location of their prey. In each hunt for prey, the three best solutions represented by *α, β* and *δ* wolves will be saved and used in each search, guiding other wolves to the possible position of the prey. The hunting formula is given by Equations (24)–(26).

→

$$
\stackrel{\rightarrow}{D}\_{\mathbf{a}} = \begin{vmatrix} \stackrel{\rightarrow}{\mathbf{C}}\_{1} \stackrel{\rightarrow}{\cdot X}\_{\mathbf{a}} - \stackrel{\rightarrow}{\mathbf{X}} \end{vmatrix}, \stackrel{\rightarrow}{D}\_{\beta} = \begin{vmatrix} \stackrel{\rightarrow}{\mathbf{C}}\_{2} \stackrel{\rightarrow}{\cdot X}\_{\beta} - \stackrel{\rightarrow}{\mathbf{X}} \end{vmatrix}, \stackrel{\rightarrow}{D}\_{\mathbf{a}} = \begin{vmatrix} \stackrel{\rightarrow}{\mathbf{C}}\_{1} \cdot \stackrel{\rightarrow}{\cdot X}\_{\mathbf{a}} - \stackrel{\rightarrow}{\mathbf{X}} \end{vmatrix} \tag{24}
$$

$$
\overrightarrow{X}\_1 = \overrightarrow{X}\_a - \overrightarrow{A}\_1 \cdot \overrightarrow{D}\_a \\
\overrightarrow{X}\_2 = \overrightarrow{X}\_\beta - \overrightarrow{A}\_2 \cdot \overrightarrow{D}\_\beta \\
\overrightarrow{X}\_3 = \overrightarrow{X}\_\delta - \overrightarrow{A}\_3 \cdot \overrightarrow{D}\_\delta \tag{25}
$$

$$
\stackrel{\rightarrow}{X}(t+1) = \frac{\stackrel{\rightarrow}{X}\_1^{\rightarrow} + \stackrel{\rightarrow}{X}\_2^{\rightarrow} + \stackrel{\rightarrow}{X}\_3^{\rightarrow}}{3} \tag{26}
$$

All in all, *GWO* starts with guiding the search process by α, β, and δ. When →*A* > 1,

they diverge and look for prey; otherwise, they find and attack the prey. Finally, if the stopping criterion is met, the optimal solution (i.e., prey) is output. In brief, in each iteration of the algorithm, individuals in the population were divided into *α*, *β*, *δ* and *ω*. The first three belong to the individuals at the decision-making level, representing the historical solution of optimal, suboptimal, and third optimal. ω corresponds to the other individuals. In the algorithm iterations, *α*, *β* and *δ* are locating prey and guiding ω to update its position, completing a sequence of actions including approaching, surrounding, and attacking the prey.

#### *4.2. Application of Multi-Objective Gray Wolf Algorithm*

The multi-objective gray wolf algorithm (MOGWO) [38] added two new components based on the gray wolf algorithm by Mirjalili in 2016. The first component is the archive, which served to store the currently acquired non-dominant Pareto optimal solutions. Then comes the leader selection strategy, which helps decision-makers to choose *α*, *β* and *δ* as the leader of the search process from the archived results. The basic flow chart is shown in Figure 3.

**Figure 3.** Multi-objective gray wolf algorithm flow chart.

#### 4.2.1. Encoding and Decoding Mechanism

Before applying the MOGWO algorithm to a specific problem, we designed an encoding and decoding scheme. It connects the solution space of the problem with the search space under consideration. Hence, designing the correct codec scheme is an important issue that affects the performance of the algorithm.

Code design is based on the types of parts/components and the number of batches. In each chromosome, a gene is represented by three or four numbers. For example, "301" represents the first sub-batch of the third parts/components, and "1003" represents the third sub-batch of the tenth parts/components. In addition, the number sequence in the chromosome represents the processing operations of each sub-batch of parts/components. As shown in Figure 4, each type of part/component needs to go through multiple processing operations. The first occurrence of "101" represents the first processing operation of the first sub-batch of parts/components *J*1, and the second occurrence represents the second processing operation, and so on. Taking an example for illustration, the optional processing

equipment of three types of workpieces and the processing time for each process are shown in Table 1.

**Figure 4.** Chromosome coding, for example.

**Table 1.** Equipment scheduling problem, for example (unit: min).

According to the coding method of the part/component arrangement, it can be assumed that the position of a gray wolf individual in this problem is [1–3], that is, the processing order of the parts/components is 1-3-2. Then, the part/component placement is decoded into a viable scheduling scheme. A Gantt chart corresponding to the first sub-batch of parts/components is shown in Figure 5. Through the decoding process, a suitable machine is selected for each process in each station for processing, and the order of each part and the start time are determined to obtain the objective function value. In Figure 5, initially, all the processes of the first sub-batch of *J*1 are arranged on the machine that can process it earliest, and then the other sub-batches of other parts/components are scheduled. The various processes of the parts/components are arranged on the machine that can complete its processing earliest; if the processing completion time on the allocated machine is less than the earliest processing start time of the scheduled parts/components, it will be arranged before the scheduled parts/components (and so on for the remaining parts/components). Each sub-batch of parts is arranged before the position of each machine; otherwise, it is arranged behind the arranged parts/components.

**Figure 5.** Gantt chart of first sub-batch.

## 4.2.2. Initialization

Initial solutions can affect the obtaining of the optimal solution to a degree. In our study, the problem can be divided into two sub-problems: equipment selection and process sequencing. Therefore, in the initialization phase, the most suitable processing equipment with the lowest processing energy consumption and processing time is selected based on the above encoding and decoding scheme. Then, the sorting plan is obtained according to the processing priority rules of the workpiece and the remaining load at most.

#### 4.2.3. Roulette Selection

In a multi-object search space, comparing solutions is usually not easy; thus, the leader selection mechanism has been designed to solve this problem. The leader provides the α, β, and δ wolves in the least crowded search places. The selection is performed by using the roulette method. The probability of each hypercube is calculated as follows: two individuals are randomly selected from *N* individuals each time, and the individuals with the lower ranking levels are selected first. If the ranking levels are the same, the crowding degree is the first large individual to generate a population of *N*/2 individuals. The process merges the two generated populations into a new progeny population (population size is *N*).

#### 4.2.4. Social Hierarchy

Due to the Pareto advantage of multi-objective planning, the optimal result is usually not a single, which is called a "trade-off solution" in multi-objective planning. According to the Pareto dominance relationship, a population can be divided into several levels. The first-level solution (non-dominant solution or compromise solution) can be expressed as a solution. If there are more than three levels in the whole, β and δ are the second- and third-level solutions. In this study, social stratification was conducted by assuming these three situations:

#### 4.2.5. Update Operator

In this study, individuals are no longer updated according to the decision level, and a hybrid search method combining local search and global search is adopted. The wolf pack generated by each iteration of the algorithm is divided into two parts: the search and tracking operations are carried out, respectively. Then, in the process of searching, the number in each group is dynamically adjusted to achieve the purpose of the individual update.

#### 4.2.6. Sorting Replacement Strategy

Different from the basic GWO algorithm, the newly generated solution is evaluated based on two fitness values: maximum time to completion and total energy consumption. The parent population and progeny individuals produced by global search and local search operations combined to a large new species, then using the non-dominated sorting method and crowded degree to sort the new species. Among them, based on the classification of the solutions, the sorting method reduces the computing complexity, and the crowded degree calculation can save with low levels of similar solution and keep the diversity of the solution space. At the same time, the distribution of individuals on the current Pareto frontier should be as broad and uniform as possible, and the introduction of an elite retention mechanism is conducive to maintaining excellent individuals and improving the overall evolution level of the population.

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

#### *5.1. Data Preparation*

Taking into account the following 10 × 10 production workshop example based on a realistic situation, it comprehensively considers the quantity required for each type of part or component, its handling energy consumption, and its level in the product structure.

For the three different types of products, their complex structures were different, resulting in different processing priorities and handling complexities. As described above, the flat product structure is a single layer, that is, all components have the same priority, and the processing and handling scheduling are relatively simple; in contrast, tall and complex products contain multi-layered parts/components structures. Taking a typical tall product as an example, the specific product structure tree is shown in Figure 6. In the figure, B represents the corresponding product, and the arrow in the figure represents the position of part *Ji*. The part *Pi* in the product structure was divided into different levels; each component or part had a certain demand; they were mixed in batches, and the production was completed. Subsequently, it was transported to the assembly workshop by a transport vehicle for assembly.

Table 2 displays the detail settings of the components, parts, and equipment. The parts and the components are set to 6 and 4, respectively. The rated capacity of the workpiece on the handling equipment is listed in Table 3. In batch scheduling, the production batch had a U-shaped relationship with the production cycle. In general, production batches that are excessively large or small will lead to a longer production cycle. In this study, each workpiece was divided into 2–3 batches, as shown in Table 4. There were three types of handling vehicles, each of which has three available equipment. The power of the handling vehicle was 20 kW, the speed was 30 m/min, and the distance between the assembly workshop and production workshop was 200 m. The information of each type of part and component and the required power is shown in Table 4. In the workshop, the distance between adjacent equipment is 5 m. All cases were simulated in MATLAB R2016b and were tested many times. The algorithm was programmed using Matlab2016b on a personal computer with an Intel(R) Core (TM) i5-930M CPU @ 2.50 GHz. The values of the algorithm parameters were determined by preliminary experiments, and the specific parameters were determined by comprehensive experiments as follows: number of iterations: 250; the number of grids per dimension:15, and population size: 20.

**Figure 6.** Product structure tree of tall.

**Table 2.** Parts/components information (10 × 10).

**Table 2.** *Cont.*

**Table 3.** The rated capacity of the workpiece on the handling equipment.

**Table 4.** Other information about workpieces.

#### *5.2. Result Analysis*

In the experiment, taking into account the total energy consumption and completion time as the objective function and comprehensively considering all the energy consumption involved in the production process, the MOGWO algorithm was used to solve the problem. As shown in Figure 7, the values of the two optimization objectives stabilize at the 189th iteration. Figure 8a shows the population distribution results after 200 iterations. There are eight solutions set in Figure 8b. The point with lower energy consumption and early completion time than others was a selection in Figure 8b. The Gantt chart of this scheme is shown in Figure 9. As the total energy consumption optimization reduces the equipment conversion and waiting time, the completion time is 38,400.86 min, and the total energy consumption is 356,750 kWh. Table 5 displays the handling process of the first sub-batch of each workpiece between the machines and the workshop, where J01 (J = 1, 2, 3, ... ,10) represents the first sub-batch of the *J*th workpiece.

**Figure 7.** Convergence process of two objectives: (**a**) convergence process of C objective; (**b**) convergence process of E objective.

**Figure 8.** Results obtained considering PEC+NPEC: (**a**) final population distribution; (**b**) Pareto frontier distribution.

**Figure 9.** Gantt chart of the optimized schedule scheme.

**Table 5.** One sub-batch handling process for each kind of workpiece.

#### *5.3. Comparison Analysis*

5.3.1. Comparison with the Traditional Model without Considering the NPEC

Different from the experiment above, only PEC is considered in the energy consumption model in this section. As shown in Figure 10, the values of the two optimization objectives stabilize at the 158th iteration. Figure 11a shows the population distribution results after 200 iterations. There are 10 solutions set in Figure 11b. The point with lower energy consumption and earlier completion time than others is selected in Figure 11b. The Gantt chart of this scheme is as shown in Figure 12, the total completion time is 40,900.43 min, and the PEC is 227,000 kWh.

**Figure 10.** Convergence process of two objectives: (**a**) convergence process of C objective; (**b**) convergence process of PEC objective.

**Figure 11.** Results obtained in comparison experiment: (**a**) final population distribution; (**b**) Pareto frontier distribution.

**Figure 12.** Gantt chart of the optimized scheme from the traditional model.

Table 6 displays the energy consumption values and completion times for the two optimization schemes. As could be seen from the table, owing to the reduced standby and conversion time, the total completion time is better than that of the comparative experiment; in terms of energy consumption, the PEC values of the two schemes are very close; the slight difference in processing times may be owing to the different preparation times required on different equipment. The optimization of the handling equipment is not considered in the comparative experiment. The handling energy consumption in the comparative experiment is estimated through the historical handling time, and the energy consumption in our experiment is optimized. As could be seen in Table 6, the optimization of the handling energy consumption can significantly reduce the NPEC, thereby reducing the total energy consumption. Compared to the comparative experiment, the standby energy consumption and handling energy consumption are reduced by 9.95% and 22.28%, respectively. The utilization rates of each machine in the two experiments are shown in Figure 13. It can be seen from the figure that the machine utilization is mostly higher than that in the comparative experiment.

**Table 6.** Comparison of energy consumption and time of two experiments. (unit: kWh).

**Figure 13.** Comparison in Cp/C for machines.

5.3.2. Algorithm Comparison with the NSGA-II

In this section, the most popular multi-objective heuristic algorithm NSGA-II [39] was selected in the experiment for comparative analysis. NSGA-II is an improvement of the NSGA algorithm. It is one of the most outstanding evolutionary multi-objective optimization algorithms so far. The parameters of the algorithm are set as follows: the population size is 100, the maximum number of iterations is set to *G* = 200, crossover probability *Pc* = 0.9, and the mutation probability *Pm* = 0.2. The parameters of MOGWO are set in Section 5.1. They use the same initialization strategy and encoding scheme.

Comparing the performance of multi-objective optimization algorithms is more based on the following criteria: (a) the degree of similarity between the solution set obtained by the operation result and the real Pareto solution set, that is, convergence; (b) the uniformity of the solution set on the Pareto frontier, namely diversity; (c) comprehensively measuring convergence and diversity. According to references [40,41], the following two indicators are used for algorithm performance evaluation: Δ metric and inverted generational distance. The real Pareto frontier of the research in this paper is the set of non-dominated solutions in the final solution set. However, the real Pareto frontier is not known. The final solution set is obtained by the calculation example through multiple independent operations of the algorithm. The specific introduction and calculation formula of these two indicators are given below.

Δ metric: Describe the uniformity of the Pareto front obtained by the algorithm. The calculation method is as follows. The smaller Δ, the more even the solution is, and the better the performance of the algorithm. When Δ is equal to 0, it indicates that the solution obtained by the algorithm is uniformly distributed in the solution set space, generally only appearing under ideal circumstances.

$$\Delta = \frac{d\_f + d\_t + \sum\_{i=1}^{n-1} \left| d\_i - \overline{d} \right|}{d\_f + d\_t + (n-1)\overline{d}} \tag{27}$$

Here, *df* and *dt* are the distances between the boundary point of the Pareto frontier obtained by the algorithm and the actual Pareto frontier boundary point of the problem to be solved; *n* represents the number of solutions in the Pareto frontier obtained by the algorithm operation, and *di* represents the value obtained by the algorithm operation. *d* represents the average value of all *dt*.

Inverted generational distance: A variant of iterative distance not only reflects the convergence of the algorithm but also shows a diversity index, which is a comprehensive evaluation index. The formula of *IGD* is as follows:

$$IGD(S, P^\*) = \frac{\sum\_{\mathbf{x} \in P^\*} dist(\mathbf{x}, S)}{|P^\*|} \tag{28}$$

Among them, *dist(x, S)* represents the individual *x* ∈ *P\** to the nearest Euclidean distance on *S*, and *|P\*|* is the cardinality of the set *P\**. The smaller the value of *IGD*, the more it can approach the entire PF. In addition, when *IGD(S, P\*) = 0*, it means *S* is a subset of *P\*.*

In this paper, two intelligent optimization algorithms are selected for comparison. The results are shown Table 7. The table counts the minimum, average, and standard deviation of the indicators.

**Table 7.** Comparison of two algorithms.

From Table 7, we can draw the following two points:

• According to the spread value ( Δ) in the table, MOGWO is better than the NSGA-II algorithm, and the MOGWO algorithm is more evenly distributed than the solution set obtained by NSGA-II. It is due to the neighborhood search mechanism of the MOGWO, which can increase the probability of obtaining the optimal solution, thereby improving the uniformity of the solution set distribution, and the algorithm has better optimization capabilities.

• According to the inverted generational distance value in the table, the solution obtained by MOGWO has better convergence and distribution than the NSGA-II algorithm. This is because of the unique hierarchical system of the MOGWO algorithm, which can be selected from different dominance levels. The optimal solution to improve the convergence and distribution of the algorithm was chosen.

The Pareto frontiers obtained by the two algorithms are shown in Figure 14. Observation shows that the solution obtained by using the MOGWO algorithm is numerically smaller than the other solution set, that is, it can dominate the solutions obtained by NSGA-II. Additionally, MOGWO proves to be effective in reducing the total energy consumption of scheduling plans. Therefore, the selected algorithm has the best solution effect.

#### *5.4. Sensitivity Analysis*

Sensitivity analysis was used to compare the influence of energy consumption of different auxiliary production operations in this section. If the objective function value and fitness value of different optimization operations change greatly, the sensitivity coefficient and the corresponding auxiliary production operation will be large and sensitive.

Sensitivity analysis is the importance of factor variables of the model to the value of the optimization objective function. The calculation is shown in Equation (29), as follows.

$$S\_{fA} = \frac{\Delta f \, / f}{\Delta A \, / \, A} \tag{29}$$

Here, *SfA* represents the sensitivity of target function value *f* to parameter *A*, and Δ*A*/*A* means the rate of change of a parameter; Δ*f/f* represents the change rate of the target function value caused by the change of factor variable Δ*A*. The sensitivity of total energy consumption to *EW*, *Es*, and *Ed* will be analyzed in our experiment. In this analysis, Δ*A/A* represents the rate of change relative to historical observations when different auxiliary operation optimizations are considered. The change rate of total energy consumption due to different auxiliary operations of optimization is expressed by Δ*f/f*. Table 8 shows the sensitivity coefficients under different optimization schemes; the larger the sensitivity coefficient, the higher the sensitivity of the target to the variation of parameters.

**Table 8.** Sensitivity coefficients under different optimization combinations.

Table 8 compares the changes in total energy consumption when considering different combinations of auxiliary optimization operations. Among them, A represents the scheduling scheme without considering auxiliary operation optimization, where *Ew*, *Es*, and *Ed* are historical forecast values. The standby time, conversion operation, and handling operation are considered separately, and the scheme is expressed by A1, A2, and A3, respectively. Data indicate the three sensitivity coefficients are all positive and the target function value has a certain sensitivity to parameter changes, proving that the optimization model proposed in our research can obtain satisfactory results in reducing energy consumption, but there have been slightly different changes in the target function to the parameter between three schemes. From Table 8, energy consumption is significantly reduced by optimizing standby time, and the target quantity changes significantly to *Ew*; second is handling energy consumption optimization, followed by conversion energy consumption. In this research, due to the large standby time and handling times, the energy consumption optimization effect is relatively significant, especially in the case of large total energy consumption. Through the sensitivity analysis of the above optimization scheme, it is reasonable and effective to comprehensively consider the NPEC in the optimization model.

#### **6. Conclusions and Future Works**

This study researched energy consumption optimization in the production process of a machining workshop, starting from the assembly relationship between the parts/ components of the multi-level product structure. By analyzing the existing energy consumption research, we conducted a systematic study on the NPEC of the workshop. The research was mainly conducted from the following two aspects.

Production scheduling considering energy-saving measures is of great significance for the realization of energy savings and emissions reduction. This study has a set of limitations: for example, we suppose the same power during the standby time and there is no interruption in the processing that limit the versatility of our method. For future works, other possibilities, such as equipment failures and emergency order insertions, will be further integrated into the optimization model.

**Author Contributions:** Conceptualization, N.L. and C.F.; methodology, N.L. and C.F.; software, C.F.; formal analysis, N.L.; investigation, C.F.; resources, N.L.; data curation, N.L. and C.F.; writing— original draft preparation, C.F.; writing—review and editing, N.L. and C.F; visualization, N.L. and

C.F.; supervision, N.L.; funding acquisition, N.L. All authors have read and agreed to the published version of the manuscript.

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

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

**Data Availability Statement:** Data is contained within the article.

**Acknowledgments:** The authors would like to thank the researchers in the Department of Industrial Engineering, China University of Mining and Technology, for their inputs.

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


*2.8. Color Response to Ammonia of Bi-Layer Films*

Each of the bi-layer films was placed into the middle-upper layer of a sealed, homemade acrylic box (500 mL). An aqueous ammonia solution was injected into the bottom of each box with 0.1 mL of different concentrations at 0–200 μM, and the color changes

were determined using a CM2300 spectrophotometer [17]. The digital values were also expressed as color changes (Δ*E*).

#### *2.9. Application in Shrimp Freshness Detection of Bi-Layer Films*

According to the determinations of the film results, the A-CBA and A-CBAL2 films were used as shrimp freshness indicators. An amount of 50 g of fresh shrimp was placed inside a sealed packing box (700 mL), whose inner surface was attached to a film at 4 ◦C for 96 h. The color of each film was obtained using a CM2300 spectrophotometer every 12 h. The total volatile basic nitrogen (TVB-N) value was determined according to the method of Zhang et al. [17].