# Heuristics for the Hybrid Flow Shop Scheduling Problem with Sequence-Dependent Setup times

Authors: Liao Yong; Li Zhantao; Li Xiang; Peng Chenfeng
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

This paper considers the hybrid flow shop scheduling problem with two types of sequence-dependent setup times, one depending on the order of jobs and the other depending on both the order of jobs and the machines where the jobs are allocated. Three constructive heuristic algorithms based on the NEH algorithm are developed with the object of minimizing the total energy consumption cost. The first algorithm, called EPRA, obtains the order of jobs using the energy consumption cost based on the processing time. The second algorithm, called ESRA, obtains the order of jobs using the energy consumption cost based on the setup time. The third algorithm, called ESPRA, takes advantage of both EPRA and ESRA by constructing fictitious jobs with ESRA and acquiring the sequence of fictitious jobs with EPRA. Due to the drawback of NEH algorithm for the FFS problem proposed in this paper, all three heuristic algorithms use EST rules to allocate the jobs to the machines. Two lower bounds are developed to evaluate the heuristic algorithms proposed in this paper. A computational experiment based on the split-plot design is developed to carry out the factor analysis, algorithm analysis, lower bounds analysis, and CPU running time analysis. The analysis reveals that all main factor interactions except the algorithm are significant; the ESPRA algorithm is generally outstanding among the three algorithms; LB1 is better than LB2; the CPU running time of ESPRA algorithm is significantly reduced compared to those of EPRA and ESST algorithms.


---

## Body

## 1. Introduction

The general flexible flow shop scheduling (FFS) (or hybrid flow shop scheduling (HFS)) problems have been found frequently in the literature. This problem is defined that a set of jobs are sequentially processed through some stages, in which each stage has some parallel machines, and each job is processed by one machine in each stage. Some assumptions are proposed: all jobs are released at time zero; the setup times of jobs are sequence-independent or included in the processing times.There are quite a number of potential variants of the HFS problem. In the beginning, most previous research studies focused on optimizing the production efficiency, e.g., minimizing makespan, the total completion time of jobs, and total tardiness. In recent years, there has been a growing interest in the development of energy saving due to the serious environmental impacts and the increment of energy costs. Thus, this paper considers minimizing the objective of energy consumption cost.Another key characteristic of the HFS problem addressed in this paper is that the setup time depends on the sequence of jobs and separated from the processing time. In recent years, a number of scheduling problems considering setup time have been found, and some reviewing papers are published until 2021 [1–3]. In the light of these reviews, there are two different variations of setup time: setup time depending on the job to be processed (sequence-independent setup time) and setup time depending on the sequence of jobs (sequence-dependent setup time). In this paper, we consider an HFS problem with two types of sequence-dependent setup time. In this research, we define two types of sequence-dependent setup time: one depending on the sequence of jobs and the other depending on both the sequence of jobs and the machine where the jobs are allocated. The setup time at each stage belongs to one of the two types. This assumption is derived from a real-life manufacturing environment faced by our cooperative company, which is one of the biggest architectural ceramic manufacturers in China. The ceramic process of this company mainly consists of five steps, namely, ball milling of original material, spraying, forming, drying, and polishing. Production switching costs plenty of time in these steps; for example, for the forming process, the setup time derived from the transformation of different colors of products depends on the sequence of the products; for the drying process, the workers need to clean the machines when the product is transformed to another. The cleaning time of each machine varies according to the residues from the processing of the last product. Thus, the setup time in the drying process depends on both the sequence of the product and the machines to which the products are allocated.


## 2. Literature Review

In this section, we review the literature from two aspects: HFS with the objective energy consumption and HFS with sequence-dependent setup time.


### 2.1. HFS with Energy Consumption

There are few of papers related to this topic which [4] address a multiobjective flexible flow shop batch scheduling problem with renewable energy. They propose a genetic algorithm to solve it and [5] investigates an energy-conscious hybrid flow shop scheduling problem with unrelated parallel machines with the energy-saving strategy of turning off and on.


### 2.2. HFS with Sequence-Dependent Setup Time (HFS-SDST)

Table1 is the literature review of the HFS-SDST problem from 2005 to 2021. For the type of setup time, Refs. [6, 14, 18, 22, 24, 30, 34] consider the family job sequence-dependent setup time; Refs. [19, 28] consider two types of setup time, but machines do not have multiple patterns; others consider the job sequence-dependent setup time. For the algorithm, most of them propose intelligent algorithms to solve this problem [6–25, 27–29, 31, 33, 35]. Fewer focus on heuristics [26, 28, 30–32, 34]. There is only one paper that focused on exact algorithm [36].Table 1
Literature review of the HFS-SDST problem.
Ref.ProblemSetup typeCriterionApproach[6]FFSMj,SijC∑LimaxSDFMakespan; total tardinessGA[7]FFSMj,SijkC∑CimaxSDMakespan; total completion timeImproved MH[8]FFSMj,SijC∑LimaxSDMakespan; total tardinessImproved SA[9]HFSMj,SijCmaxSDMakespanLeaning IG[10]HFSMj,Sij∑Lj,TSTSDTotal tardiness; total setup timesEA[11]FFSMj,Sij∑ωiCiSDThe total weighted completion timeABC[12]FFSMj,Sij∑LiSDTotal tardinessAn integrated approach based on ANN, GA, and CS[13]HFSMj,SijCmaxInner SDMakespanRegeneration GA[14]FFSMj,SijPMC,PRC,LiSDFPM cost; repair cost; tardinessSA with GA[15]HFSMj,SijC∑LimaxSDMakespan; total tardinessImproved GA[16]FFSMj,Sij∑LiSDTotal tardinessACS[17]HFSMj,SijCmaxmaxSDMakespan; tardinessEA[18]FFSMj,Sij∑CiSDFTotal completion timeMA[19]FFSMj,Sij∑ωiLi,CmaxSD/MDMakespan; the total weighted tardinessPareto approach[20]FFSMj,SijCmaxSDFMakespanGA[21]HFSMj,SijCmaxSDMakespanImproved GA[22]FFSMj,SijC∑ωiLimaxSDFMakespan; total weighted tardinessMultipopulation GA[23]HFSMj,Sij∑di−CiSDThe sum of earliness and tardinessHybrid method based on AC and SA[24]FFSMj,SijC∑ωiLimaxSDFMakespan; total weighted tardinessA multiphase approach[25]FFSMj,Sij∑ωiLiSDTotal weighted tardinessElectromagnetism algorithm[26]HFSMj,SijCmaxSDMakespanHeuristics based on SPT, LPT, and Johnson[27]HFSMj,SijC∑LimaxSDMakespan; total tardinessSA[28]FFSMj,Sij∑Ci,∑UiSD/MDTotal completion time; the number of tardy jobsHeuristics, SA, GA, TS[29]HFSMj,SijCmaxSDMakespanIA[30]FFSMj,SijCmaxSDFMakespanHeuristics[31]HFSMj,SijCmaxSDMakespanGreed method; flow line method; NEH; GA[32]FFSMj,SijCmaxSDMakespanHeuristics[33]FFSMj,SijCmaxSDMakespanTS[34]FFSMj,SijC∑LimaxSDFMakespan; total tardinessHeuristics[35]FFSMj,SijCmaxSDMakespanHybrid algorithm on GA and TS[36]FF2Mj,SijCmaxSDMakespanB&B; TS; heuristics


### 2.3. Summary of the Literature Review

Three observations are obtained according to the brief literature review. First, for the shop environments, few references focus on the HFS-SDST problem. Second, this problem addressed in this paper is one variation of the HFS-SDST problem, where the setup times depend on the sequence of jobs and the machine selected. Finally, most of the scheduling problems with setup times focus on the objective of time-related criterion and less of them have been reported on optimizing the energy efficiency.


## 2.1. HFS with Energy Consumption

There are few of papers related to this topic which [4] address a multiobjective flexible flow shop batch scheduling problem with renewable energy. They propose a genetic algorithm to solve it and [5] investigates an energy-conscious hybrid flow shop scheduling problem with unrelated parallel machines with the energy-saving strategy of turning off and on.


## 2.2. HFS with Sequence-Dependent Setup Time (HFS-SDST)

Table1 is the literature review of the HFS-SDST problem from 2005 to 2021. For the type of setup time, Refs. [6, 14, 18, 22, 24, 30, 34] consider the family job sequence-dependent setup time; Refs. [19, 28] consider two types of setup time, but machines do not have multiple patterns; others consider the job sequence-dependent setup time. For the algorithm, most of them propose intelligent algorithms to solve this problem [6–25, 27–29, 31, 33, 35]. Fewer focus on heuristics [26, 28, 30–32, 34]. There is only one paper that focused on exact algorithm [36].Table 1
Literature review of the HFS-SDST problem.
Ref.ProblemSetup typeCriterionApproach[6]FFSMj,SijC∑LimaxSDFMakespan; total tardinessGA[7]FFSMj,SijkC∑CimaxSDMakespan; total completion timeImproved MH[8]FFSMj,SijC∑LimaxSDMakespan; total tardinessImproved SA[9]HFSMj,SijCmaxSDMakespanLeaning IG[10]HFSMj,Sij∑Lj,TSTSDTotal tardiness; total setup timesEA[11]FFSMj,Sij∑ωiCiSDThe total weighted completion timeABC[12]FFSMj,Sij∑LiSDTotal tardinessAn integrated approach based on ANN, GA, and CS[13]HFSMj,SijCmaxInner SDMakespanRegeneration GA[14]FFSMj,SijPMC,PRC,LiSDFPM cost; repair cost; tardinessSA with GA[15]HFSMj,SijC∑LimaxSDMakespan; total tardinessImproved GA[16]FFSMj,Sij∑LiSDTotal tardinessACS[17]HFSMj,SijCmaxmaxSDMakespan; tardinessEA[18]FFSMj,Sij∑CiSDFTotal completion timeMA[19]FFSMj,Sij∑ωiLi,CmaxSD/MDMakespan; the total weighted tardinessPareto approach[20]FFSMj,SijCmaxSDFMakespanGA[21]HFSMj,SijCmaxSDMakespanImproved GA[22]FFSMj,SijC∑ωiLimaxSDFMakespan; total weighted tardinessMultipopulation GA[23]HFSMj,Sij∑di−CiSDThe sum of earliness and tardinessHybrid method based on AC and SA[24]FFSMj,SijC∑ωiLimaxSDFMakespan; total weighted tardinessA multiphase approach[25]FFSMj,Sij∑ωiLiSDTotal weighted tardinessElectromagnetism algorithm[26]HFSMj,SijCmaxSDMakespanHeuristics based on SPT, LPT, and Johnson[27]HFSMj,SijC∑LimaxSDMakespan; total tardinessSA[28]FFSMj,Sij∑Ci,∑UiSD/MDTotal completion time; the number of tardy jobsHeuristics, SA, GA, TS[29]HFSMj,SijCmaxSDMakespanIA[30]FFSMj,SijCmaxSDFMakespanHeuristics[31]HFSMj,SijCmaxSDMakespanGreed method; flow line method; NEH; GA[32]FFSMj,SijCmaxSDMakespanHeuristics[33]FFSMj,SijCmaxSDMakespanTS[34]FFSMj,SijC∑LimaxSDFMakespan; total tardinessHeuristics[35]FFSMj,SijCmaxSDMakespanHybrid algorithm on GA and TS[36]FF2Mj,SijCmaxSDMakespanB&B; TS; heuristics


## 2.3. Summary of the Literature Review

Three observations are obtained according to the brief literature review. First, for the shop environments, few references focus on the HFS-SDST problem. Second, this problem addressed in this paper is one variation of the HFS-SDST problem, where the setup times depend on the sequence of jobs and the machine selected. Finally, most of the scheduling problems with setup times focus on the objective of time-related criterion and less of them have been reported on optimizing the energy efficiency.


## 3. Problem Definition

### 3.1. Problem Statement

In this paper, we deal with a HFS-SDST with two types of sequence-dependent setup times. A schedule of this type is defined as follows. There aren jobs or tasks, each of which has to be processed at k stages in sequence. Each stage consists of mss=1,2,…,k unrelated parallel machines. Every machine has three modes: operational, standby, and stopping. The energy consumption of each machine consists of three parts responding to its three modes: running energy consumption with the operational mode, idle energy consumption with the standby mode, and warm-up energy consumption with mode shifting stopping to operation. An operation of job jj=1,2,3,…n at stage s should be processed on one of the ms unrelated parallel machines. Jobs are not available at zero time and each job has a release time rj. Setup time exists between two consecutive jobs allocated on the same machine. A key characteristic of this paper is that the setup times at some stages are only dependent on the sequence of jobs, while the setup times at other some stages are not only dependent on the sequence of jobs but also the machine which the jobs are allocated on. There are a few assumptions which are as follows:There is unlimited buffer or storage capacity between any two stages.(i)
All machines are simultaneously available, and there is no machine breakdown(ii)
All machines are initialized to be in the stop mode, with no energy consumption until the first job allocated to them is to be processed(iii)
Jobs have to be processed sequentially without overlapping between two successive stages(iv)
Jobs are processed without pre-emption on any machine(v)
Each machine can process at the most only one job at a time(vi)
The setup time is assumed to be asymmetric; that is, the setup time of pair of jobsi,j is different with that of pair of jobs j,i.The objective is to find a schedule of jobs to minimize the cost of energy consumptionE_min. The above problem is known as the SDST flexible flow shop scheduling problem and can be denoted by FFSrj,m1,…,mkEmin, with the three-field notation for the scheduling problem proposed by Graham et al (Appendix is an instance of FFSrj,m1,…,mkEmin.


### 3.2. Mathematical Model

The problem addressed in this paper can be expressed formally as an integer program. Before presenting this mathematical model, we introduce the notations used in the integer program.Basic notationsn: number of jobs to be scheduledk: number of stagesrj: the release time of job j, j=1,2,3,…,nms: number of unrelated parallel machines at stage s, s=1,2,3,…,kPjhs: processing time of job j on machine h at stage s, h=1,2,3,…,kBjs: starting time of job j at stage sCjs: completion time of job j at stage sQhs: set of jobs allocated on machine h at stage sths: warm-up time of machine h at stage s with the mode from stopping to operationΔths: set of idle times of machine h at stage sΔtwhs: w - th idle time of machine h at stage s, w∈Δthstijg1: setup time form job i to job j at stage g where g∈ST1 and tijg1≠tjig1tijhe2: setup time form job i to job j on machine h at stage ewheree∈ST2andtijhe2≠tjihe2Ehs1: per unite time of running energy consumption cost of machine h at stage sEhs2: per unite time of idle energy consumption cost of machine h at stage sEhs3: per unite time of warm-up energy consumption cost of machine h at stage sM: a large numberDecision variableXjhs=1 if job j is processed on machine h at stage s and 0 otherwiseYijhs=1 if job j is processed immediately behind job i on machine h at stage s and 0 otherwiseUwhs=1 if machine h at stage s with mode shifting stopping to operation in Δtwhs and 0 otherwise, w∈ΔthsThis problem can be now formulated as follows:(1)minimize Emin=∑s=1k∑h=1msEhs1pjhsXjhs+∑s=1k∑h=1msEhs3thsUwhs+∑s=1k∑h=1msEhs2Δtwhs1−Uwhsthat is subject toSubject to(2)Cjs−1≤Bjs,j=1,2,⋯,n,s=1,2,⋯k,(3)Bj1≤rj,j=1,2,⋯,n(4)∑h=1msXjhs=1,j=1,2,⋯,ns=1,2,⋯k(5)∑j,i≠jYijhs=1,Bjhs+M1-Yijhs≥Cihs.(6)i=1,…,n,j=1,…,n,s=1,2,…k,h=1,2,…,ms,Cjs=Bjs+∑h=1msXjhspjhs+tijhs,i,j∈Qhs,s=1,2,…k.(7)Bjs=maxCjhs−1,Cj−1hs;j=1,2,…,n,s=1,…k,h=1,2,…,ms(8)Δtwhs=Bjhs−Cj−1hs,Δtwhs≥ths,Ehs2⋅Δtwhs−Ehs3⋅ths≥0,Uwhs=1.(9)h=1,2,…,ms,s=1,2,…,k,w∈ΔThs,j∈Qhs,Yijhs=0Yijhs=0,i=1,…,n,j=1,…,n,i=j,s=1,2,…k,h=1,2,…,ms.(10)Xjhs,YijhsandUwhs∈1,0Xjhs,i=1,⋯,n,j=1,⋯,n,s=1,2,⋯k,h=1,2,⋯,ms.In the above formulations, the objective is to determine a schedule that minimizes the cost of energy consumption seen in equation (1). Emin consists of three parts: the cost of running energy consumption (the first term of equation (1)), standby energy consumption (the second term of equation (1)), and warm-up energy consumption (the third term of equation (1)). Constraint (2) forces that every job is processed sequentially at k stages; that is, a job cannot be started to be processed at stage j before it has not been finished at the previous one. Constraint (3) forces that the starting time of each job at the first stage is not allowed being earlier than its release time. Constraint (4) ensures that each job is processed simultaneously on one and only one machine. Constraint (5) ensures that a machine cannot process two different jobs simultaneously. If Yijhs=1, then the starting time of job j is larger than the completion time of job i where jobs i, j are allocated on the same machine. It means that for the two continuous jobs i, j allocated on the same machine h, machine h cannot start process job j until it finishes to process job i. If Yijhs=0, then constraint (5) holds where M is a large number. Constraint (6) calculates the completion time of every job at each stage, which is equal to the starting time of job j at stage s plus the processing time of job j allocated on machine h at stage s plus the setup time from its predecessor to j. Because there are two types of setup time, then the completion time of job j at stage s is calculated by the upper formulation of Constraint (6) when the setup time depends on the sequence of jobs, while it is calculated by the formulation of Constraint (6) when the setup time depends on both the sequence of jobs and the machine in which the jobs are allocated on. Constraint (7) defines the starting time of job j at stage s, which is set to the maximum between the completion time of job j at the previous stage and the available time of the machine processing job j. If s = 1, Cjhs−1=rj, and j=1, then C0hs=0. Constraint (8) determines the mode of a machine (standby mode or mode shifting, stopping on operation) in the idle time. The first formulation of Constraint (8) defines the idle time. If it satisfies the second and third formulations of Constraint (8), i.e., if the idle time of the machine is larger than its warm-up time, and the cost of warm-up energy consumption is larger than the cost of idle energy consumption in the idle time, then Uwhs=1, namely, the mode of the machine is from stopping to operation, or the mode of machine is standby Uwhs=0. Constraints (9) and (10) provide limits on the decision variables.


### 3.3. Lower Bounds

This section contains a theorem with two lower bounds for theHFS|rj,m1,…,mk|Emin problem. The notations Lg and Le are binary variables. Similar to Uwhs, Lg=1 denotes that the jobs consume the warm-up energy, while Lg=0 denotes that the jobs consume the standby energy in the setup time at stage g∈ST1 when it satisfies Constraint (8) where minj≠itjig is used instead of Δtwhs. Le has the similar indication to Lg when e∈ST2.Theorem 1.
The following are lower bounds on any feasible solution to the problemHFS|rj,m1,…,mk|Emin:(11)LB1=∑i=1n∑s=1kminh=1,⋯,mspihsEhs1+1−Lsminj≠itjihsminh=1,⋯,msEhs2+Lsminh=1,⋯,msEhs3thsLB2=∑s=1k∑h=1ms∑i=1QhspihsEhs1+∑s=1k∑h=1ms1−LsQhs−1minj≠i;i,j∈QhstjihsEhs2+LsQhs−1Ehs3ths.Proof.
LB1: This is a job-based bound. Every job i must be processed at each stage and must have setup time, both of which require at least the minimal amount of energy consumption cost. During the processing procedure of job i, it consists of three parts corresponding to the machine with three modes. The three terms of LB1, respectively, indicate the minimal of running energy consumption cost, standby energy consumption cost, and warm-up energy consumption cost.LB2: This is a machine-based bound. Each machine at every stage has three parts of energy consumption. The first term of LB2 represents the amount of running energy consumption cost, which cannot be improved. The second term of LB2 represents the minimal amount of idle energy consumption with gϵTS1. For each machine, the number of setup time is at least (Qhg−1). Therefore, the minimum of idle energy consumption cost for each machine is (Qhg−1) times of idle energy consumption cost based on the minimal setup time of jobs allocated to that machine. The third term of LB2 represents the minimal amount of idle energy consumption with eϵTSe, which is similar to the second term of LB2.


## 3.1. Problem Statement

In this paper, we deal with a HFS-SDST with two types of sequence-dependent setup times. A schedule of this type is defined as follows. There aren jobs or tasks, each of which has to be processed at k stages in sequence. Each stage consists of mss=1,2,…,k unrelated parallel machines. Every machine has three modes: operational, standby, and stopping. The energy consumption of each machine consists of three parts responding to its three modes: running energy consumption with the operational mode, idle energy consumption with the standby mode, and warm-up energy consumption with mode shifting stopping to operation. An operation of job jj=1,2,3,…n at stage s should be processed on one of the ms unrelated parallel machines. Jobs are not available at zero time and each job has a release time rj. Setup time exists between two consecutive jobs allocated on the same machine. A key characteristic of this paper is that the setup times at some stages are only dependent on the sequence of jobs, while the setup times at other some stages are not only dependent on the sequence of jobs but also the machine which the jobs are allocated on. There are a few assumptions which are as follows:There is unlimited buffer or storage capacity between any two stages.(i)
All machines are simultaneously available, and there is no machine breakdown(ii)
All machines are initialized to be in the stop mode, with no energy consumption until the first job allocated to them is to be processed(iii)
Jobs have to be processed sequentially without overlapping between two successive stages(iv)
Jobs are processed without pre-emption on any machine(v)
Each machine can process at the most only one job at a time(vi)
The setup time is assumed to be asymmetric; that is, the setup time of pair of jobsi,j is different with that of pair of jobs j,i.The objective is to find a schedule of jobs to minimize the cost of energy consumptionE_min. The above problem is known as the SDST flexible flow shop scheduling problem and can be denoted by FFSrj,m1,…,mkEmin, with the three-field notation for the scheduling problem proposed by Graham et al (Appendix is an instance of FFSrj,m1,…,mkEmin.


## 3.2. Mathematical Model

The problem addressed in this paper can be expressed formally as an integer program. Before presenting this mathematical model, we introduce the notations used in the integer program.Basic notationsn: number of jobs to be scheduledk: number of stagesrj: the release time of job j, j=1,2,3,…,nms: number of unrelated parallel machines at stage s, s=1,2,3,…,kPjhs: processing time of job j on machine h at stage s, h=1,2,3,…,kBjs: starting time of job j at stage sCjs: completion time of job j at stage sQhs: set of jobs allocated on machine h at stage sths: warm-up time of machine h at stage s with the mode from stopping to operationΔths: set of idle times of machine h at stage sΔtwhs: w - th idle time of machine h at stage s, w∈Δthstijg1: setup time form job i to job j at stage g where g∈ST1 and tijg1≠tjig1tijhe2: setup time form job i to job j on machine h at stage ewheree∈ST2andtijhe2≠tjihe2Ehs1: per unite time of running energy consumption cost of machine h at stage sEhs2: per unite time of idle energy consumption cost of machine h at stage sEhs3: per unite time of warm-up energy consumption cost of machine h at stage sM: a large numberDecision variableXjhs=1 if job j is processed on machine h at stage s and 0 otherwiseYijhs=1 if job j is processed immediately behind job i on machine h at stage s and 0 otherwiseUwhs=1 if machine h at stage s with mode shifting stopping to operation in Δtwhs and 0 otherwise, w∈ΔthsThis problem can be now formulated as follows:(1)minimize Emin=∑s=1k∑h=1msEhs1pjhsXjhs+∑s=1k∑h=1msEhs3thsUwhs+∑s=1k∑h=1msEhs2Δtwhs1−Uwhsthat is subject toSubject to(2)Cjs−1≤Bjs,j=1,2,⋯,n,s=1,2,⋯k,(3)Bj1≤rj,j=1,2,⋯,n(4)∑h=1msXjhs=1,j=1,2,⋯,ns=1,2,⋯k(5)∑j,i≠jYijhs=1,Bjhs+M1-Yijhs≥Cihs.(6)i=1,…,n,j=1,…,n,s=1,2,…k,h=1,2,…,ms,Cjs=Bjs+∑h=1msXjhspjhs+tijhs,i,j∈Qhs,s=1,2,…k.(7)Bjs=maxCjhs−1,Cj−1hs;j=1,2,…,n,s=1,…k,h=1,2,…,ms(8)Δtwhs=Bjhs−Cj−1hs,Δtwhs≥ths,Ehs2⋅Δtwhs−Ehs3⋅ths≥0,Uwhs=1.(9)h=1,2,…,ms,s=1,2,…,k,w∈ΔThs,j∈Qhs,Yijhs=0Yijhs=0,i=1,…,n,j=1,…,n,i=j,s=1,2,…k,h=1,2,…,ms.(10)Xjhs,YijhsandUwhs∈1,0Xjhs,i=1,⋯,n,j=1,⋯,n,s=1,2,⋯k,h=1,2,⋯,ms.In the above formulations, the objective is to determine a schedule that minimizes the cost of energy consumption seen in equation (1). Emin consists of three parts: the cost of running energy consumption (the first term of equation (1)), standby energy consumption (the second term of equation (1)), and warm-up energy consumption (the third term of equation (1)). Constraint (2) forces that every job is processed sequentially at k stages; that is, a job cannot be started to be processed at stage j before it has not been finished at the previous one. Constraint (3) forces that the starting time of each job at the first stage is not allowed being earlier than its release time. Constraint (4) ensures that each job is processed simultaneously on one and only one machine. Constraint (5) ensures that a machine cannot process two different jobs simultaneously. If Yijhs=1, then the starting time of job j is larger than the completion time of job i where jobs i, j are allocated on the same machine. It means that for the two continuous jobs i, j allocated on the same machine h, machine h cannot start process job j until it finishes to process job i. If Yijhs=0, then constraint (5) holds where M is a large number. Constraint (6) calculates the completion time of every job at each stage, which is equal to the starting time of job j at stage s plus the processing time of job j allocated on machine h at stage s plus the setup time from its predecessor to j. Because there are two types of setup time, then the completion time of job j at stage s is calculated by the upper formulation of Constraint (6) when the setup time depends on the sequence of jobs, while it is calculated by the formulation of Constraint (6) when the setup time depends on both the sequence of jobs and the machine in which the jobs are allocated on. Constraint (7) defines the starting time of job j at stage s, which is set to the maximum between the completion time of job j at the previous stage and the available time of the machine processing job j. If s = 1, Cjhs−1=rj, and j=1, then C0hs=0. Constraint (8) determines the mode of a machine (standby mode or mode shifting, stopping on operation) in the idle time. The first formulation of Constraint (8) defines the idle time. If it satisfies the second and third formulations of Constraint (8), i.e., if the idle time of the machine is larger than its warm-up time, and the cost of warm-up energy consumption is larger than the cost of idle energy consumption in the idle time, then Uwhs=1, namely, the mode of the machine is from stopping to operation, or the mode of machine is standby Uwhs=0. Constraints (9) and (10) provide limits on the decision variables.


## 3.3. Lower Bounds

This section contains a theorem with two lower bounds for theHFS|rj,m1,…,mk|Emin problem. The notations Lg and Le are binary variables. Similar to Uwhs, Lg=1 denotes that the jobs consume the warm-up energy, while Lg=0 denotes that the jobs consume the standby energy in the setup time at stage g∈ST1 when it satisfies Constraint (8) where minj≠itjig is used instead of Δtwhs. Le has the similar indication to Lg when e∈ST2.Theorem 1.
The following are lower bounds on any feasible solution to the problemHFS|rj,m1,…,mk|Emin:(11)LB1=∑i=1n∑s=1kminh=1,⋯,mspihsEhs1+1−Lsminj≠itjihsminh=1,⋯,msEhs2+Lsminh=1,⋯,msEhs3thsLB2=∑s=1k∑h=1ms∑i=1QhspihsEhs1+∑s=1k∑h=1ms1−LsQhs−1minj≠i;i,j∈QhstjihsEhs2+LsQhs−1Ehs3ths.Proof.
LB1: This is a job-based bound. Every job i must be processed at each stage and must have setup time, both of which require at least the minimal amount of energy consumption cost. During the processing procedure of job i, it consists of three parts corresponding to the machine with three modes. The three terms of LB1, respectively, indicate the minimal of running energy consumption cost, standby energy consumption cost, and warm-up energy consumption cost.LB2: This is a machine-based bound. Each machine at every stage has three parts of energy consumption. The first term of LB2 represents the amount of running energy consumption cost, which cannot be improved. The second term of LB2 represents the minimal amount of idle energy consumption with gϵTS1. For each machine, the number of setup time is at least (Qhg−1). Therefore, the minimum of idle energy consumption cost for each machine is (Qhg−1) times of idle energy consumption cost based on the minimal setup time of jobs allocated to that machine. The third term of LB2 represents the minimal amount of idle energy consumption with eϵTSe, which is similar to the second term of LB2.


## 4. Construction Heuristic Algorithms

### 4.1. An Improved Algorithm Based on the NEH Algorithm (EPRA)

For the HFS scheduling problem, Ref. [37] proved that it is an NP problem even with two stages where one of the two stages has at least two machines. According to the complicated theorem, the problem HFS|rj,m1,…,mk|Emin is at least having that difficulty. There is no exact algorithm in polynomial time for even small size of NP problems. Therefore, heuristic algorithms are more effective for real manufacturing. Among these heuristic algorithms, constructive heuristic algorithms are designed for the problem and have their own merit (the advantage of computational time) compared with metaheuristic algorithms. However, the performances of constructive heuristic algorithms significantly depend on the different objectives and structures of the problems. Thus, for the new problems which even have a little difference with the existing problems, the existing constructive heuristic algorithms may be not effective. In the worst case, they cannot tackle those problems. Under such circumstances, existing heuristic algorithms should be developed or new algorithms should be proposed to solve the new problems. Ref. [38] proposed a heuristic algorithm called the NEH algorithm, which is probably considered as the best constructive heuristic method for the permutation flow shop makespan problem. The core ideal of the NEH algorithm is that a job with longer total processing time on the machine has the higher priority in the processing sequence. The steps of NEH algorithm are as follows (Algorithm 1):Algorithm 1: NEH algorithm.
Step 1: Jobs are sequenced by the decreasing order of their total processing time and obtain the sequence of jobs denoted bySStep 2: Take the first two jobs ofS as the partial schedule and find the optimal scheduleStep 3: Insert the next job ofS into the available positions of the previous partial schedule and find the partial schedule with the best objective function valueStep 4: Repeat the process for the remaining jobs according toSSince this paper considers the objective of minimization of energy consumption cost, we develop a new algorithm based on the NEH algorithm (called EPRA). The EPRA algorithm is similar to the NEH algorithm except for Step 1, i.e., constructing the sequence of jobs. In the EPRA algorithm, jobs are sequenced by the decreasing order of their summed running energy consumption cost based on the total processing times (denoted byepi) which is calculated by (12)epi=∑s=1k∑h=1mspihsEhs1ms,∀i,where ∑h=1mspihsEhs1/ms denotes the average running energy consumption cost based on the processing time at each stage. Since the NEH algorithm was originally proposed for the flow shop scheduling problem where each stage has only one machine, it does not have the problem that allocates the jobs to machines. For the problem FFS|rj,m1,…,mk|Emin addressed in this paper, the original NEH algorithm does not solve that problem. Thus, we propose a rule EST (the earliest starting time) to tackle that problem in Steps 2 and 3 of the EPRA algorithm. If there is one more machine with EST rule, then allocate the job to the machine that has the lowest running energy consumption cost.The major computation of EPRA algorithm has two operations: sorting the summed running energy consumption cost of jobs and inserting jobs to the optimal position of the partial sequence. The computational efficiency of the former operation isnlogn. The computational complexity of the second operation under the worst case is by (13). Since n2+n/2>nlogn, the computational complexity of EPRA algorithm is On2Δ∑s=1kms.(13)∑p=1np=∑s=1kms1+2+⋯+n=n2+n2∑s=1kms.In this paper, we develop two more new constructive heuristic algorithms based on the NEH algorithm with the objective of minimizing the cost of energy consumption. The remained algorithms are described in the following sections.


### 4.2. Energy Consumption Cost Based on Setup Time Ranking Algorithm (ESRA)

According to the above section, the NEH algorithm considers the processing time of jobs in constructing the sequence of jobs. However, one of the characteristics of the problem addressed in this paper is that there is setup time between two consecutive jobs and takes the minimization of total energy consumption as the objective. Therefore, we propose another heuristic algorithm called the ESRA algorithm which considers the energy consumption in the setup time of consecutive jobs and not the processing time of jobs. The fundamental principle behind this algorithm is that as the setup time between jobs is reduced, the cost of energy consumption also will get reduced. The ESRA algorithm begins with the summed energy consumption matrix (denoted by ESST) of the setup times of jobs at all stages. Thus, an elementesstij of ESST is obtained by adding the cost of energy consumption in the setup time of two consecutive jobs for all stages as follows:(14)esstij=∑g=1ST1∑h=1mgEhg2mg∑i=1ntijg1+∑e=1ST2∑h=1meEhe2me∑i=1n∑h=1metijhe2me,∀i,j,g∈ST1,e∈ST2Since the setup time has two types,esstij consists of two terms where the first term of equation (14) indicates the cost of idle energy consumption with type one setup time (dependent on the sequence of jobs) and the second term of equation (14) indicates the idle energy consumption cost with type two setup time (dependent on the sequence of jobs and the machine which the jobs are allocated to). In equation (14), we take the average unit per time of idle energy consumption cost of machines at each stage as the unit per time of idle energy consumption; that is, ∑h=1mgEhg2/mg and ∑h=1meEhe2/me respectively indicate the unit time of idle energy consumption cost with type one setup time and type two setup time.The ESST matrix gives a measure of idle energy consumption cost based on the setup time between jobs. The summed cost of each row in the ESST matrix can be used for ranking the combination of jobs from the lowest sum to the maximum sum to get the order of jobs S. The order S shows the best combination of jobs in terms of idle energy consumption cost based on setup time. The first job in S is the seed of the sequence. The complete sequence is built by adding the next lowest job to the best position in the partial sequence. Since the ESRA algorithm is similar to the EPRA algorithm, the computational complexity of ESRA algorithm isOn2Δ∑s=1kms (Algorithm 2).Algorithm 2: ESRA algorithm.
Step 1: According to equation (13), calculate the summed energy consumption cost based on the setup time of each job and order the jobs in an ascending manner . Thus, obtain the sequence of jobs S.Step 2: Take the first two jobs inS and schedule them in the best objective function. Delete the two jobs in S and save the best partial sequence OP.Step 3: FromO = 3 to n, take the O-th job in S and insert it into the O possible positions in OP. Then, select the best partial schedule. Finally, obtain the job processing order OP∗.Step 4: According toOP∗, schedule the jobs on the machines at each stage with EST rule.


### 4.3. Energy Consumption Cost Based on Setup Time and Processing Ranking Algorithm (ESPRA)

ESRA is an algorithm to predominantly consider the idle energy consumption cost based on setup time, and it does not concern the processing time for constructing the sequence. Conversely, the EPRA algorithm does not give importance to the setup time while constructing the sequence. In short, ESRA is a setup-time-based algorithm, and EPRA is a processing-time-based algorithm. In this paper, we propose an algorithm (denoted by ESPRA) by considering the setup time and processing time. The ESPRA algorithm considers the energy consumption based on setup time and processing time for constructing the sequence of jobs, giving importance to both of them. Similar to the ESRA algorithm, the ESPRA algorithm starts with the summed idle energy consumption matrix ESST based on setup time and then ranks the jobs by the ascending order of the summed row of ESST matrix to get an order of jobsS. The first job in S is the seed of the sequence. In the ESPRA algorithm, the sequence is constructed by using the concept of fictitious jobs instead of constructing the sequence by adding the next job in S on the possible positions of the seed or partial sequence. The concept of fictitious jobs is explained below.


#### 4.3.1. Fictitious Jobs

The key characteristic of ESPRA algorithm is the concept of fictitious jobs. Fictitious jobs are defined by the best pair of jobs from S. It means the fictitious jobs have two jobs in them. The advantage of this concept is that it can keep the pairs which have the smallest idle energy consumption cost based on the setup time between them. The ESPRA algorithm takes this advantage which gives importance to the setup time. After identifying the fictitious jobs, the processing time of fictitious jobs at each stageFJis is calculated by equation (14); that is, the sum of the average running energy consumption cost based on the processing time of each pair of jobs at each stage:(15)FJis=∑h=1mspihsΔEhs1+∑h=1mspi+1hsΔEhs1ms,i∈S,where pn+1hs = 0. Figure 1 shows an illustration of ESPRA algorithm at stage s. J1, J2, …, J8 are individual jobs to be sequenced in the ascending order of the summed idle energy consumption based on setup time. FJ1s, FJ2s, FJ3s, and FJ4s are the fictitious jobs. For example, the fictitious job FJ1s includes the individual jobs J1 and J2. Then, the processing time of a fictitious jobs is the sum of the average running energy consumption cost of each pair of jobs instead of their processing times. After identifying the processing times of the fictitious jobs, these jobs are ranked by decreasing the order of them. The rank one job will be first added to the sequence, which forms the partial sequence. The next eligible candidate will be the job with the next highest summed running energy consumption cost based on processing time. All possible insertion positions are considered for inserting the new job in the partial sequence so that the total energy consumption cost is minimum. For example, in Figure 1, the fictitious job FJ1s with the largest summed running energy consumption cost is the partial sequence. FJ2s with the next largest summed running energy consumption cost has two possible points in the partial sequence, then inserts FJ2s into these two points, and then selects the position with less total energy consumption cost (assume that FJ2s is the front of FJ1s). The partial sequence is constructed by repeating the above operations until all fictitious jobs are completed (Algorithm 3).Figure 1
An illustration of ESPRA algorithm at stages.Algorithm 3: ESPRA algorithm.
Step 1: Set the partial sequenceOP∗ = 0.Step 2: Calculate the matrix ESST, where each entity of ESSTesstij is calculated by equation (13)Step 3: Jobs are sorted by the ascending order of the corresponding summed idle energy consumption cost based on the setup time in ESST to obtain the sequence of jobss.Step 4: According tos, identify the pairs of jobs as fictitious jobs, and calculate their running energy consumption cost based on processing time by equation (14).Step 5: Sort the fictitious jobs by the decreasing order of their running energy consumption cost and obtain the sequence of fictitious jobsσ.Step 6: Take the first two fictitious jobs ofσ and allocate them on the machines at each stage with the EST rule with a better sequence which has less total energy consumption cost. Insert the sequence of the two fictitious jobs in to OP∗ and delete them from σ.Step 7: Consider the next fictitious job ofσ.Step 8: Insert the fictitious job into all positions of the partial sequenceOP∗ and find the sequence σ' with the less total energy consumption cost. Delete this fictitious job form σ.Step 9: SetOP∗ = σ′.Step 10: Ifσ=∅, go to Step 11, else go to Step 7.Step 11: Return the total energy consumption cost.Since ESPRA algorithms use the same insertion procedure as the algorithms EPRA and ESRA, the computational complexity of ESPRA algorithm isOn2Δ∑s=1kms. Since the ESPRA algorithm uses a fictitious job concept along with the insertion procedure, the computational complexity of this operation is On2Δ∑s=1kms/4. Thus, the effective size of the problem solved by ESPRA is only half size of the original problem. However, in big O notation, the constant (1/4) is avoided, and thus, the computational complexity of ESPRA algorithm is On2Δ∑s=1kms.


## 4.1. An Improved Algorithm Based on the NEH Algorithm (EPRA)

For the HFS scheduling problem, Ref. [37] proved that it is an NP problem even with two stages where one of the two stages has at least two machines. According to the complicated theorem, the problem HFS|rj,m1,…,mk|Emin is at least having that difficulty. There is no exact algorithm in polynomial time for even small size of NP problems. Therefore, heuristic algorithms are more effective for real manufacturing. Among these heuristic algorithms, constructive heuristic algorithms are designed for the problem and have their own merit (the advantage of computational time) compared with metaheuristic algorithms. However, the performances of constructive heuristic algorithms significantly depend on the different objectives and structures of the problems. Thus, for the new problems which even have a little difference with the existing problems, the existing constructive heuristic algorithms may be not effective. In the worst case, they cannot tackle those problems. Under such circumstances, existing heuristic algorithms should be developed or new algorithms should be proposed to solve the new problems. Ref. [38] proposed a heuristic algorithm called the NEH algorithm, which is probably considered as the best constructive heuristic method for the permutation flow shop makespan problem. The core ideal of the NEH algorithm is that a job with longer total processing time on the machine has the higher priority in the processing sequence. The steps of NEH algorithm are as follows (Algorithm 1):Algorithm 1: NEH algorithm.
Step 1: Jobs are sequenced by the decreasing order of their total processing time and obtain the sequence of jobs denoted bySStep 2: Take the first two jobs ofS as the partial schedule and find the optimal scheduleStep 3: Insert the next job ofS into the available positions of the previous partial schedule and find the partial schedule with the best objective function valueStep 4: Repeat the process for the remaining jobs according toSSince this paper considers the objective of minimization of energy consumption cost, we develop a new algorithm based on the NEH algorithm (called EPRA). The EPRA algorithm is similar to the NEH algorithm except for Step 1, i.e., constructing the sequence of jobs. In the EPRA algorithm, jobs are sequenced by the decreasing order of their summed running energy consumption cost based on the total processing times (denoted byepi) which is calculated by (12)epi=∑s=1k∑h=1mspihsEhs1ms,∀i,where ∑h=1mspihsEhs1/ms denotes the average running energy consumption cost based on the processing time at each stage. Since the NEH algorithm was originally proposed for the flow shop scheduling problem where each stage has only one machine, it does not have the problem that allocates the jobs to machines. For the problem FFS|rj,m1,…,mk|Emin addressed in this paper, the original NEH algorithm does not solve that problem. Thus, we propose a rule EST (the earliest starting time) to tackle that problem in Steps 2 and 3 of the EPRA algorithm. If there is one more machine with EST rule, then allocate the job to the machine that has the lowest running energy consumption cost.The major computation of EPRA algorithm has two operations: sorting the summed running energy consumption cost of jobs and inserting jobs to the optimal position of the partial sequence. The computational efficiency of the former operation isnlogn. The computational complexity of the second operation under the worst case is by (13). Since n2+n/2>nlogn, the computational complexity of EPRA algorithm is On2Δ∑s=1kms.(13)∑p=1np=∑s=1kms1+2+⋯+n=n2+n2∑s=1kms.In this paper, we develop two more new constructive heuristic algorithms based on the NEH algorithm with the objective of minimizing the cost of energy consumption. The remained algorithms are described in the following sections.


## 4.2. Energy Consumption Cost Based on Setup Time Ranking Algorithm (ESRA)

According to the above section, the NEH algorithm considers the processing time of jobs in constructing the sequence of jobs. However, one of the characteristics of the problem addressed in this paper is that there is setup time between two consecutive jobs and takes the minimization of total energy consumption as the objective. Therefore, we propose another heuristic algorithm called the ESRA algorithm which considers the energy consumption in the setup time of consecutive jobs and not the processing time of jobs. The fundamental principle behind this algorithm is that as the setup time between jobs is reduced, the cost of energy consumption also will get reduced. The ESRA algorithm begins with the summed energy consumption matrix (denoted by ESST) of the setup times of jobs at all stages. Thus, an elementesstij of ESST is obtained by adding the cost of energy consumption in the setup time of two consecutive jobs for all stages as follows:(14)esstij=∑g=1ST1∑h=1mgEhg2mg∑i=1ntijg1+∑e=1ST2∑h=1meEhe2me∑i=1n∑h=1metijhe2me,∀i,j,g∈ST1,e∈ST2Since the setup time has two types,esstij consists of two terms where the first term of equation (14) indicates the cost of idle energy consumption with type one setup time (dependent on the sequence of jobs) and the second term of equation (14) indicates the idle energy consumption cost with type two setup time (dependent on the sequence of jobs and the machine which the jobs are allocated to). In equation (14), we take the average unit per time of idle energy consumption cost of machines at each stage as the unit per time of idle energy consumption; that is, ∑h=1mgEhg2/mg and ∑h=1meEhe2/me respectively indicate the unit time of idle energy consumption cost with type one setup time and type two setup time.The ESST matrix gives a measure of idle energy consumption cost based on the setup time between jobs. The summed cost of each row in the ESST matrix can be used for ranking the combination of jobs from the lowest sum to the maximum sum to get the order of jobs S. The order S shows the best combination of jobs in terms of idle energy consumption cost based on setup time. The first job in S is the seed of the sequence. The complete sequence is built by adding the next lowest job to the best position in the partial sequence. Since the ESRA algorithm is similar to the EPRA algorithm, the computational complexity of ESRA algorithm isOn2Δ∑s=1kms (Algorithm 2).Algorithm 2: ESRA algorithm.
Step 1: According to equation (13), calculate the summed energy consumption cost based on the setup time of each job and order the jobs in an ascending manner . Thus, obtain the sequence of jobs S.Step 2: Take the first two jobs inS and schedule them in the best objective function. Delete the two jobs in S and save the best partial sequence OP.Step 3: FromO = 3 to n, take the O-th job in S and insert it into the O possible positions in OP. Then, select the best partial schedule. Finally, obtain the job processing order OP∗.Step 4: According toOP∗, schedule the jobs on the machines at each stage with EST rule.


## 4.3. Energy Consumption Cost Based on Setup Time and Processing Ranking Algorithm (ESPRA)

ESRA is an algorithm to predominantly consider the idle energy consumption cost based on setup time, and it does not concern the processing time for constructing the sequence. Conversely, the EPRA algorithm does not give importance to the setup time while constructing the sequence. In short, ESRA is a setup-time-based algorithm, and EPRA is a processing-time-based algorithm. In this paper, we propose an algorithm (denoted by ESPRA) by considering the setup time and processing time. The ESPRA algorithm considers the energy consumption based on setup time and processing time for constructing the sequence of jobs, giving importance to both of them. Similar to the ESRA algorithm, the ESPRA algorithm starts with the summed idle energy consumption matrix ESST based on setup time and then ranks the jobs by the ascending order of the summed row of ESST matrix to get an order of jobsS. The first job in S is the seed of the sequence. In the ESPRA algorithm, the sequence is constructed by using the concept of fictitious jobs instead of constructing the sequence by adding the next job in S on the possible positions of the seed or partial sequence. The concept of fictitious jobs is explained below.


### 4.3.1. Fictitious Jobs

The key characteristic of ESPRA algorithm is the concept of fictitious jobs. Fictitious jobs are defined by the best pair of jobs from S. It means the fictitious jobs have two jobs in them. The advantage of this concept is that it can keep the pairs which have the smallest idle energy consumption cost based on the setup time between them. The ESPRA algorithm takes this advantage which gives importance to the setup time. After identifying the fictitious jobs, the processing time of fictitious jobs at each stageFJis is calculated by equation (14); that is, the sum of the average running energy consumption cost based on the processing time of each pair of jobs at each stage:(15)FJis=∑h=1mspihsΔEhs1+∑h=1mspi+1hsΔEhs1ms,i∈S,where pn+1hs = 0. Figure 1 shows an illustration of ESPRA algorithm at stage s. J1, J2, …, J8 are individual jobs to be sequenced in the ascending order of the summed idle energy consumption based on setup time. FJ1s, FJ2s, FJ3s, and FJ4s are the fictitious jobs. For example, the fictitious job FJ1s includes the individual jobs J1 and J2. Then, the processing time of a fictitious jobs is the sum of the average running energy consumption cost of each pair of jobs instead of their processing times. After identifying the processing times of the fictitious jobs, these jobs are ranked by decreasing the order of them. The rank one job will be first added to the sequence, which forms the partial sequence. The next eligible candidate will be the job with the next highest summed running energy consumption cost based on processing time. All possible insertion positions are considered for inserting the new job in the partial sequence so that the total energy consumption cost is minimum. For example, in Figure 1, the fictitious job FJ1s with the largest summed running energy consumption cost is the partial sequence. FJ2s with the next largest summed running energy consumption cost has two possible points in the partial sequence, then inserts FJ2s into these two points, and then selects the position with less total energy consumption cost (assume that FJ2s is the front of FJ1s). The partial sequence is constructed by repeating the above operations until all fictitious jobs are completed (Algorithm 3).Figure 1
An illustration of ESPRA algorithm at stages.Algorithm 3: ESPRA algorithm.
Step 1: Set the partial sequenceOP∗ = 0.Step 2: Calculate the matrix ESST, where each entity of ESSTesstij is calculated by equation (13)Step 3: Jobs are sorted by the ascending order of the corresponding summed idle energy consumption cost based on the setup time in ESST to obtain the sequence of jobss.Step 4: According tos, identify the pairs of jobs as fictitious jobs, and calculate their running energy consumption cost based on processing time by equation (14).Step 5: Sort the fictitious jobs by the decreasing order of their running energy consumption cost and obtain the sequence of fictitious jobsσ.Step 6: Take the first two fictitious jobs ofσ and allocate them on the machines at each stage with the EST rule with a better sequence which has less total energy consumption cost. Insert the sequence of the two fictitious jobs in to OP∗ and delete them from σ.Step 7: Consider the next fictitious job ofσ.Step 8: Insert the fictitious job into all positions of the partial sequenceOP∗ and find the sequence σ' with the less total energy consumption cost. Delete this fictitious job form σ.Step 9: SetOP∗ = σ′.Step 10: Ifσ=∅, go to Step 11, else go to Step 7.Step 11: Return the total energy consumption cost.Since ESPRA algorithms use the same insertion procedure as the algorithms EPRA and ESRA, the computational complexity of ESPRA algorithm isOn2Δ∑s=1kms. Since the ESPRA algorithm uses a fictitious job concept along with the insertion procedure, the computational complexity of this operation is On2Δ∑s=1kms/4. Thus, the effective size of the problem solved by ESPRA is only half size of the original problem. However, in big O notation, the constant (1/4) is avoided, and thus, the computational complexity of ESPRA algorithm is On2Δ∑s=1kms.


## 4.3.1. Fictitious Jobs

The key characteristic of ESPRA algorithm is the concept of fictitious jobs. Fictitious jobs are defined by the best pair of jobs from S. It means the fictitious jobs have two jobs in them. The advantage of this concept is that it can keep the pairs which have the smallest idle energy consumption cost based on the setup time between them. The ESPRA algorithm takes this advantage which gives importance to the setup time. After identifying the fictitious jobs, the processing time of fictitious jobs at each stageFJis is calculated by equation (14); that is, the sum of the average running energy consumption cost based on the processing time of each pair of jobs at each stage:(15)FJis=∑h=1mspihsΔEhs1+∑h=1mspi+1hsΔEhs1ms,i∈S,where pn+1hs = 0. Figure 1 shows an illustration of ESPRA algorithm at stage s. J1, J2, …, J8 are individual jobs to be sequenced in the ascending order of the summed idle energy consumption based on setup time. FJ1s, FJ2s, FJ3s, and FJ4s are the fictitious jobs. For example, the fictitious job FJ1s includes the individual jobs J1 and J2. Then, the processing time of a fictitious jobs is the sum of the average running energy consumption cost of each pair of jobs instead of their processing times. After identifying the processing times of the fictitious jobs, these jobs are ranked by decreasing the order of them. The rank one job will be first added to the sequence, which forms the partial sequence. The next eligible candidate will be the job with the next highest summed running energy consumption cost based on processing time. All possible insertion positions are considered for inserting the new job in the partial sequence so that the total energy consumption cost is minimum. For example, in Figure 1, the fictitious job FJ1s with the largest summed running energy consumption cost is the partial sequence. FJ2s with the next largest summed running energy consumption cost has two possible points in the partial sequence, then inserts FJ2s into these two points, and then selects the position with less total energy consumption cost (assume that FJ2s is the front of FJ1s). The partial sequence is constructed by repeating the above operations until all fictitious jobs are completed (Algorithm 3).Figure 1
An illustration of ESPRA algorithm at stages.Algorithm 3: ESPRA algorithm.
Step 1: Set the partial sequenceOP∗ = 0.Step 2: Calculate the matrix ESST, where each entity of ESSTesstij is calculated by equation (13)Step 3: Jobs are sorted by the ascending order of the corresponding summed idle energy consumption cost based on the setup time in ESST to obtain the sequence of jobss.Step 4: According tos, identify the pairs of jobs as fictitious jobs, and calculate their running energy consumption cost based on processing time by equation (14).Step 5: Sort the fictitious jobs by the decreasing order of their running energy consumption cost and obtain the sequence of fictitious jobsσ.Step 6: Take the first two fictitious jobs ofσ and allocate them on the machines at each stage with the EST rule with a better sequence which has less total energy consumption cost. Insert the sequence of the two fictitious jobs in to OP∗ and delete them from σ.Step 7: Consider the next fictitious job ofσ.Step 8: Insert the fictitious job into all positions of the partial sequenceOP∗ and find the sequence σ' with the less total energy consumption cost. Delete this fictitious job form σ.Step 9: SetOP∗ = σ′.Step 10: Ifσ=∅, go to Step 11, else go to Step 7.Step 11: Return the total energy consumption cost.Since ESPRA algorithms use the same insertion procedure as the algorithms EPRA and ESRA, the computational complexity of ESPRA algorithm isOn2Δ∑s=1kms. Since the ESPRA algorithm uses a fictitious job concept along with the insertion procedure, the computational complexity of this operation is On2Δ∑s=1kms/4. Thus, the effective size of the problem solved by ESPRA is only half size of the original problem. However, in big O notation, the constant (1/4) is avoided, and thus, the computational complexity of ESPRA algorithm is On2Δ∑s=1kms.


## 5. Computational Experiments

### 5.1. Data Generation

In Section1, we know that there has been no literature on this problem addressed in this paper so far. Therefore, there is no existing benchmark instance for the problem HFS|rj,m1,…,mk|Emin, and we should design the benchmark instances for this problem in this section. Because the performance of the heuristic algorithm is partly influenced by the difference of benchmark instances, we should find a method which can objectively evaluate the efficiency of the three proposed heuristic algorithms. Based on the previous research, Logendran et al. [30] proposed a split-plot experimental design for the flexible flow shop scheduling problem, where the problem size is considered as the main factor. In this paper, the problem size is defined by the number of stages and the number of jobs (namely, k and n). Three levels of each of the two parameters are selected: 4, 6, and 8 for stages and 20, 50, and 100 for jobs. These selections contribute to identifying three different sizes—−small, medium, and large. Thus, small, medium, and large problem sizes are described by the following stages and jobs, respectively: 4 and 20, 6 and 50, and 8 and 100.Flexibility (denoted byF), setup-to-processing time proportion (denoted by S/P), and the algorithms used for tackling this problem are considered as the significant important factors in this paper. Thus, they are assigned to the subplot and deemed subplot factors in the split-plot design, as opposed to the problem size, that is the only main factor. The algorithms are described in Section 5. The flexibility in a flexible flow shop is defined by the number of machines at each stage. The larger the number of stages in the flexible flow shop and the larger the number of machines at each stage, the greater flexibility is in the shop. Similar to the problem size, the flexibility introduced in this paper has three levels, namely, low, medium, and high. In this paper, we determine the low flexibility as that is represented by 1/3 of stages having more than one unrelated parallel machine, medium flexibility by 2/3 stages having more than one unrelated parallel machine, and high flexibility by all stages having more than one unrelated parallel machine. For instance, for medium size and medium flexibility, the number of stages having more than one unrelated parallel machine is (2/3) ∗ (6) = 4. If the number of stages is odd, the rounded-up integer value is used. For the same instance, if the number of stage is 7, then there are around (2/3 ∗ 7 = 4.66) = 5 stages with more than one unrelated parallel machine. Once the number of stages having more than one unrelated parallel machine in a problem instance is identified, then the challenge is to determine which stages would have how many unrelated parallel machines. Uniformly distributed integer values (U [1, k], where k represents the number of stages) are generated for the stages with unrelated parallel machines, and the number of machines at the stages with unrelated parallel machines is generated by the integer uniform distribution U [2, 8]. For example, in the five stages with more than one unrelated parallel machine case, in a seven-stage problem, 4(3)-1(2)-3(6)-6(2)-2(7) is an acceptable value of generating random numbers where the right side of “()” represents the number stage with more than one unrelated parallel machine, and the value inside of “()” represents the number of machines in that stage but not 4 (3)-1(1)-1(6)-6(2)-2(7).The setup-to-processing time proportion is another subplot factor that needs to be described. In this paper, we use eight different proportions of setup time to processing time: 5%, 10%, 25%, 50%, 75%, 100%, 125%, and 150%. Corresponding to these proportions, setup times are generated by using uniform distributionU [1, 5], U [1, 10], U [1, 25], U [1, 50], U [1, 75], U [1, 100], U [1, 125], and U [1, 150], respectively. Thus, the processing times of jobs are generated by the setup time and the setup-to-processing time proportion. The levels of the factors in the experimental design are shown in Table 2. The settings of the rest parameters are generated as follows:(i)
The per unit time of running energy consumption cost and the per unit time of idle energy consumption cost of each machine are drawn by a random number from one to ten (namely,Ths1/Ths2 = randi (1, 10)).(ii)
The per unit time of warm-up energy consumption cost of each machine is generated byThs3=randi (3, 4) ∗ Ths2, where randi (3, 4) represents a random float number from three to four. Since it is meaningless, that if the per unit time of warm-up energy consumption cost of each machine is less than the per unit time of idle energy consumption cost of each machines, then Ths3>Ths2.(iii)
The warm-up time of each machine with the mode from stopping to operation is a function with the average processing time of all jobs on that machine and is generated byths=randi0.1,0.5∗∑j=1npjhs/n, where randi (0.1, 0.5) represents a random number from 0.1 to 0.5.Table 2
Parameter settings of the factors in the experimental design.
FactorsIndexLevelsMain factorProblem sizenkSubplot factor1204250631008FlexibilityFms11/3U [2, 8]22/331Setup-to-processing time proportionS/Ptijg1 (tijhe2)15%U [1, 5]210%U [1, 10]325%U [1, 25]450%U [1, 50]575%U [1, 75]6100%U [1, 100]7125%U [1, 125]8150%U [1, 150]Algorithm1EPRA2ESRA3ESPRAIn the data generated above, all level combinations result in 3 ∗ 3 ∗ 8 ∗ 3 = 216 configurations. Five instances are generated by Table1, and each instance is run by the three proposed algorithms, respectively, which means that we have to do 216 ∗ 5 = 1080 experiments. Then record the cost of total energy consumption and CPU running time for each heuristic algorithm in every instance. All heuristic algorithms are coded with MATLAB R2012a and implemented on an Intel (R) Core i7-3370 3.4 GHZ personal computer with 4 GB RAM.


### 5.2. Analysis of the Factors

In Section 6.1, we propose a split-plot design to investigate the research problem. Here, we present a linear model representative of the split-plot design as equation (16), which consists of five factors, including one for the number of replicates (i) used in the split-plot design. The remaining four are problem size (j), flexibility (k), setup-to-processing time proportion (l), and algorithm (m). The response variable representing the total cost of energy consumption is evaluated for the flexible flow shop scheduling problem and is thus denoted by yijklm:(16)yijklm=μ+ρi+βj+ε1+γk+δl+ηm+βδjl+βγjk+βηjm+γδkl+γηkm+δηlm+βγδjkl+βδηjlm+βγηjkm+γδηklm+βγδηjklm+ε2where μ is the overall mean, ρi is the treatment effect (replicate), βj is the treatment effect (problem size), γk is the treatment effect (flexibility), δl is the treatment effect (setup-to-processing time proportion), ηm is the treatment effect (algorithm), and other terms represent two-, three-, and four-factor interactions. ε1 and ε2 are whole-plot and subplot errors, respectively.The statistical analysis is coded in Statistical Analysis System (SAS) 9.1. Table3 is the summary of the results from ANOVA. From Table 3, the effect of replicates and the main and sub-plot errors have been combined and denoted as the error term in SAS output for a total sum of squares (SS) of 8.822 × 1011 and 864 degrees of freedom. The remaining SS of 3.378 × 1013 is contributed by the main and subplot factors and factor interactions representing the model. Observe that the main factor (problem size) and the two subplot factors (flexibility and setup-to-processing time proportion) are significant as well as their interaction terms at a significance level (α) of 5%, as the P-value less than (0.0001). However, the algorithm as a subplot factor is not significant and every other interaction term that involves the algorithm is not significant. That is reasonable because three algorithms are proposed based on the NEH algorithm. However, it does not mean that there is no difference in the performance of the three algorithms, and we will discuss the performance of the three algorithms in the next section. Thus, size ∗ flexibility ∗ proportion is the highest order interaction among the factors.Table 3
Results from ANOVA.
SourceDegrees of freedomSum of squares (SS)Mean squareF-valueModel2153.378 × 10131.571 × 1011153.88Error8648.822 × 10111.021 × 109Corrected total10793.466 × 1013Sizea − 1 = 28.639 × 10124.320 × 10124230.43Flexibilityb − 1 = 27.674 × 10103.837 × 101037.58Size ∗ flexibility(a − 1) (b − 1) = 42.288 × 10115.721×1010556.03Proportionc − 1 = 71.344 × 10131.919 × 10121879.66Size ∗ proportion(a − 1) (c − 1) = 141.111 × 10137.933 × 1011776.91Flexibility ∗ proportion(b − 1) (c − 1) = 141.065 × 10117.606 × 1097.45Size ∗ flexibility ∗ proportion(a − 1) (b − 1) (c − 1) = 281.713 × 10116.117 × 1095.99Algorithmd − 1 = 22.562 × 1091.281 × 1091.25Size ∗ algorithm(a − 1) (d − 1) = 41.849 × 1094.622 × 1080.45Flexibility ∗ algorithm(b − 1) (d − 1) = 43.642 × 1089.105 × 1070.09Size ∗ flexibility ∗ algorithm(a − 1) (b − 1) (d − 1) = 82.424 × 1083.030 × 1070.03Proportion ∗ algorithm(c − 1) (d − 1) = 144.034 × 1092.881 × 1080.28Size ∗ proportion ∗ algorithm(a − 1) (c − 1) (d − 1) = 283.308 × 1091.181 × 1080.12Flexibility ∗ proportion ∗ algorithm(b − 1) (c − 1) (d − 1) = 282.755 × 1099.838 × 1070.10Size ∗ flexibility ∗ proportion ∗ algorithm(a − 1) (b − 1) (c − 1) (d − 1) = 563.378 × 1096.032 × 1070.06Error(r − 1)abcd = 8648.822 × 10111.021 × 109Corrected totalabcdr − 1 = 10793.466 × 1013As the size ∗ flexibility ∗ proportion is considered as the highest order interaction, Table4 shows the test results from effect slices for the effect corresponding to size ∗ flexibility ∗ proportion to analyze the difference within these three factors. As we can see, a statistically significant difference in the performance of the proportions is observed when medium and large problem sizes with all levels of flexibility are considered. For the remaining levels of size and flexibility, the performance of the eight proportions does not record a statistically significant difference.Table 4
Test of effect slices.
EffectSizeFlexibilityNumber of degrees of freedomDenominator degrees of freedomF-valuePr >FSize ∗ flexibility ∗ proportionSmallLow78640.011.0000Size ∗ flexibility ∗ proportionSmallMedium78640.830.5632Size ∗ flexibility ∗ proportionSmallHigh78643.300.0018Size ∗ flexibility ∗ proportionMediumLow7864107.16<0.0001Size ∗ flexibility ∗ proportionMediumMedium7864130.14<0.0001Size ∗ flexibility ∗ proportionMediumHigh786490.95<0.0001Size ∗ flexibility ∗ proportionLargeLow7864998.76<0.0001Size ∗ flexibility ∗ proportionLargeMedium78641311.04<0.0001Size ∗ flexibility ∗ proportionLargeHigh7864830.16<0.0001


### 5.3. Analysis of the Algorithms

Since the factor of the algorithm is not significant in Section5.2, we test the performance of the three algorithms in this section. From Section 5.2, we know that it is a statistically significant difference in the performance of the proportions with medium and large problem sizes in all levels of flexibility. Thus, we compare the three heuristic algorithms under the configuration of medium and large problem sizes with all levels of flexibility and setup-to-processing time proportion. To evaluate the performance of the three algorithms, we define a performance index (dev) by (17)dev=EFFS−maxLB1,LB2maxLB1,LB2×100%,where EFFS represents the value of the objective, and LB1 and LB2 represent the low bounds for the HFSri,m1,…,mkEmin problem. dev represents the deviation from the objective value to the lower bound. The smaller the dev is, the better the heuristic algorithm will be. The results of the comparison of the heuristic algorithms are shown in Tables 5 and 6.Table 5
The comparison of the heuristic algorithms (n = 50, k = 6).
FlexibilityAlgorithmProportion5%10%25%50%75%100%125%150%1/3EPRA68.6379.5287.8287.4860.9966.39154.4161.62ESRA86.0672.1481.6181.6757.2766.26147.0158.54ESPRA66.6971.7580.0480.1855.7366.39144.6558.222/3EPRA197.10231.48234.40263.35280.77330.91261.98181.16ESRA184.72221.23225.64246.83272.14317.89251.64163.08ESPRA183.76221.00301.42247.04266.09317.69246.36168.313/3EPRA226.18242.09388.69504.05433.54595.67564.66547.23ESRA218.60234.47400.94479.94405.95525.45533.70511.27ESPRA204.85229.72384.75477.50389.70578.40534.05474.46Table 6
The comparison of the heuristic algorithms (n = 100, k = 8).
FlexibilityAlgorithmProportion5%10%25%50%75%100%125%150%1/3EPRA99.69126.11117.75138.82120.58120.29119.07119.53ESRA94.67118.94111.59132.77114.96113.58115.04116.14ESPRA94.53118.50110.68134.67113.73113.31114.72115.482/3EPRA344.02308.15449.24379.20286.80455.58386.86474.91ESRA332.47296.91435.60343.17278.89457.72370.81456.32ESPRA328.68300.59431.17366.55274.60457.72370.39455.223/3EPRA436.58585.29511.22696.39736.40596.19527.52722.97ESRA423.30557.84497.17665.45695.38609.89505.36704.21ESPRA396.55550.08496.51690.63713.03592.27497.80714.56Each entity under every column in Tables5 and 6 is the average dev for the five instances run by one algorithm with the same configuration, where the best values among these algorithms are marked by bold. From Tables 5 and 6, we have the conclusions as follows:(1)
With the increase of the flexibility,dev of the three algorithms becomes larger, which means that the performance of the three algorithms deviates further from the optimal solution with the increase of the flexibility. That is because the space of feasible solutions becomes larger with the increasing flexibility, and it makes the three algorithms more difficult to find a better solution. From Tables 5 and 6, it is clearly found that dev is the least when the flexibility is low. In addition, the dev becomes larger with the increase of problem size because of the same reason.(2)
As discussed in Section5.2, the algorithm is not significant. We find that the performance of the three algorithms has no significant difference within the medium and large problem sizes, where the largest discrepancy of dev (75.61%) occurs between the ESRA algorithm and the ESPRA algorithm when the flexibility is medium, and setup-to-processing time proportion is 25%, and the problem size is medium. However, it is natural that the ESPRA algorithm is generally the best among the three algorithms (24.3% of all configurations in Tables 5 and 6) since the ESPRA algorithm simultaneously considers the setup time and the processing time. In addition, the superiority of ESPRA algorithm is more obvious with the increase of problem size. Although the performance of ESRA algorithm is very close to ESPRA algorithm’s, the running time of ESRA algorithm is greatly larger than that of the ESPRA algorithm. The details of the running time of the three algorithms will be discussed in Section 5.4. The performance of EPRA is the worst among the three algorithms.


### 5.4. Analysis of the Lower Bounds

Tables7–9 show the results of the comparisons of the lower bounds. Each entity under every column in Tables 7–9 has two lower bounds which are proposed in Section 4. Recall that LB1 is job-based bound, giving the minimal energy consumption cost over all jobs in all stages. LB2 is machine-based bound, giving the minimal value of running energy consumption over all stages to pre-emptively process and with minimal idle energy consumption based on the setup time in that stage. Each entity in every row of Tables 7–9 is the number of better ones of the lower bound under the same configuration. While LB1 is larger with 798 out of 1080 times, LB2 is larger with 282 times. Thus, the lower bound LB1 is generally the best when compared to LB2. However, LB1 is not always the best. When the problem size is small (especially flexibility is high), LB2 is better than LB1. In addition, we note that the gap performance of lower bounds between LB1 and LB2 is not large with the setup-to-processing time proportions from 5% to 50% (LB1 larges 353 times, LB2 larges 187 times), but this gap becomes larger with the setup-to-processing time proportions from 50% to 150% (LB1 larges 445 times, LB2 larges 95 times). Figure 2 is the comparison of lower bounds with different flexibilities. From Figure 2, we find that LB2 is better than LB1 when the flexibility is high.Table 7
The comparison of the lower bounds (n = 20, k = 4).
FlexibilityLower boundProportion5%10%25%50%75%100%125%150%1/3LB114109615131215LB2156902302/3LB1670810151515LB29815750003/3LB103000301LB21512151515121514Table 8
The comparison of the lower bounds (n = 20, k = 4).
FlexibilityLower boundProportion5%10%25%50%75%100%125%150%1/3LB11515151515151515LB2000000002/3LB11515151515151515LB2000000003/3LB162371011129LB29131285436Table 9
The comparison of the lower bounds (n = 100, k = 8).
FlexibilityLower boundProportion5%10%25%50%75%100%125%150%1/3LB11515151515151515LB2000000002/3LB1912151515151515LB2630000003/3LB18159910121512LB270665303Figure 2
The comparisons of the lower bounds with different flexibilities.


### 5.5. Analysis of the CPU Running Time

Figure3 shows the average CPU running time for each algorithm in different problem sizes. It is clearly found that the average CPU running time of all algorithms is larger with the increase of problem size. Especially the increase of CPU running time from medium problem size to large problem size is greatly larger than that from small problem size to medium size. As the time complexity of the algorithms mentioned in Section 5, it is an exponential function related to n. Thus, the CPU running time of all algorithms exponentially grows with the increasing problem size. In addition, the CPU running time of ESPRA algorithm is greatly less than EPRA and ESRA algorithms. That is because we define a concept of fictitious job which reduces the number of insertions in the inserting operation of ESPRA algorithm. The CPU running time of EPRA algorithm and ESRA algorithm has no significant difference. Moreover, for the ESPRA algorithm, there is an interesting phenomenon: the CPU running time has a salient point when setup-to-processing time proportion is 100% in the medium and larger problem sizes. However, the CPU running time of EPRA algorithm and ESRA algorithm has no significant difference with the setup-to-processing time proportions in all levels of problem size.Figure 3
The comparison of the three algorithms. (a)n = 20, k = 4. (b) n = 50, k = 6. (c) n = 100, k = 8.
(a)(b)(c)


## 5.1. Data Generation

In Section1, we know that there has been no literature on this problem addressed in this paper so far. Therefore, there is no existing benchmark instance for the problem HFS|rj,m1,…,mk|Emin, and we should design the benchmark instances for this problem in this section. Because the performance of the heuristic algorithm is partly influenced by the difference of benchmark instances, we should find a method which can objectively evaluate the efficiency of the three proposed heuristic algorithms. Based on the previous research, Logendran et al. [30] proposed a split-plot experimental design for the flexible flow shop scheduling problem, where the problem size is considered as the main factor. In this paper, the problem size is defined by the number of stages and the number of jobs (namely, k and n). Three levels of each of the two parameters are selected: 4, 6, and 8 for stages and 20, 50, and 100 for jobs. These selections contribute to identifying three different sizes—−small, medium, and large. Thus, small, medium, and large problem sizes are described by the following stages and jobs, respectively: 4 and 20, 6 and 50, and 8 and 100.Flexibility (denoted byF), setup-to-processing time proportion (denoted by S/P), and the algorithms used for tackling this problem are considered as the significant important factors in this paper. Thus, they are assigned to the subplot and deemed subplot factors in the split-plot design, as opposed to the problem size, that is the only main factor. The algorithms are described in Section 5. The flexibility in a flexible flow shop is defined by the number of machines at each stage. The larger the number of stages in the flexible flow shop and the larger the number of machines at each stage, the greater flexibility is in the shop. Similar to the problem size, the flexibility introduced in this paper has three levels, namely, low, medium, and high. In this paper, we determine the low flexibility as that is represented by 1/3 of stages having more than one unrelated parallel machine, medium flexibility by 2/3 stages having more than one unrelated parallel machine, and high flexibility by all stages having more than one unrelated parallel machine. For instance, for medium size and medium flexibility, the number of stages having more than one unrelated parallel machine is (2/3) ∗ (6) = 4. If the number of stages is odd, the rounded-up integer value is used. For the same instance, if the number of stage is 7, then there are around (2/3 ∗ 7 = 4.66) = 5 stages with more than one unrelated parallel machine. Once the number of stages having more than one unrelated parallel machine in a problem instance is identified, then the challenge is to determine which stages would have how many unrelated parallel machines. Uniformly distributed integer values (U [1, k], where k represents the number of stages) are generated for the stages with unrelated parallel machines, and the number of machines at the stages with unrelated parallel machines is generated by the integer uniform distribution U [2, 8]. For example, in the five stages with more than one unrelated parallel machine case, in a seven-stage problem, 4(3)-1(2)-3(6)-6(2)-2(7) is an acceptable value of generating random numbers where the right side of “()” represents the number stage with more than one unrelated parallel machine, and the value inside of “()” represents the number of machines in that stage but not 4 (3)-1(1)-1(6)-6(2)-2(7).The setup-to-processing time proportion is another subplot factor that needs to be described. In this paper, we use eight different proportions of setup time to processing time: 5%, 10%, 25%, 50%, 75%, 100%, 125%, and 150%. Corresponding to these proportions, setup times are generated by using uniform distributionU [1, 5], U [1, 10], U [1, 25], U [1, 50], U [1, 75], U [1, 100], U [1, 125], and U [1, 150], respectively. Thus, the processing times of jobs are generated by the setup time and the setup-to-processing time proportion. The levels of the factors in the experimental design are shown in Table 2. The settings of the rest parameters are generated as follows:(i)
The per unit time of running energy consumption cost and the per unit time of idle energy consumption cost of each machine are drawn by a random number from one to ten (namely,Ths1/Ths2 = randi (1, 10)).(ii)
The per unit time of warm-up energy consumption cost of each machine is generated byThs3=randi (3, 4) ∗ Ths2, where randi (3, 4) represents a random float number from three to four. Since it is meaningless, that if the per unit time of warm-up energy consumption cost of each machine is less than the per unit time of idle energy consumption cost of each machines, then Ths3>Ths2.(iii)
The warm-up time of each machine with the mode from stopping to operation is a function with the average processing time of all jobs on that machine and is generated byths=randi0.1,0.5∗∑j=1npjhs/n, where randi (0.1, 0.5) represents a random number from 0.1 to 0.5.Table 2
Parameter settings of the factors in the experimental design.
FactorsIndexLevelsMain factorProblem sizenkSubplot factor1204250631008FlexibilityFms11/3U [2, 8]22/331Setup-to-processing time proportionS/Ptijg1 (tijhe2)15%U [1, 5]210%U [1, 10]325%U [1, 25]450%U [1, 50]575%U [1, 75]6100%U [1, 100]7125%U [1, 125]8150%U [1, 150]Algorithm1EPRA2ESRA3ESPRAIn the data generated above, all level combinations result in 3 ∗ 3 ∗ 8 ∗ 3 = 216 configurations. Five instances are generated by Table1, and each instance is run by the three proposed algorithms, respectively, which means that we have to do 216 ∗ 5 = 1080 experiments. Then record the cost of total energy consumption and CPU running time for each heuristic algorithm in every instance. All heuristic algorithms are coded with MATLAB R2012a and implemented on an Intel (R) Core i7-3370 3.4 GHZ personal computer with 4 GB RAM.


## 5.2. Analysis of the Factors

In Section 6.1, we propose a split-plot design to investigate the research problem. Here, we present a linear model representative of the split-plot design as equation (16), which consists of five factors, including one for the number of replicates (i) used in the split-plot design. The remaining four are problem size (j), flexibility (k), setup-to-processing time proportion (l), and algorithm (m). The response variable representing the total cost of energy consumption is evaluated for the flexible flow shop scheduling problem and is thus denoted by yijklm:(16)yijklm=μ+ρi+βj+ε1+γk+δl+ηm+βδjl+βγjk+βηjm+γδkl+γηkm+δηlm+βγδjkl+βδηjlm+βγηjkm+γδηklm+βγδηjklm+ε2where μ is the overall mean, ρi is the treatment effect (replicate), βj is the treatment effect (problem size), γk is the treatment effect (flexibility), δl is the treatment effect (setup-to-processing time proportion), ηm is the treatment effect (algorithm), and other terms represent two-, three-, and four-factor interactions. ε1 and ε2 are whole-plot and subplot errors, respectively.The statistical analysis is coded in Statistical Analysis System (SAS) 9.1. Table3 is the summary of the results from ANOVA. From Table 3, the effect of replicates and the main and sub-plot errors have been combined and denoted as the error term in SAS output for a total sum of squares (SS) of 8.822 × 1011 and 864 degrees of freedom. The remaining SS of 3.378 × 1013 is contributed by the main and subplot factors and factor interactions representing the model. Observe that the main factor (problem size) and the two subplot factors (flexibility and setup-to-processing time proportion) are significant as well as their interaction terms at a significance level (α) of 5%, as the P-value less than (0.0001). However, the algorithm as a subplot factor is not significant and every other interaction term that involves the algorithm is not significant. That is reasonable because three algorithms are proposed based on the NEH algorithm. However, it does not mean that there is no difference in the performance of the three algorithms, and we will discuss the performance of the three algorithms in the next section. Thus, size ∗ flexibility ∗ proportion is the highest order interaction among the factors.Table 3
Results from ANOVA.
SourceDegrees of freedomSum of squares (SS)Mean squareF-valueModel2153.378 × 10131.571 × 1011153.88Error8648.822 × 10111.021 × 109Corrected total10793.466 × 1013Sizea − 1 = 28.639 × 10124.320 × 10124230.43Flexibilityb − 1 = 27.674 × 10103.837 × 101037.58Size ∗ flexibility(a − 1) (b − 1) = 42.288 × 10115.721×1010556.03Proportionc − 1 = 71.344 × 10131.919 × 10121879.66Size ∗ proportion(a − 1) (c − 1) = 141.111 × 10137.933 × 1011776.91Flexibility ∗ proportion(b − 1) (c − 1) = 141.065 × 10117.606 × 1097.45Size ∗ flexibility ∗ proportion(a − 1) (b − 1) (c − 1) = 281.713 × 10116.117 × 1095.99Algorithmd − 1 = 22.562 × 1091.281 × 1091.25Size ∗ algorithm(a − 1) (d − 1) = 41.849 × 1094.622 × 1080.45Flexibility ∗ algorithm(b − 1) (d − 1) = 43.642 × 1089.105 × 1070.09Size ∗ flexibility ∗ algorithm(a − 1) (b − 1) (d − 1) = 82.424 × 1083.030 × 1070.03Proportion ∗ algorithm(c − 1) (d − 1) = 144.034 × 1092.881 × 1080.28Size ∗ proportion ∗ algorithm(a − 1) (c − 1) (d − 1) = 283.308 × 1091.181 × 1080.12Flexibility ∗ proportion ∗ algorithm(b − 1) (c − 1) (d − 1) = 282.755 × 1099.838 × 1070.10Size ∗ flexibility ∗ proportion ∗ algorithm(a − 1) (b − 1) (c − 1) (d − 1) = 563.378 × 1096.032 × 1070.06Error(r − 1)abcd = 8648.822 × 10111.021 × 109Corrected totalabcdr − 1 = 10793.466 × 1013As the size ∗ flexibility ∗ proportion is considered as the highest order interaction, Table4 shows the test results from effect slices for the effect corresponding to size ∗ flexibility ∗ proportion to analyze the difference within these three factors. As we can see, a statistically significant difference in the performance of the proportions is observed when medium and large problem sizes with all levels of flexibility are considered. For the remaining levels of size and flexibility, the performance of the eight proportions does not record a statistically significant difference.Table 4
Test of effect slices.
EffectSizeFlexibilityNumber of degrees of freedomDenominator degrees of freedomF-valuePr >FSize ∗ flexibility ∗ proportionSmallLow78640.011.0000Size ∗ flexibility ∗ proportionSmallMedium78640.830.5632Size ∗ flexibility ∗ proportionSmallHigh78643.300.0018Size ∗ flexibility ∗ proportionMediumLow7864107.16<0.0001Size ∗ flexibility ∗ proportionMediumMedium7864130.14<0.0001Size ∗ flexibility ∗ proportionMediumHigh786490.95<0.0001Size ∗ flexibility ∗ proportionLargeLow7864998.76<0.0001Size ∗ flexibility ∗ proportionLargeMedium78641311.04<0.0001Size ∗ flexibility ∗ proportionLargeHigh7864830.16<0.0001


## 5.3. Analysis of the Algorithms

Since the factor of the algorithm is not significant in Section5.2, we test the performance of the three algorithms in this section. From Section 5.2, we know that it is a statistically significant difference in the performance of the proportions with medium and large problem sizes in all levels of flexibility. Thus, we compare the three heuristic algorithms under the configuration of medium and large problem sizes with all levels of flexibility and setup-to-processing time proportion. To evaluate the performance of the three algorithms, we define a performance index (dev) by (17)dev=EFFS−maxLB1,LB2maxLB1,LB2×100%,where EFFS represents the value of the objective, and LB1 and LB2 represent the low bounds for the HFSri,m1,…,mkEmin problem. dev represents the deviation from the objective value to the lower bound. The smaller the dev is, the better the heuristic algorithm will be. The results of the comparison of the heuristic algorithms are shown in Tables 5 and 6.Table 5
The comparison of the heuristic algorithms (n = 50, k = 6).
FlexibilityAlgorithmProportion5%10%25%50%75%100%125%150%1/3EPRA68.6379.5287.8287.4860.9966.39154.4161.62ESRA86.0672.1481.6181.6757.2766.26147.0158.54ESPRA66.6971.7580.0480.1855.7366.39144.6558.222/3EPRA197.10231.48234.40263.35280.77330.91261.98181.16ESRA184.72221.23225.64246.83272.14317.89251.64163.08ESPRA183.76221.00301.42247.04266.09317.69246.36168.313/3EPRA226.18242.09388.69504.05433.54595.67564.66547.23ESRA218.60234.47400.94479.94405.95525.45533.70511.27ESPRA204.85229.72384.75477.50389.70578.40534.05474.46Table 6
The comparison of the heuristic algorithms (n = 100, k = 8).
FlexibilityAlgorithmProportion5%10%25%50%75%100%125%150%1/3EPRA99.69126.11117.75138.82120.58120.29119.07119.53ESRA94.67118.94111.59132.77114.96113.58115.04116.14ESPRA94.53118.50110.68134.67113.73113.31114.72115.482/3EPRA344.02308.15449.24379.20286.80455.58386.86474.91ESRA332.47296.91435.60343.17278.89457.72370.81456.32ESPRA328.68300.59431.17366.55274.60457.72370.39455.223/3EPRA436.58585.29511.22696.39736.40596.19527.52722.97ESRA423.30557.84497.17665.45695.38609.89505.36704.21ESPRA396.55550.08496.51690.63713.03592.27497.80714.56Each entity under every column in Tables5 and 6 is the average dev for the five instances run by one algorithm with the same configuration, where the best values among these algorithms are marked by bold. From Tables 5 and 6, we have the conclusions as follows:(1)
With the increase of the flexibility,dev of the three algorithms becomes larger, which means that the performance of the three algorithms deviates further from the optimal solution with the increase of the flexibility. That is because the space of feasible solutions becomes larger with the increasing flexibility, and it makes the three algorithms more difficult to find a better solution. From Tables 5 and 6, it is clearly found that dev is the least when the flexibility is low. In addition, the dev becomes larger with the increase of problem size because of the same reason.(2)
As discussed in Section5.2, the algorithm is not significant. We find that the performance of the three algorithms has no significant difference within the medium and large problem sizes, where the largest discrepancy of dev (75.61%) occurs between the ESRA algorithm and the ESPRA algorithm when the flexibility is medium, and setup-to-processing time proportion is 25%, and the problem size is medium. However, it is natural that the ESPRA algorithm is generally the best among the three algorithms (24.3% of all configurations in Tables 5 and 6) since the ESPRA algorithm simultaneously considers the setup time and the processing time. In addition, the superiority of ESPRA algorithm is more obvious with the increase of problem size. Although the performance of ESRA algorithm is very close to ESPRA algorithm’s, the running time of ESRA algorithm is greatly larger than that of the ESPRA algorithm. The details of the running time of the three algorithms will be discussed in Section 5.4. The performance of EPRA is the worst among the three algorithms.


## 5.4. Analysis of the Lower Bounds

Tables7–9 show the results of the comparisons of the lower bounds. Each entity under every column in Tables 7–9 has two lower bounds which are proposed in Section 4. Recall that LB1 is job-based bound, giving the minimal energy consumption cost over all jobs in all stages. LB2 is machine-based bound, giving the minimal value of running energy consumption over all stages to pre-emptively process and with minimal idle energy consumption based on the setup time in that stage. Each entity in every row of Tables 7–9 is the number of better ones of the lower bound under the same configuration. While LB1 is larger with 798 out of 1080 times, LB2 is larger with 282 times. Thus, the lower bound LB1 is generally the best when compared to LB2. However, LB1 is not always the best. When the problem size is small (especially flexibility is high), LB2 is better than LB1. In addition, we note that the gap performance of lower bounds between LB1 and LB2 is not large with the setup-to-processing time proportions from 5% to 50% (LB1 larges 353 times, LB2 larges 187 times), but this gap becomes larger with the setup-to-processing time proportions from 50% to 150% (LB1 larges 445 times, LB2 larges 95 times). Figure 2 is the comparison of lower bounds with different flexibilities. From Figure 2, we find that LB2 is better than LB1 when the flexibility is high.Table 7
The comparison of the lower bounds (n = 20, k = 4).
FlexibilityLower boundProportion5%10%25%50%75%100%125%150%1/3LB114109615131215LB2156902302/3LB1670810151515LB29815750003/3LB103000301LB21512151515121514Table 8
The comparison of the lower bounds (n = 20, k = 4).
FlexibilityLower boundProportion5%10%25%50%75%100%125%150%1/3LB11515151515151515LB2000000002/3LB11515151515151515LB2000000003/3LB162371011129LB29131285436Table 9
The comparison of the lower bounds (n = 100, k = 8).
FlexibilityLower boundProportion5%10%25%50%75%100%125%150%1/3LB11515151515151515LB2000000002/3LB1912151515151515LB2630000003/3LB18159910121512LB270665303Figure 2
The comparisons of the lower bounds with different flexibilities.


## 5.5. Analysis of the CPU Running Time

Figure3 shows the average CPU running time for each algorithm in different problem sizes. It is clearly found that the average CPU running time of all algorithms is larger with the increase of problem size. Especially the increase of CPU running time from medium problem size to large problem size is greatly larger than that from small problem size to medium size. As the time complexity of the algorithms mentioned in Section 5, it is an exponential function related to n. Thus, the CPU running time of all algorithms exponentially grows with the increasing problem size. In addition, the CPU running time of ESPRA algorithm is greatly less than EPRA and ESRA algorithms. That is because we define a concept of fictitious job which reduces the number of insertions in the inserting operation of ESPRA algorithm. The CPU running time of EPRA algorithm and ESRA algorithm has no significant difference. Moreover, for the ESPRA algorithm, there is an interesting phenomenon: the CPU running time has a salient point when setup-to-processing time proportion is 100% in the medium and larger problem sizes. However, the CPU running time of EPRA algorithm and ESRA algorithm has no significant difference with the setup-to-processing time proportions in all levels of problem size.Figure 3
The comparison of the three algorithms. (a)n = 20, k = 4. (b) n = 50, k = 6. (c) n = 100, k = 8.
(a)(b)(c)


## 6. Results and Discussion

This paper is an important addition to the literature where the majority of scheduling studies are focused on production efficiency while ignoring vital energy-related issues. In this paper, we investigate two different types of sequence-dependent setup times, namely, setup time depending on the sequence of jobs and setup time depending on both the sequence of jobs and the machines where the jobs are allocated, in a flexible flow shop scheduling problem where one or more stages may consist of two or more unrelated parallel machines. Three new constructive heuristic algorithms (EPRA, ESRA, and ESPRA) are proposed based on the NEH algorithm. Two lower bounds have been used in the evaluation of these heuristic algorithms. Three different problem structures, ranging in size from small or medium to large, in conjunction with factors such as flexibility, setup-to-processing time, proportion, and algorithms have been considered to develop a 3 ∗ 3 ∗ 8 ∗ 3 = 216 split-plot design. Five instances are generated for each configuration. Factor analysis, algorithm analysis, lower bounds analysis, and CPU running time analysis are carried out.Factor analysis reveals that all main factor interactions except the algorithm are significant. The algorithm as a subplot factor, interactions with other subplot factors, and high-order interactions are not significant. It means that in the flexible flow shop scheduling problem addressed in this paper, the algorithm factor is not an important factor. Because all three algorithms are proposed based on the NEH algorithm, and it is reasonable that the algorithm factor has no significant difference in the results of ANOVA. However, it does not mean that there is no difference in the performance of the three algorithms.Algorithm analysis shows thatdev of the three algorithms becomes larger with the increase of flexibility. The performance of ESPRA algorithm is in general the best among the three algorithms. Although the performance of ESRA algorithm is very close to that of the ESPRA algorithm, it consumes greatly more time than the ESPRA algorithm.Lower bounds analysis shows thatLB1 is generally better than LB2. However, LB2 is better when the flexibility is high. Finally, the CPU running time of the algorithms shows that the CPU running time of ESPRA algorithm is greatly less than EPRA and ESRA algorithms, and there is no significant difference of CPU running time between the EPRA algorithm and ESRA algorithm. In addition, an interesting phenomenon is observed in the CPU running time analysis that there is a salient point when setup-to-processing time proportion is 100% in the medium and larger problem sizes.Several issues are worthy of future investigations are shown as follows: (1) it is meaningful for the theoretical analysis of the proposed heuristic algorithms, e.g., worst case analysis and (2) it is interesting to extend the proposed heuristic algorithms to deal with other performance criteria.



---

*Source: 8682203-2022-09-22.xml*