*Article* **Energy-Efficient Hybrid Flowshop Scheduling with Consistent Sublots Using an Improved Cooperative Coevolutionary Algorithm**

**Chengshuai Li 1, Biao Zhang 1,\*, Yuyan Han 1,\*, Yuting Wang 1, Junqing Li <sup>2</sup> and Kaizhou Gao <sup>3</sup>**


<sup>3</sup> Macau Institute of Systems Engineering, Macau University of Science and Technology, Taipa, Macau 999078, China

**\*** Correspondence: zhangbiao@lcu-cs.com (B.Z.); hanyuyan@lcu-cs.com (Y.H.); Tel.: +86-13863565273 (B.Z.); +86-18864974734 (Y.H.)

**Abstract:** Energy conservation, emission reduction, and green and low carbon are of great significance to sustainable development, and are also the theme of the transformation and upgrading of the manufacturing industry. This paper concentrates on studying the energy-efficient hybrid flowshop scheduling problem with consistent sublots (HFSP\_ECS) with the objective of minimizing the energy consumption. To solve the problem, the HFSP\_ECS is decomposed by the idea of "divide-andconquer", resulting in three coupled subproblems, i.e., lot sequence, machine assignment, and lot split, which can be solved by using a cooperative methodology. Thus, an improved cooperative coevolutionary algorithm (vCCEA) is proposed by integrating the variable neighborhood descent (VND) strategy. In the vCCEA, considering the problem-specific characteristics, a two-layer encoding strategy is designed to represent the essential information, and a novel collaborative model is proposed to realize the interaction between subproblems. In addition, special neighborhood structures are designed for different subproblems, and two kinds of enhanced neighborhood structures are proposed to search for potential promising solutions. A collaborative population restart mechanism is established to ensure the population diversity. The computational results show that vCCEA can coordinate and solve each subproblem of HFSP\_ECS effectively, and outperform the mathematical programming and the other state-of-the-art algorithms.

**Keywords:** hybrid flowshop scheduling; energy efficiency; consistent sublots; collaborative coevolutionary algorithm; variable neighborhood descent

**MSC:** 90B30

#### **1. Introduction**

With the changing climate and environment, green development, energy saving, and emission reduction become the themes of transformation and upgrading of the manufacturing industry. Advanced production scheduling technology can effectively improve production efficiency and reduce energy consumption in the manufacturing industry, enhancing the core competitiveness of enterprises. As a branch of scheduling problems, the hybrid flowshop scheduling problem (HFSP) [1] has a very strong industrial application background, such as microelectronics, furniture, textile, petrochemical, and pharmaceutical fields [2–5]. In HFSP, a group of jobs need to go through a series of processing stages in succession and each stage has multiple identical machines. The goal of the HFSP is to determine the job sequence and machine assignment of these jobs at each stage with considering production constraints. The problem is a very complex combinatorial optimization problem [6], and even on a very small scale, it proves to be NP-hard [1]. In the

**Citation:** Li, C.; Zhang, B.; Han, Y.; Wang, Y.; Li, J.; Gao, K. Energy-Efficient Hybrid Flowshop Scheduling with Consistent Sublots Using an Improved Cooperative Coevolutionary Algorithm. *Mathematics* **2023**, *11*, 77. https:// doi.org/10.3390/math11010077

Academic Editor: Ioannis G. Tsoulos

Received: 18 November 2022 Revised: 16 December 2022 Accepted: 22 December 2022 Published: 25 December 2022

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

most research on the HFSP [7], each job is treated as a whole, and the job cannot proceed to the next stage before the completion of the processing at a given stage [8]. In the actual production scenario [9,10], a job, called a lot in the following, usually consists of a number of identical items. When the lot is large, items already processed completely on a machine need to wait a long time in the output buffer of this machine, whereas the downstream machine may be idle. This scenario will have a negative impact on the production efficiency and lead to unnecessary energy consumption. Therefore, it is very important to develop a scheduling methodology suitable for this scenario to enhance the energy efficiency and core competitiveness of such factories.

In this paper, we introduce the technique of lot streaming into the HFSP, resulting in a novel problem, i.e., lot streaming HFSP. The lot streaming, first introduced by Reiter [11] in the context of job shop scheduling, is preferable for implementing the time-based strategy and widely adopted by top-notch companies to improve their customer service. Lot streaming is the process of splitting a large lot into several sublots and scheduling those sublots in an overlapping fashion to accelerate progress [12]. That is, the lot streaming is used to divide a lot with a large number of items into several sublots with a small number of items. Each sublot can be transported to the downstream stage for processing immediately after its completion at the upstream stage and does not have to wait for the completion of the entire lot. This method can effectively reduce the production cycle and improve production efficiency so that products can be delivered to customers faster, and more orders can be accepted within a limited time. Moreover, this method can effectively increase machine utilization, reduce machine idle time, and thus reduce energy consumption.

According to the lot streaming studies, the lot division methods [13] are equal sublots, consistent sublots, and variable sublots. With equal sublots, a lot can be divided into several sublots with equal size, i.e., each sublot contains the same number of items, and the number and size of sublots remain unchanged throughout the processing process. Consistent sublots mean that a lot is divided into several sublots that may have different sizes, and the number and size of sublots remain unchanged throughout the processing process. Equal sublots can also be understood as a specific case of consistent sublots. Unlike consistent sublots, in variable sublots [14], the number and size of sublots may change throughout the processing process. In real production, variable sublots are rarely used because their diverse nature seriously increases the difficulty of production management. Moreover, its comprehensive cost performance is not high for most enterprises. In contrast, consistent sublots are often used in most enterprises' actual production.

In sum, the energy-efficient HFSP with consistent sublots (HFSP\_ECS) is the focus of our study. To solve the problem, three coupled subproblems must be addressed, i.e., lot sequence, machine assignment, and lot split. Thus, the HFSP\_ECS is much more complex than the classical HFSP, and obviously NP-hard. With its NP-hard property, the metaheuristics are suggested to solve the problem. In addition, when using the metaheuristics, in order to obtain a globally optimal solution, the three subproblems must be coevolved and addressed simultaneously [15,16]. It is therefore natural to employ the cooperative coevolutionary algorithm (CCEA) [17]. Its design is inspired by the natural phenomenon that the coexisting species promote each other and coevolve. The algorithm adopts the strategy of "divide and conquer", which decomposes an optimization problem into several subproblems. In addition, the whole problem is optimized by a reciprocal evolutionary mechanism driven by cooperative or competitive interactions between subproblems [18]. The local search strategy also plays an important role in CCEA. This paper develops an improved cooperative coevolutionary algorithm (vCCEA) by integrating the variable neighborhood descent (VND) strategy [19]. The vCCEA can solve the whole problem by evolving the subproblems simultaneously and interacting between the subproblems. In addition, considering the problem-specific characteristics, a two-layer encoding strategy is designed to represent the solution information and a novel collaborative model is proposed to realize the interaction between subproblems. Special neighborhood structures are designed for different subproblems and two kinds of enhanced local disturbance strategies are proposed to search for potential promising solutions. This algorithm mainly contains four processes, i.e., initialization process, cooperative coevolutionary process, VND processes, and population restart processes. First, an archive that holds several complete solutions is initialized and two populations based on these solutions are built in the initialization process. Then, the two populations and archive coevolve through the collaborative model in the coevolutionary process. While in the cooperative coevolutionary process, the VND process is used to generate a new solution. With the evolution proceeding, the population restart process can be triggered to ensure the population diversity.

The main contributions of this study are as follows. (1) An energy-efficient hybrid flowshop scheduling problem with consistent sublots (HFSP\_ECS) is studied and a mathematical model is developed for it. (2) An improved cooperative coevolutionary algorithm based on the idea of "divide-and-conquer" is proposed by integrating the VND strategy. (3) A novel collaborative model suitable for the specific characteristics of HFSP\_ECS is designed to realize the interaction between the populations and the archive. (4) Two kinds of enhanced local neighborhood structures are proposed to search for potential promising solutions.

The remaining of the paper is organized as follows. A brief literature review is provided in Section 2. Section 3 describes HFSP\_ECS in detail and a linear integer programming model (MILP) is established for a better representation of this problem. Section 4 introduces the algorithm process of vCCEA and improvement strategies in detail. In Section 5, the experimental study design is presented and the results are analyzed. Finally, some conclusions are given and future research prospects are outlined in Section 6.

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

Although HFSP has been studied for several decades, little research has been carried out on energy-efficient HFSP with lot streaming. Most of the existing studies have been conducted with the objective of minimizing the production cycle to optimize HFSP with lot streaming, and little attention has been paid to the energy consumption in the production process. The following is a first review of the HFSP with lot streaming in detail, and then the existing research results on green scheduling are analyzed. Finally, the characteristics of the research problem in this paper are summarized.

With the development of a multi-species small-scale production model in recent years, more and more scholars are focusing on the HFSP with lot streaming. Depending on the number of lots, the lot streaming HFSP can be divided into two main categories, i.e., singlelot HFSP and multiple-lot HFSP. The single-lot HFSP means that only one lot needs to be processed, and how to divide lots and how to sort sublots are two major problems, i.e., sublot size and sublot sequence. Zhang et al. [12] studied a special two-stage HFSP with single-lot that the first stage has multiple identical machines and the second stage only has one machine. They first formulated the problem as an MILP considering the equal sublots, and proposed a heuristic to reach an effective solution. For the same problem, Liu [20] used linear programming and rotation methods to solve the sublot sequence and sublot size, respectively. Moreover, an effective heuristic rule is proposed for the general HFSP with equal sublots. Cheng et al. [21] studied a two-stage HFSP that the first stage only has one machine and the second stage has two parallel machines. Assuming that the number of sublots are known, the closed-form expressions are used to obtain the best sublot sizes. Then, according to the best sublot sizes, the upper bound of the sublot quantities is defined, and an algorithm combining closed-form expressions is used to obtain the global optimal solution. In addition, a heuristic is proposed for the case where the number of sublots is unknown.

Compared with the single-lot HFSP, more research focuses on multiple-lot HFSP. Potts and Baker [22] first showed how to use equal sublots in the one-job model and analyzed equal-sized sublots as a heuristic procedure. After that, they cited some difficulties in multiple-lot scheduling. Kalir and Sarin [23] studied a multiple-lot HFSP with small equal sublots, and proposed a heuristic called bottleneck minimal idleness with the objective of minimizing the maximum completion time. Naderi and Yazdani [24] studied a multiple-lot HFSP with setup time constraints. Assuming that the number of sublots were known, an MILP was established and an imperialist competitive algorithm was proposed. Zhang et al. [25] studied the HFSP with equal sublots, and a discrete fruit fly optimization algorithm was developed for solving this problem, where two main search procedures were designed to balance the exploration and exploitation abilities of the algorithm. For the same problem, Zhang et al. [26] proposed an effective migrating birds optimization algorithm with the objective of minimizing the total flow time, and a heuristic rule was introduced to address the case that the sublots from different lots have the changes to reach the downstream stage at the same time.

The multiple-lot HFSP studied above were all with equal sublots, and this means that sublots from the same lot have the equal size. When the sublots from the same lot are not equal in size, the multiple-lot HFSP is called HFSP with consistent sublots. For example, Ming Cheng and Sarin [13] studied a two-stage HFSP where the first stage only had one machine and the second stage had two identical machines. They used some conclusions from the single-lot scheduling problem, and proposed a mathematical programming-based heuristic method for this problem. Zhang et al. [27] studied a special two-stage HFSP where the first stage had multiple identical machines and the second stage only had a single machine. Additionally, two heuristic strategies were proposed to solve two subproblems, i.e., lot sequence, and lot split. Nejati et al. [28] studied a multiple-lot k-stage HFSP with a specific production scenario. They improved the genetic algorithm and simulated an annealing algorithm for this particular problem, and the effectiveness of the improved strategy was verified. Lalitha et al. [29] studied a special k-stage HFSP where the front k-1 stages only had one machine per stage and the last stage had multiple machines. An MILP was developed and some small-scale problems were solved by the optimizer. A two-stage heuristic algorithm was proposed to solve medium–large scale problems, hierarchically. Zhang et al. [30] studied an HFSP with consistent sublots and considered the setup and transportation operations. A collaboration operator was proposed and a collaborative variable neighborhood descent algorithm was developed based on this operator.

Green development, energy saving, and emission reduction are of great significance to sustainable development. Qin et al. [31] studied an HFSP with an energy-saving criterion, and considered blocking constraints. A mathematical model for HFSP with blocking constraints and energy-efficient criterion was developed and an improved iterative greedy algorithm based on the swap operator was proposed. Duan et al. [32] studied a heterogeneous multi-stage HFSP with energy-efficient for large metal component manufacturing, and an improved NSGA-II combined with the moth-flame optimization algorithm (NSGA– II–MFO) with the objective of minimizing the maximum completion time and carbon emission was proposed. Dong et al. [33] studied a distributed two-stage re-entrant green HFSP, a two-level mathematical model and an improved hybrid slap swarm and NSGA-III algorithm with the objectives of minimum completion time, total carbon emission and total energy consumption was proposed. Geng et al. [34] studied an energy-efficient re-entrant HFSP with considering customer order constraints under Time-of-Use (TOU) electricity price, and a memetic algorithm with an energy saving strategy was proposed to solve this problem.

In summary, both the lot streaming HFSP and the green HFSP have had a certain number of research results. Compared with these studies, the characteristics of our study can be summarized as follows. A k-stage energy-efficient HFSP with consistent sublots is studied in this paper, and the number of machines per stage is not limited. While Ming Cheng and Sarin [13] and Wei Zhang et al. [27] studied the special two-stage HFSP, and Lalitha et al. [29] studied a special k-stage HFSP that the first k-1 stages only have one machine at each stage and the last stage has multiple machines. Compared with these studies, the HFSP\_ECS studied in this paper has a wider scope of application. In study of Ming Cheng and Sarin [13] and Naderi and Yazdani [24], the sublots from different lots can be mixed and cross-processed, but they are prohibited in our research. This is because in real production, the machine needs to be adjusted accordingly before processing different products. Assuming that sublots are allowed to be mixed, the machine will be in a state of frequent adjustment, which has a serious impact on productivity and increases unnecessary energy consumption. Additionally, in the above studies on lot streaming, the research focuses on minimizing the completion time without considering the energy consumption in the process. However, in the actual production, energy consumption is a non-negligible factor. In our study, all machines were turned on and off uniformly, and there was a positive correlation between energy consumption and minimized completion time.

#### **3. Problem Description**

The HFSP\_ECS addressed can be described as follows. A set of lots *J* is to be consecutively processed in a series of *K* stages. Each stage *k* has *Mi* ≥ 1 identical parallel machines and at least one stage has the number of machines greater than one, i.e., *Mi* > 1. Each lot to be processed is made up of a group of identical items. The consistent sublots is employed to split a lot to several sublots with assuming that the maximum sublot quantities are limited. Each sublot contains a certain number of items, and the number of the items contained in a sublot is defined as the sublot size. The number and size of the sublots do not change during the *K* processing stages. At the same stage, different sublots from the same lot are processed continuously on the same machine. Similarly, the items from the same sublot need to be processed continuously. The sublots can proceed to the next stage immediately after its completion of the previous stage. The processing time of a sublot is the product of the sublot size and the item processing time. The processing energy consumption of the sublot is the product of the unit energy consumption and processing time. The idle energy consumption of a machine is the product of unit idle energy and idle time, the idle time, and the idle energy consumption per unit. The scheduling task of the HFSP\_ECS is to solve the three subproblems' lot sequence, machine assignment, and lot split, and its objective is to minimize the energy consumption. The assumptions are summarized as follows:


With the above description and assumptions, to better describe and solve this problem, an MILP [30] is established, the notations and constraints are described as follows:

Objective:

$$Minimize(E\_{\text{max}}) \tag{1}$$

Constraints:

$$\mathbb{C}\_{\text{max}} \ge \mathbb{C}\_{\text{K}, j, \text{L}} \quad \forall j \in \{1, 2, \dots, J\} \tag{2}$$

$$\sum\_{i=1}^{M\_k} D\_{k,j,i} = 1 \quad \forall k \in \{1, 2, \dots, M\_k\}, \forall j \in \{1, 2, \dots, J\} \tag{3}$$

$$\sum\_{x=1}^{L} N\_{j,x} = T\_j \quad \forall j \in \{1, 2, \dots, l\}, \forall c \in \{1, 2, \dots, L\} \tag{4}$$

$$N\_{\circ,\varepsilon} \ge 0 \quad \forall j \in \{1, 2, \ldots, I\}, \forall \varepsilon \in \{1, 2, \ldots, L\} \tag{5}$$

$$N\_{\mathbf{j},\varepsilon} + (1 - \mathcal{W}\_{\mathbf{j},\varepsilon}) \times G \ge 1 \quad \forall \mathbf{j} \in \{1, 2, \dots, I\}, \forall \varepsilon \in \{1, 2, \dots, L\} \tag{6}$$

$$\{N\_{\mathbf{j},\varepsilon} - \mathcal{W}\_{\mathbf{j},\varepsilon} \times G \le 0 \quad \forall \mathbf{j} \in \{1, 2, \dots, I\}, \forall e \in \{1, 2, \dots, L\} \tag{7}$$

$$\mathcal{W}\_{\mathbf{j},\varepsilon} \ge \mathcal{W}\_{\mathbf{j},\varepsilon+1} \quad \forall \mathbf{j} \in \{1, 2, \dots, I\}, \forall \varepsilon \in \{1, 2, \dots, L-1\} \tag{8}$$

$$\mathcal{S}\_{1,j,1} \ge 0 \quad \forall j \in \{1, 2, \dots, I\} \tag{9}$$

$$C\_{k,j,\varepsilon} = S\_{k,j,\varepsilon} + P\_{k,j} \times N\_{j,\varepsilon} \quad \forall k \in \{1, 2, \dots, K\}, \forall j \in \{1, 2, \dots, I\}, \forall e \in \{1, 2, \dots, L\} \tag{10}$$

$$\mathcal{S}\_{k+1,j,\varepsilon} - \mathcal{C}\_{k,j,\varepsilon} \ge 0 \quad \forall k \in \{1, 2, \dots, K - 1\}, \forall j \in \{1, 2, \dots, J\}, \forall \varepsilon \in \{1, 2, \dots, L\} \tag{11}$$

$$\mathcal{S}\_{k,j,\varepsilon+1} - \mathcal{C}\_{k,j,\varepsilon} \ge 0 \quad \forall k \in \{1, 2, \dots, K\}, \forall j \in \{1, 2, \dots, I\}, \forall \varepsilon \in \{1, 2, \dots, L-1\} \tag{12}$$

$$\begin{aligned} \forall \mathbf{y}\_{k,j,1,j} + \mathbf{y}\_{k,j1,j,1} \le D\_{k,j,i} \quad \forall k \in \{1, 2, \dots, K\}, \forall |i = j1, j \in \{1, 2, \dots, l\}, j1 \in \{1, 2, \dots, l\},\\ \forall e \in \{1, 2, \dots, L - 1\}, i \in \{1, 2, \dots, M\_k\} \end{aligned} \tag{13}$$

$$\begin{array}{ll} \Upsilon\_{k,j,1,i} + \Upsilon\_{k,j1,j,i} \le D\_{k,j1,i} \quad \forall k \in \{1, 2, \dots, K\}, \forall j! = j1, j \in \{1, 2, \dots, l\},\\ \jmath1 \in \{1, 2, \dots, l\}, i \in \{1, 2, \dots, M\_k\} \end{array} \tag{14}$$

$$\begin{aligned} \Upsilon\_{k,j,1,i} + \Upsilon\_{k,j1,i} &\geq D\_{k,j,i} + D\_{k,j1,i} - 1 \quad \forall k \in \{1, 2, \dots, K\}, \forall j! = j1, j \in \{1, 2, \dots, l\}, \\ &j1 \in \{1, 2, \dots, l\}, i \in \{1, 2, \dots, M\_k\} \end{aligned} \tag{15}$$

$$\begin{array}{ll} \mathcal{S}\_{k,j,1} - \mathcal{C}\_{k,j1,L} + \mathcal{G} \times (\mathfrak{z} - \mathcal{Y}\_{k,j1,j} - D\_{k,j1} - D\_{k,j1,i}) \ge 0 \quad \forall k \in \{1, 2, \dots, K\}, \forall j! = j1, \\\ j \in \{1, 2, \dots, l\}, j1 \in \{1, 2, \dots, l\}, i \in \{1, 2, \dots, M\_k\} \end{array} \tag{16}$$

$$E\_{\text{process}} = \sum\_{k=1}^{K} \sum\_{j=1}^{J} \sum\_{c=1}^{L} N\_{j,c} \times EP\_{k,j} \tag{17}$$

$$E\_{\rm idle} = \sum\_{k=1}^{K} \sum\_{i=1}^{M\_k} \left( \mathbb{C}\_{\rm max} - \sum\_{j=1}^{I} T\_j \times P\_{k,j} \times D\_{k,j,i} \right) \times EI\_k \tag{18}$$

$$E\_{\text{max}} = E\_{\text{process}} + E\_{\text{idle}} \tag{19}$$

Equation (1) indicates that the optimization objective minimizes the energy consumption *E*max. Equation (2) requires *C*max to be greater than or equal to the completion time of the last sublot of all lots in the last stage. Equation (3) requires that only one processing machine can be selected for each lot at the same stage. Equation (4) indicates that the sum of the items contained in all sublots from the same lot must equal the number of items contained in this lot. Equation (5) defines that the number of items contained in each sublot of lots is greater than or equal to 0. Equations (6) and (7) represent the value of *Wj*,*e*. Equation (8) shows that the nonempty sublot in the lots is expected to precede the empty sublot. Equation (9) make sure the start processing time is a non-negative number. Equation (10) shows the calculation method of completion time. In Equation (11), the sublot is required to complete the processing of the previous stage before starting the next stage of processing. Equation (12) expresses that the sublots from the same lot are processed in numbered order at each stage. Equation (13) defines *Yk*,*j*,*j*1,*<sup>i</sup>* and *Yk*,*j*1,*j*,*i*. They cannot take the value of 1 at the same time. Equations (14) and (15) are dual constraints, similar to Equation (13), emphasize that the machine can only process one lot at a time. Equation (16) indicates that at the same machine, the later lot can be processed only after the previous lot has been processed; otherwise, the equation does not work. Equations (17) and (18) give the calculation method of total process energy and total idle energy, respectively. Equation (19) shows that the total energy consumption is equal to the sum of the total process energy consumption and the total idle energy consumption.

#### **4. Improved vCCEA for Solving HFSP\_ECS Problem**

The proposed vCCEA is developed in this section. First, the design motivation and the algorithm framework are illustrated. After that, the components of vCCEA, involving solution encoding and decoding strategies, initialization strategy, cooperative coevolutionary strategy, VND strategy, and the population restart mechanism, are described in detail. Finally, the whole algorithm procedure is given.

#### *4.1. The Motivations and Framework of vCCEA*

The HFSP\_ECS is a highly complex combinatorial optimization problem. Recall that when solving the HFSP\_ECS, three subproblems must be addressed simultaneously: lot sequence, machine assignment and lot split. These subproblems are not independent but highly coupled. That is, if only one subproblem is optimized, the global optimal solution is almost impossible to be obtained. Therefore, the CCEA is employed, which uses the idea of "divide and conquer" to achieve global optimization by optimizing each subproblem as well as implementing the interaction between subproblems. Among three subproblems, the machine assignment is generally addressed by the proposed heuristic rules [35]. This paper still uses heuristic rules to solve this subproblem, and this rule is incorporated into the decoding strategy. For the other two subproblems, two populations are set up, each of which corresponds to a subproblem. These two populations are lot sequence population and lot split population, which represent lot sequence subproblem and lot split subproblem, respectively. In addition, an archive is also created, where a number of references or complete solutions are stored. It aims to establish collaborative relationships among the individuals from lot sequence and lot split populations in the collaborative coevolutionary process.

The whole algorithm consists of an initialization process, cooperative coevolutionary process of two populations, VND processes, and population restart processes. First, the archive and the two populations are initialized. Then the novel cooperative model was used to control the populations for the cooperative coevolution. In the cooperative coevolutionary process, each individual from the population collaborates with one reference randomly selected from the archive to construct a complete solution. This complete solution is perturbed by the VND process to generate new individuals for updating the population and archive. In this process, different neighborhood structures are designed for individuals from different populations. Moreover, a restart strategy is designed for the individuals who have not been updated for several generations in the population to prevent the algorithm falling into local optima. The vCCEA framework is shown in Figure 1.

#### *4.2. Ending and Decoding*

#### 4.2.1. Solution Encoding

Recall that this problem contains three subproblems: lot sequence, machine assignment, and lot split. Based on the problem specific characteristics, a two-layer encoding strategy is adopted in this paper. The first layer uses a permutation Π*<sup>J</sup>* = - *π*1, *π*<sup>2</sup> ... *π<sup>j</sup>* ... *π<sup>J</sup>* to represent the scheduling order of the lots. Where *π<sup>j</sup>* indicates the lot index and *J* represents the total number of the lots. Note that a legitimate permutation requires that each lot only appears once [36]. The lot that appears in advance in the permutation is given higher processing priority, and the scheduling order of lots in subsequent stages is determined by the heuristic rules mentioned in the solution decoding. The second layer uses a matrix *ZJ*×*<sup>L</sup>* with *J* rows and *L* columns to represent the lot split, where each row represents the segmentation information of a lot. A complete solution consists of two parts: lot permutation and lot split matrix, which can be expressed as Π*J*, *ZJ*×*<sup>L</sup>* .

**Figure 1.** The framework of vCCEA.

Here, a simple illustrative example is given for illustrating the solution encoding. In this example, there are five lots and two stages. The first stage has two parallel machines, i.e., M1 and M2, and the second stage contains three parallel machines, i.e., M3, M4, and M5. The unit idle energy consumption of machines M1 and M2 is 2, and the unit idle energy consumption of machines M3, M4, and M5 also is 2. Lot size, sublot size, item processing time and other specific production data are given in Table 1. According to the above encoding strategy, the encoding for this simple example can be expressed as Π5, *Z*5×<sup>3</sup> , where the first layer is a legal permutation-based encoding vector Π<sup>5</sup> = {3, 5, 1, 4, 2}. This permutation indicates that the current scheduling order is 3,5,1,4,2. In other words, lot 3 was processed first, followed by lot 5, lot 1, lot 4, and lot 2. The matrix *Z*5×<sup>3</sup> serves as the second layer, and is shown in Equation (20). Using lot 5 as an example, the lot 5 is divided into three sublots. The first sublot size is 1, the second sublot size is 1, and the third sublot size is 2.

$$Z\_{5 \times 3} = \begin{bmatrix} 1 & 2 & 2 \\ 2 & 3 & 3 \\ 2 & 2 & 2 \\ 1 & 2 & 2 \\ 1 & 1 & 2 \end{bmatrix} \tag{20}$$


**Table 1.** Illustrative example of HFSP\_ECS.

#### 4.2.2. Solution Decoding

The solution decoding can transform the solution into a feasible schedule, and it mainly solves two problems, lot sequence and machine assignment. For machine assignment, we give priority to idle machines, thus, the "first available machine" rule (FAM) [37] is adopted in this paper. Regarding the lot sequence, the lot scheduling order at stage is determined by the permutation Π*<sup>J</sup>* = - *π*1, *π*<sup>2</sup> ... *π<sup>j</sup>* ... *π<sup>J</sup>* , while the lot sequence of subsequent stages is determined by the "first-come–first-served" rule. That is, the lot completed earlier at the previous stage is given priority to be scheduled at the following stages. In HFSP\_ECS, the sublot of a lot can be immediately transported to the downstream stage for processing when the sublot completes the processing at the current stage. Based on this feature, the "first-come–first-served" rule based on sublot preemption is adopted, i.e., the lot whose first sublot completes the processing at the previous stage first has higher priority at the downstream stage. Under this rule, if the completion time of the first sublots of some lots at the previous stage is equal, the completion time of their second sublots is compared, and so on.

According to the above encoding and decoding strategies, the Gantt chart of the schedule for the illustrative example in Section 4.2.1 is shown in Figure 2. Here, the (a, b) represents a sublot, where a is the lot number, b is the sublot number, and then the minimum energy consumption is 555.

**Figure 2.** The schedule Gantt chart for the illustrative example.

#### *4.3. Algorithm Initialization*

At the beginning of the algorithm, an archive and two collaborative populations need to be initialized. The archive Π*R*[1] , *ZR*[1] , Π*R*[1] , *ZR*[1] ,..., Π*R*[*PS*] , *ZR*[*PS*] is made up of *PS* combinations. Each combination Π*R*[*ind*] , *ZR*[*ind*] , *ind* = 1 ... *PS* represents a complete solution, where Π*R*[*ind*] represents the lot sequence for solution *ind*. Similarly, *ZR*[*ind*] represents the lot split for solution *ind*. When the archive is initialized, the two components of each solution are initialized in two different ways. The lot sequence is

initialized by a random way, while the lot split is determined by the uniform initialization method [38]. The uniform initialization procedure is described as follows.


The two populations are the lot sequence population and the lot split population. The lot sequence population consists of *PS* individuals, i.e., Π[1] , Π[2] ,..., Π[*ps*] . That is, each individual only represents the lot sequence of a solution, this population is initialized with the lot sequences of the solutions in the archive, i.e., Π[*ind*] = Π*R*[*ind*] for *ind* = 1, 2, ... , *PS*. Obviously, the individuals in lot sequence population are indeed not the complete solutions, such that they cannot be evaluated directly. To evaluate each Π[*ind*] , a collaborator must be identified to build an evaluable solution. Here, the lot split individual *ZR*[*ind*] is determined as the collaborator, and the index of this collaborator is recorded by *Col*1[*ind*] = *ind*, where *ind* = 1, 2, ... , *PS*. The individual Π[*ind*] and its collaborator *Z*[*Col*1[*ind*]] construct a complete solution Π[*ind*] , *<sup>Z</sup>*[*Col*1[*ind*]] . The energy consumption value of the solution is also that of the individual Π[*ind*] . Similarly, the lot split population consists of *PS* lot split matrix initially, i.e., *Z*[1] , *Z*[2] ,..., *Z*[*PS*] where *Z*[*ind*] for *ind* = 1, 2, ... , *PS*. The individual <sup>Π</sup>*R*[*ind*] is determined as the collaborator for *Z*[1] , *Z*[2] ,..., *Z*[*PS*] , and the index of this collaborator is recorded by *Col*2[*ind*] = *ind*, where *ind* = 1, 2, ... , *PS*. Individual *Z*[*ind*] in this population and a lot sequence collaborator Π[*Col*2[*ind*]] comprise a new solution Π[*Col*2[*ind*]], *Z*[*ind*] . The energy consumption value of the solution is also that of the individual *Z*[*ind*] .

#### *4.4. Cooperative Coevolution Process*

The whole cooperative coevolutionary process can be divided into two parts: evolution of the lot sequence population and evolution of the lot split population. The evolution of the lot sequence population is first performed. Through this process, individuals in the lot sequence population are updated on the one hand, and certain solutions in the archive can obtain better information of lot sequences on the other hand. Then, the evolution of the lot split population is performed. This process aims to update individuals in the lot split population and ensures that solutions in the archive can obtain better lot split information. Through the above two processes, the evolution of both populations is achieved and the solutions in the archive are also updated during the evolution process. The two processes are described in detail below.

#### 4.4.1. Evolution of the Lot Sequence Population

In the evolution process of the lot sequence population, a complete solution is first constructed by individual Π[*ind*] from lot sequence population and its collaborator *Z* in the archive. To maintain the diversity of the population and to avoid premature convergence, the last collaborator *ZR*[*Col*1[*ind*]] of Π[*ind*] that is pointed by the index is not used. Instead, the lot split matrix *ZR*[*rand*] is randomly selected from the archive as the current collaborator, where *rand* is a randomly generated integer between 1 and *PS*. The combined solution here can be expressed as Π[*ind*] , *ZR*[*rand*] . Then, the VND process is performed on lot sequence Π[*ind*] of the combined solution. A new lot sequence individual Π-[*ind*] is generated by individual Π[*ind*] , and the solution Π[*ind*] , *ZR*[*rand*] comes to Π-[*ind*] , *ZR*[*rand*] . According to the VND characteristics [39], as long as a new individual Π-[*ind*] is generated, the performance of the new solution Π-[*ind*] , *Z*[*rand*] needs to be evaluated first. If the objective value of Π-[*ind*] , *Z*[*rand*] is better than Π[*ind*] , *Z*[*rand*] , then the archive and the population will be updated by a new solution Π-[*ind*] , *Z*[*rand*] . Otherwise, the VND process continues. If the whole VND process fails to find a good Π-[*ind*] , then the evolution of the individual is ended for this time.

The process of updating the archive and lot sequence population can be described as follows. If the objective value of the new solution Π-[*ind*] , *ZR*[*rand*] is better than Π[*ind*] , *<sup>Z</sup>R*[*Col*1[*ind*]] , the solution Π[*ind*] , *<sup>Z</sup>R*[*Col*1[*ind*]] will be updated by the new solution Π-[*ind*] , *ZR*[*rand*] , i.e., Π[*ind*] = Π-[*ind*] , *Col*1[*ind*] = *rand*. Note that the last collaborator *ZR*[*Col*1[*ind*]] may have been changed in the evolutionary process of the lot split population. At the same time, the archive is attempted to be updated. If the objective value of the new solution Π-[*ind*] , *ZR*[*rand*] is better than Π*R*[*rand*] , *ZR*[*rand*] , then the individual Π*R*[*rand*] will be updated by individual Π-[*ind*] , i.e., Π*R*[*rand*] = Π-[*ind*] . The above process is repeated from *ind* = 1 to *ind* = *PS*. Given the above, the coevolutionary process for the individuals Π[*ind*] in the lot sequence population is shown in Figure 3.

**Figure 3.** Evolution process of the lot sequence population.

An efficient neighborhood structure plays a key role in the whole algorithm. To design a good neighborhood structure, the problem characteristics must be exploited. When the lot sequence population evolves, the neighborhood structure only works on the lot sequence. That is, a new individual Π-[*ind*] is formed by perturbing the individual Π[*ind*] in the neighborhood structure. Therefore, to obtain a better lot sequence, three neighborhood structures are specially designed based on a solution encoding strategy, and the VND process is used to switch the neighborhoods. Two of these neighborhood structures are the insert and swap operations that are widely used in HFSP problems, referred to as lot insertion and lot swap. However, when solving the large-scale problems, a single insertion or swap may not effectively perturb the current solution. Therefore, a lot swap operation with a large search range is proposed by improving the lot swap, called the Enhanced lot swap in this instance. The lot insertion is not enhanced here since its high time complexity. The three neighborhood structures for lot sequence are described below: (1) Lot insertion. A new lot sequence is formed by taking a random lot from the lot sequence and inserting it into a randomly different position. (2) Lot swap. Take two lots at random from the lot sequence and exchange their positions in the sequence. (3) Enhanced lot swap. Perform *l* times of the lot swap, where *l* is dynamically determined by the number of lots in the lot sequence. We set *l* as *L* × *J*, where *L* is a real number between 0 and 1, and *J* is the number of lots. Additionally, the detailed process of lot sequence population evolution is shown in Algorithm 1.


1: Define a set of neighborhood structures *N*<sup>1</sup> *<sup>k</sup>*, *k* = 1, . . . , *k*max 2: for *ind* = 1 to *PS* 3: *rand* ← generate a random integer in [1 − *PS*] 4: Π[*ind*] , *ZR*[*rand*] ← constitute a complete solution 5: Define <sup>Π</sup> <sup>←</sup> <sup>Π</sup>[*ind*] 6: Let *k* ← 1 , *Count* ← 0 7: while *k* ≤ *k*max do 8: while *Count* < *C* do 9: Π-[*ind*] <sup>←</sup> *Neighborhood*(Π, *<sup>N</sup>*<sup>1</sup> *k*) 10: if Π-[*ind*] , *ZR*[*rand*] better than Π, *ZR*[*rand*] 11: *Count* ← 0 , *k* ← 1 , Π ← Π-[*ind*] 12: if Π-[*ind*] , *ZR*[*rand*] better than Π[*ind*] , *<sup>Z</sup>R*[*Col*1[*ind*]] or *ZR*[*Col*1[*ind*]] was changed 13: <sup>Π</sup>[*ind*] <sup>←</sup> <sup>Π</sup>-[*ind*] , *Col*1[*ind*] = *rand* 14: end if 15: if Π-[*ind*] , *ZR*[*rand*] better than Π[*rand*] , *ZR*[*rand*] 16: <sup>Π</sup>[*rand*] <sup>←</sup> <sup>Π</sup>-[*ind*] 17: end if 18: else 19: *Count* + + 20: end if 21: end while 22: *k* + + 23: end while 24: end for

#### 4.4.2. Evolution of the Lot Split Population

Similar to the evolution process of the lot sequence population, a complete solution is constructed by individual *Z*[*ind*] from lot split population and its collaborator Π*R*[*rand*] in the archive, where *rand* is a randomly generated integer in the range [1, *PS*]. The constructed solution here can be expressed as Π*R*[*rand*] , *Z*[*ind*] . After that, the VND procedure is executed on the solution Π*R*[*rand*] , *Z*[*ind*] , and the neighborhood structure here only works on the lot split. Through the VND process, *Z*[*ind*] becomes *Z*-[*ind*] , and thus, the solution Π*R*[*rand*] , *Z*[*ind*] comes to Π*R*[*rand*] , *Z*-[*ind*] . In the process, as long as a new solution Π*R*[*rand*] , *z*-[*ind*] is generated, the new solution Π*R*[*rand*] , *Z*-[*ind*] is evaluated. If Π*R*[*rand*] , *Z*-[*ind*] is better than Π*R*[*rand*] , *Z*[*ind*] , then the solution Π*R*[*rand*] , *Z*-[*ind*] 

is used to update the archive and lot split population. Otherwise, the VND process continues to find a potentially better solution. The process of updating the archive and lot split population can be described as follows: If the solution Π*R*[*rand*] , *Z*-[*ind*] is better than Π*R*[*Col*2[*ind*]], *Z*[*ind*] , then the *Z*[*ind*] and the *Col*2[*ind*] will be updated by *Z*-[*ind*] and Π*R*[*rand*] , i.e., *Z*[*ind*] = *Z*-[*ind*] , *Col*2[*ind*] = *rand*. It should also be noted that Π*R*[*Col*2[*ind*]] may have been changed as the lot sequence population evolves. At the same time, if the objective value of the solution Π*R*[*rand*] , *Z*-[*ind*] is better than Π*R*[*rand*] , *ZR*[*rand*] , set *ZR*[*rand*] = *Z*-[*ind*] . This process is shown in Figure 4.

**Figure 4.** Evolution process of the lot split population.

In the evolution process of the lot split population, to obtain high-quality individuals *Z*-[*ind*] , three neighborhood structures acting only on the lot split part are specially designed, and the VND process is used to switch the neighborhood. The three perturbation strategies for lot split are described below: (1) Lot split mutation. As shown in Figure 5, from the lot split matrix, a lot with two or more sublots are randomly selected. Reduce a random number in distribution U [1,5] from the size of one sublot and add the number to the size of another sublot. (2) Enhanced lot split mutation. Perform *l* times lot split mutation, where *l* is dynamically determined by the number of lots. We set *l* as *L* × *J*, where *L* is a real number between 0 and 1, and *J* is the number of lots. (3) Stochastic splits. All lots are redivided into sublots in a random manner. The procedure of the evolution of the lot split population is given in Algorithm 2.

#### **Algorithm 2** Evolution of the lot split population

1: Define a set of neighborhood structures *N*<sup>2</sup> *<sup>k</sup>*, *k* = 1, . . . , *k*max 2: for *ind* = 1 to *PS* 3: *rand* ← generate a random integer in [1 − *PS*] 4: Π*R*[*rand*] , *Z*[*ind*] ← constitute a complete solution 5: Define *<sup>Z</sup>* <sup>←</sup> *<sup>Z</sup>*[*ind*] 6: Let *k* ← 1 , *Count* ← 0 7: while *k* ≤ *k*max do 8: while *Count* < *C* do 9: *Z*-[*ind*] <sup>←</sup> *Neighborhood*(*Z*, *<sup>N</sup>*<sup>2</sup> *k*) 10: if Π*R*[*rand*] , *Z*-[*ind*] better than Π*R*[*rand*] , *Z* 11: *Count* ← 0 , *k* ← 1 , *Z* ← *Z*-[*ind*] 12: if Π*R*[*rand*] , *Z*-[*ind*] better than Π*R*[*Col*2[*ind*]], *Z*[*ind*] or Π*R*[*Col*2[*ind*]] was changed 13: *<sup>Z</sup>*[*ind*] <sup>←</sup> *<sup>Z</sup>*-[*ind*] , *Col*2[*ind*] = *rand* 14: end if 15: if Π*R*[*rand*] , *Z*-[*ind*] better than Π*R*[*rand*] , *ZR*[*rand*] 16: *<sup>Z</sup>*[*rand*] <sup>←</sup> *<sup>Z</sup>*-[*ind*] 17: end if 18: else 19: *Count* + + 20: end if 21: end while 22: *k* + + 23: end while 24: end for

**Figure 5.** Illustrations of the lot split mutation.

#### *4.5. Coevolutionary Population Restart*

With the evolving of the algorithm, the diversity of the two populations might be reduced. In this case, the efficiency of the population coevolution may be poor. To avoid the algorithm falling into the local optimality, two different restart strategies are adopted for the populations. For the lot sequence population, an individual Π[*ind*] is reinitialized if it has not been improved in a predetermined number of *R* consecutive generations. The novel individual should contain valuable information about the original individual and remain somewhat different from the original individual. For this purpose, a two-point crossover (TPX) method was used, as illustrated in Figure 6. Where two parent lot sequences are randomly selected from the archive because good solutions are stored in an archive. For the lot split population, an individual *Z*[*ind*] is also reinitialized if it has not been updated in a predetermined number of *R* consecutive generations. Due to the lot split matrix is different from the regular sequence, the classical TPX might produce infeasible schedules. Therefore, a cooperative selection operator is proposed. When determining the split information for one lot, two solutions are selected at random from the archive and compared based on their objective values, and the split information for this lot comes from the better one. The

process is repeated from the first lot to the last lot. Here, we use *Z*[*ind*] *<sup>j</sup>* to represent the lot split information for lot *j* in individual *ind*, and this process is shown in Algorithm 3.

> Π*b*, *Z<sup>b</sup>*


*j*

7: end if

6HTXHQFH

**Figure 6.** Illustration of TPX for a lot sequence.

#### *4.6. The Algorithm Procedure*

With the above description, the whole vCCEA is displayed in Algorithm 4. Where the *UpdateBestSolution*(Π, *Z* ) means update the optimal solution using solution Π, *Z* , and the *Age*(Π, *Z* ) represents the number of consecutive update failures of the solution composed of individual Π (or *Z*) and their collaborators.

**Algorithm 4** Lot split population restart


```
19: Π[ind] ← Π-
                              [ind]
                                 , Col1[ind] = rand
20: AgeΠ[ind]
                              , ZR[Col1[ind]]  ← 0
21: else
22: AgeΠ[ind]
                              , ZR[Col1[ind]]  + +
23: end if
24: if 
                   Π-
                    [ind]
                       , ZR[rand]

                              better than 
                                      Π[rand]
                                           , ZR[rand]

25: Π[rand] ← Π-
                               [ind]
26: end if
27: else
28: Count + +
29: end if
30: end while
31: k + +
32: end while
33: end for
34: for ind = 1 to PS
35: rand ← generate a random integer in [1 − PS]
36: 
        ΠR[rand]
              , Z[ind]

                   ← constitute a complete solution
37: Define Z ← Z[ind]
38: Let k ← 1 , Count ← 0
39: while k ≤ kmax do
40: while Count < C do
41: Z-
             [ind] ← Neighborhood(Z, N2
                               k)
42: if 
              ΠR[rand]
                   , Z-
                     [ind]

                         better than 
                                 ΠR[rand]
                                       , Z

43: UpdateBestSolutionΠR[rand]
                                     , Z-
                                       [ind]

44: Count ← 0 , k ← 1 , Z ← Z-
                                   [ind]
45: if 
                   ΠR[rand]
                        , Z-
                          [ind]

                              better than 
                                      ΠR[Col2[ind]], Z[ind]

46: Z[ind] ← Z-
                              [ind]
                                , Col2[ind] = rand
47: AgeΠR[Col2[ind]], Z[ind]
                                        ← 0
48: else
49: AgeΠR[Col2[ind]], Z[ind]
                                        + +
50: end if
51: if 
                   ΠR[rand]
                        , Z-
                          [ind]

                              better than 
                                      ΠR[rand]
                                            , ZR[rand]

52: Z[rand] ← Z-
                              [ind]
53: end if
54: else
55: Count + +
56: end if
57: end while
58: k + +
59: end while
60: end for
61: for ind = 1 to PS
62: if AgeΠ[ind]
                  , ZR[Col1[ind]]  > R
63: Restart1(Π[ind]
                    ), AgeΠR[Col2[ind]], Z[ind]
                                        ← 0
64: end if
65: if AgeΠR[Col2[ind]], Z[ind]
                            > R
66: Restart2(Z[ind]
                    ), AgeΠR[Col2[ind]], Z[ind]
                                        ← 0
67: end if
68: end for
69: end while
70: Output the best solution 
                   Πbest, Zbest 
                           .
```
#### **5. Experimental Analyses**

In this section, the performance of the proposed vCCEA is evaluated by experimental design and results analysis. The simulation experiment environment of this paper is a PC with 3.60 GHz Intel Core i7 processor and 32 GB RAM. The vCCEA and all compared algorithms are written in the Visual Studio 2019 C++, and run on the release x64 platform. In the algorithm test, the maximum running time is used as the algorithm termination to ensure fairness. In addition, it is considered that the algorithm has practical significance only when it can solve the problem in an acceptable time. Therefore, the termination condition is set as *t* × *J* × *K* milliseconds, where *J* indicates the number of lots and *K* represents the number of stages, respectively. Referring to the literature [30], *t* is set as 80.

#### *5.1. Experimental Dataset and Performance Indicators*

In this paper, two benchmark sets *β*<sup>1</sup> and *β*<sup>2</sup> are designed to verify the validity of the vCCEA. Where 48 small-scale instances are designed in *β*<sup>1</sup> to study the difference between the MILP and the vCCEA in solving HFSP\_ECS, and 100 medium–large scale instances solved by the metaheuristic algorithm are designed in *β*<sup>2</sup> to verify the performance of vCCEA. In *β*1, the number of lots is *J* comes from {6, 8, 10, 12, 14}, and the number of stages *S* comes from {3, 5, 8}. Thus, there are 15 different combinations of *J* × *S* that can be obtained. In *β*2, the number of lots in *J* comes from {20, 40, 60, 80, 100}, and the number of stages is *S* comes from {3, 5, 8, 10}. Similarly, there are 20 different combinations in *β*2. For *β*1, only one instance is randomly generated per combination. For *β*2, five instances are randomly generated per combination. Thus, there are 15 small scale instances and 100 medium–large scale instances in *β*<sup>1</sup> and *β*2, respectively. In *β*<sup>1</sup> and *β*2, the number of parallel machines at each stage is randomly generated from the range [1, 5]. In addition, the number of items for each lot is obtained from a uniform distribution *U*[50, 100], the processing time of items at each stage is randomly sampled from the uniform distribution *U*[1, 10], the processing energy consumption per unit time is obtained from the range [2, 5], the energy consumption per idle unit of the machine takes a value in the range [1, 3], and the maximum number of sublots is set as 5. In this study, time is measured in seconds, and energy consumption is measured in joules, and the relative percentage increase (RPI) is used as the performance metric. The RPI is calculated as in Equation (21).

$$RPI = \frac{E\_{avg} - E\_{best}}{E\_{best}} \times 100\tag{21}$$

where *Eavg* is the average energy consumption of an instance solved by the given algorithm independently performed several times, and *Ebest* is the best result obtained by all the compared algorithms. Algorithms with smaller RPI values will have better performance.

#### *5.2. Parameter Setting*

Appropriate parameters are very important to the metaheuristics, which can effectively improve their efficiency and robustness. There are four parameters in the vCCEA proposed in this paper, including the number of solutions in the archive (*PS*), the maximum number of consecutive failures in a neighborhood during the VND process (*C*), the parameters (*L*) that control the number of executions in the two enhanced neighborhood structures and the maximum number of successive generations (*R*) of updating the individual in two populations unsuccessfully. We first determine the value level of each parameter through preliminary experiments, where the details of value levels are shown in Table 2. To verify the influence of each parameter and its value at different levels, an orthogonal array *L*16 is designed using the Taguchi experimental method to determine their combinations and is displayed in Table 3. For the combinations in Table 3, five instances with different scale problems are selected from benchmark set *β*2, and the five different problem scales are 20 × 5, 40 × 5, 60 × 5, 80 × 5, and 100 × 5. Each instance is run independently 20 times, and the RPI value of each instance is calculated. Then, as shown in Table 3, the average RPI value of the five instances in each combination is collected as the response value.


**Table 2.** Parameter level factor.

**Table 3.** Orthogonal array and response values.


The trend of the parameter level is shown in Figure 7, and Table 4 gives the significance rank of each parameter of the vCCEA. According to Figure 7 and Table 4, it can be concluded that parameter *L* has the greatest impact on the algorithm among these parameters. This is because the parameter *L* is related to the VND strategy. In the process of VND, the good or bad neighborhood structure has an important influence on the algorithm. A good neighborhood structure can promote the exploration ability of the algorithm and speed up the convergence. For parameter *PS*, a larger population size can accommodate more potential solutions and help the algorithm search globally. However, it does not support longitudinal and deep search in a limited running time. Too small a population size is not conducive to global search. For parameter *C*, a too small value will not make full use of each neighborhood perturbation strategy and a too large value will waste the computational time. For parameter *R*, if the value is set too small, the good information of advanced individuals in the two populations cannot be fully utilized, and if the value is set too large, the diversity of the population cannot be guaranteed and the algorithm may converge prematurely. Therefore, the appropriate parameters are critical for vCCEA. Through the above parameter experiment and analysis, the best parameter combination we can obtain is *PS* = 10, *C* = 15, *L* = 0.3 and *R* = 150. This parameter combination is used in the following experiments.

#### *5.3. Evaluation of the Algorithm Components and Strategies*

In this subsection, the algorithm components and strategies are validated and analyzed for effectiveness. Our algorithm contains the VND process, collaborative model, and two enhanced neighborhood structures, and these three strategies are not independent but highly coupled. To verify the effect of each of the three components and the cooperation between them, three other versions of vCCEA are constructed, vCCEA\_1, vCCEA\_2, and vCCEA\_3, where the vCCEA\_1 is the vCCEA that removes the VND process. The vCCEA2 is used to verify the validity of the collaboration model. The VCCEA\_3 is the vCCEA that remove the two enhanced neighborhood structures during the VND process. And the *β*<sup>2</sup> is used to verify vCCEA and the other three versions of vCCEA. For each instance in *β*2, the four algorithms are independently run 20 times. For each algorithm, the RPI value of each instance is obtained first. Then, the average RPI of five instances from the same scale problem is calculated and represented by the average RPI values (ARPI). The results are shown in Table 5.

**Figure 7.** The trend of the parameter level.



From Table 5, we can clearly see that the whole vCCEA is the best performer. In the other three versions of the vCCEA, the vCCEA\_1 is the worst one. In the cooperative coevolutionary process, the VND process is used to generate new individuals Π(or *Z*). With the same perturbation strategy, the result of neighborhood switching using VND process is obviously better than that of traditional neighborhood perturbation. Thus, the VND process is crucial to the algorithm. Among these 20 problems with different sizes, the vCCEA\_2 is worse than the vCCEA. It can be seen that the collaborative model proposed in this paper is effective. By comparing vCCEA\_3 and vCCEA, the validity of the two enhanced neighborhood structures is proven. These two enhanced neighborhood structures can enlarge the search area of the VND process, and it is beneficial for the vCCEA to find the potential promising solution in a larger solution space.


**Table 5.** Comparison of vCCEA components.

#### *5.4. Evaluation of the vCCEA on the Small-Scale Instances*

This section focuses on the differences between vCCEA and MILP when solving the small-scale problems. We use the Gurobipy 9.1.2 optimizer to run the MILP on the instances in *β*1, and the maximum running time is limited to 3600 s. In addition, the vCCEA is used to solve the instances in *β*1, and the termination condition is set to 80 × *J* × *K*. The results are shown in Table 6.


**Table 6.** The validation results of MILP.

From Table 6, it can be concluded that for small-scale instances, both the MILP and vCCEA can find optimal solutions. The MILP and vCCEA find the same results for the first 13 instances in *β*1. In addition, for instances 14\_5 and 14\_8, the vCCEA revealed better results. As the complexity of the problem increases, the effectiveness of the MILP is gradually inferior to the vCCEA. For large-scale instances, the MILP model has difficulty in providing a good solution in a short time. As we know that time is a non-negligible factor in actual production, thus the near-optimal solutions are required in an acceptable time. Therefore, with the increasing size of the instances, the advantages of vCCEA become more and more obvious.

#### *5.5. Evaluation of vCCEA on the Medium–Large Scale Instances*

Next, the performance of the proposed vCCEA is evaluated on the medium–large scale instances in *β*2. Here, we collected five metaheuristic algorithms for comparisons, namely, CVND [30], GA [40], GAR [24], VMBO [41], and DABC [42], which are those presented for the HFSP in the literature most recently and have been proven to have excellent performance. For the HFSP\_ECS, three highly coupled subproblems need to be solved, i.e., lot sequence, machine assignment, and lot split. Due to the specificity of the problem, we retain the original characteristics of each comparison algorithm and modify it to adapt to our problem. All algorithms use the same double-layer encoding and decoding strategies as proposed in this paper, and select the corresponding lot split operator from this article. As these comparison algorithms have been partially changed for adapting to our problem, their parameters are also optimized and adjusted on the original basis by using the DOE method to ensure that these algorithms can play with better performance. For each instance in *β*2, each algorithm is run independently 20 times, and the average energy consumption and the ARPI of five instances from the same scale problem are calculated. The experimental results are given in Table 7. Additionally, to more visually demonstrate the differences between these six algorithms, the means and 95% least significant difference (LSD) confidence intervals [43] were analyzed. Figure 8 shows the confidence interval comparisons between vCCEA and each algorithm, and Figure 9 shows the confidence interval comparison among all algorithms, where the *X*-axis represents the various algorithms and the *Y*-axis is the ARPI value.

**Table 7.** Comparison results of vCCEA and other algorithms on *β*2.


**Figure 8.** Confidence interval graph. (**a**) Confidence intervals of vCCEA and CVND; (**b**) confidence intervals of vCCEA and GA; (**c**) confidence intervals of vCCEA and GAR; (**d**) confidence intervals of vCCEA and VMBO; (**e**) confidence intervals of vCCEA and DABC.

**Figure 9.** Confidence interval graph for these six algorithms.

As seen in Table 7, the vCCEA is the best one among these algorithms, which obtains best results for 19 of the 20 different problems in *β*2. The DABC algorithm finds the optimal solution for the remaining one scale problem, with the scale being 20 × 8. The last row of Table 7 gives the average energy consumption and average RPI for all instances. It is obvious that vCCEA has the best results among all the algorithms. In addition, according to the confidence intervals shown in Figures 8 and 9, it can clearly be seen that the performance of the proposed vCCEA is obviously better than that of the other five algorithms. To further evaluate the performance of the algorithm, we analyze the convergence of the algorithm, and the convergence curves of these algorithms on two examples are given. These two examples are from 40 × 5 and 80 × 10, respectively, and the convergence curves are shown in Figures 10 and 11. The *X*-axis represents the running time of the algorithm, and the *Y*-axis represents the energy consumption value.

**Figure 10.** The convergence curve for instances of 40 × 5.

**Figure 11.** The convergence curve for instances of 80 × 10.

From Figures 10 and 11, the convergence speed of vCCEA is the fastest, and the convergence degree is also better than that of the other algorithms. This is closely related to the cooperative coevolutionary strategy and VND process in the proposed vCCEA. Two populations in the algorithm evolve separately using the VND process specifically designed for them and constantly interacting with the archive. As well as with the population restart strategy, the local search capability of the vCCEA can be ensured, and it also balances the diversity and density of the populations, which further improves the algorithm performance.

Therefore, through the above analysis, the conclusion can be drawn that the vCCEA can effectively solve HFSP\_ECS and the robustness of the vCCEA can be guaranteed.

#### **6. Conclusions**

In this paper, the energy-efficient flowshop scheduling problem with consistent sublots (HFSP\_ECS) is studied, and it supports the overlap of successive operations within a multistage manufacturing system. This is a highly complex combinatorial optimization problem that consists of three highly coupled subproblems. We use the minimized energy consumption as the optimization objective. By limiting the maximum number of sublots, a linear integer programming model (MILP) of the addressed problem is established and its validity is verified by the Gurobi optimizer. An improved cooperative coevolutionary algorithm (vCCEA) is proposed by integrating the variable neighborhood decent (VND) strategy. In vCCEA, with the consideration of the problem-specific characteristics, a two-layer encoding strategy is designed, and a novel collaborative interaction model is proposed. Additionally, to ensure the local search ability of the algorithm, different neighborhood structures are designed for different subproblems, and two kinds of enhanced local neighborhood structures are proposed to search for potential promising solutions. To avoid trapping into the local optima, a population restart mechanism is designed. Moreover, through a large number of experiments on different benchmark sets, the effectiveness of the proposed strategies is proved. The experimental results show that vCCEA is significantly better than the mathematical programming and the other algorithms in solving the HFSP\_ECS.

For the HFSP\_ECS, the maximum sublot quantities is limited in this paper, so in the future, how to divide the lots and the number of sublots is a direction of our research. At the same time, we will consider more production constraints in the future, such as setup, blocking, transportation, and delivery time. In addition, the realistic manufacturing processes always have multi-objective characteristics and variability. This requires us

to consider more optimization objectives and weigh the relationship between multiple objective functions. Furthermore, the possible emergencies during production are also required and studied to derive useful dynamic and rescheduling strategies.

**Author Contributions:** C.L.: conceptualization, methodology, data curation, software, validation, writing—original draft. B.Z.: conceptualization, methodology, software, validation, writing—original draft. Y.H.: conceptualization, methodology, software, validation, writing—original draft. Y.W.: conceptualization, methodology, supervision, writing—original draft. J.L.: conceptualization, methodology, visualization, investigation. K.G.: conceptualization, methodology, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China under grant numbers 61803192, 61973203, 62106073, 62173216, and 62173356. We are grateful for Guangyue Youth Scholar Innovation Talent Program support received from Liaocheng University, the Youth Innovation Talent Introduction and Education Program support received from the Shandong Province Colleges and Universities, and the Natural Science Foundation of Shandong Province under grant numbers ZR2021QE195.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author, upon reasonable request.

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

#### **Abbreviations**

Notations


#### **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.
