Flow shop scheduling with low carbon emission and variable machining parameters

Abstract

Synergistic improvement on low carbon manufacturing is not taken full consideration in flow shop scheduling in literatures. To fill the gap, this article extends a flow shop scheduling problem by considering low carbon emission and variable machining parameters. The difference between the extended scheduling problem and the traditional one is that the former investigates the comprehensive effects of machine and scheduling levels on production throughput and environmental impact. A mathematical model of the extended problem is formulated to find the Pareto frontier, which is the set of scheduling solutions that are all Pareto efficient. Moreover, two problem-specific scheduling methods are presented and their effects on the considered problem are proved. A novel multi-objective heuristic is developed to deal with the extended problem because it is a non-deterministic polynomial acceptable problem. An actual case is studied to demonstrate the feasibility and applicability of the proposed extended problem and two scheduling methods. These finds can be employed in manufacturing execution systems and do not require additional cost.

Introduction

It is well known that global warming stems from the increasing amount of greenhouse gas due to the soaring energy consumption. In China, a National Climate Change Plan was proposed with a goal to reduce the carbon emission intensity by 40%–45% by 2020 from the 2005 level. Similar regulations were laid down in other countries and regions, including the United States, European Union, Japan, and so on.

There are two practical approaches for reducing carbon emission in manufacturing: one is to use energy-saving equipment and embodied energy framework, ¹ and the other is to employ low carbon operational strategies and optimization methods. ² The former needs a large capital investment to develop new technologies and replace old equipment, while the latter requires less cost and is more preferred by manufacturers.

Recently, extensive studies on energy saving or low carbon in manufacturing scheduling systems are found. Mouzon et al. ³ observed a fact that there were a large amount of energy savings when idle machines were turned off. In their follow-up work, a single machine scheduling model with the objectives of minimizing energy consumption and tardiness was reported. ⁴ Their model adopted an “ON–OFF” method, which switches machines off when they are idle and switches them on before new workpieces come. Fang et al. ⁵ used a speed scaling method in a flow shop scheduling problem with a peak power constraint. It is assumed that machines can run at varying speed levels and more energy is consumed if a machine works at a higher speed. Thus, it is possible to reduce energy consumption by slowing down some operations.

According to the studies,^6,7 carbon emission of any machining operation is concerned with energy consumption, resource utilization, and waste disposal. Carbon emission is a more overall criterion to indicate the environmental burden of manufacturing activities when compared with energy consumption. In the carbon efficient scheduling model proposed by Ding et al., ⁸ carbon emission was only tied to energy consumption. However, they neglected the carbon emission caused during the preparation time and tool motion time of operations. A special blocking flow shop scheduling case in turning processes for high efficiency and low carbon was investigated in Lin et al.’s ⁹ work.

Synergistic improvement can be achieved by considering the integrated optimization of multiple perspectives.^10,11 However, no much work was reported to investigate the integrated optimization. In the speed scaling method, machines only run at several prespecified levels.^5,8,12 It may not be practical in some manufacturing systems, as machine tools can run at any speed in a range. Besides, as other machining parameters (i.e. feed rate and depth of cut) also have impacts on carbon emission, ¹³ they should be considered in the integrated optimization. To fill the gap, this article addressed a flow shop scheduling problem with variable machining parameters. The contribution of this work can be summarized as follows: (1) development of a novel multi-objective mathematical model, taking into account the machining parameters and the scheduling sequence as decision variables; (2) presentation of two problem-specific scheduling methods that aim to cut down carbon emission in scheduling; and (3) demonstration of the effects of the synergistic use of methods on Pareto frontier and the importance of the order of using them in optimization procedures.

The remainder of this article is organized as follows. A carbon emission calculation method for machining operations is introduced in section “Carbon emission calculating method,” while the mathematical model of a flow shop scheduling problem with variable machining parameters is briefly described in section “Flow shop scheduling problem with variable machining parameters.” Section “Problem-specific scheduling methods” gives two scheduling methods and section “Modified multi-objective GWO” presents a modified multi-objective grey wolf optimizer (GWO). Finally, sections “Case study” and “Conclusion” present the analysis of a case study and the conclusion, respectively.

Carbon emission calculating method

As the carbon emission due to workpiece (chips) is determined at the stage of production design or process design, they are not considered during machining parameter optimization (MPO) and scheduling. Carbon emission in a machining operation, denoted by C_ms (g), is mainly caused by electricity consumption, cutting tools, and cutting fluid. ¹⁴ It can be calculated as follows

$$C_{ms}=C_{elec}+C_{tool}+C_{fluid}$$

(1)

where C_elec (g), C_tool (g), and C_fluid (g) are carbon emission due to electricity consumption, cutting tools, and cutting fluid in a machining operation, respectively.

C_elec

C_elec has a linear relationship with energy consumption, denoted by E_elec (kJ), in a machining operation, and can be determined as follows

$$C_{elec}=F_{elec}\times E_{elec}$$

(2)

where F_elec is the electricity carbon emission factor. F_elec is obtained from the study conducted by Department of National Development and Reform Commission in China in 2011. The average value of six grids in China is 0.1524 g CO₂/kJ. ¹⁵

Three states are observed in a machining operation: “basic state,”“ready state,” and “cutting state.” ¹⁶ Given a machining operation with N_step passes, E_elec can be represented as follows

$$\begin{array}{c}\multicolumn{1}{c}{E_{\mathrm{elec}}=\sum _{l=1}^{N_{\mathrm{step}}}(E_b^l+E_r^l+E_c^l)}\\ \multicolumn{1}{c}{=\sum _{l=1}^{N_{\mathrm{step}}}(\underset{0}{\overset{t_b^l}{\int }}{P}_b^l\mathrm{dt}+\underset{0}{\overset{t_r^l}{\int }}{P}_r^l\mathrm{dt}+\underset{0}{\overset{t_c^l}{\int }}{P}_c^l\mathrm{dt})}\end{array}$$

(3)

where $E_b^l$ , $E_r^l$ , and $E_c^l$ (kJ) are the electricity consumption in the basic, ready, and cutting states of the lth pass, respectively. $P_b^l$ , $P_r^l$ , and $P_c^l$ (kW) are the power consumption in the basic, ready, and cutting states of the lth pass, respectively. $t_b^l$ , $t_r^l$ , and $t_c^l$ (s) are the durations of the basic, ready, and cutting states of the lth pass, respectively.

In the basic state, a constant value of power is consumed. $E_b^l$ can be represented as follows

$$E_b^l=t_b^l\times P_b$$

(4)

where P_b is the power consumption (kW) in the basic state.

In the ready state, the power consumption, denoted by P_p (kW), consists of two parts: P_b and the operational power P_op. ¹⁶ The latter has a linear relationship with the spindle speed, denoted by n (r/min)

$$P_{\mathrm{op}}=k_{1}n+k_2$$

(5)

where k₁ and k₂ are the regression coefficients.

Then, $E_r^l$ (kJ) can be represented as follows

$$E_r^l=t_r^l\times P_r^l=t_r^l\times (P_b+k_1{n}_l+k_2)$$

(6)

Li and Kara ¹⁷ reported a specific energy consumption model to estimate $E_c^l$ . It can be applied to machining operations with both geometrically defined and undefined cutting edges

$$E_c^l=\mathrm{SE}{C}_l\times M_{\mathrm{mrv}}^l$$

(7)

where $S_{\sec }^l$ (kJ/cm³) and $M_{\mathrm{mrv}}^l$ (cm³) are the specific energy consumption and the material removal volume in the lth pass, respectively

$$\mathrm{SE}{C}_l=(r_0+\frac{r_1}{M_{\mathrm{mrr}}^l})$$

(8)

where $M_{\mathrm{mrr}}^l$ the material removal rate in the lth pass, and r₀ and r₁ are the regression coefficients.

The formulae of $M_{\mathrm{mrv}}(\mathrm{c}{\mathrm{m}}^3)$ and $M_{\mathrm{mrr}}(\mathrm{c}{\mathrm{m}}^3/\mathrm{s})$ are highly tied to machining types. Equations (9)–(11) are illustrated, for example

$$\mathrm{External}\mathrm{Turning}:M_{mrr}=\frac{vf{a}_p}{60},M_{mrv}=\frac{\pi }{4}(D_1^2-D_2^2)L$$

(9)

$$\mathrm{Internal}\mathrm{Turning}:M_{mrr}=\frac{vf{a}_p}{60},M_{mrv}=\frac{\pi }{4}(D_2^2-D_1^2)L$$

(10)

$$\mathrm{Drilling}:M_{\mathrm{mrr}}=\frac{\mathrm{vf}{d}_0}{240},M_{\mathrm{mrv}}=\pi d_0{\left(\frac{L_c}{2}\right)}^2$$

(11)

where v, f, a_p are the cutting velocity (m/min), feed rate (mm/r), and depth of cut (mm), respectively. D₁ and D₂ are the workpiece diameters before and after machining, respectively. L is the cutting length. d₀ is the tool diameter, and L_c is the depth of drilling.

Substituting C_elec in equations (3)–(7), C_elec can be re-expressed as follows

$$C_{\mathrm{elec}}=\sum _{l=1}^{N_{\mathrm{step}}}[t_b^l\times P_b+p_1{n}_l^2+p_2{n}_l+p_3+t_r^l\times (P_b+P_{\mathrm{op}})+S_{\sec }^l\times M_{\mathrm{mrv}}^l]$$

(12)

C_tool and C_fluid

Carbon emission caused by cutting tools, denoted by C_tool (g), can be calculated as follows ¹³

$$C_{\mathrm{tool}}=\sum _{l=1}^{N_{\mathrm{step}}}{\xi }_l\times F_{\mathrm{tool}}\times W_{\mathrm{tool}}$$

(13)

where ξ_l is the used tool life in the lth pass, F_tool (g CO₂/g) is the tool carbon emission factor, and W_tool (g) is the tool mass. ξ_l is given as follows

$${\xi }_l=\frac{t_c^l}{T_t^l}\times 100\%=\frac{M_{mrv}^l}{M_{mrr}^l\times T_t^l}\times 100\%$$

(14)

where $T_t^l$ is the tool life time in the lth pass. As a tool can be re-sharpened several times in its entire life time, $T_t^l(\min )$ can be calculated by

$$T_t^l=(N+1)T_{\mathrm{tool}}^l$$

(15)

where N is the times of re-sharpening and $T_{\mathrm{tool}}^l$ (min) is the tool durability in the lth pass.

Carbon emission caused by cutting fluid, denoted by C_fluid (g), can be calculated as follows ⁶

$$C_{\mathrm{fluid}}=\frac{\sum _{l=1}^{N_{\mathrm{step}}}(t_b^l+t_t^l+t_p^l+t_c^l)}{60\times T_{\mathrm{fluid}}}\times (V_{\mathrm{cc}}+V_{\mathrm{acc}})\times (F_{\mathrm{oil}}+\frac{F_{\mathrm{wc}}}{{\delta }_{\mathrm{fluid}}})$$

(16)

where F_oil (g CO₂/L) is the cutting fluid carbon emission factor, F_wc (g CO₂/L) is the carbon emission factor of cutting fluid disposal. V_cc (L) and V_acc (L) are the initial volume and additional volume of cutting fluid, respectively. δ_fluid and T_fluid (min) are the concentration and replace cycle of cutting fluid, respectively.

Flow shop scheduling problem with variable machining parameters

The flow shop scheduling problem with variable machining parameters (FVP) is an extended version of a flow shop scheduling problem (FSP). An FVP model owns a two-level structure, that is, the lower level is a MPO sub-model and the upper level is an FSP sub-model.

MPO sub-model

The optimization objectives of the MPO sub-model are machining time t_ms and carbon emission C_ms. Generally, a specific value of the depth of cut is determined by the workpiece’s process and the machining allowance, and that of the width of cut is determined by the rigidity of tools, machines, and workpieces. Therefore, the cutting velocity and feed rate are considered as decision variables in the sub-model.

Assume that operation O_ij has $n_{\mathrm{step}}^{\mathrm{ij}}$ passes. Then, $t_{\mathrm{ms}}^{\mathrm{ij}}$ is represented as follows

$$t_{\mathrm{ms}}^{\mathrm{ij}}=\sum _{l=1}^{n_{\mathrm{step}}^{\mathrm{ij}}}(t_b^{\mathrm{ijl}}+t_p^{\mathrm{ijl}}+\frac{M_{\mathrm{mrv}}^{\mathrm{ijl}}}{M_{\mathrm{mrr}}^{\mathrm{ijl}}}),i=1,2,\dots ,n_m,j=1,2,\dots ,n_{\mathrm{job}}$$

(17)

where $t_b^{\mathrm{ijl}}$ and $t_p^{\mathrm{ijl}}$ are the durations of basic and ready states of the lth pass of O_ij, respectively. $M_{\mathrm{mrv}}^{\mathrm{ijl}}$ and $M_{\mathrm{mrr}}^{\mathrm{ijl}}$ are the material removal rate and material removal volume of the lth pass, respectively. n_m is the number of machines, and n_job is the number of workpieces.

According to equations (12), (13), and (16), carbon emission of O_ij can be calculated as follows

$$\begin{array}{c}\multicolumn{1}{c}{C_{\mathrm{ms}}^{\mathrm{ij}}=F_{\mathrm{elec}}\times }\\ \multicolumn{1}{c}{\sum _{l=1}^{n_{\mathrm{step}}^{\mathrm{ij}}}[P_b^i\times t_b^{\mathrm{ijl}}+(P_r^i+P_{\mathrm{op}}^{\mathrm{ijl}})t_r^{\mathrm{ijl}}+M_{\mathrm{mrv}}^{\mathrm{ijl}}(C_0^i+\frac{C_1^i}{M_{\mathrm{mrr}}^{\mathrm{ijl}}})]}\\ \multicolumn{1}{c}{+\sum _{l=1}^{n_{\mathrm{step}}^{\mathrm{ij}}}{\xi }_{\mathrm{ijl}}\times F_{\mathrm{tool}}\times W_{\mathrm{tool}}+\frac{t_{\mathrm{ms}}^{\mathrm{ij}}}{T_{\mathrm{fluid}}^i}\times (V_{\mathrm{cc}}+V_{\mathrm{acc}})}\\ \multicolumn{1}{c}{\times (F_{\mathrm{oil}}+\frac{F_{\mathrm{wc}}}{\delta })=1,2,\dots ,n_m;j=1,2,\dots ,n_{\mathrm{job}}}\end{array}$$

(18)

It should be ensured that the machining parameters are in their upper and lower bounds, that is

$$v_{min}^{\mathrm{ij}}\le v_{\mathrm{ijl}}\le v_{max}^{\mathrm{ij}},i=1,2,\dots ,n_m;j=1,2,\dots ,n_{\mathrm{job}};l=1,2,\dots ,n_{\mathrm{step}}^{\mathrm{ij}}$$

(19)

$$f_{min}^{\mathrm{ij}}\le f_{\mathrm{ijl}}\le f_{max}^{\mathrm{ij}},i=1,2,\dots ,n_m;j=1,2,\dots ,n_{\mathrm{job}};l=1,2,\dots ,n_{\mathrm{step}}^{\mathrm{ij}}$$

(20)

where $v_{\min }^{\mathrm{ij}}(v_{max}^{\mathrm{ij}})$ is the lower (upper) bound of the cutting velocity, $f_{\min }^{\mathrm{ij}}(f_{max}^{\mathrm{ij}})$ is the lower (upper) bound of feed rate, and v_ijl (f_ijl) is the cutting velocity (feed rate) of the lth pass. The bounds are determined by referring to the machinery processing technical handbook and the technical handbook provided by the cutting tool supplier.

FSP sub-model

The optimization objectives of the FSP sub-model are makespan T_max(π) and carbon emission C_opt. The “ON–OFF” method switches machines on before their first operations and switches them off after their last operations.

Then, the details of the FSP sub-model are presented as follows:

• Optimization objectives

$$min\{\begin{array}{c}{T}_{max}\left(\pi \right)=\sum _{i=1}^{n_m-1}\sum _{j=1}^{n_{\mathrm{job}}}{x}_{j1}{t}_{\mathrm{ms}}^{\mathrm{ij}}+\sum _{j=1}^{n_{\mathrm{job}}-1}{I}_{n_m,j}\hfill \\ \multicolumn{1}{c}{C_{\mathrm{opt}}=\sum _{i=1}^{n_m}\sum _{j=1}^{n_{\mathrm{job}}}{C}_{\mathrm{ms}}^{\mathrm{ij}}+F_{\mathrm{elec}}\times \sum _{i=1}^{n_m}\sum _{j\prime =1}^{n_{\mathrm{job}}}{I}_{i,{\pi }_{j\prime }}{P}_b^i}\end{array}$$

(21)

• With constraints

$$\sum _{j=1}^{n_{\mathrm{job}}}{x}_{j,j\prime }=1,j\prime =1,2,\dots ,n_{\mathrm{job}}$$

(22)

$$\sum _{j\prime =1}^{n_{\mathrm{job}}}{x}_{j,j\prime }=1,j=1,2,\dots ,n_{\mathrm{job}}$$

(23)

$$\begin{array}{c}{I}_{i,{\pi }_{j\prime }}+\sum _{j=1}^{n_{\mathrm{job}}}{x}_{j,j\prime +1}{t}_{\mathrm{ms}}^{i,{\pi }_{j\prime }}+W_{i,{\pi }_{j\prime +1}}-W_{i,{\pi }_{j\prime }}\hfill \\ \multicolumn{1}{c}{-\sum _{j=1}^{n_{\mathrm{job}}}{x}_{j,j\prime }{t}_{\mathrm{ms}}^{i,{\pi }_{j\prime +1}}-I_{i+1,{\pi }_{j\prime }}=0}\\ j\prime =1,2,\dots ,n_{\mathrm{job}}-1;i=1,2,\dots ,n_m-1\hfill \end{array}$$

(24)

$$W_{i,{\pi }_1}=0,i=1,2,\dots ,n_m-1$$

(25)

$$I_{i,{\pi }_{j\prime }}=0,j\prime =1,2,\dots ,n_{\mathrm{job}}-1$$

(26)

where π_j′ is the j′th workpiece in π, $P_b^i$ is the power consumption of machine i in the basic state, and $I_{i,{\pi }_{j\prime }}$ is the idle time of machine i between the j′th and (j′+ 1)th workpieces in π. x_{j, j′} is a decision variable. If workpiece j is the j′th workpiece in π, x_jj′ = 1; otherwise x_jj′ = 0. W_i,j′ is the waiting time of the j′th workpiece in π between machines i and i+ 1.

It is obvious that each position in the schedule should contain just one workpiece and each workpiece should be assigned to only one position, which is supported by constraints (22) and (23). Constraint (24) is incorporated into the model to calculate the idle times.

According to the two sub-models, a feasible solution can be expressed by a schedule π and a matrix of machine parameter combinations M_mpc as follows

$$\pi =〈{\pi }_1,{\pi }_2,\dots ,{\pi }_{n_{\mathrm{job}}}〉$$

(27)

$$M_{\mathrm{mpc}}=\left[\begin{array}{cccc}{M}_{\mathrm{mpc}}(1,1)\hfill & M_{\mathrm{mpc}}(1,2)\hfill & \cdots \hfill & M_{\mathrm{mpc}}(1,n_{\mathrm{job}})\hfill \\ \multicolumn{1}{c}{M_{\mathrm{mpc}}(2,1)}& M_{\mathrm{mpc}}(2,1)\hfill & \cdots \hfill & M_{\mathrm{mpc}}(2,n_{\mathrm{job}})\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{M_{\mathrm{mpc}}(n_m,1)}& M_{\mathrm{mpc}}(n_m,2)\hfill & \cdots \hfill & M_{\mathrm{mpc}}(n_m,n_{\mathrm{job}})\hfill \end{array}\right]$$

(28)

where the elements in π are the subscripts of workpieces, and M_mpc(i, j) is a combination of machining parameters for O_ij.

Problem-specific scheduling methods

Two problem-specific scheduling methods are proposed in this article: postponing and improved two-stage optimization methods.

Postponing method

A postponing method is designed to postpone the processing of a workpiece as late as possible. Let T(π_j′, i) and S(π_j′, i) denote the completion time and start time of $O_{i,{\pi }_{j\prime }}$ , respectively. Given a specific solution x = (π, M_mpc), T(π_j′,i) can be calculated as follows

$$T({\pi }_1,i)=\sum _{l=1}^i{t}_{\mathrm{ms}}^{l,1},i=1,2,\dots ,n_m$$

(29)

$$T({\pi }_{j\prime },1)=\sum _{l=1}^{j\prime }{t}_{\mathrm{ms}}^{1,{\pi }_l},j\prime =1,2,\dots ,n_{\mathrm{job}}$$

(30)

$$T({\pi }_{j\prime },i)=max(T({\pi }_{j\prime },i-1),T({\pi }_{j\prime -1},i))+t_{\mathrm{ms}}^{i,{\pi }_{j\prime }},i=2,3,\dots ,n_m;j\prime =1,2,3,\dots ,n_{\mathrm{job}}$$

(31)

S(π_j′, i) is equal to $T({\pi }_{j\prime },i)-t_{\mathrm{ms}}^{i,{\pi }_{j\prime }}$ . The total idle time of machine i, denoted by $t_{\mathrm{idle}}^i$ , can be expressed as follows

$$t_{\mathrm{idle}}^i=T({\pi }_{n_{\mathrm{job}}},i)-S({\pi }_1,i)-\sum _{j=1}^{n_j}{t}_{\mathrm{ms}}^{\mathrm{ij}}$$

(32)

According to equation (32), $t_{\mathrm{idle}}^i$ depends on $T({\pi }_{n_j},i)$ and S(π₁,i) when M_mpc is fixed, and C_opt can be represented as follows

$$C_{\mathrm{opt}}=\sum _{i=1}^{n_m}\sum _{j=1}^{n_{\mathrm{job}}}{C}_{\mathrm{ms}}^{\mathrm{ij}}+F_{\mathrm{elec}}\times \sum _{i=1}^{n_m}{t}_{\mathrm{idle}}^i\times P_b^i$$

(33)

Let $\hat{T}({\pi }_{j\prime },i)$ and $\hat{S}({\pi }_{j\prime },i)$ denote the completion time and start time of $O_{i,{\pi }_{j\prime }}$ when the postponing method is used, respectively. Then, $\hat{T}({\pi }_{j\prime },i)$ and $\hat{S}({\pi }_{j\prime },i)$ are re-expressed as follows

$$\hat{T}({\pi }_{n_{\mathrm{job}}},i)=T({\pi }_{n_{\mathrm{job}}},i),i=n_m,n_m-1,\dots ,1$$

(34)

$$\hat{T}({\pi }_{j\prime },n_m)=\hat{T}({\pi }_{j\prime +1},n_m)-t_{\mathrm{ms}}^{n_m,{\pi }_{j\prime +1}},j\prime =n_{\mathrm{job}}-1,\dots ,1$$

(35)

$$\begin{array}{c}\hat{T}({\pi }_j,i)=min(\hat{T}({\pi }_{j\prime },i+1)-t_{\mathrm{ms}}^{i+1,{\pi }_{j\prime }},\hat{T}({\pi }_{j\prime +1},i)-t_{\mathrm{ms}}^{i,{\pi }_{j\prime +1}})\hfill \\ i=n_m-1,\dots ,1;j\prime =n_{\mathrm{job}}-1,\dots ,1\hfill \end{array}$$

(36)

$$\hat{S}({\pi }_{j\prime },i)=\hat{T}({\pi }_{j\prime },i)-t_{\mathrm{ms}}^{i,{\pi }_{j\prime }},i=1,2,\dots ,n_m;j\prime =1,2,\dots ,n_{\mathrm{job}}$$

(37)

Property 1

If the postponing method is employed on a solution, the original solution will be dominated by or equivalent to the improved solution.

Proof

The following equations are inferred from equations (29)–(31)

$$\begin{array}{c}T({\pi }_{j\prime },i)\le min\hfill \\ \multicolumn{1}{c}{(T({\pi }_{j\prime },i+1)-t_{\mathrm{ms}}^{i+1,{\pi }_{j\prime }},T({\pi }_{j\prime +1},i)-t_{\mathrm{ms}}^{i,{\pi }_{j\prime +1}})}\\ i=1,2,\dots ,n_m-1;j\prime =1,2,\dots ,n_{\mathrm{job}}-1\hfill \\ T({\pi }_{n_{\mathrm{job}}},i)\le T({\pi }_{n_{\mathrm{job}}},i+1)-t_{\mathrm{ms}}^{i+1,{\pi }_{n_{\mathrm{job}}}},\hfill \\ \multicolumn{1}{c}{i=1,2,\dots ,n_m-1}\\ T({\pi }_{j\prime },n_m)\le T({\pi }_{j\prime +1},n_m)-t_{\mathrm{ms}}^{n_m,{\pi }_{j\prime +1}},\hfill \\ \multicolumn{1}{c}{j\prime =1,2,\dots ,n_{\mathrm{job}}-1}\end{array}$$

Thus

$$\begin{array}{c}\hat{T}({\pi }_{n_{\mathrm{job}}-1},n_m)=\hat{T}({\pi }_{n_{\mathrm{job}}},n_m)-t_{\mathrm{ms}}^{n_m,{\pi }_{\mathrm{job}}}\hfill \\ =T({\pi }_{n_{\mathrm{job}}},n_m)-t_{\mathrm{ms}}^{n_m,{\pi }_{\mathrm{job}}}\hfill \\ \ge T({\pi }_{n_{\mathrm{job}}-1},n_m)\hfill \end{array}$$

Then, it can be inferred that

$$\begin{array}{c}\widehat{T}\left({\pi }_{j^{\prime }},n_m\right)=\widehat{T}\left({\pi }_{j^{\prime }+1},n_m\right)-t_{ms}^{n_m,{\pi }_{j^{\prime }+1}}\\ \ge T\left({\pi }_{j^{\prime }+1},n_m\right)-t_{ms}^{n_m,{\pi }_{j^{\prime }+1}},j^{\prime }=1,2,\dots ,n_{job}-2\\ \ge C\left({\pi }_{j^{\prime }},n_m\right)\end{array}$$

Similarly, it can be proved that $\hat{T}({\pi }_{j\prime },i)\ge T({\pi }_{j\prime },i)$ (j’ = 1, 2, …, n_job; i = 1, 2, …, n_m– 1). Thus

$$\begin{array}{c}\hat{S}({\pi }_{j\prime },i)=\hat{S}({\pi }_{j\prime },i)-t_{\mathrm{ms}}^{i,{\pi }_{j\prime }}\hfill \\ \ge T({\pi }_{j\prime },i)-t_{\mathrm{ms}}^{i,{\pi }_{j\prime }},i=1,2,\dots ,n_m;j\prime =1,2,\dots j\prime ,n_{\mathrm{job}}\hfill \\ =S({\pi }_{j\prime },i)\hfill \end{array}$$

Therefore, the following equation can be inferred

$$\begin{array}{c}{\hat{C}}_{\mathrm{opt}}=\sum _{i=1}^{n_m}\sum _{j=1}^{n_{\mathrm{job}}}{C}_{\mathrm{ms}}^{\mathrm{ij}}+F_{\mathrm{elec}}\times \sum _{i=1}^{n_m}{\hat{t}}_{\mathrm{idle}}^i\times P_b^i\hfill \\ =\sum _{i=1}^{n_m}\sum _{j=1}^{n_{\mathrm{job}}}{C}_{\mathrm{ms}}^{\mathrm{ij}}+F_{\mathrm{elec}}\hfill \\ \multicolumn{1}{c}{\times \sum _{i=1}^{n_m}(\hat{T}({\pi }_{n_{\mathrm{job}}},i)-\hat{S}({\pi }_1,i)-\sum _{j=1}^{n_j}{t}_{\mathrm{ms}}^{\mathrm{ij}})\times P_b^i}\\ \ge \sum _{i=1}^{n_m}\sum _{j=1}^{n_{\mathrm{job}}}{C}_{\mathrm{ms}}^{\mathrm{ij}}+F_{\mathrm{elec}}\hfill \\ \multicolumn{1}{c}{\times \sum _{i=1}^{n_m}(T({\pi }_{n_{\mathrm{job}}},i)-S({\pi }_1,i)-\sum _{j=1}^{n_j}{t}_{\mathrm{ms}}^{\mathrm{ij}})\times P_b^i}\\ =\sum _{i=1}^{n_m}\sum _{j=1}^{n_{\mathrm{job}}}{C}_{\mathrm{ms}}^{\mathrm{ij}}+F_{\mathrm{elec}}\times \sum _{i=1}^{n_m}{t}_{\mathrm{idle}}^i\times P_b^i\hfill \\ =C_{\mathrm{opt}}\hfill \end{array}$$

According to equation (34), ${\hat{T}}_{max}(\pi )\equiv T_{max}(\pi )$ . If ${\hat{C}}_{\mathrm{opt}}$ is smaller than C_opt, then the original solution is dominated by the improved solution. Otherwise, the former is equivalent to the latter.

Improved two-stage optimization method

As machining parameters can be set to any values in their bounds, there exist a large amount of possible combinations of them. Fortunately, all combinations are dominated by the Pareto-optimal combinations in the Pareto frontier, denoted by Γ. Moreover, Γ can be easily obtained by optimizing the MPO sub-model. In Γ, the machining time of an operation decreases but the carbon emission increases, and vice versa. Based on the fact, the following property can be proved.

Property 2

Given a solution x ₁ using the Pareto-optimal combinations. If the makespan and the starting time are fixed and a certain machining parameter combination in M_mpc is replaced with which owns a larger value of t_ms and a smaller value of C_opt, then the original solution is dominated by the new solution.

Proof

Let x ₂ and $O_{k,{\pi }_l}$ denote the new solution and the manipulated operation, respectively. The original combination is denoted by $M_{\mathrm{mpc}}(k,{\pi }_l)$ , and the replaced combination is ${\tilde{M}}_{\mathrm{mpc}}(k,{\pi }_l)$ . The carbon emission and machining time of $O_{k,{\pi }_l}$ and the idle time of machine k between the lth and (l+1)th workpieces in π when using ${\tilde{M}}_{\mathrm{mpc}}(k,{\pi }_l)$ are denoted by ${\tilde{C}}_{\mathrm{ms}}^{k,{\pi }_l}$ , ${\tilde{t}}_{\mathrm{ms}}^{k,{\pi }_l}$ , and ${\tilde{I}}_{k,{\pi }_l}$ , respectively. As ${\tilde{M}}_{\mathrm{mpc}}(k,{\pi }_l)$ owns a larger value of t_ms and a smaller value of K_ms, it follows that

$${\tilde{C}}_{\mathrm{ms}}^{k,{\pi }_l}<C_{\mathrm{ms}}^{k,{\pi }_l}\mathrm{and}{\tilde{t}}_{\mathrm{ms}}^{k,{\pi }_l}>t_{\mathrm{ms}}^{k,{\pi }_l}$$

Therefore

$$\begin{array}{l}{\tilde{I}}_{n_m,{\pi }_l}=S\left({\pi }_{l+1},n_m\right)-S\left({\pi }_l,n_m\right)-{\tilde{t}}_{ms}^{n_m,{\pi }_l}<S\left({\pi }_{l+1},n_m\right)-S\left({\pi }_l,n_m\right)-t_{ms}^{n_m,{\pi }_l}=I_{n_m,{\pi }_l}\\ l=1,2,\dots ,n_{job}-1\\ {\tilde{I}}_{k,{\pi }_l}=\min \left(S\left({\pi }_{l+1},k\right),S\left({\pi }_l,k+1\right)\right)-S\left({\pi }_l,k\right)-{\tilde{t}}_{ms}^{k,{\pi }_l}\\ <\min \left(S\left({\pi }_{l+1},k\right),S\left({\pi }_l,k+1\right)\right)-S\left({\pi }_l,k\right)-t_{ms}^{k,{\pi }_l}=I_{k,{\pi }_l}\\ k=1,2,\dots ,n_m-1;l=1,2,\dots ,n_{job}-1\end{array}$$

Thus, $C_{\mathrm{opt}}(x_2)$ is smaller than $C_{\mathrm{opt}}({\mathit{x}}_1)$ as follows

$$\begin{array}{c}{C}_{\mathrm{opt}}\left({\mathit{x}}_2\right)=\sum _{i=1}^{n_m}\sum _{j=1}^{n_j}{C}_{\mathrm{ms}}^{\mathrm{ij}}+{\tilde{C}}_{\mathrm{ms}}^{k,{\pi }_l}-C_{\mathrm{ms}}^{k,{\pi }_l}+F_{\mathrm{elec}}\hfill \\ \multicolumn{1}{c}{\times [\sum _{i=1}^{n_m}\sum _{j\prime =1}^{n_j}{I}_{i,{\pi }_{j\prime }}{P}_b^i+({\tilde{I}}_{k,{\pi }_l}-I_{k,{\pi }_l})\times P_b^k]}\\ <\sum _{i=1}^{n_m}\sum _{j=1}^{n_j}{C}_{\mathrm{ms}}^{\mathrm{ij}}+F_{\mathrm{elec}}\times \sum _{i=1}^{n_m}\sum _{j\prime =1}^{n_j}{I}_{i,{\pi }_{j\prime }}{P}_b^i\hfill \\ =C_{\mathrm{opt}}\left({\mathit{x}}_1\right)\hfill \end{array}$$

As the makespan is fixed, it can be concluded that x ₁ is dominated by x ₂ .

Based on Property 2, the steps of the improved two-stage optimization method are illustrated as follows:

• Step 1. Get the Pareto frontier of each MPO sub-model and sort the Pareto-optimal combinations in an ascending order of their machining time.

• Step 2. Calculate the maximum machining time increment of each operation without affecting the makespan, denoted by Δ. The calculation formulae are

$$\begin{array}{c}\Delta ({\pi }_{n_{\mathrm{job}}},i)=T_{max}({\pi }_{n_{\mathrm{job}}},i+1)-T({\pi }_{n_{\mathrm{job}}},i),\hfill \\ \multicolumn{1}{c}{i=1,2,\dots ,n_m-1}\\ \Delta ({\pi }_{j\prime },n_m)=S({\pi }_{j\prime +1},n_m)-T({\pi }_{j\prime },n_m),\hfill \\ \multicolumn{1}{c}{j\prime =1,2,\dots ,n_{\mathrm{job}}-1}\\ \Delta ({\pi }_{j\prime },i)=min(S({\pi }_{j\prime +1},i),S({\pi }_{j\prime },i+1))-T({\pi }_{j\prime },i),\hfill \\ \multicolumn{1}{c}{i=1,2,\dots ,n_m-1,j\prime =1,2,\dots ,n_{\mathrm{job}}-1}\end{array}$$

where $\Delta ({\pi }_{j\prime },i)$ is the Δ of $O_{i,{\pi }_{j\prime }}$ .

• Step 3. Identify the Pareto-optimal combination in the Pareto frontier which owns the minimum carbon emission and satisfies the following equation for each operation

$$t_{\mathrm{ms}}^{i,{\pi }_{j\prime }}\left({\stackrel{·}{M}}_{\mathrm{mpc}}\right)\le t_{\mathrm{ms}}^{i,{\pi }_{j\prime }}\left(M_{\mathrm{mpc}}(i,{\pi }_{j\prime })\right)+\Delta ({\pi }_{j\prime },i)$$

where $M_{\mathrm{mpc}}(i,{\pi }_{j\prime })$ and ${\stackrel{·}{M}}_{\mathrm{mpc}}$ are the current and identified combinations of $O_{i,{\pi }_{j\prime }}$ , respectively.

• Step 4. The current combinations are replaced with the identified combinations, then the machining time and carbon emission are calculated for the new solution.

Modified multi-objective GWO

GWO is a novel nature-inspired population-based meta-heuristics, which is inspired by the hunting activities of grey wolves and has been proved to be of effectiveness and efficiency. ¹⁸ Furthermore, GWO is more self-adaptive, as no pre-determined parameters are required in it except the population size and maximum number of loops. A discrete multi-objective GWO (D-MOGWO) is described as follows.

Largest order-value-based encoding scheme

As π is a sequence of discrete values and M_mpc is a matrix of continuous values, a suitable mapping should be defined between π and a continuous-value sequence. Thus, a largest order-value rule ¹⁹ is proposed to fill the gap.

Non-dominated solution archive

The optimization goal of multi-objective problems is to find the entire Pareto-optimal set or its subset. In order to reach the goal, a non-dominated solution archive ²⁰ is integrated in D-MOGWO. The archive is constructed at the beginning of each loop and involves two parts: the archive controller and adaptive grid procedure.

The archive controller is used to construct the archive. For each solution x_i in the population, the following manipulation is performed: if the archive is empty, x_i is added to it immediately; otherwise, x_i is compared with each solution x_j in the archive. If x_i dominates x_j , x_j is abandoned; if x_i cannot be dominated by any solution in the archive, it is added into the archive. As the archive has a fixed size, the adaptive grid procedure is employed when the archive is full.

The adaptive grid procedure is used to strengthen the diversity of the non-dominated solutions in the archive. The kernel idea is to divide the objective function grid of the non-dominated solutions into hyper-cubes. The procedure is re-performed when the new inserted solution lies outside the current grid. The crowding levels of hyper-cubes depend on the number of solutions in them. When the archive is full, a randomly selected solution in the most crowding hyper-cube is abandoned.

Roulette-wheel-based solution selection approach

In the basic GWO, the fittest solution in the current population is denoted by the alpha (α) wolf, and the second and third fittest solutions are denoted by the beta (β) and delta (δ) wolves, respectively. As no traditional fittest solution is available in a multi-objective optimization problem, a roulette-wheel-based solution selection approach is required in D-MOGWO. Each hyper-cube in the grid is assigned a selection probability p_i as follows

$$p_i=\frac{\frac{1}{h_i}}{\sum _{j=1}^{n_{\mathrm{hc}}}\frac{1}{h_j}}$$

(38)

where h_i is the number of solutions lie in the ith hyper-cubes and n_hc is the number of hyper-cubes.

A roulette-wheel procedure is used to choose three hyper-cubes, and then a solution is randomly selected from each hyper-cube.

The rest in the population are conceived as omega (ω) wolves, and they follow the α, β, δ wolves in the hunting.

Main loop of D-MOGWO

Following equations are used to resemble the hunting activities ¹⁸

$${\mathit{D}}_{\mathit{\alpha }}=\mathrm{abs}({\mathit{c}}_1·{\mathit{x}}_{\mathit{\alpha }}^{\mathit{t}}-{\mathit{x}}_{\mathit{i}}^{\mathit{t}}),i=1,2,\dots ,n_p$$

(39)

$${\mathit{D}}_{\mathit{\beta }}=\mathrm{abs}({\mathit{c}}_2·{\mathit{x}}_{\mathit{\beta }}^{\mathit{t}}-{\mathit{x}}_{\mathit{i}}^{\mathit{t}}),i=1,2,\dots ,n_p$$

(40)

$${\mathit{D}}_{\mathit{\delta }}=\mathrm{abs}({\mathit{c}}_3·{\mathit{x}}_{\mathit{\delta }}^{\mathit{t}}-{\mathit{x}}_{\mathit{i}}^{\mathit{t}}),i=1,2,\dots ,n_p$$

(41)

$${\mathit{x}}_{\mathit{1}}^{\mathit{t}}={\mathit{x}}_{\mathit{\alpha }}^{\mathit{t}}-{\mathit{A}}_1·{\mathit{D}}_{\alpha }$$

(42)

$${\mathit{x}}_{\mathit{2}}^{\mathit{t}}={\mathit{x}}_{\mathit{\beta }}^{\mathit{t}}-{\mathit{A}}_2·{\mathit{D}}_{\mathit{\beta }}$$

(43)

$${\mathit{x}}_{\mathit{3}}^{\mathit{t}}={\mathit{x}}_{\mathit{\delta }}^{\mathit{t}}-{\mathit{A}}_3·{\mathit{D}}_{\mathit{\delta }}$$

(44)

$${\mathit{x}}_{\mathit{i}}^{\mathit{t}\prime 1}=\frac{{\mathit{x}}_1^{\mathit{t}}+{\mathit{x}}_2^{\mathit{t}}+{\mathit{x}}_3^{\mathit{t}}}{3}$$

(45)

where ${\mathit{x}}_{\mathit{\alpha }}^{\mathit{t}}$ , ${\mathit{x}}_{\mathit{\beta }}^{\mathit{t}}$ , ${\mathit{x}}_{\mathit{\delta }}^{\mathit{t}}$ are the α, β, δ wolves in the tth loop, respectively. ${\mathit{x}}_{\mathit{i}}^{\mathit{t}}$ is the ith solution in the population in the tth loop, and ${\mathit{x}}_{\mathit{i}}^{\mathit{t}+1}$ is the (i+1)th solution in next loop. n_p is the population size. A_i and c_i (i = 1, 2, 3) are the coefficient vectors.

A ₁ and c_i can be calculated as follows

$${\mathit{A}}_{\mathit{i}}=2a\times {\mathit{r}}_1-a,i=1,2,3$$

(46)

$${\mathit{c}}_{\mathit{i}}=2·{\mathit{r}}_2,i=1,2,3$$

(47)

where the elements of r ₁ and r ₂ are random values in [0, 1]. Parameter a can be calculated as follows

$$a=2-\frac{2}{n_l}$$

(48)

where n_l is the maximum number of loops.

The flow chart of the D-MOGWO is shown in Figure 1.

Figure 1.

Flow chart of D-MOGWO.

Case study

Consider an eight-workpiece and three-machine FVP case. The diagram of the workpieces is shown in Figure 2. As shown in Table 1, eight different workpieces resulted from the combination of two different outer/inner diameters, two different lengths, and two different numbers of holes. Each workpiece requires three machining operations (i.e. external turning, internal turning, and drilling), and these operations are performed in a CNC lathe CK60 (M₁), a CNC lathe CK0628 (M₂), and a CNC drilling machine (M₃) successively.

Table 1.

Dimensions of workpieces.

Workpiece no.	∅D/∅d0 (mm)	L (mm)	∅d2 (mm)	# of holes
J 1	46/15	20	5	4
J 2	46/15	20		8
J 3	46/15	30		4
J 4	46/15	30		8
J 5	66/20	20		4
J 6	66/20	20		8
J 7	66/20	30		4
J 8	66/20	30		8

Figure 2.

Diagram of workpieces.

The formulae of $T_{\mathrm{tool}}^l$ are presented as follows

$$\mathrm{Turning}:T_{\mathrm{tool}}=60\times \sqrt[m]{\frac{C_v}{{\mathrm{va}}_p^{x_v}{f}^{y_v}}}$$

(49)

$$\mathrm{Drilling}:T_{\mathrm{tool}}=60\times \sqrt[m]{\frac{C_v{d}_0^{q_v}}{v{f}^{y_v}}}$$

(50)

where C_v, x_v, y_v, and m are the empirical coefficients related to the tool and workpiece materials.

The energy consumption regression models of these machines were established based on experimental data. The three-level orthogonal array method was used, and nine tests were need for the experiment. Then, the energy consumption regression models were established using the multiple linear regression method and are depicted in Table 2. Figure 3 gives the regression result of the SEC model of M₂. It shows that the empirical model was able to successfully predict the SEC. P_b was obtained using the sliding filter method and the values of P_b of these machines are 0.374, 0.418, and 0.473 kW, respectively. For details of the experiment, readers can refer to the authors’ previous works.^14,21

Table 2.

Coefficients of energy consumption models.

Machine no.	Regression model	Adjusted R2	RMSE
M1	Pop = 1.82 × 10−3n+ 0.037 Ssec = 2.86 + 0.40/Mmrr	0.99 0.99	0.08 0.16
M2	Pop = 3.32 × 10−4n+ 0.042 Ssec = 1.25 + 1.18/Mmrr	0.98 0.98	0.05 0.04
M3	Pop = 6.41 × 10−5n+ 0.61 Ssec = 0.85 + 0.78/Mmrr	0.95 0.98	0.10 0.09

Figure 3.

The SEC model of M_2.

Details of the operations and machining parameters are given in Table 3.^22,23 All of the machines are equipped with cemented carbide tools, the workpiece material is carbon steel, and the related constants and coefficients are listed in Table 4.^22,23 The value of F_tool is 29.6 g CO₂/g. ⁶ The period-related coefficients are given in Table 5. Both of the lathes use cutting fluid, and the drilling machine does not. The cutting fluid–related coefficients are listed in Table 6, and the values of F_oil and F_wc are 2.85 and 0.2 g CO₂/L, respectively.

Table 3.

Machining operations.

Machines	Pass no.	vmin/vmax (m/min)	fmin/fmax (mm/r)	ap (mm)
M1	1	76/135	0.3/0.5	4
	2	80/150	0.2/0.4	2
	3	84/160	0.1/0.3	1
M2	1	40/80	0.4/0.8	4
	2	50/95	0.3/0.6	2
	3	60/110	0.15/0.4	0.5
M3	\	20/40	0.1/0.3	\

Table 4.

Tool-related coefficients.

Machine no.	Wtool (g)	f (mm/r)	Cv	xv	yv	m	qv	N
M1	14	f≤ 0.30.3 < f≤ 0.7	291242	0.15	0.35	0.200.35	−	1
M2	15	f≤ 0.30.3 < f≤ 0.7	291242	0.15	0.35	0.200.35	−	1
M3	20	−	8	0	0.7	0.2	0.4	1

Table 5.

Periods of basic and ready states.

tb/tr (s)		J 1	J 2	J 3	J 4
M1	1	10/10.5	10/10.5	10/15.5	10/15.5
	2	5/5.5	5/5.5	5/10.5	5/10.5
	3	5/5.5	5/5.5	5/10.5	5/10.5
M2	1	10/12	10/12	10/17	10/17
	2	3/7	3/7	3/12	3/12
	3	3/7	3/7	3/12	3/12
M3	−	15/17	15/32	15/17	15/32
		J 5	J 6	J 7	J 8
M1	1	15/15.5	15/15.5	15/20.5	15/20.5
	2	7.5/8	7.5/8	7.5/15.5	7.5/15.5
	3	7.5/8	7.5/8	7.5/15.5	7.5/15.5
M2	1	15/17	15/17	15/22	15/22
	2	5/9.5	5/9.5	5/17	5/17
	3	5/9.5	5/9.5	5/17	5/17
M3	−	15/17	15/32	15/17	15/32

Table 6.

Cutting fluid–related coefficients.

Machine no.	Vcc (L)	Vacc (L)	δfluid	Tfluid (min)
M1	8.5	4.5	0.05	1440
M2	8.0	4.0	0.05	1440

Algorithm setting and performance metrics

D-MOGWO was coded in MATLAB R2015b in a personal computer using a 64-bit Windows 7 operating system, a Pentium Dual-Core processor, and a 4-GB RAM. The population size and the maximum number of loops are set to 500 and 200, respectively.

In order to evaluate the effectiveness of the proposed methods, following performance metrics are considered:

1. Coverage metric: ²⁴ this criterion provides a binary measure that determines which of two Pareto frontier Γ₁ and Γ₂ is better. It counts the percentage of solutions in Γ₁ that are dominated by those in Γ₂ and can be calculated as follows

$$P_{\mathrm{cm}}({\mathit{\Gamma }}_1,{\mathit{\Gamma }}_2)=\frac{|\{{\mathit{x}}_{\mathit{b}}\in {\mathit{\Gamma }}_2|\exists {\mathit{x}}_{\mathit{a}}\in {\mathit{\Gamma }}_1:{\mathit{x}}_{\mathit{a}}\succ {\mathit{x}}_{\mathit{b}}\}|}{\left|{\mathit{\Gamma }}_2\right|}$$

(51)

where x_a is a non-dominated solution in Γ₁ and x_b is a solution in Γ₂. ${\mathit{x}}_{\mathit{a}}\succ {\mathit{x}}_{\mathit{b}}$ denotes that x_a dominates x_b . Γ₁ is better than Γ₂ only if P_cm(Γ₁, Γ₂) is greater than P_cm(Γ₂, Γ₁).

• 2. Spacing metric: ²⁵ this criterion measures the extent of spread in the obtained Pareto frontier Γ. Let d_i denote Euclidean distance between solution i and its nearest solution. The spacing metric of Γ can be calculated as follows

$$P_{\mathrm{sm}}=\sqrt{\frac{\sum _{i=1}^{\left|\Gamma \right|}(d_i-\overline{d})}{\left|\Gamma \right|}}$$

(52)

where P_sm is the value of spacing metric, |Γ| is the number of solutions in Γ, and $\overline{d}$ is the average value of these distances. A smaller value of P_sm suggests that the non-dominated solutions in Γ are more evenly distributed.

D-MOGWO versus NSGA-II

To evaluate the effectiveness of the proposed heuristic, D-MOGWO is compared with NSGA-II, which is one of the most popular heuristic. The latter uses a non-dominated sorting and crowding distance and crowding distance calculating procedures. For NSGA-II, the population size is 500, the maximum number of loops is 200, the crossover rate is 0.9, and the mutation rate is 0.1. Each heuristic runs 30 times.

The performance metrics of the Pareto frontiers obtained by two heuristics are depicted in Figure 4. The part labels were described in sub-section “Algorithm settings and performance metrics”. For example, Pcm(NSGA-II, D-MOGWO) denotes the coverage metric of the Pareto frontier obtained by NSGA-II over that obtained by D-MOGWO. It is observed from the figure that the D-MOGWO outperforms NSGA-II for both performance metrics. The average of P_cm(D-MOGWO, NSGA-II) is 0.526, which means the majority of solutions in the Pareto frontier obtained by NSGA-II are dominated by D-MOGWO. The minimum, maximum, and average of P_sm(D-MOGWO) are smaller than those of P_sm(NSGA-II), respectively. The adoption of the roulette-wheel-based solution selection approach can improve the quality of population since the search process is guided by three randomly selected solutions. Therefore, it can be concluded that the D-MOGWO is effective in solving the proposed problem.

Figure 4.

Performance metrics of the Pareto frontiers obtained by two heuristics.

Analysis of problem-specific scheduling methods

For the sake of simplicity, the postponing, two-stage optimization, ⁹ and improved two-stage optimization methods are denoted by A, B, and C, respectively. Several tests are generated due to using different methods in the considered case.

Effectiveness of methods

Let Γ_B (Γ_c) denote the obtained Γ of the test using B (C). The performance metrics of the two Pareto frontiers are presented in Table 7. In the table, the minimum, maximum, and average of P_cm(Γ_B, Γ_C) are smaller than those of P_cm(Γ_C, Γ_B), respectively. This result indicates that the values of P_cm(Γ_B, Γ_C) are always smaller than those of P_cm(Γ_C, Γ_B). The average of P_cm(Γ_C, Γ_B) is larger than 0.5, whereas that of P_cm(Γ_B, Γ_C) is smaller than 0.2, which means that more than half of non-dominated solutions in Γ_B are dominated by those solutions in Γ_C in most cases. For the comparison between the spacing metrics of Γ_B and Γ_c, the results are quite interesting. Γ_C outperforms Γ_B in the average and maximum values of P_sm, while both Pareto frontiers are equal in the minimum value of P_sm. It can be inferred from this result that the solutions in Γ_C are more evenly distributed than those in Γ_B in most cases, and the best extents of spread of both Pareto frontiers are approximate to each other. These situations arise because some operations use machining parameter combinations which own smaller carbon emission and larger machining time without affecting the makespan of a schedule when using C during the optimization procedure. Therefore, it can be concluded that better Pareto frontier can be obtained when C is employed in the considered case instead of B.

Table 7.

Performance metrics of Γ_B and Γ_C.

Statistics	Pcm(ΓB, ΓC)	Pcm(ΓC, ΓB)	Psm(ΓB)	Psm(ΓC)
Minimum	0.000	0.320	0.046	0.040
Maximum	0.601	1.000	0.528	0.390
Average	0.200	0.664	0.167	0.140
Variance	0.017	0.021	0.020	0.010

Let Γ_no denote the obtained Γ of the test using no methods, and Γ_A denote the obtained Γ of the test using A. Table 8 reports the performance metrics of Γ_no, Γ_A, and Γ_C. The fact that the average of P_cm(Γ_no, Γ_A) is 0.921 shows that about 10% solutions in Γ_no are not dominated by those in Γ_A. Moreover, the minimum, maximum, and average values of P_cm(Γ_A, Γ_C) are much smaller than those of P_cm(Γ_C, Γ_A), respectively. This result indicates that the majority of solutions in Γ_C are dominated by the solutions in Γ_A. P_sm(Γ_no) (P_sm(Γ_C)) is slightly smaller than P_sm(Γ_A) (P_sm(Γ_no)). It indicates that C helps to improve the extent of spread in Pareto frontiers since the adoption of C uses the Pareto-optimal combinations instead of all combinations.

Table 8.

Performance metrics of Γ_no, Γ_A, and Γ_c.

Statistics	Pcm(Γno,ΓA)	Pcm(ΓC,ΓA)	Pcm(ΓA,ΓC)	Psm(Γno)	Psm(ΓA)	Psm(ΓC)
Minimum	0.686	0.657	0.000	0.026	0.041	0.040
Maximum	1.000	0.977	0.053	0.965	2.065	0.390
Average	0.921	0.843	0.015	0.169	0.347	0.140
Variance	0.011	0.006	0.001	0.031	0.311	0.010

Comparison with the synergistic use of methods

Two tests with both of A and C are employed are generated due to the order of using these methods. Let Γ_AC (Γ_CA) denote the obtained Γ of the test with A used before (after) C. Figure 5 shows the comparisons of the performance metrics of Γ_A, Γ_C, Γ_AC, and Γ_CA. It can be observed from the figure that no solutions in Γ_AC are dominated by those in Γ_A or Γ_C, while almost all of solutions in Γ_A and Γ_C are dominated by those in Γ_AC. Therefore, synergistic improvement on Pareto frontiers can be achieved when both of A and C are employed in the considered case. For the comparisons of the performance metrics of Γ_AC and Γ_CA, the results are quite interesting. A quite difference is observed between the mean values of P_cm(Γ_AC, Γ_CA) and P_cm(Γ_CA, Γ_AC), whereas a slight difference is observed between P_sm(Γ_AC) and P_sm(Γ_AC).

Figure 5.

Performance metrics of Γ_A, Γ_C, Γ_AC, and Γ_CA.

In order to make the analysis more statistically convincing, two-sample t tests were conducted to determine whether the performance metrics of Γ_AC and Γ_CA are significantly different from each other. The significance level used in the t tests is 5%, which means that the null hypothesis is rejected if the corresponding p value is smaller than 0.05. The two-sample t test is conducted, and the p values of t test (P_cm(Γ_AC, Γ_CA), P_cm(Γ_CA, Γ_AC)) and t test(P_sm(Γ_AC), P_sm(Γ_AC)) are 0.008 and 0.461, respectively. The results clearly show that P_cm(Γ_AC, Γ_CA) and P_cm(Γ_CA, Γ_AC) are significantly different from each other, while P_sm(Γ_AC) and P_sm(Γ_CA) are not. It can be summarized from the above experiments that the order of employing methods affects Pareto frontiers, and better Pareto frontiers can be achieved using A after C in optimization procedures.

Conclusion

This article proposes the mathematical model of a flow shop scheduling problem with variable machining parameters, which considers production efficiency and carbon emission. Two problem-specific scheduling methods (i.e. postponing and improved two-stage optimization methods) are presented and their properties are proved. Then, the steps of employing these methods based on the properties are illustrated. Moreover, a novel D-MOGWO is developed to deal with the problem because of its NP-hard nature. Experiments are designed and conducted to reveal the effectiveness of these methods on improving Pareto frontiers. Three conclusions can be inferred from the computational experiments: (1) the above-mentioned methods significantly help to improve Pareto frontiers, (2) better Pareto frontiers are achieved by the synergistic use of these methods, and (3) the order of using these methods in optimization procedures matters the optimization results.

Future works would include considering other scheduling problems, including job shop scheduling problems, hybrid flow shop scheduling problems, flexible job shop scheduling problems, and so on. Besides, more environmental friendly scheduling methods can be further investigated.