A Hybrid Grey Wolf Optimizer for Process Planning Optimization with Precedence Constraints

Process planning optimization is a well-known NP-hard combinatorial problem extensively studied in the scientific community. Its main components include operation sequencing, selection of manufacturing resources and determination of appropriate setup plans. These problems require metaheuristic-based approaches in order to be effectively and efficiently solved. Therefore, to optimize the complex process planning problem, a novel hybrid grey wolf optimizer (HGWO) is proposed. The traditional grey wolf optimizer (GWO) is improved by employing genetic strategies such as selection, crossover and mutation which enhance global search abilities and convergence of the traditional GWO. Precedence relationships among machining operations are taken into account and precedence constraints are modeled using operation precedence graphs and adjacency matrices. Constraint handling heuristic procedure is adopted to move infeasible solutions to a feasible domain. Minimization of the total weighted machining cost of a process plan is adopted as the objective and three experimental studies that consider three different prismatic parts are conducted. Comparative analysis of the obtained cost values, as well as the convergence analysis, are performed and the HGWO approach demonstrated effectiveness and flexibility in finding optimal and near-optimal process plans. On the other side, comparative analysis of computational times and execution times of certain MATLAB functions showed that the HGWO have good time efficiency but limited since it requires more time compared to considered hybrid and traditional algorithms. Potential directions to improving efficiency and performances of the proposed approach are given in conclusions.


Introduction
Process planning optimization (PPO) problems consist of two subproblems: operations selection and operations sequencing [1]. Operations selection is a task of selecting necessary machining operations for each machining feature recognized on a part or a product. In order to perform these operations, adequate machining resources are required. In that sense, the PPO problem includes selection of machine candidate, cutting tool candidate and tool approach direction (TAD) candidate for each machining operation. On the other hand, operations sequencing represents a task of finding the order of selected machining operations with respect to the predetermined precedence constraints based on relationships between machining features. Tackling with such a problem requires dealing with a number of alternatives which makes the PPO a combinatorial optimization problem. A number of alternative process plans grows in accordance with problem dimensions. In other words, as a number of machining features increases, so does a number of machining operations. This leads to exponential time growth required to find optimal or near-optimal process plans, which means that the PPO problem belongs to the class of NP-hard (non-deterministic polynomial) optimization problems.
As it is impossible to check every possible solution to an NP-hard problem, conventional non-heuristic methods have proved to be ineffective. The PPO problem requires the implementation of metaheuristic algorithms which have shown promising performance for solving difficult problems. This is especially due to their ability to achieve a trade-off between local and global search [2]. This paper proposes a novel hybrid grey wolf optimizer (HGWO) to deal with the PPO problem with precedence constraints. The manuscript is organized as follows. Section 2 provides recent studies related to the PPO. Section 3 describes the problem, focusing on manufacturing interactions and precedence constraints, representation of process plans and mathematical modelling of the PPO. Section 4 emphasizes the proposed HGWO approach starting from traditional GWO, then genetic strategies that are included, and constraint handling algorithm as part of the structure of the HGWO. Section 5 represents the experimental research that consider three experimental studies, or three prismatic parts. With detailed description of these cases, comparative analysis of cost values and convergence curves are performed to test the performances of the HGWO. Later, discussion in Section 6, emphasize the effectiveness of the proposed method and the computational efficiency. Conclusions with directions for future research are given in Section 7.

Related Research Studies
In the past few decades, many search algorithms have been employed to deal with the PPO problem. In the initial stages of research on the PPO, genetic algorithms have been widely reported among other metaheuristics. One of the earliest implementations of genetic algorithms (GAs) can be found in Zhang et al. [3]. Besides the GA, many other algorithms such as ACO, PSO, TS, SA are applied to this day.
GA it is the most popular evolutionary algorithm based on a Darwinian theory of natural selection and genetics. Traditional GAs, as well as many other algorithms, require additional modifications to improve their convergence as well as global and local search abilities. Study performed in Salehi and Bahreininejad [4] adopts hybrid GA and intelligent search for optimization of process planning in preliminary and detailed planning stages. The initial and feasible solutions were checked by order and clustering constraints. Li et al. [5] proposed a hybridization model of GA and SA where the initial process plans were generated using the GA approach while the optimal or near-optimal process plans were obtained using the SA algorithm. Kafashi [6] developed a GA approach in order to tackle integrated setup planning and operation sequencing problem in flexible manufacturing environment. With the emphasis on technological constraints, tolerance relation analysis and feature-based representation, satisfactory setup plans and operation sequences were generated. Cai et al. [7] used the GA approach based on optimization toolbox to address adaptive setup planning problem towards the process planning and scheduling integration. They considered machine availability based on tool accessibility analysis as a constraint, and selected certain scheduling criteria such as cost, makespan or machine utilization. Huang et al. [8] proposed a hybrid graph and GA approach. Combining graph theory as well as matrix theory, precisely operation precedence graphs and adjacency matrices, precedence constraints were formulated and a GA approach based on modified crossover and mutation strategies were utilized to solve the PPO problem. Two heuristic approaches are incorporated within this hybrid approach in order to adjust infeasible process plans to a feasible domain. Li et al. [9] employed a novel two phase genetic algorithm to optimize process parameters and machining sequence for two-tool parallel drilling operations with multiple blind holes distributed in a pair of parallel faces and in multiple pairs of parallel faces. Hybridization of GA and local search for the PPO on turning machine tool was proposed by Su et al. [10]. They formulated the PPO problem as a mixed 0-1 integer programming model and used a novel encoding strategy to deal with complicated precedence constraints. In the study presented by Su et al. [11], precedence constrained operation sequencing problem was formulated and edge selection based strategy was adopted to assure feasibility of solutions and improve GA's converging efficiency. Dou et al. [12] proposed the improved GA by considering fragment crossover and mutation operators with adaptive operation probabilities as well as a new elitist-based crossover strategy. Luo et al. [13] dealt with a large-sized process sequencing problem with complex association constraints. The problem with huge and complex solution space was decomposed to smaller multi-neighborhood spaces and hybrid GA and variable neighborhood search (VNS) approach were used to obtain the best solutions through all neighborhood spaces.
Furthermore, the ant colony optimization algorithm also found its application to the PPO. Liu et al. [14] used ant colony optimization for process planning optimization. They based their research on mapping the PPO problem to a weighted graph and converting to a constraint-based travelling salesman problem. Constraint and state matrices were used to ensure the feasibility of process plans. Wang et al. [15] dealt with the PPO problem by using a weighted directed graph for process plan representation. Wang et al. [16] proposed a two-stage ACO algorithm to solve the PPO problem for prismatic parts. The PPO was formulated using a directed graph and then the proposed ACO approach was used for optimizing process plans for two prismatic parts from the literature. First stage included operation selection while the second was focused on operation sequencing. Hu et al. [17] considered precedence constraint and clustering constraint relationships to ensure feasible permutations and enhanced search abilities. They used the updating method and the local search mechanism to modify the ACO algorithm.
Particle swarm optimization (PSO) is also one of the most popular metaheuristics based on swarm intelligence. This method is inspired by social behavior of living organisms that live in groups (swarms), such as birds or fish [18]. One of the recent studies concerning the PSO in process planning can be found in [19]. Here, authors proposed a modified PSO algorithm improved by adopting efficient encoding, updating and random search methods in order to tackle seven different case studies that consider prismatic parts. Petrović et al. [20] developed a new chaotic PSO algorithm for flexible process planning. By emphasizing different types of flexibilities, authors used AND/OR networks to represent process plans and then tested performances of the cPSO on several cases including cylindrical as well as prismatic parts. Miljković and Petrović [21] developed a modified multi-objective PSO algorithm for flexible process planning. By emphasizing different types of flexibilities, authors used AND/OR networks to represent feasible process plans, and tested performances of the cPSO on several cases including cylindrical as well as prismatic parts. Dou et al. [22] proposed a novel feasible sequence-oriented discrete PSO algorithm to solve the operation sequencing problems in CAPP. It incorporates novel crossover and mutation operators with adaptive probabilities to evolve particles and improve exploration ability.
Lian et al. [23] proposed the imperialist competitive algorithm to address the PPO problem. This socio-politically motivated population-based metaheuristic was inspired by imperialist competition. The approach utilizes steps of assimilation, competition, revolution and elimination to obtain optimal process plans according to predefined networks and types of flexibilities.
Wen et al. [24] proposed a new method based on honey bees mating optimization (HBMO) algorithm to optimize the PPO. The solution encoding, crossover operator and local search strategies were developed, and three experiments were carried out to demonstrate improvement of the HBMO approach.
A new heuristic method, called cross entropy approach was proposed to optimize flexible process planning by Lv and Qiao [25]. Its authors adopted AND/OR networks to represent flexible process plans and established mathematical model for the minimization of total processing time and total cost of flexible process plans. The new sample representation and probability distribution parameter were introduced and case studies were carried and discussed to indicate the performance and adaptability of the CE approach.
Wang et al. [26] adopted a hybrid bat algorithm for the PPO focusing on both crucial tasks, operations selection and operations sequencing. Encoding, decoding and initialization strategies were proposed and two local search strategies were incorporated into the standard bat algorithm (BA) to improve its local convergence. A classical simulation experiment was conducted to verify the validity of the hybrid BA.
Musharavati and Hamouda [27] investigated the possibility of implementing four different configurations of simulated annealing algorithm to solve the process planning problem in reconfigurable manufacturing systems. They used knowledge exploitation and quasi-parallelism as concepts to enhance the SA algorithms. Performances of the variant algorithms are compared and ANOVA methodology is used in computational analysis to test the means and indicate improvements towards better optimal solutions.
Mohammadi et al. [28] dealt with the multi-objective optimization of the integrated process planning and scheduling problem. They designed a slot-based mixed integer linear programming model that accounts for sequence-dependent preparation times. Minimization of manufacturing cost in process planning, and minimization of preparation times and total tardiness were considered as optimization objectives. To solve the problem, a hybrid simulated annealing approach was introduced. The SA was modified with tabu search algorithm as well as local search strategies in order to improve solution's quality and avoid premature convergence.
Xu et al. [29] introduced a novel NC process reuse-oriented flexible process planning optimization approach for prismatic parts with the objective of minimizing total manufacturing cost. A hybrid ant colony algorithm (ACA) and simulated annealing (SA) approach based on operation precedence graph (OPG) was presented to find the global optimal NC process scheme for the part.
Lian et al. [30] proposed a multi-dimensional tabu search algorithm to optimize four dimensions of a process plan, such as operation sequence, machine sequence, tool sequence, and tool approach direction sequence. Classical neighborhood strategies such as insertion and swap are used to improve the performance of MDTS.
Falih and Shammari [31] proposed a novel hybrid constrained permutation algorithm and genetic algorithm approach for process planning problem. A set of feasible operation sequences is generated using a CPA algorithm and the GA with mixed crossover operator is further employed to search for optimal or near optimal process plans.
Gao et al. [32] suggested the intelligent water drop algorithm for solving the process planning problem. Firstly, operation units were defined according to the processing characteristics, and later, the IWD algorithm was combined with the process planning problem.
Kizys et al. [33] developed an iterative local search metaheuristic and quadratic programming approach to deal with variants of the mean-variance portfolio optimization problem subject to cardinality and quantity constraints.
Sawik and Sawik [34] used stochastic programming approach to optimize cybersecurity investment in supply chains. A mixed binary optimization problem is transformed to an unconstrained binary program in order to maximize total cybersecurity value of control portfolio. It is also shown that the portfolio of security controls with maximum total cybersecurity value reduces the losses from security breaches and mitigates the impact of cyber risk.
In one of the recent studies, Milosević et al. [35] presented the intelligent process planning for the concepts of smart factory and smart manufacturing. Nature-inspired metaheuristic algorithms as smart services of artificial intelligence for intelligent process planning optimization within smart manufacturing were proposed. Three modern algorithms such as grey wolf optimizer (GWO), whale optimization algorithm (WOA) and crow search algorithm (CSA) were employed on three classical case studies from the literature.
Djurdjev et al. [36] proposed a novel genetic crow search approach (GCSA) for operation sequencing in process planning. The traditional CSA algorithm has been improved by adopting genetic strategies such as tournament selection, three-string crossover, shift mutation and resource mutation. Adaptive probabilities were introduced to improve local and global capabilities of the GCSA. Besides, the nearest mechanism strategy was added to ensure feasibility of machine, tools and TADs vectors. Lastly, repairment heuristic strategy was used to handle precedence constraints among features and machining operations.
Aforementioned studies have shown that the improvements made by authors largely increased efficiency of metaheuristic algorithms. However, as stated in [20], main drawbacks reflect in time-consuming optimization and slow convergence of TS, SA and GA approaches particularly. The key and also the most challenging task in the implementation of metaheuristic algorithms is to find a proper balance between exploration and exploitation phases [37], and thereby, maintain its simplicity and flexibility. Many different concepts have been developed and tested so far, making them available for application in different scientific fields. In this research, we focus on implementing novel metaheuristic approach to solve the PPO problem.

Process Planning Optimization Problem
Computer-aided process planning represents the link between computer-aided design and computer-aided manufacturing within the computer integrated environment. The aim of process planning is to manufacture a part or a product starting from its initial design stages to its final stage (a finished part or product). The input to process planning includes valuable geometric information obtained from a CAD file. Accordingly, process planners have a goal to map product information to a technological domain. Parts/products are generally described by design features which are characterized by geometric forms or shapes with technological attributes such as tolerances and surface finishes. Design features are transformed to manufacturing features which represent manufacturing meaning of the geometry of a part, a product or an assembly associated with certain manufacturing activities [38]. An important subset of manufacturing features consists of machining features such as holes, steps, grooves, pockets, planes, etc. Recognized machining features are most often considered as an input information for the PPO problem.
Apart from machining features that require certain operations in order to be machined on a part or a product, process planning considers other activities such as finding the sequence of machining operations, determining machining resources-machine tools, cutting tools, fixtures, determining cutting conditions and calculating machining time and cost. In the PPO, the quality of a single process plan is determined by selecting machining resources such as machine tools, cutting tools, and tool approach directions (TADs) for each machining operation, as well as generating optimal sequence of machining operations for an observed part or a product. During the decision-making process of selecting alternatives of machining resources and TADs, in order to find feasible and optimal order of machining operations, the sequence have to satisfy precedence constraints that are based on relationships between machining features and machining operations.

Representation of Process Plans
When dealing with the PPO problem, the first step is to select an appropriate representation of a process plan. In recent studies [20,21,23,25] five types of flexibilities were considered, and AND/OR networks are used to represent the PPO problem. The main advantage of these networks is in the fact that a single AND/OR network visualizes the detailed representation of all flexibilities for a considered part or a product. By using AND/OR connectors and links, a network can easily be traversed and a feasible alternative process plan can be generated.
Knowledge-based representation is the second method used for representing process plans [8,39]. Here, a process plan is represented using vectors with n bits of data information. Each bit represents a machining operation and the order of bits form an operation sequence. As mentioned, each machining operation requires a set of a machine candidate, a cutting tool candidate and a TAD candidate which are used in order to perform a given machining operation. Therefore, apart from an operation vector, this representation also consists of a machine vector, a tool vector and a TAD vector.
To form a process plan, we considered that each machining feature on a part contains a certain number of operation units [32]. Operation units represent single machining operations associated with its candidate resources and can be formulated using the following expression: where mo i represents a machining operation i, while M j , CT k and TAD l stand for a machine candidate, a cutting tool candidate and a TAD candidate for an operation i, respectively. Operation units become building elements of an operation sequence that can be represented using modified knowledge-based approach similar to the one previously mentioned. Figure 1 shows a simple example of a representation of process planning problem. In Figure 1, we have the total of 5 operation units that match the permutation of 5 numbers. They are shown in the dashed boxes in the first row of the represented process plan (left picture). Position 1 in the vector contains the operation unit 3 which further include M 1 -machine 1, CT 6 -cutting tool 6 and TAD +x -TAD "+x" as candidates for mo 3 -machining operation 3. All these machining operations form an operation sequence vector called operations[n] (n is a total number of machining operations). Candidate solutions are assigned to each machining operation and presented in the right picture. Therefore, randomly selected candidates form machine vector, cutting tool vector and TADs vector, respectively.

Manufacturing Relationships and Precedence Constraints
Process planning is a precondition for executing a machining process. According to manufacturing and geometric specifications of machining features, manufacturing relationships occur followed by potential interferences between certain consecutive operations. Therefore, in order to machine a feature, process planners are required to pay thorough attention to manufacturing relationships which strongly affect the feasibility of process plans. Consequently, generating operation sequence without considering manufacturing relationships may result in a situation where executing one operation may make it difficult or impossible to execute others. These relationships form precedence constraints in the PPO problem which affect the sequence of all machining operations required to machine a certain product or a part. With the emphasis on machining, two distinctive types of precedence constraints based on relationships among machining operations are hard and soft constraints [40].
Since the PPO problem considered in this study is formulated as an optimization problem, precedence constraints need to be respected in order to find feasible and optimal process plans. Generally, as the name implies, hard constraints are more stringent compared to soft constraints. They directly affect the feasibility of process plans and process planners must ensure the consistency with hard constraints. On the other hand, soft constraints have the influence solely on the quality, cost and efficiency of a feasible process plan and therefore can be violated in some cases.
According to [39,40], hard and soft constraints can be divided into several types. Figure 2 represents the schematic illustration of precedence constraints based on manufacturing relationships among features.

Representation and Modelling of Precedence Constraints
Precedence constraints define precedence relationships between features/machining operations whose purpose is to help in finding a feasible operation sequence. With the emphasis on hard and soft constraints, we will consider two different approaches in this research. One approach will focus only on hard precedence constraints and the other will consider both types of constraints, hard and soft. Accordingly, two different prismatic parts adopted from the literature sources are adopted to identify precedence relationships and define types of constraints based on relationships or interactions between machining operations.
The simplest way to graphically illustrate precedence relationships is with operation precedence graphs (OPGs), which are effectively proposed in [8,41]. Assuming the fact that the emphasis will be placed on two different approaches, we will adopt two different OPGs. For hard constraints only, representation will be based on classical connected OPGs (cOPGs), while for both hard and soft constraints representation of precedence relationships will be based on disconnected graphs (dOPGs).
In order to represent aforementioned constraints and relationships, we adopted the first part from Guo et al. [42]. Its 3D solid model with assigned features and operations, as well as the cOPG that depicts only hard constraints, are shown in Figure 3. According to this figure, the total number of machining operations for prismatic part 1 is 20 and can be defined as the set of machining operations [mo 1 , mo 2 , mo 3 , . . . , mo 20 ]. Each machining operation from the cOPG is represented as a node while connections among these nodes are directed edges that show precedence relationships among them. An operation sequence that satisfies these precedence relationships is considered a feasible operation sequence. The best way to manipulate such graphical data in a selected programming environment is by converting it to a matrix form. In that sense, we adopted the adjacency matrix [8]. The adjacency matrix is generally used to represent nodes of a graph. To convert graph information to an adjacency matrix, we used the following expression: where PR ij stands for a precedence relationship between machining operation i (regarding the row) and machining operation j (regarding the column). For hard constraints only, the PR ij values are binary, ones and zeroes. The first rule for this case is when PR ij = 1 and PR ji = 0. Here, a directed edge connects the nodes mo i and mo j pointing from mo i towards mo j and, therefore states that machining operation i has to be performed prior to machining operation j. Otherwise, PR ij = 0. Another rule is where PR ii = 0 for all operations from the operation set (here [mo 1 , . . . , mo 20 ]). The expression n × n means that the adjacency matrix is a square matrix where n stands for the total number of machining operations. It can be noticed from Figure 4, that the first row for mo 1 regarding the machining operation 1 has six precedence relationships that equal number 1, meaning that the machining operation mo 1 has to be performed prior to mo 2 , mo 3 , mo 5 , mo 6 , mo 11 and mo 18 . The same rule stands for each row in the matrix. To put an emphasis on both hard and soft constraints we included another prismatic part proposed by [3]. Three-dimensional (3D) solid model of the prismatic part 2 with assigned features and machining operations along with the corresponding dOPG is shown in Figure 5. Figure 5b presents the disconnected dOPG with 9 different precedence relationships among machining operations. Only node mo 5 does not have precedence relationships with other nodes. There are 4 hard constraints and 5 soft constraints depicted with red and blue edges. Accordingly, two different numbers have to be used in order to form a reliable adjacency matrix. Herewith, number 1 stands for a hard constraint and number 2 in the matrix assumes that node mo i and node mo j form a soft constraint. The formulation and rules are the same as for the cOPG. The operation precedence matrix for the prismatic part 2 is given in Figure 6. Hard constraints are the ones whose violation must be avoided and therefore have total priority over soft constraints. This means that relationships among machining operations assigned with number 2 in the matrix can in some cases be neglected if they get into conflict with hard relationships assigned with number 1.

Mathematical Model of Process Planning
So far, the most popular evaluation criteria for the PPO problem are the minimization of the total machining time and the minimization of the total machining cost. The criterion of the total machining time has been successfully used in [9,11,20,21,25]. The machining time is composed of five time factors: machine processing time, tool processing time, transportation time (between machines), total tool change time and total setup change time. On the other side, the minimization of the total machining cost has been more frequently considered in scientific studies [4,6,8,10,11,[14][15][16]24,29,31,32]. This criterion is adopted here, and consists of five cost components that are described below.
1. Total machine cost (TMC) represents the total cost of all machines that are used in a process plan and is calculated using the following equation: where j stands for operations[i].Machines (i.e., j = 1 . . . nMach, where nMach stands for total number of machines), n is the total number of machining operations while MCI j represents the machine candidate cost index for using a machine j, a constant cost value for a specific machine. 2. Total cutting tool cost (TCTC) is the total sum of cutting tool costs used in a process plan and is calculated in the following way: where j stands for operations [i].Tools (j = 1 . . . nTools, where nTools stands for total number of cutting tools), CTCI j stands for the cutting tool cost index for using a cutting tool j, a constant cost value for a specific cutting tool. 3. Total machine change cost (TMCC) and the number of machine changes (NMC) are important cost factors considered when two successive machining operations in operations[n] are performed on different machines. Number of machine changes is computed as: where MCC represents the machine candidate change cost index while operations[i].Machines stands for the machine ID used for performing the machining operation i from operations[n]. 4. Total cutting tool change cost (TCTCC) and the number of cutting tool changes (NCTC) are similarly significant cost factors compared to those of machines. Cutting tool changes do not occur only in situations when same cutting tool and same machine are used for machining two successive machining operations. In other cases, tool changes are required. Accordingly, number of cutting tool changes is computed as follows: , where CTCC is the cutting tool change cost index and operations[i].Tools stands for the cutting tool ID used for performing the operation i from operations[n]. 5. Total setup cost (TSC), number of setup changes (NSC) and the number of setups (NS) are the cost factors considered when two successive machining operations are not performed on the same machine and using the same tool approach direction. This stands for 3-axis machines which are considered in Section 4. Firstly, the NSC is calculated as follows: , The corresponding NS is calculated as: Number one assumes the starting machine setup which is then added to the corresponding number of setup changes calculated by Equation (11). Lastly, the TSC costs are computed as: where SCC represent the setup change cost index and operations[i]. TADs stands for the TAD ID used for performing the operation i from operations[n]. Similar to tool changes, setup changes do not appear only when the same TAD and the same machine are used when performing two successive machining operations. In other cases, the change is required. 6. Additional penalty costs (APC) and the number of violating constraints (NVC) are the cost factors included in the part of our study focused on soft and hard constraints. Even if soft constraints are allowed to be violated, number of violations is subject to minimization. Firstly, the number of violated soft constraints is calculated using the following equation: where operations[i] and operations[j] stand for two consecutive machining operations in operations[n] vector. The APC is thereby calculated as: where penalty assumes the fixed penalty cost index which is applied for each violated soft constraint. In this study, the penalty cost applies to prismatic part 2 that involves both hard and soft constraints. As the conflicts between these two types of precedence constraints may occur, hard constraints become the priority. In that case, soft constraints that are numerically assigned with number 2 in Figure 6, have to be violated. 7. Total weighted machining cost (TWMC) sums up all the previous cost factors in the following equation: where w 1 − w 6 are the weight coefficients used for experimental studies. The APC costs are not included in Equation (17) in those studies which consider only hard precedence constraints. 8. Fitness function is maximized in order to verify the total weighted machining costs of a process plan and is defined as follows:

Traditional Grey Wolf Optimizer
Grey wolf optimizer (GWO) is proposed by Mirjalili et al. [37]. It is a modern swarmbased algorithm whose inspiration derives from the social intelligence of grey wolves.
Grey wolves (lat. Canis lupus) are predatory animals which brings them to the top of food chain in animal world. The advantages of social life of grey wolves mostly reflect in social hunting, group care about infants and territorial defense against dangerous forces outside the pack. A strict social hierarchy among members of the pack is what brings grey wolves to the fore. From the most dominant to the most submissive wolves there are alphas, betas, deltas and omegas. According to the social dominance hierarchy, a mathematical model of GWO was designed. Three main stages of grey wolf hunting process, such as tracking, encircling and attacking a prey, are included in this model.
According to the social dominance hierarchy previously described, a mathematical model of GWO was designed. Three main stages of grey wolf hunting process, such as tracking, encircling and attacking a prey, are included in this model.
When designing the social hierarchy in the GWO, the three best solutions in a population of individuals represent alpha (α), beta (β) and delta (δ) wolves, respectively. All the other candidates are considered to be omegas (ω). The hunting (optimization) process is guided by the three fittest wolves.
The first mechanism in the GWO is encircling of prey which is defined according to the following equations: where it stands for current iteration, When modelling hunting behaviour, the fittest wolves, the alpha (α), the beta (β) and the delta (δ) have the main role. It is assumed that these wolves respectively have better knowledge about potential location of the prey. Therefore, omega wolves (ω) update their positions based on the positions of the alpha (α), the beta (β) and the delta (δ). This process is mathematically formulated by the following equations: where → X represents a position of a grey wolf over a cource of iterations it updated according to the distances from three fittest wolves, alpha, Search for prey, or divergence from prey, is another component of the GWO that emphasizes exploration. It is defined by the coefficient vector → C which ranges from 0 to 2. Its purpose is to assign random weights to a potential prey and give wolves harder or easier way to reach it. In order to use global capacities of the GWO and avoid global optima, the value of the vector → C does not decrease linearly compared to the vector → A. Grey wolf optimizer has been applied to various types of optimization problems to this day. One of its main advantages is less parameter tunning since all the vectors above mentioned are based on random values. However, main drawbacks of the GWO are low convergence rate and entrapment in local optima [20]. To overcome these issues, researchers have proposed different ways to improve or enhance the traditional GWO. Ahmed et al. [43] proposed niching GWO to deal with multi-modal optimization problems. To maintain balance between exploitation and exploration and avoid premature convergence, authors incorporated best features of the PSO algorithm and a local search technique. Yue et al. [44] developed a novel hybrid algorithm based on the GWO and fireworks algorithm (FWA) in order to combine both global abilities of FWA and local abilities of GWO. Wang and Wang [45] used quantum computing principles and operations of differential evolution with grey wolf optimizer in order to improve performance when dealing with NP-complete combinatorial 0-1 knapsack problems. Martin et al. [46] proposed a novel discrete GWO algorithm which uses update rules to distinguish between exploitation and exploration stage. Jiang and Zhang [47] applied the GWO algorithm on job shop and flexible job shop scheduling problems. They embedded crossover operation and adaptive mutation method to maintain the search within discrete domain, and avoid premature convergence and local optima. Variable neighborhood search was also introduced to improve explorative capacities of the proposed GWO. Qin et al. [48] implemented hybrid discrete multi-objective GWO to solve the casting production scheduling problem where makespan, the total production cost and the total delivery delay time were used as objective functions. A strategy based on reducing job transportation time and processing time are adopted to improve initial solutions. Additionally, improved tabu search algorithm was incorporated into the GWO to improve performances of the proposed method. Premkumar et al. [49] developed a multi-objective GWO to solve the brushless direct current (BLDC) motor design problem. Firstly, the analytical model of the BLDC motor design problem that belongs to highly non-linear electromagnetic optimization problems is presented. MOGWO was verified on standard benchmark functions, and latter applied to mono and multi objective optimization of BLDC motor design problem.
In the next section, the hybrid grey wolf optimizer based on traditional GWO and GA operations is proposed for solving NP-hard process planning combinatorial problem.

HGWO Methodology
When designing mathematical model of the HGWO method, certain adaptations of previous Equations (19), (20), (23) and (24) have to be considered. Firstly, the distances of gray wolves from the alpha (α), the beta (β) and the delta (δ) wolf for the machine vector, the tool vector and the TAD vector respectively, are formulated as follows: where → X m stands for the machine position vector of a grey wolf, X δ,tad are the alpha (α), the beta (β) and the delta (δ) wolf for the TAD vector. Furthermore, the position of a gray wolf in the iteration it + 1 is formulated on the basis of the machine vector, the tool vector and the TAD vector relative to the position vectors of the alpha (α), the beta (β) and the delta (δ) wolves. This can be expressed with the following expressions: , where → X m (it + 1), the delta (δ) wolf respectively, → X 1,t , → X 2,t and → X 3,t are the updated tool vectors according to the distance from the alpha (α), the beta (β) and the delta (δ) wolf respectively, and → X 1,tad , → X 2,tad and → X 3,tad are the updated TAD vectors according to the distance from the alpha (α), the beta (β) and the delta (δ) wolf respectively.
So far, no recent discoveries have shown the proposal of the GWO for solving process planning optimization problem. In this paper, we present a novel hybrid grey wolf optimizer (HGWO) approach to address the PPO problem with various precedence constraints focusing on prismatic parts. The pseudocode of the proposed approach can be seen in Figure 7. It includes the following components: (1) knowledge-based representation of process plans (positions of grey wolves), (2) initializing population (pack of wolves), (3) applying constraint handling heuristic, (4) fitness evaluation of positions of wolves, (5) traditional GWO steps for updating positions of wolves, (6) tournament selection, (7) crossover, (8) shift mutation, (9) resource mutation, (10) fitness re-evaluation and (11) termination criteria of the HGWO.

Initialization of Population
The first step after defining the representation of solutions is to initialize a starting population of individuals, here described as a pack of wolves. As mentioned earlier, the modified knowledge-based representation is adopted in this study. For the adopted representation of individuals, a population can be initialized. To ensure the feasibility of individuals, the search methodology has to include an appropriate approach for converting infeasible individuals to a feasible domain. Therefore, we propose two heuristic algorithms for constraint handling that will be incorporated into the proposed HGWO in order to guarantee feasible solutions. The first heuristic is adopted from [8] and utilized for dealing with hard constraints only. The other heuristic was developed in [39] to address both hard and soft constraints.
Using the information about precedence constraints and operations precedence graph let us assume that a randomly generated sequence of operations for a single individual from a population would have the following form: [mo 2 , mo 5 , mo 6 , mo 8 , mo 1 , mo 9 , mo 4 , mo 3 , mo 7 ]. To complete the initialization of a process plan, randomly selected machine, cutting tool and TAD candidate are assigned to each of these machining operations to form the operation units. Here, we will focus on the operations sequence which is crucial for elimination of infeasible solutions. Using the adjacency matrix as an input information, the constraint handling heuristic for the hard constraints is developed, Figure 8. The step-by-step procedure of the constraint handling heuristic algorithm is the following:

1.
A vector operation[n] is used for representing a single solution; its adjacency matrix is given in Equation (2); 2.
A variable pt 1 (pt 1 = n) representing an index pointer of the vector operations[n] is generated; 3.
Repeat the procedure within the for loop until pt 1 = 0; 4.
Identify the machining operation at the point pt 1 = n; 5.
For the identified machining operation calculate the sum of all number values in a corresponding row of the adjacency matrix ∑ n i=1 (PR ij ); Then, check the conditions: Leave the machining operation at the position pt 1 ; (b) Exclude the row and the column of the Adj_matrix that match the value of the machining operation at the position pt 1 ; instead of removing the row and the column, assign "not a number"; (c) Initialize n = n − 1 and move pointer pt 1 one place left; II.
Identify the machining operation at the position pt 2 ; (c) Calculate the sum ∑ n i=1 (PR ij ) of all number values in the corresponding row that matches the identified machining operation at the position pt 2 ; (d) Repeat the following steps while pt 2 > 0; (e) If ∑ n i=1 (PR ij ) = 0 Swap the machining operation at the position pt 1 with the machining operation at the position pt 2 ; Remove the row and the column of the corresponding Adj_matrix which match the value of the machining operation at the position pt 2 ; assign "not a number"; Break the loop and move to the step (3); Move to the step (d);
Reset the Adj_matrix and apply all the steps for the next individual.  4 , and select the next operation, mo 9 ; (e) Swap positions of mo 7 and mo 9 , exclude row and column of mo 9 , and select the next operation, mo 8 ; mo 1 is skipped according to the Step 5.II.f); (f) Swap positions of mo 7 and mo 8 , exclude row and column of mo 8 , and select the next operation, mo 1 ; (g) Swap positions of mo 7 and mo 1 , exclude row and column of mo 7 , and select the next operations, mo 1 and mo 6 ; (h) Leave mo 6 and mo 1 as positioned, exclude row and column of mo 6 , and select the next operation, mo 5 ; (i) mo 1 goes before mo 5 and before the last operation mo 2 ; the final sequence is generated.
The step-by step constraint handling procedure begins with the adjacency matrix given in Figure 8a. The procedure is based on identifying index pointers pt 1 and modifying the sequence within the for loop. The first selected pointer pt 1 = n = 9 denotes the last machining operation in the matrix, mo 7 . The sum of all values in the row 9 meets the second condition in the Step 5 meaning that ∑ n i=1 (PR ij ) = 0. Then, the Step 5.II.a) starts with initializing the second pointer pt 2 = pt 1 − 1 = 8 that denotes mo 3 , Figure 8b. The same condition is checked and since the sum ∑ n i=1 (PR ij ) equals zero, the operations mo 7 and mo 3 are swapped and the row and the column of the mo 3 are assigned with the non-numbering value. Operation mo 3 than takes the final position 8 in the sequence, Figure 8c. The procedure continues with the Step 6 that reduces the pointer pt 1 to 8 and repeats from the Step 5. Pointers pt 1 = 8 and pt 2 = 7 denote the operations mo 7 and mo 4 , Figure 8c. After checking conditions in Step 5.II, the row and the column of mo 4 are assigned with the non-numbering value and removed from the following steps, Figure 8d. From Figure 8d to Figure 8i, the handling procedure continues to remove the rows and the columns of machining operations for which the sum ∑ n i=1 (PR ij ) equals zero. The sequence of operations is modified accordingly. Figure 8i shows that the operation mo 5 has to swap place with mo 1 in order to complete the procedure of generating a feasible machining sequence.
The second constraint handling heuristic is adopted from [3]. Here, we will shortly focus on the most significant steps of this heuristic. In order to handle both hard and soft constraints, firstly, so called "linked list" is created. This list includes machining operations that form hard constraints, therefore leaving the positions of other operations unchanged. Herewith, a heuristic process is imposed on the linked list to ensure its consistency with the hard constraints. Such an approach is adopted for the part 2 shown in Figure 5. Afterwards, the next step of this constraint handling algorithm is to include the additional penalty costs for violated soft constraints in the objective function model, Equations (15) and (16). Considering the fact that soft constraints are not manipulated using the linked list, certain soft constraints may be subject to violation. In this step, soft constraints can be compromised and violated to achieve the minimal total weighted machining cost, Equation (17).

Genetic Components of the HGWO Approach
Genetic algorithms are techniques based on natural evolution of species that mimic the viewpoint of modern genetics, so-called "survival of the fittest". Apart from being fitnessoriented and belonging to the group of population-based optimization algorithms, their unique characterization is varying operations that mimic genetic gene changes and enable population individuals to evolve. These strategies are selection, crossover and mutation. The adopted selection strategy for the HGWO algorithm is the tournament selection which belongs to the group of proportionate-based selection schemes. This strategy is based on selecting a number of individuals (tournament size) from the population that will participate in the so called "tournament". The individual with the highest fitness value is considered as the fittest individual and therefore the winner of the tournament. The process of selecting individuals lasts until a new generation of individuals is created.
Crossover or recombination is the next strategy employed to recombine two selected parent individuals in order to obtain better offsprings. In this paper, a modified one-point crossover strategy is adopted to improve exploration capabilities of the proposed approach. Using the predefined crossover probability p c , the steps of this recombination procedure are the following:

•
Two wolf candidates selected from the tournaments are randomly chosen to be parent wolves; • By generating a random crossover point, two parent wolves are divided into two sections to create two child wolves; • Left section of the child wolf 1 is formed by copying left section of the parent wolf 1. Then, the right section of the child wolf 1 is formed in two parts. The first part is to find the remaining machining operations in the parent wolf 2 and copy them in the current order to the right section of the child wolf 1. Afterwards, the machine, cutting tool and TAD candidate are copied from these machining operations in the parent wolf 2 to the remaining machining operations in right section of the child wolf 1; • The identical, but inverse procedure using left section of the parent wolf 2 and remaining machining operations from the parent wolf 1 is performed to produce the child wolf 2.
To improve exploration capabilities of the HGWO and avoid premature convergence, two appropriate mutation strategies are introduced. Here, the shift mutation and resource mutation strategy are used. The first strategy, takes two random genes (operation units) of a randomly selected wolf from the population (pack) and exchanges them. This procedure impacts the feasibility of the mutated position of the wolf.
The second mutation strategy called "resource mutation" does minor changes to a machine, cutting tool and TAD candidate, respectively. The steps of the resource mutation strategy (example for a machine vector) are the following:

•
Randomly select a crow candidate for resource mutation; • Randomly select the mutation point, i.e., machining operation that matches certain position in the vector; • Check the availability of machines in the machine candidate list for the selected machining operation; • Randomly select one of the available machines as the current machine; • Identify all other machining operations that have the current machine in their machine candidate list; • Assign the same alternative machine candidate as the current machine for remaining machining operations; • Repeat the same steps for the CT vector and TAD vector using CT candidate list and TAD candidate list respectively.
Besides the above described genetic strategies within the HGWO approach, we also included inertia weight coefficient in order to achieve additional control of exploration and exploitation. The inertia is decreased linearly thereby emphasizing global exploration ability in the initial stages of the search process and moving towards local exploitation abilities in the latter stages. Convergence curves in Chapter 5 graphically show diversification of the search (primarily HGWO and HybGA curves) in the initial stages meaning that the exploration is at hand, whereas the latter stages clearly show less diversified search where local abilities of metaheuristics are more used. In that sense, the mathematical Equations (25)- (27) defined in the Section 4.2 will have the following form: where w represents inertia weight that is decreased linearly using the next equation: where w max and w min are maximal and minimal weight, MaxIt is the maximal number of iterations and it is the current iteration.

Experimental Studies and Results
The proposed HGWO algorithm was programmed in MATLAB environment and executed on a Windows 10 operating system by a 1.99 GHz Intel i7 processor and 8 GB RAM computer. Its performance was verified on three experimental studies that consider three different prismatic parts. The first experiment employed the prismatic part 1 taken from [42]. The 3D solid model with associated features and its cOPG are represented in Figure 3. The second experiment for testing the HECSA algorithm is the prismatic part 2 taken from [3]. Its model representation and dOPG are presented in Figure 5. The prismatic part 3 in the third experiment is adopted from [6]. Figures 9 and 10 depict the semi-transparent 3D model with its dOPG and adjacency matrix, respectively. Compared to the first and the second experimental study which considered prismatic parts with hard, and both hard and soft constraints respectively, precedence relationships in this model require certain operations to be performed in the same setup. In that sense, number 2 is added to the dOPG (e.g., mo 10 precedes mo 3 and both should be performed in the same setup). The prismatic part 3 is modeled for one condition in which all resources are available.  The detailed information about features and candidates for machines, cutting tools and TADs for prismatic part 1 is given in Table 1. Precedence relationships among machining features for prismatic part 1 are presented in Table 2. Similarly, information about features and resource candidates, as well as information about precedence relationships for prismatic part 2 are given in Tables 3 and 4, respectively. Lastly, information about features and resource candidates for prismatic part 3 are shown in Table 5. Available machines and cutting tools as well as the machining cost information about three prismatic parts are given in Tables 6 and 7, respectively.
Parameter tuning may be one of the most time-consuming challenges when testing the performances of metaheuristic algorithms. In this paper, we performed several manual tests of input parameters of the HGWO and finally we adopted the following: pack size: N = 100, maximal number of iterations: MaxIter = 700, tournament size: TourSize = 5, crossover probability: p c = 0.8, shift mutation probability p m1 = 0.4, resource mutation probability: p m2 = 0.6 and maximal and minimal inertia weights are w max = 1.2 and w min = 0.2. The adopted optimization objective is to minimize the total weighted machining cost defined in Section 3.4. The output results included minimum and maximum results achieved in 20 runs, and mean values of the obtained results.
Nine holes arranged as a replicated feature A compound hole Reaming (mo 19 ) The optimal process plans for minimal machining cost and three conditions concerning prismatic part 1 are shown in Table 8. In the first condition, all resources are available. The second condition excludes costs for cutting tools and cutting tool changes. The third condition is the same as the second, with unavailability of machine 2 and cutting tool 8.
Firstly, we will discuss the results for the prismatic part 1. For the first condition where all resources are available, the HGWO obtained the minimal TWMC of 2470 cost units in 20 runs. The total of six results have values bellow 2500 cost units. According to Table 9, the minimal TWMC of 2470 is one of the best results compared to the results obtained by other algorithms. Only the approach reported in [15] achieved better minimum. The average cost obtained by the HGWO is better than the average cost achieved by cPSO [20], PSO [42], and mACO [17]. However, the maximal value of 2805 shows a lack of consistency in achieved results compared to other metaheuristics. Considering the second condition where CT cost and CT change cost are not included, the TWMC of 1990 cost units appears to be the best minimal result. As shown in Table 9, HGWO outperforms all other metaheuristic approaches in terms of the minimal TWMC and the mean TWMC. On the other hand, the achieved maximal value was better compared to ACO [15], TSACO [16] and IWD [32] approaches. Lastly, the HGWO showed best performances for the third condition in which machine M 2 and cutting tool CT 8 are not available. Best minimal TWMC of 2490 cost units is found. Additionally, regarding the mean and the maximal values of the TWMC, the HGWO demonstrated much better consistency and effectiveness than the considered modified and hybrid approaches. Milling (mo 7 ) f7 (mo 7 ) is before f8 (mo 8 , mo 9 , mo 10 ) for the datum interactions f8 Drilling (mo 8 ) mo 8 is before (mo 9 and mo 10 ); mo 9 is before mo 10 for the fixed order of machining operations Reaming (mo 9 ) Boring (mo 10 ) f9 Milling (mo 11 ) f9 (mo 11 ) is before f10 (mo 12 , mo 13 , mo 14 ) for the datum interaction f10 Drilling (mo 12 ) mo 12 is before (mo 13 and mo 14 ); mo 13 is before mo 14 for the fixed order of machining operations; f10 (mo 12 , mo 13 , mo 14 ) is before f11 (mo 15 , mo 16 ) and mo 12

Manufacturing Interactions Precedence Relationships Type of Precedence Constraints
Tool-feature interaction mo 1 must be executed before mo 2 Hard constraints Feature-datum interactions mo 6 must be executed before mo 7 mo 10 must be executed before mo 11 mo 13 must be executed before mo 14 Thin-wall interactions mo 9 must be executed before mo 8 Soft Constraints mo 12 must be executed before mo 10 Material removal interactions mo 8 must be executed before mo 9 mo 10 must be executed before mo 12 mo 13 must be executed before mo 14 mo 3 must be executed before mo 4  Reamer -Reamer CT 10 Boring tool -Slot cutter   1  18 11  2  6  12 13 17  5  3  7  8  9  19 10 20 14  4  15 16  Machine  2  2  2  2  2  2  2  2  2  2  2  2  2  2  4  4  4  1  1  1  Cutting tool  8  7  8  8  8  3  9  8  7  8  8  3  9  9 10 10 10 2  1  3  18 17  5  11  2  6  12 13 14 19 20  7  8  9  10  4  15 16  Machine  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  1  1  1  Cutting tool  7  7  7  7  7  7  6  7  4  9   Furthermore, in Figure 11 we presented convergence curves for prismatic part 1. As can be seen, convergence of six different metaheuristic approaches is compared considering three optimization conditions. Convergence of the proposed HGWO is compared with the convergence of three traditional algorithms, PSO [42], GA [14] and GWO, as well as the hybrid and the modified metaheuristics, HybGA [8] and cPSO [20]. From the information in Figure 11a, considering the first condition, GWO, GA and HybGA showed much faster converegence compared to other algorithms. However, the GWO and GA converged towards local optima, whereas the HybGA succeded in finding much better solution. Although the HGWO showed much slower convergence than the HybGA for example, it converged towards the minimal value of 2470. This value may not be considered the global optimum, since the slightly better result in the literature has been reported so far. It could be safely argued that the HGWO converged towards the near-optimal solution. From Figure 11b and the second condition of the first study, it can be noticed that the convergence of the HGWO is similar to the one of the HybGA. After approximately 200 iterations, these approaches converged towards the optimal and near-optimal result. PSO, like CPSO, converged towards the good solutions, but with much slower convergence. GA and GWO got trapped in the local optima in the early stages. Lastly, the convergence curves for the third condition are given in Figure 11c. It can be seen that all algorithms in this analysis converged towards best solutions. The reason lies in the fact that the third condition excludes certain solution candidates as well as cost assigned with them. HGWO showed slightly slower convergence than GA, GWO and HybGA, but, on the other side, succeded in finding the new global optimum of 2490 cost units for TWMC. This result showed the effectiveness of the HGWO for the third condition. Best process plans for two conditions for prismatic part 2 are given in Table 10. After adopting similar parameters as in the previous study, we compared the results with six different approaches from the literature [8,14,20]. The comparative results are shown in Table 9. The proposed HGWO performed well by achieving optimal TWMC of 1328 cost units eight times under the first condition. The mean value was slightly greater compared to the same result achieved by TS and ACO, but much better compared to other algorithms. Similar conclusion can be made for the maximal TWMC where the HGWO showed better result compared to the HybGA, TS, SA, GA and CPSO. Under the second condition, the HGWO found the minimal TWMC of 1170 cost units which appeared 19 times in 20 runs. Besides the ACO approach, the HGWO achieved the best maximal and the best mean results compared to all other algorithms. Figure 12 presents the convergence curves for two optimization conditions of the second experimental study concerning prismatic part 2. The same six metaheuristics are considered in the convergence analysis. As shown in Figure 12a for the first condition, the HGWO and the HybGA show similar performances, with the HGWO showing slower convergence compared to the hybrid metaheuristic algorithm. Both approaches converged towards the global TWMC of 1328 cost units. GA and GWO got trapped in the local optima in early stages, whereas the CPSO and the PSO showed progress in converging towards good solutions but not as good as the ones found by the HGWO and the HybGA. Similar performances can be seen for the second condition shown in Figure 12b. Since the second condition excludes certain alternatives and costs, the HGWO showed faster convergence than for the previous condition. Like the HybGA, HGWO converged towards the global TWMC of 1170 cost units. It took HGWO approximately 150th iterations to converge towards the global optimum compared to the HybGA which was faster to some extent. However, compared to the traditional GA, GWO, PSO, as well as the CPSO, the HGWO showed more superior performance for the second condition.  The optimal process plans for prismatic part 3 are obtained and presented in Table 11. The same input parameters are used for this experimental study. Table 9 shows the comparative results for the minimal, the maximal and the mean TWMC. In 20 runs, the HGWO obtained the TWMC below 1900 cost units eight times in total, and bellow 2000 cost units eleven times. The minimal TWMC of 1815 cost units is better than the minimal values obtained by the GA, PSO, as well as the HGASA. The HGWO outperforms these algorithms in terms of the best maximal and the best mean results. Figure 13 shows the convergence curves for the third experimental study. It can be noticed that the HybGA again shows the best convergence rate and by far outperforms other algorithms in this regard. As far as the HGWO is concerned, it shows much slower convergence compared to the first and the second experimental study. The HGWO obtained the TWMC of 1815 cost units after 500 iterations. Despite the slow convergence rate, the proposed HGWO algorithm achieved near-optimal solution for this prismatic part. Traditional GWO and especially GA also showed good convergence rates but they got trapped in the local optima. CPSO and PSO performances could not match the performances of the previous algorithms. Table 11. Optimal process plans for prismatic part 3 for two conditions. 11  13  12  14  8  9  10  7  1  3  5  4  2  6  Machine  2  2  2  2  2  2  2  2  2  2  3  3  3  3  Cutting tool  5  5  5  5  5  1  6  6  6  3  8  9  2

Time Efficiency of the HGWO
According to the obtained experimental results, the HGWO approach showed very good effectiveness and flexibility by finding optimal and near-optimal solutions of the process planning optimization problem. From the perspective of the solutions quality, the comparative results in Table 9 show that the HGWO obtained the near-optimal TWMC value for the first condition, and the new optimal values for the second and the third condition of the first experimental study. In terms of the consistency, the HGWO showed good mean values for all three conditions surpassing the mean results of several algorithms reported in the literature. The optimal TWMC cost values are obtained for both conditions of the second experimental study with much better consistency that is reflected in very good mean results. In the third experimental study, the HGWO approach obtained nearoptimal results of TWMC with good consistency where mean results surpass those results obtained by traditional algorithms.
Although the proposed HGWO approach demonstrated flexibility and effectiveness in terms of the solution qualities, there are certain limitations regarding its computational time efficiency. Due to the fact a small number of researchers included computational times in their studies, a detailed comparative analysis could not be performed in this regard. In Table 12 we presented the mean computational times of six metaheuristic algorithms whose convergence curves were shown in Figures 11-13 for three experimental studies, respectively. They are expressed in seconds required to complete a single run, and the results in Table 12 are the mean computational times of 20 runs. From the perspective of the HGWO efficiency, better mean computational times are achieved against traditional PSO and modified cPSO approach. On the other side, they exceed the mean computational times of three algorithms, GWO, GA and HybGA. Similar values for the HGWO and the PSO are obtained for the prismatic part 2. Although some authors suggest less amount of time to be reasonable for large-sized NP-hard problems [50], we may argue that approximately 1 to 1.5 min of the running time is a reasonably well amount of time needed to solve different instances of an NP-hard problem such as the PPO. GA and GWO as traditional algorithms obviously show the lowest computational times in all studies, since the modified and hybrid algorithms require more time to perform additional operations. One of the reasons for lower efficiency of the HGWO compared to the HybGA is the difference in input parameters, where HGWO considered larger probabilities for crossover and mutation strategies, as well as larger tournament size, compared to HybGA. According to the convergence curves in Figures 11-13, the convergence rate of the HGWO is slower in the first stages of the search process, especially for the first conditions of prismatic parts 1 and 2, as well as prismatic part 3. For example, curve in Figure 11a shows slower rate of the HGWO in the first 300 iterations, meaning that the probabilities of crossover and mutation are much higher than in the latter iterations, which leads to the increase in computational time of the HGWO. Time efficiency of the HGWO can also be evaluated by the function execution times that are called in the programming environment. Since the MATLAB was used to implement the HGWO, we provide the execution times of several functions coded in MATLAB which are presented in Table 13. Apart from the functions that consider calling six different algorithms considered in experimental studies, fitness evaluation and constraint handling functions are also presented. Two time values, total and self-time are shown. Total time represents the total time spent in a function including all child functions, while the self-time means the total time spent in a function without any time spent in child functions. All time values in Table 13 are obtained after 20 runs and consider only the first experimental study. It can be noticed that calling the HGWO algorithm requires more time compared to almost all other algorithms except the CPSO. The big difference in total and self-time means that the HGWO requires calling many child functions, including fitness evaluation and constraint handling algorithm. According to the computational times and function execution times, we can argue that the HGWO demonstrated good time efficiency. However, this can also be considered as the main disadvantage of the proposed approach since the comparative results clearly show that the HGWO falls behind other algorithms. The fact that the constraint handling heuristic, the fitness evaluation, the tournament selection, the crossover, the shift mutation, and the resource mutation represent a separate MATLAB functions that the HGWO has to call in order to perform a single iteration, the time efficiency should be treated with enough attention. With effectiveness and flexibility as the main advantages, and time efficiency as the main disadvantage of the HGWO approach, potential improvements in terms of time efficiency and convergence rate are beyond the scope of this paper and may be a subject for future research studies.

Conclusions
A hybrid grey wolf optimizer is proposed to solve the NP-hard process planning optimization problem based on precedence constraints. The problem considers operation sequencing task along with the optimal selection of appropriate machines, cutting tools and tool approach directions that define a process plan. Modern knowledge-based representation scheme is used to represent process plans in a population of alternatives and the operation precedence graph approach with adjacency matrices is adopted to deal with precedence relationships among machining features. According to these relationships, precedence constraints are defined. To deal with the precedence constraint in an appropriate manner, we adopted constraint handling heuristic algorithm to deal with hard constraints, as well as both hard and soft constraints. In that sense, the feasibility of process plans was ensured.
To improve the performance of traditional grey wolf optimizer, the strategies of genetic algorithm were adopted. Tournament selection is utilized for selecting fittest individuals and classical two-point crossover strategy is applied on selected individuals to generate wolf offspring. Shift and resource mutation strategies preform minor changes to vector candidates and complete the evolutionary steps of the HGWO approach. To additionally achieve balance between exploration and exploitation inertia weight coefficients from the PSO are added to the mathematical model of the HGWO approach.
For the evaluation of process plans, traditional total weighted machining cost was adopted as an optimization criterion. The detailed mathematical model is presented and the criterion was formed for different conditions covered in three experimental studies which consider three prismatic parts proposed in the literature.
After manually determining input parameters of the proposed HGWO approach, popular results obtained by traditional and modern approaches were used to perform the detailed comparative analysis. Three experimental studies are conducted to test the effectiveness and efficiency of the HGWO approach. The detailed representation of optimal process plans and comparative results including the minimal, the maximal and the mean TWMC for 20 algorithm runs were given. The HGWO approach demonstrated flexibility and effectiveness by solving the process planning optimization problem and finding the optimal and the near-optimal results.
Although effectiveness and flexibility showed as the main advantages of the HGWO, time efficiency of the proposed algorithm showed certain limitations. The HGWO succeeded in finding optimal and near-optimal process plans in a reasonable time period, but the comparative analysis of computational times, as well as function execution times in MATLAB programming environment, confirm good efficiency of the HGWO. However, the HGWO falls behind certain hybrid and traditional algorithms in this regard.
One of the directions for future research may be focused on improving efficiency and convergence of the HGWO, primarily through consideration of different genetic strategies. Additionally, optimization of input parameters has already been reported in some studies and should be taken into account in order to further improve performances of the HGWO. Possible contribution to the area of process planning optimization can be directed towards the optimization of integrated process planning and scheduling. Likewise, the recent advances in the field of Industry 4.0 may consider the implementation of the HGWO as a cloud service within the Smart factory.