Minimising the machining energy consumption of a machine tool by sequencing the features of a part

non-cutting on feature sequences. to understand and characterise the MTE while machining a part. A model developed the knowledge gap, and two sub-models for speci ﬁ c energy consumption and actual cutting volume are developed. a single objective optimisation problem, minimising the MTE, is introduced. Two optimisation approaches, Depth-First Search (DFS) and Genetic Algorithm (GA), are employed to generate the optimal processing sequence. A case study is conducted, where ﬁ ve parts with 11 15 features are processed on a machining centre. By comparing the experiment results of the two algo- rithms, GA is recommended for the MTE model. The accuracy of our model achieved 96.25%. 14.13% and can be saved using DFS and GA, respectively. Moreover, the case study demonstrated a 20.69% machining time reduction.


Introduction
An enormous amount of energy (549 quadrillion Btu in 2012) is consumed annually worldwide with an estimated increase of 1.4% per year, and the global energy-related carbon dioxide (CO 2 ) emissions are expected to rise from 32 billion metric tons in 2012 to 36 billion metric tons in 2020 [1]. A large proportion (approximately 25%) is attributable to manufacturing [2,3], and reducing the manufacturing energy consumption and CO 2 emissions becomes significant for alleviating the energy crisis and environmental pollution [4e6]. Machine tools are the basic devices in manufacturing that consume considerable amounts of energy [3,7,8]. According to statistics from the U.S energy information administration, the electricity consumption of machine tools has accounted for above 50% of the total manufacturing electricity consumption [3,9]. Thus, reducing the energy consumption of machine tools has attracted a large amount of attention from both academic research and industrial applications [10,11].
The first step toward reducing the energy consumption of machine tools is to understand and characterise their energy consumption [12]. In particular, a considerable amount of energy that is consumed by machine tools is attributable to the energy consumption of a machine tool by completing a feasible processing plan for a specific part (MTE) [8], and the MTE can be divided to two types: the cutting and non-cutting energy consumption (CE and NCE) [4,13]. The NCE is defined as the energy consumption during run-time operations, including the tool path, tool change and change of spindle rotation speed [14]. The energy consumed when a part is actually cut by a machine tool can be defined as the CE [15]. It has been proved that the values of both CE and NCE are affected by the feature processing sequence of a part [8,14,15]. Finding the sequence that results in the smaller value of the NCE has been proved to be an effective energy consumption reduction approach, and modelling work for the NCE has been developed by Hu [14]. However, the potential for this approach to reduce the CE has not been well explored, and the understanding of the CE is not sufficient in existing research. Usually, the CE accounts for above 60% of the total MTE [8,16].
The CE can be modelled by multiplying the specific energy consumption (SEC) with the actual cutting volume (ACV) [15,17]. The feature processing sequence of a part can affect the value of the CE, because the ACV for a feature can vary if any of its preceding features on the processing sequence are replaced by another feature, while the SEC for the features is different [18,19]. In the existing modelling work for the CE, there are some insufficiencies. For example, the energy of the machine tool is consumed for the axial feeding, spindle rotation, coolant spray, and additional load losses while cutting the features, but these energy portions have not been added to the SEC model, which harm the model's accuracy. Moreover, the mathematic relationship between the processing sequence and the ACV of each feature has not been developed in the CE model, and it is a challenging work to dynamically reflect the relationship. Bridging important gaps and insufficiencies to model and optimise the CE within the MTE has motivated this research, and the proposed solutions are the main contributions of this paper.
Based on the above, this study aims at understanding and characterising the CE and at integrating the developed CE model with the existing NCE model to obtain the completed MTE model. To improve the accuracy, two sub-models have been developed to Nomenclature Feature sequencing problem i, k indices for the features in a part and the positions in a feature sequence n number of the actual features in a part F i i-th feature in a part F C a finite set of n features of a part, F C ¼ fF i g n i¼1 F a finite set of n þ 2 features of a part in a machining environment, F ¼ fF i g nþ1 i¼0 , F C 3F F p a specific feature in a part F 0 , F nþ1 virtual features to denote the start position and end position of a tool while machining W width of a groove w, h diameter and depth of a hole H depth of a hole with its interactive volume included S k feature at the k-th position of a sequence S a finite set to indicate all of the positions of the features in a sequence,  feature processing sequence of a part can affect the value of the  MTE, this article investigates a novel management approach to  reduce the MTE by merely adjusting the processing sequence,  without purchasing additional energy-saving and energy-recycling devices. An MTE model based on arbitrary feature sequences has been further developed. The single objective optimisation in this research is to minimise the MTE for processing a part. Depth-First Search (DFS) and Genetic Algorithm (GA) are modified and used as optimisation approaches to search for the optimal feature sequence. Based on case studies, the developed models and optimisation approaches are demonstrated, compared and discussed. In this study, it is assumed that all of the required processing for a part can be finished on a single machine tool.
In the remainder of this paper, the background and motivation are given in the next section. The problem description and the model for the MTE are presented in Section 3. In Section 4, the working procedures of the developed algorithms for solving the aforementioned optimisation problem are described. Case studies are conducted to demonstrate the developed models and optimisation approaches in Section 5. In Section 6, the results are compared and discussed. In addition, the consequence of the MTE minimisation on the machining time is analysed and discussed. Finally, a brief summary and a description of future work are given in Section 7.

Background and motivation
To understand the MTE, Hu [14] developed the mathematic models to depict the NCE, including the energy of the machine tool consumed for the tool path, tool change and change of spindle rotation speed; however the CE has not been considered. Sheng et al. [15] analysed the effects of the processing sequence of the features on the CE, and the CE was modelled by multiplying the SEC with the ACV. Further, the basic CE model was improved by imposing manufacturing constraints [18] and realising automatic identification of the features [19].
In the existing CE model, there are some insufficiencies on modelling the SEC and ACV, which harm the model's accuracy. For example, in the SEC model, the energy of the machine tool consumed for standby operation, axial feeding, spindle rotation, coolant spray and additional load losses has not been included, but these energy portions are required and account for more than 50% of the total energy consumption while cutting [20,21]. To improve the accuracy of the SEC model, much related experiment research has been conducted. For example, the standby power and coolant spray power of machine tools were accurately measured by a standardised start-up procedure [22,23]. The power models for spindle rotation and axial feeding were expressed as a piecewise linear function and a quadratic function, respectively, and the coefficients in these models were derived by linear regression [24,25]. Considering the additional load losses, the relationship between the material removal power and the process parameters was described by a generic exponential model, and coefficients of seven CNC machine tools were derived by multiple linear regressions [21]. These studies can be used as references for improving the SEC model. In addition, the research on modelling the ACV is still scarce in previous CE models, and it is necessary to develop a mathematical model to reflect the effects of the processing sequence on the ACV of each feature. Based on the previous research, there is still no accurate model for the CE within the MTE.
After developing the MTE model, optimisation approaches can be employed to reduce the MTE. It has been discovered that the processing sequence of the features can affect the machining time [26], quality [27] and energy consumption [15] while machining a specific part. Based on this discovery, many research studies have been conducted to reduce the machining time and improve the machining quality by adjusting the processing sequences of the features [27e30]. However, the application of the feature sequencing approach for reducing the MTE is still scarce, even though the manufacturing energy reduction gains more and more importance in modern manufacturing. Optimisation approaches such as deterministic algorithms and meta-heuristics have been employed to search for the optimal processing sequence of the features, which results in the minimisation of the machining time. These studies can be used as references for minimising the MTE.
Deterministic algorithms such as dynamic programming and branch-and-bound were employed to accurately find the global optimal solution [28,31,32]. However, they are normally suitable for solving small-to-medium sized problems with less than 20 features, due to the sharply increased computation time for larger problems [28]. For example, when the number of features in a part increased from 13 to 16, the computation time of a deterministic algorithm dramatically increased from 29.95 s to 1464.7 s [28]. In comparison, the computation time of a meta-heuristic increased from 0.81 s to only 1.66 s [33]. Thus, meta-heuristics have become increasingly popular because they consume much shorter computation time for large problems. Bhaskara Reddy [26] proved that Genetic Algorithm (GA) can effectively solve the time focused feature sequencing problems. In our article, GA is selected as an optimisation approach. However, finding the global optimal solution is not guaranteed when using meta-heuristics. The global optimal solution is required to be used as the benchmark for comparing the solution quality and computation time. Thus, Depth-First Search (DFS), as a deterministic algorithm, is used to deliver the global optimal solution.
According to the literature reviewed, the understanding and modelling for the MTE, in particular the CE, is still not sufficient. The existing CE model is not accurate. Accordingly, this presented paper contributes to improving the CE model by developing two accurate sub-models (SEC and ACV model). Specifically, the SEC model is improved by considering the energy of the machine tool consumed for the standby operation, axial feeding, spindle rotation, coolant spray and additional load losses. The ACV model is improved by considering the effects of the processing sequence on the ACV of each feature. Then, the developed CE model is integrated with the existing NCE model to obtain the completed MTE model. Based on the MTE model, this article investigates a novel and economical management approach to reduce the MTE by merely adjusting the processing sequence of the features in a specific part. Optimisation approaches based on DFS and GA are first adopted and compared to search for the optimal processing sequence of the features that results in the minimisation of the MTE. The proposed solutions for modelling and optimising the CE and MTE are the main contributions of this paper, and they are introduced in the following sections.

Problem statement and modelling
Considering a part that has n features, a finite set F C ¼ fF i g n i¼1 can be employed to denote these features. Because the start and end positions of the tool affect the value of the MTE, they are defined as two virtual features for the part, denoted by F 0 and F nþ1 , respectively. Hence, in the machining environment, there are n þ 2 features for an n-feature part, which are denoted as a finite set The F C is the subset of the F ðF C 3FÞ. In Fig. 1, a part with two features is used as an example to show that different processing sequences of the features can result in different values of the MTE. The two features (a groove and a hole) of this part are denoted as F 1 and F 2 , and they are milled and drilled, respectively. The start and end positions of the tool, which are virtual features, are denoted as F 0 and F 3 . Two sequences of features can be used to process this part: F 0 À F 1 À F 2 À F 3 and F 0 À F 2 À F 1 À F 3 . The cutting volume of F 1 and F 2 are labelled by solid blue lines and dashed red lines, respectively, in Fig. 1. Some of the parameters about the cutting volume of these two features are given in Fig. 2. The areas filled with purple net are the interactive volume between the F 1 and F 2 . The tool paths of the two sequences are labelled by solid black lines and dashed black lines, respectively, in Fig. 1. The power profiles of a machine tool when processing the part according to the aforementioned two sequences are shown in Fig. 3. The power profiles are developed based on the measured data and the prediction method by Jia [16] and Dahmus and Gutowski [20].
In Fig. 3, the areas filled with horizontal blue lines and vertical red lines represent the CE for milling F 1 and drilling F 2 , and the blank areas represent the NCE for the non-cutting operations such as the tool path and tool change. The total MTE to complete a feasible processing plan for a specific part is considered to be the sum of the CE and the NCE. The NCE for these two sequences are different, because both of their tool paths and tool change times are different according to Fig. 1, and the effect of the feature processing sequences on the NCE can be found in Ref. [14]. This research focuses on the effect of the sequences on the CE. Following the above example, the CE for the sequence F 0 À F 1 À F 2 À F 3 is calculated by summing the CE for milling the F 1 and drilling the F 2 [18], as follows: where SEC 1 and SEC 2 are the specific energy consumption [J/cm 3 ] for F 1 and F 2 , and v 1 and v 2 are the actual cutting volume [cm 3 ] for F 1 and F 2 under the sequence F 0 À F 1 À F 2 À F 3 , as shown in Fig. 2(c). When the sequence is changed to F 0 À F 2 À F 1 À F 3 , the milling volume for F 1 will decrease while the drilling volume for F 2 will increase. As a result, the CE for milling F 1 decreases while the CE for drilling F 2 increases, as reflected by the size change of the areas filled with horizontal and vertical lines between Fig. 3(a) and (b). More specifically, according to the dimensional parameters in Figs. 1 and 2(a), the milling volume for F 1 decreases ðH À hÞ Â p Â w 2 =4 while the drilling volume for F 2 increases ðH À hÞ Â p Â w 2 =4, and the changed cutting volume for both F 1 and F 2 is equal to the interactive volume between them, as shown in Fig. 2(b) and (d). Thus, the CE for the sequence F 0 -F 2 -F 1 -F 3 becomes the following: According to Equations (1) and (2), SEC 1 and SEC 2 should be calculated to compare the CE for these two sequences. For example, if SEC 1 < SEC 2 , the CE for the sequence F 0 -F 1 -F 2 -F 3 will be less than that for the F 0 - Following the example, the aim of this research is to determine the optimal feature processing sequence for a part that minimises the total MTE. Because there are n þ 2 features for an n-feature part in the machining environment, a finite set S ¼ fS k g nþ2 k¼1 is employed to indicate all of the positions of the features in a sequence. S k indicates the feature at the k-th position of a sequence. For example, S k ¼ F p indicates the feature at the k-th position of a sequence is the feature F p . For any part, the feature at the 1-st position and n þ 2-th position is F 0 (S 1 ¼ F 0 ) and F nþ1 (S nþ2 ¼ F nþ1 ), respectively. Then, the objective function for minimising the MTE based on a specific feature sequence can be expressed as follows: where E s is the total MTE based on a specific feature sequence; is the NCE between the feature at the k-th position and the feature at the k þ 1-th position of the sequence, and the detailed model for E ðS k ;S kþ1 Þ non can be found in Ref. [14]; E S k cut is the CE for the feature at the k-th position of the sequence. It is given that the feature at the k-th position is F i (S k ¼F i ). Then, the E S k cut is expressed as: where E Fi cut is the CE for the feature F i ; SEC i and v i are the specific energy consumption [J/cm 3 ] and the actual cutting volume [cm 3 ] for the feature F i , respectively, which are modelled in Section 3.1 and 3.2.
The constraints of the model are developed according to the precedence constraints among the features. A feasible feature sequence, namely, a feasible solution for the mathematic model, must satisfy all of the constraints. The total MTE for the corresponding feature sequence is set to infinity "∞" once any feature and its pre or post features in a sequence violate any constraint.

Modelling the specific energy consumption (SEC) for the feature
The SEC is defined as the energy consumption of the machine tool for removing 1-cm 3 material [17,34]. The SEC for the feature F i can be modelled [17,35] as: where D i and f i are the drilling diameter [mm] and feed rate [mm/r] for drilling the feature F i . In Equation (5), P i can be divided to five portions: the material removal power, axial feeding power, spindle rotation power, coolant spray power and standby power [16,25]. The standby power and coolant spray power of a machine tool are constant and remain the same while cutting different features, as labelled with P 0 and P CS , respectively, in Fig. 3. The material removal power, axial feeding power and spindle rotation power of a machine tool are variable while cutting different features. Thus, P i is expressed as follows: where P MCi , P AFi and P SRi are the material removal power, axial feeding power and spindle rotation power, respectively, for the feature F i ; and P CS and P 0 are the coolant spray power and standby power of the machine tool, which are obtained by actual measurement.
In Expression (8), P MCi is modelled according to the cutting method. For example, if the cutting method is milling, then a model in Ref. [36], which considers the additional load losses for milling and has high reliability, can be employed and is expressed as: where C M is the coefficient in the material milling power model, and v si is the cutting speed [m/min] for the feature F i ; w M , y M , x M and u M are the exponent of the cutting speed, feed rate, milling depth and milling width, respectively, in the material milling power model. The coefficient and exponents are obtained by statistical analysis based on experiment data. If the cutting method is drilling, then a model in Ref. [16], which considers the additional load losses for drilling and has high reliability, can be employed and is expressed as: where C D is the coefficient in the material drilling power model; z D and y D are the exponent of the drilling diameter and feed rate, respectively, in the material drilling power model. In Expression (8), P AFi is modelled according to the feeding speed, the type of machine tool and the feeding direction. For a 3axial machine tool in the rectilinear feeding direction, P AFi can be modelled by the sum of the feeding power of all of the axes [16], as follows: where P XFi , P YFi and P ZFi are the feeding power of the X-axis, Y-axis and Z-axis, respectively, while cutting the feature F i . By comparing several models for the axial feeding power [25], a quadratic model in Ref. [36], which has the highest accuracy, is employed. The models for P XFi , P YFi and P ZFi are similar. For example, following [36], P XFi can be modelled as: where A XF and B XF are the quadratic and monomial coefficients in the model, respectively, which can be obtained by quadratic regression based on experiment data; v fi is the feeding speed [mm/ min] for the feature F i ; and a i is the angle between the feeding direction of the tool and the X-axis while cutting the feature F i .
In Expression (8), a linear equation in Ref. [36], which has a high accuracy [25], is employed and modified to model the P SRi , as follows: where B ri SR and C ri SR are the monomial coefficient and constant in the model at the rotation speed of r i , respectively, which are obtained by linear regression based on the experiment data.

Modelling the actual cutting volume (ACV) for the feature
The ACV for the feature F i is a variable that involves the change of its position in the sequence caused by the volumetric interaction, as analysed in the aforementioned two-feature part. Even if the feature F i is located at a specific position, its ACV varies with the change of its previous features in the sequence, and this characteristic greatly increases the computation time to obtain the ACV for each feature in an arbitrary sequence. Considering the impact of the sequence of a feature on its cutting volume due to the volumetric interaction with the previous features, the inclusionexclusion principle (IEP), which is a famous and effective combinatorial method and counting technique for determining the cardinality of a set [37e39], is employed to model the ACV for the feature F i at the k-th position of the sequence, as shown in Expression (14). The basic idea of the application of the IEP is that the cutting volume for the feature F i with all of its interactive volume included is counted first, and then, the interactive volume with its previous features is excluded to obtain the ACV.
In this expression, v i is the ACV for the feature F i ; V i is the cutting volume for the feature F i with its interactive volume included, which is equal to the cutting volume for the feature F i from the blank to the objective feature, as shown in Fig. 2(c) and (d). G j is the cutting volume for the feature at the j-th position with its interactive volume included, and ∩ is the operator for calculating the interactive volume between the features. Here, j 1 , j 2 and j ðMÀ1Þ are the indices for the specific positions in the feature sequence, which are integers. The index for the position cannot be 1, because G 1 indicates the cutting volume for the feature at the 1-st position while the corresponding feature is the virtual feature F 0 , which does not have a volume. M is the maximum number of the features that interact with one another simultaneously in the part.

Solution algorithms
To alleviate the computation burden and reduce the computation time in searching for the optimal processing sequence of the features, two solution algorithms, including Depth-First Search (DFS) and Genetic Algorithm (GA), are employed, by using minimisation of the MTE as the objective. DFS is selected as one of the solutions because it can always accurately find the global optimal solution in each run. In comparison, GA normally consumes less time to reach the optimal or near-optimal solution; however, finding out the global optimum is not guaranteed. Hence, experiments are delivered in Section 5 to compare the performance of the two algorithms in solving the aforementioned problem in terms of the solution quality and computation time.

Depth-First Search (DFS)
DFS is a type of enumeration algorithm that uses the tree traversal technique, which visits a path in a search tree as far as possible prior to backtracking and visiting the next path [40]. The computation efficiency of DFS is effectively improved by a pruning operation to avoid visiting an inferior path. A flowchart of DFS can be found in Ref. [41].
At the beginning of DFS, a search tree that can cover all of the possible paths (processing sequences of the features of a part) is generated and any feasible path can be regarded as the initial upper bound. The nodes (features) of the search tree are visited based on the depth-first searching strategy [42]. If the visited node is not qualified, DFS prunes and backtracks to a previous node. If the visited node is qualified and is the end of the path, the upper bound will be updated. Then, DFS prunes and backtracks to a previous node. If the visited node is qualified but is not the end of the path, DFS continues to visit the next node. After pruning and backtracking, DFS will check whether the stopping conditions have been met or not. If the stopping conditions are met, the latest upper bound will be reported, and DFS stops; otherwise, the next node will be visited. The stopping conditions can be that all of the nodes of the search tree have been visited.

Genetic Algorithm (GA)
GA is a robust heuristic stochastic searching method that imitates the natural evolutionary process by combining "survival of the fittest" to generate offspring [43]. The weak solutions in each generation are more likely to be replaced by the offspring, and finally, the optimal or near-optimal solutions are generated to solve the optimisation problem [26]. A flowchart of GA can be found in Ref. [44].
At the beginning of GA, an initial population is randomly generated with a size of N chromosomes as the first generation. For our research, a chromosome represents a sequence of the features, encoded by integer coding [45]. Each gene in the chromosome represents a specific feature. During the generational process, chromosomes are selected to reproduce a new generation through a fitness-based evaluation by a fitness function [46]. The fitness function is defined as the genetic representation and used to evaluate the quality of each chromosome in the population [47]. In our research, the fitness function is the reciprocal of the Expression (3) as: . The chromosomes that have a higher fitness can have a higher probability to be selected through a selection operator: the proportional roulette wheel selection [48]. Then, the next generation of the chromosomes is generated from the selected chromosomes through the GA operators of (1) Crossover: an exchange of sections between chromosomes and (2) Mutation: a random modification on the chromosome [26]. Specifically, the partially mapped crossover (PMX) [49] and the swap mutation [50] are adopted as the crossover operator and the mutation operator, respectively. This generational process is repeated until a stopping condition has been met. The stopping condition can be the specified maximum generation number that is reached.

Case study
A prismatic part named part C, which has 13 actual features, is used in this case study, and part C is made of 45#Steel. In the context of this research, total 15 features are considered for this case, which are F 1 (plain), F 2 (stair), F 3 (groove), F 4 (depression), F 5 (notch), F 6 (notch), F 7 (hole), F 8 (hole), F 9 (hole), F 10 (hole), F 11 (hole), F 12 (hole), F 13 (hole) and 2 virtual features F 0 and F 14 , as shown in Fig. 4. A machining centre (XHF-714F) is employed to process part C, and the cutters for milling and drilling are W400F-FS (4 teeth) and NACHI SD8, respectively. The experiment setup for the power data collection on the XHF-714F is shown in Fig. 5. The standby power and coolant spray power of XHF-714F are shown in Table 1, which are obtained by experiment measurement. The coefficients in the power models of the spindle rotation and axial feeding are listed in Table 2, which are obtained by regression based on the experiment data. Table 3 lists the coefficients and exponents in the material removal power model under specific cutting conditions, which are also obtained by regression. The milling and drilling parameters for part C are listed in Tables 4 and 5, which are obtained from the process files. On such a basis, the MTE for each feature in part C can be calculated.
This case aims at sequencing the features in part C for the minimum MTE. According to Expression (3), the objective function is expressed as: The values of E ðS k ;S kþ1 Þ non are obtained from Table 6. Table 6 shows the NCE between the features in part C, and the detailed calculation method can be found in Ref. [14].
In the objective function, the calculation of E S6 cut is taken as an example. It is given that the feature at the 6-th position is F 2 (S 6 ¼F 2 ). According to Expression (4), the E S6 cut is: where SEC 2 and v 2 are the SEC and ACV for feature F 2 , respectively. The detailed calculation procedures for SEC 2 and v 2 are provided in Section S.1 of the Supplementary Information (SI).
In particular, the value of E S6 cut may change once the feature at the 6-th position of the sequence is not F 2 or the features at the 2-nd position, 3-rd position, 4-th position and 5-th position are not F 1 ,   Based on the data and models above and in Section S.1 of the SI, two solution algorithms, including Depth-First Search (DFS) and Genetic Algorithm (GA), are employed as optimisation approaches. DFS and GA used in this research are developed on Dev Cþþ 5.11.0 software with the programming language Cþþ. The computing platform is Intel (R) Core (TM) i7-2630 QM CPU with 2.00 GHz frequency; 4.00 GB RAM; Windows 7 (64bit) operating system. For part C, the developed DFS always returns the global minimum MTE, which is 1246.10 KJ, after 2399 iterations. A computation time of DFS is 907.7 s, and the global optimal feature sequence is F 0 À F 1 À F 2 À F 9 À F 10 À F 11 À F 12 À F 13 À F 7 À F 8 À F 5 À F 6 À F 3 À F 4 À F 14 . The searching process of DFS is provided in Section S.2 of the SI.
The parameter values used in GA are obtained by fine tuning, and their values are as follows: population size ¼ 100; crossover probability ¼ 0:9; mutation probability ¼ 0:05; and iteration Table 2 Coefficients in the power models of spindle rotation and axial feeding.

Item Coefficient
Monomial coefficient and constant in the spindle rotation power model (B ri SR , C ri SR ) (0 < r i 2200) (0.086, 14.76) Monomial coefficient and constant in the spindle rotation power model (B ri SR , C ri SR ) (2200 < r i 3000) (0.0186, 164.97) Quadratic coefficient in the feeding power model of X-axis, Y-axis, Z-axis upward and Z-axis downward (5 Â 10 À7 , À1 Â 10 À6 , À5 Â 10 À7 , À1 Â 10 À7 ) Monomial coefficient in the feeding power model of X-axis, Y-axis, Z-axis upward and Z-axis downward (0.0491, 0.043, 0.059, 0.0461) Table 3 Coefficients and exponents in the material removal power models.    Table 6 Non-cutting energy consumption (NCE) between the features in part C. ¼ 300. By running the developed GA 50 times, it sometimes returns the global minimum MTE (1246.10 KJ) after 300 iterations. A computation time of GA is 2.468 s, and the corresponding feature sequence is the same as the optimum produced by DFS. The searching process of GA for the global optimum is provided in Section S.2 of the SI. For this specific optimisation problem, GA usually stops converging within 75 iterations. However, in most of the performed runs, GA can only return the near-optimal solutions. For example, a near-minimum MTE that GA can get is 1248.69 KJ, and the corresponding feature sequence is F 0 À F 1 À F 2 À F 9 À F 4 ÀF 3 À F 6 À F 5 À F 7 À F 8 À F 11 À F 12 À F 13 À F 10 À F 14 . It also consumes less computation time (2.355 s) to get the near-optimal solution. The searching process of GA for the near-optimum is provided in Section S.2 of the SI. The comparisons of DFS and GA for the 13-feature part C in these 50 trials are summarised in Table 7.

Item
The developed models and optimisation approaches are also tested and validated on four parts with 11, 12, 14 and 15 actual features, respectively, as shown in Fig. 6. The results from using DFS and GA for optimising the feature sequence of these parts are compared and summarised in Table 7. According to the optimisation results, GA has a 30% or less possibility of finding the global optimal solution, and it frequently returns a near-optimal solution due to the nature of meta-heuristics. The deterministic algorithm, DFS, can always find the global optimal solutions in each trial, so DFS performs better than GA in terms of the solution quality for all parts. The solution quality of GA is only slightly inferior to the solution quality of DFS, but its computation time is much less than that of DFS. For example, the median MTE that was obtained by GA for part C is only 0.149% [(1247.95e1246.10)/1246.10] inferior to the solution obtained by DFS, but its computation time is approximately 99.73% [(907.7e2.426)/907.7] less than that of DFS. Moreover, when the number of features in a part is increased, the superiority of GA in terms of the computation time becomes more prominent. For example, when the number of features in the part is increased from 11 to 15, the computation time of DFS increases from 5.870 s to an intolerable 173700 s. In comparison, the computation time of GA increases from 1.807 s to only 3.182 s. In summary, GA is recommended as the better optimisation approach for our problem because GA requires much less computation time with little sacrifice in the solution quality.

Results and discussion
In this section, the developed CE model and optimisation approaches are compared with the results presented in the literature, and their superiorities are discussed and validated. In addition, the consequence of the MTE minimisation on the machining time is analysed and discussed.

Comparison and verification of the model accuracy
To compare and verify the accuracy of the developed CE model, .
Then, the experiment setup in Fig. 5 is used to measure and record the actual power and energy consumption of XHF-714F while drilling the F 1 . The measured information is shown in Figs. 7 and 8.
By analysing the fluctuations in the power values that are recorded in the database, the time interval for drilling the F 1 is (27.4s, 41.4s), as shown in Fig. 7. By querying the energy consumption data while drilling the F 1 , the values of the energy consumption are 16569.1J at the beginning time (27.4s) and 28808.7J at the end time (41.4s), as shown in Fig. 8. Thus, the measured energy consumption for drilling the F 1 is: E F1 cut ðMeasuredÞ ¼ 28808:7 À 16569:1 ¼ 12239:6J: Then, the accuracy of our CE model for the F 1 can be calculated by making a comparison between the estimated and measured values, as follows: ADM ¼ 1 À j12239:6 À 11780:3j=12239:6 ¼ 96:25%: Similarly, the accuracy of the CE model developed by Yin et al. [19] for the F 1 can be calculated as: AEM ¼ 1 À j12239:6 À 2573:7j=12239:6 ¼ 21:03%: Thus, the accuracy of our developed CE model (96.25%) is higher than that of the existing model (21.03%). However, it is more timeconsuming and labour-intensive to develop our model, because it is required in our model to collect and process the energy  consumption data of the machine tool such as axial feeding, spindle rotation, coolant spray and additional load losses. Although the CE model by Yin et al. [19] is not accurate, it is more effortless and easier to be developed without any experiments. While selecting a CE model for an industrial application, a trade-off between the model accuracy and simplification should be made.

Comparison of sequencing approaches
The proposed optimisation approaches are compared with the existing approaches for reducing the MTE, and the traditional technique Bottom-to-Top (BTT) [51] is selected as an existing approach. For example, the processing sequence of the features for part C, which is generated by the BTT, is F 0 À F 10 À F 8 À F 13 À F 5 À F 3 À F 9 À F 4 À F 11 À F 7 À F 12 À F 2 À F 6 À F 1 À F 14 . The MTE for this sequence is 1451.12 KJ. In comparison, the MTE for the sequence that is generated by our proposed optimisation approach DFS is 1246.10 KJ. Therefore, 14.13% [(1451.12e1246.10)/1451.12] MTE can be saved by using the developed DFS compared to the BTT technique. By comparing the processing sequence generated by DFS with that of BTT, the sequence by DFS has a shorter tool path and fewer tool change times, and it tends to process the features with lower SEC in priority. These contribute to saving the MTE. Further, the processing sequences generated by GA and BTT are also compared. Considering the near-optimal solutions that are generated by GA, 14.00% [(1451.12e1247.95)/1451.12] MTE can be saved by using the developed GA. In summary, both DFS and GA are effective optimisation approaches for reducing the MTE.

Machining time reductions benefit from MTE minimisation
In a real manufacturing circumstance, it is unrealistic to only reduce the MTE without controlling the machining time, which would cause a machine tool tardiness problem. Thus, the consequence of the MTE minimisation on the machining time is analysed and discussed. By referring to Expression (3), the machining time based on a specific feature sequence can be expressed as: is the non-cutting time between the feature at the k-th position and the feature at the k þ 1-th position of the sequence; and T S k cut is the cutting time for the feature at the k-th position of the sequence. Then, the machining time based on the optimal processing sequence of part C, which results in the MTE minimisation, is calculated as follows: In comparison, a processing sequence of part C without the energy-saving consideration is generated by the BTT technique, and the machining time for this sequence is calculated as: The calculation procedure for T s ðMTE minimisationÞ and T s ðBTT techniqueÞ are similar, and the measured time data and calculation for T s ðMTE minimisationÞ are taken as example in Section S.4 of the SI. In this case, 20.69% [(1499.1e1188.9)/1499.1] of the machining time reductions benefit from the MTE minimisation. However, it is not guaranteed that the MTE minimisation can always lead to the machining time reduction, because the machining time is not strictly proportional to the MTE [21]. For example, when a feature sequence with the minimum machining time is changed to a sequence with the minimum MTE, its machining time can increase. Thus, when optimising the processing sequence of the features for a specific part, the trade-off between the minimum machining time and MTE should be researched in the future.

Conclusions and future work
It has been approved that both the cutting and non-cutting energy consumption of the machine tool (CE and NCE) can be reduced by sequencing the features of a part to be processed at the process planning stage. However, sequencing features to reduce the CE has not been well explored in previous research, and the CE accounts for more than 60% of the total machining energy consumption of the machine tool (MTE). In particular, the energy of the machine tool consumed for the axial feeding, spindle rotation, coolant spray and additional load losses while cutting features has not been added to the specific energy consumption (SEC) model, which harms the CE model's accuracy. Moreover, the mathematic relationship between the processing sequence and the actual cutting volume (ACV) of each feature has not been developed in the CE model. To bridge important gaps and address insufficiencies in the CE model, more accurate sub-models for the SEC and ACV are developed. Based on the above, a mathematic model is developed for the single objective optimisation problem that minimises the MTE. Then, two algorithms, including Depth-First Search (DFS) and Genetic Algorithm (GA), are employed and compared as optimisation approaches to minimise the MTE. Because DFS can always find the global optimal solution accurately, the results obtained by GA can be validated.
In the case study, the optimal and near-optimal sequences for five parts with 11e15 actual features have been found. According to the optimisation results, GA usually returns a near-optimal solution. In comparison, the deterministic algorithm, DFS, can always find the global optimal solutions. Thus, DFS performs better than GA in terms of the solution quality for all of the parts. In these cases, the solution quality of GA is only slightly inferior to that of DFS, but its computation time is much less than that of DFS. For example, the median MTE obtained by GA for part C is only 0.149% inferior to the solution obtained by DFS, but its computation time is approximately 99.73% less than that of DFS. Moreover, when the number of features in a part is increased to 15, the superiority of GA in terms of the computation time becomes more prominent. Thus, GA is recommended between DFS and GA based on the experiment results.
A machining experiment is conducted to validate the accuracy of the developed CE model. The results show that the accuracy of our model can achieve 96.25%, which is higher than that of an existing model. To compare the effectiveness of the optimisation approaches, the proposed DFS and GA are compared with a traditional technique, Bottom-to-Top (BTT), in reducing the MTE. For part C, 14.13% and 14.00% MTE can be saved by using DFS and GA, respectively, compared to the BTT. Moreover, the consequence of the MTE minimisation on the machining time is discussed. In the case study, 20.69% of the machining time reductions benefit from the MTE minimisation. However, it is not guaranteed that the MTE minimisation can always result in the machining time reduction.
In the presented research, it is very labour-intensive to acquire and process the time and power data of the machine tool for modelling the MTE. Thus, the automation for the data acquisition and process can be improved. For the optimisation approach, DFS consumes too much time for the medium-to-large sequencing problem, although it can always find the global optimal solution. Comparatively, GA normally consumes less time, however, it usually returns a near-optimal solution. Thus, an algorithm that has a good balance between the computation speed and solution quality can be researched. The single objective is a limitation. Other optimisation objectives, including machining time, quality and cost, should be considered while optimising the MTE. In the future, research on combining the proposed energy-aware feature sequencing method with the product design software, such as SolidWorks, Pro/E, UG and CATIA, will be conducted to realise its industrial application better.