A Novel Multi-Population Artificial Bee Colony Algorithm for Energy-Efficient Hybrid Flow Shop Scheduling Problem

Considering green scheduling and sustainable manufacturing, the energy-efficient hybrid flow shop scheduling problem (EHFSP) with a variable speed constraint is investigated, and a novel multi-population artificial bee colony algorithm (MPABC) is developed to minimize makespan, total tardiness and total energy consumption (TEC), simultaneously. It is necessary for manufacturers to fully understand the notion of symmetry in balancing economic and environmental indicators. To improve the search efficiency, the population was randomly categorized into a number of subpopulations, then several groups were constructed based on the quality of subpopulations. A different search strategy was executed in each group to maintain the population diversity. The historical optimization data were also used to enhance the quality of solutions. Finally, extensive experiments were conducted. The results demonstrate that MPABC can achieve an outstanding performance on three metrics DIR, c and nd for the considered EHFSP.


Introduction
Shop scheduling is an essential subject and an effective way to improve resource utilization [1]. With the increase of energy demand and the emergence of various environmental issues, reducing energy consumption has become particularly urgent for the manufacturing industry. Many scholars devote themselves to the research of energyefficient scheduling problems, related to single machine [2], unrelated parallel machines [3], job shop [4], permutation flow shop [5] and hybrid flow shop [6].
In the typical flow shop scheduling problem, sequencing jobs are processed in predetermined orders, and each stage contains only one machine, which has been proved to be NP-hard [7]. Regarding the energy-efficient flow shop scheduling problem, Ding et al. [8] investigated a carbon-efficient scheduling of flow shops. Foumani et al. [9] considered the impact of various carbon reduction policies. Mansouri et al. [10] addressed a two-machine green flow shop scheduling. Jiang et al. [11] studied an energy-efficient permutation flow shop scheduling problem.
As an extension of flow shop, hybrid flow shop employs multiple parallel machines to enlarge the capacity of stage and eliminate the restriction of the bottleneck stage [12]. The hybrid flow shop scheduling problem (HFSP) extensively exists in various real industrial scenarios, such as steel [13], textile [14], glass [15], paper [16] and electronics [17].
EHFSP often consists of green constraints, green objectives and three sub-problems including job permutation, machine assignment and speed selection, which is more complicated than energy-efficient flow shop scheduling problems and apparently NP-hard. In recent years, EHFSP has attracted much attention from researchers and manufacturers. Bruzzone et al. [18] utilized a heuristic method to obtain an energy-aware scheduling in a sustainable manufacturing process. Dai et al. [19] designed a novel genetic-simulated annealing algorithm to compromise between makespan and TEC, where makespan is the total time length of the schedule, and TEC is the sum of energy consumption in the production process. Meng et al. [20] developed an improved genetic algorithm to solve the energy-related problem with unrelated parallel machines. Gong et al. [21] incorporated a hybrid evolutionary algorithm into an EHFSP reckoning in the worker flexibility. Wu et al. [22] applied a hybrid non-dominated sorting genetic algorithm (NSGA-II) with variable local search to EHFSP considering renewable energy. Afterwards, the optimization algorithms based on hybrid NSGA-II were also raised for EHFSP in Zeng et al. [23] and Liu et al. [24]. To investigate EHFSP in a real-life case, a multi-level optimization method [25] and problem-tailored constructive heuristic method [26] were proposed. Work conducted by Zhou and Liu [27] investigated a multi-objective algorithm based on the Nawaz-Enscore-Ham heuristic for EHFSP with fuzzy processing time. Schulz et al. [28] designed a novel iterated local search algorithm to consider peak load, total energy costs and makespan in energy-aware scheduling.
The previous work also considered EHFSP through meta-heuristic optimization algorithms. Luo et al. [29] developed an improved ant colony algorithm for HFSP considering the time-of-use electricity prices. A teaching−learning-based optimization algorithm was provided by Lin et al. [30] to minimize carbon footprint and makespan. Zhou et al. [31] utilized an imperialist competitive algorithm for minimizing makespan and TEC simultaneously. Afterwards, Zhang et al. [32] suggested a three-stage optimization framework based on decomposition to minimize makespan and TEC in EHFSP. Li et al. [33] constructed another version of MOEA/D in the lot-streaming hybrid flow shop scheduling problem. Besides, other meta-heuristic optimization algorithms were also applied in EHFSP, including particle swarm optimization algorithm [34], shuffled frog-leaping algorithm [35] and artificial fish swarm algorithm [36].
In the literature above, speed constraint in EHFSP is not considered. At present, machines are able to operate at variable processing speeds. In general, a higher speed level leads to less processing time and more energy consumption. The speed selection subproblem in EHFSP is used to choose a proper speed for the selected processing machine, which can establish a balance between productivity and energy consumption. Meanwhile, there is an available speed set for the machines of each stage, and the machine speed needs to be regarded as an independent decision variable. Three sub-problems in EHFSP should be optimized simultaneously. Therefore, it is necessary to investigate the speed constraint in EHFSP.
In recent years, only several works have been studied on EHFSP with a variable speed constraint. Liu et al. [37] introduced the speed selection into a special two-stage EHFSP in the iron and steel industry. Tang et al. [38] studied an improved particle swarm optimization method to minimize makespan and TEC. Zhang et al. [39] suggested a discrete ABC based on decomposition to address the green scheduling problem with sequence dependent setup operations. A multi-objective genetic algorithm was developed for EHFSP with lot streaming in Chen et al. [40]. Lei et al. [41] proposed a teachinglearning-based optimization algorithm for the bi-objective EHFSP. Li et al. [42] presented an imperialist competitive algorithm with the relative importance of objectives. Afterwards, Oztop et al. [43] suggested an ensemble of metaheuristics to optimize makespan and energy consumption simultaneously. In conclusion, research on EHFSP with a variable speed constraint is significant in theory and application level.
Artificial bee colony (ABC) is a promising metaheuristic algorithm motivated by the collective behavior of the bee colony [44]. Since ABC was proposed, it was mainly used to solve unconstraint and constraint continuous function optimization problems [45,46]. Compared with other metaheuristic algorithms, there are some remarkable characteristics in ABC such as its effective global exploration capability, excellent local exploitation ability and fast convergence speed. Due to its simple structure and convenient implementation, it has also been applied to solve combinatorial optimization problems, especially in the field of production scheduling problems [47][48][49][50][51][52][53][54]. Its successful applications on scheduling prove that ABC possesses an excellent search ability and convergence rate. In addition, ABC is seldom adopted in EHFSP. These factors motivated our idea to use ABC to solve EHFSP. However, ABC cannot be directly used to solve EHFSP because it is a typical discrete optimization problem and EHFSP is completely different from other scheduling problems. The algorithm needs to be modified according to the features of EHFSP. Therefore, we are committed to make improvements on ABC for solving EHFSP in this paper.
The contributions of this study can be summarized as follows: (1) EHFSP with a variable speed constraint is addressed to optimize makespan, total tardiness and TEC simultaneously, where total tardiness is a vital index for manufactures and equals the sum of tardiness of all jobs. (2) A novel multi-population artificial bee colony algorithm is proposed, in which a completely different multi-population division method is adopted. After randomly dividing the population into a few subpopulations, several groups are constructed. Each group implements different search approaches. (3) The quality of subpopulation is newly defined and mainly applied in group division. (4) The historical optimization data were stored in a memory set and fully utilized to update the weakest group.
The remainder of this paper is organized as follows. The considered problem is described in Section 2. The introduction of the original ABC is presented in Section 3, followed by the description of MPABC for EHFSP in Section 4. Extensive numerical experiments are revealed in Section 5. Conclusions and future work are provided in Section 6.

Problem Description
The notations used in this paper are shown in Table 1.

Notation Description n
The number of jobs m The number of stages s The number of available unrelated parallel machines for S j d The number of available processing speeds for M jk M jk The kth machine of stage j η ijk The basic processing requirement of job J i on machine M jk t ijkl The processing time of job J i on machine M jk at speed v l z jk (t) z jk (t) = 1, if the machine M jk is in stand-by mode at time t; otherwise, z jk (t) = 0 y jkl (t) y jkl (t) = 1, if the stage j is processed on kth machine with speed v l at time t; otherwise, y jkl (t) = 0 y ijkl y ijkl = 1, if the jth stage in job i is processed on kth machine with speed v l otherwise, y ijkl = 0 The starting processing time of job i at stage j b ij The finishing processing time of job i at stage j E jkl The energy consumption per unit time when machine M jk runs at speed v l SE jk The energy consumption per unit time when machine M jk runs in stand-by mode For M jk , variable processing speeds V = {v l } 1≤l≤d are available. Moreover, if a job is processed at one stage, an assigned machine and a selected speed must be arranged for its processing. The arranged machine and selected speed cannot be changed during the whole processing. When job J i is processed on machine M jk at speed v l , the processing time t ijkl is denoted as η ijk /v l . Several assumptions are presented as follows: • Buffer size between stages is unlimited; • Each machine can only process one operation; • Some stages can be skipped, but each job must be processed at one stage at least; • Once the process starts, interruption and preemption are not allowed.
In EHFSP, economic benefit and environmental impact are symmetrical and equally important. In this paper, the aim is to arrange jobs to the corresponding machines of each stage, decide the job sequences on each machine and choose the suitable speed for machines of each stage in order to minimize the following three objectives simultaneously when all constraints are met. min Constraints: x ii jk +x i ijk ≤ 1, ∀i, i ∈ J, j ∈ O, k ∈ S j (7) where Equation (1) is makespan, Equation (2) indicates the total tardiness (TT), and Equation (3) is total energy consumption. Constraint (4) illustrates that only one machine with a single speed level is available for each stage of jobs at a time. Constraint (5) ensures that interruption is not allowed. Constraints (6) and (7) show the sequential constraints of job processing. Constraints (8) and (9) describe the range of variables. There are some definitions in the multi-objective optimization problem with G objectives. Definition 1. Pareto dominance. x Pareto dominates y, denoted as x y, if f i (x) ≤ f i (y) for ∀i ∈ {1, 2, . . . , G}, and f j (x) < f j (y), ∃j ∈ {1, 2, . . . , G} Definition 2. Pareto optimal. Solution x * ∈ Θ is a Pareto optimal solution if and only if x ∈ Θ such that x x * . Solution x * is also called a non-dominated solution.
Definition 3. Pareto front (PF). Pareto optimal set (PS) is a set of Pareto optimal solutions. PF is a set of Pareto optimal vectors, PF ={F(x) ∈ R m |x ∈ PS }. Table 2 gives an example, including ten jobs with three stages, two machines for each stage and three processing speed levels available for each machine.

Introduction to ABC
Bionic computing is a comprehensive interdisciplinary subject in which problems are solved by discovering mechanisms in biology [55][56][57]. At present, few existing bionic computing methods can be used to simulate biological characteristics. Bees are intelligent insects. From bee dancing and optimizing honeycomb design to flying, human beings have gained a lot of enlightenment from these features. Karl von Frisch won the Nobel Prize in 1973 for his research on bee behavior. ABC is a potential bionic computing method inspired by their biological characteristics [44].
In ABC, food sources represent solutions in search space, and the nectar density corresponds to the fitness value of a feasible solution. Three main phases are included in ABC: employed bee phase, onlooker bee phase and scout bee phase. Each bee colony plays a different role in the optimization process. Employed bees fly to exploit food source information around the assigned solutions. Onlooker bees wait and choose a selected solution from the employed bees. Scout bees help to jump out of the local optima.
The basic framework of ABC is shown in Algorithm 1.

Algorithm 1
The basic framework of ABC.
scout bee phase 7. end while 8. output best solutions found so far 9. end In initialization, related parameters are set and initial population P with N solutions is randomly generated.
In the employed bee phase, food sources x m are all allocated to employed bees. A new source y m is obtained by where x n is another randomly selected solution from P, φ ∈ [−1, 1]. Greedy selection based on fitness value is conducted to determine whether x m can be updated by y m . Afterwards, employed bees dance and share information with onlooker bees. In the onlooker bee phase, the food source is chosen through the roulette selection probability where p m is the probability for x m . A new source y m is acquired and compared with x m in a similar way to that in the employed bee phase. If a solution has not been updated after Limit trails, the corresponding employed bee will become a scout bee and exploit a new source to replace the old one, where Limit is the maximum number of attempts.

MPABC for EHFSP
In previous ABCs [45][46][47][48][49][50][51][52][53], the whole population is seldom categorized into several groups, and the food sources of onlooker bees are completely from employed bees. In this study, a novel multi-population method is proposed, where the divided groups are different from the three kinds of bee colonies in the classical ABC algorithm. When the population is divided into several groups, ABC can be considered as a multi-population metaheuristic algorithm. Multi-population division is useful to maintain the population diversity and avoid premature convergence [58]; on the other hand, different search strategies can be implemented in each group simultaneously, which significantly contributes to improving the search efficiency. However, these approaches are seldom introduced into ABC. In this paper, a novel MPABC is proposed based on above ideas.
As the inverse procedure of encoding, decoding translates the solution into a feasible schedule. Generally, decoding starts from the first element of JS. For machine π i , the corresponding machine q π i j and processing speed z π i j are selected according to MS and SS. Subsequent stages are arranged in sequence until all processes are completed.

Initialization and Group Division
Initial population P with N solutions is stochastically genera

Initialization and Group Division
Initial population P with N solutions is stochastically generated. For each solution i, the fitness value F i is calculated as follows [42].
where rank i is a non-dominated sorting value [59], dist i is the crowding distance defined in [60] and θ rank i represents a set of solutions with rank i . After that, N/s solutions are randomly allocated into s subpopulations. For each subpopulation subpop j (j = 1, 2 . . . s), the quality Q j is measured as follows.
where TF j is the total fitness value of subpop j , Num Iter j represents the update quantity of solutions for subpop j in each iteration, Q j expresses the quality of subpop j , δ denotes the weight coefficient and ε is a real number. The larger Q j is, the better the subpopulation j is. δ is set to 0.6, and ε is set to 0.001 by experiments.
Unlike the previous division methods [61][62][63][64], a new group division method according to the quality of the subpopulation is described in this paper. When the population is randomly divided into s subpopulations, several groups are constructed based on Q j . The detailed division process is as follows. Sort Q j (j = 1, 2, · · ·s) in descending order, suppose that is assigned into Group 1 , the last subpopulation subpop s is allocated to Group 3 and the remaining s − p−1 subpopulations are distributed to Group 2 .Q represents the value of average quality.
Members of Group 1 are assigned to act as food sources of employed bees, solutions in Group 2 are used for sources of onlooker bees, and individuals in Group 3 are not directly involved in evolutionary operations and wait to be updated by historical optimization data.

Employed Bee Phase
Employed bees often search for new sources through crossover. In this phase, we propose a new method to execute an employed bee search. p subpopulations in Group 1 are totally regarded as the food sources of employed bees. These employed bees obtain computing resources through competition, that is, a higher quality solution can compete for more computing resources and opportunities to conduct global and multi-neighborhood searches. It is important to accelerate the search process and improve search efficiency.
The search procedures are shown in Algorithm 2.
Algorithm 2 Employed bee phase.
stochastically select a solution y ∈ subpop k (j = k) from Group 1 5.
randomly choose a solution z ∈ Ω 6. if x y then 7.
conduct global search between x and z 8.
execute multi-neighborhood search for x 9. else 10.
perform global search between y and z 11.
apply multi-neighborhood search to y 12.
end if 13. end for 14. end for 15. else 16. for each x ∈ subpop j do 17.
stochastically select y ∈ subpop k (j = k, x = y) from Group 1 , z ∈ Ω 18. execute above search strategies for the higher quality solution between x and y 19. end for 20. end if Global search is implemented based on three crossover operations [41,42], which is described as follows. For two solutions, if a random number θ < α, the crossover operation on the scheduling string is carried out; otherwise, another two crossover operations on the machine assignment string and speed selection string are performed with the same probability.
In general, the scheduling sub-problem is much more complicated than the machine assignment and speed selection sub-problems, so α > 0.5 is set to guarantee that more computing resources can be allocated to the scheduling sub-problem. In this paper, α is equal to 0.6 in experiments.
Three neighborhood structures insert, change and speed [41,42] are adopted, which are donated as N 1 , N 2 and N 3 , respectively. For solution λ, N 1 is conducted on the scheduling string. A new scheduling plan can be generated by randomly selecting a job and inserting it into a random position.
N 2 is executed on the machine assignment string and used to choose a new processing machine. The detailed steps are as follows: A machine set Ψ = {q ik ||S k | > 1, 1 ≤ i ≤ n, 1 ≤ k ≤ m} is firstly built, then randomly choose some elements from Ψ. If q ik is selected, a randomly selected machine from S k is used to replace q ik . N 3 is held on speed selection string and used for a new speed choice for the given machines. A new solution is obtained by randomly selecting some elements from the speed selection string and distributing a new speed to substitute for the selected elements stochastically. An example of three neighborhood structures is shown in Figure 2.
is firstly built, then randomly choose some elements from  . If ik q is selected, a randomly selected machine from k S is used to replace ik q .

3
N is held on speed selection string and used for a new speed choice for the given machines. A new solution is obtained by randomly selecting some elements from the speed selection string and distributing a new speed to substitute for the selected elements stochastically. An example of three neighborhood structures is shown in Figure 2. Memory set  is used to store historical optimization data, and initial  is empty. In the evolution procedure, the dominated solutions can be regarded as historical optimization data. When a solution g is used to renew  , it is directly included in  . If the number of solutions in  is greater than max  , the worst solution is removed.
is set to 20 after many experiments.
generate a new solution if g is non-dominated by  then 5.
replace  with g and updated  with g 6.
updated  with g , Job5 Multi-neighborhood search is described in Algorithm 3, where N t (λ) is a neighborhood set after the N t neighborhood search for X i , and R is the number of neighborhood searches.
The initial archive Ω consists of non-dominated solutions generated by MPABC. If a new solution g is used to update Ω, it is firstly added to Ω; then a comparison according to the Pareto dominance relationship is conducted, and all the dominated solutions are deleted.
Memory set Λ is used to store historical optimization data, and initial Λ is empty. In the evolution procedure, the dominated solutions can be regarded as historical optimization data. When a solution g is used to renew Λ, it is directly included in Λ. If the number of solutions in Λ is greater than |Λ| max , the worst solution is removed. |Λ| max is set to 20 after many experiments.
generate a new solution g ∈ N t (λ) 4.
if g is non-dominated by λ then 5.
replace λ with g and updated Ω with g 6.
end if 10. end for

Onlooker Bee Phase
Generally, sources of onlooker bees are all from employed bees. In this phase, solutions in Group 2 belong to the food sources of onlooker bees. As the elite solutions in Group 1 have better fitness than members in Group 1 , best solutions from Group 1 are chosen as the following object of most onlooker bees, which is important to generate better scheduling plans. Meanwhile, a small number of onlooker bees execute multi-neighborhood search. Diversified search strategies performed on onlooker bees are beneficial to balance exploitation and exploration.
The search process in the onlooker bee phase is shown in Algorithm 4, where µ and β are real numbers, which are set to 0.1 and 0.3 in experiments.
In the above mentioned search process, onlooker bees can own the exclusive food sources, randomly select a solution from the upper employed bees to follow and employ the multiple search strategies. Therefore, the activities of onlooker bees differ from other ABCs.
1. select µ * p * N s best solutions from Group 1 and construct a set Θ 2. for each x ∈ Group 2 do 3.
if a random number η < β then 4.
perform multi-neighborhood search on x 5. else 6.
randomly choose a solution y from set Θ 7.
execute global search between x and y 8.
end if 9. end for

Scout Bee Phase
In the original ABC, if a solution x remains unchanged after Limit trials, it will be replaced by a randomly produced solution x . In general, x is of lower quality, that is, it is probably an inefficient scheduling scheme. In this study, a solution y ∈ Λ was chosen to implement multi-neighborhood search and if a new solution z is not dominated by x, solution z replaces x.
For members in Group 3 , the worst quality is owned, and the search efficiency is lower than other ones. In this paper, the historical optimization data stored in memory set Λ are used to update Group 3 . Firstly, the solutions in Λ and Group 3 are combined, then Group 3 solutions with the best solutions are selected to replace all the ones in Group 3 . In this way, Group 3 is updated and the population quality is also enhanced.
After all groups are updated, they are regrouped according to Q j . The division process is the same as that described in Section 4.2.

Algorithm Description
The whole procedure of MPABC is shown in Algorithm 5, where the terminal criterion is the maximum elapsed time.
while the terminal criterion is not met do 4.
for each solution x ∈ Group 1 do 5.
for each solution y ∈ Group 2 do 8.
scout bee phase 11. divide all subpopulations into three groups 12. end while 13. output solutions in Ω The flow chart of MPABC is shown in Figure 3. divide all subpopulations into three groups 12.
end while 13. output solutions in  The flow chart of MPABC is shown in Figure 3. MPABC has the following characteristics. The whole population is divided into some subpopulations, and several groups are constructed according to the quality of subpopulations. Different groups execute different search strategies. The historical optimization data are used to update the worst group. Global search and multineighborhood searches are also kept balanced. These characteristics are beneficial to maintain population diversity and improve search efficiency in searching for an optimal scheduling scheme taking into account makespan, total tardiness and TEC.

Experiments and Results
Extensive computational experiments were implemented to verify the performance of MPABC. All experiments were coded in Visual Studio 2015 and run in the same environment: Intel(R) Core(TM) i5-6200U CPU @ 2.30GHz 2.40GHz, 8.00GB RAM.

Instances, Metrics and Comparative Algorithms
In this paper, we used the instances in the literature [41,42]  . Other detailed descriptions such as j S and V are described in [41]. Three performance metrics were utilized to evaluate the results of different algorithms. MPABC has the following characteristics. The whole population is divided into some subpopulations, and several groups are constructed according to the quality of subpopulations. Different groups execute different search strategies. The historical optimization data are used to update the worst group. Global search and multi-neighborhood searches are also kept balanced. These characteristics are beneficial to maintain population diversity and improve search efficiency in searching for an optimal scheduling scheme taking into account makespan, total tardiness and TEC.

Experiments and Results
Extensive computational experiments were implemented to verify the performance of MPABC. All experiments were coded in Visual Studio 2015 and run in the same environment: Intel(R) Core(TM) i5-6200U CPU @ 2.30 GHz 2.40 GHz, 8.00 GB RAM.

Instances, Metrics and Comparative Algorithms
In this paper, we used the instances in the literature [41,42] to test the performance of MPABC on EHFSP with a variable speed constraint. A total of 44 instances were adopted, consisting of jobs n ∈ {20, 30,40,50,60,70,80,90, 100, 110, 120} and stages m ∈ {2, 4, 6, 8}. Other detailed descriptions such as S j and V are described in [41]. Three performance metrics were utilized to evaluate the results of different algorithms.
Metric DI R [65] is often applied to measure the convergence performance, in which the distance of non-dominated set Ω l relative to a reference set Ω * is computed.
where d xy represents the distance between a non-dominated solution y and a reference point x in the normalized objective space. Ω l is an approximation set to PF, Ω * is a reference set that consists of the whole non-dominated solutions gained by each algorithm. The smaller the DI R value is, the better the convergence is.
Metric c [66] is often adopted to compare approximate Pareto optimal set acquired by algorithms. c(A, B) indicates that the proportion of individuals in B are dominated by members in A.
c(A, B) = |{y ∈ B: x ∈ A, x y}| |B| (16) Metric nd is usually used to denote the number of non-dominated solutions. A larger nd signifies more non-dominated solutions in the approximate PF.
Zhang et al. [39] proposed a multi-objective discrete artificial bee colony algorithm based on decomposition called MDABC for EHFSP. It can be directly applied to solve our EHFSP through adding a speed selection string and neighborhood structure speed. Thus, we selected it as the first comparative algorithm. Another algorithm used was the energyaware multi-objective optimization algorithm (EAMOA) [67], which is highly effective in solving EHFSP with makespan and energy consumption.

Parameter Settings
MPABC contains four main parameters: N, s, R and Limit. Taguchi method is conducted on instance 80 × 6. Four levels for each parameter are shown in Table 3. The orthogonal array L 16 4 4 is given in Table 4.  MPABC with each combination ran 20 times independently. The metric DI R was calculated and shown in Table 3. The trend of factor levels is demonstrated in Figure 4, and the rank of each parameter is displayed in Table 5. As shown in Figure 4, the best settings are N =100, S = 10, R = 12, Limit = 10.

Results and Discussion
MPABC was compared with MDABC, ABC and EAMOA. For each instance, al algorithms independently ran 20 times and took an average of 20 times for the final result Each algorithm ran under the same stopping criterion: the maximum CPU elapsed tim n m t   , where t was set to 100 for all instances.
The computational results about three metrics are, respectively reported in Tables 6-8. The reference set *  consists of the non-dominated solutions of 1  are non-dominated solutions gained by MPABC, ABC MDABC and EAMOA, respectively. In Table 7, "P" denotes MPABC, "D" indicate MDABC, "A" means ABC and "E" represents EAMOA. Figure 5 provides the box plot o all algorithms on metrics R DI , c and nd . Figure 6 exhibits the distribution curves o non-dominated solutions produced by four algorithms in instances 60 8  and 90 6  .
As shown in Tables 6-8, MPABC generated smaller R DI than ABC in 44 instances It also had higher   c P,A than   c A,P in 44 instances; moreover,   c P,A was equal to 1 in 43 instances, that is, the total solutions of ABC were dominated by the obtained solutions of MPABC. Similarly, the metric nd of MPABC was far greater than that o ABC in all instances, that is, ABC could not make any contribution to the reference set *  MPABC displayed a better performance than ABC. This conclusion can also be reached  All parameter settings of MDABC, ABC and EAMOA were acquired from the data of the reference.

Results and Discussion
MPABC was compared with MDABC, ABC and EAMOA. For each instance, all algorithms independently ran 20 times and took an average of 20 times for the final result. Each algorithm ran under the same stopping criterion: the maximum CPU elapsed time n × m × t, where t was set to 100 for all instances.
The computational results about three metrics are, respectively reported in Tables 6-8. The reference set Ω * consists of the non-dominated solutions of Ω 1 ∪ Ω 2 ∪ Ω 3 ∪ Ω 4 , where Ω 1 ,Ω 2 , Ω 3 and Ω 4 are non-dominated solutions gained by MPABC, ABC, MDABC and EAMOA, respectively. In Table 7, "P" denotes MPABC, "D" indicates MDABC, "A" means ABC and "E" represents EAMOA. Figure 5 provides the box plot of all algorithms on metrics DI R , c and nd. Figure 6 exhibits the distribution curves of non-dominated solutions produced by four algorithms in instances 60 × 8 and 90 × 6.   As shown in Tables 6-8, MPABC generated smaller DI R than ABC in 44 instances. It also had higher c(P, A) than c(A, P) in 44 instances; moreover, c(P, A) was equal to 1 in 43 instances, that is, the total solutions of ABC were dominated by the obtained solutions of MPABC. Similarly, the metric nd of MPABC was far greater than that of ABC in all instances, that is, ABC could not make any contribution to the reference set Ω * . MPABC displayed a better performance than ABC. This conclusion can also be reached by analyzing Figures 5 and 6. The significant difference between MPABC and ABC discloses that the multi-population scheme, different search strategies and the utilization of historical optimization data are extremely beneficial to the performance of MPABC.

Instance c(P,A) c(A,P) c(P,D) c(D,P) c(P,E) c(E,P)
In addition, results of MPABC are distinctly different from MDABC in most instances within the specified time. DI R of MPABC is smaller than that of MDABC in 38 instances and larger than that of MDABC in only 6 instances. Moreover, c(P, D) was superior to c(D, P) in 34 instances, and c(P, D) was equal to 1 in 12 instances. That is, all solutions of MDABC were dominated by those of MPABC in 12 instances. On the other hand, the nd metric produced by MDABC was much smaller than that of MPABC in 37 instances.
Compared with EAMOA, MPABC generated lower DI R in 39 instances and obtained a larger c(P, E) than c(E, P) in 31 instances. Moreover, the non-dominated solutions of MPABC dominated those of EAMOA in 12 instances. On the other hand, the nd metric generated by EAMOA was less than that of MPABC in 37 instances. Obviously, MPABC has more significant advantages over three metrics than ABC, MPABC and EAMOA.

A Real-Life Case Study
A real-life case study in a roll enterprise machining workshop was considered. There are three main processes in the crafting process: turning, milling and grinding, and three types of machines available for the processing. The first type is composed of three machines for turning (  To further verify the statistical difference in the algorithm performance, a Wilcoxon signed-rank test was conducted to compare MPABC with other algorithms. The test results are shown in Table 9. A significance level of 0.05 is assumed, which means that A is significantly different from B statistically if the asymptotically significance value is less than 0.05. This conclusion can also be drawn from Figures 5 and 6 and Table 9. In Figure 5, the median, upper and lower edges of MPABC on DI R , c and nd are superior to those of ABC, MPABC and EAMOA. Figures 5 and 6 demonstrate that solutions obtained by MPABC had better convergence and distribution. On the other hand, the asymptotically significance value was smaller than 0.05. That is to say, the statistical results of MPABC are remarkably better than other three algorithms.
MPABC performs better than comparative algorithms in solving EHFSP, which results from the novel multi-population division method, different search strategies in each group and the utilization of historical optimization data. The multi-population method is effective to maintain population diversity, different search strategies in groups are important to balance exploration and exploitation, and the utilization of historical data can rapidly and easily improve search efficiency. These features are useful to maintain the population diversity and avoid falling premature. It is concluded from above results and analysis that MPABC is more competitive in EHFSP with a variable speed constraint.

A Real-Life Case Study
A real-life case study in a roll enterprise machining workshop was considered. There are three main processes in the crafting process: turning, milling and grinding, and three types of machines available for the processing. The first type is composed of three machines for turning (M 1 ∼ M 3 ), the second type consists of three machines for milling (M 4 ∼ M 6 ), and the third class includes four machines for grinding (M 7 ∼ M 10 ). Its production process is a typical hybrid flow shop problem. The calculation method is shown in Section 4.
There are 30 jobs need to be processed. A speed set V = {1.0, 1.3, 1.55, 1.8, 2.0} is available for machines. The processing data is shown in Table 10. The non-dominated solutions obtained by MPABC are shown in Table 11. An optimal schedule is presented in Figure 7, in which the makespan is 16.35, the total tardiness is 56.92 and the total energy consumption is 780.48. From above results, the manufacturers can obtain a suitable schedule scheme to make a trade-off between makespan, TT and TEC. Therefore, optimizing the speed selection sub-problem with job scheduling and machine assignment sub-problems simultaneously is meaningful; meanwhile, MPABC is of great significance to the research of EHFSP.  7. An optimal schedule of the real-life case.

Conclusions
In this paper, we proposed a novel MPABC for EHFSP with a variable speed constraint. The proposed MPABC concerns three objectives, namely makespan, total tardiness and TEC. A new multi-population division method was embedded to achieve the group division. The division of three groups made it possible for each group to own exclusive food sources and perform different search strategies, which helps to increase the population diversity and improves the search efficiency. The weakest group is excluded from the evolution and renewed through historical optimization data, which can enhance the quality of solutions. Subpopulations are measured by the newly designed evaluation index. Finally, extensive comparative experiments with ABC, MDABC and EAMOA were conducted to verify the algorithm performance. The results demonstrate that the proposed MPABC can generate more competitive solutions in solving EHFSP. This work has not yet considered the feedback in evolution or the cooperation and competition relationship between groups. In the future research, these factors will be fully discussed, and more effective search strategies such as Q-learning can be explored to enhance the search performance. Accordingly, more intelligent optimization algorithms can be introduced to solve EHFSP. In addition, we are interested in investigating the capability of MPABC in the energy-efficient distributed scheduling problem, fuzzy flow shop scheduling problem and other engineering optimization problems.
Author Contributions: Y.Z. contributed to the conception of the study, the background research, method design and analysis of the experimental results, and wrote the manuscript; Z.F. provided the methodological guidance, writing review and funding support; T.Z. reviewed and revised the  Figure 7. An optimal schedule of the real-life case.

Conclusions
In this paper, we proposed a novel MPABC for EHFSP with a variable speed constraint. The proposed MPABC concerns three objectives, namely makespan, total tardiness and TEC. A new multi-population division method was embedded to achieve the group division. The division of three groups made it possible for each group to own exclusive food sources and perform different search strategies, which helps to increase the population diversity and improves the search efficiency. The weakest group is excluded from the evolution and renewed through historical optimization data, which can enhance the quality of solutions. Subpopulations are measured by the newly designed evaluation index. Finally, extensive comparative experiments with ABC, MDABC and EAMOA were conducted to verify the algorithm performance. The results demonstrate that the proposed MPABC can generate more competitive solutions in solving EHFSP. This work has not yet considered the feedback in evolution or the cooperation and competition relationship between groups. In the future research, these factors will be fully discussed, and more effective search strategies such as Q-learning can be explored to enhance the search performance. Accordingly, more intelligent optimization algorithms can be introduced to solve EHFSP. In addition, we are interested in investigating the capability of MPABC in the energy-efficient distributed scheduling problem, fuzzy flow shop scheduling problem and other engineering optimization problems.