A memetic algorithm for energy-efficient distributed re-entrant hybrid flow shop scheduling problem

Abstract

Facing the worsening environmental problems, green manufacturing and sustainable development have attracted much attention. Aiming at the energy-efficient distributed re-entrant hybrid flow shop scheduling problem considering the customer order constraints (EDORHFSP) under Time-of-Use (TOU) electricity price, a mathematical model is established to minimize the maximum completion time and total consumption energy cost. In the study, some customer orders require production in multiple factories and jobs belonging to the same customer order must be processed in one factory. Firstly, a memetic algorithm (MA) was proposed to solve the problem. To improve the performance of the algorithm, encoding and decoding methods, energy cost saving procedure, three heuristic rules about the population initialization and some neighborhood search methods are designed. Then, Taguchi method is adopted to research the influence of parameters setting. Lastly, numerical experiments demonstrate the effectiveness and superiority of MA for the EDORHFSP.

Keywords

Energy-efficient memetic algorithm Time-of-Use electricity price distributed re-entrant hybrid flow shop scheduling customer order constraints

1 Introduction

The development of manufacturing industry has brought great convenience to people’s life, but also led to some new problems, such as environmental pollution, global warming, energy depletion and so on. As rapid industrialization of the world, the constraint of resources and environment has become the bottleneck that restricts the sustainable development of industry. Meanwhile, under the background of economic globalization and cooperative production, manufacturing has changed from the traditional single factory production to multi-location factoriesproduction and distributed manufacturing has become a common production mode. Distributed manufacturing can make full use of the resources, reduce the cost and improve the competitiveness of enterprises through reasonable allocation, optimal combination and resource sharing.

Re-entrant hybrid flow shop scheduling problem (RHFSP), which has been proved to be a NP-hard problem [1], is an extension of classical hybrid flow shop scheduling problem (HFSP) and is widely used in chemical, textile, steel, semiconductor [2] and other industries. Therefore, the study of RHFSP has important theoretical significance and practical value. Since Graves and Meal et al. [3] proposed the re-entrant scheduling problem in 1983, more and more scholars have studied RHFSP in recent years. Kim and Lee [4] considered the unrelated parallel machine based on the general RHFSP with the objectives of minimizing the maximum completion time and the total tardiness time. But in fact, the tardiness time objective is treated as a constraint condition. Dugardin and Yalaoui et al. [5] proposed a new multi-objective genetic algorithm to solve the RHFSP, aiming at minimizing the maximum completion time and maximizing utilization of bottleneck machine, simultaneously. Cho and Bae et al. [6], Ying and Lin et al. [7], and Shen and Wang et al. [8, 9] proposed the pareto genetic algorithm, iterative pareto greedy algorithm, improved teaching and learning optimization algorithm and pareto based discrete harmony search algorithm respectively to solve the bi-objective RHFSP with the objectives of minimizing the maximum completion time and the total tardiness. Lin and Liu et al. [10] designed a hybrid harmony search and genetic algorithm for the RHFSP considering limited buffer capacity, aiming at minimizing the maximum completion time and mean flowtime. Zhou and Hu et al. [101] proposed a hybrid differential evolution algorithm to solve the RHFSP considering inspection and repair operations. The objective is minimizing total weighted completion time. To sum up, the optimization objectives are mainly time-related in previous researches, and rarely consider the environmental and energy cost related objectives [8].

According to statistics, by 2050, global energy demand will increase by 37%, and CO2 emissions will increase by about 80% in the next 30–40 years [11]. As we all know, the industrial sector consumes the most energy around the word [12, 13]. Thus, how to save energy has become an urgent problem in modern industrial production. Developing energy-saving equipment is a direct way to save energy, but it is difficult to achieve the target of energy saving and economic benefits in a short time due to large investment and long cycle [14]. In addition, green scheduling is another means to achieve energy saving and consumption reduction. It can improve energy utilization at a very small cost through reasonable optimization of resource allocation, operation sequencing and operation mode [15]. In recent years, more and more attention has been paid to green scheduling problems, in which several mechanisms are mainly considered, such as power off mechanism [12, 16], dynamic machine speed [17 –19] and TOU electricity price mechanism [19 –21]. To reduce the burden on the grid in peak periods, TOU electricity price has been implemented in many countries to encourage consumers to shift the usage from peak periods to off-peak or mid-peak periods. Sin et al. [22] investigated the scheduling problem considering preventive maintenance under TOU electricity price. Ding and Schulz et al. [22] developed a hybrid particle swarm optimization algorithm to solve the flexible flow shop scheduling considering TOU electricity price. Zhou and Jin et al. [23] used a hybrid multi-objective meta-heuristic algorithm to study the energy-efficient scheduling problem considering dynamic job arrival times under TOU electricity price. Zheng and Zhou [24] proposed a hybrid ant colony optimization algorithm to solve the two-stage blocking green permutation flow shop scheduling problem under TOU electricity price. In the paper, we study the EDORHFSP under the TOU electricity price mechanism.

Most of the previous studies on scheduling problem assume that all the jobs are processed in a single factory. Nevertheless, in today’s highly competitive marketplace, distributed scheduling has become increasingly popular owing to its production efficiency and efficient resource utilization [25]. Distributed scheduling problem (DSP) mainly involves two sub-problem: assigning jobs to factories and determining jobs sequence. Gao et al. [26] used an enhanced tabu search algorithm for the distributed permutation flow shop scheduling problem (DPFSP). Lin et al. [27] adopted an improved iterated greedy algorithm to solve the DPFSP. Wang et al. [14] prosed a knowledge-based cooperative algorithm to solve the energy-efficient DPFSP. It is not difficult to find that in the classic DSP, jobs are processed as independent entities. However, several jobs belong to the same customer order in the actual production scenario. In that case, it is wise to handle the jobs from the same customer order in the one factory to reduce production costs and enhance customer satisfaction [28, 29]. Obviously, if jobs belonging to the same customer order are processed in different factories, it is necessary to transport partially completed jobs among the factories in order to package the entire order together [30]. As far as we know, there is no literature on energy-saving RHFSP considering both distributed features and customer order constraints simultaneously. In this paper, we tackle the EDORHFSP under TOU electricity price with the objectives of minimizing the maximum completion time and total energy consumption cost.

The remainder of the paper is organized as follows. In Section 2, the EDORHFSP is detailed. The MA is designed in detail to solve the EDORHFSP in Section 3. In Section 4, numerical experiments are given. Finally, some conclusions and future works.

2 Problem description and formulation

2.1 Problem description

EDORHFSP can be described as follows: r customer orders need to be processed in f factories. Each customer order l contains λ_l jobs, and each customer order can be processed by any one of the f factories. Each factory contains s stages. Each stage has one or more parallel machines and some jobs may need to visit the same stage several times. Besides the classic flow shop constraints, there are some additional constraints according to the distributed characteristics and customer order constraints. The jobs belong to the same customer order must be processed in one factory. Jobs from different orders can be intermingled in the factory. Once a job is assigned to a factory, all its operations must be processed in the same factory. Thus, the major task of EDORHFSP is to assign customer orders to the factories, assign jobs to machines, determine the set of operations on all machines, so as to minimize maximum completion time and total energy cost. Figure 1 shows the diagram of EDORHFSP.

Fig.1

Diagram of the EDORHFSP.

2.2 Symbols definition

f: number of factories

s: number of stages in each factory

r: number of customer orders

h: number of re-entrances

g: index for factories, g = 1, 2, . . . , f

l, w: index for customer orders, l, w = 1, 2, . . . , r

λ_l: number of jobs in customer order l

J_l: jobs set belonging to customer order l

n: number of jobs, namely $n = \sum_{l = 1}^{r} λ_{l}$

k: index of operations for the job

j,v: index for jobs, j, v = 1, 2, . . . , n

C_j: completion time of job j

N_j: operations set of job j

s: number of stages in each factory

i: index for stages, i = 1, 2, . . . , s

m_ig: number of machines at stage i in factory g

a: index for machines at stage i in factory g, a = 1, 2, . . . , m_ig

M_total: number of machines at all stages in all factories g, $M_{total} = \sum_{g = 1}^{f} \sum_{i = 1}^{s} m_{ig}$

q: index for machines

R_q: number of operations on machine q

O_jk: the k-th operation of job j

M: a large positive number

U = {e₀ + 1, e₀ + 2, . . . , e₁, e₁ + 1, e₁ + 2, . . . , e₂, . . . , e_r-1 + 1, e_r-1 + 2, . . . , e_r}, in which, {e_l-1 + 1, e_l-1 + 2, . . . , e_l} are the jobs belonging to customer order l. $e_{l} = \sum_{w = 0}^{l} λ_{w}$ , l = 1, 2, . . . , r, and we can get e_r = n. Note that if the original jobs order does not fit this setting, it must be renumbered.

U_ig: operations set that are processed at stage i in factory g

S_jkg: starting time of O_jk in factory g

E_jkg: ending time of O_jk in factory g

p_jkiag: processing time of O_jk in factory g at stage i on machine a

T_q: completion time of machine q

S_q: start-up time of machine q

PI_q: unit energy consumption of machine q in the idle state

PW_q: unit energy consumption of machine q in the processing state

f (t): TOU electricity price function

r_ijkag: binary variable, 1 if O_jk is processed at stage i on machine a in factory g, and 0 otherwise

Y_jg: binary variable, 1 if the job j is processed in factory g, and 0 otherwise

Z_jkj′k′: binary variable, 1 if O_jk precedes O_j′k′, and 0 otherwise

$y_{q}^{t}$ : binary variable, 1 if the machine q is in processing state at time t, and 0 otherwise

2.3 Objectives functions and main constraints

Two objective functions C_max and TEC are listed as follows based on the existing literatures [7, 31].

Objective functions: $minimize C_{\max} = \max (C_{j})$ (1) $minimizeTEC = PC + IC$ (2)

In which, C_max and TEC represent the maximum completion time and the total energy consumption cost respectively. In addition, TEC is made up of two parts, one is the energy consumption cost PC when machines are in the processing state, the other one is the energy consumption cost IC when machines are in the idle state. $PC = \sum_{q = 1}^{M_{total}} \sum_{t = S_{q}}^{T_{q}} {PW}_{q} y_{q}^{t} f (t)$ (3) $IC = \sum_{q = 1}^{M_{total}} \sum_{t = S_{q}}^{T_{q}} {PI}_{q} (1 - y_{q}^{t}) f (t)$ (4)

s.t. $\sum_{g = 1}^{f} Y_{jg} = 1, \forall j$ (5)

$\begin{matrix} Y_{jg} = & Y_{vg}, j, v = e_{l - 1} + 1, e_{l - 1} + 2, . . ., e_{l}, \\ j \neq v, l = 1, 2, . . ., r \forall g \end{matrix}$ (6) $\sum_{a = 1}^{m_{ig}} r_{ijkag} = 1, \forall i, j, k, g, O_{jk} \in U_{ig}$ (7)

$\begin{matrix} \sum_{a = 1}^{m_{i^{' g}}} r_{i^{'} jk + 1 ag} S_{jk + 1 g} \\ ⩾ \sum_{a = 1}^{m_{ig}} r_{ijkag} (S_{jkg} + p_{jkiag}), \forall i, i^{'}, j, k, g \end{matrix}$ (8)

$\begin{matrix} M (1 - Z_{{jkj}^{'} k^{'}}) + M (1 - r_{ijkag}) + M (1 - r_{{ij}^{'} k^{'} ag}) \\ + (S_{j^{'} k^{'} g} - S_{jkg}) ⩾ p_{jkiag} \forall i, j < j^{'}, a, O_{jk} \\ \in U_{i}, O_{j^{'} k^{'}} \in U_{ig} \end{matrix}$ (9)

$\begin{matrix} M (1 - r_{ijkag}) + M (1 - r_{{ij}^{'} k^{'} ag}) + {MZ}_{{jkj}^{'} k^{'}} \\ + (S_{jkg} - S_{j^{'} k^{'} g}) ⩾ p_{j^{'} k^{'} iag} \forall i, j < j^{'}, a, O_{jk} \\ \in U_{i}, O_{j^{'} k^{'}} \in U_{i g} \end{matrix}$ (10)

$\begin{matrix} M (1 - r_{ijkag}) + M (1 - r_{{ijk}^{'} ag}) + (S_{{jk}^{'} g} - S_{jkg}) \\ ⩾ p_{jkiag} \forall i, a, j, k < k^{'}, O_{jk} \in U_{i}, O_{{jk}^{'}} \in U_{ig} \end{matrix}$ (11) $C_{j} = \sum_{i = 1}^{s} \sum_{a = 1}^{m_{ig}} r_{{ijN}_{j} ag} (S_{{jN}_{j} g} + P_{{jN}_{j} iag}), \forall j, g$ (12) $C_{j} ⩽ C_{\max}, \forall j$ (13) $C_{j} ⩾ 0, \forall j$ (14) $S_{jkg} ⩾ 0, \forall j, k, g$ (15) $E_{jkg} = S_{jkg} + p_{jkiag}, \forall j, k, g$ (16)

Formulas (1) and (2) are two objective functions. Constraint (5) means that each job must be assigned to one and only one factory. Constraint (6) ensures that jobs belonging to the same customer order should be processed in the same factory. Constraint (7) ensures that each operation can only be processed on one machine at the corresponding stage. Constraint (8) ensures that the starting time is not earlier than the ending time of the previous operation. Constraints (9), (10) and (11) mean that one machine processes at most one operation at a time. Constraints (12), (13) and (14) define the C_j. Constraints (15) and (16) define the relationship between S_jkg andE_jkg.

In the bi-objective optimization problem, solution X₁ dominates X₂ (denotes as X₁ ≺ X₂), if and only if ∀i ∈ {1, 2} , f_i (X₁) ⩽ f_i (X₂) ∧ ∃ i ∈ {1, 2}, f_i (X₁) < f_i (X₂). A solution X is the non-dominant solution when there is no solution that can dominate it. EDORHFSP usually involves a large number of jobs, machines and operations, and due to the re-entrant characteristics, it is difficult to solve using the mathematical model directly. Given this situation, a memetic algorithm is proposed to solve the problem with objectives of minimizing the maximum completion time and total energy consumption cost simultaneously.

3 The Memetic Algorithm for the EDORHFSP

Memetic algorithm (MA) is an algorithm framework based on the combination of global search and local search. Different search strategies can form different cultural MA. For example, genetic algorithm can be used for global search strategy, tabu search and simulated annealing can be used for local search strategy. By establishing the balance between global search and local search, MA is widely used in various multi-objective optimization problems and shows good performance. For the EDORHFSP, an Memetic algorithm is designed. Moreover, to enhance the performance of MA, encoding and decoding methods, energy cost saving procedure, three heuristic rules about the population initialization and some neighborhood search methods are designed. Figure 2 shows the flow chart of MA for the EDORHFSP.

Fig.2

The flow chart of MA for the EDORHFSP.

3.1 Encoding and decoding methods

In the paper, a solution contains two segments: factories assignment for all customer orders and jobs sequence in each factory. To be specific, each solution consists of f vectors and each vector represents a sequence of jobs in the corresponding factory, for instance, n = 6, f = 2, r = 3, s = 3, J₁ = [1, 2, 3], J₂ = [4, 5] and J₃ = [6]. In the example, 3 customer orders are processed in 2 identical factories and there are 3, 2 and 2 machines at each stage, respectively. Specifically, customer order 1 includes job 1, 2 and 3, customer order 2 includes job 4 and 5 and customer order 3 just includes of job 6. Table 1 lists the processing times of all operations. For convenience, let us suppose that PW = 8 and PI = 1 in thepaper.

Table 1
Processing times

Processing times First-pass Re-entrant1

Stage 1 Stage 2 Stage 3 Stage 1 Stage 2 Stage 3

Job1 2 2 1 1 2 1

Job2 1 3 2 3 2 1

Job3 3 1 2 2 1 2

Job4 2 1 1 2 3 2

Job5 1 2 3 1 1 2

Job6 2 2 3 2 3 2

Processing times	First-pass	Re-entrant1
Job1	2	2	1	1	2	1
Job2	1	3	2	3	2	1
Job3	3	1	2	2	1	2
Job4	2	1	1	2	3	2
Job5	1	2	3	1	1	2
Job6	2	2	3	2	3	2

In different places or countries, TOU electricity price schemes vary greatly. Ttaking the electricity price scheme of Shanghai in 2019 as an example, which is shown in Table 2. Such pricing structure can encourage consumers to transfer the electricity consumption from peak periods to mid-peak or off-peak periods, because the price of electricity in peak periods is higher.

Table 2

TOU electricity price scheme

Type	Peak	Mid-peak	Off-peak
Periods	8 : 00–11 : 00	6 : 00–8 : 00	0 : 00–6 : 00
	13 : 00–15 : 00	11 : 00–13 : 00	22 : 00–24 : 00
		15 : 00–22 : 00
Price (CNY/Kwh)	1.202	0.749	0.285

The decoding method is detailed as follows:

Step 1: Select π^g form the solution π, here π^g refers to jobs sequence processed in factory g, initialize g = 1.

Step 2: Select the first job j from π^g, arrange all the operations of job j on the machines that can finish it as early as possible.

Step 3: Select the next job v from π^g and arrange operation O_vk to the appropriate machine from the available machine set (MS) for it.

Step 4: Select machine a from the MS, and traverse all idle periods [Starts, Ends] of it. The earliest starting time S_vkg of operation O_vk is shown in formula (17).

$S_{vkg} = \max {E_{vk - 1 g}, Starts}$ (17)

Step 5: According to Equation (18), find an appropriate insertion point for O_vk. If more than one qualified S_vkg is found, select the smallest one as the starting time. In addition, if the insertion condition is not met, then setting S_vkg = Lm_a, where Lm_a signifies the ending time of machine a at the current state.

$S_{vkg} + p_{vkiag} ⩽ Ends$ (18)

Step 6: Repeat step 4 and step 5 until all machines in MS are traversed. If more than one machine with the same S_vkg, we choose the machine with less idle time between insertion point and the previousoperation.

Step 7: Calculate E_vkgof operation O_vk, according to formula (19).

$E_{vkg} = S_{vkg} + p_{vkiag}$ (19)

Step 8: Traverse all jobs in the π^g and repeat step 2-step 7.

Step 9: Update g = g + 1, repeat step2–step8 until all factories have been decoded.

Step 10: According to the scheduling results of all operations, two objective function values C_max and TEC are calculated.

With above mentioned encoding method, a solution can be denoted as π = {π¹, π², . . . , π^f}. To make it easier to understand, let us suppose that π = {, [4, 5, 6]} where π¹ = [2, 1, 3] and π² = [4, 5, 6]. Take factory 1 for example, firstly, all of the operations of job 2 are scheduled on the machine that can finish them as early as possible. Then, each operation of job 1 will be arranged on the machine by finding the appropriate insertion points. By analogy, all operations of job 3 are arranged on the appropriate machines. Figures 3 and 4 are the gantt charts, in which operation “302” denotes the second operation of job 3. Two objective function values of thescheduling scheme are C_max= 15 hours and TEC = 360.0280 CNY.

Fig.3

Schedule of factory 1.

Fig.4

Schedule of factory 2.

3.2 Energy cost saving procedure

Because the waiting state of machines and jobs inevitably appears in the RHFSP, especially in the case of large-scale jobs, we can make full use of these waiting times to adjust the processing periods of the operations. Based on machine allocation and operation sequencing, the energy cost saving procedure is added to select the appropriate starting processing time for all operations, and try to avoid the operations being processed at high electricity price, so as to achieve the purpose of reducing the TEC without deteriorating the C_max. According to the decoding rules, all operations can only move to the right. Moreover, the shift of the next operation will affect the previous one, so the energy cost saving procedure should follow the rule from back to front. Specifically, energy cost saving procedure mainly includes two parts: right-shift (RS) procedure and local adjustment (LA) procedure. The detailed RS procedure is as follows, and the diagram is shown inFig. 5.

Fig.5

RS procedure.

Step 1: Select one factory and all operations are arranged according to the completion time.

Step 2: Under the premise of not affecting the maximum completion time of the selected factory, each operation should move to the right to avoid operations being processed at high electricity price as far as possible according the TOU electricity price.

Step 3: Repeat step 2 until all operations have been traversed.

Step 4: Repeat step 1 to step 3 until all factories have been traversed.

Step 5: Update the pareto archive with the newly explored solutions.

Taking Table 1 as an example, after adding the RS procedure to Figs. 3 and 4, the new gantt charts are shown in Figs. 6 and 7. The comparison shows that operations “106”, “301”, “406”, “506” and “601” have moved to the right. Investigate the reason we can from following aspects thinking: (1) The electricity price of operations “106” and “601” remains the same before and after the procedure. Although the RS procedure cannot reduce the energy consumption, it can strive for a certain right-shift space for the previous operations. (2) The right-shift of operation “406” can reduce the energy consumption cost by moving the processing period from high electricity price 1.202 yuan/Kwh to low electricity price 0.749 yuan/Kwh. Meanwhile, the idle time of machine 6 did not increase in factory 2. (3) Although the shift of operation “506” from high electricity price 1.202 yuan/Kwh to low electricity price 0.749 yuan/Kwh also resulted in an increase of the idle time, taking the two factors into consideration, the total energy cost is reduced. (4) Although the electricity price of operation “301” before and after the RS procedure remains unchanged, energy cost can be reduced by shortening the idle time in the factory. To sum up, two objective function values of the scheduling scheme are C_max= 15 hours and TEC = 349.7712 CNY after the RS procedure, and TEC is reduced by 10.2568 CNY by comparison.

Fig.6

Schedule of factory 1 after RS procedure.

Fig.7

Schedule of factory 2 after RS procedure.

After the RS procedure, some operations may be transferred to other machines at the same stage, which we call it LA procedure. Without deteriorating the maximum completion time, some machines can be shut-down in advance or power-on later by the LA procedure, which is help to reduce the idle energy consumption cost (IC). Figures 8 and 9 are the gantt charts after the LA procedure. The processing period of operation “301” is at 1–4 hours in Fig. 8. Machine 1, which is at the same stage as machine 3, is idle at 1–4 hours. Therefore, operation “301” can be transferred to machine 1 for processing, and machine 3 can be powered on later, thus reducing idle energy consumption effectively. In summary, after the LA procedure, two objective function values of the scheduling scheme are C_max= 15 hours and TEC = 349.2089 CNY, and TEC is reduced by 0.5623 CNY by comparison.

Fig.8

Schedule of factory 1 after LA procedure.

Fig.9

Schedule of factory 2 after LA procedure.

The pseudo code for LA procedure is as follows, in which, Idle(i) denotes the total idle time on all machines at stage i.

The pseudo code for LA procedure:
1:	for g = 1:f
2:	for i = 1:s
3:	Calculate the Idle(i) before adjustment for stage i
4:	for each O operation inU_i
5:	for a = 1:m_ig
6:	if operation O can be arranged on machine a
7:	Then calculate the Idle(i)’
8:	if Idle(i)>Idle(i)’ then
9:	Change machine for operation O
10:	Idle(i)=Idle(i)’
11:	end if
12:	end if
13:	end for
14:	Find the suitable schedule scheme for operation O
with minimum Idle(i)
15:	end for
16:	end for
17:	end for
18:	Update scheduling scheme and calculate two objective
function values

The energy consumption curves (EC curves) of all machines along with time in each factory are plotted in Figs. 10 and 11. Obviously, the electricity load shifts in some periods. If there is a difference between the dotted line and the solid line, it means that the electricity load has been transferred. Specifically, the electricity is charged at high (low) price about 8–9 (3–4) hours in factory 1. After the energy cost saving procedure, total energy consumption is decline in the former period, on the contrary, total energy consumption increased in the latter period. It’s the same situation in factory 2. Therefore, the TEC can be saved by transferring the usage from peak period to off-peak or mi-peak period after the energy cost saving procedure.

Fig.10

Energy consumption curves of factory 1.

Fig.11

Energy consumption curves of factory 2.

Take operation O_jk as an example, supposing operation O_jk is processed on machine a in factory g and S_j′k′g is the starting time of the next operation after O_jk on machine a. Meanwhile, the range of starting time of O_jk is [t_min, t_max]. The pseudo code for the energy cost saving procedure is as follows.

The pseudo code for the energy cost saving procedure:
1:	Get the starting time range of operation O_jk,
t_min = S_jk, t_max = min {S_j′k′g, S_jk+1g} - P_jkiag
2:	Initialize goodresult = TEC, when Start=t_min
3:	Divide the possible starting time, the number is (t_max - t_min)
4:	for i = 1: (t_max - t_min)
5:	Calculate the TEC, when Start = t_min + i
6:	if TEC≤good result
7:	goodresult = TEC
8:	goodstart=Start
9:	end if
10:	for
11:	The appropriate starting time of O _jk after RS procedure is goodstart
12:	LA procedure
13:	Update non-dominated solutions set

3.3 Population initialization

For the EDORHFSP, three heuristic rules are designed for assigning the customer orders to the appropriate factories and sequencing jobs in each factory. Rule1, Rule2 and Rule3 begin with constructing a job sequence θ randomly, and each job appears once in it.

3.3.1 Rule1

The specific steps of Rule1 are as follows:

Step 1: Generate a jobs sequence θ randomly.

Step 2: Remove the first job from θ and insert it into the first stage of the first factory.

Step 3: If θ is not empty, remove the next job j of customer order l. If no jobs of the customer order l have been extracted from θ, then randomly select a factory g and place job j at the end of sequence π^g, where π^g is the jobs sequence to be processed in factory g. if one or more jobs of customer order l have already been assigned to the factory (for example factory G), then place job j at the end of sequence π^G.

Step4: Repeat Step3 until θ is empty.

Step5: Get π = {π¹, π², . . . , π^f}.

3.3.2 Rule2

In Rule2, consider C_max when arranging jobs. The details of Rule2 is given as follows.

Step 1: Generate a jobs sequence θ randomly.

Step 2: Remove the first job from θ and insert it into the first stage of the first factory.

Step 3: If θ is not empty, remove the next job j of customer order l. If no jobs of the customer order l have been extracted from θ, then tentatively insert job j into all possible positions of π^g where g = 1, 2, . . . , f. Through calculation, we choose to put the job j in the position with minimum C_max of all factories, then the job sequence π^g^′ is obtained after arranging job j. if there are some jobs belonging to customer order l have been assigned to the factory (for example factory G), then job j must be assigned to the factory G and place job j in the position with minimum C_max, then the job sequence π^G is obtained after arranging job j.

Step4: Repeat Step3 until θ is empty.

Step5: Get π = {π¹, π², . . . , π^f}.

3.3.3 Rule3

Rule3 is similar to Rule2, the difference is that when assigning jobs, we take minimizing TEC as the criterion. Here, we will not elaborate on the specific steps of Rule3.

Figure 12 shows the hybrid initialization strategy. P_init is the percentage of initialization population to execute Rule2, and Rule3 strategies, and the remaining (1-2*P_init) will execute the Rule1 strategy.

Fig.12

Population initialization based on hybrid strategy.

3.4 Neighborhood search

According to the characteristics of EDORHFSP, five neighborhood search methods are designed in the paper:

Orders_Swap_among_Factories, Orders_Insert_among_Factories, Job_Best_Insert_within_Factory, Job_Best_Swap_within_Factory and

Jobs_Inverse_within_Factory.

3.4.1 Orders_Swap_among_Factories

Firstly, find the critical factory f_c with maximum completion time among all factories. Then, try to exchange the customer order in factory f_c with other customer orders in the remaining factories sequentially. Next, the jobs of two chosen orders are reinserted into the appropriate position of the target factories one by one. Accept the exchange if the C_max can be reduced until all possible exchanges are completed.

3.4.2 Orders_Random_Insert_other_Factory

At first, find the critical factory f_c with maximum completion time among all factories. Then, select a customer order randomly form π^fc and insert it into another factory (select a factory randomly different from f_c). Obviously, the procedure does not guarantee to find a better solution, but it can greatly increase the diversity of individuals.

3.4.3 Job_Best_Insert_within_Factory

A job j is randomly extracted from π^g to generate jobs sequence π^g^′ and inserted it into all possible positions of π^g^′. The job will be re-assigned, if a better solution is gotten. The pseudo code for the Job_Best_Insert_within_Factory is as follows.

The pseudo code for the Job_Best_Insert_within_Factory
(factory g, individual x)
1:	Get the sequence π^g^′ after extracting n-th job randomly for π^g, n∈ (1, \|π^g\|)
% \|π^g\| indicates the number of jobs processed in factory g
2:	% Try to insert job into π^g^′ with the number of \|π^g^′\|+1 possible positions
3	for i = 1: \|π^g^′\|+1
4	x′ = x, temp = zeros(1, \|π^g\|), temp (i) = n
5	if i==1
6	temp(2:end)=π^g^′
7	elseif i==\|π^g^′\|+1
8	temp(1:i-1)=π^g^′
9	Else
10	temp(1:i-1)=π^g^′(1:i-1), temp(i + 1:end)=π^g^′(i:end)
11	end if
12:	Update the π^g in individual x′ using temp
13:	% $\bar{x'}$ indicates the objective function values [C_max, TEC] of the individual x′.
14:	if $\bar{x^{'}} ≺ \bar{x}$
15:	x = x′
16:	end if
17:	end for
18:	Get the individual x

3.4.4 Job_Best_Swap_within_Factory

This procedure attempts to swap the job j with any one of the rest jobs within the factory. If a better solution is obtained, the jobs sequence is updated. The pseudo code for the Job_Best_Swap_within_Factory is as follows.

The pseudo code for the Job_Best_Swap_within_Factory
(factory g, individual x)
1:	Get the sequence π^g^′ after extracting n-th job randomly for π^g, n∈ (1, \|π^g\|)
% \|π^g\| indicates the number of jobs processed in factory g
2:	% Swap the n-th job with the rest jobs in π^g^′ one by one
3:	for i = 1: \|π^g\|
4:	x′ = x, temp = π^g
5:	tt=temp(i), temp(i)=temp(n), temp(n)=tt
6:	Update the π^g in individual x’ using temp
7:	% $\bar{x'}$ indicates the objective function values [C_max, TEC] of the individual x′.
8:	if $\bar{x^{'}} ≺ \bar{x}$
9:	x = x′
10:	end if
11:	end for
12:	Get the individual x

3.4.5 Jobs_Inverse_within_Factory

Invert the sub-sequence randomly between two different positions within the critical factory f_c.

Due to a large number of jobs, machines and operations are involved in this EDORHFSP, the neighborhood local search is conducted mainly within the critical factory f_c. In addition, to improve the diversity of the population, we will do some exploration in the other factories. The pseudo code for Neighborhood_Search_methods is as follows.

The pseudo code for Neighborhood_Search_methods (individual x):
1:	Find the critical factory f_c among all factories with the maximum completion time
2:	Job_Best_Insert_within_Factory(factoryf_c, individual x)
3:	Job_Best_Swap_within_Factory(factoryf_c, individual x)
4:	Orders_Swap_among_Factories(factoryf_c, individual x)
5:	if the solution cannot be improved
6:	Select a factory f_other randomly different from f_c
7:	Switch N % N=randperm(3,1)
8:	Case 1: Job_Best_Insert_within_Factory(factoryf_other, individual x)
9:	Case 2: Job_Best_Swap_within_Factory(factoryf_other, individual x)
10:	Case 3: Jobs_Inverse_within_Factory(factoryf_c, individual x)
11:	if rand < P_Insert
12:	Orders_Random_Insert_other_Factory(factoryf_c, factoryf_other)
13:	end if
14:	end if
15:	Update the individual x

4 Numerical experiments and comparisons

4.1 Generating test data

As a new problem, there is no ready-made test data in related research. Therefore, we will generate the test data based on the famous benchmark sets created by Cho et al. [6] for the RHFSP. For each instance, the number of jobs, stages, machines, re-entrances, processing times, customer orders and factories are random integers from the uniform distributions U[10, 50], U[5, 25], U[1, 6], U[1,6, 1,6], U[1, 30], U[5, 8] and U[2, 4], respectively.

In this paper, the TOU electricity price scheme consists of seven pricing periods one day. Figure 13 shows the piecewise function diagram of electricity price. In the figure, the price (1.202 CNY/Kwh) in peak periods is more than four times as high as that (0.285 CNY/Kwh) in off-peak periods. Such a huge change in electricity prices provides a huge opportunity for energy intensive manufacturing companies to save costs.

Fig.13

TOU electricity price over 24 hours.

All programs are coded in Matlab on the computer with i7-10510U CPU@1.80 GHz and 16 G RAM.

4.2 Performance metrics

The following three metrics (IGD, C, SP) are used to evaluate the performance of MA. In addition, since the optimal pareto fronts (OPF) of the instances are unknown, we use an approximate pareto front (APF) instead of OPF. Specifically, for each instance, APF can be obtained by combining the non-dominated solutions obtained by all the algorithms into one set and then removing the dominated solutions from the set. In addition, due to different dimensions, it is necessary to standardize before calculating them.

IGD (Inverted Generational Distance) is a comprehensive metric which is used to evaluate evaluates the convergence and distribution performance of the algorithm. The calculation formula is as shown in (20). $IGD = \frac{\sum_{x \in N^{*}} dist (x, N)}{| N^{*} |}$ (20)

In which, |N*| denotes the number of non-dominated solutions on the OPF. N is the non-dominated solutions set obtained by one algorithm. dist (x, N) represents the Euclidean distance between the solution x ∈ N^* and the nearest solution on N. Clearly, the smaller the IGD, the better the performance. In an extreme case, IGD = 0 means that N is a subset of N*.

C metric is used to measure the dominating relationship between two non-dominated solution sets. C(A, B) means the percentage of the solutions in B that are dominated by at least one solution in A. C (A, B) =0 indicates that no solution in B is dominated by any solution in A. C (A, B) =1 indicates that all solutions in B are dominated by some solutions in A. $C (A, B) = \frac{| {u \in B | \exists v \in A : v ≺ u} |}{| B |}$ (21)

SP is used to measure the distribution uniformity of the non-dominated solutions on the pareto front. In an extreme case, SP = 0 indicates that the non-dominated solutions on the pareto front are evenly distributed. $SP = \sqrt{\frac{1}{| N | - 1} \sum_{i = 1}^{| N |} (\bar{D} - D_{i})^{2}}$ (22)

Where, $D_{i} = min (\sum_{k = 1}^{M} | f_{k}^{i} - f_{k}^{j} |)$ , i ≠ j, $\bar{D} = \frac{1}{n} \sum_{i = 1}^{n} D_{i}$ and |N| is the number of non-dominated solutions.

4.3 Parameters setting

MA mainly involves four key parameters: PopSize (population size), P_init (the percentage of initialization population to execute Rule2 and Rule3 in population initialization), P_Insert (the probability of performing Orders_Random_Insert_other_Factory), G_max (the maximum number of generations). Taguchi method of design-of-experiments (DOE) is used to investigate the effect of parameters setting on the performance of MA. As shown in Table 3, each parameter contains four factor levels. Moreover, the orthogonal array L₁₆(4⁴) is used and the MA will run 20 times independently for every parameter combination on each instance. After all the experiments are done, the average of dominant metric Ω of each parameter combination is regarded as the response variable (RV) value, as shown in Table 4. The dominant metric Ω represents the percentage of non-dominated solutions obtained by an algorithm that are also non-dominated solutions on the OPF.

Table 3
Parameters combination

Parameters Factor level

1 2 3 4

PopSize 20 40 60 80

P _init 0.1 0.2 0.3 0.4

P _Insert 0.3 0.4 0.5 0.6

G _max 60 80 100 120

Parameters	Factor level
PopSize	20	40	60	80
P _init	0.1	0.2	0.3	0.4
P _Insert	0.3	0.4	0.5	0.6
G _max	60	80	100	120

Table 4

Orthogonal array and ARV values

Experiment number	Parameter				ARV
	PopSize	P _init	P _Insert	G _max
1	20	0.1	0.10	60	0.0642
2	20	0.2	0.15	80	0.0305
3	20	0.3	0.20	100	0.0364
4	20	0.4	0.25	120	0.0546
5	40	0.1	0.15	100	0.0503
6	40	0.2	0.10	120	0.0807
7	40	0.3	0.25	60	0.0877
8	40	0.4	0.20	80	0.0612
9	60	0.1	0.20	120	0.0858
10	60	0.2	0.25	100	0.1041
11	60	0.3	0.10	80	0.1113
12	60	0.4	0.15	60	0.0325
13	80	0.1	0.25	80	0.0758
14	80	0.2	0.20	60	0.0347
15	80	0.3	0.15	120	0.0617
16	80	0.4	0.10	100	0.1050

The average RV and the significance rank are shown in Table 5. The larger the delta is, the more significant the parameter is. Table 5 shows that P_Insert is the most significant parameter, PopSize ranks the second, G_max and P_init ranks the third and the last, respectively. To be specific, a proper P_Insert is help to prevent the algorithm from falling into local optimum. For PopSize, a small PopSize will lead to insufficient sampling of the solution space, while a large PopSize could lead to inadequate evolutionary processes. As for P_init, a suitable value of P_init is helpful to balance the diversity and quality of the initial solutions. The trend of each factor is shown in Fig. 14. Based on the above analysis, it is suggested to select the appropriate parameter combination as PopSize = 60, P_init = 0.3, P_Insert = 0.1 and G_max = 100 in the following numerical experiments.

Table 5

RV and rank of each parameter

Level	PopSize	P _init	P _Insert	G _max
1	0.04643	0.06903	0.09030	0.05477
2	0.06997	0.06250	0.04375	0.06970
3	0.08342	0.07427	0.05453	0.07395
4	0.06930	0.06332	0.08055	0.07070
Delta	0.03700	0.01177	0.04655	0.01918
Rank	2	4	1	3

Fig.14

Factor level trend.

4.4 Effectiveness of population initialization

In this paper, three heuristic rules are designed for population initialization. It is necessary to compare MA with MA1 (random initialization using Rule1 only) to verify the effectiveness of population initialization. Each instance runs for 20 times independently. The average C metric of 8 instances and the p values obtained by Wilcoxon signed-rank test are listed in Table 6 and the boxplots of C metric are shown in Fig. 15. When p value is less than 0.05, there is a significant difference between the two algorithms.

Table 6
Comparison of C metric between MA and MA1

Scale C(MA,MA1) C(MA1,MA) p C(MA,MA1) C(MA1,MA) p C(MA,MA1) C(MA1,MA) p

(n,s,h) r = 4f = 2 r = 6f = 3 r = 8f = 4

(20,10,1) 0.6078 0.1713 0.001 0.6625 0.1518 0 0.5975 0.19852 0.002

(14,10,1) 0.8394 0.1369 0 0.6690 0.1879 0 0.6043 0.2280 0.003

(19,9,1) 0.7336 0.1934 0 0.7750 0.1809 0 0.6833 0.2036 0.001

(16,6,1) 0.6933 0.1685 0 0.5796 0.1645 0.002 0.7754 0.2198 0.001

(45,13,5) 0.6104 0.3735 0.022 0.7663 0.1936 0 0.7052 0.1638 0

(10,21,5) 0.7599 0.1566 0 0.7485 0.2090 0.001 0.5877 0.2514 0.008

(29,7,5) 0.6482 0.2061 0 0.6760 0.1625 0 0.5011 0.2236 0.011

(29,14,5) 0.6275 0.1998 0.002 0.6749 0.1663 0.001 0.4866 0.2508 0.033

Scale	C(MA,MA1)	C(MA1,MA)	p	C(MA,MA1)	C(MA1,MA)	p	C(MA,MA1)	C(MA1,MA)	p
(20,10,1)	0.6078	0.1713	0.001	0.6625	0.1518	0	0.5975	0.19852	0.002
(14,10,1)	0.8394	0.1369	0	0.6690	0.1879	0	0.6043	0.2280	0.003
(19,9,1)	0.7336	0.1934	0	0.7750	0.1809	0	0.6833	0.2036	0.001
(16,6,1)	0.6933	0.1685	0	0.5796	0.1645	0.002	0.7754	0.2198	0.001
(45,13,5)	0.6104	0.3735	0.022	0.7663	0.1936	0	0.7052	0.1638	0
(10,21,5)	0.7599	0.1566	0	0.7485	0.2090	0.001	0.5877	0.2514	0.008
(29,7,5)	0.6482	0.2061	0	0.6760	0.1625	0	0.5011	0.2236	0.011
(29,14,5)	0.6275	0.1998	0.002	0.6749	0.1663	0.001	0.4866	0.2508	0.033

Fig.15

Boxplot of C metric between MA and MA1.

Table 6 shows that C(MA, MA1) is larger than C(MA1, MA) on all instances, which implies that most non-dominated solutions obtained by MA1 are dominated by those obtained by MA. In addition, the p values of all the instances are less than 0.05, which indicates that the designed population initialization procedure is significantly effective to solve the EDORHFSP.

4.5 Effectiveness of energy cost saving

To demonstrate the effectiveness of energy cost saving, we compare the MA with MA without energy cost saving, denoted as MA2. Each instance runs 20 times independently. The average C metric of 8 instances and the p values obtained by Wilcoxon signed-rank test are recorded in Table 7 and Fig. 16 is the boxplots of C metric.

Table 7
Comparison of C metric between MA and MA2

Scale C(MA,MA2) C(MA2,MA) p C(MA,MA2) C(MA2,MA) p C(MA,MA2) C(MA2,MA) p

(n,s,h) r=4 f = 2 r=6 f = 3 r=8 f = 4

(20,10,1) 0.7762 0.1123 0 0.7598 0.0745 0 0.8432 0.0639 0

(14,10,1) 0.8160 0.1052 0 0.8394 0.1004 0 0.8274 0.0973 0

(19,9,1) 0.8882 0.0503 0 0.8710 0.0462 0 0.8604 0.0583 0

(16,6,1) 0.6929 0.1112 0 0.7298 0.0926 0 0.7744 0.0762 0

(45,13,5) 0.4348 0.1992 0.024 0.6064 0.2313 0.013 0.7727 0.1710 0

(10,21,5) 0.5468 0.1321 0.003 0.6796 0.1709 0.001 0.7533 0.1969 0.001

(29,7,5) 0.6650 0.1418 0 0.6936 0.1040 0 0.7224 0.0852 0

(29,14,5) 0.5687 0.2049 0.011 0.7078 0.2329 0.001 0.6385 0.2185 0.005

Scale	C(MA,MA2)	C(MA2,MA)	p	C(MA,MA2)	C(MA2,MA)	p	C(MA,MA2)	C(MA2,MA)	p
(20,10,1)	0.7762	0.1123	0	0.7598	0.0745	0	0.8432	0.0639	0
(14,10,1)	0.8160	0.1052	0	0.8394	0.1004	0	0.8274	0.0973	0
(19,9,1)	0.8882	0.0503	0	0.8710	0.0462	0	0.8604	0.0583	0
(16,6,1)	0.6929	0.1112	0	0.7298	0.0926	0	0.7744	0.0762	0
(45,13,5)	0.4348	0.1992	0.024	0.6064	0.2313	0.013	0.7727	0.1710	0
(10,21,5)	0.5468	0.1321	0.003	0.6796	0.1709	0.001	0.7533	0.1969	0.001
(29,7,5)	0.6650	0.1418	0	0.6936	0.1040	0	0.7224	0.0852	0
(29,14,5)	0.5687	0.2049	0.011	0.7078	0.2329	0.001	0.6385	0.2185	0.005

Fig.16

Boxplot of C metric between MA and MA2.

Table 7 shows that C(MA, MA2) is larger than C(MA2, MA) on all instances, which implies that most non-dominated solutions obtained by MA2 are dominated by those obtained by MA. In addition, p values of all the instances are less than 0.05, which means that the energy cost saving procedure is significantly effective to solve the EDORHFSP. The pareto fronts on the scale (20,10,1) instance with r = 8, f = 4 and the scale (16,6,1) instance with r = 4, f = 2 are shown in Figs. 17 and 18, which further verify the above conclusion.

Fig.17

Pareto front on the scale (20,10,1) instance with r = 8 and f = 4.

Fig.18

Pareto front on the scale (16,6,1) instance with r = 4 and f = 2.

4.6 Comparison with other algorithms

To validate the performance of the proposed MA, SPEAII, MODE and NSGAII are selected for comparative study. 8 instances of different sizes are randomly selected from the 240 instances with three different [r, f] combinations. Each instance runs 20 times independently. To make an overall comparison, the average values of three performance metrics are listed in the Tables 8, 9 and 10. The boxplots of three metrics are shown in Figs. 19 and 20. As we all know, the average values of three metrics can only represent the performance from the macroscopic perspective, so it is necessary to carry out the Wilcoxon signed-rank test to show the difference is significant or not.

Table 8
Performance metrics of three algorithms (r = 4, f = 2, A: MA, B: MODE, C: SPEA-II, D: NSGA-II)

Scale IGD/SP C metric

(n,s,h) MA MODE SPEA-II NSGA-II C(A,B) C(B,A) p C(A,C) C(C,A) p C(A,D) C(D,A) p

(20,10,1) 0.0747/0.1492 0.2821/0.2378 0.4053/0.5194 0.2705/0.382 0.5472 0.2059 0.001 0.6476 0.1083 0 0.4995 0.3452 0.041

(14,10,1) 0.0238/0.0617 0.2783/0.2377 0.5137/0.4823 0.2275/0.3654 0.6777 0.2668 0.001 0.7756 0.2094 0 0.6775 0.3201 0.001

(19,9,1) 0.0661/0.0092 0.1921/0.1909 0.3872/0.6623 0.2647/0.3341 0.5885 0.1821 0 0.6107 0.0764 0 0.4796 0.3845 0.036

(16,6,1) 0.0737/0.1502 0.1897/0.2923 0.5451/0.4447 0.1859/0.3058 0.5157 0.0799 0 0.7745 0.1675 0 0.6547 0.3075 0.002

(45,13,5) 0.0111/0.1637 0.3385/0.4395 0.5931/0.4455 0.168/0.2967 0.4918 0.2011 0.032 0.5657 0.2198 0.003 0.5800 0.2021 0.004

(10,21,5) 0.0258/0.0484 0.3634/0.2562 0.3411/0.4660 0.1831/0.2933 0.6661 0.3246 0.021 0.6474 0.1826 0 0.699 0.3082 0

(29,7,5) 0.0617/0.1575 0.3254/0.3684 0.3217/0.3897 0.2294/0.3744 0.8042 0.1901 0 0.6173 0.2087 0 0.4603 0.3467 0.037

(29,14,5) 0.0933/0.0407 0.2932/0.2502 0.4094/0.4497 0.2711/0.3656 0.7727 0.2676 0 0.8167 0.1827 0 0.5551 0.3775 0.021

Average 0.0537/0.0975 0.2828/0.2841 0.4395/0.4825 0.2250/0.3396 0.6329 0.2147 0.6819 0.1694 0.5757 0.3240

Scale	IGD/SP	C metric
(20,10,1)	0.0747/0.1492	0.2821/0.2378	0.4053/0.5194	0.2705/0.382	0.5472	0.2059	0.001	0.6476	0.1083	0	0.4995	0.3452	0.041
(14,10,1)	0.0238/0.0617	0.2783/0.2377	0.5137/0.4823	0.2275/0.3654	0.6777	0.2668	0.001	0.7756	0.2094	0	0.6775	0.3201	0.001
(19,9,1)	0.0661/0.0092	0.1921/0.1909	0.3872/0.6623	0.2647/0.3341	0.5885	0.1821	0	0.6107	0.0764	0	0.4796	0.3845	0.036
(16,6,1)	0.0737/0.1502	0.1897/0.2923	0.5451/0.4447	0.1859/0.3058	0.5157	0.0799	0	0.7745	0.1675	0	0.6547	0.3075	0.002
(45,13,5)	0.0111/0.1637	0.3385/0.4395	0.5931/0.4455	0.168/0.2967	0.4918	0.2011	0.032	0.5657	0.2198	0.003	0.5800	0.2021	0.004
(10,21,5)	0.0258/0.0484	0.3634/0.2562	0.3411/0.4660	0.1831/0.2933	0.6661	0.3246	0.021	0.6474	0.1826	0	0.699	0.3082	0
(29,7,5)	0.0617/0.1575	0.3254/0.3684	0.3217/0.3897	0.2294/0.3744	0.8042	0.1901	0	0.6173	0.2087	0	0.4603	0.3467	0.037
(29,14,5)	0.0933/0.0407	0.2932/0.2502	0.4094/0.4497	0.2711/0.3656	0.7727	0.2676	0	0.8167	0.1827	0	0.5551	0.3775	0.021
Average	0.0537/0.0975	0.2828/0.2841	0.4395/0.4825	0.2250/0.3396	0.6329	0.2147		0.6819	0.1694		0.5757	0.3240

Table 9

Performance metrics of three algorithms (r = 6, f = 3, A: MA, B: MODE, C: SPEA-II, D: NSGA-II)

Scale	IGD/SP				C metric
(n,s,h)	MA	MODE	SPEA-II	NSGA-II	C(A,B)	C(B,A)	p	C(A,C)	C(C,A)	p	C(A,D)	C(D,A)	p
(20,10,1)	0.0107/0.1047	0.1467/0.2497	0.5449/0.2661	0.1794/0.2496	0.5962	0.3106	0.042	0.8472	0.0922	0	0.4677	0.3563	0.028
(14,10,1)	0.0567/0.0880	0.1841/0.1816	0.5882/0.6004	0.1984/0.3339	0.6683	0.2215	0.001	1	0	0	0.6901	0.2944	0
(19,9,1)	0.0118/0.1154	0.1588/0.2634	0.1936/0.2574	0.2765/0.2731	0.6413	0.2980	0.004	0.8268	0.0303	0	0.4412	0.3256	0.016
(16,6,1)	0.0623/0.0881	0.1490/0.1997	0.5184/0.3417	0.2345/0.3992	0.7401	0.2011	0	0.7924	0.1885	0	0.6415	0.3034	0.033
(45,13,5)	0.0787/0.1369	0.3056/0.1702	0.2664/0.5034	0.2985/0.3231	0.6864	0.2375	0	0.6273	0.2336	0.001	0.6525	0.3149	0.007
(10,21,5)	0.0421/0.0109	0.3024/0.2396	0.4833/0.6009	0.1646/0.3306	0.6987	0.2682	0.001	0.6938	0.1305	0	0.6659	0.3256	0.012
(29,7,5)	0.0234/0.1145	0.2970/0.1793	0.5119/0.6726	0.2706/0.2480	0.7772	0.2163	0	0.8013	0.1251	0	0.4620	0.2593	0
(29,14,5)	0.0628/0.1391	0.2563/0.2603	0.3740/0.2782	0.2449/0.4411	0.6335	0.1925	0	0.8293	0.1031	0	0.5439	0.2894	0
Average	0.0374/0.0996	0.2249/0.2179	0.4350/0.4400	0.2334/0.3248	0.6802	0.2432		0.8022	0.1129		0.5706	0.3086

Table 10

Performance metrics of three algorithms (r = 8, f = 4, A: MA, B: MODE, C: SPEA-II, D: NSGA-II)

Scale	IGD/SP				C metric
(n,s,h)	MA	MODE	SPEA-II	NSGA-II	C(A,B)	C(B,A)	p	C(A,C)	C(C,A)	p	C(A,D)	C(D,A)	p
(20,10,1)	0.0303/0.1149	0.2162/0.1901	0.3028/0.2745	0.1723/0.3499	0.7438	0.2218	0	0.8212	0.1401	0	0.5076	0.2515	0
(14,10,1)	0.0664/0.1677	0.2045/0.3073	0.3052/0.2057	0.2109/0.3453	0.7385	0.2008	0	1	0	0	0.6480	0.3095	0.018
(19,9,1)	0.0092/0.1303	0.1012/0.2593	0.3299/0.2660	0.1978/0.4402	0.7262	0.2787	0	0.7439	0.2129	0	0.5522	0.3733	0.021
(16,6,1)	0.0491/0.0844	0.1850/0.1773	0.2935/0.1974	0.2319/0.2610	0.7614	0.1725	0	1	0	0	0.6700	0.3300	0
(45,13,5)	0.0389/0.1767	0.2255/0.3647	0.3505/0.2662	0.2661/0.2654	0.6535	0.2884	0.001	0.7582	0.2166	0	0.4873	0.3442	0.032
(10,21,5)	0.0399/0.1057	0.2503/0.2679	0.3646/0.3109	0.1479/0.3634	0.8393	0.1045	0	0.5765	0.2573	0.004	0.5086	0.2510	0.001
(29,7,5)	0.0404/0.1226	0.2385/0.1937	0.4003/0.1762	0.2696/0.2880	0.7201	0.1691	0	0.7418	0.1567	0	0.4778	0.3502	0.037
(29,14,5)	0.0298/0.1458	0.1991/0.2671	0.2052/0.3576	0.2730/0.4273	0.6470	0.2787	0.001	0.6685	0.1142	0	0.4754	0.3011	0.032
Average	0.0380/0.1310	0.2025/0.2534	0.3190/0.2568	0.2211/0.3425	0.7287	0.2143		0.7887	0.1372		0.5409	0.3139

Fig.19

Boxplot of IGD and SP metrics for all algorithms.

Fig.20

Boxplot of C metric for all algorithms.

The above three tables show that IGD and SP of the MA is smaller than those of the MODE, SPEA-II and NSGA-II on every instance, which implies that the MA is superior to the MODE, SPEA-II and NSGA-II in terms of IGD and SP metrics. Meanwhile, C(MA, MODE), C(MA, SPEA-II) and C(MA, NSGA-II) are greater than C(MODE, MA), C(SPEA-II, MA) and C(NSGA-II, MA) respectively on all instances, which implies that most non-dominated solutions obtained by MODE, SPEA-II and NSGA-II are dominated by those obtained by MA. In terms of C metric, since most of p values are smaller than 0.05, therefore, MA is clearly superior to that of MODE, SPEA-II and NSGA-II. Figures 19 and 20 further verify the above conclusions.

Take the instance with scale (19,9,1) as an example, in which there are 19 jobs, 9 stages, 1 re-entrance, 2 factories, 5 customer orders and J₁ = {1, 2} , J₂ = {3} , J₃ = {4, 5, 6, 7, 8, 9} , J₄ = {10, 11} and J₅ = {12, 13, 14, 15, 16, 17, 18, 19}. In addition, there are 2, 1, 2, 2, 1, 2, 1, 2 and 1 machines at each stage, respectively. Take a non-dominated solution π = {, linebreak [19, 11, 12, 14, 18, 10, 16, 3, 17, 13, 15]} as an example. Obviously, Customer orders J₁ and J₃ are processed in factory 1 and customer orders J₂, J₄ and J₅ are processed in factory 2. The Gantt charts without considering TOU are illustrated in Figs. 21 and 22 with C_max= 100 hours and TEC = 5730.9650 CNY. The Gantt charts considering TOU are illustrated in Figs. 23 and 24 with C_max= 100 hours and TEC = 5492.1233 CNY. Comparisons show that the total TEC is reduced by 4.17% to 238.8417 CNY without affecting the C_max. Exploring its reason, it is not difficult to discover that processing periods of some operations have been shifted to right. For example, the processing periods of operations “102”, “201”, “202”, “6016”, “706”, “304”, “305”, “1301”, “1302” and “1607” move from higher price to lower price. Meanwhile, the total energy consumption curves in the two factories are shown in Fig. 25. The figure shows that electricity load shifts obviously in some periods, such as, about 25 hours, 80 hours and 90 hours of factories 1, about 12 hours, 40 hours and 70 hours of factories 2. To study the change trend of the energy consumption cost before and after considering TOU, taking the corresponding scheduling scheme in Figs. 21–24 as an example, the energy consumption cost histogram is shown in Fig. 26. It suggests that in factory 1 and factory 2, the energy consumption cost in idle state increases slightly after considering TOU, but the energy consumption cost is significantly reduced in processing state, which leads to the reduction of TEC. The reason is that the unit energy consumption of machines in the processing state is far greater than that in the idle state. Therefore, it is necessary to reduce the processing period at the high electricity price as far as possible, but this process may increase the energy consumption cost in the idle period. To sum up, the energy consumption cost can be saved by avoiding the periods at the high electricity price, meanwhile, the MA algorithm can solve the EDORHFSP effectively.

Fig.21

Gantt charts without considering TOU of factory 1.

Fig.22

Gantt charts without considering TOU of factory 2.

Fig.23

Gantt charts considering TOU of factory 1.

Fig.24

Gantt charts considering TOU of factory 2.

Fig.25

Energy consumption curves.

Fig.26

The comparison about the energy consumption cost before and after considering TOU.

5 Conclusions and future works

In this study, a memetic algorithm is proposed to solve the energy-efficient distributed re-entrant hybrid flow shop scheduling problem considering the customer order constraints under TOU electricity price. some customer orders require production in multiple factories and jobs belonging to the same customer order must be processed in the same factory. The main contributions of this paper are as follows: (1) Order constraint is introduced to scheduling problem, which is closer to the actual production. (2) A mathematical model is established for the EDORHFSP aiming at minimizing the C_max and TEC. (3) Based on the characteristics of the problem, the encoding scheme, decoding methods and three heuristic rules of population initialization are designed. (4) In order to avoid local optimum and premature convergence, energy cost saving procedure and some local neighborhood search methods are designed to improve the performance of the algorithm. (5) The superiority and effectiveness of the algorithm are verified by a large number of simulation experiments.

Numerical experiments show that the proposed algorithm has better performances than the MODE, SPEA-II and NSGA-II algorithms and TEC can be reduced by transferring the operations from peak periods to mid-peak or off-peak periods. In the future, EDORHFSP will be further studied, such as designing better algorithms, introducing dynamic speed of machines, uncertain scheduling, andso on.

Funding

The work is supported by National Natural Science Foundation of China (NO. 71840003).

References

Wang

M.Y.

, Sethi

S.P.

, van de Velde

S.L.

, Minimizing Makespan in a Class of Reentrant Shops[J], Operations Research 45(5) (1997), 702–712.

Choi

H.S.

, Kim

J.S.

, Lee

D.H.

, Real-time scheduling for reentrant hybrid flow shops: A decision tree based mechanism and its application to a TFT-LCD line[J], Expert Systems with Applications 38(4) (2011), 3514–3521.

Graves

S.C.

, Meal

H.C.

, Stefek

, Zeghmi

A.H.

, Scheduling of re-entrant flow shops[J], Journal of Operations Management 3(4) (1983), 197–207.

Kim

H.W.

, Lee

D.H.

, Heuristic algorithms for re-entrant hybrid flow shop scheduling with unrelated parallel machines[J], Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture 223(4) (2009), 433–442.

Dugardin

, Yalaoui

, Amodeo

, New multi-objective method to solve reentrant hybrid flow shop scheduling problem[J], European Journal of Operational Research 203(1) (2010), 22–31.

Cho

H.M.

, Bae

S.J.

, Kim

, Jeong

I.J.

, Bi-objective scheduling for reentrant hybrid flow shop using Pareto genetic algorithm[J], Computers & Industrial Engineering 61(3) (2011), 529–541.

Ying

K.C.

, Lin

S.W.

, Wan

S.Y.

, Bi-objective reentrant hybrid flowshop scheduling: an iterated Pareto greedy algorithm[J], International Journal of Production Research 52(19) (2014), 5735–5747.

Shen

, Wang

, Deng

, Zheng

, A Pareto-Based Discrete Harmony Search Algorithm for Bi-objective Reentrant Hybrid Flowshop Scheduling Problem, in: Kim

J.H.

and Geem

Z.W.

(Eds.), Harmony Search Algorithm, Springer Berlin Heidelberg, 2016, pp. 435–445.

Shen

, Wang

, Zheng

, A modified teaching–learning-based optimisation algorithm for bi-objective re-entrant hybrid flowshop scheduling[J], International Journal of Production Research 54(12) (2016), 3622–3639.

10.

Lin

C.C.

, Liu

W.Y.

, Chen

Y.H.

, Considering stockers in reentrant hybrid flow shop scheduling with limited buffer capacity[J], Computers & Industrial Engineering 139 (2020), 106154.

11.

International Energy Outlook 2016 [EB/OL], Energy Information Administration. https://www.eia.gov/ (accessed December 31, 2020.

12.

, Gao

, Li

, Pan

, Wang

, Energy-efficient permutation flow shop scheduling problem using a hybrid multi-objective backtracking search algorithm[J], Journal of Cleaner Production 144 (2017), 228–238.

13.

Jovane

, Yoshikawa

, Alting

, Boër

C.R.

, Westkamper

, Williams

, Tseng

, Seliger

, Paci

A.M.

, The incoming global technological and industrial revolution towards competitive sustainable manufacturing[J], CIRP Annals 57(2) (2008), 641–659.

14.

Wang

J.J.

, Wang

, A Knowledge-Based Cooperative Algorithm for Energy-Efficient Scheduling of Distributed Flow-Shop[J], IEEE Transactions on Systems, Man, and Cybernetics: Systems 50(5) (2020), 1805–1819.

15.

Chen

, Zhang

, Arinez

, Biller

, Energy-Efficient Production Systems Through Schedule-Based Operations[J], IEEE Transactions on Automation Science and Engineering 10(1) (2013), 27–37.

16.

Che

, Wu

, Peng

, Yan

, Energy-efficient bi-objective single-machine scheduling with power-down mechanism[J], Computers & Operations Research 85 (2017), 172–183.

17.

, Che

, Energy-efficient no-wait permutation flow shop scheduling by adaptive multi-objective variable neighborhood search[J], Omega 94 (2020), 102–117.

18.

Fang

, Uhan

N.A.

, Zhao

, Sutherland

J.W.

, Scheduling on a single machine under time-of-use electricity tariffs[J], Annals of Operations Research 238(1) (2016), 199–227.

19.

, Che

, A memetic differential evolution algorithm for energy-efficient parallel machine scheduling[J], Omega 82 (2019), 155–165.

20.

Ding

J.Y.

, Song

S.J.

, Zhang

, Chiong

, Wu

, Parallel Machine Scheduling Under Time-of-Use Electricity Prices: New Models and Optimization Approaches[J], IEEE Transactions on Automation Science and Engineering 13(2) (2016), 1138–1154.

21.

Geng

K.F.

, Ye

C.M.

, Dai

Z.H.

, Liu

, Bi-Objective Re-Entrant Hybrid Flow Shop Scheduling considering Energy Consumption Cost under Time-of-Use Electricity Tariffs[J], Complexity 2020 (2020), 1–17.

22.

Ding

, Schulz

, Shen

, Buscher

, Lü

, Energy aware scheduling in flexible flow shops with hybrid particle swarm optimization[J], Computers & Operations Research 125 (2021), 1–17.

23.

Zhou

, Jin

, Du

, Energy-efficient scheduling of a single batch processing machine with dynamic job arrival times[J], Energy 209 (2020), 118420.

24.

Zheng

, Zhou

, Xu

, Chen

, Energy-efficient scheduling for multi-objective two-stage flow shop using a hybrid ant colony optimisation algorithm[J], International Journal of Production Research 58(13) (2020), 4103–4120.

25.

Toptal

, Sabuncuoglu

, Distributed scheduling: a review of concepts and applications[J], International Journal of Production Research 48(18) (2010), 5235–5262.

26.

Gao

, Chen

, Deng

, An efficient tabu search algorithm for the distributed permutation flowshop scheduling problem[J], International Journal of Production Research 51(3) (2013), 641–651.

27.

Lin

S.W.

, Ying

K.C.

, Huang

C.Y.

, Minimising makespan in distributed permutation flowshops using a modified iterated greedy algorithm[J], International Journal of Production Research 51(16) (2013), 5029–5038.

28.

Meng

, Pan

Q.K.

, Wang

, A distributed permutation flowshop scheduling problem with the customer order constraint[J], Knowledge-Based Systems 184 (2019), 104894.

29.

Framinan

J.M.

, Perez-Gonzalez

, New approximate algorithms for the customer order scheduling problem with total completion time objective[J], Computers & Operations Research 78 (2017), 181–192.

30.

Framinan

J.M.

, Perez-Gonzalez

, Order scheduling with tardiness objective: Improved approximate solutions[J], European Journal of Operational Research 266(3) (2018), 840–850.

31.

Pan

J.C.H.

, Chen

J.S.

, Mixed binary integer programming formulations for the reentrant job shop scheduling problem[J], Computers & Operations Research 32(5) (2005), 1197–1212.

A memetic algorithm for energy-efficient distributed re-entrant hybrid flow shop scheduling problem

Abstract

Keywords

1 Introduction

2 Problem description and formulation

2.1 Problem description

2.3 Objectives functions and main constraints

Table 1 Processing times Processing times First-pass Re-entrant1 Stage 1 Stage 2 Stage 3 Stage 1 Stage 2 Stage 3 Job1 2 2 1 1 2 1 Job2 1 3 2 3 2 1 Job3 3 1 2 2 1 2 Job4 2 1 1 2 3 2 Job5 1 2 3 1 1 2 Job6 2 2 3 2 3 2

3.3.1 Rule1

3.3.2 Rule2

3.3.3 Rule3

3.4.1 Orders_Swap_among_Factories

3.4.2 Orders_Random_Insert_other_Factory

3.4.3 Job_Best_Insert_within_Factory

3.4.4 Job_Best_Swap_within_Factory

3.4.5 Jobs_Inverse_within_Factory

4 Numerical experiments and comparisons

4.1 Generating test data

Table 3 Parameters combination Parameters Factor level 1 2 3 4 PopSize 20 40 60 80 P init 0.1 0.2 0.3 0.4 P Insert 0.3 0.4 0.5 0.6 G max 60 80 100 120

Funding

References

Table 1
Processing times

Processing times First-pass Re-entrant1

Stage 1 Stage 2 Stage 3 Stage 1 Stage 2 Stage 3

Job1 2 2 1 1 2 1

Job2 1 3 2 3 2 1

Job3 3 1 2 2 1 2

Job4 2 1 1 2 3 2

Job5 1 2 3 1 1 2

Job6 2 2 3 2 3 2

Table 3
Parameters combination

Parameters Factor level

1 2 3 4

PopSize 20 40 60 80

P _init 0.1 0.2 0.3 0.4

P _Insert 0.3 0.4 0.5 0.6

G _max 60 80 100 120