A Knowledge-Based Cooperative Differential Evolution Algorithm for Energy-Efficient Distributed Hybrid Flow-Shop Rescheduling Problem

Abstract

Due to the increasing level of customization and globalization of competition, rescheduling for distributed manufacturing is receiving more attention. In the meantime, environmentally friendly production is becoming a force to be reckoned with in intelligent manufacturing industries. In this paper, the energy-efficient distributed hybrid flow-shop rescheduling problem (EDHFRP) is addressed and a knowledge-based cooperative differential evolution (KCDE) algorithm is proposed to minimize the makespan of both original and newly arrived orders and total energy consumption (simultaneously). First, two heuristics were designed and used cooperatively for initialization. Next, a three-dimensional knowledge base was employed to record the information carried out by elite individuals. A novel DE with three different mutation strategies is proposed to generate the offspring. A local intensification strategy was used for further enhancement of the exploitation ability. The effects of major parameters were investigated and extensive experiments were carried out. The numerical results prove the effectiveness of each specially-designed strategy, while the comparisons with four existing algorithms demonstrate the efficiency of KCDE in solving EDHFRP.

1. Introduction

Productivity has become one of the core competencies in the manufacturing industry, particularly with the integration of the world economy. However, with the ongoing trend of mass customization and the accelerated variety of products, a growing number of manufacturing companies have realized the importance of environmentally friendly manufacturing. As the most energy-consuming and ${\mathrm{CO}}_2$ emission entity, China’s energy structure highly depends on coal. More than half of the electrical production is fueled by coal, while over 70% of power plants are coal-fired [1]. Energy-efficient scheduling of manufacturing industries is an effective way to reduce carbon emissions [2]. Therefore, it is of extraordinary significance to develop production schedules for intelligent manufacturing industries with energy-related constraints and objectives.

By extending the production model into the distributed environment and establishing multiple factories in various remote geographical locations, distributed manufacturing has been widely applied in various fields, including automotive, steel-making, and chemical processing [3]. In the past few decades, various evolutionary algorithms, such as artificial bee colony [4], iterated greedy algorithm [5], Jaya algorithm [6], and memetic algorithm [7], have been proposed to solve distributed shop scheduling problems with different constraints. The authors of [8] focused on the distributed parallel machine-scheduling problem and proposed a knowledge-based two-population optimization algorithm to minimize the energy consumption and total tardiness simultaneously. A distributed flow-shop group scheduling problem was investigated and solved by a cooperative co-evolutionary algorithm with a novel collaboration model and a reinitialization scheme in [9]. The authors of [10] studied a distributed heterogeneous hybrid flow-shop scheduling problem (FSP) under nonidentical time-of-use electricity tariffs, of which, the makespan and total electricity charges were considered as optimization objectives. Distributed hybrid differentiation FSP was investigated in [11], for which a distributed co-evolutionary memetic algorithm was proposed to minimize the makespan. Recently, the combination of distributed flow-shop and parallel machines, namely the distributed hybrid flow-shop scheduling problem (DHFSP), has emerged as a practical production model. In order to provide an efficient productive schedule for DHFSP, ref. [7] proposed a bi-population cooperative memetic algorithm. For energy-aware DHFSP with multiprocessor tasks, [12] proposed a multi-objective evolutionary algorithm based on decomposition while [13] designed a knowledge-based multi-objective algorithm.

Uncertainties always exist in real-life manufacturing. Unlike determinate processing, rescheduling aims to make schedules under dynamic environments, dealing with unpredicted events, such as new job arrivals and machine breakdowns. Rescheduling in FSP with variable processing times by using real-time information was studied in [14]. Rescheduling jobs in flexible job-shop scheduling problem (FJSP) is receiving more attention. FJSP (with newly arrived priority jobs) was investigated in [15], while machine breakdowns and FJSP recovery were considered in [16,17]. The authors of [18,19] took both real-time order acceptance and machine maintenance into consideration when rescheduling the FJSP. Furthermore, [20] provided energy-efficient schedules for FJSP with dynamic variations occurring frequently during production. The authors of [21] focused on rescheduling for the hybrid flow-shop scheduling problem (HFSP), in which energy consumption was regarded as an objective and the real-time order acceptance was considered. Moreover, [22] focused on steelmaking systems and proposed a hybrid fruit fly optimization algorithm to solve the rescheduling problem in HFSP with machine breakdowns and processing variations, while the release variation was considered in [23]. Moreover, rescheduling for a single-machine problem was studied in [24], in which the total waiting time was set as an objective and the newly arrived reworked jobs were considered. Rescheduling (in several identical parallel-machine problems) was investigated in [25]; the sequence-dependent setup times, release times, and newly arrived jobs were considered.

Inspired by the natural evolution of species, differential evolution (DE) is a simple (yet efficient) heuristic for global optimization over continuous spaces [26]; it has been improved in the past few decades and widely applied as a state-of-the-art global optimization technique in a variety of scheduling problems [27]. In [28], DE was combined with a local search mechanism to minimize the makespan and tardiness of FSP simultaneously. In order to solve permutation FSP while minimizing the makespan, optimization algorithms based on discrete DE were designed in [29], and a DE with the algebraic differential mutation was designed in [30], considering the total flowtime criterion. Time-of-use electricity tariffs were considered in [31], where a multi-objective discrete DE with novel mutations and crossover operations was developed. The authors of [32] combined DE with heuristic strategies in order to solve FJSP with outsourcing operations and job priority constraints, while a new mutation strategy based on searching was designed for DE in [33] to minimize the makespan in FJSP with reconfigurable machine tools. Moreover, DE was used for scheduling hybrid FSP [34] and JSP with fuzzy processing times [35,36]. In [8], DE performed cooperatively with NSGA-II on the population in parallel and provided efficient solutions for the distributed energy-efficient parallel machine-scheduling problem. Moreover, DE is used for solving other scheduling problems, such as robotic cell scheduling with batch-processing machines [37] and a single batch-processing machine with arbitrary job sizes and release times [38].

This study aims to solve the energy-efficient distributed hybrid flow-shop rescheduling problem (EDHFRP) by minimizing the makespan of both the first-order and the newly arrived order, along with the total energy consumption. EDHFRP can be considered as an extension of the existing DHFSP, which is studied in [7,12], by adding constraints, including energy-efficient and newly-arrived jobs. The mathematical model of EDHFRP is provided in Section 2.1, appending the calculation of energy consumption and the description of rescheduling on the basis of the model of DHFSP defined in [7,12]. In order to solve this problem, a knowledge-based cooperative differential evolution (KCDE) algorithm is proposed. In KCDE, two heuristics are designed for initialization according to the characteristics of EDHFRP. Then, a three-dimensional knowledge base is proposed to describe the information carried out by elite individuals. Moreover, three novel differential mutation strategies along with a crossover method are employed to the knowledge base for evolution. Additionally, a local intensification strategy was designed for elite individuals in order to further improve the quality of solutions. Finally, extensive comparisons and tests were carried out to verify the efficiency and effectiveness of the proposed KCDE in solving EDHFRP.

The remainder of this paper is divided into four parts. Section 2 describes EDHFRP in detail and provides the mathematical model. Every strategy and design of KCDE is introduced in Section 3. The numerical results and discussion are shown in Section 4. Finally, the conclusion provides a brief summary and critique of the findings in Section 5.

2. Distributed Hybrid Flow-Shop Rescheduling Problem

2.1. Problem Description

The indices, parameters, variables and decision variables used in the following description is shown as follows.

Indices	Description
F	Number of factories.
f	Index of factories, $f\in \{1,2,\cdots ,F\}$.
$O_1$	The first order.
$O_2$	The newly arrived order during the processing of $O_1$.
$n_1$	Number of jobs in $O_1$.
$n_2$	Number of jobs in $O_2$.
$n_p$	Number of jobs under processing at time $t_a$ in $O_1$.
$n_u$	Number of unprocessed jobs at time $t_a$ in $O_1$.
$j_1$	Index of jobs in $O_1$, $j_1\in \{1,2,\cdots ,n_1\}$.
$j_2$	Index of jobs in $O_2$, $j_2\in \{1,2,\cdots ,n_2\}$.
$j_p$	Index of jobs under processing at time $t_a$, $j_p\in \{1,2,\cdots ,n_p\}$.
$j_u$	Index of unprocessed jobs at time $t_a$, $j_u\in \{1,2,\cdots ,n_u\}$.
j	Index of all jobs, $j\in \{1,2,\cdots ,n_1+n_2\}$.
s	Number of stages.
k	Index of stages, $k\in \{1,2,\cdots ,s\}$.
$m_{f,k}$	Number of machines for stage k in factory f.
i	Index of machines, $i\in \{1,2,\cdots ,m_{f,k}\}$.

Parameters	Description
$t_a$	The arrival time of $O_2$.
$P_{j,k}$	Processing time of job j at stage k.
$E{P}_{f,k,i}$	Energy consumption of machine i per unit time at stage k in factory f in
	processing mode.
$EI$	Energy consumption per unit time in idle mode.

Variables	Description
$EC{P}_{f,k,i}$	Energy consumption of machine i in processing mode at stage k in factory f.
$EC{I}_{f,k,i}$	Energy consumption of machine i in idle mode at stage k in factory f.
$S_{j,k}$	Begin time of processing job j at stage k.
$C_{j,k}$	Completion time of processing job j at stage k.
$M{S}_1$	Makespan of $O_1$.
$M{S}_2$	Makespan of $O_2$.
$TEC$	Total energy consumption.

Decision Variables	Description
$x_{j,f}$	Binary variable whose value equals 1 when job j is assigned to
	factory f or 0 otherwise.
$y_{j,f,k,i}$	Binary variable whose value equals 1 when job j is assigned to
	machine i at stage k in factory f or 0 otherwise.
$z_{j,j^{\prime },f,k,i}$	Binary variable whose value equals 1 when job j is processed
	before $j^{\prime }$ on machine i at stage k in factory f or 0 otherwise.

In a typical DHFSP, the total number of n jobs need to be processed on a set of machines. Each job includes s stages and there are one or more machines for each stage. Moreover, the machines belong to F-geographically dispersed factories. Each job needs to be assigned to one of the factories before processing, and should not be removed from the certain factory until completion. Therefore, in order to schedule DHFSP, three sub-problems demand solving, namely factory assignment, machine selection, and job sequencing.

In EDHFRP, during the processing of $n_1$ jobs in the first-order $O_1$, another order $O_2$ with $n_2$ jobs arrives at time $t_a$ and needs to be completed as soon as possible. At time $t_a$, the possible states of jobs in $O_1$ are completed, under processing, and unprocessed. For jobs in the first two states, the original schedule is kept, while the unprocessed jobs in $O_1$ need rescheduling along with the new jobs in $O_2$. However, during the rescheduling, the available time of each machine changes from 0 to the completion time of jobs under processing at time $t_a$. Furthermore, the unprocessed jobs in $O_1$ have to remain in the factories to which they are assigned, while the newly arrived jobs in $O_2$ can be assigned to any factory. Therefore, the complexity of scheduling EDHFRP is much higher than the typical DHFSP.

The Gantt chart of a solution of an instance with $n_1=8$, $n_2=5$, $s=3$, $F=2$, and $t_a=30$ is illustrated in Figure 1. The original schedule for $O_1$ is shown in Figure 1a; it can be seen that job 1 and job 7 in factory 1, as well as job 8 in factory 2 are unprocessed when new order $O_2$ arrives at $t_a$. Thus, the three unprocessed jobs along with jobs in $O_2$ are rescheduled, while job 1 and job 7 still remain in factory 1 and job 8 in factory 2. The processing plan after rescheduling is shown in Figure 1b, in which the rectangles labeled II-1∼II-5 denote the jobs in $O_2$. The makespan of $O_1$ and $O_2$ are 163 and 195 min, respectively, and the total energy consumption can be calculated as 134.08 kJ.

2.2. Mathematical Model

The mixed integer linear programming model for EDHFRP minimizing $M{S}_1$, $M{S}_2$ and $TEC$ is formulated as follows:

$$min(M{S}_1,M{S}_2,TEC)$$

(1)

Subject to:

$$\sum _{f=1}^F{x}_{j,f}=1 \forall j$$

(2)

$$\sum _{i=1}^{m_{f,k}}{y}_{j,f,k,i}=1 \forall f,j,k$$

(3)

$$S_{j_1,1}\ge 0 \forall j_1$$

(4)

$$n_p+n_u=n_1$$

(5)

$$C_{j,k}=S_{j,k}+P_{j,k} \forall j,k$$

(6)

$$S_{j,k+1}\ge C_{j,k} \forall j \forall k\in \{1,2,\cdots ,s-1\}$$

(7)

$$z_{j,j^{\prime },f,k,i}+z_{j^{\prime },j,f,k,i}\le 1 \forall j\ne j^{\prime },f,k,i\in \{1,2,\cdots ,m_{f,k}\}$$

(8)

$$z_{j,j^{\prime },f,k,i}+z_{j^{\prime },j,f,k,i}\ge y_{j,f,k,i}+y_{j^{\prime },f,k,i}-1 \forall j\ne j^{\prime },f,k,i$$

(9)

$$S_{j^{\prime },k}-C_{j,k}+{10}^6\times (3-y_{j,f,k,i}-y_{j^{\prime },f,k,i}-z_{j,j^{\prime },f,k,i})\ge 0 \forall j\ne j^{\prime },f,k,i$$

(10)

$$M{S}_1\ge C_{j_1,s} \forall j_1$$

(11)

$$M{S}_2\ge C_{j_2,s} \forall j_2$$

(12)

$$EC{P}_{f,k,i}=\sum _{j=1}^n{y}_{j,f,k,i}\times E{P}_{f,k,i}\times P_{j,k} \forall f,k,i\in \{1,2,\cdots ,m_{f,k}\}$$

(13)

$$\begin{array}{cc}\hfill EC{I}_{f,k,i}=EI\times [\underset{j}{max}(C_{j,k}& \times y_{j,f,k,i})-\underset{j}{min}(S_{j,k+1}\times y_{j,f,k,i})-\sum _{j=1}^n{y}_{j,f,k,i}\times P_{j,k}]\hfill \\ \hfill & \forall f,k,i\in \{1,2,\cdots ,m_{f,k}\}\hfill \end{array}$$

(14)

$$TEC=\sum _{f=1}^F\sum _{k=1}^s\sum _{i=1}^{m_{f,k}}(EC{P}_{f,k,i}+EC{I}_{f,k,i})$$

(15)

$$x_{j,f}\in \{0,1\} \forall j,f$$

(16)

$$y_{j,f,k,i}\in \{0,1\} \forall j,f,k,i\in \{1,2,\cdots ,m_{f,k}\}$$

(17)

$$z_{j,j^{\prime },f,k,i}\in \{0,1\} \forall j,j^{\prime },f,k,i$$

(18)

where Equation (1) describes the three optimization objectives, Equation (2) assures that only one factory must be selected for each job, Equation (3) ensures that the processing of each stage of each job must happen on only one machine, Equation (4) guarantees all jobs in $O_1$ must be processed after time 0, Equation (5) makes sure that each job in $O_1$ is either under-processed or unprocessed at time $t_a$, Equation (6) defines the calculation of the completion time, Equation (7) ensures that each stage can only be processed after the completion of the previous one, Equations (8)–(10) guarantee that each machine can process—at most—one job at any time, Equations (11)–(12) define the makespan of two orders, Equations (13)–(15) calculate the total energy consumption. Equations (16)–(18) indicate the binary decision variables.

3. KCDE for EDHFRP

The classical DE refers to an efficient and powerful population-based stochastic search technique for solving optimization problems over a continuous space, aiming to evolve a population of D-dimensional parameter vectors, which encode the candidate solutions toward the global optimum. However, as a discrete optimization problem, the parameter vectors built for job-shop scheduling problems are limited to being positive integers according to the constraints; thus, it is difficult for DE to reach its full potential when directly used in solving job-shop scheduling problems. The knowledge base proposed in this paper is capable of making an explicit connection between continuous space and discrete parameter vectors. In this section, the representation of the solution is first introduced briefly, following by the details of four vital components of KCDE, including hybrid initialization, knowledge base, cooperative differential mutation, and local intensification.

3.1. Solution Representation

For each job in EDHFRP, the factory assignment should be determined first, then a machine is selected for each stage, while the processing sequence on each machine should be determined at the same time. In this paper, the efficient and widely used permutation-encoding method is employed to represent the job sequence in each factory. Machine assignments are obtained through decoding operations according to the job sequence and factory assignment. Therefore, the solutions of EDHFRP can be encoded as $\{N_1,\cdots ,N_F,{\Pi }_1,\cdots ,{\Pi }_F\}$, where $N_f,f\in \{1,2,\cdots ,F\}$ denotes the number of jobs assigned to factories f, ${\Pi }_f,f\in \{1,2,\cdots ,F\}$ represents the processing sequence of jobs at the first stage in factory f. Considering the feasibility of the solutions, each job should occur in the job sequence of any factory once and only once.

To be specific, the corresponding solution of the Gantt chart in Figure 1a—before new order arrival—is shown in Figure 2a, from which it can be seen that jobs 3, 5, 7, and 1 are processed in factory 1 sequentially, while the first job processed in factory 2 is job 4, followed by jobs 2, 6, and 8. When $O_2$ arrives at time $t_a$, jobs 3, 5, 4, and 2 are being processed, while jobs 7, 1, 6, and 8 remain unprocessed. Therefore, the 4 unprocessed jobs along with the 5 new jobs are rescheduled as in Figure 2b, in which jobs 7, 1, 6, and 8 are still in the factories they were assigned to at first.

3.2. Hybrid Initialization

In order to strike a balance between the quality and variety of the initial population, three different initialization strategies, including two EDHFRP characteristic-based heuristics and a random method, were designed and employed in a cooperative manner.

It has been conclusively shown that the Nawaz–Enscore–Ham (NEH) heuristic [39] is one of the most effective heuristics for solving FSP [40]. Recently, the variants of NEH heuristics have been widely used in DFSP and DHFSP [5,6,41,42]. In this paper, a greedy NEH heuristic (GNEH) and a NEH heuristic based on lower bounds (NEHLB) are proposed.

For both heuristics, the initial job sequences are generated randomly and the jobs are arranged sequentially. The main method of GNEH is described as Algorithm 1 in detail. First, each job is inserted into all possible position in each factories. Next, the values of $M{S}_1$ and $TEC$ of each insertion are calculated. Then, the insertion with the best result is selected. Moreover, $NP/3$ individuals are generated through GNEH.

Algorithm 1: Greedy NEH heuristic

The NEHLB is designed as Algorithm 2 and implemented to generate $NP/3$ individuals. For each job, the current $L{B}_{M{S}_1}$ and $L{B}_{TEC}$ of each factory are calculated first. Next, the factory with the smallest value is chosen. Then, the job is inserted into one of the possible positions in a certain factory (randomly). In order to increase the diversity of the initial population, the rest of the population is generated randomly.

Algorithm 2: NEH heuristic with biased optimization

3.3. 3D Knowledge Base

Taking the features of DE and EDHFRP into consideration, a three-dimensional knowledge base is designed to store the information of elite individuals. In the proposed knowledge base, both x-axis and y-axis represent the jobs, while the z-axis represents the factories. Therefore, the $xy$-plane is a square, whose side lengths are equal to the number of jobs n. The value of the knowledge base at position $(x_0,y_0,z_0)$, where $x_0,y_0\in 1,2,\cdots ,n,x_0\ne y_0$ and $z_0\in \{1,2,\cdots ,F\}$, denotes the probability that both job $x_0$ and $y_0$ are assigned to the same factory $z_0$ while job $x_0$ is processed before job $y_0$, while the value of position $(x_0,x_0,z_0)$ means the probability that job $x_0$ is assigned to factory $z_0$.

The knowledge base is initialized as Equation (19), where $j_1,j_2\in \{1,2,\cdots ,n_1\}$ and $f\in \{1,2,\cdots ,F\}$.

$$K{B}_0(j_1,j_2,f)=\left\{\begin{array}{c}0.5, \mathrm{if} j_1\ne j_2\hfill \\ 1, \mathrm{if} j_1=j_2\hfill \end{array}\right.$$

(19)

At the end of each generation, the knowledge base is updated according to the non-dominated solutions. For each elite individual $po{p}_e$ in $E\left(g\right)$, $e\in 1,2,\cdots ,n_E$, an incremental matrix ${\delta }_e$ is calculated according to Equations (20) and (21) where $j_1,j_2\in \{1,2,\cdots ,n_F\}$, $j_1\ne j_2$, and $n_E$ is the number of non-dominated solutions.

$${\delta }_e(j_1,j_2,f)=\left\{\begin{array}{c}1, \mathrm{if} j_1 \mathrm{is} \mathrm{processed} \mathrm{before} j_2 \mathrm{in} \mathrm{factory} f\hfill \\ 0, \mathrm{otherwise}\hfill \end{array}\right.$$

(20)

$${\delta }_e(j,j,f)=\left\{\begin{array}{c}1, \mathrm{if} j \mathrm{is} \mathrm{assigned} \mathrm{to} \mathrm{factory} f\hfill \\ 0, \mathrm{otherwise}\hfill \end{array}\right.$$

(21)

Finally, the knowledge base is updated via Equation (22) and normalized as Equation (23).

$$K{B}_g=(1-\alpha )\times K{B}_{g-1}+\alpha \times \frac{{\sum }_{e=1}^{n_E}{\delta }_e}{n_E}$$

(22)

$$K{B}_g(x,y,z)=\left\{\begin{array}{c}\frac{K{B}_g(x,y,z)}{{\sum }_{z=1}^{F}KB(x,y,z)}, \mathrm{if} x=y\hfill \\ \frac{K{B}_g(x,y,z)}{{\sum }_{x=1}^n{\sum }_{y=1}^{n}KB(x,y,z)}, \mathrm{if} x\ne y\hfill \end{array}\right.$$

(23)

Through the operations above, the useful information of better-performed solutions is transformed from the discrete model into the continuous space, which is more suitable to work DE to its advantage. Moreover, the processing priority between different jobs and the information of the factory assignment are linked together.

3.4. Cooperative DE

As the most important step of DE, the mutation operator provides new individuals for the population and considerable research studies have proposed modifications in mutation strategies in the past two decades [43]. A total of six most popular types of mutation operators of DE and their characteristics are represented in [44]. As a result of the performance insufficiency of each single mutation strategy, the cooperation among several different strategies has been used to solve complicated scheduling problems [33,45,46]. Since the mutation operator DE/rand/1 has the most exploration properties but very limited exploitation abilities, DE/best/1 has the most exploitation ability but may end up in premature convergence, while DE/current-to-best/1 has both the exploration and exploitation properties according to its utilization of the memory set of individuals, the cooperation among them may avoid being trapped in local optima and ensure the balance between the diversification and intensification when searching for the optimal solution. In this paper, three differential mutation strategies, i.e., DE/rand/1, DE/best/1, and DE/current-to-best/1, are used cooperatively with the knowledge base according to the weight vector $\mathbf{w}=[w_1,w_2,w_3]$. The main method of cooperative DE is shown in Algorithm 3. The individuals for each DE strategy are calculated as $N{P}_1=\mathrm{round}(NP\times w_1)$, $N{P}_2=\mathrm{round}(NP\times w_2)$, and $N{P}_3=NP-N{P}_1-N{P}_2$.

Algorithm 3: Cooperative Differential Evolution

3.4.1. Mutation

For each individual in $Po{p}_1\left(g\right)$, the DE/rand/1 strategy is carried out on $K{B}_g$, for which three individuals are selected from $Po{p}_1\left(g\right)$ randomly. Next, the best of them is recorded as ${\mathbf{x}}_{r_0}$ according to the non-dominated sort, while the worst individual is named ${\mathbf{x}}_{r_2}$, and the rest are denoted as ${\mathbf{x}}_{r_1}$. Then, a variant matrix can be obtained through differential mutation, as described in Equations (24) and (25), where $x,y\in \{1,2,\cdots ,n\}$, $x\ne y$, $z\in \{1,2,\cdots ,F\}$ and $F_m$ denote the mutation factors.

$${\Delta }_1(x,y,z)=\left\{\begin{array}{c}1, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{processed} \mathrm{before} \mathrm{job} y \mathrm{in} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_{r_0}\hfill \\ F_m, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{processed} \mathrm{before} \mathrm{job} y \mathrm{in} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_{r_1}\hfill \\ -F_m, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{processed} \mathrm{before}\mathrm{job} y \mathrm{in} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_{r_2}\hfill \end{array}\right.$$

(24)

$${\Delta }_1(x,x,z)=\left\{\begin{array}{c}1, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{assigned} \mathrm{to} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_{r_0}\hfill \\ F_m, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{assigned} \mathrm{to} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_{r_1}\hfill \\ -F_m, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{assigned} \mathrm{to} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_{r_2}\hfill \end{array}\right.$$

(25)

The DE/best/1 strategy is carried out on $K{B}_g$ for each individual in $Po{p}_2\left(g\right)$. A total of three individuals, including the best one in $Po{p}_2\left(g\right)$ and two other ones randomly selected, are selected and denoted as ${\mathbf{x}}_{best}$, ${\mathbf{x}}_{r_1}$ and ${\mathbf{x}}_{r_2}$. Then, the differential mutation is carried out as Equations (24) and (25), in which the ${\mathbf{x}}_{r_0}$ is replaced by ${\mathbf{x}}_{best}$, thus the variant matrix ${\Delta }_2$ is received.

As for the DE/current-to-best/1, the current individual along with the best one in $Po{p}_3\left(g\right)$ and two randomly selected ones are chosen and recorded as ${\mathbf{x}}_i$, ${\mathbf{x}}_{best}$, ${\mathbf{x}}_{r_1}$ and ${\mathbf{x}}_{r_2}$ respectively. Then, the variant matrix is generated through Equations (26) and (27).

$${\Delta }_3(x,y,z)=\left\{\begin{array}{c}2\times F_m, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{processed} \mathrm{before} \mathrm{job} y \mathrm{in} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_{best}\hfill \\ 1-2\times F_m, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{processed} \mathrm{before} \mathrm{job} y \mathrm{in} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_i\hfill \\ F_m, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{processed} \mathrm{before} \mathrm{job} y \mathrm{in} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_{r_1}\hfill \\ -F_m, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{processed} \mathrm{before} \mathrm{job} y \mathrm{in} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_{r_2}\hfill \end{array}\right.$$

(26)

$${\Delta }_3(x,x,z)=\left\{\begin{array}{c}2\times F_m, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{assigned} \mathrm{to} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_{best}\hfill \\ 1-2\times F_m, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{assigned} \mathrm{to} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_i\hfill \\ F_m, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{assigned} \mathrm{to} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_{r_1}\hfill \\ -F_m, \mathrm{if} \mathrm{job} x \mathrm{is} \mathrm{assigned} \mathrm{to} \mathrm{factory} z \mathrm{in} {\mathbf{x}}_{r_2}\hfill \end{array}\right.$$

(27)

3.4.2. Crossover

The crossover operator constructs a new trial individual under the guidance of a target individual and a mutant individual [43]. Most of the existing crossover operators of DE are binomial variant, which is essentially suitable for solving continuous optimization problems. However, as a discrete problem, job-shop scheduling problems have various limits on structuring individuals. Most of the existing research using DE to solve scheduling problems have paid attention to the modification of crossover strategies or the structure of individuals instead of alteration in both of them [47], which may cause mismatching between the problem and the algorithm. Thus, on the basis of the knowledge base, a problem-specific crossover operator is proposed as follows.

With ${\Delta }_d$ and $K{B}_g$, an offspring matrix can be obtained through crossover, which is shown in Equation (28).

$$K{B}_{gd}=(1-CR)\times K{B}_g+CR\times {\Delta }_d$$

(28)

where $CR$ is the crossover probability, ${\Delta }_d,d\in \{1,2,3\}$ indicates the three variant matrices provided by three DE strategies, and $K{B}_{gd}$ is the corresponding offspring matrix.

Finally, an offspring individual is sampled from $K{B}_{g1}$ according to the Roulette wheel selection strategy. The sampling operation mainly consists of two steps, choosing the factory for every job and sorting the jobs in each factory. In the first step, the certain factory f with ${max}_{f}K{B}_{gd}(j,j,f)$ is selected for each job j. Then, jobs assigned to the same factory are chosen one by one via the roulette wheel selection strategy according to the value of $K{B}_{gd}(j,j^{\prime },f)$; thus, jobs with higher values are more likely to be processed earlier.

3.4.3. Selection

The offspring generated through DE/rand/1, DE/best/1, and DE/current-to-best/1, denoted as $Po{p}_{New1}\left(g\right)$, $Po{p}_{New2}\left(g\right)$, and $Po{p}_{New3}\left(g\right)$, respectively, are gathered together as $Po{p}_{New}\left(g\right)$, with a size of $NP$. Next, the newly generated population $Po{p}_{New}\left(g\right)$ and the current population $Pop\left(g\right)$ are united and the best $NP$ individuals are selected according to non-dominated sorting. During the selection, the chosen individuals belonging to $Po{p}_{New1}\left(g\right)$, $Po{p}_{New2}\left(g\right)$, and $Po{p}_{New3}\left(g\right)$ are recorded as $Nu{m}_1$, $Nu{m}_2$, and $Nu{m}_3$, respectively. Then, the weight vector is updated as Equation (29).

$$\mathbf{w}=[\frac{Nu{m}_1}{Nu{m}_1+Nu{m}_2+Nu{m}_3},\frac{Nu{m}_2}{Nu{m}_1+Nu{m}_2+Nu{m}_3},\frac{Nu{m}_3}{Nu{m}_1+Nu{m}_2+Nu{m}_3}]$$

(29)

According to the descriptions of the three DE strategies, the minimum number of individuals for each group is four. In order to guarantee the feasibility of the cooperative DE, the lower bound of $w_1$, $w_2$, and $w_3$ in the weight vector is set as 0.1. In other words, the value of elements in $\mathbf{w}$ smaller than 0.1 will be improved to 0.1, while the other elements will be decreased in proportion.

3.5. Local Intensification

In order to further enhance the quality of solutions, a local intensification is designed and carried out on non-dominated individuals at the end of each generation. The main method of local intensification is described in Algorithm 4. For each non-dominated individual, the factories with the maximum and minimum completion times are first selected, denoted as the critical factory $f_c$ and the easiest factory $f_e$, respectively. Next, a job assigned to $f_c$ is chosen randomly, which is then inserted into all possible positions in $f_e$. After that, the objective values of each insertion are calculated and the best one is kept as a neighbor. Finally, the $NP$ best individuals are selected from the union of the current population and neighbors.

Algorithm 4: Local intensification

3.6. Rescheduling Strategy

When the new order $O_2$ arrives, the following tasks are completed before carrying out the KCDE.

• The state of jobs in $O_1$ at time $t_a$ are recorded, and the unprocessed jobs are counted.
• The unprocessed jobs are then put with jobs in $O_2$, together forming the new order. The factory assignments of unprocessed jobs are also recorded, which are used as the constraints during the rescheduling stage since the jobs cannot be transformed into other factories once assigned.
• The available machine times are updated into the completion times of jobs, which are processed at time $t_a$.

3.7. Framework of KCDE

The framework of the proposed KCDE is depicted in Figure 3. Three differential evolution strategies are employed cooperatively on the knowledge base, and a self-adaptive weight vector is designed to control the proportion of the application of three strategies according to the quality of the offspring generated by them. Moreover, several problem-specific strategies, such as hybrid initialization, update of the knowledge base, and local intensification, are implemented. The performance of such a properly designed algorithm in solving EDHFRP is highly expected.

4. Numerical Results and Comparisons

4.1. Experimental Settings

A considerable amount of instances are generated and employed to test the performance of KCDE in solving EDHFRP. Each instance has four decisive factors for its scale, including the number of factories F, the number of jobs in the first order $n_1$, the number of jobs in the new order $n_2$, and the number of stages s. In this paper, the value of F is collected from $\{3,4,5\}$, the value of $n_1$ is from $\{30,50,80\}$, while $n_2$ is selected from $\{20,30,50\}$ and s is from $\{4,5,6\}$. Moreover, the processing time and energy consumption of the per-unit time is sampled from the instances provided by early research for similar problems [12,41]. Moreover, 10 instances are generated for each combination of F, $n_1$, $n_2$, and s; thus, a total of $3\times 3\times 3\times 3=81$ instances in 81 different scales are used in the following tests and comparisons. The arrival time of $O_2$ is set randomly, while instances with the same scales own the same $t_a$.

Generally, more running time is needed in each iteration for a scheduling instance with higher complexity, while the operation is faster when it comes to an easy instance. In order to guarantee the reasonableness and fairness of the experiments on instances with different scales, the maximum running time is set according to the four decisive factors for the scale. [d=R1.]Therefore, for Foreach compared algorithm and testing algorithm, the maximum running time for the first order is set as $0.1\times F\times n_1\times s$ seconds, while it is set as $0.1\times F\times n_2\times s$ seconds after $O_2$’s arrival. To be fair, all algorithms were implemented through MATLAB2021b and conducted in the same computer environment with 12th Gen Intel Core i7-12700K, 3.61GHz, 32GB RAM, and Windows 10 OS.

Three commonly used metrics in the multi-optimization field, i.e., hyper-volume (HV), C metric (CM), and generational distance (GD) [48] are utilized to evaluate the performance of the algorithms, which focus on measuring the convergence and diversity of the Pareto set approximations.

4.2. Model Validation

In order to validate the mathematical model provided in Section 2.2, 10 small-scale instances are generated randomly, and the exact solver IBM ILOG CPLEX 12.9 is used to compare with the proposed algorithm. For each comparison, both CPLEX and KCDE run 10 times independently. In each run, the CPU time limit of CPLEX is set to 1 h while the stopping criterion of KCDE is set as mentioned above. Since the rescheduling solution is obtained on the basis of the schedule of the original order, and depends on the producing state at time $t_a$, the comparison on each instance is divided into two parts.

In the first part, the original schedule $O_1$ is scheduled by both CPLEX and KCDE, the results ($M{S}_1$ and $TEC$) are compared. In the second part, the scheduling plan of the original order obtained by CPLEX is used for both the solver and the algorithm to solve the rescheduling problem; the results $M{S}_1$, $M{S}_2$, and $TEC$ are compared again.

The results of the first comparison are shown in the left column of

Table 1. Comparisons between CPLEX and KCDE.

Instance	Original Order				Rescheduling Plan
	CPLEX		KCDE		CPLEX			KCDE
f2n4+3s2	117.70	33.67	117.70	33.27	119.10	107.40	68.34	118.20	107.40	65.14
f2n4+3s2	118.50	33.02	118.20	32.66	143.90	143.00	59.62	143.80	142.50	58.70
f2n4+3s3	123.00	55.35	123.00	55.27	166.00	179.70	103.55	164.10	178.20	99.95
f2n4+3s3	146.20	61.90	145.60	61.43	153.00	195.10	96.44	153.00	195.00	94.55
f2n5+3s2	82.30	44.07	82.00	43.23	121.40	119.00	76.44	121.30	118.10	74.04
f2n5+3s2	150.20	43.42	150.20	43.39	149.60	151.50	79.13	149.60	151.50	76.13
f2n5+3s3	126.90	58.47	126.90	58.48	146.00	102.90	115.95	145.70	102.40	113.76
f2n5+3s3	200.60	74.07	200.60	73.04	228.90	176.80	121.70	228.60	176.00	120.34
f3n5+3s2	96.30	40.88	95.90	40.71	99.10	105.00	99.84	98.50	103.00	98.47
f3n5+3s2	139.50	40.53	139.50	40.51	140.00	142.30	96.76	138.80	140.10	96.32

The better value of each instance is marked in bold.

, from which it can be seen that small differences exist between the results obtained by CPLEX and KCDE, while KCDE achieves most of the better results. The right column in Table 1 provides the results of the second comparison, in which the objective values given by KCDE are better than those provided by CPLEX in all instances. Therefore, the mathematical model of EDHFRP is validated.

4.3. Parameter Setting

There exist three important parameters in the proposed KCDE: (1) the mutation factor in the differential mutation operation $\left(F\right)$, (2) the probability of crossover in DE $\left(CR\right)$, (3) and the weight factor used for updating the knowledge base $\left(\alpha \right)$. In order to inspect the effects of these parameters on the performance of KCDE and choose the best combination accordingly, the famous method, the design of the experiment (DOE), is carried out as follows.

In order to be more accurate, four levels are set for each parameter. According to the empirical conclusions provided by early studies on DE, the values of F and $CR$ are set as $F\in \{0.5,0.55,0.6,0.65\}$ and $CR\in \{0.2,0.25,0.3,0.35\}$. On the basis of the characteristics of EDHFRP and the design of the knowledge base, the value of $\alpha$ is set as $\alpha \in \{0.75,0.8,0.85,0.9\}$. Then, 16 different combinations of $(F,CR,\alpha )$ are tested on a group of 27 typical instances, including three small-sized instances, three moderate-sized instances, and three large-sized instances for each $f\in \{3,4,5\}$, according to the orthogonal array $L_{16}\left(4^3\right)$.

To be persuasive, the HV indicator is used to represent the performance of the algorithm, and KCDE runs 10 times on each instance independently with each combination. The average HV value is recorded and shown in

Table 2. Orthogonal Arrays and HV Values.

Experiment	Factor Level			HV
Number	F	$\mathbf{CR}$	$\mathit{\alpha }$
1	1	1	1	0.77799
2	1	2	2	0.76874
3	1	3	3	0.78359
4	1	4	4	0.78363
5	2	1	2	0.77442
6	2	2	1	0.77722
7	2	3	4	0.78411
8	2	4	3	0.78266
9	3	1	3	0.77947
10	3	2	4	0.78020
11	3	3	1	0.78163
12	3	4	2	0.78488
13	4	1	4	0.77749
14	4	2	3	0.77721
15	4	3	2	0.77504
16	4	4	1	0.78175

. The combination of the minimum values of the objectives among all results for each instance is regarded as the lower bound, while the maximum ones are combined and regarded as the higher bound when calculating HV.

Next, the average value of HV is calculated for each level of each parameter, and the certain level with the largest HV value is chosen. The average HVs and the rank of the parameters are shown in

Table 3. Average HVs and the Rank of Each Parameter.

Level	F	$\mathbf{CR}$	$\mathit{\alpha }$
1	0.77849	0.77734	0.77965
2	0.77960	0.77584	0.77577
3	0.78154	0.78109	0.78073
4	0.77787	0.78323	0.78136
Delta	0.00367	0.00738	0.00559
Rank	3	1	2

The best value of each level is marked in bold.

, where Delta is the difference between the largest and the smallest HV value for each parameter. The parameter with a larger Delta has a more significant effect on the performance of KCDE. Meanwhile, the main effects of the parameters are illustrated in Figure 4.

It can be seen from Table 3 that $CR$ has the greatest influence on KCDE when solving EDHFRP in different scales, followed by $\alpha$ and $CR$. As illustrated in Figure 4, $CR$ denotes the proportion that the offspring learns from the trial individual rather than its parent, and a small value of $CR$ may result in a rather slow evolutionary process. $\alpha$ determines how much the knowledge base learns from the elite individuals in each generation; from the trend depicted in Figure 4, it can be seen that a small $\alpha$ results in lower HV, which is caused by the loss of useful information. F controls the influence of parent vectors on the trial vector during differential mutation; a small F is more likely to decrease the diversity of the population, while the F of a too-large value will lead the population into slower convergence, as shown in Figure 4. Therefore, according to the analysis above and the results in Table 3, $F=0.6$, $CR=0.35$, and $\alpha =0.9$ are used in the following tests and experiments.

4.4. Effect of Hybrid Initialization

The proposed KCDE is compared to its variant KCDE-Ran, in which the initial population is generated randomly in order to investigate the effectiveness of hybrid initialization. For each instance, each algorithm runs 10 times independently. The average C(KCDE, KCDE-Ran) and C(KCDE-Ran, KCDE) values of the instances with the same scale are calculated and provided in

Table 4. Average C Metric Values of KCDE and KCDE-Ran.

Dataset	${\mathit{C}}_{1,2}$1	${\mathit{C}}_{2,1}$2	p	Dataset	${\mathit{C}}_{1,2}$	${\mathit{C}}_{2,1}$	p	Dataset	${\mathit{C}}_{1,2}$	${\mathit{C}}_{2,1}$	p
f3n30+20s4	0.43	0.31	0.02	f4n30+20s4	0.53	0.24	0.02	f5n30+20s4	0.50	0.23	0.02
f3n30+20s5	0.47	0.34	0.02	f4n30+20s5	0.44	0.32	0.02	f5n30+20s5	0.45	0.29	0.02
f3n30+20s6	0.43	0.32	0.02	f4n30+20s6	0.43	0.29	0.02	f5n30+20s6	0.39	0.31	0.02
f3n30+30s4	0.47	0.30	0.02	f4n30+30s4	0.50	0.30	0.02	f5n30+30s4	0.45	0.28	0.02
f3n30+30s5	0.52	0.27	0.02	f4n30+30s5	0.53	0.24	0.00	f5n30+30s5	0.56	0.24	0.00
f3n30+30s6	0.53	0.28	0.02	f4n30+30s6	0.46	0.36	0.02	f5n30+30s6	0.53	0.25	0.02
f3n30+50s4	0.61	0.19	0.00	f4n30+50s4	0.53	0.28	0.02	f5n30+50s4	0.58	0.23	0.02
f3n30+50s5	0.54	0.26	0.02	f4n30+50s5	0.53	0.27	0.02	f5n30+50s5	0.54	0.23	0.02
f3n30+50s6	0.49	0.35	0.02	f4n30+50s6	0.48	0.31	0.02	f5n30+50s6	0.61	0.19	0.00
f3n50+20s4	0.42	0.28	0.02	f4n50+20s4	0.47	0.26	0.02	f5n50+20s4	0.44	0.24	0.00
f3n50+20s5	0.53	0.28	0.02	f4n50+20s5	0.41	0.28	0.02	f5n50+20s5	0.46	0.32	0.02
f3n50+20s6	0.48	0.29	0.02	f4n50+20s6	0.50	0.24	0.00	f5n50+20s6	0.46	0.32	0.02
f3n50+30s4	0.54	0.19	0.02	f4n50+30s4	0.55	0.24	0.02	f5n50+30s4	0.48	0.29	0.02
f3n50+30s5	0.57	0.27	0.02	f4n50+30s5	0.46	0.33	0.02	f5n50+30s5	0.54	0.26	0.02
f3n50+30s6	0.51	0.29	0.02	f4n50+30s6	0.51	0.26	0.02	f5n50+30s6	0.49	0.30	0.02
f3n50+50s4	0.60	0.22	0.02	f4n50+50s4	0.58	0.23	0.00	f5n50+50s4	0.57	0.22	0.00
f3n50+50s5	0.63	0.20	0.02	f4n50+50s5	0.56	0.26	0.02	f5n50+50s5	0.62	0.19	0.00
f3n50+50s6	0.59	0.24	0.02	f4n50+50s6	0.57	0.23	0.00	f5n50+50s6	0.59	0.23	0.00
f3n80+20s4	0.47	0.19	0.02	f4n80+20s4	0.45	0.26	0.02	f5n80+20s4	0.40	0.22	0.02
f3n80+20s5	0.50	0.28	0.02	f4n80+20s5	0.43	0.22	0.02	f5n80+20s5	0.52	0.22	0.00
f3n80+20s6	0.45	0.29	0.02	f4n80+20s6	0.49	0.24	0.02	f5n80+20s6	0.36	0.23	0.02
f3n80+30s4	0.48	0.18	0.02	f4n80+30s4	0.47	0.22	0.02	f5n80+30s4	0.48	0.15	0.00
f3n80+30s5	0.55	0.24	0.02	f4n80+30s5	0.46	0.27	0.02	f5n80+30s5	0.50	0.29	0.02
f3n80+30s6	0.47	0.24	0.02	f4n80+30s6	0.56	0.20	0.00	f5n80+30s6	0.44	0.28	0.02
f3n80+50s4	0.52	0.23	0.02	f4n80+50s4	0.54	0.24	0.02	f5n80+50s4	0.62	0.16	0.00
f3n80+50s5	0.57	0.19	0.02	f4n80+50s5	0.50	0.26	0.02	f5n80+50s5	0.60	0.23	0.02
f3n80+50s6	0.52	0.21	0.02	f4n80+50s6	0.55	0.19	0.02	f5n80+50s6	0.60	0.22	0.02

¹C_1,2 = C(KCDE,KCDE-Ran); ² C_2,1 = C(KCDE-Ran,KCDE).

. Furthermore, the p-value of the nonparametric test with a 95% confidence level is calculated and summarized for each scale (to analyze the significance of the difference between the two compared algorithms). The boxplots of all $CM$ values are illustrated in Figure 5a.

It can be seen from Table 4 and Figure 5a that for each instance, C(KCDE, KCDE-Ran) is larger than C(KCDE-Ran, KCDE) while all p-values are smaller than 0.05. Therefore, it can be concluded that most of the non-dominated solutions obtained by KCDE are better than those provided by KCDE-Ran and the hybrid initialization considerably enhances the quality of the initial population (and, thus, the performance of KCDE).

4.5. Effect of Knowledge Base

In order to examine the effectiveness of the knowledge base, the comparisons between KCDE and its variant KCDE-nKB were carried out. In KCDE-nKB, each offspring is generated with the information carried out by the corresponding selected parents. Each algorithm runs 10 times independently, the average $CM$ values of the instances with the same scale are calculated. Moreover, to analyze the significance of the difference between KCDE and KCDE-nKB, the sign test with a 95% confidence level was carried out, and the p-value of each scale is calculated and summarized in

Table 5. Average C Metric Values of KCDE and KCDE-nKB.

Dataset	${\mathit{C}}_{1,3}$1	${\mathit{C}}_{3,1}$2	p	Dataset	${\mathit{C}}_{1,3}$	${\mathit{C}}_{3,1}$	p	Dataset	${\mathit{C}}_{1,3}$	${\mathit{C}}_{3,1}$	p
f3n30+20s4	0.68	0.11	0.00	f4n30+20s4	0.68	0.12	0.00	f5n30+20s4	0.64	0.09	0.00
f3n30+20s5	0.64	0.15	0.00	f4n30+20s5	0.53	0.18	0.02	f5n30+20s5	0.70	0.07	0.00
f3n30+20s6	0.61	0.13	0.00	f4n30+20s6	0.65	0.11	0.00	f5n30+20s6	0.63	0.10	0.00
f3n30+30s4	0.67	0.10	0.00	f4n30+30s4	0.62	0.12	0.02	f5n30+30s4	0.70	0.06	0.00
f3n30+30s5	0.68	0.12	0.00	f4n30+30s5	0.67	0.09	0.00	f5n30+30s5	0.71	0.08	0.00
f3n30+30s6	0.65	0.09	0.00	f4n30+30s6	0.70	0.12	0.00	f5n30+30s6	0.69	0.12	0.00
f3n30+50s4	0.64	0.09	0.00	f4n30+50s4	0.65	0.11	0.02	f5n30+50s4	0.76	0.06	0.00
f3n30+50s5	0.61	0.19	0.02	f4n30+50s5	0.68	0.10	0.00	f5n30+50s5	0.67	0.09	0.00
f3n30+50s6	0.58	0.12	0.00	f4n30+50s6	0.59	0.11	0.00	f5n30+50s6	0.65	0.09	0.00
f3n50+20s4	0.55	0.11	0.00	f4n50+20s4	0.60	0.10	0.00	f5n50+20s4	0.59	0.09	0.00
f3n50+20s5	0.58	0.16	0.00	f4n50+20s5	0.55	0.15	0.00	f5n50+20s5	0.67	0.10	0.00
f3n50+20s6	0.61	0.09	0.00	f4n50+20s6	0.54	0.12	0.00	f5n50+20s6	0.67	0.10	0.00
f3n50+30s4	0.59	0.07	0.00	f4n50+30s4	0.68	0.09	0.00	f5n50+30s4	0.65	0.07	0.00
f3n50+30s5	0.66	0.09	0.00	f4n50+30s5	0.59	0.10	0.00	f5n50+30s5	0.66	0.11	0.00
f3n50+30s6	0.62	0.07	0.00	f4n50+30s6	0.59	0.12	0.02	f5n50+30s6	0.62	0.09	0.00
f3n50+50s4	0.65	0.07	0.00	f4n50+50s4	0.61	0.10	0.00	f5n50+50s4	0.61	0.10	0.00
f3n50+50s5	0.60	0.11	0.00	f4n50+50s5	0.62	0.12	0.00	f5n50+50s5	0.62	0.12	0.00
f3n50+50s6	0.66	0.09	0.00	f4n50+50s6	0.59	0.12	0.00	f5n50+50s6	0.56	0.13	0.00
f3n80+20s4	0.50	0.08	0.00	f4n80+20s4	0.63	0.12	0.00	f5n80+20s4	0.49	0.10	0.00
f3n80+20s5	0.59	0.08	0.00	f4n80+20s5	0.60	0.11	0.00	f5n80+20s5	0.57	0.14	0.02
f3n80+20s6	0.55	0.09	0.00	f4n80+20s6	0.60	0.12	0.02	f5n80+20s6	0.47	0.12	0.00
f3n80+30s4	0.58	0.06	0.00	f4n80+30s4	0.65	0.07	0.00	f5n80+30s4	0.52	0.11	0.00
f3n80+30s5	0.65	0.03	0.00	f4n80+30s5	0.61	0.06	0.00	f5n80+30s5	0.65	0.08	0.02
f3n80+30s6	0.59	0.05	0.00	f4n80+30s6	0.61	0.08	0.00	f5n80+30s6	0.58	0.10	0.00
f3n80+50s4	0.56	0.06	0.00	f4n80+50s4	0.64	0.05	0.00	f5n80+50s4	0.59	0.05	0.00
f3n80+50s5	0.58	0.07	0.00	f4n80+50s5	0.67	0.04	0.00	f5n80+50s5	0.64	0.08	0.00
f3n80+50s6	0.58	0.08	0.00	f4n80+50s6	0.63	0.07	0.00	f5n80+50s6	0.61	0.09	0.00

¹C_1,3 = C(KCDE,KCDE-nKB); ² C_3,1 = C(KCDE-nKB,KCDE).

. The boxplots of all $CM$ values are provided in Figure 5b.

From Table 5 and Figure 5b, it is obvious that C(KCDE,KCDE-nKB) is larger than C(KCDE-nKB,KCDE) for each instance and all p-values are lower than 0.05, which demonstrates that most solutions in the Pareto set approximation obtained by KCDE are better than those yielded by KCDE-nKB. Thus, it can be concluded that the utilization of the knowledge base can efficiently support KCDE by providing information carried out by elite individuals.

4.6. Effect of Local Intensification

The effectiveness of local intensification is investigated by comparing the proposed KCDE and KCDE-nLI, in which the local intensification part is absent. The average C(KCDE,KCDE-nLI) and C(KCDE-nLI, KCDE) values of the instances with the same scales were calculated and summarized in

Table 6. Average C Metric Values of KCDE and KCDE-nLI.

Dataset	${\mathit{C}}_{1,4}$1	${\mathit{C}}_{4,1}$2	p	Dataset	${\mathit{C}}_{1,4}$	${\mathit{C}}_{4,1}$	p	Dataset	${\mathit{C}}_{1,4}$	${\mathit{C}}_{4,1}$	p
f3n30+20s4	0.67	0.03	0.00	f4n30+20s4	0.73	0.02	0.00	f5n30+20s4	0.66	0.03	0.00
f3n30+20s5	0.55	0.03	0.00	f4n30+20s5	0.64	0.03	0.02	f5n30+20s5	0.67	0.03	0.00
f3n30+20s6	0.55	0.05	0.00	f4n30+20s6	0.72	0.02	0.00	f5n30+20s6	0.63	0.03	0.00
f3n30+30s4	0.59	0.02	0.00	f4n30+30s4	0.75	0.02	0.00	f5n30+30s4	0.65	0.03	0.00
f3n30+30s5	0.56	0.03	0.00	f4n30+30s5	0.72	0.02	0.00	f5n30+30s5	0.71	0.01	0.00
f3n30+30s6	0.65	0.03	0.00	f4n30+30s6	0.78	0.01	0.00	f5n30+30s6	0.69	0.02	0.00
f3n30+50s4	0.65	0.03	0.00	f4n30+50s4	0.71	0.02	0.00	f5n30+50s4	0.80	0.01	0.00
f3n30+50s5	0.53	0.02	0.00	f4n30+50s5	0.79	0.01	0.00	f5n30+50s5	0.72	0.01	0.00
f3n30+50s6	0.62	0.03	0.00	f4n30+50s6	0.74	0.01	0.00	f5n30+50s6	0.79	0.01	0.00
f3n50+20s4	0.56	0.07	0.00	f4n50+20s4	0.69	0.04	0.00	f5n50+20s4	0.63	0.03	0.00
f3n50+20s5	0.51	0.07	0.02	f4n50+20s5	0.70	0.04	0.00	f5n50+20s5	0.61	0.05	0.00
f3n50+20s6	0.48	0.06	0.00	f4n50+20s6	0.58	0.05	0.00	f5n50+20s6	0.61	0.05	0.00
f3n50+30s4	0.63	0.05	0.00	f4n50+30s4	0.75	0.03	0.00	f5n50+30s4	0.71	0.02	0.00
f3n50+30s5	0.52	0.06	0.00	f4n50+30s5	0.68	0.03	0.00	f5n50+30s5	0.62	0.02	0.00
f3n50+30s6	0.52	0.05	0.00	f4n50+30s6	0.58	0.05	0.00	f5n50+30s6	0.66	0.03	0.00
f3n50+50s4	0.62	0.04	0.00	f4n50+50s4	0.78	0.03	0.00	f5n50+50s4	0.74	0.02	0.00
f3n50+50s5	0.51	0.04	0.00	f4n50+50s5	0.72	0.02	0.00	f5n50+50s5	0.67	0.01	0.00
f3n50+50s6	0.47	0.04	0.00	f4n50+50s6	0.53	0.04	0.00	f5n50+50s6	0.66	0.02	0.00
f3n80+20s4	0.55	0.09	0.02	f4n80+20s4	0.59	0.05	0.00	f5n80+20s4	0.60	0.04	0.00
f3n80+20s5	0.46	0.09	0.02	f4n80+20s5	0.56	0.05	0.02	f5n80+20s5	0.57	0.04	0.00
f3n80+20s6	0.48	0.08	0.00	f4n80+20s6	0.50	0.07	0.00	f5n80+20s6	0.58	0.04	0.00
f3n80+30s4	0.52	0.07	0.00	f4n80+30s4	0.66	0.04	0.00	f5n80+30s4	0.59	0.05	0.00
f3n80+30s5	0.44	0.08	0.02	f4n80+30s5	0.57	0.06	0.02	f5n80+30s5	0.52	0.04	0.02
f3n80+30s6	0.48	0.08	0.02	f4n80+30s6	0.42	0.06	0.02	f5n80+30s6	0.71	0.04	0.02
f3n80+50s4	0.49	0.05	0.00	f4n80+50s4	0.71	0.04	0.00	f5n80+50s4	0.71	0.03	0.00
f3n80+50s5	0.43	0.06	0.02	f4n80+50s5	0.59	0.03	0.00	f5n80+50s5	0.60	0.02	0.00
f3n80+50s6	0.46	0.05	0.00	f4n80+50s6	0.44	0.05	0.00	f5n80+50s6	0.70	0.02	0.00

¹C_1,4 = C(KCDE,KCDE-nLI); ² C_4,1 = C(KCDE-nLI,KCDE).

. Furthermore, to analyze the significance of the difference between the two algorithms, the sign test with the 95% confidence level was carried out and the p-values were calculated. The boxplots of all $CM$ values are provided in Figure 5c.

As shown in Table 6 and Figure 5c, for each instance, C(KCDE,KCDE-nLI) is larger than C(KCDE-nLI, KCDE). Moreover, all of the p-values are lower than 0.05; it can be concluded that most elite solutions yielded by KCDE can dominate those obtained by KCDE-nLI. Thus, it can be concluded that the proposed local intensification operator has significant influence on enhancing the performance of KCDE.

4.7. Comparisons to Other Algorithms

To our knowledge, there is currently no research on EDHFRP. Therefore, two state-of-the-art algorithms designed for problems similar to EDHFRP, i.e., DCMA [11] and BCMA [7], along with two classical multi-objective optimization algorithms, namely NSGA-II [49] and MOEA/D [50], were adapted for comparison. The parameters of DCMA and BCMA were set as suggestions in [7,11], respectively. The stopping criterion and running times are set as mentioned before.

The performances of all of the algorithms on instances of different scales are first evaluated by the C metric, which describes the dominant relationship between Pareto set approximations obtained by different algorithms. The average $CM$ between KCDE and the other 4 algorithms in 10 instances of each scale are indicated in

Table 7. Average C Metric Values of KCDE and NSGA-II, MOEA/D, BCMA, DCMA.

Dataset	${\mathit{C}}_{1,5}$1	${\mathit{C}}_{5,1}$2	p	${\mathit{C}}_{1,6}$3	${\mathit{C}}_{6,1}$4	p	${\mathit{C}}_{1,7}$5	${\mathit{C}}_{7,1}$6	p	${\mathit{C}}_{1,8}$7	${\mathit{C}}_{8,1}$8	p
f3n30+20	0.735	0.018	0.00	0.796	0.014	0.00	0.517	0.013	0.00	0.705	0.000	0.00
f3n30+30	0.763	0.013	0.00	0.792	0.011	0.00	0.500	0.010	0.00	0.595	0.000	0.00
f3n30+50	0.731	0.012	0.00	0.752	0.012	0.00	0.489	0.008	0.00	0.381	0.000	0.00
f3n50+20	0.676	0.034	0.00	0.768	0.022	0.00	0.521	0.022	0.00	0.898	0.000	0.00
f3n50+30	0.749	0.021	0.00	0.792	0.019	0.00	0.561	0.018	0.00	0.820	0.000	0.00
f3n50+50	0.751	0.014	0.00	0.814	0.014	0.00	0.595	0.009	0.00	0.705	0.000	0.00
f3n80+20	0.674	0.035	0.01	0.710	0.032	0.00	0.446	0.036	0.01	0.952	0.000	0.00
f3n80+30	0.675	0.031	0.00	0.731	0.027	0.00	0.458	0.026	0.00	0.933	0.000	0.00
f3n80+50	0.740	0.024	0.01	0.768	0.020	0.00	0.547	0.018	0.01	0.903	0.000	0.00
f4n30+20	0.784	0.011	0.00	0.842	0.011	0.00	0.499	0.008	0.00	0.713	0.000	0.00
f4n30+30	0.827	0.008	0.00	0.851	0.005	0.00	0.532	0.006	0.00	0.634	0.000	0.00
f4n30+50	0.848	0.004	0.00	0.869	0.003	0.00	0.468	0.004	0.00	0.549	0.000	0.00
f4n50+20	0.770	0.020	0.00	0.821	0.013	0.00	0.567	0.013	0.00	0.860	0.000	0.00
f4n50+30	0.817	0.011	0.00	0.860	0.012	0.00	0.569	0.011	0.00	0.773	0.000	0.00
f4n50+50	0.847	0.008	0.00	0.890	0.009	0.00	0.578	0.009	0.00	0.653	0.000	0.00
f4n80+20	0.735	0.026	0.00	0.782	0.015	0.00	0.508	0.018	0.01	0.935	0.000	0.00
f4n80+30	0.772	0.020	0.00	0.813	0.019	0.00	0.563	0.012	0.00	0.895	0.000	0.00
f4n80+50	0.861	0.010	0.00	0.883	0.010	0.00	0.594	0.009	0.00	0.852	0.000	0.00
f5n30+20	0.769	0.013	0.00	0.802	0.010	0.00	0.519	0.006	0.00	0.645	0.003	0.00
f5n30+30	0.838	0.006	0.00	0.864	0.004	0.00	0.574	0.004	0.00	0.480	0.001	0.00
f5n30+50	0.831	0.005	0.00	0.863	0.003	0.00	0.523	0.003	0.00	0.656	0.001	0.00
f5n50+20	0.754	0.011	0.00	0.858	0.006	0.00	0.571	0.007	0.00	0.855	0.000	0.00
f5n50+30	0.818	0.010	0.00	0.910	0.008	0.00	0.643	0.006	0.00	0.764	0.001	0.00
f5n50+50	0.822	0.007	0.00	0.914	0.004	0.00	0.653	0.005	0.00	0.828	0.000	0.00
f5n80+20	0.703	0.027	0.01	0.824	0.018	0.00	0.521	0.020	0.01	0.946	0.002	0.00
f5n80+30	0.763	0.017	0.01	0.869	0.010	0.00	0.559	0.011	0.00	0.917	0.001	0.00
f5n80+50	0.823	0.011	0.00	0.913	0.006	0.00	0.562	0.013	0.00	0.798	0.000	0.00

¹C_1,5 = C(KCDE,NSGA-II); ² C_5,1 = C(NSGA-II,KCDE); ³ C_1,6 = C(KCDE,MOEA/D) ; ⁴ C_6,1 = C(MOEA/D,KCDE); ⁵ C_1,7 = C(KCDE,BCMA); ⁶ C_7,1 = C(MOEA/D,BCMA); ⁷ C_1,8 = C(KCDE,DCMA); ⁸ C_8,1 = C(MOEA/D,DCMA).

, while Figure 6 represents the boxplots of all $CM$. Moreover, the nonparametric test with a confidence level of 95% is carried out to show if there are significant differences between different algorithms; the p values are also recorded in Table 7. It is noticeable that C(KCDE, NSGA-II), C(KCDE, MOEA/D), C(KCDE, BCMA), and C(KCDE, DCMA) are much larger than C(NSGA-II, KCDE), C(MOEA/D, KCDE), C(BCMA, KCDE), and C(DCMA, KCDE) in all instances of different scales, especially large-scale instances. Moreover, the p values are all lower than 0.05, which indicates that most of the solutions obtained by KCDE can dominate those of the other algorithms. Therefore, it can be demonstrated that KCDE performs significantly better than NSGA-II, MOEA/D, BCMA, and DCMA.

The average HV and GD values of all algorithms are provided in

Table 8. HVs and GDs of KCDE and NSGA-II, MOEA/D, BCMA, DCMA.

Dataset	HV					GD
	KCDE	NSGA	MOEA/D	BCMA	DCMA	KCDE	NSGA	MOEA/D	BCMA	DCMA
f3n30+20	0.85	0.10	0.10	0.08	0.04	0.13	1.06	1.02	1.00	1.40
f3n30+30	0.86	0.11	0.11	0.09	0.05	0.10	1.02	0.96	0.99	1.35
f3n30+50	0.86	0.12	0.11	0.10	0.06	0.09	0.95	0.91	1.00	1.26
f3n50+20	0.86	0.08	0.07	0.07	0.03	0.15	1.11	1.04	1.03	1.50
f3n50+30	0.88	0.09	0.09	0.07	0.03	0.11	1.06	0.99	1.03	1.45
f3n50+50	0.88	0.10	0.10	0.08	0.05	0.07	0.98	0.92	1.01	1.39
f3n80+20	0.87	0.06	0.06	0.05	0.02	0.18	1.16	1.08	1.07	1.56
f3n80+30	0.87	0.06	0.07	0.06	0.02	0.15	1.14	1.03	1.05	1.54
f3n80+50	0.89	0.08	0.08	0.07	0.03	0.09	1.02	0.98	1.03	1.51
f4n30+20	0.87	0.10	0.10	0.07	0.04	0.10	1.03	0.99	1.07	1.27
f4n30+30	0.90	0.11	0.11	0.08	0.05	0.07	0.98	0.95	1.04	1.22
f4n30+50	0.91	0.13	0.12	0.10	0.06	0.06	0.95	0.90	1.01	1.16
f4n50+20	0.88	0.07	0.07	0.05	0.02	0.12	1.09	1.02	1.13	1.46
f4n50+30	0.90	0.09	0.08	0.06	0.03	0.09	1.06	0.97	1.12	1.37
f4n50+50	0.90	0.10	0.10	0.08	0.05	0.06	0.99	0.92	1.09	1.27
f4n80+20	0.87	0.06	0.06	0.05	0.02	0.16	1.14	1.07	1.14	1.51
f4n80+30	0.89	0.06	0.06	0.05	0.02	0.12	1.09	1.02	1.14	1.47
f4n80+50	0.91	0.08	0.08	0.06	0.03	0.08	1.03	0.97	1.14	1.41
f5n30+20	0.84	0.09	0.09	0.06	0.05	0.17	1.09	1.05	1.07	1.36
f5n30+30	0.90	0.11	0.11	0.08	0.05	0.08	1.03	0.97	1.05	1.32
f5n30+50	0.91	0.13	0.12	0.09	0.04	0.07	0.97	0.92	1.04	1.40
f5n50+20	0.87	0.07	0.07	0.05	0.03	0.15	1.12	1.08	1.15	1.47
f5n50+30	0.91	0.08	0.08	0.06	0.04	0.08	1.06	1.01	1.13	1.42
f5n50+50	0.91	0.10	0.09	0.07	0.03	0.05	1.00	0.96	1.11	1.49
f5n80+20	0.85	0.05	0.05	0.04	0.02	0.20	1.20	1.14	1.17	1.54
f5n80+30	0.89	0.06	0.06	0.04	0.03	0.11	1.12	1.08	1.15	1.51
f5n80+50	0.91	0.07	0.07	0.05	0.04	0.07	1.05	1.00	1.12	1.36

. Figure 7 and Figure 8 illustrate the boxplots of two evaluation indicators on instances of different scales and the distribution of Pareto set approximations respectively. For each instance, the objective values are normalized with the largest value and the lowest value obtained by all the algorithms as higher bounds and lower bounds, respectively.

From Table 8 and Figure 7, one can see that the KCDE results in higher HV and lower GD in all instances. Thus, KCDE performs better than NSGA-II, MOEA/D, BCMA, and DCMA in solving EDHFRP since higher HV values indicate better convergence and diversity of the Pareto set approximations while lower GD values represent greater convergence properties. The Pareto set approximations of five algorithms in instances f3n30+30s5 and f3n80+20s5 are displayed in Figure 8, from which one can observe that the non-dominated solutions achieved by KCDE distribute more evenly than those obtained by the other four algorithms. As a result, it can be concluded that KCDE is more able to provide better solutions with shorter makespans and lower total energy consumption for EDHFRP.

5. Conclusions

This paper focuses on solving the energy-efficient distributed hybrid flow-shop rescheduling problem, minimizing the makespans of original and newly arrived orders as well as the total energy consumption, simultaneously. To this end, a knowledge-based cooperative differential evolution algorithm was designed. A group of broad instances with stochastic arrival times for the new order were generated. The effects of major parameters on the proposed algorithm were investigated and the importance of each problem-specific strategy was verified through extensive tests and experiments. The comprehensive comparisons and analysis demonstrate that KCDE is able to perform in a superior way in each instance.

In the future, the proposed algorithm can be extended to solve distributed rescheduling problems with more complicated constraints, such as fuzzy processing time, sequence-dependent setup time, and limited buffers. Moreover, this work can be used to solve scheduling problems in other kinds of environments, including no-idle, no-wait, lot-streaming, and heterogeneous factories, by simply changing the model. Moreover, optimization objectives with more practical significance, including tardiness, time-of-use electricity costs, and transportation can be considered in the algorithm.